Getobjective does not return value of my objective function

jump
#1
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!

#2

Please read: PSA: make it easier to help you

What is the error? getObjectiveValue should probably be getobjectivevalue.

1 Like
#3

Thank you for your reply! I change it to getobjectivevalue but the result returned is still not the value of my objective function.

#4

There are two problems:

@NLobjective(m, Min, sum{gp[i]^2*generators[i].pi1 + gp[i]*generators[i].pi2+generators[i].pi3, i=numgenerators})

doesn’t do what you think it does.
i=numgenerators will only add the last generator. You need i=1:numgenerators.

It looks like you are using an old version of JuMP (indicated by sum{}). Update to the latest JuMP and use

@NLobjective(m, Min, sum(generators[i[.pi1 * gp[i]^2 + generators[i].pi2 * gp[i] + generators[i].pi3 for i in 1:numgenerators))
#5

I change it to i=1:numgenerators and it returns my objective function value this time. Thank you so much for your time!