Skip to content

Complex matrices #376

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Closed
lytex wants to merge 12 commits into from
Closed

Complex matrices #376

Show file tree
Hide file tree
Changes from 5 commits
Commits
Show all changes
12 commits
Select commit Hold shift + click to select a range
e6f2f60
allow _ssmatrix to accept complex values
lytex Feb 9, 2020
53dbfe4
add new tests with complex state matrices
lytex Feb 29, 2020
05d730a
rename to test_evalfr_complex, use self.sysC332
lytex Feb 29, 2020
a55165a
missed one sys.evalfr in previous commit
lytex Feb 29, 2020
2d231a8
missed one sys.evalfr in previous commit
lytex Feb 29, 2020
f1c9b9e
add complex to test_freq_resp
lytex Mar 3, 2020
e7a93ea
remove warnings settings for debugging
lytex Mar 3, 2020
bb44827
ease precision of reference values in test
lytex Mar 3, 2020
fc5c76b
Merge branch 'master' into complex
lytex Mar 30, 2020
5f58723
build an auxiliary system to compare against
lytex Mar 30, 2020
1803d38
build an auxiliary system to compare against
lytex Mar 30, 2020
84bc8ef
add code to use slycot complex function ab08nz if any matrix is complex
lytex Apr 1, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
10 changes: 8 additions & 2 deletions control/statesp.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -93,9 +93,15 @@ def _ssmatrix(data, axis=1):
# Convert the data into an array or matrix, as configured
# If data is passed as a string, use (deprecated?) matrix constructor
if config.defaults['statesp.use_numpy_matrix'] or isinstance(data, str):
arr = np.matrix(data, dtype=float)
if np.isrealobj(data):
arr = np.matrix(data, dtype=float)
elif np.iscomplexobj(data):
arr = np.matrix(data, dtype=complex)
else:
arr = np.array(data, dtype=float)
if np.isrealobj(data):
arr = np.array(data, dtype=float)
elif np.iscomplexobj(data):
arr = np.array(data, dtype=complex)
Comment on lines -96 to +104

@bnavigator bnavigator Jul 15, 2020

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is this really necessary? Why not just use matrix() and asarray() without dtype and figure out the dtype later?

@bnavigator bnavigator Jul 29, 2020 •
edited
Loading

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Turns out if dtype is omitted completely, it could still be int and trip operations later on.

The test suite already has complex casting warnings from Scipy and Slycot functions returning complex arrays (with 0 imaginary part), which are passed to _ssmatrix. Thus, my approach in #438 is to use dtype=np.result_type(np.float, arr), which ensures to cast from int to at least float and retains complex if arr is already complex.

ndim = arr.ndim
shape = arr.shape

Expand Down
215 changes: 215 additions & 0 deletions control/tests/statesp_test.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -63,6 +63,50 @@ def setUp(self):
D623 = np.zeros((3, 2))
self.sys623 = StateSpace(A623, B623, C623, D623)

# Systems with complex matrices, sysC322, sysC222, sysC623
# These systems have the same shape as the previous ones

A = np.array([[6.0 + 9.0j, 8.0 + 9.0j, -4.0 - 7.0j],
[8.0 - 7.0j, 3.0, 1.0 - 1.0j],
[-7.0 + 9.0j, -8.0 + 6.0j, 9.0 + 8.0j]])
B = np.array([[6.0 + 3.0j, -9.0 - 2.0j],
[9.0 + 5.0j, 7.0 + 3.0j],
[3.0 + 5.0j, 8.0 - 6.0j]])
C = np.array([[4.0 + 4.0j, -4.0 + 9.0j, -8.0 - 1.0j],
[-9.0 - 3.0j, -9.0 - 9.0j, 6.0 - 2.0j]])
D = np.array([[5.0 - 1.0j, -6.0 + 4.0j],
[6.0 + 3.0j, 5.0j]])
self.sysC322 = StateSpace(A, B, C, D)

A = np.array([[-4.0 - 7.0j, 3.0 + 9.0j],
[3.0, -6.0 - 3.0j]])
B = np.array([[2.0, 5.0 + 7.0j],
[-5.0 + 4.0j, -5.0 + 9.0j]])
C = np.array([[1.0 + 6.0j, -7.0 + 6.0j],
[-7.0 - 5.0j, -5.0 - 5.0j]])
D = np.array([[8.0 + 2.0j, -6.0 - 3.0j],
[-3.0 - 1.0j, -5.0 + 6.0j]])
self.sysC222 = StateSpace(A, B, C, D)

A = np.array(
[[2.0 - 8.0j, -2.0 + 6.0j, 8.0 - 1.0j, -6.0 + 7.0j, -5.0 - 3.0j, -5.0 - 6.0j],
[1.0 - 1.0j, 1.0 + 7.0j, -7.0 + 8.0j, 6.0 + 2.0j, 3.0, 8.0 - 5.0j],
[8.0 - 7.0j, -8.0 - 8.0j, 1.0 - 6.0j, -4.0 + 1.0j, 4.0 - 2.0j, -7.0 - 2.0j],
[-4.0 + 9.0j, -8.0 - 2.0j, -1.0 - 4.0j, 1.0 - 7.0j, 5.0 - 8.0j, 6.0 - 9.0j],
[5.0 - 9.0j, 1.0 - 5.0j, -9.0 - 7.0j, -6.0 + 7.0j, -1.0 - 5.0j, 1.0 + 8.0j],
[5.0 + 5.0j, 5.0 + 6.0j, -3.0 - 7.0j, 2.0 + 2.0j, -8.0 - 7.0j, 9.0 + 8.0j]])
B = np.array([[6.0j, 5.0 + 3.0j, 8.0 + 4.0j],
[-9.0j, -2.0 - 1.0j, 9.0 - 6.0j],
[-3.0 - 9.0j, -5.0 + 1.0j, 1.0 - 2.0j],
[8.0 - 6.0j, -2.0 - 4.0j, -8.0 + 2.0j],
[-2.0 + 3.0j, -8.0 + 5.0j, -5.0 + 5.0j],
[-7.0 + 4.0j, -7.0 - 6.0j, -3.0 - 8.0j]])
C = np.array([[8.0 + 6.0j, -3.0j, -1.0 + 7.0j, 2.0j, 6.0 - 6.0j, 3.0 - 1.0j],
[5.0 + 1.0j, -1.0 + 8.0j, -4.0 + 1.0j, 2.0j, 6.0 - 4.0j, -2.0 - 5.0j]])
D = np.array([[7.0 - 4.0j, -5.0 - 1.0j, -5.0 + 8.0j],
[-6.0 + 8.0j, -6.0 - 6.0j, -1.0 + 9.0j]])
self.sysC623 = StateSpace(A, B, C, D)

def test_D_broadcast(self):
"""Test broadcast of D=0 to the right shape"""
# Giving D as a scalar 0 should broadcast to the right shape
Expand Down Expand Up @@ -519,6 +563,177 @@ def test_lft(self):
np.testing.assert_allclose(np.array(pk.C).reshape(-1), Cmatlab)
np.testing.assert_allclose(np.array(pk.D).reshape(-1), Dmatlab)

def test_pole_complex(self):
"""Evaluate the poles of a complex MIMO system."""

p = np.sort(self.sysC322.pole())
true_p = np.sort([21.7521962567499187 +7.8381752267389437j
-6.4620734598457208 +2.8701578198630169j,
2.7098772030958163 +6.2916669533980274j])

np.testing.assert_array_almost_equal(p, true_p)


@unittest.skipIf(not slycot_check(), "slycot not installed")
def test_zero_siso_complex(self):
"""Evaluate the zeros of a complex SISO system."""
# extract only first input / first output system of sysC222. This system is denoted sysC111
# or tfC111
tfC111 = ss2tf(self.sysC222)
sysC111 = tf2ss(tfC111[0, 0])

# compute zeros as root of the characteristic polynomial at the numerator of tfC111
# this method is simple and assumed as valid in this test
true_z = np.sort(tfC111[0, 0].zero())
# Compute the zeros through ab08nd, which is tested here
z = np.sort(sysC111.zero())

np.testing.assert_almost_equal(true_z, z)

@unittest.skipIf(not slycot_check(), "slycot not installed")
def test_zero_mimo_sysC322_square_complex(self):
"""Evaluate the zeros of a square complex MIMO system."""

z = np.sort(self.sysC322.zero())
true_z = np.sort([36.4937595620286217 +8.0738640861708575j,
-4.7860079612333388 +29.3266582804945379j,
-7.4509692664104161 +5.5262915134608006j])
np.testing.assert_array_almost_equal(z, true_z)

@unittest.skipIf(not slycot_check(), "slycot not installed")
def test_zero_mimo_sysC222_square_complex(self):
"""Evaluate the zeros of a square complex MIMO system."""

z = np.sort(self.sysC222.zero())
true_z = np.sort([-10.5685005560737366,
3.3685005560737391])
np.testing.assert_array_almost_equal(z, true_z)

@unittest.skipIf(not slycot_check(), "slycot not installed")
def test_zero_mimo_sysC623_non_square_complex(self):
"""Evaluate the zeros of a non square complex MIMO system."""

z = np.sort(self.sysC623.zero())
true_z = np.sort([]) # System has no transmission zeros, not sure what slycot returns
np.testing.assert_array_almost_equal(z, true_z)

def test_add_ss_complex(self):
"""Add two complex MIMO systems."""

A = np.array([[6.0 +9.0j, 8.0 +9.0j, -4.0 -7.0j, 0.0, 0.0],
[8.0 -7.0j, 3.0, 1.0 -1.0j, 0.0, 0.0],
[-7.0 +9.0j, -8.0 +6.0j, 9.0 +8.0j, 0.0, 0.0],
[0.0, 0.0, 0.0, -4.0 -7.0j, 3.0 +9.0j],
[0.0, 0.0, 0.0, 3.0, -6.0 -3.0j]])
B = np.array([[6.0 +3.0j, -9.0 -2.0j],
[9.0 +5.0j, 7.0 +3.0j],
[3.0 +5.0j, 8.0 -6.0j],
[2.0, 5.0 +7.0j],
[-5.0 +4.0j, -5.0 +9.0j]])
C = np.array([[4.0 +4.0j, -4.0 +9.0j, -8.0 -1.0j, 1.0 +6.0j, -7.0 +6.0j],
[-9.0 -3.0j, -9.0 -9.0j, 6.0 -2.0j, -7.0 -5.0j, -5.0 -5.0j]])
D = np.array([[13.0 +1.0j, -12.0 +1.0j],
[3.0 +2.0j, -5.0 +11.0j]])

sys = self.sysC322 + self.sysC222

np.testing.assert_array_almost_equal(sys.A, A)
np.testing.assert_array_almost_equal(sys.B, B)
np.testing.assert_array_almost_equal(sys.C, C)
np.testing.assert_array_almost_equal(sys.D, D)

def test_subtract_ss_complex(self):
"""Subtract two complex MIMO systems."""

A = np.array([[6.0 +9.0j, 8.0 +9.0j, -4.0 -7.0j, 0.0, 0.0],
[8.0 -7.0j, 3.0, 1.0 -1.0j, 0.0, 0.0],
[-7.0 +9.0j, -8.0 +6.0j, 9.0 +8.0j, 0.0, 0.0],
[0.0, 0.0, 0.0, -4.0 -7.0j, 3.0 +9.0j],
[0.0, 0.0, 0.0, 3.0, -6.0 -3.0j]])
B = np.array([[6.0 +3.0j, -9.0 -2.0j],
[9.0 +5.0j, 7.0 +3.0j],
[3.0 +5.0j, 8.0 -6.0j],
[2.0, 5.0 +7.0j],
[-5.0 +4.0j, -5.0 +9.0j]])
C = np.array([[4.0 +4.0j, -4.0 +9.0j, -8.0 -1.0j, -1.0 -6.0j, 7.0 -6.0j],
[-9.0 -3.0j, -9.0 -9.0j, 6.0 -2.0j, 7.0 +5.0j, 5.0 +5.0j]])
D = np.array([[-3.0 -3.0j, 7.0j],
[9.0 +4.0j, 5.0 -1.0j]])

sys = self.sysC322 - self.sysC222

np.testing.assert_array_almost_equal(sys.A, A)
np.testing.assert_array_almost_equal(sys.B, B)
np.testing.assert_array_almost_equal(sys.C, C)
np.testing.assert_array_almost_equal(sys.D, D)

def test_multiply_ss_complex(self):
"""Multiply two complex MIMO systems."""

A = np.array([[6.0 +9.0j, 8.0 +9.0j, -4.0 -7.0j, 41.0 +98.0j, -25.0 +70.0j],
[8.0 -7.0j, 3.0, 1.0 -1.0j, -55.0 +3.0j, -113.0 -31.0j],
[-7.0 +9.0j, -8.0 +6.0j, 9.0 +8.0j, -113.0 +25.0j, -121.0 -27.0j],
[0.0, 0.0, 0.0, -4.0 -7.0j, 3.0 +9.0j],
[0.0, 0.0, 0.0, 3.0, -6.0 -3.0j]])
B = np.array([[67.0 +51.0j, 30.0 -80.0j],
[44.0 +42.0j, -92.0 -30.0j],
[-16.0 +56.0j, -7.0 +39.0j],
[2.0, 5.0 +7.0j],
[-5.0 +4.0j, -5.0 +9.0j]])
C = np.array([[4.0 +4.0j, -4.0 +9.0j, -8.0 -1.0j, 73.0 +31.0j, 21.0 +47.0j],
[-9.0 -3.0j, -9.0 -9.0j, 6.0 -2.0j, 13.0 +4.0j, -35.0 -10.0j]])
D = np.array([[64.0 -4.0j, -27.0 -65.0j],
[47.0 +21.0j, -57.0 -61.0j]])

sys = self.sysC322 * self.sysC222

np.testing.assert_array_almost_equal(sys.A, A)
np.testing.assert_array_almost_equal(sys.B, B)
np.testing.assert_array_almost_equal(sys.C, C)
np.testing.assert_array_almost_equal(sys.D, D)

def test_evalfr_complex(self):
"""Evaluate the frequency response at one frequency."""

resp = [[4.6799374736968655 -34.9854626345217383j,

@bnavigator bnavigator Mar 2, 2020

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Please be easier on the precision of your reference values. This will almost certainly fail on some build systems due to rounding errors or depending on the used BLAS/LAPACK implementation.

-10.8392352552155344 -10.3778031623880267j],
[28.8313352973005479 +17.1145433776227947j,
5.6628990560933081 +8.8759694583057787j]]

# Correct versions of the call
np.testing.assert_almost_equal(evalfr(self.sysC322, 1j), resp)
np.testing.assert_almost_equal(self.sysC322._evalfr(1.), resp)

# Deprecated version of the call (should generate warning)
import warnings
with warnings.catch_warnings(record=True) as w:
# Set up warnings filter to only show warnings in control module
warnings.filterwarnings("ignore")
warnings.filterwarnings("always", module="control")

# Make sure that we get a pending deprecation warning
self.sysC322.evalfr(1.)
assert len(w) == 1
assert issubclass(w[-1].category, PendingDeprecationWarning)

@unittest.skipIf(not slycot_check(), "slycot not installed")
def test_freq_resp(self):
"""Evaluate the frequency response at multiple frequencies."""

true_mag = [[15.2721178549527039, 3.9176691825112484, 20.5865790875032246],
[24.5384389050919864, 2.8374330975514015, 18.2268344283306227]]

true_phase = [[1.0345533469994428, 2.2133291186570987, 2.6715324185062164],
[1.8217663044282106, -2.8266936088743044, 2.2694910839768196]]
true_omega = [0.1, 10, 0.01j, -1j];

mag, phase, omega = self.sysC623.freqresp(true_omega)

np.testing.assert_almost_equal(mag, true_mag)
np.testing.assert_almost_equal(phase, true_phase)
np.testing.assert_equal(omega, true_omega)


class TestRss(unittest.TestCase):
"""These are tests for the proper functionality of statesp.rss."""

Expand Down