Python constrained non-linear optimization

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 npM = 10N = 3Q = 10C = 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)]) - vbounds = [(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.0from mystic.solvers import diffev2from mystic.monitors import VerboseMonitormon = 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.434718Generation 10 has Chi-Squared: -1.733787Generation 20 has Chi-Squared: -1.859787Generation 30 has Chi-Squared: -1.860533Generation 40 has Chi-Squared: -1.860533Generation 50 has Chi-Squared: -1.860533Generation 60 has Chi-Squared: -1.860533Generation 70 has Chi-Squared: -1.860533Generation 80 has Chi-Squared: -1.860533Generation 90 has Chi-Squared: -1.860533Generation 100 has Chi-Squared: -1.860533Generation 110 has Chi-Squared: -1.860533Generation 120 has Chi-Squared: -1.860533Generation 130 has Chi-Squared: -1.860533Generation 140 has Chi-Squared: -1.860533Generation 150 has Chi-Squared: -1.860533Generation 160 has Chi-Squared: -1.860533Generation 170 has Chi-Squared: -1.860533Generation 180 has Chi-Squared: -1.860533Generation 190 has Chi-Squared: -1.860533Generation 200 has Chi-Squared: -1.860533Generation 210 has Chi-Squared: -1.860533STOP("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*x2bounds = [(0,1e6)]*4# with penalty='penalty' applied, solution is:xs = [0.72759093, 0.35964857, 37.69901188, 240.0]ys = 5804.3762083from mystic.symbolic import generate_constraint, generate_solvers, simplifyfrom mystic.symbolic import generate_penalty, generate_conditionsequations = """-x0 + 0.0193*x2 <= 0.0-x1 + 0.00954*x2 <= 0.0-pi*x2**2*x3 - (4/3.)*pi*x2**3 + 1296000.0 <= 0.0x3 - 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.

scipy has a spectacular package for constrained non-linear optimization.

You can get started by reading the optimize doc, but here's an example with SLSQP:

minimize(func, [-1.0,1.0], args=(-1.0,), jac=func_deriv, constraints=cons, method='SLSQP', options={'disp': True})

As others have commented as well, the SciPy minimize package is a good place to start. We also have a review of many other optimization packages in the Python Gekko paper (see Section 4). I've included an example below (Hock Schittkowski #71 benchmark) that includes an objective function, equality constraint, and inequality constraint in Scipy.optimize.minimize.

import numpy as npfrom scipy.optimize import minimizedef objective(x):    return x[0]*x[3]*(x[0]+x[1]+x[2])+x[2]def constraint1(x):    return x[0]*x[1]*x[2]*x[3]-25.0def constraint2(x):    sum_eq = 40.0    for i in range(4):        sum_eq = sum_eq - x[i]**2    return sum_eq# initial guessesn = 4x0 = np.zeros(n)x0[0] = 1.0x0[1] = 5.0x0[2] = 5.0x0[3] = 1.0# show initial objectiveprint('Initial SSE Objective: ' + str(objective(x0)))# optimizeb = (1.0,5.0)bnds = (b, b, b, b)con1 = {'type': 'ineq', 'fun': constraint1} con2 = {'type': 'eq', 'fun': constraint2}cons = ([con1,con2])solution = minimize(objective,x0,method='SLSQP',\                    bounds=bnds,constraints=cons)x = solution.x# show final objectiveprint('Final SSE Objective: ' + str(objective(x)))# print solutionprint('Solution')print('x1 = ' + str(x[0]))print('x2 = ' + str(x[1]))print('x3 = ' + str(x[2]))print('x4 = ' + str(x[3]))

Here is the same problem with Python Gekko:

from gekko import GEKKOm = GEKKO()x1,x2,x3,x4 = m.Array(m.Var,4,lb=1,ub=5)x1.value = 1; x2.value = 5; x3.value = 5; x4.value = 1m.Equation(x1*x2*x3*x4>=25)m.Equation(x1**2+x2**2+x3**2+x4**2==40)m.Minimize(x1*x4*(x1+x2+x3)+x3)m.solve(disp=False)print(x1.value,x2.value,x3.value,x4.value)

There is also a more comprehensive discussion thread on nonlinear programming solvers for Python if SLSQP can't solve your problem. My course material on Engineering Design Optimization is available if you need additional information on the solver methods.