Getobjectivevalue does not return value of my objective function

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

0 Likes