Skip to content

Commit 00e4052

Browse files
repaghrepagh
authored
Merge branch 'master' into modal_form
2 parents 136a4ec + 077a8bc commit 00e4052

10 files changed

Lines changed: 269 additions & 189 deletions

File tree

‎.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

‎control/canonical.py‎

Lines changed: 4 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -135,7 +135,10 @@ def observable_form(xsys):
135135
Wrz = obsv(zsys.A, zsys.C)
136136

137137
# Transformation from one form to another
138-
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.")
139142

140143
# Finally, compute the output matrix
141144
zsys.B = Tzx * xsys.B

‎control/statesp.py‎

Lines changed: 43 additions & 24 deletions
Original file line numberDiff line numberDiff line change
@@ -54,8 +54,7 @@
5454
import math
5555
import numpy as np
5656
from numpy import all, angle, any, array, asarray, concatenate, cos, delete, \
57-
dot, empty, exp, eye, matrix, ones, poly, poly1d, roots, shape, sin, \
58-
zeros, squeeze
57+
dot, empty, exp, eye, isinf, matrix, ones, pad, shape, sin, zeros, squeeze
5958
from numpy.random import rand, randn
6059
from numpy.linalg import solve, eigvals, matrix_rank
6160
from numpy.linalg.linalg import LinAlgError
@@ -516,17 +515,41 @@ def pole(self):
516515
def zero(self):
517516
"""Compute the zeros of a state space system."""
518517

519-
if self.inputs > 1 or self.outputs > 1:
520-
raise NotImplementedError("StateSpace.zeros is currently \
521-
implemented only for SISO systems.")
518+
if not self.states:
519+
return np.array([])
522520

523-
den = poly1d(poly(self.A))
524-
# Compute the numerator based on zeros
525-
#! TODO: This is currently limited to SISO systems
526-
num = poly1d(poly(self.A - dot(self.B, self.C)) + ((self.D[0, 0] - 1) *
527-
den))
528-
529-
return roots(num)
521+
# Use AB08ND from Slycot if it's available, otherwise use
522+
# scipy.lingalg.eigvals().
523+
try:
524+
from slycot import ab08nd
525+
526+
out = ab08nd(self.A.shape[0], self.B.shape[1], self.C.shape[0],
527+
self.A, self.B, self.C, self.D)
528+
nu = out[0]
529+
return sp.linalg.eigvals(out[8][0:nu,0:nu], out[9][0:nu,0:nu])
530+
except ImportError: # Slycot unavailable. Fall back to scipy.
531+
if self.C.shape[0] != self.D.shape[1]:
532+
raise NotImplementedError("StateSpace.zero only supports "
533+
"systems with the same number of "
534+
"inputs as outputs.")
535+
536+
# This implements the QZ algorithm for finding transmission zeros
537+
# from
538+
# https://dspace.mit.edu/bitstream/handle/1721.1/841/P-0802-06587335.pdf.
539+
# The QZ algorithm solves the generalized eigenvalue problem: given
540+
# `L = [A, B; C, D]` and `M = [I_nxn 0]`, find all finite lambda
541+
# for which there exist nontrivial solutions of the equation
542+
# `Lz - lamba Mz`.
543+
#
544+
# The generalized eigenvalue problem is only solvable if its
545+
# arguments are square matrices.
546+
L = concatenate((concatenate((self.A, self.B), axis=1),
547+
concatenate((self.C, self.D), axis=1)), axis=0)
548+
M = pad(eye(self.A.shape[0]), ((0, self.C.shape[0]),
549+
(0, self.B.shape[1])), "constant")
550+
return np.array([x for x in sp.linalg.eigvals(L, M,
551+
overwrite_a=True)
552+
if not isinf(x)])
530553

531554
# Feedback around a state space system
532555
def feedback(self, other=1, sign=-1):
@@ -757,7 +780,6 @@ def _convertToStateSpace(sys, **kw):
757780
[1., 1., 1.]].
758781
759782
"""
760-
761783
from .xferfcn import TransferFunction
762784
import itertools
763785
if isinstance(sys, StateSpace):
@@ -771,20 +793,17 @@ def _convertToStateSpace(sys, **kw):
771793
try:
772794
from slycot import td04ad
773795
if len(kw):
774-
raise TypeError("If sys is a TransferFunction, _convertToStateSpace \
775-
cannot take keywords.")
796+
raise TypeError("If sys is a TransferFunction, "
797+
"_convertToStateSpace cannot take keywords.")
776798

777799
# Change the numerator and denominator arrays so that the transfer
778800
# function matrix has a common denominator.
779-
num, den = sys._common_den()
780-
# Make a list of the orders of the denominator polynomials.
781-
index = [len(den) - 1 for i in range(sys.outputs)]
782-
# Repeat the common denominator along the rows.
783-
den = array([den for i in range(sys.outputs)])
784-
#! TODO: transfer function to state space conversion is still buggy!
785-
#print num
786-
#print shape(num)
787-
ssout = td04ad('R',sys.inputs, sys.outputs, index, den, num,tol=0.0)
801+
# matrices are also sized/padded to fit td04ad
802+
num, den, denorder = sys.minreal()._common_den()
803+
804+
# transfer function to state space conversion now should work!
805+
ssout = td04ad('C', sys.inputs, sys.outputs,
806+
denorder, den, num, tol=0)
788807

789808
states = ssout[0]
790809
return StateSpace(ssout[1][:states, :states],

‎control/tests/canonical_test.py‎

Lines changed: 44 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -82,3 +82,47 @@ def test_modal_form(self):
8282
#np.testing.assert_array_almost_equal(sys_check.C, C_true)
8383
np.testing.assert_array_almost_equal(sys_check.D, D_true)
8484
#np.testing.assert_array_almost_equal(T_check, T_true)
85+
86+
def test_observable_form(self):
87+
"""Test the observable canonical form"""
88+
89+
# Create a system in the observable canonical form
90+
coeffs = [1.0, 2.0, 3.0, 4.0, 1.0]
91+
A_true = np.polynomial.polynomial.polycompanion(coeffs)
92+
A_true = np.fliplr(np.flipud(A_true))
93+
B_true = np.matrix("1.0 1.0 1.0 1.0").T
94+
C_true = np.matrix("1.0 0.0 0.0 0.0")
95+
D_true = 42.0
96+
97+
# Perform a coordinate transform with a random invertible matrix
98+
T_true = np.matrix([[-0.27144004, -0.39933167, 0.75634684, 0.44135471],
99+
[-0.74855725, -0.39136285, -0.18142339, -0.50356997],
100+
[-0.40688007, 0.81416369, 0.38002113, -0.16483334],
101+
[-0.44769516, 0.15654653, -0.50060858, 0.72419146]])
102+
A = np.linalg.solve(T_true, A_true)*T_true
103+
B = np.linalg.solve(T_true, B_true)
104+
C = C_true*T_true
105+
D = D_true
106+
107+
# Create a state space system and convert it to the observable canonical form
108+
sys_check, T_check = canonical_form(ss(A, B, C, D), "observable")
109+
110+
# Check against the true values
111+
np.testing.assert_array_almost_equal(sys_check.A, A_true)
112+
np.testing.assert_array_almost_equal(sys_check.B, B_true)
113+
np.testing.assert_array_almost_equal(sys_check.C, C_true)
114+
np.testing.assert_array_almost_equal(sys_check.D, D_true)
115+
np.testing.assert_array_almost_equal(T_check, T_true)
116+
117+
def test_unobservable_system(self):
118+
"""Test observable canonical form with an unobservable system"""
119+
120+
# Create an unobservable system
121+
A = np.matrix("1.0 2.0 2.0; 4.0 5.0 5.0; 7.0 8.0 8.0")
122+
B = np.matrix("1.0 1.0 1.0").T
123+
C = np.matrix("1.0 1.0 1.0")
124+
D = 42.0
125+
sys = ss(A, B, C, D)
126+
127+
# Check if an exception is raised
128+
np.testing.assert_raises(ValueError, canonical_form, sys, "observable")

‎control/tests/convert_test.py‎

Lines changed: 21 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -40,7 +40,7 @@ def setUp(self):
4040
# Set to True to print systems to the output.
4141
self.debug = False
4242
# get consistent results
43-
np.random.seed(9)
43+
np.random.seed(7)
4444

4545
def printSys(self, sys, ind):
4646
"""Print system to the standard output."""
@@ -141,10 +141,9 @@ def testConvert(self):
141141
ssxfrm_real = ssxfrm_mag * np.cos(ssxfrm_phase)
142142
ssxfrm_imag = ssxfrm_mag * np.sin(ssxfrm_phase)
143143
np.testing.assert_array_almost_equal( \
144-
ssorig_real, ssxfrm_real)
144+
ssorig_real, ssxfrm_real)
145145
np.testing.assert_array_almost_equal( \
146-
ssorig_imag, ssxfrm_imag)
147-
146+
ssorig_imag, ssxfrm_imag)
148147
#
149148
# Make sure xform'd TF has same frequency response
150149
#
@@ -198,8 +197,9 @@ def testTf2ssStaticMimo(self):
198197
"""Regression: tf2ss for MIMO static gain"""
199198
import control
200199
# 2x3 TFM
201-
gmimo = control.tf2ss(control.tf([[ [23], [3], [5] ], [ [-1], [0.125], [101.3] ]],
202-
[[ [46], [0.1], [80] ], [ [2], [-0.1], [1] ]]))
200+
gmimo = control.tf2ss(control.tf(
201+
[[ [23], [3], [5] ], [ [-1], [0.125], [101.3] ]],
202+
[[ [46], [0.1], [80] ], [ [2], [-0.1], [1] ]]))
203203
self.assertEqual(0, gmimo.states)
204204
self.assertEqual(3, gmimo.inputs)
205205
self.assertEqual(2, gmimo.outputs)
@@ -229,6 +229,21 @@ def testSs2tfStaticMimo(self):
229229
numref = np.asarray(d)[...,np.newaxis]
230230
np.testing.assert_array_equal(numref, np.array(gtf.num) / np.array(gtf.den))
231231

232+
def testTf2SsDuplicatePoles(self):
233+
"""Tests for "too few poles for MIMO tf #111" """
234+
import control
235+
try:
236+
import slycot
237+
num = [ [ [1], [0] ],
238+
[ [0], [1] ] ]
239+
240+
den = [ [ [1,0], [1] ],
241+
[ [1], [1,0] ] ]
242+
g = control.tf(num, den)
243+
s = control.ss(g)
244+
np.testing.assert_array_equal(g.pole(), s.pole())
245+
except ImportError:
246+
print("Slycot not present, skipping")
232247

233248
def suite():
234249
return unittest.TestLoader().loadTestsFromTestCase(TestConvert)

‎control/tests/discrete_test.py‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -6,6 +6,7 @@
66
import unittest
77
import numpy as np
88
from control import *
9+
from control import matlab
910

1011
class TestDiscrete(unittest.TestCase):
1112
"""Tests for the DiscreteStateSpace class."""

‎control/tests/matlab_test.py‎

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -396,9 +396,9 @@ def testBalred(self):
396396
@unittest.skipIf(not slycot_check(), "slycot not installed")
397397
def testModred(self):
398398
modred(self.siso_ss1, [1])
399-
modred(self.siso_ss2 * self.siso_ss3, [0, 1])
400-
modred(self.siso_ss3, [1], 'matchdc')
401-
modred(self.siso_ss3, [1], 'truncate')
399+
modred(self.siso_ss2 * self.siso_ss1, [0, 1])
400+
modred(self.siso_ss1, [1], 'matchdc')
401+
modred(self.siso_ss1, [1], 'truncate')
402402

403403
@unittest.skipIf(not slycot_check(), "slycot not installed")
404404
def testPlace_varga(self):

‎control/tests/statesp_test.py‎

Lines changed: 29 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -42,15 +42,37 @@ def testPole(self):
4242

4343
np.testing.assert_array_almost_equal(p, true_p)
4444

45-
def testZero(self):
46-
"""Evaluate the zeros of a SISO system."""
45+
@unittest.skipIf(not slycot_check(), "slycot not installed")
46+
def testMIMOZero_nonsquare(self):
47+
"""Evaluate the zeros of a MIMO system."""
48+
49+
z = np.sort(self.sys1.zero())
50+
true_z = np.sort([44.41465, -0.490252, -5.924398])
51+
52+
np.testing.assert_array_almost_equal(z, true_z)
53+
54+
A = np.array([[1, 0, 0, 0, 0, 0],
55+
[0, 1, 0, 0, 0, 0],
56+
[0, 0, 3, 0, 0, 0],
57+
[0, 0, 0,-4, 0, 0],
58+
[0, 0, 0, 0,-1, 0],
59+
[0, 0, 0, 0, 0, 3]])
60+
B = np.array([[0,-1],
61+
[-1,0],
62+
[1,-1],
63+
[0, 0],
64+
[0, 1],
65+
[-1,-1]])
66+
C = np.array([[1, 0, 0, 1, 0, 0],
67+
[0, 1, 0, 1, 0, 1],
68+
[0, 0, 1, 0, 0, 1]])
69+
D = np.zeros((3,2))
70+
sys = StateSpace(A, B, C, D)
4771

48-
sys = StateSpace(self.sys1.A, [[3.], [-2.], [4.]], [[-1., 3., 2.]], [[-4.]])
49-
z = sys.zero()
72+
z = np.sort(sys.zero())
73+
true_z = np.sort([2., -1.])
5074

51-
np.testing.assert_array_almost_equal(z, [4.26864638637134,
52-
-3.75932319318567 + 1.10087776649554j,
53-
-3.75932319318567 - 1.10087776649554j])
75+
np.testing.assert_array_almost_equal(z, true_z)
5476

5577
def testAdd(self):
5678
"""Add two MIMO systems."""

‎control/tests/xferfcn_test.py‎

Lines changed: 9 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -425,10 +425,16 @@ def testPoleMIMO(self):
425425
[[[1., 2.], [1., 3.]], [[1., 4., 4.], [1., 9., 14.]]])
426426
p = sys.pole()
427427

428-
np.testing.assert_array_almost_equal(p, [-7., -3., -2., -2.])
429-
430-
# Tests for TransferFunction.feedback.
428+
np.testing.assert_array_almost_equal(p, [-2., -2., -7., -3., -2.])
431429

430+
@unittest.skipIf(not slycot_check(), "slycot not installed")
431+
def testDoubleCancelingPoleSiso(self):
432+
433+
H = TransferFunction([1,1],[1,2,1])
434+
p = H.pole()
435+
np.testing.assert_array_almost_equal(p, [-1, -1])
436+
437+
# Tests for TransferFunction.feedback
432438
def testFeedbackSISO(self):
433439
"""Test for correct SISO transfer function feedback."""
434440

0 commit comments

Comments
 (0)