# Modeling the amplitude and energy decay of a weakly damped harmonic
oscillator using the energy dissipation rate and a simple trick

Karlo Lelas<EMAIL_ADDRESS>Faculty of Textile Technology, University of
Zagreb, Prilaz baruna Filipovića 28a, 10000 Zagreb, Croatia  Robert Pezer
Department of Physical Metallurgy, Faculty of Metallurgy, University of
Zagreb, Aleja narodnih heroja 3, 44000 Sisak, Croatia<EMAIL_ADDRESS>
###### Abstract

We demonstrate how to derive the exponential decrease of amplitude and an
excellent approximation of the energy decay of a weakly damped harmonic
oscillator. This is achieved using a basic understanding of the undamped
harmonic oscillator and the connection between the damping force’s power and
the energy dissipation rate. The trick is to add the energy dissipation rates
corresponding to two specific pairs of initial conditions with the same
initial energy. In this way, we obtain a first-order differential equation
from which we quickly determine the time-dependent amplitude and the energies
corresponding to each pair of considered initial conditions. Comparing the
results of our model to the exact solutions and energies yielded an excellent
agreement. The physical concepts and mathematical tools we utilize are
familiar to first-year undergraduates.

##  I Introduction

The damped harmonic oscillator is a topic regularly covered in physics
textbooks for the first year of undergraduate studies, e.g., see [1, 2, 3].
Sometimes, only figures with a graphic representation of solutions of the
corresponding equation of motion are given [1], and occasionally, analytical
expressions of these solutions are given along with such figures [2, 3].
Still, none of these books show the derivation of these solutions. The reason
for this is that the equation of motion of the damped harmonic oscillator is a
second-order differential equation, the solution of which requires
mathematical tools that students at that level are not yet familiar with.
Thus, solutions are just given, and usually, only the exponential
approximation of the energy decay in the weak damping limit is discussed [2,
3].

On the other hand, in these same books [1, 2, 3], before getting to the part
with the damped harmonic oscillator, the simple harmonic oscillator (i.e.,
undamped harmonic oscillator) is treated in great detail. Students know the
solutions of an undamped harmonic oscillator and how to relate the amplitude
and phase of the undamped solution to the initial conditions. Another
important concept that students learn even before learning about undamped
harmonic motion is the concept of instantaneous power and its connection to
the rate of change of the system’s energy [2, 3].

In this paper, we show how to model the solutions and complete behavior of the
energy of a weakly damped harmonic oscillator, using only prior knowledge
about the undamped harmonic oscillator and the connection between the power of
the damping force and the energy dissipation rate. Of course, we must also
have some physical insight into the behavior of the system we are modeling.
Therefore, it would be best to start by providing students with figures of the
results of carefully conducted experiments on systems described by the damped
harmonic oscillator model, such as experiments with mechanical oscillatory
systems [4, 5] or with an oscillating RLC circuits [6], in a weakly damped
regime. Another option we choose in this paper is to show graphically the
exact solutions of the equation of motion of a weakly damped harmonic
oscillator as a substitute for the experimental results to gain insight into
its behavior and the approximations we can use.

This paper is organized into five sections. In section II, we present the
theory needed for modeling a weakly damped harmonic oscillator, and we state
the problem we are dealing with in this paper. In section III, we explain the
approximations we use in our model and derive the exponential decrease of the
amplitude and the corresponding expressions for the energies for the two pairs
of initial conditions we consider. In section IV, we compare the results of
our model with exact solutions and energies. In section V, we comment on a
simpler version of our model and derivation, which are more suitable for high
school students.

##  II Basic theory and the statement of the problem

As an example of a damped harmonic oscillator, we consider a block of mass $m$
that oscillates under the influence of the restoring force of an ideal spring
of stiffness $k$ and the damping force $F_{d}(t)=-bv(t)$, where $b$ is the
damping constant and $v(t)$ is the velocity of the block. The equation of
motion of this system is [2]

$ma(t)=-bv(t)-kx(t)\,,$ (1)

where $x(t)$ is the block’s displacement from the equilibrium position, and
$a(t)$ is its acceleration. Since $v(t)=dx(t)/dt$ and $a(t)=d^{2}x(t)/dt^{2}$,
(1) is a second-order linear differential equation. The methods for finding
its solutions are not given in physics textbooks for the first year of
undergraduate studies, e.g., see [1, 2, 3]. For any $b\geq 0$, the energy
(potential plus kinetic) of the block-spring system is given by

$E(t)=\frac{kx^{2}(t)}{2}+\frac{mv^{2}(t)}{2}\,.$ (2)

If we put $b=0$ in (1) we get the equation of the undamped harmonic oscillator
[2], with general solution $x(t)=A_{0}\cos(\omega_{0}t+\phi_{0})$, where
$\omega_{0}=\sqrt{k/m}$ is the angular frequency of the undamped system, while
$A_{0}$ and $\phi_{0}$ are constants (amplitude and phase) that are determined
from the initial conditions. The undamped system oscillates with conserved,
i.e., constant, energy [2]. For the damped systems, i.e., for $b>0$, the
energy is not conserved due to the power of the damping force
$P_{d}(t)=F_{d}(t)v(t)$, and the energy dissipation rate is given by [3]

$\frac{dE(t)}{dt}=P_{d}(t)=-bv^{2}(t)\,.$ (3)

Figure 1: (a) The red curve shows the displacement $\tilde{x}_{1}(t)$ of
weakly damped system (with some $m$, $k$ and $b$) for initial conditions
$\tilde{x}_{1}(0)=x_{0}$ and $\tilde{v}_{1}(0)=0$. For comparison, the black
curve shows the undamped solution with the same initial conditions, i.e.,
$x(t)=x_{0}\cos(\omega_{0}t)$. (b) The blue curve shows the displacement
$\tilde{x}_{2}(t)$ of the same weakly damped system for initial conditions
$\tilde{x}_{2}(0)=0$ and $\tilde{v}_{2}(0)=v_{0}=\omega_{0}x_{0}$, while the
black curve shows the undamped solution with these initial conditions, i.e.,
$x(t)=x_{0}\sin(\omega_{0}t)$, for comparison. In both figures, the time is in
$T_{0}=2\pi/\omega_{0}$ units, and the displacement is in $x_{0}$ units.

Let’s imagine that the instructor gives his students Fig. 1 with graphic
representations of the exact solutions of equation (1) for a weakly damped
system and an undamped system, for two pairs of initial conditions, which have
the same initial energy $E_{0}$. One pair of initial conditions has a purely
potential initial energy, while the other pair of initial conditions has a
purely kinetic initial energy. The analytical expressions of the damped
solutions used in Fig. 1 are given in section IV. Now comes the statement of
the problem. The instructor says to the students: ”Find the expression for the
energy decay of weakly damped system, using the insights provided by Fig. 1
and the energy dissipation rate (3). Furthermore, find what condition $b$ must
satisfy, compared to $m$ and $k$, for the system to be weakly damped.”

##  III Modeling of the amplitude decrease and the energy decay of weakly
damped harmonic oscillator

In Fig. 1(a) and (b), we see that the frequency and phase of
$\tilde{x}_{1}(t)$ and $\tilde{x}_{2}(t)$ remain unchanged compared to the
undamped solutions (at least as far as we can see from Fig. 1(a) and (b)),
while the amplitudes of $\tilde{x}_{1}(t)$ and $\tilde{x}_{2}(t)$ slowly
decrease with time in the same fashion. It is reasonable to describe the slow
decay of the amplitudes of $\tilde{x}_{1}(t)$ and $\tilde{x}_{2}(t)$ with the
same function, even though we can notice that the first few extremes of
$\tilde{x}_{2}(t)$ are slightly smaller in magnitude than the extremes of
$\tilde{x}_{1}(t)$ (for later times, the difference in the magnitude of
extremes is hardly visible), since this can be explained by the fact that
$\tilde{x}_{2}(t)$ reaches each extreme $T_{0}/4$ time later than
$\tilde{x}_{1}(t)$, so $\tilde{x}_{2}(t)$ is exposed to damping for more
extended time up to the moment of reaching the extremes compared to
$\tilde{x}_{1}(t)$. These insights suggest that the weakly damped
displacements shown in Fig. 1(a) and (b) can be modeled with

$x_{1}(t)=f(t)x_{0}\cos(\omega_{0}t)\,$ (4)

and

$x_{2}(t)=f(t)x_{0}\sin(\omega_{0}t)\,$ (5)

where function $f(t)$ describes the amplitude decrease in time. From (4) and
(5) we calculate the corresponding velocities

$v_{1}(t)=\frac{dx_{1}(t)}{dt}=-f(t)\omega_{0}x_{0}\sin(\omega_{0}t)+\dot{f}(t)x_{0}\cos(\omega_{0}t)$
(6)

and

$v_{2}(t)=\frac{dx_{2}(t)}{dt}=f(t)\omega_{0}x_{0}\cos(\omega_{0}t)+\dot{f}(t)x_{0}\sin(\omega_{0}t)\,,$
(7)

where we used notation $\dot{f}(t)=df(t)/dt$. Since $\tilde{x}_{1}(0)=x_{0}$
and $\tilde{v}_{2}(0)=\omega_{0}x_{0}$ we can easily see that $f(0)=1$ must
hold in order that our modeled solutions $x_{1}(t)$ and $x_{2}(t)$ satisfy
these initial conditions. Furthermore, our modeled solution $x_{1}(t)$ has a
initial velocity, i.e. $v_{1}(0)=\dot{f}(0)x_{0}$, while $\tilde{x}_{1}(t)$
has zero initial velocity. We deal with this issue in the following two
paragraphs.

It is clear from Fig. 1(a) and (b) that the amplitudes of $\tilde{x}_{1}(t)$
and $\tilde{x}_{2}(t)$ change very slowly in time in comparison to the
oscillating parts. Therefore, the rate of change of function $f(t)$ over time
has to be much slower than the rate of change of $\cos(\omega_{0}t)$ and
$\sin(\omega_{0}t)$, which is of the order $\omega_{0}$. We can conclude that
$|\dot{f}(t)|\ll\omega_{0}$ holds for weak damping for all time, i.e., for
$t\geq 0$ (we have put absolute value because we can expect the negative
derivative of $f(t)$ since the amplitude decreases with time).

Regarding the discrepancy of our modeled solution $x_{1}(t)$, i.e., (4), and
$\tilde{x}_{1}(t)$ in initial velocity, we can easily see that adding a phase
$\phi$ to cosine in $x_{1}(t)$ could resolve this issue since it would lead to
the initial velocity of the form

$v_{1}(0)=-f(0)\omega_{0}x_{0}\sin(\phi)+\dot{f}(0)x_{0}\cos(\phi)$ (8)

and, by setting $v_{1}(0)=0$, one easily gets that
$\phi=\arctan(\dot{f}(0)/\omega_{0}))$ resolves this issue, but, due to the
condition $|\dot{f}(0)|\ll\omega_{0}$, we can see that $\phi\ll 1$ holds. This
leads us to the conclusion that $\tilde{x}_{1}(t)$ has a phase, but which is
so small that it is unobservable in Fig. 1(a). Thus, adding this phase into
our model would be an unnecessary complication, i.e., (4) is a valid model of
$\tilde{x}_{1}(t)$.

By using our modeled displacements and velocities in (2), we get the
corresponding energies

$E_{1}(t)=E_{0}\left[f^{2}(t)-f(t)\left(\frac{\dot{f}(t)}{\omega_{0}}\right)\sin(2\omega_{0}t)+\left(\frac{\dot{f}(t)}{\omega_{0}}\right)^{2}\cos^{2}(\omega_{0}t)\right]$
(9)

and

$E_{2}(t)=E_{0}\left[f^{2}(t)+f(t)\left(\frac{\dot{f}(t)}{\omega_{0}}\right)\sin(2\omega_{0}t)+\left(\frac{\dot{f}(t)}{\omega_{0}}\right)^{2}\sin^{2}(\omega_{0}t)\right]\,,$
(10)

where $E_{0}=m\omega_{0}^{2}x_{0}^{2}/2$ is the initial energy, and we used
well know trigonometric identities
$\sin^{2}(\omega_{0}t)+\cos^{2}(\omega_{0}t)=1$ and
$\sin(2\omega_{0}t)=2\sin(\omega_{0}t)\cos(\omega_{0}t)$ in deriving (9) and
(10). We can see that our modeled solution $x_{2}(t)$ has the right initial
energy, i.e. $E_{2}(0)=E_{0}$, while $x_{1}(t)$ has initial energy
$E_{1}(0)=E_{0}\left(1+(\dot{f}(0)/\omega_{0})^{2}\right)$. Since
$|\dot{f}(t)|\ll\omega_{0}$ holds for weak damping, we can safely neglect the
terms of the order $(\dot{f}(t)/\omega_{0})^{2}$ in (9) and (10). Thus, to the
first order in $(\dot{f}(t)/\omega_{0})$, our modeled expressions for the
energies are

$E_{1}(t)=E_{0}\left[f^{2}(t)-f(t)\left(\frac{\dot{f}(t)}{\omega_{0}}\right)\sin(2\omega_{0}t)\right]$
(11)

and

$E_{2}(t)=E_{0}\left[f^{2}(t)+f(t)\left(\frac{\dot{f}(t)}{\omega_{0}}\right)\sin(2\omega_{0}t)\right]\,.$
(12)

Both (11) and (12) have the correct initial energy. Now we focus on
determining the $f(t)$.

If one would put, e.g., the energy (11) into the equation (3), i.e.,
$dE_{1}(t)/dt=-bv_{1}(t)^{2}$, it would lead to a very complicated
differential equation for $f(t)$, but instead of that we can use the following
trick. We add the dissipation rates corresponding to two pairs of initial
conditions, i.e.

$\frac{d}{dt}\left(E_{1}(t)+E_{2}(t)\right)=-b\left(v_{1}^{2}(t)+v_{2}^{2}(t)\right)\,.$
(13)

Since $E_{1}(t)+E_{2}(t)=2E_{0}f(t)^{2}$, and
$v_{1}^{2}(t)+v_{2}^{2}(t)=\omega_{0}^{2}x_{0}^{2}f^{2}(t)$ (in adding the
squared velocities we again neglected the terms of the order
$(\dot{f}(t)/\omega_{0})^{2}$) we get the differential equation

$\frac{d}{dt}f^{2}(t)=-\frac{b}{m}f^{2}(t)\,.$ (14)

Using substitutions $F(t)=f^{2}(t)$ and $\alpha=-b/m$, (14) becomes
$dF(t)/dt=\alpha F(t)$, and first year undergraduates know that the function
$F(t)=e^{\alpha t}$ satisfies this differential equation [3]. Thus, we obtain

$f(t)=\sqrt{F(t)}=e^{-\frac{b}{2m}t}$ (15)

as the slowly decaying amplitude of our modeled solutions (4) and (5). Using
(15), and its time derivative, in (11) and (12) we get

$E_{1}(t)=E_{0}e^{-\frac{b}{m}t}\left[1+\frac{b}{2m\omega_{0}}\sin(2\omega_{0}t)\right]$
(16)

and

$E_{2}(t)=E_{0}e^{-\frac{b}{m}t}\left[1-\frac{b}{2m\omega_{0}}\sin(2\omega_{0}t)\right]\,$
(17)

as final expressions for the energies in our model.

We can now easily find what condition $b$ must satisfy, compared to $m$ and
$k$, for the system to be weakly damped. As we stated before, for weak damping
condition $|\dot{f}(t)|\ll\omega_{0}$ holds. This condition holds for all $t$,
thus, if we consider $t=0$, for which $\dot{f}(0)=-b/(2m)$, we get

$\frac{b}{2m}\ll\omega_{0}$ (18)

as the condition of weak damping. Since $\omega_{0}=\sqrt{k/m}$, (18) can be
rewritten as $b\ll 2\sqrt{mk}$, which is a common form of the weak damping
condition in physics textbooks [3].

##  IV Comparison of our model with exact results

The exact solution of differential equation (1), for general initial
conditions, is given, e.g., in [6]. For a pair of initial conditions, which we
considered in Fig. 1(a) and (b), the exact solutions are

$\tilde{x}_{1}(t)=\frac{\omega_{0}x_{0}}{\omega}e^{-\frac{b}{2m}t}\cos\left(\omega
t-\arctan\left(\frac{b}{2m\omega}\right)\right)\,$ (19)

and

$\tilde{x}_{2}(t)=\frac{\omega_{0}x_{0}}{\omega}e^{-\frac{b}{2m}t}\sin\left(\omega
t\right)\,$ (20)

where $\omega=\sqrt{\omega_{0}^{2}-\left(\frac{b}{2m}\right)^{2}}$ is the
damped angular frequency. In Fig. 1(a) and (b) we took
$\frac{b}{2m}=0.02\omega_{0}$, which is in good agreement with weak damping
condition (18). In this case, (19) and (20) become

$\tilde{x}_{1}(t)=1.0002\,e^{-0.02\omega_{0}t}x_{0}\cos\left(0.9998\,\omega_{0}t-0.0020\right)\,$
(21)

and

$\tilde{x}_{1}(t)=1.0002\,e^{-0.02\omega_{0}t}x_{0}\sin\left(0.9998\,\omega_{0}t\right)\,,$
(22)

where we rounded numbers to four decimals, thus, even without showing a
graphical comparison, it is clear that our modeled solution (4) and (5), with
$f(t)=e^{-0.02\omega_{0}t}$, approximate excellently the exact solutions (21)
and (22). To further test the validity, but also to gain insight into the
limits of our model, in Fig. 2 we show the modeled solutions (4) and (5) for
$\frac{b}{2m}=0.1\omega_{0}$, and the corresponding exact solutions (19) and
(20). We see only a slight mismatch between the modeled and the exact
solutions. Thus, our model describes well the exact solutions (at least) for
$0<\frac{b}{2m}\leq 0.1\omega_{0}$.

Figure 2: Modeled solutions $x_{1}(t)$ (red curve) and $x_{2}(t)$ (blue curve)
given by (4) and (5) for $\frac{b}{2m}=0.1\omega_{0}$. The corresponding exact
solutions $\tilde{x}_{1}(t)$ and $\tilde{x}_{2}(t)$ given by (19) and (20) are
shown as dotted curves.

The exact expression for the energy of the damped harmonic oscillator, for
general initial conditions, is given, e.g., in [7]. For a pair of initial
conditions, which we considered in Fig. 1(a) and (b), the exact energies are

$\tilde{E}_{1,2}(t)=E_{0}e^{-\frac{b}{m}t}\left[\cos^{2}(\omega
t)\pm\frac{b}{2m\omega}\sin(2\omega
t)+\frac{\omega_{0}^{2}+\left(\frac{b}{2m}\right)^{2}}{\omega^{2}}\sin^{2}(\omega
t)\right]\,,$ (23)

where the upper sign in $\pm$ corresponds to $\tilde{E}_{1}(t)$, and lower
sign to $\tilde{E}_{2}(t)$. In Fig. 3 we show the modeled energies and the
exact energies for
$\frac{b}{2m}=\\{0.02\omega_{0},0.05\omega_{0},0.1\omega_{0}\\}$. We can see
excellent agreement between our modeled energies and the exact energies.

Figure 3: (a) Modeled energy $E_{1}(t)$, i.e. (16). (b) Modeled energy
$E_{2}(t)$, i.e. (17). Both figures show the behavior for
$\frac{b}{2m}=\\{0.02\omega_{0},0.05\omega_{0},0.1\omega_{0}\\}$. The
corresponding exact energies $\tilde{E}_{1}(t)$ and $\tilde{E}_{2}(t)$, given
by (23), are shown as dotted curves on both figures.

##  V A simpler version of our model

Exponential amplitude decay and purely exponential approximation of the energy
decay of weakly damped harmonic oscillator have recently been derived using
the energy dissipation rate averaged over one period of the corresponding
undamped system, i.e. $T_{0}$, with the approximation that the amplitude
remains constant over time intervals $\Delta t\sim T_{0}$ [8]. In our model,
this approximation corresponds to neglecting terms with $\dot{f}(t)$ in the
velocities (6) and (7), while the displacements (4) and (5) remain the same.
In this case, it is easy to show that the modeled energies for both pairs of
initial conditions are the same, i.e.

$E_{1}(t)=E_{2}(t)=\frac{m\omega_{0}^{2}x_{0}^{2}}{2}f^{2}(t)\,,$ (24)

and the corresponding dissipation rates are

$\frac{dE_{1}(t)}{dt}=-bf^{2}(t)\omega_{0}^{2}x_{0}^{2}\sin^{2}(\omega_{0}t)\,$
(25)

and

$\frac{dE_{2}(t)}{dt}=-bf^{2}(t)\omega_{0}^{2}x_{0}^{2}\cos^{2}(\omega_{0}t)\,.$
(26)

Now, it is even easier to see that the sum of (25) and (26) leads to equation
(14). Thus, we obtain $f(t)=e^{-\frac{b}{2m}t}$ again, but in this case, we
only get the exponential energy decay, i.e.,
$E_{1}(t)=E_{2}(t)=E_{0}e^{-\frac{b}{m}t}$, which is a good approximation of
the exact energies (23), with $\frac{b}{2m}\ll\omega_{0}$, only if one is
interested in behavior of the energy on time scales greater than $T_{0}$. This
simplified version of our derivation is suitable for high school students,
i.e., from a mathematical point of view, it is at a slightly lower level of
difficulty than the method presented in [8] since there is no need to
calculate the time integral of the energy dissipation rate if our trick is
used.

##  VI Conclusion

We showed how to obtain an excellent approximation of the solutions and the
energy of a weakly damped harmonic oscillator without solving the differential
equation of the second order, i.e., equation (1), and without any previous
insights into the analytical form of the solutions of that equation. The model
we introduced and the trick we used are suitable for first-year
undergraduates. In addition, we commented on a simplified version of our model
and derivation, which are more suitable for high school students.

##  VII Acknowledgments

This work was supported by the QuantiXLie Center of Excellence, a project co-
financed by the Croatian Government and European Union through the European
Regional Development Fund, the Competitiveness and Cohesion Operational
Programme (Grant No. KK.01.1.1.01.0004).

## References

  * [1] J.D. Cutnell and K.W. Johnson.  Physics.  John Wiley & Sons, 2009.
  * [2] David Halliday, Robert Resnick, and Jearl Walker.  Fundamentals of Physics.  John Wiley & Sons, 2013.
  * [3] Hugh D. Young and Roger A. Freedman.  University Physics with Modern Physics.  Pearson, 2020.
  * [4] Peter F. Hinrichsen.  Acceleration, Velocity, and Displacement for Magnetically Damped Oscillations.  The Physics Teacher, 57(4):250–253, 04 2019.
  * [5] Manuel I González and Alfredo Bol.  Controlled damping of a physical pendulum: experiments near critical conditions.  European Journal of Physics, 27(2):257, jan 2006.
  * [6] Karlo Lelas, Nikola Poljak, and Dario Jukić.  Damped harmonic oscillator revisited: The fastest route to equilibrium.  American Journal of Physics, 91(10):767–775, 10 2023.
  * [7] K. Lelas and I. Nakić.  Optimal damping of vibrating systems: Dependence on initial conditions.  Journal of Sound and Vibration, 576:118303, 2024.
  * [8] Stylianos-Vasileios Kontomaris and Anna Malamou.  An elementary proof of the amplitude’s exponential decrease in damped oscillations.  Physics Education, 59(2):025030, feb 2024.


# Flexible-step MPC for Switched Linear Systems with No Quadratic Common
Lyapunov Function

Annika Fürnsinn  Department of Mathematics and Statistics
Queen’s University, Kingston, ON K7L 3N6, Canada<EMAIL_ADDRESS>,
Christian Ebenbauer  Chair of Intelligent Control Systems
RWTH Aachen University, 52062 Aachen, Germany<EMAIL_ADDRESS>aachen.de and  Bahman Gharesifard  Department of Electrical and Computer
Engineering
University of California at Los Angeles, Los Angeles, CA 90095
<EMAIL_ADDRESS>
###### Abstract.

In this paper, we develop a systematic method for constructing a generalized
discrete-time control Lyapunov function for the flexible-step Model Predictive
Control (MPC) scheme, recently introduced in [3], when restricted to the class
of linear systems. Specifically, we show that a set of Linear Matrix
Inequalities (LMIs) can be used for this purpose, demonstrating its
tractability. The main consequence of this LMI formulation is that, when
combined with flexible-step MPC, we can effectively stabilize switched control
systems, for which no quadratic common Lyapunov function exists.

##  1\. Introduction

Model predictive control (MPC) is a widely utilized method to approximate the
solution of an infinite-horizon optimal control problem, where the underlying
dynamics can be of nonlinear or linear nature. In this paper, we focus on the
latter. The method itself consists of considering a sequence of finite-horizon
optimal control problems implemented in a receding horizon fashion. The
performance and stability properties of MPC schemes heavily rely on terminal
conditions [7, 8, 9, 10]; in fact, in the standard scheme, the optimal value
function is chosen as the Lyapunov function and stability is ensured by
enforcing a one-step descent of this optimal value function.

In our recent work [3], we proposed the flexible-step MPC scheme, which
notably allows for implementing a flexible number of control inputs in each
iteration in a computationally attractive manner. In order to provide
stability guarantees, we introduced the notion of generalized discrete-time
control Lyapunov functions (g-dclfs). As opposed to the strict descent
required of classical Lyapunov functions, this notion only requires a descent
of the averaged values of these functions over a maximal allowable
implementation window. We incorporate this average decrease condition as a
constraint in our MPC formulation. As demonstrated in our prior work [3, 4],
the flexible-step implementation of our scheme has many advantages; for
example, it can be used for stabilizing nonholonomic systems, where it is
known that standard MPC may lack asymptotic convergence [13].

A characteristic of Lyapunov theory is its reliance on a suitable Lyapunov
function. In spite of its powerful properties, the same is true for flexible-
step MPC; here, we need to find a suitable g-dclf. The purpose of the current
paper is to provide a systematic way of constructing such functions for the
class of linear control systems. In particular, we show that a simple
quadratic function can be selected as the g-dclf. In correspondence with the
related non-monotonic Lyapunov functions [2], we can cast the average decrease
constraint (adc) as an LMI, which can then be efficiently solved. More
importantly, we show that by choosing a large enough time step, this LMI
always admits a solution.

One advantage of this LMI characterization is that it can be utilized without
further knowledge about the cost function. This is due to the fact that in
flexible-step MPC, the two tasks of optimization and stabilization are
decoupled, i.e. the g-dclf can generally be a different function than the
optimal cost function. In standard MPC on the other hand, the tasks of
optimization and stabilization go hand in hand. The connection is particularly
apparent, when working with a linear system and a quadratic cost function. In
this case, the terminal cost consists of the solution to a Riccati or Lyapunov
equation. For non-quadratic cost functions, there is no straightforward
procedure for the design in standard MPC, however, as mentioned earlier, the
LMI approach can still be utilized for flexible-step MPC.

More importantly, the primary advantage of this LMI characterization is that
it allows us to use the flexible-step MPC scheme in the context of switched
linear systems. It is well-known that there are globally asymptotically stable
switched systems for which no classical quadratic common Lyapunov function
exists [11, 5]. In this paper, we provide an example of such a system, for
which, in contrast, a quadratic common g-dclf exists. In particular, we find
this function by coupling the LMI conditions for the individual linear
components. Our result demonstrates that similar to nonholonomic systems [3],
flexible-step MPC is also a useful tool for stabilizing switched linear
systems.

The remainder of the paper is organized as follows. We give a brief summary of
the flexible-step MPC scheme in Section 2. Section 3 contains our main result,
where we explain the translation of the adc constraint into an LMI and discuss
a numerical MPC example with a non-quadratic cost function. In Section 4, we
first apply the LMI approach and then flexible-step MPC to a linear switched
control system and present the numerical results. We discuss our future
directions in Section 5.

###  1.1. Notation

We denote the set of non-negative (positive) integers by $\mathbb{N}$
($\mathbb{N}_{>0}$), the set of (non-negative) reals by ($\mathbb{R}_{\geq
0}$), and the interior of a subset $S\subseteq\mathbb{R}^{n}$ by
$\operatorname{int}S$. The set of real matrices of the size $p\times q$ is
denoted by p,q. We make use of boldface when considering a sequence of finite
vectors, e.g. $\mathbf{u}=[u_{0},\dots,u_{N-1}]\in{}^{p,N}$, we refer to $j$th
component by $u_{j}$ and to the subsequence going from component $i$ to $j$ by
$\mathbf{u}_{[i:j]}$. Similarly,
$\mathbf{U}_{[0:N-1]}\subseteq\mathbb{R}^{p,N}$ denotes a set of (ordered)
sequences of vectors with components zero to $N-1$. When such a set depends on
the initial state $x$, it is expressed by $\mathbf{U}_{[0:N-1]}(x)$. As usual,
$\mathbf{u}^{*}$ denotes the solution to the optimal control problem solved in
the iteration of an MPC scheme, the symbol $\mathbf{u}^{*-}$ denotes the
optimal control sequence of the _previous_ iteration. For later use, we recall
that a function $V:\mathbb{R}^{n}\times\mathbb{R}^{p,q}\to\mathbb{R}$ is
called positive definite if $V(x,\mathbf{u})=0$ is equivalent to
$(x,\mathbf{u})=(0,\mathbf{0})$ and $V(x,\mathbf{u})>0$ for all
$(x,\mathbf{u})\in\mathbb{R}^{n}\times\mathbb{R}^{p,q}\setminus\\{(0,\mathbf{0})\\}$.
Furthermore, a function
$\alpha:\mathbb{R}^{n}\times\mathbb{R}^{p,q}\to\mathbb{R}$ is called radially
unbounded if $\|(x,\mathbf{u})\|\to\infty$ implies
$\alpha(x,\mathbf{u})\to\infty$, where $\|\cdot\|$ is the Euclidean norm. Note
that $x$ is a vector and $\mathbf{u}$ is a matrix here, so with
$\|(x,\mathbf{u})\|$ we implicitly refer to the norm of
$[x^{T},\operatorname{vec}(\mathbf{u})^{T}]$, where
$\operatorname{vec}(\mathbf{u})$ is the usual vectorization of a matrix into a
column vector.

In the setting of MPC, it is important to distinguish between predictions and
the actual implementations at a given time index $k\in\mathbb{N}$. In
particular, we use $x^{k},u^{k}$ to refer to _predictions_ , whereas we use
$x(k),u(k)$ to refer to the _actual_ states and _implemented_ inputs,
respectively.

##  2\. Flexible-step MPC

To begin with, let us consider nonlinear discrete-time control systems of the
form

$\displaystyle x^{k+1}$ $\displaystyle=f(x^{k},u^{k}),$ (1)

where $k\in\mathbb{N}$ denotes the time index, $x^{k}\in
X\subseteq\mathbb{R}^{n}$ is the state with the initial state $x^{0}\in X$ and
$u^{k}\in U\subseteq\mathbb{R}^{p}$ is an input. The state and input
constraints satisfy $0\in\operatorname{int}X$, $0\in\operatorname{int}U$ and
$f:X\times U\to\mathbb{R}^{n}$ is continuous with $f(0,0)=0$. The core of the
flexible-step MPC scheme consists of g-dclfs, which are defined as follows:

######  Definition 2.1 (Set of Feasible Controls).

For $x\in X$, we define a set of feasible controls as

$\displaystyle\mathbf{U}_{[0:N-1]}(x):=\\{\mathbf{u}_{[0:N-1]}\in\mathbb{R}^{p,N}$
$\displaystyle:\text{ with $x^{0}=x$},u_{j}\in U,x^{j+1}=f(x^{j},u_{j})\in X,$
(2) $\displaystyle j=0,\dots,N-1\\}.$

An infinite sequence of control inputs is called feasible when equation (2) is
fulfilled for all times $j\in\mathbb{N}$.

######  Definition 2.2 (g-dclf).

Consider the control system (1). Let $m\in\mathbb{N}_{>0}$ and
$q\in\mathbb{N}$. We call
$V:\mathbb{R}^{n}\times\mathbb{R}^{p,q}\to\mathbb{R}$ a generalized discrete-
time control Lyapunov function of order $m$ (g-dclf) for system (1) if $V$ is
continuous, positive definite and additionally:

If $q=0$:

  1. i)

there exists a continuous, radially unbounded and positive definite function
${\alpha:\mathbb{R}^{n}\to\mathbb{R}}$ such that for any $x^{0}\in$
$\mathbb{R}^{n}$ we have $V(x^{0})-\alpha(x^{0})\geq 0;$

  2. ii)

for any $x^{0}\in\mathbb{R}^{n}$ there exists
$\boldsymbol{\nu}_{[0:m-1]}\in\mathbf{U}_{[0:m-1]}(x^{0})$, which steers
$x^{0}$ to some $x^{m}$, such that

$\tfrac{1}{m}(\sigma_{m}V(x^{m})+\dots+\sigma_{1}V(x^{1}))-V(x^{0})\leq-\alpha(x^{0}),$
(3)

where $\sigma_{m},\ldots,\sigma_{1}\in{}_{\geq 0}$ and

$\tfrac{1}{m}(\sigma_{m}+\sigma_{m-1}+\ldots+\sigma_{1})-1\geq 0.$ (4)

If $q\neq 0$:

  1. i’)

there exists a continuous, radially unbounded and positive definite function
${\alpha:\mathbb{R}^{n}\times\mathbb{R}^{p,q}\to\mathbb{R}}$ such that for any
$(x^{0},\mathbf{u}_{[0:q-1]})\in$
${\mathbb{R}^{n}\times\mathbf{U}_{[0:q-1]}(x^{0})}$ we have

$\displaystyle
V(x^{0},\mathbf{u}_{[0:q-1]})-\alpha(x^{0},\mathbf{u}_{[0:q-1]})\geq 0;$ (5)

  2. ii’)

for any
$(x^{0},\mathbf{u}_{[0:q-1]})\in\mathbb{R}^{n}\times\mathbf{U}_{[0:q-1]}(x^{0})$
there exists $\boldsymbol{\nu}_{[0:q+m-1]}\in\mathbb{R}^{p,q+m}$ with
$\boldsymbol{\nu}_{[l:q+l-1]}\in\mathbf{U}_{[0:q-1]}(x^{l})$ for every
$l\in\\{0,1,\dots,m\\}$, which steers $x^{0}$ to some $x^{m}$, such that

$\displaystyle\tfrac{1}{m}(\sigma_{m}V(x^{m},\boldsymbol{\nu}_{[m:q+m-1]})+\dots+\sigma_{1}V(x^{1},\boldsymbol{\nu}_{[1:q]}))-V(x^{0},\mathbf{u}_{[0:q-1]})\leq-\alpha(x^{0},\mathbf{u}_{[0:q-1]}),$
(6)

where $\sigma_{m},\ldots,\sigma_{1}\in{}_{\geq 0}$ satisfy (4).

The interpretation of Definition 2.2 is that the sequence of Lyapunov function
values decreases on average every $m$ steps, and not at every step as in the
classical Lyapunov function. This is why we will refer to (6) (or (3) in the
case $q=0$) as the average decrease constraint (adc). This notion is key in
the flexible-step MPC scheme proposed in [3], which aims to solve the
following optimal control problem:

$\displaystyle\begin{aligned}
&\min\sum_{j=0}^{N_{p}-1}f_{0}(x^{j},u^{j})+\phi(x^{N_{p}})\\\ \ &\text{s.t.
}x^{j+1}=f(x^{j},u^{j}),\ j=0,\dots,N_{p}-1,\\\ &\hphantom{\text{s.t.
}}x^{0}=x,\\\ &\hphantom{\text{s.t. }}u^{j}\in U,\,x^{j}\in X,\,x^{N_{p}}\in
X^{N_{p}},\ j=0,\dots,N_{p}-1\\\ &\hphantom{\text{s.t.
}}[u^{l},\dots,u^{l+q-1}]\in\mathbf{U}_{[0:q-1]}(x^{l})\text{ for
}l=0,1,\dots,m\\\ &\hphantom{\text{s.t.
}}\tfrac{1}{m}\sum_{j=1}^{m}\sigma_{j}V(x^{j},[u^{j},\dots,u^{j-1+q}])-V(x^{0},\mathbf{u}^{*-}_{[0:q-1]})\leq-\alpha(x^{0},\mathbf{u}^{*-}_{[0:q-1]}),\end{aligned}$
(7)

where $N_{p}\in\mathbb{N},q\leq N_{p}$ with $q\in\mathbb{N},m\leq N_{p}$ with
$m\in\mathbb{N}_{>0},X^{N_{p}}\subseteq X\subseteq\mathbb{R}^{n}$ with
$0\in\operatorname{int}X^{N_{p}}$, $U\subseteq\mathbb{R}^{p}$ with
$0\in\operatorname{int}U$, $\sigma_{m},\ldots,\sigma_{1}\in{}_{\geq 0}$
satisfying (4), $\phi:\mathbb{R}^{n}\to\mathbb{R}$ is a positive semi-definite
function, the function $f_{0}:\mathbb{R}^{n}\times\mathbb{R}^{p}\to\mathbb{R}$
is positive definite and $V:\mathbb{R}^{n}\times\mathbb{R}^{p,q}\to\mathbb{R}$
is a g-dclf of order $m$. Note that condition (6) was added as a constraint to
the above optimal control problem. The flexible-step MPC scheme is given in
Algorithm 1.

Algorithm 1 Flexible-step MPC scheme

1: set $k=k_{0}$, measure the initial state $x(k_{0})$ and choose an arbitrary
$\mathbf{u}^{*-}_{[0:q-1]}\in\mathbf{U}_{[0:q-1]}(x(k_{0}))$

2: measure the current state $x(k)$ of (1)

3: solve problem (7) with $x=x(k)$ and obtain the optimal input
${\mathbf{u}^{*}_{[0:N-1]}}$, where $N=\max\\{q+m,N_{p}\\}$

4: choose an index $1\leq\ell_{\textup{decr}}\leq m$ for which
$V(x^{*\ell_{\textup{decr}}},\mathbf{u}^{*}_{[\ell_{\textup{decr}}:q+\ell_{\textup{decr}}-1]})-V(x,\mathbf{u}^{*-}_{[0:q-1]})\leq-\alpha(x,\mathbf{u}^{*-}_{[0:q-1]})$

5: implement
$\mathbf{u}^{*}_{[0:\ell_{\textup{decr}}-1]}=:[c^{k}_{mpc},\dots,c^{k+\ell_{decr}-1}_{mpc}]$
and redefine
$\mathbf{u}^{*-}_{[0:q-1]}:=\mathbf{u}^{*}_{[\ell_{\textup{decr}}:q+\ell_{\textup{decr}}-1]}$

6: increase $k:=k+\ell_{\textup{decr}}$ and go to 2

We refer the reader to our prior work in [3, 4] for further details. This
paper is concerned with finding a systematic way to design the g-dclf $V$, so
that it can be utilized for adc inside of (7). Similar to Lyapunov theory, in
general settings, this is quite a challenging task; for this reason, we
restrict ourselves to linear systems in what follows.

##  3\. g-dclfs for linear systems

From now on, we consider linear discrete-time control systems

$\displaystyle x^{k+1}=Ax^{k}+Bu^{k}$ (8)

with some initial state $x^{0}\in{}^{n}$,
$x^{k}\in{}^{n},\,u^{k}\in{}^{p},\,A\in{}^{n,n}$ and $B\in{}^{n,p}$, where
$(A,B)$ is stabilizable. Since we are now in the linear setting, let us focus
on the quadratic g-dclf

$V(x)=\|x\|^{2}.$

In this case, adc reads:

$\displaystyle\frac{\sigma_{m}\|x^{m}\|^{2}+\dots+\sigma_{1}\|x^{1}\|^{2}}{m}-\|x^{0}\|^{2}\leq-\alpha(x^{0}).$
(9)

Suppose we utilize state feedback $u=Kx$ with some matrix $K\in{}^{p,n}$.
Clearly, we have that

$x^{j}=(A+BK)^{j}x^{0},$

for all $j\in\\{1,\ldots m\\}$. Therefore,

$\mathbf{x}=\Phi x^{0},$

where $\mathbf{x}=((x^{1})^{T},\dots,(x^{m})^{T})^{T}$ and

$\displaystyle\Phi$ $\displaystyle=\begin{bmatrix}A+BK\\\ (A+BK)^{2}\\\
\vdots\\\ (A+BK)^{m}\end{bmatrix}\in{}^{nm,n}.$

By defining

$\displaystyle\Lambda:=\begin{bmatrix}\dfrac{\sigma_{1}}{m}I&&\\\ &\ddots&\\\
&&\dfrac{\sigma_{m}}{m}I\end{bmatrix}\in{}^{nm,nm},$

and also $\alpha(x^{0}):=\varepsilon\|x^{0}\|^{2}$, where $1>\varepsilon>0$,
we can rewrite (9) as a matrix inequality:

$\displaystyle\mathbf{x}^{T}\Lambda\mathbf{x}-(x^{0})^{T}x^{0}$
$\displaystyle=(\Phi x^{0})^{T}\Lambda\Phi x^{0}-(x^{0})^{T}x^{0}$
$\displaystyle\leq-\alpha(x^{0})=-\varepsilon\|x^{0}\|^{2}.$ By rearranging
further, we obtain $\displaystyle-(x^{0})^{T}\Phi^{T}\Lambda\Phi
x^{0}+(x^{0})^{T}Ix^{0}$
$\displaystyle\geq\varepsilon(x^{0})^{T}x^{0}\quad\forall x^{0}\in{}^{n}$
$\displaystyle\Rightarrow-\Phi^{T}\Lambda\Phi+I$
$\displaystyle\succeq\varepsilon I,$

which is a Linear Matrix Inequality (LMI) with the decision variable
$\Lambda$. If we pose no constraint on this LMI, the solution will just be the
zero matrix, as $I\succeq\varepsilon I$ for small $\varepsilon$. We do,
however, have two constraints on the weights: First, they need to be non-
negative; second, they need to satisfy (4). Altogether, our LMI consists of
three parts:

$\displaystyle-\Phi^{T}\Lambda\Phi+I-\varepsilon I$ $\displaystyle\succeq 0$
(10) $\displaystyle\Lambda$ $\displaystyle\succeq 0$ (11)
$\displaystyle\tfrac{1}{m}(\sigma_{1}+\dots+\sigma_{m})$ $\displaystyle\geq
1\text{ or }\operatorname{trace}(\Lambda)\geq n.$ (12)

This leads us to the first result of this paper.

######  Theorem 3.1.

Consider system (8) and a stabilizing matrix $K\in\mathbb{R}^{p,n}$, i.e. the
absolute value of each eigenvalue of $A+BK$ is strictly less than one. Then
there exists a sufficiently large $m\in\mathbb{N}_{>0}$ such that
$\sigma_{1},\dots,\sigma_{m}\in\mathbb{R}_{\geq 0}$ satisfy the LMI described
by (10), (11), (12) with $0<\varepsilon<1$.

###### Proof.

Suppose $K\in\mathbb{R}^{p,n}$ satisfies $\rho(A+BK)<1$. Our goal is to find a
solution to the LMI (10), (11), (12). Let us first consider the setting where
$\sigma_{1},\dots,\sigma_{m-1}=0$ and $\sigma_{m}=m$. With this choice, (11)
and (12) are automatically satisfied. We show that we can find
$m\in\mathbb{N}_{>0}$ such that the average descent property described by (10)
is satisfied. With the above choices, (10) reads

$\displaystyle-((A+BK)^{m})^{T}(A+BK)^{m}+(1-\varepsilon)I\succeq 0.$ (13)

Note that

$\displaystyle\,\|((A+BK)^{m})^{T}(A+BK)^{m}\|$
$\displaystyle\leq\|((A+BK)^{m})^{T}\|\|(A+BK)^{m}\|$
$\displaystyle=\|(A+BK)^{m}\|\|(A+BK)^{m}\|.$

Since $\rho(A+BK)<1$, it is a well-known fact [12] that
$\lim_{m\to\infty}\|(A+BK)^{m}\|=0$; therefore,

$\displaystyle\lim_{m\to\infty}((A+BK)^{m})^{T}(A+BK)^{m}$ $\displaystyle=0.$

Let us denote the eigenvalues of ${((A+BK)^{m})^{T}(A+BK)^{m}}$ by
$\\{\lambda_{m,1},\dots,\lambda_{m,n}\\}$. This in turn implies that each of
the eigenvalues $\lambda_{m,i}$ converges to the eigenvalues of the zero
matrix, which are all zero of course, as $m$ goes to infinity. Let us now
choose $\tilde{\varepsilon}$ with $0<\tilde{\varepsilon}\leq 1-\varepsilon$.
Then there exists $\tilde{m}\in\mathbb{N}_{>0}$ such that for all eigenvalues
$\lambda_{\tilde{m},i}$ we have $|\lambda_{\tilde{m},i}|<\tilde{\varepsilon}$.
Since the eigenvalues of

$-((A+BK)^{\tilde{m}})^{T}(A+BK)^{\tilde{m}}+(1-\varepsilon)I$ (14)

are equal to $(1-\varepsilon)-\lambda_{\tilde{m},i}$, we have that

$\displaystyle(1-\varepsilon)-\lambda_{\tilde{m},i}$
$\displaystyle\geq(1-\varepsilon)-|\lambda_{\tilde{m},i}|$
$\displaystyle>(1-\varepsilon)-\tilde{\varepsilon}\geq 0.$

With all the eigenvalues of the matrix in (14) being non-negative, we conclude
that the matrix itself is positive (semi) definite. In other words, we have
found $\tilde{m}\in\mathbb{N}_{>0}$ such that (13) is satisfied. In the more
general case, where all the weights $\sigma_{i}$ are positive, we have
$\tilde{m}-1$ more matrices to consider in (10). Clearly, we can find small
enough weights $\sigma_{1},\dots,\sigma_{\tilde{m}-1}$ such that
$-\Phi^{T}\Lambda\Phi$ is still dominated by
$\frac{\sigma_{\tilde{m}}}{\tilde{m}}(-((A+BK)^{\tilde{m}})^{T}(A+BK)^{\tilde{m}})$
and hence, stays positive definite. This finishes the proof. ∎

###  3.1. Numerical Example

To illustrate our results, we will use the introduced LMI approach for an
example. In particular, we consider the following problem:

######  Problem 3.2.

$\displaystyle\min$
$\displaystyle\sum_{j=0}^{N_{p}-1}\|x^{j}\|_{1}+5\|u^{j}\|^{2}$
$\displaystyle\mathrm{s.t.}$
$\displaystyle\,x^{j+1}=\begin{bmatrix}2.13&1&1\\\ 0&1&0.3\\\
0&0&0.5\end{bmatrix}x^{j}+\begin{bmatrix}0\\\ 0\\\ 1\end{bmatrix}u^{j}$
$\displaystyle\,x^{0}=x$
$\displaystyle\,\frac{1}{N_{p}}(\sigma_{N_{p}}\|x^{N_{p}}\|^{2}+\dots+\sigma_{1}\|x^{1}\|^{2})-\|x^{0}\|^{2}\leq-\varepsilon\|x^{0}\|^{2}.$

The initial state is $[4\,\,12\,\,15]^{T}$ and the prediction horizon is
chosen to be $N_{p}=10$. Note that this problem has a non-quadratic running
cost, which means that it is not directly clear which terminal cost should be
utilized in the standard MPC. To apply our methodology, we first select a
stabilizing feedback $K$ for the control system in Problem 3.2

$K=\begin{bmatrix}-3.5507&-2.6749&-2.4633\end{bmatrix}.$

Next, we find a sufficiently large $m\in\mathbb{N}_{>0}$ such that the LMI
(10), (11), (12) admits a solution; note that the existence of such an $m$ is
guaranteed by Theorem 3.1. Using numerical evaluations, one can observe that
selecting $m\geq 6$ suffices. We have chosen $m=10$ and have solved the LMI
(10), (11), (12) to obtain:

$\displaystyle\sigma_{1}=0.0055,\sigma_{2}=0.0524,\sigma_{3}=0.0660,\sigma_{4}=0.0655,\sigma_{5}=0.0762,$
$\displaystyle\sigma_{6}=0.0952,\sigma_{7}=0.1201,\sigma_{8}=0.1479,\sigma_{9}=0.1745,\sigma_{10}=0.1967.$

With the weights for adc at hand, we solve Problem 3.2 with the flexible-step
MPC scheme given in Algorithm 1 by using fmincon from MATLAB. In step four of
Algorithm 1, we find the index $\ell_{decr}\in\\{1,\dots,N_{p}\\}$ with the
_maximal_ descent for the g-dclf $V(x)=\|x\|^{2}$. We then choose
$\varepsilon=10^{-10}$ for the right-hand side of adc. Figure 1 shows the
state trajectories according to the solution of Problem 3.2. After a transient
phase, the state is successfully stabilized to the origin. A summary of the
implemented steps at each optimization instance of Problem 3.2 is given in
Figure 2, where notably, a flexible number of steps is implemented in each
optimization instance.

Figure 1. State trajectories according to the solution of Problem 3.2 Figure
2. In each optimization instance, Problem 3.2 is solved. The solution yields
the number of implemented steps, depicted on the vertical axis. Note that the
sum of implemented steps of the optimization instances corresponds to the time
step of the implementation.

##  4\. Quadratic Common g-dclfs

We now demonstrate how the introduced LMI approach can be beneficial in the
context of stabilization of switched linear systems. Let us start by recalling
some terminology [11]. We consider a family of nonlinear control systems
described by

$\displaystyle x^{k+1}=f_{\uptheta}(x^{k},u^{k}),\quad\uptheta\in\Uptheta,$
(15)

where $\Uptheta$ is some index set and
$f_{\uptheta}:\mathbb{R}^{n}\times\mathbb{R}^{p}\to\mathbb{R}^{n}$ is
continuous for each $\uptheta\in\Uptheta$. With abuse of notation, we denote
the switching signal by $\uptheta:\mathbb{N}\to\Uptheta$, which specifies, at
each time instant $k\in\mathbb{N}$, the index $\uptheta(k)\in\Uptheta$ of the
active subsystem, i.e. the system from the family (15) that is currently being
followed. Hence, a switched system with time-dependent switching can be
described by the equation

$x^{k+1}=f_{\uptheta(k)}(x^{k},u^{k}).$ (16)

In what follows, we discuss the asymptotic stability of the switched system
(16) under arbitrary switching. Let us precisely define what that entails.

######  Definition 4.1.

A switched system (16) is uniformly asymptotically stable if for all
$\varepsilon>0$ there exists a $\delta>0$ such that for all switching signals
$\uptheta\in\Uptheta$ the trajectories of (16) with $0<\|x^{0}\|<\delta$
satisfy $\|x^{k}\|<\varepsilon$.

We assume that the individual subsystems have the origin as the common
equilibrium point: $f_{\uptheta}(0,0)=0$ for all $\uptheta\in\Uptheta$. A
useful tool when investigating the stability of a _family_ of control systems
is the notion of _common_ Lyapunov functions, which we define next.

######  Definition 4.2.

A Lyapunov function $V:\mathbb{R}^{n}\to\mathbb{R}_{\geq 0}$ is called common
with respect to (15) if it is a continuous, radially unbounded, positive
definite function and satisfies the Lyapunov condition, i.e. there exists a
continuous, positive definite function $W:\mathbb{R}^{n}\to\mathbb{R}_{\geq
0}$ such that for all $x\in\mathbb{R}^{n}$ and all $\uptheta\in\Uptheta$ there
exists $u\in\mathbb{R}^{p}$ with

$V(f_{\uptheta}(x,u))-V(x)\leq-W(x).$

In other words, the single Lyapunov function $V$ decreases along trajectories
of all subsystems of (15). For simplicity, we assume in what follows that the
index set $\Uptheta$ is finite. Then the existence of a common Lyapunov
function abiding the slightly weaker condition

$V(f_{\uptheta}(x,u))-V(x)<0$ (17)

for all $x\in\mathbb{R}^{n}\setminus\\{0\\}$ and all $\uptheta\in\Uptheta$
implies global uniform (w.r.t. $\uptheta$) asymptotic stability, compare with
the continuous time result [11, Theorem 2.1].111 More generally, condition
(17) is sufficient for global uniform asymptotic stability under a compactness
assumption on the index set $\Uptheta$ and a continuity assumption on
$f_{\uptheta}$ with respect to $\uptheta$. Both of these assumptions trivially
hold when $\Uptheta$ is finite [11].

Similar to the common Lyapunov function, a common g-dclf for (15) can be
defined. Essentially, it must hold that for all $\uptheta\in\Uptheta$ there
exists a feasible control such that adc is satisfied.

In the context of linear systems, a natural question is whether it is
sufficient to work with _quadratic_ common Lyapunov functions, i.e. Lyapunov
functions $V:\mathbb{R}^{n}\to\mathbb{R}_{\geq 0}$ in the form of
$V(x)=x^{T}Px$, where $P\in\mathbb{R}^{n,n}$ is positive definite. The next
main result addresses precisely this question, with the first part being
classical [11].

######  Theorem 4.3.

The switched system

$x^{k+1}=x^{k}+\uptheta\,h\begin{bmatrix}0&1\\\
-(1+\alpha)&0\end{bmatrix}x^{k}+h\begin{bmatrix}0\\\ 1\end{bmatrix}u^{k},$
(18)

where $\uptheta\in\\{-1,1\\},h>0$ and $\alpha=0.5$, does not admit a quadratic
common control Lyapunov function. It does, however, admit a quadratic common
g-dclf for $h=0.1$.

We prove this result through two lemmata.

######  Lemma 4.4.

The switched system (18) does not admit a quadratic common control Lyapunov
function.

###### Proof.

The linear system (18) can be written more compactly as

$\displaystyle x^{k+1}$
$\displaystyle=A_{\uptheta}x^{k}+Bu^{k},\quad\uptheta\in\\{-1,1\\}$ (19) with
$\displaystyle A_{1}$ $\displaystyle=\begin{bmatrix}1&h\\\
-0.5h&1\end{bmatrix},\,A_{-1}=\begin{bmatrix}1&-h\\\
0.5h&1\end{bmatrix},\,B=h\begin{bmatrix}0\\\ 1\end{bmatrix}.$

For later use, we also introduce the matrices

$\displaystyle\mathcal{A}_{1}$ $\displaystyle=\begin{bmatrix}0&1\\\
-0.5&0\end{bmatrix},\,\mathcal{A}_{-1}=\begin{bmatrix}0&-1\\\
0.5&0\end{bmatrix},\,\mathcal{B}=\begin{bmatrix}0\\\ 1\end{bmatrix}.$

By way of contradiction, let us assume that a quadratic common control
Lyapunov function $V(x)=x^{T}Px$, with a positive definite matrix
$P\in\mathbb{R}^{n,n}$, exists. Without loss of generality, the matrix $P$ has
the form

$\displaystyle P=\begin{bmatrix}1&q\\\ q&r\end{bmatrix}.$

Since $\mathcal{B}=\begin{bmatrix}0&1\end{bmatrix}^{T}$, we have that

$\mathcal{B}^{T}Py=0,\quad\quad\mathrm{for}\
y=\begin{bmatrix}-r&q\end{bmatrix}^{T}.$

The function $V$ needs to satisfy the descent property (17)

$\displaystyle(A_{\uptheta}x+Bu)^{T}P(A_{\uptheta}x+Bu)-x^{T}Px$
$\displaystyle<0$

with one common $P$ for both $\uptheta=-1$ and $\uptheta=1$. Expanding the
above inequality yields

$\displaystyle
x^{T}A_{\uptheta}^{T}PA_{\uptheta}x+2uB^{T}PA_{\uptheta}x+uB^{T}PBu-x^{T}Px$
$\displaystyle<0.$

We rewrite the above inequality by using
$A_{\uptheta}=I+h\mathcal{A_{\uptheta}}$ and $B=h\mathcal{B}$

$\displaystyle
x^{T}(I+h\mathcal{A_{\uptheta}})^{T}P(I+h\mathcal{A_{\uptheta}})+2hu\mathcal{B}^{T}P(I+h\mathcal{A_{\uptheta}})x+h^{2}u\mathcal{B}^{T}P\mathcal{B}u-x^{T}Px$
$\displaystyle<0$ $\displaystyle\Rightarrow
x^{T}(P+h\mathcal{A_{\uptheta}}^{T}P+hP\mathcal{A_{\uptheta}}+h^{2}\mathcal{A_{\uptheta}}^{T}P\mathcal{A_{\uptheta}}-P)x+2hu\mathcal{B}^{T}P(I+h\mathcal{A_{\uptheta}})x+h^{2}u\mathcal{B}^{T}P\mathcal{B}u$
$\displaystyle<0$ $\displaystyle\Rightarrow
hx^{T}(\mathcal{A_{\uptheta}}^{T}P+P\mathcal{A_{\uptheta}})x+h^{2}x^{T}\mathcal{A_{\uptheta}}^{T}P\mathcal{A_{\uptheta}}x+2hu\mathcal{B}^{T}P(I+h\mathcal{A}_{\uptheta})x+h^{2}u\mathcal{B}^{T}P\mathcal{B}u$
$\displaystyle<0.$

By dividing by $h>0$, this is equivalent to

$\displaystyle
x^{T}(\mathcal{A}_{\uptheta}^{T}P+P\mathcal{A}_{\uptheta})x+hx^{T}\mathcal{A}_{\uptheta}^{T}P\mathcal{A}_{\uptheta}x+2u\mathcal{B}^{T}P(I+h\mathcal{A}_{\uptheta})x+hu\mathcal{B}^{T}P\mathcal{B}u<0.$

Since $P$ is positive definite, we have
$hu^{T}\mathcal{B}^{T}P\mathcal{B}u\geq 0$ and
$hx^{T}\mathcal{A}^{T}_{\uptheta}P\mathcal{A}_{\uptheta}x\geq 0$. Thus, in
order for the above expression to be negative, we must have that

$\displaystyle
x^{T}(\mathcal{A}_{\uptheta}^{T}P+P\mathcal{A}_{\uptheta})x+2u\mathcal{B}^{T}P(I+h\mathcal{A}_{\uptheta})x<0.$

To be precise, the following two inequalities need to hold

$\displaystyle
x^{T}(\mathcal{A}_{1}^{T}P+P\mathcal{A}_{1})x+2u\mathcal{B}^{T}P(I+h\mathcal{A}_{1})x$
$\displaystyle<0$ $\displaystyle
x^{T}(\mathcal{A}_{-1}^{T}P+P\mathcal{A}_{-1})x+2u\mathcal{B}^{T}P(I+h\mathcal{A}_{-1})x$
$\displaystyle<0.$

We will now evaluate these inequalities for two different $x$ values. We
evaluate the first equation at $x_{1}=A_{1}^{-1}y=(I+h\mathcal{A}_{1})^{-1}y$
and the second equation at $x_{-1}=A_{-1}^{-1}y=(I+h\mathcal{A}_{-1})^{-1}y$.
This yields

$\displaystyle x_{1}^{T}(\mathcal{A}_{1}^{T}P+P\mathcal{A}_{1})x_{1}+0$
$\displaystyle<0$ $\displaystyle
x_{-1}^{T}(\mathcal{A}_{-1}^{T}P+P\mathcal{A}_{-1})x_{-1}+0$ $\displaystyle<0$

because $\mathcal{B}^{T}Py=0$. Since

$\displaystyle A_{1}^{-1}$ $\displaystyle=\frac{1}{\det A_{1}}A_{-1},\quad
A_{-1}^{-1}=\frac{1}{\det A_{1}}A_{1}$ $\displaystyle\mathcal{A}_{-1}A_{1}$
$\displaystyle=0.5hI+\mathcal{A}_{-1},\quad\mathcal{A}_{1}A_{-1}=0.5hI+\mathcal{A}_{1},$

we obtain for the first inequality

$\displaystyle x_{1}^{T}(\mathcal{A}_{1}^{T}P+P\mathcal{A}_{1})x_{1}$
$\displaystyle<0$ $\displaystyle\frac{1}{(\det
A_{1})^{2}}y^{T}(A_{-1}^{T}\mathcal{A}_{1}^{T}P+P\mathcal{A}_{1}A_{-1})y$
$\displaystyle<0$
$\displaystyle\Leftrightarrow\,y^{T}((0.5hI+\mathcal{A}_{1})^{T}P+P(0.5hI+\mathcal{A}_{1}))y$
$\displaystyle<0$ $\displaystyle
y^{T}(\mathcal{A}_{1}^{T}P+P\mathcal{A}_{1})y+hy^{T}Py$ $\displaystyle<0.$

For the second inequality we obtain

$\displaystyle x_{-1}^{T}(\mathcal{A}_{-1}^{T}P+P\mathcal{A}_{-1})x_{-1}$
$\displaystyle<0$ $\displaystyle\frac{1}{(\det
A_{1})^{2}}y^{T}(A_{1}^{T}\mathcal{A}_{-1}^{T}P+P\mathcal{A}_{-1}A_{1})y$
$\displaystyle<0$
$\displaystyle\Leftrightarrow\,y^{T}((0.5hI+\mathcal{A}_{-1})^{T}P+P(0.5hI+\mathcal{A}_{-1}))y$
$\displaystyle<0$
$\displaystyle\,y^{T}(\mathcal{A}_{-1}^{T}P+P\mathcal{A}_{-1})y+hy^{T}Py$
$\displaystyle<0.$

Again, since $P$ is positive definite, we have $hy^{T}Py\geq 0$. Hence, in
order for the above inequalities to be negative, it must hold

$\displaystyle y^{T}(\mathcal{A}_{1}^{T}P+P\mathcal{A}_{1})y$
$\displaystyle<0$ $\displaystyle
y^{T}(\mathcal{A}_{-1}^{T}P+P\mathcal{A}_{-1})y$ $\displaystyle<0.$

Finally, by using the identity $\mathcal{A}_{-1}=-\mathcal{A}_{1}$, we see
that the above inequalities become

$\displaystyle y^{T}(\mathcal{A}_{1}^{T}P+P\mathcal{A}_{1})y$
$\displaystyle<0$ $\displaystyle
y^{T}(\mathcal{A}_{1}^{T}P+P\mathcal{A}_{1})y$ $\displaystyle>0,$

which is a contradictory statement. To summarize, the existence of a quadratic
common control Lyapunov function has led us to a contradiction, which means
that there cannot exist such a Lyapunov function. ∎

######  Lemma 4.5.

The switched system (18) admits a quadratic common g-dclf for $h=0.1$.

###### Proof.

It can be readily verified that for $m=10$ the weights222The weights can be
found through MATLAB and the sedumi package by solving the coupled LMI
$\displaystyle-\Phi_{1}^{T}\Lambda\Phi_{1}+I-\varepsilon I$
$\displaystyle\succeq 0$
$\displaystyle-\Phi_{-1}^{T}\Lambda\Phi_{-1}+I-\varepsilon I$
$\displaystyle\succeq 0$ $\displaystyle\Lambda$ $\displaystyle\succeq 0$
$\displaystyle\operatorname{trace}(\Lambda)$ $\displaystyle\geq n=2.$

$\displaystyle\sigma_{1}=0.0644,\sigma_{2}=0.0570,\sigma_{3}=0.0589,\sigma_{4}=0.0655,\sigma_{5}=0.0775,$
$\displaystyle\sigma_{6}=0.0959,\sigma_{7}=0.1227,\sigma_{8}=0.1646,\sigma_{9}=0.2488,\sigma_{10}=0.5447.$
(20)

satisfy (10), (11), (12) simultaneously for system (18) with $\uptheta=1$ and
system (18) with $\uptheta=-1$. This can be seen by using the stabilizing
matrices

$\displaystyle K_{1}=[-5.4017\,\,\,-7.0985],\,\,K_{-1}=[5.4017\,\,\,-7.0985]$

for (18) with $\uptheta=1$ and $\uptheta=-1$ and building $\Phi_{1}$ and
$\Phi_{-1}$ consisting of powers of $(A_{1}+BK_{1})$ and $(A_{-1}+BK_{-1})$,
respectively. Since the above weights satisfy (10), (11), (12) for the
switched system, i.e. system (18) with $\uptheta=1$ and $\uptheta=-1$, it
admits a common g-dclf, namely $V(x)=\|x\|^{2}$. ∎

We now illustrate how the quadratic common g-dclf facilitates stabilization by
solving the optimal control Problem 4.6 within the flexible-step MPC scheme
(Algorithm 1). The weights for adc are given by (4), the initial state is
$\begin{bmatrix}4&5\end{bmatrix}^{T}$, the prediction horizon is $N_{p}=10$
and we utilize the solver fmincon from MATLAB.

######  Problem 4.6.

$\displaystyle\min\sum_{j=0}^{N_{p}-1}\|x^{j}\|^{2}+5\|u^{j}\|^{2}$
$\displaystyle\,\mathrm{s.t.}\,x^{j+1}=x^{j}+h\begin{bmatrix}0\\\
1\end{bmatrix}u^{j}+\left\\{\begin{array}[]{l}h\begin{bmatrix}0&1\\\
-0.5&0\end{bmatrix}x^{j},\,\text{if }k+j\ \mathrm{even}\\\
-h\begin{bmatrix}0&1\\\
-0.5&0\end{bmatrix}x^{j}\hfill\mathrm{otherwise}\end{array}\right.$
$\displaystyle\hskip 22.76228pt\,x^{0}=x$ $\displaystyle\hskip
22.76228pt\,\frac{1}{N_{p}}(\sigma_{N_{p}}\|x^{N_{p}}\|^{2}+\dots+\sigma_{1}\|x^{1}\|^{2})-\|x^{0}\|^{2}\leq-\varepsilon\|x^{0}\|^{2}$

Here, $k$ is the time step of the closed-loop system. In step four of
Algorithm 1, we find the _first_ index $\ell_{decr}\in\\{1,\dots,N_{p}\\}$
which achieves a descent for the g-dclf $V(x)=\|x\|^{2}$. To achieve a
descent, we choose $\varepsilon=10^{-5}$ for the right-hand side of adc. The
state trajectories of the closed-loop system are displayed in Figure 3. It is
evident, that after a transient phase, the closed-loop system was successfully
stabilized to the origin. In Figure 4, the number of implemented steps in each
optimization instance is displayed. We see that in multiple optimization
instances, the solver makes use of more than one step of the computed control
sequence.

It is worth pointing out that our scheme works for any arbitrary switching
pattern (which must be known for the next $m$ time steps so that our scheme
can be implemented).

Figure 3. State trajectories according to the solution of Problem 4.6 Figure
4. In each optimization instance, Problem 4.6 is solved. The solution yields
the number of implemented steps, depicted on the vertical axis. Note that the
sum of implemented steps of the optimization instances corresponds to the time
step of the implementation.  Figure 5. State trajectories according to
standard MPC with the terminal cost
$\phi(x^{N_{p}})=480(x^{N_{p}})^{T}x^{N_{p}}$

We compare the results from flexible-step MPC with stan- dard MPC. In standard
MPC, the cost

$\displaystyle\sum_{j=0}^{N_{p}-1}\|x^{j}\|^{2}+5\|u^{j}\|^{2}+\phi(x^{N_{p}})$

is optimized in each iteration. In other words, the stabilizing ingredient of
a quadratic terminal cost $\phi:\mathbb{R}^{n}\to\mathbb{R}_{\geq 0}$ is used.
As it is unclear how to choose the terminal cost systematically, we have tried
three different terminal costs, i.e.
$\phi(x^{N_{p}})=(x^{N_{p}})^{T}P_{1}x^{N_{p}},\phi(x^{N_{p}})=(x^{N_{p}})^{T}P_{-1}x^{N_{p}}$
and $\phi(x^{N_{p}})=480(x^{N_{p}})^{T}x^{N_{p}}$, where $P_{1}$ and $P_{-1}$
are the solutions of the Riccati equation for the $\uptheta=1$ and
$\uptheta=-1$ dynamics (18), respectively. As expected, standard MPC fails to
stabilize the state with all of these terminal costs. This is depicted for one
of the terminal costs, e.g. $\phi(x^{N_{p}})=480(x^{N_{p}})^{T}x^{N_{p}}$, in
Figure 5, where the state trajectories are shown according to standard MPC.

##  5\. Conclusion

In this paper, we showed that the average decrease constraint of the recently
introduced framework of flexible-step MPC can be verified using an LMI in the
case of linear systems. As a key consequence, we were able to find a quadratic
common _generalized_ Lyapunov function for a switched linear control system,
where the existence of quadratic common Lyapunov functions is not guaranteed.
Using this, we showcased an exemplary optimal control problem with a switched
system, where flexible-step MPC can overcome challenges faced in standard MPC.
In the future, we will explore applications to robotics and power systems,
where switched systems naturally appear and create challenges that we believe
can be overcome by our method.

## References

  * [1] E. Sontag, _Stability and Stabilization: Discontinuities and the Effect of Disturbances_ , Nonlinear Analysis, Differential Equations and Control, 1999, pp. 551–598, Springer.
  * [2] A. Ahmadi and P. Parrilo, _Non-monotonic Lyapunov Functions for Stability of Discrete Time Nonlinear and Switched Systems_ , 2008 47th IEEE Conference on Decision and Control, pp. 614–621, 2008.
  * [3] A. Fürnsinn, C. Ebenbauer and B. Gharesifard, _Flexible-step Model Predictive Control based on Generalized Lyapunov Functions_ , Automatica (provisionally accepted), 2022.
  * [4] A. Fürnsinn, C. Ebenbauer and B. Gharesifard, _Relaxed Feasibility and Stability Criteria for Flexible-Step MPC_ , IEEE Control Systems Letters, vol. 7, pp. 2851–2856, 2023.
  * [5] J. Hespanha, D. Liberzon, D. Angeli, E. Sontag, _Nonlinear norm-observability notions and stability of switched systems_ , IEEE Transactions on Automatic Control, vol. 50, no. 2, pp. 154–168, 2005.
  * [6] Z. Artstein, _Stabilization with relaxed controls_ , Nonlinear Analysis: Theory, Methods & Applications, vol. 7, no. 11, pp. 1163–1173, 1983.
  * [7] D. Mayne and J. Rawlings and C. Rao and P. Scokaert, _Constrained model predictive control: Stability and optimality_ , Automatica, vol. 7, no. 6, pp. 789–814, 2000.
  * [8] D. Mayne, _Model predictive control: Recent developments and future promise_ , Automatica, vol. 50, no. 12, pp. 2967–2986, 2014.
  * [9] J. Rawlings and D. Mayne and M. Diehl, _Model Predictive Control: Theory, Computation, and Design_ , vol. 2, Nob Hill Publishing Madison, WI, 2017.
  * [10] L. Grüne and J. Pannek, _Nonlinear Model Predictive Control_ , Springer, 2017.
  * [11] D. Liberzon, _Switching in Systems and Control_ , vol. 190, Springer, 2003.
  * [12] E. Isaacson and H. Keller, _Analysis of Numerical Methods_ , Courier Corporation, 1994.
  * [13] M. Müller and K. Worthmann, _Quadratic costs do not always work in MPC_ , Automatica, vol. 82, pp. 269–277, 2017.