Python solvers for mixed-integer nonlinear constrained optimization

Question

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).

Answer 1

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.

Answer 2

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.

Answer 3

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.