Distributed Algorithms for Linearly-Solvable Optimal Control in Networked Multi-Agent Systems

For a MAS consisting of $N\in\mathbb{N}$ agents, $\mathcal{D}=\{1,2,\cdots,N\}$ denotes the index set of agents, and the communication network among agents is described by an undirected graph $\mathcal{G}=\{\mathcal{V},\mathcal{E}\}$, which implies that the communication channel between any two agents is bilateral. The communication network $\mathcal{G}$ is assumed to be connected. An agent $i\in\mathcal{D}$ in the network is denoted as a vertex $v_{i}\in\mathcal{V}=\{v_{1},v_{2},\cdots,v_{N}\}$, and an undirected edge $(v_{i},v_{j})\in\mathcal{E}\subset\mathcal{V}\times\mathcal{V}$ in graph $\mathcal{G}$ implies that the agents $i$ and $j$ can measure the states of each other, which is denoted by $x_{i}=[x_{i(1)},x_{i(2)},\cdots,x_{i(M)}]^{\top}\in\mathbb{R}^{M}$ for agent $i\in\mathcal{D}$. Agents $i$ and $j$ are neighboring or adjacent agents, if there exists a communication channel between them, and the index set of agents neighboring to agent $i$ is denoted by $\mathcal{N}_{i}$ with $\bar{\mathcal{N}}_{i}=\mathcal{N}_{i}\cup\{i\}$. We use column vectors $\bar{x}_{i}=[x_{i}^{\top},x^{\top}_{j\in\mathcal{N}_{i}}]^{\top}\in\mathbb{R}^{M\cdot|\mathcal{\bar{N}}_{i}|}$ to denote the joint state of agent $i$ and its adjacent agents $\mathcal{N}_{i}$, which together group up the factorial subsystem $\bar{\mathcal{N}}_{i}$ of agent $i$, and $x=[x_{1}^{\top},x_{2}^{\top},\cdots,x^{\top}_{n}]^{\top}\in\mathbb{R}^{M\cdot N}$ denotes the global states of MAS. Figure 1 shows a MAS and all its factorial subsystems.

To optimize the correlations between neighboring agents while not intensifying the communication and computational complexities of the network or each agent, our paper proposes a distributed planning and execution framework, as a trade-off scheme between the easiness of implementation and optimality of performance, for multi-agent LSOC problems. Under this distributed framework, the local control action $u_{i}$ of agent $i\in N$ is computed by solving the local LSOC problem defined in factorial subsystem $\bar{\mathcal{N}}_{i}$. Instead of requiring the cost functions fully factorized over agents [31] or the knowledge of global states [30, 25], joint cost functions are permitted in every factorial subsystem, which captures the correlations and cooperation between neighboring agents, and the local control action $u_{i}$ only relies on the local observation $\bar{x}_{i}$ of agent $i$, which also simplifies the structure of communication network. Meanwhile, the global computational complexity is no longer exponential with respect to the total amount of agents in network $|\mathcal{D}|$, but becomes linear with respect to it in general, and the computational complexity of local agent $i$ is only related to the number of agents in factorial subsystem $\bar{\mathcal{N}}_{i}$. However, since the local control actions are computed from the local observations with only partial information, the distributed LSOC law obtained under this framework is usually a sub-optimal solution, despite that it coincides with the global optimal solution when the communication network is fully connected. More explanations and discussions on this will be given at the end of this section. Before that, we first reformulate the discrete-time and continuous-time LSOC problems from the single-agent scenario to a multi-agent setting that is compatible with our distributed framework.

Discrete-time SOC for single-agent systems, also known as the single-agent MDP, is briefly reviewed and then generalized to the networked MAS scenario. We consider the single-agent MDPs with finite state space and continuous control space in the first-exit setting. For a single agent $i\in\mathcal{D}$, the state variable $x_{i}$ belongs to a finite set $\mathcal{S}_{i}=\{s^{i}_{1},s^{i}_{2},\cdots\}=\mathcal{I}_{i}\cup\mathcal{B}_{i}$, which may be generated from an infinite-dimensional state problem by an appropriate coding scheme [34]. $\mathcal{I}_{i}\subset\mathcal{S}_{i}$ denotes the set of interior states of agent $i$, and $\mathcal{B}_{i}\subset\mathcal{S}_{i}$ denotes the set of boundary states. Without communication and interference from other agents, the passive dynamics of agent $i$ follows the probability distribution

$$x^{\prime}_{i}\sim p_{i}(\cdot|x_{i}),$$

where $x_{i},x^{\prime}_{i}\in\mathcal{S}_{i}$, and $p_{i}(x^{\prime}_{i}|x_{i})$ denotes the transition probability from state $x_{i}$ to $x^{\prime}_{i}$. When taking control action $u_{i}$ at state $x_{i}$, the controlled dynamics of agent $i$ is described by the distribution mapping

$$x_{i}^{\prime}\sim u_{i}(\cdot|x_{i})=p_{i}(\cdot|x_{i},u_{i}),$$

(1)

where $u_{i}(x^{\prime}_{i}|x_{i})$ or $p_{i}(x^{\prime}_{i}|x_{i},u_{i})$ denotes the transition probability from state $x_{i}$ to state $x^{\prime}_{i}$ subject to control $u_{i}$ that belongs to a continuous space. We require that $u_{i}(x^{\prime}_{i}|x_{i})=0$, whenever $p_{i}(x^{\prime}_{i}|x_{i})=0$, to prevent the direct transitions to the goal states. When $x_{i}\in\mathcal{I}_{i}\subset\mathcal{S}_{i}$, the immediate or running cost function of LSMDPs is designed as:

$$c_{i}(x_{i},u_{i})=q_{i}(x_{i})+\textrm{KL}(u_{i}(\cdot|x_{i})\parallel p_{i}(\cdot|x_{i})),$$

(2)

where the state cost $q_{i}(x_{i})$ can be an arbitrary function encoding how (un)desirable different states are, and the KL-divergence^***The KL-divergence (relative entropy) between two discrete probability mass functions $p(x)$ and $q(x)$ is defined as $$\textrm{KL}(p\parallel q)=\sum_{x\in\mathcal{X}}p(x)\log[{p(x)}/{q(x)}],$$ which has an absolute minimum $0$ when $p(x)=q(x),\forall x\in\mathcal{X}$. For two continuous probability density functions $p(x)$ and $q(x)$, the KL-divergence is defined as $$\textrm{KL}(p\parallel q)=\int_{x\in\chi}p(x)\log[p(x)/q(x)]dx.$$ measures the cost of control actions. When $x_{i}\in\mathcal{B}_{i}\subset\mathcal{S}_{i}$, the final cost function is defined as $\phi_{i}(x_{i})\geq 0$. The cost-to-go function of first-exit problem starting at state-time pair $(x_{i}^{t_{0}},t_{0})$ is defined as

$$J^{u_{i}}_{i}(x^{t_{0}}_{i},t_{0})=\mathbb{E}^{u_{i}}\bigg{[}\phi_{i}(x_{i}^{t_{f}})+\sum_{\tau=t_{0}}^{t_{f}-1}c_{i}(x_{i}^{\tau},u^{\tau}_{i})\bigg{]},$$

(3)

where $(x^{\tau}_{i},u_{i}^{\tau})$ is the state-action pair of agent $i$ at time step $\tau$, $x^{t_{f}}_{i}$ is the terminal or exit state, and the expectation $\mathbb{E}^{u_{i}}$ is taken with respect to the probability measure under which $x_{i}$ satisfies (1) given the control law $u_{i}=(u^{t_{0}}_{i},u^{t_{1}}_{i},\cdots,u^{t_{f}-1}_{i})$ and initial condition $x_{i}^{t_{0}}$. The objective of discrete-time stochastic optimal control problem is to find the optimal policy $u_{i}^{*}$ and value functions $V_{i}(x_{i})$ by solving the Bellman equation

$$V_{i}(x_{i})=\min_{u_{i}}\left\{c_{i}(x_{i},u_{i})+\mathbb{E}_{x^{\prime}_{i}\sim u_{i}(\cdot|x_{i})}[V_{i}(x^{\prime}_{i})]\right\},$$

(4)

where the value function $V_{i}(x_{i})$ is defined as the expected cumulative cost for starting at state $x_{i}$ and acting optimally thereafter, i.e. $V_{i}(x_{i})=\min_{u_{i}}J^{u_{i}}_{i}(x^{t_{0}}_{i},t_{0})$.

Based on the formulations of single-agent LSMDP and augmented dynamics, we introduce a multi-agent LSMDP formulation subject to the distributed framework of the factorial subsystem. For simplicity, we assume that the passive dynamics of agents in MAS are homogeneous and mutually independent, i.e. agents without control are governed by identical dynamics and do not interfere or collide with each other. This assumption is also posited in many previous papers on multi-agent LSOC or distributed control [35, 31, 11, 25]. Since agent $i\in\mathcal{D}$ can only observe the states of neighboring agents $\mathcal{N}_{i}$, we are interested in the subsystem $\bar{\mathcal{N}}_{i}=\mathcal{N}_{i}\cup\{i\}$ when computing the control law of agent $i$. Hence, the autonomous dynamics of subsystem $\bar{\mathcal{N}}_{i}$ follow the distribution mapping

$$\bar{x}^{\prime}_{i}\sim\bar{p}_{i}(\cdot|\bar{x}_{i})=\prod_{j\in{\mathcal{\bar{N}}}_{i}}p_{j}(\cdot|x_{j}),$$

(5)

where the joint state $\bar{x}_{i},\bar{x}^{\prime}_{i}\in\prod_{j\in\mathcal{\bar{N}}_{i}}S_{j}=\bar{\mathcal{S}}_{i}=\{\bar{s}^{i}_{1},\bar{s}^{i}_{2},\cdots\}=\bar{\mathcal{I}}_{i}\cup\bar{\mathcal{B}}_{i}$, and distribution information $\bar{p}_{i}(\cdot|\bar{x}_{i})$ are generally only accessible to agent $i$. Similarly, the global state of all agents $x^{\prime}\sim p(\cdot|x)=\prod_{i=1}^{N}p_{i}(\cdot|x_{i})$, which is usually not available to a local agent $i$ unless it can receive the information from all other agents $\mathcal{D}\backslash\{i\}$, e.g. agent 3 in Figure 1. Since the local control action $u_{i}$ only relies on the local observation of agent $i$, we assume that the joint posterior state $\bar{x}^{\prime}_{i}$ of subsystem $\bar{\mathcal{N}}_{i}$ is exclusively determined by the joint prior state $\bar{x}_{i}$ and joint control $\bar{u}_{i}$ when computing the optimal control actions in subsystem $\bar{\mathcal{N}}_{i}$. More intuitively, under this assumption, the local LSOC algorithm in subsystem $\bar{\mathcal{N}}_{i}$ only requires the measurement of joint state $\bar{x}_{i}$ and treats subsystem $\bar{\mathcal{N}}_{i}$ as a complete connected network, as shown in Figure 1. When each agent in $\mathcal{\bar{N}}_{i}$ samples their control action independently, the joint controlled dynamics of the factorial subsystem $\bar{\mathcal{N}}_{i}$ satisfies

$$\bar{x}^{\prime}_{i}\sim\bar{u}_{i}(\cdot|\bar{x}_{i})=\prod_{j\in\mathcal{\bar{N}}_{i}}u_{j}(\cdot|\bar{x}_{i})=\prod_{j\in\mathcal{\bar{N}}_{i}}p_{j}(\cdot|x_{i},x_{j\in\mathcal{N}_{i}},u_{j}),$$

(6)

where the joint state $\bar{x}_{i}$ and joint distribution $\bar{u}_{i}(\cdot|\bar{x}_{i})$ are only accessible to agent $i$ in general. Once we figure out the joint control distribution $\bar{u}_{i}(\cdot|\bar{x}_{i})$ for subsystem $\bar{\mathcal{N}}_{i}$, the local control distribution $u_{i}(\cdot|\bar{x}_{i})$ of agent $i$ can be retrieved by calculating the marginal distribution. The joint immediate cost function for subsystem $\bar{\mathcal{N}}_{i}$ when $\bar{x}_{i}\in\mathcal{\bar{I}}_{i}$ is defined as follows

$$c_{i}(\bar{x}_{i},\bar{u}_{i})=q_{i}(\bar{x}_{i})+\textrm{KL}(\bar{u}_{i}(\cdot|\bar{x}_{i})\parallel\bar{p}_{i}(\cdot|\bar{x}_{i}))=q_{i}(\bar{x}_{i})+\sum_{j\in\bar{\mathcal{N}}_{i}}\textrm{KL}(u_{j}(\cdot|\bar{x}_{i})\parallel p_{j}(\cdot|x_{j})),$$

(7)

where the state cost $q_{i}(\bar{x}_{i})$ can be an arbitrary function of joint state $\bar{x}_{i}$, i.e. a constant or the norm of disagreement vector, and the second equality follows from (5) and (6), which implies that the joint control cost is the cumulative sum of the local control costs. When $\bar{x}_{i}\in\mathcal{\bar{B}}_{i}$, the exit cost function $\phi_{i}(\bar{x}_{i})=\sum_{j\in\mathcal{\bar{N}}_{i}}\omega^{i}_{j}\cdot\phi_{j}(x_{j})$, where $\omega^{i}_{j}>0$ is a weight measuring the priority of assignment on agent $j$. In order to improve the success rate in application, it is preferable to assign $\omega_{i}^{i}$ as the largest weight when computing the control distribution $\bar{u}_{i}$ in subsystem $\mathcal{\bar{N}}_{i}$. Subsequently, the joint cost-to-go function of first-exit problem in subsystem $\bar{\mathcal{N}}_{i}$ becomes

$$J^{\bar{u}_{i}}_{i}(\bar{x}^{t_{0}}_{i},t_{0})=\mathbb{E}^{\bar{u}_{i}}\bigg{[}\phi_{i}(\bar{x}_{i}^{t_{f}})+\sum_{\tau=t_{0}}^{t_{f}-1}c_{i}(\bar{x}_{i}^{\tau},\bar{u}^{\tau}_{i})\bigg{]}.$$

(8)

Some abuses of notations occur when we define the cost functions $c_{i}$ and $\phi_{i}$, cost-to-go function $J_{i}$, and value function $V_{i}$ in single-agent setting and the factorial subsystem; one can differentiate the different settings from the arguments of these functions. Derived from the single-agent Bellman equation (4), the joint optimal control action $\bar{u}_{i}^{*}(\cdot|\bar{x}_{i})$ subject to the joint cost function (7) can be solved from the following joint Bellman equation in subsystem $\bar{\mathcal{N}}_{i}$

$$V_{i}(\bar{x}_{i})=\min_{\bar{u}_{i}}\left\{c_{i}(\bar{x}_{i},\bar{u}_{i})+\mathbb{E}_{\bar{x}^{\prime}_{i}\sim\bar{u}_{i}(\cdot|\bar{x}_{i})}[V_{i}(\bar{x}^{\prime}_{i})]\right\},$$

(9)

where $V_{i}(\bar{x}_{i})$ is the (joint) value function of joint state $\bar{x}_{i}$. A linearization method as well as a parallel programming method for solving (9) will be discussed in Section 3.

For continuous-time LSOC problems, we first consider the dynamics of single agent $i$ described by the following Itô diffusion process

$$dx_{i}=f_{i}(x_{i},t)dt+B_{i}(x_{i})[u_{i}(x_{i},t)dt+\sigma_{i}dw_{i}],$$

(10)

where $x_{i}\in\mathbb{R}^{M}$ is agent $i$’s state vector from an uncountable state space; $f_{i}(x_{i},t)+B_{i}(x_{i})\cdot u_{i}(x_{i},t)\in\mathbb{R}^{M}$ is the deterministic drift term with passive dynamics $f_{i}(x_{i},t)$, control matrix $B_{i}(x_{i})\in\mathbb{R}^{M\times P}$ and control action $u_{i}(x_{i},t)\in\mathbb{R}^{P}$; noise $dw_{i}\in\mathbb{R}^{P}$ is a vector of possibly correlated^†††When the components of $d{\tilde{w}}_{i}=[d{\tilde{w}}_{i,(1)},\cdots,d{\tilde{w}}_{i,(P)}]^{\top}$ are correlated and satisfy a multi-variate normal distribution $N(0,\Sigma_{i})$, by using the Cholesky decomposition $\Sigma_{i}=\sigma_{i}\sigma^{\top}_{i}$, we can rewrite $d\tilde{w}_{i}=\sigma_{i}dw_{i}$, where $dw_{i}$ is a vector of Brownian components with zero drift and unit-variance rate. Brownian components with zero mean and unit rate of variance, and the positive semi-definite matrix $\sigma_{i}\in\mathbb{R}^{P\times P}$ denotes the covariance of noise $dw_{i}$. When $x_{i}\in\mathcal{I}_{i}$, the running cost function is defined as

$$c_{i}(x_{i},u_{i})=q_{i}(x_{i})+\dfrac{1}{2}u_{i}(x_{i},t)^{\top}R_{i}u_{i}(x_{i},t),$$

(11)

where $q_{i}(x_{i})\geq 0$ is the state-related cost, and $u_{i}^{\top}R_{i}u_{i}$ is the control-quadratic term with matrix $R_{i}\in\mathbb{R}^{P\times P}$ being positive definite. When $x^{t_{f}}_{i}\in\mathcal{B}_{i}$, the terminal cost function is $\phi_{i}(x^{t_{f}}_{i})$, where $t_{f}$ is the exit time. Hence, the cost-to-go function of first-exit problem is defined as

$$J^{u_{i}}_{i}(x_{i}^{t},t)=\mathbb{E}^{u_{i}}_{x_{i}^{t},t}\left[\phi_{i}(x_{i}^{t_{f}})+\int_{t}^{t_{f}}c_{i}(x_{i}(\tau),u_{i}(\tau))\ d\tau\right],$$

(12)

where the expectation is taken with respect to the probability measure under which $x_{i}$ is the solution to (10) given the control law $u_{i}$ and initial condition $x_{i}(t)$. The value function is defined as the minimal cost-to-go function $V_{i}(x_{i},t)=\min_{u_{i}}J_{i}^{u_{i}}(x_{i}^{t},t)$. Subject to the dynamics (10) and running cost function (11), the optimal control action $u_{i}^{*}$ can be solved from the following single-agent stochastic Hamilton–Jacobi–Bellman (HJB) equation:

	$\displaystyle-\partial_{t}V_{i}(x_{i},t)=\min_{u_{i}}\Big{\{}c_{i}(x_{i},u_{i})+[f_{i}(x_{i},t)$	$\displaystyle+B_{i}(x_{i})u_{i}(x_{i},t)]^{\top}\cdot\nabla_{x_{i}}V_{i}(x_{i},t)$		(13)
		$\displaystyle+\frac{1}{2}\textrm{tr}\left[B_{i}(x_{i})\sigma_{i}\sigma^{\top}_{i}B_{i}(x_{i})^{\top}\cdot\nabla^{2}_{x_{i}x_{i}}V_{i}(x_{i},t)\right]\Big{\}},$

where $\nabla_{x_{i}}$ and $\nabla^{2}_{x_{i}x_{i}}$ respectively refer to the gradient and Hessian matrix with $\nabla_{x_{i}}V_{i}=[\partial V_{i}/\partial x_{i(1)},\break\cdots,\partial V_{i}/\partial x_{i(M)}]^{\top}$ and elements $[\nabla^{2}_{x_{i}x_{i}}V_{i}]_{m,n}=\partial^{2}V_{i}/\partial x_{i(m)}\partial x_{i(n)}$. A few methods have been proposed to solve the stochastic HJB in (13), such as the approximation methods via discrete-time MDPs or eigenfunction [36] and path integral approaches [37, 31, 22].

Similar to the extension of discrete-time LSMDP from single-agent setting to MAS, the joint continuous-time dynamics for factorial subsystem $\bar{\mathcal{N}}_{i}$ is described by

$$d\bar{x}_{i}=\bar{f}_{i}(\bar{x}_{i},t)dt+\bar{B}_{i}(\bar{x}_{i})\left[\bar{u}_{i}(\bar{x}_{i},t)dt+\bar{\sigma}_{i}d\bar{w}_{i}\right],$$

(14)

where the joint passive dynamics vector is denoted by $\bar{f}_{i}(\bar{x}_{i},t)=[f_{i}(x_{i},t)^{\top},f_{j\in\mathcal{N}_{i}}(x_{j},t)^{\top}]^{\top}\in\mathbb{R}^{M\cdot|\mathcal{\bar{N}}_{i}|}$, the joint control matrix is denoted by $\bar{B}_{i}(\bar{x}_{i})=\textrm{diag}\{B_{i}(x_{i}),B_{j\in\mathcal{N}_{i}}(x_{j})\}\in\break\mathbb{R}^{M\cdot|\bar{\mathcal{N}}_{i}|\times P\cdot|\bar{\mathcal{N}}_{i}|}$, $\bar{u}_{i}(\bar{x}_{i},t)=[u_{i}(\bar{x}_{i},t)^{\top},u_{j\in\mathcal{N}_{i}}(\bar{x}_{i},t)^{\top}]^{\top}\in\mathbb{R}^{P\cdot|\bar{\mathcal{N}}_{i}|}$ is the joint control action, $d\bar{w}_{i}=[dw^{\top}_{i},dw^{\top}_{j\in\mathcal{N}_{i}}]^{\top}\in\mathbb{R}^{P\cdot|\bar{\mathcal{N}}_{i}|}$ is the joint noise vector, and the joint covariance matrix is denoted by $\bar{\sigma}_{i}=\textrm{diag}\{\sigma_{i},\sigma_{j\in\mathcal{N}_{i}}\}\in\mathbb{R}^{P\cdot|\bar{\mathcal{N}}_{i}|\times P\cdot|\bar{\mathcal{N}}_{i}|}$. Analogous to the discrete-time scenario, we assume that the passive dynamics of agents are homogeneous and mutually independent, and for the local planning algorithm on agent $i$ or subsystem $\bar{\mathcal{N}}_{i}$, which computes the local control action $u_{i}(\bar{x}_{i},t)$ for agent $i$ and joint control action $\bar{u}_{i}(\bar{x}_{i},t)$ in subsystem $\bar{\mathcal{N}}_{i}$, the evolution of joint state $\bar{x}_{i}$ only depends on the current values of $\bar{x}_{i}$ and joint control $\bar{u}_{i}(\bar{x}_{i},t)$. When $\bar{x}_{i}\in\mathcal{\bar{I}}_{i}$, the joint immediate cost function for subsystem $\bar{\mathcal{N}}_{i}$ is defined as

$$c_{i}(\bar{x}_{i},\bar{u}_{i})=q_{i}(\bar{x}_{i})+\frac{1}{2}\bar{u}_{i}(\bar{x}_{i},t)^{\top}\bar{R}_{i}\bar{u}_{i}(\bar{x}_{i},t),$$

(15)

where the state-related cost $q_{i}(\bar{x}_{i})$ can be an arbitrary function measuring the (un)desirability of different joint states $\bar{x}_{i}\in\mathcal{\bar{S}}_{i}$, and $\bar{u}_{i}^{\top}\bar{R}_{i}\bar{u}_{i}$ is the control-quadratic term with matrix $\bar{R}_{i}\in\mathbb{R}^{P\cdot|\bar{\mathcal{N}}_{i}|\times P\cdot|\bar{\mathcal{N}}_{i}|}$ being positive definite. When $\bar{R}_{i}=\textrm{diag}\{R_{i},R_{j\in\mathcal{N}_{i}}\}$ with $R_{i}$ and $R_{j}$ defined in (11), the joint control cost term in (15) satisfies $\bar{u}_{i}^{\top}\bar{R}_{i}\bar{u}_{i}=\sum_{j\in\mathcal{\bar{N}}_{i}}u_{j}^{\top}R_{j}u_{j}$, which is symmetric with respect to the relationship of discrete-time control costs in (7). When $\bar{x}_{i}\in\bar{\mathcal{B}}_{i}$, the terminal cost function is defined as $\phi_{i}(\bar{x}_{i})=\sum_{j\in\mathcal{\bar{N}}_{i}}\omega_{j}^{i}\cdot\phi_{j}(x_{j})$, where $\omega_{j}^{i}>0$ is the weight measuring the priority of assignment on agent $j$, and we let the weight $\omega_{i}^{i}$ dominate other weights $\omega_{j\in\mathcal{N}_{i}}^{i}$ to improve the success rate. Compared with the cost functions fully factorized over agents, the joint cost functions in (15) can gauge and facilitate the correlation and cooperation between neighboring agents. Subsequently, the joint cost-to-go function of first-exit problem in subsystem $\bar{\mathcal{N}}_{i}$ is defined as

$$J^{\bar{u}_{i}}_{i}(\bar{x}_{i}^{t},t)=\mathbb{E}^{\bar{u}_{i}}_{\bar{x}_{i}^{t},t}\left[\phi_{i}(\bar{x}_{i}^{t_{f}})+\int_{t}^{t_{f}}c_{i}(\bar{x}_{i}(\tau),\bar{u}_{i}(\tau))\ d\tau\right].$$

Let the (joint) value function $V_{i}(\bar{x}_{i},t)$ be the minimal cost-to-go function, i.e. $V_{i}(\bar{x}_{i},t)=\min_{\bar{u}_{i}}J_{i}^{\bar{u}_{i}}(\bar{x}_{i}^{t},t)$. We can then compute the joint optimal control action $\bar{u}_{i}^{*}$ of subsystem $\mathcal{\bar{N}}_{i}$ by solving the following joint optimality equation

$$V_{i}(\bar{x}_{i},t)=\min_{\bar{u}_{i}}\mathbb{E}^{\bar{u}_{i}}_{\bar{x}_{i}^{t},t}\left[\phi_{i}(\bar{x}_{i}^{t_{f}})+\int_{t}^{t_{f}}c_{i}(\bar{x}_{i}(\tau),\bar{u}_{i}(\tau))\ d\tau\right].$$

(16)

A distributed path integral control algorithm and a distributed REPS algorithm for solving (16) will be respectively discussed in Section 3. Discussion on the relationship between discrete-time and continuous-time LSOC dynamics can be found in [16, 38].

We first consider the discrete-time MAS with dynamics (6) and immediate cost function (7), which give the joint Bellman equation (9) in subsystem $\mathcal{\bar{N}}_{i}$

$$V_{i}(\bar{x}_{i})=\min_{\bar{u}_{i}}\left\{c_{i}(\bar{x}_{i},\bar{u}_{i})+\mathbb{E}_{\bar{x}^{\prime}_{i}\sim\bar{u}_{i}(\cdot|\bar{x}_{i})}[V_{i}(\bar{x}^{\prime}_{i})]\right\}.$$

In order to compute the value function $V_{i}(\bar{x}_{i})$ and optimal control action $\bar{u}_{i}^{*}(\cdot|\bar{x}_{i})$ from equation (9), an exponential or Cole-Hopf transformation is employed to linearize (9) into a system of linear equations, which can be cast into a decentralized programming and solved in parallel. Local optimal control distribution $u_{i}^{*}(\cdot|\bar{x}_{i})$ that depends on the local observation of agent $i$ is then derived by marginalizing the joint control distribution $\bar{u}_{i}^{*}(\cdot|\bar{x}_{i})$.

We now consider the LSOC for continuous-time MASs subject to joint dynamics (14) and joint immediate cost function (15), which can be formulated into the joint optimality equation (16) as follows

$$V_{i}(\bar{x}_{i},t)=\min_{\bar{u}_{i}}\mathbb{E}^{\bar{u}_{i}}_{\bar{x}_{i}^{t},t}\left[\phi_{i}(\bar{x}_{i}^{t_{f}})+\int_{t}^{t_{f}}c_{i}(\bar{x}_{i}(\tau),\bar{u}_{i}(\tau))\ d\tau\right].$$

In order to solve (16) and derive the local optimal control action $u^{*}_{i}(\bar{x}_{i},t)$ for agent $i$, we first cast the joint optimality equation (16) into a joint stochastic HJB equation that gives an analytic form for joint optimal control action $\bar{u}^{*}_{i}(\bar{x}_{i},t)$. To solve for the value function, by resorting to the Cole-Hopf transformation, the stochastic HJB equation is linearized into a partial differential equation (PDE) with respect to desirability function. Feynman-Kac formula is then invoked to formulate the solution of the linearized PDE and joint optimal control action as the path integral formulae forward in time, which are later approximated respectively by a distributed Monte Carlo (MC) sampling method and a sample-efficient distributed REPS algorithm.