# On the Solution of the Eigenvalue Assignment Problem for Discrete-Time Systems

Authors: El-Sayed M. E. Mostafa; Abdallah W. Aboutahoun; Fatma F. S. Omar
Journal: Journal of Applied Mathematics
 (2017)

---

## Abstract

The output feedback eigenvalue assignment problem for discrete-time systems is considered. The problem is formulated first as an unconstrained minimization problem, where a three-term nonlinear conjugate gradient method is proposed to find a local solution. In addition, a cut to the objective function is included, yielding an inequality constrained minimization problem, where a logarithmic barrier method is proposed for finding the local solution. The conjugate gradient method is further extended to tackle the eigenvalue assignment problem for the two cases of decentralized control systems and control systems with time delay. The performance of the methods is illustrated through various test examples.


---

## Body

## 1. Introduction

In this work, we consider the following unconstrained minimization problem:(1)minK∈Rp×rfK=ρAK,where f:Rp×r→R+ is generally semismooth and nonconvex and ρ(·) is the spectral radius of the matrix function A(K) that will be defined later on.It is desirable to have a local solutionK∗ of (1) such that ρ(A(K∗))<1. Therefore, we include such a constraint as a cut to the objective function implying the following inequality constrained minimization problem: (2)minfK=ρAKs.t.gK=ρAK-1+τ≤0,where τ∈(0,1) is a given constant.It is well known that, for any matrix norm, it holds thatρ(M)≤M, where M is any square matrix. Then, we can replace the eigenvalue constraint of problem (2) by a regular inequality constraint which yields the following relaxed minimization problem: (3)minKfK=ρAKs.t.g^K=AK-1+τ≤0,where τ is as defined above. This problem can be tackled easily by any constrained optimization solver. However, there is no guarantee in general that a feasible solution exists for this problem.The three problems (1)–(3) concern the well-known eigenvalue assignment problem (EAP) for discrete-time systems. These problems are generally nonconvex and semismooth optimization problems. Over the past decades, a considerable amount of attention has been given to the EAP where a lot of research can be found in systems and control literature in particular for continuous-time systems (see, e.g., [1–3] and the references therein). In the framework of linear discrete-time systems, a set of eigenvalue assignment algorithms have been developed (see, e.g., [2–11] and the references therein).A related problem in output feedback control design is the linear-quadratic control problem in which the goal is to design an output feedback gain matrixK that minimizes a certain performance cost function while all eigenvalues of the closed-loop system matrix A(K) must be within the unit circle (see, e.g., [12]). Such a controller K can be calculated to this problem in case of continuous-time systems by available public software packages (e.g., HIFOO) [13]. However, for discrete-time systems, to the best of our knowledge, public software has not yet been developed.In this work, we focus on the two problems (1) and (2) where the attempt is to minimize the spectral radius of the nonsymmetric real matrix A(K). In this regard, we apply a three-term nonlinear conjugate gradient (CG) method [14] to find a local solution to the minimization problem (1) which is attempting to stabilize the associated control system; see the next section. In addition, a logarithmic barrier interior-point method is proposed to tackle the constrained minimization problem (2) for the same purpose.Nonlinear conjugate gradient methods are widely studied and comprise a class of unconstrained optimization algorithms which are characterized by low memory requirements and strong global convergence properties (see the survey [15] and later references [3, 14]). We focus on a three-term nonlinear CG method which has a nice performance (see [14]). The CG methods are descent direction methods, which means that, starting from a given point K0∈Rp×r, these methods generate a sequence {Kk}⊂Rp×r according to the following relation: (4)Kk+1=Kk+αkΔKk,where the step size αk>0 satisfies the line search rule and ΔKk is a descent direction for f at Kk. The update of the new search direction varies from one CG method to another.The logarithmic barrier interior-point method is one of the standard methods for solving constrained optimization problems (see, e.g., [16]). This method is employed to tackle the constrained problem (2) which converts it into a sequence of unconstrained minimization problems.This article is organized as follows. In the next section, we state the formulation of the eigenvalue assignment problem and introduce some basic concepts which are needed in the subsequent analysis. In Section3, we evaluate the required derivatives of the objective function. In Section 4, we introduce the proposed three-term CG method that finds a local solution of the unconstrained minimization problem. In Section 5, we extend the CG method to tackle the output feedback EAP problem for decentralized control systems. In Section 6, we reformulate the discrete-time control system with time delay as an augmented system without any delay so that the EAP can be tackled by the proposed CG method. In Section 7, we introduce a logarithmic barrier method for finding a local solution to problem (2). In Section 8, we test the proposed methods on different test examples from the literature. Then, we end with a conclusion.Notations. For vectors, the symbol · is the 2-norm, while for matrices · denotes the Frobenius norm defined by M=M,M, where ·,· is the inner product given by M1,M2=TrM2∗TM1 and Tr· is the trace operator. The eigenvalues of a matrix M∈Rn×n are denoted by λiM,  i=1,…,n. The Greek letter ρM denotes the spectral radius of a square matrix M and ΛM≔diag(λ1,…,λn) is a diagonal matrix with eigenvalues on its main diagonal. Moreover, Re(·) and Im(·) denote the real and imaginary parts of a complex number, respectively. Sometimes and for the sake of simplicity, we skip the arguments of the considered functions; for example, we use fk to denote f(Kk).


## 2. Problem Formulation and Preliminary

The static output feedback eigenvalue assignment problem for discrete-time systems can be stated as follows (see, e.g., the abovementioned citations and the references therein). Consider the linear time-invariant discrete-time system with the following state space realization:(5)xk+1=Axk+Buk,yk=Cxk,k=0,1,2,…,where x(k)∈Rn, u(k)∈Rp, and y(k)∈Rr are the state, the control input, and the measured output vectors, respectively. Moreover, A∈Rn×n, B∈Rn×p, and C∈Rr×n are given constant matrices. Such a system is often closed by the control law u(k)=Ky(k) which yields (6)xk+1=A+BKCxk=AKxk,k=0,1,2,…,where A(K)≔A+BKC and K∈Rp×r is the output feedback gain matrix which represents the unknown.The following definitions are needed for later use.Definition 1.
The spectral radius of a matrixA∈Cn×n with eigenvalues λ1,…,λn is defined as (7)ρA=maxλi:i∈1,2,…,n.Definition 2.
The discrete-time control system (6) is asymptotically stable (i.e., x(k)→0 as k→∞ for any initial x(0)) if and only if ρ(A(K))<1.The eigenvalue assignment problem is to design an output feedback gain matrixK providing a closed-loop system in a satisfactory stage by shifting controllable eigenvalues to desirable locations in the complex plane. In particular, the EAP requires the spectral radius of the closed-loop matrix A(K) to be strictly within the unit circle in the complex plane.A necessary condition for the EAP by constant output feedback is given in [9]. Fu [7] has also shown that the EAP via static output feedback is NP-hard. Systems with a symmetric state space realization, namely, A=AT, C=BT, occur in different applications such as RC networks [10]. The symmetric EAP might be stated as follows: find a matrix K such that (8)ρA+BKBT<1.It is well known that the eigenvalues of a real symmetric matrix are not everywhere differentiable. A classical theorem (see, e.g., [17]) states that each eigenvalue of a symmetric matrix is the difference of two convex functions, which implies that the eigenvalues are semismooth functions. This fact allows us to use the theory of semismoothness to establish convergence results for the proposed methods.Definition 3.
A functionalψ:Rp×r→R is said to be locally Lipschitz continuous at K∈Rp×r with a constant L>0 if and only if there exists a number σ>0 such that (9)ψK1-ψK2≤LK1-K2,∀K1,K2∈BσK,where Bσ(K) is an open ball with center at K and radius σ.Letψ:Rp×r→R be a locally Lipschitz continuous function. Radmacher’s theorem (see, e.g., [18]) implies that such a mapping ψ is differentiable almost everywhere. Let Dψ be the set of all K at which ψ is differentiable and let ∇ψT be its Jacobian whenever it exists. Let D⊆Rp×r be the set of all K at which f is differentiable. Such a set is open. Therefore, it is convenient to replace D by the following level set: (10)Ω≔K∈D:fK≤c,where c>0 is a given constant. This set is assumed to be bounded.


## 3. Derivative of the Objective Function

The eigenvalues of the matrixA(K) of (1) are not in general differentiable at a point K where A(K) has repeated eigenvalues. Therefore, let us consider the following assumption on eigenvalues of A(K).Assumption 4.
Assume thatA(K) has no multiple eigenvalues for all K∈Rp×r.LetA(K) be diagonalizable and let U(K)∈Cn×n and V(K)∈Cn×n be a couple of matrices whose columns are the left and right eigenvectors of A(K); that is, U(K), V(K), and Λ(A(K)) satisfy (11)UKTAK=ΛUKT,AKVK=VKΛ.The columns ofU(K) and V(K) satisfy (12)Ui∗Vi≠0,i=1,…,n.The following lemma provides the first-order derivatives of the objective functionf(K) of the minimization problem (1) required by the CG methods.Lemma 5.
Suppose thatA(K) satisfies Assumption 4 and f(K) is differentiable at K∈Rp×r. Let λ^=Re(λ^)+iIm⁡(λ^) be the largest in magnitude eigenvalue of A(K). Then, the entries of the gradient of the objective function f are given by (13)∂fK∂Kkl=Reλ^λ^u1TBEklCv1+u2TBEklCv2+Imλ^λ^u1TBEklCv2-u2TBEklCv1,k=1,…,p,l=1,…,r,where (u1+iu2) and (v1+iv2) are the left and right eigenvectors associated with λ^ and Ekl∈Rp×r is a matrix with zero entries except at the (k,l) position where its value is one.Proof (see [19, Lemma 6.3.10]).
The two eigenvectors are normalized such that(14)u1-iu2Tv1+iv2=1,which by differentiation yields (15)u1′-iu2′Tv1+iv2+u1-iu2Tv1′+iv2′=0,where the dash denotes the first-order derivative with respect to the entries of K. The eigenvalue λ^ also satisfies that (16)λ^=u1-iu2TAKv1+iv2,which by differentiation gives (17)λ^′=u1-iu2TAK′v1+iv2=u1TAK′v1+u2TAK′v2+iu1TAK′v2-u2TAK′v1.Suppose that λ^=Re(λ^)+iIm⁡(λ^) is the largest in magnitude eigenvalue of A(K). Then, (18)ρAK=λ^=Reλ^2+Imλ^2,and by the chain rule we have (19)∂ρAK∂Kkl=∂ρAK∂Reλ^∂Reλ^∂Kkl+∂ρAK∂Imλ^∂Imλ^∂Kkl=Reλ^λ^u1T∂AK∂Kklv1+u2T∂AK∂Kklv2+Imλ^λ^u1T∂AK∂Kklv2-u2T∂AK∂Kklv1,k=1,…,p,l=1,…,r.Since ∂(A(K))/∂Kkl=BEklC, then the result follows.


## 4. Three-Term CG Method for the EAP

In this section, we analyze and study a three-term nonlinear CG method for computing the local solution of the minimization problem (1) (see [14]). For a given starting K0∈Rp×r, the CG method generates a sequence of iterates according to the relation (20)Kk+1=Kk+αkΔKk,where ΔKk∈Rp×r is supposed to be a descent direction for f at Kk and αk>0 is the step size. The new search direction for the CG method is given by the following relation: (21)ΔKk+1=-∇fk+1-δkSk-ηkYk,ΔK0=-∇f0,where ∇fk+1 is the gradient of f at Kk+1, Sk=Kk+1-Kk, Yk=∇fk+1-∇fk, and (22)δk=1+2Yk2TrYkTSkTrSkT∇fk+1TrYkTSk-TrYkT∇fk+1TrYkTSk,ηk=TrSkT∇fk+1TrYkTSk.In order to globalize this CG method, we recall Wolfe conditions (see, e.g., [16]) to update a suitable step size αk for the new iterate (20): (23)fKk+αkΔKk≤fKk+γαkTr∇fkTΔKk,(24)Tr∇fk+1TΔKk≥γ^Tr∇fkTΔKk,where 0<γ<γ^<1. The strong Wolfe conditions replace (24) by the following condition: (25)Tr∇fk+1TΔKk≤γ^Tr∇fkTΔKk.The following theorem shows that the search directionΔKk is a descent direction to the objective function.Theorem 6.
LetKk∈D be generated by (20) and let the step size αk satisfy Wolfe conditions (23)-(24). Then, ΔKk+1 evaluated by (21)-(22) is a descent direction to f for all k=0,1,2,….Proof (see also [14, Proposition 3.1]).
From Wolfe’s second condition (24), we obtain the curvature condition Tr(YkTSk)>0. Moreover, from (21)-(22), it follows that (26)Tr∇fk+1TΔKk+1=-∇fk+12-1+2Yk2TrYkTSkTrSkT∇fk+12TrYkTSk≤0.The three-term nonlinear CG algorithm is stated in the following lines.Algorithm 7 (three-term CG method for the output feedback EAP).
(0)
LetK0∈Rp×r, 0<γ<γ^<1, 0<τ<1 be given constants and ϵ∈(0,1) be the tolerance. Moreover, let A,B,C be given constant matrices. Choose α0>0, compute ∇f0, and set ΔK0=-∇f0. If ∇f0≤ϵ or f0<1-τ, stop; otherwise, set k←0 and go to the next step.While
∇
f
k
>
ϵ or fk≥1-τ, do(1)
Computeαk>0 that satisfies Wolfe conditions (23)-(24); set Kk+1=Kk+αkΔKk and then calculate the gradient ∇fk+1.(2)
If∇fk+1≤ϵ or fk+1<1-τ, stop; otherwise, go to the next step.(3)
Calculate a new search directionΔKk+1 by (21)-(22).(4)
Setk←k+1 and repeat.End (do)Remark 8.
For the considered linear control system, one of the major tasks is to achieve a stabilizing output feedback controller, that is. to calculateK such that ρ(A(K)) is strictly less than one. However, there is no relationship between achieving such a controller and finding K, a local minimum of problem (1). Therefore, it is reasonable to stop the CG method as soon as a local solution is achieved or a stabilizing K is reached with a sufficient stability margin.To prove global convergence of the CG algorithm, we assume that∇fk≠0 for all k; otherwise, a stationary point is found.Assumption 9.
The following assumptions are assumed to hold:(a)
The level setΩ as defined in (10) is bounded; that is, there exists a constant β>0 such that, for all K∈Ω, K<β.(b)
In some neighborhoodN of the level set (10), the gradient ∇f(K) of the objective function f is Lipschitz continuous.From Assumption9(a), one has Sk≤Kk+1+Kk≤2β and from Assumption 9(b) we deduce that there exists a constant Γ>0 such that (27)∇fK≤Γ∀K∈Ω.The following lemma provides a lower bound to the step sizeαk; see [14, Proposition 5.1].Lemma 10.
Let{Kk} be generated by Algorithm 7 and suppose that ΔKk is a descent direction for fk. In addition, let Assumption 9 hold. Then, (28)αk≥1-γ^Tr∇fkTΔKkLΔKk2.Proof.
By subtractingTr(∇fkTΔKk) from both sides of (24) and using Lipschitz condition, (29)γ^-1Tr∇fkTΔKk≤Tr∇fk+1-∇fkΔKk=TrYkTΔKk≤TrYkTΔKk≤YkΔKk≤αkLΔKk.Since ΔKk is a descent direction for f and γ^ is less than one, then (28) follows.The result of the following lemma is used in proving the main global convergence theorem.Lemma 11.
Let{Kk} be generated by Algorithm 7 and assume that ΔKk is a descent direction. Furthermore, let Assumption 9 hold. Then, (30)∑k=1∞Tr∇fkTΔKk2ΔKk2<+∞.Proof.
From Wolfe condition (23) and Lemma 10, one has (31)fKk-fKk+αkΔKk≥-γαkTr∇fkTΔKk≥-γ1-γ^Tr∇fkTΔKk2LΔKk2.Then, Assumption 9 implies condition (30).Next, we have the following result for the three-term CG method (see, e.g., [14, Proposition 5.3]).Lemma 12.
Let Assumption9 hold. Consider {Kk} to be generated by Algorithm 7 where ΔKk is a descent direction for f and αk satisfies Wolfe conditions (23) and (24). If (32)∑k>11ΔKk2=∞,then(33)liminfk→∞∇fk=0.Proof (see also [20, Lemma 2.2]).
Suppose (33) is not true. Then, there exists a constant ε>0 such that (34)∇fk≥ε∀k.Therefore, from (32) and Lemma 11, we have (35)∑k=0∞ε4ΔKk2≤∑k=0∞∇fk4ΔKk2≤∑k=0∞Tr∇fkTΔKk2ΔKk2=∞.This contradicts Lemma 12 which completes the proof.Under condition (30), the following global convergence result is obtained.Theorem 13.
LetΔKk be generated by Algorithm 7. Assume that Assumption 9 holds. Assume further that there exists a constant Γ1>0 such that Tr(YkTSk)≥Γ1 for any k≥1. Then, it holds that (36)liminfk→∞∇fk=0.Proof (see also [14, Theorem 5.2]).
SinceTr(∇fkTSk)<0 for all k, then (37)TrSkT∇fk+1=TrYkTSk+Tr∇fkTSk<TrYkTSk.By direct computation, we have (38)ΔKk+1=-∇fk+1+TrYkT∇fk+1TrYkTSk-1+2Yk2TrYkTSkTrSkT∇fk+1TrYkTSkSk-TrSkT∇fk+1TrYkTSkYk.Consequently, one has (39)ΔKk+1≤∇fk+1+TrYkT∇fk+1TrYkTSkSk+1+2Yk2TrYkTSkSk+Yk<Γ+1+2LΓβ+8L2β2Γ12β+2Lβ.Then, from Lemma 12, (36) follows.


## 5. The EAP for Decentralized Control Systems

Consider the linear time-invariant decentralized control system withν control stations: (40)xk+1=Axk+∑i=1νBiuik,yik=Cixk,i=0,1,2,…,ν,where x(k)∈Rn, ui(k)∈Rpi, and yi(k)∈Rri are the state, the control input, and the measured output vectors, respectively. A∈Rn×n, Bi∈Rn×pi, and Ci∈Rri×n are given constant matrices, i=1,…,ν.The output feedback EAP for the decentralized system (40) is to find output feedback gain matrices Ki∈Rpi×ri that place, by using the control law (41)uik=Kiyik,i=1,…,ν,the eigenvalues of the closed-loop system matrix A(K1,…,Kν)=A+∑i=1νBiKiCi to strictly lie within the unit disk.By introducing the following augmented matrices, it is straightforward to rewrite the decentralized control system (40) in the original structure (5): (42)B=B1⋯Bν,C=C1T⋯CνTT.The corresponding closed-loop system matrix is A(K)=A+BKC, where the output feedback gain matrix K is given by the block-diagonal matrix (43)K=diagK1,…,Kν,and the corresponding unconstrained minimization problem takes the form (44)minK1,…,KνfK1,…,Kν=ρAK.The three-term CG method of the last section can be applied to tackle the minimization problem (44). Let D^⊆Rp1×r1×⋯×Rpν×rν be the set of all (45)K=diagK1,…,Kν∈Rp1×r1×⋯×Rpν×rνat which f is differentiable. For a given starting K0, the three-term CG method generates a sequence of the form (46)Kk+1=Kk+αkΔKk,k=0,1,2,…,where αk>0 is the step size that must satisfy Wolfe conditions (23)-(24), while the search direction ΔKk=diag(ΔKk1,…,ΔKkν) is a descent direction for fk at Kk. The new search direction ΔKk+1 is updated by using (21)-(22).Algorithm7 is restated in the following lines to tackle the output feedback EAP problem for decentralized control systems.Algorithm 14 (three-term CG method for the decentralized output feedback EAP).
(0)
LetK0∈Rp1×r1×⋯×Rpν×rν, ϵ∈(0,1), 0<τ<1, and 0<γ<γ^<1 be given constants. Moreover, let A,B1,…,Bν,C1,…,Cν be given constant matrices. Choose α0>0. Compute ∇f0 and set ΔK0=-∇f0. If ∇f0≤ϵ or f0<1-τ, stop; otherwise, set k≔0 and go to the next step.While
∇
f
k
>
ϵ or fk<1-τ, do(1)
Calculate a step sizeαk>0 that satisfies Wolfe conditions (23)-(24); set Kk+1=Kk+αkΔKk and then calculate the gradient ∇fk+1.(2)
If∇fk+1≤ϵ or fk+1<1-τ, stop; otherwise, go to the next step.(3)
Calculate the new search directionΔKk+1 according to (21)-(22).(4)
Setk←k+1 and repeat.End (do)Remark 15.
The way in which the CG method is designed allows us to maintain the block-diagonal structure of the unknown matrixK efficiently without skipping any information of the data matrices of the control system.


## 6. The EAP for Time-Delay Systems

Consider the following linear discrete-time and time-delay system (see, e.g., [21, 22]): (47)x^k+1=∑i=0dAix^k-i+B^uk,yk=∑i=0dCix^k-i,k=0,1,2,…,where Ai∈Rn×n, B^∈Rn×p, and Ci∈Rr×n are given constant matrices and d∈N is the time delay in the state vector.The above time-delay system has a specific feature that it can be converted into an augmented linear system without any delay (see, e.g., [22]). This is achieved by introducing the augmented state vector (48)xTk=x^Tk⋯x^Tk-d.The time-delay system (47) is equivalently rewritten as (49)xk+1=Axk+Buk,yk=Cxk,where (50)A=A0A1⋯Ad-1AdIn0⋯000In⋱00⋮⋱⋱⋮⋮00⋯In0,B=B^00⋮0,C=C0,C1,…,Cd-1,Cd.By using the control lawu(k)=Ky(k), we close system (49) as (51)xk+1=A+BKCxk≕AKxk,where A(K)=A+BKC. For such a system, one can apply Algorithm 7 to compute an output feedback gain matrix K∗, a local solution to the minimization problem (1). It is known that the eigenvalues of the time-delay system coincide with the eigenvalues of the augmented linear system (49).We emphasize that the same approach can be considered in case of control systems with input delay.


## 7. Logarithmic Barrier Method

The attempt of this section is to improve formulation (1), where we aim to achieve a stabilizing K together with minimizing the spectral radius of the closed-loop system matrix. Therefore, let us consider the inequality constrained minimization problem (2). By following the idea of the logarithmic barrier method, we consider the following unconstrained minimization problem: (52)minK∈Rp×rϕμK=fK-μlog⁡-gK,where μ>0 is the barrier parameter and the functions f and g are as defined in (2). According to the theory of barrier methods, the minimizer of ϕμ approaches a local solution of the original problem (2) as μ→0 under certain conditions.Let us define the following strict feasible region:(53)F≔K∈Rp×r:gK≤0.We assume that F≠∅.The Lagrangian functions for problem (2) are defined as (54)lK,υ=fK+υgK,where f and g are as defined in (2) and υ≥0 is the associated Lagrange multiplier.First-order derivative ofϕμ with respect to the entries of K is obtained in the following lemma.Lemma 16.
Suppose thatA(K) satisfies Assumption 4 and f(K) is differentiable at K∈Rp×r. Let λ^=Re(λ^)+iIm⁡(λ^) be the largest in magnitude eigenvalue of A(K). Then, the first derivative of ϕμ is given by (55)∂ϕμK∂Kkl=1+μgKReλ^λ^u1TBEklCv1+u2TBEklCv2+1+μgKImλ^λ^u1TBEklCv2-u2TBEklCv1,k=1,…,p,l=1,…,r,where (u1+iu2) and (v1+iv2) are the left and right eigenvectors associated with λ^ and Ekl∈Rp×r is a matrix with zero entries except at the (k,l) position where it has a value of one.Proof.
By differentiatingϕμ with respect to the entries of K and utilizing the derivative of the spectral radius ρ(A(K)) obtained in Lemma 5, the result follows.As can be seen in (55), all terms of the first derivative of ϕμ depend on μ, which leads to an ill-conditioned Hessian matrix of the barrier function ϕμ at the solution. Therefore, second-order methods are not recommended to compute a local solution of (52), but rather first-order methods such as nonlinear CG methods of Section 4 are recommended.The gradient ofl takes the form (56)∇KlK,υ=∇fK+υ∇gK=1+υ∇ρAK.From (54) and (55), we see that if υ satisfies (57)υ∗=μ∗gK∗,then the solution obtained by the proposed log-barrier method satisfies the stationary requirement of the Karush–Kuhn–Tucker conditions.The logarithmic barrier interior-point method is stated in the following lines.Algorithm 17 (logarithmic barrier method for the output feedback EAP).
(0)
Choose a starting barrier parameterμ0>0, outer and inner-loop tolerances ϵ,ϵin>0, a starting feasible point K0s∈Rp×r. Let a,b,τ∈(0,1) be given constants.For
j
=
0,1
,
2
,
…
,
d
o
(1)
Find an approximate local minimizerKj of ϕμj(K) starting from Kjs; terminate if ∇ϕμj≤ϵjin.
(2)
If the final stopping test∇ϕμj≤ϵ is satisfied or ρ(A(Kj))<1-τ, stop; otherwise, go to the next step.
(3)
Choose a new barrier parameterμj+1∈(0,aμj) and a new inner-loop tolerance ϵj+1in∈(0,bϵjin).
(4)
Choose a new starting pointKj+1s and set j←j+1.
End (do)Remark 18.
The major task of Algorithm17 lies in Step (1) where an unconstrained minimization method has to be employed to obtain an approximate local solution up to the prescribed accuracy represented by the inner-loop tolerance ϵin. In the implementation, we consider the three-term CG method of Section 4 to calculate such a local solution.


## 8. Illustrative Numerical Examples

In this section, various test examples are provided to illustrate the performance of the proposed methods. Among the considered test problems are two examples of decentralized control systems and two examples for time-delay systems. A starting feasible point is required for the log-barrier method which might be obtained by the CG method applied on problem (1). The log-barrier method aims to obtain a local solution for the constrained problem (2) or at least achieves stabilizing output feedback that strictly lies within the unit circle.The methods are implemented using Matlab and all computations were carried out on a Laptop with 3.07 GHz and 1 GB RAM. Some of the considered test problems are for  continuous-time systems. The Matlab functionc2d from the control toolbox is employed to provide the corresponding discrete-time data matrices.As mentioned in Remark8, we stop the method as soon as the objective function fk is strictly less than one or a local solution is achieved. The CG method uses the sufficient decrease condition (23) with γ=10-4 for the backtracking line search.Example 1.
This test problem is borrowed from [23]. The constant data matrices for the corresponding discrete-time model are (58)A=1.00000.10000.005001.00000.1000001.0000,B=0.10520.00020.10500.00500.10000.1000,CT=100100,where the spectral radius of the system matrix A is 1.0. Starting from K0, a matrix of ones, the proposed CG method successfully converges to a local minimizer for the unconstrained problem (1). The achieved output feedback gain matrix is (59)K∗=-7.7969-7.9723-3.3475-2.7390,and the corresponding objective function value is f(K∗)=0.9387.
Table1 shows the convergence behavior of the considered CG method to a local solution of problem (1), which is also a stabilizing output feedback gain to the corresponding control system (5).Table 1
Convergence behavior of the CG method to a local solution for Example1.
k
f
(
K
k
)
∇
f
(
K
k
)
1
1.2783
e
+
000
1.8255
e
-
001
2
9.4060
e
-
001
5.9548
e
-
002
3
9.4044
e
-
001
5.1829
e
-
002
⋮
⋮
⋮
26
9.3874
e
-
001
1.5049
e
-
002
27
9.3874
e
-
001
1.0949
e
-
005
28
9.3874
e
-
001
6.4275
e
-
006Example 2.
This test problem is the aircraft model in cruise flight conditions [24, AC1] for a continuous-time control system. The functionc2d is used to have the discrete-time counterpart with the following data matrices: (60)A=1.00000.00140.11320.0005-0.096700.9945-0.0171-0.00050.006800.00031.00000.0957-0.004800.0060-0.00000.9131-0.09360-0.02770.00020.09730.9287,B=-0.00760.00000.0003-0.01150.09970.00000.02120.0000-0.00810.41520.0003-0.15890.1742-0.0014-0.0154,C=100000100000100.The spectral radius of the system matrix A is 1. Starting from K0, a matrix of ones, the CG method achieves the following stabilizing output feedback gain matrix after 12 iterations: (61)K∗=0.74730.74900.0293-0.2442-1.2308-2.29481.09501.13581.2428.The corresponding objective function value at both points is f(K0)=1.0841 and f(K∗)=0.9836.Example 3.
This test problem represents a chemical reactor model [24, REA1]. The data matrices of the corresponding discrete-time system are as follows: (62)A=1.17820.00150.5116-0.4033-0.05150.6619-0.01100.06130.07620.33510.56060.3824-0.00060.33530.08930.8494,B=0.0045-0.08760.46720.00120.2132-0.23530.2131-0.0161,C=101-10100001-1.The spectral radius of the system matrix A is 1.2203. The system is clearly Schur unstable. By starting with K0, the zero matrix, the CG method converges to a local solution K∗ after 25 iterations, where (63)K∗=0.3847-0.28280.33152.26670.44461.0987,and the corresponding objective function value is f(K∗)=0.9382.


### 8.1. The CG Method for Decentralized Control Systems

Example 4.
This test problem is borrowed from [25] of a decentralized control system with two control stations; each station has one input and one measured output. The given constant data matrices are (64)A=1.00-1.6-1.01.0-0.300.41.0,B1=100,B2=011,C1=100,C2=011.The uncontrolled system is discrete-time Schur unstable, where ρ(A)=1.8154>1. Starting with the following K0, the CG method requires only 6 iterations to converge to a stabilizing output feedback gain matrix K∗, where (65)K0=0.1000.1,K∗=-1.006000-0.4585.The objective function value at the two points is f(K0)=1.9924 and f(K∗)=0.6996<1, respectively. Although f(K∗) seems to be relatively small, a stationary point fails to be achieved.Example 5.
This is a fifth-order decentralized control system with two control stations. The data matrices are randomly generated which are as follows:(66)A=0.8755-0.1324-0.01050.17620.04840.02890.9512-0.0056-0.02570.0048-0.0516-0.01531.0470-0.09000.0207-0.09150.2109-0.02921.0087-0.09150.13050.06080.0939-0.13450.8591,B1=-0.12510.0472-0.05690.0821-0.03220.1360-0.11950.0421-0.02670.1732,B2T=0.05830.0432-0.1143-0.03400.01790.11860.01530.1403-0.1515-0.0926,C1=0100000100,C2=0000100000.The uncontrolled system is discrete-time Schur unstable, where ρ(A)=1.1114>1. Starting from the following randomly generated K0, the CG method reaches a stabilizing K∗ after 8 iterations, where (67)K0=-0.9513-1.8038000.60901.267100000.38700.3936000.26480.4900,K∗=9.4329-12.435200-23.61421.260600001.10311.1097001.27481.5000.The corresponding objective function value at both points is f(K0)=1.2056 and f(K∗)=0.9669, respectively.


### 8.2. The CG Method for Time-Delay Systems

Example 6.
This test problem represents a time-delay system (see [22]) which has the following constant data matrices: (68)A0=0.950.78-1.76-0.87,A1=0.120.89-0.95-2.97,B^=0.5-0.95,C0T=0.95-2.98,C1=00,where d=1. The data matrices for this system after converting it from continuous to discrete are as follows: (69)A=1.09280.08230.00870.0811-0.18050.8957-0.0914-0.29030.10470.00401.00050.0042-0.00900.0951-0.00460.9854,B=0.0485-0.09480.0025-0.0047,CT=0.9500-2.980000.The spectral radius of the system matrix A is 1.0098. Starting from the following K0, the three-term CG method successfully converges to a local solution of the minimization problem (1). The starting and achieved local solutions are (70)K0=-0.4838,K∗=-1.3741.The objective function value at both points is f(K0)=1.0008 and f(K∗)=0.9887, respectively.Example 7.
This test problem represents a time-delay system (see [22]) which has the following constant data matrices: (71)A0=1-0.60.40.5,A1=0.50.20.60.4,B^=0.1100.1,C0=1001,where d=1. By converting the system from continuous to discrete-time, the corresponding data matrices are (72)A=1.1065-0.06370.05070.01980.04631.05200.06260.04150.1052-0.00311.00250.00100.00220.10260.00311.0020,B=0.01050.02070.00020.01070.00050.00100.00000.0005,CT=10000100.The spectral radius of the system matrix A is 1.1183. Starting from K0, the three-term CG method successfully converges to a local solution for the minimization problem (1). The starting and achieved local solutions are as follows: (73)K0=-147.214216.259572.9190-39.6464,K∗=-162.54464.720752.8133-71.3982.The objective function value at both points is f(K0)=1.1136 and f(K∗)=0.9989, respectively.


### 8.3. The Logarithmic Barrier Method

The following two examples quite show the performance of the logarithmic barrier interior-point method for finding a local solution to problem (2). A starting feasible point for the method might be obtained by executing the CG method with some iterations on problem (1) until a stabilizing K is obtained.Example 8.
This test problem is borrowed from [23]. The constant data matrices for the corresponding discrete-time model are as follows: (74)A=1.00000.10000.005001.00000.1000001.0000,B=0.10520.00020.10500.00500.10000.1000,CT=100100,where the spectral radius of the system matrix A is 1.0. Starting from the following feasible K0, the interior-point method successfully converges to a local solution K∗. The two matrices are as follows: (75)K0=-7.79-7.97-3.34-2.73,K∗=-7.7592-7.9516-3.3145-2.7147,where the corresponding objective function values are f(K0)=0.9430 and f(K∗)=0.9388, respectively.Example 9.
This test problem represents a chemical reactor model [24, REA2]. The data matrices for the corresponding discrete-time system are as follows: (76)A=1.18050.00140.5122-0.4038-0.05150.6619-0.01100.06130.07630.33510.56060.3823-0.00060.33530.08930.8494,B=0.0045-0.08760.46720.00120.2132-0.23530.2131-0.0161,C=101-10100.The spectral radius of the system matrix A is 1.2227. Starting from the following feasible K0, the log-barrier method converges after 4 outer-loop iterations and 77 inner-loop iterations to a local solution K∗. The starting and achieved local solutions are as follows: (77)K0=0.07-0.491.53-0.09,K∗=-0.5680-0.01315.8580-0.8531.Table 2 shows the convergence behavior of the log-barrier interior-point method while converging to the local solution.Table 2
Convergence behavior of the log–barrier interior–point method to a local solution for Example9.
j
ϵ
j
i
n
μ
j
k
f
(
K
k
)
∇
ϕ
μ
j
K
k
0
9.0
e
-
001
1.0
e
-
001
0
9.4856
e
-
001
2.9173
e
-
001
1
9.4063
e
-
001
3.6122
e
-
001
1
9.0
e
-
002
4.0
e
-
003
0
9.4063
e
-
001
3.6122
e
-
001
1
9.3240
e
-
001
1.0758
e
-
001
2
9.3240
e
-
001
5.3765
e
-
001
⋮
⋮
⋮
61
8.2291
e
-
001
3.5946
e
-
001
62
8.2283
e
-
001
1.0708
e
-
002
2
9.0
e
-
003
6.4
e
-
006
0
8.2283
e
-
001
1.0708
e
-
002
⋮
⋮
⋮
5
8.2282
e
-
001
6.5932
e
-
003
3
9.0
e
-
004
4.1
e
-
010
0
8.2282
e
-
001
6.5932
e
-
003
⋮
⋮
⋮
6
8.2281
e
-
001
7.3226
e
-
004
4
9.0
e
-
005
1.0
e
-
015
0
8.2281
e
-
001
7.3226
e
-
004
⋮
⋮
⋮
3
8.2281
e
-
001
4.5061
e
-
005The above results quite show the convergence behavior of the proposed three-term CG method and the logarithmic barrier interior-point method. For all the considered test problems, the CPU time is less than 3 sec. If a local solution fails to be achieved, the CG method at least finds a stabilizing output feedback controller for the closed-loop system.


## 8.1. The CG Method for Decentralized Control Systems

Example 4.
This test problem is borrowed from [25] of a decentralized control system with two control stations; each station has one input and one measured output. The given constant data matrices are (64)A=1.00-1.6-1.01.0-0.300.41.0,B1=100,B2=011,C1=100,C2=011.The uncontrolled system is discrete-time Schur unstable, where ρ(A)=1.8154>1. Starting with the following K0, the CG method requires only 6 iterations to converge to a stabilizing output feedback gain matrix K∗, where (65)K0=0.1000.1,K∗=-1.006000-0.4585.The objective function value at the two points is f(K0)=1.9924 and f(K∗)=0.6996<1, respectively. Although f(K∗) seems to be relatively small, a stationary point fails to be achieved.Example 5.
This is a fifth-order decentralized control system with two control stations. The data matrices are randomly generated which are as follows:(66)A=0.8755-0.1324-0.01050.17620.04840.02890.9512-0.0056-0.02570.0048-0.0516-0.01531.0470-0.09000.0207-0.09150.2109-0.02921.0087-0.09150.13050.06080.0939-0.13450.8591,B1=-0.12510.0472-0.05690.0821-0.03220.1360-0.11950.0421-0.02670.1732,B2T=0.05830.0432-0.1143-0.03400.01790.11860.01530.1403-0.1515-0.0926,C1=0100000100,C2=0000100000.The uncontrolled system is discrete-time Schur unstable, where ρ(A)=1.1114>1. Starting from the following randomly generated K0, the CG method reaches a stabilizing K∗ after 8 iterations, where (67)K0=-0.9513-1.8038000.60901.267100000.38700.3936000.26480.4900,K∗=9.4329-12.435200-23.61421.260600001.10311.1097001.27481.5000.The corresponding objective function value at both points is f(K0)=1.2056 and f(K∗)=0.9669, respectively.


## 8.2. The CG Method for Time-Delay Systems

Example 6.
This test problem represents a time-delay system (see [22]) which has the following constant data matrices: (68)A0=0.950.78-1.76-0.87,A1=0.120.89-0.95-2.97,B^=0.5-0.95,C0T=0.95-2.98,C1=00,where d=1. The data matrices for this system after converting it from continuous to discrete are as follows: (69)A=1.09280.08230.00870.0811-0.18050.8957-0.0914-0.29030.10470.00401.00050.0042-0.00900.0951-0.00460.9854,B=0.0485-0.09480.0025-0.0047,CT=0.9500-2.980000.The spectral radius of the system matrix A is 1.0098. Starting from the following K0, the three-term CG method successfully converges to a local solution of the minimization problem (1). The starting and achieved local solutions are (70)K0=-0.4838,K∗=-1.3741.The objective function value at both points is f(K0)=1.0008 and f(K∗)=0.9887, respectively.Example 7.
This test problem represents a time-delay system (see [22]) which has the following constant data matrices: (71)A0=1-0.60.40.5,A1=0.50.20.60.4,B^=0.1100.1,C0=1001,where d=1. By converting the system from continuous to discrete-time, the corresponding data matrices are (72)A=1.1065-0.06370.05070.01980.04631.05200.06260.04150.1052-0.00311.00250.00100.00220.10260.00311.0020,B=0.01050.02070.00020.01070.00050.00100.00000.0005,CT=10000100.The spectral radius of the system matrix A is 1.1183. Starting from K0, the three-term CG method successfully converges to a local solution for the minimization problem (1). The starting and achieved local solutions are as follows: (73)K0=-147.214216.259572.9190-39.6464,K∗=-162.54464.720752.8133-71.3982.The objective function value at both points is f(K0)=1.1136 and f(K∗)=0.9989, respectively.


## 8.3. The Logarithmic Barrier Method

The following two examples quite show the performance of the logarithmic barrier interior-point method for finding a local solution to problem (2). A starting feasible point for the method might be obtained by executing the CG method with some iterations on problem (1) until a stabilizing K is obtained.Example 8.
This test problem is borrowed from [23]. The constant data matrices for the corresponding discrete-time model are as follows: (74)A=1.00000.10000.005001.00000.1000001.0000,B=0.10520.00020.10500.00500.10000.1000,CT=100100,where the spectral radius of the system matrix A is 1.0. Starting from the following feasible K0, the interior-point method successfully converges to a local solution K∗. The two matrices are as follows: (75)K0=-7.79-7.97-3.34-2.73,K∗=-7.7592-7.9516-3.3145-2.7147,where the corresponding objective function values are f(K0)=0.9430 and f(K∗)=0.9388, respectively.Example 9.
This test problem represents a chemical reactor model [24, REA2]. The data matrices for the corresponding discrete-time system are as follows: (76)A=1.18050.00140.5122-0.4038-0.05150.6619-0.01100.06130.07630.33510.56060.3823-0.00060.33530.08930.8494,B=0.0045-0.08760.46720.00120.2132-0.23530.2131-0.0161,C=101-10100.The spectral radius of the system matrix A is 1.2227. Starting from the following feasible K0, the log-barrier method converges after 4 outer-loop iterations and 77 inner-loop iterations to a local solution K∗. The starting and achieved local solutions are as follows: (77)K0=0.07-0.491.53-0.09,K∗=-0.5680-0.01315.8580-0.8531.Table 2 shows the convergence behavior of the log-barrier interior-point method while converging to the local solution.Table 2
Convergence behavior of the log–barrier interior–point method to a local solution for Example9.
j
ϵ
j
i
n
μ
j
k
f
(
K
k
)
∇
ϕ
μ
j
K
k
0
9.0
e
-
001
1.0
e
-
001
0
9.4856
e
-
001
2.9173
e
-
001
1
9.4063
e
-
001
3.6122
e
-
001
1
9.0
e
-
002
4.0
e
-
003
0
9.4063
e
-
001
3.6122
e
-
001
1
9.3240
e
-
001
1.0758
e
-
001
2
9.3240
e
-
001
5.3765
e
-
001
⋮
⋮
⋮
61
8.2291
e
-
001
3.5946
e
-
001
62
8.2283
e
-
001
1.0708
e
-
002
2
9.0
e
-
003
6.4
e
-
006
0
8.2283
e
-
001
1.0708
e
-
002
⋮
⋮
⋮
5
8.2282
e
-
001
6.5932
e
-
003
3
9.0
e
-
004
4.1
e
-
010
0
8.2282
e
-
001
6.5932
e
-
003
⋮
⋮
⋮
6
8.2281
e
-
001
7.3226
e
-
004
4
9.0
e
-
005
1.0
e
-
015
0
8.2281
e
-
001
7.3226
e
-
004
⋮
⋮
⋮
3
8.2281
e
-
001
4.5061
e
-
005The above results quite show the convergence behavior of the proposed three-term CG method and the logarithmic barrier interior-point method. For all the considered test problems, the CPU time is less than 3 sec. If a local solution fails to be achieved, the CG method at least finds a stabilizing output feedback controller for the closed-loop system.


## 9. Conclusion

The eigenvalue assignment problem for discrete-time systems formulated as an unconstrained optimization is considered. A three-term CG method is proposed to find a local solution of the minimization problem or at least to achieve a stabilizing output feedback gain matrix.The conjugate gradient method is extended to tackle the eigenvalue assignment problem for the case of decentralized control systems and the case of control systems with time delay. Global convergence is established for the conjugate gradient method.In order to maintain a selected stability margin for the control system, the optimization problem is reformulated as an inequality constraint problem by adding a cut to the objective function. A logarithmic barrier method is proposed for tackling this problem. The performance of the CG method as well as the logarithmic barrier method is shown on various test problems. It is quite interesting to see the performance of some variants of quasi-Newton and inexact Newton methods on the eigenvalue assignment problem for discrete-time systems.



---

*Source: 7256769-2017-08-10.xml*