# Closed-form expressions for the magnetic field of permanent magnets in three
dimensions

V. Ziemann, Uppsala University

###### Abstract

We derive a closed-form expression of the magnetic field of a finite-size
current sheet and use it to calculate the field of permanent magnets, which
are modeled through their surface current densities. We illustrate the method
by determining the multipoles and the effective length due to fringe fields of
a finite-length dipole constructed of magnetic cubes.

##  1 Introduction

Today permanent magnets are widely used in the construction of efficient
electrical motors, for example those used in battery-driven automobiles, but
also for scientific applications, such as ion traps to store and analyze rare
atoms, or for undulators to generate synchrotron radiation, and for multipole
magnets to guide and confine charged-particle beams. These magnets do not need
to be powered and are therefore very energy efficient, which makes them
attractive regarding today’s focus on sustainability. In the 1970s Klaus
Halbach showed [1] how to design multipole magnets and undulators using paper
and pencil only. His approach made use of variables in the complex plane and
is essentially limited to two spatial dimensions. In order to analyze three-
dimensional features we typically use finite-element codes [2]. In this
report, we offer a complementary approach based on solving Biot-Savart’s
equation for a finite-size rectangular current sheet and modeling the magnets
as an assembly of such sheets. In iron-free geometries this allows us to
calculate the three components of the magnetic field $\vec{B}(z)$ at any point
in space as the superposition of the contributions from the current-carrying
sheets. Since all integrations are done analytically this results in
expressions that make calculating the fields very rapid. Moreover, the field
values are exact, no approximations due to finite meshes are involved.

A key feature of permanent magnet material is that the relative magnetic
permeability $\mu_{r}$ is very close to unity, which makes them transparent to
other, externally generated fields. As a consequence, the field at a given
point $\vec{z}$ is the superposition of the fields due to all current sources.
And here the second key feature enters; the magnetic field of permanent
magnets can be calculated from the equivalent surface current density
$I^{\prime}=dI/ds=B_{r}/\mu_{0}$ that is proportional to the remanent magnetic
field $B_{r}$ of the permanent magnet material. Here $\mu_{0}=4\pi\times
10^{-7}\,$A/Tm is the permeability of free space. Since $B_{r}$ is on the
order of one Tesla, the corresponding current densities are enormous.

These observations now guide us through the remainder of this report. We first
calculate the field of a straight current-carrying filament, followed by
spreading the filament laterally to obtain a current sheet and calculating the
field it generates. In the subsequent section we briefly introduce software,
based on MATLAB, to assemble magnets that are composed of these sheets and
calculate their fields. We then use the software to illustrate the field of a
single sheet, a magnet cube assembled of four sheets, and a dipole magnet,
made of four cubes, whose field quality and fringe fields we analyze. Finally,
we consider solenoid magnets which help us to validate the calculations and
the software.

##  2 Field of a filament

The magnetic field at point $\vec{r}_{2}$ due to a filament that carries a
current $I$, starts at position $\vec{r}_{a}$, and ends at $\vec{r}_{b}$ is
given by Biot-Savart’s law

$\vec{B}(\vec{r}_{2})=\frac{\mu_{0}I}{4\pi}\int\frac{d\vec{r}_{1}\times(\vec{r}_{2}-\vec{r}_{1})}{|\vec{r}_{2}-\vec{r}_{1}|^{3}}=\frac{\mu_{0}I}{4\pi}\int_{0}^{1}\frac{[\vec{r}_{b}-\vec{r}_{a}]\times[\vec{r}_{2}-t\vec{r}_{b}-(1-t)\vec{r}_{a}]}{|\vec{r}_{2}-t\vec{r}_{b}-(1-t)\vec{r}_{a})|^{3}}dt\
,$ (1)

where we parameterize a point on the filament $\vec{r}_{1}$ by
$\vec{r}_{1}=(1-t)\vec{r}_{a}+t\vec{r}_{b}$ with $0<t<1$. First we simplify
the numerator in the last equality and find

$(\vec{r}_{b}-\vec{r}_{a})\times(\vec{r}_{2}-t\vec{r}_{b}-(1-t)\vec{r}_{a})=(\vec{r}_{b}-\vec{r}_{a})\times(\vec{r}_{2}-\vec{r}_{a})$
(2)

whereas for the denominator we obtain

$\displaystyle|\vec{r}_{2}-t\vec{r}_{b}-(1-t)\vec{r}_{a})|^{2}$
$\displaystyle=$
$\displaystyle(\vec{r}_{2}-t\vec{r}_{b}-(1-t)\vec{r}_{a})\cdot(\vec{r}_{2}-t\vec{r}_{b}-(1-t)\vec{r}_{a})$
(3) $\displaystyle=$
$\displaystyle(\vec{r}_{2}-\vec{r}_{a})^{2}-2t(\vec{r}_{2}-\vec{r}_{a})\cdot(\vec{r}_{b}-\vec{r}_{a})+t^{2}(\vec{r}_{b}-\vec{r}_{a})^{2}$
$\displaystyle=$ $\displaystyle A-2Bt+Ct^{2}$

with

$A=(\vec{r}_{2}-\vec{r}_{a})^{2},\quad
B=(\vec{r}_{2}-\vec{r}_{a})\cdot(\vec{r}_{b}-\vec{r}_{a}),\quad\,\mathrm{and}\quad
C=(\vec{r}_{b}-\vec{r}_{a})^{2}\ .$ (4)

Note that these parameters have simple geometric interpretations: $A$
describes the squared distance $d_{P}$ between the observation point $P$ at
$\vec{r}_{2}$ and the starting point of the filament $\vec{r}_{a}$. The
variable $B$ is the product of $d_{P}$, the length $L$ of the filament, and
the cosine of the enclosed angle, whereas $C$ is the length of the filament
squared. Inserting these expressions into the equation for the magnetic field
$\vec{B}$ we get

$\vec{B}(\vec{r}_{2})=\frac{\mu_{0}I}{4\pi}(\vec{r}_{b}-\vec{r}_{a})\times(\vec{r}_{2}-\vec{r}_{a})\int_{0}^{1}\frac{dt}{(A-2Bt+Ct^{2})^{3/2}}\
.$ (5)

This integral can be found in integral tables [3] and is given by

$\int\frac{dt}{(A-2Bt+Ct^{2})^{3/2}}=\frac{Ct-B}{(AC-B^{2})\sqrt{A-2Bt+Ct^{2}}}$
(6)

such that we write $\vec{B}(\vec{r}_{2})$ as

$\vec{B}(\vec{r}_{2})=\frac{\mu_{0}I}{4\pi}\frac{(\vec{r}_{b}-\vec{r}_{a})\times(\vec{r}_{2}-\vec{r}_{a})}{AC-B^{2}}\left[\frac{C-B}{\sqrt{A-2B+C}}+\frac{B}{\sqrt{A}}\right]$
(7)

where $A,B$ and $C$ are given in Equation 4. Equation 7 is already quite
useful to determine the fields of systems containing a finite number of
current-carrying wires. We previously used it to design earth-field
compensation coils [4].

##  3 Field of a finite current sheet

Figure 1: Current sheet.

In order to calculate the field of a sheet, we first consider the field that
is generated by a filament that is laterally displaced by the vector

$\vec{p}=\alpha(\vec{r}_{b}-\vec{r}_{a})\times\vec{n}\ ,$ (8)

where $\vec{n}$ is the vector normal to the sheet and $\alpha$ is the aspect
ratio of the sheet, such that
$w=\alpha|\vec{r}_{b}-\vec{r}_{a}|=\alpha\sqrt{C}$ is the width of the sheet.
The construction is illustrated in Figure 1. Note that $\vec{p}$ is
perpendicular to the filament; the scalar product of $\vec{r}_{b}-\vec{r}_{a}$
and $\vec{p}$ is zero. This construction permits us to parameterize the
parallel-displaced filament by its end points $\vec{r}_{a}+s\vec{p}$ and
$\vec{r}_{b}+s\vec{p}$, where $1/2\leq s\leq 1/2$ describes by how far the
filament is displaced, as shown in Figure 1 by the thin dotted line that is
parallel to $\vec{r}_{b}-\vec{r}_{a}$. The field from the displaced filament
is then given by replacing the parameters $A,B$ and $C$ by their $s$-dependent
counterparts given by

$\displaystyle A(s)$ $\displaystyle=$
$\displaystyle(\vec{r}_{2}-\vec{r}_{a}-s\vec{p})^{2}$ $\displaystyle B(s)$
$\displaystyle=$
$\displaystyle(\vec{r}_{2}-\vec{r}_{a}-s\vec{p})\cdot(\vec{r}_{b}-\vec{r}_{a})$
(9) $\displaystyle C(s)$ $\displaystyle=$
$\displaystyle(\vec{r}_{b}+s\vec{p}-\vec{r}_{a}-s\vec{p})^{2}=C\ .$

Moreover, the cross-product before the integral becomes
$(\vec{r}_{b}-\vec{r}_{a})\times(\vec{r}_{2}-\vec{r}_{a}-s\vec{p})$. We now
“spread out” the current $I=I^{\prime}w$ evenly across the width
$w=\alpha\sqrt{C}$ of the sheet with constant current density $I^{\prime}$ and
sum up the contributions of all displaced filaments, which corresponds to
integrating over $s$. We also note that the current density $I^{\prime}$ is
related to the remanent magnetic field $B_{r}$ through $\mu I^{\prime}=B_{r}$.
Combining these concepts we find that the the field at point $\vec{r}_{2}$ is
given by the following integral

$\displaystyle\vec{B}(\vec{r}_{2})$ $\displaystyle=$
$\displaystyle\frac{B_{r}\alpha\sqrt{C}}{4\pi}\int_{-1/2}^{1/2}\frac{(\vec{r}_{b}-\vec{r}_{a})\times(\vec{r}_{2}-\vec{r}_{a}-s\vec{p})}{A(s)C-B(s)^{2}}$
$\displaystyle\qquad\qquad\quad\times\left[\frac{C-B(s)}{\sqrt{A(s)-2B(s)+C}}+\frac{B(s)}{\sqrt{A(s)}}\right]ds$

where $A(s)$ and $B(s)$ are at most quadratic functions of $s$. The evaluation
of the integral, which is rather lengthy and therefore deferred to Appendix A,
leads to

$\displaystyle\vec{B}(\vec{r}_{2})$ $\displaystyle=$
$\displaystyle\frac{B_{r}\alpha
C}{4\pi}\left[(C-B)J_{1}+BJ_{3}\right](\vec{r}_{b}-\vec{r}_{a})\times(\vec{r}_{2}-\vec{r}_{a})$
$\displaystyle\qquad+\frac{B_{r}\alpha
C^{2}}{4\pi}\left[(C-B)J_{2}+BJ_{4}\right]\vec{n}\ ,$

where $J_{1}$ to $J_{4}$ are defined in Appendix A in Equations 22 to 24. We
point out that $\alpha C$ is the area of the sheet, such that $B_{r}\alpha C$
describes the integrated “strength” of the sheet to excited magnetic fields.

##  4 Software implementation

Equation 3 describes the magnetic field vector $\vec{B}(\vec{r}_{2})$ that is
caused by one sheet that is characterized by
$\vec{r}_{a},\vec{r}_{b},\vec{n},\alpha,B_{r}$ in a coordinate-free way and
only requires function evaluations; no integrations are necessary which makes
the calculation of $\vec{B}(\vec{r}_{2})$ very fast. In order to calculate the
field at a given point $\vec{r}_{2}$ due to several sheets we only have to sum
over the contributions from each of the sheets.

We therefore represent each sheet by a one-dimensional array that contains the
three spatial components of $\vec{r}_{a},\vec{r}_{b},\vec{n}$ as well as
$\alpha$ and $B_{r}$ for that sheet; thus there are eleven numbers to store
for each sheet. The entire geometry is then given as a collection of, say $M$,
such sheets that is stored in an array of dimension $M\times 11$.

In order to place the sheets within the geometry, we prepare functions to
translate sheets by a vector $\vec{a}$; the function simply adds $\vec{a}$ to
$\vec{r}_{a}$ and $\vec{r}_{b}$ in the array with the sheets. Likewise, three
functions to rotate a sheet by an angle $\theta$ around each of the three
axes, just multiply $\vec{r}_{a},\vec{r}_{b}$, and $\vec{n}$ in a sheet by the
appropriate rotation matrix. Using these functions, it is easy to create any
geometry, very similar to the way that three-dimensional modeling software
composes complex geometries from geometric primitives, often triangles, but
also from cubes or spheres.

Finally functions to calculate the field at a point or along a path just loop
over the sheets and add their contributions. We implemented these functions in
MATLAB, briefly explain the functions in Appendix B, and illustrate their use
with a number of simple examples in the next section.

##  5 Examples

Figure 2: Field pattern of a current sheet (left) and a four sheets assembled
to generate the field of a permanent magnet cube (right).

As a first example we consider the field of a single sheet, which is shown on
the left-hand image in Figure 2. The sheet is shown as the blue rectangle and
the current as the green line stretching in the center of the sheet. The
vertical green arrow indicates the normal vector $\vec{n}$ of the sheet. We
observe that the magnetic field vector $\vec{B}$ curls around the sheet,
perpendicular to the current’s direction, pointing to the right above and to
the left below the sheet. Making the sheets very narrow causes the field to
approach a circular form, just as expected from Ampere’s law.

In the next step we assemble four quadratic current sheets to represent a
permanent magnet cube, which is shown on the right-hand image in Figure 2. We
see that the current, again shown as a green line, flows around the surface of
the cube. We use the convention that the normal vectors, again shown as green
arrows, point towards the outside of the magnet. The easy axis of the magnet
is indicated by the yellow arrow on the sheets, that is oriented perpendicular
to the current flow and points upwards. The magnetic field created by the four
sheets representing the permanent magnet cube is illustrated by the red
arrows. We observe that the field points upwards at the top of the cube, then
turns over such that it points down on the outside of the cube, then turns
over again to point upwards below the cube. We point out that the four current
sheets form a square solenoid and the field therefore corresponds to that
geometry.

Figure 3:  The geometry of a permanent magnet solenoid (top left) and the
vertical field in the midplane in the $x$-direction (top right). A 3D
rendition of the magnet (middle left) and the field along the $z$-axis (middle
right). The bottom left shows the field along the circle shown in cyan and the
multipole coefficients (absolute and relative) in the midplane. The bottom
right shows the corresponding integrated values.

We reach the next level of complexity by assembling eight 5 mm-cubes with
$B_{r}=1.47\,$T to form an approximation of a Halbach dipole, discussed in
Appendix A.3 in [5]. The top left image in Figure 3 shows the view onto the
assembly from the top. The easy axis of the magnets, shown by the yellow
arrows, rotates twice when moving from one magnet to the next; in every fourth
of the eight magnets it points in the same direction. In the center of the
assembly, the magnetic field points upwards, as is expected for a dipole. The
plot on the top right shows the vertical field component $B_{y}$ along the
black dotted line shown on the left-hand image. We observe that the field is
on the order of 0.35 T and has a small dip in the center. The image on the
left-hand side in the middle row shows a perspective view of the same geometry
which makes the three-dimensional character obvious. The vertical red line in
the center of the assembly indicates the direction of the $z$-axis and the
plot on its right shows $B_{y}$ along this direction, but extending over $\pm
20\,$mm. We see that the peak field of about 0.35 T appears in the center of
the magnet, but the fringe field extends significantly outside the magnet,
which only has a a height of $\pm 2.5\,$ mm.

We determine the multipolarity of the assembly at the radius $r$ shown by the
cyan circle on he top-left image. A tangent to the circle, illustrated by the
cyan straight line, is also shown. In order to determine the multipoles we
calculate the scalar product of the tangent vector with the magnetic field at
every point on the circle. This yields the tangential field component $B_{t}$
which we display in the top panel of the image at the bottom left plot as a
function of the azimuthal angle $\phi$. The multipolarity, shown on the middle
panel then follows from a Fourier-transformation, as discussed in Appendix C.
We find that mainly the dipole component with multipolarity $m=1$ is present
and that its amplitude is $B_{max}=0.353\,$T. Moreover, from the real and
imaginary part of the transform, we determine the roll angle of the magnet,
which here turns out to be close to zero. On the bottom panel we show the
amplitudes of the higher multipoles, normalized to $B_{max}$. We find the
decapole component ($m=5$) with a relative amplitude on the order of a percent
to give the dominant contribution. Finally, we calculate the integral of
$B_{t}$ along a line in the range $z=\pm 20\,$mm with $x$ and $y$-positions
determined by the cyan circle and display it on the upper panel in the image
on the bottom-right of Figure 3. Again, mostly the dipole component with $m=1$
is present. Moreover the field integral turns out to be $3.63\,$Tmm, which
allows us to determine the effective length of the magnet to be about 10.3 mm,
a value that significantly exceeds the physical length of the magnet, which is
only 5 mm. Again the decapole component with a relative magnitude of about
$6\times 10^{-3}$ gives the dominant contribution to the integrated field.

Figure 4:  The geometry of a permanent magnet solenoid (left) and the field on
axis (right).

As last example we consider an axial permanent magnet solenoid with
$B_{r}=1.47\,$T, inner radius 3 mm, and outer radius 10 mm. It is modeled by
circular currents flowing in one direction at the inner radius and in the
opposite direction at the outer radius. The left-hand image in Figure 4 shows
the assembly, here the circular magnet is approximated by a hollow polygon
with 18 sides. The field along the central axis of two long concentric
circular solenoids is known analytically [6]. We therefore compare it to our
numerical solution, based on Equation 3, and show both on the right-hand plot
in Figure 4. We observe that the agreement is very good.

##  6 Conclusions

We determined an analytic expression for the magnetic field generated by a
rectangular current sheet and used it to model permanent magnets in iron-free
geometries. We illustrated the method by calculating the fields of simple
geometries such as a single sheet, cubes, and solenoids assembled from
multiple sheets.

We expect this algorithm to be useful for rapid prototyping of three-
dimensional permanent-magnet structures. It is rather fast because all
integrations are done analytically and only fields at the points of interest
need to be calculated. This is in contrast to finite-element methods. They
first need to determine the fields on a grid in the whole geometry and then
interpolate the fields to the points of interest. On the other hand is the
present method limited to iron-free geometries and to magnets that can be
modeled with rectangular current sheets.

Discussions with Roger Ruber are gratefully acknowledged. In part, funding is
provided through the project Disseminating technology for cold magnets to
provide access to a wider international market that is supported by the
European Regional Development Fund (ERDF) and Region Kronoberg.

## References

  * [1] K. Halbach, Design of permanent multipole magnets with oriented rare earth cobald magnets, Nuclear Instruments and Methods 169 (1980) 1.
  * [2] Finite-element codes, for example, COMSOL https://www.comcol.com.
  * [3] I. Gradstein, I. Ryshik, Tables of series, products, and integrals, Verlag Harri Deutsch, Thun, 1981.
  * [4] V. Ziemann, R. Wedberg, T. Peterson, A. Wiren, Retro-Fitting Earth-Field Compensation Coils to the Vertical Cryostat GERSEMI in FREIA, Instruments 2020, 4(1), 8.
  * [5] V. Ziemann, Hands-On Accelerator Physics Using MATLAB, CRC Press, Boca Raton, 2019.
  * [6] Q. Peng, S. McMurry, J. Coey, Axial magnetic field produced by axially and radially magnetized permanent rings, Journal of Magnetism and Magnetic Materials 268 (2004) 165.

##  Appendix A The integral

In this appendix we evaluate the integral, given in Equation 3. First we
consider the cross product before the integral

$\displaystyle(\vec{r}_{b}-\vec{r}_{a})\times(\vec{r}_{2}-\vec{r}_{a}-s\vec{p})$
$\displaystyle\qquad\qquad=(\vec{r}_{b}-\vec{r}_{a})\times(\vec{r}_{2}-\vec{r}_{a})-\alpha
s(\vec{r}_{b}-\vec{r}_{a})\times\left[(\vec{r}_{b}-\vec{r}_{a})\times\vec{n}\right]$
$\displaystyle\qquad\qquad=(\vec{r}_{b}-\vec{r}_{a})\times(\vec{r}_{2}-\vec{r}_{a})+s\alpha
C\vec{n}$ (12)

where we used the identity
$\vec{a}\times[\vec{b}\times\vec{c}]=\vec{b}(\vec{a}\cdot\vec{c})-\vec{c}(\vec{a}\cdot\vec{b})$
to evaluate the double cross product in the second line. Next, we calculate

$A(s)=(\vec{r}_{2}-\vec{r}_{a})^{2}-2s(\vec{r}_{2}-\vec{r}_{a})\cdot\vec{p}+s^{2}\vec{p}^{2}=A-2Ds+s^{2}\vec{p}^{2}$
(13)

where we introduce $D=(\vec{r}_{2}-\vec{r}_{a})\cdot\vec{p}$. For
$\vec{p}^{2}$ in the last term we find

$\displaystyle\vec{p}^{2}$ $\displaystyle=$
$\displaystyle\alpha^{2}\left[(\vec{r}_{b}-\vec{r}_{a})\times\vec{n}\right]\cdot\left[(\vec{r}_{b}-\vec{r}_{a})\times\vec{n}\right]$
(14) $\displaystyle=$
$\displaystyle\alpha^{2}\left[(\vec{r}_{b}-\vec{r}_{a})^{2}\vec{n}^{2}-\left((\vec{r}_{b}-\vec{r}_{a})\cdot\vec{n}\right)^{2}\right]=\alpha^{2}C\
,$

which allows us to write

$A(s)=A-2Ds+\alpha^{2}Cs^{2}$ (15)

as a second-order polynomial in the variable $s$. For $B(s)$ we arrive at

$\displaystyle B(s)$ $\displaystyle=$
$\displaystyle(\vec{r}_{2}-\vec{r}_{a}-s\vec{p})\cdot(\vec{r}_{b}-\vec{r}_{a})$
$\displaystyle=$ $\displaystyle
B-s\alpha\left[(\vec{r}_{b}-\vec{r}_{a})\times\vec{n}\right]\cdot(\vec{r}_{b}-\vec{r}_{a})=B\
,$

which is constant. Inserting these parameters into Equation 3 leads to

$\displaystyle\vec{B}(\vec{r}_{2})$ $\displaystyle=$
$\displaystyle\frac{\mu_{0}I^{\prime}\alpha\sqrt{C}}{4\pi}\int_{-1/2}^{1/2}ds\frac{(\vec{r}_{b}-\vec{r}_{a})\times(\vec{r}_{2}-\vec{r}_{a})+\alpha
sC\vec{n}}{AC-B^{2}+C^{2}-2sCD+\alpha^{2}C^{2}s}$
$\displaystyle\left[\frac{\sqrt{C}(C-B)}{\sqrt{AC-2BC+C^{2}-2sCD+\alpha^{2}C^{2}s^{2}}}+\frac{\sqrt{C}(B)}{\sqrt{AC-2sCD+\alpha^{2}C^{2}s^{2}}}\right]\
.$

Close inspection of the integrand shows that the second-order polynomial in
the denominator of the first line only differs by a constant term from the
polynomials under the roots in the second line. Introducing the abbreviations

$R_{1}=a_{1}+b_{1}s+c_{1}s^{2}\quad\mathrm{with}\quad a_{1}=AC-2BC+C^{2},\
b_{1}=-2CD,\ c_{1}=\alpha^{2}C^{2}$ (18)

and

$R_{2}=a_{2}+b_{2}s+c_{2}s^{2}\quad\mathrm{with}\quad a_{2}=AC,\
b_{2}=-2CD=b_{1},\ c_{2}=\alpha^{2}C^{2}=c_{1}$ (19)

we can write

$\displaystyle\vec{B}(\vec{r}_{2})$ $\displaystyle=$
$\displaystyle\frac{\mu_{0}I^{\prime}\alpha
C}{4\pi}\int_{-1/2}^{1/2}ds\left[\frac{(\vec{r}_{b}-\vec{r}_{a})\times(\vec{r}_{2}-\vec{r}_{a})+\alpha
sC\vec{n}}{p_{1}+R_{1}}\frac{C-B}{\sqrt{R_{1}}}\right.$ (20)
$\displaystyle\qquad\qquad\qquad\quad\left.+\frac{(\vec{r}_{b}-\vec{r}_{a})\times(\vec{r}_{2}-\vec{r}_{a})+\alpha
sC\vec{n}}{p_{2}+R_{2}}\frac{B}{\sqrt{R_{2}}}\right]$ $\displaystyle=$
$\displaystyle\frac{\mu_{0}I^{\prime}\alpha
C}{4\pi}\left[(\vec{r}_{b}-\vec{r}_{a})\times(\vec{r}_{2}-\vec{r}_{a})(C-B)\int_{-1/2}^{1/2}\frac{ds}{(p_{1}+R_{1})\sqrt{R_{1}}}\right.$
$\displaystyle\qquad+\vec{n}\alpha
C(C-B)\int_{-1/2}^{1/2}\frac{sds}{(p_{1}+R_{1})\sqrt{R_{1}}}$
$\displaystyle\qquad+(\vec{r}_{b}-\vec{r}_{a})\times(\vec{r}_{2}-\vec{r}_{a})B\int_{-1/2}^{1/2}\frac{ds}{(p_{2}+R_{2})\sqrt{R_{2}}}$
$\displaystyle\qquad\left.+\vec{n}\alpha
CB\int_{-1/2}^{1/2}\frac{sds}{(p_{2}+R_{2})\sqrt{R_{2}}}\right]$

with $p_{1}=-(B-C)^{2}$ and $p_{2}=-B^{2}.$ We denote the integrals by

$\displaystyle J_{1}$ $\displaystyle=$
$\displaystyle\int_{-1/2}^{1/2}\frac{ds}{(p_{1}+R_{1})\sqrt{R_{1}}}\ ,\qquad
J_{2}=\int_{-1/2}^{1/2}\frac{sds}{(p_{1}+R_{1})\sqrt{R_{1}}}\ ,$
$\displaystyle J_{3}$ $\displaystyle=$
$\displaystyle\int_{-1/2}^{1/2}\frac{ds}{(p_{2}+R_{2})\sqrt{R_{2}}}\ ,\qquad
J_{4}=\int_{-1/2}^{1/2}\frac{sds}{(p_{2}+R_{2})\sqrt{R_{2}}}\ .$ (21)

We note that $J_{3}$ and $J_{4}$ resemble $J_{1}$ and $J_{2}$, only the
constant parameters $a_{i},b_{i},c_{i}$, and $p_{i}$ for $i=1,2$ differ. We
therefore only consider $J_{1}$ and $J_{2}$ in the following. As a matter of
fact, closed expressions for these integrals are given as number 2.284 in [3].
In particular $J_{1}$ and $J_{2}$ can be determined from

$UJ_{2}+VJ_{1}=\int\frac{Us+V}{(p_{1}+R_{1})\sqrt{R_{1}}}ds=\frac{U}{c_{1}}I_{1}(s)-\frac{2Vc_{1}-Ub_{1}}{\sqrt{c_{1}^{2}p_{1}\left[b_{1}^{2}-4(a_{1}+p_{1})c_{1}\right]}}I_{2}(s)$
(22)

with $I_{1}(s)$ given by [3]

$I_{1}(s)=\left\\{\begin{array}[]{ll}\frac{1}{\sqrt{p_{1}}}\arctan\left(\sqrt{\frac{R_{1}}{p_{1}}}\right)&\quad\mathrm{for}\
p_{1}>0\\\
\frac{1}{2\sqrt{-p_{1}}}\ln\left(\frac{\sqrt{-p_{1}}-\sqrt{R_{1}}}{\sqrt{-p_{1}}+\sqrt{R_{1}}}\right)&\quad\mathrm{for}\
p_{1}<0\ .\end{array}\right.$ (23)

With $d_{1}=p_{1}[b_{1}^{2}-4(a_{1}+p_{1})c_{1}]$ the second contribution
$I_{2}(s)$ is given by [3]

$I_{2}(s)=\left\\{\begin{array}[]{ll}-\arctan\left(\sqrt{\frac{p_{1}}{d_{1}}}\right)\frac{b_{1}+2c_{1}s}{\sqrt{R_{1}}}&\quad\mathrm{for}\
p_{1}d_{1}>0\\\
\frac{1}{2i}\ln\left(\frac{\sqrt{-d_{1}}\sqrt{R_{1}}+\sqrt{p_{1}}(b_{1}+2c_{1}s)}{\sqrt{-d_{1}}\sqrt{R_{1}}-\sqrt{p_{1}}(b_{1}+2c_{1}s)}\right)&\quad\mathrm{for}\
p_{1}d_{1}<0\ \mathrm{and}\ p_{1}>0\\\
\frac{1}{2i}\ln\left(\frac{\sqrt{d_{1}}\sqrt{R_{1}}+\sqrt{-p_{1}}(b_{1}+2c_{1}s)}{\sqrt{d_{1}}\sqrt{R_{1}}-\sqrt{-p_{1}}(b_{1}+2c_{1}s)}\right)&\quad\mathrm{for}\
p_{1}d_{1}<0\ \mathrm{and}\ p_{1}<0\\\ \end{array}\right.$ (24)

where $R_{1}$ is a function of $s$, given in Equation 18. We point out that
the sign of $I_{2}(s)$ in Equation 22 is reversed with respect to [3] to make
the numerically evaluated integral consistent with the value given in Equation
24. Inserting $J_{1}$ through $J_{4}$ in Equation 20 and reordering terms then
leads to the result stated in Equation 3 in the main body of this report.

##  Appendix B MATLAB functions

The following MATLAB functions implement the functionality described in
Section 4 and are used to generate the plots shown in the main body of this
report. Note that sheets is a $M\times 11$ array that describes $M$ sheets.

  * •

B=Bsheets(sheets,r2): returns $\vec{B}(\vec{r}_{2})$ for the sheets defined by
the input argument.

  * •

sheets=make_cubez(a,alpha,Br): returns the four sheets for a square block with
width $a$ and height $\alpha a$ with the easy axis pointing in the
$z$-direction.

  * •

sheets=make_brick(a,b,c,Br): returns the four sheets for a rectangular block
with sizes $a,b,$ and $c$ in the $x,y$, and $z$ direction, respectively. The
easy axis points in the $z$-direction.

  * •

sheets=make_polygon(n,r,h,Br): returns the $n$ sheets for a polygon with $n$
sides, where $r$ specifies the distance from the center to the corners and $h$
is the height in the $z$ direction. The easy axis also points in the
$z$-direction.

  * •

sheets=make_polygon_hollow(n,ri,ro,h,Br): returns the $2n$ sheets for a hollow
polygon with $n$ sides, where $r_{i}$ specifies the distance from the center
to the corners of the inner polygon and $r_{o}$ to the outer polygon; $h$ is
the height in the $z$ direction. The easy axis also points in the
$z$-direction. This function creates “solenoidal” permanent magnets with
single call.

  * •

sheets=sheets_translate(sheets,a): returns the sheets provided as input after
adding a vector $\vec{a}$ to $\vec{r}_{a}$ and $\vec{r}_{b}$ of the sheets.

  * •

sheets=sheets_rotate_x(sheets,theta): returns the sheets provided as input
after the corresponding vectors $\vec{r}_{a},\vec{r}_{b}$, and $\vec{n}$ of
the sheets are rotated by the angle $\theta$. There are corresponding
functions for rotations around the $y$ and $z$-axis.

  * •

draw_sheet(sheet): draws one sheet as a blue surface and adds a green line to
indicate the direction of the current flow. A yellow arrow indicates the easy
axis.

  * •

B=field_along_line(sheets,line): returns an array of $\vec{B}(\vec{r}_{2})$ at
the sequence of points $\vec{r}_{2}$ provided in the array line.

The use of these functions is best illustrated in the code to prepare a
Halbach multipole. The code below first prepares a 5 mm cube and the rotates
it such that the easy axis points in the vertical direction. After defining
the number of cubes and the multipolarity we loop over the cubes. Inside the
loop we first rotate the cube as needed to generate the desired multipolarity
$m$, then translate it radially along the $x$ axis by 8 mm before rotating
into place a second time. Last, we add the sheets of the currently handled
cube to the sheets which, after the loop completes, contains all the sheets
that define the multipole.

      cube=make_cubez(5,1,Br);
      cubey=sheets_rotate_x(cube,-90);
      MM=8; % number of cubes
      m=1;  % multipolarity (m=1: dipole, 2: quadrupole)
      for k=0:MM-1
        tmp=sheets_rotate_z(cubey,-k*m*360/MM);
        tmp=sheets_translate(tmp,[8;0;0]);
        tmp=sheets_rotate_z(tmp,-k*360/MM);
        if k==0
          sheets=tmp;
        else
          sheets=[sheets;tmp];
        end
      end

##  Appendix C Determining the multipolarity

On magnet test benches the multipolarity of magnets is often determined from
the voltage that is induced in a coil that rotates inside the magnet. In a
dipole the polarity changes once per revolution and the voltage shows a
spectral component at the revolution frequency. In a quadrupole, the polarity
changes twice per revolution and the spectrum shows a peak at the second
revolution harmonic, whose amplitude is proportional to the quadrupole
gradient.

Here, we use a simplified version of this procedure and calculate the
projection of the magnetic field onto a (normalized) vector that is tangent to
a circle with radius $r$ inside the magnet. This is illustrated by the cyan
circle and tangent vector visible on the top-left image in Figure 3. In much
the same way as on a measurement bench, the number of polarity changes when
traveling around the circle is related to the multipolarity $m$ of the magnet.
The amplitude $\hat{B}$ of the Fourier harmonic directly gives the magnetic
field harmonic on the circle. For example, in a quadrupole, $\hat{B}$ is
related to the gradient $g$ by $g=\hat{B}/r$.

Moreover, if we denote the real and imaginary part of the largest Fourier
harmonic by $a_{m}+ib_{m}$, the roll angle $\hat{\phi}$ of the magnet around
the $z$-axis (pointing upwards on the top-left image in Figure 3) is given by
$\hat{\phi}=(1/m)\arctan(b_{m}/a_{m})$. In the legend of the images in the
bottom row this value is also reported. It is, however, very small, because
the magnet was not rotated.

If we use the field on the cyan circle, we measure the harmonic at a given
$z$-position, similar to what a short rotating coil would do. If, on the other
hand, we integrate the field along a line extending in the $z$-direction and
perpendicular to positions on the circle, we determine the harmonics of the
integrated field, similar to what a long coil would do. This information is
reported on the bottom left image in Figure 3.


# Self-Tuning Tube-based Model Predictive Control

Damianos Tranos, Alessio Russo, and Alexandre Proutiere This work was
supported by the Wallenberg AI, Autonomous Systems and Software Program (WASP)
funded by the Knut and Alice Wallenberg Foundation.D. Tranos, A. Russo, and A.
Proutiere are with the Division of Decision and Control Systems, School of
Electrical Engineering and Computer Science, Royal Institute of Technology
(KTH), Stockholm, Sweden. Emails<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>}.

###### Abstract

We present Self-Tuning Tube-based Model Predictive Control (STT-MPC), an
adaptive robust control algorithm for uncertain linear systems with additive
disturbances based on the least-squares estimator and polytopic tubes. Our
algorithm leverages concentration results to bound the system uncertainty set
with prescribed confidence, and guarantees robust constraint satisfaction for
this set, along with recursive feasibility and input-to-state stability.
Persistence of excitation is ensured without compromising the algorithm’s
asymptotic performance or increasing its computational complexity. We
demonstrate the performance of our algorithm using numerical experiments.

##  I Introduction

Model Predictive Control (MPC) [1] addresses the infinite horizon optimal
control problem in the presence of input and state constraints by
approximating it as a sequence of finite horizon optimization problems. When
the dynamics of the system are uncertain, robust MPC methods [2] can be
employed to ensure constraint satisfaction for pre-specified sets of system
parameters and disturbances. This robustness comes at the cost of a reduced
closed-loop performance. To mitigate this performance loss, one may leverage
adaptive control techniques [3] to learn the system dynamics in an online
manner, and in turn reduce the uncertainty causing this loss.

Adaptive control schemes have been developed and studied mainly for
unconstrained control problems. The first results [4] concerned their
asymptotic convergence properties and established conditions under which the
controller derived from these schemes actually approaches the optimal feedback
controller obtained by assuming the full knowledge on the system dynamics.
More recently, see e.g. [5, 6, 7], researchers managed to quantify the
convergence rate of some classical adaptive control schemes, such the
celebrated self-tuning regulators, as well as the price that one has to pay
(in terms of cumulative losses – often captured through the notion of regret)
to learn the system dynamics. These recent important results are however
restricted to the LQR problems in unconstrained linear systems.

In this paper, we investigate the design and the performance analysis of MPC
schemes handling system dynamics uncertainties. We propose to combine adaptive
and robust control methods. The adaptive control component of our schemes
allows to rapidly reduce over time the uncertainties in a controlled and
quantifiable manner, whereas the robust control component ensures constraint
satisfaction. More precisely, our contributions are as follows.

Contributions. We address the problem of controlling an uncertain linear
system with parametric and additive disturbances, subject to deterministic
constraints. We present Self-Tuning Tube-based Model Predictive Control (STT-
MPC), an algorithm combining adaptive and robust control techniques. STT-MPC
uses a simple Least Squares Estimator (LSE) for estimating the dynamics and a
parameter set compatible with the observations with prescribed level of
certainty. We derive tight concentration results (similar to those presented
in [7] for the LQR problem) for this set, and exploit these results to
construct fixed complexity polyhedral approximations of it. In STT-MPC, these
approximations are used to build a polytopic tube MPC scheme [8] to ensure
robust constraint satisfaction. We establish the recursive feasibility and the
input-to-state stability of the proposed scheme.

In contrast to previously proposed robust MPC schemes (refer to §II for
details), STT-MPC enjoys the following properties. (i) STT-MPC uses a
probabilistic rather than robust estimation scheme, leading to notably faster
adaptation rates and smaller parameter sets for which the constraints need to
be satisfied. (ii) Persistence of Excitation (PE), required to get performance
guarantees for the LSE, is achieved without modifying either the cost or
objective function of the MPC. The added excitation is treated as an additive
disturbance and can be chosen to decay to zero with time (e.g., at a rate
$1/\sqrt{t}$). This allows us to asymptotically recover the performance of the
standard MPC algorithm with full knowledge of the dynamics. (iii) An input-to-
state stability analysis of STT-MPC is possible even when the PE condition is
active and without imposing any additional restrictions on the choice of the
estimate of the nominal model parameter.

Notations. For a time dependent vector $x_{t}$, we denote by $x_{k|t}$ its
prediction at time $k+t$ given information at time $t$. For any two sets
$\mathcal{A}$ and $\mathcal{B}$, we define their Minkowski sum as the set
$\mathcal{A}\oplus\mathcal{B}:=\\{a+b:a\in\mathcal{A},\ b\in\mathcal{B}\\}$.
We also define, for any constant $\lambda\geq 0$, the scaled set
$\lambda\mathcal{A}:=\\{\lambda a,\ a\in\mathcal{A}\\}.$ For any
$d\in{\mathbb{N}}$, $x\in{\mathbb{R}}^{d}$, and $\epsilon>0$ let
$\mathcal{B}(x,\epsilon)$ denote the $\epsilon$ ball of the spectral norm
centered on $x$. The unit ball centered at the origin is denoted simply as
$\mathcal{B}$. For any set $\mathcal{S}$, and any $\varepsilon>0$ there exists
a polytope $\mathcal{P}$ that is an outer approximation of $S$, i.e.
$\mathcal{S}\subset\mathcal{P}\oplus\varepsilon\mathcal{B}$. We refer to this
polytope as the outer polyhedral approximation of $\mathcal{S}$. Finally, a
function $\kappa:{\mathbb{R}}_{+}\to{\mathbb{R}}_{+}$ is of class
$\mathcal{K}$ if it is strictly increasing and $\kappa(0)=0$ and is of class
$\mathcal{K}_{\infty}$ if in addition $\kappa(x)\to\infty$ as $x\to\infty$.

##  II Related Work

An essential ingredient of adaptive control is Persistent Excitation (PE) [9]
(which is equivalent to the notion of required exploration in reinforcement
learning [10]). In the context of MPC, different solutions for PE have been
investigated: [11] introduces a sufficient PE condition as an additional
constraint in their MPC optimization. On the other hand, [12], [13], and later
[14], propose dual-MPC schemes where the predicted parameter covariance matrix
is included in the cost function, though these dual-control methods lack
feasibility and stability guarantees.

In a different direction, the works of [15] and [16] introduce robust MPC
schemes with guaranteed recursive feasibility and constraint satisfaction.
These schemes come with the cost of only considering adaptation for the
nominal model, where PE is not ensured, and robustness is guaranteed for a
fixed parameter set that is not updated online. In contrast, our scheme, STT-
MPC, updates the parameter set each step, and the latter rapidly concentrates
around the true paraemter.

In [17], the authors combine set-membership identification, which involves
updating a set of parameters compatible with the observed state-trajectory,
and robust constraint tightening for FIR models. This approach was later
extended in [18] to the case of linear state-space models where online set-
membership is combined with homothetic tube MPC [19]. Subsequently, [20]
proposed an adaptive robust MPC scheme using again set-membership
identification and a less conservative polytopic tube MPC scheme [8]. As in
[11], they ensure PE in the form of added convex constraints and a
modification of the cost function.

##  III Model, Objective, and Approach

###  III-A Model and assumptions

We consider the following discrete time, linear, time-invariant system:

$\displaystyle x_{t+1}=A(\theta^{\star})x_{t}+B(\theta^{\star})u_{t}+w_{t},$
(1)

where $x_{t},w_{t}\in\mathbb{R}^{d_{x}}$ and $u_{t}\in\mathbb{R}^{d_{u}}$. The
state transition and state-action transition matrices $A(\theta^{\star})$ and
$B(\theta^{\star})$ are initially unknown. The set of possible such matrices
is parametrized by $\theta\in\mathbb{R}^{d_{\theta}}$ (here $\theta$ could
well parametrize each entry of the matrices, in which case
$d_{\theta}=d_{x}(d_{x}+d_{u})$). To simplify the notations, for two possible
parameters $\theta_{1},\theta_{2}$, we define
$\lVert\theta_{1}-\theta_{2}\rVert:=\max(\lVert
A(\theta_{1})-A(\theta_{2})\rVert_{2},\lVert
B(\theta_{1})-B(\theta_{2})\rVert_{2})$. We make the following assumptions.

######  Assumption 1 (Parameter uncertainty).

For some $\epsilon_{0}>0$, ${\cal
B}(\theta^{*},\epsilon_{0})\subset\Theta_{0}$ where $\Theta_{0}$ is a known
convex polytope.

######  Assumption 2 (Additive disturbance).

The sequence $(w_{t})_{t\geq 0}$ is i.i.d, and for each $t\geq 0$, $w_{t}$ is
zero-mean, isotropic, with support in the ball $\mathcal{B}(0,3\sigma)$.
Hence, $w_{t}$ is $\sigma^{2}$-sub-gaussian. Further define $\mathcal{W}$, a
convex polytope providing a conservative approximation of
$\mathcal{B}(0,3\sigma)$, i.e., $\mathcal{B}(0,3\sigma)\subset\mathcal{W}$.

The system also needs to obey the following state and input constraints for
all $t\geq 0$,

$\displaystyle Fx_{t}+Gu_{t}\leq\textbf{1},$ (2)

where $F\in{\mathbb{R}}^{d_{c}\times d_{x}}$ and
$G\in{\mathbb{R}}^{d_{c}\times d_{u}}$ define the state and input constraints
respectively.

######  Assumption 3 (State and input constraints).

The set

$\displaystyle\mathcal{C}=\\{(x,u)\in{\mathbb{R}}^{d_{x}}\times{\mathbb{R}}^{d_{u}}:Fx+Gu\leq\textbf{1}.\\}$

is compact and contains the origin in its interior.

Finally, we assume that we have access to a stabilizer $K$:

######  Assumption 4 (Stabilizing Controller).

There exists a known, robustly stabilizing feedback gain $K$ such that
$A(\theta)+B(\theta)K$ is stable (i.e., $\rho(A(\theta)+B(\theta)K)<1$) for
all $\theta\in\Theta_{0}$.

###  III-B Objective and MPC

We wish to minimize the long-term cost defined as
$\lim\sup_{T\to\infty}\frac{1}{T}\sum_{t=0}^{T-1}x_{t}^{\top}Qx_{t}+u_{t}^{\top}Ru_{t}$,
through some positive semi-definite matrices $Q,R$. To this aim, we use MPC,
with a receding horizon $N$. Specifically, at time $t$, given the current
system state $x_{t}$ and the past observations used to derive an estimator
$\theta_{t}$ of $\theta^{\star}$, we will identify a control policy
$(u_{k|t})_{k=0,\ldots,N-1}$ minimizing the cost along a predicted system
trajectory $(x_{k|t})_{k=0,\ldots,N}$. We use the well-known dual mode
prediction paradigm [21] with the following predicted control sequence,

$\displaystyle u_{k|t}=\begin{cases}Kx_{k|t}+v_{k|t}\quad\forall
k\in\\{0,\dots,N-1\\},\\\ Kx_{k|t}\quad\forall k\geq N,\end{cases}$ (3)

where $\\{v_{0|t},\dots v_{N-1|t}\\}$ are the optimization variables to be
determined by the MPC. The resulting prediction dynamics will be for
$k\in\\{0,\dots,N-1\\}$,

$\displaystyle x_{0|t}$ $\displaystyle=x_{t},$ (4a) $\displaystyle x_{k+1|t}$
$\displaystyle=\Phi(\theta_{t})x_{k|t}+B(\theta_{t})v_{k|t},$ (4b)

where $\Phi(\theta):=A(\theta)+B(\theta)K$ for any $\theta\in\Theta_{0}$.

###  III-C General approach

To handle the uncertainty due to both the noise and the fact that
$\theta^{\star}$ is unknown, we apply a tube-based MPC approach. The tube used
in step $t$ is essentially constructed from a polytope approximating a ball
centered at the LSE $\hat{\theta}_{t}$ of $\theta^{\star}$ and whose radius
corresponds to a prescribed level of confidence $\delta$ that the user wishes
to guarantee. Persistent excitation is achieved by adding (bounded) noise to
the input. The resulting algorithm is presented in Algorithm 1, and its
ingredients are detailed in the next section. Its analysis is given in Section
V.

Input: Initial state $x_{0}$; confidence $\delta$; estimate $\theta_{0}$;
uncertain parameter set $\Theta_{0}$

1 Find stabilizing matrix $K$ for all $\theta\in\Theta_{0}$

2 Compute matrices $T$ and $H_{c}$ using (10b)

3 for  _$t=1,\dots,\mathcal{T}$_ do

4 if  _$t <t^{\star}(\delta)$_ then

5 Set $\theta_{t}\leftarrow\theta_{t-1}$ and
$\Theta_{t}\leftarrow\Theta_{t-1}$

6

7 else

8 Update $\theta_{t}\leftarrow\hat{\theta}_{t}$ using (5)

9 Compute $\Delta_{t}$ and set
$\Theta_{t}\leftarrow\Theta_{t-1}\cap\Delta_{t}$

10

11 Compute matrices $H_{t}^{(j)}$ using (10a)

12 Obtain $P(\theta_{t})$ by solving (20)

13 Solve
$\mathcal{P}_{N}(x_{t},\theta_{\rho(t)},\Theta_{\rho(t)},\bar{\mathcal{W}}_{t})$

14 Apply $u_{t}=Kx_{t}+v^{*}_{0|t}+\zeta_{t}$

15

Algorithm 1 STT-MPC

##  IV Self-tuning Tube-based MPC

###  IV-A LSE and persistent excitation

Our algorithm starts with an initial parameter $\theta_{0}\in\Theta_{0}$,
which is then updated using the LSE. We let $\theta_{1}=\theta_{0}$. For
$t\geq 2$, the LSE enjoys the following explicit expression:

$\displaystyle\hat{\theta}_{t}={\left(\sum\limits_{k=0}^{t-2}x_{k+1}\begin{bmatrix}x_{k}\\\
u_{k}\end{bmatrix}^{\top}\right)\left(\sum\limits_{k=0}^{t-2}\begin{bmatrix}x_{k}\\\
u_{k}\end{bmatrix}\begin{bmatrix}x_{k}\\\
u_{k}\end{bmatrix}^{\top}\right)^{\dagger}}.$ (5)

For $t\geq t^{\star}(\delta)$, we align our prediction parameter $\theta_{t}$
to $\hat{\theta}_{t}$. As shown in [7], a finite-time analysis of the
performance of the LSE is rather intricate but possible even if the feedback
controller varies over time. The performance is tightly related to the minimal
eigenvalue of the cumulative covariate matrix
$\lambda_{\min}(\sum_{s=0}^{t-2}y_{s}y_{s}^{\top})$ where
$y_{s}=\begin{bmatrix}x_{s}\\\ u_{s}\end{bmatrix}$. More precisely, for the
LSE to lead to a good approximation of $\theta^{\star}$, we need to ensure
that this eigenvalue grows with time. To this aim, we add an isotropic and
bounded noise to the control input. This noise is represented by the random
vector $\zeta_{t}$ taken to be the projection of $\xi_{t}$ on
$\mathcal{B}(0,3\sigma_{t})$, where $\xi_{t}$ is i.i.d. according to a normal
distribution, i.e., $\xi_{t}\sim\mathcal{N}(0,\sigma_{t}^{2}I_{d_{u}})$. The
choice of $\sigma_{t}$ directly impacts the performance of the LSE but also
the overall performance of the controller (a higher $\sigma_{t}$ means higher
excitation and hence better LSE, but at the expense of a worse overall
controller, refer to Section VII). Here we set
$\sigma_{t}^{2}=\sqrt{d_{x}}\sigma^{2}t^{-\alpha}$ for some $\alpha\in(0,1)$.
This ensures that (i) the LSE $\theta_{t}$ converges to $\theta^{\star}$
(refer to Section V for a precise statement) and (ii) the controller converges
to that obtained through a classical tube-based MPC framework with known
$\theta^{\star}$. When analyzing the performance of the LSE, we will show that
the good event ${\cal G}$ holds with probability at least $1-\delta$:

${\cal G}=\left(\|\hat{\theta}_{t}-\theta^{\star}\|\leq\epsilon_{t},\ \forall
t\geq t^{\star}(\delta)\right),$

where $t^{\star}(\delta)=c_{1}+c_{2}\log(1/\delta)$ and
$\epsilon_{t}^{2}=c_{3}\log(t/\delta)/t^{1-\alpha}$ for some positive
constants $c_{1},c_{2},c_{3}$. From the above result, we define $\Delta_{t}$
as an outer polyhedral approximation of
$\mathcal{B}(\theta_{t},2\epsilon_{t})$. We further recursively define the
uncertainty sets as follows: $\Theta_{t}=\Theta_{t-1}\cap\Delta_{t}$ for all
$t\geq t^{\star}(\delta)$ and $\Theta_{t}=\Theta_{0}$ for
$t<t^{\star}(\delta)$. By construction, the true parameter $\theta^{\star}$
belongs to the interior of $\Theta_{t}$ with high probability in the following
sense:

######  Lemma 1.

Under event ${\cal G}$,
$\mathcal{B}(\theta^{\star},\epsilon_{t})\subset\Theta_{t}$ for all $t\geq 1$.

###### Proof.

For $t<t^{\star}(\delta)$ we have $\Theta_{t}=\Theta_{0}$ and the result holds
by Assumption 1. Let $t\geq t^{\star}(\delta)$. We show that
$\mathcal{B}(\theta^{*},\epsilon_{t})\subset\Delta_{t}$ and
$\mathcal{B}(\theta^{*},\epsilon_{t})\subset\Theta_{t-1}$. For the first, we
have for all $\theta\in\mathcal{B}(\theta^{*},\epsilon_{t})$, using the
triangle inequality:
$\lVert\theta-\theta_{t}\rVert\leq\lVert\theta-\theta^{*}\rVert+\lVert\theta_{t}-\theta^{*}\rVert\leq
2\epsilon_{t}$, where the second inequality holds under ${\cal G}$. Thus
$\theta\in\mathcal{B}(\theta_{t},2\epsilon_{t})$ and so
$\mathcal{B}(\theta^{*},\epsilon_{t})\subset\mathcal{B}(\theta_{t},2\epsilon_{t})\subset\Delta_{t}$.
We prove $\mathcal{B}(\theta^{*},\epsilon_{t})\subset\Theta_{t-1}$ by
induction. Assume that
$\mathcal{B}(\theta^{*},\epsilon_{t})\subset\Theta_{t-1}$. Then we show that
$\mathcal{B}(\theta^{*},\epsilon_{t+1})\subset\Theta_{t}$. Let $\theta$ such
that $\|\theta-\theta^{\star}\|\leq\epsilon_{t+1}$. Then
$\|\theta-\theta_{t}\|\leq\lVert\theta-\theta^{*}\rVert+\lVert\theta_{t}-\theta^{*}\rVert\leq\epsilon_{t+1}+\epsilon_{t}$.
This implies that $\theta\in\Delta_{t}$. In addition, $\theta\in\Theta_{t-1}$.
Indeed, $\|\theta-\theta^{\star}\|\leq\epsilon_{t+1}\leq\epsilon_{t}$, and we
conclude using the induction assumption. Hence $\theta\in\Theta_{t}$. $\Box$

###  IV-B Polytopic tubes and associated linear constraints

With the considered control inputs, the system dynamics becomes:

$x_{t+1}=\Phi(\theta^{\star})x_{t}+B(\theta^{\star})(v_{0|t}+\zeta_{t})+w_{t}.$
(6)

We apply a tube-based approach, and at time $t$, we build a polytopic tube
based on:

  * (i)

$\Theta_{t}$, encoding the uncertainty about $\theta^{\star}$. We denote by
$m_{t}$ the number of vertices of $\Theta_{t}$, and the vertices themselves by
$\theta_{t}^{(j)}$, $j=1,\ldots,m_{t}$.

  * (ii)

A polytope $\bar{\mathcal{W}}_{t}$, handling the uncertainty due to the noise
$B(\theta^{\star})\zeta_{t}+w_{t}$, including that due to the persistent
excitation. To define $\bar{\mathcal{W}}_{t}$, let
$\bar{B}_{t}=\max_{\theta\in\Theta_{t}}\lVert B(\theta)\rVert_{2}$. We define
$\mathcal{Z}_{t}$ as the outer polyhedral approximation to
$\mathcal{B}(0,3\sigma_{t})$, and
$\bar{\mathcal{W}}_{t}=\mathcal{W}\oplus\bar{B}_{t}\mathcal{Z}_{t}$.

We define the state tube cross sections as the sets: for $k=0,\ldots,N$,

$\displaystyle X_{k|t}=\\{x:Tx\leq\alpha_{k|t}\\},$ (7)

where the matrix $T\in{\mathbb{R}}^{d_{\alpha}\times d_{x}}$ is chosen such
that, for some $\lambda\in[0,1)$, the set $\\{x:Tx\leq\mathbf{1}\\}$ is
$\lambda$-contractive with respect to the system $x_{t+1}=\Phi(\theta)x_{t}$
for all $\theta\in\Theta_{0}$. This property is needed to ensure the robust
positive invariance of $X_{N|t}$ (see Lemma 5.7 in [21]). To derive the
associated linear constraints, we apply a standard result to ensure inclusion
of polyhedra (see Proposition 3.31 in [22]). More precisely, for any
$j=1,\ldots,m_{t}$, we have
$X_{k|t}\subseteq\\{x:\Phi(\theta_{t}^{(j)})x+B(\theta_{t}^{(j)})v_{k|t}+w\in
X_{k+1|t}\\}$ for all $w\in\bar{\mathcal{W}}_{t}$ if there exists
$H_{t}^{(j)}\geq 0$ such that:

$\displaystyle H_{t}^{(j)}T=T\Phi(\theta_{t}^{(j)}),$ (8a) $\displaystyle
H_{t}^{(j)}\alpha_{k|t}+TB(\theta_{t}^{(j)})v_{k|t}+\bar{w}_{t}\leq\alpha_{k+1|t},$
(8b)

where $\bar{w}_{t}$ is such that
$[\bar{w}_{t}]_{i}=\max_{w\in\bar{\mathcal{W}}_{t}}[Tw]_{i}$ for
$i\in\\{1,\dots,d_{\alpha}\\}$. Similarly, we have
$X_{k|t}\subseteq\\{x:(F+GK)x+G(v_{k|t}+\zeta)\leq\textbf{1}\\}$ for all
$\zeta\in\mathcal{Z}_{t}$, if there exists $H_{c}\geq 0$, such that:

$\displaystyle H_{c}T=F+GK$ (9a) $\displaystyle
H_{c}\alpha_{k|t}+Gv_{k|t}+\bar{\zeta}_{t}\leq\textbf{1},$ (9b)

where $\bar{\zeta}_{t}$ is such that
$[\bar{\zeta}_{t}]_{i}=\max_{\zeta\in\mathcal{Z}_{t}}[G\zeta]_{i}$ for
$i\in\\{1,\dots,d_{c}\\}$.

######  Remark 1.

For simplicity, we have chosen $\bar{w}_{t}$ and $\bar{\zeta}_{t}$ so as to be
conservative with respect to the noise we apply. However, notice that we could
well pre-sample the noise sequence $(\zeta_{t})_{t\geq 0}$ and use it as part
of our predictions.

The matrices $H_{t}^{(j)}$, and $H_{c}$ are chosen such that for all
$i\in\\{1,\dots,d_{\alpha}\\}$ and for all $j\in\\{1,\dots,m_{t}\\}$,

$\displaystyle\big{(}H_{t}^{(j)}\big{)}_{i}$
$\displaystyle=\arg\min_{h}\\{\textbf{1}^{\top}h:h^{\top}T=(T)_{i}\Phi(\theta_{t}^{(j)}),\
h\geq 0\\},$ (10a) and for all $i\in\\{1,\dots,d_{c}\\}$,
$\displaystyle\big{(}H_{c}\big{)}_{i}$
$\displaystyle=\arg\min_{h}\\{\textbf{1}^{\top}h:h^{\top}T=(F+GK)_{i},\ h\geq
0\\}.$ (10b)

Finally, we have the terminal conditions

$\displaystyle H_{t}^{(j)}\alpha_{N}+\bar{w}_{t}$
$\displaystyle\leq\alpha_{N},$ (11) $\displaystyle
H_{c}\alpha_{N}+\bar{\zeta}_{t}$ $\displaystyle\leq\mathbf{1}.$ (12)

###  IV-C The tube MPC problem

Fix a step $t$. Let $\boldsymbol{v}_{t}=\\{v_{0|t},\dots,v_{N-1|t}\\}$ and
$\boldsymbol{\alpha}_{t}=\\{\alpha_{0|t},\dots,\alpha_{N|t}\\}$. The resulting
tube MPC problem, denoted as
$\mathcal{P}_{N}(x_{t},\theta_{t},\Theta_{t},\bar{\mathcal{W}}_{t})$, and
evaluated at every $t$, is

$\displaystyle\underset{\boldsymbol{v}_{t},\boldsymbol{\alpha}_{t}}{\textnormal{minimize}}\sum_{k=0}^{N-1}\left(\
x_{k|t}^{\top}Qx_{k|t}+v_{k|t}^{\top}Rv_{k|t}\right)+x_{N|t}^{\top}P(\theta_{t})x_{N|t}$
$\displaystyle\textnormal{subject to},\textnormal{for all
}j=1,\dots,m_{t},\textnormal{and }k=0,\dots,N-1:$ initial constraints:
$\displaystyle\quad x_{0|t}=x_{t},$ (13) $\displaystyle\quad
Tx_{0|t}\leq\alpha_{0|t},$ (14) $\displaystyle\textnormal{system
constraints},$ $\displaystyle\quad
x_{k+1|t}=\Phi(\theta_{t})x_{k|t}+B(\theta_{t})v_{k|t},$ (15)
$\displaystyle\textnormal{tube constraints},$ $\displaystyle\quad
H_{t}^{(j)}\alpha_{k|t}+TBv_{k|t}+\bar{w}_{t}\leq\alpha_{k+1|t},$ (16)
$\displaystyle\quad H_{c}\alpha_{k|t}+Gv_{k|t}+\bar{\zeta}_{t}\leq\mathbf{1},$
(17) terminal constraints: $\displaystyle\quad
H_{t}^{(j)}\alpha_{N|t}+\bar{w}_{t}\leq\alpha_{N|t},$ (18) $\displaystyle\quad
H_{c}\alpha_{N|t}+\bar{\zeta}_{t}\leq\mathbf{1}.$ (19)

where $P(\theta_{t})$ is obtained by solving the Lyapunov equation:

$\displaystyle
P(\theta_{t})-\Phi(\theta_{t})^{\top}P(\theta_{t})\Phi(\theta_{t})=Q+K^{\top}RK$
(20)

The feasibility of the above problem depends on whether
$\theta^{\star}\in\Theta_{t}$, i.e., on the event ${\cal G}$. However, for the
sake of the analysis, we would like to ensure that the problem the algorithm
solves in each step is always feasible (with probability 1). In turn, this
will allows us to analyze the performance of the LSE and to establish that
indeed ${\cal G}$ holds with probability at least $1-\delta$ (refer to the
next section for details). To ensure that the tube MPC problem is always
feasible, i.e., even when $\mathcal{G}$ does not occur, we instead solve
$\mathcal{P}_{N}(x_{t},\theta_{\rho(t)},\Theta_{\rho(t)},\bar{\mathcal{W}}_{t})$111More
precisely, the problem is obtained by replacing $\theta_{t}$ by
$\theta_{\rho(t)}$, $H_{t}^{(j)}$ by $H_{\rho(t)}^{(j)}$ for all $j$ in
$\mathcal{P}_{N}(x_{t},\theta_{t},\Theta_{t},\bar{\mathcal{W}}_{t})$. in
Algorithm 1 where,

$\displaystyle\rho(t):=\max\\{\tau\leq
t:\mathcal{P}_{N}(x_{t},\theta_{\tau},\Theta_{\tau},\bar{\mathcal{W}}_{t})\
\textnormal{is feasible}\\}.$

Essentially, in the unlikely event that
$\mathcal{P}_{N}(x_{t},\theta_{t},\Theta_{t},\bar{\mathcal{W}}_{t})$ is not
feasible, we instead solve
$\mathcal{P}_{N}(x_{t},\theta_{\tau},\Theta_{\tau},\bar{\mathcal{W}}_{t})$
using the latest estimates $\theta_{\tau}$ and $\Theta_{\tau}$ for which the
problem is feasible. This modification is of little practical consequence as
we can (and typically do) choose $\delta$ to be very small.

##  V Analysis

This section is devoted to the analysis of the Self-tuning Tube-based MPC
(STT-MPC) algorithm. We first provides guarantees for the LSE, and then
establish the recursive feasibility and the input-to-state stability of the
algorithm.

###  V-A Performance of the LSE

######  Theorem 1.

Assume that the optimization problem
$\mathcal{P}_{N}(x_{0},\theta_{0},\Theta_{0},\bar{\mathcal{W}}_{0})$ is
feasible for initial state $x_{0}\in{\mathbb{R}}^{d_{x}}$ and parameter
$\theta_{0}\in\Theta_{0}$. Under the STT-MPC algorithm, we have: for all
$\delta\in(0,1)$222Note that the result holds for all $\delta\in(0,1)$ and in
particular for the parameter $\delta$ chosen as input of the STT-MPC
algorithm.,

$\mathbb{P}[{\cal
G}]=\mathbb{P}[\|\hat{\theta}_{t}-\theta^{\star}\|\leq\epsilon_{t},\forall
t\geq t^{\star}(\delta)]\geq 1-\delta.$

###### Proof.

The proof follows a similar path as that used to prove Proposition 1 in [7],
and we provide a sketch of the arguments below. We first decompose the mean
square error as (see the proof of Proposition 1 in [7]):

$\displaystyle\|\hat{\theta}_{t}-\theta^{\star}\|^{2}\leq\left\|\left(\sum_{s=0}^{t-2}y_{s}y_{s}^{\top}\right)^{-1/2}\right.$
$\displaystyle\left.\left(\sum_{s=0}^{t-2}y_{s}w_{s}^{\top}\right)\right\|^{2}$
$\displaystyle\times\frac{1}{\lambda_{\min}(\sum_{s=0}^{t-2}y_{s}y_{s}^{\top})}.$

The first term in the r.h.s. of the above inequality is referred to as the
self-normalizing term, and it can be analyzed using Hanson-Wright inequality.
Specifically, we can show as in the proof of Proposition 1 in [7] that it is
smaller than $c_{4}\log(t/\delta)$ with probability at least $1-\delta/2$ for
some constant $c_{4}>0$ and for any $\delta\in(0,1)$. To prove this upper
bound, we leverage the fact that the tube MPC problem solved under STT-MPC is
feasible in each step (see Theorem 2 (ii)), which in turn ensures that $y_{t}$
is bounded, and simplifies the analysis of the self-normalizing term. Next, we
can derive a lower bound of the minimal eigenvalue of the cumulative covariate
matrix using the techniques developed in Appendix G in [7]. As in Proposition
11 in [7], we can show that for any $\delta\in(0,1)$,
$\lambda_{\min}(\sum_{s=0}^{t-2}y_{s}y_{s}^{\top})\geq c_{5}t^{1-\alpha}$ with
probability at least $1-\delta$, provided that $t^{1-\alpha}\geq
c_{6}\log(e/\delta)$, for some constants $c_{5},c_{6}>0$. We conclude that for
any $\delta\in(0,1)$, $\|\hat{\theta}_{t}-\theta^{\star}\|\leq\epsilon_{t}$
with probability at least $1-\delta$ when $t^{1-\alpha}\geq
c_{6}\log(e/\delta)$. The upper bound of $\|\hat{\theta}_{t}-\theta^{\star}\|$
can be made uniform over time $t\geq t^{\star}(\delta)$ applying Lemma 23 in
[7], which concludes the proof. $\Box$

###  V-B Recursive feasibility

######  Theorem 2.

If the optimization problem
$\mathcal{P}_{N}(x_{0},\theta_{0},\Theta_{0},\bar{\mathcal{W}}_{0})$ is
feasible for initial state $x_{0}\in{\mathbb{R}}^{d_{x}}$ and parameter
$\theta_{0}\in\Theta_{0}$, then, for all $t>0$,

  1. (i)

under event $\mathcal{G}$, the problem
$\mathcal{P}_{N}(x_{t},\theta_{t},\Theta_{t},\bar{\mathcal{W}}_{t})$ is
feasible;

  2. (ii)

the problem
$\mathcal{P}_{N}(x_{t},\theta_{\rho(t)},\Theta_{\rho(t)},\bar{\mathcal{W}}_{t})$
is feasible.

###### Proof.

For (i), we have that $\mathcal{G}$ holds which, together with Assumption 1,
ensures that $\theta^{*}\in\Theta_{t}$ for all $t>0$. Let
$\mathcal{P}_{N}(x_{t-1},\theta_{t-1},\Theta_{t-1},\bar{\mathcal{W}}_{t-1})$
be feasible with solution
$(\boldsymbol{v}_{t-1}^{*},\boldsymbol{\alpha}_{t-1}^{*})$. We define the tail
point $(\hat{\boldsymbol{v}}_{t},\hat{\boldsymbol{\alpha}}_{t})$ as

$\displaystyle\hat{\boldsymbol{v}}_{t}=\\{v^{*}_{1|t-1},\dots,v^{*}_{N-1|t-1},0\\},$
$\displaystyle\hat{\boldsymbol{\alpha}}_{t}=\\{\alpha_{1|t-1},\dots,\alpha_{N|t-1},\alpha_{N|t-1}\\},$

and we show that it is a feasible point for
$\mathcal{P}_{N}(x_{t},\theta_{t},\Theta_{t},\bar{\mathcal{W}}_{t})$, i.e.,
that it satisfies the constraints (14)-(19) at time $t$.

The initial constraint (14) is satisfied since
$(\boldsymbol{v}_{t-1}^{*},\boldsymbol{\alpha}_{t-1}^{*})$ satisfies (16), and
$\theta^{*}\in\Theta_{t-1}$, which together imply that
$x_{t}\in\mathcal{X}_{1|t-1}$. The system dynamics constraint (15) holds
automatically. We now consider only the vertices of $\Theta_{t-1}$ and note
that $\bar{w}_{t-1}\geq\bar{w}_{t}$ and
$\bar{\zeta}_{t-1}\geq\bar{\zeta}_{t}$. At time $t$, the tube constraints (16)
and (17) for $k\in\\{0,\dots,N-2\\}$ are satisfied whenever the tube
constraints for $k\in\\{1,\dots,N-1\\}$ at time $t-1$ are satisfied.
Furthermore, noting that $v_{N-1|t}=0$ and $\alpha_{N|t}=\alpha^{*}_{N|t-1}$,
by simple substitution, we will have that (16) and (17) for $k=N-1$ at time
$t$ hold when the terminal constraints (18) and (19) at time $t-1$ hold. If
the terminal constraints hold at time $t-1$, they will also hold at time $t$.
Finally, since $\Theta_{t}\subseteq\Theta_{t-1}$, these constraints will also
hold for every vertex of $\Theta_{t}$.

For (ii) note that, due to Assumption 1, $\theta^{*}\in\Theta_{0}$ so that
there will always exist $\tau\leq t$ such that $\theta^{*}\in\Theta_{\tau}$
and such that
$\mathcal{P}_{N}(x_{\tau},\theta_{\tau},\Theta_{\tau},\bar{\mathcal{W}}_{\tau})$
is feasible. Then, using a similar argument as in (i) we will have that
$\mathcal{P}_{N}(x_{\tau+1},\theta_{\tau},\Theta_{\tau},\bar{\mathcal{W}}_{\tau+1})$
is feasible and thus
$\mathcal{P}_{N}(x_{t+1},\theta_{\tau},\Theta_{\tau},\bar{\mathcal{W}}_{t+1})$
is feasible. $\Box$

###  V-C Input-to-State Stability (ISS)

Let $X_{\mathcal{P}}$ be the set of states $x_{0}$ such that problem
$\mathcal{P}_{N}(x_{0},\theta_{0},\Theta_{0},\bar{\mathcal{W}}_{0})$ is
feasible. Additionally, let $J(x_{t},\boldsymbol{v}_{t},\theta_{t})$ be the
objective of problem
$\mathcal{P}_{N}(x_{t},\theta_{t},\Theta_{t},\bar{\mathcal{W}}_{t})$, let
$\boldsymbol{v}_{t}^{*}$ be its optimal solution and let
$V^{*}(x_{t},\theta_{t})$ be its associated value. To simplify notations, we
refer to the closed loop system (6) as,

$\displaystyle x_{t+1}=f(x_{t},v_{t},w_{t},\theta_{t})$ (21)

where, we refer to $B(\theta_{t})\zeta_{t}+w_{t}$ as $w_{t}$. We further
define the stage cost $L(x,v)=x^{\top}(Q+K^{\top}RK)x+v^{\top}Rv$, and
terminal cost $\phi(x,\theta)=x^{\top}P(\theta)x$ so that

$\displaystyle
J(x_{t},\boldsymbol{v}_{t},\theta_{t})=\sum_{k=0}^{N-1}L(x_{k|t},v_{k|t})+\phi(x_{N|t},\theta_{t}).$

######  Lemma 2 (Theorem 3 in [23]).

The system (21) with control law $v_{t}(x_{t},\theta_{t})$ is ISS with region
of attraction $\mathcal{R}\subseteq{\mathbb{R}}^{d_{x}}$ if

  1. (i)

$\mathcal{R}$ is a compact set containing the origin, and is robust positive
invariant for (21).

  2. (ii)

There exist $\mathcal{K}_{\infty}$-functions
$\alpha_{1},\alpha_{2},\alpha_{3}$, $\mathcal{K}$-functions
$\sigma_{1},\sigma_{2}$, and a continuous Lyapunov function
$V:\mathcal{R}\to{\mathbb{R}}_{+}$ such that for all states
$x_{t}\in\mathcal{R}$ and disturbances $w_{t}\in\bar{\mathcal{W}}_{t}$ we
have,

$\displaystyle\alpha_{1}(\lVert x_{t}\rVert)\leq V(x_{t})$
$\displaystyle\leq\alpha_{2}(\lVert x_{t}\rVert)$ (22) $\displaystyle
V(x_{t+1})-V(x_{t})$ $\displaystyle\leq-\alpha_{3}(\lVert
x_{t}\rVert)+\sigma_{1}(\lVert w_{t}\rVert)$
$\displaystyle\quad+\sigma_{2}(\lVert\theta_{t}-\theta^{*}\rVert).$ (23)

######  Theorem 3.

Under event $\mathcal{G}$ and under the STT-MPC algorithm, for all initial
conditions $x_{0}\in X_{\mathcal{P}}$ and all parameters
$\theta^{*}\in\Theta_{0}$, the system (1) is ISS with region of attraction
$X_{\mathcal{P}}$.

It suffices to show that $X_{\mathcal{P}}$ satisfies condition (i) of Lemma 2
and $V^{*}(x_{t},\theta_{t})$ is a Lyapunov function satisfying condition
(ii).

For condition (i), suppose that
$\mathcal{P}_{N}(x_{0},\theta_{0},\Theta_{0},\bar{\mathcal{W}}_{0})$ is
feasible, i.e. $X_{\mathcal{P}}$ is not empty. Then the terminal constraints
(18) and (19) ensure that there exists an $\alpha_{N|0}$ such that the
terminal cross-section set $X_{N|0}=\\{x:Tx\leq\alpha_{N|0}\\}$ satisfies the
constraints (2) and $f(x,0,\theta,w)\in X_{N|0}$ for all $x\in X_{N|0}$,
$w\in\bar{\mathcal{W}}_{0}$, and $\theta\in\Theta_{0}$, i.e. it is robustly
invariant. For every $x\in X_{N|0}$, the point
$(\boldsymbol{v},\boldsymbol{\alpha})=(\\{0,\dots,0\\},\\{\alpha_{N|0},\dots,\alpha_{N|0}\\})$
is feasible for
$\mathcal{P}_{N}(x,\theta_{0},\Theta_{0},\bar{\mathcal{W}}_{0})$, thus
$X_{N|0}\subseteq X_{\mathcal{P}}$. The robust positive invariance of
$X_{\mathcal{P}}$ will then follow directly from Theorem 2. Since
$X_{N|t}\subseteq X_{\mathcal{P}}$ and $0\in\bar{\mathcal{W}}_{t}$ by
Assumption 2, we have $0\in X_{\mathcal{P}}$. Finally, $X_{\mathcal{P}}$ is
compact due to Assumption 3.

For condition (ii) we first show that the bound (22) holds for
$V^{\star}(x_{t},\theta_{t})$. For any given state $x_{t}$, and parameters
$\theta_{t}\in\Theta_{t}$, and for any $Q,R\succ 0$, problem
$\mathcal{P}_{N}(x_{t},\theta_{t},\Theta_{t},\bar{\mathcal{W}}_{t})$ will be a
quadratic program and it follows from Theorem 4 in [24] that
$V^{\star}(x_{t},\theta_{t})$ is a continuous function of $x_{t}$. Since
$\Theta_{t}$ is compact by Assumption 1, it follows from Lemma 4.2 in [25]
that there exist $\mathcal{K}_{\infty}$-functions $\alpha_{1}$, $\alpha_{2}$
such that for any $x_{t}\in X_{\mathcal{P}}$,

$\displaystyle\alpha_{1}(\lVert x_{t}\rVert)\leq
V^{\star}(x_{t},\theta_{t})\leq\alpha_{2}(\lVert x_{t}\rVert).$

We now show the bound (23). First, note that the continuity of
$f(x,u,w,\theta)$, $L(x,u)$, and $\phi(x,\theta)$ with respect to $x,u,w$, and
$\theta$, implies that there exist $\mathcal{K}_{\infty}$-functions
$\sigma_{x},\sigma_{w},\sigma_{\theta},\sigma_{L},\sigma_{\phi}$, so that for
all $x,x_{1},x_{2}\in X_{\mathcal{P}}$, $v\in{\mathbb{R}}^{d_{u}}$,
$w,w_{1},w_{2}\in\bar{\mathcal{W}}_{0}$, and
$\theta,\theta_{1},\theta_{2}\in\Theta_{t}$,

$\displaystyle\lVert
f(x_{1},v,w_{1},\theta_{1})-f(x_{2},v,w_{2},\theta_{2})\rVert$
$\displaystyle\leq\sigma_{x}(\lVert x_{1}-x_{2}\rVert)$
$\displaystyle+\sigma_{w}(\lVert w_{1}-w_{2}\rVert)$
$\displaystyle+\sigma_{\theta}(\lVert\theta_{1}-\theta_{2}\rVert),$
$\displaystyle|L(x_{1},v)-L(x_{2},v)|$ $\displaystyle\leq\sigma_{L}(\lVert
x_{1}-x_{2}\rVert),$
$\displaystyle\lVert\phi(x_{1},\theta)-\phi(x_{2},\theta)\rVert$
$\displaystyle\leq\sigma_{\phi}(\lVert x_{1}-x_{2}\rVert).$

Recall that $\boldsymbol{v}_{t}^{*}=\\{v_{0|t}^{*},\dots,v_{N-1|t}^{*}\\}$ is
the predicted control sequence that arises from the solution of
$\mathcal{P}_{N}(x_{t},\theta_{t},\Theta_{t},\bar{\mathcal{W}}_{t})$. We
denote by $\boldsymbol{x}_{t}^{*}=\\{x_{0|t}^{*},\dots,x_{N|t}^{*}\\}$, the
associated predicted state trajectory (using (4), with parameter $\theta_{t}$,
when applying $\boldsymbol{v}_{t}^{*}$). Similarly, recall the tail control
sequence
$\hat{\boldsymbol{v}}_{t+1}:=\\{\hat{v}_{0|t+1},\dots,\hat{v}_{N-1|t+1}\\}=\\{v_{1|t}^{*},\dots,v_{N-1|t}^{*},0\\}$
for which we define, an associated predicted state trajectory
$\hat{\boldsymbol{x}}_{t+1}=\\{\hat{x}_{0|t+1},\dots,\hat{x}_{N|t+1}\\}$
(using (4) with parameter $\theta_{t+1}$, when applying
$\hat{\boldsymbol{v}}_{t+1}$) with the property that
$\hat{x}_{0|t+1}=x_{1|t}^{*}$. To simplify the presentation, we omit
references to the subscript $t$ when there is no ambiguity. Noting the sub-
optimality of $\hat{\boldsymbol{v}}_{t+1}$, we will have

$\displaystyle J($ $\displaystyle
x_{t+1},\hat{\boldsymbol{v}}_{t+1},\theta_{t+1})-J(x_{t},\boldsymbol{v}^{*}_{t},\theta_{t})$
$\displaystyle\leq V^{*}(x_{t+1},\theta_{t+1})-V^{*}(x_{t},\theta_{t})$
$\displaystyle\leq\sum_{k=0}^{N-1}L(\hat{x}_{k},\hat{v}_{k})-L(x_{k}^{*},v_{k}^{*})+\phi(\hat{x}_{N},\theta_{t+1})-\phi(x_{N}^{*},\theta_{t})$
$\displaystyle=-L(x_{0|t},v^{*}_{0})+\sum_{k=0}^{N-2}L(\hat{x}_{k},v_{k+1}^{*})-L(x_{k+1}^{*},v_{k+1}^{*})$
$\displaystyle\quad+L(\hat{x}_{N-1},0)+\phi(\hat{x}_{N},\theta_{t+1})-\phi(\hat{x}_{N-1},\theta_{t+1})$
$\displaystyle\quad+\phi(\hat{x}_{N-1},\theta_{t+1})-\phi(x_{N}^{*},\theta_{t}).$

We will proceed to upper bound the above quantities. We first note that, for
some $\mathcal{K}_{\infty}$ function $\kappa_{\epsilon}$ we will have,

$\displaystyle\sigma_{\theta}(\lVert\theta_{t+1}-\theta_{t}\rVert)\leq\sigma_{\theta}(\lVert\theta_{t}-\theta^{*}\rVert)+\kappa_{\epsilon}(\epsilon_{t})$

and for all $k\in\\{1,\dots,N-1\\}$,

$\displaystyle\lVert\hat{x}_{k}-x^{*}_{k+1}\rVert$
$\displaystyle\leq\sigma_{x}(\lVert\hat{x}_{k-1}-x^{*}_{k}\rVert)+\sigma_{\theta}(\lVert\theta_{t+1}-\theta_{t}\rVert)$
$\displaystyle\leq\sigma_{x}(\lVert\hat{x}_{k-1}-x^{*}_{k}\rVert)+\sigma_{\theta}(\lVert\theta_{t}-\theta^{*}\rVert)+\kappa_{\epsilon}(\epsilon_{t}).$

Here we the triangle inequality, the fact that event $\mathcal{G}$ holds, and
the monotonicity of $\epsilon_{t}$. Applying this relation recursively leads
to,

$\displaystyle\lVert\hat{x}_{k}-x^{*}_{k+1}\rVert\leq$
$\displaystyle\frac{1}{2}(2\sigma_{x})^{k}(\lVert\hat{x}_{0}-x^{*}_{1}\rVert)$
$\displaystyle+\sum_{i=0}^{k-1}\frac{1}{2}(2\sigma_{x})^{i}\circ
2(\sigma_{\theta}(\lVert\theta_{t}-\theta^{*}\rVert+\kappa_{\epsilon}(\epsilon_{t})),$

where
$\lVert\hat{x}_{0}-x^{*}_{1}\rVert\leq\sigma_{w}(w_{t})+\sigma_{\theta}(\lVert\theta_{t}-\theta^{*}\rVert)+\kappa_{\epsilon}(\epsilon_{t})$.

Now, since $Q,R\geq 0$, there exists a $K_{\infty}$-function $\alpha_{3}$ such
that $-L(x_{0|t},v^{*}_{0})\leq-\alpha_{3}(\lVert x_{0|t}\rVert)$. We will
also have, for all $k\in\\{0,\dots,N-2\\}$,

$\displaystyle|L(\hat{x}_{k},v_{k+1}^{*})-L(x_{k+1}^{*},v_{k+1}^{*})|\leq\sigma_{L}(\lVert\hat{x}_{k}-x^{*}_{k+1}\rVert).$

Then, from (20) it follows that,

$\displaystyle
L(\hat{x}_{N-1},0)+\phi(\hat{x}_{N,\theta_{t+1}})-\phi(\hat{x}_{N-1},\theta_{t+1})=0.$

Finally, recall that (20) is linear in $P(\theta)$ and note that its solution
is unique for all $\theta\in\Theta_{0}$ (Theorem 3.7 [25] using that
$\Phi(\theta)$ is stable). It follows then that $P(\theta)$ is Lipschitz
continuous and that there will exist a $\mathcal{K}_{\infty}$-function
$\kappa_{\phi}$ such that,

$\displaystyle|\phi(\hat{x}_{N-1},\theta_{t+1})-\phi(x^{*}_{N},\theta_{t})|$
$\displaystyle\leq\sigma_{\phi}(\lVert\hat{x}_{N-1}-x^{*}_{N}\rVert)$
$\displaystyle+\kappa_{\phi}(\lVert\theta_{t}-\theta^{*}\rVert)+\kappa_{\epsilon}(\epsilon_{t}).$

We now group the above terms together to get,

$\displaystyle V^{*}(x_{t+1},\theta_{t+1})-V^{*}(x_{t},\theta_{t})\leq$
$\displaystyle-\alpha_{3}(\lVert x_{0|t}\rVert)+\bar{\sigma}_{w}(\lVert
w_{t}\rVert)$
$\displaystyle+\bar{\sigma}_{\theta}(\lVert\theta_{t}-\theta^{*}\rVert)+\bar{\kappa}_{\epsilon}(e_{t})$

where $\bar{\sigma}_{w},\bar{\sigma}_{\theta}$, and $\bar{\kappa}_{\epsilon}$
are $\mathcal{K}$-functions. Since $e_{t}\to 0$ as $t\to\infty$, we will have
that condition (iii) is satisfied. $\Box$

##  VI Numerical Experiments

Figure 1: Average volume of $\textrm{Vol}(\Theta_{t})$ for different values of
$\delta$. Results are averaged over 10 different runs (the small shadowed
areas depicts 99.99% confidence interval).

We illustrate the performance of STT-MPC on a second-order linear systems with
true parameters

$A=\begin{bmatrix}0.6&0.2\\\ -0.1&0.4\end{bmatrix},\quad B=\begin{bmatrix}1\\\
0.6\end{bmatrix},\quad\sigma=0.01,$ (24)

where $B$ is assumed to be known, while $A$ is the uncertain parameter. The
initial state is $x_{0}=(6,3)$, and we consider $\Theta_{0}$ to be a
$4$-dimensional hypercube centered at
$\theta_{0}=\begin{bmatrix}0.57&0.17&-0.12&0.42\end{bmatrix}$ with side length
$0.14$. Consequently, the resulting stabilizing feedback gain is
$K=\begin{bmatrix}-0.426&-0.290\end{bmatrix}$. The least-squares confidence
region $\Delta_{t}$ is over-approximated by an hypercube centered in
$\theta_{t}$ with side length $4\varepsilon_{t}$. We consider the following
state and input constraints: $[x_{t}]_{1}\geq-0.15,[x_{t}]_{2}\geq-1.1$ and
$u_{t}\leq 0.5$. The matrix $T$ is computed according to the relation (5.98)
in [21], with $\lambda=0.999$. We consider $Q=I_{d_{x}}$, $R=I_{d_{u}}$, and,
for simplicity, we let the worst case noise realization be
$\bar{w}_{t}=\bar{w}_{0}$, for all $t\geq 1$. All the simulations were
performed in Python 3.9333Find the code here:
https://github.com/rssalessio/SelfTuningTubeBasedMPC, using the CVXPY library
[26] and the MOSEK solver.

TABLE I: Average volume of $\textrm{Vol}(\Theta_{t})$ with respect to $\textrm{Vol}(\Theta_{0})$ Time step $t$ | 1 | 5 | 15 | 50 | 100
---|---|---|---|---|---
$\delta=10^{-1}$ | $100\%$ | $1.312\%$ | $0.401\%$ | $0.071\%$ | $0.023\%$
$\delta=10^{-2}$ | $100\%$ | $3.354\%$ | $0.827\%$ | $0.125\%$ | $0.038\%$
$\delta=10^{-3}$ | $100\%$ | $7.058\%$ | $1.525\%$ | $0.205\%$ | $0.059\%$
$\delta=10^{-4}$ | $100\%$ | $12.651\%$ | $2.593\%$ | $0.319\%$ | $0.088\%$

Figure 1 shows the average volume of $\Theta_{t}$ over 10 different runs. We
note that the reduction in volume is significant, and well below $1\%$ of the
volume of $\Theta_{0}$ after $50$ steps (see also the details in table I;
after 100 steps the volume is reduced below $0.1\%$ of the original volume).
In Figure 2 we present the trajectory of a single run of STT-MPC. In the plot
are shown the tubes computed by the algorithm. We note that the constraints
are robustly satisfied. Moreover, we also note a sensible reduction of the
tubes over time.

Figure 2: Illustration of the state trajectory using STT-MPC, and its tubes,
over a single run. The gray areas depict the state constraints
$[x_{t}]_{1}\geq-0.15,[x_{t}]_{2}\geq-1.1$. As a comparison, the plot of the
state trajectory under $u_{t}=Kx_{t}$ is also shown.

##  VII Conclusions

We proposed STT-MPC which combines least-squares estimaton with a polytopic
tube-based MPC method to ensure robust constraint satisfaction while learning
the system dynamics. In contrast with existing conservative set-based
identification approaches, we construct parameter sets that hold for all time
with a confidence level $\delta$ chosen as a design parameter. Persistent
excitation is ensured by injecting a truncated noise signal which decays at a
rate $t^{-\alpha}$, with $\alpha>0$ controlling the trade-off between
(transient) estimation accuracy and controller performance. Importantly, we
asymptotically recover the performance of the equivalent tube-based MPC that
has full knowledge of the dynamics.

We analysed our proposed algorithm by deriving performance guarantees for the
LSE, as well as by establishing recursive feasibility and input-to-state
stability guarantees. We demonstrated the performance of our adaptation
scheme, its dependence on $\delta$, and the robust constraint satisfaction of
STT-MPC via a numerical example.

Future work will focus on the analysis of the performance of our algorithm and
quantifying the additional objective cost that must be paid for learning the
dynamics. In particular, we will investigate whether a PE noise decaying as
$1/\sqrt{t}$ could ensure a regret scaling at most as $\sqrt{t}$ (as this was
established in [7] for LQR problems).

## References

  * [1] D. Q. Mayne, J. B. Rawlings, C. V. Rao, and P. O. Scokaert, “Constrained model predictive control: Stability and optimality,” _Automatica_ , vol. 36, no. 6, pp. 789–814, 2000.
  * [2] D. Q. Mayne, M. M. Seron, and S. Raković, “Robust model predictive control of constrained linear systems with bounded disturbances,” _Automatica_ , vol. 41, no. 2, pp. 219–224, 2005.
  * [3] D. Q. Mayne, “Model predictive control: Recent developments and future promise,” _Automatica_ , vol. 50, no. 12, pp. 2967–2986, 2014.
  * [4] K. J. Åström and B. Wittenmark, “On self tuning regulators,” _Automatica_ , vol. 9, no. 2, pp. 185–199, 1973.
  * [5] Y. Abbasi-Yadkori and C. Szepesvári, “Regret bounds for the adaptive control of linear quadratic systems,” in _Proceedings of the 24th Annual Conference on Learning Theory_. JMLR Workshop and Conference Proceedings, 2011, pp. 1–26.
  * [6] H. Mania, S. Tu, and B. Recht, “Certainty equivalence is efficient for linear quadratic control,” _Advances in Neural Information Processing Systems_ , vol. 32, 2019.
  * [7] Y. Jedra and A. Proutiere, “Minimal expected regret in linear quadratic control,” in _International Conference on Artificial Intelligence and Statistics_. PMLR, 2022, pp. 10 234–10 321.
  * [8] J. Fleming, B. Kouvaritakis, and M. Cannon, “Robust tube mpc for linear systems with multiplicative uncertainty,” _IEEE Transactions on Automatic Control_ , vol. 60, no. 4, pp. 1087–1092, 2014.
  * [9] K. S. Narendra and A. M. Annaswamy, “Persistent excitation in adaptive systems,” _International Journal of Control_ , vol. 45, no. 1, pp. 127–160, 1987.
  * [10] D. Tranos, “Regret minimization in structured reinforcement learning,” Ph.D. dissertation, KTH Royal Institute of Technology, 2021.
  * [11] G. Marafioti, R. R. Bitmead, and M. Hovd, “Persistently exciting model predictive control,” _International Journal of Adaptive Control and Signal Processing_ , vol. 28, no. 6, pp. 536–552, 2014.
  * [12] T. A. N. Heirung, B. Foss, and B. E. Ydstie, “Mpc-based dual control with online experiment design,” _Journal of Process Control_ , vol. 32, pp. 64–76, 2015.
  * [13] T. A. N. Heirung, B. E. Ydstie, and B. Foss, “Dual adaptive model predictive control,” _Automatica_ , vol. 80, pp. 340–348, 2017.
  * [14] Z. Cai and B. E. Ydstie, “Dual adaptive model predictive control with disturbances,” _IFAC-PapersOnLine_ , vol. 54, no. 3, pp. 206–211, 2021.
  * [15] A. Aswani, H. Gonzalez, S. S. Sastry, and C. Tomlin, “Provably safe and robust learning-based model predictive control,” _Automatica_ , vol. 49, no. 5, pp. 1216–1226, 2013.
  * [16] S. Di Cairano, “Indirect adaptive model predictive control for linear systems with polytopic uncertainty,” in _2016 American Control Conference (ACC)_. IEEE, 2016, pp. 3570–3575.
  * [17] M. Tanaskovic, L. Fagiano, R. Smith, and M. Morari, “Adaptive receding horizon control for constrained mimo systems,” _Automatica_ , vol. 50, no. 12, pp. 3019–3029, 2014.
  * [18] M. Lorenzen, M. Cannon, and F. Allgöwer, “Robust mpc with recursive model update,” _Automatica_ , vol. 103, pp. 461–471, 2019.
  * [19] S. V. Raković, B. Kouvaritakis, R. Findeisen, and M. Cannon, “Homothetic tube model predictive control,” _Automatica_ , vol. 48, no. 8, pp. 1631–1638, 2012.
  * [20] X. Lu and M. Cannon, “Robust adaptive tube model predictive control,” in _2019 American Control Conference (ACC)_. IEEE, 2019, pp. 3695–3701.
  * [21] B. Kouvaritakis and M. Cannon, “Model predictive control,” _Switzerland: Springer International Publishing_ , vol. 38, 2016.
  * [22] F. Blanchini and S. Miani, _Set-theoretic methods in control_. Springer, 2008, vol. 78.
  * [23] D. Limon, T. Alamo, D. M. Raimondo, D. Pena, J. M. Bravo, A. Ferramosca, and E. F. Camacho, “Input-to-state stability: a unifying framework for robust model predictive control,” in _Nonlinear model predictive control_. Springer, 2009, pp. 1–26.
  * [24] A. Bemporad, M. Morari, V. Dua, and E. N. Pistikopoulos, “The explicit linear quadratic regulator for constrained systems,” _Automatica_ , vol. 38, no. 1, pp. 3–20, 2002.
  * [25] H. K. Khalil, _Nonlinear control_. Pearson New York, 2015, vol. 406.
  * [26] S. Diamond and S. Boyd, “Cvxpy: A python-embedded modeling language for convex optimization,” _The Journal of Machine Learning Research_ , vol. 17, no. 1, pp. 2909–2913, 2016.