Python constrained non-linear optimization

Question

What's the recommended package for constrained non-linear optimization in python ?

The specific problem I'm trying to solve is this:

I have an unknown X (Nx1), I have M (Nx1) u vectors and M (NxN) s matrices.

max [5th percentile of (ui_T*X), i in 1 to M]
st 
0<=X<=1 and
[95th percentile of (X_T*si*X), i in 1 to M]<= constant

When I started out the problem I only had one point estimate for u and s and I was able to solve the problem above with cvxpy.

I realized that instead of one estimate for u and s, I had the entire distribution of values so I wanted to change my objective function so that I could use the entire distribution. The problem description above is my attempt to include that information in a meaningful way.

cvxpy cannot be used to solve this, I've tried scipy.optimize.anneal, but I can't seem to set bounds on the unknown values. I've looked at pulp too but it doesnt allow nonlinear constraints.

Answer 1

While the SLSQP algorithm in scipy.optimize.minimize is good, it has a bunch of limitations. The first of which is it's a QP solver, so it works will for equations that fit well into a quadratic programming paradigm. But what happens if you have functional constraints? Also, scipy.optimize.minimize is not a global optimizer, so you often need to start very close to the final results.

There is a constrained nonlinear optimization package (called mystic) that has been around for nearly as long as scipy.optimize itself -- I'd suggest it as the go-to for handling any general constrained nonlinear optimization.

For example, your problem, if I understand your pseudo-code, looks something like this:

import numpy as np

M = 10
N = 3
Q = 10
C = 10

# let's be lazy, and generate s and u randomly...
s = np.random.randint(-Q,Q, size=(M,N,N))
u = np.random.randint(-Q,Q, size=(M,N))

def percentile(p, x):
    x = np.sort(x)
    p = 0.01 * p * len(x)
    if int(p) != p:
        return x[int(np.floor(p))]
    p = int(p)
    return x[p:p+2].mean()

def objective(x, p=5): # inverted objective, to find the max
    return -1*percentile(p, [np.dot(np.atleast_2d(u[i]), x)[0] for i in range(0,M-1)])


def constraint(x, p=95, v=C): # 95%(xTsx) - v <= 0
    x = np.atleast_2d(x)
    return percentile(p, [np.dot(np.dot(x,s[i]),x.T)[0,0] for i in range(0,M-1)]) - v

bounds = [(0,1) for i in range(0,N)]

So, to handle your problem in mystic, you just need to specify the bounds and the constraints.

from mystic.penalty import quadratic_inequality
@quadratic_inequality(constraint, k=1e4)
def penalty(x):
  return 0.0

from mystic.solvers import diffev2
from mystic.monitors import VerboseMonitor
mon = VerboseMonitor(10)

result = diffev2(objective, x0=bounds, penalty=penalty, npop=10, gtol=200, \
                 disp=False, full_output=True, itermon=mon, maxiter=M*N*100)

print result[0]
print result[1]

The result looks something like this:

Generation 0 has Chi-Squared: -0.434718
Generation 10 has Chi-Squared: -1.733787
Generation 20 has Chi-Squared: -1.859787
Generation 30 has Chi-Squared: -1.860533
Generation 40 has Chi-Squared: -1.860533
Generation 50 has Chi-Squared: -1.860533
Generation 60 has Chi-Squared: -1.860533
Generation 70 has Chi-Squared: -1.860533
Generation 80 has Chi-Squared: -1.860533
Generation 90 has Chi-Squared: -1.860533
Generation 100 has Chi-Squared: -1.860533
Generation 110 has Chi-Squared: -1.860533
Generation 120 has Chi-Squared: -1.860533
Generation 130 has Chi-Squared: -1.860533
Generation 140 has Chi-Squared: -1.860533
Generation 150 has Chi-Squared: -1.860533
Generation 160 has Chi-Squared: -1.860533
Generation 170 has Chi-Squared: -1.860533
Generation 180 has Chi-Squared: -1.860533
Generation 190 has Chi-Squared: -1.860533
Generation 200 has Chi-Squared: -1.860533
Generation 210 has Chi-Squared: -1.860533
STOP("ChangeOverGeneration with {'tolerance': 0.005, 'generations': 200}")
[-0.17207128  0.73183465 -0.28218955]
-1.86053344078

mystic is very flexible, and can handle any type of constraints (e.g. equalities, inequalities) including symbolic and functional constraints. I specified the constraints as "penalties" above, which is the traditional way, in that they apply a penalty to the objective when the constraint is violated. mystic also provides nonlinear kernel transformations, which constrain solution space by reducing the space of valid solutions (i.e. by a spatial mapping or kernel transformation).

As an example, here's mystic solving a problem that breaks a lot of QP solvers, since the constraints are not in the form of a constraints matrix. It's optimizing the design of a pressure vessel.

"Pressure Vessel Design"

def objective(x):
    x0,x1,x2,x3 = x
    return 0.6224*x0*x2*x3 + 1.7781*x1*x2**2 + 3.1661*x0**2*x3 + 19.84*x0**2*x2

bounds = [(0,1e6)]*4
# with penalty='penalty' applied, solution is:
xs = [0.72759093, 0.35964857, 37.69901188, 240.0]
ys = 5804.3762083

from mystic.symbolic import generate_constraint, generate_solvers, simplify
from mystic.symbolic import generate_penalty, generate_conditions

equations = """
-x0 + 0.0193*x2 <= 0.0
-x1 + 0.00954*x2 <= 0.0
-pi*x2**2*x3 - (4/3.)*pi*x2**3 + 1296000.0 <= 0.0
x3 - 240.0 <= 0.0
"""
cf = generate_constraint(generate_solvers(simplify(equations)))
pf = generate_penalty(generate_conditions(equations), k=1e12)


if __name__ == '__main__':

    from mystic.solvers import diffev2
    from mystic.math import almostEqual
    from mystic.monitors import VerboseMonitor
    mon = VerboseMonitor(10)

    result = diffev2(objective, x0=bounds, bounds=bounds, constraints=cf, penalty=pf, \ 
                     npop=40, gtol=50, disp=False, full_output=True, itermon=mon)

    assert almostEqual(result[0], xs, rel=1e-2)
    assert almostEqual(result[1], ys, rel=1e-2)

Find this, and roughly 100 examples like it, here: https://github.com/uqfoundation/mystic.

I'm the author, so I am slightly biased. However, the bias is very slight. mystic is both mature and well-supported, and is unparalleled in capacity to solve hard constrained nonlinear optimization problems.