# Is there a high quality nonlinear programming solver for Python?

I have several challenging non-convex global optimization problems to solve. Currently I use MATLAB's Optimization Toolbox (specifically, fmincon() with algorithm='sqp'), which is quite effective. However, most of my code is in Python, and I'd love to do the optimization in Python as well. Is there a NLP solver with Python bindings that can compete with fmincon()? It must

• be able to handle nonlinear equality and inequality constraints
• not require the user to provide a Jacobian.

It's okay if it doesn't guarantee a global optimum (fmincon() does not). I'm looking for something that robustly converges to a local optimum even for challenging problems, and even if it's slightly slower than fmincon().

I have tried several of the solvers available through OpenOpt and found them to be inferior to MATLAB's fmincon/sqp.

Just for emphasis I already have a tractable formulation and a good solver. My goal is merely to change languages in order to have a more streamlined workflow.

Geoff points out that some characteristics of the problem may be relevant. They are:

• 10-400 decision variables
• 4-100 polynomial equality constraints (polynomial degree ranges from 1 to about 8)
• A number of rational inequality constraints equal to about twice the number of decision variables
• The objective function is one of the decision variables

The Jacobian of the equality constraints is dense, as is the Jacobian of the inequality constraints.

-
David, this is now unfortunately a completely different question :) The difference between local minimum and global is the subject of a potential infinite number of PhDs, and by the No Free Lunch Theorem, any solver that is good for one general global optimization problem is provably bad for another. I might suggest that you start by considering formulation options (Is there a mixed integer form? Does a convex approximation exist?) – Aron Ahmadia Nov 30 '11 at 19:34
David, Aron makes a good point. Formulation is definitely key in terms of obtaining numerical solutions of non-convex NLPs, let alone obtaining good solutions quickly. It may be worth considering alternative formulations, and then using the structure of those formulations to guide your choice of solver. Using a solver that exploits any structure (such as sparsity, multi-stage stochastic programming, using constraints to generate cuts) that you can induce in your problem is key to getting good solutions. – Geoff Oxberry Dec 11 '11 at 5:22
@DavidKetcheson: Since you have a formulation that you want to use, could you at least comment on the characteristics of your formulation? Is the Jacobian of the Lagrangian dense or sparse? Roughly how many variables does it have? It does you no good for us to recommend software that implements solution methods that are ill-suited for your problem, and that's the only reason people are talking about formulations in the first place. – Geoff Oxberry Dec 12 '11 at 17:56

fmincon(), as you mentioned, employs several strategies that are well-known in nonlinear optimization that attempt to find a local minimum without much regard for whether the global optimum has been found. If you're okay with this, then I think you have phrased the question correctly (nonlinear optimization).

The best package I'm aware of for general nonlinear optimization is IPOPT[1]. Apparently Matthew Xu maintains a set of Python bindings to IPOPT, so this might be somewhere to start.

[1]: Andreas Wachter is a personal friend, so I may be a bit biased.

-
 Andreas does good work, but his solver also requires Jacobian matrix information (or at the very least, sparsity information for the Jacobian matrix). When you say that you want a solver that does not require a Jacobian matrix, do you mean that you want a solver that does not require you to provide the Jacobian matrix analytically (so that a finite-difference calculation would suffice) or do you want a solver that does not require Jacobian matrix information at all (which would limit you to derivative-free optimization methods)? – Geoff Oxberry♦ Dec 11 '11 at 5:17 Good catch. I mean the former; I've updated the question. – David Ketcheson Dec 11 '11 at 8:10 I was finally able to apply IPOPT to my problem using sage.openopt.org. It's great! – David Ketcheson Mar 16 '12 at 19:46

I work in a lab that does global optimization of mixed-integer and non-convex problems. My experience with open source optimization solvers has been that the better ones are typically written in a compiled language, and they fare poorly compared to commercial optimization packages.

If you can formulate your problem as an explicit system of equations and need a free solver, your best bet is probably IPOPT, as Aron said. Other free solvers can be found on the COIN-OR web site. To my knowledge, the nonlinear solvers do not have Python bindings provided by the developers; any bindings you find would be third-party. In order to obtain good solutions, you would also have to wrap any nonlinear, convex solver you found in appropriate stochastic global optimization heuristics, or in a deterministic global optimization algorithm such as branch-and-bound. Alternatively, you could use Bonmin or Couenne, both of which are deterministic non-convex optimization solvers that perform serviceably well compared to the state-of-the-art solver, BARON.

If you can purchase a commercial optimization solver, you might consider looking at the GAMS modeling language, which includes several nonlinear optimization solvers. Of particular mention are the interfaces to the solvers CONOPT, SNOPT, and BARON. (CONOPT and SNOPT are convex solvers.) A kludgey solution that I've used in the past is to use the Fortran (or Matlab) language bindings to GAMS to write a GAMS file and call GAMS from Fortran (or Matlab) to calculate the solution of an optimization problem. GAMS has Python language bindings, and a very responsive support staff willing to help out if there's any trouble. (Disclaimer: I have no affiliation with GAMS, but my lab does own a GAMS license.) The commercial solvers should be no worse than fmincon; in fact, I'd be surprised if they weren't a lot better. If your problems are sufficiently small in size, then you may not even need to purchase a GAMS license and licenses to solvers, because an evaluation copy of GAMS may be downloaded from their web site. Otherwise, you would probably want to decide which solvers to purchase in conjunction with a GAMS license. It's worth noting that BARON requires a mixed-integer linear programming solver, and that licenses for the two best mixed-integer linear programming solvers CPLEX and GUROBI are free for academics, so you might be able to get away with just purchasing the GAMS interfaces rather than the interfaces and the solver licenses, which can save you quite a bit of money.

This point bears repeating: for any of the deterministic non-convex optimization solvers I've mentioned above, you need to be able to formulate the model as an explicit set of equations. Otherwise, the non-convex optimization algorithms won't work, because all of them rely on symbolic analysis to construct convex relaxations for branch-and-bound-like algorithms.

UPDATE: One thought that hadn't occurred to me at first was that you could also call the Toolkit for Advanced Optimization (TAO) and PETSc using tao4py and petsc4py, which would have the potential added benefit of easier parallelization, and leveraging familiarity with PETSc and the ACTS tools.

UPDATE #2: Based on the additional information you mentioned, sequential quadratic programming (SQP) methods are going to be your best bet. SQP methods are generally considered more robust than interior point methods, but have the drawback of requiring dense linear solves. Since you care more about robustness than speed, SQP is going to be your best bet. I can't find a good SQP solver out there written in Python (and apparently, neither could Sven Leyffer at Argonne in this technical report). I'm guessing that the algorithms implemented in packages like SciPy and OpenOpt have the basic skeleton of some SQP algorithms implemented, but without the specialized heuristics that more advanced codes use to overcome convergence issues. You could try NLopt, written by Steven Johnson at MIT. I don't have high hopes for it because it doesn't have any reputation that I know of, but Steven Johnson is a brilliant guy who writes good software (after all, he did co-write FFTW). It does implement a version of SQP; if it's good software, let me know.

I was hoping that TAO would have something in the way of a constrained optimization solver, but it doesn't. You could certainly use what they have to build one up; they have a lot of the components there. As you pointed out, though, it'd be much more work for you to do that, and if you're going to that sort of trouble, you might as well be a TAO developer.

With that additional information, you are more likely to get better results calling GAMS from Python (if that's an option at all), or trying to patch up the IPOPT Python interface. Since IPOPT uses an interior point method, it won't be as robust, but maybe Andreas' implementation of an interior point method is considerably better than Matlab's implementation of SQP, in which case, you may not be sacrificing robustness at all. You'd have to run some case studies to know for sure.

You're already aware of the trick to reformulate the rational inequality constraints as polynomial inequality constraints (it's in your book); the reason this would help BARON and some other nonconvex solvers is that it can use term analysis to generate additional valid inequalities that it can use as cuts to improve and speed up solver convergence.

Excluding the GAMS Python bindings and the Python interface to IPOPT, the answer is no, there aren't any high quality nonlinear programming solvers for Python yet. Maybe @Dominique will change that with NLPy.

UPDATE #3: More wild stabs at finding a Python-based solver yielded PyGMO, which is a set of Python bindings to PaGMO, a C++ based global multiobjective optimization solver. Although it was created for multiobjective optimization, it can also be used to single objective nonlinear programming, and has Python interfaces to IPOPT and SNOPT, among other solvers. It was developed within the European Space Agency, so hopefully there's a community behind it. It was also released relatively recently (November 24, 2011).

-

APM Python is a free optimization toolbox that has interfaces to APOPT, BPOPT, IPOPT, and other solvers. It provides first (Jacobian) and second (Hessian) information to the solvers and provides an optional web-interface to view results.

We've done a couple benchmark tests and found that the combination of APOPT (active set method) and IPOPT (interior point method) can solve a large percentage of benchmark problems. There are a number of example problems that are included with the download zip file. The one that you'll probably want to start with is the Hock Schittkowski #71 problem. It is the simplest example and demonstrates how to solve constrained optimization problems.

There is a browser interface and an API to Python / MATLAB. The API to Python is a single script (apm.py) that is available for download from the apmonitor.com homepage. Once the script is loaded into a Python code, it gives the ability to solve problems of:

• Nonlinear equations
• Differential and algebraic equations
• Least squares model fitting
• Moving horizon estimation
• Nonlinear model predictive control
• etc.

For the new user, the APM Python software has a Google Groups forum where a user can post questions. There are bi-weekly webinars that showcase optimization problems in operations research and engineering.

Below is an example of an optimization problem (hs71.apm).

Model
Variables
x[1] = 1, >=1, <=5
x[2] = 5, >=1, <=5
x[3] = 5, >=1, <=5
x[4] = 1, >=1, <=5
End Variables

Equations
x[1] * x[2] * x[3] * x[4] > 25
x[1]^2 + x[2]^2 + x[3]^2 + x[4]^2 = 40

minimize  x[1] * x[4] * (x[1]+x[2]+x[3]) + x[3]
End Equations
End Model


The optimization problem is solved with the following Python script:

from apm import *
server = 'http://xps.apmonitor.com'

# Application name
app = 'eqn'

# Clear previous application
apm(server,app,'clear all')

# Load model file

# Option to select solver (1=APOPT, 2=BPOPT, 3=IPOPT)
apm_option(server,app,'nlc.solver',3)

# Solve on APM server
solver_output = apm(server,app,'solve')

# Display solver output
print solver_output

# Retrieve results
results = apm_sol(server,app)

# Display results
print '--- Results of the Optimization Problem ---'
print results

# Display Results in Web Viewer
url = apm_var(server,app)
print "Opened Web Viewer: " + url


APM Python is a free web-service for optimization. The optimization problems are solved on remote servers and results are returned to the local Python script. We recently added parallel processing support for both MATLAB and Python.

-
John, I see that APM Python is freely available, but I can't figure out from looking at the package whether it contains solvers that it uses locally or it requires a connection to the AP Monitor website to do the computations. I'm curious as to which. – Aron Ahmadia Feb 6 '12 at 5:19
Aron, the MATLAB or Python scripts require an internet connection to the APM servers to solve the optimization problems. This has a number of advantages and disadvantages. On the positive side, a web-service for optimization allows for cross-platform compatibility, free access to some commercial solvers, and software upgrades that are transparent to the user. On the downside, APM is not as flexible as some open-source alternatives but is designed for industrial users who favor a turn-key solution for optimization applications. – John Hedengren Feb 21 '12 at 19:29

Though this does not entirely answer your question, I author a Python package for nonlinear programming named NLPy. The most recent version may be retrieved from https://github.com/dpo/nlpy

I must stress that NLPy is research-grade and the solvers included are by no means as robust as more seasoned codes like IPOPT. Moreover, they currently require that Jacobians be provided. That being said, the point of NLPy is to provide the tools needed for researchers to assemble custom solvers if they need to. At any rate, I'll be interested to hear your comments offline if you do give it a try. You may also find the related packages https://github.com/dpo/pykrylov and https://github.com/dpo/pyorder useful. Currently, the documentation of NLPy is definitely lacking. The other two should be reasonable.

-

There's cvxmod, a Python wrapper around Stephen Boyd's convex optimization software. It's part of the Sage package.

-
 But the OP is asking about a non-convex optimization problem. – Alejandro Nov 30 '11 at 23:10 The OP is asking about a non-convex optimization problem, but all of the solvers mentioned so far are only guaranteed to find epsilon-optimal solutions to convex optimization problems without additional metaheuristics (multistart, or other stochastic global optimization algorithms that call on deterministic, nonlinear, convex optimization solvers) or branch-and-bound-like algorithms (such as branch-and-bound, branch-and-cut, and branch-and-reduce) that require relaxations of the objective function and constraints. This answer is no worse than any of the others mentioned as of Dec 11th. – Geoff Oxberry♦ Dec 11 '11 at 5:07 Geoff, how can I apply cvxmod to a non-convex problem? – David Ketcheson Dec 11 '11 at 8:08 I haven't used the software, but in theory, like any other convex solver, you'd use it to find locally optimal solutions much like you're using fmincon now (which is also a convex solver). One way to use it would be multistart. Generate a list of points to be used as initial guesses for your convex solver. For each point used as a guess, record the solution returned by the solver. The point that corresponds to the minimum objective function value over all solutions returned is the best approximation to the global optimum. – Geoff Oxberry♦ Dec 11 '11 at 8:23 Another approach would be to generate successive convex relaxations that approximate your non-convex problem, and solve each subproblem with a convex solver. Solutions to relaxations will give lower bounds on the minimum objective function value for your problem, and any solution returned by the convex solver on the non-convex problem will yield an upper bound on the minimum objective function value. Such a strategy forms the basis for many deterministic non-convex optimization algorithms. – Geoff Oxberry♦ Dec 11 '11 at 8:26
show 3 more comments

pyomo is a full GAMS/AMPL-like modeling environment for optimization in python. It is extremely powerful, has interfaces to all solvers that are supported by AMPL, and generates Jacobians etc. automatically. However, due to it running in a 'virtual python environment', it might not be trivial to link it to existing code.

-

What about scipy.fmin_slsqp ?

http://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.fmin_slsqp.html

-
 Thanks, but that is one that I tried (through OpenOpt, which provides an additional interface to it). It was never better than fmincon/sqp and failed in many cases where the latter succeeded. – David Ketcheson Dec 1 '11 at 5:48 Update: I tried this one directly from SciPy. It fails even on problems where fmincon is able to consistently find the global optimum in a few seconds. – David Ketcheson Dec 3 '11 at 18:15

Since MATLAB has a JIT compiler while CPython does not yet (at least, until pypy gets full numpy support). It seems like you want a free solver that outperforms commercially produced fmincon. Isn't it too much?

IIRC among commercial NLP solvers, only snopt has provided a Python API until now (although it's rather ugly).

Which OpenOpt solvers have you tried? How many variables and constraints do you have in your nonconvex task?

You could try IPOPT through OpenOpt / Funcdesigner API without installation on OpenOpt Sage server (pay attention to the "switch from sage to python" picture).

Algencan, also connected to OpenOpt, says it solve nonconvex probs rather good (unfortunately it's absent in OO sage server). On the other hand, there is no gradient-based solver capable of yielding a guaranteed solution for nonconvex problems. Even some convex NLPs with 2 variables can be constructed where any solver will fail due to machine roundoff errors, e.g. $10^{300}(x-0.1)^2 + 10^{-300}(y-0.2)^2$ with starting guess like $(x,y) = (1,1)$.

-
If you read carefully, I'm just asking for something with similar robustness to fmincon. It doesn't need to be better, and it can even be slower. – David Ketcheson Dec 12 '11 at 11:30

for global problems you could be interested in http://openopt.org/interalg and other openopt global solvers (http://openopt.org/GLP ) for local optimization openopt also provides variety of solvers: http://openopt.org/NLP

-
 Yes, I tried some of those, but none measured up to fmincon. – David Ketcheson Jun 21 '12 at 16:25

PyGMO, http://pagmo.sourceforge.net/pygmo/index.html, contains several solvers, providing the same interface to them. IPOPT and scipy slsqp are included in case you compile the code and download/install the third party code indipendently.

As a bonus, parallel use of the solver is made really easy (multistart) via the archipelago class!!

-