# An Augmented Lagrangian Method for the OptimalH∞ Model Order Reduction Problem

Authors: Hongli Yang; Maoxian Zhao
Journal: Mathematical Problems in Engineering
 (2017)

---

## Abstract

This paper treats the computational method of the optimalH∞ model order reduction (MOR) problem of linear time-invariant (LTI) systems. Optimal solution of MOR problem of LTI systems can be obtained by solving the LMIs feasibility coupling with a rank inequality constraint, which makes the solutions much harder to be obtained. In this paper, we show that the rank inequality constraint can be formulated as a linear rank function equality constraint. Properties of the linear rank function are discussed. We present an iterative algorithm based on augmented Lagrangian method by replacing the rank inequality with the linear rank function. Convergence analysis of the algorithm is given, which is distinct to the now available heuristic methods. Numerical experiments for the MOR problems of continuous LTI system illustrate the practicality of our method.


---

## Body

## 1. Introduction

Model order reduction (MOR) problem has received considerable attention since it was put up. In physical or engineering system, the mathematical modeling of the system often results in the high-order controllers, while the simulation and physical implementation of the higher order controllers are more difficult to realize due to the high order. The purpose of the MOR is to obtain a stable lower-order system that approximates the higher one according to some given criterion. Many techniques and results on MOR problems have been reported, and there are various approaches such as the aggregation method [1], balanced truncation approach [2], and optimal Hankel norm approximations methods [3, 4]. It is a difficult work to list all the references since we have many references in this field. The H∞ MOR problem has been studied by many researchers such as [4–13] and references therein. It mainly consists of finding a lower-order model G^ of the original system G such that the H∞ error G-G^ is small. It has received quite a few attentions since the H∞ norm of the difference between two systems is one of the most meaningful measures of approximations. Linear matrix inequality (LMI) is a well-known convex feasibility problem which can be solved by the semidefinite programming (SDP) interior algorithms [14]. LMI is also a powerful tool in control and systems theory and was widely used in output feedback stabilization, reduced-order H∞ synthesis, μ synthesis with constant scaling [15], and so on. In [6], an LMI approach was followed for the H∞ MOR problem, utilizing the Bounded Real Lemma [16], and an iterative two-step algorithm was proposed but with no guaranteed convergence to the global optimum. More recently, [5] studied the deterministic H∞ MOR problem for deterministic case, necessarily and sufficient conditions are derived for γ- suboptimal MOR problems in terms of LMIs, numerical techniques based on alternating projections are proposed to solve the MOR problems, and [7] extends the results to the stochastic MOR problem.The optimal indexγ of H∞ MOR problem is computed by minimizing the model given in [5], and an algorithm to solve the related LMIs feasibility coupling with a rank constraint is also proposed in that paper. The rank constrained LMI problem can be considered as a generalization of LMI problem, but the lack of convexity makes the rank constrained LMI problem much harder to solve. Now available algorithms for the rank constrained LMI problems are largely heuristic in nature and most of them are not convergent. The existing methods for LMI coupling with rank constrained are those based on linearization [17], alternating projections [18, 19], trace minimization methods that solve the problem by solving a related convex problem [20], augmented Lagrangian methods (ALF) [21], sequential semidefinite programming [22], and so on. Aside from [21, 22], these methods did not give the proof of convergence and the challenge still remains to find the numerical schemes with verifiable local/global convergence. More recently, nuclear norm minimizations method is a good relaxation approximation method to this problem [23].The contribution of this paper is as follows: an algorithm is presented using the augmented Lagrangian method by means of the equivalent condition to the rank constraints; ALF method has more advantages compared with the PF method. For PF method, the penalty item becomes more and more possibly ill-conditioned when the penalty parameter become decreasing, which makes the solution more difficult to find. One of the modification methods to solve this problem for PF is the ALF method, which has better performance in the aspect of robustness and convergence rate. So we proposed an ALF method for Problem10. Convergence of the algorithm is also investigated in this paper. Numerical experiments showed that our algorithm is effective.The rest of this paper is structured as follows. Section2 gives an introduction to the H∞ MOR problem and the feasible set of convex LMI constraints coupling with rank constraint. Section 3 presents an equivalent condition of the rank constraints and reformulates the MOR problem as an NLP problem using a cost function which combines H∞ gains index γ and a penalty term accounting for the rank constraint. An augmented Lagrangian method solving the NLP problem is proposed in Section 4, as well as convergence results. Section 5 gives some computational details about the algorithms. Numerical experiments are given in Section 6 for continuous LTI MOR problems. Section 7 concludes this paper.The standard notation used throughout this paper is as follows.Sn denotes the set of all symmetric matrix in Rn, and In denotes the n identity matrix. For a matrix A, its transpose, rank, inverse, and trace are denoted by AT, rank(A), A-1, and trace(A), respectively. If A is a symmetric matrix, A>(≥) 0 means that A is positive (semidefinite) definite. Symbols ∇ and ∇2 mean the gradient and Hessian matrix of a function, σ(A) denotes the maximum singular value of a matrix, ∥·∥ denotes the H∞ norm of a rational transfer function, and A,B denotes the inner product of two matrices.


## 2. Problem Formulation

Consider an asymptotically stablenth-order linear, time-invariant continuous system G with a state space representation (1)x˙=Ax+Bu,y=Cx+Du,where A∈Rn×n, B∈Rn×m, C∈Rp×n, and D∈Rp×m. The optimal H∞ MOR problem is to find a stable n^th order system G^(n^≪n) with a state space representation (2)x^˙=A^x^+B^u^,y^=C^x^+D^u^,where A^∈Rn^×n^, B^∈Rn^×m, C^∈Rp×n^, and D∈Rp×m, such that the H∞ gain index of G-G^∞ with respect to G^ is minimized. The balanced truncation method for MOR problem [2] and the optimal Hankel norm MOR method [3] are the most popular among the now available methods. Methods based on the LMIs are the ones which can reformulate the MOR problem as LMIs by means of the Bounded real lemma. The optimal MOR solution can be obtained by solving a feasible point of LMIs feasible set coupling with a rank constraint. The method is proposed in [6] in terms of time domain response of the error system and [5] extends the result to the H∞ MOR problem in the state space. Reference [7] extends the results to the stochastic case. More recently, basing on LMIs method, [24] presents a sufficient and necessary condition for positivity-preserving H∞ MOR of LTI positive system.Theorem 1 (see [5]).
There exists a stable systemG^ such that the H∞ MOR problem of the nth-order continuous-time system in the form of (1) is solvable for a given γ>0 if and only if there exists positive definite matrix X>0, Y>0 satisfying (3)minγ,X,Yγ(4)AX+XAT+BBT<0,(5)YA+ATY+CTC<0,(6)XγIγIY≥0,rankXγIγIY≤n+n^.The optimal H∞ error and the corresponding feasible solution for X,Y can be obtained by solving the following optimization Problem 2, and all γ- suboptimal of n^th- order models that correspond to a feasible matrix pair (X,Y) are given by (7)D^C^A^B^=G1^+G2^LG3^,where L∈R(p+n^)×(m+n^) is any matrix such that L<1, and (8)G1^=M1-Q12Q22-1M2TM2Q22-1M2T-1,(9)G2^=-Q11+Q12Q22-1Q12T-G1G32G1T1/2,(10)G3^=-M2Q22-1M2T1/2,where(11)M1=000Rx2,M2=0IRxLxT0,Q11=-γ2ICLxRxRxLxTCT0,Q12=CXDRxLxTAT0,Q22=AX+XATBBTIand Rx∈Rn^×n^ is an arbitrary positive definite matrix, Lx∈Rn×n^ is an arbitrary matrix factor such that LxLxT=X-γ2Y-1.Problem 2.
According to Theorem1, the solution of the optimal H∞ MOR problem can be obtained by solving the following LMIs constrained minimization problem coupling with rank constraints.(12)minγ,X,Yγ(13)AX+XAT+BBT<0,(14)YA+ATY+CTC<0,(15)XγIγIY≥0,rankXγIγIY≤n+n^.Assumption 3.
Throughout the paper, we assume that Problem2 has at least one feasible solution.Theorem 4 (see [12]).
The rank of a symmetric positive semidefinite matrixW satisfies rank(W)≤r if and only if there exist matrices V∈Rn×(n-r) with VTV=In-r such that tr(VTWV)=0.Theorem4 is a lemma in [12]; we give a detailed different proof here.Proof.
Given a symmetric positive semidefinite matrixW≥0, rank(W)≤r, r>0, there exists an orthogonal matrix U such that(16)UTWU=λ1⋱λn,where the eigenvalues are ordered λ1≥λ2≥⋯≥λn≥0, since the rank of matrix W is less than or equal to r; there must be at most r eigenvalues that are not 0. The columns of matrix U are the corresponding eigenvectors of the related eigenvalues; take V=(ur+1,…,un)∈Rn×(n-r); the columns of V consist of the eigenvectors corresponding to the n-r smallest eigenvalues of W, and VTV=In-r; this ends the necessity part.
If there exists a matrixV∈Rn×(n-r), with VTV=In-r satisfying tr(VTWV)=0, suppose the eigenvalues of the symmetric matrix VTWV are λ1,λ2,…,λn-r. They are also the eigenvalues of the matrix W since matrix V is full rank which means λ1,λ2,…,λn-r≥0 because W is a symmetric positive semidefinite matrix. We have (17)λ1+λ2+⋯+λn-r=trVTWV=0,from the nonnegativity of the eigenvalues λi, i=1,2,…,n-r, we have λi=0, i=1,2,…,n-r, which means rank(W)≤r; this ends the proof of sufficiency part.For simplification, we note the closure of the semipositive feasible constraints set in Problem2 as Ω={(γ,X,Y)∣γ≥0,∃X,Y>0, such that AX+XAT+BBT≤0, YA+ATY+CTC≤0, -XγIγIY≤0}.Proposition 5.
SetΩ as a convex in variable γ.
The set is the intersection of some polynomial inequalities constrained set. The proof of Proposition5 can be done by the positive definiteness (positive semidefiniteness) of matrix.Proposition 6.
Denoteh(γ;V)=tr(VTWV), where W=XγIγIY; then h(γ) is a linear function in matrix variable X, Y for a given fixed matrix V.Proof.
From Theorem1, matrices X, Y are functions in variable γ; according to the definition of function h(γ) and the inner product of matrix, we have(18)hγ=trVTWV=V,WV,W=XγIγIY.Its obvious that h(γ) is a linear function in X, Y for a given V and W=XγIγIY but keep in mind that X, Y are related to variable γ, and V varies with matrix W. Matrix V can be obtained by the singular decomposition of the specific matrix W.Proposition 7.
For a given matrixW with a fixed rank n^, function h(γ;V) is continuous differentiable in variable γ for fixed X, Y. Suppose γ∗ is the solution of equation h(γ)=0; then when γ→γ∗, h(γ)→0.Proof.
From Theorem4, matrix V is a submatrix of the orthogonal transformation matrix, so VTWV has the same eigenvalues as matrix W. The value of h(γ;V) is the addition of some eigenvalues of matrix W. From the definition of eigenvalues of a matrix and Ostrowski Theorem, the eigenvalues of matrix W are a continuous and differentiable function in the entries of the matrix W (while matrix V is not the case). The proof of this lemma is obvious. Conclusions can be drawn according to the continuity of h(γ) and Proposition 6.Proposition7 suggests that |h(γ;V)|≤ε can be used as the iteration termination criterion, although the rank condition for matrix W in Theorem 1 does not hold.Proposition 8.
If the rank of matrixW is full, function h(γ;V)>0 for any γ≥0 and X, Y satisfies the LMIs; if the rank of matrix W is less than or equal to n+n^(n^≪n), there must exists a γ¯>0 such that h(γ¯;V)=0.Proof.
W
=
X
γ
I
γ
I
Y, it is obvious that rank(W)≥n since W has a block matrix X which is positive definite. If rank(W)≤n+n^, then(19)rankW=rankX00γ2I-XY=rankX+rankγ2I-XY≤n+n^.Since rank(X)=n, we have rank(γ2I-XY)≤n^.
MatrixY=S2, where S is also a symmetric positive definite matrix for the reason that Y is a symmetric positive definite matrix; then we have(20)rankγ2I-XY=rankS-1γ2I-SXSS;there must exist an orthogonal matrix Q such that(21)QTSXSQ=μ1⋱μn,so we have(22)rankγ2I-XY=rankS-1γ2I-SXSS=rankS-1Qγ2-μ1⋱γ2-μnQTS,so there must exist at least n-n^ eigenvalues which are zero; that is, γ¯2=μij, j=1,2,…,n-n^. The property h(0)>0 comes from the positive definite of matrix X, Y. Propositions 5–8 show that the rank constraint in (15) can be replaced by the equality constraint h(γ;V)=0.


## 3. Reformulation of Problem2

By means of the functionh(γ,X,Y) and Theorem 1, Problem 2 can be reformulated as follows.Problem 9.
(23)
m
i
n
γ
,
X
,
Y
γ
A
X
+
X
A
T
+
B
B
T
<
0
,
Y
A
+
A
T
Y
+
C
T
C
<
0
,
X
γ
I
γ
I
Y
≥
0
,
h
γ
;
V
=
0
.The above minimization Problem 9 can be written into the following standard optimization Problem 10, but keep in mind that X, Y is related to γ. We make a notation x=(γ,X,Y) (we omit the dimension of matrix here), which is the decision variable.Problem 10.
(24)
m
i
n
γ
,
X
,
Y
f
x
(25)
h
x
;
V
=
0
,
x
=
γ
,
X
,
Y
∈
Ω
.
Problem10 is an optimization problem with equality constraint and closed feasible set; standard optimization techniques and algorithms can be used to solve this problem. Recently, [12] presented a penalty function (PF) algorithm by means of the rank function given in this paper. Numerical experiments showed that the algorithm is effective while there is no convergence analysis given in that paper. Also there are some disadvantages such as ill conditions for penalty function methods. This motivates us to present a convergent algorithm with good numerical performance for Problem 10 (also Problem 2) by means of the function h(γ;V).


## 4. Augmented Lagrangian Function Algorithm

### 4.1. ALF Algorithm

One of the effective convergent algorithms for Problem10 with equality constraints is the augmented Lagrangian function (ALF) method. ALF method has more advantages compared with PF method. For PF method, the penalty term becomes more and more possibly ill-conditioned when the penalty parameter is decreasing, which makes the solution more difficult to find. One of the modification methods to overcome this problem for PF is the ALF method, which has better performance in the aspect of robustness and convergence rate. So we proposed an ALF method for Problem 10 in this paper. For more detailed information about the ALF methods, refer to [25–28].Problem10 can be considered as a standard equality constrained nonlinear programming problem (NLP). Augmented Lagrangian function method is an effective method to solve an optimization with equality constraints. We reformulate Problem 10 as the form of augmented Lagrangian function optimization Problem 11, where λ is the multiplier and c is the penalty parameter.Problem 11.
(26)
m
i
n
γ
,
X
,
Y
Φ
c
=
f
x
+
λ
h
x
;
V
+
c
2
h
2
x
;
V
(27)
x
=
γ
,
X
,
Y
∈
Ω
,
λ
∈
Λ
.
Suppose thekth iteration solution of Problem 11 is xk=(γk,Xk,Yk). The rationale for the ALF method is based on the fact that when λk is bounded and ck→∞, then the term λh(x;V)+c/2h2(x;V) tends to infinity if h(x;V)≠0 and is equal to 0 if h(x;V)=0. The minimum of Problem 10 can be obtained by solving the unconstrained Problem 11 under some conditions.ALF Algorithm for NLP Problem10Step 1 (initialization).
Initialize the algorithm by determining a feasible of the LMI constraintsX0, Y0, and γ0. Then initialize the Lagrangian multiplier λ0 and penalty parameter c0>0, k≔0,Step 2 (inner optimization).
Fork=0,1,2,…, solving the NLP problem for fixed ck, λk(28)minγ,X,YΦck=fx+λkhx;V+ck2h2x;V(29)x=γ,X,Y∈Ω,λ∈Λ.Using algorithm solving LMI feasibility (such as SeDuMi [29] or YALMIP [29]) to solve the subproblem, an eigenvalue decomposition of matrix W needs to be done at kth iteration, and let (30)xk=γk,Xk,Ykbe the solution of xk=argminΦckf(x)+λkh(x;V)+ck/2h2(x;V). If xk satisfies the termination criterion in Step 4, stop!Step 3 (update penalty and multiplier).
(31)
λ
k
+
1
=
λ
k
+
c
k
h
x
k
+
1
;
V
,
c
k
+
1
=
ρ
c
k
,
i
f
h
x
k
+
1
>
μ
h
x
k
;
V
;
c
k
,
i
f
h
x
k
+
1
;
V
≤
μ
h
x
k
;
Vfor given ρ>1, 0<μ<1. To ensure that λk→λ∗ during the iteration, we use the above update of the multiplier. This method is discussed in [25, 26]. Let k=k+1; go to Step 2.Step 4 (terminating phase).
Due to the linearity of the constraints functionh(x), the terminate criterion of the algorithm can be chosen with (γ,X,Y)∈Ω and |h(x)|≤ε.Proposition 12.
ALF algorithm is well defined.Proof.
From Assumption3, the feasible set of Problem 11 is nonempty; Propositions 5 and 6 showed that there must exist x¯ such that h(x¯)=0; it is obvious that the algorithm is well defined.


### 4.2. The Optimality Conditions for Problem11

We note the Lagrangian part ofΦc(x,λ) as Φ(x,λ); that is, (32)Φcx,λ=Φx,λ+c2h2x;V.Letx∗ be a local minimizer of Problem 10; suppose there exists a Lagrangian multiplier λ∗ such that the first- and second-order optimality conditions for Problem 10 are satisfied [25]. Then there exists a large enough c¯>0 for any c>c¯, at the point (x∗,λ∗); we have(33)Φcx∗,λ∗=Φx∗,λ∗,∇Φcx∗,λ∗=∇Φx∗,λ∗,∇2Φcx∗,λ∗≤∇2Φx∗,λ∗so we can obtain the optimality conditions for Problem 11,(34)∇Φcx∗,λ∗Tx-x∗≥0,for anyx∈Ω,(First order condition)x-x∗T∇2Φcx∗,λ∗x∗-λ∗,for anyx≠x∗≥0.(Second order optimality condition).Then from the Taylor expansion and the optimality conditions, we have(35)Φcx,λ∗≥Φcx∗,λ∗+η2x-x∗2,for any x∈Ω, x-x∗<ε, and certain η, ϵ>0, c>c¯.Proposition 13.
From the above inequality (35), we conclude that there must exist a local minimum of Φc(x,λ) that is close to x∗ for every c>c¯ when λ is close to λ∗.


### 4.3. The Convergence of ALF Algorithm

Theorem 14 (local convergence of the algorithm).
Suppose the feasible set{x∣x=(γ,X,Y)∈Ω} is nonempty. For k=1,2,…, let xk be global minima of the objective function in Problem 10 over the feasible set. If the Lagrangian multiplier sequence xk is bounded, for all k, 0<ck<ck+1, and ck→0, then every convergent point of the sequence xk is a global optimal solution of Problem 10.Proof.
From Propositions5 and 7, f(x), h(x) are continuous in closed set Ω; the proof of the theorem is similar to that of Bertsekas in [25].Remark 15.
The constraint set of Problem11 is a closure of the open set in Problem 2. According to Theorem 14 and Weierstrass Theorem, Problem 11 must have a global minimum. But whether the ALF algorithm will obtain a unique global minimum point or not depends on the concrete problem given. Usually the global point is not unique for the LMIs problem. But the feasible set for Problem 2 is not closed, so the convergent point of Problem 11 may be a local solution to Problem 2.Remark 16.
The convergence rate of the ALF algorithm presented in this paper depends on the optimization method in the subproblem. If the inner algorithm in the subproblem is superlinear convergent, and the ALF algorithm will be superlinear convergent.


## 4.1. ALF Algorithm

One of the effective convergent algorithms for Problem10 with equality constraints is the augmented Lagrangian function (ALF) method. ALF method has more advantages compared with PF method. For PF method, the penalty term becomes more and more possibly ill-conditioned when the penalty parameter is decreasing, which makes the solution more difficult to find. One of the modification methods to overcome this problem for PF is the ALF method, which has better performance in the aspect of robustness and convergence rate. So we proposed an ALF method for Problem 10 in this paper. For more detailed information about the ALF methods, refer to [25–28].Problem10 can be considered as a standard equality constrained nonlinear programming problem (NLP). Augmented Lagrangian function method is an effective method to solve an optimization with equality constraints. We reformulate Problem 10 as the form of augmented Lagrangian function optimization Problem 11, where λ is the multiplier and c is the penalty parameter.Problem 11.
(26)
m
i
n
γ
,
X
,
Y
Φ
c
=
f
x
+
λ
h
x
;
V
+
c
2
h
2
x
;
V
(27)
x
=
γ
,
X
,
Y
∈
Ω
,
λ
∈
Λ
.
Suppose thekth iteration solution of Problem 11 is xk=(γk,Xk,Yk). The rationale for the ALF method is based on the fact that when λk is bounded and ck→∞, then the term λh(x;V)+c/2h2(x;V) tends to infinity if h(x;V)≠0 and is equal to 0 if h(x;V)=0. The minimum of Problem 10 can be obtained by solving the unconstrained Problem 11 under some conditions.ALF Algorithm for NLP Problem10Step 1 (initialization).
Initialize the algorithm by determining a feasible of the LMI constraintsX0, Y0, and γ0. Then initialize the Lagrangian multiplier λ0 and penalty parameter c0>0, k≔0,Step 2 (inner optimization).
Fork=0,1,2,…, solving the NLP problem for fixed ck, λk(28)minγ,X,YΦck=fx+λkhx;V+ck2h2x;V(29)x=γ,X,Y∈Ω,λ∈Λ.Using algorithm solving LMI feasibility (such as SeDuMi [29] or YALMIP [29]) to solve the subproblem, an eigenvalue decomposition of matrix W needs to be done at kth iteration, and let (30)xk=γk,Xk,Ykbe the solution of xk=argminΦckf(x)+λkh(x;V)+ck/2h2(x;V). If xk satisfies the termination criterion in Step 4, stop!Step 3 (update penalty and multiplier).
(31)
λ
k
+
1
=
λ
k
+
c
k
h
x
k
+
1
;
V
,
c
k
+
1
=
ρ
c
k
,
i
f
h
x
k
+
1
>
μ
h
x
k
;
V
;
c
k
,
i
f
h
x
k
+
1
;
V
≤
μ
h
x
k
;
Vfor given ρ>1, 0<μ<1. To ensure that λk→λ∗ during the iteration, we use the above update of the multiplier. This method is discussed in [25, 26]. Let k=k+1; go to Step 2.Step 4 (terminating phase).
Due to the linearity of the constraints functionh(x), the terminate criterion of the algorithm can be chosen with (γ,X,Y)∈Ω and |h(x)|≤ε.Proposition 12.
ALF algorithm is well defined.Proof.
From Assumption3, the feasible set of Problem 11 is nonempty; Propositions 5 and 6 showed that there must exist x¯ such that h(x¯)=0; it is obvious that the algorithm is well defined.


## 4.2. The Optimality Conditions for Problem11

We note the Lagrangian part ofΦc(x,λ) as Φ(x,λ); that is, (32)Φcx,λ=Φx,λ+c2h2x;V.Letx∗ be a local minimizer of Problem 10; suppose there exists a Lagrangian multiplier λ∗ such that the first- and second-order optimality conditions for Problem 10 are satisfied [25]. Then there exists a large enough c¯>0 for any c>c¯, at the point (x∗,λ∗); we have(33)Φcx∗,λ∗=Φx∗,λ∗,∇Φcx∗,λ∗=∇Φx∗,λ∗,∇2Φcx∗,λ∗≤∇2Φx∗,λ∗so we can obtain the optimality conditions for Problem 11,(34)∇Φcx∗,λ∗Tx-x∗≥0,for anyx∈Ω,(First order condition)x-x∗T∇2Φcx∗,λ∗x∗-λ∗,for anyx≠x∗≥0.(Second order optimality condition).Then from the Taylor expansion and the optimality conditions, we have(35)Φcx,λ∗≥Φcx∗,λ∗+η2x-x∗2,for any x∈Ω, x-x∗<ε, and certain η, ϵ>0, c>c¯.Proposition 13.
From the above inequality (35), we conclude that there must exist a local minimum of Φc(x,λ) that is close to x∗ for every c>c¯ when λ is close to λ∗.


## 4.3. The Convergence of ALF Algorithm

Theorem 14 (local convergence of the algorithm).
Suppose the feasible set{x∣x=(γ,X,Y)∈Ω} is nonempty. For k=1,2,…, let xk be global minima of the objective function in Problem 10 over the feasible set. If the Lagrangian multiplier sequence xk is bounded, for all k, 0<ck<ck+1, and ck→0, then every convergent point of the sequence xk is a global optimal solution of Problem 10.Proof.
From Propositions5 and 7, f(x), h(x) are continuous in closed set Ω; the proof of the theorem is similar to that of Bertsekas in [25].Remark 15.
The constraint set of Problem11 is a closure of the open set in Problem 2. According to Theorem 14 and Weierstrass Theorem, Problem 11 must have a global minimum. But whether the ALF algorithm will obtain a unique global minimum point or not depends on the concrete problem given. Usually the global point is not unique for the LMIs problem. But the feasible set for Problem 2 is not closed, so the convergent point of Problem 11 may be a local solution to Problem 2.Remark 16.
The convergence rate of the ALF algorithm presented in this paper depends on the optimization method in the subproblem. If the inner algorithm in the subproblem is superlinear convergent, and the ALF algorithm will be superlinear convergent.


## 5. Computational Issues of the ALF Algorithm and Choice of Parameter

### 5.1. Choice of Initial Point

The initial point of the algorithm can be found by solving the following LMIs feasibility problem using YALMIP [30]:(36)mint∣AX+XAT+BBT<tI,YA+ATY+CTC<tI,-XγIγIY<tI;then determine X0, Y0, γ0 such that h(x) is as close as possible to 0. The initial value γ0 can be chosen as the square root of minimum eigenvalue of matrix SXS, which shows in Theorem 1 that γ→γ¯ as k→∞.


### 5.2. The Inner Optimization for the Subproblem in Step2

About the inner optimization for the subproblem in Step2, we use the interior point algorithm for LMI feasibility such as SeDuMi in [29] and YALMIP in [30], respectively, in our numerical experiments. At the kth iteration, we replace h(xk;V) with h~(xk;V)=tr(Vk-1TW(x)Vk-1). Since h(x;V) is continuous, when x→x∗, we have (37)hx;V-h~x;V=trVTWxV-trVk-1TWxVk-1=trVTWxV-Vk-1TWxVk-1⟶0.


### 5.3. The Choice of the Penalty Parameter and Lagragian Multiplier

#### 5.3.1. Penalty Parameterc0

The initial value ofc0 should not be too large, and ck should not be increased too fast to a point where the subproblem (27) and (28) becomes ill-conditioned, also ck should not be increased too slowly in that the first-order update of the Lagrange multiplier has a poor convergence rate otherwise.


#### 5.3.2. Lagrangian Multiplierλ0

Any prior knowledge should be exploited to selectλ0 close to λ∗, but this is generally difficult. Many trials need to be done in our experiments.


## 5.1. Choice of Initial Point

The initial point of the algorithm can be found by solving the following LMIs feasibility problem using YALMIP [30]:(36)mint∣AX+XAT+BBT<tI,YA+ATY+CTC<tI,-XγIγIY<tI;then determine X0, Y0, γ0 such that h(x) is as close as possible to 0. The initial value γ0 can be chosen as the square root of minimum eigenvalue of matrix SXS, which shows in Theorem 1 that γ→γ¯ as k→∞.


## 5.2. The Inner Optimization for the Subproblem in Step2

About the inner optimization for the subproblem in Step2, we use the interior point algorithm for LMI feasibility such as SeDuMi in [29] and YALMIP in [30], respectively, in our numerical experiments. At the kth iteration, we replace h(xk;V) with h~(xk;V)=tr(Vk-1TW(x)Vk-1). Since h(x;V) is continuous, when x→x∗, we have (37)hx;V-h~x;V=trVTWxV-trVk-1TWxVk-1=trVTWxV-Vk-1TWxVk-1⟶0.


## 5.3. The Choice of the Penalty Parameter and Lagragian Multiplier

### 5.3.1. Penalty Parameterc0

The initial value ofc0 should not be too large, and ck should not be increased too fast to a point where the subproblem (27) and (28) becomes ill-conditioned, also ck should not be increased too slowly in that the first-order update of the Lagrange multiplier has a poor convergence rate otherwise.


### 5.3.2. Lagrangian Multiplierλ0

Any prior knowledge should be exploited to selectλ0 close to λ∗, but this is generally difficult. Many trials need to be done in our experiments.


## 5.3.1. Penalty Parameterc0

The initial value ofc0 should not be too large, and ck should not be increased too fast to a point where the subproblem (27) and (28) becomes ill-conditioned, also ck should not be increased too slowly in that the first-order update of the Lagrange multiplier has a poor convergence rate otherwise.


## 5.3.2. Lagrangian Multiplierλ0

Any prior knowledge should be exploited to selectλ0 close to λ∗, but this is generally difficult. Many trials need to be done in our experiments.


## 6. Numerical Experiments

We provide two examples for the applicability of the above method for determining the suboptimal model errors and the corresponding reduced state space. We use the benchmark model reduction examples: Wilson’s example and Boiler system. Both examples are continuous-time LTI systemH∞ model reduction problems.Example 1.
Consider the fourth-order system reduces to a second system, Wilson example in [31].(38)x˙t=000-150100-245010-113001-19xt+4100ut,yt=0001xt.Example 2.
Consider the ninth-order system reduce to a third-order system, Boiler System in [32].(39)x˙t=-0.910000000000-4.449000000000-10.262571.4790000000-571.479-10.262000000000-10.987000000000-15.21411.6220000000-11.622-15.214000000000-89.874000000000-502.665xt+-4.336-3.69110.141-1.61216.629-242.476-14.26113.67282.187ut,yt=-0.422-0.736-0.04160.232-0.816-0.7150.546-0.235-0.080xt.
Tables1 and 2 list the optimal indexes obtained by our algorithm. Figures 1 and 2 indicate the error and the convergence by the ALF algorithm. The computation parameters used in the ALF algorithms were ϵ=1e-12, μ=0.85, ρ=1.1, c0=1000, and λ0=1000. Our experience showed that this algorithm works quite well for both problems compared with the now available algorithms [31, 32] under the frame of feasibility coupled with a rank constraint. The iterations were both less than 100 times when the termination were satisfied; this showed that our algorithm has good convergence performance and also the error is less than the results obtained by other methods, especially for Boiler System; our algorithm performed quite better than that in [32]. The suboptimal errors can be obtained as well as the corresponding positive matrices X, Y by means of our algorithm. The reduced state space can be evaluated according to (6)–(8) (Theorem 1 in [5]), and the obtained reduced system is stable.Table 1
Wilson’s example.
Method
γ
Wilson’s model reduction
2.525e-4
H
∞ model reduction
2.729e-4
ALF model reduction
1.638e-4Table 2
Boiler system.
Method
γ
H
∞ model reduction
1.127e-1
ALF model reduction
0.043e-1Figure 1
Gamma convergence.Figure 2
Gamma convergence.


## 7. Conclusion

We presented an ALF algorithm for optimalH∞ MOR problem of the LTI system by means of an augmented Lagrangian method. First, we give a rank function which can reformulate the rank constraint of the optimal H∞ MOR problem as a linear function equality constraint; then we reformulate the original problem as an optimization problem with LMIs constraints and linear equality constraints; an ALF algorithm is presented in the following. Compared with the now available algorithms, the goodness of our algorithm is the better convergence, which is opposite to the heuristic algorithm; our algorithm can also overcome the ill condition of PF algorithm in the process iteration. The algorithm is applicable to both continuous system and discrete system. Numerical experiments showed that the algorithm is effective for deterministic H∞ MOR problem. MOR is an interesting and practical topic; not only is it applied in physics or engineering systems, but also it can be used in some models in biological systems [33–36].



---

*Source: 8230978-2017-11-09.xml*