# Evaluating Data-driven Performances of Mixed Integer Bilinear Formulations
for Book Placement Planning

Xuan Lin1, Gabriel I. Fernandez1, and Dennis W. Hong1 1X. Lin, Gabriel I.
Fernandez, and Dennis W. Hong are with the Robotics and Mechanisms Laboratory,
Department of Mechanical and Aerospace Engineering, University of California,
Los Angeles, CA 90095, USA<EMAIL_ADDRESS>
###### Abstract

Mixed integer bilinear programs (MIBLPs) offer tools to resolve robotics
motion planning problems with orthogonal rotation matrices or static moment
balance, but require long solving times. Recent work utilizing data-driven
methods has shown potential to overcome this issue allowing for applications
on larger scale problems. To solve mixed-integer bilinear programs online with
data-driven methods, several re-formulations exist including mathematical
programming with complementary constraints (MPCC), and mixed-integer
programming (MIP). In this work, we compare the data-driven performances of
various MIBLP reformulations using a book placement problem that has discrete
configuration switches and bilinear constraints. The success rate, cost, and
solving time are compared along with non-data-driven methods. Our results
demonstrate the advantage of using data-driven methods to accelerate the
solving speed of MIBLPs, and provide references for users to choose the
suitable re-formulation.

##  I Introduction

Optimization-based methods are useful tools for solving robotic motion
planning problems. Typical approaches such as mixed-integer convex programs
(MICPs) [1, 2], nonlinear or nonconvex programs (NLPs) [3, 4, 5] and mixed-
integer nonlinear programs (MINLPs) [6] offer powerful tools to formulate
these problems. However, each has its drawbacks. NLPs tend to converge to
local optimal solutions. In practice, local optimal solutions can have bad
properties, such as inconsistent behavior as they depend on initial guesses.
To introduce discrete variables, NLPs require complementary constraints [7],
which is numerically challenging as it violates the majority of Constraint
Qualifications established for standard nonlinear optimization [8]. Mixed-
integer programs (MIPs) explicitly treat the discrete variables, which is
usually solved with Branch and Bound [9], or Benders Decomposition [10, 11].
MIP solvers seek global optimal solutions with more consistent behavior than
NLP solvers, but can require impractically long solving times for problems
with a large number of integer variables [12]. MINLPs, such as mixed-integer
bilinear programs (MIBLPs), incorporate both integer variables and nonlinear
(bilinear) constraints, hence, very expressive. Unfortunately, we lack
efficient algorithms to tackle MINLPs. As a result, it is difficult to
implement most of the optimization schemes mentioned above online.

Figure 1: Two snapshots of a book placement trajectory.

Recently, researchers have started to investigate machine learning methods to
gather problem-specific heuristics and speed up the discrete or nonlinear
optimization solvers. Standard algorithms to solve MIPs such as branch-and-
bound rely on heuristics to quickly remove infeasible regions. Learning
methods can be used to acquire better heuristics. For example, [13] used graph
neural networks to learn heuristics. [14] used reinforcement learning to
discover efficient cutting planes. On the other hand, data can be collected to
learn and solve specific problems [15, 16]. In [16], the authors proposed CoCo
which turns prediction of binary variables into classification problems, where
each set of binary variables has a unique class label. For online solving, the
neural network proposes candidate binary solutions, reducing the problems into
convex programs. Similar strategy is used by [17] which used K-nearest-
neighbor to predict initial guesses for NLPs, followed by a searching scheme
to satisfy the nonlinear constraints.

For multi-agent motion planning problems, in addition to bilinear constraints
such as orthogonal rotation matrices, agents can operate under various modes
or configurations encoded with discrete variables. The discrete configuration
switches with bilinear constraints lead to challenging MIBLPs. There exist
several approaches to re-formulate the MIBLP into either an MICP or an NLP,
such that the aforementioned data-driven methods can directly apply to speed
up the solver, including:

  1. 1.

Turning binary variables into continuous variables with complementary
constraints, hence the problem becomes an MPCC [18], and data-driven methods
such as [17, 19] can directly apply.

  2. 2.

Turning bilinear constraints into mixed-integer linear constraints using
McCormick envelopes [20], hence the problem becomes MICP, and data-driven
methods such as [15, 16] can directly apply.

In this paper, we ask the following question:

Given a certain amount of problems solved offline as warm-start data, which
re-formulation can solve the original MIBLP faster and more reliably? How much
data is needed to achieve such performances?

To get a comprehensive answer of this question, we introduce the book
placement problem for a benchmark. Given a bookshelf with several books on
top, an additional book needs to be placed on the shelf with minimal
disturbance to the existing books. Fig. 1 shows two snapshots of the
trajectory to insert one book into the shelf while interacting with existing
books. The bookshelf problem works well as a good benchmark because it: 1) is
an MINLP that can be converted to an MICP problem with hundreds of integer
variables, 2) can easily be scaled to push algorithms to their limit, and 3)
has practical significance where data can reasonably be gathered, such as in
the logistics industry.

To summarize, our contributions are as follows:

  1. 1.

Formulate the book placement problem as an MIBLP and solve it with data-driven
methods.

  2. 2.

Benchmark the data-driven performance of different MIBLP re-formulations on
the bookshelf problem.

##  II Related Works

###  II-A Parametric Programming

Parametric programming is a technique to build a function that maps to the
optimization solutions from varying parameters of the optimization problem
[21, 22]. Previous works have investigated parametric programming on linear
programming [23], quadratic programming [24], mixed-integer nonlinear
programming [25]. As an implementation of parametric programming on controller
design, explicit MPC [24, 26] tried to solve the model-predictive control
problem offline and stored the active sets. When computed online, the problem
parameter can be used to retrieve the active sets. [24] solved a constrained
linear quadratic regulator problem with explicit MPC and proved that the
active sets are polyhedrons. Therefore, the parameter space can be partitioned
using an algorithm proposed by [27]. However, when the problem is non-convex,
the active sets do not in general form polyhedrons. For non-convex problems,
previous works [17, 28, 15] directly stored the data points and picked out
warm-start using non-parametric learning such as K-nearest neighbor (KNN), and
solved the online formulation. On the other hand, modern learning techniques
such as neural networks can learn an embedding of a larger set of parameters
that maps to the solutions [16]. One advantage of using neural network methods
is the capability to generalize to out-of-distribution situations that are not
included in the training set [29].

###  II-B Mixed-integer Programming

On the other hand, mathematical programs with complementary constraints (MPCC)
models discrete modes through continuous variables with complementary
constraints. Complementary constraints enforce a pair of variables such that
if one of them is non-zero, the other one should be zero. This constraint is
traditionally hard to solve and is sensitive to initial guesses. Algorithms
such as time-stepping [30], pivoting [31], central path methods [32] are
proposed to resolve complementarity. In the robotics community, complementary
constraints are typically used to optimize over gaits for trajectory
optimization [7, 33] or control with implicit contacts [34] where contact
forces and distance to the ground are complementary with each other. On the
other hand, complementary constraints can also be used to model binary
variables [18].

###  II-C Mathematical Programming with Complementary Constraints

Mathematical programs with complementary constraints (MPCC) model discrete
variables by turning them into continuous variables with an additional
complementary constraint. Complementary constraints enforce a pair of
variables such that if one of them is non-zero, the other one should be zero.
This constraint is challenging to solve as it violates the majority of
Constraint Qualifications established for standard nonlinear optimization [8].
Algorithms such as time-stepping [30], pivoting [31], and central path methods
[32] were proposed to resolve this challenge. In the robotics community,
complementary constraints have been used to optimize over gaits for trajectory
optimization [7, 33] or control with implicit contacts [34] where ground
contact forces and distance to the ground are complementary with each other.
On the other hand, complementary constraints were also used to model binary
variables [18].

###  II-D Item Arrangement Planning

Optimal object arrangement and re-arrangement problems have been studied by a
large amount of previous literature. Using the task and motion planning
framework, previous works [35, 36, 37] planned a sequence of motion primitives
to arrange the objects into the target configurations. [38] used MIP models to
optimize the item placements on a shelf to minimize the access time. [39] also
used MIP models to generate book relocation plans in a library. [40] proposed
an algorithm to determine where to relocate objects during the rearrangement
process. Those works focused on item arrangements while ignoring their
interactions during the manipulation process. Along with the stow project
published by Amazon Science [41], several recent works investigated how items
interact with each other within the limited stowing space, and how to generate
a motion plan considering object configurations and contact forces [19, 42].
However, this task remains challenging due to complex item interactions such
as friction, and the variety of item configurations.

##  III Book Placement Problem Setup

Figure 2: Complete formulation of the book placement problem.

Assume a 2D bookshelf with limited width $W\in\mathbb{R}$ and height
$H\in\mathbb{R}$ contains rectangular books where book $i$ has width
$W_{i}\in\mathbb{R}$ and height $H_{i}\in\mathbb{R}$ for $i=1,...,N-1$. A new
book, $i=N$, is to be inserted into the shelf. The bookshelf contains enough
books in various orientations, where in order to insert book $N$, the other
$N-1$ books may need to be moved, i.e., optimize for minimal movement of $N-1$
books. This paper focuses on inserting one book using a single robot motion
during which the book is always grasped. The problem settings can be extended
to a sequence of motions with re-grasping.

Fig. 2 shows the constraints, variables and objective function for the
bookshelf problem. The variables that characterize book $i$ are: position
$\boldsymbol{x}_{i}\in\mathbb{R}^{2}=[x_{i},y_{i}]$ and angle
$\theta_{i}\in\mathbb{R}$ about its centroid. $\theta_{i}=0$ when a book
stands upright. The rotation matrix is: $\textbf{R}_{i}\in\mathbb{R}^{2\times
2}=[cos(\theta_{i}),\ -sin(\theta_{i});\ sin(\theta_{i}),\ cos(\theta_{i})]$.
Let the 4 vertices of book $i$ be $\boldsymbol{v}_{i,k}\in\mathbb{R}^{2}$,
$k=1,2,3,4$. The constraint A in Fig. 2 shows the linear relationship between
$\boldsymbol{x}_{i}$ and $\boldsymbol{v}_{i,k}$, where
$\boldsymbol{h}_{i,k}\in\mathbb{R}^{2}$ is the constant offset vector from its
centroid to vertices. Constraint B enforces that all vertices of all books
stay within the bookshelf, a linear constraint. Constraint C enforces the
orthogonality of the rotation matrix, a bilinear (non-convex) constraint.
Constraint D enforces that the angle $\theta_{i}$ stays within
$[-90\degree,90\degree]$, storing books right side up.

To ensure that the final book positions and orientations do not overlap with
each other, separating plane constraints are enforced. For convex shapes, the
two shapes do not overlap with each other if and only if there exists a
separating hyperplane $\boldsymbol{a}^{T}\boldsymbol{x}=b$ in between for
$\boldsymbol{a}\in\mathbb{R}^{2}$ representing the plane normal direction and
$b\in\mathbb{R}$ representing plane offset [43]. That is, for any point
$\boldsymbol{p}_{1}\in\mathbb{R}^{2}$ inside shape 1 then
$\boldsymbol{a}^{T}\boldsymbol{p}_{1}\leq b$, and for any point
$\boldsymbol{p}_{2}\in\mathbb{R}^{2}$ inside shape 2 then
$\boldsymbol{a}^{T}\boldsymbol{p}_{2}\geq b$. This is represented by
constraint E. Constraint F enforces $\boldsymbol{a}$ to be a normal vector.
Both E and F are bilinear constraints.

Finally, we need to assign configurations to each book $i$. For each book, it
can be standing straight up, laying down on its left or right, or leaning
towards left or right against some other book, as shown in the far right
column in Fig. 2. For each book $i$, we assign $n_{z}=3$ binary variables
$\boldsymbol{z}_{i}\in\\{0,1\\}^{n_{z}}$. If book $i$ stands upright
($z_{i,2}=1$) or lays flat on its left ($z_{i,1}=1$) or right ($z_{i,3}=1$),
constraints I1 or J1 or H1 are enforced, respectively. If book $i$ leans
against another book (or left/right wall which can be treated as static books)
on the left or right, constraints in K and L are enforced, respectively. To
this end checks need to indicate the contact between books. By looking at the
right column in Fig. 2, we can reasonably assume that the separating plane
$\boldsymbol{a}^{T}\boldsymbol{x}=b$ always crosses vertex 1 of the book on
the left and vertex 4 of the book on the right. This is represented by
bilinear constraints, K1 and L1. In addition, the books need to remain stable
given gravity. Constraints K2 and L2 enforce that a book is stable if its $x$
position stays between the supporting point of itself (vertex 2 if leaning
rightward and vertex 3 if leaning leftward) and the $x$ position of the book
that it is leaning onto. Lastly, constraint K3 and L3 enforce that the books
have contact with the _ground_. For practical reasons, we assume that books
cannot stack onto each other, i.e, each book has to touch the _ground_ of
bookshelf at at least one point. We note that constraints in H, I, J, K, and L
can be formulated as mixed-integer linear constraints using big-M formulation
[44]. Any contact conditions between pairs of books may also be added into
this problem as long as it can be formulated as mixed-integer convex
constraint. Overall, this is a problem with integer variables,
$\boldsymbol{z}_{i}$, and bilinear constraints C, E, F, K1, and L1, hence, an
MIBLP problem.

Practically, this problem presents challenges for retrieving high-quality
solutions. If robots were used to store books, the permissible solving time
would be several seconds. Several non-data-driven approaches have been
benchmarked in this paper to solve this mixed-integer bilinear problem: fix
one set of bilinear variables and solve MICP [45], nonlinear ADMM approach
[46, 47], or directly apply off-the-shelf MINLP solvers such as BONMIN [48].
These non-data-driven approaches struggle to meet the requirements.

##  IV Data-driven Algorithm

Assume that we are given a set of problems parametrized by
$\boldsymbol{\Theta}\in\mathbb{R}^{p}$ that is drawn from a distribution
$D(\boldsymbol{\Theta})$. We seek to solve a parameterized version of MIBLP
[49]:

$\displaystyle\underset{\boldsymbol{x},\ \boldsymbol{z}}{\text{minimize}}\
||\boldsymbol{x}-\boldsymbol{x}_{g}(\boldsymbol{\Theta})||^{2}_{\boldsymbol{Q}}$
(1a) s.t.
$\displaystyle\boldsymbol{A}(\boldsymbol{\Theta})\boldsymbol{x}\leq\boldsymbol{h}(\boldsymbol{\Theta})$
(1b) $\displaystyle l_{i}\leq x[i]\leq u_{i}\ \ \ i=1,...,dim(x)$ (1c)
$\displaystyle x[j]\in\\{0,1\\}\ \ \ \ j\in\mathcal{J}$ (1d) $\displaystyle
x[r_{s}]=x[p_{s}]x[q_{s}]\ s=1,...,S$ (1e)

Where $\boldsymbol{x}\in\mathbb{R}^{dim(x)}$ denotes both continuous and
binary variables. $\boldsymbol{x}_{g}$ denotes the target values of
$\boldsymbol{x}$. $\mathcal{J}$ gives an index set for binary variables.
Constraints (1b) is mixed-integer linear, meaning if the binary variables are
relaxed into continuous variables $x[j]\in[0,1]\forall j\in\mathcal{J}$, (1b)
becomes linear. Equality constraints in (1b) are omitted as they can be turned
into two inequality constraints from opposite directions. Constraints (1e) are
bilinear.

As mentioned in section I, we compare the data-driven performance of two re-
formulations: converting a mixed-integer bilinear program to either an MIQP or
an MPCC, and solve it with data-driven methods. If the learning agent with the
problem parameter $\boldsymbol{\Theta}$ can provide a list of binary
variables, the problem reduces to a convex optimization permitting convex
solvers. On the other hand, if the learning agent provides a good initial
guess for the variables involved in complementary constraints, the feasibility
of solving the problem significantly increases. We first collect the dataset.
Let the dataset of size $D$ be a set of tuples
$\\{(\boldsymbol{\Theta}_{i},\boldsymbol{x}_{i})\\}$, $i=1,...,D$. We
construct $\\{\boldsymbol{\Theta}_{i}\\}$ from uniform samples, then solve
each $\boldsymbol{\Theta}_{i}$ with a non-data-driven approach to get
$\boldsymbol{x}_{i}$. The dataset is then used to learn candidate binary
variables or initial guesses for continuous variables. In the following
sections, we describe the implementation details for each online formulation.

###  IV-A MPCC Re-formulation

One re-formulation of (1) is to turn all the binary variables into continuous
variables with complementary constraints, such that the MIBLP formulation
becomes MPCC and can be solved through NLP solvers. We remove constraint (1d),
and impose the equivalent complementary constraint from [18]:

$x[j](1-x[j])\leq\epsilon\ \ j\in\mathcal{J}$ (2)

Where $\epsilon>0$ is a small constant to relax the constraint for better
solver performance. We have also tried other implementations in [50] and get
similar or worse performances.

With the dataset, the learner samples multiple initial guesses
$\boldsymbol{x}$ for the online problem $\boldsymbol{\Theta}$, using the
K-nearest-neighbor approach. The MPCC is then solved using the sampled
$\boldsymbol{x}$ as initial guesses, followed by a search scheme until
nonlinear constraints are satisfied. This algorithm is given by Algorithm 1.

Input: $\\{(\boldsymbol{\Theta}_{i},\boldsymbol{x}_{i})\\}$,
$\boldsymbol{\Theta}$

1 Find the K-nearest-neighbors $\boldsymbol{\Theta}_{1}$, …,
$\boldsymbol{\Theta}_{K}$ in $\\{\boldsymbol{\Theta}_{i}\\}$, sorted in order
of increasing Euclidean distance
$d(\boldsymbol{\Theta},\boldsymbol{\Theta}_{k})$. Retrieve
$\boldsymbol{x}_{1}$, …, $\boldsymbol{x}_{K}$.

2 for  _$k=1,...,K$_ do

3 Solve the MPCC re-formulation of (1) using applicable NLP solvers with
$\boldsymbol{x}_{k}$ as initial guess.

4 if  _successful_ then

5 return _$\boldsymbol{x}$_

6

return _nil_

Algorithm 1 MPCC Online Solving Procedure

###  IV-B MIQP Re-formulation

Another re-formulatoin is to convert bilinear constraints into mixed-integer
linear constraints. The key idea is for each constraint (1e), we partition the
range of $x[p_{s}]$ and $x[q_{s}]$ into small intervals, then approximate the
constraint within each interval using the linear McCormick envelope relaxation
[51]. For $x[r_{s}]=x[p_{s}]x[q_{s}]$, we set in total $n_{p_{s}}$ small
intervals for $x[p_{s}]$ as
$[l_{p_{s}},x^{1}_{p_{s}},,...,x^{n_{p_{s}}-1}_{p_{s}},u_{p_{s}}]$, and in
total $n_{q_{s}}$ small intervals for $x[q_{s}]$ as
$[l_{q_{s}},x^{1}_{q_{s}},,...,x^{n_{q_{s}}-1}_{q_{s}},u_{q_{s}}]$. We then
adopt the special ordered set of type 2 (sos2) formulation introduced by [52,
53] that used additional continuous variables
$\boldsymbol{\alpha}_{s}\in\mathbb{R}^{n_{p_{s}}+1}$,
$\boldsymbol{\beta}_{s}\in\mathbb{R}^{n_{q_{s}}+1}$,
$\boldsymbol{\gamma}_{s}\in\mathbb{R}^{(n_{p_{s}}+1)\times(n_{q_{s}}+1)}$, and
$n_{\delta}$ additional binary variables
$\boldsymbol{\delta}\in\\{0,1\\}^{n_{\delta}}$, $n_{\delta}=\lceil
log_{2}((n_{p_{s}}+1)\times(n_{q_{s}}+1))\rceil$, to formulate
$x[r_{s}]=x[p_{s}]x[q_{s}]$ into mixed-integer linear constraints:

$\displaystyle x[p_{s}]=\sum_{i}{\alpha_{s}[i]x^{i}_{p_{s}}},\ \
x[q_{s}]=\sum_{i}{\beta_{s}[i]x^{i}_{q_{s}}}$ (3) $\displaystyle
x[r_{s}]=\sum_{i}\sum_{j}\gamma_{s}[i,j]x^{i}_{p_{s}}x^{j}_{q_{s}}$
$\displaystyle\sum_{j}\gamma_{s}[i,j]=\alpha_{s}[i],\ \
\sum_{i}\gamma_{s}[i,j]=\beta_{s}[j],\ \ \gamma_{s}[i,j]\geq 0$
$\displaystyle\boldsymbol{\alpha}_{s},\boldsymbol{\beta}_{s}\ \ \text{are in
sos2 (enforced by $n_{\delta}$ binary variables)}$

We refer readers to [52] for more details of sos2 enforcement using binary
variables. Despite the relatively smaller number of binary variables due to
$log_{2}(\cdot)$ formulation, approximating nonlinear constraints with mixed-
integer linear constraints still generates a large number of integer variables
when high approximation accuracy is desired. As a result, the formulation
suffers from extended solving time without good warm-starts. We used offline
solved problems to train the learner. Online, the learner provides a) complete
set of binary variables to fix the binary variables in (1d); b) complete set
of continuous variables, such that one can locate the intervals for
$x[p_{s}]$, $x[q_{s}]$, and replace constraints in (1e) by linear constraints.
Therefore, the problem becomes QP permitting fast solvers. For the learner, we
have implemented both K-nearest neighbor approach [15], and the neural-network
approach [16]. This KNN approach is detailed Algorithm 2. For the neural-
network approach, we refer readers to [16] for details of network structure,
training, and online sampling. After the neural network samples
$\boldsymbol{x}_{1},...,\boldsymbol{x}_{K}$, the method to turn MIQPs into QPs
is identical to Algorithm 2.

Input: $\\{(\boldsymbol{\Theta}_{i},\boldsymbol{x}_{i})\\}$,
$\boldsymbol{\Theta}$

1 Find the K-nearest-neighbors $\boldsymbol{\Theta}_{1}$, …,
$\boldsymbol{\Theta}_{K}$ in $\\{\boldsymbol{\Theta}_{i}\\}$, sorted in order
of increasing Euclidean distance
$d(\boldsymbol{\Theta},\boldsymbol{\Theta}_{k})$. Retrieve
$\boldsymbol{x}_{1}$, …, $\boldsymbol{x}_{K}$.

2 for  _$k=1,...,K$_ do

3 for  _$s=1,...,S$_ do

4 Find the interval $[x^{i}_{p_{s}},x^{i+1}_{p_{s}}]$ and
$[x^{j}_{q_{s}},x^{j+1}_{q_{s}}]$ that $x_{k}[p_{s}]$, $x_{k}[q_{s}]$ stays
within, replace (1e) with linear constraints:

5 $\displaystyle
x[p_{s}]=\alpha_{s}[i]x^{i}_{p_{s}}+\alpha_{s}[i+1]x^{i+1}_{p_{s}}$
$\displaystyle
x[q_{s}]=\beta_{s}[j]x^{j}_{q_{s}}+\beta_{s}[j+1]x^{j+1}_{q_{s}}$
$\displaystyle
x[r_{s}]=\gamma_{s}[i,j]x^{i}_{p_{s}}x^{j}_{q_{s}}+\gamma_{s}[i+1,j]x^{i+1}_{p_{s}}x^{j}_{q_{s}}$
$\displaystyle+\gamma_{s}[i,j+1]x^{i}_{p_{s}}x^{j+1}_{q_{s}}+\gamma_{s}[i+1,j+1]x^{i+1}_{p_{s}}x^{j+1}_{q_{s}}$
$\displaystyle\gamma_{s}[i,j]+\gamma_{s}[i,j+1]=\alpha_{s}[i]$
$\displaystyle\gamma_{s}[i,j]+\gamma_{s}[i+1,j]=\beta_{s}[j]$
$\displaystyle\alpha_{s}[i]+\alpha_{s}[i+1]=1$
$\displaystyle\beta_{s}[j]+\beta_{s}[j+1]=1$
$\displaystyle\alpha_{s}[i],\alpha_{s}[i+1],\beta_{s}[j],\beta_{s}[j+1]\geq 0$
$\displaystyle\gamma_{s}[i,j],\gamma_{s}[i+1,j],\gamma_{s}[i,j+1],\gamma_{s}[i+1,j+1]\geq
0$

6 Fix $x[j]=x_{k}[j]\ \forall\ j\in\mathcal{J}$

7 Solve the resultant QP using applicable solvers.

8 if  _successful_ then

9 return _$\boldsymbol{x}$_

10

return _nil_

Algorithm 2 MIQP Online Solving Procedure

###  IV-C ADMM Algorithm for MINLPs

For benchmark, we implemented a non-data-driven ADMM algorithm designed for
MINLPs [46]. The idea is to separate the MIBLP into a MIQP formulation and an
NLP formulation with the exact copies of variables but different constraints,
and iterate between them until the consensus is reached. The MIQP formulation
contains constraints
$\textbf{A},\textbf{B},\textbf{D},\textbf{G},\textbf{H1},\textbf{I1},\textbf{J1},\textbf{K2},\textbf{K3},\textbf{L2},\textbf{L3}$,
while missing bilinear constraints
$\textbf{C},\textbf{E},\textbf{F},\textbf{K1},\textbf{L1}$. The NLP
formulation contains constraints
$\textbf{A},\textbf{B},\textbf{C},\textbf{D},\textbf{E},\textbf{F},\textbf{H},\textbf{I},\textbf{J},\textbf{K},\textbf{L}$,
while missing the mode constraint G saying $\boldsymbol{z}$s are integers. The
formulation is elaborated in Eqn. 4.

$\displaystyle\underset{\boldsymbol{x}_{i},\ \boldsymbol{R}_{i},\
\boldsymbol{v}_{i,k},\ \boldsymbol{a}_{j},\ b_{j},\
z_{i}}{\text{minimize}}\quad f_{obj}$ (4) subject to Mixed integer linear
constraints: A, B, D, G, H1, I1, J1, K2, K3, L2, L3 Nonlinear constraints: A,
B, C, D, E, F, H, I, J, K, L

The ADMM iterates the following 3 steps: first it solves the MIQP formulation
in (4), then solves the NLP formulation in (4). Finally it updates the
residual variables. The process continues till convergence. Let
$\boldsymbol{x}_{1}\in\mathbb{R}^{n_{x}}$ and
$\boldsymbol{x}_{2}\in\mathbb{R}^{n_{x}}$ be the copies of the variables in
formulation (4). Let $\boldsymbol{w}\in\mathbb{R}^{n_{x}}$ be the residual
variables to achieve consensus. $\gamma$ be a constant. $\boldsymbol{G}$ are
the weights. The first step is:

$\displaystyle\underset{\boldsymbol{x}_{1}}{\text{minimize}}\quad||\boldsymbol{x}_{1}^{i}-\boldsymbol{x}_{2}^{i}+\boldsymbol{w}^{i}||_{\boldsymbol{G}_{k}}$
s.t. Mixed-integer constraints in (4)

The next step solves the NLP:

$\displaystyle\underset{\boldsymbol{x}_{2}}{\text{minimize}}\quad||\boldsymbol{x}_{2}^{i}-(\boldsymbol{x}_{1}^{i+1}+\boldsymbol{w}^{i}))||_{\boldsymbol{G}_{k}}$
s.t. Nonlinear constraints in (4)

Finally, we update the residual variables and weights:

$\displaystyle\boldsymbol{w}^{i+1}=\boldsymbol{w}^{i}+\boldsymbol{x}_{1}^{i}-\boldsymbol{x}_{2}^{i}$
(5) $\displaystyle\boldsymbol{G}_{k+1}=\gamma\boldsymbol{G}_{k}$
$\displaystyle\boldsymbol{w}^{i+1}=\boldsymbol{w}^{i}/\gamma$

##  V Experiment

###  V-A Experiment Setup

We place 3 books inside the shelf where 1 additional book is to be inserted.
Grids are assigned to the variables involved in the bilinear constraint C, E,
F, K1, L1: $\boldsymbol{R}_{i}(\theta_{i})$, $\boldsymbol{a}_{j}$ and
$\boldsymbol{v}_{i,k}$. These variables span a 48-dim space.
$\boldsymbol{\Theta}\in\mathbb{R}^{17}$ includes the center positions, angles,
heights and widths of stored books and height and width of the book to be
inserted.

For MIQP re-formulation, Table I lists the number of grids for each variable
involved in the bilinear constraints. The rotation angles $\theta_{i}$, which
includes $\boldsymbol{R}_{i}$, are gridded at $\frac{\pi}{8}$ intervals.
Elements in $\boldsymbol{a}_{j}$ are gridded on 0.25 intervals. Elements in
$\boldsymbol{v}_{i,k}$ are gridded at intervals $\frac{1}{4}$ the shelf width
$W$ and height $H$. Our MIQP re-formulation results in 130 integer variables
in total, out of which 18 variables represent discrete configurations, and the
rest comes from bilinear constraints. For MPCC re-formulations, a small
relaxation of $\epsilon=10^{-8}$ on the right hand side of the complementary
constraints is used to increase the chance of getting a feasible solution.

TABLE I: Segmentations of Non-convex Variables variable | range | # of segmentations
---|---|---
Item orientation $\theta$ (rad) | [-pi/2, pi/2] | 8
Separating plane normal $a$ | [-1, 1] | 8
Vertex x position $v_{x}$ (cm) | [-9, 9] | 4
Vertex y position $v_{y}$ (cm) | [0, 11] | 4

The problem instances are generated using a 2-dim simulated environment of
books on a shelf. Initially, 4 randomly sized books are arbitrarily placed on
the shelf, and then 1 is randomly removed and regarded as the book to be
inserted. Contrary to the sequence, the initial state with 4 books represents
one feasible (not necessarily optimal) solution to the problem of placing a
book on a shelf with 3 existing books. This guarantees that all problem
instances are feasible. Since this problem can be viewed as high-level
planning for robotic systems, the simulated data is sufficient. For
applications outside of the scope of this paper real-world data may be
preferable in this pipeline.

All methods are tested on 400 randomly sampled bookshelf problems using the
same method mentioned above. The solving process is done on a Core i7-7800X
3.5GHz × 12 machine. For NLPs, we used off-the-shelf solver KNITRO with its
default interior point method. For MIQPs, we used off-the-shelf solver Gurobi.
For QPs, we used off-the-shelf solver OSQP that runs ADMM. We have also tried
other solvers such as IPOPT for NLPs, and Gurobi for QPs and got similar or
worse results. All results are listed in table V-A. We omit the error bars as
they are consistent among the results: success rates are within $\pm 3\%$,
averaged solving times are within $\pm 10ms$, and averaged objective values
are within $\pm 400$.

TABLE II: Benchmark Results

Formulation & |  |  |  |  |  |  |  |  |
---|---|---|---|---|---|---|---|---|---
initial guess | Zero initial guess |  |  |  |  |  |  |  |
MPCC | Manual initial guess |  |  |  |  |  |  |  |
MPCC | ADMM |  Data-driven |  |  |  |  |  |  |
initial guess MPCC |  |  | Data-driven |  |  |  |  |  |
initial guess MIQP |  |  |  |  |  |  |  |  |
Note | Zero |  |  |  |  |  |  |  |
initial guess | One Manual |  |  |  |  |  |  |  |
initial guess | ADMM for |  |  |  |  |  |  |  |
MINLPs | KNN selects |  |  |  |  |  |  |  |
top 3 guesses |  |  | CoCo select |  |  |  |  |  |
top 10 guesses |  |  |  |  |  |  |  |  |
Amount of data | 0 |  0 | 0 | 100 | 500 |  1000 | 15000 | 30000 | 80000
Success Rate | <1% |  78.25% | 96.5% | 93.25% | 97.75% |  99.40% | 94.50% | 95.50% | 98.75%
Avg. Solve Time | / |  29ms | 260ms | 61ms | 30ms |  28ms | 65ms | 45ms | 16ms
Max. Solve Time | / |  198ms | 1.29sec | 257ms | 128ms |  110ms | 190ms | 182ms | 174ms
Avg. Objective | / |  -8620 | -7250 | -7982 | -8492 |  -8506 | -8637 | -8631 | -8607
Avg. # trials | / |  1 | 4.73 | 1.27 | 1.05 |  1.02 | 4.85 | 3.85 | 2.57
Solver | Knitro |  Knitro | Gurobi+Knitro | Knitro |  |  | OSQP |  |

###  V-B Baselines

  1. 1.

Zero Initial Guesses As a baseline, we solve the MPCC re-formulation with the
default initial guess from the solvers (all zeros). The results are shown in
the column named “Zero initial guess MPCC” of Table V-A.

  2. 2.

Manually Designed Initial Guesses Another baseline is to solve the MPCC re-
formulation with manually designed initial guesses. The book placement
requires inserting one book minimizing the movements of the existing books.
Therefore, a good initial guess is to directly use the continuous and discrete
variables from the original scene. We solve the separating planes using the
original scene as initial guesses for $\boldsymbol{a}$ and $b$. The results
are shown in column named “Manual initial guess MPCC” of Table V-A.

  3. 3.

ADMM for MINLPs For the last non-data-driven baseline, we implemented ADMM
method designed for MINLP described in Section Section IV-C. Manually designed
initial guesses described above are used for the nonlinear sub-problem to
increase the feasibility rate. Efforts are devoted to tuning the weights
$\boldsymbol{G}$, $\gamma$, and scaling the variables properly to ensure the
performance is on the better end. The results are shown in the column named
“ADMM” of Table V-A.

Despite the $\epsilon$ relaxation of complementary constraints, the rate of
solving the problem successfully with zero warm-start is almost 0, revealing
the challenge of solving MPCCs without good initial guesses. The performance
is improved significantly with one manually designed initial guess, yet
remains unsatisfactory. We also noticed that the success rates for manually
designed initial guesses drop as the number of books increases. When there are
four books on the shelf, inserting the 5th book has less than $60\%$ success
rate. Therefore, this approach may not scale to problems with much larger
sizes.

ADMM works fairly well for getting feasible solutions. However, it usually
takes 4-5 iterations (1 iteration includes 1 MIQP and 1 NLP) to get a feasible
solution, hence the solving speed is impaired. The nonlinear formulation takes
more than 70% of the solving time, due to the separating plane constraints. In
addition, the average optimal cost from ADMM is harmed by the extra consensus
terms in the objective function to guide the convergence.

###  V-C MPCC Re-formulation

We used the K-nearest neighbor algorithm to pick out the top three candidate
solutions from the dataset of sizes 100, 500, and 1000, respectively. The
results are shown in the column named “Data-driven initial guess MPCC” of
Table V-A. A simple KNN algorithm with a few candidate solutions can
significantly increase the rate of getting a feasible solution to more than
$98\%$. The average trial of solving the problem is only slightly more than 1
showing that most problem instances are solved with only one trial. The
average solving speed is also one order of magnitude faster than the non-data-
driven ADMM algorithm.

###  V-D MIQP Re-formulation

We trained the CoCo network to pick out the first 10 candidate binary
solutions from datasets of sizes 15000, 30000, and 80000, respectively. The
resultant QPs were then solved by OSQP, which used the previous solutions to
warm-start the next problem for faster solving speeds. The results are shown
in the column named “Data-driven initial guess MIQP” of Table V-A. Our results
show that for this specific problem, the MIQP re-formulation achieves the
fastest solving speed among the formulations tested. On the other hand, CoCo
takes on more trials to find a feasible solution than the MPCC re-formulation.
We also tested the KNN approach to select the first 10 candidate binary
solutions and got similar results as using the CoCo network.

###  V-E Analysis of Results

According to our test, data-driven methods can simultaneously achieve better
optimality, faster solving speed, and higher success rates than non-data-
driven methods for the book placement problem. The MPCC re-formulation and
MIQP re-formulation differs mostly at the amount of data used to achieve
similar solving speeds and success rates. The solving speeds heavily rely on
the solver used. Since we do not have access to the details of numerical
schemes inside KNITRO or Gurobi, we avoid reasoning in this aspect.

Our MIQP re-formulation achieves $>95\%$ success rate using a dataset about
size 30000. For datasets less than 10000, CoCo has success rates generally
$<90\%$. On the other hand, MPCC re-formulation achieves $>90\%$ success rates
with the number of data on the magnitude of hundreds. We noticed that KNITRO
rarely changes the initial binary guesses for $x[j],j\in\mathcal{J}$.
Therefore, this difference comes from the bilinear constraints (1e) converted
into multiple linear envelopes with additional binary variables
$\boldsymbol{\delta}$. For MIQP re-formulation, the learner directly provide
$\boldsymbol{\delta}$ that fixes the intervals
$[x^{i}_{p_{s}},x^{i+1}_{p_{s}}]$ and $[x^{j}_{q_{s}},x^{j+1}_{q_{s}}]$ for
$x[p_{s}]$ and $x[q_{s}]$ by adding linear constraints. The QP solver finds
$x[p_{s}]$ and $x[q_{s}]$ inside the intervals. Since $x[p_{s}]$ and
$x[q_{s}]$ within the predicted intervals, if the learner provides the wrong
intervals, the QP is infeasible, and the learner would need to try again.
Sufficiently large dataset is required to train the learner to predict the
correct intervals. However, for MPCC re-formulation, KNN provides initial
guesses for $x[p_{s}]$, $x[q_{s}]$, then the NLP solver uses numerical schemes
to satisfy (1e) and other constraints. Since NLP solvers are designed to
search and satisfy nonlinear constraints, even restore from infeasible initial
guesses, the learner does not have to guarantee the feasibility of initial
guesses as seen by smaller numbers of trials.

Roughly speaking, if one has a limited amount of data, our results support
keeping the bilinear constraints (1e) and letting the learner provide initial
guesses of $x[p_{s}]$, $x[q_{s}]$ to the NLP solvers, instead of decomposing
(1e) into multiple linear envelopes like (3) does and letting the learner
select the correct intervals for $x[p_{s}]$, $x[q_{s}]$ to reside in.

##  VI Conclusion, Discussion and Future Work

This paper compares data-driven performances of MPCC or MIP re-formulations
applied to the book placement problem formulated as a MIBLP.

Data-driven methods can simultaneously achieve better optimality, faster
solving speed, and higher success rates than non-data-driven methods. The MPCC
re-formulation has a higher chance of finding feasible solutions than the MIQP
re-formulation from less amount of data. The MIQP re-formulation requires
significantly more data to achieve similar performances as the MPCC re-
formulation.

An obvious challenge with the data-driven approach is how to generalize the
training to cases outside the dataset. To achieve this goal, simple learners
such as KNN do not perform well, and learning methods such as auto-encoders
are to be explored.

## References

  * [1] R. Deits and R. Tedrake, “Footstep planning on uneven terrain with mixed-integer convex optimization,” in _2014 IEEE-RAS international conference on humanoid robots_. IEEE, 2014, pp. 279–286.
  * [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] H. Dai, A. Valenzuela, and R. Tedrake, “Whole-body motion planning with centroidal dynamics and full kinematics,” in _IEEE RAS Int. Conf. on Humanoid Robots_ , 2014.
  * [4] A. W. Winkler, C. D. Bellicoso, M. Hutter, and J. Buchli, “Gait and trajectory optimization for legged systems through phase-based end-effector parameterization,” _IEEE Robotics and Automation Letters_ , vol. 3, no. 3, pp. 1560–1567, 2018.
  * [5] Y. Shirai, X. Lin, Y. Tanaka, A. Mehta, and D. Hong, “Risk-aware motion planning for a limbed robot with stochastic gripping forces using nonlinear programming,” _IEEE Robotics and Automation Letters_ , vol. 5, no. 4, pp. 4994–5001, 2020.
  * [6] M. Soler, A. Olivares, P. Bonami, and E. Staffetti, “En-route optimal flight planning constrained to pass through waypoints using minlp,” 2011.
  * [7] M. Posa, C. Cantu, and R. Tedrake, “A direct method for trajectory optimization of rigid bodies through contact,” _Int. Journal of Robotics Research_ , vol. 33, no. 1, pp. 69–81, 2014.
  * [8] Z.-Q. Luo, J.-S. Pang, and D. Ralph, _Mathematical programs with equilibrium constraints_. Cambridge University Press, 1996.
  * [9] S. Boyd and J. Mattingley, “Branch and bound methods,” _Notes for EE364b, Stanford University_ , pp. 2006–07, 2007.
  * [10] A. M. Geoffrion, “Generalized benders decomposition,” _Journal of optimization theory and applications_ , vol. 10, pp. 237–260, 1972.
  * [11] X. Lin, “Generalized benders decomposition with continual learning for hybrid model predictive control in dynamic environment,” _arXiv preprint arXiv:2310.03344_ , 2023.
  * [12] X. Lin, M. S. Ahn, and D. Hong, “Designing multi-stage coupled convex programming with data-driven mccormick envelope relaxations for motion planning,” in _2021 IEEE/RSJ International Conference on Robotics and Automation (ICRA)_. IEEE, 2021.
  * [13] V. Nair, S. Bartunov, F. Gimeno, I. von Glehn, P. Lichocki, I. Lobov, B. O’Donoghue, N. Sonnerat, C. Tjandraatmadja, P. Wang _et al._ , “Solving mixed integer programs using neural networks,” _arXiv preprint arXiv:2012.13349_ , 2020.
  * [14] Y. Tang, S. Agrawal, and Y. Faenza, “Reinforcement learning for integer programming: Learning to cut,” in _International Conference on Machine Learning_. PMLR, 2020, pp. 9367–9376.
  * [15] J.-J. Zhu and G. Martius, “Fast non-parametric learning to accelerate mixed-integer programming for hybrid model predictive control,” _IFAC-Papers Online_ , vol. 53, no. 2, pp. 5239 – 5245, 2020.
  * [16] 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, 2022.
  * [17] K. Hauser, “Learning the problem-optimum map: Analysis and application to global optimization in robotics,” _IEEE Transactions on Robotics_ , vol. 33, no. 1, pp. 141–152, 2016.
  * [18] J. Park and S. Boyd, “A semidefinite programming method for integer convex quadratic minimization,” _Optimization Letters_ , vol. 12, no. 3, pp. 499–518, 2018.
  * [19] 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.
  * [20] A. Gupte, S. Ahmed, M. S. Cheon, and S. Dey, “Solving mixed integer bilinear problems using milp formulations,” _SIAM Journal on Optimization_ , vol. 23, no. 2, pp. 721–744, 2013.
  * [21] A. V. Fiacco, _Introduction to sensitivity and stability analysis in nonlinear programming_. Academic press, 1983, vol. 165.
  * [22] E. N. Pistikopoulos, V. Dua, N. A. Bozinis, A. Bemporad, and M. Morari, “On-line optimization via off-line parametric optimization tools,” _Computers & Chemical Engineering_, vol. 26, no. 2, pp. 175–185, 2002.
  * [23] T. Gal, _Postoptimal Analyses, Parametric Programming, and Related Topics: degeneracy, multicriteria decision making, redundancy_. Walter de Gruyter, 2010.
  * [24] 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.
  * [25] V. Dua and E. N. Pistikopoulos, “Algorithms for the solution of multiparametric mixed-integer nonlinear optimization problems,” _Industrial & engineering chemistry research_, vol. 38, no. 10, pp. 3976–3987, 1999.
  * [26] P. Tøndel, T. A. Johansen, and A. Bemporad, “An algorithm for multi-parametric quadratic programming and explicit mpc solutions,” _Automatica_ , vol. 39, no. 3, pp. 489–497, 2003.
  * [27] V. Dua and E. N. Pistikopoulos, “An algorithm for the solution of multiparametric mixed integer linear programming problems,” _Annals of operations research_ , vol. 99, no. 1, pp. 123–139, 2000.
  * [28] G. Tang and K. Hauser, “A data-driven indirect method for nonlinear optimal control,” _Astrodynamics_ , vol. 3, no. 4, pp. 345–359, 2019.
  * [29] T. Power and D. Berenson, “Variational inference mpc using normalizing flows and out-of-distribution projection,” _arXiv preprint arXiv:2205.04667_ , 2022.
  * [30] M. Anitescu and F. A. Potra, “Formulating dynamic multi-rigid-body contact problems with friction as solvable linear complementarity problems,” _Nonlinear Dynamics_ , vol. 14, no. 3, pp. 231–247, 1997.
  * [31] E. Drumwright, “Rapidly computable viscous friction and no-slip rigid contact models,” _arXiv preprint arXiv:1504.00719_ , 2015.
  * [32] M. Kojima, N. Megiddo, T. Noma, and A. Yoshise, “A unified approach to interior point algorithms for linear complementarity problems: A summary,” _Operations Research Letters_ , vol. 10, no. 5, pp. 247–254, 1991.
  * [33] 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.
  * [34] S. L. Cleac’h, T. Howell, M. Schwager, and Z. Manchester, “Fast contact-implicit model-predictive control,” _arXiv preprint arXiv:2107.05616_ , 2021.
  * [35] A. Krontiris and K. E. Bekris, “Dealing with difficult instances of object rearrangement.” in _Robotics: Science and Systems_ , vol. 1123, 2015.
  * [36] C. Nam, J. Lee, S. H. Cheong, B. Y. Cho, and C. Kim, “Fast and resilient manipulation planning for target retrieval in clutter,” in _2020 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2020, pp. 3777–3783.
  * [37] R. Wang, K. Gao, D. Nakhimovich, J. Yu, and K. E. Bekris, “Uniform object rearrangement: From complete monotone primitives to efficient non-monotone informed search,” in _2021 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2021, pp. 6621–6627.
  * [38] L. Y. Chen, H. Huang, and K. Goldberg, “Optimal arrangement and rearrangement of objects on shelves to minimize robot retrieval cost,” _IEEE Transactions on Automation Science and Engineering_ , 2023.
  * [39] A. Ohnishi, S. Koichi, and A. Suzuki, “A mip model to generate an optimal book relocation plan for improving the management of shelf occupancy–a case study at nanzan university library–,” _Journal of Advanced Mechanical Design, Systems, and Manufacturing_ , vol. 16, no. 4, pp. JAMDSM0036–JAMDSM0036, 2022.
  * [40] S. H. Cheong, B. Y. Cho, J. Lee, C. Kim, and C. Nam, “Where to relocate?: Object rearrangement inside cluttered and confined environments for robotic manipulation,” in _2020 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2020, pp. 7791–7797.
  * [41] S. O’Neill. (2022) How Amazon Robotics researchers are solving a “beautiful problem”. Teaching robots to stow items presents a challenge so large it was previously considered impossible — until now. [Online]. Available: https://www.amazon.science/latest-news/how-amazon-robotics-researchers-are-solving-a-beautiful-problem
  * [42] H. Chen, Y. Niu, K. Hong, S. Liu, Y. Wang, Y. Li, and K. R. Driggs-Campbell, “Predicting object interactions with behavior primitives: An application in stowing tasks,” in _7th Annual Conference on Robot Learning_ , 2023.
  * [43] S. Boyd, S. P. Boyd, and L. Vandenberghe, _Convex optimization_. Cambridge university press, 2004.
  * [44] J. P. Vielma, “Mixed integer linear programming formulation techniques,” _Siam Review_ , vol. 57, no. 1, pp. 3–57, 2015.
  * [45] M. Posa, M. Tobenkin, and R. Tedrake, “Stability analysis and control of rigid-body systems with impacts and friction,” _IEEE Transactions on Automatic Control_ , vol. 61, no. 6, pp. 1423–1437, 2015.
  * [46] Y. Shirai, X. Lin, A. Schperberg, Y. Tanaka, H. Kato, V. Vichathorn, and D. Hong, “Simultaneous contact-rich grasping and locomotion via distributed optimization enabling free-climbing for multi-limbed robots,” in _2022 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. IEEE, 2022, pp. 13 563–13 570.
  * [47] X. Lin, G. I. Fernandez, Y. Liu, T. Zhu, Y. Shirai, and D. Hong, “Multi-modal multi-agent optimization for limms, a modular robotics approach to delivery automation,” in _2022 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. IEEE, 2022, pp. 12 674–12 681.
  * [48] P. Bonami and J. Lee, “Bonmin user’s manual,” _Numer Math_ , vol. 4, pp. 1–32, 2007.
  * [49] M. Fischetti and M. Monaci, “A branch-and-cut algorithm for mixed-integer bilinear programming,” _European Journal of Operational Research_ , vol. 282, no. 2, pp. 506–514, 2020.
  * [50] O. Stein, J. Oldenburg, and W. Marquardt, “Continuous reformulations of discrete–continuous optimization problems,” _Computers & chemical engineering_, vol. 28, no. 10, pp. 1951–1966, 2004.
  * [51] G. P. McCormick, “Computability of global solutions to factorable nonconvex programs: Part i—convex underestimating problems,” _Mathematical programming_ , vol. 10, no. 1, pp. 147–175, 1976.
  * [52] J. P. Vielma and G. L. Nemhauser, “Modeling disjunctive constraints with a logarithmic number of binary variables and constraints,” _Mathematical Programming_ , vol. 128, pp. 49–72, 2011.
  * [53] H. Dai, G. Izatt, and R. Tedrake, “Global inverse kinematics via mixed-integer convex optimization,” _Int. Journal of Robotics Research_ , vol. 38, no. 12-13, pp. 1420–1441, 2019.