Investigate constrained optimization for multiple springpot models

[moustachio-belvedere]

We should consider building in some clever functionality which embeds model constraints (particularly for multi-springpot fractional models) into the optimisation constraints.

This issue is related to https://github.com/JuliaRheology/RHEOS.jl/issues/62 as it should help us get rid of the modified Mittag-Leffler function dependency.

Whereas before, the goal was limited to just automatically notice we had springpots, and thus placing 0<a<1 bounds on the springpot parameters, I was reading through the NLopt documentation again today, and the library does allow for constraint relationships between parameters:

https://github.com/JuliaOpt/NLopt.jl

It would be very cool if we had a system that automatically embedded all the constraint information from the models into the optimization process.

As a concrete example, for the fractional Kelvin-Voigt model, we could automatically enforce:

(a<1) & (a>0)
(β<1) & (β>0)
-a+β < 0

without having to type anything out explicitly! The first two should come from regular upper and lower bounds, and up till now were doable manually in RHEOS, the last one could be done automatically using an explicit constraint relationship as shown in the NLopt readme example!

[moustachio-belvedere]

This is quite involved on two fronts:

1) Finding the most sensible way of programatically embedding model constraints into the optimisation object (may need changes in model struct constraints field)

2) Our current local gradient-free algorithm does not accept constraints, only bounds. NLopt docs state that COBYLA is the only local gradient-free algorithm that supports constraints. If COBYLA is not as effective as SUBPLEX then we would need to consider other optimisation packages. Supplying gradients is probably not an option for the more complicated fractional models and user supplied models.

[moustachio-belvedere]

If NLopt -> COBYLA is a reasonable direction, the simplest way to address point 1 is to embed our constraints as functions rather than expressions. If we do this we should be able to simply loop through constraints and add them to the optimiser object before beginning optimisation. I've made a small example script to demonstrate an applied constraint, would be easy to extend to incorporate the above if we want to:

https://github.com/JuliaRheology/RheoBenchAndTest/blob/master/20200607_constraintsinvestigate.jl

For NLopt to understand them, the constraints would have to be in terms of the parameter index of the fitting parameters vector. If we want to keep the constraints written in their functions in terms of actual keyword parameters then we could probably use a thin wrapper function, need to think carefully about how this would update when parameters are 'frozen' though.

[moustachio-belvedere]

Update to this which I investigated but forgot to document some time ago, automatic differentiation may be the way forward to open up our options WRT optimisation algorithms. Its ability to handle the mittag peddler function is TBC.