Model-based reinforcement learning in differential graphical games

TL;DR

This paper integrates differential game theory with actor-critic-identifier reinforcement learning to develop cooperative control strategies for multi-agent formation tracking with uncertain nonlinear dynamics, demonstrated through simulations.

Contribution

It introduces a novel model-based reinforcement learning approach combining differential game theory for formation control of multi-agent systems with uncertain dynamics.

Findings

Effective formation tracking achieved in simulations

Robustness to uncertain nonlinear dynamics demonstrated

Communication topology with spanning tree facilitates control

Abstract

This paper seeks to combine differential game theory with the actor-critic-identifier architecture to determine forward-in-time, approximate optimal controllers for formation tracking in multi-agent systems, where the agents have uncertain heterogeneous nonlinear dynamics. A continuous control strategy is proposed, using communication feedback from extended neighbors on a communication topology that has a spanning tree. A model-based reinforcement learning technique is developed to cooperatively control a group of agents to track a trajectory in a desired formation. Simulation results are presented to demonstrate the performance of the developed technique.

Equations

$${\dot{x}}_i=f_i\left(x_i\right)+g_i\left(x_i\right)u_i,$$

$$e_i=\sum _{j\in \left\{0\right\}\cup {\mathcal{N}}_{-i}}{a}_{ij}\left(\left(x_i-x_j\right)-x_{dij}\right).$$

$$e_i=\sum _{j\in \left\{0\right\}\cup {\mathcal{N}}_{-i}}{a}_{ij}\left(\left(x_i-x_{di0}\right)-\left(x_j-x_{dj0}\right)\right).$$

$$\mathcal{E}=\left(\left(\mathcal{L}+{\mathcal{A}}_0\right)\otimes I_n\right)\left(\mathcal{X}-{\mathcal{X}}_d-{\mathcal{X}}_0\right),$$

$$f_i\left(x_j+x_{dij}\right)+g_i\left(x_j+x_{dij}\right)u_{ij}={\dot{x}}_j.$$

$$\dot{x}=g\left(x\right)u,g\left(x\right)=\left[\begin{array}{cc}\cos \left(x_3\right)& 0\\ \sin \left(x_3\right)& 0\\ 0& 1\end{array}\right].$$

$${\mu }_i\triangleq \sum _{j\in {\mathcal{N}}_{-i}\cup \left\{0\right\}}{a}_{ij}\left(u_i-u_{ij}\right).$$

$${\mu }_{{\mathcal{S}}_i}={\mathcal{L}}_{gi}\left({\mathcal{E}}_i\right)u_{{\mathcal{S}}_i}-F_i\left({\mathcal{E}}_i\right),$$

$${\left[{\mathcal{L}}_{gi}\left({\mathcal{E}}_i\right)\right]}_{kl}=\{\begin{array}{ll}-a_{{\lambda }_i^k{\lambda }_i^l}{g}_{{\lambda }_i^k{\lambda }_i^l}\left(x_{{\lambda }_i^l}\right),& \forall l\ne k,\\ \text{sideset}{}^{{\lambda }_i^k}\sum a_{{\lambda }_i^{k}j}{I}_{m_{{\lambda }_i^k}},& \forall l=k,\end{array}$$

$$F_i\left({\mathcal{E}}_i\right)\triangleq {\left[\text{sideset}{}^i\sum a_{{\lambda }_i^{1}j}{f}_{{\lambda }_i^{1}j}^T\left(x_j\right),\cdots ,\text{sideset}{}^{{\lambda }_i^{s_i}}\sum a_{{\lambda }_i^{s_i}j}{f}_{{\lambda }_i^{s_i}j}^T\left(x_j\right)\right]}^T.$$

$${\left[{\mathcal{L}}_{gi}\left({\mathcal{E}}_i\right)\right]}_{kl}=\{\begin{array}{ll}-a_{{\lambda }_i^k{\lambda }_i^l}{I}_2,& \forall l\ne k,\\ \text{sideset}{}^{{\lambda }_i^k}\sum a_{{\lambda }_i^{k}j}{I}_2,& \forall l=k,\end{array}$$

$$u_{{\mathcal{S}}_i}={\mathcal{L}}_{gi}^{-1}\left({\mathcal{E}}_i\right){\mu }_{{\mathcal{S}}_i}+{\mathcal{L}}_{gi}^{-1}\left({\mathcal{E}}_i\right)F_i\left({\mathcal{E}}_i\right).$$

$$u_i={\mathcal{L}}_{gi}^i\left({\mathcal{E}}_i\right){\mu }_{{\mathcal{S}}_i}+{\mathcal{L}}_{gi}^i\left({\mathcal{E}}_i\right)F_i\left({\mathcal{E}}_i\right),$$

$$u_j={\mathcal{L}}_{gi}^j\left({\mathcal{E}}_i\right){\mu }_{{\mathcal{S}}_i}+{\mathcal{L}}_{gi}^j\left({\mathcal{E}}_i\right)F_i\left({\mathcal{E}}_i\right).$$

$${\dot{e}}_i={\mathcal{F}}_i\left({\mathcal{E}}_i\right)+{\mathcal{G}}_i\left({\mathcal{E}}_i\right){\mu }_{{\mathcal{S}}_i},$$

$${\dot{x}}_i={\mathcal{F}}_i\left({\mathcal{E}}_i\right)+{\mathcal{G}}_i\left({\mathcal{E}}_i\right){\mu }_{{\mathcal{S}}_i},$$

$$\begin{gathered} {\mathcal{F}}_i\left({\mathcal{E}}_i\right)\triangleq \text{sideset}{}^i\sum a_{ij}{g}_i\left(x_i\right){\mathcal{L}}_{gi}^i\left({\mathcal{E}}_i\right)F_i\left({\mathcal{E}}_i\right)-\text{sideset}{}^i\sum a_{ij}{f}_j\left(x_j\right) \\ -\text{sideset}{}^i\sum a_{ij}{g}_j\left(x_j\right){\mathcal{L}}_{gi}^j\left({\mathcal{E}}_i\right)F_i\left({\mathcal{E}}_i\right)+\text{sideset}{}^i\sum a_{ij}{f}_i\left(x_i\right), \end{gathered}$$

$${\mathcal{G}}_i\left({\mathcal{E}}_i\right)\triangleq \text{sideset}{}^i\sum a_{ij}\left(g_i\left(x_i\right){\mathcal{L}}_{gi}^i\left({\mathcal{E}}_i\right)-g_j\left(x_j\right){\mathcal{L}}_{gi}^j\left({\mathcal{E}}_i\right)\right),$$

$${\mathcal{F}}_i\left({\mathcal{E}}_i\right)\triangleq f_i\left(x_i\right)+g_i\left(x_i\right){\mathcal{L}}_{gi}^i\left({\mathcal{E}}_i\right)F_i\left({\mathcal{E}}_i\right),$$

$$J_i\left(e_i\left(\cdot \right),{\mu }_i\left(\cdot \right)\right)\triangleq \underset{0}{\overset{\infty }{\int }}{r}_i\left(e_i\left(\sigma \right),{\mu }_i\left(\sigma \right)\right)\text{d}\sigma$$

$$\begin{gathered} V_i^{{\overline{\mu }}_i,{\overline{\mu }}_{{\mathcal{S}}_{-i}}}\left({\mathcal{E}}_i\right)\triangleq \\ \underset{t}{\overset{\infty }{\int }}{r}_i\left(h_{ei}^{{\overline{\mu }}_i,{\overline{\mu }}_{{\mathcal{S}}_{-i}}}\left(\sigma ,t,{\mathcal{E}}_i\right),{\overline{\mu }}_i\left({\mathcal{H}}_i^{{\overline{\mu }}_i,{\overline{\mu }}_{{\mathcal{S}}_{-i}}}\left(\sigma ,t,{\mathcal{E}}_i\right)\right)\right)\text{d}\sigma , \end{gathered}$$

$$V_j^{\ast }\left({\mathcal{E}}_j\right)\triangleq V_j^{{\mu }_j^{\ast },{\mu }_{{\mathcal{S}}_{-j}}^{\ast }}\left({\mathcal{E}}_j\right)\le V_j^{{\overline{\mu }}_j,{\mu }_{{\mathcal{S}}_{-j}}^{\ast }}\left({\mathcal{E}}_j\right),$$

$$\sum_{j\in\mathcal{S}_{i}}\nabla_{e_{j}}V_{i}^{*}\left(\mathcal{E}_{i}\right)\left(\mathscr{F}_{j}\left(\mathcal{E}_{i}\right)+\mathscr{G}_{j}\left(\mathcal{E}_{i}\right)\mu_{\mathcal{S}_{j}}^{*}\left(\mathcal{E}_{i}\right)\right)\\ +\nabla_{x_{i}}V_{i}^{*}\left(\mathcal{E}_{i}\right)\left(\mathcal{F}_{i}\left(\mathcal{E}_{i}\right)+\mathcal{G}_{i}\left(\mathcal{E}_{i}\right)\mu_{\mathcal{S}_{i}}^{*}\left(\mathcal{E}_{i}\right)\right)\\ +\overline{Q}_{i}\left(\mathcal{E}_{i}\right)+\mu_{i}^{*T}\left(\mathcal{E}_{i}\right)R_{i}\mu_{i}^{*}\left(\mathcal{E}_{i}\right)=0,\>\forall\mbox{$\mathcal{E}$}_{i}\in\mathbb{R}^{n\left(s_{i}+1\right)},$$

$$\begin{gathered} V_i^{{\overline{\mu }}_i,{\mu }_{{\mathcal{S}}_{-i}}^{\ast }}\left({\left[e_i^T,{\mathcal{E}}_{-i}^T\right]}^T\right)= \\ \underset{t}{\overset{\infty }{\int }}{r}_i\left(h_{ei}^{{\overline{\mu }}_i,{\mu }_{{\mathcal{S}}_{-i}}^{\ast }}\left(\sigma ,t,{\mathcal{E}}_i\right),{\overline{\mu }}_i\left({\mathcal{H}}_i^{{\overline{\mu }}_i,{\mu }_{{\mathcal{S}}_{-i}}^{\ast }}\left(\sigma ,t,{\mathcal{E}}_i\right)\right)\right)\text{d}\sigma . \end{gathered}$$

$${\overline{V}}_i^{{\overline{\mu }}_i,{\mu }_{{\mathcal{S}}_{-i}}^{\ast }}\left(e_i,t\right)\triangleq V_i^{{\overline{\mu }}_i,{\mu }_{{\mathcal{S}}_{-i}}^{\ast }}\left({\left[e_i^T,{\mathcal{E}}_{-i}^T\left(t\right)\right]}^T\right),$$

$$\begin{gathered} \frac{\partial {\overline{V}}_i^{\ast }\left(e_i,t\right)}{\partial e_i}\left({\mathcal{F}}_i\left({\mathcal{E}}_i\right)+{\mathcal{G}}_i\left({\mathcal{E}}_i\right){\mu }_{{\mathcal{S}}_i}^{\ast }\left({\mathcal{E}}_i\right)\right)+\frac{\partial {\overline{V}}_i^{\ast }\left(e_i,t\right)}{\partial t} \\ +Q_i\left(e_i\right)+{\mu }_i^{\ast T}\left({\mathcal{E}}_i\right)R_i{\mu }_i^{\ast }\left({\mathcal{E}}_i\right)=0. \end{gathered}$$

$$\frac{\partial {\overline{V}}_i^{\ast }\left(e_i,t\right)}{\partial e_i}=\frac{\partial V_i^{\ast }\left({\mathcal{E}}_i\right)}{\partial e_i},$$

$$\begin{gathered} \frac{\partial {\overline{V}}_i^{\ast }\left(e_i,t\right)}{\partial t}=\frac{\partial V_i^{\ast }\left({\mathcal{E}}_i\right)}{\partial x_i}\left({\mathcal{F}}_i\left({\mathcal{E}}_i\right)+{\mathcal{G}}_i\left({\mathcal{E}}_i\right){\mu }_{{\mathcal{S}}_i}^{\ast }\left({\mathcal{E}}_i\right)\right) \\ +\sum _{j\in {\mathcal{S}}_{-i}}\frac{\partial V_i^{\ast }\left({\mathcal{E}}_i\right)}{\partial e_j}\left({\mathcal{F}}_j\left({\mathcal{E}}_i\right)+{\mathcal{G}}_j\left({\mathcal{E}}_i\right){\mu }_{{\mathcal{S}}_j}^{\ast }\left({\mathcal{E}}_i\right)\right), \end{gathered}$$

$$\begin{gathered} {\mu }_i^{\ast }\left({\mathcal{E}}_i\right)=-\frac{1}{2}{R}_i^{-1}\sum _{j\in {\mathcal{S}}_i}{\left({\mathcal{G}}_j^i\left({\mathcal{E}}_i\right)\right)}^T{\left({\nabla }_{e_j}{V}_i^{\ast }\left({\mathcal{E}}_i\right)\right)}^T \\ -\frac{1}{2}{R}_i^{-1}{\left({\mathcal{G}}_i^i\left({\mathcal{E}}_i\right)\right)}^T{\left({\nabla }_{x_i}{V}_i^{\ast }\left({\mathcal{E}}_i\right)\right)}^T, \end{gathered}$$

$$\begin{gathered} {\delta }_i\left({\mathcal{E}}_i,{\hat{W}}_{ci},{\left({\hat{W}}_a\right)}_{{\mathcal{S}}_i}\right)\triangleq {\hat{\mu }}_i^T\left({\mathcal{E}}_i,{\hat{W}}_{ai}\right)R{\hat{\mu }}_i\left({\mathcal{E}}_i,{\hat{W}}_{ai}\right) \\ +\sum _{j\in {\mathcal{S}}_i}{\nabla }_{e_j}{\hat{V}}_i\left({\mathcal{E}}_i,{\hat{W}}_{ci}\right){\mathcal{G}}_j\left({\mathcal{E}}_j\right){\hat{\mu }}_{{\mathcal{S}}_j}\left({\mathcal{E}}_j,{\left({\hat{W}}_a\right)}_{{\mathcal{S}}_j}\right) \\ +{\nabla }_{x_i}{\hat{V}}_i\left({\mathcal{E}}_i,{\hat{W}}_{ci}\right)\left({\mathcal{F}}_i\left({\mathcal{E}}_i\right)+{\mathcal{G}}_i\left({\mathcal{E}}_i\right){\hat{\mu }}_{{\mathcal{S}}_i}\left({\mathcal{E}}_i,{\left({\hat{W}}_a\right)}_{{\mathcal{S}}_i}\right)\right) \\ +\sum _{j\in {\mathcal{S}}_i}{\nabla }_{e_j}{\hat{V}}_i\left({\mathcal{E}}_i,{\hat{W}}_{ci}\right){\mathcal{F}}_j\left({\mathcal{E}}_j\right)+Q_i\left(e_i\right). \end{gathered}$$

Full text

I Introduction

In the past few decades, reinforcement learning (RL)-based techniques have been established as primary tools for online real-time optimization [1, 2, 3, 4, 5, 6, 7]. RL techniques are valuable not only for optimization but also for control synthesis in complex systems such as a distributed network of cognitive agents. Combined efforts from multiple autonomous agents can yield tactical advantages including: improved munitions effects; distributed sensing, detection, and threat response; and distributed communication pipelines [8, 9]. While coordinating behaviors among autonomous agents is a challenging problem that has received mainstream focus, unique challenges arise when seeking optimal autonomous collaborative behaviors. For example, most collaborative control literature focuses on centralized approaches that require all nodes to continuously communicate with a central agent, yielding a heavy communication demand that is subject to failure due to delays, and missing information[10]. Furthermore, the central agent is required to carry enough on-board computational resources to process the data and to generate command signals. These challenges motivate the need to minimize communication for guidance, navigation and control tasks, and to distribute the computational burden among the agents.

Since all the agents in a network have independent collaborative or competitive objectives, the resulting optimization problem is a multi-objective optimization problem. Differential game theory is often used to define optimality in multi-objective optimization problems [11, 12, 13, 14, 15, 16]. For example, a Nash equilibrium solution to a multi-objective optimization problem is said to be achieved if none of the players can benefit from a unilateral deviation from the equilibrium [17]. Thus, Nash equilibrium solutions provide a secure set of strategies in the sense that none of the players have an incentive to diverge from their equilibrium policy. Hence, Nash equilibrium has been a widely used solution concept in differential game-based control techniques. Online real-time solutions to differential games with centralized objectives are presented in results such as [18, 19, 20, 21, 22]; however, since these results solve problems with centralized objectives (i.e., each agent minimizes or maximizes a cost function that penalizes the states of all the agents in the network), they are not applicable for a network of agents with independent decentralized objectives (i.e., each agent minimizes or maximizes a cost function that penalizes only the error states corresponding to itself).

In this paper, the objective is to obtain an online forward-in-time feedback-Nash equilibrium solution (cf. [23, 24, 25, 26, 27, 28]) to an infinite-horizon formation tracking problem, where each agent desires to follow a mobile leader while the group maintains a desired formation. The agents try to minimize cost functions that penalize their own formation tracking errors and their own control efforts.

Various methods have been developed to solve optimal tracking problems for linear systems. In [29, 30, 31, 32], optimal controllers are developed to cooperatively control agents with linear dynamics. In [33], a differential game-based approach is developed for unmanned aerial vehicles to achieve distributed Nash strategies. In [34], an optimal consensus algorithm is developed for a cooperative team of agents with linear dynamics using only partial information.

For nonlinear systems, a MPC-based approach is presented in [35]; however, no stability or convergence analysis is presented. A stable distributed MPC-based approach is presented in [36] for nonlinear discrete-time systems with known nominal dynamics. Asymptotic stability is proved without any interaction between the nodes; however, a nonlinear optimal control problem needs to be solved at every iteration to implement the controller. An optimal tracking approach for formation control is presented in [37] using single network adaptive critics where the value function is learned offline. Recently, a leader-based consensus algorithm is developed in [38] where exact model of the system dynamics is utilized, and convergence to optimality is obtained under a persistence of excitation condition.

For multi-agent problems with decentralized objectives, the desired action by an individual agent depends on the actions and the resulting trajectories of its neighbors; hence, the error system for each agent is a complex nonautonomous dynamical system. Nonautonomous systems, in general, have non-stationary value functions. Since non-stationary functions are difficult to approximate using parameterized function approximation schemes such as neural networks (NNs), designing optimal policies for nonautonomous systems is challenging.

Since the external influence from neighbors renders the dynamics of each agent nonautonomous, optimization in a network of agents presents challenges similar to optimal tracking problems. Using insights gained from the authors’ previous work on optimal tracking problems [39], this paper develops a model-based RL technique to generate feedback-Nash equilibrium policies online, for agents in a network with cooperative or competitive objectives. In particular, the network of agents is separated into autonomous subgraphs, and the differential game is solved separately on each subgraph.

The primary contribution of this paper is the formulation and online approximate feedback-Nash equilibrium solution of an optimal network formation tracking problem. A relative control error minimization technique is introduced to facilitate the formulation of a feasible infinite-horizon total-cost differential graphical game. Dynamic programming-based feedback-Nash equilibrium solution of the differential graphical game is facilitated via the development of a set of coupled Hamilton-Jacobi (HJ) equations. The developed approximate feedback-Nash equilibrium solution is analyzed using a Lyapunov-based stability analysis to demonstrate ultimately bounded formation tracking in the presence of uncertainties.

II Notation

Throughout the paper, $\mathbb{R}^{n}$ denotes $n-$dimensional Euclidean space, $\mathbb{R}_{>a}$ denotes the set of real numbers strictly greater than $a\in\mathbb{R}$, and $\mathbb{R}_{\geq a}$ denotes the set of real numbers greater than or equal to $a\in\mathbb{R}$. Unless otherwise specified, the domain of all the functions is assumed to be $\mathbb{R}_{\geq 0}$. Functions with domain $\mathbb{R}_{\geq 0}$ are defined by abuse of notation using only their image. For example, the function $x:\mathbb{R}_{\geq 0}\to\mathbb{R}^{n}$ is defined by abuse of notation as $x\in\mathbb{R}^{n}$. By abuse of notation, the state variables are also used to denote state trajectories. For example, the state variable $x$ in the equation $\dot{x}=f\left(x\right)+u$ is also used as $x\left(t\right)$ to denote the state trajectory, i.e., the general solution $x:\mathbb{R}_{\geq 0}\to\mathbb{R}^{n}$ to $\dot{x}=f\left(x\right)+u$ evaluated at time $t$. Unless otherwise specified, all the mathematical quantities are assumed to be time-varying. Unless otherwise specified, an equation of the form $g\left(x\right)=f+h\left(y,t\right)$ is interpreted as $g\left(x\left(t\right)\right)=f\left(t\right)+h\left(y\left(t\right),t\right)$ for all $t\in\mathbb{R}_{\geq 0}$, and a definition of the form $g\left(x,y\right)\triangleq f\left(y\right)+h\left(x\right)$ for functions $g:A\times B\to C$, $f:B\to C$ and $h:A\to C$ is interpreted as $g\left(x,y\right)\triangleq f\left(y\right)+h\left(x\right),\>\forall\left(x,y\right)\in A\times B$. The total derivative $\frac{\partial f\left(x\right)}{\partial x}$ is denoted by $\nabla f$ and the partial derivative $\frac{\partial f\left(x,y\right)}{\partial x}$ is denoted by $\nabla_{x}f\left(x,y\right)$. An $n\times n$ identity matrix is denoted by $I_{n}$, $n\times m$ matrices of zeros and ones are denoted by $\mathbf{0}_{n\times m}$ and $\mathbf{1}_{n\times m}$, respectively, and $\mathbf{1}_{S}$ denotes the indicator function of the set $S$.

III Graph Theory Preliminaries

Consider a set of $N$ autonomous agents moving in the state space $\mathbb{R}^{n}$. The control objective is for the agents to maintain a desired formation with respect to a leader. The state of the leader is denoted by $x_{0}\in\mathbb{R}^{n}$. The agents are assumed to be on a network with a fixed communication topology modeled as a static directed graph (i.e. digraph).

Each agent forms a node in the digraph. The set of all nodes excluding the leader is denoted by $\mathcal{N}=\left\{1,\cdots N\right\}$ and the leader is denoted by node 0. If node $i$ can receive information from node $j$ then there exists a directed edge from the $j\textsuperscript{th}$ to the $i$th node of the digraph, denoted by the ordered pair $\left(j,i\right)$. Let $E$ denote the set of all edges. Let there be a positive weight $a_{ij}\in\mathbb{R}$ associated with each edge $\left(j,i\right)$. Note that $a_{ij}\neq 0$ if and only if $\left(j,i\right)\in E.$ The digraph is assumed to have no repeated edges, i.e., $\left(i,i\right)\notin E,\forall i$, which implies $a_{ii}=0,\forall i$. The neighborhood sets of node $i$ are denoted by $\mathcal{N}_{-i}$ and $\mathcal{N}_{i}$, defined as $\mathcal{N}_{-i}\triangleq\left\{j\in\mathcal{N}\mid\left(j,i\right)\in E\right\}$ and $\mathcal{N}_{i}\triangleq\mathcal{N}_{-i}\cup\left\{i\right\}$.

To streamline the analysis, an adjacency matrix $\mathcal{A}\in\mathbb{R}^{N\times N}$ is defined as $\mathcal{A}\triangleq\left[a_{ij}\mid i,j\in\mathcal{N}\right]$, a diagonal pinning gain matrix $\mathcal{A}_{0}\in\mathbb{R}^{N\times N}$ is defined as $\mathcal{A}_{0}\triangleq\mbox{diag}\left(\left[a_{10},\cdots,a_{N0}\right]\right)$, an in-degree matrix $\mathcal{D}\in\mathbb{R}^{N\times N}$ is defined as $\mathcal{D}\triangleq\mbox{diag}\left(d_{i}\right),$ where $d_{i}\triangleq\sum_{j\in\mathcal{N}_{i}}a_{ij}$, and a graph Laplacian matrix $\mathcal{L}\in\mathbb{R}^{N\times N}$ is defined as $\mathcal{L}\triangleq\mathcal{D}-\mathcal{A}$. The graph is assumed to have a spanning tree, i.e., given any node $i$, there exists a directed path from the leader [math] to node $i$. A node $j$ is said to be an extended neighbor of node $i$ if there exists a directed path from node $j$ to node $i$. The extended neighborhood set of node $i$, denoted by $\mathcal{S}_{-i}$, is defined as the set of all extended neighbors of node $i.$ Formally, $\mathcal{S}_{-i}\triangleq\{j\in\mathcal{N}\mid j\neq i\wedge\exists\kappa\leq N,\>\left\{j_{1},\cdots j_{\kappa}\right\}\subset\mathcal{N}\mid\left\{\left(j,j_{1}\right),\left(j_{1},j_{2}\right),\cdots,\left(j_{\kappa},i\right)\right\}\subset 2^{E}\}$. Let $\mathcal{S}_{i}\triangleq\mathcal{S}_{-i}\cup\left\{i\right\}$, and let the edge weights be normalized such that $\sum_{j}a_{ij}=1$ for all $i\in\mathcal{N}$. Note that the sub-graphs are nested in the sense that $\mathcal{S}_{j}\subseteq\mathcal{S}_{i}$ for all $j\in\mathcal{S}_{i}$.

IV Problem Formulation

The state $x_{i}\in\mathbb{R}^{n}$ of each agent evolves according to the control affine dynamics

[TABLE]

where $u_{i}\in\mathbb{R}^{m_{i}}$ denotes the control input, and $f_{i}:\mathbb{R}^{n}\to\mathbb{R}^{n}$ and $g_{i}:\mathbb{R}^{n}\to\mathbb{R}^{n\times m_{i}}$ are locally Lipschitz continuous functions.

Assumption 1**.**

The dynamics of the leader are described by $\dot{x}_{0}=f_{0}\left(x_{0}\right),$ where $f_{0}:\mathbb{R}^{n}\to\mathbb{R}^{n}$ is a locally Lipschitz continuous function. The function $f_{0}$, and the initial condition $x_{0}\left(t_{0}\right)$ are selected such that the trajectory $x_{0}\left(t\right)$ is uniformly bounded for all $t\in\mathbb{R}_{\geq t_{0}}$.

The control objective is for the agents to maintain a predetermined formation (with respect to an inertial reference frame) around the leader while minimizing their own cost functions. For all $i\in\mathcal{N}$, the $i\textsuperscript{th}$ agent is aware of its constant desired relative position $x_{dij}\in\mathbb{R}^{n}$ with respect to all its neighbors $j\in\mathcal{N}_{-i}$, such that the desired formation is realized when $x_{i}-x_{j}\to x_{dij}$ for all $i,j\in\mathcal{N}$.111The vectors $x_{dij}$ are assumed to be fixed in an inertial reference frame, i.e., the final desired formation is rigid and its motion in an inertial reference frame can be described as pure translation. To facilitate the control design, the formation is expressed in terms of a set of constant vectors $\left\{x_{di0}\in\mathbb{R}^{n}\right\}_{i\in\mathcal{N}}$ where each $x_{di0}$ denotes the constant final desired position of agent $i$ with respect to the leader. The vectors $\left\{x_{di0}\right\}_{i\in\mathcal{N}}$ are unknown to the agents not connected to the leader, and the known desired inter agent relative position can be expressed in terms of $\left\{x_{di0}\right\}_{i\in\mathcal{N}}$ as $x_{dij}=x_{di0}-x_{dj0}$. The control objective is thus satisfied when $x_{i}\to x_{di0}+x_{0}$ for all $i\in\mathcal{N}$. To quantify the objective, local neighborhood tracking error signals are defined as

[TABLE]

To facilitate the analysis, the error signals in (2) are expressed in terms of the unknown leader-relative desired positions as

[TABLE]

Stacking the error signals in a vector \mathcal{E}\triangleq\left[\begin{array}[t]{cccc}e_{1}^{T},&e_{2}^{T},&\cdots,&e_{N}^{T}\end{array}\right]^{T}\in\mathbb{R}^{nN} the equation in (3) can be expressed in a matrix form

[TABLE]

where $\mathcal{X}=$ $\left[x_{1}^{T},\right.$ $x_{2}^{T},$ $\cdots,$ \left.x_{N}^{T}\right]^{T}$$\in\mathbb{R}^{nN}, $\mathcal{X}_{d}=$ $\left[x_{d10}^{T},\right.$ $x_{d20}^{T},$ $\cdots,$ \left.x_{dN0}^{T}\right]^{T}$$\in\mathbb{R}^{nN}, $\mathcal{X}_{0}=$ $\left[x_{0}^{T},\right.$ $x_{0}^{T},$ $\cdots,$ \left.x_{0}^{T}\right]^{T}$$\in\mathbb{R}^{nN}, and $\otimes$ denotes the Kronecker product. Using (4), it can be concluded that provided the matrix $\left(\left(\mathcal{L}+\mathcal{A}_{0}\right)\otimes I_{n}\right)\in\mathbb{R}^{nN\times nN}$ is nonsingular, $\left\|\mathcal{E}\right\|\to 0$ implies $x_{i}\to x_{di0}+x_{0}$ for all $i\in\mathcal{N}$, and hence, the satisfaction of control objective. The matrix $\left(\left(\mathcal{L}+\mathcal{A}_{0}\right)\otimes I_{n}\right)$ is nonsingular provided the graph has a spanning tree with the leader at the root [40]. To facilitate the formulation of an optimization problem, the following section explores the functional dependence of the state-value functions for the network of agents.

IV-A *Elements of the value

function*

The dynamics for the open-loop neighborhood tracking error are \dot{e}_{i}=\sum_{j\in\left\{0\right\}\cup\mathcal{N}_{-i}}a_{ij}\Bigl{(}f_{i}\left(x_{i}\right)+g_{i}\left(x_{i}\right)u_{i}-f_{j}\left(x_{j}\right)-g_{j}\left(x_{j}\right)u_{j}\Bigr{)}. Under the temporary assumption that each controller $u_{i}$ is an error-feedback controller, i.e. $u_{i}\left(t\right)=\hat{u}_{i}\left(e_{i}\left(t\right),t\right)$, the error dynamics are expressed as \dot{e}_{i}=\sum_{j\in\left\{0\right\}\cup\mathcal{N}_{-i}}a_{ij}\Bigl{(}f_{i}\left(x_{i}\right)+g_{i}\left(x_{i}\right)\hat{u}_{i}\left(e_{i},t\right)-f_{j}\left(x_{j}\right)-g_{j}\left(x_{j}\right)\hat{u}_{j}\left(e_{j},t\right)\Bigr{)}. Thus, the error trajectory $\left\{e_{i}\left(t\right)\right\}_{t=t_{0}}^{\infty}$ , where $t_{0}$ denotes the initial time, depends on $\hat{u}_{j}\left(e_{j}\left(t\right),t\right)$, $\forall j\in\mathcal{N}_{i}$. Similarly, the error trajectory $\left\{e_{j}\left(t\right)\right\}_{t=t_{0}}^{\infty}$ depends on $\hat{u}_{k}\left(e_{k}\left(t\right),t\right),\forall k\in\mathcal{N}_{j}$. Recursively, the trajectory $\left\{e_{i}\left(t\right)\right\}_{t=t_{0}}^{\infty}$ depends on $\hat{u}_{j}\left(e_{j}\left(t\right),t\right)$, and hence, on $e_{j}\left(t\right),\forall j\in\mathcal{S}_{i}$. Thus, even if the controller for each agent is restricted to use local error feedback, the resulting error trajectories are interdependent. In particular, a change in the initial condition of one agent in the extended neighborhood causes a change in the error trajectories corresponding to all the extended neighbors. Consequently, the value function corresponding to an infinite-horizon optimal control problem where each agent tries to minimize $\intop_{t_{0}}^{\infty}\left(Q\left(e_{i}\left(\tau\right)\right)+R\left(u_{i}\left(\tau\right)\right)\right)d\tau$, where $Q:\mathbb{R}^{n}\to\mathbb{R}$ and $R:\mathbb{R}^{m_{i}}\to\mathbb{R}$ are positive definite functions, is dependent on the error states of all the extended neighbors.

Since the steady-state controllers required for formation tracking are generally nonzero, quadratic total-cost optimal control problems result in infinite costs, and hence, are infeasible. In the following section, relative steady-state controllers are derived to facilitate the formulation of a feasible optimal control problem.

IV-B Optimal formation tracking problem

When the agents are perfectly tracking the desired trajectory in the desired formation, even though the states of all the agents are different, the time-derivatives of the states of all the agents are identical. Hence, in steady state, the control signal applied by each agent must be such that the time derivatives of the states corresponding to the set of extended neighbors are identical. In particular, the relative control signal $u_{ij}\in\mathbb{R}^{m_{i}}$ that will keep node $i$ in its desired relative position with respect to node $j\in\mathcal{S}_{-i}$, i.e., $x_{i}=x_{j}+x_{dij}$, must be such that the time derivative of $x_{i}$ is the same as the time derivative of $x_{j}$. Using the dynamics of the agents from (1), and substituting the desired relative positions $x_{j}+x_{dij}$ for the states $x_{i}$, the relative control signals $u_{ij}$ must satisfy

[TABLE]

The relative steady-state control signals can be expressed in an explicit form provided the following assumption is satisfied.

Assumption 2**.**

The matrix $g_{i}\left(x\right)$ is full rank for all $i\in\mathcal{N}$ and for all $x\in\mathbb{R}^{n}$; furthermore, the relative steady-state control signal expressed as $u_{ij}=f_{ij}\left(x_{j}\right)+g_{ij}\left(x_{j}\right)u_{j},$ satisfies (5) along the desired trajectory, where $f_{ij}\left(x_{j}\right)\triangleq g_{i}^{+}\left(x_{j}+x_{dij}\right)\left(f_{j}\left(x_{j}\right)-f_{i}\left(x_{j}+x_{dij}\right)\right)\in\mathbb{R}^{m_{i}}$, $g_{ij}\left(x_{j}\right)\triangleq g_{i}^{+}\left(x_{j}+x_{dij}\right)g_{j}\left(x_{j}\right)\in\mathbb{R}^{m_{i}\times m_{j}}$, $g_{0}\left(x\right)\triangleq 0$ for all $x\in\mathbb{R}^{n}$, $u_{i0}\equiv 0$ for all $i\in\mathcal{N}$, and $g_{i}^{+}\left(x\right)$ denotes a pseudoinverse of the matrix $g_{i}\left(x\right)$ for all $x\in\mathbb{R}^{n}$ and for all $i\in\mathcal{N}$.

Assumption 2 places restrictions on the control-effectiveness matrices. The matrices $g_{i}\left(x\right)$ are full rank for a large class of systems including, but not limited to, kinematic wheels and fully actuated Euler-Lagrange systems with invertible inertia matrices. The second part of Assumption 2 requires the existence of a feedback controller that can keep the system on the desired trajectory if the system starts on the desired trajectory. This assumption depends on the systems, the network, the desired formation, and the desired trajectory; hence, insights into its satisfaction are hard to obtain in general. The satisfaction of this assumption needs to be verified on a case-by-case basis. For example, consider a kinematic wheel modeled as

[TABLE]

In this case, provided the formation satisfies $x_{dij}\left(3\right)=0$, that is, the target formation is such that all the kinematic wheels have the same steering angle, the functions $f_{ij}$ and $g_{ij}$ can be computed as $f_{ij}=0,$ and $g_{ij}=I_{2}$. The relative steady-state control is then $u_{ij}=u_{j}$, which satisfies $g\left(x_{j}+x_{dij}\right)u_{j}=\dot{x}_{j}$, and hence, Assumption 2 holds.

To facilitate the formulation of an optimal formation tracking problem, define the control errors $\mu_{i}\in\mathbb{R}^{m_{i}}$ as

[TABLE]

The control errors $\left\{\mu_{i}\right\}$ are treated as the design variables in the remainder of this paper. Since the control errors $\left\{\mu_{i}\right\}$ are designed and the controllers $\left\{u_{i}\right\}$ are implemented in practice, it is essential to invert the relationship in (7). To facilitate the inversion, let $\mathcal{S}_{i}^{o}\triangleq\left\{1,\cdots,s_{i}\right\}$, where $s_{i}\triangleq\left|\mathcal{S}_{i}\right|$. Let $\lambda_{i}:\mathcal{S}_{i}^{o}\to\mathcal{S}_{i}$ be a bijective map such that $\lambda_{i}\left(1\right)=i$. For notational brevity, let $\left(\cdot\right)_{\mathcal{S}_{i}}$ denote the concatenated vector $\left[\left(\cdot\right)_{\lambda_{i}^{1}}^{T},\left(\cdot\right)_{\lambda_{i}^{2}}^{T},\cdots,\left(\cdot\right)_{\lambda_{i}^{s_{i}}}^{T}\right]^{T}$, let $\left(\cdot\right)_{\mathcal{S}_{-i}}$ denote the concatenated vector $\left[\left(\cdot\right)_{\lambda_{i}^{2}}^{T},\cdots,\left(\cdot\right)_{\lambda_{i}^{s_{i}}}^{T}\right]^{T}$, let $\sum^{i}$ denote $\sum_{j\in\mathcal{N}_{-i}\cup\left\{0\right\}}$, let $\lambda_{i}^{j}$ denote $\lambda_{i}\left(j\right)$, let $\mathcal{E}_{i}\triangleq\left[e_{\mathcal{S}_{i}}^{T},x_{\lambda_{i}^{1}}^{T}\right]^{T}\in\mathbb{R}^{n\left(s_{i}+1\right)}$, and let $\mathcal{E}_{-i}\triangleq\left[e_{\mathcal{S}_{-i}}^{T},x_{\lambda_{i}^{1}}^{T}\right]^{T}\in\mathbb{R}^{ns_{i}}$. Then, the control error vectors $\mu_{\mathcal{S}_{i}}\in\mathbb{R}^{\sum_{k\in\mathcal{S}_{i}}m_{k}}$ can be expressed as

[TABLE]

where the matrices $\mathscr{L}_{gi}:\mathbb{R}^{n\left(s_{i}+1\right)}\to\mathbb{R}^{\sum_{k\in\mathcal{S}_{i}}m_{k}\times\sum_{k\in\mathcal{S}_{i}}m_{k}}$ are defined by

[TABLE]

where $k,l=1,2,\cdots,s_{i}$, and $F_{i}:\mathbb{R}^{n\left(s_{i}+1\right)}\to\mathbb{R}^{\sum_{k\in\mathcal{S}_{i}}m_{k}}$ are defined as

[TABLE]

Assumption 3**.**

The matrix $\mathscr{L}_{gi}\left(\mathcal{E}_{i}\left(t\right)\right)$ is invertible for all $t\in\mathbb{R}$ and for all $i\in\mathcal{N}$.

Assumption 3 is a controllability-like condition. Intuitively, Assumption 3 requires the control effectiveness matrices to be compatible to ensure the existence of relative control inputs that allow the agents to follow the desired trajectory in the desired formation. Assumption 3 depends on the systems, the network, the desired formation, and the desired trajectory; hence, insights into its satisfaction are hard to obtain in general. The satisfaction of this assumption needs to be verified on a case-by-case basis. For example, consider the kinematic wheel in (6). Provided the formation satisfies $x_{dij}\left(3\right)=0$, that is, the target formation is such that all the kinematic wheels have the same steering angle, we have $g_{ij}=I_{2},$ and hence, the matrices $\mathscr{L}_{gi}$ are given by

[TABLE]

It can be shown that $\mathscr{L}_{gi}=\mathcal{L}_{\mathcal{S}_{i}}\otimes I_{2}$, where $\mathcal{L}_{\mathcal{S}_{i}}$ denotes the Laplacian matrix corresponding to the subgraph $\mathcal{S}_{i}$. Hence, the graph connectivity condition ensures that the matrices $\mathscr{L}_{gi}$ are invertible, and in this specific case, Assumption 3 holds.

Using Assumption 3, the control vectors can be expressed as

[TABLE]

Let $\mathscr{L}_{gi}^{k}$ denote the $\left(\lambda_{i}^{-1}\left(k\right)\right)$th block row of $\mathscr{L}_{gi}^{-1}$. Then, the controllers $u_{i}$ can be implemented as

[TABLE]

and for any $j\in\mathcal{N}_{-i}$,

[TABLE]

Using (10) and (11), the error and the state dynamics for the agents can be represented as

[TABLE]

and

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

and $\mathcal{G}_{i}\left(\mathcal{E}_{i}\right)\triangleq g_{i}\left(x_{i}\right)\mathscr{L}_{gi}^{i}\left(\mathcal{E}_{i}\right)$.

Let $h_{ei}^{\overline{\mu}_{i},\overline{\mu}_{\mathcal{S}_{-i}}}\left(t,t_{0},\mathcal{E}_{i0}\right)$ and $h_{xi}^{\overline{\mu}_{i},\overline{\mu}_{\mathcal{S}_{-i}}}\left(t,t_{0},\mathcal{E}_{i0}\right)$ denote the trajectories of (12) and (13), respectively, with the initial time $t_{0}$, initial condition $\mathcal{E}_{i}\left(t_{0}\right)=\mathcal{E}_{i0}$, and policies $\overline{\mu}_{j}:\mathbb{R}^{n\left(s_{i}+1\right)}\to\mathbb{R}^{m_{i}},\>j\in\mathcal{S}_{i}$, and let $\mathcal{H}_{i}\triangleq\left[\left(h_{e}\right)_{\mathcal{S}_{i}}^{T},h_{x\lambda_{i}^{1}}^{T}\right]^{T}$. Define the cost functionals

[TABLE]

where $r_{i}:\mathbb{R}^{n}\times\mathbb{R}^{m_{i}}\to\mathbb{R}_{\geq 0}$ denote the local costs defined as $r_{i}\left(e_{i},\mu_{i}\right)\triangleq Q_{i}\left(e_{i}\right)+\mu_{i}^{T}R_{i}\mu_{i},$ where $Q_{i}:\mathbb{R}^{n}\to\mathbb{R}_{\geq 0}$ are positive definite functions, and $R_{i}\in\mathbb{R}^{m_{i}\times m_{i}}$ are constant positive definite matrices. The objective of each agent is to minimize the cost functional in (14). To facilitate the definition of a feedback-Nash equilibrium solution, define the value functions $V_{i}:\mathbb{R}^{n\left(s_{i}+1\right)}\to\mathbb{R}_{\geq 0}$ as

[TABLE]

where $V_{i}^{\overline{\mu}_{i},\overline{\mu}_{\mathcal{S}_{-i}}}\left(\mathcal{E}_{i}\right)$ denotes the total cost-to-go for Agent $i$ under the policies $\overline{\mu}_{\mathcal{S}_{i}}$, when the sub-graph $\mathcal{S}_{i}$ starts from the state $\mathcal{E}_{i}$. Note that the value functions in (15) are time-invariant because the dynamical systems $\left\{\dot{e}_{j}=\mathscr{F}_{j}\left(\mathcal{E}_{i}\right)+\mathscr{G}_{j}\left(\mathcal{E}_{i}\right)\mu_{\mathcal{S}_{j}}\right\}_{j\in\mathcal{S}_{i}}$ and $\dot{x}_{i}=\mathcal{F}_{i}\left(\mathcal{E}_{i}\right)+\mathcal{G}_{i}\left(\mathcal{E}_{i}\right)\mu_{\mathcal{S}_{i}}$ together form an autonomous dynamical system.

A graphical feedback-Nash equilibrium solution within the subgraph $\mathcal{S}_{i}$ is defined as the tuple of policies $\left\{\mu_{j}^{*}:\mathbb{R}^{n\left(s_{j}+1\right)}\to\mathbb{R}^{m_{j}}\right\}_{j\in\mathcal{S}_{i}}$ such that the value functions in (15) satisfy

[TABLE]

for all $j\in\mathcal{S}_{i}$, for all $\mathcal{E}_{i}\in\mathbb{R}^{n\left(s_{i}+1\right)}$ and for all admissible policies $\overline{\mu}_{j}$. Provided a feedback-Nash equilibrium solution exists and the value functions (15) are continuously differentiable for all $i\in\mathcal{N}$, the feedback-Nash equilibrium value functions can be characterized in terms of the following system of HJ equations:

[TABLE]

where $\overline{Q}_{i}:\mathbb{R}^{n\left(s_{i}+1\right)}\to\mathbb{R}$ is defined as $\overline{Q}_{i}\left(\mathcal{E}_{i}\right)\triangleq Q_{i}\left(e_{i}\right)$.

Theorem 1**.**

Provided a feedback-Nash equilibrium solution exists and that the value functions in (15) are continuously differentiable, the system of HJ equations in (16) constitutes a necessary and sufficient condition for $\left\{\mu_{j}^{*}:\mathbb{R}^{n\left(s_{j}+1\right)}\to\mathbb{R}^{m_{j}}\right\}_{j\in\mathcal{S}_{i}}$ to be a feedback-Nash equilibrium solution within the subgraph $\mathcal{S}_{i}$.

Proof:

Consider the cost functional in (14), and assume that all the extended neighbors of the $i\textsuperscript{th}$ agent follow their feedback-Nash equilibrium policies. The value function corresponding to any admissible policy $\overline{\mu}_{i}$ can be expressed as

[TABLE]

Treating the dependence on $\mathcal{E}_{-i}$ as explicit time dependence define

[TABLE]

for all $e_{i}\in\mathbb{R}^{n}$ and for all $t\in\mathbb{R}_{\geq 0}$. Assuming that the optimal controller that minimizes (14) when all the extended neighbors follow their feedback-Nash equilibrium policies exists, and that the optimal value function $\overline{V}_{i}^{*}\triangleq\overline{V}_{i}^{\mu_{i}^{*},\mu_{\mathcal{S}_{-i}}^{*}}$ exists and is continuously differentiable, optimal control theory for single objective optimization problems (cf. [41]) can be used to derive the following necessary and sufficient condition

[TABLE]

Using (17), the partial derivative with respect to the state can be expressed as

[TABLE]

for all $e_{i}\in\mathbb{R}^{n}$ and for all $t\in\mathbb{R}_{\geq 0}$, and the partial derivative with respect to time can be expressed as

[TABLE]

for all $e_{i}\in\mathbb{R}^{n}$ and for all $t\in\mathbb{R}_{\geq 0}$. Substituting (19) and (20) into (18) and repeating the process for each $i$, the system of HJ equations in (16) is obtained. ∎

Minimizing the HJ equations using the stationary condition, the feedback-Nash equilibrium solution is expressed in the explicit form

[TABLE]

for all $\mathcal{E}_{i}\in\mathbb{R}^{n\left(s_{i}+1\right)}$, where $\mathscr{G}_{j}^{i}\triangleq\mathscr{G}_{j}\frac{\partial\mu_{\mathcal{S}_{j}}^{*}}{\partial\mu_{i}^{*}}$, and $\mathcal{G}_{i}^{i}\triangleq\mathcal{G}_{i}\frac{\partial\mu_{\mathcal{S}_{i}}^{*}}{\partial\mu_{i}^{*}}$. Since an analytical solution of system of HJ equations in (16) is generally infeasible to obtain, the feedback-Nash value functions and the feedback-Nash policies are approximated using parametric approximation schemes $\hat{V}_{i}\left(\mathcal{E}_{i},\hat{W}_{ci}\right)$ and $\hat{\mu}_{i}\left(\mathcal{E}_{i},\hat{W}_{ai}\right)$, respectively, where $\hat{W}_{ci}\in\mathbb{R}^{L_{i}}$ and $\hat{W}_{ai}\in\mathbb{R}^{L_{i}}$ are parameter estimates. Substitution of the approximations $\hat{V}_{i}$ and $\hat{\mu}_{i}$ in (16) leads to a set of Bellman errors (BEs) $\delta_{i}$ defined as

[TABLE]

Approximation of the feedback-Nash equilibrium policies is realized by tuning the estimates $\hat{V}_{i}$ and $\hat{\mu}_{i}$ so as to minimize the BEs $\delta_{i}$. However, computation of $\delta_{i}$ in (22) and $u_{ij}$ in (7) requires exact model knowledge. In the following, a CL-based system identifier is developed to relax the exact model knowledge requirement and to facilitate the implementation of model-based RL via BE extrapolation (cf. [39]). In particular, the developed controllers do not require the knowledge of the system drift functions $f_{i}$.

V System Identification

On any compact set $\chi\subset\mathbb{R}^{n}$ the function $f_{i}$ can be represented using a NN as

[TABLE]

for all $x\in\mathbb{R}^{n}$, where $\theta_{i}\in\mathbb{R}^{P_{i}+1\times n}$ denote the unknown output-layer NN weights, $\sigma_{\theta i}:\mathbb{R}^{n}\to\mathbb{R}^{P_{i}+1}$ denotes a bounded NN basis function, $\epsilon_{\theta i}:\mathbb{R}^{n}\to\mathbb{R}^{n}$ denotes the function reconstruction error, and $P_{i}\in\mathbb{N}$ denotes the number of NN neurons. Using the universal function approximation property of single layer NNs, provided the rows of $\sigma_{\theta i}\left(x\right)$ form a proper basis, there exist constant ideal weights $\theta_{i}$ and positive constants $\overline{\theta_{i}}\in\mathbb{R}$ and $\overline{\epsilon_{\theta i}}\in\mathbb{R}$ such that $\left\|\theta_{i}\right\|_{F}\leq\overline{\theta_{i}}<\infty$ and $\sup_{x\in\chi}\left\|\epsilon_{\theta i}\left(x\right)\right\|\leq\overline{\epsilon_{\theta i}}$, where $\left\|\cdot\right\|_{F}$ denotes the Frobenius norm, i.e., $\left\|\theta\right\|_{F}\triangleq\sqrt{\mbox{tr}\left(\theta^{T}\theta\right)}$.

Assumption 4**.**

The bounds $\overline{\theta_{i}}$ and $\overline{\epsilon_{\theta i}}$ are known for all $i\in\mathcal{N}.$

Using an estimate $\hat{\theta}_{i}\in\mathbb{R}^{P_{i}+1\times n}$ of the weight matrix $\theta_{i},$ the function $f_{i}$ can be approximated by the function $\hat{f}_{i}:\mathbb{R}^{n}\times\mathbb{R}^{P_{i}+1\times n}\to\mathbb{R}^{n}$ defined by $\hat{f}_{i}\left(x,\hat{\theta}\right)\triangleq\hat{\theta}^{T}\sigma_{\theta i}\left(x\right).$ Based on (23), an estimator for online identification of the drift dynamics is developed as

[TABLE]

where $\tilde{x}_{i}\triangleq x_{i}-\hat{x}_{i}$, and $k_{i}\in\mathbb{R}$ is a positive constant learning gain. The following assumption facilitates concurrent learning (CL)-based system identification.

Assumption 5**.**

[42, 43] A history stack containing recorded state-action pairs $\left\{x_{i}^{k},u_{i}^{k}\right\}_{k=1}^{M_{\theta i}}$ along with numerically computed state derivatives $\left\{\dot{\bar{x}}_{i}^{k}\right\}_{k=1}^{M_{\theta i}}$ that satisfies

[TABLE]

is available a priori. In (25), $\sigma_{\theta i}^{k}\triangleq\sigma_{\theta i}\left(x_{i}^{k}\right)$, $\overline{d_{i}},\underline{\sigma_{\theta i}}\in\mathbb{R}$ are known positive constants, and $\lambda_{\min}\left(\cdot\right)$ denotes the minimum eigenvalue.

The weight estimates $\hat{\theta}_{i}$ are updated using the following CL-based update law:

[TABLE]

where $g_{i}^{k}\triangleq g_{i}\left(x_{i}^{k}\right)$, $k_{\theta i}\in\mathbb{R}$ is a constant positive CL gain, and $\Gamma_{\theta i}\in\mathbb{R}^{P_{i}+1\times P_{i}+1}$ is a constant, diagonal, and positive definite adaptation gain matrix.

To facilitate the subsequent stability analysis, a candidate Lyapunov function $V_{0i}:\mathbb{R}^{n}\times\mathbb{R}^{P_{i}+1\times n}\to\mathbb{R}$ is selected as

[TABLE]

where $\tilde{\theta}_{i}\triangleq\theta_{i}-\hat{\theta}_{i}$ and $\mbox{tr}\left(\cdot\right)$ denotes the trace of a matrix. Using (24)-(26), the identity $\mbox{tr}\left(\tilde{\theta}^{T}\left(\sum_{j=1}^{M_{\theta i}}\sigma_{\theta i}^{j}\sigma_{\theta i}^{j}\right)\tilde{\theta}\right)=\left(\mbox{vec}\left(\tilde{\theta}_{i}\right)\right)^{T}\left(\left(\sum_{j=1}^{M_{\theta i}}\sigma_{\theta i}^{j}\sigma_{\theta i}^{j}\right)\otimes I_{p+1}\right)\left(\mbox{vec}\left(\tilde{\theta}_{i}\right)\right)$, and the facts that $\lambda_{\min}\left\{\left(\left(\sum_{j=1}^{M_{\theta i}}\sigma_{\theta i}^{j}\sigma_{\theta i}^{j}\right)\otimes I_{p+1}\right)\right\}$ $=$ $\lambda_{\min}\left\{\sum_{j=1}^{M_{\theta i}}\sigma_{\theta i}^{j}\sigma_{\theta i}^{j}\right\}$ and $\lambda_{\max}\left\{\left(\left(\sum_{j=1}^{M_{\theta i}}\sigma_{\theta i}^{j}\sigma_{\theta i}^{j}\right)\otimes I_{p+1}\right)\right\}$ $=\lambda_{\max}\left\{\sum_{j=1}^{M_{\theta i}}\sigma_{\theta i}^{j}\sigma_{\theta i}^{j}\right\}$ (cf. [44, Theorem 4.2.12]), the following bound on the time derivative of $V_{0i}$ is established:

[TABLE]

where $\overline{d_{\theta i}}\triangleq\overline{d}_{i}\sum_{k=1}^{M_{\theta i}}\left\|\sigma_{\theta i}^{k}\right\|+\sum_{k=1}^{M_{\theta i}}\left(\left\|\epsilon_{\theta i}^{k}\right\|\left\|\sigma_{\theta i}^{k}\right\|\right)$. Using (27) and (28), a Lyapunov-based stability analysis can be used to show that $\hat{\theta}_{i}$ converges exponentially to a neighborhood around $\theta_{i}$.

VI Approximation of the BE and the relative steady-state controller

Using the approximations $\hat{f}_{i}$ for the functions $f_{i}$, the BEs in (22) can be approximated as

[TABLE]

In (29),

[TABLE]

The approximations $\hat{F}_{i}$, $\hat{\mathscr{F}}_{i}$, and $\hat{\mathcal{F}}_{i}$ are related to the original unknown functions as $\hat{F}_{i}\left(\mathcal{E}_{i},\theta_{\mathcal{S}_{i}}\right)+B_{i}\left(\mathcal{E}_{i}\right)=F_{i}\left(\mathcal{E}_{i}\right)$, $\hat{\mathscr{F}}_{i}\left(\mathcal{E}_{i},\theta_{\mathcal{S}_{i}}\right)+\mathscr{B}_{i}\left(\mathcal{E}_{i}\right)=\mathscr{F}_{i}\left(\mathcal{E}_{i}\right)$, and $\hat{\mathcal{F}}_{i}\left(\mathcal{E}_{i},\theta_{\mathcal{S}_{i}}\right)+\mathcal{B}_{i}\left(\mathcal{E}_{i}\right)=\mathcal{F}_{i}\left(\mathcal{E}_{i}\right)$, where $B_{i}$, $\mathscr{B}_{i}$, and $\mathcal{B}_{i}$ are $O\left(\left(\overline{\epsilon_{\theta}}\right)_{\mathcal{S}_{i}}\right)$ terms that denote bounded function approximation errors.

Using the approximations $\hat{f}_{i}$, an implementable form of the controllers in (9) is expressed as

[TABLE]

Using (8) and (30), an unmeasurable form of the virtual controllers implemented on the systems (12) and (13) is given by

[TABLE]

VII Value function approximation

On any compact set $\chi\in\mathbb{R}^{n\left(s_{i}+1\right)}$, the value functions can be represented as

[TABLE]

where $W_{i}\in\mathbb{R}^{L_{i}}$ are ideal NN weights, $\sigma_{i}:\mathbb{R}^{n\left(s_{i}+1\right)}\to\mathbb{R}^{L_{i}}$ are NN basis functions and $\epsilon_{i}:\mathbb{R}^{n\left(s_{i}+1\right)}\to\mathbb{R}$ are function approximation errors. Using the universal function approximation property of single layer NNs, provided $\sigma_{i}\left(\mathcal{E}_{i}\right)$ forms a proper basis, there exist constant ideal weights $W_{i}$ and positive constants $\overline{W_{i}}\in\mathbb{R}$ and $\overline{\epsilon_{i}},\overline{\nabla\epsilon_{i}}\in\mathbb{R}$ such that $\left\|W_{i}\right\|\leq\overline{W_{i}}<\infty$, $\sup_{\mathcal{E}_{i}\in\chi}\left\|\epsilon_{i}\left(\mathcal{E}_{i}\right)\right\|\leq\overline{\epsilon_{i}}$, and $\sup_{\mathcal{E}_{i}\in\chi}\left\|\nabla\epsilon_{i}\left(\mathcal{E}_{i}\right)\right\|\leq\overline{\nabla\epsilon_{i}}$.

Assumption 6**.**

The constants $\overline{\epsilon_{i}},$ $\overline{\nabla\epsilon_{i}},$ and $\overline{W_{i}}$ are known for all $i\in\mathcal{N}$.

Using (21) and (32), the feedback-Nash equilibrium policies are

[TABLE]

for all $\mathcal{E}_{i}\in\mathbb{R}^{n\left(s_{i}+1\right)},$ where $G_{\sigma i}\left(\mathcal{E}_{i}\right)\triangleq\sum_{j\in\mathcal{S}_{i}}\left(\mathscr{G}_{j}^{i}\left(\mathcal{E}_{i}\right)\right)^{T}\left(\nabla_{e_{j}}\sigma_{i}\left(\mathcal{E}_{i}\right)\right)^{T}+\left(\mathcal{G}_{i}^{i}\left(\mathcal{E}_{i}\right)\right)^{T}\left(\nabla_{x_{i}}\sigma_{i}\left(\mathcal{E}_{i}\right)\right)^{T}$ and $G_{\epsilon i}\left(\mathcal{E}_{i}\right)\triangleq\sum_{j\in\mathcal{S}_{i}}\left(\mathscr{G}_{j}^{i}\left(\mathcal{E}_{i}\right)\right)^{T}\left(\nabla_{e_{j}}\epsilon_{i}\left(\mathcal{E}_{i}\right)\right)^{T}+\left(\mathcal{G}_{i}^{i}\left(\mathcal{E}_{i}\right)\right)^{T}\left(\nabla_{x_{i}}\epsilon_{i}\left(\mathcal{E}_{i}\right)\right)^{T}.$ The value functions and the policies are approximated using NNs as

[TABLE]

where $\hat{W}_{ci}$ and $\hat{W}_{ai}$ are estimates of the ideal weights $W_{i}$, introduced in (22).

VIII Simulation of experience via BE extrapolation

A consequence of Theorem 1 is that the BE provides an indirect measure of how close the estimates $\hat{W}_{ci}$ and $\hat{W}_{ai}$ are to the ideal weights $W_{i}$. From a reinforcement learning perspective, each evaluation of the BE along the system trajectory can be interpreted as experience gained by the critic, and each evaluation of the BE at points not yet visited can be interpreted as simulated experience. In previous results such as [21, 20, 29, 45, 4], the critic is restricted to the experience gained (in other words BEs evaluated) along the system state trajectory. The development in [21, 20, 29, 45] can be extended to employ simulated experience; however, the extension requires exact model knowledge. In results such as [4], the formulation of the BE does not allow for simulation of experience. The formulation in (29) employs the system identifier developed in Section V to facilitate approximate evaluation of the BE at off-trajectory points.

To simulate experience, a set of points $\left\{\mathcal{E}_{i}^{k}\right\}_{k=1}^{M_{i}}$ is selected corresponding to each agent $i$ , and the instantaneous BE in (22) is approximated at the current state and at the selected points using (37). The approximation at the current state is denoted by $\hat{\delta}_{ti}$ and the approximation at the selected points is denoted by $\hat{\delta}_{ti}^{k}$, where $\hat{\delta}_{ti}$ and $\hat{\delta}_{ti}^{k}$ are defined as

[TABLE]

Note that once $\left\{e_{j}\right\}_{j\in\mathcal{S}_{i}}$ and $x_{i}$ are selected, the $i\textsuperscript{th}$ agent can compute the states of all the remaining agents in the sub-graph. For notational brevity, the arguments to the functions $\sigma_{i}$, $\hat{\mathscr{F}}_{i}$, $\mathscr{G}_{i}$, $\mathcal{G}_{i}$, $\hat{\mathcal{F}}_{i}$, $\hat{\mu}_{i}$, $G_{\sigma i}$, $G_{\epsilon i}$, and $\epsilon_{i}$ are suppressed hereafter.

The critic uses simulated experience to update the value function weights using a least squares-based update law

[TABLE]

where $\rho_{i}\triangleq 1+\nu_{i}\omega_{i}^{T}\Gamma_{i}\omega_{i}$, $\Gamma_{i}\in\mathbb{R}^{L_{i}\times L_{i}}$ denotes the time-varying least-squares learning gain, $\overline{\Gamma}_{i}\in\mathbb{R}$ denotes the saturation constant, $\left\|\Gamma_{i}\left(t_{0}\right)\right\|\leq\overline{\Gamma}_{i},$ and $\eta_{c1i},\eta_{c2i},\beta_{i},\nu_{i}\in\mathbb{R}$ are constant positive learning gains. In (34),

[TABLE]

where for a function $\phi_{i}\left(\mathcal{E}_{i},\left(\cdot\right)\right)$, the notation $\phi_{i}^{k}$ indicates evaluation at $\mathcal{E}_{i}=\mathcal{E}_{i}^{k}$; i.e., $\phi_{i}^{k}\triangleq\phi_{i}\left(\mathcal{E}_{i}^{k},\left(\cdot\right)\right)$. The actor updates the policy weights using the following update law derived based on the Lyapunov-based stability analysis in section IX.

[TABLE]

where $\eta_{a1i},\eta_{a2i}\in\mathbb{R}$ are constant positive learning gains. The following assumption facilitates simulation of experience.

Assumption 7**.**

[43] For each $i\in\mathcal{N}$, there exists a finite set of points $\left\{\mathcal{E}_{i}^{k}\right\}_{k=1}^{M_{i}}$ such that

[TABLE]

where $\lambda_{\min}$ denotes the minimum eigenvalue, and $\underline{\rho_{i}}\in\mathbb{R}$ is a positive constant.

IX Stability analysis

To facilitate the stability analysis, the left hand side of (16) is subtracted from (29) to express the BEs in terms of the weight estimation errors as

[TABLE]

where $\tilde{\left(\cdot\right)}\triangleq\left(\cdot\right)-\hat{\left(\cdot\right)}$, $\Delta_{i}=O\left(\left(\overline{\epsilon}\right)_{\mathcal{S}_{i}},\left(\overline{\nabla\epsilon}\right)_{\mathcal{S}_{i}},\left(\overline{\epsilon_{\theta}}\right)_{\mathcal{S}_{i}}\right)$, and $\mathcal{R}_{\mathcal{S}_{j}}\triangleq\mbox{diag}\left(\left[R_{\lambda_{j}^{1}}^{-1}G_{\sigma\lambda_{j}^{1}}^{T},\cdots,R_{\lambda_{j}^{s_{j}}}^{-1}G_{\sigma\lambda_{j}^{s_{j}}}^{T}\right]\right)$ are block diagonal matrices. Consider a set of extended neighbors $\mathcal{S}_{p}$ corresponding to the $p$th agent. To analyze asymptotic properties of the agents in $\mathcal{S}_{p},$ consider the following candidate Lyapunov function

[TABLE]

where $Z_{p}\in\mathbb{R}^{\left(2ns_{i}+2L_{i}s_{i}+n\left(P_{i}+1\right)s_{i}\right)}$ is defined as

[TABLE]

$\mbox{vec}\left(\cdot\right)$ denotes the vectorization operator, and $V_{ti}:\mathbb{R}^{ns_{i}}\times\mathbb{R}\to\mathbb{R}$ is defined as

[TABLE]

for all $e_{\mathcal{S}_{i}}\in\mathbb{R}^{ns_{i}}$ and for all $t\in\mathbb{R}_{\geq t_{0}}$. Since $V_{ti}^{*}$ depends on $t$ only through uniformly bounded leader trajectories, Lemma 1 from [46] can be used to show that $V_{ti}$ is a positive definite and decrescent function.222Since the graph has a spanning tree, the mapping between the errors and the states is invertible. Hence, the state of an agent can be expressed as $x_{i}=h_{i}\left(e_{\mathcal{S}_{i}},x_{0}\right)$ for some function $h_{i}$. Thus, the value function can be expressed as $V_{i}^{*}\left(e_{\mathcal{S}_{i}},x_{0}\right)=V_{i}^{*}\left(e_{\mathcal{S}_{i}},h\left(e_{\mathcal{S}_{i}},x_{0}\right)\right)$. Then, $V_{ti}^{*}$ can be alternatively defined as $V_{ti}\left(e_{\mathcal{S}_{i}},t\right)\triangleq V_{i}^{*}\left(\begin{bmatrix}e_{\mathcal{S}_{i}}\\ x_{0}\left(t\right)\end{bmatrix}\right).$ Since $x_{0}$ is a uniformly bounded function of $t$ by assumption, Lemma 1 from [46] can be used to conclude that $V_{ti}$ is a positive definite and decrescent function. Thus, using Lemma 4.3 from [47], the following bounds on the candidate Lyapunov function in (38) are established

[TABLE]

for all $Z_{p}\in\mathbb{R}^{\left(2ns_{i}+2L_{i}s_{i}+n\left(P_{i}+1\right)s_{i}\right)}$ and for all $t\in\mathbb{R}_{\geq t_{0}}$, where $\underline{v_{lp}},\overline{v_{lp}}:\mathbb{R}\to\mathbb{R}$ are class $\mathcal{K}$ functions.

To facilitate the stability analysis, given any compact ball $\chi_{p}\subset\mathbb{R}^{2ns_{i}+2L_{i}s_{i}+n\left(P_{i}+1\right)s_{i}}$ of radius $r_{p}\in\mathbb{R}$ centered at the origin, a positive constant $\iota_{p}\in\mathbb{R}$ is defined as

[TABLE]

where for any function $\varpi:\mathbb{R}^{l}\to\mathbb{R}$, $l\in\mathbb{N}$, the notation $\overline{\left\|\varpi\right\|}$ denotes $\sup_{y\in\mathcal{\chi}_{p}\cap\mathbb{R}^{l}}\left\|\varpi\left(y\right)\right\|$ and $A_{i}^{\theta}$, $B_{i}^{\theta}$, and $A_{i}^{a1}$ are uniformly bounded state-dependent terms. The following sufficient gain conditions facilitate the subsequent stability analysis.

[TABLE]

[TABLE]

[TABLE]

where $A_{ij}^{1a\theta}$, $B_{ij}^{1a\theta}$, and $A_{ji}^{1ac}$ are uniformly bounded state-dependent terms.

Theorem 2**.**

Provided Assumptions 1 - 7 hold and the sufficient gain conditions in (41)-(43) are satisfied, the controller in (33) along with the actor and critic update laws in (34) and (35), and the system identifier in (24) along with the weight update laws in (26) ensure that the local neighborhood tracking errors $e_{i}$ are ultimately bounded and that the policies $\hat{\mu}_{i}$ converge to a neighborhood around the feedback-Nash policies $\mu_{i}^{*}$ for all $i\in\mathcal{N}$.

Proof:

The time derivative of the candidate Lyapunov function in (38) is given by

[TABLE]

Using (16), (28), (31), and (37), the update laws in (34) and (35), and the definition of $V_{ti}$ in (39), the derivative in (44) can be bounded as333For a detailed derivation of the bound, see [48].

[TABLE]

Let $v_{lp}:\mathbb{R}\to\mathbb{R}$ be a class $\mathcal{K}$ function such that

[TABLE]

where $\underline{q_{i}}:\mathbb{R}\to\mathbb{R}$ are class $\mathcal{K}$ functions such that $\underline{q_{i}}\left(\left\|e\right\|\right)\leq Q_{i}\left(e\right),\>\forall e\in\mathbb{R}^{n},\>\forall i\in\mathcal{N}$. Then, the Lyapunov derivative can be bounded as

[TABLE]

for all $Z_{p}$ such that $Z_{p}\in\chi_{p}$ and $\left\|Z_{p}\right\|\geq v_{lp}^{-1}\left(\iota_{p}\right)$. Using the bounds in (40), the sufficient conditions in (41)-(43), and the inequality in (46), Theorem 4.18 in [47] can be invoked to conclude that every trajectory $Z_{p}\left(t\right)$ satisfying $\left\|Z_{p}\left(t_{0}\right)\right\|\leq\overline{v_{lp}}^{-1}\left(\underline{v_{lp}}\left(r_{p}\right)\right)$, is bounded for all $t\in\mathbb{R}_{\geq t_{0}}$ and satisfies

[TABLE]

Since the choice of the subgraph $\mathcal{S}_{p}$ was arbitrary, the neighborhood tracking errors $e_{i}$ are ultimately bounded for all $i\in\mathcal{N}$. Furthermore, the weight estimates $\hat{W}_{ai}$ converge to a neighborhood of the ideal weights $W_{i}$; hence, invoking Theorem 1, the policies $\hat{\mu}_{i}$ converge to a neighborhood of the feedback-Nash equilibrium policies $\mu_{i}^{*}$ for all $i\in\mathcal{N}$. ∎

X Simulations

This section provides a simulation example to demonstrate the applicability of the developed technique. The agents are assumed to have the communication topology as shown in Figure 1 with unit pinning gains and edge weights. The motion of the agents is described by identical nonlinear one-dimensional dynamics of the form (1) where $f_{i}\left(x_{i}\right)=\theta_{i1}x_{i}+\theta_{i2}x_{i}^{2}$, and $g_{i}\left(x_{i}\right)=\left(\cos(2x_{i1})+2\right)$ for all $i=1,\cdots,5.$ The ideal values of the unknown parameters are selected to be $\theta_{i1}=0,$ $0,$ $0.1,$ $0.5,$ and $0.2$, and $\theta_{i2}=1$, $0.5,$ $1,$ $1,$ and $1,$ for $i=1,\cdots,5$, respectively. The agents start at $x_{i}=2$ for all $i$, and their final desired locations with respect to each other are given by $x_{d12}=0.5,$ $x_{d21}=-0.5,$ $x_{d43}=-0.5,$ and $x_{d53}=-0.5$. The leader traverses an exponentially decaying trajectory $x_{0}\left(t\right)=e^{-0.1t}$. The desired positions of agents 1 and 3 with respect to the leader are $x_{d10}=0.75$ and $x_{d30}=1$, respectively.444The optimal control problem parameters, basis functions, and adaptation gains for all the agents and the plots for weight estimates corresponding to agents 1-5 are available in [48]

For each agent $i,$ five values of $e_{i}$, three values of $x_{i}$, and three values of errors corresponding to all the extended neighbors are selected for BE extrapolation, resulting in $5\times 3^{s_{i}}$ total values of $\mathcal{E}_{i}$. All agents estimate the unknown drift parameters using history stacks containing thirty points recorded online using a singular value maximizing algorithm (cf. [49]), and compute the required state derivatives using a fifth order Savitzky-Golay smoothing filter (cf. [50]).

Figures 2 - 4 show the tracking error, the state trajectories compared with the desired trajectories, and the control inputs for all the agents demonstrating convergence to the desired formation and the desired trajectory. Note that Agents 2, 4, and 5 do not have a communication link to the leader, nor do they know their desired relative position with respect to the leader. The convergence to the desired formation is achieved via cooperative control based on decentralized objectives. Figure 5 shows the evolution and convergence of the value function weights and the parameters estimates for the drift dynamics for Agent 1. The errors between the ideal drift parameters and their respective estimates are large, however, as demonstrated by Figure 3, the resulting dynamics are sufficiently close to the actual dynamics for the developed technique to generate stabilizing policies. It is unclear whether the value function and the policy weights converge to their ideal values. Since an alternative method to solve this problem is not available to the best of the author’s knowledge, a comparison between value function and policy weight estimates and their corresponding ideal values is infeasible.

XI Concluding Remarks

A simulation-based actor-critic-identifier architecture is developed to obtain feedback-Nash equilibrium solutions to a class of differential graphical games*. *It is established that in a cooperative game based on minimization of the local neighborhood tracking errors, the value function corresponding to an agent depends on information obtained from all their extended neighbors. A set of coupled HJ equations are developed that serve as necessary and sufficient conditions for feedback-Nash equilibrium, and closed-form expressions for the feedback-Nash equilibrium policies are developed based on the HJ equations. The fact that the developed technique requires each agent to communicate with all of its extended neighbors motivates the search for a decentralized method to generate feedback-Nash equilibrium policies.