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

python-control / python-control / 13107996232

03 Feb 2025 06:53AM UTC coverage: 94.731% (+0.02%) from 94.709%
13107996232

push

github

web-flow
Merge pull request #1094 from murrayrm/userguide-22Dec2024

Updated user documentation (User Guide, Reference Manual)

9673 of 10211 relevant lines covered (94.73%)

8.28 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
# Initial authors: Steve Brunton, Kevin Chen, Lauren Padilla
4
# Creation date: 30 Nov 2010
5

6
"""Tools for model simplification.
7

8
This module contains routines for obtaining reduced order models for state
9
space systems.
10

11
"""
12

13
import warnings
9✔
14

15
# External packages and modules
16
import numpy as np
9✔
17

18
from .exception import ControlArgument, ControlDimension, ControlSlycot
9✔
19
from .iosys import isctime, isdtime
9✔
20
from .statefbk import gram
9✔
21
from .statesp import StateSpace
9✔
22
from .timeresp import TimeResponseData
9✔
23

24
__all__ = ['hankel_singular_values', 'balanced_reduction', 'model_reduction',
9✔
25
           'minimal_realization', 'eigensys_realization', 'markov', 'hsvd',
26
           'balred', 'modred', 'minreal', 'era']
27

28

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

37
    Parameters
38
    ----------
39
    sys : `StateSpace`
40
        State space system.
41

42
    Returns
43
    -------
44
    H : array
45
        List of Hankel singular values.
46

47
    See Also
48
    --------
49
    gram
50

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

59
    Examples
60
    --------
61
    >>> G = ct.tf2ss([1], [1, 2])
62
    >>> H = ct.hsvd(G)
63
    >>> H[0]
64
    np.float64(0.25)
65

66
    """
67
    # TODO: implement for discrete-time systems
68
    if (isdtime(sys, strict=True)):
5✔
69
        raise NotImplementedError("Function not implemented in discrete time")
70

71
    Wc = gram(sys, 'c')
5✔
72
    Wo = gram(sys, 'o')
5✔
73
    WoWc = Wo @ Wc
5✔
74
    w, v = np.linalg.eig(WoWc)
5✔
75

76
    hsv = np.sqrt(w)
5✔
77
    hsv = np.array(hsv)
5✔
78
    hsv = np.sort(hsv)
5✔
79
    # Return the Hankel singular values, high to low
80
    return hsv[::-1]
5✔
81

82

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

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

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

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

115
    Returns
116
    -------
117
    rsys : `StateSpace`
118
        Reduced order model.
119

120
    Raises
121
    ------
122
    ValueError
123
        If `method` is not either 'matchdc' or 'truncate'.
124
    NotImplementedError
125
        If the 'matchdc' method is used for a discrete-time system.
126

127
    Warns
128
    -----
129
    UserWarning
130
        If eigenvalues of `sys.A` are not all stable.
131

132
    Examples
133
    --------
134
    >>> G = ct.rss(4)
135
    >>> Gr = ct.model_reduction(G, [0, 2], method='matchdc')
136
    >>> Gr.nstates
137
    2
138

139
    See Also
140
    --------
141
    balanced_reduction, minimal_realization
142

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

151
    States, inputs, and outputs can be specified using integer offsets or
152
    using signal names.  Slices can also be specified, but must use the
153
    Python `slice` function.
154

155
    """
156
    if not isinstance(sys, StateSpace):
9✔
157
        raise TypeError("system must be a StateSpace system")
×
158

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

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

179
        elim = np.atleast_1d(_expand_key(elim))
9✔
180
        keep = np.atleast_1d(_expand_key(keep))
9✔
181

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

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

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

207
    B1 = sys.B[keep_states, :]
9✔
208
    B2 = sys.B[elim_states, :]
9✔
209

210
    C1 = sys.C[:, keep_states]
9✔
211
    C2 = sys.C[:, elim_states]
9✔
212

213
    # Figure out the new state space system
214
    if method == 'matchdc' and A22.size > 0:
9✔
215
        if sys.isdtime(strict=True):
9✔
216
            raise NotImplementedError(
217
                "'matchdc' not (yet) supported for discrete-time systems")
218

219
        # if matchdc, residualize
220
        # Check if the matrix A22 is invertible
221
        if np.linalg.matrix_rank(A22) != len(elim_states):
9✔
222
            raise ValueError("Matrix A22 is singular to working precision.")
×
223

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

233
        Ar = A11 - A12 @ A22I_A21
9✔
234
        Br = B1 - A12 @ A22I_B2
9✔
235
        Cr = C1 - C2 @ A22I_A21
9✔
236
        Dr = sys.D - C2 @ A22I_B2
9✔
237

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

245
    else:
246
        raise ValueError("Oops, method is not supported!")
×
247

248
    # Get rid of additional inputs and outputs
249
    Br = Br[:, keep_inputs]
9✔
250
    Cr = Cr[keep_outputs, :]
9✔
251
    Dr = Dr[keep_outputs, :][:, keep_inputs]
9✔
252

253
    rsys = StateSpace(Ar, Br, Cr, Dr)
9✔
254
    return rsys
9✔
255

256

257
def balanced_reduction(sys, orders, method='truncate', alpha=None):
9✔
258
    """Balanced reduced order model of system of a given order.
259

260
    States are eliminated based on Hankel singular value.  If `sys` has
261
    unstable modes, they are removed, the balanced realization is done on
262
    the stable part, then reinserted in accordance with [1]_.
263

264
    References
265
    ----------
266
    .. [1] C. S. Hsu and D. Hou, "Reducing unstable linear control
267
       systems via real Schur transformation".  Electronics Letters,
268
       27, 984-986, 1991.
269

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

287
    Returns
288
    -------
289
    rsys : `StateSpace`
290
        A reduced order model or a list of reduced order models if orders is
291
        a list.
292

293
    Raises
294
    ------
295
    ValueError
296
        If `method` is not 'truncate' or 'matchdc'.
297
    ImportError
298
        If slycot routine ab09ad, ab09md, or ab09nd is not found.
299
    ValueError
300
        If there are more unstable modes than any value in orders.
301

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

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

324
    # Check for ss system object, need a utility for this?
325

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

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

342
    rsys = []                   # empty list for reduced systems
5✔
343

344
    # check if orders is a list or a scalar
345
    try:
5✔
346
        iter(orders)
5✔
347
    except TypeError:           # if orders is a scalar
5✔
348
        orders = [orders]
5✔
349

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

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

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

381

382
def minimal_realization(sys, tol=None, verbose=True):
9✔
383
    """Eliminate uncontrollable or unobservable states.
384

385
    Eliminates uncontrollable or unobservable states in state-space
386
    models or canceling pole-zero pairs in transfer functions. The
387
    output `sysr` has minimal order and the same response
388
    characteristics as the original model `sys`.
389

390
    Parameters
391
    ----------
392
    sys : `StateSpace` or `TransferFunction`
393
        Original system.
394
    tol : real
395
        Tolerance.
396
    verbose : bool
397
        Print results if True.
398

399
    Returns
400
    -------
401
    rsys : `StateSpace` or `TransferFunction`
402
        Cleaned model.
403

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

411

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

417
    H = np.zeros((q*m, p*n))
9✔
418

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

424
    return H
9✔
425

426

427
def eigensys_realization(arg, r, m=None, n=None, dt=True, transpose=False):
9✔
428
    r"""eigensys_realization(YY, r)
429

430
    Calculate ERA model based on impulse-response data.
431

432
    This function computes a discrete-time system
433

434
    .. math::
435

436
        x[k+1] &= A x[k] + B u[k] \\\\
437
        y[k] &= C x[k] + D u[k]
438

439
    of order :math:`r` for a given impulse-response data (see [1]_).
440

441
    The function can be called with 2 arguments:
442

443
    * ``sysd, S = eigensys_realization(data, r)``
444
    * ``sysd, S = eigensys_realization(YY, r)``
445

446
    where `data` is a `TimeResponseData` object, `YY` is a 1D or 3D
447
    array, and r is an integer.
448

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

472
    Returns
473
    -------
474
    sys : `StateSpace`
475
        State space model of the specified order.
476
    S : array
477
        Singular values of Hankel matrix. Can be used to choose a good `r`
478
        value.
479

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

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

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

504
    q, p, l = YY.shape
9✔
505

506
    if m is None:
9✔
507
        m = 2*r
9✔
508
    if n is None:
9✔
509
        n = 2*r
9✔
510

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

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

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

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

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

532
    return StateSpace(Ar, Br, Cr, Dr, dt), S
9✔
533

534

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

538
    Calculate Markov parameters [D CB CAB ...] from data.
539

540
    This function computes the the first `m` Markov parameters [D CB CAB
541
    ...] for a discrete-time system.
542

543
    .. math::
544

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

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

553
    The function can be called with either 1, 2 or 3 arguments:
554

555
    * ``H = markov(data)``
556
    * ``H = markov(data, m)``
557
    * ``H = markov(Y, U)``
558
    * ``H = markov(Y, U, m)``
559

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

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

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

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

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

610
    """
611

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

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

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

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

657
    t = 0
9✔
658
    if truncate:
9✔
659
        t = m
9✔
660

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

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

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

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

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

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

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

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

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

727
# Function aliases
728
hsvd = hankel_singular_values
9✔
729
balred = balanced_reduction
9✔
730
modred = model_reduction
9✔
731
minreal = minimal_realization
9✔
732
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

© 2026 Coveralls, Inc