Movatterモバイル変換

[0]ホーム
URL:

Skip to content

Navigation Menu

Sign in
Appearance settings

Search code, repositories, users, issues, pull requests...

Provide feedback

We read every piece of feedback, and take your input very seriously.

Saved searches

Use saved searches to filter your results more quickly

 Sign up
Appearance settings

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.

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

Already on GitHub?Sign in to your account

Closed
lytex wants to merge12 commits intopython-control:mainfromlytex:complex
Closed
Show file tree
Hide file tree
Changes fromall commits
Commits
Show all changes
12 commits
Select commitHold shift + click to select a range
e6f2f60
allow _ssmatrix to accept complex values
lytexFeb 9, 2020
53dbfe4
add new tests with complex state matrices
lytexFeb 29, 2020
05d730a
rename to test_evalfr_complex, use self.sysC332
lytexFeb 29, 2020
a55165a
missed one sys.evalfr in previous commit
lytexFeb 29, 2020
2d231a8
missed one sys.evalfr in previous commit
lytexFeb 29, 2020
f1c9b9e
add complex to test_freq_resp
lytexMar 3, 2020
e7a93ea
remove warnings settings for debugging
lytexMar 3, 2020
bb44827
ease precision of reference values in test
lytexMar 3, 2020
fc5c76b
Merge branch 'master' into complex
lytexMar 30, 2020
5f58723
build an auxiliary system to compare against
lytexMar 30, 2020
1803d38
build an auxiliary system to compare against
lytexMar 30, 2020
84bc8ef
add code to use slycot complex function ab08nz if any matrix is complex
lytexApr 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
23 changes: 17 additions & 6 deletionscontrol/statesp.py
View file
Open in desktop
Original file line numberDiff line numberDiff 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
Copy link
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 usematrix() andasarray() withoutdtype and figure out thedtype later?

Copy link
Contributor

@bnavigatorbnavigatorJul 29, 2020•
edited
Loading

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 beint 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 usedtype=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 DownExpand Up@@ -572,10 +578,15 @@ def zero(self):
# Use AB08ND from Slycot if it's available, otherwise use
# scipy.lingalg.eigvals().
try:
from slycot import ab08nd

out = ab08nd(self.A.shape[0], self.B.shape[1], self.C.shape[0],
self.A, self.B, self.C, self.D)
if np.isrealobj(self.A) and np.isrealobj(self.B) and\
np.isrealobj(self.C) and np.isrealobj(self.D):
from slycot import ab08nd
out = ab08nd(self.A.shape[0], self.B.shape[1], self.C.shape[0],
self.A, self.B, self.C, self.D)
else:
from slycot import ab08nz
out = ab08nz(self.A.shape[0], self.B.shape[1], self.C.shape[0],
self.A, self.B, self.C, self.D)
nu = out[0]
if nu == 0:
return np.array([])
Expand Down
220 changes: 218 additions & 2 deletionscontrol/tests/statesp_test.py
View file
Open in desktop
Original file line numberDiff line numberDiff line change
Expand Up@@ -7,13 +7,12 @@
import numpy as np
from numpy.linalg import solve
from scipy.linalg import eigvals, block_diag
from control import matlab
from control import matlab, series
from control.statesp import StateSpace, _convertToStateSpace, tf2ss
from control.xferfcn import TransferFunction, ss2tf
from control.lti import evalfr
from control.exception import slycot_check


class TestStateSpace(unittest.TestCase):
"""Tests for the StateSpace class."""

Expand DownExpand Up@@ -63,6 +62,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 DownExpand Up@@ -519,6 +562,179 @@ 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.493759 + 8.073864j,
-4.7860079 + 29.3266582j,
-7.4509692 +5.5262915j])
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.56850055,
3.3685005])
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_aux = StateSpace(A, B, C, D)

sys = series(self.sysC322, self.sysC222)

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

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

resp = [[4.6799374 - 34.9854626j,
-10.8392352 - 10.3778031j],
[28.8313352 + 17.1145433j,
5.6628990 + 8.8759694j]]

# 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_complex(self):
"""Evaluate the frequency response at multiple frequencies."""

true_mag = [[15.2721178, 3.9176691, 20.5865790],
[24.5384389, 2.8374330, 18.2268344]]

true_phase = [[1.03455334, 2.2133291, 2.6715324],
[1.8217663, -2.8266936, 2.2694910]]
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
[8]ページ先頭
©2009-2025 Movatter.jp