8000 Change API to work with TimeResponseData and ndarray as input, add py… · python-control/python-control@614a080 · GitHub
[go: up one dir, main page]
Skip to content

Commit 614a080

Browse files
KybernetikJoKybernetikJo
committed
Change API to work with TimeResponseData and ndarray as input, add pytest, fix small things
1 parent 6bbad5f commit 614a080

File tree

3 files changed

+128
-16
lines changed

3 files changed

+128
-16
lines changed

‎control/modelsimp.py

Lines changed: 38 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -48,6 +48,7 @@
4848
from .iosys import isdtime, isctime
4949
from .statesp import StateSpace
5050
from .statefbk import gram
51+
from .timeresp import TimeResponseData
5152

5253
__all__ = ['hsvd', 'balred', 'modred', 'era', 'markov', 'minreal']
5354

@@ -368,8 +369,10 @@ def minreal(sys, tol=None, verbose=True):
368369
return sysr
369370

370371

371-
def era(data, r, m=None, n=None, dt=True):
372-
"""Calculate an ERA model of order `r` based on the impulse-response data.
372+
def era(arg, r, m=None, n=None, dt=True, transpose=False):
373+
r"""era(YY, r)
374+
375+
Calculate an ERA model of order `r` based on the impulse-response data.
373376
374377
This function computes a discrete time system
375378
@@ -380,8 +383,19 @@ def era(data, r, m=None, n=None, dt=True):
380383
381384
for a given impulse-response data (see [1]_).
382385
386+
The function can be called with 2 arguments:
387+
388+
* ``sysd, S = era(data, r)``
389+
* ``sysd, S = era(YY, r)``
390+
391+
where `response` is an `TimeResponseData` object, and `YY` is 1D or 3D
392+
array and r is an integer.
393+
383394
Parameters
384395
----------
396+
YY : array_like
397+
impulse-response data from which the StateSpace model is estimated,
398+
1D or 3D array.
385399
data : TimeResponseData
386400
impulse-response data from which the StateSpace model is estimated.
387401
r : integer
@@ -398,11 +412,16 @@ def era(data, r, m=None, n=None, dt=True):
398412
It can be used to scale the StateSpace model in order to match
399413
the impulse response of this library.
400414
Default values is True.
415+
transpose : bool, optional
416+
Assume that input data is transposed relative to the standard
417+
:ref:`time-series-convention`. For TimeResponseData this parameter
418+
is ignored.
419+
Default value is False.
401420
402421
Returns
403422
-------
404423
sys : StateSpace
405-
A reduced order model sys=StateSpace(Ar,Br,Cr,Dr,dt)
424+
A reduced order model sys=StateSpace(Ar,Br,Cr,Dr,dt).
406425
S : array
407426
Singular values of Hankel matrix.
408427
Can be used to choose a good r value.
@@ -416,6 +435,10 @@ def era(data, r, m=None, n=None, dt=True):
416435
417436
Examples
418437
--------
438+
>>> T = np.linspace(0, 10, 100)
439+
>>> _, YY = ct.impulse_response(ct.tf([1], [1, 0.5], True), T)
440+
>>> sysd, _ = ct.era(YY, r=1)
441+
419442
>>> T = np.linspace(0, 10, 100)
420443
>>> response = ct.impulse_response(ct.tf([1], [1, 0.5], True), T)
421444
>>> sysd, _ = ct.era(response, r=1)
@@ -434,10 +457,16 @@ def block_hankel_matrix(Y, m, n):
434457

435458
return H
436459

437-
Y = np.array(data.outputs, ndmin=3)
438-
if data.transpose:
439-
Y = np.transpose(Y)
440-
q, p, l = Y.shape
460+
if isinstance(arg, TimeResponseData):
461+
YY = np.array(arg.outputs, ndmin=3)
462+
if arg.transpose:
463+
YY = np.transpose(YY)
464+
else:
465+
YY = np.array(arg, ndmin=3)
466+
if transpose:
467+
YY = np.transpose(YY)
468+
469+
q, p, l = YY.shape
441470

442471
if m is None:
443472
m = 2*r
@@ -450,7 +479,7 @@ def block_hankel_matrix(Y, m, n):
450479
if (l-1) < m+n:
451480
raise ValueError("Not enough data for requested number of parameters")
452481

453-
H = block_hankel_matrix(Y[:,:,1:], m, n+1) # Hankel matrix (q*m, p*(n+1))
482+
H = block_hankel_matrix(YY[:,:,1:], m, n+1) # Hankel matrix (q*m, p*(n+1))
454483
Hf = H[:,:-p] # first p*n columns of H
455484
Hl = H[:,p:] # last p*n columns of H
456485

@@ -463,7 +492,7 @@ def block_hankel_matrix(Y, m, n):
463492
Ar = Sigma_inv @ Ur.T @ Hl @ Vhr.T @ Sigma_inv
464493
Br = Sigma_inv @ Ur.T @ Hf[:,0:p]*dt
465494
Cr = Hf[0:q,:] @ Vhr.T @ Sigma_inv
466-
Dr = Y[:,:,0]
495+
Dr = YY[:,:,0]
467496

468497
return StateSpace(Ar,Br,Cr,Dr,dt), S
469498

‎control/tests/modelsimp_test.py

Lines changed: 81 additions & 2 deletions
166
Original file line numberDiff line numberDiff line change
@@ -7,10 +7,10 @@
77
import pytest
88

99

10-
from control import StateSpace, forced_response, tf, rss, c2d
10+
from control import StateSpace, impulse_response, step_response, forced_response, tf, rss, c2d
1111
from control.exception import ControlMIMONotImplemented
1212
from control.tests.conftest import slycotonly
13-
from control.modelsimp import balred, hsvd, markov, modred
13+
from control.modelsimp import balred, hsvd, markov, modred, era
1414

1515

1616
class TestModelsimp:
@@ -111,6 +111,85 @@ def testMarkovResults(self, k, m, n):
111111
# for k=5, m=n=10: 0.015 %
112112
np.testing.assert_allclose(Mtrue, Mcomp, rtol=1e-6, atol=1e-8)
113113

114+
def testERASignature(self):
115+
116+
# test siso
117+
# Katayama, Subspace Methods for System Identification
118+
# Example 6.1, Fibonacci sequence
119+
H_true = np.array([0.,1.,1.,2.,3.,5.,8.,13.,21.,34.])
120+
121+
# A realization of fibonacci impulse response
122+
A = np.array([[0., 1.],[1., 1.,]])
123+
B = np.array([[1.],[1.,]])
124+
C = np.array([[1., 0.,]])
125+
D = np.array([[0.,]])
126+
127+
T = np.arange(0,10,1)
128+
sysd_true = StateSpace(A,B,C,D,True)
129+
ir_true = impulse_response(sysd_true,T=T)
130+
131+
# test TimeResponseData
132+
sysd_est, _ = era(ir_true,r=2)
133+
ir_est = impulse_response(sysd_est, T=T)
134+
_, H_est = ir_est
135+
136+
np.testing.assert_allclose(H_true, H_est, rtol=1e-6, atol=1e-8)
137+
138+
# test ndarray
139+
_, YY_true = ir_true
140+
sysd_est, _ = era(YY_true,r=2)
141+
ir_est = impulse_response(sysd_est, T=T)
142+
_, H_est = ir_est
143+
144+
np.testing.assert_allclose(H_true, H_est, rtol=1e-6, atol=1e-8)
145+
146+
# test mimo
147+
# Mechanical Vibrations: Theory and Application, SI Edition, 1st ed.
148+
# Figure 6.5 / Example 6.7
149+
# m q_dd + c q_d + k q = f
150+
m1, k1, c1 = 1., 4., 1.
151+
m2, k2, c2 = 2., 2., 1.
152+
k3, c3 = 6., 2.
153+
154+
A = np.array([
155+
[0., 0., 1., 0.],
156+
[0., 0., 0., 1.],
157+
[-(k1+k2)/m1, (k2)/m1, -(c1+c2)/m1, c2/m1],
158+
[(k2)/m2, -(k2+k3)/m2, c2/m2, -(c2+c3)/m2]
159+
])
160+
B = np.array([[0.,0.],[0.,0.],[1/m1,0.],[0.,1/m2]])
161+
C = np.array([[1.0, 0.0, 0.0, 0.0],[0.0, 1.0, 0.0, 0.0]])
162+
D = np.zeros((2,2))
163+
164+
sys = StateSpace(A, B, C, D)
165+
+
dt = 0.1
167+
T = np.arange(0,10,dt)
168+
sysd_true = sys.sample(dt, method='zoh')
169+
ir_true = impulse_response(sysd_true, T=T)
170+
171+
# test TimeResponseData
172+
sysd_est, _ = era(ir_true,r=4,dt=dt)
173+
174+
step_true = step_response(sysd_true)
175+
step_est = step_response(sysd_est)
176+
177+
np.testing.assert_allclose(step_true.outputs,
178+
step_est.outputs,
179+
rtol=1e-6, atol=1e-8)
180+
181+
# test ndarray
182+
_, YY_true = ir_true
183+
sysd_est, _ = era(YY_true,r=4,dt=dt)
184+
185+
step_true = step_response(sysd_true, T=T)
186+
step_est = step_response(sysd_est, T=T)
187+
188+
np.testing.assert_allclose(step_true.outputs,
189+
step_est.outputs,
190+
rtol=1e-6, atol=1e-8)
191+
192+
114193
def testModredMatchDC(self):
115194
#balanced realization computed in matlab for the transfer function:
116195
# num = [1 11 45 32], den = [1 15 60 200 60]

‎examples/era_msd.py

Lines changed: 9 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -11,10 +11,12 @@
1111
import control as ct
1212

1313
# set up a mass spring damper system (2dof, MIMO case)
14-
# m q_dd + c q_d + k q = u
15-
m1, k1, c1 = 1., 1., .1
16-
m2, k2, c2 = 2., .5, .1
17-
k3, c3 = .5, .1
14+
# Mechanical Vibrations: Theory and Application, SI Edition, 1st ed.
15+
# Figure 6.5 / Example 6.7
16+
# m q_dd + c q_d + k q = f
17+
m1, k1, c1 = 1., 4., 1.
18+
m2, k2, c2 = 2., 2., 1.
19+
k3, c3 = 6., 2.
1820

1921
A = np.array([
2022
[0., 0., 1., 0.],
@@ -39,7 +41,7 @@
3941
sys = ct.StateSpace(A, B, C, D)
4042

4143

42-
dt = 0.5
44+
dt = 0.1
4345
sysd = sys.sample(dt, method='zoh')
4446
response = ct.impulse_response(sysd)
4547
response.plot()
@@ -48,7 +50,9 @@
4850
sysd_est, _ = ct.era(response,r=4,dt=dt)
4951

5052
step_true = ct.step_response(sysd)
53+
step_true.sysname="H_true"
5154
step_est = ct.step_response(sysd_est)
55+
step_est.sysname="H_est"
5256

5357
step_true.plot(title=xixo)
5458
step_est.plot(color='orange',linestyle='dashed')

0 commit comments

Comments
 (0)
0