Solution of Algebraic Riccati Equations Arising in Control of Partial Differential Equations
Carmeliza Navasca, Kirsten Morris · Lecture notes in pure and applied mathematics · 2005
Algebraic Riccati equations of large dimension arise when using approximations to design controllers for systems modelled by partial differential equations. For large model order direct solution methods based on eigenvector calculation fail. In this paper we describe an iterative method that takes advantage of several special features of these problems: (1) sparsity of the matrices (2) much fewer controls than approximation order and (3) convergence of the control with increasing model order. The algorithm is straightforward to code. Performance is illustrated with a number of standard examples. Introduction We consider the problem of calculating feedback controls for systems modelled by partial differential or delay differential equations. In these systems the state x(t) lies in an infinite-dimensional space. A classical controller design objective is to find a control u(t) so that the objective function ∫ ∞ 0 〈Cx(t), Cx(t)〉+ u∗(t)Ru(t)dt (1) is minimized where R is a positive definite matrix and the observation C ∈ L(X,Rp). The theoretical solution to this problem for many infinite-dimensional systems parallels the theory for finite-dimensional systems [9, 16, 17, e.g.]. In practice, the control is calculated through approximation. This leads to solving an algebraic Riccati equation A∗P + PA− PBR−1B∗P = −C∗C. (2) for a feedback operator K = −R−1B′P. (3) The matrices A, B, C arise in a finite dimensional approximation of the infinite dimensional system. Let n indicate the order of the approximation, m the number of control inputs and p the number of observations. Thus, A is n×n, B is n×m and C is p×n. There have been many papers written describing conditions under which approximations lead to approximating controls that converge to 2 the control for the original infinite-dimensional system [3, 10, 13, 16, 17, e.g.]. In this paper we will assume that an approximation has been chosen so that a solution to the Riccati equation (2) exists for sufficiently large n and also that the approximating feedback operators converge. For problems where the model order is small, n < 50, a direct method based on calculating the eigenvectors of the associated Hamiltonian works well [18]. Due the limitations of the calculation of eigenvectors for large non-symmetric matrices, this method is not suitable for problems where n becomes large. Unfortunately, many infinite-dimensional control problems lead to Riccati equations of large order. This is particularly evident in control of systems modelled by partial differential equations with more than one space dimension. For such problems, iterative methods are more appropriate. There are two methods that may be used: Chrandrasekhar and Newton-Kleinman iterations. In Chrandrasekhar iterations, the Riccati equation is not itself solved directly [2, 6]. A system of 2 differential equations K(t) = −B∗L∗(t)L(t), K(0) = 0, L(t) = L(t)(A−BK(t)), L(0) = C, is solved for K ∈ Rm×n, L ∈ Rp×n. The feedback operator is obtained as limt→−∞K(t). The advantage to this approach is that the number of controls m and number of observations p is typically much less than the approximation model order n. This leads to significant savings in storage. Furthermore, the matrices arising in approximation are typically sparse and this can be used in implementation of this algorithm. Unfortunately, the convergence of K(t) can be very slow and a very accurate algorithm suitable for stiff systems must be used. This can lead to very large computation times. Another approach to solving large Riccati equations is the Newton-Kleinman method [15]. The Riccati equation (2) can be rewritten as (A−BK)∗P + P (A−BK) = −C∗C −K∗RK. (4) We say a matrix Ao is Hurwitz if σ(Ao) ⊂ C− If A − BK is Hurwitz, then the above equation is a Lyapunov equation. An initial feedback K0 must be chosen so A − BK0 is Hurwitz. Define Si = A−BKi, and solve the Lyapunov equation S∗ i Xi + XiSi = −C∗C −K∗ i RKi (5) for Xi and then update the feedback as Ki+1 = −RBXi. If A − BK0 is Hurwitz, then Xi converges quadratically to P [15]. For an arbitrary large Riccati equation, this condition may be difficult to satisfy. However, this condition is not restrictive for Riccati equations arising in control of infinite-dimensional systems. First, many of these systems are stable even when uncontrolled and so the initial iterate K0 may be chosen as zero. Second, if the approximation procedure is valid then convergence of the feedback gains is obtained with increasing model order. Thus, a gain obtained from a lower order approximation, perhaps using a direct solution, may be used as an initial estimate, or ansatz, for a higher order approximation. This technique was used successfully in [12, 24]and later in this paper. In this paper we use a modified Newton-Kleinman iteration first proposed by Banks and Ito [2]as a refinement for a partial solution to the Chandraskehar equation. In that paper, they partially solve the Chandrasekhar equations and then use the resulting feedback K as a stabilizing initial guess for a modified Newton-Kleinman method. Instead of the standard Newton-Kleinman form (5) above, Banks and Ito rewrote the Riccati equation in the form (A−BKi)Xi + Xi(A−BKi) = −D∗ i Di (6) where Xi = Pi−1−Pi, Ki+1 = Ki−B Xi, and Di = Ki−Ki−1. The resulting Lyapunov equation is solved for Xi. Equation (6) has fewer inhomogeneous terms than the equation in the standard Solution of Algebraic Riccati Equations Arising in Control of Partial Differential Equations 3 Newton-Kleinman method (5). Also, the non-homogeneous term D depends on m inputs, not the observation C. In [2]a Smith’s method was used to solve the Lyapunov equations. Although convergent, this method is slow. Solution of the Lyapunov equation is a key step in implementing either modified or standard Newton-Kleinman. The Lyapunov equations arising in the Newton-Kleinman method have several special features: (1) the model order n is generally much larger than m or p and (2) the matrices are often sparse. We use a recently developed method [19, 23]that uses these features, leading to an efficient algorithm. In the next section we describe the implementation of this Lyapunov solver. We then use this Lyapunov solver with both standard and modified Newton-Kleinman to solve a number of standard control examples, including one with several space variables. Our results indicate that modified Newton-Kleinman achieves considerable savings in computation time over standard Newton-Kleinman. We also found that using the solution from a lower-order approximation as an ansatz for a higher-order approximation significantly reduced the computation time. 1. Solution of Lyapunov Equation Solution of a Lyapunov equation is a key step in each iteration of the Newton-Kleinman method. Thus, it is imperative to use a good Lyapunov algorithm. As for the Riccati equation, direct methods such as Bartels-Stewart [4]are only appropriate for low model order and do not take advantage of sparsity in the matrices. The Alternating Direction Implicit (ADI) and Smith methods are two well-known iterative schemes. These will be briefly described before describing a modification that leads to reduced memory requirements and faster computation. Consider the Lyapunov equation XAo + AoX = −DD∗ (7) where Ao ∈ Rn×n and D ∈ Rn×r. In the case of standard Newton-Kleinman, r = m + p while for modified Newton-Kleinman, r is only m. If Ao is Hurwitz, then the Lyapunov equation has a symmetric positive semidefinite solution X. For p < 0, define U = (Ao − pI)(Ao + pI)−1 and V = −2p(Ao + pI)DD(Ao + pI)−1. In Smith’s method [26], equation (7) is rewritten. The solution X is found by using successive substitutions: X = limi→∞Xi where Xi = UXi−1U + V (8) with X0 = 0. Convergence of the iterations can be improved by careful choice of the parameter p e.g. [25, pg. 197]. This method of successive substitition is unconditionally convergent, but has only linear convergence. The ADI method [20, 28]improves Smith’s method by using a different parameter pi at each step. Two alternating linear systems, (Ao + piI)Xi− 1 2 = −DD ∗ −Xi−1(Ao − piI) (9) (Ao + piI)X ∗ i = −DD∗ −X∗ i− 12 (Ao − piI) (10) are solved recursively starting with X0 = 0 ∈ Rn×n and parameters pi < 0. If all parameters pi = p then equations (9,10) reduce to Smith’s method. If the ADI parameters pi are chosen appropriately, then convergence is obtained in J iterations where J ' n. Choice of the ADI parameters is discussed below. If Ao is sparse, then the linear systems (9,10) can be solved efficiently. However, full calculation of the dense iterates Xi is required at each step. Setting X0 = 0, it can be easily shown that Xi is symmetric and positive semidefinite for all i, and so we can write X = ZZ∗ where Z is a Cholesky factor of X [19, 23]. (A Cholesky factor does not need to be square or be lower triangular.)