hi there, I have a non-linear problem where the decision variable is a matrix of dimensions m * n * k. The constraints may also be somewhat complex. When I set m to a small subset of the dataset, e.g., 10, the program seems to work normally. However, when I set m to a slightly larger value, e.g., 50, the program does not produce any output, even though I set print_level to 12. There is no output at all. I am a bit confused and don’t know how to debug this. Is the program running normally (even though it hasn’t completed a single iteration in 12 hours)? How should I determine where the problem is? I am quite new to Julia and optimization problems. I am willing to provide more information. Any help would be greatly appreciated.
Hi @Chi_Jan,
Can you provide a (small) reproducible example?
You’re probably creating a very large problem by mistake, or you are creating deeply nested expression graphs that JuMP’s automatic differentiation problem is struggling with. Hard to say without an example.
To debug, try commenting out constraints until you have something that runs quickly. That should give you a hint about the problematic line.
Thank you for responding so quickly. I provide the simplified code here. Although there are many additional variables ignored in this gist, it shows all the constraints and objective functions in my program. Thanks again for your help.
I can’t run your code because of various errors.
You can improve your reproducible example by making it a single script that I can copy-paste into the REPL and run. You should also remove any unnecessary features (like using Parquet???).
As a comment, this is a very unique way of writing JuMP code. Once we have a reproducible example, I can suggest something that should work much better.
As some initial tips, Julia is not MATLAB. You’re welcome to use for-loops and sums. Vectorizing everything does not make the code faster.
I would write:
# function c1()
# m = dropdims(sum(decision_matrix,dims=1),dims=1)
# @constraint(model,0.99 .<=m .<=1.01)
# end
for r in 1:length(R), i in 1:length(I)
@constraint(model, 0.99 <= sum(decision_matrix[:, r, i]) <= 1.01)
end
The other issue, which is probably the cause of your slow-down, is that JuMP does not exploit common subexpressions. This is particularly problematic when you have a large summation in the denominator of a constraint fraction.
I would write:
# function c3()
# total = 0
# for f in F
# i = F_idx[f]
# total += sum(dropdims(sum(decision_matrix[i,:,:] .* pm1[:,:] .* pm2[i,:,:],dims=1),dims=1))
# end
# for f in F
# i = F_idx[f]
# n = sum(dropdims(sum(decision_matrix[f_idx,:,:] .* pm1[:,:] .* pm2[i,:,:],dims=1),dims=1))
# @constraint(model, lower[f]<= n / total <=upper[f])
# end
# end
@variable(model, x_pm1_pm2[F])
for f in F
@constraints(model, begin
x_pm1_pm2[f] == sum(decision_matrix[F_idx[f],:,:] .* (pm1 .* pm2[F_idx[f],:,:]))
lower[f] <= x_pm1_pm2[f] / sum(x_pm1_pm2) <= upper[f]
# better yet, assuming the denominator is >= 0
x_pm1_pm2[f] <= upper[f] * sum(x_pm1_pm2)
x_pm1_pm2[f] >= lower[f] * sum(x_pm1_pm2)
end)
end
In the better yet case, JuMP can recognize and exploit the fact that these constraints are linear.
Edit: Actually, since things are linear, you could even do:
@variable(model, x_pm1_pm2[F])
@expression(
model,
x_pm1_pm2[f in F],
sum(decision_matrix[F_idx[f],:,:] .* (pm1 .* pm2[F_idx[f],:,:])),
)
total = sum(x_pm1_pm2)
@constraints(model, begin
[f in F], x_pm1_pm2[f] <= upper[f] * total
[f in F], x_pm1_pm2[f] >= lower[f] * total
end)
1 Like
Hi, I have updated the code to ensure it is runnable. Would you mind taking some time to check it 
? The code in this example represents the general logic of my code work. The actual dataset used will be larger.
Can you post it here as a single code-block that I can copy-and-paste into the REPL?
hi, I put the code here for your information.
using JuMP,Ipopt
F = ["f1","f2","f3"]
R = ["r1","r2","r3"]
I = ["i1","i2","i3"]
model = Model(Ipopt.Optimizer)
decision_matrix = fill(0.1,length(F),length(R),length(I))
# Parameter matrix
pm1 = rand(length(R),length(I))
pm2 = rand(length(F),length(R),length(I))
pm3 = rand(length(F),length(R),length(I))
pm4 = rand(length(F),length(R),length(I))
println(size(pm1),size(pm2),size(pm3),size(pm4))
# Indexes for F,R,I
idx_arr = [i for i in 1:length(F)]
F_idx = Dict( first(row)=> getindex(row,2) for row in collect(zip(F,idx_arr)))
idx_arr = [i for i in 1:length(R)]
R_idx = Dict( first(row)=> getindex(row,2) for row in collect(zip(R,idx_arr)))
idx_arr = [i for i in 1:length(I)]
I_idx = Dict( first(row)=> getindex(row,2) for row in collect(zip(I,idx_arr)))
F_revese_idx = Dict([val => key for (key, val) in F_idx])
R_reverse_idx = Dict([val => key for (key, val) in R_idx])
I_reverse_idx = Dict([val => key for (key, val) in I_idx])
# bounds
c2_lower = Dict( f=> 0.5 * rand() for f in F)
c2_upper = Dict( f => 0.5 + 0.5 * rand() for f in F)
c3_lower = Dict( f=> 0.5 * rand() for f in F)
c3_upper = Dict( f => 0.5 + 0.5 * rand() for f in F)
c4_lower = Dict( (f,r,i)=> 0.5 * rand() for f in F for r in R for i in I )
c4_upper = Dict( (f,r,i) => 0.5 + 0.5 * rand() for f in F for r in R for i in I)
c5_lower = Dict( (f,r)=> 0.5 * rand() for f in F for r in R )
c5_upper = Dict( (f,r) => 0.5 + 0.5 * rand() for f in F for r in R )
target = Dict(f=> rand() for f in F)
# Constraints
function c1()
m = dropdims(sum(decision_matrix,dims=1),dims=1)
@constraint(model,0.99 .<=m .<=1.01)
end
function c2()
for f in F
f_idx = F_idx[f]
r1 = sum(dropdims(sum(decision_matrix[f_idx,:,:] .* pm1[:,:] .* pm2[f_idx,:,:],dims=1),dims=1)) # n*k .* n*k .* n *k => n*k
r2 = sum(dropdims(sum(decision_matrix[f_idx,:,:] .* pm1[:,:] .* pm2[f_idx,:,:] .* pm3[f_idx,:,:],dims=1),dims=1))
r3 = sum(dropdims(sum(decision_matrix[f_idx,:,:] .* pm1[:,:] .* pm2[f_idx,:,:] .* pm4[f_idx,:,:],dims=1),dims=1))
rr1 = r2 / r1
rr2 = r3 / r1
c2_expr = rr1 * (1-rr2) - rr2
@constraint(model, c2_lower[f]<=c2_expr<=c2_upper[f])
end
end
function c3()
total = 0
for f in F
f_idx = F_idx[f]
total += sum(dropdims(sum(decision_matrix[f_idx,:,:] .* pm1[:,:] .* pm2[f_idx,:,:],dims=1),dims=1))
end
for f in F
f_idx = F_idx[f]
n = sum(dropdims(sum(decision_matrix[f_idx,:,:] .* pm1[:,:] .* pm2[f_idx,:,:],dims=1),dims=1))
@constraint(model, c3_lower[f]<= n / total <=c3_upper[f])
end
end
function c4()
for f in F
for i in I
r1_total = 0
r2_total = 0
r3_total = 0
total = 0
for r in R
if occursin("r1",r)
r1_total += decision_matrix[F_idx[f],R_idx[r],I_idx[i]] * pm1[R_idx[r],I_idx[i]]
total += decision_matrix[F_idx[f],R_idx[r],I_idx[i]] * pm1[R_idx[r],I_idx[i]]
elseif occursin("r2",r)
r2_total += decision_matrix[F_idx[f],R_idx[r],I_idx[i]] * pm1[R_idx[r],I_idx[i]]
total += decision_matrix[F_idx[f],R_idx[r],I_idx[i]] * pm1[R_idx[r],I_idx[i]]
elseif occursin("r3",r)
r3_total += decision_matrix[F_idx[f],R_idx[r],I_idx[i]] * pm1[R_idx[r],I_idx[i]]
total += decision_matrix[F_idx[f],R_idx[r],I_idx[i]] * pm1[R_idx[r],I_idx[i]]
end
end
@constraint(model, get(c4_lower,(f,"r1",i),0)<=r1_total / total <= get(c4_upper,(f,"r1",i),1))
@constraint(model, get(c4_lower,(f,"r2",i),0)<=r2_total / total <=get(c4_upper,(f,"r2",i) ,1))
@constraint(model, get(c4_lower,(f,"r3",i),0)<=r3_total / total <=get(c4_upper,(f,"r3",i),1))
end
end
end
function c5()
for f in F
total = 0
r1_total = 0
r2_total = 0
r3_total = 0
for r in R
for i in I
if occursin("r1",r)
r1_total += decision_matrix[F_idx[f],R_idx[r],I_idx[i]] * pm1[R_idx[r],I_idx[i]]
total += decision_matrix[F_idx[f],R_idx[r],I_idx[i]] * pm1[R_idx[r],I_idx[i]]
elseif occursin("r2",r)
r2_total += decision_matrix[F_idx[f],R_idx[r],I_idx[i]] * pm1[R_idx[r],I_idx[i]]
total += decision_matrix[F_idx[f],R_idx[r],I_idx[i]] * pm1[R_idx[r],I_idx[i]]
elseif occursin("r3",r)
r3_total += decision_matrix[F_idx[f],R_idx[r],I_idx[i]] * pm1[R_idx[r],I_idx[i]]
total += decision_matrix[F_idx[f],R_idx[r],I_idx[i]] * pm1[R_idx[r],I_idx[i]]
end
end
end
@constraint(model, get(c5_lower,(f,"r1"),0)<=r1_total / total <= get(c5_upper,(f,"r1"),1))
@constraint(model, get(c5_lower,(f,"r2"),0)<=r2_total / total <=get(c5_upper,(f,"r2") ,1))
@constraint(model, get(c5_lower,(f,"r3"),0)<=r3_total / total <=get(c5_upper,(f,"r3"),1))
end
end
function obj_fn()
total = sum([ sum(decision_matrix[F_idx[f],:,:] .* pm1[:,:] .* pm2[F_idx[f],:,:]) for f in F])
err =0
for f in F
f_idx = F_idx[f]
r1 = sum(decision_matrix[f_idx,:,:] .* pm1[:,:] .* pm2[f_idx,:,:])
r2 = sum(decision_matrix[f_idx,:,:] .* pm1[:,:] .* pm2[f_idx,:,:] .* pm3[f_idx,:,:])
r3 = sum(decision_matrix[f_idx,:,:] .* pm1[:,:] .* pm2[f_idx,:,:] .* pm4[f_idx,:,:])
r3 = r3 / r1
r2 = r2 / r1
rr = r3 * (1.0 - r2 ) - r2
err += (r1/ total )*abs(rr - target[f])
end
return err
end
@variable(model,decision_matrix[1:length(F),1:length(R),1:length(I)],lower_bound=0.0,upper_bound=1.0)
obj = obj_fn()
c1()
c2()
c3()
c4()
c5()
@objective(model,Min,obj)
println(model)
optimize!(model)
I’d start by:
using JuMP, Ipopt
F = ["f1", "f2", "f3"]
R = ["r1", "r2", "r3"]
I = ["i1", "i2", "i3"]
model = Model(Ipopt.Optimizer)
@variable(model, 0 <= decision_matrix[F, R, I] <= 1)
# Parameter matrix
pm1 = rand(length(R), length(I))
pm2 = rand(length(F), length(R), length(I))
pm3 = rand(length(F), length(R), length(I))
pm4 = rand(length(F), length(R), length(I))
println(size(pm1), size(pm2), size(pm3), size(pm4))
# Indexes for F,R,I
F_idx = Dict(zip(F, 1:length(F)))
R_idx = Dict(zip(R, 1:length(R)))
I_idx = Dict(zip(I, 1:length(I)))
F_revese_idx = Dict(v => k for (k, v) in F_idx)
R_reverse_idx = Dict(v => k for (k, v) in R_idx)
I_reverse_idx = Dict(v => k for (k, v) in I_idx)
# bounds
c2_lower = Dict(f => 0.5 * rand() for f in F)
c2_upper = Dict(f => 0.5 + 0.5 * rand() for f in F)
c3_lower = Dict(f => 0.5 * rand() for f in F)
c3_upper = Dict(f => 0.5 + 0.5 * rand() for f in F)
c4_lower = Dict((f, r, i) => 0.5 * rand() for f in F for r in R for i in I)
c4_upper = Dict((f, r, i) => 0.5 + 0.5 * rand() for f in F for r in R for i in I)
c5_lower = Dict((f, r) => 0.5 * rand() for f in F for r in R )
c5_upper = Dict((f, r) => 0.5 + 0.5 * rand() for f in F for r in R )
target = Dict(f => rand() for f in F)
# Constraints
# function c1()
# m = dropdims(sum(decision_matrix,dims=1),dims=1)
# @constraint(model,0.99 .<=m .<=1.01)
# end
@constraint(model, [r in R, i in I], 0.99 <= sum(decision_matrix[:, r, i]) <= 1.01)
# function c3()
# total = 0
# for f in F
# f_idx = F_idx[f]
# total += sum(dropdims(sum(decision_matrix[f_idx,:,:] .* pm1[:,:] .* pm2[f_idx,:,:],dims=1),dims=1))
# end
# for f in F
# f_idx = F_idx[f]
# n = sum(dropdims(sum(decision_matrix[f_idx,:,:] .* pm1[:,:] .* pm2[f_idx,:,:],dims=1),dims=1))
# @constraint(model, c3_lower[f]<= n / total <=c3_upper[f])
# end
# end
@expression(model, x_pm1_pm2[f in F], sum(decision_matrix[f, :, :] .* (pm1 .* pm2[F_idx[f],:,:])))
@constraints(model, begin
[f in F], x_pm1_pm2[f] <= c3_upper[f] * sum(x_pm1_pm2)
[f in F], x_pm1_pm2[f] >= c3_lower[f] * sum(x_pm1_pm2)
end)
I’ll take a look at the rest tomorrow. Your main issue are terms like x[f] / sum(x). If the terms are linear, multiply out the denominator. If the terms are nonlinear, add a new decision variable and set y == sum(x) and use x[f] / y.
Thank u for your kindly help. According to your suggestion, I have changed my code to this. The program can now produce output.
# function c4()
# for f in F
# for i in I
# r1_total = 0
# r2_total = 0
# r3_total = 0
# total = 0
# for r in R
# if occursin("r1",r)
# r1_total += decision_matrix[F_idx[f],R_idx[r],I_idx[i]] * pm1[R_idx[r],I_idx[i]]
# total += decision_matrix[F_idx[f],R_idx[r],I_idx[i]] * pm1[R_idx[r],I_idx[i]]
# elseif occursin("r2",r)
# r2_total += decision_matrix[F_idx[f],R_idx[r],I_idx[i]] * pm1[R_idx[r],I_idx[i]]
# total += decision_matrix[F_idx[f],R_idx[r],I_idx[i]] * pm1[R_idx[r],I_idx[i]]
# elseif occursin("r3",r)
# r3_total += decision_matrix[F_idx[f],R_idx[r],I_idx[i]] * pm1[R_idx[r],I_idx[i]]
# total += decision_matrix[F_idx[f],R_idx[r],I_idx[i]] * pm1[R_idx[r],I_idx[i]]
# end
# end
# @constraint(model, get(c4_lower,(f,"r1",i),0)<=r1_total / total <= get(c4_upper,(f,"r1",i),1))
# @constraint(model, get(c4_lower,(f,"r2",i),0)<=r2_total / total <=get(c4_upper,(f,"r2",i) ,1))
# @constraint(model, get(c4_lower,(f,"r3",i),0)<=r3_total / total <=get(c4_upper,(f,"r3",i),1))
# end
# end
# end
function c4v2()
@variable(model,Z[F,["r1","r2","r3"],I])
for f in F
for i in I
@constraints(model,begin
Z[f,"r1",i] == sum( [ decision_matrix[F_idx[f],R_idx[r],I_idx[i]] * pm1[R_idx[r],I_idx[i]] for r in R if occursin("r1",r)] )
Z[f,"r2",i] == sum( [ decision_matrix[F_idx[f],R_idx[r],I_idx[i]] * pm1[R_idx[r],I_idx[i]] for r in R if occursin("r2",r)] )
Z[f,"r3",i] == sum( [ decision_matrix[F_idx[f],R_idx[r],I_idx[i]] * pm1[R_idx[r],I_idx[i]] for r in R if occursin("r3",r)] )
Z[f,"r1",i] <= sum( Z[f,:,i]) * get(c4_upper,(f,"r1",i),1)
Z[f,"r1",i] >= sum( Z[f,:,i]) * get(c4_lower,(f,"r1",i),1)
Z[f,"r2",i] <= sum( Z[f,:,i]) * get(c4_upper,(f,"r2",i),1)
Z[f,"r2",i] >= sum( Z[f,:,i]) * get(c4_lower,(f,"r2",i),1)
Z[f,"r3",i] <= sum( Z[f,:,i]) * get(c4_upper,(f,"r3",i),1)
Z[f,"r3",i] >= sum( Z[f,:,i]) * get(c4_lower,(f,"r3",i),1)
end)
end
end
end
# function c5()
# for f in F
# total = 0
# r1_total = 0
# r2_total = 0
# r3_total = 0
# for r in R
# for i in I
# if occursin("r1",r)
# r1_total += decision_matrix[F_idx[f],R_idx[r],I_idx[i]] * pm1[R_idx[r],I_idx[i]]
# total += decision_matrix[F_idx[f],R_idx[r],I_idx[i]] * pm1[R_idx[r],I_idx[i]]
# elseif occursin("r2",r)
# r2_total += decision_matrix[F_idx[f],R_idx[r],I_idx[i]] * pm1[R_idx[r],I_idx[i]]
# total += decision_matrix[F_idx[f],R_idx[r],I_idx[i]] * pm1[R_idx[r],I_idx[i]]
# elseif occursin("r3",r)
# r3_total += decision_matrix[F_idx[f],R_idx[r],I_idx[i]] * pm1[R_idx[r],I_idx[i]]
# total += decision_matrix[F_idx[f],R_idx[r],I_idx[i]] * pm1[R_idx[r],I_idx[i]]
# end
# end
# end
# @constraint(model, get(c5_lower,(f,"r1"),0)<=r1_total / total <= get(c5_upper,(f,"r1"),1))
# @constraint(model, get(c5_lower,(f,"r2"),0)<=r2_total / total <=get(c5_upper,(f,"r2") ,1))
# @constraint(model, get(c5_lower,(f,"r3"),0)<=r3_total / total <=get(c5_upper,(f,"r3"),1))
# end
# end
function c5v2()
@variable(model,Y[F,["r1","r2","r3"]])
for f in F
for r in R
if occursin("r1",r)
@constraints(model,begin
Y[f,"r1"] == sum([decision_matrix[F_idx[f],R_idx[r],I_idx[i]] * pm1[R_idx[r],I_idx[i]] for i in I ] )
Y[f,"r1"] <= sum(Y[f,:]) * get(c5_upper,(f,"r1"),1)
Y[f,"r1"] >= sum(Y[f,:]) * get(c5_lower,(f,"r1"),0)
end)
elseif occursin("r2",r)
@constraints(model,begin
Y[f,"r2"] == sum([decision_matrix[F_idx[f],R_idx[r],I_idx[i]] * pm1[R_idx[r],I_idx[i]] for i in I ] )
Y[f,"r2"] <= sum(Y[f,:]) * get(c5_upper,(f,"r1"),1)
Y[f,"r2"] >= sum(Y[f,:]) * get(c5_lower,(f,"r1"),0)
end)
elseif occursin("r3",r)
@constraints(model,begin
Y[f,"r2"] == sum([decision_matrix[F_idx[f],R_idx[r],I_idx[i]] * pm1[R_idx[r],I_idx[i]] for i in I ] )
Y[f,"r2"] <= sum(Y[f,:]) * get(c5_upper,(f,"r1"),1)
Y[f,"r2"] >= sum(Y[f,:]) * get(c5_lower,(f,"r1"),0)
end)
end
end
end
end
function obj_fnv2()
total = sum([ sum(decision_matrix[F_idx[f],:,:] .* pm1[:,:] .* pm2[F_idx[f],:,:]) for f in F])
err =0
for f in F
f_idx = F_idx[f]
r1 = sum(decision_matrix[f_idx,:,:] .* pm1[:,:] .* pm2[f_idx,:,:])
r2 = sum(decision_matrix[f_idx,:,:] .* pm1[:,:] .* pm2[f_idx,:,:] .* pm3[f_idx,:,:])
r3 = sum(decision_matrix[f_idx,:,:] .* pm1[:,:] .* pm2[f_idx,:,:] .* pm4[f_idx,:,:])
r3 = r3 / r1
r2 = r2 / r1
rr = r3 * (1.0 - r2 ) - r2
err += (r1/ total )*abs(rr - target[f])
end
return err
end
If possible, please allow me to ask a few silly questions.
- How can I ensure the consistency of two version of constraints. Can I verify it programmatically?
- Compared to Python with Pyomo, does Julia with JuMP have performance advantages for the same problems and constraints?
- Perhaps you have some suggestions for using Julia with JuMP to solve optimization problems, such as some dos or don’ts or some recipes.
1 Like
Your code still “looks” not very JuMP-like.
I didn’t test this, but I’d write your model like this. Hopefully it gives you some ideas. You want to re-use expressions where possible.
R_subset = ["r1", "r2", "r3"]
function rr_f(f)
r2 = x_pm1_pm2_pm3[f] / x_pm1_pm2[f]
r3 = x_pm1_pm2_pm4[f] / x_pm1_pm2[f]
return r3 * (1.0 - r2) - r2
end
model = Model(Ipopt.Optimizer)
@variables(model, begin
0 <= decision_matrix[F, R, I] <= 1
x_pm1_pm2[f in F]
x_pm1_pm2_pm3[f in F]
x_pm1_pm2_pm4[f in F]
Z[F, R_subset, I]
Y[F, R_subset]
end)
@constraints(model, begin
[f in F],
x_pm1_pm2[f] == sum(decision_matrix[f, :, :] .* (pm1 .* pm2[F_idx[f],:,:]))
[f in F],
x_pm1_pm2_pm3[f] == sum(decision_matrix[f, :, :] .* (pm1 .* pm2[F_idx[f],:,:] .* pm3[F_idx[f],:,:]))
[f in F],
x_pm1_pm2_pm4[f] == sum(decision_matrix[f, :, :] .* (pm1 .* pm2[F_idx[f],:,:] .* pm4[F_idx[f],:,:]))
# c1
[r in R, i in I], 0.99 <= sum(decision_matrix[:, r, i]) <= 1.01
# c2
[f in F], c2_lower[f] <= rr_f(f) <= c2_upper[f]
# c3
[f in F], x_pm1_pm2[f] <= c3_upper[f] * sum(x_pm1_pm2)
[f in F], x_pm1_pm2[f] >= c3_lower[f] * sum(x_pm1_pm2)
# c4
[f in F, r in R_subset, i in I],
Z[f, r, i] == sum(decision_matrix[f, rr, i] * pm1[rr, i] for rr in R if occursin(r, rr))
[f in F, r in R_subset, i in I],
Z[f, r, i] >= sum(Z[f, :, i]) * get(c4_lower, (f, r, i), 0)
[f in F, r in R_subset, i in I],
Z[f, r, i] <= sum(Z[f, :, i]) * get(c4_upper, (f, r, i), 1)
# c5
[f in F, r in R_subset],
Y[f, r] == sum(decision_matrix[f, rr, i] * pm1[rr, i] for i in I for rr in R if occursin(r, rr))
[f in F, r in R_subset],
Y[f, r] <= sum(Y[f, :]) * get(c5_upper, (f, r), 1)
[f in F, r in R_subset],
Y[f, r] >= sum(Y[f, :]) * get(c5_lower, (f, r), 0)
end)
@objective(
model,
Min,
sum(x_pm1_pm2[f] * abs(rr_f(f) - target[f]) for f in F) / sum(x_pm1_pm2),
)
optimize!(model)
For your questions:
- Nothing helpful. All you can really do is get a solution and verify that it meets your constraints. We can’t check if two JuMP models have the same feasibile set.
- For problems like this, not really. Pyomo uses AMPL’s solver library for derivatives, so for some models it will be faster and some it will be slower.
- There are two things I would say about your model:
i) there’s a math formulation problem. You have lots of normalizing denominators. Some solvers don’t like that. You’re repeating the construction of large expressions multiple times, so thedecisio_variable_pm1_pm2appears in c2 and the objective, etc. You’re also usingabsin the objective, which is non-smooth. Ipopt assumes functions are smooth and twice differentiable, so this can cause convergence issues.
ii) the coding style. I think this is just a “new to Julia” thing. Directly translating code that might work well in MATLAB or Python is not always the best thing to do in Julia. You also had a fewsum(dropdims(sum(...)))things, where only a single sum was needed, etc. I don’t have any helpful tips other than keep writing code, reading this forum, looking at the JuMP documentation, etc. There’s definitely a learning curve for writing Julia code, but this forum can be pretty helpful at pointing out different ways to do things.
please allow me to ask a few silly questions.
No question is silly. These threads are often quite helpful for lurking readers who have identical questions 
2 Likes