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];
teta_esp_harm_rad = teta_esp_harm*pi/180;
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] = teta_esp_harm_rad[c,h]+h*(I1_fase[i]-teta_esp_harm_rad[c,1]);
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).