# Multimodeling Control via System Balancing

Authors: Nada Ratković Kovačević; Dobrila Škatarić
Journal: Mathematical Problems in Engineering
 (2010)

---

## Abstract

A new approach in multimodeling strategy is proposed. Multimodel strategies in which control agents use different simplified models of the same system are being developed using balancing transformation and the corresponding order reduction concepts. Traditionally, the multimodeling concept was studied using the ideas of multitime scales (singular perturbations) and weak subsystem coupling. For all reduced-order models obtained, a Linear Quadratic Gaussian (LQG) control problem was solved. Different order reduction techniques were compared based on the values of the optimized criteria for the closed-loop case where the full-order balanced model utilizes regulators calculated to be the optimal for various reduced-order models. The results obtained were demonstrated on a real-world example: a multiarea power system consisting of two identical areas, that is, two identical power plants.


---

## Body

## 1. Introduction

Large-scale systems have been the subject of research work for several decades [1–20]. Some order reduction techniques were developed for the singularly perturbed class of systems, based on different mathematical procedures, such as graph metric [21]. The concept of multimodel strategies for large-scale systems is originated from [17]. According to that concept, a large-scale system may be controlled by several independent agents using several simplified models of the system. In the singular perturbation [17], methodology has been used to develop multimodel strategies by exploiting the nature of the system that has two fast subsystems mutually weakly coupled and both strongly connected to the slow subsystem. The basic contribution of [17] has been to establish a set of conditions under which a multimodel strategy is well posed in the sense that the performance of the control system is close to the performance that would have been obtained had the control strategy been designed by a single controller knowing the exact model of the overall system. Multimodel strategies [17] have been studied and employed in several papers, either within the context of singular perturbations or in different mathematical set-ups [14]. Conditions for a multimodel strategy to be well-posed are investigated in [14] for a linear quadratic Gaussian (LQG) optimal control problem [22].The multimodeling structure used in the classical approach via singular perturbations and weak coupling is defined by a linear dynamic system that has one slow and two fast subsystems [11]. The fast subsystems are strongly connected to the slow subsystem and weakly connected (or not connected) among themselves. This structure describes well the dynamics of several real-world systems, for example, power systems [15, 17] and automobiles [2]. The corresponding multimodeling representation [11, 17] is defined by (1.1):
(1.1)[ẋ0(t)ε1ẋ1(t)ε2ẋ2(t)]=[A00A01A02A10A11ε3A12A20ε3A21A22][x0(t)x1(t)x2(t)]+[B01B02B11ε3B12ε3B21B22][u1(t)u2(t)],(1.2)y(t)=C[x0(t)x1(t)x2(t)]+D[u1(t)u2(t)],
where x0∈ℜn0 are slow state variables, x1∈ℜn1, x2∈ℜn2 are fast state variables, u1∈ℜm1, u2∈ℜm2 are control inputs, y∈ℜp are outputs, ε1, ε2 are small positive singular perturbation parameters, and ε3 is small weak-coupling parameter. (For ε3=0, (1.1) describes multiparameter singularly perturbed system (MSPS) studied in detail in literature [5–11].)For the purpose of deterministic optimal control of the above multimodeling structure, the quadratic performance criterion has to be minimized by the proper choice of the control variablesu1(t) and u2(t). The performance criterion for the linear quadratic Gaussian optimization [22] is given by
(1.3)J=12∫0+∞[xT(t)Qx(t)+uT(t)Ru(t)]dt,Q=QT≥0,R=RT>0,
where
(1.4)x(t)=[x0(t)x1(t)x2(t)],u(t)=[u1(t)u2(t)],Q=[Q00Q01Q02Q01TQ110Q02T0Q22],R=[R100R2],Q=qTq=[q01q110q020q22]T[q01q110q020q22]=[q01Tq01+q02Tq02q01Tq11q02Tq22q11Tq01q11Tq110q22Tq020q22Tq22].
In the general multimodeling case, all zero-elements in matrices R and Q can be replaced by O(ε) elements. (O(εi) is defined by O(εi)<cεi, where c is a bounded constant, i is a real number, and ε=∥[ε1ε2]∥.)To solve the multimodeling problem one proceeds with constructing two different models of (1.1), obtained by setting ε1=0, which leads to the first model for the first controller, and by setting ε2=0, which produces the second model for the second controller. In order to simplify equations, without loss of generality, small coupling parameter ε3 is set to zero.The fast dynamics of the other subsystem is approximated by an algebraic equation (the correspondingεi is set to zero).Two approximations of the original model derived from (1.1) are
(1.5)[ẋ0(t)ε1ẋ1(t)0]=[A00A01A02A10A110A200A22][x0(t)x1(t)x2(t)]+[B01B02B1100B22][u1(t)u2(t)],[ẋ0(t)0ε2ẋ2(t)]=[A00A01A02A10A110A200A22][x0(t)x1(t)x2(t)]+[B01B02B1100B22][u1(t)u2(t)].
The above equations can be rewritten as:
(1.6)[ẋ0(t)ε1ẋ1(t)]=[(A00-A02A22-1A20)A01A10A11][x0(t)x1(t)]+[B01(B02-A02A22-1B22)B110][u1(t)u2(t)],[ẋ0(t)ε2ẋ2(t)]=[(A00-A01A11-1A10)A02A20A22][x0(t)x2(t)]+[(B01-A01A11-1B11)B020B22][u1(t)u2(t)],
leading to two different models of the original system. The algebraic equations defined in (1.5) are used in (1.3) by each controller to form their own performance criterion. Such simplified criteria are optimized by each controller via the corresponding reduced-order model and the obtained control strategies form the multimodeling strategy. The multimodeling strategy is well posed if it is O(ε) close to the global optimal control strategy obtained by performing direct optimization on the original system and the original performance criterion, as shown in [14].


## 2. The Use of Balancing Transformation for  the System Order Reduction

Robust order reduction based on the use of balancing transformation has been described in [20]. Concisely, consideration is given to a linear, time invariant system as in (2.1):
(2.1)dx(t)dt=Ax(t)+Bu(t),x(t0)=x0,y(t)=Cx(t)+Du(t),
where x(t) is an n-dimensional state vector, u(t) is an m-dimensional input vector, and y(t) is a p-dimensional output vector.For a linear, time invariant system (2.1) the corresponding transfer function for the open-loop system is given by (2.2):
(2.2)G(s)=C(sI-A)-1B+D.It is assumed that the system (2.1) is asymptotically stable and that a G(s) is of minimal realization.Assumption.
A system is asymptotically stable, a pair(A,B) is controllable, and a pair (A,C) is observable.The controllability and observability Gramians of the original system (2.1) satisfy the algebraic equations of Lyapunov [15, 16]:
(2.3)PAT+AP+BBT=0,QA+ATQ+CTC=0.For a system that is controllable and observable, both controllability Gramian and observability Gramian are positive definite matrices,P>0, Q>0.The balancing transformation is applied on the space vector in order to achieve that the controllability and the observability Gramians are identical and diagonal, that is,(2.4)xb(t)=Tx(t),det(T)≠0⇒dxb(t)dt=Abxb(t)+Bbu(t),yb(t)=Cbxb(t)+Dbu(t)=y(t),(2.5)Ab=TAT-1,Bb=TB,Cb=CT-1,Db=D,(2.6)Pb=Qb=Σ=diag{σ1,σ2,…,σn},σ1≥σ2≥⋯≥σn>0,
where σi are known as the Hankel singular values (HSVs).Assuming that the original system is controllable and observable, a balanced system will also be both controllable and observable, since the similarity transformation preserves controllability and observability of the system that was transformed [15, 16]. Hence all σi are positive. Furthermore, both original and balanced systems are of minimal realization. The transfer function of the balanced system given by
(2.7)Gb(s)=Cb(sI-Ab)-1Bb+D=G(s)
stays unchanged thanks to a coordinate change through a nonsingular transformation. The balanced controllability and observability Gramians satisfy the following algebraic Lyapunov equations:
(2.8)ΣAbT+AbΣ+BbBbT=0⇔PbAbT+AbPb+BbBbT=0,ΣAb+AbTΣ+CbTCb=0⇔QbAb+AbTQb+CbTCb=0.The idea of the order reduction through balancing transformation can be linked with the canonical system decomposition. It was shown that the system’s modes that were either uncontrollable or unobservable did not appear in the system transfer function. In [15, 16] it is shown that the system modes that are both weakly controllable and weakly observable have little influence on the system dynamics; so they can be neglected.However, it was noticed that those modes which are weakly controllable and well observable cannot be neglected as can neither be neglected those modes that are well controllable and weakly observable. Let us assume that the balanced system (2.4)–(2.6) is partitioned in the following way:
(2.9)Ab=[A11A12A21A22],Bb=[B11B22],Cb=[C11C22],Db=D,Σ=[Σ100Σ2],Σ1=diag{σ1,σ2,…,σr},Σ2=diag{σr+1,σr+2,…,σn}.
Assuming that σr>σr+1, balanced truncationproduces a system of lower order, r, defined by
(2.10)dx1(t)dt=A11x1(t)+B11u(t),y(t)=C11x1(t)+Du(t),
and the corresponding transfer function of the reduced-order system is
(2.11)G11(s)=C11(sI-A11)-1B11+D.The reduced-order system attained in this way is both controllable and observable since all corresponding HSVs are positive. Furthermore, the reduced-order system is balanced and asymptotically stable. It was shown in literature, for example, [20] that the H∞ norm for the reduced-order system, obtained through the truncation procedure given above, satisfies the condition:
(2.12)∥G(s)-G11(s)∥∞≤2(σr+1+σr+2+⋯+σn).It was noticed that the reduced-order system obtained through the balanced truncation procedure gives very good approximation of the original system in the case of pulse input for both control signals (good spectra approximation on higher frequencies) but shows a considerable, steady-state error in the case of step input (poor spectra approximation on lower frequencies). This error is due to the fact that the original system and the reduced order system have different DC gains. Actually, after the above-described truncation through balancing transformation, most of the spectra on lower frequencies are kept and also some of the spectra on higher frequencies, but some of the spectra on lower frequencies are lost as well as most of the spectra on higher frequencies. By eliminating part of the spectra on lower frequencies (which occur in the neglected part of the system—state variablesx2(t)) we have produced gain that differs from the gain of the original system that was balanced. This discrepancy was eliminated in [16] where a technique of balanced residualization was proposed that produced an accurate (exact) DC gain and a very good spectra approximation on lower frequencies and sometimes even on middle frequencies. It should be noted that in [20] a residualisation technique was also used. An improved truncation method that preserves the exact DC gain value as in the original system is given in [16] and applied in [19].


## 3. Multimodeling via System Balancing

Two multimodeling order reduction concepts are mentioned above: the first one uses the small-parameter idea and the second one is based on balancing transformation, presented through two specific methods—balanced truncation and balanced residualization. Here, we will present the idea of multimodeling via system balancing, which is more general than multimodeling via singular perturbation since the latter requires special structures of the original model. Multimodeling via singular perturbation is performed on the assumption that two fast weakly connected subsystems, whose system matrices are invertible, are both strongly coupled to the slow subsystem while no such assumption is necessary for the multimodeling via system balancing.Let us consider a time invariant linear system represented by (3.1):
(3.1)[ẋ0(t)ẋ1(t)ẋ2(t)]=A[x0(t)x1(t)x2(t)]+B[u1(t)u2(t)],y(t)=C[x0(t)x1(t)x2(t)]+D[u1(t)u2(t)],
where
(3.2)A=[A00A01A02A10A11A12A20A21A22],B=[B01B02B11B12B21B22].
The balanced model is derived as in (2.4)-(2.5), that is,
(3.3)xb(t)=Tx(t),det(T)≠0⇒dxb(t)dt=Abxb(t)+Bbu(t),yb(t)=Cbxb(t)+Dbu(t)=y(t),Ab=TAT-1,Bb=TB,Cb=CT-1,Db=D,Pb=Qb=Σ=diag{σ1,σ2,…,σn},σ1≥σ2≥⋯≥σn>0.
The HSVs are the same as the HSVs of the original model (1.1) and are sorted in a descending order so that they can help in determining how to partition the model: the first few state space variables for the slow subsystem and the following two for the two fast parts, as in (3.4).The balanced model will be partitioned in the following manner (3.4)-(3.5):
(3.4)[ẋb0(t)ẋb1(t)ẋb2(t)]=[Ab00Ab01Ab02Ab10Ab11Ab12Ab20Ab21Ab22][xb0(t)xb1(t)xb2(t)]+[Bb01Bb02Bb11Bb12Bb21Bb22][u1(t)u2(t)],(3.5)yb(t)=[Cb0Cb1Cb2][xb0(t)xb1(t)xb2(t)]+[Db1Db2][u1(t)u2(t)]=y(t),
where xb0∈ℜn0 are slow or the common-core state variables of the balanced model; xb1∈ℜn1, xb2∈ℜn2 are state variables of the balanced model, having a small amount of total energy (n0+n1+n2=n); u1∈ℜm1, u2∈ℜm2 are control inputs of the balanced and original model (m1+m2=m); yb∈ℜp, y∈ℜp are outputs of the balanced and original model; respectively, and Abij,Bbik,Cbj, and Dbk (i,j=0,1,2 and k=1,2) are submatrices of the balanced system matrices Ab,Bb,Cb, and Db, having the corresponding dimensions.Instead of zero-submatrices in [12] a more general case is considered here, where submatrices Ab12, Ab21, Bb12, and Bb21 are nonzero matrices, as it is usually the case in the real-world systems.The corresponding quadratic performance criterion which has to be minimized is(3.6)Jb=12∫0+∞[xbT(t)Qbxb(t)+uT(t)Rbu(t)]dt,Qb=QbT≥0,Rb=RbT>0,
where, without loss of generality, we assume that Qb and Rb are block-diagonal matrices:
(3.7)Qb=[Qb00000Qb11000Qb22],Rb=[Rb1100Rb22].Now the decomposition of the balanced model (3.4) is performed into two reduced-order models, both of them composed of the common core—a subsystem corresponding to slow modes of the balanced model, and one of the two different remaining subsystems, corresponding to the modes having a small amount of total energy, namely, the first reduced-order model and the second reduced-order model. Cross-coupling matrices Ab12,Ab21 will be neglected though they have nonzero values. The first reduced-order model obtained is
(3.8)[ẋb0(t)ẋb1(t)]=A1r[xb0(t)xb1(t)]+B1r[u1(t)u2(t)],yb1(t)=[Cb0Cb1][xb0(t)xb1(t)]+Db1u1(t),(3.9)A1r=[Ab00Ab01Ab10Ab11],B1r=[Bb01Bb02Bb11Bb12].The second reduced-order model is
(3.10)[ẋb0(t)ẋb2(t)]=A2r[xb0(t)xb2(t)]+B2r[u1(t)u2(t)],yb2(t)=[Cb0Cb2][xb0(t)xb2(t)]+Db2u2(t),
where
(3.11)A2r=[Ab00Ab02Ab20Ab22],B2r=[Bb01Bb02Bb12Bb22].It is interesting to note thatthe first and the second reduced-order models in (3.8) and (3.10) both have HSVs corresponding to the common core modes of the original system and two complement sets of the HSVs corresponding to the modes of the original system that have a small amount of total energy, if we keep all parts of the Bb matrix.If, however, coupling submatricesBb12 and Bb21 are set to zero-matrices of the corresponding dimensions, as suggested in [11], the HSVs obtained for the first and the second subsystems will differ from those of the original full-order model, and the reduced-order models will be as in (3.12)–(3.13) and (3.14). These two types of reduced-order models would still contain all control signals of the original and balanced models.Thefirst reduced-order model is obtained as in [11, 17]:(3.12)[ẋb0(t)ẋb1(t)]=A1r[xb0(t)xb1(t)]+B1r[u1(t)u2(t)],yb1(t)=[Cb0Cb1][xb0(t)xb1(t)]+Db1u1(t),
where
(3.13)A1r=[Ab00Ab01Ab10Ab11],B1r=[Bb01Bb02Bb110].Thesecond reduced-order model is obtained as in [11, 17]:(3.14)[ẋb0(t)ẋb2(t)]=A2r[xb0(t)xb2(t)]+B2r[u1(t)u2(t)],yb2(t)=[Cb0Cb2][xb0(t)xb2(t)]+Db2u2(t).A2r=[Ab00Ab02Ab20Ab22],B2r=[Bb01Bb020bb22].The same approximation is done for the performance criterion (1.3); hence two performance criteria are needed, which leads to multicriteria optimization problem. Depending on the actual problem setup, very often described by differential games, the two controllers find their own optimal strategies and apply such strategies to the global system defined by (1.1). In such a way, the multimodeling strategy is well posed if the performance criterion under the multimodeling strategy is O(σr+1/σr) close to the global optimal control strategy obtained by performing direct optimization of the original performance criterion for the original system.The following two criteria to be optimized now are(3.15)J1=12∫0+∞{[xb0(t)xb1(t)]Q1[xb0(t)xb1(t)]+[u1(t)u2(t)]R1[u1(t)u2(t)]}dt,Q1=Q1T≥0,R1=R1T>0,(3.16)J2=12∫0+∞{[xb0(t)xb2(t)]Q2[xb0(t)xb2(t)]+[u1(t)u2(t)]R2[u1(t)u2(t)]}dt,Q2=Q2T≥0,R2=R2T>0.Q1 and Q2 are block-diagonal matrices and R1 and R2 are submatrices of Rb having corresponding dimensions:
(3.17)Q1=Q1T=[Qb0000Qb11]≥0,R1=R1T=Rb11>0,Q2=Q2T=[Qb0000Qb22]≥0,R2=R2T=Rb22>0.The optimal control problem we refer to, following the suggestions in [2, 14, 17], is to minimize a convex sum of J1 and J2 given in (3.15) and (3.16), that is, for some γ1 and γ2(3.18)JMMC=γ1J1+γ2J2,γ1+γ2=1,γ1>0,γ2>0.
This corresponds to a Pareto optimal cooperative strategy [6].In Pareto optimal strategy a situation is considered in which decision makers should decide on their strategies through mutual cooperation [14, 17]. The essence of this is that no variation from a Pareto optimal strategy can decrease the costs of either of the decision makers. Let each decision maker have a quadratic cost functional as in (3.15) and (3.16). A Pareto solution is a pair uopt1(t),uopt2(t) which minimizes (3.18) for some γ1 and γ2 [6]. The optimal feedback solution to (3.15) and (3.16) is given by (3.19):
(3.19)uopt(t)=[uopt1(t)uopt2(t)]=-Rb-1BbTPbxb(t).
Here Pb is the positive semidefinite stabilizing solution of the algebraic Riccati equation:
(3.20)AbTPb+PbAb+Qb-PbSbPb=0,Sb=BbRb-1BbT.The optimal state regulator is a special case of this problem where the decision makers agree on a choice ofγ1 and γ2 as weighting factors. Without prejudice to the generality, we chose γ1=γ2=0.5. We expect that the approximation of the optimization criteria would be
(3.21)12(J1+J2)≈Jb,
and according to [17] it is to be expected that:
(3.22)JMMC=12J1+12J2+O(σr+1σr)≈J.The original result is obtained in [9], where instead of O(σr+1/σr), O(∥μ∥2) was used, where μ stood for a norm of [ε1ε2ε3], in specific ∥μ∥=ε12+ε22+ε32. Still adhering to this idea but including criteria for balancing truncation, instead of O(∥μ∥2) we can use O(σr+1/σr) having in mind that the order r is chosen where σr+1 is of order of percent of the σr. In [6] a recursive algorithm is developed for solving multiparameter Riccati equations with the rate of convergence O(∥μ∥i) and the rate of accuracy for the near-optimal strategy is O(∥μ∥i+1) where i is the iteration number.The required solution of the algebraic Riccati equation (3.20) is based on the standard assumption [5, 12].Assumption.
The triples (A1r,B1r,chol(Q1)) and (A2r,B2r,chol(Q2)) are stabilizable-detectable. (Here chol(Q) is the Cholesky decomposition of a matrix Q.)Assumption.
The pairs(A1r,chol(Q1)) and (A2r,chol(Q2)) are detectable-observable.The matrices(Air,Bir), i=1,2 are given by (3.23):
(3.23)A1r=[Ab00Ab01Ab10Ab11],B1r=[Bb01Bb02Bb110],A2r=[Ab00Ab02Ab20Ab22],B2r=[Bb01Bb020Bb22].
The optimal feedback solution to (3.12)–(3.13) and (3.15) is given by
(3.24)uthefirstopt(t)=[ufirstopt1(t)ufirstopt2(t)]=-R1-1B1rTP1[xb0(t)xb1(t)],
where P1 is the positive semidefinite stabilizing solution of the corresponding algebraic Riccati equation (3.25):
(3.25)A1rTP1+P1A1r+Q1-P1S1P1=0,S1=B1rR1-1B1rT.
The optimal feedback solution to (3.14) and (3.16) is given by
(3.26)uthesecondopt(t)=[usecondopt1(t)usecondopt2(t)]=-R2-1B2rTP2[xb0(t)xb2(t)],
where P2 is the positive semidefinite stabilizing solution of the corresponding algebraic Riccati equation:
(3.27)A2rTP2+P2A2r+Q2-P2S2P2=0,S2=B2rR2-1B2rT.It is intuitively clear that a further reduction could be made for some models. If the control signals in the original model are weakly coupled, then it is possible to form two subsystems as described above, howeverthe first one having only inputs u1 and the second subsystem having only inputs u2. This would mean that the neglected modes having a small amount of total energy are considered to have reached their steady-state values, rather than as changing variables. Models of the type (1.1) are sound candidates for this type of reduction. This would be as if the submatrices Bb02 in (3.13) and Bb01 in (3.14) were set to zero-matrices of the corresponding dimensions, resulting in (3.28)–(3.31) and (3.32)–(3.35), respectively. The first reduced-order model obtained by neglecting the control variable u2 in (3.12)–(3.13) is
(3.28)[ẋb0(t)ẋb1(t)]=[Ab00Ab01Ab10Ab11][xb0(t)xb1(t)]+[Bb010Bb110][u1(t)u2(t)],
or
(3.29)[ẋb0(t)ẋb1(t)]=A1r[xb0(t)xb1(t)]+B1ru1(t),(3.30)yb1(t)=[Cb0Cb1][xb0(t)xb1(t)]+Db1u1(t),(3.31)A1r=[Ab00Ab01Ab10Ab11],B1r=[Bb01Bb11].The second reduced-order model obtained by neglecting the control variables u1 in (3.14) is presented in (3.32)–(3.35):
(3.32)[ẋb0(t)ẋb2(t)]=[Ab00Ab02Ab20Ab22][xb0(t)xb2(t)]+[0Bb020Bb22][u1(t)u2(t)],
or
(3.33)[ẋb0(t)ẋb2(t)]=A2r[xb0(t)xb2(t)]+B2ru2(t),(3.34)yb2(t)=[Cb0Cb2][xb0(t)xb2(t)]+Db2u2(t),(3.35)A2r=[Ab00Ab02Ab20Ab22],B2r=[Bb02Bb22].The same approximation is done for the performance criterion (3.15)-(3.16); hence two performance criteria are obtained, as in (3.36)-(3.37) and (3.38)-(3.39), each of them having its own control signal. Taking into consideration multicriteria optimization, both of these criteria would be regarded as having equal importance. Again, the multimodeling strategy will be considered to be well posed if the performance criterion under the multimodeling strategy is O(σr+1/σr) close to the global optimal control strategy obtained by performing direct optimization on the original system and the original performance criterion.The two criteria to be optimized now are(3.36)Jthefirst=12∫0+∞{[xb0(t)xb1(t)]Q1[xb0(t)xb1(t)]+u1(t)TR1u1(t)}dt,(3.37)Q1=Q1T=[Qb0000Qb11]≥0,R1=R1T=Rb11>0,(3.38)Jthesecond=12∫0+∞{[xb0(t)xb2(t)]Q2[xb0(t)xb2(t)]+u2(t)TR2u2(t)}dt,(3.39)Q2=Q2T=[Qb0000Qb22]≥0,R2=R2T=Rb22>0.In the simulation example it has been chosen thatQ1 and Q2 are unity matrices of the corresponding dimensions, while R1 and R2 are unity matrices of the corresponding dimensions multiplied by some weighing factor which is in accordance with the suggestions made in [17]:
(3.40)Q1=Q1T=In0+n1≥0,R1=R1T=Const·Im1>0,(3.41)Q2=Q2T=In0+n2≥0,R2=R2T=Const·Im2>0.
The optimal feedback solution to (3.36)-(3.37) and (3.40) is given by (3.42):
(3.42)uthefirstopt(t)=-R1-1B1rTPthefirst[xb0(t)xb1(t)],B1r=[Bb01Bb11],
where Pthefirst is the positive semidefinite stabilizing solution of the algebraic Riccati equation:
(3.43)A1rTPthefirst+PthefirstA1r+Q1-PthefirstSthefirstPthefirst=0,Sthefirst=B1rR1-1B1rT,
while the optimal feedback solution to (3.38)-(3.39) and (3.41) is given by:
(3.44)uthesecondopt(t)=-R2-1B2rTPthesecond[xb0(t)xb2(t)],B2r=[Bb02Bb22].
Here Pthesecond is the positive semidefinite stabilizing solution of the algebraic Riccati equation:
(3.45)A2rTPthesecond+PthesecondA2r+Q2-PthesecondSthesecondPthesecond=0,Sthesecond=B2rR2-1B2rT.
It could be explored at a later stage what the sufficient and necessary conditions should be for the existence of solutions to Riccati equations (3.43) and (3.45). Some useful ideas and efficient iterative solutions with the existence conditions could be found in the work of Mukaidani et al. [5–11].The next step is to implement optimal control from (3.42) and (3.44) to the full-order model. Optimal regulators (3.42) and (3.44) are easier to design than an optimal regulator for the full-order model (3.6). Since both control agents are of equal importance to the full-order model, approximation of the optimal control strategy could be
(3.46)uapproxopt(t)=[uthefirstopt(t)uthesecondopt(t)]=[-R1-1B1rTPthefirst[xb0(t)xb1(t)]-R2-1B2rTPthesecond[xb0(t)xb2(t)]],
as was shown in the simulation example.The optimal control problem we refer to, following the suggestions in [2, 14, 17], is to minimize a convex sum of Jthefirst and Jthesecond given in (3.36) and (3.38) for some γ1 and γ2:
(3.47)Japproximated=γ1Jthefirst+γ2Jthesecond,γ1+γ2=1,γ1>0,γ2>0.
This corresponds to a Pareto optimal cooperative strategy [6].In Pareto optimal strategy we have considered a situation in which decision makers decide on their strategies through mutual cooperation [14, 17]. The essence of this approach is that no variation from a Pareto optimal strategy can decrease the costs of any of the decision makers. Let each decision maker have a quadratic cost functional as in (3.36) and (3.38). A Pareto solution is a pair uthefirstopt(t),uthesecondopt(t) which minimized (3.47) for some γ1 and γ2 [6]. The optimal state regulator is a special aspect of this problem where the decision makers agree on the choice of γ1 and γ2 as weighting factors.Without prejudice to the generality of the above considerations, we choseγ1=γ2=0.5. According to [17] it is to be expected that the approximation of the optimization criteria would be
(3.48)Japproximated=Jbal.IandII=12(Jthefirst+Jthesecond)≈Jb.


## 4. Example

The methods for the order reduction displayed here were tested on an example known from literature and taken from [17] with some modification. Such example is a state space model of a power system consisting of two interconnected identical areas, where each area consists of one plant. The model behavior was simulated in the open-loop as well as in the closed-loop.For the model example a system was chosen having two inputs and one output. Based on HSVs a decision was made on what should be taken for the reduced order while the order of the slow or common-core subsystem and of the two fast parts (i.e., parts having a small amount of total energy) was chosen as suggested in [17]. Several reduced-order models were produced, using the methods as mentioned and described above.The efficiency of these approximations was compared first in the open-loop case for typical input functions: impulse, step, ramp, and sine. In the open-loop case a comparison was also made with the balanced model for all frequency characteristics (both magnitude and phase spectra) of all available transfer functions. Control strategy in this paper was different from the one used in [17]. Optimal linear quadratic regulators (LQRs) were designed for original and balanced models as well as for the ones of reduced order. Gain matrices of the reduced-order models were expanded by zeros to obtain the full order, and these reduced-order regulators were used to close the loop on the balanced model. Different approximations were compared with respect to the values of optimized LQG criteria.The model in the state-space representation is of the ordern=9 with the astatism of the second order therefore being unstable and is described by means of the matrices:(4.1)A=[004.50100000004.5-1000000-0.050-0.10.1000000-0.050.1000.100032.7-32.70000000000-550000-40000-10000000000-55000-400000-10],BT=[000.100010000000.1000010].Inspired by [1, page 227, Case 3] and [17] we have chosen
(4.2)C=[110001010],D=[00].For the model chosen for simulation, state-space variables as in [17] are as follows: x1 and v1 the integral of the area control error for the area 1, x2 and v2 the integral of the area control error for the area 2, x3 and Δf1 frequency variation for the area 1, x4 and Δf2 frequency variation for the area 2, x5 and ΔP12 tie-line power flow variation, x6 and ΔPG1 turbine output variation for the plant 1 or the area 1, x7 and Δa1 turbine valve position variation for the plant 1 or the area 1, x8 and ΔPG2 turbine output variation for the plant 2 or the area 2, and x9 and Δa2 turbine valve position variation for the plant 2 or the area 2, and control signals are as follows: u1 and ΔPc1 speed changer variation for the plant 1 or the area 1, and u2 and ΔPc1 speed changer variation for the plant 2 or the area 2.It is significant to note that the astatism is not inherent to the system—it is induced by including the integrals of the area control error (ACE)v1 and v2 into state vector. So the simulation could be performed with the modification of this model where these two variables would be omitted from the state space vector but retained in the output or outputs, resulting in the model having full-order n=7.Here we have chosen another modification. To make the model in the simulation example stabilized, matrixA is changed according to the prescribed degree of stability:
(4.3)A=A+αI9,
where α is chosen to be -10 and I9 denotes unity matrix of dimension 9. Other system matrices remained unchanged. Matrix A which is changed to have the prescribed degree of stability is
(4.4)A=[-1004.50100000-1004.5-1000000-10.050-0.10.1000000-10.050.1000.100032.7-32.7-10000000000-1550000-40000-20000000000-155000-400000-20].However, the system obtained in this way is still unobservable—only four of its modes are observable. From the matrix dimensions in (4.4) it is obvious that n=9 is the order of the system.Table1 shows HSVs of the original and balanced models. On the basis of the HSVs a decision was made that the reduced order could be (and so it was chosen to be) r=2, since the third Hankel singular value is more than 400 (410.4444) times smaller than the second one. For the reduced-order model of order r=2 good quality approximation will be achieved, both in the open- and in closed-loop cases. However, in this manner we would omit state space variables that are important to be retained.Table 1
Hankel singular values.
For the original and balanced modelsFor thefirst ROM obtained through multimodelingFor thesecond ROM obtained through multimodeling0.140060.0990370.0990370.0239210.0169150.0169150.58281·10-40.41211·10-40.41211·10-40.41813·10-50.29566·10-50.29566·10-50.25878·10-80.17367·10-210.13413·10-200.61066·10-90.1615·10-230.19797·10-210.52878·10-101.4047·10-285.6210·10-260.1103·10-100.58426·10-13Multimodeling allows us to reduce the order of the model while keeping all state space variables accounted for, only decoupled into two intersecting subsets.As for the multimodeling, the reduced-order system which consisted of a common core subsystem and subsystem corresponding to modes having a small amount of total energy was calledthe firstwhile the second referred to the one which consisted of a common core subsystem and subsystem corresponding to the remaining modes, having a smaller amount of total energy.According to the suggestions made in [17] the slow or common core subsystem order was chosen to be of order n0=5, and the orders of the two fast subsystems, that is, subsystems corresponding to modes with a small amount of total energy were n1=2, and n2=2, respectively; so the first and the second subsystems were of the same order, n0+n1=n0+n2=7. Forthe first and the second reduced-order models a multimodeling was performed through balanced truncation. In this way two subsystems were obtained having one input, with u1 being input to the first subsystem, as in (3.28)–(3.31), and u2 being input to the second one, as in (3.32)–(3.35).Table1 contains also the HSVs for the first and the second reduced-order models (ROMs) obtained through multimodeling using balanced truncation. It seems that the largest four HSVs for the firstsubsystem and the secondone are identical. However, the difference between them is of the order 10-15 or less.Each of the model approximations as well as the original one was tested on the open-loop case for typical input functions: impulse, step, ramp, and sine. In the time domain all of the approximations produce good and similar performances, except for both thefirstand the second subsystems obtained through multimodeling, that exhibit a considerable degree of difference in behavior when compared to models in the cases of ramp and sine input signals. However, the behavior of the first modeland the second oneis almost identical.The optimal gain matrices were subsequently calculated for the original model, its balanced equivalent, and the four reduced-order models: one obtained throughbalanced truncation, one obtained through balanced residualization, whereas two reduced-order models were comprised of the slow subsystem and one fast subsystem—called the first and the second—consisted of the slow subsystem and the rest of the fast subsystem. Matrices needed for optimization criteria, chosen in accordance with those from [17], and used for the original and balanced nonreduced model and for various reduced-order models are in (4.5):
(4.5)Q=I9,R=20I2,Qtruncated=Qrezidualized=I2,Rtruncated=Rrezidualized=20I2,Qthefirst=Qthesecond=I7,Rthefirst=Rthesecond=20.
Here Ik denotes unity matrix of dimension k.Simulations were carried out using MATLAB. Proportional regulator gain matrices were computed for the original and the balanced models of the full-order and the truncated and residualized reduced-order models of order 2 as well as for the first and the second reduced models of order 7.Corresponding matrices for the first and the second subsystems,Kthefirst and Kthesecond, are close to submatrices of Kbalanced and of the necessary dimensions:(4.6)Kbalanced=[-3.28·10-3-1.561·10-30.34·10-4-0.15·10-4-0.33·10-30.89·10-3-0.32·10-5-0.111·10-3-2.65·10-3-3.28·10-3-1.561·10-30.34·10-4-0.15·10-40.586·10-3-0.89·10-30.32·10-50.111·10-32.65·10-3],Ktruncated=[-3.2778·10-3-1.5597·10-3-3.2778·10-3-1.5597·10-3]Krezidualized=[-3.2785·10-3-1.5635·10-3-3.2785·10-3-1.5635·10-3],Kthefirst=[-3.28·10-3-1.561·10-30.34·10-4-0.15·10-4-0.397·10-3-0.78·10-4-0.247·10-3],Kthesecond=[-3.280·10-3-1.561·10-30.34·10-4-0.15·10-40.586·10-3-0.111·10-32.234·10-3].Finally, a regulator is composed of theKthefirst and Kthesecond as in (4.7)Kbal.IandII=[-3.28·10-3-1.561·10-30.34·10-4-0.15·10-4-0.397·10-3-0.78·10-4-0.247·10-300-3.28·10-3-1.561·10-30.34·10-4-0.15·10-40.586·10-300-0.111·10-32.234·10-3].The closed-loop systems were also tested for various input functions: impulse, step, ramp and sine. Regulators, calculated on reduced order models were subsequently used to form the closed loop on the balanced full-order model.The corresponding Lyapunov equations were solved and the traces of their solutions were compared to derive which of the closed-loop models obtained in this manner is the closest to the optimal case with a regulator having gain matrixKbalanced. The values of the optimal criteria for the balanced model with its optimal P regulator and the values of suboptimal criteria for the balanced model with an optimal regulator for the various reduced-order models, are seemingly identical, as in CriteriaB:
(4.8)CriteriaB=[0.45840.45840.45840.45840.45840.4584].
From the values of CriteriaB it appears that we have reached the desired approximation of the optimization criteria:
(4.9)Japproximated=Jbal.IandII=12Jthefirst+12Jthesecond≈0.4584.
Absolute errors with respect to balanced full-order optimal case are(4.10)CriteriaErrorB=[00.768·10-50.772·10-50.1402·10-40.1030·10-40.469·10-5].However, the difference in the optimal criteria value is in theCriteriaErrorB and for the balanced model with optimal regulator for the truncated model it is 0.768·10-5.For the balanced model with optimal regulator for the residualized the difference in the optimal criteria value is 0.772·10-5. The difference in the optimal criteria value for the balanced model with optimal regulator for the first subsystem is 0.1402·10-4. For the balanced model with optimal regulator for the second subsystem the difference is 0.1030·10-4.For the regulator composed of theKthefirst and Kthesecond gain matrices as described above, the difference between the corresponding suboptimal criteria and the optimal one is the smallest of all tested reduced-order models: 0.469·10-5. This concurs with the quality approximation expected from the results of [17].The corresponding relative errors with respect to the optimal criteria value for the balanced system are(4.11)CriteriaRelativeErrorB=[0%0.0017%0.0017%0.0031%0.0022%0.0010%].
Apparently, a combination of the first and the second optimal gain matrices produces performance of the closed-loop model that is the closest to the balanced full-order model. The relative difference in the optimal criteria value is in the Criteria Relative Error B and, respectively, for the balanced model using regulator optimal for the truncated model it is 0.0017%. As for the balanced model with optimal regulator for the residualized, it is identical to the above value, that is, 0.0017%. For the balanced model using regulator optimal for the first subsystemrelative criteria error is 0.0031%, and finally, for the balanced model with optimal regulator for the second subsystemit is 0.0022%.For the regulator composed of theKthefirst and Kthesecond as described above, the relative difference between the corresponding suboptimal criteria and the optimal one is the smallest of all tested reduced-order models: 0.0010%. This result confirms that the approximation achieved in this way is of a good quality.


## 5. Conclusion

Many order reduction techniques have been developed throughout the past few decades, and this problem will stay under consideration for as long as the engineering practice keeps developing. Here are performed and combined two already established techniques for order reduction by means of multimodeling and balancing. The novel approach has given very good results in the open-loop as well as in the closed-loop model approximation for real engineering purposes. In those cases where physical interpretation of the state space variables requires that not a single variable should be lost, multimodeling provides a way to reduce order and simplify the regulator design. The cases of this kind are the real-world models of power systems and cars, for example. The choice of order reduction technique can be made more appropriate provided that all kinds of information about the original system and the related constraints are incorporated.



---

*Source: 841830-2010-07-18.xml*