(LABEL:)Eq.

11institutetext:  D. Yang  ✉22institutetext: University of California at Los
Angeles, Los Angeles, USA
22email<EMAIL_ADDRESS>33institutetext: P. Balaprakash
44institutetext: Argonne National Laboratory, Lemont, USA
44email<EMAIL_ADDRESS>55institutetext: S. Leyffer  66institutetext:
Argonne National Laboratory, Lemont, USA
66email<EMAIL_ADDRESS>
# Modeling Design and Control Problems Involving Neural Network Surrogates

Dominic Yang  Prasanna Balaprakash  Sven Leyffer

(Received: / Accepted: date)

###### Abstract

We consider nonlinear optimization problems that involve surrogate models
represented by neural networks. We demonstrate first how to directly embed
neural network evaluation into optimization models, highlight a difficulty
with this approach that can prevent convergence, and then characterize
stationarity of such models. We then present two alternative formulations of
these problems in the specific case of feedforward neural networks with ReLU
activation: as a mixed-integer optimization problem and as a mathematical
program with complementarity constraints. For the latter formulation we prove
that stationarity at a point for this problem corresponds to stationarity of
the embedded formulation. Each of these formulations may be solved with state-
of-the-art optimization methods, and we show how to obtain good initial
feasible solutions for these methods. We compare our formulations on three
practical applications arising in the design and control of combustion
engines, in the generation of adversarial attacks on classifier networks, and
in the determination of optimal flows in an oil well network.

###### Keywords:

Mixed-integer programming nonlinear programming complementarity constraints
machine learning neural networks

###### MSC:

MSC codes

###### Acknowledgements.

This material is based upon work supported by the U.S. Department of Energy,
Office of Science, Office of Advanced Scientific Computing Research, under
Contract DE-AC02-06CH11357. This work was also supported by the U.S.
Department of Energy through grant DE-FG02-05ER25694. The first author was
also supported through an NSF-MSGI fellowship.

##  1 Introduction and Background

We are interested in solving general optimization problems that include deep
neural networks (DNNs) that are used as surrogate models of complex functions
(e.g., physical processes aithal2019maltese , classification schemes
krizhevsky2012imagenet ; simonyan2014very ; he2016deep ). In particular, we
consider optimization problems of the form

$\mathop{\mbox{minimize}}_{x}\;f(\text{DNN}(x),x)\quad\mbox{subject
to}\;c(\text{DNN}(x),x)\leq 0,\;x\in X,$ (1)

where $f$ and $c$ are smooth functions representing the objective function and
constraints, respectively; $\text{DNN}(x)$ is the output of a DNN at $x$
(which we assume to have been previously trained on suitable data); $x$ are
the optimization variables; and $X\subset\mathbb{R}^{n}$ is a compact set that
may include integer restrictions. We consider feedforward neural networks in
this paper that are composed of a sequence of multiple layers of neurons. The
values of neurons in layer $\ell$, $x^{\ell}$ are a nonlinear function (the
activation function) applied to a linear transformation of the values in the
prior layer:

$x^{\ell}=\sigma(W^{{\ell}^{T}}x^{\ell-1}+b^{\ell}),\;\ell=1,\ldots,L.$ (2)

We have $x^{0}=x$, $x^{L}=\text{DNN}(x)$, $L$ is the number of layers, $W$ and
$b$ are weights determined by a training procedure, and $\sigma$ is the
activation function applied componentwise. Unless otherwise specified, we will
take $\sigma$ to be the ReLU function: $\text{ReLU(x)}=\max\\{x,0\\}$.

An example of an optimization problem we wish to solve is to minimize the
output of a neural network regressor that predicts the quantity of emissions
from automobile engine specifications. In this case $X$ would be the set of
existing automobile specifications, the constraints $c$ would be constraints
to ensure the engine is realistic, and the objective would be a function of
the emissions or engine performance. We will discuss this problem in detail in
Section 2.1.

To arrive at a tractable form of (1), we make the following assumptions.

###### Assumption 1

We assume that the following conditions are satisfied for problem (1):

  1. 1.

The objective function, $f:\mathbb{R}^{n}\to\mathbb{R}$, and the constraint
functions, $c:\mathbb{R}^{n}\to\mathbb{R}^{m}$, are twice continuously
differentiable convex functions.

  2. 2.

The set $X\subset\mathbb{R}^{n}$ is convex and compact.

The most restrictive assumption is the convexity assumption. We can relax this
assumption by leveraging standard global optimization techniques (see, e.g.,
sahinidis:96 ; tawarmalani.sahinidis:02 ; belotti2009couenne ), at the expense
of making the reformulated problem harder to solve. The smoothness assumption
on $f$ and $c$ can be relaxed to Lipschitz continuity by using subgradients,
and the compactness assumption is typically satisfied as long as $x$ are
constrained by bounds or the constraint function has compact level sets.
Below, we consider deep neural networks with ReLU activation functions that
satisfy this assumption.

##### Outline and Contributions.

We start by discussing three applications that employ DNN surrogates within an
optimization problem such as (1): (a) the minimization of emissions in an
engine design problem (as in aithal2019maltese ), (b) the generation of
optimal adversarial examples to “fool” a given classifier (as in
fischetti2018deep ), and (c) the optimization of the pump configuration for
oil (as in grimstad2019relu ). We develop a warmstart heuristic for both
approaches that generates good initial guesses from the training data and
helps us overcome the challenges of the nonconvex formulation. We also show
how to add constraints that restrict the optimization problem to search only
near where existing training points can be added.

Next, we demonstrate empirically that simply including a ReLU DNN within an
optimization problem can lead to poor convergence results. In particular, we
have developed a new nonlinear constraint for JuMP JuMP  that allows us to
directly include DNNs within an optimization model specified in JuMP. We then
give a compact characterization of stationarity of the embedded formulation of
our optimization problem.

Next, we consider alternative formulations and show that DNNs that use purely
ReLU activation functions can be formulated as mixed-integer sets, building on
fischetti2018deep ; anderson2019strong . We can then formulate (1)
equivalently as a convex mixed-integer nonlinear program (MINLP), which we
refer to as the mixed-integer program (MIP) formulation. We also introduce a
new formulation that results in a nonconvex nonlinear program (NLP) with
complementarity constraints. We refer to this as the mathematical program with
complementary constraints (MPCC) formulation. We prove this formulation has
stationarity conditions equivalent to the embedded formulation. Our
reformulations involve a lifting into a higher-dimensional space in which the
MINLP problem is convex.

We show empirically for each of our applications that for moderately sized
machine learning (ML) models the resulting programs can be solved by using
state-of-the-art MINLP and NLP solvers, making optimization problems with ML
models tractable in practice. We demonstrate that using the MIP formulation,
we can find provably optimal solutions for small problems. Using the embedded
and MPCC formulation, we show that we can address significantly larger
networks at the cost of guaranteeing only locally optimal solutions, and we
showcase scenarios where the MPCC formulation outperforms the embedded
formulation in terms of both optimal value found and consistency in
convergence to a solution.

Throughout this paper we assume that the deep neural network has been trained
and is fixed for the optimization, and we do not consider the question of
updating the neural network weights during the optimization loop. One
limitation of our approach is that we use standard MINLP solvers to tackle the
reformulated MINLP, which limits the size of the neural network,
$\text{DNN}(x)$, that can be used in the optimization.

##### Related Work.

Prior approaches to optimization over neural networks using MIP formulations
include (fischetti2018deep, ; cheng2017maximum, ; dutta2018output, ;
khalil2018combinatorial, ; serra2020empirical, ; tjeng2017evaluating, ). These
approaches primarily model the piecewise ReLU constraint using standard big-M
modeling tricks. They generally use the same basic formulation, but each
augments the solve by adding methods to tighten the big-M constraints
(fischetti2018deep, ; tjeng2017evaluating, ; grimstad2019relu, ), decomposing
the problem into smaller problems (khalil2018combinatorial, ), or adding local
search routines (dutta2018output, ). Anderson et al. (anderson2019strong, )
provide an in-depth overview of how to strengthen these models to an ideal
formulation with exponentially many constraints, as well as a method to
separate in linear time.

Some other approaches to these problems using methods from MINLP have been
tried. Cheon (cheon2020outer, ) solves an inverse problem over ReLU
constraints using an outer-approximation-inspired method, but without global
optimality guarantees. Katz et al. (katz2017reluplex, ) use an approach from
satisfiability modulo theory to address an optimization problem over neural
networks. Scheweidtmann and Mitsos (schweidtmann2019deterministic, ) use a
MINLP approach involving McCormick relaxations.

One major focus of these optimization problems is on testing the resilience of
neural networks against adversarial attack (carlini2016evaluating, ). This
involves either maximizing a notion of resilience (cheng2017maximum, ) or
finding minimal perturbations needed to misclassify an image
(fischetti2018deep, ). Some work has also been done in using optimization to
visualize features corresponding to neurons (fischetti2018deep, ) and for
surrogate optimization in the context of optimizing the production of a set of
oil wells (grimstad2019relu, ). Other applications involve the use of neural
networks as surrogates in the context of policy design for reinforcement
learning ryu2019caql ; delarue2020reinforcement . A recent paper by
Papalexopoulos et al. (papalexopoulos2021constrained, ) discusses the use of
ReLU neural network surrogates for the purposes of black-box optimization.

##  2 Modeling Optimization Applications Involving Neural Network Surrogates

In this section we describe three optimization models that make use of neural-
network surrogates, and we discuss some of the challenges that arise.

###  2.1 Optimal Design of Combustion Engine

Automobile engine operation is typically modeled by using highly intensive
physics-based simulation code, which even on powerful computers can take hours
just to model even a short drive. Hence, modeling the evaluation of this
simulation code by using a neural network surrogate model can produce
significant time savings, at the cost of producing slightly less accurate
results. A study examining this approach is documented in (aithal2019maltese,
).

Given a trained neural network, we could formulate several optimization
problems answering questions related to the operation of an engine on a given
commute. The problem we consider is that of minimizing emissions over the
course of a given drive. The resulting solution will then be both the engine
type and the drive style (e.g., RPM at all times in the drive) that produces
the most environmentally efficient commute. The following section demonstrates
how such a problem can be formulated.

####  2.1.1 Simple Engine Design and Control Problem

Suppose we have a trained DNN that predicts engine behavior based on engine
specifications and driving parameters, given in Table 1. The DNN is trained by
using 64 trips each split into 1,500 observations (25 minutes observed at
second intervals).

Input Parameter | Output Parameter
---|---
fuel injection (g) /s | nitrogen oxide, NO /s
engine RPM /s | carbon monoxide, CO /s
compression ratio | torque /s
Table 1: Input/output parameters of engine DNN model. Time-dependent
parameters are shown with time per second (/s).

The problem we propose to solve is the optimal design and control of an engine
for a given 25-minute trip. We will use a prescribed torque profile as a
surrogate for the trip characteristics and minimize a weighted sum of NO and
CO output. We have the following design variables that are input to the DNN:
fuel injection, $f$; engine RPM, $r$; and compression ratio, $c$. Note that
$f$ and $r$ are characteristics that change over time whereas $c$ is an engine
parameter that is fixed for the full drive. We also have variables NO and CO
that represent nitrogen oxide and carbon monoxide emissions produced by the
engine, as well as a variable torque that indicates engine torque. Each of
these quantities is predicted by the neural network for each time interval of
the drive. Formally, we state the following optimal design and control problem
over the time horizon $t=1,\ldots,T$ (where $T$ is the number of time
intervals):

$\begin{array}[]{lll}\displaystyle\mathop{\mbox{minimize}}_{f,r,c,\text{CO},\text{NO},\text{Torque}}&\displaystyle\sum_{t=1}^{T}(\text{NO}_{t}+\lambda\text{CO}_{t})\Delta
t&\text{min. NO and CO emissions}\\\ \mbox{subject
to}&\text{Torque}_{t}\geq\text{PrescribedTorque}_{t},\;\forall
t=1,\ldots,T&\text{trip profile}\\\
&\left(\text{NO}_{t},\text{CO}_{t},\text{Torque}_{t}\right)=\text{DNN}\left(f_{t},r_{t},c\right)&\text{DNN
constraints}\\\
&f_{t}\in[f_{\min},f_{\max}],\;r_{t}\in[r_{\min},r_{\max}],\;c\in[c_{\min},c_{\max}]&\text{bounds
on controls}.\end{array}$ (3)

We are using 1,500 time intervals (corresponding to $T=1500$ seconds in a
25-minute drive). As written, we will have $2T$ \+ 1 continuous control
variables in addition to the $3T$ output variables predicted by the DNN for a
total of $5T+1$ continuous variables in this formulation.

Our model has a separate DNN evaluation for each time interval $t$, meaning
each evaluation engenders different neural network activations. Because of the
presence of the design variable $c$, which is independent of $t$, this problem
does not decompose into individual time steps. One may construct a bilevel
optimization problem wherein on the upper level we decide $c$ and other
engine-level variables of interest and on the lower level we determine the
drive-specific variables that change over time by solving $T$ separate
optimization problems of much smaller size. This approach will not be
addressed in this paper.

####  2.1.2 Convex Hull Constraints

In addition to the constraints that encode the evaluation of a neural network,
constraints must be added to ensure that the solver does not extrapolate
significantly from the training data. Without these constraints, the
optimization routine may find that the optimal solution resides in an area for
which the neural network has not learned the behavior of the modeled function,
leading to a solution that, while optimal, does not reflect the true function
behavior and may be nonsensical. In fact, in an earlier implementation without
the convex hull constraint, we observed that the optimal design was obtained
for an engine that produced negative emissions. We rectify this error by
introducing constraints that constrain the input data to our model to be
within the convex hull of the input training data.

In our experiments we examined how the optimization models operated with
simple box constraints that bounded the input by the extremal values of the
training data, as well as with the convex hull constraints. Figure 1 shows the
location of fuel mass and RPM for solutions computed for a sample neural
network trained on our data set compared with the actual training data. The
corresponding solutions are displayed in Figure 2.

Figure 1: Location of training data and optimal solution for an engine design
problem on a neural network with 1 hidden layer of 16 nodes. On the left the
input is constrained by box constraints and on the right by convex hull
constraints.

Without the convex hull constraints, the optimal solution of neural network
evaluation appears to wander off to an area for which no training data exists.
This situation immediately creates problems in the computed solutions, which
are obviously nonsensical because they involve negative emissions. This
suggests that our neural network does not generalize well to data that has not
yet been seen by the network. On the other hand, constraining the input data
to be within the convex hull has the effect of producing solutions that look
like the training data.

Figure 2: Optimal configuration for each time step when using box constraints
(left) and using convex hull constraints (right). The red lines indicate the
drive cycles used to train the neural network, and the blue lines indicate the
optimal solution.

Unfortunately, these constraints have the effect of restricting possible
solutions to those that are close to the training data, when the function
being modeled may have a minimum that appears far from the training data. This
situation speaks to the limitations of using solely a surrogate model for
optimization. To find a global optimum of the original function would
necessarily involve an alternating approach wherein the surrogate model was
optimized and then this solution was queried against the modeled function for
new data to be added to the surrogate. This type of approach is addressed in
(queipo2005surrogate, ) and a similar approach is implemented in
papalexopoulos2021constrained .

If we fix the compression ratio, we can plot the contours of the objective
function as a function of fuel mass and RPM. Figure 3 shows the contours of
the objective function for a 5-layer DNN. We observe that this objective
function is highly nonconvex in the reduced space of the original variables
but it is convex in the lifted MINLP space according to Proposition 4. We also
observe that the convex hull constraint fulfills a second function by acting
as a mild convexifier of the problem by restricting the variables to a small
sliver of the feasible set.

Figure 3: Contour plot for the NO output of a 5-layer DNN.

####  2.1.3 Warmstarting the Problem

Seeding a MILP or NLP solver with a feasible solution of high quality can
significantly reduce the solve time. For MILP solvers, a feasible solution
with sufficiently low value can be used to prune subproblems with higher
objective value in a branch-and-bound tree search. For NLP solvers, beginning
with a feasible solution avoids the need to search for a feasible solution,
and the choice of a good solution can ensure quicker convergence to better
locally optimal solutions.

For the optimal engine design problem, we can use our collection of training
data to find parameterizations of the engine that exceed the desired torque
value but also have a low amount of emissions produced. We do not directly
work with the output data from the training set but rather with the output
from the neural network applied to the input data since this is what is
constrained.

To produce a high-quality solution is then a matter of fixing the parameters
that must remain constant for the entire drive (i.e., the compression ratio)
and then, from the training data with these fixed parameters, choosing the
remaining controls so as to minimize emissions while still exceeding the
desired torque for each time step. The remaining variable values correspond to
activation of each neuron in the neural network when applied to the input data
and hence must be set to exactly those values. This procedure is summarized in
Algorithm 1.

Data: Input DNN, compression ratio $c_{0}$, training data $\mathcal{X}$,
torque profile $\text{TorqueProfile}_{t}$, $T$

1 Select subset of training data $\tilde{\mathcal{X}}$ with compression ratio
$c_{0}$;

2 Let $\tilde{\mathcal{Y}}=\text{DNN}(\tilde{\mathcal{X}})$;

3 for  _$t=1,\ldots,T$_ do

4 Let $\text{NO},\text{CO},\text{Torque}=\tilde{\mathcal{Y}}$;

5 Let $\mathcal{T}$ be the set of times $t$ where
$\text{Torque}_{t}\geq\text{TorqueProfile}_{t}$;

6 Let $\tau=\text{argmin}_{t\in\mathcal{T}}(NO_{t}+\lambda CO_{t})$;

7 Let $x_{t}=\tilde{\mathcal{X}}_{\tau}$;

8

9 end for

10Set remaining variable values from activation of DNN when applied to $x_{t}$
for all $t$.;

Algorithm 1 Algorithm for finding a warmstart solution to (3)

###  2.2 Adversarial Attack Generation

Our next problem involves the resilience of neural networks to incorrect
classifications. Deep neural networks have the immensely useful property of
being able to uniformly approximate any function in a certain general class,
but at the cost of being fairly opaque in terms of how the underlying
machinery works. This opacity may mask unpredictable behavior that makes these
neural networks susceptible to attacks that disrupt the correct
classifications observed on training and testing data. Szegedy et al.
(szegedy2013intriguing, ) first demonstrated that almost imperceptible
perturbations of image data may lead to misclassifications, in essence
demonstrating key instabilities in neural networks. This work has led to a
number of papers (carlini2016evaluating, ; cheng2017maximum, ;
fischetti2018deep, ) developing algorithms that produce estimates of the
resilience of neural networks, estimating how close an incorrectly classified
image can be to a correctly classified one.

The central problem of interest as introduced by (szegedy2013intriguing, ) is
described as follows. Given a classifier DNN, an image $x\in[0,1]^{m}$, and a
desired classification label $l$, we have the following problem:

$\begin{array}[]{lll}\displaystyle\mathop{\mbox{minimize}}_{z}&\|x-z\|&\\\
\mbox{subject to}&\text{DNN}(z)=l&\\\ &z\in[0,1]^{m},\end{array}$ (4)

where $\|\cdot\|$ is a suitable norm. Essentially, we ask for an image of
minimal distance to image $x$ that is classified with a different label.

We add additional constraints (see, e.g., (fischetti2018deep, )) to ensure
that this is a very confident classification. Specifically, we replace the
constraint $\text{DNN}(z)=l$ with a constraint asking that the activation for
label $l$ is some factor larger than the activation for each of the other
labels; in other words, if $y_{L}$ is the final layer, we have
$y_{L,l}\geq\alpha y_{L,i}$ for $i\neq l$. Often the final layer is given by
applying a softmax function, namely,

$y_{L,i}=\sigma(z_{1},\ldots,z_{n})_{i}=\exp(z_{i})/\sum_{j=1}^{n}\exp(z_{j}),$
(5)

where the produced values $y_{L,i}$ effectively represent probabilities that
the image is a given label $i$. In this case we can represent this constraint
in terms of $z_{i}$ as

$\sigma(z)_{l}/\sigma(z)_{i}\geq\alpha\Leftrightarrow\exp(z_{l}-z_{i})\geq\alpha\Leftrightarrow
z_{l}\geq z_{i}+\log(\alpha),$

which is a linear constraint. Then if we write $y_{L}=\text{DNN}(z)$ as the
final layer prior to the softmax layer, we will have the following problem:

$\begin{array}[]{lll}\displaystyle\mathop{\mbox{minimize}}_{z}&\|x-z\|&\\\
\mbox{subject to}&y_{L}=\text{DNN}(z)&\\\ &y_{L,l}\geq
y_{L,i}+\log(\alpha),\qquad i\neq l\\\ &z\in[0,1]^{m}.\end{array}$ (6)

We observe that the formulation in (6) is easily extendable and that
constraints on the allowable perturbations can augment the formulation. For
example, we may restrict the magnitude of the perturbation to any given pixel
by restricting $|x_{i}-z_{i}|\leq\epsilon$ for all $i$. Alternatively, we may
be interested in continuous perturbations and therefore restrict
$|(x_{i}-z_{i})-(x_{j}-z_{j})|<\epsilon$ for all pixels $j$ adjacent to pixel
$i$. Constraints of this sort have been used to extend adversarial attack
optimization problems in fischetti2018deep .

###  2.3 Surrogate Modeling of Oil Well Networks

The next example of an optimization problem with an embedded neural network
involves the operation of an offshore oil platform and is taken from
(grimstad2019relu, ). The full optimization problem is reproduced in Eqns.
(7), and the associated sets are given in Table 2.

Set | Description
---|---
$N$ | Set of nodes in the network.
$N^{w}$ | Subset of well nodes in network.
$N^{m}$ | Subset of manifold nodes in the network.
$N^{s}$ | Subset of separator nodes in the network.
$E$ | Set of edges in the network.
$E^{d}$ | Subset of discrete edges that can be turned on or off.
$E^{r}$ | Subset of riser edges.
$E_{i}^{\text{in}}$ | Subset of edges entering node $i$.
$E_{i}^{\text{out}}$ | Subset of edges leaving node $i$.
$C$ | Oil, gas, and water.
Table 2: Sets used in oil well optimization problem

In this problem we have a network comprising three sets of nodes: the wells
$N^{w}$ that are the sources of oil, water, and gas; the manifolds $N^{m}$
that are connected to wells and mix incoming flows of oil, water, and gas; and
the separators $N^{s}$ that are the sinks of all the flow. Each well is
connected to each manifold by a pipeline in $E^{d}$ that may be turned on and
off. Each manifold is then connected to a unique separator by a riser in
$E^{r}$. An example network with 8 wells, 2 manifolds, and 2 separators taken
from grimstad2019relu  is depicted in Figure 4.

$\displaystyle\displaystyle\mathop{\mbox{maximize}}_{y,q,p,\Delta p}$
$\displaystyle\displaystyle\sum_{e\in E^{r}}q_{e,\text{oil}}$ (7a) subject to
$\displaystyle\sum_{e\in E_{i}^{\text{in}}}q_{e,c}=\sum_{e\in
E_{i}^{\text{out}}}q_{e,c},\;\forall c\in C,i\in N^{m}$ (7b) $\displaystyle
p_{j}=\text{DNN}_{e}(q_{e,\text{oil}},q_{e,\text{gas}},q_{e,\text{wat}},p_{i}),$
$\displaystyle\forall e=(i,j)\in E^{r}$ (7c) $\displaystyle\Delta
p_{e}=p_{i}-p_{j},$ $\displaystyle\forall e=(i,j)\in E^{r}$ (7d)
$\displaystyle-M_{e}(1-y_{e})\leq p_{i}-p_{j}-\Delta p_{e},$
$\displaystyle\forall e=(i,j)\in E^{d}$ (7e) $\displaystyle p_{i}-p_{j}-\Delta
p_{e}\leq M_{e}(1-y_{e}),$ $\displaystyle\forall e=(i,j)\in E^{d}$ (7f)
$\displaystyle\sum_{e\in E_{i}^{out}}y_{e}\leq 1,$ $\displaystyle\forall i\in
N^{w}$ (7g) $\displaystyle y_{e}q^{L}_{e,c}\leq q_{e,c}\leq y_{e}q_{e,c}^{U}$
$\displaystyle\forall c\in C,e\in E^{d}$ (7h) $\displaystyle p_{i}^{L}\leq
p_{i}\leq p_{i}^{U},$ $\displaystyle\forall i\in N$ (7i)
$\displaystyle\sum_{e\in
E_{i}^{\text{out}}}q_{e,\text{oil}}=\text{DNN}_{i}(p_{i}),$
$\displaystyle\forall i\in N^{w}$ (7j) $\displaystyle\sum_{e\in
E_{i}^{\text{out}}}q_{e,\text{gas}}=c_{e,\text{gor}}\sum_{e\in
E_{i}^{\text{out}}}q_{e,\text{oil}},$ $\displaystyle\forall i\in N^{w}$ (7k)
$\displaystyle\sum_{e\in
E_{i}^{\text{out}}}q_{e,\text{wat}}=c_{e,\text{wor}}\sum_{e\in
E_{i}^{\text{out}}}q_{e,\text{oil}},$ $\displaystyle\forall i\in N^{w}$ (7l)
$\displaystyle p_{i}=p_{i}^{s}=\text{const.},$ $\displaystyle\forall i\in
N^{w}$ (7m) $\displaystyle y_{e}\in\\{0,1\\},\forall e\in E^{d}$ (7n)

Our goal in this problem is then to optimize the flow rate of oil to each of
the sinks. Various physical and logical constraints are used to encode the
operation of this flow network. (7b) ensures that the flow into each node
matches the flow out of each node. (7e) and (7f) ensure that if a pipeline is
open, then the difference in pressure between the two nodes is actually
represented by $\Delta p_{e}$. (7g) ensures that each well routes its flow
only to a single manifold. (7h) bounds the flow rate of each material when
pipelines are active and forces the rate to zero when inactive, whereas (7i)
bounds the pressure at each node. (7k) and (7l) establish the expected ratio
of flow rates of each material.

The neural networks $\text{DNN}_{e}$ and $\text{DNN}_{i}$ appear in the
problem in (7c) and (7j) and represent nonlinear functions that predict the
separator pressure and outgoing oil flow rate based on incoming flow rate and
pressure, respectively. There is one neural network $\text{DNN}_{e}$ for each
riser edge $e$ and one neural network $\text{DNN}_{i}$ for each well $i$. Each
is trained separately; and for our particular configuration of the network, we
have 2 risers and 8 wells for a total of 10 separate neural networks that are
encoded into our problem.

The authors in grimstad2019relu  considered two different configurations of
neural networks for this problem: a “shallow” network with few layers of many
nodes and a “deep” network with many layers of fewer nodes. The shallow
configuration had two hidden layers of 20 nodes for the well networks $f_{i}$
and two hidden layers of 50 nodes for the riser networks $g_{e}$. The deep
configurations had four hidden layers of 10 nodes for the well networks and 5
layers of 20 nodes for the riser networks. All in all, this corresponds to 520
ReLU nodes for each configuration.

Figure 4: Example oil well network with 8 wells, 2 manifolds, and 2
separators, from grimstad2019relu .

##  3 Embedded Neural Network Formulation

Given that many deep learning libraries (e.g.,, TensorFlow abadi2016tensorflow
and PyTorch paszke2017automatic ) have well-developed built-in automatic
differentiation capabilities, we naturally want to see whether we can directly
embed the evaluation of the neural network into an NLP:

$\mathop{\mbox{minimize}}_{x}\;f(\text{DNN}(x),x)\quad\mbox{subject
to}\;c(\text{DNN}(x),x)=0,\;x\in X.$ (8)

Unlike the formulations in Section 4, this formulation has the advantage of
not requiring auxiliary variables for each of the internal nodes of the neural
network. That is, this formulation should scale significantly better as the
number of nodes in the neural network increases.

The modeling package JuMP DunningHuchetteLubin2017  in the programming
language Julia is a library that establishes a general framework for
representing generic optimization problems and interfacing with solvers. Of
particular use in this problem is its ability to handle nonlinear functions
with user-provided gradients.

The current release of JuMP supports only univariate user-defined nonlinear
functions, but we can represent vector-valued functions by listing each output
component separately. A single nonlinear constraint representing
$y_{i}=\text{DNN}(x)_{i}$ can be written as

        register(model, DNN_i, n, DNN_i, DNN_i_prime)))
        @NLconstraint(model, y[i] == DNN_i(x...)))),

where model is the JuMP optimization model, n is the dimension of the input,
and DNN_i and DNN_i_prime represent function evaluations of the neural
network’s $i$th output neuron and its gradient that can be provided by
TensorFlow.

For networks with many outputs, explicitly listing these commands can quickly
become unwieldy, but we can enumerate these constraints using macros. The
following code demonstrates how one can encode $y=M(x)$, where
$x\in\mathbb{R}^{n}$ and $y\in\mathbb{R}^{m}$, and f and f_p are evaluation of
the neural network and the Jacobian, respectively.

    macro DNNConstraints_grad(model, x, y, n, m, f, f_p)
      ex = Expr(:block)
      for i = 1:m
      dnn = gensym("DNN")
      push!(ex.args, :($dnn_s(z...) = $f(z...)[$i]))
        p = gensym("DNN_prime")
        push!(ex.args,:($prime(g, z...) = begin g .= $f_p(z...)[$i,:] end))
        push!(ex.args,:(register($model, $(QuoteNode(dnn_s)), $n, $dnn, $p)))
        push!(ex.args,:(@NLconstraint($model, $y[$i] == $dnn($x...))))
      end
      ex
    end

This macro can then be called to add constraints on a neural network that
takes in values from $\mathbb{R}^{20}$ and outputs values in $\mathbb{R}^{10}$
using the following command,

    @DNNConstraints_grad(model, x, y, 20, 10, DNN, DNN_prime),

where DNN and DNN_prime are now vector- and matrix-valued functions that
return the output vector and Jacobian for neural network evaluation,
respectively.

With this macro, we now have the ability to directly embed evaluation of a
neural network and its derivatives within a mathematical program. This allows
us to treat neural network evaluation as simply another function that appears
in functions and constraints so that we can use any NLP solver for our
optimization problems.

###  3.1 Convergence Behavior

We have observed in our experiments that state-of-the-art solvers express real
difficulties with convergence when applied to this formulation. We believe
that the reason has to do with the choice of activation function,
$\text{ReLU(a)}=\max(a,0)$. This activation function results in nonconvex
constraints and objective in (1); moreover, the function is nonsmooth whenever
$a=0$.

The nonsmoothness of the ReLU function has generally been considered to be of
little concern since, in practice, neural networks involving ReLU neurons have
found great success especially in problems of classification he2016deep ;
simonyan2014very ; krizhevsky2012imagenet  so that they have become an
industry standard in deep learning. Part of their success is attributed to
their ease in implementation as well as their tendency to produce sparse
activation patterns and avoid the vanishing gradient problem in training
glorot2011deep . Furthermore, some theoretical results have affirmed some
convergence guarantees under mild conditions li2017convergence ;
du2018gradient  and goodfellow2014qualitatively  have demonstrated
qualitatively that the training problem often experiences few of the issues
common in nonconvex and nonsmooth optimization.

Figure 5: Objective value, primal infeasibility, and dual infeasibility
plotted against iteration number for a sample ReLU network (left) and a sample
swish network (right) over an Ipopt solve.

These results pertain primarily to the optimization problems solved in the
training process. We have observed relatively little work in the literature
regarding optimization problems involving trained networks. Our preliminary
experiments suggest that, in this context, the expectation of good behavior
may be unfounded. Initial solves of Problem (3) using the state-of-the-art NLP
solver Ipopt wachter.biegler:06  experienced serious issues with convergence.
Figure 5 illustrates the objective value and primal and dual infeasibility for
one example solve of the problem. Qualitatively, we observe that instead of
terminating, the solver bounces about among objective values of about the same
magnitude. The reason for the lack of convergence can be explained by the fact
that one measure of convergence, the dual infeasibility, remains high
throughout the duration of the solve. This type of behavior makes this
formulation difficult for general use because there are not clear conditions
in general for checking whether the current solution at hand is “good.” We
suspect that these convergence issues can arise in part due to the
nondifferentiability of the ReLU neurons. Classical optimization algorithms
such as quasi-Newton algorithms powell1969method  often have basked in the
assumption of differentiability, and convergence results typically depend on
them.

An alternative type of activation function that has seen some success in
classifier networks ramachandran2017swish  is the swish function, which is
defined as follows:

$\text{swish}(a)=\frac{a}{1+e^{-\beta a}}.$ (9)

$\beta$ is a hyperparameter that may be learned but is typically set to 1.
Note that for $\beta=0$ this is simply a linear function and that for
$\beta\to\infty$ this approaches the ReLU function. With the activation
function in (9), all functions in (1) are twice continuously differentiable
under Assumption 1.

To ensure that failure to converge when using the ReLU network is not due to
choice of NLP solver, we perform the same solve of Problem (3) using instead
swish neurons. We plot the objective and infeasibilities in Figure 5. Instead
of failing to terminate after 3,000 iterations, Ipopt converges successfully
after only 25 iterations; and both primal infeasibility and dual infeasibility
converge rapidly to zero. As a differentiable network, the swish formulation
does not express the convergence issues that the nondifferentiable ReLU
networks appear to exhibit.

These initial results confirm the observation in (bolte2020conservative, ,
Fig. 1) that showed that ReLU networks can be nondifferentiable almost
everywhere for certain types of networks. This observation makes the use of
ReLU networks questionable as embedded constraints within optimization solvers
that rely on differentiable problem functions. Hence, in Section 4 we develop
alternative formulations that have better numerical properties.

###  3.2 Stationarity in the Embedded Formulation

In this section we will give a condition for the stationarity of the embedded
formulation as well as prove a specific form for the Jacobian of a ReLU neural
network that holds under mild genericity conditions. To handle the points of
nondifferentiability in the embedded formulation, we consider generalized
Jacobians first introduced in clarke.et.al:96 , which can be defined for all
$x$ in the domain of a function $g$ as follows:

$\partial g(x)=\text{conv}\\{\lim_{k}\nabla g(x_{k}):x_{k}\rightarrow
x,x_{k}\notin S\\},$ (10)

where $S$ is the set of points of nondifferentiability of $g$ and
$\text{conv}(A)$ is the convex hull of the set $A$. Since ReLU networks are
piecewise linear functions by construction, the generalized Jacobian will be
given by a convex hull of the Jacobians of linear functions on finitely many
regions in the neighborhood of a given point. Stationarity of $g$ at a point
$x$ is then equivalent to having $0\in\partial g(x)$.

Figure 6: Regions of a randomly instantiated neural network with two inputs
followed by 1, 2, and 3 layers of 5 ReLU neurons each.

We elaborate on what we mean by region in the context of feedforward ReLU
networks. As these neural networks are created by composing sequences of ReLU
activation functions with affine maps, we can specify the region of the
function by deciding which ReLU neurons will be active; each selection
corresponds to a different affine function. Given a set of active neurons $N$
(we will refer to these as activation patterns), we will define the associated
region $R(N)$ to be the set of points $x$ such that the value at a neuron in
$N$ prior to activation is positive and this value for any other neuron is
negative. The value prior to activation being zero corresponds to being on the
boundary of another region. $R(N)$ will be empty if no point $x$ engenders the
activation pattern $N$. The regions will tile the input space, and adjacent
regions are separated by the 0-level sets of a particular neuron.

Depictions of the regions of neural networks with 1, 2, and 3 layers can be
seen in Figure 6. We observe that the regions of a single-layer network arise
from arrangements of hyperplanes but thst as more layers are added, the
regions become more complex as they are subdivided. The total amount of
nonempty regions of both shallow and deep neural networks has been studied in
detail in montufar2014number  and pascanu2013number . However, their
definition of region as maximally connected linear regions of the piecewise
function slightly differs from ours because we distinguish adjacent regions
that have the different activation patterns but ultimately get mapped to the
same linear map.

Before we can discuss stationarity of the embedded formulation, we will need
to introduce some notation. We will identify any neuron by its layer index
$\ell$ and the index $i$ within the layer. Given a set of neurons in the
neural network $N=\\{(\ell_{1},i_{1}),\ldots,(\ell_{n},i_{n})\\}$, we define
$F_{N}$ to be the affine form given when the set of neurons in $N$ are active
and all other neurons are inactive. Put precisely, given the affine transform
in layer $\ell$, $A^{\ell}x=W^{\ell^{T}}x+b^{\ell}$, if neuron $(\ell,i)$ is
inactive, we replace column $w_{i}$ and component $b_{i}$ with zeros. We
define our new affine form as

$\hat{A}^{\ell}_{N}(x)=\hat{W}^{\ell^{T}}x+\hat{b}^{\ell},\;\hat{W}^{\ell}=W\text{diag}(\kappa^{\ell}),\;\hat{b}=\text{diag}(\kappa^{\ell})b.\,$

where $\kappa^{\ell}$ is a vector with $\kappa^{\ell}_{i}$ equal to 1 if
neuron $(\ell,i)\in B$, and zero otherwise. Then
$F_{N}(x)=\hat{A}^{L}_{N}\hat{A}^{L-1}_{N}\cdots\hat{A}^{1}_{N}x$.

We can partition the neurons of layer $\ell$ into three sets given $x$:
strictly inactive neurons $I^{\ell,-}(x)$, nonstrictly inactive neurons
$I^{\ell,0}(x),$ and active neurons $I^{\ell,+}(x)$. We define these to be the
indices $i$ of the neurons in layer $\ell$ for which
$A^{\ell}_{i}\hat{A}^{\ell-1}\cdots\hat{A}^{1}x$ is negative, zero, and
positive, respectively. The sets $I^{+}(x),I^{0}(x),I^{-}(x)$ are then defined
to be all layer-neuron pairs $(\ell,i)$, where $i\in
I^{\ell,+}(x),I^{\ell,0}(x),$ or $I^{\ell,-}(x)$, respectively. With these
definitions we observe that for any $x$, $\text{DNN}(x)=F_{I^{+}(x)}(x)$.

Note that $F_{N}$ is now smooth as a composition of affine functions, and we
can compute the Jacobian $J_{F_{N}}$ simply by the chain rule:

$J_{F_{N}}=\hat{W}^{1}\hat{W}^{2}\cdots\hat{W}^{L}.$ (11)

In light of (10), the generalized gradient at $x$ can then be written as the
convex hull of these Jacobians associated with the activation pattern for each
region in the neighborhood of $x$. Any such activation pattern must
necessarily include $I^{+}(x)$ and exclude $I^{-}(x)$ but may only include a
subset of $I^{0}(x)$. If there is a region corresponding to each set of the
form $I^{+}(x)\cup N$, where $N\subset I^{0}(x)$, we can write the generalized
gradient in a nice form given by the following proposition.

###### Proposition 1

Suppose at a point $x$ that the region $R(I^{+}(x)\cup M)$ is nonempty for
each choice of $M\subset I^{0}(x)$. Let $S$ be the set of all matrices of the
form

$W^{1}\text{diag}(\kappa^{1})\cdots W^{L}\text{diag}(\kappa^{L}),$

where

$\begin{cases}\kappa^{\ell}_{i}=1&(\ell,i)\in I^{+}(x),\\\
\kappa^{\ell}_{i}\in[0,1]&(\ell,i)\in I^{0}(x),\\\
\kappa^{\ell}_{i}=0&(\ell,i)\in I^{-}(x).\end{cases}$

Then the generalized gradient of $\text{DNN}(x)$ at $x$,
$\partial\text{DNN}(x)=S$.

###### Proof

First, we remark that any Jacobian $J_{F_{N}}$, where $N=I^{+}(x)\cup M$, is
clearly in $S$ since we can take $\kappa^{\ell}_{i}=1$ for all $(\ell,i)\in M$
and take it equal to 0 for $(\ell,i)\in I^{0}(x)\setminus M$, and the two
expressions agree. Furthermore any convex combination of the $J_{F_{B}}$
should be contained in $S$ since $S$ is a convex set. Hence
$\partial\text{DNN}(x)\subset S$.

Now given a matrix $V=W^{1}\text{diag}(\kappa^{1})\cdots
W^{L}\text{diag}(\kappa^{L})$ in $S$, we show that it is in
$\partial\text{DNN}(x)$. We enumerate
$I^{0}(x)=\\{(\ell_{1},i_{1}),\ldots,(\ell_{n},i_{n})\\}$. Consider first any
alteration $V_{1}$ of $V$ by fixing $\kappa^{\ell_{1}}_{i_{1}}$ and replacing
$\kappa^{\ell_{j}}_{i_{j}}$ by $\hat{\kappa}^{\ell_{j}}\in\\{0,1\\}$ for
$j>1$. Note that $V_{1}$ is a convex combination of the two matrices given by
replacing $\kappa^{\ell_{1}}_{i_{1}}$ with 1 and by replacing it with 0. Both
these matrices are in $\partial\text{DNN}(x)$ since this set includes all
Jacobians associated with any activation pattern on $I^{0}(x)$. Hence we must
have that $V_{1}\in\partial\text{DNN}(x)$. If we then consider fixing the
first two components and replacing the remaining such components to produce
$V_{2}$, in the same fashion we can find two matrices of the form $V_{1}$ of
which $V_{2}$ is a convex combination and hence
$V_{2}\in\partial\text{DNN}(x)$. Repeating until we have fixed all components,
we will have that $V\in\partial\text{DNN}(x)$.

Unless the region is nonempty for each of these activation patterns, the
generalized Jacobian will not necessarily take this form. As a counterexample,
consider a zero-bias two-layer network with one input variable and
$W^{1}=(1,1)$, $W^{2}=(1,-1)^{T}$. This corresponds to having
$\text{DNN}(x)=\text{ReLU}(x)-\text{ReLU}(x)=0$, and so the generalized
gradient is $\\{0\\}$ everywhere, whereas the formulation given above would
introduce many more possible gradients.

Figure 7: Regions of a zero-bias, single-layer network with two input
variables and three output variables.

For a slightly more interesting example, consider a one-layer network with two
input variables, no bias, and three output variables, namely,
$\text{DNN}(x)=\text{ReLU}(W^{T}x)=\text{ReLU}(w_{1}^{T}x)+\text{ReLU}(w_{2}^{T}x)+\text{ReLU}(w_{3}^{T}x)$.
The input domain $\mathbb{R}^{2}$ is then divided into regions by the lines
$w_{i}^{T}x=0$ for $i=1,2,3$ each intersecting at the origin, as can be seen
in Figure 7. Observe that in the neighborhood of the origin, there are only
six neighboring regions and two regions are missing, one where the output
would be $(w_{1}+w_{2})^{T}x$ and another where the output is $w_{3}^{T}$,
corresponding to activation patterns where only the first two neurons are
active and only the third neuron is active, respectively. In this scenario the
generalized gradient is similarly not as expansive as we desire.

Requiring each region $R(I^{+}(x)\cup N)$ to be nonempty seems to be an
onerous requirement, especially given there are $2^{|I^{0}(x)|}$ such regions.
However, it turns out that under some mild genericity conditions, we can
expect this to be true. These conditions arise from the theory of hyperplane
arrangements. Suppose we have a collection of hyperplanes
$H_{1},\ldots,H_{n}\subset\mathbb{R}^{d}$, where $H_{i}$ is given by
$\\{x\in\mathbb{R}^{d}:a_{i}^{T}x=b_{i}\\}$ for some vectors
$a_{i}\in\mathbb{R}^{d}$ and constants $b_{i}\in\mathbb{R}$. We say the
collection is in general position if the intersection of $m$ hyperplanes has
dimension $d-m$ if $m\leq d$ and is empty if $m>d$. If the intersection of all
hyperplanes is nonempty, it is called a central arrangement. In particular, if
$d=2$, the hyperplanes are lines in $\mathbb{R}^{2}$, and they are in general
position if no three lines intersect at a point and no two lines are parallel.
A well-known result from Zaslavsky zaslavsky1975facing  is that if we have $m$
hyperplanes in general position, the number of regions they divide
$\mathbb{R}^{d}$ into is given by $\sum_{i=0}^{d}\binom{m}{i}$. Note that if
$m\leq d$, this quantity is $2^{m}$, and we can identify each region by
choosing whether $a_{i}^{T}x<b$ or $a_{i}^{T}x>b$ for all $i$.

To ensure that the necessary regions are not empty, we will require a specific
set of hyperplanes to be in general position. Before we describe the
hyperplanes, we will prove a useful lemma.

###### Lemma 1

Suppose a central arrangement of hyperplanes $H_{1},\ldots,H_{n}$, where
$H_{i}=\\{x\in\mathbb{R}^{n}:a_{i}^{T}x=b_{i}\\}$. Then the set
$H_{1},\ldots,H_{n-1},H_{n}^{A}$, where $H_{n}^{A}=\\{x:(\sum_{i\in
A}w_{i}a_{i}+a_{n})^{T}x=\sum_{i\in A}w_{i}b_{i}+b_{n}\\}$ for any choice of
constants $w_{i}\in\mathbb{R}$, is also in general position for any
$A\subset\\{1,\ldots,n-1\\}$.

###### Proof

Because the hyperplanes intersect, we can change coordinates so that each
hyperplane intersects the origin; thus, without loss of generality, we can
write $b_{i}=0$ for all $i$. In this setting, hyperplanes being in general
position is equivalent to the normal vectors being linearly independent.
Obviously if $a_{1},\ldots,a_{n}$ are linearly independent,
$a_{1},\ldots,a_{n-1},a_{n}+\sum_{i\in A}w_{i}a_{i}$ are linearly independent.

Using this lemma, in the following proposition we can now show exactly which
hyperplanes need to be in general position.

###### Proposition 2

Suppose that the set of hyperplanes
$\\{H_{\ell,i}=\\{x:A^{\ell}_{i}\hat{A}^{\ell-1}_{I^{+}(x)}\cdots\hat{A}^{1}_{I^{+}(x)}x=0\\}:\ell=1,\ldots,L,i\in
I^{0,\ell}(x)\\}$ is in general position. Then for any choice of $N\subset
I^{0}(x)$, the set of hyperplanes
$\\{H_{\ell,i}=\\{x:A^{\ell}_{i}\hat{A}^{\ell-1}_{M}\cdots\hat{A}^{1}_{M}x=0\\}:\ell=1,\ldots,L,i\in
I^{0,\ell}(x)\\}$ is also in general position with $M=I^{+}(x)\cup N$.

###### Proof

Observe for any hyperplane
$A^{\ell}_{i}\hat{A}^{\ell-1}_{I^{+}(x)}\cdots\hat{A}^{1}_{I^{+}(x)}x=0$ and
any $j$ in $1,\ldots,\ell-1$, we consider the associated affine form
$A^{\ell}_{i}\hat{A}^{\ell-1}_{I^{+}(x)}\cdots\hat{A}^{j+1}_{I^{+}(x)}(\hat{A}^{j}_{I^{+}(x)}+\hat{A}^{j}_{N})\hat{A}^{j-1}_{I^{+}(x)}\cdots\hat{A}^{1}_{I^{+}(x)}x=0$.
Let $b^{T}x+c$ denote the affine form associated with
$A^{\ell}_{i}\hat{A}^{\ell-1}_{I^{+}(x)}\cdots\hat{A}^{j+1}_{I^{+}(x)}$, and
for each neuron $i$ in layer $j$ let $d_{i}^{T}x+e_{i}$ denote the affine form
associated with
$\hat{A}^{j}_{i}\hat{A}^{j-1}_{I^{+}(x)}\cdots\hat{A}^{1}_{I^{+}(x)}x$. Then
we can write

$A^{\ell}_{i}\hat{A}^{\ell-1}_{I^{+}(x)}\cdots\hat{A}^{j+1}_{I^{+}(x)}(\hat{A}^{j}_{I^{+}(x)}+\hat{A}^{j}_{N})\hat{A}^{j-1}_{I^{+}(x)}\cdots\hat{A}^{1}_{I^{+}(x)}x=\sum_{i\in
I^{+,j}(x)\cup N_{j}}b_{i}(d_{i}^{T}x+e_{i})+c_{i},$

and by assumption the hyperplane defined by $H=\\{x:\sum_{i\in
I^{+,j}(x)}b_{i}(d_{i}^{T}x+e_{i})+c_{i}=0\\}$ is in general position with the
hyperplanes defined by $H_{i}=\\{x:d_{i}^{T}x+e_{i}=0\\}$ for $i\in N_{j}$.
The new hyperplane is created by adding to the normal vector of $H$ a linear
combination of the normal vectors of each $H_{i}$, and so it remains in
general position with the $H_{i}$. Hence, we have replaced one $I^{+}(x)$ with
$M$. Repeating in this fashion, we can replace each $I^{+}(x)$ in each affine
form with $M$ and still remain in general position.

Now that for any choice of $N\subset I^{0}(x)$ the associated hyperplanes will
be in general position, we can be assured that the associated region is
nonempty.

###### Proposition 3

Suppose the hypothesis of Proposition 2 holds at a given $x$. Then for each
subset $N\subset I^{0}(x)$, the region $R(M)$ with $M=I^{+}(x)\cup N$ is
nonempty.

###### Proof

From Proposition 2, the hyperplanes
$\\{H_{\ell,i}=\\{x:A^{\ell}_{i}\hat{A}^{\ell-1}_{M}\cdots\hat{A}^{1}_{M}x=0\\}:\ell=1,\ldots,L,i\in
I^{0,\ell}(x)\\}$ are in general position; and since these form a central
arrangement (centered at $x$), there must be a nonempty region with
$A^{\ell}_{i}\hat{A}^{\ell-1}_{M}\cdots\hat{A}^{1}_{M}x>0$ if $(\ell,i)\in M$
and $A^{\ell}_{i}\hat{A}^{\ell-1}_{M}\cdots\hat{A}^{1}_{M}x<0$ otherwise.

The condition that the hyperplanes enumerated in Proposition 2 be in general
position is fairly mild. Certain types of neural networks will violate this
assumption: for example, zero-bias neural networks will violate this
assumption since each hyperplane will go through the origin. Additionally, we
cannot have $|I^{0}(x)|$ be larger than the dimension of the input space since
this also will have too many intersecting hyperplanes. That being said, we
generally expect a collection of randomly generated hyperplanes to be in
general position with probability 1 under most natural probability
distributions.

With these tools in hand, we can now state in full the stationarity conditions
for the embedded problem:

###### Theorem 3.1

Suppose at a given $x^{*}$, the conditions of Proposition 2 hold. Then $x^{*}$
is a stationary point of (8) if there exist nonnegative multipliers $\mu^{*}$
and matrices $\hat{W}^{1},\ldots,\hat{W}^{L}$ with
$\hat{W}^{1}\cdots\hat{W}^{L}\in\partial\text{DNN}(x)$ such that

  1. 1.

$c(\text{DNN}(x^{*}),x^{*})\leq 0$

  2. 2.

$\mu^{*^{T}}c(\text{DNN}(x^{*}),x^{*})=0$

  3. 3.

$\nabla_{x}f^{*}+\nabla_{x}c^{*}\mu^{*}+\hat{W}^{1}\cdots\hat{W}^{L}(\nabla_{y}f^{*}+\nabla_{y}c^{*}\mu^{*})=0$,

where $f^{*},c^{*}$ indicate evaluation of $f,c$ at $x^{*}$.

###### Proof

The first two conditions are standard primal feasibility and complementary
slackness conditions. To arrive at the third, we note that the gradient of the
Lagrangian $L(x,\mu)=f(\text{DNN}(x),x)+\mu c(\text{DNN}(x),x)$ on each region
(with activation pattern $B$) nearby is given by
$\nabla_{x}f+\nabla_{x}c\mu+J_{F_{B}}(\nabla_{y}f+\nabla_{y}c\mu)$ evaluated
at $x$. As $x\to x^{*}$ in a given region, and since $f$ and $c$ are smooth,
this approaches
$\nabla_{x}f^{*}+\nabla_{x}c^{*}\mu+J_{F_{B}}(\nabla_{y}f^{*}+\nabla_{y}c^{*}\mu)$.
As the generalized gradient at $x^{*}$ is given by the convex hull of all such
gradients and by Proposition 3 there is one for region $I^{+}(x)\cup A$ with
$A\subset I^{0}(x)$, it follows in a similar fashion as in the proof of
Proposition 1 that any element of the generalized gradient can be written
$\nabla_{x}f^{*}+\nabla_{x}c^{*}\mu^{*}+\hat{W}^{1}\cdots\hat{W}^{L}(\nabla_{y}f^{*}+\nabla_{y}c^{*}\mu^{*})$.
Thus condition 3 holds.

##  4 Formulating DNNs as Optimization Models

We provide two alternative formulations of ReLU DNNs in terms of optimization
models that avoid the pitfalls of the embedded formulation in Section 3.1. The
first formulation uses binary variables to model the $\max$-functions,
resulting in a mixed-integer program, building on grimstad2019relu ;
fischetti2018deep ; anderson2019strong . Our formulation differs from
cheon2020outer , which considered only inverse problems and did not exploit
the convex structure of the DNNs that arises when we lift the DNN constraint,
resulting in a nonconvex mode. The second formulation uses complementarity
constraints that can be solved as systems of nonlinear inequalities. We derive
theoretical properties of both formulations.

###  4.1 Formulating DNNs with Mixed-Integer Sets

In this section we show how general optimization problems involving DNNs, such
as (1), can be equivalently formulated as convex MIPs, extending
fischetti2018deep .

We assume that the DNN is a deep neural network with ReLU activation
functions, and we rewrite the nonconvex problem (1) as a constrained problem:

$\mathop{\mbox{minimize}}_{x,y^{L}}\;f(x,y^{L})\quad\mbox{subject
to}\;y^{L}=\text{DNN}(x),\;c(x,y^{L})\leq 0,\;x\in X,$ (12)

where $L$ is the number of layers of the DNN. Fischetti and Jo
(fischetti2018deep, ) have shown that the nonconvex constraint,
$y^{L}=\text{DNN}(x)$, can be formulated as a mixed-integer linear set in the
case of DNNs with ReLU activation functions. We let $w^{\ell}_{i}$ denote the
weights of neuron $i=1,\ldots,N_{\ell}$ at level $\ell=1,\ldots,L$ and
$b^{\ell}_{i}$ its corresponding bias. Then, the levels $\ell=1,\ldots,L$ are
computed as

$y^{\ell}_{i}=\text{ReLU}\left(w^{\ell^{T}}_{i}y^{\ell-1}+b^{\ell}_{i}\right)=\max\left(0,w^{\ell^{T}}_{i}y^{\ell-1}+b^{\ell}_{i}\right).$
(13)

We lift this constraint by introducing slack variables, $s^{\ell}_{i}$, and
binary variables, $z^{\ell}_{i}\in\\{0,1\\}$, and observe that (13) is
equivalent to the mixed-integer linear constraints

$y^{\ell}_{i}-s^{\ell}_{i}=w^{\ell^{T}}_{i}y^{\ell-1}+b^{\ell}_{i},\;0\leq
y^{\ell}_{i}\leq M^{\ell,y}_{i}(1-z^{\ell}_{i}),\;0\leq s^{\ell}_{i}\leq
M^{\ell,s}_{i}z^{\ell}_{i},\;z^{\ell}_{i}\in\\{0,1\\},$ (14)

where $M^{\ell,y}_{i},M^{\ell,s}_{i}>0$ are sufficiently large constants (if
$z^{\ell}_{i}=1$, we are on the $0$-branch of ReLU, and if $z^{\ell}_{i}=0$,
we are on the positive branch). By substituting (14) into (12) we obtain a
convex MINLP that is equivalent to the problem (1):

$\begin{array}[]{lll}\displaystyle\mathop{\mbox{minimize}}_{x,y,s,z}&f(x,y_{L})&\\\
\mbox{subject to}&c(x,y_{L})\leq 0&\\\
&y^{\ell}_{i}-s^{\ell}_{k}=w^{\ell^{T}}_{i}y^{\ell-1}+b^{\ell}_{i},&\ell=1,\ldots,L\\\
&0\leq y^{\ell}_{i}\leq M^{\ell,y}_{i}(1-z^{\ell}_{i}),\;0\leq
s^{\ell}_{i}\leq
M^{\ell,s}_{i}z^{\ell}_{i},&i=1,\ldots,N_{\ell},\ell=1,\ldots,L\\\
&y_{0}=x,\;z^{\ell}_{i}\in\\{0,1\\},x\in X.&\end{array}$ (15)

We observe that we have as many binary variables in this problem as we have
ReLU nodes. Next we show that the resulting lifted formulation results in a
tractable convex MINLP.

###### Proposition 4

Let Assumption 1 hold. Then it follows that (15) is a convex MINLP in the
sense that the continuous relaxation of (15) is a convex NLP.

###### Proof

The result follows from the convexity of $f$, $c$, and $X$ and the fact that
the remaining constraints are affine (with the exception of the integrality
restriction on $z$).

The choice of the big-$M$ constants $M^{\ell}_{i}$ in (14) can have a great
impact on the solution time of our MIP. If $M^{\ell}_{i}$ is too small, the
problem excludes solutions that should be feasible; but if it is too large,
the space that must be searched by the solver may become so large as to be too
computationally intractable for most computers.

These upper bounds are common in mixed-integer programs, and we follow the
approach in (fischetti2018deep, ) where we consider each neuron in our neural
network individually, removing all constraints on other neurons in the same or
subsequent layers. Then we set as our objective function to maximize
$y^{\ell}_{i}$ in one iteration and $s^{\ell}_{i}$ in another iteration. These
computed optimal values will serve as upper bounds $M^{\ell,y}_{i}$ and
$M^{\ell,s}_{i}$, respectively, and can be used in our model. This process of
obtaining bounds is related to the optimality-based bound-tightening technique
used in global optimization; see, for example, gleixner2017three . Since these
constants depend solely on the weights for the neural network as well as
inputs to the neural network, we need to compute these bounds only once and
may reuse them in any optimization problem involving this neural network.
Alternative methods for handling the big-$M$ constraints are considered in
(grimstad2019relu, ), and an alternative formulation with exponentially many
constraints alleviating the above concerns is presented in
(anderson2019strong, ).

With this formulation established, we may then easily pass the problem as is
to any standard MINLP solver such as Gurobi (gurobi5, ), CPLEX (cplex12, ),
Bonmin bonami2007bonmin , MINOTAUR MINOTAUR , or Baron sahinidis:96  to
compute the solution. Because MINLPs are NP-complete, however, these problems
do not scale well, and only problems involving modestly sized networks (on the
order of 100 hidden nodes) may be tractably solved. Neural networks used in
commercial settings generally involve at least thousands of hidden nodes,
resulting in thousands of binary variables, a number that typically is well
beyond the reach of commercial solvers.

###  4.2 Formulating DNNs with Complementarity Constraints

In this section we discuss an alternative formulation of (1) as a nonconvex
nonlinear program using complementarity constraints. This approach has the
advantage of scaling significantly better for problems with larger neural
networks, but with the caveat that solutions produced can be guaranteed only
to be locally optimal.

Our approaches are based on the following observations. We can rewrite the
ReLU activation function in (13) equivalently as a complementarity constraint
(using vector notation):

$y^{\ell}=\max(W^{\ell^{T}}y^{\ell-1}+b^{\ell},0)\;\Leftrightarrow\;0\leq
y^{\ell}\;\perp\;y^{\ell}\geq W^{\ell^{T}}y^{\ell-1}+b^{\ell},$ (16)

where $\perp$ means that for each component, $i$, both inequalities
$y^{\ell}_{i}\geq 0$ and
$y^{\ell}_{i}\geq[W^{\ell^{T}}y^{\ell-1}+b^{\ell}]_{i}$ are satisfied and at
least one is satisfied at equality. By replacing the ReLU function with these
complementarity constraints, we obtain a mathematical program with
complementarity constraints, which we can solve using standard NLP solvers;
see, for example, (Leyf03a, ; LeyLopNoc:06, ; FLRS:06, ; RaghBieg:05, ). The
mathematical program modeling language AMPL (fourer.gay.kernighan:93, ) allows
the modeling of complementarity constraints: for example, the above constraint
can be written as follows:

    ReLUCompl{l in Level, i in Neuron[l]}: 0 <= y[l,i] complements
    y[l,i] >= sum{j in Neuron[l-1]} W[l,i,j] y[l-1,j] + b[l,i] .

The most successful NLP solvers handle MPCCs by reformulating the
complementarity constraints in (16) by first introducing the same slack
variables as for the MINLP and lifting the formulation

$s^{\ell}=y^{\ell}-w^{\ell^{T}}_{i}y^{\ell-1}+b^{\ell},$

and then rewriting (16) equivalently as

$s^{\ell}=y^{\ell}-w^{\ell^{T}}_{i}y^{\ell-1}+b^{\ell}\;\text{and}\;0\leq
y^{\ell}\;\perp\;s^{\ell}\geq 0.$

A nonlinear optimization formulation is then given as

$s^{\ell}=y^{\ell}-w^{\ell^{T}}_{i}y^{\ell-1}+b^{\ell},\quad y^{\ell}\geq
0,\quad s^{\ell}\geq 0,\;\text{and}\;y^{\ell^{T}}s^{\ell}\leq 0,$

where the lower bound, $y^{\ell^{T}}s^{\ell}\geq 0$, is implied by the
nonnegative bounds on $y^{\ell},s^{\ell}\geq 0$, and again omitted for
numerical reasons; see (FLRS:06, ).

Because the last constraint involves the nonconvex term
$y^{\ell^{T}}s^{\ell}$, the problem as a whole is nonconvex, and standard NLP
solvers will produce only locally optimal solutions instead of globally
optimal ones. Altogether, we then have the following NLP:

$\begin{array}[]{lll}\displaystyle\mathop{\mbox{minimize}}_{x,y,s,z}&f(x,y_{L})&\\\
\mbox{subject to}&c(x,y_{L})\leq 0&\\\
&y^{\ell}_{i}-s^{\ell}_{i}=w^{\ell^{T}}_{i}y^{\ell-1}+b^{\ell}{i},&\ell=1,\ldots,L\\\
&0\leq y^{\ell}_{i},\;0\leq
s^{\ell}_{i},&i=1,\ldots,N_{\ell},\ell=1,\ldots,L\\\
&\displaystyle\sum_{\ell=1}^{L}\sum_{i=1}^{N_{\ell}}y^{\ell}_{i}s^{\ell}_{i}\leq
0,\\\ &y_{0}=x,\;x\in X.&\end{array}$ (17)

Note that we could have used the constraint $y^{\ell^{T}}s^{\ell}\leq 0$ for
each $\ell$ separately. We prefer the formulation in (17) because it has
better convergence behavior in practice; see FletLeyf:04 . Again we have a
continuous variable for every node in our neural network as well as a slack
variable for each ReLU node. The absence of integer variables, however,
produces a much more scalable problem.

We now demonstrate that the stationarity conditions for the MPCC formulation
coincide with those of the embedded formulation. In MPCC form, the
optimization problem is given by the following:

$\begin{array}[]{ll}\displaystyle\mathop{\mbox{minimize}}_{x}&f\left(x^{0},x^{L})\right)\\\
\mbox{subject to}&c\left(x^{0},x^{L})\right)\leq 0\\\ &0\leq
x^{\ell}\;\perp\;x^{\ell}\geq
W^{\ell^{T}}x^{\ell-1}+b^{\ell},\;\ell=1,\ldots,L.\end{array}$ (18)

First, we state the definition of strong stationarity for the MPCCs (18); see,
for example, SchSch:00  for its general form. We then show that the two
conditions are equivalent in our case.

######  Definition 1 (Scheel and Scholtes, SchSch:00 )

We say that $(x^{0,*},x^{1,*},\ldots,x^{L,*})$ is a strongly stationary point
of (18) if there exist multipliers $\mu^{*}\geq 0$ and
$\nu_{1}^{*},\nu_{2}^{*}$ such that the following conditions are satisfied:

$\displaystyle c^{*}\leq 0\;\text{and}\;0\leq
x^{\ell,*}\;\perp\;x_{\ell,*}\geq
W^{\ell^{T}}x^{\ell-1,*}+b^{\ell},\;\ell=1,\ldots,L$ (19a)
$\displaystyle\mu^{*,T}c^{*}=0$ (19b) $\displaystyle
0=\nabla_{x^{0}}f^{*}+\nabla_{x^{0}}c^{*}\mu^{*}+W^{1}\nu_{2}^{1,*}$ (19c)
$\displaystyle-\nu_{1}^{\ell,*}-\nu_{2}^{\ell,*}+W^{\ell+1}\nu_{2}^{\ell+1,*}=0,\;\ell=1,\ldots,L-1$
(19d)
$\displaystyle\nabla_{x^{L}}f^{*}+\nabla_{x^{L}}c^{*}\mu^{*}-\nu_{1}^{L,*}-\nu_{2}^{L,*}=0$
(19e) $\displaystyle
x^{\ell,*}_{j}>0\Rightarrow\nu^{\ell,*}_{1j}=0\;\text{and}\;x^{\ell,*}_{j}>w^{\ell^{T}}_{j}x^{\ell-1,*}+b^{\ell}_{j}\Rightarrow\nu^{\ell,*}_{2j}=0,\;\ell=1,\ldots,L$
(19f) $\displaystyle
0=x^{\ell,*}_{j}=w^{\ell,T}_{j}x^{\ell-1,*}+b_{j}\Rightarrow\nu^{\ell,*}_{1j}\geq
0,\nu_{2j}^{*}\geq 0,\;\ell=1,\ldots,L,$ (19g)

where all functions and gradients are evaluated at $(x^{*},y^{*})$, that is,
$f^{*}:=f(x^{*},y^{*})$.

Using these stationarity conditions, we can arrive at the following conditions
on the gradients of our optimization problem.

###### Proposition 5

Given a strongly stationary point satisfying (19), there exist
$\hat{W}^{1},\ldots,\hat{W}^{L}$, column scalings of $W^{1},\ldots,W^{L}$ such
that
$\nabla_{x^{0}}f^{*}+\nabla_{x^{0}}c^{*}\mu^{*}+\hat{W}^{1}\cdots\hat{W}^{L}(\nabla_{x^{L}}f^{*}+\nabla_{x^{L}}c^{*}\mu^{*})=0$.

###### Proof

To construct such matrices, we work backward starting with layer $L$ and
examine each component $j=1,\ldots,N_{L}$. Backsubstituting $\nu^{L,*}_{2}$
using (19d) and (19e), we have

$-\nu^{L-1,*}_{1}-\nu^{L-1,*}_{2}+W^{L}(\nabla_{x^{L}}f^{*}+\nabla_{x^{L}}c^{*}\mu^{*}-\nu^{L,*})=0.$

For a given component $j$ we can determine the appropriate scaling factor
$\kappa^{L}_{j}$ in three cases:

  1. 1.

$w^{L^{T}}_{j}x^{L-1,*}+b^{L}_{j}<0$: Then it follows that $x^{L,*}_{j}=0$,
$\nu^{L,*}_{2j}=0$ and $\nu_{1j}^{*}$ is not sign-constrained. It follows from
(19e) that
$\nu^{L,*}_{1j}=[\nabla_{x^{L}}f^{*}-\nabla_{x^{L}}c^{*}\mu^{*}]_{j}$ and
hence that $\kappa^{L}_{j}=0$.

  2. 2.

$w^{L^{T}}_{j}x^{L-1,*}+b^{L}_{j}>0$: Then it follows that
$x^{L,*}_{j}=w^{L,T}_{j}x^{*}+b_{j}>0$ and hence that $\nu^{L,*}_{1j}=0$.
Thus, we can set $\kappa^{L}_{j}=1$.

  3. 3.

$w^{L^{T}}_{j}x^{L,*}+b^{L}_{j}=0=x^{L,*}_{j}$: Then it follows that
$\nu^{L,*}_{1j},\nu^{L,*}_{2j}\geq 0$ from (19g). From (19e), we obtain that
$[\nabla_{x^{L}}f^{*}+\nabla_{x^{L}}c^{*}\mu^{*}]_{j}\geq\nu^{L,*}_{1j},\nu^{L,*}_{2j}\geq
0$; and given
$\nu_{1j}^{*}\in\left[0,[\nabla_{y}f^{*}+\nabla_{y}c^{*}\mu^{*}]_{j}\right]$,
we choose
$\kappa^{L}_{j}=\frac{[\nabla_{x^{L}}f^{*}+\nabla_{x^{L}}c^{*}\mu^{*}]_{j}-\nu^{L,*}_{j}}{[\nabla_{x^{L}}f^{*}+\nabla_{x^{L}}c^{*}\mu^{*}]_{j}}$,
which is in $[0,1]$.

Hence we can write
$W^{L}(\nabla_{x^{L}}f^{*}+\nabla_{x^{L}}c^{*}\mu^{*}-\nu^{L,*})=\hat{W}^{L}(\nabla_{x^{L}}f^{*}+\nabla_{x^{L}}c^{*}\mu^{*})$,
where $\hat{W}^{L}=W\text{diag}(\kappa^{L})$. Repeating this process for each
$\ell$, we arrive at
$\nabla_{x^{0}}f^{*}+\nabla_{x^{0}}c^{*}\mu^{*}+\hat{W}^{1}\cdots\hat{W}^{L}(\nabla_{x^{L}}f^{*}+\nabla_{x^{L}}c^{*}\mu^{*})=0$.

We can now prove the following theorem.

###### Theorem 4.1

Given $x^{*}$, suppose the assumptions of Proposition 2 hold. Then $x^{*}$ is
a stationary point of (8) if and only if $(x^{0,*},x^{1,*},\ldots,x^{L,*})$ is
a strongly stationary point of (18) with $x^{\ell,*}$ being the activations of
layer $\ell$ in the neural network.

###### Proof

Since the given expression lies in the generalized gradient for (8) and equals
zero, strong stationarity of the MPCC formulation clearly implies stationarity
in the embedded formulation. On the other hand, if $0$ is in the generalized
gradient for the embedded formulation, as we know the form of the generalized
gradient, we can find $\hat{W}^{1},\ldots,\hat{W}^{L}$ such that
$\nabla_{x^{0}}f^{*}+\nabla_{x^{0}}c^{*}\mu^{*}+\hat{W}^{1}\cdots\hat{W}^{L}(\nabla_{x^{L}}f^{*}+\nabla_{x^{L}}c^{*}\mu^{*})=0$.
The choice of scaling factors $\kappa$ for the matrices will then determine
the values of $\nu$ as in Proposition 5 so that the conditions for strong
stationarity hold.

##  5 Numerical Experiments

Here we present our numerical experiments with neural-network surrogates for
the three sample applications introduced in Section 2.

We assess the performance of solvers for the MIP, MPCC, and embedded
formulations of (3). We use the commercial solver CPLEX (cplex12, ) to solve
the MIP formulation of the problem and the solver Ipopt (wachter.biegler:06, )
to solve the MPCC and embedded formulations of the problem. Both solvers are
run with the default options. All experiments are performed on a single thread
on an Intel Xeon Gold 6130 CPU and 188 GB of memory.

###  5.1 Numerical Experiments with Engine Design Optimization

We have applied our models and algorithms to a collection of neural networks
with varying numbers of layers to study how well each formulation scales. The
architectures we tested all have the simple structure of an input layer with
the three input variables, followed by $n$ hidden layers of 16 nodes, and then
the output layer with the three output variables. We consider networks with 1,
3, and 5 hidden layers and train each of these networks for 20 epochs on the
simulation data produced by the simulator of (aithal2019maltese, ). We use the
adam solver in TensorFlow abadi2016tensorflow  to train the neural network.

The number of auxiliary variables for the MIP and MPCC formulations scales
with the number of ReLU neurons as well as with the number of time steps. The
last layer has no ReLU neurons since we want the final layer to be able to
take all real values, but each of the hidden layers uses ReLU neurons. We then
have $T\times(\\#\text{Hidden Layers})\times(\\#\text{Nodes Per Layer})$
additional sets of auxiliary variables.

Solving the full integer program with 1,500 time steps to optimality is
generally intractable (even for a single-layer network, this amounts to
$1500\times 16\times 1=24000$ binary variables, beyond the scope of any
general-purpose state-of-the-art solver), so we consider instead a coarser
discretization using larger time steps.

Instead of using data at each second as is presented in the original data, we
consider time intervals of 750, 500, 250, 150, 30, 15, 10, 5, and 1 seconds.
Since we have 1,500 seconds of data, this corresponds to 3, 6, 10, 50, 100,
150, 300, and 1,500 time steps for which we are evaluating our neural network.
At each of the larger time steps, our prescribed torque profile will be the
average of the prescribed torques in that interval. This approach gives us
separate problems for each choice of coarseness of the discretization and for
each choice of neural network architecture. As the number of time steps
increases, we have observed that the computed solutions converge to the
solution of the full problem.

For each instance, we run until convergence to a solution, 3 hours have
passed, or Ipopt has performed 3,000 iterations. If the solver has not
converged at the end of the experiment, we record the best solution discovered
that is feasible and that we determined to be having a constraint violation
under $10^{-6}$.

Figure 8 shows the amount of time needed until a solver found its best
solution with a time limit of 3 hours on each architecture with each time step
size. The percentage gap between the computed objective and the best-known
objective for each solver is shown in Figure 9.

Figure 8: Time until a solver found its best solution within 3 hours (10,800
seconds) on each architecture for varying values of $T$. If a runtime of
10,800 seconds is recorded, this indicates the solver failed to find a
feasible solution in the full 3 hours (if a warmstart is used, this means it
failed to find a feasible solution better than the warmstart). Figure 9:
Percentage gap in objective between final objective value and best-known
objective value for each solver on each architecture for varying values of
$T$. A value of 0 indicates the solver found the best-known objective value,
and a value of $\infty$ indicates that no solution was found.

Except for the smallest problems, CPLEX fails to converge to the optimal
solution of the MIP within the allotted time limit of 3 hours. Even worse, it
fails to even find a feasible solution for most problems with architectures
with 3 or 5 hidden layers. This result is to be expected because, except for
the smallest cases, these problems can involve thousands of binary variables.
With the warmstart, CPLEX is actually able to find solutions with
significantly better objective values for each of the configurations. It still
fails to prove optimality for any of these solutions, however, and times out
for the same set of problems as without the warmstart.

For the MPCC, on the other hand, Ipopt always finds locally optimal solutions
within the time limit, often in seconds for problems with a small number of
time steps. These solutions may be suboptimal, however, worse by up to 10% in
some instances. With the warmstart solution, Ipopt solves the MPCC in terms of
speed generally by a factor of 2. Furthermore, the objective solution tends to
have a significantly smaller objective, although there are a few exceptions.
Solving the MIP without warmstart finds a better solution than does the MPCC
formulation in only three problem instances, and in all three cases the MPCC
solution is within 0.1% of the MIP’s globally optimal solution.

In almost all cases the embedded network formulation is able to find its best
solution more quickly than does either the MIP or the MPCC formulation.
Without the warmstart, however, the embedded ReLU network formulation performs
worse in all but a few cases, with an objective value that is up to 10% larger
than the MPCC formulation objective values. With the warmstart, the objective
values are comparable except for the largest network, where in most cases
Ipopt failed to find any solutions better than the warmstart solution in the
allotted time, and so the computed objectives are worse by a small margin.
Ipopt finds the best solution to the embedded formulation faster than does the
MIP formulation, but it is often slower than the MPCC formulation. Each
iteration is faster since the formulations involve significantly fewer
variables; but because of issues with convergence, the solver often takes
significantly more steps. This situation is noticeable for the 5-layer network
when warmstarting, since Ipopt often fails to converge to any solution better
than the warmstart.

Overall, these experiments confirm our prior suspicions that using the MILP
formulation without warmstart quickly becomes computationally intractable as
the size of the neural network increases past modest architectures. Switching
to the MPCC formulation of the problem offers significant speedup at the cost
of losing global optimality, although when the MPCC formulation obtains a
worse objective value, it tends to be by only a marginal amount, and when it
does better, the improvement can be by a significant amount. The embedded
formulation, on the other hand, can provide speed and scalability but often
encounters difficulty with convergence and finds slightly worse objective
values. For all formulations, significant gains can be realized in terms of
both solving time and solution quality by providing the solver with a high-
quality warmstart solution.

###  5.2 Numerical Experiments with Adversarial Attack Generation

Next we consider each of the formulations of the adversarial attack generation
problem in (6) by considering neural networks trained on the MNIST handwritten
digit recognition data set (lecun1998mnist, ). We consider 10 different
architectures corresponding to having an input layer with $28\times 28$ nodes
followed by either 1 or 2 fully connected layers with 20, 40, 60, 80, or 100
hidden ReLU nodes each and then the output layer with 10 nodes and a softmax
activation function.

If we use the mixed-integer formulation from (15) for the DNN constraints,
then depending on our choice of norm, it becomes a mixed-integer linear
program (for $L^{1}$ or $L^{\infty}$ norms) or a mixed-integer quadratic
program (for $L^{2}$ norm) that we can solve using CPLEX. We will use the
$L^{2}$ norm in experiments.

For each neural network, we train the model by minimizing the categorical
cross-entropy loss function for 10 epochs using the 60,000 digits of training
data. Each neural network attains a test accuracy of 96–98% when tested on the
10,000 digits of testing data.

We set $\alpha=1.2$, meaning the probability for the given classification $l$
is 1.2 times higher than for any other. On each architecture, we solve the
problem for 100 digits from the training data, with the goal of finding the
closest image to the given digit that will be classified as a zero. As a
warmstart, we initialize the solution in each solve to be a digit that is
classified as the desired digit, so that the solver starts at a feasible
point. We run each iteration until convergence to the globally optimal
solution (for the MIP formulation) or to a locally optimal solution (for the
MPCC and embedded network formulations) or until one hour has passed or 3,000
iterations have occurred in Ipopt, at which point we terminate with the best
feasible solution seen so far.

In Table 3 we tabulate the results of the experiments for all three
formulations. For each neural network architecture we record the average solve
times and objectives for CPLEX and Ipopt over the 100 solves. We also record
how many times CPLEX finds the optimal solution to the MIP as well as how many
times in each formulation a feasible solution that is better than the initial
solution is found. For the embedded formulation we also include the average
objective excluding the infeasible problems since we observed a significant
difference between these cases.

Architecture | MIP | MPCC | Embedded ReLU
---|---|---|---
|  2cmAvg. |  |  |  |  |  |  |  |  |  |
Time |  2cmAvg. |  |  |  |  |  |  |  |  |  |
Obj. |  2cmNum. |  |  |  |  |  |  |  |  |  |
Opt. |  2cmNum. |  |  |  |  |  |  |  |  |  |
Feas. |  2cmAvg. |  |  |  |  |  |  |  |  |  |
Time |  2cmAvg. |  |  |  |  |  |  |  |  |  |
Obj. |  2cmNum. |  |  |  |  |  |  |  |  |  |
Feas. |  2cmAvg. |  |  |  |  |  |  |  |  |  |
Time |  2cmAvg. |  |  |  |  |  |  |  |  |  |
Obj. |  2cmAvg. |  |  |  |  |  |  |  |  |  |
Feas. |  |  |  |  |  |  |  |  |  |  |
Obj. |  2cmNum. |  |  |  |  |  |  |  |  |  |
Feas. |  |  |  |  |  |  |  |  |  |  |
784,20,10 | 1.8 | 1.79 | 100 | 100 | 0.3 | 1.80 | 100 | 35.5 | 14.67 | 1.80 | 88
784,40,10 | 76.9 | 1.75 | 100 | 100 | 1.5 | 1.77 | 100 | 28.4 | 8.64 | 1.77 | 94
784,60,10 | 2302 | 1.85 | 52 | 100 | 4.8 | 1.87 | 100 | 31.8 | 1.91 | 1.91 | 100
784,80,10 | 2571 | 2.04 | 38 | 100 | 14.2 | 2.07 | 100 | 48.1 | 17.35 | 2.08 | 86
784,100,10 | 3078 | 2.16 | 16 | 100 | 22.9 | 2.18 | 100 | 30.1 | 10.66 | 2.09 | 92
784,20,20,10 | 106.1 | 1.20 | 100 | 100 | 0.8 | 1.22 | 100 | 26.6 | 4.04 | 4.04 | 100
784,40,40,10 | 3100 | 1.90 | 16 | 100 | 4.4 | 1.77 | 100 | 38.8 | 19.97 | 1.68 | 84
784,60,60,10 | 3513 | 6.85 | 5 | 100 | 14.3 | 1.88 | 100 | 29.8 | 12.61 | 12.61 | 100
784,80,80,10 | 3587 | 27.59 | 1 | 76 | 36.1 | 1.46 | 100 | 27.8 | 3.83 | 3.83 | 100
784,100,100,10 | 3600 | 54.83 | 0 | 42 | 57.0 | 1.78 | 100 | 35.1 | 1.83 | 1.83 | 100
Table 3: Comprehensive results for solves of the MIP formulation using CPLEX,
the MPCC formulation using Ipopt, and the embedded ReLU formulation also using
Ipopt. Times are recorded in seconds.

The average time to convergence for the CPLEX solves of the MIP formulation is
significantly higher than that of the Ipopt solves of the MPCC formulation.
Most CPLEX solves for the larger networks time out at 1 hour before finding an
optimal solution; for the smaller networks, it can take 10 to 100 times longer
to converge to a solution. Ipopt, on the other hand, takes less than a minute
on average to converge to a locally optimal solution of the MPCC formulation.
The embedded formulations are fairly uniform in how long they take to find
their best solution, which makes sense given that the problem does not change
in size with the number of neurons in the network.

The longer solve times of the MIP formulation have the advantage of eventually
uncovering better solutions the majority of the time even if they are not
provably optimal. For all the single hidden layer architectures and the two
smallest double hidden layer architectures, the MIP formulation almost always
produces a better solution than the MPCC formulation does. These solutions
generally offer only a marginal improvement on the order of 1–2% difference in
the objective. For perturbations of this magnitude, the differences are
essentially imperceptible.

For the largest networks the MIP formulation fares poorly and cannot find a
better feasible solution than the initial solution in all cases for the two
largest networks. The feasible solutions it does find are of noticeably lower
quality: of all instances from the three largest networks, CPLEX found the
best solution in only 13 cases out of 300. The MIP solutions for the larger
networks are visually worse, as can be seen in a comparison between the
perturbed images generated from CPLEX and Ipopt in Figure 11 with respect to
the original images in Figure 10.

Figure 10: Suite of 10 digits from the MNIST training data for which
adversarial attacks are generated.

Figure 11: Perturbed images produced after 1 hour of computation when using
CPLEX to solve the MIP formulation (left) and Ipopt to solve the NLP
formulation (right) with given digits in Figure 10 for the neural network with
2 hidden layers of 100 nodes.

When compared with the MIP and MPCC formulations, the embedded formulation,
while often faster and more scalable, can be more unstable in the sense that
while using it, Ipopt can fail to find feasible solutions even though for this
problem they certainly exist. The ReLU networks are more prone to failure and
exhibit the poor convergence behavior as in Figure 5 leading to failure to
find solutions in about 5% of instances. In a few cases the objective value of
the found solutions is also noticeably worse.

###  5.3 Numerical Experiments with Oil Well Networks

Next we consider solving the oil well problem using each of our formulations.
We again used CPLEX for the MIP formulation and solved the problem with a time
limit of 1 hour, recording the best time for each neural network
configuration. This problem involves binary variables even in the
complementary constraint formulation, so we needed to use an MINLP solver and
elected to use the solver Bonmin bonami2007bonmin . We used the default
branch-and-bound scheme with Ipopt as the NLP solver. Because no MINLP solver
is currently supported on JuMP for handling both integer variables and user-
defined nonlinear functions, we could not solve the problem as stated using
the embedded neural network formulations.

The results for the solves on the full problem are presented in Table 4. We
observe that in this small set of problems, the MPCC formulation outperforms
the MIP formulation in both network configurations in terms of time taken as
well as objective value produced. CPLEX fails to even produce a solution for
the deep neural network configuration. These are hard problems; indeed, in 3
of the 4 solves performed, the solutions found are not proven optimal.

| Shallow | Deep
---|---|---
Solver | Time (s) | Objective | Time (s) | Objective
CPLEX | $3388^{*}$ | 1.271 | $3600^{*}$ | NaN
bonmin | $1865^{*}$ | 1.283 | 1275.3 | 1.304
Table 4: Time to find the best solution for each instance of the full oil well
problem with a time limit of 1 hour. The * indicates that the solver could not
prove the solution optimal.

To obtain a problem that did not have binary variables, we considered fixing
each of the binary variables randomly to either 0 or 1 while still ensuring
feasibility. If we represent the neural networks using complementary
constraints or using the embedded formulation, the problems become standard
NLPs that we can directly solve using Ipopt as before. In our setup we
considered 10 different configurations of the binary variables, and we solved
them using CPLEX for the MIP formulation and Ipopt for the rest. We ran each
instance until convergence or 1 hour had passed, and we recorded the best
solution found.

| Shallow | Deep
---|---|---
Instance | MIP | MPCC | Embedded | MIP | MPCC | Embedded
1 | 17.5 | 1.8 | 10.3 | $3600^{*}$ | 1.6 | 16.7
2 | 22.0 | 1.2 | 8.2 | $3600^{*}$ | 1.4 | 19.4
3 | 23.9 | 1.5 | 7.5 | $3600^{*}$ | 1.7 | 14.2
4 | 40.3 | 1.5 | 6.2 | $3600^{*}$ | 1.4 | 7.2
5 | 19.0 | 1.4 | 4.0 | $3600^{*}$ | 1.6 | 13.0
6 | 21.0 | 2.1 | 3.0 | $3600^{*}$ | 1.7 | 7.1
7 | 169.4 | 1.7 | 5.4 | $3600^{*}$ | 1.2 | 14.1
8 | 14.9 | 1.5 | 5.7 | $3600^{*}$ | 1.4 | 7.3
9 | 86.3 | 1.9 | 4.3 | $3600^{*}$ | 1.6 | 12.1
10 | 15.5 | 1.7 | 5.8 | $3600^{*}$ | 1.5 | 5.8
Table 5: Time until a solver finds the best feasible solution for each
formulation and for both shallow and deep neural networks on each particular
fixing of the oil well problem. The * indicates that the solver could not
prove the found solution optimal.

The results for the experiments performed with the binary variables fixed to
arbitrary configurations are presented in Table 5, which details the solve
times, and in Table 6, which details the objective values produced. We observe
that, overall, the MPCC formulation clearly outperforms all the other
formulations in terms of speed in that it is able to find solutions in about 2
seconds where the embedded network formulations take 10 to 20 seconds and the
MIP formulation can take much longer. In fact, the MIP formulation fails to
find a single solution for any of the instances using the deep neural network
configuration, which suggests that optimization problems with deeper networks
are more difficult than networks with shallow networks even if the number of
nodes is the same.

| Shallow | Deep
---|---|---
Instance | MIP | MPCC | ReLU NN | MIP | MPCC | Embedded
1 | 1.231 | 1.231 | 1.231 | NaN | 1.248 | 1.183
2 | 1.228 | 1.228 | 1.228 | NaN | 1.228 | 1.155
3 | 1.219 | 1.219 | 1.219 | NaN | 1.229 | 1.274
4 | 1.264 | 1.264 | 1.264 | NaN | 1.274 | 1.255
5 | 1.182 | 1.182 | 1.182 | NaN | 1.192 | 1.253
6 | 1.252 | 1.252 | 1.252 | NaN | 1.256 | 1.266
7 | 1.229 | 1.229 | 1.229 | NaN | 1.244 | 1.215
8 | 1.207 | 1.207 | 1.207 | NaN | 1.201 | 1.225
9 | 1.239 | 1.239 | 1.239 | NaN | 1.249 | 1.282
10 | 1.263 | 1.263 | 1.263 | NaN | 1.276 | 1.243
Table 6: Objective value of the best feasible solution for each formulation
for both shallow and deep neural networks on each particular fixing of the oil
well problem. NaN indicates that no feasible solution was found.

In terms of objective value, for the shallow network all formulations were
able to find the global optimum in every instance relatively quickly. For the
deep networks the different formulations outperformed each other on different
instances, leaving no clear winner. The MPCC formulation solves each instance
the quickest, because of the relatively small size of each of the neural
networks.

Unlike the other two problems, we observed for these instances that the
embedded ReLU network formulation did not have trouble with convergence and
that, in spite of nondifferentiability, the dual infeasibility was brought
down to zero, signifying convergence in all instances. We postulate that the
reason for the success here as compared with the other instances may be due to
the simplicity of the networks: 8 of the networks are single input and single
output, and so the difficulties that might emerge in multiple dimensions do
not appear.

##  6 Conclusions and Future Research

We have presented three alternative formulations of ReLU deep-neural network
constraints as a mixed-integer problem, an optimization problem with
complementarity constraints, and a problem with the neural network directly
embedded. The MIP and MPCC formulations can be viewed as lifted formulations,
and we have shown that the lifting convexifies optimization problems with deep
neural network constraints in the case of the mixed-integer formulation. We
have also presented a warmstart technique that uses training data of the
neural network to construct good initial solutions. We have compared the three
formulations on three examples arising in the design of engines, the design of
images that “fool” a given classifier, and the assignment of flow in an oil
well network. Each formulation has its advantages. We have shown the MIP
formulation to be useful in finding an optimal solution, but at the cost of a
particularly long solve time. We observed that the new complementarity
constraint formulation generally outperforms the mixed-integer formulation in
terms of solution time but may not find the optimal solution, although it
often comes close. We also observed that the embedded neural network
formulation has the advantage of being scalable and quick to solve but has
difficulties with convergence related to the nondifferentiability of the ReLU
activation function (which could be rectified by using a smooth activation
function like the swish function). The experiments reported in this report
portray the versatility of each of the different formulations in handling the
breadth of optimization problems that may involve neural networks.

## References

  * [1] M. Abadi, P. Barham, J. Chen, Z. Chen, A. Davis, J. Dean, M. Devin, S. Ghemawat, G. Irving, M. Isard, et al.  TensorFlow: A system for large-scale machine learning.  In 12th $\\{$USENIX$\\}$ Symposium on Operating Systems Design and Implementation ($\\{$OSDI$\\}$ 16), pages 265–283, 2016.
  * [2] S. M. Aithal and P. Balaprakash.  MaLTESE: Large-scale simulation-driven machine learning for transient driving cycles.  In High Performance Computing, pages 186–205, Cham, 2019. Springer International Publishing.
  * [3] R. Anderson, J. Huchette, C. Tjandraatmadja, and J. P. Vielma.  Strong mixed-integer programming formulations for trained neural networks.  In International Conference on Integer Programming and Combinatorial Optimization, pages 27–42, 2019.
  * [4] P. Belotti.  Couenne: A user’s manual.  Technical report, FICO, 2020.
  * [5] J. Bolte and E. Pauwels.  Conservative set valued fields, automatic differentiation, stochastic gradient methods and deep learning.  Mathematical Programming, pages 1–33, 2020.
  * [6] P. Bonami and J. Lee.  BONMIN user’s manual.  Numer Math, 4:1–32, 2007.
  * [7] N. Carlini and D. Wagner.  Towards evaluating the robustness of neural networks.  arXiv preprint arXiv:1608.04644, 2016.
  * [8] C.-H. Cheng, G. Nührenberg, and H. Ruess.  Maximum resilience of artificial neural networks.  arXiv preprint arXiv:1705.01040, 2017.
  * [9] M.-S. Cheon.  An outer-approximation guided optimization approach for constrained neural network inverse problems.  arXiv preprint arXiv:2002.10404, 2020.
  * [10] L. Clarke, J. Linderoth, E. Johnson, G. Nemhauser, R. Bhagavan, and M. Jordan.  Using OSL to improve the computational results of a MIP logistics model.  EKKNEWS, 16, 1996.
  * [11] A. Delarue, R. Anderson, and C. Tjandraatmadja.  Reinforcement learning with combinatorial actions: An application to vehicle routing.  Advances in Neural Information Processing Systems, 33, 2020.
  * [12] S. S. Du, X. Zhai, B. Poczos, and A. Singh.  Gradient descent provably optimizes over-parameterized neural networks.  In International Conference on Learning Representations, 2018.
  * [13] I. Dunning, J. Huchette, and M. Lubin.  Jump: A modeling language for mathematical optimization.  SIAM review, 59(2):295–320, 2017.
  * [14] I. Dunning, J. Huchette, and M. Lubin.  Jump: A modeling language for mathematical optimization.  SIAM Review, 59(2):295–320, 2017.
  * [15] S. Dutta, S. Jha, S. Sankaranarayanan, and A. Tiwari.  Output range analysis for deep feedforward neural networks.  In NASA Formal Methods Symposium, pages 121–138. Springer, 2018\.
  * [16] M. Fischetti and J. Jo.  Deep neural networks and mixed integer linear optimization.  Constraints, 23(3):296–309, 2018.
  * [17] R. Fletcher and S. Leyffer.  Solving mathematical program with complementarity constraints as nonlinear programs.  Optimization Methods and Software, 19(1):15–40, 2004.
  * [18] R. Fletcher, S. Leyffer, D. Ralph, and S. Scholtes.  Local convergence of SQP methods for mathematical programs with equilibrium constraints.  SIAM Journal on Optimization, 17(1):259––286, 2006.
  * [19] R. Fourer, D. M. Gay, and B. W. Kernighan.  AMPL: A Modeling Language for Mathematical Programming.  The Scientific Press, 1993.
  * [20] A. M. Gleixner, T. Berthold, B. Müller, and S. Weltge.  Three enhancements for optimization-based bound tightening.  Journal of Global Optimization, 67(4):731–757, 2017.
  * [21] X. Glorot, A. Bordes, and Y. Bengio.  Deep sparse rectifier neural networks.  In Proceedings of the fourteenth international conference on artificial intelligence and statistics, pages 315–323. JMLR Workshop and Conference Proceedings, 2011.
  * [22] I. J. Goodfellow, O. Vinyals, and A. M. Saxe.  Qualitatively characterizing neural network optimization problems.  arXiv preprint arXiv:1412.6544, 2014.
  * [23] B. Grimstad and H. Andersson.  ReLU networks as surrogate models in mixed-integer linear programs.  Computers & Chemical Engineering, 131:106580, 2019.
  * [24] Gurobi Optimization, Inc.  Gurobi Optimizer Reference Manual, Version 5.0, 2012.
  * [25] K. He, X. Zhang, S. Ren, and J. Sun.  Deep residual learning for image recognition.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 770–778, 2016.
  * [26] IBM Corp.  IBM Ilog CPLEX V12.1: User’s Manual for CPLEX, 2009.
  * [27] G. Katz, C. Barrett, D. L. Dill, K. Julian, and M. J. Kochenderfer.  Reluplex: An efficient SMT solver for verifying deep neural networks.  In International Conference on Computer Aided Verification, pages 97–117. Springer, 2017.
  * [28] E. B. Khalil, A. Gupta, and B. Dilkina.  Combinatorial attacks on binarized neural networks.  arXiv preprint arXiv:1810.03538, 2018.
  * [29] A. Krizhevsky, I. Sutskever, and G. E. Hinton.  Imagenet classification with deep convolutional neural networks.  Advances in neural information processing systems, 25:1097–1105, 2012.
  * [30] Y. LeCun.  The MNIST database of handwritten digits.  http://yann. lecun. com/exdb/mnist/, 1998.
  * [31] S. Leyffer.  Mathematical programs with complementarity constraints.  SIAG/OPT Views-and-News, 14(1):15–18, 2003.
  * [32] S. Leyffer, G. Lopez-Calva, and J. Nocedal.  Interior methods for mathematical programs with complementarity constraints.  SIAM Journal on Optimization, 17(1):52–77, 2006.
  * [33] Y. Li and Y. Yuan.  Convergence analysis of two-layer neural networks with ReLU activation.  Advances in Neural Information Processing Systems, 30:597–607, 2017\.
  * [34] A. Mahajan, S. Leyffer, J. Linderoth, J. Luedtke, and T. Munson.  MINOTAUR: a toolkit for solving mixed-integer nonlinear optimization.  wiki-page, 2011.  http://wiki.mcs.anl.gov/minotaur.
  * [35] G. F. Montufar, R. Pascanu, K. Cho, and Y. Bengio.  On the number of linear regions of deep neural networks.  Advances in Neural Information Processing Systems, 27:2924–2932, 2014.
  * [36] T. Papalexopoulos, C. Tjandraatmadja, R. Anderson, J. P. Vielma, and D. Belanger.  Constrained discrete black-box optimization using mixed-integer programming.  arXiv preprint arXiv:2110.09569, 2021.
  * [37] R. Pascanu, G. Montufar, and Y. Bengio.  On the number of response regions of deep feed forward networks with piece-wise linear activations.  arXiv preprint arXiv:1312.6098, 2013.
  * [38] A. Paszke, S. Gross, S. Chintala, G. Chanan, E. Yang, Z. DeVito, Z. Lin, A. Desmaison, L. Antiga, and A. Lerer.  Automatic differentiation in PyTorch.  2017\.
  * [39] M. Powell.  A method for nonlinear constraints in minimization problems in optimization.  In R. Fletcher, editor, Optimization. Academic Press, 1969.
  * [40] N. V. Queipo, R. T. Haftka, W. Shyy, T. Goel, R. Vaidyanathan, and P. K. Tucker.  Surrogate-based analysis and optimization.  Progress in aerospace sciences, 41(1):1–28, 2005.
  * [41] A. Raghunathan and L. T. Biegler.  An interior point method for mathematical programs with complementarity constraints (MPCCs).  SIAM Journal on Optimization, 15(3):720–750, 2005.
  * [42] P. Ramachandran, B. Zoph, and Q. V. Le.  Swish: a self-gated activation function.  arXiv preprint arXiv:1710.05941, 7:1, 2017.
  * [43] M. Ryu, Y. Chow, R. Anderson, C. Tjandraatmadja, and C. Boutilier.  Caql: Continuous action q-learning.  In International Conference on Learning Representations, 2019.
  * [44] N. V. Sahinidis.  BARON: A general purpose global optimization software package.  Journal of Global Optimization, 8(2):201–205, 1996.
  * [45] H. Scheel and S. Scholtes.  Mathematical program with complementarity constraints: Stationarity, optimality and sensitivity.  Mathematics of Operations Research, 25:1–22, 2000.
  * [46] A. M. Schweidtmann and A. Mitsos.  Deterministic global optimization with artificial neural networks embedded.  Journal of Optimization Theory and Applications, 180(3):925–948, 2019.
  * [47] T. Serra and S. Ramalingam.  Empirical bounds on linear regions of deep rectifier networks.  In AAAI, pages 5628–5635, 2020.
  * [48] K. Simonyan and A. Zisserman.  Very deep convolutional networks for large-scale image recognition.  arXiv preprint arXiv:1409.1556, 2014.
  * [49] C. Szegedy, W. Zaremba, I. Sutskever, J. Bruna, D. Erhan, I. Goodfellow, and R. Fergus.  Intriguing properties of neural networks.  arXiv preprint arXiv:1312.6199, 2013.
  * [50] M. Tawarmalani and N. V. Sahinidis.  Convexification and Global Optimization in Continuous and Mixed-Integer Nonlinear Programming: Theory, Algorithms, Software, and Applications.  Kluwer Academic Publishers, Boston MA, 2002.
  * [51] V. Tjeng, K. Xiao, and R. Tedrake.  Evaluating robustness of neural networks with mixed integer programming.  arXiv preprint arXiv:1711.07356, 2017.
  * [52] A. Wächter and L. T. Biegler.  On the implementation of a primal-dual interior point filter line search algorithm for large-scale nonlinear programming.  Mathematical Programming, 106(1):25–57, 2006.
  * [53] T. Zaslavsky.  Facing up to arrangements: Face-count formulas for partitions of space by hyperplanes: Face-count formulas for partitions of space by hyperplanes, volume 154.  American Mathematical Soc., 1975.

The submitted manuscript has been created by UChicago Argonne, LLC, Operator
of Argonne National Laboratory (“Argonne”). Argonne, a U.S. Department of
Energy Office of Science laboratory, is operated under Contract No. DE-
AC02-06CH11357. The U.S. Government retains for itself, and others acting on
its behalf, a paid-up nonexclusive, irrevocable worldwide license in said
article to reproduce, prepare derivative works, distribute copies to the
public, and perform publicly and display publicly, by or on behalf of the
Government. The Department of Energy will provide public access to these
results of federally sponsored research in accordance with the DOE Public
Access Plan. http://energy.gov/downloads/doe-public-access-plan.