8000 StateSpace.zero() now supports MIMO systems · MaxGaukler/python-control@6e65070 · GitHub
[go: up one dir, main page]
Skip to content

Commit 6e65070

Browse files
calcmogulcalcmogul
committed
StateSpace.zero() now supports MIMO systems
scipy.linalg.eigvals() is used to solve the generalized eigenvalue problem. The zero locations used in the unit test were obtained via the zero() function in MATLAB R2017b with two more digits of precision added based on the results of StateSpace.zero().
1 parent 601b581 commit 6e65070

File tree

2 files changed

+61
-18
lines changed

2 files changed

+61
-18
lines changed

‎control/statesp.py

Lines changed: 34 additions & 12 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,40 @@ 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 λ for
541+
# which there exist nontrivial solutions of the equation `Lz - λMz`.
542+
#
543+
# The generalized eigenvalue problem is only solvable if its
544+
# arguments are square matrices.
545+
L = concatenate((concatenate((self.A, self.B), axis=1),
546+
concatenate((self.C, self.D), axis=1)), axis=0)
547+
M = pad(eye(self.A.shape[0]), ((0, self.C.shape[0]),
548+
(0, self.B.shape[1])), "constant")
549+
return np.array([x for x in sp.linalg.eigvals(L, M,
550+
overwrite_a=True)
551+
if not isinf(x)])
530552

531553
# Feedback around a state space system
532554
def feedback(self, other=1, sign=-1):

‎control/tests/statesp_test.py

Lines changed: 27 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -43,14 +43,35 @@ def testPole(self):
4343
np.testing.assert_array_almost_equal(p, true_p)
4444

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

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

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

5576
def testAdd(self):
5677
"""Add two MIMO systems."""

0 commit comments

Comments
 (0)
0