• Home
  • Features
  • Pricing
  • Docs
  • Announcements
  • Sign In

python-control / python-control / 12266904606

11 Dec 2024 12:03AM UTC coverage: 94.73% (+0.009%) from 94.721%
12266904606

Pull #1078

github

web-flow
Merge 0b9be7a84 into e2c0ff85f
Pull Request #1078: fix issue with multiplying MIMO LTI system by scalar

9330 of 9849 relevant lines covered (94.73%)

8.27 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

92.73
control/modelsimp.py
1
# modelsimp.py - tools for model simplification
2
#
3
# Author: Steve Brunton, Kevin Chen, Lauren Padilla
4
# Date: 30 Nov 2010
5
#
6
# This file contains routines for obtaining reduced order models
7
#
8

9
import warnings
9✔
10

11
# External packages and modules
12
import numpy as np
9✔
13

14
from .exception import ControlArgument, ControlDimension, ControlSlycot
9✔
15
from .iosys import isctime, isdtime
9✔
16
from .statefbk import gram
9✔
17
from .statesp import StateSpace
9✔
18
from .timeresp import TimeResponseData
9✔
19

20
__all__ = ['hankel_singular_values', 'balanced_reduction', 'model_reduction',
9✔
21
           'minimal_realization', 'eigensys_realization', 'markov', 'hsvd',
22
           'balred', 'modred', 'minreal', 'era']
23

24

25
# Hankel Singular Value Decomposition
26
#
27
# The following returns the Hankel singular values, which are singular values
28
# of the matrix formed by multiplying the controllability and observability
29
# Gramians
30
def hankel_singular_values(sys):
9✔
31
    """Calculate the Hankel singular values.
32

33
    Parameters
34
    ----------
35
    sys : StateSpace
36
        A state space system
37

38
    Returns
39
    -------
40
    H : array
41
        A list of Hankel singular values
42

43
    See Also
44
    --------
45
    gram
46

47
    Notes
48
    -----
49
    The Hankel singular values are the singular values of the Hankel operator.
50
    In practice, we compute the square root of the eigenvalues of the matrix
51
    formed by taking the product of the observability and controllability
52
    gramians.  There are other (more efficient) methods based on solving the
53
    Lyapunov equation in a particular way (more details soon).
54

55
    Examples
56
    --------
57
    >>> G = ct.tf2ss([1], [1, 2])
58
    >>> H = ct.hsvd(G)
59
    >>> H[0]
60
    np.float64(0.25)
61

62
    """
63
    # TODO: implement for discrete time systems
64
    if (isdtime(sys, strict=True)):
5✔
65
        raise NotImplementedError("Function not implemented in discrete time")
66

67
    Wc = gram(sys, 'c')
5✔
68
    Wo = gram(sys, 'o')
5✔
69
    WoWc = Wo @ Wc
5✔
70
    w, v = np.linalg.eig(WoWc)
5✔
71

72
    hsv = np.sqrt(w)
5✔
73
    hsv = np.array(hsv)
5✔
74
    hsv = np.sort(hsv)
5✔
75
    # Return the Hankel singular values, high to low
76
    return hsv[::-1]
5✔
77

78

79
def model_reduction(
9✔
80
        sys, elim_states=None, method='matchdc', elim_inputs=None,
81
        elim_outputs=None, keep_states=None, keep_inputs=None,
82
        keep_outputs=None, warn_unstable=True):
83
    """Model reduction by input, output, or state elimination.
84

85
    This function produces a reduced-order model of a system by eliminating
86
    specified inputs, outputs, and/or states from the original system.  The
87
    specific states, inputs, or outputs that are eliminated can be
88
    specified by either listing the states, inputs, or outputs to be
89
    eliminated or those to be kept.
90

91
    Two methods of state reduction are possible: 'truncate' removes the
92
    states marked for elimination, while 'matchdc' replaces the eliminated
93
    states with their equilibrium values (thereby keeping the input/output
94
    gain unchanged at zero frequency ["DC"]).
95

96
    Parameters
97
    ----------
98
    sys : StateSpace
99
        Original system to reduce.
100
    elim_inputs, elim_outputs, elim_states : array of int or str, optional
101
        Vector of inputs, outputs, or states to eliminate.  Can be specified
102
        either as an offset into the appropriate vector or as a signal name.
103
    keep_inputs, keep_outputs, keep_states : array, optional
104
        Vector of inputs, outputs, or states to keep.  Can be specified
105
        either as an offset into the appropriate vector or as a signal name.
106
    method : string
107
        Method of removing states: either 'truncate' or 'matchdc' (default).
108
    warn_unstable : bool, option
109
        If `False`, don't warn if system is unstable.
110

111
    Returns
112
    -------
113
    rsys : StateSpace
114
        Reduced order model.
115

116
    Raises
117
    ------
118
    ValueError
119
        If `method` is not either 'matchdc' or 'truncate'.
120
    NotImplementedError
121
        If the 'matchdc' method is used for a discrete time system.
122

123
    Warns
124
    -----
125
    UserWarning
126
        If eigenvalues of `sys.A` are not all stable.
127

128
    Examples
129
    --------
130
    >>> G = ct.rss(4)
131
    >>> Gr = ct.model_reduction(G, [0, 2], method='matchdc')
132
    >>> Gr.nstates
133
    2
134

135
    See Also
136
    --------
137
    balanced_reduction : Eliminate states using Hankel singular values.
138
    minimal_realization : Eliminate unreachable or unobservable states.
139

140
    Notes
141
    -----
142
    The model_reduction function issues a warning if the system has
143
    unstable eigenvalues, since in those situations the stability of the
144
    reduced order model may be different than the stability of the full
145
    model.  No other checking is done, so users must to be careful not to
146
    render a system unobservable or unreachable.
147

148
    States, inputs, and outputs can be specified using integer offsets or
149
    using signal names.  Slices can also be specified, but must use the
150
    Python ``slice()`` function.
151

152
    """
153
    if not isinstance(sys, StateSpace):
9✔
154
        raise TypeError("system must be a StateSpace system")
×
155

156
    # Check system is stable
157
    if warn_unstable:
9✔
158
        if isctime(sys) and np.any(np.linalg.eigvals(sys.A).real >= 0.0) or \
9✔
159
           isdtime(sys) and np.any(np.abs(np.linalg.eigvals(sys.A)) >= 1):
160
            warnings.warn("System is unstable; reduction may be meaningless")
9✔
161

162
    # Utility function to process keep/elim keywords
163
    def _process_elim_or_keep(elim, keep, labels):
9✔
164
        def _expand_key(key):
9✔
165
            if key is None:
9✔
166
                return []
9✔
167
            elif isinstance(key, str):
9✔
168
                return labels.index(key)
9✔
169
            elif isinstance(key, list):
9✔
170
                return [_expand_key(k) for k in key]
9✔
171
            elif isinstance(key, slice):
9✔
172
                return range(len(labels))[key]
9✔
173
            else:
174
                return key
9✔
175

176
        elim = np.atleast_1d(_expand_key(elim))
9✔
177
        keep = np.atleast_1d(_expand_key(keep))
9✔
178
            
179
        if len(elim) > 0 and len(keep) > 0:
9✔
180
            raise ValueError(
9✔
181
                "can't provide both 'keep' and 'elim' for same variables")
182
        elif len(keep) > 0:
9✔
183
            keep = np.sort(keep).tolist()
9✔
184
            elim = [i for i in range(len(labels)) if i not in keep]
9✔
185
        else:
186
            elim = [] if elim is None else np.sort(elim).tolist()
9✔
187
            keep = [i for i in range(len(labels)) if i not in elim]
9✔
188
        return elim, keep
9✔
189

190
    # Determine which states to keep
191
    elim_states, keep_states = _process_elim_or_keep(
9✔
192
        elim_states, keep_states, sys.state_labels)
193
    elim_inputs, keep_inputs = _process_elim_or_keep(
9✔
194
        elim_inputs, keep_inputs, sys.input_labels)
195
    elim_outputs, keep_outputs = _process_elim_or_keep(
9✔
196
        elim_outputs, keep_outputs, sys.output_labels)
197

198
    # Create submatrix of states we are keeping
199
    A11 = sys.A[:, keep_states][keep_states, :]     # states we are keeping
9✔
200
    A12 = sys.A[:, elim_states][keep_states, :]     # needed for 'matchdc'
9✔
201
    A21 = sys.A[:, keep_states][elim_states, :]
9✔
202
    A22 = sys.A[:, elim_states][elim_states, :]
9✔
203

204
    B1 = sys.B[keep_states, :]
9✔
205
    B2 = sys.B[elim_states, :]
9✔
206

207
    C1 = sys.C[:, keep_states]
9✔
208
    C2 = sys.C[:, elim_states]
9✔
209

210
    # Figure out the new state space system
211
    if method == 'matchdc' and A22.size > 0:
9✔
212
        if sys.isdtime(strict=True):
9✔
213
            raise NotImplementedError(
214
                "'matchdc' not (yet) supported for discrete time systems")
215
        
216
        # if matchdc, residualize
217
        # Check if the matrix A22 is invertible
218
        if np.linalg.matrix_rank(A22) != len(elim_states):
9✔
219
            raise ValueError("Matrix A22 is singular to working precision.")
×
220

221
        # Now precompute A22\A21 and A22\B2 (A22I = inv(A22))
222
        # We can solve two linear systems in one pass, since the
223
        # coefficients matrix A22 is the same. Thus, we perform the LU
224
        # decomposition (cubic runtime complexity) of A22 only once!
225
        # The remaining back substitutions are only quadratic in runtime.
226
        A22I_A21_B2 = np.linalg.solve(A22, np.concatenate((A21, B2), axis=1))
9✔
227
        A22I_A21 = A22I_A21_B2[:, :A21.shape[1]]
9✔
228
        A22I_B2 = A22I_A21_B2[:, A21.shape[1]:]
9✔
229

230
        Ar = A11 - A12 @ A22I_A21
9✔
231
        Br = B1 - A12 @ A22I_B2
9✔
232
        Cr = C1 - C2 @ A22I_A21
9✔
233
        Dr = sys.D - C2 @ A22I_B2
9✔
234

235
    elif method == 'truncate' or A22.size == 0:
9✔
236
        # Get rid of unwanted states
237
        Ar = A11
9✔
238
        Br = B1
9✔
239
        Cr = C1
9✔
240
        Dr = sys.D
9✔
241

242
    else:
243
        raise ValueError("Oops, method is not supported!")
×
244

245
    # Get rid of additional inputs and outputs
246
    Br = Br[:, keep_inputs]
9✔
247
    Cr = Cr[keep_outputs, :]
9✔
248
    Dr = Dr[keep_outputs, :][:, keep_inputs]
9✔
249

250
    rsys = StateSpace(Ar, Br, Cr, Dr)
9✔
251
    return rsys
9✔
252

253

254
def balanced_reduction(sys, orders, method='truncate', alpha=None):
9✔
255
    """Balanced reduced order model of sys of a given order.
256

257
    States are eliminated based on Hankel singular value.
258
    If sys has unstable modes, they are removed, the
259
    balanced realization is done on the stable part, then
260
    reinserted in accordance with the reference below.
261

262
    Reference: Hsu,C.S., and Hou,D., 1991,
263
    Reducing unstable linear control systems via real Schur transformation.
264
    Electronics Letters, 27, 984-986.
265

266
    Parameters
267
    ----------
268
    sys : StateSpace
269
        Original system to reduce.
270
    orders : integer or array of integer
271
        Desired order of reduced order model (if a vector, returns a vector
272
        of systems).
273
    method : string
274
        Method of removing states, either ``'truncate'`` or ``'matchdc'``..
275
    alpha : float
276
        Redefines the stability boundary for eigenvalues of the system
277
        matrix A.  By default for continuous-time systems, alpha <= 0
278
        defines the stability boundary for the real part of A's eigenvalues
279
        and for discrete-time systems, 0 <= alpha <= 1 defines the stability
280
        boundary for the modulus of A's eigenvalues. See SLICOT routines
281
        AB09MD and AB09ND for more information.
282

283
    Returns
284
    -------
285
    rsys : StateSpace
286
        A reduced order model or a list of reduced order models if orders is
287
        a list.
288

289
    Raises
290
    ------
291
    ValueError
292
        If `method` is not ``'truncate'`` or ``'matchdc'``
293
    ImportError
294
        if slycot routine ab09ad, ab09md, or ab09nd is not found
295

296
    ValueError
297
        if there are more unstable modes than any value in orders
298

299
    Examples
300
    --------
301
    >>> G = ct.rss(4)
302
    >>> Gr = ct.balred(G, orders=2, method='matchdc')
303
    >>> Gr.nstates
304
    2
305

306
    """
307
    if method != 'truncate' and method != 'matchdc':
5✔
308
        raise ValueError("supported methods are 'truncate' or 'matchdc'")
×
309
    elif method == 'truncate':
5✔
310
        try:
5✔
311
            from slycot import ab09ad, ab09md
5✔
312
        except ImportError:
×
313
            raise ControlSlycot(
×
314
                "can't find slycot subroutine ab09md or ab09ad")
315
    elif method == 'matchdc':
5✔
316
        try:
5✔
317
            from slycot import ab09nd
5✔
318
        except ImportError:
×
319
            raise ControlSlycot("can't find slycot subroutine ab09nd")
×
320

321
    # Check for ss system object, need a utility for this?
322

323
    # TODO: Check for continous or discrete, only continuous supported for now
324
    #   if isCont():
325
    #       dico = 'C'
326
    #   elif isDisc():
327
    #       dico = 'D'
328
    #   else:
329
    dico = 'C'
5✔
330

331
    job = 'B'                   # balanced (B) or not (N)
5✔
332
    equil = 'N'                 # scale (S) or not (N)
5✔
333
    if alpha is None:
5✔
334
        if dico == 'C':
5✔
335
            alpha = 0.
5✔
336
        elif dico == 'D':
×
337
            alpha = 1.
×
338

339
    rsys = []                   # empty list for reduced systems
5✔
340

341
    # check if orders is a list or a scalar
342
    try:
5✔
343
        order = iter(orders)
5✔
344
    except TypeError:           # if orders is a scalar
5✔
345
        orders = [orders]
5✔
346

347
    for i in orders:
5✔
348
        n = np.size(sys.A, 0)
5✔
349
        m = np.size(sys.B, 1)
5✔
350
        p = np.size(sys.C, 0)
5✔
351
        if method == 'truncate':
5✔
352
            # check system stability
353
            if np.any(np.linalg.eigvals(sys.A).real >= 0.0):
5✔
354
                # unstable branch
355
                Nr, Ar, Br, Cr, Ns, hsv = ab09md(
×
356
                    dico, job, equil, n, m, p, sys.A, sys.B, sys.C,
357
                    alpha=alpha, nr=i, tol=0.0)
358
            else:
359
                # stable branch
360
                Nr, Ar, Br, Cr, hsv = ab09ad(
5✔
361
                    dico, job, equil, n, m, p, sys.A, sys.B, sys.C,
362
                    nr=i, tol=0.0)
363
            rsys.append(StateSpace(Ar, Br, Cr, sys.D))
5✔
364

365
        elif method == 'matchdc':
5✔
366
            Nr, Ar, Br, Cr, Dr, Ns, hsv = ab09nd(
5✔
367
                dico, job, equil, n, m, p, sys.A, sys.B, sys.C, sys.D,
368
                alpha=alpha, nr=i, tol1=0.0, tol2=0.0)
369
            rsys.append(StateSpace(Ar, Br, Cr, Dr))
5✔
370

371
    # if orders was a scalar, just return the single reduced model, not a list
372
    if len(orders) == 1:
5✔
373
        return rsys[0]
5✔
374
    # if orders was a list/vector, return a list/vector of systems
375
    else:
376
        return rsys
5✔
377

378

379
def minimal_realization(sys, tol=None, verbose=True):
9✔
380
    """ Eliminate uncontrollable or unobservable states.
381

382
    Eliminates uncontrollable or unobservable states in state-space
383
    models or cancelling pole-zero pairs in transfer functions. The
384
    output sysr has minimal order and the same response
385
    characteristics as the original model sys.
386

387
    Parameters
388
    ----------
389
    sys : StateSpace or TransferFunction
390
        Original system.
391
    tol : real
392
        Tolerance.
393
    verbose : bool
394
        Print results if True.
395

396
    Returns
397
    -------
398
    rsys : StateSpace or TransferFunction
399
        Cleaned model.
400

401
    """
402
    sysr = sys.minreal(tol)
5✔
403
    if verbose:
5✔
404
        print("{nstates} states have been removed from the model".format(
5✔
405
                nstates=len(sys.poles()) - len(sysr.poles())))
406
    return sysr
5✔
407

408

409
def _block_hankel(Y, m, n):
9✔
410
    """Create a block Hankel matrix from impulse response."""
411
    q, p, _ = Y.shape
9✔
412
    YY = Y.transpose(0, 2, 1) # transpose for reshape
9✔
413

414
    H = np.zeros((q*m, p*n))
9✔
415

416
    for r in range(m):
9✔
417
        # shift and add row to Hankel matrix
418
        new_row = YY[:, r:r+n, :]
9✔
419
        H[q*r:q*(r+1), :] = new_row.reshape((q, p*n))
9✔
420

421
    return H
9✔
422

423

424
def eigensys_realization(arg, r, m=None, n=None, dt=True, transpose=False):
9✔
425
    r"""eigensys_realization(YY, r)
426

427
    Calculate ERA model of order `r` based on impulse-response data `YY`.
428

429
    This function computes a discrete time system
430

431
    .. math::
432

433
        x[k+1] &= A x[k] + B u[k] \\\\
434
        y[k] &= C x[k] + D u[k]
435

436
    for a given impulse-response data (see [1]_).
437

438
    The function can be called with 2 arguments:
439

440
    * ``sysd, S = eigensys_realization(data, r)``
441
    * ``sysd, S = eigensys_realization(YY, r)``
442

443
    where `data` is a `TimeResponseData` object, `YY` is a 1D or 3D
444
    array, and r is an integer.
445

446
    Parameters
447
    ----------
448
    YY : array_like
449
        Impulse response from which the StateSpace model is estimated, 1D
450
        or 3D array.
451
    data : TimeResponseData
452
        Impulse response from which the StateSpace model is estimated.
453
    r : integer
454
        Order of model.
455
    m : integer, optional
456
        Number of rows in Hankel matrix. Default is 2*r.
457
    n : integer, optional
458
        Number of columns in Hankel matrix. Default is 2*r.
459
    dt : True or float, optional
460
        True indicates discrete time with unspecified sampling time and a
461
        positive float is discrete time with the specified sampling time.
462
        It can be used to scale the StateSpace model in order to match the
463
        unit-area impulse response of python-control. Default is True.
464
    transpose : bool, optional
465
        Assume that input data is transposed relative to the standard
466
        :ref:`time-series-convention`. For TimeResponseData this parameter
467
        is ignored. Default is False.
468

469
    Returns
470
    -------
471
    sys : StateSpace
472
        A reduced order model sys=StateSpace(Ar,Br,Cr,Dr,dt).
473
    S : array
474
        Singular values of Hankel matrix. Can be used to choose a good r
475
        value.
476

477
    References
478
    ----------
479
    .. [1] Samet Oymak and Necmiye Ozay, Non-asymptotic Identification of
480
       LTI Systems from a Single Trajectory.
481
       https://arxiv.org/abs/1806.05722
482

483
    Examples
484
    --------
485
    >>> T = np.linspace(0, 10, 100)
486
    >>> _, YY = ct.impulse_response(ct.tf([1], [1, 0.5], True), T)
487
    >>> sysd, _ = ct.eigensys_realization(YY, r=1)
488

489
    >>> T = np.linspace(0, 10, 100)
490
    >>> response = ct.impulse_response(ct.tf([1], [1, 0.5], True), T)
491
    >>> sysd, _ = ct.eigensys_realization(response, r=1)
492
    """
493
    if isinstance(arg, TimeResponseData):
9✔
494
        YY = np.array(arg.outputs, ndmin=3)
9✔
495
        if arg.transpose:
9✔
496
            YY = np.transpose(YY)
×
497
    else:
498
        YY = np.array(arg, ndmin=3)
9✔
499
        if transpose:
9✔
500
            YY = np.transpose(YY)
×
501

502
    q, p, l = YY.shape
9✔
503

504
    if m is None:
9✔
505
        m = 2*r
9✔
506
    if n is None:
9✔
507
        n = 2*r
9✔
508

509
    if m*q < r or n*p < r:
9✔
510
        raise ValueError("Hankel parameters are to small")
×
511

512
    if (l-1) < m+n:
9✔
513
        raise ValueError("not enough data for requested number of parameters")
×
514

515
    H = _block_hankel(YY[:, :, 1:], m, n+1) # Hankel matrix (q*m, p*(n+1))
9✔
516
    Hf = H[:, :-p] # first p*n columns of H
9✔
517
    Hl = H[:, p:] # last p*n columns of H
9✔
518

519
    U,S,Vh = np.linalg.svd(Hf, True)
9✔
520
    Ur =U[:, 0:r]
9✔
521
    Vhr =Vh[0:r, :]
9✔
522

523
    # balanced realizations
524
    Sigma_inv = np.diag(1./np.sqrt(S[0:r]))
9✔
525
    Ar = Sigma_inv @ Ur.T @ Hl @ Vhr.T @ Sigma_inv
9✔
526
    Br = Sigma_inv @ Ur.T @ Hf[:, 0:p]*dt # dt scaling for unit-area impulse
9✔
527
    Cr = Hf[0:q, :] @ Vhr.T @ Sigma_inv
9✔
528
    Dr = YY[:, :, 0]
9✔
529

530
    return StateSpace(Ar, Br, Cr, Dr, dt), S
9✔
531

532

533
def markov(*args, m=None, transpose=False, dt=None, truncate=False):
9✔
534
    """markov(Y, U, [, m])
535

536
    Calculate the first `m` Markov parameters [D CB CAB ...] from data.
537

538
    This function computes the Markov parameters for a discrete time
539
    system
540

541
    .. math::
542

543
        x[k+1] &= A x[k] + B u[k] \\\\
544
        y[k] &= C x[k] + D u[k]
545

546
    given data for u and y.  The algorithm assumes that that C A^k B = 0
547
    for k > m-2 (see [1]_).  Note that the problem is ill-posed if the
548
    length of the input data is less than the desired number of Markov
549
    parameters (a warning message is generated in this case).
550

551
    The function can be called with either 1, 2 or 3 arguments:
552

553
    * ``H = markov(data)``
554
    * ``H = markov(data, m)``
555
    * ``H = markov(Y, U)``
556
    * ``H = markov(Y, U, m)``
557

558
    where `data` is a `TimeResponseData` object, `YY` is a 1D or 3D
559
    array, and r is an integer.
560

561
    Parameters
562
    ----------
563
    Y : array_like
564
        Output data. If the array is 1D, the system is assumed to be
565
        single input. If the array is 2D and transpose=False, the columns
566
        of `Y` are taken as time points, otherwise the rows of `Y` are
567
        taken as time points.
568
    U : array_like
569
        Input data, arranged in the same way as `Y`.
570
    data : TimeResponseData
571
        Response data from which the Markov parameters where estimated.
572
        Input and output data must be 1D or 2D array.
573
    m : int, optional
574
        Number of Markov parameters to output. Defaults to len(U).
575
    dt : True of float, optional
576
        True indicates discrete time with unspecified sampling time and a
577
        positive float is discrete time with the specified sampling time.
578
        It can be used to scale the Markov parameters in order to match
579
        the unit-area impulse response of python-control. Default is True
580
        for array_like and dt=data.time[1]-data.time[0] for
581
        TimeResponseData as input.
582
    truncate : bool, optional
583
        Do not use first m equation for least squares. Default is False.
584
    transpose : bool, optional
585
        Assume that input data is transposed relative to the standard
586
        :ref:`time-series-convention`. For TimeResponseData this parameter
587
        is ignored. Default is False.
588

589
    Returns
590
    -------
591
    H : ndarray
592
        First m Markov parameters, [D CB CAB ...].
593

594
    References
595
    ----------
596
    .. [1] J.-N. Juang, M. Phan, L. G.  Horta, and R. W. Longman,
597
       Identification of observer/Kalman filter Markov parameters - Theory
598
       and experiments. Journal of Guidance Control and Dynamics, 16(2),
599
       320-329, 2012. http://doi.org/10.2514/3.21006
600

601
    Examples
602
    --------
603
    >>> T = np.linspace(0, 10, 100)
604
    >>> U = np.ones((1, 100))
605
    >>> T, Y = ct.forced_response(ct.tf([1], [1, 0.5], True), T, U)
606
    >>> H = ct.markov(Y, U, 3, transpose=False)
607

608
    """
609

610
    # Convert input parameters to 2D arrays (if they aren't already)
611
    # Get the system description
612
    if len(args) < 1:
9✔
613
        raise ControlArgument("not enough input arguments")
9✔
614

615
    if isinstance(args[0], TimeResponseData):
9✔
616
        data = args[0]
9✔
617
        Umat = np.array(data.inputs, ndmin=2)
9✔
618
        Ymat = np.array(data.outputs, ndmin=2)
9✔
619
        if dt is None:
9✔
620
            dt = data.time[1] - data.time[0]
9✔
621
            if not np.allclose(np.diff(data.time), dt):
9✔
622
                raise ValueError("response time values must be equally "
×
623
                                 "spaced.")
624
        transpose = data.transpose
9✔
625
        if data.transpose and not data.issiso:
9✔
626
            Umat, Ymat = np.transpose(Umat), np.transpose(Ymat)
9✔
627
        if len(args) == 2:
9✔
628
            m = args[1]
9✔
629
        elif len(args) > 2:
9✔
630
            raise ControlArgument("too many positional arguments")
9✔
631
    else:
632
        if len(args) < 2:
9✔
633
            raise ControlArgument("not enough input arguments")
9✔
634
        Umat = np.array(args[1], ndmin=2)
9✔
635
        Ymat = np.array(args[0], ndmin=2)
9✔
636
        if dt is None:
9✔
637
            dt = True
9✔
638
        if transpose:
9✔
639
            Umat, Ymat = np.transpose(Umat), np.transpose(Ymat)
9✔
640
        if len(args) == 3:
9✔
641
            m = args[2]
9✔
642
        elif len(args) > 3:
9✔
643
            raise ControlArgument("too many positional arguments")
9✔
644

645
    # Make sure the number of time points match
646
    if Umat.shape[1] != Ymat.shape[1]:
9✔
647
        raise ControlDimension(
9✔
648
            "Input and output data are of differnent lengths")
649
    l = Umat.shape[1]
9✔
650

651
    # If number of desired parameters was not given, set to size of input data
652
    if m is None:
9✔
653
        m = l
9✔
654

655
    t = 0
9✔
656
    if truncate:
9✔
657
        t = m
9✔
658

659
    q = Ymat.shape[0] # number of outputs
9✔
660
    p = Umat.shape[0] # number of inputs
9✔
661

662
    # Make sure there is enough data to compute parameters
663
    if m*p > (l-t):
9✔
664
        warnings.warn("Not enough data for requested number of parameters")
9✔
665

666
    # the algorithm - Construct a matrix of control inputs to invert
667
    #
668
    # (q,l)   = (q,p*m) @ (p*m,l)
669
    # YY.T    = H @ UU.T
670
    #
671
    # This algorithm sets up the following problem and solves it for
672
    # the Markov parameters
673
    #
674
    # (l,q)   = (l,p*m) @ (p*m,q)
675
    # YY      = UU @ H.T
676
    #
677
    # [ y(0)   ]   [ u(0)    0       0                 ] [ D           ]
678
    # [ y(1)   ]   [ u(1)    u(0)    0                 ] [ C B         ]
679
    # [ y(2)   ] = [ u(2)    u(1)    u(0)              ] [ C A B       ]
680
    # [  :     ]   [  :      :        :          :     ] [  :          ]
681
    # [ y(l-1) ]   [ u(l-1)  u(l-2)  u(l-3) ... u(l-m) ] [ C A^{m-2} B ]
682
    #
683
    # truncated version t=m, do not use first m equation
684
    #
685
    # [ y(t)   ]   [ u(t)    u(t-1)  u(t-2)    u(t-m)  ] [ D           ]
686
    # [ y(t+1) ]   [ u(t+1)  u(t)    u(t-1)    u(t-m+1)] [ C B         ]
687
    # [ y(t+2) ] = [ u(t+2)  u(t+1)  u(t)      u(t-m+2)] [ C B         ]
688
    # [  :     ]   [  :      :        :          :     ] [  :          ]
689
    # [ y(l-1) ]   [ u(l-1)  u(l-2)  u(l-3) ... u(l-m) ] [ C A^{m-2} B ]
690
    #
691
    # Note: This algorithm assumes C A^{j} B = 0
692
    # for j > m-2.  See equation (3) in
693
    #
694
    #   J.-N. Juang, M. Phan, L. G. Horta, and R. W. Longman, Identification
695
    #   of observer/Kalman filter Markov parameters - Theory and
696
    #   experiments. Journal of Guidance Control and Dynamics, 16(2),
697
    #   320-329, 2012. http://doi.org/10.2514/3.21006
698
    #
699

700
    # Set up the full problem
701
    # Create matrix of (shifted) inputs
702
    UUT = np.zeros((p*m, l))
9✔
703
    for i in range(m):
9✔
704
        # Shift previous column down and keep zeros at the top
705
        UUT[i*p:(i+1)*p, i:] = Umat[:, :l-i]
9✔
706

707
    # Truncate first t=0 or t=m time steps, transpose the problem for lsq
708
    YY = Ymat[:, t:].T
9✔
709
    UU = UUT[:, t:].T
9✔
710

711
    # Solve for the Markov parameters from  YY = UU @ H.T
712
    HT, _, _, _ = np.linalg.lstsq(UU, YY, rcond=None)
9✔
713
    H = HT.T/dt # scaling
9✔
714

715
    H = H.reshape(q, m, p) # output, time*input -> output, time, input
9✔
716
    H = H.transpose(0, 2, 1) # output, input, time
9✔
717

718
    # for siso return a 1D array instead of a 3D array
719
    if q == 1 and p == 1:
9✔
720
        H = np.squeeze(H)
9✔
721

722
    # Return the first m Markov parameters
723
    return H if not transpose else np.transpose(H)
9✔
724

725
# Function aliases
726
hsvd = hankel_singular_values
9✔
727
balred = balanced_reduction
9✔
728
modred = model_reduction
9✔
729
minreal = minimal_realization
9✔
730
era = eigensys_realization
9✔
STATUS · Troubleshooting · Open an Issue · Sales · Support · CAREERS · ENTERPRISE · START FREE · SCHEDULE DEMO
ANNOUNCEMENTS · TWITTER · TOS & SLA · Supported CI Services · What's a CI service? · Automated Testing

© 2025 Coveralls, Inc