# Forecasting the changes between endemic and epidemic phases of a contagious
disease, with the example of COVID-19

Jacques Demongeot 1, Pierre Magal2111Corresponding author. e-mail:
<EMAIL_ADDRESS>and Kayode Oshnubi 1
1 AGEIS laboratory, UGA, 38700 La Tronche, France
2 Univ. Bordeaux, IMB, UMR 5251, F-33400 Talence, France.
CNRS, IMB, UMR 5251, F-33400 Talence, France.

###### Abstract

Predicting the endemic/epidemic transition during the temporal evolution of a
contagious disease.
Methods:  Defining indicators for detecting the transition endemic/epidemic,
with four scalars to be compared, calculated from the daily reported news
cases: coefficient of variation, skewness, kurtosis, and entropy. The
indicators selected are related to the shape of the empirical distribution of
the new cases observed over 14 days. This duration has been chosen to smooth
out the effect of weekends when fewer new cases are registered. For finding a
forecasting variable, we have used the PCA (principal component analysis),
whose first principal component (a linear combination of the selected
indicators) explains a large part of the observed variance and can then be
used as a predictor of the phenomenon studied (here the occurrence of an
epidemic wave).
Results: A score has been built from the four proposed indicators using a
Principal Component Analysis (PCA), which allows an acceptable level of
forecasting performance by giving a realistic retro-predicted date for the
rupture of the stationary endemic model corresponding to the entrance in the
epidemic exponential growth phase. This score is applied to the retro-
prediction of the limits of the different phases of the COVID-19 outbreak in
successive endemic/epidemic transitions and three countries, France, India,
and Japan.
Conclusion: We provided a new forecasting method for predicting an epidemic
wave occurring after an endemic phase for a contagious disease.

Keywords: Contagious disease; Endemic phase; Epidemic wave; Endemic/epidemic
transition forecasting; COVID-19 epidemic wave prediction

The paper is dedicated to James D. Murray, whose pioneering work in
mathematical biology we admire.

##  1 Introduction

Finding a reliable prediction method of the frontiers between different
stationary and non-stationary periods of a time series is a challenging
problem. Since the seminal work by Deshayes and Picard on the stationarity
rupture in time series [1, 2, 3], many works have dealt with the break in
stationarity [4, 5, 6, 7, 8, 9], the most recent using the concepts of
functional statistics [10, 11, 12, 13, 14, 15, 16]. Indeed, stationarity is
crucial as many forecasting models of time series rely on stationarity for
easy modeling and obtaining reliable results. A stationary time series
presents statistical properties which do not change over time, as the
empirical distribution of the random variable observed in the series, with its
main characteristic parameters, mean, coefficient of variation, moments, and
entropy. In the event of a break in stationarity, there may be a sudden
transition with a sudden change in the values of these parameters and the
appearance of a non-constant trend. The problem of the existence of this
transition arises with particular acuity in the case of contagious diseases,
which alternate stationary endemic periods and epidemic peaks having an
initial exponential trend, which must be predicted to prevent the spread of
the disease does not give rise to a pandemic.

The term endemic phase is understood to mean a period in which there is an
equilibrium in a model whose parameters have changed in value following an
epidemic phase, due to mitigation measures or a change in the virulence of the
infectious agent. In the case of chronic diseases observed outside epidemic
phases (bacterial meningitis, rabies, smallpox before its eradication, etc.),
the definition of endemicity corresponds more to the sporadic appearance (by
birth or human displacement) of individuals much more susceptible than the
general population to the infectious agent [17, 18, 19].

We will propose in this article a method to estimate the breakdown of endemic
stationarity based on four parameters linked to the empirical distribution of
the number of daily reported new cases of COVID-19 in several countries,
parameters whose isolated or joint predictive power will be analyzed. These
parameters are the coefficient of variation , the skewness, the kurtosis, and
the entropy of the stationary empirical measure calculated in a moving window.

The pair represented by the succession of an endemic phase and an epidemic
wave in COVID-19 outbreak can be considered as a functional unit, the break
between the two phases having to be found [20, 21, 22, 23, 24, 25]. The
endemic phase is characterized by a low average level of new cases, with low
variance. At the start of the epidemic phase, the average number of cases will
grow exponentially, and the standard deviation will grow proportionally at the
beginning, then saturate, like that of an additive noise in part independent
of growth, which explains the increase then decrease in the coefficient of
variation and kurtosis, therefore those of the first principal component at
the endemic/epidemic boundary.

##  2 Materials and Methods

We use in the following a moving window of length 14 days for calculating the
empirical distribution of the random variable equal to the number of daily
reported new cases. The empirical distribution $N_{t}$ on day $t$ is obtained
from the daily number of reported new cases considered as a random variable
$N_{t}=\left(N(t-13),N(t-12),\ldots,N(t)\right)$.

The length of the window has been chosen for eliminating the effect of the
week periodicity observed in data due to the lack of reporting each weekend.
Indeed, daily new infection cases are highly affected by weekends, such that
new case numbers are lowest at the start of the week, and increase afterwards
[26].

###  2.1 Empirical first four moments

In the following we will use the terminology endemic period, to describe a
period during which the daily new cases occurs randomly around some mean
value. An epidemic wave is a period during which the daily new cases occurs by
contact between susceptible and infected.

Our goal in this paper is to explore the transition between the endemic period
and the following epidemic wave which will be studied by calculating several
parameters in a moving window around the frontier of 14 days on which we
suspect this transition occurred.

We consider the first four moments of $N_{t}$. We start with the mean

$\mu=E(N_{t})=\dfrac{\displaystyle\sum_{i=0}^{13}N(t-i)}{14},$ (2.1)

where $E$ is the expectation operator, with the standard deviation

$\sigma=E\left(\left(N_{t}-\mu\right)^{2}\right)^{1/2}=\sqrt{\dfrac{\displaystyle\sum_{i=0}^{13}\left(N(t-i)-\mu\right)^{2}}{14}}.$
(2.2)

From these two first parameters, we can compute the coefficient of variation

$CV(N_{t})=\dfrac{\sigma}{\mu}.$ (2.3)

The skewness of the random variable $N_{t}$ is the third standardized moment,
defined as

$Skew(N_{t})=E\left(\left(\dfrac{N_{t}-\mu}{\sigma}\right)^{3}\right).$ (2.4)

Recall that the skewness verifies

$Skew(N_{t})=\dfrac{E(N_{t}^{3})-3\mu\sigma^{2}-\mu^{3}}{\sigma^{3}}=\dfrac{E(N_{t}^{3})}{\sigma^{3}}-\left(3\dfrac{1}{CV}+\dfrac{1}{CV^{3}}\right).$
(2.5)

The kurtosis is the fourth standardized moment, defined as

$Kurt(N_{t})=E\left(\left(\dfrac{N_{t}-\mu}{\sigma}\right)^{4}\right).$ (2.6)

The empirical entropy $\mathcal{E}$ of the empirical distribution is defined
as follows:

$\mathcal{E}(N_{t})=-\sum_{i=1:d\text{ with }p_{i}>0}p_{i}\log p_{i},$ (2.7)

where the $p_{i}$ are the weights of a histogram on d value intervals of
$N_{t}$. In the Results’ section, we use the approximate entropy. We refer to
[36] for more details.

###  2.2 Phenomenological model used for multiple epidemic waves

To represent the data, we used a phenomenological model to fit the curve of
cumulative reported cases. Such an idea is not new since it was already
proposed by Bernoulli [17] in 1760 in the context of the smallpox epidemic.
Here we used the so-called Bernoulli–Verhulst [37] model to describe the
epidemic phase. Bernoulli [17] investigated an epidemic phase followed by an
endemic phase. This appears clearly in Figures 9 and 10 of the paper by Dietz,
and Heesterbeek [38] who revisited the original article of Bernoulli. We also
refer to Blower [18] for another article revisiting the original work of
Bernoulli. Several works comparing cumulative reported cases data and the
Bernoulli–Verhulst model appear in the literature (see [39, 40, 41]). The
Bernoulli–Verhulst model is sometimes called Richard’s model, although
Richard’s work came much later in 1959.

The phenomenological model deals with data series of new infectious cases
decomposed into two successive phases: 1) endemic phases followed by 2)
epidemic phases.

Endemic phase: During the endemic phase, the dynamics of new cases appears to
fluctuate around an average value independently of the number of cases.
Therefore the average cumulative number of cases is given by

${\rm CR}(t)=N_{0}+(t-t_{0})\times a,\text{ for }t\in[t_{0},t_{1}],$ (2.8)

where $t_{0}$ denotes the beginning of the endemic phase, $N_{0}$ is the
number of new cases at time $t_{0}$, and $a$ is the average value of the daily
number of new cases.

We assume that the average daily number of new cases is constant. Therefore
the daily number of new cases is given by

${\rm CR}^{\prime}(t)=a.$ (2.9)

Epidemic phase: In the epidemic phase, the new cases are contributing to
produce secondary cases. Therefore the daily number of new cases is no longer
constant, but varies with time as follows

${\rm
CR}(t)=N_{\rm{base}}+\dfrac{\mathrm{e}^{\chi(t-t_{0})}N_{0}}{\left[1+\dfrac{N_{0}^{\theta}}{N_{\infty}^{\theta}}\left(\mathrm{e}^{\chi\theta\left(t-t_{0}\right)}-1\right)\right]^{1/\theta}},\text{
for }t\in[t_{0},t_{1}].$ (2.10)

In other words, the daily number of new cases follows the Bernoulli–Verhulst
[17, 37] equation. Namely, by setting

$N(t)={\rm CR}(t)-N_{\rm{base}},$ (2.11)

we obtain

$N^{\prime}(t)=\chi\,N(t)\,\left[1-\left(\dfrac{N(t)}{N_{\infty}}\right)^{\theta}\right],$
(2.12)

completed with the initial value

$N(t_{0})=N_{0}.$

In the model, $N_{\rm{base}}+N_{0}$ corresponds to the value ${\rm CR}(t_{0})$
of the cumulative number of cases at time $t=t_{0}$. The parameter
$N_{\infty}+N_{\rm{base}}$ is the maximal value of the cumulative reported
cases after the time $t=t_{0}$. $\chi>0$ is a Malthusian growth parameter, and
$\theta$ regulates the speed at which ${\rm CR}(t)$ increases to
$N_{\infty}+N_{\rm{base}}$.

##  3 Results

Here we use cumulative numbers of reported new cases for COVID-19 in France,
India, and Japan taken from WHO [42]. Data shows a succession of endemic
periods (yellow background color regions) followed by epidemic waves (blue
background color regions).

###  3.1 Data for France

Figure 1:  In this figure we plot in blue the phenomenological model and in
black the data. Data is the cumulative reported number of new cases with a
14-day rolling average.

In Figure 2, each colored segment corresponds to a single endemic or epidemic
phenomenological model. This change of color may occur due to a change of
dynamic inside an epidemic phase when a second wave comes before the end of
the previous one.

Figure 2:  In this figure we plot with multiple colors the phenomenological
models obtained for each period.  Figure 3:  In this figure we plot in blue
the first derivative of the phenomenological model and in black the data. Data
is the daily reported number of new cases with a 14-day rolling average.

By making a principal component analysis (i.e. the matlab function PCA)
between the standardized variables
$CV_{s}(N_{t}),Skew_{s}(N_{t}),Kurt_{s}(N_{t}),\mathcal{E}_{s}(N_{t})$, we
obtain the percentage of the variance explained by each principal component

$Explain=\left(\begin{array}[]{cccc}70.71\\\ 21.92\\\ 4.94\\\ 2.43\\\
\end{array}\right)$

and the matrix giving the projection coefficients of the principal components

$coeff=\left(\begin{array}[]{cccc}0.5527&-0.1480&0.7163&-0.3995\\\
0.5631&-0.2162&-0.0348&0.7968\\\ 0.5577&-0.0795&-0.6955&-0.4461\\\
0.2577&0.9618&0.0449&0.0808\end{array}\right).$

By using the first column of the above matrix, we deduce the first principal
component

$0.55CV_{s}(N_{t})+0.56Skew_{s}(N_{t})+0.56Kurt_{s}(N_{t})+0.26\mathcal{E}_{s}(N_{t})$
(3.1)

which explains $70.71\%$ of the variability.

We deduce that Kurtosis, Skewness, and the coefficient of variation (in
decreasing order of importance) best explain the variability.

Figure 4:  In this figure we plot the first principal component for France
(see formula (3.1)). The horizontal green lines correspond to $\pm 1$.

###  3.2 Data for India

In this subsection, we consider the data for India. We present the same curves
as for France, describing successively raw data and simulated results by the
phenomenological models.

Figure 5:  In this figure we plot in blue the phenomenological model and in
black the data. Data is the cumulative reported number of new cases with a
14-day rolling average.

In Figure 6, each colored segment corresponds to a single endemic or epidemic
phenomenological model. This change of color may occur due to a change of
dynamic inside an epidemic phase when a second wave comes before the end of
the previous one.

Figure 6:  In this figure we plot with multiple colors the phenomenological
models obtained for each period.  Figure 7:  In this figure we plot in blue
the first derivative of the phenomenological model and in black the data. Data
is the daily reported number of new cases with a 14-day rolling average.

By making a principal component analysis (i.e. the matlab function PCA)
between the standardized variables
$CV_{s}(N_{t}),Skew_{s}(N_{t}),Kurt_{s}(N_{t})$, and $\mathcal{E}_{s}(N_{t})$,
we obtain the percentage of the variance explained by each principal component

$Explain=\left(\begin{array}[]{cccc}53.13\\\ 22.55\\\ 14.22\\\ 10.09\\\
\end{array}\right)$

and the matrix giving the projection coefficients of the principal components

$coeff=\left(\begin{array}[]{cccc}0.5161&-0.2296&0.8212&-0.0810\\\
0.5714&-0.0587&-0.4434&-0.6881\\\ 0.5660&-0.2235&-0.3478&0.7133\\\
0.2946&0.9455&0.0896&0.1062\end{array}\right).$

By using the first column of the above matrix, we deduce the first principal
component

$0.52CV_{s}(N_{t})+0.57Skew_{s}(N_{t})+0.57Kurt_{s}(N_{t})+0.29\mathcal{E}_{s}(N_{t})$
(3.2)

which explains $53\%$ of the variability.

Figure 8:  In this figure we plot the first principal component for India (see
formula (3.2)). The horizontal green lines correspond to the values $\pm 1$.

###  3.3 Data for Japan

In this subsection, we consider the data for Japan. We present the same curves
as for France and India, describing successively raw data and simulated
results by the phenomenological models.

Figure 9:  In this figure we plot in blue the phenomenological model and in
black the data. Data is the cumulative reported number of new cases with a
14-day rolling average.  Figure 10:  In this figure we plot with multiple
colors the phenomenological models obtained for each period.  Figure 11:  In
this figure we plot in blue the first derivative of the phenomenological model
and in black the data. Data is the daily reported number of new cases with a
14-day rolling average.

By making a principal component analysis (i.e. the matlab function PCA)
between the standardized variables
$CV_{s}(N_{t}),Skew_{s}(N_{t}),Kurt_{s}(N_{t})$, and $\mathcal{E}_{s}(N_{t})$,
we obtain the percentage of the variance explained by each principal component

$Explain=\left(\begin{array}[]{cccc}71.62\\\ 17.54\\\ 6.87\\\
3.97\end{array}\right)$

and the matrix giving the projection coefficients of the principal components

$coeff=\left(\begin{array}[]{cccc}0.5234&-0.2243&0.7969&-0.2018\\\
0.5452&-0.2406&-0.2310&0.7691\\\ 0.5401&-0.1760&-0.5575&-0.6054\\\
0.3703&0.9278&0.0270&0.0358\par\end{array}\right).$

By using the first column of the above matrix, we deduce the first principal
component

$0.52CV_{s}(N_{t})+0.55Skew_{s}(N_{t})+0.54Kurt_{s}(N_{t})+0.37\mathcal{E}_{s}(N_{t})$
(3.3)

which explains $71\%$ of the variability.

We deduce that Kurtosis, Skewness, and the coefficient of variation (in
decreasing order of importance) best explain the variability.

Figure 12:  In this figure we plot the first principal component for Japan
(see formula (3.3)). The horizontal green lines correspond to the values $\pm
1$.

##  4 Discussion

The forecasting of the epidemic waves of the COVID-19 outbreak is based on a
change in the nature of the time series dynamics related to the number of
daily new reported cases of this contagious disease. This change can concern
the moments or the entropy of the empirical distribution of the stationary
component at the end of the endemic phase, which disrupts when a not constant
trend occurs, marking the start of an epidemic wave.

From a careful examination of Figures 4, 8 and 12, we can conclude that there
are not constant but frequent patterns for $C_{1}(t)$ identifiable in the
three studied countries and for a majority of their endemic/epidemic
transitions.

The predictive power of the first principal component $C_{1}(t)$ can be
quantified by its performance ratio, that is, by the percentage of correct
retro-predictions obtained by fixing variation thresholds to forecast the
occurrence of an epidemic wave. For France, if we fix the threshold to the
value $1$, $C_{1}(t)$ predicts correctly an epidemic outbreak a weak after a
decrease from this threshold value not reached elsewhere in the endemic phase.
This prediction is correct at $53\%$ only for India. For Japan, the
performance of the retro-prediction is $71\%$ and $70\%$ for France.

It is clear that the level of prediction is not very high (71% in the best
case), but, in the absence of a currently reliable predictor, we can consider
that it is sufficient to trigger mitigation measures at level of a population.
A more systematic study of the changing shape of the empirical distribution is
needed, looking at many epidemic waves in many countries. A parameter
measuring the deviation from classical laws (such as those linked to the
Kolmogorov-Smirnov, chi-square or Shapiro-Wilk tests in the case of the normal
law) could thus be added in subsequent works.

As previously noticed in [43], a classical epidemic peak with a near-symmetric
growth and decline may be preceded or followed by a shoulder-like behavior
corresponding to a prematurely stopped wave followed by another. Shoulder-like
behavior for epidemic waves can be explained by using multiple sub-group
epidemic models. This idea was first explored for SARS-CoV-1 in [44] and
reconsidered for SARS-CoV-2 by [45]. But the changes between endemic and
epidemic are still challenging to model.

##  5 Conclusions

We have studied in this article the evolution of four parameters related to
the dynamics of the number of the daily reported new cases $N_{t}$ at day $t$
of a contagious disease, which can serve as early indicators of the appearance
of epidemic waves from a previous endemic state. By applying this parameter
calculation to COVID-19, we showed that a score obtained by PCA based on the
linear combination of the four chosen parameters with specific coefficients
for each one could forecast the variations of the empirical distribution of
the daily reported number of new cases $N_{t}$, then could be considered often
as a good predictor (for the countries and the epidemic waves in these
countries) of the endemic-epidemic transition. A systematic study of
contagious diseases other than COVID-19 is necessary to confirm this
forecasting property’s existence. Still, we can already propose this score as
a realistic indicator of the next occurrence of an epidemic outbreak from a
change in the dynamics of the observed daily new cases during the endemic
periods.

## References

  * [1] Deshayes, J, & Picard, D. Application aux tests de rupture de régression. Astérisque 1979, 68, 73-98.
  * [2] Deshayes, J., & Picard, D. Lois asymptotiques des tests et estimateurs de rupture dans un modèle statistique classique. Annales de l’I.H.P. Probabilités et statistiques 1984, 20, 309-327.
  * [3] Picard, D. Testing and estimating change-points in time series. Advances in applied probability, 1985, 17(4), 841-867.
  * [4] Vogt, M., & Dette, H. Detecting gradual changes in locally stationary processes. The Annals of Statistics 2015, 43, 713-740.
  * [5] van Delft, A., & Eichler, M. Locally stationary functional time series. Electronic Journal of Statistics 2018, 12, 107-170.
  * [6] Palachy, S. Detecting stationarity in time series data. Medium Towards Data Sci. (2019), 9, 53.
  * [7] Ting, K.M., Liu, Z., Zhang, H., & Zhu, Y. A new distributional treatment for time series and an anomaly detection investigation. Proceedings of the VLDB Endowment (2022), 15, 2321-2333.
  * [8] Hauber, A.L., Sigloch, C., & Timmer, J. Detecting frequency modulation in stochastic time-series data. Physical Review E (2022), 106, 024204.
  * [9] Bawdekar, A.A., Prusty, B.R., & Bingi, K. Sensitivity Analysis of Stationarity Tests’ Outcome to Time Series Facets and Test Parameters. Mathematical Problems in Engineering 2022, 2022, 2402989.
  * [10] Demongeot, J., Laksaci, A., Madani, F., & Rachdi, M. Functional data: local linear estimation of the conditional density and its application. Statistics (2013), 47, 26-44.
  * [11] Rachdi, M., Laksaci, A., Demongeot, J., Abdali, A., & Madani, F. Theoretical and practical aspects on the quadratic error in the local linear estimation of the conditional density for functional data. Comp. Stat. Data Anal. (2014), 73, 53-68.
  * [12] Demongeot, J., Laksaci, A., Rachdi, M., & Rahmani, S. On the local linear modelization of the conditional distribution for functional data. Sankhya A (2014), 76, 328-355.
  * [13] Demongeot, J., Hamie, A., Laksaci, A., & Rachdi, M. Relative-Error Prediction in Nonparametric Functional Statistics: Theory and Practice. Journal of Multivariate Analysis (2016), 146, 261-268.
  * [14] Belkis, A., Demongeot, J., Laksaci, A., & Rachdi, M. Functional data analysis: Estimation of the relative error in functional re-gression under random left-truncation. J. Nonparametric Statistics (2018), 30, 472-490.
  * [15] Oshinubi, K., Ibrahim, F., Rachdi, M., & Demongeot, J. Functional Data Analysis: Application to Daily Observation of COVID-19 Prevalence in France. AIMS Mathematics (2022), 7, 5347-5385.
  * [16] Jelassi, M., Oshinubi, K., Rachdi, M., & Demongeot, J. Epidemic Dynamics on Social Interaction Networks. AIMS Bioengineering (2022), 9, 348-361.
  * [17] Bernoulli, D. Essai d’une nouvelle analyse de la petite Vérole, & des avantages de l’Inoculation pour la prévenir, Mémoire Académie Royale des Sciences, Paris, 1760.
  * [18] Blower S, & Bernoulli D. An attempt at a new analysis of the mortality caused by smallpox and of the advantages of inoculation to prevent it. (1766). Rev Med Virol. (2004) Sep-Oct;14(5):275-88.
  * [19] Murray J.D., Stanley E.A., Brown D.L. On the spatial spread of rabies among foxes. Proc R Soc Lond B Biol Sci. (1986) Nov 22;229 (1255):111-50.
  * [20] Cinaglia, P., & Cannataro, M. Forecasting COVID-19 epidemic trends by combining a neural network with $R_{t}$ estimation. Entropy (2022), 24(7), 929.
  * [21] Srivastava, S. R., Meena, Y. K., & Singh, G. Forecasting on Covid-19 infection waves using a rough set filter driven moving average models. Appl. Soft. Comput. (2022) Dec;131:109750.
  * [22] Lynch, C. J., & Gore, R. Short-range forecasting of COVID-19 during early onset at county, health district, and state geographic levels using seven methods: comparative forecasting study. Journal of medical Internet research, 23(3) (2021), e24925.
  * [23] Lynch C.J., & Gore R. Application of one-, three-, and seven-day forecasts during early onset on the COVID-19 epidemic dataset using moving average, autoregressive, autoregressive moving average, autoregressive integrated moving average, and naïve forecasting methods. Data Brief. (2021) Apr;35:106759.
  * [24] Yu CS, Chang SS, Chang TH, Wu JL, Lin YJ, Chien HF, & Chen RJ. A COVID-19 Pandemic Artificial Intelligence-Based System With Deep Learning Forecasting and Automatic Statistical Data Acquisition: Development and Implementation Study. J Med Internet Res. (2021) May 20;23(5):e27806.
  * [25] ArunKumar KE, Kalaga DV, Sai Kumar CM, Chilkoor G, Kawaji M, & Brenza TM. Forecasting the dynamics of cumulative COVID-19 cases (confirmed, recovered and deaths) for top-16 countries using statistical machine learning models: Auto-Regressive Integrated Moving Average (ARIMA) and Seasonal Auto-Regressive Integrated Moving Average (SARIMA). Appl Soft Comput. 2021 May;103:107161.
  * [26] Demongeot, J., Oshinubi, K., Rachdi, M., Seligmann, H., Thuderoz, F., & Waku, J. Estimation of daily reproduction rates in COVID-19 outbreak, Computation, 9, 109 (2021).
  * [27] Demongeot, J., & Magal, P. Spectral method in epidemic time series. Biology (2022), 11, 1825.
  * [28] Worldometer. Available online: https://www.worldometers.info/coronavirus (Accessed on 6 January 2023).
  * [29] Caswell, H. Matrix population models: construction, analysis and interpretation, 2nd ed. Sinauer: Sunderland, (2001).
  * [30] Demongeot, J., & Demetrius, L. Complexity and Stability in Biological Systems. Int. J. Bifurcation and Chaos (2015), 25, 40013.
  * [31] Nicol-Harper, A., Dooley, C., Packman, D., Mueller, M., Bijak, J., Hodgson, D., Townley, S., & Ezard, T. Inferring transient dynamics of human populations from matrix non-normality. Popul Ecol. (2018), 60, 185-196.
  * [32] Xu, Z., Yang, D., Wang, J., & Demongeot, J. Statistical Analysis Supports UTR (untranslated region) Deletion Theory in SARS-CoV-2. Virulence (2022), 13, 1772-1789.
  * [33] Demongeot, J., Griette, Q., Maday, Y., & Magal, P. A Kermack-McKendrick model with age of infection starting from a single or multiple cohorts of infected patients. Proc. Royal Society A (2023), 479, 2022.0381.
  * [34] Kammegne, B., Oshinubi, K., Babasola, T., Peter, O.J., Longe, O.B., Ogunrinde, R.B., Titiloye, E.O., & Demongeot, J. Mathematical modelling of spatial distribution of COVID-19 outbreak using diffusion equation. Pathogens (2023), 12, 88.
  * [35] Xu, Z., Wei, D., Zeng, Q., Zhang, H., Sun, Y., & Demongeot, J. More or less deadly? A mathematical model that predicts 1 SARS-CoV-2 evolutionary direction. Computers in Biology & Medicine 2023, 153, 106510.
  * [36] Pincus, S. M. Approximate entropy as a measure of system complexity. Proceedings of the National Academy of Sciences, 88(6) (1991), 2297-2301.
  * [37] Verhulst, P. F. Notice sur la loi que la population poursuit dans son accroissement, Corresp. Mathématique Phys., (1938),113-121.
  * [38] Dietz, K., & Heesterbeek, J. A. P. Daniel Bernoulli’s epidemiological model revisited, Math. Biosci., 180 (2002), 1-21.
  * [39] Hsieh, Y. H. Richards model: a simple procedure for real-time prediction of outbreak severity, in Modeling and dynamics of infectious diseases, (2009), 216–236.
  * [40] Wang, X.S., Wu, J., & Yang, Y. Richards model revisited: validation by and application to infection dynamics, J. Theoret. Biol., 313 (2012), 12-19.
  * [41] Zhou, G., & Yan, G. Severe acute respiratory syndrome epidemic in Asia, Emerg. Infect. Dis., 9 (2003), 1608-1610.
  * [42] Data from WHO. Accessed: July 20, 2022. https://COVID19.who.int/WHO-COVID-19-global-data.csv
  * [43] Weitz, J. S., Park, S. W., Eksin, C., & Dushoff, J. Awareness-driven behavior changes can shift the shape of epidemics away from peaks and toward plateaus, shoulders, and oscillations. Proceedings of the National Academy of Sciences, 117(51) (2020), 32764-32771.
  * [44] Magal, P., Seydi, O., & Webb, G. Final size of a multi-group SIR epidemic model: Irreducible and non-irreducible modes of transmission. Mathematical biosciences, 301 (2018), 59-67.
  * [45] Chowell, G., Dahal, S., Tariq, A., Roosa, K., Hyman, J. M., & Luo, R. An ensemble n-sub-epidemic modeling framework for short-term forecasting epidemic trajectories: Application to the COVID-19 pandemic in the USA. PLoS Computational Biology, 18(10) (2022), e1010602.


# Regret and Conservatism of Constrained Stochastic Model Predictive Control

Maik Pfefferkorn<EMAIL_ADDRESS>
Control and Cyber-Physical Systems Laboratory
Technical University of Darmstadt, Germany Venkatraman Renganathan
<EMAIL_ADDRESS>
Department of Automatic Control - LTH
Lund University, Sweden Rolf Findeisen<EMAIL_ADDRESS>
Control and Cyber-Physical Systems Laboratory
Technical University of Darmstadt, Germany

###### Abstract

We analyse conservatism and regret of stochastic model predictive control
(SMPC) when using moment-based ambiguity sets for modeling unknown
uncertainties. To quantify the conservatism, we compare the deterministic
constraint tightening while taking a distributionally robust approach against
the optimal tightening when the exact distributions of the stochastic
uncertainties are known. Furthermore, we quantify the regret by comparing the
performance when the distributions of the stochastic uncertainties are known
and unknown. Analysing the accumulated sub-optimality of SMPC due to the lack
of knowledge about the true distributions of the uncertainties marks the novel
contribution of this work.

##  1 Introduction

Recently, there has been a surge of interest in analyzing control algorithms
that operate subject to unknown entities through the lens of regret analysis.
Entities of interest, also referred to as _quantities of interest_ (QIs), are
numerous, e.g., disturbance sequences, uncertainty distributions or cost
function parameterizations. Lack of knowledge about a QI induces regret for
the controller, i.e., a performance loss compared to when the QI was known.
Analysing regret enables to better understand the performance trade-off of
control algorithms due to the lack of knowledge about the entities. Notably,
many robust control algorithms incorporate caution mechanisms into their
decision making scheme in order to account for the lack of knowledge about the
QIs. The cautiousness in the decision making can be quantified as the
conservatism of the control algorithm. As a result of cautious decision
making, robust control algorithms are known to incur regret. For instance, the
regret associated with $H_{\infty}$ control compared against an optimal
hindsight controller that knows the disturbance sequences exactly beforehand
was examined in Karapetyan et al. (2022). Further works investigated regret in
different variants of the finite- and infinite-horizon linear-quadratic
regulator problem, see e.g. Li et al. (2019); Chen et al. (2023) and
references therein.

We are interested in comparing controllers that are robust with respect to a
moment-based ambiguity set of distributions against a fully informed
counterpart that knows the true uncertainty distribution. Specifically, we
focus on stochastic model predictive control (SMPC) Rawlings et al. (2019).
For SMPC, we analyze conservatism and regret of its distributionally robust
formulation for unknown uncertainty distributions compared against its fully
informed counterpart. We provide a detailed performance and constraint
tightening analysis for SMPC. Our main contributions are as follows:

  1. 1.

We define constraint conservatism and distributional regret of
distributionally robust SMPC with uncertainties modeled using moment-based
ambiguity sets.

  2. 2.

We develop a framework for quantifying conservatism and regret and for
analyzing their behaviors. Specifically, we derive analytic expressions for
conservatism and regret to identify and analyze the effects that lead to
regret. In doing so, we show how regret analysis can open the door for
leveraging the power of dynamically setting state and input constraints in
SMPC to achieve better performance.

  3. 3.

We underline our theoretical findings in simulations.

The paper is organised as follows: The SMPC problem is introduced in Section
2. Regret and conservatism associated with SMPC are presented in Section 3,
while some of its features are demonstrated using numerical simulation in
Section 4 before concluding in Section 5.

Notation and Preliminaries. The set of real numbers, integers and the natural
numbers are denoted by $\mathbb{R},\mathbb{Z},\mathbb{N}$ respectively and the
subset of real numbers greater than a given constant $a\in\mathbb{R}$ is
denoted by $\mathbb{R}_{>a}$. The subset of natural numbers between two
constants $a,b\in\mathbb{N}$ with $a<b$ is denoted by $[a:b]$. Given two sets
$A\subset\mathbb{R}^{n},B\subset\mathbb{R}^{n}$, their Pontryagin difference
is denoted by $A\ominus B:=\\{a\in A\mid a+b\in A,\forall b\in
B\\}\subset\mathbb{R}^{n}$. For a matrix $A\in\mathbb{R}^{n\times n}$, we
denote its transpose and its trace by $A^{\top}$ and $\mathbf{Tr}(A)$
respectively. An identity matrix of dimension $n$ is denoted by $I_{n}$. We
denote by $\mathbb{S}^{n}$ the set of symmetric matrices in
$\mathbb{R}^{n\times n}$. For $A\in\mathbb{S}^{n}$, we denote by $A\succ
0(A\succeq 0)$ to mean that $A$ is positive definite (positive semi-definite).
Given $x\in\mathbb{R}^{n}$, and $A\in\mathbb{R}^{n\times n}$, the notation
${\left\|x\right\|}^{2}_{A}$ denotes $x^{\top}Ax$.

##  2 Stochastic Model Predictive Control (SMPC) Problem Formulation

We consider SMPC with joint chance constraints similar to Paulson et al.
(2020) using moment-based ambiguity set modeling to describe the stochastic
system uncertainties. The detailed set-up is as follows.

###  2.1 System Dynamics & Constraints

We consider a stochastic discrete-time linear time-invariant system given by

$x_{k+1}=Ax_{k}+Bu_{k}+w_{k},\quad\forall k\in\mathbb{N},$ (1)

where $x_{k}\in\mathbb{R}^{n}$ and $u_{k}\in\mathbb{R}^{m}$ are the system
state and input at time step $k$, respectively. The matrices
$A\in\mathbb{R}^{n\times n}$ and $B\in\mathbb{R}^{n\times m}$ denote the
dynamics matrix and the input matrix respectively. For the ease of exposition,
we will assume that $x_{0}\in\mathbb{R}^{n}$ is known. The process noise
$w_{k}\in\mathbb{R}^{n}$ is a zero-mean random vector that is independent and
identically distributed over time. The distribution of $w_{k}$, namely
$\mathbb{P}^{w}$, is unknown. However, it is known to belong to a moment-based
ambiguity set of distributions given by

$\mathcal{P}^{w}=\left\\{\mathbb{P}^{w}\mid\mathbb{E}[w_{k}]=0,\mathbb{E}[w_{k}w_{k}^{\top}]={\Sigma}_{w}\right\\}.$
(2)

With $N\in\mathbb{N}$, we denote the prediction horizon considered in the
predictive control problem to be defined shortly. The states, inputs and the
disturbances over the prediction horizon are given by

$\displaystyle\mathbf{x}_{k}$ $\displaystyle:=\begin{bmatrix}x^{\top}_{0\mid
k}&x^{\top}_{1\mid k}&\cdots&x^{\top}_{N\mid
k}\end{bmatrix}^{\top}\in\mathbb{R}^{(N+1)n},$ (3)
$\displaystyle\mathbf{u}_{k}$ $\displaystyle:=\begin{bmatrix}u^{\top}_{0\mid
k}&u^{\top}_{1\mid k}&\cdots&u^{\top}_{N-1\mid
k}\end{bmatrix}^{\top}\in\mathbb{R}^{Nm},$ (4) $\displaystyle\mathbf{w}_{k}$
$\displaystyle:=\begin{bmatrix}w^{\top}_{0\mid k}&w^{\top}_{1\mid
k}&\cdots&w^{\top}_{N-1\mid k}\end{bmatrix}^{\top}\in\mathbb{R}^{Nn}.$ (5)

Here, the subscript $i\\!\mid\\!k$ indicates $i$ time steps ahead of $k$ and
$x_{0\mid k}=x_{k}$. Starting at $x_{k}$, we compactly write the evolution of
(1) over the prediction horizon as

$\displaystyle\mathbf{x}_{k}=\mathbf{A}x_{0\mid
k}+\mathbf{B}\mathbf{u}_{k}+\mathbf{D}\mathbf{w}_{k},$ (6)

with matrices $\mathbf{A},\mathbf{B}$, and $\mathbf{D}$ being of appropriate
dimension. The mean and covariance of (6) evolve as

$\displaystyle\bar{\mathbf{x}}_{k}$
$\displaystyle:=\mathbb{E}[\mathbf{x}_{k}]=\mathbf{A}x_{0\mid
k}+\mathbf{B}\mathbf{u}_{k},$ (7) $\displaystyle\mathbf{\Sigma}_{\mathbf{x}}$
$\displaystyle:=\mathbb{E}[(\mathbf{x}_{k}-\bar{\mathbf{x}}_{k})(\mathbf{x}_{k}-\bar{\mathbf{x}}_{k})^{\top}]=\mathbf{D}\mathbf{\Sigma}_{\mathbf{w}}\mathbf{D}^{\top},$
(8)

where $\mathbf{\Sigma}_{\mathbf{w}}$ is a block diagonal matrix with $N$
blocks and each block being equal to $\Sigma_{w}$. Since the support of
$\mathbb{P}^{w}$ can be unbounded, we cannot guarantee the satisfaction of
hard state constraints for all $k$. Hence, we enforce a distributionally
robust joint risk constraint on the states as

$\displaystyle\sup_{\mathbb{P}_{k}^{\mathbf{x}}\in\mathcal{P}^{\mathbf{x}}_{k}}\mathbb{P}_{k}^{\mathbf{x}}\left[\mathbf{x}_{k}\notin\mathcal{X}\right]\leq\Delta,$
(9)

where $\mathcal{P}^{\mathbf{x}}_{k}$ denotes the moment-based ambiguity set
for the compact state $\mathbf{x}_{k}$ and is given by

$\displaystyle\mathcal{P}^{\mathbf{x}}_{k}=\left\\{\mathbb{P}^{\mathbf{x}}_{k}\mid\mathbb{E}[\mathbf{x}_{k}]=\bar{\mathbf{x}}_{k},\mathbb{E}[(\mathbf{x}_{k}\\!-\\!\bar{\mathbf{x}}_{k})(\mathbf{x}_{k}-\bar{\mathbf{x}}_{k})^{\top}]=\mathbf{\Sigma}_{\mathbf{x}}\right\\}.$

The set $\mathcal{X}$ is assumed to be a convex polytope defined by

$\displaystyle\mathcal{X}$
$\displaystyle:=\bigcap^{n_{x}}_{i=1}\left\\{\mathbf{x}\mid
f^{\top}_{i}\mathbf{x}\leq
g_{i}\right\\}=\left\\{\mathbf{x}\mid\mathbf{F}\mathbf{x}\leq\mathbf{g}\right\\},$
(10)

where
$f_{i}\in\mathbb{R}^{(N+1)n},g_{i}\in\mathbb{R},\mathbf{F}\in\mathbb{R}^{n_{x}\times(N+1)n}$
and $\mathbf{g}\in\mathbb{R}^{n_{x}}$. Further, $\Delta\in(0,0.5]$ represents
a user-prescribed total risk budget for the worst-case probability of
constraint violation over the entire prediction horizon. A similar
probabilistic treatment of control constraints can also be formulated. For the
ease of exposition, we consider the hard input constraint formulation
$u_{k}\in\mathcal{U},\forall k$, where $\mathcal{U}$ is a convex polytope

$\displaystyle\mathcal{U}$
$\displaystyle:=\bigcap^{n_{u}}_{j=1}\left\\{\mathbf{u}\mid
c^{\top}_{j}\mathbf{u}\leq
d_{j}\right\\}=\left\\{\mathbf{u}\mid\mathbf{C}\mathbf{u}\leq\mathbf{d}\right\\}$
(11)

with
$c_{j}\in\mathbb{R}^{Nm},d_{i}\in\mathbb{R},\mathbf{C}\in\mathbb{R}^{n_{u}\times
Nm}$ and $\mathbf{d}\in\mathbb{R}^{n_{u}}$. The objective of the SMPC is to
minimize the cost

$\displaystyle J_{\mathrm{SMPC}}(\mathbf{u}_{k},x_{0\mid
k}):=\mathbb{E}\left[\left\lVert{\mathbf{x}_{k}}\right\rVert^{2}_{\mathbf{Q}}+\left\lVert{\mathbf{u}_{k}}\right\rVert^{2}_{\mathbf{R}}\right].$
(12)

Note that $\mathbf{Q}$ is a block diagonal matrix with $N+1$ blocks and each
block being equal to $Q\succeq 0$. Similarly, $\mathbf{R}$ is a block diagonal
matrix with $N$ blocks and each block being equal to $R\succ 0$. Then,
$J_{\mathrm{SMPC}}(\mathbf{u}_{k},x_{0\mid k})$ can be rewritten as

$\displaystyle J_{\mathrm{SMPC}}(\mathbf{u}_{k},x_{0\mid
k})=\left\lVert{\bar{\mathbf{x}}_{k}}\right\rVert^{2}_{\mathbf{Q}}+\left\lVert{\mathbf{u}_{k}}\right\rVert^{2}_{\mathbf{R}}+\mathbf{Tr}(\mathbf{Q}\mathbf{\Sigma}_{\mathbf{x}}).$
(13)

The following problem formally establishes the SMPC optimization problem along
with the state and input constraints.

###  2.2 Stochastic Model Predictive Control (SMPC)

######  Problem 1.

Given an initial state $x_{0|k}=x_{k}\in\mathbb{R}^{n}$, process noise
$w_{i\mid k}\sim\mathbb{P}^{w}\in\mathcal{P}^{w},\forall i\in[0:N-1]$, and the
penalty matrices $Q\succeq 0,R\succ 0$, we seek to find an input sequence
$\mathbf{u}_{k}$ that optimizes the following optimal control problem (OCP):

$\displaystyle\\!\underset{\mathbf{u}_{k}}{\mathrm{min}}$ $\displaystyle
J_{\mathrm{SMPC}}(\mathbf{u}_{k},x_{0\mid k})$ (14a)
$\displaystyle\mathrm{s.}\,\mathrm{t.}$
$\displaystyle\mathbf{x}_{k}=\mathbf{A}x_{0\mid
k}+\mathbf{B}\mathbf{u}_{k}+\mathbf{D}\mathbf{w}_{k},\quad x_{0\mid k}=x_{k},$
(14b) $\displaystyle\mathbf{u}_{k}\in\mathcal{U},w_{i\mid
k}\sim\mathbb{P}^{w}\in\mathcal{P}^{w},$ (14c)
$\displaystyle\sup_{\mathbb{P}_{k}^{\mathbf{x}}\in\mathcal{P}^{\mathbf{x}}_{k}}\mathbb{P}_{k}^{\mathbf{x}}\left[\mathbf{x}_{k}\notin\mathcal{X}\right]\leq\Delta.$
(14d)

Problem 1 is solved online in a receding horizon fashion at each time instance
$k$. The first element $u^{\star}_{0\mid k}$ of the optimal input sequence
$\mathbf{u}^{\star}_{k}$ is applied to the system given by (1) until the next
sampling instant $k+1$111We refer the reader to Rawlings et al. (2019) for a
comprehensive introduction to and overview of model predictive control.. Note
that the OCP (14) is in general computationally intractable as (14d) is an
infinite-dimensional distributionally robust joint risk constraint. Using
Boole’s inequality, (14d) can be decomposed into individual risk constraints
across each time step along the prediction horizon. Suppose that
$\sum^{n_{x}}_{i=1}\delta_{i}\leq\Delta$ and

$\displaystyle\sup_{\mathbb{P}_{k}^{\mathbf{x}}\in\mathcal{P}^{\mathbf{x}}_{k}}\mathbb{P}_{k}^{\mathbf{x}}\left[f^{\top}_{i}\mathbf{x}_{k}>g_{i}\right]$
$\displaystyle\leq\delta_{i},\quad\forall i=1,\dots,n_{x}.$ (15)

Then, (15) implies (14d). That is,

$\displaystyle\sup_{\mathbb{P}_{k}^{\mathbf{x}}\in\mathcal{P}^{\mathbf{x}}_{k}}\mathbb{P}_{k}^{\mathbf{x}}\\!\left[\mathbf{x}_{k}\\!\notin\\!\mathcal{X}\right]$
$\displaystyle=\sup_{\mathbb{P}_{k}^{\mathbf{x}}\in\mathcal{P}^{\mathbf{x}}_{k}}\mathbb{P}_{k}^{\mathbf{x}}\left[\mathbf{x}_{k}\\!\in\\!\bigcup^{n_{x}}_{i=1}\left\\{\mathbf{x}_{k}\\!\mid\\!f^{\top}_{i}\mathbf{x}_{k}\\!>\\!g_{i}\right\\}\right]$
$\displaystyle\leq\sum^{n_{x}}_{i=1}\sup_{\mathbb{P}_{k}^{\mathbf{x}}\in\mathcal{P}^{\mathbf{x}}_{k}}\mathbb{P}_{k}^{\mathbf{x}}\left[f^{\top}_{i}\mathbf{x}_{k}>g_{i}\right]$
$\displaystyle\leq\sum^{n_{x}}_{i=1}\delta_{i}$ $\displaystyle\leq\Delta.$

Though (15) is an infinite-dimensional risk constraint, it can be equivalently
formulated as a second order cone constraint on the state mean through
deterministic constraint tightening using Lemma 1 or Lemma 2 given below
without proofs.

######  Lemma 1.

(From Ono & Williams (2008)) Let $\mathbf{x}\in\mathbb{R}^{n}$ be a random
vector with known distribution $\mathbb{P}^{\mathbf{x}}$ defined using mean
$\bar{\mathbf{x}}$ and covariance $\mathbf{\Sigma}^{\mathbf{x}}$. Then,
$\forall\delta_{i}\\!\in\\!(0,0.5]$, and $i\\!\in\\![1:n_{x}]$, we have

$\footnotesize\mathbb{P}^{\mathbf{x}}\left[f^{\top}_{i}\mathbf{x}>g_{i}\right]\leq\delta_{i}\Leftrightarrow
f^{\top}_{i}\bar{\mathbf{x}}\leq
g_{i}-\Phi^{-1}_{\mathbb{P}^{\mathbf{x}}}(1-\delta_{i})\left\lVert{\mathbf{\Sigma}_{\mathbf{x}}^{\frac{1}{2}}f_{i}}\right\rVert_{2},$

where $\Phi_{\mathbb{P}^{\mathbf{x}}}$ denotes the Cumulative Distribution
Function (CDF) of the known distribution $(\mathbb{P}^{\mathbf{x}})$ when
normalized.

######  Lemma 2.

(From Calafiore & El Ghaoui (2006)) Let $\mathbf{x}\in\mathbb{R}^{n}$ be
random with mean $\bar{\mathbf{x}}$ and covariance
$\mathbf{\Sigma}_{\mathbf{x}}$. Then,
$\forall\delta_{i}\in(0,0.5],i\in[1:n_{x}]$,

$\sup_{\mathbb{P}^{\mathbf{x}}\in\mathcal{P}^{\mathbf{x}}}\mathbb{P}^{\mathbf{x}}\left[f^{\top}_{i}\mathbf{x}\\!>\\!g_{i}\right]\leq\delta_{i}\Leftrightarrow
f^{\top}_{i}\bar{\mathbf{x}}\\!\leq\\!g_{i}-\sqrt{\frac{1\\!-\\!\delta_{i}}{\delta_{i}}}\left\lVert{\mathbf{\Sigma}_{\mathbf{x}}^{\frac{1}{2}}f_{i}}\right\rVert_{2}.$

The surrogate of problem 1, called the surrogate SMPC, can be written as

$\displaystyle\\!\underset{\mathbf{u}_{k}}{\mathrm{min}}$ $\displaystyle
J_{\mathrm{SMPC}}(\mathbf{u}_{k},x_{0\mid k})$ (16a)
$\displaystyle\mathrm{s.}\,\mathrm{t.}$
$\displaystyle\mathbf{x}_{k}=\mathbf{A}x_{0\mid
k}+\mathbf{B}\mathbf{u}_{k}+\mathbf{D}\mathbf{w}_{k},\quad x_{0\mid k}=x_{k},$
(16b)
$\displaystyle\mathbf{u}_{k}\in\mathcal{U},w_{k}\sim\mathbb{P}^{w}\in\mathcal{P}^{w},$
(16c) $\displaystyle f^{\top}_{i}\bar{\mathbf{x}}_{k}\leq
g_{i}-\psi_{i}\left\lVert{\mathbf{\Sigma}_{\mathbf{x}}^{\frac{1}{2}}f_{i}}\right\rVert_{2},i\in[1:n_{x}],$
(16d)

Here, the tightening constant $\psi_{i},\forall i=1,\dots,n_{x}$ is given by

$\displaystyle\psi_{i}:=\begin{cases}\Phi^{-1}_{\mathbb{P}_{k}^{\mathbf{x}}}(1-\delta_{i}),&\text{
when $\mathbb{P}_{k}^{\mathbf{x}}$ is known},\\\
\sqrt{\frac{1-\delta_{i}}{\delta_{i}}},&\text{ when
$\mathbb{P}_{k}^{\mathbf{x}}$ is unknown}.\end{cases}$ (17)

Note that (16) has finite-dimensional constraints unlike its original
counterpart in (14). Let the control input sequences that minimize (16) with
exact constraint tightening and distributionally robust constraint tightening
according to (17) be denoted by $\mathbf{u}^{*}_{k}$ and
$\mathbf{u}^{\dagger}_{k}$ respectively. Their corresponding optimal value
functions are given by

$\displaystyle J^{*}_{\mathrm{SMPC}}(x_{k})$
$\displaystyle:=J_{\mathrm{SMPC}}(\mathbf{u}^{*}_{k},x_{k}).$ (18)
$\displaystyle J^{\dagger}_{\mathrm{SMPC}}(x_{k})$
$\displaystyle:=J_{\mathrm{SMPC}}(\mathbf{u}^{\dagger}_{k},x_{k}).$ (19)

##  3 Conservatism & Regret Analyses

In this section, we introduce and define the concepts of conservatism and
regret for the previously introduced SMPC algorithm.

###  3.1 Conservatism of SMPC

We would like to study the difference in constraint tightening when the true
distributions of the stochastic uncertainties are known and when they are
unknown. To do so, we define the constraint conservatism associated with the
SMPC problem with joint chance constraint formulation as follows.

######  Definition 1.

The constraint conservatism, denoted by
$\mathfrak{C}_{\mathrm{SMPC}}\in\mathbb{R}$, associated with the SMPC problem
is defined as the difference in volume of the deterministically tightened
state constraint set with and without the knowledge of $\mathbb{P}^{w}$ and
$\mathbb{P}_{k}^{\mathbf{x}}$ respectively. That is,

$\displaystyle\mathfrak{C}_{\mathrm{SMPC}}$
$\displaystyle:=\int_{\underbrace{\left(\mathcal{X}_{\mathrm{True}}\ominus\mathcal{X}_{\mathrm{DR}}\right)}_{:=\mathcal{X}_{\mathrm{Diff}}}}dx,\quad\text{where},$
(20) $\displaystyle\mathcal{X}_{\mathrm{True}}$
$\displaystyle:=\bigcap^{n_{x}}_{i=1}\left\\{\mathbf{x}\mid
f^{\top}_{i}{\mathbf{x}}\leq
g_{i}-\Phi^{-1}_{\mathbb{P}_{k}^{\mathbf{x}}}(1-\delta_{i})\left\lVert{\mathbf{\Sigma}_{\mathbf{x}}^{\frac{1}{2}}f_{i}}\right\rVert_{2}\right\\},$
(21) $\displaystyle\mathcal{X}_{\mathrm{DR}}$
$\displaystyle:=\bigcap^{n_{x}}_{i=1}\left\\{\mathbf{x}\mid
f^{\top}_{i}{\mathbf{x}}\leq
g_{i}-\sqrt{\frac{1-\delta_{i}}{\delta_{i}}}\left\lVert{\mathbf{\Sigma}_{\mathbf{x}}^{\frac{1}{2}}f_{i}}\right\rVert_{2}\right\\}.$
(22)

Note that $\mathcal{X}_{\mathrm{True}}$ and $\mathcal{X}_{\mathrm{DR}}$ can be
equivalently written as

$\displaystyle\mathcal{X}_{\mathrm{True}}$
$\displaystyle=\left\\{\mathbf{x}\mid\mathbf{F}{\mathbf{x}}\leq\mathbf{g}_{\mathrm{True}}\right\\},$
(23) $\displaystyle\mathcal{X}_{\mathrm{DR}}$
$\displaystyle=\left\\{\mathbf{x}\mid\mathbf{F}{\mathbf{x}}\leq\mathbf{g}_{\mathrm{DR}}\right\\},\quad\text{where}\quad\forall
i\in[1:n_{x}]$ (24)
$\displaystyle\begin{bmatrix}\mathbf{g}_{\mathrm{True}}\end{bmatrix}_{i}$
$\displaystyle:=\begin{bmatrix}g_{i}-\Phi^{-1}_{\mathbb{P}_{k}^{\mathbf{x}}}(1-\delta_{i})\left\lVert{\mathbf{\Sigma}_{\mathbf{x}}^{\frac{1}{2}}f_{i}}\right\rVert_{2}\end{bmatrix},\quad\text{and},$
(25) $\displaystyle\begin{bmatrix}\mathbf{g}_{\mathrm{DR}}\end{bmatrix}_{i}$
$\displaystyle:=\begin{bmatrix}g_{i}-\sqrt{\frac{1-\delta_{i}}{\delta_{i}}}\left\lVert{\mathbf{\Sigma}_{\mathbf{x}}^{\frac{1}{2}}f_{i}}\right\rVert_{2}\end{bmatrix}.$
(26)

Let
$\mathbf{z}:=\left[\begin{array}[]{c|c}\mathbf{x}^{\top}&\mathbf{x}^{\top}\end{array}\right]^{\top}\in\mathbb{R}^{2n(N+1)}$
and

$\displaystyle\mathbf{h}:=\left[\begin{array}[]{c}\mathbf{g}_{\mathrm{True}}\\\
\hline\cr\mathbf{g}_{\mathrm{DR}}\end{array}\right]\in\mathbb{R}^{2n_{x}},\quad\mathbf{H}:=\left[\begin{array}[]{c|c}\mathbf{F}&\mathbf{F}\\\
\hline\cr\mathbf{0}_{n_{x}\times n(N+1)}&\mathbf{F}\end{array}\right].$ (31)

Then, using Theorem 3.3 in Gabidullina (2019), we can express
$\mathcal{X}_{\mathrm{Diff}}$ as

$\displaystyle\mathcal{X}_{\mathrm{Diff}}:=\left\\{\mathbf{z}\mid\mathbf{H}{\mathbf{z}}\leq\mathbf{h}\right\\}.$
(32)

Although analyzing (32) is in general hard, its volume (i.e., conservatism)
can be computed numerically, see e.g. Chevallier et al. (2022).

Remark 1: Note that (20) uses the Pontryagin difference of the two
deterministically tightened convex polyhedrons $\mathcal{X}_{\mathrm{True}}$
and $\mathcal{X}_{\mathrm{DR}}$. The constraint tightening in both (21) and
(22) depends upon the risk222When a non-uniform risk allocation is performed
by solving an optimization problem as in Ono & Williams (2008), both
$\mathcal{X}_{\mathrm{Diff}}$ and $\mathfrak{C}_{\mathrm{SMPC}}$ will have
interesting observations. Analysing this aspect is out of the scope of this
paper and is left for future. $\delta_{i}$ corresponding to the $i$th
constraint. Only for $\delta_{i},i=1,\dots,n_{x}$ close to $0$, the
distributionally robust constraint tightening will be significantly stricter
than the exact tightening.

###  3.2 Regret of SMPC

We start by noting that for SMPC formulation (16), regret is introduced by
conservatism: Since the deterministic constraint tightening will be different
for the case when the true distributions $\mathbb{P}_{w}$ and
$\mathbb{P}_{\mathbf{x}}$ are known and unknown respectively, the resulting
optimal cost from solving the Surrogate SMPC problem given by (16) will be
different when tightened state constraints become active. This difference of
the optimal costs is essentially referred to as regret.

######  Assumption 1.

Both the systems under the distributionally robust and fully informed SMPC
respectively encounter the same disturbance realization $w_{k}$ at all time
steps $k\in\mathbb{N}$.

######  Definition 2.

Given Assumption 1, the distributional regret associated with the SMPC problem
at time step $k$ denoted by $\mathfrak{R}^{\mathrm{SMPC}}_{k}\in\mathbb{R}$ is
defined as the difference in the cost associated while solving the Surrogate
SMPC problem (16) with and without the knowledge of $\mathbb{P}^{w}$ and
$\mathbb{P}_{k}^{\mathbf{x}}$ respectively. That is,

$\displaystyle\mathfrak{R}^{\mathrm{SMPC}}_{k}:=J^{\dagger}_{\mathrm{SMPC}}(x_{k}^{\dagger})-J^{*}_{\mathrm{SMPC}}(x_{k}^{*}).$
(33)

To analyze the regret, we derive a closed-form expression for the optimal
input sequences $\mathbf{u}^{*}_{k}$ and $\mathbf{u}^{\dagger}_{k}$. We
reformulate OCP (16) as a quadratic program (QP) of the form

$\displaystyle\\!\underset{\mathbf{u_{k}}}{\mathrm{min}}$
$\displaystyle\frac{1}{2}\left\lVert{\mathbf{u}_{k}}\right\rVert^{2}_{H}+h_{k}^{\top}\mathbf{u}_{k}+r_{k}$
(34a) $\displaystyle\mathrm{s.}\,\mathrm{t.}$
$\displaystyle\mathbf{M}\mathbf{u}_{k}-\mathbf{b}_{k}\leq\mathbf{0}.$ (34b)

We start by reformulating the cost function (13) of OCP (16). Substituting (7)
and (8) into (13), we obtain

$\displaystyle
J_{\mathrm{SMPC}}(\mathbf{u}_{k},x_{k})=\frac{1}{2}\left\lVert{\mathbf{u}_{k}}\right\rVert^{2}_{H}+h_{k}^{\top}\mathbf{u}_{k}+r_{k},$
(35)

where $H=2(\mathbf{B}^{\top}\mathbf{Q}\mathbf{B}+\mathbf{R})\succ
0,h_{k}^{\top}=2x^{\top}_{0\mid k}\mathbf{A}^{\top}\mathbf{Q}\mathbf{B}$ and
$r_{k}=\mathbf{Tr}\left(\mathbf{Q}\mathbf{D}\mathbf{\Sigma}_{w}\mathbf{D}^{\top}\right)+\left\lVert{x_{0\mid
k}}\right\rVert^{2}_{\mathbf{A}^{\top}\mathbf{Q}\mathbf{A}}$. Next, we
substitute (7) and (8) into (16d) to reformulate the tightened state
constraints in terms of only the input as

$\displaystyle\underbrace{f_{i}^{\top}\mathbf{B}}_{=:\bar{f}_{i}^{\top}}\mathbf{u}_{k}\leq\underbrace{g_{i}-f^{\top}_{i}\mathbf{A}x_{0\mid
k}}_{=:\bar{g}_{k,i}}-\psi_{i}\underbrace{\left\lVert{(\mathbf{D}\mathbf{\Sigma}_{\mathbf{w}}\mathbf{D}^{\top})^{\frac{1}{2}}f_{i}}\right\rVert_{2}}_{=:v_{i}}.$
(36)

Then, $\forall i\in[1:n_{x}]$, writing (36) equivalently in vectorized form as
$\mathbf{\bar{F}}\mathbf{u}_{k}\leq\mathbf{\bar{g}}_{k}-\mathbf{\psi}^{\top}\mathbf{v}$
enables us to define (34b) via

$\displaystyle\mathbf{M}=\begin{bmatrix}\mathbf{C}\\\
\mathbf{\bar{F}}\end{bmatrix},\quad\mathbf{b}_{k}=\begin{bmatrix}\mathbf{d}\\\
\mathbf{\bar{g}}_{k}-\mathrm{diag}({\bm{\psi}})\mathbf{v}\end{bmatrix}.$ (37)

######  Definition 3.

An inequality constraint is said to be active if
$\mathbf{M}_{i:}\mathbf{u}_{k}-\mathbf{b}_{k,i}=0$ and inactive if
$\mathbf{M}_{i:}\mathbf{u}_{k}-\mathbf{b}_{k,i}<0$, where $\mathbf{M}_{i:}$
denotes the i${}^{\text{th}}$ row of $\mathbf{M}$ and $\mathbf{b}_{k,i}$ the
i${}^{\text{th}}$ entry of $\mathbf{b}_{k}$. The active set
$\mathcal{A}_{k}\subseteq\mathcal{I}:=[1:n_{x}]$ is the index set of active
inequality constraints.

######  Assumption 2.

The active set $\mathcal{A}_{k}^{\diamond}$ of QP (34) is known $\forall k$.

######  Assumption 3.

QP (34) satisfies the linear independence constraint qualification (LICQ)
criterion, i.e., the gradients of the active inequality constraints are
linearly independent.

######  Proposition 1.

Under Assumptions 2 and 3, the unique and global solution to QP (34) at time
step $k$ is given by

$\mathbf{u}_{k}^{\diamond}=\left(V_{k}^{\diamond}\mathbf{\tilde{M}}^{\diamond}_{k}H^{-1}-H^{-1}\right)h_{k}+V_{k}^{\diamond}\mathbf{\tilde{b}}_{k}^{\diamond},$
(38)

where
$V_{k}^{\diamond}=H^{-1}\mathbf{\tilde{M}}_{k}^{\diamond^{\top}}\left(\mathbf{\tilde{M}}_{k}^{\diamond}H^{-1}\mathbf{\tilde{M}}_{k}^{\diamond^{\top}}\right)^{-1}$,
$\mathbf{\tilde{M}_{k}^{\diamond}}=[\mathbf{M}_{i:}]_{i\in\mathcal{A}_{k}^{\diamond}}$,
and
$\mathbf{\tilde{b}}_{k}^{\diamond}=[\mathbf{b}_{k,i}]_{i\in\mathcal{A}^{\diamond}_{k}}$.

###### Proof.

The result directly follows from applying the method of Lagrange multipliers
Ghojogh et al. (2021); Boyd & Vandenberghe (2004) to QP (34), exploiting
sufficiency of the active set Arnström (2023). We first formulate the
Lagrangian as

$\displaystyle\mathcal{L}(\mathbf{u}_{k},\mathbf{\mu}_{k})=\frac{1}{2}\mathbf{u}_{k}^{\top}H\mathbf{u}_{k}\\!+\\!h_{k}^{\top}\mathbf{u}_{k}\\!+\\!r_{k}\\!+\\!\bm{\mu}_{k}^{\top}(\mathbf{M}\mathbf{u}_{k}\\!-\\!\mathbf{b}_{k}),$
(39)

where $\mathbf{\mu}_{k}$ is the vector of Lagrange multipliers. Applying the
Karush-Kuhn-Tucker (KKT) conditions yields

$\displaystyle H\mathbf{u}_{k}^{\diamond}+h_{k}+\mathbf{M}^{\top}\bm{\mu}_{k}$
$\displaystyle=\mathbf{0}$ (40a)
$\displaystyle\mathbf{M}\mathbf{u}_{k}^{\diamond}-\mathbf{b}_{k}$
$\displaystyle\leq\mathbf{0}$ (40b) $\displaystyle\bm{\mu}_{k}$
$\displaystyle\geq\mathbf{0}$ (40c)
$\displaystyle\mu_{k,i}(\mathbf{M}_{i:}\mathbf{u}_{k}^{\diamond}-\mathbf{b}_{k,i})$
$\displaystyle=0,\quad\forall i\in[0:n_{x}].$ (40d)

According to the complementary slackness condition (40d), we have
$\mu_{k,i}=0,\forall i\in\mathcal{I}\setminus\mathcal{A}_{k}^{\diamond}$. This
enables to remove inactive inequality constraints from the problem formulation
(c.f., sufficiency of the active set Arnström (2023)). Hence, the KKT system
(40) is reduced to a linear system of equations given by

$\begin{bmatrix}H&\mathbf{\tilde{M}}_{k}^{\diamond^{\top}}\\\
\mathbf{\tilde{M}}_{k}^{\diamond}&\mathbf{0}\end{bmatrix}\begin{bmatrix}\mathbf{u}_{k}^{\diamond}\\\
\bm{\tilde{\mu}}_{k}\end{bmatrix}=\begin{bmatrix}-h_{k}\\\
\mathbf{\tilde{b}}_{k}^{\diamond}\end{bmatrix},$ (41)

where
$\bm{\tilde{\mu}}_{k}=[\mathbf{\mu}_{k,i}]_{i\in\mathcal{A}_{k}^{\diamond}},\mathbf{\tilde{M}_{k}^{\diamond}}=[\mathbf{M}_{i:}]_{i\in\mathcal{A}_{k}^{\diamond}}$,
$\mathbf{\tilde{b}}_{k}^{\diamond}=[\mathbf{b}_{k,i}]_{i\in\mathcal{A}^{\diamond}_{k}}$
and the dual feasibility condition (40c) is trivially fulfilled if
$\mathcal{A}_{k}^{\diamond}$ is the correct active set. Exploiting
invertibility of the coefficient matrix of (41) yields (38). Expression (38)
is guaranteed to be the optimal solution to QP (34) as the KKT conditions are
necessary and sufficient under Assumption 3. The solution is unique and global
as $H\succ 0$. ∎

The optimal input sequences $\mathbf{u}^{*}_{k}$ and
$\mathbf{u}^{\dagger}_{k}$ are obtained using (38) (Proposition 1) when exact
and distributionally robust constraint tightening are employed, respectively.

Remark 2: The active set $\mathcal{A}_{k}^{\diamond}$ can be constructed in an
iterative manner through additions and removals of constraints to a working
set $\mathcal{A}_{k}$ until primal feasibility (40b) and dual feasibility
(40c) are satisfied by solution (38). This idea is exploited in so-called
active-set methods for solving QPs Arnström (2023). Hence, Assumption 2 is not
restrictive and only used to avoid computations that are out of the scope of
this work. As for Assumption 3, the LICQ criterion is commonly employed in
practice as it is a rather weak condition. In the considered controller set-
up, LICQ can be established during control design and is naturally given in
many cases. However, other constraint qualifications might be used, see e.g.
Bergmann & Herzog (2019).

Remark 3: Next, we investigate the regret as defined in (33). Given the
optimal input sequences $\mathbf{u}_{k}^{\star}$ and
$\mathbf{u}_{k}^{\dagger}$ according to Proposition 1 for appropriate
constraint tightening, it can be seen that regret is induced through three
main effects:

  1. (i)

The active sets $\mathcal{A}_{k}^{\star}$ and $\mathcal{A}_{k}^{\dagger}$
differ as the distributionally robust controller is more cautious and
constraints might become active or inactive at earlier or later times when
compared to the fully-informed controller.

  2. (ii)

The initial states $x_{0\mid k}^{\star}$ and $x_{0\mid k}^{\dagger}$ differ as
the feasible set of states of the distributionally robust controller is
smaller than that of the fully-informed controller.

  3. (iii)

The tightening factors (17) of the active state constraints differ between the
distributionally robust and the fully-informed controller.

While regret can be easily computed using the previously stated results,
analyzing its general behavior is challenging. Since the active sets enter the
respective control sequences in a highly nonlinear way, their influence on the
regret is hardly analyzable. Therefore, we focus in the following on a special
case defined in the following assumption.

######  Assumption 4.

The active sets $\mathcal{A}_{k}^{*}$ and $\mathcal{A}_{k}^{\dagger}$ are
identical at time step $k$, such that
$\mathbf{\tilde{M}}_{k}^{\star}=\mathbf{\tilde{M}}_{k}^{\dagger}=\mathbf{\tilde{M}}_{k}$,
$\mathbf{V}_{k}^{\star}=\mathbf{V}_{k}^{\dagger}=\mathbf{V}_{k}$ and the
difference between $\mathbf{\tilde{b}_{k}^{\star}}$ and
$\mathbf{\tilde{b}_{k}^{\dagger}}$ is only due to the different tightening
factors $\bm{\tilde{\psi}}^{\star}$ and $\bm{\tilde{\psi}}^{\dagger}$
corresponding to the active state constraints.

######  Theorem 1.

Let Assumption 4 hold and write
$\mathbf{V}_{k}=\begin{bmatrix}\mathbf{V}_{1,k}&\mathbf{V}_{2,k}\end{bmatrix}$
as block-matrix composed of $\mathbf{V}_{1,k}$ and $\mathbf{V}_{2,k}$. Then,
the regret at time step $k$ is given by

$\displaystyle\begin{split}\mathfrak{R}^{\mathrm{SMPC}}_{k}\\!=\\!&-(x_{0\mid
k}^{\star}\\!-\\!x_{0\mid k}^{\dagger})^{\top}\Lambda_{1,k}(x_{0\mid
k}^{\star}\\!+\\!x_{0\mid k}^{\dagger})\\\
&-\mathbf{\tilde{v}}_{k}^{\top}\mathrm{diag}(\bm{\tilde{\psi}}_{k}^{\star}\\!-\\!\bm{\tilde{\psi}}_{k}^{\dagger})\Lambda_{2,k}\mathrm{diag}(\bm{\tilde{\psi}}_{k}^{\star}\\!+\\!\bm{\tilde{\psi}}_{k}^{\dagger})\mathbf{\tilde{v}}_{k}\\\
&+(x_{0\mid k}^{\star}\\!-\\!x_{0\mid
k}^{\dagger})^{\top}\Lambda_{3,k}\mathrm{diag}(\bm{\tilde{\psi}}^{\star}_{k}\\!+\\!\bm{\tilde{\psi}}^{\dagger}_{k})\mathbf{\tilde{v}}_{k}\\\
&+(x_{0\mid k}^{\star}\\!+\\!x_{0\mid
k}^{\dagger})^{\top}\Lambda_{3,k}\mathrm{diag}(\bm{\tilde{\psi}}^{\star}_{k}\\!-\\!\bm{\tilde{\psi}}^{\dagger}_{k})\mathbf{\tilde{v}}_{k}\\\
&-(x_{0\mid k}^{\star}\\!-\\!x_{0\mid k}^{\dagger})^{\top}\Lambda_{4,k}\\\
&+\mathbf{\tilde{v}}^{\top}_{k}\mathrm{diag}(\bm{\tilde{\psi}}^{\star}_{k}\\!-\\!\bm{\tilde{\psi}}^{\dagger}_{k})\Lambda_{5,k},\end{split}$
(42)

where

$\displaystyle\Lambda_{1,k}=\frac{1}{2}\alpha_{k}^{\top}H\alpha_{k}+\frac{1}{2}(\tilde{h}^{\top}\alpha_{k}+\alpha_{k}^{\top}\tilde{h})+\mathbf{A}^{\top}\mathbf{QA},$
$\displaystyle\Lambda_{2,k}=\frac{1}{2}\mathbf{V}_{2,k}^{\top}H\mathbf{V}_{2,k},\quad\Lambda_{3,k}=\frac{1}{2}(\tilde{h}^{\top}\mathbf{V}_{2,k}+\alpha_{k}^{\top}H\mathbf{V}_{2,k})$
$\displaystyle\Lambda_{4,k}=\tilde{h}^{\top}\gamma_{k}+\alpha_{k}^{\top}H\gamma_{k},\quad\Lambda_{5,k}=\mathbf{V}_{2,k}^{\top}H\gamma_{k}$

with
$\alpha_{k}=(V_{k}\mathbf{\tilde{M}}_{k}H^{-1}-H^{-1})\tilde{h}-\mathbf{V}_{2,k}\mathbf{\tilde{F}}_{k}\mathbf{A}$,
$\tilde{h}=2\mathbf{B}^{\top}\mathbf{Q}\mathbf{A}$ and
$\gamma_{k}=\mathbf{V}_{1,k}\mathbf{\tilde{d}}_{k}+\mathbf{V}_{2,k}\mathbf{\tilde{g}}_{k}$.

###### Proof.

Expression (42) is obtained from (33) using the cost function representation
(35), substituting the optimal input sequences $\mathbf{u}^{*}_{k}$ and
$\mathbf{u}^{\dagger}_{k}$ from (38) and the definitions of
$\mathbf{\tilde{M}}_{k}$ and $\mathbf{\tilde{b}}_{k}$ and performing a series
of algebraic operations (not shown here for the brevity of presentation). ∎

Note that (42) is quadratic in the initial conditions $x_{0\mid k}^{\star}$,
$x_{0\mid k}^{\dagger}$ and the tightening factors
$\bm{\tilde{\psi}}_{k}^{\star}$, $\bm{\tilde{\psi}}_{k}^{\dagger}$ but not
directly in their respective differences.

######  Corollary 1.

If no (state and input) constraints are active, then (42) simplifies to

$\displaystyle\begin{split}\mathfrak{R}^{\mathrm{SMPC}}_{k}&=-(x_{0\mid
k}^{\star}\\!-\\!x_{0\mid k}^{\dagger})^{\top}\Lambda_{1}(x_{0\mid
k}^{\star}\\!+\\!x_{0\mid k}^{\dagger}),\\\
\Lambda_{1}&=\mathbf{A}^{\top}\mathbf{QA}-\frac{1}{2}\tilde{h}^{\top}H^{-1}\tilde{h}\end{split}$
(43)

and the optimal input sequences are given by the linear feedback laws
$\mathbf{u}_{k}^{\star}=-H^{-1}\tilde{h}x_{0\mid k}^{\star}$ and
$\mathbf{u}_{k}^{\dagger}=-H^{-1}\tilde{h}x_{0\mid k}^{\dagger}$.

We continue analyzing the general behavior of regret when no constraints are
active, given that the fully informed and the distributionally robust SMPC are
stabilizing in the sense of the following definitions.

######  Definition 4.

(From Culbertson et al. (2023)) A random variable $w\sim\mathbb{P}_{w}$
belongs to $L^{p}$, denoted by $w\\!\in\\!L^{p}$, for $p\\!>\\!0$ if $\lVert
w\rVert_{L^{p}}:=\mathbb{E}\left[\lVert\xi\rVert^{p}\right]^{\frac{1}{p}}<\infty$,
where $\lVert\cdot\rVert$ defines a typical norm on $\mathbb{R}^{n}$.

######  Definition 5.

(From Culbertson et al. (2023)) System (1) is input-to-state stable in
probability (ISSp) with respect to $L^{p}$ if $\forall\epsilon\\!\in\\!(0,1)$,
$K\\!\in\\!\mathbb{N}$ and $w_{k}\\!\in\\!L_{p}$,
$\exists\beta\\!\in\\!\mathcal{KL}$ and $\exists\varrho\\!\in\\!\mathcal{K}$
such that

$\mathbb{P}_{x}\left[\lVert x_{k+i}\rVert\leq\beta(\lVert
x_{k}\rVert,k\\!+\\!i)\\!+\\!\varrho(\lVert w_{k+i}\rVert_{L^{p}}),\forall
i\\!\leq\\!K\right]\\!\geq\\!1\\!-\\!\epsilon.$ (44)

If this holds for $\beta(\lVert x_{k}\rVert,k+i)=M\nu^{k+i}\lVert x_{k}\rVert$
with $M>0$ and $\nu\in(0,1)$, system (1) is exponentially ISSp (eISSp).

######  Corollary 2.

(From Culbertson et al. (2023)) The nominal (undisturbed) system is
asymptotically stable according to

$\mathbb{P}[\lVert\bar{x}_{k+i}\rVert\leq\beta(\lVert x_{k}\rVert,k+i),\forall
i\leq K)]=1.$ (45)

######  Definition 6.

(From Culbertson et al. (2023)) A continuous function
$V:\mathbb{R}^{n}\rightarrow\mathbb{R}_{\geq 0}$ is an ISSp Lyapunov function
for system (1) if there exist functions
$\kappa_{1},\kappa_{2},\kappa_{3}\\!\in\\!\mathcal{K}_{\infty}$ and
$\kappa_{4}\\!\in\\!\mathcal{K}$ such that

$\displaystyle\kappa_{1}(\lVert x_{k}\rVert)\leq V(x_{k})\leq\kappa_{2}(\lVert
x_{k}\rVert)$
$\displaystyle\mathbb{E}[V(f(x_{k},u_{k})\\!-\\!V(x_{k})]\leq\kappa_{3}(V(x_{k}))\\!+\\!\kappa_{4}(\lVert
w_{k}\rVert_{L^{p}})$

hold for all $x_{k}\in\mathcal{X}$ and $w_{k}\in L_{p}$. If $\kappa_{1}(\lVert
x_{k}\rVert)=a\lVert x_{k}\rVert^{c}$, $\kappa_{2}(\lVert x_{k}\rVert)=b\lVert
x_{k}\rVert^{c}$ and $\kappa_{3}(V(x))=\nu V(x)$ with
$a,b,c\in\mathbb{R}_{>0}$ and $\nu\in(0,1)$, then $V$ is an eISSp Lyapunov
function.

######  Assumption 5.

System (1) is eISSp under SMPC for both the exact and the distributionally
robust cases. That is, $\exists\beta^{*},\beta^{\dagger}\in\mathcal{KL}$,
given by
$\beta^{*}(s,t)=M^{*}\nu^{*^{t}}s,\beta^{\dagger}(s,t)=M^{\dagger}\nu^{\dagger^{t}}s$
with $M^{*},M^{\dagger}>0$ and $\nu^{*},\nu^{\dagger}\in(0,1)$ such that (44)
holds for $w_{k}\in L_{p},w_{k}\sim\mathbb{P}^{w}$ with $p>0$, given
$\varepsilon\\!\in\\!(0,1)$ and $\bar{K}=\min\\{K^{*},K^{\dagger}\\}$. The
optimal value functions
$J_{\mathrm{SMPC}}^{*}(x_{k}),J_{\mathrm{SMPC}}^{*}(x_{k})$ are eISSp Lyapunov
functions.

Note that there is an interplay between the probability $\varepsilon$ and the
horizon lengths $K^{*},K^{\dagger}$. For fixed $\varepsilon$, the (maximum)
horizon lengths $K^{*},K^{\dagger}$ for which (44) holds true in the
respective case is determined by the problem formulation and vice versa. We
consider the common horizon length $\bar{K}=\min\\{K^{*},K^{\dagger}\\}$ to
compare the fully informed and the distributionally robust case as (44)
remains valid for $\bar{K}$ in both cases.

######  Assumption 6.

There exists a set $\Phi:=\\{x\in\mathcal{X}\mid\lVert x\rVert\leq
r\\}\subset\mathcal{X}_{\mathrm{DR}}\subset\mathcal{X}_{\mathrm{True}}\subset\mathcal{X}$,
where $r\in\mathbb{R}_{>0}$ is chosen such that $\Phi$ is the largest set such
that no input or (tightened) state constraints are active for all $x\in\Phi$
for both the fully informed and the distributionally robust controller.

######  Definition 7.

(From Culbertson et al. (2023)) For a bounded set $\mathcal{S}\in\mathcal{X}$,
the hitting time is defined as
$\tau_{\mathcal{S}}(x):=\inf\\{k\in\mathbb{N}\mid
x_{k}\in\mathcal{S},x_{0}=x\\}$. The set $\mathcal{S}$ is recurrent if
$\forall x\in\mathcal{X},\mathbb{P}[\tau_{\mathcal{S}}(x)<\infty]=1$.

######  Lemma 3.

The set $\Phi$ is recurrent for system (1) under both the fully informed and
the distributionally robust controller. Furthermore, if $x_{k}\in\Phi$, then
system (1) satisfies

$\mathbb{P}[x_{k+i}\in\Phi,\forall i\leq\bar{K}]\geq 1-\varepsilon.$ (47)

###### Proof.

If $V$ is an eISSp Lyapunov functions for system (1), any sublevel set
$\mathcal{V}_{\gamma}\\!:=\\!\\{x\\!\in\\!\mathcal{X}\\!\mid\\!V(x)\\!\leq\\!\gamma\\}$
is recurrent, see Culbertson et al. (2023). As the optimal value functions are
eISSp Lyapunov functions for the respective closed-loop systems according to
Assumption 5 and $r>0$ by Assumption 6, there exist
$\gamma^{*},\gamma^{\dagger}$ such that
$\mathcal{V}_{\gamma^{*}}\subset\Phi,\mathcal{V}_{\gamma^{\dagger}}\subset\Phi$.
Hence, $\Phi$ is recurrent under both the fully informed and the
distributionally robust controller. Expression (47) follows from (44) as
$\beta(\lVert x_{k}\rVert,k\\!+\\!i)\\!+\\!\rho(\lVert
w_{k+i}\rVert_{L_{p}})\leq\beta(\lVert x_{k}\rVert,k)\\!+\\!\rho(\lVert
w_{k}\rVert_{L_{p}})=\lVert x_{k}\rVert\leq r$ for $x_{k}\\!\in\\!\Phi$. ∎

By Lemma 3, $\Phi$ takes the role of a probabilistic invariant set.

######  Lemma 4.

There exists $\tilde{\beta}\in\mathcal{KL}$ such that if
$x_{k}^{*},x_{k}^{\dagger}\in\Phi$,

$\mathbb{P}[\lVert x_{k+i}^{*}-x_{k+i}^{\dagger}\rVert\leq\tilde{\beta}(\lVert
x_{k}^{*}\rVert+\lVert x_{k}^{\dagger}\rVert,k+i),\forall i\leq\bar{K}]\geq
1-\varepsilon.$ (48)

###### Proof.

If $x_{k}^{*},x_{k}^{\dagger}\in\Phi$, both systems remain in $\Phi$ for
horizon $\bar{K}$ with probability at least $1-\varepsilon$ according to (47).
While $x_{k}^{*},x_{k}^{\dagger}\in\Phi$, the unconstrained error dynamics
between the fully informed and the distributionally robust closed-loop system
is given by
$e_{k+i}:=x_{k+i}^{*}-x_{k+i}^{\dagger}=\bar{x}_{k+i}^{*}-\bar{x}_{k+i}^{\dagger}$
as both systems encounter the same disturbance $w_{k+i}$ according to
Assumption 1. Hence, exploiting (45), we find that $\forall i\leq\bar{K}$

$\displaystyle\mathbb{P}[e_{k+i}\\!=\\!\bar{x}_{k+i}^{*}\\!-\\!\bar{x}_{k+i}^{\dagger}]\\!\geq\\!1\\!-\\!\varepsilon$
(49) $\displaystyle\Leftrightarrow~{}$ $\displaystyle\mathbb{P}[\lVert
e_{k+i}\rVert\\!\leq\\!\lVert x_{k+i}^{*}\rVert\\!+\\!\lVert
x_{k+i}^{\dagger}\rVert]\\!\geq\\!1\\!-\\!\varepsilon$ (50)
$\displaystyle\Leftrightarrow~{}$ $\displaystyle\mathbb{P}[\lVert
e_{k+i}\rVert\\!\leq\\!\beta^{*}(\lVert
x_{k}^{*}\rVert,k\\!+\\!i)\\!+\\!\beta^{\dagger}(\lVert
x_{k}^{\dagger}\rVert,k\\!+\\!i)]\\!\geq\\!1\\!-\\!\varepsilon$ (51)

with $\beta^{*},\beta^{\dagger}$ from Assumption 5. We can now construct
$\tilde{\beta}(\lVert x_{k}^{*}\rVert\\!+\\!\lVert
x_{k}^{\dagger}\rVert,k+i)=(M^{*}\\!+\\!M^{\dagger})(\nu^{*}\\!+\\!\nu^{\dagger})^{k+i}(\lVert
x_{k}^{*}\rVert\\!+\\!\lVert x_{k}^{\dagger}\rVert)\\!\geq\\!\beta^{*}(\lVert
x_{k}^{*}\rVert,k\\!+\\!i)\\!+\\!\beta^{\dagger}(\lVert
x_{k}^{\dagger}\rVert,k\\!+\\!i)$, where clearly
$\tilde{\beta}\in\mathcal{KL}$, which concludes the proof. ∎

######  Theorem 2.

If $x_{k}^{*},x_{k}^{\dagger}\in\Phi$, then $\exists\sigma\in\mathcal{KL}$
such that

$\mathbb{P}[\lVert\mathfrak{R}_{k+i}^{\mathrm{SMPC}}\rVert\\!\leq\\!\sigma(\lVert
x_{k}^{*}\rVert\\!+\\!\lVert x_{k}^{\dagger}\rVert,k\\!+\\!i),\forall
i\\!\leq\\!\bar{K}]\\!\geq\\!1\\!-\\!\varepsilon.$ (52)

###### Proof.

If $x_{k}^{*},x_{k}^{\dagger}\in\Phi$, both systems remain in $\Phi$ for
horizon $\bar{K}$ with probability at least $1-\varepsilon$ according to (47).
As no constraints are active for $x\in\Phi$, we find that $\forall
i\leq\bar{K}$

$\displaystyle\mathbb{P}[\mathfrak{R}_{k+i}^{\mathrm{SMPC}}=\eqref{eq:regret_unconstr}]\\!\geq\\!1\\!-\\!\varepsilon$
$\displaystyle\Leftrightarrow~{}$
$\displaystyle\mathbb{P}[\lVert\mathfrak{R}_{k+i}^{\mathrm{SMPC}}\rVert\\!\leq\\!\lVert
x_{k+i}^{*}\\!-\\!x_{k+i}^{\dagger}\rVert\lVert\Lambda_{1}\rVert\lVert
x_{k+i}^{*}\\!+\\!x_{k+i}^{\dagger}\rVert]\\!\geq\\!1\\!-\\!\varepsilon$
$\displaystyle\Leftrightarrow~{}$
$\displaystyle\mathbb{P}[\lVert\mathfrak{R}_{k+i}^{\mathrm{SMPC}}\rVert\\!\leq\\!\tilde{\beta}(\lVert
x_{k}^{*}\rVert\\!+\\!\lVert
x_{k}^{\dagger}\rVert,k\\!+\\!i)2\lVert\Lambda_{1}\rVert
r]\\!\geq\\!1\\!-\\!\varepsilon,$

where $\lVert x_{k+i}^{*}\\!+\\!x_{k+i}^{\dagger}\rVert\leq\lVert
x_{k}^{*}\rVert\\!+\\!\lVert x_{k}^{\dagger}\rVert\leq 2r$ as
$x_{k+i}^{*},x_{k+i}^{\dagger}\in\Phi$. Clearly, $2\lVert\Lambda_{1}\rVert
r>0$ and hence, $\sigma(\lVert x_{k}^{*}\rVert\\!+\\!\lVert
x_{k}^{\dagger}\rVert,k\\!+\\!i)=2\lVert\Lambda_{1}\rVert
r\tilde{\beta}(\lVert x_{k}^{*}\rVert\\!+\\!\lVert
x_{k}^{\dagger}\rVert,k\\!+\\!i)$ is a class $\mathcal{KL}$ function. ∎

Theorem 2 implies that the cumulative regret
${\mathfrak{cR}}_{k}^{\mathrm{SMPC}}=\sum_{i=0}^{k}\mathfrak{R}^{\mathrm{SMPC}}_{i}$
will have sublinear growth over time. Furthermore, since regret is only
observed when state constraints become active, those results indicate the
potential of dynamically setting the constraints by leveraging the insights
obtained from regret analysis.

##  4 Numerical Simulation

(a) Control inputs ${u}^{\star}_{k}$ & ${u}^{\dagger}_{k}$

(b) System State: $[x]_{1}$

(c) System State: $[x]_{2}$

Figure 1: Simulation results for the distributionally robust and the fully
informed case. The fully informed, distributionally robust, upper bound and
lower bound quantities for all the respective parameters in all three subplots
are plotted in blue, red, yellow and purple colors respectively. The green,
dashed line indicates the time $\tau_{\Phi}$ when both systems entered the set
$\Phi$. It is evident that the distributionally robust system makes cautious
decisions and incurs regret.

(a) Instantaneous Distributional Regret

(b) Cumulative Distributional Regret

Figure 2: Distributional regret associated with the distributionally robust &
the fully informed controller cases are shown here.

We consider a discretized double integrator with discretization time step
$dt=0.05$ given by

$\displaystyle A=\begin{bmatrix}1&dt\\\ 0&1\end{bmatrix},B=\begin{bmatrix}0\\\
dt\end{bmatrix},C=I_{2},D=0_{2\times 1}.$ (53)

The model was simulated for a total of $T=200$ time steps with a total risk
budget of $\Delta=0.1$. The prediction horizon was set to be $N=5$ time steps.
The process noise was sampled from a multivariate Laplacian distribution with
zero mean and covariance of $\mathbf{\Sigma}^{\mathbf{w}}={0.01}^{2}I_{2}$.
The state and control penalty matrices are $Q=\mathrm{diag}(1,0.1)$ and
$R=0.1$ respectively. The control and the state constraints are given by
$u\in[-20,2]$ and $x\in[-4,15]\times[-4,1.5]$. The tightening constants for
both the distributionally robust case and the fully informed case was obtained
using (17). Further, the regret was computed using (33) for all time
steps333The code is available under: https://github.com/mpfefferk/Regret-and-
Conservatism-of-SMPC..

The results of simulating the system given by (53) using both the fully
informed and the distributionally robust SMPC controller are shown in Figures
1 and 2. Specifically, the control inputs and the states are plotted in
Figures 1(a) and 1(b-c) respectively. Similarly, the instantaneous and the
cumulative distributional regret associated with the distributionally robust
SMPC controller are plotted in Figures 2(a) and 2(b) respectively. It is
evident that in general the distributionally robust system exhibits a cautious
decision making and that we observe regret according to the three effects
identified in Remark 3. Precisely, at $t=4.1s$, all constraints become
inactive for both the distributionally robust and the fully informed
controller, indicating that both systems have entered the (here implicitly
defined) set $\Phi$. The distributionally robust controller being cautious to
avoid constraint violation is at that time closer to the reference than the
fully informed controller as seen in Figures 1(b-c). Since its remaining
unconstrained cost-to-go is lower than its counterpart, this explains the
observed, negative instantaneous regret in Figure 2(a).

##  5 Conclusion

A distributionally robust SMPC formulation exploiting moment-based ambiguity
sets has been analysed in regard of conservatism and regret. We have performed
an analysis of both conservatism and regret associated with distributionally
robust SMPC for linear stochastic systems. We observed a sublinear growth of
regret over time when state constraints became active. These findings were
underlined by simulations. Future research will be dedicated towards applying
the established framework to extended SMPC formulations that may include
online estimation or learning of unknown quantities. Furthermore, the effect
of non-uniform risk allocation on conservatism and regret is to be further
investigated.

## Acknowledgment

This project has received funding from the European Research Council (ERC)
under the European Union’s Horizon 2020 research and innovation program under
grant agreement No 834142 (Scalable Control) and from the German Research
Foundation (Research Training Group 2297). V. Renganathan is a member of the
ELLIIT Strategic Research Area in Lund University.

## References

  * Arnström (2023) Daniel Arnström.  _Real-Time Certified MPC: Reliable Active-Set QP Solvers_.  PhD thesis, Division of Automatic Control, Department of Electrical Engineering, Linköping University (Sweden), 2023.
  * Bergmann & Herzog (2019) Ronny Bergmann and Roland Herzog.  Intrinsic Formulation of KKT Conditions and Constraint Qualifications on Smooth Manifolds.  _SIAM Journal on Optimization_ , 29(4):2423–2444, 2019.
  * Boyd & Vandenberghe (2004) S. Boyd and L. Vandenberghe.  _Convex Optimization_.  Cambridge University Press, 2004.
  * Calafiore & El Ghaoui (2006) G. C. Calafiore and L. El Ghaoui.  On distributionally robust chance-constrained linear programs.  _Journal of Optimization Theory and Applications_ , 130(1), December 2006.
  * Chen et al. (2023) Yitian Chen, Timothy L Molloy, Tyler Summers, and Iman Shames.  Regret analysis of online lqr control via trajectory prediction and tracking.  In _Proceedings of The 5th Annual Learning for Dynamics and Control Conference_ , volume 211 of _Proceedings of Machine Learning Research_ , pp. 248–258. PMLR, 15–16 Jun 2023.
  * Chevallier et al. (2022) Augustin Chevallier, Frédéric Cazals, and Paul Fearnhead.  Efficient computation of the volume of a polytope in high-dimensions using piecewise deterministic Markov processes.  In _25th International Conference on Artificial Intelligence and Statistics_ , 2022.
  * Culbertson et al. (2023) Preston Culbertson, Ryan K. Cosner, Maegan Tucker, and Aaron D. Ames.  Input-to-state stability in probability.  _arXiv:2304.14578_ , 2023.
  * Gabidullina (2019) ZR Gabidullina.  The Minkowski difference for convex polyhedra and some its applications.  _arXiv:1903.03590_ , 2019.
  * Ghojogh et al. (2021) Benyamin Ghojogh, Ali Ghodsi, Fakhri Karray, and Mark Crowley.  KKT Contitions, First-Order and Second-Order Optimization, and Distributed Optimization: Tutorial and Survey.  _arXiv:2110.01858_ , 2021.
  * Karapetyan et al. (2022) Aren Karapetyan, Andrea Iannelli, and John Lygeros.  On the regret of H∞ control.  In _2022 IEEE 61st Conference on Decision and Control (CDC)_ , pp. 6181–6186. IEEE, 2022.
  * Li et al. (2019) Yingying Li, Xin Chen, and Na Li.  Online Optimal Control with Linear Dynamics and Predictions: Algorithms and Regret Analysis.  In _33rd Conference on Neural Information Processing Systems_ , 2019\.
  * Ono & Williams (2008) Masahiro Ono and Brian C Williams.  Iterative risk allocation: A new approach to robust model predictive control with a joint chance constraint.  In _2008 47th IEEE Conference on Decision and Control_ , pp. 3427–3432. IEEE, 2008.
  * Paulson et al. (2020) Joel A Paulson, Edward A Buehler, Richard D Braatz, and Ali Mesbah.  Stochastic model predictive control with joint chance constraints.  _International Journal of Control_ , 93(1):126–139, 2020.
  * Rawlings et al. (2019) James B. Rawlings, David Q. Mayne, and Moritz M. Diehl.  _Model Predictive Control: Theory, Computation, and Design_.  Nob Hill Publishing, 2 edition, 2019.