# Decentralized Signal Temporal Logic Control for Perturbed Interconnected
Systems via Assume-Guarantee Contract Optimization

Kasra Ghasemi, Sadra Sadraddini, and Calin Belta This work was partially
supported by the NSF under grants IIS-2024606 and IIS-1723995.K. Ghasemi and
C. Belta are with the Division of System Engineering, Boston University,
Boston, MA 02215, USA<EMAIL_ADDRESS><EMAIL_ADDRESS>Sadraddini is with the
Computer Science and Artificial Intelligence Laboratory, Massachusetts
Institute of Technology, Cambridge, MA 02139, USA<EMAIL_ADDRESS>
###### Abstract

We develop a novel decentralized control method for a network of perturbed
linear systems with dynamical couplings subject to Signal Temporal Logic (STL)
specifications. We first transform the STL requirements into set containment
problems and then we develop controllers to solve these problems. Our approach
is based on treating the couplings between subsystems as disturbances, which
are bounded sets that the subsystems negotiate in the form of parametric
assume-guarantee contracts. The set containment requirements and parameterized
contracts are added to the subsystems’ constraints. We introduce a centralized
optimization problem to derive the contracts, reachability tubes, and
decentralized closed-loop control laws. We show that, when the STL formula is
separable with respect to the subsystems, the centralized optimization problem
can be solved in a distributed way, which scales to large systems. We present
formal theoretical guarantees on robustness of STL satisfaction. The
effectiveness of the proposed method is demonstrated via a power network case
study.

##  I INTRODUCTION

Multi agent systems benefit from decentralized control laws that require only
local information. Online computations lessen as each system implements its
own control law. Furthermore, communication issues are also mitigated as
agents do not need to constantly share information. Multi agent control drew a
lot of interest in control theory, particularly in the last two decades [22,
6, 9]. For multi agent systems with hard constraints and uncertainties,
rigorous mathematical tools are required to reason about the closed-loop
performance of the aggregate system.

Formal methods provide mathematical guarantees for the behavior of control
systems. Formal languages, such as temporal logics [1], can be used to
describe system specifications. With particular relevance to this work, Signal
Temporal Logic (STL) [18] can describe a broad range of temporally bounded
constraints. For example, the STL formula
$\psi={\textbf{F}}_{[0,10]}(x>10)\vee{\textbf{G}}_{[5,20]}(x<0)$ reads in
plain English “ _eventually_ between time 0 and 10 the value of $x$ exceeds 10
or _always_ between 5 and 20 the value of $x$ remains below 0”. The use of
formal methods in multi-agent systems has also been investigated [23, 16, 15].
But, with only one exception [15], they only studied dynamically decoupled
agents, and none of them took into account the presence of additive
disturbances. A related approach in formal methods is based on set-valued
dynamics. Analyzing such systems enables characterizing all the possible
responses in the presence of bounded uncertainties. Reachability analysis and
correct-by-design control synthesis, which guarantee correctness without the
need for system testing, received a lot of attention in recent years [11, 17,
8].

Formal methods come with a high computational cost, which makes it challenging
to apply them to multi agent systems. That is especially true when we are
considering systems with disturbances, and want to guarantee the satisfaction
of temporal logic specification under all allowed disturbances. Divide and
conquer techniques are a natural way to break the problem into smaller pieces.
They can be applied to interconnected systems, where the dynamics of the
agents are coupled. Assume-guarantee contracts [4] formalize the promises that
systems make and provide over dynamical couplings. For instance, assume-
guarantee contracts were used to describe vehicular flow between neighborhoods
of a traffic network [13], aircraft power distributions, [21], and dynamics of
an aerial robot tethered to a ground one [19].

In this paper, we study the problem of decentralized control design for
interconnected perturbed linear systems subject to STL constraints. Unlike
approaches that assume that feasible assume-guarantee contracts are given
a-priori [20, 5], we parameterize the contracts and search for feasibility.
Unlike the search methods in [13, 14], our parameterization, which is based on
our prior work [10], has a special convexity property that leads to a
tractable solution. The approach in [19] also parameterized contracts and
found them using convex optimization, but was limited to polytopic invariant
sets. Here we include complex, non-convex STL constraints, and retain the
parameterization from [10]. The main contributions of this paper are as
follows:

  1. 1.

By fixing the “logical behavior” through solving a mixed-integer program, we
are able to convert the STL specifications into set containment problems.
Then, a linear program is proposed to jointly optimize assume-guarantee
contracts, set-valued trajectories, and decentralized closed loop control
laws. This allows steering the aggregate system in a way that the global STL
formulae is satisfied, while disturbances are rejected in a decentralized
manner. The resulting bounds are computed using assume-guarantee contracts and
are connected to the STL robustness score, a signed distance to satisfaction.

  2. 2.

When the given STL formula is separable with respect to the subsystems, we
provide a method to make the contribution above computationally more tractable
for large networks by making it compositional. We use the convexity properties
in [10] to optimize contracts, reachability sets, and controllers in a
distributed way.

The rest of the paper is organized as follows. We first provide the notation
and the necessary background in Section II. We state the problem in Section
III. The solution is provided in Sections IV and V. Finally, an illustrative
example is shown in Section VI.

##  II Notations and Preliminaries

###  II-A Notation

$\mathbb{R}$, $\mathbb{R}_{+}$ and $\mathbb{N}$ stand for the sets of real,
non-negative real, and non-negative integers, respectively; $\mathbb{N}_{h}$
represents the set of non-negative numbers up to $h\in\mathbb{N}$. An
$h$-dimensional box is defined as
$\mathbb{B}_{h}:=\\{b\in\mathbb{R}^{h}|||b||_{\infty}\leq 1\\}$.
$\mathbb{S}_{1}\oplus\mathbb{S}_{2}:=\\{s_{1}+s_{2}|s_{1}\in\mathbb{S}_{1},s_{2}\in\mathbb{S}_{2}\\}$
is the Minkowski sum of two sets $\mathbb{S}_{1}$ and $\mathbb{S}_{2}$. The
Directed Hausdorff distance $d_{DH}(\mathbb{S}_{1},\mathbb{S}_{2})$ is a
quantitative measure of how far $\mathbb{S}_{2}$ is from being a subset of
$\mathbb{S}_{1}$, and it can be computed as:

$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 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 sets
$\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}$ by $\prod_{i=1}^{N}\mathbb{S}_{i}$.
$I_{n}$, $0_{n}$, and $[A_{1},A_{2}]$ represent the $n\times n$ identity
matrix, the $n$-dimensional zero vector, and the horizontal concatenation of
matrices $A_{1}$, $A_{2}$ with the same number of rows, respectively.

###  II-B Zonotopes

A zonotope is a symmetric shape set representation defined as
$\mathcal{Z}(c,G):=\\{c+Gb|\forall b\in\mathbb{B}_{q}\\}$, where
$c\in\mathbb{R}^{n}$ and $G\in\mathbb{R}^{n\times q}$ $(n,q\in\mathbb{N})$
denote the zonotope’s center and generator, respectively. The order of the
zonotope is equal to $\dfrac{q}{n}$. Zonotopes are convenient for set
calculations, such as Minkowski sums and linear transformations. Given two
sets $\mathbb{S}_{1}=\mathcal{Z}(c_{1},G_{1})$ and
$\mathbb{S}_{2}=\mathcal{Z}(c_{2},G_{2})$, a matrix $A\in\mathbb{R}^{m\times
n}$, and a vector $b\in\mathbb{R}^{n}$, where $c_{1},c_{2}\in\mathbb{R}^{n}$
and $G_{1}\in\mathbb{R}^{n\times q_{1}}$ , $G_{2}\in\mathbb{R}^{n\times
q_{2}}$, we have
$\mathbb{S}_{1}\oplus\mathbb{S}_{2}=\mathcal{Z}(c_{1}+c_{2},[G_{1},G_{2}])$
and $A\mathbb{S}_{1}+b=\mathcal{Z}(Ac_{1}+b,AG_{1})$.

###  II-C Specifications

Signal Temporal Logic (STL) was introduced in [18] to specify Boolean and
temporal properties of real-valued, time signals. A discrete-time signal is a
function $s:\mathbb{N}\rightarrow\mathbb{R}^{q}$. We use $(s,[t_{1},t_{2}])$
to denote the sequence $s(t_{1}),...,s(t_{2})$ and $(s,t)$ for
$(s,[t,\infty])$. An STL formula is defined with the following recursive
grammar:

$\varphi::=\pi|\neg\varphi|\varphi\wedge\psi|\varphi|\varphi\vee\psi|\textbf{F}_{[t_{1},t_{2}]}\varphi|\textbf{G}_{[t_{1},t_{2}]}\varphi|\varphi\textbf{U}_{[t_{1},t_{2}]}\psi$
(2)

where $\pi$ is a predicate. All predicates are assumed to be linear in the
form $p(s)\leq c$ or $p(s)\geq c$, with $c$ being a scalar and
$p:\mathbb{R}^{q}\rightarrow\mathbb{R}$ being a linear function. Symbols
$\neg$ , $\wedge$ , and $\vee$ denote Boolean negation, conjunction, and
disjunction, respectively; $\textbf{F}_{[t_{1},t_{2}]}$,
$\textbf{G}_{[t_{1},t_{2}]}$, and $\textbf{U}_{[t_{1},t_{2}]}$ are temporal
operators for “eventually”,“always”, and “until”, respectively. Also,
$(s,t)\models\varphi$ denotes that signal $s$ satisfies formula $\varphi$ at
time $t$, and $(s,t)\nvDash\varphi$ if this is not the case.

###### Definition 1

The satisfaction of a formula by a signal $s$ at time $t$ is defined as
follows:

  * •

$(s,t)\models(p(s)\geq c)\Leftrightarrow p(s(t))\geq c$ ,

  * •

$(s,t)\models(p(s)\leq c)\Leftrightarrow p(s(t))\leq c$ ,

  * •

$(s,t)\models\neg\varphi\Leftrightarrow(s,t)\nvDash\varphi$ ,

  * •

$(s,t)\models\varphi_{1}\wedge\varphi_{2}\Leftrightarrow(s,t)\models\varphi_{1}\wedge(s,t)\models\varphi_{2}$
,

  * •

$(s,t)\models\varphi_{1}\vee\varphi_{2}\Leftrightarrow(s,t)\models\varphi_{1}\vee(s,t)\models\varphi_{2}$
,

  * •

$(s,t)\models\textbf{G}_{[t_{1},t_{2}]}\varphi\Leftrightarrow\forall
t^{\prime}\in[t_{1},t_{2}],(s,t)\models\varphi$ ,

  * •

$(s,t)\models\textbf{F}_{[t_{1},t_{2}]}\varphi\Leftrightarrow\exists
t^{\prime}\in[t_{1},t_{2}],(s,t^{\prime})\models\varphi$ ,

  * •

$(s,t)\models\varphi_{1}\textbf{U}_{[t_{1},t_{2}]}\varphi_{2}\Leftrightarrow\exists
t^{\prime}\in[t_{1},t_{2}],(s,t^{\prime})\models\varphi_{2}\wedge\forall
t^{\prime\prime}\in[t_{1},t^{\prime\prime}](s,t^{\prime\prime})\models\varphi_{1}$.

For simplicity, $(s,0)\models\varphi$ is denoted by $s\models\varphi$. The
horizon of a formula is the shortest amount of time required to determine
whether a formula $\varphi$ is satisfied, and it is denoted by $hrz(\varphi)$
[2]. The robustness [7] of an STL formula with respect to a signal determines
how strongly the signal satisfies / violates the formula. Robustness is a real
function that produces a score, where larger scores mean stronger
satisfaction. The robustness of the formula $\varphi$ with respect to the
signal $s$ at time $t$ is denoted by $\rho(s,\varphi,t)$ and can be computed
recursively. The procedure begins with the predicates, where each predicate’s
robustness is defined as $\rho(s,(p(s)\geq c),t)=\rho=p(s)-c$. Without loss of
generality, we only consider negation free formulas in this paper. This is not
restrictive, as any STL formula can be made negation-free. It is also worth
noting that, while predicates with inequalities are used in the semantics
definition, strict inequalities and equalities can be formed using the Boolean
operators.

##  III Problem Definition and Approach

Consider the following network of coupled time-variant linear subsystems:

$x_{i}(t+1)=A_{ii}(t)x_{i}(t)+B_{ii}(t)u_{i}(t)+\sum_{j\neq
i}A_{ij}(t)x_{j}(t)\\\ +\sum_{j\neq
i}B_{ij}(t)u_{j}(t)+w_{i}(t),\;i\in\mathcal{I},$ (3)

where $\mathcal{I}$ is an index set for the subsystems;
$A_{ii}(t)\in\mathbb{R}^{n_{i}\times n_{i}}$,
$A_{ij}(t)\in\mathbb{R}^{n_{i}\times n_{j}}$,
$B_{ii}(t)\in\mathbb{R}^{n_{i}\times m_{i}}$, and
$B_{ij}(t)\in\mathbb{R}^{n_{i}\times m_{j}}$ are given, time-variant matrices
for subsystem $i$. Let $\eta=|\mathcal{I}|$ denote the number of subsystems in
the network. The state, control input, and disturbance for subsystem $i$ at
time step $t$ are represented by $x_{i}(t)\in\mathbb{R}^{n_{i}}$,
$u_{i}(t)\in\mathbb{R}^{m_{i}}$, and $w_{i}(t)\in\mathbb{R}^{n_{i}}$, which
are bounded by given polytopic sets $x_{i}(t)\in
X_{i}(t)\subseteq\mathbb{R}^{n_{i}}$, $u_{i}(t)\in
U_{i}(t)\subseteq\mathbb{R}^{m_{i}}$, and $w_{i}(t)\in
W_{i}(t)\subset\mathbb{R}^{n_{i}}$, respectively. A decentralized controller
$\mu_{i}(.,t):X_{i}(t)\rightarrow U_{i}(t)$ is a function that maps the
current state of subsystem $i$ into a control input in the control space of
the same subsystem. System (3) with no disturbances is called a nominal
system.

######  Definition 2 (Decentralized Finite-Time Viable Sets)

Given $h\in\mathbb{N}$, the sequences of sets
$\Omega_{i}(0),\Omega_{i}(1),...,\Omega_{i}(h)$, $i\in\mathcal{I}$ for the
interconnected system in (3) are called decentralized _viable_ sets, if for
all $t\in\mathbb{N}_{h},\forall i\in\mathcal{I}$, $\Omega_{i}(t)\subseteq
X_{i}(t)$ and there exists a set of policies $\mu_{i}(.,t)$ such that
$\Theta_{i}(t)\subseteq U_{i}(t)$ and $\forall t\in\mathbb{N}_{h-1},\forall
x_{i}(t)\in\Omega_{i}(t),\forall w_{i}(t)\in W_{i}(t)\Rightarrow
x_{i}(t+1)\in\Omega_{i}(t+1)$, where $\Theta_{i}(t):=\mu(\Omega_{i}(t),t)$ is
called action set.

A signal $s:\mathbb{N}\rightarrow X\times U\subset\mathbb{R}^{n+m}$ is a
trajectory where $s(t)$ represents a vector stacking the state and control of
the aggregated system at time step $t$, which is represented by
$s(t)=(x(t),u(t))$, where
$x(t)=[x^{T}_{1}(t),\cdots,x^{T}_{\eta}(t)]^{T}\in\mathbb{R}^{n}$ and
$u(t)=[u^{T}_{1}(t),\cdots,u^{T}_{\eta}(t)]^{T}\in\mathbb{R}^{m}$ and
$n=\sum_{i\in\mathcal{I}}n_{i}$ and $m=\sum_{i\in\mathcal{I}}m_{i}$.

In this paper, we consider the following problem:

###### Problem 1

Given a network of perturbed linear systems in the form (3), the initial
states $x_{i}^{initial}(0)\in X_{i}(0),\forall i\in\mathcal{I}$, a bounded STL
formula $\varphi$ with linear predicates in the states and / or controls, and
a quadratic cost $J:\mathbb{S}\rightarrow\mathbb{R}_{+}$, find the optimal
decentralized controllers $\mu_{i}(x_{i}(t),t),\forall i\in\mathcal{I}$ and
their corresponding sequence of viable sets $\Omega_{i}(t)$ such that $J$ is
minimized, $x_{i}(t)\in X_{i}(t)$, $u_{i}(t)\in U_{i}(t)$, and
$s\models\varphi$. If such a signal does not exist, find $\Omega_{i}(t)$
corresponding to the maximum possible value of the robustness, i.e, find a
signal with the least amount of violation.

To solve Problem 1, a two-step optimization-based approach is proposed. We
begin by solving a mixed-integer program for the aggregated nominal system,
which is constrained by the STL formula $\varphi$ [2]. It allows us to
determine active predicates at each time and convert the STL formula
satisfaction into a set containment problems, which is shown to be a convex
programming problem [25]. In the second step, we take into account the
additive disturbance, along with the set containment constraints, and we find
a set of decentralized closed-loop controllers and viable sets. The technical
details are explained in the next sections.

##  IV Converting STL formulas into Set Containment Problems

In this section, the method from [2] is used to encode an STL formula into a
mixed-integer linear program. Then, the set of predicates whose satisfaction
corresponds to the maximum robustness for the nominal system is identified and
transformed to a set containment problem.

###  IV-A Encoding the STL Formulas

Following the predicate-based encoding from [2], a binary variable
$z_{t}^{\pi}\in\\{0,1\\}$ is dedicated to each predicate $\pi=(y\geq 0)$,
which must be assigned to $1$ if the predicate is true, and to 0 otherwise.
The relation between $z_{t}^{\pi}$, the robustness $\rho$, and $y_{t}$ is
encoded as

$y_{t}+M(1-z_{t}^{\pi})\geq\rho\quad,\quad y_{t}-Mz_{t}^{\pi}<\rho,$ (4)

where $M$ is a sufficiently large number such that for all time steps,
$M\geq\text{max\color[rgb]{1,0,0} }y_{i},i\in\mathbb{N}_{n_{y}}$. The
equations in (4) enforce the binary variable $z_{t}^{\pi}$ to be equal to $1$
when $y_{t}\geq\rho$ and equal to $0$ when $y_{t}<\rho$. Disjunctions and
conjunctions are captured by the following constraints:

$z=\bigwedge_{i=0}^{n_{z}}z_{i}\Rightarrow z\leq
z_{i},i\in\mathbb{N}_{n_{z}},z=\bigvee_{i=0}^{n_{z}}z_{i}\Rightarrow
z\leq\sum_{i=0}^{n_{z}}z_{i},$ (5)

where $n_{z}\in\mathbb{N}$ and $z\in[0,1]$ is declared as a continuous
variable. However, as the above equation shows, it can only take binary
values. In [12], [24] upper-bounding constraints are added to create a
necessary and sufficient condition:

$z=\bigwedge_{i=0}^{n_{z}}z_{i}\Leftrightarrow
z\geq\sum_{i=0}^{n_{z}}z_{i}-n_{z}+1,z\leq z_{i},i\in\mathbb{N}_{n_{z}}$ (6a)
$z=\bigvee_{i=0}^{n_{z}}z_{i}\Leftrightarrow z\geq
z_{i},i\in\mathbb{N}_{n_{z}},z\leq\sum_{i=0}^{n_{z}}z_{i}.$ (6b)

The upper-bound constraints are necessary when the specification does not
include negation. $z_{t}^{\varphi}\in[0,1]$ is the variable that indicates
whether $(s,t)\models\varphi$. A recursive translation of an STL formula is as
follows:

$\varphi=\bigwedge_{i=1}^{n_{\varphi}}\varphi_{i}\Rightarrow
z_{t}^{\varphi}=\bigwedge_{i=1}^{n_{\varphi}}z_{k}^{\varphi_{i}};\varphi=\bigvee_{i=1}^{n_{\varphi}}\varphi_{i}\Rightarrow
z_{t}^{\varphi}=\bigvee_{i=1}^{n_{\varphi}}z_{t}^{\varphi_{i}};\\\
\varphi=G_{I}\psi\Rightarrow z_{t}^{\varphi}=\bigwedge_{t^{\prime}\in
I}z_{t^{\prime}}^{\psi};\varphi=F_{I}\psi\Rightarrow
z_{t}^{\varphi}=\bigvee_{t^{\prime}\in I}z_{t^{\prime}}^{\psi};\\\
\varphi=\psi_{1}U_{I}\psi_{2}\Rightarrow
z_{t}^{\varphi}=\bigvee_{t^{\prime}\in
I}(z_{t^{\prime}}^{\psi_{2}}\wedge\bigwedge_{t^{\prime\prime}\in[t,t^{\prime}]}z_{t^{\prime\prime}}^{\psi_{1}}),$
(7)

where $n_{\varphi}\in\mathbb{N}$. Given a formula $\varphi$, the set of
constraints recursively constructed by equations (4), (6), and (7) is denoted
by $\mathcal{C}_{\varphi}$.

######  Theorem 1 (Adapted from [2])

The following properties hold for the above mixed-integer linear program
encoding:(i) $(s,t)\models\varphi$, if adding $z_{t}^{\varphi}=1$ and
$\rho\geq 0$ to the constraints makes $\mathcal{C}_{\varphi}$ feasible, (ii)
$(s,t)\nvDash\varphi$, if adding $z_{t}^{\varphi}=1$ and $\rho\geq 0$ makes
$\mathcal{C}_{\varphi}$ infeasible, (iii) the largest $\rho$ such that
$z_{t}^{\varphi}=1$ and $\mathcal{C}_{\varphi}$ is feasible is equal to the
robustness.

It is shown in [2] that when the STL formulas are negation free, $\rho$ equals
robustness. As a result, it can be used as an objective function to maximize
robustness.

###  IV-B Set Containment for STL Formula Satisfaction

The objective of this subsection is to get the set of $z_{t}^{\pi}$s equal to
$1$, which are called active predicates, for the maximum robustness while
considering the nominal system. If the disturbance bound is small enough, it
can be assumed that the perturbed and nominal systems have the same set of
active predicates, and a closed-loop controller can be found to ensure that
the system’s reachability set still satisfies those predicates. We can do the
synthesis for the aggregate nominal system [2] rewritten as
$x(t+1)=A(t)x(t)+Bu(t)$ from (3), by using the STL satisfaction constraints
introduced before:

$\displaystyle\max_{x(t),u(t),z_{t}^{\pi},\rho}\quad-J(s[0,hrz(\varphi)])+\mathcal{M}(|\rho|-\rho)$
(8) s.t. $\displaystyle x(t+1)=A(t)x(t)+Bu(t),t\in\mathbb{N}_{hrz(\varphi)-1}$
$\displaystyle x(0)=[x^{initial}_{1},...,x^{initial}_{\eta}],$
$\displaystyle\mathcal{C}_{\varphi},z_{0}^{\varphi}=1.$

As long as robustness is positive, the proposed objective function minimizes
the user defined cost function $J(.)$, which can be a regular quadratic
function in the form of
$\sum_{t=0}^{hrz(\varphi)}x(t)^{T}Qx(t)+\sum_{t=0}^{hrz(\varphi)}u(t)^{T}Ru(t)$.
Otherwise, it maximizes robustness due to the effect of the large scalar
$\mathcal{M}$ and finds the nominal trajectory with the least violation.

Each active predicate is actually a set, $y_{t}\geq\rho,\forall y_{t}$, which
must hold for all possible signals at time $t$. By definition, we have
$s(t)\in\prod_{i}\Omega_{i}(t)\times\prod_{i}\Theta_{i}(t)$. Assuming the set
$\prod_{i}\Omega_{i}(t)\times\prod_{i}\Theta_{i}(t)$ is represented by a
zonotopic set $\mathcal{Z}(c,G)$ (notation $t$ is removed for readability),
then any possible signal must satisfy
$\mathrm{e}\geq\rho,\quad\forall\mathrm{e}\in\mathcal{Z}(p(c),p(G))$. Also, by
definition, the zonotope $\mathcal{Z}(c,G)$ has the following upper and lower
bounds $c-\sum_{i}|g_{i}|\leq\mathcal{Z}(c,G)\leq c+\sum_{i}|g_{i}|$, where
$g_{i}$ is the $i$th column of $G$. Using these bounds, the satisfaction
constraint for an active predicate would be:

$-p(c)+\sum_{i}|p(G)_{i}|\leq-\rho$ (9)

where $p(G)_{i}$ is the $i$th element of $p(G)$.

###### Theorem 2

The constraint in (9) can be written as a set of linear constraints as
follows:

$-p(c)+\sum_{i}p^{\prime}_{i}\leq-\rho,p^{\prime}_{i}\geq
p(G)_{i},p^{\prime}_{i}\geq-p(G)_{i}.$ (10)

###### Proof:

It can be easily seen that if such $p^{\prime}_{i}$s exist, the following
relation holds:

$-p(c)+\sum_{i}|p(G)_{i}|\leq-p(c)+\sum_{i}p^{\prime}_{i}\leq-\rho,$ (11)

which also satisfies the original constraint (9). Also, because (10) is the
relaxed form of the original problem, if such $p^{\prime}_{i}$s do not exist,
the original problem is also infeasible. ∎

Finally, the set of linear constraints that guarantees any possible
trajectories in viable and action sets satisfies the STL formula $\varphi$ is
denoted by $\mathcal{G}_{\varphi}$.

##  V Computation of Viable Sets under Additive Disturbance

The original problem has been transformed into a decentralized control
synthesis problem with zonotopic set containment constraints. The latter
problem was considered in [10], where a compositional approach using assume-
guarantee contracts is proposed. In this section, we give a brief overview of
[10] and incorporate the linear constraints $\mathcal{G}_{\varphi}$ into its
formulation.

###  V-A Decentralized Synthesis

First, the subsystems are decoupled from each other by considering the effects
of other subsystems as disturbances, and by making some assumptions on the
operational sets of each subsystem, as follows:

$x_{i}(t+1)=A_{ii}(t)x_{i}(t)+B_{ii}(t)u_{i}(t)+w_{i}^{aug}(t),$ (12)

where $w_{i}^{aug}(t)$ is the augmented disturbance set on subsystem $i$,
which belongs to:

$w_{i}^{aug}(t)\in\bigoplus_{j\neq
i}A_{ij}(t)\mathcal{X}_{j}(t)\oplus\bigoplus_{j\neq
i}B_{ij}(t)\mathcal{U}_{j}(t)\oplus W_{i}(t),$ (13)

where $\mathcal{X}_{j}(t)$ and $\mathcal{U}_{j}(t)$ are assumed operational
sets for the state and the control input of subsystem $j\in\mathcal{I}$. It
can be seen that the performance of each subsystem affects the assumptions of
the other subsystems. This give and take contracts are called assume-guarantee
contracts.

######  Definition 3 (Assume-Guarantee Contracts)

An assume-guarantee contract for subsystem $i\in\mathcal{I}$ is a pair
$\mathcal{C}_{i}=(\mathcal{A}_{i},\mathcal{G}_{i})$, where:

  * •

The assumption $\mathcal{A}_{i}$ is the assumption set over the disturbance
$w_{i}^{aug}(t)\in\mathcal{W}_{i}(t)$,

  * •

The guarantee $\mathcal{G}_{i}$ is the promise of subsystem $i$ over its state
and control input
$x_{i}(t)\in\mathcal{X}_{i}(t),u_{i}(t)\in\mathcal{U}_{i}(t)$.

As seen in (13), the following relation holds between the guarantee of other
subsystems $\mathcal{X}_{j},\mathcal{U}_{j},j\neq i$ and the assumption of
subsystem $i$, $\mathcal{A}_{i}$:

$\mathcal{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 W_{i}(t)$ (14)

The above zonotopic set is represented by
$W_{i}(t)=\mathcal{Z}(d^{w}_{i}(t),G^{w}_{i}(t))$, where
$d^{w}_{i}(t)\in\mathbb{R}^{n_{i}}$ and
$G^{w}_{i}(t)\in\mathbb{R}^{n_{i}\times l(t)}$. Next, we define a parametric
assume-guarantee contract, which is similar to the regular contract except
that the sets $\mathcal{X}_{i}(t)$, $\mathcal{U}_{i}(t)$ are replaced with the
parametric sets below:

$\mathcal{X}_{i}(t,\alpha_{i}^{x}(t)):=\mathcal{Z}(c_{i}^{x}(t),G_{i}^{x}\text{Diag}(\alpha_{i}^{x}(t))),$
(15a)
$\mathcal{U}_{i}(t,\alpha_{i}^{u}(t)):=\mathcal{Z}(c_{i}^{u}(t),G_{i}^{u}\text{Diag}(\alpha_{i}^{u}(t))),$
(15b)

where $G_{i}^{x}\in\mathbb{R}^{n_{i}\times f_{i}}$ and
$G_{i}^{u}\in\mathbb{R}^{m_{i}\times g_{i}}$ ($f_{i},g_{i}\in\mathbb{N}$) are
given matrices defined by the user, and the vectors
$c_{i}^{x}(t)\in\mathbb{R}^{n_{i}}$ ,
$\alpha_{i}^{x}(t)\in\mathbb{R}^{f_{i}}$, $c_{i}^{u}(t)\in\mathbb{R}^{m_{i}}$,
and $\alpha_{i}^{u}(t)\in\mathbb{R}^{g_{i}}$ are parameters. Also, the
parametric assumption set $\mathcal{W}_{i}(t,\alpha^{ext})$ is derived by
replacing the above parametric sets into equation (14), where $\alpha^{ext}$
denotes the set of all parameters. To deal with the mismatch between the
assumed and real operational disturbance sets, we introduce the notion of
correctness:

######  Definition 4 (Correctness)

A set of parametric contracts $\mathcal{C}_{i}$ is correct if

$\mathcal{W}_{i}(t)\subseteq\bigoplus_{j\neq i}A_{ij}(t)\Omega_{j}(t)\oplus
B_{ij}(t)\Theta_{j}(t)\oplus W_{i}(t),\forall i,t.$ (16)

The preceding definition is required to resolve the circularity problem of
assumption-guarantee contracts. It was shown in [10] that the following
sufficient constraints imply (16):

$\mathcal{X}_{i}(t,\alpha_{i}^{x}(t))\subseteq\Omega_{i}(t),\mathcal{U}_{i}(t,\alpha_{i}^{u}(t))\subseteq\Theta_{i}(t).$
(17)

The next step is to design a robust controller for each subsystem. The
following decentralized controller structure is proposed for each subsystem:

$x_{i}(t)=\bar{x}^{i}_{t}+T^{i}_{t}\zeta,u_{i}(t)=\bar{u}^{i}_{t}+M^{i}_{t}\zeta,\zeta\in\mathbb{B}_{k},$
(18)

where $\bar{x}^{i}_{t}\in\mathbb{R}^{n_{i}}$,
$\bar{u}^{i}_{t}\in\mathbb{R}^{m_{i}}$, $T^{i}_{t}\in\mathbb{R}^{n_{i}\times
q_{i}(t)}$, and $M^{i}_{t}\in\mathbb{R}^{m_{i}\times q_{i}(t)}$ are unknowns
that need to be tuned and $q(t)=k+\sum_{i=0}^{i=t}{l(i)}$, where
$k\in\mathbb{N}$ is a hyper-parameter. Then, for subsystem $i$, it can be
shown that the following linear constraints are sufficient for tuning the
control parameters:

$[A_{ii}(t)T^{i}_{t}+B_{ii}(t)M^{i}_{t},G^{aug}_{i}(t)]=[T^{i}_{t+1}],t\in\mathbb{N}_{h-1}$
(19a)
$A_{ii}(t)\bar{x}^{i}_{t}+B_{ii}(t)\bar{u}^{i}_{t}+d^{aug}_{i}(t)=\bar{x}^{i}_{t+1},t\in\mathbb{N}_{h-1}.$
(19b)

If such parameters exist,
$\Omega_{i}(t)=\mathcal{Z}(\bar{x}^{i}_{t},T^{i}_{t})$ is the viable set,
$\Theta_{i}(t)=\mathcal{Z}(\bar{u}^{i}_{t},M^{i}_{t})$ is the action set, and
(18) is the controller. Intuitively, constraints (19a) and (19b) are set
containment constraints; (19b) adjusts the centers of the viable sets and
(19a) takes care of the set expansion at each step, such that all the possible
trajectories are contained within the tube
$\Omega_{i}(0),\Omega_{i}(1),\cdots,\Omega_{i}(h)$, where $h$ is the horizon,
which is set to $hrz(\varphi)$ in our problem. Additionally, the following
constraints are proposed to impose hard constraints:

$\mathcal{Z}(\bar{x}^{i}_{t},T^{i}_{t})\subseteq
X_{i}(t),\mathcal{Z}(\bar{u}^{i}_{t},M^{i}_{t})\subseteq
U_{i}(t),t\in\mathbb{N}_{h}.$ (20)

It was demonstrated in [25] that zonotope and polytope containment problems
can be encoded into linear constraints. Thus, all of the suggested constraints
(16), (19), (20), and $\mathcal{G}_{\varphi}$ for all subsystems and time
steps may be merged to build a centralized linear program to solve Problem 1.
The objective function is ad-hoc, but we recommend the mean square error
between the center line of viable/action sets and the nominal
trajectory/controllers generated by (8).

###  V-B Compositional Computation of Decentralized Viable Sets

Despite the fact that the centralized solution presented at the end of the
preceding subsection is a linear program, it still suffers from curse of
dimensionality in high dimensions. Nevertheless, it is demonstrated in [10]
that the suggested parameterization (15) allows for compositional computation
of viable sets in a time-efficient manner by transforming a single, large
linear program into a group of smaller linear programs. We show that if the
STL formula in Problem 1 is separable by subsystems, we can also use the
parameterization to solve the same centralized approach in the previous
section in a compositional manner. Additionally, convergence is ensured due to
the convexity of the problem set.

###### Assumption 1

The STL formula in Problem 1 is separable by the subsystems, meaning it can
take the form $\varphi=\varphi_{1}\wedge...\wedge\varphi_{\eta}$, where $i$ is
the subsystem’s index.

In [10], we proposed a parametric potential function that quantifies how far a
set of contracts is from correctness. This comes in contrast to the previously
introduced correctness property, which was either true or false. The larger
the parametric potential function, the farther the set of contracts is from
correctness, so the goal is to minimize the proposed potential function. Here,
the parametric potential function is modified by including the containment
constraints coming from the STL formulas, as well as adding the sum of the
directed Hausdorff distances between the hard constraints and the
viable/action sets in (20) into the potential function.

######  Definition 5 (Parametric potential function)

The parametric potential function $\mathcal{V}(\alpha^{ext})$ is defined as
$\mathcal{V}(\alpha^{ext})=\sum_{i\in\mathcal{I}}\mathcal{V}_{i}(\alpha^{ext})$,
where

$\mathcal{V}_{i}(\alpha^{ext}):=\\\
\sum_{t\in\mathbb{N}_{hrz(\varphi)}}[d_{DH}(\mathcal{X}_{i}(t),\Omega_{i}(t))+d_{DH}(\mathcal{U}_{i}(t),\Theta_{i}(t))\\\
+d_{DH}(X_{i}(t),\Omega_{i}(t))+d_{DH}(U_{i}(t),\Theta_{i}(t))].$ (21)

Using the technique explained in Subsection IV-B, the satisfaction of the STL
formula $\varphi_{i}$ for subsystem $i$ can be encoded as a set of linear
constraints denoted by $\mathcal{G}_{\varphi_{i}}$. Each component of the
parametric potential function $\mathcal{V}_{i}(\alpha^{ext})$ can be computed
using these constraints and (19) by solving the following linear program:

$\displaystyle\mathcal{V}_{i}(\alpha)=\min_{\underset{,d_{t}^{x},d_{t}^{u},\bar{d}_{t}^{x},\bar{d}_{t}^{u}}{\mathrm{x}^{i},T^{i},\mathrm{u}^{i},M^{i}}}\begin{aligned}
&\sum_{t\in\mathbb{N}_{hrz(\varphi)}}[d^{x}_{t}+\bar{d}^{x}_{t}]+\sum_{t\in\mathbb{N}_{hrz(\varphi)-1}}[d^{u}_{t}+\bar{d}^{u}_{t}]\end{aligned}$
subject to
$\displaystyle[A_{ii}(t)T^{i}_{t}+B_{ii}(t)M^{i}_{t},G_{i}^{w}(t)]=[T^{i}_{t+1}],\forall
t\in\mathbb{N}_{hrz(\varphi)-1}$ (22a) $\displaystyle
A_{ii}(t)\bar{\mathrm{x}}^{i}_{t}+B_{ii}(t)\bar{\mathrm{u}}^{i}_{t}+d_{i}^{w}(t)=\bar{\mathrm{x}}^{i}_{t+1},\forall
t\in\mathbb{N}_{hrz(\varphi)-1}$ (22b)
$\displaystyle\mathcal{Z}(\bar{\mathrm{x}}^{i}_{t},T^{i}_{t})\subseteq\
\mathcal{X}_{i}(t,\alpha^{x}_{i}(t))\oplus\mathcal{Z}(0,d^{x}_{t}I_{n_{i}}),\forall
t\in\mathbb{N}_{hrz(\varphi)}$ (22c)
$\displaystyle\mathcal{Z}(\bar{\mathrm{u}}^{i}_{t},M^{i}_{t})\subseteq\
\mathcal{U}_{i}(t,\alpha^{u}_{i}(t))\oplus\mathcal{Z}(0,d^{u}_{t}I_{m_{i}}),\forall
t\in\mathbb{N}_{hrz(\varphi)-1}$ (22d)
$\displaystyle\mathcal{Z}(\bar{\mathrm{x}}^{i}_{t},T^{i}_{t})\subseteq
X_{i}(t)\oplus\mathcal{Z}(0,\bar{d}^{x}_{t}I_{n_{i}}),\forall
t\in\mathbb{N}_{hrz(\varphi)}$ (22e)
$\displaystyle\mathcal{Z}(\bar{\mathrm{u}}^{i}_{t},M^{i}_{t})\subseteq
U_{i}(t)\oplus\mathcal{Z}(0,\bar{d}^{u}_{t}I_{m_{i}}),\forall
t\in\mathbb{N}_{hrz(\varphi)-1}$ (22f)
$\displaystyle\mathcal{G}_{\varphi_{i}},\bar{\mathrm{x}}^{i}_{0}=x_{i}^{initial}(0)$
(22g) $\displaystyle d_{t}^{x},\bar{d}_{t}^{x}\geq 0,\hskip 5.69054pt\forall
t\in\mathbb{N}_{hrz(\varphi)},$ (22h) $\displaystyle
d_{t}^{u},\bar{d}_{t}^{u}\geq 0,\hskip 5.69054pt\forall
t\in\mathbb{N}_{hrz(\varphi)-1}.$ (22i)

Constraints (22a) and (22b) originate from (19). Also, the STL satisfaction
constraints and the initial state constraint are added in (22g). The remaining
constraints along with the objective function compute an over-approximation
for the summation of the directed Hausdorff distance between sets
$\Omega_{i}(t)$ and $\mathcal{X}_{i}(t)/X_{i}(t)$ and $\Theta_{i}(t)$ and
$\mathcal{U}_{i}(t)$/$U_{i}(t)$ over all time steps. This approach of
computing the Directed Hausdorff distance is inspired from [25].

######  Theorem 3 (Convexity of the potential function)

The potential function proposed above is convex with respect to the
parameters. The set of acceptable parameters (correct and valid) is also a
convex set.

###### Proof:

As seen in (22), each component of the potential function is a linear program,
which makes $\mathcal{V}_{i}(\alpha^{ext})$ a convex and piecewise affine
function (a sum of convex functions is convex). Also, it is a well-known fact
that the level set of a convex function is a convex set, thus, the set of
acceptable parameters, which is equal to the zero level set of the potential
function, is also a convex set. ∎

The idea is to minimize the potential function using gradient descent and
iteratively update the parameters:

$\alpha^{ext}\leftarrow\alpha^{ext}-\sum_{i\in\mathcal{I}}\nabla_{\alpha^{ext}}\mathcal{V}_{i}(\alpha^{ext}),$
(23)

Convergence to the global minimum is guaranteed because the proposed potential
function is convex. Each subsystem can find the direction that is best for it
($\nabla_{\alpha^{ext}}\mathcal{V}_{i}(\alpha^{ext})$), using its own local
information and the common knowledge parameters, breaking the problem down
into many smaller linear programs. If the minimum of the potential function is
zero (by definition, the potential function is always larger than zero), it
indicates that both the set of derived parametric contracts are correct and
the viable and actions sets are within hard constraints. Thus, the desired
control policies and viable sets are determined. Also, the nominal
trajectories and controllers derived from (8) can be used as initial values
for the center parameters in our parameterized sets to give the gradient
descent a warm start.

##  VI Case Study

We apply the method developed in this paper to address the load-frequency
problem in power networks [3]. A network is made up of several areas, each
with its own power generator and demands, and some of them can be connected to
each other to interchange power as needed, depending on the network
architecture. Each area’s state is represented by a 2-dimensional vector
$[\delta_{i}(t),f_{i}(t)]^{T}$, where $\delta_{i}(t)\in\mathbb{R}$ is the
deviation of the phase angle and $f_{i}(t)\in\mathbb{R}$ is the deviation of
the frequency at time $t$ for area $i\in\mathbb{N}$. Also,
$u_{i}(t)\in\mathbb{R}$ is the control input, which is the amount of change
from its nominal value in the power generated by the generator at the area $i$
and time $t$. The dynamics for each area is given by:

$\dot{\delta}_{i}(t)=2\pi f_{i}(t)\hskip 5.69054pt,\hskip
5.69054pt\dot{f}_{i}(t)=-\dfrac{f_{i}(t)}{T_{p_{i}}}+\dfrac{K_{p_{i}}u_{i}(t)}{T_{p_{i}}}-\\\
\dfrac{K_{p_{i}}}{2\pi
T_{p_{i}}}(\sum_{j\in\mathcal{N}_{i}}K_{s_{ij}}[\delta_{i}(t)-\delta_{j}(t)])-\dfrac{K_{p_{i}}\omega_{i}(t)}{T_{p_{i}}},$
(24)

where $K_{p_{i}},K_{s_{ij}},T_{p_{i}}$ are the system gain, synchronizing
coefficient between area $i$ and $j$, and system model time constant. In this
case study, they are set to 110, 0.5, and 25, respectively, for all areas.
Also, $\omega_{i}(t)$ is the load disturbance for area $i$ at time $t$, which
is bounded by $|\omega_{i}(t)|\leq 0.001$. In addition, $\mathcal{N}_{i}$
denotes the neighbours of area $i$. Here, we consider the ring network
architecture consisted of $20$ areas. Also, the control input is bounded by
$|u_{i}(t)|\leq 0.1$. We use the Euler method to discretize the dynamics for
every $0.1$ unit of time. For all areas, the initial state is $[0.1,0.1]^{T}$
and the STL specification is
$\varphi_{i}=\textbf{F}_{[0,6]}\textbf{G}_{[0,2]}\psi_{1}\wedge\textbf{F}_{[0,8]}\psi_{2}$,
where $\psi_{1}=[\delta_{i}\leq 0.26]\wedge[\delta_{i}\geq
0.14]\wedge[f_{i}\leq-0.04]\wedge[f_{i}\geq-0.16]$ and
$\psi_{2}=[\delta_{i}\leq 0.01]\wedge[\delta_{i}\geq-0.01]\wedge[f_{i}\leq
0.01]\wedge[f_{i}\geq-0.01]$. The goal is to synthesize decentralized
controllers for each area subject to the specifications. We set the horizon to
nine and synthesize the controllers using our approach. The baseline
parametric sets are selected to be the viable and action sets generated from
(19) while couplings to other areas are ignored. The initial value of all
parameters in the distributed algorithm is one. We used Gurobi on a MacBook
Pro with 2.6 GHz 6-Core Intel Core i7 and 16 GB memory to run the algorithm.
The results are shown in Fig.1(a), and Fig.1(b). It can be seen that any
possible trajectory that passes through the viable sets satisfies the STL
specification and at the same time all the implemented controllers satisfy the
hard constraint on the control input.

To demonstrate the approach’s scalability, we experimented with various number
of areas in the ring network and reported the running time in Fig 1(c). The
stated time period includes only the time spent on second step, but not on
solving the MILP. That is because both distributed and centralized approaches
share the first step. Additionally, to ensure that the solution exists for
high-dimensional state spaces, we consider a large bound for the
controller(i.e. $|u_{i}(t)|\leq 10$). As predicted, the distributed approach
has a slower growth rate, making it more appropriate for the state spaces
larger than $40$. Moreover, one of the primary benefits of the distributed
technique is that it may be calculated in parallel. While we handled
everything sequentially here, if multiprocessing is employed, the stated time
could be reduced more depending on the number of cores used.

(a)

(b)

(c)

Figure 1: (a) The green sets illustrate the viable sets for one of the areas
in the case study. The blue and orange sets are the set of states satisfying
$\psi_{1}$ and $\psi_{2}$, respectively. The red sets show the parameterized
sets defined on the state space at different time steps for this specific area
(some of them are tightly close to the viable sets and are not visible). The
black line represents the trajectory traveled by this area. (b) Controllers
for each area as time series. (c) Reported time in seconds for the distributed
and centralized approach for different state space dimensions.

##  VII CONCLUSIONS

Control synthesis subject to both a STL formula and a bounded disturbance is a
computationally challenging problem. To overcome this challenge, we propose a
solution which consists of two steps: First, we convert satisfaction of the
STL formula into a set containment problem. To handle it, we consider the
nominal system and use a centralized MILP. We claim that for small enough
disturbances, both systems would have the same set of active predicates, which
are seen as bounds. Second, we synthesize controllers subject to these bounds.
Since the second step needs a set-based calculation, it has a relatively
higher computational cost and thus creates a bottleneck for large scale
systems. We show that this step can be achieved in a compositional fashion
when the STL formula is separable by subsystems.

In the future, we will investigate the possibility of replacing the MILP in
the first step, which remains a barrier due to its computational cost.
Sampling methods are a promising direction. Additionally, we are considering
employing a distributed control architecture instead of the decentralized
architecture proposed here.

## References

  * [1] Christel Baier and Joost-Pieter Katoen.  Principles of model checking.  MIT press, 2008.
  * [2] Calin Belta and Sadra Sadraddini.  Formal Methods for Control Synthesis: An Optimization Perspective.  Robotics, and Autonomous Systems Annu. Rev. Control Robot. Auton. Syst, 28:12–13, 2018.
  * [3] E. Camponogara, D. Jia, B.H. Krogh, and S. Talukdar.  Distributed model predictive control.  IEEE Control Systems Magazine, 22(1):44–52, 2002.
  * [4] Krishnendu Chatterjee and Thomas A Henzinger.  Assume-guarantee synthesis.  In International Conference on Tools and Algorithms for the Construction and Analysis of Systems, pages 261–275. Springer, 2007.
  * [5] Yuxiao Chen, James Anderson, Karan Kalsi, Steven H Low, and Aaron D Ames.  Compositional set invariance in network systems with assume-guarantee contracts.  In 2019 American Control Conference (ACC), pages 1027–1034. IEEE, 2019.
  * [6] Jorge Cortes, Sonia Martinez, Timur Karatas, and Francesco Bullo.  Coverage control for mobile sensing networks.  IEEE Transactions on robotics and Automation, 20(2):243–255, 2004\.
  * [7] Alexandre Donzé and Oded Maler.  Robust satisfaction of temporal logic over real-valued signals.  In International Conference on Formal Modeling and Analysis of Timed Systems, pages 92–106. Springer, 2010.
  * [8] Souradeep Dutta, Xin Chen, and Sriram Sankaranarayanan.  Reachability analysis for neural feedback systems using regressive polynomial rule inference.  In Proceedings of the 22nd ACM International Conference on Hybrid Systems: Computation and Control, pages 157–168, 2019.
  * [9] Magnus Egerstedt and Xiaoming Hu.  Formation constrained multi-agent control.  IEEE transactions on robotics and automation, 17(6):947–951, 2001\.
  * [10] 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, pages 1–10, 2020.
  * [11] Antoine Girard.  Reachability of uncertain linear systems using zonotopes.  In International Workshop on Hybrid Systems: Computation and Control, pages 291–305. Springer, 2005.
  * [12] Sertac Karaman, Ricardo G Sanfelice, and Emilio Frazzoli.  Optimal control of mixed logical dynamical systems with linear temporal logic specifications.  In 2008 47th IEEE Conference on Decision and Control, pages 2117–2122. IEEE, 2008.
  * [13] 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.
  * [14] Weixuan Lin and Eilyan Bitar.  Decentralized control of constrained linear systems via assume-guarantee contracts.  In 2020 American Control Conference (ACC), pages 917–924. IEEE, 2020.
  * [15] Lars Lindemann and Dimos V Dimarogonas.  Control barrier functions for multi-agent systems under conflicting local signal temporal logic tasks.  IEEE control systems letters, 3(3):757–762, 2019.
  * [16] Zhiyu Liu, Bo Wu, Jin Dai, and Hai Lin.  Distributed communication-aware motion planning for multi-agent systems from stl and spatel specifications.  In 2017 IEEE 56th Annual Conference on Decision and Control (CDC), pages 4452–4457. IEEE, 2017.
  * [17] Anirudha Majumdar and Russ Tedrake.  Funnel libraries for real-time robust feedback motion planning.  The International Journal of Robotics Research, 36(8):947–982, 2017\.
  * [18] Oded Maler and Dejan Nickovic.  Monitoring temporal properties of continuous signals.  In Formal Techniques, Modelling and Analysis of Timed and Fault-Tolerant Systems, pages 152–166. Springer, 2004.
  * [19] 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.
  * [20] Pierluigi Nuzzo, Alberto L Sangiovanni-Vincentelli, Davide Bresolin, Luca Geretti, and Tiziano Villa.  A platform-based design methodology with contracts and related tools for the design of cyber-physical systems.  Proceedings of the IEEE, 103(11):2104–2132, 2015.
  * [21] Chanwook Oh, Eunsuk Kang, Shinichi Shiraishi, and Pierluigi Nuzzo.  Optimizing assume-guarantee contracts for cyber-physical system design.  In 2019 Design, Automation & Test in Europe Conference & Exhibition (DATE), pages 246–251. IEEE, 2019.
  * [22] Reza Olfati-Saber, J Alex Fax, and Richard M Murray.  Consensus and cooperation in networked multi-agent systems.  Proceedings of the IEEE, 95(1):215–233, 2007.
  * [23] Yash Vardhan Pant, Houssam Abbas, Rhudii A Quaye, and Rahul Mangharam.  Fly-by-logic: Control of multi-drone fleets with temporal logic objectives.  In 2018 ACM/IEEE 9th International Conference on Cyber-Physical Systems (ICCPS), pages 186–197. IEEE, 2018.
  * [24] Vasumathi Raman, Alexandre Donzé, Mehdi Maasoumy, Richard M Murray, Alberto Sangiovanni-Vincentelli, and Sanjit A Seshia.  Model predictive control with signal temporal logic specifications.  In 53rd IEEE Conference on Decision and Control, pages 81–87. IEEE, 2014.
  * [25] Sadra Sadraddini and Russ Tedrake.  Linear encodings for polytope containment problems.  In 2019 IEEE 58th Conference on Decision and Control (CDC), pages 4367–4372. IEEE, 2019.


# To the synthesis and characterization of layered metal phosphorus
triselenides proposed for electrochemical sensing and energy applications
(Manuscript was initially submitted to ACS Catalysis as a Comment for Ref. 22;
however was rejected after appeal)

Yuriy Dedkov<EMAIL_ADDRESS>Mouhui Yan  Elena Voloshina
<EMAIL_ADDRESS>Department of Physics, Shanghai University, 200444
Shanghai, China  Institute of Physical and Organic Chemistry, Southern Federal
University,
344090 Rostov on Don, Russia

###### Abstract

Recent studies reported on the synthesis and characterization of several bulk
crystals of layered metal triselenophosphites MPSe3 (M = transition metals).
In these works characterization was performed via a combination of different
bulk- and surface-sensitive experimental methods accompanied by DFT
calculations. However, the critical examination of the available experimental
and theoretical data demonstrates that these results do not support the
conclusions on the electrochemical sensing and energy applications of studied
triselenophosphites. These conclusions are made without any relation to the
age of discussed data and possible recent progress in experimental and
theoretical approaches.

###### keywords:

trichalcogenides, DFT, XRD, EDS, XPS

##  1 Introduction

The discovery of the unique transport properties of graphene in in 2004 [1, 2]
lead to the increased attention to graphene-based systems [3, 4, 5, 6] and
also to other classes of 2D materials. Among them are h-BN [7], two-
dimensional dichalcogenides [8, 9], and many others, as well as heterosystems
on the basis of these materials [10]. These studies led to the discovery of
many exciting properties, which can be used in different applications, like,
touch screens [11, 12], gas sensors [13, 14, 15], very good thermal conductors
[16, 17], etc.

Recently, a new class of low-dimensional materials, namely transition metal
phosphorus trichalcogenides received a lot of attention, because they are
considered as a new class of 2D materials, which can find application in
electronics, sensing, catalysis, or energy conversion. Many experimental and
theoretical works as well as several respective review articles on the studies
of the transition metal phosphorus trichalcogenides appeared in literature
[18, 19, 20, 21, 22, 23, 24]. However, the available in literature
experimental and theoretical works on the studies of, e. g. electrocatalytic
properties of these materials, contain sometimes misleading information on the
synthesis and characterization. These inaccurate results obtained during
experiment/sample preparation and treatment of experimental data lead to the
fact that the discussed effects and main conclusions are not supported by the
presented experimental data.

Recently a work of Gusmão et al. [22] was published, which is devoted to the
studies of the electrochemical performance (hydrogen and oxygen evolution
reactions) of transition metal phosphorus trichalcogenides MPSe3 (M = Cd, Cr,
Fe, Mn, Sn, Zn). It was found that among all synthesised samples, FePSe3 and
MnPSe3 have the highest efficiency for the hydrogen evolution reaction with
good stability. A the same time it was found that MnPSe3 holds the lowest
oxidation potential, although it was assigned to the presence of MnO2 in the
sample. In the same work MPSe3 samples were synthesised using standard
chemical vapour transport method. Structural characterization was performed by
means of x-ray diffraction (XRD) in the $\theta-2\theta$ geometry using Cu
K$\alpha$ line. Morphology and chemical analysis were studied using scanning
electron microscopy (SEM) combined with energy-dispersive x-ray spectroscopy
(EDS) applying the gentle electron beam with the energy of $2$ keV. X-ray
photoelectron spectroscopy (XPS) was used for the surface analysis of the
studied samples. These data obtained during samples characterization were used
to support the observed effects in further experiments on the electrochemical
performance of these materials (hydrogen and oxygen evolution reactions).
However, as shown below, the presented claims in the discussed manuscript have
to be reconsidered, because experimental as well as theoretical data on the
sample characterization contain several flaws and errors and these results
have to be fully reanalyzed. Therefore the main conclusions on the
electrochemical performance of transition metal phosphorus trichalcogenides
have to be reconsidered.

##  2 Experimental and computational details

Samples synthesis. Manganese ($99.9$%), phosphorus ($99.999$%), and selenium
($99.999$%) from Shanghai Macklin Biochemical Co., Ltd. and Alfa Aesar were
used during synthesis. A stoichiometric amount of high-purity elements (mole
ratio $\mathrm{Mn}:\mathrm{P}:\mathrm{Se}=1:1:3$, $1$ g in total) and iodine
(about $20$ mg) as a transport agent were sealed into a quartz ampule (length
$17$ cm, external diameter approximately $15$ mm) and kept in a two-zone
furnace ($650-600^{\circ}$ C). The pressure inside the ampule was pumped down
to $1\times 10^{-3}$ Torr. After 10 days of heating, the ampule was cooled
down to room temperature with bulk crystals in the colder edge.

Characterization. XRD patterns were collected with a Bruker D2 Phaser
diffractometer using Cu K$\alpha$ ($1.54178$ Å) radiation at room temperature.
Optical images were collected with optical microscope at different
magnifications. Presented reference spectra of WSe2 were collected using SPECS
PHOIBOS 150 energy analyzer and Al K$\alpha$ monochromatized x-ray source
($h\nu=1486.6$ eV).

DFT calculations. Spin-polarised DFT calculations based on plane-wave basis
sets of $500$ eV cutoff energy were performed with the Vienna ab initio
simulation package (VASP) [25, 26]. The Perdew-Burke-Ernzerhof (PBE) exchange-
correlation functional [27] was employed. The electron-ion interaction was
described within the projector augmented wave (PAW) method [28] with Mn ($3p$,
$3d$, $4s$), Cr ($3p$, $3d$, $4s$), Fe ($3p$, $3d$, $4s$), P ($3s$, $3p$), and
Se ($4s$, $4p$) states treated as valence states. The Brillouin-zone
integration was performed on $\Gamma$-centred symmetry reduced Monkhorst-Pack
meshes using a Gaussian smearing with $\sigma=0.1$ eV, except for the
calculation of density of states (DOS). For these calculations, the
tetrahedron method with Blöchl corrections [29] was employed. A $12\times
12\times 4$ $k$-mesh was used in the case of ionic relaxations and $24\times
24\times 8$ for the DOS calculations, respectively. The DFT+$\,U$ scheme [30,
31] was adopted for the treatment of Mn-, Cr- and Fe-$3d$ orbitals. Dispersion
interactions were considered adding a $1/r^{6}$ atom-atom term as
parameterised by Grimme (“D2” parameterisation) [32]. During structure
optimisation, the convergence criteria for energy and force were set equal to
$10^{-5}$ eV and $1\times 10^{-2}$ eV/Å, respectively.

##  3 Results and discussions

The studied 3D MPSe3 crystals have either $C2/m$ (M = Cr) or $R\bar{3}$ (M =
Mn, Zn, Cd, Fe) space group symmetry. The both cases can be viewed as layered
structures, where each stacked MPSe3 layer has a $D_{3d}$ symmetry (see Fig.
1). Therefore, one might expect the hexagonal-like shape for the grown
crystals with angles between crystallite edges either $120^{\circ}$ or
$60^{\circ}$ [20]. Analysis of the presented in Ref. 22 SEM images of the
grown crystals demonstrate the absence of this representative feature,
indicating the low quality of the studied crystals (see Fig. 2(a) for MnPSe3;
see also Figs. 2 and S1 in Ref. 22 for other MPSe3 crystals). At the same
time, we present in Fig. 2(b,c) optical images of MnPSe3 crystal grown
recently in our laboratory. One can clearly see that this sample demonstrates
very high quality with well ordered planes and steps oriented with respect to
each other either by $120^{\circ}$ or $60^{\circ}$. The measured XRD patterns
of our crystals demonstrate extraordinary quality without additional phases
(Fig. 2(d)) [20], whereas the data presented in Ref. 22 (Fig. S3) demonstrate
the existence of other undesired crystal phases in the studied samples.

The SEM/EDS combination was used in Ref. 22 to study the morphology as well as
the composition of the studied MPSe3 samples. Despite the low energy of the
used electron beam during this analysis, the EDS method cannot be considered
as a surface sensitive as, for example, at $5$ keV of beam energy the probing
depth can reach $0.5\,\mu\mathrm{m}$. Therefore, the results presented in Fig.
S2 and Tab S1 of Ref. 22 demonstrate extremely poor quality of the studied
samples in bulk as well as at the surface: some samples demonstrate the
absence of phosphorus in the sample and the level of the C- and
O-contamination is very high (it is varied between $\approx 43$ at.% and
$\approx 65$ at.% in total).

The electronic structure of the studied MPSe3 crystals was investigated in
Ref. 22 using density functional theory (DFT) within the PBE+$U$ approach with
$U=3$ eV for the proper treatment of electron correlations in the valence
band. Here we would like to point out, that from the presented in Ref. 22 data
it is absolutely unclear if bulk 3D or monolayer 2D phases were considered in
the DFT calculations. This is crucial for the description of the band gap
around the Fermi level as its value strongly depends on the system
dimensionality as shown in Ref. 33. Thus, only comparison of experimental
results with theoretical data obtained for the same phase (here, bulk 3D) will
make sense. The calculated values for the band gap are summarized in Table 1,
where we compare data from Ref. 22 and our own data obtained in the framework
of the present study. Fig. 3 shows density of states (DOS) plots for a series
of MPSe3 compounds calculated in the present study for the 3D phase and
different values of $U$.

As can be seen, the values of band gap for MPSe3 in Ref. 22 are approximately
by factor of $2$ smaller compared to the experimental values, which appeared
in earlier publications [20] and to the values obtained with DFT in the
present study. These values are also not consistent with the recently
published values for MnPSe3 [33], which are in good agreement with available
experimental data. It is also very confusing, that in Ref. 22, the band gap of
$0.5$ eV is claimed for FePSe3, that contradicts to the metallic state for
this compound as can be seen in Fig. 3 of this reference. Moreover, high
density of states at the Fermi level indicates an unstable situation, which is
not energetically favourable.

When looking at the DFT results, i. e. DOS calculated for MPSe3, one can note
a significant variation of a band gap depending on the nature of M. Comparison
of the calculated band gap with, e.g., the free energy of water splitting
($1.23$ eV) may be used as an indication that MnPSe3 (with $\Delta E_{g}=1.8$
eV) is a good candidate for water-splitting catalyst. Such a consideration,
however, is oversimplified, because for the strait forward conclusion one has
to take into account some important factors. First of all, any catalytic
process with MPX3 will take place at the surface and it is known, that the
electronic properties of these compounds strongly depend on the system
dimensionality [33]. Secondly, one has to consider presence of defects, which
on the one hand are expected to enhance the catalytic activity [34], on the
other hand will influence (reduce) the band gap [33]. In the case of
electrochemical water splitting, the band edges of the catalysts must straddle
the redox potentials of water, which in their turn depend on the $pH$ value
[21]. Thus, in order to discuss the activity of MPSe3 for the electrochemical
applications, one has to perform very detailed, accurate, and systematic
study, where the calculations on 3D MPSe3 are just the initial step. At the
same timeIt is necessary to note, that the scientific weight of the drawn
conclusions will depend on the reliability of this initial step descriptions.

The most confusing part of Ref. 22 regards to the XPS characterization of the
studied MPSe3 samples and their interpretation. Because these XPS data are
used to support the respective electrochemical performance results for the
studied material, we believe that the main conclusions of the discussed
manuscript have to be reconsidered and the respective data have to be fully
reanalized.

To remind, here we present the basic consideration for the analysis of XPS
emission lines. Core levels in XPS use the notation $nl_{j}$, where $n$ is the
principal quantum number, $l$ is the angular momentum quantum number and
$j=l+s$ (where $s=\pm 1/2$ is the spin angular momentum number). All orbital
levels, except the $s$ levels ($l=0$), give rise to a doublet with the two
possible states having different binding energies. This is known as spin-orbit
splitting. The peaks will also have specific area ratios based on the
degeneracy of each spin state, i.e. the number of different spin combinations
that can give rise to the total $j$. For example, for the $2p$ spectra, where
$n=2$ and $l=1$, $j$ will be $1/2$ and $3/2$. The area ratio for the two spin
orbit peaks ($2p_{1/2}:2p_{3/2}$) will be $1:2$ (corresponding to $2$
electrons in the $2p_{1/2}$ level and $4$ electrons in the $2p_{3/2}$ level).
The similar consideration is valid for $d$ and $f$ levels, where the area
ratios for the respective spin-orbit split components are $2:3$ and $3:4$,
respectively. These ratios must be taken into account when analyzing spectra
of the $p$, $d$ and $f$ core levels. Spin-orbit splitting values (in eV) can
be found in different databases (e. g.:
https://xpssimplified.com/periodictable.php). Here, we present two
representative examples, demonstrating the application of the fit procedure
for the Se $3d$ spectra, where the above described peaks’ ratios are used. The
Se $3d$ spectrum for a WSe2 single crystal measured using monochromotized Al
K$\alpha$ x-ray source and shown in Fig. 4(a) demonstrates clear spin-orbit
split doublet with components having equal full width at half maximum (FWHM),
peak separation of $0.85$ eV and intensities ratio of
$I(3d_{3/2})/I(3d_{5/2})=0.67$. The Se $3d$ spectrum for $\alpha$-P4Se3
measured using non-monochromatized Mg K$\alpha$ x-ray source shows that all
spin-orbit doublets must be fit in order to properly identify the species
present in the sample. The $3d_{3/2}$ and $3d_{5/2}$ doublet for each chemical
specie is constrained to have $\approx 2:3$ peak area ratios, equal FWHM, and
a peak separation of $\approx 1$ eV (Fig. 4(b)) [35]. The discussed examples
demonstrate the universality of the discussed approach without any connection
to the specific sample.

Coming back to the XPS data presented in Ref. 22 we have to mention the high
level of the C- and O-contamination that does not allow to carefully study the
chemical states of elements in the studied compounds (see Fig. S5 of the
discussed manuscript). Therefore, it is very difficult to draw the clear
conclusions, which might support the further data obtained in the studies of
the electrochemical performance of these materials. Moreover, the performed
analysis and interpretation of the respective core levels presented in this
manuscript does not correspond to the high standards accepted in the
scientific publications. As an example, we present here several XPS spectra
for MnPSe3, FePSe3, and ZnPSe3 extracted from Fig. 5 of Ref. 22. The following
list marks the serious deficiencies in the spectra interpretation (Fig.
4(c-e)):

  * 1.

The respective fit of the Se $3d$ XPS lines is not correct for all data
presented in the discussed manuscript. As discussed before, the ratio of
integral intensities for $3d_{3/2}$ and $3d_{5/2}$ lines have to be $\approx
0.67$ and cannot be more than $1$ as it is presented for the case of ZnPSe3.
If necessary, several spin-orbit split lines corresponding to different
species have to be introduced for the proper fit of the XPS spectra.

  * 2.

The respective fit of P $2p$ XPS line is not correct for several data sets
presented in the discussed manuscript (see the respective spectra for FePSe3
and ZnPSe3). As discussed before, the ratio of integral intensities for
$2p_{1/2}$ and $2p_{3/2}$ lines have to be $\approx 0.5$ and cannot be more
than $1$ as it is presented for the case of ZnPSe3. If necessary, several
spin-orbit split lines corresponding to different species have to be
introduced. It has to be also noted that the same broad peak in the P $2p$
spectra for different MPSe3 compounds, which is located at $\approx 134$ eV of
binding energy, is assigned to different emission lines: P4O10 or P
$2p_{1/2}$, although it is obvious that this component corresponds to the same
chemical state of P atoms, namely to phosphorus oxide (PxOy). Ironically, on
the basis of this inaccurate fit, authors of Ref. 22 made the conclusion on
the low concentration of PxOy for FePSe3 and on the absence of the oxidation
state of P for ZnPSe3, which is, of course, wrong.

  * 3.

For the Mn $2p_{3/2}$ XPS line, two components in the spectra were assigned to
the emission from Mn atoms either in MnPSe3 or in MnO2 (Fig. 4(c)). However,
this statement is not supported by any additional reference XPS measurements.
The quality of the XPS spectra presented in Ref. 22 is very poor, not allowing
to perform clear identification of the components and several Mn-oxide phases
can be assigned to this “MnO2”-peak: MnO, Mn2O3, or MnO2 [36].

  * 4.

The presented fit of the Fe $2p$ XPS spectra is not correct. One can clearly
see a well resolved shoulder at the lower binding energies for Fe $2p_{3/2}$
(marked by the arrow in Fig. 4(d)). This component is missed in the analysis;
however, it might be assigned to another chemical state of Fe atoms in FePSe3,
or the main fitted peak can be due to FexOy in the studied sample [37, 38] and
the small not-assigned peak is due to the FePSe3 phase.

  * 5.

The previous consideration is also valid for the analysis of the Zn $2p$ XPS
spectra. One can clearly see a well resolved shoulder at the lower binding
energies for Zn $2p_{3/2}$ (marked by the arrow in Fig. 4(e)), which is also
is not included in the respective analysis [39, 40].

  * 6.

Analysis of Table S2 of Ref. 22 shows that binding energies for the spin-orbit
split components of the Se $3d$ XPS line change the order for different
compounds, which is, of course, wrong - the binding energy of the $3d_{5/2}$
has to be always smaller. However, this mistake can be related to the typing
error.

##  4 Conclusions

In conclusion, we have demonstrated that the recent work of Gusmão et al. [22]
contains several serious flaws and errors, which are related to the
characterization of the synthesised MPSe3 samples (M = Cd, Cr, Fe, Mn, Sn,
Zn). In this work the structural characterization was performed using XRD and
sample morphology and chemical analysis was studied using the SEM/EDS
combinations. XPS was used for the surface analysis of the studied samples.
Because these data are considered as the main results for the MPSe3 samples
characterization and they are used to support the respective results from the
further studies of the electrochemical performance of the studied objects
(hydrogen and oxygen evolution reactions), we believe that the main
conclusions of the discussed manuscript are not valid and the respective
experimental data, which are criticised here, have to be carefully
reconsidered and reanalyzed. Here we would like to point out that the recently
published data have no limitation period because all presently used
experimental and theoretical approaches were and are available for the last 10
years (at least), as demonstrated by the respective examples.

## Acknowledgement

This contribution was supported by the National Natural Science Foundation of
China (Grant No. 21973059). Y.D. and E.V. thank the support by the Ministry of
Science and Higher Education of the Russian Federation (State assignment in
the field of scientific activity, Southern Federal University, 2020).

## References

## References

  * [1] K. Novoselov, A. Geim, S. Morozov, D. Jiang, M. Katsnelson, I. Grigorieva, S. Dubonos, A. Firsov, Two-dimensional gas of massless Dirac fermions in graphene, Nature 438 (2005) 197–200.
  * [2] Y. Zhang, Y. Tan, H. Stormer, P. Kim, Experimental observation of the quantum Hall effect and Berry’s phase in graphene, Nature 438 (2005) 201–204.
  * [3] Y. Dedkov, E. Voloshina, Graphene growth and properties on metal substrates, J. Phys.: Condens. Matter 27 (2015) 303002.
  * [4] P. Janthon, F. Viñes, S. M. Kozlov, J. Limtrakul, F. Illas, Theoretical assessment of graphene-metal contacts, J. Chem. Phys. 138 (2013) 244701\.
  * [5] R. Roy, R. Thapa, S. Chakrabarty, A. Jha, P. R. Midya, E. M. Kumar, K. K. Chattopadhyay, Role of oxygen functionality on the band structure evolution and conductance of reduced graphene oxide, Chem. Phys. Lett. 677 (2017) 80–86.
  * [6] Y. Dedkov, E. Voloshina, Epitaxial graphene/Ge interfaces: a minireview, Nanoscale, accepted (2020), doi: 10.1039/D0NR00185F.
  * [7] C. Oshima, A. Nagashima, Ultra-thin epitaxial films of graphite and hexagonal boron nitride on solid surfaces, J. Phys.: Condens. Matter 9 (1997) 1–20.
  * [8] S. Manzeli, D. Ovchinnikov, D. Pasquier, O. V. Yazyev, A. Kis, 2D transition metal dichalcogenides, Nature Rev. Mater. 2 (2017) 147.
  * [9] L. Meng, S. Hu, W. Yan, J. Feng, H. Li, X. Yan, Controlled synthesis of large scale continuous monolayer WS2 film by atmospheric pressure chemical vapor deposition, Chem. Phys. Lett. 739 (2020) 136945.
  * [10] A. K. Geim, I. V. Grigorieva, Van der Waals heterostructures, Nature 499 (2014) 419–425.
  * [11] S. Bae, H. Kim, Y. Lee, X. Xu, J.-S. Park, Y. Zheng, J. Balakrishnan, T. Lei, H. R. Kim, Y. I. Song, Y.-J. Kim, K. S. Kim, B. Ozyilmaz, J.-H. Ahn, B. H. Hong, S. Iijima, Roll-to-roll production of 30-inch graphene films for transparent electrodes, Nat. Nanotech. 5 (2010) 574–578.
  * [12] J. Ryu, Y. Kim, D. Won, N. Kim, J. S. Park, E.-K. Lee, D. Cho, S.-P. Cho, S. J. Kim, G. H. Ryu, H.-A.-S. Shin, Z. Lee, B. H. Hong, S. Cho, Fast synthesis of high-performance graphene films by hydrogen-free rapid thermal chemical vapor deposition, ACS Nano 8 (2014) 950–956.
  * [13] F. Schedin, A. K. Geim, S. V. Morozov, E. W. Hill, P. Blake, M. I. Katsnelson, K. S. Novoselov, Detection of individual gas molecules adsorbed on graphene, Nature Mater. 6 (2007) 652–655.
  * [14] T. O. Wehling, M. I. Katsnelson, A. I. Lichtenstein, Adsorbates on graphene: Impurity states and electron scattering, Chem. Phys. Lett. 476 (2009) 125–134.
  * [15] V. Nagarajan, R. Chandiramouli, Adsorption behavior of NH3 and NO2 molecules on stanene and stanane nanosheets - A density functional theory study, Chem. Phys. Lett. 695 (2018) 162–169.
  * [16] N. Wang, M. K. Samani, H. Li, L. Dong, Z. Zhang, P. Su, S. Chen, J. Chen, S. Huang, G. Yuan, X. Xu, B. Li, K. Leifer, L. Ye, J. Liu, Tailoring the thermal and mechanical properties of graphene film by structural engineering, Small 14 (2018) 1801346.
  * [17] X. Wu, V. Varshney, J. Lee, Y. Pang, A. K. Roy, T. Luo, How to characterize thermal transport capability of 2D materials fairly? - Sheet thermal conductance and the choice of thickness, Chem. Phys. Lett. 669 (2017) 233–237.
  * [18] X. Li, T. Cao, Q. Niu, J. Shi, J. Feng, Coupling the valley degree of freedom to antiferromagnetic order, Proc. Natl. Acad. Sci. USA 110 (2013) 3738–3742.
  * [19] X. Li, X. Wu, J. Yang, Half-Metallicity in MnPSe3 exfoliated nanosheet with carrier doping, J. Am. Chem. Soc. 136 (2014) 11065–11069.
  * [20] K.-z. Du, X.-z. Wang, Y. Liu, P. Hu, M. I. B. Utama, C. K. Gan, Q. Xiong, C. Kloc, Weak Van der Waals stacking, Wwide-range band gap, and Raman study on ultrathin layers of metal phosphorus trichalcogenides, ACS Nano 10 (2016) 1738–1743.
  * [21] X. Zhang, X. Zhao, D. Wu, Y. Jing, Z. Zhou, MnPSe3 monolayer: A promising 2D visible-light photohydrolytic catalyst with high carrier mobility, Adv. Sci. 3 (2016) 1600062.
  * [22] R. Gusmão, Z. Sofer, D. Sedmidubský, Š. Huber, M. Pumera, The role of the metal element in layered metal phosphorus triselenides upon their electrochemical sensing and energy applications, ACS Catal. 7 (2017) 8159–8170.
  * [23] M. A. Susner, M. Chyasnavichyus, M. A. McGuire, P. Ganesh, P. Maksymovych, Metal thio- and selenophosphates as multifunctional van der Waals layered materials, Adv. Mater. 29 (2017) 1602852.
  * [24] F. Wang, T. A. Shifa, P. Yu, P. He, Y. Liu, F. Wang, Z. Wang, X. Zhan, X. Lou, F. Xia, J. He, New frontiers on van der Waals layered metal phosphorous trichalcogenides, Adv. Funct. Mater. 28 (2018) 1802151.
  * [25] G. Kresse, J. Furthmuller, Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set, Phys. Rev. B 54 (1996) 11169–11186.
  * [26] G. Kresse, J. Hafner, Norm-conserving and ultrasoft pseudopotentials for first-row and transition elements, J. Phys.: Condens. Matter 6 (1994) 8245–8257.
  * [27] J. Perdew, K. Burke, M. Ernzerhof, Generalized gradient approximation made simple, Phys. Rev. Lett. 77 (1996) 3865–3868.
  * [28] P. E. Blöchl, Projector augmented-wave method, Phys. Rev. B 50 (1994) 17953–17979.
  * [29] P. E. Blöchl, O. Jepsen, O. Andersen, Improved tetrahedron method for Brillouin-zone integrations., Phys. Rev. B 49 (1994) 16223–16233.
  * [30] V. I. Anisimov, A. I. Poteryaev, M. A. Korotin, A. O. Anokhin, G. Kotliar, First-principles calculations of the electronic structure and spectra of strongly correlated systems: dynamical mean-field theory, J. Phys.: Condens. Matter 9 (1997) 7359–7367.
  * [31] S. L. Dudarev, G. A. Botton, S. Y. Savrasov, C. J. Humphreys, A. P. Sutton, Electron-energy-loss spectra and the structural stability of nickel oxide: An LSDA+U study, Phys. Rev. B 57 (1998) 1505–1509.
  * [32] S. Grimme, Semiempirical GGA-type density functional constructed with a long-range dispersion correction, J. Comput. Chem. 27 (2006) 1787–1799.
  * [33] J. Yang, Y. Zhou, Q. Guo, Y. Dedkov, E. Voloshina, Electronic, magnetic and optical properties of MnPX3 (X = S, Se) monolayers with and without chalcogen defects: a first-principles study, RSC Adv. 10 (2020) 851–864.
  * [34] X. Li, Y. Fang, J. Wang, B. Wei, K. Qi, H. Y. Hoh, Q. Hao, T. Sun, Z. Wang, Z. Yin, Y. Zhang, J. Lu, Q. Bao, C. Su, High-yield electrochemical production of large-sized and thinly layered NiPS3 flakes for overall water splitting, Small 15 (2019) 1902427.
  * [35] J. R. Rollo, G. R. Burns, W. T. Robinson, R. J. H. Clark, H. M. Dawes, M. B. Hursthouse, A new polymorph of tetraphosphorus triselenide, $\alpha^{\prime}$-P4Se3: an x-ray, Raman, and XPS study of the normal crystalline phases and a DSC study of the crystalline and the orientationally disordered phases of P4Se3, Inorg. Chem. 29 (1990) 2889–2894.
  * [36] M. C. Biesinger, B. P. Payne, A. P. Grosvenor, L. W. M. Lau, A. R. Gerson, R. S. C. Smart, Resolving surface chemical states in XPS analysis of first row transition metals, oxides and hydroxides: Cr, Mn, Fe, Co and Ni, Appl. Surf. Sci. 257 (2011) 2717–2730.
  * [37] Y. S. Dedkov, M. Fonin, U. Ruediger, C. Laubschat, Graphene-protected iron layer on Ni(111), Appl. Phys. Lett. 93 (2008) 022509.
  * [38] Y. S. Dedkov, A. Generalov, E. N. Voloshina, M. Fonin, Structural and electronic properties of Fe3O4/graphene/Ni(111) junctions, Phys. Status Solidi RRL 5 (2011) 226–228.
  * [39] M. Wang, L. Jiang, E. J. Kim, S. H. Hahn, Electronic structure and optical properties of Zn(OH)2: LDA+U calculations and intense yellow luminescence, RSC Adv. 5 (2015) 87496–87503.
  * [40] D. Kundu, B. D. Adams, V. Duffort, S. H. Vajargah, L. F. Nazar, A high-capacity and long-life aqueous rechargeable zinc battery using a metal oxide intercalation cathode, Nature Energy 1 (2016) 928.

Table 1: Comparison of band gaps obtained in different works using theoretical (DFT) and experimental methods. Compound | PBE+$U$ | PBE+$U$+D2 | HSE06+D2 | Experiment
---|---|---|---|---
|  $U$= 3 eV (Ref. 22) |  $U$= 3 eV |  $U$= 4 eV |  $U$= 5 eV | (Ref. 33) | (Ref. 20)
MnPSe3 |  $1.40$ eV |  $1.60$ eV |  $1.72$ eV |  $1.80$ eV |  $2.50$ eV |  $2.25$ eV
SnPSe3 |  $1.20$ eV |  |  |  |  |
CrPSe3 |  $0.25$ eV |  $0.39$ eV |  $0.70$ eV |  $0.96$ eV |  |
FePSe3 |  $0.50$ eV |  $0.95$ eV |  $1.04$ eV |  $1.08$ eV |  |  $1.20$ eV

Figure 1: Top (a) and side (b) views of a single layer of MPSe3. Spheres of
different size/colour represent ions of different type.

Figure 2: (a) SEM image of MnPSe3 extracted from Ref. 22 (white bar
corresponds to $10\,\mu\mathrm{m}$). The MnPSe3 crystals obtained in our
experiments: (b) general view of the crystal (grid size is
$1\,\mathrm{mm}\times 1\,\mathrm{mm}$), (c) optical microscopy image, (d)
respective XRD patterns.

Figure 3: DOS of bulk MPSe3 (M = Mn, Fe, Cr) in the antiferromagnetic state
calculated in the present work for $U=3$ eV and $U=5$ eV.

Figure 4: (a) The Se $3d$ XPS spectrum measured for WSe2 using monochromatized
Al K$\alpha$ line. Respective spin-orbit split $3d_{3/2}$ and $3d_{5/2}$
emission lines are marked. (b) The Se $3d$ XPS spectrum for $\alpha$-P4Se3
measured using non-monochromatized Mg K$\alpha$ x-ray source [35]. Respective
spin-orbit split emission lines for two kinds of P species are marked. (c-e)
Examples of XPS spectra for MnPSe3, FePSe3, and ZnPSe3 extracted from Ref. 22
and discussed in the present commentary.