Hi, I have problem while running optimization of our biophysical model (which is using Markov chain methodology): some parameter values during JuMP optimisation procedure are set to “0” by autodiff, but this is out of bounds, that causes NaN values of optimized function.
I use JuMP autodiff. Optimization variables values passed by ForwardDiff displayed using (x->x.value).(par) command.
Parameters are stored in 4x8 array. “Out ouf bounds” illustration:
Creating model..
Par lower bound:
4×8 Array{Float64,2}:
0.01 5.0 2.0e-8 1.0e-9 1.0e10 1.0e10 0.5 -1.0
0.01 5.0 2.0e-8 1.0e-9 1.0e10 1.0e10 0.5 1.0
0.01 5.0 2.0e-8 1.0e-15 1.0e10 1.0e10 0.5 -1.0
0.01 5.0 2.0e-8 1.0e-15 1.0e10 1.0e10 0.5 1.0
Par upper bound:
4×8 Array{Float64,2}:
0.3 40.0 5.0e-8 3.0e-9 1.0e10 1.0e10 0.5 -1.0
0.3 40.0 4.0e-8 3.0e-9 1.0e10 1.0e10 0.5 1.0
0.3 40.0 5.0e-8 1.0e-15 1.0e10 1.0e10 0.5 -1.0
0.3 40.0 4.0e-8 1.0e-15 1.0e10 1.0e10 0.5 1.0
Starting optimization..
This is Ipopt version 3.12.10, running with linear solver mumps.
NOTE: Other linear solvers might be more efficient (see Ipopt documentation).
Number of nonzeros in equality constraint Jacobian...: 0
Number of nonzeros in inequality constraint Jacobian.: 0
Number of nonzeros in Lagrangian Hessian.............: 0
par=
4×8 Array{Float64,2}:
0.0 0.0 0.0 0.0 1.0e10 1.0e10 0.5 -1.0
0.0 0.0 0.0 0.0 1.0e10 1.0e10 0.5 1.0
0.0 0.0 0.0 1.0e-15 1.0e10 1.0e10 0.5 -1.0
0.0 0.0 0.0 1.0e-15 1.0e10 1.0e10 0.5 1.0
...
MWE of our model optimization:
using JuMP, Ipopt, PyPlot, LinearAlgebra, DualNumbers, ForwardDiff
function MC36SM_Mindaugo_SS(vj::Float64, par::Array{T,2}, p::Array{T}, pc1c2::Float64, pc2c1::Float64) where {T <: Real}
PPP = Array{T}(undef, 36, 36); PPP .= 0;
p1 = 0;
p2 = 0;
p3 = 0;
p4 = 0;
A = par[:,1];
v0 = par[:,2];
Go = par[:,3];
Gc = par[:,4];
Ro = par[:,5];
Rc = par[:,6];
Pt = par[:,7];
P = par[:,8];
K = Pt[1];
g = Array{T}(undef, 36, 4); g .= 0;
gg = Array{T}(undef, 36); gg .= 0; #buffer for the sums
v = Array{T}(undef, 36, 4); v .= 0;
k = Array{T}(undef, 36, 4); k .= 0;
R = Array{T}(undef, 36, 4); R .= 0;
# //states conductances
for i1 = 1:2
for i2 = 1:3
for i3 = 1:3
for i4 = 1:2
if (i1 == 1)
g[(i1 - 1) * 18 + (i2 - 1) * 6 + (i3 - 1) * 2 + (i4 - 1) * 1 + 1,1] = Go[1];
R[(i1 - 1) * 18 + (i2 - 1) * 6 + (i3 - 1) * 2 + (i4 - 1) * 1 + 1,1] = Ro[1];
else
g[(i1 - 1) * 18 + (i2 - 1) * 6 + (i3 - 1) * 2 + (i4 - 1) * 1 + 1,1] = Gc[1];
R[(i1 - 1) * 18 + (i2 - 1) * 6 + (i3 - 1) * 2 + (i4 - 1) * 1 + 1,1] = Rc[1];
end
if (i2 == 1)
g[(i1 - 1) * 18 + (i2 - 1) * 6 + (i3 - 1) * 2 + (i4 - 1) * 1 + 1,2] = Go[2];
R[(i1 - 1) * 18 + (i2 - 1) * 6 + (i3 - 1) * 2 + (i4 - 1) * 1 + 1,2] = Ro[2];
else
g[(i1 - 1) * 18 + (i2 - 1) * 6 + (i3 - 1) * 2 + (i4 - 1) * 1 + 1,2] = Gc[2];
R[(i1 - 1) * 18 + (i2 - 1) * 6 + (i3 - 1) * 2 + (i4 - 1) * 1 + 1,2] = Rc[2];
end
if (i3 == 1)
g[(i1 - 1) * 18 + (i2 - 1) * 6 + (i3 - 1) * 2 + (i4 - 1) * 1 + 1,3] = Go[3];
R[(i1 - 1) * 18 + (i2 - 1) * 6 + (i3 - 1) * 2 + (i4 - 1) * 1 + 1,3] = Ro[3];
else
g[(i1 - 1) * 18 + (i2 - 1) * 6 + (i3 - 1) * 2 + (i4 - 1) * 1 + 1,3] = Gc[3];
R[(i1 - 1) * 18 + (i2 - 1) * 6 + (i3 - 1) * 2 + (i4 - 1) * 1 + 1,3] = Rc[3];
end
if (i4 == 1)
g[(i1 - 1) * 18 + (i2 - 1) * 6 + (i3 - 1) * 2 + (i4 - 1) * 1 + 1,4] = Go[4];
R[(i1 - 1) * 18 + (i2 - 1) * 6 + (i3 - 1) * 2 + (i4 - 1) * 1 + 1,4] = Ro[4];
else
g[(i1 - 1) * 18 + (i2 - 1) * 6 + (i3 - 1) * 2 + (i4 - 1) * 1 + 1,4] = Gc[4];
R[(i1 - 1) * 18 + (i2 - 1) * 6 + (i3 - 1) * 2 + (i4 - 1) * 1 + 1,4] = Rc[4];
end
end
end
end
end
if any(isnan.(g))
display(g)
# display(i)
global db = g
error("breakpoint (not error)")
end
for a = 1:4 #// voltages and rectification
# gg = 1 ./ sum( (1 ./ g), dims = 2 ); #sums rows
fill!(gg, 0.0)
row_sum!(gg, g) #happens twice so took it out into function
# v=vj*[gg gg gg gg]./g;
for j = 1:4, i = 1:36 #loop so no allocations
@inbounds v[i, j] = vj * gg[i] / g[i, j];
end
g = g .* exp.(v ./ R);
end
for i = 1:36
for j = 1:4
@inbounds k[i,j] = exp(A[j] * (v[i,j] * P[j] - v0[j]));
end
end
# //------------Matrica P-------------------------------------------------
for i1 = 1:2 for i2 = 1:3 for i3 = 1:3 for i4 = 1:2
for j1 = 1:2 for j2 = 1:3 for j3 = 1:3 for j4 = 1:2
i = (i1 - 1) * 18 + (i2 - 1) * 6 + (i3 - 1) * 2 + (i4 - 1) * 1 + 1;
j = (j1 - 1) * 18 + (j2 - 1) * 6 + (j3 - 1) * 2 + (j4 - 1) * 1 + 1;
# %-------be Pt---------------------
# % p1
if (i1 == 1)
if (j1 == 1)
p1 = (1 - K * k[i,1] / (1 + k[i,1]));
else
p1 = K * k[i,1] / (1 + k[i,1]);
end
else
if (j1 == 2)
p1 = (1 - K / (1 + k[i,1]));
else
p1 = K / (1 + k[i,1]);
end
end
# %------------------------------
# % p2
if (i2 == 1)
if (j2 == 1)
p2 = (1 - K * k[i,2] / (1 + k[i,2]));
end
if (j2 == 2)
p2 = K * k[i,2] / (1 + k[i,2]);
end
if (j2 == 3)
p2 = 0;
end
end
if (i2 == 2)
if (j2 == 1)
p2 = (K / (1 + k[i,2])) * (1 - pc1c2);
end
if (j2 == 2)
p2 = (1 - K / (1 + k[i,2])) * (1 - pc1c2);
end
if (j2 == 3)
p2 = pc1c2;
end
end
if (i2 == 3)
if (j2 == 1)
p2 = 0;
end
if (j2 == 2)
p2 = pc2c1;
end
if (j2 == 3)
p2 = 1 - pc2c1;
end
end
# %----------------------------------------
# % p3
if (i3 == 1)
if (j3 == 1)
p3 = (1 - K * k[i,3] / (1 + k[i,3]));
end
if (j3 == 2)
p3 = K * k[i,3] / (1 + k[i,3]);
end
if (j3 == 3)
p3 = 0;
end
end
if (i3 == 2)
if (j3 == 1)
p3 = (K / (1 + k[i,3])) * (1 - pc1c2);
end
if (j3 == 2)
p3 = (1 - K / (1 + k[i,3])) * (1 - pc1c2);
end
if (j3 == 3)
p3 = pc1c2;
end
end
if (i3 == 3)
if (j3 == 1)
p3 = 0;
end
if (j3 == 2)
p3 = pc2c1;
end
if (j3 == 3)
p3 = 1 - pc2c1;
end
end
# %----------------------------
# %p4
if (i4 == 1)
if (j4 == 1)
p4 = (1 - K * k[i,4] / (1 + k[i,4]));
else
p4 = K * k[i,4] / (1 + k[i,4]);
end
else
if (j4 == 2)
p4 = (1 - K / (1 + k[i,4]));
else
p4 = K / (1 + k[i,4]);
end
end
PPP[j,i] = p1 * p2 * p3 * p4;
end
end
end
end
end
end
end
end
d_g = 100000;
g_old = 100000;
g_final = 0;
i = 0;
tmp = similar(p);
# ggg = zeros(MVector{36,Float64});
# ggg = 1 ./ sum( (1 ./ g), dims=2 ); #%sums rows
fill!(gg, 0.0);
row_sum!(gg, g)
while (d_g / (g_old + g_final) > .001e-5)
i = i + 1;
# p=p*PPP;
mul!(tmp, PPP, p)
p, tmp = tmp, p;
g_final = dot(gg, p); # (p*ggg)[1];
d_g = abs(g_old - g_final);
g_old = g_final;
# print(" ",i,", ", d_g / (g_old + g_final), " ")
end
# print(i, " ");
# println(i);
return g_final, i;
end
@inline function row_sum!(buffer, g) #fill the buffer in a loop -> no allocations
for j = 1:size(g, 2), i = 1:size(g, 1)
@inbounds buffer[i] += 1 / g[i, j];
end
buffer .= 1 ./ buffer;
end
# end
function getGj(N, vj, par::Array{T,2}) where {T <: Real}
println("getGJ()")
# global pp = par
# if readline() == "s"
# error("---breakpoint---")
# end
gj = Array{T}(undef, N); gj .= 0;
pc1c2 = 0;
pc2c1 = 0;
pc1c2 = 0.01;
pc2c1 = 0.001;
ppp = Array{T}(undef, 36); ppp .= 0;
ppp[1] = 1;
iter = 0;
for i = 1:N
gj[i], .. = MC36SM_Mindaugo_SS(vj[i], par, ppp, pc1c2, pc2c1);
end
println("gj=")
display(gj)
return gj
end
function my_model(N, vj, par_tuple)
println("my_model():")
par = reshape(collect(par_tuple), 4, 8);
y = getGj(N, vj, par)
if typeof(par) == Array{ForwardDiff.Dual{ForwardDiff.Tag{getfield(JuMP, Symbol("##85#87")){typeof(min_f)},Float64},Float64,11},2}
println("Forw.Diff detected, (x->x.value).(par) =")
display((x->x.value).(par))
# global pp = par
# error("aaaaaa")
plot(vj, (x->x.value).(y)) #using broadcast to get Dual.value of Dual tuple (because we can't do elementwise ops on tuple)
else
println("par=")
display(par);
plot(vj, y)
end
return y
end
function min_f(par_tuple...) # N , x, y <= args
println("minf()")
global ii += 1
println("iter=",ii)
s = sum(((my_model(N, vj_e, par_tuple) .- gj_e) ./ gj_e).^2) #TODO: sqrt(..), abs()
println("min_f value:", s)
return s
end
println("Starting.. ploting test function..")
global ii = 0; #iter counter
parAllFast = [ .15 20 20e-9 2e-9 1e10 1e10 .5 -1
.15 20 20e-9 2e-9 1e10 1e10 .5 1 ];
parAllSlow = [ .1 30 40e-9 1e-15 1e10 1e10 .5 -1
.1 30 40e-9 1e-15 1e10 1e10 .5 1 ];
par_test = [parAllFast; parAllSlow];
N = 40;
vj_e = collect(LinRange(0, 120, N));
gj_e = getGj(N, vj_e, par_test); # test gj
pygui(true)
plot(vj_e, gj_e, "b+--", markersize = 12)
autoscale(false)
println("Creating model..")
model = Model(with_optimizer(Ipopt.Optimizer))
register(model, :min_f, 32, min_f, autodiff = true) # vector argyments are not supported in user-defined functions, so using 1D tuple, length 32; par array(8x4) automaticaly converts to tuple(32) at "register(...)"
lb_fast_A = [0.01 5 20e-9 1e-9 1e10 1e10 0.5 -1];
lb_fast_B = [0.01 5 20e-9 1e-9 1e10 1e10 0.5 1];
lb_slow_A = [0.01 5 20e-9 1e-15 1e10 1e10 0.5 -1];
lb_slow_B = [0.01 5 20e-9 1e-15 1e10 1e10 0.5 1];
ub_fast_A = [ .30 40 50e-9 3e-9 1e10 1e10 0.5 -1];
ub_fast_B = [ .30 40 40e-9 3e-9 1e10 1e10 0.5 1];
ub_slow_A = [ .30 40 50e-9 1e-15 1e10 1e10 0.5 -1];
ub_slow_B = [ .30 40 40e-9 1e-15 1e10 1e10 0.5 1];
lb_AB = [lb_fast_A; lb_fast_B; lb_slow_A; lb_slow_B];
ub_AB = [ub_fast_A; ub_fast_B; ub_slow_A; ub_slow_B];
@variable(model, lb_AB[i,j] <= par_optim[i = 1:4, j = 1:8] <= ub_AB[i,j])
println("Par lower bound:")
display(lower_bound.(par_optim))
println("Par upper bound:")
display(upper_bound.(par_optim))
@NLobjective(model, Min, min_f(par_optim...))
println("Starting optimization..")
optimize!(model)
println("values = ")
display(reshape(value.(par_optim), 4, 8))
println("Total calls:")
display(ii)
Maybe somebody could suggest what cause such “out of bounds” behavior?