11institutetext: Dipartimento di Fisica, Università di Ferrara, via Saragat 1,
44122 Ferrara, Italy
11email<EMAIL_ADDRESS>22institutetext: Keldysh Institute of Applied
Mathematics, Russian Academy of Sciences. 4 Miusskaya sq., Moscow, 125047,
Russia
22email<EMAIL_ADDRESS>
# On a discontinuity at the base of the transition layer located between the
Keplerian accretion disk and the compact object

L. Titarchuk 11 I. Kalashnikov 22

(Received; accepted.)

###### Abstract

Context. We study the geometry of the transition layer (TL) between the
classical Keplerian accretion disk (the TL outer boundary) and the compact
object at the TL inner boundary.

Aims.  Our goal is to demonstrate using the hydrodynamical formalism that the
TL is created along with a shock due to a discontinuity and to an adjustment
of the Keplerian disk motion to a central object.

Methods. We apply hydrodynamical equations to describe a plasma motion near a
central object in the TL.

Results.  We point out that before matter accretes to a central object the TL
cloud is formed between an adjustment radius and the TL inner boundary which
is probably a site where the emergent Compton spectrum comes from. Using a
generalization of the Randkine-Hugoniot relation and a solution of the
azimutal force balance equation we have reproduced the geometric
characteristics of TL.

###### Key Words.:

accretion, accretion disks–radiation mechanisms: general– black hole
physics–neutron stars–white dwarfs

##  1 Introduction

Accretion onto compact objects as white dwarf (WD), neutron star (NS), and
black hole (BH) binaries are very similar. The accreting matter with
relatively large angular momentum can form a disk or go to outflow (see
Titarchuk et al. (2007) and Shaposhnikov & Titarchuk (2009), hereafter TSA07,
ST09, respectively). On the other hand, the plasma with low angular momentum
proceeds towards the compact object almost in a free-fall manner until the
centrifugal force becomes strong to halt the flow (see, e.g., (Chakrabarti &
Titarchuk 1995), hereafter CT95) Thus, we can suggest that three definitive
regions a (Keplerian) disk, the shock, and a transition layer (TL) between the
last stable orbit near WD, NS or BH and the shock, which are probably a place
where the emergent X-ray spectrum is formed. In this paper, we study the
characteristics of this adjustment of the Keplerian disk flow to an WD, an NS
or an BH.

The disk starts to deviate from a Keplerian motion at a certain radius to
adapt itself to the boundary conditions at WD or NS surface (or at the last
stable orbit, at $r_{1}=3r_{\rm S}$, where $r_{\rm S}$ is the Schwarzchild
radius, in the case of a BH). Titarchuk et al. (1998) (hereafter TLM98)
explained the millisecond variability detected by Rossi X-Ray Timing Explorer
(RXTE) in the X-ray emission from a number of low-mass X-ray binary systems
(XRBs). Later, Seifina & Titarchuk (2010), Farinelli & Titarchuk (2011), and
so on analyzed X-ray data from Sco X-1, 4U 1728$-$34, 4U 1608$-$522, 4U
1636$-$536, 4U 0614$-$091, 4U 1735$-$44, 4U 1820$-$30, GX 5$-$1) in terms of
dynamics of the centrifugal barrier, in a hot boundary region (the TL)
surrounding a neutron star (NS). They demonstrated that this region may
experience relaxation oscillations and that the displacements of a gas element
both in radial and vertical directions occur at the same main frequency, of an
order of the local Keplerian frequency.

Observations of black hole X-ray binaries and active galactic nuclei indicate
that the accretion flows around black holes are composed of hot and cold gas,
which have been theoretically described in terms of a hot corona next to an
optically thick relatively cold disk. Liu & Qiao (2022) (hereafter LQ22)
reviews the accretion flows around black holes, with an emphasis on the
physics that determines the configuration of hot and cold accreting gas, and
how the configuration varies with the accretion rate and thereby produces
various luminosity and spectra. They provide references to the famous
solutions of the standard disk model such as, for example, Shakura & Sunyaev
(1973) (hereafter SS73). The most successful applications of this model are
the steady state and time-dependent nature of thin disks in dwarf novae and
soft X-ray transients. Application of the four accretion models to black holes
presented by LQ22 is constrained by both the theoretical assumption of
specific models and observational luminosity and spectrum.

The soft photons from the innermost part of the disk are scattered off the hot
electrons thus forming the Comptonized X-ray spectrum [see Sunyaev & Titarchuk
(1980), hereafter ST80]. The electron temperature is regulated by the supply
of soft photons from a disk, which depends on the ratio of the energy release
(accretion rate) in a disk and the energy release in the TL. For example, the
electron temperature is higher for lower accretion rates (see, e.g., CT95),
while for a high accretion rate (of an order of the Eddington one) the TL is
cooled down very efficiently due to the soft photon flux and the
Comptonization. This is a possible mechanism for the low-frequency QPOs
(quasi-periodic oscillations), discovered by RXTE in a number of low-mass
X-ray binary systems (LMXBs; see Strohmayer et al. (1996); van der Klis et al.
(1996), hereafter S96 and VK96a, respectively; Zhang et al. (1996), TLM98,
Titarchuk et al. (2007) hereafter TSA07). These observations reveal a wealth
of high and low-frequency X-ray variabilities that are believed to be due to
the processes occurring in the very vicinity of an accreting NS and BH.

Lančová et al. (2019) discovered a new class of black hole accretion disk
solutions through 3D radiative magnetohydrodynamic simulations applying the
general relativity formalism. These solutions combine features of thin, slim,
and thick disk models. They provide a realistic description of black hole
disks and have a thermally stable, optically thick, Keplerian region supported
by magnetic fields. Mishra et al. (2020) have conducted a detailed analysis of
two-dimensional viscous, radiation hydrodynamic numerical simulations of the
Shakura-Sunyaev thin disks around a stellar mass black hole. They found
multiple robust, coherent oscillations occurring in the disk, including a
trapped fundamental g-mode and standing-wave p-modes. The study suggests that
these findings could be of astrophysical importance in observed twin peak,
high-frequency QPOs. Large scale 3D magnetohydrodynamic simulations of
accretion onto magnetized stars with tilted magnetic and rotational axes were
done by Romanova et al. (2021).It was shown that the inner parts of the disc
become warped, tilted, and precess due to magnetic interaction between the
magnetosphere and the disc. According to results of the numerical simulations
by Lugovskii & Chechetkin (2012), large-scale instability in accretion disks,
reveals changes in a disk flow structure because of t the formation of large
vortexes. Over time, this leads to the development of asymmetric spiral
structures and results in angular-momentum redistribution within the disk.

Suková et al. (2017) have conducted simulations of accretion flows with a low
angular momentum, filling the gap between spherically symmetric Bondi
accretion and disc-like accretion flows. They identify ranges of parameters
for which the shock after formation moves towards or outwards the central
black hole or the long-lasting oscillating shock is observed. The results are
scalable with a central black hole mass and can be compared to QPOs of
selected microquasars and supermassive black holes in the centres of weakly
active galaxies. In order to explain NICER telescope evidences for non-dipolar
magnetic field structures in rotation-powered millisecond pulsars, Das et al.
(2022) have conducted a suite of general relativistic magnetohydrodynamic
simulations of accreting neutron stars for dipole, quadrupole, and
quadrudipolar stellar field geometries. The study found that the location and
a size of the accretion columns resulting in hotspots changes significantly
depending on initial stellar field strength and geometry, providing a viable
mechanism to explain the radio power in observed neutron star jets.

Despite the increasing role of numerical simulations of accretion, the
consideration of relatively simple analytical models allows us to identify the
response of the system to changes in a set of parameters, which in the case of
numerical simulations becomes a very cumbersome task. In the present paper we,
along with other authors (e.g. Abramowicz & Fragile (2013); Ajay et al.
(2022))), propose a model of accretion on the central object. The point is
that existing models of shock waves (see e.g. Chakrabarti (1989)) cannot
correctly describe the abrupt transition between the accretion disk and the
TL. Our goal is to construct a satisfactory model of the shock wave lying at
the base of TL, which can correctly describe observational data for
temperature and density (Seifina & Titarchuk (2010)). In Sect. 2 we describe
our approach to the TL model. In Sect. 3. we formulate equations to determine
the vertical dependencies of density and pressure in the TL. In Secs. 4-5 we
introduce the generalized Rankine-Hugoniot relation and take into account the
rotation effect in the TL, respectively. We discuss our setup and final
results in Sect. 6. Finally, we summarize our conclusions in Sect. 7.

##  2 Motivation and the model

Figure 1: The schematic representation of the considered model.

Consideration of the two-dimensional structure of accretion disks is
associated with obvious mathematical difficulties. Typically, a two-
dimensional structure is only accounted for by an non-homogeneous distribution
of matter within a layer of constant thickness $2h$. Consideration of a disk
of constant thickness is made possible by ignoring the vertical velocity
$v_{z}$. The more general case $h=h(r)$ requires solving equation
$dh/dr=(v_{z}/u)|_{z=h(r)}$, where $u$ is the radial velocity, to describe the
behavior of the height as the radius changes.

According to the observations and their interpretations (TLM98, ST09) the disk
height may achieve very significant values near the TL, while at larger
distances it is pretty small. The thermodynamic characteristics of the plasma
change considerably: in the TL, its temperature can reach tens of keV and its
density drops drastically (see e.g. ST09). This change is widely believed to
be presumably due to the presence of a shock wave (see e.g. CT95, TLM98 and
ST09). Such a region near the central object (CO) is associated with the TL.

In an attempt to get closer to a two-dimensional consideration, we have chosen
the following model describing the shock (Fig. 1). The flow velocity is
supposed to have only radial and azimuthal components in the both sides of the
shock, which are independent of the height. It is assumed that the height of
the flow changes sharply by discontinuity from $h_{2}$ (the disk thickness) to
$h_{1}$ (the TL thickness).

We require the angular velocity $\omega$ of the TL to match to the CO at its
surface and coincide with disk velocity at the shock position $r_{2}$. This
angular velocity of the disk is supposed to equal to Keplerian one. Only
$r\phi$ component of the viscous stress tensor was considered because it is
responsible for taking off the angular momentum. We considered the stress
tensor form $\eta r\omega\omega^{\prime}$, where $\eta$ is the turbulent
viscosity.

The thermodynamic quantities of the flow have vertical dependencies which are
determined by the CO gravity and an equation of state. In order to identify
them, we did not assume the smallness of the vertical coordinate compared to
the radial one.

##  3 The vertical dependences

Let us start with the hydrostatic balance equation in the vertical direction:

$\displaystyle\frac{\partial p}{\partial z}=-\rho\frac{\partial\Phi}{\partial
z},$ (1)

where $p$, $\rho$ are pressure and density and $\Phi$ is the gravitational
potential creating by the central object. We supposed that the vertical
distribution of matter is given by polytropic relation $p=K\rho^{1+1/n}$,
where $K$ is a constant related to the entropy and $n$ is the polytropic
index. Then the formal solution of (1) may be written as:

$\displaystyle K(n+1)\rho^{1/n}=\widetilde{f}(r)-\Phi(r,z),$ (2)

where $\widetilde{f}$ is an arbitrary function. We require the density and
pressure to be zero at the inner boundary $z=\pm h$ and equal to some
functions $\rho_{e}$, $p_{e}$ correspondingly at the equator $z=0$. Then we
have:

$\displaystyle\rho=\rho_{e}(r)\left(\frac{\Phi(r,h)-\Phi(r,z)}{\Phi(r,h)-\Phi(r,0)}\right)^{n},$
(3) $\displaystyle
p=p_{e}(r)\left(\frac{\Phi(r,h)-\Phi(r,z)}{\Phi(r,h)-\Phi(r,0)}\right)^{n+1}.$
(4)

Assuming the height is small compared to the radius $|z|\ll r$ we may get the
usually used (Hōshi (1977); Matsumoto et al. (1984)) dependence
$\rho=\rho_{e}(r)(1-z^{2}/h^{2})^{n}$. However, we did not make such an
assumption later on.

We will need integrals from such vertical dependencies:

$\displaystyle\zeta_{n}(r,h)=\int_{-h}^{h}\left(\frac{\Phi(r,h)-\Phi(r,z)}{\Phi(r,h)-\Phi(r,0)}\right)^{n}dz.$
(5)

Apparently there is no general expression for such an integral for arbitrary
$n$, but it may be calculated for particular values.

##  4 Generalised Rankine-Hugoniot relation

Let us denote by index 1 the values in TL and by index 2 the values in the
disc. Then the mass continuity condition leads to:

$\displaystyle\rho_{e1}u_{1}\zeta_{n,1}=\rho_{e2}u_{2}\zeta_{n,2}=j,$ (6)

where, for the sake of brevity, we have redefined (5) as
$\zeta_{k,i}=\zeta_{k}(r_{2},h_{i})$. The momentum continuity:

$\displaystyle\rho_{e1}u_{1}^{2}\zeta_{n,1}+p_{e1}\zeta_{(n+1),1}=\rho_{e2}u_{2}^{2}\zeta_{n,2}+p_{e2}\zeta_{(n+1),2}$
(7)

From (6)-(7) we may get the following expression for the mass flux:

$\displaystyle
j^{2}=\frac{p_{e2}\zeta_{(n+1),2}-p_{e1}\zeta_{(n+1),1}}{[\rho_{e1}\zeta_{n,1}]^{-1}-[\rho_{e2}\zeta_{n,2}]^{-1}}.$
(8)

The integrated over height energy continuity condition has the following form:

$\displaystyle
h_{1}u_{1}^{2}+\zeta_{1,1}w_{e1}=h_{2}u_{2}^{2}+\zeta_{1,2}w_{e2},$ (9)

where $w_{e}$ – enthalpy at the equator. From (9) taking into account (6)-(8)
we may derive the generalization of the Rankine-Hugoniot relation:

$\displaystyle\frac{h_{1}[\rho_{e1}\zeta_{n,1}]^{-2}-h_{2}[\rho_{e2}\zeta_{n,2}]^{-2}}{[\rho_{e1}\zeta_{n,1}]^{-1}-[\rho_{e2}\zeta_{n,2}]^{-1}}=\frac{\zeta_{1,2}w_{e2}-\zeta_{1,1}w_{e1}}{p_{e2}\zeta_{(n+1),2}-p_{e1}\zeta_{(n+1),1}}.$
(10)

If we set homogeneous vertical distribution
$\zeta_{k,i}=\zeta_{k}(r_{2},h_{i})=2h_{i}$ and $h_{1}=h_{2}$ then the usual
expression may be gotten:

$\displaystyle\frac{\rho_{e2}^{-1}+\rho_{e1}^{-1}}{2}=\frac{w_{e2}-w_{e1}}{p_{e2}-p_{e1}}.$
(11)

For a given equation of state and known values of thermodynamic quantities on
both sides of the discontinuity, (10) may be considered as an implicit
equation relating the heights $h_{1}$, $h_{2}$ and the radius $r_{2}$.

Figure 2: The dependence of the ratio of the TL height to that of disk
$h_{1}/h_{2}$ on the Reynolds number $\gamma_{2}={\dot{M}}/{4\pi
h_{2}\eta_{2}}$ and temperature in the TL $T_{1}$ for the ratio of turbulent
viscosities $\eta_{2}/\eta_{1}=45$. All other calculations were carried out
along the diagonal line, which reflects the transition between the high/soft
and low/hard states.

Figure 3: The dependence of the TL length $L$ (left) and the ratio (right) of
the TL height to disk height $h_{1}/h_{2}$ on the Reynolds number,
$\gamma_{2}={\dot{M}}/{4\pi h_{2}\eta_{2}}$ for different ratios of turbulent
viscosities $\eta_{2}/\eta_{1}$: 60 (solid line), 45 (dashed) and 30 (dotted).

##  5 Rotation

The continuity equation, integrated over the height has the following
solution:

$\displaystyle r\rho_{e}u\zeta_{n}=-q,$ (12)

where $q=\dot{M}/2\pi>0$ is a constant. The equation of force balance in
$\phi$ direction for the TL is the following:

$\displaystyle\frac{\partial}{\partial
r}(r^{2}\rho_{1}u_{1}r\omega_{1}-\eta_{1}r^{3}\omega^{\prime}_{1})=0.$ (13)

Integrated over the height it has the solution
$\omega_{1}=C_{1}r^{-\gamma_{1}}+C_{2}r^{-2}$, where we denoted

$\displaystyle\gamma_{1}=\frac{q}{2h_{1}\eta_{1}}=\frac{\dot{M}}{4\pi
h_{1}\eta_{1}},$ (14)

which is nothing more than the Reynolds number. Before defining the unknown
constants, let us discuss the boundary conditions.

In our model, we do not take into account the gradual broadening of the flow
after the shock wave, nor its gradual narrowing near the CO. Therefore,
neglecting this, we require equality of angular velocities to the disk one at
the outer boundary and the CO one at the inner boundary:

$\displaystyle\omega_{1}(r_{1})=\omega_{s},$ (15)
$\displaystyle\omega_{1}(r_{2})=\omega_{K}(r_{2}).$ (16)

In addition to the conservation of radial momentum at the discontinuity, there
must also be conservation of momentum in the $\phi$-direction. Given the
viscosity, the law of $\phi$-momentum continuity, integrated over height,
reads as follows:

$\displaystyle
j\,\omega_{1}-2h_{1}\eta_{1}\omega^{\prime}_{1}=j\,\omega_{2}-2h_{2}\eta_{2}\omega^{\prime}_{2},$
(17)

where $j$ is the mass flux from (6).

Since we suppose that at the outer boundary $r=r_{2}$ both angular velocities
are equal to the Keplerian one (16), then we have the following condition for
the angular velocity derivative:

$\displaystyle\omega^{\prime}_{1}(r_{2})=\frac{\eta_{2}}{\eta_{1}}\frac{h_{2}}{h_{1}}\omega^{\prime}_{K}(r_{2}).$
(18)

Thus, for the second-order differential equation (13) we have three boundary
conditions (15)(16)(18), i.e. the problem is overdetermined. Using the
conditions (15)(16) we may write down the solution:

$\displaystyle\omega_{1}=\frac{1}{r_{1}^{2}r_{2}^{\gamma_{1}}-r_{1}^{\gamma_{1}}r_{2}^{2}}\left(\frac{r_{1}^{2}r_{2}^{2}}{r^{2}}(\omega_{\text{K}}r_{2}^{\gamma_{1}}-\omega_{s}r_{1}^{\gamma_{1}})-\frac{r_{1}^{\gamma_{1}}r_{2}^{\gamma_{1}}}{r^{\gamma_{1}}}(\omega_{\text{K}}r_{2}^{2}-\omega_{s}r_{1}^{2})\right),$
(19)

and (18) may be used to get another implicit equation for $h_{1}$, $h_{2}$ and
$r_{2}$ with given $\omega_{s}$, $\omega_{K}(r_{2})$ and $r_{1}$. Thus,
knowing the CO radius $r_{1}$ and height of the accretion disk $h_{2}$, as
well as the thermodynamic characteristics of the plasma $\rho_{e}$, $p_{e}$,
$w_{e}$ on both sides of the discontinuity, the length, and height of the
transition layer can be calculated using (10)(18)(19).

##  6 Setup and results

Since we aim to calculate the geometric characteristics $h_{1}$, $r_{2}$ of
the TL, in (19) we have switched to the Reynolds number for the disk
$\gamma_{2}$, which is independent of $h_{1}$:

$\displaystyle\gamma_{2}=\gamma_{1}\left(\frac{\eta_{1}}{\eta_{2}}\right)\left(\frac{h_{1}}{h_{2}}\right)=\frac{\dot{M}}{4\pi
h_{2}\eta_{2}}.$ (20)

As the CO we considered an NS with a mass of $M=1.5M_{\odot}$, a radius of
$r_{1}=15\text{ km}$ and an angular velocity $\omega_{s}=10^{-2}\text{
s}^{-1}$. The gravitation potential of the NS is the Newtonian,
$\Phi=-GM/(r^{2}+z^{2})^{1/2}$.

Although in deriving the vertical distributions (3)-(4) it was assumed that
the gas is polytropic, the equation of state for the quantities at the equator
can, generally speaking, be given arbitrarily; but we considered the same
equation of state for vertical and radial dependences. For both TL and
accretion disk we took the polytropic index $n=3$, which corresponds to the
photon gas equation of state with $w=4p/\rho$ and $p=\alpha T^{4}/3$, where
$\alpha$ – radiation density constant, $T$ – temperature.

As a disk height, we took three Schwarzchild radii $h_{2}=3r_{\rm
S}=13.3\text{ km}$. The maximum (at the equator) density and temperature in
the disk were taken as $\rho_{e2}=1.2\cdot 10^{-5}\text{ g cm}^{-3}$ and
$T_{e2}=0.5\text{ keV}$ see e.g. ST09.

The density for the TL was chosen to correspond to the observed optical depth
relative to Thompson scattering $\tau\simeq 2$ (see e.g. ST09). This value can
be achieved by selecting the concentration $10^{17}\text{ cm}^{-3}$ which
gives $\rho_{e2}=2\cdot 10^{-7}\text{ g cm}^{-3}$.

For neutron stars, the TL temperature is around $25$ keV in the low/hard state
and evolves to approximately $2$ keV in the high/soft state (e.g. Seifina &
Titarchuk (2010)). In our analysis, there is no way that we can establish a
law of temperature change as a function of the mass accretion rate. Therefore,
to begin with, we solved (10)(18)(19) for the entire temperature range. As can
be seen from Fig. 2 the height of the TL is almost independent of the
temperature in it. The same picture is observed for the TL length. However,
the solution of (10) with some fixed $r_{2}$ alone shows considerable
sensitivity to the temperature. The independence of the geometry of the TL on
its temperature obtained in our model is achieved only by solving the
equations for height and length together. Thus, we can conclude that the
geometric characteristics of the TL depend mainly on the accretion rate
(Reynolds number $\gamma$). Hereafter, these geometric characteristics have
been calculated assuming a simple linear relation, $T_{1}(\gamma_{2})$,
despite the insignificance of this correction. The relation,
$T_{1}(\gamma_{2})$, represented in Fig. 2, reflects the described above
transition between the high/soft and low/hard states.

The turbulent viscosity may be estimated in different ways. SS73 proposed its
value as $\eta=\alpha\rho c_{s}l_{\text{turb}}$, where $\alpha$ is a
dimensionless constant, $c_{s}$ is the sound speed, and $l_{\text{turb}}$ is a
turbulent length scale, which is equal to the height $h$ or, in addition, to
the value $(h^{-2}+h_{r}^{-2})^{-1/2}$, where $h_{r}$ is the radial pressure
scale $h_{r}=|p/p^{\prime}|$ (Popham & Sunyaev (2001), PS01 hereafter). Also
in radiation-pressure dominated conditions it may be estimated (TLM98) as
$\eta=m_{p}n_{\text{ph}}cl/3$, where $m_{p}$ is the proton mass,
$n_{\text{ph}}$ is the photon number density, $c$ is the speed of light, and
$l$ is the photon mean free path. So the turbulent viscosity may depend on the
conditions in the flow as well as its geometry. We avoided further reduction
and took several values of the ratio of turbulent viscosities
$\eta_{2}/\eta_{1}$ to make the resulting TL lengths $L$ to be near the
observable ones as well as $\gamma_{2}$ corresponds to $\alpha^{-1}$ parameter
from the SS73 model.

We calculated (10)(18)(19) with parameters describes above for three different
values of the ratio $\eta_{2}/\eta_{1}$: 30, 45, and 60. When this value is
increased, there is no significant change, whereas, at values quite below 30,
the length of TL becomes absurdly large.

As can be seen from Fig. 3, the TL length decreases with increasing
$\gamma_{2}$. At extremely low accretion rates, the length of the TL
considerably exceeds the radius of the CO, and as it increases, it can drop to
tenths of it. The motion of a gas element is a superposition of rotation with
an angular velocity $\omega$ and falling on the CO with the velocity $u$, i.e.
it is a spiral motion. A smaller radial gas flow (i.e. smaller $\dot{M}$ and
$\gamma_{2}$) increases the distance between the coils of such a spiral,
thereby making it larger. Therefore, with a fixed viscosity responsible for
decreasing the angular momentum, increasing the mass accretion rate leads to
decreasing the TL length.

The behavior of the TL height $h_{1}$ is similar – it is about 80 times the
height of the disk at a low accretion rate and drops to about 45 when the
accretion rate increases (see Fig. 3). Thus, the height of the TL of an NS can
be up to a thousand kilometers in the low/hard state and half that in the
high/soft state. In the frame of our model, this result is quite explainable.
The fact is that an increase in the mass accretion rate at a fixed density
$\rho_{2}$ implies an increase in radial velocity. In this way, more matter
flows through a certain area per unit of time, making it smaller.

##  7 Conclusion

We have solved the problem of the transition layer geometry near a central
object (CO). It was supposed to be weakly magnetized, therefore we study only
the hydrodynamical properties of such configurations. We have derived the
generalized Rankine-Hugoniot relation taking into account different the
heights of the accretion disk and TL and heterogeneous distributions along
them. Then we solved the equations of continuity and force balance in $\phi$
direction and have shown that the last one ought to have three boundary
conditions. In order to satisfy them, a TL must have certain geometric
characteristics.

As a CO we considered NS, setting the gas parameters on both sides of the
shock wave according to observations (see ST09). Then the length and height of
the transition layer were calculated. We suggest that the turbulent viscosity
is constant on each side of the shock, which values have been chosen as such
that the Reynolds number $\gamma_{2}$ corresponds to $\alpha^{-1}$ parameter
from the SS73 model. As can be seen from Fig. 3, both the TL length and height
decrease with increasing $\gamma_{2}$. Therefore we may conclude that:

1\. We clarify the nature of the region where X-ray spectra are formed.

2\. When the source evolves to a softer state, the Compton corona region
becomes more compact (see Fig. 3 and also ST06, ST09). We should emphasize
that the integrated power $P_{x}$ of the resulting power density spectra
rapidly declines toward soft states (TSA07).

3\. We should point out that a ratio of the TL height to the disk height
$h_{1}/h_{2}$ on the Reynolds number is quite large (see Fig. 3).

4\. The calculations were done for a NS, although the same behavior should be
expected for other types of COs (a WD and a BH).

5\. According to the presented calculations and observed manifestations (ST06,
ST09) the TL height $h_{1}$ sufficiently exceeds the disk height $h_{2}$.
Therefore, as also noted in PS01, one-dimensional constant-height models need
to be improved to correctly describe the transition layer.

## References

  * Abramowicz & Fragile (2013) Abramowicz, M. A. & Fragile, P. C. 2013, Living Reviews in Relativity, 16, 1
  * Ajay et al. (2022) Ajay, A., Rajesh, S. R., & Singh, N. K. 2022, arXiv e-prints, arXiv:2210.04997
  * Chakrabarti & Titarchuk (1995) Chakrabarti, S. & Titarchuk, L. G. 1995, ApJ, 455, 623
  * Chakrabarti (1989) Chakrabarti, S. K. 1989, ApJ, 337, L89
  * Das et al. (2022) Das, P., Porth, O., & Watts, A. L. 2022, MNRAS, 515, 3144
  * Farinelli & Titarchuk (2011) Farinelli, R. & Titarchuk, L. 2011, A&A, 525, A102
  * Hōshi (1977) Hōshi, R. 1977, Progress of Theoretical Physics, 58, 1191
  * Lančová et al. (2019) Lančová, D., Abarca, D., Kluźniak, W., et al. 2019, ApJ, 884, L37
  * Liu & Qiao (2022) Liu, B. F. & Qiao, E. 2022, iScience, 25, 103544
  * Lugovskii & Chechetkin (2012) Lugovskii, A. Y. & Chechetkin, V. M. 2012, Astronomy Reports, 56, 96
  * Matsumoto et al. (1984) Matsumoto, R., Kato, S., Fukue, J., & Okazaki, A. T. 1984, PASJ, 36, 71
  * Mishra et al. (2020) Mishra, B., Kluźniak, W., & Fragile, P. C. 2020, MNRAS, 497, 1066
  * Popham & Sunyaev (2001) Popham, R. & Sunyaev, R. 2001, The Astrophysical Journal, 547, 355
  * Romanova et al. (2021) Romanova, M. M., Koldoba, A. V., Ustyugova, G. V., et al. 2021, MNRAS, 506, 372
  * Seifina & Titarchuk (2010) Seifina, E. & Titarchuk, L. 2010, ApJ, 722, 586
  * Shakura & Sunyaev (1973) Shakura, N. I. & Sunyaev, R. A. 1973, A&A, 500, 33
  * Shaposhnikov & Titarchuk (2009) Shaposhnikov, N. & Titarchuk, L. 2009, ApJ, 699, 453
  * Strohmayer et al. (1996) Strohmayer, T. E., Zhang, W., Swank, J. H., et al. 1996, ApJ, 469, L9
  * Suková et al. (2017) Suková, P., Charzyński, S., & Janiuk, A. 2017, MNRAS, 472, 4327
  * Sunyaev & Titarchuk (1980) Sunyaev, R. A. & Titarchuk, L. G. 1980, A&A, 86, 121
  * Titarchuk et al. (1998) Titarchuk, L., Lapidus, I., & Muslimov, A. 1998, ApJ, 499, 315
  * Titarchuk et al. (2007) Titarchuk, L., Shaposhnikov, N., & Arefiev, V. 2007, ApJ, 660, 556
  * van der Klis et al. (1996) van der Klis, M., Swank, J. H., Zhang, W., et al. 1996, ApJ, 469, L1
  * Zhang et al. (1996) Zhang, W., Lapidus, I., White, N. E., & Titarchuk, L. 1996, ApJ, 473, L135

Make science, not war


# LGR-MPC: A user-friendly software based on Legendre-Gauss-Radau pseudo
spectral method for solving Model Predictive Control problems

Saeid Bayat  James T.Allison

###### Abstract

Active components, such as actuators, constitute a fundamental aspect of
engineering systems, affording the freedom to shape system behavior as
desired. However, this capability necessitates energy consumption, primarily
in the form of electricity. Thus, a trade-off emerges between energy usage and
desired outcomes. While open-loop optimal control methods strive for
efficiency, practical implementation is hampered by disturbances and model
discrepancies, underscoring the need for closed-loop controllers. The
Proportional-Integral-Derivative (PID) controller is widely favored in
industry due to its simplicity, despite sub-optimal responses in many cases.
To bridge this gap, Model Predictive Control (MPC) offers a solution, yet its
complexity limits its broad applicability. This paper introduces user-friendly
Python-based MPC software, enabling easy access to MPC. The effectiveness of
this software is demonstrated through multiple examples, including those with
a known analytical solution

###### keywords:

Model Predictive Control (MPC) , Pseudo Spectral Method , Optimization ,
Software

[author-ise] department = Department of Industrial and Enterprise Systems
Engineering,, organization = University of Illinois at Urbana-Champaign,,
city=Urbana, state=IL, country=USA

##  1 Introduction

Actuators, known as active components, are commonly employed in engineering
systems to fulfill the objectives set by the designer. By manipulating the
actuators, the system’s response can be fine-tuned to achieve optimal outcomes
and stabilize unstable systems. However, the operation of these components
necessitates energy consumption, resulting in a trade-off between the actions
performed and energy usage. Designers strive to attain the best possible
results while minimizing energy consumption. To accomplish this objective,
Open Loop Control (OLC) problems are formulated, where no structural
assumptions are made, enabling the optimization algorithm to generate the
optimal control signal freely. While OLC offers advantages, it faces a
significant limitation: it cannot be practically implemented due to the
disparities between the real system and the model used for OLC design [1].
Thus, the adoption of Closed Loop Control (CLC) becomes necessary, utilizing
system state feedback. CLC controllers excel at handling disturbances, model
mismatches, measurement noise, and other factors, making them widely applied
across various applications.

The industry widely recognizes Proportional-Integral-Derivative (PID) control
as the most renowned CLC method, thanks to its favorable response and
straightforward design [2]. However, the PID response may significantly differ
from that of OLC because PID has a specific structure that prevents it from
replicating an OLC-like response, even with the assistance of optimization
algorithms. This disparity between OLC and CLC can be bridged by Model
Predictive Control (MPC). MPC offers varying levels of information, allowing
for the inclusion of controllers ranging from those with no knowledge of
future events to those equipped with complete information. Despite its
advantages, MPC faces a significant drawback in terms of computational cost
[3, 4]. However, recent advancements in transistor and processor design have
expanded the applicability of MPC in engineering systems [5, 4, 6]. Another
challenge lies in the design of MPC itself, particularly for industry
professionals whose expertise may not primarily lie in control theory. To
address this, user-friendly software becomes essential in making MPC a more
versatile controller. This paper introduces a solution where users are not
required to possess in-depth knowledge of the MPC design process. By
leveraging MPC hyper parameters, users can effectively design an MPC for their
specific study system without delving into intricate details.

The concept of developing closed-loop control based on an open-loop controller
was first introduced in literature more than 50 years ago. Notably, Lee and
Markus [7] presented the essence of Model Predictive Control (MPC). In the
present day, with advancements in transistor technology and the development of
high-frequency processors, MPC finds applications in a wide range of areas.
These include suspension systems [8], automotive powertrains [9], power
converters [10], vehicle traction control [11], and various other domains.

Several frameworks have been introduced in the literature to address Model
Predictive Control (MPC) challenges. Notable examples include DO-MPC [12],
ACADO [13], and MUSCOD II [14]. While these tools offer solutions, many of
them demand specific software dependencies or necessitate defining problems
using particular syntax, which might not be user-friendly. For instance, in
the case of MUSCOD II [14], users are required to submit files without access
to the underlying code, which is also not publicly available. On the other
hand, DO-MPC [12] provides an open-access modular codebase, allowing users to
define and modify different modules as needed. However, the usage of CASADI
[15] for automatic differentiation poses a challenge. Users must formulate
problems using CASADI’s syntax, which has a steep learning curve initially.
Furthermore, DO-MPC [12] assumes the control horizon to be the same as the
prediction horizon, limiting flexibility. Additionally, due to the CASADI
dependency, working with black-box functions becomes difficult or even
impractical. This is problematic because many engineering systems involve
black-box functions in dynamics, constraints, or objective functions.

Several tools, including DO-MPC [12], have publicly available codebases.
However, the complexity of these codes hinders their educational utility,
making it challenging to comprehend the functioning of each component for
teaching purposes. Moreover, certain techniques like single shooting employed
in tools like SS-MPC [4] suffer from issues such as ill-conditioning.
Furthermore, specific MPC toolsets are tailored for particular tasks like
building operations. Examples of such task-specific software include MShoot
[16], MPCPY [17], and TACO [18], all designed exclusively for this purpose.

To address the aforementioned challenges, a novel MPC toolset is introduced in
this paper, characterized by the following key attributes:

  * •

Open Source: The toolset is made publicly available, with access to the
underlying source code. This encourages broader utilization, including
educational purposes, and provides transparency.

  * •

Independence from Third-Party Software: This toolset operates autonomously,
without necessitating the installation of additional third-party software.
This alleviates concerns about compatibility across different machines.

  * •

Simplified Syntax: Unlike tools relying on specific syntax like CASADI [15],
this toolset is developed in Python, requiring only a basic understanding of
Python for effective usage.

  * •

Black Box Function Compatibility: Given that this paper does not rely on
automatic differentiation, it becomes possible to employ any black box
function seamlessly. Therefore, users are not compelled to either rephrase
their black box functions using a predefined syntax or modify their problem to
convert the black box function into a symbolic form.

  * •

Stability through Pseudospectral Methods: Employing pseudospectral methods,
the toolset transforms the problem into a nonlinear program, mitigating the
ill-conditioning problems associated with certain methods like Shooting [4].

  * •

Comprehensibility: The paper elaborates on the essential steps involved in
formulating MPC as a Nonlinear Programming (NLP) problem. This offers a
detailed understanding of the toolset’s inner workings.

  * •

Flexibility: Each aspect, including dynamics, objectives, constraints, MPC
parameters, and optimization settings, can be modified easily without
necessitating the complete redefinition of the problem. This flexibility
empowers users to assess the impact of each parameter on the solution.
Additionally, control horizon and prediction horizon can be defined
independently, offering greater control in MPC design compared to tools like
DO-MPC [12].

  * •

Simplified Design: The toolset eliminates the need for users to manually
implement the MPC loop and individual function calls for tasks such as
simulation and result retrieval. The code automates these steps, providing
users with a vectorized result across the entire time horizon.

  * •

Active GitHub Repository: The toolset is backed by an active GitHub
repository, ensuring ongoing updates and support for addressing reported
issues.

  * •

Accessible Complexity: Although toolsets like DO-MPC [12] offer modular
frameworks applicable in various scenarios, such as robust MPC or state
estimation, their intricacy often obscures the understanding of how distinct
code components operate. In contrast, this toolset aims to strike a balance.
It endeavors to provide software that remains relatively uncomplicated while
establishing a robust foundation built on advanced techniques, such as
pseudospectral methods, for designing MPC problems.

It’s essential to note that this toolset is not intended to replace existing
powerful MPC toolsets. Instead, its purpose is to specifically address the
outlined goals, offering a solution that aligns with the aforementioned
attributes. The remainder of this article is as follows: In Sec. 1.1 through
1.2, you will find a comprehensive overview of MPC parameters and Pseudo-
Spectral methods. In Sec. 2 the MPC implementation is discussed. Software
structure is discussed in Sec. 3 and different examples are solved in Sec. 4.
Finally, Sec. 5 provides a concluding summary. In A, a variety of methods for
resolving MPC problems are discussed, and B contains three distinct examples
that have been solved using this software.

###  1.1 MPC

To design an MPC controller, three hyper parameters are needed:

  * •

Control Sampling Time ($T_{s}$): is the rate at which the controller executes
the control algorithm. At the time interval between each sampling point, the
control is constant and is equal to the previous time step. It also determines
how fast the controller can react to disturbances and setpoint changes. In
addition, it determines how advanced processors are needed because as $T_{s}$
decreases, the computation time increase, but with the benefit of responding
faster to disturbances and setpoint changes. The rule of thumb to determine
$T_{s}$ is to put 10 to 20 samples in system rise time, which is defined as
the time that takes for the system step response to rise from $10\%$ to $90\%$
of the steady-state response [19].

  * •

Prediction Horizon ($p$): indicates the number of time steps utilized for
making predictions using the system model. The product $pT_{s}$ represents the
duration for which the prediction is made. Increasing the value of $p$ leads
to more accurate MPC results as it provides more information about the system.
However, this comes at the expense of higher computation costs. To strike the
right balance, it is advisable to include sufficient samples to cover the
system’s settling time. The settling time refers to the duration when the
error between the system’s step response and its steady-state becomes less
than or equal to $2\%$ [20].

  * •

Control Horizon ($m$): $m\leq p$ shows the number of time steps for which MPC
calculates optimal actions. The maximum value for $m$ is $p$, and the length
of the time window in which control is not constant is $mT_{s}$. In the time
interval $[mT_{s}\,\,pT_{s}]$ control is constant and is equal to control
action at time $mT_{s}$. A general rule to set $m$ is to use a value between
$0.1p$ and $0.2p$ with a minimum of 3 [21].

MPC addresses a problem that is similar to optimal control in an open-loop
setting but with an iterative approach. In each iteration of MPC, the
controller solves an optimization problem to determine the optimal control
sequence for a finite time horizon. This includes predicting the system’s
future behavior based on a model and optimizing the control inputs to minimize
a specified cost function. The iterative nature of MPC allows for the
incorporation of real-time measurements and the ability to adapt the control
actions based on the evolving system dynamics, disturbances, and changing
setpoints. By continuously updating the control inputs, MPC effectively
addresses the dynamic nature of the system, leading to improved performance
and robustness compared to open-loop optimal control. The general optimal
control problem is shown in Eq. 1 [22, 23]. Here $t$ is time, $t_{0}$ is
initial time, $t_{f}$ is final time, $u$ is action, $\bm{\xi}$ shows state,
$\Phi$ is Mayer cost, $\mathcal{L}$ is Lagrange cost, $\mathcal{\bm{C}}$ shows
path constraint and $\mathcal{\bm{\phi}}$ shows boundary constraint. Dynamis
is also shown by $\bm{f}$. It should be noted that depending on the method
used to solve the optimal control problem, we have a different notation for
dynamic in the OLC formula. This will be discussed in the next section.

$\displaystyle\min_{\bm{u}(t)}\,\,\,$
$\displaystyle\Phi\left(\bm{\xi}(t_{0}),t_{0},\bm{\xi}(t_{f}),t_{f}\right)+\int^{t_{f}}_{t_{0}}\mathcal{L}\left(t,\bm{\xi}(t),\bm{u}(t)\right)dt$
(1) $\displaystyle\mathrm{Subject\,\,to:}$
$\displaystyle\dot{\bm{\xi}}-\bm{f}\left(t,\bm{\xi}(t),\bm{u}(t)\right)=0$
$\displaystyle\mathcal{\bm{C}}\left(t,\bm{\xi}(t),\bm{u}(t)\right)\leq 0$
$\displaystyle\mathcal{\bm{\phi}}\left(t_{0},\bm{\xi}(t_{0}),t_{f},\bm{\xi}(t_{f})\right)=0$

The distinction between MPC and open-loop optimal control problems lies in the
iterative nature of MPC. In each iteration, MPC solves a problem defined by
Eq. 1, where $t_{0}$ and $t_{f}$ represent the starting and final times for
that particular iteration. After each iteration, the starting time moves
forward by one sampling time $T_{s}$, while the final time is adjusted to
ensure that the time interval between the initial and final times covers the
entire prediction horizon. The MPC loop, illustrated in Fig. 1, demonstrates
how each iteration minimizes the objective over the corresponding time window
while satisfying dynamic equations, path constraints, and boundary conditions.
Once convergence is achieved, the resulting control actions are applied to the
real system, denoted as $\bm{f}$, and the iteration count is increased by 1.
It is important to note that the Mayer cost and boundary conditions are only
considered when the initial or final time of the system falls within the time
window being solved by MPC. During each iteration, an optimization algorithm
is employed to solve the MPC problem. The model used in the optimization
process, denoted as $\hat{\bm{f}}$, can differ from the actual system model.
For instance, a simplified model might be used for faster computation, or
disturbances may be unknown, resulting in the optimization model not
accounting for them.

Figure 1: Schematic illustrating MPC loop. At each iteration, the optimization
problem is solved using a nonlinear programming software while satisfying
constraints and system model dynamics ($\hat{\bm{f}}$). When it converges, the
first time step of the controller is applied to the real plant ($\bm{f}$).

MPC problem shares similarities with the open-loop optimal control problem. A
delves into various methods employed to address open-loop optimal control
problems. This paper uses the pseudo-spectral method to solve the MPC
problems. Pseudo-spectral methods can be divided into different categories
depending on the nodes used for collocation and regression. This method is
discussed in detail in Sec. 1.2.

###  1.2 Pseudo spectral Methods

In pseudospectral methods, as discussed in A, the states and controls are
discretized. High-order polynomials are employed to approximate the continuous
function of states, while no assumptions are made for control signals at
points other than the discretized ones. Collocation is used in the
pseudospectral method to satisfy the dynamic equation. Specifically, a set of
points, denoted as $\tau={\tau_{1},\tau_{2},\cdots,\tau_{N}}$, is employed to
equate the left-hand side of the dynamic equation to the right-hand side. This
results in the generation of defect constraints, represented by $\zeta$:

$\displaystyle\dot{\bm{\xi}}(\tau_{j})-\bm{f}(\tau_{j},\bm{\xi}(\tau_{j}),\bm{u}(\tau_{j}))=0\,\,,(j=1,\cdots,N)$
(2)

In pseudo-spectral methods, the collocation points ($\tau$) includes 3 main
categories: Gauss methods, Radau method and Lobatto methods. In a gauss
method, neither of the end points ($\tau_{1},\,\tau_{N}$) are collocation
points. In Radao method at least on of the end points
($\tau_{1},\,\mathrm{or},\,\tau_{N}$) is a collocation point. In a Lobatto
method, both end points ($\tau_{1},\,\mathrm{and},\,\tau_{N}$) are collocation
points.

In the pseudo-spectral method, the set of points used for collocation is
obtained from the root of the Chebyshev or Legendre polynomial [24]. These
points provide orthogonal collocation property. The benefit of this is that
the quadrature approximation to a definite integral is pretty accurate [24].
When roots of Legendre polynomial are used, these methods are called Legendre-
Gauss(LG), Legendre-Gauss-Radau(LGR), and Legendre-Gauss-Lobatto (LGL). All
these points are defined on the interval $[-1,+1]$, So the time domain
$[t_{0},t_{f}]$ is mapped to $[-1,1]$ by using the following equation:

$\displaystyle t=\frac{t_{f}-f_{0}}{2}\tau+\frac{t_{f}+t_{0}}{2}$ (3)

In these method, the $N$ collocation points are obtained from Eq. 4. Figure 2
also shown the location of these points. As we see, LGL contains both
endpoints, LGR contains just starting point, Flipped LGR includes end point,
and LG contains no end point.

$\displaystyle\mathrm{LG}:\,\,\mathrm{Roots\,\,of\,\,}P_{N}(\tau)$ (4)
$\displaystyle\mathrm{LGR}:\,\,\mathrm{Roots\,\,of\,\,}P_{N-1}(\tau)+P_{N}(\tau)$
$\displaystyle\mathrm{LGL}:\,\,\mathrm{Roots\,\,of\,\,}\dot{P}_{N-1}(\tau)+\\{-1\\}$
$\displaystyle\mathrm{Flipped\,\,LGL}:\,\,\mathrm{Roots\,\,of}\dot{P}_{N-1}(\tau)+\\{+1\\}$

Figure 2: Collocation nodes of different pseudo-spectral methods

As mentioned earlier, the open-loop optimal control problem will be
transcribed to an NLP using direct methods. As a result, The equivalence of
the NLP and the OLC is essential. In Rao [24], it is shown that the LGL
methods do not yield a complete equivalence, but LGR and LG methods do. In
addition, because the LGR method has a collocation point at the start, it is
better than the LG method because after solving the NLP, the optimal control
value can be obtained at time 0. In contrast, if the LG method is used, there
is no value at time 0 for state or control. Therefore, the LGR method is used
in this paper.

The procedure of LGR method is as follows: Consider N LGR points,
$\tau=\\{\tau_{1},\cdots,\tau_{N}\\}$, where $\tau_{1}=-1$, and $\tau_{N}<+1$.
We define a new point at $\tau_{N+1}=+1$. Then, Lagrange polynomials with
degree N are used to interpolate states at LGR points plus $\tau_{N+1}=+1$, so
we have:

$\displaystyle L_{i}(t)=\prod_{k=1,\,k\neq i}^{N}\frac{t-t_{k}}{t_{i}-t_{k}}$
(5)
$\displaystyle\bm{\xi}(\tau)\approx\bm{\Xi}(\tau)=\sum_{i=1}^{N+1}\bm{\Xi}_{i}L_{i}(\tau)$
(6)

where $\bm{\Xi}_{i}=\bm{\Xi}(\tau_{i})$. As a result, each state is
approximated with a Lagrange polynomial with degree N. The collocation points
are the N LGR points and are used to provide defect constraints that make the
left-hand side of the dynamic equal to the right-hand side. By taking the
derivative of the above equation, we have [25]:

$\displaystyle\bm{\dot{}}{\xi}(\tau)\approx\bm{\dot{}}{\Xi}(\tau)=\sum_{i=1}^{N+1}\bm{\Xi}_{i}\dot{L}_{i}(\tau)$
(7)

By equating the derivative of state to the dynamic for the LGR points, we have
[25]:

$\displaystyle\sum_{i=1}^{N+1}\bm{\Xi}_{i}\dot{L}_{i}(\tau_{k})$
$\displaystyle=\frac{t_{f}-t_{0}}{2}\bm{f}\left(\bm{\Xi}_{k},\bm{U}_{k},\tau,t_{0},t_{f}\right),\,\,(k=1,\cdots,N)$
(8) $\displaystyle\sum_{i=1}^{N+1}D_{ki}\bm{\Xi}_{i}$
$\displaystyle=\frac{t_{f}-t_{0}}{2}\bm{f}\left(\bm{\Xi}_{k},\bm{U}_{k},\tau,t_{0},t_{f}\right),\,\,(D_{ki}=\dot{L}_{i}(\tau_{k}))$
(9)

Where $\bm{U}_{k}=\bm{U}(\tau_{k})$, and the $D_{ki}$ ,$(1\leq k\leq
N)\,,\,(1\leq i\leq N+1)$ is a $N\times(N+1)$ non-square matrix and is called
differentiation matrix. The D matrix is non-square because approximation of
the states is used at $N+1$ points, but only the LGR points ($N$ points) are
used for collocation. Let $\bm{\Xi}^{\mathrm{LGR}}$ and
$\bm{U}^{\mathrm{LGR}}$ be defined as:

$\displaystyle\bm{\Xi}^{\mathrm{LGR}}$
$\displaystyle=\begin{bmatrix}\bm{\Xi}_{1}\\\ \vdots\\\
\bm{\Xi}_{N+1}\end{bmatrix}_{(N+1)\times n_{\xi}}$ (10)
$\displaystyle\bm{U}^{\mathrm{LGR}}$
$\displaystyle=\begin{bmatrix}\bm{U}_{1}\\\ \vdots\\\
\bm{U}_{N}\end{bmatrix}_{(N)\times n_{u}}$ (11)

Figure 3: Applied constraint to satisfy continuity of states through segments

where $n_{\xi}$ shows number of states and .$n_{u}$ shows number of controls.
As a result, we have:

$\displaystyle\bm{D}_{k}\bm{\Xi}^{\mathrm{LGR}}=\frac{t_{f}-t_{0}}{2}\bm{f}\left(\bm{\Xi}_{k},\bm{U}_{k},\tau,t_{o},t_{f}\right)\,,\,(k=1,\cdots,N)$
(12)

Furthermore, the Lagrange part of the objective can be obtained through
numerical quadrature. As a result, the equivalent NLP form of the OLC is:

$\displaystyle\min_{\bm{\Xi}^{\mathrm{LGR}},\,\bm{U}^{\mathrm{LGR}}}J$
$\displaystyle=\Phi(\bm{\Xi}(\tau_{1}),\tau_{1},\bm{\Xi}(\tau_{N+1}),\tau_{N+1})$
(13)
$\displaystyle\quad+\frac{t_{f}-t_{0}}{2}\sum_{k=1}^{N}w_{k}L(\bm{\Xi}_{k},\bm{U}_{k},\tau,t_{0},t_{f})$
$\displaystyle\,\,\,\,\,\,\mathrm{Subject\,\,to}:$
$\displaystyle\bm{D}_{k}\bm{\Xi}^{\mathrm{LGR}}-$
$\displaystyle\frac{t_{f}-t_{0}}{2}\bm{f}\left(\bm{\Xi}_{k},\bm{U}_{k},\tau,t_{0},t_{f}\right)\,=\,0$
$\displaystyle\bm{\phi}\left(\bm{\Xi}(\tau_{1}),\tau_{1},\bm{\Xi}(\tau_{N+1}),\tau_{N+1}\right)\,\,=\,\,0$
$\displaystyle\frac{t_{f}-t_{0}}{2}\bm{\mathcal{C}}\left(\bm{\Xi}_{k},\bm{U}_{k},\tau,t_{0},t_{f}\right)\leq$
$\displaystyle\mathrm{for\,\,}k=1,\cdots,N$

In cases where the time window $[t_{0},t_{f}]$ is large, it can lead to a
large number of collocation points and a correspondingly high-order
polynomial, which can introduce difficulties. To address this, a strategy
involves partitioning the time window into $S$ segments and employing an
$N^{th}$ degree polynomial within each segment. Let there be $S$ segments,
denoted by $s=1,\cdots,S$, with $N_{s}$ collocation points allocated to each
segment. Additionally, a regression point (at $\tau^{s}_{N+1}=1$) is placed at
the end of each segment. This segmentation introduces additional constraints:
The terminal state value at the current segment must be equal to the initial
state value of the succeeding segment. It’s noteworthy that while the terminal
value doesn’t coincide with a collocation point, the initial value does, as
depicted in Figure 3.

##  2 Implementation

One important difference between MPC and OLC is the time window. In OLC, the
optimization problem is solved through the entire time window $[t_{0},t_{f}]$.
However, in MPC, the problem is solved many time. At each iteration, the time
interval is $[t^{\mathrm{iter}}_{0}\,\,t^{\mathrm{iter}}_{f}]$ . Where,
$t^{\mathrm{iter}}_{0}=(\mathrm{iter-1})T_{s}$, and
$t^{\mathrm{iter}}_{f}=\min(t_{f},t^{\mathrm{iter}}_{0}+pT_{s})$. In addition,
the time window of controller is
$[t^{\mathrm{iter}}_{0c}\,\,t^{\mathrm{iter}}_{fc}]$ . Where,
$t^{\mathrm{iter}}_{0c}=(\mathrm{iter-1})T_{s}$, and
$t^{\mathrm{iter}}_{fc}=\min(t_{f},t^{\mathrm{iter}}_{0c}+mT_{s})$. Figure 4
shows prediction horizon and control horizon time window at each iteration.
Here, it is assumed $m=2$, and $P=3$. As we see, the length of time window may
decrease as it reaches the final time.

Figure 4: Prediction horizon and control horizon through mpc iterations

Another important difference between MPC and OLC is the sample rate, which is
fixed between each sample point at $T_{s}$. As was mentioned before, the
collocation points are obtained by LGR nodes, and the difference between two
consecutive LGR nodes may be smaller than $T_{S}$. Therefore, the control at
those points should not be considered as a design variable but should be fixed
and equal to previos node that has satisfied these criteria. This concept is
shown in Fig. 5. The final time of the control horizon is
$(t^{\mathrm{iter}}_{0c}+mT_{s})$, so as we see, the control at a time beyond
this point is fixed and is equal to the control signal at this final point. In
times smaller than $(t^{\mathrm{iter}}_{0c}+mT_{s})$, this criteria is checked
on collocation points. If the length of time window between two LGR points is
less than $T_{s}$, it is not considered a control variable but is fixed. These
points are shown by hollow shapes. Other points are shown by filled shapes. In
addition, the control signal between these control points is constant and
equal to the control variable at that point.

Figure 5: Control variables at each iteration that satisfy $T_{s}$ sampling
rate

At each iteration, the optimization variables include states and controls.
States are defined at LGR points and end points of each section, which are
used for regression. However, for the control signal, those collocation points
in the control time window that satisfy $T_{s}$ criteria are considered as
control variables. Fig. 6 shows the optimization variable at each iteration.
$\bm{X}^{\mathrm{iter}}$ consists of many blocks. Each block corresponds to
state and control variables at each segment.

Figure 6: MPC variables at each iteration which includes both state and
control values

As a similar problem is solved at each iteration, the current result can be
used as a guess for the next. This is shown in Fig. 7. If the success flag of
the optimization algorithm is “True”, then the result is considered as the
initial guess. However, if it is not true, then a signal that varies linearly
from lower bound to upper bound is considered for the initial guess.
Additionally, the number of design variables may vary from iteration to
iteration. The code also takes these into account, but they are not shown in
this flowchart for simplicity.

Figure 7: Method to provide initial guess

One thing that is important for optimization algorithms is scaling. Here, all
state and control variables are scaled to $[-1,1]$. To do these, 4 matrices
are considered in the code that should be defined by the user :
$A_{\mathrm{scaled}},B_{\mathrm{scaled}},C_{\mathrm{scaled}},\,\mathrm{and}\,D_{\mathrm{scaled}}$.
$A_{\mathrm{scaled}}\,\mathrm{and}\,B_{\mathrm{scaled}}$ are arrays with
length $n_{\xi}$, and $C_{\mathrm{scaled}}\,\mathrm{and}\,D_{\mathrm{scaled}}$
are arrays with length $n_{u}$.
$A_{\mathrm{scaled}}(i)=\frac{\mathrm{ub}^{i}+\mathrm{lb}^{i}}{2}$,
$B_{\mathrm{scaled}}(i)=\frac{\mathrm{ub}^{i}-\mathrm{lb}^{i}}{2}$,
$C_{\mathrm{scaled}}(j)=\frac{\mathrm{ub}^{j}+\mathrm{lb}^{j}}{2}$,
$D_{\mathrm{scaled}}(j)=\frac{\mathrm{ub}^{j}-\mathrm{lb}^{j}}{2}$, for
$i=1,\cdots\,n_{\xi}$, and $j=1,\cdots,n_{u}$, where “lb” is the lower bound
and “ub” is upper bound.

One more thing that is used in the code is ODE simulation. When pseudo-
spectral methods are used for optimal control, feasibility is one of the
biggest concerns. I.e., the initial mesh may not be adequate to satisfy
dynamic equations and regridding is needed. Here, after each iteration, the
dynamic is simulated in that time interval using obtained control. As a
result, the result is always feasible.

##  3 Software Structure

In this section, the different structures that need to be defined by the user
are defined, and different parts of the code are shown in Fig. 8. Here we go
through them one by one:

(a) (b) (c) (d) (e)

Figure 8: Different structures used in MPC software

  * •

Model Data:
$t_{0}$ is initial time and $t_{f}$ is final time. Scaling matrices are A, B,
C, and D. $\bm{\xi}_{t_{0}}$ and $\bm{\xi}_{t_{f}}$ are states at the initial
time and final time and are used as boundary constraints.
$\mathrm{Flag}_{\bm{\xi}_{0}}$ and $\mathrm{Flag}_{\bm{\xi}_{f}}$ are flags,
which shows whether boundary constraints should be used at the initial and
final point or not. $\bm{\xi}_{\mathrm{lb}}$ and $\bm{\xi}_{\mathrm{ub}}$ are
lower bound and upper bound of states, and $\bm{u}_{\mathrm{lb}}$ and
$\bm{u}_{\mathrm{ub}}$ are lower bound and upper bound of controls.
$\mathrm{Flag}_{\mathrm{Plot}}$ is a flag that shows the plot at each
iteration if the flag is “True”. However, this increases the computation time.
$\mathrm{File\,\,Name}$ is an address to save the result. All other data that
are needed in system dynamics, constraints, or objectives should also be
defined in the Model Data Structure. The Model Data is input to dynamics,
constraints, and objective functions so that the user can get access to these
data through the Model Data structure.

  * •

MPC Parameters:
$T_{s}$ is control sampling time, $m$ is control horizon, and $p$ is
prediction horizon.

  * •

Mesh Parameters:
As was mentioned before, at each iteration, MPC is divided into several
sections and each section has $N$ nodes. Here $\Delta_{t_{\mathrm{segment}}}$
shows time interval for each section. And $\\#\mathrm{Nodes}_{\mathrm{seg}}$
shows number of nodes at each section. As a result, at each iteration the time
interval is divided into some sections base on the division of that time
window to $\Delta_{t_{\mathrm{segment}}}$ and
$\\#\mathrm{Nodes}_{\mathrm{seg}}$ shows number of nodes at each segment.
$n_{h_{\mathrm{sol}}}$ is an integer value and time step that is used to store
ODE result is $\frac{T_{s}}{n_{h_{\mathrm{sol}}}}$.

  * •

NLP Options:
$\mathrm{Iter}_{\mathrm{max}}$ is maximum number of iterations.
$\mathrm{Disp}$ is a flag and if it is “True”, optimization algorithm displays
optimization messages. Here Scipy is used; therefore, Scipy messages will be
displaced.

  * •

Functions:
This structure includes 4 functions: Lagrange cost, Mayer cost, system
dynamics, and path constraints. The input of all these functions are time,
state, control, and Model Data. The output of Lagrange cost is $\mathcal{L}$.
the output of Mayer cost is $\mathcal{M}$, the output of system dynamics is
$\bm{f}$ and the output of path constraint is $\mathcal{C(\cdot)}>0$ which is
a vector and if all elements are equal or grater than zero then the path
constraint is satisfied.

Different constraints used in this code is shown in Fig. 9. These are defined
as:

  * •

Defect Constraint:
This satisfies the dynamic equation at collocation points. The user needs to
define dynamics, and this code itself used that function to generate defect
constraints at all collocation points.

  * •

Path Constraint:
This uses path constraint function to provide path constraint.

  * •

Boundary Constraint:
This uses $\bm{\xi}_{t_{0}}$, $\bm{\xi}_{t_{f}}$,
$\mathrm{Flag}_{\xi_{t_{0}}}$ and $\mathrm{Flag}_{\xi_{t_{f}}}$ to generate
boundary constraints.

  * •

$\Psi_{\mathrm{iter}}$ provide continuity equations over iterations. In other
words, the starting value of current iteration is equal to end value of last
iteration.

  * •

$\Psi_{\mathrm{seg}}$ provide continuity equations over segments. In other
words, the regression point of current segment is equal to first collocation
point of next segment.

Figure 9: Constraint Flowchart

Different objective functions used in this code are shown in Fig. 10. These
are defined as:

  * •

Lagrange cost: This function computes Lagrange cost based on user Lagrange
cost defined in Functions.

  * •

Mayer cost: This function uses user Mayer cost defined in Functions to compute
Mayer cost.

Figure 10: Objective Flowchart

The whole algorithm is shown in the Algorithm. 1. First, it reads user-defined
data and functions. Then it computes all LGR and regression points, goes
through each iteration, and solves the MPC problem. Then it uses the ODE
solver and stores the data, and based on that, it provides an initial point
for states for the next iteration. This is posed as a constraint. In the end,
all data are stored in a pickle file whose name is provided by the user, and
the results are plotted.

Train;

$\mathrm{Define\,Model\,\,Data}$;

$\mathrm{Define\,MPC\,\,Parameters}$;

$\mathrm{Define\,Mesh\,\,Parameters}$;

$\mathrm{Define\,dynamics,\,constraint,\,and\,obj\,funcs}$;

$\mathrm{Define\,Scipy\,Options}$;

for  _$\mathrm{iter}$ $\mathrm{in}$ range($n_{\mathrm{iter}}$)_ do

$t_{0}^{\mathrm{iter}}=(\mathrm{iter}-1)T_{s}$;

$t_{f}^{\mathrm{iter}}=\min(t_{f},t_{0}^{\mathrm{iter}}+pT_{s})$;

$t_{f_{c}}^{\mathrm{iter}}=\min(t_{f},t_{0}^{\mathrm{iter}}+mT_{s})$;

if  _$\mathrm{Flag}_{\mathrm{Plot}}==\mathrm{True}$_ then

$\mathrm{Plot\,\,Results\,\,in\,\,[t_{0}^{\mathrm{iter}},\,t_{f}^{\mathrm{iter}}]}$

end if

$\xi^{\ast},u^{\ast}\leftarrow\mathrm{Solve\,\,MPC\,\,in\,\,interval\,[t_{0}^{\mathrm{iter}},\,t_{f}^{\mathrm{iter}}]}$;

$\mathrm{Store\,\,Data}\leftarrow\mathrm{Run\,ODE\,\,Dynamics\,\,in\,\,[t_{0}^{\mathrm{iter}},\,t_{f_{c}}^{\mathrm{iter}}]}$;

$\mathrm{Update\,\xi_{0}\,for\,the\,next\,iterartion\,uinsg\,ODE}$

end for

$\mathrm{load\,Data\,\,in\,\,pickle\,\,file}$;

$\mathrm{Plot\,Result\,\,in\,\,[t_{0},\,t_{f}]}$;

Algorithm 1 Software Procedure

##  4 Example

In this section, an example with an analytical solution will be explored, and
a comparison will be drawn between the analytical solution and the one derived
through MPC. More exmaples are studied in .B. In these examples the effect of
control sampling time, control horizon, and prediction horizon are
investigated. Also the code in the GitHub repository provides a good examples
how the user should define problems.

The problem studied here is obtained from Bryson and Ho [26] pp.166-167:

$\displaystyle\min_{\bm{u}(t)}\frac{1}{2}\int^{t_{f}}_{t_{0}}u^{2}dt$ (14)
$\displaystyle\mathrm{where:}$ (15)
$\displaystyle\dot{\bm{\xi}}=\begin{bmatrix}0&1\\\
-1&0\end{bmatrix}\bm{\xi}+\begin{bmatrix}0\\\ 1\end{bmatrix}u$ (16)
$\displaystyle\xi_{1}(0)=x_{0},\,\,\xi_{2}(0)=v_{0},\,\,\xi_{1}(t_{f})=0,\,\,\xi_{2}(t_{f})=0$
(17)

The exact open-loop optimal control when $x_{0}=-1/2$ and $v_{0}=1$ is:

$\displaystyle\scriptsize
u^{\ast}(t)=-\frac{2}{t_{f}^{2}-\sin^{2}(t_{f})}\scriptsize\begin{bmatrix}x_{0}\\\
v_{0}\end{bmatrix}^{T}\begin{bmatrix}\sin(t_{f}-t)\sin(t_{f})-t_{f}\sin(t)\\\
-\cos(t_{f}-t)\sin(t_{f})+t_{f}\cos(t)\end{bmatrix}$ (18)

The analytical solution and MPC results are shown in Fig. 11. In the first two
MPC cases, the prediction horizon is the same as the final time ($t_{f}$), so
it has access to the full time-window, and the only difference between these
two cases is in their sampling time: In the first one, $T_{s}=0.2$, and in the
second one $T_{s}=1.0$. In the third case, the time window is $1$ s which is
smaller than the final time ($2$ s), so the MPC has no information about
boundary constraints at the final point until the start time of the MPC is 1
s. As it is shown in Fig. 11(c), the control signal is 0 from start to $1$ s,
and after that, by having access to the boundary constraint, the control
changes to satisfy that constraint.

(a) (b) (c)

Figure 11: Examples 1 response through different scenarios. Here “exact” is
obtained by analytical solution, and other legends show the corresponding MPC
parameters in each scenario.

Using this software, it is possible to modify MPC parameters such as $T_{s}$,
$t_{p}$, and $t_{m}$, enabling an exploration of their impact on the outcomes.
In Fig. 12(a), the variation of the objective value across the 2D space
defined by $T_{s}$ and $t_{p}$ is presented. In this case, the assumption is
that $t_{m}$ is equal to $t_{p}$. Evidently, shifting towards the upper left
corner of the plot leads to diminished objective values, signifying an
improved outcome. This outcome arises because a greater authority is vested in
the MPC controller in this region. Additionally, Fig. 12(b) displays the
objective value across the 2D domain of $T_{s}$ and $t_{m}$, assuming that
$t_{p}=20$. Once again, moving towards the upper left corner results in
enhanced performance. This enhancement is due to the reduction in $T_{s}$ and
the increase in $t_{m}$ along this trajectory. These analyses play a crucial
role in guiding decisions related to actuator, hardware, and sensor selection.
For instance, if there’s no need for excessively small $T_{s}$, more cost-
effective actuators and hardware components can be adopted. Alternatively, if
variations in $t_{m}$ have minimal influence on the outcomes, economical
sensors can be preferred.

(a) (b)

Figure 12: (a) Exploring objective values across a 2D design space of $T_{s}$
and $t_{p}$, assuming $t_{m}=t_{p}$. (b) Exploring objective values across a
2D design space of $T_{s}$ and $t_{m}$, assuming $t_{p}=20.0\mathrm{s}$

##  5 Conclusion

In this article, a user-friendly software is developed to solve MPC problems.
This problem is converted to a nonlinear program using the pseudo-spectral
method, and the optimization problem is solved using Scipy. This code is
written in python, and the corresponding Github address is provided to get
access to the source code and examples. This code was tested under different
examples with exact solutions, and the MPC results were compared with the
exact one. This paper demonstrated how this code works and how different parts
should be defined. Due to its minimal requirement for in-depth MPC
understanding, engineers from diverse disciplines can employ it effectively.

## Declaration of Competing Interest

The authors declare that they have no competing interest.

## Data Availability

All data required to replicate the results can be generated by the Python
optimization code. The Python optimization codes for all the problems
demonstrated in the manuscript are available online using the mentioned GitHub
repository.

##  Appendix A Methods to solvey MPC

  1. 1

Indirect Method: This method is also called “optimize then discretize” where
Pontryagin’s Maximum Principle (PMP) or Calculus of Variation (CV) are used to
derive optimality equations [24] , and then numerical methods are used to
discretize the problem and solve it. The Strength and shortcomings of this
method are:

    1. $+$

: Derived optimality equations gives us insight into the structure of the
solution, ex: Linear Quadratic Regulator (LQR)

    2. $-$

: Obtained boundary value problem is difficult to solve

    3. $-$

: Path constraints are difficult to handle

    4. $-$

: Initial guess for co-state is needed, and this is not straightforward
because co-states do not represent physical states

    5. $-$

: Problem is ill-conditioned and is not robust

    6. $-$

: Cannot be used for black-box functions

  2. 2

Direct Methods: These methods are also called “discretize then optimize”
because the problem is first discretized and then solved using an NLP solver
[25]. The Strength and shortcomings of this method are:

    1. $+$

: Much more effective than direct methods

    2. $+$

: More effective in handling constraints

    3. $+$

: Powerful NLP solver can be used

    4. $+$

: There exist well-established methods

    5. $-$

: Mathematics is not as elegant as indirect methods

    6. $-$

: It gives us no insight into the structure of the problem

Direct methods includes 2 sub-methods: “Sequanrtial” and “Simultaneous”:

    1. 2.1

Sequential: In this method, just the control signal is discretized and dynamic
is satisfied through simulation, and the result is always feasible [23]. As a
result, the dynamic equation in the optimal control problem is not under the
“Subject to” section but will be in the “where” section because it is no
longer satisfied by the optimization algorithm but is obtained through
simulation. The strength and shortcomings of this method are:

      1. $+$

: Smaller number of defect constraints than simultaneous because states are
not discretized

      2. $+$

: Feasible solution can be easily find

      3. $+$

: Powerful differential algebraic solvers can be used

      4. $-$

: Need to perform full simulation for each parameter perturbation

      5. $-$

: Is not effective for highly nonlinear problems or unstable problems

      6. $-$

: It is computationally inefficient

The sequential method itself includes two methods: Single Shooting and
Multiple Shooting [22]:

      1. 2.1.1

Single Shooting: In this method, the whole prediction horizon time window is
simulated in one section. The strength and shortcomings of this method are:

        1. $+$

: Simple

        2. $+$

: Dynamics are always feasible

        3. $+$

: No defect constraint is needed

        4. $-$

: NLP inherits the ill-conditioning of the problem

      2. 2.1.2

Multiple Shooting: In this method, the time window is divided into several
sections, and each section is optimized independently. The states at these
section points must be identical and satisfied through defect constraints.
I.e., the state and final time of the previous section will be the same as the
state at the initial time of the current section. The strength and
shortcomings of this method are:

        1. $+$

: Better convergence than single shooting

        2. $+$

: problem is well-conditioned

        3. $+$

: it has a good structure for parallel computation

        4. $-$

: More complex than single shooting

        5. $-$

: Defect constraint is needed so make the NLP more complex

    2. 2.2

Simultaneous: In this method, both control and states are discretized. In this
approach, the dynamic equation is satisfied through defect constraints. The
strength and shortcomings of this method are:

      1. $+$

: The resulting NLP is sparse, so powerful sparse NLP solvers like IPOPR and
SNOPT can be used

      2. $+$

: Knowledge of the state trajectory can be used as a guess in the
initialization

      3. $+$

: It has fast local convergence

      4. $+$

: Effective for unstable systems

      5. $+$

: Effective for path and boundary constraints

      6. $+$

: Work well even for ill-conditioned problems because the NLP does not inherit
the ill-conditioning of the problem

      7. $-$

: a large number of defect constraints

      8. $-$

: Mesh regrinding is needed and changes the NLP dimension

      9. $-$

: More challenging to find a feasible solution than Sequential metho

      10. $-$

: State derivative function ($\bm{f}$) is needed, so it does not work with
data-driven models

The simultaneous method itself consists of two methods: Single Step and Pseudo
Spectral method:

      1. 2.2.1

Single Step: In this method, the location of discretized points is selected
evenly, and low-order polynomials are used in each section. The strength and
shortcomings of this method are:

        1. $+$

Less difficult to implement than the Pseudospectral method

        2. $-$

Finer mesh than pseudospectral needs to be sued

        3. $-$

more computation cost than pseudospectral

      2. 2.2.2

Pseudospectral: In this method, the location of discretized points are not
placed evenly but are obtained through roots of some polynomials like Legendre
polynomials. Furthermore, high-order polynomials like Lagrange polynomials are
used to interpolate states. The strength and shortcomings of this method are:

        1. $+$

Coarser mesh than single step can be used

        2. $+$

Lead to an easier problem for NLP than single step

##  Appendix B Examples

###  B.1 Example 2

The second problem is obtained from Bryson and Ho [26] pp.120-122.

$\displaystyle\min_{\bm{u}(t)}\int^{t_{f}}_{t_{0}}\frac{1}{2}u^{2}dt$ (19)
$\displaystyle\mathrm{where:}$
$\displaystyle\dot{\bm{\xi}}=\begin{bmatrix}0&1\\\
0&0\end{bmatrix}\bm{\xi}+\begin{bmatrix}0\\\ 1\end{bmatrix}u$
$\displaystyle\xi_{1}(0)=0,\,\,\xi_{2}(0)=v_{0},\,\,\xi_{1}(t_{f})=1,\,\,\xi_{2}(t_{f})=-1$
$\displaystyle\xi_{1}(t)\leq l$

In this example, there are two states and one control. Based on the bounds and
constraints shown in Eq. 19, the scaled matrices can be defined as:

$\displaystyle A_{\mathrm{scaled}}$ $\displaystyle=\begin{bmatrix}(l-l)/2\\\
(2-2)/2\end{bmatrix}$ (20) $\displaystyle B_{\mathrm{scaled}}$
$\displaystyle=\begin{bmatrix}(l+l)/2\\\ (2+2)/2\end{bmatrix}$ (21)
$\displaystyle C_{\mathrm{scaled}}$
$\displaystyle=\begin{bmatrix}(20-20)/2\end{bmatrix}$ (22) $\displaystyle
D_{\mathrm{scaled}}$ $\displaystyle=\begin{bmatrix}(20+20)/2\end{bmatrix}$
(23)

The exact open-loop optimal control when $l=1/9$ is:

$u^{\ast}=\left\\{\begin{array}[]{ll}-\frac{2}{3l}\left(1-\frac{t}{3l}\right)&\mbox{if
$0\leq t\leq 3l$};\\\ 0&\mbox{if $3l\leq t\leq 1-3l$};\\\
-\frac{2}{3l}\left(1-\frac{1-t}{3l}\right)&\mbox{if $1-3l\leq
t$}.\end{array}\right.$ (24)

The analytical solution and MPC results are shown in Fig. 13. In all 3 MPC
cases, the prediction horizon is the same as the final time, and only the
control sampling time is different. The result is closer to the analytical
solution when the control sampling time is shorter. When it is bigger, the
result will deviate from the exact solution obtained by open-loop optimal
control.

(a) (b) (c)

Figure 13: Examples 2 response through different scenarios. Here “exact” is
obtained by analytical solution, and other legends show the corresponding MPC
parameters in each scenario.

###  B.2 Example 3

The third problem is obtained from Bryson and Ho [26] pp.109-110.

$\displaystyle\min_{\bm{u}(t)}\frac{a^{2}}{2}\xi^{2}_{t_{f}}+$
$\displaystyle\int^{t_{f}}_{t_{0}}\frac{1}{2}u^{2}dt$ (25)
$\displaystyle\mathrm{where:}$ (26) $\displaystyle\dot{\bm{\xi}}=b(t)u$ (27)
$\displaystyle\xi(0)=\xi_{0}$ (28) $\displaystyle|u(t)|\leq 1$ (29)

The exact open-loop optimal control when is:

$u^{\ast}=-\mathrm{sat}\left[a^{2}b(t)\xi(t_{f})\right]$ (30)

The analytical solution and MPC results are shown in Fig. 14 for the case
where $t_{f}=1$, $\xi_{0}=1$, $a=1$ and $b(t)=t\cos(20\pi t)-1/4$. Here
control sampling time is the same for both MPC cases, and the only difference
is in their prediction horizon. In the first MPC case, the prediction horizon
covers the full time-window, but in the second case, it only covers $0.2$ s.,
which is smaller than the final time ($1$ s); as a result, in this case, MPC
has no access to the Mayer cost, until it reaches $0.8$ s. As we see, the
control is zero for this case from $t=0$ to $t=0.8$, while for the first case,
the result is close to the exact solution.

(a) (b)

Figure 14: Examples 3 response through different scenarios. Here “exact” is
obtained by analytical solution, and other legends show the corresponding MPC
parameters in each scenario.

###  B.3 Example 4

Fig. 15 shows a simple model of vehicle suspension system defined in Herber
and Allison [27], Bayat and Allison [6]. Here $\delta$ is road profile, $U$ is
unsprung mass, $S$ is sprung mass, $z_{\mathrm{U}}$ is the displacement of
unsprung mass, and $z_{\mathrm{S}}$ is the displacement of the sprung mass.
Furthermore, $F$ shows control force, $m$ shows mass, $k$ shows spring, and
$b$ shows damper. The goal is to provide a smooth ride for passengers by
transferring forces through different components like spring, damper, mass,
and actuator. To reach this goal, the objectives is defined in Eq. 31 [27].
Here $w_{1}(z_{U}-\delta)^{2}$ represents handling performance,
$w_{2}\ddot{z}_{S}^{2}$ shows passenger comfort, and $w_{3}F^{2}$ is a penalty
function for control effort. In this example, the initial value of states is
set to zero, and the displacement between sprung and unsprung mass is bound by
$r_{\mathrm{max}}$. The parameters used in this example are whin in Table. 1.
For the MPC control sampling time, prediction horizon and control horizon are:
$T_{s}=0.01$, $p=50$, and $m=50$, so the length of time window interval is
$0.5$ s.

Figure. 16 shows the results obtained by MPC and OLC. Here, the entire time
horizon is from $0$ to $3$ seconds. Here $\xi_{1}=z_{U}-\delta$,
$\xi_{2}=\dot{z}_{U}$, $\xi_{3}=z_{S}-z_{U}$, $\xi_{4}=\dot{z}_{S}$. In the
first MPC case the prediction horizon is $0.4\,\mathrm{[s]}$ and control
sampling time is $0.01\,\mathrm{[s]}$, and in the second MPC case, the
prediction horizon is $3\,\mathrm{[s]}$ and control sampling time is
$0.1\,\mathrm{[s]}$. As we see, the result of the first case is much closer to
the OLC because of having a shorter sampling time. The OLC response has a high
frequency, so here having a shorter sampling time is more important than
having a bigger prediction horizon.

$\displaystyle\Uppsi_{d}=\int_{t_{0}}^{t_{F}}\left(w_{1}(z_{U}-\delta)^{2}+w_{2}\ddot{z}_{S}^{2}+w_{3}F^{2}\right)\mathrm{d}t$
(31)

Figure 15: Vehicle suspension system Table 1: Vehicle suepnsion parameters

Parameter | Value | Parameter | Value
---|---|---|---
$t_{0}$ | 0 s | $t_{f}$ | 3 s
$b$ |  $100.8$ Ns/m | $k$ |  $2.15*10^{4}$ N/m
$r_{\mathrm{max}}$ |  $0.04$ m | $s_{\mathrm{max}}$ |  $0.04$ m
$m_{\mathrm{U}}$ |  $65$ kg | $m_{\mathrm{S}}$ |  $325$ kg
$w_{1}$ | $10^{5}\,s^{-1}m^{-2}$ | $k_{t}$ |  $232.5\times 10^{3}$ N/m
$w_{2}$ | $0.5\,s^{3}m^{-2}$ | $b_{t}$ |  $0$ Ns/m
$w_{3}$ | $10^{-5}\,s^{-1}N^{-2}$ |  |

(a) (b) (c) (d) (e)

Figure 16: Vehicle suspension results

## References

  * Chen et al. [2020] S. Chen, H. Wang, M. Morari, V. M. Preciado, N. Matni,  Robust closed-loop model predictive control via system level synthesis,  in: 2020 59th IEEE Conference on Decision and Control (CDC), IEEE, 2020, pp. 2152–2159.
  * Foley et al. [2005] M. W. Foley, R. H. Julien, B. R. Copeland,  A comparison of pid controller tuning methods,  The Canadian Journal of Chemical Engineering 83 (2005) 712–722.
  * Leuer and Böcker [2014] M. Leuer, J. Böcker,  Real-time implementation of an online model predictive control for ipmsm using parallel computing on fpga,  in: 2014 International Power Electronics Conference (IPEC-Hiroshima 2014-ECCE ASIA), IEEE, 2014, pp. 346–350.
  * Bayat and Allison [2023] S. Bayat, J. T. Allison,  Ss-mpc: A user-friendly software based on single shooting optimization to solve model predictive control problems,  Software Impacts (2023) 100566.
  * Frison and Jørgensen [2015] G. Frison, J. B. Jørgensen,  Mpc related computational capabilities of armv7a processors,  in: 2015 European Control Conference (ECC), IEEE, 2015, pp. 3414–3421.
  * Bayat and Allison [2023] S. Bayat, J. T. Allison,  Control Co-Design with varying available information applied to vehicle suspensions,  in: International Design Engineering Technical Conferences and Computers and Information in Engineering Conference, American Society of Mechanical Engineers, 2023.
  * Lee and Markus [1967] E. B. Lee, L. Markus, Foundations of optimal control theory, Technical Report, Minnesota Univ Minneapolis Center For Control Sciences, 1967.
  * Giorgetti et al. [2006] N. Giorgetti, A. Bemporad, H. E. Tseng, D. Hrovat,  Hybrid model predictive control application towards optimal semi-active suspension,  International Journal of Control 79 (2006) 521–533.
  * Saerens et al. [2008] B. Saerens, M. Diehl, J. Swevers, E. Van den Bulck,  Model predictive control of automotive powertrains-first experimental results,  in: 2008 47th IEEE Conference on Decision and Control, IEEE, 2008, pp. 5692–5697.
  * Richter et al. [2010] S. Richter, S. Mariethoz, M. Morari,  High-speed online mpc based on a fast gradient method applied to power converter control,  in: Proceedings of the 2010 American Control Conference, IEEE, 2010, pp. 4737–4743.
  * Borrelli et al. [2006] F. Borrelli, A. Bemporad, M. Fodor, D. Hrovat,  An mpc/hybrid system approach to traction control,  IEEE Transactions on Control Systems Technology 14 (2006) 541–552.
  * Lucia et al. [2017] S. Lucia, A. Tătulea-Codrean, C. Schoppmeyer, S. Engell,  Rapid development of modular and sustainable nonlinear model predictive control solutions,  Control Engineering Practice 60 (2017) 51–62.
  * Houska et al. [2011] B. Houska, H. J. Ferreau, M. Diehl,  Acado toolkit—an open-source framework for automatic control and dynamic optimization,  Optimal Control Applications and Methods 32 (2011) 298–312.
  * Kühl et al. [2007] P. Kühl, J. Ferreau, J. Albersmeyer, C. Kirches, L. Wirsching, S. Sager, A. Potschka, G. Schulz, M. Diehl, D. B. Leineweber, et al.,  Muscod-ii users manual,  University of Heidelberg (2007).
  * Andersson et al. [2019] J. A. Andersson, J. Gillis, G. Horn, J. B. Rawlings, M. Diehl,  Casadi: a software framework for nonlinear optimization and optimal control,  Mathematical Programming Computation 11 (2019) 1–36.
  * Arendt and Veje [2019] K. Arendt, C. T. Veje, Mshoot: an open source framework for multiple shooting mpc in buildings, 2019.
  * Blum and Wetter [2019] D. Blum, M. Wetter, MPCPy: An open-source software platform for model predictive control in buildings, Technical Report, Lawrence Berkeley National Lab.(LBNL), Berkeley, CA (United States), 2019.
  * Jorissen et al. [2019] F. Jorissen, W. Boydens, L. Helsen,  Taco, an automated toolchain for model predictive control of building systems: implementation and verification,  Journal of building performance simulation 12 (2019) 180–192.
  * Bordons et al. [2020] C. Bordons, F. Garcia-Torres, M. A. Ridao, C. Bordons, F. Garcia-Torres, M. A. Ridao,  Model predictive control fundamentals,  Model Predictive Control of Microgrids (2020) 25–44.
  * Wang [2009] L. Wang, Model predictive control system design and implementation using MATLAB®, Springer Science & Business Media, 2009.
  * Allgöwer and Zheng [2012] F. Allgöwer, A. Zheng, Nonlinear model predictive control, volume 26, Birkhäuser, 2012.
  * Herber [2017] D. R. Herber, Advances in combined architecture, plant, and control design, Ph.D. thesis, University of Illinois at Urbana-Champaign, 2017.
  * Herber [2014] D. R. Herber,  Dynamic system design optimization of wave energy converters utilizing direct transcription (2014).
  * Rao [2009] A. V. Rao,  A survey of numerical methods for optimal control,  Advances in the Astronautical Sciences 135 (2009) 497–528.
  * Fahroo and Ross [2008] F. Fahroo, I. M. Ross,  Advances in pseudospectral methods for optimal control,  in: AIAA guidance, navigation and control conference and exhibit, 2008, p. 7309.
  * Bryson and Ho [2018] A. E. Bryson, Y.-C. Ho, Applied optimal control: optimization, estimation, and control, Routledge, 2018.
  * Herber and Allison [2019] D. R. Herber, J. T. Allison,  A problem class with combined architecture, plant, and control design applied to vehicle suspensions,  Journal of Mechanical Design 141 (2019).