8000 Improve markov function, add mimo support, change api to TimeResponse… · python-control/python-control@ec38f2e · GitHub
[go: up one dir, main page]
Skip to content

Commit ec38f2e

Browse files
KybernetikJoKybernetikJo
committed
Improve markov function, add mimo support, change api to TimeResponseData
1 parent 4acc78b commit ec38f2e

File tree

2 files changed

+133
-103
lines changed

2 files changed

+133
-103
lines changed

‎control/modelsimp.py

Lines changed: 103 additions & 83 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 .st 10000 atefbk import gram
51+
from .timeresp import TimeResponseData
5152

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

@@ -402,9 +403,9 @@ def era(YY, m, n, nin, nout, r):
402403
raise NotImplementedError('This function is not implemented yet.')
403404

404405

405-
def markov(Y, U, m=None, transpose=False):
406+
def markov(data, m=None, dt=True, truncate=False):
406407
"""Calculate the first `m` Markov parameters [D CB CAB ...]
407-
from input `U`, output `Y`.
408+
from data
408409
409410
This function computes the Markov parameters for a discrete time system
410411
@@ -420,23 +421,45 @@ def markov(Y, U, m=None, transpose=False):
420421
421422
Parameters
422423
----------
423-
Y : array_like
424-
Output data. If the array is 1D, the system is assumed to be single
425-
input. If the array is 2D and transpose=False, the columns of `Y`
426-
are taken as time points, otherwise the rows of `Y` are taken as
427-
time points.
428-
U : array_like
429-
Input data, arranged in the same way as `Y`.
424+
data : TimeResponseData
425+
Response data from which the Markov parameters where estimated.
426+
Input and output data must be 1D or 2D array.
430427
m : int, optional
431428
Number of Markov parameters to output. Defaults to len(U).
432-
transpose : bool, optional
433-
Assume that input data is transposed relative to the standard
434-
:ref:`time-series-convention`. Default value is False.
429+
dt : (True of float, optional)
430+
True indicates discrete time with unspecified sampling time,
431+
positive number is discrete time with specified sampling time.
432+
It can be used to scale the markov parameters in order to match
433+
the impulse response of this library.
434+
Default values is True.
435+
truncate : bool, optional
436+
Do not use first m equation for least least squares.
437+
Default value is False.
435438
< 10000 span class=pl-s>
436439
Returns
437440
-------
438-
H : ndarray
439-
First m Markov parameters, [D CB CAB ...]
441+
H : TimeResponseData
442+
Markov parameters / impulse response [D CB CAB ...] represented as
443+
a :class:`TimeResponseData` object containing the following properties:
444+
445+
* time (array): Time values of the output.
446+
447+
* outputs (array): Response of the system. If the system is SISO,
448+
the array is 1D (indexed by time). If the
449+
system is not SISO, the array is 3D (indexed
450+
by the output, trace, and time).
451+
452+
* inputs (array): Inputs of the system. If the system is SISO,
453+
the array is 1D (indexed by time). If the
454+
system is not SISO, the array is 3D (indexed
455+
by the output, trace, and time).
456+
457+
Notes
458+
-----
459+
It works for SISO and MIMO systems.
460+
461+
This function does comply with the Python Control Library
462+
:ref:`time-series-convention` for representation of time series data.
440463
441464
References
442465
----------
@@ -445,95 +468,69 @@ def markov(Y, U, m=None, transpose=False):
445468
and experiments. Journal of Guidance Control and Dynamics, 16(2),
446469
320-329, 2012. http://doi.org/10.2514/3.21006
447470
448-
Notes
449-
-----
450-
Currently only works for SISO systems.
451-
452-
This function does not currently comply with the Python Control Library
453-
:ref:`time-series-convention` for representation of time series data.
454-
Use `transpose=False` to make use of the standard convention (this
455-
will be updated in a future release).
456-
457471
Examples
458472
--------
459473
>>> T = np.linspace(0, 10, 100)
460474
>>> U = np.ones((1, 100))
461-
>>> T, Y = ct.forced_response(ct.tf([1], [1, 0.5], True), T, U)
462-
>>> H = ct.markov(Y, U, 3, transpose=False)
475+
>>> response = ct.forced_response(ct.tf([1], [1, 0.5], True), T, U)
476+
>>> H = ct.markov(response, 3)
463477
464478
"""
465479
# Convert input parameters to 2D arrays (if they aren't already)
466-
Umat = np.array(U, ndmin=2)
467-
Ymat = np.array(Y, ndmin=2)
480+
Umat = np.array(data.inputs, ndmin=2)
481+
Ymat = np.array(data.outputs, ndmin=2)
468482

469483
# If data is in transposed format, switch it around
470-
if transpose:
484+
if data.transpose and not data.issiso:
471485
Umat, Ymat = np.transpose(Umat), np.transpose(Ymat)
472486

473-
# Make sure the system is a SISO system
474-
if Umat.shape[0] != 1 or Ymat.shape[0] != 1:
475-
raise ControlMIMONotImplemented
476-
477487
# Make sure the number of time points match
478488
if Umat.shape[1] != Ymat.shape[1]:
479489
raise ControlDimension(
480490
"Input and output data are of differnent lengths")
481-
n = Umat.shape[1]
491+
l = Umat.shape[1]
482492

483493
# If number of desired parameters was not given, set to size of input data
484494
if m is None:
485-
m = Umat.shape[1]
495+
m = l
496+
497+
t = 0
498+
if truncate:
499+
t = m
500+
501+
q = Ymat.shape[0] # number of outputs
502+
p = Umat.shape[0] # number of inputs
486503

487504
# Make sure there is enough data to compute parameters
488-
if m > n:
505+
if m*p > (l-t):
489506
warnings.warn("Not enough data for requested number of parameters")
490507

508+
# the algorithm - Construct a matrix of control inputs to invert
491509
#
492-
# Original algorithm (with mapping to standard order)
493-
#
494-
# RMM note, 24 Dec 2020: This algorithm sets the problem up correctly
495-
# until the final column of the UU matrix is created, at which point it
496-
# makes some modifications that I don't understand. This version of the
497-
# algorithm does not seem to return the actual Markov parameters for a
498-
# system.
499-
#
500-
# # Create the matrix of (shifted) inputs
501-
# UU = np.transpose(Umat)
502-
# for i in range(1, m-1):
503-
# # Shift previous column down and add a zero at the top
504-
# newCol = np.vstack((0, np.reshape(UU[0:n-1, i-1], (-1, 1))))
505-
# UU = np.hstack((UU, newCol))
506-
#
507-
# # Shift previous column down and add a zero at the top
508-
# Ulast = np.vstack((0, np.reshape(UU[0:n-1, m-2], (-1, 1))))
509-
#
510-
# # Replace the elements of the last column new values (?)
511-
# # Each row gets the sum of the rows above it (?)
512-
# for i in range(n-1, 0, -1):
513-
# Ulast[i] = np.sum(Ulast[0:i-1])
514-
# UU = np.hstack((UU, Ulast))
515-
#
516-
# # Solve for the Markov parameters from Y = H @ UU
517-
# # H = [[D], [CB], [CAB], ..., [C A^{m-3} B], [???]]
518-
# H = np.linalg.lstsq(UU, np.transpose(Ymat))[0]
519-
#
520-
# # Markov parameters are in rows => transpose if needed
521-
# return H if transpose else np.transpose(H)
522-
523-
#
524-
# New algorithm - Construct a matrix of control inputs to invert
510+
# (q,l) = (q,p*m) @ (p*m,l)
511+
# YY.T = H @ UU.T
525512
#
526513
# This algorithm sets up the following problem and solves it for
527514
# the Markov parameters
528515
#
516+
# (l,q) = (l,p*m) @ (p*m,q)
517+
# YY = UU @ H.T
518+
#
529519
# [ y(0) ] [ u(0) 0 0 ] [ D ]
530520
# [ y(1) ] [ u(1) u(0) 0 ] [ C B ]
531521
# [ y(2) ] = [ u(2) u(1) u(0) ] [ C A B ]
532522
# [ : ] [ : : : : ] [ : ]
533-
# [ y(n-1) ] [ u(n-1) u(n-2) u(n-3) ... u(n-m) ] [ C A^{m-2} B ]
523+
# [ y(l-1) ] [ u(l-1) u(l-2) u(l-3) ... u(l-m) ] [ C A^{m-2} B ]
534524
#
535-
# Note: if the number of Markov parameters (m) is less than the size of
536-
# the input/output data (n), then this algorithm assumes C A^{j} B = 0
525+
# truncated version t=m, do not use first m equation
526+
#
527+
# [ y(t) ] [ u(t) u(t-1) u(t-2) u(t-m) ] [ D ]
528+
# [ y(t+1) ] [ u(t+1) u(t) u(t-1) u(t-m+1)] [ C B ]
529+
# [ y(t+2) ] = [ u(t+2) u(t+1) u(t) u(t-m+2)] [ C B ]
530+
# [ : ] [ : : : : ] [ : ]
531+
# [ y(l-1) ] [ u(l-1) u(l-2) u(l-3) ... u(l-m) ] [ C A^{m-2} B ]
532+
#
533+
# Note: This algorithm assumes C A^{j} B = 0
537534
# for j > m-2. See equation (3) in
538535
#
539536
# J.-N. Juang, M. Phan, L. G. Horta, and R. W. Longman, Identification
@@ -542,17 +539,40 @@ def markov(Y, U, m=None, transpose=False):
542539
# 320-329, 2012. http://doi.org/10.2514/3.21006
543540
#
544541

542+
# Set up the full problem
545543
# Create matrix of (shifted) inputs
546-
UU = Umat
547-
for i in range(1, m):
548-
# Shift previous column down and add a zero at the top
549-
new_row = np.hstack((0, UU[i-1, 0:-1]))
550-
UU = np.vstack((UU, new_row))
551-
UU = np.transpose(UU)
552-
553-
# Invert and solve for Markov parameters
554-
YY = np.transpose(Ymat)
555-
H, _, _, _ = np.linalg.lstsq(UU, YY, rcond=None)
556-
544+
UUT = np.zeros((p*m,(l)))
545+
for i in range(m):
546+
# Shift previous column down and keep zeros at the top
547+
UUT[i*p:(i+1)*p,i:] = Umat[:,:l-i]
548+
549+
# Truncate first t=0 or t=m time steps, transpose the problem for lsq
550+
YY = Ymat[:,t:].T
551+
UU = UUT[:,t:].T
552+
553+
# Solve for the Markov parameters from YY = UU @ H.T
554+
HT, _, _, _ = np.linalg.lstsq(UU, YY, rcond=None)
555+
H = HT.T/dt # scaling
556+
557+
H = H.reshape(q,m,p) # output, time*input -> output, time, input
558+
H = H.transpose(0,2,1) # output, input, time
559+
560+
# Create unit area impulse inputs
561+
inputs = np.zeros((q,p,m))
562+
trace_labels, trace_types = [], []
563+
for i in range(p):
564+
inputs[i,i,0] = 1/dt # unit area impulse
565+
trace_labels.append(f"From {data.input_labels[i]}")
566+
trace_types.append('impulse')
567+
568+
# Markov parameters as TimeResponseData with unit area impulse inputs
557569
# Return the first m Markov parameters
558-
return H if transpose else np.transpose(H)
570+
return TimeResponseData(time=data.time[:m],
571+
outputs=H,
572+
output_labels=data.output_labels,
573+
inputs=inputs,
574+
input_labels=data.input_labels,
575+
trace_labels=trace_labels,
576+
trace_types=trace_types,
577+
transpose=data.transpose,
578+
issiso=data.issiso)

‎control/tests/modelsimp_test.py

Lines changed: 30 additions & 20 deletions
Original file line numberDiff line numberDiff line change
@@ -7,7 +7,7 @@
77
import pytest
88

99

10-
from control import StateSpace, forced_response, tf, rss, c2d
10+
from control import StateSpace, forced_response, tf, rss, c2d, TimeResponseData
1111
from control.exception import ControlMIMONotImplemented
1212
from control.tests.conftest import slycotonly
1313
from control.modelsimp import balred, hsvd, markov, modred
@@ -33,36 +33,44 @@ def testHSVD(self):
3333
assert not isinstance(hsv, np.matrix)
3434

3535
def testMarkovSignature(self):
36-
U = np.array([[1., 1., 1., 1., 1.]])
36+
U = np.array([1., 1., 1., 1., 1.])
3737
Y = U
38+
response = TimeResponseData(time=np.arange(U.shape[-1]),
39+
outputs=Y,
40+
output_labels='y',
41+
inputs=U,
42+
input_labels='u',
43+
)
3844
m = 3
39-
H = markov(Y, U, m, transpose=False)
40-
Htrue = np.array([[1., 0., 0.]])
41-
np.testing.assert_array_almost_equal(H, Htrue)
45+
H = markov(response, m)
46+
Htrue = np.array([1., 0., 0.])
47+
np.testing.assert_array_almost_equal(H.outputs, Htrue)
4248

4349
# Make sure that transposed data also works
44-
H = markov(np.transpose(Y), np.transpose(U), m, transpose=True)
45-
np.testing.assert_array_almost_equal(H, np.transpose(Htrue))
50+
response.transpose=True
51+
H = markov(response, m)
52+
np.testing.assert_array_almost_equal(H.outputs, np.transpose(Htrue))
4653

4754
# Generate Markov parameters without any arguments
48-
H = markov(Y, U, m)
49-
np.testing.assert_array_almost_equal(H, Htrue)
55+
response.transpose=False
56+
H = markov(response, m)
57+
np.testing.assert_array_almost_equal(H.outputs, Htrue)
5058

5159
# Test example from docstring
5260
T = np.linspace(0, 10, 100)
5361
U = np.ones((1, 100))
54-
T, Y = forced_response(tf([1], [1, 0.5], True), T, U)
55-
H = markov(Y, U, 3, transpose=False)
62+
response = forced_response(tf([1], [1, 0.5], True), T, U)
63+
H = markov(response, 3)
5664

5765
# Test example from issue #395
58-
inp = np.array([1, 2])
59-
outp = np.array([2, 4])
60-
mrk = markov(outp, inp, 1, transpose=False)
66+
#inp = np.array([1, 2])
67+
#outp = np.array([2, 4])
68+
#mrk = markov(outp, inp, 1, transpose=False)
6169

6270
# Make sure MIMO generates an error
63-
U = np.ones((2, 100)) # 2 inputs (Y unchanged, with 1 output)
64-
with pytest.raises(ControlMIMONotImplemented):
65-
markov(Y, U, m)
71+
#U = np.ones((2, 100)) # 2 inputs (Y unchanged, with 1 output)
72+
#with pytest.raises(ControlMIMONotImplemented):
73+
# markov(Y, U, m)
6674

6775
# Make sure markov() returns the right answer
6876
@pytest.mark.parametrize("k, m, n",
@@ -98,18 +106,20 @@ def testMarkovResults(self, k, m, n):
98106
Mtrue = np.hstack([Hd.D] + [
99107
Hd.C @ np.linalg.matrix_power(Hd.A, i) @ Hd.B
100108
for i in range(m-1)])
109+
110+
Mtrue = np.squeeze(Mtrue)
101111

102112
# Generate input/output data
103113
T = np.array(range(n)) * Ts
104114
U = np.cos(T) + np.sin(T/np.pi)
105-
_, Y = forced_response(Hd, T, U, squeeze=True)
106-
Mcomp = markov(Y, U, m)
115+
response = forced_response(Hd, T, U, squeeze=True)
116+
Mcomp = markov(response, m)
107117

108118
# Compare to results from markov()
109119
# experimentally determined probability to get non matching results
110120
# with rtot=1e-6 and atol=1e-8 due to numerical errors
111121
# for k=5, m=n=10: 0.015 %
112-
np.testing.assert_allclose(Mtrue, Mcomp, rtol=1e-6, atol=1e-8)
122+
np.testing.assert_allclose(Mtrue, Mcomp.outputs, rtol=1e-6, atol=1e-8)
113123

114124
def testModredMatchDC(self):
115125
#balanced realization computed in matlab for the transfer function:

0 commit comments

Comments
 (0)
0