# Python solvers for mixed-integer nonlinear constrained optimization

I want to minimize a black box function $f(x)$, which takes a 8$\times$3 matrix of non-negative integers as input. Each row specifies a variable, whereas each column specifies a certain time period so that $x_{ij}$ is the $i$th variable in the $j$th time period.

$f(x)$ is the activation function for a convolutional filter in the 2. convolutional layer of a convolutional neural network. It is therefore expected to be nonlinear and nonconvex. The input is subject to constraints as seen below, where $\text{sgn}(x)$ is the sign function as defined at Wikipedia.

\begin{aligned} \text{Objective:} \hspace{1cm} & \text{minimize} \hspace{0.2cm} f(x)\\ \text{Constraints:} \hspace{1cm} & \sum_{i=3}^{6}\sum_{j=1}^{3} x_{ij} = 10 & \\ & x_{3j} + x_{5j} \geq x_{1j} &\forall j = 1,2,3\\ & x_{4j} + x_{6j} \geq x_{2j} &\forall j = 1,2,3\\ & x_{7j} \leq 15 &\forall j = 1,2,3\\ & x_{8j} \leq 15 &\forall j = 1,2,3\\ & \text{sgn}(x_{7j}) = \text{sgn}(x_{3j}) &\forall j = 1,2,3 \\ & \text{sgn}(x_{8j}) = \text{sgn}(x_{4j}) &\forall j = 1,2,3 \\ & x_{ij} \in \mathbb{N_0} &\forall i = 1,2,\cdots,8 \hspace{0.2cm} \forall j = 1,2,3 \end{aligned}

Can anyone recommend any Python packages that would be able to solve this problem? Any commercial software with an interface to Python and a free academic license/evaluation period would also be great.

EDIT: It should be noted that the optimization does not have to find a global minimum (although that is, of course, preferred).

• I might not have explained the problem well enough (see the updated question). It is nonlinear and I am certain that the optimal solution will have D's unequal to 0.
– pir
Commented Jun 7, 2015 at 8:01
• @felbo: I'm guess that your nonlinearity consists of the restrictions on signs of certain pairs of unknowns. If you are asking how to solve this problem with off-the-shelf software packages, one approach would be breaking each of those pair-sign constraints into pairs of inequalities, e.g. $x_{7,j} \gt 0$ and $x_{3,j} \gt 0$ OR $x_{7,j} \lt 0$ and $x_{3,j} \lt 0$ (trusting you to fill in the details about how signs for zero values should be treated). But perhaps you are not as much interested in this one problem as in a broader survey of integer optimization software.
– hardmath
Commented Jun 7, 2015 at 14:00
• @hardmath: I'm not interested in a broader survey of integer optimization software, but in solving this exact problem. Breaking the pair-sign constraints into pairs of inequalities sounds like a great idea if that can be used as input for off-the-shelf packages! The nonlinearity I mentioned in the question is not regarding the constraints, but regarding the function $f(x)$. Do you have any recommendation on off-the-shelf software for black box optimization after your reformulation of the constraints?
– pir
Commented Jun 7, 2015 at 14:12
• @felbo: I see. Given the difficult (nonlinear, nonconvex) behavior of your function $f(x)$, I would focus on effciently generating all feasible points $(x_{i,j})$ in your domain and evaluating the function at those points. Off-the-shelf software would not be likely to do anything more intelligent than such exhaustive search.
– hardmath
Commented Jun 7, 2015 at 14:19
• What if I could start the optimization with values $x_0$, where it is known that $f(x_0)$ is close to a local minimum? So it would just have to find closest minimum. Perhaps Ipopt as suggested by @DavidKetcheson would be a good bet?
– pir
Commented Jun 7, 2015 at 15:00

I want to minimize a black box function $f(x)$, which takes a 8×3 matrix of non-negative integers as input.

In general, this sort of problem should be solved with a derivative-free MINLP (or if the function is linear, an MILP) solver. A very cursory glance at the literature suggests that the algorithm DFL, presented in a JOTA paper (also see preprint) could work; there might be other algorithms that also solve derivative-free MINLPs. The authors' implementation is in Fortran 90, so you would need to write some wrappers.

Generally speaking, though, you need some solver that can solve:

• black-box/derivative-free problems (I assume here that derivatives are not available, otherwise it would not be "black box")
• that are also mixed-integer

Since your problem contains no continuous decision variables, exhaustive sampling, as proposed by @hardmath, is another option that is probably easier to implement if you'd rather not write Python wrappers to a Fortran package (I wouldn't blame you).

However, IPOPT would not work, because it requires at least gradient information, and solves continuous problems, not mixed-integer problems. Even ignoring that estimating both gradient and Hessian information with finite-difference-type approaches is likely to yield poor outcomes due to numerical errors in these approximations, IPOPT will only solve the continuous relaxation of the problem. IPOPT (like any other continuous optimization solver) would have to be augmented with branch-and-cut or branch-and-bound-type methods, which is a great deal of work. If you're going to forge ahead and use a solver that assumes you have gradient information (even though you don't), you would be much better served using an MINLP solver (e.g., BARON, Couenne, Bonmin) that would at least give you integer-feasible solutions.

• Why is it mixed-integer? All the parameters are integers.
– pir
Commented Jun 18, 2015 at 15:11
• Purely integer programs are a subset of mixed-integer programs. Since methods for solving integer programs include approaches like relaxing the integer variables to continuous variables, then solving using a branch-and-bound (or branch-and-cut, or other) scheme, from an implementation perspective, implementing an integer programming solver frequently means mixed-integer programs can also be solved. Commented Jun 18, 2015 at 22:15
• Can you recommend an MINLP solver that has an interface to something more accessible than Fortran and C? A solver with a Python interface would be amazing.
– pir
Commented Mar 7, 2016 at 13:45
• @pir: I know there are black-box MINLP solvers implemented in MATLAB, but I don't know if they've been released. If you're just looking for MINLP solvers in the standard case where you know your objective and constraints (and these are not black box, plus you can compute derivatives), you might try interfacing to them via modeling languages like GAMS, AMPL, JuMP (Julia), or Pyomo (Python). There should be a MATLAB interface as well, I just don't remember what it is off the top of my head. Commented Mar 7, 2016 at 20:51
• Thanks! I am able to compute the first-order derivatives, but the function would still be black box in the sense that it cannot be formulated in a concise mathematical way (output of a neural network). Do you have any recommendations for which solver to use of the ones available in Pyomo? (software.sandia.gov/downloads/pub/pyomo/…)
– pir
Commented Mar 7, 2016 at 22:19

I agree with all of the answers provided here but I wanted to supplement with a Python implementation. We developed the Python GEKKO package for solving similar problems. We're also working on machine learning functions that may be able to combine a convolutional neural network with this constrained mixed-integer problem as a single optimization. Here is a potential solution with Python GEKKO (>0.2rc4).

import numpy as np
from gekko import GEKKO
m = GEKKO()
ni = 8
nj = 3
x = [[m.Var(lb=0,integer=True) for j in range(nj)] for i in range(ni)]
s = 0
for i in range(ni):
for j in range(nj):
s += x[i][j]
m.Equation(s==10)
m.Equations([x[2][j]+x[4][j]>=x[0][j] for j in range(nj)])
m.Equations([x[3][j]+x[5][j]>=x[1][j] for j in range(nj)])
for j in range(nj):
x[6][j].upper=15
x[7][j].upper=15
m.Equations([(m.sign3(x[6][j])==m.sign3(x[2][j])) for j in range(nj)])
m.Equations([(m.sign3(x[7][j])==m.sign3(x[3][j])) for j in range(nj)])
m.Obj(sum([m.tanh(x[i][0]) for i in range(ni)]))
m.options.SOLVER=1
m.solve()
print(x)


When you switch to IPOPT with m.options.SOLVER=3, it gives a non-integer solution, as expected:

---------------------------------------------------
Solver         :  IPOPT (v3.12)
Solution time  :   0.109899999995832      sec
Objective      :   7.933869530630442E-009
Successful solution
---------------------------------------------------

[[0.0], [0.2838769309], [0.2838769309]]
[[0.0], [0.2838769309], [0.2838769309]]
[[3.9669348007e-09], [0.99734110398], [0.99734110398]]
[[3.96693473e-09], [0.99734110398], [0.99734110398]]
[[0.0], [0.52439022369], [0.52439022369]]
[[0.0], [0.52439022369], [0.52439022369]]
[[0.0], [0.69439174379], [0.69439174379]]
[[0.0], [0.69439174379], [0.69439174379]]


However, when you switch to the APOPT MINLP solver with m.options.SOLVER = 1, it gives an integer solution:

---------------------------------------------------
Solver         :  APOPT (v1.0)
Solution time  :   3.900000000430737E-002 sec
Objective      :   0.000000000000000E+000
Successful solution
---------------------------------------------------

[[0.0], [2.0], [1.0]]
[[0.0], [1.0], [1.0]]
[[0.0], [0.0], [0.0]]
[[0.0], [0.0], [0.0]]
[[0.0], [2.0], [1.0]]
[[0.0], [1.0], [1.0]]
[[0.0], [0.0], [0.0]]
[[0.0], [0.0], [0.0]]


As you can see, the solution is fast but it is likely because the problem size is small and the objective function is not correct. For small problems, an exhaustive search may be the right solution. However, if you scale-up to larger problems with more integer variables or include the convolutional neural network model, each NLP branch and bound sub-problem will take longer to solve or may not converge.

• Branch and Bound Commented Oct 20, 2023 at 5:05
• still Gekko returns float representation of result.x, thus can convert it to integers with res_milp.x=[int(xi) for xi in res_milp.x] Commented Oct 20, 2023 at 5:07
• That is correct. Gekko has an integer tolerance where it can find a solution within 1e-2 of the integer value. This can be adjusted with the solver option minlp_integer_tol. Solutions are always reported as floating point values. Solver options are available here: apmonitor.com/wiki/index.php/Main/OptionApmSolver Commented Oct 20, 2023 at 12:53

The feasible region can be broken up into several "convex" integer lattice subdomains by setting out the specific combinations of sign-agreement pairs. [That is, each sign agreement between $x_{7,j}$ and $x_{3,j}$, resp. between $x_{8,j}$ and $x_{4,j}$ can occur in one of two ways: both positive or both zero. Given the three "time" indices $j=1,2,3$, there are $2^6 = 64$ "components" which are the integer lattice solutions of linear constraints.]

When the region described by linear constraints is too complicated or large for direct generation of all feasible points, a popular technique is to generate a random sample of the points, e.g. by a Monte Carlo Markov Chain (MCMC) algorithm.

There is a Python implementation of MCMC samplers as pymc by Christopher Fonnesbeck, Anand Patil and David Huard. Although targeted to Bayesian estimation, the sampling portion of the toolkit could be used to generate "random" feasible points from each subdomain described above.

The minimum would then be approximated by the least value of $f(x)$ attained over all the sampled points.

Here is a sketch of how one might generate all the feasible points for these particular constraints.

Note that the "odd" row and "even" row indices $i$ are linked only by the sum-to-ten constraint at top. We can rearrange the unknowns from the $8\times 3$ matrix $x$ into a $4\times 6$ matrix:

$$\begin{pmatrix} x_{1,1} & x_{1,2} & x_{1,3} & x_{2,1} & x_{2,2} & x_{2,3} \\ x_{3,1} & x_{3,2} & x_{3,3} & x_{4,1} & x_{4,2} & x_{4,3} \\ x_{5,1} & x_{5,2} & x_{5,3} & x_{6,1} & x_{6,2} & x_{6,3} \\ x_{7,1} & x_{7,2} & x_{7,3} & x_{8,1} & x_{8,2} & x_{8,3} \end{pmatrix}$$

Each column of non-negative integers $(a,b,c,d)^T$ in this form obeys the same constraints:

$$a \le b + c$$

$${b = 0} \Leftrightarrow {d = 0}$$

$$d \le 15$$

We can add the constraint $b+c \le 10$ deduced from the combined middle two row sums in the new form being exactly ten, and we get a finite set of possibilities $S \subset \mathbb{N}_0^4$ from which all six columns are chosen.

Furthermore the sum-to-ten constraint for the two middle rows allows us to assign parts that add to $10$ for the six columns, thereby generating all the feasible points for this problem. There are $35$ partitions of $10$ in not more than six parts, and $\binom{15}{5} = 3003$ weak compositions of $10$ in exactly six summands, so this appears to be a tractable task computationally.

Added: A program written to carry out these calculations found some 8 trillion feasible solutions. More precisely there are 8,054,848,375,116 black box function evaluations needed to exhaust the search space.

• Sounds interesting! Note that $x_{i,j} \in \mathbb{N_0}$ $\forall i = 1,2,..,8$ $\forall j = 1,2,3$, allowing only two ways of sign agreement. Based on my understanding there would then be $2^{2 \cdot 3} = 64$ components.
– pir
Commented Jun 7, 2015 at 19:33
• How do you know that each integer lattice subdomain is "convex"?
– pir
Commented Jun 7, 2015 at 19:51
• The mention of "convex" was perhaps too terse/cryptic. Any nonempty region of $\mathbb{R}^n$ defined by linear equalities and inequalities is necessarily convex. What I meant by "convex" is the nonmepty intersection of the integer lattice $\mathbb{Z}^n$ with a convex subset of $\mathbb{R}^n$. Unfortunately there seems to be no succinct name for this. Because these sets are disconnected if they contain more than one point, generally these are not truly convex sets (though their theory and the application of MCMC algorithms is tied to convexity).
– hardmath
Commented Jun 8, 2015 at 0:09
• The constraint is not across all $x_{ij}$, but only for $i = 3,\cdots,6$.
– pir
Commented Jun 8, 2015 at 16:24
• @hardmath: The terminology I'd use is that the continuous relaxation of the problem is convex. Commented Jun 17, 2015 at 1:07