c Pleiades Publishing, Ltd., 2017. ISSN 0005-1179, Automation and Remote Control, 2017, Vol. 78, No. 7, pp. 1203–1217.  c M.M. Tchaikovsky, V.N. Timin, 2017, published in Avtomatika i Telemekhanika, 2017, No. 7, pp. 39–56. Original Russian Text 

LINEAR SYSTEMS

Synthesis of Anisotropic Suboptimal Control for Linear Time-Varying Systems on Finite Time Horizon M. M. Tchaikovsky∗,∗∗,a and V. N. Timin∗∗,b ∗

∗∗

Academician Pilyugin Center, Moscow, Russia Institute of Control Sciences, Russian Academy of Sciences, Moscow, Russia e-mail: a [email protected], b [email protected] Received February 24, 2016

Abstract—This paper presents the statement and solution to the problem of synthesis of the anisotropic control that guarantees some prescribed level of attenuation of the uncertain stochastic disturbances affecting the linear discrete time-varying system on finite time horizon. The anisotropy of a random vector is considered as a measure of the statistical uncertainty of the disturbance. The closed-loop system capabilities to attenuate the external disturbances are characterized by its anisotropic norm. The synthesis problem solution is formulated in form of sufficient conditions of existence of a controller that guarantees the anisotropic norm of the closed-loop system to be bounded by some given threshold level. The controller synthesis algorithm is based on solving the system of matrix inequalities recursively. Keywords: linear discrete time-varying system, finite time horizon, attenuation of stochastic disturbances, statistic uncertainty, anisotropy, norm, static output feedback, matrix inequalities, recursive solution. DOI: 10.1134/S0005117917070037

1. INTRODUCTION The paper considers the problem of control of a linear discrete-time system with time-varying parameters affected by correlated random disturbances with the unknown distribution. The timevarying parameters of the control object are the known functions of time. The goal of the feedback control consists in attenuation of the random disturbances affecting the closed-loop system. Application of the classical approach of optimal linear-quadratic Gaussian (LQG) control for solving this problem supposes that covariance matrices of the random disturbances are known precisely [1]. From the other hand, the H∞ control [2] designed for a worst-case scenario of the deterministic square-summable disturbance that guarantees the 2 -induced operator norm of the closed-loop system to be bound can be rather conservative and require excessive energy to solve this problem. In general case, when the probability distribution of the correlated disturbances is not known precisely, one can consider its anisotropy, i.e., the relative entropy of its distribution w.r.t. a family of Gaussian distributions with zero mean and scalar covariance matrices of a kind λI, λ > 0 [3], as a measure of its uncertainty. In this case the worst-case performance of the system is characterized by its anisotropic norm defined as the maximum mean-square gain of the system w.r.t. the disturbances with the anisotropy bounded above by some known number a > 0 [3]. In two limiting cases, when a → 0 and a → ∞, the anisotropic norm of the system becomes its scaled Frobenius norm and 2 -induced operator norm, respectively. The anisotropy-based theory of stochastic robust control, which is mainly set forth in [4–6] as applied to the linear discrete time-invariant systems on infinite time horizon, is aimed at synthesis 1203

1204

TCHAIKOVSKY, TIMIN

of the control that is more robust with respect to the correlated disturbances than LQG and H2 control, and at the same time less conservative than H∞ control. Application of the anisotropic norm as a performance criterion in the optimal control problem led to synthesis of the optimal anisotropic controller based on solving a system of cross-coupled nonlinear matrix algebraic equations [6]. Formulation of a criterion for the anisotropic norm to be bounded by a given threshold value in terms of linear matrix inequalities (LMIs) [7] allowed to apply the methods of convex optimization and LMI technique for solving the suboptimal anisotropic control and filtering problems. Definitions of the anisotropy of a random vector and anisotropic norm of a linear discrete timevarying (LDTV) system introduced in [3] founded the basis for solving the anisotropic control and filtering problems on finite time horizon. Necessary and sufficient conditions for verifying whether the anisotropic norm of the system with time-varying parameters is bounded above by a given threshold value were presented in [8]. These conditions require solving the difference Riccati equation and inequality with respect to the determinant of some positive definite matrix. Based on this criterion, a solution to the anisotropic filtering problem was proposed in [9]. In this paper we solve the problem of synthesis of the anisotropic controller for the LDTV system that guarantees the anisotropic norm of the closed-loop system to be bounded by some given threshold value on finite time horizon. The paper has the following structure. In Section 2, the problem statement is set forth. In Section 3, the control problem solution is derived. Section 4 considers a numerical example. Concluding remarks are given in Section 5. The following notations are used through the paper. Rn denotes the set of the real n-dimensional vectors, Rn×m stands for the set of the real (n × m) matrices. For the real matrix M = [mij ], M T denotes its transpose, M T := [mji ]. For the real symmetric matrices M = M T , N = N T , M  N (or M  N ) means that the matrix M − N is positive definite (or semi-definite). In the block-symmetric matrices, the symbol ∗ substitutes the blocks defined by symmetry. The (n × n) identity matrix  is denoted by In . The maximum singular value of the matrix M is denoted by σ(M ) := λmax (M T M ), where λmax (M T M ) stands for the maximum eigenvalue of the ma trix M T M . The trace of the matrix M = [mij ] is denoted by tr M := k mkk , its determinant is denoted by det M . The Euclidian norm of a vector or the absolute value of a scalar value is denoted by | · |. The symbol E stands for the expectation of a random value. 2. PROBLEM STATEMENT Let us consider a linear discrete time-invariant system P described on the finite horizon [0, N ] of the discrete time k by the state-space model with the nx -dimensional state X = (xk ), mw -dimensional disturbance input W = (wk ), mu -dimensional control input U = (uk ), pz -dimensional controlled output Z = (zk ), and py -dimensional measured output Y = (yk ) as ⎡

⎤

⎡

⎤⎡

⎤

xk+1 Ak Bw,k Bu,k xk ⎢ ⎥ ⎢ ⎥⎢ ⎥ P : ⎣ zk ⎦ = ⎣ Cz,k Dzw,k Dzu,k ⎦ ⎣ wk ⎦ , yk Cy,k Dyw,k 0 uk

k = 0, . . . , N,

(2.1)

with the initial state x0 = 0, where Ak , Bw,k , Bu,k , Cz,k , Dzw,k , Dzu,k , Cy,k , Dyw,k are the known properly-dimensioned matrix-valued functions of the discrete time k. The signals in this system are understood as the sequences of vectors. It is supposed that the system P is controllable in full, i.e., for each time point i ∈ [0, N ] there exists such time point j  i + 1 that the symmetric positive semi-definite matrix Πij =

j−1

k=i

T Φj,k+1Bu,k Bu,k ΦT j,k+1

AUTOMATION AND REMOTE CONTROL

Vol. 78

No. 7

2017

SYNTHESIS OF ANISOTROPIC SUBOPTIMAL CONTROL

1205

is nonsingular [1]. Here Φjk := Aj−1 Aj−2 . . . Ak when j > k, Φkk := Inx when j = k is the state transition matrix of the undisturbed system (2.1) from xk to xj . The control goal is to ensure some desirable quality of attenuation of the external random disturbances. Let us find the control law as the static output feedback with the time-varying matrix gain that depends on the discrete time uk = Kk yk ,

k = 0, . . . , N.

(2.2)

The control law (2.2) represents rather wide class of the controllers that easily allows to turn to a number of particular cases. For example, if in the model (2.1) of the control object the matrices Cy,k = Inx , Dyw,k = 0, we get the time-varying state-feedback controller uk = Kk xk . The control law (2.2) also includes the problem of designing the controller in form of the time-varying fixed-order dynamic compensator 

xck+1





=

uk

Ack Bkc



Ckc Dkc

xck



,

yk

k = 0, . . . , N,

(2.3)

where xck ∈ Rnc , xc0 = 0, by expanding the state vector of the control object (2.1) with the state of the controller (2.3) and considering the augmented model with the matrices ⎡ ⎢

⎡

Ak

0

Bw,k

0

Bu,k

⎤ ⎢ 0 0 0 In c 0 ⎢ Ak Bw,k Bu,k ⎢ ⎢ ⎥ ⎣ Cz,k Dzw,k Dzu,k ⎦ := ⎢ ⎢ Cz,k 0 Dzw,k 0 Dzu,k ⎢ Cy,k Dyw,k 0 ⎢ 0 In c 0 0 0 ⎣

0

Cy,k

Dyw,k

0

⎤ ⎥ ⎥ ⎥ ⎥ ⎥, ⎥ ⎥ ⎥ ⎦



Kk :=

Ack Bkc



Ckc Dkc

0

as the control object and controller. The equations of the closed-loop system F consisting of the control object (2.1) and controller (2.2) are as follows 

F : where



Ak Bk Ck Dk

xk+1 zk





:=





=

Ak Bk Ck Dk



xk wk



,

k = 0, . . . , N,

Ak + Bu,k Kk Cy,k Bw,k + Bu,k Kk Dyw,k Cz,k + Dzu,k Kk Cy,k Dzw,k + Dzu,k Kk Dyw,k

(2.4)

.

(2.5)

Gathering the elements of the input and output signals W and Z of the closed-loop system on the time interval [0, N ] to the column vectors ⎡

W0:N

⎤

w0 ⎢ .. ⎥ := ⎣ . ⎦ , wN

⎡

Z0:N

⎤

z0 ⎢ .. ⎥ := ⎣ . ⎦ zN

(2.6)

let us identify the input-output operator of the closed-loop system F with the matrix F0:N being the lower triangular matrix F0:N := block [fjk ] 0j,kN

AUTOMATION AND REMOTE CONTROL

Vol. 78

No. 7

2017

(2.7)

1206

TCHAIKOVSKY, TIMIN

composed of the blocks fjk ∈ Rpz ×mw defined by

⎧ ⎪ ⎨ Cj Φj,k+1 Bk when j > k

Dk ⎪ ⎩ 0

fjk :=

when j = k when j < k,

where Φjk := Aj−1 Aj−2 . . . Ak when j > k, Φkk := Inx when j = k is the state transition matrix of the closed-loop system (2.4) from xk to xj . Since the closed-loop system F has zero initial state, its vectorized output on the finite time horizon [0, N ] is given by Z0:N = F0:N W0:N .

(2.8)

In this connection, any norm of the system F is understood as the respective norm of the matrix F0:N . In particular, the analogues of the H2 and H∞ norms of the DLTV system F on the finite time horizon [0, N ] are the Frobenius and induced norms of the matrix F0:N given by [3] F 2 :=



T F tr F0:N 0:N ,

F ∞ := σ(F0:N ).

(2.9)

Denote the dimension of the vector W0:N by  := mw (N + 1). It is assumed that the column vector W0:N belongs to the class L2 of square-integrable R -valued random vectors distributed absolutely continuously relative to the -dimensional Lebesgue measure mes  . Recall that the anisotropy A(W0:N ) of the vector W0:N ∈ L2 with the probability distribution function f : R → R+ is defined in [3] as the minimum value of the relative entropy D(f p,λ ) with respect to Gaussian distributions p,λ in R with the zero mean and scalar covariance matrices λI : 

λ>0



 2πe ln E|W0:N |2 − h(W0:N ), 2 

A(W0:N ) := min D(f p,λ ) = 

(2.10)

where h(W0:N ) := −E ln f (W0:N ) = − R f (w) ln f (w)dw is the differential entropy of W0:N [10]. Hereinafter it is assumed that the disturbance anisotropy is bounded above by some given parameter a > 0: A(W0:N )  a.

(2.11)

To achieve the declared control objective, let us quantitatively characterize the ability of the closed-loop system F to attenuate the external disturbance W by its anisotropic norm defined in [3] as 



|||F |||a = sup R(F0:N , W0:N ) : W0:N ∈ L2 , A(W0:N )  a , where R(F0:N , W0:N ) :=

(2.12)



E (|Z0:N |2 ) /E (|W0:N |2 )

is the mean-square gain of the operator F with respect to W. From [3] it is known that for any fixed system F its anisotropic norm is a nondecreasing function of the parameter a  0 that satisfies the relations √ (2.13) F 2 /  = |||F |||0  |||F |||a  lim |||F |||a = F ∞ . a→+∞

The statement of the anisotropic suboptimal control synthesis problem for LDTV system on the finite time horizon is formulated as follows. AUTOMATION AND REMOTE CONTROL

Vol. 78

No. 7

2017

SYNTHESIS OF ANISOTROPIC SUBOPTIMAL CONTROL

1207

Problem. Given the control plant P with the state-space model (2.1), anisotropy level a > 0 of its input W, and some desirable threshold value γ > 0, find the matrix-valued gain function K = (Kk ), k = 0, . . . , N , of the static output-feedback controller (2.2) that guarantees that the anisotropic norm (2.12) of the closed-loop system F does not exceed the threshold value γ, that is |||F |||a < γ.

(2.14)

3. SOLUTION OF ANISOTROPIC CONTROLLER SYNTHESIS PROBLEM To solve the problem stated above, let us make use of the following criterion for the anisotropic norm of LDTV system to be bounded by some given value formulated and proved in [8]: Lemma. Given γ > 0 and a  0, the anisotropic norm (2.12) of the linear discrete time-varying system (2.4) satisfies the condition |||F |||a  γ on the finite time horizon [0, N ] if and only if there k = R  T  0, exists the number q  0 such that for the symmetric positive semi-definite matrices R k k = 0, . . . , N , satisfying the difference Riccati equation T T   T R k+1 = Ak Rk Ak + qBk Bk + Mk Sk Mk ,

(3.1)

 CT + qB DT )S −1 , Mk = −(Ak R k k k k k  k CT − qDk DT S k = Ip − C k R k

z

(3.2) (3.3)

k

 0 = 0 the matrices S0 , . . . , SN are all positive definite and satisfy the with the initial condition R inequality N

ln det Sk   ln(1 − qγ 2 ) + 2a,

(3.4)

k=0

where  = mw (N + 1). In [8] it is noted that the matrices Sk defined by (3.3) are positive definite if and only if q < F −2 ∞ . For any such q the left part of the inequality (3.4) is nonpositive since Sk  Ipz , hence any q that satisfies the conditions of the Lemma must also satisfy the inequalities 



γ −2 1 − e−2a/  q < γ −2 .

(3.5)

Direct application of the Lemma conditions to the realization of the closed-loop system (2.4), (2.5) doesn’t lead to a relatively simple constructive solution of the synthesis problem. Therefore, let us proceed from the Lemma conditions in form of the Eq. (3.1)–(3.3) and inequality (3.4) to conditions in form of the matrix inequalities. From the monotonicity property of a solution of the difference Riccati equation w.r.t. the initial condition and the equation coefficient matrices it fol k of the Eqs. (3.1)–(3.3), k = 0, . . . , N , with the initial condition R 0 lows [11, 12] that the solutions R T exist, if there exist the symmetric positive semi-definite matrices Rk = Rk  0, k = 0, . . . , N , that satisfy the difference Riccati inequality T Rk+1  Ak Rk AT k + qBk Bk



T + Ak Rk CT k + qBk Dk



T Ipz − Ck Rk CT k − qDk Dk

−1 

T Ak Rk CT k + qBk Dk

T

(3.6)

 0 for the same value of q, at that Rk  R  k , k = 0, . . . , N . It is with the initial condition R0  R obvious that for such matrices Rk the inequality N

k=0









T 2 ln det Ipz − Ck Rk CT + 2a k − qDk Dk   ln 1 − qγ

AUTOMATION AND REMOTE CONTROL

Vol. 78

No. 7

2017

(3.7)

1208

TCHAIKOVSKY, TIMIN

T holds true since Ipz − Ck Rk CT k − qDk Dk  Sk ∀k ∈ [0, N ]. After some elementary transformations let us rewrite (3.7) as N  k=0

T 2a 2 mw (N +1) det (Ipz − Ck Rk CT . k − qDk Dk )  e (1 − qγ )

(3.8)

In the synthesis problem statement the anisotropy level of the disturbance is supposed to be strictly positive (a > 0). It is equivalent to the strict inequality q > 0 holding true, which implies 0 ≺ Sk ≺ Ipz and strict inequality in the localization condition (3.5), as well as in (3.8). Let us note that q = 0 if and only if a = 0 that corresponds to one of the limiting cases of the anisotropic norm, namely the Frobenius norm in the inequalities (2.13). The obvious sufficient condition for the inequality (3.8) to hold true is simultaneous satisfying the conditions 





T 2a/(N +1) 1 − qγ 2 det Ipz − Ck Rk CT k − qDk Dk > e

mw

,

k = 0, . . . , N.

(3.9)

 0 = 0, we obtain the inequality If for k = 0 in (3.9) we choose R0 = R 





2a/(N +1) 1 − qγ 2 det Ipz − qD0 DT 0 >e

mw

,

from which, taking into account the corollary of Schur’s Lemma (see e.g., [13]) 







T det (ηImw ) det Ipz − D0 (ηIm )−1 DT 0 = det ηImw − D0 D0 ,

where η := q −1 , we get 

1/mw

η − e−2a/(N +1) det(ηImw − DT 0 D0 )

< γ2.

(3.10)

Introduce the slack variable, (mw × mw ) matrix Θ0 = ΘT 0 such that 

ηImw −

DT 0 D0

 Θ0  0

⇔

ηImw − Θ0 DT 0 D0 Ip z



 0, Θ0  0.

(3.11)

The inequality (3.10) holds true if the inequalities (3.11) and 

η − e−2a/(N +1) det Θ0

1/mw

< γ2

(3.12)

hold. Now consider the Riccati inequality (3.6) when k = 0. If we choose the initial condition R0 =  0 = 0, then the inequality (3.6) becomes R 

2 T T R1 − qB0 BT 0 − q B0 D0 Ipz − qD0 D0

−1 

B0 DT 0

T

 0.

(3.13)

Let us demand the condition R1  R0 = 0 to hold true when q > 0. Applying Schur’s Lemma to the inequality (3.13) subject to η = q −1 > 0 we get ⎡

⎤

R1 0 B0 ⎢ ⎥ D0 ⎦  0, ⎣ 0 Ip z T BT 0 D0 ηImw

R1  0.

AUTOMATION AND REMOTE CONTROL

(3.14)

Vol. 78

No. 7

2017

SYNTHESIS OF ANISOTROPIC SUBOPTIMAL CONTROL

1209

The inequalities (3.11), (3.12), (3.14) define the sufficient conditions for the anisotropic norm of the system F to be bounded by γ at the zero point of the discrete time k. Now let us derive the similar sufficient conditions for the time interval [1, N ]. The inequalities (3.9) written w.r.t. the new variables η = q −1 and Pk = ηRk become 





T 2a/(N +1) pz −mw η η − γ2 det ηIpz − Ck Pk CT k − Dk Dk > e

mw

,

k = 1, . . . , N.

(3.15)

Introduce the slack variable, (pz × pz ) matrix Ψk = ΨT k such that T ηIpz − Ck Pk CT k − Dk Dk  Ψk  0

or, by Schur’s Lemma, ⎡ ⎢ ⎢ ⎢ ⎣

ηIpz − Ψk

Ck

CT k

Pk−1

DT k

0

⎤

Dk

⎥ ⎥ ⎦

0 ⎥  0,

Ψk  0,

k = 1, . . . , N.

(3.16)

Imw

The inequalities (3.15) hold true if the inequalities (3.16) and det Ψk > e2a/(N +1) η pz −mw (η − γ 2 )mw ,

k = 1, . . . , N,

(3.17)

hold. Denoting Pk := ηRk and applying Schur’s Lemma to the inequality (3.6) we obtain the inequalities similar to (3.6) for k = 1, . . . , N : ⎡

Pk+1

⎢ ⎢ 0 ⎢ ⎢ ⎢ ⎢ AT k ⎣

BT k

⎤

0

Ak

Bk

ηIpz

Ck

Dk ⎥ ⎥

CT k

Pk−1

0 ⎥

DT k

0

⎥ ⎥  0, ⎥ ⎦

Pk  0.

(3.18)

Imw

Multiplying (3.18) and (3.16) from the both sides by the symmetric positive definite block-diagonal matrices blockdiag(Inx , Ipz , Pk , Imw ), blockdiag(Ipz , Pk , Imw ), respectively, for k = 1, . . . , N we have ⎡

Pk+1

0

Ak Pk

⎢ ⎢ 0 ηIpz Ck Pk ⎢ ⎢ ⎢ T ⎢ Pk AT Pk k Pk Ck ⎣

BT k

DT k

0

Bk

⎤

⎥ Dk ⎥ ⎥ ⎥  0, ⎥ 0 ⎥ ⎦

Imw

⎡ ⎢ ⎢ ⎢ ⎣

ηIpz − Ψk Ck Pk Pk CT k

Pk

DT k

0

Dk

⎤ ⎥ ⎥ ⎦

0 ⎥  0, Imw

Pk  0, Ψk  0.

(3.19)

The inequalities (3.17) and (3.19) establish the sufficient conditions for the anisotropic norm of the system F to be bounded by given threshold value γ on the finite time horizon [1, N ]. Joining up these inequalities with the inequalities (3.12), (3.14) for k = 0 and application to the closed-loop system realization (2.4), (2.5) leads to the solution of the synthesis problem formulated as the following theorem. Theorem. Given a > 0 and γ > 0, the static output-feedback controller (2.2) with the feedback gain coefficient matrix K = (Kk ), k = 0, . . . , N , solves the synthesis problem for the control obT ject (2.1) if there exist a scalar η > 0, (mw × mw )-matrix Θ0 = ΘT 0 , and (nx × nx )-matrix R1 = R1 AUTOMATION AND REMOTE CONTROL

Vol. 78

No. 7

2017

1210

TCHAIKOVSKY, TIMIN

such that 

η − e−2a/(N +1) det Θ0 

1/mw

< γ 2,

T T K TDT ηImw − Θ0 Dzw,0 + Dyw,0 0 zu,0

∗ ⎡

∗

 0,

Ip z

R1

⎢ ⎢ ∗ ⎣

(3.20)

0

Bw,0 + Bu,0 K0 Dyw,0

(3.21)

⎤ ⎥

Ipz Dzw,0 + Dzu,0 K0 Dyw,0 ⎥ ⎦  0, ∗

(3.22)

ηImw R1  0,

Θ0  0,

(3.23)

T as well as (pz × pz )-matrices Ψk = ΨT k and (nx × nx )-matrices Pk = Pk , k = 1, . . . , N, that satisfy the system of inequalities



e−2a/(N +1) det Ψk

⎡ ⎢ ⎢ ⎣ ⎡ ⎢ ⎢ ⎢ ⎢ ⎣



1 pz



> η 1 − γ 2 /η

 mw pz

,

ηIpz − Ψk (Cz,k + Dzu,k Kk Cy,k )Pk Dzw,k + Dzu,k Kk Dyw,k

Pk+1 ∗

∗

Pk

0

∗

∗

Imw

(Ak + Bu,k Kk Cy,k )Pk

Bw,k + Bu,k Kk Dyw,k

0

(3.24) ⎤ ⎥ ⎥  0, ⎦

(3.25)

⎤ ⎥ ⎥  0, ⎥ ⎦

ηIpz (Cz,k + Dzu,k Kk Cy,k )Pk Dzw,k + Dzu,k Kk Dyw,k ⎥

∗

∗

Pk

0

∗

∗

∗

Imw

Pk  0,

Ψk  0,

(3.26)

k = 1, . . . , N,

(3.27)

with the initial condition P1 = ηR1 . 1

From [14] it is known that the function −(det Ψk ) pz of the real (pz × pz )-matrix Ψk = ΨT k  0 is convex w.r.t. its argument and representable as an intersection of finite number of second-order cones, i.e., as the system of linear matrix inequalities (LMIs). The sufficient conditions of the Theorem divide the anisotropic control synthesis procedure into two stages. First, from solving the system of inequalities (3.20)–(3.23) for k = 0 one should compute the scalar η  = (q  )−1 > 0 and matrix R1  0 that define the initial condition P1 = η  R1 of the system of difference inequalities (3.24)–(3.27), as well as the matrix of the controller parameters K0 . Then the inequalities (3.24)–(3.27) are solved recursively for k = 1, 2, . . . , N with fixed η = η  w.r.t. the unknown variables Pk+1  0, Ψk  0, and Kk . As it was noted before, from the monotonicity property of the solution of the difference Riccati  k  0, equation w.r.t. the initial conditions and equation matrix coefficients it follows that Rk  R  k = 0, . . . , N, where Rk is the solution of the Riccati equation (3.1)–(3.3) with the initial condition  0 = 0, and Rk is the solution of the respective Riccati inequality (3.6) with the initial condition R  0 . Therefore, for the inequalities (3.24)–(3.27) of the Theorem with some fixed value of R0  R   0 = 0, the η = η computed from solving the system of inequalities (3.20)–(3.23) derived for R0 = R    holds true for all k = 1, 2, . . . , N . It is known that for the positive inequality Pk = η Rk  η R k  from R  R  it follows that tr R  tr R  [13]. It can be intersemi-definite matrices Rk and R k k k k k min esting to find the minimum solution Rk , which is the most close to the solution of the Eq. (3.1) AUTOMATION AND REMOTE CONTROL

Vol. 78

No. 7

2017

SYNTHESIS OF ANISOTROPIC SUBOPTIMAL CONTROL

1211

over the set of all solutions of the inequality (3.6) that satisfy the inequality (3.7) for some fixed q  = (η  )−1 and synthesize the control that corresponds to this minimum solution. The sufficient conditions of the Theorem allow to do so by means of recursive solution of the convex optimization problems. Corollary. In conditions of the Theorem, the anisotropic suboptimal controller K = (Kk ) that corresponds to the minimum solution of the Riccati inequality (3.6) can be obtained from solving the convex optimization problems consecutively. First, solve the problem tr R1 → inf η, Φ0 , R1 , K0

over

satisfying

(3.28)

(3.20)–(3.23)

to obtain η  , R1min , and K0 . After the matrix P1 = η  R1min is computed, choose it as the initial condition and solve the convex optimization problems tr Pk+1 → inf Ψk , Pk+1 , Kk

over

satisfying

(3.29)

(3.24)–(3.27)

for k = 1, 2, . . . , N recursively. Remark 1. Simultaneous holding true of (3.20), (3.24) for all k is the sufficient condition for the inequality N 

(det Ψk )

1 pz

>e

2a pz

 

2

η 1 − γ /η

 mw N +1 pz

(3.30)

k=0

to hold true, but not necessary. It brings certain conservatism to the conditions of the Theorem. From the other hand, application of (3.30) in the conditions of the Theorem instead of the inequalities (3.20), (3.24) complicates the computational algorithm of the synthesis problem solution, because in this case one has to apply the recursive procedure

(det Ψk )

1 pz

>e

2a pz

 

η 1 − γ 2 /η

 mw N +1 pz

⎛ ⎝

k−1 

⎞−1

(det Ψj )1/pz ⎠

(3.31)

j=0

to construct the sequence of its solutions Ψ = (Ψk ), k = 0, . . . , N . First, solving the inequalities (3.31) recursively doesn’t allow to get η  immediately as a value of the variable that enters the system of inequalities linearly, so the inequalities (3.31) are not convex w.r.t. η. Second, even relatively small levels of the external disturbance anisotropy a can lead to appearance of big numbers in the right part of the inequality (3.31) and make difficult a practical implementation of the computational algorithm of the controller synthesis. Remark 2. The matrices of the controller parameters Kk , k = 0, . . . , N , enter the Theorem inequalities affinely making it possible to impose some additional structural constraints on them to design decentralized or fixed-structure controllers. Remark 3. If the external disturbance anisotropy level a → +∞, then from the inequalities (3.5) it follows that q → γ −2 and η → γ 2 . In this case the inequalities (3.20), (3.21), and (3.24), (3.25) AUTOMATION AND REMOTE CONTROL

Vol. 78

No. 7

2017

1212

TCHAIKOVSKY, TIMIN

hold true for any Θ0 , Φk  0, and the LMIs (3.22), (3.26) become ⎡ ⎢ ⎢ ⎣ ⎡ ⎢ ⎢ ⎢ ⎢ ⎣

γ 2 Ipz (Cz,k + Dzu,k Kk Cy,k )Pk Dzw,k + Dzu,k Kk Dyw,k ∗

Pk

0

∗

∗

Imw

Pk+1

0

(Ak + Bu,k Kk Cy,k )Pk

⎤ ⎥ ⎥  0, ⎦

Bw,k + Bu,k Kk Dyw,k

(3.32)

⎤ ⎥ ⎥  0, ⎥ ⎦

γ 2 Ipz (Cz,k + Dzu,k Kk Cy,k )Pk Dzw,k + Dzu,k Kk Dyw,k ⎥

∗ ∗

∗

Pk

0

∗

∗

∗

Imw Pk  0,

P0 = 0,

(3.33)

k = 0, . . . , N.

(3.34)

Existence of the solution of the system of recursive LMIs (3.32)–(3.34) is the necessary and sufficient condition for the inequality F ∞ < γ to hold true, where the induced norm of the closed-loop system is defined by the expression (2.9). 4. NUMERICAL EXAMPLE Let us consider the problem of control for the 2nd order linear oscillator 

x˙ 1 x˙ 2





=

0



1

−ω02 (t) −2ξ(t) ⎡

⎤

⎡

x1 x2





+

⎤

0 0 1 0



⎡

w1 w2





+

0 1



u,

(4.1)

⎤

z1 1 0  0 ⎢ ⎥ ⎢ ⎥ x1 ⎢ ⎥ + ⎣ 0 ⎦ u, ⎣ z2 ⎦ = ⎣ 0 1 ⎦ x2 z3 0 0 1 

y=



0 1



x1 x2





+



0 0, 1



w1 w2

(4.2)

(4.3)

with the time-varying parameters, namely, the natural frequency ω0 (t) and the damping factor ξ(t): ω0 (t) = 1 + 155| sin(3πt/T )|,

ξ(t) = −3t/T

(4.4)

on the finite time horizon [0, T ], T = 2 s. The model (4.1)–(4.4) was transformed to the discretetime form with the time step 0.004 s and written in the standard form (2.1). It is supposed that the anisotropy of the disturbance W is bounded by the number a > 0, A(W0:N ) < a. The typical diagrams of the disturbance signals for a = 500, 1000 are presented at Fig. 1. The random disturbances drive the open-loop system without control (4.1)–(4.4) to oscillations with variable frequency and increasing magnitude, see Fig. 2b. Based on recursive solving the system of inequalities (3.20)–(3.27) for the discrete-time system (4.1)–(4.4) we designed the control laws (2.2) in form of the static output feedback guaranteeing that |||F |||a < γ holds true for various given a, γ > 0. For the purpose of comparison, from solving the recursive LMIs (3.32)–(3.34) we also found the control laws (2.2) guaranteeing F ∞ < γ holds. While computing the controllers, the minimum solutions Pkmin were found according to the Corollary. Computations were realized in Matlab by means of YALMIP interface and SDPT3 solver [15]. Some results of the synthesis problem solution for various a, γ are presented in the table. For the closed-loop systems with the controllers computed for a√number of values of a and γ listed in the table, the values of the scaled Frobenius norm F0:N 2 /  and induced norm F0:N ∞ AUTOMATION AND REMOTE CONTROL

Vol. 78

No. 7

2017

SYNTHESIS OF ANISOTROPIC SUBOPTIMAL CONTROL

1213

Fig. 1. The random disturbances with different anisotropy level: (a) a = 500; (b) a = 1000.

Fig. 2. (a) The time-varying parameters of the control object: natural frequency (top plot), damping factor (bottom plot); (b) the trajectories of the open-loop system influenced by the disturbances (a = 1000).

were calculated by the formulas (2.9). For the considered control object and the values presented in the table these norms of the closed-loop system increase as a, γ, increase. Their maximum values correspond to the limiting case a → +∞ (see the first two rows of the table). The next three rows of the table present the values of η  found from solving the inequalities (3.20)–(3.23), as well as the left and right boundaries of the intervals of allowed values of η for the listed a and γ. These boundaries are defined by the inequalities (3.5) and, taking into account η = q −1 , are given by γ 2 < η  γ 2 (1 − e−2a/ )−1 , where  = mw (N + 1). All values of η  presented in the table belong to the allowed intervals. These intervals become narrower as a increases, and the computed values of η  become closer to their right boundaries. The maximum and minimum values of the function (det Ψk )1/pz over k = 1, . . . , N presented in the table show that the range of its changing grows as a increases and e−2a/(pz (N +1)) decreases. AUTOMATION AND REMOTE CONTROL

Vol. 78

No. 7

2017

1214

TCHAIKOVSKY, TIMIN

Fig. 3. The values of function (det Ψk )1/pz (top plots) and eigenvalues of the matrix Pk (middle and bottom plots) computed for the different anisotropy levels: (a) a = 500, (b) a = 1000, for γ = 32.

The last amount plays a role of the constant coefficient in the inequality (3.24). From the table data one can see that the function (e−2a/(N +1) det Ψk )1/pz for the listed values of a and γ is bounded !m /p below by η  1 − γ 2 /η  w z as consistent with (3.24). The illustration of this fact can be observed at the plots of Fig. 3, where the diagrams of the function (det Ψk )1/pz are represented for a = 500, a = 1000, and γ = 32 together with the diagrams of the eigenvalues of the matrix Pk denoted by the solid lines and symbols . For the purpose of comparison, these plots also represent the diagrams of the eigenvalues of the matrix Pk computed as a → +∞ and denoted by the dashed lines and symComparison of the closed-loop systems with the static-feedback controllers computed for various values of a and γ a 10 50 100 500 1000 → +∞ γ 16 18 20 32 32 32 Solution results: √ F0:N 2 /  F0:N ∞

0.5959 4.4534

0.8001 4.2397

0.814 4.2331

0.8789 4.2154

1.0462 4.3056

1.1807 5.5906

γ2 η 2 γ /(1 − e−2a/ ) maxk (det Ψk )1/pz mink (det Ψk )1/pz e−2a/(pz (N +1)) !mw /pz η  1 − γ 2 /η  maxk λ1 (Pk ) maxk λ2 (Pk ) maxk |Kk |

256 12934.95 12954.03 12934.86 12934.69 0.9868 12763.71 0.0067 1.94×10−6 8.615

Simulation results: maxk |uk | maxk |x1,k | maxk |x2,k |

3.022 4.864 5.051 4.323 5.092 6.647 18.97×10−4 16.2×10−4 14.9×10−4 10.2×10−4 13.6×10−4 14.3×10−4 0.0984 0.082 0.0796 0.0553 0.0734 0.0913

324 400 1024 1024 1024 3403.262 2206.661 1621.149 1184.537 γ2 3411.174 2210.649 1621.83 1185.016 1024 3403.066 2206.453 1620.875 1183.82 — 3402.708 2206.074 1620.372 1182.505 — 0.9356 0.8754 0.5141 0.2643 0 3183.491 1931.181 833.0324 312.5356 0 0.0038 0.0037 0.0033 0.0027 0.0022 1.69×10−6 1.69×10−6 1.67×10−6 1.67×10−6 1.66×10−6 12.56 12.94 14.88 23.6 48.96

AUTOMATION AND REMOTE CONTROL

Vol. 78

No. 7

2017

SYNTHESIS OF ANISOTROPIC SUBOPTIMAL CONTROL

1215

Fig. 4. The time-varying feedback gain Kk , computed for different anisotropy levels: (a) a = 500, (b) a = 1000, for γ = 32.

Fig. 5. The control signals uk = Kk yk computed for different anisotropy levels: (a) a = 500, (b) a = 1000, for γ = 32.

bols ♦. The diagrams of λ1,2 (Pk ) that correspond to the finite values of the disturbance anisotropy a practically coincide with the diagrams of the eigenvalues of the matrix Pk computed as a → +∞ for that points of the discrete time k, where the function (e−2a/(N +1) det Ψk )1/pz doesn’t achieve its minimum value conditioned by the constraint (3.24). For that k, when the constraint (3.24) is active, one can observe the growth of λ1,2 (Pk ) computed for the finite values of a in comparison with the eigenvalues of Pk as a → +∞. As it follows from the table data, these distinctions become more apparent as a becomes less. At Fig. 4 one can find the diagrams of the time-varying output-feedback gain coefficients Kk computed for a = 500, a = 1000, and a → +∞ with γ = 32. The diagrams of the gains Kk for the finite values of a practically coincide with the diagrams of Kk computed as a → +∞ for that points of the discrete time k, when the constraint (3.24) is not active (see Fig. 3)). Activation of this constraint gives rise to boundedness of the feedback gain amplitude. From the table data, where the maximum values of |Kk | over k = 0, . . . , N are given, one can see that the maximum amplitude of Kk is less as the disturbance anisotropy level a is less. In comparison with the ∞-controller, the AUTOMATION AND REMOTE CONTROL

Vol. 78

No. 7

2017

1216

TCHAIKOVSKY, TIMIN

Fig. 6. The phase trajectories x1,k , x2,k of the systems closed by the control laws computed for different anisotropy levels: (a) a = 500, (b) a = 1000, for γ = 32.

less values of the gains Kk of the anisotropic controllers lead to the less magnitudes of the control signals uk = Kk yk in the systems closed by these control laws, see Fig. 5 and the table. The trajectories of the closed-loop systems with the anisotropic controllers computed for a = 500, a = 1000, a → +∞ and γ = 32 represented at Fig. 6 are computed for feeding the system input by the random disturbance signals presented at Fig. 1. One can see that the magnitudes of oscillations in the closed-loop systems do not increase in time as distinct from the trajectories of the open-loop system without control (Fig. 2b). The maximum values of these oscillations are almost 40 times lesser than the same in the open-loop system. By the table data, in the closed-loop systems with the anisotropic controllers the magnitudes of the phase trajectories are less in comparison with the system closed by the ∞-controller computed as a → +∞. This fact is also confirmed by the diagrams of Fig. 6.

5. CONCLUSION In this paper, we derived the sufficient existence conditions for the anisotropic controller with time-varying parameters that guarantees the anisotropic norm of the closed-loop linear discrete time-varying system to be bounded by some given threshold value on finite time horizon. These conditions define the algorithm for computing the matrix of the controller parameters from recursive solving of the systems of the matrix inequalities or the sequence of the optimization problems. The control goal is to decrease the influence of the correlated random disturbances with unknown distribution law on the closed-loop system. The statistical uncertainty of the disturbance is measured by its anisotropy, which is the entropy-based measure of deviation of the disturbance distribution from the distribution of the nominal Gaussian white noise with scalar covariance matrix. The sensitivity of the closed-loop system is characterized by the anisotropic norm defined as the maximum mean-square gain of the system over the set of all disturbances with the anisotropy bounded by some given value. The ∞-controller that guarantees boundedness of the 2 -induced norm of the closed-loop linear time-varying system is the limiting case of the anisotropic controller. The procedure for computing the parameters of such controller comes to solving the system of linear matrix inequalities recursively. AUTOMATION AND REMOTE CONTROL

Vol. 78

No. 7

2017

SYNTHESIS OF ANISOTROPIC SUBOPTIMAL CONTROL

1217

The developed algorithm of synthesis allows to design the fixed-order, static state- and outputfeedback, fixed-structure, and decentralized controllers. Some drawbacks of the algorithm consist in its redundancy w.r.t. the given threshold value bounding the anisotropic norm of the closed-loop system, impossibility to minimize this threshold value, as well as certain computational complexity conditioned by need for recursive solving the system of matrix inequalities or convex optimization problems for each points of the discrete time on finite horizon. ACKNOWLEDGMENTS This work is partially supported by Program 15 of Division of Power Engineering, MachineBuilding, Mechanics and Control Processes of Russian Academy of Sciences, and Russian Foundation for Basic Research, project no. 17-08-00185. REFERENCES 1. Kwakernaak, H. and Sivan, R., Linear Optimal Control Systems, New York: Wiley, 1972. Translated under the title Lineinye optimal’nye sistemy upravleniya, Moscow: Mir, 1977. 2. Green, M. and Limebeer, D.J.N., Linear Robust Control, Englewood Cliffs: Prentice Hall, 1995. 3. Vladimirov, I.G., Diamond, F., and Kloeden, P., Anisotropy-Based Robust Performance Analysis of Finite Horizon Linear Discrete Time Varying Systems, Autom. Remote Control , 2006, vol. 67, no. 8, pp. 1265–1282. 4. Vladimirov, I.G., Kurdjukov, A.P., and Semyonov, A.V., Anisotropy of Signals and Entropy of Linear Stationary Systems, Dokl. Math., 1995, vol. 74, no. 3, pp. 388–390. 5. Diamond, P., Vladimirov, I.G., Kurdyukov, A.P., and Semyonov, A.V., Anisotropy-Based Performance Analysis of Linear Discrete Time Invariant Control Systems, Int. J. Control , 2001, no. 74, pp. 28–42. 6. Vladimirov, I.G., Kurdjukov, A.P., and Semyonov, A.V., State-Space Solution to Anisotropy-Based Stochastic H∞ -Optimization Problem, Proc. 13th IFAC World Congr., San-Francisco, USA, 1996, pp. 427–432. 7. Tchaikovsky, M.M. and Kurdyukov, A.P., Strict Anisotropic Norm Bounded Real Lemma in Terms of Matrix Inequalities, Dokl. Math., 2011, vol. 84, no. 3, pp. 895–898. 8. Maximov, E.A., Kurdyukov, A.P., and Vladimirov, I.G., Anisotropic Norm Bounded Real Lemma for Linear Discrete Time Varying Systems, Proc. 18th IFAC World Congr., Milano, Italy, 2011, pp. 4701–4706. 9. Yaesh, I. and Stoica, A.-M., Linear Time-Varying Anisotropic Filtering and its Application to Nonlinear Systems State Estimation, Proc. 2014 European Control Conf., Strasbourg, France, 2014, pp. 975–980. 10. Cover, T.M. and Thomas, J.A., Elements of Information Theory, New York: Wiley, 1991. 11. Freiling, G. and Ionescu, V., Time-Varying Discrete Riccati Equation: Some Monotonicity Results, Lin. Algebra Appl., 1999, vol. 286 (1), pp. 135–148. 12. Freiling, G. and Ionescu, V., Monotonicity and Convexity Properties of Matrix Riccati Equations, IMA J. Math. Control & Inf., 2001, vol. 18, pp. 61–72. 13. Bernstein, D.S., Matrix Mathematix: Theory, Facts, and Formulas with Application to Linear Systems Theory, New Jersey: Princeton Univ. Press, 2005. 14. Ben-Tal, A. and Nemirovskii, A., Lectures on Modern Convex Optimization, Haifa: Technion, 2000. 15. L¨ofberg, J., YALMIP: A Toolbox for Modeling and Optimization in Matlab, Proc. CACSD Conf., Taipei, Taiwan, 2004.

This paper was recommended for publication by A.P. Kurdyukov, a member of the Editorial Board AUTOMATION AND REMOTE CONTROL

Vol. 78

No. 7

2017