Skip to content

Commit 32b0c2b

Browse files
murrayrmmurrayrm
authored
Merge pull request #193 from Sup3rGeo/zgrid-cherry
Discrete omega-damping plot and tweaks. Took a look at the revisions and tested out on MacOS 10.14, Python 3.6 + verified that ReadTheDocs and TravisCI both pass. All looks good. Merging into master.
2 parents 72f427d + 89dc937 commit 32b0c2b

File tree

3 files changed

+216
-105
lines changed

3 files changed

+216
-105
lines changed

‎control/grid.py‎

Lines changed: 182 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,182 @@
1+
import numpy as np
2+
from numpy import cos, sin, sqrt, linspace, pi, exp
3+
import matplotlib.pyplot as plt
4+
from mpl_toolkits.axisartist import SubplotHost
5+
from mpl_toolkits.axisartist.grid_helper_curvelinear import GridHelperCurveLinear
6+
import mpl_toolkits.axisartist.angle_helper as angle_helper
7+
from matplotlib.projections import PolarAxes
8+
from matplotlib.transforms import Affine2D
9+
10+
class FormatterDMS(object):
11+
'''Transforms angle ticks to damping ratios'''
12+
def __call__(self,direction,factor,values):
13+
angles_deg = values/factor
14+
damping_ratios = np.cos((180-angles_deg)*np.pi/180)
15+
ret = ["%.2f"%val for val in damping_ratios]
16+
return ret
17+
18+
class ModifiedExtremeFinderCycle(angle_helper.ExtremeFinderCycle):
19+
'''Changed to allow only left hand-side polar grid'''
20+
def __call__(self, transform_xy, x1, y1, x2, y2):
21+
x_, y_ = np.linspace(x1, x2, self.nx), np.linspace(y1, y2, self.ny)
22+
x, y = np.meshgrid(x_, y_)
23+
lon, lat = transform_xy(np.ravel(x), np.ravel(y))
24+
25+
with np.errstate(invalid='ignore'):
26+
if self.lon_cycle is not None:
27+
lon0 = np.nanmin(lon)
28+
lon -= 360. * ((lon - lon0) > 360.) # Changed from 180 to 360 to be able to span only 90-270 (left hand side)
29+
if self.lat_cycle is not None:
30+
lat0 = np.nanmin(lat)
31+
lat -= 360. * ((lat - lat0) > 360.) # Changed from 180 to 360 to be able to span only 90-270 (left hand side)
32+
33+
lon_min, lon_max = np.nanmin(lon), np.nanmax(lon)
34+
lat_min, lat_max = np.nanmin(lat), np.nanmax(lat)
35+
36+
lon_min, lon_max, lat_min, lat_max = \
37+
self._adjust_extremes(lon_min, lon_max, lat_min, lat_max)
38+
39+
return lon_min, lon_max, lat_min, lat_max
40+
41+
def sgrid():
42+
# From matplotlib demos:
43+
# https://matplotlib.org/gallery/axisartist/demo_curvelinear_grid.html
44+
# https://matplotlib.org/gallery/axisartist/demo_floating_axis.html
45+
46+
# PolarAxes.PolarTransform takes radian. However, we want our coordinate
47+
# system in degree
48+
tr = Affine2D().scale(np.pi/180., 1.) + PolarAxes.PolarTransform()
49+
# polar projection, which involves cycle, and also has limits in
50+
# its coordinates, needs a special method to find the extremes
51+
# (min, max of the coordinate within the view).
52+
53+
# 20, 20 : number of sampling points along x, y direction
54+
sampling_points = 20
55+
extreme_finder = ModifiedExtremeFinderCycle(sampling_points, sampling_points,
56+
lon_cycle=360,
57+
lat_cycle=None,
58+
lon_minmax=(90,270),
59+
lat_minmax=(0, np.inf),)
60+
61+
grid_locator1 = angle_helper.LocatorDMS(15)
62+
tick_formatter1 = FormatterDMS()
63+
grid_helper = GridHelperCurveLinear(tr,
64+
extreme_finder=extreme_finder,
65+
grid_locator1=grid_locator1,
66+
tick_formatter1=tick_formatter1
67+
)
68+
69+
fig = plt.figure()
70+
ax = SubplotHost(fig, 1, 1, 1, grid_helper=grid_helper)
71+
72+
# make ticklabels of right invisible, and top axis visible.
73+
visible = True
74+
ax.axis[:].major_ticklabels.set_visible(visible)
75+
ax.axis[:].major_ticks.set_visible(False)
76+
ax.axis[:].invert_ticklabel_direction()
77+
78+
ax.axis["wnxneg"] = axis = ax.new_floating_axis(0, 180)
79+
axis.set_ticklabel_direction("-")
80+
axis.label.set_visible(False)
81+
ax.axis["wnxpos"] = axis = ax.new_floating_axis(0, 0)
82+
axis.label.set_visible(False)
83+
ax.axis["wnypos"] = axis = ax.new_floating_axis(0, 90)
84+
axis.label.set_visible(False)
85+
axis.set_axis_direction("left")
86+
ax.axis["wnyneg"] = axis = ax.new_floating_axis(0, 270)
87+
axis.label.set_visible(False)
88+
axis.set_axis_direction("left")
89+
axis.invert_ticklabel_direction()
90+
axis.set_ticklabel_direction("-")
91+
92+
# let left axis shows ticklabels for 1st coordinate (angle)
93+
ax.axis["left"].get_helper().nth_coord_ticks = 0
94+
ax.axis["right"].get_helper().nth_coord_ticks = 0
95+
ax.axis["left"].get_helper().nth_coord_ticks = 0
96+
ax.axis["bottom"].get_helper().nth_coord_ticks = 0
97+
98+
fig.add_subplot(ax)
99+
100+
### RECTANGULAR X Y AXES WITH SCALE
101+
#par2 = ax.twiny()
102+
#par2.axis["top"].toggle(all=False)
103+
#par2.axis["right"].toggle(all=False)
104+
#new_fixed_axis = par2.get_grid_helper().new_fixed_axis
105+
#par2.axis["left"] = new_fixed_axis(loc="left",
106+
# axes=par2,
107+
# offset=(0, 0))
108+
#par2.axis["bottom"] = new_fixed_axis(loc="bottom",
109+
# axes=par2,
110+
# offset=(0, 0))
111+
### FINISH RECTANGULAR
112+
113+
ax.grid(True, zorder=0,linestyle='dotted')
114+
115+
_final_setup(ax)
116+
return ax, fig
117+
118+
def _final_setup(ax):
119+
ax.set_xlabel('Real')
120+
ax.set_ylabel('Imaginary')
121+
ax.axhline(y=0, color='black', lw=1)
122+
ax.axvline(x=0, color='black', lw=1)
123+
plt.axis('equal')
124+
125+
def nogrid():
126+
f = plt.figure()
127+
ax = plt.axes()
128+
129+
_final_setup(ax)
130+
return ax, f
131+
132+
def zgrid(zetas=None, wns=None):
133+
'''Draws discrete damping and frequency grid'''
134+
135+
fig = plt.figure()
136+
ax = fig.gca()
137+
138+
# Constant damping lines
139+
if zetas is None:
140+
zetas = linspace(0, 0.9, 10)
141+
for zeta in zetas:
142+
# Calculate in polar coordinates
143+
factor = zeta/sqrt(1-zeta**2)
144+
x = linspace(0, sqrt(1-zeta**2),200)
145+
ang = pi*x
146+
mag = exp(-pi*factor*x)
147+
# Draw upper part in retangular coordinates
148+
xret = mag*cos(ang)
149+
yret = mag*sin(ang)
150+
ax.plot(xret,yret, 'k:', lw=1)
151+
# Draw lower part in retangular coordinates
152+
xret = mag*cos(-ang)
153+
yret = mag*sin(-ang)
154+
ax.plot(xret,yret,'k:', lw=1)
155+
# Annotation
156+
an_i = int(len(xret)/2.5)
157+
an_x = xret[an_i]
158+
an_y = yret[an_i]
159+
ax.annotate(str(round(zeta,2)), xy=(an_x, an_y), xytext=(an_x, an_y), size=7)
160+
161+
# Constant natural frequency lines
162+
if wns is None:
163+
wns = linspace(0, 1, 10)
164+
for a in wns:
165+
# Calculate in polar coordinates
166+
x = linspace(-pi/2,pi/2,200)
167+
ang = pi*a*sin(x)
168+
mag = exp(-pi*a*cos(x))
169+
# Draw in retangular coordinates
170+
xret = mag*cos(ang)
171+
yret = mag*sin(ang)
172+
ax.plot(xret,yret,'k:', lw=1)
173+
# Annotation
174+
an_i = -1
175+
an_x = xret[an_i]
176+
an_y = yret[an_i]
177+
num = '{:1.1f}'.format(a)
178+
ax.annotate("$\\frac{"+num+"\pi}{T}$", xy=(an_x, an_y), xytext=(an_x, an_y), size=9)
179+
180+
_final_setup(ax)
181+
return ax, fig
182+

‎control/pzmap.py‎

Lines changed: 20 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -40,15 +40,17 @@
4040
#
4141
# $Id:pzmap.py 819 2009-05-29 21:28:07Z murray $
4242

43-
from numpy import real, imag
44-
from .lti import LTI
43+
from numpy import real, imag, linspace, exp, cos, sin, sqrt
44+
from math import pi
45+
from .lti import LTI, isdtime, isctime
46+
from .grid import sgrid, zgrid, nogrid
4547

4648
__all__ = ['pzmap']
4749

4850
# TODO: Implement more elegant cross-style axes. See:
4951
# http://matplotlib.sourceforge.net/examples/axes_grid/demo_axisline_style.html
5052
# http://matplotlib.sourceforge.net/examples/axes_grid/demo_curvelinear_grid.html
51-
def pzmap(sys, Plot=True, title='Pole Zero Map'):
53+
def pzmap(sys, Plot=True, grid=False, title='Pole Zero Map'):
5254
"""
5355
Plot a pole/zero map for a linear system.
5456
@@ -59,6 +61,8 @@ def pzmap(sys, Plot=True, title='Pole Zero Map'):
5961
Plot: bool
6062
If ``True`` a graph is generated with Matplotlib,
6163
otherwise the poles and zeros are only computed and returned.
64+
grid: boolean (default = False)
65+
If True plot omega-damping grid.
6266
6367
Returns
6468
-------
@@ -75,18 +79,22 @@ def pzmap(sys, Plot=True, title='Pole Zero Map'):
7579

7680
if (Plot):
7781
import matplotlib.pyplot as plt
82+
83+
if grid:
84+
if isdtime(sys, strict=True):
85+
ax, fig = zgrid()
86+
else:
87+
ax, fig = sgrid()
88+
else:
89+
ax, fig = nogrid()
90+
7891
# Plot the locations of the poles and zeros
7992
if len(poles) > 0:
80-
plt.scatter(real(poles), imag(poles), s=50, marker='x')
93+
ax.scatter(real(poles), imag(poles), s=50, marker='x', facecolors='k')
8194
if len(zeros) > 0:
82-
plt.scatter(real(zeros), imag(zeros), s=50, marker='o',
83-
facecolors='none')
84-
# Add axes
85-
#Somewhat silly workaround
86-
plt.axhline(y=0, color='black')
87-
plt.axvline(x=0, color='black')
88-
plt.xlabel('Re')
89-
plt.ylabel('Im')
95+
ax.scatter(real(zeros), imag(zeros), s=50, marker='o',
96+
facecolors='none', edgecolors='k')
97+
9098

9199
plt.title(title)
92100

‎control/rlocus.py‎

Lines changed: 14 additions & 93 deletions
Original file line numberDiff line numberDiff line change
@@ -47,12 +47,15 @@
4747

4848
# Packages used by this module
4949
import numpy as np
50-
from scipy import array, poly1d, row_stack, zeros_like, real, imag
50+
from scipy import array, poly1d, row_stack, zeros_like, real, imag, exp, sin, cos, linspace, sqrt
51+
from math import pi
5152
import scipy.signal # signal processing toolbox
5253
import pylab # plotting routines
5354
from .xferfcn import _convertToTransferFunction
5455
from .exception import ControlMIMONotImplemented
5556
from functools import partial
57+
from .lti import isdtime
58+
from .grid import sgrid, zgrid, nogrid
5659

5760
__all__ = ['root_locus', 'rlocus']
5861

@@ -82,7 +85,7 @@ def root_locus(sys, kvect=None, xlim=None, ylim=None, plotstr='-', Plot=True,
8285
If True, report mouse clicks when close to the root-locus branches,
8386
calculate gain, damping and print
8487
grid: boolean (default = False)
85-
If True plot s-plane grid.
88+
If True plot omega-damping grid.
8689
8790
Returns
8891
-------
@@ -110,12 +113,18 @@ def root_locus(sys, kvect=None, xlim=None, ylim=None, plotstr='-', Plot=True,
110113
while new_figure_name in figure_title:
111114
new_figure_name = "Root Locus " + str(rloc_num)
112115
rloc_num += 1
113-
f = pylab.figure(new_figure_name)
116+
if grid:
117+
if isdtime(sys, strict=True):
118+
ax, f = zgrid()
119+
else:
120+
ax, f = sgrid()
121+
else:
122+
ax, f = nogrid()
123+
pylab.title(new_figure_name)
114124

115125
if PrintGain:
116126
f.canvas.mpl_connect(
117127
'button_release_event', partial(_RLFeedbackClicks, sys=sys))
118-
ax = pylab.axes()
119128

120129
# plot open loop poles
121130
poles = array(denp.r)
@@ -128,17 +137,13 @@ def root_locus(sys, kvect=None, xlim=None, ylim=None, plotstr='-', Plot=True,
128137

129138
# Now plot the loci
130139
for col in mymat.T:
131-
ax.plot(real(col), imag(col), plotstr)
140+
ax.plot(real(col), imag(col), plotstr, lw=3)
132141

133142
# Set up plot axes and labels
134143
if xlim:
135144
ax.set_xlim(xlim)
136145
if ylim:
137146
ax.set_ylim(ylim)
138-
ax.set_xlabel('Real')
139-
ax.set_ylabel('Imaginary')
140-
if grid:
141-
_sgrid_func()
142147
return mymat, kvect
143148

144149

@@ -350,88 +355,4 @@ def _RLFeedbackClicks(event, sys):
350355
print("Clicked at %10.4g%+10.4gj gain %10.4g damp %10.4g" %
351356
(s.real, s.imag, K.real, -1 * s.real / abs(s)))
352357

353-
354-
def _sgrid_func(fig=None, zeta=None, wn=None):
355-
if fig is None:
356-
fig = pylab.gcf()
357-
ax = fig.gca()
358-
xlocator = ax.get_xaxis().get_major_locator()
359-
360-
ylim = ax.get_ylim()
361-
ytext_pos_lim = ylim[1] - (ylim[1] - ylim[0]) * 0.03
362-
xlim = ax.get_xlim()
363-
xtext_pos_lim = xlim[0] + (xlim[1] - xlim[0]) * 0.0
364-
365-
if zeta is None:
366-
zeta = _default_zetas(xlim, ylim)
367-
368-
angules = []
369-
for z in zeta:
370-
if (z >= 1e-4) and (z <= 1):
371-
angules.append(np.pi/2 + np.arcsin(z))
372-
else:
373-
zeta.remove(z)
374-
y_over_x = np.tan(angules)
375-
376-
# zeta-constant lines
377-
378-
index = 0
379-
380-
for yp in y_over_x:
381-
ax.plot([0, xlocator()[0]], [0, yp*xlocator()[0]], color='gray',
382-
linestyle='dashed', linewidth=0.5)
383-
ax.plot([0, xlocator()[0]], [0, -yp * xlocator()[0]], color='gray',
384-
linestyle='dashed', linewidth=0.5)
385-
an = "%.2f" % zeta[index]
386-
if yp < 0:
387-
xtext_pos = 1/yp * ylim[1]
388-
ytext_pos = yp * xtext_pos_lim
389-
if np.abs(xtext_pos) > np.abs(xtext_pos_lim):
390-
xtext_pos = xtext_pos_lim
391-
else:
392-
ytext_pos = ytext_pos_lim
393-
ax.annotate(an, textcoords='data', xy=[xtext_pos, ytext_pos], fontsize=8)
394-
index += 1
395-
ax.plot([0, 0], [ylim[0], ylim[1]], color='gray', linestyle='dashed', linewidth=0.5)
396-
397-
angules = np.linspace(-90, 90, 20)*np.pi/180
398-
if wn is None:
399-
wn = _default_wn(xlocator(), ylim)
400-
401-
for om in wn:
402-
if om < 0:
403-
yp = np.sin(angules)*np.abs(om)
404-
xp = -np.cos(angules)*np.abs(om)
405-
ax.plot(xp, yp, color='gray',
406-
linestyle='dashed', linewidth=0.5)
407-
an = "%.2f" % -om
408-
ax.annotate(an, textcoords='data', xy=[om, 0], fontsize=8)
409-
410-
411-
def _default_zetas(xlim, ylim):
412-
"""Return default list of dumps coefficients"""
413-
sep1 = -xlim[0]/4
414-
ang1 = [np.arctan((sep1*i)/ylim[1]) for i in np.arange(1, 4, 1)]
415-
sep2 = ylim[1] / 3
416-
ang2 = [np.arctan(-xlim[0]/(ylim[1]-sep2*i)) for i in np.arange(1, 3, 1)]
417-
418-
angules = np.concatenate((ang1, ang2))
419-
angules = np.insert(angules, len(angules), np.pi/2)
420-
zeta = np.sin(angules)
421-
return zeta.tolist()
422-
423-
424-
def _default_wn(xloc, ylim):
425-
"""Return default wn for root locus plot"""
426-
427-
wn = xloc
428-
sep = xloc[1]-xloc[0]
429-
while np.abs(wn[0]) < ylim[1]:
430-
wn = np.insert(wn, 0, wn[0]-sep)
431-
432-
while len(wn) > 7:
433-
wn = wn[0:-1:2]
434-
435-
return wn
436-
437358
rlocus = root_locus

0 commit comments

Comments
 (0)