GalactimOptim + AutoForwardDiff + IPNewton problem

Specific Domains Optimization (Mathematical)
1

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)
2

That was a bad error message. Fixed: Better error message for constrained optimizers by ChrisRackauckas · Pull Request #198 · SciML/Optimization.jl · GitHub .

NLopt BFGS acts similarly, so I don’t think it’s a problem of Optim as much of BFGS here. I found that NLopt.LD_TNEWTON_PRECOND_RESTART works best here.

using NLopt
sol = solve(prob, NLopt.LD_LBFGS(), cb=mycallback, maxiters = 10000).minimum # 4.156774187867866
sol = solve(prob, NLopt.LD_TNEWTON_PRECOND_RESTART(), cb=mycallback, maxiters = 10000).minimum # 3.2253551599824806

Note that one thing you can do to help optimizations is help it exit correctly when the solver fails, i.e.

    if state.retcode != :Success
        return Inf
    end
3

Hi Chris. Thanks for the fix and advices. BTW, the display callback silently fails (there is no display), is the keyword different for NLopt solvers ? Another question: how do you get the final gradient ? Maybe using GalacticOptim is too high level for me (what would you advise as an alternative to it, still having easy access to AD) ?

S.

4

Open an issue so it shows up in my emails. It’s like a one line fix to fix that but I just haven’t gotten to it and right now it’ll get lost in the v1.7 update stuff, but an email could probably put it in queue for tonight.

5

OK, where do I file the issue ? NLopt.jl ?

S.

6

GalacticOptim.jl

7

https://github.com/SciML/GalacticOptim.jl/issues/199

1 Like