Skip to content

Commit 1ecf93a

Browse files
murrayrmmurrayrm
committed
* Moved nichols() code to separate file, added unit test; add'l updates pending
* Moved slycot imports inside functions that need them * Updated pvtol-nested example to work with v0.4 * See ChangeLog for detailed listing of changes
1 parent 5c72d5e commit 1ecf93a

9 files changed

Lines changed: 341 additions & 224 deletions

File tree

‎ChangeLog‎

Lines changed: 24 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,3 +1,27 @@
1+
2011-03-31 Richard Murray <murray@malabar.local>
2+
3+
* examples/pvtol-nested.py: updated stability margin plot to use
4+
proper calling format for bode().
5+
6+
* src/statesp.py (_convertToStateSpace): moved slycot import
7+
to the location where it is actually needed (allows running some
8+
commands without slycot installed)
9+
10+
* src/xferfcn.py (_convertToTransferFunction): moved slycot import
11+
to the location where it is actually needed (allows running some
12+
commands without slycot installed)
13+
14+
* src/nichols.py: new file for Nichols plot routines; move
15+
nichols(), nichols_grid(), closed_loop_contours(), m_circles(),
16+
n_circles()
17+
18+
* src/__init__.py, src/freqresp.py, src/matlab.py: updated to match
19+
new file structure for Nichols charts
20+
21+
* src/nichols.py (nichols): updated processing of freqresp to take
22+
into account the fact that return arguments are now a matrix of
23+
results (even for a SISO system)
24+
125
2011-03-30 Richard Murray <murray@malabar.local>
226

327
* tests/: added top level subdirectory, to be used for unit tests.

‎examples/pvtol-nested.py‎

Lines changed: 7 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -84,22 +84,22 @@
8484
# (gm, pm, wgc, wpc) = margin(L);
8585

8686
#! TODO: this figure has something wrong; axis limits mismatch
87-
figure(6); clf; subplot(221);
88-
(magh, phaseh) = bode(L);
87+
figure(6); clf;
88+
bode(L);
8989

9090
# Add crossover line
91-
subplot(magh); hold(True);
92-
loglog([10^-4, 10^3], [1, 1], 'k-')
91+
subplot(211); hold(True);
92+
loglog([1e-4, 1e3], [1, 1], 'k-')
9393

9494
# Replot phase starting at -90 degrees
9595
bode(L, logspace(-4, 3));
9696
(mag, phase, w) = freqresp(L, logspace(-4, 3));
9797
phase = phase - 360;
98-
subplot(phaseh);
99-
semilogx([10^-4, 10^3], [-180, -180], 'k-')
98+
subplot(212);
99+
semilogx([1e-4, 1e3], [-180, -180], 'k-')
100100
hold(True);
101101
semilogx(w, np.squeeze(phase), 'b-')
102-
axis([10^-4, 10^3, -360, 0]);
102+
axis([1e-4, 1e3, -360, 0]);
103103
xlabel('Frequency [deg]'); ylabel('Phase [deg]');
104104
# set(gca, 'YTick', [-360, -270, -180, -90, 0]);
105105
# set(gca, 'XTick', [10^-4, 10^-2, 1, 100]);

‎src/__init__.py‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -59,6 +59,7 @@
5959
from xferfcn import *
6060
from statesp import *
6161
from freqplot import *
62+
from nichols import *
6263
from bdalg import *
6364
from statefbk import *
6465
from delay import *

‎src/freqplot.py‎

Lines changed: 0 additions & 213 deletions
Original file line numberDiff line numberDiff line change
@@ -212,138 +212,6 @@ def nyquist(syslist, omega=None, Plot=True):
212212
plt.plot([-1], [0], 'r+')
213213

214214
return x, y, omega
215-
# Nichols plot
216-
# Contributed by Allan McInnes <Allan.McInnes@canterbury.ac.nz>
217-
#! TODO: need unit test code
218-
def nichols(syslist, omega=None):
219-
"""Nichols plot for a system
220-
221-
Usage
222-
=====
223-
magh = nichols(sys, omega=None)
224-
225-
Plots a Nichols plot for the system over a (optional) frequency range.
226-
227-
Parameters
228-
----------
229-
syslist : linsys
230-
List of linear input/output systems (single system is OK)
231-
omega : freq_range
232-
Range of frequencies (list or bounds) in rad/sec
233-
234-
Return values
235-
-------------
236-
None
237-
"""
238-
239-
# If argument was a singleton, turn it into a list
240-
if (not getattr(syslist, '__iter__', False)):
241-
syslist = (syslist,)
242-
243-
# Select a default range if none is provided
244-
if (omega == None):
245-
omega = default_frequency_range(syslist)
246-
247-
for sys in syslist:
248-
# Get the magnitude and phase of the system
249-
mag, phase, omega = sys.freqresp(omega)
250-
251-
# Convert to Nichols-plot format (phase in degrees,
252-
# and magnitude in dB)
253-
x = unwrap(sp.degrees(phase), 360)
254-
y = 20*sp.log10(mag)
255-
256-
# Generate the plot
257-
plt.plot(x, y)
258-
259-
plt.xlabel('Phase (deg)')
260-
plt.ylabel('Magnitude (dB)')
261-
plt.title('Nichols Plot')
262-
263-
# Mark the -180 point
264-
plt.plot([-180], [0], 'r+')
265-
266-
# Nichols grid
267-
def nichols_grid():
268-
"""Nichols chart grid
269-
270-
Usage
271-
=====
272-
nichols_grid()
273-
274-
Plots a Nichols chart grid on the current axis.
275-
276-
Parameters
277-
----------
278-
None
279-
280-
Return values
281-
-------------
282-
None
283-
"""
284-
mag_min_default = -40.0 # dB
285-
mag_step = 20.0 # dB
286-
287-
# Chart defaults
288-
phase_min, phase_max, mag_min, mag_max = -360.0, 0.0, mag_min_default, 40.0
289-
290-
# Set actual chart bounds based on current plot
291-
if plt.gcf().gca().has_data():
292-
phase_min, phase_max, mag_min, mag_max = plt.axis()
293-
294-
# M-circle magnitudes.
295-
# The "fixed" set are always generated, since this guarantees a recognizable
296-
# Nichols chart grid.
297-
mags_fixed = np.array([-40.0, -20.0, -12.0, -6.0, -3.0, -1.0, -0.5, 0.0,
298-
0.25, 0.5, 1.0, 3.0, 6.0, 12.0])
299-
300-
if mag_min < mag_min_default:
301-
# Outside the "fixed" set of magnitudes, the generated M-circles
302-
# are extended in steps of 'mag_step' dB to cover anything made
303-
# visible by the range of the existing plot
304-
mags_adjust = np.arange(mag_step*np.floor(mag_min/mag_step),
305-
mag_min_default, mag_step)
306-
mags = np.concatenate((mags_adjust, mags_fixed))
307-
else:
308-
mags = mags_fixed
309-
310-
# N-circle phases (should be in the range -360 to 0)
311-
phases = np.array([-0.25, -10.0, -20.0, -30.0, -45.0, -60.0, -90.0,
312-
-120.0, -150.0, -180.0, -210.0, -240.0, -270.0,
313-
-310.0, -325.0, -340.0, -350.0, -359.75])
314-
315-
# Find the M-contours
316-
m = m_circles(mags, phase_min=np.min(phases), phase_max=np.max(phases))
317-
m_mag = 20*sp.log10(np.abs(m))
318-
m_phase = sp.mod(sp.degrees(sp.angle(m)), -360.0) # Unwrap
319-
320-
# Find the N-contours
321-
n = n_circles(phases, mag_min=np.min(mags), mag_max=np.max(mags))
322-
n_mag = 20*sp.log10(np.abs(n))
323-
n_phase = sp.mod(sp.degrees(sp.angle(n)), -360.0) # Unwrap
324-
325-
# Plot the contours behind other plot elements.
326-
# The "phase offset" is used to produce copies of the chart that cover
327-
# the entire range of the plotted data, starting from a base chart computed
328-
# over the range -360 < phase < 0 (see above). Given the range
329-
# the base chart is computed over, the phase offset should be 0
330-
# for -360 < phase_min < 0.
331-
phase_offset_min = 360.0*np.ceil(phase_min/360.0)
332-
phase_offset_max = 360.0*np.ceil(phase_max/360.0) + 360.0
333-
phase_offsets = np.arange(phase_offset_min, phase_offset_max, 360.0)
334-
for phase_offset in phase_offsets:
335-
plt.plot(m_phase + phase_offset, m_mag, color='gray',
336-
linestyle='dashed', zorder=0)
337-
plt.plot(n_phase + phase_offset, n_mag, color='gray',
338-
linestyle='dashed', zorder=0)
339-
340-
# Add magnitude labels
341-
for x, y, m in zip(m_phase[:][-1], m_mag[:][-1], mags):
342-
align = 'right' if m < 0.0 else 'left'
343-
plt.text(x, y, str(m) + ' dB', size='small', ha=align)
344-
345-
# Make sure axes conform to any pre-existing plot.
346-
plt.axis([phase_min, phase_max, mag_min, mag_max])
347215

348216
# Gang of Four
349217
#! TODO: think about how (and whether) to handle lists of systems
@@ -462,84 +330,3 @@ def default_frequency_range(syslist):
462330

463331
return omega
464332

465-
# Compute contours of a closed-loop transfer function
466-
def closed_loop_contours(Hmag, Hphase):
467-
"""Contours of the function H = G/(1+G).
468-
469-
Usage
470-
=====
471-
contours = closed_loop_contours(mags, phases)
472-
473-
Parameters
474-
----------
475-
mags : array-like
476-
Meshgrid array of magnitudes of the contours
477-
phases : array-like
478-
Meshgrid array of phases in radians of the contours
479-
480-
Return values
481-
-------------
482-
contours : complex array
483-
Array of complex numbers corresponding to the contours.
484-
"""
485-
# Compute the contours in H-space
486-
H = Hmag*sp.exp(1.j*Hphase)
487-
488-
# Invert H = G/(1+G) to get an expression for the contours in G-space
489-
return H/(1.0 - H)
490-
491-
# M-circle
492-
def m_circles(mags, phase_min=-359.75, phase_max=-0.25):
493-
"""Constant-magnitude contours of the function H = G/(1+G).
494-
495-
Usage
496-
=====
497-
contours = m_circles(mags)
498-
499-
Parameters
500-
----------
501-
mags : array-like
502-
Array of magnitudes in dB of the M-circles
503-
phase_min : degrees
504-
Minimum phase in degrees of the N-circles
505-
phase_max : degrees
506-
Maximum phase in degrees of the N-circles
507-
508-
Return values
509-
-------------
510-
contours : complex array
511-
Array of complex numbers corresponding to the contours.
512-
"""
513-
# Convert magnitudes and phase range into a grid suitable for
514-
# building contours
515-
phases = sp.radians(sp.linspace(phase_min, phase_max, 500))
516-
Hmag, Hphase = sp.meshgrid(10.0**(mags/20.0), phases)
517-
return closed_loop_contours(Hmag, Hphase)
518-
519-
# N-circle
520-
def n_circles(phases, mag_min=-40.0, mag_max=12.0):
521-
"""Constant-phase contours of the function H = G/(1+G).
522-
523-
Usage
524-
=====
525-
contour = n_circles(angles)
526-
527-
Parameters
528-
----------
529-
phases : array-like
530-
Array of phases in degrees of the N-circles
531-
mag_min : dB
532-
Minimum magnitude in dB of the N-circles
533-
mag_max : dB
534-
Maximum magnitude in dB of the N-circles
535-
536-
Return values
537-
-------------
538-
contours : complex array
539-
Array of complex numbers corresponding to the contours.
540-
"""
541-
# Convert phases and magnitude range into a grid suitable for
542-
# building contours
543-
mags = sp.linspace(10**(mag_min/20.0), 10**(mag_max/20.0), 2000)
544-
Hphase, Hmag = sp.meshgrid(sp.radians(phases), mags)
545-
return closed_loop_contours(Hmag, Hphase)

‎src/matlab.py‎

Lines changed: 3 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -74,7 +74,8 @@
7474

7575
# Import MATLAB-like functions that can be used as-is
7676
from ctrlutil import unwrap
77-
from freqplot import nyquist, nichols, gangof4
77+
from freqplot import nyquist, gangof4
78+
from nichols import nichols, nichols_grid
7879
from bdalg import series, parallel, negate, feedback
7980
from pzmap import pzmap
8081
from statefbk import ctrb, obsv, gram, place, lqr
@@ -752,7 +753,7 @@ def ngrid():
752753
=====
753754
ngrid()
754755
"""
755-
freqplot.nichols_grid()
756+
nichols_grid()
756757

757758
#
758759
# Modifications to scipy.signal functions

0 commit comments

Comments
 (0)