Precise Integration Method for Solving Noncooperative LQ Differential Game

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”.

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 𝑁-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 with 𝑁-players, the controlled equation can be defined for the LQ differential game problem as follows: 𝐱̇(𝑡)=𝐀𝐱(𝑡)+^𝑁_𝑖=1𝐁_𝑖𝐮_𝑖(𝑡),𝐱(0)=𝐱₀. (1) The quadratic cost function associated with the player 𝑖 is given by 𝐽_𝑖𝐮₁,…,𝐮_{𝑁𝑖=1,…,𝑁}=12𝐱^𝑇(𝑇)𝐐_𝑖𝑇𝐱(𝑇)+12^𝑇₀𝐱^𝑇(𝑡)𝐐_𝑖𝐱(𝑡)+^𝑁_𝑗=1𝐮^𝑇_𝑗(𝑡)𝐑_𝑖𝑗𝐮_𝑗(𝑡)d𝑡, (2) where 𝑇 denotes the fixed prescribed duration of the game, 𝐱₀ is the initial state, 𝐱∈ℝ^𝑛, and 𝐮_𝑖∈ℝ^𝑚𝑖(𝑖=1,…,𝑁). The weighting matrices 𝐐_𝑖𝑇, 𝐐_𝑖, and 𝐑_𝑖𝑗(𝑖≠𝑗) are assumed to be nonnegative; 𝐑_𝑖𝑖 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 𝐽_𝑖𝐮₁,…,𝐮_{𝑁𝑖=1,…,𝑁}→min𝐮_𝑖∈ℝ^𝑚𝑖. (3)

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 𝐮^∗(𝐮^∗₁,…,𝐮^∗_𝑁) is said to be Nash-equilibrium strategy if for any admissible 𝐮_𝑖∈ℝ^𝑚𝑖(𝑖=1,…,𝑁): 𝐽_𝑖𝐮^∗₁,…,𝐮^∗_𝑁≤𝐽_𝑖𝐮^∗₁,…,𝐮^∗_𝑖−1,𝐮_𝑖,𝐮^∗_𝑖+1,…,𝐮^∗_𝑁. (4)

This definition means that if all players are applying Nash-equilibrium strategies 𝐮^∗(𝐮^∗₁,…,𝐮^∗_𝑁), 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 the 𝑁-players differential game, there exists a unique solution set 𝐏_𝑖(𝑡)(𝑖=1,…,𝑁) to the coupled matrix Riccati differential equations: −𝐏̇_𝑖(𝑡)=𝐐_𝑖+𝐀^𝑇𝐏_𝑖(𝑡)+𝐏_𝑖(𝑡)𝐀−^𝑁_𝑗=1𝐏_𝑖(𝑡)𝐒_𝑗𝐏_𝑗(𝑡),𝐏_𝑖(𝑇)=𝐐_𝑖𝑇, (5) where 𝐒_𝑗=𝐁_𝑗𝐑⁻¹_𝑗𝑗𝐁^𝑇_𝑗. Then, the LQ differential game admits a unique Nash-equilibrium solution given by 𝐮^∗_𝑖(𝑡)=−𝐑⁻¹_𝑖𝑖𝐁^𝑇_𝑖𝐏_𝑖(𝑡)𝐱(𝑡). (6) 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 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 𝐏_𝑖(𝑡)(𝑖=1,…,𝑁) into a single whole matrix, the coefficient matrices of (5) can be rewritten and composed as follows: 𝐖(𝑡)=𝐏₁(𝑡)⋮𝐏_𝑁(𝑡)∈ℝ^{𝑁𝑛×𝑛},𝐃=−𝐒₁⋮−𝐒_𝑁^𝑇∈ℝ^{𝑛×𝑁𝑛},𝐁=−𝐐₁⋮−𝐐_𝑁∈ℝ^{𝑁𝑛×𝑁},𝐂=−𝐀^𝑇𝟎⋯𝟎𝟎−𝐀^𝑇⋯𝟎⋮⋮⋱⋮𝟎𝟎⋯−𝐀^𝑇∈ℝ^{𝑁𝑛×𝑁𝑛}. (7)(8)(9) Then the coupled matrix Riccati differential equation (5) can be rewritten as −𝐖̇(𝑡)=−𝐁−𝐂𝐖(𝑡)+𝐖(𝑡)𝐀+𝐖(𝑡)𝐃𝐖(𝑡). (10)

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 𝐲(𝑡)∈ℝ^{(𝑁+1)𝑛×𝑛} be the solution of the following linear ODE with the terminal boundary condition: 𝐲̇(𝑡)=𝐲(𝑡)(𝑡),𝐲(𝑇)=𝐔(𝑇)𝐕𝐀𝐃𝐁𝐂=𝐇𝐲₁(𝑇)⋮𝐕_𝑁(𝑇)=𝐈𝐐_1𝑇⋮𝐐_𝑁𝑇. (11) By solving linear ODE (11), the solution of nonlinear asymmetric Riccati differential equation can be obtained by 𝐖(𝑡)=𝐏₁(𝑡)⋮𝐏_𝑁(𝑡)=𝐕₁(𝑡)⋮𝐕_𝑁(𝑡)𝐔⁻¹(𝑡). (12) 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 𝐔(𝑡) or 𝐕_𝑖(𝑡) 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)

To find the solutions for the general ARDE, the relationship between the state vectors 𝐪_𝑎,𝐩_𝑎 at time 𝑡=𝑡_𝑎 and the state vectors 𝐪_𝑏,𝐩_𝑏 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 𝐪_𝑏=𝐅𝐪_𝑎−𝐆𝐩_𝑏,𝐩_𝑎=𝐐𝐪_𝑎+𝐄𝐩_𝑏. (14)

It can be proved that the interval matrices 𝐅,𝐆,𝐐, and 𝐄 satisfy the following differential equations and boundary conditions [9]: ̇̇()(),̇()=()̇𝐆=−𝐃+𝐀𝐆−𝐆𝐂+𝐆𝐁𝐆=−𝐅𝐃𝐄,𝐅=𝐆𝐁+𝐀𝐅=𝐅𝐀+𝐃𝐐𝐄=𝐄𝐁𝐆−𝐂𝐐𝐃−𝐂𝐄,𝐐=−𝐄𝐁𝐅=−𝐁+𝐐𝐀−𝐂𝐐+𝐐𝐃𝐐,𝐆=𝟎,𝐐=𝟎,𝐄=𝐈_𝑝,𝐅=𝐈_𝑞, (15)(16) where 𝐈_𝑞 and 𝐈_𝑝 are 𝑞 and 𝑝 dimension identity matrix, respectively.

Comparing (15) with (10), it shows that both the matrix 𝐐 and the matrix 𝐖 satisfy the same differential equation. Therefore, if the matrix 𝐐 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,…,𝑡_𝑘=𝑘𝜂,…,𝑡_𝑁=𝑀𝜂=𝑇⋅(𝜂=𝑇/𝑀). According to the PIM, the subinterval 𝜂 is further divided into 2^𝑁 extremely small intervals with the equal length 𝜏=𝜂/2^𝑁. The matrices 𝐅(𝜏),𝐆(𝜏),𝐐(𝜏), and 𝐄(𝜏) can be expressed as the function of the time interval 𝜏 by using a fourth-order Taylor series, and then the matrices 𝐅(𝜂),𝐆(𝜂),𝐐(𝜂), and 𝐄(𝜂) can be combined 𝑁 times with the following formulae: 𝐆_𝑐=𝐆+𝐈_𝑞+𝐅𝐈_𝑞+𝐆𝐐⁻¹𝐆𝐈_𝑝+𝐄,𝐐_𝑐=𝐐+𝐈_𝑝+𝐄𝐐𝐈_𝑞+𝐆𝐐⁻¹𝐈_𝑞+𝐅,𝐅_𝑐=2𝐅+𝐅𝐅−𝐈_𝑞+𝐅𝐆𝐐𝐈_𝑞+𝐆𝐐⁻¹+𝐈_𝑞+𝐆𝐐⁻¹𝐆𝐐𝐈_𝑞+𝐅2,𝐄_𝑐=2𝐄+𝐄𝐄−𝐈_𝑝+𝐄𝐐𝐆𝐈_𝑝+𝐐𝐆⁻¹+𝐈_𝑝+𝐐𝐆⁻¹𝐐𝐆𝐈_𝑝+𝐄2. (17)

Then the matrices 𝐅(𝜂),𝐄(𝜂) can be obtained by using the following formulae: 𝐅𝜂=𝐈_𝑞+𝐅𝜂,𝐄𝜂=𝐈_𝑝+𝐄𝜂. (18)

Because the matrices 𝐅(𝜂),𝐆(𝜂),𝐐(𝜂), and 𝐄(𝜂) for the subinterval 𝜂 have been computed, the matrices 𝐅(𝑘𝜂),𝐆(𝑘𝜂),𝐐(𝑘𝜂), and 𝐄(𝑘𝜂) at time 𝑘𝜂 can be obtained by implementing the following formulae: 𝐆𝑘×𝐈𝜂+1=𝐆𝑘𝜂+𝐅𝑘𝜂_𝑞+𝐆𝜂𝐐𝑘𝜂⁻¹𝐆𝜂𝐄,𝐐𝑘𝜂𝜂×𝐈𝑘𝜂𝜂+1=𝐐+𝐄_𝑝+𝐐𝑘𝜂𝐆𝜂⁻¹𝐐𝐅𝜂,𝐅𝑘𝐈𝑘𝜂𝜂+1=𝐅𝑘𝜂_𝑞+𝐆𝜂𝐐𝑘𝜂⁻¹𝐅𝜂,𝐄𝑘𝜂+1=𝐄𝜂𝐈_𝑝+𝐐𝑘𝜂𝐆𝜂⁻¹𝐄𝑘𝜂. (19)

According to the theorem of the computational structural mechanics [13], when the length of subinterval is zero, the subinterval matrices 𝐆, 𝐐 will become zero matrices. When the length of subinterval is going to infinity, the subinterval matrices 𝐆, 𝐐 will become constant matrices. For the linear system, the inversion of the matrices 𝐈_𝑞+𝐆𝐐 and 𝐈_𝑝+𝐐𝐆 can always be guaranteed [13]. Furthermore, because the matrix 𝐐 and the matrix 𝐖 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 𝐅(𝜏),𝐆(𝜏),𝐐(𝜏), and 𝐄(𝜏). 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 𝐅(𝜏),𝐆(𝜏),𝐐(𝜏), and 𝐄(𝜏) 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 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 𝐪_𝑎,𝐩_𝑎 and the state vectors 𝐪_𝑏,𝐩_𝑏 in the time interval [𝑡_𝑎,𝑡_𝑏] has been given in (14) by employing the interval matrices 𝐅,𝐆,𝐐,  and 𝐄. Furthermore, this relationship can also be expressed with the solution of linear ODE (11) or (13) as follows: 𝐪_𝑏𝐩_𝑏=Φ₁₁Φ₁₂Φ₂₁Φ₂₂𝐪_𝑎𝐩_𝑎𝐪=Φ^𝑎𝐩_𝑎, (21) where Φ denotes the exponential of the matrix 𝐇. Comparing (14) with (21), we can get the following relationship: 𝐅=Φ₁₁−Φ₁₂Φ⁻¹₂₂Φ₂₁,𝐆=−Φ₁₂Φ⁻¹₂₂,𝐐=−Φ⁻¹₂₂Φ₂₁,𝐄=Φ⁻¹₂₂. (22)

In fact, on the one hand, the length of interval [𝑡_𝑎,𝑡_𝑏] cannot be made arbitrarily large. This is because the columns of Φ₂₂ 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 Φ₂₂ always exists. On the other hand, the computing of the exponential of the matrix 𝐇 and then transforming it into interval matrices 𝐅,𝐆,𝐐, and  𝐄 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: Φ=exp𝐇𝜂=exp(𝐇𝜏)^2𝑁. (23) 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, Φ=exp(𝐇𝜏)≈𝐈+𝐓_𝑎, (24) where 𝐓_𝑎𝐇=𝐇𝜏+²𝜏²+𝐇2!³𝜏³𝐇3!+⋯+^𝑞𝜂^𝑞.𝑞! (25) Substituting (24) into (23) gives Φ=𝐈+𝐓_𝑎^2𝑁=𝐈+𝐓_𝑎^2(𝑁−1)×𝐈+𝐓_𝑎^2(𝑁−1). (26) Meanwhile, the addition theorem is applied only to the incremental part 𝐓_𝑎, and it gives 𝐓_𝑎⟵2𝐓_𝑎+𝐓_𝑎×𝐓_𝑎. (27) Finally, by performing (27) 𝑁 times, the matrix exponential Φ is obtained by Φ=𝐈+𝐓_𝑎. (28)

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.

Theorem 5. If ‖𝐗‖/2^𝑁≤1/2, then [𝐓_𝑞(‖𝐗‖/2^𝑁)]^2𝑁=exp(𝐗+𝐄), where ‖𝐄‖‖𝐗‖‖𝐗‖≤82^𝑁𝑞1≤1𝑞+12^𝑞−31𝑞+1 (29) in which ‖𝐗‖ denotes the norm of the matrix 𝐗 and 𝐓_𝑞 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 𝐀, 𝐁, 𝐂 and 𝐃 of the ADRE by (7)~(10), then the assembled matrix 𝐇 can be obtained by (11);
(2) Give the duration of the game time 𝑇, then the time step is given by 𝜂=𝑇/𝑀;
(3) Let 𝐇=𝐇𝜂/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 𝐓_𝑎 by (25);
(5) Perform the following steps For iter=1∶𝑁 𝐓_𝑎=𝐓_𝑎(𝐓_𝑎+2𝐈); End
(6) Get the matrix exponential Φ(𝜂) by (28) and then transform it into the interval matrices 𝐅(𝜂),𝐆(𝜂),𝐐(𝜂),𝐄(𝜂) by (22);
(7) Compute interval matrices 𝐅(𝑘𝜂),𝐆(𝑘𝜂),𝐐(𝑘𝜂),𝐄(𝑘𝜂) 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).

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⁻¹² for the “direct integration method” and the “linear transformation method”.

Example 7. Considering a simple pursuit-evasion problem in space [7], the symbol 𝑢₁(𝑡) denotes the control acceleration of the purser, and the symbol 𝑢₂(𝑡) 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 𝑟̇(𝑡)(𝑡),𝑣̇(𝑡)=𝑣=𝑢₁(𝑡)−𝑢₂(𝑡). (30) If we denote the state vector 𝐱(𝑡)=[𝑟(𝑡),𝑣(𝑡)]^𝑇, then (30) can be rewritten as the following: 𝐱̇(𝑡)=𝐱(𝑡)+01𝑢0100¹(𝑡)+0𝑢−1²(𝑡). (31) The cost functions are selected by 𝐽₁=12𝐱^𝑇(1)𝐱(1)[phone]𝑐𝑢^𝑇₁(𝑡)𝑢₁(𝑡)𝐽d𝑡,₂=12𝐱^𝑇(1)𝐱(1)+−101002¹₀𝑐𝑢^𝑇₂(𝑡)𝑢₂(𝑡)d𝑡. (32) The exact solutions of the coupled matrix Riccati differential equation are 𝐏₁(𝑡)=1𝜔(𝑡)()11−𝑡1−𝑡1−𝑡²,𝜔(𝑡)1=1+31𝑐−𝑐(1−𝑡)³,𝐏₂(𝑡)=−𝐏₁(𝑡). (33)

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.

(a)

(b)

(c)

Figure 1 

The numerical solutions of the coupled matrix Riccati differential equation for 𝑐=2.

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.

(a)

(b)

(c)

Figure 2 

The numerical solutions of the coupled matrix Riccati differential equation for 𝑐=0.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.

Figure 3 

The values of the function 𝜔(𝑡) for different values of parameter 𝑐.

Example 8. Considering another LQ differential game problem [5], the controlled differential equation can be expressed by 𝐱̇(𝑡)=𝐱(𝑡)+0𝑢−0.03660.02710.0188−0.45550.0482−1.0100.0024−4.02080.10020.2855−0.70701.322900100.44223.0447−5.52¹(𝑡)+0𝑢0.1761−7.59224.99²(𝑡). (34) The weighting matrices in the cost function are given by 𝐐₁=diag(3.5,2,4,5),𝐐₂(),𝐐=diag1.5,6,3,1_1𝑇=𝐈₄,𝐐_2𝑇=𝐈₄,𝐑₁₁=1,𝐑₂₂=2. (35) 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”.

(a)

(b)

(c)

(d)

Figure 4 

The numerical solutions of the coupled matrix Riccati differential equation.

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 𝐔(𝑡) exists in (12), if the matrices 𝐔(𝑡) 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 𝐔(𝑡) and 𝐕(𝑡) in Figure 5. We can see that the determinants of the matrix functions 𝐔(𝑡) and 𝐕(𝑡) from time 𝑡=10 to time 𝑡=8 increase much linearly. Therefore, the matrices functions 𝐔(𝑡) and 𝐕(𝑡) are both ill-conditioned matrices, and these causes led to the bad results.

Figure 5 

The determinants of the matrix functions 𝐔(𝑡) and 𝐕(𝑡).

5. 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.

Acknowledgments

The project is supported by the National Science Foundation of China (11102031, 11072044) and the Fundamental Research Funds for the Central Universities (DUT13LK25).