# Good languages/packages for interior point optimization with non-linear constraints?

I'm currently using Python's scipy.optimize package to perform parameter estimation for a system of 10 ODEs. I have some observed data, and I'm trying to find the set of parameters which makes the ODE solution best fit the data in the least squares sense. (The ODE, by the way, doesn't have an analytic solution, so I'm not supplying a gradient.) Each parameter is constrained to an interval, and there are also a couple inequality constraints. There are about 15-30 parameters in total that I'd like to optimize. I tend to get decent fits with SLSQP and so-so results with 'trust-constr' in scipy.optimize.minimize. The former method tends to be faster (30-90 seconds vs 5-10 minutes for the latter) and also tends to give a lower squared error. However, it is preferred for the problem I'm working on to find a minimum (or minima) which are not on the boundaries. Since Python's 'trust-constr' method hasn't been working all that well for me, I'm looking for alternatives.

Can anyone provide a recommendation for a language and/or package that would have a suitable interior-point (or similar) method? I've heard some good things about Julia, though in briefly reading the JuliaOpt and JuMP pages, I'm not sure how good the support is for problems with inequality constraints. Performance is something that is semi-important as I plan to run hundreds of optimizations and I'd like to be able to do this in a reasonable amount of time. Any ideas or recommendations would be greatly appreciated.

• You should know that some ODE solvers have "sensitivity analysis" features that allow you to compute your gradients along with the solution to the ODE. This feature is available, for example, in the Sundials package and in Julia's differential equations package. Nov 9 '20 at 22:55
• @Brian Borchers Thanks for your comment, but it's not clear to me how this relates to my question... Nov 12 '20 at 21:03
• You wrote "The ODE, by the way, doesn't have an analytic solution, so I'm not supplying a gradient." My comment explained how you could provide a gradient rather than depending on the optimization routine to use finite difference approximation. This would be helpful with whatever optimization solver you might use. Nov 12 '20 at 22:25
• If you could do the nonlinear constraints as penalty terms / augmented Lagrangian, try scipy least_squares. Feb 19 at 9:34

See also this ugly plot of 3 NLopt algorithms $$\times$$ a dozen synthetic 8d test functions.