$\mathsf{S}_{min}\subset\mathbb{Z}_{f}\oplus\mathsf{S}$ holds.

Proof. This follows directly from the proof of Theorem 2. If the minimal RPI
set $\mathsf{S}_{min}$ is contained in $\mathbb{Z}_{f}\oplus\mathsf{S}$ and
$\mathsf{S}_{min}$ is robustly asymptotically stable, the state reaches
$\mathbb{Z}_{f}\oplus\mathsf{S}$ in finite time steps.∎

In the proposed approach, a multi-stage MPC solution is computed on the
scenario tree for the large uncertainties with recourse, i.e. a tree of future
inputs depends on the realization of the uncertainty. The affine feedback is
added “on top” to robustify the solution against the small disturbances. The
feedback gain $\mathchoice{\mbox{$\displaystyle K$}}{\mbox{$\textstyle
K$}}{\mbox{$\scriptstyle K$}}{\mbox{$\scriptscriptstyle K$}}_{\mathrm{inv}}$
is fixed only for small disturbances and the degrees of freedom are increased
using the multi-stage approach for large uncertainties resulting in an
improved trade-off between optimality and complexity. Also, if the robust
horizon is chosen as $N_{r}\geq 1$, we have different feed-forward terms at
each stage in the predictions. This results in a scheme with the following
advantages when compared to multi-stage MPC and tube-based MPC independently:

  1. 1.

The growth in problem complexity is reduced when compared to a pure multi-
stage approach because the small uncertainties are not considered in the
scenario tree.

  2. 2.

The structurally relaxed recourse which is modeled in the prediction for the
realizations of the large uncertainties until robust horizon reduces the
conservatism compared to pure tube-based MPC.

  3. 3.

The choice of a robust horizon on the one hand limits the rapid growth of the
scenario tree and on the other hand, provides increased degrees of freedom
when compared to a standard tube-based MPC resulting in an improved trade-off.

  4. 4.

The use of low complexity tubes is possible for less conservatism because of
increased degrees of freedom in the form of feed-forward terms beyond the
robust horizon. This enables the application of the approach to high
dimensional systems.

##  6 Case study

The example considered in this paper is a continuous stirred-tank reactor
(CSTR) with a reaction scheme that is adapted from [31]. Two chemical
reactions take place in the reactor:

$\displaystyle\mathrm{A}\leavevmode\nobreak\ $
$\displaystyle\to\,\mathrm{B}\leavevmode\nobreak\ \to\,\mathrm{C}$
$\displaystyle 2\mathrm{A}\leavevmode\nobreak\ $
$\displaystyle\to\,\mathrm{D}$

The linearized discrete-time model has the form given in (1), where
$\mathchoice{\mbox{$\displaystyle x$}}{\mbox{$\textstyle
x$}}{\mbox{$\scriptstyle x$}}{\mbox{$\scriptscriptstyle x$}}=[\Delta
C_{\mathrm{a}},\,\Delta C_{\mathrm{b}},\,\Delta T_{\mathrm{R}},\Delta
T_{\mathrm{J}}]^{T}$ and $\mathchoice{\mbox{$\displaystyle
u$}}{\mbox{$\textstyle u$}}{\mbox{$\scriptstyle u$}}{\mbox{$\scriptscriptstyle
u$}}=\Delta F$. $\Delta C_{\mathrm{a}}$ and $\Delta C_{\mathrm{b}}$ denotes
the deviations of concentration of component A and B in mol/l ,
$\Delta{T}_{\mathrm{R}}$ and $\Delta{T}_{\mathrm{J}}$ denote the change in the
temperature of the reactor and in the jacket temperature in ∘C with respect to
the equilibrium point. The input is the deviation of the feed from the
equilibrium in l/h. The uncertain model is given as:

$\mathchoice{\mbox{$\displaystyle A$}}{\mbox{$\textstyle
A$}}{\mbox{$\scriptstyle A$}}{\mbox{$\scriptscriptstyle
A$}}_{\mathrm{unc}}(d_{1},d_{2},d_{3},d_{4})=\begin{pmatrix}0.3+d_{1}&-0.09&-0.01&0\\\
0.2&0.29+d_{2}&0.002&0\\\ d_{3}&d_{4}&1.10&0.15\\\
0.05&0.07&0.13&0.68\end{pmatrix}.$

The system vertex matrices are

$\displaystyle\mathchoice{\mbox{$\displaystyle A$}}{\mbox{$\textstyle
A$}}{\mbox{$\scriptstyle A$}}{\mbox{$\scriptscriptstyle A$}}_{1}$
$\displaystyle=\mathchoice{\mbox{$\displaystyle A$}}{\mbox{$\textstyle
A$}}{\mbox{$\scriptstyle A$}}{\mbox{$\scriptscriptstyle
A$}}_{\mathrm{unc}}(0.1,0.1,0.33,0.26),$
$\displaystyle\mathchoice{\mbox{$\displaystyle A$}}{\mbox{$\textstyle
A$}}{\mbox{$\scriptstyle A$}}{\mbox{$\scriptscriptstyle A$}}_{2}$
$\displaystyle=\mathchoice{\mbox{$\displaystyle A$}}{\mbox{$\textstyle
A$}}{\mbox{$\scriptstyle A$}}{\mbox{$\scriptscriptstyle
A$}}_{\mathrm{unc}}(-0.1,-0.1,-0.33,-0.26),$
$\displaystyle\mathchoice{\mbox{$\displaystyle A$}}{\mbox{$\textstyle
A$}}{\mbox{$\scriptstyle A$}}{\mbox{$\scriptscriptstyle A$}}_{3}$
$\displaystyle=\mathchoice{\mbox{$\displaystyle A$}}{\mbox{$\textstyle
A$}}{\mbox{$\scriptstyle A$}}{\mbox{$\scriptscriptstyle
A$}}_{\mathrm{unc}}(0.1,-0.1,0.33,-0.26),$
$\displaystyle\mathchoice{\mbox{$\displaystyle A$}}{\mbox{$\textstyle
A$}}{\mbox{$\scriptstyle A$}}{\mbox{$\scriptscriptstyle A$}}_{4}$
$\displaystyle=\mathchoice{\mbox{$\displaystyle A$}}{\mbox{$\textstyle
A$}}{\mbox{$\scriptstyle A$}}{\mbox{$\scriptscriptstyle
A$}}_{\mathrm{unc}}(-0.1,0.1,-0.33,0.26).$

The input vertex matrices are $\mathchoice{\mbox{$\displaystyle
B$}}{\mbox{$\textstyle B$}}{\mbox{$\scriptstyle B$}}{\mbox{$\scriptscriptstyle
B$}}_{1}=\mathchoice{\mbox{$\displaystyle B$}}{\mbox{$\textstyle
B$}}{\mbox{$\scriptstyle B$}}{\mbox{$\scriptscriptstyle
B$}}_{2}=\mathchoice{\mbox{$\displaystyle B$}}{\mbox{$\textstyle
B$}}{\mbox{$\scriptstyle B$}}{\mbox{$\scriptscriptstyle
B$}}_{3}=\mathchoice{\mbox{$\displaystyle B$}}{\mbox{$\textstyle
B$}}{\mbox{$\scriptstyle B$}}{\mbox{$\scriptscriptstyle
B$}}_{4}=[0.1,-0.05,0.8,0.1]^{T}$. The additive uncertainty bounds on all
states are $\mathchoice{\mbox{$\displaystyle w$}}{\mbox{$\textstyle
w$}}{\mbox{$\scriptstyle w$}}{\mbox{$\scriptscriptstyle
w$}}\in\mathchoice{\mbox{$\displaystyle\mathsf{W}$}}{\mbox{$\textstyle\mathsf{W}$}}{\mbox{$\scriptstyle\mathsf{W}$}}{\mbox{$\scriptscriptstyle\mathsf{W}$}}:=\\{\mathchoice{\mbox{$\displaystyle
w$}}{\mbox{$\textstyle w$}}{\mbox{$\scriptstyle w$}}{\mbox{$\scriptscriptstyle
w$}}|\,||w||_{\infty}\leq 0.1\\}$. The state constraints are
$\mathchoice{\mbox{$\displaystyle x$}}{\mbox{$\textstyle
x$}}{\mbox{$\scriptstyle x$}}{\mbox{$\scriptscriptstyle
x$}}\in\mathbb{X}:=\\{x|[-5,-5,-3,-5]^{T}\leq x\leq[5,5,3,5]^{T}\\}$ and the
input bounds are given by $\mathchoice{\mbox{$\displaystyle
u$}}{\mbox{$\textstyle u$}}{\mbox{$\scriptstyle u$}}{\mbox{$\scriptscriptstyle
u$}}\in\mathbb{U}:=\\{\mathchoice{\mbox{$\displaystyle u$}}{\mbox{$\textstyle
u$}}{\mbox{$\scriptstyle u$}}{\mbox{$\scriptscriptstyle u$}}|\,|u|\leq 2\\}$.
The control task is to take the system to a bounded set around the origin
while respecting the state and input constraints. The $\textstyle Q$ matrix of
the stage cost is chosen as an identity matrix and the $\textstyle R$ matrix
is chosen as $0.01$. The length of the prediction horizon of is chosen as
$N_{p}=5$.

###  6.1 Details of the simulation study

The tube-based part of the proposed scheme was implemented with three
different types of tubes and different simulation studies were performed. The
types of the tubes studied are given in Table 1.

Table 1: Details of the different types of tubes studied in this work. Tube type |  Implemented Optimization problem |  Tube complexity |  Recursive feasibility |  Stability
---|---|---|---|---
general complexity tube tube |  (18) |  18 inequalities and 44 vertices |  YES |  YES
Homothetic tube |  (23) |  18 inequalities and 44 vertices |  YES |  YES
Low complexity tube |  (23) with $8$ row $T$ matirx |  8 inequalities and 16 vertices |  YES |  YES

The additive disturbances are considered as small disturbances and an offline
invariant tube $\mathsf{S}$ was obtained. The set $\textstyle\mathsf{S}$ is
obtained by using (15) using the $\lambda-$contractive set for the value of
$\lambda=0.68$ obtained for the system (1) and it is contained in the box
given by
$\mathchoice{\mbox{$\displaystyle\mathsf{S}$}}{\mbox{$\textstyle\mathsf{S}$}}{\mbox{$\scriptstyle\mathsf{S}$}}{\mbox{$\scriptscriptstyle\mathsf{S}$}}\subseteq\\{x|[-0.4088,-0.5670,-0.3936,-0.3518]^{T}\leq
x\leq[0.4088,0.5670,0.3936,0.3518]^{T}\\}$ for the LQ-optimal feedback gain
$\mathchoice{\mbox{$\displaystyle K$}}{\mbox{$\textstyle
K$}}{\mbox{$\scriptstyle K$}}{\mbox{$\scriptscriptstyle
K$}}_{\mathrm{inv}}=\mathchoice{\mbox{$\displaystyle K$}}{\mbox{$\textstyle
K$}}{\mbox{$\scriptstyle K$}}{\mbox{$\scriptscriptstyle
K$}}_{\mathrm{pred}}=\mathchoice{\mbox{$\displaystyle K$}}{\mbox{$\textstyle
K$}}{\mbox{$\scriptstyle K$}}{\mbox{$\scriptscriptstyle
K$}}_{f}=K=[-0.0493,-0.0004,-1.3330,-0.3485]$. In Table 2, the volumes and the
computation times of the proposed approach for different robust horizons and
the tube-based MPC approach are given. In the following, we investigate all
aspects of the proposed approach.

Table 2: Quantitative comparison of feasible domains obtained using the proposed tube-enhanced multi-stage (TEMS) MPC scheme for varying robust horizons and the tube-based MPC. Tube type |  Para-meters |  Tube MPC (without the invariant tube) |  TEMS MPC with $N_{r}=0$ |  TEMS MPC with $N_{r}=1$ |  TEMS MPC with $N_{r}=2$ |  TEMS MPC with $N_{r}=3$ |  TEMS MPC with $N_{r}=4$ |  TEMS MPC with $N_{r}=5$
---|---|---|---|---|---|---|---|---
|  Volume |  1197.1 |  1110.7 |  4007.6 |  4392.7 |  4570.9 |  4574.6 |  4574.6
General complexity |  Comp. time [s] |  6.55 |  0.15 |  0.45 |  1.27 |  3.54 |  8.65 |  1.2
|  Volume |  1065.2 |  1001.0 |  3820.3 |  4319.5 |  4536.6 |  4574.2 |  as
Homo- thetic |  Comp. time [s] |  5.04 |  0.20 |  0.73 |  2.3 |  7.03 |  14.18 |  above
|  Volume |  96.02 |  96.02 |  1415.9 |  3563.2 |  4411.3 |  4545.9 |  as
Low complexity |  Comp. time [s] |  1.1 |  0.07 |  0.28 |  0.82 |  2.8 |  6.5 |  above

###  6.2 The effect of the invariant tube

If the pure multi-stage MPC is applied, it gives rise to 64 branches per node
in the scenario tree resulting in more than 100 million scenarios. The
additive disturbances therefore are removed from the multi-stage part and are
formulated in the (invariant) tube-based part of the scheme. Hence,
$\overline{\mathsf{W}}=\\{0\\}$ and $\underline{\mathsf{W}}=\mathsf{W}$ for
the studied example. This results in four branches at every node which is a
dramatic reduction when compared to the $64$ branches required in the case of
a full scenario tree for all the uncertainties.

For comparison purposes, tube-based MPC was implented for the three different
types of tubes considered in the proposed framework. The $\lambda-$contractive
set resulted in $32$ inequalities for the system (1) compared to $18$
inequalities for the system (3). The contractive sets were computed using the
multi-parametric toolbox [27]. Since the additive disturbances are not
considered in (3), the resulting complexity of the $\lambda-$contractive set
is different. The $3$D projections of the $\lambda-$ contractive sets and the
small disturbances invariant set $\mathsf{S}$ are shown in Figure 4.

(a)  Projection on $\Delta C_{a}$, $\Delta C_{b}$ and $\Delta T_{\mathrm{R}}$.
(b)  Projection on $\Delta C_{b}$, $\Delta T_{\mathrm{R}}$ and $\Delta
T_{\mathrm{J}}$.

Figure 4: Plot of $\lambda-$ contractive sets and the disturbance invariant
set $\mathsf{S}$

(a)  Projection on $\Delta C_{a}$, $\Delta C_{b}$ and $\Delta T_{\mathrm{R}}$.
(b)  Projection on $\Delta C_{b}$, $\Delta T_{\mathrm{R}}$ and $\Delta
T_{\mathrm{J}}$.

Figure 5: Comparison of the feasible domains obtained by the tube-based NMPC
(without the invariant tube $\mathsf{S}$) and the proposed scheme with
$N_{r}=0$ .

If the tube-based MPC [16] is employed without the classification of
uncertainties proposed in this paper, the tube at every time step is
characterized by $32$ inequalities in the case of homothetic and general
complexity tube tubes and the number of constraints increases with the number
of additive and parametric uncertainties. In the studied example, we must
consider $4$ vertex matrices and $16$ vertices of the additive disturbance
set. Hence, to formulate the propagation of tubes, $32\times 4\times 16=2048$
constraints are required per prediction step. In contrast, for the proposed
scheme with robust horizon $N_{r}=0$, only $18\times 4=72$ constraints are
required to formulate the propagation of the tubes. This is because the small
uncertainties are not considered both in the computation of contractive sets
and in the online problem. Instead, a suitable back-off is obtained by making
use of the disturbance invariant set $\mathsf{S}$. If the low complexity tube
is employed, the complexity of the tube-based MPC scheme without the invariant
tube is $16$ times larger than the proposed scheme with $N_{r}=0$ because of
the $16$ vertices of the additive disturbance bounds.

The feasible domains of the tube-based MPC scheme (without the invariant tube)
and the proposed tube-enhanced multi-stage MPC scheme with $N_{r}=0$ obtained
using the general complexity tube predicted tubes are shown in Figure 5. It
can be seen that the feasible region of the tube-based MPC scheme that does
not employ the invariant tube is larger than the feasible region of the
proposed scheme with $N_{r}=0$. The volume of the feasible domain of tube-
based MPC scheme (without the invariant tube) is 1197.1 whereas the proposed
scheme with $N_{r}=0$ results in a volume of 1110.7 (approx. 7% smaller).
However, the computation time of the proposed approach with $N_{r}=0$ is only
0.15 seconds compared to the other scheme which has a computation time of 6.55
seconds. If the homothetic tubes are employed, the feasible domain of the
proposed scheme is $6\%$ smaller while the computation time of the proposed
approach is approximately $96\%$ smaller (for exact values refer to Table 2).
In the case of low complexity tube, there is no difference in the volumes of
feasible domains observed.

With this example, it can be seen that the invariant set can introduce a
certain conservatism in the closed-loop. However, it reduces significantly the
computational complexity and improves the computation time of the approach
largely (up to 98% reduction). Thanks to the use an invariant tube for small
uncertainties, an important reduction in computation time can be expected at
the cost of only minor additional conservativeness.

###  6.3 The effect of the robust horizon

(a)  Projection on $\Delta C_{a}$, $\Delta C_{b}$ and $\Delta T_{\mathrm{R}}$.
(b)  Projection on $\Delta C_{b}$, $\Delta T_{\mathrm{R}}$ and $\Delta
T_{\mathrm{J}}$.

Figure 6: Comparison of the feasible domains obtained using the proposed
scheme that uses general complexity tubes for varying robust horizons.

(a)  Projection on $\Delta C_{a}$, $\Delta C_{b}$ and $\Delta T_{\mathrm{R}}$.
(b)  Projection on $\Delta C_{b}$, $\Delta T_{\mathrm{R}}$ and $\Delta
T_{\mathrm{J}}$.

Figure 7: Comparison of the feasible domains obtained using the proposed
scheme that uses homethetic tubes for varying robust horizons.

(a)  Projection on $\Delta C_{a}$, $\Delta C_{b}$ and $\Delta T_{\mathrm{R}}$.
(b)  Projection on $\Delta C_{b}$, $\Delta T_{\mathrm{R}}$ and $\Delta
T_{\mathrm{J}}$.

Figure 8: Comparison of the feasible domains obtained using the proposed
scheme that uses low complexity homethetic tubes for varying robust horizons.

If the proposed tube-enhanced multi-stage approach is employed to handle large
uncertainties, the number of scenarios considered by the optimization problem
increases with the length of robust horizon. Total number of scenarios in the
problem is determined by $N_{r}$. For $N_{r}=1$, the problem has only four
scenarios and for $N_{r}=5$, the problem has $4^{5}=1024$ scenarios. First, we
discuss the effect of the types of the predicted tubes that are employed and
then summarize the observations.

####  6.3.1 General complexity tube

The scheme with a robust horizon $N_{r}=4$ has the same feasible domain as the
scheme with full robust horizon. When the robust horizon decreases further,
the volume of the feasible domain decreases monotonically. The comparison of
the feasible domains of the proposed scheme with different robust horizons
that uses general complexity tubes are shown in Figure 6. The proposed scheme
with $N_{r}=1$ has a volume that is approximately 12% smaller than the volume
of the full horizon case. The volumes of the feasible domains of the tube-
based MPC implemented without the invariant tube and the proposed scheme with
$N_{r}=0$ are significantly smaller (approx. 74% smaller and 76% smaller). In
the tube-based MPC scheme, the optimization problem has one feed-forward term
per prediction step as degrees of freedom. Whereas, in the cased of $N_{r}=1$,
there are four feed-forward terms optimized at every stage. This improves the
degrees of freedom of the controller and results in an improved performance.
The polytopes are plotted with the help of the multi-parametric toolbox [27].

The computation times of the scheme with different robust horizons however
does not show a uniform trend. The computation times of robust horizons
$N_{r}=0$ and $N_{r}=1$ are smaller than the full robust horizon. However, the
schemes with the robust horizons 2 to 4 have computation times larger than the
scheme with full robust horizon. This is because of the difference in
complexity associated with the tree and the tube. The scheme with $N_{r}=2$
has $16$ scenarios. However, from the second prediction step, the propagation
of tubes is characterized by $72$ inequalities and this is formulated for all
$16$ scenarios requiring $72\times 16=1152$ inequalities after the second
prediction step. This leads to an increased computational effort when compared
to a full tree with $4\times 4=16$ equality constraints per node (though
exponentially increasing every stage).

####  6.3.2 Homothetic tube

The trends in the volumes of the feasible domain and the computation times are
similar to the general complexity tube case, if the homothetic tubes are
employed in the predictions. However, the scheme results in an increased
conservatism compared to the general complexity tube case. This is expected
because the shape of the predicted sets is restricted. The proposed scheme
with $N_{r}=1$ has a volume of the feasible domain that is smaller by
approximately $17\%$ in this case compared to the full robust horizon case.
The feasible domains of the proposed scheme implemented with homothetic tubes
are given in Figure 7. There is also an increase in computation times when
compared to the general complexity tube case for a fixed horizon. This is
because of the increase in the number of constraints due to $44$ vertices of
the tube considered to obtain the tight upper bound of the stage cost. This
leads to a proportional increase in computational cost.

####  6.3.3 Low complexity tube

The computation times are smaller in this case as expected because the tube is
represented using the minimal number of inequalities. There is however large
conservatism as a result. When $N_{r}=0$ is employed, the volume of the
feasible region is only $96.02$ which is more than $10$ times smaller than for
the homothetic tube case. Similar reductions in the feasible domains can be
observed across all robust horizons. The feasible domains of the proposed
scheme implemented with the Low complexity tube are given in Figure 8. An
interesting point to note here is that when $N_{r}=1$ is applied, the volume
of the feasible domain is much smaller for the same robust horizon than when
general complexity/homothetic tubes are employed. Despite this reduction, the
volume is larger than the best feasible domain observed for the tube-based
MPC. This clearly demonstrates the advantages of the proposed approach. For
robust horizon $N_{r}=2$, the feasible domain is larger and the computational
cost is still lower than tube-based MPC. Hence the scheme can be used to
control high dimensional systems, where the pure tube-based scheme is either
highly conservative (in the case of using a Low complexity tube) or
intractable (in the case of a higher complexity tube). The proposed scheme
offers an alternative with an improved trade-off between optimality and
complexity.

The computation time and the feasible domain for the case $N_{r}=N_{p}=5$ does
not depend on the type of tube employed because there is no tube-based part in
the predictions of the optimization problems (18) and (23) for this case.

####  6.3.4 Summary

Below are the summary of the observations:

  1. 1.

The volumes of the feasible domains of the proposed scheme increases
monotonically with the increase in the robust horizon. For a given robust
horizon, the volume of the feasible domain of the proposed scheme that employs
the Low complexity tube is lower and that of the general complexity tube is
higher. The volume of the scheme that employs a homothetic tube is in between
but it is closer to the general complexity tube than to the Low complexity
tube.

  2. 2.

The proposed scheme with $N_{r}=1$ has a larger volume of the feasible domain
and a smaller computational effort than the tube-based MPC that is employed
without the invariant tube. Even the proposed scheme with $N_{r}=1$ that
employs the Low complexity tube has a larger volume than the tube-based MPC
that employs a general complexity tube.

  3. 3.

Though the proposed scheme with Low complexity tube shows conservatism with
respect to the scheme with full robust horizon, it gives the best computation
times. In addition, the proposed scheme employed with Low complexity tubes
with $N_{r}\geq 1$ shows better performance than the tube-based MPC scheme
with complex/general complexity tubes.

  4. 4.

The scheme offers flexibility with respect to the choice of tube and robust
horizon. For example, the proposed scheme that employs a homothetic tube with
$N_{r}=1$ has a feasible volume and computation times comparable to that of
the proposed scheme with Low complexity tube with $N_{r}=2$. The scheme offers
a wide variety of options to achieve a desired trade-off between optimality
and complexity than the existing robust MPC schemes.

Figure 9: Trajectories of $\Delta C_{\mathrm{a}},\,\Delta
C_{\mathrm{b}},\,\Delta T_{\mathrm{R}},\,\Delta F$ obtained using 100
simulation runs.

For this example, $N_{r}=1$ using a homothetic/general complexity tube,
$N_{r}=2$ using Low complexity tube and $N_{r}=5$ are possible choices to
obtain good trade-offs overall in optimality and computational complexity. The
scheme with $N_{r}=5$ gives the best performance for a reasonable
computational time. Though the volume of the feasible domain of $N_{r}=1$
using general complexity tube is approximately 12% smaller, the computation
time is reduced by 62%. Hence $N_{r}=1$ using general complexity/homothetic
tube is a good choice for this example. The closed-loop state and input
trajectories of the proposed scheme employed using a general complexity tube
with $N_{r}=1$ for random initial conditions and random realizations of the
uncertainties for $100$ simulation runs is shown in Figure 9.

####  6.3.5 The effect of the prediction horizon

(a)  Y-axis in Linear scale  (b)  Y-axis in Logarithmic scale

Figure 10: Trend of computation times of the proposed scheme for varying
different robust horizons

The growth of the problem complexity with respect to the prediction horizon is
analyzed by comparing the average computation times of the schemes with
different robust horizons. The results are plotted in Figure 10. It can be
seen that the computational cost increases exponentially if a full robust
horizon is used and is significantly larger than for the other robust horizons
considered for $N_{p}=8$. The computation times of the schemes with
$N_{r}<N_{p}$ grows linearly in complexity with respect to the prediction
horizon. However, the slope is seen increasing when the robust horizon
increases. The proposed scheme with robust horizon $N_{r}=1$ has a computation
time of less than one second and the scheme with $N_{r}=2$ has a computation
time of less than $3$ seconds. The variations in the times are much smaller
than in the full robust horizon case. As the horizon grows larger, the
computational advantages of the proposed scheme increase.

##  7 Conclusion

In this paper, we have shown that the combination of multi-stage and tube-
based model predictive control schemes offers a flexible framework to manage
the trade-off between the performance and computational complexity of the
robust scheme. The proposed method uses a tube-based method to handle
uncertainties that are small or occur far in the prediction (after the robust
horizon), while the multi-stage approach handles the significant and immediate
uncertainties to increase the performance. The stability of the scheme for
linear systems with parametric and additive disturbances was demonstrated for
any choice of robust horizon, including the pure multi-stage case. Simulation
results show that the proposed method provides flexibility to obtain a good
trade-off between complexity and performance.

## References

  * Campo and Morari [1987] Campo PJ, Morari M.  Robust model predictive control.  In: Proc. of the American Control Conference; 1987. p. 1021–1026.
  * Lee and Yu [1997] Lee JH, Yu ZH.  Worst-case formulations of model predictive control for systems with bounded parameters.  Automatica 1997;33(5):763–781.
  * Scokaert and Mayne [1998] Scokaert POM, Mayne DQ.  Min-max feedback model predictive control for constrained linear systems.  IEEE Transactions on Automatic Control 1998;43(8):1136–1142.
  * Bernardini and Bemporad [2009] Bernardini D, Bemporad A.  Scenario-based model predictive control of stochastic constrained linear systems.  In: Proc. of the 48th IEEE Conference on Decision and Control, 2009.; 2009\. p. 6333–6338.
  * Lucia et al. [2013] Lucia S, Finkler T, Engell S.  Multi-stage Nonlinear Model Predictive Control Applied to a Semi-batch Polymerization Reactor under Uncertainty.  Journal of Process Control 2013;23:1306–1319.
  * Kothare et al. [1996] Kothare MV, Balakrishnan V, Morari M.  Robust constrained model predictive control using linear matrix inequalities.  Automatica 1996;32(10):1361 – 1379.
  * Lee and Kouvaritakis [2000] Lee YI, Kouvaritakis B.  Robust receding horizon predictive control for systems with uncertain dynamics and input saturation.  Automatica 2000;36(10):1497–1504.
  * Chisci et al. [2001] Chisci L, Rossiter JA, Zappa G.  Systems with persistent disturbances: predictive control with restricted constraints.  Automatica 2001;37(7):1019 – 1028.
  * Löfberg [2003] Löfberg J.  Approximations of closed-loop minimax MPC.  In: Proc. of the 42nd IEEE Conference on Decision and Control; 2003. p. 1438–1442.
  * Langson et al. [2004] Langson W, Chryssochoos I, Raković SV, Mayne DQ.  Robust model predictive control using tubes.  Automatica 2004;40(1):125 – 133.
  * Mayne et al. [2005] Mayne DQ, Seron MM, Rakovic SV.  Robust model predictive control of constrained linear systems with bounded disturbances.  Automatica 2005;41:219 – 224.
  * Rakovic et al. [2012] Rakovic SV, Kouvaritakis B, Cannon M, Panos C, Findeisen R.  Parameterized tube model predictive control.  IEEE Transactions on Automatic Control 2012;57(11):2746–2761.
  * Raković et al. [2012] Raković SV, Kouvaritakis B, Findeisen R, Cannon M.  Homothetic tube model predictive control.  Automatica 2012;48(8):1631–1638.
  * Raković et al. [2016] Raković SV, Levine WS, Açıkmeşe B.  Elastic tube model predictive control.  In: American Control Conference (ACC), 2016 IEEE; 2016. p. 3594–3599.
  * Villanueva et al. [2017] Villanueva ME, Quirynen R, Diehl M, Chachuat B, Houska B.  Robust MPC via min–max differential inequalities.  Automatica 2017;77:311 – 321.
  * Fleming et al. [2015] Fleming J, Kouvaritakis B, Cannon M.  Robust Tube MPC for Linear Systems With Multiplicative Uncertainty.  IEEE Transactions on Automatic Control 2015;60(4):1087–1092.
  * Muñoz-Carpintero et al. [2015] Muñoz-Carpintero D, Cannon M, Kouvaritakis B.  Robust MPC strategy with optimized polytopic dynamics for linear systems with additive and multiplicative uncertainty.  Systems & Control Letters 2015;81(Supplement C):34 – 41.
  * Lee and Kouvaritakis [2000] Lee YI, Kouvaritakis B.  A linear programming approach to constrained robust predictive control.  IEEE Transactions on Automatic Control 2000 Sep;45(9):1765–1770.
  * Blanco et al. [2010] Blanco TB, Cannon M, Moor BD.  On efficient computation of low-complexity controlled invariant sets for uncertain linear systems.  International Journal of Control 2010;83(7):1339–1346.
  * Subramanian et al. [2018] Subramanian S, Lucia S, Birjandi SAB, Paulen R, Engell S.  A Combined Multi-Stage and Tube-Based MPC Scheme for Constrained Linear Systems.  6th IFAC Conference on Nonlinear Model Predictive Control 2018;p. 577–582.
  * Rawlings and Mayne [2009] Rawlings JB, Mayne DQ.  Model Predictive Control Theory and Design.  Nob Hill Pub; 2009.
  * Blanchini and Miani [2008] Blanchini F, Miani S.  Set-theoretic methods in control.  Springer; 2008.
  * Hadjiyiannis et al. [2011] Hadjiyiannis MJ, Goulart PJ, Kuhn D.  An Efficient Method to Estimate the Suboptimality of Affine Controllers.  IEEE Transactions on Automatic Control 2011;56:2841– 2853.
  * Kouramas et al. [2005] Kouramas KI, Rakovic SV, Kerrigan EC, Allwright JC, Mayne DQ.  On the Minimal Robust Positively Invariant Set for Linear Difference Inclusions.  In: Proceedings of the 44th IEEE Conference on Decision and Control; 2005\. p. 2296–2301.
  * Lucia [2014] Lucia S.  Robust Multi-stage Nonlinear Model Predictive Control.  Shaker; 2014.
  * Kerrigan and Maciejowski [2004] Kerrigan EC, Maciejowski JM.  Feedback min-max model predictive control using a single linear program: robust stability and the explicit solution.  International Journal of Robust and Nonlinear Control 2004;14(4):395–413.
  * Herceg et al. [2013] Herceg M, Kvasnica M, Jones CN, Morari M.  Multi-Parametric Toolbox 3.0.  In: Proc. of the European Control Conference Zürich, Switzerland; 2013\. p. 502–510.  http://control.ee.ethz.ch/~mpt.
  * Scokaert and Mayne [1998] Scokaert POM, Mayne DQ.  Min-max feedback model predictive control for constrained linear systems.  IEEE Transactions on Automatic Control 1998;43(8):1136–1142.
  * Lu and Cannon [2019] Lu X, Cannon M.  Robust adaptive tube model predictive control.  In: 2019 American Control Conference (ACC) IEEE; 2019. p. 3695–3701.
  * Kouvaritakis and Cannon [2015] Kouvaritakis B, Cannon M.  Model Predictive Control: Classical, Robust and Stochastic.  Advanced Textbooks in Control and Signal Processing, Springer International Publishing; 2015.
  * Klatt and Engell [1998] Klatt KU, Engell S.  Gain-scheduling trajectory control of a continuous stirred tank reactor.  Computers & Chemical Engineering 1998;22:491–502.


|  |  | 0.00 | 0.00 | 0.00 | 0.04 | 0.46 | 0.00
|  |  | 0.00 | 0.00 | 0.00 | 0.04 | 0.46 | 0.00
|  | 0 | 0.00 | 0.00 | 0.00 | 0.07 | 0.08 | 0.00
|  |  | 0.00 | 0.00 | 0.00 | 0.07 | 0.08 | 0.00
|  |  | 0.00 | 0.00 | 0.00 | 0.07 | 0.08 | 0.00
|  | 1 | 0.00 | 0.00 | 0.00 | 0.41 | 0.50 | 0.00
|  |  | 0.00 | 0.00 | 0.00 | 0.41 | 0.50 | 0.00
|  |  | 0.00 | 0.00 | 0.00 | 0.41 | 0.50 | 0.00
| 1000 | -1 | 0.00 | 0.00 | 0.00 | 0.00 |  |
|  |  | 0.00 | 0.00 | 0.00 | 0.00 |  |
|  |  | 0.00 | 0.00 | 0.00 | 0.00 |  |
|  | 0 | 0.00 | 0.00 | 0.00 | 0.00 |  |
|  |  | 0.00 | 0.00 | 0.00 | 0.00 |  |
|  |  | 0.00 | 0.00 | 0.00 | 0.00 |  |
|  | 1 | 0.00 | 0.00 | 0.00 | 0.03 |  |
|  |  | 0.00 | 0.00 | 0.00 | 0.03 |  |
|  |  | 0.00 | 0.00 | 0.00 | 0.03 |  |
| 5000 | -1 | 0.00 | 0.00 | 0.00 | 0.00 |  |
|  |  | 0.00 | 0.00 | 0.00 | 0.00 |  |
|  |  | 0.00 | 0.00 | 0.00 | 0.00 |  |
|  | 0 | 0.00 | 0.00 | 0.00 | 0.00 |  |
|  |  | 0.00 | 0.00 | 0.00 | 0.00 |  |
|  |  | 0.00 | 0.00 | 0.00 | 0.00 |  |
|  | 1 | 0.00 | 0.00 | 0.00 | 0.00 |  |
|  |  | 0.00 | 0.00 | 0.00 | 0.00 |  |
|  |  | 0.00 | 0.00 | 0.00 | 0.00 |  |

Table 10: Data aggregation by treatment arm and outcome status when generating
one dataset using the data generating algorithm in Table 2 for small sample
sizes ($n=200$) under low event rate.

size | conf | treat_prev | assigned_treatment | outcome | followup_time | 0 | 1 | 2 | 3 | 4
---|---|---|---|---|---|---|---|---|---
200 | 0.1 | -1 | 0 | 0 | 134 | 86 | 58 | 35 | 24
|  |  |  | 1 | 8 | 7 | 2 | 1 | 0
|  |  | 1 | 0 | 57 | 19 | 4 | 1 | 0
|  |  |  | 1 | 1 | 1 | 1 | 0 | 0
|  | 0 | 0 | 0 | 87 | 42 | 19 | 8 | 3
|  |  |  | 1 | 6 | 1 | 2 | 1 | 1
|  |  | 1 | 0 | 102 | 56 | 27 | 16 | 9
|  |  |  | 1 | 5 | 2 | 1 | 1 | 0
|  | 1 | 0 | 0 | 51 | 16 | 6 | 1 | 1
|  |  |  | 1 | 0 | 0 | 0 | 0 | 0
|  |  | 1 | 0 | 140 | 103 | 73 | 50 | 33
|  |  |  | 1 | 9 | 2 | 3 | 1 | 5
| 0.5 | -1 | 0 | 0 | 146 | 101 | 69 | 45 | 29
|  |  |  | 1 | 9 | 7 | 3 | 3 | 1
|  |  | 1 | 0 | 44 | 16 | 2 | 0 | 0
|  |  |  | 1 | 1 | 1 | 0 | 0 | 0
|  | 0 | 0 | 0 | 97 | 40 | 27 | 10 | 4
|  |  |  | 1 | 7 | 3 | 0 | 1 | 0
|  |  | 1 | 0 | 95 | 48 | 24 | 12 | 7
|  |  |  | 1 | 1 | 2 | 1 | 2 | 2
|  | 1 | 0 | 0 | 58 | 12 | 3 | 0 | 0
|  |  |  | 1 | 3 | 1 | 0 | 0 | 0
|  |  | 1 | 0 | 135 | 90 | 61 | 39 | 28
|  |  |  | 1 | 4 | 3 | 5 | 0 | 1
| 0.9 | -1 | 0 | 0 | 123 | 84 | 60 | 38 | 22
|  |  |  | 1 | 11 | 4 | 4 | 4 | 3
|  |  | 1 | 0 | 62 | 18 | 6 | 2 | 0
|  |  |  | 1 | 4 | 1 | 0 | 0 | 0
|  | 0 | 0 | 0 | 92 | 45 | 27 | 11 | 3
|  |  |  | 1 | 5 | 1 | 0 | 1 | 1
|  |  | 1 | 0 | 97 | 51 | 22 | 10 | 3
|  |  |  | 1 | 6 | 1 | 0 | 0 | 0
|  | 1 | 0 | 0 | 59 | 16 | 7 | 3 | 0
|  |  |  | 1 | 5 | 0 | 0 | 0 | 0
|  |  | 1 | 0 | 127 | 71 | 48 | 34 | 16
|  |  |  | 1 | 9 | 2 | 3 | 0 | 2

Table 11: Data aggregation by treatment arm and outcome status when generating
one dataset using the data generating algorithm in Table 2 for medium sample
sizes ($n=1000$) under low event rate.

size | conf | treat_prev | assigned_treatment | outcome | followup_time | 0 | 1 | 2 | 3 | 4
---|---|---|---|---|---|---|---|---|---
1000 | 0.1 | -1 | 0 | 0 | 664 | 472 | 311 | 212 | 151
|  |  |  | 1 | 59 | 32 | 18 | 13 | 4
|  |  | 1 | 0 | 263 | 57 | 15 | 2 | 0
|  |  |  | 1 | 14 | 3 | 0 | 0 | 0
|  | 0 | 0 | 0 | 492 | 235 | 111 | 55 | 32
|  |  |  | 1 | 25 | 22 | 4 | 2 | 1
|  |  | 1 | 0 | 453 | 221 | 106 | 47 | 21
|  |  |  | 1 | 30 | 13 | 5 | 2 | 0
|  | 1 | 0 | 0 | 246 | 76 | 23 | 7 | 2
|  |  |  | 1 | 21 | 2 | 0 | 0 | 0
|  |  | 1 | 0 | 691 | 470 | 342 | 242 | 159
|  |  |  | 1 | 42 | 28 | 6 | 9 | 5
| 0.5 | -1 | 0 | 0 | 662 | 456 | 317 | 223 | 154
|  |  |  | 1 | 51 | 21 | 11 | 11 | 3
|  |  | 1 | 0 | 269 | 68 | 13 | 1 | 0
|  |  |  | 1 | 18 | 1 | 1 | 0 | 0
|  | 0 | 0 | 0 | 474 | 233 | 111 | 56 | 28
|  |  |  | 1 | 26 | 15 | 4 | 6 | 1
|  |  | 1 | 0 | 478 | 216 | 88 | 44 | 20
|  |  |  | 1 | 22 | 13 | 6 | 5 | 1
|  | 1 | 0 | 0 | 277 | 74 | 24 | 7 | 1
|  |  |  | 1 | 15 | 6 | 1 | 0 | 2
|  |  | 1 | 0 | 679 | 457 | 310 | 212 | 149
|  |  |  | 1 | 29 | 19 | 10 | 12 | 5
| 0.9 | -1 | 0 | 0 | 626 | 433 | 316 | 206 | 149
|  |  |  | 1 | 47 | 23 | 17 | 13 | 1
|  |  | 1 | 0 | 291 | 69 | 22 | 8 | 3
|  |  |  | 1 | 36 | 4 | 1 | 0 | 0
|  | 0 | 0 | 0 | 449 | 202 | 106 | 54 | 23
|  |  |  | 1 | 29 | 14 | 5 | 2 | 0
|  |  | 1 | 0 | 486 | 197 | 90 | 37 | 16
|  |  |  | 1 | 36 | 17 | 5 | 3 | 1
|  | 1 | 0 | 0 | 257 | 74 | 21 | 7 | 2
|  |  |  | 1 | 17 | 5 | 2 | 1 | 0
|  |  | 1 | 0 | 673 | 415 | 273 | 163 | 99
|  |  |  | 1 | 53 | 19 | 7 | 10 | 6

Table 12: Data aggregation by treatment arm and outcome status when generating
one dataset using the data generating algorithm in Table 2 for large sample
sizes ($n=5000$) under low event rate.

size | conf | treat_prev | assigned_treatment | outcome | followup_time | 0 | 1 | 2 | 3 | 4
---|---|---|---|---|---|---|---|---|---
5000 | 0.1 | -1 | 0 | 0 | 3374 | 2352 | 1639 | 1170 | 811
|  |  |  | 1 | 255 | 145 | 84 | 65 | 41
|  |  | 1 | 0 | 1308 | 329 | 91 | 32 | 15
|  |  |  | 1 | 63 | 19 | 7 | 4 | 0
|  | 0 | 0 | 0 | 2327 | 1101 | 507 | 246 | 120
|  |  |  | 1 | 157 | 68 | 24 | 12 | 7
|  |  | 1 | 0 | 2373 | 1126 | 549 | 262 | 131
|  |  |  | 1 | 143 | 70 | 20 | 14 | 3
|  | 1 | 0 | 0 | 1281 | 322 | 92 | 23 | 6
|  |  |  | 1 | 82 | 19 | 4 | 0 | 1
|  |  | 1 | 0 | 3461 | 2398 | 1673 | 1124 | 794
|  |  |  | 1 | 176 | 125 | 95 | 49 | 27
| 0.5 | -1 | 0 | 0 | 3356 | 2254 | 1599 | 1130 | 801
|  |  |  | 1 | 260 | 153 | 76 | 42 | 33
|  |  | 1 | 0 | 1290 | 347 | 93 | 22 | 3
|  |  |  | 1 | 94 | 19 | 9 | 3 | 0
|  | 0 | 0 | 0 | 2314 | 1079 | 545 | 263 | 140
|  |  |  | 1 | 150 | 73 | 25 | 12 | 5
|  |  | 1 | 0 | 2411 | 1108 | 512 | 227 | 97
|  |  |  | 1 | 125 | 57 | 22 | 11 | 8
|  | 1 | 0 | 0 | 1324 | 345 | 95 | 29 | 6
|  |  |  | 1 | 67 | 21 | 7 | 1 | 0
|  |  | 1 | 0 | 3421 | 2317 | 1551 | 1030 | 688
|  |  |  | 1 | 188 | 95 | 79 | 40 | 35
| 0.9 | -1 | 0 | 0 | 3243 | 2118 | 1413 | 982 | 644
|  |  |  | 1 | 235 | 153 | 87 | 46 | 41
|  |  | 1 | 0 | 1386 | 325 | 75 | 18 | 1
|  |  |  | 1 | 136 | 26 | 7 | 3 | 0
|  | 0 | 0 | 0 | 2369 | 1142 | 571 | 266 | 136
|  |  |  | 1 | 161 | 67 | 36 | 16 | 6
|  |  | 1 | 0 | 2267 | 959 | 406 | 152 | 69
|  |  |  | 1 | 203 | 61 | 36 | 12 | 5
|  | 1 | 0 | 0 | 1453 | 405 | 123 | 40 | 11
|  |  |  | 1 | 70 | 24 | 5 | 1 | 0
|  |  | 1 | 0 | 3227 | 2045 | 1256 | 757 | 477
|  |  |  | 1 | 250 | 112 | 75 | 51 | 28

Table 13: Summary statistics of IPTWs distribution in various simulation scenarios when $n=200$ Outcome event rate | Treatment prevalence | Confounding strength | Mean | Standard deviation | Max
---|---|---|---|---|---
Low | -1 | 0.10 | 0.99 | 0.11 | 1.55
Low | -1 | 0.50 | 1 | 0.28 | 5.58
Low | -1 | 0.90 | 1 | 0.35 | 4.36
Low | 0.00 | 0.10 | 1 | 0.18 | 3.35
Low | 0.00 | 0.50 | 0.99 | 0.16 | 2.07
Low | 0.00 | 0.90 | 0.98 | 0.47 | 7.23
Low | 1 | 0.10 | 1 | 0.12 | 1.57
Low | 1 | 0.50 | 1 | 0.28 | 3.26
Low | 1 | 0.90 | 1 | 0.32 | 6.14
Medium | -1 | 0.10 | 1 | 0.05 | 1.26
Medium | -1 | 0.50 | 1 | 0.20 | 3.82
Medium | -1 | 0.90 | 1 | 0.26 | 4.14
Medium | 0.00 | 0.10 | 1 | 0.07 | 1.69
Medium | 0.00 | 0.50 | 0.99 | 0.30 | 4.33
Medium | 0.00 | 0.90 | 1 | 0.51 | 9.54
Medium | 1 | 0.10 | 1.01 | 0.17 | 2.17
Medium | 1 | 0.50 | 1.03 | 0.64 | 11.43
Medium | 1 | 0.90 | 1 | 0.34 | 5.06
High | -1 | 0.10 | 1 | 0.08 | 1.91
High | -1 | 0.50 | 1 | 0.22 | 2.79
High | -1 | 0.90 | 1 | 0.30 | 4.94
High | 0.00 | 0.10 | 1.01 | 0.19 | 3.80
High | 0.00 | 0.50 | 1 | 0.38 | 5.60
High | 0.00 | 0.90 | 0.95 | 0.34 | 5.89
High | 1 | 0.10 | 1 | 0.08 | 1.59
High | 1 | 0.50 | 0.99 | 0.17 | 2.29
High | 1 | 0.90 | 0.99 | 0.36 | 4.77
Table 14: Summary statistics of inverse probability of IPTWs distribution in various simulation scenarios when $n=1000$ Outcome_prevalence | Treatment_prevalence | Confounding | Mean | standard deviation | Max
---|---|---|---|---|---
Low | -1 | 0.10 | 1 | 0.12 | 4.79
Low | -1 | 0.50 | 1 | 0.21 | 4.80
Low | -1 | 0.90 | 1 | 0.36 | 11.74
Low | 0.00 | 0.10 | 1 | 0.13 | 3.36
Low | 0.00 | 0.50 | 0.99 | 0.24 | 4.21
Low | 0.00 | 0.90 | 0.98 | 0.42 | 7.85
Low | 1 | 0.10 | 1 | 0.10 | 2.22
Low | 1 | 0.50 | 1.01 | 0.22 | 3.73
Low | 1 | 0.90 | 1.01 | 0.50 | 14.10
Medium | -1 | 0.10 | 0.99 | 0.10 | 1.82
Medium | -1 | 0.50 | 1 | 0.20 | 3.45
Medium | -1 | 0.90 | 0.99 | 0.29 | 6.51
Medium | 0.00 | 0.10 | 0.99 | 0.14 | 2.25
Medium | 0.00 | 0.50 | 1 | 0.25 | 3.53
Medium | 0.00 | 0.90 | 1.02 | 0.58 | 15.86
Medium | 1 | 0.10 | 1 | 0.09 | 1.54
Medium | 1 | 0.50 | 1 | 0.22 | 4.43
Medium | 1 | 0.90 | 0.99 | 0.35 | 8.52
High | -1 | 0.10 | 0.99 | 0.09 | 2.43
High | -1 | 0.50 | 1 | 0.16 | 2.86
High | -1 | 0.90 | 0.99 | 0.32 | 9.25
High | 0.00 | 0.10 | 0.99 | 0.14 | 2.62
High | 0.00 | 0.50 | 0.99 | 0.19 | 3.11
High | 0.00 | 0.90 | 1.02 | 0.58 | 20.46
High | 1 | 0.10 | 1 | 0.11 | 2.51
High | 1 | 0.50 | 1.01 | 0.19 | 3.74
High | 1 | 0.90 | 1.01 | 0.43 | 12.29
Table 15: Summary statistics of inverse probability of IPTWs distribution in various simulation scenarios when $n=5000$ Outcome_prevalence | Treatment_prevalence | Confounding | Mean | standard deviation | Max
---|---|---|---|---|---
Low | -1 | 0.10 | 1 | 0.10 | 2.52
Low | -1 | 0.50 | 1 | 0.19 | 5.06
Low | -1 | 0.90 | 1 | 0.38 | 22.74
Low | 0.00 | 0.10 | 0.99 | 0.13 | 4.05
Low | 0.00 | 0.50 | 0.99 | 0.24 | 4.69
Low | 0.00 | 0.90 | 0.99 | 0.45 | 15.12
Low | 1 | 0.10 | 1 | 0.10 | 2.49
Low | 1 | 0.50 | 1 | 0.20 | 3.54
Low | 1 | 0.90 | 1 | 0.39 | 13.06
Medium | -1 | 0.10 | 1 | 0.09 | 2.13
Medium | -1 | 0.50 | 1 | 0.17 | 4.61
Medium | -1 | 0.90 | 0.99 | 0.29 | 9.82
Medium | 0.00 | 0.10 | 0.99 | 0.13 | 3.17
Medium | 0.00 | 0.50 | 0.99 | 0.24 | 5.97
Medium | 0.00 | 0.90 | 1 | 0.43 | 11.59
Medium | 1 | 0.10 | 1 | 0.10 | 3.02
Medium | 1 | 0.50 | 1 | 0.20 | 3.55
Medium | 1 | 0.90 | 0.99 | 0.37 | 11.27
High | -1 | 0.10 | 1 | 0.07 | 2.51
High | -1 | 0.50 | 1 | 0.19 | 6.27
High | -1 | 0.90 | 0.99 | 0.36 | 29.88
High | 0.00 | 0.10 | 0.99 | 0.12 | 2.43
High | 0.00 | 0.50 | 0.99 | 0.23 | 7.53
High | 0.00 | 0.90 | 0.99 | 0.39 | 8.35
High | 1 | 0.10 | 1 | 0.08 | 2.46
High | 1 | 0.50 | 1 | 0.22 | 4.78
High | 1 | 0.90 | 0.99 | 0.34 | 6.38

Figure 13: Coverage of the CIs under medium event rates. Bootstrap: CIs
constructed by nonparametric bootstrap; LEF both: CIs constructed by applying
Approach 2 of LEF bootstrap; LEF outcome: CIs constructed by applying Approach
1 of LEF bootstrap; Jackknife Wald: CIs constructed by applying Approach 1 of
Jackknife resampling; Jackknife MVN: CIs constructed by applying Approach 2 of
Jackknife resampling; Sandwich: CIs based on the sandwich variance estimator

Figure 14: Coverage of the CIs under high event rates. Bootstrap: CIs
constructed by nonparametric bootstrap; LEF both: CIs constructed by applying
Approach 2 of LEF bootstrap; LEF outcome: CIs constructed by applying Approach
1 of LEF bootstrap; Jackknife Wald: CIs constructed by applying Approach 1 of
Jackknife resampling; Jackknife MVN: CIs constructed by applying Approach 2 of
Jackknife resampling; Sandwich: CIs based on the sandwich variance estimator

Figure 15: Ratio of mean estimated standard error to empirical standard
deviation of MRD sampling distribution under medium event rates. Bootstrap:
CIs constructed by nonparametric bootstrap; LEF both: CIs constructed by
applying Approach 2 of LEF bootstrap; LEF outcome: CIs constructed by applying
Approach 1 of LEF bootstrap; Jackknife Wald: CIs constructed by applying
Approach 1 of Jackknife resampling; Jackknife MVN: CIs constructed by applying
Approach 2 of Jackknife resampling; Sandwich: CIs based on the sandwich
variance estimator

Figure 16: Ratio of mean estimated standard error to empirical standard
deviation of MRD sampling distribution under high event rates. Bootstrap: CIs
constructed by nonparametric bootstrap; LEF both: CIs constructed by applying
Approach 2 of LEF bootstrap; LEF outcome: CIs constructed by applying Approach
1 of LEF bootstrap; Jackknife Wald: CIs constructed by applying Approach 1 of
Jackknife resampling; Jackknife MVN: CIs constructed by applying Approach 2 of
Jackknife resampling; Sandwich: CIs based on the sandwich variance estimator

Figure 17: Ratio of mean estimated standard error to empirical standard
deviation of MRD sampling distribution under high event rates. Bootstrap: CIs
constructed by nonparametric bootstrap; LEF both: CIs constructed by applying
Approach 2 of LEF bootstrap; LEF outcome: CIs constructed by applying Approach
1 of LEF bootstrap; Jackknife Wald: CIs constructed by applying Approach 1 of
Jackknife resampling; Jackknife MVN: CIs constructed by applying Approach 2 of
Jackknife resampling; Sandwich: CIs based on the sandwich variance estimator

Figure 18: Ratio of mean estimated standard error to empirical standard
deviation of MRD sampling distribution under high event rates. Bootstrap: CIs
constructed by nonparametric bootstrap; LEF both: CIs constructed by applying
Approach 2 of LEF bootstrap; LEF outcome: CIs constructed by applying Approach
1 of LEF bootstrap; Jackknife Wald: CIs constructed by applying Approach 1 of
Jackknife resampling; Jackknife MVN: CIs constructed by applying Approach 2 of
Jackknife resampling; Sandwich: CIs based on the sandwich variance estimator

Figure 19: Bias-eliminated coverage of the CIs under low event rates.
Bootstrap: CIs constructed by nonparametric bootstrap; LEF both: CIs
constructed by applying Approach 2 of LEF bootstrap; LEF outcome: CIs
constructed by applying Approach 1 of LEF bootstrap; Jackknife Wald: CIs
constructed by applying Approach 1 of Jackknife resampling; Jackknife MVN: CIs
constructed by applying Approach 2 of Jackknife resampling; Sandwich: CIs
based on the sandwich variance estimator

Figure 20: Bias-eliminated coverage of the CIs under medium event rates.
Bootstrap: CIs constructed by nonparametric bootstrap; LEF both: CIs
constructed by applying Approach 2 of LEF bootstrap; LEF outcome: CIs
constructed by applying Approach 1 of LEF bootstrap; Jackknife Wald: CIs
constructed by applying Approach 1 of Jackknife resampling; Jackknife MVN: CIs
constructed by applying Approach 2 of Jackknife resampling; Sandwich: CIs
based on the sandwich variance estimator

Figure 21: Bias-eliminated coverage of the CIs under high event rates.
Bootstrap: CIs constructed by nonparametric bootstrap; LEF both: CIs
constructed by applying Approach 2 of LEF bootstrap; LEF outcome: CIs
constructed by applying Approach 1 of LEF bootstrap; Jackknife Wald: CIs
constructed by applying Approach 1 of Jackknife resampling; Jackknife MVN: CIs
constructed by applying Approach 2 of Jackknife resampling; Sandwich: CIs
based on the sandwich variance estimator