Novel Results on Output-Feedback LQR Design

This article provides novel developments in output-feedback stabilization for linear time-invariant systems within the linear quadratic regulator (LQR) framework. First, we derive the necessary and sufficient conditions for output-feedback stabilizability in connection with the LQR framework. Then, we propose a novel iterative Newton's method for output-feedback LQR design and a computationally efficient modified approach that requires solving only a Lyapunov equation at each iteration step. We show that the proposed modified approach guarantees convergence from a stabilizing state feedback to a stabilizing output-feedback solution and succeeds in solving high-dimensional problems where other, state-of-the-art methods, fail. Finally, numerical examples illustrate the effectiveness of the proposed methods.


I. INTRODUCTION
O NE of the most fundamental problems in the control theory is the linear quadratic regulator (LQR) design problem [1]. The so-called infinite horizon linear quadratic problem of finding a control function u(t) = −Kx(t) for x 0 ∈ R n x that minimizes the cost functional with R > 0, Q − NR −1 N T ≥ 0 subject to x(0) = x 0 , anḋ has been studied by many authors [1], [2], [3], [4], where x(t) ∈ R n x , y(t) ∈ R n y , and u(t) ∈ R n u denote the state, measurable output, and the control input vectors, respectively. Furthermore, A ∈ R n x ×n x , B ∈ R n x ×n u , and C ∈ R n y ×n x are constant known matrices.
Often it is not possible or economically feasible to measure all the state variables. In this case, an output-feedback control law defined as would be more beneficial. However, finding an optimal outputfeedback control law in the form (3), which minimizes (1), is still one of the most important open questions in control engineering, despite the availability of many approaches and numerical algorithms, as it is pointed out in [5] and [6]. This is mainly due to the lack of testable necessary and sufficient conditions for output-feedback stabilizability. Furthermore, the majority of algorithms for the outputfeedback LQR design are formulated in terms of linear matrix inequalities (LMIs) [7], [8], [9], [10], [11], [12], [13], [14] or bilinear matrix inequalities (BMIs) [15], [16], [17], [18], [19], [20]. These algorithms are dependent on the used BMI/LMI solvers and could work well for small/medium-sized problems, but may fail to converge to a solution or become computationally too heavy as the problem size increases [21]. In addition, available iterative numerical algorithms with guaranteed convergence, such as [22] and [23], or algorithms using nonlinear programming (NLP) such as [24] and [25], as well as the recently introduced ray-shooting-method-based approaches, e.g., [21] and [26], unfortunately require a selection of an initial stabilizing outputfeedback gain. However, a direct procedure for finding such a gain is unknown and could be hard to get, as highlighted in [5]. The author in [27] has proposed a state-feedback projection theory to bypass the need of a stabilizing output-feedback gain. However, the introduced iterative controller design problem results in a coupled nonlinear matrix equations, and conditions for the existence and global uniqueness are neither introduced nor discussed. Furthermore, the proposed Newton approach ensures only sufficient conditions for output-feedback stabilizability. Finally, the authors in [28] proposed an algorithm that iterates a Riccati equation from an initial state-feedback solution, but it applies to a restrictive problem description and its convergence has not been proven.
In general, finding stabilizing static output-feedbacks (SOFs) is suspected to be nondeterministic polynomial-time hard (NPhard), as it is discussed in [5], [21], and [26]. The problem is known to be NP-hard if structural constraints or bounds are imposed on the entries of the controllers, see, e.g., [29] and This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see http://creativecommons.org/licenses/by/4.0/ [30]. Furthermore, minimal-norm SOFs with bounded entries, pole-placement, and simultaneous stabilization via SOFs are also considered to be NP-hard (see [21], [29], and [31], respectively). Moreover, the authors in [5] have reviewed results from the computational complexity theory to suggest that "such hope that someone can come up with an algorithm that can solve most of the SOFs problems in practice may not be realistic, at least for moderate and large-size problems." This prediction/prognosis from 20 years afar has been more or less proven since, as described previously, after years of extensive research in this filed, there are still unsolved problems, especially if we consider large-size problems.
Even though most of the output-feedback problems are considered to be NP-hard, we have shown in our recent paper [32] that within the LQR framework, it is possible to find SOFs in a reasonable time even for large-scale systems. In this article, we expand and complete our results from [32]. First, we derive the necessary and sufficient conditions for output-feedback stabilizability in connection with the LQR framework. Then, we propose a novel iterative output-feedback LQR design approach for linear time-invariant (LTI) systems, using Newton's method. Afterwards, we show that with a simple modification, a new iterative algorithm can be obtained which has a guaranteed convergence to an optimal output-feedback solution from any stabilizing state-feedback solution. In addition, the proposed modified algorithm requires solving only a Lyapunov equation at each iteration step, which is computationally much more tractable than algorithms in the literature, including approaches based on LMIs, BMIs, NLP, and ray-shooting methods. Finally, we propose/review some simple and useful modifications/extensions.
The mathematical notation of this article is as follows. The set of real and complex numbers are denoted by R and C, respectively. Given a matrix C ∈ R n y ×n x , its pseudoinverse is denoted by C + . For matrices A, B ∈ R n x ×n x , their Hadamard (Schur) and Kronecker products are denoted by A • B and A ⊗ B, respectively. Matrices, if not explicitly stated, are assumed to have compatible dimensions. The real part of a complex number z is denoted by (z). Finally, for any positive integer n x , the n x × n x identity and zero matrices are denoted by I n x , 0 n x ∈ R n x ×n x , respectively. For matrices, . means any matrix norm, consequently, . F , and . 2 means the Frobenius and induced 2-norm, respectively.

II. NECESSARY AND SUFFICIENT CONDITIONS FOR OUTPUT-FEEDBACK STABILIZABILITY
This section formulates the necessary and sufficient conditions for output-feedback stabilizability in the LQR framework, essential for the main results.
Considering the system (2) and the output-feedback control law (3), let us first recall some related terminology.
Definition 1: A square matrix A ∈ R n x ×n x is said to be stable if and only if for every eigenvalue λ i of A, (λ i ) ≤ 0.
Definition 2: The pair (A, B) is said to be stabilizable if and only if there exists a real matrix K ∈ R n u ×n x such that A − BK is stable.

Definition 3:
The pair (A, C) is said to be detectable if and only if there exists a real matrix L ∈ R n x ×n y such that A − LC is stable.
Definition 4: The system (2) is said to be SOF stabilizable if and only if there exists a real matrix F ∈ R n u ×n y such that A − BF C is stable.
Then, the novel necessary and sufficient stability conditions for output-feedback stabilizability in the LQR framework can be formulated as follows.
Theorem 1: The system (2) is SOF stabilizable if and only if the pair (A, B) is stabilizable, the pair (A, C) is detectable and there exist real matrices F ∈ R n u ×n y and G ∈ R n u ×n x such that where P ∈ R n x ×n x is the real symmetric positive semidefinite solution of for given Q ∈ R n x ×n x , N ∈ R n x ×n u , and R ∈ R n u ×n u matrices satisfying Proof: We will first prove the necessity of Theorem 1. Assume that A − BF C is stable for some F , i.e., the system (2) is output-feedback stabilizable. Then, the pair (A, B) is stabilizable since A − BK is stable for K = F C, and consequently, the pair (A, C) is detectable, since A − LC is stable for L = BF . Furthermore, because A − BF C is stable, there exists a unique symmetric positive semidefinite matrix P (see Appendix A for details), such that Rearranging (7), one can obtain Hence, setting G = F C − R −1 (B T P + N T ) implies the necessity of Theorem 1. Now assume that the pair (A, B) is stabilizable, the pair (A, C) is detectable and there exist real matrices F and G satisfying (4). From (4) and (5), it follows that (7) is satisfied. From the second condition, it follows that A − LC is stable for some L. Noting that it follows that the pair is detectable as well. Since P is symmetric and positive semidefinite, we conclude from (7) that A − BF C is stable, and hence, the sufficiency of Theorem 1 is proved as well. Remark 1: Similar conditions for output-feedback stabilizability have been obtained in [28,Th. 1], but for a restricted problem formulation with Q = C T C, R = I, and N = 0.
From Theorem 1 follows that if the system (2) is outputfeedback stabilizable, then there exists a state feedback gain K = F C such that A − BK is stable. If C is a square and nonsingular matrix, then we can easily express the output-feedback gain as F = KC −1 . However, for most of the output-feedback problems, the matrix C is nonsquare, i.e., noninvertible. Therefore, by expressing the output-feedback gain using a pseudoinverse as F = KC + , a so-called pseudoinverse error appears that can be calculated as which is identical to (4), since K = R −1 (B T P + N T ). Hence, from the aforementioned and the Theorem 1, it follows that if the system (2) is output-feedback stabilizable, then for given Q, R, and N matrices satisfying (6), there exists a real positive semidefinite matrix P such that (5) is fulfilled for The next identity is straightforward and used later to obtain the main results.
Identity 1: Suppose that F = R −1 (B T P + N T )C + and G = F C − R −1 (B T P + N T ). Then, the following statements are identical: Proof: The identity can be proved by substituting back all the denotations.

III. INFINITE HORIZON OUTPUT-FEEDBACK LQR DESIGN
Equations (14)- (16) are algebraic Riccati-like equations. In general, Newton's method and its modifications are widely used to solve algebraic Riccati equations [33], [34], [35]. Inspired by [33] and [34], in this section, we first propose a Newton'smethod-based algorithm to design stabilizing SOF controllers. Then, we show that by a simple modification, a computationally similarly tractable stabilizing output-feedback controller design approach can be obtained, while guaranteeing convergence from any initial state-feedback LQR solution. Finally, after a short sensitivity analysis, we show the relation of these approaches to the infinite horizon output-feedback LQR problem (i.e., to find a control law in the form (3), minimizing the cost function (1), subject to system dynamics (2) and initial state x 0 ).

A. Newton's Method for Stabilizing SOF Controller Design
The Fréchet derivative of a matrix function where the norm is any matrix norm and F = R or C [36], [37]. The Fréchet derivative, if it exists, can be shown to be unique [38]. Consider R defined by the Riccati like matrix equation (15). Then, its Fréchet derivative at the matrix P is given by (see Appendix C) where Z = C + C, and Now, we can formulate the Newton's method in Banach space (see [35] and [39]) for the solution of (15) as follows: Furthermore, we can compute P j+1 directly from (21) as where X j is solved from Equation (23) is a generalized Sylvester equation, which can be, based on Identity 1, rewritten to the form where Lemma 1: Suppose thatĀ ∈ F m×n ,B ∈ F p×q , and X ∈ F n×p . Then,  (24) gives The generalized Sylvester equation (24) has a unique solution if and only if K L is nonsingular. In this case, the solution can be obtained analytically as or can be approximated either by gradient-based iterative methods (such as [41], [42], and [43]), or by any other methods in the literature. The proposed Newton's method for the SOF controller design using (22) and (23) is summarized in Algorithm 1.
Remark 2: It follows from (24) that if C = I, then Z = I and the generalized Sylvester equation (24) reduces to Hence, the Algorithm 1 becomes equivalent to [34, Algorithm 1.1] for the state-feedback LQR design.
Remark 3: Algorithm 1 has a termination condition that depends on a constant > 0, → 0, which describes the expected tolerance on the numerical solution. For example = 10 −d means d digit desired accuracy in the numerical solution.
The results from this subsection are used only as an intermediate step to obtain the main results, the modified Newton's method. Therefore, global convergence and existence of a stabilizing solution remains to be proven. Although, standard local q-quadratic convergence results for Newton's method apply [ , Algorithm 1 requires an initial guess P 1 = P T 1 , which is close enough to a solvent for which the Fréchet derivative (i.e., K L ) is nonsingular. However, a direct procedure to get such initial guess is out of the topic of this article, since the results from this subsection are only used as an intermediate step to obtain the main results, the modified Newton's method.
In the next subsection, we show that with a simple modification, a new iterative algorithm can be obtained that has a guaranteed convergence from any stabilizing state-feedback solution to an output-feedback solution.

B. Modified Newton's Method for Stabilizing SOF Controller Design
In order to calculate the Newton step in the Algorithm 1, we need to solve the generalized Sylvester equation (24). In this subsection, we show that with a simple modification, we can approximate the Newton step and converge to a solution with similar computational effort, but with a guaranteed convergence from any state-feedback solution.
By freezing the matrix G in (15), the term G T RG becomes a constant during an iteration step and the Fréchet derivative reduces tô and the Newton's method to Equation (30) is a Lyapunov equation, which can be solved efficiently and with much less computational effort than solving (24) with (27) or with other iterative methods. The Algorithm 2 summarizes the proposed modified Newton's method for the SOF controller design using (29)-(31).

1) Convergence:
In this subsubsection, we show that under certain assumptions, Algorithm 2 has a guaranteed convergence from a stabilizing starting guess P 1 (i.e.,Ã −SP 1 is stable for someQ ≥ 0), to a stabilizing output-feedback solution.
Let us recall some results relating to the convergence proof.
Proof: For proof, see [44,Lemma 8,p. 7]. The next Proposition shows that if the conditions described in Theorem 1 hold, then with a stabilizing starting guess (P 1 ) the Algorithm 2 cannot fail due to a singular Lyapunov operator.
Proposition 1: Suppose that the conditions in Theorem 1 hold, and the pair (Ã,C q ) is detectable, whereQ =C T qCq is a full-rank factorization ofQ. If P 1 is stabilizing, and Algorithm 2 is applied to (15), then the Lyapunov operator of the Lyapunov equation in step 7 from Algorithm 2 is nonsingular for all j and the sequence of approximate solutions X j is well defined.
Proof: Suppose that the pair (Ã,C q ) is detectable. From step 7 from Algorithm 2 applied to (15), we can get sinceQ andS are positive semi-definite, due to Q − NR −1 N T ≥ 0 and R > 0. From (32) follows that ifÃ −SP j is stable, thenÃ −S(P j + X j ) is also stable. Furthermore, Lemma 3 implies that P j + X j is positive semidefinite. The Lyapunov operator corresponding to the Lyapunov equation in step 7 from Algorithm 2 is well defined, precisely as for X j ∈ R n x ×n x and j = 1, 2, . . .. Let us recall the following Lemma.  (15), then the Newton step (of Algorithm 2) is a descent direction of R(P j + X j ) F . It follows that we have R(P j + X j ) F ≤ R(P j ) F and R(P j + X j ) F = R(P j ) F if and only if R(P j ) = 0. That is, the residual decreases as long as P j is not a solution of (15).
Remark 6: It is important to note that P 1 ≥ P 2 , where P 1 is the initial guess, is not true in general. This is one of the drawbacks of Newton's methods. In [33] and [34], the authors have introduced a step-size control (for state-feedback LQR design), which can efficiently solve the problem of a potentially disastrous first Newton step.
Collecting the results so far, we have the following convergence result for the modified Newton's method.
Theorem 2: Suppose that the pair (Ã, B) is stabilizable, the pair (Ã,C q ) is detectable, and there exist real matrices F and G such that F C − R −1 (B T P + N T ) = G. If Algorithm 2 is applied to (15) with a stabilizing starting guess P 1 (i.e.,Ã − BK 1 is stable for someQ ≥ 0), then P * = lim j→∞ P j exists and is the stabilizing solution of the generalized Riccati-like equation (15).

Remark 7:
The assumption that the pair (Ã,C q ) is detectable, whereQ =C T qCq is a full-rank factorization ofQ, is a requirement even for the standard state-feedback LQR design.  Ã, B) and observability of the pair (Ã,C q )).
Remark 9: From Theorem 2 follows that the convergence rate of Algorithm 2 is at least sublinear. We have observed from the examples studied later in Section V that the convergence rate is in fact linear, ifÃ −SP * has no eigenvalues on the imaginary axis, although further investigation is needed for a formal proof. IfÃ −SP * has eigenvalues on the imaginary axis, the convergence behavior remains an open problem (as it is still an open problem even for the standard state-feedback LQR design; see, for example, [34, Remark 1.1]).
Remark 10: If the system (2) is stabilizable and (Ã,C q ) is detectable, then the standard state-feedback LQR solution for (2) for someQ ≥ 0 always gives a P 1 for whichÃ −SP 1 is stable.

C. Sensitivity Analysis
It is well known that the Newton's-method-based approaches, in general, are highly sensitive to ill-conditioning. Condition numbers measure the sensitivity of a problem to perturbation in the data. The unstructured absolute condition number cond(R(P )) can be expressed in terms of the Fréchet derivative of R(P ) in (15), evaluated at P By applying a Frobenius norm cond(R(P )) = max the problem of computing cond(R(P )) reduces to finding the 2-norm of K L . The relative condition number of R(P ) at P , denoted by rcond(R(P )), can be written in terms of the absolute condition number cond(R(P )) (see, [46, Sec. 2, p. 776]) as rcond(R(P )) = cond(R(P )) P R(P ) .
For structured condition numbers of R(P ) at P , as well as for level-2 condition numbers using higher order Fréchet derivatives, see [46] and [47] and references therein. The effect of condition numbers on convergence will be investigated later in Section V.

D. Connection to Infinite-Horizon LQR With Output Feedback
This subsection describes the relation of Algorithms 1 and 2 to infinite-horizon LQR with output feedback. First, let us recall the necessary conditions for the solution of the LQR problem with output feedback, i.e., the existence of a control law in the form (3) minimizing (1) subject to (2) with R > 0, Q − NR −1 N T ≥ 0 and x(0) = x 0 .
Lemma 5: The necessary conditions for the solution of the LQR problem with output feedback are given by with The dependence of X x 0 in (38) in the initial states x 0 makes the optimal gain dependent on the initial state through (38). In many applications, x 0 may not be known (which is typical for the output-feedback design, as it is pointed out in [48]). It is usual (see, for example, [49]), to sidestep this problem by replacing where E{X x 0 } = E{x 0 x T 0 } is the initial autocorrelation of the state. Usually, it is assumed that nothing is known of x 0 except that it is uniformly distributed on a surface described by X x 0 . Most of the papers on the output-feedback LQR design assume that the initial states are uniformly distributed on the unit sphere, i.e., X x 0 = I (e.g., [14], [18], [48], and [50]).
The next theorem describes how Algorithms 1 and 2 are connected to the initial condition problem described previously.
By assuming that Y = I, (41) reduces to which is identical to the step 3 in Algorithms 1 and 2. Furthermore, from (38) for Y = I, it follows that Finally, setting G = F C − R −1 (B T P + N T ) and by rearranging (37), we can get (14), which is equivalent (see Identity 1) to the step 5 in Algorithms 1 and 2. Hence, the proof is completed. Remark 11: The initial state x 0 is generally free and so is Y , which is a function of X x 0 . Hence, instead of guessing X x 0 = I, we may guess for Y = I, and thus, the nonlinearity in Y in (37)-(39) disappears. So, one can get a simple Riccati-like equation (37), which can be solved easily using Algorithm 1 or 2.
A direct comparison of setting X x 0 = I versus Y = I will be investigated in Section V.

E. Output-Feedback LQR Problem With Known Initial Conditions
In the previous Section III-D, we have shown the relation of Algorithms 1 and 2 to output-feedback LQR problem, and that the proposed algorithms involve less nonlinearities compared to other approaches in the literature when the initial conditions are not given priory, i.e., x 0 is unknown. In this subsection, we show (see Algorithm 3) how the Algorithm 2 can be extended if the initial conditions are known.
The numerical examples in Section V suggest that Algorithm 3 converges to a solution if R(P ) is well-conditioned. But at this writing, we are not aware of a proof for this conjecture.
and if all uncontrollable state variables of the system (2) are asymptotically stable, then A c is negative definite. Hence, Y i+1 exists and is symmetric and positive definite. It follows that Y i+1 has a full rank and if CY i+1 C T is nonsingular, the Lyapunov operator corresponding to the Lyapunov equation in step 7 from Algorithm 3 is well defined, precisely as for Y i+1 ∈ R n x ×n x and i = 1, 2 . . .. Hence, Algorithm 3 cannot fail due to a singular Lyapunov operator in step 7. Therefore, if in step 3 of Algorithm 3, the Algorithm 2 succeeds in finding F i at each step, then Algorithm 3 produces a sequence of symmetric where Y * is the solution satisfying (37)- (39).

IV. USEFUL TECHNIQUES, EXTENSIONS, AND MODIFICATIONS
The control law (3) is defined in an SOF form. Many different controller structures can be transformed to this SOF form, like proportional-integral (PI), realizable proportional-integralderivative (PID f ), realizable proportional-derivative (PD f ), realizable derivative (D f ), even full/reduced order dynamic output-feedback controllers (DOF), dynamic output-feedback with integral and realizable derivative part (DOFID f ), or dynamic output-feedback with realizable derivative part (DOFD f ), by augmenting the system with additional state variables. For more information, see [14].
Since the proposed algorithms (Algorithms 1-3) belong to the LQR framework, all the well-known techniques, modifications, and extensions of the standard LQR design can be applied here as well. Therefore, one can apply the following: 1) Bryson's rules [51,Sec. 5.2] for selecting the weighting matrices Q and R; 2) methods/techniques in [48] for damping, decoupling, tracking, disturbance rejection, etc., controller design; 3) techniques in [52] for different eigenvalue placements (pole-placement techniques in LQ) and guaranteed convergence rate; 4) techniques in [53] and [54] for frequency weighting (frequency shaped LQ); 5) some other methods/techniques in the LQR framework; see, e.g., [48] and [51] and references therein.

V. NUMERICAL EXAMPLES
In order to show the viability of the previous proposed algorithms (Algorithms 1-3 [55]. As algorithms to be compared, the iterative LMI (iLMI) method from [14] and the BMI formulation of the outputfeedback LQR (OFLQR) problem (see Appendix B, Lemma [57] and by Mosek LMI solver v. 9.0.103 [58] using YALMIP R20190425 [59]. Finally, the proposed algorithms (Algorithms 1-3) have been implemented in MATLAB programming language (see , where for the Algorithms 2 and 3, for the step 7, the built-in MATLAB lyap subrutin has been used, and for the Algorithm 1, for the step 7, (27) is used. Furthermore, in order to simplify the code of the Algorithm 2, the Identity 1 has been used and (15) has been replaced with (16).
MATLAB implementations and examples are fully provided in Listing 1-3. The first set of examples, can be downloaded from repository, 1 while the second set of examples are all the SOF stabilizable plants from the COMPl e ib library, which is freely available (see [55]).

A. First Set of Examples
The first set of examples contains 1000 SOF stabilizable examples in ten groups generated by MATLAB's rss subrutin. 2 Each group represents 100 examples with different size of system order (see Table I). Hence, we can test the behavior and effectiveness of the proposed algorithms and compare them to other output-feedback algorithms in the LQR framework for increasing number of states.
The weighting matrices have been chosen as Q = C T C, R = I, and N = 0. The initial Lyapunov matrix for Algorithms 1-3 is the optimal Lyapunov matrix from the standard state-feedback LQR design. The stopping criterion and maximal iteration number for Algorithms 1-3 have been chosen as = 10 −12 and maxIteration = 9 × 10 6 . Finally, for Algorithm 3 and for the iLMI and BMI methods, the initial state matrix has been chosen as X x 0 = I, i.e., it has been assumed that the initial states are uniformly distributed on the unit sphere.
The effect of increasing the number of states on the running time of one iteration and on the number of iterations of Algorithms 1 and 2 are shown in Figs. 1 and 2. The effect of increasing the number of states on the rcond(R(P )) and on the number of iterations of Algorithms 1 and 2 is shown in Fig. 3. It can be observed that the number of iterations, and therefore, the overall running times of Algorithms 1 and 2 are sensitive to the relative condition number of (14) (rcond(R(P ))). That was to be expected, since it is well known that the Newton's-method-based approaches, in general, are sensitive to ill-conditioning. The effect of increasing the number of states on the average running times and on solved examples of Algorithms 1-3, and of the iLMI and BMI methods are shown in Fig. 4. It can be observed  that the Algorithm 2 outperformed all the other algorithms and methods since it has solved all the examples in this test set, while the running time was very close or sometimes better than the running time of Algorithm 1. Furthermore, it can be observed that even though for the Algorithm 1 we do not have a convergence proof from a state-feedback solution (as for Algorithm 2), it has solved many more examples than the iLMI or the BMI methods. Furthermore, Algorithm 3 for n x ≥ 30 failed to converge to a solution for few examples. It should be noted that for all those examples the rcond(R(P )) > 10 10 , and even the Algorithm 2 has struggled, since the number of iterations for those examples was higher than 10 4 , while for the rest of the examples in those groups was smaller with almost 1 or 2, sometimes with 3-4 orders of magnitude. This is the reason for that large trajectory of number of iterations of Algorithm 2 in Fig. 2 for n x > 10.
Finally, Fig. 5 compares how far the actual linear quadratic cost is for some randomly generated initial state conditions (x 0 within a unit sphere) for X x 0 = I and for Y = I from the optimal output-feedback cost (minimizing the linear quadratic cost for the given x 0 ). It can be observed that the distribution of distances from the optimal cost for different initial conditions for Y = I, i.e., for Algorithm 2, is comparable with the choice of setting X x 0 = I. Hence, Algorithm 2 is a viable approach for the outputfeedback LQR design with unknown initial conditions.

B. Second Set of Examples
The second set includes all the continuous-time SOF stabilizable examples from the COMPl e ib library [55] (see Table II), i.e., all the examples expect the ten reduced-order control (ROC) instances, and the four examples pointed out in [60], which are not continuous-time stabilizable (REA4, NN3, NN10, and NN12, indicated with the a superscript in Table II). This rich-full library contains benchmark examples from a wide spectrum of real-world applications and academic problems even with n x > 4000. Therefore, we can test the behavior and effectiveness of our proposed algorithms on large-scale stable/unstable The weighting matrices for all examples have been chosen as Q = C T C + αI, R = I, and N = 0, where α = 10 −9 . The αI has been introduced to ensure the positive-definiteness of the Q matrix, as most of the examples in the COMPl e ib library are ill-conditioned causing the C T C to became negative definite due to numerical errors. The initial Lyapunov matrix for the Algorithm 2 is the optimal Lyapunov matrix from the standard  state-feedback LQR design. Furthermore, the stopping criterion and maximal iteration number for Algorithm 2 have been chosen as = 10 −12 and maxIteration = 9 × 10 6 . For the iLMI and BMI methods, the X x 0 has been chosen as X x 0 = I, and all other solver related parameters for the Mosek LMI and Penlab BMI solvers have been kept as default.
The results summarized in Table II indicate that the proposed approach is superior compared to BMI and iLMI formulations. While the proposed Algorithm 2 has solved 92% (101/110) of the examples, the iLMI formulation 63% (59/110) and the BMI formulation only 17% (19/110). That is, Algorithm 2 solved 71% more examples than iLMI, and 432% more than the BMI method. In addition, even with the built-in MATLAB lyap subrutin, which is not well-suited for large-scale problems, we were able to solve examples with order higher than 4000 within minutes. The LAH example, see Table II, well demonstrates that the proposed approach is computationally much more tractable than approaches based on LMIs and/or BMIs . While the Algorithm 2 converged to a solution in 1. 23  . The COMPl e ib library well demonstrates that without proper regularization or preconditioning, the proposed algorithms may fail to converge due to numerical issues. The same is true for the iLMI and BMI methods. Table II also indicates that with system balancing (in our case with the built-in MATLAB balreal subrutin), we are able to solve 26% more examples with the Algorithm 2 than without any system balancing. This number can be further increased by preconditioning the Lyapunov/Sylvester equation within the Newton's method similarly as in [35]. Furthermore, the proposed approach can be easily extended with exact line-search, similarly as it is done in [33] and [34] to speed up the convergence and to reduce the overall running time even further.
In Remark 9, we have discussed that the convergence rate of the Algorithm 2 is at least sublinear. However, we have observed from the aforementioned examples that the convergence rate is in fact linear, ifÃ −SP * has no eigenvalues on the imaginary axis. Convergence rates of the proposed Algorithms 1 and 2 on different COMPl e ib plants (AC3, AC4, and DIS2), for initial Lyapunov matrices obtained from the standard state-feedback LQR design, are shown in Figs. 6 and 8. Convergence rates of the Algorithms 1 and 2 on the COMPl e ib plant AC3 for random initial Lyapunov matrices are shown in Figs. 7 and 9. From the figures, it can be observed that the convergence rate is quadratic for Algorithm 1 and is linear for Algorithm 2, if the initial Lyapunov matrix is calculated by the standard LQR design, and that it becomes quadratic/linear in the neighborhood of the solution when the Lyapunov matrix is randomly initialized.
In summary, the proposed algorithms, unlike other methods based on linear/bilinear matrix inequalities, NLPs, and rayshooting methods, can solve almost all the SOF examples in the COMPl e ib library (most of them within milliseconds). Until now it has been achieved only by some multivariate direct search methods applied to SOF stabilization in [60]. However, the author's attention in that publication was restricted to SOF Fig. 6. Convergence rate of the Algorithm 2 on COMPl e ib plants AC3, AC4, and DIS2. The initial Lyapunov matrix is obtained by the standard state-feedback LQR design. It can be observed that the convergence rate is linear. Fig. 7. Convergence rate of the Algorithm 2 on COMPl e ib plant AC3 for random initial Lyapunov matrices. It can be observed that the convergence rate becomes linear in the neighborhood of the solution. Fig. 8. Convergence rate of the Algorithm 1 on COMPl e ib plants AC3, AC4, and DIS2. The initial Lyapunov matrix is obtained by the standard state-feedback LQR design. It can be observed that the convergence rate is quadratic. Fig. 9. Convergence rate of the Algorithm 1 on COMPl e ib plant AC3 for random initial Lyapunov matrices. stabilization only, i.e., no attempt was made to optimize closedloop performance criteria relevant to control engineering. For more comparison, the readers are referred to [61], where the authors have evaluated different controller design approaches on the COMPl e ib library (including frequency-domain approaches minimizing H ∞ and/or H 2 norms as well).

VI. CONCLUSION
This article provides novel results on the SOF controller design for LTI systems in the LQR framework. Even though most of the output-feedback control problems are considered to be NP-hard, we show that within the LQR framework, it is possible to find SOFs in sublinear (linear) time even for large-scale systems. The proposed framework, with novel necessary and sufficient conditions for output-feedback stabilizability, opens the possibility to use well-known methods, such as the Newton's methods, to design SOFs with guaranteed convergence from a stabilizing state-feedback solution to a stabilizing outputfeedback solution. Hence, we can get computationally efficient approaches that succeed in solving high-dimensional problems where other, state-of-the-art methods fail.
The usability, tractability, and effectiveness is also verified on more than 1000 numerical examples in addition to all the SOF stabilizable plants from the COMPl e ib library. The proposed Algorithm 2, unlike other methods based on linear/bilinear matrix inequalities, NLPs, can solve almost all the SOF examples in the COMPl e ib library (most of them within milliseconds). Along this line, numerical results also indicate that the proposed algorithms suffer from the well-known drawbacks of Newton's methods. Therefore, regularization and proper scaling is needed to improve the usability of the proposed approaches for illconditioned problems.
In terms of future works, the Lyapunov equation can be preconditioned within the Newton's method similarly as in [35], and the proposed approaches can be easily extended with exact line search, similarly as it is done in [33] and [34].
APPENDIX A EXISTENCE OF P ≥ 0 Lemma 6: Let F ∈ R n u ×n y be given such that A − BF C is stable. Substitution of u(t) = −F y(t) = −F Cx(t) into the cost function (1) gives Since Q + C T F T RF C − C T F T N T − NF C ≥ 0, because R > 0 and Q − NR −1 N T ≥ 0, and since A − BF C is stable, it follows that the Lyapunov equation: (A − BF C) T P + P (A − BF C) has a unique solution P ≥ 0.

APPENDIX B BMI FORMULATION OF THE OFLQR DESIGN PROBLEM
Lemma 7: The SOF LQR design problem can be formulated as the following optimization problem: subject to BMI and LMI constraints Proof: Assume that the Lyapunov candidate is positive semidefinite. Then, from the Bellman-Lyapunov inequalitẏ V (x(t)) + J(x(t)) ≤ 0 →V (x(t)) ≤ −J(x(t)) where which indicates that the closed-loop system is stable. Integrating both sides from 0 to ∞, we can obtain the upper bound of the cost function which completes the proof.

APPENDIX C FRÉCHET DERIVATIVE OF (15)
By substituting back F = R −1 (B T P + N T )C + , G = F C − R −1 (B T P + N T ),S Q = NR −1 N T , andS A = BR −1 N T to (15), and perturbing with X, we can get R(P + X) = +Q +Ã T (P + X) + (P + X)Ã − (P + X)S(P + X) + Z T (P + X)S(P + X)Z + Z T (P + X)S A Z − Z T (P + X)S(P + X) we get R(P + X) = R(P ) + L(P, X) + E o (X), or R(P + X) − R(P ) − L(P, X) = E o (X). If P is the solution of (15), then R(P + X) F = R(P ) F = 0, and consequently, X F = 0. From this, follows that lim X F →0 E o (X) F = 0, and L(P, X) is the Fréchet derivative of (15) at P .