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.

Hi, I am still observing that JuMP passes optimization parameter value of 0 during first call to objective function, while bounds are @variable(model, 0.2 <= td <= 0.3). That causes NaN value of objective function and solver stops on message EXIT: Invalid number in NLP function or derivative detected.
Also, before JuMP stops, there was a single call to the objective function by ForwadDiff (probably for calculation of gradient), and in that case parameter (td) value was correctly withing the bounds.
What could be wrong?

You can try providing a starting value: @variable(model, 0.2 <= td <= 0.3, start = 0.2). See https://www.juliaopt.org/JuMP.jl/v0.20.0/variables/#Start-values-1

Thanks, that worked.