Closed-Form Solutions to Continuous-Time Algebraic Riccati Equation for Second-Order Systems

Continuous-time algebraic Riccati equations (CAREs) often appear in control and optimization problems such as linear quadratic regulator (LQR) [1], linear quadratic Gaussian [2], Kalman filter [2], $H2$ and $H∞$ methods [1], the gap (⁠$δ$⁠) and $ν$-gap (⁠$δν$⁠) [1] computation, etc. As a result, they are employed in applications such as control design for helicopters [3], unmanned aerial vehicles [4], spacecraft formation [5], space structures [6], induction motors [7], and others.

These CAREs are generally solved using the Schur eigenvector method [8], which is formulated for first-order systems. The first-order state-space form, with the state dimension $n¯$⁠, is used to create a Hamiltonian matrix (⁠$2n¯×2n¯$⁠), whose rearranged eigenvectors are used to compute the CARE solution. Since the computational effort is proportional to $n¯3$ [9], the method in Ref. [8] works well only for moderate size systems (⁠$n¯≤100$⁠) [8]. It becomes inefficient and prohibitive for larger systems [9–11], particularly when the system is poorly conditioned. Thus, it becomes a limiting factor for applications that require solving CARE for large-size systems such as the control of large space structures (⁠$n¯=128$⁠; [12]), finite element techniques [13], etc. Solving CARE in real-time is challenging even for small-size systems. For example, Zhao et al. [14] mention that solving CARE in real-time for only six states is considerably difficult for a processor onboard aerial robot.

Furthermore, the linear dynamics of many systems can be expressed as second-order linear time-invariant (LTI) ordinary differential equations (ODEs), described by the mass (⁠$M$⁠), stiffness (⁠$K$⁠), damping (⁠$C$⁠), and input $B$ matrices. The standard approach to solve CARE for such a system, of size $n×n$⁠, is to convert it to an equivalent first-order system, of size $n¯×n¯=2n×2n$⁠, which doubles the system size. Thus, the problem in solving CARE is exacerbated for second-order systems, such as in Refs. [6,12], reducing the limit of $n$ to $n≤50$⁠, beyond which the method in Ref. [8] becomes progressively less reliable and inefficient. Moreover, conversion of second-order systems to their first-order form leads to the loss of desirable properties such as symmetry and sparsity of the second-order matrices [15].

However, second-order systems have a special structure that can be used to solve CAREs in a closed-form manner, i.e., expressing the solution explicitly as an analytical formula. This special structure has been exploited in Refs. [15,16] to find simple, specialized conditions for the controllability and observability of second-order systems. In this article, we use a similar strategy to obtain simple closed-form solutions to CAREs for second-order systems. Prior works [17–20] that attempted to find analytical or closed-form solutions to CARE for second-order systems are briefly described next, revealing their limitations.

Prior works are evaluated based on whether they actually solve the CARE and if certain requirements are met. Specifically, the solution matrix $P$ and design matrix $Q$ must be positive semidefinite, and another design matrix $R$ must be positive definite. Also, prior works that claim a unique positive semidefinite solution should ensure the satisfaction of the detectability and stabilizability conditions.

In Ref. [17], a method that could simplify solving CARE for systems where $M$ and $K$ are simultaneously congruently diagonalizable is suggested. To that end, the second-order system is transformed into a first-order one, and a modal transformation is applied to diagonalize $M$ and $K$ (see Ref. [21] for an extensive analysis of simultaneous diagonalizability of matrices). The obtained CARE is broken down into four equations involving the transformed matrices. One of these equations is treated as a CARE by restrictively assuming the off-diagonal blocks of $P$ to be symmetric and solved using the eigenvector-based approach [8]. However, this assumption, while mathematically convenient, is not generally true and can lead to an incorrect solution. Also, the approach in Ref. [17] can lead to a non-symmetric, and hence, non-positive semidefinite $P$⁠. Therefore, the method in Ref. [17] does not always solve the CARE.

In Ref. [18], a proposed closed-form solution for the CARE for undamped fully actuated second-order systems is expressed in terms of scalar weights $α$ and $β$⁠. The assumption of no damping significantly limits the application of Ref. [18].

Furthermore, Hanks and Skelton [19] propose a method to solve CARE for second-order systems but only if the sensors and actuators are collocated. The method includes the design of the matrix $Q$ (or $Q^)$⁠, which is intrinsic to the CARE formulation, and the CARE solution $P$ in terms of scalar constants $a$ and $b$⁠, and another matrix $θ^$⁠. However, there is no discussion about how $a$ and $b$ should be selected to ensure that $Q$ (or $Q^$⁠) is positive semidefinite, as required by the CARE. Moreover, to prove that $P$ is positive definite, certain constraints are imposed which, if solved, could lead to $a$ and $b$ such that $P≻0$⁠. These constraints might be unsolvable, or be very expensive to solve, and levy structural impositions on the system limiting its application to very particular cases.

In Ref. [20], a closed-form solution to the CARE for second-order systems is presented in terms of scalars $α$ and $β$⁠. This solution satisfies the CARE and the authors discuss how to find $α$ and $β$⁠. However, the work in Ref. [20] also assumes collocated actuators and sensors. Under this assumption, the authors prove that the obtained $P$ always exponentially stabilizes the system using LQR. By additionally proving the observability condition, they invoke the argument of a unique stabilizing solution of Ref. [8] and conclude that their solution is unique and positive semidefinite. However, the proof is valid only for systems with collocated actuators and sensors. The authors do not address the more realistic and practical non-collocated actuators and sensors systems [22].

Finally, the icare method [23] in Matlab uses a modified version of the method of Ref. [8] for first-order systems, by employing a proprietary scaling routine. This extends the limit of $n¯≤100$ (or $n≤50$⁠) without much loss in accuracy, albeit with a very significant increase in computation time. In comparison, the unscaled icare of Matlab produces poor results as $n¯$ increases beyond $100$ (or $n>50$⁠).

Main Contribution: We present two novel closed-form solutions, in Theorems 2 and 3, to CARE for second-order systems. The first applies when the mass and stiffness matrices {$M,K$} can be simultaneously congruently diagonalized. The second is valid for any positive definite $M,C$⁠, and $K$⁠. The solutions are unique, positive semidefinite, and optimal for the cost in (3). They are expressed in terms of the second-order system matrices $M,C,K$⁠, and offer physical insights. The proposed methods work very well for $n>>50$ and achieve better accuracy than the Matlab’s unscaled icare [23] and similar accuracy to Matlab’s scaled icare, but are significantly faster than both of them.

Section 2 discusses the problem statement. The novel methods to solve CARE are presented in Secs. 3 and 4, followed by results in Sec. 5. Some applications are presented in Sec. 6 and the conclusions follow in Sec. 7.

Moreover, if the conditions of stabilizability and detectability are, respectively, strengthened to controllability and observability, then the solution becomes positive definite. Note that the conditions of stabilizability and detectability are critical. Thus, any method that claims to solve the CARE for a unique positive semidefinite, stabilizing solution that is optimal to the cost in (3) should prove these conditions.

Moreover, $P11$⁠, from (15), depends on terms that may be non-symmetric. For example, $C^$ is not necessarily diagonal, the term $C^P12T$ might be non-symmetric because even if we assume $P12$ to be symmetric, the product of symmetric matrices is not generally symmetric. Similarly, $P22L^P12T$ is not necessarily symmetric. This would result in $P$ not being necessarily symmetric, in contradiction with a key requirement for the correct CARE solution. Furthermore, it should be shown that the solution (⁠$P11,P12,P22$⁠) and the design matrices (⁠$Q11,Q12,Q22$⁠) from (13) to (16) would result in positive semidefinite $P$ and $Q$⁠. The work in Ref. [17] does not address these issues in finding a positive semidefinite solution to CARE. Also, the proposed solution is not closed-form and still depends on the eigenvector method [8].

We now present novel closed-form methods that solve CARE for second-order systems.

Furthermore, since $L$ and $Q11$ are positive definite and semidefinite matrices, respectively, and $X$ is full-rank, the matrix $J=XTL−1X+Q11$ is positive definite. Thus, all of the conditions for the existence of a solution to (17), as per Ref. [24], are satisfied.

Thus, if we can show that the proposed $P$ and $Q$ are positive semidefinite along with the stabilizability and detectability conditions, then this $P$ would be the same as that obtained from Ref. [8] by virtue of uniqueness. Therefore, the proposed $P$ would be the unique positive semidefinite and stabilizing solution to (4) that is optimal for the LQR cost in (3) if the choice of ${Q11,Q12,Q22}$ ensures that matrices $Q$ and $P$ are positive semidefinite, and the conditions of stabilizability and detectability are met.

We now modify the presented method by transforming the system to its modal form. Furthermore, the assumption made on $Z$ being invertible, after (25), is no longer needed in the modified method.

Since $C^$ is a real symmetric matrix, $C^1/2$ can be directly obtained by orthogonal diagonalization. Note that all the matrix products inside the square bracket in (48) are symmetric. Thus, their positively weighted sum is also symmetric with real eigenvalues. Hence, the inequality in (49) is valid. Thus, the matrix $Q$ is positive semidefinite for $a$ given by (49). This choice of $a$ also ensures that $Q22$⁠, from (43), is positive semidefinite.

Now, we present another novel closed-form solution to (4). Here, we do not require to transform the system to its modal form, and only require the matrices $M,C,K$ to be positive definite. The solution is in terms of $M,C,K$⁠, as shown in Theorem 3.

The formulas for $Q,P$ were designed by choosing candidate elements of $Q$ (⁠$Q11,Q12,Q22$⁠) and $P$ (⁠$P11,P12,P22$⁠) such that Eqs. (17)–(19) were satisfied. Additionally, these elements should also satisfy $P⪰0$ and $Q⪰0$⁠. It can be easily verified that the given choice of $P,Q,R$ exactly satisfy the associated CARE. In the following, we prove that the matrices $P$ and $Q$ are positive semidefinite.

Note that the terms I and II are both symmetric and so is their positively weighted sum. Hence, the combined matrix on the right-hand side of (67) has only real eigenvalues. Thus, we have shown that $Q⪰0$ for any $γ>0$ and for $a$ as per (67).

For $P$ and $Q$ to be simultaneously positive semidefinite, the scalar $a$ can be chosen as greater than or equal to the maximum of the values given by (67) and (71). If the inequality in (67) is considered strict, then $Q≻0$⁠. Then, based on the arguments given toward the end of method-I (after (54)), $P$ is the unique positive definite and stabilizing solution to (4) that is optimal for the LQR cost in (3).

Now, we present simulation results that numerically verify the proposed method. Since method-II works for any positive definite $M,C,K$⁠, it is more general than method-I, which requires modal decomposition. Thus, we only present results from method-II. We compare method-II’s results with the scaled and unscaled versions of matlab’s icare on the basis of the residual and computation time.

Figure 1 shows the $2$-norm of the residual left behind, using scattered points, by different methods when solving the CARE. The median of the residuals is also shown for an interval of $50$ by solid lines. It is clear that method-II and the scaled and unscaled icare work almost equally well when $n≤50$ for second-order systems. However, as $n$ increases, Matlab’s unscaled icare starts giving progressively large residuals. The poor performance of the unscaled icare could be attributed to the errors generated while computing the eigenvectors for large-size systems (⁠$n>50$⁠). In contrast, method-II avoids eigenvector calculations and its associated errors.

The scaled icare produces much better results than the unscaled one because of the proprietary scaling routine. However, it does so at the cost of significant time and computation penalty, as shown in Fig. 2. The proposed method-II achieves similar accuracy to the scaled icare, and at a cost even smaller than the unscaled icare, see Fig. 2.

Figure 2 shows the time taken to compute the CARE solution with respect to the size of the second-order system, clearly illustrating the advantage of method-II from the computation cost perspective. Since Matlab’s unscaled icare requires the computation and reordering of the Hamiltonian matrix’s eigenvectors, it becomes computationally expensive for large-size systems. In addition, the scaled icare scales the input matrices according to a proprietary scaling routine, which further increases the computation time and cost. In contrast, method-II is closed-form and avoids eigenvector computation. Consequently, the time taken to compute the solution is significantly less for method-II, as evident from Fig. 2. Thus, method-II would perform much better than other methods in solving CARE for large-size systems in real-time.

In this section, we discuss the advantages of our work in the context of some specific engineering problems such as vibration control.

Many large space structures require solving CARE for second-order systems to achieve vibration control [6,12]. Consider Ref. [6], where the authors improve the active vibration control of large space structures through structural modifications. To that end, the authors of Ref. [6] first solve a CARE to find the LQR solution. Then, with this solution, the authors obtain the closed-loop system, compute its eigenvalues, and the damping factors corresponding to these eigenvalues. Since these damping factors (and the eigenvalues) are dependent on the system’s matrices $M,C,andK$⁠, any structural modification will change these matrices as well as the damping factors. Thus, required changes in the damping factors can be obtained through modification of these matrices by changing properties such as the area, length, etc. of the constituent structural elements.