Help with my MINLP


#1

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:

pic2 pic3

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


#2

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

For example instead of

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

#3

Thank you! I’ll try those things!


#4

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?


#5

Take a look at some of the syntax for expressions in the JuMP documentation: http://www.juliaopt.org/JuMP.jl/0.18/refexpr.html#methods

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

Instead you need

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

#6

Thanks for the fast reply.

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?


#7

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