Skip to content

Commit 43d8b95

Browse files
murrayrmmurrayrm
authored
Merge branch 'master' into sisotool_final
2 parents 3e7cad6 + 1254ccb commit 43d8b95

31 files changed

+11676
-9826
lines changed

‎.travis.yml‎

Lines changed: 18 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -51,8 +51,10 @@ before_install:
5151
- conda create -q -n test-environment python="$TRAVIS_PYTHON_VERSION" pip coverage
5252
- source activate test-environment
5353
# Install openblas if slycot is being used
54+
# also install scikit-build for the build process
5455
- if [[ "$SLYCOT" != "" ]]; then
55-
conda install openblas;
56+
conda install openblas;
57+
conda install -c conda-forge scikit-build;
5658
fi
5759
# Make sure to look in the right place for python libraries (for slycot)
5860
- export LIBRARY_PATH="$HOME/miniconda/envs/test-environment/lib"
@@ -65,16 +67,30 @@ install:
6567
- conda install $SCIPY matplotlib
6668
# Build slycot from source
6769
# For python 3, need to provide pointer to python library
70+
# Use "Unix Makefiles" as generator, because Ninja cannot handle Fortran
6871
#! git clone https://github.com/repagh/Slycot.git slycot;
6972
- if [[ "$SLYCOT" != "" ]]; then
7073
git clone https://github.com/python-control/Slycot.git slycot;
71-
cd slycot; python setup.py install; cd ..;
74+
cd slycot; python setup.py install -G "Unix Makefiles"; cd ..;
7275
fi
7376

7477
# command to run tests
7578
script:
7679
- 'if [ $SLYCOT != "" ]; then python -c "import slycot"; fi'
7780
- coverage run setup.py test
7881

82+
# only run examples if Slycot is install
83+
# set PYTHONPATH for examples
84+
# pmw needed for examples/tfvis.py
85+
# future is needed for Python 2, also for examples/tfvis.py
86+
87+
- if [[ "$SLYCOT" != "" ]]; then
88+
export PYTHONPATH=$PWD;
89+
conda install -c conda-forge pmw future;
90+
cd examples; bash run_examples.sh; cd ..;
91+
fi
92+
93+
# arbitrary change to try to trigger travis build
94+
7995
after_success:
8096
- coveralls

‎control/bdalg.py‎

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -239,14 +239,14 @@ def feedback(sys1, sys2=1, sign=-1):
239239
# its feedback member function.
240240
if isinstance(sys1, (int, float, complex, np.number)):
241241
if isinstance(sys2, tf.TransferFunction):
242-
sys1 = tf._convertToTransferFunction(sys1)
242+
sys1 = tf._convert_to_transfer_function(sys1)
243243
elif isinstance(sys2, ss.StateSpace):
244244
sys1 = ss._convertToStateSpace(sys1)
245245
elif isinstance(sys2, frd.FRD):
246246
sys1 = ss._convertToFRD(sys1)
247247
else: # sys2 is a scalar.
248-
sys1 = tf._convertToTransferFunction(sys1)
249-
sys2 = tf._convertToTransferFunction(sys2)
248+
sys1 = tf._convert_to_transfer_function(sys1)
249+
sys2 = tf._convert_to_transfer_function(sys2)
250250

251251
return sys1.feedback(sys2, sign)
252252

‎control/freqplot.py‎

Lines changed: 160 additions & 106 deletions
Large diffs are not rendered by default.

‎control/margins.py‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -144,7 +144,7 @@ def stability_margins(sysdata, returnall=False, epsw=0.0):
144144
sys = frdata.FRD(mag * np.exp(1j * phase * math.pi/180),
145145
omega, smooth=True)
146146
else:
147-
sys = xferfcn._convertToTransferFunction(sysdata)
147+
sys = xferfcn._convert_to_transfer_function(sysdata)
148148
except Exception as e:
149149
print (e)
150150
raise ValueError("Margin sysdata must be either a linear system or "
@@ -299,7 +299,7 @@ def phase_crossover_frequencies(sys):
299299
"""
300300

301301
# Convert to a transfer function
302-
tf = xferfcn._convertToTransferFunction(sys)
302+
tf = xferfcn._convert_to_transfer_function(sys)
303303

304304
# if not siso, fall back to (0,0) element
305305
#! TODO: should add a check and warning here

‎control/mateqn.py‎

Lines changed: 28 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -411,7 +411,7 @@ def dlyap(A,Q,C=None,E=None):
411411

412412
#### Riccati equation solvers care and dare
413413

414-
def care(A,B,Q,R=None,S=None,E=None):
414+
def care(A, B, Q, R=None, S=None, E=None, stabilizing=True):
415415
""" (X,L,G) = care(A,B,Q,R=None) solves the continuous-time algebraic Riccati
416416
equation
417417
@@ -527,7 +527,11 @@ def care(A,B,Q,R=None,S=None,E=None):
527527
raise e
528528

529529
try:
530-
X,rcond,w,S_o,U,A_inv = sb02md(n,A,G,Q,'C')
530+
if stabilizing:
531+
sort = 'S'
532+
else:
533+
sort = 'U'
534+
X, rcond, w, S_o, U, A_inv = sb02md(n, A, G, Q, 'C', sort=sort)
531535
except ValueError as ve:
532536
if ve.info < 0 or ve.info > 5:
533537
e = ValueError(ve.message)
@@ -613,8 +617,12 @@ def care(A,B,Q,R=None,S=None,E=None):
613617
# Solve the generalized algebraic Riccati equation by calling the
614618
# Slycot function sg02ad
615619
try:
616-
rcondu,X,alfar,alfai,beta,S_o,T,U,iwarn = \
617-
sg02ad('C','B','N','U','N','N','S','R',n,m,0,A,E,B,Q,R,S)
620+
if stabilizing:
621+
sort = 'S'
622+
else:
623+
sort = 'U'
624+
rcondu, X, alfar, alfai, beta, S_o, T, U, iwarn = \
625+
sg02ad('C', 'B', 'N', 'U', 'N', 'N', sort, 'R', n, m, 0, A, E, B, Q, R, S)
618626
except ValueError as ve:
619627
if ve.info < 0 or ve.info > 7:
620628
e = ValueError(ve.message)
@@ -671,7 +679,7 @@ def care(A,B,Q,R=None,S=None,E=None):
671679
else:
672680
raise ControlArgument("Invalid set of input parameters.")
673681

674-
def dare(A,B,Q,R,S=None,E=None):
682+
def dare(A, B, Q, R, S=None, E=None, stabilizing=True):
675683
""" (X,L,G) = dare(A,B,Q,R) solves the discrete-time algebraic Riccati
676684
equation
677685
@@ -692,8 +700,8 @@ def dare(A,B,Q,R,S=None,E=None):
692700
matrix :math:`G = (B^T X B + R)^{-1} (B^T X A + S^T)` and the closed loop
693701
eigenvalues L, i.e., the eigenvalues of A - B G , E.
694702
"""
695-
if S is not None or E is not None:
696-
return dare_old(A, B, Q, R, S, E)
703+
if S is not None or E is not None or not stabilizing:
704+
return dare_old(A, B, Q, R, S, E, stabilizing)
697705
else:
698706
Rmat = asmatrix(R)
699707
Qmat = asmatrix(Q)
@@ -702,7 +710,7 @@ def dare(A,B,Q,R,S=None,E=None):
702710
L = eigvals(A - B.dot(G))
703711
return X, L, G
704712

705-
def dare_old(A,B,Q,R,S=None,E=None):
713+
def dare_old(A, B, Q, R, S=None, E=None, stabilizing=True):
706714
# Make sure we can import required slycot routine
707715
try:
708716
from slycot import sb02md
@@ -795,7 +803,12 @@ def dare_old(A,B,Q,R,S=None,E=None):
795803
raise e
796804

797805
try:
798-
X,rcond,w,S,U,A_inv = sb02md(n,A,G,Q,'D')
806+
if stabilizing:
807+
sort = 'S'
808+
else:
809+
sort = 'U'
810+
811+
X, rcond, w, S, U, A_inv = sb02md(n, A, G, Q, 'D', sort=sort)
799812
except ValueError as ve:
800813
if ve.info < 0 or ve.info > 5:
801814
e = ValueError(ve.message)
@@ -884,8 +897,12 @@ def dare_old(A,B,Q,R,S=None,E=None):
884897
# Solve the generalized algebraic Riccati equation by calling the
885898
# Slycot function sg02ad
886899
try:
887-
rcondu,X,alfar,alfai,beta,S_o,T,U,iwarn = \
888-
sg02ad('D','B','N','U','N','N','S','R',n,m,0,A,E,B,Q,R,S)
900+
if stabilizing:
901+
sort = 'S'
902+
else:
903+
sort = 'U'
904+
rcondu, X, alfar, alfai, beta, S_o, T, U, iwarn = \
905+
sg02ad('D', 'B', 'N', 'U', 'N', 'N', sort, 'R', n, m, 0, A, E, B, Q, R, S)
889906
except ValueError as ve:
890907
if ve.info < 0 or ve.info > 7:
891908
e = ValueError(ve.message)

‎control/matlab/timeresp.py‎

Lines changed: 47 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -4,7 +4,7 @@
44
Note that the return arguments are different than in the standard control package.
55
"""
66

7-
__all__ = ['step', 'impulse', 'initial', 'lsim']
7+
__all__ = ['step', 'stepinfo', 'impulse', 'initial', 'lsim']
88

99
def step(sys, T=None, X0=0., input=0, output=None, return_x=False):
1010
'''
@@ -66,6 +66,52 @@ def step(sys, T=None, X0=0., input=0, output=None, return_x=False):
6666

6767
return yout, T
6868

69+
def stepinfo(sys, T=None, SettlingTimeThreshold=0.02, RiseTimeLimits=(0.1,0.9)):
70+
'''
71+
Step response characteristics (Rise time, Settling Time, Peak and others).
72+
73+
Parameters
74+
----------
75+
sys: StateSpace, or TransferFunction
76+
LTI system to simulate
77+
78+
T: array-like object, optional
79+
Time vector (argument is autocomputed if not given)
80+
81+
SettlingTimeThreshold: float value, optional
82+
Defines the error to compute settling time (default = 0.02)
83+
84+
RiseTimeLimits: tuple (lower_threshold, upper_theshold)
85+
Defines the lower and upper threshold for RiseTime computation
86+
87+
Returns
88+
-------
89+
S: a dictionary containing:
90+
RiseTime: Time from 10% to 90% of the steady-state value.
91+
SettlingTime: Time to enter inside a default error of 2%
92+
SettlingMin: Minimum value after RiseTime
93+
SettlingMax: Maximum value after RiseTime
94+
Overshoot: Percentage of the Peak relative to steady value
95+
Undershoot: Percentage of undershoot
96+
Peak: Absolute peak value
97+
PeakTime: time of the Peak
98+
SteadyStateValue: Steady-state value
99+
100+
101+
See Also
102+
--------
103+
step, lsim, initial, impulse
104+
105+
Examples
106+
--------
107+
>>> S = stepinfo(sys, T)
108+
'''
109+
from ..timeresp import step_info
110+
111+
S = step_info(sys, T, SettlingTimeThreshold, RiseTimeLimits)
112+
113+
return S
114+
69115
def impulse(sys, T=None, X0=0., input=0, output=None, return_x=False):
70116
'''
71117
Impulse response of a linear system

‎control/nichols.py‎

Lines changed: 34 additions & 22 deletions
Original file line numberDiff line numberDiff line change
@@ -1,3 +1,13 @@
1+
"""nichols.py
2+
3+
Functions for plotting Black-Nichols charts.
4+
5+
Routines in this module:
6+
7+
nichols.nichols_plot aliased as nichols.nichols
8+
nichols.nichols_grid
9+
"""
10+
111
# nichols.py - Nichols plot
212
#
313
# Contributed by Allan McInnes <Allan.McInnes@canterbury.ac.nz>
@@ -45,17 +55,17 @@
4555
from .ctrlutil import unwrap
4656
from .freqplot import default_frequency_range
4757

48-
__all__ = ['nichols_plot', 'nichols']
58+
__all__ = ['nichols_plot', 'nichols', 'nichols_grid']
59+
4960

50-
# Nichols plot
51-
def nichols_plot(syslist, omega=None, grid=True):
61+
def nichols_plot(sys_list, omega=None, grid=True):
5262
"""Nichols plot for a system
5363
5464
Plots a Nichols plot for the system over a (optional) frequency range.
5565
5666
Parameters
5767
----------
58-
syslist : list of LTI, or LTI
68+
sys_list : list of LTI, or LTI
5969
List of linear input/output systems (single system is OK)
6070
omega : array_like
6171
Range of frequencies (list or bounds) in rad/sec
@@ -68,14 +78,14 @@ def nichols_plot(syslist, omega=None, grid=True):
6878
"""
6979

7080
# If argument was a singleton, turn it into a list
71-
if (not getattr(syslist, '__iter__', False)):
72-
syslist = (syslist,)
81+
if not getattr(sys_list, '__iter__', False):
82+
sys_list = (sys_list,)
7383

7484
# Select a default range if none is provided
7585
if omega is None:
76-
omega = default_frequency_range(syslist)
86+
omega = default_frequency_range(sys_list)
7787

78-
for sys in syslist:
88+
for sys in sys_list:
7989
# Get the magnitude and phase of the system
8090
mag_tmp, phase_tmp, omega = sys.freqresp(omega)
8191
mag = np.squeeze(mag_tmp)
@@ -100,9 +110,8 @@ def nichols_plot(syslist, omega=None, grid=True):
100110
if grid:
101111
nichols_grid()
102112

103-
# Nichols grid
104-
#! TODO: Consider making linestyle configurable
105-
def nichols_grid(cl_mags=None, cl_phases=None):
113+
114+
def nichols_grid(cl_mags=None, cl_phases=None, line_style='dotted'):
106115
"""Nichols chart grid
107116
108117
Plots a Nichols chart grid on the current axis, or creates a new chart
@@ -116,6 +125,8 @@ def nichols_grid(cl_mags=None, cl_phases=None):
116125
cl_phases : array-like (degrees), optional
117126
Array of closed-loop phases defining the iso-phase lines on a custom
118127
Nichols chart. Must be in the range -360 < cl_phases < 0
128+
line_style : string, optional
129+
.. seealso:: https://matplotlib.org/gallery/lines_bars_and_markers/linestyles.html
119130
120131
Returns
121132
-------
@@ -137,12 +148,12 @@ def nichols_grid(cl_mags=None, cl_phases=None):
137148
# The key set of magnitudes are always generated, since this
138149
# guarantees a recognizable Nichols chart grid.
139150
key_cl_mags = np.array([-40.0, -20.0, -12.0, -6.0, -3.0, -1.0, -0.5, 0.0,
140-
0.25, 0.5, 1.0, 3.0, 6.0, 12.0])
151+
0.25, 0.5, 1.0, 3.0, 6.0, 12.0])
141152
# Extend the range of magnitudes if necessary. The extended arange
142153
# will end up empty if no extension is required. Assumes that closed-loop
143154
# magnitudes are approximately aligned with open-loop magnitudes beyond
144155
# the value of np.min(key_cl_mags)
145-
cl_mag_step = -20.0 # dB
156+
cl_mag_step = -20.0 # dB
146157
extended_cl_mags = np.arange(np.min(key_cl_mags),
147158
ol_mag_min + cl_mag_step, cl_mag_step)
148159
cl_mags = np.concatenate((extended_cl_mags, key_cl_mags))
@@ -163,12 +174,12 @@ def nichols_grid(cl_mags=None, cl_phases=None):
163174
# Find the M-contours
164175
m = m_circles(cl_mags, phase_min=np.min(cl_phases), phase_max=np.max(cl_phases))
165176
m_mag = 20*sp.log10(np.abs(m))
166-
m_phase = sp.mod(sp.degrees(sp.angle(m)), -360.0) # Unwrap
177+
m_phase = sp.mod(sp.degrees(sp.angle(m)), -360.0) # Unwrap
167178

168179
# Find the N-contours
169180
n = n_circles(cl_phases, mag_min=np.min(cl_mags), mag_max=np.max(cl_mags))
170181
n_mag = 20*sp.log10(np.abs(n))
171-
n_phase = sp.mod(sp.degrees(sp.angle(n)), -360.0) # Unwrap
182+
n_phase = sp.mod(sp.degrees(sp.angle(n)), -360.0) # Unwrap
172183

173184
# Plot the contours behind other plot elements.
174185
# The "phase offset" is used to produce copies of the chart that cover
@@ -182,10 +193,10 @@ def nichols_grid(cl_mags=None, cl_phases=None):
182193

183194
for phase_offset in phase_offsets:
184195
# Draw M and N contours
185-
plt.plot(m_phase + phase_offset, m_mag, color='gray',
186-
linestyle='dotted', zorder=0)
187-
plt.plot(n_phase + phase_offset, n_mag, color='gray',
188-
linestyle='dotted', zorder=0)
196+
plt.plot(m_phase + phase_offset, m_mag, color='lightgray',
197+
linestyle=line_style, zorder=0)
198+
plt.plot(n_phase + phase_offset, n_mag, color='lightgray',
199+
linestyle=line_style, zorder=0)
189200

190201
# Add magnitude labels
191202
for x, y, m in zip(m_phase[:][-1] + phase_offset, m_mag[:][-1], cl_mags):
@@ -203,7 +214,7 @@ def nichols_grid(cl_mags=None, cl_phases=None):
203214
# generating Nichols plots
204215
#
205216

206-
# Compute contours of a closed-loop transfer function
217+
207218
def closed_loop_contours(Gcl_mags, Gcl_phases):
208219
"""Contours of the function Gcl = Gol/(1+Gol), where
209220
Gol is an open-loop transfer function, and Gcl is a corresponding
@@ -229,7 +240,7 @@ def closed_loop_contours(Gcl_mags, Gcl_phases):
229240
# Invert Gcl = Gol/(1+Gol) to map the contours into the open-loop space
230241
return Gcl/(1.0 - Gcl)
231242

232-
# M-circle
243+
233244
def m_circles(mags, phase_min=-359.75, phase_max=-0.25):
234245
"""Constant-magnitude contours of the function Gcl = Gol/(1+Gol), where
235246
Gol is an open-loop transfer function, and Gcl is a corresponding
@@ -255,7 +266,7 @@ def m_circles(mags, phase_min=-359.75, phase_max=-0.25):
255266
Gcl_mags, Gcl_phases = sp.meshgrid(10.0**(mags/20.0), phases)
256267
return closed_loop_contours(Gcl_mags, Gcl_phases)
257268

258-
# N-circle
269+
259270
def n_circles(phases, mag_min=-40.0, mag_max=12.0):
260271
"""Constant-phase contours of the function Gcl = Gol/(1+Gol), where
261272
Gol is an open-loop transfer function, and Gcl is a corresponding
@@ -281,5 +292,6 @@ def n_circles(phases, mag_min=-40.0, mag_max=12.0):
281292
Gcl_phases, Gcl_mags = sp.meshgrid(sp.radians(phases), mags)
282293
return closed_loop_contours(Gcl_mags, Gcl_phases)
283294

295+
284296
# Function aliases
285297
nichols = nichols_plot

0 commit comments

Comments
 (0)