# 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
#     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.