# Precise Integration Method for Solving Noncooperative LQ Differential Game

Authors: Hai-Jun Peng; Sheng Zhang; Zhi-Gang Wu; Biao-Song Chen
Journal: Mathematical Problems in Engineering
 (2013)

---

## Abstract

The key of solving the noncooperative linear quadratic (LQ) differential game is to solve the coupled matrix Riccati differential equation. The precise integration method based on the adaptive choosing of the two parameters is expanded from the traditional symmetric Riccati differential equation to the coupled asymmetric Riccati differential equation in this paper. The proposed expanded precise integration method can overcome the difficulty of the singularity point and the ill-conditioned matrix in the solving of coupled asymmetric Riccati differential equation. The numerical examples show that the expanded precise integration method gives  more stable and accurate numerical results than the “direct integration method” and the “linear transformation method”.


---

## Body

## 1. Introduction

The research for the differential game is derived from the military field, such as the pursuit evasion problem, the spacecraft interception problem, and the cooperative or noncooperative problem. The emergence of the Isaacs’s monograph [1] gives birth to the differential games. Particularly, the differential game has been introduced into the economic field by Nash [2], and the equilibrium solution of the LQ differential game had been named by Nash-equilibrium solution. So far, many works for the LQ differential games have been done by researchers.The qualitative analysis about the LQ differential game, that is, the existence, stability of solutions, has been reported in many references. Reference [3] gives the series Nash solution of two-person nonzero-sum LQ differential game. Reference [4] discusses the Nash-equilibrium point and the existence, uniqueness of the algebraic Riccati equation for the closed LQ differential game. Reference [5] proposes the global existence of solutions to coupled Riccati differential equation (CRDE) in closed-loop Nash games. In contrast to the qualitative analysis, the numerical method for solving the LQ differential game problem to obtain the high performance solutions is still a challenging problem. Reference [6] founds that a special type of LQ differential game problem has the character of Hamiltonian system and applies the symplectic Runge-Kutta algorithm for solving this differential game problem. Reference [7] employs the Magnus integrators for solving N-player LQ differential games, and the CRDE which is reformulated as a linear ordinary equation is solved by Magnus integrators.From the above discussions, the key for solving LQ differential game is to solve the CRDE. Only a few problems can obtain their analytical solutions. It is necessary to employ numerical methods for solving the CRDE. The traditional numerical method is to employ an ordinary differential equation (ODE) solver for backward integration of the nonlinear CRDE. This kind of numerical methods is effective in most of cases. However, it may fail when the singularity point is encountered. Besides, the CRDE is also usually reformulated as a linear ODE such as [6, 7]. The advantage of this transformation method is that the linear ODE can be easier solved by comparing with the nonlinear ODE. However, the inverse matrix exists for the expression of the CRDE. Therefore, it restricts the application of this numerical method for many LQ differential game problems.In the field of LQ optimal control, the uncoupled Riccati differential equation (URDE) has been successfully solved by the precise integration method (PIM) [8–10]. Particularly, the PIM can obtain the machine accuracy solutions even for stiff problems. However, the difference between URDE and CRDE is obvious, and the CRDE is much difficult to be solved. There is no PIM for solving this CRDE from the LQ differential game problem. Therefore, it is partly inspired by the high accuracy of the PIM for solving the URDE, and we employ the adaptive parameter criterion to expand the PIM and apply it for solving the CRDE in this paper. Finally, the noncooperative LQ differential game problem can be solved by the proposed expanded PIM.


## 2. Noncooperative LQ Differential Game

### 2.1. Problem Description

For the general case withN-players, the controlled equation can be defined for the LQ differential game problem as follows:
(1)x˙(t)=Ax(t)+∑i=1NBiui(t),x(0)=x0.
The quadratic cost function associated with the player i is given by
(2)Ji(u1,…,uN)i=1,…,N=12xT(T)QiTx(T)+12∫0T[xT(t)Qix(t)+∑j=1NujT(t)Rijuj(t)]dt,
where T denotes the fixed prescribed duration of the game, x0 is the initial state, x∈ℝn, and ui∈ℝmi(i=1,…,N). The weighting matrices QiT, Qi, and Rij(i≠j) are assumed to be nonnegative; Rii is assumed to be positive definite.Definition 1.
The noncooperativeN-players LQ differential game is that the purpose of each player is to minimize the cost function Ji in  (2) under the controlled differential equation (1). It means that
(3)Ji(u1,…,uN)i=1,…,N→minui∈ℝmi.If each player employs strategy (3), then a conflict situation will be provoked. Thus a definition of an equilibrium behavior which in some sense is suitable for all players is demanded.Definition 2.
The system is described by (1), (2), and (3); a strategy u*(u1*,…,uN*) is said to be Nash-equilibrium strategy if for any admissible ui∈ℝmi(i=1,…,N):
(4)Ji(u1*,…,uN*)≤Ji(u1*,…,ui-1*,ui,ui+1*,…,uN*).This definition means that if all players are applying Nash-equilibrium strategiesu*(u1*,…,uN*), then there is no sense to change them unilaterally; that is, a player may only increase its losses.


### 2.2. Coupled Matrix Riccati Differential Equation

From the above description, the standard solution for a noncooperative and nonzero-sum LQ differential game problem can be obtained by calculus of variations. The theorem can be extracted from [11] as follows.Theorem 3.
For theN-players differential game, there exists a unique solution set Pi(t)(i=1,…,N) to the coupled matrix Riccati differential equations:
(5)-P˙i(t)=Qi+ATPi(t)+Pi(t)A-∑j=1NPi(t)SjPj(t),Pi(T)=QiT,
where Sj=BjRjj-1BjT. Then, the LQ differential game admits a unique Nash-equilibrium solution given by
(6)ui*(t)=-Rii-1BiTPi(t)x(t).
The proof can be found at [11].From the above theorem, the key for obtaining the Nash-equilibrium strategy of LQ differential game is to find the solution of the coupled matrix Riccati differential equation (5). However, it is not an easy task to find the solution of (5). It is comparatively rare that an exact analysis can be made. Numerical methods must be employed to solve this coupled matrix Riccati differential equation in most cases.The traditional or the most obvious method is to employ a numerical integration method. Then, the direct integration from the end timeT to the initial time 0 can be implemented by using an ordinary differential solver, such as ODE45 function in MATLAB software. However, this direct integration method has some difficulty in practical implementation. Because the coupled Riccati differential equations are high nonlinear differential equations, another method is to transfer the coupled Riccati differential equations into the linear ODE. Then it is solved by an ordinary differential solver.If we take all thePi(t)(i=1,…,N) into a single whole matrix, the coefficient matrices of (5) can be rewritten and composed as follows:
(7)W(t)=[P1(t)⋮PN(t)]∈ℝNn×n,(8)D=[-S1⋮-SN]T∈ℝn×Nn,(9)B=[-Q1⋮-QN]∈ℝNn×N,C=[-AT0⋯00-AT⋯0⋮⋮⋱⋮00⋯-AT]∈ℝNn×Nn.
Then the coupled matrix Riccati differential equation (5) can be rewritten as
(10)-W˙(t)=-B-CW(t)+W(t)A+W(t)DW(t).At last, the coupled matrix Riccati differential equations are transferred into a single whole asymmetric Riccati differential equation.Theorem 4.
The solution of (10) can be obtained by solving a linear ODE. Let y(t)∈ℝ(N+1)n×n be the solution of the following linear ODE with the terminal boundary condition:
(11)y˙(t)=[ADBC]y(t)=Hy(t),y(T)=[U(T)V1(T)⋮VN(T)]=[IQ1T⋮QNT].
By solving linear ODE (11), the solution of nonlinear asymmetric Riccati differential equation can be obtained by
(12)W(t)=[P1(t)⋮PN(t)]=[V1(t)⋮VN(t)]U-1(t).
The proof can be found at [12].From the above process, we can see that the solving of linear ODE is easier than the solving of nonlinear ODE. However, if the matricesU(t) or Vi(t) are ill-conditioned, then the solution of matrix Riccati differential equation (5) or asymmetric Riccati differential equation (10) cannot be obtained easily.


## 2.1. Problem Description

For the general case withN-players, the controlled equation can be defined for the LQ differential game problem as follows:
(1)x˙(t)=Ax(t)+∑i=1NBiui(t),x(0)=x0.
The quadratic cost function associated with the player i is given by
(2)Ji(u1,…,uN)i=1,…,N=12xT(T)QiTx(T)+12∫0T[xT(t)Qix(t)+∑j=1NujT(t)Rijuj(t)]dt,
where T denotes the fixed prescribed duration of the game, x0 is the initial state, x∈ℝn, and ui∈ℝmi(i=1,…,N). The weighting matrices QiT, Qi, and Rij(i≠j) are assumed to be nonnegative; Rii is assumed to be positive definite.Definition 1.
The noncooperativeN-players LQ differential game is that the purpose of each player is to minimize the cost function Ji in  (2) under the controlled differential equation (1). It means that
(3)Ji(u1,…,uN)i=1,…,N→minui∈ℝmi.If each player employs strategy (3), then a conflict situation will be provoked. Thus a definition of an equilibrium behavior which in some sense is suitable for all players is demanded.Definition 2.
The system is described by (1), (2), and (3); a strategy u*(u1*,…,uN*) is said to be Nash-equilibrium strategy if for any admissible ui∈ℝmi(i=1,…,N):
(4)Ji(u1*,…,uN*)≤Ji(u1*,…,ui-1*,ui,ui+1*,…,uN*).This definition means that if all players are applying Nash-equilibrium strategiesu*(u1*,…,uN*), then there is no sense to change them unilaterally; that is, a player may only increase its losses.


## 2.2. Coupled Matrix Riccati Differential Equation

From the above description, the standard solution for a noncooperative and nonzero-sum LQ differential game problem can be obtained by calculus of variations. The theorem can be extracted from [11] as follows.Theorem 3.
For theN-players differential game, there exists a unique solution set Pi(t)(i=1,…,N) to the coupled matrix Riccati differential equations:
(5)-P˙i(t)=Qi+ATPi(t)+Pi(t)A-∑j=1NPi(t)SjPj(t),Pi(T)=QiT,
where Sj=BjRjj-1BjT. Then, the LQ differential game admits a unique Nash-equilibrium solution given by
(6)ui*(t)=-Rii-1BiTPi(t)x(t).
The proof can be found at [11].From the above theorem, the key for obtaining the Nash-equilibrium strategy of LQ differential game is to find the solution of the coupled matrix Riccati differential equation (5). However, it is not an easy task to find the solution of (5). It is comparatively rare that an exact analysis can be made. Numerical methods must be employed to solve this coupled matrix Riccati differential equation in most cases.The traditional or the most obvious method is to employ a numerical integration method. Then, the direct integration from the end timeT to the initial time 0 can be implemented by using an ordinary differential solver, such as ODE45 function in MATLAB software. However, this direct integration method has some difficulty in practical implementation. Because the coupled Riccati differential equations are high nonlinear differential equations, another method is to transfer the coupled Riccati differential equations into the linear ODE. Then it is solved by an ordinary differential solver.If we take all thePi(t)(i=1,…,N) into a single whole matrix, the coefficient matrices of (5) can be rewritten and composed as follows:
(7)W(t)=[P1(t)⋮PN(t)]∈ℝNn×n,(8)D=[-S1⋮-SN]T∈ℝn×Nn,(9)B=[-Q1⋮-QN]∈ℝNn×N,C=[-AT0⋯00-AT⋯0⋮⋮⋱⋮00⋯-AT]∈ℝNn×Nn.
Then the coupled matrix Riccati differential equation (5) can be rewritten as
(10)-W˙(t)=-B-CW(t)+W(t)A+W(t)DW(t).At last, the coupled matrix Riccati differential equations are transferred into a single whole asymmetric Riccati differential equation.Theorem 4.
The solution of (10) can be obtained by solving a linear ODE. Let y(t)∈ℝ(N+1)n×n be the solution of the following linear ODE with the terminal boundary condition:
(11)y˙(t)=[ADBC]y(t)=Hy(t),y(T)=[U(T)V1(T)⋮VN(T)]=[IQ1T⋮QNT].
By solving linear ODE (11), the solution of nonlinear asymmetric Riccati differential equation can be obtained by
(12)W(t)=[P1(t)⋮PN(t)]=[V1(t)⋮VN(t)]U-1(t).
The proof can be found at [12].From the above process, we can see that the solving of linear ODE is easier than the solving of nonlinear ODE. However, if the matricesU(t) or Vi(t) are ill-conditioned, then the solution of matrix Riccati differential equation (5) or asymmetric Riccati differential equation (10) cannot be obtained easily.


## 3. Improved Precise Integration Method

In this section, we will firstly give the fundamental PIM for solving the asymmetric Riccati differential equation. In order to increase the efficiency of the fundamental PIM, then an improved precise integration method which is proposed by [10] for solving the symmetric Riccati differential equation will be expanded for solving the asymmetric Riccati differential equation.


### 3.1. PIM for Solving Asymmetric Riccati Differential Equation

Zhong and Zhu [8] firstly proposed the PIM for solving the symmetric Riccati differential equation (SRDE) by using the relationship between the solutions of SRDE and the fundamental solutions of the linear Hamiltonian system. Then this idea is developed for solving the asymmetric Riccati differential equation (ARDE) [9]. In this section, a brief introduction of PIM for solving ARDE is given.The PIM connects the coefficient matrices of the ARDE with the linear ODE as follows:(13)q˙=Aq+Dp,p˙=Bq+Cp.To find the solutions for the general ARDE, the relationship between the state vectorsqa,pa at time t=ta and the state vectors qb,pb at time t=tb is introduced. If the time interval [ta,tb](0≤ta<tb≤T) is considered as a subinterval within the whole integration interval [0,T], this relationship can be formulated as
(14)qb=Fqa-Gpb,pa=Qqa+Epb.It can be proved that the interval matricesF,G,Q, and E satisfy the following differential equations and boundary conditions [9]:
(15)G˙=-D+AG-GC+GBG=-FDE,F˙=(GB+A)F=F(A+DQ),E˙=E(BG-C)=(QD-C)E,Q˙=-EBF=-B+QA-CQ+QDQ,(16)G=0,Q=0,E=Ip,F=Iq,
where Iq and Ip are q and p dimension identity matrix, respectively.Comparing (15) with (10), it shows that both the matrix Q and the matrix W satisfy the same differential equation. Therefore, if the matrix Q defined by (14) can be obtained, then the solutions of the ARDE (10) can be given. In order to solve the ARDE, the continuous time should be discreted into a series of time intervals; that is, t0=0,…,tk=kη,…,tN=Mη=T·(η=T/M). According to the PIM, the subinterval η is further divided into 2N extremely small intervals with the equal length τ=η/2N. The matrices F-(τ),G(τ),Q(τ), and E-(τ) can be expressed as the function of the time interval τ by using a fourth-order Taylor series, and then the matrices F-(η),G(η),Q(η), and E-(η) can be combined N times with the following formulae:
(17)Gc=G+(Iq+F-)(Iq+GQ)-1G(Ip+E-),Qc=Q+(Ip+E-)Q(Iq+GQ)-1(Iq+F-),F-c=2F-+F-F--(Iq+F-)[GQ(Iq+GQ)-1+(Iq+GQ)-1GQ](Iq+F-)2,E-c=2E-+E-E--(Ip+E-)[QG(Ip+QG)-1+(Ip+QG)-1QG](Ip+E-)2.Then the matricesF(η),E(η) can be obtained by using the following formulae:
(18)F(η)=Iq+F-(η),E(η)=Ip+E-(η).Because the matricesF(η),G(η),Q(η), and E(η) for the subinterval η have been computed, the matrices F(kη),G(kη),Q(kη), and E(kη) at time kη can be obtained by implementing the following formulae:
(19)G(k(η+1))=G(kη)+F(kη)×(Iq+G(η)Q(kη))-1G(η)E(kη),Q(k(η+1))=Q(η)+E(η)×(Ip+Q(kη)G(η))-1Q(kη)F(η),F(k(η+1))=F(kη)(Iq+G(η)Q(kη))-1F(η),E(k(η+1))=E(η)(Ip+Q(kη)G(η))-1E(kη).According to the theorem of the computational structural mechanics [13], when the length of subinterval is zero, the subinterval matrices G, Q will become zero matrices. When the length of subinterval is going to infinity, the subinterval matrices G, Q will become constant matrices. For the linear system, the inversion of the matrices Iq+GQ and Ip+QG can always be guaranteed [13]. Furthermore, because the matrix Q and the matrix W satisfy the same differential equation, the solution for the ARDE (10) can be obtained according to (19):
(20)W(kη)=Q(T-kη)+E(T-kη)×(Ip+W(T)G(T-kη))-1W(T)F(T-kη).Thus, the PIM for solving the ARDE has been given by implementing the formulae (20) at last. From the above procedure of the PIM, the key of the PIM is the computation of the interval matrices F-(τ),G(τ),Q(τ), and E-(τ). The original PIM employs the Taylor series method with the fixed number of intervals N. With the increasing N, the precision of the PIM will increase. However, the computational work will increase much, especially for a large degree of freedom problem. Therefore, a main issue is how to compute the interval matrices F-(τ),G(τ),Q(τ), and E-(τ) with a higher precision and a lower computational work.


### 3.2. Improved Parameter Computation for PIM

Since the accuracy and efficiency of the PIM are highly dependent on the choices of parameterN, [10] has proposed an improved precise integration method (IPIM) for solving the symmetric Riccati differential equation. Therefore, in this paper, we expand the idea of the IPIM and propose the adaptive PIM for solving the asymmetric Riccati differential equation.The relationship between the adjacent state vectorsqa,pa and the state vectors qb,pb in the time interval [ta,tb] has been given in (14) by employing the interval matrices F,G,Q,  and E. Furthermore, this relationship can also be expressed with the solution of linear ODE (11) or (13) as follows:
(21)[qbpb]=[Φ11Φ12Φ21Φ22][qapa]=Φ[qapa],
where Φ denotes the exponential of the matrix H. Comparing (14) with (21), we can get the following relationship:
(22)F=Φ11-Φ12Φ22-1Φ21,G=-Φ12Φ22-1,Q=-Φ22-1Φ21,E=Φ22-1.In fact, on the one hand, the length of interval[ta,tb] cannot be made arbitrarily large. This is because the columns of Φ22 become more linearly dependent on the length of intervals [ta,tb] [14]. Consequently, if the choice of length of interval is not infinite and appropriate,  the inversion of the matrix Φ22 always exists. On the other hand, the computing of the exponential of the matrix H and then transforming it into interval matrices F,G,Q, and  E using (22) are feasible and stable for a small interval.The implementation of the PIM for computing the matrix exponentialΦ in one time step η can be given as follows:
(23)Φ=exp(Hη)=(exp(Hτ))2N.
Since the time step η is not very large and τ is extremely small, then the Taylor series approximation method can be used to evaluate Φ; that is,
(24)Φ=exp(Hτ)≈I+Ta,
where
(25)Ta=Hτ+H2τ22!+H3τ33!+⋯+Hqηqq!.
Substituting (24) into (23) gives
(26)Φ=(I+Ta)2N=(I+Ta)2(N-1)×(I+Ta)2(N-1).
Meanwhile, the addition theorem is applied only to the incremental part Ta, and it gives
(27)Ta⟵2Ta+Ta×Ta.
Finally, by performing (27) N times, the matrix exponential Φ is obtained by
(28)Φ=I+Ta.In order to choose the appropriate parametersN and q, a theorem has been given by [15] for the scaling and squaring method. Actually, the PIM considers the roundoff errors in the computational procedure. However, the scaling and squaring method does not. Therefore, the PIM is more stable. In this paper, we employ Theorem 5 to choose the parameters N and q adaptively.Theorem 5.
If∥X∥/2N≤1/2, then [Tq(∥X∥/2N)]2N=exp(X+E), where
(29)∥E∥∥X∥≤8(∥X∥2N)q1q+1≤(12)q-31q+1
in which ∥X∥ denotes the norm of the matrix X and Tq denotes the Taylor series of order q.According to the analysis given above, the detailed implementation of the adaptive PIM for the ARDE is listed as follows.Algorithm 6.
(
1
) Form the coefficient matrices A, B, C and D of the ADRE by (7)~(10), then the assembled matrix H can be obtained by (11);
(
2
) Give the duration of the game time T, then the time step is given by η=T/M;
(
3
) Let H-=Hη/2N and then determine the two parameters N and q according to the adaptive error Theorem 5 for a given error tolerance;
(
4
) Compute the q-order of the Taylor series approximation formulation Ta by (25);
(
5
) Perform the following stepsForiter=1:N
Ta=Ta(Ta+2I);End
(
6
) Get the matrix exponential Φ(η) by (28) and then transform it into the interval matrices F(η),G(η),Q(η),E(η) by (22);
(
7
) Compute interval matrices F(kη),G(kη),Q(kη),E(kη) at time kη by (19), and then the solution of the ARDE (10) can be obtained from the recursive formula (20);
(
8
) The Nash equilibrium solution of the LQ differential game can be given by (6).


## 3.1. PIM for Solving Asymmetric Riccati Differential Equation

Zhong and Zhu [8] firstly proposed the PIM for solving the symmetric Riccati differential equation (SRDE) by using the relationship between the solutions of SRDE and the fundamental solutions of the linear Hamiltonian system. Then this idea is developed for solving the asymmetric Riccati differential equation (ARDE) [9]. In this section, a brief introduction of PIM for solving ARDE is given.The PIM connects the coefficient matrices of the ARDE with the linear ODE as follows:(13)q˙=Aq+Dp,p˙=Bq+Cp.To find the solutions for the general ARDE, the relationship between the state vectorsqa,pa at time t=ta and the state vectors qb,pb at time t=tb is introduced. If the time interval [ta,tb](0≤ta<tb≤T) is considered as a subinterval within the whole integration interval [0,T], this relationship can be formulated as
(14)qb=Fqa-Gpb,pa=Qqa+Epb.It can be proved that the interval matricesF,G,Q, and E satisfy the following differential equations and boundary conditions [9]:
(15)G˙=-D+AG-GC+GBG=-FDE,F˙=(GB+A)F=F(A+DQ),E˙=E(BG-C)=(QD-C)E,Q˙=-EBF=-B+QA-CQ+QDQ,(16)G=0,Q=0,E=Ip,F=Iq,
where Iq and Ip are q and p dimension identity matrix, respectively.Comparing (15) with (10), it shows that both the matrix Q and the matrix W satisfy the same differential equation. Therefore, if the matrix Q defined by (14) can be obtained, then the solutions of the ARDE (10) can be given. In order to solve the ARDE, the continuous time should be discreted into a series of time intervals; that is, t0=0,…,tk=kη,…,tN=Mη=T·(η=T/M). According to the PIM, the subinterval η is further divided into 2N extremely small intervals with the equal length τ=η/2N. The matrices F-(τ),G(τ),Q(τ), and E-(τ) can be expressed as the function of the time interval τ by using a fourth-order Taylor series, and then the matrices F-(η),G(η),Q(η), and E-(η) can be combined N times with the following formulae:
(17)Gc=G+(Iq+F-)(Iq+GQ)-1G(Ip+E-),Qc=Q+(Ip+E-)Q(Iq+GQ)-1(Iq+F-),F-c=2F-+F-F--(Iq+F-)[GQ(Iq+GQ)-1+(Iq+GQ)-1GQ](Iq+F-)2,E-c=2E-+E-E--(Ip+E-)[QG(Ip+QG)-1+(Ip+QG)-1QG](Ip+E-)2.Then the matricesF(η),E(η) can be obtained by using the following formulae:
(18)F(η)=Iq+F-(η),E(η)=Ip+E-(η).Because the matricesF(η),G(η),Q(η), and E(η) for the subinterval η have been computed, the matrices F(kη),G(kη),Q(kη), and E(kη) at time kη can be obtained by implementing the following formulae:
(19)G(k(η+1))=G(kη)+F(kη)×(Iq+G(η)Q(kη))-1G(η)E(kη),Q(k(η+1))=Q(η)+E(η)×(Ip+Q(kη)G(η))-1Q(kη)F(η),F(k(η+1))=F(kη)(Iq+G(η)Q(kη))-1F(η),E(k(η+1))=E(η)(Ip+Q(kη)G(η))-1E(kη).According to the theorem of the computational structural mechanics [13], when the length of subinterval is zero, the subinterval matrices G, Q will become zero matrices. When the length of subinterval is going to infinity, the subinterval matrices G, Q will become constant matrices. For the linear system, the inversion of the matrices Iq+GQ and Ip+QG can always be guaranteed [13]. Furthermore, because the matrix Q and the matrix W satisfy the same differential equation, the solution for the ARDE (10) can be obtained according to (19):
(20)W(kη)=Q(T-kη)+E(T-kη)×(Ip+W(T)G(T-kη))-1W(T)F(T-kη).Thus, the PIM for solving the ARDE has been given by implementing the formulae (20) at last. From the above procedure of the PIM, the key of the PIM is the computation of the interval matrices F-(τ),G(τ),Q(τ), and E-(τ). The original PIM employs the Taylor series method with the fixed number of intervals N. With the increasing N, the precision of the PIM will increase. However, the computational work will increase much, especially for a large degree of freedom problem. Therefore, a main issue is how to compute the interval matrices F-(τ),G(τ),Q(τ), and E-(τ) with a higher precision and a lower computational work.


## 3.2. Improved Parameter Computation for PIM

Since the accuracy and efficiency of the PIM are highly dependent on the choices of parameterN, [10] has proposed an improved precise integration method (IPIM) for solving the symmetric Riccati differential equation. Therefore, in this paper, we expand the idea of the IPIM and propose the adaptive PIM for solving the asymmetric Riccati differential equation.The relationship between the adjacent state vectorsqa,pa and the state vectors qb,pb in the time interval [ta,tb] has been given in (14) by employing the interval matrices F,G,Q,  and E. Furthermore, this relationship can also be expressed with the solution of linear ODE (11) or (13) as follows:
(21)[qbpb]=[Φ11Φ12Φ21Φ22][qapa]=Φ[qapa],
where Φ denotes the exponential of the matrix H. Comparing (14) with (21), we can get the following relationship:
(22)F=Φ11-Φ12Φ22-1Φ21,G=-Φ12Φ22-1,Q=-Φ22-1Φ21,E=Φ22-1.In fact, on the one hand, the length of interval[ta,tb] cannot be made arbitrarily large. This is because the columns of Φ22 become more linearly dependent on the length of intervals [ta,tb] [14]. Consequently, if the choice of length of interval is not infinite and appropriate,  the inversion of the matrix Φ22 always exists. On the other hand, the computing of the exponential of the matrix H and then transforming it into interval matrices F,G,Q, and  E using (22) are feasible and stable for a small interval.The implementation of the PIM for computing the matrix exponentialΦ in one time step η can be given as follows:
(23)Φ=exp(Hη)=(exp(Hτ))2N.
Since the time step η is not very large and τ is extremely small, then the Taylor series approximation method can be used to evaluate Φ; that is,
(24)Φ=exp(Hτ)≈I+Ta,
where
(25)Ta=Hτ+H2τ22!+H3τ33!+⋯+Hqηqq!.
Substituting (24) into (23) gives
(26)Φ=(I+Ta)2N=(I+Ta)2(N-1)×(I+Ta)2(N-1).
Meanwhile, the addition theorem is applied only to the incremental part Ta, and it gives
(27)Ta⟵2Ta+Ta×Ta.
Finally, by performing (27) N times, the matrix exponential Φ is obtained by
(28)Φ=I+Ta.In order to choose the appropriate parametersN and q, a theorem has been given by [15] for the scaling and squaring method. Actually, the PIM considers the roundoff errors in the computational procedure. However, the scaling and squaring method does not. Therefore, the PIM is more stable. In this paper, we employ Theorem 5 to choose the parameters N and q adaptively.Theorem 5.
If∥X∥/2N≤1/2, then [Tq(∥X∥/2N)]2N=exp(X+E), where
(29)∥E∥∥X∥≤8(∥X∥2N)q1q+1≤(12)q-31q+1
in which ∥X∥ denotes the norm of the matrix X and Tq denotes the Taylor series of order q.According to the analysis given above, the detailed implementation of the adaptive PIM for the ARDE is listed as follows.Algorithm 6.
(
1
) Form the coefficient matrices A, B, C and D of the ADRE by (7)~(10), then the assembled matrix H can be obtained by (11);
(
2
) Give the duration of the game time T, then the time step is given by η=T/M;
(
3
) Let H-=Hη/2N and then determine the two parameters N and q according to the adaptive error Theorem 5 for a given error tolerance;
(
4
) Compute the q-order of the Taylor series approximation formulation Ta by (25);
(
5
) Perform the following stepsForiter=1:N
Ta=Ta(Ta+2I);End
(
6
) Get the matrix exponential Φ(η) by (28) and then transform it into the interval matrices F(η),G(η),Q(η),E(η) by (22);
(
7
) Compute interval matrices F(kη),G(kη),Q(kη),E(kη) at time kη by (19), and then the solution of the ARDE (10) can be obtained from the recursive formula (20);
(
8
) The Nash equilibrium solution of the LQ differential game can be given by (6).


## 4. Numerical Examples

In this section, two examples are given to show the validity of the PIM for solving noncooperative LQ differential game problem. Meanwhile, two other numerical methods are also employed to be compared with the PIM. The first numerical method is to directly integrate the coupled matrix Riccati differential equation (5) by ODE45 solver in MATLAB software. We call this first method “direct integration method”. The second numerical method is to integrate the linear ODE  (11) by ODE45 solver in MATLAB software. We call this second method “linear transformation method”. Because an adaptive step Runge-Kutta is implemented in ODE45 solver, the relative error and the absolute error are both set as 1×10-12 for the “direct integration method” and the “linear transformation method”.Example 7.
Considering a simple pursuit-evasion problem in space [7], the symbol u1(t) denotes the control acceleration of the purser, and the symbol u2(t) denotes the control acceleration of the evader. The symbol r(t) denotes the relative position of the purser with respect to the evader. The differential equation can be expressed by
(30)r˙(t)=v(t),v˙(t)=u1(t)-u2(t).
If we denote the state vector x(t)=[r(t),v(t)]T, then (30) can be rewritten as the following:
(31)x˙(t)=[0100]x(t)+[01]u1(t)+[0-1]u2(t).
The cost functions are selected by
(32)J1=12xT(1)[1000]x(1)+12∫011cu1T(t)u1(t)dt,J2=12xT(1)[-1000]x(1)+12∫01cu2T(t)u2(t)dt.
The exact solutions of the coupled matrix Riccati differential equation are
(33)P1(t)=1ω(t)[11-t1-t(1-t)2],ω(t)=1+13(c-1c)(1-t)3,P2(t)=-P1(t).From the exact solutions  (33), it is shown that the parameter c has important influence on the coupled matrix Riccati differential equation. Hence, two cases are considered, that is, Case 1“c=2” and Case 2  “c=0.2”. In both cases, the different numerical results are denoted by different symbols. The symbol “⋄” denotes the solution from the PIM, the symbol “×” denotes the solution from the “direct integration method”, and the symbol “°” denotes the solution from the “linear transformation method”.Case 1 (c=2).
The length of step for the PIM is given asη=0.05. Figure 1 shows the numerical solutions of the coupled matrix Riccati differential equation obtained from the PIM, the “direct integration method” and the “linear transformation method” for Case 1. It can be easily seen that the symbols “⋄”, “×”, and “°” overlap each other, and they are all nearly on the line of the exact solutions. Therefore, the above three methods all have enough precisions for Case 1.The numerical solutions of the coupled matrix Riccati differential equation forc=2.
(a)
(b)
(c)Case 2 (c=0.2).
The length of step for the PIM is also given asη=0.05. Figure 2 shows the numerical solutions of the coupled matrix Riccati differential equation obtained from the above three methods for Case 2. The three methods overlap each other after about t=0.2. Before t=0.2, the PIM and the “linear transformation method” are both nearly on the line of the exact solutions. However, the results from the “direct integration method” disappear in Figure 2. Therefore, the “direct integration method” fails for Case 2.The numerical solutions of the coupled matrix Riccati differential equation forc=0.2.
(a)
(b)
(c)From the exact solutions (33), it can be seen that the function ω(t) will change much for different parameter c. We plot the curves of the function ω(t) for Cases 1 and 2 in Figure 3. The line corresponding to Case 1 has no zero point. However, the line corresponding to Case 2 has a zero point. Consequently, we find the reason for the failure of the “direct integration method” in Case 2. Because the “direct integration method” which integrates the coupled matrix Riccati different equation from the end time T to the initial time 0 cannot go through the singularity point, the integration procedure stops at this singularity point.Figure 3
The values of the functionω(t) for different values of parameter c.Example 8.
Considering another LQ differential game problem [5], the controlled differential equation can be expressed by
(34)x˙(t)=[-0.03660.02710.0188-0.45550.0482-1.0100.0024-4.02080.10020.2855-0.70701.32290010]x(t)+[0.44223.0447-5.520]u1(t)+[0.1761-7.59224.990]u2(t).
The weighting matrices in the cost function are given by
(35)Q1=diag(3.5,2,4,5),Q2=diag(1.5,6,3,1),Q1T=I4,Q2T=I4,R11=1,R22=2.
The fixed prescribed duration of the game is T=10.This example has no analytical solutions, and we also employ three different numerical methods for solving this example. In Figure4, the solid lines denote the solution from the PIM, and the time step length is η=0.2; the symbol “×” denotes the solution from the “direct integration method”, and the symbol “°” denotes the solution from the “linear transformation method”.The numerical solutions of the coupled matrix Riccati differential equation.
(a)
(b)
(c)
(d)From Figure4, we can see that the numerical results from the “direct integration method” are nearly the same with the numerical results from the PIM in the all range of time T. However, the numerical results from the “linear transformation method” have much difference with the above two methods. Particularly, in the range of time interval [0,8], the “linear transformation method” gives the wrong results. Therefore, the “linear transformation method” fails for this example.The key point of the “linear transformation method” can be easily found in (12) in Theorem 4. Because the inversion matrix of U(t) exists in (12), if the matrices U(t) is ill-conditioned, then the solution of the coupled matrix Riccati differential equation (5) cannot be obtained. In order to find the failure reason of the “linear transformation method” for this example, we have given the determinants of the matrix functions U(t) and V(t) in Figure 5. We can see that the determinants of the matrix functions U(t) and V(t) from time t=10 to time t=8 increase much linearly. Therefore, the matrices functions U(t) and V(t) are both ill-conditioned matrices, and these causes led to the bad results.Figure 5
The determinants of the matrix functionsU(t) and V(t).


## 5. Conclusions

In order to obtain the Nash-equilibrium strategy of theN-players noncooperative LQ differential game problem, the key is to solve the coupled matrix Riccati differential equation. The precise integration method is expanded and applied for solving the coupled matrix Riccati differential equation in this paper. Furthermore, a comparison of the algorithm performance between the precise integration method, the “direct integration method”, and the “linear transformation method” is also given. At last, the numerical examples show that the “direct integration method” cannot go through the singularity point, and the integration procedure stops at this singularity point. The “linear transformation method” fails in the context of the ill-conditioned or singular matrix. However, the expanded precise integration method can always work and overcome the above difficulties. Therefore, the expanded precise integration method gives the more stable and accurate numerical results.



---

*Source: 713725-2013-06-10.xml*