# Generalized Benders Decomposition with Continual Learning for Hybrid Model
Predictive Control in Dynamic Environment

Xuan Lin† †Xuan Lin conducts this research as an independent researcher
<EMAIL_ADDRESS>
###### Abstract

Hybrid model predictive control (MPC) with both continuous and discrete
variables is widely applicable to robotic control tasks, especially those
involving contact with the environment. Due to the combinatorial complexity,
the solving speed of hybrid MPC can be insufficient for real-time
applications. In this paper, we proposed a hybrid MPC solver based on
Generalized Benders Decomposition (GBD) with continual learning. The algorithm
accumulates cutting planes from the invariant dual space of the subproblems.
After a short cold-start phase, the accumulated cuts provide warm-starts for
the new problem instances to increase the solving speed. Despite the randomly
changing environment that the control is unprepared for, the solving speed
maintains. We verified our solver on controlling a cart-pole system with
randomly moving soft contact walls and show that the solving speed is 2-3
times faster than the off-the-shelf solver Gurobi.

##  I Introduction

Hybrid model predictive control (MPC) with both continuous and discrete
variables is widely applicable to robotic control tasks, especially those
involving contact with the environment. However, discontinuous variables are
oftentimes computed offline for Hybrid MPC [1, 2, 3, 4] due to their
combinatorial complexities. These include gaits for legged robots and contact
sequences for manipulation tasks. Several models with mixed discrete-
continuous variables were proposed including mixed-logic dynamic systems
(MLDs) [5], linear complementary models (LCs) [6], and piece-wise affine
systems (PWAs) [7]. Their conditional equivalences were established in [8]
(for example, LCs are equivalent to MLDs provided that the complementary
variables are bounded). Despite several recent works that solve MPC on these
systems [9, 10, 11], the problems addressed in those papers are demonstrate in
static environments. In real robotic applications, it is beneficial to further
increase the solving speed to reduce the control error since the models are
never accurate.

In this paper, we propose a novel hybrid MPC solver based on Generalized
Benders decomposition (GBD) [12] to solve problems including MLD constraints
under changing environments. Benders decomposition separates the problem into
a master problem which solves part of the variables named complicating
variables, and a subproblem which solves the rest of the variables. It uses a
delayed constraint generation technique that builds up representations of the
feasible region and optimal cost function for complicating variables inside
the master problem. These representations are constructed as cutting planes
based on the dual solutions of the subproblem. The key idea we propose is to
accumulate the cuts as more problem instances are solved since the dual
feasible set is invariant under the changing environments. The accumulated
cuts then provide warm-starts for new problem instances. As a result, the
solving speed keeps on increasing. Under the fastest situation, GBD only needs
to solve one master problem and one subproblem to get a globally optimal
solution. The proposed solver is compared against the recent works on warm-
started Branch-Bound solvers and the commercialized solver Gurobi.

We list the contributions below:

  1. 1.

We propose a novel solver based on Generalized Benders Decomposition for
hybrid MPCs, where the accumulated cuts provide warm-starts for the new
problem instances leading to faster solving speeds despite changing
environments.

  2. 2.

We tested our solver on controlling a cart-pole system with randomly moving
soft contact walls, a more challenging test than cart-pole with static walls
prevail in previous literature [11, 13, 9, 14]. We show that our GBD solver
runs faster on average than warm-started Branch and Bound solvers and the off-
the-shelf solver Gurobi.

_Notations_ Vectors are bold lowercase; matrices are bold uppercase; sets are
script or italicized uppercase. The real number set is $\mathbb{R}$. For
$\boldsymbol{x},\boldsymbol{y}\in\mathbb{R}^{n}$,
$\boldsymbol{x}\leq\boldsymbol{y}$ indicates element-wise inequality. For
$\boldsymbol{A}\in\mathbb{R}^{n\times n}$ and
$\boldsymbol{B}\in\mathbb{R}^{m\times m}$,
$\text{diag}(\boldsymbol{A},\boldsymbol{B})\in\mathbb{R}^{(n+m)\times(n+m)}$
denotes the block diagonal matrix with diagonal blocks $\boldsymbol{A}$ and
$\boldsymbol{B}$, and zeros otherwise. $\boldsymbol{I}_{n}$ denotes an
identity matrix of dimension $n$. The open ball
$\textit{B}_{\epsilon}(\boldsymbol{p})$ denotes
$\\{\boldsymbol{q}:||\boldsymbol{p}-\boldsymbol{q}||<\epsilon\\}$.

##  II Related Works

###  II-A Mixed-Logic Dynamic Models (MLDs)

In [5], the authors proposed mixed-logic dynamics systems as a general
modeling tool for control systems incorporating physical laws, logic rules and
operating constraints. MLD is rich enough to incorporate dynamic systems such
as finite state machines, nonlinear systems that can be written as piece-wise
linear form. MLDs have been widely used to model energy storage systems [15],
transportation systems [16], temperature management systems [17], to name a
few. Recently, MLDs and its equivalent models such as linear complementary
models are introduced into the robotics locomotion and manipulation community
to model real-time control involving contacts [11, 10, 9, 18].

MLDs incorporate states and inputs that can be a mixture of continuous and
discrete variables, and quadratic objective functions. Since MLDs incorporate
a mixture of discrete and continuous variables, solving it online for fast
motion planning and model-predictive control demands an efficient MIQP solver.
Several methods have been proposed including explicit MPC [19], Branch-and-
Bound [11], ADMM [9], Lagrange relaxation [20], elastic modes [21], and
generalized Benders decomposition [22]. We briefly go over them below.

Explicit MPC solves the problem completely or partially offline, such that the
solver only picks out solutions online. [23] building polyhedron regions of
invariant optimal control function for LQR problems. Explicitly building
polytope regions is computationally expensive hence limited to simple
problems. [24] stored all solutions and used a K-nearest-neighbor approach to
pick out binary solutions. [19] explored a partial explicit approach that
combines offline solved binary solutions with online convex computation.
Related to this work is the library building approach where the problem is
solved offline and recorded into a dataset. The solver then picks out data
points and uses them as warm-starts. The online data selection approach can be
a K-nearest neighbor classifier [25, 26], or a learned neural-network [18].
However, this approach generally has difficulty facing out-of-distribution
scenarios from the dataset. Except for [11], the works mentioned above only
use offline optimal solutions. The infeasible or suboptimal solutions are not
used.

Branch and bound is a common approach to solve mixed-integer programs used by
off-the-shelf solvers such as Gurobi. This approach first relaxes the integer
programming problem into a convex programming on the root node where all
binaries are extended to continuous variables between zero and one. It then
fixes the binary variables one-by-one, and generate branches from the root to
a solution on the leaf node where all binary variables are fixed. Previous
studies tried to use information from previous solve to warm-start the new
solve online. [27] studied propagating the path from root to leaf from the
previous iteration to the next one for warm-start, such that a number of
parent nodes do not need to be resolved. [11] further explored propagating
complete B&B tree to warm-start the next problem. Even with the proposed
techniques, B&B can still be slow as it has too many subproblems to keep track
of, particularly under noise and model inaccuracies.

Another approach to solve mixed-integer programs is through alternating
direction method of multipliers (ADMM). ADMM solves two or multiple problems
iteratively until they reach a consensus through a penalty in the objective
function. Computer vision and operation research community has used ADMM to
solve large scale MIP problems [28]. In the robotics community, ADMM has been
implemented for solving complementary control problems [9] at a fast speed.
Despite [9] does not discuss it, ADMM allows for easy warm-start [29], using
previous solution to accelerate the solving of the next solution. On the other
hand, ADMM does not have convergence guarantee for MIP problems unless special
assumptions are made such as [28].

###  II-B Benders Decomposition

Benders decomposition [30] can be regarded as Danzig-Wolfe decomposition [31]
applied to the dual. Both of them use delayed column or constraint generation
techniques. Benders decomposition identifies the complicating variables and
defines a subproblem such that those variables are fixed. For this technique
to work well, the subproblem should be much easier to solve than the complete
problem. For mixed-integer programming, the subproblem is the convex part with
complicating variables being the discrete variables [22]. As subproblems are
solved, cutting planes are added to the master problem to build the feasible
set and optimal cost function for the subproblem.

Benders decomposition was originally proposed to solve linear duals e.g.
MILPs. In [12], the author proposed Generalized Benders decomposition (GBD)
that extends the theory to nonlinear duals. Several authors have investigated
solving MIQPs using GBD [22, 32, 33, 34]. In [35, 36], the authors propose
logic-based Benders decomposition which further generalized the theory to so-
called inference dual, which is a logic combination of propositions. This
method extends the application of BD to planning and scheduling problems such
as satisfiability of 0-1 programming problems whose dual is not a traditional
linear or nonlinear programming problem. Using this idea, [37] proposed a
formulation of combinatorial Benders feasibility cut for MILPs that does not
depend on the big-M constant.

As Benders decomposition involves master-subproblem structure, it suits the
large-scale distributed problems, or problems with a large number of possible
scenarios like stochastic programs [38]. For applications such as distributed
control [39], the subproblems can be decoupled into multiple smaller-scale
problems and solved in parallel to reduce the computation demand. As pointed
out by the review paper [40], many authors report over 90% solving time spent
on the master problem. Therefore, a number of previous work investigated on
how to speed up the master problem, or use its results more efficiently.
Examples include local branching heuristics [41], heuristic master problem
solutions [42], generating pareto-optimal cuts [43], cut initialization [44],
valid inequalities [45], etc. See [40] for a comprehensive review of these
methods. [46] points out that classic Benders feasibility cuts do not carry
objective function value leading to convergence issues. They proposed
additional feasibility cuts to resolve this issue.

GBD can also be used to learn objective functions. This has been applied to
dual dynamic programming for MPC over long-term problems [47, 48]. Previous
work [49, 50] uses Benders cuts to construct lower bounds for infinitely long
objective functions using Bellman operators for both nonlinear and mixed-
integer linear systems. Through learning Benders cuts from offline dataset,
one avoids hand-tuning terminal cost of objective function. Despite a more
optimal objective being learned, the online solving speed of MIP is invariant
of objective functions.

##  III Problem Model

We develop MPC control laws for Mixed Logic Dynamic (MLD) systems as proposed
by [5]:

$\displaystyle\boldsymbol{x}_{k+1}$
$\displaystyle=\boldsymbol{E}\boldsymbol{x}_{k}+\boldsymbol{F}\boldsymbol{u}_{k}+\boldsymbol{G}\boldsymbol{\delta}_{k}+\boldsymbol{n}_{k}$
(1a) $\displaystyle\boldsymbol{H}_{1}\boldsymbol{x}_{k}$
$\displaystyle+\boldsymbol{H}_{2}\boldsymbol{u}_{k}+\boldsymbol{H}_{3}\boldsymbol{\delta}_{k}\leq\boldsymbol{h}(\boldsymbol{\theta})$
(1b)

At time $k$, $\boldsymbol{x}_{k}\in\mathbb{R}^{n_{x}}$ is the continuous
state. $\boldsymbol{u}_{k}\in\mathbb{R}^{n_{u}}$ denotes the continuous input.
$\boldsymbol{\delta}_{k}\in\\{0,1\\}^{n_{\delta}}$ is the binary input.
$\boldsymbol{n}_{k}\in\mathbb{R}^{n_{x}}$ is the disturbance input. Matrices
representing system dynamics are $\boldsymbol{E}\in\mathbb{R}^{n_{x}\times
n_{x}}$, $\boldsymbol{F}\in\mathbb{R}^{n_{x}\times n_{u}}$,
$\boldsymbol{G}\in\mathbb{R}^{n_{x}\times n_{\delta}}$.
$\boldsymbol{H}_{1}\in\mathbb{R}^{n_{c}\times n_{x}}$.
$\boldsymbol{H}_{2}\in\mathbb{R}^{n_{c}\times n_{u}}$.
$\boldsymbol{H}_{3}\in\mathbb{R}^{n_{c}\times n_{\delta}}$. The submatrices
are of appropriate dimensions. The right-hand side of the constraint (1b) is
$\boldsymbol{h}\in\mathbb{R}^{n_{c}}$ where $n_{c}$ is the number of
inequality constraints. $\boldsymbol{\theta}$ parameterizes $\boldsymbol{h}$
to represent the changing environments. We assume that matrices
$\boldsymbol{E}$, $\boldsymbol{F}$, $\boldsymbol{G}$, $\boldsymbol{H}_{1}$,
$\boldsymbol{H}_{2}$, $\boldsymbol{H}_{3}$ are independent of
$\boldsymbol{\delta}_{k}$ and $\boldsymbol{\theta}$. This makes
$\boldsymbol{\delta}_{k}$ and $\boldsymbol{\theta}$ as inputs to the system
while the inherent physics of the system is invariant.

######  Remark.

The goal of parameter $\boldsymbol{\theta}$ is to represent a sudden change in
the environment that the controller is uninformed of and cannot prepare for it
down the MPC horizon. Note this is different from the time-varying system
investigated by [11] where the controller is well-informed of the change in
advance ($\boldsymbol{\theta}_{t},t=0,...,T$ is known).

We formulate a hybrid MPC for this system. The MPC formulation solves an
optimization problem to get a sequence of control inputs. However, only the
first one is used. It then takes the sensor feedback and resolves the problem.
If this could be done fast enough on the hardware, the robot can reject
disturbances. The MPC formulation is:

$\displaystyle\underset{\boldsymbol{x}_{k}\in X_{k},\ \boldsymbol{u}_{k},\
\boldsymbol{\delta}_{k}}{\text{minimize}}$ $\displaystyle\ \
\sum_{k=0}^{N-1}\boldsymbol{x}_{k}^{T}\boldsymbol{Q}_{k}\boldsymbol{x}_{k}+\boldsymbol{u}_{k}^{T}\boldsymbol{R}_{k}\boldsymbol{u}_{k}+\boldsymbol{x}_{N}^{T}\boldsymbol{Q}_{N}\boldsymbol{x}_{N}$
s.t. $\displaystyle\boldsymbol{x}_{0}=\boldsymbol{x}_{ini}$
$\displaystyle\boldsymbol{x}_{k+1}=\boldsymbol{E}\boldsymbol{x}_{k}+\boldsymbol{F}\boldsymbol{u}_{k}+\boldsymbol{G}\boldsymbol{\delta}_{k}$
$\displaystyle\boldsymbol{H}_{1}\boldsymbol{x}_{k}+\boldsymbol{H}_{2}\boldsymbol{u}_{k}+\boldsymbol{H}_{3}\boldsymbol{\delta}_{k}\leq\boldsymbol{h}(\boldsymbol{\theta})$
$\displaystyle\boldsymbol{\delta}_{k}\in\\{0,1\\}^{n_{\delta}},\ k=0,...,N-1$

The matrices $\boldsymbol{Q}_{k}$ and $\boldsymbol{Q}_{N}$ are positive
definite matrices. $X_{k}$ is the domain of $\boldsymbol{x}_{k}$. The system
is written into a more compact form:

$\displaystyle\underset{\boldsymbol{x}\in X,\
\boldsymbol{\delta}}{\text{minimize}}$
$\displaystyle\boldsymbol{x}^{T}\boldsymbol{Q}\boldsymbol{x}$ (2) s.t.
$\displaystyle\boldsymbol{A}\boldsymbol{x}=\boldsymbol{b}(\boldsymbol{x}_{ini},\boldsymbol{\delta})$
$\displaystyle\boldsymbol{C}\boldsymbol{x}\leq\boldsymbol{d}(\boldsymbol{\theta},\boldsymbol{\delta})$
$\displaystyle\boldsymbol{\delta}_{k}\in\\{0,1\\}^{n_{\delta}}$

Let $n_{xu}=n_{x}+n_{u}$. The matrices and vectors have structures:

$\displaystyle\boldsymbol{x}$
$\displaystyle=\begin{bmatrix}\boldsymbol{x}^{T}_{0}&\boldsymbol{u}^{T}_{0}&\cdots&\boldsymbol{x}^{T}_{N-1}&\boldsymbol{u}^{T}_{N-1}&\boldsymbol{x}^{T}_{N}\end{bmatrix}^{T}\in\mathbb{R}^{Nn_{xu}+n_{x}}$
(3a) $\displaystyle\boldsymbol{\delta}$
$\displaystyle=\begin{bmatrix}\boldsymbol{\delta}^{T}_{0}&\cdots&\boldsymbol{\delta}^{T}_{N}\end{bmatrix}^{T}\in\mathbb{R}^{(N+1)n_{\delta}}$
(3b)

Hence domain of $\boldsymbol{x}$ is
$X=X_{k}\times\mathbb{R}^{n_{u}}\times\cdots\times\mathbb{R}^{n_{u}}\times
X_{k}$.

$\displaystyle\boldsymbol{A}=$
$\displaystyle\begin{bmatrix}\boldsymbol{I}_{n_{x}}\ \ \boldsymbol{0}\\\
-\boldsymbol{E}\ -\boldsymbol{F}&\boldsymbol{I}_{n_{x}}\ \ \boldsymbol{0}\\\
&-\boldsymbol{E}\ -\boldsymbol{F}&\ddots\\\ &&\ddots&\boldsymbol{I}_{n_{x}}\ \
\boldsymbol{0}\\\ &&&-\boldsymbol{E}\
-\boldsymbol{F}&\boldsymbol{I}_{n_{x}}\end{bmatrix}$ (4)
$\displaystyle\in\mathbb{R}^{(N+1)n_{x}\times(Nn_{xu}+n_{x})}$

$\displaystyle\boldsymbol{b}(\boldsymbol{x}_{ini},\boldsymbol{\delta})$
$\displaystyle=\begin{bmatrix}\boldsymbol{x}_{ini}^{T}&(\boldsymbol{G}\boldsymbol{\delta}_{0})^{T}&\cdots&(\boldsymbol{G}\boldsymbol{\delta}_{N-1})^{T}\end{bmatrix}^{T}$
(5) $\displaystyle\in\mathbb{R}^{(N+1)n_{x}}$

$\displaystyle\boldsymbol{C}=$
$\displaystyle\begin{bmatrix}\boldsymbol{H}_{1}\ \boldsymbol{H}_{2}\\\
&\boldsymbol{H}_{1}\ \boldsymbol{H}_{2}\\\ &&\ddots\\\ &&&\boldsymbol{H}_{1}\
\boldsymbol{H}_{2}&\boldsymbol{0}\end{bmatrix}$ (6)
$\displaystyle\in\mathbb{R}^{Nn_{c}\times(Nn_{xu}+n_{x})}$

$\displaystyle\boldsymbol{d}(\boldsymbol{\theta},\boldsymbol{\delta})$
$\displaystyle=\begin{bmatrix}(\boldsymbol{h}(\boldsymbol{\theta})-\boldsymbol{H}_{3}\boldsymbol{\delta}_{0})^{T}&\cdots&(\boldsymbol{h}(\boldsymbol{\theta})-\boldsymbol{H}_{3}\boldsymbol{\delta}_{N-1})^{T}\end{bmatrix}^{T}$
(7) $\displaystyle\in\mathbb{R}^{Nn_{c}}$

$\boldsymbol{Q}=\text{diag}(\boldsymbol{Q}_{k},\boldsymbol{R}_{k})\in\mathbb{R}^{(Nn_{xu}+n_{x})\times(Nn_{xu}+n_{x})}$
(8)

Problem (2) is an MIQP and can be solved through an off-the-shelf mixed-
integer convex programming solver based on Branch and Bound, such as Gurobi.
However, Branch and Bound algorithms keep track of a large number of
subproblems that relax the binary constraints in different ways. Despite the
MPC warm-start scheme such as shifting contact sequence can be used [11], many
subproblems still need to be solved for a new problem instance. For
applications that require extremely fast solving speed, this can be
insufficient. In this paper, we propose to use generalized Benders
decomposition to solve the problem several times faster than Gurobi.

##  IV Benders decomposition formulation

In this section, we apply Benders decomposition to our hybrid MPC problem
(III). Benders decomposition deals with the problem of the following form:

$\displaystyle\underset{\boldsymbol{x},\boldsymbol{y}}{\text{minimize}}\ \
f(\boldsymbol{x},\boldsymbol{y})$ (9) $\displaystyle\text{s.t.}\ \
\boldsymbol{G}(\boldsymbol{x},\boldsymbol{y})\leq 0$
$\displaystyle\boldsymbol{x}\in X,\boldsymbol{y}\in Y$

where $\boldsymbol{y}$ is a vector of complicating variables. If
$\boldsymbol{y}$ is fixed, the optimization problem is much easier to solve.
Benders decomposition partitions the problem into a master problem by
projecting onto the $\boldsymbol{y}$ space:

$\displaystyle\underset{\boldsymbol{y}}{\text{minimize}}\ \ v(\boldsymbol{y})$
(10) $\displaystyle\text{s.t.}\ \ \boldsymbol{y}\in Y\cap V$

The function $v(\boldsymbol{y})$ is defined to provide the best objective
function with fixed complicating variable $\boldsymbol{y}$:

$\displaystyle v(\boldsymbol{y})=\underset{\boldsymbol{x}}{\text{infimum}}\ \
f(\boldsymbol{x},\boldsymbol{y})$ (11) $\displaystyle\text{s.t.}\ \
\boldsymbol{G}(\boldsymbol{x},\boldsymbol{y})\leq 0$
$\displaystyle\boldsymbol{x}\in X$

$V$ contains all $\boldsymbol{y}$’s such that problem (11) is feasible:

$V=\\{\boldsymbol{y}:\boldsymbol{G}(\boldsymbol{x},\boldsymbol{y})\leq 0,\ \
\exists\boldsymbol{x}\in X\\}$ (12)

For our hybrid MPC, we define the complicating variable $\boldsymbol{y}$ as
the binary variable $\boldsymbol{\delta}$, the initial condition
$\boldsymbol{x}_{ini}$, and the parameter $\boldsymbol{\theta}$. The
subproblem is:

$\displaystyle
v(\boldsymbol{x}_{ini},\boldsymbol{\theta},\boldsymbol{\delta})=\underset{\boldsymbol{x}\in
X}{\text{minimize}}$
$\displaystyle\boldsymbol{x}^{T}\boldsymbol{Q}\boldsymbol{x}$ (13) s.t.
$\displaystyle\boldsymbol{A}\boldsymbol{x}=\boldsymbol{b}(\boldsymbol{x}_{ini},\boldsymbol{\delta})$
$\displaystyle\boldsymbol{C}\boldsymbol{x}\leq\boldsymbol{d}(\boldsymbol{\theta},\boldsymbol{\delta})$

Given fixed $\boldsymbol{x}_{ini}$, $\boldsymbol{\theta}$,
$\boldsymbol{\delta}$, (13) is a quadratic programming and can be solved
through off-the-shelf QP solvers. The master problem is:

$\displaystyle\underset{\boldsymbol{\delta}}{\text{minimize}}$ $\displaystyle
v(\boldsymbol{x}_{ini},\boldsymbol{\theta},\boldsymbol{\delta})$ (14) s.t.
$\displaystyle\boldsymbol{\delta}_{k}\in\\{0,1\\}$
$\displaystyle\boldsymbol{\delta}\in
V\coloneqq\\{\boldsymbol{\delta}:\boldsymbol{A}\boldsymbol{x}=\boldsymbol{b}(\boldsymbol{x}_{ini},\boldsymbol{\delta})$
$\displaystyle\boldsymbol{C}\boldsymbol{x}\leq\boldsymbol{d}(\boldsymbol{\theta},\boldsymbol{\delta}),\exists\boldsymbol{x}\in
X\\}$

The essential issue with solving (14) is that function
$v(\boldsymbol{x}_{ini},\boldsymbol{\theta},\boldsymbol{\delta})$ and set $V$
are only implicitly known through their definitions. Benders decomposition is
a process that iteratively solves problem (10) and (11) to build
approximations of $v$ and $V$ in the problem (10).

We will constantly work with the dual of problem (13), given the advantage
that the dual is invariant with respect to the complicating variables. We
derive the dual for reference. Recall the definition of Lagrangian for the
subproblem (13):

$\displaystyle\mathcal{L}(\boldsymbol{x},\boldsymbol{\nu},\boldsymbol{\lambda};\boldsymbol{x}_{ini},\boldsymbol{\theta},\boldsymbol{\delta})=f_{obj}(\boldsymbol{x})$
$\displaystyle+\boldsymbol{\nu}^{T}(\boldsymbol{A}\boldsymbol{x}-\boldsymbol{b}(\boldsymbol{x}_{ini},\boldsymbol{\delta}))$
(15)
$\displaystyle+\boldsymbol{\lambda}^{T}(\boldsymbol{C}\boldsymbol{x}-\boldsymbol{d}(\boldsymbol{\theta},\boldsymbol{\delta}))$

where $\boldsymbol{\nu}\in\mathbb{R}^{(N+1)n_{x}}$,
$\boldsymbol{\lambda}\in\mathbb{R}^{Nn_{c}}$ are the dual variables associated
with
$\boldsymbol{A}\boldsymbol{x}=\boldsymbol{b}(\boldsymbol{x}_{ini},\boldsymbol{\delta})$,
$\boldsymbol{C}\boldsymbol{x}\leq\boldsymbol{d}(\boldsymbol{\theta},\boldsymbol{\delta})$,
respectively. The Lagrange dual function $g$ is:

$g(\boldsymbol{\nu},\boldsymbol{\lambda};\boldsymbol{x}_{ini},\boldsymbol{\theta},\boldsymbol{\delta})=\underset{\boldsymbol{x}\in
X}{\text{minimize}}\ \
\mathcal{L}(\boldsymbol{x},\boldsymbol{\nu},\boldsymbol{\lambda};\boldsymbol{x}_{ini},\boldsymbol{\theta},\boldsymbol{\delta})$
(16)

where $f_{obj}=\boldsymbol{x}^{T}\boldsymbol{Q}\boldsymbol{x}$. Let
$\boldsymbol{x}^{0}$ be the unconstrained minimizer of $\mathcal{L}$ (in
general, $\boldsymbol{x}^{0}$ different from the opitmal primal solution
$\boldsymbol{x}^{*}$). By taking derivative, we have:

$\boldsymbol{x}^{0}=-\frac{1}{2}\boldsymbol{Q}^{-1}(\boldsymbol{A}^{T}\boldsymbol{\nu}+\boldsymbol{C}^{T}\boldsymbol{\lambda})$
(17)

Hence the Lagrange dual problem is:

$\displaystyle\underset{\boldsymbol{\nu},\
\boldsymbol{\lambda}}{\text{maximize}}$ $\displaystyle\ \
-\frac{1}{4}||\boldsymbol{A}^{T}\boldsymbol{\nu}+\boldsymbol{C}^{T}\boldsymbol{\lambda}||^{2}_{\boldsymbol{Q}^{-1}}$
(18) $\displaystyle\ \
-\boldsymbol{b}(\boldsymbol{x}_{ini},\boldsymbol{\delta})^{T}\boldsymbol{\nu}-\boldsymbol{d}(\boldsymbol{\theta},\boldsymbol{\delta})^{T}\boldsymbol{\lambda}$
s.t. $\displaystyle\ \ \boldsymbol{\lambda}\geq\boldsymbol{0}$

As the feasibility of (13) is independent of the objective function, (13) is
feasible if and only if the following problem is feasible:

$\displaystyle\underset{\boldsymbol{x}\in X}{\text{minimize}}$
$\displaystyle\boldsymbol{0}$ (19) s.t.
$\displaystyle\boldsymbol{A}\boldsymbol{x}=\boldsymbol{b}(\boldsymbol{x}_{ini},\boldsymbol{\delta})$
$\displaystyle\boldsymbol{C}\boldsymbol{x}\leq\boldsymbol{d}(\boldsymbol{\theta},\boldsymbol{\delta})$

Problem (19) has the dual:

$\displaystyle\underset{\boldsymbol{\nu},\
\boldsymbol{\lambda}}{\text{maximize}}$ $\displaystyle\ \
-\boldsymbol{b}(\boldsymbol{x}_{ini},\boldsymbol{\delta})^{T}\boldsymbol{\nu}-\boldsymbol{d}(\boldsymbol{\theta},\boldsymbol{\delta})^{T}\boldsymbol{\lambda}$
(20) s.t.
$\displaystyle\boldsymbol{A}^{T}\boldsymbol{\nu}+\boldsymbol{C}^{T}\boldsymbol{\lambda}=\boldsymbol{0}$
$\displaystyle\boldsymbol{\lambda}\geq\boldsymbol{0}$

###  IV-A Feasibility cuts

If at iteration $p$, the subproblem is infeasible under the given
$\boldsymbol{\delta}_{p}$, this $\boldsymbol{\delta}_{p}$ needs to be removed
from the master problem. This can be achieved by adding a cutting plane. Since
the problem (III) is linearly constrained, the Farkas certificates can be used
to add feasibility cuts. They can be discovered by solving (19) with a dual
simplex solver ([51], Chapter 6.5). The theorem of alternatives for (19) is:

######  Lemma 1.

Given $\boldsymbol{A}\in\mathbb{R}^{l\times n}$,
$\boldsymbol{b}\in\mathbb{R}^{l}$, $\boldsymbol{C}\in\mathbb{R}^{m\times n}$,
$\boldsymbol{d}\in\mathbb{R}^{m}$, exactly one of the following statements is
true:

  1. 1.

There exists an $\boldsymbol{x}\in\mathbb{R}^{n}$ that satisfies
$\boldsymbol{A}\boldsymbol{x}=\boldsymbol{b},\boldsymbol{C}\boldsymbol{x}\leq\boldsymbol{d}$.

  2. 2.

There exist $\boldsymbol{y}\in\mathbb{R}^{l}$,
$\boldsymbol{z}\in\mathbb{R}^{m}$ that satisfy
$\boldsymbol{z}\geq\boldsymbol{0}$,
$\boldsymbol{A}^{T}\boldsymbol{y}+\boldsymbol{C}^{T}\boldsymbol{z}=\boldsymbol{0}$,
$\boldsymbol{b}^{T}\boldsymbol{y}+\boldsymbol{d}^{T}\boldsymbol{z}<0$.

###### Proof.

See Appendix A. ∎

If (19) is infeasible for $\boldsymbol{\delta}_{p}$. Then we can add a cutting
plane to the master problem to remove a set of $\boldsymbol{\delta}$’s
including $\boldsymbol{\delta}_{p}$. Farkas lemma guarantees the existence of
$\tilde{\boldsymbol{\nu}}_{p}\in\mathbb{R}^{(N+1)n_{x}}$,
$\tilde{\boldsymbol{\lambda}}_{p}\in\mathbb{R}^{Nn_{c}}$ such that:

$\displaystyle\tilde{\boldsymbol{\lambda}}_{p}\geq\boldsymbol{0}$ (21)
$\displaystyle\boldsymbol{A}^{T}\tilde{\boldsymbol{\nu}}_{p}+\boldsymbol{C}^{T}\tilde{\boldsymbol{\lambda}}_{p}=\boldsymbol{0}$
$\displaystyle\boldsymbol{b}(\boldsymbol{x}_{ini},\boldsymbol{\delta}_{p})^{T}\tilde{\boldsymbol{\nu}}_{p}+\boldsymbol{d}(\boldsymbol{\theta},\boldsymbol{\delta}_{p})^{T}\tilde{\boldsymbol{\lambda}}_{p}<0$

To prevent the master problem from giving this $\boldsymbol{\delta}_{p}$, a
constraint to defeat the Farkas infeasible proof is added to the master
problem:

$\boldsymbol{b}(\boldsymbol{x}_{ini},\boldsymbol{\delta})^{T}\tilde{\boldsymbol{\nu}}_{p}+\boldsymbol{d}(\boldsymbol{\theta},\boldsymbol{\delta})^{T}\tilde{\boldsymbol{\lambda}}_{p}\geq
0$ (22)

We note that this cutting plane will not remove any feasible
$\boldsymbol{\delta}$ from the subproblem. We state this as a lemma.

######  Lemma 2.

For given $\boldsymbol{x}_{ini}$ and $\boldsymbol{\theta}$, any
$\boldsymbol{\delta}$ that contradicts (22) proves infeasibility for (19).

###### Proof.

As $\tilde{\boldsymbol{\nu}}_{p}$, $\tilde{\boldsymbol{\lambda}}_{p}$
discovered by the dual simplex solver satisfy the first two conditions of
(21), they are feasible for the dual problem (20). Let $a\in\mathbb{R}^{+}$ be
an arbitrary positive value,
$(a\tilde{\boldsymbol{\nu}}_{p},a\tilde{\boldsymbol{\lambda}}_{p})$ are also
feasible for (20). Let $\boldsymbol{\delta}$ be any value that contradicts
(22), we have
$-a\tilde{\boldsymbol{\nu}}_{p}^{T}\boldsymbol{b}(\boldsymbol{x}_{ini},\boldsymbol{\delta})-a\tilde{\boldsymbol{\lambda}}_{p}^{T}\boldsymbol{d}(\boldsymbol{\theta},\boldsymbol{\delta})\rightarrow+\infty$
as $a\rightarrow+\infty$, hence the dual problem is unbounded which proves
that the primal problem (19) is infeasible (from Corollary 4.1 of [51]). ∎

Since hybrid MPC needs to be solved fast online, it is important to maximize
the usage of computations so the number of iteration to find a feasible
solution is reduced. Many previous works added one feasibility cut each
iteration. Some previous works [43, 52, 53] propose adding multiple cutting
planes each iteration, or re-formulate the problem such that stronger cuts can
be generated. However, the subproblem structure has not been explored by those
papers. We propose an innovative technique to add multiple feasibility cuts to
the master problem via subproblem recursive structure. The online computation
time prevents us to solve any additional optimization problems (even convex
ones), but those cutting planes can be retrieved without any additional
computation given the planes we already have.

Define $\tilde{\boldsymbol{\nu}}_{p}^{m}$,
$\tilde{\boldsymbol{\lambda}}_{p}^{m}$, $m=1,...,N-1$ such that:

$\displaystyle\tilde{\boldsymbol{\nu}}_{p,k}^{m}$
$\displaystyle=\begin{cases}\tilde{\boldsymbol{\nu}}_{p,k+m}\ \ &\forall
k+m\leq N\\\ \boldsymbol{0}\ \ &\forall k+m>N\\\ \end{cases}$ (23)
$\displaystyle\tilde{\boldsymbol{\lambda}}_{p,k}^{m}$
$\displaystyle=\begin{cases}\tilde{\boldsymbol{\lambda}}_{p,k+m}\ \ &\forall
k+m\leq N-1\\\ \boldsymbol{0}\ \ &\forall k+m>N-1\\\ \end{cases}$

For each $m$, we add an additional cutting planes:

$\boldsymbol{b}(\boldsymbol{x}_{ini},\boldsymbol{\delta})^{T}\tilde{\boldsymbol{\nu}}_{p}^{m}+\boldsymbol{d}(\boldsymbol{\theta},\boldsymbol{\delta})^{T}\tilde{\boldsymbol{\lambda}}_{p}^{m}\geq
0,\ \ m=1,...,N-1$ (24)

The addition of cuts (24) works in two ways. First, the solutions that
contradicts (24) has good optimality as predicted by the optimality cuts,
hence may be selected by the master problem as the next trial solution.
Therefore, addition of (24) eliminates those trials and accelerates the master
problem to find a feasible solution, especially in the cold start stage when
the master problem is almost empty. Second, cuts (24) predicts the future
infeasible cased by shifting the current infeasible cases into the future,
such that the future solves do not need to re-discover them. Similar to
$\tilde{\boldsymbol{\nu}}_{p}$ and $\tilde{\boldsymbol{\lambda}}_{p}$, we
present:

######  Corollary 2.1.

Any $\boldsymbol{\delta}$ that contradicts (24) proves infeasibility for (19)
with given $\boldsymbol{x}_{ini}$ and $\boldsymbol{\theta}$.

###### Proof.

We can verify that $\tilde{\boldsymbol{\nu}}_{p}^{m}$,
$\tilde{\boldsymbol{\lambda}}_{p}^{m}$ are dual feasible for any $m$ if
$\tilde{\boldsymbol{\nu}}_{p}$, $\tilde{\boldsymbol{\lambda}}_{p}$ are dual
feasible by simply plugging them into the dual feasibility constraints. This
is a simple extension of Lemma 2 . ∎

For our problem, $\boldsymbol{b}(\boldsymbol{x}_{ini},\boldsymbol{\delta})$,
$\boldsymbol{d}(\boldsymbol{\theta},\boldsymbol{\delta})$ depend on
$\boldsymbol{\delta}$ linearly, it is interesting to realize that from an
infeasible subproblem with one $\boldsymbol{\delta}$, we construct a plane
that may remove a set of infeasible $\boldsymbol{\delta}$’s. This contributes
to the efficacy of Benders decomposition as it takes usage of infeasible
samples which are usually thrown away by the methods that learn binary
solutions offline [19, 18, 26].

###  IV-B Optimality cuts

If at iteration $q$, the sub-problem is solved to optimal under given
$\boldsymbol{\delta}_{q}$, we want to add a cutting plane as a lower bound
that approaches
$v(\boldsymbol{x}_{ini},\boldsymbol{\theta},\boldsymbol{\delta})$ from below.
This can be realized through duality theory. For any $\boldsymbol{\nu}$ and
$\boldsymbol{\lambda}\geq\boldsymbol{0}$,
$g(\boldsymbol{\nu},\boldsymbol{\lambda};\boldsymbol{x}_{ini},\boldsymbol{\theta},\boldsymbol{\delta}_{q})\leq
v(\boldsymbol{x}_{ini},\boldsymbol{\theta},\boldsymbol{\delta}_{q})$. Since
the subproblem is convex and we assume there exists a stricly feasible
solution (Slater’s), strong duality is achieved and the best lower bound is
tight:

$v(\boldsymbol{x}_{ini},\boldsymbol{\theta},\boldsymbol{\delta}_{q})=\underset{\boldsymbol{\nu},\boldsymbol{\lambda}\geq\boldsymbol{0}}{\text{maximize}}\
\
g(\boldsymbol{\nu},\boldsymbol{\lambda};\boldsymbol{x}_{ini},\boldsymbol{\theta},\boldsymbol{\delta}_{q})$
(25)

Therefore, we add the best lower bound as a cutting plane to the master
problem. Let $\boldsymbol{x}_{q}^{0}$ be the unconstrained minimizer of
$\mathcal{L}$ at iteration $q$, and
$\boldsymbol{\nu}_{q}^{*},\boldsymbol{\lambda}_{q}^{*}$ be the maximizer of
$g(\boldsymbol{\nu},\boldsymbol{\lambda};\boldsymbol{x}_{ini},\boldsymbol{\theta},\boldsymbol{\delta}_{q})$,
the cutting plane takes the form:

$\displaystyle
v(\boldsymbol{x}_{ini},\boldsymbol{\theta},\boldsymbol{\delta})$
$\displaystyle\geq\mathcal{L}(\boldsymbol{x}_{q}^{0},\boldsymbol{\nu}_{q}^{*},\boldsymbol{\lambda}_{q}^{*};\boldsymbol{x}_{ini},\boldsymbol{\theta},\boldsymbol{\delta})$
(26)
$\displaystyle=f_{obj}(\boldsymbol{x}_{q}^{0})+\boldsymbol{\nu}_{q}^{*T}\boldsymbol{A}\boldsymbol{x}_{q}^{0}+\boldsymbol{\lambda}_{q}^{*T}\boldsymbol{C}\boldsymbol{x}_{q}^{0}$
$\displaystyle-\boldsymbol{\nu}_{q}^{*T}\boldsymbol{b}(\boldsymbol{x}_{ini},\boldsymbol{\delta})-\boldsymbol{\lambda}_{q}^{*T}\boldsymbol{d}(\boldsymbol{\theta},\boldsymbol{\delta})\triangleq\mathcal{L}^{*}(\boldsymbol{x}_{ini},\boldsymbol{\theta},\boldsymbol{\delta})$

Note that $\boldsymbol{\nu}_{q}^{*},\boldsymbol{\lambda}_{q}^{*}$ depend on
$\boldsymbol{\delta}_{q},\boldsymbol{x}_{ini}$. We make one key observation:

######  Proposition 3.

$\boldsymbol{x}_{q}^{0}$ depends on
$\boldsymbol{\nu}_{q}^{*},\boldsymbol{\lambda}_{q}^{*}$ but does not have
explicit dependency on $\boldsymbol{\delta}_{q},\boldsymbol{x}_{ini}$,
$\boldsymbol{\theta}$.

###### Proof.

This is true given (17). ∎

With Proposition 3, when $\boldsymbol{\delta}_{q},\boldsymbol{x}_{ini}$,
$\boldsymbol{\theta}$ change, $\boldsymbol{x}_{q}^{0}$ is still the
unconstrained minimizer of $\mathcal{L}$ as long as we do not swap
$\boldsymbol{\nu}_{q}^{*},\boldsymbol{\lambda}_{q}^{*}$. However,
$\boldsymbol{\nu}_{q}^{*},\boldsymbol{\lambda}_{q}^{*}$ are no longer the
maximizer of $g$. Hence,
$\mathcal{L}^{*}(\boldsymbol{x}_{ini},\boldsymbol{\theta},\boldsymbol{\delta})$
only provides a loose lower bound for $\boldsymbol{\delta}$ other than the
current $\boldsymbol{\delta}_{q}$ used to generate this optimality cut. The
subscript $q$ is dropped in (26) indicating that the inequality is valid for
general $\boldsymbol{\delta}$. As $\boldsymbol{x}_{ini}$ and
$\boldsymbol{\theta}$ take the same position as $\boldsymbol{\delta}_{q}$ in
$\mathcal{L}$, the same argument applies when $\boldsymbol{x}_{ini}$ and
$\boldsymbol{\theta}$ are updated. This will be used to construct warm-starts
for hybrid MPC.

If the solver used for the subproblem does not return the unconstrained
optimizer $\boldsymbol{x}_{q}^{0}$, we can leverage on strong duality to avoid
computing $\boldsymbol{x}_{q}^{0}$. Since
$v(\boldsymbol{x}_{ini},\boldsymbol{\theta},\boldsymbol{\delta}_{q})=f_{obj}(\boldsymbol{x}_{q}^{*})\triangleq
f_{obj,q}^{*}=\mathcal{L}^{*}(\boldsymbol{x}_{ini},\boldsymbol{\theta},\boldsymbol{\delta}_{q})$,
the cutting plane takes the form:

$\displaystyle
v(\boldsymbol{x}_{ini},\boldsymbol{\theta},\boldsymbol{\delta})\geq
f_{obj,q}^{*}$
$\displaystyle+\boldsymbol{\nu}_{q}^{*T}(\boldsymbol{b}(\boldsymbol{x}_{ini},\boldsymbol{\delta}_{q})-\boldsymbol{b}(\boldsymbol{x}_{ini},\boldsymbol{\delta}))$
(27)
$\displaystyle+\boldsymbol{\lambda}_{q}^{*T}(\boldsymbol{d}(\boldsymbol{\theta},\boldsymbol{\delta}_{q})-\boldsymbol{d}(\boldsymbol{\theta},\boldsymbol{\delta}))$

###  IV-C The Benders master problem

The final form of the master problem (14) is:

$\displaystyle\underset{\boldsymbol{\delta}}{\text{minimize}}\ \ z_{0}$ (28)
s.t. $\displaystyle\boldsymbol{\delta}_{k}\in\\{0,1\\}$ $\displaystyle for\
p=1,...,\text{current \\# of infeasible subproblem:}$
$\displaystyle\boldsymbol{b}(\boldsymbol{x}_{ini},\boldsymbol{\delta})^{T}\tilde{\boldsymbol{\nu}}_{p}+\boldsymbol{d}(\boldsymbol{\theta},\boldsymbol{\delta})^{T}\tilde{\boldsymbol{\lambda}}_{p}\geq
0$
$\displaystyle\boldsymbol{b}(\boldsymbol{x}_{ini},\boldsymbol{\delta})^{T}\tilde{\boldsymbol{\nu}}_{p}^{m}+\boldsymbol{d}(\boldsymbol{\theta},\boldsymbol{\delta})^{T}\tilde{\boldsymbol{\lambda}}_{p}^{m}\geq
0,\ \ m=1,...,N-1$ $\displaystyle for\ q=1,...,\text{current \\# of optimal
subproblem:}$ $\displaystyle z_{0}\geq
f_{obj,q}^{*}+\boldsymbol{\nu}_{q}^{*T}(\boldsymbol{b}(\boldsymbol{x}_{ini},\boldsymbol{\delta}_{q})-\boldsymbol{b}(\boldsymbol{x}_{ini},\boldsymbol{\delta}))+$
$\displaystyle\ \ \ \ \ \ \
\boldsymbol{\lambda}_{q}^{*T}(\boldsymbol{d}(\boldsymbol{\theta},\boldsymbol{\delta}_{q})-\boldsymbol{d}(\boldsymbol{\theta},\boldsymbol{\delta}))$

We define the epigraph variable $z_{0}$ and make $z_{0}\geq
v(\boldsymbol{x}_{ini},\boldsymbol{\theta},\boldsymbol{\delta}_{q})\ \forall
q$ to find the smallest value of optimality cuts corresponding to binary mode
$\boldsymbol{\delta}$. For our MPC problem with relatively small $N$ ($10\sim
30$), this problem is a small-scale MIP problem that contains $2N$ binary
variables and one continuous variable. As the algorithm proceeds, constraints
will be added to the master problem. For our test, the number of feasibility
cuts and optimality cuts has a scale of hundreds or thousands. The master
problem is solved within the presolve stage by an off-the-shelf MIP solver.

###  IV-D Upper bound and lower bound of the cost

Benders decomposition is an iterative procedure that adds feasibility cuts and
optimality cuts to improve the approximation of $v$ and $V$ in the master
problem. Since the optimality cuts are lower bounds of the actual cost $v$
according to (26), the master problem will "underestimate the difficulty" of
the subproblem when few cuts are generated. This means it will propose
$\boldsymbol{\delta}$ that is oftentimes infeasible, or with an optimal cost
cannot be reached by the subproblem. We define the master problem cost
$z_{0,i}$ at the latest iteration $i$ as the lower bound of the global optimal
cost. As more cuts are added, the lower bound increases to approach the actual
cost. Since at each iteration $i$, the MIP solver finds the global optimal
value of the master problem, the lower bound does not decrease with respect to
$i$ (if at next iteration $i+1$, $z_{0,i+1}<z_{0,i}$, $z_{0,i+1}$ should be
achievable at iteration $i$ since the problem has one additional constraint at
iteration $i+1$ than $i$, contradicting the fact that MIP solvers find the
global optimal solution).

On the other hand, the best subproblem cost till the current iteration
provides a upper bound of the global optimal cost. Since the accumulated cuts
in the master problem underestimate the actual subproblem cost, the cost of
the iteration $i+1$, $f^{*}_{obj,i+1}$, may not decrease comparing to the cost
$f^{*}_{obj,i}$ at iteration $i$, if the estimated cost is far from the actual
one. The upper bound is introduced to keep track of the best cost.

As the algorithm proceeds, the lower bound and upper bound will approach each
other. Eventually, they converge to each other and the algorithm terminates (a
proof of convergence can be found in [12]). We introduce the $MIPgap$ for the
termination condition. For fair comparison, we use the same definition as
Gurobi [54]:

$g_{a}=|z_{P}-z_{D}|/|z_{P}|$ (29)

Where $z_{P}$ is the primal bound, or upper bound of the global optimal cost,
and $z_{D}$ is the dual bound, or the lower bound of the cost. When $g_{a}$ is
smaller than a predefined threshold $G_{a}$, the inner loop algorithm
terminates with an optimal $\boldsymbol{\delta}^{*}$ and control signal
$\boldsymbol{u}^{*}$. We present the Benders MPC algorithm in Algorithm 1.

Input: $G_{a}$

1 Initialization $LB\coloneqq-\infty$, $UB\coloneqq\infty$, iteration
$i\coloneqq 0$

2 while  _$|UB-LB|/|UB|\geq G_{a}$_ do

3 Solve master problem (28) to get $\boldsymbol{\delta}$ and $m^{*}_{obj,i}$

4 Let $LB\coloneqq m^{*}_{obj,i}$

5 Solve problem (19) with $\boldsymbol{\delta}$ using dual simplex

6 if  _Feasible_ then

7 Solve (13) with solutions from (19) as warm-starts

8 Let the optimal cost be $f^{*}_{obj,i}$

9 if  _$f^{*}_{obj,i}$ <UB_ then

10 Let $UB\coloneqq f^{*}_{obj,i}$,
$\boldsymbol{u}^{*}\coloneqq\boldsymbol{u}$

11 Add constraint (26) to master problem (28)

12 else Infeasible

13 Add constraint (22) to master problem (28)

14 if  _$\tilde{\boldsymbol{\nu}}_{i}^{m}$ ,
$\tilde{\boldsymbol{\lambda}}_{i}^{m}$ is not equivalent to any existing
$\tilde{\boldsymbol{\nu}}$, $\tilde{\boldsymbol{\lambda}}$_ then

15 Add constraint (24) to master problem (28)

16

17 $i\coloneqq i+1$

return _$\boldsymbol{u}^{*}$_

Algorithm 1 Benders MPC

##  V Continual learning for warm-start

The idea behind the original Benders decomposition for linear programming
subproblems [30] is to express the solutions to the dual of the subproblem by
extreme rays and extreme points, and construct feasibility cuts and optimality
cuts based on them. If all extreme rays and extreme points are found for the
given subproblem, the master problem is guaranteed to find the global optimal
$\boldsymbol{\delta}$. For the QP subproblem in this paper, the optimal
solutions are not necessarily extreme points, but not far from them. In fact,
$||\boldsymbol{A}^{T}\boldsymbol{\nu}+\boldsymbol{C}^{T}\boldsymbol{\lambda}||$
is minimized in (18). If this norm is zero at the optimality, it is equivalent
to the linear dual feasible set of (20). As long as this norm is bounded, the
dual solutions $(\boldsymbol{\nu}^{*},\boldsymbol{\lambda}^{*})$ come from a
bounded (hence totally bounded) region around the extreme points of
$\boldsymbol{A}^{T}\boldsymbol{\nu}+\boldsymbol{C}^{T}\boldsymbol{\lambda}=0$,
which can be finitely covered by open balls
$\textit{B}_{\epsilon}(\boldsymbol{\nu}^{*},\boldsymbol{\lambda}^{*})$.
However, the number of extreme rays and covers is astronomical, hence an
iterative procedure is used to only add necessary cuts in searching for
optimal solutions. This avoids completely constructing the problem-solution
mappings such as explicit MPC [23], especially when the matrices
$\boldsymbol{A}$, $\boldsymbol{C}$ are undetermined before the solver begins
(for example, the robot may have an unknown payload until it is hand over).

In this paper, we extend this idea to continual learning with dynamic
environment represented by $\boldsymbol{\theta}$ shifting online. Only a small
number of extreme rays and covers are added for a given $\boldsymbol{\theta}$.
The solver continuously generates more rays and covers as
$\boldsymbol{\theta}$ shifts. The new rays and covers are retained, while the
duplicated ones are discarded. Once the solution is retained, it is never
removed. The retained solutions provide increasingly better warm-starts for
the incoming new problem instances. Since the number of extreme rays and
covers is finite, this provess will terminate and the master problem does not
grow infinitely large. This approach bares similarity of the continual
learning framework111A large number of literature continual learning is based
on tasks ([55, 56, 57], to name a few). On the other hand, this work does not
define tasks, hence more in-line with the task-free continual learning such as
[58]..

As MPC proceeds online, problem (III) needs to be constantly resolved with
different initial conditions $\boldsymbol{x}_{ini}$ and parameter
$\boldsymbol{\theta}$. Since $\boldsymbol{x}_{ini}$ and $\boldsymbol{\theta}$
take the same position in the subproblem (13) as $\boldsymbol{\delta}$, all
the optimality cuts (26) that construct lower bounds for $\boldsymbol{\delta}$
are also valid lower bounds for changing $\boldsymbol{x}_{ini}$ and
$\boldsymbol{\theta}$. In addition, the feasibility cuts (22) can also be used
for new initial condition as $\tilde{\boldsymbol{\nu}}_{p}$,
$\tilde{\boldsymbol{\lambda}}_{p}$ are independent of $\boldsymbol{x}_{ini}$
and $\boldsymbol{\theta}$. Assume we have feasible and optimality cuts as
listed in (28). When new initial condition $\boldsymbol{x}_{ini}^{\prime}$ and
different parameter $\boldsymbol{\theta}^{\prime}$ come in, we update the
cutting planes:

$\displaystyle for\ p=1,...,\text{current \\# of infeasible subproblem:}$ (30)
$\displaystyle\boldsymbol{b}(\boldsymbol{x}_{ini}^{\prime},\boldsymbol{\delta})^{T}\tilde{\boldsymbol{\nu}}_{p}+\boldsymbol{d}(\boldsymbol{\theta}^{\prime},\boldsymbol{\delta})^{T}\tilde{\boldsymbol{\lambda}}_{p}\geq
0$

$\displaystyle for\ q=1,...,$ current # of optimal subproblem (31)
$\displaystyle z_{0}\geq f_{obj,q}^{*}$
$\displaystyle+\boldsymbol{\nu}_{q}^{*T}(\boldsymbol{b}(\boldsymbol{x}_{ini},\boldsymbol{\delta}_{q})-\boldsymbol{b}(\boldsymbol{x}_{ini}^{\prime},\boldsymbol{\delta}))$
$\displaystyle+\boldsymbol{\lambda}_{q}^{*T}(\boldsymbol{d}(\boldsymbol{\theta},\boldsymbol{\delta}_{q})-\boldsymbol{d}(\boldsymbol{\theta}^{\prime},\boldsymbol{\delta}))$

######  Corollary 3.1.

For given $\boldsymbol{x}_{ini}^{\prime}$ and $\boldsymbol{\theta}^{\prime}$,
any $\boldsymbol{\delta}$ that contradicts (30) proves infeasibility for (19).

###### Proof.

A simple result from Lemma (2) given $\tilde{\boldsymbol{\nu}}_{p}$ and
$\tilde{\boldsymbol{\lambda}}_{p}$ are independent of $\boldsymbol{x}_{ini}$
and $\boldsymbol{\theta}$. ∎

######  Corollary 3.2.

$f_{obj,q}^{*}+\boldsymbol{\nu}_{q}^{*T}(\boldsymbol{b}(\boldsymbol{x}_{ini},\boldsymbol{\delta}_{q})-\boldsymbol{b}(\boldsymbol{x}_{ini}^{\prime},\boldsymbol{\delta}))+\boldsymbol{\lambda}_{q}^{*T}(\boldsymbol{d}(\boldsymbol{\theta},\boldsymbol{\delta}_{q})-\boldsymbol{d}(\boldsymbol{\theta}^{\prime},\boldsymbol{\delta}))$

gives a lower bound of
$v(\boldsymbol{x}_{ini}^{\prime},\boldsymbol{\theta}^{\prime},\boldsymbol{\delta})$.

###### Proof.

A simple result from that $\boldsymbol{x}_{ini}$, $\boldsymbol{\theta}$ and
$\boldsymbol{\delta}$ take the same position in (26). ∎

With this technique, when a new initial condition is received, we run the
updated master problem first. The master problem automatically provides a good
warm-start using knowledge of previously accumulated cuts, reducing the number
of iterations. The modified MPC algorithm with continual learning is provided
in Algorithm 2.

Input: $\boldsymbol{x}_{ini},\ \boldsymbol{\theta},\ G_{a},\ \epsilon$

1 Initialization $LB\coloneqq-\infty$, $UB\coloneqq\infty$, iteration
$i\coloneqq 0$

2 Update $\boldsymbol{x}_{ini}$, $\boldsymbol{\theta}$ in all the existing
feasibility and optimality cuts in master problem (28)

3 while  _$|UB-LB|/|UB|\geq G_{a}$_ do

4 Solve master problem (28) to get $\boldsymbol{\delta}$ and $m^{*}_{obj,i}$

5 Let $LB\coloneqq m^{*}_{obj,i}$

6 Solve problem (19) with $\boldsymbol{\delta}$ using dual simplex

7 if  _Feasible_ then

8 Solve (13) with solutions from (19) as warm-starts

9 Let the optimal cost be $f^{*}_{obj,i}$, optimal dual variables be
$\boldsymbol{\nu}_{i}^{*},\boldsymbol{\lambda}_{i}^{*}$

10 if
_$(\boldsymbol{\nu}_{i}^{*},\boldsymbol{\lambda}_{i}^{*})\notin\textit{B}_{\epsilon}(\boldsymbol{\nu}_{i}^{*},\boldsymbol{\lambda}_{i}^{*})$_
then

11 Add constraint (26) to master problem (28)

12 if  _$f^{*}_{obj,i}$ <UB_ then

13 Let $UB\coloneqq f^{*}_{obj,i}$,
$\boldsymbol{u}^{*}\coloneqq\boldsymbol{u}$

14

15

16 else Infeasible

17 Add constraint (22) to master problem (28)

18 if  _$\tilde{\boldsymbol{\nu}}_{i}^{m}$ ,
$\tilde{\boldsymbol{\lambda}}_{i}^{m}$ is not equivalent to any existing
$\tilde{\boldsymbol{\nu}}$, $\tilde{\boldsymbol{\lambda}}$_ then

19 Add constraint (24) to master problem (28)

20

21 $i\coloneqq i+1$

return _$\boldsymbol{u}^{*}$_

Algorithm 2 Benders MPC with continual learning

##  VI Experiment

We test our Benders MPC to control the inverted pendulum with moving soft
walls. This is also presented as a verification problem in previous works [9,
11, 18], except that we additionally randomize the wall motion.

Figure 1: Cart-pole system with moving soft contact walls.

####  VI-1 Problem setup

The setup is shown in Fig. 1. Let the nonlinear pendulum dynamics be
$\dot{\boldsymbol{x}}=f(\boldsymbol{x},\boldsymbol{u})+\boldsymbol{n}$.
$x_{1}$ is the position of the cart, $x_{2}$ is the angle of the pole, and
$x_{3}$, $x_{4}$ are their derivatives. The control input $\boldsymbol{u}$ is
the horizontal actuation force to push the cart. $\boldsymbol{n}$ is a random
disturbance torque acting on the pole. The moving elastic pads are located to
the right of the origin at a distance of $d_{1}$, and to the left at a
distance of $d_{2}$. Let $l$ be the length of the pole. When the pole
penetrates ($x_{1}-lcos(x_{2})\geq d_{1}$ or $x_{1}-lcos(x_{2})\leq-d_{2}$),
additional contact force is generated at the tip of the pole. Let the
parameter $\boldsymbol{\theta}=\begin{bmatrix}d_{1}&d_{2}\end{bmatrix}$.

We linearize the pendulum model around $x_{2}=0$ and use a linear elastic law
for the wall contact. The linear model is:

$\displaystyle\dot{x}_{1}$ $\displaystyle=x_{3}$ (32)
$\displaystyle\dot{x}_{2}$ $\displaystyle=x_{4}$ $\displaystyle\dot{x}_{3}$
$\displaystyle=\frac{gm_{p}}{m_{c}}x_{2}+\frac{u}{m_{c}}$
$\displaystyle\dot{x}_{4}$
$\displaystyle=\frac{g(m_{c}+m_{p})}{lm_{c}}x_{2}+\frac{u}{lm_{c}}+\frac{\lambda_{1}}{lm_{p}}-\frac{\lambda_{2}}{lm_{p}}$

Where $m_{p}$ is the mass of the pole, $m_{c}$ is the mass of the cart,
$\lambda_{1}$, $\lambda_{2}$ are contact forces from the right and the left
walls, respectively. They are both assumed to be positive. $g$ is the
gravitational acceleration. We define the control input
$\boldsymbol{u}=\begin{bmatrix}u&\lambda_{1}&\lambda_{2}\end{bmatrix}^{T}$. If
penetration happens, there is a non-zero contact force. This can be modeled as
mixed-integer linear constraints:

$\displaystyle\delta_{i}=0$ $\displaystyle\Rightarrow\lambda_{i}=0,\
a_{i}(lx_{2}-x_{1})+\frac{\lambda_{i}}{k_{i}}+d_{i}\geq 0$ (33)
$\displaystyle\delta_{i}=1$ $\displaystyle\Rightarrow\lambda_{i}\geq 0,\
a_{i}(lx_{2}-x_{1})+\frac{\lambda_{i}}{k_{i}}+d_{i}=0$

Where $i=1,2$. $a_{1}=1$ and $a_{2}=-1$. $k_{1}$ and $k_{2}$ are elastic
coefficients to model the right and left wall contacts. These logic laws are
enforced as mixed-integer linear constraints using the standard big-M approach
[59], where the maximal distance possible from pole to wall, $D_{max}$, and
maximal possible contact force, $\lambda_{max}$, are used as big-M constants.
We also define the maximal and minimal cart position limits $d_{min}$ and
$d_{max}$, and angle limits to be $\pm\frac{\pi}{2}$. The velocity limits of
the cart, angular velocity limits of the pole, and control limits $u_{max}$
are also defined accordingly. This problem has $n_{x}=4$, $n_{u}=3$,
$n_{z}=2$, $n_{c}=20$ (including variable limits). The matrices (variable
limits are excluded) after discretization are defined such that:

$\boldsymbol{E}=\boldsymbol{I}_{4}+dt\begin{bmatrix}0&0&1&0\\\ 0&0&0&1\\\
0&gm_{p}/m_{c}&0&0\\\ 0&g(m_{c}+m_{p})/(lm_{c})&0&0\end{bmatrix}$ (34)
$\boldsymbol{F}=dt\begin{bmatrix}0&0&0\\\ 0&0&0\\\ 1/m_{c}&0&0\\\
1/(lm_{c})&1/(lm_{p})&-1/(lm_{p})\end{bmatrix}$ (35)
$\boldsymbol{H}_{1}=\begin{bmatrix}0&0&0&0\\\ 0&0&0&0\\\ -1&l&0&0\\\
1&-l&0&0\\\ 1&-l&0&0\\\ -1&l&0&0\end{bmatrix}\
\boldsymbol{H}_{2}=\begin{bmatrix}0&1&0\\\ 0&0&1\\\ 0&1/k_{1}&0\\\
0&-1/k_{1}&0\\\ 0&0&1/k_{2}\\\ 0&0&-1/k_{2}\end{bmatrix}$ (36)
$\boldsymbol{H}_{3}=\begin{bmatrix}-\lambda_{max}&0\\\ 0&-\lambda_{max}\\\
D_{max}&0\\\ 0&0\\\ 0&D_{max}\\\ 0&0\end{bmatrix}\
\boldsymbol{h}(\boldsymbol{\theta})=\begin{bmatrix}0\\\ 0\\\ -d_{1}+D_{max}\\\
d_{1}\\\ -d_{2}+D_{max}\\\ d_{2}\end{bmatrix}$ (37)

Other matrices are zeros. The objective function penalizes the control
efforts, the velocities, and tries to regulate the pole to the zero position.
We choose $\boldsymbol{Q}_{k}=\boldsymbol{I}_{4}$,
$\boldsymbol{R}_{k}=\boldsymbol{I}_{2}$. The terminal cost
$\boldsymbol{Q}_{N}$ is obtained by solving a discrete algebraic Ricatti
equation. In the actual experiment, we choose $m_{c}=1.0kg$, $m_{p}=0.4kg$,
$l=0.6m$, $k_{1}=k_{2}=50N/m$, $u_{max}=20N$. The discretization step size
$dt=0.02s$ and planning horizon $N=10$.

####  VI-2 Monte-Carlo experiment

We implement Algorithm 2 to solve this problem. We choose $gap=0.1$, which is
identical among all the benchmark methods. The $\epsilon$ is chosen properly
to reduce the number of optimality cuts. We use off-the-shelf solver Gurobi to
solve both the master problems (MIPs) and the subproblems (QPs).

The controller is coded in Python, and tested inside a pybullet environment
[60] on a 12th Gen Intel Core i7-12800H × 20 laptop with 16GB memory. At the
beginning of each test episode, the pendulum begins from a perturbed angle of
$x_{2}=10\degree$ such that it will bump into the wall to regain balance. For
the rest of each episode, the persistent random disturbance torque
$\boldsymbol{n}$ is generated from a Gaussian distribution with zero mean and
a standard deviation $\boldsymbol{\sigma}=8Nm$. The system is constantly
disturbed and has to frequently touch the wall for re-balance. The wall motion
is generated by a random walk on top of a base sinusoidal motion with a
constant offset $d_{\text{off},i}$:

$d_{i}=d_{\text{off},i}+Asin(wt+\theta_{1})+m_{i},\ i=1,2$ (38)

Where $m_{i}$ is the integration of a Gaussian noise. We conduct statistical
analysis for 10 feasible trajectories under different disturbance torque
$\boldsymbol{n}$ and wall motion $m_{i}$. The data is collected from solved
problems where at least one contact is planned, removing the cases when
contact is not involved. The disturbance torque is unknown to the controller.
The wall motions are provide to controller only at run-time.

The following methods are also used for benchmark:

  1. (i)

Warm-started Branch and Bound. We implemented the Branch and Bound algorithm
in Python with warm-start as described by [11].

  2. (ii)

Off-the-shelf solver. We implemented the off-the-shelf solver Gurobi. The
problem is setup only once and solved iteratively such that warm-starts are
automatically used to minimize the solving time. The default setting is used
to optimize the performance.

  3. (iii)

GBD without warm-start. We implemented Algorithm 1 such that previous cuts are
not used to warm-start the next problem.

####  VI-3 Results

Fig. 2 gives the histogram result showing the number of iterations to solve
the problem and their frequencies. Thanks to the continual learning, 99.2% of
problem instances are solved within 5 iterations by GBD, except for a few
problems during the cold-start phase taking more than 10 iterations. This
99.2% problem instances have an average solving speed of $500-600Hz$. This
solving time represents a complete procedure of Algorithm 2, where the time
spend in solving master and subproblems account for 60% of the total time.
Note that previous work [40] reported over 90% solving time spent on the
master problem. On the contrary, we report that the master problem is
oftentimes solved within the presolve stage, since hundreds of cuts are
accumulated after a few problem instances. This accounts for less than 30% of
the total solving time.

On the other hand, Branch and Bound algorithm relies on subproblems. The warm-
start scheme reduces the number of solved subproblems by 50%. However, the BB
solver still goes through more than 10$\times$ subproblems to converge
compared to the GBD solver, from our averaged data. The gist of this warm-
start scheme is to shift the covers in time. Ideally, the computations up till
$k=N-2$ can be reused and the algorithm only needs to compute the binary input
at $k=N-1$. However, this scheme becomes less effective as the amplitude of
wall motion $A$ increases. The reason is that the current contact sequence
cannot simply be shifted and has to be recomputed before $k=N-1$.
Consequently, additional covers need to be refined.

Figure 2: Comparison of number of solver iterations for different problems of
$(\boldsymbol{x}_{ini},\boldsymbol{\theta})$. x-axis is the range of solver
iterations. y axis is the count of problem instances from the collected
trajectories. Left: The proposed GBD with continual learning. Right: Branch
and Bound with warm-start [11]. Figure 3: A case of solving procedure from
cold-start when the pole bumps into the moving elastic wall. x-axis is time
and y-axis is the solving speed in Hz. Left: Comparison of solving speed
between GBD with continual learning for warm-start, GBD without any warm-
start, and off-the-shelf solver Gurobi. Right: The number of cuts accumulated
during the solving procedure.

Fig. 3 shows the solving speed in Hz during the beginning of an episode from
our data. The solver begins from cold-start but has to plan contact ever since
$t=0$. After one iteration of from cold-start (taking $200ms$ in Fig. 3), the
cuts accumulate to provide warm-start for the next iteration, and solving
speed increases over Gurobi ($200-300Hz$). Without warm-start, the solving
speed remains on average $25Hz$.

Due to the fast cold-start, even if the system dynamics are only known at run-
time or completely change, the time cost to learn new dynamics for our problem
is at the scale of hundreds of milliseconds. This is much faster than training
neural-network-based policies [18]. If global optimal solutions are not
required in the beginning, the learning time can be further reduced.

##  VII Conclusion, Discussion and Future Work

In this paper, we proposed a hybrid MPC solver based on Generalized Benders
decomposition with continual learning. The algorithm accumulates cutting
planes from the invariant dual space of the subproblems under than randomly
changing environment. After a cold-start phase at the scale of hundreds of
milliseconds, the accumulated cuts provide warm-starts for the new problem
instances to increase the solving speed. This leads to solving speeds that are
2-3 times faster than the commercialized solver Gurobi in controlling the
cart-pole system with randomly moving soft walls.

There are several theoretical analyses and hardware experiments that can make
this preliminary results more thorough. For example, analyzing the
generalizability of the learned cuts to the new problem instances, or testing
the scalability of this algorithm to more complex problems. We can also
combine Branch and Bound with Benders cuts to leverage both of their
strengths. Although we already have results to some of the questions above,
they do not fit into the current paper and will be included in the future
journal version.

##  Appendix A Proof of lemma 1

We present a proof of Lemma 1. Recall the Farkas’ lemma (Theorem 4.6 in [51]):

######  Theorem 4.

Let $\tilde{\boldsymbol{A}}\in\mathbb{R}^{m\times n}$ and
$\tilde{\boldsymbol{b}}\in\mathbb{R}^{m}$. Then, exactly one of the two
following alternatives holds:

  1. 1.

There exists some $\tilde{\boldsymbol{x}}\geq\boldsymbol{0}$ such that
$\tilde{\boldsymbol{A}}\tilde{\boldsymbol{x}}=\tilde{\boldsymbol{b}}$.

  2. 2.

There exists some vector $\tilde{\boldsymbol{p}}$ such that
$\tilde{\boldsymbol{p}}^{T}\tilde{\boldsymbol{A}}\geq\boldsymbol{0}^{T}$ and
$\tilde{\boldsymbol{p}}^{T}\tilde{\boldsymbol{b}}<0$.

For any vector $\boldsymbol{x}$, there exists
$\boldsymbol{y}\geq\boldsymbol{0}$, $\boldsymbol{z}\geq\boldsymbol{0}$ such
that $\boldsymbol{x}=\boldsymbol{y}-\boldsymbol{z}$. The inequality constraint
$\boldsymbol{C}\boldsymbol{x}\leq\boldsymbol{d}$ is equivalent to
$\boldsymbol{C}\boldsymbol{x}+\boldsymbol{\delta}=\boldsymbol{d},\exists\boldsymbol{\delta}\geq\boldsymbol{0}$.
Hence the first condition of Lemma 1 is equivalent to the existence of
$\tilde{\boldsymbol{x}}=\begin{bmatrix}\boldsymbol{y}^{T}&\boldsymbol{z}^{T}&\boldsymbol{\delta}^{T}\end{bmatrix}^{T}\geq
0$ such that:

$\underbrace{\begin{bmatrix}\boldsymbol{A}&-\boldsymbol{A}&\boldsymbol{0}\\\
\boldsymbol{C}&-\boldsymbol{C}&\boldsymbol{I}\end{bmatrix}}_{\tilde{\boldsymbol{A}}}\underbrace{\begin{bmatrix}\boldsymbol{y}\\\
\boldsymbol{z}\\\
\boldsymbol{\delta}\end{bmatrix}}_{\tilde{\boldsymbol{x}}}=\underbrace{\begin{bmatrix}\boldsymbol{b}\\\
\boldsymbol{d}\end{bmatrix}}_{\tilde{\boldsymbol{b}}}$ (39)

By Theorem 4, this condition is alternative to the existence of
$\tilde{\boldsymbol{p}}=\begin{bmatrix}\boldsymbol{y}^{T}&\boldsymbol{z}^{T}\end{bmatrix}^{T}$
such that
$\tilde{\boldsymbol{p}}^{T}\tilde{\boldsymbol{A}}\geq\boldsymbol{0}^{T}$ and
$\tilde{\boldsymbol{p}}^{T}\tilde{\boldsymbol{b}}<0$, which gives the second
condition of Lemma 1.

Acknowledgements The author would like to thank Zehui Lu, Shaoshuai Mou, and
Yan Gu for their helpful discussions and suggestions.

## References

  * [1] S. Kuindersma, R. Deits, M. Fallon, A. Valenzuela, H. Dai, F. Permenter, T. Koolen, P. Marion, and R. Tedrake, “Optimization-based locomotion planning, estimation, and control design for the atlas humanoid robot,” _Autonomous robots_ , vol. 40, pp. 429–455, 2016.
  * [2] X. Lin, J. Zhang, J. Shen, G. Fernandez, and D. W. Hong, “Optimization based motion planning for multi-limbed vertical climbing robots,” in _2019 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. IEEE, 2019, pp. 1918–1925.
  * [3] F. R. Hogan and A. Rodriguez, “Reactive planar non-prehensile manipulation with hybrid model predictive control,” _The International Journal of Robotics Research_ , vol. 39, no. 7, pp. 755–773, 2020.
  * [4] J. Zhang, X. Lin, and D. W. Hong, “Transition motion planning for multi-limbed vertical climbing robots using complementarity constraints,” in _2021 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2021, pp. 2033–2039.
  * [5] A. Bemporad and M. Morari, “Control of systems integrating logic, dynamics, and constraints,” _Automatica_ , vol. 35, no. 3, pp. 407–427, 1999.
  * [6] W. Heemels, J. M. Schumacher, and S. Weiland, “Linear complementarity systems,” _SIAM journal on applied mathematics_ , vol. 60, no. 4, pp. 1234–1269, 2000.
  * [7] E. Sontag, “Nonlinear regulation: The piecewise linear approach,” _IEEE Transactions on automatic control_ , vol. 26, no. 2, pp. 346–358, 1981.
  * [8] W. P. Heemels, B. De Schutter, and A. Bemporad, “Equivalence of hybrid dynamical models,” _Automatica_ , vol. 37, no. 7, pp. 1085–1091, 2001.
  * [9] A. Aydinoglu, A. Wei, and M. Posa, “Consensus complementarity control for multi-contact mpc,” _arXiv preprint arXiv:2304.11259_ , 2023.
  * [10] S. L. Cleac’h, T. Howell, M. Schwager, and Z. Manchester, “Fast contact-implicit model-predictive control,” _arXiv preprint arXiv:2107.05616_ , 2021.
  * [11] T. Marcucci and R. Tedrake, “Warm start of mixed-integer programs for model predictive control of hybrid systems,” _IEEE Transactions on Automatic Control_ , vol. 66, no. 6, pp. 2433–2448, 2020.
  * [12] A. M. Geoffrion, “Generalized benders decomposition,” _Journal of optimization theory and applications_ , vol. 10, pp. 237–260, 1972.
  * [13] R. Deits, T. Koolen, and R. Tedrake, “Lvis: Learning from value function intervals for contact-aware robot controllers,” in _2019 International Conference on Robotics and Automation (ICRA)_. IEEE, 2019, pp. 7762–7768.
  * [14] A. Cauligi, P. Culbertson, B. Stellato, D. Bertsimas, M. Schwager, and M. Pavone, “Learning mixed-integer convex optimization strategies for robot planning and control,” in _2020 59th IEEE Conference on Decision and Control (CDC)_. IEEE, 2020, pp. 1698–1705.
  * [15] J. Tobajas, F. Garcia-Torres, P. Roncero-Sánchez, J. Vázquez, L. Bellatreche, and E. Nieto, “Resilience-oriented schedule of microgrids with hybrid energy storage system using model predictive control,” _Applied Energy_ , vol. 306, p. 118092, 2022.
  * [16] Y. Cao, Z. Zhang, F. Cheng, and S. Su, “Trajectory optimization for high-speed trains via a mixed integer linear programming approach,” _IEEE Transactions on Intelligent Transportation Systems_ , vol. 23, no. 10, pp. 17 666–17 676, 2022.
  * [17] M. Klaučo, J. Drgoňa, M. Kvasnica, and S. Di Cairano, “Building temperature control by simple mpc-like feedback laws learned from closed-loop data,” _IFAC Proceedings Volumes_ , vol. 47, no. 3, pp. 581–586, 2014.
  * [18] A. Cauligi, P. Culbertson, E. Schmerling, M. Schwager, B. Stellato, and M. Pavone, “Coco: Online mixed-integer control via supervised learning,” _IEEE Robotics and Automation Letters_ , vol. 7, no. 2, pp. 1447–1454, 2021\.
  * [19] T. Marcucci, R. Deits, M. Gabiccini, A. Bicchi, and R. Tedrake, “Approximate hybrid model predictive control for multi-contact push recovery in complex environments,” in _2017 IEEE-RAS 17th International Conference on Humanoid Robotics (Humanoids)_. IEEE, 2017, pp. 31–38.
  * [20] A. M. Geoffrion, “Lagrangian relaxation for integer programming,” _50 Years of Integer Programming 1958-2008: From the Early Years to the State-of-the-Art_ , pp. 243–281, 2010.
  * [21] M. Anitescu, “On using the elastic mode in nonlinear programming approaches to mathematical programs with complementarity constraints,” _SIAM Journal on Optimization_ , vol. 15, no. 4, pp. 1203–1236, 2005.
  * [22] R. Lazimy, “Mixed-integer quadratic programming,” _Mathematical Programming_ , vol. 22, pp. 332–349, 1982.
  * [23] A. Bemporad, M. Morari, V. Dua, and E. N. Pistikopoulos, “The explicit linear quadratic regulator for constrained systems,” _Automatica_ , vol. 38, no. 1, pp. 3–20, 2002.
  * [24] J.-J. Zhu and G. Martius, “Fast non-parametric learning to accelerate mixed-integer programming for hybrid model predictive control,” _IFAC-PapersOnLine_ , vol. 53, no. 2, pp. 5239–5245, 2020.
  * [25] X. Lin, G. I. Fernandez, and D. W. Hong, “Reduce: Reformulation of mixed integer programs using data from unsupervised clusters for learning efficient strategies,” in _2022 International Conference on Robotics and Automation (ICRA)_. IEEE, 2022, pp. 4459–4465.
  * [26] X. Lin, F. Xu, A. Schperberg, and D. Hong, “Learning near-global-optimal strategies for hybrid non-convex model predictive control of single rigid body locomotion,” _arXiv preprint arXiv:2207.07846_ , 2022.
  * [27] P. Hespanhol, R. Quirynen, and S. Di Cairano, “A structure exploiting branch-and-bound algorithm for mixed-integer model predictive control,” in _2019 18th European Control Conference (ECC)_. IEEE, 2019, pp. 2763–2768.
  * [28] B. Wu and B. Ghanem, “$\ell_{p}$-box admm: A versatile framework for integer programming,” _IEEE transactions on pattern analysis and machine intelligence_ , vol. 41, no. 7, pp. 1695–1708, 2018.
  * [29] B. Stellato, G. Banjac, P. Goulart, A. Bemporad, and S. Boyd, “Osqp: An operator splitting solver for quadratic programs,” _Mathematical Programming Computation_ , vol. 12, no. 4, pp. 637–672, 2020.
  * [30] J. Benders, “Partitioning procedures for solving mixed-variables programming problems ‘,” _Numerische mathematik_ , vol. 4, no. 1, pp. 238–252, 1962\.
  * [31] G. B. Dantzig and P. Wolfe, “Decomposition principle for linear programs,” _Operations research_ , vol. 8, no. 1, pp. 101–111, 1960.
  * [32] L. J. Watters, “Reduction of integer polynomial programming problems to zero-one linear programming problems,” _Operations Research_ , vol. 15, no. 6, pp. 1171–1174, 1967.
  * [33] F. Glover, “Improved linear integer programming formulations of nonlinear integer problems,” _Management science_ , vol. 22, no. 4, pp. 455–460, 1975\.
  * [34] R. McBride and J. Yormark, “An implicit enumeration algorithm for quadratic integer programming,” _Management Science_ , vol. 26, no. 3, pp. 282–296, 1980.
  * [35] J. N. Hooker and H. Yan, “Logic circuit verification by benders decomposition,” _Principles and practice of constraint programming: the newport papers_ , pp. 267–288, 1995.
  * [36] J. N. Hooker and G. Ottosson, “Logic-based benders decomposition,” _Mathematical Programming_ , vol. 96, no. 1, pp. 33–60, 2003.
  * [37] G. Codato and M. Fischetti, “Combinatorial benders’ cuts for mixed-integer linear programming,” _Operations Research_ , vol. 54, no. 4, pp. 756–766, 2006.
  * [38] J. R. Birge and F. Louveaux, _Introduction to stochastic programming_. Springer Science & Business Media, 2011.
  * [39] P.-D. Moroşan, R. Bourdais, D. Dumur, and J. Buisson, “A distributed mpc strategy based on benders’ decomposition applied to multi-source multi-zone temperature regulation,” _Journal of Process Control_ , vol. 21, no. 5, pp. 729–737, 2011.
  * [40] R. Rahmaniani, T. G. Crainic, M. Gendreau, and W. Rei, “The benders decomposition algorithm: A literature review,” _European Journal of Operational Research_ , vol. 259, no. 3, pp. 801–817, 2017.
  * [41] W. Rei, J.-F. Cordeau, M. Gendreau, and P. Soriano, “Accelerating benders decomposition by local branching,” _INFORMS Journal on Computing_ , vol. 21, no. 2, pp. 333–345, 2009.
  * [42] A. M. Costa, J.-F. Cordeau, B. Gendron, and G. Laporte, “Accelerating benders decomposition with heuristicmaster problem solutions,” _Pesquisa Operacional_ , vol. 32, pp. 03–20, 2012.
  * [43] T. L. Magnanti and R. T. Wong, “Accelerating benders decomposition: Algorithmic enhancement and model selection criteria,” _Operations research_ , vol. 29, no. 3, pp. 464–484, 1981.
  * [44] J.-F. Cordeau, G. Stojković, F. Soumis, and J. Desrosiers, “Benders decomposition for simultaneous aircraft routing and crew scheduling,” _Transportation science_ , vol. 35, no. 4, pp. 375–388, 2001.
  * [45] J. A. Rodríguez, M. F. Anjos, P. Côté, and G. Desaulniers, “Accelerating benders decomposition for short-term hydropower maintenance scheduling,” _European Journal of Operational Research_ , vol. 289, no. 1, pp. 240–253, 2021.
  * [46] A. Grothey, S. Leyffer, and K. McKinnon, “A note on feasibility in benders decomposition,” _Numerical Analysis Report NA/188, Dundee University_ , 1999\.
  * [47] M. V. Pereira and L. M. Pinto, “Stochastic optimization of a multireservoir hydroelectric system: A decomposition approach,” _Water resources research_ , vol. 21, no. 6, pp. 779–792, 1985.
  * [48] ——, “Multi-stage stochastic optimization applied to energy planning,” _Mathematical programming_ , vol. 52, pp. 359–375, 1991.
  * [49] J. Warrington, P. N. Beuchat, and J. Lygeros, “Generalized dual dynamic programming for infinite horizon problems in continuous state and action spaces,” _IEEE Transactions on Automatic Control_ , vol. 64, no. 12, pp. 5012–5023, 2019.
  * [50] S. Menta, J. Warrington, J. Lygeros, and M. Morari, “Learning q-function approximations for hybrid control problems,” _IEEE Control Systems Letters_ , vol. 6, pp. 1364–1369, 2021.
  * [51] D. Bertsimas and J. N. Tsitsiklis, _Introduction to linear optimization_. Athena scientific Belmont, MA, 1997, vol. 6.
  * [52] L. Wu and M. Shahidehpour, “Accelerating the benders decomposition for network-constrained unit commitment problems,” _Energy Systems_ , vol. 1, no. 3, pp. 339–376, 2010.
  * [53] N. Beheshti Asl and S. MirHassani, “Accelerating benders decomposition: multiple cuts via multiple solutions,” _Journal of Combinatorial Optimization_ , vol. 37, pp. 806–826, 2019.
  * [54] G. Optimization. (2023) Mipgap. [Online]. Available: https://www.gurobi.com/documentation/9.5/refman/mipgap2.html
  * [55] S. Thrun, “Is learning the n-th thing any easier than learning the first?” _Advances in neural information processing systems_ , vol. 8, 1995.
  * [56] S. Thrun and T. M. Mitchell, “Lifelong robot learning,” _Robotics and autonomous systems_ , vol. 15, no. 1-2, pp. 25–46, 1995.
  * [57] P. Ruvolo and E. Eaton, “Ella: An efficient lifelong learning algorithm,” in _International conference on machine learning_. PMLR, 2013, pp. 507–515.
  * [58] R. Aljundi, K. Kelchtermans, and T. Tuytelaars, “Task-free continual learning,” in _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , 2019, pp. 11 254–11 263.
  * [59] J. P. Vielma, “Mixed integer linear programming formulation techniques,” _Siam Review_ , vol. 57, no. 1, pp. 3–57, 2015.
  * [60] E. Coumans and Y. Bai, “Pybullet, a python module for physics simulation for games, robotics and machine learning,” 2016.