# On the Solution of a Nonlinear Semidefinite Program Arising in Discrete-Time Feedback Control Design

Authors: El-Sayed M. E. Mostafa
Journal: Journal of Applied Mathematics
 (2014)

---

## Abstract

A sequential quadratic programming method with line search is analyzed and studied for finding the local solution of a nonlinear semidefinite programming problem resulting from the discrete-time output feedback problem. The method requires an initial feasible point with respect to two positive definite constraints. By parameterizing the optimization problem we ease that requirement. The method is tested numerically on several test problems chosen from the benchmark
collection (Leibfritz, 2004).


---

## Body

## 1. Introduction

In this paper, the following nonlinear semidefinite programming (NSDP) problem is considered:(1)
min
⁡
J
(
X
)
s
.
t
.
H
(
X
)
=
0
,
G
(
X
)
≻
0
,
where J
:
ℝ
n
×
n
×
ℝ
p
×
r
→
ℝ, H
:
ℝ
n
×
n
×
ℝ
p
×
r
→
ℝ
n
×
n, G
:
ℝ
n
×
n
×
ℝ
p
×
r
→
ℝ
n
×
n are sufficiently smooth matrix functions, where G
≻
0 means that G is positive definite. This problem is assumed to be nonlinear and generally nonconvex. The origin of the considered NSDP problem is the static output feedback (SOF) design problem in which the unknown X is partitioned as X
=
(
L
,
F
), where L
∈
ℝ
n
×
n represents the state variable and F
∈
ℝ
p
×
r represents the control variable.In the last decade NSDP has attracted the attention of many authors in the optimization community. For instance, Jarre [1] introduced an interior-point method for nonconvex semi-definite programs. Leibfritz and Mostafa [2] proposed an interior-point trust region method for a special class of NSDP problems resulting from the continuous-time SOF problem. Kočvara et al. [3] considered an augmented Lagrangian method for a similar NSDP problem. Sun et al. [4] investigated the rate of convergence of the augmented Lagrangian approach for NSDP. Correa and Ramirez [5] proposed sequential semi-definite programming method for nonlinear semi-definite programming. Yamashita and Yabe [6] study local and superlinear convergence of a primal-dual interior point method. Freund et al. [7] proposed a sequential semidefinite programming approach for a nonlinear program with nonlinear semidefinite constraints.The design problem of optimal output feedback controllers is one of the most studied problems over the last four decades in the community of systems and control. Clearly, this is due to the existence of numerous practical applications in system and control engineering, finance, and statistics. The benchmark collection [8], for instance, contains over 160 applications from system and control engineering only. Mäkilä and Toivonen in the survey [9] summarized several special purpose algorithms such as the Levine-Athans method, the dual Levine-Athans method, and the descent Anderson-Moore method as well as Newton’s method for finding the local solution of a matrix optimization problem that corresponds to the discrete-time SOF design problem. Rautert and Sachs [10] proposed a structured quasi-Newton method for the continuous-time SOF problem. Mäkilä [11] introduced a matrix optimization problem related to the discrete-time SOF design problem, which is an extension of the former problem considered in [9]. However, it has not yet been tested numerically.Several constrained optimization techniques have been developed for constrained optimization problems related to the continuous-time SOF design problem. For instance, Apkarian et al. [12] introduced a linear matrix inequality approach for the H
2 synthesis problem. Leibfritz and Mostafa [13] proposed two trust region methods for two different formulations of the continuous-time SOF problem. Mostafa [14, 15] developed a trust region method for the decentralized SOF problem and an augmented Lagrangian SQP method for the continuous-time SOF problem, respectively. More recently, Mostafa [16] introduced a trust region method for a particular constrained optimization problem that originates from the discrete-time SOF problem.In the current work a different formulation of the design problem is considered that underlies the NSDP problem formulation  (1). In fact such an NSDP problem can be seen as a generalization of the matrix optimization problem introduced by Mäkilä [11]. The goal is to study and analyze an SQP method with line search for finding a local solution of such an NSDP problem. The reduced Hessian with a step decomposition approach (see, e.g., [17]) is employed with the SQP method. Moreover, we compare the SQP method numerically versus Newton’s method proposed by Mäkilä and Toivonen [9] but applied on the matrix optimization problem [11]. As is well known, SQP methods are among the most effective techniques for nonlinear constrained optimization problems; see, for example, [17]. In particular, the problem structure is exploited with the considered method.The existence of an initial feasible pointX
0
=
(
L
0
,
F
0
) with respect to the positive definite constraints is necessary to start the iteration sequence of the SQP method. Such an initial feasible point can be easy to find in simple problems but can be difficult or very time consuming for moderate and complex problems. A parametrization approach is employed on the NSDP problem that overcomes this difficulty.This paper is organized as follows. In the next section the discrete-time SOF problem is introduced and the related NSDP problem is stated. In Section3 the SQP method that finds a local solution of the NSDP is presented. Section 4 is devoted to overcome the difficulty of finding an initial feasible point (
L
0
,
F
0
) to start the iteration sequence of the SQP method. In order to demonstrate the considered approach we have tested in Section 5 the SQP method on various test problems derived from benchmark collection [8].Notations. For a matrix G
∈
ℝ
n
×
n the notation G
≻
0 (G
≽
0) denotes that G is strictly positive (semi)definite. The symbol ⊗ denotes the Kronecker product of two matrices and 〈
·
,
·
〉 denotes the inner product operator. Throughout the paper, ∥
·
∥ is the Frobenius norm; otherwise; it will be mentioned. Moreover, I
n denotes the identity matrix of order n. We omit the argument when it is known from the context; for example, we use J
k to denote J
(
X
k
).


## 2. The Static Output Feedback Control Problem

The problem of designing an optimal output feedback controller for discrete-time control systems is an important and extensively studied topic in system and control literature; see, for example, the two survey papers by Mäkilä and Toivonen [9] and Syrmos et al. [18] and numerous references therein, including Garcia et al. [19], Lee and Khargonekar [20], Mäkilä [11], Mostafa [21, 22], and Sulikowski et al. [23].Consider the linear time-invariant control system described by the following discrete time state-space realization:(2)
x
k
+
1
=
A
x
k
+
B
u
k
+
v
k
,
y
k
=
C
x
k
+
e
k
,
k
≥
0
,
where x
k
∈
ℝ
n is the state vector, u
k
∈
ℝ
p is the control (input) vector, y
k
∈
ℝ
r is the measured output vector, and v
k
∈
ℝ
n and e
k
∈
ℝ
r are disturbance (noise) terms. Moreover A, B, and C are given constant matrices of appropriate dimensions.We will consider the following quadratic cost function as an objective function to be minimized:(3)
J
=
𝔼
[
∑
k
=
0
∞
〈
x
k
,
Q
x
k
〉
+
〈
u
k
,
R
u
k
〉
‍
]
,
where 𝔼
[
·
] is the expected value and Q
∈
ℝ
n
×
n and R
∈
ℝ
p
×
p are given constant weight matrices.The following control law is often used to close the above control system:(4)
u
k
=
F
y
k
,
where F
∈
ℝ
p
×
r is the output feedback gain matrix representing the unknown.In order to remove the dependency of the problem on the initial state vectorx
0 the following assumption is often imposed.Assumption 1.
Assume thatx
0 is a random variable uniformly distributed over the unit ball such that 𝔼
[
x
0
]
=
0.
Under this assumption the optimal control problem  (2)–(4) can be transformed into the matrix optimization problem (see [11]):
(5)
min
⁡
F
J
^
(
F
)
≔
J
(
L
(
F
)
,
F
)
=
Tr
⁡
(
L
(
F
)
(
Q
+
C
T
F
T
R
F
C
)
)
+
Tr
⁡
(
F
T
R
F
R
e
)
,
where Tr
⁡
(
·
) is the trace operator and L
(
F
) solves the discrete Lyapunov equation:
(6)
L
(
F
)
=
(
A
+
B
F
C
)
L
(
F
)
(
A
+
B
F
C
)
T
+
V
+
B
F
R
e
F
T
B
T
.Assumption 2.
Throughout the paper we assume thatA
∈
ℝ
n
×
n, B
∈
ℝ
n
×
p, and C
∈
ℝ
r
×
n are given constant matrices. Furthermore, we assume that Q
∈
ℝ
n
×
n and R
∈
ℝ
p
×
p are given constant weight matrices and R
e
∈
ℝ
r
×
r and V
∈
ℝ
n
×
n are given covariance matrices. We assume that Q
≽
0 while R, R
e and V
≻
0.From Lyapunov stability theory the unknownF must be chosen from the following set:
(7)
𝒮
F
=
{
F
∈
ℝ
p
×
r
:
ρ
(
A
+
B
F
C
)
<
1
}
,
where ρ
(
·
) is the spectral radius. This restriction is required ultimately to guarantee the existence of a unique solution of the discrete Lyapunov equation  (6); see, for example, [24]. Moreover, from the Lyapunov stability theory the next lemma gives an equivalence between the asymptotic stability of the system  (2) and fulfilling two positive definite constraints.Lemma 3 ([24, Lemma 5.7
.
19]).
Consider the SOF problem  (2)–(4). An F
∈
𝒮
F exists if and only if there exists (
L
,
F
)
∈
𝒟
s, where
(8)
𝒟
s
=
{
(
L
,
F
)
∈
ℝ
n
×
n
×
ℝ
p
×
r
:
L
≻
0
,
Y
(
L
,
F
)
≻
0
}
,
where
(9)
Y
(
L
,
F
)
=
L
-
(
A
+
B
F
C
)
L
(
A
+
B
F
C
)
T
.According to Lemma3 and by considering both of L and F as independent variables, where we skip the imposed dependency of L on F that was considered in  (5)-(6), then the minimization problem can be stated as the following NSDP:
(10)
min
⁡
L
,
F
J
(
L
,
F
)
=
Tr
⁡
(
L
(
Q
+
C
T
F
T
R
F
C
)
)
+
Tr
⁡
(
F
T
R
F
R
e
)
,
(11)
s
.
t
.
H
(
L
,
F
)
=
-
L
+
(
A
+
B
F
C
)
L
(
A
+
B
F
C
)
T
+
V
+
B
F
R
e
F
T
B
T
=
0
,
(12)
Y
(
L
,
F
)
≻
0
,
L
≻
0
.Instead of considering  (5)–(7) we would rather consider the constrained minimization problem  (10)–(12) which can be seen as a generalization of the former problem. The equivalence between the stability condition  (7) and the positive definite constraints of  (8) provides a straightforward idea to fulfill such a constraint within the numerical algorithm explicitly. There are, however, different alternatives to handle the positive definite constraints within a numerical algorithm. First, one might incorporate those positive definite terms within the objective function and apply an interior-point technique on the corresponding constrained problem; see [2]. We might also use slack matrix variables for the inequality constraints. Then one solves the problem with only equality constraints.We end this section by noting that the two sets𝒮
F and 𝒟
s are open. Therefore, it is convenient to replace them by the two subsets Ω
^ and Ω, respectively, where both subsets are assumed to be compact. For example, the level set
(13)
Ω
^
=
{
F
∈
𝒮
F
:
J
^
(
F
)
≤
J
^
(
F
0
)
}
⊆
𝒮
F
is compact; see [9].


## 3. An SQP Method for the Design Problem

Let us consider the minimization problem  (10)–(12). The Lagrangian function associated with the equality constraint is defined as
(14)
ℓ
(
L
,
F
,
K
)
=
J
(
L
,
F
)
+
〈
K
,
H
(
L
,
F
)
〉
,
where K
∈
ℝ
n
×
n is the Lagrange multiplier. The next lemma includes first- and second-order derivatives of ℓ and first derivatives of the function H. Similar results can be found in [16, Lemma 2.1] but on a simpler problem. In the following discussion we denote A
(
F
)
=
A
+
B
F
C and Q
(
F
)
=
Q
+
C
T
F
T
R
F
C.Lemma 4.
Consider the NSDP problem  (10)–(12). The Lagrangian function ℓ
(
L
,
F
,
K
) is twice continuously differentiable and its first- and second-order directional derivatives as well as first-order derivatives of the constraint function H are the following:
(15)
ℓ
L
(
·
)
Δ
L
=
Tr
⁡
(
(
A
(
F
)
T
K
A
(
F
)
+
Q
(
F
)
-
K
)
Δ
L
)
,
ℓ
F
(
·
)
Δ
F
=
Tr
⁡
(
2
K
F
L
C
T
Δ
F
T
+
2
(
B
T
K
B
+
R
)
F
R
e
Δ
F
T
)
,
ℓ
L
L
(
·
)
(
Δ
L
,
Δ
L
)
=
0
,
ℓ
L
F
(
·
)
(
Δ
L
,
Δ
F
)
=
ℓ
F
L
(
L
,
F
,
K
)
(
Δ
F
,
Δ
L
)
=
Tr
⁡
(
2
K
F
Δ
L
C
T
Δ
F
T
)
,
ℓ
F
F
(
·
)
(
Δ
F
,
Δ
F
)
=
Tr
⁡
(
2
(
B
T
K
B

+

R
)
Δ
F
(
C
L
C
T

+

R
e
)
Δ
F
T
)
,
H
L
(
·
)
Δ
L
≔
〈
Δ
L
,
∇
L
H
(
·
)
〉
=
Δ
L
-
A
(
F
)
Δ
L
A
(
F
)
T
,
H
F
(
·
)
Δ
F
≔
〈
Δ
F
,
∇
F
H
(
·
)
〉
=
B
Δ
F
L
F
T
+
L
F
Δ
F
T
B
T
,
where K
F
=
B
T
K
A
(
F
)
+
R
F
C and L
F
=
A
(
F
)
L
C
T
+
B
F
R
e.Proof.
By using the directional derivatives with both ofℓ and H the above derivatives are obtained.The Karush-Kuhn-Tucker (KKT) optimality conditions yield a nonlinear system of matrix equations as given in the next lemma.Lemma 5 (see, e.g., [7, Section 5]).
Let(
L
,
F
)
∈
𝒟
s be a local minimizer of the problem  (10)–(12) and let  (10)–(12) be regular at (
L
,
F
). Then there exists a Lagrange multiplier K
∈
ℝ
n
×
n satisfying with L and F the following KKT conditions:
(16)
ℳ
(
U
)
=
(
2
K
F
L
C
T
+
2
(
B
T
K
B
+
R
)
F
R
e
-
L
+
A
(
F
)
L
A
(
F
)
T
+
V
+
B
F
R
e
F
T
B
T
-
K
+
A
(
F
)
T
K
A
(
F
)
+
Q
(
F
)
)
=
0
,
where K
F is as defined in Lemma 4 and U
=
(
L
,
F
,
K
)
∈
𝒟
s
×
ℝ
n
×
n.Clearly, the system  (16) is nonlinear in the unknowns L, F, and K. However, the last two equations of  (16) might be viewed as discrete Lyapunov equations of the following form
(17)
L
=
A
(
F
)
L
A
(
F
)
T
+
V
+
B
F
R
e
F
T
B
T
,
(18)
K
=
A
(
F
)
T
K
A
(
F
)
+
Q
(
F
)
.The SQP method is based on minimizing successfully a quadratic program of the form:(19)
min
⁡
(
Δ
L
,
Δ
F
)
q
(
Δ
L
,
Δ
F
)
s
.
t
.
h
(
Δ
L
,
Δ
F
)
=
0
,
where
(20)
q
(
Δ
L
,
Δ
F
)
=
ℓ
L
Δ
L
+
ℓ
F
Δ
F
+
1
2
(
ℓ
L
L
(
Δ
L
,
Δ
L
)
+
2
ℓ
L
F
(
Δ
L
,
Δ
F
)
+
ℓ
F
F
(
Δ
F
,
Δ
F
)
)
,
h
(
Δ
L
,
Δ
F
)
=
H
L
(
L
,
F
)
Δ
L
+
H
F
(
L
,
F
)
Δ
F
+
H
(
L
,
F
)
.
All derivatives required for constructing this QP are given in Lemma 4.An SQP method known also as Newton-Lagrange method generates a sequence of iterates by solving the augmented linear system:(21)
(
∇
L
L
2
ℓ
k
∇
L
F
2
ℓ
k
∇
L
H
k
*
∇
F
L
2
ℓ
k
∇
F
F
2
ℓ
k
∇
F
H
k
*
∇
L
H
k
∇
F
H
k
0
)
(
Δ
L
Δ
F
Δ
K
)
=
-
(
∇
L
ℓ
k
∇
F
ℓ
k
H
k
)
and updates the step by
(22)
(
L
k
+
1
,
F
k
+
1
,
K
k
+
1
)
=
(
L
k
+
α
k
Δ
L
k
,
F
k
+
α
k
Δ
F
k
,
K
k
+
α
k
Δ
K
k
)
,
where α
∈
(
0,1
] is the step size.Assumption 6.
Throughout this paper we make the following assumptions.(H1)
The mappingsJ and H are at least twice continuously differentiable.
(H2)
For a given(
L
,
F
)
∈
𝒟
s the mapping ∇
L
H
(
L
,
F
) is invertible.
(H3)
There exists an(
L
0
,
F
0
)
∈
𝒟
s. Moreover, assume that for a given (
L
k
,
F
k
)
∈
𝒟
s there always exists an α
k
>
0 such that (
L
k
+
α
k
Δ
L
k
,
F
k
+
α
k
Δ
F
k
)
∈
𝒟
s.The next theorem provides us with the linearization of the nonlinear system  (16).Theorem 7.
For a givenU
:
=
(
L
,
F
,
K
)
∈
𝒟
s
×
ℝ
n
×
n the linearized equations of the KKT system (16) are the following:
(23)
(
(
B
T
K
B
+
R
)
Δ
F
(
C
L
C
T
+
R
e
)
+
K
F
Δ
L
C
T
+
B
T
Δ
K
L
F
B
Δ
F
L
F
T
+
L
F
Δ
F
T
B
T
-
Δ
L
+
A
(
F
)
Δ
L
A
(
F
)
T
C
T
Δ
F
T
K
F
+
K
F
T
Δ
F
C
-
Δ
K
+
A
(
F
)
T
Δ
K
A
(
F
)
)
=
-
(
K
F
L
C
T
+
(
B
T
K
B
+
R
)
F
R
e
-
L
+
A
(
F
)
L
A
(
F
)
T
+
V
+
B
F
R
e
F
T
B
T
-
K
+
A
(
F
)
T
K
A
(
F
)
+
Q
(
F
)
)
,
where K
F and L
F are as defined in Lemma 4.Proof.
The linearization of the nonlinear systemℳ
(
U
)
=
0 about the point U
=
(
L
,
F
,
K
) in the direction of Δ
U
=
(
Δ
L
,
Δ
F
,
Δ
K
) takes the form
(24)
∇
ℳ
(
U
)
T
Δ
U
=
-
ℳ
(
U
)
,
where ∇
ℳ
(
U
)
T is the Jacobian matrix, with which we directly obtain the system  (23).Under certain assumptions it is possible to show that the SQP directionΔ
U given by  (23) is well defined. The difficulty that one faces when solving  (23) is due to the fact that we do not have an explicit representation of the Jacobian of ℳ
(
U
). However, if we write U
=
(
L
,
F
,
K
) as a long column vector vec
⁡
(
U
)
∈
ℝ
p
r
+
2
n
2 such that f
(
vec
⁡
(
U
)
)
=
ℳ
(
U
), then we are able to evaluate the Jacobian of f using the basic properties of the Kronecker product which is an (
2
n
2
+
p
r
)
×
(
2
n
2
+
p
r
) matrix.Theorem 8.
LetU
=
(
L
,
F
,
K
)
∈
𝒟
s
×
ℝ
n
×
n be given. The Jacobian matrix of ℳ
(
U
) is given by:
(25)
∇
ℳ
(
U
)
T
=
(
(
C
L
C
T
+
R
e
)
⊗
(
B
T
K
B
+
R
)
C
⊗
K
F
L
F
T
⊗
B
T
[
L
F
⊗
B
+
K
n
n
(
L
F
⊗
B
)
]
S
F
0
[
K
r
n
(
C
T

⊗

K
F
T
)

+

(
C
T

⊗

K
F
T
)
]
0
S
F
T
)
,
where K
F and L
F are as defined in Lemma 4, S
F
=
-
I
n
2
+
A
(
F
)
⊗
A
(
F
), and K
m
n is an m
n
×
m
n commutation matrix defined as vec
⁡
(
W
T
)
=
K
r
n
vec
⁡
(
W
) for any m
×
n matrix W. Furthermore, the linear system (23) is explicitly written as
(26)
∇
ℳ
(
U
)
T
vec
⁡
(
U
)
=
-
vec
⁡
(
ℳ
)
,
where
(27)
vec
⁡
(
ℳ
)
=
(
vec
⁡
(
K
F
L
C
T
+
(
B
T
K
B
+
R
)
F
R
e
)
vec
⁡
(
-
L
+
A
(
F
)
L
A
(
F
)
T
+
V
+
B
F
R
e
F
T
B
T
)
vec
⁡
(
-
K
+
A
(
F
)
T
K
A
(
F
)
+
Q
(
F
)
)
)
.Proof.
For any given three matricesA, B, and C of appropriate dimensions the linear matrix equation
(28)
A
X
B
=
C
is equivalent to the explicit linear system
(29)
(
B
T
⊗
A
)
vec
⁡
(
X
)
=
vec
⁡
(
C
)
.
If we apply this result on the above linear system the entries of the Jacobian matrix can be easily obtained.Indeed, the last result shows that the dimension of the Jacobian matrix  (25) can be very large because of the Kronecker products. Moreover, the evaluation of the Jacobian matrix explicitly in the SQP method is unnecessary and makes the method impractical. By taking this fact in mind the reduced Hessian approach is rather considered as an ideal candidate for finding the local solution of the NSDP  (10)–(12); see, for example, [2, 16].


### 3.1. The Reduced Hessian Approach

The reduced Hessian approach relies on the existence of the null space and the right inverse of the Jacobian∇
H
T. In optimal control such operators are not available explicitly. However, it is possible to construct them from the problem itself; for more details see, for example, [2, 16, 25]. The two operators are defined as
(30)
Z
(
L
,
F
)
=
(
-
H
L
-
1
(
L
,
F
)
H
F
(
L
,
F
)
,
I
p
)
,
R
^
(
L
,
F
)
=
(
-
H
L
-
1
(
L
,
F
)
,
0
)
,
where I
p and 0 are the identity and the zero operators, respectively. We note that if the Jacobian ∇
H
T is surjective then such a right inverse R
^ exists. We also note that the range space of Z coincides with the null space of the Jacobian ∇
H
T; see [16, Lemma 2].Next, with the help ofZ and R
^ the following result shows that the linearized equality constraint is decomposed into two discrete Lyapunov equations.Lemma 9.
Let(
L
,
F
)
∈
𝒟
s be given and let (
Δ
L
,
Δ
F
)
∈
ℝ
n
×
n
+
p
×
r be the solution of the QP  (19). Then the linearized equality constraint of  (19) is decomposed into the following discrete Lyapunov equations:
(31)
Δ
L
n
=
A
(
F
)
Δ
L
n
(
A
F
)
T
+
H
(
L
,
F
)
,
(32)
Δ
L
t
=
A
(
F
)
Δ
L
t
A
(
F
)
T
+
B
Δ
F
L
F
T
+
L
F
Δ
F
T
B
T
,
where L
F is as given in Lemma 5.Proof.
(see also [16, Lemma 3]). The linearized equality constraint of (19) can be rewritten as
(33)
∇
L
H
(
L
,
F
)
Δ
L
t
+
∇
F
H
(
L
,
F
)
Δ
F
+
∇
L
H
(
L
,
F
)
Δ
L
n
+
H
(
L
,
F
)
=
0
,
where Δ
L
=
Δ
L
n
+
Δ
L
t. According to the definition of Z and R
^ the above equation can be decomposed as
(34)
∇
L
H
(
L
,
F
)
Δ
L
n
+
H
(
L
,
F
)
=
0
∀
Δ
L
n
,
H
L
(
L
,
F
)
Δ
L
t
(
Δ
F
)
+
H
F
(
L
,
F
)
Δ
F
=
0
∀
Δ
F
,
which yield the discrete Lyapunov equations  (31) and (32), respectively.Following the proof of Lemma9 the step (
Δ
L
,
Δ
F
) that solves the QP (19) with the help of Z and R
^ is decomposed as
(35)
(
Δ
L
,
Δ
F
)
=
Z
(
L
,
F
)
Δ
F
+
R
^
(
L
,
F
)
H
(
L
,
F
)
=
:
(
Δ
L
t
(
Δ
F
)
,
Δ
F
)
+
(
Δ
L
n
,
0
)
,
where (
Δ
L
t
(
Δ
F
)
,
Δ
F
) and (
Δ
L
n
,
0
) are the tangential and the normal steps, respectively.Consequently, one can obtain the normal step(
Δ
L
n
,
0
) by solving the discrete Lyapunov equation  (31) for Δ
L
n. Having obtained such a step component we solve the following tangential subproblem for Δ
F:
(36)
min
⁡
Δ
F
∈
ℝ
p
×
r
ψ
(
Δ
F
)
:
=
q
(
Δ
L
t
(
Δ
F
)
,
Δ
F
)
=
〈
Δ
L
t
(
Δ
F
)
,
∇
L
ℓ
(
·
)
〉
+
〈
Δ
F
,
∇
F
ℓ
(
·
)
〉
+
〈
Δ
F
,
∇
F
L
2
ℓ
(
·
)
Δ
L
n
〉
+
〈
Δ
L
t
(
Δ
F
)
,
∇
L
L
2
ℓ
(
·
)
Δ
L
n
〉
+
1
2
〈
Δ
F
,
Z
T
∇
2
ℓ
(
·
)
Z
Δ
F
〉
,
where
(37)
〈
Δ
F
,
Z
T
∇
2
ℓ
(
·
)
Z
Δ
F
〉
=
〈
Δ
L
t
(
Δ
F
)
,
∇
L
L
2
ℓ
(
·
)
Δ
L
t
(
Δ
F
)
〉
+
〈
Δ
F
,
∇
F
F
2
ℓ
(
·
)
Δ
F
〉
+
2
〈
Δ
F
,
∇
F
L
2
ℓ
(
·
)
Δ
L
t
(
Δ
F
)
〉
and Δ
L
t
(
Δ
F
) is the corresponding solution of the discrete Lyapunov equation  (32).In order to find the local solution of the minimization problem  (36) let us apply the necessary optimality conditions, which give the following result.Lemma 10.
Let(
L
,
F
,
K
)
∈
𝒟
s
×
ℝ
n
×
n and Δ
L
n
∈
ℝ
n
×
n be given. Moreover, let Δ
F
∈
ℝ
p
×
r be a local minimum for the problem (36). Then Δ
F solves the following linear matrix equation:
(38)
(
B
T
K
B
+
R
)
Δ
F
(
C
L
C
T
+
R
e
)
+
K
F
Δ
L
t
(
Δ
F
)
C
T
+
B
T
Δ
K
(
Δ
F
)
L
F
=
-
K
F
(
L
+
Δ
L
n
)
C
T
-
(
B
T
K
B
+
R
)
F
R
e
-
B
T
S
L
F
,
where Δ
L
t
(
Δ
F
), Δ
K
(
Δ
F
), and S, respectively, solve the discrete Lyapunov equations  (32) and
(39)
Δ
K
=
A
(
F
)
T
Δ
K
A
(
F
)
+
K
F
T
Δ
F
C
+
C
T
Δ
F
T
K
F
,
(40)
S
=
A
(
F
)
T
S
A
(
F
)
+
M
+
M
T
,
where K
F and L
F are as defined in Lemma 4 and
(41)
M
=
-
K
+
A
(
F
)
T
K
A
(
F
)
+
Q
(
F
)
.Proof.
We use the derivatives ofℓ given in Lemma 4 to construct ψ
(
Δ
F
). If we differentiate ψ with respect to Δ
F and apply the necessary optimality conditions, we obtain
(42)
Z
T
∇
2
ℓ
(
·
)
Z
Δ
F
=
-
∇
F
ℓ
(
·
)
-
∇
F
L
2
Δ
L
n
-
∂
∂
Δ
F
〈
Δ
L
t
(
Δ
F
)
,
∇
L
ℓ
(
·
)
〉
,
which is another representation of  (38). As an example of evaluating these derivatives we obtain the last term of  (42) as follows:
(43)
∂
∂
Δ
F
〈
Δ
L
t
(
Δ
F
)
,
∇
L
ℓ
(
·
)
〉
=
∂
∂
Δ
F
Tr
⁡
(
M
Δ
L
t
(
Δ
F
)
)
=
∂
∂
Δ
F
Tr
⁡
(
M
∑
k
=
0
∞
A
(
F
)
k
(
B
Δ
F
L
F
T
∑
k
=
0
∞
A
(
F
)
k
(
B
Δ
F
L
F
T
+
L
F
Δ
F
T
B
T
)
(
A
(
F
)
T
)
k
)
=
B
T
(
∑
k
=
0
∞
(
A
(
F
)
T
)
k
(
M
+
M
T
)
A
(
F
)
k
‍
)
L
F
=
B
T
S
L
F
,
where
(44)
Δ
L
t
(
Δ
F
)
=
∑
k
=
0
∞
‍
A
(
F
)
k
(
B
Δ
F
L
F
T
+
L
F
Δ
F
T
B
T
)
(
A
(
F
)
T
)
k
,
S
=
∑
k
=
0
∞
‍
(
A
(
F
)
T
)
k
(
M
+
M
T
)
A
(
F
)
k
are the exact solutions of the discrete Lyapunov equations  (32) and  (40), respectively.Because of using directional derivatives the resulting linear system  (38)–(40) is not given explicitly. However, one can modify the conjugate gradient (CG) method [17] and solve this system iteratively. The problem structure is exploited, where the partition of the unknown X
=
(
L
,
F
) is utilized.For a givenΔ
F calculated by the CG method both of Δ
L
t
(
Δ
F
) and Δ
K
(
Δ
F
) are uniquely determined by solving the discrete Lyapunov equations  (32) and  (39), respectively. The CG algorithm for solving the system  (38)–(40) is stated in the following lines, where we use T
i, D
i, and E
i at the ith CG inner iteration to denote the approximation of the step component Δ
F, the CG direction, and the residual, respectively.Algorithm 11 (modified CG method).
LetU
=
(
L
,
F
,
K
)
∈
𝒟
s
×
ℝ
n
×
n and Δ
L
n be given. Moreover, let A, B, C, Q, R, R
e, and V be given constant matrices. Set T
0
=
0, E
0
=
-
K
F
L
C
T
-
(
B
T
K
B
+
R
)
F
R
e, and D
0
=
E
0. Choose η and ϵ
∈
(
0,1
).
Fori
=
0,1
,
…, do(1)
Solve forΔ
L
t
(
D
i
) and Δ
K
(
D
i
) the discrete Lyapunov equations  (32) and  (39), respectively.
(2)
Calculate(45)
N
(
D
i
)
≔
(
B
T
K
B
+
R
)
D
i
(
C
L
C
T
+
R
e
)
+
K
F
Δ
L
t
(
D
i
)
C
T
+
B
T
Δ
K
(
D
i
)
L
F
.
(3)
IfTr
⁡
(
D
i
T
N
(
D
i
)
)
≤
ϵ
∥
D
i
∥
2, update the step Δ
F as follows: if i
=
0, set Δ
F
=
D
0; else, set Δ
F
=
T
i and exit.
(4)
Calculate the ratioθ
i
=
Tr
⁡
(
E
i
T
E
i
)
/
Tr
⁡
(
D
i
T
N
(
D
i
)
). Set T
i
+
1
=
T
i
+
θ
i
D
i and E
i
+
1
=
E
i
-
θ
i
N
(
D
i
).
(5)
If∥
E
i
+
1
∥
≤
η
∥
E
0
∥ set Δ
F
=
T
i
+
1 and exit.
(6)
Calculate the ratio(46)
β
i
=
Tr
⁡
(
E
i
+
1
T
E
i
+
1
)
Tr
⁡
(
E
i
T
E
i
)
.
SetD
i
+
1
=
E
i
+
1
+
β
i
D
i, i
:
=
i
+
1, and go to (
1
).
End (do)There are two exits in this CG method. The first one takes place when a negative or zero curvature is encountered by the calculated search direction. The second exit occurs when the norm of the new residual is sufficiently small.Having evaluated both ofΔ
L
k
n and (
Δ
L
t
(
Δ
F
k
)
,
Δ
F
k
), then according to the step decomposition  (35) the total step (
Δ
L
k
,
Δ
F
k
) is obtained. Consequently, the new iterate takes the form
(47)
(
L
k
+
1
,
F
k
+
1
)
=
(
L
k
+
α
k
Δ
L
k
,
F
k
+
α
k
Δ
F
k
)
,
where α
k
>
0 is the step size. Starting from α
0
=
1, we calculate the step size α
k by a backtracking rule (e.g., by halving α) until we satisfy the following sufficient decrease condition:
(48)
Φ
μ
k
(
L
k
+
α
k
Δ
L
,
F
k
+
α
k
Δ
F
)
≤
Φ
k
μ
k
+
τ
α
k
[
Φ
L
μ
k
(
L
k
,
F
k
)
Δ
L
+
Φ
F
μ
k
(
L
k
,
F
k
)
Δ
F
]
,
where τ
∈
(
0,1
/
2
) is a given parameter. Let α
~
k be the achieved step size. Starting from that α
~
k we satisfy the following stability condition by backtracking in the same way until an α
k
>
0 is reached such that
(49)
L
k
+
α
k
Δ
L
k
≻
0
,
Y
(
L
k
+
α
k
Δ
L
k
,
F
k
+
α
k
Δ
F
k
)
≻
0
.We assume that there always exists anα
k
>
0 such that  (48) is satisfied. If α is decreased and reached a lower bound α
min
⁡
>
0 without fulfilling  (49), then we might follow one of the ideas often used in modifying the Hessian matrix of Newton-based methods; see, for example, [17, Section 6.3]. Let W
k denote any of the two matrices L
k
+
α
min
⁡
Δ
L
k or Y
(
L
k
+
α
min
⁡
Δ
L
k
,
F
k
+
α
min
⁡
Δ
F
k
) with a spectral decomposition W
k
=
Q
~
Λ
Q
~
T. A correction matrix Δ
W
k of minimum Frobenius norm that ensures λ
min
⁡
(
W
k
+
Δ
W
k
)
≥
δ
>
0 is given by
(50)
Δ
W
k
=
Q
~
diag
⁡
(
t
i
)
Q
~
T
,
where
(51)
t
i
=
{
0
,
λ
i
>
δ
,
δ
-
λ
i
,
otherwise
.
The modified matrix is
(52)
W
k
+
Δ
W
k
=
Q
~
(
Λ
+
diag
⁡
(
t
i
)
)
Q
~
T
,
where λ
min
⁡
(
W
k
) means the smallest eigenvalue of W
k. Note that the nature of the stability constraint is an inactive constraint for most test problems.The Matlab function[
~
,
j
]
=
cholinc
(
M
) of the incomplete Cholesky factorization can be used for checking the positive definiteness, where the indicator j is used for that purpose. On the other hand, as mentioned at the end of Section 2 we can check the stability condition via the spectral radius of the closed-loop system matrix; namely, ρ
(
A
(
F
k
+
α
k
Δ
F
k
)
)
<
1.We use as a merit function the followingl
1-penalty function for globalizing the SQP method:
(53)
Φ
μ
(
L
,
F
)
=
J
(
L
,
F
)
+
μ
∥
H
(
L
,
F
)
∥
1
,
where μ
>
0 is the penalty parameter. This function is directionally differentiable.The following update rule for the penalty parameterμ is considered (see, e.g., [17, page 547]):
(54)
μ
k
=
{
μ
k
-
1
,
if
μ
k
-
1
≥
∥
Δ
K
k
∥
∞
+
γ
,
∥
Δ
K
k
∥
∞
+
2
γ
,
otherwise
,
where γ
>
0 is a given constant and Δ
K
k is obtained by solving the discrete Lyapunov equation  (39).The following algorithm demonstrates the SQP method that finds the local solution of the NSDP  (10)–(12).Algorithm 12 (the method DSQP).
Let(
L
0
,
F
0
,
K
0
)
∈
𝒟
s
×
ℝ
n
×
n be given. Moreover, let A, B, C, Q, R, R
e, and V be given constant matrices. Choose τ
∈
(
0,1
/
2
), γ
>
0, and ϵ
tol
∈
(
0,1
).
While∥
∇
U
ℓ
k
∥
>
ϵ
tol, do(1)
CalculateΔ
L
k
n solution of the discrete Lyapunov equation  (31).
(2)
For a givenΔ
L
k
n calculate the local solution (
Δ
L
t
(
Δ
F
k
)
,
Δ
F
k
) of the problem (36) by Algorithm 11. Then set (
Δ
L
k
,
Δ
F
k
)
=
(
Δ
L
k
n
,
0
)
+
(
Δ
L
t
(
Δ
F
k
)
,
Δ
F
k
).
(3)
CalculateΔ
K
k by solving (39) and then update the penalty parameter μ
k using  (54).
(4)
Starting fromα
0
=
1 apply a backtracking step-size rule that satisfies the sufficient decrease conditions  (48).
(5)
If(
L
k
+
α
k
Δ
L
k
,
F
k
+
α
k
Δ
F
k
)
)
∉
𝒟
s, then maintain the positive definiteness of  (49) as described above.
(6)
SetF
k
+
1
=
F
k
+
α
k
Δ
F
k, L
k
+
1
=
L
k
+
α
k
Δ
L
k, K
k
+
1
=
K
k
+
α
k
Δ
K
k, and k
:
=
k
+
1; then go to step (
1
). End (While)Step(
5
) of Algorithm 12 can also be modified as follows. Set F
k
+
1
=
F
k
+
α
k
Δ
F
k. Then obtain L
k
+
1
=
L
(
F
k
+
1
) and K
k
+
1
=
K
(
F
k
+
1
) by solving the discrete Lyapunov equations  (17) and  (18), respectively. This means we project the search direction from the space 𝒟
s to the control-variable space ℝ
p
×
r.The global convergence behavior of the method DSQP is beyond the scope of this work. However, we state some of the main global convergence results, which can be shown under the following general assumptions.Assumption 13.
In addition to Assumption6 we further assume that the functions J and H together with their derivatives are in Ω
⊆
𝒟
s. We further assume that the operator H
L
-
1 is uniformly bounded in 𝒟
s; furthermore, we assume that the sequences {
∇
2
ℓ
k
}, {
L
k
}, {
F
k
}, {
K
k
}, and {
Z
(
L
k
,
F
k
)
} are all bounded in Ω.The next lemma provides an upper bound for the directional derivative of the merit function by the computed search direction.Lemma 14.
Let(
Δ
L
,
Δ
F
,
Δ
K
) be generated by the SQP iteration  (21). Then the directional derivative of Φ
μ
(
L
,
F
) in the direction (
Δ
L
,
Δ
F
) satisfies
(55)
Φ
L
μ
Δ
L
+
Φ
F
μ
Δ
F
≤
-
[
ℓ
L
L
(
Δ
L
,
Δ
L
)
+
2
ℓ
L
F
(
Δ
L
,
Δ
F
)
+
ℓ
F
F
(
Δ
F
,
Δ
F
)
]
-
[
μ
-
∥
Δ
K
∥
∞
]
∥
H
∥
1
,
where the directional derivatives of Φ
μ
(
L
,
F
) are given by
(56)
Φ
L
μ
Δ
L
+
Φ
F
μ
Δ
F
=
J
L
(
L
,
F
)
Δ
L
+
J
F
(
L
,
F
)
Δ
F
-
μ
∥
H
(
L
,
F
)
∥
1
=
Tr
⁡
(
Δ
L
Q
(
F
)
)
+
2
Tr
⁡
(
Δ
F
T
R
F
(
C
L
C
T
+
R
e
)
)
-
μ
∥
H
(
L
,
F
)
∥
1
.Proof.
See [17, Lemma 18.2].Following [17, Theorem 18.3] the search direction computed by the method DSQP is descent for the merit function Φ.Theorem 15.
Suppose that(
L
k
,
F
k
)
∈
𝒟
s is not a stationary point of the NSDP  (10)–(12) and the reduced Hessian Z
k
∇
2
ℓ
k
Z
k is positive definite. Then the search direction (
Δ
L
k
,
Δ
F
k
) by the method DSQP is a descent direction for the merit function Φ if μ is updated by  (54).


## 3.1. The Reduced Hessian Approach

The reduced Hessian approach relies on the existence of the null space and the right inverse of the Jacobian∇
H
T. In optimal control such operators are not available explicitly. However, it is possible to construct them from the problem itself; for more details see, for example, [2, 16, 25]. The two operators are defined as
(30)
Z
(
L
,
F
)
=
(
-
H
L
-
1
(
L
,
F
)
H
F
(
L
,
F
)
,
I
p
)
,
R
^
(
L
,
F
)
=
(
-
H
L
-
1
(
L
,
F
)
,
0
)
,
where I
p and 0 are the identity and the zero operators, respectively. We note that if the Jacobian ∇
H
T is surjective then such a right inverse R
^ exists. We also note that the range space of Z coincides with the null space of the Jacobian ∇
H
T; see [16, Lemma 2].Next, with the help ofZ and R
^ the following result shows that the linearized equality constraint is decomposed into two discrete Lyapunov equations.Lemma 9.
Let(
L
,
F
)
∈
𝒟
s be given and let (
Δ
L
,
Δ
F
)
∈
ℝ
n
×
n
+
p
×
r be the solution of the QP  (19). Then the linearized equality constraint of  (19) is decomposed into the following discrete Lyapunov equations:
(31)
Δ
L
n
=
A
(
F
)
Δ
L
n
(
A
F
)
T
+
H
(
L
,
F
)
,
(32)
Δ
L
t
=
A
(
F
)
Δ
L
t
A
(
F
)
T
+
B
Δ
F
L
F
T
+
L
F
Δ
F
T
B
T
,
where L
F is as given in Lemma 5.Proof.
(see also [16, Lemma 3]). The linearized equality constraint of (19) can be rewritten as
(33)
∇
L
H
(
L
,
F
)
Δ
L
t
+
∇
F
H
(
L
,
F
)
Δ
F
+
∇
L
H
(
L
,
F
)
Δ
L
n
+
H
(
L
,
F
)
=
0
,
where Δ
L
=
Δ
L
n
+
Δ
L
t. According to the definition of Z and R
^ the above equation can be decomposed as
(34)
∇
L
H
(
L
,
F
)
Δ
L
n
+
H
(
L
,
F
)
=
0
∀
Δ
L
n
,
H
L
(
L
,
F
)
Δ
L
t
(
Δ
F
)
+
H
F
(
L
,
F
)
Δ
F
=
0
∀
Δ
F
,
which yield the discrete Lyapunov equations  (31) and (32), respectively.Following the proof of Lemma9 the step (
Δ
L
,
Δ
F
) that solves the QP (19) with the help of Z and R
^ is decomposed as
(35)
(
Δ
L
,
Δ
F
)
=
Z
(
L
,
F
)
Δ
F
+
R
^
(
L
,
F
)
H
(
L
,
F
)
=
:
(
Δ
L
t
(
Δ
F
)
,
Δ
F
)
+
(
Δ
L
n
,
0
)
,
where (
Δ
L
t
(
Δ
F
)
,
Δ
F
) and (
Δ
L
n
,
0
) are the tangential and the normal steps, respectively.Consequently, one can obtain the normal step(
Δ
L
n
,
0
) by solving the discrete Lyapunov equation  (31) for Δ
L
n. Having obtained such a step component we solve the following tangential subproblem for Δ
F:
(36)
min
⁡
Δ
F
∈
ℝ
p
×
r
ψ
(
Δ
F
)
:
=
q
(
Δ
L
t
(
Δ
F
)
,
Δ
F
)
=
〈
Δ
L
t
(
Δ
F
)
,
∇
L
ℓ
(
·
)
〉
+
〈
Δ
F
,
∇
F
ℓ
(
·
)
〉
+
〈
Δ
F
,
∇
F
L
2
ℓ
(
·
)
Δ
L
n
〉
+
〈
Δ
L
t
(
Δ
F
)
,
∇
L
L
2
ℓ
(
·
)
Δ
L
n
〉
+
1
2
〈
Δ
F
,
Z
T
∇
2
ℓ
(
·
)
Z
Δ
F
〉
,
where
(37)
〈
Δ
F
,
Z
T
∇
2
ℓ
(
·
)
Z
Δ
F
〉
=
〈
Δ
L
t
(
Δ
F
)
,
∇
L
L
2
ℓ
(
·
)
Δ
L
t
(
Δ
F
)
〉
+
〈
Δ
F
,
∇
F
F
2
ℓ
(
·
)
Δ
F
〉
+
2
〈
Δ
F
,
∇
F
L
2
ℓ
(
·
)
Δ
L
t
(
Δ
F
)
〉
and Δ
L
t
(
Δ
F
) is the corresponding solution of the discrete Lyapunov equation  (32).In order to find the local solution of the minimization problem  (36) let us apply the necessary optimality conditions, which give the following result.Lemma 10.
Let(
L
,
F
,
K
)
∈
𝒟
s
×
ℝ
n
×
n and Δ
L
n
∈
ℝ
n
×
n be given. Moreover, let Δ
F
∈
ℝ
p
×
r be a local minimum for the problem (36). Then Δ
F solves the following linear matrix equation:
(38)
(
B
T
K
B
+
R
)
Δ
F
(
C
L
C
T
+
R
e
)
+
K
F
Δ
L
t
(
Δ
F
)
C
T
+
B
T
Δ
K
(
Δ
F
)
L
F
=
-
K
F
(
L
+
Δ
L
n
)
C
T
-
(
B
T
K
B
+
R
)
F
R
e
-
B
T
S
L
F
,
where Δ
L
t
(
Δ
F
), Δ
K
(
Δ
F
), and S, respectively, solve the discrete Lyapunov equations  (32) and
(39)
Δ
K
=
A
(
F
)
T
Δ
K
A
(
F
)
+
K
F
T
Δ
F
C
+
C
T
Δ
F
T
K
F
,
(40)
S
=
A
(
F
)
T
S
A
(
F
)
+
M
+
M
T
,
where K
F and L
F are as defined in Lemma 4 and
(41)
M
=
-
K
+
A
(
F
)
T
K
A
(
F
)
+
Q
(
F
)
.Proof.
We use the derivatives ofℓ given in Lemma 4 to construct ψ
(
Δ
F
). If we differentiate ψ with respect to Δ
F and apply the necessary optimality conditions, we obtain
(42)
Z
T
∇
2
ℓ
(
·
)
Z
Δ
F
=
-
∇
F
ℓ
(
·
)
-
∇
F
L
2
Δ
L
n
-
∂
∂
Δ
F
〈
Δ
L
t
(
Δ
F
)
,
∇
L
ℓ
(
·
)
〉
,
which is another representation of  (38). As an example of evaluating these derivatives we obtain the last term of  (42) as follows:
(43)
∂
∂
Δ
F
〈
Δ
L
t
(
Δ
F
)
,
∇
L
ℓ
(
·
)
〉
=
∂
∂
Δ
F
Tr
⁡
(
M
Δ
L
t
(
Δ
F
)
)
=
∂
∂
Δ
F
Tr
⁡
(
M
∑
k
=
0
∞
A
(
F
)
k
(
B
Δ
F
L
F
T
∑
k
=
0
∞
A
(
F
)
k
(
B
Δ
F
L
F
T
+
L
F
Δ
F
T
B
T
)
(
A
(
F
)
T
)
k
)
=
B
T
(
∑
k
=
0
∞
(
A
(
F
)
T
)
k
(
M
+
M
T
)
A
(
F
)
k
‍
)
L
F
=
B
T
S
L
F
,
where
(44)
Δ
L
t
(
Δ
F
)
=
∑
k
=
0
∞
‍
A
(
F
)
k
(
B
Δ
F
L
F
T
+
L
F
Δ
F
T
B
T
)
(
A
(
F
)
T
)
k
,
S
=
∑
k
=
0
∞
‍
(
A
(
F
)
T
)
k
(
M
+
M
T
)
A
(
F
)
k
are the exact solutions of the discrete Lyapunov equations  (32) and  (40), respectively.Because of using directional derivatives the resulting linear system  (38)–(40) is not given explicitly. However, one can modify the conjugate gradient (CG) method [17] and solve this system iteratively. The problem structure is exploited, where the partition of the unknown X
=
(
L
,
F
) is utilized.For a givenΔ
F calculated by the CG method both of Δ
L
t
(
Δ
F
) and Δ
K
(
Δ
F
) are uniquely determined by solving the discrete Lyapunov equations  (32) and  (39), respectively. The CG algorithm for solving the system  (38)–(40) is stated in the following lines, where we use T
i, D
i, and E
i at the ith CG inner iteration to denote the approximation of the step component Δ
F, the CG direction, and the residual, respectively.Algorithm 11 (modified CG method).
LetU
=
(
L
,
F
,
K
)
∈
𝒟
s
×
ℝ
n
×
n and Δ
L
n be given. Moreover, let A, B, C, Q, R, R
e, and V be given constant matrices. Set T
0
=
0, E
0
=
-
K
F
L
C
T
-
(
B
T
K
B
+
R
)
F
R
e, and D
0
=
E
0. Choose η and ϵ
∈
(
0,1
).
Fori
=
0,1
,
…, do(1)
Solve forΔ
L
t
(
D
i
) and Δ
K
(
D
i
) the discrete Lyapunov equations  (32) and  (39), respectively.
(2)
Calculate(45)
N
(
D
i
)
≔
(
B
T
K
B
+
R
)
D
i
(
C
L
C
T
+
R
e
)
+
K
F
Δ
L
t
(
D
i
)
C
T
+
B
T
Δ
K
(
D
i
)
L
F
.
(3)
IfTr
⁡
(
D
i
T
N
(
D
i
)
)
≤
ϵ
∥
D
i
∥
2, update the step Δ
F as follows: if i
=
0, set Δ
F
=
D
0; else, set Δ
F
=
T
i and exit.
(4)
Calculate the ratioθ
i
=
Tr
⁡
(
E
i
T
E
i
)
/
Tr
⁡
(
D
i
T
N
(
D
i
)
). Set T
i
+
1
=
T
i
+
θ
i
D
i and E
i
+
1
=
E
i
-
θ
i
N
(
D
i
).
(5)
If∥
E
i
+
1
∥
≤
η
∥
E
0
∥ set Δ
F
=
T
i
+
1 and exit.
(6)
Calculate the ratio(46)
β
i
=
Tr
⁡
(
E
i
+
1
T
E
i
+
1
)
Tr
⁡
(
E
i
T
E
i
)
.
SetD
i
+
1
=
E
i
+
1
+
β
i
D
i, i
:
=
i
+
1, and go to (
1
).
End (do)There are two exits in this CG method. The first one takes place when a negative or zero curvature is encountered by the calculated search direction. The second exit occurs when the norm of the new residual is sufficiently small.Having evaluated both ofΔ
L
k
n and (
Δ
L
t
(
Δ
F
k
)
,
Δ
F
k
), then according to the step decomposition  (35) the total step (
Δ
L
k
,
Δ
F
k
) is obtained. Consequently, the new iterate takes the form
(47)
(
L
k
+
1
,
F
k
+
1
)
=
(
L
k
+
α
k
Δ
L
k
,
F
k
+
α
k
Δ
F
k
)
,
where α
k
>
0 is the step size. Starting from α
0
=
1, we calculate the step size α
k by a backtracking rule (e.g., by halving α) until we satisfy the following sufficient decrease condition:
(48)
Φ
μ
k
(
L
k
+
α
k
Δ
L
,
F
k
+
α
k
Δ
F
)
≤
Φ
k
μ
k
+
τ
α
k
[
Φ
L
μ
k
(
L
k
,
F
k
)
Δ
L
+
Φ
F
μ
k
(
L
k
,
F
k
)
Δ
F
]
,
where τ
∈
(
0,1
/
2
) is a given parameter. Let α
~
k be the achieved step size. Starting from that α
~
k we satisfy the following stability condition by backtracking in the same way until an α
k
>
0 is reached such that
(49)
L
k
+
α
k
Δ
L
k
≻
0
,
Y
(
L
k
+
α
k
Δ
L
k
,
F
k
+
α
k
Δ
F
k
)
≻
0
.We assume that there always exists anα
k
>
0 such that  (48) is satisfied. If α is decreased and reached a lower bound α
min
⁡
>
0 without fulfilling  (49), then we might follow one of the ideas often used in modifying the Hessian matrix of Newton-based methods; see, for example, [17, Section 6.3]. Let W
k denote any of the two matrices L
k
+
α
min
⁡
Δ
L
k or Y
(
L
k
+
α
min
⁡
Δ
L
k
,
F
k
+
α
min
⁡
Δ
F
k
) with a spectral decomposition W
k
=
Q
~
Λ
Q
~
T. A correction matrix Δ
W
k of minimum Frobenius norm that ensures λ
min
⁡
(
W
k
+
Δ
W
k
)
≥
δ
>
0 is given by
(50)
Δ
W
k
=
Q
~
diag
⁡
(
t
i
)
Q
~
T
,
where
(51)
t
i
=
{
0
,
λ
i
>
δ
,
δ
-
λ
i
,
otherwise
.
The modified matrix is
(52)
W
k
+
Δ
W
k
=
Q
~
(
Λ
+
diag
⁡
(
t
i
)
)
Q
~
T
,
where λ
min
⁡
(
W
k
) means the smallest eigenvalue of W
k. Note that the nature of the stability constraint is an inactive constraint for most test problems.The Matlab function[
~
,
j
]
=
cholinc
(
M
) of the incomplete Cholesky factorization can be used for checking the positive definiteness, where the indicator j is used for that purpose. On the other hand, as mentioned at the end of Section 2 we can check the stability condition via the spectral radius of the closed-loop system matrix; namely, ρ
(
A
(
F
k
+
α
k
Δ
F
k
)
)
<
1.We use as a merit function the followingl
1-penalty function for globalizing the SQP method:
(53)
Φ
μ
(
L
,
F
)
=
J
(
L
,
F
)
+
μ
∥
H
(
L
,
F
)
∥
1
,
where μ
>
0 is the penalty parameter. This function is directionally differentiable.The following update rule for the penalty parameterμ is considered (see, e.g., [17, page 547]):
(54)
μ
k
=
{
μ
k
-
1
,
if
μ
k
-
1
≥
∥
Δ
K
k
∥
∞
+
γ
,
∥
Δ
K
k
∥
∞
+
2
γ
,
otherwise
,
where γ
>
0 is a given constant and Δ
K
k is obtained by solving the discrete Lyapunov equation  (39).The following algorithm demonstrates the SQP method that finds the local solution of the NSDP  (10)–(12).Algorithm 12 (the method DSQP).
Let(
L
0
,
F
0
,
K
0
)
∈
𝒟
s
×
ℝ
n
×
n be given. Moreover, let A, B, C, Q, R, R
e, and V be given constant matrices. Choose τ
∈
(
0,1
/
2
), γ
>
0, and ϵ
tol
∈
(
0,1
).
While∥
∇
U
ℓ
k
∥
>
ϵ
tol, do(1)
CalculateΔ
L
k
n solution of the discrete Lyapunov equation  (31).
(2)
For a givenΔ
L
k
n calculate the local solution (
Δ
L
t
(
Δ
F
k
)
,
Δ
F
k
) of the problem (36) by Algorithm 11. Then set (
Δ
L
k
,
Δ
F
k
)
=
(
Δ
L
k
n
,
0
)
+
(
Δ
L
t
(
Δ
F
k
)
,
Δ
F
k
).
(3)
CalculateΔ
K
k by solving (39) and then update the penalty parameter μ
k using  (54).
(4)
Starting fromα
0
=
1 apply a backtracking step-size rule that satisfies the sufficient decrease conditions  (48).
(5)
If(
L
k
+
α
k
Δ
L
k
,
F
k
+
α
k
Δ
F
k
)
)
∉
𝒟
s, then maintain the positive definiteness of  (49) as described above.
(6)
SetF
k
+
1
=
F
k
+
α
k
Δ
F
k, L
k
+
1
=
L
k
+
α
k
Δ
L
k, K
k
+
1
=
K
k
+
α
k
Δ
K
k, and k
:
=
k
+
1; then go to step (
1
). End (While)Step(
5
) of Algorithm 12 can also be modified as follows. Set F
k
+
1
=
F
k
+
α
k
Δ
F
k. Then obtain L
k
+
1
=
L
(
F
k
+
1
) and K
k
+
1
=
K
(
F
k
+
1
) by solving the discrete Lyapunov equations  (17) and  (18), respectively. This means we project the search direction from the space 𝒟
s to the control-variable space ℝ
p
×
r.The global convergence behavior of the method DSQP is beyond the scope of this work. However, we state some of the main global convergence results, which can be shown under the following general assumptions.Assumption 13.
In addition to Assumption6 we further assume that the functions J and H together with their derivatives are in Ω
⊆
𝒟
s. We further assume that the operator H
L
-
1 is uniformly bounded in 𝒟
s; furthermore, we assume that the sequences {
∇
2
ℓ
k
}, {
L
k
}, {
F
k
}, {
K
k
}, and {
Z
(
L
k
,
F
k
)
} are all bounded in Ω.The next lemma provides an upper bound for the directional derivative of the merit function by the computed search direction.Lemma 14.
Let(
Δ
L
,
Δ
F
,
Δ
K
) be generated by the SQP iteration  (21). Then the directional derivative of Φ
μ
(
L
,
F
) in the direction (
Δ
L
,
Δ
F
) satisfies
(55)
Φ
L
μ
Δ
L
+
Φ
F
μ
Δ
F
≤
-
[
ℓ
L
L
(
Δ
L
,
Δ
L
)
+
2
ℓ
L
F
(
Δ
L
,
Δ
F
)
+
ℓ
F
F
(
Δ
F
,
Δ
F
)
]
-
[
μ
-
∥
Δ
K
∥
∞
]
∥
H
∥
1
,
where the directional derivatives of Φ
μ
(
L
,
F
) are given by
(56)
Φ
L
μ
Δ
L
+
Φ
F
μ
Δ
F
=
J
L
(
L
,
F
)
Δ
L
+
J
F
(
L
,
F
)
Δ
F
-
μ
∥
H
(
L
,
F
)
∥
1
=
Tr
⁡
(
Δ
L
Q
(
F
)
)
+
2
Tr
⁡
(
Δ
F
T
R
F
(
C
L
C
T
+
R
e
)
)
-
μ
∥
H
(
L
,
F
)
∥
1
.Proof.
See [17, Lemma 18.2].Following [17, Theorem 18.3] the search direction computed by the method DSQP is descent for the merit function Φ.Theorem 15.
Suppose that(
L
k
,
F
k
)
∈
𝒟
s is not a stationary point of the NSDP  (10)–(12) and the reduced Hessian Z
k
∇
2
ℓ
k
Z
k is positive definite. Then the search direction (
Δ
L
k
,
Δ
F
k
) by the method DSQP is a descent direction for the merit function Φ if μ is updated by  (54).


## 4. Initial Feasible Solution

As can be seen in Algorithm12 an initial (
L
0
,
F
0
)
∈
𝒟
s is required to start the method DSQP. Such a matrix pair can be obtained by using any method that calculates suboptimal solutions (see, e.g., [21, Algorithm Fstab]) or by a pole assignment method; see, for example, [18]. Another alternative is to parameterize the system matrix A such that ρ
(
A
)
<
1 and therefore F
0 can be chosen as the zero matrix. Hence, according to Lemma 3 one can obtain L
0
≻
0 by solving the discrete Lyapunov equation  (17). This idea is described in [16] and is stated herein for completeness. In the parametrization approach one replaces the system matrix A by A
t
=
(
1
-
t
)
A, where 0
≤
t
<
1 is a parameter. Then instead of solving the NSDP problem  (10)–(12) we solve the following parameterized problem:
(57)
min
⁡
L
,
F
J
(
L
,
F
)
=
Tr
⁡
(
L
(
Q
+
C
T
F
T
R
F
C
)
)
+
Tr
⁡
(
F
T
R
F
R
e
)
,
s
.
t
.
H
t
(
L
,
F
)
=
-
L
+
A
t
(
F
)
L
A
t
(
F
)
T
+
V
+
B
F
R
e
F
T
B
T

=

0
,
Y
t
(
L
,
F
)

=

L

-

(
A
t

+

B
F
C
)
L
(
A
t

+

B
F
C
)
T

≻

0
,
L

≻

0
,
where A
t
(
F
)
=
A
t
+
B
F
C and the constant matrices A, B, C, Q, R
e, and V are as defined in Section 2. The matrix pair (
L
,
F
) must be chosen from the following parameterized set:
(58)
𝒟
s
t
=
{
(
L
,
F
)
∈
ℝ
n
×
n
×
ℝ
p
×
r
:
L
≻
0
,
Y
t
(
L
,
F
)
≻
0
}
.
Obviously as t
→
0 the set 𝒟
s
t approaches the set 𝒟
s.The initialization of the method DSQP can be as follows. Starting fromt
0
=
0, where t
0
∈
[
0,1
), we choose t
0 as the first candidate such that ρ
(
A
t
0
)
<
1 and set F
0
=
0. Having (
t
0
,
F
0
), then initial L
0 and K
0 can be obtained by solving the discrete Lyapunov equations:
(59)
L
=
A
t
0
(
F
0
)
L
A
t
0
(
F
0
)
T
+
V
,
K
=
A
t
0
(
F
0
)
T
K
A
t
0
(
F
0
)
+
Q
.
Then Algorithm 12 can be applied directly on the parameterized problem (57) but after replacing A
(
F
) by A
t
(
F
). The parameter t is updated in the main iteration loop according to the decreasing sequence {
t
k
}:
(60)
t
k
+
1
=
c
k
t
k
,
where {
c
k
}
↓
0. As mentioned above we choose t
0
∈
[
0,1
) as the first candidate such that ρ
(
A
t
0
)
<
1, starting from t
0
=
0.On the other hand, we might test the method DSQP for finding feasible points of the problem  (57), that is, finding points that satisfy the equality constraints as well as the positive definite constraints. This can be achieved if we solve the problem with zero objective function. For that purpose we set the objective function and all its derivatives to be zero in Algorithm 12. As a result and with the exception of  (15), the simplification of all formulas has no drawbacks on Algorithm 12. In that case the zero term of the objective function in  (15) implies  (18) to reduce to
(61)
K
=
A
t
(
F
)
T
K
A
t
(
F
)
,
which cannot be used in its current form to update the multiplier K as suggested in Algorithm 12, because it yields the trivial solution K
=
0. In order to overcome this difficulty the multiplier is updated instead by the modified equation:
(62)
K
=
A
t
(
F
)
T
K
A
t
(
F
)
+
I
n
.


## 5. Numerical Results

This section includes an implementation of the method DSQP using Matlab. Moreover, the method DSQP is compared versus Newton’s method (NM) [9] applied on the minimization problem  (5)-(6). Four test problems from the benchmark collection COMPlib [8] are considered in detail in addition to an application of a hydraulic turbine [26]. The benchmark COMPlib includes mainly continuous-time test problems. Therefore, the function c2d from the control system toolbox of Matlab is considered to convert the continuous-time model to the discrete-time counterpart.For those test problems withρ
(
A
)
<
1 according to Lemma 3 it is possible to start the method DSQP with F
0
=
0 and to obtain L
0 by solving the discrete Lyapunov equation  (17). However, for test problems with ρ
(
A
)
≥
1 an initial (
L
0
,
F
0
)
∈
𝒟
s is required; otherwise, the parametrization approach described in Section 4 can be applied. The function dlyap from the control system toolbox of Matlab is used for calculating the approximate solutions of the discrete Lyapunov equations. For all test problems the tolerance is ϵ
tol
=
1
×
1
0
-
6. The constant weight matrices Q and R and the covariance matrices V and R
e are taken as the unit matrix.Example 16.
The first test problem is a supersonic transport aircraft under Mach 2.7 flight condition; see [8, AC16]. The discrete-time data matrices are as follows:
(63)
A
=
[
0.9901
0.0012
0.0000
-
0.0994
-
0.0315
1.0000
0.0989
0.0014
-
0.6269
-
0.0004
0.9773
0.0376
0.1237
0.0001
0.0016
0.9892
]
,
B
=
[
0.0005
0.0004
0.0004
0.0040
0.0078
0.0794
-
0.0086
-
0.0066
]
,
C
=
I
4
.
Given the data matrices A, B, C, Q, R, R
e, and V the goal is to calculate a local solution (
L
*
,
F
*
)
∈
𝒟
s of the NSDP problem  (10)–(12) using the method DSQP. The corresponding F
* is considered as an approximate solution of the design problem  (2)–(4). The open-loop system is discrete-time Schur stable, because ρ
(
A
)
=
0.9989
<
1. Hence, according to Lemma 3 one can choose F
0
=
0 and obtain L
0
≻
0 the corresponding solution of the discrete Lyapunov equation  (17). The two methods NM and DSQP converge to the same local minimum after 25 and 21 iterations, respectively. Moreover, we have also tested the method DSQP for finding feasible points as explained at the end of Section 4. The achieved KKT and feasible points are, respectively
(64)
F
*
=
[
-
1.6109
0.1684
0.6795
6.3050
4.0166
-
0.8769
-
1.4994
-
2.9913
]
,
F
feas
=
[
-
107.5749
8.8006
10.7506
129.7845
18.3855
-
12.8539
-
13.9596
-
13.1650
]
,
where L
* and L
feas are the corresponding solutions of the discrete Lyapunov equation  (17). Table 1 shows the convergence behavior of the method DSQP to the local solution.
The method DSQP converges to a local minimum in less number of iterations than the method NM for most of the considered test problems. The next examples have also such a behavior.Table 1
Convergence behavior of the method DSQP for Example16.
k
J
k
∥
∇
U
ℓ
k
∥
ρ
(
A
(
F
k
)
)
i
cg
1
3.1135
e
+
005
8.4992
e
+
007
9.9838
e
-
001
0
2
2.1007
e
+
005
3.7773
e
+
007
9.9745
e
-
001
1
3
1.4256
e
+
005
1.6786
e
+
007
9.9581
e
-
001
1
⋮
⋮
⋮
⋮
⋮
19
1.5151
e
+
003
2.8400
e
-
003
9.6853
e
-
001
6
20
1.5151
e
+
003
7.6411
e
-
005
9.6853
e
-
001
6
21
1.5151
e
+
003
9.0844
e
-
009
9.6853
e
-
001
8Example 17.
This test problem is the discrete version of a transport aircraft model [8, AC8]. The problem has the dimensions n
=
9, p
=
1, and r
=
5. Therefore, we do not list the data matrices A, B, and C.
The open-loop system of this example is not discrete-time Schur stable, whereρ
(
A
)
=
1.0012
>
1. The two methods NM and DSQP are considered with the parametrization approach, where t
0
=
0.002 and c
k
=
(
0.8
)
k. Both methods successfully converge to the same local minimum after 21 and 18 iterations, respectively. The achieved KKT and feasible points are, respectively
(65)
F
*
=
[
2.1458
-
0.5334
-
1.4506
0.0386
1.0690
]
,
F
feas
=
[
2.5640
-
0.6294
-
1.6525
0.0469
1.2302
]
,
where L
* and L
feas are the corresponding solutions of the discrete Lyapunov equation  (17).
Table2 shows the behavior of the two methods while converging to the local solution. The objective functions and the convergence criteria are listed columnwise in that table. The two methods show fast local convergence rate starting from F
0
=
0.
Figure1 shows the state variables for the open- and closed-loop systems of the current example, where the achieved output feedback controller matrix F
* enforces all state variables to decay to the zero state.Table 2
Convergence behavior of the methods NM and DSQP for Example17.
k
J
^
k
∥
∇
J
^
k
∥
k
J
k
∥
∇
U
ℓ
k
∥
t
k
1
2.7105
e
+
003
1.3120
e
+
005
1
1.9178
e
+
003
5.4776
e
+
004
2.00
e
-
003
2
4.0079
e
+
003
2.2196
e
+
005
2
1.9541
e
+
003
3.9033
e
+
004
1.60
e
-
003
3
6.0553
e
+
003
3.4516
e
+
005
3
2.4908
e
+
003
3.4364
e
+
004
1.02
e
-
003
⋮
⋮
⋮
⋮
⋮
⋮
⋮
19
1.5413
e
+
003
2.3402
e
-
001
16
1.5413
e
+
003
1.9083
e
-
002
1.32
e
-
016
20
1.5413
e
+
003
4.6373
e
-
003
17
1.5413
e
+
003
2.2029
e
-
006
2.98
e
-
018
21
1.5413
e
+
003
7.4319
e
-
007
18
1.5413
e
+
003
3.8485
e
-
011
5.36
e
-
020Open- and closed-loop systems for the aircraft model of Example17.
(a)
(b)Example 18.
This test problem is the decentralized interconnected system [8, DIS3]. The problem has the dimensions n
=
8, p
=
4, and r
=
4. Therefore, we do not list the data matrices A, B, and C.
The open-loop system for this example is discrete-time Schur stable, whereρ
(
A
)
=
0.9620
<
1. Starting from F
0
=
0 where L
0 is the corresponding solution of  (17) the method DSQP successfully converges to a KKT point (
L
*
,
F
*
) in 14 iterations, while the method NM requires 16 iterations to converge to the same local minimum. Furthermore, the method DSQP requires 17 iterations to reach a feasible point. Both of the KKT and the feasible points are, respectively,
(66)
F
*
=
[
-
1.7344
-
0.5988
-
0.1937
0.0938
0.2451
-
0.1286
0.4894
-
0.5011
0.0037
-
0.0705
-
0.3206
0.0853
0.0606
-
0.2429
-
0.1231
-
0.1806
]
,
F
feas
=
[
-
10.5946
-
0.7911
-
0.4047
0.7671
5.5324
-
2.1539
3.5120
-
3.1956
0.2289
-
1.8917
-
12.0235
1.9402
0.4861
-
2.2963
-
0.3051
-
2.0438
]
,
where L
* and L
feas are the corresponding solutions of  (17). Furthermore, Table 3 shows the convergence behavior of the method DSQP to the local solution (
L
*
,
F
*
)
∈
𝒟
s.Table 3
Convergence behavior of the method DSQP for Example18.
k
J
k
∥
∇
U
ℓ
k
∥
ρ
(
A
(
F
k
)
)
i
cg
1
1.3612
e
+
003
1.7167
e
+
004
9.6199
e
-
001
1
2
9.6941
e
+
002
7.6729
e
+
003
9.4539
e
-
001
1
3
6.9504
e
+
002
3.4950
e
+
003
9.6356
e
-
001
1
⋮
⋮
⋮
⋮
⋮
12
6.7653
e
+
001
1.2410
e
-
002
9.0019
e
-
001
5
13
6.7653
e
+
001
7.0258
e
-
005
9.0021
e
-
001
8
14
6.7653
e
+
001
1.3558
e
-
008
9.0021
e
-
001
13Example 19.
This test problem represents a nuclear reactor model [8, REA3]. The dimensions of this example are n
=
12, p
=
1, and r
=
3. Therefore, we skip listing the data matrices A, B, and C.
The open-loop system for this example is not discrete-time Schur stable. By applying the parametrization approach, the methods NM and DSQP converge to the same local minimizer after 31 and 27 iterations, respectively. The achieved KKT and feasible points are, respectively,(67)
F
*
=
[
-
0.0350
-
0.3937
-
8.3675
]
,
F
feas
=
[
-
0.1269
-
1.4385
-
3.0157
]
,
where L
* and L
feas are the corresponding solutions of  (17).
Table4 shows the convergence behavior of the method DSQP to the KKT point.Table 4
Convergence behavior of the method DSQP for Example19.
k
J
k
∥
∇
U
ℓ
k
∥
ρ
(
A
t
k
(
F
k
)
)
t
k
i
cg
0
3.0141
e
+
008
1.5819
e
+
010
9.9779
e
-
001
2.00
e
-
003
1
1
2.7360
e
+
008
1.2031
e
+
010
9.9819
e
-
001
1.60
e
-
003
1
2
2.7696
e
+
008
1.0624
e
+
010
9.9877
e
-
001
1.02
e
-
003
1
⋮
⋮
⋮
⋮
⋮
⋮
25
1.8154
e
+
004
8.4501
e
-
001
9.9979
e
-
001
6.39
e
-
035
3
26
1.8154
e
+
004
5.7828
e
-
004
9.9979
e
-
001
1.93
e
-
037
3
27
1.8154
e
+
004
7.7346
e
-
007
9.9979
e
-
001
4.67
e
-
040
4Example 20 (hydraulic turbine monitoring).
This model was considered by Wang and Daley [26], where their objective was to study a fault detection method for that model. However, we aim to apply the SQP method for calculating the required SOF controller F
*. Such a controller is supposed to simultaneously stabilize the control system (2) and minimize the quadratic cost function (3) as well.
The hydraulic turbine generating set consists of a reservoir, an inlet pipe, a hydraulic turbine, a generator, a transmission line, and an infinite bus of the produced current. At the end of the inlet pipe a gate is located to control the water flow rate inside the pipe. A change in the gate position yields a change in the speed of the turbine and the produced power.
The dynamical behavior of the hydraulic turbine generating set is represented by the following equations (see [26]):
(68)
T
a
N
˙
(
t
)
=
M
(
N
(
t
)
,
H
(
t
)
,
G
(
t
)
)
-
M
d
(
t
)
,
H
(
t
)
=
-
T
w
Q
˙
(
t
)
,
Q
(
t
)
=
Q
(
N
(
t
)
,
H
(
t
)
,
G
(
t
)
)
,
where N
(
t
), H
(
t
), and G
(
t
) are the speed of the turbine, water pressure at the gate, and the gate opening position, respectively. Moreover, M
(
N
(
t
)
,
H
(
t
)
,
G
(
t
)
), Q
(
N
(
t
)
,
H
(
t
)
,
G
(
t
)
), and M
d
(
t
) are the driving torque on the gate, water flow inside the pipe, and the load into the infinite bus, respectively; T
a and T
w are total inertia and water starting time constant, respectively. Both of M and Q are assumed to be nonlinear and unknown functions in the variables N
(
t
), H
(
t
), and G
(
t
).
LetP
n
(
N
n
,
H
n
,
G
n
) be a nominated reference point. By linearizing the system  (68) about P
n one obtains the following continuous-time linear control system (see [26] for details):
(69)
x
˙
(
t
)
=
A
c
x
(
t
)
+
B
c
u
(
t
)
+
E
c
v
(
t
)
,
where
(70)
A
c
=
[
-
50.0
0
0
0.2500
-
0.1250
0.3570
1000.05
-
0.0250
-
2.2450
]
,
B
c
=
[
50
0
-
100
]
,
E
c
=
[
0
-
0.2500
-
0.0050
]
.
The corresponding discrete-time state space model takes the following form:
(71)
x
k
+
1
=
A
x
k
+
B
u
k
+
E
v
k
,
y
k
=
C
x
k
,
where
(72)
A
=
e
A
c
T
^
,
B
=
∫
0
T
^
‍
e
A
c
τ
B
c
d
τ
,
E
=
∫
0
T
^
‍
e
A
c
τ
E
c
d
τ
.
By choosing the sample time period T
^
=
0.1sec., the following data matrices are obtained:
(73)
A
=
[
0.0067
0
0
0.0590
0.9875
0.0331
1.6359
-
0.0022
0.7846
]
,
B
=
[
0.9933
-
0.0341
-
1.6315
]
,
(74)
E
=
[
0
-
0.250
-
0.005
]
.
The discrete-time system  (71) is Schur stable, where ρ
(
A
)
=
0.9871
<
1. Similar to the above examples we run the methods NM and DSQP on this model using the above data matrices. The two methods converge to the same stationary point after 12 and 11 iterations, respectively, where
(75)
F
*
=
[
-
0.0083
-
0.0117
]
.Figure2 shows the state variables of the open- and closed-loop systems. Although the open-loop system is Schur stable, the computed SOF controller F
* accelerates the decay of the state variables to the zero state.Open- and closed-loop systems for the hydraulic turbine model of Example20.
(a)
(b)In Table5 we further compare the SQP method for finding the local solution of the NSDP  (10)–(12) versus Newton’s method applied on the special formulation of the optimization problem  (5)-(6). The comparison is performed on 50 test problems from the benchmark collection [8]. In that table we list the problem name, the problem dimensions, the stability indicators for the open- and closed-loop control systems, and the number of iterations of the two methods to reach the stationary points. A dash “–” in the last two columns indicates that the method could not achieve the prescribed accuracy. The obtained results quite show that the SQP method outperforms Newton’s method with respect to number of iterations.Table 5
Comparison between the methods NM and DSQP for finding the local solutions of the problems   (5)-(6) and (10)–(12), respectively; test problems are from [8].
Problem
Problem dimension
Stability indicator
No. of iterations
n
p
r
ρ
(
A
)
ρ
(
A
(
F
*
)
)
NM
DSQP
AC1
5
3
3
1.0000
9.6958
e
-
001
11
7
AC3
5
2
4
0.9991
9.5419
e
-
001
20
19
AC4
4
1
2
1.2942
9.9501
e
-
001
11
10
AC6
7
2
4
0.9992
9.1586
e
-
001
19
17
AC7
9
1
2
1.0174
9.9693
e
-
001
7
5
AC8
9
1
5
1.0012
9.6072
e
-
001
21
18
AC15
4
2
3
0.9990
9.6497
e
-
001
23
22
AC16
4
2
4
0.9989
9.6853
e
-
001
25
21
AC17
4
1
2
0.9723
9.4295
e
-
001
13
11
HE1
4
2
1
1.0280
9.9116
e
-
001
5
4
HE2
4
2
2
0.9971
9.6999
e
-
001
16
13
HE3
8
4
6
1.0088
9.6678
e
-
001
9
6
REA1
4
2
2
1.2203
8.9332
e
-
001
8
4
REA2
4
2
3
1.2227
8.9821
e
-
001
7
4
REA3
12
1
3
1.0000
9.9979
e
-
001
31
27
DIS1
8
4
4
0.9912
9.5066
e
-
001
13
11
DIS2
3
2
2
1.1824
7.7117
e
-
001
8
4
DIS3
6
4
4
0.9620
9.0021
e
-
001
16
14
DIS4
6
4
6
1.1551
8.7595
e
-
001
9
6
AGS
12
2
2
0.9786
9.7976
e
-
001
7
5
TG1
10
2
2
0.9768
9.6791
e
-
001
19
17
UWV
8
2
2
0.9989
3.0749
e
-
001
20
19
IH
21
11
10
1.0000
9.7497
e
-
001
5
3
EB1
10
1
1
0.9990
9.9370
e
-
001
14
8
EB2
10
1
1
0.9990
9.9370
e
-
001
9
8
TF1
7
2
4
1
9.9241
e
-
001
—
8
PSM
7
2
3
0.9495
9.1393
e
-
001
11
9
NN2
2
1
1
1.0000
9.4185
e
-
001
3
2
NN4
4
2
3
0.9969
9.3285
e
-
001
14
12
NN8
3
2
2
0.9971
9.5459
e
-
001
14
13
NN11
16
3
5
0.9048
9.0783
e
-
001
8
6
NN13
6
2
2
1.2147
8.0133
e
-
001
7
4
NN15
3
2
2
1
9.9880
e
-
001
7
6
NN16
8
4
4
1.0000
9.8135
e
-
001
10
6
HF2D10
5
2
3
1.0133
9.2187
e
-
001
4
3
HF2D11
5
2
3
1.0253
8.0074
e
-
001
4
2
HF2D12
5
2
4
0.9830
9.4527
e
-
001
13
—
HF2D13
5
2
4
0.9756
8.0712
e
-
001
13
12
HF2D14
5
2
4
1.0227
9.6188
e
-
001
14
27
HF2D15
5
2
4
1.1689
8.7364
e
-
001
8
6
HF2D17
5
2
4
1.0557
8.0215
e
-
001
9
5
HF2D18
5
2
2
1.0285
9.9607
e
-
001
8
6
MFP
4
3
2
0.9979
9.9567
e
-
001
—
11
TMD
6
2
4
1
9.8667
e
-
001
8
5
LAH
48
1
1
0.9742
9.7416
e
-
001
3
2
DLR1
10
2
2
0.9995
9.9902
e
-
001
9
6
DLR2
40
2
2
0.9995
9.9943
e
-
001
6
5
DLR3
40
2
2
0.9995
9.9943
e
-
001
7
5
HF1
130
1
2
0.9981
9.9573
e
-
001
7
6
ISS1
270
3
3
0.9997
9.9962
e
-
001
11
7


## 6. Conclusion

In this paper, an SQP method with line search is introduced for finding the local solution of some NSDP problem resulting from the discrete-time static output feedback design problem. One of the two unknowns of the KKT point obtained by the numerical method is the output feedback controller matrixF that represents the approximate solution of the design problem.In order to initialize the SQP method a straightforward strategy was considered to create a starting feasible point with respect to the positive definite constraints  (12). The method DSQP outperformed Newton’s method [9] applied on a special formulation of the optimization problem. Furthermore, the obtained results quite show that the method DSQP converges from remote starting points to a KKT point at a fast local rate, which is typical for SQP methods. Feasible points of the NSDP problem have been also calculated by the SQP method. Such feasible points can be utilized, for example, by methods that look for the optimal solutions. Finally, the considered numerical experiments show that the proposed algorithm is effective in practical applications related to this problem class.



---

*Source: 683797-2014-01-02.xml*