Skip to content

Commit

Permalink
tests passing, closes #778, closes #829
Browse files Browse the repository at this point in the history
  • Loading branch information
@ccoffrin
ccoffrin committed Aug 7, 2022
1 parent a893bd3 commit 3757e31
Show file tree
Hide file tree
Showing 5 changed files with 231 additions and 106 deletions.
2 changes: 2 additions & 0 deletions src/core/base.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -218,6 +218,8 @@ function _check_missing_keys(dict, keys, type)
end
end

# enables support for v[1]
Base.getindex(v::JuMP.VariableRef, i::Int) = v

nw_ids(pm::AbstractPowerModel) = _IM.nw_ids(pm, pm_it_sym)
nws(pm::AbstractPowerModel) = _IM.nws(pm, pm_it_sym)
Expand Down
250 changes: 172 additions & 78 deletions src/core/objective.jl
View file
Original file line number Diff line number Diff line change
@@ -1,6 +1,3 @@
# enables support for v[1]
Base.getindex(v::JuMP.VariableRef, i::Int) = v

"""
Checks if any generator cost model will require a JuMP nonlinear expression
"""
Expand Down Expand Up @@ -39,7 +36,7 @@ function objective_min_fuel_and_flow_cost(pm::AbstractPowerModel; kwargs...)
nl_gen = check_nl_gen_cost_models(pm)
nl_dc = check_nl_dcline_cost_models(pm)

nl = nl_gen || nl_dc
nl = nl_gen || nl_dc || typeof(pm) <: AbstractIVRModel

expression_pg_cost(pm; nonlinear=nl, kwargs...)
expression_p_dc_cost(pm; nonlinear=nl, kwargs...)
Expand All @@ -52,10 +49,21 @@ function objective_min_fuel_and_flow_cost(pm::AbstractPowerModel; kwargs...)
for (n, nw_ref) in nws(pm))
)
else
pg_cost = Dict()
p_dc_cost = Dict()
for (n, nw_ref) in nws(pm)
for (i,gen) in nw_ref[:gen]
pg_cost[(n,i)] = var(pm, n, :pg_cost, i)
end
for (i,dcline) in nw_ref[:dcline]
p_dc_cost[(n,i)] = var(pm, n, :p_dc_cost, i)
end
end

return JuMP.@NLobjective(pm.model, Min,
sum(
sum( var(pm, n, :pg_cost, i) for (i,gen) in nw_ref[:gen]) +
sum( var(pm, n, :p_dc_cost, i) for (i,dcline) in nw_ref[:dcline])
sum( pg_cost[n,i] for (i,gen) in nw_ref[:gen]) +
sum( p_dc_cost[n,i] for (i,dcline) in nw_ref[:dcline])
for (n, nw_ref) in nws(pm))
)
end
Expand All @@ -64,20 +72,29 @@ end

""
function objective_min_fuel_cost(pm::AbstractPowerModel; kwargs...)
nl = check_nl_gen_cost_models(pm)
nl_gen = check_nl_gen_cost_models(pm)

nl = nl_gen || typeof(pm) <: AbstractIVRModel

expression_pg_cost(pm; nonlinear=nl, kwargs...)

if !nl
return JuMP.@objective(pm.model, Min,
sum(
sum( var(pm, n, :pg_cost, i) for (i,gen) in nw_ref[:gen])
sum( var(pm, n, :pg_cost, i) for (i,gen) in nw_ref[:gen])
for (n, nw_ref) in nws(pm))
)
else
pg_cost = Dict()
for (n, nw_ref) in nws(pm)
for (i,gen) in nw_ref[:gen]
pg_cost[(n,i)] = var(pm, n, :pg_cost, i)
end
end

return JuMP.@NLobjective(pm.model, Min,
sum(
sum( var(pm, n, :pg_cost, i) for (i,gen) in nw_ref[:gen])
sum( pg_cost[n,i] for (i,gen) in nw_ref[:gen])
for (n, nw_ref) in nws(pm))
)
end
Expand All @@ -88,7 +105,7 @@ end
cleans up raw pwl cost points in preparation for building a mathamatical model.
The key mathematical properties,
- the first and last points are strickly outside of the pmin-to-pmax range
- the first and last points are strictly outside of the pmin-to-pmax range
- pmin and pmax occur in the first and last line segments.
"""
function calc_pwl_points(ncost::Int, cost::Vector{<:Real}, pmin::Real, pmax::Real; tolerance=1e-2)
Expand Down Expand Up @@ -179,46 +196,28 @@ function expression_pg_cost(pm::AbstractPowerModel; nonlinear::Bool=false, repor
pg_cost = var(pm, n)[:pg_cost] = Dict{Int,Any}()

for (i,gen) in ref(pm, n, :gen)
if isa(var(pm, n, :pg, i), Array)
pg_terms = [var(pm, n, :pg, i)[c] for c in conductor_ids(pm, n)]
else
pg_terms = [var(pm, n, :pg, i)]
end

if gen["model"] == 1
pg_vars = [var(pm, n, :pg, i)[c] for c in conductor_ids(pm, n)]
pmin = sum(JuMP.lower_bound.(pg_vars))
pmax = sum(JuMP.upper_bound.(pg_vars))
if isa(pg_terms, Array{JuMP.VariableRef})
pmin = sum(JuMP.lower_bound.(pg_terms))
pmax = sum(JuMP.upper_bound.(pg_terms))
else
pmin = sum(gen["pmin"][c] for c in conductor_ids(pm, n))
pmax = sum(gen["pmax"][c] for c in conductor_ids(pm, n))
end

# note pmin/pmax may be different from gen["pmin"]/gen["pmax"] in the on/off case
points = calc_pwl_points(gen["ncost"], gen["cost"], pmin, pmax)

pg_cost_lambda = JuMP.@variable(pm.model,
[i in 1:length(points)], base_name="$(n)_pg_cost_lambda",
lower_bound = 0.0,
upper_bound = 1.0
)
JuMP.@constraint(pm.model, sum(pg_cost_lambda) == 1.0)

pg_expr = 0.0
pg_cost_expr = 0.0
for (i,point) in enumerate(points)
pg_expr += point.mw*pg_cost_lambda[i]
pg_cost_expr += point.cost*pg_cost_lambda[i]
end
JuMP.@constraint(pm.model, pg_expr == sum(pg_vars))
pg_cost[i] = pg_cost_expr
pg_cost[i] = _pwl_cost_expression(pm, pg_terms, points, nw=n, id=i, var_name="pg")

elseif gen["model"] == 2
pg = sum( var(pm, n, :pg, i)[c] for c in conductor_ids(pm, n))

cost_rev = reverse(gen["cost"])
if length(cost_rev) == 0
pg_cost[i] = 0.0
elseif length(cost_rev) == 1
pg_cost[i] = cost_rev[1]
elseif length(cost_rev) == 2
pg_cost[i] = cost_rev[1] + cost_rev[2]*pg
elseif length(cost_rev) == 3
pg_cost[i] = cost_rev[1] + cost_rev[2]*pg + cost_rev[3]*pg^2
else # length(cost_rev) >= 4
cost_rev_nl = cost_rev[4:end]
pg_cost[i] = JuMP.@NLexpression(pm.model, cost_rev[1] + cost_rev[2]*pg + cost_rev[3]*pg^2 + sum( v*pg^(d+3) for (d,v) in enumerate(cost_rev_nl)) )
end

pg_cost[i] = _polynomial_cost_expression(pm, pg_terms, cost_rev, nw=n, id=i, var_name="pg")
else
Memento.error(_LOGGER, "Only cost models of types 1 and 2 are supported at this time, given cost model type of $(model) on generator $(i)")
end
Expand All @@ -237,47 +236,28 @@ function expression_p_dc_cost(pm::AbstractPowerModel; nonlinear::Bool=false, rep
for (i,dcline) in ref(pm, n, :dcline)
arc = (i, dcline["f_bus"], dcline["t_bus"])

if isa(var(pm, n, :pg, i), Array)
p_dc_terms = [var(pm, n, :p_dc, arc)[c] for c in conductor_ids(pm, n)]
else
p_dc_terms = [var(pm, n, :p_dc, arc)]
end

if dcline["model"] == 1
p_dc_vars = [var(pm, n, :p_dc, arc)[c] for c in conductor_ids(pm, n)]
pmin = sum(JuMP.lower_bound.(p_dc_vars))
pmax = sum(JuMP.upper_bound.(p_dc_vars))
if isa(p_dc_terms, Array{JuMP.VariableRef})
pmin = sum(JuMP.lower_bound.(p_dc_terms))
pmax = sum(JuMP.upper_bound.(p_dc_terms))
else
pmin = sum(dcline["pminf"][c] for c in conductor_ids(pm, n))
pmax = sum(dcline["pmaxf"][c] for c in conductor_ids(pm, n))
end

# note pmin/pmax may be different from dcline["pminf"]/dcline["pmaxf"] in the on/off case
points = calc_pwl_points(dcline["ncost"], dcline["cost"], pmin, pmax)

dc_p_cost_lambda = JuMP.@variable(pm.model,
[i in 1:length(points)], base_name="$(n)_dc_p_cost_lambda",
lower_bound = 0.0,
upper_bound = 1.0
)
JuMP.@constraint(pm.model, sum(dc_p_cost_lambda) == 1.0)

dc_p_expr = 0.0
dc_p_cost_expr = 0.0
for (i,point) in enumerate(points)
dc_p_expr += point.mw*dc_p_cost_lambda[i]
dc_p_cost_expr += point.cost*dc_p_cost_lambda[i]
end

JuMP.@constraint(pm.model, dc_p_expr == sum(p_dc_vars))
p_dc_cost[i] = dc_p_cost_expr
p_dc_cost[i] = _pwl_cost_expression(pm, p_dc_terms, points, nw=n, id=i, var_name="dc_p")

elseif dcline["model"] == 2
p_dc = sum( var(pm, n, :p_dc, arc) for c in conductor_ids(pm, n))

cost_rev = reverse(dcline["cost"])
if length(cost_rev) == 0
p_dc_cost[i] = 0.0
elseif length(cost_rev) == 1
p_dc_cost[i] = cost_rev[1]
elseif length(cost_rev) == 2
p_dc_cost[i] = cost_rev[1] + cost_rev[2]*p_dc
elseif length(cost_rev) == 3
p_dc_cost[i] = cost_rev[1] + cost_rev[2]*p_dc + cost_rev[3]*p_dc^2
else # length(cost_rev) >= 4
cost_rev_nl = cost_rev[4:end]
p_dc_cost[i] = JuMP.@NLexpression(pm.model, cost_rev[1] + cost_rev[2]*p_dc + cost_rev[3]*p_dc^2 + sum( v*p_dc^(d+3) for (d,v) in enumerate(cost_rev_nl)) )
end
p_dc_cost[i] = _polynomial_cost_expression(pm, p_dc_terms, cost_rev, nw=n, id=i, var_name="dc_p")
else
Memento.error(_LOGGER, "only cost models of types 1 and 2 are supported at this time, given cost model type of $(model) on dcline $(i)")
end
Expand All @@ -288,6 +268,120 @@ function expression_p_dc_cost(pm::AbstractPowerModel; nonlinear::Bool=false, rep
end


function _pwl_cost_expression(pm::AbstractPowerModel, x_list::Array{JuMP.VariableRef}, points; nw=0, id=1, var_name="x")
cost_lambda = JuMP.@variable(pm.model,
[i in 1:length(points)], base_name="$(nw)_$(var_name)_cost_lambda_$(id)",
lower_bound = 0.0,
upper_bound = 1.0
)
JuMP.@constraint(pm.model, sum(cost_lambda) == 1.0)

expr = 0.0
cost_expr = 0.0
for (i,point) in enumerate(points)
expr += point.mw*cost_lambda[i]
cost_expr += point.cost*cost_lambda[i]
end
JuMP.@constraint(pm.model, expr == sum(x_list))

return cost_expr
end

function _pwl_cost_expression(pm::AbstractPowerModel, x_list, points; nw=0, id=1, var_name="x")
cost_lambda = JuMP.@variable(pm.model,
[i in 1:length(points)], base_name="$(nw)_$(var_name)_cost_lambda_$(id)",
lower_bound = 0.0,
upper_bound = 1.0
)
JuMP.@constraint(pm.model, sum(cost_lambda) == 1.0)

expr = 0.0
cost_expr = 0.0
for (i,point) in enumerate(points)
expr += point.mw*cost_lambda[i]
cost_expr += point.cost*cost_lambda[i]
end
JuMP.@NLconstraint(pm.model, expr == sum(x for x in x_list))

return cost_expr
end



# note that `cost_terms` should be providing in ascending order (the reverse of the Matpower spec.)
function _polynomial_cost_expression(pm::AbstractPowerModel, x_list::Array{JuMP.VariableRef}, cost_terms; nw=0, id=1, var_name="x")
x = sum(x_list)
if length(cost_terms) == 0
return 0.0
elseif length(cost_terms) == 1
return cost_terms[1]
elseif length(cost_terms) == 2
return cost_terms[1] + cost_terms[2]*x
elseif length(cost_terms) == 3
return cost_terms[1] + cost_terms[2]*x + cost_terms[3]*x^2
else # length(cost_terms) >= 4
cost_nl = cost_terms[4:end]
return JuMP.@NLexpression(pm.model, cost_terms[1] + cost_terms[2]*x + cost_terms[3]*x^2 + sum( v*x^(d+3) for (d,v) in enumerate(cost_nl)) )
end
end

# note that `cost_terms` should be providing in ascending order (the reverse of the Matpower spec.)
function _polynomial_cost_expression(pm::AbstractConicModels, x_list::Array{JuMP.VariableRef}, cost_terms; nw=0, id=1, var_name="x")
x = sum(x_list)
if length(cost_terms) == 0
return 0.0
elseif length(cost_terms) == 1
return cost_terms[1]
elseif length(cost_terms) == 2
return cost_terms[1] + cost_terms[2]*x
elseif length(cost_terms) == 3
x_lb = sum(JuMP.lower_bound.(x_list))
x_ub = sum(JuMP.upper_bound.(x_list))

x_sqr_lb = 0.0
x_sqr_ub = max(x_lb^2, x_ub^2)
if x_lb > 0.0
x_sqr_lb = x_lb^2
end
if x_ub < 0.0
x_sqr_lb = x_ub^2
end

x_sqr = JuMP.@variable(pm.model,
base_name="$(nw)_$(var_name)_sqr_$(id)",
lower_bound = x_sqr_lb,
upper_bound = x_sqr_ub,
start = 0.0
)
JuMP.@constraint(pm.model, [0.5, x_sqr, x] in JuMP.RotatedSecondOrderCone())

return cost_terms[1] + cost_terms[2]*x + cost_terms[3]*x_sqr
else # length(cost_terms) >= 4
Memento.error(_LOGGER, "the network cost data features a polynomial cost function that is not compatible with conic mathematical programs.")
end
end

# note that `cost_terms` should be providing in ascending order (the reverse of the Matpower spec.)
function _polynomial_cost_expression(pm::AbstractPowerModel, x_list, cost_terms; nw=0, id=1, var_name="x")
x = JuMP.@NLexpression(pm.model, sum(x for x in x_list))
if length(cost_terms) == 0
return 0.0
elseif length(cost_terms) == 1
return cost_terms[1]
elseif length(cost_terms) == 2
return JuMP.@NLexpression(pm.model, cost_terms[1] + cost_terms[2]*x)
elseif length(cost_terms) == 3
return JuMP.@NLexpression(pm.model, cost_terms[1] + cost_terms[2]*x + cost_terms[3]*x^2)
else # length(cost_terms) >= 4
cost_nl = cost_terms[4:end]
return JuMP.@NLexpression(pm.model, cost_terms[1] + cost_terms[2]*x + cost_terms[3]*x^2 + sum( v*x^(d+3) for (d,v) in enumerate(cost_nl)) )
end
end






function objective_max_loadability(pm::AbstractPowerModel)
nws = nw_ids(pm)
Expand Down
32 changes: 31 additions & 1 deletion src/util/obbt.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -406,7 +406,37 @@ end


function _constraint_obj_bound(pm::AbstractPowerModel, bound)
model = PowerModels.check_cost_models(pm)

model = nothing
for (n, nw_ref) in nws(pm)
for (i,gen) in nw_ref[:gen]
if haskey(gen, "cost")
if model == nothing
model = gen["model"]
else
if gen["model"] != model
Memento.error(_LOGGER, "cost models are inconsistent, the typical model is $(model) however model $(gen["model"]) is given on generator $(i)")
end
end
else
Memento.error(_LOGGER, "no cost given for generator $(i)")
end
end
for (i,dcline) in nw_ref[:dcline]
if haskey(dcline, "model")
if model == nothing
model = dcline["model"]
else
if dcline["model"] != model
Memento.error(_LOGGER, "cost models are inconsistent, the typical model is $(model) however model $(dcline["model"]) is given on dcline $(i)")
end
end
else
Memento.error(_LOGGER, "no cost given for dcline $(i)")
end
end
end

if model != 2
Memento.error(_LOGGER, "Only cost models of type 2 is supported at this time, given cost model type $(model)")
end
Expand Down
8 changes: 4 additions & 4 deletions test/opf-var.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -127,8 +127,8 @@ end
result = PowerModels._solve_opf_cl(data, SDPWRMPowerModel, sdp_solver)

@test result["termination_status"] == OPTIMAL
@test isapprox(result["objective"], 5728.62; atol = 1e0)
#@test isapprox(result["objective"], 5747.63; atol = 1e0)
#@test isapprox(result["objective"], 5728.62; atol = 1e0)
@test isapprox(result["objective"], 5747.63; atol = 1e0)
end
@testset "5-bus case" begin
data = build_current_data("../test/data/matpower/case5.m")
Expand All @@ -144,8 +144,8 @@ end
result = PowerModels._solve_opf_cl(data, SDPWRMPowerModel, sdp_solver)

@test result["termination_status"] == OPTIMAL || result["termination_status"] == ALMOST_OPTIMAL
@test isapprox(result["objective"], 7505.33; atol = 1e0)
#@test isapprox(result["objective"], 7637.95; atol = 1e0)
#@test isapprox(result["objective"], 7505.33; atol = 1e0)
@test isapprox(result["objective"], 7637.95; atol = 1e0)
end
end

Expand Down
Loading

0 comments on commit 3757e31

Please sign in to comment.