# Some JuMP variables values are out of bounds

**URL:** <https://discourse.julialang.org/t/some-jump-variables-values-are-out-of-bounds/26008>\
**Category:** Optimization (Mathematical)\
**Created:** [July 4, 2019, 11:32am UTC](https://discourse.julialang.org/t/some-jump-variables-values-are-out-of-bounds/26008 "2019-07-04T11:32:21Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![KestutisMa](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kestutisma/32/7835_2.png) [@KestutisMa](https://discourse.julialang.org/u/KestutisMa)\
**Post date:** [July 4, 2019, 11:32am UTC](https://discourse.julialang.org/t/some-jump-variables-values-are-out-of-bounds/26008/1 "2019-07-04T11:32:21Z")

</div>

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:

```julia
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:

```julia
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?

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [July 5, 2019, 1:38pm UTC](https://discourse.julialang.org/t/some-jump-variables-values-are-out-of-bounds/26008/2 "2019-07-05T13:38:11Z")

</div>

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](https://www.coin-or.org/Ipopt/documentation/node42.html).

You can pass parameters like so:

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

```

---

<div class="post-metadata">

**Author:** ![KestutisMa](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kestutisma/32/7835_2.png) [@KestutisMa](https://discourse.julialang.org/u/KestutisMa)\
**Post date:** [July 5, 2019, 5:28pm UTC](https://discourse.julialang.org/t/some-jump-variables-values-are-out-of-bounds/26008/3 "2019-07-05T17:28:35Z")

</div>

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

```julia
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):

```julia
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:

```nohighlight
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.

```

---

<div class="post-metadata">

**Author:** ![miles.lubin](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/miles.lubin/32/279_2.png) [@miles.lubin](https://discourse.julialang.org/u/miles.lubin)\
**Post date:** [July 5, 2019, 6:27pm UTC](https://discourse.julialang.org/t/some-jump-variables-values-are-out-of-bounds/26008/4 "2019-07-05T18:27:26Z")

</div>

See also Ipopt’s [`bound_relax_factor`](https://www.coin-or.org/Ipopt/documentation/node44.html#SECTION000114010000000000000) option. Zero is a valid setting, but it may inhibit convergence.

---

<div class="post-metadata">

**Author:** ![KestutisMa](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kestutisma/32/7835_2.png) [@KestutisMa](https://discourse.julialang.org/u/KestutisMa)\
**Post date:** [July 10, 2019, 4:32pm UTC](https://discourse.julialang.org/t/some-jump-variables-values-are-out-of-bounds/26008/5 "2019-07-10T16:32:53Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![KestutisMa](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kestutisma/32/7835_2.png) [@KestutisMa](https://discourse.julialang.org/u/KestutisMa)\
**Post date:** [November 4, 2019, 12:58pm UTC](https://discourse.julialang.org/t/some-jump-variables-values-are-out-of-bounds/26008/6 "2019-11-04T12:58:27Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [November 4, 2019, 5:22pm UTC](https://discourse.julialang.org/t/some-jump-variables-values-are-out-of-bounds/26008/7 "2019-11-04T17:22:47Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![KestutisMa](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kestutisma/32/7835_2.png) [@KestutisMa](https://discourse.julialang.org/u/KestutisMa)\
**Post date:** [November 6, 2019, 2:20pm UTC](https://discourse.julialang.org/t/some-jump-variables-values-are-out-of-bounds/26008/8 "2019-11-06T14:20:35Z")

</div>

Thanks, that worked.