# An Algorithm for Solving a Class of Multiplayer Feedback-Nash Differential Games

Authors: Jorge Herrera de la Cruz; Benjamin Ivorra; Ángel M. Ramos
Journal: Mathematical Problems in Engineering
 (2019)

---

## Abstract

In this work, we introduce a novel numerical algorithm, called RaBVItG (Radial Basis Value Iteration Game) to approximate feedback-Nash equilibria for deterministic differential games. More precisely, RaBVItG is an algorithm based on value iteration schemes in a meshfree context. It is used to approximate optimal feedback Nash policies for multiplayer, trying to tackle the dimensionality that involves, in general, this type of problems. Moreover, RaBVItG also implements a game iteration structure that computes the game equilibrium at every value iteration step, in order to increase the accuracy of the solutions. Finally, with the purpose of validating our method, we apply this algorithm to a set of benchmark problems and compare the obtained results with the ones returned by another algorithm found in the literature. When comparing the numerical solutions, we observe that our algorithm is less computationally expensive and, in general, reports lower errors.


---

## Body

## 1. Introduction

From a general point of view, differential games (DG) is a topic involving game theory and controlled systems of differential equations. Here, we focus on conflict problems in whichN players control a system with an associated cost function per player. There exist two main types of control systems for these problems:closed-loop(or feedback), in which the controls (optimal solutions) are functions of the state variables of the system; andopen-loopsolutions, in which controls are only functions of time. Focusing on differential games, feedback controls are more robust than open-loop controls due to the possibility of the players to react to sudden changes in the values of the state variables of their opponents while they are playing the game. Based on the behavior of the N players, it is possible to classify the game and its corresponding equilibrium solutions ascooperativeandnoncooperative.The first class is used when the N players make coordinated decisions trying to optimize a joint criterium, while the second class is used when the N players compete each other trying to optimize their own criterium, taking into account the strategies from the rest N-1 players. We note that there also exist intermediate possibilities.In this work, we focus on Nash (i.e.,noncooperative)feedback solutions (i.e.,closed-loop). Our aim is twofold. Firstly, we develop an algorithm to solve numerically a deterministic multiplayer feedback-Nash differential game (FNDG). The algorithm is called RaBVIt (Radial Basis Value Iteration Game) and is based on radial basis, value iteration, and game iteration. Secondly, we apply the algorithm to some benchmark problems found in the literature of differential games (generally, linear or linear quadratic examples with explicit solutions) in order to compare the performance of RaBVItG to a previous algorithm (see [1], for more details).Focusing on recent works, we highlight the fact that, in a recent survey (see, [2]), the authors report that some existing numerical (computational) methods in the literature of DG are focused in some feasible subclass of games: two-person zero-sum games (see, [3]). These games do not strictly match our purpose (of dealing with N players nonzero games) but are interesting, from an analytical point of view, because they tackle the lack of differentiability of the value function by using the viscosity solution, as this scheme can be reduced to a one-dimensional game. Viscosity solutions in N dimensional games are, in general, complex to obtain. Additionally, other advances are in the line of linear quadratic models (see, [4]), but this framework omits to deal with interesting nonlinearities in the dynamic system or running cost functions. According to the authors, computational differential games is an area that needs to grow up.In this context, RaBVItG is an algorithm based on some recent reinforcement learning schemes (see, e.g., [5]) which are discrete-time methods closely related to dynamic programming. More precisely, we consider methods from reinforcement learning algorithms that are, mainly, used to simulate and approximate the behavior (in their usual terminology) of a set ofagents that takeactions in order to get a cumulativereward. In particular, we focus on value iteration (VI) which is a general technique used by reinforcement learning to iterate over the value function of the problem in order to obtain a fixed point solution. In our case, we have a coupled system of value functions (one per player) involved in a FNDG. We also use in RaBVItG the concept of game iteration (GI), meaning that, at eachvalue iteration step, we iterate again to find the corresponding game equilibrium (Nash, in this case) associated with the set of value functions. This algorithm is used to simulate the way in which the players make decisions (for instance, in a Nash context, fixing the opponent’s strategies to obtain the optimal strategy for a player). Finally, we also use function approximation techniques that allow us to simplify the model in order to approximate, using mesh-free techniques, the value function for each player (see, for instance, [6, 7]).This numerical method has been designed for solving a coupled system ofN-player Hamilton-Jacobi-Bellman equations (HJB). HJB equations provide the value function which gives optimal policies for a given dynamic system with an associated cost function (see, for example, [8]). They are the key for finding a feedback solution. In order to discretize our problem (in time and space), we use the techniques developed in [9], which introduces a semi-Lagrangian discretization framework for HJB equations and proves the convergence of the scheme based on the viscosity solution (see, e.g., [10]).In order to validate our algorithm, we apply it to solve a set of benchmark problems found in the literature (see, [1, 4]). We compare the obtained CPU time and error with the ones returned by another numerical algorithm found in [1].This paper is organized as follows. In Section2, we present the theoretical model by describing the relevant variables involved in the game, the coupled optimization problems, and the basic Nash equilibrium concepts. In Section 3, we explain the numerical implementation of our method, based on a semi-Lagrangian discretization, value iteration, and radial basis interpolation. In Section 4, we show the performance of the method by solving some benchmark problems.


## 2. Materials and Methods

This section deals with the explanation of the considered deterministic theoretical model. Firstly, we introduce the differential game of interest. Secondly, we define the considered feedback Nash equilibrium and the Hamilton-Jacobi-Bellman equations.


### 2.1. Deterministic Differential Game

We first define the class of differential games we are dealing with in this work.Let us consider a set ofN players. Each player, i∈1,2,…,N, has a payoff functional, Jiy,u, given by(1)Jiy,u=∫0∞e-ρitfixt,utdt,with y∈X, u∈U, and(2)x˙t=gxt,utt∈0,∞,and(3)x0=y.Defining the following, (1)
u:0,∞→U is the control function with U⊂Rm, so that U=U1×U2×⋯×UN. Let us define ui:0,∞→Ui, as the control associated with thei- th player, with Ui⊂Rmi a given subset of admissible controls. We also denote by Ui=ui:0,∞→Ui,suchthatuiismeasurable, U=U1×U2⋯×UN, and u(t)=u1t,…,uNt∈U.(2)
x:0,∞→X is a function such that x(t) denotes the state of the system at time t, where X⊂RP is the set of admissible states (p≠N). The evolution of those state variables is driven by (2), called thestate equation. We note that y∈X and g:D→RP, with D⊂RP×Rm, and assume that g is continuous and it exists Lg such that gx1,u-gx2,u≤Lgx1-x2, for all x1,x2∈X, t∈R+, and u∈U. Those conditions ensure that the system has a unique solution (this can be proved by using the Carathèodory theorem, see [9]).(3)
Ji:X×U→R is the payoff functional of player i. We assume fi:D→R, where fixt,u(t) represents the instantaneous payoff of player i for the choice uit. Note that the integral is affected by an exponential discounting parameter ρi>0 that actualizes the value of the payoff (see, for instance, [4]). The presence of the discount factor ensures that integral (1) is finite whenever fi is bounded.


### 2.2. Feedback Optimization Problem

In this Section, we recall thesynthesis proceduredetailed in, e.g., [11].Letvi:X→R be a continuously differentiable mapping calledvalue functionand defined by(4)viy=supu∈U⁡Jiy,u.For eachy∈X, we assume the existence of, at least, one optimal control uy∗∈U such that(5)viy=maxu∈U⁡Jiy,u=Jiy,uy∗=∫0∞e-ρitfixyt,uy∗,uy∗tdt,where xy(t,uy) is an admissible trajectory satisfying (2)-(3). We denote by xy∗(t)=xy(t,uy∗) the optimal trajectory.According to [11], the function hy,uyi0,∞→R given by(6)hy,uyiτ≔∫0τfi(xyt;uy,uyte-ρitdt+vixyτ;uye-ρiτis constant and nonincreasing for all t>0 if and only if uy∈U is an optimal control for the initial condition y. Thus, considering uy∗∈U, and derivating hy,uyi with respect to τ, we obtain(7)ρivixy∗τ-fixy∗τ,uy∗τ-∇vixy∗τTgxy∗τ,uy∗τ=0.Furthermore, forτ=0, we have that(8)ρiviy-fiy,uy∗0-∇viyTgy,uy∗0=0.We note that, for any other admissible constant control of the formu(t)=ξ∈U, for all t≥0, as hy,uyi(τ) is not increasing, we have that(9)ddτ∫0τfixyτ,ξe-ρitdt+vxyτ,ξe-ρtτ≤0. Again, for τ=0, we obtain(10)ρiviy-fiy,ξ-∇viyTgy,ξ≤0. So, under previous hypothesis, we conclude that for all y∈X we satisfy the following HJB equation:(11)ρiviy=maxu∈U⁡fiy,u+∇viyTgy,u,anduy∗ is such that(12)uy∗0∈arg⁡maxξ∈U⁡fiy,ξ+∇viyTgy,ξ.We remark that, in general, ifvi is only continuous, such equation needs to be interpreted in terms of viscosity solution of Problem (11). However, when it applies to differential games, there is, so far, no general theorem on the existence or uniqueness of solutions (see [12]).Now, we defineSi:X→Ui calledfeedback-mapper player i=1,…,N, such that(13)Siy∈arg⁡maxu∈U⁡fiy,u+∇ViyTgy,u,wherey∈Rp. We consider an optimal control defined by(14)uy∗t=Sixy∗t,for almost everyt>0.The abovementioned synthesis procedure consists in obtaining, using thefeedback-map,an optimal decision uy∗(t)=Si(xy∗(t)) related to the corresponding optimal trajectory, by solving(15)x˙t=gxt,S1xt,…,SNxtt∈0,∞,and(16)x0=y.So, provided an initial positiony, we say that uy∗,xy∗∈U×C[0,∞),X is an optimal paircontrol-trajectory for every initial condition, and it corresponds to the optimal feedback policy we are trying to estimate.


### 2.3. Feedback-Nash Equilibrium

Now, we adapt the previous expressions in order to get a feedback-Nash equilibrium.To do so, we define afeedback-Nash mapSi:X→Ui, per player i=1,…,N, such that, for all y∈Rp,(17)Siy∈arg⁡maxui∈Ui⁡fiy,ui;u-i+∇ViyTgy,ui;u-i,where the array ui;u-i denotes the pair of the control for player i and the controls associated with the rest of players (denoted by -i).Considering this feedback-Nash map, we apply the synthesis procedure described previously. To do so, we define a feedback N-tupleS≡S1,…,Si-1,Si,Si+1,…,SN and a feedback (N-1)-tuple S-i≡S1,…,Si-1,Si+1,…,SN. Then,S∗≡S∗,1,…,S∗,i-1,S∗,i,S∗,i+1,…,S∗,N is a feedback-Nash equilibrium (FNE) if(18)Jiy;S∗≥Jiy;ui,S∗,-iforeachy∈X,Si:X→Uiandi∈1,…,N,where Ji(y,Si) denotes Ji(y,Si(y)) and ui,S∗,-i is the N-vector of controls obtained by replacing the ith component in Si∗ by ui.Assumingy and S∗,-i are fixed, finding ui by maximizing Jiy;ui,S∗.-i (i.e., solving (18)) can be solved by using dynamic programming (see, for instance, [7, 13]).Finally, we define(19)viy;S∗,-i=maxui∈Ui⁡Jiy;ui,S∗,-i,which is the solution of the following HJB equation (see [3]):(20)ρiviy;S∗,-i=maxui∈Ui⁡fiy,uiS∗,-i+∇Viy;S∗,-iTgy,ui,S∗,-i,where the notationvi.;S∗,-i means that we are finding a Nash equilibrium for player i, by fixing the optimal strategies of the N-1 players. We note that, regarding current literature, apart from zero-sum games, to find relevant cases where (20) is well posed, we need to focus on games in one spatial dimension. For instance, in [14], an existence theorem of Nash equilibria in feedback form, valid for one-space dimensionnoncooperativegames, is given. However, as far as we know, there are no general existence theorems for feedback-Nash equilibrium (for spaces in dimension greater than 1).


## 2.1. Deterministic Differential Game

We first define the class of differential games we are dealing with in this work.Let us consider a set ofN players. Each player, i∈1,2,…,N, has a payoff functional, Jiy,u, given by(1)Jiy,u=∫0∞e-ρitfixt,utdt,with y∈X, u∈U, and(2)x˙t=gxt,utt∈0,∞,and(3)x0=y.Defining the following, (1)
u:0,∞→U is the control function with U⊂Rm, so that U=U1×U2×⋯×UN. Let us define ui:0,∞→Ui, as the control associated with thei- th player, with Ui⊂Rmi a given subset of admissible controls. We also denote by Ui=ui:0,∞→Ui,suchthatuiismeasurable, U=U1×U2⋯×UN, and u(t)=u1t,…,uNt∈U.(2)
x:0,∞→X is a function such that x(t) denotes the state of the system at time t, where X⊂RP is the set of admissible states (p≠N). The evolution of those state variables is driven by (2), called thestate equation. We note that y∈X and g:D→RP, with D⊂RP×Rm, and assume that g is continuous and it exists Lg such that gx1,u-gx2,u≤Lgx1-x2, for all x1,x2∈X, t∈R+, and u∈U. Those conditions ensure that the system has a unique solution (this can be proved by using the Carathèodory theorem, see [9]).(3)
Ji:X×U→R is the payoff functional of player i. We assume fi:D→R, where fixt,u(t) represents the instantaneous payoff of player i for the choice uit. Note that the integral is affected by an exponential discounting parameter ρi>0 that actualizes the value of the payoff (see, for instance, [4]). The presence of the discount factor ensures that integral (1) is finite whenever fi is bounded.


## 2.2. Feedback Optimization Problem

In this Section, we recall thesynthesis proceduredetailed in, e.g., [11].Letvi:X→R be a continuously differentiable mapping calledvalue functionand defined by(4)viy=supu∈U⁡Jiy,u.For eachy∈X, we assume the existence of, at least, one optimal control uy∗∈U such that(5)viy=maxu∈U⁡Jiy,u=Jiy,uy∗=∫0∞e-ρitfixyt,uy∗,uy∗tdt,where xy(t,uy) is an admissible trajectory satisfying (2)-(3). We denote by xy∗(t)=xy(t,uy∗) the optimal trajectory.According to [11], the function hy,uyi0,∞→R given by(6)hy,uyiτ≔∫0τfi(xyt;uy,uyte-ρitdt+vixyτ;uye-ρiτis constant and nonincreasing for all t>0 if and only if uy∈U is an optimal control for the initial condition y. Thus, considering uy∗∈U, and derivating hy,uyi with respect to τ, we obtain(7)ρivixy∗τ-fixy∗τ,uy∗τ-∇vixy∗τTgxy∗τ,uy∗τ=0.Furthermore, forτ=0, we have that(8)ρiviy-fiy,uy∗0-∇viyTgy,uy∗0=0.We note that, for any other admissible constant control of the formu(t)=ξ∈U, for all t≥0, as hy,uyi(τ) is not increasing, we have that(9)ddτ∫0τfixyτ,ξe-ρitdt+vxyτ,ξe-ρtτ≤0. Again, for τ=0, we obtain(10)ρiviy-fiy,ξ-∇viyTgy,ξ≤0. So, under previous hypothesis, we conclude that for all y∈X we satisfy the following HJB equation:(11)ρiviy=maxu∈U⁡fiy,u+∇viyTgy,u,anduy∗ is such that(12)uy∗0∈arg⁡maxξ∈U⁡fiy,ξ+∇viyTgy,ξ.We remark that, in general, ifvi is only continuous, such equation needs to be interpreted in terms of viscosity solution of Problem (11). However, when it applies to differential games, there is, so far, no general theorem on the existence or uniqueness of solutions (see [12]).Now, we defineSi:X→Ui calledfeedback-mapper player i=1,…,N, such that(13)Siy∈arg⁡maxu∈U⁡fiy,u+∇ViyTgy,u,wherey∈Rp. We consider an optimal control defined by(14)uy∗t=Sixy∗t,for almost everyt>0.The abovementioned synthesis procedure consists in obtaining, using thefeedback-map,an optimal decision uy∗(t)=Si(xy∗(t)) related to the corresponding optimal trajectory, by solving(15)x˙t=gxt,S1xt,…,SNxtt∈0,∞,and(16)x0=y.So, provided an initial positiony, we say that uy∗,xy∗∈U×C[0,∞),X is an optimal paircontrol-trajectory for every initial condition, and it corresponds to the optimal feedback policy we are trying to estimate.


## 2.3. Feedback-Nash Equilibrium

Now, we adapt the previous expressions in order to get a feedback-Nash equilibrium.To do so, we define afeedback-Nash mapSi:X→Ui, per player i=1,…,N, such that, for all y∈Rp,(17)Siy∈arg⁡maxui∈Ui⁡fiy,ui;u-i+∇ViyTgy,ui;u-i,where the array ui;u-i denotes the pair of the control for player i and the controls associated with the rest of players (denoted by -i).Considering this feedback-Nash map, we apply the synthesis procedure described previously. To do so, we define a feedback N-tupleS≡S1,…,Si-1,Si,Si+1,…,SN and a feedback (N-1)-tuple S-i≡S1,…,Si-1,Si+1,…,SN. Then,S∗≡S∗,1,…,S∗,i-1,S∗,i,S∗,i+1,…,S∗,N is a feedback-Nash equilibrium (FNE) if(18)Jiy;S∗≥Jiy;ui,S∗,-iforeachy∈X,Si:X→Uiandi∈1,…,N,where Ji(y,Si) denotes Ji(y,Si(y)) and ui,S∗,-i is the N-vector of controls obtained by replacing the ith component in Si∗ by ui.Assumingy and S∗,-i are fixed, finding ui by maximizing Jiy;ui,S∗.-i (i.e., solving (18)) can be solved by using dynamic programming (see, for instance, [7, 13]).Finally, we define(19)viy;S∗,-i=maxui∈Ui⁡Jiy;ui,S∗,-i,which is the solution of the following HJB equation (see [3]):(20)ρiviy;S∗,-i=maxui∈Ui⁡fiy,uiS∗,-i+∇Viy;S∗,-iTgy,ui,S∗,-i,where the notationvi.;S∗,-i means that we are finding a Nash equilibrium for player i, by fixing the optimal strategies of the N-1 players. We note that, regarding current literature, apart from zero-sum games, to find relevant cases where (20) is well posed, we need to focus on games in one spatial dimension. For instance, in [14], an existence theorem of Nash equilibria in feedback form, valid for one-space dimensionnoncooperativegames, is given. However, as far as we know, there are no general existence theorems for feedback-Nash equilibrium (for spaces in dimension greater than 1).


## 3. The RaBVItG Algorithm for Solving Deterministic Differential Games

In this section, we describe the numerical implementation of the algorithm RaBVItG used to solve the considered deterministic differential game presented in Section2.2. To do so, we propose a semi-Lagrangian discretization scheme of the HJB equation (see [9]). Then, we describe the general structure of the algorithm used to solve this problem. Finally, we introduce a particular implementation for the case of a feedback-Nash N-players differential game.


### 3.1. Model Discretization

Here, we propose a particular discrete version of (1)-(3).Leth∈(0,1) be a time step and tk=kh, k∈N∪0. First, we aim to approximate xk≈x(tk). To do so, given y∈X, we consider the following discrete approximation of (1)-(2):(21)Jihy,uh=h∑k=0∞e-ρikfixk,uk,where uk=u1,k,…,uN,k∈Uh, with ui,k∈Uih, being Uh=U1h×U2h×⋯×UNh and with the assumption that the controls are constant on tk,tk+1; uh=u1h,…,uNh∈Uh, with uih∈Uih the control for player i, and Uh=U1h×U2h×⋯×UNh, being Uih=uih:uih=ui,0,ui,1,… with ui,k∈Uih.Then, starting fromx0=y, we use a first-order Euler scheme for the state equations:(22)xk+1=xk+hgxk,uk.Next, we discretize the HJB equation (20). To this aim, given y∈RN, we approximate (19) by considering(23)vihy=maxuih∈Uih⁡Jihy,uih,u-ih.Following [9], we obtain a first-order discrete-time HJB equation for player i (i.e., a discrete version of (20)) of the form:(24)vihy=maxuih∈Ui⁡hfiy,uih,u-ih+1-ρihvihy+hgy,uih,u-ih,and we point out that 1-ρih corresponds to the first-order Taylor expansion of e-ρih which is typically used as a discounting factor in discrete dynamic programming (see, i.e., [15]).Now, focusing on the synthesis procedure, we define the following discrete-time versions ofdiscrete feedback N-tuple asSh≡S1h,…,Si-1h,Sih,Si+1h,…,SNh and thediscrete feedback (N-1)-tuple asS-ih≡S1h,…,Si-1h,Si+1h,…,SNh.Note that thediscrete feedback-NashS∗h≡S∗1h,…,S∗Nh satisfies(25)Jihy;S∗h≥Jiy;uih,S∗,-ih,for each y∈X,Sih:X→Ui and i∈1,…,N.Furthermore, according to the definition of thefeedback-map,(26)Sihy;S-i,h∈arg⁡maxuih∈Ui⁡hfiy,uih,S-ih+1-ρihvihz+hgy,uih,S-ih;S-ih.Thus, we obtain(27)vihy=h∑k=0∞e-ρikfixk∗,Shxk∗,where x0∗=y for k∈N and(28)xk+1∗=xk∗+h·xk∗,Shxk∗.However, to determinevih, satisfying (24) is still not always feasible. Thus, we aim to obtain an approximation, denoted by vi,jh, by considering a spatial discretization of (20). To do so, we consider a set of Q∈N arbitrary points yjj=1Q∈X~, with X~ being a closed subset of X. Next, we approximate vi,jhyj;S∗-ih for each i∈1,…,N and j∈1,…,Q. To this aim, let y♯=yj+h·gyj,ui,S∗-ih∈X~ (which are not necessarily in X~). For those points, the HJB equation (24) can be approximated as(29)vi,jhyj;S∗-ih=maxui∈Ui⁡hfiyj,ui,S∗-ih+1-ρihIViy♯;S∗-ih,where IVi· is the approximation of vi,jhy♯;S∗-ih at points y♯ computed by a collocation algorithm (see, for instance, [6]) using a “mesh-free” method based on scattered nodes.More precisely,IVi· is of the form(30)IViy♯=∑j=1PλjΦy♯-yj2,y♯∈X~,where λj∈R, Δ2 is the Euclidean norm and Φr is a real-valued radial basis function (see, for instance, [16]). Here, we use the Gaussian RBF given by Φ(r)=exp⁡-y♯-yj2/2σ2 with σ2>0. In order to determine λj, for j=1,…,Q, we consider that(31)Aλ¯i=Vi,where Ajl=Φyj-yl, for j,l=1,…,Q, λ¯i=λ1,..,λQT and Vi=vihy1,…,vihyQT for i=1,…,N.


### 3.2. Algorithm Structure

In this section, we present the general structure of the algorithm used to solve the problem defined by (1)-(3) and (18) and using the discrete HJB equation (29). This algorithm is based on two main nested loops. It combines a process of the main loop called nvalue iteration(see [13]) with an inner loop calledgame iteration (GI), consisting in a relaxation algorithm to find the proper (convergent) Nash equilibrium for anapproximated value (until reaching convergence of this value).Firstly, before presenting the algorithm, we introduce some useful notations. LetV be the array of values for all the players evaluated in all the original set of points X~ and given by(32)V=v1hy1,…,v1hyQ,…,vNhy1,…,vNhyQT∈RQ·N. We also define Λ as a matrix that stores the controls of each player as follows:(33)Λ=Λ1,…,Λi,…,ΛN∈MQ·N×m, where Λi=ui1,…,uiQT∈MQ×m and Mn×m denotes the set of all real-valued matrices of dimension n×m. We also define(34)Λ-i=Λ1,…,Λi-1,Λi+1,…,ΛN∈MQN-1×m. Additionally, let F be a vector that quantifies the cost for each player at every data point, given by(35)FΛ=f1x,Λ1,…,f1x,Λp,…,fNx,Λ1,…,fNx,ΛpT∈RQ·N.Next, we introduce two operators, T=T1,…,TQ·N:RQ·N→RQ·N and G=G1,…,GQ·N:RQ·N→RQ·N, with(36)TjV;Λ=hFΛj+1-ρihIVΛj,j=1,…,Q·N,and(37)Gi,jΛi;Vi=arg⁡maxΛi∈Ui⁡Jiy,Λi,j,Λ-i,ji=1,…,N,j=1,…,Q·N. In the previous expression, IVΛ=A·λ1¯⋯⋯A·λN¯∈RQ·N is a Q·N interpolation block vector, where λi¯ and A are defined in (31).Secondly, we pretend to solve the following fixed-point schemes:(38)V=TV;Λ,△andΛ=GΛ;V,△,where △ denotesspatiotemporal discretization parameters. More precisely, let s=0,1,… and r=0,1,…,. We define the following processes: (i)
Game iteration: first, for each playeri, we generate Λis+1, a candidate to optimal policy at steps+1 as follows:(39)Λis+1j=θΛisj+1-θarg⁡maxΛi∈Ui⁡hFΛi,Λi-1sj+1-ρihIVirΛi,Λi-1sjj=1,…,Q·Nwithy♯=yj+hg(yj,Λi,Λ-isj), 0<θ<1 a weight coefficient (see, e.g., [17]), and Vir described below.We iterate this process until reaching convergence (i.e., consideringΛs+1-Λs<ϵ1, with ϵ1>0). Once the convergence is reached, we consider the following candidate for the feedback-Nash:(40)Λ∗=Λ1s+1,…,ΛNs+1.Doing so, we obtain atrue (in the sense of convergence) Nash equilibrium for afalse value function (until our value function also converges).(ii)
Value iteration: once we have obtained a candidate for the Nash equilibrium (i.e., the previous game iteration loop ended), we update the value function at stepr+1 as follows:(41)Vir+1j=hFΛ∗j+1-ρihIVirΛ∗j.From a general point of view, our numerical scheme consists in a coupled system of the formUs+1=G(Us;Vr,△), Vr+1=T(Vr;Us+1,△). It can be summarized as(42)Whileerror>ϵ1r=0.Fors=0,1,…Λs+1=GΛs;VrIfΛs+1-Λs<ϵ2BreakForEndforVr+1=TVr;Λs+1error=Vr+1-Vrs=0r=r+1EndWhile.Remark 1.
As done in Section4, we compare our algorithm with a competitive one found in the literature (see [1]), and we briefly recall the former algorithm main structure:(43)Whileerror>ϵ1r=0Vr+1=TVr;Λerror=Vr+1-Vrr=r+1EndWhile.The main difference between both methods is the existence of the step calledgame iteration.


### 3.3. Algorithm Pseudocode

Now, we present in detail a pseudocode version of our algorithm presented in the previous section.(i)
Initialize all parameters and counters:
(a)
Set:N the number of players; 0<h<1 the time step; ρi the discount parameter; τi=1-ρih; 0<θ<1 the update parameter for the relaxation algorithm ingame iteration. Fix It1max∈N and It2max∈N, the maximum number of iterations for the two while loops defined below.(b)
Initialize:r=0 thevalue iterationstep; s=0 thegame iterationstep;the tolerance value for the stopping criteria: tol≈0 and tolG≈0, for Value Iteration and Game Iteration, respectively.(c)
Define bounds for controls and state variables according to the experimental data.(ii)
Build an initial approximation:
(a)
Generate a setX of Q random scattered points, with a uniform distribution, for each player from X~:(44)X=X11,…,X1Q,…,XN1,…,X1Q∈RQ·N(b)
Generate a set of convenient initial values and controls:(45)Vr=V11,…,V1Q,…,VN1,…,VNQ∈RQ·N.If there is no a priori information, they are set to zero.(46)Λs=Λ1,…,Λi,…,ΛN∈MQ·N×m(iii)
We defineA, using radial basis functions (see (30))(47)Ajl=ΦXj-Xl,forj=1,…,Qandl=1,…,Q,
(a)
Determineλ¯it (for i=1,…,N) by collocation (see (31)) satisfying:(48)Aλ¯ir=VirWHILE (dif>tol) and (r<It1max)WHILE (difG>tolG) and (s<It2max)(iv)
Forj=1,…,Q·N, perform:
(a)
Sety=Xj(b)
Set(49)Λsj=Λ1s,…,Λis,…,ΛNsj(c)
Apply the game iteration process, for each playeri=1,…,N:
(1)
Obtain a state variable:y♯=y+hg(y,Λi,Λ-isj).(2)
GetAi=ΦXj-y♯, j=1,…,Q·N,(3)
ComputeΛis+1j=θΛisj+(1-θ)argmaxΛi∈Ui⁡hF(Λi,Λsj)+(1-ρih)λi¯r·Aij(4)
Set(50)Λs+1j=Λ1s+1,…,Λis+1,…,ΛNs+1j(d)
ComputedifG=Λs+1-Λs.
(1)
Sets=s+1END WHILE (v)
DefineΛ∗j=Λ1s+1,…,Λis+1,…,ΛNs+1j (vi)
Sety♯=y+h·g(y,Λ∗j). (vii)
GetAij∗=ΦXj-y♯, j=1,…,Q·N, (viii)
Approximate value function:IΛi∗Virj=λ¯ir·Ai∗j (ix)
Actualize the value function (Value Iteration process):Vit+1j=hF(Λ∗)j+(1-ρih)(λ¯ir·Ai∗)j (x)
Actualizeλ¯ir+1 by collocation(51)Aiλ¯ir+1=Vir+1 (xi)
Computedif=Vr+1-Vr (xii)
r = r + 1 (xiii)
s = 0 .END WHILE (xiv)
In cases for which the algorithm converges, we return the equilibrium policiesΛ1∗,…,ΛN∗ and their associated optimal values: V1∗,…,VN∗, both depending on X. (xv)
By using RBF, we approximate the obtained feedback optimal controls to all values of the state variables which are in the bounded state space defined during the initialization step.Remark 2.
We note that, in the context of the FNDG problems tackled in this work, it is not a trivial task to prove the convergence of the considered semi-Lagrangian discretization to the viscosity solution of the HJB equations (which can be seen as a weak solution). As shown in [11], the definition of viscosity solution is mainly based on comparison principles with upper and lower solutions. However, those techniques are not adapted to the case of vector-valued functions contemplated here (see, e.g., [18]). Additionally, there exist some results in the literature for systems of Hamilton-Jacobi equations (see, for instance, [19]) studying the limits of vanishing viscosity solutions. In practice, it is a quite complex task to prove the existence and uniqueness of those limits. Therefore, in our case, the analysis of the convergence of RaBVItG is only studied through the convergence of the considered numerical experiments. A theoretical analysis of the convergence properties of the proposed method should be performed as a future work.


## 3.1. Model Discretization

Here, we propose a particular discrete version of (1)-(3).Leth∈(0,1) be a time step and tk=kh, k∈N∪0. First, we aim to approximate xk≈x(tk). To do so, given y∈X, we consider the following discrete approximation of (1)-(2):(21)Jihy,uh=h∑k=0∞e-ρikfixk,uk,where uk=u1,k,…,uN,k∈Uh, with ui,k∈Uih, being Uh=U1h×U2h×⋯×UNh and with the assumption that the controls are constant on tk,tk+1; uh=u1h,…,uNh∈Uh, with uih∈Uih the control for player i, and Uh=U1h×U2h×⋯×UNh, being Uih=uih:uih=ui,0,ui,1,… with ui,k∈Uih.Then, starting fromx0=y, we use a first-order Euler scheme for the state equations:(22)xk+1=xk+hgxk,uk.Next, we discretize the HJB equation (20). To this aim, given y∈RN, we approximate (19) by considering(23)vihy=maxuih∈Uih⁡Jihy,uih,u-ih.Following [9], we obtain a first-order discrete-time HJB equation for player i (i.e., a discrete version of (20)) of the form:(24)vihy=maxuih∈Ui⁡hfiy,uih,u-ih+1-ρihvihy+hgy,uih,u-ih,and we point out that 1-ρih corresponds to the first-order Taylor expansion of e-ρih which is typically used as a discounting factor in discrete dynamic programming (see, i.e., [15]).Now, focusing on the synthesis procedure, we define the following discrete-time versions ofdiscrete feedback N-tuple asSh≡S1h,…,Si-1h,Sih,Si+1h,…,SNh and thediscrete feedback (N-1)-tuple asS-ih≡S1h,…,Si-1h,Si+1h,…,SNh.Note that thediscrete feedback-NashS∗h≡S∗1h,…,S∗Nh satisfies(25)Jihy;S∗h≥Jiy;uih,S∗,-ih,for each y∈X,Sih:X→Ui and i∈1,…,N.Furthermore, according to the definition of thefeedback-map,(26)Sihy;S-i,h∈arg⁡maxuih∈Ui⁡hfiy,uih,S-ih+1-ρihvihz+hgy,uih,S-ih;S-ih.Thus, we obtain(27)vihy=h∑k=0∞e-ρikfixk∗,Shxk∗,where x0∗=y for k∈N and(28)xk+1∗=xk∗+h·xk∗,Shxk∗.However, to determinevih, satisfying (24) is still not always feasible. Thus, we aim to obtain an approximation, denoted by vi,jh, by considering a spatial discretization of (20). To do so, we consider a set of Q∈N arbitrary points yjj=1Q∈X~, with X~ being a closed subset of X. Next, we approximate vi,jhyj;S∗-ih for each i∈1,…,N and j∈1,…,Q. To this aim, let y♯=yj+h·gyj,ui,S∗-ih∈X~ (which are not necessarily in X~). For those points, the HJB equation (24) can be approximated as(29)vi,jhyj;S∗-ih=maxui∈Ui⁡hfiyj,ui,S∗-ih+1-ρihIViy♯;S∗-ih,where IVi· is the approximation of vi,jhy♯;S∗-ih at points y♯ computed by a collocation algorithm (see, for instance, [6]) using a “mesh-free” method based on scattered nodes.More precisely,IVi· is of the form(30)IViy♯=∑j=1PλjΦy♯-yj2,y♯∈X~,where λj∈R, Δ2 is the Euclidean norm and Φr is a real-valued radial basis function (see, for instance, [16]). Here, we use the Gaussian RBF given by Φ(r)=exp⁡-y♯-yj2/2σ2 with σ2>0. In order to determine λj, for j=1,…,Q, we consider that(31)Aλ¯i=Vi,where Ajl=Φyj-yl, for j,l=1,…,Q, λ¯i=λ1,..,λQT and Vi=vihy1,…,vihyQT for i=1,…,N.


## 3.2. Algorithm Structure

In this section, we present the general structure of the algorithm used to solve the problem defined by (1)-(3) and (18) and using the discrete HJB equation (29). This algorithm is based on two main nested loops. It combines a process of the main loop called nvalue iteration(see [13]) with an inner loop calledgame iteration (GI), consisting in a relaxation algorithm to find the proper (convergent) Nash equilibrium for anapproximated value (until reaching convergence of this value).Firstly, before presenting the algorithm, we introduce some useful notations. LetV be the array of values for all the players evaluated in all the original set of points X~ and given by(32)V=v1hy1,…,v1hyQ,…,vNhy1,…,vNhyQT∈RQ·N. We also define Λ as a matrix that stores the controls of each player as follows:(33)Λ=Λ1,…,Λi,…,ΛN∈MQ·N×m, where Λi=ui1,…,uiQT∈MQ×m and Mn×m denotes the set of all real-valued matrices of dimension n×m. We also define(34)Λ-i=Λ1,…,Λi-1,Λi+1,…,ΛN∈MQN-1×m. Additionally, let F be a vector that quantifies the cost for each player at every data point, given by(35)FΛ=f1x,Λ1,…,f1x,Λp,…,fNx,Λ1,…,fNx,ΛpT∈RQ·N.Next, we introduce two operators, T=T1,…,TQ·N:RQ·N→RQ·N and G=G1,…,GQ·N:RQ·N→RQ·N, with(36)TjV;Λ=hFΛj+1-ρihIVΛj,j=1,…,Q·N,and(37)Gi,jΛi;Vi=arg⁡maxΛi∈Ui⁡Jiy,Λi,j,Λ-i,ji=1,…,N,j=1,…,Q·N. In the previous expression, IVΛ=A·λ1¯⋯⋯A·λN¯∈RQ·N is a Q·N interpolation block vector, where λi¯ and A are defined in (31).Secondly, we pretend to solve the following fixed-point schemes:(38)V=TV;Λ,△andΛ=GΛ;V,△,where △ denotesspatiotemporal discretization parameters. More precisely, let s=0,1,… and r=0,1,…,. We define the following processes: (i)
Game iteration: first, for each playeri, we generate Λis+1, a candidate to optimal policy at steps+1 as follows:(39)Λis+1j=θΛisj+1-θarg⁡maxΛi∈Ui⁡hFΛi,Λi-1sj+1-ρihIVirΛi,Λi-1sjj=1,…,Q·Nwithy♯=yj+hg(yj,Λi,Λ-isj), 0<θ<1 a weight coefficient (see, e.g., [17]), and Vir described below.We iterate this process until reaching convergence (i.e., consideringΛs+1-Λs<ϵ1, with ϵ1>0). Once the convergence is reached, we consider the following candidate for the feedback-Nash:(40)Λ∗=Λ1s+1,…,ΛNs+1.Doing so, we obtain atrue (in the sense of convergence) Nash equilibrium for afalse value function (until our value function also converges).(ii)
Value iteration: once we have obtained a candidate for the Nash equilibrium (i.e., the previous game iteration loop ended), we update the value function at stepr+1 as follows:(41)Vir+1j=hFΛ∗j+1-ρihIVirΛ∗j.From a general point of view, our numerical scheme consists in a coupled system of the formUs+1=G(Us;Vr,△), Vr+1=T(Vr;Us+1,△). It can be summarized as(42)Whileerror>ϵ1r=0.Fors=0,1,…Λs+1=GΛs;VrIfΛs+1-Λs<ϵ2BreakForEndforVr+1=TVr;Λs+1error=Vr+1-Vrs=0r=r+1EndWhile.Remark 1.
As done in Section4, we compare our algorithm with a competitive one found in the literature (see [1]), and we briefly recall the former algorithm main structure:(43)Whileerror>ϵ1r=0Vr+1=TVr;Λerror=Vr+1-Vrr=r+1EndWhile.The main difference between both methods is the existence of the step calledgame iteration.


## 3.3. Algorithm Pseudocode

Now, we present in detail a pseudocode version of our algorithm presented in the previous section.(i)
Initialize all parameters and counters:
(a)
Set:N the number of players; 0<h<1 the time step; ρi the discount parameter; τi=1-ρih; 0<θ<1 the update parameter for the relaxation algorithm ingame iteration. Fix It1max∈N and It2max∈N, the maximum number of iterations for the two while loops defined below.(b)
Initialize:r=0 thevalue iterationstep; s=0 thegame iterationstep;the tolerance value for the stopping criteria: tol≈0 and tolG≈0, for Value Iteration and Game Iteration, respectively.(c)
Define bounds for controls and state variables according to the experimental data.(ii)
Build an initial approximation:
(a)
Generate a setX of Q random scattered points, with a uniform distribution, for each player from X~:(44)X=X11,…,X1Q,…,XN1,…,X1Q∈RQ·N(b)
Generate a set of convenient initial values and controls:(45)Vr=V11,…,V1Q,…,VN1,…,VNQ∈RQ·N.If there is no a priori information, they are set to zero.(46)Λs=Λ1,…,Λi,…,ΛN∈MQ·N×m(iii)
We defineA, using radial basis functions (see (30))(47)Ajl=ΦXj-Xl,forj=1,…,Qandl=1,…,Q,
(a)
Determineλ¯it (for i=1,…,N) by collocation (see (31)) satisfying:(48)Aλ¯ir=VirWHILE (dif>tol) and (r<It1max)WHILE (difG>tolG) and (s<It2max)(iv)
Forj=1,…,Q·N, perform:
(a)
Sety=Xj(b)
Set(49)Λsj=Λ1s,…,Λis,…,ΛNsj(c)
Apply the game iteration process, for each playeri=1,…,N:
(1)
Obtain a state variable:y♯=y+hg(y,Λi,Λ-isj).(2)
GetAi=ΦXj-y♯, j=1,…,Q·N,(3)
ComputeΛis+1j=θΛisj+(1-θ)argmaxΛi∈Ui⁡hF(Λi,Λsj)+(1-ρih)λi¯r·Aij(4)
Set(50)Λs+1j=Λ1s+1,…,Λis+1,…,ΛNs+1j(d)
ComputedifG=Λs+1-Λs.
(1)
Sets=s+1END WHILE (v)
DefineΛ∗j=Λ1s+1,…,Λis+1,…,ΛNs+1j (vi)
Sety♯=y+h·g(y,Λ∗j). (vii)
GetAij∗=ΦXj-y♯, j=1,…,Q·N, (viii)
Approximate value function:IΛi∗Virj=λ¯ir·Ai∗j (ix)
Actualize the value function (Value Iteration process):Vit+1j=hF(Λ∗)j+(1-ρih)(λ¯ir·Ai∗)j (x)
Actualizeλ¯ir+1 by collocation(51)Aiλ¯ir+1=Vir+1 (xi)
Computedif=Vr+1-Vr (xii)
r = r + 1 (xiii)
s = 0 .END WHILE (xiv)
In cases for which the algorithm converges, we return the equilibrium policiesΛ1∗,…,ΛN∗ and their associated optimal values: V1∗,…,VN∗, both depending on X. (xv)
By using RBF, we approximate the obtained feedback optimal controls to all values of the state variables which are in the bounded state space defined during the initialization step.Remark 2.
We note that, in the context of the FNDG problems tackled in this work, it is not a trivial task to prove the convergence of the considered semi-Lagrangian discretization to the viscosity solution of the HJB equations (which can be seen as a weak solution). As shown in [11], the definition of viscosity solution is mainly based on comparison principles with upper and lower solutions. However, those techniques are not adapted to the case of vector-valued functions contemplated here (see, e.g., [18]). Additionally, there exist some results in the literature for systems of Hamilton-Jacobi equations (see, for instance, [19]) studying the limits of vanishing viscosity solutions. In practice, it is a quite complex task to prove the existence and uniqueness of those limits. Therefore, in our case, the analysis of the convergence of RaBVItG is only studied through the convergence of the considered numerical experiments. A theoretical analysis of the convergence properties of the proposed method should be performed as a future work.


## 4. Numerical Experiments

In Sections1 and 2, we have presented and developed the RaBVItG algorithm to solve differential games involving N-players pursuing a feedback-Nash. Now, in the current section, we apply the algorithm to several benchmark problems in order to evaluate its main properties.Remark 3.
In a recent survey (see [2]), the existing numerical methods in the literature focus on some feasible subclass of games, regarding the ones we are dealing with in this work: two-person zero-sum games (see [3]) which can deal with the lack of differentiability of the value function by using the viscosity solution since this game can be reduced to a one-dimensional game. Additionally, the other advances are in line with linear quadratic models (see [4]), avoiding the interesting nonlinearity in the dynamic system or running cost. According to the authors, computational differential games is an area that needs to grow up.Here, we solve different test problems, selected from [1, 4].In order to check the efficiency of our approach, we compute (i) the error of the method with respect to the analytic solution; and (ii) the CPU time and the number of iterations. Furthermore, we also compare our algorithm performances by comparing its results to the one returned by the method proposed in [1], which is called in this work CCF (i.e., from the initials of the authors). To our knowledge, at this moment, it is the only algorithm found in the literature designed to solve similar problems. Both algorithms are based on the semi-Lagrangian discretization of the problem and the design of value iterations as the main loop. However, they exhibit several relevant differences, such as the following:(i)
RaBVItG is a meshfree algorithm and, thus,
(1)
only evaluated at some points of the state space,(2)
uses a RBF method to approximate functions,(3)
does not require a discretization of the control space, so that we use a gradient algorithm to find the critical points at each algorithm step.(ii)
RaBVItG incorporates an iteration inner loop, called game iteration, where players alternately update their strategies by averaging their current strategy with the best response to the other player’s current strategy. The main advantage is that it finds the proper Nash equilibrium for an approximated value until reaching convergence in the value iteration loop. Our hypothesis is that this step increases the accuracy and, thus, reduces the approximation error.


### 4.1. Numerical Test 1 (LQ, Scalar)

This example is introduced in [4]. In this case, we consider a two-player (N=2) scalar (X⊆R) differential game, where m=2. Let the following quadratic cost functional to be minimized, with i=1,2:(52)Jiy,u1,u2=∫0∞qix2t+riuit2dt,subject to the following linear dynamic system:(53)x˙t=axt+b1u1t+b2u2t,x0=y,y∈X.According to [4], the analytic feedback-Nash equilibria for this game are defined by the following linear expressions:(54)ui∗t=fi∗x∗t≡-birikix∗t, where k1,k2 are the equilibria for the following system of Ricatti ordinary differential equations, derived from the solution of the corresponding HJB equations: (55)k˙1=-b12r1k12-2b22r2k1k2+2ak1+q1,k˙2=-b22r2k22-2b12r1k1k2+2ak2+q2.Furthermore, an additional condition to be a stable equilibrium is to satisfy the following inequality:(56)a-b12r1k1-b22r2k2<0,where (56) divides the plane into “stable” and “unstable” regions. So, all the feedback-Nash equilibria are obtained as the intersection points of both hyperbolas in the “stable” region.Now, we consider the following particular parameter values:a=3,q1=q2=1,b1=b2=2,r1=r2=1 (see [4], p. 396). Equation (55) becomes(57)k˙1=-4k12-8k1k2+6k1+1,k˙2=-4k22-8k1k2+6k2+1.We obtain those equilibria values:(58)k11=1/2k21=1k12=1k22=1/2k13=0.7287k23=0.7287k14=-0.2287k24=-0.2287. However, according to Theorem  8.22 in [4], since we have here three feedback-Nash equilibria, then two of them are stable nodes and one is a saddle point. We note that the pair (k1(4),k2(4)) does not satisfy condition (56). In Figure 1, we show the phase diagram with the three stable equilibria, where ki(1) and ki(3) are stable nodes and ki(2) is a saddle point. Finally, ki(4) is unstable.Figure 1
Numerical Test 1: phase diagram of System (55).Next, we analyze the performance of our algorithm to solve this problem and compare it with CCF, mentioned previously. To do so, we choose three different starting points which are in the basin of attraction of the different stable feedback-Nash equlibria. As it can be seen in Figure1, equilibria k(2) have a symmetric behavior with respect to k(1) and, so, the results will remain the same. Thus, we only analyze two of the three equilibria. In order to analyze the performance of both algorithms, we run them three times starting from the following initial points: the corresponding equilibrium and points in the attraction domain of each equilibrium (see Figure 1).We define our state variable as0≤x≤1 and we discretize it by using Q=9 points, so X=[0,0.125,0.25,0.375,0.5,0.625,0.75,0.875,1]. Regarding the control space, it is bounded by considering values in [-5,5]. RaBVItG does not need any grid in the control space. In the case of CCF, we have generated a uniform grid of 20 elements from the defined space.In Table1, we summarize the results obtained with both algorithms (i.e., RaBVItG and CCF) and their performances to reach the two equilibria of interest. For the three different starting equilibria for each player, and for each starting point, we report the errors values: (59)LerrorA∞=max⁡-2x-uA1,-x-uA2,where uAi is the vector of solutions (for i=1,2) obtained by the algorithms trying to achieve the equilibrium A=-2x,-x, and(60)LerrorB∞=max⁡-1.4574x-uB1,-1.4574x-uB2, where uBi are the numerical solutions obtained for the equilibrium B=-1.4574x,-1574x, starting from the different initial values.Table 1
Numerical Test 1 LQ, Scalar: results obtained with RaBVItG and CCF. In each cell, we report two values: (i) the upper number representsLerrorA∞ or LerrorB∞ depending on the case; (ii) and the bottom number represents the number of iterations the algorithm requires for convergence.
A = - 2 x , - x B = - 1.4574 x , - 1.4574 x Time Step Algorithm - 2.75 x , - 0.5 x - 2 x , - x - 4 x , - 0.25 x - x , - x - 1.4574 x , - 1.4574 x - 4 x , - 4 x  h=0.1 RaBVitG 0.1427 14 0.3212 10 0.2227 26 0.1039 10 0.1025 10 0.0927 5 CCF 0.2832 22 0.2832 10 0.2331 24 0.1161 9 0.1161 6 0.1161 12  h=0.05 RaBVitG 0.0809 14 0.0159 4 0.0793 74 0.0109 19 0.0104 19 0.0089 7 CCF 0.1955 67 0.2206 22 0.1704 52 0.0535 20 0.0535 8 0.0535 12  h=0.01 RaBVitG 0.0267 847 0.0139 136 0.0267 912 0.009 57 0.009 56 0.0068 23 CCF 0.0451 440 0.0450 68 0.0840 584 0.0284 56 0.0284 25 0.0284 19As it can be seen in Table1, both algorithms show a similar qualitative pattern: the smaller the value of h, the smaller the error and the higher the number of iterations needed. It is interesting to remark that our algorithm, RaBVitG, obtains lower errors but requires, in general, more iterations than CCF. However, as we point out in the next experiments, the computational time of RaBVitG is dramatically lower than the one of the CCF algorithm. This tends to show that RaBVitG seems to be more efficient and accurate than CCF.On the other hand, Figure2 shows how both algorithms converge to different equilibrium, depending on the starting point. To perform this experiment, we have run both algorithms using a mesh of 162 initial conditions from the set -4,-3.5,…,-0.5×-4,-3.5,…,-0.5. As we can see in Figure 1, for each plot, the white region represents the initial points that converge to the corresponding equilibrium (i.e., they satisfy the tolerance criterium of convergence with respect to the point represented by a circle). Additionally, the grey region is composed by the rest of initial points that do not satisfy this tolerance criterium. Thus, the white region can be interpreted as a numerical approximation of the basin of attraction of each equilibrium.Figure 2
Numerical Test 1 (LQ, Scalar): Convergence of the algorithm towards Nash equilibria: (top) RaBVItG and (bottom) CCF).f1,f2 represent the starting points ui0(t)=fi0x(t), i=1,2. The white region means that the algorithm starting from a point inside this area tends to the considered equilibrium (represented by a circle). The grey region means that the algorithm starting from this area tends to another equilibrium.Comparing both algorithms, RaBVItG has less problems than CCF to identify the correct equilibrium in the area(f1,f2)∈[0,-1]×[0,-1]. Focusing on the saddle equilibrium, the stable manifold is situated on the bisector, so any trajectory starting out of this bisector should converge to the closest stable node (i.e., (-2,-1) or (-1,-2)). As it can be seen in Figure 2, CCF exhibits more wrong equilibria than RaBVItG. When starting from the attraction region of a stable equilibrium, for some initial points, both algorithms converge to the saddle point. However, RaBVItG performs better than CCF, since the convergence band around this equilibrium is finer.


### 4.2. Numerical Test 2 (Non-LQ, Scalar)

The next experiment was proposed in [1]. Again, we present a two-player scalar differential game, with the following cost functional to minimize, with i=1,2:(61)Jiy,u1,u2=∫0∞sixt+uit22dt,subject to the dynamic system(62)x˙t=u1t+u2t,x0=y,with y∈R.According to [1], we obtain the analytical value functions per player attained to their feedback-Nash equilibrium,(63)v1∗x=s1x-s1s2-s122and(64)v2∗x=s2x-s1s2-s222.Here, we considers1=-1, s2=2. It is straightforward to deduce, once we know the true value function and using the corresponding HJB equations, that the optimal feedback-Nash policies remain constant for all x(t) values:(65)u1∗=1,u2∗=-2.Focusing on the RaBVItG and CCF configurations, we use the following specifications: the admissible values for the state variable are in[-50,50]. We discretize it using Q=8 equispaced points: x=-50,-37.5,-25,-12.5,0,12.5,25,37.5,50. The set of admissible controls is in [-5,5]. Additionally, only in the case of CCF, we discretize the control variable using a grid of M equispaced points. We choose M=50,150,400 in order to compare the effect of improving the control space discretization. We use an interpolator based on splines to approximate the value function in the corresponding algorithm step.In Figure3, we depict the different errors in infinity-norm between the numerical feedback-Nash policies and the real ones (65). As we can see in this figure, the greater M is, the better the approximation to the analytical solution is. Indeed, in the case of M=400, CCF produces slightly better results than RaBVItG. However, there is an important trade-off between improving the error measures and the number of points that CCF needs. Finally, reported in Table 2, our algorithm is clearly more time efficient and requires, approximately, the same quantity of iterations in value than CCF.Table 2
Numerical Test 2 (non LQ, scalar): performance comparison between RaBVitG and CCF. The upper number reports the CPU time (s) and the bottom number the number of iterations.
h = 0.7 h = 0.25 h = 0.1 h = 0.05 h = 0.01 RaBVitG - 0.048 3 0.052 8 0.054 27 0.062 68 0.075 501  CCF M = 50 1.049 3 6.036 8 24.408 27 65.3135 69 523.53 507 M = 150 3.30 3 14.87 8 61.03 27 200.8 68 1702.1 503 M = 400 3.37 3 19.77 8 81.2 27 214.9 68 1763.3 502Figure 3
Numerical Test 2 (non-LQ, scalar):L∞ errors between numerical and analytical feedback-Nash policies obtained considering RaBVItG and CCF with M=50,150,400.


### 4.3. Numerical Test 3:(LQ, Scalar, Three Players)

We consider the following three-player minimization differential game, fori={1,2,3}:(66)Jiy,u1,u2,u3=∫0∞qix2t+riuit2dt,subject to the following dynamic system:(67)x˙t=axt+b1u1t+b2u2t+b3u3t,x0=y, with y∈R. Additionally, the considered parameter values are a=-1, b1=1, b2=1, b3=0.5, q1=2, r1=1, q2=2, r2=2, q3=1, and r3=1.According to [4], this model has only one feedback-Nash equilibrium, which consists in (68)u1∗t=0.6202xtu2∗t=-0.2805xtu3∗t=0.0436xt.We set our admissible values for the state variable in[0,1] and we discretize it by using Q=9 points, so X=[0,0.125,0.25,0.375,0.5,0.625,0.75,0.875,1]. Regarding the control space, we consider values in [-2,2]. For the CCF algorithm, again we discretize the control variable using a grid of M equally spaced points. As previously, we choose M=50,150,400, and we use an interpolator based on splines.As we can see in Table3, RaBVItG again exhibits better computational times, showing the greatest efficiency for all the different time discretization values. Furthermore, in Figure 4, we observe that, even with M=400, RaBVItG produces errors smaller than using CCF. The difference with respect to Test 2 is that the solution, in this case, is not a constant one, being a linear function.Table 3
Numerical Test 3 (LQ, Scalar, Three Players): performance comparison between RaBVitG and CCF. The upper number reports the CPU time (s) and the bottom number the number of iterations.
h = 0.7 h = 0.25 h = 0.1 h = 0.05 h = 0.01 RaBVitG - 0.0049 8 0.0050 17 0.0058 36 0.0048 55 0.0081 215  CCF M = 50 0.8453 4 1.0181 5 3.2296 14 9.003 36 177.2065 668 M = 150 1.5780 4 2.2935 5 8.4177 13 22.1266 31 180.2918 238 M = 400 4.30 4 6.65 5 22.72 13 61.06 31 477.56 234Figure 4
Numerical Test 3 (LQ, scalar, three players):L∞ errors between numerical and analytical feedback-Nash policies obtained considering RaBVItG and CCF with M={50,150,400}.


### 4.4. Numerical Test 4: (Nonscalar, Two Players)

This last minimization differential game is interesting as is not a scalar one, because its state variable has two dimensions. The running-cost fori={1,2} is(69)Jiy,u1,u2=∫0∞e-tx12t+x2t2dtifx12t+x2t2>10otherwisesubject to the dynamic system(70)x˙1t=u2tx˙2t=u1tx0=y.This test case was introduced in [1] but, unfortunately, the authors do not provide any analytical solution to this differential game in order to compare with. However, they show a figure with a value function comparable to the ones reported in Figure 5. We run CCF and RaBVItG using different mesh sizes (in the case of CCF, M1=202,302,402,502) and meshfree points (M2=20,30,40,50) in the case of RaBVItG. In Figure 5, we represent the different value functions for Player 1 (Player 2 has the same result since it is a symmetric game). The fact that CCF requires a mesh for the state variables (and controls) clearly affects the CPU time required to achieve convergence. Indeed, in order to compare the solution of RaBVItG with CCF algorithms, once we obtain the value function per player (as a M2- dimensional array), we interpolate our solution into the mesh defined in CCF using, again, RBF methods. Both algorithms converge to the same solution as the mesh size increases. However, as shown in Table 4, CCF method needs, in general, less value iterations but is much computationally expensive than RaBVItG costly (around 10 times greater).Table 4
Numerical Test 4 (non Scalar, Two Players): performance comparison between RaBVitG and CCF: mesh sizeM, CPU time (s) and number of iterations (Its).
CCF RaBVitG size 2 0 2 3 0 2 4 0 2 5 0 2 20 30 40 50  CPU time 2140 5140.3 9623.6 12854 479.11 603.08 734.8 1272.8  Its 333 287 266 248 542 457 397 275Figure 5
Numerical Test 4 (nonscalar, two players): value functions for player 1 returned by RaBVItG and CCF.


## 4.1. Numerical Test 1 (LQ, Scalar)

This example is introduced in [4]. In this case, we consider a two-player (N=2) scalar (X⊆R) differential game, where m=2. Let the following quadratic cost functional to be minimized, with i=1,2:(52)Jiy,u1,u2=∫0∞qix2t+riuit2dt,subject to the following linear dynamic system:(53)x˙t=axt+b1u1t+b2u2t,x0=y,y∈X.According to [4], the analytic feedback-Nash equilibria for this game are defined by the following linear expressions:(54)ui∗t=fi∗x∗t≡-birikix∗t, where k1,k2 are the equilibria for the following system of Ricatti ordinary differential equations, derived from the solution of the corresponding HJB equations: (55)k˙1=-b12r1k12-2b22r2k1k2+2ak1+q1,k˙2=-b22r2k22-2b12r1k1k2+2ak2+q2.Furthermore, an additional condition to be a stable equilibrium is to satisfy the following inequality:(56)a-b12r1k1-b22r2k2<0,where (56) divides the plane into “stable” and “unstable” regions. So, all the feedback-Nash equilibria are obtained as the intersection points of both hyperbolas in the “stable” region.Now, we consider the following particular parameter values:a=3,q1=q2=1,b1=b2=2,r1=r2=1 (see [4], p. 396). Equation (55) becomes(57)k˙1=-4k12-8k1k2+6k1+1,k˙2=-4k22-8k1k2+6k2+1.We obtain those equilibria values:(58)k11=1/2k21=1k12=1k22=1/2k13=0.7287k23=0.7287k14=-0.2287k24=-0.2287. However, according to Theorem  8.22 in [4], since we have here three feedback-Nash equilibria, then two of them are stable nodes and one is a saddle point. We note that the pair (k1(4),k2(4)) does not satisfy condition (56). In Figure 1, we show the phase diagram with the three stable equilibria, where ki(1) and ki(3) are stable nodes and ki(2) is a saddle point. Finally, ki(4) is unstable.Figure 1
Numerical Test 1: phase diagram of System (55).Next, we analyze the performance of our algorithm to solve this problem and compare it with CCF, mentioned previously. To do so, we choose three different starting points which are in the basin of attraction of the different stable feedback-Nash equlibria. As it can be seen in Figure1, equilibria k(2) have a symmetric behavior with respect to k(1) and, so, the results will remain the same. Thus, we only analyze two of the three equilibria. In order to analyze the performance of both algorithms, we run them three times starting from the following initial points: the corresponding equilibrium and points in the attraction domain of each equilibrium (see Figure 1).We define our state variable as0≤x≤1 and we discretize it by using Q=9 points, so X=[0,0.125,0.25,0.375,0.5,0.625,0.75,0.875,1]. Regarding the control space, it is bounded by considering values in [-5,5]. RaBVItG does not need any grid in the control space. In the case of CCF, we have generated a uniform grid of 20 elements from the defined space.In Table1, we summarize the results obtained with both algorithms (i.e., RaBVItG and CCF) and their performances to reach the two equilibria of interest. For the three different starting equilibria for each player, and for each starting point, we report the errors values: (59)LerrorA∞=max⁡-2x-uA1,-x-uA2,where uAi is the vector of solutions (for i=1,2) obtained by the algorithms trying to achieve the equilibrium A=-2x,-x, and(60)LerrorB∞=max⁡-1.4574x-uB1,-1.4574x-uB2, where uBi are the numerical solutions obtained for the equilibrium B=-1.4574x,-1574x, starting from the different initial values.Table 1
Numerical Test 1 LQ, Scalar: results obtained with RaBVItG and CCF. In each cell, we report two values: (i) the upper number representsLerrorA∞ or LerrorB∞ depending on the case; (ii) and the bottom number represents the number of iterations the algorithm requires for convergence.
A = - 2 x , - x B = - 1.4574 x , - 1.4574 x Time Step Algorithm - 2.75 x , - 0.5 x - 2 x , - x - 4 x , - 0.25 x - x , - x - 1.4574 x , - 1.4574 x - 4 x , - 4 x  h=0.1 RaBVitG 0.1427 14 0.3212 10 0.2227 26 0.1039 10 0.1025 10 0.0927 5 CCF 0.2832 22 0.2832 10 0.2331 24 0.1161 9 0.1161 6 0.1161 12  h=0.05 RaBVitG 0.0809 14 0.0159 4 0.0793 74 0.0109 19 0.0104 19 0.0089 7 CCF 0.1955 67 0.2206 22 0.1704 52 0.0535 20 0.0535 8 0.0535 12  h=0.01 RaBVitG 0.0267 847 0.0139 136 0.0267 912 0.009 57 0.009 56 0.0068 23 CCF 0.0451 440 0.0450 68 0.0840 584 0.0284 56 0.0284 25 0.0284 19As it can be seen in Table1, both algorithms show a similar qualitative pattern: the smaller the value of h, the smaller the error and the higher the number of iterations needed. It is interesting to remark that our algorithm, RaBVitG, obtains lower errors but requires, in general, more iterations than CCF. However, as we point out in the next experiments, the computational time of RaBVitG is dramatically lower than the one of the CCF algorithm. This tends to show that RaBVitG seems to be more efficient and accurate than CCF.On the other hand, Figure2 shows how both algorithms converge to different equilibrium, depending on the starting point. To perform this experiment, we have run both algorithms using a mesh of 162 initial conditions from the set -4,-3.5,…,-0.5×-4,-3.5,…,-0.5. As we can see in Figure 1, for each plot, the white region represents the initial points that converge to the corresponding equilibrium (i.e., they satisfy the tolerance criterium of convergence with respect to the point represented by a circle). Additionally, the grey region is composed by the rest of initial points that do not satisfy this tolerance criterium. Thus, the white region can be interpreted as a numerical approximation of the basin of attraction of each equilibrium.Figure 2
Numerical Test 1 (LQ, Scalar): Convergence of the algorithm towards Nash equilibria: (top) RaBVItG and (bottom) CCF).f1,f2 represent the starting points ui0(t)=fi0x(t), i=1,2. The white region means that the algorithm starting from a point inside this area tends to the considered equilibrium (represented by a circle). The grey region means that the algorithm starting from this area tends to another equilibrium.Comparing both algorithms, RaBVItG has less problems than CCF to identify the correct equilibrium in the area(f1,f2)∈[0,-1]×[0,-1]. Focusing on the saddle equilibrium, the stable manifold is situated on the bisector, so any trajectory starting out of this bisector should converge to the closest stable node (i.e., (-2,-1) or (-1,-2)). As it can be seen in Figure 2, CCF exhibits more wrong equilibria than RaBVItG. When starting from the attraction region of a stable equilibrium, for some initial points, both algorithms converge to the saddle point. However, RaBVItG performs better than CCF, since the convergence band around this equilibrium is finer.


## 4.2. Numerical Test 2 (Non-LQ, Scalar)

The next experiment was proposed in [1]. Again, we present a two-player scalar differential game, with the following cost functional to minimize, with i=1,2:(61)Jiy,u1,u2=∫0∞sixt+uit22dt,subject to the dynamic system(62)x˙t=u1t+u2t,x0=y,with y∈R.According to [1], we obtain the analytical value functions per player attained to their feedback-Nash equilibrium,(63)v1∗x=s1x-s1s2-s122and(64)v2∗x=s2x-s1s2-s222.Here, we considers1=-1, s2=2. It is straightforward to deduce, once we know the true value function and using the corresponding HJB equations, that the optimal feedback-Nash policies remain constant for all x(t) values:(65)u1∗=1,u2∗=-2.Focusing on the RaBVItG and CCF configurations, we use the following specifications: the admissible values for the state variable are in[-50,50]. We discretize it using Q=8 equispaced points: x=-50,-37.5,-25,-12.5,0,12.5,25,37.5,50. The set of admissible controls is in [-5,5]. Additionally, only in the case of CCF, we discretize the control variable using a grid of M equispaced points. We choose M=50,150,400 in order to compare the effect of improving the control space discretization. We use an interpolator based on splines to approximate the value function in the corresponding algorithm step.In Figure3, we depict the different errors in infinity-norm between the numerical feedback-Nash policies and the real ones (65). As we can see in this figure, the greater M is, the better the approximation to the analytical solution is. Indeed, in the case of M=400, CCF produces slightly better results than RaBVItG. However, there is an important trade-off between improving the error measures and the number of points that CCF needs. Finally, reported in Table 2, our algorithm is clearly more time efficient and requires, approximately, the same quantity of iterations in value than CCF.Table 2
Numerical Test 2 (non LQ, scalar): performance comparison between RaBVitG and CCF. The upper number reports the CPU time (s) and the bottom number the number of iterations.
h = 0.7 h = 0.25 h = 0.1 h = 0.05 h = 0.01 RaBVitG - 0.048 3 0.052 8 0.054 27 0.062 68 0.075 501  CCF M = 50 1.049 3 6.036 8 24.408 27 65.3135 69 523.53 507 M = 150 3.30 3 14.87 8 61.03 27 200.8 68 1702.1 503 M = 400 3.37 3 19.77 8 81.2 27 214.9 68 1763.3 502Figure 3
Numerical Test 2 (non-LQ, scalar):L∞ errors between numerical and analytical feedback-Nash policies obtained considering RaBVItG and CCF with M=50,150,400.


## 4.3. Numerical Test 3:(LQ, Scalar, Three Players)

We consider the following three-player minimization differential game, fori={1,2,3}:(66)Jiy,u1,u2,u3=∫0∞qix2t+riuit2dt,subject to the following dynamic system:(67)x˙t=axt+b1u1t+b2u2t+b3u3t,x0=y, with y∈R. Additionally, the considered parameter values are a=-1, b1=1, b2=1, b3=0.5, q1=2, r1=1, q2=2, r2=2, q3=1, and r3=1.According to [4], this model has only one feedback-Nash equilibrium, which consists in (68)u1∗t=0.6202xtu2∗t=-0.2805xtu3∗t=0.0436xt.We set our admissible values for the state variable in[0,1] and we discretize it by using Q=9 points, so X=[0,0.125,0.25,0.375,0.5,0.625,0.75,0.875,1]. Regarding the control space, we consider values in [-2,2]. For the CCF algorithm, again we discretize the control variable using a grid of M equally spaced points. As previously, we choose M=50,150,400, and we use an interpolator based on splines.As we can see in Table3, RaBVItG again exhibits better computational times, showing the greatest efficiency for all the different time discretization values. Furthermore, in Figure 4, we observe that, even with M=400, RaBVItG produces errors smaller than using CCF. The difference with respect to Test 2 is that the solution, in this case, is not a constant one, being a linear function.Table 3
Numerical Test 3 (LQ, Scalar, Three Players): performance comparison between RaBVitG and CCF. The upper number reports the CPU time (s) and the bottom number the number of iterations.
h = 0.7 h = 0.25 h = 0.1 h = 0.05 h = 0.01 RaBVitG - 0.0049 8 0.0050 17 0.0058 36 0.0048 55 0.0081 215  CCF M = 50 0.8453 4 1.0181 5 3.2296 14 9.003 36 177.2065 668 M = 150 1.5780 4 2.2935 5 8.4177 13 22.1266 31 180.2918 238 M = 400 4.30 4 6.65 5 22.72 13 61.06 31 477.56 234Figure 4
Numerical Test 3 (LQ, scalar, three players):L∞ errors between numerical and analytical feedback-Nash policies obtained considering RaBVItG and CCF with M={50,150,400}.


## 4.4. Numerical Test 4: (Nonscalar, Two Players)

This last minimization differential game is interesting as is not a scalar one, because its state variable has two dimensions. The running-cost fori={1,2} is(69)Jiy,u1,u2=∫0∞e-tx12t+x2t2dtifx12t+x2t2>10otherwisesubject to the dynamic system(70)x˙1t=u2tx˙2t=u1tx0=y.This test case was introduced in [1] but, unfortunately, the authors do not provide any analytical solution to this differential game in order to compare with. However, they show a figure with a value function comparable to the ones reported in Figure 5. We run CCF and RaBVItG using different mesh sizes (in the case of CCF, M1=202,302,402,502) and meshfree points (M2=20,30,40,50) in the case of RaBVItG. In Figure 5, we represent the different value functions for Player 1 (Player 2 has the same result since it is a symmetric game). The fact that CCF requires a mesh for the state variables (and controls) clearly affects the CPU time required to achieve convergence. Indeed, in order to compare the solution of RaBVItG with CCF algorithms, once we obtain the value function per player (as a M2- dimensional array), we interpolate our solution into the mesh defined in CCF using, again, RBF methods. Both algorithms converge to the same solution as the mesh size increases. However, as shown in Table 4, CCF method needs, in general, less value iterations but is much computationally expensive than RaBVItG costly (around 10 times greater).Table 4
Numerical Test 4 (non Scalar, Two Players): performance comparison between RaBVitG and CCF: mesh sizeM, CPU time (s) and number of iterations (Its).
CCF RaBVitG size 2 0 2 3 0 2 4 0 2 5 0 2 20 30 40 50  CPU time 2140 5140.3 9623.6 12854 479.11 603.08 734.8 1272.8  Its 333 287 266 248 542 457 397 275Figure 5
Numerical Test 4 (nonscalar, two players): value functions for player 1 returned by RaBVItG and CCF.


## 5. Conclusions

In this work, we have developed a novel numerical algorithm (RaBVItG) for solving multiplayer feedback-Nash differential games. RaBVItG is an algorithm based on radial basis function approximators (known to work properly without mesh), value iteration (a classical approach in dynamic programming), reinforcement learning (to solve a class of Hamilton-Jacobi-Bellman equations), and game iteration (to obtain, at each value iteration step, the corresponding feedback-Nash equilibrium).The general purpose of this algorithm is to deal with multiplayer problems, with two or more players, and allow the dimension of the control space for each player to be greater than one. For this reason, we have designed an algorithm which is meshfree for both state and control spaces. Additionally, we have also selected the semi-Lagrangian discretization for the HJB equation due to its fine results generally reported in the literature for problems similar to the differential games studied in this paper.We have validated RaBVItG by comparing the obtained results with the ones returned by another algorithm published in the literature (here, called CCF, see [1]). This particular algorithm is based on meshes in the state and control spaces and the way to get the Nash solution is relatively different from our approach. Indeed, CCF performs value iteration steps for a particular player fixing at every step the controls for the remaining players, until reaching convergence. In our case, RaBVItG alternates value iteration steps with a new step called game iteration. During this step, fixing the value for each player, the algorithm iterates until reaching convergence to a Nash equilibrium. So, we find the true (in the sense of convergence) Nash equilibrium for an approximated value (again, in the sense of convergence). As shown in our numerical experiments, our approach seems to exhibit a better accuracy when approximating the true equilibrium. Indeed, the different experiments performed in order to compare both algorithms tend to show that RaBVItG is, on average, 10 times faster than CCF and, in some of the experiments, with a smaller quantity of data points, RaBVitG obtains smaller errors.Additionally, RaBVItG can be implemented in real problems with more than two players and more than two controls per player, which are problems difficult to solve with CCF.Future lines of work on RaBVItG include the study of a stochastic version (e.g., including a multivariate diffusion term in the state equations) and different game iterations schemes in order to find other typical equilibria studied in the literature (e.g., Stackelberg, or Pareto). Furthermore, as discussed in Remark2, a study of the theoretical convergence properties of the proposed algorithm should be performed.



---

*Source: 1417275-2019-05-28.xml*