Skip to content

Commit 77b25d5

Browse files
murrayrmmurrayrm
committed
initial refactoring of model_reduction (no changes in functionality)
1 parent e6824d7 commit 77b25d5

1 file changed

Lines changed: 88 additions & 102 deletions

File tree

‎control/modelsimp.py‎

Lines changed: 88 additions & 102 deletions
Original file line numberDiff line numberDiff line change
@@ -1,52 +1,20 @@
1-
#! TODO: add module docstring
21
# modelsimp.py - tools for model simplification
32
#
43
# Author: Steve Brunton, Kevin Chen, Lauren Padilla
54
# Date: 30 Nov 2010
65
#
76
# This file contains routines for obtaining reduced order models
87
#
9-
# Copyright (c) 2010 by California Institute of Technology
10-
# All rights reserved.
11-
#
12-
# Redistribution and use in source and binary forms, with or without
13-
# modification, are permitted provided that the following conditions
14-
# are met:
15-
#
16-
# 1. Redistributions of source code must retain the above copyright
17-
# notice, this list of conditions and the following disclaimer.
18-
#
19-
# 2. Redistributions in binary form must reproduce the above copyright
20-
# notice, this list of conditions and the following disclaimer in the
21-
# documentation and/or other materials provided with the distribution.
22-
#
23-
# 3. Neither the name of the California Institute of Technology nor
24-
# the names of its contributors may be used to endorse or promote
25-
# products derived from this software without specific prior
26-
# written permission.
27-
#
28-
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
29-
# "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
30-
# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
31-
# FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL CALTECH
32-
# OR THE CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
33-
# SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
34-
# LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF
35-
# USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
36-
# ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
37-
# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
38-
# OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
39-
# SUCH DAMAGE.
40-
#
41-
# $Id$
8+
9+
import warnings
4210

4311
# External packages and modules
4412
import numpy as np
45-
import warnings
46-
from .exception import ControlSlycot, ControlArgument, ControlDimension
47-
from .iosys import isdtime, isctime
48-
from .statesp import StateSpace
13+
14+
from .exception import ControlArgument, ControlDimension, ControlSlycot
15+
from .iosys import isctime, isdtime
4916
from .statefbk import gram
17+
from .statesp import StateSpace
5018
from .timeresp import TimeResponseData
5119

5220
__all__ = ['hankel_singular_values', 'balanced_reduction', 'model_reduction',
@@ -108,35 +76,48 @@ def hankel_singular_values(sys):
10876
return hsv[::-1]
10977

11078

111-
def model_reduction(sys, ELIM, method='matchdc', warn_unstable=True):
112-
"""Model reduction by state elimination.
113-
114-
Model reduction of `sys` by eliminating the states in `ELIM` using a given
115-
method.
79+
def model_reduction(
80+
sys, elim_states=None, method='matchdc', elim_inputs=None,
81+
elim_outputs=None, keep_states=None, keep_inputs=None,
82+
keep_outputs=None, remove_hidden_states='unobs', 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, or states from the original system. The
87+
specific states, inputs, or outputs that are eliminated can be
88+
specified by listing either the states, inputs, or outputs to be
89+
eliminated or those to be kept. In addition, unused states (those that
90+
do not affect the inputs or outputs of the reduced system) can also be
91+
eliminated.
11692
11793
Parameters
11894
----------
11995
sys : StateSpace
12096
Original system to reduce.
121-
ELIM : array
122-
Vector of states to eliminate.
97+
elim_inputs, elim_outputs, elim_states : array of int or str, optional
98+
Vector of inputs, outputs, or states to eliminate. Can be specified
99+
either as an offset into the appropriate vector or as a signal name.
100+
keep_inputs, keep_outputs, keep_states : array, optional
101+
Vector of inputs, outputs, or states to keep. Can be specified
102+
either as an offset into the appropriate vector or as a signal name.
103+
remove_hidden_states : str, optional
104+
If `unobs` (default), then eliminate any states that are unobservable
105+
via the reduced-order system outputs.
123106
method : string
124-
Method of removing states in `ELIM`: either 'truncate' or
107+
Method of removing states in `elim`: either 'truncate' or
125108
'matchdc' (default).
126109
warn_unstable: bool, option
127110
If `False`, don't warn if system is unstable.
128111
129112
Returns
130113
-------
131114
rsys : StateSpace
132-
A reduced order model.
115+
Reduced order model.
133116
134117
Raises
135118
------
136119
ValueError
137120
If `method` is not either ``'matchdc'`` or ``'truncate'``.
138-
NotImplementedError
139-
If `sys` is a discrete time system.
140121
141122
Warns
142123
-----
@@ -146,55 +127,53 @@ def model_reduction(sys, ELIM, method='matchdc', warn_unstable=True):
146127
Examples
147128
--------
148129
>>> G = ct.rss(4)
149-
>>> Gr = ct.modred(G, [0, 2], method='matchdc')
130+
>>> Gr = ct.model_reduction(G, [0, 2], method='matchdc')
150131
>>> Gr.nstates
151132
2
152133
153134
"""
135+
if not isinstance(sys, StateSpace):
136+
raise TypeError("system must be a a StateSpace system")
154137

155-
# Check for ss system object, need a utility for this?
138+
# Check system is stable
139+
if warn_unstable:
140+
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):
142+
warnings.warn("System is unstable; reduction may be meaningless")
156143

157-
# TODO: Check for continous or discrete, only continuous supported for now
158-
# if isCont():
159-
# dico = 'C'
160-
# elif isDisc():
161-
# dico = 'D'
162-
# else:
163-
if (isctime(sys)):
164-
dico = 'C'
165-
else:
166-
raise NotImplementedError("Function not implemented in discrete time")
144+
# 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]
167148

168-
# Check system is stable
169-
if np.any(np.linalg.eigvals(sys.A).real >= 0.0) and warn_unstable:
170-
warnings.warn("System is unstable; reduction may be meaningless")
171-
172-
ELIM = np.sort(ELIM)
173-
# Create list of elements not to eliminate (NELIM)
174-
NELIM = [i for i in range(len(sys.A)) if i not in ELIM]
175-
# A1 is a matrix of all columns of sys.A not to eliminate
176-
A1 = sys.A[:, NELIM[0]].reshape(-1, 1)
177-
for i in NELIM[1:]:
178-
A1 = np.hstack((A1, sys.A[:, i].reshape(-1, 1)))
179-
A11 = A1[NELIM, :]
180-
A21 = A1[ELIM, :]
181-
# A2 is a matrix of all columns of sys.A to eliminate
182-
A2 = sys.A[:, ELIM[0]].reshape(-1, 1)
183-
for i in ELIM[1:]:
184-
A2 = np.hstack((A2, sys.A[:, i].reshape(-1, 1)))
185-
A12 = A2[NELIM, :]
186-
A22 = A2[ELIM, :]
187-
188-
C1 = sys.C[:, NELIM]
189-
C2 = sys.C[:, ELIM]
190-
B1 = sys.B[NELIM, :]
191-
B2 = sys.B[ELIM, :]
149+
# Determine which states to keep
150+
elim_states, keep_states = _process_elim_or_keep(
151+
elim_states, keep_states, sys.state_labels)
152+
153+
keep_inputs = slice(None, None)
154+
keep_outputs = slice(None, None)
155+
156+
# Create submatrix of states we are keeping
157+
A11 = sys.A[:, keep_states][keep_states, :] # states we are keeping
158+
A12 = sys.A[:, elim_states][keep_states, :] # needed for 'matchdc'
159+
A21 = sys.A[:, keep_states][elim_states, :]
160+
A22 = sys.A[:, elim_states][elim_states, :]
192161

162+
B1 = sys.B[keep_states, :]
163+
B2 = sys.B[elim_states, :]
164+
165+
C1 = sys.C[:, keep_states]
166+
C2 = sys.C[:, elim_states]
167+
168+
# Figure out the new state space system
193169
if method == 'matchdc':
194-
# if matchdc, residualize
170+
if sys.isdtime(strict=True):
171+
raise NotImplementedError(
172+
"'matchdc' not (yet) supported for discrete time systems")
195173

174+
# if matchdc, residualize
196175
# Check if the matrix A22 is invertible
197-
if np.linalg.matrix_rank(A22) != len(ELIM):
176+
if np.linalg.matrix_rank(A22) != len(elim_states):
198177
raise ValueError("Matrix A22 is singular to working precision.")
199178

200179
# Now precompute A22\A21 and A22\B2 (A22I = inv(A22))
@@ -210,22 +189,29 @@ def model_reduction(sys, ELIM, method='matchdc', warn_unstable=True):
210189
Br = B1 - A12 @ A22I_B2
211190
Cr = C1 - C2 @ A22I_A21
212191
Dr = sys.D - C2 @ A22I_B2
192+
213193
elif method == 'truncate':
214194
# if truncate, simply discard state x2
215195
Ar = A11
216196
Br = B1
217197
Cr = C1
218198
Dr = sys.D
199+
219200
else:
220201
raise ValueError("Oops, method is not supported!")
221202

203+
# Get rid of additional inputs and outputs
204+
Br = Br[:, keep_inputs]
205+
Cr = Cr[keep_outputs, :]
206+
Dr = Dr[keep_outputs, :][:, keep_inputs]
207+
222208
rsys = StateSpace(Ar, Br, Cr, Dr)
223209
return rsys
224210

225211

226212
def balanced_reduction(sys, orders, method='truncate', alpha=None):
227213
"""Balanced reduced order model of sys of a given order.
228-
214+
229215
States are eliminated based on Hankel singular value.
230216
If sys has unstable modes, they are removed, the
231217
balanced realization is done on the stable part, then
@@ -280,7 +266,7 @@ def balanced_reduction(sys, orders, method='truncate', alpha=None):
280266
raise ValueError("supported methods are 'truncate' or 'matchdc'")
281267
elif method == 'truncate':
282268
try:
283-
from slycot import ab09md, ab09ad
269+
from slycot import ab09ad, ab09md
284270
except ImportError:
285271
raise ControlSlycot(
286272
"can't find slycot subroutine ab09md or ab09ad")
@@ -350,7 +336,7 @@ def balanced_reduction(sys, orders, method='truncate', alpha=None):
350336

351337
def minimal_realization(sys, tol=None, verbose=True):
352338
""" Eliminate uncontrollable or unobservable states.
353-
339+
354340
Eliminates uncontrollable or unobservable states in state-space
355341
models or cancelling pole-zero pairs in transfer functions. The
356342
output sysr has minimal order and the same response
@@ -382,14 +368,14 @@ def _block_hankel(Y, m, n):
382368
"""Create a block Hankel matrix from impulse response."""
383369
q, p, _ = Y.shape
384370
YY = Y.transpose(0, 2, 1) # transpose for reshape
385-
371+
386372
H = np.zeros((q*m, p*n))
387-
373+
388374
for r in range(m):
389375
# shift and add row to Hankel matrix
390376
new_row = YY[:, r:r+n, :]
391377
H[q*r:q*(r+1), :] = new_row.reshape((q, p*n))
392-
378+
393379
return H
394380

395381

@@ -435,7 +421,7 @@ def eigensys_realization(arg, r, m=None, n=None, dt=True, transpose=False):
435421
unit-area impulse response of python-control. Default is True.
436422
transpose : bool, optional
437423
Assume that input data is transposed relative to the standard
438-
:ref:`time-series-convention`. For TimeResponseData this parameter
424+
:ref:`time-series-convention`. For TimeResponseData this parameter
439425
is ignored. Default is False.
440426
441427
Returns
@@ -470,7 +456,7 @@ def eigensys_realization(arg, r, m=None, n=None, dt=True, transpose=False):
470456
YY = np.array(arg, ndmin=3)
471457
if transpose:
472458
YY = np.transpose(YY)
473-
459+
474460
q, p, l = YY.shape
475461

476462
if m is None:
@@ -480,14 +466,14 @@ def eigensys_realization(arg, r, m=None, n=None, dt=True, transpose=False):
480466

481467
if m*q < r or n*p < r:
482468
raise ValueError("Hankel parameters are to small")
483-
469+
484470
if (l-1) < m+n:
485471
raise ValueError("not enough data for requested number of parameters")
486-
472+
487473
H = _block_hankel(YY[:, :, 1:], m, n+1) # Hankel matrix (q*m, p*(n+1))
488474
Hf = H[:, :-p] # first p*n columns of H
489475
Hl = H[:, p:] # last p*n columns of H
490-
476+
491477
U,S,Vh = np.linalg.svd(Hf, True)
492478
Ur =U[:, 0:r]
493479
Vhr =Vh[0:r, :]
@@ -504,7 +490,7 @@ def eigensys_realization(arg, r, m=None, n=None, dt=True, transpose=False):
504490

505491
def markov(*args, m=None, transpose=False, dt=None, truncate=False):
506492
"""markov(Y, U, [, m])
507-
493+
508494
Calculate the first `m` Markov parameters [D CB CAB ...] from data.
509495
510496
This function computes the Markov parameters for a discrete time
@@ -583,7 +569,7 @@ def markov(*args, m=None, transpose=False, dt=None, truncate=False):
583569
# Get the system description
584570
if len(args) < 1:
585571
raise ControlArgument("not enough input arguments")
586-
572+
587573
if isinstance(args[0], TimeResponseData):
588574
data = args[0]
589575
Umat = np.array(data.inputs, ndmin=2)
@@ -643,7 +629,7 @@ def markov(*args, m=None, transpose=False, dt=None, truncate=False):
643629
# This algorithm sets up the following problem and solves it for
644630
# the Markov parameters
645631
#
646-
# (l,q) = (l,p*m) @ (p*m,q)
632+
# (l,q) = (l,p*m) @ (p*m,q)
647633
# YY = UU @ H.T
648634
#
649635
# [ y(0) ] [ u(0) 0 0 ] [ D ]
@@ -679,7 +665,7 @@ def markov(*args, m=None, transpose=False, dt=None, truncate=False):
679665
# Truncate first t=0 or t=m time steps, transpose the problem for lsq
680666
YY = Ymat[:, t:].T
681667
UU = UUT[:, t:].T
682-
668+
683669
# Solve for the Markov parameters from YY = UU @ H.T
684670
HT, _, _, _ = np.linalg.lstsq(UU, YY, rcond=None)
685671
H = HT.T/dt # scaling

0 commit comments

Comments
 (0)