Precise Integration Method for Solving Noncooperative LQ Differential Game

The key of solving the noncooperative linear quadratic (LQ) differential game is to solve the coupledmatrix 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”.


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 -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][9][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.

Noncooperative LQ Differential Game
2.1.Problem Description.For the general case with -players, the controlled equation can be defined for the LQ differential game problem as follows: The quadratic cost function associated with the player  is given by where  denotes the fixed prescribed duration of the game, x 0 is the initial state, x ∈ R  , and u  ∈ R   ( = 1, . . ., ).The weighting matrices Q  , Q  , and R  ( ̸ = ) are assumed to be nonnegative; R  is assumed to be positive definite.Definition 1.The noncooperative -players LQ differential game is that the purpose of each player is to minimize the cost function   in (2) under the controlled differential equation (1).It means that 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) ), then there is no sense to change them unilaterally; that is, a player may only increase its losses.

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 the -players differential game, there exists a unique solution set P  () ( = 1, . . ., ) to the coupled matrix Riccati differential equations: where Then, the LQ differential game admits a unique Nash-equilibrium solution given by The proof can be found at [11].
From the above theorem, the key for obtaining the Nashequilibrium 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 time  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 the P  () ( = 1, . . ., ) into a single whole matrix, the coefficient matrices of ( 5) can be rewritten and composed as follows: Then the coupled matrix Riccati differential equation ( 5) can be rewritten as 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() ∈ R (+1)× be the solution of the following linear ODE with the terminal boundary condition: By solving linear ODE (11), the solution of nonlinear asymmetric Riccati differential equation can be obtained by 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 matrices U() or V  () are ill-conditioned, then the solution of matrix Riccati differential equation (5) or asymmetric Riccati differential equation ( 10) cannot be obtained easily.

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.

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: To find the solutions for the general ARDE, the relationship between the state vectors q  , p  at time  =   and the state vectors q  , p  at time  =   is introduced.If the time interval [  ,   ] (0 ≤   <   ≤ ) is considered as a subinterval within the whole integration interval [0, ], this relationship can be formulated as It can be proved that the interval matrices F, G, Q, and E satisfy the following differential equations and boundary conditions [9]: where I  and I  are  and  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,  0 = 0, . . .,   = , . . .,   =  =  ⋅ ( = /).According to the PIM, the subinterval  is further divided into 2  extremely small intervals with the equal length  = /2  .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  times with the following formulae: 2 , Then the matrices F(), E() can be obtained by using the following formulae: Because the matrices F(), G(), Q(), and E() for the subinterval  have been computed, the matrices F(), G(), Q(), and E() at time  can be obtained by implementing the following formulae: 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 I  + GQ and I  + 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) 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 .With the increasing , 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.

Improved Parameter Computation for PIM.
Since the accuracy and efficiency of the PIM are highly dependent on the choices of parameter , [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 vectors q  , p  and the state vectors q  , p  in the time interval [  ,   ] 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: where Φ denotes the exponential of the matrix H. Comparing ( 14) with (21), we can get the following relationship: In fact, on the one hand, the length of interval [  ,   ] cannot be made arbitrarily large.This is because the columns of Φ 22 become more linearly dependent on the length of intervals [  ,   ] [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: 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, where Substituting ( 24) into (23) gives × (I + T  ) Meanwhile, the addition theorem is applied only to the incremental part T  , and it gives Finally, by performing (27)  times, the matrix exponential Φ is obtained by In order to choose the appropriate parameters  and , 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  and  adaptively.
in which ‖X‖ denotes the norm of the matrix X and T  denotes the Taylor series of order .
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 , then the time step is given by  = /; (3) Let H = H/2  and then determine the two parameters  and  according to the adaptive error Theorem 5 for a given error tolerance; (4) Compute the -order of the Taylor series approximation formulation T  by (25); (5) Perform the following steps For iter = 1 :  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(), G(), Q(), E() at time  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).

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  1 () denotes the control acceleration of the purser, and the symbol  2 () denotes the control acceleration of the evader.The symbol () denotes the relative position of the purser with respect to the evader.The differential equation can be expressed by If we denote the state vector x() = [(), V()]  , then (30) can be rewritten as the following: The cost functions are selected by The exact solutions of the coupled matrix Riccati differential equation are From the exact solutions (33), it is shown that the parameter  has important influence on the coupled matrix Riccati differential equation.Hence, two cases are considered, that is, Case 1 " = 2" and Case 2 " = 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 ( = 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.
Case 2 ( = 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  = 0.2.Before  = 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.
From the exact solutions (33), it can be seen that the function () will change much for different parameter .We plot the curves of the function () 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  to the initial time 0 cannot go through the singularity point, the integration procedure stops at this singularity point.
The weighting matrices in the cost function are given by The fixed prescribed duration of the game is  = 10.
This example has no analytical solutions, and we also employ three different numerical methods for solving this example.In Figure 4, 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".
From Figure 4, 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 .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() exists in (12), if the matrices U() is illconditioned, 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() and V() in Figure 5.We can see that the determinants of the matrix functions U() and V() from time  = 10 to time  = 8 increase much linearly.Therefore, the matrices functions U() and V() are both ill-conditioned matrices, and these causes led to the bad results.

Conclusions
In order to obtain the Nash-equilibrium strategy of the 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.

Figure 1 :
Figure 1: The numerical solutions of the coupled matrix Riccati differential equation for  = 2.

Figure 2 :Figure 3 :
Figure 2: The numerical solutions of the coupled matrix Riccati differential equation for  = 0.2.

Figure 4 :Example 8 .
Figure 4: The numerical solutions of the coupled matrix Riccati differential equation.

Figure 5 :
Figure 5: The determinants of the matrix functions U() and V().