Skip to content

MIMO impulse and step response #514

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

Merged
merged 3 commits into from
Jan 20, 2021
Merged

MIMO impulse and step response #514

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
e838b4b
initial implementation of MIMO step, impulse response
murrayrm Jan 18, 2021
85a480c
finalize MIMO implementation: remove duplicate code, unit test updates
murrayrm Jan 18, 2021
28bbe7f
fixes for @sawyerbfuller review
murrayrm Jan 20, 2021
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
1 change: 0 additions & 1 deletion control/config.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -17,7 +17,6 @@
_control_defaults = {
'control.default_dt': 0,
'control.squeeze_frequency_response': None,
'control.squeeze_time_response': True,
'control.squeeze_time_response': None,
'forced_response.return_x': False,
}
Expand Down
5 changes: 3 additions & 2 deletions control/iosys.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1428,8 +1428,9 @@ def input_output_response(sys, T, U=0., X0=0, params={}, method='RK45',
for i in range(len(T)):
u = U[i] if len(U.shape) == 1 else U[:, i]
y[:, i] = sys._out(T[i], [], u)
return _process_time_response(sys, T, y, [], transpose=transpose,
return_x=return_x, squeeze=squeeze)
return _process_time_response(
sys, T, y, np.array((0, 0, np.asarray(T).size)),
transpose=transpose, return_x=return_x, squeeze=squeeze)

# create X0 if not given, test if X0 has correct shape
X0 = _check_convert_array(X0, [(nstates,), (nstates, 1)],
Expand Down
14 changes: 7 additions & 7 deletions control/tests/matlab2_test.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -121,25 +121,25 @@ def test_step(self, SISO_mats, MIMO_mats, mplcleanup):
#print("gain:", dcgain(sys))

subplot2grid(plot_shape, (0, 0))
t, y = step(sys)
y, t = step(sys)
plot(t, y)

subplot2grid(plot_shape, (0, 1))
T = linspace(0, 2, 100)
X0 = array([1, 1])
t, y = step(sys, T, X0)
y, t = step(sys, T, X0)
plot(t, y)

# Test output of state vector
t, y, x = step(sys, return_x=True)
y, t, x = step(sys, return_x=True)

#Test MIMO system
A, B, C, D = MIMO_mats
sys = ss(A, B, C, D)

subplot2grid(plot_shape, (0, 2))
t, y = step(sys)
plot(t, y)
y, t = step(sys)
plot(t, y[:, 0, 0])

def test_impulse(self, SISO_mats, mplcleanup):
A, B, C, D = SISO_mats
Expand Down Expand Up @@ -168,8 +168,8 @@ def test_impulse_mimo(self, MIMO_mats, mplcleanup):
#Test MIMO system
A, B, C, D = MIMO_mats
sys = ss(A, B, C, D)
t, y = impulse(sys)
plot(t, y, label='MIMO System')
y, t = impulse(sys)
plot(t, y[:, :, 0], label='MIMO System')

legend(loc='best')
#show()
Expand Down
151 changes: 109 additions & 42 deletions control/tests/timeresp_test.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -347,7 +347,7 @@ def test_impulse_response_mimo(self, mimo_ss2):
yref_notrim = np.zeros((2, len(t)))
yref_notrim[:1, :] = yref
_t, yy = impulse_response(sys, T=t, input=0)
np.testing.assert_array_almost_equal(yy, yref_notrim, decimal=4)
np.testing.assert_array_almost_equal(yy[:,0,:], yref_notrim, decimal=4)

@pytest.mark.skipif(StrictVersion(sp.__version__) < "1.3",
reason="requires SciPy 1.3 or greater")
Expand Down Expand Up @@ -639,9 +639,10 @@ def test_time_vector(self, tsystem, fun, squeeze, matarrayout):
if hasattr(tsystem, 't'):
# tout should always match t, which has shape (n, )
np.testing.assert_allclose(tout, tsystem.t)
if squeeze is False or sys.outputs > 1:

if squeeze is False or not sys.issiso():
assert yout.shape[0] == sys.outputs
assert yout.shape[1] == tout.shape[0]
assert yout.shape[-1] == tout.shape[0]
else:
assert yout.shape == tout.shape

Expand Down Expand Up @@ -725,21 +726,22 @@ def test_time_series_data_convention_2D(self, siso_ss1):

@pytest.mark.usefixtures("editsdefaults")
@pytest.mark.parametrize("fcn", [ct.ss, ct.tf, ct.ss2io])
@pytest.mark.parametrize("nstate, nout, ninp, squeeze, shape", [
[1, 1, 1, None, (8,)],
[2, 1, 1, True, (8,)],
[3, 1, 1, False, (1, 8)],
[3, 2, 1, None, (2, 8)],
[4, 2, 1, True, (2, 8)],
[5, 2, 1, False, (2, 8)],
[3, 1, 2, None, (1, 8)],
[4, 1, 2, True, (8,)],
[5, 1, 2, False, (1, 8)],
[4, 2, 2, None, (2, 8)],
[5, 2, 2, True, (2, 8)],
[6, 2, 2, False, (2, 8)],
@pytest.mark.parametrize("nstate, nout, ninp, squeeze, shape1, shape2", [
# state out in squeeze in/out out-only
[1, 1, 1, None, (8,), (8,)],
[2, 1, 1, True, (8,), (8,)],
[3, 1, 1, False, (1, 1, 8), (1, 8)],
[3, 2, 1, None, (2, 1, 8), (2, 8)],
[4, 2, 1, True, (2, 8), (2, 8)],
[5, 2, 1, False, (2, 1, 8), (2, 8)],
[3, 1, 2, None, (1, 2, 8), (1, 8)],
[4, 1, 2, True, (2, 8), (8,)],
[5, 1, 2, False, (1, 2, 8), (1, 8)],
[4, 2, 2, None, (2, 2, 8), (2, 8)],
[5, 2, 2, True, (2, 2, 8), (2, 8)],
[6, 2, 2, False, (2, 2, 8), (2, 8)],
])
def test_squeeze(self, fcn, nstate, nout, ninp, squeeze, shape):
def test_squeeze(self, fcn, nstate, nout, ninp, squeeze, shape1, shape2):
# Figure out if we have SciPy 1+
scipy0 = StrictVersion(sp.__version__) < '1.0'

Expand All @@ -750,27 +752,56 @@ def test_squeeze(self, fcn, nstate, nout, ninp, squeeze, shape):
else:
sys = fcn(ct.rss(nstate, nout, ninp, strictly_proper=True))

# Keep track of expect users warnings
warntype = UserWarning if sys.inputs > 1 else None

# Generate the time and input vectors
tvec = np.linspace(0, 1, 8)
uvec = np.dot(
np.ones((sys.inputs, 1)),
np.reshape(np.sin(tvec), (1, 8)))

#
# Pass squeeze argument and make sure the shape is correct
with pytest.warns(warntype, match="Converting MIMO system"):
_, yvec = ct.impulse_response(sys, tvec, squeeze=squeeze)
assert yvec.shape == shape
#
# For responses that are indexed by the input, check against shape1
# For responses that have no/fixed input, check against shape2
#

_, yvec = ct.initial_response(sys, tvec, 1, squeeze=squeeze)
assert yvec.shape == shape
# Impulse response
if isinstance(sys, StateSpace):
# Check the states as well
_, yvec, xvec = ct.impulse_response(
sys, tvec, squeeze=squeeze, return_x=True)
if sys.issiso():
assert xvec.shape == (sys.states, 8)
else:
assert xvec.shape == (sys.states, sys.inputs, 8)
else:
_, yvec = ct.impulse_response(sys, tvec, squeeze=squeeze)
assert yvec.shape == shape1

with pytest.warns(warntype, match="Converting MIMO system"):
# Step response
if isinstance(sys, StateSpace):
# Check the states as well
_, yvec, xvec = ct.step_response(
sys, tvec, squeeze=squeeze, return_x=True)
if sys.issiso():
assert xvec.shape == (sys.states, 8)
else:
assert xvec.shape == (sys.states, sys.inputs, 8)
else:
_, yvec = ct.step_response(sys, tvec, squeeze=squeeze)
assert yvec.shape == shape
assert yvec.shape == shape1

# Initial response (only indexed by output)
if isinstance(sys, StateSpace):
# Check the states as well
_, yvec, xvec = ct.initial_response(
sys, tvec, 1, squeeze=squeeze, return_x=True)
assert xvec.shape == (sys.states, 8)
else:
_, yvec = ct.initial_response(sys, tvec, 1, squeeze=squeeze)
assert yvec.shape == shape2

# Forced response (only indexed by output)
if isinstance(sys, StateSpace):
# Check the states as well
_, yvec, xvec = ct.forced_response(
Expand All @@ -779,52 +810,54 @@ def test_squeeze(self, fcn, nstate, nout, ninp, squeeze, shape):
else:
# Just check the input/output response
_, yvec = ct.forced_response(sys, tvec, uvec, 0, squeeze=squeeze)
assert yvec.shape == shape
assert yvec.shape == shape2

# Test cases where we choose a subset of inputs and outputs
_, yvec = ct.step_response(
sys, tvec, input=ninp-1, output=nout-1, squeeze=squeeze)
# Possible code if we implemenet a squeeze='siso' option
if squeeze is False:
# Shape should be unsqueezed
assert yvec.shape == (1, 8)
assert yvec.shape == (1, 1, 8)
else:
# Shape should be squeezed
assert yvec.shape == (8, )

# For InputOutputSystems, also test input_output_response
# For InputOutputSystems, also test input/output response
if isinstance(sys, ct.InputOutputSystem) and not scipy0:
_, yvec = ct.input_output_response(sys, tvec, uvec, squeeze=squeeze)
assert yvec.shape == shape
assert yvec.shape == shape2

#
# Changing config.default to False should return 3D frequency response
#
ct.config.set_defaults('control', squeeze_time_response=False)

with pytest.warns(warntype, match="Converting MIMO system"):
_, yvec = ct.impulse_response(sys, tvec)
assert yvec.shape == (sys.outputs, 8)
_, yvec = ct.impulse_response(sys, tvec)
if squeeze is not True or sys.inputs > 1 or sys.outputs > 1:
assert yvec.shape == (sys.outputs, sys.inputs, 8)

_, yvec = ct.initial_response(sys, tvec, 1)
assert yvec.shape == (sys.outputs, 8)
_, yvec = ct.step_response(sys, tvec)
if squeeze is not True or sys.inputs > 1 or sys.outputs > 1:
assert yvec.shape == (sys.outputs, sys.inputs, 8)

with pytest.warns(warntype, match="Converting MIMO system"):
_, yvec = ct.step_response(sys, tvec)
assert yvec.shape == (sys.outputs, 8)
_, yvec = ct.initial_response(sys, tvec, 1)
if squeeze is not True or sys.outputs > 1:
assert yvec.shape == (sys.outputs, 8)

if isinstance(sys, ct.StateSpace):
_, yvec, xvec = ct.forced_response(
sys, tvec, uvec, 0, return_x=True)
assert xvec.shape == (sys.states, 8)
else:
_, yvec = ct.forced_response(sys, tvec, uvec, 0)
assert yvec.shape == (sys.outputs, 8)
if squeeze is not True or sys.outputs > 1:
assert yvec.shape == (sys.outputs, 8)

# For InputOutputSystems, also test input_output_response
if isinstance(sys, ct.InputOutputSystem) and not scipy0:
_, yvec = ct.input_output_response(sys, tvec, uvec)
assert yvec.shape == (sys.noutputs, 8)
if squeeze is not True or sys.outputs > 1:
assert yvec.shape == (sys.outputs, 8)

@pytest.mark.parametrize("fcn", [ct.ss, ct.tf, ct.ss2io])
def test_squeeze_exception(self, fcn):
Expand Down Expand Up @@ -861,3 +894,37 @@ def test_squeeze_0_8_4(self, nstate, nout, ninp, squeeze, shape):

_, yvec = ct.initial_response(sys, tvec, 1, squeeze=squeeze)
assert yvec.shape == shape

@pytest.mark.parametrize(
"nstate, nout, ninp, squeeze, ysh_in, ysh_no, xsh_in", [
[4, 1, 1, None, (8,), (8,), (8, 4)],
[4, 1, 1, True, (8,), (8,), (8, 4)],
[4, 1, 1, False, (8, 1, 1), (8, 1), (8, 4)],
[4, 2, 1, None, (8, 2, 1), (8, 2), (8, 4, 1)],
[4, 2, 1, True, (8, 2), (8, 2), (8, 4, 1)],
[4, 2, 1, False, (8, 2, 1), (8, 2), (8, 4, 1)],
[4, 1, 2, None, (8, 1, 2), (8, 1), (8, 4, 2)],
[4, 1, 2, True, (8, 2), (8,), (8, 4, 2)],
[4, 1, 2, False, (8, 1, 2), (8, 1), (8, 4, 2)],
[4, 2, 2, None, (8, 2, 2), (8, 2), (8, 4, 2)],
[4, 2, 2, True, (8, 2, 2), (8, 2), (8, 4, 2)],
[4, 2, 2, False, (8, 2, 2), (8, 2), (8, 4, 2)],
])
def test_response_transpose(
self, nstate, nout, ninp, squeeze, ysh_in, ysh_no, xsh_in):
sys = ct.rss(nstate, nout, ninp)
T = np.linspace(0, 1, 8)

# Step response - input indexed
t, y, x = ct.step_response(
sys, T, transpose=True, return_x=True, squeeze=squeeze)
assert t.shape == (T.size, )
assert y.shape == ysh_in
assert x.shape == xsh_in

# Initial response - no input indexing
t, y, x = ct.initial_response(
sys, T, 1, transpose=True, return_x=True, squeeze=squeeze)
assert t.shape == (T.size, )
assert y.shape == ysh_no
assert x.shape == (T.size, sys.states)
Loading