# Help with my MINLP

Hello, sorry for this long post and problem… but if any help is trully appreciated.

I need help with my MINLP, since the results i’m getting is not what i was expecting… This is an electrical engineering problem, where i need to allocate harmonic filters in some buses of a radial distribution feeder (maybe someone from the area should know what this is).

Here’s a WE (sorry, not so Minimal, it might take a while to run it). I’ve tried with a small similar problem, and its working fine, but if this one, i’m having issues. It has a bunch of precalculated stuff, so if you want, you can just skip the the optimization part.

``````using JuMP, CouenneNL

# Some previous calculation, don't bother with these. Go to line 94.

Nt = 17;
hmax = 25;
de = [0;1;2;3;4;5;6;7;8;3;2;11;12;12;14;14;16];
para = [1;2;3;4;5;6;7;8;9;10;11;12;13;14;15;16;17];
kj = 1;
fq = 100;
Qmod = 100000.0;
u = [13100.5; 12870.6; 12853.2; 12820.3; 12811.1; 12752.1; 12749.1; 12735.9; 12694.2; 12844.7; 12754.0; 12703.1; 12676.4; 12547.0; 12543.5; 12483.8; 12478.2];
DHTv3 = [0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0];
DHTvI = [0.240953; 4.67092; 4.77751; 4.92099; 4.93973; 5.04262; 5.06766; 5.2084; 5.20804; 4.77749; 5.80687; 6.31057; 6.31013; 7.50116; 7.9857; 8.37866; 8.37864];
R = [0.010485; 0.0601575; 0.0214154; 0.0416004; 0.0110106; 0.09254; 0.0100485; 0.122281; 0.361213; 0.12022; 0.349655; 0.203799; 0.952345; 0.655657; 0.349655; 0.573609; 0.201135];
X = [0.0721325; 1.30207; 0.0598253; 0.116079; 0.0307316; 0.258411; 0.0280675; 0.150719; 0.270953; 0.155667; 0.452749; 0.263874; 1.23291; 0.848981; 0.452749; 0.706895; 0.247775];
Qshunt = [0.0; 0.0; 600000.0; 600000.0; 600000.0; 1.8e6; 0.0; 600000.0; 0.0; 0.0; 600000.0; 600000.0; 0.0; 0.0; 600000.0; 600000.0; 0.0];
PLlin = [0.0; 0.0; 200000.0; 400000.0; 1.5e6; 0.0; 800000.0; 200000.0; 1.0e6; 500000.0; 1.0e6; 300000.0; 200000.0; 800000.0; 500000.0; 1.0e6; 200000.0];
QLlin = [0.0; 0.0; 120000.0; 250000.0; 930000.0; 0.0; 500000.0; 120000.0; 620000.0; 310000.0; 620000.0; 190000.0; 120000.0; 500000.0; 310000.0; 620000.0; 120000.0];
PLnlin = [0.0; 0.0; 0.0; 0.0; 0.0; 3.0e6; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0];
QLnlin = [0.0; 0.0; 0.0; 0.0; 0.0; 2.26e6; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0];
I1_fase = [0.159969; -2.81295; 1.08462; 0.877378; 0.712251; 0.471352; 0.251137; 0.546692; 0.578663; 1.2784; 1.1575; 0.708211; 0.516676; 0.227036; -0.157791; -0.23005; -0.405142]
DHTv3_max = 5;
DHTvI_max = 6;
y1 = zeros(Complex128,Nt,hmax);
y2 = zeros(Complex128,Nt,hmax);
y3 = zeros(Complex128,Nt,hmax);
i_harm_1 = zeros(Complex128,Nt,1);

for h = 3:2:hmax
for i = 1:Nt
y1[i,h] = 1/(R[i]+im*h*X[i]);
y2[i,h] = im*h*(Qshunt[i]/(u[i]^2));
y3[i,h] = (PLlin[i]/u[i]^2)-im*(QLlin[i]/(h*u[i]^2));
end
end

for i = 1:Nt
i_harm_1[i] = conj((PLnlin[i]+im*QLnlin[i])/u[i]);
end

perc_esp_harm = [100.0 0.0 0.0 0.0  20.0 0.0 14.3 0.0 0.0 0.0  9.10 0.0 7.70 0.0 0.0 0.0  5.9 0.0   5.3 0.0 0.0 0.0  4.3 0.0   4.0];
teta_esp_harm = [  0.0 0.0 0.0 0.0 -67.8 0.0 11.9 0.0 0.0 0.0 -7.13 0.0 68.6 0.0 0.0 0.0 46.5 0.0 116.4 0.0 0.0 0.0 87.5 0.0 159.3];
c = 1;
Ih = zeros(Complex128,Nt,hmax);
Ih_mod = zeros(Nt,hmax);
Ih_fase = zeros(Nt,hmax);

for h = 3:2:hmax
for i = 1:Nt
if PLnlin[i] != 0 || QLnlin[i] != 0
Ih_mod[i,h] = (perc_esp_harm[c,h]/100)*abs(i_harm_1[i]);
Ih_fase[i,h] = Ih_fase[i,h]*180/pi;
Ih[i,h] = Ih_mod[i,h]*cos(Ih_fase[i,h])+im*Ih_mod[i,h]*sin(Ih_fase[i,h]);
c = c+1;
end
end
c = 1;
end

Ih_real = real(Ih);
Ih_imag = imag(Ih);

Yh = zeros(Complex128,Nt,Nt,hmax);

for h = 3:2:hmax
for L = 1:Nt
for C = 1:Nt
for i = 1:Nt
if L == C && (de[i] == L || para[i] == L)
if i == L
Yh[L,C,h] = Yh[L,C,h]+y1[i,h]+y2[i,h]+y3[i,h];
end
if i != L
Yh[L,C,h] = Yh[L,C,h]+y1[i,h];
end
end
if L != C && (de[i] == L && para[i] == C)
Yh[L,C,h] = Yh[L,C,h]-y1[i,h];
Yh[C,L,h] = Yh[L,C,h];
end
end
end
end
end

Yh_real = real(Yh);
Yh_imag = imag(Yh);

# End of the previous calculations

# My problem starts here

barramin = 12;
barramax = 15;
hsmin = 3;
hsmax = 25;
w1 = 0;
w2 = 0;
w3 = 1;

Xcf = zeros(Nt,1);
for i = 1:Nt
if i >= barramin && i <= barramax
Xcf[i] = u[i]^2/(kj*Qmod);
end
end

m = Model(solver=CouenneNLSolver())

@variable(m, a[i = barramin:barramax], Bin)
@variable(m, hsmin <= hs[i=barramin:barramax] <= hsmax, Int)
@variable(m, hs_aux[i = barramin:barramax], Int) # Auxiliar variable to make hs an odd number
@variable(m, y5_real[i = barramin:barramax,h = 3:2:hmax])
@variable(m, y5_imag[i = barramin:barramax,h = 3:2:hmax])
@variable(m, Yh_filtro_real[L = 1:Nt,C = 1:Nt, h = 3:2:hmax])
@variable(m, Yh_filtro_imag[L = 1:Nt,C = 1:Nt, h = 3:2:hmax])
@variable(m, Vh_filtro_real[i = 1:Nt,h = 3:2:hmax])
@variable(m, Vh_filtro_imag[i = 1:Nt,h = 3:2:hmax])
@variable(m, Vh_filtro[i = 1:Nt,h = 3:2:hmax])
@variable(m, DHTv3_filtro[i = 1:Nt])
@variable(m, DHTvI_filtro[i = 1:Nt])

@objective(m, Min, sum(w1*DHTv3_filtro[i]+w2*DHTvI_filtro[i] for i = 1:Nt)+sum(w3*a[i] for i = barramin:barramax))

@constraint(m, con1[i = barramin:barramax], hs[i] == 2*hs_aux[i]+1) # Making hs an odd number
@NLconstraint(m , con2[i = barramin:barramax, h = 3:2:hmax], y5_real[i,h] == (Xcf[i]/(fq*hs[i]))/((Xcf[i]/(fq*hs[i]))^2+(h*Xcf[i]/(hs[i]^2)-Xcf[i]/h)^2))
@NLconstraint(m , con3[i = barramin:barramax, h = 3:2:hmax], y5_imag[i,h] == (-h*Xcf[i]/(hs[i]^2)+Xcf[i]/h)/((Xcf[i]/(fq*hs[i]))^2+(h*Xcf[i]/(hs[i]^2)-Xcf[i]/h)^2))

for h = 3:2:hmax, L = 1:Nt, C = 1:Nt
if L == C && (L >= barramin && L <= barramax)
@NLconstraint(m, Yh_filtro_real[L, C, h] == Yh_real[L,C,h]+a[L]*y5_real[L,h]);
@NLconstraint(m, Yh_filtro_imag[L, C, h] == Yh_imag[L,C,h]+a[L]*y5_imag[L,h]);
else
@constraint(m, Yh_filtro_real[L, C, h] == Yh_real[L,C,h]);
@constraint(m, Yh_filtro_imag[L, C, h] == Yh_imag[L,C,h]);
end
end

for h = 3:2:hmax, L = 1:Nt
@constraint(m, sum(Yh_filtro_real[L,C,h]*Vh_filtro_real[C,h]-Yh_filtro_imag[L,C,h]*Vh_filtro_imag[C,h] for C = 1:Nt) == Ih_real[L,h])
@constraint(m, sum(Yh_filtro_imag[L,C,h]*Vh_filtro_real[C,h]+Yh_filtro_real[L,C,h]*Vh_filtro_imag[C,h] for C = 1:Nt) == Ih_imag[L,h])
end

@NLconstraint(m, con4[i = 1:Nt, h = 3:2:hmax], Vh_filtro[i,h] == sqrt(1e-20+Vh_filtro_real[i,h]^2+Vh_filtro_imag[i,h]^2))

@NLconstraint(m, con5[i = 1:Nt], DHTv3_filtro[i] == sqrt(1e-20+sum(Vh_filtro[i,h]^2 for h = 3:2:hmax if rem(h,3) == 0))/u[i]*100)
@NLconstraint(m, con6[i = 1:Nt], DHTvI_filtro[i] == sqrt(1e-20+sum(Vh_filtro[i,h]^2 for h = 3:2:hmax if rem(h,3) != 0))/u[i]*100)
@constraint(m, con7[i = 1:Nt], DHTv3_filtro[i] <= DHTv3_max)
@constraint(m, con8[i = 1:Nt], DHTvI_filtro[i] <= DHTvI_max)

status = solve(m)

DHTv3_filtro_sol = getvalue(DHTv3_filtro);
DHTvI_filtro_sol = getvalue(DHTvI_filtro);
println("Objective value: ", getobjectivevalue(m))
println("a = ", getvalue(a))
println("hs = ", getvalue(hs))
println("DHTv3_filtro = ", getvalue(DHTv3_filtro))
println("DHTvI_filtro = ", getvalue(DHTvI_filtro))
``````

I’ve attached some pictures showing the model and explaning a bit more of my problem:

As result, i got this:

Wich doesn’t make sense to me, because i know theres a better answer than that. For example, if only `a[12] = 1` and `hs[12] = 5`, i have a solution to my problem (all those `DHTv3[i]` and `DHTvI[i]` are under the `5` and `6` limits).

Any help is really appreciated, even for programming tips, to make my code run faster (as you can see, its taking a lot of time to run).

If you know a better solution, you should use the `start` keyword in `@variable`. For example

``````@variable(model, x >= 0, start = 1.234)
``````

MINLP’s are very difficult to solve. So it isn’t surprising that it takes a long time to solve. To speed things up, you could change constraints 7 and 8 into variable bounds. You could also look at reducing the number of variables and constraints by using `@expression` and `@NLexpression` (some docs).

``````@variable(m, x)
@variable(m, y)
@NLconstraint(m, y == sin(x) + cos(x))
@objective(m, Max, y)
``````

go

``````@variable(m, x)
@NLexpression(m, y_expr, sin(x) + cos(x)
@NLobjective(m, Max, y_expr)
``````

Thank you! I’ll try those things!

I was trying to rewrite the following contraints/NLconstraints as expression/NLexpressions, since they are adding too many constraints to my problem.

``````   Nt = 17;
for h = 3:2:hmax, L = 1:Nt, C = 1:Nt
if L == C && (L >= barramin && L <= barramax)
@NLconstraint(m, Yh_filtro_real[L, C, h] == Yh_real[L,C,h]+a[L]*y5_real[L,h]);
@NLconstraint(m, Yh_filtro_imag[L, C, h] == Yh_imag[L,C,h]+a[L]*y5_imag[L,h]);
else
@constraint(m, Yh_filtro_real[L, C, h] == Yh_real[L,C,h]);
@constraint(m, Yh_filtro_imag[L, C, h] == Yh_imag[L,C,h]);
end
end
``````

I tried replacing for this:

``````  for h = 3:2:hmax, L = 1:Nt, C = 1:Nt
if L == C && (L >= barramin && L <= barramax)
@NLexpression(m, Yh_filtro_real[L, C, h], Yh_real[L,C,h]+a[L]*y5_real[L,h]);
@NLexpression(m, Yh_filtro_imag[L, C, h], Yh_imag[L,C,h]+a[L]*y5_imag[L,h]);
else
@expression(m, Yh_filtro_real[L, C, h], Yh_real[L,C,h]);
@expression(m, Yh_filtro_imag[L, C, h], Yh_imag[L,C,h]);
end
end
``````

But it didn’t worked, i got the following error:

`ERROR: LoadError: Failed attempt to index JuMPArray along dimension 1: 1 ∉ 17`

In my code `Nt = 17`.

Any tips in how can i replace those constraints for expressions?

Take a look at some of the syntax for expressions in the JuMP documentation: Expressions and Constraints — JuMP -- Julia for Mathematical Optimization 0.18 documentation

For example, you have something like the following, which won’t work:

``````model = Model()
@variable(model, x[1:2])
for i in 1:2
@expression(model, my_ex[i], x[i]^2)
end
``````

``````model = Model()
@variable(model, x[1:2])
@expression(model, my_ex[i in 1:2], x[i]^2)
``````

Another option is

``````model = Model()
@variable(model, x[1:2])
my_ex = Dict()
for i in 1:2
my_ex[i] = @expression(model, x[i]^2)
end
``````

So, i assume this option won’t work for me, since i got an `if` condition in my model:

``````model = Model()
@variable(model, x[1:2])
@expression(model, my_ex[i in 1:2], x[i]^2)
``````

I guess the second option would work, i just need to put the conditions inside that `for` loop… But what does `Dict()` means?

`Dict()` creates a dictionary.

For loops also work:

``````model = Model()
@variable(model, x[1:2])
my_ex = Dict()
for i in 1:2
if isodd(i)
my_ex[i] = @expression(model, x[i]^2)
else
my_ex[i] = @expression(model, x[i])
end
end
``````
1 Like

Hello again,

Say i got this integer variable:

`@variable(m, x_min <= x[i=1:3] <= x_max, Int)`

Is there anyway i can write a constraint that every `x[i]` should be different from the other?

There are a variety of formulations for all different. A simple one would be `1 <= abs(x[i] - x[j])` or `1 <= (x[i] - x[j])^2` for each pair of variables.

2 Likes