A Structure Exploiting Branch-and-Bound Algorithm for Mixed-Integer Model Predictive Control

TL;DR

This paper introduces a specialized branch-and-bound algorithm for mixed-integer model predictive control that leverages problem structure and warm-starting to improve computational efficiency in real-time applications.

Contribution

It presents a novel structure-exploiting branch-and-bound algorithm tailored for MI-MPC, incorporating warm-start strategies based on previous solutions.

Findings

Enhanced computational performance over existing solvers

Effective exploitation of problem structure and warm-starting

Successful application to automotive and aerospace case studies

Abstract

Mixed-integer model predictive control (MI-MPC) requires the solution of a mixed-integer quadratic program (MIQP) at each sampling instant under strict timing constraints, where part of the state and control variables can only assume a discrete set of values. Several applications in automotive, aerospace and hybrid systems are practical examples of how such discrete-valued variables arise. We utilize the sequential nature and the problem structure of MI-MPC in order to provide a branch-and-bound algorithm that can exploit not only the block-sparse optimal control structure of the problem but that can also be warm started by propagating information from branch-and-bound trees and solution paths at previous time steps. We illustrate the computational performance of the proposed algorithm and compare against current state-of-the-art solvers for multiple MPC case studies, based on a preliminary implementation in MATLAB and C code.

Figures

Tables

Table 1. TABLE I: Timing results (ms) per sampling step of hybrid MPC test problem for different horizon lengths N 𝑁 N . Computation times for MI-NNLS solver are taken directly from  [ 14 ] .

N	BB-PRESAS	GUROBI	MOSEK	MI-NNLS
	(mean/max)	(mean/max)	(mean/max)	(mean/max)
2	0.1/0.2	0.7/1.4	2.1/4.0	2.0/2.6
3	0.2/0.3	1.0/2.3	15.1/24.7	2.5/4.8
4	0.4/0.9	1.7/4.6	21.7/35.5	3.1/6.9
5	0.9/1.7	2.5/4.9	28.7/39.3	3.9/13.0
6	1.5/3.5	3.2/7.5	36.8/58.8	5.1/18.3
7	2.3/4.9	4.0/6.9	51.8/109.3	6.4/30.2
8	3.5/7.6	5.1/10.0	70.4/185.8	8.1/43.4
9	5.1/10.3	6.6/12.5	98.7/347.1	11.1/69.8
10	6.8/14.3	8.4/16.1	126.7/465.3	14.4/103.2
11	8.8/22.1	9.8/17.2	168.2/587.8	20.6/179.1
12	11.3/23.7	11.6/20.5	219.2/765.0	26.9/263.4
13	15.0/31.6	14.3/29.5	276.3/996.0	35.5/384.9
14	17.8/35.1	16.4/44.6	334.1/1241.9	46.3/562.4
15	21.0/41.6	21.9/71.6	450.8/1606.8	61.7/766.9

Equations

$$\underset{X,U}{\min }$$

$$x_0-{\hat{x}}_0=0,$$

$$A_i{x}_i+B_i{u}_i+a_i=x_{i+1},$$

$$l_i^c\le C_i{x}_i+D_i{u}_i\le u_i^c,$$

$$F_i{u}_i\in \{0,1\},$$

$$l_N^c\le C_N{x}_N\le u_N^c,$$

$$S_{i,k}=\text{score}({\Delta }_{i,k}^{-},{\Delta }_{i,k}^{+})=\text{max}({\Delta }_{i,k}^{+},ϵ)\cdot \text{max}({\Delta }_{i,k}^{-},ϵ),$$

$$S_{i,k}=\text{score}({\bar{u}}_{i,k}{\varphi }_{i,k}^{-},(1-{\bar{u}}_{i,k}){\varphi }_{i,k}^{+}),$$

$$\bar{A}X=\bar{B}U+b+E_0{\hat{x}}_0,$$

$$\bar{A}$$

$$\bar{B}$$

$$X={\bar{A}}^{-1}\bar{B}U+{\bar{A}}^{-1}(b+E_0{\hat{x}}_0).$$

$$\underset{U}{\min }$$

$${\bar{l}}^c\le D_{\mathrm{c}}U\le {\bar{u}}^c$$

$$F_i{u}_i\in \{0,1\},$$

$$H_{\mathrm{c}}$$

$$h_{\mathrm{c}}$$

$${\bar{l}}^c$$

$$d_i{u}_i\le u_b-\sum _{j\ne i}{d}_j{u}_j\le \underset{=:{\bar{u}}_{b,i}}{\underbrace{u_b-\sum _{j\ne i,d_j>0}{d}_j{\underline{u}}_j-\sum _{j\ne i,d_j<0}{d}_j{\bar{u}}_j}},$$

$$u_i\le \frac{{\bar{u}}_{b,i}}{d_i},\text{ if }{d}_i>0\text{or}{u}_i\ge \frac{{\bar{u}}_{b,i}}{d_i},\text{ if }{d}_i<0.$$

$${\bar{u}}_i=\min ({\bar{u}}_i,\frac{{\bar{u}}_{b,i}}{d_i}),\text{or}{\underline{u}}_i=\max ({\underline{u}}_i,\frac{{\bar{u}}_{b,i}}{d_i}),$$

$${\bar{u}}_i=\min ({\bar{u}}_i,\lfloor \frac{{\bar{u}}_{b,i}}{d_i}\rfloor ),\text{or}{\underline{u}}_i=\max ({\underline{u}}_i,\lceil \frac{{\bar{u}}_{b,i}}{d_i}\rceil ).$$

Full text

I INTRODUCTION

Optimization based control and estimation techniques, such as model predictive control (MPC) and moving horizon estimation (MHE), allow a model-based design framework in which the system dynamics and constraints can directly be taken into account [1]. This framework can be further extended to hybrid systems [2], providing a powerful technique to model a large range of problems, e.g., including dynamical systems with mode switchings or quantized control, problems with logic rules or no-go zone constraints. However, the resulting optimization problems are highly non-convex because they contain variables that only take integer values. When using a quadratic objective in combination with linear system dynamics and linear inequality constraints, the resulting optimal control problem (OCP) can be formulated as a mixed-integer quadratic program (MIQP).

We aim to solve MIQP problems of the following form:

[TABLE]

where the optimization variables are the state $X=[x_{0}^{\top},\ldots,x_{N}^{\top}]^{\top}$ and control trajectory $U=[u_{0}^{\top},\ldots,u_{N-1}^{\top}]^{\top}$. The set of constraints (1e) are binary equality constraints, since the left-hand side needs to be equal to either [math] or $1$. For simplicity of notation, we further consider only binary control variables instead of more general integer constraints for an affine function of both state and control variables. MPC for several classes of hybrid systems can be straightforwardly formulated as in (1). Notable examples are mixed logical systems [2], where auxiliary continuous and discrete variables can be added to the input vector. Moreover, in combination with the binary constraints (1e), the affine inequalities (1d) can model various complicated but practical restrictions on the feasible region, such as no-go zones and disjoint polyhedral constraints for states and inputs.

A hybrid MPC controller aims to solve the MIQP (1) at every sampling time instant. This is a difficult task, given that mixed-integer programming is $\mathcal{NP}$-hard in general, and several methods for solving such a sequence of MIQPs have been explored in the literature. These approaches can be divided into heuristic techniques, which seek to efficiently find sub-optimal solutions to the problem, and optimization algorithms which attempt to solve the MIQPs to optimality. Examples of the former include rounding and pumping schemes [3, 4], approximate optimization algorithms [5, 6], and approximate dynamic programming [7]. The downside of fast heuristic approaches is often the lack of guarantees for finding an optimal or even an integer-feasible solution. Heuristic rounding-based approaches to mixed-integer nonlinear OCPs can be found, e.g., in [8, 9].

As for solving these problems to optimality, most of the optimization algorithms for MIQPs are based on the classical branch-and-bound (B&B) technique [10]. For the purpose of mixed-integer MPC, the standard B&B strategy has been combined with various methods for solving the relaxed convex QPs. For example, a B&B algorithm for mixed-integer MPC (MI-MPC) has been proposed in combination with a dual active-set solver in [11], with an interior point algorithm in [12], dual projected gradient methods in [13, 6], a nonnegative least squares solver in [14], and the alternating direction method of multipliers (ADMM) in [15]. Branch-and-bound methods for solving mixed-integer nonlinear OCPs have also been studied, e.g., in [16].

Another important research topic focuses on general pre-processing and modeling techniques to reduce the size and strengthen the mixed-integer problem formulations [17]. These presolve techniques are vital to the good performance of current state-of-the-art mixed-integer solvers [18], such that these methods can often solve seemingly intractable problems in practice. Lastly, the branch-and-bound method itself has been extensively studied with several improvements in branching and variable selection techniques [19, 20], including recent developments in applying machine learning techniques in order to learn “better” branching rules [21]. Finally, the branch-and-bound strategy has been generalized further, e.g., using cutting planes to tighten the convex problem relaxations, resulting in branch-and-cut or branch-and-price variants of the algorithm [10, 17]. Unlike state-of-the-art mixed-integer solvers, e.g., GUROBI [22] and MOSEK [23], our aim is to propose a tailored algorithm and its solver implementation for fast embedded MI-MPC applications, i.e., running on microprocessors with considerably less computational resources and available memory. The optimization algorithm should be relatively simple to code with a moderate use of resources, while the software implementation is preferably compact and library independent.

In this paper, our first contribution is to propose a branch-and-bound based MPC algorithm, which exploits the features of a recently proposed structure-exploiting primal active-set solver called PRESAS [24]. The latter algorithm is tailored to efficiently solve QPs with a block-sparse optimal control structure. Our second contribution is to bring various mixed-integer programming techniques, such as bound strengthening, domain propagation, and advanced branching rules, to the context of MI-MPC. In particular, we present an algorithm that exploits the sequential nature of MPC, in order to warm-start the branch-and-bound search tree and to re-use information gathered at previous time steps. A similar type of approach was proposed recently by [14], but in this work we provide not only a warm-start procedure for the integer variables but we also show how to improve the branching strategy by warm starting and how to efficiently combine this with presolving techniques for MI-MPC. Finally, the computational performance of the proposed algorithm, for a preliminary implementation in Matlab and C code, is illustrated and compared against current state-of-the-art solvers for multiple MPC case studies.

The paper is organized as follows. Section II presents the basic idea of branch-and-bound methods for mixed-integer programming. Then, Section III presents presolve techniques in the context of mixed-integer optimal control. The resulting MI-MPC algorithm and its tailored warm-starting strategies are discussed in Section IV and its performance is illustrated based on multiple numerical case studies in Section V. Finally, Section VI concludes the paper.

II MIXED-INTEGER QUADRATIC PROGRAMMING

We first introduce some of the basic concepts in mixed-integer programming based on branch-and-bound methods, such as convex relaxations and branching strategies.

II-A Convex Quadratic Program Relaxations

A standard approach to solve the MIQP (1) is to create convex relaxations of this problem (by either dropping some constraints or by re-formulating the problem and providing an approximation scheme) and then solve the relaxations in order to approach the solution to the original MIQP. A straightforward idea is to obtain convex QP relaxations by dropping the binary equality constraints (1e) and instead enforcing the affine inequality constraints $0\leq F_{i}u_{i}\leq 1$. Other convex relaxations for MIQPs have been studied in the literature such as moment or SDP relaxations that are often tighter than QP relaxations [25, 26], but they can be relatively expensive to solve for larger problems.

For the purpose of this paper, we will focus our attention on QP relaxations where we allow the binary variables to take on real values. The main reason for choosing this relaxation is that we utilize a tailored structure exploiting active-set solver, called PRESAS [24], proposed recently for efficiently solving the convex QP relaxations. The latter solver has been shown to be competitive with state-of-the-art QP solvers for embedded MPC, and it benefits strongly from warm-starting, which can be exploited when solving the sequence of QPs within the branch-and-bound strategy. Note that the relaxations need to be convex, i.e., the weight matrices $Q_{i}$, $R_{i}$ and $P$ need to be positive (semi-) definite in (1a) such that each solution to a QP relaxation is globally optimal.

II-B Branch-and-Bound Algorithm

The main idea of the branch-and-bound (B&B) algorithm is to sequentially create partitions of the original problem and then attempt to solve those partitions. While solving each partition may still be challenging, it is fairly efficient to obtain local lower bounds on the optimal objective value, by solving relaxations of the mixed-integer program or by using duality. If we happen to obtain an integer-feasible solution while solving a relaxation, we can then use it to obtain a global upper bound for the solution to the original problem. This may help to avoid solving or branching certain partitions that were already created, i.e., these partitions or nodes can be pruned. The general algorithmic idea of partitioning is better illustrated as a binary search tree, see Figure 1.

A key step in this approach is how to create the partitions, i.e., which node to choose and which binary variable to select for branching. Since we solve a QP relaxation at every node of the tree, it is natural to branch on one of the binary variables with fractional values in the optimal solution of the QP relaxation. Therefore, if a variable, e.g., $u_{i,k}\in\{0,1\}$ has a fractional value in a given QP relaxation, then we create two partitions where we respectively add the equality constraint $u_{i,k}=0$ and $u_{i,k}=1$. Another key step is how to choose the order in which the created subproblems are solved. These two steps have been extensively explored in the literature and various heuristics are implemented in state-of-the-art tools [19]. We provide next a brief description of strategies that we implemented in our B&B solver.

II-C Tree Search: Node Selection Strategies

A common implementation of the branch-and-bound method is based on a depth-first node selection strategy, which can be readily implemented using a last-in-first-out (LIFO) buffer. The next node to be solved is selected as one of the children of the current node and this process is repeated until a node is pruned, i.e., the node is either infeasible, optimal or dominated by the upper bound, which is followed by a backtracking procedure. Instead, a best-first strategy selects the node with the lowest local lower bound so far. In what follows, we will employ a combination of the depth-first and best-first node selection approach. This idea is motivated by aiming to find an integer-feasible solution quickly at the start of the branch-and-bound procedure (depth-first) to allow for early pruning, followed by a more greedy search for better feasible solutions (best-first).

II-D Reliability Branching for Variable Selection

The idea of reliability branching is to combine two powerful concepts for variable selection: strong branching and pseudo-costs [19]. Strong branching relies on temporarily branching, both up (to higher integer) and down (to lower integer), for every binary variable that has a fractional value in the solution of a QP relaxation in a given node, before committing to the variable that provides the highest value for a particular score function. The increase in objective values $\Delta_{i,k}^{+}$, $\Delta_{i,k}^{-}$ are computed when branching the binary variable $u_{i,k}$, respectively, up and down. Given these quantities, a simple scoring function $\text{score}(\cdot,\cdot)$ is computed for each binary variable. For instance, based on the product [20]:

[TABLE]

given a small positive value $\epsilon>0$. This branching rule has been empirically shown to provide smaller search trees in practice [19]. The downside is that this procedure is relatively expensive since several QP relaxations are solved in order to select one variable to branch on.

The idea of pseudo-costs aims at approximating the increase of the objective function to decide which variable to branch on, without having to solve additional QP relaxations. This can be done by keeping statistic information for each binary variable, i.e., the pseudo-costs that represent the average increase in the objective value per unit change in that particular binary variable when branching. Every time that a given variable is chosen to be branched on, and the resulting relaxation is feasible, then we update each corresponding pseudo-cost with the observed increase in the objective, divided by the distance of the real to the binary value, in the form of a cumulative average. Therefore, each variable has two pseudo-costs, $\phi_{i,k}^{-}$ when the variable was branched “down” and $\phi_{i,k}^{+}$ when it was branched “up”. Given the solution to a QP relaxation, one can then use the pseudo-costs to select the binary variable with the highest score value to be branched on next:

[TABLE]

given a fractional value $\bar{u}_{i,k}$ in the QP relaxation.

This way, we select variables based on their past behavior throughout the branch-and-bound tree. However, at the beginning of the algorithm, the pseudo-costs are not yet initialized, which is when branching decisions typically impact the tree size the most. Reliability branching uses strong branching to initialize the pseudo-costs until a certain condition of reliability is satisfied, e.g., one switches to using pseudo-costs only once that particular variable has been branched on a specified number $\eta_{rel}$ of times [19]. The resulting branching rule is summarized in Algorithm 1. Note that reliability branching coincides with pseudo-cost branching if $\eta_{rel}=0$, with strong branching if $\eta_{rel}=\infty$, but typically a value $1\leq\eta_{rel}\leq 4$ is chosen.

This rule can be further augmented by implementing a look ahead limit in the number of candidates, as well as a limit in the number of iterations for each QP relaxation in the strong branching step. Note that many other branching rules exist such as, e.g., “most infeasible” branching which selects the binary variable with fractional part that is closest to $0.5$. Even though the latter rule is used quite often, e.g., in [14], it generally does not perform very well in practice [19]. Extensive empirical experiments with different branching strategies are beyond the scope of this paper.

III PRESOLVE TECHNIQUES FOR MIXED-INTEGER OPTIMAL CONTROL

As mentioned earlier, presolve techniques are often crucial in making convex relaxations tighter such that typically fewer nodes need to be explored, sometimes to such an extent that seemingly intractable problems become tractable. Next, we briefly describe some of these concepts with a focus on domain propagation for bound strengthening and its implementation for mixed-integer optimal control.

III-A Domain Propagation for Condensed QP Subproblem

Several strengthening techniques are implemented as part of “presolve” routines in commercial solvers [18]. One particular technique that is suitable to mixed-integer optimal control is based on domain propagation, in which the goal is to strengthen bound values based on the inequality constraints (1d)-(1f) in the problem. However, the results of such a strategy are rather weak when directly applied to the block-sparse QP in (1), because the stage-wise coupling of the state variables (1c) needs to be taken into account. Therefore, we use instead the equivalent dense QP formulation in which the state variables are numerically eliminated, such that stronger bounds can be obtained for the control variables. Hence, we can use the block-structured sparsity to efficiently solve the QP relaxations, while we use the equivalent but dense format to effectively perform domain propagation.

Let us concatenate all state variables in a vector $X$ and all control variables in the vector $U$, such that Eqs. (1b)-(1c) can be written more compactly as

[TABLE]

where we define the block-sparse matrices

[TABLE]

The matrix $\bar{A}$ is invertible such that we can write:

[TABLE]

Now, we can substitute the latter expression for the state vector in OCP (1) to obtain the condensed form

[TABLE]

where the condensed matrices and vectors read as

[TABLE]

where $\bar{b}:=b+E_{0}\hat{x}_{0}$ is defined and given $Q=\text{blkdiag}(Q_{1},\ldots,Q_{N-1},P)$ , $R=\text{blkdiag}(R_{0},\ldots,R_{N-1})$, and $l^{c}=[l^{c^{\top}}_{1},\ldots,l^{c^{\top}}_{N}]^{\top}$ and $u^{c}=[u^{c^{\top}}_{1},\ldots,u^{c^{\top}}_{N}]^{\top}$.

Given the condensed problem formulation, which can be computed offline and which is parametric in the current state value $\hat{x}_{0}$, we can then apply the following bound strengthening procedure, which is explained next for a single affine constraint $l_{b}\leq\sum_{i}d_{i}u_{i}\leq u_{b}$ in (7b). This constraint can be used to try and tighten bound values for all control variables $u_{i}$ for which $d_{i}\neq 0$, where $u_{i}$ denotes a single control variable in the vector $U$. Let $\bar{u}_{i},\underline{u}_{i}$ be the current upper/lower bounds for $u_{i}$ such that

[TABLE]

in which we divide by $d_{i}$ in order to obtain

[TABLE]

This results, respectively, in the updated bound values

[TABLE]

or, in case $u_{i}$ is an integer or binary variable,

[TABLE]

where $\lfloor\cdot\rfloor$ and $\lceil\cdot\rceil$ are the floor and ceiling operations, respectively. Thus, this can result in strengthening of bound values for both continuous and integer/binary control variables. The procedure can be executed for each control variable and each inequality constraint in an iterative manner, see Algorithm 2, since bound strengthening for one variable can lead to strengthening for other variables [18]. The process is typically stopped when the bound values do not sufficiently change or a certain limit on the computation time is met.

Domain propagation can lead to considerable reductions in the amount of explored nodes, e.g., because variables are fixed, when $\bar{u}_{i}=\underline{u}_{i}$, or because of infeasibility detection, when $\bar{u}_{i}<\underline{u}_{i}$, without the need to solve any QP relaxations. In addition, the updated bound values for all control variables can be used to strengthen QP relaxations in the future. Lastly, we can use domain propagation in order to improve and generalize Hessian-based fixing strategies, such as the one proposed in [27]. Hessian-based fixing typically can only be applied to unconstrained problems, since it fixes the variables solely based on the objective. Here, we propose to use domain propagation to compute the feasibility impact of certain variable fixings. More specifically, a particular variable can be fixed based on optimality, if and only if this fixing does not induce feasibility-based fixings.

III-B Probing Strategies and Cutting Planes

Probing [28] is a classical technique that can be incorporated in any branch-and-bound method to derive stronger inequalities or better bounds. It consists of tentatively trying to fix some variables and to derive potential logical implications on other variables. We do not further describe probing strategies in detail, but we refer to [18] for an overview. The computational cost and performance of probing can be greatly improved by relying on some of the other techniques that were discussed earlier. For example, the pseudo-costs can be used in order to choose the bound value for each binary variable that is likely to result in a low objective value. In turn, the QP relaxations that are solved in the probing procedure can be used to update the pseudo-cost values. In addition, domain propagation and other variable fixing strategies can be used to reduce the amount of QP relaxations that need to be solved.

Other presolving techniques such as cut generation can be applied using the condensed problem, and can be fully transferred to the original OCP formulation. In the present paper, we refrain from using cut generation techniques as they produce inequalities that potentially couple variables across stages. Such coupling between stages is not desirable as we rely on a block-sparsity exploiting QP solver.

III-C Resulting MIQP Algorithm for Optimal Control

Algorithm 3 describes the most important steps in our proposed B&B method for solving the MIQP in (1). It solves a block-structured QP relaxation using PRESAS [24] at every node and utilizes reliability branching (Algorithm 1) to decide the branching variables. As discussed earlier, the node selection strategy is based on a depth-first search followed by a best-first search as soon as an integer-feasible solution has been found. Note that the upper bound value UB provided to Alg. 3 can be based on an integer-feasible solution guess or it can be set to $+\infty$. Because of space limitations, the present paper will not further discuss all parameter choices in the algorithm such as, e.g., the reliability branching parameters, presolving frequency, memory usage, etc.

IV MIXED-INTEGER MPC ALGORITHM

In embedded applications of mixed-integer MPC, one needs to solve an MIQP (1) at each sampling instant under strict timing constraints. We can leverage the fact that we solve a sequence of similar problems (parametrized by the initial condition $\hat{x}_{0}$), in order to warm-start the B&B-algorithm. We refer to our proposed warm-starting procedure as tree propagation, because the main goal is to “propagate” the B&B tree forward by one time step. We describe this process in detail below. Then, we present the resulting mixed-integer MPC algorithm.

IV-A Warm Starting based on Tree Propagation

The warm-starting procedure aims to use knowledge of one MIQP, i.e., the search tree after solving the problem, in order to improve the B&B search for the next MIQP. Our idea is to store the path from the root to the leaf node where the optimal solution to the MIQP was found, as well as the branching order of the variables. We can then perform a shifting of this path in order to obtain a “warm-started tree” to start our search to solve the MIQP at the next time step. We illustrate this procedure in Figure 2, where the optimal path at the current time step is denoted by the sequence of nodes $P_{1}\rightarrow P_{2}\rightarrow P_{4}\rightarrow P_{6}\rightarrow P_{7}$. Let us consider a corresponding sequence of variables $u_{2}\rightarrow u_{3}\rightarrow u_{0}\rightarrow u_{1}$ that we branched on in order to create such optimal path. After shifting by one time step, all branched variables in the first control interval can be ignored, e.g., resulting in a shifted and shorter path of variables $u_{1}\rightarrow u_{2}\rightarrow u_{0}$.

At the subsequent time step, after obtaining the new state estimate, we execute all presolving techniques and we solve the QP relaxation corresponding to the root node. After removing from the warm-started tree the nodes that correspond to branched variables which are already integer feasible in the relaxed solution at the root node, we proceed by solving all the leaf nodes on the warm-started path. As we solve both children of a node on this path, we do not have to solve the parent node itself and therefore reduce computations by solving less QP relaxations. Hence, we go over the tree in the order depicted by the index of each node in Fig. 2. After the warm-started branch has been explored, we resume normal procedure of the B&B method. Algorithm 4 summarizes the proposed tree propagation technique.

The sequential nature of the problem also allows to shift and re-use the pseudo-cost information from one MPC time step to the next. This idea has the potential of producing smaller search trees as the MPC progresses, without the need to perform strong branching at every MPC step. The propagation of pseudo-costs can be coupled with an update of the reliability parameters to improve the overall performance. For example, the reliability number should be reduced for each variable from one time step to the next, in order to force strong branching for variables that have not been branched on in a sufficiently long time. In addition, nodes can be removed from the warm-started path in case they correspond to branched variables for which there is no pseudo-cost information or it is not sufficiently reliable, in an attempt to avoid bad branching decisions. Finally, these warm-started pseudo-costs can also be used to re-order the warm-started tree, in order to result in smaller search tree sizes.

The proposed tree propagation technique, with the additional re-use of pseudo-cost information, has been summarized in Algorithm 4. This procedure can improve the overall performance of the B&B method in multiple ways. First of all, the optimal path and pseudo-cost information is re-used to make better branching decisions for the mixed-integer program at the next time step, because the search trees are often similar for two subsequent problems. Also, the computational cost can be reduced by solving less QP relaxations to explore the warm-started tree. In addition, the shifted optimal path can be used in an attempt to efficiently obtain an integer-feasible solution, and therefore an important upper bound in the B&B algorithm, for the MPC problem at the next time step. Lastly, one can store the relaxed QP solutions on the optimal path, shift them by one time step and use them to warm-start the QP solver for nodes on the shifted optimal path.

IV-B MI-MPC Algorithm Implementation

Algorithm 5 summarizes the proposed MI-MPC algorithm. It solves a sequence of MIQPs where the branch-and-bound tree is warm-started at every time step, as well as the pseudo-cost and QP condensing information. As mentioned earlier, the B&B strategy and the additional presolve, warm-start and heuristic branching techniques have been implemented in Matlab, based on a C code implementation of the PRESAS algorithm [24] to solve each QP relaxation. In Section V, we illustrate the computational performance of the presented MI-MPC algorithm, including these presolving and warm-starting techniques, for two numerical case studies of mixed-integer MPC. A self-contained C code implementation is part of ongoing work, in order to illustrate the computational efficiency of the proposed algorithmic techniques.

IV-C Real-Time Embedded Applications of MI-MPC

Note that, in practice, the proposed warm-starting strategies often allow one to obtain an integer-feasible solution in a computationally efficient manner. However, even if the tree propagation immediately provides the globally optimal solution to the MIQP (1), a branch-and-bound algorithm still needs to perform relatively many iterations to prove optimality by pruning remaining nodes in the search tree. This motivates the use of a maximum number of B&B iterations in order to meet strict timing requirements of the embedded control application. Even if the algorithm does not terminate within this specified number of iterations, a feasible or even optimal solution may be available. This and other heuristic strategies for real-time implementation of MI-MPC are straightforward but outside the scope of this paper.

V CASE STUDIES: MIXED-INTEGER MPC

We report two numerical case studies to illustrate the computational performance of our MIQP-based MPC algorithm: a hybrid MPC test example and a satellite orbit re-centering application with a no-go zone in the orbital path. Our branch-and-bound algorithm has been implemented in Matlab in conjunction with the PRESAS active-set solver in C. To evaluate the performance, we compare our algorithm with the state-of-the-art GUROBI [22] and MOSEK [23] solvers for mixed-integer programming. It is important to emphasize that all advanced presolve and heuristic options have been activated for both software tools, resulting in fair computational comparisons.

V-A Hybrid MPC: Benchmark Example

The first case study is a hybrid MPC problem from [2], with the default settings as in bm99sim.m, which is a part of the Hybrid Toolbox for Matlab. This demo example has been used also more recently for numerical comparisons in [14]. The system is modeled using the HYSDEL toolbox [29] to obtain the mixed logical dynamical (MLD) system formulation. Figure 3 illustrates the average and worst-case CPU times taken by our algorithm, GUROBI and MOSEK for a range of control horizon lengths $N$.

Table I presents a detailed comparison for this test example, including additional timing results for the MI-NNLS solver that are taken directly from [14]. The latter computational results can serve only as a reference since they have been obtained on a different computer, with respect to the one used here with a 2.80 GHz Intel Xeon E3-1505M v5 processor and 32 GB of RAM. An important feature of our method is that its worst-case computation time is often rather close to the average performance in closed-loop MI-MPC simulations. This highlights the effectiveness of our tree propagation warm-starting procedure, such that consecutive branch-and-bound trees have approximately the same size. In addition, it can be observed from Table I that our proposed BB-PRESAS solver is either competitive with, or is a factor $2$ or $3$ times faster than GUROBI. The computational speedup is much larger when compared with other state-of-the-art tools such as MOSEK, our solver can be more than $10$ times faster in this particular MI-MPC test example. It shall be noted that GUROBI is a heavily optimized and fairly large software, which is unlikely to be amenable for embedded microprocessors, due to its code size, memory requirements, and software library dependencies.

V-B Satellite Station Keeping with No-Go Zones

The second case study is motivated by a real-world application, namely, orbit control of a satellite in a circular low earth orbit, 400km from earth surface. The satellite propulsion system is composed of two on/off thrusters, one on each of the in-track faces of the satellite, with gimbals rotating along the vertical axis and subject to angle constraints [30]. Thus, the propulsion system is controlled by two binary and two continuous control signals. The satellite dynamics are formulated by relative motion equations (HCW) with respect to the target position along the orbit, and the cone constraints of the thrust forces are approximated as simplexes [30]. Here, we consider a re-centering maneuver in which the satellite, previously drifting, is re-centered close to the target position along the orbit. Furthermore, the error coordinates from the target position are constrained in a station keeping window ($-300\leq X\leq 300,-150\leq Y\leq 150$).

Thus, our problem is simplified from [30], by considering only the orbital dynamics in the orbital plane, i.e., ignoring the out-of-orbital-plane and attitude dynamics, and as a consequence using a simpler propulsion system with only two thrusters. To better highlight the potential of the MI-MPC method, we add an exclusion zone in the station keeping window, i.e., an area that must be avoided, which makes the allowed region of positions to be non-convex. This additional constraint is modeled using standard integer programming techniques (see, e.g., [17]), resulting in three additional binary variables for each prediction step of the mixed-integer OCP to implement the logical exclusion zone constraints. In Figure 4, we show the trajectory of the satellite in relative coordinates, where the origin is the desired satellite position along the orbit, for the simulation of the satellite controlled by the mixed-integer MPC. The depicted area in the figure corresponds to the station keeping window, in which the satellite should be kept, and the shaded area is the exclusion zone that must be avoided, at least pointwise in time. The computational timing results for this particular closed-loop MPC simulation can be found in Figure 5. One can observe that our proposed algorithm has a very competitive runtime at every MPC time step, when compared to the commercial GUROBI solver. Most importantly, the BB-PRESAS algorithm appears to perform at least as good for this particular case study in terms of worst-case computation times.

VI CONCLUSIONS & OUTLOOK

In this paper, we proposed a branch-and-bound algorithm for mixed-integer MPC that exploits the optimal control problem structure to strengthen variable bounds, re-use pseudo-costs and warm-start the search tree at every MPC time step. More specifically, tailored domain propagation and tree propagation strategies have been presented. We showed preliminary results that illustrate the computational performance of our algorithm for two different MI-MPC case studies. A compact, efficient, but self-contained C code implementation of the proposed algorithm is under development to enable real-time embedded applications of hybrid MPC.