Skip to content

Commit b6b92d0

Browse files
murrayrmmurrayrm
authored
Merge branch 'master' into rlocus_zoom
2 parents c65c100 + 4286329 commit b6b92d0

39 files changed

+1421
-824
lines changed

‎.travis.yml

Lines changed: 5 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -28,10 +28,9 @@ env:
2828
before_install:
2929
# Install gfortran for testing slycot; use apt-get instead of conda in
3030
# order to include the proper CXXABI dependency (updated in GCC 4.9)
31-
# Also need to include liblapack here, to make sure paths are right
3231
- if [[ "$SLYCOT" != "" ]]; then
3332
sudo apt-get update -qq;
34-
sudo apt-get install gfortran liblapack-dev;
33+
sudo apt-get install gfortran;
3534
fi
3635
# Install display manager to allow testing of plotting functions
3736
- export DISPLAY=:99.0
@@ -51,6 +50,10 @@ before_install:
5150
- conda info -a
5251
- conda create -q -n test-environment python="$TRAVIS_PYTHON_VERSION" pip coverage
5352
- source activate test-environment
53+
# Install openblas if slycot is being used
54+
- if [[ "$SLYCOT" != "" ]]; then
55+
conda install openblas;
56+
fi
5457
# Make sure to look in the right place for python libraries (for slycot)
5558
- export LIBRARY_PATH="$HOME/miniconda/envs/test-environment/lib"
5659
# coveralls not in conda repos => install via pip instead

‎README.rst

Lines changed: 10 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -24,7 +24,7 @@ Features
2424
Links
2525
=====
2626

27-
- Project home page: http://python-control.sourceforge.net
27+
- Project home page: http://python-control.org
2828
- Source code repository: https://github.com/python-control/python-control
2929
- Documentation: http://python-control.readthedocs.org/
3030
- Issue tracker: https://github.com/python-control/python-control/issues
@@ -46,7 +46,7 @@ https://github.com/python-control/Slycot
4646
Installation
4747
============
4848

49-
The package may be installed using pip or distutils.
49+
The package may be installed using pip, conda, or distutils.
5050

5151
Pip
5252
---
@@ -56,6 +56,14 @@ To install using pip::
5656
pip install slycot # optional
5757
pip install control
5858

59+
conda-forge
60+
-----------
61+
62+
Binaries are available from conda-forge for selected platforms (Linux and
63+
MacOS). Install using
64+
65+
conda install -c conda-forge control
66+
5967
Distutils
6068
---------
6169

‎control/bdalg.py

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -62,12 +62,12 @@
6262
__all__ = ['series', 'parallel', 'negate', 'feedback', 'append', 'connect']
6363

6464
def series(sys1, *sysn):
65-
"""Return the series connection (... * sys3 *) sys2 * sys1
65+
"""Return the series connection (... \* sys3 \*) sys2 \* sys1
6666
6767
Parameters
6868
----------
6969
sys1: scalar, StateSpace, TransferFunction, or FRD
70-
*sysn: other scalars, StateSpaces, TransferFunctions, or FRDs
70+
sysn: other scalars, StateSpaces, TransferFunctions, or FRDs
7171
7272
Returns
7373
-------

‎control/canonical.py

Lines changed: 75 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -6,8 +6,8 @@
66
from .statesp import StateSpace
77
from .statefbk import ctrb, obsv
88

9-
from numpy import zeros, shape, poly
10-
from numpy.linalg import solve, matrix_rank
9+
from numpy import zeros, shape, poly, iscomplex, hstack
10+
from numpy.linalg import solve, matrix_rank, eig
1111

1212
__all__ = ['canonical_form', 'reachable_form', 'observable_form']
1313

@@ -22,7 +22,7 @@ def canonical_form(xsys, form='reachable'):
2222
Canonical form for transformation. Chosen from:
2323
* 'reachable' - reachable canonical form
2424
* 'observable' - observable canonical form
25-
* 'modal' - modal canonical form [not implemented]
25+
* 'modal' - modal canonical form
2626
2727
Returns
2828
-------
@@ -37,6 +37,8 @@ def canonical_form(xsys, form='reachable'):
3737
return reachable_form(xsys)
3838
elif form == 'observable':
3939
return observable_form(xsys)
40+
elif form == 'modal':
41+
return modal_form(xsys)
4042
else:
4143
raise ControlNotImplemented(
4244
"Canonical form '%s' not yet implemented" % form)
@@ -133,9 +135,78 @@ def observable_form(xsys):
133135
Wrz = obsv(zsys.A, zsys.C)
134136

135137
# Transformation from one form to another
136-
Tzx = inv(Wrz) * Wrx
138+
Tzx = solve(Wrz, Wrx) # matrix left division, Tzx = inv(Wrz) * Wrx
139+
140+
if matrix_rank(Tzx) != xsys.states:
141+
raise ValueError("Transformation matrix singular to working precision.")
137142

138143
# Finally, compute the output matrix
139144
zsys.B = Tzx * xsys.B
140145

141146
return zsys, Tzx
147+
148+
def modal_form(xsys):
149+
"""Convert a system into modal canonical form
150+
151+
Parameters
152+
----------
153+
xsys : StateSpace object
154+
System to be transformed, with state `x`
155+
156+
Returns
157+
-------
158+
zsys : StateSpace object
159+
System in modal canonical form, with state `z`
160+
T : matrix
161+
Coordinate transformation: z = T * x
162+
"""
163+
# Check to make sure we have a SISO system
164+
if not issiso(xsys):
165+
raise ControlNotImplemented(
166+
"Canonical forms for MIMO systems not yet supported")
167+
168+
# Create a new system, starting with a copy of the old one
169+
zsys = StateSpace(xsys)
170+
171+
# Calculate eigenvalues and matrix of eigenvectors Tzx,
172+
eigval, eigvec = eig(xsys.A)
173+
174+
# Eigenvalues and according eigenvectors are not sorted,
175+
# thus modal transformation is ambiguous
176+
# Sorting eigenvalues and respective vectors by largest to smallest eigenvalue
177+
idx = eigval.argsort()[::-1]
178+
eigval = eigval[idx]
179+
eigvec = eigvec[:,idx]
180+
181+
# If all eigenvalues are real, the matrix of eigenvectors is Tzx directly
182+
if not iscomplex(eigval).any():
183+
Tzx = eigvec
184+
else:
185+
# A is an arbitrary semisimple matrix
186+
187+
# Keep track of complex conjugates (need only one)
188+
lst_conjugates = []
189+
Tzx = None
190+
for val, vec in zip(eigval, eigvec.T):
191+
if iscomplex(val):
192+
if val not in lst_conjugates:
193+
lst_conjugates.append(val.conjugate())
194+
if Tzx is not None:
195+
Tzx = hstack((Tzx, hstack((vec.real.T, vec.imag.T))))
196+
else:
197+
Tzx = hstack((vec.real.T, vec.imag.T))
198+
else:
199+
# if conjugate has already been seen, skip this eigenvalue
200+
lst_conjugates.remove(val)
201+
else:
202+
if Tzx is not None:
203+
Tzx = hstack((Tzx, vec.real.T))
204+
else:
205+
Tzx = vec.real.T
206+
207+
# Generate the system matrices for the desired canonical form
208+
zsys.A = solve(Tzx, xsys.A).dot(Tzx)
209+
zsys.B = solve(Tzx, xsys.B)
210+
zsys.C = xsys.C.dot(Tzx)
211+
212+
return zsys, Tzx

‎control/config.py

Lines changed: 5 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -18,6 +18,8 @@
1818
def use_matlab_defaults():
1919
"""
2020
Use MATLAB compatible configuration settings
21+
22+
The following conventions are used:
2123
* Bode plots plot gain in dB, phase in degrees, frequency in Hertz
2224
"""
2325
# Bode plot defaults
@@ -28,7 +30,9 @@ def use_matlab_defaults():
2830
# Set defaults to match FBS (Astrom and Murray)
2931
def use_fbs_defaults():
3032
"""
31-
Use `Astrom and Murray <http://fbsbook.org>`_ compatible settings
33+
Use `Feedback Systems <http://fbsbook.org>`_ (FBS) compatible settings
34+
35+
The following conventions are used:
3236
* Bode plots plot gain in powers of ten, phase in degrees,
3337
frequency in Hertz
3438
"""

‎control/frdata.py

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -452,8 +452,8 @@ def _convertToFRD(sys, omega, inputs=1, outputs=1):
452452
scalar, then the number of inputs and outputs can be specified
453453
manually, as in:
454454
455-
>>> sys = _convertToFRD(3.) # Assumes inputs = outputs = 1
456-
>>> sys = _convertToFRD(1., inputs=3, outputs=2)
455+
>>> frd = _convertToFRD(3., omega) # Assumes inputs = outputs = 1
456+
>>> frd = _convertToFRD(1., omegs, inputs=3, outputs=2)
457457
458458
In the latter example, sys's matrix transfer function is [[1., 1., 1.]
459459
[1., 1., 1.]].

‎control/freqplot.py

Lines changed: 60 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -47,6 +47,7 @@
4747
import math
4848
from .ctrlutil import unwrap
4949
from .bdalg import feedback
50+
from .margins import stability_margins
5051

5152
__all__ = ['bode_plot', 'nyquist_plot', 'gangof4_plot',
5253
'bode', 'nyquist', 'gangof4']
@@ -60,17 +61,18 @@
6061

6162
# Bode plot
6263
def bode_plot(syslist, omega=None, dB=None, Hz=None, deg=None,
63-
Plot=True, omega_limits=None, omega_num=None, *args, **kwargs):
64-
"""Bode plot for a system
64+
Plot=True, omega_limits=None, omega_num=None,margins=None, *args, **kwargs):
65+
"""
66+
Bode plot for a system
6567
6668
Plots a Bode plot for the system over a (optional) frequency range.
6769
6870
Parameters
6971
----------
7072
syslist : linsys
7173
List of linear input/output systems (single system is OK)
72-
omega : freq_range
73-
Range of frequencies in rad/sec
74+
omega : list
75+
List of frequencies in rad/sec to be used for frequency response
7476
dB : boolean
7577
If True, plot result in dB
7678
Hz : boolean
@@ -84,7 +86,9 @@ def bode_plot(syslist, omega=None, dB=None, Hz=None, deg=None,
8486
If Hz=True the limits are in Hz otherwise in rad/s.
8587
omega_num: int
8688
number of samples
87-
*args, **kwargs:
89+
margins : boolean
90+
if True, plot gain and phase margin
91+
\*args, \**kwargs:
8892
Additional options to matplotlib (color, linestyle, etc)
8993
9094
Returns
@@ -105,7 +109,7 @@ def bode_plot(syslist, omega=None, dB=None, Hz=None, deg=None,
105109
2. If a discrete time model is given, the frequency response is plotted
106110
along the upper branch of the unit circle, using the mapping z = exp(j
107111
\omega dt) where omega ranges from 0 to pi/dt and dt is the discrete
108-
time base. If not timebase is specified (dt = True), dt is set to 1.
112+
timebase. If not timebase is specified (dt = True), dt is set to 1.
109113
110114
Examples
111115
--------
@@ -208,7 +212,7 @@ def bode_plot(syslist, omega=None, dB=None, Hz=None, deg=None,
208212
color=pltline[0].get_color())
209213

210214
# Add a grid to the plot + labeling
211-
ax_mag.grid(True, which='both')
215+
ax_mag.grid(False if margins else True, which='both')
212216
ax_mag.set_ylabel("Magnitude (dB)" if dB else "Magnitude")
213217

214218
# Phase plot
@@ -218,6 +222,50 @@ def bode_plot(syslist, omega=None, dB=None, Hz=None, deg=None,
218222
phase_plot = phase
219223
ax_phase.semilogx(omega_plot, phase_plot, *args, **kwargs)
220224

225+
# Show the phase and gain margins in the plot
226+
if margins:
227+
margin = stability_margins(sys)
228+
gm, pm, Wcg, Wcp = margin[0], margin[1], margin[3], margin[4]
229+
if pm >= 0.:
230+
phase_limit = -180.
231+
else:
232+
phase_limit = 180.
233+
234+
ax_mag.axhline(y=0 if dB else 1, color='k', linestyle=':')
235+
ax_phase.axhline(y=phase_limit if deg else math.radians(phase_limit), color='k', linestyle=':')
236+
mag_ylim = ax_mag.get_ylim()
237+
phase_ylim = ax_phase.get_ylim()
238+
239+
if pm != float('inf') and Wcp != float('nan'):
240+
if dB:
241+
ax_mag.semilogx([Wcp, Wcp], [0.,-1e5],color='k', linestyle=':')
242+
else:
243+
ax_mag.loglog([Wcp,Wcp], [1.,1e-8],color='k',linestyle=':')
244+
245+
if deg:
246+
ax_phase.semilogx([Wcp, Wcp], [1e5, phase_limit+pm],color='k', linestyle=':')
247+
ax_phase.semilogx([Wcp, Wcp], [phase_limit + pm, phase_limit],color='k')
248+
else:
249+
ax_phase.semilogx([Wcp, Wcp], [1e5, math.radians(phase_limit)+math.radians(pm)],color='k', linestyle=':')
250+
ax_phase.semilogx([Wcp, Wcp], [math.radians(phase_limit) +math.radians(pm), math.radians(phase_limit)],color='k')
251+
252+
if gm != float('inf') and Wcg != float('nan'):
253+
if dB:
254+
ax_mag.semilogx([Wcg, Wcg], [-20.*np.log10(gm), -1e5],color='k', linestyle=':')
255+
ax_mag.semilogx([Wcg, Wcg], [0,-20*np.log10(gm)],color='k')
256+
else:
257+
ax_mag.loglog([Wcg, Wcg], [1./gm,1e-8],color='k', linestyle=':')
258+
ax_mag.loglog([Wcg, Wcg], [1.,1./gm],color='k')
259+
260+
if deg:
261+
ax_phase.semilogx([Wcg, Wcg], [1e-8, phase_limit],color='k', linestyle=':')
262+
else:
263+
ax_phase.semilogx([Wcg, Wcg], [1e-8, math.radians(phase_limit)],color='k', linestyle=':')
264+
265+
ax_mag.set_ylim(mag_ylim)
266+
ax_phase.set_ylim(phase_ylim)
267+
plt.suptitle('Gm = %.2f %s(at %.2f rad/s), Pm = %.2f %s (at %.2f rad/s)'%(20*np.log10(gm) if dB else gm,'dB ' if dB else '\b',Wcg,pm if deg else math.radians(pm),'deg' if deg else 'rad',Wcp))
268+
221269
if nyquistfrq_plot:
222270
ax_phase.axvline(nyquistfrq_plot, color=pltline[0].get_color())
223271

@@ -236,7 +284,7 @@ def genZeroCenteredSeries(val_min, val_max, period):
236284
ylim = ax_phase.get_ylim()
237285
ax_phase.set_yticks(genZeroCenteredSeries(ylim[0], ylim[1], math.pi / 4.))
238286
ax_phase.set_yticks(genZeroCenteredSeries(ylim[0], ylim[1], math.pi / 12.), minor=True)
239-
ax_phase.grid(True, which='both')
287+
ax_phase.grid(False if margins else True, which='both')
240288
# ax_mag.grid(which='minor', alpha=0.3)
241289
# ax_mag.grid(which='major', alpha=0.9)
242290
# ax_phase.grid(which='minor', alpha=0.3)
@@ -253,7 +301,8 @@ def genZeroCenteredSeries(val_min, val_max, period):
253301
# Nyquist plot
254302
def nyquist_plot(syslist, omega=None, Plot=True, color='b',
255303
labelFreq=0, *args, **kwargs):
256-
"""Nyquist plot for a system
304+
"""
305+
Nyquist plot for a system
257306
258307
Plots a Nyquist plot for the system over a (optional) frequency range.
259308
@@ -267,7 +316,7 @@ def nyquist_plot(syslist, omega=None, Plot=True, color='b',
267316
If True, plot magnitude
268317
labelFreq : int
269318
Label every nth frequency on the plot
270-
*args, **kwargs:
319+
\*args, \**kwargs:
271320
Additional options to matplotlib (color, linestyle, etc)
272321
273322
Returns
@@ -283,6 +332,7 @@ def nyquist_plot(syslist, omega=None, Plot=True, color='b',
283332
--------
284333
>>> sys = ss("1. -2; 3. -4", "5.; 7", "6. 8", "9.")
285334
>>> real, imag, freq = nyquist_plot(sys)
335+
286336
"""
287337
# If argument was a singleton, turn it into a list
288338
if (not getattr(syslist, '__iter__', False)):

0 commit comments

Comments
 (0)