beginalgorithmic 2024 1

Pierer v. Esch et al. Sensitivity-Based Distributed Model Predictive Control
for Nonlinear Systems under Inexact Optimization

Maximilian Pierer von Esch, Chair of Automatic Control, Friedrich-Alexander-
Universität Erlangen-Nürnberg, 91058 Erlangen, Germany.

Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under
project no.Gr 3870/6-1.

# Sensitivity-Based Distributed Model Predictive Control for Nonlinear Systems
under Inexact Optimization

Maximilian Pierer von Esch  Andreas Völz  Knut Graichen  Chair of Automatic
Control, Friedrich-Alexander-Universität Erlangen-Nürnberg, Erlangen, Germany
<EMAIL_ADDRESS>M. Pierer von Esch  A. Völz  K. Graichen

(Date Month Year; Date Month Year; Date Month Year)

###### Abstract

[Abstract]This paper presents a distributed model predictive control (DMPC)
scheme for nonlinear continuous-time systems. The underlying distributed
optimal control problem is cooperatively solved in parallel via a sensitivity-
based algorithm. The algorithm is fully distributed in the sense that only one
neighbor-to-neighbor communication step per iteration is necessary and that
all computations are performed locally. Sufficient conditions are derived for
the algorithm to converge towards the central solution. Based on this result,
stability is shown for the suboptimal DMPC scheme under inexact minimization
with the sensitivity-based algorithm and verified with numerical simulations.
In particular, stability can be guaranteed with either a suitable stopping
criterion or a fixed number of algorithm iterations in each MPC sampling step
which allows for a real-time capable implementation.

, , Sensitivity-Based Distributed Model Predictive Control for Nonlinear
Systems under inexact optimization

###### keywords:

Distributed model predictive control, sensitivities, sensitivity-based
distributed optimal control, nonlinear systems, networked systems, stability,
suboptimality

††articletype: Research Article††journal: Journal††volume: 00††footnotetext:
Abbreviations: ADMM, alternating direction multiplier method; CLF, control
Lyapunov function; DMPC, distributed model predictive control; MPC, model
predictive control; OCP, optimal control problem

##  1 Introduction

Model predictive control (MPC) has emerged as a powerful strategy for
controlling both linear and nonlinear systems, demonstrating its versatility
in handling constraints and optimizing cost functions.1, 2, 3 However,
confronted with large-scale systems, a centralized MPC approach for the entire
system is often impractical or undesirable. This limitation may be due to the
unavailability of a centralized control entity with sufficient computing
resources, scalability concerns, or communication restrictions imposed by a
given communication topology. These restrictions in combination with a shift
towards networked and decentralized control architectures have led to the
concept of distributed MPC (DMPC) in which the central MPC controller is
replaced by local MPCs, which control the single subsystems of the global
system.4, 5, 6, 7 The subsystems together with the local MPCs then form a
network of so-called agents.

The development of a DMPC scheme heavily depends on the considered problem
setting since various sorts of couplings between subsystems and communication
topologies give rise to different requirements for the DMPC controller.6 For
instance, the subsystems can be dynamically coupled in the sense that the
dynamics of the individual subsystems depend on other subsystems’ states
and/or inputs. Such couplings typically arise from physically connected
subsystems found in large-scale processing plants 8, 9, 10, infrastructure
networks such as smart grids11, 12, 13 and water distribution networks14, 15,
16 or coupled mechanical systems, for example, found in cooperative payload
scenarios17, 18. A majority of the DMPC literature considers linear discrete-
time systems, in which the subsystems are assumed to be linearly coupled to
the state and/or control of their neighbors.19, 20, 21, 22, 23 In recent
years, however, also dynamically coupled nonlinear continuous-time systems
have been the focus of research. 24, 25, 26, 27, 28

In other settings, the global system is composed of a group of dynamically
independent subsystems whose dynamics only depend on their own states and
controls. Typically, this is the case in multi-agent networks in which the
agents need to solve different cooperative tasks like formation control29, 30,
platooning31, synchronization and consensus32, flocking33 or coverage
problems34. These cooperative tasks are then often formulated in terms of
coupled cost functions and/or constraints. The difficulty in this regard is
that the control task usually differs from the classical setpoint
stabilization scenario.

Finally, the available communication topology has a major influence on the
DMPC scheme since the communication topology in general does not need to
correspond to the coupling structure which results in different design
requirements. A large share of DMPC schemes considers either neighbor-to-
neighbor communication where the communication and coupling topology are
identical21, 26, 35, 20, 36, 37, star networks in which the agents communicate
only with a central node 38, 39 or system-wide communication which allows the
agents to communicate with all other agents or a subset thereof40.
Consequently, different combinations of these communication structures are
possible. Reducing communication is a vital aspect of DMPC, as previous
studies have reported that communication can be responsible for a major part
of the overall computation time 26, 41, 37, 42. Neighbor-to-neighbor
communication is favorable in that regard as it lowers the system-wide
communication load and provides a more flexible and scalable communication
structure without a single point of failure.

In addition to different problem settings and communication requirements, the
inherent characteristics of the DMPC controller can vary drastically in how
neighboring subsystems are considered. Since a comprehensive and comparative
overview of different DMPC structures is out of scope of this paper, we refer
the reader to the survey papers 4, 6, 43. However, one of the most promising
approaches is distributed optimization, where the central MPC problem is
solved iteratively and in parallel by each agent. The idea is that after
enough algorithm iterations in each MPC step, the local solutions will
converge to the centralized one leading to the same performance as the
centralized MPC scheme which is known as cooperative DMPC.43

Distributed optimization algorithms require suitable decomposition schemes
such that the central problem can be decomposed into smaller subproblems which
are then solved locally and in parallel by the individual agents. For an
introductory overview of different optimization methods used in DMPC, we refer
to the surveys 44, 45. A majority of the distributed optimization methods in
DMPC rely on a decomposition of the dual problem, where coupling constraints
are taken into account by Lagrangian multipliers. Notable algorithms that fall
into this category are dual decomposition22, the alternating direction method
of multipliers (ADMM) 46, decentralized interior point methods 40, the
augmented Lagrangian-based alternating direction inexact Newton method 47 and
its bi-level distributed variant48. A large share of DMPC schemes employ these
dual decomposition-based optimization methods.49, 24, 26, 25, 22, 30, 23, 50,
19 However, the disadvantage of dual optimization methods such as ADMM lies in
the fact that it achieves primal feasibility only in the limit, strong duality
is required to hold, and slow convergence for highly coupled systems has been
observed in practice.37, 36

Another class of decomposition methods concerns primal decomposition. These
approaches are characterized by directly exchanging required quantities
between neighboring subsystems rather than using Lagrangian multipliers to
append the coupling constraints.51, 13, 8, 52, 53, 54 A comparison of primal
and dual decomposition methods in DMPC is for example given in references 44,
45, 55, 56. An advantage of primal approaches is that strong duality is not
required to hold and that the iterates are primal feasible under certain
conditions. However, it needs to be ensured that the distributed solution
corresponds with the centralized one. Most notable, sensitivity-based
approaches have recently emerged as viable distributed optimization algorithms
for DMPC.8, 57, 58, 35, 59, 37 All these works follow the same ideas which can
be traced back to the 1970s 60, 61. In particular, sensitivity-based
approaches have been utilized for automatic building control 62, 59, water
distribution networks 58, or process control 8. Furthermore, the
contributions37, 63 show that the sensitivities can be calculated locally by
each agent in a computationally efficient manner using optimal control theory
resulting in a fully distributed algorithm with only neighbor-to-neighbor
communication which overcomes the major drawback that the agents usually
require knowledge of the overall system dynamics, additional communication
steps or a central coordination step. In addition, convergence of the
sensitivity-based algorithm in the context of DMPC has been investigated for
the linear discrete-time case 8, an algorithm variant considering the complete
system dynamic at agent level 59 and most recently for nonlinear continuous-
time systems in a fully distributed setting 63 with an inexact optimization of
the underlying OCP. This makes the sensitivity-based approach an interesting
option for nonlinear DMPC.

Common approaches to ensure the stability of DMPC schemes utilizing
distributed optimization algorithms, especially for dynamically coupled
systems, are to use classical MPC terminal ingredients such as local terminal
costs and set constraints which can be extended to the distributed setting.64,
65, 20, 66 Compared to the central MPC case, synthesis of the terminal costs
and sets is non-trivial in the distributed case and usually requires
optimization-based approaches 20, 66, although iterative algorithms exist 67.
However, terminal set constraints that are designed for stability requirements
are often restrictive and unfavorable from a numerical viewpoint as they lead
to an increased computational load compared to an MPC formulation without
terminal constraints 68 which extends to DMPC 24. Another approach is relaxed
dynamic programming, in which the terminal conditions are dropped and
stability is guaranteed via a relaxed descent condition 69 which has found
application in DMPC.49 An important aspect of a (D)MPC scheme is to ensure
stability even after a finite number of optimization iterations68, 24, often
referred to as stability under inexact optimization 44, as this limits the
time-consuming communication steps. Typically, a suitable stopping criterion
or a lower bound on the required iterations is used to guarantee stability.
70, 24, 19, 39, 71

This contribution presents a fully distributed DMPC scheme for nonlinear
continuous-time systems that are coupled in dynamics and/or cost functions via
the states of their respective neighbors. The scheme is fully distributed in
the sense that only neighbor-to-neighbor communication and local computations
are required. The centralized MPC problem is formulated without terminal set
constraints which allows for a computationally efficient solution.72, 68
Stability is ensured via terminal cost functions which act as a control-
Lyaponuv function (CLF) and are synthesized offline via an optimization-based
approach such that a separable structure is obtained. The DMPC scheme is
realized via distributed optimization with a sensitivity-based primal
decomposition approach. In order to limit the communication between
subsystems, the algorithm is terminated either after a suitable stopping
criterion or a fixed number of iterations are reached. Prematurely stopping
the algorithm leads to a suboptimal control solution and subsequently a
mismatch between the optimal centralized and distributed solution. Although
related schemes have been investigated in 24, 70, 19, the approach presented
in this paper requires far less assumptions and considers nonlinear systems.
In particular, by establishing linear convergence of the sensitivity-based
algorithm, exponential stability and incremental improvement of the DMPC
scheme are shown if either the stopping criterion is sufficiently small or if
a certain number of fixed algorithm iterations are performed in each MPC step.

The paper is structured as follows: In Section 2 the problem statement and
system class are introduced for which the central optimal control problem
(OCP) is stated. Subsequently, the optimal MPC and DMPC control strategy is
reviewed and it is shown how separable terminal costs in the continuous-time
nonlinear case can be synthesized. Following the central considerations, the
distributed solution via the sensitivity-based approach is discussed in
Section 3. In particular, the notion of sensitivities is derived via classical
optimal control theory, and the sensitivity-based DMPC scheme is presented.
Sufficient conditions are derived such that the algorithm converges towards
the central optimal solution. The stability of the closed-loop system under
the DMPC control law is analyzed in Section 4. The algorithm is evaluated via
numerical simulations and compared to ADMM in Section 5 before the paper is
summarized in Section 6.

Several conventions are used throughout this text. The standard 2-norm
$\|\boldsymbol{q}\|:=\|\boldsymbol{q}\|_{2}=(|q_{1}|^{2}+...+|q_{n}|^{2})^{\frac{1}{2}}$
is used for vectors $\boldsymbol{q}\in\mathds{R}^{n}$, while for time
functions $\boldsymbol{q}:[0,T]\rightarrow\mathds{R}^{n}$, the vector-valued
$L_{\infty}$-norm with
$\|\boldsymbol{q}\|_{L_{\infty}}=\max_{t\in[0,T]}\|\boldsymbol{q}(t)\|$ is
utilized along with the corresponding function space $L_{\infty}(0,T;P)$ such
that $\boldsymbol{q}\in L_{\infty}(0,T;P)$ implies
$\|\boldsymbol{q}\|_{L_{\infty}}<\infty$ on $t\in[0,T]$ and
${\boldsymbol{q}(t)\in P\subseteq\mathds{R}^{n},\,t\in[0,T]}$. An
$r$-neighborhood of a point $\boldsymbol{v}_{0}\in\mathds{R}^{v}$ is denoted
as
$\mathcal{B}(\boldsymbol{v}_{0},r):=\\{\boldsymbol{v}\in\mathds{R}^{v}\,|\,\|\boldsymbol{v}-\boldsymbol{v}_{0}\|\leq
r\\}$, while an $r$-neighborhood to a set $\mathcal{S}\subset\mathds{R}^{v}$
is defined as
$\mathcal{S}^{r}:=\bigcup_{\boldsymbol{s}_{0}\in\mathcal{S}}\mathcal{B}(\boldsymbol{s}_{0},r)$.
The stacking and reordering of individual vectors
$\boldsymbol{v}_{i}\in\mathds{R}^{v_{i}},\,i\in\mathcal{V}$ from a set
$\mathcal{V}$ is defined as $[\boldsymbol{v}_{i}]_{i\in\mathcal{V}}$. The
partial derivative of a function $f(\boldsymbol{x},\boldsymbol{y})$ w.r.t. to
one of its arguments $\boldsymbol{x}$ is denoted as
$\partial_{\boldsymbol{x}}\,f(\boldsymbol{x},\boldsymbol{y})$. For given
iterates $(\boldsymbol{x}^{k},\boldsymbol{y}^{k})$ at step $k$ of an arbitrary
algorithm, the short-hand notation
$\partial_{\boldsymbol{x}}\,f^{k}=\partial_{\boldsymbol{x}}\,f(\boldsymbol{x},\boldsymbol{y})\big{|}_{\boldsymbol{x}=\boldsymbol{x}^{k},\boldsymbol{y}=\boldsymbol{y}^{k}}$
is used when applicable. The argument of the time functions is omitted in the
paper when it is convenient. Finally, system variables are underlined (e.g.
$\underaccent{\bar}{x}$) to distinguish them from MPC-internal variables.

##  2 Problem statement and (D)MPC strategy

The structure of multi-agent systems is conveniently described by a graph
$\mathcal{G}=(\mathcal{V},\mathcal{E})$ in which the vertices
$\mathcal{V}=\\{1,\ldots,N\\}$ represent single dynamic subsystems referred to
as agents and the edges $\mathcal{E}\subset\mathcal{V}\times\mathcal{V}$
reflect a coupling between two agents. In this paper, the dynamics of an agent
$i\in\mathcal{V}$ with states
$\boldsymbol{\underaccent{\bar}{x}}_{i}\in\mathds{R}^{n_{i}}$ and controls
$\boldsymbol{\underaccent{\bar}{u}}_{i}\in\mathds{R}^{m_{i}}$ are described by
the following nonlinear neighbor-affine system

$\displaystyle\boldsymbol{\dot{\underaccent{\bar}{x}}}_{i}$
$\displaystyle=\boldsymbol{f}_{ii}(\boldsymbol{\underaccent{\bar}{x}}_{i},\boldsymbol{\underaccent{\bar}{u}}_{i})+\sum_{j\in\mathcal{N}_{i}^{\leftarrow}}\boldsymbol{f}_{ij}(\boldsymbol{\underaccent{\bar}{x}}_{i},\boldsymbol{\underaccent{\bar}{x}}_{j})=:\boldsymbol{f}_{i}(\boldsymbol{\underaccent{\bar}{x}}_{i},\boldsymbol{\underaccent{\bar}{u}}_{i},\boldsymbol{\underaccent{\bar}{x}}_{\mathcal{N}_{i}^{\leftarrow}})\,,\quad\boldsymbol{\underaccent{\bar}{x}}_{i}(0)=\boldsymbol{\underaccent{\bar}{x}}_{i,0}\,.$
(1)

The coupling to the neighbors is given by the states
$\boldsymbol{\underaccent{\bar}{x}}_{j}\in\mathds{R}^{n_{j}}$ using the
stacked notation
$\boldsymbol{\underaccent{\bar}{x}}_{\mathcal{N}_{i}^{\leftarrow}}:=[\boldsymbol{\underaccent{\bar}{x}}_{j}]_{j\in\mathcal{N}_{i}^{\leftarrow}}\in\mathds{R}^{p_{i}}$
with $p_{i}=\sum_{j\in\mathcal{N}_{i}^{\leftarrow}}n_{j}$ and
$p=\sum_{i\in\mathcal{V}}p_{i}$ where the set of sending neighbors
$\mathcal{N}^{\leftarrow}_{i}=\\{j\in\mathcal{V}\,:\,(j,i)\in\mathcal{E},\,j\neq
i\\}$ describes all the agents influencing agent $i\in\mathcal{V}$ while the
set of receiving neighbors
$\mathcal{N}_{i}^{\rightarrow}=\\{j\in\mathcal{V}\,:\,(i,j)\in\mathcal{E},\,i\neq
j\\}$ refers to the agents being influenced by agent $i\in\mathcal{V}$. The
union of both sets defines the neighborhood
$\mathcal{N}_{i}=\mathcal{N}^{\leftarrow}_{i}\cup\mathcal{N}_{i}^{\rightarrow}$.
In addition, it is allowed that the agents are able to communicate bi-
directionally with all neighbors $j\in\mathcal{N}_{i}$. Neighbor-affinity
means that in addition to the local dynamic functions
$\boldsymbol{f}_{ii}:\mathds{R}^{n_{i}}\times\mathds{R}^{m_{i}}\rightarrow\mathds{R}^{n_{i}}$,
the neighbors’ dynamic functions
$\boldsymbol{f}_{ij}:\mathds{R}^{n_{i}}\times\mathds{R}^{n_{j}}\rightarrow\mathds{R}^{n_{i}}$
enter additively into the dynamics (1) and only depend on exactly one other
state $\boldsymbol{\underaccent{\bar}{x}}_{j}$,
$j\in\mathcal{N}_{i}^{\leftarrow}$ .26 The individual controls
$\boldsymbol{\underaccent{\bar}{u}}_{i}$ are constrained to the compact and
convex constraint sets
$\boldsymbol{\underaccent{\bar}{u}}_{i}(t)\in\mathbb{U}_{i}\subset\mathds{R}^{m_{i}}$,
$t\geq 0$ which contain the origin $\boldsymbol{0}\in\mathds{R}^{m_{i}}$ in
their interior. The coupled subsystems (1) can equivalently be written in a
centralized form

$\displaystyle\boldsymbol{\dot{\underaccent{\bar}{x}}}=\boldsymbol{f}(\boldsymbol{\underaccent{\bar}{x}},\boldsymbol{\underaccent{\bar}{u}})\,,\quad\boldsymbol{\underaccent{\bar}{x}}(0)=\boldsymbol{\underaccent{\bar}{x}}_{0}$
(2)

with the central dynamics
$\boldsymbol{f}=[\boldsymbol{f}_{i}]_{i\in\mathcal{V}}$ as well as stacked
state
$\boldsymbol{\underaccent{\bar}{x}}=[\boldsymbol{\underaccent{\bar}{x}}_{i}]_{i\in\mathcal{V}}\in\mathds{R}^{n}$
and control
$\boldsymbol{\underaccent{\bar}{u}}=[\boldsymbol{\underaccent{\bar}{u}}_{i}]_{i\in\mathcal{V}}\in\mathds{R}^{m}$
vectors of dimension $n=\sum_{i\in\mathcal{V}}n_{i}$ and
$m=\sum_{i\in\mathcal{V}}m_{i}$, respectively. The stacked initial state is
given by
$\boldsymbol{\underaccent{\bar}{x}}_{0}=[\boldsymbol{\underaccent{\bar}{x}}_{i,0}]_{i\in\mathcal{V}}\in\mathds{R}^{n}$.
The local input constraints are summarized as
$\boldsymbol{\underaccent{\bar}{u}}(t)\in\mathbb{U}\subset\mathds{R}^{m}$ with
the set $\mathbb{U}$ defined as the Cartesian product
$\mathbb{U}=\prod_{i\in\mathcal{V}}\mathbb{U}_{i}$. Furthermore, it is assumed
that the local dynamics (1) and consequently the central dynamics (2) as well
as cost functions (4) are at least twice continuously differentiable w.r.t.
their arguments. In addition, system (2) must yield a bounded solution for any
initial condition $\boldsymbol{\underaccent{\bar}{x}}_{0}\in\mathds{R}^{n}$
and input $\boldsymbol{\underaccent{\bar}{u}}(t)\in\mathbb{U}$, $t\in[0,T]$
for some $0<T<\infty$. Without loss of generality, the control tasks consist
in controlling each subsystem (1) to the origin, i.e.
$\boldsymbol{f}_{i}(\boldsymbol{0},\boldsymbol{0},\boldsymbol{0})=\boldsymbol{0}$
or equivalently $\boldsymbol{f}(\boldsymbol{0},\boldsymbol{0})=\boldsymbol{0}$
in the centralized form (2).

###  2.1 Central optimal control problem and MPC stability

In this section, well known centralized MPC results are summarized which will
form the basis for the following investigations in the distributed case. 1,
68, 73, 24 In particular, the central MPC problem relies on the repeated
online solution of the following OCP at each sampling point $t_{k}=k\Delta t$,
$k\in\mathds{N}_{0}$

$\displaystyle\min_{\boldsymbol{u}}\quad$ $\displaystyle
J(\boldsymbol{u};\boldsymbol{\underaccent{\bar}{x}}_{k})=V(\boldsymbol{x}(T))+\int_{0}^{T}l(\boldsymbol{x}(\tau),\boldsymbol{u}(\tau))\,\mathrm{d}\tau$
(3a) $\displaystyle\leavevmode\nobreak\ \operatorname{s.t.}\quad$
$\displaystyle\boldsymbol{\dot{x}}=\boldsymbol{f}(\boldsymbol{x},\boldsymbol{u})\,,\quad\boldsymbol{x}(0)=\boldsymbol{\underaccent{\bar}{x}}_{k}$
(3b) $\displaystyle\boldsymbol{u}(\tau)\in\mathbb{U}\,,\quad\tau\in[0,\,T]$
(3c)

where
$\boldsymbol{\underaccent{\bar}{x}}(t_{k})=\boldsymbol{\underaccent{\bar}{x}}_{k}=[\boldsymbol{\underaccent{\bar}{x}}_{i,k}]_{i\in\mathcal{V}}$
is the state of system (2) at time $t=t_{k}$ and $\Delta t>0$ is the sampling
time. The cost function (3a) with horizon length $T>0$ consists of the
separable terminal and integral costs

$\displaystyle
V(\boldsymbol{x}(T))=\sum_{i\in\mathcal{V}}V_{i}(\boldsymbol{x}_{i}(T))\,,\quad
l(\boldsymbol{x},\boldsymbol{u})=\sum_{i\in\mathcal{V}}l_{i}(\boldsymbol{x}_{i},\boldsymbol{u}_{i},\boldsymbol{x}_{\mathcal{N}_{i}^{\leftarrow}})$
(4)

with $V_{i}:\mathds{R}^{n_{i}}\rightarrow\mathds{R}^{+}_{0}$ and
$l_{i}:\mathds{R}^{n_{i}}\times\mathds{R}^{m_{i}}\times\prod_{j\in\mathcal{N}_{i}}\mathds{R}^{n_{j}}\rightarrow\mathds{R}^{+}_{0}$.
Similar to (1), the coupled integral cost function (4) is given in the
neighbor-affine form

$l_{i}(\boldsymbol{x}_{i},\boldsymbol{u}_{i},\boldsymbol{x}_{\mathcal{N}_{i}^{\leftarrow}}):=l_{ii}(\boldsymbol{x}_{i},\boldsymbol{u}_{i})+\sum_{j\in\mathcal{N}_{i}^{\leftarrow}}l_{ij}(\boldsymbol{x}_{i},\boldsymbol{x}_{j})$
(5)

with the local part
$l_{ii}:\mathds{R}^{n_{i}}\times\mathds{R}^{m_{i}}\rightarrow\mathds{R}_{0}^{+}$
and the coupled part
$l_{ij}:\mathds{R}^{n_{i}}\times\mathds{R}^{n_{j}}\rightarrow\mathds{R}_{0}$.
The terminal and integral cost additionally satisfy the quadratic bounds

$\displaystyle m_{l}(\|\boldsymbol{x}\|^{2}+\|\boldsymbol{u}\|^{2})\leq
l(\boldsymbol{x},\boldsymbol{u})\leq
M_{l}(\|\boldsymbol{x}\|^{2}+\|\boldsymbol{u}\|^{2})\,,\quad
m_{V}\|\boldsymbol{x}\|^{2}\leq V(\boldsymbol{x})\leq
M_{V}\|\boldsymbol{x}\|^{2}$ (6)

for some constants $M_{l}\geq m_{l}>0$ and $M_{V}\geq m_{V}>0$. The admissible
input space to OCP (3) then follows as
$\mathcal{U}=L_{\infty}(0,T;\mathbb{U})$. Furthermore, there exists a non-
empty and open set $\Gamma\subset\mathds{R}^{n}$ with
$\boldsymbol{0}\in\Gamma$ such that for all
$\boldsymbol{\underaccent{\bar}{x}}_{k}\in\Gamma$, a minimizing and unique
solution $(\boldsymbol{x}^{*}(\tau;\boldsymbol{\underaccent{\bar}{x}}_{k})$,
$\boldsymbol{u}^{*}(\tau;\boldsymbol{\underaccent{\bar}{x}}_{k}))$,
$\tau\in[0,\,T]$ of (3) exists. The existence of an optimal solution is not
too restrictive as no terminal constraints are considered and all functions
are assumed to be continuously differentiable. Throughout this paper, the
central necessary optimality conditions to OCP (3) are needed. To this end,
define the central Hamiltonian as

$H(\boldsymbol{x},\boldsymbol{u},\boldsymbol{\lambda})=l(\boldsymbol{x},\boldsymbol{u})+\boldsymbol{\lambda}^{\mathrm{T}}\boldsymbol{f}(\boldsymbol{x},\boldsymbol{u})$
(7)

with the adjoint state $\boldsymbol{\lambda}\in\mathds{R}^{n}$. Then, the
first-order optimality conditions require that there exist optimal adjoint
states
$\boldsymbol{\lambda}^{*}(\tau;\boldsymbol{\underaccent{\bar}{x}}_{k})$,
$\tau\in[0,\,T]$ such that the canonical boundary value problem

$\displaystyle\boldsymbol{\dot{x}}^{*}(\tau;\boldsymbol{\underaccent{\bar}{x}}_{k})$
$\displaystyle=\boldsymbol{f}(\boldsymbol{x}^{*}(\tau;\boldsymbol{\underaccent{\bar}{x}}_{k}),\boldsymbol{u}^{*}(\tau;\boldsymbol{\underaccent{\bar}{x}}_{k}))\,,$
$\displaystyle\quad\boldsymbol{x}^{*}(0;\boldsymbol{\underaccent{\bar}{x}}_{k})=\boldsymbol{\underaccent{\bar}{x}}_{k}$
(8a)
$\displaystyle\boldsymbol{\dot{\lambda}}^{*}(\tau;\boldsymbol{\underaccent{\bar}{x}}_{k})$
$\displaystyle=-\partial_{\boldsymbol{x}}H(\boldsymbol{x}^{*}(\tau;\boldsymbol{\underaccent{\bar}{x}}_{k}),\boldsymbol{u}^{*}(\tau;\boldsymbol{\underaccent{\bar}{x}}_{k}),\boldsymbol{\lambda}^{*}(\tau;\boldsymbol{\underaccent{\bar}{x}}_{k}))=:\boldsymbol{G}(\boldsymbol{x}^{*},\boldsymbol{u}^{*},\boldsymbol{\lambda}^{*})\,,$
$\displaystyle\quad\boldsymbol{\lambda}^{*}(T;\boldsymbol{\underaccent{\bar}{x}}_{k})=\partial_{\boldsymbol{x}}V(\boldsymbol{x}^{*}(T;\boldsymbol{\underaccent{\bar}{x}}_{k}))$
(8b) is satisfied and that
$\boldsymbol{u}^{*}(\tau;\boldsymbol{\underaccent{\bar}{x}}_{k})$ minimizes
the Hamiltonian (7)
$\displaystyle\min_{\boldsymbol{u}\in\mathbb{U}}H(\boldsymbol{x}^{*}(\tau;\boldsymbol{\underaccent{\bar}{x}}_{k}),\boldsymbol{u},\boldsymbol{\lambda}^{*}(\tau;\boldsymbol{\underaccent{\bar}{x}}_{k}))\,,\quad\tau\in[0,\,T]\,.$
(8c)

The corresponding optimal (minimal) cost is denoted as

$\displaystyle
J^{*}(\boldsymbol{\underaccent{\bar}{x}}_{k}):=J(\boldsymbol{u}^{*}(\cdot,\boldsymbol{\underaccent{\bar}{x}}_{k});\boldsymbol{\underaccent{\bar}{x}}_{k})=\sum_{i\in\mathcal{V}}J_{i}(\boldsymbol{u}_{i}^{*}(\cdot;\boldsymbol{\underaccent{\bar}{x}}_{k});\boldsymbol{\underaccent{\bar}{x}}_{k})\,,$
(9)

where
$J_{i}(\boldsymbol{u}_{i}^{*}(\cdot;\boldsymbol{\underaccent{\bar}{x}}_{k});\boldsymbol{\underaccent{\bar}{x}}_{k})=:J_{i}^{*}(\boldsymbol{\underaccent{\bar}{x}}_{k})$
are the optimal agent costs. MPC strategies usually assume that this optimal
solution to OCP (3) is exactly known at each sampling point $t_{k}$. Then, the
first part of the optimal control trajectory on the sampling interval
$[t_{k},t_{k+1})$ is applied to the centralized system (2) which can be
interpreted as a nonlinear control law of the form

$\displaystyle\boldsymbol{\underaccent{\bar}{u}}(t_{k}+\tau)=\boldsymbol{u}^{*}(\tau;\boldsymbol{\underaccent{\bar}{x}}_{k})=:\boldsymbol{\kappa}^{*}(\boldsymbol{x}^{*}(\tau;\boldsymbol{\underaccent{\bar}{x}}_{k});\boldsymbol{\underaccent{\bar}{x}}_{k})\,,\quad\tau\in[0,\,\Delta
t)$ (10)

with sampling time $\Delta t<T$ and
$\boldsymbol{\kappa}^{*}(\boldsymbol{0};\boldsymbol{\underaccent{\bar}{x}}_{k})=\boldsymbol{0}$.
In the the next MPC time step $t_{k+1}$, the process of solving (3) is
repeated again with
$\boldsymbol{x}(0)=\boldsymbol{\underaccent{\bar}{x}}_{k+1}$ that (in the
nominal case) is given by
$\boldsymbol{\underaccent{\bar}{x}}_{k+1}=\boldsymbol{x}^{*}(\Delta
t;\boldsymbol{\underaccent{\bar}{x}}_{k})$. Since the central MPC problem is
formulated without terminal constraints, it is often assumed that the terminal
cost $V(\boldsymbol{x})$ in (6) represents a local CLF on an invariant set
$\Omega_{\beta}\subset\Gamma$ containing the origin at its center.68, 74 There
exists a feedback law
$\boldsymbol{u}=\boldsymbol{r}(\boldsymbol{x})\in\mathbb{U}$ and a non-empty
compact set
$\Omega_{\beta}:=\\{\boldsymbol{x}\in\Gamma:V(\boldsymbol{x})\leq\beta\\}\subset\Gamma$
containing the origin at its center such that for all
$\boldsymbol{x}\in\Omega_{\beta}$ the CLF inequality

$\dot{V}(\boldsymbol{x},\boldsymbol{r}(\boldsymbol{x}))+l(\boldsymbol{x},\boldsymbol{r}(\boldsymbol{x}))\leq
0$ (11)

with $\dot{V}(\boldsymbol{x},\boldsymbol{r}(\boldsymbol{x}))=\frac{\partial
V}{\partial\boldsymbol{x}}\boldsymbol{f}(\boldsymbol{x},\boldsymbol{r}(\boldsymbol{x}))$
is satisfied. A classical approach in MPC is to choose $V(\boldsymbol{x})$ as
the quadratic function
$V(\boldsymbol{x})=\boldsymbol{x}^{\mathrm{T}}\boldsymbol{P}\boldsymbol{x}$
where $\boldsymbol{P}\succ 0$ follows from the solution of a Lyapunov or
Riccati equation, given that the system (3) is stabilizable around the origin.
This results in a linear feedback law
$\boldsymbol{r}(\boldsymbol{x})=\boldsymbol{K}\boldsymbol{x}$ which stabilizes
the nonlinear system (2) on the (possibly small) invariant set
$\Omega_{\beta}$.75, 76, 77 However, due to the structural constraint (4)
concerning the terminal cost $V(\boldsymbol{x})$ and the requirement that only
direct neighbors are allowed to communicate, the design is more involved in
the distributed setting.20, 66 Therefore, a simple optimization-based approach
to design $V(\boldsymbol{x})$ which considers the separability constraint (4)
and to synthesize a structured control law, i.e.
$\boldsymbol{r}(\boldsymbol{x})=[\boldsymbol{r}_{i}(\boldsymbol{x}_{i},\boldsymbol{x}_{\mathcal{N}_{i}})]_{i\in\mathcal{V}}$
with
$\boldsymbol{x}_{\mathcal{N}_{i}}:=[\boldsymbol{x}_{j}]_{j\in\mathcal{N}_{i}}$,
will be discussed in Section 2.2. Throughout this paper, the following compact
level set of the optimal cost

$\displaystyle\Gamma_{\alpha}=\\{\boldsymbol{x}\in\Gamma:J^{*}(\boldsymbol{x})\leq\alpha\\}\,,\quad\alpha:=\beta\big{(}1+\frac{m_{l}}{m_{V}}T\big{)}\,,$
(12)

which characterizes the domain of attraction of the MPC scheme, will be
needed. Based on the preceding assumptions, the following stability results
for the MPC scheme without terminal constraints under the centralized control
law (10) can be shown.68, 78

######  Theorem 2.1.

Suppose that the Assumption 2.1 is satisfied. Then, for
$\boldsymbol{\underaccent{\bar}{x}}_{0}\in\Gamma_{\alpha}$ it holds that
$\boldsymbol{\underaccent{\bar}{x}}_{k}\in\Gamma_{\alpha}$ for all MPC steps
and
$\boldsymbol{x}^{*}(T;\boldsymbol{\underaccent{\bar}{x}}_{k})\in\Omega_{\beta}$
with $\Omega_{\beta}\subset\Gamma_{\alpha}$. Furthermore, the optimal cost in
the next MPC step decreases according to

$J^{*}(\boldsymbol{x}^{*}(\Delta
t;\boldsymbol{\underaccent{\bar}{x}}_{k}))\leq
J^{*}(\boldsymbol{\underaccent{\bar}{x}}_{k})-\int_{0}^{\Delta
t}l\big{(}\boldsymbol{x}^{*}(\tau;\boldsymbol{\underaccent{\bar}{x}}_{k}),\boldsymbol{u}^{*}(\tau;\boldsymbol{\underaccent{\bar}{x}}_{k})\big{)}\,\mathrm{d}\tau$
(13)

for all $\boldsymbol{\underaccent{\bar}{x}}_{k}\in\Gamma_{\alpha}$ and the
origin of the system (2) under the centralized MPC control law (10) is
asymptotically stable in the sense that the closed-loop trajectories satisfy
$\lim_{t\rightarrow\infty}\|\boldsymbol{\underaccent{\bar}{x}}(t)\|=0$.

Theorem 2.1 and the corresponding proofs can be found in references68, 73, 78.
Asymptotic stability of the closed-loop trajectories follows from (13) using
Barbalat’s Lemma and can be strengthened to exponential stability if
additionally the optimal cost (9) is continuously differentiable. Theorem 2.1
will serve as the basis for further stability considerations in the
distributed case.

###  2.2 Terminal cost design for distributed MPC

In this section it is discussed how the central MPC stability considerations
of the previous section can be transferred to the distributed case and in
particular how the separable terminal costs
$V(\boldsymbol{x})=\sum_{i\in\mathcal{V}}V_{i}(\boldsymbol{x}_{i})$ and
structured terminal controllers
$\boldsymbol{r}(\boldsymbol{x})=[\boldsymbol{r}_{i}(\boldsymbol{x}_{i},\boldsymbol{x}_{\mathcal{N}_{i}})]_{i\in\mathcal{V}}$
can be synthesized in the considered nonlinear continuous-time setting such
that stability is ensured. In cooperative distributed MPC, the central OCP (3)
is decomposed into smaller subproblems which are solved in parallel by the
agents via a distributed optimization algorithm. After an optimal input has
been found, each agent locally applies the controls

$\boldsymbol{\underaccent{\bar}{u}}_{i}(t_{k}+\tau)=\boldsymbol{u}_{i}^{*}(\tau;\boldsymbol{\underaccent{\bar}{x}}_{k})\,,\quad\tau\in[0,\,\Delta
t)\,,\quad i\in\mathcal{V}$ (14)

to the subsystems (1). Compared to the central case, the terminal feedback law
$\boldsymbol{r}(\boldsymbol{x})$ is not allowed to consider the complete
system state in the distributed setting, i.e.
$\boldsymbol{r}_{i}(\boldsymbol{x})$ for all $i\in\mathcal{V}$, since this
violates the requirement of only neighbor-to-neighbor communication. Rather
$\boldsymbol{r}_{i}$ can only involve neighboring states
$\boldsymbol{x}_{j},\,j\in\mathcal{N}_{i}$ which can be communicated, i.e.
$\boldsymbol{r}_{i}(\boldsymbol{x}_{i},\boldsymbol{x}_{\mathcal{N}_{i}})$.
However, classical central MPC stability proofs, e.g. the works 68, 75, 73,
74, rely on a full state feedback law. In addition, the required separability
in (4) imposes an additional structural constraint on the CLFs. These two
additional requirements need to be accounted for when transferring Theorem 2.1
to the DMPC setting. Therefore, it is necessary to replace Assumption 2.1 with
the following strengthened assumption. There exists a structured feedback law
$\boldsymbol{u}=[\boldsymbol{r}_{i}(\boldsymbol{x}_{i},\boldsymbol{x}_{\mathcal{N}_{i}})]_{i\in\mathcal{V}}\in\mathbb{U}$
and a structured CLF, i.e.
$V(\boldsymbol{x})=\sum_{i\in\mathcal{V}}V_{i}(\boldsymbol{x}_{i})$ such that
the CLF inequality (11) holds for all
$[\boldsymbol{x}_{i}]_{i\in\mathcal{V}}\in\Omega_{\beta}=\\{\boldsymbol{x}\in\Gamma:V(\boldsymbol{x})\leq\beta\\}$.
The assumption concerning the existence of such a structured feedback
controller is typical in the DMPC setting 20, 66, 67, 19 and has recently been
shown to be satisfied by a broad class of systems 79. Furthermore, this
assumption is not as restrictive as requiring local stabilizing controllers,
i.e.
$\boldsymbol{u}=[\boldsymbol{r}_{i}(\boldsymbol{x}_{i})]_{i\in\mathcal{V}}$
28. The derivation of general conditions for the existence of such (linear)
structured feedback controllers is beyond the scope of this paper and poses a
difficult problem in itself. Related literature includes the reference80,
where conditions on the decentralized stabilizability of the linearized
dynamics (2) are discussed, or the work81 where conditions on the coupling
strength between agents are derived such that a structured controller exists.
Furthermore, algorithms for designing structured feedback controllers for
continuous-time systems can be found e.g. in the contributions 82, 83. The
following corollary reveals that Theorem 2.1 holds in the same manner for the
optimal distributed setting under the strengthened Assumption 2.2.

######  Corollary 2.2.

Suppose that Assumption 2.2 holds. Then, the origin of the closed-loop system
under the DMPC control law (14) is asymptotically stable for all
$\boldsymbol{\underaccent{\bar}{x}}_{0}\in\Gamma_{\alpha}$ and the optimal
cost (9) decays asymptotically.

The proof of Corollary 2.2 is analogous to the proof of Theorem 2.1 as
Assumption 2.2 represents a special case of Assumption 2.1. It shows that the
classical MPC result of Theorem 2.1 extends to the distributed setting if it
is additionally assumed that the feedback control
$\boldsymbol{r}(\boldsymbol{x})$ and terminal cost $V(\boldsymbol{x})$ are
distributed as well. Note that this result is independent of the underlying
distributed optimization algorithm. The case in which only a suboptimal
solution with a limited amount of iterations is found by the considered
sensitivity-based distributed optimal control algorithm will be investigated
in Section 4.

In the following, a simple optimization-based approach is proposed to obtain a
structured terminal cost $V(\boldsymbol{x})$ and controller
$\boldsymbol{u}=\boldsymbol{r}(\boldsymbol{x})$ offline such that Assumption
2.2 is satisfied. For the remainder of this section, we will consider the
specific quadratic terminal and integral cost functions

$V(\boldsymbol{x})=\sum_{i\in\mathcal{V}}\boldsymbol{x}_{i}^{\mathrm{T}}\boldsymbol{P}_{i}\boldsymbol{x}_{i}\,,\quad
l(\boldsymbol{x},\boldsymbol{u})=\sum_{i\in\mathcal{V}}\boldsymbol{x}_{i}^{\mathrm{T}}\boldsymbol{Q}_{i}\boldsymbol{x}_{i}+\boldsymbol{u}_{i}^{\mathrm{T}}\boldsymbol{R}_{i}\boldsymbol{u}_{i}$
(15)

with weighting matrices $\boldsymbol{P}_{i},\,\boldsymbol{Q}_{i}\succ 0$ and
$\boldsymbol{R}_{i}\succ 0$, $i\in\mathcal{V}$. The main difficulty and
difference to classical MPC are that the quadratic terminal cost
$V(\boldsymbol{x})=\boldsymbol{x}^{\mathrm{T}}\boldsymbol{P}\boldsymbol{x}$
and linear feedback law
$\boldsymbol{r}(\boldsymbol{x})=\boldsymbol{K}\boldsymbol{x}$ cannot be
obtained by solving a Lyapunov or Riccati equation.20, 66 This is due to the
fact that in general the solution matrix $\boldsymbol{P}$ and feedback
controller $\boldsymbol{K}$ will not exhibit the separable structure as
required by Assumption 2.2. Furthermore, consider the linearized system of (2)
at the origin which can be derived as

$\boldsymbol{\dot{x}}=\boldsymbol{A}\boldsymbol{x}+\boldsymbol{B}\boldsymbol{u}\,,\quad\boldsymbol{x}(0)=\boldsymbol{\underaccent{\bar}{x}}_{0}\,,$
(16)

where
$\boldsymbol{A}=\partial_{\boldsymbol{x}}\,\boldsymbol{f}(\boldsymbol{0},\boldsymbol{0})$
and
$\boldsymbol{B}:=\partial_{\boldsymbol{u}}\,\boldsymbol{f}(\boldsymbol{0},\boldsymbol{0})$.
Note that similar to (1), the linearized system exhibits a sparse structure
with
$\partial_{\boldsymbol{x}_{j}}\,\boldsymbol{f}_{i}(\boldsymbol{0},\boldsymbol{0},\boldsymbol{0})=:\boldsymbol{A}_{ij}\neq\boldsymbol{0}$
only if $j\in\mathcal{N}_{i}^{\leftarrow}$ and
$\boldsymbol{B}=\mathrm{diag}(\boldsymbol{B}_{1},\boldsymbol{B}_{2},...,\boldsymbol{B}_{N})$
with
$\boldsymbol{B}_{i}=\partial_{\boldsymbol{u}_{i}}\,\boldsymbol{f}_{i}(\boldsymbol{0},\boldsymbol{0},\boldsymbol{0})$.
By Assumption 2.1 there exists a stabilizing controller
$\boldsymbol{u}=\boldsymbol{K}\boldsymbol{x}$ such that the controlled linear
system $\boldsymbol{A}+\boldsymbol{B}\boldsymbol{K}$ is stable. It is well
known from MPC literature that if the inequality

$(\boldsymbol{A}+\boldsymbol{B}\boldsymbol{K})^{\mathrm{T}}\boldsymbol{P}+\boldsymbol{P}(\boldsymbol{A}+\boldsymbol{B}\boldsymbol{K})+\gamma(\boldsymbol{Q}+\boldsymbol{K}^{\mathrm{T}}\boldsymbol{R}\boldsymbol{K})\leq
0$ (17)

with $\gamma\in(1,\infty)$ is satisfied, then there exists a non-empty set
$\Omega_{\beta}=\\{\boldsymbol{x}\in\Gamma:\boldsymbol{x}^{\mathrm{T}}\boldsymbol{P}\boldsymbol{x}\leq\beta\\}$
for some $\beta>0$ such that $\Omega_{\beta}$ is control invariant with the
control law
$\boldsymbol{\underaccent{\bar}{u}}=\boldsymbol{K}\boldsymbol{\underaccent{\bar}{x}}$,
i.e. any initial state
$\boldsymbol{\underaccent{\bar}{x}}_{0}\in\Omega_{\beta}$ implies that
$\boldsymbol{\underaccent{\bar}{x}}(t)\in\Omega_{\beta}$ and
$\boldsymbol{\underaccent{\bar}{u}}(t)\in\mathbb{U}$ for all $t$ and that for
any $\boldsymbol{\underaccent{\bar}{x}}(t)\in\Omega_{\beta}$, the CLF
inequality
$\frac{\mathrm{d}}{\mathrm{d}t}(\boldsymbol{x}^{\mathrm{T}}\boldsymbol{P}\boldsymbol{x})+\boldsymbol{x}^{\mathrm{T}}(\boldsymbol{\boldsymbol{Q}+\boldsymbol{K}^{\mathrm{T}}\boldsymbol{R}\boldsymbol{K}})\boldsymbol{x}\leq
0$ holds 75, 2. The requirements of Assumption 2.2 and condition (17) can be
reformulated as the following convex semi-definite optimization problem in
$\boldsymbol{P}$ and $\boldsymbol{K}$

$\displaystyle\min_{\boldsymbol{E}\succeq\boldsymbol{0},\,\boldsymbol{Y}}\quad$
$\displaystyle-\log(\det(\boldsymbol{E}))$ (18a)
$\displaystyle\operatorname{s.t.}\quad$
$\displaystyle\begin{bmatrix}\boldsymbol{A}\boldsymbol{E}+\boldsymbol{E}\boldsymbol{A}^{\mathrm{T}}+\boldsymbol{B}\boldsymbol{Y}+\boldsymbol{Y}^{\mathrm{T}}\boldsymbol{B}^{\mathrm{T}}&\boldsymbol{E}\boldsymbol{Q}^{\frac{1}{2}}&\boldsymbol{Y}^{\mathrm{T}}\boldsymbol{R}^{\frac{1}{2}}\\\
\boldsymbol{Q}^{\frac{1}{2}}\boldsymbol{E}&-\frac{1}{\gamma}\boldsymbol{I}&\boldsymbol{0}\\\
\boldsymbol{R}^{\frac{1}{2}}\boldsymbol{Y}&\boldsymbol{0}&-\frac{1}{\gamma}\boldsymbol{I}\end{bmatrix}\preceq\boldsymbol{0}$
(18b)
$\displaystyle\boldsymbol{E}=\mathrm{diag}\\{\boldsymbol{E}_{i}\,:\,i\in\mathcal{V}\\}$
(18c) $\displaystyle\boldsymbol{Y}_{ij}=\boldsymbol{0}\,,\quad\quad
i\in\mathcal{V},\,j\notin\mathcal{N}_{i}\,,$ (18d)

with the substitutions $\boldsymbol{E}:=\boldsymbol{P}^{-1}$ and
$\boldsymbol{Y}:=\boldsymbol{K}\boldsymbol{E}$. The equivalence of conditions
(17) and (18b) can be shown via Schur compliment techniques84. For the
considered quadratic cost functions, the constraint (18c) ensures that the
terminal cost $V(\boldsymbol{x})$ has the required separable structure, while
the constraint (18d) guarantees that the controller fulfills Assumption 2.2,
i.e.
$\boldsymbol{r}_{i}(\boldsymbol{x}_{i},\boldsymbol{x}_{\mathcal{N}_{i}})=\boldsymbol{K}_{ii}\boldsymbol{x}_{i}+\sum_{j\in\mathcal{N}_{i}}\boldsymbol{K}_{ij}\boldsymbol{x}_{j}$.
The cost objective requires maximizing $\log(\det\boldsymbol{P}^{-1})$ which
is favorable in regard of enlarging the size of the terminal set as it
maximizes the volume of the 1-level set ellipsoid
$\Omega_{1}=\\{\boldsymbol{x}\in\Gamma:\boldsymbol{x}^{\mathrm{T}}\boldsymbol{P}\boldsymbol{x}\leq
1\\}$ 20. Since problem (18) is convex, global minimizerers $\boldsymbol{P}$
and $\boldsymbol{K}$ with the required structure exists as long as (18b) as
well as $\boldsymbol{E}\succeq\boldsymbol{0}$ are feasible. In this way, the
optimization procedure (18) provides a constructive test for the existence of
a separable terminal cost acting as a CLF and a structured terminal controller
as both required by Assumption 2.2. Conservatism, however, is introduced by
the choice of a block diagonal terminal cost and the structured linear
controller when compared to the classical MPC approach without structural
constraints. The size of the terminal region, i.e. $\beta$, can for example be
determined similar to the procedure in 75.

##  3 Distributed solution via sensitivities

The main idea of iterative and cooperative DMPC schemes is to solve (3) in a
distributed fashion in each MPC step $k$. To this end, a sensitivity-based
approach is used in which the local cost function (5) of each agent is
extended by so-called sensitivities of its neighbors such that the agents
cooperatively solve the central OCP (3).

###  3.1 Extended local optimal control problem

The idea of utilizing first-order sensitivities to account for the influence
of neighbors has been utilized before in the DMPC setting 35, 59, 8, 57, 37
and is derived for the problem at hand. By extending the individual cost
functions (4) of each agent $i\in\mathcal{V}$ with the first-order
sensitivities of its neighbors $j\in\mathcal{N}_{i}$, the agents not only
consider their local cost objective but also a first-order approximation of
their neighbors’ costs. The sensitivity represents the information about the
expected change in the overall cost objective (3a) that the states of an agent
$i\in\mathcal{V}$ induce via the local cost of agent
$j\in\mathcal{N}_{i}^{\rightarrow}$. The local problem of agent
$i\in\mathcal{V}$ to be solved at step $q$ of the sensitivity-based algorithm
therefore reads

$\displaystyle\min_{\boldsymbol{u}_{i}}\quad$
$\displaystyle\bar{J}_{i}(\boldsymbol{u}_{i};\boldsymbol{\underaccent{\bar}{x}}_{k}):=V_{i}(\boldsymbol{x}_{i}(T))+\int_{0}^{T}l_{i}(\boldsymbol{x}_{i},\boldsymbol{u}_{i},\boldsymbol{x}_{\mathcal{N}_{i}^{\leftarrow}}^{q-1})\,\mathrm{d}\tau+\sum_{j\in\mathcal{N}_{i}^{\rightarrow}}\delta\bar{J}_{j}(\boldsymbol{u}_{j}^{q-1};\boldsymbol{\underaccent{\bar}{x}}_{k})(\delta\boldsymbol{x}_{i})$
(19a) $\displaystyle\leavevmode\nobreak\ \operatorname{s.t.}\quad$
$\displaystyle\boldsymbol{\dot{x}}_{i}=\boldsymbol{f}_{i}(\boldsymbol{x}_{i},\boldsymbol{u}_{i},\boldsymbol{x}_{\mathcal{N}_{i}^{\leftarrow}}^{q-1})\,,\quad\boldsymbol{x}(0)=\boldsymbol{\underaccent{\bar}{x}}_{i,k}$
(19b)
$\displaystyle\boldsymbol{u}_{i}(\tau)\in\mathbb{U}_{i}\,,\quad\tau\in[0,\,T]$
(19c)

with the neighbor’s cost sensitivity represented by the corresponding Gâteaux
derivative
$\delta\bar{J}_{j}(\boldsymbol{u}_{j}^{q-1};\boldsymbol{\underaccent{\bar}{x}}_{k})(\delta\boldsymbol{x}_{i})$,
$j\in\mathcal{N}_{i}^{\rightarrow}$ in direction of
$\delta\boldsymbol{x}_{i}:=\boldsymbol{x}_{i}-\boldsymbol{x}_{i}^{q-1}$
evaluated for the trajectories of the previous iteration $q-1$. The dynamics
(1) and costs (5) are decoupled by using the state trajectories
$\boldsymbol{x}_{\mathcal{N}_{i}^{\leftarrow}}^{q-1}$ of the previous
iteration which need to be communicated by each neighbor
$j\in\mathcal{N}_{i}^{\leftarrow}$ to the agent $i\in\mathcal{V}$. The Gâteaux
derivative $\delta\bar{J}_{j}$ in (19a) can be expressed as (see Appendix A
for an explicit derivation)

$\displaystyle\delta\bar{J}_{j}(\boldsymbol{u}_{j}^{q-1};\boldsymbol{\underaccent{\bar}{x}}_{k})(\delta\boldsymbol{x}_{i})=\int_{0}^{T}\big{(}\underbrace{\partial_{\boldsymbol{x}_{i}}l_{ji}(\boldsymbol{x}_{j}^{q-1},\boldsymbol{x}_{i}^{q-1})+\partial_{\boldsymbol{x}_{i}}\,\boldsymbol{f}_{ji}(\boldsymbol{x}_{j}^{q-1},\boldsymbol{x}_{i}^{q-1})^{\mathrm{T}}\boldsymbol{\lambda}_{j}^{q-1}}_{\boldsymbol{g}_{ji}^{q-1}(\tau)}\big{)}^{\mathrm{T}}\,\delta\boldsymbol{x}_{i}\,\mathrm{d}\tau\,,\quad
j\in\mathcal{N}_{i}^{\rightarrow}$ (20)

where $\boldsymbol{\lambda}_{j}\in\mathds{R}^{n_{j}}$ denotes the adjoint
state of the neighbors $j\in\mathcal{N}_{i}^{\rightarrow}$. Defining the local
Hamiltonian for the local OCP (19)

$H_{i}(\boldsymbol{x}_{i},\boldsymbol{u}_{i},\boldsymbol{\lambda}_{i}):=l_{i}(\boldsymbol{x}_{i},\boldsymbol{u}_{i},\boldsymbol{x}_{\mathcal{N}_{i}^{\leftarrow}}^{q-1})+\boldsymbol{\lambda}_{i}^{\mathrm{T}}\boldsymbol{f}_{i}(\boldsymbol{x}_{i},\boldsymbol{u}_{i},\boldsymbol{x}_{\mathcal{N}_{i}^{\leftarrow}}^{q-1})+\sum_{j\in\mathcal{N}_{i}^{\rightarrow}}(\,\boldsymbol{g}_{ji}^{q-1})^{\mathrm{T}}(\boldsymbol{x}_{i}-\boldsymbol{x}_{i}^{q-1})\,,$
(21)

the adjoint state can be calculated via the backward integration of

$\boldsymbol{\dot{\lambda}}_{i}(\tau)=-\partial_{\boldsymbol{x}_{i}}H_{i}(\boldsymbol{x}_{i},\boldsymbol{u}_{i},\boldsymbol{\lambda}_{i})\,,\quad\boldsymbol{\lambda}_{i}(T)=\partial_{\boldsymbol{x}_{i}}V_{i}(\boldsymbol{x}_{i}(T;\boldsymbol{\underaccent{\bar}{x}}_{k}))\,,\quad\tau\in[0,\,T]\,.$
(22)

and needs to be communicated by each agent $i\in\mathcal{V}$ to its neighbors
$j\in\mathcal{N}_{i}^{\leftarrow}$ such that (20) can be evaluated. Note that
for sending neighbors $j\in\mathcal{N}_{i}^{\leftarrow}$ the sensitivity (20)
vanishes since $l_{ji}=0$ and $\boldsymbol{f}_{ji}=\boldsymbol{0}$. The
bracketed term in (20) can be interpreted as the (time-dependent) gradient
$\boldsymbol{g}_{ji}(\tau)\in\mathds{R}^{n_{i}},\tau\in[0,\,T]$ of a
neighbor’s $j\in\mathcal{N}_{i}^{\rightarrow}$ cost functional w.r.t. the
(time-dependent) external trajectories $\boldsymbol{x}_{i}(\tau)$ of agent
$i\in\mathcal{V}$. The gradient $\boldsymbol{g}_{ji}^{q-1}$ can be computed
locally by each agent for its neighbors $j\in\mathcal{N}_{i}^{\rightarrow}$
due to the neighbor-affine structure of (1) and (5) as long as the agent has
access to the coupling functions
$l_{ji}(\boldsymbol{x}_{j},\boldsymbol{x}_{i})$ and
$\boldsymbol{f}_{ji}(\boldsymbol{x}_{j},\boldsymbol{x}_{i})$,
$j\in\mathcal{N}_{i}^{\rightarrow}$.

###  3.2 Sensitivity-based algorithm

The separable structure of OCP (19) is exploited by solving the individual
problems in parallel on the agent level. The procedure is summarized in
Algorithm 1 which shows the distributed sensitivity-based solution of the
original central OCP (3). The algorithm is executed in each MPC step with the
current system state
$\boldsymbol{\underaccent{\bar}{x}}_{k}=[\boldsymbol{\underaccent{\bar}{x}}_{i,k}]_{i\in\mathcal{V}}=\boldsymbol{\underaccent{\bar}{x}}(t_{k})$
in parallel and locally by each agent $i\in\mathcal{V}$. Algorithm 1 consists
of two local computations of which the first calculates the gradient
$\boldsymbol{g}_{ji}$ for all neighbors $j\in\mathcal{N}_{i}^{\rightarrow}$
and the second requires the solution of the local OCP (19). In addition, one
communication step is needed in which the state and adjoint trajectories are
exchanged between neighboring agents. Thus, Algorithm 1 constitutes a fully
distributed algorithm with only local computation steps and neighbor-to-
neighbor communication. For an efficient solution of OCP (27) the fixed-point
iteration scheme63 or the projected gradient method73, 85 can be used which
also conveniently calculate the adjoint state $\boldsymbol{\lambda}_{i}$.
Alternatively, a standard OCP solver can be used to compute
$(\boldsymbol{u}_{i}^{k},\boldsymbol{x}_{i}^{k})$ in (27), but then the
adjoint dynamics (22) must be solved backward in time to obtain
$\boldsymbol{\lambda}_{i}^{k}$.

The stopping criterion in Step 4 terminates the algorithm either if a maximum
number of iterations $q_{\mathrm{max}}$ has been reached or if a convergence
criterion is satisfied at a certain iteration number $q_{k}$ in the current
MPC step. The criterion (28) evaluates the progress of the state and adjoint
state trajectories
($\boldsymbol{x}_{i}^{q}(\tau),\boldsymbol{\lambda}_{i}^{q}(\tau)$) between
two algorithm iterations w.r.t. the norm of the current system state
$\boldsymbol{\underaccent{\bar}{x}}_{k}$ weighted by the constant $d>0$. This
leads to a contraction of the stopping criterion during the stabilization of
the system to the origin. Although the stopping criterion is only evaluated
locally, it needs to be satisfied by all agents simultaneously in each MPC
step which requires either a central node or a global communication step. In a
practical implementation, this convergence monitoring is usually done by a
central node which is also responsible for triggering the individual algorithm
steps in a synchronized manner.25

This central convergence check can be avoided if instead a maximum number of
iterations $q_{\mathrm{max}}$ is performed without evaluating the stopping
criterion. In this way, the convergence criterion (28) can be used if a
certain accuracy is required, while stopping after a maximum number of
iterations is beneficial for real-time applications with a fixed-sampling time
$\Delta t$. In the next section, conditions on the stopping criterion constant
$d$ and the maximum number of iterations $q_{max}$ are derived such that the
resulting DMPC controller stabilizes the system exponentially.

In the first MPC step $k=0$, Algorithm 1 is locally initialized with
appropriate trajectories, e.g.
$\boldsymbol{x}_{i}^{0}(\tau;\boldsymbol{x}_{0})=\boldsymbol{\underaccent{\bar}{x}}_{i,0}$
and
$\boldsymbol{\lambda}_{i}^{0}(\tau,\boldsymbol{x}_{0})=\partial_{\boldsymbol{x}_{i}}V_{i}(\boldsymbol{\underaccent{\bar}{x}}_{i,0})$.
In the following MPC steps $k\geq 1$, the algorithm is warm-started by the
solution trajectories of the last iteration $q_{k}$ in the previous MPC step,
i.e.

$\displaystyle\boldsymbol{x}_{i}^{0}(\tau;\boldsymbol{\underaccent{\bar}{x}}_{k+1})=\boldsymbol{x}_{i}^{q_{k}}(\tau;\boldsymbol{\underaccent{\bar}{x}}_{k})\,,\quad\boldsymbol{\lambda}_{i}^{0}(\tau;\boldsymbol{\underaccent{\bar}{x}}_{k+1})=\boldsymbol{\lambda}_{i}^{q_{k}}(\tau;\boldsymbol{\underaccent{\bar}{x}}_{k})\,,\quad\tau\in[0,\,T]\,,\quad
i\in\mathcal{V}$ (23)

for all agents $i\in\mathcal{V}$. As soon as the stopping criterion (28) is
fulfilled for all agents $i\in\mathcal{V}$ or the maximum number of iterations
$q_{\mathrm{max}}$ are reached, the calculated control trajectory
$\boldsymbol{u}_{i}^{q_{k}}(\tau;\boldsymbol{\underaccent{\bar}{x}}_{k})$ from
the last iteration $q_{k}$ is taken as the input for the actual subsystems
(1), i.e.

$\boldsymbol{\underaccent{\bar}{u}}_{i}(t_{k}+\tau)=\boldsymbol{u}_{i}^{q_{k}}(\tau;\boldsymbol{\underaccent{\bar}{x}}_{k})\,,\quad\tau\in[0,\,\Delta
t)\,,\quad i\in\mathcal{V}.$ (24)

Similar to (10), the current controls
$\boldsymbol{u}^{q_{k}}(\tau;\boldsymbol{\underaccent{\bar}{x}}_{k})=[\boldsymbol{u}_{i}]_{i\in\mathcal{V}}$
can interpreted as the nonlinear sampled control law

$\boldsymbol{u}^{q_{k}}(\tau;\boldsymbol{\underaccent{\bar}{x}}_{k})=:\boldsymbol{\kappa}(\boldsymbol{x}^{q_{k}}(\tau;\boldsymbol{\underaccent{\bar}{x}}_{k});\boldsymbol{x}^{q_{k}-1}(\tau;\boldsymbol{\underaccent{\bar}{x}}_{k}),\boldsymbol{\lambda}^{q_{k}-1}(\tau;\boldsymbol{\underaccent{\bar}{x}}_{k}),\boldsymbol{\underaccent{\bar}{x}}_{k})$
(25)

with $\tau\in[0,\,\Delta t)$. Note, however, in contrast to the centralized
MPC case (10), the control law (25) is parameterized by the state and adjoint
trajectories
$\boldsymbol{x}^{q_{k}-1}(\tau;\boldsymbol{\underaccent{\bar}{x}}_{k})$ and
$\boldsymbol{\lambda}^{q_{k}-1}(\tau;\boldsymbol{\underaccent{\bar}{x}}_{k})$
of the previous iteration of Algorithm 1.

Algorithm 1 Distributed sensitivity-based solution of OCP (3)

1:Initialize $(\boldsymbol{x}_{i}^{0}$, $\boldsymbol{\lambda}_{i}^{0})$, set
$q=1$, choose $d>0$ or $q_{\mathrm{max}}$, send $\boldsymbol{x}_{i}^{0}$ to
$j\in\mathcal{N}_{i}$ and $\boldsymbol{\lambda}_{i}^{0}$ to
$j\in\mathcal{N}_{i}^{\leftarrow}$

2:Compute the gradients $\boldsymbol{g}_{ji}^{q-1}$ for all
$j\in\mathcal{N}_{i}^{\rightarrow}$ as
$\displaystyle\boldsymbol{g}_{ji}^{q-1}(\tau)=\partial_{\boldsymbol{x}_{i}}l_{ji}(\boldsymbol{x}_{j}^{q-1},\boldsymbol{x}_{i}^{q-1})+\partial_{\boldsymbol{x}_{i}}\,\boldsymbol{f}_{ji}(\boldsymbol{x}_{j}^{q-1},\boldsymbol{x}_{i}^{q-1})^{\mathrm{T}}\boldsymbol{\lambda}_{j}^{q-1}\,,\quad\tau\in[0,\,T]$
(26)

3:Compute
$(\boldsymbol{u}_{i}^{q},\boldsymbol{x}_{i}^{q},\boldsymbol{\lambda}_{i}^{q})$
by solving  $\displaystyle\min_{\boldsymbol{u}_{i}}\quad$ $\displaystyle
V_{i}(\boldsymbol{x}_{i}(T))+\int_{0}^{T}l_{i}(\boldsymbol{x}_{i},\boldsymbol{u}_{i},\boldsymbol{x}_{\mathcal{N}_{i}^{\leftarrow}}^{q-1})+\\!\\!\sum_{j\in\mathcal{N}_{i}^{\rightarrow}}(\boldsymbol{g}_{ji}^{q-1})^{\mathrm{T}}(\boldsymbol{x}_{i}-\boldsymbol{x}_{i}^{q-1})\,\mathrm{d}\tau$
(27a) $\displaystyle\leavevmode\nobreak\ \operatorname{s.t.}\quad$
$\displaystyle\boldsymbol{\dot{x}}_{i}=\boldsymbol{f}_{i}(\boldsymbol{x}_{i},\boldsymbol{u}_{i},\boldsymbol{x}_{\mathcal{N}_{i}^{\leftarrow}}^{q-1})\,,\quad\boldsymbol{x}_{i}(0)=\boldsymbol{\underaccent{\bar}{x}}_{i,k}$
(27b)
$\displaystyle\boldsymbol{u}_{i}(\tau)\in\mathbb{U}_{i}\,,\quad\tau\in[0,\,T]$
(27c)

4:Send state trajectories $\boldsymbol{x}_{i}^{q}$ to $j\in\mathcal{N}_{i}$
and adjoint trajectories $\boldsymbol{\lambda}_{i}^{q}$ to
$j\in\mathcal{N}_{i}^{\leftarrow}$

5:if $q=q_{\mathrm{max}}$ or
$\bigg{\|}\begin{matrix}\boldsymbol{x}_{i}^{q}-\boldsymbol{x}_{i}^{q-1}\\\
\boldsymbol{\lambda}_{i}^{q}-\boldsymbol{\lambda}_{i}^{q-1}\\\
\end{matrix}\bigg{\|}_{L_{\infty}}\leq d\|\underaccent{\bar}{\vm
x}_{i,k}\|\,,\quad\forall i\in\mathcal{V}$ (28)

6:then Quit with $q_{k}:=q$

7:else Increment $q\leftarrow q+1$ and go to Step 1

##  4 DMPC stability analysis under inexact optimization

By construction, Algorithm 1 iterates until a suboptimal input
$\boldsymbol{u}_{i}^{q_{k}}(\tau;\boldsymbol{\underaccent{\bar}{x}}_{k})$,
$\tau\in[0,\,\Delta t)$ is found, either when the stopping criterion (28) is
fulfilled or a maximum number of iterations $q_{\mathrm{max}}$ is reached.
While Algorithm 1 is consequently suitable for a real-time capable
implementation, stability of the DMPC scheme is not guaranteed since Corollary
2.2 does not apply anymore. That is why stability results for the DMPC scheme
under an inexact minimization of the central OCP (3) by Algorithm 1 are
analyzed in this section. Two scenarios are considered: First an upper bound
on the constant $d$ in (28) is derived in Section 4.2 and then the bound is
reformulated in Section 4.3 in terms of the required number of iterations
$q_{\mathrm{max}}$ that must be executed in each MPC step. The proof follows
along the lines of the work24, where the (D)MPC stability of a similar scheme,
based on ADMM, was investigated. However, we require less assumptions and in
particular do not need to assume convergence and boundedness of the iterates
of the employed distributed optimization algorithm but rather are able to show
linear convergence of Algorithm 1 which turns out to be critical for the
subsequent stability analysis of the DMPC scheme. In contrast the trade-off is
that convergence of Algorithm 1 can only be guaranteed for some maximum
horizon length which is investigated in more detail in Lemma 4.1 at the end of
Section 4.1.

###  4.1 Preliminaries for the suboptimal case

Compared to the optimal (D)MPC feedback laws (10) and (14), the control (25)
represents a suboptimal feedback law. This is caused by prematurely stopping
Algorithm 1 after a certain number of iterations. As a consequence, the
solution trajectories of the sensitivity-based DMPC algorithm are not
consistent with the optimal MPC solution and could lead to a possibly
destabilizing solution as Theorem 2.1 or Corollary 2.2 do not hold anymore.

While stopping Algorithm 1 prematurely limits the number of required
iterations and thus the time-expensive communication steps, it also leads to
suboptimal control trajectories
$\boldsymbol{u}^{q_{k}}(\tau;\boldsymbol{\underaccent{\bar}{x}}_{k})$ that
differ from the optimal ones
$\boldsymbol{u}^{*}(\tau;\boldsymbol{\underaccent{\bar}{x}}_{k})$. That is why
in general, the state trajectories of the local MPC predictions as part of the
solution of the local OCP (27), the actual realizations from applying the
suboptimal control (25) of each agent to the actual system (2), and the
optimal ones resulting from the centralized MPC solution (3) will not be
identical. In particular, we have to distinguish between the following
trajectories:

  * •

The individual predicted state trajectories
$\boldsymbol{x}^{q_{k}}(\cdot;\boldsymbol{\underaccent{\bar}{x}}_{k})=[\boldsymbol{x}_{i}^{q_{k}}(\cdot;\boldsymbol{\underaccent{\bar}{x}}_{k})]_{i\in\mathcal{V}}$
resulting from the solution of OCP (27) in the last iteration $q_{k}$, i.e.

$\boldsymbol{\dot{x}}^{q_{k}}(\tau;\boldsymbol{\underaccent{\bar}{x}}_{k})=\boldsymbol{\hat{f}}(\boldsymbol{x}^{q_{k}}(\tau;\boldsymbol{\underaccent{\bar}{x}}_{k}),\boldsymbol{u}^{q_{k}}(\tau;\boldsymbol{\underaccent{\bar}{x}}_{k}),\boldsymbol{\hat{x}}^{q_{k}-1}(\tau;\boldsymbol{\underaccent{\bar}{x}}_{k}))\,,\quad\boldsymbol{x}^{q_{k}}(0;\boldsymbol{\underaccent{\bar}{x}}_{k})=\boldsymbol{\underaccent{\bar}{x}}_{k}$
(29)

with the stacked notation
$\boldsymbol{\hat{x}}^{q_{k}-1}(\tau;\boldsymbol{\underaccent{\bar}{x}}_{k})=[\boldsymbol{x}_{\mathcal{N}_{i}^{\leftarrow}}^{q_{k}-1}]_{i\in\mathcal{V}}\in\mathds{R}^{p}$
and $\boldsymbol{\hat{f}}=[\boldsymbol{f}_{i}]_{i\in\mathcal{V}}$, where the
notation $\boldsymbol{\hat{f}}$ explicitly captures the dependency on
$\boldsymbol{\hat{x}}^{q_{k}-1}(\tau;\boldsymbol{\underaccent{\bar}{x}}_{k})$.

  * •

The actual state trajectories
$\boldsymbol{x}_{c}^{q_{k}}(\cdot;\boldsymbol{\underaccent{\bar}{x}}_{k})$
resulting from applying the suboptimal controls
$\boldsymbol{u}^{q_{k}}(\cdot;\boldsymbol{\underaccent{\bar}{x}}_{k})=[\boldsymbol{u}_{i}^{q_{k}}(\cdot;\boldsymbol{\underaccent{\bar}{x}}_{k})]_{i\in\mathcal{V}}$
to the central system (2), i.e.

$\boldsymbol{\dot{x}}_{c}^{q_{k}}(\tau;\boldsymbol{\underaccent{\bar}{x}}_{k})=\boldsymbol{f}(\boldsymbol{x}_{c}^{q_{k}}(\tau;\boldsymbol{\underaccent{\bar}{x}}_{k}),\boldsymbol{u}^{q_{k}}(\tau;\boldsymbol{\underaccent{\bar}{x}}_{k}))\,,\quad\boldsymbol{x}_{c}^{q_{k}}(0;\boldsymbol{\underaccent{\bar}{x}}_{k})=\boldsymbol{\underaccent{\bar}{x}}_{k}\,.$
(30)

  * •

The optimal trajectories
$\boldsymbol{x}^{*}(\cdot;\boldsymbol{\underaccent{\bar}{x}}_{k})=[\boldsymbol{x}_{i}^{*}(\cdot;\boldsymbol{\underaccent{\bar}{x}}_{k})]_{i\in\mathcal{V}}$
following from solving the central OCP (3) cf. (8a), i.e.

$\boldsymbol{\dot{x}}^{*}(\tau;\boldsymbol{\underaccent{\bar}{x}}_{k})=\boldsymbol{f}(\boldsymbol{x}^{*}(\tau;\boldsymbol{\underaccent{\bar}{x}}_{k}),\boldsymbol{u}^{*}(\tau;\boldsymbol{\underaccent{\bar}{x}}_{k}))\,,\quad\boldsymbol{x}^{*}(0;\boldsymbol{\underaccent{\bar}{x}}_{k})=\boldsymbol{\underaccent{\bar}{x}}_{k}\,.$
(31)

In the nominal case, this implies that the system state of (2) in the next MPC
step, i.e.
$\boldsymbol{\underaccent{\bar}{x}}_{k+1}=\boldsymbol{\underaccent{\bar}{x}}(t_{k+1})$,
will lie on the actual state trajectory

$\boldsymbol{\underaccent{\bar}{x}}_{k+1}=\boldsymbol{x}_{c}^{q_{k}}(\Delta
t;\boldsymbol{\underaccent{\bar}{x}}_{k})$ (32)

and not on the individual state trajectory $\boldsymbol{x}^{q_{k}}(\Delta
t;\boldsymbol{\underaccent{\bar}{x}}_{k})$ or the optimal one
$\boldsymbol{x}^{*}(\Delta t;\boldsymbol{\underaccent{\bar}{x}}_{k})$. This
difference between optimal and actual state trajectory is expressed by the
error between actual state trajectory and optimal state trajectory

$\Delta\boldsymbol{x}_{c}^{q_{k}}(\tau;\boldsymbol{\underaccent{\bar}{x}}_{k}):=\boldsymbol{x}_{c}^{q_{k}}(\tau;\boldsymbol{\underaccent{\bar}{x}}_{k})-\boldsymbol{x}^{*}(\tau;\boldsymbol{\underaccent{\bar}{x}}_{k})\,,\quad\tau\in[0,\,T]$
(33)

that can be interpreted as a suboptimality measure in each MPC step $k$.
Several assumptions and intermediate results concerning the optimal solution
of the central MPC scheme and the properties of Algorithm 1 are necessary to
proceed. The optimal cost (9) is twice continuously differentiable and the
feedback laws (10) and (25) are locally Lipschitz w.r.t. their arguments.
These continuity assumptions are typical in (D)MPC 72, 68 and are required to
derive certain bounds in the stability analysis. A consequence of Assumption
4.1 is that there exist constants $m_{J}$, $M_{J}>0$ such that the optimal
cost satisfies (cf. (Sensitivity-Based Distributed Model Predictive Control
for Nonlinear Systems under Inexact Optimization) and (104) in Appendix E)

$m_{J}\|\boldsymbol{\underaccent{\bar}{x}}_{k}\|^{2}\leq
J^{*}(\boldsymbol{\underaccent{\bar}{x}}_{k})\leq
M_{J}\|\boldsymbol{\underaccent{\bar}{x}}_{k}\|^{2}\,,\quad\forall\boldsymbol{\underaccent{\bar}{x}}_{k}\in\Gamma_{\alpha}\,.$
(34)

In order to proceed with the stability analysis, the rate of convergence of
Algorithm 1 needs to be characterized. The next lemma shows linear convergence
of the sensitivity-based DMPC Algorithm.

######  Lemma 4.1.

Suppose that Assumption 4.1 holds. Then, there exists an upper bound on the
horizon length $T_{\max}>0$ such that for $T<T_{\max}$, Algorithm 1 converges
linearly, i.e.

$\displaystyle\bigg{\|}\begin{matrix}\boldsymbol{x}^{*}(\cdot;\boldsymbol{\underaccent{\bar}{x}}_{k})-\boldsymbol{x}^{q}(\cdot;\boldsymbol{\underaccent{\bar}{x}}_{k})\\\
\boldsymbol{\lambda}^{*}(\cdot;\boldsymbol{\underaccent{\bar}{x}}_{k})-\boldsymbol{\lambda}^{q}(\cdot;\boldsymbol{\underaccent{\bar}{x}}_{k})\\\
\end{matrix}\bigg{\|}_{L_{\infty}}\leq
p\bigg{\|}\begin{matrix}\boldsymbol{x}^{*}(\cdot;\boldsymbol{\underaccent{\bar}{x}}_{k})-\boldsymbol{x}^{q-1}(\cdot;\boldsymbol{\underaccent{\bar}{x}}_{k})\\\
\boldsymbol{\lambda}^{*}(\cdot;\boldsymbol{\underaccent{\bar}{x}}_{k})-\boldsymbol{\lambda}^{q-1}(\cdot;\boldsymbol{\underaccent{\bar}{x}}_{k})\\\
\end{matrix}\bigg{\|}_{L_{\infty}}$ (35)

for some $p\in(0,1)$.

######  Proof 4.2.

See Appendix B.

Lemma 4.1 reveals the fact that convergence can be guaranteed with a
sufficiently small prediction horizon $T$. This shows that there exists a
trade-off between the size of the domain of attraction $\Gamma_{\alpha}$ which
can be enlarged by increasing $T$, cf. Theorem 2.1, and the stability of the
distributed algorithm which can be achieved by decreasing $T$. It should be
pointed out that no convexity assumptions on the original OCP (3) are
necessary and convergence is explicitly shown without critical assumptions
which sets it apart from, e.g. ADMM.24

######  Remark 4.3.

In order to improve convergence and enlarge the maximum horizon
$T_{\mathrm{max}}$ the iterates of Algorithm 1 can be damped as first proposed
in 72, 63, i.e. by performing

$\displaystyle\boldsymbol{x}_{i}^{q}(\tau)\leftarrow(1-\epsilon)\boldsymbol{x}_{i}^{q}(\tau)+\epsilon\boldsymbol{x}_{i}^{q-1}(\tau)\,,\quad\boldsymbol{\lambda}_{i}^{q}(\tau)\leftarrow(1-\epsilon)\boldsymbol{\lambda}_{i}^{q}(\tau)+\epsilon\boldsymbol{\lambda}_{i}^{q-1}(\tau)\,,\quad\tau\in[0,\,T]$
(36)

as an intermediate step between Step 2 and 3 of Algorithm 1 with damping
factor $\epsilon\in[0,1)$. This measure seeks to prevent drastic changes
during the iterations and is therefore beneficial as the sensitivities can be
interpreted as a first-order approximation of the change in costs and has a
strong impact on the admissible horizon 63.

###  4.2 Convergence properties with stopping criterion

The investigation of stability of the DMPC scheme requires to look at the
error
$\Delta\boldsymbol{x}_{c}^{q_{k}}(\tau;\boldsymbol{\underaccent{\bar}{x}}_{k})$
in (33) as the difference between the suboptimal sensitivity-based and the
optimal (D)MPC solution. The following lemma bounds this error in terms of the
stopping criterion (28) for all agents $i\in\mathcal{V}$.

######  Lemma 4.4.

Suppose that Assumption 4.1 and $T<T_{\mathrm{max}}$ hold. Then, there exists
a constant $D>0$ such that the error (33) satisfies

$\|\Delta\boldsymbol{x}_{c}^{q_{k}}(\cdot;\boldsymbol{\underaccent{\bar}{x}}_{k})\|_{L_{\infty}}\leq
Dd\|\boldsymbol{\underaccent{\bar}{x}}_{k}\|\,,\quad\forall\boldsymbol{\underaccent{\bar}{x}}_{k}\in\Gamma_{\alpha}\,.$
(37)

######  Proof 4.5.

See Appendix C.

In the centralized optimal MPC case, the next sampling point lies on the
optimal state trajectory, i.e.
$\boldsymbol{\underaccent{\bar}{x}}_{k+1}=\boldsymbol{x}^{*}(\Delta
t;\boldsymbol{\underaccent{\bar}{x}}_{k})$ and the cost
$J^{*}(\boldsymbol{x}^{*}(\Delta t;\boldsymbol{\underaccent{\bar}{x}}_{k}))$
decreases according to (13). As indicated by (32), this is not the case for
the suboptimal distributed MPC scheme. The next lemma relates the cost
decrease in the suboptimal case to the optimal cost
$J^{*}(\boldsymbol{\underaccent{\bar}{x}}_{k})$ and the error (33).

######  Lemma 4.6.

Suppose that Assumptions 2.1 to 4.1 hold. Then, there exists some
$\bar{\alpha}<\alpha$ and correspondingly a subset
$\Gamma_{\bar{\alpha}}\subset\Gamma_{\alpha}$ such that
$\boldsymbol{x}_{c}^{q_{k}}(\Delta
t;\boldsymbol{\underaccent{\bar}{x}}_{k})\in\Gamma_{\alpha}$ for all
$\boldsymbol{\underaccent{\bar}{x}}_{k}\in\Gamma_{\bar{\alpha}}$. Moreover,
there exist constants $0<a\leq 1$ and $b,c>0$ such that the optimal cost at
the next sampling point $t_{k+1}$ satisfies

$J^{*}(\boldsymbol{x}_{c}^{q_{k}}(\Delta
t;\boldsymbol{\underaccent{\bar}{x}}_{k}))\leq(1-a)J^{*}(\boldsymbol{\underaccent{\bar}{x}}_{k})+b\sqrt{J^{*}(\boldsymbol{\underaccent{\bar}{x}}_{k})}\|\Delta\boldsymbol{x}_{c}^{q_{k}}(\cdot;\boldsymbol{\underaccent{\bar}{x}}_{k})\|_{L_{\infty}}+c\|\Delta\boldsymbol{x}_{c}^{q_{k}}(\cdot;\boldsymbol{\underaccent{\bar}{x}}_{k})\|_{L_{\infty}}^{2}\,,\quad\forall\boldsymbol{\underaccent{\bar}{x}}_{k}\in\Gamma_{\bar{\alpha}}\,.$
(38)

######  Proof 4.7.

See Appendix D.

Lemma 4.6 reveals that the cost decrease and overall MPC performance is
limited by the error
$\|\Delta\boldsymbol{x}_{c}^{q_{k}}(\cdot;\boldsymbol{\underaccent{\bar}{x}}_{k})\|_{L_{\infty}}$
which opposes the contraction term
$(1-a)J^{*}(\boldsymbol{\underaccent{\bar}{x}}_{k})$. Moreover, the original
domain of attraction $\Gamma_{\alpha}$ is reduced to the smaller set
$\Gamma_{\bar{\alpha}}$. Based on Lemmas 4.4 and 4.6, the stability of the
DMPC scheme under inexact minimization can be shown.

######  Theorem 4.8.

Suppose that Assumptions 2.1 through 4.1 as well as $T<T_{\mathrm{max}}$ hold
and let the constant $d$ in the stopping criterion (28) satisfy

$d<\frac{\sqrt{m_{J}}}{2cD}\bigg{(}\sqrt{b^{2}+4ac}-b\bigg{)}\,.$ (39)

Then, the optimal cost (9) and the error (33) decay exponentially for all
$\boldsymbol{\underaccent{\bar}{x}}_{0}\in\Gamma_{\bar{\alpha}}$ and the
origin of the closed-loop system under the DMPC control law (25) is
exponentially stable.

######  Proof 4.9.

The result of Lemma 4.4 can be expressed in terms of the optimal cost
$J^{*}(\boldsymbol{\underaccent{\bar}{x}}_{k})$ using the quadratic bound
(34), i.e.

$\|\Delta\boldsymbol{x}_{c}^{q_{k}}(\cdot;\boldsymbol{\underaccent{\bar}{x}}_{k})\|_{L_{\infty}}\leq\frac{Dd}{\sqrt{m_{J}}}\sqrt{J^{*}(\boldsymbol{\underaccent{\bar}{x}}_{k})}$
(40)

which, inserted in the relation (38) of Lemma 4.6, results in a bound on the
cost in the next MPC step

$\displaystyle J^{*}(\boldsymbol{x}_{c}^{q_{k}}(\Delta
t;\boldsymbol{\underaccent{\bar}{x}}_{k})\leq
p_{J}J^{*}(\boldsymbol{\underaccent{\bar}{x}}_{k})$ (41)

with $p_{J}:=(1-a)+\frac{bDd}{\sqrt{m_{J}}}+\frac{cD^{2}d^{2}}{m_{J}}$. By
(39), the contraction ratio satisfies $p_{J}<1$ which in return implies the
exponential decay of the optimal cost
$J^{*}(\boldsymbol{\underaccent{\bar}{x}}_{k})=J^{*}(\boldsymbol{x}_{c}^{q_{k}}(\Delta
t;\boldsymbol{x}_{k-1}))$

$\displaystyle J^{*}(\boldsymbol{\underaccent{\bar}{x}}_{k})\leq
p_{J}J^{*}(\boldsymbol{\underaccent{\bar}{x}}_{k-1})\leq(p_{J})^{k}J^{*}(\boldsymbol{\underaccent{\bar}{x}}_{0})\leq(p_{J})^{k}\bar{\alpha}$
(42)

and of the error
$\|\Delta\boldsymbol{x}_{c}^{q_{k}}(\cdot;\boldsymbol{\underaccent{\bar}{x}}_{k})\|_{L_{\infty}}$
by (40)

$\|\Delta\boldsymbol{x}_{c}^{q_{k}}(\cdot;\boldsymbol{\underaccent{\bar}{x}}_{k})\|_{L_{\infty}}\leq\frac{Dd}{\sqrt{m_{J}}}\sqrt{(p_{J})^{k}\bar{\alpha}}\,.$
(43)

To show exponential stability in continuous time, the state trajectory of the
closed loop, i.e.
$\boldsymbol{\underaccent{\bar}{x}}(t)=\boldsymbol{\underaccent{\bar}{x}}(t_{k}+\tau)=\boldsymbol{x}_{c}^{q_{k}}(\tau;\boldsymbol{\underaccent{\bar}{x}}_{k})$,
$\tau\in[0,\,\Delta t)$, $k\in\mathds{N}_{0}$ needs to be bounded by an
exponential envelope function. Using the triangle inequality, the bound (98),
and the bound (43), results in

$\displaystyle\|\boldsymbol{x}_{c}^{q_{k}}(\tau;\boldsymbol{\underaccent{\bar}{x}}_{k})\|$
$\displaystyle\leq\|\boldsymbol{x}^{*}(\tau;\boldsymbol{\underaccent{\bar}{x}}_{k})\|+\|\Delta\boldsymbol{x}^{q_{k}}_{c}(\tau;\boldsymbol{\underaccent{\bar}{x}}_{k})\|\leq\|\boldsymbol{\underaccent{\bar}{x}}_{k}\|\mathrm{e}^{\hat{L}\tau}+\|\Delta\boldsymbol{x}^{q_{k}}_{c}(\cdot;\boldsymbol{\underaccent{\bar}{x}}_{k})\|_{L_{\infty}}\leq(\mathrm{e}^{\hat{L}\tau}+Dd)\|\boldsymbol{\underaccent{\bar}{x}}_{k}\|$
$\displaystyle\leq(\mathrm{e}^{\hat{L}\tau}+Dd)\frac{1}{\sqrt{m_{J}}}\sqrt{J^{*}(\boldsymbol{\underaccent{\bar}{x}}_{k})}\leq(\mathrm{e}^{\hat{L}\tau}+Dd)\frac{1}{\sqrt{m_{J}}}(p_{J})^{\frac{1}{2}k}\sqrt{J^{*}(\boldsymbol{\underaccent{\bar}{x}}_{0})}$
(44)

which bounds the state trajectory in each MPC step $k$ for $\tau\in[0,\,\Delta
t)$. The exponential decay of (4.9) implies that there exists constants
$\gamma_{1},\,\gamma_{2}>0$ such that
$\|\boldsymbol{x}(t)\|\leq\gamma_{1}\mathrm{e}^{-\gamma_{2}t}$.

Theorem 4.8 shows that an explicit value for $d$ can be computed offline
ensuring stability and incremental improvement of the DMPC scheme. In
particular, it is evident that the value of $d>0$ controls the contraction
rate $p_{J}$ in (41) and the optimal (D)MPC case, i.e.
$p_{J}\rightarrow(1-a)$, is recovered for $d\rightarrow 0$. In fact, it is
possible to explicitly control the worst-case contraction rate $p_{J}$ in the
interval $(1-a,1)$ via $d$. However, calculating the value of $d$ via (39) is
usually too conservative for design purposes due to the various involved
Lipschitz and continuity estimates. Nevertheless, Theorem 4.8 states that a
stabilizing $d$ can always be found.

######  Remark 4.10.

Compared to the scheme given in the work24, we do not require to assume
boundedness and linear convergence of the underlying distributed optimization
algorithm to show stability of the MPC scheme, but explicitly prove that there
exists an upper bound on the horizon such that the iterates are bounded and
linear convergence of Algorithm 1 is ensured. The drawback is that the horizon
cannot be chosen arbitrarily long which can possibly decrease the performance
of the (D)MPC-controller. However, this can be mitigated by damping the
trajectories, cf. (36).

###  4.3 Convergence properties with a fixed number of iterations

In the following, it is discussed how the stopping criterion (28) can be
equivalently fulfilled with a finite number of iterations $q$ without
explicitly evaluating it. The following theorem precisely states a lower bound
on the required iterations $q_{\mathrm{max}}$.

######  Theorem 4.11.

Suppose that Assumptions 2.1 through 4.1 as well as $T<T_{\mathrm{max}}$ hold
and $d$ is chosen according to Theorem 4.8. Then, for all
$\boldsymbol{\underaccent{\bar}{x}}_{0}\in\Gamma_{\bar{\alpha}}$ there exists
a constant $e>0$ such that if the number of iterations per MPC step and the
initial optimization error satisfy

$q_{\mathrm{max}}>1+\log_{p}\bigg{(}\frac{d}{e(1+p)}\bigg{)}\,,\quad\bigg{\|}\begin{matrix}\boldsymbol{x}^{1}(\cdot;\boldsymbol{\underaccent{\bar}{x}}_{0})-\boldsymbol{x}^{0}(\cdot;\boldsymbol{\underaccent{\bar}{x}}_{0})\\\
\boldsymbol{\lambda}^{1}(\cdot;\boldsymbol{\underaccent{\bar}{x}}_{0})-\boldsymbol{\lambda}^{0}(\cdot;\boldsymbol{\underaccent{\bar}{x}}_{0})\\\
\end{matrix}\bigg{\|}_{L_{\infty}}\leq
dN\sqrt{\frac{\bar{\alpha}}{m_{J}}}p^{1-q_{\mathrm{max}}}\,,$ (45)

then the origin of the closed-loop system under the DMPC control law (25) is
exponentially stable and the optimization error (33) decays exponentially.

######  Proof 4.12.

The proof consists of showing by induction that the stopping criterion, i.e.

$\bigg{\|}\begin{matrix}\delta\boldsymbol{x}^{q}(\cdot;\boldsymbol{\underaccent{\bar}{x}}_{k})\\\
\delta\boldsymbol{\lambda}^{q}(\cdot;\boldsymbol{\underaccent{\bar}{x}}_{k})\end{matrix}\bigg{\|}_{L_{\infty}}\leq
d\|\boldsymbol{\underaccent{\bar}{x}}_{k}\|\,,$ (46)

with
$\delta\boldsymbol{x}^{q}(\tau;\boldsymbol{\underaccent{\bar}{x}}_{k}):=\boldsymbol{x}^{q}(\tau;\boldsymbol{\underaccent{\bar}{x}}_{k})-\boldsymbol{x}^{q-1}(\tau;\boldsymbol{\underaccent{\bar}{x}}_{k})$
and
$\delta\boldsymbol{\lambda}^{q}(\tau;\boldsymbol{\underaccent{\bar}{x}}_{k}):=\boldsymbol{\lambda}^{q}(\tau;\boldsymbol{\underaccent{\bar}{x}}_{k})-\boldsymbol{\lambda}^{q-1}(\tau;\boldsymbol{\underaccent{\bar}{x}}_{k})$
holds in each MPC step $k$ for the given conditions (45). Consequently, the
optimal cost $J^{*}(\boldsymbol{\underaccent{\bar}{x}}_{k})$ and optimization
error
$\|\Delta\boldsymbol{x}_{c}^{q_{k}}(\cdot;\boldsymbol{\underaccent{\bar}{x}}_{k})\|_{L_{\infty}}$
decrease according to (42) and (43) as stated by Theorem 4.8. Note that the
linear convergence property (35) can equivalently be expressed as

$\bigg{\|}\begin{matrix}\delta\boldsymbol{x}^{q}(\cdot;\boldsymbol{\underaccent{\bar}{x}}_{k})\\\
\delta\boldsymbol{\lambda}^{q}(\cdot;\boldsymbol{\underaccent{\bar}{x}}_{k})\end{matrix}\bigg{\|}_{L_{\infty}}\leq
p^{q-1}\bigg{\|}\begin{matrix}\delta\boldsymbol{x}^{1}(\cdot;\boldsymbol{\underaccent{\bar}{x}}_{k})\\\
\delta\boldsymbol{\lambda}^{1}(\cdot;\boldsymbol{\underaccent{\bar}{x}}_{k})\end{matrix}\bigg{\|}_{L_{\infty}}.$
(47)

The relation (46) must hold at $k=0$ (induction start), i.e.

$\bigg{\|}\begin{matrix}\delta\boldsymbol{x}^{q}(\cdot;\boldsymbol{\underaccent{\bar}{x}}_{0})\\\
\delta\boldsymbol{\lambda}^{q}(\cdot;\boldsymbol{\underaccent{\bar}{x}}_{0})\end{matrix}\bigg{\|}_{L_{\infty}}\leq
p^{q-1}\bigg{\|}\begin{matrix}\delta\boldsymbol{x}^{1}(\cdot;\boldsymbol{\underaccent{\bar}{x}}_{0})\\\
\delta\boldsymbol{\lambda}^{1}(\cdot;\boldsymbol{\underaccent{\bar}{x}}_{0})\end{matrix}\bigg{\|}_{L_{\infty}}\leq
Nd\|\boldsymbol{\underaccent{\bar}{x}}_{0}\|$ (48)

which in return shows that the initial optimization error must satisfy

$\bigg{\|}\begin{matrix}\delta\boldsymbol{x}^{1}(\cdot;\boldsymbol{\underaccent{\bar}{x}}_{0})\\\
\delta\boldsymbol{\lambda}^{1}(\cdot;\boldsymbol{\underaccent{\bar}{x}}_{0})\end{matrix}\bigg{\|}_{L_{\infty}}\leq
dN\sqrt{\frac{\bar{\alpha}}{m_{J}}}p^{1-q_{\mathrm{max}}}$ (49)

which is true by assumption in Theorem 4.11. Following the lines of Theorem
4.8, this implies that
$\|\Delta\boldsymbol{x}_{c}^{q_{k}}(\cdot;\boldsymbol{\underaccent{\bar}{x}}_{0})\|_{L_{\infty}}\leq
Dd\|\boldsymbol{\underaccent{\bar}{x}}_{0}\|\leq
d\frac{D}{\sqrt{m_{J}}}\sqrt{\bar{\alpha}}$ which in return shows that
$J^{*}(\boldsymbol{\underaccent{\bar}{x}}_{1})\leq
p_{J}J^{*}(\boldsymbol{\underaccent{\bar}{x}}_{0})$ completing the induction
start. We continue with the induction hypothesis that (46) holds at MPC step
$k$. At MPC step $k+1$ the following bound can be given via the linear
convergence property (47) and Minkowski’s inequality

$\displaystyle\bigg{\|}\begin{matrix}\delta\boldsymbol{x}^{q}(\cdot;\boldsymbol{\underaccent{\bar}{x}}_{k+1})\\\
\delta\boldsymbol{\lambda}^{q}(\cdot;\boldsymbol{\underaccent{\bar}{x}}_{k+1})\end{matrix}\bigg{\|}_{L_{\infty}}\leq\bigg{\|}\begin{matrix}\boldsymbol{x}^{q}(\cdot;\boldsymbol{\underaccent{\bar}{x}}_{k+1})-\boldsymbol{x}^{*}(\cdot;\boldsymbol{\underaccent{\bar}{x}}_{k+1})\\\
\boldsymbol{\lambda}^{q}(\cdot;\boldsymbol{\underaccent{\bar}{x}}_{k+1})-\boldsymbol{\lambda}^{*}(\cdot;\boldsymbol{\underaccent{\bar}{x}}_{k+1})\end{matrix}\bigg{\|}_{L_{\infty}}+\bigg{\|}\begin{matrix}\boldsymbol{x}^{q-1}(\cdot;\boldsymbol{\underaccent{\bar}{x}}_{k+1})-\boldsymbol{x}^{*}(\cdot;\boldsymbol{\underaccent{\bar}{x}}_{k+1})\\\
\boldsymbol{\lambda}^{q-1}(\cdot;\boldsymbol{\underaccent{\bar}{x}}_{k+1})-\boldsymbol{\lambda}^{*}(\cdot;\boldsymbol{\underaccent{\bar}{x}}_{k+1})\end{matrix}\bigg{\|}_{L_{\infty}}$
$\displaystyle\leq(1+p)p^{q-1}\bigg{\|}\begin{matrix}\boldsymbol{x}^{0}(\cdot;\boldsymbol{\underaccent{\bar}{x}}_{k+1})-\boldsymbol{x}^{*}(\cdot;\boldsymbol{\underaccent{\bar}{x}}_{k+1})\\\
\boldsymbol{\lambda}^{0}(\cdot;\boldsymbol{\underaccent{\bar}{x}}_{k+1})-\boldsymbol{\lambda}^{*}(\cdot;\boldsymbol{\underaccent{\bar}{x}}_{k+1})\end{matrix}\bigg{\|}_{L_{\infty}}\leq(1+p)p^{q-1}\bigg{(}\bigg{\|}\begin{matrix}\boldsymbol{x}^{0}(\cdot;\boldsymbol{\underaccent{\bar}{x}}_{k+1})-\boldsymbol{x}^{*}(\cdot;\boldsymbol{\underaccent{\bar}{x}}_{k})\\\
\boldsymbol{\lambda}^{0}(\cdot;\boldsymbol{\underaccent{\bar}{x}}_{k+1})-\boldsymbol{\lambda}^{*}(\cdot;\boldsymbol{\underaccent{\bar}{x}}_{k})\end{matrix}\bigg{\|}_{L_{\infty}}+\bigg{\|}\begin{matrix}\boldsymbol{x}^{*}(\cdot;\boldsymbol{\underaccent{\bar}{x}}_{k})-\boldsymbol{x}^{*}(\cdot;\boldsymbol{\underaccent{\bar}{x}}_{k+1})\\\
\boldsymbol{\lambda}^{*}(\cdot;\boldsymbol{\underaccent{\bar}{x}}_{k})-\boldsymbol{\lambda}^{*}(\cdot;\boldsymbol{\underaccent{\bar}{x}}_{k+1})\end{matrix}\bigg{\|}_{L_{\infty}}\bigg{)}.$
(50)

We now bound the first norm on the right hand side of (4.12) by using the re-
initialization for the adjoint state
$\boldsymbol{\lambda}^{0}(\tau;\boldsymbol{\underaccent{\bar}{x}}_{k+1})=\boldsymbol{\lambda}^{q}(\tau;\boldsymbol{\underaccent{\bar}{x}}_{k})$
and state
$\boldsymbol{x}^{0}(\tau;\boldsymbol{\underaccent{\bar}{x}}_{k+1})=\boldsymbol{x}^{q}(\tau;\boldsymbol{\underaccent{\bar}{x}}_{k})$
as follows

$\displaystyle\bigg{\|}\begin{matrix}\boldsymbol{x}^{0}(\cdot;\boldsymbol{\underaccent{\bar}{x}}_{k+1})-\boldsymbol{x}^{*}(\cdot;\boldsymbol{\underaccent{\bar}{x}}_{k})\\\
\boldsymbol{\lambda}^{0}(\cdot;\boldsymbol{\underaccent{\bar}{x}}_{k+1})-\boldsymbol{\lambda}^{*}(\cdot;\boldsymbol{\underaccent{\bar}{x}}_{k})\end{matrix}\bigg{\|}_{L_{\infty}}=\bigg{\|}\begin{matrix}\boldsymbol{x}^{q}(\cdot;\boldsymbol{\underaccent{\bar}{x}}_{k})-\boldsymbol{x}^{*}(\cdot;\boldsymbol{\underaccent{\bar}{x}}_{k})\\\
\boldsymbol{\lambda}^{q}(\cdot;\boldsymbol{\underaccent{\bar}{x}}_{k})-\boldsymbol{\lambda}^{*}(\cdot;\boldsymbol{\underaccent{\bar}{x}}_{k})\end{matrix}\bigg{\|}_{L_{\infty}}\leq\frac{p}{1-p}\bigg{\|}\begin{matrix}\boldsymbol{x}^{q}(\cdot;\boldsymbol{\underaccent{\bar}{x}}_{k})-\boldsymbol{x}^{q-1}(\cdot;\boldsymbol{\underaccent{\bar}{x}}_{k})\\\
\boldsymbol{\lambda}^{q}(\cdot;\boldsymbol{\underaccent{\bar}{x}}_{k})-\boldsymbol{\lambda}^{q-1}(\cdot;\boldsymbol{\underaccent{\bar}{x}}_{k})\end{matrix}\bigg{\|}_{L_{\infty}}\leq\frac{p}{1-p}Nd\|\boldsymbol{\underaccent{\bar}{x}}_{k}\|\,,$
(51)

where the induction hypothesis (46) was used in the last inequality of (51).
Regarding the second norm on the right hand side of (4.12), Equations (101)
and (Sensitivity-Based Distributed Model Predictive Control for Nonlinear
Systems under Inexact Optimization) in Appendix E show that there exists
Lipschitz constants $L_{x},L_{\lambda}>0$ such that
$\|\boldsymbol{x}^{*}(\cdot;\boldsymbol{\underaccent{\bar}{x}}_{k+1})-\boldsymbol{x}^{*}(\cdot;\boldsymbol{\underaccent{\bar}{x}}_{k})\|_{L_{\infty}}\leq
L_{x}\|\boldsymbol{\underaccent{\bar}{x}}_{k+1}-\boldsymbol{\underaccent{\bar}{x}}_{k}\|$
and
$\|\boldsymbol{\lambda}^{*}(\cdot;\boldsymbol{\underaccent{\bar}{x}}_{k+1})-\boldsymbol{\lambda}^{*}(\cdot;\boldsymbol{\underaccent{\bar}{x}}_{k})\|_{L_{\infty}}\leq
L_{\lambda}\|\boldsymbol{\underaccent{\bar}{x}}_{k+1}-\boldsymbol{\underaccent{\bar}{x}}_{k}\|$
for all $\boldsymbol{\underaccent{\bar}{x}}_{k}\in\Gamma_{\alpha}$,
$k\in\mathbb{N}_{0}$ (recall that
$\boldsymbol{\underaccent{\bar}{x}}_{k+1}\in\Gamma_{\alpha}$ due to Lemma 4.6)
with which the norm can be bounded as

$\displaystyle\bigg{\|}\begin{matrix}\boldsymbol{x}^{*}(\cdot;\boldsymbol{\underaccent{\bar}{x}}_{k})-\boldsymbol{x}^{*}(\cdot;\boldsymbol{\underaccent{\bar}{x}}_{k+1})\\\
\boldsymbol{\lambda}^{*}(\cdot;\boldsymbol{\underaccent{\bar}{x}}_{k})-\boldsymbol{\lambda}^{*}(\cdot;\boldsymbol{\underaccent{\bar}{x}}_{k+1})\end{matrix}\bigg{\|}_{L_{\infty}}\leq\|\boldsymbol{x}^{*}(\cdot;\boldsymbol{\underaccent{\bar}{x}}_{k})-\boldsymbol{x}^{*}(\cdot;\boldsymbol{\underaccent{\bar}{x}}_{k+1})\|_{L_{\infty}}+\|\boldsymbol{\lambda}^{*}(\cdot;\boldsymbol{\underaccent{\bar}{x}}_{k})-\boldsymbol{\lambda}^{*}(\cdot;\boldsymbol{\underaccent{\bar}{x}}_{k+1})\|_{L_{\infty}}\leq(L_{\lambda}+L_{x})\|\boldsymbol{\underaccent{\bar}{x}}_{k+1}-\boldsymbol{\underaccent{\bar}{x}}_{k}\|.$
(52)

The norm of the current system state $\boldsymbol{\underaccent{\bar}{x}}_{k}$
can be related to the next state
$\boldsymbol{\underaccent{\bar}{x}}_{k+1}=\boldsymbol{x}_{c}^{q_{k}}(\Delta
t;\boldsymbol{\underaccent{\bar}{x}}_{k})$ via the error (33)

$\displaystyle\|\boldsymbol{x}_{c}^{q_{k}}(\Delta
t;\boldsymbol{\underaccent{\bar}{x}}_{k})\|=\|\boldsymbol{x}^{*}(\Delta
t;\boldsymbol{\underaccent{\bar}{x}}_{k})+\Delta\boldsymbol{x}_{c}^{q_{k}}(\Delta
t;\boldsymbol{\underaccent{\bar}{x}}_{k})\|\geq\|\boldsymbol{x}^{*}(\Delta
t;\boldsymbol{\underaccent{\bar}{x}}_{k})\|-\|\Delta\boldsymbol{x}_{c}^{q_{k}}(\Delta
t;\boldsymbol{\underaccent{\bar}{x}}_{k})\|\geq(\mathrm{e}^{-\hat{L}\Delta
t}-Dd)\|\boldsymbol{\underaccent{\bar}{x}}_{k}\|\,,$ (53)

where the reverse triangle inequality together with (99) and (37) was used. If
$d\leq\mathrm{e}^{-\hat{L}\Delta t}/D$, then $(\mathrm{e}^{-\hat{L}\Delta
t}-Dd)>0$ and we can upper bound the current state
$\|\boldsymbol{\underaccent{\bar}{x}}_{k}\|$ by the next state
$\|\boldsymbol{\underaccent{\bar}{x}}_{k+1}\|$ according to

$\|\boldsymbol{\underaccent{\bar}{x}}_{k}\|\leq\frac{1}{\mathrm{e}^{-\hat{L}\Delta
t}-Dd}\|\boldsymbol{\underaccent{\bar}{x}}_{k+1}\|.$ (54)

Utilizing (54) to bound the current state
$\|\boldsymbol{\underaccent{\bar}{x}}_{k}\|$ in (52) and (51) via the next
state $\boldsymbol{\underaccent{\bar}{x}}_{k+1}$ leads to

$\displaystyle\bigg{\|}\begin{matrix}\boldsymbol{x}^{q}(\cdot;\boldsymbol{\underaccent{\bar}{x}}_{k})-\boldsymbol{x}^{*}(\cdot;\boldsymbol{\underaccent{\bar}{x}}_{k})\\\
\boldsymbol{\lambda}^{q}(\cdot;\boldsymbol{\underaccent{\bar}{x}}_{k})-\boldsymbol{\lambda}^{*}(\cdot;\boldsymbol{\underaccent{\bar}{x}}_{k})\end{matrix}\bigg{\|}_{L_{\infty}}\leq
c_{1}\|\boldsymbol{\underaccent{\bar}{x}}_{k+1}\|\,,\quad\bigg{\|}\begin{matrix}\boldsymbol{x}^{*}(\cdot;\boldsymbol{\underaccent{\bar}{x}}_{k})-\boldsymbol{x}^{*}(\cdot;\boldsymbol{\underaccent{\bar}{x}}_{k+1})\\\
\boldsymbol{\lambda}^{*}(\cdot;\boldsymbol{\underaccent{\bar}{x}}_{k})-\boldsymbol{\lambda}^{*}(\cdot;\boldsymbol{\underaccent{\bar}{x}}_{k+1})\end{matrix}\bigg{\|}_{L_{\infty}}\leq
c_{2}\|\boldsymbol{\underaccent{\bar}{x}}_{k+1}\|$ (55)

with $c_{1}:=\frac{Ndp}{(1-p)(\mathrm{e}^{-\hat{L}\Delta t}-Dd)}>0$ and
$c_{2}:=(L_{\lambda}+\bar{L}\mathrm{e}^{\hat{L}T})(1+\frac{1}{\mathrm{e}^{-\hat{L}\Delta
t}-Dd})>0$. Inserting these bounds into (4.12) results in

$\displaystyle\bigg{\|}\begin{matrix}\delta\boldsymbol{x}^{q}(\cdot;\boldsymbol{\underaccent{\bar}{x}}_{k+1})\\\
\delta\boldsymbol{\lambda}^{q}(\cdot;\boldsymbol{\underaccent{\bar}{x}}_{k+1})\end{matrix}\bigg{\|}_{L_{\infty}}\leq(1+p)(c_{1}+c_{2})p^{q-1}\|\boldsymbol{\underaccent{\bar}{x}}_{k+1}\|.$
(56)

Utilizing the condition on the iteration number (45) with $e:=c_{1}+c_{2}>0$,
the bound (4.12) eventually becomes

$\displaystyle\bigg{\|}\begin{matrix}\delta\boldsymbol{x}^{q}(\cdot;\boldsymbol{\underaccent{\bar}{x}}_{k+1})\\\
\delta\boldsymbol{\lambda}^{q}(\cdot;\boldsymbol{\underaccent{\bar}{x}}_{k+1})\end{matrix}\bigg{\|}_{L_{\infty}}\leq
d\|\boldsymbol{\underaccent{\bar}{x}}_{k+1}\|\,,$ (57)

which completes the induction step. Following the lines of Theorem 4.8, this
implies the exponential decay of the error, i.e.
$\|\Delta\boldsymbol{x}_{c}^{q_{k}}(\cdot;\boldsymbol{\underaccent{\bar}{x}}_{k})\|_{L_{\infty}}\leq
d\frac{D}{\sqrt{m_{J}}}\sqrt{p_{J}^{k}\bar{\alpha}}$, which in return shows
the exponential decay of the cost, i.e.
$J^{*}(\boldsymbol{\underaccent{\bar}{x}}_{k})\leq p_{J}^{k}\bar{\alpha}$.

Theorem 4.11 shows that there exists a sufficiently large number of iterations
$q_{\mathrm{max}}$ such that the stopping criterion (28) is fulfilled in every
MPC iteration, thus, guaranteeing exponential stability of the DMPC scheme.
Executing Algorithm 1 with a fixed number of iterations per MPC step has the
advantage that the computation time on the agent level and the communication
steps are kept constant throughout the MPC steps which is critical in real-
time applications.

######  Remark 4.13.

The number of iterations $q_{k}$ for satisfying the stopping criterion (28)
will be lower in practice since (45) is a worst-case estimate for the required
number of iterations. Typically, the required number of iterations is highest
in the first MPC step $k=0$ due to the "cold" start with some initial
trajectories and decreases with the progression of the MPC scheme to the warm
start (23), see also the numerical example in Section 5.1.

##  5 Numerical Evaluation

Two numerical examples are utilized to demonstrate the performance of
Algorithm 1 and its theoretical properties. At first, a nonlinear system of
three agents is utilized to compare the distributed to the central solution
and to verify the linear convergence property. In addition, a comparison to
the ADMM algorithm46 is presented. The next example concerns the distributed
end region and the effect of the suboptimal predicted trajectories on the
closed-loop behavior. Algorithm 1 is implemented in C++ within the modular
DMPC framework GRAMPC-D26, 37 in which the local OCPs (19) are solved with the
toolbox GRAMPC85 via the projected gradient method.

###  5.1 Benchmark system

In this section, a typical benchmark system which is often used in nonlinear
DMPC 28, 24, 86, 87, 88 is considered for the application of the sensitivity-
based DMPC scheme. The dynamics of the system describing a bipedal locomotor
experiment are given by the following nonlinear differential equations

$\displaystyle\ddot{\theta}_{1}$
$\displaystyle=0.1(1-5.25\theta_{1}^{2})\dot{\theta}_{1}-\theta_{1}+u_{1}$
$\displaystyle\ddot{\theta}_{2}$
$\displaystyle=0.001(1-6070\theta_{2}^{2})\dot{\theta}_{2}-4\theta_{2}+0.057\theta_{1}\dot{\theta}_{1}+0.1(\dot{\theta}_{2}-\dot{\theta}_{3})+u_{2}$
(58) $\displaystyle\ddot{\theta}_{3}$
$\displaystyle=0.001(1-192\theta_{3}^{2})\dot{\theta}_{3}-4\theta_{3}+0.057\theta_{1}\dot{\theta}_{1}+0.1(\dot{\theta}_{3}-\dot{\theta}_{2})+u_{3}$

with the states $\boldsymbol{x}_{i}=[\theta_{i},\dot{\theta}_{i}]$ and
controls $u_{i}$ of each agent $i\in\mathcal{V}=\\{1,2,3\\}$ for which the
system (58) can be represented in the neighbor-affine form (1). The coupling
structure of the system follows as $\mathcal{N}_{1}^{\rightarrow}=\\{2,3\\}$,
$\mathcal{N}_{2}^{\rightarrow}=\\{3\\}$,
$\mathcal{N}_{3}^{\rightarrow}=\\{2\\}$ and
$\mathcal{N}_{1}^{\leftarrow}=\\{\\}$,
$\mathcal{N}_{2}^{\leftarrow}=\\{1,3\\}$,
$\mathcal{N}_{3}^{\leftarrow}=\\{1,2\\}$.

The controls are constrained to the set $u_{i}(t)\in[-1,1]$, $\forall
i\in\mathcal{V}$ and the control task is to drive the system from its stable
limit cycle to the origin. To this end, the quadratic cost functions
$l_{i}(\boldsymbol{x}_{i},u_{i})=\boldsymbol{x}_{i}^{\mathrm{T}}\boldsymbol{Q}_{i}\boldsymbol{x}_{i}+R_{i}u_{i}^{2},\,V_{i}(\boldsymbol{x}_{i})=\boldsymbol{x}_{i}^{\mathrm{T}}\boldsymbol{P}_{i}\boldsymbol{x}_{i}$
with the weighting matrices $\boldsymbol{Q}_{i}=\mathrm{diag}(30,30)$ and
$R_{i}=0.1$ are employed for each agent. We follow the approach of Section 2.2
and linearize the system around the origin to obtain the linear system with
$\boldsymbol{A}=\partial_{\boldsymbol{x}}\boldsymbol{f}(\boldsymbol{x},\boldsymbol{u})$
and
$\boldsymbol{B}=\partial_{\boldsymbol{u}}\boldsymbol{f}(\boldsymbol{x},\boldsymbol{u})$
which is linearly controllable but unstable. The optimization procedure (18)
with $\gamma=1.2$ is utilized offline to find the terminal weights

$\boldsymbol{P}_{1}=\begin{bmatrix}37.4&2.0\\\
2.0&2.2\end{bmatrix}\,,\quad\boldsymbol{P}_{2}=\boldsymbol{P}_{3}=\begin{bmatrix}38.8&1.7\\\
1.7&2.2\end{bmatrix}$ (59)

and structured feedback controller $\boldsymbol{K}_{11}=[-16.3,-18.3]$,
$\boldsymbol{K}_{12}=\boldsymbol{K}_{13}=\boldsymbol{0}$ and
$\boldsymbol{K}_{22}=\boldsymbol{K}_{33}=[-13.8,\,-15.3]$,
$\boldsymbol{K}_{23}=K_{32}=[0.02,\,0.04]$ and
$\boldsymbol{K}_{21}=\boldsymbol{K}_{31}=\boldsymbol{0}$. Both the conditions
that the linear controller satisfies the input constraints, i.e.
$\underaccent{\bar}{u}_{i}(t)\in[-1,\,1]$, and that the CLF inequality (11)
are fulfilled in the terminal region
$\Omega_{\beta}=\\{\boldsymbol{x}^{\mathrm{T}}\boldsymbol{P}\boldsymbol{x}\leq
0.9\\}$. Finally, the prediction horizon is set to $T=3.0\,\unit{}$ and the
sampling time to $\Delta t=0.05\unit{}$, respectively.

The initial values for the (D)MPC simulation illustrated in Figure 1 are
chosen as $\theta_{1,0}=0.7$, $\theta_{2,0}=0.28$, $\theta_{3,0}=-0.61$ and
$\dot{\theta}_{i,0}=0$ according to reference28, for all $i\in\mathcal{V}$ and
the system is simulated for a total of $6\unit{}$. The top two figures show a
comparison of the states and controls for the centralized optimal solution and
the solution obtained by Algorithm 1 with $d=0.1$ in the stopping criterion
(28). Furthermore, the lower left figure shows the exponential decay of the
cost function
$J(\boldsymbol{\underaccent{\bar}{x}}_{k})=\sum_{i\in\mathcal{V}}J_{i}(\boldsymbol{u}_{i};\boldsymbol{\underaccent{\bar}{x}}_{k})$
in each MPC step which is in line with Theorems 2.1 and 4.8. The required
iterations $q_{k}$ such that the stopping criterion (28) is fulfilled are
depicted in the lower right figure for different values of $d$. Clearly, a
tighter stopping bound leads to more required iterations. As indicated by
Theorem 4.11 the required iterations are highest in the initial MPC step due
to the large initial optimization error and then converge to a stationary
value. This is especially prevalent in the case of $d=0.1$ where four
algorithm iterations are required to overcome the initial optimization error
and then subsequently drop to two iterations for the remaining MPC steps due
to the warm-start (23).

The overall communication effort can be characterized by the number of
trajectories sent between an agent $i$ and its neighbor $j\in\mathcal{N}_{i}$
in each MPC sampling step. Consequently, Algorithm 1 requires a total number
of $q_{k}n_{i}(|\mathcal{N}_{i}^{\leftarrow}|+|\mathcal{N}_{i}|)$ trajectories
to be communicated by each agent. In this example, the first MPC step for
$d=0.1$ thus requires $80$ trajectories to be exchanged in the communication
network and is halved to $40$ trajectories in the following steps. In
practice, the trajectories are transmitted in discretized form and the actual
data amount depends on the number of discretization points (21 in this
example). The major advantage however is that only one local communication
step is needed in each iteration $q$ which greatly reduces the overall
communication time in the network.37

Figure 1: Closed-loop optimal and distributed solution of System (58).

Figure 2 shows the evolution of the envelope and mean norm between optimal and
suboptimal trajectories in each step $q$ of Algorithm 1 in the first MPC step
for different initial conditions. According to Lemma 4.1, linear convergence
in each step is guaranteed as long as the prediction horizon is chosen
sufficiently small which is confirmed by this example system. It is also
apparent that the adjoint states are the limiting factor for convergence as
their progress is slower than that of the states which also can be seen in the
proof of Lemma 4.1. This is mainly due to the fact that initialization is
further away from the optimal solution than in the case of the states and that
the "negotiation" between agents takes place via the adjoint states.

Figure 2: Envelope over all sampled trajectories together with the median
convergence curve of Algorithm 1 for random initial conditions in the first
DMPC step $k=0$.

In addition, a comparison with the popular ADMM algorithm 46 is given. In
particular, we utilize the ADMM implementation of GRAMPC-D26 with penalty
parameter adaption and employ the stabilizing stopping criterion from 24 which
is similar to (28) with $d=5\times 10^{-3}$ and yields approximately the same
control performance. The goal of stabilizing the coupled oscillators (58)
remains the same. However, to ensure maximum comparability, we use the same
cost matrices and initial values as in 24 to arrive at an identical setup. The
convergence behavior is evaluated by observing the normalized cost progression
w.r.t. the number of algorithm iterations in the first MPC step $k=0$ which is
visualized in the left of Figure 3. The sensitivity-based algorithm takes
about 4 iterations to convergence while the ADMM algorithm requires
approximately 30 iterations to converge towards the central solution.
Furthermore, the number algorithm iterations $q_{k}$ to achieve a stabilizing
solution is depicted on the right of Figure 3. Clearly, the required
iterations $q_{k}$ in each MPC-step of the proposed algorithm are lower
compared to ADMM underpinning the promising application possibilities of the
sensitivity-based algorithm to DMPC.

Figure 3: Comparison of sensitivity-based and ADMM-based DMPC in terms of the
convergence behavior in the first DMPC step (left) as well as the required
iterations for a stabilizing solution (right).

###  5.2 Distributed end region

The next example concerns the distributed end region and investigates how the
suboptimality of the predicted trajectories affect the region of attraction.
As an example we use a similar system as in the work75

$\dot{x}_{i}=(\mu_{i}+(1-\mu_{i})x_{i})u_{i}+\sum_{j\in\mathcal{N}_{i}^{\leftarrow}}\epsilon_{ij}x_{j}\,,$
(60)

where a number of $N=2$ agents are considered which are coupled bi-
directionally resulting in $\mathcal{N}_{1}=\\{2\\}$ and
$\mathcal{N}_{2}=\\{1\\}$. The controls are constrained to the set
$u_{i}\in[-2,\,2]$ for all $i\in\mathcal{V}$. Again the quadratic functions
(15) are chosen with the weighting matrices
$\boldsymbol{Q}=\mathrm{diag}(10,10)$ and $\boldsymbol{R}=\mathrm{diag}(1,1)$.
For now set $\epsilon=\epsilon_{12}=\epsilon_{21}=2$ and
$\mu_{1}=\mu_{2}=\mu=0.5$. The optimization problem (18) is solved with
$\gamma=1.1$ with and without the structural constraint (18c) to obtain a
separable
$V(\boldsymbol{x})=\boldsymbol{x}^{\mathrm{T}}\boldsymbol{P}_{d}\boldsymbol{x}$
and non-separable terminal cost function
$V(\boldsymbol{x})=\boldsymbol{x}^{\mathrm{T}}\boldsymbol{P}_{c}\boldsymbol{x}$.
Consequently, the invariant end region
$\Omega_{\beta}=\\{\boldsymbol{x}^{\mathrm{T}}\boldsymbol{P}_{c/d}\boldsymbol{x}\leq
0.9\\}$ is obtained for both terminal cost functions since the dynamics are
identical for both agents. The left plot of Figure 4 visualizes the two
terminal regions. Clearly, the size of $\Omega_{\beta}$ is reduced in the
distributed setting and the separable terminal region is contained within the
non-separable end region. The right plot shows the change in area of the
terminal region which is proportional to $\det(\boldsymbol{P}^{-1})$ with
respect to an increasing coupling strength $\epsilon$ for different $\mu$. The
intuitive result that the separable terminal region is reduced compared to the
non-separable one with an increase in coupling strength is visible. In
addition, an increased nonlinearity (for $\mu=1$ system (60) is linear)
worsens this effect and reduces the terminal region even further.

Figure 4: Comparison of the end regions following from separable and non-
separable terminal cost functions.

The unstable system (60) is now stabilized via the DMPC control scheme for
different initial conditions. The coupling strengths are set to
$\epsilon_{12}=0.5$ and $\epsilon_{21}=2.0$, $\mu_{1}=1.0$ and $\mu_{2}=0.5$,
the prediction horizon to $T=0.5\,\unit{}$ and the sampling to $\Delta
t=0.05\,\unit{}$. The terminal cost is again designed via the procedure of
Section 2.2 to obtain the end region
$\Omega_{\beta}=\\{\boldsymbol{x}^{\mathrm{T}}\boldsymbol{P}\boldsymbol{x}\leq
1.05\\}$. Figure 5 shows the optimal predictions obtained by solving OCP (3),
the suboptimal prediction of Algorithm 1 after $q=1$ iteration in which OCP
(19) is solved once and the resulting closed-loop trajectory of the suboptimal
DMPC controller with $q_{\mathrm{max}}=1$ as well as the terminal region
$\Omega_{\beta}$. All initial conditions where chosen to be located within the
region of attraction $\Gamma_{\alpha}$ since all endpoints of the predicted
trajectories reach the terminal region, i.e.
$\boldsymbol{x}^{*}(T;\boldsymbol{x}_{0})\in\Omega_{\beta}$. Thus, with the
chosen initial conditions the optimal MPC controller is guaranteed to be
stable according to Theorem 2.1. This, however, is not necessarily the case
for the suboptimal distributed controller.

Note the trajectories for the upper left initial condition at
$\boldsymbol{x}_{0}=[-1.3\,1.4]^{\mathrm{T}}$. They visualize the difference
between the original domain of attraction $\Gamma_{\alpha}$ and the reduced
domain of attraction $\Gamma_{\bar{\alpha}}$ as the next sampling point of the
suboptimal closed-loop trajectory $\boldsymbol{x}_{c}(\Delta
t;\underaccent{\bar}{\vm x}_{0})$ does not lie within the region of attraction
$\Gamma_{\alpha}$ such that Theorems 4.8 and 4.11 hold. This can be seen from
the fact that the endpoint of the optimal predicted trajectory at the next
sampling point $\underaccent{\bar}{\vm x}_{1}=\boldsymbol{x}_{c}(\Delta
t;\underaccent{\bar}{\vm x}_{0})$ does not lie within the terminal region,
i.e. $\boldsymbol{x}^{*}(T;\underaccent{\bar}{\vm
x}_{1})\notin\Omega_{\beta}$.

Thus, the DMPC controller with this number of algorithm iterations is not
guaranteed to be stable for this particular initial condition. This could be
circumvented by either performing more iterations or choosing an initial
condition closer to the origin such that it is located within
$\Gamma_{\bar{\alpha}}$. However, the scheme is still robust enough to
stabilize the system in this particular case. Moreover, the dependency on the
sampling time can also be seen in this particular trajectory as a higher
sampling time would take the next sampling point $\boldsymbol{x}_{c}(\Delta
t;\underaccent{\bar}{\vm x}_{0})$ farther away from the region of attraction
and thus possibly destabilize the scheme.

Figure 5: Predicted and closed-loop trajectories of the controlled system (60)
for different initial conditions.

##  6 Conclusion

This paper presented a DMPC scheme for continuous-time nonlinear systems which
relies on sensitivities to cooperatively solve the underlying OCP in each MPC
sampling step in a distributed fashion. The agents are dynamically coupled in
a neighbor-affine form which has the consequence that the sensitivities can be
evaluated locally. The algorithm only requires local computations and one
neighbor-to-neighbor communication step per iteration and thus constitutes a
scalable and fully distributed DMPC scheme. Furthermore, it was shown that
linear convergence can be guaranteed as long as the prediction horizon is
sufficiently small which represents a compromise between the domain of
attraction and convergence speed of the algorithm. The computation and
communication load is limited by either a contracting stopping criterion or a
fixed number of algorithm iterations. For both scenarios exponential stability
of the DMPC controller is shown. Numerical evaluations have demonstrated the
effectiveness of the presented scheme as only a few iterations per MPC step
are necessary to achieve a nearly optimal performance.

Future research concerns the experimental validation of the presented DMPC
scheme. In addition, convergence of the algorithm will be analyzed in the
presence of packet loss and in the asynchronous setting to further reduce the
execution time.

## References

  * 1 Rawlings J, Mayne D, Diehl M. Model predictive control: theory, computation, and design.  Nob Hill Publishing Madison.  2nd ed., 2017.
  * 2 Mayne D, Rawlings J, Rao C, Scokaert P. Constrained model predictive control: stability and optimality. Automatica. 2000;36(6):789–814.
  * 3 Grüne L, Pannek J. Nonlinear model predictive control: Theory and Algorithms.  Springer, 2011.
  * 4 Negenborn R, Maestre J. Distributed model predictive control: An overview and roadmap of future research opportunities. Control Systems Magazine. 2014;34(4):87–97.
  * 5 Camponogara E, Jia D, Krogh BH, Talukdar S. Distributed model predictive control. Control Systems Magazine. 2002;22(1):44–52.
  * 6 Christofides P, Scattolini R, La Muñoz de Peña D, Liu J. Distributed model predictive control: A tutorial review and future research directions. Computers & Chemical Engineering. 2013;51:21–41.
  * 7 Al-Gherwi W, Budman H, Elkamel A. Robust distributed model predictive control: A review and recent developments. Canadian Journal of Chemical Engineering. 2011;89(5):1176–1190.
  * 8 Scheu H, Marquardt W. Sensitivity-based coordination in distributed model predictive control. Journal of Process Control. 2011;21(5):715–728.
  * 9 Liu J, Peña M. d. lD, Christofides PD. Distributed model predictive control of nonlinear process systems. AIChE Journal. 2009;55(5):1171–1184.
  * 10 Venkat AN, Rawlings JB, Wright SJ. Distributed model predictive control of large-scale systems. In: Findeisen R, Allgöwer F, Biegler L. , eds. Workshop on Assessment and Future Directions of Nonlinear Model Predictive Control, 1st ed., Lecture Notes in Control and Information Sciences. Springer, 2007:591–605.
  * 11 Jiang Y, Sauerteig P, Houska B, Worthmann K. Distributed optimization using ALADIN for MPC in smart grids. IEEE Transactions on Control Systems Technology. 2020;29(5):2142–2152.
  * 12 Braun P, Grüne L, Kellett CM, Weller SR, Worthmann K. A distributed optimization algorithm for the predictive control of smart grids. IEEE Transactions on Automatic Control. 2016;61(12):3898–3911.
  * 13 Venkat AN, Hiskens IA, Rawlings JB, Wright SJ. Distributed MPC strategies with application to power system automatic generation control. IEEE Transactions on control systems technology. 2008;16(6):1192–1206.
  * 14 Leirens S, Zamora C, Negenborn R, De Schutter B. Coordination in urban water supply networks using distributed model predictive control. In: Proceedings of the American Control Conference (ACC). 2010:3957–3962.
  * 15 Trnka P, Pekař J, Havlena V. Application of distributed MPC to Barcelona water distribution network. In: Proceedings of the IFAC World Congress. 2011; Italy: Mailand:9139–9144.
  * 16 Hentzelt S, Klingler A, Graichen K. Experimental results for distributed model predictive control applied to a water distribution system. In: International Symposium on Intelligent Control (ISIC). 2014; France: Juan Les Pins:1100–1106.
  * 17 Wehbeh J, Rahman S, Sharf I. Distributed model predictive control for UAVs collaborative payload transport. In: Proceedings of the Conference on Intelligent Robots and Systems (IROS). 2020; USA: Las Vegas:11666–11672.
  * 18 Mulagaleti SK, Van Parys R, Pipeleers G. Distributed model predictive control of multiple vehicles transporting a flexible payload. In: Proceedings of the European Control Conference (ECC). 2018; Cyprus: Limassol:1417–1422.
  * 19 Köhler J, Müller MA, Allgöwer F. Distributed model predictive control—recursive feasibility under inexact dual optimization. Automatica. 2019;102:1–9.
  * 20 Conte C, Jones CN, Morari M, Zeilinger MN. Distributed synthesis and stability of cooperative distributed model predictive control for linear systems. Automatica. 2016;69:117–125.
  * 21 Farina M, Scattolini R. Distributed predictive control: A non-cooperative algorithm with neighbor-to-neighbor communication for linear systems. Automatica. 2012;48(6):1088-1096.
  * 22 Giselsson P, Doan M, Keviczky T, De Schutter B, Rantzer A. Accelerated gradient methods and dual decomposition in distributed model predictive control. Automatica. 2013;49(3):829–833.
  * 23 Rostami R, Costantini G, Görges D. ADMM-based distributed model predictive control: Primal and dual approaches. In: Proceedings of the Conference on Decision and Control (CDC). 2017; Australia: Melbourne:6598–6603.
  * 24 Bestler A, Graichen K. Distributed model predictive control for continuous-time nonlinear systems based on suboptimal ADMM. Optimal Control Applications and Methods. 2019;40(1):1–23.
  * 25 Burk D, Völz A, Graichen K. Towards a modular framework for distributed model predictive control of nonlinear neighbor-affine systems. In: Proceedings of the Conference on Decision and Control (CDC). 2019; France: Nice:5279–5284.
  * 26 Burk D, Völz A, Graichen K. A modular framework for distributed model predictive control of nonlinear continuous-time systems (GRAMPC-D). Optimization and Engineering. 2021;23(2):771–795.
  * 27 Mirzaei A, Ramezani A. Cooperative optimization-based distributed model predictive control for constrained nonlinear large-scale systems with stability and feasibility guarantees. ISA Transactions. 2021;116:81–96.
  * 28 Dunbar WB. Distributed receding horizon control of dynamically coupled nonlinear systems. IEEE Transactions on Automatic Control. 2007;52(7):1249–1263.
  * 29 Zhihao C, Longhong W, Jiang Z, Kun W, Yingxun W. Virtual target guidance-based distributed model predictive control for formation control of multiple UAVs. Chinese Journal of Aeronautics. 2020;33(3):1037–1056.
  * 30 Zhao R, Miao M, Lu J, Wang Y, Li D. Formation control of multiple underwater robots based on ADMM distributed model predictive control. Ocean Engineering. 2022;257.
  * 31 Liu P, Kurt A, Ozguner U. Distributed model predictive control for cooperative and flexible vehicle platooning. IEEE Transactions on Control Systems Technology. 2018;27(3):1115–1128.
  * 32 Müller MA, Reble M, Allgöwer F. A general distributed MPC framework for cooperative control. In: Proceedengs of the IFAC World Congress. 2011; Italy: Milan:7987–7992.
  * 33 Lyu Y, Hu J, Chen BM, Zhao C, Pan Q. Multivehicle flocking with collision avoidance via distributed model predictive control. IEEE Transactions on Cybernetics. 2019;51(5):2651–2662.
  * 34 Farina M, Perizzato A, Scattolini R. Application of distributed predictive control to motion and coordination problems for unicycle autonomous robots. Robotics and Autonomous Systems. 2015;72:248–260.
  * 35 Huber H, Graichen K. A sensitivity-based distributed model predictive control algorithm for nonlinear continuous-time systems. In: Proceedings of the Conference on Control Technology and Applications (CCTA). 2021; USA: San Diego:195–201.
  * 36 Costantini G, Rostami R, Görges D. Decomposition methods for distributed quadratic programming with application to distributed model predictive control. In: Proceedings of the Conference on Communication, Control, and Computing. 2018; USA: Allerton:943–950.
  * 37 Pierer von Esch M, Völz A, Graichen K. Sensitivity-based distributed model predictive control: synchronous and asynchronous execution compared to ADMM. at-Automatisierungstechnik. 2024;72(2):91–106.
  * 38 Pierer von Esch M, Völz A, Graichen K. Asynchronous ADMM for nonlinear continuous-time systems. Optimal Control Applications and Methods. 2023\.
  * 39 Belgioioso G, Liao-McPherson D, Badyn dMH, Pelzmann N, Lygeros J, Dörfler F. Stability and robustness of distributed suboptimal model predictive control. In: Proceedings of the IFAC World Congress. 2023; Japan: Yokohama:5115–5120.
  * 40 Engelmann A, Stomberg G, Faulwasser T. An essentially decentralized interior point method for control. In: Proceedings of the Conference on Decision and Control (CDC). 2021; USA: Austin:2414-2420.
  * 41 Burk D, Völz A, Graichen K. Towards asynchronous ADMM for distributed model predictive control of nonlinear systems. In: Proceedings of the European Control Conference (ECC). 2021; Netherlands: Rotterdam:1957–1962.
  * 42 Rostami R, Costantini G, Görges D. Communication reduction in stabilizing distributed model predictive control. In: Proceedings of the Conference on Decision and Control (CDC). 2018; USA: Miami:5164–5169.
  * 43 Müller MA, Allgöwer F. Economic and distributed model predictive control: Recent developments in optimization-based control. SICE Journal of Control, Measurement, and System Integration. 2017;10(2):39–52.
  * 44 Stomberg G, Engelmann A, Faulwasser T. A compendium of optimization algorithms for distributed linear-quadratic MPC. at-Automatisierungstechnik. 2022;70(4):317–330.
  * 45 Braun P, Grüne L. Verteilte Optimierung: Anwendungen in der Modellprädiktiven Regelung (German). at-Automatisierungstechnik. 2018;66(11):939–949.
  * 46 Boyd S, Parikh N, Chu E, Peleato B, Eckstein J. Distributed optimization and statistical learning via the alternating direction method of multipliers. Foundations and Trends in Machine Learning. 2010;3(1):1–122.
  * 47 Houska B, Frasch J, Diehl M. An augmented Lagrangian based algorithm for distributed nonconvex optimization. SIAM Journal on Optimization. 2016;26(2):1101–1127.
  * 48 Engelmann A, Jiang Y, Houska B, Faulwasser T. Decomposition of nonconvex optimization via bi-level distributed ALADIN. IEEE Transactions on Control of Network Systems. 2020;7(4):1848–1858.
  * 49 Farokhi F, Shames I, Johansson KH. Distributed MPC via dual decomposition and alternative direction method of multipliers. In: Maestre JM, Negenbor RR. , eds. Distributed model predictive control made easy, 1st ed., Intelligent Systems, Control and Automation: Science and Egineering. Springer, 2013:115–131.
  * 50 Kögel M, Findeisen R. Cooperative distributed MPC using the alternating direction multiplier method. In: Proceedings of the Symposium on Advanced Control of Chemical Processes (ADCHEM). 2012; Singapur: Singapur:445–450.
  * 51 Venkat AN, Rawlings JB, Wright SJ. Stability and optimality of distributed model predictive control. In: Proceedings of the Conference on Decision and Control (CDC). 2005; Spain: Seville:6680–6685.
  * 52 Doan M, Diehl M, Keviczky T, De Schutter B. A Jacobi decomposition algorithm for distributed convex optimization in distributed model predictive control. In: Proceedings of the IFAC World Congress. 2017; France: Toulouse:4905–4911.
  * 53 Stewart BT, Venkat AN, Rawlings JB, Wright SJ, Pannocchia G. Cooperative distributed model predictive control. Systems & Control Letters. 2010;59(8):460–469.
  * 54 Stewart BT, Wright SJ, Rawlings JB. Cooperative distributed model predictive control for nonlinear systems. Journal of Process Control. 2011;21(5):698–704.
  * 55 Paulen R, Nazari S, Shahidi SA, Sonntag C, Engell S. Primal and dual decomposition for distributed MPC—Theory, implementation, and comparison in a SoS simulation framework. In: Proceedings of the Mediterranean Conference on Control and Automation (MED). 2016; Greece: Athens:286–291.
  * 56 Pflaum P, Alamir M, Lamoudi MY. Comparison of a primal and a dual decomposition for distributed MPC in smart districts. In: Proceedings of the Conference on Smart Grid Communications (SmartGridComm). 2014:55–60.
  * 57 Scheu H, Busch J, Marquardt W. Nonlinear distributed dynamic optimization based on first order sensitivities. In: Proceedings of the American Control Conference (ACC). 2010; USA: Baltimore:1574–1579.
  * 58 Alvarado I, Limon D, De La Peña DM, et al. A comparative analysis of distributed MPC techniques applied to the HD-MPC four-tank benchmark. Journal of Process Control. 2011;21(5):800–815.
  * 59 Huber H, Burk D, Graichen K. Comparison of Sensitivity-Based and ADMM-Based DMPC Applied to Building Automation. In: Proceedings of the Conference on Control Technology and Applications (CCTA). 2022; Italy: Trieste:546-553.
  * 60 Cohen G. On an algorithm of decentralized optimal control. Journal of Mathematical Analysis and Applications. 1977;59(2):242–259.
  * 61 Mesarovic MD, Macko D, Takahara Y. Theory of hierarchical, multilevel, systems.  Elsevier, 2000.
  * 62 Darure T, Puig V, Yamé J, Hamelin F, Wang Y. Distributed Model Predictive Control applied to a VAV based HVAC system based on Sensitivity Analysis. In: Proceedings of the Conference on Nonlinear Model Predictive Control. 2018; USA: Wisconsin:259–264.
  * 63 Pierer von Esch M, Völz A, Graichen K. A fixed-point iteration scheme for sensitivity-based distributed optimal control. IEEE Transactions on Automatic Control (submitted). 2023\.
  * 64 Müller MA, Reble M, Allgöwer F. Cooperative control of dynamically decoupled systems via distributed model predictive control. International Journal of Robust and Nonlinear Control. 2012;22(12):1376–1397.
  * 65 Müller MA, Schürmann B, Allgöwer F. Robust cooperative control of dynamically decoupled systems via distributed MPC. Proceedings of the Conference on Nonlinear Model Predictive Control. 2012;45(17):412–417.
  * 66 Costantini G, Rostami R, Görges D. Distributed linear quadratic regulator for the synthesis of a separable terminal cost for distributed model predictive control. In: Proceedings of the Conference on Decision and Control (CDC). 2018; USA: Miami:5170–5175.
  * 67 Costantini G, Rostami R, Görges D. Fast separable terminal cost synthesis for distributed MPC. In: Proceedings of the European Control Conference (ECC). 2019; Italy: Naples:313–318.
  * 68 Graichen K, Kugi A. Stability and incremental improvement of suboptimal MPC without terminal constraints. IEEE Transactions on Automatic Control. 2010;55(11):2576–2580.
  * 69 Grune L, Rantzer A. On the Infinite Horizon Performance of Receding Horizon Controllers. IEEE Transactions on Automatic Control. 2008;53(9):2100-2111.
  * 70 Giselsson P, Rantzer A. On feasibility, stability and performance in distributed model predictive control. IEEE Transactions on Automatic Control. 2013;59(4):1031–1036.
  * 71 Jané-Soneira P, Gießler A, Pfeifer M, Hohmann S. Distributed Suboptimal Model Predictive Control with Minimal Information Exchange. In: Proceedings of the Conference on Decision and Control (CDC). 2022; Mexico: Cancun:3939–3946.
  * 72 Graichen K. A fixed-point iteration scheme for real-time model predictive control. Automatica. 2012;48(7):1300–1305.
  * 73 Graichen K, Käpernick B. A Real-Time Gradient Method for Nonlinear Model Predictive Control. In: Zheng T. , ed. Frontiers of Model Predictive Control, 1st ed., IntechOpen. , 2012.
  * 74 Jadbabaie A, Yu J, Hauser J. Unconstrained receding-horizon control of nonlinear systems. IEEE Transactions on Automatic Control. 2001;46(5):776–783.
  * 75 Chen H, Allgöwer F. A quasi-infinite horizon nonlinear model predictive control scheme with guaranteed stability. Automatica. 1998;34(10):1205–1217.
  * 76 Parisini T, Zoppoli R. A receding-horizon regulator for nonlinear systems and a neural approximation. Automatica. 1995;31(10):1443–1451.
  * 77 Michalska H, Mayne DQ. Robust receding horizon control of constrained nonlinear systems. IEEE Transactions on Automatic Control. 1993;38(11):1623–1633.
  * 78 Limón D, Alamo T, Salas F, Camacho EF. On the stability of constrained MPC without terminal constraint. IEEE Transactions on Automatic Control. 2006;51(5):832–836.
  * 79 Zhang RC, Li W, Li N. On the optimal control of network LQR with spatially-exponential decaying structure. In: Proceedings of the American Control Conference (ACC). 2023; USA: San Diego:1775–1780.
  * 80 Siljak DD. Decentralized control of complex systems.  Courier Corporation, 2011.
  * 81 Xu L, Guo B, Galimberti C, Farina M, Carli R, Trecate GF. Suboptimal distributed LQR design for physically coupled systems. Proceedings of the IFAC World Congress. 2020;53(2):11032–11037.
  * 82 Deroo F, Ulbrich M, Anderson BD, Hirche S. Accelerated iterative distributed controller synthesis with a Barzilai-Borwein step size. In: Proceedings of the Conference on Decision and Control (CDC). 2012; USA: Maui:4864–4870.
  * 83 Mårtensson K, Rantzer A. A scalable method for continuous-time distributed control synthesis. In: Proceedings of the American Control Conference (ACC). 2012; Canada: Montreal:6308–6313.
  * 84 Boyd S, El Ghaoui L, Feron E, Balakrishnan V. Linear matrix inequalities in system and control theory.  SIAM, 1994.
  * 85 Englert T, Völz A, Mesmer F, Rhein S, Graichen K. A software framework for embedded nonlinear model predictive control using a gradient-based augmented Lagrangian approach (GRAMPC). Optimization and Engineering. 2019;20(3):769–809.
  * 86 Yuanqiang Z, Dewei L, Yugeng X, Lihui C, Yuli X, Zhongxue G. Distributed model predictive control for nonlinear networked systems with asynchronous communication. In: Proceedings of the Chinese Control Conference (CCC). 2017; China: Dalian.
  * 87 Zhou Y, Li D, Lu J, Xi Y, Cen L. Networked and distributed predictive control of non-linear systems subject to asynchronous communication. IET Control Theory & Applications. 2018;12(4):504–514.
  * 88 Zhou Y, Li D, Xi Y, Gan Z. Periodic event-triggered control for distributed networked multiagents with asynchronous communication: A predictive control approach. International Journal of Robust and Nonlinear Control. 2019;29(1):43–66.

Computation of Sensitivities In order to derive the first-order sensitivities
for OCP (19), the OCP of an individual neighbor
$j\in\mathcal{N}_{i}^{\rightarrow}$ of the central OCP (3) from the
perspective of the agent $i\in\mathcal{V}$ is considered

$\displaystyle\min_{\boldsymbol{u}_{j}}\quad$ $\displaystyle
J_{j}(\boldsymbol{x}_{j},\boldsymbol{u}_{j},\boldsymbol{x}_{\mathcal{N}_{j}^{\leftarrow}};\boldsymbol{\underaccent{\bar}{x}}_{k}):=V_{j}(\boldsymbol{x}_{j}(T))+\int_{0}^{T}l_{j}(\boldsymbol{x}_{j},\boldsymbol{u}_{j},\boldsymbol{x}_{\mathcal{N}_{j}^{\leftarrow}})\,\mathrm{d}\tau$
(61a) $\displaystyle\leavevmode\nobreak\ \operatorname{s.t.}\quad$
$\displaystyle\boldsymbol{\dot{x}}_{j}=\boldsymbol{f}_{j}(\boldsymbol{x}_{j},\boldsymbol{u}_{j},\boldsymbol{x}_{\mathcal{N}_{j}^{\leftarrow}})\,,\quad\boldsymbol{x}_{j}(0)=\boldsymbol{\underaccent{\bar}{x}}_{j,k}$
(61b)
$\displaystyle\boldsymbol{u}_{j}(\tau)\in\mathbb{U}_{j}\,,\quad\tau\in[0,\,T]\,.$
(61c)

Now the first-order sensitivity is formulated as the Gâteaux derivative of OCP
(61)

$\delta
J_{j}(\boldsymbol{x}_{j}^{q-1},\boldsymbol{u}_{j}^{q-1},\boldsymbol{x}_{\mathcal{N}_{j}^{\leftarrow}}^{q-1};\boldsymbol{\underaccent{\bar}{x}}_{k})(\delta\boldsymbol{x}_{i})=\frac{\mathrm{d}J_{j}(\boldsymbol{x}_{j}^{q-1},\boldsymbol{u}_{j}^{q-1},\boldsymbol{x}_{\mathcal{N}_{j}^{\leftarrow}}^{q-1}+\epsilon\delta\boldsymbol{x}_{ji};\boldsymbol{\underaccent{\bar}{x}}_{k})}{\mathrm{d}\epsilon}\bigg{|}_{\epsilon=0}$
(62)

at some point
$(\boldsymbol{x}_{j}^{q-1},\boldsymbol{u}_{j}^{q-1},\boldsymbol{x}_{\mathcal{N}_{j}^{\leftarrow}}^{q-1})$
w.r.t. the admissible direction
$\delta\boldsymbol{x}_{ji}=[\boldsymbol{0}^{\mathrm{T}}\dots\boldsymbol{0}^{\mathrm{T}}\,\delta\boldsymbol{x}_{i}^{\mathrm{T}}\,\boldsymbol{0}^{\mathrm{T}}\dots\boldsymbol{0}^{\mathrm{T}}]^{\mathrm{T}}\in\mathds{R}^{p_{j}}$
where $\delta\boldsymbol{x}_{i}=\boldsymbol{x}_{i}-\boldsymbol{x}_{i}^{q-1}$
shows up at the position in $\delta\boldsymbol{x}_{ji}$ corresponding to
$\boldsymbol{x}_{i}$ in $\boldsymbol{x}_{\mathcal{N}_{j}^{\leftarrow}}$.
Adjoining the dynamics to the cost via (time-dependent) Lagrange multipliers
$\boldsymbol{\lambda}_{j}\in\mathds{R}^{n_{j}}$ results in

$\displaystyle\delta
J_{j}(\boldsymbol{x}_{j}^{q-1},\boldsymbol{u}_{j}^{q-1},\boldsymbol{x}_{\mathcal{N}_{j}^{\leftarrow}}^{q-1};\boldsymbol{\underaccent{\bar}{x}}_{k})(\delta\boldsymbol{x}_{i})$
$\displaystyle=\int_{0}^{T}\frac{\mathrm{d}}{\mathrm{d}\epsilon}l_{j}(\boldsymbol{x}_{j}^{q-1},\boldsymbol{u}_{j}^{q-1},\boldsymbol{x}_{\mathcal{N}_{j}^{\leftarrow}}^{q-1}+\epsilon\delta\boldsymbol{x}_{ji})+\frac{\mathrm{d}}{\mathrm{d}\epsilon}\big{(}(\boldsymbol{\lambda}_{j}^{q-1})^{\mathrm{T}}(\,\boldsymbol{f}_{j}(\boldsymbol{x}_{j}^{q-1},\boldsymbol{u}_{j}^{q-1},\boldsymbol{x}_{\mathcal{N}_{j}^{\leftarrow}}^{q-1}+\epsilon\delta\boldsymbol{x}_{ji})-\boldsymbol{\dot{x}}_{j}^{q-1})\big{)}\,\mathrm{d}t\,\bigg{|}_{\epsilon=0}$
$\displaystyle=\int_{0}^{T}\partial_{\boldsymbol{x}_{ji}}\big{(}l_{j}(\boldsymbol{x}_{j}^{q-1},\boldsymbol{u}_{j}^{q-1},\boldsymbol{x}_{\mathcal{N}_{j}^{\leftarrow}}^{q-1})+\boldsymbol{f}_{j}(\boldsymbol{x}_{j}^{q-1},\boldsymbol{u}_{j}^{q-1},\boldsymbol{x}_{\mathcal{N}_{j}^{\leftarrow}}^{q-1})^{\mathrm{T}}\boldsymbol{\lambda}_{j}^{q-1}\big{)}^{\mathrm{T}}\,\delta\boldsymbol{x}_{ji}\,\mathrm{d}\tau.$
(63)

Considering the particular formulation of dynamics (1) and cost terms (5) in
neighbor-affine form, the Gâteaux derivative further simplifies to

$\displaystyle\delta
J_{j}(\boldsymbol{x}_{j}^{q-1},\boldsymbol{u}_{j}^{q-1},\boldsymbol{x}_{\mathcal{N}_{j}^{\leftarrow}}^{q-1};\boldsymbol{\underaccent{\bar}{x}}_{k})(\delta\boldsymbol{x}_{i})=\int_{0}^{T}\big{(}\partial_{\boldsymbol{x}_{i}}l_{ji}(\boldsymbol{x}_{j}^{q-1},\boldsymbol{x}_{i}^{q-1})+\partial_{\boldsymbol{x}_{i}}\,\boldsymbol{f}_{ji}(\boldsymbol{x}_{j}^{q-1},\boldsymbol{x}_{i}^{q-1})^{\mathrm{T}}\boldsymbol{\lambda}_{j}^{q-1}\big{)}^{\mathrm{T}}\,\delta\boldsymbol{x}_{i}\,\mathrm{d}\tau.$
(64)

In OCP (19), the sensitivities are calculated recursively for the already
augmented cost functional. Recalling the definition of
$\delta\boldsymbol{x}_{i}=\boldsymbol{x}_{i}-\boldsymbol{x}_{i}^{q-1}$, it is
clear that (64) is linear w.r.t. $\boldsymbol{x}_{i}$ and can be incorporated
as a linear term into the local cost function term
$l_{ii}(\boldsymbol{x}_{i},\boldsymbol{u}_{i})$ in (5). Thus, a repeated
application of the Gâteaux derivative (62) still results in (64), i.e. the
sensitivity (20)
$\delta\bar{J}_{j}(\boldsymbol{u}_{j}^{q-1};\boldsymbol{\underaccent{\bar}{x}}_{k})$
of the augmented cost function in (19a). Proof of Lemma 1

######  Proof .1.

The following considerations require several Lipschitz estimates. Based on the
continuous differentiability of the dynamics (29), there exists a finite
(local) Lipschitz constant $L_{f}<\infty$ for some $r_{x}>0$ such that

$\displaystyle\|\boldsymbol{\hat{f}}(\boldsymbol{x},\boldsymbol{u},\boldsymbol{\hat{x}})-\boldsymbol{\hat{f}}(\boldsymbol{y},\boldsymbol{v},\boldsymbol{\hat{y}})\|$
$\displaystyle\leq
L_{f}(\|\boldsymbol{x}-\boldsymbol{y}\|+\|\boldsymbol{u}-\boldsymbol{v}\|+\|\boldsymbol{\hat{x}}-\boldsymbol{\hat{y}}\|)$
(65)

for all
$\boldsymbol{x},\boldsymbol{y},\boldsymbol{\hat{x}},\boldsymbol{\hat{y}}\in\Gamma_{\alpha}^{r_{x}}$
and $\boldsymbol{u},\boldsymbol{v}\in\mathbb{U}$, where
$\Gamma_{\alpha}^{r_{x}}:=\bigcup_{\boldsymbol{x}_{k}\in\Gamma_{\alpha}}\mathcal{B}(\boldsymbol{x}_{k},r_{x})$
is the $r_{x}$-neighborhood to the (compact) domain of attraction
$\Gamma_{\alpha}$. Similarly, consider the stacked adjoint dynamics (22) in
each iteration $q$ of Algorithm 1, i.e.

$\boldsymbol{\dot{\lambda}}^{q}(\tau)=-\partial_{\boldsymbol{x}}\sum_{i\in\mathcal{V}}H_{i}(\boldsymbol{x}_{i},\boldsymbol{u}_{i},\boldsymbol{\lambda}_{i})=:\boldsymbol{G}_{d}(\boldsymbol{x}^{q},\boldsymbol{u}^{q},\boldsymbol{\lambda}^{q};\boldsymbol{\bar{x}}^{q-1},\boldsymbol{\bar{\lambda}}^{q-1})\,,\quad\boldsymbol{\lambda}^{q}(T)=\partial_{\boldsymbol{x}}V(\boldsymbol{x}^{q}(T))\,,$
(66)

where the notation
$\boldsymbol{\bar{x}}=[\boldsymbol{x}_{i}^{\mathrm{T}}\,\boldsymbol{x}_{\mathcal{N}_{i}}^{\mathrm{T}}]_{i\in\mathcal{V}}\in\mathds{R}^{p_{x}}$,
$p_{x}=\sum_{i\in\mathcal{V}}(n_{i}+\sum_{j\in\mathcal{N}_{i}}n_{j})$ and
$\boldsymbol{\bar{\lambda}}=[\boldsymbol{\lambda}_{\mathcal{N}_{i}^{\rightarrow}}]_{i\in\mathcal{V}}\in\mathds{R}^{p_{\lambda}}$,
$p_{\lambda}=\sum_{i\in\mathcal{V}}\sum_{j\in\mathcal{N}_{i}^{\rightarrow}}n_{j}$
explicitly captures the dependency on the trajectories of iteration $q-1$. Due
to the assumed differentiability of all dynamics and cost functions, there
exists a finite Lipschitz constant $L_{G}<\infty$ for some $r_{\lambda}>0$
such that

$\displaystyle\|\boldsymbol{G}_{d}(\boldsymbol{x},\boldsymbol{u},\boldsymbol{\lambda};\boldsymbol{\bar{x}},\boldsymbol{\bar{\lambda}})-\boldsymbol{G}_{d}(\boldsymbol{y},\boldsymbol{v},\boldsymbol{\mu};\boldsymbol{\bar{y}},\boldsymbol{\bar{\mu}})\|$
$\displaystyle\leq
L_{G}(\|\boldsymbol{x}-\boldsymbol{y}\|+\|\boldsymbol{u}-\boldsymbol{v}\|+\|\boldsymbol{\lambda}-\boldsymbol{\mu}\|+\|\boldsymbol{\bar{x}}-\boldsymbol{\bar{y}}\|+\|\boldsymbol{\bar{\lambda}}-\boldsymbol{\bar{\mu}}\|)$
(67)

for all for all
$\boldsymbol{x},\boldsymbol{y},\boldsymbol{\bar{x}},\boldsymbol{\bar{y}}\in\Gamma_{\alpha}^{r_{x}}$
and $\boldsymbol{u},\boldsymbol{v}\in\mathbb{U}$ as well as
$\boldsymbol{\lambda},\boldsymbol{\mu},\boldsymbol{\bar{\lambda}},\boldsymbol{\bar{\mu}}\in\mathcal{S}^{r_{\lambda}}$,
where and $\mathcal{S}^{r_{\lambda}}$ is the $r_{\lambda}$-neighborhood to the
compact set
$\mathcal{S}=\\{\partial_{\boldsymbol{x}}V(\boldsymbol{x})|\boldsymbol{x}\in\Gamma_{\alpha}^{r_{x}}\\}$.

At first it is shown that the iterates are bounded for each
$\boldsymbol{\underaccent{\bar}{x}}_{k}\in\Gamma_{\alpha}$ for a sufficiently
short prediction horizon $T$, i.e.
$\boldsymbol{x}^{q}(\tau)\in\Gamma_{\alpha}^{r_{x}}$ and
$\boldsymbol{\lambda}^{q}(\tau)\in\mathcal{S}^{r_{\lambda}}$, $\tau\in[0,\,T]$
for $q=1,2,\dots,q_{\mathrm{max}}$. We proceed by induction. To this end,
assume that $\boldsymbol{x}^{q-1}(\tau)\in\Gamma_{\alpha}^{r_{x}}$ and
$\boldsymbol{\lambda}^{q-1}(\tau)\in\mathcal{S}^{r_{\lambda}}$,
$\tau\in[0,\,T]$ and consider the integral form of the dynamics (29) for
$\boldsymbol{\underaccent{\bar}{x}}_{k}\in\Gamma_{\alpha}$. By adding and
subtracting
$\boldsymbol{\hat{f}}(\boldsymbol{\underaccent{\bar}{x}}_{k},\boldsymbol{0},\boldsymbol{\underaccent{\bar}{\hat
x}}_{k})$, with $\boldsymbol{\underaccent{\bar}{\hat
x}}_{k}=[[\boldsymbol{\underaccent{\bar}{x}}_{j,k}]_{j\in\mathcal{N}_{i}^{\leftarrow}}]_{i\in\mathcal{V}}$
as well as using the Lipschitz property (65), we get (omitting time arguments)

$\displaystyle\|\boldsymbol{x}^{q}(\tau;\boldsymbol{\underaccent{\bar}{x}}_{k})-\boldsymbol{\underaccent{\bar}{x}}_{k}\|$