METHOD, DEVICE, AND STORAGE MEDIUM FOR DECENTRALIZED OPTIMAL CONTROL FOR LARGE-SCALE MULTIAGENT SYSTEMS

Abstract

The present disclosure provides a method, a device, and a storage medium for decentralized optimal control for a large-scale multi-agent system. The method includes initializing errors to obtain an initialized error of each of an actor NN, a critic NN and a mass NN; initializing error thresholds to obtain an initialized error threshold of each of the actor NN, the critic NN and the mass NN; and if the initialized error of each of the actor NN, the critic NN and the mass NN is greater than or equal to a corresponding initialized error threshold, calculating weights of each of the actor NN, the critic NN and the mass NN, and updating the actor NN, the critic NN, and the mass NN; and calculating errors of each of the actor NN, the critic NN and the mass NN, and updating the actor NN, the critic NN, and the mass NN.

Description

FIELD OF THE DISCLOSURE

The present disclosure generally relates to the field of hierarchical heterogeneous planning and scheduling technology and, more particularly, relates to a method, a device, and a storage medium for decentralized optimal control for large-scale multi-agent systems.

BACKGROUND

In recent years, large-scale multi-agent systems (LS-MAS) have attracted significant interest from both research societies and industrial communities due to its capability of upgrading conventional multi-agent system performance by using its diversity gain. For instance, the tracking control problem in the LS-MAS has been studied. However, It is extremely difficult to directly utilize conventional control into LS-MAS due to three challenges. The first challenge is notorious “curse of dimensionality”. Since conventional cooperative control needs each agent to know other agents' states, the computational complexity of distributed control is exponentially increased along with increased number of agents. The second challenge is lacking a realistic reliable communication network that can timely support information exchange among LS-MAS. Due to the limitation of communication capability in practice, conventional distributed cooperative control techniques are extremely difficult to be applied. The last challenge is that the constraints from physical system limitation and practical environment may cause difficulty in LS-MAS optimal control design. Therefore, there is a need to overcome these challenges simultaneously and lead to an intelligent, reliable and applicable control for LS-MAS.

BRIEF SUMMARY OF THE DISCLOSURE

One aspect or embodiment of the present disclosure provides a method for decentralized optimal control for a large-scale multi-agent system. The large-scale multi-agent system includes multiple agents, and each agent includes three neural networks (NNs) including an actor NN, a critic NN, and a mass NN. The method includes initializing errors to obtain an initialized error of the actor NN, an initialized error of the critic NN, and an initialized error of the mass NN; initializing error thresholds to obtain an initialized error threshold of the actor NN, an initialized error threshold of the critic NN, and an initialized error threshold of the mass NN; if the initialized error of the actor NN is greater than or equal to the initialized error threshold of the actor NN, if the initialized error of the critic NN is greater than or equal to the initialized error threshold of the critic NN, and if the initialized error of the mass NN is greater than or equal to the initialized error threshold of the mass NN: calculating NN weights of the actor NN, the critic NN, and the mass NN, respectively; and updating the actor NN, the critic NN, and the mass NN using corresponding calculated NN weights, respectively; and calculating NN errors of the actor NN, the critic NN, and the mass NN, respectively; and updating the actor NN, the critic NN, and the mass NN using corresponding calculated NN errors, respectively.

Another aspect or embodiment of the present disclosure provides a device for decentralized optimal control for a large-scale multi-agent system. The large-scale multi-agent system includes multiple agents, and each agent includes three neural networks (NNs) including an actor NN, a critic NN, and a mass NN. The device includes a memory, configured to store program instructions for performing a method for decentralized optimal control for the large-scale multi-agent system; and a processor, coupled with the memory and, when executing the program instructions, configured for: initializing errors to obtain an initialized error of the actor NN, an initialized error of the critic NN, and an initialized error of the mass NN; initializing error thresholds to obtain an initialized error threshold of the actor NN, an initialized error threshold of the critic NN, and an initialized error threshold of the mass NN; if the initialized error of the actor NN is greater than or equal to the initialized error threshold of the actor NN, if the initialized error of the critic NN is greater than or equal to the initialized error threshold of the critic NN, and if the initialized error of the mass NN is greater than or equal to the initialized error threshold of the mass NN: calculating NN weights of the actor NN, the critic NN, and the mass NN, respectively; and updating the actor NN, the critic NN, and the mass NN using corresponding calculated NN weights, respectively; and calculating NN errors of the actor NN, the critic NN, and the mass NN, respectively; and updating the actor NN, the critic NN, and the mass NN using corresponding calculated NN errors, respectively.

Another aspect or embodiment of the present disclosure provides a non-transitory computer-readable storage medium, containing program instructions for, when being executed by a processor, performing a method for decentralized optimal control for a large-scale multi-agent system. The large-scale multi-agent system includes multiple agents, and each agent includes three neural networks (NNs) including an actor NN, a critic NN, and a mass NN. The method includes initializing errors to obtain an initialized error of the actor NN, an initialized error of the critic NN, and an initialized error of the mass NN; initializing error thresholds to obtain an initialized error threshold of the actor NN, an initialized error threshold of the critic NN, and an initialized error threshold of the mass NN; if the initialized error of the actor NN is greater than or equal to the initialized error threshold of the actor NN, if the initialized error of the critic NN is greater than or equal to the initialized error threshold of the critic NN, and if the initialized error of the mass NN is greater than or equal to the initialized error threshold of the mass NN: calculating NN weights of the actor NN, the critic NN, and the mass NN, respectively; and updating the actor NN, the critic NN, and the mass NN using corresponding calculated NN weights, respectively; and calculating NN errors of the actor NN, the critic NN, and the mass NN, respectively; and updating the actor NN, the critic NN, and the mass NN using corresponding calculated NN errors, respectively.

Other aspects or embodiments of the present disclosure may be understood by those skilled in the art in light of the description, the claims, and the drawings of the present disclosure.

BRIEF DESCRIPTION OF THE DRAWINGS

The following drawings are merely examples for illustrative purposes according to various disclosed embodiments and are not intended to limit the scope of the present disclosure.

FIG. 1 depicts a flowchart of an exemplary method for decentralized optimal control for a large-scale multi-agent system according to various disclosed embodiments of the present disclosure.

FIG. 2 depicts an exemplary barrier-actor-critic-mass (BACM) algorithm according to various disclosed embodiments of the present disclosure.

FIG. 3 depicts an exemplary structure of a barrier-actor-critic-mass system according to various disclosed embodiments of the present disclosure.

FIG. 4 depicts an exemplary overall trajectory schematic according to various disclosed embodiments of the present disclosure.

FIG. 5 depicts an exemplary tracking error plot of all agents in an x axis according to various disclosed embodiments of the present disclosure.

FIG. 6 depicts an exemplary tracking error plot of all agents in a y axis according to various disclosed embodiments of the present disclosure.

FIG. 7 depicts an exemplary HJB (Hamiltonian-Jacobi-Bellman) equation error plot according to various disclosed embodiments of the present disclosure.

FIG. 8 depicts an exemplary FPK (Fokker-Planck-Kolmogorov) equation error plot according to various disclosed embodiments of the present disclosure.

DETAILED DESCRIPTION

References may be made in detail to exemplary embodiments of the disclosure, which may be illustrated in the accompanying drawings. Wherever possible, same reference numbers may be used throughout the accompanying drawings to refer to same or similar parts.

Mean field game theory (MFG) may be adopted to address the “curse of dimensionality” in LS-MAS. In MFG, individual agents may use a probability density function (PDF) (i.e. “mass”) of all agents to observe the behavior of entire population without requiring their states and control inputs. Then, infinity players' non-cooperative game may be shifted into a two-players game that includes a single agent versus entire population. Meanwhile, practical physical system limitations as well as complex environment may cause constraints into the control design for LS-MAS. For example, both state and density constraints may be considered in MFG based control for LS-MAS, respectively. To better integrate those constraints into the MFG-based LS-MAS optimal control problem formulation, barrier function may be adopted for handling individual agent state constraint and mass function's density constraint. With the barrier function and MFG, the constrained LS-MAS optimal control problem may be formulated. However, to obtain optimal control, a pair of forward and backward partial differential equation (PDE), called Fokker-Planck-Kolmogorov (FPK) equation and Hamiltonian-Jacobi-Bellman (HJB) equation, may need to be solved. It is extremely difficult and even impossible to directly solve these PDEs since these two PDEs are closely coupled with each other. To address such difficulty, adaptive dynamic programming and reinforcement learning technique may be adopted. Furthermore, a barrier-actor-critic-mass (BACM) learning algorithm may be developed with mass NN (neural network) for learning behaviors of large population via estimating the solution of FPK equation with barrier function, critic NN for obtaining optimal cost function by learning the solution of the HJB equation with barrier function, and actor NN for solving decentralized optimal tracking control based on the information provided by the mass NN and the critic NN. The key contributions of such configuration may be the following: the boundary and density constraints may be integrated into conventional MFG based LS-MAS optimization through a barrier function based system transformation; and the barrier-actor-critic-mass algorithm may be developed to solve the constrained HJB and FPK equations simultaneously and further obtain the optimal control for LS-MAS in real-time.

According to various embodiments of the present disclosure, LS-MAS tracking optimal control is described hereinafter. N may represent the number of homogeneous agents moving in a l dimensional configuration space, which is enclosed by an upper and lower boundary. An agent i may be controlled by the stochastic differential equation with their states being constrained as follows:

$\begin{array}{cc}{\mathrm{dx}}_i=\left[f\left(x_i\right)+g\left(x_i\right)u_i\right]\mathrm{dt}+\sqrt{2v}{\mathrm{dB}}_i& \left(1\right)\end{array}$

where f(x_i) and g(x_i) may be nonlinear functions, x_i may be an agent state which includes the position and velocity of the agent, u_i may be a control input, B_i may be standard Brownian motion which represents the process noise; and v may be a non-negative parameter.

A predefined time varying trajectory x_r(t) may be given to all agents, where t is time. The objective of individual agent may be to track the reference trajectory by minimizing the tracking error which is defined as the following:

$\begin{array}{cc}{\tilde{x}}_i\left(t\right)=x_i\left(t\right)-x_r\left(t\right)& \left(2\right)\end{array}$

Moreover, the tracking error dynamics may be derived as follows:

$\begin{array}{cc}\begin{array}{c}d{\tilde{x}}_i\left(t\right)={\mathrm{dx}}_i\left(t\right)-{\mathrm{dx}}_r\left(t\right)\\ =\left[f\left(x_i\right)+g\left(x_i\right)u_i-\frac{\text{?}}{\mathrm{dt}}\right]\mathrm{dt}+\sqrt{2\upsilon }{\mathrm{dB}}_i\\ =\left[f^{\prime }\left({\tilde{x}}_i\right)+g^{\prime }\left({\tilde{x}}_i\right)u_i\right]\mathrm{dt}+\sqrt{2\upsilon }{\mathrm{dB}}_i\end{array}& \left(3\right)\end{array}$ $\mathrm{where},$ $f^{\prime }\left({\tilde{x}}_i\right)=f\left({\tilde{x}}_i+x_r\right)-\left({\mathrm{dx}}_r/\mathrm{dt}\right)$ $\mathrm{and}$ $g^{\prime }\left({\tilde{x}}_i\right)=g\left({\tilde{x}}_i+x_r\right)$ $\text{?}\text{indicates text missing or illegible when filed}$

The optimal objective of each agent may be to track the reference trajectory by minimizing the following cost function:

$\begin{array}{cc}\text{?}\left({\tilde{x}}_i,m\right)=E\left\{{\int }_0^{\infty }\left[L\left({\tilde{x}}_i,u_i\right)+C\left({\tilde{x}}_i,m\right)\right]\mathrm{dt}\right\}& \left(4\right)\end{array}$ $\text{?}\text{indicates text missing or illegible when filed}$

where m({tilde over (x)}_i, t) may denote the probability density function (mass) of the population's tracking error at time t. Also, C({tilde over (x)}_i, m) may be the mean field coupling function which represents the interaction between agent i and the whole population of other agents. Since the dimension of the PDF and each agent state are same, the mean field coupling function can greatly reduce the computational complexity problem. Moreover, L({tilde over (x)}_i, u_i)=∥{tilde over (x)}_i∥+∥u_i∥_R², where Q and R have compatible dimensions.

Next, a barrier-function based system transformation may be applied to the original system to ensure both the tracking error state and density constraints. Let the Barrier function B(.):→ is defined on (l_{{tilde over (x)},i}), u_{{tilde over (x)},i}), then the tracking error state {tilde over (x)}_i of the system may be represented as follows:

$\begin{array}{cc}{s}_i=B_i\left({\tilde{x}}_i;l_{\tilde{x},i},u_{\tilde{x},i}\right)=\ln \frac{u_{\tilde{x},i}\left(l_{\tilde{x},i}-{\tilde{x}}_i\right)}{l_{\tilde{x},i}\left(u_{\tilde{x},i}-{\tilde{x}}_i\right)}& \left(5\right)\end{array}$

• • where l_{{tilde over (x)},i} and u_{{tilde over (x)},i} may satisfy l_{{tilde over (x)},i}<{tilde over (x)}i<u_{{tilde over (x)},i} and s_i may be the tracking error state of the transformed system. Also, the Barrier function may be invertible on interval (l_{{tilde over (x)},i}, u_{{tilde over (x)},i}), for example:

$\begin{array}{cc}{\tilde{x}}_i=B_i^{-1}\left(s_i;l_{\tilde{x},i},u_{\tilde{x},i}\right)=l_{\tilde{x},i},u_{\tilde{x},i}\frac{\text{?}-\text{?}}{l_{\tilde{x},i}\text{?}-u_{\tilde{x},i}\text{?}}& \left(6\right)\end{array}$ $\text{?}\text{indicates text missing or illegible when filed}$

Similarly, barrier function may be generated for ensuring density constraint as follows:

$\begin{array}{cc}\rho =B_m\left(m;{\rho }_1,p_2\right)=\ln \frac{{\rho }_2\left({\rho }_1-m\left({\tilde{x}}_i,t\right)\right)}{{\rho }_1\left({\rho }_2-m\left({\tilde{x}}_i,t\right)\right)}& \left(7\right)\end{array}$

• • where ρ₁ and ρ₂ may be two constants satisfying ρ₁<ρ₂. The inverse of the barrier function may be represented as follows:

$\begin{array}{cc}m=B_m^{-1}\left(\rho ;{\rho }_1,{\rho }_2\right)={\rho }_1{\rho }_2\frac{\text{?}-\text{?}}{{\rho }_1\text{?}-{\rho }_2\text{?}}& \left(8\right)\end{array}$ $\text{?}\text{indicates text missing or illegible when filed}$

• • where ρ may be the density of the transformed system.

In one embodiment of the present disclosure, the barrier functions B(.) may take finite value when the arguments are within above defined region and approach to infinity as the state and density approach the boundary of the defined region, respectively.

The dynamics of the transformed state s_i may be obtained by using following chain rule:

$\begin{array}{cc}\begin{array}{c}{\mathrm{ds}}_i=\frac{\frac{d{\tilde{x}}_i}{\mathrm{dt}}}{\frac{d{\tilde{x}}_i}{{\mathrm{ds}}_i}}\mathrm{dt}\\ =\frac{\left[f^{\prime }\left({\tilde{x}}_i\right)+g^{\prime }\left({\tilde{x}}_i\right)u_i\right]\mathrm{dt}}{\frac{d\left(\text{?}\frac{\text{?}}{\text{?}}\right)}{{\mathrm{ds}}_i}}\\ =\left[f^{\prime }\left({\tilde{x}}_i\right)+g^{\prime }\left({\tilde{x}}_i\right)u_i\right]\frac{u_{\tilde{x},i}^2\text{?}-l_{\tilde{x},i}{u}_{\tilde{x},i}+l_{\tilde{x},i}^2\text{?}}{u_{\tilde{x},i}{l}_{\tilde{x},i}^2-l_{\tilde{x},i}{u}_{\tilde{x},i}^2}\mathrm{dt}\\ =\left[F\left(s_i\right)+G\left(s_i\right)u_i\right]\mathrm{dt}+\sqrt{2\upsilon }{\mathrm{dB}}_i\end{array}& \left(9\right)\end{array}$ $\mathrm{where}$ $F\left(s_i\right)=f\left({\tilde{x}}_i\right)\frac{\text{?}}{\text{?}}$ $\mathrm{and}$ $G\left(s_i\right)=g\left({\tilde{x}}_i\right)\frac{\text{?}}{\text{?}}$ $\text{?}\text{indicates text missing or illegible when filed}$

F(s_i) may be Lipschitz, and there may exist a constant a_f such that for s_i∈Ω, ∥F(s_i)∥≤a_f∥s_i∥, where Ω may be a compact set containing the origin. In addition, G(s_i) may be bounded on Ω, i.e., there may exist a constant a_g such that ∥(s_i)∥≤a_g. Moreover, the system in equation (1) may be controllable over the compact set Ω.

Next, a new cost function of the transformed state may be represented as follows:

$\begin{array}{cc}{V}_i\left(s_i,\rho \right)=E{\int }_0^{\infty }\left[L\left(s_i,u_i\right)+C\left(s_i,\rho \right)\right]\mathrm{dt}& \left(10\right)\end{array}$ $\mathrm{where}$ $L\left(s_i,u_i\right)=s_i\text{?}+{u_i}_R^2.$ $\text{?}\text{indicates text missing or illegible when filed}$

Then, a Hamiltonian may be defined as follows:

$\begin{array}{cc}H\left[s_i,{\mathrm{DV}}_i\left(s_i,\rho ,t\right)\right]=L\left(s_i,u_i\right)+{{\mathrm{DV}}_i\left(s_i,\rho ,t\right)}^T\left[F\left(s_i\right)+G\left(s_i\right)u_i\right]& \left(11\right)\end{array}$

Next, the following HJB equation may be obtained by substituting the optimal evaluation function into the Hamiltonian as follows:

$\begin{array}{cc}-{\partial }_t{V}_i^{*}\left(s_i,\rho ,t\right)-{\mathrm{\upsilon \Delta V}}_i^{*}\left(s_i,\rho ,t\right)+H\left[s_i,{\mathrm{DV}}_i^{*}\left(s_i,\rho ,t\right)\right]=C\left(s_i,\rho \right)& \left(12\right)\end{array}$

Then, the optimal control for each agent may be derived as follows:

$\begin{array}{cc}{u}_i^{*}\left(s_i\right)=-\frac{\text{?}}{2}{R}^{-1}{g}^T\left(s_i\right){\mathrm{DV}}_i^{*}\left(s_i,\rho ,t\right)& \left(13\right)\end{array}$ $\text{?}\text{indicates text missing or illegible when filed}$

To obtain the HJB equation in equation (12), the practical probability density function (PDF) (i.e., mass function p) may be required. The mass function may be obtained by solving the FPK equation, where the FPK equation with density constraint may be obtained as follows:

$\begin{array}{cc}\begin{array}{c}{\partial }_t\rho \left(s_i,t\right)=\frac{{\partial }_{t}m\left(\tilde{x},t\right)}{\frac{\partial m\left(\tilde{x},t\right)}{\text{?}}}\\ =\frac{\mathrm{\upsilon \Delta }m\left(\tilde{x},t\right)+\mathrm{di}\text{?}\left(m{D}_{p}H\left[\tilde{x},{\mathrm{DV}}_i\left(\tilde{x},m,t\right)\right]\right)}{\frac{\text{?}\frac{\text{?}}{\text{?}}}{\text{?}}}\\ =\left(\upsilon \Delta \rho \left(s_i,t\right)+\mathrm{di}\text{?}\left(\rho D_{p}H\left[s_i,{\mathrm{DV}}_i\left(s_i,\rho ,t\right)\right]\right)\right)\\ \frac{{\rho }_1^2\text{?}-2{\rho }_1{\rho }_2+{\rho }_1^2\text{?}}{{\rho }_2{\rho }_1^2-{\rho }_1{\rho }_2^2}\end{array}& \left(14\right)\end{array}$ $\text{?}\text{indicates text missing or illegible when filed}$

Next, the FPK equation with the optimal cost function may be obtained as follows:

$\begin{array}{c}\left(15\right)\end{array}$ ${\partial }_t\rho \left(s_i,t\right)-[\upsilon \Delta \rho \left(s_i,t\right)+\mathrm{di}\text{?}\left(\rho D_{p}H\left[s_i,{\mathrm{DV}}_i^{*}\left(s_i,\rho ,t\right)\right]\right))]\frac{{\rho }_1^2\text{?}-2{\rho }_1{\rho }_2+{\rho }_1^2\text{?}}{{\rho }_2{\rho }_1^2-{\rho }_1{\rho }_2^2}=0$ $\text{?}\text{indicates text missing or illegible when filed}$

According to various embodiments of the present disclosure, to obtain the optimal control policy, the coupled HJB-FPK equation may need to be solved in real time. However, the HJB and FPK equations may be multi-dimensional nonlinear PDEs whose solution may be difficult to achieve with state and density constraints. Therefore, in the present disclosure, the barrier-actor-critic-mass based NNs may be developed to learn the solution of coupled HJB-FPK equations.

According to various embodiments of the present disclosure, a method for decentralized optimal control for a large-scale multi-agent system is described hereinafter.

FIG. 1 depicts a flowchart of an exemplary method for decentralized optimal control for a large-scale multi-agent system according to various disclosed embodiments of the present disclosure. FIG. 2 depicts an exemplary BACM algorithm according to various disclosed embodiments of the present disclosure. FIG. 3 depicts an exemplary structure of a barrier-actor-critic-mass system according to various disclosed embodiments of the present disclosure.

The large-scale multi-agent system includes multiple agents; and each agent includes three neural networks (NNs) including an actor NN, a critic NN, and a mass NN. Referring to FIGS. 1-3, the method includes initializing errors to obtain an initialized error of the actor NN, an initialized error of the critic NN, and an initialized error of the mass NN; initializing error thresholds to obtain an initialized error threshold of the actor NN, an initialized error threshold of the critic NN, and an initialized error threshold of the mass NN; if the initialized error of the actor NN is greater than or equal to the initialized error threshold of the actor NN, if the initialized error of the critic NN is greater than or equal to the initialized error threshold of the critic NN, and if the initialized error of the mass NN is greater than or equal to the initialized error threshold of the mass NN: calculating NN weights of the actor NN, the critic NN, and the mass NN, respectively; and updating the actor NN, the critic NN, and the mass NN using corresponding calculated NN weights, respectively; and calculating NN errors of the actor NN, the critic NN, and the mass NN, respectively; and updating the actor NN, the critic NN, and the mass NN using corresponding calculated NN errors, respectively.

In one embodiment, the method further includes, if the initialized error of the actor NN is less than the initialized error threshold of the actor NN, obtaining previous calculated NN weights of the actor NN; or if the initialized error of the critic NN is less than the initialized error threshold of the critic NN, obtaining previous calculated NN weights of the critic NN; or if the initialized error of the mass NN is less than the initialized error threshold of the mass NN, obtaining previous calculated NN weights of the mass NN.

In one embodiment, the method further includes using the previous calculated NN weights of the actor NN to calculate a control; and executing the calculated control.

In one embodiment, the method further includes, before initializing the errors, initializing a state and a density of the agent, where the state of the agent includes a position and a velocity; and calculating an error of the agent using the state of the agent and a predefined trajectory.

In one embodiment, the method further includes, before initializing the errors and after calculating the error of the agent, performing a barrier-function based system transformation on the error and the density of the agent to obtain to a transformed error state and a transformed density state, respectively.

In one embodiment, the transformed error state and the transformed density state are configured to calculate corresponding NN weights and errors.

In one embodiment, the method further includes, before initializing the errors, randomly initializing the NN weights of the actor NN, the critic NN, and the mass NN.

In one embodiment, the critic NN is configured to estimate a cost function; and the mass NN is configured to estimate a probability density function.

In one embodiment, the agent includes an unmanned aerial vehicle.

In one embodiment, referring to FIG. 2, there may be no need for using two levels of while-loops; the outside loop may be only for the purpose of debugging; and optionally, lines 7 and 15 may be removed or omitted.

According to various embodiments of the present disclosure, the barrier-actor-critic-mass algorithm is described hereinafter. Referring to FIGS. 1-3, in the BACM, each agent may maintain three neural networks (NN). The actor NN may approximate the optimal control policy, the critic NN may approximate the optimal evaluation function and the mass NN may estimate the density of the entire population. Meanwhile, the barrier function may be applied into three NNs to ensure both tracking error and density constraints being satisfied during the learning process.

According to various embodiments of the present disclosure, critic learning is described in the following. The optimal value function may be represented as follows:

$\begin{array}{cc}\text{?}\left(s_i,\rho ,t\right)=W_{V,i}^T{\varphi }_{V,i}+{\epsilon }_{V,i}& \left(16\right)\end{array}$ $\text{?}\text{indicates text missing or illegible when filed}$

where W_v,i may be an ideal critic NN weight and ∅_v,i may be the critic NN activation function. In addition, e may represent the reconstruction error of critic neural network. Next, the optimal cost function may be approximated as follows:

$\begin{array}{cc}{\hat{V}}_i\left(s_i,{\hat{\rho }}_i,t\right)={\hat{W}}_{V,i}^T{\hat{\varphi }}_{V,i}& \left(17\right)\end{array}$

where Ŵ_V,i may be the approximated NN weights.

By substituting equation (17) to equation (12), a residual error used to tune the weight of the critic NN may be obtained as follows:

$\begin{array}{cc}{e}_{\mathrm{HJBi}}=C\left(s_i,{\hat{\rho }}_i\right)+{\hat{W}}_{V,i}^T\left[{\partial }_t{\hat{\varphi }}_{V,i}+\mathrm{\nu \Delta }{\hat{\varphi }}_{V,i}-{\hat{H}}_W\right]\mathrm{where},\hat{H}=H\left[s_i,{\mathrm{DV}}_i^{*}\left(s_i,{\rho }_i,t\right)\right]\mathrm{and}\hat{H}={\hat{W}}_{V,i}^T{\hat{H}}_W.& \left(18\right)\end{array}$

Next, the equation (18) may be simplified as follows:

$\begin{array}{cc}{e}_{\mathrm{HJBi}}=C\left(s_i,{\hat{\rho }}_i\right)+{\hat{W}}_{V,i}^T{\psi }_{V,i}\left(s_i,{\hat{\rho }}_i,t\right)& \left(19\right)\end{array}$ $\begin{array}{cc}\mathrm{where}{\psi }_{V,i}\left(s_i,{\hat{\rho }}_i,t\right)=\left[{\partial }_t{\hat{\varphi }}_{V,i}+\mathrm{\nu \Delta }{\hat{\varphi }}_{V,i}-{\hat{H}}_W\right]& \left(20\right)\end{array}$

By substituting the optimal cost function from equation (16) to equation (12), it may obtain:

$\begin{array}{cc}C\left(s_i,{\rho }_i\right)+W_{V,i}^T\left[{\partial }_t{\varphi }_{V,i}+{\mathrm{\nu \Delta \varphi }}_{V,i}-H_W\right]+{\epsilon }_{\mathrm{HJBi}}=0& \left(21\right)\end{array}$

where H=W_V,i^TH_W and ε_HJBi may be an error caused by the reconstruction error.

After the simplification, the equation (21) may be written as follows:

$\begin{array}{cc}C\left(s_i,{\rho }_i\right)+W_{V,i}^T{\psi }_{V,i}\left(s_i,{\rho }_i,t\right)+{\epsilon }_{\mathrm{HJBi}}=0& \left(22\right)\end{array}$

The approximation error of the coupling function may be derived as follows:

$\begin{array}{cc}\tilde{C}\left(s_i,{\tilde{\rho }}_i\right)=C\left(s_i,{\tilde{\rho }}_i\right)-C\left(s_i,{\rho }_i\right)& \left(23\right)\end{array}$

By substituting equation (23) to equation (22), it may obtain:

$\begin{array}{cc}C\left(s_i,{\tilde{\rho }}_i\right)-\tilde{C}\left(s_i,{\tilde{\rho }}_i\right)+W_{V,i}^T{\psi }_{V,i}\left(s_i,{\rho }_i,t\right)+{\epsilon }_{\mathrm{HJBi}}=0& \left(24\right)\end{array}$

Next, by substituting equation (24) into equation (19), it may obtain:

$\begin{array}{cc}{e}_{\mathrm{HJBi}}=\tilde{C}\left(s_i,{\tilde{\rho }}_i\right)-W_{V,i}^T{\psi }_{V,i}\left(s_i,{\rho }_i,t\right)-{\epsilon }_{\mathrm{HJBi}}+{\hat{W}}_{V,i}^T{\psi }_{V,i}\left(s_i,{\hat{\rho }}_i,t\right)& \left(25\right)\end{array}$

Next, the critic NN weight approximation error and HJB equation approximation error may be respectively defined as follows:

$\begin{array}{cc}{\hat{W}}_{V,i}=W_{V,i}-{\hat{W}}_{V,i}& \left(26\right)\end{array}$ $\begin{array}{cc}{\tilde{\psi }}_{V,i}\left(s_i,{\tilde{\rho }}_i,t\right)={\psi }_{V,i}\left(s_i,{\rho }_i,t\right)-{\psi }_{V,i}\left(s_i,{\hat{\rho }}_i,t\right)& \left(27\right)\end{array}$

By substituting equation (26) and equation (27) into equation (25), it may obtain:

$\begin{array}{cc}\begin{array}{c}{e}_{\mathrm{HJBi}}=\tilde{C}\left(s_i,{\tilde{\rho }}_i\right)-W_{V,i}^T\left({\tilde{\psi }}_{V,i}\left(s_i,{\tilde{\rho }}_i,t\right)+{\psi }_{V,i}\left(s_i,{\hat{\rho }}_i,t\right)\right)+\\ \left(W_{V,i}^T-{\tilde{W}}_{V,i}^T\right){\psi }_{V,i}\left(s_i,{\hat{\rho }}_i,t\right)-{\epsilon }_{\mathrm{HJBi}}\\ =\tilde{C}\left(s_i,{\tilde{\rho }}_i\right)-W_{V,i}^T{\tilde{\psi }}_{V,i}\left(s_i,{\tilde{\rho }}_i,t\right)-\\ {\tilde{W}}_{V,i}^T{\psi }_{V,i}\left(s_i,{\hat{\rho }}_i,t\right)-{\epsilon }_{\mathrm{HJBi}}\end{array}& \left(28\right)\end{array}$

Next, the update law for critic NN may be obtained by using the gradient descent along with the HJB approximation error as follows:

$\begin{array}{cc}{\stackrel{.}{\hat{W}}}_{V,i}=-{\alpha }_{V,i}\frac{{\psi }_{V,i}\left(s_i,{\hat{\rho }}_i,t\right)e_{\mathrm{HJBi}}^T}{1+{{\psi }_{V,i}\left(s_i,{\hat{\rho }}_i,t\right)}^2}& \left(29\right)\end{array}$

where a_V,i may be the learning rate.

According to various embodiments of the present disclosure, mass learning is described as the following. The mass function be represented as follows:

$\begin{array}{cc}\rho \left(s_i,t\right)=W_{\rho ,i}^T{\varphi }_{\rho ,i}+{\epsilon }_{\rho ,i}& \left(30\right)\end{array}$

• • where Ŵ_ρ,i and ϕ_ρ,i may be the ideal mass NN weight and activation function, respectively. ε_ρ,i may be the reconstruction error of the mass NN.

Then, the mass distribution may be estimated as follows:

$\begin{array}{cc}\hat{\rho }\left(s_i,{\stackrel{\_}{\rho }}_i,t\right)={\hat{W}}_{\rho ,i}^T{\hat{\varphi }}_{\rho ,i}& \left(31\right)\end{array}$

• • where Ŵ_ρ,i may be the approximated mass NN weight. Moreover, ρ may be the averaged historical density defined as

$\stackrel{\_}{\rho }=\frac{1}{T}{\int }_t^{t-T}\rho d\rho ,$

and T may be a constant historical window.

The residual error for the mass NN may be defined by substituting equation (31) to equation (15) as follows:

$\begin{array}{cc}{e}_{{\mathrm{FPK}}_i}={\hat{W}}_{\rho ,i}^T\left[{\partial }_t{\hat{\varphi }}_{\rho ,i}-\left[\mathrm{\nu \Delta }{\hat{\varphi }}_{\rho ,i}+\mathrm{div}\left({\hat{\varphi }}_{\rho ,i}{D}_p\hat{H}\right)\right]\frac{{\rho }_1^2{e}^{-{\stackrel{\_}{\rho }}_i}-2{\rho }_1{\rho }_2+{\rho }_1^2{e}^{{\stackrel{\_}{\rho }}_i}}{{\rho }_2{\rho }_1^2-{\rho }_1{\rho }_2^2}\right]& \left(32\right)\end{array}$ $\mathrm{where}\hat{H}=H\left[s_i,D_s{\hat{\varphi }}_{V,i}\right]$

Equation (32) may be simplified as follows:

$\begin{array}{cc}{e}_{\mathrm{FPKi}}={\hat{W}}_{\rho ,i}^T{\psi }_{\rho ,i}\left(s_i,{\stackrel{\_}{\rho }}_i,{\hat{V}}_i,t\right)& \left(33\right)\end{array}$ $\begin{array}{cc}{\psi }_{\rho ,i}\left(s_i,\rho {\text{?}}_i,{\hat{V}}_i,\text{⁠}t\right)=\left[\text{⁠}{\partial }_t{\hat{\varphi }}_{\rho ,i}-\left[v\Delta {\hat{\varphi }}_{\rho ,i}+\mathrm{div}\left({\hat{\varphi }}_{\rho ,i}{D}_{\rho }\hat{H}\right)\right]\frac{{\rho }_1^2{e}^{-{\rho }_i}\text{?}-2{\rho }_1{\rho }_2+{\rho }_1^2{e}^{-\rho }\text{?}}{{\rho }_2{\rho }_1^2-{\rho }_1{\rho }_2^2}\right]& \left(34\right)\end{array}$ $\text{?}\text{indicates text missing or illegible when filed}$

Next, by substituting the mass function from equation (30) to equation (15), it may obtain:

$\begin{array}{cc}{W}_{\rho ,i}^T[{\partial }_t{\varphi }_{\rho ,i}-\left[v\Delta {\varphi }_{\rho ,i}+\mathrm{div}\left({\varphi }_{\rho ,i}{D}_{\rho }H\left[s_i,D_s{V}_i\left(s_i,{\rho }_i,t\right)\right]\right)\right]& \left(35\right)\end{array}$ $\frac{{\rho }_1^2{e}^{-{\rho }_i}\text{?}-2{\rho }_1{\rho }_2+{\rho }_1^2{e}^{-{\rho }_i}\text{?}}{{\rho }_2{\rho }_1^2-{\rho }_1{\rho }_2^2}]+{\text{?}}_{\mathrm{FPKi}}=0$ $\text{?}\text{indicates text missing or illegible when filed}$

• • which may be simplified as follows:

$\begin{array}{cc}{W}_{\rho ,i}^T{\psi }_{\rho ,i}\left(s_i,\rho {\text{?}}_i,V_i,t\right)+{\text{?}}_{\mathrm{FPKi}}=0& \left(36\right)\end{array}$ $\text{?}\text{indicates text missing or illegible when filed}$

The mass NN weight approximation error and FPK equation approximation error may be defined as follows:

$\begin{array}{cc}{\tilde{W}}_{\rho ,i}=W_{\rho ,i}={\hat{W}}_{\rho ,i}& \left(37\right)\end{array}$ $\begin{array}{cc}\psi {\text{?}}_{\rho ,i}\left(s_i,\rho {\text{?}}_i,{\tilde{V}}_i,t\right)={\psi }_{\rho ,i}\left(s_i,\rho {\text{?}}_i,V_i,t\right)-{\psi }_{\rho ,i}\left(s_i,\rho {\text{?}}_i,{\hat{V}}_i,t\right)& \left(38\right)\end{array}$ $\text{?}\text{indicates text missing or illegible when filed}$

Next, be substituting equation (36) into equation (33), it may obtain:

$\begin{array}{cc}{\text{?}}_{\mathrm{FPKi}}=-W{\text{?}}_{\rho ,i}^T{\psi }_{\rho ,i}\left(s_i,{\stackrel{\_}{\rho }}_i,{\hat{V}}_i,t\right)-W_{\rho ,i}^T\psi {\text{?}}_{\rho ,i}\left(s_i,\rho {\text{?}}_i,V{\text{?}}_i,t\right)={\text{?}}_{\mathrm{FPKi}}& \left(39\right)\end{array}$ $\text{?}\text{indicates text missing or illegible when filed}$

Then, by applying the gradient descent along with FPK estimation error, the update law for mass NN may be generated as follows:

$\begin{array}{cc}{\stackrel{.}{\hat{W}}}_{\rho ,i}={\alpha }_{\rho .i}\frac{{\psi }_{\rho ,i}\left(s_i,{\stackrel{\_}{\rho }}_i,{\hat{V}}_i,t\right)e_{\mathrm{FPKi}}^T}{1+{{\psi }_{\rho ,i}\left(s_i,{\stackrel{\_}{\rho }}_i,{\hat{V}}_i,t\right)}^2}& \left(40\right)\end{array}$

where α_ρ,i may be the mass NN learning rate.

According to various embodiments of the present disclosure, actor learning is described as the following. The optimal control may be represented as follows:

$\begin{array}{cc}{u}_i\left(s_i,\rho ,t\right)=W_{u,i}^T{\varphi }_{u,i}+{\epsilon }_{u,i}& \left(41\right)\end{array}$

where W_u,i and ϕ_u,i may be the ideal actor NN weight and activation function, respectively. ε_u,i may be the reconstruction error of the Actor NN.

Then, the optimal control may be estimated as follows:

$\begin{array}{cc}{\hat{u}}_i\left(s_i,{\hat{\rho }}_i,t\right)={\hat{W}}_{u,i}^T{\hat{\varphi }}_{u,i}& \left(42\right)\end{array}$

where Ŵ_u,i is the approximated actor NN weight.

The residual error after substituting equation (42) into equation (13) may be represented as follows:

$\begin{array}{cc}{e}_{u,i}={\hat{W}}_{u,i}^T{\hat{\varphi }}_{u,i}+\frac{1}{2}{R}^{-1}{g}^T\left(s_i\right)D_s{\hat{V}}_i\left(s_i,{\hat{\rho }}_i,t\right)& \left(43\right)\end{array}$

Furthermore, the update law for actor NN may be designed as follows:

$\begin{array}{cc}{W}_{u,i}\text{?}=-{\alpha }_{u,i}\frac{{\varphi }_{u,i}\left(s_i,{\hat{\rho }}_i,t\right)e_{\mathrm{ui}}^T}{1+{{\varphi }_{u,i}\left(s_i,{\hat{\rho }}_i,t\right)}^2}& \left(44\right)\end{array}$ $\text{?}\text{indicates text missing or illegible when filed}$

The designed BACM algorithm has been implemented into the large-scale multi-UAV (unmanned aerial vehicle) system to address the decentralized mean field based optimal tracking control problem. In one embodiment, a total of 3000 agents (e.g., UAV) may be deployed with system dynamics under physical limitation and uncertain environment. A reference trajectory may have been given ahead of the mission planning. The goal of each agent may be to track the reference trajectory while avoiding the obstacle during the mission. Therefore, the movements of all agents may be limited to a fixed area with specific boundary and density constraint. The initial positions of all agents may be generated randomly following a normal distribution with mean 0.5 and variance 0.16. The initial velocities of all agents may be set to zero. In one embodiment, the reference trajectory may be given as follows:

$x_r\left(t\right)=\left[\begin{array}{c}0.2\sin \left(2t\right)+0.002t^2+0.5\\ 0.2t\\ 0.2\cos \left(2t\right)+0.004t\\ 0.2\end{array}\right]$

In one embodiment, the agent intrinsic dynamics may be given as follows:

$f^{\prime }\left(x\right)=\left[\begin{array}{c}{x}_2-x_1\\ x_4-x_3\\ \frac{x_2}{2}\left[\left({\cos \left(2x_1+2\right)}^2-1\right)\right]-\frac{x_1}{2}\\ \frac{x_4}{2}\left[\left({\cos \left(2x_3+2\right)}^2-1\right)\right]-\frac{x_3}{2}\end{array}\right]$ $g^{\prime }\left(x\right)=\left[\begin{array}{c}0\\ 0\\ \cos \left(2x_1\right)+2\\ \cos \left(2x_3\right)+2\end{array}\right]$

The non-negative parameter v may be selected as 0.02. The mean field cost function may be selected as C(s_i, m)=∥s_i−(ρ)∥, which represents the difference between current tracking error of the agent i and current average tracking error of the whole population. In addition, the state and density constraints may be considered as follows:

$\left[\begin{array}{c}{l}_{\tilde{x},i}\\ u_{\tilde{x},i}\end{array}\right]=\left[\begin{array}{c}{x}_r\left(t\right)+1\\ -\left(x_r\left(t\right)+1\right)\end{array}\right]$

where l_{{tilde over (x)},i} and u_{{tilde over (x)},i} may be the lower and upper bound of the state constraint, respectively.

Furthermore, the lower and upper bound of the density constraint may be defined as follows:

$\left[\begin{array}{c}{\rho }_1\\ {\rho }_2\end{array}\right]=\left[\begin{array}{c}0.2\\ 1\end{array}\right]$

• • where ρ(s)=1 may denote that the tracking error of all agents are same.

The barrier function-based system transformation may have been employed for state constraint. The new dynamics of the transformed system may be given as follows:

$\begin{array}{cc}F\left(s_i\right)=f\left(B_i^{-1}\left(s_i\right)\right)\frac{u_{\stackrel{\_}{x},i}^2{e}^{-s_i}-2l_{\stackrel{\_}{x},i}{u}_{\stackrel{\_}{x},i}+l_{\stackrel{\_}{x},i}^2{e}^{s_i}}{u_{\stackrel{\_}{x},i}{l}_{\stackrel{\_}{x},i}^2-l_{\stackrel{\_}{x},i}{u}_{\stackrel{\_}{x},i}^2}& \left(45\right)\end{array}$ $\begin{array}{cc}G\left(s_i\right)=g\left(B_i^{-1}\left(s_i\right)\right)\frac{u_{\stackrel{\_}{x},i}^2{e}^{-s_i}-2l_{\stackrel{\_}{x},i}{u}_{\stackrel{\_}{x},i}+l_{\stackrel{\_}{x},i}^2{e}^{s_i}}{u_{\stackrel{\_}{x},i}{l}_{\stackrel{\_}{x},i}^2-l_{\stackrel{\_}{x},i}{u}_{\stackrel{\_}{x},i}^2}& \left(46\right)\end{array}$ $\mathrm{where}f\left(B_i^{-1}\left(s_i\right)\right)=f^{\prime }\left(x\right)-\left({\mathrm{dx}}_r/\mathrm{dt}\right)\mathrm{and}g\left(B_i^{-1}\left(s_i\right)\right)=g^{\prime }\left(x\right).$

In one embodiment, the coefficients to evaluate the cost of actions and tracking errors may be selected as R=1, and Q=1. The learning rate of the neural network may be defined as α_u,i=2×10⁻⁴, α_V,i=2×10⁻⁶, α_ρ,i=1×10⁻³. Furthermore, the thresholds may be defined as δ_u=1×10⁻³, δ_FPK=1×10⁻³, and δ_HJB=1×10⁻⁴.

According to various embodiments of the present disclosure, the overall performance schematic of developed BACM based decentralized optimal tracking control is shown in FIG. 4. Referring to FIG. 4, the black curve may mark the reference trajectory and grey curves may represent the boundary constraints. It should be noted that the developed algorithm may force all the agents to track the reference trajectory while satisfying the given state constraints (i.e., boundary).

The tracking errors of all agents has been analyzed in various embodiments of the present disclosure. FIG. 5 depicts an exemplary tracking error plot of all agents in an x axis according to various disclosed embodiments of the present disclosure. FIG. 6 depicts an exemplary tracking error plot of all agents in a y axis according to various disclosed embodiments of the present disclosure. FIGS. 5-6 illustrate the tracking errors of all agents in the x-axis and y-axis, respectively. Both figures show that the tracking errors may converge to near zero along with time, which may indicate that the designed algorithm may track the reference trajectory in real time.

According to various embodiments of the present disclosure, the neural networks performance may be demonstrated by analyzing the HJB equation error along with the FPK equation error of agents. FIG. 7 depicts an exemplary HJB equation error plot according to various disclosed embodiments of the present disclosure. FIG. 8 depicts an exemplary FPK equation error plot according to various disclosed embodiments of the present disclosure. Without loss of generality, the optimality for, for example, agent 1, may be evaluated. Referring to FIGS. 7-8, the mean field equations error for agent 1 may converge to near zero, which may indicate that the solution of the HJB-FPK coupled equation system may be successfully approximated, such that the e-Nash equilibrium may have been reached.

According to various embodiments of the present disclosure, the BACM framework may have been developed based on mean field game theory. The decentralized optimal control for LS-MAS may have been obtained by solving the coupled HJB-FPK equations under the state and density constraints that is ensured through appropriate barrier functions. Three neural networks may be employed to solve the barrier function based mean field game, where the actor NN is for learning optimal control, the critic NN is for estimating optimal cost function, and the mass NN is for approximating the LS-MAS's probability density function (i.e., mass). Furthermore, a series of numerical simulations may have demonstrated the effectiveness of the developed method in embodiments of the present disclosure.

Various embodiments of the present disclosure further provide a device for decentralized optimal control for a large-scale multi-agent system. The large-scale multi-agent system includes multiple agents. Each agent includes three neural networks (NNs) including an actor NN, a critic NN, and a mass NN. The device includes a memory, configured to store program instructions for performing a method for decentralized optimal control for the large-scale multi-agent system; and a processor, coupled with the memory and, when executing the program instructions, configured for: initializing errors to obtain an initialized error of the actor NN, an initialized error of the critic NN, and an initialized error of the mass NN; initializing error thresholds to obtain an initialized error threshold of the actor NN, an initialized error threshold of the critic NN, and an initialized error threshold of the mass NN; if the initialized error of the actor NN is greater than or equal to the initialized error threshold of the actor NN, if the initialized error of the critic NN is greater than or equal to the initialized error threshold of the critic NN, and if the initialized error of the mass NN is greater than or equal to the initialized error threshold of the mass NN: calculating NN weights of the actor NN, the critic NN, and the mass NN, respectively; and updating the actor NN, the critic NN, and the mass NN using corresponding calculated NN weights, respectively; and calculating NN errors of the actor NN, the critic NN, and the mass NN, respectively; and updating the actor NN, the critic NN, and the mass NN using corresponding calculated NN errors, respectively.

Various embodiments of the present disclosure further provide a non-transitory computer-readable storage medium, containing program instructions for, when being executed by a processor, performing a method for decentralized optimal control for a large-scale multi-agent system. The large-scale multi-agent system includes multiple agents. Each agent includes three neural networks (NNs) including an actor NN, a critic NN, and a mass NN. The method includes initializing errors to obtain an initialized error of the actor NN, an initialized error of the critic NN, and an initialized error of the mass NN; initializing error thresholds to obtain an initialized error threshold of the actor NN, an initialized error threshold of the critic NN, and an initialized error threshold of the mass NN; if the initialized error of the actor NN is greater than or equal to the initialized error threshold of the actor NN, if the initialized error of the critic NN is greater than or equal to the initialized error threshold of the critic NN, and if the initialized error of the mass NN is greater than or equal to the initialized error threshold of the mass NN: calculating NN weights of the actor NN, the critic NN, and the mass NN, respectively; and updating the actor NN, the critic NN, and the mass NN using corresponding calculated NN weights, respectively; and calculating NN errors of the actor NN, the critic NN, and the mass NN, respectively; and updating the actor NN, the critic NN, and the mass NN using corresponding calculated NN errors, respectively.

The embodiments disclosed herein may be exemplary only. Other applications, advantages, alternations, modifications, or equivalents to the disclosed embodiments may be obvious to those skilled in the art and be intended to be encompassed within the scope of the present disclosure.

Claims

1. A method for decentralized optimal control for a large-scale multi-agent system, the large-scale multi-agent system including multiple agents each including three neural networks (NNs) including an actor NN, a critic NN, and a mass NN, the method comprising: initializing errors to obtain an initialized error of the actor NN, an initialized error of the critic NN, and an initialized error of the mass NN;initializing error thresholds to obtain an initialized error threshold of the actor NN, an initialized error threshold of the critic NN, and an initialized error threshold of the mass NN; andif the initialized error of the actor NN is greater than or equal to the initialized error threshold of the actor NN, if the initialized error of the critic NN is greater than or equal to the initialized error threshold of the critic NN, and if the initialized error of the mass NN is greater than or equal to the initialized error threshold of the mass NN: calculating NN weights of the actor NN, the critic NN, and the mass NN, respectively; and updating the actor NN, the critic NN, and the mass NN using corresponding calculated NN weights, respectively; andcalculating NN errors of the actor NN, the critic NN, and the mass NN, respectively; and updating the actor NN, the critic NN, and the mass NN using corresponding calculated NN errors, respectively.

2. The method according to claim 1, further including: if the initialized error of the actor NN is less than the initialized error threshold of the actor NN, obtaining previous calculated NN weights of the actor NN; orif the initialized error of the critic NN is less than the initialized error threshold of the critic NN, obtaining previous calculated NN weights of the critic NN; orif the initialized error of the mass NN is less than the initialized error threshold of the mass NN, obtaining previous calculated NN weights of the mass NN.

3. The method according to claim 2, further including: using the previous calculated NN weights of the actor NN to calculate a control; and

4. The method according to claim 1, wherein before initializing the errors, the method further includes: initializing a state and a density of the agent, wherein the state of the agent includes a position and a velocity; andcalculating an error of the agent using the state of the agent and a predefined trajectory.

5. The method according to claim 4, wherein before initializing the errors and after calculating the error of the agent, the method further includes: performing a barrier-function based system transformation on the error and the density of the agent to obtain to a transformed error state and a transformed density state, respectively.

6. The method according to claim 5, wherein: the transformed error state and the transformed density state are configured to calculate corresponding NN weights and errors.

7. The method according to claim 1, wherein before initializing the errors, the method further includes: randomly initializing the NN weights of the actor NN, the critic NN, and the mass NN.

8. The method according to claim 1, wherein: the critic NN is configured to estimate a cost function; andthe mass NN is configured to estimate a probability density function.

9. The method according to claim 1, wherein: the agent includes an unmanned aerial vehicle.

10. A device for decentralized optimal control for a large-scale multi-agent system, the large-scale multi-agent system including multiple agents each including three neural networks (NNs) including an actor NN, a critic NN, and a mass NN, the device comprising: a memory, configured to store program instructions for performing a method for decentralized optimal control for the large-scale multi-agent system; anda processor, coupled with the memory and, when executing the program instructions, configured for: initializing errors to obtain an initialized error of the actor NN, an initialized error of the critic NN, and an initialized error of the mass NN;initializing error thresholds to obtain an initialized error threshold of the actor NN, an initialized error threshold of the critic NN, and an initialized error threshold of the mass NN;if the initialized error of the actor NN is greater than or equal to the initialized error threshold of the actor NN, if the initialized error of the critic NN is greater than or equal to the initialized error threshold of the critic NN, and if the initialized error of the mass NN is greater than or equal to the initialized error threshold of the mass NN:calculating NN weights of the actor NN, the critic NN, and the mass NN, respectively; and updating the actor NN, the critic NN, and the mass NN using corresponding calculated NN weights, respectively; andcalculating NN errors of the actor NN, the critic NN, and the mass NN, respectively; and updating the actor NN, the critic NN, and the mass NN using corresponding calculated NN errors, respectively.

11. The device according to claim 10, wherein the processor is further configured for: if the initialized error of the actor NN is less than the initialized error threshold of the actor NN, obtaining previous calculated NN weights of the actor NN; orif the initialized error of the critic NN is less than the initialized error threshold of the critic NN, obtaining previous calculated NN weights of the critic NN; orif the initialized error of the mass NN is less than the initialized error threshold of the mass NN, obtaining previous calculated NN weights of the mass NN.

12. The device according to claim 11, wherein the processor is further configured for: using the previous calculated NN weights of the actor NN to calculate a control; andexecuting the calculated control.

13. The device according to claim 10, wherein before initializing the errors, the processor is further configured for: initializing a state and a density of the agent, wherein the state of the agent includes a position and a velocity; andcalculating an error of the agent using the state of the agent and a predefined trajectory.

14. The device according to claim 13, wherein before initializing the errors and after calculating the error of the agent, the processor is further configured for: performing a barrier-function based system transformation on the error and the density of the agent to obtain to a transformed error state and a transformed density state, respectively.

15. The device according to claim 14, wherein: the transformed error state and the transformed density state are configured to calculate corresponding NN weights and errors.

16. The device according to claim 10, wherein before initializing the errors, the processor is further configured for: randomly initializing the NN weights of the actor NN, the critic NN, and the mass NN.

17. A non-transitory computer-readable storage medium, containing program instructions for, when being executed by a processor, performing a method for decentralized optimal control for a large-scale multi-agent system which includes multiple agents each including three neural networks (NNs) including an actor NN, a critic NN, and a mass NN, the method comprising: initializing errors to obtain an initialized error of the actor NN, an initialized error of the critic NN, and an initialized error of the mass NN;initializing error thresholds to obtain an initialized error threshold of the actor NN, an initialized error threshold of the critic NN, and an initialized error threshold of the mass NN;if the initialized error of the actor NN is greater than or equal to the initialized error threshold of the actor NN, if the initialized error of the critic NN is greater than or equal to the initialized error threshold of the critic NN, and if the initialized error of the mass NN is greater than or equal to the initialized error threshold of the mass NN: calculating NN weights of the actor NN, the critic NN, and the mass NN, respectively; and updating the actor NN, the critic NN, and the mass NN using corresponding calculated NN weights, respectively; andcalculating NN errors of the actor NN, the critic NN, and the mass NN, respectively; and updating the actor NN, the critic NN, and the mass NN using corresponding calculated NN errors, respectively.

18. The storage medium according to claim 17, wherein the method further includes: if the initialized error of the actor NN is less than the initialized error threshold of the actor NN, obtaining previous calculated NN weights of the actor NN; orif the initialized error of the critic NN is less than the initialized error threshold of the critic NN, obtaining previous calculated NN weights of the critic NN; orif the initialized error of the mass NN is less than the initialized error threshold of the mass NN, obtaining previous calculated NN weights of the mass NN.

19. The storage medium according to claim 18, wherein the method further includes: using the previous calculated NN weights of the actor NN to calculate a control; andexecuting the calculated control.

20. The storage medium according to claim 17, wherein before initializing the errors, the method further includes: initializing a state and a density of the agent, wherein the state of the agent includes a position and a velocity; andcalculating an error of the agent using the state of the agent and a predefined trajectory.