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).

1. Introduction

In this paper, the following nonlinear semidefinite programming (NSDP) problem is considered: where , , are sufficiently smooth matrix functions, where means that 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 is partitioned as , where represents the state variable and 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 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 point 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 Section 3 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 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 the notation () denotes that 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, denotes the identity matrix of order . We omit the argument when it is known from the context; for example, we use to denote .

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: where is the state vector, is the control (input) vector, is the measured output vector, and and are disturbance (noise) terms. Moreover , , and are given constant matrices of appropriate dimensions.

We will consider the following quadratic cost function as an objective function to be minimized: where is the expected value and and are given constant weight matrices.

The following control law is often used to close the above control system: where is the output feedback gain matrix representing the unknown.

In order to remove the dependency of the problem on the initial state vector the following assumption is often imposed.

Assumption 1. Assume that is a random variable uniformly distributed over the unit ball such that .
Under this assumption the optimal control problem (2)–(4) can be transformed into the matrix optimization problem (see [11]): where is the trace operator and solves the discrete Lyapunov equation:

Assumption 2. Throughout the paper we assume that , , and are given constant matrices. Furthermore, we assume that and are given constant weight matrices and and are given covariance matrices. We assume that while , and .

From Lyapunov stability theory the unknown must be chosen from the following set: 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 ]). Consider the SOF problem (2)–(4). An exists if and only if there exists , where where

According to Lemma 3 and by considering both of and as independent variables, where we skip the imposed dependency of on that was considered in (5)-(6), then the minimization problem can be stated as the following NSDP:

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 and 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 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 where is the Lagrange multiplier. The next lemma includes first- and second-order derivatives of and first derivatives of the function . Similar results can be found in [16, Lemma 2.1] but on a simpler problem. In the following discussion we denote and .

Lemma 4. Consider the NSDP problem (10)–(12). The Lagrangian function is twice continuously differentiable and its first- and second-order directional derivatives as well as first-order derivatives of the constraint function are the following: where and .

Proof. By using the directional derivatives with both of and 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 be a local minimizer of the problem (10)–(12) and let (10)–(12) be regular at . Then there exists a Lagrange multiplier satisfying with and the following KKT conditions: where is as defined in Lemma 4 and .

Clearly, the system (16) is nonlinear in the unknowns , , and . However, the last two equations of (16) might be viewed as discrete Lyapunov equations of the following form

The SQP method is based on minimizing successfully a quadratic program of the form: where 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: and updates the step by where is the step size.

Assumption 6. Throughout this paper we make the following assumptions.(H1)The mappings and are at least twice continuously differentiable.(H2)For a given the mapping is invertible.(H3)There exists an . Moreover, assume that for a given there always exists an such that .

The next theorem provides us with the linearization of the nonlinear system (16).

Theorem 7. For a given the linearized equations of the KKT system (16) are the following: where and are as defined in Lemma 4.

Proof. The linearization of the nonlinear system about the point in the direction of takes the form where is the Jacobian matrix, with which we directly obtain the system (23).

Under certain assumptions it is possible to show that the SQP direction 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 . However, if we write as a long column vector such that , then we are able to evaluate the Jacobian of using the basic properties of the Kronecker product which is an matrix.

Theorem 8. Let be given. The Jacobian matrix of is given by: where and are as defined in Lemma 4, , and is an commutation matrix defined as for any matrix . Furthermore, the linear system (23) is explicitly written as where

Proof. For any given three matrices , , and of appropriate dimensions the linear matrix equation is equivalent to the explicit linear system 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 . 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 where and are the identity and the zero operators, respectively. We note that if the Jacobian is surjective then such a right inverse exists. We also note that the range space of coincides with the null space of the Jacobian ; see [16, Lemma 2].

Next, with the help of and the following result shows that the linearized equality constraint is decomposed into two discrete Lyapunov equations.

Lemma 9. Let be given and let be the solution of the QP (19). Then the linearized equality constraint of (19) is decomposed into the following discrete Lyapunov equations: where is as given in Lemma 5.

Proof. (see also [16, Lemma 3]). The linearized equality constraint of (19) can be rewritten as where . According to the definition of and the above equation can be decomposed as which yield the discrete Lyapunov equations (31) and (32), respectively.

Following the proof of Lemma 9 the step that solves the QP (19) with the help of and is decomposed as where and are the tangential and the normal steps, respectively.

Consequently, one can obtain the normal step by solving the discrete Lyapunov equation (31) for . Having obtained such a step component we solve the following tangential subproblem for : where and 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 and be given. Moreover, let be a local minimum for the problem (36). Then solves the following linear matrix equation: where , , and , respectively, solve the discrete Lyapunov equations (32) and where and are as defined in Lemma 4 and

Proof. We use the derivatives of given in Lemma 4 to construct . If we differentiate with respect to and apply the necessary optimality conditions, we obtain which is another representation of (38). As an example of evaluating these derivatives we obtain the last term of (42) as follows: where 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 is utilized.

For a given calculated by the CG method both of and 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 , , and at the th CG inner iteration to denote the approximation of the step component , the CG direction, and the residual, respectively.

Algorithm 11 (modified CG method). Let and be given. Moreover, let , , , , , , and be given constant matrices. Set , , and . Choose and .
For , do(1)Solve for and the discrete Lyapunov equations (32) and (39), respectively.(2)Calculate (3)If , update the step as follows: if , set ; else, set and exit.(4)Calculate the ratio . Set and .(5)If set and exit.(6)Calculate the ratio
Set , , and go to .
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 and , then according to the step decomposition (35) the total step is obtained. Consequently, the new iterate takes the form where is the step size. Starting from , we calculate the step size by a backtracking rule (e.g., by halving ) until we satisfy the following sufficient decrease condition: where is a given parameter. Let be the achieved step size. Starting from that we satisfy the following stability condition by backtracking in the same way until an is reached such that

We assume that there always exists an such that (48) is satisfied. If is decreased and reached a lower bound 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 denote any of the two matrices or with a spectral decomposition . A correction matrix of minimum Frobenius norm that ensures is given by where The modified matrix is where means the smallest eigenvalue of . Note that the nature of the stability constraint is an inactive constraint for most test problems.

The Matlab function of the incomplete Cholesky factorization can be used for checking the positive definiteness, where the indicator 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, .

We use as a merit function the following -penalty function for globalizing the SQP method: where 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]): where is a given constant and 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 be given. Moreover, let , , , , , , and be given constant matrices. Choose , , and .
While , do(1)Calculate solution of the discrete Lyapunov equation (31).(2)For a given calculate the local solution of the problem (36) by Algorithm 11. Then set .(3)Calculate by solving (39) and then update the penalty parameter using (54).(4)Starting from apply a backtracking step-size rule that satisfies the sufficient decrease conditions (48).(5)If , then maintain the positive definiteness of (49) as described above.(6)Set , , , and ; then go to step . End (While)

Step of Algorithm 12 can also be modified as follows. Set . Then obtain and by solving the discrete Lyapunov equations (17) and (18), respectively. This means we project the search direction from the space to the control-variable space .

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 Assumption 6 we further assume that the functions and together with their derivatives are in . We further assume that the operator is uniformly bounded in ; furthermore, we assume that the sequences , , , , and 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 be generated by the SQP iteration (21). Then the directional derivative of in the direction satisfies where the directional derivatives of are given by

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 is not a stationary point of the NSDP (10)–(12) and the reduced Hessian is positive definite. Then the search direction 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 Algorithm 12 an initial 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 such that and therefore can be chosen as the zero matrix. Hence, according to Lemma 3 one can obtain 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 by , where is a parameter. Then instead of solving the NSDP problem (10)–(12) we solve the following parameterized problem: where and the constant matrices , , , , , and are as defined in Section 2. The matrix pair must be chosen from the following parameterized set: Obviously as the set approaches the set .

The initialization of the method DSQP can be as follows. Starting from , where , we choose as the first candidate such that and set . Having , then initial and can be obtained by solving the discrete Lyapunov equations: Then Algorithm 12 can be applied directly on the parameterized problem (57) but after replacing by . The parameter is updated in the main iteration loop according to the decreasing sequence : where . As mentioned above we choose as the first candidate such that , starting from .

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 which cannot be used in its current form to update the multiplier as suggested in Algorithm 12, because it yields the trivial solution . In order to overcome this difficulty the multiplier is updated instead by the modified equation:

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 according to Lemma 3 it is possible to start the method DSQP with and to obtain by solving the discrete Lyapunov equation (17). However, for test problems with an initial 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 . The constant weight matrices and and the covariance matrices and 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: Given the data matrices , , , , , , and the goal is to calculate a local solution of the NSDP problem (10)–(12) using the method DSQP. The corresponding is considered as an approximate solution of the design problem (2)–(4). The open-loop system is discrete-time Schur stable, because . Hence, according to Lemma 3 one can choose and obtain 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 where and 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.

Example 17. This test problem is the discrete version of a transport aircraft model [8, AC8]. The problem has the dimensions , , and . Therefore, we do not list the data matrices , , and .
The open-loop system of this example is not discrete-time Schur stable, where . The two methods NM and DSQP are considered with the parametrization approach, where and . Both methods successfully converge to the same local minimum after 21 and 18 iterations, respectively. The achieved KKT and feasible points are, respectively where and are the corresponding solutions of the discrete Lyapunov equation (17).
Table 2 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 .
Figure 1 shows the state variables for the open- and closed-loop systems of the current example, where the achieved output feedback controller matrix enforces all state variables to decay to the zero state.

Example 18. This test problem is the decentralized interconnected system [8, DIS3]. The problem has the dimensions , , and . Therefore, we do not list the data matrices , , and .
The open-loop system for this example is discrete-time Schur stable, where . Starting from where is the corresponding solution of (17) the method DSQP successfully converges to a KKT point 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, where and are the corresponding solutions of (17). Furthermore, Table 3 shows the convergence behavior of the method DSQP to the local solution .

Example 19. This test problem represents a nuclear reactor model [8, REA3]. The dimensions of this example are , , and . Therefore, we skip listing the data matrices , , and .
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, where and are the corresponding solutions of (17).
Table 4 shows the convergence behavior of the method DSQP to the KKT point.

Example 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 . 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]): where , , and are the speed of the turbine, water pressure at the gate, and the gate opening position, respectively. Moreover, , , and are the driving torque on the gate, water flow inside the pipe, and the load into the infinite bus, respectively; and are total inertia and water starting time constant, respectively. Both of and are assumed to be nonlinear and unknown functions in the variables , , and .
Let be a nominated reference point. By linearizing the system (68) about one obtains the following continuous-time linear control system (see [26] for details): where The corresponding discrete-time state space model takes the following form: where By choosing the sample time period sec., the following data matrices are obtained: The discrete-time system (71) is Schur stable, where . 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

Figure 2 shows the state variables of the open- and closed-loop systems. Although the open-loop system is Schur stable, the computed SOF controller accelerates the decay of the state variables to the zero state.

In Table 5 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.

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 matrix 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.

Conflict of Interests

The author declares that there is no conflict of interests regarding the publication of this paper.