Filetype PDF | Posted on 27 Jan 2023 | 2 years ago
Authentication
digital resources
136x Filetype PDF File size 0.39 MB Source: ww3.math.ucla.edu
Iterative Solution of Algebraic Riccati Equations for Damped Systems
Kirsten Morris and Carmeliza Navasca
Abstract—Algebraic Riccati equations (ARE) of large dimen- common iterative methods are Chandrasekhar [3], [9] and
sion arise when using approximations to design controllers for Newton-Kleinman iterations [14], [17], [28]. Here we use a
systems modelled by partial differential equations. We use a modification to the Newton-Kleinman method first proposed
modified Newton method to solve the ARE. Since the modified by by Banks and Ito [3] as a refinement for a partial solution
Newton method leads to a right-hand side of rank equal to to the Chandraskehar equation.
the number of inputs, regardless of the weights, the resulting
Lyapunov equation can be more efficiently solved. A low-rank Solution of the Lyapunov equation is a key step in im-
Cholesky-ADI algorithm is used to solve the Lyapunov equation plementing either modified or standard Newton-Kleinman.
resulting at each step. The algorithm is straightforward to code. The Lyapunov equations arising in the Newton-Kleinman
Performance is illustrated with an example of a beam, with method have several special features: (1) the model order
different levels of damping. Results indicate that for weakly n is generally much larger than the number of controls m
damped problems a low rank solution to the ARE may not
exist. Further analysis supports this point. or number of observations p and (2) the matrices are often
INTRODUCTION sparse. We use a recently developed method [20], [26] that
uses these features, leading to an efficient algorithm. Use
Aclassical controller design objective is to find a control of this Lyapunov solver with a Newton-Kleinman method is
u(t) so that the objective function described in [27], [6], [7], [8].
! ∞!x(t),Qx(t)"+u∗(t)Ru(t)dt (1) In the next section we describe the algorithm to solve the
Lyapunov equations. We use this Lyapunov solver with both
0 standard and modified Newton-Kleinman to solve a number
is minimized where R is a symmetric positive definite matrix of standard control examples, including one with several
and Q is symmetric positive semi-definite. It is well-known space variables. Our results indicate that modified Newton-
that the solution to this problem is found by solving an Kleinman achieves considerable savings in computation time
algebraic Riccati equation over standard Newton-Kleinman. We also found that using
A∗P +PA−PBR−1B∗P =−Q. (2) the solution from a lower-order approximation as an ansatz
for a higher-order approximation significantly reduced the
for a feedback operator computation time.
K=−R−1B′P. (3) I. DESCRIPTION OF ALGORITHM
If the system is modelled by partial differential or delay One approach to solving large Riccati equations is the
differential equations then the state x(t) lies in an infinite- Newton-Kleinman method [17]. The Riccati equation (2) can
dimensional space. The theoretical solution to this problem be rewritten as
for many infinite-dimensional systems parallels the theory for ∗ ∗
finite-dimensional systems [11], [12], [18], [19]. In practice, (A−BK) P +P(A−BK)=−Q−K RK. (4)
the control is calculated through approximation. The matrices We say a matrix A is Hurwitz if σ(A ) ⊂ C .
A, B, and C arise in a finite dimensional approximation of o o −
Theorem 1.1: [17] Consider a stabilizable pair (A,B)
the infinite dimensional system. Let n indicate the order of with a feedback K0 so that A − BK0 is Hurwitz. Define
the approximation, m the number of control inputs and p S =A−BK,andsolve the Lyapunov equation
i i
the number of observations. Thus, A ∈ Rn×n, B ∈ Rn×m S∗P +P S =−Q−K∗RK (5)
and C ∈ Rp×n. There have been many papers written i i i i i i
describing conditions under which approximations lead to for P and then update the feedback as K =−R−1B∗P.
approximating controls that converge to the control for the i i+1 i
Then
original infinite-dimensional system [4], [12], [15], [18], lim P = P
[19]. In this paper we will assume that an approximation i→∞ i
has been chosen so that a solution to the Riccati equation (2) with quadratic convergence.
exists for sufficiently large n and also that the approximating The key assumption in the above theorem is that an initial
feedback operators converge. Unfortunately, particularly in estimate, or ansatz, K , such that A − BK is Hurwitz
0 0
partial differential equation models with more than one space is available. For an arbitrary large Riccati equation, this
dimension, many infinite-dimensional control problems lead condition may be difficult to satisfy. However, this condition
to Riccati equations of large order. A survey of currently is not restrictive for Riccati equations arising in control of
used methods to solve large ARE can be found in [8]. Two infinite-dimensional systems. First, many of these systems
are stable even when uncontrolled and so the initial iterate TABLE I
K may be chosen as zero. Second, if the approximation COMPLEXITY OF CF-ADI AND ADI [20]
0
procedure is valid then convergence of the feedback gains is
obtained with increasing model order. Thus, a gain obtained CF-ADI ADI
2
Sparse A O(Jrn) O(Jrn )
from a lower order approximation, perhaps using a direct 2 3
method, may be used as an ansatz for a higher order Full A O(Jrn ) O(Jrn )
approximation. This technique was used successfully in [14],
[24], [28]. TABLE II
Amodification to this method was proposed by Banks and CHOLESKY-ADI METHOD
Ito [3]. In that paper, they partially solve the Chandrasekhar Given A and D
equations and then use the resulting feedback K as a o
Choose ADI parameters {p ,...,p } with ℜ(p ) < 0
√ 1 J i
∗ −1
stabilizing initial guess. Instead of the standard Newton- Define z = −2p (A +p I) D
1 1 o 1
and Z =[z ]
Kleinman iteration (5) above, Banks and Ito rewrote the 1 1
For i = 2,...,J
Riccati equation in the form √−2pi+1 ∗ −1
Define W =( √ )[I − (p −p)(A +p I) ]
i −2pi i+1 i o i+1
∗ ∗ (1) z = W z
(A−BK) X +X(A−BK) = −D RD, i i i−1
i i i i i i (2) If $z$ > tol
Di = Ki −Ki−1, Ki+1 = Ki−B∗Xi. Z =[Z z ]
i i−1 i
Solution of a Lyapunov equation is a key step in each Else, stop.
iteration of the Newton-Kleinman method, both standard and
modified. Consider the Lyapunov equation
XAo+A∗X=−DD∗ (6) symmetric A the optimal parameters can be easily calculated.
o A heuristic procedure to calculate the parameters in the
where Ao ∈ Rn×n is Hurwitz and D ∈ Rn×r. Factor general case is in [26]. If the spectrum of A is complex,
Q into C∗C and let nQ indicate the rank of C. In the we estimated the complex ADI parameters as in [10]. As
case of standard Newton-Kleinman, r = m + nQ while for the spectrum of A flattens to the real axis the ADI parameters
modified Newton-Kleinman, r = m. If Ao is Hurwitz, then are closer to optimal.
the Lyapunov equation has a symmetric positive semidefinite If A is sparse, then sparse arithmetic can be used in
solution X. As for the Riccati equation, direct methods such o
as Bartels-Stewart [5] are only appropriate for low model calculation of Xi. However, full calculation of the dense
order and do not take advantage of sparsity in the matrices. iterates Xi is required at each step. By setting X0 = 0,
In [3] Smith’s method was used to solve the Lyapunov it can be easily shown that Xi is symmetric and positive
−1 semidefinite for all i, and so we can write X = ZZ∗
equations. For p < 0, define U = (A − pI)(A + pI)
o o where Z is a Cholesky factor of X [20], [26]. (A Cholesky
∗ −1 ∗ −1
and V = −2p(Ao+pI) DD (Ao+pI) . The Lyapunov factor does not need to be square or be lower triangular.)
equation (6) can be rewritten as Substituting Z Z∗ for X , setting X = 0, and defining
i i i 0
∗ ∗ −1
X=U XU+V. M =(A +pI) weobtain the following iterates
i o i
In Smith’s method [30], the solution X is found by using Z = "−2p M D
successive substitutions: X = limi→∞Xi where 1 " 1 1
Z = [ −2pMD, M(A∗−pI)Z ].
Xi = U∗Xi−1U +V (7) i i i i o i i−1
with X0 = 0. Convergence of the iterations can be improved Note that Z ∈ Rn×r, Z ∈ Rn×2r, and Z ∈ Rn×ir.
by careful choice of the parameter p e.g. [29, pg. 297]. 1 2 i
This method of successive substitition is unconditionally The algorithm is stopped when the Cholesky iterations
convergent, but has only linear convergence. converge within some tolerance. In [20] these iterates are
The ADI method [21], [32] improves Smith’s method by reformulated in a more efficient form, using the observation
using a different parameter p at each step. Two alternating that the order in which the ADI parameters are used is
i irrelevant. This leads to the algorithm shown in Table II.
linear systems, This solution method results in considerable savings in
(A∗ +piI)Xi−1 = −DD∗−Xi−1(Ao−piI) computation time and memory. In calculation of the feedback
o 2
∗ ∗ ∗ ∗ K, the full solution X never needs to be constructed. Also,
(A +piI)X = −DD −X 1(Ao−piI)
o i i−2 the complexity of this method (CF-ADI) is also considerably
are solved recursively starting with X0 = 0 ∈ Rn×n and less than that of standard ADI as shown in Table I. Recall
parameters pi < 0. If all parameters pi = p then the above that for standard Newton-Kleinman, r is the sum of the rank
equations reduce to Smith’s method. If the ADI parameters of the state weight nQ and controls m. For modified Newton-
pi are chosen appropriately, then convergence is obtained Kleinman r is equal to m. The reduction in complexity
in J iterations where J ≪ n. The parameter selection of CF-ADI method over ADI is marked for the modified
problem has been studied extensively [10], [21], [31]. For Newton-Kleinman method.
10 2 2 2
E 2.68 ×10 N/m where M = EI d2φ+C I d2ψ ∈H2(0,1).
I 1.64 ×10−9 m4 bdr d bdr
b We use R = 1 and define C by the tip position:
ρ 1.02087 kg/m
2
Cv 2 Ns/m
Cd 2.5 ×108 Ns/m2 w(1,t) = C[w(x,t) w˙(x,t)].
L 1 m
I 121.9748 kg m2 Then Q = C∗C. We also solved the control problem with
h
d .041 kg−1 Q=I.
TABLE III Let H2N ⊂ H be a sequence of finite-dimensional
TABLE OF PHYSICAL PARAMETERS. subspaces spanned by the standard cubic B-splines with a
uniform partition of [0,1] into N subintervals. This yields
an approximation in H2N × H2N [16, e.g.] of dimension
II. CONTROL OF EULER-BERNOULLI BEAM n = 4N. This approximation method yields a sequence
of solutions to the algebraic Riccati equation that converge
In this section, we consider a Euler-Bernoulli beam strongly to the solution to the infinite-dimensonal Riccati
clamped at one end (r = 0) and free to vibrate at the other equation corresponding to the original partial differential
end (r = 1). Let w(r,t) denote the deflection of the beam equation description [4], [22].
from its rigid body motion at time t and position r. The
deflection is controlled by applying a torque u(t) at the III. LOW RANK APPROXIMATIONS TO LYAPUNOV
clamped end (r = 0). We assume that the hub inertia I FUNCTION
¨ h
is much larger than the beam inertia I so that I θ ≈ u(t).
b h Tables IV,V show the effect of varying C on the number
Thepartial differential equation model with Kelvin-Voigt and d
of iterations required for convergence. Larger values of C
viscous damping is d
(i.e. smaller values of α) leads to a decreasing number of
iterations. Small values of C lead to a large number of
ρw (r,t)+C w (r,t)+... d
tt v t required iterations in each solution of a Lyapunov equation.
∂2 ρr
[C I w (x,t) +EI w (r,t)] = u(t), Recall that the CF-ADI algorithm used here starts with a
∂r2 d b rrt b rr I
h rank 1 initial estimate of the Cholesky factor and the rank of
with boundary conditions the solution is increased at each step. The efficiency of the
w(0,t) = 0 Cholesky-ADI method relies on the existence of a low-rank
w (1,t) = 0. approximation of the solution X to the Lyapunov equation.
r This is true of many other iterative algorithms to solve large
EI w (1,t)+C I w (1,t) = 0
b rr d b rrt Lyapunov equations.
∂ [EI w (r,t)+C I w (r,t)] = 0. Theorem 3.1: For any symmetric, positive semi-definite
∂r b rr d b rrt r=1 matrix X of rank n let µ ≥ µ ... ≥ µ be its eigenvalues.
1 2 n
The values of the physical parameters in Table II are as in For any rank k < n matrix Xk,
[2]. *X−X * µ
Let k 2 ≥ k+1.
# $ *X* µ
dw 2 1
H= w∈H2(0,1):w(0)= (0) = 0 Proof: See, for example, [13, Thm. 2.5.3]. !
dr Thus, the relative size of the largest and smallest eigen-
be the closed linear subspace of the Sobolev space H2(0,1) values determines the existence of a low rank approximation
2 Xk that is close to X, regardless of how this approximation
and define the state-space to be X = H×L (0,1) with state is obtained.
z(t) = (w(·,t), ∂ w(·,t)). A state-space formulation of the
∂t The ratio µ /µ is plotted for several values of C in
above partial differential equation problem is k+1 1 d
Figure 1. For larger values of C the solution X is closer to
d
d a low rank matrix than it is for smaller values of C .
x(t) = Ax(t)+Bu(t), d
dt A bound on the rank of the solution to a Lyapunov
where equation where A is symmetric is given in [25]. A tighter
0 I bound on the error in a low-rank approximation has been
A= obtained [1] in terms of the eigenvalues and eigenvectors of
EI 4 C I 4 C A. This bound is applicable to all diagonalizable matrices
− bd − d b d − v
ρ dr4 ρ dr4 ρ A. The bound for the case where the right-hand-side D has
and 0 rank one is as follows.
B= Theorem 3.2: [1, Thm 3.1] Let V be the matrix of right
r eigenvectors of A, and denote its condition number by κ(X).
Ih Denote the eigenvalues of A by λ1,...λn. There is a rank
with domain k < n approximation to X satisfying
d 2 2
dom(A)={(φ,ψ)∈X :ψ∈H, M(L)= M(L)=0} *X−X * ≤(n−k) δ (κ(V)*D* ) (8)
dr k 2 k+1 2
TABLE IV
where BEAM: : EFFECT OF CHANGING C (STANDARD NEWTON-KLEINMAN)
* * d
k 2
−1 )*λk+1−λj*
δk+1 = * ∗ * . Cv Cd α N.-K. It’s Lyapunov Itn’s CPU time
2Real(λk+1) j=1*λk+1 +λj* 2 1E04 1.5699 – – –
3E05 1.5661 – – –
The eigenvalues λi are ordered so that δi are decreasing. 4E05 1.5654 3 1620;1620;1620 63.14
In order to calculate this bound, all the eigenvalues and 1E07 1.5370 3 1316;1316;1316 42.91
eigenvectors of A are needed. 1E08 1.4852 3 744;744;744 18.01
If the condition number κ(V ) is not too large, for instance 5E08 1.3102 3 301;301;301 5.32
if A is normal, then δk+1/δ1 gives a relative error estimate of
the error µ /µ in the approximation X . This estimate TABLE V
k+1 1 k BEAM: EFFECT OF CHANGING C (MODIFIED NEWTON-KLEINMAN)
is observed in [1] to estimate the actual decay rate of the d
eigenvalues of X quite well, even in cases where the decay C C α N.-K. It’s Lyapunov It’s CPU time
is very slow. v d
2 1E04 1.5699 2 – –
Consider a parabolic partial differential equation, with A- 3E05 1.5661 2 – –
matrix of the approximation is symmetric and negative defi- 4E05 1.5654 2 1620;1 24.83
nite. Then all the eigenvalues are real. A routine calculation 1E07 1.5370 2 1316;1 16.79
yields that 1E08 1.4852 2 744;1 7.49
5E08 1.3102 2 301;1 2.32
δk+1 ≈ λ1
δ1 λk
and so the rate of growth of the eigenvalues determines the [23]. A simple calculation yields that for all k,
accuracy of the low rank approximation. The accuracy of
the low-rank approximant can be quite different if A is non- k−1
symmetric. The solution X to the Lyapunov equation could δk = ) 1 2 ≈1.
4c
δ1 1+ v 2
be, for instance, the identity matrix [1], [25] in which case j=1 (wk−wj)
the eigenvalues of X do not decay. This indicates that the eigenvalues of the Lyapunov solution
It is observed through several numerical examples in [1] do not decay.
that it is not just complex eigenvalues that cause the decay
rate to slow down, but the dominance of the imaginary parts IV. CONCLUSIONS
of the eigenvalues as compared to the real parts, together with The results indicate a problem with solving the algebraic
absence of clustering in the spectrum of A. This effect was Riccati equation for systems with light damping where the
observed in the beam equation. Figure 2 shows the change eigenvalues are not decaying rapidly. Although better choice
in the spectrum of A as C is varied. Essentially, increasing
d of the ADI parameters might help convergence, the fact
C increases the angle that the spectrum makes with the
d that the low rank approximations to the solution converge
imaginary axis. Note that the eigenvalues do not cluster. As very slowly when the damping is small limits the achievable
the damping is decreased, the dominance of the imaginary improvement. This may have consequences for control of
parts of the eigenvalues of A increases and the decay rate of coupled acoustic-structure problems where the spectra are
the eigenvalues of X slows. The Cholesky-ADI calculations closer to those of hyperbolic systems than those of parabolic
for the beam with C = 0 did not terminate on a solution
d systems.
Xk with k < n.
These observations are supported by results in the theory
of control of partial differential equations. If C > 0 the
d -./0123450+637.8+98:+;<=)+>.7?+23:@.1/+A
B
! +
semigroup for the original partial differential equation is #!
parabolic and the solution to the Riccati equation converges !#!
uniformly in operator norm [18, chap.4]. However, if C = 0, #!
d
!$!
the partial differential equation is hyperbolic and only strong #!
convergence of the solution is obtained [19]. Thus, one might !%!
#!
expect a greater number of iterations in the Lyapunov loop #
! !&!
+,+*#!
to be required as C is decreased. Consider the bound (8) !
d
!"!
for the beam with C = 0, and an eigenfunction basis for #!
d
the approximation so that the eigenvalues of A are the exact !)!
2 #!
eigenvalues. Defining c = C /ρ, c = EI /ρandindicating '
v v b "C#!
!(! '
the roots of 1 + cos(a)cosh(a) by a , the eigenvalues are #! #C#!
k
(
#C#!
!'!+
λ =−c /2±iω #!
k v k ! " #! #" $! $" %! %" &!
*
where +
2 2 Fig. 1. Eigenvalue ratio for solution to beam ARE for varying C
ωk = ca −c /4 d
k v
no reviews yet
Please Login to review.
