ACOPF Problem

I am trying to run an ACOPF (Alternating current optimal power flow) in the two-node network for a variable load in time “t”. The result of the model is just the zero values and nothing else. Can someone suggest to me a possible solution to it? Highly appreciated…

Code is attached as:

#Generator Limits
 
Pmin = [0.1,0.1];
Pmax = [100,100];
Qmin = [0.000001,0.000001];
Qmax = [200.0,200.0];

# Generation costs
c = [0.4,1.6];
# Line Limits
SLmax = 1.0;

# Voltage limits
Vmin = 0.95;
Vmax = 1.10;

end 
 
bb= Model(with_optimizer(Ipopt.Optimizer, print_level =1)) 
 
# Decision variables 
@variable(bb, p1[t = LoadData.Hour]); 
@variable(bb, p2[t = LoadData.Hour]);
@variable(bb, q1[t = LoadData.Hour]);
@variable(bb, q2[t = LoadData.Hour]);
@variable(bb, pl12[t = LoadData.Hour]);
@variable(bb, pl21[t = LoadData.Hour]);
@variable(bb, ql12[t = LoadData.Hour]);
@variable(bb, ql21[t = LoadData.Hour]);
@variable(bb, rc[i in Nodes, t in LoadData.Hour]);
@variable(bb, rd[i in Nodes, t in LoadData.Hour]);
@variable(bb, vi[i in Nodes, t in LoadData.Hour]);
@variable(bb, θi[i in Nodes, t in LoadData.Hour]); 
 
@constraint(bb, PL1[i = 1, j = 2,t in LoadData.Hour],p1[t] == pl12[t]) 
@constraint(bb, PL2[i = 1, j = 2,t in LoadData.Hour],p2[t] + LoadData.PV1[t]/100 -LoadData.Demand[t]/100== pl21[t])
@constraint(bb, QL1[i = 1, j = 2,t in LoadData.Hour],q1[t] == ql12[t]) 
@constraint(bb, QL2[i = 1, j = 2,t in LoadData.Hour],q2[t] == ql21[t])
@constraint(bb, PL1[i = 1, j = 2,t in LoadData.Hour],p1[t] == pl12[t])  
 
 
@NLconstraint(bb, pnn12[i in Nodes, j in Nodes, t in LoadData.Hour], pl12[t] == -(vi[1,t]^2)* Gbus[1,2] + vi[1,t]*vi[2,t]*(Gbus[1,2]*cos(θi[1,t]-θi[2,t])+ Bbus[1,2]*sin(θi[1,t]-θi[2,t])));
 
@NLconstraint(bb, pnn21[i in Nodes, j in Nodes, t in LoadData.Hour], pl21[t] == -(vi[2,t]^2)* Gbus[2,1] + vi[1,t]*vi[2,t]*(Gbus[2,1]*cos(θi[2,t]-θi[1,t])+ Bbus[2,1]*sin(θi[2,t]-θi[1,t])));

@NLconstraint(bb,qnn12[i in Nodes, j in Nodes, t in LoadData.Hour], ql12[t] == (vi[1,t]^2)* Bbus[1,2] + vi[1,t]*vi[2,t]*(Gbus[1,2]*sin(θi[1,t]-θi[2,t])- Bbus[1,2]*cos(θi[1,t]-θi[2,t]))); 

@NLconstraint(bb,qnn21[i in Nodes, j in Nodes, t in LoadData.Hour], ql21[t] == (vi[2,t]^2)* Bbus[2,1] + vi[1,t]*vi[2,t]*(Gbus[2,1]*sin(θi[2,t]-θi[1,t])- Bbus[2,1]*cos(θi[2,t]-θi[1,t])));  
 
 
#Technical generation limits 

@constraint(bb, glimit1[i in Nodes, t in LoadData.Hour],Pmin[1] <= p1[t] <= Pmax[1]);
@constraint(bb, glimit2[i in Nodes, t in LoadData.Hour],Pmin[1] <= p2[t] <= Pmax[1]);
@constraint(bb, qlimit1[i in Nodes, t in LoadData.Hour],Pmin[1] <= q1[t] <= Pmax[1]);
@constraint(bb, qlimit2[i in Nodes, t in LoadData.Hour],Pmin[1] <= q2[t] <= Pmax[1]); 
 
@constraint(bb, voltagelimit[i in Nodes, t in LoadData.Hour], vi[1,t] == 1.00); 
@constraint(bb, voltagelimit2[i in Nodes, t in LoadData.Hour], 0.95 <= vi[2,t]<= 1.05);
@constraint(bb, tetali[i in Nodes, t in LoadData.Hour], θi[1,t] == 0);
@constraint(bb, tetali2[i in Nodes, t in LoadData.Hour], -pi <= θi[2,t]<= pi);
 
 
@NLconstraint(bb, lc1[i in Nodes, t in LoadData.Hour], (pl12[t])^2+ (ql12[t])^2 <= (SLmax[1])^2);
@NLconstraint(bb, lc2[i in Nodes, t in LoadData.Hour], (pl21[t])^2+ (ql21[t])^2 <= (SLmax[1])^2); 
 
@objective(bb, Min, sum(c[1]*p1[t] + c[2]*p2[t] for t in LoadData.Hour))





From a quick glance at your code nothing stands out as a clear bug to me. There are some details that are unclear, for example if the Gbus and Bbus are defined correctly. I would recommend the usual mathematical programming debugging procedures to isolate and identify the problem.

If a reference implementation would he helpful in the debugging process you can have a look a this one here,

1 Like

Don’t add variable bounds via @constraint:

@variable(bb, p1[t = LoadData.Hour]);
# ...
@constraint(bb, glimit1[i in Nodes, t in LoadData.Hour],Pmin[1] <= p1[t] <= Pmax[1])

do instead:

@variable(bb, Pmin[1] <= p1[t = LoadData.Hour] <= Pmax[1]);

That might help quite a bit.