A method to construct reduced-order parameter ... - Semantic Scholar

INTERNATIONAL JOURNAL OF ROBUST AND NONLINEAR CONTROL Int. J. Robust. Nonlinear Control 2015; 00:1–30 Published online in Wiley InterScience (www.interscience.wiley.com). DOI: 10.1002/rnc

A method to construct reduced-order parameter varying models Jennifer Annoni∗ and Peter Seiler Aerospace Engineering and Mechanics, University of Minnesota

SUMMARY This paper describes a method to construct reduced-order models for high dimensional nonlinear systems. It is assumed that the nonlinear system has a collection of equilibrium operating points parameterized by a scheduling parameter. First, a reduced-order linear system is constructed at each equilibrium point using input/output data. This step combines techniques from dynamic mode decomposition and direct subspace system identification. This yields discrete-time models that are linear from input to output but whose state matrices are functions of the scheduling parameter. Second, a parameter varying linearization is used to connect these linear models across the various operating points. The key technical issue in this second step is to ensure the reduced-order model has a consistent state definition across all operating points. These two steps yield a reduced-order Linear Parameter Varying (LPV) system that approximates the nonlinear system c 2015 John Wiley & Sons, Ltd. even when the operating point changes in time. Copyright Received . . .

KEY WORDS: dynamic mode decomposition; reduced-order modeling; parameter varying

1. INTRODUCTION This paper describes a method to construct a reduced-order model for a high dimensional nonlinear system. The work is motivated by the control of systems that involve fluid and/or structural dynamics. One specific example is wind farm control. The overall performance of a wind farm can be improved through proper coordination of the turbines [1]. High fidelity computational fluid dynamic models have been developed for wind farms [2, 3]. These high fidelity models are accurate but are not suitable for control law design due to their computational complexity. A second example involves the control of flexible aircraft. More fuel efficient aircraft can be designed by reducing structural weight thus leading to increased flexibility [4]. The vibrational modes can significantly degrade the performance and can even lead to aeroelastic instabilities (flutter) [4, 5]. High fidelity, computational fluid/structural models also exist for this application [6,7] but these are also too complex for control design. Simplified, control-oriented models are needed in both of these examples. There are several relevant reduced-order modeling techniques developed by the fluids and controls communities. Proper Orthogonal Decomposition (POD), summarized in Section 2.1, is a standard method in the fluids community [8–10]. Briefly, the state is projected onto a low dimensional subspace of POD modes constructed using data from the high-order system. Dynamic mode decomposition (DMD) [11–13] is a more recent approach that fits time-domain data with linear dynamics on a reduced-order subspace. Both approaches are, in their basic form, for autonomous (unforced) systems although there are variations for input/output systems, e.g DMD with control [14] or IOROMs [7]. In the controls community, balanced truncation [15, 16] is a standard input/output model reduction technique. This requires the solution of Lyapunov equations ∗ Correspondence

to: University of Minnesota 110 Union St SE, Minneapolis, MN 55455 . E-mail: [email protected]

c 2015 John Wiley & Sons, Ltd. Copyright Prepared using rncauth.cls [Version: 2010/03/27 v2.00]

for controllability and observability Gramians. This is computationally costly for systems with extremely large state dimension. Balanced POD [17, 18] addresses this computational issue by forming approximate Gramians using POD modes from simulations of the system and a related adjoint. Finally, system identification techniques, e.g. subspace methods [19–21], can be used to identify reduced-order models from input/output data. This paper describes a new method to construct reduced-order linear parameter varying (LPV) models that approximate a high-order nonlinear model. The nonlinear system is assumed to have a parameterized collection of equilibrium operating points. For example, the free-stream wind speed parameterizes the equilibrium condition in a wind farm. For flexible aircraft, the dynamics depend on the aircraft airspeed and altitude. The proposed approach involves two steps. First, DMD is extended to handle inputs/outputs (Section 3). This new method, termed Input/Output DMD (IODMD), uses direct N4SID subspace identification [19] on a low-dimensional subspace to construct a reduced-order, linear model at one operating condition. Second, the reduced-order models constructed at fixed operating conditions are “stitched” together using a parameter varying linearization (Section 4). The key technical issue is that the states of the reduced-order model must have a consistent meaning across all operating conditions (as described in Section 4.3). The proposed approach handles this issue by constructing a single reduced-order subspace that is used at all operating conditions. The proposed IODMD and LPV linearization methods are demonstrated on a nonlinear mass-spring-damper example (Section 5). There are several benefits of the proposed reduced-order LPV models. The main benefit is that the models can be used for standard, gain-scheduled control. In addition, more formal control analysis and synthesis tools have been developed for LPV systems [22–25]. This LPV modeling approach can yield models that are accurate over a wide range of operating conditions for a nonlinear system. This is in contrast to the existing reduced-order modeling approaches described above which, for the most part, are used to construct a single linear model. In addition, the IODMD step relies on input/output data from a forced response and does not require linearization of the system or the construction/simulation of the system adjoint as in Balanced POD. Finally, the reduced-order model involves a projection onto a well-defined reduced-order subspace. This retains a physical meaning in the reduced-order states. In other words, the reduced-order state can be used to approximate the fullorder state of the system. This provides insight into the key spatial modes of fluid/structure systems and is not simply a black-box identification technique, often seen in subspace identification. 2. BACKGROUND This section summarizes the basic linear algebra results used in this paper. Details can be found in linear algebra textbooks, e.g. [26, 27]. First, every matrix M ∈ Rn×m has a singular value T n×n decomposition (SVD) given by M = U ΣV and V ∈ Rm×m are orthogonal. If h where iU ∈ R ¯ Σ n×m ¯ = diag(σ1 , . . . , σm ). The diagonal n ≥ m then Σ ∈ R has the form Σ = 0(n−m)×m where Σ ¯ are the singular values of M and can be ordered as σ1 ≥ · · · ≥ σm ≥ 0. The rank of M entries of Σ is equal to the number of nonzero singular values of M . If m > n , then the SVD takes a similar form with Σ = [ Σ¯ 0n×(m−n) ]. qP P n m 2 Next, the Frobenius norm for a matrix M ∈ Rn×m is defined as kM kF := i=1 j=1 |Mij | . The Frobenius norm can also be expressed as kM k2F = tr(M T M ) = tr(M M T ).  M1 using the trace 2 Note that if M is partitioned as M = M2 then kM kF = kM1 k2F + kM2 k2F . Similarly, if M = [ M3 M4 ] then kM k2F = kM3 k2F + kM4 k2F . It can also be shown that kM kF = kQM W kF for any orthogonal matrices Q ∈ Rn×n and W ∈ Rm×m . In other words, the Frobenius norm is invariant to multiplication by orthogonal matrices. This invariance property of the Frobenius norm will be applied to solve several matrix optimization problems. The first result stated below is a matrix least squares problem. Lemma 1 Let M ∈ Rn×m and N ∈ Rr×m be given matrices. Assume m ≥ r and N has rank r so that the

SVD has the form N = U ΣV T = U [ Σ¯ 0r×(m−r) ]

h

V1T V2T

i

¯ ∈ Rr×r is nonsingular. Then where Σ

min kM − Y N k2F = kM V2 k2F

(1)

Y ∈Rn×r

The unique solution is Yopt = M N † where N † = V

h

¯ −1 Σ

0(m−r)×r

i

U T is the Moore-Penrose inverse.

Proof For any choice of Y , the orthogonal invariance of the Frobenius norm implies

 ¯ kM − Y N k2F = k(M − Y N )V k2F = M V1 − Y U Σ

 2 M V2 F

(2)

¯ 2 + kM V2 k2 . The Frobenius norm can be split into two terms as: kM − Y N k2F = kM V1 − Y U Σk F F The second term is unaffected by the choice of Y . The first term can be made equal to zero by ¯ −1 U T = M N † . This optimal choice has approximation error the (unique) choice Yopt = M V1 Σ 2 kM V2 kF .

2.1. Proper Orthogonal Decomposition (POD) This section briefly reviews the use of POD to construct reduced-order dynamic models for autonomous (unforced) systems. The summary is written to align with the results developed in the remainder of the paper. More complete details on POD can be found in [8–10]. Consider an autonomous system modeled by the following discrete-time, nonlinear dynamics: xk+1 = f (xk )

(3)

where x ∈ Rnx is the state vector. POD is typically applied when the state dimension nx is extremely s large (> 1, 000). A collection of ns snapshot measurements {xk }nk=1 ⊂ Rnx is obtained from the system. The objective is to approximate the system on a low dimensional subspace. This requires the selection of the low dimensional subspace as well as the approximation of the snapshots on the chosen subspace. First, assume that the low-dimensional subspace has been chosen. Formally, an r-dimensional subspace is defined by a set of linearly independent basis vectors {qj }rj=1 ⊂ Rnx . Further assume the basis vectors are orthonormal: qjT qi = 0 if j 6= i and qiT qi = 1.† The snapshots are approximated by linear combinations of the basis vectors. That is, coefficients cjk ∈ R are computed such Pr that xk ≈ j=1 cjk qj . To simplify notation, define the matrices of snapshots, basis vectors, and coefficients as X := [x1 , . . . , xns ] ∈ Rnx ×ns , Q := [q1 , . . . , qr ] ∈ Rnx ×r , and C ∈ Rr×ns . Orthonormality of the basis vectors implies QT Q = Ir . In matrix notation, given a subspace defined by Q the objective is to find a coefficient matrix C such that X ≈ QC . The next lemma provides the optimal least-squares solution for C . Lemma 2 Let X ∈ Rnx ×ns be given snapshot data. Let Q ∈ Rnx ×r with QT Q = Ir span a given rdimensional subspace. Then min

C∈Rr×ns

kX − QCk2F = kX − QQT Xk2F

(4)

with unique minimizer given by Copt := QT X .

† The assumption that the basis vectors are orthonormal is without loss of generality. Specifically, given any r linearly independent vectors, Gram-Schmidt orthogonalization can be used to obtain r orthonormal vectors that span the same subspace.

Proof   Use Gram-Schmidt orthogonalization to construct a matrix Q⊥ ∈ Rnx ×(nx −r) such that Q Q⊥ is an nx × nx orthogonal matrix. The orthogonal invariance of the Frobenius norm implies

h T i

2

h T i 2

Q

kX − QCk2F = Q (X − QC) = QQX−C (5)

T T X F

⊥

⊥

F

T

The second equality follows from Q Q = Ir and QT⊥ Q = 0(nx −r)×r . The error can be split into two terms: kX − QCk2F = kQT X − Ck2F + kQT⊥ Xk2F . The second term is unaffected by the choice of C . The first term can be made equal to zero by the unique choice Copt = QT X . This yields an approximation error kQT⊥ Xk2F which corresponds to the projection of the snapshots onto the subspace orthogonal to that spanned by Q. This minimal error can also be expressed as kX − QCopt k2F = kX − QQT Xk2F . The optimal least squares approximation Copt = QT X is simply the projection of the snapshots onto the basis vectors. Lemma 2 assumes that the basis vectors Q are given. This raises the question as to which r-dimensional subspace yields the minimal projection error. This question is solved by the Proper Orthogonal Decomposition as formally stated in the next theorem. Theorem 1 Let X ∈ Rnx ×ns be given snapshot data with SVD X = U ΣV T . Let r be any non-negative integer such that r ≤ rank(X) and σr > σr+1 . Then rank(X)

min

Q∈R

nx ×r

T

, Q Q=Ir

C∈Rr×ns

kX −

QCk2F

=

X

σk2

(6)

k=r+1

The (unique) optimal r-dimensional subspace is spanned by Qopt = Ur where Ur contains the first r columns of U . The corresponding optimal projection coefficients are given by Copt = Σr VrT where Vr contains the first r columns of V and Σr := diag(σ1 , . . . , σr ). Proof This is a consequence of the Eckart-Young-Mirsky Theorem [28, 29]. Briefly, the orthogonal ˆ := U T QCV invariance of the Frobenius norm implies kX − QCkF = kΣ − U T QCV kF . Define Σ ˆ ˆ 2 ≥ and note that rank(Σ) ≤ r for any choice of Q and C . It can be shown that kΣ − Σk F Prank(X) 2 T ˆ k=r+1 σk for any rank(Σ) ≤ r . The choice of Qopt = Ur and Copt = Σr Vr achieves this lower bound. The columns of U are the POD modes of the snapshots in X . The energy lost by projecting the Prank(X) state onto the first r POD modes is k=r+1 σk2 . The lost energy relative to the total energy is given by: Prank(X) 2 σk (7) P OD := Pk=r+1 rank(X) 2 σk k=1 The number of modes r is typically selected to retain a certain percentage of the energy in the snapshot data, e.g. r is selected to achieve P OD < 0.01 or some similar tolerance. A reduced-order nonlinear model can be constructed by projecting the state onto the selected POD modes. Specifically, let Q := Ur ∈ Rnx ×r contain the r selected POD modes. Define the reduced state as z := QT x ∈ Rr . A reduced-order nonlinear model that approximates Equation 3 is given by zk+1 = QT f (Qzk ). Solutions of the reduced-order model zk can be used to construct approximate solutions of the full order model as xk ≈ Qzk . This procedure constructs a reduced-order model by using the Galerkin projection onto the POD modes of the snapshot data. Additional details on the Galerkin projection can be found in [30] with specific application to the Navier-Stokes equations for fluid flows.

3. INPUT/OUTPUT DYNAMIC MODE DECOMPOSITION This section describes a method to construct a reduced-order linear model for nonlinear systems with inputs and outputs. The approach projects data to a low-dimensional subspace and then applies the direct N4SID subspace identification method [19]. This differs from POD in that it uses snapshot data to construct both the low dimensional subspace as well as the model state-space matrices. For autonomous systems (no inputs) this proposed method reduces to the dynamic mode decomposition [11–13]. Hence the procedure is termed Input/Output Dynamic Mode Decomposition (IODMD). A related method is the DMD with control described in [14]. DMD with control uses a projection on both the state and input snapshot measurements while IODMD only uses a projection on the state measurements. The approaches are similar but it may be more intuitive to avoid projection on the inputs/outputs (as done here) if these dimensions are relatively small. Consider a discrete-time nonlinear system: xk+1 = f (xk , uk )

(8)

yk = h(xk , uk )

where x ∈ Rnx , u ∈ Rnu , and y ∈ Rny are the state, input, and output vectors. A linear model will be constructed using data collected near some equilibrium point. Specifically, assume that the system has an equilibrium condition described by (¯ x, u ¯, y¯) such that x ¯ = f (¯ x, u ¯)

(9)

y¯ = h(¯ x, u ¯)

If the state is initialized at x0 = x¯ and the input is held fixed at uk = u¯ for k = 0, 1, . . . then the state and output remain in equilibrium at xk = x¯ and yk = y¯ for k = 0, 1, . . .. A collection of ns snapshots are obtained by exciting the system near this equilibrium point. The snapshots include s s s measurements of the state {xk }nk=0 ⊂ Rnx , input {uk }nk=0 ⊂ Rnu , and output {yk }nk=0 ⊂ R ny . Define perturbations of these measurements from the equilibrium condition: δxk := xk − x ¯,

δuk := uk − u ¯,

δyk := yk − y¯

(10)

The proposed IODMD attempts to fit the snapshots with a linear time-invariant model of the form: δxk+1 = A δxk + B δuk δyk = C δxk + D δuk

(11)

The state matrices (A, B, C, D) are real matrices with dimensions compatible to those of (x, u, y). The constructed model is simply a best linear fit to the collected data and is not necessarily a linearization. Moreover, the equilibrium point assumption is not essential to the development here but it extends naturally in the next section to LPV models. The intent is to apply the method for systems with extremely large state dimension (nx > 1, 000) but with moderate input/output dimensions (nu , ny < 100). The method projects the state onto a low dimensional subspace in order to make the computations tractable. The problem is first simplified by assuming that an r-dimensional subspace of Rnx has been selected. An orthonormal basis for this subspace is specified by the columns of a matrix Q ∈ Rnx ×r with QT Q = Ir . It is possible to construct an additional nx − r orthonormal vectors that span the remaining dimensions of the state space, e.g. by using Gram-Schmidt orthogonalization. This yields a matrix Q⊥ ∈ Rnx ×(nx −r) such that Q Q⊥ is an nx × nx orthogonal matrix. The matrix   Q Q⊥ defines a (nonsingular) state-space coordinate transformation:     z = Qz + Q⊥ z⊥ (12) δx := Q Q⊥ z⊥ where z ∈ Rr and z⊥ ∈ Rnx −r . This decomposes the state as projections z := QT x and z⊥ := QT⊥ x onto the subspaces defined by Q and Q⊥ , respectively. The discrete-time dynamics in Equation 11

can be expressed in the new coordinates as: h T i h T i Q AQ QT AQ⊥ Q B z ] + δuk [ zz⊥ ]k+1 = Q [ T T z ⊥ k QT ⊥ AQ Q⊥ AQ⊥ ⊥B   z δyk = CQ CQ⊥ [ z⊥ ]k + D δuk

(13)

It is assumed that the dynamics can be well-approximated on the subspace defined by Q. As a result, the states z⊥ can be truncated with negligible error. The truncated, reduced-order model takes the form: zk+1 = (QT AQ)zk + (QT B)δuk := F zk + Gδuk

(14) (15)

δyk = (CQ)zk + Dδuk := Hzk + Dδuk

The state matrices of the reduced-order system have dimensions F ∈ Rr×r , G ∈ Rr×nu , and H ∈ Rny ×r . This is equivalent to the following low-rank approximation for the full-order state matrices:         A B QF QT QG Q 0 F G QT 0 ≈ = (16) C D 0 I ny H D 0 Inu HQT D The next theorem provides the optimal choice for the reduced-order state matrices (F, G, H, D) given a specific subspace spanned by Q. Theorem 2 Let the following snapshot data be given: X0 := [x0 , . . . , xns −1 ] ∈ Rnx ×ns , U0 := [u0 , . . . , uns −1 ] ∈ R

nu ×ns

,

X1 := [x1 , . . . , xns ] ∈ Rnx ×ns Y0 := [y0 , . . . , yns −1 ] ∈ R

(17)

ny ×ns

Let Q ∈ Rnx ×r with QT Q = Ir span a given r-dimensional subspace. Assume R(r+nu )×ns has rank equal to (r + nu ). Then

  h i

X1 Q 0 F

Y 0 − 0 In y H h i min F G ∈R(r+ny )×(r+nu ) H D

has a unique minimizer given by 

F H

G D

 opt

G D

 h QT

0

0 In u

i

X0 U0



2

(18) i QT X0 ∈ U0

(19)

F

 T   T † Q X1 Q X0 = Y0 U0

Proof The orthogonal invariance of the Frobenius norm implies

 T 

  h i

Q 0  X  h Q 0 i   h QT 0 i  X0 

X1

2 Q 0 F G 1 F

QT⊥ 0 =

Y0 − 0 I n y H D

U0 Y0 − 0 I n y H 0 Inu

0 I F ny

h

(20)

G D

2  h QT X i 0

U0

F

(21) Use QT Q = Ir and QT⊥ Q = 0(nx −r)×r in order to split the error into two terms:

h T i   h QT X i

Q X1

2 F G 0 − H

Y0

+ kQT⊥ X1 k2F D U0

(22)

F

The second term is unaffected by the choice of (F, G, H, D). By Lemmah1, the first i term is uniquely QT X0 minimized by the state matrices in Equation 20. The assumption that has rank equal to U0 (r + nu ) can be removed although the optimal state matrices will no longer be unique in this case.

Theorem 2 projects the state onto a subspace defined by Q. The optimal reduced-order state matrices (Equation 20) are then obtained by least squares. This is the direct N4SID subspace method [19] for estimating state matrices given data for the input, output, and reduced-order state. The error introduced by this approximation consists of two terms. The second term of Equation 22 is the error due to projecting the state snapshots onto the subspace defined by Q. The first term of Equation 22 evaluated with the optimal state matrices gives the error associated with fitting the projected snapshots with a linear system. This raises the question as to which r-dimensional subspace yields the minimum error. The optimal choice for Q must balance the two sources of error due to projection and linear fitting. A sub-optimal, but useful, choice for the projection subspace is given by the POD modes of X0 . Specifically, let the SVD of X0 be given by X0 = U ΣV T . Approximate the state on a subspace defined by the first r POD modes of X0 , i.e. Q := Ur . By Theorem 2, the optimal reduced-order state-matrices for this choice is:    T  † F G Ur X1 Σr Vr = (23) H D opt U0 Y0 This reduces to the standard DMD method if the system contains no inputs. Specifically, this yields the following low-rank approximation (assuming no inputs): A ≈ Ur Fopt UrT

(24)

The low-rank model provides information about the dynamic modes of the system. Specifically, let Fopt have an eigenvalue decomposition T ΛT −1 where T ∈ Cr×r and Λ = diag(λ1 , . . . , λr ) has the (possibly complex) eigenvalues of Fopt . The matrix T := [t1 , . . . , tr ] contains the corresponding eigenvectors, i.e. Fopt tj = λj tj . The DMD modes are given by mj = Ur tj for j = 1, . . . , r. By Equation 24, these modes satisfy Amj = λj mj . Thus if δx0 = mj then δxk = λkj mj for k = 0, 1, . . .. In other words, each mode mj provides spatial information regarding a specific temporal frequency λj for the system. The IODMD procedure constructs state matrices (Equation 23) for systems with inputs. As with standard DMD, an eigenvalue decomposition of Fopt can be used to construct DMD modes that provide spatial modes associated with a specific temporal frequency. IODMD also yields input/output information for the model. This proposed method is a tractable implementation of the direct N4SID subspace method and can be applied for very large systems. In particular the pseudoinverse and matrix multiplications in Equation 23 involve matrices with dimensions (r, nu , ny , ns ). This is not simply a black-box (input-output) approach because the reduced-order state zk can be used to approximately reconstruct the full order state by δxk ≈ Ur zk . The approach only requires input/output/state data from the model. Construction and simulation of an adjoint system, as in Balanced POD, is not required.

4. PARAMETER VARYING DMD 4.1. Parameter Varying Linearizations Consider a discrete-time nonlinear system: xk+1 = f (xk , uk , ρk ) yk = h(xk , uk , ρk )

(25)

where x ∈ Rnx , u ∈ Rnu , and y ∈ Rny are the state, input, and output vectors. The model also depends on a time-varying parameter ρ ∈ Rnρ . The parameter describes the operating condition. The IODMD approach can be extended to this case using a parameter varying linearization similar to the one in [31]. Assume there is a collection of equilibrium points (¯ x(ρ), u ¯(ρ), y¯(ρ)) such that x ¯(ρ) = f (¯ x(ρ), u ¯(ρ), ρ) y¯(ρ) = h(¯ x(ρ), u ¯(ρ), ρ)

(26)

If the state is initialized at x0 = x¯(ρ), the input is held fixed at uk = u¯(ρ), and the operating condition is frozen at ρk = ρ then the state and output will remain in equilibrium at xk = x¯(ρ) and yk = y¯(ρ) for k = 0, 1, . . .. Thus (¯ x(ρ), u ¯(ρ), y¯(ρ)) defines an equilibrium condition for each fixed value of ρ. The nonlinear dynamics can be linearized around the equilibrium points defined for each fixed value of ρ. Define perturbations from the equilibrium condition as: δxk := xk − x ¯(ρ),

δuk := uk − u ¯(ρ),

δyk := yk − y¯(ρ)

(27)

For a fixed operating condition (ρk = ρ for k = 0, 1, . . .), a standard linearization yields: δxk+1 = A(ρ)δxk + B(ρ)δuk δyk = C(ρ)δxk + D(ρ)δuk

where the linearized state matrices are defined by ∂f ∂f , B(ρ) := A(ρ) := ∂x (¯x(ρ),¯u(ρ),ρ) ∂u (¯x(ρ),¯u(ρ),ρ) ∂h ∂h C(ρ) := , D(ρ) := ∂x (¯x(ρ),¯u(ρ),ρ) ∂u (¯x(ρ),¯u(ρ),ρ)

(28)

(29)

Next consider the case where the operating condition, specified by the parameter ρk , varies in time. In general (xk , uk , yk ) := (¯ x(ρk ), u ¯(ρk ), y¯(ρk )) is not a valid solution of the nonlinear system. In other words, the parameterized state/input/output values only define an equilibrium condition for fixed values of ρ. Despite this fact, it is still possible to construct a time varying linearization around the parameterized values (¯ x(ρ), u ¯(ρ), y¯(ρ)). Re-define perturbation variables as follows for the case where ρ varies in time: δxk := xk − x ¯(ρk ),

δuk := uk − u ¯(ρk ),

δyk := yk − y¯(ρk )

(30)

A valid (time-varying) linearization can be obtained when ρk is time-varying as follows: δxk+1 = f (xk , uk , ρk ) − x ¯(ρk+1 )

(31)

A Taylor series expansion of f (xk , uk , ρk ) yields: f (¯ xk (ρk ) + δxk , u ¯k (ρk ) + δuk , ρk ) ≈ x ¯(ρk ) + A(ρk )δxk + B(ρk )δuk

(32)

where A and B are as defined in Equation 29. A similar Taylor series approximation can be performed to linearize the output function h in terms of the matrices C and D defined in Equation 29. Combining Equations 31 and 32 leads to the following parameter varying linearization: δxk+1 = A(ρk )δxk + B(ρk )δuk + (¯ x(ρk ) − x ¯(ρk+1 )) δyk = C(ρk )δxk + D(ρk )δuk

(33)

This differs from the standard linearization at a single fixed operating point in two respects. First, ρk varies in time and hence this is a time-varying system. More precisely, the time variations in the state matrices (A, B, C, D) arise due to ρk and hence this is called a linear parameter varying (LPV) system. There are many tools in the controls literature that address this class of systems [22–25]. Second, the equilibrium state varies in time due to the changing operating condition. This effect is retained by the term x¯(ρk ) − x¯(ρk+1 ) which provides a forcing term for the dynamics. This model linearizes the dependence on the state and the input. The linearization approximation is accurate as long as the state, input, and output remain near the manifold of equilibrium points (¯ x(ρk ), u ¯(ρk ), y¯(ρk )). It is important that the nonlinear dependence on the operating condition (specified by ρk ) is retained.

4.2. Example: Nonlinear Mass-Spring-Damper (1-Block) A simple nonlinear mass-spring-damper example is presented to illustrate the parameter varying linearization method. A single mass, shown in Figure 1, is connected to a wall by a linear damper and nonlinear spring. Two forces F and u are applied to the block. Newton’s second law yields the following continuous-time, nonlinear model:     d p v (34) = 1 dt v m (F + u − bv − k(p)p) where p is the position of the block (units of m) and v is the velocity (m/s). b is the (linear) damping constant (in N/(m/s)), and k(p) is a nonlinear spring constant (in N/m). The spring constant has the form k(p) = k1 + k2 p2 and is assumed to be a stiffening spring k2 > 0. The corresponding values for these system constants are given in Figure 1. Finally, the block has two forces: an exogenous (disturbance) force F and a controllable input force u (both in N).

Parameter m k1 k2 b

Value 1 0.5 1 1

Unit kg N/m N/m3 N/(m/s)

Figure 1. Nonlinear mass-spring-damper (Left) and system constants (Right).

A continuous-time parameter varying linearization can be performed. Instead, the model is discretized to align with the discrete-time development used in this the paper. A simple Euler approximation with sample time ∆t yields a two-state discrete-time model: xk+1 = f (xk , uk , Fk )

where x := [p, v]T is the state and the discretized dynamics are given by:   vk f (xk , uk , Fk ) := xk + ∆t 1 m (Fk + uk − bvk − k(pk )pk )

(35)

(36)

For simplicity, the outputs for this model will be the full state, i.e. y := x. The disturbance force Fk is treated as the varying parameter in the dynamics: ρk := Fk . If u ≡ 0 and the disturbance force is held constant at F ≡ ρ then the mass moves to an equilibrium position p¯(ρ). This equilibrium position can be determined by solving the cubic equation ρ = k (¯ p) p¯ for p¯. Thus a parameterized collection of equilibrium conditions for the mass spring damper system is: ¯ 
(¯ x(ρ), u ¯(ρ)) := p(ρ) ,0 (37) 0 The parameter varying linearization around this collection of equilibrium points is given by: δxk+1 = A(ρk )δxk + B(ρk )δuk + (¯ x(ρk ) − x ¯(ρk+1 )) (38) i h i h 0 1 0 where: A(ρ) := I2 + ∆t − klin (ρ) − b and B(ρ) := ∆t m1 . The linearized spring constant at m

m

the operating condition ρ is klin (ρ) := k1 + 3k2 p¯(ρ)2 . The spring stiffens as it is stretched and, as a result, the dynamics of the system change when the exogenous disturbance force F is applied. Figure 2 shows responses of the nonlinear system and the corresponding parameter varying (LPV) linearization. The parameter varying linearization is in discrete-time with a sample time of

∆t = 0.01sec. The external force is ρk := Fk = 0.8 cos(0.1k∆t) N. The initial condition is x(0) = [¯ p(ρ0 ), 0]T = [0.751, 0]T which corresponds to δx0 = 0. The controllable input is set to δuk = 0 for the entire simulation. Thus the only forcing arises due to the changing external disturbance force. This appears as a forcing term in the parameter varying linearization (Equation 38) due to the varying equilibrium term. The nonlinear responses were simulated in continuous-time. The left subplot of Figure 2 shows the position and velocity vs. time. The parameter varying linearization agrees closely with the true nonlinear response. There are some small discrepancies between the two velocity trajectories. The parameter varying linearization is accurate even though the exogenous disturbance force, position, and velocity all vary over a large range. The right subplot shows the position and velocity vs. exogenous force. This subplot also shows the equilibrium (trim) relation x ¯(ρ) given in Equation 37. The key point is that the position and velocity stay near this equilibrium manifold throughout the simulation. As a consequence, (δxk , δuk ) remain small and the Taylor series approximation is accurate. If the frequency of ρk is increased then the trajectories will deviate more significantly from the equilibrium manifold. As a consequence, faster changes in the exogenous force lead to larger linearization errors.

Figure 2. System response vs. time (left) and vs. parameter (right)

4.3. State Consistency Issue Parameter varying linearizations can be constructed using data from steady operating conditions specified by constant values of ρ. Specifically, the linearization only requires state matrices (A(ρ), B(ρ), C(ρ), D(ρ)) and equilibrium conditions (¯ x(ρ), u ¯(ρ), y¯(ρ)) at each fixed value of ρ. In principal, the proposed IODMD or another model reduction method can be used to identify

reduced-order state matrices at each ρ. One key technical issue is state consistency. To clarify this issue, consider an autonomous system without inputs and outputs. Let A(ρ) denote the state matrix that appears in the parameter varying linearization. DMD can be used at each fixed value of ρ to identify a subspace spanned by Q(ρ) ∈ Rnx ×r(ρ) and a reduced-order matrix F (ρ) ∈ Rr(ρ)×r(ρ) such that A(ρ) ≈ Q(ρ)F (ρ)Q(ρ)T . The reduced-order state at the operating point ρ is defined as z := Q(ρ)T δx ∈ Rr(ρ) . This reduced-order state z will lack consistency if the projection subspace Q(ρ) depends on the parameter. In other words, the meaning of z and dimension r(ρ) will depend on the value of ρ. Hence the state z at ρ1 will not, in general, be consistent in either meaning or dimension with the state z at ρ2 6= ρ1 . To circumvent this issue, a single consistent subspace Q ∈ Rnx ×r should be used at all operating conditions. The reduced-order state z := QT δx ∈ Rr then has a consistent meaning for all parameter values. Moreover, a reduced-order matrix F (ρ) ∈ Rr×r can be identified at each value of ρ such that A(ρ) ≈ QF (ρ)QT . This would lead to a reduced-order parameter varying linearization of the form: zk+1 = F (ρk )zk + (¯ z (ρk ) − z¯(ρk+1 ))

(39)

where z¯(ρ) := QT x¯(ρ) is the equilibrium state projected onto the reduced subspace. The next section presents a method to address this state consistency issue. 4.4. Reduced-Order Parameter Varying Linearizations The proposed approach, shown in Algorithm 1, constructs a reduced-order model for a discrete-time nonlinear system (Equation 25) in four steps. First, a collection of parameter values are selected (Line 1). Second, data is collected from the nonlinear system at these selected parameter values (Lines 2-5). The collected data includes the equilibrium conditions as well as state/input/output snapshots obtained by exciting the nonlinear system. The algorithm, as written, assumes that the same number of snapshots ns are obtained at each grid point. However, the number of snapshots can easily change with each grid point. Third, a single r-dimensional subspace of the state space is constructed (Lines 6-9). The subspace, defined by an orthogonal matrix Q ∈ Rnx ×r , is constructed from the POD modes of the snapshots obtained from all parameter values. Fourth, reduced-order state matrices are computed at each parameter vector using the IODMD approach described in Section 3 (Lines 10-11). The outcome of Algorithm 1 is a single r-dimensional subspace Q as well as equilibrium conditions and reduced-order state matrices computed at the selected grid points. This yields a reduced-order parameter varying linearization of the nonlinear system of the form: zk+1 = F (ρk )zk + G(ρk )δuk + (¯ z (ρk ) − z¯(ρk+1 )) δyk = H(ρk )zk + D(ρk )δuk

(44)

where z := QT δx ∈ Rr is the reduced-order state and z¯(ρ) := QT x¯(ρ) ∈ Rr is the reduced-order equilibrium state. The single subspace defined by Q is used to construct state matrices at all parameter values. Hence the reduced order state z retains a consistent meaning across the parameter domain. Note that Algorithm 1 only computes the state matrices and equilibrium conditions on a grid of specified parameter values. Interpolation (e.g. linear, spline, etc.) must be used to evaluate the state matrices and equilibrium conditions at any parameter value not contained in this grid. It is assumed that the grid of parameter values is sufficiently dense that this interpolation is accurate. The parameter varying linearization can be used to approximate the response of the nonlinear system for an initial condition x0 , input uk , and parameter trajectory ρk . Specifically, the initial condition and input for the nonlinear system define a corresponding initial condition z0 = QT (x0 − x ¯(ρ0 )) and input δuk = uk − u ¯(ρk ) for the parameter varying system. The reduced-order, parameter varying linearization (Equation 44) can be simulated to obtain the state response zk and output δyk . These correspond to the the state response xk = Qzk + x¯(ρk ) and output yk = δyk + y¯(ρk ) for the full-order nonlinear system. The subspace construction step of Algorithm 1 (Lines 6-9) requires the SVD of the matrix X0 that contains the snapshot data from all the operating points. This matrix has ns ng columns and hence the

Algorithm 1 Reduced-Order Parameter Varying Linearization 1:

n

g Given: Collection of parameter grid points {pj }j=1 .

Data Collection: At each grid point j = 1, . . . , ng do the following: Equilibrium: Compute the equilibrium condition (¯ x(pj ), u ¯(pj ), y¯(pj )) at ρ = pj . Experiment: Excite the nonlinear system (Equation 25) with fixed ρk = pj . The initial condition x0 and input uk should be near the equilibrium condition (¯ x(pj ), u ¯(pj )). 5: Snapshots: Define the matrices of snapshot deviations from the equilibrium conditions at pj :

2: 3: 4:

X0 (pj ) := [x0 − x ¯(pj ), . . . , xns −1 − x ¯(pj )] ∈ Rnx ×ns j

j

j

j

j

j

j

X1 (p ) := [x1 − x ¯(p ), . . . , xns − x ¯(p )] ∈ R

nx ×ns

j

U0 (p ) := [u0 − u ¯(p ), . . . , uns −1 − u ¯(p )] ∈ R j

Y0 (p ) := [y0 − y¯(p ), . . . , yns −1 − y¯(p )] ∈ R 6: 7: 8: 9: 10: 11:

nu ×ns

ny ×ns

(40) (41) (42) (43)

Construct Subspace for State Reduction: Stack Data: Define matrix of all state data: X0 := [X0 (p1 ), . . . , X0 (png )] ∈ Rnx ×(ns ng ) . POD Modes: Compute POD modes of X0 (Theorem 1). Subspace: Choose r modes, denoted Q ∈ Rnx ×r , to capture sufficient energy in X0 . Reduced-Order State Matrices: At each grid point j = 1, . . . , ng do the following: IODMD: Use snapshot data (X0 (pj ), X1 (pj ), U0 (pj ), Y0 (pj )) and subspace Q to compute state matrices (F (pj ), G(pj ), H(pj ), D(pj )) via IODMD (Theorem 2).

Algorithm 2 Iterative Method to Construct Subspace n

g Given: Collection of parameter grid points {pj }j=1 and snapshots X0 (pj ) from each grid point. 2: Initial Point: Use standard POD (Theorem 1) to compute r1 modes Q(p1 ) ∈ Rnx ×r1 to capture sufficient energy in X0 (p1 ). 3: Iterative Processing: For j = 2, . . . , ng , iteratively compute additional modes at each grid point. Given Q1 := [Q(p1 ), . . . Q(pj−1 )], use iterative POD (Theorem 3) to compute rj additional modes Q(pj ) ∈ Rnx ×rj to capture sufficient energy in X0 (pj ). 4: Subspace: Stack modes to form a single subspace defined by Q := [Q(p1 ), . . . , Q(png )].

1:

SVD of X0 may be computationally intractable if there are many grid points. A suboptimal, but more computationally efficient approach, is to iteratively process the snapshot data from each grid point. The basic idea is to determine a modes Q(p1 ) that capture the energy in the snapshots at the first grid point. Next, additional modes Q(p2 ) are computed so that the combined set [Q(p1 ), Q(p2 )] captures the energy in the snapshots at the second grid point. The procedure continues iteratively computing new modes Q2 := Q(pj ) to combine with previously computed modes Q1 := [Q(p1 ), . . . , Q(pj−1 )]. The benefit is that only snapshots obtained from one grid point are required for the calculations. This iterative procedure requires a method to compute the optimal (new) modes Q2 that should be added to some given modes Q1 . Theorem 3 in the appendix provides a POD-type result to perform this iterative calculation. Algorithm 2 gives the detailed steps for the iterative subspace construction. This iterative method can replace the single-step method described in Lines 6-9 of Algorithm 1. The remaining steps of Algorithm 1 are unchanged even when combined with the iterative subspace calculation.

5. EXAMPLE: NONLINEAR MASS-SPRING-DAMPER (M-BLOCKS) This section provides an example using both the IODMD and parameter-varying DMD approaches to construct low-order models. This extends the nonlinear mass-spring-damper (Section 4.2) to consist of M blocks. Each block has mass m and is connected to its neighboring block by a linear damper and nonlinear spring. In addition, each block is connected to “ground” (fixed wall) by a linear damper and nonlinear spring. All dampers have constant b and all springs are given by k(p) = k1 + k2 p2 . Finally, the M th block has an exogenous (disturbance) force F and a controllable input force u (both in N). Let (pi , vi ) denote the position and velocity of block i. The continuous-time dynamics for this M -block mass-spring-damper system are described by: mv˙ i = −Fi − Fi,i−1 − Fi,i+1 mv˙ M = F + u − FM − FM,M −1

(i = 1, . . . , M − 1)

(45)

Fi denotes the force due to a spring and damper connecting block i to the fixed wall. Fi,j denotes the force due to a spring and damper connecting blocks i and j . These forces are defined as Fi := bvi + k(pi )pi and Fi,j := b(vi − vj ) + k(pi − pj ) (pi − pj ). The convention for the first block is F1,0 ≡ 0. The position of the M th block is the output for the system, i.e. y = pM . All mass, spring, and damping constants are the same as those given in Figure 1. The result is a nonlinear system with state x := [pT , v T ]T ∈ R2M , a single controllable input u and single output y . The exogenous force F is again considered as a parameter in this model (ρ := F ). All results are given for M = 100 blocks. This is small enough that the reduced and full order models can be directly compared to assess the feasibility of the method.

5.1. Example: IODMD The IODMD method described in Section 3 is used to construct a single time-invariant model with the exogenous force held constant at ρ ≡ 1N. The system is simulated with ρ ≡ 1 and u = 0 to determine the equilibrium point x¯(ρ) for this parameter value. Next, the nonlinear system is excited in continuous-time for Tf = 20 sec with a chirp excitation input u(t) = 0.1 sin(ω(t)t). The chirp frequency is ω(t) = 0.1 (20)t/Tf which excites the nonlinear system with frequencies from 0.1 to 2 rad/sec. The nonlinear (continuous-time) system is simulated and snapshots of the state/output are collected every ∆t = 0.1 seconds. This yields 200 snapshots over the Tf = 20 second simulation. The left subplot of Figure 3 shows the POD and least-squares fitting errors obtained by the IODMD procedure. The horizontal axis is the number of POD modes selected from the snapshots X0 . The POD error is the relative energy lost in X0 by selecting the specified number of POD (as defined in Equation 7). The model relative fitting error is obtained  evaluating the IODMD  by 1 least-squares cost (Equation 19) with the optimal fit and dividing by k X Y0 kF . Both relative errors decrease rapidly with increasing number of POD modes. The fitting error levels off around 4 modes hence this was selected for the model fit. The relative fit error does not go to zero because the discrete-time linear system cannot perfectly fit the snapshots from the continuous-time nonlinear system. The right subplot shows a Bode plot for the four-state, reduced-order, discrete time model. The plot also shows a 200-state, full order, continuous time linearization. The reduced-order model agrees quite closely with the full order model in the frequency band of input excitation. The mismatch between the two models at high frequencies can be reduced by reducing the snapshot sampling time ∆t and increasing the chirp excitation frequency. A time-domain step response was used as another validation of the reduced order model. Specifically, the left subplot of Figure 4 shows the output of the discrete-time, reduced-order fourstate model with a step input δuk = 0.1. The exogenous force remains constant at F = ρ = 1. The vertical axis shows the position of the M th block relative to the trim position. The plot also shows the response of the full order, continuous time nonlinear system. The nonlinear system is initialized at the trim condition x(0) = x¯(ρ). The response of the reduced-order model is in close agreement with the nonlinear model. Moreover, the states of the reduced-order model zk retain a physical significance. These can be “lifted” to estimate the response of the full order system using the chosen

POD modes (δxk ≈ Qzk ). The right subplot of Figure 4 compares a subset of these state estimates (position of blocks 98-100) for the reduced-order system with the true response of the nonlinear system. Again there is close agreement although some discrepancies appear in the response for block 98. This error is due to the lack of input excitation. Specifically, the input force u is applied to block 100 and the chosen chirp amplitude provides less excitation for blocks located further from block 100. Finally, it is important to note that the quality of the reduced order model degrades if too many POD modes (i.e. much more than 4) are selected for the fit. This is due to overfitting a linear model to a nonlinear response. Inverse modeling techniques may be used to select the number of modes to accurately represent the system without overfitting [32]. This will be addressed in future work. The POD and model fit errors shown in Figure 3 provide useful criteria to avoid such overfitting.

Figure 3. IODMD Results: POD/Fit Relative Errors (left) and Bode Plot (right)

5.2. Example: Parameter Varying DMD Next, the parameter varying DMD method described in Section 4 is used to construct parameter varying reduced order model. A grid of eleven parameter values was selected for the exogenous forces as {0, 0.2, . . . , 2}. At each fixed value of ρ = pj in this grid (j = 1, . . . 11) the system is simulated with u = 0 to determine the corresponding equilibrium point x¯(pj ). Next, the nonlinear system is excited for Tf = 20 sec at each pj with a chirp excitation input u(t) = 0.1 sin(ω(t)t). The exogenous force F = ρ stretches the springs and this increases the system speed of response due to the spring stiffening. Thus the chirp frequency was chosen to depend on the parameter value as t/T ω(t) = 0.1 (10 + 10ρ) f . This excites the nonlinear system with frequencies from 0.1 to (1 + ρ)

Figure 4. IODMD Results: Output vs. Time (left) and States vs. Time (right)

rad/sec. The nonlinear (continuous-time) system is simulated and snapshots of the state/output are collected every ∆t = 0.1 seconds. This yields 200 snapshots over the Tf = 20 second simulation for each parameter grid point. The basic single-step procedure (Algorithm 1) is used to construct the subspace modes Q. This example is sufficiently small such that it was possible to compute the SVD on snapshots obtained at all grid points. Five modes were selected and reduced order models were constructed at each grid point using the IODMD procedure. The model relative fitting error was less than 0.05 at each grid point. This indicates that the reduced order models at each “fixed” grid point accurately match the recorded snapshot data. The iterative method (Algorithm 2) was also implemented. This also yielded a five-state subspace Q although with slightly larger fitting errors. The single-step and iterative approach gave similar results and hence the remainder of the example focuses on the single-step algorithm. Figure 5 compares a time-domain step response of the full order nonlinear system and the reduced order LPV model. The exogenous disturbance varies as F (t) = ρ(t) = − cos(0.5t) + 1. The controllable input is a step u(t) = 0 for t < 25sec and u(t) = 0.1 for t ≥ 25sec. The nonlinear system was simulated in continuous-time and the reduced order parameter varying model was simulated in discrete time. Linear interpolation was used to compute the state matrices and equilibrium point appearing in the LPV model for parameter values not contained in the 11 point grid. The left subplot of Figure 5 shows the position of the M th block. This is obtained from the parameter varying model as yk = y¯(ρk ) + δyk . The response of the reduced order model is in close agreement with the nonlinear model. The effect of the step input u is apparent after t = 25sec. Again, the states of the reduced model zk can be “lifted” to estimate the response of the full order system using the chosen POD modes as xk = x¯(ρk ) + δxk ≈ x¯(ρk ) + Qzk . The right subplot of Figure 5

compares a subset of these state estimates (position of blocks 98-100) for the reduced order system with the true response of the nonlinear system. Again, there is close agreement.

Figure 5. LPV DMD Results: Output vs. Time (left) and States vs. Time (right)

6. CONCLUSIONS This paper described a method to construct reduced-order models for high dimensional nonlinear systems. It is assumed that the nonlinear system has a collection of equilibrium operating points. The method has two main components. First, a reduced-order linear system is constructed at each equilibrium point using input/output data. Second, a parameter varying linearization is used to connect these linear models. A nonlinear mass-spring-damper example was used to demonstrate the method. Future work will apply this method to construct reduced order models for wind farm control.

7. ACKNOWLEDGMENTS This work was partially supported by the National Science Foundation under Grant No. NSFCMMI-1254129 entitled “CAREER: Probabilistic Tools for High Reliability Monitoring and Control of Wind Farms”. The work was also supported by the University of Minnesota Institute on the Environment, IREE Grant No. RS-0039-09.

REFERENCES 1. Johnson K, Thomas N. Wind farm control: Addressing the aerodynamic interaction among wind turbines. American Control Conference, 2009; 2104–2109. 2. Churchfield M, Lee S. NWTC design codes-SOWFA 2012. Http://wind.nrel.gov/designcodes/simulators/SOWFA. 3. Yang X, Sotiropoulos F, Conzemius R, Wachtler J, Strong M. Large-eddy simulation of turbulent flow past wind turbines/farms: the virtual wind simulator (VWiS). Wind Energy 2014; . 4. Beranek J, Nicolai L, Buonanno M, Burnett E, Atkinson C, Holm-Hansen B, Flick P. Conceptual design of a multiutility aeroelastic demonstrator. 13th AIAA/ISSMO Multidisciplinary Analysis Optimization Conference, 2010; AIAA–2010–9350. 5. Lind R, Brenner M. Robust aeroservoelastic stability analysis. Springer-Verlag, 1999. 6. MSC Software Corporation. MSC/NASTRAN: Structural and Multidiscipline Finite Element Analysis 2012. 7. Danowsky B, Thompson P, Farhat C, Lieu T, Harris C, Lechniak J. Incorporation of feedback control into a highfidelity aeroservoelastic fighter aircraft model. Journal of Aircraft 2010; 47(4):1274–1282. 8. Lo´eve M. Probability Theory. Van Nostrand, 1955. 9. Berkooz G, Holmes P, Lumley J. The proper orthogonal de-composition in the analysis of turbulent flows. Annual Review of Fluid Mechanics 1993; 25:539–575. 10. Holmes P, Lumley J, Berkooz G. Turbulence, Coherent Structures, Dynamical Systems and Symmetry. Cambridge University Press, 1988. 11. Schmid P. Dynamic mode decomposition of numerical and experimental data. Journal of Fluid Mechanics 2010; 656:5–28. 12. Schmid P, Li L, Juniper M, Pust O. Applications of the dynamic mode decomposition. Theoretical and Computational Fluid Dynamics 2010; 25:249–259. 13. Tu J, Rowley C, Luchtenburg D, Brunton S, Kutz J. On dynamic mode decomposition: Theory and applications. submitted to the Journal of Computational Dynamics 2013; . 14. Proctor J, Brunton S, Kutz J. Dynamic mode decomposition with control. arXiv:1409.6358 2014; . 15. Moore B. Principal component analysis in linear systems: Controllability, observability, and model reduction. IEEE Transactions on Automatic Control 1981; 26(1):17–32. 16. Pernebo L, Silverman L. Model reduction via balanced state space representations. IEEE Transactions on Automatic Control 1982; 27(2):382–387. 17. Willcox K, Peraire J. Balanced model reduction via the proper orthogonal decomposition. AIAA journal 2002; 40(11):2323–2330. 18. Rowley C. Model reduction for fluids, using balanced proper orthogonal decomposition. International Journal of Bifurcation and Chaos 2005; 15(03):997–1013. 19. Viberg M. Subspace-based methods for identification of linear time-invariant systems. Automatica 1995; 31(12):1835–1851. 20. van Overschee P, De Moor B. Subspace Identification for Linear Systems: Theory, Implementation, and Applications. Kluwer, 1996. 21. Ljung L. System Identification: Theory for the User. Second edn., Prentice Hall, 1999. 22. Packard A. Gain scheduling via linear fractional transformations. Systems and Control Letters 1994; 22:79–92. 23. Wu F. A generalized LPV system analysis and control synthesis framework. International Journal of Control 2001; 74:745–759. 24. Scherer C. Advances in linear matrix inequality methods in control, chap. Robust mixed control and linear parameter-varying control with full-block scalings. SIAM, 2000; 187–207. 25. Pfifer H, Seiler P. Robustness analysis of linear parameter varying systems using integral quadratic constraints. International Journal of Robust and Nonlinear Control 2014; doi:10.1002/rnc.3240. 26. Horn R, Johnson C. Matrix Analysis. Cambridge University Press, 1990. 27. Golub G, Loan CV. Matrix Computations. Johns Hopkins University Press, 1996. 28. Eckart G, Young G. The approximation of one matrix by another of lower rank. Psychometrika 1936; 1:211–218. 29. Markovsky I. Low-Rank Approximation: Algorithms, Implementation, Applications. Springer, 2012. 30. Holmes P, Lumley J, Berkooz G. Turbulence, coherent structures, dynamical systems and symmetry. Cambridge university press, 1998. 31. Takarics B, Seiler P. Gain scheduling for nonlinear systems via integral quadratic constraints. accepted to the American Control Conference, 2015. 32. Aster R, Borchers B, Thurber C. Parameter estimation and inverse problems. Academic Press, 2011.

Theorem 3 Let Q1 ∈ Rnx ×r1 be a given matrix with QT1 Q1 = Ir1 . Let X ∈ Rnx ×ns be given snapshot data. Define the SVD of the projected snapshot matrix (Inx − Q1 QT1 )X = U ΣV T . Let r2 be any nonnegative integer such that r2 ≤ rank((Inx − Q1 QT1 )X) and σr2 > σr2 +1 . Then rank(X)

min

Q2 ∈Rnx ×r2 , QT 2 Q2 =Ir2 C1 ∈Rr1 ×ns r2 ×ns C2 ∈R

kX − Q1 C1 −

Q2 C2 k2F

=

X

σk2

(46)

k=r2 +1

An optimal solution is given by C1,opt = QT1 X , Q2,opt = Ur , and C2,opt = Σr VrT where Σr , Ur , and Vr are associated with the first r singular values and vectors of (Inx − Q1 QT1 )X .

Proof orthogonalization to construct a matrix Q1,⊥ ∈ Rnx ×(nx −r1 ) such that Use Gram-Schmidt  Q1 Q1,⊥ is orthogonal. The orthogonal invariance of the Frobenius norm thus implies

h T i

2

h T i

Q

Q X−C1 −QT1 Q2 C2 2 kX − Q1 C1 − Q2 C2 k2F = QT1 (X − Q1 C1 − Q2 C2 ) = Q1T X−Q

T Q C 1,⊥

F

1,⊥

1,⊥

2

2

(47)

F

The second equality follows from QT1 Q1 = Ir1 and QT1,⊥ Q1 = 0(nx −r1 )×r1 . The error can be split as: kX − Q1 C1 − Q2 C2 k2F = kQT1 X − C1 − QT1 Q2 C2 k2F + kQT1,⊥ X − QT1,⊥ Q2 C2 k2F

(48)

The second term is unaffected by the choice of C1 . Moreover, for any (Q2 , C2 ) the first term can be made equal to zero by the choice C1,opt = QT1 X − QT1 Q2 C2 . In fact, QT1 Q2 = 0 may be assumed without loss of generality. Specifically, the choice of Q2 only affects the second term of Equation 48 (assuming the optimal choice for C1 just specified). a change of h ˜Perform i   h ˜2 i Q2 variables Q2 = Q1 Q1,⊥ Q . This change of variables from Q to is invertible since 2 ˆ2 ˆ2 Q Q   Q1 Q1,⊥ is orthogonal. Substitute this change of variables into the second error term to obtain ˆ 2 C2 k2 . Thus Q ˜ 2 has no effect on the second term and kQT1,⊥ X − QT1,⊥ Q2 C2 k2F = kQT1,⊥ X − Q F ˆ 2 and hence QT Q2 = 0. In this case, may be set to zero. Q2 can be selected to have the form Q1,⊥ Q 1 T the optimal choice for C1 simplifies to C1,opt = Q1 X . Retaining the assumption that QT1 Q2 = 0 as well as C1,opt = QT1 X , the total error is given by: kX − Q1 C1 − Q2 C2 k2F = k(Inx − Q1 QT1 )X − Q2 C2 k2F

(49)

By the standard POD result (Theorem 1) this cost is minimized by the choice Q2,opt = Ur , and C2,opt = Σr VrT where Σr , Ur , and Vr are associated with the first r singular values and vectors of (Inx − Q1 QT1 )X . It can be shown that this construction satisfies the assumption QT1 Q2,opt = 0.