Some JuMP variables values are out of bounds

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?

All solvers use a tolerance for comparing bounds. When you write x >= 0.0, the computer interprets it as x >= -tolerance. By default, Ipopt uses a tolerance of 1e-8. So you should expect to observe small values like 1.0e-15.

Here are the Ipopt parameters if you want to change the tolerance: https://www.coin-or.org/Ipopt/documentation/node42.html.

You can pass parameters like so:

model = Model(with_optimizer(Ipopt.Optimizer, tol=1e-9))

I tried set tol values 1e-9 and 1e-16 but that didn’t change result.

It is strange that only first five objective function calls are with some parameters provided as zeros. I did detecting if parameters type are ForwardDiff Duals, and if they are converted to normal numbers (just for output):

4th call (many parameters are 0 (that is out of bounds), non zeros are only those for which I set same upper and lower bounds):

Forw.Diff detected, (x->x.value).(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

5th call was not from ForwadDiff (just two parameters are 0):

par=
4×8 Array{Float64,2}:
 0.0129  5.05  1.05e-8  0.0      1.0e10  1.0e10  0.5  -1.0
 0.0129  5.05  1.04e-8  0.0      1.0e10  1.0e10  0.5   1.0
 0.0129  5.05  1.05e-8  1.0e-15  1.0e10  1.0e10  0.5  -1.0
 0.0129  5.05  1.04e-8  1.0e-15  1.0e10  1.0e10  0.5   1.0

And terminated with:

Number of objective gradient evaluations             = 1
Number of equality constraint evaluations            = 0
Number of inequality constraint evaluations          = 0
Number of equality constraint Jacobian evaluations   = 0
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations             = 0
Total CPU secs in IPOPT (w/o function evaluations)   =      2.542
Total CPU secs in NLP function evaluations           =      0.196

EXIT: Invalid number in NLP function or derivative detected.

See also Ipopt’s bound_relax_factor option. Zero is a valid setting, but it may inhibit convergence.

Thank you, I decreased bound_relax_factor and error and invalid values are gone. It is strange that at the beginning first 4 calls to objective function parameters still are 0 (maybe some solver “warm-up” code?), but that doesn’t give error at end as before.