Skip to content

Commit b05492c

Browse files
sawyerbfullersawyerbfuller
committed
fixed remaining failing unit tests and cleanup of docfiles
1 parent 1a6d806 commit b05492c

9 files changed

Lines changed: 147 additions & 101 deletions

File tree

‎control/frdata.py‎

Lines changed: 7 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -47,7 +47,7 @@
4747
# External function declarations
4848
from warnings import warn
4949
import numpy as np
50-
from numpy import angle, array, empty, ones, isin, \
50+
from numpy import angle, array, empty, ones, \
5151
real, imag, absolute, eye, linalg, where, dot, sort
5252
from scipy.interpolate import splprep, splev
5353
from .lti import LTI
@@ -354,23 +354,24 @@ def eval(self, omega, squeeze=True):
354354
355355
squeeze: bool, optional (default=True)
356356
If True and sys is single input, single output (SISO), return a
357-
1D array or scalar depending on omega's length.
357+
1D array or scalar the same length as omega.
358358
359359
"""
360360
omega_array = np.array(omega, ndmin=1) # array-like version of omega
361361
if any(omega_array.imag > 0):
362362
raise ValueError("FRD.eval can only accept real-valued omega")
363363

364364
if self.ifunc is None:
365-
elements = isin(self.omega, omega)
365+
elements = np.isin(self.omega, omega) # binary array
366366
if sum(elements) < len(omega_array):
367367
raise ValueError(
368-
"not all frequencies omega are in frequency list of FRD "
369-
"system. Try an interpolating FRD for additional points.")
368+
'''not all frequencies omega are in frequency list of FRD
369+
system. Try an interpolating FRD for additional points.''')
370370
else:
371371
out = self.fresp[:, :, elements]
372372
else:
373-
out = empty((self.outputs, self.inputs, len(omega_array)))
373+
out = empty((self.outputs, self.inputs, len(omega_array)),
374+
dtype=complex)
374375
for i in range(self.outputs):
375376
for j in range(self.inputs):
376377
for k, w in enumerate(omega_array):

‎control/lti.py‎

Lines changed: 7 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -465,9 +465,14 @@ def evalfr(sys, x, squeeze=True):
465465
466466
.. todo:: Add example with MIMO system
467467
"""
468+
out = sys.horner(x)
469+
if not hasattr(x, '__len__'):
470+
# received a scalar x, squeeze down the array along last dim
471+
out = np.squeeze(out, axis=2)
468472
if squeeze and issiso(sys):
469-
return sys.horner(x)[0][0]
470-
return sys.horner(x)
473+
return out[0][0]
474+
else:
475+
return out
471476

472477

473478
def freqresp(sys, omega, squeeze=True):

‎control/margins.py‎

Lines changed: 9 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -213,15 +213,15 @@ def stability_margins(sysdata, returnall=False, epsw=0.0):
213213
# a bit coarse, have the interpolated frd evaluated again
214214
def _mod(w):
215215
"""Calculate |G(jw)| - 1"""
216-
return np.abs(sys(1j * w)[0][0]) - 1
216+
return np.abs(sys(1j * w)) - 1
217217

218218
def _arg(w):
219219
"""Calculate the phase angle at -180 deg"""
220-
return np.angle(-sys(1j * w)[0][0])
220+
return np.angle(-sys(1j * w))
221221

222222
def _dstab(w):
223223
"""Calculate the distance from -1 point"""
224-
return np.abs(sys(1j * w)[0][0] + 1.)
224+
return np.abs(sys(1j * w) + 1.)
225225

226226
# Find all crossings, note that this depends on omega having
227227
# a correct range
@@ -232,7 +232,7 @@ def _dstab(w):
232232

233233
# find the phase crossings ang(H(jw) == -180
234234
widx = np.where(np.diff(np.sign(_arg(sys.omega))))[0]
235-
widx = widx[np.real(sys(1j * sys.omega[widx])[0][0]) <= 0]
235+
widx = widx[np.real(sys(1j * sys.omega[widx])) <= 0]
236236
w_180 = np.array(
237237
[sp.optimize.brentq(_arg, sys.omega[i], sys.omega[i+1])
238238
for i in widx])
@@ -296,11 +296,10 @@ def phase_crossover_frequencies(sys):
296296
(array([ 1.73205081, 0. ]), array([-0.5 , 0.25]))
297297
"""
298298

299+
if not sys.issiso():
300+
raise ValueError("MIMO systems not yet implemented.")
299301
# Convert to a transfer function
300302
tf = xferfcn._convert_to_transfer_function(sys)
301-
302-
# if not siso, fall back to (0,0) element
303-
#! TODO: should add a check and warning here
304303
num = tf.num[0][0]
305304
den = tf.den[0][0]
306305

@@ -313,7 +312,9 @@ def phase_crossover_frequencies(sys):
313312

314313
# using real() to avoid rounding errors and results like 1+0j
315314
# it would be nice to have a vectorized version of self.evalfr here
316-
gain = np.real(np.asarray([tf(1j * f)[0][0] for f in realposfreq]))
315+
# update Sawyer B. Fuller 2020.08.15: your wish is my command.
316+
#gain = np.real(np.asarray([tf(1j * f) for f in realposfreq]))
317+
gain = np.real(tf(1j * realposfreq))
317318

318319
return realposfreq, gain
319320

‎control/rlocus.py‎

Lines changed: 5 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -493,7 +493,7 @@ def _RLFindRoots(nump, denp, kvect):
493493
"""Find the roots for the root locus."""
494494
# Convert numerator and denominator to polynomials if they aren't
495495
roots = []
496-
for k in kvect:
496+
for k in np.array(kvect, ndmin=1):
497497
curpoly = denp + k * nump
498498
curroots = curpoly.r
499499
if len(curroots) < denp.order:
@@ -577,10 +577,10 @@ def _RLFeedbackClicksPoint(event, sys, fig, ax_rlocus, sisotool=False):
577577
# Catch type error when event click is in the figure but not in an axis
578578
try:
579579
s = complex(event.xdata, event.ydata)
580-
K = -1. / sys.horner(s)
581-
K_xlim = -1. / sys.horner(
580+
K = -1. / sys(s)
581+
K_xlim = -1. / sys(
582582
complex(event.xdata + 0.05 * abs(xlim[1] - xlim[0]), event.ydata))
583-
K_ylim = -1. / sys.horner(
583+
K_ylim = -1. / sys(
584584
complex(event.xdata, event.ydata + 0.05 * abs(ylim[1] - ylim[0])))
585585

586586
except TypeError:
@@ -621,7 +621,7 @@ def _RLFeedbackClicksPoint(event, sys, fig, ax_rlocus, sisotool=False):
621621
ax_rlocus.plot(s.real, s.imag, 'k.', marker='s', markersize=8,
622622
zorder=20, label='gain_point')
623623

624-
return K.real[0][0]
624+
return K.real
625625

626626

627627
def _removeLine(label, ax):

‎control/statesp.py‎

Lines changed: 75 additions & 47 deletions
Original file line numberDiff line numberDiff line change
@@ -437,60 +437,91 @@ def __rdiv__(self, other):
437437

438438
raise NotImplementedError("StateSpace.__rdiv__ is not implemented yet.")
439439

440-
def __call__(self, s, squeeze=True):
441-
"""Evaluate system's transfer function at complex frequencies s/z.
440+
def __call__(self, x, squeeze=True):
441+
"""Evaluate system's transfer function at complex frequencies.
442+
443+
Evaluates at complex frequency x, where x is s or z dependong on
444+
whether the system is continuous or discrete-time.
442445
443446
If squeeze is True (default) and sys is single input, single output
444-
(SISO), return a 1D array or scalar depending on s.
447+
(SISO), return a 1D array or scalar depending on the size of x.
448+
For a MIMO system, returns an (n_outputs, n_inputs, n_x) array.
445449
446450
"""
447-
# Use TB05AD from Slycot if available
451+
# Use Slycot if available
448452
try:
449-
from slycot import tb05ad
450-
451-
# preallocate
452-
s_arr = np.array(s, ndmin=1) # array-like version of s
453-
out = np.empty((self.outputs, self.inputs, len(s_arr)),
454-
dtype=complex)
455-
n = self.states
456-
m = self.inputs
457-
p = self.outputs
458-
# The first call both evaluates C(sI-A)^-1 B and also returns
459-
# Hessenberg transformed matrices at, bt, ct.
460-
result = tb05ad(n, m, p, s_arr[0], self.A,
461-
self.B, self.C, job='NG')
462-
# When job='NG', result = (at, bt, ct, g_i, hinvb, info)
463-
at = result[0]
464-
bt = result[1]
465-
ct = result[2]
466-
467-
# TB05AD frequency evaluation does not include direct feedthrough.
468-
out[:, :, 0] = result[3] + self.D
469-
470-
# Now, iterate through the remaining frequencies using the
471-
# transformed state matrices, at, bt, ct.
472-
473-
# Start at the second frequency, already have the first.
474-
for kk, s_kk in enumerate(s_arr[1:len(s_arr)]):
475-
result = tb05ad(n, m, p, s_kk, at, bt, ct, job='NH')
476-
# When job='NH', result = (g_i, hinvb, info)
477-
478-
# kk+1 because enumerate starts at kk = 0.
479-
# but zero-th spot is already filled.
480-
out[:, :, kk+1] = result[0] + self.D
481-
482-
if not hasattr(s, '__len__'):
483-
# received a scalar s, squeeze down the array
484-
out = np.squeeze(out, axis=2)
485-
except ImportError: # Slycot unavailable. Fall back to horner.
486-
out = self.horner(s)
453+
out = self.slycot_horner(x)
454+
except ImportError: # Slycot unavailable. use built-in horner.
455+
out = self.horner(x)
456+
if not hasattr(x, '__len__'):
457+
# received a scalar x, squeeze down the array along last dim
458+
out = np.squeeze(out, axis=2)
487459
if squeeze and self.issiso():
488460
return out[0][0]
489461
else:
490462
return out
491463

464+
def slycot_horner(self, s):
465+
"""Evaluate system's transfer function at complex frequencies s
466+
using Horner's method from Slycot.
467+
468+
Expects inputs and outputs to be formatted correctly. Use __call__
469+
for a more user-friendly interface.
470+
471+
Parameters
472+
s : array-like
473+
474+
Returns
475+
output : array of size (outputs, inputs, len(s))
476+
477+
"""
478+
from slycot import tb05ad
479+
480+
# preallocate
481+
s_arr = np.array(s, ndmin=1) # array-like version of s
482+
out = np.empty((self.outputs, self.inputs, len(s_arr)),
483+
dtype=complex)
484+
n = self.states
485+
m = self.inputs
486+
p = self.outputs
487+
# The first call both evaluates C(sI-A)^-1 B and also returns
488+
# Hessenberg transformed matrices at, bt, ct.
489+
result = tb05ad(n, m, p, s_arr[0], self.A,
490+
self.B, self.C, job='NG')
491+
# When job='NG', result = (at, bt, ct, g_i, hinvb, info)
492+
at = result[0]
493+
bt = result[1]
494+
ct = result[2]
495+
496+
# TB05AD frequency evaluation does not include direct feedthrough.
497+
out[:, :, 0] = result[3] + self.D
498+
499+
# Now, iterate through the remaining frequencies using the
500+
# transformed state matrices, at, bt, ct.
501+
502+
# Start at the second frequency, already have the first.
503+
for kk, s_kk in enumerate(s_arr[1:len(s_arr)]):
504+
result = tb05ad(n, m, p, s_kk, at, bt, ct, job='NH')
505+
# When job='NH', result = (g_i, hinvb, info)
506+
507+
# kk+1 because enumerate starts at kk = 0.
508+
# but zero-th spot is already filled.
509+
out[:, :, kk+1] = result[0] + self.D
510+
return out
511+
492512
def horner(self, s):
493-
"""Evaluate systems's transfer function at complex frequencies s.
513+
"""Evaluate system's transfer function at complex frequencies s
514+
using Horner's method.
515+
516+
Expects inputs and outputs to be formatted correctly. Use __call__
517+
for a more user-friendly interface.
518+
519+
Parameters
520+
s : array-like
521+
522+
Returns
523+
output : array of size (outputs, inputs, len(s))
524+
494525
"""
495526
s_arr = np.array(s, ndmin=1) # force to be an array
496527
# Preallocate
@@ -501,9 +532,6 @@ def horner(self, s):
501532
np.dot(self.C,
502533
solve(s_idx * eye(self.states) - self.A, self.B)) \
503534
+ self.D
504-
if not hasattr(s, '__len__'):
505-
# received a scalar s, squeeze down the array along last dim
506-
out = np.squeeze(out, axis=2)
507535
return out
508536

509537
def freqresp(self, omega):
@@ -879,7 +907,7 @@ def dcgain(self):
879907
if self.isctime():
880908
gain = np.asarray(self.D-self.C.dot(np.linalg.solve(self.A, self.B)))
881909
else:
882-
gain = self.horner(1)
910+
gain = np.squeeze(self.horner(1))
883911
except LinAlgError:
884912
# eigenvalue at DC
885913
gain = np.tile(np.nan, (self.outputs, self.inputs))

‎control/tests/frd_test.py‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -387,9 +387,9 @@ def test_eval(self):
387387
np.testing.assert_almost_equal(sys_tf(1j), frd_tf(1j))
388388

389389
# Should get an error if we evaluate at an unknown frequency
390-
self.assertRaises(Exception, frd_tf.eval(2))
390+
self.assertRaises(ValueError, frd_tf.eval, 2)
391391
# Should get an error if we use __call__ at real-valued frequency
392-
self.assertRaises(ValueError, frd_tf(2))
392+
self.assertRaises(ValueError, frd_tf, 2)
393393

394394
def test_repr_str(self):
395395
# repr printing

‎control/tests/margin_test.py‎

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -76,7 +76,7 @@ def test_stability_margins_omega(sys, refout, refoutall):
7676
def test_stability_margins_3input(sys, refout, refoutall):
7777
"""Test stability_margins() function with mag, phase, omega input"""
7878
omega = np.logspace(-2, 2, 2000)
79-
mag, phase, omega_ = sys.freqresp(omega)
79+
mag, phase, omega_ = sys.frequency_response(omega)
8080
out = stability_margins((mag, phase*180/np.pi, omega_))
8181
assert_allclose(out, refout, atol=1.5e-3)
8282

@@ -92,7 +92,7 @@ def test_margin_sys(sys, refout, refoutall):
9292
def test_margin_3input(sys, refout, refoutall):
9393
"""Test margin() function with mag, phase, omega input"""
9494
omega = np.logspace(-2, 2, 2000)
95-
mag, phase, omega_ = sys.freqresp(omega)
95+
mag, phase, omega_ = sys.frequency_response(omega)
9696
out = margin((mag, phase*180/np.pi, omega_))
9797
assert_allclose(out, np.array(refout)[[0, 1, 3, 4]], atol=1.5e-3)
9898

@@ -113,7 +113,7 @@ def test_phase_crossover_frequencies():
113113
[[3], [4]]],
114114
[[[1, 2, 3, 4], [1, 1]],
115115
[[1, 1], [1, 1]]])
116-
omega, gain = phase_crossover_frequencies(tf)
116+
omega, gain = phase_crossover_frequencies(tf[0,0])
117117
assert_allclose(omega, [1.73205, 0.], atol=1.5e-3)
118118
assert_allclose(gain, [-0.5, 0.25], atol=1.5e-3)
119119

@@ -123,7 +123,7 @@ def test_mag_phase_omega():
123123
sys = TransferFunction(15, [1, 6, 11, 6])
124124
out = stability_margins(sys)
125125
omega = np.logspace(-2, 2, 1000)
126-
mag, phase, omega = sys.freqresp(omega)
126+
mag, phase, omega = sys.frequency_response(omega)
127127
out2 = stability_margins((mag, phase*180/np.pi, omega))
128128
ind = [0, 1, 3, 4] # indices of gm, pm, wg, wp -- ignore sm
129129
marg1 = np.array(out)[ind]

‎control/tests/statesp_array_test.py‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -505,7 +505,7 @@ def test_horner(self):
505505
# Make sure result agrees with frequency response
506506
mag, phase, omega = self.sys322.frequency_response([1])
507507
np.testing.assert_array_almost_equal(
508-
self.sys322.horner(1.j),
508+
np.squeeze(self.sys322.horner(1.j)),
509509
mag[:,:,0] * np.exp(1.j * phase[:,:,0]))
510510

511511
def tearDown(self):

0 commit comments

Comments
 (0)