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