# Robust data-driven control design for linear systems subject to input
saturations

A. Seuret, S. Tarbouriech This work has been partially funded by the European
Union - NextGenerationEU. Research by Sophie Tarbouriech is partially funded
by ANR via project HANDY, number ANR-18-CE40-0010.A: Seuret is with the Univ.
de Sevilla, Camino de los Descubrimientos, s/n 41092 Sevilla, Spain. (e-mail:
aseuret@us.es).S. Tarbouriech is with LAAS-CNRS, Université de Toulouse, CNRS,
Toulouse, France. (e-mail: tarbour@laas.fr)..

###### Abstract

This paper deals with the problem of providing a data-driven solution to the
local stabilization of linear systems subject to input saturation. After
presenting a model-based solution to this well-studied problem, a systematic
method to transform model-driven into data-driven LMI conditions is presented.
This technical solution is demonstrated to be equivalent to the recent
advanced results on LMI formulations based on S-procedure or Peterson Lemmas.
However, the advantage of the proposed method relies on its simplicity and its
potential to be applicable to a wide class of problems of stabilization of
(non)linear discrete-time systems. The method is then illustrated on an
academic example.

###### Index Terms:

Linear systems, Saturation, LMI, Data-driven control design

##  I Introduction

The robustness properties for dynamical control systems have been studied in
several works from an analysis or design context [7, 12, 13]. Then, robustness
conditions to face uncertainty or presence of additive disturbances have been
proposed allowing to ensure the stability and a certain level of performance
for the systems under consideration. Most of these conditions are formulated
in the form of linear matrix inequalities (LMIs), due to the powerful
numerical and optimization procedures, as semi-definite programming, which can
be handled. The general common feature of this kind of methods is that they
are model-based, possibly modeling also the presence of uncertainties as norm-
bounded or polytopic uncertainties. In front of complex systems to model, we
need to consider very imprecise or even unknown mathematical models of the
dynamic evolution, including for examples nonlinearities. The consequence is
then that the model-based methods for analysis and controller design may
reveal to be difficult or even impossible to apply. Adapting the tools issued
from the control theory, some work have emerged based on some data
information: see, for example, [1, 2, 9, 19] and the references therein.

Hence, some studies addressing the stability or stabilization criteria from a
data-driven point of view, have been published. In particular several problems
related to the design of state/output feedback controllers for linear systems
have been revisited [6, 18, 19, 5]. In these works considering the case of
exact data experiments for linear time-invariant systems, i.e. without noise
or uncertainties, equivalent formulations between model-based and data-based
criteria both in stability analysis and control design have been exhibited.

Another important feature when dealing with dynamical control systems pertains
the control input saturation due to limitations on the actuators [16]. At our
knowledge, few works in the literature addressed the problem to deal with
constrained control (see, for example, [11] and the references therein) and
more especially with saturating control input. The objective of the current
paper is to bring some preliminary bricks in considering this aspect. To this
end, we consider the design of a state-feedback saturated control law allowing
to ensure the regional (local) asymptotic stability of the origin when the
system is noisy-free and to ensure the convergence to an attractor when the
system is affected by noise. Indeed, the characterization of an inner-
approximation of the basin of attraction of the origin together with an outer-
approximation of the attractor are proposed first by model-based techniques
and then by data-driven techniques. Taking inspiration from the model-based
conditions (see Theorem 1 in Section III), which are written in a friendly-
data form, the data-driven local stabilization is formulated through matrix
inequalities conditions (See Theorem 2 in Section IV-C). The technique
allowing to exhibit the sufficient conditions is based on the use of Lyapunov
arguments, generalized sector-bounded conditions to deal with the saturation
and S-procedure to handle the presence of noise. A first attempt in this
direction was proposed in [15], to revisit the well-know problem of
$\mathcal{L}_{2}$ stabilization but using a data-driven approach. Differently
from this paper, the implicit objectives are here to maximize the inner-
approximation of the basin of attraction of the origin and to minimize the
outer-approximation of the attractor. The main rational behind the matrix
inequalities formulation is due to the friendly form of the model-based
conditions which allows to directly derive the data-driven conditions thanks
to a matrix-constrained relaxation (see Lemma 2 in Section IV-B). Note that
the matrix inequalities conditions are quasi-LMIs in the sense that there is
the product between matrices and scalars.

The paper is organized as follows. Section II presents the system under
consideration and formally state the control problem. In Section III,
sufficient conditions to solve the control problem are formulated as quasi-
LMIs, in the sense that there is a product between a matrix and a scalar.
Section IV deals with the local stabilization from a data-driven point of
view. Section V illustrates the theoretical results and proposes some insights
regarding the influence of the tuning parameters. Comparison between model-
based and data-driven approaches is proposed for different collections of data
to depict the trade-off between the estimate of the basin of attraction and of
the attractor. Some concluding remarks end the paper in Section VI.

Notation. Throughout the paper, $\mathbb{N}$ denotes the set of natural
numbers, $\mathbb{R}$ the real numbers, $\mathbb{R}^{n}$ the $n$-dimensional
Euclidean space, $\mathbb{R}^{n\times m}$ the set of all real $n\times m$
matrices and $\mathbb{S}^{n}$ ($\mathbb{S}^{n}_{+}$) the set of symmetric
(positive definite) matrices in $\mathbb{R}^{n\times n}$. For any $n$ and $m$
in $\mathbb{N}$, matrices $I_{n}$ and $\boldsymbol{0}_{n,m}$
($\boldsymbol{0}_{n}=\boldsymbol{0}_{n,n}$) denote the identity matrix of
$\mathbb{R}^{n\times n}$ and the null matrix of $\mathbb{R}^{n\times m}$,
respectively. When no confusion is possible, the subscripts of these matrices
that precise the dimension, will be omitted. For any matrix $M$ of
$\mathbb{R}^{n\times n}$, the notation $M\succ 0$, ($M\prec 0$) means that $M$
is in $\mathbb{S}^{n}_{+}$. For any matrices $A=A^{\\!\top},B,C=C^{\\!\top}$
of appropriate dimensions, matrix $\left[\begin{smallmatrix}A&B\\\
\ast&C\end{smallmatrix}\right]$ denotes the symmetric matrix
$\left[\begin{smallmatrix}A&B\\\ B^{\\!\top}&C\end{smallmatrix}\right]$. For
any matrix $N\in\mathbb{R}^{n\times m}$, notation $N_{(i)}$, for any
$i=1,\dots,n$, stands for the $i^{th}$ row of $N$. $\|x\|$ denotes the
Euclidean norm of $x$. For a symmetric matrix, $\lambda_{m}(M)$ and
$\lambda_{M}(M)$ denote its minimal and maximal eigenvalues respectively of a
square matrix $M$. For a matrix $M\in\mathbb{S}^{n}$, $M\succ 0$ and any
positive scalar $\alpha\in\mathbb{R}_{>0}$, we denote the ellipsoid
$\mathcal{E}(M,\alpha)=\left\\{x\in\mathbb{R}^{n},~{}x^{\\!\top}Mx\leq\alpha^{-1}\right\\}$.

##  II Problem formulation

###  II-A System data

In this paper, we consider the discrete-time linear system subject to an input
saturation and affected by an external perturbation. Such a system is
described by the following equations

$\displaystyle\left\\{\begin{array}[]{lcl}x^{+}&=&Ax+B\operatorname{sat}(u)+w,\\\
x_{0}&\in&\mathbb{R}^{n_{x}},\end{array}\right.$ (3)

where $x\in\mathbb{R}^{n_{x}}$ is the state vector, which adopts the following
notation $x^{+}=x_{k+1}$ and $x=x_{k}$, $u=u_{k}\in\mathbb{R}^{n_{u}}$ is the
control input. The system is perturbed by the unknown noise signal $w=w_{k}$,
assumed to be bounded leading to the following assumption.

###### Assumption 1

There exists $\lambda>0$ such that the norm of the disturbance verifies
$w^{\\!\top}_{k}w_{k}<\lambda$, for all $k\geq 0$. In other words, the
disturbance belongs to the following set $\Omega_{\lambda}$

$\Omega_{\lambda}:=\\{v\in\mathbb{R}^{n},\quad v^{\\!\top}v\leq\lambda\\}.$
(4)

The model $\mathcal{G}$ is defined from the matrices of the system as
$\mathcal{G}:=\begin{pmatrix}A,B\end{pmatrix}\in\mathbb{R}^{{n_{x}}\times{n_{x}}}\times\mathbb{R}^{{n_{x}}\times{n_{u}}}$.

The saturation function $\operatorname{sat}(u)$ is the classical decentralized
vector-valued saturation map from $R^{{n_{u}}}$ to $R^{{n_{u}}}$, whose the
components are defined by:

$\operatorname{sat}(u_{i})=\operatorname{sign}(u_{i})\min(|u_{i}|,\bar{u}_{i}),\quad\forall
i=1,\dots,{n_{u}}$ (5)

where $u_{i}$ refers to the $i^{th}$ control input and $\bar{u}_{i}$ is the
$i^{th}$ entry of the vector $\bar{u}\in\mathbb{R}^{{n_{u}}}$ that is the
level vector of the saturation.

In this paper, we focus on the design of a static state-feedback control law
of the form

$u=Kx,$ (6)

with $K\in\mathbb{R}^{n_{u}\times n_{x}}$, that regionally (locally)
stabilizes the trajectories of the closed-loop system (3) with (6) in the
absence of perturbation.

###  II-B Preliminaries on generalized sector conditions

Due to the presence of the saturation map in system (3), one has to take care
of the notion of stability we can ensure for the closed-loop system (3)-(6)
[16]. In order to address this problem and therefore to handle the presence of
the saturation, we follow the approach proposed in [8] and [16], which
consists in rewriting (3) as

$x^{+}=(A+BK)x+B\phi(u)+w,$ (7)

where $\phi(u)$ is the dead-zone function defined by

$\phi(u)=\operatorname{sat}(u)-u.$ (8)

Following Remark 7.3 in [16], we use the following Generalized sector
condition lemma to handle the dead-zone (8).

###### Lemma 1

([16]) Consider a matrix $G\in\mathbb{R}^{n_{u}\times n_{x}}$. The following
relation holds

$\phi(u)^{\\!\top}T[\operatorname{sat}(u)+Gx]\leq 0,$ (9)

with any diagonal positive definite matrix
$T\in\mathbb{R}^{{n_{u}}\times{n_{u}}}$ provided that $x\in\mathcal{S}(G)$
where

$\mathcal{S}(G)=\\{x\in\mathbb{R}^{{n_{x}}}:|G_{(i)}x|\leq\bar{u}_{i},\quad\forall
i=1,\dots,n_{u}\\}.$ (10)

###  II-C Control objectives

It is well-known (see [10], [16] and the references therein) that ensuring
global asymptotic stability of the origin for saturated systems (3) (without
perturbation) is in general impossible. In other words, one cannot ensure the
global asymptotic stability of the origin for any initial condition except if
the open loop is not exponentially unstable, and in particular if the open
loop is Hurwitz. That means that the basin of attraction of the origin is not
the whole state space and needs to be determined, that turns out to be a
complex task. In this case, the objective is to characterize an approximation
of the basin of attraction of the origin. Considering level sets built from
Lyapunov functions leads to inner-approximation of the basin of attraction of
the origin. Furthermore, in presence of perturbation satisfying Assumption 1
(that is bounded but not vanishing noise $w$ in the dynamics), the closed-loop
trajectories will not converge to the origin but to some equilibrium points
and therefore one can be interested by characterizing the attractor toward
which will converge the closed-loop trajectories. Using also level sets built
from Lyapunov functions leads to outer-approximation of the real attractor.

Hence, this paper deals with the following problems:

  * (P1)

In the case $w=0$, characterize an approximation of the basin of attraction,
expressed as a level of the Lyapunov function given by $\mathcal{E}(P,1)$,
where $P$ is a symmetric positive definite matrix to be designed. Thus, for
any initial condition belonging to $\mathcal{E}(P,1)$, the closed-loop
trajectories will asymptotically converge to the origin.

  * (P2)

In the case $w\neq 0$, characterize an approximation of the attractor
capturing the closed-loop trajectories expressed as a level set of of the same
Lyapunov function given by $\mathcal{E}(P,\varepsilon)$, where $P$ is a
symmetric positive definite matrix and $\varepsilon>0$ are to be designed.

  * (P3)

In addition to the previous problems, maximize the approximation of the basin
of attraction $\mathcal{E}(P,)$, as well as minimize the approximation of the
attractor $\mathcal{E}(P,\varepsilon)$.

Moreover, as we consider that the model is not assumed to be perfectly known
but only approached via finite set of data experiments, the final objective is
to solve Problems (P1), (P2) and (P3) but from a data-driven point of view.
Hence, we want to revisit the solutions obtained from a model-based approach
to provide a new data-driven practical stabilization criterion, which arises
from matrix manipulations of a suitable model-based preliminary result. This
transformation is made possible thanks to a new lemma, which can be seen as a
particular case of those from the literature, but which is particularly well-
suited for data-driven design.

##  III Model-based local practical stabilization

The following theorem follows the methods presented in Chapter 3 in [16] and
deals with the model-based design of stabilizing controllers (6) for system
(3) under Assumption 1. Therefore, the following theorem states a solution to
solve Problems (P1), (P2) and (P3) through a model-based approach.

###### Theorem 1

Under Assumption 1, and for a given $\mu\in(0,1)$, $\alpha_{1}>0$ and
$\alpha_{2}>0$, assume that there exists

$\begin{array}[]{lcl}\mathcal{D}_{\mathcal{V}}^{1}&:=&\\{\varepsilon,\mu,W,S,Y,Z\\}\\\
&\in&\mathbb{R}_{>0}\\!\times\\!\mathbb{R}_{>0}\\!\times\\!\mathbb{S}_{n_{x}}\\!\times\\!\mathbb{D}_{n_{x}}\\!\times\\!{\mathbb{R}^{{n_{x}}\\!\times\\!{n_{u}}}}\\!\times\\!{\mathbb{R}^{{n_{x}}\\!\times\\!{n_{u}}}}\end{array}$

solution to the following optimization problem

$\begin{array}[]{lcl}\max_{\mathcal{D}_{\mathcal{V}}^{1}}&\alpha_{1}\varepsilon+\alpha_{2}Tr(W)\\\
s.t.&\varepsilon>1,\ \Phi(\mathcal{G})\succ 0,\
\begin{bmatrix}W&Z_{(i)}^{\\!\top}\\\ \ast&\bar{u}_{i}^{2}\\\
\end{bmatrix}\succ 0,\end{array}$ (11)

for all $i=1,\dots,{n_{u}}$, where

$\begin{array}[]{lcl}\Phi(\mathcal{G})&=&\begin{bmatrix}(1-\mu)W&Y^{\\!\top}\\!\\!+\\!Z^{\\!\top}&WA^{\\!\top}\\!\\!+\\!Y^{\\!\top}B^{\\!\top}\\\
\ast&2S&SB^{\\!\top}\\\
\ast&\ast&W-\frac{\lambda\varepsilon}{\mu}I_{n_{x}}\end{bmatrix}.\end{array}$

Then, the control law (6) with $K=YW^{-1}$ ensures that Problems (P1), (P2)
and (P3) are solved, i.e.

  * •

When $w=0$, the ellipsoid $\mathcal{E}(W^{-1},1)$ is an approximation of the
basin of attraction of the origin for the closed-loop system (3)-(6).

  * •

When $w\in\Omega_{\lambda}$ and $w\neq 0$, the solutions to the closed-loop
system (3)-(6) initialized in
$\mathcal{E}(W^{-1},1)\setminus\mathcal{E}(W^{-1},\varepsilon)$ converge to
the attractor $\mathcal{E}(W^{-1},\varepsilon)$.

  * •

Sets $\mathcal{E}(W^{-1},1)$ and $\mathcal{E}(W^{-1},\varepsilon)$ are optimal
with respect to the cost function $\alpha_{1}\varepsilon+\alpha_{2}Tr(W)$,
which aims deriving a compromise between enlarging $\mathcal{E}(W^{-1},1)$ and
minimizing $\mathcal{E}(W^{-1},\varepsilon)$.

###### Proof:

Consider the Lyapunov function given by $V(x)=x^{\\!\top}Px$, for a given
$P\in\mathbb{S}_{+}^{n_{x}}$. In order to address Problems (P1) and (P2), the
objective is to ensure that the forward increment of $V$, i.e. $\Delta
V(x)={x^{+}}^{\\!\top}Px^{+}-x^{\\!\top}Px$ verifies:

$\Delta V(x)\\!\preceq\\!0,\ \forall(x,w)\mbox{
s.t.}\left\\{\\!\\!\begin{array}[]{ll}x\\!\in\\!\mathcal{E}(P,1),&\\!(\mbox{i.e.
}x^{\\!\top}Px\leq 1)\\\
x\\!\notin\\!\mathcal{E}(P,\varepsilon),&\\!(\mbox{i.e.
}x^{\\!\top}Px\geq\varepsilon^{-1})\\\
w\\!\in\\!\Omega_{\lambda},&\\!(\mbox{i.e.
}w^{\\!\top}w\leq\lambda)\end{array}\right.$ (12)

which is implied, thanks to the use of S-procedure [4], by the following
condition

$\begin{array}[]{lcl}\mathcal{L}(x,w)&\\!\\!\\!\\!\\!=&\\!\\!\\!\\!\Delta
V(x)\\!+\\!\bar{\mu}_{1}(x^{\\!\top}Px\\!-\\!\varepsilon^{-1})\\!+\\!\bar{\mu}_{2}(\lambda\\!-\\!w^{\\!\top}w)\\\
&&\\!\\!\\!+\bar{\mu}_{3}(1-x^{\\!\top}Px)\preceq 0\end{array}$ (13)

with $\bar{\mu}_{i}\geq 0$, $i=1,2,3$. Let us introduce the following
augmented vector $\xi=(x,\phi(Kx),w)$, then $\Delta V(x)$ reads:

$\Delta V(x)=\xi^{\\!\top}\left(\begin{bmatrix}(A+BK)^{\\!\top}\\\
B^{\\!\top}\\\
I_{n_{x}}\end{bmatrix}\\!\\!\\!P\begin{bmatrix}(A+BK)^{\\!\top}\\\
B^{\\!\top}\\\ I_{n_{x}}\end{bmatrix}^{\\!\top}-\begin{bmatrix}P&0&0\\\
0&0&0\\\ 0&0&0\end{bmatrix}\right)\xi\\\ $ (14)

The satisfaction of the last inequalities in (11) means that the ellipsoid
$\mathcal{E}(W,1)$ is included in the set $\mathcal{S}(G)$. From Lemma 1, for
$x\in\mathcal{E}(W,1)$ inequality (9) holds and reads using $\xi$ as follows:

$\xi^{\\!\top}\begin{bmatrix}0&(K+G)^{\\!\top}T&0\\\ T(K+G)&2T&0\\\
0&0&0\end{bmatrix}\xi\leq 0.$ (15)

Then, merging (15) and (13) yields

$\mathcal{L}_{1}(x,w):=\mathcal{L}(x,w)-\xi^{\\!\top}\begin{bmatrix}0&(K+G)^{\\!\top}T&0\\\
T(K+G)&2T&0\\\ 0&0&0\end{bmatrix}\xi\geq\mathcal{L}(x,w).$

Therefore, having $\mathcal{L}_{1}(x,w)$ negative implies that
$\mathcal{L}(x,w)$ is also negative for all $(x,w)$. Let us rewrite the
expression of $\mathcal{L}_{1}(x,w)$ in a more compact form.

$\mathcal{L}_{1}(x,w):=-\xi^{\\!\top}\Phi_{1}(\mathcal{G})\xi-\bar{\Theta}_{\lambda}$

where

$\begin{array}[]{rcl}\Phi_{1}(\mathcal{G})&\\!\\!\\!\\!=&\\!\\!\\!\\!\\!\begin{bmatrix}\\!\\!(1\\!-\\!\bar{\mu}_{1}\\!+\\!\bar{\mu}_{3})P&(K\\!+\\!G)^{\\!\top}T&\\!0\\\
T(K\\!+\\!G)&2T&\\!0\\\
0&0&\\!\\!\\!\\!\bar{\mu}_{2}I_{n_{x}}\end{bmatrix}\\!-\\!\begin{bmatrix}\\!(A\\!+\\!BK)^{\\!\top}\\!\\!\\\
B^{\\!\top}\\\
I_{n_{x}}\end{bmatrix}\\!P\begin{bmatrix}\\!(A\\!+\\!BK)^{\\!\top}\\!\\!\\\
B^{\\!\top}\\\ I_{n_{x}}\end{bmatrix}^{\\!\\!\top}\\\
\Theta_{\lambda}&\\!\\!\\!\\!=&\\!\\!\\!\\!\bar{\mu}_{1}\varepsilon^{-1}-\bar{\mu}_{2}\lambda-\bar{\mu}_{3}\end{array}$

In the next developments, we will show that the existence of a solution to
condition $\Phi(\mathcal{G})\succ 0$ ensures the existence of
$\bar{\mu}_{1}>0$, $\bar{\mu}_{2}>0$ and $\bar{\mu}_{3}>0$ such that
$\mathcal{L}_{1}(x,w)\leq 0$.

Applying the Schur complement to $\Phi_{1}(\mathcal{G})$, pre- and post-
multiplying it by $\operatorname{diag}(W,S,I_{2n_{x}})$, with $W=P^{-1}$ and
$S=T^{-1}$, and selecting $Y=KW$ and $Z=GW$ yield

$\begin{bmatrix}(1\\!-\\!\bar{\mu}_{1}\\!+\\!\bar{\mu}_{3})W&Y^{\\!\top}\\!\\!+\\!Z^{\\!\top}&0&WA^{\\!\top}\\!+\\!Y^{\\!\top}B^{\\!\top}\\\
\ast&2S&0&SB^{\\!\top}\\\ \ast&\ast&\bar{\mu}_{2}I_{n_{x}}&I_{n_{x}}\\\
\ast&\ast&\ast&W\end{bmatrix}\succ 0,$

which is equivalent to

$\Phi_{2}(\mathcal{G}):=\begin{bmatrix}(1\\!-\\!\bar{\mu}_{1}\\!+\\!\bar{\mu}_{3})W&Y^{\\!\top}\\!\\!\\!+\\!\\!Z^{\\!\top}&WA^{\\!\top}\\!+\\!Y^{\\!\top}\\!B^{\\!\top}\\\
\ast&2S&SB^{\\!\top}\\\
\ast&\ast&W-\bar{\mu}_{2}^{-1}I_{n_{x}}\end{bmatrix}\succ 0$ (16)

We have shown so far that having $\mathcal{L}_{1}\leq 0$ is equivalent to
$\Phi_{2}(\mathcal{G})\succ 0$ and $\Theta_{\lambda}>0$. Enforcing the
introduction of $\Phi(\mathcal{G})$, we get the following expression

$\Phi_{2}(\mathcal{G})\\!=\\!\Phi(\mathcal{G})\\!+\\!\begin{bmatrix}(\mu\\!-\\!\bar{\mu}_{1}\\!+\\!\bar{\mu}_{3})W&0&0\\\
\ast&0&0\\\
\ast&\ast&\left(\frac{\lambda\varepsilon}{\mu}\\!-\\!\frac{1}{\bar{\mu}_{2}}\right)I_{n_{x}}\end{bmatrix}.$
(17)

Then, consider the following selection

$\begin{array}[]{lcl}\bar{\mu}_{1}:=\mu\in(0,1),\quad\bar{\mu}_{2}:=\frac{\mu}{\lambda\varepsilon}\left(1-\frac{\varepsilon\gamma}{\mu}\right),\quad\bar{\mu}_{3}<\gamma,\end{array}$

for sufficiently small $\gamma>0$ such that $\bar{\mu}_{2}>0$. This selection
ensures that $\Theta_{\lambda}=\gamma-\bar{\mu}_{3}>0$, which is required.
Using this selection for $\bar{\mu}_{2}$, we note that

$\frac{\lambda\varepsilon}{\mu}\\!-\\!\frac{1}{\bar{\mu}_{2}}\\!=\\!\frac{\lambda\varepsilon}{\mu}\left(1\\!-\\!\frac{1}{1\\!-\\!\gamma\varepsilon/\mu}\right)\\!=\\!\frac{-\gamma\lambda\varepsilon^{2}}{\mu^{2}\\!-\\!\gamma\varepsilon\mu}>\frac{-\bar{\mu}_{3}\lambda\varepsilon^{2}}{\mu^{2}\\!-\\!\gamma\varepsilon\mu},$

so that the following inequality holds

$\begin{array}[]{lcl}\Phi_{2}(\mathcal{G})\succeq\Phi(\mathcal{G})+\bar{\mu}_{3}\begin{bmatrix}W&0&0\\\
\ast&0&0\\\
\ast&\ast&-\frac{\lambda\varepsilon^{2}}{\mu^{2}\\!-\\!\gamma\varepsilon\mu}I_{n_{x}}\end{bmatrix}.\end{array}$
(18)

Hence, if condition $\Phi(\mathcal{G})\succ 0$ holds, then there exists a
sufficiently small $\bar{\mu}_{3}$ such that $\Phi_{2}(\mathcal{G})\succ 0$.
Therefore, condition $\Phi(\mathcal{G})\succ 0$ ensures that there exist
$\bar{\mu}_{1}>0,\bar{\mu}_{2}>0,\bar{\mu}_{3}>0$ such that
$\mathcal{L}_{1}(x,w)$ is negative, which ensures (12).

Therefore, the satisfaction of relation $\Phi(\mathcal{G})\succ 0$ ensures the
local stability of the closed-loop system with the control gain $K=YW^{-1}$,
for any initial condition in $x\in\mathcal{E}(W^{-1},1)$.

Since $\varepsilon>1$, one gets $\mathcal{E}(W^{-1},\varepsilon)\subseteq$
$\mathcal{E}(W^{-1},1)$. In order to prove that
$\mathcal{E}(W^{-1},\varepsilon)$ is an attractor for the closed-loop system,
it remains to demonstrate that $\mathcal{E}(W^{-1},\varepsilon)$ is invariant,
i.e.

$x\in\mathcal{E}(W^{-1},\varepsilon)\quad\Rightarrow\quad
x^{+}\in\mathcal{E}(W^{-1},\varepsilon).$

To do so, the satisfaction of $\Phi(\mathcal{G})\succ 0$ ensures that

$\begin{array}[]{lcl}V(x^{+})&\\!\\!\\!\\!\\!=&\displaystyle\\!\\!V(x)+\underbrace{\mathcal{L}_{1}(x,w)}_{\leq
0}+\underbrace{\xi^{\\!\top}\begin{bmatrix}0&(K+G)^{\\!\top}T&0\\\
\ast&2T&0\\\ \ast&\ast&0\end{bmatrix}\xi}_{\leq 0}\\\
&\\!\\!\\!\\!\\!-&\\!\\!\\!\\!\\!\bar{\mu}_{1}(V(x)\\!-\\!\varepsilon^{-1})\\!\underbrace{-\bar{\mu}_{2}(\lambda\\!-\\!w^{\\!\top}w)}_{\leq
0}\\!-\bar{\mu}_{3}(1\\!-\\!V(x))\\\
&\\!\\!\\!\\!\\!\leq&\\!\\!\displaystyle(1-\bar{\mu}_{1}+\bar{\mu}_{3})V(x)+\\!\bar{\mu}_{1}\varepsilon^{-1}-\\!\bar{\mu}_{3}\\\
\end{array}$

Since $\bar{\mu}_{1}=\mu$ is in $(0,1)$, there exists a sufficiently small
$\bar{\mu}_{3}$ such that $1-\bar{\mu}_{1}+\bar{\mu}_{3}>0$ and, consequently,
as $V(x)\leq\varepsilon^{-1}$, we have

$\begin{array}[]{lcl}V(x^{+})\\!\leq\\!(1\\!-\\!\bar{\mu}_{1}\\!+\\!\bar{\mu}_{3})\varepsilon^{-1}\\!+\\!\bar{\mu}_{1}\varepsilon^{-1}\\!-\\!\bar{\mu}_{3}\\!=\\!\varepsilon^{-1}\\!-\\!\bar{\mu}_{3}(1\\!-\\!\varepsilon^{-1}),\end{array}$

which implies that $V(x^{+})\\!\leq\\!\varepsilon^{-1}$, since
$\varepsilon>1$.

The last step of the proof addresses the particular case when $w=0$ or
equivalently $\lambda=0$. In this situation, $\varepsilon$ disappears from the
condition and can be selected as large as possible, so that
$\mathcal{E}(W^{-1},\varepsilon)$ shrinks to $\\{0\\}$. Therefore, the
solution to the system, initialized in $\mathcal{E}(W^{-1},1)$ converges to
the origin. ∎

###### Remark 1

Note that $\Phi(\mathcal{G})\succ 0$ is not an LMI due to the terms $\mu W$
and $\lambda\varepsilon\mu^{-1}$. Nevertheless, once the parameter $\mu$ is
fixed a priori, the problem becomes an LMI and can easily be solved. A
gridding on the parameter $\mu$ should be included to the LMI solver.

##  IV Local data-driven control design

The objective of this section is to transform the previous model-based theorem
into a data-based one. The basic idea is to use the inherent robustness
property of the previous theorem with respect to $(A,B)$, to derive such a
transformation. After presenting the main features and assumptions on the data
experiments, a new lemma is presented to achieve this transformation in a
systematic manner. An application of this lemma to the problem in hand finally
allows us to exhibit the main result of this paper.

###  IV-A Data collections and assumption

Unlike the usual situation where the system model is available, our objective
is here to provide a data-driven-based design result. To do this, let us first
specify the notion of data, we are considering here. We define the following
data matrices that collect the available measurement for the control design.

$\begin{array}[]{lclcllcr}\mathcal{X}^{+}&\\!\\!:=&\\![&\\!x^{+}_{1}&x^{+}_{2}&\dots&x^{+}_{p}&\\!\\!],\\\
\mathcal{X}&\\!\\!:=&\\![&\\!x_{1}&x_{2}&\dots&x_{p}&\\!\\!],\\\
\mathcal{U}&\\!\\!:=&\\![&\\!\operatorname{sat}(u_{1})&\operatorname{sat}(u_{2})&\dots&\operatorname{sat}(u_{p})&\\!\\!],\\\
\omega&\\!\\!:=&\\![&\\!w_{1}&w_{2}&\dots&w_{p}&\\!\\!].\\\ \end{array}$ (19)

We define the data collection $\mathcal{D}$ as follows

$\mathcal{D}:=(\mathcal{X}^{+},\mathcal{X},\mathcal{U})\in\mathbb{R}^{n_{x}\times
p}\times\mathbb{R}^{n_{x}\times p}\times\mathbb{R}^{n_{u}\times p}.$ (20)

According the systems dynamics, these data verify

$\mathcal{X}^{+}=A\mathcal{X}+B\mathcal{U}+\omega\quad\in\mathbb{R}^{n_{x}\times
p}.$ (21)

for some matrices $(A,B)$. The noise samples $w(0),w(1),$ $\dots,w(p)$
collected in the matrix $\omega$ are unknown. However, the following
assumption is made on the “energy” of the noise, following the presentation of
[3].

###### Assumption 2

Assume that there exists a known matrix $\Delta_{\omega}$ in
$\mathbb{S}^{n_{x}}_{+}$ such that the noise samples matrix $\omega$ verify

$\omega\in\Omega:=\\{v\in\mathbb{R}^{n_{x}\times p},\ vv^{\\!\top}\leq
p\lambda\Delta_{\omega}\\}.$ (22)

Observe in this assumption that matrix $\Delta_{\omega}$ is independent of the
noise amplitude $\lambda$ and of the number of experiments $p$. This
assumption can be related to the co-variance of the noise $\omega$ as
mentioned in [3] or [18].

The authors of [3] have interestingly presented an alternative formulation of
$\Omega$ using the data collection $\mathcal{D}$. This method that was also
developed in many other papers such as, for example, [18], [19], consists in
redefining the following set

$\mathcal{C}\\!:=\\!\\{[A\ B],\
\mathcal{X}^{+}\\!\\!=\\!A\mathcal{X}\\!+\\!B\mathcal{U}\\!+\\!\omega,\
\forall\omega\in\Omega\\}\subset\mathbb{R}^{n_{x}\times(n_{x}+n_{u})}.$ (23)

Thus, an alternative representation of $\mathcal{C}$ can be derived as noticed
in [3], as follows

$\mathcal{C}\\!\\!:=\\!\left\\{\\![A\ B],\begin{bmatrix}I\\\ A^{\\!\top}\\\
B^{\\!\top}\end{bmatrix}^{\\!\\!\top}\\!\\!\begin{bmatrix}\mathcal{X}^{+}{\mathcal{X}^{+}}^{\\!\top}\\!\\!-p\lambda\Delta_{\omega}&-\mathcal{X}^{+}\mathcal{X}^{\\!\top}&-\mathcal{X}^{+}\mathcal{U}^{\\!\top}\\\
\ast&\mathcal{X}\mathcal{X}^{\\!\top}&\mathcal{X}\mathcal{U}^{\\!\top}\\\
\ast&\ast&\mathcal{U}\mathcal{U}^{\\!\top}\\\
\end{bmatrix}\\!\begin{bmatrix}I\\\ A^{\\!\top}\\\
B^{\\!\top}\end{bmatrix}\\!\\!\preceq\\!0\\!\right\\}$ (24)

We classically consider that the data are informative as defined below.

###### Definition 1

The data collection $\mathcal{D}$ in (19) is said informative if matrix
$\begin{bmatrix}\mathcal{X}\\\ \mathcal{U}\\\
\end{bmatrix}\begin{bmatrix}\mathcal{X}\\\ \mathcal{U}\\\
\end{bmatrix}^{\\!\top}$ is non singular.

###  IV-B Matrix-constrained Relaxation

The following lemma, which is the main brick of the paper, presents a generic
method to transform a problem of a particular matrix inequality which depends
on parameters verifying a quadratic constraint into a formulation that is
independent of these parameters. It is stated below.

###### Lemma 2

For given positive integers $n_{1},n_{2},n_{3}$, consider matrices
$(\mathcal{M}_{1},\mathcal{M}_{2},\mathcal{M}_{3})$ in
$\mathbb{S}^{n_{1}}_{+}\times\mathbb{R}^{n_{1}\times
n_{2}}\times\mathbb{S}^{n_{3}}_{+}$ and
$(\mathcal{N}_{1},\mathcal{N}_{2},\mathcal{N}_{3})$ in $\
S^{n_{3}}\times\mathbb{R}^{n_{3}\times n_{2}}\times\mathbb{S}^{n_{3}}_{+}$,
i.e. with $\mathcal{N}_{3}\succ 0$.

Then, the following statements are equivalent

  * (i)

Inequality

$\mathcal{M}(\mathcal{A})=\begin{bmatrix}\mathcal{M}_{1}&\mathcal{M}_{2}\mathcal{A}\\\
\ast&\mathcal{M}_{3}\end{bmatrix}\succ
0,\quad\forall\mathcal{A}\in\Sigma_{\mathcal{N}},$ (25)

holds true where $\Sigma_{\mathcal{N}}$ represents the set of allowable
uncertain matrices $\mathcal{A}$ characterized by a quadratic constraint
defined as follows

$\\!\\!\\!\\!\Sigma_{\mathcal{N}}\\!:=\\!\left\\{\mathcal{A}\in\mathbb{R}^{n_{2}\times
n_{3}},\left[\begin{matrix}I_{n_{3}}\\\ \mathcal{A}\\\
\end{matrix}\right]^{\\!\top}\\!\\!\underbrace{\begin{bmatrix}\mathcal{N}_{1}&\\!\\!\mathcal{N}_{2}\\\
\ast&\\!\\!\mathcal{N}_{3}\\\
\end{bmatrix}}_{\mathcal{N}}\\!\left[\begin{matrix}I_{n_{3}}\\\ \mathcal{A}\\\
\end{matrix}\right]\\!\preceq\\!0.\right\\}$ (26)

  * (ii)

There exists $\eta>0$ such that

$\begin{bmatrix}\mathcal{M}_{1}&\boldsymbol{0}_{n_{1},{n_{3}}}&\mathcal{M}_{2}\\\
\ast&\mathcal{M}_{3}+\eta\mathcal{N}_{1}&\eta\mathcal{N}_{2}\\\
\ast&\ast&\eta\mathcal{N}_{3}\end{bmatrix}\succ 0.$ (27)

The proof of this lemma is postponed to the appendix.

Lemma 2 provides an alternative formulation in robust analysis for uncertain
matrices subject to quadratic constraints of the form (26) compared to the one
presented in [18, 19]. For the sake of consistency, the S-Lemma provided in
[18] is recalled in the following lemma.

###### Lemma 3

[18, Th.9] Let
$\mathcal{M}_{s},\mathcal{N}_{s}\in\mathbb{R}^{(n_{s}+m_{s})\times(n_{s}+m_{s})}$
be symmetric matrices and assume that inequality $\begin{bmatrix}I_{n_{s}}\\\
\mathcal{A}_{s}^{\\!\top}\end{bmatrix}^{\\!\top}\mathcal{N}_{s}\begin{bmatrix}I_{n_{s}}\\\
\mathcal{A}_{s}^{\\!\top}\end{bmatrix}\succeq 0$ for at least one matrix
$\mathcal{A}_{s}\in\mathbb{R}^{n_{s}\times m_{s}}$. Then the next statements
are equivalent:

  1. (i)

$\begin{bmatrix}I_{n_{s}}\\\
\mathcal{A}_{s}^{\\!\top}\end{bmatrix}^{\\!\top}\mathcal{M}_{s}\begin{bmatrix}I_{n_{s}}\\\
\mathcal{A}_{s}^{\\!\top}\end{bmatrix}\succ 0,\
\forall\mathcal{A}_{s}\in\Sigma_{\mathcal{N}_{s}}$,
where set $\Sigma_{\mathcal{N}_{s}}$ has the same definition as in (26) but
replacing $\mathcal{A}$ and $\mathcal{N}$ by $\mathcal{A}_{s}$ and
$\mathcal{N}_{s}$, respectively.

  2. (ii)

There exists $\eta>0$ such that $\mathcal{M}_{s}-\eta\mathcal{N}_{s}\succ 0$.

Note that both lemmas address the problem of the satisfaction of an inequality
subject to uncertain matrices characterized by a quadratic constraint. The
main interest of both lemmas is to derive equivalent inequalities that are
independent of the uncertain matrix $\mathcal{A}$ (or $\mathcal{A}_{s}$).
Finally, both lemmas can be seen as application of the usual manipulations on
LMI such as Schur Complement, Finsler’s lemma and S-procedure. Apart from
presenting these similarities, both lemmas have substantial differences.
First, Lemma 3 requires that matrix $\mathcal{M}_{s}$ has the same size as
$\mathcal{N}_{s}$, the matrix that characterizes the quadratic constraint on
$\mathcal{A}$. Lemma 2 is more flexible in this sense, as there is no
relationship between matrices $\mathcal{M}_{1}$ and $\mathcal{N}$, which are
independent. This flexibility has the benefit of reducing the initial
manipulations to derive, from usual stability or control problems, the
appropriate expression of $\mathcal{M}_{i}$ and $\mathcal{N}_{i}$, for
$i=1,2,3$, to fit the framework of Lemma 3. In fact, the relationship between
both lemmas can be seen by selecting $n_{s}=n_{1}+n_{3}$, $m_{s}=n_{2}$,
$\mathcal{A}_{s}^{\\!\top}=\begin{bmatrix}0_{n_{2}\times
n_{3}}&\mathcal{A}^{\\!\top}\end{bmatrix}$ and

$\begin{array}[]{c}\mathcal{M}_{s}=\begin{bmatrix}\mathcal{M}_{1}&0&\mathcal{M}_{2}\\\
\ast&\mathcal{M}_{3}&0\\\ \ast&\ast&0\\\
\end{bmatrix},\quad\mathcal{N}_{s}=-\begin{bmatrix}0&0&0\\\ \ast&\ \
\mathcal{N}_{1}&\mathcal{N}_{2}\\\ \ast&\ast&\mathcal{N}_{3}\\\
\end{bmatrix}.\end{array}$

From this selection, it is clear that item (ii) of Lemma 3 is equivalent to
item (ii) of Lemma 2, showing that Lemma 2 is a particular case of Lemma 3.
That being said, the main advantage of Lemma 2 is that the structure of
$\mathcal{M}(\mathcal{A})$ arises in many LMI problems in control as it will
be shown in the next section. Therefore it avoids one from enforcing an
initial LMI problem to fit with the structure of Lemma 3, which is not an easy
task in general.

###### Remark 2

A similar lemma was already presented in [14], where block
$\mathcal{M}_{2}\mathcal{A}$ is extended to
$\bar{\mathcal{M}}_{2}+\mathcal{M}_{2}\mathcal{A}$, where
$\bar{\mathcal{M}}_{2}$ is a term that is independent of the uncertain matrix.
In addition, a deeper discussion on the similarities with the existing lemmas
from the literature has been proposed therein.

###  IV-C Data-driven local stabilization of saturated systems

This section provides a new contribution on the data-based design of
stabilizing control law for linear systems subject to input saturation. The
method is highly inspired from [18] but has been reformulated to get a simpler
and user-friendly formulation. In this paragraph, we will demonstrate how
Lemma 2 can be easily applied to the stabilization problem of saturated
systems. The robust local stabilization of system (3) can be formalized in the
following theorem.

###### Theorem 2

Under Assumptions 1 and 2, i.e. for a given matrix $\Delta_{w}>0$, and for
given $\mu\in(0,1)$, $\alpha_{1}>0$ and $\alpha_{2}>0$, assume that there
exist

$\begin{array}[]{lcl}\mathcal{D}_{\mathcal{V}}^{2}&:=&\\{\varepsilon,\eta,W,S,Y,Z\\}\\\
&\in&\mathbb{R}_{>0}\\!\times\\!\mathbb{R}_{>0}\\!\times\\!\mathbb{S}_{n_{x}}\\!\times\\!\mathbb{D}_{n_{x}}\\!\times\\!{\mathbb{R}^{{n_{x}}\\!\times\\!{n_{u}}}}\\!\times\\!{\mathbb{R}^{{n_{x}}\\!\times\\!{n_{u}}}}\end{array}$

that are solution to the following optimization problem

$\begin{array}[]{lcl}\max_{\mathcal{D}_{\mathcal{V}}^{2}}&\alpha_{1}\varepsilon+\alpha_{2}Tr(W)\\\
s.t.&\varepsilon>1,\ \eta>0,\ \Psi(\mathcal{D})\succ 0,\
\begin{bmatrix}W&Z_{(i)}^{\\!\top}\\\ \ast&\bar{u}_{i}^{2}\\\
\end{bmatrix}\succ 0,\end{array}$ (28)

for all $i=1,\dots,{n_{u}}$, where

$\begin{array}[]{rcl}\Psi(\mathcal{D})&\\!\\!\\!\\!=&\\!\\!\\!\\!\begin{bmatrix}(1\\!-\\!\mu)W&\\!\\!\\!Y^{\\!\top}\\!\\!+\\!Z^{\\!\top}&\\!\\!\\!\\!\\!0&\\!\\!\\!\\!W&\\!\\!\\!\\!Y^{\\!\top}\\\
\ast&\\!\\!\\!2S&\\!\\!\\!\\!\\!0&\\!\\!\\!\\!0&\\!\\!\\!\\!\\!S\\\
\ast&\\!\\!\\!\ast&\\!\\!\\!\\!\\!{\Psi}_{3}&\\!\\!-\eta\mathcal{X}^{+}\mathcal{X}^{\\!\top}&\\!\\!-\eta\mathcal{X}^{+}\mathcal{U}^{\\!\top}\\\
\ast&\\!\\!\\!\ast&\\!\\!\\!\\!\\!\ast&\\!\\!\eta\mathcal{X}\mathcal{X}^{\\!\top}&\\!\\!\eta\mathcal{X}\mathcal{U}^{\\!\top}\\\
\ast&\\!\\!\\!\ast&\\!\\!\\!\\!\ast&\\!\\!\ast&\\!\\!\eta\mathcal{U}\mathcal{U}^{\\!\top}\end{bmatrix}\\\
{\Psi}_{3}&\\!\\!\\!\\!\\!=&\\!W-\frac{\lambda\varepsilon}{\mu}I_{n_{x}}+\eta\mathcal{X}^{+}{\mathcal{X}^{+}}^{\\!\top}-\eta
p\lambda\Delta_{\omega}\\\ \end{array}$ (29)

Then, the control law (6) with $K=YW^{-1}$ ensures that Problems (P1) and (P2)
are solved, or equivalently, the following statements hold:

  * •

When $w=0$, the ellipsoid $\mathcal{E}(W^{-1},1)$ is an approximation of the
basin of attraction of the origin for the closed-loop system (3)-(6);

  * •

When $w\in\Omega_{\lambda}$ and $w\neq 0$, the solutions to the closed-loop
system (3)-(6) initialized in
$\mathcal{E}(W^{-1},1)\setminus\mathcal{E}(W^{-1},\varepsilon)$ converge to
the attractor $\mathcal{E}(W^{-1},\varepsilon)$.

  * •

The estimation of the basin of attraction and of the attractor are optimal
with respect to the cost function $\alpha_{1}\varepsilon+\alpha_{2}Tr(W)$,
leading to a compromise between enlarging $\mathcal{E}(W^{-1},1)$ and
minimizing $\mathcal{E}(W^{-1},\varepsilon)$.

###### Proof:

Let us first note that, in Theorem 1, conditions $\varepsilon>1$ and
$\begin{bmatrix}W&Z_{(i)}^{\\!\top}\\\ \ast&\bar{u}_{i}\\\ \end{bmatrix}\succ
0$ do not involve the matrices of the system and thus remains unchanged. On
the other side, condition $\Phi(\mathcal{G})\succ 0$ has to be modified in
order to remove the system’s matrices. To do so, let us note that
$\Phi(\mathcal{G})$ is rewritten as follows

$\Phi(\mathcal{G})=\begin{bmatrix}\begin{bmatrix}(1-\mu)W&X^{\\!\top}+Y^{\\!\top}\\\
\ast&2S\end{bmatrix}&\begin{bmatrix}W&Y^{\\!\top}\\\
0&S\end{bmatrix}\left[\begin{matrix}A^{\\!\top}\\\ B^{\\!\top}\\\
\end{matrix}\right]\\\
\ast&W-\frac{\lambda\varepsilon}{\mu}I_{n_{x}}\end{bmatrix},$ (30)

which has the same structure as the one of matrix $\mathcal{M}(\mathcal{A})$
introduced in Lemma 2. Indeed, selecting $n_{1}=n_{2}={n_{x}}+{n_{u}}$,
$n_{3}={n_{x}}$ and

$\begin{array}[]{l}\mathcal{M}_{1}\\!=\\!\begin{bmatrix}\\!(1\\!-\\!\mu)W&\\!\\!\\!\\!X^{\\!\top}+Y^{\\!\top}\\\
\ast&\\!\\!\\!\\!2S\end{bmatrix}\\!\\!,\
\mathcal{M}_{2}\\!=\\!\begin{bmatrix}W&\\!\\!\\!\\!Y^{\\!\top}\\\
0&\\!\\!\\!\\!S\end{bmatrix}\\!,\\\
\mathcal{M}_{3}\\!=\\!W\\!\\!-\\!\frac{\lambda\varepsilon}{\mu}I_{n_{x}},\\\
\mathcal{N}_{1}=\mathcal{X}^{+}{\mathcal{X}^{+}}^{\\!\top}\\!\\!-\\!p\lambda\Delta_{\omega},\
\mathcal{N}_{2}\\!=\\!-\mathcal{X}^{+}\begin{bmatrix}\mathcal{X}\\\
\mathcal{U}\\\ \end{bmatrix}^{\\!\top}\\!\\!\\!,\
\mathcal{N}_{3}\\!=\\!\begin{bmatrix}\mathcal{X}\\\ \mathcal{U}\\\
\end{bmatrix}\begin{bmatrix}\mathcal{X}\\\ \mathcal{U}\\\
\end{bmatrix}^{\\!\top}\\!\\!\\!.\end{array}$

with the uncertainty matrix

$\mathcal{A}=\left[\begin{matrix}A^{\\!\top}\\\ B^{\\!\top}\\\
\end{matrix}\right]\in\mathcal{C}\subset\mathbb{R}^{(n_{x}+n_{u})\times
n_{x}}.$

Altogether, the problem can be expressed as in (25) with this set of matrices.

Then, the problem resumes to the satisfaction of condition
$\Phi(\mathcal{G})\succ 0$ for all matrices $[A\
B]^{\\!\top}=\mathcal{A}\subset\mathcal{C}$. Note that $\mathcal{N}_{3}\succ
0$ is required in Lemma 2, which refers to the informativity of the data and
is also required for the satisfaction of condition $\Psi(\mathcal{D})\succ 0$.
We are thus in position to apply Lemma 2, which states that having
$\Phi(\mathcal{G})\succ 0$ for all $\mathcal{A}\in\mathcal{C}$ is equivalent
to the existence of $\eta>0$ such that $\Psi(\mathcal{D})\succ 0$. ∎

###### Remark 3

Note that the informativity of the data is a necessary condition for the
solvability of the LMI optimization problem. In other words, if matrix
$\begin{bmatrix}\mathcal{X}\\\ \mathcal{U}\\\
\end{bmatrix}\begin{bmatrix}\mathcal{X}\\\ \mathcal{U}\\\
\end{bmatrix}^{\\!\top}$ is singular, then the LMI problem cannot be solved.

###### Remark 4

Note that, contrary to [1] or [18], our method does not require to study the
dual system involving $(A+BK)^{\\!\top}$. First, this trick is not possible
for nonlinear systems. Second, our method directly treats the initial problem
without using this artifice of calculus to avoid a technical problem.

###### Remark 5

Differently from [3], the additional variable $\eta$ cannot be “absorbed” by
the other decision variables. This due to the problem of local stabilization,
for which the scalability of the Lyapunov matrix $W$ is not permitted since it
appears multiplied by the scalar $\mu$.

###### Remark 6

It is worth noting that the resulting LMI condition is very similar to the
ones presented in [3] based on the application of the Pertersen’s lemma for
the case of linear systems. Indeed, selecting $S=0$, $Y=Z=0$ and
$\varepsilon=0$, the same LMI condition as the one of Th.2 in [3] is
retrieved. This makes us think that Lemma 2 can be related to the Petersen’s
lemma as well as we have showed the link with the S-Procedure of [18]. Again,
the advantage of Lemma 2 over both these lemmas is that the structure of many
LMI problems in control fits with the structure of Lemma 2, so that no
additional manipulation is required.

##  V Numerical application

Consider the discrete-time systems (7) borrowed from [17, 20], with the
following matrices

$A=\begin{bmatrix}0.8&0.5\\\ -0.4&1.2\end{bmatrix},\quad B=\begin{bmatrix}0\\\
1\end{bmatrix},\quad\bar{u}_{1}=5.$ (31)

###  V-A Model-based solution

This section aims at illustrating the model-based stabilization result of
Theorem 1. Solving the optimization problem (11) with $\alpha_{1}=1$ and
$\alpha_{2}=10^{-3}$ for $\lambda=0.01,0.05$ and $0.1$ and $\mu=0.08,0.3$ and
$0.6$. Figure 1 shows the estimation of the basin of attraction (dashed blue
line) and of the attractor (black line) for all combinations of
$(\lambda,\mu)$. In addition, each sub-figure also depicts $40$ simulations
that are initiated at the boundary of $\mathcal{E}(W^{-1},1)$, which all
converge the attractor $\mathcal{E}(W^{-1},\varepsilon)$. For the particular
case $(\lambda,\mu)=(0.05,0.4)$, we have obtained
$W=\begin{bmatrix}78.67&-14.16\\\ -14.16&27.09\end{bmatrix}$ and
$\varepsilon=79.54$.

Figure 1: Estimations of the basin of attraction and of the attractor obtained
by solving the model-based optimization problem (11) for various $\lambda$ and
$\mu$.

One can observe on Figure 1 that increasing the magnitude of the noise (i.e.,
$\lambda$), increases the size of the attractor, while the size of
approximation of the basin of attraction remains almost the same.
Interestingly, Figure 1 shows that the tuning parameter $\mu$ has an important
effect on the solution. Indeed, for small values of $\mu$, the approximation
of the basin of attraction is larger, but the approximation of the attractor
is rather poor, since the black ellipsoids are very large compared to the
chattering around the origin. Reversely, when $\mu$ is large, the optimization
process provides a more accurate approximation of the attractor but at the
price of drastically reducing the size of the approximation of the basin of
attraction. This shows that the tuning parameter $\mu$ plays the role of the
balance between the optimization of estimations of the basin of attraction and
of the attractor, which cannot be performed simultaneously. This compromise
can be expected from the LMI condition, since having $(1-\mu)$ small, implies
that $W$ has to be large (block $(1,1)$ of $\Phi(\mathcal{G})$) and having
$\mu$ small implies that $W\varepsilon^{-1}$ has to be large (block $(3,3)$ of
$\Phi(\mathcal{G})$).

(a) Simulations obtained with a collection of $p=5$ data.

(b) Simulations obtained with a collection of $p=20$ data.

Figure 2: Estimations of the basin of attraction and of the attractor obtained
by solving the data-based optimization problem (11) for various $\lambda$ and
$\mu$.

###  V-B Data-driven local stabilization

This section illustrates the impact of the data-driven approach compared to
the model-driven one. Let us first present the construction of the data.
Following the dynamics of the system, we have selected $p=5$ and $p=20$ random
values of $x_{i}$ and $u_{i}$ in (19) where each component of $x_{i}$ lies in
$[-1,1]$. Then, we have created the random noise signal $w_{i}$ so that the
collection of noise vector $\omega$ belongs to $\Omega$ as described in
Assumption 2 with $\Delta_{\omega}=0.05I$. Altogether, vector
$\mathcal{X}^{+}$ has been computed using equation (21).

The results are illustrated in Figure 2. More precisely, Figures 2(a) and 2(b)
illustrate the results obtained by solving the data-driven stabilization
conditions of Theorem 2 with $p=5$ and $p=20$ experiments, respectively. As
for the model-based condition of Theorem 1, one can see that increasing the
magnitude of the noise leads to a reduction of the size of the estimation of
the basin of attraction (dashed blue ellipsoids), and also an increase of the
size of the attractor (black ellipsoids). The same effect of the tuning
parameters can be also detected.

Compared to the model-driven results in Figure 1, it occurs a notable
reduction of both ellipsoids for each case. This is due to the fact that the
model-based solution has no uncertainties in the matrices $A$ and $B$, while
the data-driven condition includes the inherent uncertainties due to the noise
affecting the data.

Comparing now both data-driven result, one can see that each estimation of the
basin of attraction obtained for $p=5$ are smaller than the one derived with
$p=20$. This can be interpreted by the fact that increasing the number of data
experiments give more information on the system so that the uncertainties due
to the noise are reduced. In other words the set of allowable matrices in
$\mathcal{C}$ becomes “smaller” as the number of data increases.

As a last comment, it seems surprising that the second column of Figure2(a)
(corresponding to $\lambda=0.05$) shows smallest estimation of the basin of
attraction than the third column (corresponding to $\lambda=0.1$). In light of
the previous paragraph, a possible explanation is again related to the set of
allowable matrices $\mathcal{C}$. As the data have been generated randomly, a
possible explanation is that set $\mathcal{C}$ obtained with a random noise
signal in $\Omega$ with $\lambda=0.05$ is “smaller” than the one obtained with
$\lambda=0.1$.

##  VI Conclusion

This paper addressed the problem of providing a data-driven solution to the
local stabilization of linear systems subject to input saturation. Using
model-based solution to this well-studied problem, a systematic method to
transform model-driven into data-driven LMI conditions is presented thanks to
some adequate rewriting of the conditions. Although this technical solution is
shown to be equivalent to some recent advanced results, its main advantage
relies on its simplicity and its potential to be applicable to a broad class
of problems of stabilization of (non)linear discrete-time systems.

The proposed results pave the way for future work, as in particular the
possibility to consider more complex dynamics or isolated nonlinearities

## References

  * [1] J. Berberich, A. Koch, C.W. Scherer, and F. Allgöwer.  Robust data-driven state-feedback design.  In Proceedings of the IEEE American Control Conference, pages 1532–1538. IEEE, 2020.
  * [2] J. Berberich, J. Köhler, M. A Müller, and F. Allgöwer.  Data-driven model predictive control with stability and robustness guarantees.  IEEE Trans. on Automatic Control, 66(4):1702–1717, 2020.
  * [3] A. Bisoffi, C. De Persis, and P. Tesi.  Data-driven control via Petersen’s lemma.  Automatica, 145:110537, 2022.
  * [4] S. Boyd, L. El Ghaoui, E. Feron, and V. Balakrishnan.  Linear Matrix Inequalities in System and Control Theory.  Society for Industrial and Applied Mathematics, 1997.
  * [5] V. Breschi, C. De Persis, S. Formentin, and P. Tesi.  Direct data-driven model-reference control with Lyapunov stability guarantees.  In Proceedings of the IEEE Conference on Decision and Control, pages 1456–1461. IEEE, 2021.
  * [6] C. De Persis and P. Tesi.  Formulas for data-driven control: Stabilization, optimality, and robustness.  IEEE Trans. on Automatic Control, 65(3):909–924, 2019.
  * [7] Y. Ebihara, D. Peaucelle, and D. Arzelier.  S-variable approach to LMI-based robust control, volume 6.  Springer, 2015.
  * [8] J.M. Gomes Da Silva and S. Tarbouriech.  Antiwindup design with guaranteed regions of stability: an LMI-based approach.  IEEE Trans. on Automatic Control, 50(1):106–111, 2005.
  * [9] Z.S. Hou and Z. Wang.  From model-based control to data-driven control: Survey, classification and perspective.  Information Sciences, 235:3–35, 2013.
  * [10] T. Hu and Z. Lin.  Control systems with actuator saturation: analysis and design.  Birkhauser, Boston, 2001.
  * [11] D. Piga, S. Formentin, and A. Bemporad.  Direct data-driven control of constrained systems.  IEEE Transactions on Control Systems Technology, 26:1422–1429, 2018\.
  * [12] I. Postlethwaite, M.C. Turner, and G. Herrmann.  Robust control applications.  Annual Reviews in Control, 31(1):27–39, 2007.
  * [13] C.W. Scherer.  Theory of robust control.  Citeseer, 2001.
  * [14] A. Seuret, C. Albea, and F. Gordillo.  LMI relaxations and its application to data-driven control design for switched affine systems.  Submitted to Int. Journal of Robust and Nonlinear Control, 2022\.
  * [15] A. Seuret and S. Tarbouriech.  A data-driven approach to the $\mathcal{l}_{2}$ stabilization of linear systems subject to input saturations.  IEEE Control Systems Letters, under review, 2023.
  * [16] S. Tarbouriech, G. Garcia, J.M. Gomes da Silva Jr., and I. Queinnec.  Stability and Stabilization of Linear Systems with Saturating Actuators.  Springer, 2011.
  * [17] S. Tarbouriech and J.M. Gomes Da Silva Jr.  Admissible polyhedra for discrete-time linear systems with saturating controls.  In Proceedings of the 1997 American Control Conference (Cat. No. 97CH36041), volume 6, pages 3915–3919. IEEE, 1997.
  * [18] H.J. Van Waarde, M.K. Camlibel, and M. Mesbahi.  From noisy data to feedback controllers: non-conservative design via a matrix S-lemma.  IEEE Trans. on Automatic Control, 67(1):162 – 175, 2022.
  * [19] H.J. Van Waarde, J. Eising, H.L. Trentelman, and M.K. Camlibel.  Data informativity: a new perspective on data-driven analysis and control.  IEEE Trans. on Automatic Control, 65(11):4753–4768, 2020.
  * [20] M. Vassilaki, J.C. Hennet, and G. Bitsoris.  Feedback control of linear discrete-time systems under state and control constraints.  International Journal of control, 47(6):1727–1735, 1988.

## Appendix-Proof of Lemma 2

###### Proof:

The proof is divided into two steps.
(i)$\Rightarrow$(ii): The first step of the proof is to find an appropriate
expression of matrices $\mathcal{A}$ that belong to $\Sigma_{\mathcal{N}}$.
For any matrix $\mathcal{A}$ in $\Sigma_{\mathcal{N}}$, it holds

$\begin{array}[]{lcl}0&\\!\\!\preceq&\\!\\!\mathcal{R}-\begin{bmatrix}I_{n_{3}}\\\
\mathcal{A}\\\
\end{bmatrix}^{\\!\top}\begin{bmatrix}\mathcal{R}+\mathcal{N}_{1}&\\!\\!\mathcal{N}_{2}\\\
\ast&\\!\\!\mathcal{N}_{3}\\\ \end{bmatrix}\left[\begin{matrix}I_{n_{3}}\\\
\mathcal{A}\\\ \end{matrix}\right]\\\ \end{array}$

where $\mathcal{R}$ is any matrix in $\mathbb{S}^{n_{3}}$. In addition,
$\mathcal{N}_{3}\succ 0$ ensures that there exists a matrix $\mathcal{R}$ such
that $\begin{bmatrix}\mathcal{R}+\mathcal{N}_{1}&\\!\\!\mathcal{N}_{2}\\\
\ast&\\!\\!\mathcal{N}_{3}\end{bmatrix}$ is positive definite. That allows
applying the Schur complement as follows to obtain

$\begin{bmatrix}\begin{bmatrix}\mathcal{R}+\mathcal{N}_{1}&\\!\\!\mathcal{N}_{2}\\\
\ast&\\!\\!\mathcal{N}_{3}\\\ \end{bmatrix}^{-1}&\
\left[\begin{matrix}I_{n_{3}}\\\ \mathcal{A}\\\ \end{matrix}\right]\\\
\ast&\mathcal{R}\end{bmatrix}\succ 0$

Note that it is not the usual way to apply the Schur Complement, but this dual
way has been considered to keep block $(1,2)$ (resp. $(2,1)$) with
$\left[\begin{smallmatrix}I_{n_{3}}\\\ \mathcal{A}\\\
\end{smallmatrix}\right]$ (resp. its transpose). Next, pre- and post-multiply
the previous inequality by

$\begin{bmatrix}\boldsymbol{0}_{n_{1},{n_{3}}}&\mathcal{Z}&\boldsymbol{0}_{{n_{1}},{n_{3}}}\\\
\boldsymbol{0}_{n_{3}}&\boldsymbol{0}_{{n_{3}},{n_{2}}}&I_{n_{3}}\end{bmatrix}$

and its transpose, respectively, where $\mathcal{Z}$ is any matrix in
$\mathbb{R}^{n_{1}\times{n_{2}}}$, i.e. of the same dimensions as
$\mathcal{M}_{2}$. This yields that having
$\mathcal{A}\in\Sigma_{\mathcal{N}}$ implies

$\begin{array}[]{l}\begin{bmatrix}\begin{bmatrix}\boldsymbol{0}_{{n_{1}},{n_{3}}}\\\
\mathcal{Z}^{\\!\top}\end{bmatrix}^{\\!\\!\top}\\!\\!\begin{bmatrix}\mathcal{R}+\mathcal{N}_{1}&\\!\\!\mathcal{N}_{2}\\\
\ast&\\!\\!\mathcal{N}_{3}\\\
\end{bmatrix}^{-1}\\!\\!\begin{bmatrix}\boldsymbol{0}_{n_{1},{n_{3}}}\\\
\mathcal{Z}^{\\!\\!\top}\end{bmatrix}^{\\!\\!\top}&\\!\\!\\!\\!\begin{bmatrix}\boldsymbol{0}_{{n_{1}},{n_{3}}}\\\
\mathcal{Z}^{\\!\top}\end{bmatrix}^{\\!\\!\top}\\!\\!\left[\begin{matrix}I_{n_{3}}\\\
\mathcal{A}\\\ \end{matrix}\right]\\\
\ast&\\!\\!\\!\\!\mathcal{R}\end{bmatrix}\\\
=\begin{bmatrix}\begin{bmatrix}\boldsymbol{0}_{n_{1},{n_{3}}}\\\
\mathcal{Z}^{\\!\top}\end{bmatrix}^{\\!\\!\top}\begin{bmatrix}\mathcal{R}+\mathcal{N}_{1}&\\!\\!\mathcal{N}_{2}\\\
\ast&\\!\\!\mathcal{N}_{3}\\\
\end{bmatrix}^{-1}\begin{bmatrix}\boldsymbol{0}_{{n_{1}},{n_{3}}}\\\
\mathcal{Z}^{\\!\top}\end{bmatrix}^{\\!\top}&\mathcal{Z}\mathcal{A}\\\
\ast&\mathcal{R}\end{bmatrix}\succ 0.\end{array}$ (32)

The previous calculations ensure that inequality (25) can be rewritten as
$\mathcal{M}(\mathcal{A})\succ 0$, for all $\mathcal{A}$ such that (32) holds
true. Using an S-procedure, this statement is equivalent to the existence of a
positive scalar $\eta>0$ such that

$\begin{array}[]{l}\begin{bmatrix}\mathcal{M}_{1}&\mathcal{M}_{2}\mathcal{A}\\\
\ast&\mathcal{M}_{3}\end{bmatrix}\\\
-\eta\begin{bmatrix}\begin{bmatrix}\boldsymbol{0}_{{n_{1}},{n_{3}}}\\\
\mathcal{Z}^{\\!\top}\end{bmatrix}^{\\!\top}\begin{bmatrix}\mathcal{R}+\mathcal{N}_{1}&\\!\\!\mathcal{N}_{2}\\\
\ast&\\!\\!\mathcal{N}_{3}\\\
\end{bmatrix}^{-1}\begin{bmatrix}\boldsymbol{0}_{n_{1},{n_{3}}}\\\
\mathcal{Z}^{\\!\top}\end{bmatrix}^{\\!\top}&\mathcal{Z}\mathcal{A}\\\
\ast&\mathcal{R}\end{bmatrix}\succ 0,\end{array}$

which together with the Schur complement can be written

$\begin{bmatrix}\mathcal{M}_{1}&(\mathcal{M}_{2}-\eta\mathcal{Z})\mathcal{A}&\boldsymbol{0}_{n_{1},{n_{3}}}&\eta\mathcal{Z}\\\
\ast&\mathcal{M}_{3}-\eta\mathcal{R}&\boldsymbol{0}_{n_{3}}&\boldsymbol{0}_{{n_{3}},{n_{2}}}\\\
\ast&\ast&\mathcal{R}+\mathcal{N}_{1}&\\!\\!\mathcal{N}_{2}\\\
\ast&\ast&\ast&\\!\\!\mathcal{N}_{3}\\\ \end{bmatrix}\succ 0$

Selecting $\eta\mathcal{R}=\mathcal{M}_{3}-\epsilon I_{n_{3}}$, for any
arbitrarily small $\epsilon>0$ and $\eta\mathcal{Z}=\mathcal{M}_{2}$ so that
the second (block) row and columns are zero and the uncertain matrix
$\mathcal{A}$ disappears. The previous matrix becomes

$\begin{bmatrix}\mathcal{M}_{1}&\boldsymbol{0}_{n,p}&\boldsymbol{0}_{n_{1},{n_{3}}}&\mathcal{M}_{2}\\\
\ast&\epsilon
I_{n_{3}}&\boldsymbol{0}_{n_{3}}&\boldsymbol{0}_{{n_{3}},{n_{2}}}\\\
\ast&\ast&\mathcal{M}_{3}+\eta\mathcal{N}_{1}-\epsilon
I_{n_{3}}&\eta\mathcal{N}_{2}\\\
\ast&\ast&\ast&\eta\mathcal{N}_{3}\end{bmatrix}\succ 0$

which is equivalent to (27) knowing that $\epsilon$ is an arbitrarily small
scalar.
(ii)$\Rightarrow$(i): Pre- and post-multiply (27) by
$\begin{bmatrix}I_{n_{1}}&\boldsymbol{0}_{{n_{1}},{n_{3}}}&\boldsymbol{0}_{{n_{1}},{n_{3}}}\\\
\boldsymbol{0}_{{n_{3}},{n_{1}}}&I_{n_{3}}&\mathcal{A}^{\\!\top}\end{bmatrix}$
and its transpose, respectively, leads to

$\begin{array}[]{lcl}0&\\!\\!\prec&\\!\\!\begin{bmatrix}\mathcal{M}_{1}&\mathcal{M}_{2}\mathcal{A}\\\
\ast&\mathcal{M}_{3}\end{bmatrix}+\eta\begin{bmatrix}\boldsymbol{0}_{n_{1}}&\boldsymbol{0}_{n_{1},{n_{3}}}\\\
\ast&\left[\begin{matrix}I_{n_{3}}\\\ \mathcal{A}\\\
\end{matrix}\right]^{\\!\\!\top}\\!\\!\begin{bmatrix}\mathcal{N}_{1}&\\!\\!\mathcal{N}_{2}\\\
\ast&\\!\\!\mathcal{N}_{3}\\\ \end{bmatrix}\left[\begin{matrix}I_{n_{3}}\\\
\mathcal{A}\\\ \end{matrix}\right]\end{bmatrix}\\\
&\\!\\!\prec&\\!\\!\begin{bmatrix}\mathcal{M}_{1}&\mathcal{M}_{2}\mathcal{A}\\\
\ast&\mathcal{M}_{3}\end{bmatrix}\end{array}$

where the last inequality holds only for matrices $\mathcal{A}$ that belongs
to $\Sigma_{\mathcal{N}}$, which concludes the proof. ∎


#
Beyond QUBO and HOBO formulations, solving the Travelling Salesman Problem on
a quantum boson sampler

Daniel Goldsmith<EMAIL_ADDRESS>Digital Catapult, 101
Euston Road, London NW1 2RA  Joe Day-Evans  Digital Catapult, 101 Euston Road,
London NW1 2RA

###### Abstract

The Travelling Salesman Problem (TSP) is an important combinatorial
optimisation problem, and is usually solved on a quantum computer using a
Quadratic Unconstrained Binary Optimisation (QUBO) formulation or a Higher
Order Binary Optimisation (HOBO) formulation. In these formulations, penalty
terms are added to the objective function for outputs that don’t map to valid
routes.

We present a novel formulation which needs fewer binary variables, and where,
by design, there are no penalty terms because all outputs from the quantum
device are mapped to valid routes. Simulations of a quantum boson sampler were
carried out which demonstrate that larger networks can be solved with this
penalty-free formulation than with formulations with penalties. Simulations
were successfully translated to hardware by running a non-QUBO formulation
with penalties on an early experimental prototype ORCA PT-1 boson sampler.
Although we worked with a boson sampler, we believe that this novel
formulation is relevant to other quantum devices.

This work shows that a good embedding for combinatorial optimisation problems
can solve larger problems with the same quantum computing resource. The
flexibility of boson sampling quantum devices is a powerful asset in solving
combinatorial optimisation problem, because it enables formulations where the
output string is always mapped to a valid solution, avoiding the need for
penalties.

##  1 Introduction

We are in the ‘Noisy Intermediate Scale era of quantum computing’ (NISQ), a
phrase coined by Preskill [1], and contemporary quantum devices have high
error rates, limited connectivity and lack scale. Combinatorial optimization,
where an objective function with binary arguments is minimised, has been
identified as a potential early use case for quantum computing, along with
quantum chemistry and quantum machine learning [2].

An important combinatorial optimisation problem is the well known NP-Hard
Travelling Salesman Problem (TSP), which aims to find the shortest Hamiltonian
cycle through a graph $G=(V,E)$ of $N$ locations represented as vertices $V$
in the graph, and edges $E$ between the vertices representing the connections.
Each edge between vertices $u$ and $v$ has a distance $D_{u,v}$. The TSP
problem consists of finding a permutation $\pi$ of the $N$ locations, that
minimises the distance $C(\pi)$ of visiting each of the $N$ locations once,
and returning to the start. The distance $C(\pi)$ is shown in Equation 1 where
we note that $\pi(N+1)=\pi(1)$ because the salesman must return to the
starting location.

$C(\pi)=\sum_{i=1}^{N}D_{\pi(i),\pi(i+1)}$ (1)

As well as classical methods, quantum methods to solve TSP include quantum
inspired algorithms [3], adiabatic quantum computing [4, 5], quantum heuristic
based on a Grover search [6], quantum phase estimation [7], quantum
approximate optimization algorithm (QAOA)[8, 9, 10, 11] and using qudits
rather than qubits [12].

The standard approach to solving this problem with a quantum computer is the
Quadratic Unconstrained Binary Optimisation (QUBO) formulation as outlined by
Lucas [13], using either quantum annealing devices or the QAOA algorithm for
circuit model devices. In 2024 Ramezani and others [8] designed and
implemented a Higher Order Binary Optimisation (HOBO) algorithm for QAOA where
the number of qubits scales as $\mathcal{O}(N\log_{2}{}N)$, rather than
$\mathcal{O}(N^{2})$ for the standard approach, but at the cost of increasing
the number of 2-qubit gates and hence the circuit depth. In these QUBO and
HOBO formulations, outputs that don’t map to a valid solution attract penalty
terms, and these invalid solutions may lead to a flat objective function that
is difficult to optimise. The QUBO and HOBO formulations are described in the
Literature Review.

In the Methods section we present a novel formulation where, by design, there
are no penalty terms because all outputs are mapped to valid routes, and
explain how this formulation works with a boson sampler. This penalty-free
formulation scales as $\mathcal{O}(N\log_{2}{}N)$ and requires slightly fewer
binary variables than Ramezani’s formulation [8], and we believe that it may
be relevant for other quantum devices where the output is a binary string with
a probability distribution depending on parameters that can be modified by a
classical optimiser.

In the Experiments, Data and Results section we describe simulations of a
quantum boson sampler that demonstrate that larger networks can be solved with
this penalty-free formulation than with formulations with penalties, such as
the QUBO formulation, or a formulation where the location reference is encoded
in binary, similar to [8]. We also describe how we translated simulations to
hardware by running a non-QUBO formulation with penalty terms on an early
experimental prototype ORCA PT-1 boson sampler.

##  2 Literature review

###  2.1 QUBO formulation of TSP

In the Quadratic Unconstrained Binary Optimisation (QUBO) formulation of a
combinatorial optimisation problem the problem Hamiltonian $H_{P}$ is written
as a function of binary variables $x_{i}$ of no higher than quadratic order,
as shown in Equation 2, and is minimised. The problem Hamiltonian contains the
objective function to be minimised as well as penalty terms for any
constraints, multiplied by arbitary constants to ensure that it is not
energetically favourable to break the constraints:

$H_{P}=\sum_{i,j}Q_{i,j}x_{i}x_{j}$ (2)

If the binary bit string is represented as a vector $\boldsymbol{x}$ this is
equivalent to the matrix formulation shown in Equation 3:

$H_{P}=\boldsymbol{x}^{T}\boldsymbol{Q}\boldsymbol{x}$ (3)

Lucas [13] describes a formulation of TSP as a QUBO. Labelling the vertices
$0,1,...,N-1$, where $N$ is the number of vertices, or locations to be
visited, $N^{2}$ binary variables $x_{v,i}$ are contemplated, where $v$
represents the vertex and $i$ its order in a prospective cycle. Since, without
loss of generality, the cycle can always start at vertex $0$, the number of
binary variables required is reduced to $\left(N-1\right)^{2}$.

As each edge between vertices $u$ and $v$ has a distance $D_{u,v}$ then the
total distance is found as the sum of the distances for each edge:

$H_{obj}=\sum_{(uv)\in E}D_{u,v}\sum_{j=0}^{N-1}x_{u,j}x_{v,j+1}$ (4)

Noting that the $v^{th}$ vertex can only appear once in the cycle, and that
there must be a $j^{th}$ node in the cycle for each $j$, these constraints are
encoded in the penalty Hamiltonian:

$H_{pen}=A\sum_{v=0}^{N-1}\left(1-\sum_{j=0}^{N-1}x_{v,j}\right)^{2}+A\sum_{j=0}^{N-1}\left(1-\sum_{v=0}^{N-1}x_{v,j}\right)^{2}$
(5)

with $A$ large enough so that it is never favourable to violate the
constraints of $H_{pen}$.

In summary, in the QUBO formulation of TSP, we minimise the problem
Hamiltonian:

$H_{P}=H_{obj}+H_{pen}$ (6)

####  2.1.1 Solving the QUBO formulation of TSP on a quantum annealing device

The QUBO formulation in 6 can be solved with a quantum annealing device using
adiabatic quantum computing because the quadratic problem Hamiltonian $H_{P}$
in Equation 6 can be mapped to the Hamiltonian for a quantum annealing device
as shown in Equation 7 which is also quadratic:

$H_{P}=\sum_{i}h_{i}\sigma_{z}^{(i)}+\sum_{i,j}J_{i,j}\sigma_{z}^{(i)}\sigma_{z}^{(j)}$
(7)

where $\sigma_{z}^{(i)}$ denotes the Pauli-Z matrix acting on the $i^{th}$
qubit, the $h_{i}$ biases are individually set for each qubit, and the
$J_{i,j}$ biases control the strength of interactions between the $i^{th}$ and
$j^{th}$ qubit. Since each binary variable in the problem is mapped to one
qubit, and more qubits may be needed if one qubit can’t be connected to
another, the problem requires at least $\mathcal{O}(N^{2})$ qubits

Only small TSP networks can be solved at present using Quantum Annealing. For
example, in 2021 Jain [4] struggled to solve networks of more than 10
locations with the D-Wave quantum annealer.

####  2.1.2 Solving the QUBO formulation of TSP on a circuit model device

In the QAOA algorithm [9] the time evolution of the adiabatic quantum
computing Hamiltonian is Trotterised using a paramaterized quantum circuit
which encodes a Unitary operator $U(\bm{\gamma},\bm{\beta})$, where
$\bm{\gamma}$ and $\bm{\beta}$ are sets of variational parameters.

The Unitary operator takes the form shown in Equation 8

$U(\bm{\gamma},\bm{\beta})=\prod_{l=1}^{L}\exp\left(-i\beta_{l}H_{M}\right)\exp\left(-i\gamma_{l}H_{P}\right)$
(8)

where $L$ denotes the number of optimization layers. The problem Hamiltonian
$H_{P}$ is as described in Equation 7 and the $H_{M}$ is a mixing Hamiltonian.

###  2.2 HOBO formulation of TSP

The QUBO formulation above uses a One-Hot encoding scheme. Only a few
combinations of binary variables are valid solutions because if the binary
variables $x_{v,i}$ are represented as a matrix, few matrices are valid
permutation matrices. It is reasonable to ask if we can do better with a
Higher Order Binary Optimisation (HOBO) formulation.

In 2021 Vargas-Calderon and others [12] formulated TSP using qudits rather
than qubits. They encoded a valid route, for example, $0\to 3\to 2\to 1\to 0$,
using an ordered set of four 4-level qudits as a quantum state
$\ket{0}\otimes\ket{3}\otimes\ket{2}\otimes\ket{1}$. They devised a
Hamiltonian that returns the total distance for valid cycles, and a penalty
term otherwise, and is a sum of quadratic terms involving no more than two
qudits. The team simulated valid solutions for up to 40 locations, although
they could not yet run this formulation on an existing device, since
contemporary devices use qubits, not qudits.

In 2024 Ramezani and others [8] built on Vargas-Calderon’s work by introduced
a HOBO formulation. Routes are encoded using the binary representation of the
location number. For example, a valid cycle $0\to 3\to 2\to 1\to 0$ is encoded
as a binary string of ${00111001}_{2}$ = ${00}_{2}$ ${11}_{2}$ ${10}_{2}$
${01}_{2}$.

Since each binary variable is mapped to a qubit, this reduces the qubits
required to $N\left\lceil\log_{2}{}N\right\rceil$, since there are $N$
locations, and each location reference requires
$\left\lceil\log_{2}{}N\right\rceil$ binary variables. In the problem
Hamiltonian the reduction in the number of binary variables comes at the
expense of introducing higher order binary terms. Because of this the number
of two-qubit gates scales as $\mathcal{O}(N^{4}\log_{2}{}N)$ rather than
$\mathcal{O}(N^{3})$ for the conventional QAOA algorithm. Simulation by
Ramezani’s team showed that the HOBO formulation returned better solutions
than standard methods, but only for small networks with four locations.

##  3 Methods

###  3.1 Formulation of the TSP without penalty terms

A novel formulation is presented, which unlike the QUBO and HOBO formulations
above, does not use penalty terms. Although this formulation was implemented
on a quantum boson sampler, we believe the penalty free formulation may be
relevant for other quantum devices where the output is a binary string with a
probability distribution depending on parameters that can be modified by a
classical optimiser.

An initial sequential list of location numbers is created and items are moved
from the initial list to a re-ordered list until the initial list is empty and
the re-ordered list holds a valid cycle. To achieve this the output bit string
from the quantum device is split into smaller strings, which are interpreted
as the binary encoding of indices pointing to the next item of the ordered
list to be moved. Use of modulo arithmetic prevents errors where the index is
too high to reference a valid element of the ordered list.

The required bit string length $L_{1}$ is found by summing the length of these
binary index strings, noting that the length of the string is
$\left\lceil\log_{2}(i)\right\rceil$ where $i$ is the number of elements
remaining in the ordered list. $L_{1}$ is given by Equation 9 where $N$ is the
number of locations:

$L_{1}=\sum_{i=1}^{N-1}\left\lceil\log_{2}(i)\right\rceil$ (9)

The number of binary variables $L_{1}$ required is less than the
$N\left\lceil\log_{2}{}N\right\rceil$ for Ramezani’s formulation [8], because
this encoding is more efficient than Ramezani’s binary conversion of each
location code, and substantially less than the $\left(N-1\right)^{2}$ needed
for the standard QUBO formulation [13].

Pseudo code for implementation of this algorithm on a boson sampler is shown
as Algorithm 1 in the Appendix. This formulation was simulated on a boson
sampling device, as described below:

###  3.2 Optimisation on a boson sampling device

In a photonic boson sampler indistinguishable photons are injected into a
linear interferometer and the output distribution is measured. The output
photon probability distribution can be shown to be hard to simulate
classically [14], a necessary condition for quantum advantage. Bradler and
Wallner [15] explain how a photonic boson sampling device can carry out
optimisation using Variational Quantum Algorithms (VQAs) [2], and illustrate
this by solving a Mobius graph problem in a QUBO formulation, and a portfolio
optimization problem with binary investment problem in a non-QUBO formulation.

An example of a photonic boson sampler is the ORCA Computing PT-Series quantum
computer, which produces of a train of optical pulses that interfere with each
other via programmable beam splitters located within one or more fibre loops,
with a parameter vector given by the beam splitter angles
$\boldsymbol{\theta}=\theta^{(i)}$, as in Figure 1.

Figure 1: The ORCA PT-Series following [15] with a train of optical pulses
that interfere with each other via programmable beam splitters parameterised
by the beam splitter angles $\theta^{(i)}$

The output of this boson sampler, a quantum Fock state, with an integer number
of photons in each output position, is measured and then mapped to a bit
string using component wise parity functions. An example of a component-wise
parity function mapping from the optical circuit outputs to bit strings is
show in Figure 2 with even or odd photon counts mapped to $0$ or $1$ [15]. It
can shown that the device needs to be run in four configurations,
corresponding to sampling steps for odd and even parity counts and for $M$ and
$M-1$ photons, so that all bit strings can be sampled.

Figure 2: Component wise parity function from optical circuit outputs to bit
strings based on [15].

The objective function is evaluated classically for each output bit string and
then optimised classically by tuning the parameter vector of beam splitter
angles $\boldsymbol{\theta}=\theta^{(i)}$ in a feedback loop, which is
repeated until the objective function converges, at which point the lowest
energy found and the corresponding bit string found are recorded.

###  3.3 Comparison with other formulations

The performance of the penalty-term free formulation on the bosom sampler was
assessed by comparing it to a non-QUBO and a QUBO formulation, both with
penalty terms

####  3.3.1 Non-QUBO formulation of the TSP with penalty terms on a boson
sampler

Using the flexibility of the boson sampler, Ramezani’s [8] technique can be
used to encode the boson sampler output binary string as a location list,
using a binary representation of the location label, but without explicitly
defining the full HOBO formulation. If the location list contains duplicate
location labels, or location labels higher than $(N-1)$ then a constant
penalty term $H_{pen}$, as defined in Equation 10, is applied with a value of
$\rho$ high enough so that these terms are not selected in preference to a
valid cycle. In the simulations a value of $\rho=5$ was used. Otherwise the
location list defines a valid cycle, and the cycle distance is calculated.

$H_{pen}=\rho\sum_{i=0}^{N-1}D_{0,i}$ (10)

Pseudo code is shown as Algorithm 2 in the Appendix. This formulation was
simulated on a boson sampling device and also run on the ORCA PT-1 device.

Whilst this formulation was able to find solutions for up to 15 locations, as
discussed in the results below, the solutions were not high quality. The
reason seems likely to be that many of the bit strings result in the penalty
function being applied, so the cost landscape is flat, and the optimiser
struggles to find the minimum.

####  3.3.2 QUBO formulation of the TSP with penalty terms on a boson sampler

A boson sampler output binary string of length $L_{3}=(N-1)^{2}$ was reshaped
as a matrix with elements $x_{v,j}$, to reproduce the QUBO formulation.
Distances were calculated for binary strings representing valid routes, and
penalty terms were applied if the route was not valid. Pseudo code is shown as
Algorithm 3 in the Appendix. This formulation was simulated on a boson
sampling device. More binary variables were required, and only small networks
could be solved, as discussed in the results below.

##  4 Experiments, data and results

###  4.1 Experiments carried out

The three boson sampler TSP formulations described above were simulated on a
boson sampler using the ORCA SDK with networks of between four and 48
locations. The solution quality ($Q_{sol}$) was determined as the shortest
distance found experimentally around a valid cycle ($D_{min\\_exp}$) as a
percentage of the best known minimum distance ($D_{min\\_best}$) as shown in
Equation 11

$Q_{sol}=100\frac{D_{min\\_exp}}{D_{min\\_best}}$ (11)

A solution quality of $Q_{sol}=100\%$ represents a solution where the shortest
distance found is optimal, and a lower percentage represents a sub-optimal
solution.

Key training graphs are presented showing how the objective function varies
with the number of training iterations. Four different plots are shown on each
graph represent the four sampling configurations for odd and even parity
counts and for $M$ and $M-1$ photons, as explained in the description of the
boson sampling device.

The simulations were run on a laptop with an Intel(R) Core(TM) i5-10210U CPU
clocked at 1.60GHz with 8.00 GB (7.76 GB usable) RAM running Windows 11. We
used the latest ORCA SDK which can be obtained from ORCA Computing. A link to
an earlier version of the SDK can be found in [15]. We have placed the
relevant Jupyter Notebooks and Python modules on a Digital Catapult GitHub
repository [16].

To assess the impact of hyper-parameters, simulations were run with different
numbers of iterations and learning rates. The optimisation method was either
the parameter shift rule, or Simultaneous Perturbation Stochastic
Approximation (SPSA) [17], where instead of estimating the gradient for each
parameter by taking two function estimates per parameter, the gradient was
estimated with two function estimation steps based on a stochastic variation.

The non-QUBO based formulation with penalty terms was also run on an
experimental ORCA PT-1 boson sampler. Since this early device could only
accommodate six binary variables it was run for four locations.

###  4.2 Data

Data sets with $5$, $15$, $26$, $42$ and $48$ locations were obtained from
[18, 19], and a data set with $4$ locations was used for Unit Testing, where
the results of the simulations were checked manually.

###  4.3 Results

The non-QUBO formulation with penalties was run on a network of 4 locations on
an experimental ORCA-PT1 boson sampler, and the correct solution was found.

Results of simulation with different formulations, numbers of locations,
iterations, learning rates and optimisation methods are summarised in the
table below:

---
Table 1: Summary of simulation results

####  4.3.1 Simulation of small networks

Each of the three formulations simulated found the correct solution for
networks of $4$ or $5$ locations.

####  4.3.2 Simulations of a network with 15 locations

The QUBO formulation could not find a valid solution for 15 locations. We
believe this is because almost all binary strings were outside the valid
solution space, and attracted the same penalty term shown in Equation 10,
leading to a flat cost landscape, where it is difficult for the optimiser to
find a minimum.

Both non-QUBO formulations found valid solution for 15 locations. However, the
training graph with SPSA and a learning rate of 0.01 for the non-QUBO
formulation with penalties in Figure 3 shows that most bit strings have a cost
equal to the penalty term, apart from a few strings with a valid solution. The
average objective function does not reduce with the number of iterations, and
there is no convergence to an optimum solution.

Figure 3: Objective function against iterations for non-QUBO formulation with
penalties showing that because almost all bit strings attract the same penalty
term, no convergence is possible on a network with 15 locations. The four
plots represent the four sampling configurations.

The penalty free formulation solution quality was higher for 15 locations that
the formulation with penalties, with a solution quality of up to $90.7\%$.

####  4.3.3 Simulations of network with over 26 locations

The non-QUBO formulation with penalties could not find a valid solution for 26
locations. The penalty-free formulation found solutions with up to 48
locations, even though the solution quality declined with larger networks. The
penalty-free formulation showed evidence of training with the value of the
objective function decreasing during training. This is clearly seen with 26
locations, SPSA and a learning rate of 0.01, as shown in Figure 4, although it
was less pronounced in other network sizes.

Figure 4: Objective function against iterations for the penalty-free
formulation showing evidence of training for a network of 26 locations. The
four plots represent the four sampling configurations.

####  4.3.4 Impact of optimisation method

SPSA was generally much quicker, presumably because only two estimates were
made, and was noiser, usually with a lower solution quality.

##  5 Discussion

Simulations of a novel, penalty free formulation of TSP on a boson sampler has
shown useful results, with valid solutions for networks of up to 48 locations,
even though the solution quality declined with larger networks. The use of
penalties may be associated with flat, difficult to optimise, training graphs.

Whilst the problem size and solution quality is not yet comparable to best of
breed classical methods, the simulation carried out was of an early boson
sampler with only one fibre loop to store photons, and more advanced devices
are likely soon. Also, much information from the measurement of the Fock space
is discarded, since integer values are mapped to binary parities. Accessing
these integers may lead to an improved formulation.

Future work is needed to investigate the performance of this penalty free
formulation by simulating and/or running it on larger and more complex boson
samplers, or other devices.

##  6 Conclusions

This work shows that:

  1. 1.

It is important to use a good embedding for TSP, and by inference, other
combinatorial problems.

  2. 2.

An embedding designed to always map the binary output string of a quantum
device to a valid solution, avoiding the need for penalty terms, was able to
solve larger routes than formulations with penalties.

##  7 Acknowledgments

We are very grateful to Dr Clive Emary from Frazer-Nash and Dr William
Clements from ORCA for their invaluable contribution to this work. Many thanks
to the ORCA Computing team for allowing access to an ORCA PT-1 experimental
prototype to test the algorithm on a real device, and allowing access to the
latest version of their software development kit (SDK). Many thanks to Sheevan
Jayasinghe from Bahut, a QTAP participant, who contributed the TSP use case,
and an example coded in the QUBO formulation.

This work was carried out during the use case dissemination phase of the
Quantum Technology Access Programme (QTAP), part of the Quantum Data Centre of
the Future, project number 10004793, funded by the Innovate UK Industrial
Strategy Challenge Fund (ISCF).

## References

  * [1] John Preskill.  Fault-tolerant quantum computation.  https://arxiv.org/abs/quant-ph/9712048, 12 1997.
  * [2] Kishor Bharti, Alba Cervera-Lierta, Thi Ha Kyaw, Tobias Haug, Sumner Alperin-Lea, Abhinav Anand, Matthias Degroote, Hermanni Heimonen, Jakob S. Kottmann, Tim Menke, Wai-Keong Mok, Sukin Sim, Leong-Chuan Kwek, and Alán Aspuru-Guzik.  Noisy intermediate-scale quantum (nisq) algorithms.  Reviews of Modern Physics, 94, 1 2021.
  * [3] Vicente P. Soloviev, Concha Bielza, and Pedro Larrañaga.  Quantum-inspired estimation of distribution algorithm to solve the travelling salesman problem.  2021 IEEE Congress on Evolutionary Computation (CEC), pages 416–425, 2021.
  * [4] Siddhartha Jain.  Solving the traveling salesman problem on the d-wave quantum computer.  In Frontiers of Physics, 2021.
  * [5] Tien D. Kieu.  The travelling salesman problem and adiabatic quantum computation: an algorithm.  Quantum Information Processing, 18:1–19, 2018.
  * [6] Jeongho Bang, Junghee Ryu, Changhyoup Lee, Seokwon Yoo, James Lim, and Jinhyoung Lee.  A quantum heuristic algorithm for the traveling salesman problem.  Journal of the Korean Physical Society, 61:1944 – 1949, 2010.
  * [7] Karthik Srinivasan, Saipriya Satyajit, Bikash K. Behera, and Prasanta K. Panigrahi.  Efficient quantum algorithm for solving travelling salesman problem: An ibm quantum experience.  https://arxiv.org/abs/1805.10928, 2018\.
  * [8] Mehdi Ramezani, Sadegh Salami, Mehdi Shokhmkar, Morteza Moradi, and Alireza Bahrampour.  Reducing the number of qubits from $n^{2}$ to $n\log_{2}(n)$ to solve the traveling salesman problem with quantum computers: A proposal for demonstrating quantum supremacy in the nisq era.  https://arxiv.org/abs/2402.18530, 2024\.
  * [9] Edward Farhi, Jeffrey Goldstone, and Sam Gutmann.  A quantum approximate optimization algorithm.  http://arxiv.org/abs/1411.4028, 11 2014\.
  * [10] Yue Ruan, Samuel Marsh, Xiling Xue, Zhihao Liu, and Jingbo B. Wang.  The quantum approximate algorithm for solving traveling salesman problem.  Computers, Materials & Continua, 2020.
  * [11] Matthew Radzihovsky, Joey Murphy, and Mason Swofford.  A qaoa solution to the traveling salesman problem using pyquil.  https://cs269q.stanford.edu/projects2019/radzihovsky_murphy_swofford_Y.pdf, 2019\.
  * [12] Vladimir Vargas-Calderón, Nicolas Parra-A., Herbert Vinck-Posada, and Fabio A. González.  Many-qudit representation for the travelling salesman problem optimisation.  Journal of the Physical Society of Japan, 90(11):114002, November 2021.
  * [13] Andrew Lucas.  Ising formulations of many np problems.  https://arxiv.org/abs/1302.5843, 2 2013\.
  * [14] Scott Aaronson and Alexei Y. Arkhipov.  The computational complexity of linear optics.  https://arxiv.org/abs/1011.3245, 2010\.
  * [15] Kamil Bradler and Hugo Wallner.  Certain properties and applications of shallow bosonic circuits.  https://arxiv.org/abs/2112.09766, 12 2021\.
  * [16] Daniel Goldsmith.  Quantum technology access programme software repository.  https://github.com/digicatapult/QTAP, 2023\.
  * [17] James C. Spall.  Simultaneous perturbation stochastic approximation: Spsa - a method for system optimisation.  https://www.jhuapl.edu/SPSA/, 2001\.
  * [18] Tsp - data for the traveling salesperson problem.  https://people.sc.fsu.edu/ jburkardt/datasets/tsp/tsp.html, 2019\.  Reviewed 13 February 2024.
  * [19] Gerhard Reinelt.  TSPLIB–a traveling salesman problem library.  ORSA Journal on Computing, 3(4):376–384, 1991.

## Appendix

### Algorithms

In these algorithms we follow the convention from Python that the first item
of a list is indexed by $0$.

Algorithm 1 Formulation of TSP without penalty terms on a boson sampler

Input:

Number of locations $N$

Distance matrix $D_{u,v}$

Number of samples $N_{s}$

Output:

Minimum distance $D_{min}$

Cycle with minimum distance $LIST_{min}$

Calculate the bit string length $L_{1}$ from Equation 9

Configure the interferometer to produce bit strings $\boldsymbol{b}$ of length
$L_{1}$

for $i\leftarrow 1$ to $N_{s}$ do

Measure the output pattern generated by the interferometer.

Classical map the output pattern to a bit string $\boldsymbol{b}$ of length
$L_{1}$

Initialise an ordered list of locations $LIST_{ordered}=[0,1,...N-1]$

Initialise an empty list of locations $LIST_{cycle}=[\>]$

Move the $0^{th}$ element $0$ from $LIST_{ordered}$ to $LIST_{cycle}$ since,
without loss of generality, the cycle can start at location $0$.

for $j\leftarrow 1$ to $N-2$ do

Calculate the number of items $m_{j}$ left in $LIST_{ordered}$ as $m_{j}=N-j$

Calculate the length $l_{j}$ of bit string $\boldsymbol{b}_{j}$ as
$l_{j}=\left\lceil log_{2}(m_{j})\right\rceil$

Move $l_{j}$ bits from $\boldsymbol{b}$ to $\boldsymbol{b}_{j}$

Evaluate $\boldsymbol{b}_{j}$ as an integer $p_{j}$

Calculate an index $q_{j}=p_{j}mod(m_{j})$

Move the $q_{j}^{th}$ element from $LIST_{ordered}$ to $LIST_{cycle}$

end for

Move the final element $LIST_{ordered}$ to $LIST_{cycle}$

Calculate the distance $D_{i}$ for cycle $LIST_{cycle}$ using the distance
matrix $D_{u,v}$

if $D_{min}$ is initial then

$D_{min}\leftarrow D_{i}$

else if $D_{i}<D_{min}$ then

$D_{min}\leftarrow D_{i}$

$LIST_{min}\leftarrow LIST_{cycle}$

end if

end for

Iterate changing the parameter vector $\boldsymbol{\theta}$ until convergence
is reached.

Algorithm 2 Non-QUBO formulation of the TSP with penalty terms on a boson
sampler

Input:

Number of locations $N$

Distance matrix $D_{u,v}$

Number of samples $N_{s}$

Output:

Bit string for which a minimum distance is found $\boldsymbol{b}_{min}$

Minimum problem Hamiltonian $H_{min}$

Minimum distance $D_{min}$

Cycle with minimum distance $LIST_{min}$

Calculate the bit string length $L_{2}=\left\lceil
log_{2}(N)\right\rceil(N-1)$

Configure the interferometer to produce bit strings $\boldsymbol{b}$ of length
$L_{2}$

for $i\leftarrow 1$ to $N_{s}$ do

Measure the output pattern generated by the interferometer.

Initialise an empty list of locations $LIST_{cycle}=[\>]$

for $j\leftarrow 1$ to $N-2$ do

Move $\left\lceil log_{2}(N)\right\rceil$ bits from $\boldsymbol{b}$ to
$\boldsymbol{b}_{j}$

Evaluate $\boldsymbol{b}_{j}$ as an integer $p_{i}$

Append $p_{i}$ to $LIST_{cycle}$

end for

Augment $LIST_{cycle}$ by adding the missing location by elimination

if $LIST_{cycle}$ is valid then

Calculate the distance $D_{i}$ for cycle $LIST_{cycle}$ using the distance
matrix $D_{u,v}$

$H_{problem}\leftarrow D_{i}$

else

$H_{problem}\leftarrow H_{pen}$

end if

if $H_{min}$ is initial then

$H_{min}\leftarrow H_{problem}$

else if $H_{problem}<H_{min}$ then

$H_{min}\leftarrow H_{problem}$

$\boldsymbol{b}_{min}\leftarrow\boldsymbol{b}$

end if

end for

Iterate, changing the parameter vector $\boldsymbol{\theta}$ until convergence
is reached.

Calculate $LIST_{min}$ from $\boldsymbol{b}_{min}$

Calculate $D_{min}$ from $LIST_{min}$ and $D_{u,v}$

Algorithm 3 QUBO formulation of TSP on a boson sampler.

Input:

Number of locations $N$

Distance matrix $D_{u,v}$

Number of samples $N_{s}$

Output:

Bit string for which a minimum distance is found $\boldsymbol{b}_{min}$

Minimum problem Hamiltonian $H_{min}$

Minimum distance $D_{min}$

Cycle with minimum distance $LIST_{min}$

Calculate the bit string length $L_{3}$ as $(N-1)^{2}$

Configure the interferometer to produce bit strings $\boldsymbol{b}$ of length
$L_{3}$

for $i\leftarrow 1$ to $N_{s}$ do

Measure the output pattern generated by the interferometer.

Classical map the output pattern to a bit string $\boldsymbol{b}$ of length
$L_{3}$

Reshape $\boldsymbol{b}$ to a matrix with elements $x_{i,j}$

Calculate $H_{problem}$ classically using Equation 6 and $D_{u,v}$

if $H_{min}$ is initial then

$H_{min}\leftarrow H_{problem}$

else if $H_{problem}<H_{min}$ then

$H_{min}\leftarrow H_{problem}$

$\boldsymbol{b}_{min}\leftarrow\boldsymbol{b}$

end if

end for

Iterate changing the parameter vector $\boldsymbol{\theta}$ until convergence
is reached.

Calculate $LIST_{min}$ from $\boldsymbol{b}_{min}$

Calculate $D_{min}$ from $LIST_{min}$ and $D_{u,v}$