using MatpowerCases
using JuMP
using Ipopt
include("acopf_input.jl")
path = pwd()"/case.dat"
casefilename, refbus, loadscale, thermalLimitscale, extras, mvaBase = readconfig(path)
generators, generatorlist, buses, lines = readcase(casefilename, extras, loadscale, thermalLimitscale, mvaBase)
function solve_acopf(thermalLimitscale, generators, buses, lines)
numbuses = length(buses)
numlines = length(lines)
numgenerators = length(generators)
theta_u=deg2rad(float(extras["theta_u"]))
m=Model(solver=IpoptSolver())
PV = find(map(b->b.kind==:PV, buses))
PQ = find(map(b->b.kind==:PQ, buses))
REF = find(map(b->b.kind==:Ref, buses))
@assert length(REF) == 1
@variable(m,ir[1:numbuses]) #real current injected at node
@variable(m,ij[1:numbuses]) #immaginary current injected at node
@variable(m,irij[1:numlines]) #real current flow from i to j
@variable(m,ijij[1:numlines]) #immaginary current flow from i to j
@variable(m,irji[1:numlines])# real current flow from j to i
@variable(m,ijji[1:numlines])#reactive current flow from j to i
@variable(m,theta[1:numbuses])#bus angle
@variable(m,theta_ij[1:numlines])
@variable(m, gp[1:numgenerators])#Active power generated at certain generator
@variable(m, gq[1:numgenerators])#Reactive power generated
@variable(m, v[1:numbuses] >=0)#Bus voltage
@variable(m, vr[1:numbuses] >=0)#Bus voltage real part
@variable(m, vj[1:numbuses])#Bus voltage imanginary part
#Objective Functions
@NLobjective(m, Min, sum{gp[i]^2*generators[i].pi1 + gp[i]*generators[i].pi2+generators[i].pi3, i=numgenerators})
#AC transmission constraint
irij_constr=[]
ijij_constr=[]
irji_constr=[]
ijji_constr=[]
for i=1:numlines
@constraint(m, theta_ij[i]==(theta[lines[i].head]-theta[lines[i].tail]))
@constraint(m,-6.28 <=theta_ij[i]<=6.28 )
#current leaving i to j
@NLconstraint(m,irij[i] == lines[i].γ * v[lines[i].head] * cos(theta[lines[i].head])-v[lines[i].head] *sin(theta[lines[i].head])*lines[i].β-lines[i].γ * v[lines[i].tail] *cos(theta[lines[i].tail])+lines[i].β * v[lines[i].tail] *sin(theta[lines[i].tail]))
@NLconstraint(m,ijij[i] == (lines[i].β-lines[i].b_charge/2) * v[lines[i].head] * cos(theta[lines[i].head])+v[lines[i].head] *sin(theta[lines[i].head])*lines[i].γ-lines[i].β * v[lines[i].tail] *cos(theta[lines[i].tail])-lines[i].γ * v[lines[i].tail] *sin(theta[lines[i].tail]))
#Current leaving j towards i
@NLconstraint(m,irji[i] == lines[i].γ * v[lines[i].tail] * cos(theta[lines[i].tail])-v[lines[i].tail] *sin(theta[lines[i].tail])*lines[i].β-lines[i].γ * v[lines[i].head] *cos(theta[lines[i].head])+lines[i].β * v[lines[i].head] *sin(theta[lines[i].head]))
@NLconstraint(m,ijji[i] == (lines[i].β-lines[i].b_charge/2) * v[lines[i].tail] * cos(theta[lines[i].tail])+v[lines[i].tail] *sin(theta[lines[i].tail])*lines[i].γ-lines[i].β * v[lines[i].head] *cos(theta[lines[i].head])-lines[i].γ * v[lines[i].head] *sin(theta[lines[i].head]))
#Thermal Limits
if lines[i].Imax > 0
@constraint(m, -lines[i].Imax<=irij[i] <= lines[i].Imax)
@constraint(m, -lines[i].Imax<=ijij[i] <= lines[i].Imax)
@NLconstraint(m,irij[i]^2+ijij[i]^2<= lines[i].Imax^2)
@NLconstraint(m,irji[i]^2+ijji[i]^2<= lines[i].Imax^2)
end
end
#Nodal Current Balancing & bus voltage limit &nodal power balanceing
vr_constr=[]
vj_constr=[]
ir_constr=[]
ij_constr=[]
for i=1:numbuses
@constraint(m,-theta_u <=theta[i]<= theta_u)
@NLconstraint(m,vr[i]==v[i]*cos(theta[i]))
@NLconstraint(m,vj[i]==v[i]*sin(theta[i]))
@constraint(m,-buses[i].Vmax <= vr[i] <= buses[i].Vmax)
@constraint(m,-buses[i].Vmax <= vj[i] <= buses[i].Vmax)
@NLconstraint(m,(buses[i].Vmin)^2 <= vr[i]^2+vj[i]^2 <= (buses[i].Vmax)^2)#@1 to linearize
if buses[i].kind == :Ref
@constraint(m,theta[i] == 0)
elseif buses[i].kind == :PV #why this process
end
@NLconstraint(m, ir[i]==-sum{irij[k], k in buses[i].outlist}+sum{irij[k], k in buses[i].inlist}+vr[i]*buses[i].Gs-vr[i]*buses[i].Bs)
@NLconstraint(m, ij[i]==-sum{ijij[k], k in buses[i].outlist}+sum{ijij[k], k in buses[i].inlist}+vj[i]*buses[i].Gs-vj[i]*buses[i].Bs)
@NLconstraint(m, sum{gp[k], k in buses[i].genids}+vr[i]*ir[i]+vj[i]*ij[i]==buses[i].Pd)
@NLconstraint(m, sum{gq[k], k in buses[i].genids}+vj[i]*ir[i]-vr[i]*ij[i]==buses[i].Qd)
end
for i=1:numgenerators
@constraint(m, 0<=gp[i]<=generators[i].Pgmax)
@constraint(m, -generators[i].Qgmax<=gq[i]<=generators[i].Qgmax)
end
solve(m)
return getvalue(gp),getvalue(gq),getObjectiveValue(m),getvalue(vr),getvalue(vj), getvalue(irij),getvalue(ijij),getvalue(theta)
end
(gp_opt, gq_opt,obj,vr_opt,vj_opt, irij_opt,ijij_opt,theta_opt)=solve_acopf(thermalLimitscale, generators, buses, lines);
println("********************************************************************")
println("The optimal cost is ", obj, gp_opt,vr_opt,"Current flow is ",irij_opt)
println("********************************************************************")
So after I run the code the solver status is optimal and everything is fine. Compare the result with test case, they are also very similar. But the only problem is that getobjectivevalue() does not return the value of my objective function “sum{gp[i]^2*generators[i].pi1 + gp[i]*generators[i].pi2+generators[i].pi3, i=numgenerators}” rather, it returns the value of sum{gp[i], i=numgenerators}. May I ask what is the reason for this?
Thank you!