- Email: [email protected]
MACSYMA -Generated closed-form solutions to some matrix Riccati equations
wvzww-Generated Closed-Form Solutions to Some Matrix Riccati Equations Ralph M. Wilcox* Hughes Aircraft Company Electra-Optical Data Systems Group El Segundo, Califwia 90245 and Leo P. Harten Massachusetts Znstitute of Technology Cambridge, Massachusetts 02139
Transmitted by Harriet Kagiwada
ABSTRACT Closed-form solutions to some continuous-time matrix Riccati equations have been generated by employing MACSYMA, the symbolic manipulation program developed at MIT. These equations, whose solution by hand would be prohibitive, occur in a description of the performance of imaging infrared trackers. Solutions have been obtained for both secondorder and third-order matrices (corresponding, respectively, to the application to the target of a white-noise or a Poisson-type stochastic force), and for the cases in which the variance matrix is initially either zero or singular.
1.
INTRODUCTION
Matrix Riccati equations occur in many engineering and scientific applications [l-11], particularly in control, estimation, and scattering theories. Thus
*This work has been supported, in part, by the Night Vision and Electrc+Optics Laboratory through U.S. Army Research Office contract DAAG2976Da100, awarded to Hughes Aircraft Company on 26 June 1981. APPLIED MATHEMATZCS AND COMPUTATION 14:149-166 Q R. M. Wilcox et al., 1984
(1984)
149
150
RALPH M. WILCOX AND LEO P. HARTEN
solutions or techniques found useful in one problem may be significant for other applications as well. Because of the complexity of the coupled nonlinear differential equations, closed-form solutions have been for the most part confined to the scalar case. A notable exception is the work of Nishimura and Kano [ll], where by setting M = 0 and choosing A and N in a particular way in Equation (1) below, these authors obtained a periodic solution for the 2 x 2 matrix case. The matrix Riccati equations solved in this paper have occurred as part of an effort to analyze the performance of imaging R trackers [12]. Four problems are defined in Section 2, and solved in Sections 3, 4, 5, and 6. All problems correspond to a model where the image of a maneuvering target is imaged upon a rectangular array of sensor pixels. Although in this model the target motion is confined to one dimension, the model may be extended to two dimensions. When this is done, it is found that, for the initial conditions employed, the corresponding Riccati equation decouples into two simpler Riccati equations of the form solved in this paper. Problems I and II are for the case where the target is acted upon by a white-noise random force, while problems III and IV are for the case where the target is subject to random step accelerations (or decelerations) occurring at random times whose distribution is Poisson.’ Problems I and III correspond to zero initial conditions (the target’s initial position is known exactly), while problems II and IV correspond to singular initial conditions (the target’s initial position is completely unknown). In all problems the measurement model has been linearized, thus yielding the Riccati form of the variance matrix equations. The difficulty of obtaining analytical solutions to the matrix Riccati equations has been recognized by Sage and White who, after obtaining a solution for the scalar case, state: 2 As is readily apparent, the effort involved in solving this simple Riccati equation is such as to suggest that the analytical solution of a Riccati equation of higher order than the first would be prohibitive.
These authors evidently were not considering the possibility of employing a computer to perform the required manipulations. This paper presents results which were obtained and verified by MACSYMA,the general-purpose user-extensible algebraic manipulator developed by the Mathlab Group at the MIT Laboratory for Computer Science.3 ‘See p. 84 of Reference [8]. ‘See p. 217 of Reference [7]. 3See Reference [ 131 and Reference [ 141
MACSYMA-GCT~WT~~~~ Closed-Form Solutions 2.
151
GENERAL DISCUSSION The basic Riccati equation solved in this work is of the form
dV dt
-=AV+VAT+M-VNV,
(1)
where V is the variance matrix to be solved for as a function of time t; A, M, and iV are constant matrices; and AT is the transpose of the matrix A. We present the solutions to Equation (1) for the four problems labeled I, II, III, and IV in the following four sections. Problem I is defined mathematically by the second-order matrices
(2) (3)
(4 together with the initial condition
v(0) = 0.
(5)
Physically these equations arise for the optimalestimation problem of a rectangular target moving with one degree of freedom viewed by an imaging tracker [12]. It is assumed that the target is subjected to a white-noise force, and that the sensor pixel measurements are also subject to additive independent white noise. In obtaining the Riccati form (l), it has been necessary to linearize the measurement equations. Also, the special form of the matrices M and N in terms of the parameter f is due to a change of scale of the time variable t which makes time dimensionless. Problem II is the same as problem I except that the initial condition is specified by
v-l(o)
= 0,
where V- '(t ) is the inverse of the matrix V( t ).
(6)
152
RALPH M. WILCOX AND LEO P. HARTEN
Equations (5) and (6) are two extreme conditions, neither of which occurs in practice; Equation (5) corresponds to the case where the target’s position is known exactly at t = 0, while Equation (6) corresponds to the case where nothing is known about the target’s position at t = 0. Equation (6) implies that V( t ) is singular at time t = 0. The inverse variance V-‘(t) satisfies the following Riccati equation: dV-’ -= dt
-V-‘A-ATV-‘+N-V-‘MV-l,
0)
as is readily verified. Note that the system specified by Equations (7) and (6) may be obtained from the system specified by (1) and (5) by means of the following transpositions:
v-v-‘, A-
-AT,
(8)
Problem III is the same as problem I except that it is defined by the third-order matrices
(10)
(11)
The only physical difference between this problem and problem I is that instead of being subjected to a white-noise force, the target is acted upon by a Poisson process consisting of random step accelerations (or decelerations) at random times. The time scaling defining the parameter f in Equations (10) and (11) is also different from that used for Equations (3) and (4). Problem IV is the same as problem III except that we use the initial condition (6) instead of (5).
153
M.MzsYMA-GeneratedClosed-Form Solutions
The solution to these problems is based upon a well-known theorem which converts the problem of solving a matrix Riccati equation into the problem of solving a linear matrix differential equation in twice as many variables.4 Thus let
v=
XY-‘,
02)
where X and Y are square matrices of the same dimension as V satisfying the linear system e=AX+MY, dt 03)
subject to the initial conditions
x(0) = 0,
(14)
Y(O)=Z. Then it is easily verified that V satisfies Equation (1) subject to the initial condition (5). Equations (13) and (14) may also be written in the single matrix form dZ
dt=BZ,
z(o,=(;)~
(16)
where Z and B are extended matrices defined by
ZE
(x1 Y ’
‘-($
_;T)*
4See e.g. p. 159 of Reference [2] or p. 262 of Reference [9].
07)
(1%
RALPH M. WILCOX AND LEO P. HARTEN
154 3.
SOLUTION
TO PROBLEM
I
Solution of this problem, though straightforward, would be extremely tedious for a human to perform by hand. Consequently we have written a batch program in MACSYMA which performs all steps required to obtain the solution following the procedure outlined above, and which verifies the result by simplifying to zero the difference of the right- and left-hand sides of the Riccati equation (1) using the solution. The differential equation for this problem assumes the form
dXl1 dt=
dXE2 dt=
X 21> X
22’
dX *=fy,,,
dX22
7
=
fr223
(19) dY,l
x,1
x=7'
dY,,
dt d%
dt-
x,2
=f’ - _ y
dY22 -=
dt
11,
-y 12
subject to the initial conditions Y,,=Y,=l, Xl, = x,2 = x2,=x22
=
Y,, = Y,, = 0.
(20)
MACSYMA solves such a system of coupled differential equations by using a routine DESOLVE (contributed by Richard Bogen and documented in Reference [13]) which uses the Laplace-transform method and then inverts the transforms. LISP
MAcsYMA-kwrated
155
Closed-Form !?ok&nas
After the X and Y matrices are separated out from the soIution of Equations (19) and (20), the Y matrix is inverted and multiplied from the left by the X matrix in accordance with Equation (12), and the result is simplified, yielding the following elements of V:
v,&)=v,l(t)=~ v
(tJ
_
f
22
f
cosh(@t)
- co+&) W
Wd
’
sinh(&t)+sin(Gt) w ’ J3T
Gw
where D(t)=d[2+cosh(@t)+cos(&t)].
(22)
Computer plots of these curves are presented in Figure 1 for the case f = 1;
1
2
5
4
5
FIG. 1. Normalized solution ( f = 1) for problem I. The solid line is V,,(t); V12(t); the dotted line is V,(t).
6
7
the dashed line is
RALPH M. WILCOX AND LEO P. HARTEN
156
the linear scaling of the components of V with fallows the results to be found for other values of 5 From the relative simplicity of the resulting solution the reader may be deceived into believing that a solution by hand would not really be very tedious to perform. This is decidedly not the case, since much simplification of long expressions was required before obtaining the highly economical form of the final result. The simplicity of the result obtained does, however, suggest the possibility that more elegant methods of solution may exist. The MACSYMA batch program verifies that V indeed satisfies the coupled nonlinear differential equations specified by Equations (l)-(4):
dKl=2V I2
dt
-- wJ2 f
-=V,--
dvn
v,,v,,
dt
f
’
’
cw
It is readily shown that the solution vanishes at t = 0, as required, while for large t it approaches the steady-state solution V&)=&f,
(244
vlz(~)=vzlw
V,(m)=fif.
4.
=
f9
@4
(24c)
SOLUTION TO PROBLEM II
To obtain the solution for V-i satisfying Equations (7) and (6), we note that the transpositions (8) are equivalent to the same initial-value problem (13), (14) as before, but with the role of X and Y interchanged so that instead of Equation (17) we now have
z=
(f;,1
w
MACSYMA-GSWW~~ Closed-Form Solutions
157
instead of Equation (14) we have
x(0) = I, (26)
Y(0) = 0, and instead of Equation (12) we have v-1
=yx-1.
(27)
Thus the solution for V-i is obtained by minor modifications to the batch program used to solve problem I. However, an even easier method of solution is possible by comparing in detail the differential equations satisfied by the elements of V and V-‘. Defining
dull
-
dt
dU,2
- dt
1 =- -
f
= - f
f(Q2, Ul,U,
dU2s = - f(Q2 dt
(294 Ull,
@W
- 2U12.
(294
-
Since Equation (29) may be obtained from Equation (23) by making the substitutions
(30)
the solution for U may be found by making these same substitutions in
RALPH M. WILCOX AND LEO P. HARTEN
158
‘0
I 1
I 2
FIG.2. Normalized solution ( f= is VIZ(t); the dotted line is V,(t).
-i3-
1) for problem
-
I 4
II. The solid line is V,,(t); the dashed line
Equations (21) and (22). When this is done, and MACSYMA inverts U to obtain V, it is found, after considerable further simplification, that V again is represented by Equation (21) and that instead of Equation (22), o(t) is now to be defined by D(t)=g[
-2+
cosh(&t)+cos(\/;Zt)].
(31)
Plots of these solutions for the case f = 1 are given in Figure 2. The same program used to verify solution I also verifies that this solution satisfies Equation (23). For small values of t, MACSYMA performed a Taylor-Laurent expansion about t = 0 to yield the leading behavior of V:
v,, = 4f/t
+ *. * ,
V,,=V,,=6f/t2+
v, = 12f/t3
+ ** - .
(32a) ..a,
WW (324
MACSYMA-~~~~~~Closed-Form Solutions
159
For large values of t, it is again easily found, either by using MACSYMA or by hand calculation, that the steady-state limiting behavior is again given by Equations (24).
5.
SOLUTION
TO PROBLEM
III
The solution of problem III proceeds by modifying the batch programs used on problem I; the main differences are the replacement of the matrix definitions (2), (3), and (4) by Equations (9), (lo), and (ll), and the use of different trigonometric, hyperbolic, and exponential identities for the simplification of results. We obtain the following results:
(334 (3W (33c)
(33d)
@W
where the quantities a(t),
b(t), g(t), h(t), k(t), m(t), and d(t)
are defined
by a(t) = sinh2t + [4+cos(fit)]sinht
+4cos(ifit)sinhit,
b(t)=fisin(fit)cosht
+46sin($fit)cosh+t,
g(t)=fisin(fit)sinht
+2&sin(+Gt)sinhgt,
h(t)= k(t) m(t)= d(t)
cosh2t + [2 - cos(&)]cosht
- 2cos(+&t)cosh$,
(34)
= sinh2t - 2cos(fit)sinht, 4sinht
- 8cos($fit)sinhbt,
= cosh2t + [8+2cos(&)]coshl+16cos(+&)cosh+t
+9.
RALPH M. WILCOX AND LEO P. HARTEN
160
0
1
2
3
4
5
6
FIG. 3. Normalizedsolution(f= 1) for problem III. The solid line is VII(t);the long-dashed line is V,,( t ); the shortdashed line is V,,( t ); the dotted line is V,( t ); the dash-dot line is V,(t); the dashdotdot line is V,(t).
Computer plots of the six distinct matrix elements are given in Figure 3 for thecasef= 1. The MACSYMA batch program verifies that V satisfies the following coupled nonlinear differential equations defined by Equations (l), (9), (lo), and (11):
dV,, -_=2v,,_~, dt VllVl2 f ’
dV,,
-=v13+vv22--
dt
(354
dh.3 -=V,--.--
v,2v,3
f'
dt
d&s
dt
-f_
W3)"
f.
(354
Wf)
~csnm-Generated
161
Closed-Form Solutiot~
The amount
of work involved in performing this verification is almost as great as that involved in deriving the solution. It is readily shown that the solution vanishes at t = 0, as required, while for large t it approaches the steady-state solution
v,,(4 = 2f >
(364
v,2(4 =%W = 2fT
(36b)
v,,b) =Gb4 = f?
(364
%&4
6.
SOLUTION
=
3f,
(36d)
&3W=YB(4=2f,
(364
v,(co)
(360
= 2f.
TO PROBLEM
IV
The solution of problem IV is similar to that Equations (7), (9), (lo), (ll), and (28) it follows that
of problem
II. From
dv,, = -1 - f(Q2,
-
dt
dun =
- q, - f
q$J,
9
(3%)
dU,,=
- q, - f
q&J,
3
(374
d%2 -=
-
fwi3)2~
dU23 -=
- q,
d&3 -=
-2u,-
-
-
dt
dt
dt
dt
dt
Comparing
(374
f
these equations
w2
-
- u, -
with Equation
f q&$Jz,,
f(UJ2. (35), it is seen that Equation
(374 (374
Wf) (37)
RALPH M. WILCOX AND LEO P. HARTEN
162
may be obtained from Equation (35) by making the following substitutions:
The solution U is thus found by making these substitutions into Equation (33). When this is done, and MACSYMA inverts U to yield V, it is found, after very considerable further simplification, that the elements of V are given by
v
(Q= 11
fm+Ptt) d(t)
(394 ’
w$
tw (394 (394 v
(t)= 33
f4(t)-Po) d(t)
’
(39f)
~cwm-Generated
163
Closed-Form Solutions
where the quantities 9(t), p(t), w(t), q(t) = 4cosh2t
x(t), y(t), x(t), and d(t) are defined by
- [4cos(&)+16]cosht
+16cos($&)cosh&t, (40a)
p(t) = 4fisinhtsin(&t)
- 16fisinh&in($&t),
(4Ob)
w(t)=4sinh2t+[4cos(&t)-8]sinht--8sinhdtcos(~fit),
(4Oc)
x(t)=&coshtsin(fit)-8ficoshbtsin(+fit), y(t) = 2cosh2t
+4coshtcos(&)
- 6,
z(t) = 24 - 16cosh$cos(;&) d(t)
(4W (40(5)
- 8cosh(t),
= 32sinhbtcos($fit)+2sinh2t
(4of)
- sinht[4cos(&t)+16]. (4Og)
In performing the simplifications, the determinant obtained while inverting U was found to contain the determinant of the original problem as a double factor. The batch program written to verify solution III also verifies this one. For small values of t, MACSYMA performed a Taylor-Laurent expansion about t = 0 to yield the leading behavior of V:
v,, = 9f/t
+ *. *
(414
)
V,,=V,,=36f/t2+
....
(4W
V,,=V,,=eof/t”+
**.,
(41c)
v,
= 192f/P
+ **. )
V,=V3,=360f/t4+
v, = 720f/P
(414 ..*,
+ *. . .
(41c) @1f)
Note that in this case d(t) has a Taylor expansion starting with the tg term; a hand calculation of the first nonvanishing derivative about t = 0
RALPH M. WILCOX AND LEO P. HARTEN
164
-. ‘0
I 1
I 2
I 3
-I --a
4
FIG. 4. Normalized solution (f= 1) for problem IV. The solid line is V,,(t); the long-dashed line is V,,(t); the shortdashed line is V,,(t); the dotted line is V.(t); the dash-dot line is V,(t); the dashdotdot line is V,(t). As explained in Section 6, we plotted the Taylor-Laurent series for t < 0.3 and the analytic solutions for t > 0.3 to obtain accurate results.
would involve a great deal of time; and owing to the high order zero of d(t), it would be unreliable to numerically establish the Taylor series by differences, because the form of d(t) given in Equation (40g) is unsuitable for computation in single-precision arithmetic when t < 0.3, despite its relatively trivial form. MACSYMA retains the rational numbers generated in the Taylor series to infinite precision, and thus we can be certain that the answer is correct in having perfect cancellation of the coefficients of the to through t8 terms. We thus plotted the first two terms of the Taylor-Laurent series for t.< 0.3 and then matched onto the analytical form for t > 0.3 to get the accurate results in Figure 4, and we again took f = 1. For large values of t, it is again easily found, either by using MACSYMA or by hand calculation, that the steady-state limiting behavior is again given by Equations (36). 7.
FURTHER
REMARKS
One of the advantages of computer calculations is that programs written for one purpose can by small modification be used for other purposes as well.
r+rAcsym-Generated Closed-Form Solutions
165
Consequently, we will not be surprised if modifications of our programs obtain additional closed-form solutions as well. The use of general and powerful computer algebra programs, such as MACSYMA, goes far beyond the derivation of matrix Riccati closed-form solutions. Since the MACSYMAmanual [ 131 details the available capabilities, we merely mention the basic process by which we generate such programs: We use the interactive line-at-a-time mode in which MACSYMAevaluates and then simplifies expressions as they are typed in. When we see that the operations perform in the desired way (to get through steps of a larger organized routine), we start combining the intermediate steps in simple ways and continue to develop more powerful abilities. We reach a stage where the whole program prints out only the results which we wish to use (since many of the intermediate steps yield large results which we are not really interested in viewing), or we get stuck and try another approach. There may be several hours of such interactive use, followed by a batch program which runs the desired case in a few minutes without further interaction. In combination with MACSYMA’Swide-ranging routines, we can also examine a variety of techniques to apply to particular problems, and thus select the ones which appear to be the most general, or the most efficient, depending on the scope of the task that we wish to carry out. We also note the utility of exact versus numerical solutions, particularly in the context of singular problems; the singular initial condition is well known mathematically, but is hard to pose numerically without making some transformation to a reciprocal quantity which is zero initially. Yet analytic results can be obtained in either case with the aid of MACSYMA. One of the authors (RMW) thanks his colleagues Larry Rubin, Carleton Nealy, and Richard Frey of Hughes Aircraft Company, and Stan Rodak and Frank Shields of the Night Vision Electra-Optical Laboratory, for their interest and encouragement of this work. The research presented here was carried out with the aid of MACSYMA, a large symbolic manipulation program developed at the MIT Laboratory for Computer Science and supported by the National Aeronautics and Space Administration under grant NSG 1323, by the Office of Naval Research under grant NOOO14-77-C-0641, by the U.S. Department of Energy under grant ET-78-C-02-4687, and by the U.S. Air Force under grant F49620- 79-C-020. REFERENCES 1 2
W. T. Reid, Riccati Differential Equations, Academic, New York, 1972. R. Bellman, Introduction to the Mathematical Theory of Control Processes, Vol. 1, Academic, New York, 1967.
166 3 4 5 6 7 8 9 10 11 12
13 14
RALPH M. WILCOX
AND LEO P. HARTEN
Ft. Bellman, Methods of Non-linear Analysis, Vol. 2, Academic, New York, 1967. It. S. Bucy and P. D. Joseph, Filtering fo* Stochastic Processes with Applications to Guidance, Interscience, New York, 1968. K. J. Astrom, Introduction to Stochastic Control Theory, Academic, New York, 1970. S. Bamett, Matrices in Control Theory, Van Nostrand Reinhold, London, 1971. A. P. Sage and C. C. White, III, Optimum Systems Control, 2nd ed., Prentice-Hall, Englewood Cliffs, N.J., 1977. T. P. McCarty, Stochustic Systems and State Estimation, Wiley, New York, 1974. Z. Schuss, Theory and Applications of Stochastic Dffmential Equations, Wiley, New York, 1978. R. M. Redheffer, On the relation of transmission line theory to scattering and transfer, 1. Math. and Phys. 41: 1-41 (1962). T. Nishimura and H. Kano, Periodic oscillations of matrix Riccati equations in time-invariant systems, IEEE Trans. Automat. Control AC-25: 749-755 (1980). R. L. Frey, C. D. Nealy, L. M. Rubin, and R. M. Wilcox, Analytic performance modeling of imaging trackers, Preliminary Report submitted to Night Vision ElectroOpticaI Laboratories, 5 March 1982. R. Bogen et al., MACSYMAReferatce Manual, The Mathlab Group, Lab. for Computer Sci., Mass. Inst. of Tech., Cambridge, MA 02139, 1977. R. Pavelle, M. Rothstein, and J. Fitch, Computer algebra, Sci. Amer. 245(6): 136-152 (1981).
Transmitted by Harriet Kagiwada
ABSTRACT Closed-form solutions to some continuous-time matrix Riccati equations have been generated by employing MACSYMA, the symbolic manipulation program developed at MIT. These equations, whose solution by hand would be prohibitive, occur in a description of the performance of imaging infrared trackers. Solutions have been obtained for both secondorder and third-order matrices (corresponding, respectively, to the application to the target of a white-noise or a Poisson-type stochastic force), and for the cases in which the variance matrix is initially either zero or singular.
1.
INTRODUCTION
Matrix Riccati equations occur in many engineering and scientific applications [l-11], particularly in control, estimation, and scattering theories. Thus
*This work has been supported, in part, by the Night Vision and Electrc+Optics Laboratory through U.S. Army Research Office contract DAAG2976Da100, awarded to Hughes Aircraft Company on 26 June 1981. APPLIED MATHEMATZCS AND COMPUTATION 14:149-166 Q R. M. Wilcox et al., 1984
(1984)
149
150
RALPH M. WILCOX AND LEO P. HARTEN
solutions or techniques found useful in one problem may be significant for other applications as well. Because of the complexity of the coupled nonlinear differential equations, closed-form solutions have been for the most part confined to the scalar case. A notable exception is the work of Nishimura and Kano [ll], where by setting M = 0 and choosing A and N in a particular way in Equation (1) below, these authors obtained a periodic solution for the 2 x 2 matrix case. The matrix Riccati equations solved in this paper have occurred as part of an effort to analyze the performance of imaging R trackers [12]. Four problems are defined in Section 2, and solved in Sections 3, 4, 5, and 6. All problems correspond to a model where the image of a maneuvering target is imaged upon a rectangular array of sensor pixels. Although in this model the target motion is confined to one dimension, the model may be extended to two dimensions. When this is done, it is found that, for the initial conditions employed, the corresponding Riccati equation decouples into two simpler Riccati equations of the form solved in this paper. Problems I and II are for the case where the target is acted upon by a white-noise random force, while problems III and IV are for the case where the target is subject to random step accelerations (or decelerations) occurring at random times whose distribution is Poisson.’ Problems I and III correspond to zero initial conditions (the target’s initial position is known exactly), while problems II and IV correspond to singular initial conditions (the target’s initial position is completely unknown). In all problems the measurement model has been linearized, thus yielding the Riccati form of the variance matrix equations. The difficulty of obtaining analytical solutions to the matrix Riccati equations has been recognized by Sage and White who, after obtaining a solution for the scalar case, state: 2 As is readily apparent, the effort involved in solving this simple Riccati equation is such as to suggest that the analytical solution of a Riccati equation of higher order than the first would be prohibitive.
These authors evidently were not considering the possibility of employing a computer to perform the required manipulations. This paper presents results which were obtained and verified by MACSYMA,the general-purpose user-extensible algebraic manipulator developed by the Mathlab Group at the MIT Laboratory for Computer Science.3 ‘See p. 84 of Reference [8]. ‘See p. 217 of Reference [7]. 3See Reference [ 131 and Reference [ 141
MACSYMA-GCT~WT~~~~ Closed-Form Solutions 2.
151
GENERAL DISCUSSION The basic Riccati equation solved in this work is of the form
dV dt
-=AV+VAT+M-VNV,
(1)
where V is the variance matrix to be solved for as a function of time t; A, M, and iV are constant matrices; and AT is the transpose of the matrix A. We present the solutions to Equation (1) for the four problems labeled I, II, III, and IV in the following four sections. Problem I is defined mathematically by the second-order matrices
(2) (3)
(4 together with the initial condition
v(0) = 0.
(5)
Physically these equations arise for the optimalestimation problem of a rectangular target moving with one degree of freedom viewed by an imaging tracker [12]. It is assumed that the target is subjected to a white-noise force, and that the sensor pixel measurements are also subject to additive independent white noise. In obtaining the Riccati form (l), it has been necessary to linearize the measurement equations. Also, the special form of the matrices M and N in terms of the parameter f is due to a change of scale of the time variable t which makes time dimensionless. Problem II is the same as problem I except that the initial condition is specified by
v-l(o)
= 0,
where V- '(t ) is the inverse of the matrix V( t ).
(6)
152
RALPH M. WILCOX AND LEO P. HARTEN
Equations (5) and (6) are two extreme conditions, neither of which occurs in practice; Equation (5) corresponds to the case where the target’s position is known exactly at t = 0, while Equation (6) corresponds to the case where nothing is known about the target’s position at t = 0. Equation (6) implies that V( t ) is singular at time t = 0. The inverse variance V-‘(t) satisfies the following Riccati equation: dV-’ -= dt
-V-‘A-ATV-‘+N-V-‘MV-l,
0)
as is readily verified. Note that the system specified by Equations (7) and (6) may be obtained from the system specified by (1) and (5) by means of the following transpositions:
v-v-‘, A-
-AT,
(8)
Problem III is the same as problem I except that it is defined by the third-order matrices
(10)
(11)
The only physical difference between this problem and problem I is that instead of being subjected to a white-noise force, the target is acted upon by a Poisson process consisting of random step accelerations (or decelerations) at random times. The time scaling defining the parameter f in Equations (10) and (11) is also different from that used for Equations (3) and (4). Problem IV is the same as problem III except that we use the initial condition (6) instead of (5).
153
M.MzsYMA-GeneratedClosed-Form Solutions
The solution to these problems is based upon a well-known theorem which converts the problem of solving a matrix Riccati equation into the problem of solving a linear matrix differential equation in twice as many variables.4 Thus let
v=
XY-‘,
02)
where X and Y are square matrices of the same dimension as V satisfying the linear system e=AX+MY, dt 03)
subject to the initial conditions
x(0) = 0,
(14)
Y(O)=Z. Then it is easily verified that V satisfies Equation (1) subject to the initial condition (5). Equations (13) and (14) may also be written in the single matrix form dZ
dt=BZ,
z(o,=(;)~
(16)
where Z and B are extended matrices defined by
ZE
(x1 Y ’
‘-($
_;T)*
4See e.g. p. 159 of Reference [2] or p. 262 of Reference [9].
07)
(1%
RALPH M. WILCOX AND LEO P. HARTEN
154 3.
SOLUTION
TO PROBLEM
I
Solution of this problem, though straightforward, would be extremely tedious for a human to perform by hand. Consequently we have written a batch program in MACSYMA which performs all steps required to obtain the solution following the procedure outlined above, and which verifies the result by simplifying to zero the difference of the right- and left-hand sides of the Riccati equation (1) using the solution. The differential equation for this problem assumes the form
dXl1 dt=
dXE2 dt=
X 21> X
22’
dX *=fy,,,
dX22
7
=
fr223
(19) dY,l
x,1
x=7'
dY,,
dt d%
dt-
x,2
=f’ - _ y
dY22 -=
dt
11,
-y 12
subject to the initial conditions Y,,=Y,=l, Xl, = x,2 = x2,=x22
=
Y,, = Y,, = 0.
(20)
MACSYMA solves such a system of coupled differential equations by using a routine DESOLVE (contributed by Richard Bogen and documented in Reference [13]) which uses the Laplace-transform method and then inverts the transforms. LISP
MAcsYMA-kwrated
155
Closed-Form !?ok&nas
After the X and Y matrices are separated out from the soIution of Equations (19) and (20), the Y matrix is inverted and multiplied from the left by the X matrix in accordance with Equation (12), and the result is simplified, yielding the following elements of V:
v,&)=v,l(t)=~ v
(tJ
_
f
22
f
cosh(@t)
- co+&) W
Wd
’
sinh(&t)+sin(Gt) w ’ J3T
Gw
where D(t)=d[2+cosh(@t)+cos(&t)].
(22)
Computer plots of these curves are presented in Figure 1 for the case f = 1;
1
2
5
4
5
FIG. 1. Normalized solution ( f = 1) for problem I. The solid line is V,,(t); V12(t); the dotted line is V,(t).
6
7
the dashed line is
RALPH M. WILCOX AND LEO P. HARTEN
156
the linear scaling of the components of V with fallows the results to be found for other values of 5 From the relative simplicity of the resulting solution the reader may be deceived into believing that a solution by hand would not really be very tedious to perform. This is decidedly not the case, since much simplification of long expressions was required before obtaining the highly economical form of the final result. The simplicity of the result obtained does, however, suggest the possibility that more elegant methods of solution may exist. The MACSYMA batch program verifies that V indeed satisfies the coupled nonlinear differential equations specified by Equations (l)-(4):
dKl=2V I2
dt
-- wJ2 f
-=V,--
dvn
v,,v,,
dt
f
’
’
cw
It is readily shown that the solution vanishes at t = 0, as required, while for large t it approaches the steady-state solution V&)=&f,
(244
vlz(~)=vzlw
V,(m)=fif.
4.
=
f9
@4
(24c)
SOLUTION TO PROBLEM II
To obtain the solution for V-i satisfying Equations (7) and (6), we note that the transpositions (8) are equivalent to the same initial-value problem (13), (14) as before, but with the role of X and Y interchanged so that instead of Equation (17) we now have
z=
(f;,1
w
MACSYMA-GSWW~~ Closed-Form Solutions
157
instead of Equation (14) we have
x(0) = I, (26)
Y(0) = 0, and instead of Equation (12) we have v-1
=yx-1.
(27)
Thus the solution for V-i is obtained by minor modifications to the batch program used to solve problem I. However, an even easier method of solution is possible by comparing in detail the differential equations satisfied by the elements of V and V-‘. Defining
dull
-
dt
dU,2
- dt
1 =- -
f
= - f
f(Q2, Ul,U,
dU2s = - f(Q2 dt
(294 Ull,
@W
- 2U12.
(294
-
Since Equation (29) may be obtained from Equation (23) by making the substitutions
(30)
the solution for U may be found by making these same substitutions in
RALPH M. WILCOX AND LEO P. HARTEN
158
‘0
I 1
I 2
FIG.2. Normalized solution ( f= is VIZ(t); the dotted line is V,(t).
-i3-
1) for problem
-
I 4
II. The solid line is V,,(t); the dashed line
Equations (21) and (22). When this is done, and MACSYMA inverts U to obtain V, it is found, after considerable further simplification, that V again is represented by Equation (21) and that instead of Equation (22), o(t) is now to be defined by D(t)=g[
-2+
cosh(&t)+cos(\/;Zt)].
(31)
Plots of these solutions for the case f = 1 are given in Figure 2. The same program used to verify solution I also verifies that this solution satisfies Equation (23). For small values of t, MACSYMA performed a Taylor-Laurent expansion about t = 0 to yield the leading behavior of V:
v,, = 4f/t
+ *. * ,
V,,=V,,=6f/t2+
v, = 12f/t3
+ ** - .
(32a) ..a,
WW (324
MACSYMA-~~~~~~Closed-Form Solutions
159
For large values of t, it is again easily found, either by using MACSYMA or by hand calculation, that the steady-state limiting behavior is again given by Equations (24).
5.
SOLUTION
TO PROBLEM
III
The solution of problem III proceeds by modifying the batch programs used on problem I; the main differences are the replacement of the matrix definitions (2), (3), and (4) by Equations (9), (lo), and (ll), and the use of different trigonometric, hyperbolic, and exponential identities for the simplification of results. We obtain the following results:
(334 (3W (33c)
(33d)
@W
where the quantities a(t),
b(t), g(t), h(t), k(t), m(t), and d(t)
are defined
by a(t) = sinh2t + [4+cos(fit)]sinht
+4cos(ifit)sinhit,
b(t)=fisin(fit)cosht
+46sin($fit)cosh+t,
g(t)=fisin(fit)sinht
+2&sin(+Gt)sinhgt,
h(t)= k(t) m(t)= d(t)
cosh2t + [2 - cos(&)]cosht
- 2cos(+&t)cosh$,
(34)
= sinh2t - 2cos(fit)sinht, 4sinht
- 8cos($fit)sinhbt,
= cosh2t + [8+2cos(&)]coshl+16cos(+&)cosh+t
+9.
RALPH M. WILCOX AND LEO P. HARTEN
160
0
1
2
3
4
5
6
FIG. 3. Normalizedsolution(f= 1) for problem III. The solid line is VII(t);the long-dashed line is V,,( t ); the shortdashed line is V,,( t ); the dotted line is V,( t ); the dash-dot line is V,(t); the dashdotdot line is V,(t).
Computer plots of the six distinct matrix elements are given in Figure 3 for thecasef= 1. The MACSYMA batch program verifies that V satisfies the following coupled nonlinear differential equations defined by Equations (l), (9), (lo), and (11):
dV,, -_=2v,,_~, dt VllVl2 f ’
dV,,
-=v13+vv22--
dt
(354
dh.3 -=V,--.--
v,2v,3
f'
dt
d&s
dt
-f_
W3)"
f.
(354
Wf)
~csnm-Generated
161
Closed-Form Solutiot~
The amount
of work involved in performing this verification is almost as great as that involved in deriving the solution. It is readily shown that the solution vanishes at t = 0, as required, while for large t it approaches the steady-state solution
v,,(4 = 2f >
(364
v,2(4 =%W = 2fT
(36b)
v,,b) =Gb4 = f?
(364
%&4
6.
SOLUTION
=
3f,
(36d)
&3W=YB(4=2f,
(364
v,(co)
(360
= 2f.
TO PROBLEM
IV
The solution of problem IV is similar to that Equations (7), (9), (lo), (ll), and (28) it follows that
of problem
II. From
dv,, = -1 - f(Q2,
-
dt
dun =
- q, - f
q$J,
9
(3%)
dU,,=
- q, - f
q&J,
3
(374
d%2 -=
-
fwi3)2~
dU23 -=
- q,
d&3 -=
-2u,-
-
-
dt
dt
dt
dt
dt
Comparing
(374
f
these equations
w2
-
- u, -
with Equation
f q&$Jz,,
f(UJ2. (35), it is seen that Equation
(374 (374
Wf) (37)
RALPH M. WILCOX AND LEO P. HARTEN
162
may be obtained from Equation (35) by making the following substitutions:
The solution U is thus found by making these substitutions into Equation (33). When this is done, and MACSYMA inverts U to yield V, it is found, after very considerable further simplification, that the elements of V are given by
v
(Q= 11
fm+Ptt) d(t)
(394 ’
w$
tw (394 (394 v
(t)= 33
f4(t)-Po) d(t)
’
(39f)
~cwm-Generated
163
Closed-Form Solutions
where the quantities 9(t), p(t), w(t), q(t) = 4cosh2t
x(t), y(t), x(t), and d(t) are defined by
- [4cos(&)+16]cosht
+16cos($&)cosh&t, (40a)
p(t) = 4fisinhtsin(&t)
- 16fisinh&in($&t),
(4Ob)
w(t)=4sinh2t+[4cos(&t)-8]sinht--8sinhdtcos(~fit),
(4Oc)
x(t)=&coshtsin(fit)-8ficoshbtsin(+fit), y(t) = 2cosh2t
+4coshtcos(&)
- 6,
z(t) = 24 - 16cosh$cos(;&) d(t)
(4W (40(5)
- 8cosh(t),
= 32sinhbtcos($fit)+2sinh2t
(4of)
- sinht[4cos(&t)+16]. (4Og)
In performing the simplifications, the determinant obtained while inverting U was found to contain the determinant of the original problem as a double factor. The batch program written to verify solution III also verifies this one. For small values of t, MACSYMA performed a Taylor-Laurent expansion about t = 0 to yield the leading behavior of V:
v,, = 9f/t
+ *. *
(414
)
V,,=V,,=36f/t2+
....
(4W
V,,=V,,=eof/t”+
**.,
(41c)
v,
= 192f/P
+ **. )
V,=V3,=360f/t4+
v, = 720f/P
(414 ..*,
+ *. . .
(41c) @1f)
Note that in this case d(t) has a Taylor expansion starting with the tg term; a hand calculation of the first nonvanishing derivative about t = 0
RALPH M. WILCOX AND LEO P. HARTEN
164
-. ‘0
I 1
I 2
I 3
-I --a
4
FIG. 4. Normalized solution (f= 1) for problem IV. The solid line is V,,(t); the long-dashed line is V,,(t); the shortdashed line is V,,(t); the dotted line is V.(t); the dash-dot line is V,(t); the dashdotdot line is V,(t). As explained in Section 6, we plotted the Taylor-Laurent series for t < 0.3 and the analytic solutions for t > 0.3 to obtain accurate results.
would involve a great deal of time; and owing to the high order zero of d(t), it would be unreliable to numerically establish the Taylor series by differences, because the form of d(t) given in Equation (40g) is unsuitable for computation in single-precision arithmetic when t < 0.3, despite its relatively trivial form. MACSYMA retains the rational numbers generated in the Taylor series to infinite precision, and thus we can be certain that the answer is correct in having perfect cancellation of the coefficients of the to through t8 terms. We thus plotted the first two terms of the Taylor-Laurent series for t.< 0.3 and then matched onto the analytical form for t > 0.3 to get the accurate results in Figure 4, and we again took f = 1. For large values of t, it is again easily found, either by using MACSYMA or by hand calculation, that the steady-state limiting behavior is again given by Equations (36). 7.
FURTHER
REMARKS
One of the advantages of computer calculations is that programs written for one purpose can by small modification be used for other purposes as well.
r+rAcsym-Generated Closed-Form Solutions
165
Consequently, we will not be surprised if modifications of our programs obtain additional closed-form solutions as well. The use of general and powerful computer algebra programs, such as MACSYMA, goes far beyond the derivation of matrix Riccati closed-form solutions. Since the MACSYMAmanual [ 131 details the available capabilities, we merely mention the basic process by which we generate such programs: We use the interactive line-at-a-time mode in which MACSYMAevaluates and then simplifies expressions as they are typed in. When we see that the operations perform in the desired way (to get through steps of a larger organized routine), we start combining the intermediate steps in simple ways and continue to develop more powerful abilities. We reach a stage where the whole program prints out only the results which we wish to use (since many of the intermediate steps yield large results which we are not really interested in viewing), or we get stuck and try another approach. There may be several hours of such interactive use, followed by a batch program which runs the desired case in a few minutes without further interaction. In combination with MACSYMA’Swide-ranging routines, we can also examine a variety of techniques to apply to particular problems, and thus select the ones which appear to be the most general, or the most efficient, depending on the scope of the task that we wish to carry out. We also note the utility of exact versus numerical solutions, particularly in the context of singular problems; the singular initial condition is well known mathematically, but is hard to pose numerically without making some transformation to a reciprocal quantity which is zero initially. Yet analytic results can be obtained in either case with the aid of MACSYMA. One of the authors (RMW) thanks his colleagues Larry Rubin, Carleton Nealy, and Richard Frey of Hughes Aircraft Company, and Stan Rodak and Frank Shields of the Night Vision Electra-Optical Laboratory, for their interest and encouragement of this work. The research presented here was carried out with the aid of MACSYMA, a large symbolic manipulation program developed at the MIT Laboratory for Computer Science and supported by the National Aeronautics and Space Administration under grant NSG 1323, by the Office of Naval Research under grant NOOO14-77-C-0641, by the U.S. Department of Energy under grant ET-78-C-02-4687, and by the U.S. Air Force under grant F49620- 79-C-020. REFERENCES 1 2
W. T. Reid, Riccati Differential Equations, Academic, New York, 1972. R. Bellman, Introduction to the Mathematical Theory of Control Processes, Vol. 1, Academic, New York, 1967.
166 3 4 5 6 7 8 9 10 11 12
13 14
RALPH M. WILCOX
AND LEO P. HARTEN
Ft. Bellman, Methods of Non-linear Analysis, Vol. 2, Academic, New York, 1967. It. S. Bucy and P. D. Joseph, Filtering fo* Stochastic Processes with Applications to Guidance, Interscience, New York, 1968. K. J. Astrom, Introduction to Stochastic Control Theory, Academic, New York, 1970. S. Bamett, Matrices in Control Theory, Van Nostrand Reinhold, London, 1971. A. P. Sage and C. C. White, III, Optimum Systems Control, 2nd ed., Prentice-Hall, Englewood Cliffs, N.J., 1977. T. P. McCarty, Stochustic Systems and State Estimation, Wiley, New York, 1974. Z. Schuss, Theory and Applications of Stochastic Dffmential Equations, Wiley, New York, 1978. R. M. Redheffer, On the relation of transmission line theory to scattering and transfer, 1. Math. and Phys. 41: 1-41 (1962). T. Nishimura and H. Kano, Periodic oscillations of matrix Riccati equations in time-invariant systems, IEEE Trans. Automat. Control AC-25: 749-755 (1980). R. L. Frey, C. D. Nealy, L. M. Rubin, and R. M. Wilcox, Analytic performance modeling of imaging trackers, Preliminary Report submitted to Night Vision ElectroOpticaI Laboratories, 5 March 1982. R. Bogen et al., MACSYMAReferatce Manual, The Mathlab Group, Lab. for Computer Sci., Mass. Inst. of Tech., Cambridge, MA 02139, 1977. R. Pavelle, M. Rothstein, and J. Fitch, Computer algebra, Sci. Amer. 245(6): 136-152 (1981).