Skip to content

Commit ae43212

Browse files
murrayrmmurrayrm
committed
add input/output removal + unit tests
1 parent 77b25d5 commit ae43212

3 files changed

Lines changed: 77 additions & 29 deletions

File tree

‎control/modelsimp.py‎

Lines changed: 24 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -138,20 +138,33 @@ def model_reduction(
138138
# Check system is stable
139139
if warn_unstable:
140140
if isctime(sys) and np.any(np.linalg.eigvals(sys.A).real >= 0.0) or \
141-
np.any(np.abs(np.linalg.eigvals(sys.A)) >= 1):
141+
isdtime(sys) and np.any(np.abs(np.linalg.eigvals(sys.A)) >= 1):
142142
warnings.warn("System is unstable; reduction may be meaningless")
143143

144144
# Utility function to process keep/elim keywords
145-
def _process_elim_or_keep(elim, keep, labels):
146-
elim = np.sort(elim).tolist()
147-
return elim, [i for i in range(len(labels)) if i not in elim]
145+
def _process_elim_or_keep(elim, keep, labels, allow_both=False):
146+
elim = [] if elim is None else np.atleast_1d(elim)
147+
keep = [] if keep is None else np.atleast_1d(keep)
148+
149+
if len(elim) > 0 and len(keep) > 0:
150+
if not allow_both:
151+
raise ValueError(
152+
"can't provide both 'keep' and 'elim' for same variables")
153+
elif len(keep) > 0:
154+
keep = np.sort(keep).tolist()
155+
elim = [i for i in range(len(labels)) if i not in keep]
156+
else:
157+
elim = [] if elim is None else np.sort(elim).tolist()
158+
keep = [i for i in range(len(labels)) if i not in elim]
159+
return elim, keep
148160

149161
# Determine which states to keep
150162
elim_states, keep_states = _process_elim_or_keep(
151163
elim_states, keep_states, sys.state_labels)
152-
153-
keep_inputs = slice(None, None)
154-
keep_outputs = slice(None, None)
164+
elim_inputs, keep_inputs = _process_elim_or_keep(
165+
elim_inputs, keep_inputs, sys.input_labels)
166+
elim_outputs, keep_outputs = _process_elim_or_keep(
167+
elim_outputs, keep_outputs, sys.output_labels)
155168

156169
# Create submatrix of states we are keeping
157170
A11 = sys.A[:, keep_states][keep_states, :] # states we are keeping
@@ -166,11 +179,11 @@ def _process_elim_or_keep(elim, keep, labels):
166179
C2 = sys.C[:, elim_states]
167180

168181
# Figure out the new state space system
169-
if method == 'matchdc':
182+
if method == 'matchdc' and A22.size > 0:
170183
if sys.isdtime(strict=True):
171184
raise NotImplementedError(
172185
"'matchdc' not (yet) supported for discrete time systems")
173-
186+
174187
# if matchdc, residualize
175188
# Check if the matrix A22 is invertible
176189
if np.linalg.matrix_rank(A22) != len(elim_states):
@@ -190,8 +203,8 @@ def _process_elim_or_keep(elim, keep, labels):
190203
Cr = C1 - C2 @ A22I_A21
191204
Dr = sys.D - C2 @ A22I_B2
192205

193-
elif method == 'truncate':
194-
# if truncate, simply discard state x2
206+
elif method == 'truncate' or A22.size == 0:
207+
# Get rid of unwanted states
195208
Ar = A11
196209
Br = B1
197210
Cr = C1

‎control/tests/docstrings_test.py‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -39,7 +39,7 @@
3939
control.ss2tf: '48ff25d22d28e7b396e686dd5eb58831',
4040
control.tf: '53a13f4a7f75a31c81800e10c88730ef',
4141
control.tf2ss: '086a3692659b7321c2af126f79f4bc11',
42-
control.markov: 'eda7c4635bbb863ae6659e574285d356',
42+
control.markov: 'a4199c54cb50f07c0163d3790739eafe',
4343
control.gangof4: '0e52eb6cf7ce024f9a41f3ae3ebf04f7',
4444
}
4545

‎control/tests/modelsimp_test.py‎

Lines changed: 52 additions & 17 deletions
Original file line numberDiff line numberDiff line change
@@ -3,6 +3,7 @@
33
RMM, 30 Mar 2011 (based on TestModelSimp from v0.4a)
44
"""
55

6+
import math
67
import warnings
78

89
import numpy as np
@@ -45,7 +46,7 @@ def testMarkovSignature(self):
4546
inputs=U,
4647
input_labels='u',
4748
)
48-
49+
4950
# setup
5051
m = 3
5152
Htrue = np.array([1., 0., 0.])
@@ -104,7 +105,6 @@ def testMarkovSignature(self):
104105
HT = markov(response, m)
105106
np.testing.assert_array_almost_equal(HT, np.transpose(Htrue))
106107
response.transpose=False
107-
108108

109109
# Test example from docstring
110110
# TODO: There is a problem here, last markov parameter does not fit
@@ -215,7 +215,7 @@ def testMarkovResults(self, k, m, n):
215215
Mtrue = np.hstack([Hd.D] + [
216216
Hd.C @ np.linalg.matrix_power(Hd.A, i) @ Hd.B
217217
for i in range(m-1)])
218-
218+
219219
Mtrue = np.squeeze(Mtrue)
220220

221221
# Generate input/output data
@@ -241,8 +241,8 @@ def testMarkovResults(self, k, m, n):
241241
Mcomp_scaled = markov(response, m, dt=Ts)
242242

243243
np.testing.assert_allclose(Mtrue, Mcomp, rtol=1e-6, atol=1e-8)
244-
np.testing.assert_allclose(Mtrue_scaled, Mcomp_scaled, rtol=1e-6, atol=1e-8)
245-
244+
np.testing.assert_allclose(
245+
Mtrue_scaled, Mcomp_scaled, rtol=1e-6, atol=1e-8)
246246

247247
def testERASignature(self):
248248

@@ -256,7 +256,7 @@ def testERASignature(self):
256256
B = np.array([[1.],[1.,]])
257257
C = np.array([[1., 0.,]])
258258
D = np.array([[0.,]])
259-
259+
260260
T = np.arange(0,10,1)
261261
sysd_true = StateSpace(A,B,C,D,True)
262262
ir_true = impulse_response(sysd_true,T=T)
@@ -265,15 +265,15 @@ def testERASignature(self):
265265
sysd_est, _ = eigensys_realization(ir_true,r=2)
266266
ir_est = impulse_response(sysd_est, T=T)
267267
_, H_est = ir_est
268-
268+
269269
np.testing.assert_allclose(H_true, H_est, rtol=1e-6, atol=1e-8)
270270

271271
# test ndarray
272272
_, YY_true = ir_true
273273
sysd_est, _ = eigensys_realization(YY_true,r=2)
274274
ir_est = impulse_response(sysd_est, T=T)
275275
_, H_est = ir_est
276-
276+
277277
np.testing.assert_allclose(H_true, H_est, rtol=1e-6, atol=1e-8)
278278

279279
# test mimo
@@ -307,18 +307,18 @@ def testERASignature(self):
307307
step_true = step_response(sysd_true)
308308
step_est = step_response(sysd_est)
309309

310-
np.testing.assert_allclose(step_true.outputs,
310+
np.testing.assert_allclose(step_true.outputs,
311311
step_est.outputs,
312312
rtol=1e-6, atol=1e-8)
313-
313+
314314
# test ndarray
315315
_, YY_true = ir_true
316316
sysd_est, _ = eigensys_realization(YY_true,r=4,dt=dt)
317317

318318
step_true = step_response(sysd_true, T=T)
319319
step_est = step_response(sysd_est, T=T)
320320

321-
np.testing.assert_allclose(step_true.outputs,
321+
np.testing.assert_allclose(step_true.outputs,
322322
step_est.outputs,
323323
rtol=1e-6, atol=1e-8)
324324

@@ -403,7 +403,7 @@ def testBalredTruncate(self):
403403
B = np.array([[2.], [0.], [0.], [0.]])
404404
C = np.array([[0.5, 0.6875, 0.7031, 0.5]])
405405
D = np.array([[0.]])
406-
406+
407407
sys = StateSpace(A, B, C, D)
408408
orders = 2
409409
rsys = balred(sys, orders, method='truncate')
@@ -424,7 +424,7 @@ def testBalredTruncate(self):
424424
# Apply a similarity transformation
425425
Ar, Br, Cr = T @ Ar @ T, T @ Br, Cr @ T
426426
break
427-
427+
428428
# Make sure we got the correct answer
429429
np.testing.assert_array_almost_equal(Ar, Artrue, decimal=2)
430430
np.testing.assert_array_almost_equal(Br, Brtrue, decimal=4)
@@ -444,20 +444,20 @@ def testBalredMatchDC(self):
444444
B = np.array([[2.], [0.], [0.], [0.]])
445445
C = np.array([[0.5, 0.6875, 0.7031, 0.5]])
446446
D = np.array([[0.]])
447-
447+
448448
sys = StateSpace(A, B, C, D)
449449
orders = 2
450450
rsys = balred(sys,orders,method='matchdc')
451451
Ar, Br, Cr, Dr = rsys.A, rsys.B, rsys.C, rsys.D
452-
452+
453453
# Result from MATLAB
454454
Artrue = np.array(
455455
[[-4.43094773, -4.55232904],
456456
[-4.55232904, -5.36195206]])
457457
Brtrue = np.array([[1.36235673], [1.03114388]])
458458
Crtrue = np.array([[1.36235673, 1.03114388]])
459459
Drtrue = np.array([[-0.08383902]])
460-
460+
461461
# Look for possible changes in state in slycot
462462
T1 = np.array([[1, 0], [0, -1]])
463463
T2 = np.array([[-1, 0], [0, 1]])
@@ -467,9 +467,44 @@ def testBalredMatchDC(self):
467467
# Apply a similarity transformation
468468
Ar, Br, Cr = T @ Ar @ T, T @ Br, Cr @ T
469469
break
470-
470+
471471
# Make sure we got the correct answer
472472
np.testing.assert_array_almost_equal(Ar, Artrue, decimal=2)
473473
np.testing.assert_array_almost_equal(Br, Brtrue, decimal=4)
474474
np.testing.assert_array_almost_equal(Cr, Crtrue, decimal=4)
475475
np.testing.assert_array_almost_equal(Dr, Drtrue, decimal=4)
476+
477+
478+
@pytest.mark.parametrize("kwargs, nstates, noutputs, ninputs", [
479+
({'elim_states': [1, 3]}, 3, 3, 3),
480+
({'elim_inputs': [1, 2], 'keep_states': [1, 3]}, 2, 3, 1),
481+
({'elim_outputs': [1, 2], 'keep_inputs': [0, 1],}, 5, 1, 2),
482+
({'keep_states': [2, 0], 'keep_outputs': [0, 1]}, 2, 2, 3),
483+
({'elim_inputs': [0, 1, 2]}, 5, 3, 0), # no inputs
484+
({'elim_outputs': [0, 1, 2]}, 5, 0, 3), # no outputs
485+
({'elim_states': [0, 1, 2, 3, 4]}, 0, 3, 3), # no states
486+
({'elim_states': [0, 1], 'keep_states': [1, 2]}, None, None, None)
487+
])
488+
@pytest.mark.parametrize("method", ['truncate', 'matchdc'])
489+
def test_model_reduction(method, kwargs, nstates, noutputs, ninputs):
490+
sys = ct.rss(5, 3, 3)
491+
492+
if nstates is None:
493+
# Arguments should generate an error
494+
with pytest.raises(ValueError, match="can't provide both"):
495+
red = ct.model_reduction(sys, **kwargs, method=method)
496+
return
497+
else:
498+
red = ct.model_reduction(sys, **kwargs, method=method)
499+
500+
assert red.nstates == nstates
501+
assert red.ninputs == ninputs
502+
assert red.noutputs == noutputs
503+
504+
if method == 'matchdc':
505+
# Define a new system with truncated inputs and outputs
506+
# (assumes we always keep the initial inputs and outputs)
507+
chk = ct.ss(
508+
sys.A, sys.B[:, :ninputs], sys.C[:noutputs, :],
509+
sys.D[:noutputs, :][:, :ninputs])
510+
np.testing.assert_allclose(red(0), chk(0))

0 commit comments

Comments
 (0)