# Compositional Synthesis for Linear Systems via Convex Optimization of
Assume-Guarantee Contracts

Kasra Ghasemi<EMAIL_ADDRESS>Sadra Sadraddini<EMAIL_ADDRESS>Calin Belta
<EMAIL_ADDRESS>Boston University, Boston, MA  Massachusetts Institute of
Technology, Cambridge, MA

###### Abstract

We take a divide and conquer approach to design controllers for reachability
problems given large-scale linear systems with polyhedral constraints on
states, controls, and disturbances. Such systems are made of small subsystems
with coupled dynamics. We treat the couplings as additional disturbances and
use assume-guarantee (AG) contracts to characterize these disturbance sets.
For each subsystem, we design and implement a robust controller locally,
subject to its own constraints and contracts. The main contribution of this
paper is a method to derive the contracts via a novel parameterization and a
corresponding potential function that characterizes the distance to the
correct composition of controllers and contracts, where all contracts are
held. We show that the potential function is convex in the contract
parameters. This enables the subsystems to negotiate the contracts with the
gradient information from the dual of their local synthesis optimization
problems in a distributed way, facilitating compositional control synthesis
that scales to large systems. We present numerical examples, including a
scalability study on a system with tens of thousands of dimensions, and a case
study on applying our method to a distributed Model Predictive Control (MPC)
problem in a power system.

###### keywords:

Assume-guarantee contract; compositional correct-by-construction synthesis;
distributed robust model predictive control; viable set; invariant set;
zonotope; interconnected linear systems.

, ,

##  1 Introduction

Verification and synthesis of controllers for large-scale systems with hard
constraints are encountered in many applications such as resource management,
biology, and robotics [3, 2]. The large-scale nature of such systems poses
certain computational challenges for analysis and controller design. First,
verification and synthesis problems usually reduce to optimization problems
with prohibitive sizes. This particularly concerns the applicability of real-
time algorithms such as Model Predictive Control (MPC). Second, large-scale
systems are typically networks of many compartments with physically separated
actuators. A central command unit would lead to a complex communication
architecture during implementation, commonly referred to as online mode. Thus
decentralized computations of the controllers are preferred such that each
subsystem computes its own decisions based on local information.

Roughly, decentralized control can be categorized into two approaches. In the
first, the synthesis problem is solved in a centralized fashion, but the
controllers are deployed in a decentralized manner in the online mode. In this
case, one has to pose structural constraints and costs on the sparsity of the
control law and solve a large optimization problem subject to those
constraints (in offline mode) before the deployment phase. The works in [34,
14, 26, 29] follow this approach. In the second, the synthesis problem itself
is solved in a distributed manner. In other words, the control laws are
determined by iteratively solving local synthesis problems. In this case,
there is no single large optimization problem, and the corresponding
computation time is significantly reduced.

In this paper, we are interested in the latter case, which is more
challenging, as one has to formulate and solve a distributed optimization
problem. We consider linear discrete-time systems with additive set-valued
bounded disturbances and constrained states and control inputs. Couplings
across subsystems are allowed over states, controls, and constraints. The
problem is to design decentralized feedback controllers such that all closed-
loop trajectories of the aggregate system satisfy all the constraints under
all disturbances, while a cost function is optimized. The cost can be chosen
trivially as a constant if the goal is to only find feasible controllers. If
this is not the case, we assume that the cost function is separable over the
subsystems. Both finite-time (with possibly time-varying system dynamics and
constraints) and infinite-time (in the form of robust set-invariance) problems
are considered. The subsystems can communicate offline to synthesize the
decentralized controllers and identify their closed-loop set-valued
trajectories, also known as reachability sets. Thus, online computations use
only local information. If online communication is possible, controllers and
set-valued trajectories can be recomputed online in an MPC fashion with
recursive feasibility guarantee.

The system interacts with its external world, which we refer to as its
environment. Exogenous disturbances are examples of such interactions. Central
to our approach is a divide and conquer approach based on assume-guarantee
(AG) reasoning [20]. We use the notion of AG contracts to formalize the
promises a system takes from and makes to its environment. One use of AG
contracts is to decouple subsystems within a network by constructing one
contract for each subsystem. By doing so, we are seeing the effect of
couplings as exogenous, which can lead to mismatch between the assumptions
over the environment and what happens in reality to some subsystems. Correct
composition happens when there is no mismatch in a set of AG contracts.
Achieving correct composition is a challenging problem. While AG contracts
have been used in control synthesis [19, 30, 6], the focus has been on how to
use AG contracts given a priori. In contrast, the focus of this paper is on
synthesizing contracts. The main contributions of this paper are as follows:

  * •

We provide a novel parameterization of AG contracts in the form of zonotopes
and a corresponding potential function that quantifies how far a set of
contracts is from correct composition. We show that the potential function and
the set of parameters leading to correct composition are convex (Section 4 and
Section 5).

  * •

We develop an algorithm for subsystems to coordinate AG contracts and solve
local optimization problems such that all subsystems’ closed-loop responses
are guaranteed to satisfy the constraints over states and control inputs. We
also show that our distributed parameter synthesis method is not more
conservative than the centralized version of contract optimization (Section
6).

  * •

We include examples showing that our method can handle in a reasonable amount
of time systems of orders of magnitude larger than what centralized methods
can handle. We also provide an application to MPC for a power network (Section
7 and Section 8).

###  1.1 Related Work

Distributed MPC  Distributed optimization for MPC has a long history. Earlier
works such as [40, 42] focus on systems with decoupled dynamics. This way, the
optimization problem needs to be only decomposed through the cost function.
The work in [12] handles coupled dynamics and formulates MPC as a distributed
optimization problem that can be solved using methods such as Alternating
Direction Method of Multipliers (ADMM) [7]. However, the authors of [12]
consider deterministic systems only and there are no set-valued couplings. In
this paper, we deal with set-valued uncertainties and couplings through
assume-guarantee contracts.

Contracts for Control  As mentioned earlier, the majority of works assume
contracts given a priori and focus on how to use them for control synthesis
[30, 37]. In [38], the whole feasible sets of the individual subsystems were
considered in the contracts. This approach is very conservative and the
problem becomes infeasible in case the feasible sets are unbounded as it makes
the couplings also unbounded. The work in [22] studied parametric assume-
guarantee contracts and developed a notion of small gain theorem for systems
with dynamical couplings. In [35, 13, 27] the contracts were synthesized in a
centralized optimization problem, but controllers were deployed in a
decentralized way. The work in [23] found contracts through binary search over
box-shaped contract sets for monotone systems and specifications. The
parameterization of contracts in [27] leads to a non-convex optimization
problem that is conservatively approximated by a convex semi definite program.
In contrast to the works mentioned above, contract synthesis as well as the
control synthesis are carried out in a distributed way in our paper. In
addition, our parameterization has convexity properties, which combined with
the distributability mentioned above, lead to scalability and convergence of
our proposed method over a finite number of iterations.

Controller Parameterization Synthesizing controllers using convex programs is
a well studied subject [7]. A seminal result is that linear feedback over past
disturbances instead of the current state leads to a convex parameterization
of system’s predicted behaviours, which facilitates optimization subject to
convex constraints [16]. The same idea can be applied to an output feedback
setting [1]. The authors in [32] also use a similar technique but for infinite
time set invariance. The zonotope-based controllers in this paper are loosely
based on disturbance feedback policies. By leveraging the properties of
zonotopes and of zonotope operations, such as zonotope order reduction, we are
able to provide a distributed set-based computation of robust controllers
subject to bounded disturbance sets, which is correct-by-construction and
scalable.

Conference Version  Some of the results from this paper appeared in the
conference version [15], where we introduced the convex parameterization and
the synthesis method. The parameterization defined here is richer than the one
in [15] due to the inclusion of contract set centroids, which is particularly
useful for reachability problems. This paper also considers optimality in
control synthesis as part of MPC, while [15] only focused on reach and avoid
constraints. In this paper, we also treat coupled constraints and provide a
semidefinite programming approach to decouple them with the least amount of
conservativesness. Finally, complete proofs, updated examples, and an MPC case
study are included in addition to the material from [15].

##  2 Preliminaries and Definitions

###  2.1 Notation

The sets of real, non-negative real, integers, and non-negative integers are
represented by $\mathbb{R}$, $\mathbb{R}_{+}$, $\mathbb{N}$, and
$\mathbb{N}_{+}$, respectively. We use
$\mathbb{N}_{h_{1},h_{2}}:=\\{h_{1},h_{1}+1,\cdots,h_{2}\\}$,
$h_{1},h_{2}\in\mathbb{N},h_{1}<h_{2}$, and
$\mathbb{N}_{h}:=\mathbb{N}_{0:h}$. Additionally, $I_{n}$, $0_{n}$, and
$\mathbb{1}_{n}$ represent the $n$-dimensional identity matrix, vector of
zeros, and vector of ones, respectively. Given matrices
$A_{1},A_{2},\cdots,A_{N}$ with the same number of rows,
$[A_{1},A_{2},\cdots,A_{N}]$ is their horizontal concatenation. Given a matrix
$A$, $|A|$ is its element-wise absolute value. We use $\text{Blk}(.)$ to
denote block diagonal concatenation of a matrix and sum() to denote the sum of
all its elements. Given $\alpha\in\mathbb{R}^{n}$,
$\operatorname{{\text{Diag}}}(\alpha)$ is the diagonal matrix composed of the
entries of $\alpha$. The transpose of matrix $A$ is denoted by $A^{T}$.
Inequality relations denoted by $\geq$ and $\leq$ over matrices are
interpreted element-wise, and $A\succeq 0$ indicates that $A$ is a positive
semi-definite matrix. The infinity-norm of matrix $A\in\mathbb{R}^{n_{1}\times
n_{2}},n_{1},n_{2}\in\mathbb{N}_{+}$ is defined as
$\|A\|_{\infty}:=\text{max}(|A|\mathbb{1}_{n_{1}})$, where the max function
returns the largest element.

Given $\mathbb{S}_{1},\mathbb{S}_{2}\subseteq\mathbb{R}^{n}$, their _Minkowski
sum_ is denoted by
$\mathbb{S}_{1}\oplus\mathbb{S}_{2}:=\\{s_{1}+s_{2}|s_{1}\in\mathbb{S}_{1},s_{2}\in\mathbb{S}_{2}\\}$.
We interpret $s+\mathbb{S}$ as $\\{s\\}\oplus\mathbb{S}$, and
$\mu(\mathbb{S})$ as $\\{\mu(s)|s\in\mathbb{S}\\}$, where $\mu$ is a function
with $\mathbb{S}$ within its domain. The Directed Hausdorff distance
$d_{DH}(\mathbb{S}_{1},\mathbb{S}_{2})$ quantifies how distant is
$\mathbb{S}_{2}$ from being a subset of $\mathbb{S}_{1}$:

$d_{DH}(\mathbb{S}_{1},\mathbb{S}_{2}):=\sup_{s_{2}\in\mathbb{S}_{2}}\inf_{s_{1}\in\mathbb{S}_{1}}d(s_{1},s_{2}),$
(1)

where $d:\mathbb{R}^{n}\times\mathbb{R}^{n}\rightarrow\mathbb{R}_{+}$ is a
metric. For closed compact sets, $d_{DH}(\mathbb{S}_{1},\mathbb{S}_{2})=0$, if
and only if $\mathbb{S}_{2}\subseteq\mathbb{S}_{1}$. The Cartesian product of
$\mathbb{S}_{1}$ and $\mathbb{S}_{2}$ is denoted by
$\mathbb{S}_{1}\times\mathbb{S}_{2}$ and the Cartesian product of
$\mathbb{S}_{1},\cdots,\mathbb{S}_{N}$ is denoted by
$\prod_{i=1}^{N}\mathbb{S}_{i}$.

###  2.2 Zonotopes

A zonotope $\mathcal{Z}(c,G)$ is defined as $c\oplus
G\mathbb{B}_{p}\subset\mathbb{R}^{n}$, where $c\in\mathbb{R}^{n}$ is the
_center_ , the columns of $G\in\mathbb{R}^{n\times p}$ are the _generators_ ,
and $\mathbb{B}_{p}:=\\{b\in\mathbb{R}^{p}|||b||_{\infty}\leq 1\\}$. The order
of a zonotope is defined as $\frac{p}{n}$. Zonotopes are convenient to
manipulate with affine transformations and Minkowski sums:

$A\mathcal{Z}(c,G)+b=\mathcal{Z}(Ac+b,AG),$ (2a)
$\mathcal{Z}(c_{1},G_{1})\oplus\mathcal{Z}(c_{2},G_{2})=\mathcal{Z}(c_{1}+c_{2},[G_{1},G_{2}]).$
(2b)

The volume of a full-dimensional zonotope of order 1, also known as
parallelotope, is given as [17]:

$\text{Vol}(\mathcal{Z}(c,G))=2^{n}\text{det}(G^{T}G).$ (3)

The Cartesian product of zonotopes is:

$\prod_{i=1}^{N}\mathcal{Z}\left(c_{i},G_{i})=\mathcal{Z}([c_{1}^{T},\cdots,c_{N}^{T}]^{T},\text{Blk}(G_{1},\cdots,G_{N})\right).$

Zonotope order reduction methods over-approximate a zonotope by another one
with a smaller order. Several order reduction methods were reported in [25]
and [43]. In this paper, the Boxing method [8, 21] with order 1 is used to
over-approximate a given zonotope $\mathcal{Z}(c,G)$ by a hyper-box:

$\operatorname{{\text{Red}}}(\mathcal{Z}(c,G)):=\mathcal{Z}(c,\operatorname{{\text{Diag}}}(\sum_{i}{|g_{i}|})),$
(4)

where $g_{i}$ represents the $i$th column of $G$. Equation (4) is based on the
following properties of zonotopes:

$c-\sum_{i=1}^{p}{|g_{i}|}\leq\mathcal{Z}(c,G)\leq c+\sum_{i=1}^{p}{|g_{i}|}.$
(5)

######  Lemma 1 (Zonotope Containment [36]).

Given two zonotopes $\mathcal{Z}(c_{1},G_{1})$ and $\mathcal{Z}(c_{2},G_{2})$,
where $c_{1},c_{2}\in\mathbb{R}^{q}$ and $G_{1}\in\mathbb{R}^{q\times r}$,
$G_{2}\in\mathbb{R}^{q\times s}$, we have
$\mathcal{Z}(c_{1},G_{1})\subseteq\mathcal{Z}(c_{2},G_{2})$, if
$\exists\Gamma\in\mathbb{R}^{s\times r}$ and $\gamma\in\mathbb{R}^{s}$ s.t.

$G_{1}=G_{2}\Gamma,$ (6a) $c_{2}-c_{1}=G_{2}\gamma,$ (6b)
$\left\|[\Gamma,\gamma]\right\|_{\infty}\leq 1.$ (6c)

While Lemma 1 provides a sufficiency condition, it was shown in [36] that its
necessity gap, in a particular quantitative sense, is often small. The
following modified version of Lemma 1 enables us to independently scale each
generator of $G_{2}$.

######  Lemma 2 (Weighted Zonotope Containment).

Given two zonotopes $\mathcal{Z}(\bar{c}_{1},G_{1})$ and
$\mathcal{Z}(\bar{c}_{2},G_{2})$ (in the same format as Lemma 1), and a vector
$\alpha\in\mathbb{R}^{s}_{+}$, we have
$\mathcal{Z}(\bar{c}_{1},G_{1})\subseteq\mathcal{Z}(\bar{c}_{2},G_{2}\operatorname{{\text{Diag}}}(\alpha))$,
if the conditions in Lemma 1 hold while constraint (6c) changes to the
following element-wise inequality:

$[|\Gamma|,|\gamma|]\mathbb{1}_{s}\leq\alpha.$ (7)

Proof: see Appendix

######  Lemma 3.

(Directed Hausdorff Distance Computation [36]) Given two closed sets
$\mathbb{S}_{1}$ and $\mathbb{S}_{2}$ and $\|.\|_{\infty}$ as the underlying
metric, $d_{DH}(\mathbb{S}_{1},\mathbb{S}_{2})$ is the optimal value for the
following linear program:

$\displaystyle\underset{d}{\text{argmin}}\quad$ $\displaystyle d$ (8a) s.t.
$\displaystyle\mathbb{S}_{2}\subseteq\mathbb{S}_{1}\oplus
d\mathcal{Z}(0_{n},I_{n}),$ (8b) $\displaystyle d\geq 0.$ (8c)

###  2.3 Constrained Linear Discrete-Time Systems

A discrete-time (time-varying) linear system is characterized by the following
equation:

$x_{t+1}=A_{t}x_{t}+B_{t}u_{t}+d_{t},$ (9)

where $x_{t}\in\mathbb{X}_{t}\subseteq\mathbb{R}^{n}$,
$u_{t}\in\mathbb{U}_{t}\subseteq\mathbb{R}^{m}$, and
$d_{t}\in\mathbb{D}_{t}\subseteq\mathbb{R}^{n}$ are the state, control, and
disturbance at time $t\in\mathbb{N}$, respectively. The polytopic sets
$\mathbb{X}_{t}$, $\mathbb{U}_{t}$, and $\mathbb{D}_{t}$ define time varying
bounds and are assumed to be given. The matrices $A_{t}\in\mathbb{R}^{n\times
n}$ and $B_{t}\in\mathbb{R}^{n\times m}$ may be time dependent. If $A_{t}$,
$B_{t}$ and the sets $\mathbb{X}_{t},\mathbb{U}_{t}$, and $\,\mathbb{D}_{t}$
are all time-invariant, the system is linear time-invariant (LTI). If at least
one of them is time-variant, the system is linear time-variant (LTV) and we
are mainly interested in a time-limited response over a finite horizon
$h\in\mathbb{N}_{+}$ (i.e. $t\in\mathbb{N}_{h}$). Due to space limits, the
optimization problems and algorithms shown here are for the LTV class of
problems. The LTI versions can be seen as particular cases and can be
formulated easily.

A control policy $\mu$ is characterized by a set of functions
$\mu_{t}:\mathbb{X}_{t}\rightarrow\mathbb{U}_{t},t\in\mathbb{N}_{h}$, where
$h\in\mathbb{N}_{+}$ is a finite horizon. For infinite horizon, the policy
$\mu(.):\mathbb{X}\rightarrow\mathbb{U}$ is not a function of time.

######  Definition 1 (Finite-time Viable Sets).

For a given horizon $h\in\mathbb{N}_{+}$, a sequence of sets
$\Omega_{0},\Omega_{1},...,\Omega_{h}$ for system (9) is a sequence of
_viable_ sets, if for all $t\in\mathbb{N}_{h}$, $\Omega_{t}\subseteq X_{t}$
and there exists a policy $\mu_{t}$ such that $\Theta_{t}\subseteq U_{t}$,
where $\Theta_{t}:=\mu_{t}(\Omega_{t})$, and

$\forall t\in\mathbb{N}_{h-1},\forall x_{t}\in\Omega_{t},\forall d_{t}\in
D_{t}\Rightarrow x_{t+1}\in\Omega_{t+1}.$ (10)

$\Theta_{t}$ is called action set.

######  Definition 2 (Infinite-time Viable Set [5]).

A set $\Omega\subseteq X$ for an LTI system is an infinite-time viable set,
also known as robust control invariant (RCI) set, if there exists a control
policy $\mu$ such that $\Theta:=\mu(\Omega)\subseteq U$ and

$\forall t\in\mathbb{N},\forall x_{t}\in\Omega,\forall d_{t}\in D\Rightarrow
x_{t+1}\in\Omega.$ (11)

##  3 Problem Statement and Approach

Consider a network of linear subsystems that are dynamically coupled in the
following form:

$\begin{array}[]{ll}x_{i,t+1}=&A_{ii,t}x_{i,t}+B_{ii,t}u_{i,t}+d_{i,t}\\\
&+\displaystyle\sum_{j\neq i}{A_{ij,t}x_{j,t}}+\displaystyle\sum_{j\neq
i}{B_{ij,t}u_{j,t}},\end{array}$ (12)

where $x_{i,t}\in\mathbb{R}^{n_{i}}$, $u_{i,t}\in\mathbb{R}^{m_{i}}$, and
$d_{i,t}\in D_{i,t},D_{i,t}\subseteq\mathbb{R}^{n_{i}}$ are the state,
control, and the disturbance for subsystem $i\in\mathcal{I}$,
$\mathcal{I}=\mathbb{N}_{1:\eta}$, and $\eta$ is the number of subsystems.
Also, $D_{i,t}=\mathcal{Z}(\bar{d}_{i,t},G_{i,t}^{d})$, where
$\bar{d}_{i,t}\in\mathbb{R}^{n_{i}}$ and
$G_{i,t}^{d}\in\mathbb{R}^{n_{i}\times p_{i}^{d_{t}}}$ is the given
disturbance set for subsystem $i$. We denote the disturbance set for the
aggregate system as $D_{t}=\prod_{i\in\mathcal{I}}D_{i,t}$. The matrices
$A_{ii,t}\in\mathbb{R}^{n_{i}\times n_{i}}$ and
$B_{ii,t}\in\mathbb{R}^{n_{i}\times m_{i}}$ characterize the (time-variant)
internal dynamics of subsystem $i$. Also, $A_{ij,t}\in\mathbb{R}^{n_{i}\times
n_{j}}$ and $B_{ij,t}\in\mathbb{R}^{n_{i}\times m_{j}}$ characterize the
coupling effects of subsystem $j$ on subsystem $i$.

######  Remark 1.

If the sets $\mathcal{D}_{i,t}$ are not zonotopic, they can be over-
approximated by zonotopic sets. Furthermore, if the disturbance set is given
for the aggregated system, it may be decomposed into sets for each subsystem
$i$, such that the aggregated disturbance set is the subset of the Cartesian
product of the new sets for each subsystem.

We consider coupled constraints over state $(X_{t})$ and control input
$(U_{t})$. These sets can be polytopic sets or zonotopes. For the network of
subsystems in (12), the constraints are

$[x_{1,t}^{T},\cdots,x^{T}_{\eta,t}]^{T}\in X_{t},\hskip
8.53581pt[u^{T}_{1,t},\cdots,u^{T}_{\eta,t}]^{T}\in U_{t}.$ (13)

We are now ready to formulate the problem of finding decentralized viable sets
and their corresponding decentralized controllers $\mu_{i}(.)$ with respect to
the constraints.

######  Problem 1.

(Decentralized Finite-time Viable Sets) Given a network of interconnected
subsystems in the form of (12) and the horizon $h\in\mathbb{N}_{+}$, find sets
$\Omega_{i,t}$ and decentralized controllers $\mu_{i,t}(.),\forall
i\in\mathcal{I}$ and $t\in\mathbb{N}_{h}$ such that
$\prod_{i\in\mathcal{I}}\Omega_{i,t}\subseteq X_{t}$,
$\prod_{i\in\mathcal{I}}\Theta_{i,t}\subseteq U_{t}$, and

$\begin{array}[]{ll}&\forall x_{i,t}\in\Omega_{i,t},\forall
u_{j,t}\in\Theta_{j,t},(j\neq i),\forall d_{i,t}\in D_{i,t}\\\ &\Rightarrow
x_{i,t+1}\in\Omega_{i,t+1},\end{array}$ (14)

where $\Theta_{i,t}=\mu_{i,t}(\Omega_{i,t})$.

In essence, Problem 1 is a tube-based trajectory optimization problem with
reachability constraints for $h$ steps ahead, while taking into account all
possible disturbances to ensure guaranteed reachability. This is similar to a
robust MPC formulation with horizon $h$ when solved recursively, but without
an objective function because the goal can be achieved by reachability set
containment constraints. Thus, the objective function is optional. A terminal
condition can also be introduced to the MPC formulation for the purpose of
ensuring stability and recursive feasibility.

In addition to the previous problem, which focuses on the time-limited
response of an LTV system, we also consider the infinite-time response of an
LTI system:

######  Problem 2.

(Decentralized Infinite-time Viable Sets) For the particular case when each
subsystem in (12) is time invariant, find sets $\Omega_{i}$ and decentralized
controllers $\mu_{i}$(.) for $\forall i\in\mathcal{I}$, such that
$\prod_{i\in\mathcal{I}}\Omega_{i}\subseteq X$,
$\prod_{i\in\mathcal{I}}\Theta_{i}\subseteq U$, and

$\begin{array}[]{ll}&\forall x_{i,t}\in\Omega_{i},\forall
x_{j,t}\in\Omega_{j},\forall u_{j,t}\in\Theta_{j}(j\neq i),\forall d_{i,t}\in
D_{i}\\\ &\Rightarrow x_{i,t+1}\in\Omega_{i},\end{array}$ (15)

where $\Theta_{i}=\mu_{i}(\Omega_{i})$.

It is worth noting that the concept of the infinite-time viable sets and
Problem 2 can be extended to $T$-periodic systems where
$A_{t+T}=A_{t},B_{t+T}=B_{t},X_{t+T}=X_{t},U_{t+T}=U_{t},D_{t+T}=D_{t},\forall
t\in\mathbb{N}$. We omit studying this class of systems in this paper.

We wish to handle the coupled constraints (13) separately. To this goal, we
consider an additional step of decoupling the constraints. This allows each
subsystem to impose constraints independently of the others, while still
satisfying the coupled constraints when the subsystems are combined, therefore
allowing for distributed computations. Consider the following domain and range
for each controller:

$\mu_{i,t}(.):X_{i,t}\rightarrow U_{i,t},$ (16)

where $X_{i,t}$ and $U_{i,t}$ are admissible sets in the state space and
control space for subsystem $i$ at time $t$, respectively. However, since the
constraints are given in coupled forms, $X_{i,t}$ and $U_{i,t}$ are unknown
and need to be found by decomposing $X_{t}$ and $U_{t}$, respectively. We
address this in the following subproblem:

######  Subproblem 1 (Set Decomposition).

Given a set $X\subset\mathbb{R}^{n}$ in the form of a zonotope and a set of
integers $n_{i}$ ,$\forall i\in\mathcal{I}$ $(\sum_{i}n_{i}=n)$, find sets
$X_{i}=\mathcal{Z}(c_{i},G^{x}_{i})$ where $c_{i}\in\mathbb{R}^{n_{i}}$ and
$G^{x}_{i}\in\mathbb{R}^{n_{i}\times p^{x}_{i}},p^{x}_{i}\in\mathbb{N}_{+}$
with the maximum volume for the set $\prod_{i\in\mathcal{I}}X_{i}$, such that

$\forall x_{i}\in
X_{i}\Rightarrow[x_{1}^{T},x_{2}^{T},...,x_{\eta}^{T}]^{T}\in X.$

Solution: See Appendix.

Defining this subproblem for all time steps of $X_{t}$ and $U_{t}$ will lead
to full decomposition of constraints.

##  4 Assume-Guarantee Contracts

In this section, we formalize assume-guarantee contracts for one system and
provide details on the convex parameterization of contracts and controllers.

###  4.1 Definitions

######  Definition 3 (Assume-Guarantee Contract).

An assume-guarantee contract for system (9) is a pair
$\mathcal{C}=(\mathcal{A},\mathcal{G})$, where:

  * •

$\mathcal{A}$ is the assumption, which is a sequence of disturbance sets
$\mathcal{D}_{t},t\in\mathbb{N}_{h-1}$;

  * •

$\mathcal{G}$ is the guarantee, which is a sequence of tuples
$(\mathcal{X}_{t},\mathcal{U}_{t})$, where $\mathcal{X}_{t}$ and
$\mathcal{U}_{t}$ are two sets at time $t\in\mathbb{N}_{h}$ in the state space
and control space, respectively.

######  Definition 4 (Contract Validity).

A contract is _valid_ if its guarantee respects the system’s constraints
$\mathcal{X}_{t}\subseteq X_{t},\forall
t\in\mathbb{N}_{h},\mathcal{U}_{t}\subseteq U_{t},\forall
t\in\mathbb{N}_{h-1}$.

######  Definition 5 (Contract Satisfiability).

A valid contract is _satisfiable_ if it is possible to find a control policy
and viable sets such that $\Omega_{t}\subseteq\mathcal{X}_{t},\forall
t\in\mathbb{N}_{h},\Theta_{t}\subseteq\mathcal{U}_{t},\forall
t\in\mathbb{N}_{h-1}$.

###  4.2 Finite Horizon Contract Satisfiability

For a single system with a given disturbance bound, we show that a satisfiable
contract can be found using convex programs, which encode a specific form of
control policies.

######  Theorem 1.

Given an LTV system in the form (9) with the bounded disturbance set
$D_{t}=\mathcal{Z}(\bar{d}_{t},G^{d}_{t})$, where
$\bar{d}_{t}\in\mathbb{R}^{n}$ and $G^{d}_{t}\in\mathbb{R}^{n\times p_{t}}$, a
finite horizon contract is satisfiable if $\exists k\in\mathbb{N}$, vectors
$\bar{\mathrm{x}}_{t}\in\mathbb{R}^{n},\bar{\mathrm{u}}_{t}\in\mathbb{R}^{m}$,
and matrices $T_{t}\in\mathbb{R}^{n\times l_{t}}$ and
$M_{t}\in\mathbb{R}^{m\times l_{t}}$, where $l_{0}=k$ and $l_{t\neq
0}=k+\sum_{\hat{t}=0}^{t-1}{p_{\hat{t}}}$ such that the following relations
hold:

$[A_{t}T_{t}+B_{t}M_{t},G^{d}_{t}]=T_{t+1},\quad\forall t\in\mathbb{N}_{h-1},$
(17a)
$A_{t}\bar{\mathrm{x}}_{t}+B_{t}\bar{\mathrm{u}}_{t}+\bar{d}_{t}=\bar{\mathrm{x}}_{t+1},\quad\forall
t\in\mathbb{N}_{h-1},$ (17b) $\mathcal{Z}(\bar{\mathrm{x}}_{t},T_{t})\subseteq
X_{t},\quad\forall t\in\mathbb{N}_{h},$ (17c)
$\mathcal{Z}(\bar{\mathrm{u}}_{t},M_{t})\subseteq U_{t},\quad\forall
t\in\mathbb{N}_{h-1}.$ (17d)

Then $\Omega_{t}=\mathcal{Z}(\bar{\mathrm{x}}_{t},T_{t}),t\in\mathbb{N}_{h}$
is the sequence of viable sets for horizon $h$ and
$\Theta_{t}=\mathcal{Z}(\bar{\mathrm{u}}_{t},M_{t})$ is the sequence of action
sets. Moreover, the controller $\mu_{t}(x_{t})$ can be computed by:

$\mu_{t}(x_{t})=\bar{\mathrm{u}}_{t}+M_{t}\zeta(x),$ (18a)
$x_{t}=\bar{\mathrm{x}}_{t}+T_{t}\zeta(x),$ (18b)
$\zeta\in\mathbb{B}_{l_{t}},$ (18c)

######  Proof 1.

The proof is by construction. Substituting (18) in (9) yields:

$x_{t+1}\in
A_{t}(T_{t}\zeta+\bar{\mathrm{x}}_{t})+B_{t}(M_{t}\zeta+\bar{\mathrm{u}}_{t})+\bar{d}_{t}\oplus
G^{d}_{t}\mathbb{B}_{p_{t}},$ (19)

where the right hand side set is equal to

$\\{A_{t}\bar{\mathrm{x}}_{t}+B_{t}\bar{\mathrm{u}}_{t}+\bar{d}_{t}\\}\oplus[A_{t}T_{t}+B_{t}M_{t},G^{d}_{t}]\mathbb{B}_{l_{t}+p_{t}}.$
(20)

By definition, we know that

$x_{t+1}\in\mathcal{Z}(\bar{\mathrm{x}}_{t+1},T_{t+1}).$ (21)

By using the Minkowski sum property of zonotopes shown in (2b) and by assuming
that the two sets shown in (20) and (21) are identical, it is straightforward
to reach (17a) and (17b). The next two constraints (17c) and (17d) are
imposing the constraints at each time step.

Using the linear encoding for the set containment problem proposed in [36],
containment constraints (17c) and (17d) can be encoded as linear constraints
with respect to $T_{t},M_{t},\bar{x}_{t},$ and $\bar{u}_{t}$. The cost
function can be chosen depending on the application. We typically choose to
minimize the summation of Frobenious norms of $T_{t}$ for $t\in\mathbb{N}_{h}$
as a heuristic to minimize the size of the viable sets.

######  Remark 2.

Note that the order of zonotope $\Omega_{t}$ is increasing at each time step.
This makes the number of variables and constraints in the program grow
quadratically with $h$. The complexity can be decreased by fixing the number
of columns in matrices $T_{t}$ and $M_{t}$ at $k$ and by changing equation
(17a) to

$[A_{t}T_{t}+B_{t}M_{t},G^{d}_{t}]=[0_{n\times p_{t}},T_{t+1}].$ (22)

However, this modification leads to a more conservative computation and may
cause infeasibility.

###  4.3 Infinite Horizon Contract Satisfiability

Inspired by the method in [32], we provide a linear programming approach to
compute RCI sets.

######  Theorem 2.

Given an LTI system in the form (9) with the bouned disturbance set
$D=\mathcal{Z}(\bar{d},G^{d})$, where $\bar{d}\in\mathbb{R}^{n}$ and
$G^{d}\in\mathbb{R}^{n\times p}$, an infinite horizon contract is satisfiable
if $\exists k\in\mathbb{N},\beta\in[0,1)$, vectors
$\bar{\mathrm{x}}\in\mathbb{R}^{n},\bar{\mathrm{u}}\in\mathbb{R}^{m}$, and
matrices $T\in\mathbb{R}^{n\times k}$, $M\in\mathbb{R}^{m\times k}$, and
$E\in\mathbb{R}^{n\times p}$, such that the following relations hold:

$[AT+BM,G^{d}]=[E,T],$ (23a) $\mathcal{Z}(0,E)\subseteq\mathcal{Z}(0,\beta
G^{d}),$ (23b) $A\bar{\mathrm{x}}+B\bar{\mathrm{u}}+\bar{d}=\bar{\mathrm{x}},$
(23c) $\mathcal{Z}(\bar{\mathrm{x}},\dfrac{1}{1-\beta}T)\subseteq X,$ (23d)
$\mathcal{Z}(\bar{\mathrm{u}},\dfrac{1}{1-\beta}M)\subseteq U.$ (23e)

Then $\Omega=\mathcal{Z}(\bar{\mathrm{x}},(1-\beta)^{-1}T)$ is a RCI set and
the action set is $\Theta=\mathcal{Z}(\bar{\mathrm{u}},(1-\beta)^{-1}M)$.
Furthermore, the controller can be computed by

$\mu(x)=\bar{\mathrm{u}}+(1-\beta)^{-1}M\zeta(x),$ (24a)
$x=\bar{\mathrm{x}}+(1-\beta)^{-1}T\zeta(x),$ (24b) $\zeta\in\mathbb{B}_{k}.$
(24c)

######  Proof 2.

Substituting policy (24) in (9), we obtain relation (23a). In order to prove
invariance, we observe that:

$(AT+BM)\mathbb{B}_{k}\oplus\mathcal{Z}(0,G^{d})\subseteq
T\mathbb{B}_{k}\oplus\mathcal{Z}(0,\beta G^{d}),$ (25)

and

$A\bar{\mathrm{x}}+B\bar{\mathrm{u}}+\bar{d}=\bar{\mathrm{x}}.$ (26)

We subtract $\mathcal{Z}(0,\beta G^{d})$ from both sides in the Pontryagin
difference sense [24]. This is a valid operation when both sides are convex
polytopes. We omit the proof as it is based on the properties of support
functions [33] of convex sets. We arrive at:

$(AT+BM)\mathbb{B}_{k}\oplus(1-\beta)\mathcal{Z}(0,G^{d})\subseteq
T\mathbb{B}_{k}.$ (27)

By multiplying both sides of (27) by $\dfrac{1}{1-\beta}$, we reach the
conclusion that $\mathcal{Z}(\bar{\mathrm{x}},(1-\beta)^{-1}T)$ is a RCI set
with $\Theta=\mathcal{Z}(\bar{\mathrm{u}},(1-\beta)^{-1}M)$, and the proof is
complete.

Similar to (17), there exists a sufficient linear encoding for (23) to find
the viable sets and action sets for fixed $k$ and $\beta$. The feasibility of
the linear program implies the satisfiability of the contract.

######  Remark 3.

We can simplify Theorem 2 by assuming $E=0_{n\times p}$ and $\beta=0$. As a
result, there is no need for constraint (23b). However, this assumption
increases conservativeness and may lead to infeasibility.

Note that, in both Theorem 1 and Theorem 2, the structures of matrices $T_{t}$
and $M_{t}$ (or $T$ and $M$) depend on a constant $k$ that acts like a hyper-
parameter for our proposed linear program. Therefore, changing $k$ leads to a
different result. This allows iterations over different $k$ to find a feasible
solution for $T_{t}$ and $M_{t}$ (or $T$ and $M$). However, finding a feasible
solution with the smallest possible value for $k$ is preferable, since it
decreases the number of variables in the optimization problem.

##  5 Composition of Parametric Assume-Guarantee Contracts

In this section, we focus on the network of coupled systems (12) and provide
the first steps for the solutions to Problems 1 and 2. First, we decouple the
subsystems by viewing the coupling effects of other subsystems as
disturbances, which are called augmented disturbances $d_{i,t}^{aug}$ and are
equal to:

$d_{i,t}^{aug}:=\sum_{j\neq i}A_{ij,t}x_{j,t}+\sum_{j\neq
i}B_{ij,t}u_{j,t}+d_{i,t},$ (28)

where $i$ is the subsystem’s index and $t$ is the time step. Next, for each
subsystem $i\in\mathcal{I}$, we define one AG contract denoted by
$\mathcal{C}_{i}$. The collection of these contracts forms a set of AG
contracts $\mathcal{C}=\\{\mathcal{C}_{i}|i\in\mathcal{I}\\}$ for a network of
coupled subsystems. From (28), it can be seen that, for subsystem $i$, the
assumption set over the disturbance space, which is denoted by $W_{i,t}$, can
be defined as a function of the guarantees of other subsystems
$(\mathcal{X}_{j,t},\mathcal{U}_{j,t})(j\neq i)$ as follows:

$W_{i,t}:=\bigoplus_{j\neq i}A_{ij,t}\mathcal{X}_{j,t}\oplus\bigoplus_{j\neq
i}B_{ij,t}\mathcal{U}_{j,t}\oplus D_{i,t}.$ (29)

However, unlike the single-system case, in which by using Theorem 1 (or
Theorem 2 for the LTI case), we could directly compute the guarantee for a
given disturbance set, in the case of interconnected subsystems, finding a set
of satisfiable contracts is not straightforward. This is due to the fact that
the guarantee of one subsystem affects the assumptions of other subsystems as
a result of treating the coupling effects as disturbance. This is known as the
circularity problem of AG contracts.

Also, it is worth noting that the contracts are common knowledge among all
subsystems because each subsystem has to be able to compute its assumed
disturbance using (29). This, in fact, helps in the compositional computation
of the synthesis problem.

Next, we introduce a correctness criterion to break the circularity of AG
contracts and parametric contracts to search over contract sets to find a set
of satisfiable contracts. Also, we propose a potential function as a
quantitative indicator of how far a set of given contracts is from correct
composition. Finally, we introduce a set of sufficient constraints for the
validity of parametric contracts.

###  5.1 Composition Correctness

From (28), it can be shown that after finding the viable sets $\Omega_{i,t}$
and action sets $\Theta_{i,t}$ for all $i\in\mathcal{I}$, the actual augmented
disturbance set for subsystem $i$ at time $t$, denoted by $D_{i,t}^{aug}$, can
be obtained as:

$D^{aug}_{i,t}:=\bigoplus_{j\neq i}A_{ij,t}\Omega_{j,t}\oplus\bigoplus_{j\neq
i}B_{ij,t}\Theta_{j,t}\oplus D_{i,t}.$ (30)

If the actual experienced disturbance set for every subsystem becomes a subset
of its assumed disturbance set, the circularity of AG contracts is broken. We
say that the composition of the contracts is correct. The formal definition of
correctness is given as follows:

######  Definition 6 (Composition Correctness).

Consider a set of locally valid assume-guarantee contracts
$\mathcal{C}_{i}=(\mathcal{A}_{i},\mathcal{G}_{i})$, $i\in\mathcal{I}$. The
composition is correct if the following relation holds:

$D_{i,t}^{aug}\subseteq W_{i,t},\forall i\in\mathcal{I},\forall
t\in\mathbb{N}_{h-1}$ (31)

However, in order to make compositional computations easier, relation (31) is
replaced with the following more conservative version:

$\Omega_{i,t}\subseteq\mathcal{X}_{i,t},\Theta_{i,t}\subseteq\mathcal{U}_{i,t},\forall
i\in\mathcal{I},\forall t\in\mathbb{N}_{h-1},$ (32)

When (32) holds, from (29) and (30), it is evident that (31) holds as well.
Our characterization of correctness is Boolean, i.e., it indicates whether or
not the contract composition is correct. Nevertheless, having a quantitative
measure that describes how far a set of given contracts is from correct
composition is desirable. Such measurement can be very helpful in directing us
toward correct contract composition. The following “potential function” does
exactly that by assigning a score to a set of contracts.

######  Definition 7 (Potential Function).

Given a set of contracts $\mathcal{C}=\\{\mathcal{C}_{i}|i\in\mathcal{I}\\}$,
its potential function is

$\mathcal{V}(\mathcal{C})=\sum_{i\in\mathcal{I}}\mathcal{V}_{i}(\mathcal{C}),$
(33)

where $\mathcal{V}_{i}(\mathcal{C})$ is defined as follows:

$\mathcal{V}_{i}(\mathcal{C}):=\sum_{t\in\mathbb{N}_{h}}d_{DH}(\mathcal{X}_{i,t},\Omega_{i,t})+\sum_{t\in\mathbb{N}_{h-1}}d_{DH}(\mathcal{U}_{i,t},\Theta_{i,t}).$
(34)

The potential function is equal to the sum of the directed Hausdorff distances
between the sets in (32). Given that (32) leads to (31) and the potential
function is zero if and only if (32) holds, it can be inferred that when the
potential function is zero, the composition of the contracts is correct,
implying that each system anticipates a larger set of disturbance than it will
actually experience. Thus, this quantitative characterization of correct
composition is sound. Fig 1 shows an example of the sets in the state space of
subsystem $i$ at a given time $t$. A similar diagram can be drawn for the
control space to show the other component of the potential function.

Figure 1: The admissible set $X_{i,t}$ (red), guarantee set in the state space
$\mathcal{X}_{i,t}$ (blue), and the viable set $\Omega_{i,t}$ (green) for a
2-dimensional example subsystem $i$ at time $t$. The first component of the
potential function in (34) at time step $t$ (i.e.
$d_{DH}(\mathcal{X}_{i,t},\Omega_{i,t})$) is denoted by $\nu_{i,t}^{x}$. If
$\nu_{i,t}^{x}=0$, then $\Omega_{i,t}\subseteq\mathcal{X}_{i,t}$. Our goal is
to find a set of contracts such that
$\Omega_{i,t}\subseteq\mathcal{X}_{i,t}\subseteq X_{i,t}$ for all $i,t$.

###  5.2 Parametric Contracts

To be able to search over the contract sets, we introduce parametric sets. The
baseline sets are zonotopes $\mathcal{Z}(\bar{c}^{x}_{i,t},C^{x}_{i,t})$ and
$\mathcal{Z}(\bar{c}^{u}_{i,t},C^{u}_{i,t})$, for all $i\in\mathcal{I}$ and
$t\in\mathbb{N}_{h}$, where $\bar{c}^{x}_{i,t}\in\mathbb{R}^{n_{i}}$,
$\bar{c}^{u}_{i,t}\in\mathbb{R}^{m_{i}}$, and
$C^{x}_{i,t}\in\mathbb{R}^{n_{i}\times\zeta^{x}_{i,t}}$,
$C^{u}_{i,t}\in\mathbb{R}^{m_{i}\times\zeta^{u}_{i,t}}$. The baseline sets
need to be defined by the user in advance. One option is to select them as the
viable sets and action sets for each subsystem, respectively, while ignoring
the couplings to other subsystems. The parametric sets are denoted by
$\mathcal{X}_{i,t}(\alpha^{x}_{i,t})$ and
$\mathcal{U}_{i,t}(\alpha^{u}_{i,t})$, and are defined as follows:

$\mathcal{X}_{i,t}(\alpha^{x}_{i,t}):=\mathcal{Z}(\bar{c}^{x}_{i,t},C^{x}_{i,t}\operatorname{{\text{Diag}}}(\alpha^{x}_{i,t})),$
(35a)
$\mathcal{U}_{i,t}(\alpha^{u}_{i,t}):=\mathcal{Z}(\bar{c}^{u}_{i,t},C^{u}_{i,t}\operatorname{{\text{Diag}}}(\alpha^{u}_{i,t})),$
(35b)

where parameters $\alpha^{x}_{i,t}\in\mathbb{R}_{+}^{\zeta^{x}_{i,t}}$ and
$\alpha^{u}_{i,t}\in\mathbb{R}_{+}^{\zeta^{u}_{i,t}}$ are vectors with non-
negative real entries that scale the generators’ columns of the baseline sets.
Also, the set of all parameters is denoted by
$\alpha=\\{\alpha^{x}_{i,t},\alpha^{u}_{i,t}\\}_{\forall i,t}$. Finally, a
parametric contract $\mathcal{C}_{i}$ has its guarantee in the form of (35) in
the state and control space, respectively, and its assumption can be directly
derived by replacing (35) in (29), which is shown by $W_{i,t}(\alpha)$.

###  5.3 Parametric Potential Function

Here we introduce parametric potential functions, which quantify how close a
set of contracts/parameters is from composition correctness when the contracts
are given in parametric form.

######  Definition 8 (Parametric Potential Function).

Given a set of parameters $\alpha$, the parametric potential function is
defined as:

$\mathcal{V}(\alpha)=\sum_{i\in\mathcal{I}}\mathcal{V}_{i}(\alpha),$ (36)

where

$\begin{array}[]{ll}\mathcal{V}_{i}(\alpha):=&\displaystyle\sum_{t\in\mathbb{N}_{h}}d_{DH}(\mathcal{X}_{i,t}(\alpha^{x}_{i,t}),\Omega_{i,t})+\\\
&\displaystyle\sum_{t\in\mathbb{N}_{h-1}}d_{DH}(\mathcal{U}_{i,t}(\alpha^{u}_{i,t}),\Theta_{i,t}).\end{array}$
(37)

Note that the parametric potential function is zero if and only if

$\Omega_{i,t}\subseteq\mathcal{X}_{i,t}(\alpha^{x}_{i,t})\quad\text{and}\quad\Theta_{i,t}\subseteq\mathcal{U}_{i,t}(\alpha^{u}_{i,t}).$
(38)

Thus the level set of the potential function at zero equals the set of
parameters that correspond to the correct contract composition. These
parameters are referred to as correct parameters. Using the linear program
(8), the optimization problem for computing $\mathcal{V}_{i}(\alpha)$ is as
follows:

$\displaystyle\mathcal{V}_{i}(\alpha)=\min_{\bar{\mathrm{x}}^{i},T^{i},\bar{\mathrm{u}}^{i},M^{i},d_{t}^{x},d_{t}^{u}}\begin{aligned}
&\sum_{t\in\mathbb{N}_{h}}d^{x}_{t}+\sum_{t\in\mathbb{N}_{h-1}}d^{u}_{t}\end{aligned}$
subject to
$\displaystyle[A_{ii,t}T^{i}_{t}+B_{ii,t}M^{i}_{t},\mathcal{W}_{i,t}^{aug}]=T^{i}_{t+1},\quad\forall
t\in\mathbb{N}_{h-1}$ (39a) $\displaystyle
A_{ii,t}\bar{\mathrm{x}}^{i}_{t}+B_{ii,t}\bar{\mathrm{u}}^{i}_{t}+\bar{d}_{i,t}^{aug}=\bar{\mathrm{x}}^{i}_{t+1},\quad\forall
t\in\mathbb{N}_{h-1}$ (39b)
$\displaystyle\mathcal{Z}(\bar{\mathrm{x}}^{i}_{t},T^{i}_{t})\subseteq\
\mathcal{X}_{i,t}(\alpha^{x}_{i,t})\oplus\mathcal{Z}(0_{n_{i}},d^{x}_{t}I_{n_{i}}),\forall
t\in\mathbb{N}_{h}$ (39c)
$\displaystyle\mathcal{Z}(\bar{\mathrm{u}}^{i}_{t},M^{i}_{t})\subseteq\
\mathcal{U}_{i}(t,\alpha^{u}_{i,t})\oplus\mathcal{Z}(0_{m_{i}},d^{u}_{t}I_{m_{i}}),\quad\forall
t\in\mathbb{N}_{h-1}$ (39d) $\displaystyle d_{t}^{x}\geq 0,\hskip
5.69054pt\forall t\in\mathbb{N}_{h},$ (39e) $\displaystyle d_{t}^{u}\geq
0,\hskip 5.69054pt\forall t\in\mathbb{N}_{h-1},$ (39f)

where $\bar{\mathrm{x}}^{i},T^{i},\bar{\mathrm{u}}^{i}$, and $M^{i}$ are sets
containing all $\bar{\mathrm{x}}^{i}_{t},T^{i}_{t},\bar{\mathrm{u}}^{i}_{t}$,
and $M^{i}_{t}$, respectively. Constraints (39a) and (39b) come from Theorem
1, which enforce viability conditions (must be replaced by the constraints in
Theorem 2 in the LTI case), where $\bar{d}_{i,t}^{aug}$ and
$\mathcal{W}_{i,t}^{aug}$ are the center and generator of the assumed
disturbance set for the subsystem $i$, respectively:

$\mathcal{Z}(\bar{d}_{i,t}^{aug},\mathcal{W}_{i,t}^{aug}):=\operatorname{{\text{Red}}}(W_{i,t}(\alpha)),\forall
t\in\mathbb{N}_{h-1}.$ (40)

Since the assumed disturbance set for each subsystem is known for a given
$\alpha$, zonotope order reduction methods may be used to over-approximate the
disturbance set in order to reduce computational complexity. Also, scalar
variables $d^{x}_{t}$ and $d^{u}_{t}$, constraints (39c)-(39f) and the
objective function originate from Lemma 3 in order to compute the sum of
Directed Hausdorff distances over all time steps on both state and control
spaces. The original problem of computing all of the directed Hausdorff
distances is a multi-objective optimization problem, i.e. one objective
function for each Directed Hausdorff distance computation. However, since the
optimal points of the optimization problems are not conflicting with each
other, we may add the objective functions and ensure that the outcome is the
same as when each optimization problem is solved independently. Note that
constraints (39c) and (39d) can be encoded into linear constraints using Lemma
2. The following theorem is the main result of this section:

######  Theorem 3.

(Convexity of parametric potential function) Using our parameterization (35)
and the linear encoding with containment introduced in Lemma 1 and Lemma 2,
the parametric potential function $\mathcal{V}(\alpha)$ is a convex piecewise
affine function. The set of correct parameters, which is the level set of
$\mathcal{V}(\alpha)$ at zero, is also a convex set.

######  Proof 3.

As shown in (39), each $\mathcal{V}_{i}(\alpha)$ is formulated in a linear
program, implying that each $\mathcal{V}_{i}(\alpha)$ is a convex piecewise
affine function [4]. Also, because the summation of convex piecewise affine
functions remains convex and piecewise affine, $\mathcal{V}(\alpha)$ is also a
convex piecewise affine function. Since the level set of a convex function is
a convex set, the set of correct parameters is a convex set.

###  5.4 Validity

A set of desirable contracts needs to be _compositionally correct_ and
individually _satisfiable_ and _valid_ , as indicated in Section 4. The set of
parameters that results in composition correctness and satisfiability does not
necessarily lead to validity as well. In other words, in Fig. 1, it may be
possible that $\Omega_{i,t}\subseteq\mathcal{X}_{i,t}$ but
$\Omega_{i,t}\not\subseteq X_{i,t}$. In this subsection, the focus is on
finding the set of valid parameters, which are the set of parameters that
corresponds to a valid contract for each subsystem. A contract is valid if and
only if:

$\mathcal{X}_{i,t}(\alpha^{x}_{i,t})\subseteq X_{i,t},\forall
t\in\mathbb{N}_{h},$ (41a) $\mathcal{U}_{i,t}(\alpha^{u}_{i,t})\subseteq
U_{i,t},\forall t\in\mathbb{N}_{h-1}.$ (41b)

######  Theorem 4 (Convexity of valid parameters set).

Using parameterization (35) and the linear encoding with containment
introduced in Lemma 2, the set of valid parameters derived from (41) is a
convex set.

######  Proof 4.

The set of valid parameters is the feasible region of the linear encoding of
(41). Therefore, it is a convex set [4].

Notably, in a centralized method, we may avoid decoupling constraints via the
solution for Subproblem 1 and instead impose them by
$\prod_{i\in\mathcal{I}}\mathcal{X}_{i,t}(\alpha^{x}_{i,t})\subseteq
X_{t},\forall t\in\mathbb{N}_{h}$ and
$\prod_{i\in\mathcal{I}}\mathcal{U}_{i,t}(\alpha^{u}_{i,t})\subseteq
U_{t},\forall t\in\mathbb{N}_{h-1}$.

##  6 Compositional Synthesis and Computations

Two approaches based on parametric assume-guarantee contracts are proposed in
this section to solve Problems 1 and 2. The first provides a single
centralized optimization to identify a set of desirable contracts,
decentralized viable sets, and decentralized controllers. The second method
uses the notion of parametric potential function to achieve the same task in a
compositional fashion.

###  6.1 Single Convex Program

By combining all of the encodings presented so far, we obtain the following
centralized linear program:

$\displaystyle\Omega,\Theta=\underset{\bar{\mathrm{x}}^{i}_{t},T^{i}_{t},\bar{\mathrm{u}}^{i}_{t},M^{i}_{t},\alpha,\mathrm{d}_{i},\mathrm{G}_{i}}{\text{argmin}}\begin{aligned}
&\sum_{t\in\mathbb{N}_{h-1}}\sum_{i\in\mathcal{I}}{\text{sum}(\alpha_{i,t}^{x})}\end{aligned}$
subject to
$\displaystyle[A_{ii,t}T^{i}_{t}+B_{ii,t}M^{i}_{t},\mathrm{G}_{i,t}]=T^{i}_{t+1},\forall
t\in\mathbb{N}_{h-1},\forall i\in\mathcal{I}$ (42a) $\displaystyle
A_{ii,t}\bar{\mathrm{x}}^{i}_{t}+B_{ii,t}\bar{\mathrm{u}}^{i}_{t}+\mathrm{d}_{i,t}=\bar{\mathrm{x}}^{i}_{t+1},\forall
t\in\mathbb{N}_{h-1},\forall i\in\mathcal{I}$ (42b)
$\displaystyle\mathcal{Z}(\mathrm{d}_{i,t},\mathrm{G}_{i,t})=W_{i,t}(\alpha),\forall
t\in\mathbb{N}_{h-1},\forall i\in\mathcal{I}$ (42c)
$\displaystyle\prod_{i\in\mathcal{I}}\mathcal{Z}(\bar{\mathrm{x}}^{i}_{t},T_{t}^{i})\subseteq
X_{t},\quad\forall t\in\mathbb{N}_{h}$ (42d)
$\displaystyle\prod_{i\in\mathcal{I}}\mathcal{Z}(\bar{\mathrm{u}}^{i}_{t},M_{t}^{i})\subseteq
U_{t},\quad\forall t\in\mathbb{N}_{h-1}$ (42e)
$\displaystyle\mathcal{Z}(\bar{\mathrm{x}}_{t}^{i},T^{i}_{t})\subseteq\mathcal{X}_{i,t}(\alpha^{x}_{i,t}),\quad\forall
t\in\mathbb{N}_{h-1},\forall i\in\mathcal{I}$ (42f)
$\displaystyle\mathcal{Z}(\bar{\mathrm{u}}^{i}_{t},M^{i}_{t})\subseteq\mathcal{U}_{i,t}(\alpha^{u}_{i,t}),\quad\forall
t\in\mathbb{N}_{h-1},\forall i\in\mathcal{I}$ (42g)
$\displaystyle\alpha^{x}_{i,t},\alpha^{u}_{i,t}\geq 0,\hskip 11.38109pt\forall
t\in\mathbb{N}_{h-1},\forall i\in\mathcal{I},$ (42h)

where
$\Omega=\\{\Omega_{i}|\Omega_{i}=\mathcal{Z}(\bar{\mathrm{x}}^{i}_{0},T^{i}_{0}),...,\mathcal{Z}(\bar{\mathrm{x}}^{i}_{h},T^{i}_{h}),\forall
i\in\mathcal{I}\\}$ and
$\Theta=\\{\Theta_{i}|\Theta_{i}=\mathcal{Z}(\bar{\mathrm{u}}^{i}_{0},M^{i}_{0}),...,\mathcal{Z}(\bar{\mathrm{u}}^{i}_{h-1},M^{i}_{h-1}),\forall
i\in\mathcal{I}\\}$ are decentralized viable sets and action sets,
respectively. The constraints (42a) and (42b) are sufficient constraints for
viable sets introduced in Theorem 1 (must be replaced with the constraints in
Theorem 2 for Problem 2). The constraints on the state and control input are
imposed in (42d) and (42e). They are applied over the aggregated system in
order to reduce the conservatism caused by decoupling the coupled constraints.
The constraint (42c) calculates the assumed disturbance set from (29). We
cannot use zonotope order reduction techniques on this formulation since doing
so will result in a non-convex encoding. The constraints (42f) and (42g) are
also the proposed requirements suggested by the correctness criterion (38).
Lemma 2 may be used to encode these two constraints into linear constraints.
It is worth noting that there is no need to provide validity constraints for
parametric AG contracts because constraints (42d) and (42e) suffice.
Furthermore, the objective function is ad-hoc and it may be selected as shown
in (42), which is a heuristic strategy for decreasing the volume of the viable
sets. The constant $k$ (introduced in Theorem 1) must be established before
solving the linear program in (42). Our strategy is to start with a small
initial $k$ and solve (42), then increase it by one unit until feasibility is
attained.

Note that this method is sound, because the correctness criterion is enforced
in the process by using zonotope containment constraints. As a result, the
output is correct-by-construction. Also, as $k$ goes to infinity, our approach
is theoretically able to find a set of correct parameters, if one exists.

###  6.2 Compositional Approach

The centralized method proposed in the previous section suffers from the curse
of dimensionality. Despite the fact that the control is decentralized, solving
a single, very large linear program is impractical for large-scale systems.
The computational complexity is driven by two main factors: (i) the large
number of variables and constraints in the single optimization problem (42);
(ii) the order of the zonotope
$\mathcal{Z}(\bar{d}_{i,t}^{aug},D_{i,t}^{aug})$ in the constraint (42c),
which becomes very large when the number of neighboring subsystems is large
due to the Minkowski sum computations in (29), which makes (42a) infeasible
for thin $T_{t}^{i}$ matrices. Also, because of the existing recursivity in
(42a), the widths of the matrices $T_{t}^{i}$ increase with time steps.

The main contribution of this paper is a compositional method for the
computation of the viable sets. In this subsection, we show that, by using our
proposed parameterized sets (35) and the convex potential function (36), the
single optimization problem in the previous subsection can be transformed into
a number of iterative small linear optimization problems with convergence
guarantees. Also, following (39), we use zonotope order reduction methods to
lower the order of the augmented disturbance sets, which substantially speeds
up the calculation.

As explained before, when the potential function is zero, the parameters are
inside the correct set of parameters and we have a set of compositionally
correct cont,racts. Thus, the goal is to determine $\alpha^{*}$ such that:

$\mathcal{V}^{*}=\underset{\alpha}{\text{min
}}\sum_{i\in\mathcal{I}}\mathcal{V}_{i}(\alpha),$ (43)

where $\mathcal{V}_{i}(\alpha)$ is the potential function that corresponds to
subsystem $i$. Each subsystem concurrently computes its own portion of the
potential function by (39). When the optimal value of the above optimization
problem is greater than zero $(\mathcal{V}^{*}>0)$, the optimal parameters are
not in the set of correct parameters. In this case, the algorithm raises the
hyper-parameter $k$ (defined in Theorem 1) by one unit and repeats the same
procedure, until it either hits a user-defined maximum threshold for $k$ (no
solution found) or it finds parameters with zero value of the potential
function. We solve the aforementioned optimization problem using gradient
descent. We start with a random initial guess for $\alpha$, and then update
the parameters using:

$\alpha\leftarrow\alpha-\delta\sum_{i\in\mathcal{I}}\nabla_{\alpha}\mathcal{V}_{i}(\alpha),$
(44)

where $\delta$ is the step size. We use duality to calculate the gradient of
the potential function. The optimal dual variable of a constraint in a linear
program is equal to the derivative of the objective function with respect to
the right hand side of the corresponding constraint, according to sensitivity
analysis of linear programs [4]. Using this property and the chain rule, we
can compute $\nabla_{\alpha}\mathcal{V}_{i}(\alpha)$ after solving the linear
program proposed for $\mathcal{V}_{i}(\alpha)$. In our implementation, we used
the built-in function in Gurobi [18] to find the optimal dual variables of the
constraints that contain the parameters. It should be noted that in (40),
additional tighter zonotope order reduction approaches can also be used to
further minimize conservatism. In our implementation, we used the Boxing
method to keep the chain rule simpler.

There are $\eta$ gradients for each point in the parameter space, where $\eta$
is the number of subsystems. In other words, each subsystem selfishly proposes
a direction that best suits it at each point in the parameters space. By
calculating the sum of these directions, we can determine which direction is
the best for the aggregated system. This is because the suggested objective
function is in the form of $\sum_{i\in\mathcal{I}}\mathcal{V}_{i}(\alpha)$,
and its derivative with respect to $\alpha$ is equal to the sum of the
gradients computed by each subsystem. Finally, we stop updating the parameters
when the value of the potential function reaches a plateau and there is no
further improvement in the objective function after the set of parameters is
updated.

As stated in Subsection 5.4, the existence of a zero potential function does
not imply that the set of parameters is also valid. There are two ways to
address this issue. First, we can add the directed Hausdorff distance between
parameterized sets and admissible sets to the objective function in (39), as
we did for correctness. Then, we can determine a set of parameters that will
cause the new objective function to reach zero. Second, we can limiting the
domain of the potential function to the set of valid parameters. To avoid
leaving the valid set of parameters after each update, the parameters must be
projected to the set of valid parameters. In this case, the potential function
remains convex because the set of valid parameters is convex. The second
method projects the parameters onto the target set in a single step, whereas
the first requires many iterations. Here, we pick the second option due to its
speed. The projection is performed via the following optimization problem:

$\displaystyle\alpha^{x^{p}}_{i,t}=\underset{\alpha^{p}}{\text{argmin }}$
$\displaystyle||\alpha^{x^{p}}_{i,t}-\alpha^{x}_{i,t}||_{2}$ (45) s.t.
$\displaystyle\mathcal{X}_{i,t}(\alpha^{x^{p}}_{i,t})\subseteq X_{i,t},$

which projects $\alpha_{i,t}^{x}$ to $\alpha^{x^{p}}_{i,t}$ for parameters
corresponding to the state space of subsystem $i$ at time step $t$. There is a
similar QP for the parameters corresponding to the control space for each
subsystem as well. The algorithm is shown in Fig. 2 as a flowchart.

$\alpha$ initializationprojection $\forall i,\forall t$
(45)sys1$\rightarrow$(39)sys2$\rightarrow$(39)sys$\eta\rightarrow$(39)$+$Gradient
Descent (44)return $\alpha,\Omega,\Theta$$\alpha^{p}$$\alpha^{p}$$\alpha^{p}$.
. .$\nu_{1}(\alpha)$ $\Omega_{1}$ $\Theta_{1}$
$\nabla_{\alpha}\nu_{1}$$\nu_{2}(\alpha)$ $\Omega_{2}$ $\Theta_{2}$
$\nabla_{\alpha}\nu_{2}$$\nu_{\eta}(\alpha)$ $\Omega_{\eta}$ $\Theta_{\eta}$
$\nabla_{\alpha}\nu_{\eta}$ $\sum_{i}\nabla_{\alpha}\nu_{i}$$\alpha$if
$\sum_{i}\nu_{i}=0$ Figure 2: Flowchart representation of the algorithm from
Sec. 6.2. The parameters are initialized and then projected to the set of
valid parameters. Then each subsystem independently determines the AG
contracts and its compartment of the potential function. One of the following
three scenarios occurs once the potential function has been determined: (i) If
it is zero, the result is returned; (ii) If it hits a plateau, $k$ is increase
by one unit and the algorithm is restarted with the most recently updated
parameters (this step is not shown); (iii) If (i) or (ii) do not apply, the
gradient information is used to update the parameters and the the algorithm
reiterates.

##  7 Distributed Robust MPC With Stability Guarantee

In this section, we show that our solution for Problem 1 can be used for the
problem of Distributed Robust MPC when solved recursively. Each MPC needs to
have: (i) a given initial state $x_{i,0},\forall i\in\mathcal{I}$, (ii) a
quadratic cost function
$J_{i}:=\sum_{t}{\mathrm{x}^{i}_{t}}^{T}Q\mathrm{x}^{i}_{t}+{\mathrm{u}^{i}_{t}}^{T}R\mathrm{u}^{i}_{t}$
for each subsystem $i$, and (ii) the RCI sets $(\Omega_{i})$ as terminal sets.
In [10, 39], it is shown that considering a terminal condition for an MPC
ensures the stability and recursive feasibility of the approach. The terminal
condition specifies a target set for the MPC’s last time step, which is
encoded using set containment constraints similar to those described in the
preceding sections. Using the provided solution for problem 2, the RCI sets
can be determined in offline mode. This section incorporates these adjustments
into the solution of Problem 1 presented above.

Parameters: The potential function is defined in the same way as before, but
the optimization in (39) and the parameter set are modified. Since the state
of each subsystem changes in the corresponding state space, the centers of the
baseline sets $\bar{c}^{x}_{i,t}$ and $\bar{c}^{u}_{i,t}$ in the parametric AG
contracts can no longer be set in advance, and are included in the set of
parameters. The extended set of parameters is denoted by
$\alpha^{ext}:=\\{\alpha^{x}_{i},\alpha^{u}_{i},\bar{c}_{i}^{x},\bar{c}_{i}^{u}\\}_{i\in\mathcal{I}}$,
where
$\alpha^{x}_{i}=\alpha^{x}_{i,1},...,\alpha^{x}_{i,h-1},\alpha^{u}_{i}=\alpha^{u}_{i,1},...,\alpha^{u}_{i,h-1},\bar{c}_{i}^{x}=\bar{c}_{i,1}^{x},...,\bar{c}_{i,h-1}^{x},\bar{c}_{i}^{u}=\bar{c}_{i,1}^{u},...,\bar{c}_{i,h-1}^{u}$.
The convexity property of the potential function is not affected by this
change in the parameter set. However, the initial time step has no parameters
anymore because the initial states are known. Therefore, there is no need to
construct parameterized sets for it (i.e., $T^{i}_{0}=0_{n_{i}}$ and
$M^{i}_{0}=0_{m_{i}}$). Also, since the goal sets, which are the terminal
sets, are already known, there is no need to specify any parameterized sets
for the last time step. Finally, the potential function is reconstructed as
the solution to the following optimization problem:

$\displaystyle\mathcal{V}_{i}(\alpha^{ext})=\underset{\bar{\mathrm{x}}^{i}_{t},T^{i}_{t},\bar{\mathrm{u}}^{i}_{t},M^{i}_{t},d_{t}^{x},d_{t}^{u}}{\text{min}}\begin{aligned}
&\sum_{t\in\mathbb{N}_{1,h}}d^{x}_{t}+\sum_{t\in\mathbb{N}_{1,h-1}}d^{u}_{t}+\omega
J_{i}\end{aligned}$ subject to
$\displaystyle[0_{k-n_{i}},D_{i,0}^{aug}]=T^{i}_{1},$ (46a)
$\displaystyle\mathcal{Z}(\bar{d}_{i,0}^{aug},D_{i,0}^{aug})=\operatorname{{\text{Red}}}(\sum_{j\neq
i}A_{ij}\bar{\mathrm{x}}^{i}_{0}+B_{ij}\bar{\mathrm{u}}^{i}_{0}\oplus D_{i}),$
(46b)
$\displaystyle[A_{ii}T^{i}_{t}+B_{ii}M^{i}_{t},D_{i,t}^{aug}]=[T^{i}_{t+1}],\hskip
5.69054pt\forall t\in\mathbb{N}_{1,h-1}$ (46c) $\displaystyle
A_{ii}\bar{\mathrm{x}}^{i}_{t}+B_{ii}\bar{\mathrm{u}}^{i}_{t}+\bar{d}_{i,t}^{aug}=\bar{\mathrm{x}}^{i}_{t+1},\hskip
5.69054pt\forall t\in\mathbb{N}_{h-1}$ (46d)
$\displaystyle\mathcal{Z}(\bar{d}_{i,t}^{aug},D_{i,t}^{aug})=\operatorname{{\text{Red}}}(\bigoplus_{j\neq
i}A_{ij}\mathcal{X}_{j,t}(\bar{c}_{i,t}^{x},\alpha^{x}_{j,t})$
$\displaystyle\hskip 22.76228pt\oplus\bigoplus_{j\neq
i}B_{ij}\mathcal{U}_{j,t}(\bar{c}_{i,t}^{u},\alpha^{u}_{j,t})\oplus
D_{i}),\forall t\in\mathbb{N}_{1,h-1}$ (46e)
$\displaystyle\mathcal{Z}(\bar{\mathrm{x}}^{i}_{t},T^{i}_{t})\subseteq\
\mathcal{X}_{i,t}(\bar{c}_{i,t}^{x},\alpha^{x}_{i,t})\oplus\mathcal{Z}(0,d^{x}_{t}I_{n_{i}}),\forall
t\in\mathbb{N}_{1,h-1}$ (46f)
$\displaystyle\mathcal{Z}(\bar{\mathrm{x}}^{i}_{h},T^{i}_{h})\subseteq\
\Omega_{i}\oplus\mathcal{Z}(0,d^{x}_{h}I_{n_{i}}),\hskip 5.69054pt$ (46g)
$\displaystyle\mathcal{Z}(\bar{\mathrm{u}}^{i}_{t},M^{i}_{t})\subseteq\
\mathcal{U}_{i,t}(\bar{c}_{i,t}^{u},\alpha^{u}_{i,t})\oplus\mathcal{Z}(0,d^{u}_{t}I_{m_{i}}),\forall
t\in\mathbb{N}_{1,h-1}$ (46h) $\displaystyle\bar{\mathrm{x}}^{i}_{0}=x_{i,0},$
(46i) $\displaystyle d_{t}^{x},d_{t}^{u}\geq 0,\hskip 5.69054pt\forall
t\in\mathbb{N}_{1,h},$ (46j)

where constraints (46a) and (46b) are originated from Theorem 1, while the
generators are assigned to zero matrices, where $k$ is the number of columns
in $T_{1}^{i}$. Also, (46i) specifies the current state and (46g) imposes the
terminal condition when $d_{h}^{x}=0$, which happens when $\mathcal{V}^{*}=0$.

Objective function: In our formulation, we capture the goal sets using the set
inclusion constraints. However, one might be interested in the optimality of
the MPC formulation, instead of only a feasible solution to a goal set. In
this case, we wish to minimize the potential function as well as a user-
defined cost function $J_{i}$. The main challenge is to satisfy all the set
containment requirements while minimizing $J_{i}$. A very little protrusion of
one of the sets, on the other hand, may cause a cascade of instability among
the subsystems. The objective function is the weighted sum shown in (46),
where $\omega$ is a scalar weight. Given the user-defined cost function, we
start the iterations with $\omega$ set to zero. Following the discovery of
correct contracts, where the objective functions with zero weight for all
subsystems are zero, the weight $\omega$ is increased as long as
$\sum_{t\in\mathbb{N}_{1,h}}d_{t}^{x^{*}}+\sum_{t\in\mathbb{N}_{1,h-1}}d_{t}^{u^{*}}$
stays zero for all subsystems (i.e. the contracts remain correct), where
$d_{t}^{x^{*}}$ and $d_{t}^{u^{*}}$ are optimal values for $d_{t}^{x}$ and
$d_{t}^{u}$, respectively. Since the weight is increased to the point that one
of the subsystems disregards correctness, the number of iterations is limited.

Communication: The method proposed here can be seen as an _iterative_
distributed MPC [11] with disturbance rejection, since neighboring subsystems
share information, including parameters and generated gradients, as many times
as needed until convergence, at each time step. Following convergence, each
subsystem implements the optimal value calculated for
$\bar{\mathrm{u}}^{i}_{0}$, and the cycle continues. Note that after each
convergence, the decentralized output controllers give a guaranteed solution
to the goal sets, making our method resilient to connection loss issues. By
repeatedly solving the problem, on the other hand, the most up-to-date
information is used, which leads to a more optimal solution. In addition,
because we may reuse the optimal variables from the previous iteration, there
is no need to randomly initialize the parameters each time the MPC is solved.
The warm start speeds up the convergence process.

##  8 Case Studies

This section contains four case studies. The first is a simple example of
finding decentralized RCI sets for a connected LTI system with coupled state
constraints. The second example shows the viable sets of a connected network
of LTV systems. The goal is to demonstrate sequences of decentralized time-
limited viable sets. The third illustrates the scalability of the proposed
approach by comparing its computation time with those obtained using a series
of benchmark methods. The final example demonstrates the capability of our
proposed MPC approach in a power network. The source code is available on
GitHub111https://github.com/Kasraghasemi/parsi. A MacBook Pro 2.6 GHz Intel
Core i7 with Gurobi [18] as the optimizer was used for the implementations.

### Case Study 1

Consider an LTI system with

$A=\left[\begin{array}[]{c c| c c |c c}0.1&0.1&0.1&0.02&0.04&-0.02\\\
0&0.1&-0.1&0.06&0&-0.04\\\ \hline\cr-0.08&0&0.1&0.1&0.04&0.1\\\
0.02&-0.06&0&0.1&0.08&0\\\ \hline\cr 0&0&0.04&0.02&0.1&0.1\\\
0&0&0.02&0.1&0&0.1\end{array}\right]$

as the aggregate matrix. We assume that the dynamics of the subsystems are
only coupled through states, hence $B$ is a block diagonal matrix. The state
constraints are expressed in the following coupled form:

$X=\mathcal{Z}(0_{6},\left[\begin{array}[]{c c| c c |c c}1&0.1&0.2&0&0&0\\\
0.1&1&0.02&0&0&0.1\\\ \hline\cr 0&0.01&1&0&0.1&0.1\\\ 0.2&0&0.03&1&0&0.1\\\
\hline\cr 0&0.1&0.1&0&1&0.2\\\ -0.1&-0.02&0.1&0&0&1\end{array}\right]).$

The subsystems are defined by (12) with $B_{ii}=\begin{bmatrix}0\\\
0.1\end{bmatrix},U_{i}=\mathcal{Z}(0_{1},[1]),D_{i}=\mathcal{Z}(0_{2},\begin{bmatrix}0.3&0\\\
0&0.3\end{bmatrix})$. Our goal is to identify decentralized RCI sets
$\Omega_{1},\Omega_{2},\Omega_{3}$ by breaking the 6-dimensional system into
three 2-dimensional subsystems. The first step is to decouple the state
admissible set $X$ into $X_{1},X_{2}$, and $X_{3}$. We accomplish this by
using the semi-definite program in (52). The outcomes are illustrated in Fig
3. The decentralized RCI sets and the correct set of parameters are determined
using the compositional method provided in Sec. 6.2. Fig 4 shows the
projection of potential function and the correct set of parameters on
$\alpha^{x}_{1}[1]-\alpha^{x}_{1}[2]$ space, which are two of the parameters.

Figure 3: Case Study 1: The projections of the 6-dimensional admissible sets
to the corresponding plane for each subsystem are shown in green. The
decomposed sets are depicted in red. As expected, all the red sets are subsets
of the green sets.

Figure 4: Case Study 1: [Left] The green polytope shows the projection of the
correct set of parameters onto the plane of two of the parameters. The blue
trajectory represents the path obtained using the compositional method to get
to the final correct parameters, starting from random parameters. The three
arrows show three preferred directions ($\nabla_{\alpha}\nu_{i}$) for each of
the three subsystems at the start point. [Right] The potential function as a
function of the same two parameters, with all other parameters set at their
final correct values. Its level set at zero corresponds to the correct set of
parameters. The potential function is a piece-wise affine function.

### Case Study 2

Our goal in this case study is to determine the sequences of decentralized
limited-time viable sets for a 8-dimensional LTV system with

$0.01\left[\begin{array}[]{c c| c c |c c | c c}t&100&0.1&0.1&0&0&-\sin(t)&t\\\ 0&100&0.1&0.1&0&0&\log(t+1)&\cos(t)\\\ \hline\cr 0.1&0.1&100&100&0.1&0.1&0&0\\\ 0.1&0.1&0&100&0.1&0.1&0&-1\\\ \hline\cr 0&0&0.1&0.1&100&100&0&1\\\ 0.1t^{2}&0&0.1&0.1&0&100&0&0\\\ \hline\cr t&0&0&0&0&0&100&100\\\ -t&1&0&1&0&0&0&t\end{array}\right],$

as the aggregate time-variant $A_{t}$ matrix. There are four subsystems in the
form of (12) that are only coupled through states in their dynamics. We assume
$B_{ii,t}=[0,0.1]^{T},D_{i,t}=\mathcal{Z}(0_{2},\operatorname{{\text{Diag}}}([0.4,0.4])$
,$X_{1,t}=\mathcal{Z}(0_{2},\operatorname{{\text{Diag}}}([5-\sin\dfrac{\pi
t}{15},6-5.5\sin\dfrac{\pi t}{12}])$
,$X_{2,t}=\mathcal{Z}(0_{2},\operatorname{{\text{Diag}}}([5-2\sin\dfrac{\pi
t}{8},6-5.5\sin\dfrac{\pi t}{20}])$
,$X_{3,t}=\mathcal{Z}(0_{2},\operatorname{{\text{Diag}}}([5-\cos\dfrac{\pi
t}{15},6-5.5\cos\dfrac{\pi t}{12}])$
,$X_{4,t}=\mathcal{Z}(0_{2},\operatorname{{\text{Diag}}}([5-\dfrac{t}{5},5-\dfrac{t}{5}]),U_{i}=\mathcal{Z}(0,[10])$
The horizon is $15$ ($t\in\mathbb{N}_{15}$). The results, which are displayed
in Fig 5, were found using the centralized method described in Sec. 6.1. The
sizes of the viable sets tend to increase in time due to the additive
disturbances. Furthermore, since the objective function is designed to reduce
the area of the viable sets, the viable sets at the first time step are
points, and there are no sets for the first time step in Fig. 5.

Figure 5: The decentralized viable sets for the four coupled LTV subsystems
from Case Study 2. To avoid sets from overlapping in each subsystem, each set
is shifted to the right of the preceding set. In fact, they all have a common
center. The gray areas are state constraints $X_{i,t}$, which correctly
contain the viable sets entirely.

### Case Study 3

This example is adopted from [28], which creates a vast random network of
connected LTI subsystems by randomly spreading a predetermined number of
points in a $100$ by $100$ square area. Each point represents a subsystem with
$A_{ii}=\begin{bmatrix}1&0.2\\\ 0&1\end{bmatrix},B_{ii}=\begin{bmatrix}0\\\
0.2\end{bmatrix},X_{i}=\mathcal{Z}(0,5I_{2}),U_{i}=\mathcal{Z}(0,5I_{1}),D_{i}=\mathcal{Z}(0,0.1I_{2}).$
If the Euclidean distance $\operatorname{{\text{dist}}}(i,j)$ between two
points representing two subsystems $i$ and $j$ is less than $10$, they are
dynamically coupled with the following
$A_{ij}=\dfrac{\lambda}{1+\operatorname{{\text{dist}}}(i,j)}\begin{bmatrix}1&1\\\
1&1\end{bmatrix}$, where $\lambda$ is a constant scalar that limits the effect
of couplings when the number of subsystems is large. We compare three
approaches to obtain the RCI sets: (i) Cen uses Theorem 2 to construct a
centralized RCI set as an output of centralized optimization problem, (ii)
DecCen uses the centralized optimization presented in Sec. 6.1, with a set of
decentralized RCI sets as the result, and (iii) Compose uses the compositional
method described in Sec. 6.2, which likewise produces decentralized RCI sets.
Note that in the LTI case, we need to apply the approach from Remark 3 to make
all constraints linear.

Table 1 shows the average synthesis times for the three approaches and the
average number of iterations for the compositional method for 10 sample runs
for systems with increasing state dimensions $n$. The reported times include
only the time spent solving the optimization problems, not the time spent
constructing them. Approach (i) failed to handle large dimensions and timed
out. Given the polynomial computational complexity of linear programs, this
was expected. While having a larger number of variables, approach (ii)
provides better results, which is mainly due to the sparsity in the control
matrices. The sparsity originates from the fact that we are considering
separable viable and action sets for different subsystems. The second strategy
is the quickest of the three in smaller dimensions. Approach (iii), on the
other hand, demonstrates its scalability by expanding to massive dimensions.
The synthesis time could be reduced further by utilizing parallel computation
techniques using multi-core computers, which we will address in future work.

Table 1: Average synthesis times (in seconds) for the three approaches from
Case Study 3 and the average number of iterations for the compositional method

$n$ | $\lambda$ | Cen | DecCen | Compose | Avg. Iter.
---|---|---|---|---|---
10 | 0.1 | 0.26 | 0.01 | 0.03 | 3.4
20 | 0.1 | 3.23 | 0.10 | 0.05 | 38.3
30 | 0.1 | 19.60 | 0.26 | 0.24 | 99.0
50 | 0.01 | 193.69 | 1.03 | 0.56 | 24.8
100 | 0.01 | — | 4.34 | 1.66 | 33.8
200 | 0.01 | — | 37.03 | 6.46 | 73.6
500 | 0.01 | — | 576.72 | 10.63 | 47.0
1000 | 0.001 | — | — | 12.40 | 30.3
10000 | 0.0001 | — | — | 153.00 | 30.0

### Case Study 4

The load-frequency control (LFC) problem in power networks [9] is addressed in
this example. A power network is a collection of areas, where each area has
its own generator and users, and can transfer its excess power to other areas
through a set of links. The phase angle $\delta_{i,t}$ and frequency $f_{i,t}$
comprise the state of area/subsystem $i$ at time $t$. The control input is the
amount of power $u_{i,t}$ produced by the generator. The goals are to keep the
system in a state that is close to its nominal value and to keep the power
exchanges between areas on schedule. Since only the departure from the nominal
state matters, the state vector for subsystem $i$ is written as
$x_{i,t}=[\Delta\delta_{i,t},\Delta f_{i,t}]^{T}$, with $\Delta$ denoting the
divergence from the nominal value. The system dynamics from [9] is discretized
using the Euler method, resulting in $A_{ii}=\begin{bmatrix}1,&2\pi\Delta t\\\
\dfrac{-\Delta tk_{p_{i}}}{2\pi
T_{p_{i}}}\displaystyle\sum_{j\in\mathcal{N}_{i}}K_{s_{ij}},&1-\dfrac{\Delta
t}{T_{p_{i}}}\end{bmatrix},A_{ij}=\begin{bmatrix}0,&0\\\ \dfrac{\Delta
tk_{p_{i}}K_{s_{ij}}}{2\pi
T_{p_{i}}},&0\end{bmatrix},B_{ii}=\begin{bmatrix}0\\\ \dfrac{K_{p_{i}}\Delta
t}{T_{p_{i}}}\end{bmatrix}$, where $\mathcal{N}_{i}$ is a set containing the
neighbours of subsystem $i$ and $\Delta t,K_{p_{i}},K_{s_{ij}},T_{p_{i}}$ are
the discretized time interval, system gain, synchronizing coefficient between
area $i$ and $j$, and system model time constant, respectively, which are set
to 0.1s, 110, 0.5, and 25s for all areas. Consider a network of four fully
connected areas with bidirectional connections between each pair, with the
following bounds over control input and load disturbance for each subsystem:
$|u_{i,t}|\leq 0.1$ and
$D_{i}=\mathcal{Z}(0_{2},\operatorname{{\text{Diag}}}([\epsilon,\dfrac{-\Delta
tK_{p_{i}}\Delta P_{d_{i}}}{T_{p_{i}}}]))$, where $\epsilon=10^{-10}$ and
$\Delta P_{d_{i}}=0.01p.u$ is the load disturbance for area $i$. The initial
state for the aggregate system is
$[-0.2,0.1,0.02,0.01,0.1,-0.05,-0.03,-0.01]^{T}$ and the goal is to go to the
set where $|\Delta\delta_{i,h}|\leq 0.01,|f_{i,h}|\leq 0.01,\forall i$. Using
the solution to Problem 2, we first find a set of decentralized RCI sets
inside the goal set. Then, we proceed with the strategy from Sec. 7 with
horizon $h=5$. Also, in the cost function, the $R$ and $Q$ matrices are set to
zero and the identity matrix, respectively. The output viable sets after
solving the problem are illustrated in green for each subsystem separately in
Fig. 6, where as expected, the final viable set in each subsystem ends up
inside its goal set. If the goal is just to reach the goal set, there is
actually no need to solve the problem again in the next time step for the new
initial states. Each subsystem can implement the controller in the derived
sequence, which by construction. has been shown it maintains the states inside
the derived viable sets.

Figure 6: The goal sets and sequences of viable sets for each area of the
power network from Case Study 4 are shown in blue and green, respectively. The
red dashed lines begin at the initial state and link the centers of the viable
sets, i.e., the nominal predicted trajectory.

##  9 Conclusion and Future Work

We proposed a convex parameterization of assume-guarantee contracts that
facilitates compositional control synthesis of large-scale linear systems with
uncertainties. We achieved linear time complexity, while other state-of-the-
art algorithms offer polynomial time complexity at best. Extensions to
nonlinear and (partially) unknown dynamics, investigating probabilistic
assume-guarantee contracts, and implementations in robotic experimental
platforms are directions of future research.

This work was partially supported by the National Science Foundation under
grants IIS-2024606 and IIS-1723995.

## References

  * [1] James Anderson, John C Doyle, Steven H Low, and Nikolai Matni.  System level synthesis.  Annual Reviews in Control, 47:364–393, 2019.
  * [2] Gianluca Antonelli.  Interconnected dynamic systems: An overview on distributed control.  IEEE Control Systems Magazine, 33(1):76–88, 2013.
  * [3] Lubomir Bakule.  Decentralized control: An overview.  Annual reviews in control, 32(1):87–98, 2008.
  * [4] Dimitris Bertsimas and John N Tsitsiklis.  Introduction to linear optimization, volume 6.  Athena Scientific Belmont, MA, 1997.
  * [5] F. Blanchini.  Set invariance in control.  Automatica, 35(11):1747 – 1767, 1999.
  * [6] Sergiy Bogomolov, Goran Frehse, Marius Greitschus, Radu Grosu, Corina Pasareanu, Andreas Podelski, and Thomas Strump.  Assume-guarantee abstraction refinement meets hybrid systems.  In Haifa verification conference, pages 116–131. Springer, 2014\.
  * [7] Stephen Boyd, Neal Parikh, and Eric Chu.  Distributed optimization and statistical learning via the alternating direction method of multipliers.  Now Publishers Inc, 2011.
  * [8] C. Combastel.  A state bounding observer based on zonotopes. in proc. of the european control conference.  page 2589–2594, 2003.
  * [9] Eduardo F. Camacho and Carlos Bordons.  Distributed model predictive control.  Optimal Control Applications and Methods, 36(3):269–271, 2015.
  * [10] Luigi Chisci, A Lombardi, and Edoardo Mosca.  Dual-receding horizon control of constrained discrete time systems.  European Journal of Control, 4(2):278–285, 1996.
  * [11] Panagiotis D Christofides, Riccardo Scattolini, David Munoz de la Pena, and Jinfeng Liu.  Distributed model predictive control: A tutorial review and future research directions.  Computers & Chemical Engineering, 51:21–41, 2013.
  * [12] Christian Conte, Colin N Jones, Manfred Morari, and Melanie N Zeilinger.  Distributed synthesis and stability of cooperative distributed model predictive control for linear systems.  Automatica, 69:117–125, 2016.
  * [13] Georgios Darivianakis, Angelos Georghiou, and John Lygeros.  Decentralized decision making for networks of uncertain systems.  arXiv preprint arXiv:1803.07660, 2018.
  * [14] Makan Fardad, Fu Lin, and Mihailo R Jovanović.  Sparsity-promoting optimal control for a class of distributed systems.  In Proceedings of the 2011 American Control Conference, pages 2050–2055. IEEE, 2011.
  * [15] Kasra Ghasemi, Sadra Sadraddini, and Calin Belta.  Compositional synthesis via a convex parameterization of assume-guarantee contracts.  In Proceedings of the 23rd International Conference on Hybrid Systems: Computation and Control, HSCC ’20, New York, NY, USA, 2020. Association for Computing Machinery.
  * [16] Paul J Goulart, Eric C Kerrigan, and Jan M Maciejowski.  Optimization over state feedback policies for robust control with constraints.  Automatica, 42(4):523–533, 2006.
  * [17] Eugene Gover and Nishan Krikorian.  Determinants and the volumes of parallelotopes and zonotopes.  Linear Algebra and its Applications, 433(1):28–40, jul 2010.
  * [18] LLC Gurobi Optimization.  Gurobi optimizer reference manual, 2020.
  * [19] Thomas A Henzinger, Marius Minea, and Vinayak Prabhu.  Assume-guarantee reasoning for hierarchical hybrid systems.  In International Workshop on Hybrid Systems: Computation and Control, pages 275–290. Springer, 2001.
  * [20] Thomas A Henzinger, Shaz Qadeer, and Sriram K Rajamani.  You assume, we guarantee: Methodology and case studies.  In International Conference on Computer Aided Verification, pages 440–451. Springer, 1998.
  * [21] W Kiihn.  Rigorously Computed Orbits of Dynamical Systems without the Wrapping Effect.  Technical report, 1998.
  * [22] Eric S Kim, Murat Arcak, and Sanjit A Seshia.  A Small Gain Theorem for Parametric Assume-Guarantee Contracts *.
  * [23] Eric S Kim, Murat Arcak, and Sanjit A Seshia.  Compositional controller synthesis for vehicular traffic networks.  In 2015 54th IEEE Conference on Decision and Control (CDC), pages 6165–6171. IEEE, 2015.
  * [24] Ilya Kolmanovsky and Elmer G Gilbert.  Theory and computation of disturbance invariant sets for discrete-time linear systems.  Mathematical problems in engineering, 4, 1998.
  * [25] Anna Kathrin Kopetzki, Bastian Schurmann, and Matthias Althoff.  Methods for order reduction of zonotopes.  2017 IEEE 56th Annual Conference on Decision and Control, CDC 2017, 2018-Janua(Cdc):5626–5633, 2018.
  * [26] Fu Lin, Makan Fardad, and Mihailo R Jovanović.  Design of optimal sparse feedback gains via the alternating direction method of multipliers.  IEEE Transactions on Automatic Control, 58(9):2426–2431, 2013.
  * [27] Weixuan Lin and Eilyan Bitar.  Design of Robust Decentralized Controllers via Assume-Guarantee Contracts.  dec 2020.
  * [28] N. Motee and A. Jadbabaie.  Optimal control of spatially distributed systems.  IEEE Transactions on Automatic Control, 53(7):1616–1629, 2008.
  * [29] Petter Nilsson and Necmiye Ozay.  Synthesis of separable controlled invariant sets for modular local control design.  In 2016 American Control Conference (ACC), pages 5656–5663. IEEE, 2016.
  * [30] Pierluigi Nuzzo, Huan Xu, Necmiye Ozay, John B Finn, Alberto L Sangiovanni-Vincentelli, Richard M Murray, Alexandre Donzé, and Sanjit A Seshia.  A contract-based methodology for aircraft electric power system design.  IEEE Access, 2:1–25, 2013.
  * [31] B. O’Donoghue, E. Chu, N. Parikh, and S. Boyd.  Conic optimization via operator splitting and homogeneous self-dual embedding.  Journal of Optimization Theory and Applications, 169(3):1042–1068, June 2016.
  * [32] S. V. Raković, E. C. Kerrigan, D. Q. Mayne, and K. I. Kouramas.  Optimized robust control invariance for linear discrete-time systems: Theoretical foundations.  Automatica, 43(5):831–841, 2007.
  * [33] Ralph Tyrell Rockafellar.  Convex analysis.  Princeton university press, 2015.
  * [34] Michael Rotkowitz and Sanjay Lall.  A characterization of convex problems in decentralized control.  IEEE transactions on Automatic Control, 50(12):1984–1996, 2005\.
  * [35] Sadra Sadraddini, János Rudan, and Calin Belta.  Formal synthesis of distributed optimal traffic control policies.  In Proceedings of the 8th International Conference on Cyber-Physical Systems, pages 15–24, 2017.
  * [36] Sadra Sadraddini and Russ Tedrake.  Linear encodings for polytope containment problems.  arXiv preprint arXiv:1903.05214, 2019.
  * [37] Adnane Saoud, Antoine Girard, and Laurent Fribourg.  Contract-based design of symbolic controllers for safety in distributed multiperiodic sampled-data systems.  IEEE Transactions on Automatic Control, 2020.
  * [38] Sheila Scialanga and Konstantinos Ampountolas.  Robust constrained interpolating control of interconnected systems.  In 2018 IEEE Conference on Decision and Control (CDC), pages 7016–7021. IEEE, 2018.
  * [39] Pierre OM Scokaert, David Q Mayne, and James B Rawlings.  Suboptimal model predictive control (feasibility implies stability).  IEEE Transactions on Automatic Control, 44(3):648–654, 1999.
  * [40] Tyler H Summers and John Lygeros.  Distributed model predictive consensus via the alternating direction method of multipliers.  In 2012 50th annual Allerton conference on communication, control, and computing (Allerton), pages 79–84. IEEE, 2012.
  * [41] Russ Tedrake and the Drake Development Team.  Drake: Model-based design and verification for robotics, 2019.
  * [42] Paul Trodden and Arthur Richards.  Cooperative distributed mpc of linear systems with coupled constraints.  Automatica, 49(2):479–487, 2013.
  * [43] Xuejiao Yang and Joseph K. Scott.  A comparison of zonotope order reduction techniques.  Automatica, 2018.

## Appendix

Proof of Lemma 2 Using constraints (6) for
$\mathcal{Z}(\bar{c}_{1},G_{1})\subseteq\mathcal{Z}(\bar{c}_{2},G_{2}\operatorname{{\text{Diag}}}(\alpha))$,
we have:

$\begin{array}[]{c}G_{1}=G_{2}\operatorname{{\text{Diag}}}(\alpha)\Gamma,\\\
c_{2}-c_{1}=G_{2}\operatorname{{\text{Diag}}}(\alpha)\gamma,||[\Gamma,\gamma]||_{\infty}\leq
1.\end{array}$ (47)

We can replace $\operatorname{{\text{Diag}}}(\alpha)\Gamma$ and
$\operatorname{{\text{Diag}}}(\alpha)\gamma$ with $\Gamma^{new}$ and
$\gamma^{new}$, respectively. We have:

$G_{1}=G_{2}\Gamma^{new},$ (48a) $c_{2}-c_{1}=G_{2}\gamma^{new},$ (48b)
$||[\operatorname{{\text{Diag}}}(\alpha^{-1}),\Gamma^{new},\operatorname{{\text{Diag}}}(\alpha^{-1})\gamma^{new}]||_{\infty}\leq
1,$ (48c)

where $\alpha^{-1}$ is element-wise. In (48c), each row of matrix
$[\Gamma^{new},\gamma^{new}]$ is divided by the corresponding element in
vector $\alpha$. Because all the elements of $\alpha$ are positive, we can
multiply the inequality by $\operatorname{{\text{Diag}}}(\alpha)$ and have:

$[|\Gamma^{new}|,|\gamma^{new}|]\mathbb{1}_{s}\leq\alpha$ (49)

## Solution to Subproblem 1

By definition, we have $\\{[x_{1}^{T},x_{2}^{T},...,x_{\eta}^{T}]^{T}|\forall
x_{i}\in X_{i},i\in\mathcal{I}\\}=\prod_{i\in\mathcal{I}}X_{i}$. Thus, it can
be seen that for all $x_{i}\in
X_{i},[x_{1}^{T},x_{2}^{T},...,x_{\eta}^{T}]^{T}\in X$ holds, if

$\prod_{i\in\mathcal{I}}X_{i}\subseteq X,$ (50)

where $X_{i}=\mathcal{Z}(c_{i},G_{i})$ and

$\prod_{i\in\mathcal{I}}X_{i}=\mathcal{Z}([c_{1}^{T},\cdots,c_{\eta}^{T}]^{T},\text{BlockDiag}([G^{x}_{1},\cdots,G^{x}_{\eta}])).$

Since $c_{i}$s and $G_{i}^{x}$s are unknown. we encode (50) into a linear
program with the help of the linear containment encoding proposed in [36].
However, we are interested in finding the maximum volume for
$\prod_{i\in\mathcal{I}}X_{i}$.

######  Assumption 1.

Each $G^{x}_{i},i\in\mathcal{I}$ is a symmetric square positive semi-definite
matrix.

The assumption that the matrices are square makes the calculation of the
volume of $\prod_{i\in\mathcal{I}}X_{i}$ easier by the formulation in (3).
However, the function in (3) is still non-concave, which is not desirable for
an optimization. It can be shown that by further assuming that all $G_{i}^{x}$
are symmetric positive semi-definite $((G^{x}_{i})^{T}=G^{x}_{i}\text{ and
}G^{x}_{i}\succeq 0)$, the logarithm of the volume in (3), which is equal to

$\text{log}(\text{Volume}(\prod_{i\in\mathcal{I}}X_{i}))=\\\
n\text{log}(2)+2\text{log}(\text{det}(\text{BlockDiag}([G^{x}_{1},\cdots,G^{x}_{\eta}])))$
(51)

becomes concave because $log(det(G))$ when $G$ is a positive semi-definite
matrix is a concave function. Considering the assumptions and removing the
constants in the objective function, the final optimization problem for
decomposing $X$ is as follows:

$\displaystyle\min_{c_{i},G^{x}_{i}}\quad$
$\displaystyle\text{log}(\text{det}(\text{Blk}([G^{x}_{1},\cdots,G^{x}_{\eta}])))$
(52a) s.t.
$\displaystyle\prod_{i\in\mathcal{I}}{\mathcal{Z}(c_{i},G^{x}_{i})}\subseteq
X,\hskip 5.69054pt$ (52b)
$\displaystyle\text{Blk}(G_{1}^{x},...,G_{\eta}^{x})\succeq 0,$ (52c)
$\displaystyle G_{i}^{x}\text{ is a symmetric square matrix},\forall
i\in\mathcal{I}$ (52d)

This can be solved with any semi-definite programming solver. In this paper,
we solved it by SCS [31] using drake [41].