Skip to content

Commit 6b66e28

Browse files
bnavigatorbnavigator
authored
Merge branch 'master' into zero-poles-time-vector
2 parents ddef939 + 76ec2de commit 6b66e28

30 files changed

+250
-616
lines changed

‎.travis.yml

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -49,6 +49,11 @@ jobs:
4949
services: xvfb
5050
python: "3.8"
5151
env: SCIPY=scipy SLYCOT=source
52+
- name: "use numpy matrix"
53+
dist: xenial
54+
services: xvfb
55+
python: "3.8"
56+
env: SCIPY=scipy SLYCOT=source PYTHON_CONTROL_STATESPACE_ARRAY=1
5257

5358
# Exclude combinations that are very unlikely (and don't work)
5459
exclude:

‎README.rst

Lines changed: 9 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -99,10 +99,16 @@ You can check out the latest version of the source code with the command::
9999
Testing
100100
-------
101101

102-
You can run a set of unit tests to make sure that everything is working
103-
correctly. After installation, run::
102+
You can run the unit tests with `pytest`_ to make sure that everything is
103+
working correctly. Inside the source directory, run::
104104

105-
python setup.py test
105+
pytest -v
106+
107+
or to test the installed package::
108+
109+
pytest --pyargs control -v
110+
111+
.. _pytest: https://docs.pytest.org/
106112

107113
License
108114
-------

‎control/canonical.py

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

9-
from numpy import zeros, shape, poly, iscomplex, hstack, dot, transpose
9+
from numpy import zeros, zeros_like, shape, poly, iscomplex, vstack, hstack, dot, \
10+
transpose, empty
1011
from numpy.linalg import solve, matrix_rank, eig
1112

1213
__all__ = ['canonical_form', 'reachable_form', 'observable_form', 'modal_form',
@@ -70,9 +71,9 @@ def reachable_form(xsys):
7071
zsys = StateSpace(xsys)
7172

7273
# Generate the system matrices for the desired canonical form
73-
zsys.B = zeros(shape(xsys.B))
74+
zsys.B = zeros_like(xsys.B)
7475
zsys.B[0, 0] = 1.0
75-
zsys.A = zeros(shape(xsys.A))
76+
zsys.A = zeros_like(xsys.A)
7677
Apoly = poly(xsys.A) # characteristic polynomial
7778
for i in range(0, xsys.states):
7879
zsys.A[0, i] = -Apoly[i+1] / Apoly[0]
@@ -124,9 +125,9 @@ def observable_form(xsys):
124125
zsys = StateSpace(xsys)
125126

126127
# Generate the system matrices for the desired canonical form
127-
zsys.C = zeros(shape(xsys.C))
128+
zsys.C = zeros_like(xsys.C)
128129
zsys.C[0, 0] = 1
129-
zsys.A = zeros(shape(xsys.A))
130+
zsys.A = zeros_like(xsys.A)
130131
Apoly = poly(xsys.A) # characteristic polynomial
131132
for i in range(0, xsys.states):
132133
zsys.A[i, 0] = -Apoly[i+1] / Apoly[0]
@@ -144,7 +145,7 @@ def observable_form(xsys):
144145
raise ValueError("Transformation matrix singular to working precision.")
145146

146147
# Finally, compute the output matrix
147-
zsys.B = Tzx * xsys.B
148+
zsys.B = Tzx.dot(xsys.B)
148149

149150
return zsys, Tzx
150151

@@ -174,9 +175,9 @@ def modal_form(xsys):
174175
# Calculate eigenvalues and matrix of eigenvectors Tzx,
175176
eigval, eigvec = eig(xsys.A)
176177

177-
# Eigenvalues and according eigenvectors are not sorted,
178+
# Eigenvalues and corresponding eigenvectors are not sorted,
178179
# thus modal transformation is ambiguous
179-
# Sorting eigenvalues and respective vectors by largest to smallest eigenvalue
180+
# Sort eigenvalues and vectors from largest to smallest eigenvalue
180181
idx = eigval.argsort()[::-1]
181182
eigval = eigval[idx]
182183
eigvec = eigvec[:,idx]
@@ -189,23 +190,18 @@ def modal_form(xsys):
189190

190191
# Keep track of complex conjugates (need only one)
191192
lst_conjugates = []
192-
Tzx = None
193+
Tzx = empty((0, xsys.A.shape[0])) # empty zero-height row matrix
193194
for val, vec in zip(eigval, eigvec.T):
194195
if iscomplex(val):
195196
if val not in lst_conjugates:
196197
lst_conjugates.append(val.conjugate())
197-
if Tzx is not None:
198-
Tzx = hstack((Tzx, hstack((vec.real.T, vec.imag.T))))
199-
else:
200-
Tzx = hstack((vec.real.T, vec.imag.T))
198+
Tzx = vstack((Tzx, vec.real, vec.imag))
201199
else:
202200
# if conjugate has already been seen, skip this eigenvalue
203201
lst_conjugates.remove(val)
204202
else:
205-
if Tzx is not None:
206-
Tzx = hstack((Tzx, vec.real.T))
207-
else:
208-
Tzx = vec.real.T
203+
Tzx = vstack((Tzx, vec.real))
204+
Tzx = Tzx.T
209205

210206
# Generate the system matrices for the desired canonical form
211207
zsys.A = solve(Tzx, xsys.A).dot(Tzx)

‎control/freqplot.py

Lines changed: 24 additions & 21 deletions
Original file line numberDiff line numberDiff line change
@@ -40,11 +40,13 @@
4040
# SUCH DAMAGE.
4141
#
4242
# $Id$
43+
44+
import math
45+
4346
import matplotlib as mpl
4447
import matplotlib.pyplot as plt
45-
import scipy as sp
4648
import numpy as np
47-
import math
49+
4850
from .ctrlutil import unwrap
4951
from .bdalg import feedback
5052
from .margins import stability_margins
@@ -110,9 +112,9 @@ def bode_plot(syslist, omega=None,
110112
config.defaults['freqplot.number_of_samples'].
111113
margins : bool
112114
If True, plot gain and phase margin.
113-
*args : `matplotlib` plot positional properties, optional
115+
*args : :func:`matplotlib.pyplot.plot` positional properties, optional
114116
Additional arguments for `matplotlib` plots (color, linestyle, etc)
115-
**kwargs : `matplotlib` plot keyword properties, optional
117+
**kwargs : :func:`matplotlib.pyplot.plot` keyword properties, optional
116118
Additional keywords (passed to `matplotlib`)
117119
118120
Returns
@@ -128,21 +130,22 @@ def bode_plot(syslist, omega=None,
128130
----------------
129131
grid : bool
130132
If True, plot grid lines on gain and phase plots. Default is set by
131-
config.defaults['bode.grid'].
133+
`config.defaults['bode.grid']`.
134+
132135
133136
The default values for Bode plot configuration parameters can be reset
134137
using the `config.defaults` dictionary, with module name 'bode'.
135138
136139
Notes
137140
-----
138-
1. Alternatively, you may use the lower-level method (mag, phase, freq)
139-
= sys.freqresp(freq) to generate the frequency response for a system,
140-
but it returns a MIMO response.
141+
1. Alternatively, you may use the lower-level method
142+
``(mag, phase, freq) = sys.freqresp(freq)`` to generate the frequency
143+
response for a system, but it returns a MIMO response.
141144
142145
2. If a discrete time model is given, the frequency response is plotted
143-
along the upper branch of the unit circle, using the mapping z = exp(j
144-
\\omega dt) where omega ranges from 0 to pi/dt and dt is the discrete
145-
timebase. If not timebase is specified (dt = True), dt is set to 1.
146+
along the upper branch of the unit circle, using the mapping z = exp(j
147+
\\omega dt) where omega ranges from 0 to pi/dt and dt is the discrete
148+
timebase. If not timebase is specified (dt = True), dt is set to 1.
146149
147150
Examples
148151
--------
@@ -183,12 +186,12 @@ def bode_plot(syslist, omega=None,
183186
if Hz:
184187
omega_limits *= 2. * math.pi
185188
if omega_num:
186-
omega = sp.logspace(np.log10(omega_limits[0]),
189+
omega = np.logspace(np.log10(omega_limits[0]),
187190
np.log10(omega_limits[1]),
188191
num=omega_num,
189192
endpoint=True)
190193
else:
191-
omega = sp.logspace(np.log10(omega_limits[0]),
194+
omega = np.logspace(np.log10(omega_limits[0]),
192195
np.log10(omega_limits[1]),
193196
endpoint=True)
194197

@@ -464,9 +467,9 @@ def nyquist_plot(syslist, omega=None, plot=True, label_freq=0,
464467
Label every nth frequency on the plot
465468
arrowhead_width : arrow head width
466469
arrowhead_length : arrow head length
467-
*args : `matplotlib` plot positional properties, optional
470+
*args : :func:`matplotlib.pyplot.plot` positional properties, optional
468471
Additional arguments for `matplotlib` plots (color, linestyle, etc)
469-
**kwargs : `matplotlib` plot keyword properties, optional
472+
**kwargs : :func:`matplotlib.pyplot.plot` keyword properties, optional
470473
Additional keywords (passed to `matplotlib`)
471474
472475
Returns
@@ -529,8 +532,8 @@ def nyquist_plot(syslist, omega=None, plot=True, label_freq=0,
529532
phase = np.squeeze(phase_tmp)
530533

531534
# Compute the primary curve
532-
x = sp.multiply(mag, sp.cos(phase))
533-
y = sp.multiply(mag, sp.sin(phase))
535+
x = np.multiply(mag, np.cos(phase))
536+
y = np.multiply(mag, np.sin(phase))
534537

535538
if plot:
536539
# Plot the primary curve and mirror image
@@ -539,13 +542,13 @@ def nyquist_plot(syslist, omega=None, plot=True, label_freq=0,
539542
ax = plt.gca()
540543
# Plot arrow to indicate Nyquist encirclement orientation
541544
ax.arrow(x[0], y[0], (x[1]-x[0])/2, (y[1]-y[0])/2, fc=c, ec=c,
542-
head_width=arrowhead_width,
545+
head_width=arrowhead_width,
543546
head_length=arrowhead_length)
544547

545548
plt.plot(x, -y, '-', color=c, *args, **kwargs)
546549
ax.arrow(
547550
x[-1], -y[-1], (x[-1]-x[-2])/2, (y[-1]-y[-2])/2,
548-
fc=c, ec=c, head_width=arrowhead_width,
551+
fc=c, ec=c, head_width=arrowhead_width,
549552
head_length=arrowhead_length)
550553

551554
# Mark the -1 point
@@ -556,7 +559,7 @@ def nyquist_plot(syslist, omega=None, plot=True, label_freq=0,
556559
ind = slice(None, None, label_freq)
557560
for xpt, ypt, omegapt in zip(x[ind], y[ind], omega[ind]):
558561
# Convert to Hz
559-
f = omegapt / (2 * sp.pi)
562+
f = omegapt / (2 * np.pi)
560563

561564
# Factor out multiples of 1000 and limit the
562565
# result to the range [-8, 8].
@@ -601,7 +604,7 @@ def gangof4_plot(P, C, omega=None, **kwargs):
601604
Linear input/output systems (process and control)
602605
omega : array
603606
Range of frequencies (list or bounds) in rad/sec
604-
**kwargs : `matplotlib` plot keyword properties, optional
607+
**kwargs : :func:`matplotlib.pyplot.plot` keyword properties, optional
605608
Additional keywords (passed to `matplotlib`)
606609
607610
Returns

‎control/iosys.py

Lines changed: 6 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -1514,8 +1514,9 @@ def find_eqpt(sys, x0, u0=[], y0=None, t=0, params={},
15141514
return_y : bool, optional
15151515
If True, return the value of output at the equilibrium point.
15161516
return_result : bool, optional
1517-
If True, return the `result` option from the scipy root function used
1518-
to compute the equilibrium point.
1517+
If True, return the `result` option from the
1518+
:func:`scipy.optimize.root` function used to compute the equilibrium
1519+
point.
15191520
15201521
Returns
15211522
-------
@@ -1529,9 +1530,9 @@ def find_eqpt(sys, x0, u0=[], y0=None, t=0, params={},
15291530
If `return_y` is True, returns the value of the outputs at the
15301531
equilibrium point, or `None` if no equilibrium point was found and
15311532
`return_result` was False.
1532-
result : scipy root() result object, optional
1533-
If `return_result` is True, returns the `result` from the scipy root
1534-
function.
1533+
result : :class:`scipy.optimize.OptimizeResult`, optional
1534+
If `return_result` is True, returns the `result` from the
1535+
:func:`scipy.optimize.root` function.
15351536
15361537
"""
15371538
from scipy.optimize import root

‎control/matlab/timeresp.py

Lines changed: 8 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -22,7 +22,7 @@ def step(sys, T=None, X0=0., input=0, output=None, return_x=False):
2222
LTI system to simulate
2323
2424
T: array-like or number, optional
25-
Time vector, or simulation time duration if a number (time vector is
25+
Time vector, or simulation time duration if a number (time vector is
2626
autocomputed if not given)
2727
2828
X0: array-like or number, optional
@@ -67,7 +67,7 @@ def step(sys, T=None, X0=0., input=0, output=None, return_x=False):
6767

6868
return yout, T
6969

70-
def stepinfo(sys, T=None, SettlingTimeThreshold=0.02, RiseTimeLimits=(0.1,0.9)):
70+
def stepinfo(sys, T=None, SettlingTimeThreshold=0.02, RiseTimeLimits=(0.1, 0.9)):
7171
'''
7272
Step response characteristics (Rise time, Settling Time, Peak and others).
7373
@@ -77,7 +77,7 @@ def stepinfo(sys, T=None, SettlingTimeThreshold=0.02, RiseTimeLimits=(0.1,0.9)):
7777
LTI system to simulate
7878
7979
T: array-like or number, optional
80-
Time vector, or simulation time duration if a number (time vector is
80+
Time vector, or simulation time duration if a number (time vector is
8181
autocomputed if not given)
8282
8383
SettlingTimeThreshold: float value, optional
@@ -110,7 +110,7 @@ def stepinfo(sys, T=None, SettlingTimeThreshold=0.02, RiseTimeLimits=(0.1,0.9)):
110110
'''
111111
from ..timeresp import step_info
112112

113-
S = step_info(sys, T, SettlingTimeThreshold, RiseTimeLimits)
113+
S = step_info(sys, T, None, SettlingTimeThreshold, RiseTimeLimits)
114114

115115
return S
116116

@@ -130,9 +130,9 @@ def impulse(sys, T=None, X0=0., input=0, output=None, return_x=False):
130130
LTI system to simulate
131131
132132
T: array-like or number, optional
133-
Time vector, or simulation time duration if a number (time vector is
133+
Time vector, or simulation time duration if a number (time vector is
134134
autocomputed if not given)
135-
135+
136136
X0: array-like or number, optional
137137
Initial condition (default = 0)
138138
@@ -186,9 +186,9 @@ def initial(sys, T=None, X0=0., input=None, output=None, return_x=False):
186186
LTI system to simulate
187187
188188
T: array-like or number, optional
189-
Time vector, or simulation time duration if a number (time vector is
189+
Time vector, or simulation time duration if a number (time vector is
190190
autocomputed if not given)
191-
191+
192192
X0: array-like object or number, optional
193193
Initial condition (default = 0)
194194

0 commit comments

Comments
 (0)