Hi, I am trying to solve a parameter fitting problem in a stiff biological ODE. The problem solves with correct precision in Matlab (with fmincon/IP) using the complex-step technique to compute the gradient, but since I am migrating a lot of stuff to Julia this example is very important to me. I made a first port and used Flux but constraints are not handled and the choice of solvers is small. I adapted the code for GalacticOptim, but the Optim.IPNewton() solver does not seem to like by combination:

```
ERROR: LoadError: Use OptimizationFunction to pass the derivatives or automatically generate them with one of the autodiff backends
```

Optim.BFGS() likes it, but fails to converge correctly.

Really guys, this is very HARD to migrate stuff to Julia, even for motivated guys like me (Julia code below).

S.

```
function syst_AM2(state,theta,t)
Y_x1,Y_x2,mu1_max,mu2_max,Ks1,Ks2,KI,mu0,S00,S10,S20,coeff = theta;
X1,X2,S0,S1,S2,CH4 = state
r0=mu0*S0;
r1=mu1_max*S1*X1/(S1+Ks1);
r2=mu2_max*S2*X2/(S2+Ks2+(S2^2/KI));
[Y_x1*r1, Y_x2*r2, -r0, -r1+r0, (1-Y_x1)*r1-r2, (1-Y_x2)*r2]
end
function loss(theta,p)
state0 = [X10;X20;theta[9:11];0];
coeff = theta[12];
prob = ODEProblem(syst_AM2,state0,(t_CH4[1],t_CH4[end]),theta)
state = solve(prob, Rodas5(),saveat=t_CH4)
S2 = state[5,:]
CH4 = state[6,:]
deb_CH4 = CH4[2:end]-CH4[1:end-1]
e1 = CH4[1:end-1]-Exp_CH4*1
e2 = (deb_CH4-Exp_deb_CH4)*4
e3 = (S2[1:end-1]-coeff*Exp_AGV)*1
(sum(e1.^2)+sum(e2.^2)+sum(e3.^2))/length(t_CH4);
end
function mycallback(p, l)
display(loss(p))
return false
end
using Plots, DifferentialEquations, GalacticOptim, Optim
H=0.39;
R=0.195;
Surf=pi*R^2;
Exp_CH4 = [0,0.146695751,0.452562719,0.85245185,1.390508363,1.88574014,2.333176518,2.79329662,3.417909766,4.299654989,6.515952412,10.29002283,13.99900971,16.46767139,18.77271153,21.23159231,23.82660177,26.81067853,30.00135498,32.59923757,34.18509462,35.5158513,36.74890559,37.86754592,38.79872368,39.69511074,40.56481051,41.3865614,42.13182338,42.84323204,43.49045536];
t_CH4=0.0:size(Exp_CH4)[1]; # temps (jours)
Exp_AGV = [0.6307493,1.8081079,3.0002554,3.4371393,4.1710245,4.306228,4.4414315,4.576635,5.0225553,5.2529314,4.6852442,3.6354049,3.1327115,2.6300181,2.1273247,1.9948686,1.8574543,1.5407059,1.2099629,1.1136956,1.0174283,0.9211824,0.9121245,0.9195297,0.8298519,0.9236508,0.9041827,0.8846932,0.8652251,0.8372571,0.808326];
Exp_deb_CH4 = 0*Exp_CH4;
Exp_deb_CH4[2:end]= Exp_CH4[2:end]-Exp_CH4[1:end-1];
X10 = 0.283861214/(H*Surf);
X20 = 0.094620405/(H*Surf);
S00 = 5.180333444/(H*Surf);
S10 = 0.01/(H*Surf);
S20 = 0.029385961/(H*Surf);
theta_c0 = [0.25, 0.43, 0.5, 3.1, 21.0, 23.0, 4.0, 0.04, S00, S10, S20, 1.0]
cstr = [
0 1; #Y_X1
0 1; #Y_X2
0 Inf; #mu1max
0 Inf; #mu2max
1e-3 Inf; #Ks1
1e-3 Inf; #Ks2
1e-3 Inf; #KI
0 Inf; # mu0
0 Inf; # S00
0 Inf; # S10
0 Inf; # S20
0 Inf]; # AGV coefficient
fun = OptimizationFunction(loss, GalacticOptim.AutoForwardDiff())
prob = OptimizationProblem(fun, theta_c0, [], lb=cstr[:,1], ub=cstr[:,2])
sol = solve(prob, Optim.IPNewton(), cb=mycallback)
```