Optimal Powerflow Nonlinear Optimization Benchmark
Data Load and Setup Code
This is generic setup code usable for all solver setups. Basically removing some unnecessary untyped dictionaries before getting to the benchmarks.
PRINT_LEVEL = 0
# This is a soft upper limit to the time of each optimization.
# If times go above this, they will halt early
MAX_CPU_TIME = 100.0
# Maximum number of variables in an optimization problem for the benchmark
# Anything with more variables is rejected and not run
# This is for testing. 100 is a good size for running a test of changes
# Should be set to typemax(Int) to run the whole benchmark
SIZE_LIMIT = 10001000import PowerModels
import ConcreteStructs
using BenchmarkTools
ConcreteStructs.@concrete struct DataRepresentation
data
ref
var_lookup
var_init
var_lb
var_ub
ref_gen_idxs
lookup_pg
lookup_qg
lookup_va
lookup_vm
lookup_lij
lookup_p_lij
lookup_q_lij
cost_arrs
f_bus
t_bus
ref_bus_idxs
ref_buses_idxs
ref_bus_gens
ref_bus_arcs
ref_branch_idxs
ref_arcs_from
ref_arcs_to
p_idxmap
q_idxmap
bus_pd
bus_qd
bus_gs
bus_bs
br_g
br_b
br_tr
br_ti
br_ttm
br_g_fr
br_b_fr
br_g_to
br_b_to
end
function load_and_setup_data(file_name)
data = PowerModels.parse_file(file_name)
PowerModels.standardize_cost_terms!(data, order=2)
PowerModels.calc_thermal_limits!(data)
ref = PowerModels.build_ref(data)[:it][:pm][:nw][0]
# Some data munging to type-stable forms
var_lookup = Dict{String,Int}()
var_init = Float64[]
var_lb = Float64[]
var_ub = Float64[]
var_idx = 1
for (i,bus) in ref[:bus]
push!(var_init, 0.0) #va
push!(var_lb, -Inf)
push!(var_ub, Inf)
var_lookup["va_$(i)"] = var_idx
var_idx += 1
push!(var_init, 1.0) #vm
push!(var_lb, bus["vmin"])
push!(var_ub, bus["vmax"])
var_lookup["vm_$(i)"] = var_idx
var_idx += 1
end
for (i,gen) in ref[:gen]
push!(var_init, 0.0) #pg
push!(var_lb, gen["pmin"])
push!(var_ub, gen["pmax"])
var_lookup["pg_$(i)"] = var_idx
var_idx += 1
push!(var_init, 0.0) #qg
push!(var_lb, gen["qmin"])
push!(var_ub, gen["qmax"])
var_lookup["qg_$(i)"] = var_idx
var_idx += 1
end
for (l,i,j) in ref[:arcs]
branch = ref[:branch][l]
push!(var_init, 0.0) #p
push!(var_lb, -branch["rate_a"])
push!(var_ub, branch["rate_a"])
var_lookup["p_$(l)_$(i)_$(j)"] = var_idx
var_idx += 1
push!(var_init, 0.0) #q
push!(var_lb, -branch["rate_a"])
push!(var_ub, branch["rate_a"])
var_lookup["q_$(l)_$(i)_$(j)"] = var_idx
var_idx += 1
end
@assert var_idx == length(var_init)+1
ref_gen_idxs = [i for i in keys(ref[:gen])]
lookup_pg = Dict{Int,Int}()
lookup_qg = Dict{Int,Int}()
lookup_va = Dict{Int,Int}()
lookup_vm = Dict{Int,Int}()
lookup_lij = Tuple{Int,Int,Int}[]
lookup_p_lij = Int[]
lookup_q_lij = Int[]
cost_arrs = Dict{Int,Vector{Float64}}()
for (i,gen) in ref[:gen]
lookup_pg[i] = var_lookup["pg_$(i)"]
lookup_qg[i] = var_lookup["qg_$(i)"]
cost_arrs[i] = gen["cost"]
end
for (i,bus) in ref[:bus]
lookup_va[i] = var_lookup["va_$(i)"]
lookup_vm[i] = var_lookup["vm_$(i)"]
end
for (l,i,j) in ref[:arcs]
push!(lookup_lij, (l,i,j))
push!(lookup_p_lij,var_lookup["p_$(l)_$(i)_$(j)"])
push!(lookup_q_lij,var_lookup["q_$(l)_$(i)_$(j)"])
end
f_bus = Dict{Int,Int}()
t_bus = Dict{Int,Int}()
for (l,branch) in ref[:branch]
f_bus[l] = branch["f_bus"]
t_bus[l] = branch["t_bus"]
end
ref_bus_idxs = [i for i in keys(ref[:bus])]
ref_buses_idxs = [i for i in keys(ref[:ref_buses])]
ref_bus_gens = ref[:bus_gens]
ref_bus_arcs = ref[:bus_arcs]
ref_branch_idxs = [i for i in keys(ref[:branch])]
ref_arcs_from = ref[:arcs_from]
ref_arcs_to = ref[:arcs_to]
p_idxmap = Dict(lookup_lij[i] => lookup_p_lij[i] for i in 1:length(lookup_lij))
q_idxmap = Dict(lookup_lij[i] => lookup_q_lij[i] for i in 1:length(lookup_lij))
bus_pd = Dict(i => 0.0 for (i,bus) in ref[:bus])
bus_qd = Dict(i => 0.0 for (i,bus) in ref[:bus])
bus_gs = Dict(i => 0.0 for (i,bus) in ref[:bus])
bus_bs = Dict(i => 0.0 for (i,bus) in ref[:bus])
for (i,bus) in ref[:bus]
if length(ref[:bus_loads][i]) > 0
bus_pd[i] = sum(ref[:load][l]["pd"] for l in ref[:bus_loads][i])
bus_qd[i] = sum(ref[:load][l]["qd"] for l in ref[:bus_loads][i])
end
if length(ref[:bus_shunts][i]) > 0
bus_gs[i] = sum(ref[:shunt][s]["gs"] for s in ref[:bus_shunts][i])
bus_bs[i] = sum(ref[:shunt][s]["bs"] for s in ref[:bus_shunts][i])
end
end
br_g = Dict(i => 0.0 for (i,branch) in ref[:branch])
br_b = Dict(i => 0.0 for (i,branch) in ref[:branch])
br_tr = Dict(i => 0.0 for (i,branch) in ref[:branch])
br_ti = Dict(i => 0.0 for (i,branch) in ref[:branch])
br_ttm = Dict(i => 0.0 for (i,branch) in ref[:branch])
br_g_fr = Dict(i => 0.0 for (i,branch) in ref[:branch])
br_b_fr = Dict(i => 0.0 for (i,branch) in ref[:branch])
br_g_to = Dict(i => 0.0 for (i,branch) in ref[:branch])
br_b_to = Dict(i => 0.0 for (i,branch) in ref[:branch])
for (i,branch) in ref[:branch]
g, b = PowerModels.calc_branch_y(branch)
tr, ti = PowerModels.calc_branch_t(branch)
br_g[i] = g
br_b[i] = b
br_tr[i] = tr
br_ti[i] = ti
br_ttm[i] = tr^2 + ti^2
br_g_fr[i] = branch["g_fr"]
br_b_fr[i] = branch["b_fr"]
br_g_to[i] = branch["g_to"]
br_b_to[i] = branch["b_to"]
end
DataRepresentation(
data,
ref,
var_lookup,
var_init,
var_lb,
var_ub,
ref_gen_idxs,
lookup_pg,
lookup_qg,
lookup_va,
lookup_vm,
lookup_lij,
lookup_p_lij,
lookup_q_lij,
cost_arrs,
f_bus,
t_bus,
ref_bus_idxs,
ref_buses_idxs,
ref_bus_gens,
ref_bus_arcs,
ref_branch_idxs,
ref_arcs_from,
ref_arcs_to,
p_idxmap,
q_idxmap,
bus_pd,
bus_qd,
bus_gs,
bus_bs,
br_g,
br_b,
br_tr,
br_ti,
br_ttm,
br_g_fr,
br_b_fr,
br_g_to,
br_b_to)
end
file_name = "../../benchmarks/OptimizationFrameworks/opf_data/pglib_opf_case5_pjm.m"
dataset = load_and_setup_data(file_name);[info | PowerModels]: extending matpower format with data: areas 1x3
[info | PowerModels]: removing 1 cost terms from generator 4: [4000.0, 0.0]
[info | PowerModels]: removing 1 cost terms from generator 1: [1400.0, 0.0]
[info | PowerModels]: removing 1 cost terms from generator 5: [1000.0, 0.0]
[info | PowerModels]: removing 1 cost terms from generator 2: [1500.0, 0.0]
[info | PowerModels]: removing 1 cost terms from generator 3: [3000.0, 0.0]
[info | PowerModels]: updated generator 4 cost function with order 2 to a f
unction of order 3: [0.0, 4000.0, 0.0]
[info | PowerModels]: updated generator 1 cost function with order 2 to a f
unction of order 3: [0.0, 1400.0, 0.0]
[info | PowerModels]: updated generator 5 cost function with order 2 to a f
unction of order 3: [0.0, 1000.0, 0.0]
[info | PowerModels]: updated generator 2 cost function with order 2 to a f
unction of order 3: [0.0, 1500.0, 0.0]
[info | PowerModels]: updated generator 3 cost function with order 2 to a f
unction of order 3: [0.0, 3000.0, 0.0]Test Setup
Ensure that all objectives and constraints evaluate to the same value on a random vector on the same dataset
test_u0 = [0.6292298794022337, 0.30740951571225206, 0.0215258802699263, 0.38457509230779996, 0.9419186480931858, 0.34961116773074874, 0.875763562401991, 0.3203478635827923, 0.6354060958226175, 0.45537545721771266, 0.3120599359696674, 0.2421238802331842, 0.886455177641366, 0.49797378087768696, 0.652913329799645, 0.03590201299300255, 0.5618806749518928, 0.8142146688533769, 0.3973557130434364, 0.27827135011662674, 0.16456134856048643, 0.7465018431665373, 0.4898329811551083, 0.6966035226583556, 0.7419662648518377, 0.8505905798503723, 0.27102126066405097, 0.1988238097281576, 0.09684601934490256, 0.49238142828542797, 0.1366594202307445, 0.6337080281764231, 0.28814906958008235, 0.5404996094640431, 0.015153517398975858, 0.6338449294034381, 0.5165464961007717, 0.572879113636733, 0.9652420600585092, 0.26535868365228543, 0.865686920119479, 0.38426996353892773, 0.007412077949221274, 0.3889835001514599]
test_obj = 7079.190664351089
test_cons = [0.0215258802699263, -1.0701734802505833, -5.108902216849063, -3.49724505910433, -2.617834191007569, 0.5457423426033834, -0.7150251969424766, -2.473175092089014, -2.071687022809815, -1.5522321037165985, -1.0107399030803794, 3.0047739260369246, 0.2849522377447594, -2.8227966798520674, 3.2236954017592256, 1.0793383525116511, -1.633412293595111, -3.1618224299953224, -0.7775962590542184, 1.7252573527333024, -4.23535583005632, -1.7030832394691608, 1.5810450617647889, -0.33289810365419437, 0.19476447251065077, 1.0688558672739048, 1.563372246165339, 9.915310272572729, 1.4932615291788414, 2.0016715378998793, -1.4038702698147258, -0.8834081057449231, 0.21730536348839036, -7.40879932706212, -1.6000837514115611, 0.8542376821320647, 0.06615508569119477, -0.6077039991323074, 0.6138802155526912, 0.0061762164203837955, -0.3065125522705683, 0.5843454392910835, 0.7251928172073308, 1.2740182727083802, 0.11298343104675009, 0.2518186223833513, 0.4202616621130535, 0.3751697141306502, 0.4019890236200105, 0.5950107614751935, 1.0021074654956683, 0.897077248544158, 0.15136310228960612]53-element Vector{Float64}:
0.0215258802699263
-1.0701734802505833
-5.108902216849063
-3.49724505910433
-2.617834191007569
0.5457423426033834
-0.7150251969424766
-2.473175092089014
-2.071687022809815
-1.5522321037165985
⋮
0.11298343104675009
0.2518186223833513
0.4202616621130535
0.3751697141306502
0.4019890236200105
0.5950107614751935
1.0021074654956683
0.897077248544158
0.15136310228960612Setup and Validations
Now is the setup code for each optimization framework, along with the validation runs on the test case. Any test which fails the validation case, i.e. has x_test_res[1] !≈ test_obj or x_test_res[2] !≈ test_cons should be considered invalidated as this means that the model in that modeling platform does not evaluate to give the same results
Optimization.jl
Constraint optimization implementation reference: https://github.com/SciML/Optimization.jl/blob/master/lib/OptimizationMOI/test/runtests.jl Other AD libraries can be considered: https://docs.sciml.ai/dev/modules/Optimization/API/optimization_function/
import Optimization
import OptimizationMOI
import ModelingToolkit
import Ipopt
import Enzyme
function build_opf_optimization_prob(dataset; adchoice = Optimization.AutoEnzyme())
(;data,
ref,
var_lookup,
var_init,
var_lb,
var_ub,
ref_gen_idxs,
lookup_pg,
lookup_qg,
lookup_va,
lookup_vm,
lookup_lij,
lookup_p_lij,
lookup_q_lij,
cost_arrs,
f_bus,
t_bus,
ref_bus_idxs,
ref_buses_idxs,
ref_bus_gens,
ref_bus_arcs,
ref_branch_idxs,
ref_arcs_from,
ref_arcs_to,
p_idxmap,
q_idxmap,
bus_pd,
bus_qd,
bus_gs,
bus_bs,
br_g,
br_b,
br_tr,
br_ti,
br_ttm,
br_g_fr,
br_b_fr,
br_g_to,
br_b_to) = dataset
#total_callback_time = 0.0
function opf_objective(x, param)
#start = time()
cost = 0.0
for i in ref_gen_idxs
pg = x[lookup_pg[i]]
_cost_arr = cost_arrs[i]
cost += _cost_arr[1]*pg^2 + _cost_arr[2]*pg + _cost_arr[3]
end
#total_callback_time += time() - start
return cost
end
function opf_constraints(ret, x, param)
offsetidx = 0
# va_con
for (reti,i) in enumerate(ref_buses_idxs)
ret[reti + offsetidx] = x[lookup_va[i]]
end
offsetidx += length(ref_buses_idxs)
# @constraint(model,
# sum(p[a] for a in ref[:bus_arcs][i]) ==
# sum(pg[g] for g in ref_bus_gens[i]) -
# sum(load["pd"] for load in bus_loads) -
# sum(shunt["gs"] for shunt in bus_shunts)*x[lookup_vm[i]]^2
# )
# power_balance_p_con
for (reti,i) in enumerate(ref_bus_idxs)
ret[reti + offsetidx] = sum(x[lookup_pg[j]] for j in ref_bus_gens[i]; init=0.0) -
bus_pd[i] -
bus_gs[i]*x[lookup_vm[i]]^2 -
sum(x[p_idxmap[a]] for a in ref_bus_arcs[i])
end
offsetidx += length(ref_bus_idxs)
# @constraint(model,
# sum(q[a] for a in ref[:bus_arcs][i]) ==
# sum(x[lookup_qg[g]] for g in ref_bus_gens[i]) -
# sum(load["qd"] for load in bus_loads) +
# sum(shunt["bs"] for shunt in bus_shunts)*x[lookup_vm[i]]^2
# )
# power_balance_q_con
for (reti,i) in enumerate(ref_bus_idxs)
ret[reti + offsetidx] = sum(x[lookup_qg[j]] for j in ref_bus_gens[i]; init=0.0) -
bus_qd[i] +
bus_bs[i]*x[lookup_vm[i]]^2 -
sum(x[q_idxmap[a]] for a in ref_bus_arcs[i])
end
offsetidx += length(ref_bus_idxs)
# @NLconstraint(model, p_fr == (g+g_fr)/ttm*vm_fr^2 + (-g*tr+b*ti)/ttm*(vm_fr*vm_to*cos(va_fr-va_to)) + (-b*tr-g*ti)/ttm*(vm_fr*vm_to*sin(va_fr-va_to)) )
# power_flow_p_from_con =
for (reti,(l,i,j)) in enumerate(ref_arcs_from)
ret[reti + offsetidx] = (br_g[l]+br_g_fr[l])/br_ttm[l]*x[lookup_vm[f_bus[l]]]^2 +
(-br_g[l]*br_tr[l]+br_b[l]*br_ti[l])/br_ttm[l]*(x[lookup_vm[f_bus[l]]]*x[lookup_vm[t_bus[l]]]*cos(x[lookup_va[f_bus[l]]]-x[lookup_va[t_bus[l]]])) +
(-br_b[l]*br_tr[l]-br_g[l]*br_ti[l])/br_ttm[l]*(x[lookup_vm[f_bus[l]]]*x[lookup_vm[t_bus[l]]]*sin(x[lookup_va[f_bus[l]]]-x[lookup_va[t_bus[l]]])) -
x[p_idxmap[(l,i,j)]]
end
offsetidx += length(ref_arcs_from)
# @NLconstraint(model, p_to == (g+g_to)*vm_to^2 + (-g*tr-b*ti)/ttm*(vm_to*vm_fr*cos(va_to-va_fr)) + (-b*tr+g*ti)/ttm*(vm_to*vm_fr*sin(va_to-va_fr)) )
# power_flow_p_to_con
for (reti,(l,i,j)) in enumerate(ref_arcs_to)
ret[reti + offsetidx] = (br_g[l]+br_g_to[l])*x[lookup_vm[t_bus[l]]]^2 +
(-br_g[l]*br_tr[l]-br_b[l]*br_ti[l])/br_ttm[l]*(x[lookup_vm[t_bus[l]]]*x[lookup_vm[f_bus[l]]]*cos(x[lookup_va[t_bus[l]]]-x[lookup_va[f_bus[l]]])) +
(-br_b[l]*br_tr[l]+br_g[l]*br_ti[l])/br_ttm[l]*(x[lookup_vm[t_bus[l]]]*x[lookup_vm[f_bus[l]]]*sin(x[lookup_va[t_bus[l]]]-x[lookup_va[f_bus[l]]])) -
x[p_idxmap[(l,i,j)]]
end
offsetidx += length(ref_arcs_to)
# @NLconstraint(model, q_fr == -(b+b_fr)/ttm*vm_fr^2 - (-b*tr-g*ti)/ttm*(vm_fr*vm_to*cos(va_fr-va_to)) + (-g*tr+b*ti)/ttm*(vm_fr*vm_to*sin(va_fr-va_to)) )
# power_flow_q_from_con
for (reti,(l,i,j)) in enumerate(ref_arcs_from)
ret[reti + offsetidx] = -(br_b[l]+br_b_fr[l])/br_ttm[l]*x[lookup_vm[f_bus[l]]]^2 -
(-br_b[l]*br_tr[l]-br_g[l]*br_ti[l])/br_ttm[l]*(x[lookup_vm[f_bus[l]]]*x[lookup_vm[t_bus[l]]]*cos(x[lookup_va[f_bus[l]]]-x[lookup_va[t_bus[l]]])) +
(-br_g[l]*br_tr[l]+br_b[l]*br_ti[l])/br_ttm[l]*(x[lookup_vm[f_bus[l]]]*x[lookup_vm[t_bus[l]]]*sin(x[lookup_va[f_bus[l]]]-x[lookup_va[t_bus[l]]])) -
x[q_idxmap[(l,i,j)]]
end
offsetidx += length(ref_arcs_from)
# @NLconstraint(model, q_to == -(b+b_to)*vm_to^2 - (-b*tr+g*ti)/ttm*(vm_to*vm_fr*cos(va_to-va_fr)) + (-g*tr-b*ti)/ttm*(vm_to*vm_fr*sin(va_to-va_fr)) )
# power_flow_q_to_con
for (reti,(l,i,j)) in enumerate(ref_arcs_to)
ret[reti + offsetidx] = -(br_b[l]+br_b_to[l])*x[lookup_vm[t_bus[l]]]^2 -
(-br_b[l]*br_tr[l]+br_g[l]*br_ti[l])/br_ttm[l]*(x[lookup_vm[t_bus[l]]]*x[lookup_vm[f_bus[l]]]*cos(x[lookup_va[t_bus[l]]]-x[lookup_va[f_bus[l]]])) +
(-br_g[l]*br_tr[l]-br_b[l]*br_ti[l])/br_ttm[l]*(x[lookup_vm[t_bus[l]]]*x[lookup_vm[f_bus[l]]]*sin(x[lookup_va[t_bus[l]]]-x[lookup_va[f_bus[l]]])) -
x[q_idxmap[(l,i,j)]]
end
offsetidx += length(ref_arcs_to)
# @constraint(model, va_fr - va_to <= branch["angmax"])
# @constraint(model, va_fr - va_to >= branch["angmin"])
# power_flow_vad_con
for (reti,(l,i,j)) in enumerate(ref_arcs_from)
ret[reti + offsetidx] = x[lookup_va[f_bus[l]]] - x[lookup_va[t_bus[l]]]
end
offsetidx += length(ref_arcs_from)
# @constraint(model, p_fr^2 + q_fr^2 <= branch["rate_a"]^2)
# power_flow_mva_from_con
for (reti,(l,i,j)) in enumerate(ref_arcs_from)
ret[reti + offsetidx] = x[p_idxmap[(l,i,j)]]^2 + x[q_idxmap[(l,i,j)]]^2
end
offsetidx += length(ref_arcs_from)
# @constraint(model, p_to^2 + q_to^2 <= branch["rate_a"]^2)
# power_flow_mva_to_con
for (reti,(l,i,j)) in enumerate(ref_arcs_to)
ret[reti + offsetidx] = x[p_idxmap[(l,i,j)]]^2 + x[q_idxmap[(l,i,j)]]^2
end
offsetidx += length(ref_arcs_to)
@assert offsetidx == length(ret)
nothing
end
con_lbs = Float64[]
con_ubs = Float64[]
#@constraint(model, va[i] == 0)
for (i,bus) in ref[:ref_buses]
push!(con_lbs, 0.0)
push!(con_ubs, 0.0)
end
#power_balance_p_con
for (i,bus) in ref[:bus]
push!(con_lbs, 0.0)
push!(con_ubs, 0.0)
end
#power_balance_q_con
for (i,bus) in ref[:bus]
push!(con_lbs, 0.0)
push!(con_ubs, 0.0)
end
#power_flow_p_from_con
for (l,i,j) in ref[:arcs_from]
push!(con_lbs, 0.0)
push!(con_ubs, 0.0)
end
#power_flow_p_to_con
for (l,i,j) in ref[:arcs_to]
push!(con_lbs, 0.0)
push!(con_ubs, 0.0)
end
#power_flow_q_from_con
for (l,i,j) in ref[:arcs_from]
push!(con_lbs, 0.0)
push!(con_ubs, 0.0)
end
#power_flow_q_to_con
for (l,i,j) in ref[:arcs_to]
push!(con_lbs, 0.0)
push!(con_ubs, 0.0)
end
#power_flow_vad_con
for (l,i,j) in ref[:arcs_from]
branch = ref[:branch][l]
push!(con_lbs, branch["angmin"])
push!(con_ubs, branch["angmax"])
end
#power_flow_mva_from_con
for (l,i,j) in ref[:arcs_from]
branch = ref[:branch][l]
push!(con_lbs, -Inf)
push!(con_ubs, branch["rate_a"]^2)
end
#power_flow_mva_to_con
for (l,i,j) in ref[:arcs_to]
branch = ref[:branch][l]
push!(con_lbs, -Inf)
push!(con_ubs, branch["rate_a"]^2)
end
model_variables = length(var_init)
ret = Array{Float64}(undef, length(con_lbs))
model_constraints = length(con_lbs)
optf = Optimization.OptimizationFunction(opf_objective, adchoice; cons=opf_constraints)
optprob = Optimization.instantiate_function(optf, var_init, adchoice, nothing, length(con_lbs))
optprob = Optimization.OptimizationFunction(opf_objective, adchoice; cons=opf_constraints,
grad = optprob.grad, hess = optprob.hess, cons_j = optprob.cons_j, cons_h = optprob.cons_h,
hess_prototype = optprob.hess_prototype, cons_jac_prototype = optprob.cons_jac_prototype, cons_hess_prototype = optprob.cons_hess_prototype,
hess_colorvec = optprob.hess_colorvec, cons_jac_colorvec = optprob.cons_jac_colorvec, cons_hess_colorvec = optprob.cons_hess_colorvec)
prob = Optimization.OptimizationProblem(optprob, var_init; lb=var_lb, ub=var_ub, lcons=con_lbs, ucons=con_ubs)
end
function solve_opf_optimization(dataset; adchoice = Optimization.AutoSparseReverseDiff(true))
model_build_time = @elapsed prob = build_opf_optimization_prob(dataset; adchoice)
# Correctness tests
ret = zeros(length(prob.lcons))
prob.f.cons(ret, prob.u0, nothing)
@allocated prob.f(prob.u0, nothing) == 0
@allocated prob.f.cons(ret, prob.u0, nothing) == 0
solve_time_with_compilation = @elapsed sol = Optimization.solve(prob, Ipopt.Optimizer(), print_level = PRINT_LEVEL, max_cpu_time = MAX_CPU_TIME)
cost = sol.minimum
feasible = (sol.retcode == Optimization.SciMLBase.ReturnCode.Success)
#println(sol.u) # solution vector
solve_time_without_compilation = @elapsed sol = Optimization.solve(prob, Ipopt.Optimizer(), print_level = PRINT_LEVEL, max_cpu_time = MAX_CPU_TIME)
return (prob,sol),Dict(
"case" => file_name,
"variables" => length(prob.u0),
"constraints" => length(prob.lcons),
"feasible" => feasible,
"cost" => cost,
"time_build" => model_build_time,
"time_solve" => solve_time_without_compilation,
"time_solve_compilation" => solve_time_with_compilation,
)
end
function test_optimization_prob(dataset, test_u0)
prob = build_opf_optimization_prob(dataset)
ret = zeros(length(prob.lcons))
prob.f.cons(ret, test_u0, nothing)
obj = prob.f(test_u0, nothing)
obj, ret
endtest_optimization_prob (generic function with 1 method)optimization_test_res = test_optimization_prob(dataset, test_u0)(7079.190664351089, [0.0215258802699263, -1.0701734802505833, -5.1089022168
49063, -3.49724505910433, -2.617834191007569, 0.5457423426033834, -0.715025
1969424766, -2.473175092089014, -2.071687022809815, -1.5522321037165985 …
1.2740182727083802, 0.11298343104675009, 0.2518186223833513, 0.42026166211
30535, 0.3751697141306502, 0.4019890236200105, 0.5950107614751935, 1.002107
4654956683, 0.897077248544158, 0.15136310228960612])@assert optimization_test_res[1] == test_obj@assert optimization_test_res[2] == test_consModelingToolkit.jl
Showcases symbolic interface to Optimization.jl, through ModelingToolkit.jl. Equivalent to using AutoModelingToolkit as the AD backend in OptimizationFunction, the values for the objectives don't match exactly because of structural simplification.
import PowerModels
import Ipopt
using ModelingToolkit, Optimization, OptimizationMOI
import ModelingToolkit: ≲
function build_opf_mtk_prob(dataset)
(;data, ref) = dataset
vars = Num[]
lb = Float64[]
ub = Float64[]
ModelingToolkit.@variables va[1:maximum(keys(ref[:bus]))]
for i in keys(ref[:bus])
push!(lb, -Inf)
push!(ub, Inf)
end
ModelingToolkit.@variables vm[1:maximum(keys(ref[:bus]))]
for i in keys(ref[:bus])
push!(lb, ref[:bus][i]["vmin"])
push!(ub, ref[:bus][i]["vmax"])
end
vars = vcat(vars, [va[i] for i in keys(ref[:bus])], [vm[i] for i in keys(ref[:bus])])
ModelingToolkit.@variables pg[1:maximum(keys(ref[:gen]))]
for i in keys(ref[:gen])
push!(lb, ref[:gen][i]["pmin"])
push!(ub, ref[:gen][i]["pmax"])
end
ModelingToolkit.@variables qg[1:maximum(keys(ref[:gen]))]
for i in keys(ref[:gen])
push!(lb, ref[:gen][i]["qmin"])
push!(ub, ref[:gen][i]["qmax"])
end
vars = vcat(vars, [pg[i] for i in keys(ref[:gen])], [qg[i] for i in keys(ref[:gen])])
i_inds, j_inds, l_inds = maximum(first.(ref[:arcs])), maximum(getindex.(ref[:arcs], Ref(2))), maximum(last.(ref[:arcs]))
ModelingToolkit.@variables p[1:i_inds, 1:j_inds, 1:l_inds]
ModelingToolkit.@variables q[1:i_inds, 1:j_inds, 1:l_inds]
for (l, i, j) in ref[:arcs]
push!(vars, p[l, i, j])
push!(lb, -ref[:branch][l]["rate_a"])
push!(ub, ref[:branch][l]["rate_a"])
end
for (l, i, j) in ref[:arcs]
push!(vars, q[l, i, j])
push!(lb, -ref[:branch][l]["rate_a"])
push!(ub, ref[:branch][l]["rate_a"])
end
loss = sum(gen["cost"][1] * pg[i]^2 + gen["cost"][2] * pg[i] + gen["cost"][3] for (i, gen) in ref[:gen])
cons = Array{Union{ModelingToolkit.Equation,ModelingToolkit.Inequality}}([])
for (i, bus) in ref[:ref_buses]
push!(cons, va[i] ~ 0)
end
for (i, bus) in ref[:bus]
bus_loads = [ref[:load][l] for l in ref[:bus_loads][i]]
bus_shunts = [ref[:shunt][s] for s in ref[:bus_shunts][i]]
push!(cons,
sum(p[a...] for a in ref[:bus_arcs][i]) ~
(sum(pg[g] for g in ref[:bus_gens][i]; init = 0.0)) -
(sum(load["pd"] for load in bus_loads; init = 0.0)) -
sum(shunt["gs"] for shunt in bus_shunts; init = 0.0)*vm[i]^2
)
push!(cons,
sum(q[a...] for a in ref[:bus_arcs][i]) ~
(sum(qg[g] for g in ref[:bus_gens][i]; init = 0.0)) -
(sum(load["qd"] for load in bus_loads; init = 0.0))
+ sum(shunt["bs"] for shunt in bus_shunts; init = 0.0)*vm[i]^2
)
end
# Branch power flow physics and limit constraints
for (i, branch) in ref[:branch]
f_idx = (i, branch["f_bus"], branch["t_bus"])
t_idx = (i, branch["t_bus"], branch["f_bus"])
p_fr = p[f_idx...]
q_fr = q[f_idx...]
p_to = p[t_idx...]
q_to = q[t_idx...]
vm_fr = vm[branch["f_bus"]]
vm_to = vm[branch["t_bus"]]
va_fr = va[branch["f_bus"]]
va_to = va[branch["t_bus"]]
g, b = PowerModels.calc_branch_y(branch)
tr, ti = PowerModels.calc_branch_t(branch)
ttm = tr^2 + ti^2
g_fr = branch["g_fr"]
b_fr = branch["b_fr"]
g_to = branch["g_to"]
b_to = branch["b_to"]
# From side of the branch flow
push!(cons, p_fr ~ (g + g_fr) / ttm * vm_fr^2 + (-g * tr + b * ti) / ttm * (vm_fr * vm_to * cos(va_fr - va_to)) + (-b * tr - g * ti) / ttm * (vm_fr * vm_to * sin(va_fr - va_to)))
push!(cons, q_fr ~ -(b + b_fr) / ttm * vm_fr^2 - (-b * tr - g * ti) / ttm * (vm_fr * vm_to * cos(va_fr - va_to)) + (-g * tr + b * ti) / ttm * (vm_fr * vm_to * sin(va_fr - va_to)))
# To side of the branch flow
push!(cons, p_to ~ (g + g_to) * vm_to^2 + (-g * tr - b * ti) / ttm * (vm_to * vm_fr * cos(va_to - va_fr)) + (-b * tr + g * ti) / ttm * (vm_to * vm_fr * sin(va_to - va_fr)))
push!(cons, q_to ~ -(b + b_to) * vm_to^2 - (-b * tr + g * ti) / ttm * (vm_to * vm_fr * cos(va_to - va_fr)) + (-g * tr - b * ti) / ttm * (vm_to * vm_fr * sin(va_to - va_fr)))
# Voltage angle difference limit
push!(cons, va_fr - va_to ≲ branch["angmax"])
push!(cons, branch["angmin"] ≲ va_fr - va_to)
# Apparent power limit, from side and to side
push!(cons, p_fr^2 + q_fr^2 ≲ branch["rate_a"]^2)
push!(cons, p_to^2 + q_to^2 ≲ branch["rate_a"]^2)
end
optsys = ModelingToolkit.OptimizationSystem(loss, vars, [], constraints=cons, name=:rosetta)
optsys = ModelingToolkit.structural_simplify(optsys)
u0map = Dict([k => 0.0 for k in collect(optsys.states)])
for key in keys(ref[:bus])
if vm[key] in keys(u0map)
# println(vm[key])
u0map[vm[key]] = 1.0
end
end
inds = Int[]
for k in collect(optsys.states)
push!(inds, findall(x -> isequal(x, k), vars)[1])
end
prob = Optimization.OptimizationProblem(optsys, u0map, lb = lb[inds], ub = ub[inds], grad=true, hess=true, cons_j=true, cons_h=true, cons_sparse=true, sparse=true)
end
function solve_opf_mtk(dataset)
model_build_time = @elapsed prob = build_opf_mtk_prob(dataset)
# @assert prob.f(prob.u0, nothing) == 0.0 #MTK with simplification doesn't evaluate the same
ret = zeros(length(prob.lcons))
prob.f.cons(ret, prob.u0, nothing, )
@allocated prob.f(prob.u0, nothing) == 0
@allocated prob.f.cons(ret, prob.u0, nothing) == 0
solve_time_with_compilation = @elapsed sol = OptimizationMOI.solve(prob, Ipopt.Optimizer())
solve_time_without_compilation = @elapsed sol = OptimizationMOI.solve(prob, Ipopt.Optimizer())
cost = sol.minimum
feasible = (sol.retcode == Optimization.SciMLBase.ReturnCode.Success)
return (prob,sol),Dict(
"case" => file_name,
"variables" => length(prob.u0),
"constraints" => length(prob.lcons),
"feasible" => feasible,
"cost" => cost,
"time_build" => model_build_time,
"time_solve" => solve_time_without_compilation,
"time_solve_compilation" => solve_time_with_compilation,
)
end
function test_mtk_prob(dataset, test_u0)
prob = build_opf_mtk_prob(dataset)
ret = zeros(length(prob.lcons))
prob.f.cons(ret, test_u0, nothing)
obj = prob.f(test_u0, nothing)
obj, ret
endtest_mtk_prob (generic function with 1 method)mtk_test_res = test_optimization_prob(dataset, test_u0)(7079.190664351089, [0.0215258802699263, -1.0701734802505833, -5.1089022168
49063, -3.49724505910433, -2.617834191007569, 0.5457423426033834, -0.715025
1969424766, -2.473175092089014, -2.071687022809815, -1.5522321037165985 …
1.2740182727083802, 0.11298343104675009, 0.2518186223833513, 0.42026166211
30535, 0.3751697141306502, 0.4019890236200105, 0.5950107614751935, 1.002107
4654956683, 0.897077248544158, 0.15136310228960612])@assert mtk_test_res[1] == test_obj@assert mtk_test_res[2] == test_consJuMP.jl
Implementation reference: https://github.com/lanl-ansi/PowerModelsAnnex.jl/blob/master/src/model/ac-opf.jl Only the built-in AD library is supported
import PowerModels
import Ipopt
import JuMP
function build_opf_jump_prob(dataset)
(;data, ref) = dataset
constraints = Any[]
model = JuMP.Model(Ipopt.Optimizer)
vars = [JuMP.@variable(model, va[i in keys(ref[:bus])]),
JuMP.@variable(model, ref[:bus][i]["vmin"] <= vm[i in keys(ref[:bus])] <= ref[:bus][i]["vmax"], start=1.0),
JuMP.@variable(model, ref[:gen][i]["pmin"] <= pg[i in keys(ref[:gen])] <= ref[:gen][i]["pmax"]),
JuMP.@variable(model, ref[:gen][i]["qmin"] <= qg[i in keys(ref[:gen])] <= ref[:gen][i]["qmax"]),
JuMP.@variable(model, -ref[:branch][l]["rate_a"] <= p[(l,i,j) in ref[:arcs]] <= ref[:branch][l]["rate_a"]),
JuMP.@variable(model, -ref[:branch][l]["rate_a"] <= q[(l,i,j) in ref[:arcs]] <= ref[:branch][l]["rate_a"])]
JuMP.@objective(model, Min, sum(gen["cost"][1]*pg[i]^2 + gen["cost"][2]*pg[i] + gen["cost"][3] for (i,gen) in ref[:gen]))
for (i,bus) in ref[:ref_buses]
push!(constraints,JuMP.@constraint(model, va[i] == 0))
end
for (i,bus) in ref[:bus]
bus_loads = [ref[:load][l] for l in ref[:bus_loads][i]]
bus_shunts = [ref[:shunt][s] for s in ref[:bus_shunts][i]]
push!(constraints,JuMP.@constraint(model,
sum(p[a] for a in ref[:bus_arcs][i]) ==
sum(pg[g] for g in ref[:bus_gens][i]) -
sum(load["pd"] for load in bus_loads) -
sum(shunt["gs"] for shunt in bus_shunts)*vm[i]^2
))
push!(constraints,JuMP.@constraint(model,
sum(q[a] for a in ref[:bus_arcs][i]) ==
sum(qg[g] for g in ref[:bus_gens][i]) -
sum(load["qd"] for load in bus_loads) +
sum(shunt["bs"] for shunt in bus_shunts)*vm[i]^2
))
end
# Branch power flow physics and limit constraints
for (i,branch) in ref[:branch]
f_idx = (i, branch["f_bus"], branch["t_bus"])
t_idx = (i, branch["t_bus"], branch["f_bus"])
p_fr = p[f_idx]
q_fr = q[f_idx]
p_to = p[t_idx]
q_to = q[t_idx]
vm_fr = vm[branch["f_bus"]]
vm_to = vm[branch["t_bus"]]
va_fr = va[branch["f_bus"]]
va_to = va[branch["t_bus"]]
g, b = PowerModels.calc_branch_y(branch)
tr, ti = PowerModels.calc_branch_t(branch)
ttm = tr^2 + ti^2
g_fr = branch["g_fr"]
b_fr = branch["b_fr"]
g_to = branch["g_to"]
b_to = branch["b_to"]
# From side of the branch flow
push!(constraints,JuMP.@NLconstraint(model, p_fr == (g+g_fr)/ttm*vm_fr^2 + (-g*tr+b*ti)/ttm*(vm_fr*vm_to*cos(va_fr-va_to)) + (-b*tr-g*ti)/ttm*(vm_fr*vm_to*sin(va_fr-va_to)) ))
push!(constraints,JuMP.@NLconstraint(model, q_fr == -(b+b_fr)/ttm*vm_fr^2 - (-b*tr-g*ti)/ttm*(vm_fr*vm_to*cos(va_fr-va_to)) + (-g*tr+b*ti)/ttm*(vm_fr*vm_to*sin(va_fr-va_to)) ))
# To side of the branch flow
push!(constraints,JuMP.@NLconstraint(model, p_to == (g+g_to)*vm_to^2 + (-g*tr-b*ti)/ttm*(vm_to*vm_fr*cos(va_to-va_fr)) + (-b*tr+g*ti)/ttm*(vm_to*vm_fr*sin(va_to-va_fr)) ))
push!(constraints,JuMP.@NLconstraint(model, q_to == -(b+b_to)*vm_to^2 - (-b*tr+g*ti)/ttm*(vm_to*vm_fr*cos(va_to-va_fr)) + (-g*tr-b*ti)/ttm*(vm_to*vm_fr*sin(va_to-va_fr)) ))
# Voltage angle difference limit
push!(constraints,JuMP.@constraint(model, branch["angmin"] <= va_fr - va_to <= branch["angmax"]))
# Apparent power limit, from side and to side
push!(constraints,JuMP.@constraint(model, p_fr^2 + q_fr^2 <= branch["rate_a"]^2))
push!(constraints,JuMP.@constraint(model, p_to^2 + q_to^2 <= branch["rate_a"]^2))
end
model_variables = JuMP.num_variables(model)
# for consistency with other solvers, skip the variable bounds in the constraint count
non_nl_constraints = sum(JuMP.num_constraints(model, ft, st) for (ft, st) in JuMP.list_of_constraint_types(model) if ft != JuMP.VariableRef)
model_constraints = JuMP.num_nonlinear_constraints(model) + non_nl_constraints
model, vars, constraints
end
function solve_opf_jump(dataset)
model_build_time = @elapsed model = build_opf_jump_prob(dataset)[1]
JuMP.set_attribute(model, "max_cpu_time", MAX_CPU_TIME)
JuMP.set_attribute(model, "print_level", PRINT_LEVEL)
solve_time_with_compilation = @elapsed JuMP.optimize!(model)
solve_time_without_compilation = @elapsed JuMP.optimize!(model)
cost = JuMP.objective_value(model)
feasible = (JuMP.termination_status(model) == JuMP.LOCALLY_SOLVED)
nlp_block = JuMP.MOI.get(model, JuMP.MOI.NLPBlock())
total_callback_time =
nlp_block.evaluator.eval_objective_timer +
nlp_block.evaluator.eval_objective_gradient_timer +
nlp_block.evaluator.eval_constraint_timer +
nlp_block.evaluator.eval_constraint_jacobian_timer +
nlp_block.evaluator.eval_hessian_lagrangian_timer
model_variables = JuMP.num_variables(model)
non_nl_constraints = sum(JuMP.num_constraints(model, ft, st) for (ft, st) in JuMP.list_of_constraint_types(model) if ft != JuMP.VariableRef)
model_constraints = JuMP.num_nonlinear_constraints(model) + non_nl_constraints
return model, Dict(
"case" => file_name,
"variables" => model_variables,
"constraints" => model_constraints,
"feasible" => feasible,
"cost" => cost,
"time_build" => model_build_time,
"time_solve" => solve_time_without_compilation,
"time_solve_compilation" => solve_time_with_compilation,
)
end
function test_jump_prob(dataset, test_u0)
model, vars, constraints = build_opf_jump_prob(dataset)
(;
lookup_pg,
lookup_qg,
lookup_va,
lookup_vm,
lookup_lij,
lookup_p_lij,
lookup_q_lij) = dataset
f = JuMP.objective_function(model)
flatvars = reduce(vcat,[reduce(vcat,vars[i]) for i in 1:length(vars)])
point = Dict()
for v in flatvars
varname, varint = split(JuMP.name(v), "[")
idx = if varint[1] == '('
varint = (parse(Int, varint[2]), parse(Int, varint[5]), parse(Int, varint[8]))
if varname == "p"
lookup_p_lij[findfirst(x->x==varint,lookup_lij)]
elseif varname == "q"
lookup_q_lij[findfirst(x->x==varint,lookup_lij)]
else
error("Invalid $varname, $varint")
end
else
varint = parse(Int, varint[1:end-1])
if varname == "va"
lookup_va[varint]
elseif varname == "pg"
lookup_pg[varint]
elseif varname == "qg"
lookup_qg[varint]
elseif varname == "vm"
lookup_vm[varint]
else
error("Invalid $varname, $varint")
end
end
point[v] = test_u0[idx]
end
obj = JuMP.value(x->point[x], f)
# The JuMP assertion error is because JuMP and optimization.jl build different problems. JuMP builds f(x) == a and optimization.jl builds f(x) - a == 0
# Workaround this for consistent evaluation
# It's not a general purpose approach because only some of the Optimization.jl constraints are written as f(x) - a = 0 .
# Others are written as f(x) <= a, like the p_fr^2 + q_fr^2 <= branch["rate_a"]^2 constraints
primal_value(set::JuMP.MOI.EqualTo) = JuMP.MOI.constant(set)
primal_value(set) = 0.0
function primal_value(f, constraint)
object = JuMP.constraint_object(constraint)
fx = JuMP.value(f, object.func)
return fx - primal_value(object.set)
end
function primal_value(f, constraint::JuMP.NonlinearConstraintRef)
return JuMP.value(f, constraint)
end
obj = JuMP.value(x->point[x], f)
cons = [primal_value(x->point[x], c) for c in constraints]
obj, cons
endtest_jump_prob (generic function with 1 method)jump_test_res = test_jump_prob(dataset, test_u0)(7079.190664351089, [0.0215258802699263, 1.0701734802505833, 0.715025196942
4766, 5.108902216849064, 2.473175092089014, 3.49724505910433, 2.07168702280
9815, 2.6178341910075695, 1.5522321037165985, -0.5457423426033834 … 0.006
1762164203837955, 0.2518186223833513, 0.897077248544158, 1.633412293595111,
-1.4932615291788414, -1.5810450617647889, 1.6000837514115611, -0.306512552
2705683, 0.4202616621130535, 0.15136310228960612])@assert jump_test_res[1] ≈ test_obj@assert sort(abs.(jump_test_res[2])) ≈ sort(abs.(test_cons))NLPModels.jl
Implementation reference: https://juliasmoothoptimizers.github.io/ADNLPModels.jl/stable/tutorial/ Other AD libraries can be considered: https://juliasmoothoptimizers.github.io/ADNLPModels.jl/stable/
import ADNLPModels
import NLPModelsIpopt
function build_opf_nlpmodels_prob(dataset)
(;data, ref) = dataset
bus_pd = Dict(i => 0.0 for (i,bus) in ref[:bus])
bus_qd = Dict(i => 0.0 for (i,bus) in ref[:bus])
bus_gs = Dict(i => 0.0 for (i,bus) in ref[:bus])
bus_bs = Dict(i => 0.0 for (i,bus) in ref[:bus])
for (i,bus) in ref[:bus]
if length(ref[:bus_loads][i]) > 0
bus_pd[i] = sum(ref[:load][l]["pd"] for l in ref[:bus_loads][i])
bus_qd[i] = sum(ref[:load][l]["qd"] for l in ref[:bus_loads][i])
end
if length(ref[:bus_shunts][i]) > 0
bus_gs[i] = sum(ref[:shunt][s]["gs"] for s in ref[:bus_shunts][i])
bus_bs[i] = sum(ref[:shunt][s]["bs"] for s in ref[:bus_shunts][i])
end
end
br_g = Dict(i => 0.0 for (i,branch) in ref[:branch])
br_b = Dict(i => 0.0 for (i,branch) in ref[:branch])
br_tr = Dict(i => 0.0 for (i,branch) in ref[:branch])
br_ti = Dict(i => 0.0 for (i,branch) in ref[:branch])
br_ttm = Dict(i => 0.0 for (i,branch) in ref[:branch])
br_g_fr = Dict(i => 0.0 for (i,branch) in ref[:branch])
br_b_fr = Dict(i => 0.0 for (i,branch) in ref[:branch])
br_g_to = Dict(i => 0.0 for (i,branch) in ref[:branch])
br_b_to = Dict(i => 0.0 for (i,branch) in ref[:branch])
for (i,branch) in ref[:branch]
g, b = PowerModels.calc_branch_y(branch)
tr, ti = PowerModels.calc_branch_t(branch)
br_g[i] = g
br_b[i] = b
br_tr[i] = tr
br_ti[i] = ti
br_ttm[i] = tr^2 + ti^2
br_g_fr[i] = branch["g_fr"]
br_b_fr[i] = branch["b_fr"]
br_g_to[i] = branch["g_to"]
br_b_to[i] = branch["b_to"]
end
var_lookup = Dict{String,Int}()
var_init = Float64[]
var_lb = Float64[]
var_ub = Float64[]
var_idx = 1
for (i,bus) in ref[:bus]
push!(var_init, 0.0) #va
push!(var_lb, -Inf)
push!(var_ub, Inf)
var_lookup["va_$(i)"] = var_idx
var_idx += 1
push!(var_init, 1.0) #vm
push!(var_lb, bus["vmin"])
push!(var_ub, bus["vmax"])
var_lookup["vm_$(i)"] = var_idx
var_idx += 1
end
for (i,gen) in ref[:gen]
push!(var_init, 0.0) #pg
push!(var_lb, gen["pmin"])
push!(var_ub, gen["pmax"])
var_lookup["pg_$(i)"] = var_idx
var_idx += 1
push!(var_init, 0.0) #qg
push!(var_lb, gen["qmin"])
push!(var_ub, gen["qmax"])
var_lookup["qg_$(i)"] = var_idx
var_idx += 1
end
for (l,i,j) in ref[:arcs]
branch = ref[:branch][l]
push!(var_init, 0.0) #p
push!(var_lb, -branch["rate_a"])
push!(var_ub, branch["rate_a"])
var_lookup["p_$(l)_$(i)_$(j)"] = var_idx
var_idx += 1
push!(var_init, 0.0) #q
push!(var_lb, -branch["rate_a"])
push!(var_ub, branch["rate_a"])
var_lookup["q_$(l)_$(i)_$(j)"] = var_idx
var_idx += 1
end
@assert var_idx == length(var_init)+1
#total_callback_time = 0.0
function opf_objective(x)
#start = time()
cost = 0.0
for (i,gen) in ref[:gen]
pg = x[var_lookup["pg_$(i)"]]
cost += gen["cost"][1]*pg^2 + gen["cost"][2]*pg + gen["cost"][3]
end
#total_callback_time += time() - start
return cost
end
function opf_constraints!(cx, x)
#start = time()
va = Dict(i => x[var_lookup["va_$(i)"]] for (i,bus) in ref[:bus])
vm = Dict(i => x[var_lookup["vm_$(i)"]] for (i,bus) in ref[:bus])
pg = Dict(i => x[var_lookup["pg_$(i)"]] for (i,gen) in ref[:gen])
qg = Dict(i => x[var_lookup["qg_$(i)"]] for (i,gen) in ref[:gen])
p = Dict((l,i,j) => x[var_lookup["p_$(l)_$(i)_$(j)"]] for (l,i,j) in ref[:arcs])
q = Dict((l,i,j) => x[var_lookup["q_$(l)_$(i)_$(j)"]] for (l,i,j) in ref[:arcs])
vm_fr = Dict(l => vm[branch["f_bus"]] for (l,branch) in ref[:branch])
vm_to = Dict(l => vm[branch["t_bus"]] for (l,branch) in ref[:branch])
va_fr = Dict(l => va[branch["f_bus"]] for (l,branch) in ref[:branch])
va_to = Dict(l => va[branch["t_bus"]] for (l,branch) in ref[:branch])
# va_con = [va[i] for (i,bus) in ref[:ref_buses]]
k = 0
for (i,bus) in ref[:ref_buses]
k += 1
cx[k] = va[i]
end
# @constraint(model,
# sum(p[a] for a in ref[:bus_arcs][i]) ==
# sum(pg[g] for g in ref[:bus_gens][i]) -
# sum(load["pd"] for load in bus_loads) -
# sum(shunt["gs"] for shunt in bus_shunts)*vm[i]^2
# )
for (i, bus) in ref[:bus]
k += 1
cx[k] = sum(pg[j] for j in ref[:bus_gens][i]; init=0.0) - bus_pd[i] - bus_gs[i]*vm[i]^2 - sum(p[a] for a in ref[:bus_arcs][i])
end
# @constraint(model,
# sum(q[a] for a in ref[:bus_arcs][i]) ==
# sum(qg[g] for g in ref[:bus_gens][i]) -
# sum(load["qd"] for load in bus_loads) +
# sum(shunt["bs"] for shunt in bus_shunts)*vm[i]^2
# )
for (i, bus) in ref[:bus]
k += 1
cx[k] = sum(qg[j] for j in ref[:bus_gens][i]; init=0.0) - bus_qd[i] + bus_bs[i]*vm[i]^2 - sum(q[a] for a in ref[:bus_arcs][i])
end
# @NLconstraint(model, p_fr == (g+g_fr)/ttm*vm_fr^2 + (-g*tr+b*ti)/ttm*(vm_fr*vm_to*cos(va_fr-va_to)) + (-b*tr-g*ti)/ttm*(vm_fr*vm_to*sin(va_fr-va_to)) )
for (l,i,j) in ref[:arcs_from]
k += 1
cx[k] = (br_g[l]+br_g_fr[l])/br_ttm[l]*vm_fr[l]^2 +
(-br_g[l]*br_tr[l]+br_b[l]*br_ti[l])/br_ttm[l]*(vm_fr[l]*vm_to[l]*cos(va_fr[l]-va_to[l])) +
(-br_b[l]*br_tr[l]-br_g[l]*br_ti[l])/br_ttm[l]*(vm_fr[l]*vm_to[l]*sin(va_fr[l]-va_to[l])) -
p[(l,i,j)]
end
# @NLconstraint(model, p_to == (g+g_to)*vm_to^2 + (-g*tr-b*ti)/ttm*(vm_to*vm_fr*cos(va_to-va_fr)) + (-b*tr+g*ti)/ttm*(vm_to*vm_fr*sin(va_to-va_fr)) )
for (l,i,j) in ref[:arcs_to]
k += 1
cx[k] = (br_g[l]+br_g_to[l])*vm_to[l]^2 +
(-br_g[l]*br_tr[l]-br_b[l]*br_ti[l])/br_ttm[l]*(vm_to[l]*vm_fr[l]*cos(va_to[l]-va_fr[l])) +
(-br_b[l]*br_tr[l]+br_g[l]*br_ti[l])/br_ttm[l]*(vm_to[l]*vm_fr[l]*sin(va_to[l]-va_fr[l])) -
p[(l,i,j)]
end
# @NLconstraint(model, q_fr == -(b+b_fr)/ttm*vm_fr^2 - (-b*tr-g*ti)/ttm*(vm_fr*vm_to*cos(va_fr-va_to)) + (-g*tr+b*ti)/ttm*(vm_fr*vm_to*sin(va_fr-va_to)) )
for (l,i,j) in ref[:arcs_from]
k += 1
cx[k] = -(br_b[l]+br_b_fr[l])/br_ttm[l]*vm_fr[l]^2 -
(-br_b[l]*br_tr[l]-br_g[l]*br_ti[l])/br_ttm[l]*(vm_fr[l]*vm_to[l]*cos(va_fr[l]-va_to[l])) +
(-br_g[l]*br_tr[l]+br_b[l]*br_ti[l])/br_ttm[l]*(vm_fr[l]*vm_to[l]*sin(va_fr[l]-va_to[l])) -
q[(l,i,j)]
end
# @NLconstraint(model, q_to == -(b+b_to)*vm_to^2 - (-b*tr+g*ti)/ttm*(vm_to*vm_fr*cos(va_to-va_fr)) + (-g*tr-b*ti)/ttm*(vm_to*vm_fr*sin(va_to-va_fr)) )
for (l,i,j) in ref[:arcs_to]
k += 1
cx[k] = -(br_b[l]+br_b_to[l])*vm_to[l]^2 -
(-br_b[l]*br_tr[l]+br_g[l]*br_ti[l])/br_ttm[l]*(vm_to[l]*vm_fr[l]*cos(va_to[l]-va_fr[l])) +
(-br_g[l]*br_tr[l]-br_b[l]*br_ti[l])/br_ttm[l]*(vm_to[l]*vm_fr[l]*sin(va_to[l]-va_fr[l])) -
q[(l,i,j)]
end
# @constraint(model, va_fr - va_to <= branch["angmax"])
# @constraint(model, va_fr - va_to >= branch["angmin"])
for (l,i,j) in ref[:arcs_from]
k += 1
cx[k] = va_fr[l] - va_to[l]
end
# @constraint(model, p_fr^2 + q_fr^2 <= branch["rate_a"]^2)
for (l,i,j) in ref[:arcs_from]
k += 1
cx[k] = p[(l,i,j)]^2 + q[(l,i,j)]^2
end
# @constraint(model, p_to^2 + q_to^2 <= branch["rate_a"]^2)
for (l,i,j) in ref[:arcs_to]
k += 1
cx[k] = p[(l,i,j)]^2 + q[(l,i,j)]^2
end
#total_callback_time += time() - start
return cx
end
con_lbs = Float64[]
con_ubs = Float64[]
#@constraint(model, va[i] == 0)
for (i,bus) in ref[:ref_buses]
push!(con_lbs, 0.0)
push!(con_ubs, 0.0)
end
#power_balance_p_con
for (i,bus) in ref[:bus]
push!(con_lbs, 0.0)
push!(con_ubs, 0.0)
end
#power_balance_q_con
for (i,bus) in ref[:bus]
push!(con_lbs, 0.0)
push!(con_ubs, 0.0)
end
#power_flow_p_from_con
for (l,i,j) in ref[:arcs_from]
push!(con_lbs, 0.0)
push!(con_ubs, 0.0)
end
#power_flow_p_to_con
for (l,i,j) in ref[:arcs_to]
push!(con_lbs, 0.0)
push!(con_ubs, 0.0)
end
#power_flow_q_from_con
for (l,i,j) in ref[:arcs_from]
push!(con_lbs, 0.0)
push!(con_ubs, 0.0)
end
#power_flow_q_to_con
for (l,i,j) in ref[:arcs_to]
push!(con_lbs, 0.0)
push!(con_ubs, 0.0)
end
#power_flow_vad_con
for (l,i,j) in ref[:arcs_from]
branch = ref[:branch][l]
push!(con_lbs, branch["angmin"])
push!(con_ubs, branch["angmax"])
end
#power_flow_mva_from_con
for (l,i,j) in ref[:arcs_from]
branch = ref[:branch][l]
push!(con_lbs, -Inf)
push!(con_ubs, branch["rate_a"]^2)
end
#power_flow_mva_to_con
for (l,i,j) in ref[:arcs_to]
branch = ref[:branch][l]
push!(con_lbs, -Inf)
push!(con_ubs, branch["rate_a"]^2)
end
model_variables = length(var_init)
model_constraints = length(opf_constraints!(similar(con_lbs), var_init))
#=
backend = ADNLPModels.ADModelBackend(model_variables, opf_objective, model_constraints, opf_constraints!;
gradient_backend = ADNLPModels.ReverseDiffADGradient,
hprod_backend = ADNLPModels.SDTForwardDiffADHvprod,
jprod_backend = ADNLPModels.ForwardDiffADJprod,
jtprod_backend = ADNLPModels.ReverseDiffADJtprod,
jacobian_backend = ADNLPModels.ForwardDiffADJacobian, # SDTSparseADJacobian,
hessian_backend = ADNLPModels.ForwardDiffADHessian, # SparseADJacobian,
ghjvprod_backend = ADNLPModels.ForwardDiffADGHjvprod,
hprod_residual_backend = ADNLPModels.ReverseDiffADHvprod,
jprod_residual_backend = ADNLPModels.ForwardDiffADJprod,
jtprod_residual_backend = ADNLPModels.ReverseDiffADJtprod,
jacobian_residual_backend = ADNLPModels.ForwardDiffADHessian, # SparseADJacobian,
hessian_residual_backend = ADNLPModels.ForwardDiffADHessian
)
=#
nlp = ADNLPModels.ADNLPModel!(opf_objective, var_init, var_lb, var_ub, opf_constraints!, con_lbs, con_ubs, backend = :optimized)
end
function solve_opf_nlpmodels(dataset)
model_build_time = @elapsed nlp = build_opf_nlpmodels_prob(dataset)
solve_time_with_compilation = @elapsed output = NLPModelsIpopt.ipopt(nlp, print_level = PRINT_LEVEL, max_cpu_time = MAX_CPU_TIME)
solve_time_without_compilation = @elapsed output = NLPModelsIpopt.ipopt(nlp, print_level = PRINT_LEVEL, max_cpu_time = MAX_CPU_TIME)
cost = output.objective
feasible = (output.primal_feas <= 1e-6)
model_variables = nlp.meta.nvar
model_constraints = nlp.meta.ncon
return (nlp, output), Dict(
"case" => file_name,
"variables" => model_variables,
"constraints" => model_constraints,
"feasible" => feasible,
"cost" => cost,
"time_build" => model_build_time,
"time_solve" => solve_time_without_compilation,
"time_solve_compilation" => solve_time_with_compilation,
)
end
function test_nlpmodels_prob(dataset, test_u0)
nlp = build_opf_nlpmodels_prob(dataset)
ret = zeros(nlp.meta.ncon)
nlp.c!(ret, test_u0)
obj = nlp.f(test_u0)
obj, ret
endtest_nlpmodels_prob (generic function with 1 method)nlpmodels_test_res = test_nlpmodels_prob(dataset, test_u0)(7079.190664351089, [0.0215258802699263, -1.0701734802505833, -5.1089022168
49063, -3.49724505910433, -2.617834191007569, 0.5457423426033834, -0.715025
1969424766, -2.473175092089014, -2.071687022809815, -1.5522321037165985 …
1.2740182727083802, 0.11298343104675009, 0.2518186223833513, 0.42026166211
30535, 0.3751697141306502, 0.4019890236200105, 0.5950107614751935, 1.002107
4654956683, 0.897077248544158, 0.15136310228960612])@assert nlpmodels_test_res[1] == test_obj@assert nlpmodels_test_res[2] == test_consNonconvex
Implementation reference: https://julianonconvex.github.io/Nonconvex.jl/stable/problem/ Currently does not converge due to an upstream issue with the AD backend Zygote: https://github.com/JuliaNonconvex/Nonconvex.jl/issues/130
import Nonconvex
Nonconvex.@load Ipopt
function build_opf_nonconvex_prob(dataset)
(;data, ref) = dataset
time_model_start = time()
model = Nonconvex.DictModel()
bus_pd = Dict(i => 0.0 for (i,bus) in ref[:bus])
bus_qd = Dict(i => 0.0 for (i,bus) in ref[:bus])
bus_gs = Dict(i => 0.0 for (i,bus) in ref[:bus])
bus_bs = Dict(i => 0.0 for (i,bus) in ref[:bus])
for (i,bus) in ref[:bus]
if length(ref[:bus_loads][i]) > 0
bus_pd[i] = sum(ref[:load][l]["pd"] for l in ref[:bus_loads][i])
bus_qd[i] = sum(ref[:load][l]["qd"] for l in ref[:bus_loads][i])
end
if length(ref[:bus_shunts][i]) > 0
bus_gs[i] = sum(ref[:shunt][s]["gs"] for s in ref[:bus_shunts][i])
bus_bs[i] = sum(ref[:shunt][s]["bs"] for s in ref[:bus_shunts][i])
end
end
br_g = Dict(i => 0.0 for (i,branch) in ref[:branch])
br_b = Dict(i => 0.0 for (i,branch) in ref[:branch])
br_tr = Dict(i => 0.0 for (i,branch) in ref[:branch])
br_ti = Dict(i => 0.0 for (i,branch) in ref[:branch])
br_ttm = Dict(i => 0.0 for (i,branch) in ref[:branch])
br_g_fr = Dict(i => 0.0 for (i,branch) in ref[:branch])
br_b_fr = Dict(i => 0.0 for (i,branch) in ref[:branch])
br_g_to = Dict(i => 0.0 for (i,branch) in ref[:branch])
br_b_to = Dict(i => 0.0 for (i,branch) in ref[:branch])
br_rate_a = Dict(i => 0.0 for (i,branch) in ref[:branch])
br_angmin = Dict(i => 0.0 for (i,branch) in ref[:branch])
br_angmax = Dict(i => 0.0 for (i,branch) in ref[:branch])
for (i,branch) in ref[:branch]
g, b = PowerModels.calc_branch_y(branch)
tr, ti = PowerModels.calc_branch_t(branch)
br_g[i] = g
br_b[i] = b
br_tr[i] = tr
br_ti[i] = ti
br_ttm[i] = tr^2 + ti^2
br_g_fr[i] = branch["g_fr"]
br_b_fr[i] = branch["b_fr"]
br_g_to[i] = branch["g_to"]
br_b_to[i] = branch["b_to"]
br_rate_a[i] = branch["rate_a"]
br_angmin[i] = branch["angmin"]
br_angmax[i] = branch["angmax"]
end
for (i,bus) in ref[:bus]
addvar!(model, "va_$(i)", -Inf, Inf, init=0.0) #va
addvar!(model, "vm_$(i)", bus["vmin"], bus["vmax"], init=1.0) #vm
end
for (i,gen) in ref[:gen]
addvar!(model, "pg_$(i)", gen["pmin"], gen["pmax"], init=0.0) #pg
addvar!(model, "qg_$(i)", gen["qmin"], gen["qmax"], init=0.0) #qg
end
for (l,i,j) in ref[:arcs]
branch = ref[:branch][l]
addvar!(model, "p_$(l)_$(i)_$(j)", -branch["rate_a"], branch["rate_a"], init=0.0) #p
addvar!(model, "q_$(l)_$(i)_$(j)", -branch["rate_a"], branch["rate_a"], init=0.0) #q
end
# JuMP.@objective(model, Min, sum(gen["cost"][1]*pg[i]^2 + gen["cost"][2]*pg[i] + gen["cost"][3] for (i,gen) in ref[:gen]))
function opf_objective(x::OrderedDict)
cost = 0.0
for (i,gen) in ref[:gen]
pg = x["pg_$(i)"]
cost += gen["cost"][1]*pg^2 + gen["cost"][2]*pg + gen["cost"][3]
end
return cost
end
Nonconvex.set_objective!(model, opf_objective)
# JuMP.@constraint(model, va[i] == 0)
function const_ref_bus(x::OrderedDict, i)
return x["va_$(i)"]
end
for (i,bus) in ref[:ref_buses]
add_eq_constraint!(model, x -> const_ref_bus(x,i))
end
# @constraint(model,
# sum(p[a] for a in ref[:bus_arcs][i]) ==
# sum(pg[g] for g in ref[:bus_gens][i]) -
# sum(load["pd"] for load in bus_loads) -
# sum(shunt["gs"] for shunt in bus_shunts)*vm[i]^2
# )
function const_power_balance_p(x::OrderedDict, b)
balance = - bus_pd[b] - bus_gs[b]*x["vm_$(b)"]^2
for (l,i,j) in ref[:bus_arcs][b]
balance -= x["p_$(l)_$(i)_$(j)"]
end
for j in ref[:bus_gens][b]
balance += x["pg_$(j)"]
end
return balance
end
# @constraint(model,
# sum(q[a] for a in ref[:bus_arcs][i]) ==
# sum(qg[g] for g in ref[:bus_gens][i]) -
# sum(load["qd"] for load in bus_loads) +
# sum(shunt["bs"] for shunt in bus_shunts)*vm[i]^2
# )
function const_power_balance_q(x::OrderedDict, b)
balance = - bus_qd[b] + bus_bs[b]*x["vm_$(b)"]^2
for (l,i,j) in ref[:bus_arcs][b]
balance -= x["q_$(l)_$(i)_$(j)"]
end
for j in ref[:bus_gens][b]
balance += x["qg_$(j)"]
end
return balance
end
for (i,bus) in ref[:bus]
add_eq_constraint!(model, x -> const_power_balance_p(x,i))
add_eq_constraint!(model, x -> const_power_balance_q(x,i))
end
# @NLconstraint(model, p_fr == (g+g_fr)/ttm*vm_fr^2 + (-g*tr+b*ti)/ttm*(vm_fr*vm_to*cos(va_fr-va_to)) + (-b*tr-g*ti)/ttm*(vm_fr*vm_to*sin(va_fr-va_to)) )
function const_flow_p_from(x::OrderedDict,l,i,j)
return (br_g[l]+br_g_fr[l])/br_ttm[l]*x["vm_$(i)"]^2 +
(-br_g[l]*br_tr[l]+br_b[l]*br_ti[l])/br_ttm[l]*(x["vm_$(i)"]*x["vm_$(j)"]*cos(x["va_$(i)"]-x["va_$(j)"])) +
(-br_b[l]*br_tr[l]-br_g[l]*br_ti[l])/br_ttm[l]*(x["vm_$(i)"]*x["vm_$(j)"]*sin(x["va_$(i)"]-x["va_$(j)"])) -
x["p_$(l)_$(i)_$(j)"]
end
# @NLconstraint(model, q_fr == -(b+b_fr)/ttm*vm_fr^2 - (-b*tr-g*ti)/ttm*(vm_fr*vm_to*cos(va_fr-va_to)) + (-g*tr+b*ti)/ttm*(vm_fr*vm_to*sin(va_fr-va_to)) )
function const_flow_q_from(x::OrderedDict,l,i,j)
return -(br_b[l]+br_b_fr[l])/br_ttm[l]*x["vm_$(i)"]^2 -
(-br_b[l]*br_tr[l]-br_g[l]*br_ti[l])/br_ttm[l]*(x["vm_$(i)"]*x["vm_$(j)"]*cos(x["va_$(i)"]-x["va_$(j)"])) +
(-br_g[l]*br_tr[l]+br_b[l]*br_ti[l])/br_ttm[l]*(x["vm_$(i)"]*x["vm_$(j)"]*sin(x["va_$(i)"]-x["va_$(j)"])) -
x["q_$(l)_$(i)_$(j)"]
end
# @NLconstraint(model, p_to == (g+g_to)*vm_to^2 + (-g*tr-b*ti)/ttm*(vm_to*vm_fr*cos(va_to-va_fr)) + (-b*tr+g*ti)/ttm*(vm_to*vm_fr*sin(va_to-va_fr)) )
function const_flow_p_to(x::OrderedDict,l,i,j)
return (br_g[l]+br_g_to[l])*x["vm_$(j)"]^2 +
(-br_g[l]*br_tr[l]-br_b[l]*br_ti[l])/br_ttm[l]*(x["vm_$(j)"]*x["vm_$(i)"]*cos(x["va_$(j)"]-x["va_$(i)"])) +
(-br_b[l]*br_tr[l]+br_g[l]*br_ti[l])/br_ttm[l]*(x["vm_$(j)"]*x["vm_$(i)"]*sin(x["va_$(j)"]-x["va_$(i)"])) -
x["p_$(l)_$(j)_$(i)"]
end
# @NLconstraint(model, q_to == -(b+b_to)*vm_to^2 - (-b*tr+g*ti)/ttm*(vm_to*vm_fr*cos(va_to-va_fr)) + (-g*tr-b*ti)/ttm*(vm_to*vm_fr*sin(va_to-va_fr)) )
function const_flow_q_to(x::OrderedDict,l,i,j)
return -(br_b[l]+br_b_to[l])*x["vm_$(j)"]^2 -
(-br_b[l]*br_tr[l]+br_g[l]*br_ti[l])/br_ttm[l]*(x["vm_$(j)"]*x["vm_$(i)"]*cos(x["va_$(j)"]-x["va_$(i)"])) +
(-br_g[l]*br_tr[l]-br_b[l]*br_ti[l])/br_ttm[l]*(x["vm_$(j)"]*x["vm_$(i)"]*sin(x["va_$(j)"]-x["va_$(i)"])) -
x["q_$(l)_$(j)_$(i)"]
end
function const_thermal_limit(x::OrderedDict,l,i,j)
return x["p_$(l)_$(i)_$(j)"]^2 + x["q_$(l)_$(i)_$(j)"]^2 - br_rate_a[l]^2
end
function const_voltage_angle_difference_lb(x::OrderedDict,l,i,j)
return br_angmin[l] - x["va_$(i)"] + x["va_$(j)"]
end
function const_voltage_angle_difference_ub(x::OrderedDict,l,i,j)
return x["va_$(i)"] - x["va_$(j)"] - br_angmax[l]
end
for (l,i,j) in ref[:arcs_from]
add_eq_constraint!(model, x -> const_flow_p_from(x,l,i,j))
add_eq_constraint!(model, x -> const_flow_q_from(x,l,i,j))
add_eq_constraint!(model, x -> const_flow_p_to(x,l,i,j))
add_eq_constraint!(model, x -> const_flow_q_to(x,l,i,j))
add_ineq_constraint!(model, x -> const_thermal_limit(x,l,i,j))
add_ineq_constraint!(model, x -> const_thermal_limit(x,l,j,i))
add_ineq_constraint!(model, x -> const_voltage_angle_difference_lb(x,l,i,j))
add_ineq_constraint!(model, x -> const_voltage_angle_difference_ub(x,l,i,j))
end
model
end
function solve_opf_nonconvex(dataset)
model_build_time = @elapsed model = build_opf_nonconvex_prob(dataset)
solve_time_with_compilation = @elapsed result = Nonconvex.optimize(
model,
IpoptAlg(),
NonconvexCore.getinit(model);
options = IpoptOptions(; first_order=false, symbolic=false, sparse=true, print_level = PRINT_LEVEL, max_cpu_time = MAX_CPU_TIME),
)
solve_time_without_compilation = @elapsed result = Nonconvex.optimize(
model,
IpoptAlg(),
NonconvexCore.getinit(model);
options = IpoptOptions(; first_order=false, symbolic=false, sparse=true, print_level = PRINT_LEVEL, max_cpu_time = MAX_CPU_TIME),
)
cost = result.minimum
feasible = result.status == 0 # just guessing this is correct for Ipopt
model_variables = Nonconvex.NonconvexCore.getnvars(model)
model_constraints = Nonconvex.NonconvexCore.getnconstraints(model)
return (model, result), Dict(
"case" => file_name,
"variables" => model_variables,
"constraints" => model_constraints,
"feasible" => feasible,
"cost" => cost,
"time_build" => model_build_time,
"time_solve" => solve_time_without_compilation,
"time_solve_compilation" => solve_time_with_compilation,
)
end
function test_nonconvex_prob(dataset, test_u0)
model = build_opf_nonconvex_prob(dataset)
(;
lookup_pg,
lookup_qg,
lookup_va,
lookup_vm,
lookup_lij,
lookup_p_lij,
lookup_q_lij) = dataset
point = Dict()
for v in keys(model.init)
varsplit = split(v, "_")
varname = varsplit[1]
varint = parse.(Int, varsplit[2:end])
idx = if varname == "p"
lookup_p_lij[findfirst(x->x==Tuple(varint),lookup_lij)]
elseif varname == "q"
lookup_q_lij[findfirst(x->x==Tuple(varint),lookup_lij)]
elseif varname == "va"
lookup_va[varint[1]]
elseif varname == "pg"
lookup_pg[varint[1]]
elseif varname == "qg"
lookup_qg[varint[1]]
elseif varname == "vm"
lookup_vm[varint[1]]
else
error("Invalid $varname, $varint")
end
point[v] = test_u0[idx]
end
u0 = OrderedDict(keys(model.init) .=> getindex.((point,),keys(model.init)))
obj = model.objective(u0)
cons = vcat(model.eq_constraints(u0), model.ineq_constraints(u0))
obj, cons
endtest_nonconvex_prob (generic function with 1 method)nonconvex_test_res = test_nonconvex_prob(dataset, test_u0)(7079.190664351089, [0.0215258802699263, -1.0701734802505833, -0.7150251969
424766, -5.108902216849064, -2.4731750920890136, -3.4972450591043294, -2.07
1687022809815, -2.617834191007569, -1.5522321037165985, 0.5457423426033834
… -1.13747899115099, 0.09028143995439242, -17.895781377616647, -17.250522
75145584, -0.5297749920186825, -0.517422559177915, -15.579738337886946, -15
.848636897710394, -0.21708622332773042, -0.8301113278688671])@assert nonconvex_test_res[1] ≈ test_obj@assert sort(abs.(nonconvex_test_res[2])) ≈ sort(abs.(test_cons))Error: DimensionMismatch: dimensions must match: a has dims (Base.OneTo(59)
,), b has dims (Base.OneTo(53),), mismatch at 1println(sort(abs.(nonconvex_test_res[2])))[0.0215258802699263, 0.08410522353400862, 0.09028143995439242, 0.1947644725
1065077, 0.21708622332773042, 0.21730536348839036, 0.2849522377447594, 0.33
063890653376593, 0.33289810365419437, 0.45744368990710405, 0.51742255917791
5, 0.5297749920186825, 0.5457423426033834, 0.5897538612894937, 0.7150251969
424766, 0.7775962590542184, 0.8301113278688671, 0.8834081057449231, 1.01073
99030803794, 1.0688558672739048, 1.0701734802505833, 1.0793383525116511, 1.
1313027747306061, 1.13747899115099, 1.3778364577303637, 1.4038702698147258,
1.4932615291788414, 1.5522321037165985, 1.563372246165339, 1.5810450617647
889, 1.6000837514115611, 1.633412293595111, 1.7030832394691608, 1.725257352
7333024, 2.0016715378998793, 2.071687022809815, 2.4731750920890136, 2.61783
4191007569, 2.8227966798520674, 3.0047739260369246, 3.1618224299953224, 3.2
236954017592256, 3.4972450591043294, 4.23535583005632, 4.48598172729162, 5.
108902216849064, 5.164989238524806, 7.40879932706212, 9.915310272572729, 15
.579738337886946, 15.848636897710394, 17.145492534504328, 17.25052275145584
, 17.422407182792668, 17.563254560708913, 17.745610976379986, 17.7724302858
69348, 17.895781377616647, 18.034616568953247]println(sort(abs.(test_cons)))[0.0061762164203837955, 0.0215258802699263, 0.06615508569119477, 0.11298343
104675009, 0.15136310228960612, 0.19476447251065077, 0.21730536348839036, 0
.2518186223833513, 0.2849522377447594, 0.3065125522705683, 0.33289810365419
437, 0.3751697141306502, 0.4019890236200105, 0.4202616621130535, 0.54574234
26033834, 0.5843454392910835, 0.5950107614751935, 0.6077039991323074, 0.613
8802155526912, 0.7150251969424766, 0.7251928172073308, 0.7775962590542184,
0.8542376821320647, 0.8834081057449231, 0.897077248544158, 1.00210746549566
83, 1.0107399030803794, 1.0688558672739048, 1.0701734802505833, 1.079338352
5116511, 1.2740182727083802, 1.4038702698147258, 1.4932615291788414, 1.5522
321037165985, 1.563372246165339, 1.5810450617647889, 1.6000837514115611, 1.
633412293595111, 1.7030832394691608, 1.7252573527333024, 2.0016715378998793
, 2.071687022809815, 2.473175092089014, 2.617834191007569, 2.82279667985206
74, 3.0047739260369246, 3.1618224299953224, 3.2236954017592256, 3.497245059
10433, 4.23535583005632, 5.108902216849063, 7.40879932706212, 9.91531027257
2729]Optim.jl
Implementation reference: https://julianlsolvers.github.io/Optim.jl/stable/#examples/generated/ipnewton_basics/ Currently does not converge to a feasible point, root cause in unclear debug/optim-debug.jl can be used to confirm it will converge if given a suitable starting point
import Optim
function build_opf_optim_prob(dataset)
(;data, ref) = dataset
bus_pd = Dict(i => 0.0 for (i,bus) in ref[:bus])
bus_qd = Dict(i => 0.0 for (i,bus) in ref[:bus])
bus_gs = Dict(i => 0.0 for (i,bus) in ref[:bus])
bus_bs = Dict(i => 0.0 for (i,bus) in ref[:bus])
for (i,bus) in ref[:bus]
if length(ref[:bus_loads][i]) > 0
bus_pd[i] = sum(ref[:load][l]["pd"] for l in ref[:bus_loads][i])
bus_qd[i] = sum(ref[:load][l]["qd"] for l in ref[:bus_loads][i])
end
if length(ref[:bus_shunts][i]) > 0
bus_gs[i] = sum(ref[:shunt][s]["gs"] for s in ref[:bus_shunts][i])
bus_bs[i] = sum(ref[:shunt][s]["bs"] for s in ref[:bus_shunts][i])
end
end
br_g = Dict(i => 0.0 for (i,branch) in ref[:branch])
br_b = Dict(i => 0.0 for (i,branch) in ref[:branch])
br_tr = Dict(i => 0.0 for (i,branch) in ref[:branch])
br_ti = Dict(i => 0.0 for (i,branch) in ref[:branch])
br_ttm = Dict(i => 0.0 for (i,branch) in ref[:branch])
br_g_fr = Dict(i => 0.0 for (i,branch) in ref[:branch])
br_b_fr = Dict(i => 0.0 for (i,branch) in ref[:branch])
br_g_to = Dict(i => 0.0 for (i,branch) in ref[:branch])
br_b_to = Dict(i => 0.0 for (i,branch) in ref[:branch])
for (i,branch) in ref[:branch]
g, b = PowerModels.calc_branch_y(branch)
tr, ti = PowerModels.calc_branch_t(branch)
br_g[i] = g
br_b[i] = b
br_tr[i] = tr
br_ti[i] = ti
br_ttm[i] = tr^2 + ti^2
br_g_fr[i] = branch["g_fr"]
br_b_fr[i] = branch["b_fr"]
br_g_to[i] = branch["g_to"]
br_b_to[i] = branch["b_to"]
end
var_lookup = Dict{String,Int}()
var_init = Float64[]
var_lb = Float64[]
var_ub = Float64[]
var_idx = 1
for (i,bus) in ref[:bus]
push!(var_init, 0.0) #va
push!(var_lb, -Inf)
push!(var_ub, Inf)
var_lookup["va_$(i)"] = var_idx
var_idx += 1
push!(var_init, 1.0) #vm
push!(var_lb, bus["vmin"])
push!(var_ub, bus["vmax"])
var_lookup["vm_$(i)"] = var_idx
var_idx += 1
end
for (i,gen) in ref[:gen]
#push!(var_init, 0.0) #pg
push!(var_init, (gen["pmax"]+gen["pmin"])/2) # non-standard start
push!(var_lb, gen["pmin"])
push!(var_ub, gen["pmax"])
var_lookup["pg_$(i)"] = var_idx
var_idx += 1
#push!(var_init, 0.0) #qg
push!(var_init, (gen["qmax"]+gen["qmin"])/2) # non-standard start
push!(var_lb, gen["qmin"])
push!(var_ub, gen["qmax"])
var_lookup["qg_$(i)"] = var_idx
var_idx += 1
end
for (l,i,j) in ref[:arcs]
branch = ref[:branch][l]
push!(var_init, 0.0) #p
push!(var_lb, -branch["rate_a"])
push!(var_ub, branch["rate_a"])
var_lookup["p_$(l)_$(i)_$(j)"] = var_idx
var_idx += 1
push!(var_init, 0.0) #q
push!(var_lb, -branch["rate_a"])
push!(var_ub, branch["rate_a"])
var_lookup["q_$(l)_$(i)_$(j)"] = var_idx
var_idx += 1
end
@assert var_idx == length(var_init)+1
#total_callback_time = 0.0
function opf_objective(x)
#start = time()
cost = 0.0
for (i,gen) in ref[:gen]
pg = x[var_lookup["pg_$(i)"]]
cost += gen["cost"][1]*pg^2 + gen["cost"][2]*pg + gen["cost"][3]
end
#total_callback_time += time() - start
return cost
end
function opf_constraints(c,x)
#start = time()
va = Dict(i => x[var_lookup["va_$(i)"]] for (i,bus) in ref[:bus])
vm = Dict(i => x[var_lookup["vm_$(i)"]] for (i,bus) in ref[:bus])
pg = Dict(i => x[var_lookup["pg_$(i)"]] for (i,gen) in ref[:gen])
qg = Dict(i => x[var_lookup["qg_$(i)"]] for (i,gen) in ref[:gen])
p = Dict((l,i,j) => x[var_lookup["p_$(l)_$(i)_$(j)"]] for (l,i,j) in ref[:arcs])
q = Dict((l,i,j) => x[var_lookup["q_$(l)_$(i)_$(j)"]] for (l,i,j) in ref[:arcs])
vm_fr = Dict(l => vm[branch["f_bus"]] for (l,branch) in ref[:branch])
vm_to = Dict(l => vm[branch["t_bus"]] for (l,branch) in ref[:branch])
va_fr = Dict(l => va[branch["f_bus"]] for (l,branch) in ref[:branch])
va_to = Dict(l => va[branch["t_bus"]] for (l,branch) in ref[:branch])
va_con = [va[i] for (i,bus) in ref[:ref_buses]]
# @constraint(model,
# sum(p[a] for a in ref[:bus_arcs][i]) ==
# sum(pg[g] for g in ref[:bus_gens][i]) -
# sum(load["pd"] for load in bus_loads) -
# sum(shunt["gs"] for shunt in bus_shunts)*vm[i]^2
# )
power_balance_p_con = [
sum(pg[j] for j in ref[:bus_gens][i]; init=0.0) -
bus_pd[i] -
bus_gs[i]*vm[i]^2 -
sum(p[a] for a in ref[:bus_arcs][i])
for (i,bus) in ref[:bus]
]
# @constraint(model,
# sum(q[a] for a in ref[:bus_arcs][i]) ==
# sum(qg[g] for g in ref[:bus_gens][i]) -
# sum(load["qd"] for load in bus_loads) +
# sum(shunt["bs"] for shunt in bus_shunts)*vm[i]^2
# )
power_balance_q_con = [
sum(qg[j] for j in ref[:bus_gens][i]; init=0.0) -
bus_qd[i] +
bus_bs[i]*vm[i]^2 -
sum(q[a] for a in ref[:bus_arcs][i])
for (i,bus) in ref[:bus]
]
# @NLconstraint(model, p_fr == (g+g_fr)/ttm*vm_fr^2 + (-g*tr+b*ti)/ttm*(vm_fr*vm_to*cos(va_fr-va_to)) + (-b*tr-g*ti)/ttm*(vm_fr*vm_to*sin(va_fr-va_to)) )
power_flow_p_from_con = [
(br_g[l]+br_g_fr[l])/br_ttm[l]*vm_fr[l]^2 +
(-br_g[l]*br_tr[l]+br_b[l]*br_ti[l])/br_ttm[l]*(vm_fr[l]*vm_to[l]*cos(va_fr[l]-va_to[l])) +
(-br_b[l]*br_tr[l]-br_g[l]*br_ti[l])/br_ttm[l]*(vm_fr[l]*vm_to[l]*sin(va_fr[l]-va_to[l])) -
p[(l,i,j)]
for (l,i,j) in ref[:arcs_from]
]
# @NLconstraint(model, p_to == (g+g_to)*vm_to^2 + (-g*tr-b*ti)/ttm*(vm_to*vm_fr*cos(va_to-va_fr)) + (-b*tr+g*ti)/ttm*(vm_to*vm_fr*sin(va_to-va_fr)) )
power_flow_p_to_con = [
(br_g[l]+br_g_to[l])*vm_to[l]^2 +
(-br_g[l]*br_tr[l]-br_b[l]*br_ti[l])/br_ttm[l]*(vm_to[l]*vm_fr[l]*cos(va_to[l]-va_fr[l])) +
(-br_b[l]*br_tr[l]+br_g[l]*br_ti[l])/br_ttm[l]*(vm_to[l]*vm_fr[l]*sin(va_to[l]-va_fr[l])) -
p[(l,i,j)]
for (l,i,j) in ref[:arcs_to]
]
# @NLconstraint(model, q_fr == -(b+b_fr)/ttm*vm_fr^2 - (-b*tr-g*ti)/ttm*(vm_fr*vm_to*cos(va_fr-va_to)) + (-g*tr+b*ti)/ttm*(vm_fr*vm_to*sin(va_fr-va_to)) )
power_flow_q_from_con = [
-(br_b[l]+br_b_fr[l])/br_ttm[l]*vm_fr[l]^2 -
(-br_b[l]*br_tr[l]-br_g[l]*br_ti[l])/br_ttm[l]*(vm_fr[l]*vm_to[l]*cos(va_fr[l]-va_to[l])) +
(-br_g[l]*br_tr[l]+br_b[l]*br_ti[l])/br_ttm[l]*(vm_fr[l]*vm_to[l]*sin(va_fr[l]-va_to[l])) -
q[(l,i,j)]
for (l,i,j) in ref[:arcs_from]
]
# @NLconstraint(model, q_to == -(b+b_to)*vm_to^2 - (-b*tr+g*ti)/ttm*(vm_to*vm_fr*cos(va_to-va_fr)) + (-g*tr-b*ti)/ttm*(vm_to*vm_fr*sin(va_to-va_fr)) )
power_flow_q_to_con = [
-(br_b[l]+br_b_to[l])*vm_to[l]^2 -
(-br_b[l]*br_tr[l]+br_g[l]*br_ti[l])/br_ttm[l]*(vm_to[l]*vm_fr[l]*cos(va_to[l]-va_fr[l])) +
(-br_g[l]*br_tr[l]-br_b[l]*br_ti[l])/br_ttm[l]*(vm_to[l]*vm_fr[l]*sin(va_to[l]-va_fr[l])) -
q[(l,i,j)]
for (l,i,j) in ref[:arcs_to]
]
# @constraint(model, va_fr - va_to <= branch["angmax"])
# @constraint(model, va_fr - va_to >= branch["angmin"])
power_flow_vad_con = [
va_fr[l] - va_to[l]
for (l,i,j) in ref[:arcs_from]
]
# @constraint(model, p_fr^2 + q_fr^2 <= branch["rate_a"]^2)
power_flow_mva_from_con = [
p[(l,i,j)]^2 + q[(l,i,j)]^2
for (l,i,j) in ref[:arcs_from]
]
# @constraint(model, p_to^2 + q_to^2 <= branch["rate_a"]^2)
power_flow_mva_to_con = [
p[(l,i,j)]^2 + q[(l,i,j)]^2
for (l,i,j) in ref[:arcs_to]
]
c .= [
va_con...,
power_balance_p_con...,
power_balance_q_con...,
power_flow_p_from_con...,
power_flow_p_to_con...,
power_flow_q_from_con...,
power_flow_q_to_con...,
power_flow_vad_con...,
power_flow_mva_from_con...,
power_flow_mva_to_con...,
]
#total_callback_time += time() - start
return c
end
con_lbs = Float64[]
con_ubs = Float64[]
#@constraint(model, va[i] == 0)
for (i,bus) in ref[:ref_buses]
push!(con_lbs, 0.0)
push!(con_ubs, 0.0)
end
#power_balance_p_con
for (i,bus) in ref[:bus]
push!(con_lbs, 0.0)
push!(con_ubs, 0.0)
#push!(con_lbs, -Inf)
#push!(con_ubs, Inf)
end
#power_balance_q_con
for (i,bus) in ref[:bus]
push!(con_lbs, 0.0)
push!(con_ubs, 0.0)
#push!(con_lbs, -Inf)
#push!(con_ubs, Inf)
end
#power_flow_p_from_con
for (l,i,j) in ref[:arcs_from]
push!(con_lbs, 0.0)
push!(con_ubs, 0.0)
end
#power_flow_p_to_con
for (l,i,j) in ref[:arcs_to]
push!(con_lbs, 0.0)
push!(con_ubs, 0.0)
end
#power_flow_q_from_con
for (l,i,j) in ref[:arcs_from]
push!(con_lbs, 0.0)
push!(con_ubs, 0.0)
end
#power_flow_q_to_con
for (l,i,j) in ref[:arcs_to]
push!(con_lbs, 0.0)
push!(con_ubs, 0.0)
end
#power_flow_vad_con
for (l,i,j) in ref[:arcs_from]
branch = ref[:branch][l]
push!(con_lbs, branch["angmin"])
push!(con_ubs, branch["angmax"])
end
#power_flow_mva_from_con
for (l,i,j) in ref[:arcs_from]
branch = ref[:branch][l]
push!(con_lbs, -Inf)
push!(con_ubs, branch["rate_a"]^2)
end
#power_flow_mva_to_con
for (l,i,j) in ref[:arcs_to]
branch = ref[:branch][l]
push!(con_lbs, -Inf)
push!(con_ubs, branch["rate_a"]^2)
end
df = Optim.TwiceDifferentiable(opf_objective, var_init)
dfc = Optim.TwiceDifferentiableConstraints(opf_constraints, var_lb, var_ub, con_lbs, con_ubs)
df, dfc, var_init, con_lbs, con_ubs
end
function solve_opf_optim(dataset)
model_build_time = @elapsed df, dfc, var_init, con_lbs, con_ubs = build_opf_optim_prob(dataset)
options = Optim.Options(show_trace=PRINT_LEVEL != 0,time_limit=MAX_CPU_TIME)
solve_time_with_compilation = @elapsed res = Optim.optimize(df, dfc, var_init, Optim.IPNewton(), options)
solve_time_without_compilation = @elapsed res = Optim.optimize(df, dfc, var_init, Optim.IPNewton(), options)
sol = res.minimizer
cost = res.minimum
# NOTE: confirmed these constraint violations can be eliminated
# if a better starting point is used
sol_eval = dfc.c!(zeros(dfc.bounds.nc), sol)
vio_lb = [max(v,0) for v in (con_lbs .- sol_eval)]
vio_ub = [max(v,0) for v in (sol_eval .- con_ubs)]
const_vio = vio_lb .+ vio_ub
constraint_tol = 1e-6
feasible = (sum(const_vio) <= constraint_tol)
return (res,), Dict(
"case" => file_name,
"variables" => length(var_init),
"constraints" => dfc.bounds.nc,
"feasible" => feasible,
"cost" => cost,
"time_build" => model_build_time,
"time_solve" => solve_time_without_compilation,
"time_solve_compilation" => solve_time_with_compilation,
)
end
function test_optim_prob(dataset, test_u0)
df, dfc, var_init, con_lbs, con_ubs = build_opf_optim_prob(dataset)
obj = df.f(test_u0)
cons = dfc.c!(zeros(dfc.bounds.nc), test_u0)
obj, cons
endtest_optim_prob (generic function with 1 method)optim_test_res = test_optim_prob(dataset, test_u0)(7079.190664351089, [0.0215258802699263, -1.0701734802505833, -5.1089022168
49063, -3.49724505910433, -2.617834191007569, 0.5457423426033834, -0.715025
1969424766, -2.473175092089014, -2.071687022809815, -1.5522321037165985 …
1.2740182727083802, 0.11298343104675009, 0.2518186223833513, 0.42026166211
30535, 0.3751697141306502, 0.4019890236200105, 0.5950107614751935, 1.002107
4654956683, 0.897077248544158, 0.15136310228960612])@assert optim_test_res[1] == test_obj@assert optim_test_res[2] == test_consCASADI
Implementation reference: https://github.com/lanl-ansi/PowerModelsAnnex.jl/blob/master/src/model/ac-opf.jl
CASADI Segfaults so removed for now.
import PowerModels
import PythonCall
import CondaPkg
CondaPkg.add("casadi")
function solve_opf_casadi(dataset)
(;data, ref) = dataset
time_model_start = time()
casadi = PythonCall.pyimport("casadi")
x, x0, lbx, ubx, cons, lbg, ubg = [], [], [], [], [], [], []
va, vm = Dict{Int,Any}(), Dict{Int,Any}()
for (k, _) in ref[:bus]
va[k] = casadi.SX.sym("va$k")
push!(x, va[k])
push!(x0, 0.0)
push!(lbx, -casadi.inf)
push!(ubx, casadi.inf)
vm[k] = casadi.SX.sym("vm$k")
push!(x, vm[k])
push!(x0, 1.0)
push!(lbx, ref[:bus][k]["vmin"])
push!(ubx, ref[:bus][k]["vmax"])
end
pg, qg = Dict{Int,Any}(), Dict{Int,Any}()
for (k, ) in ref[:gen]
pg[k] = casadi.SX.sym("pg$k")
push!(x, pg[k])
push!(x0, 0.0)
push!(lbx, ref[:gen][k]["pmin"])
push!(ubx, ref[:gen][k]["pmax"])
qg[k] = casadi.SX.sym("qg$k")
push!(x, qg[k])
push!(x0, 0.0)
push!(lbx, ref[:gen][k]["qmin"])
push!(ubx, ref[:gen][k]["qmax"])
end
p, q = Dict{NTuple{3,Int},Any}(), Dict{NTuple{3,Int},Any}()
for k in ref[:arcs]
a = ref[:branch][k[1]]["rate_a"]
p[k] = casadi.SX.sym("p$k")
push!(x, p[k])
push!(x0, 0.0)
push!(lbx, -a)
push!(ubx, a)
q[k] = casadi.SX.sym("q$k")
push!(x, q[k])
push!(x0, 0.0)
push!(lbx, -a)
push!(ubx, a)
end
f = sum(
cost["cost"][1] * pg[k]^2 +
cost["cost"][2] * pg[k] +
cost["cost"][3] for (k, cost) in ref[:gen]
)
for (k, _) in ref[:ref_buses]
push!(cons, va[k])
push!(lbg, 0)
push!(ubg, 0)
end
for (i, _) in ref[:bus]
bus_loads = [ref[:load][l] for l in ref[:bus_loads][i]]
bus_shunts = [ref[:shunt][s] for s in ref[:bus_shunts][i]]
push!(
cons,
sum(p[k] for k in ref[:bus_arcs][i]) -
sum(pg[g] for g in ref[:bus_gens][i]; init = 0) +
sum(load["pd"] for load in bus_loads; init = 0) +
sum(shunt["gs"] for shunt in bus_shunts; init = 0) * vm[i]^2
)
push!(lbg, 0)
push!(ubg, 0)
push!(
cons,
sum(q[k] for k in ref[:bus_arcs][i]) -
sum(qg[g] for g in ref[:bus_gens][i]; init = 0) +
sum(load["qd"] for load in bus_loads; init = 0) -
sum(shunt["bs"] for shunt in bus_shunts; init = 0) * vm[i]^2
)
push!(lbg, 0)
push!(ubg, 0)
end
for (i, branch) in ref[:branch]
f_idx = (i, branch["f_bus"], branch["t_bus"])
t_idx = (i, branch["t_bus"], branch["f_bus"])
p_fr = p[f_idx]
q_fr = q[f_idx]
p_to = p[t_idx]
q_to = q[t_idx]
vm_fr = vm[branch["f_bus"]]
vm_to = vm[branch["t_bus"]]
va_fr = va[branch["f_bus"]]
va_to = va[branch["t_bus"]]
g, b = PowerModels.calc_branch_y(branch)
tr, ti = PowerModels.calc_branch_t(branch)
ttm = tr^2 + ti^2
g_fr = branch["g_fr"]
b_fr = branch["b_fr"]
g_to = branch["g_to"]
b_to = branch["b_to"]
push!(
cons,
(g+g_fr)/ttm*vm_fr^2 + (-g*tr+b*ti)/ttm*(vm_fr*vm_to*casadi.cos(va_fr-va_to)) + (-b*tr-g*ti)/ttm*(vm_fr*vm_to*casadi.sin(va_fr-va_to)) - p_fr
)
push!(
cons,
-(b+b_fr)/ttm*vm_fr^2 - (-b*tr-g*ti)/ttm*(vm_fr*vm_to*casadi.cos(va_fr-va_to)) + (-g*tr+b*ti)/ttm*(vm_fr*vm_to*casadi.sin(va_fr-va_to)) - q_fr
)
push!(
cons,
(g+g_to)*vm_to^2 + (-g*tr-b*ti)/ttm*(vm_to*vm_fr*casadi.cos(va_to-va_fr)) + (-b*tr+g*ti)/ttm*(vm_to*vm_fr*casadi.sin(va_to-va_fr)) - p_to
)
push!(
cons,
-(b+b_to)*vm_to^2 - (-b*tr+g*ti)/ttm*(vm_to*vm_fr*casadi.cos(va_to-va_fr)) + (-g*tr-b*ti)/ttm*(vm_to*vm_fr*casadi.sin(va_to-va_fr)) - q_to
)
for i in 1:4
push!(lbg, 0)
push!(ubg, 0)
end
push!(cons, va_fr - va_to)
push!(lbg, branch["angmin"])
push!(ubg, branch["angmax"])
push!(cons, p_fr^2 + q_fr^2)
push!(lbg, -casadi.inf)
push!(ubg, branch["rate_a"]^2)
push!(cons, p_to^2 + q_to^2)
push!(lbg, -casadi.inf)
push!(ubg, branch["rate_a"]^2)
end
nlp = Dict("x" => casadi.vcat(x), "f" => f, "g" => casadi.vcat(cons))
options = PythonCall.pydict(Dict("error_on_fail" => true))
model = casadi.nlpsol("model", "ipopt", PythonCall.pydict(nlp), options)
model_variables = length(x)
model_constraints = length(lbg)
model_build_time = time() - time_model_start
time_solve_start = time()
solution = model(; lbx = lbx, ubx = ubx, lbg = lbg, ubg = ubg, x0 = x0)
cost = PythonCall.pyconvert(Float64, (PythonCall.pyfloat(solution["f"])))
feasible = true # error if not feasible
solve_time = time() - time_solve_start
total_time = time() - time_data_start
println("")
println("\033[1mSummary\033[0m")
println(" case........: $(file_name)")
println(" variables...: $(model_variables)")
println(" constraints.: $(model_constraints)")
println(" feasible....: $(feasible)")
println(" cost........: $(round(Int, cost))")
println(" total time..: $(total_time)")
println(" data time.: $(data_load_time)")
println(" build time: $(model_build_time)")
println(" solve time: $(solve_time)")
# println(" callbacks: $(total_callback_time)")
println("")
return Dict(
"case" => file_name,
"variables" => model_variables,
"constraints" => model_constraints,
"feasible" => feasible,
"cost" => cost,
"time_total" => total_time,
"time_data" => data_load_time,
"time_build" => model_build_time,
"time_solve" => solve_time,
#"time_callbacks" => TBD,
)
end
solve_opf_casadi("../../benchmarks/OptimizationFrameworks/opf_data/pglib_opf_case5_pjm.m")Test the Benchmarking
file_name = "../../benchmarks/OptimizationFrameworks/opf_data/pglib_opf_case5_pjm.m"
dataset = load_and_setup_data(file_name);Error: IOError: stream is closed or unusablemodel, res = solve_opf_optimization(dataset);
res***************************************************************************
***
This program contains Ipopt, a library for large-scale nonlinear optimizati
on.
Ipopt is released as open source code under the Eclipse Public License (EP
L).
For more information visit https://github.com/coin-or/Ipopt
***************************************************************************
***
Dict{String, Any} with 8 entries:
"cost" => 17551.9
"variables" => 44
"constraints" => 53
"case" => "../../benchmarks/OptimizationFrameworks/opf_
data…
"time_build" => 9.85608
"time_solve_compilation" => 7.83109
"time_solve" => 0.0995637
"feasible" => truemodel, res = solve_opf_jump(dataset);
resDict{String, Any} with 8 entries:
"cost" => 17551.9
"variables" => 44
"constraints" => 53
"case" => "../../benchmarks/OptimizationFrameworks/opf_
data…
"time_build" => 0.00204758
"time_solve_compilation" => 0.954648
"time_solve" => 0.0160452
"feasible" => truemodel, res = solve_opf_nlpmodels(dataset);
resDict{String, Any} with 8 entries:
"cost" => 17551.9
"variables" => 44
"constraints" => 53
"case" => "../../benchmarks/OptimizationFrameworks/opf_
data…
"time_build" => 0.234455
"time_solve_compilation" => 2.61176
"time_solve" => 0.0376862
"feasible" => truemodel, res = solve_opf_nonconvex(dataset);
resError: MethodError: no method matching Float64(::ForwardDiff.Dual{ForwardDi
ff.Tag{NonconvexUtils.var"#101#108"{NonconvexUtils.var"#100#107"{NonconvexU
tils.var"#97#104"{NonconvexCore.VectorOfFunctions{Tuple{NonconvexCore.IneqC
onstraint{NonconvexCore.var"#80#82"{NonconvexCore.FunctionWrapper{Main.var"
##WeaveSandBox#225".var"#215#257"{Main.var"##WeaveSandBox#225".var"#const_t
hermal_limit#250"{Dict{Int64, Float64}}, Int64, Int64, Int64}}, Differentia
bleFlatten.Unflatten{Tuple{OrderedCollections.OrderedDict{String, Float64}}
, DifferentiableFlatten.var"#unflatten_to_Tuple#11"{Tuple{Int64}, Tuple{Int
64}, Tuple{DifferentiableFlatten.var"#unflatten_to_Dict#16"{typeof(identity
), OrderedCollections.OrderedDict{String, Float64}}}}}}, Float64}, Nonconve
xCore.IneqConstraint{NonconvexCore.var"#80#82"{NonconvexCore.FunctionWrappe
r{Main.var"##WeaveSandBox#225".var"#216#258"{Main.var"##WeaveSandBox#225".v
ar"#const_thermal_limit#250"{Dict{Int64, Float64}}, Int64, Int64, Int64}},
DifferentiableFlatten.Unflatten{Tuple{OrderedCollections.OrderedDict{String
, Float64}}, DifferentiableFlatten.var"#unflatten_to_Tuple#11"{Tuple{Int64}
, Tuple{Int64}, Tuple{DifferentiableFlatten.var"#unflatten_to_Dict#16"{type
of(identity), OrderedCollections.OrderedDict{String, Float64}}}}}}, Float64
}, NonconvexCore.IneqConstraint{NonconvexCore.var"#80#82"{NonconvexCore.Fun
ctionWrapper{Main.var"##WeaveSandBox#225".var"#217#259"{Main.var"##WeaveSan
dBox#225".var"#const_voltage_angle_difference_lb#251"{Dict{Int64, Float64}}
, Int64, Int64, Int64}}, DifferentiableFlatten.Unflatten{Tuple{OrderedColle
ctions.OrderedDict{String, Float64}}, DifferentiableFlatten.var"#unflatten_
to_Tuple#11"{Tuple{Int64}, Tuple{Int64}, Tuple{DifferentiableFlatten.var"#u
nflatten_to_Dict#16"{typeof(identity), OrderedCollections.OrderedDict{Strin
g, Float64}}}}}}, Float64}, NonconvexCore.IneqConstraint{NonconvexCore.var"
#80#82"{NonconvexCore.FunctionWrapper{Main.var"##WeaveSandBox#225".var"#218
#260"{Main.var"##WeaveSandBox#225".var"#const_voltage_angle_difference_ub#2
52"{Dict{Int64, Float64}}, Int64, Int64, Int64}}, DifferentiableFlatten.Unf
latten{Tuple{OrderedCollections.OrderedDict{String, Float64}}, Differentiab
leFlatten.var"#unflatten_to_Tuple#11"{Tuple{Int64}, Tuple{Int64}, Tuple{Dif
ferentiableFlatten.var"#unflatten_to_Dict#16"{typeof(identity), OrderedColl
ections.OrderedDict{String, Float64}}}}}}, Float64}, NonconvexCore.IneqCons
traint{NonconvexCore.var"#80#82"{NonconvexCore.FunctionWrapper{Main.var"##W
eaveSandBox#225".var"#215#257"{Main.var"##WeaveSandBox#225".var"#const_ther
mal_limit#250"{Dict{Int64, Float64}}, Int64, Int64, Int64}}, Differentiable
Flatten.Unflatten{Tuple{OrderedCollections.OrderedDict{String, Float64}}, D
ifferentiableFlatten.var"#unflatten_to_Tuple#11"{Tuple{Int64}, Tuple{Int64}
, Tuple{DifferentiableFlatten.var"#unflatten_to_Dict#16"{typeof(identity),
OrderedCollections.OrderedDict{String, Float64}}}}}}, Float64}, NonconvexCo
re.IneqConstraint{NonconvexCore.var"#80#82"{NonconvexCore.FunctionWrapper{M
ain.var"##WeaveSandBox#225".var"#216#258"{Main.var"##WeaveSandBox#225".var"
#const_thermal_limit#250"{Dict{Int64, Float64}}, Int64, Int64, Int64}}, Dif
ferentiableFlatten.Unflatten{Tuple{OrderedCollections.OrderedDict{String, F
loat64}}, DifferentiableFlatten.var"#unflatten_to_Tuple#11"{Tuple{Int64}, T
uple{Int64}, Tuple{DifferentiableFlatten.var"#unflatten_to_Dict#16"{typeof(
identity), OrderedCollections.OrderedDict{String, Float64}}}}}}, Float64},
NonconvexCore.IneqConstraint{NonconvexCore.var"#80#82"{NonconvexCore.Functi
onWrapper{Main.var"##WeaveSandBox#225".var"#217#259"{Main.var"##WeaveSandBo
x#225".var"#const_voltage_angle_difference_lb#251"{Dict{Int64, Float64}}, I
nt64, Int64, Int64}}, DifferentiableFlatten.Unflatten{Tuple{OrderedCollecti
ons.OrderedDict{String, Float64}}, DifferentiableFlatten.var"#unflatten_to_
Tuple#11"{Tuple{Int64}, Tuple{Int64}, Tuple{DifferentiableFlatten.var"#unfl
atten_to_Dict#16"{typeof(identity), OrderedCollections.OrderedDict{String,
Float64}}}}}}, Float64}, NonconvexCore.IneqConstraint{NonconvexCore.var"#80
#82"{NonconvexCore.FunctionWrapper{Main.var"##WeaveSandBox#225".var"#218#26
0"{Main.var"##WeaveSandBox#225".var"#const_voltage_angle_difference_ub#252"
{Dict{Int64, Float64}}, Int64, Int64, Int64}}, DifferentiableFlatten.Unflat
ten{Tuple{OrderedCollections.OrderedDict{String, Float64}}, DifferentiableF
latten.var"#unflatten_to_Tuple#11"{Tuple{Int64}, Tuple{Int64}, Tuple{Differ
entiableFlatten.var"#unflatten_to_Dict#16"{typeof(identity), OrderedCollect
ions.OrderedDict{String, Float64}}}}}}, Float64}, NonconvexCore.IneqConstra
int{NonconvexCore.var"#80#82"{NonconvexCore.FunctionWrapper{Main.var"##Weav
eSandBox#225".var"#215#257"{Main.var"##WeaveSandBox#225".var"#const_thermal
_limit#250"{Dict{Int64, Float64}}, Int64, Int64, Int64}}, DifferentiableFla
tten.Unflatten{Tuple{OrderedCollections.OrderedDict{String, Float64}}, Diff
erentiableFlatten.var"#unflatten_to_Tuple#11"{Tuple{Int64}, Tuple{Int64}, T
uple{DifferentiableFlatten.var"#unflatten_to_Dict#16"{typeof(identity), Ord
eredCollections.OrderedDict{String, Float64}}}}}}, Float64}, NonconvexCore.
IneqConstraint{NonconvexCore.var"#80#82"{NonconvexCore.FunctionWrapper{Main
.var"##WeaveSandBox#225".var"#216#258"{Main.var"##WeaveSandBox#225".var"#co
nst_thermal_limit#250"{Dict{Int64, Float64}}, Int64, Int64, Int64}}, Differ
entiableFlatten.Unflatten{Tuple{OrderedCollections.OrderedDict{String, Floa
t64}}, DifferentiableFlatten.var"#unflatten_to_Tuple#11"{Tuple{Int64}, Tupl
e{Int64}, Tuple{DifferentiableFlatten.var"#unflatten_to_Dict#16"{typeof(ide
ntity), OrderedCollections.OrderedDict{String, Float64}}}}}}, Float64}, Non
convexCore.IneqConstraint{NonconvexCore.var"#80#82"{NonconvexCore.FunctionW
rapper{Main.var"##WeaveSandBox#225".var"#217#259"{Main.var"##WeaveSandBox#2
25".var"#const_voltage_angle_difference_lb#251"{Dict{Int64, Float64}}, Int6
4, Int64, Int64}}, DifferentiableFlatten.Unflatten{Tuple{OrderedCollections
.OrderedDict{String, Float64}}, DifferentiableFlatten.var"#unflatten_to_Tup
le#11"{Tuple{Int64}, Tuple{Int64}, Tuple{DifferentiableFlatten.var"#unflatt
en_to_Dict#16"{typeof(identity), OrderedCollections.OrderedDict{String, Flo
at64}}}}}}, Float64}, NonconvexCore.IneqConstraint{NonconvexCore.var"#80#82
"{NonconvexCore.FunctionWrapper{Main.var"##WeaveSandBox#225".var"#218#260"{
Main.var"##WeaveSandBox#225".var"#const_voltage_angle_difference_ub#252"{Di
ct{Int64, Float64}}, Int64, Int64, Int64}}, DifferentiableFlatten.Unflatten
{Tuple{OrderedCollections.OrderedDict{String, Float64}}, DifferentiableFlat
ten.var"#unflatten_to_Tuple#11"{Tuple{Int64}, Tuple{Int64}, Tuple{Different
iableFlatten.var"#unflatten_to_Dict#16"{typeof(identity), OrderedCollection
s.OrderedDict{String, Float64}}}}}}, Float64}, NonconvexCore.IneqConstraint
{NonconvexCore.var"#80#82"{NonconvexCore.FunctionWrapper{Main.var"##WeaveSa
ndBox#225".var"#215#257"{Main.var"##WeaveSandBox#225".var"#const_thermal_li
mit#250"{Dict{Int64, Float64}}, Int64, Int64, Int64}}, DifferentiableFlatte
n.Unflatten{Tuple{OrderedCollections.OrderedDict{String, Float64}}, Differe
ntiableFlatten.var"#unflatten_to_Tuple#11"{Tuple{Int64}, Tuple{Int64}, Tupl
e{DifferentiableFlatten.var"#unflatten_to_Dict#16"{typeof(identity), Ordere
dCollections.OrderedDict{String, Float64}}}}}}, Float64}, NonconvexCore.Ine
qConstraint{NonconvexCore.var"#80#82"{NonconvexCore.FunctionWrapper{Main.va
r"##WeaveSandBox#225".var"#216#258"{Main.var"##WeaveSandBox#225".var"#const
_thermal_limit#250"{Dict{Int64, Float64}}, Int64, Int64, Int64}}, Different
iableFlatten.Unflatten{Tuple{OrderedCollections.OrderedDict{String, Float64
}}, DifferentiableFlatten.var"#unflatten_to_Tuple#11"{Tuple{Int64}, Tuple{I
nt64}, Tuple{DifferentiableFlatten.var"#unflatten_to_Dict#16"{typeof(identi
ty), OrderedCollections.OrderedDict{String, Float64}}}}}}, Float64}, Noncon
vexCore.IneqConstraint{NonconvexCore.var"#80#82"{NonconvexCore.FunctionWrap
per{Main.var"##WeaveSandBox#225".var"#217#259"{Main.var"##WeaveSandBox#225"
.var"#const_voltage_angle_difference_lb#251"{Dict{Int64, Float64}}, Int64,
Int64, Int64}}, DifferentiableFlatten.Unflatten{Tuple{OrderedCollections.Or
deredDict{String, Float64}}, DifferentiableFlatten.var"#unflatten_to_Tuple#
11"{Tuple{Int64}, Tuple{Int64}, Tuple{DifferentiableFlatten.var"#unflatten_
to_Dict#16"{typeof(identity), OrderedCollections.OrderedDict{String, Float6
4}}}}}}, Float64}, NonconvexCore.IneqConstraint{NonconvexCore.var"#80#82"{N
onconvexCore.FunctionWrapper{Main.var"##WeaveSandBox#225".var"#218#260"{Mai
n.var"##WeaveSandBox#225".var"#const_voltage_angle_difference_ub#252"{Dict{
Int64, Float64}}, Int64, Int64, Int64}}, DifferentiableFlatten.Unflatten{Tu
ple{OrderedCollections.OrderedDict{String, Float64}}, DifferentiableFlatten
.var"#unflatten_to_Tuple#11"{Tuple{Int64}, Tuple{Int64}, Tuple{Differentiab
leFlatten.var"#unflatten_to_Dict#16"{typeof(identity), OrderedCollections.O
rderedDict{String, Float64}}}}}}, Float64}, NonconvexCore.IneqConstraint{No
nconvexCore.var"#80#82"{NonconvexCore.FunctionWrapper{Main.var"##WeaveSandB
ox#225".var"#215#257"{Main.var"##WeaveSandBox#225".var"#const_thermal_limit
#250"{Dict{Int64, Float64}}, Int64, Int64, Int64}}, DifferentiableFlatten.U
nflatten{Tuple{OrderedCollections.OrderedDict{String, Float64}}, Differenti
ableFlatten.var"#unflatten_to_Tuple#11"{Tuple{Int64}, Tuple{Int64}, Tuple{D
ifferentiableFlatten.var"#unflatten_to_Dict#16"{typeof(identity), OrderedCo
llections.OrderedDict{String, Float64}}}}}}, Float64}, NonconvexCore.IneqCo
nstraint{NonconvexCore.var"#80#82"{NonconvexCore.FunctionWrapper{Main.var"#
#WeaveSandBox#225".var"#216#258"{Main.var"##WeaveSandBox#225".var"#const_th
ermal_limit#250"{Dict{Int64, Float64}}, Int64, Int64, Int64}}, Differentiab
leFlatten.Unflatten{Tuple{OrderedCollections.OrderedDict{String, Float64}},
DifferentiableFlatten.var"#unflatten_to_Tuple#11"{Tuple{Int64}, Tuple{Int6
4}, Tuple{DifferentiableFlatten.var"#unflatten_to_Dict#16"{typeof(identity)
, OrderedCollections.OrderedDict{String, Float64}}}}}}, Float64}, Nonconvex
Core.IneqConstraint{NonconvexCore.var"#80#82"{NonconvexCore.FunctionWrapper
{Main.var"##WeaveSandBox#225".var"#217#259"{Main.var"##WeaveSandBox#225".va
r"#const_voltage_angle_difference_lb#251"{Dict{Int64, Float64}}, Int64, Int
64, Int64}}, DifferentiableFlatten.Unflatten{Tuple{OrderedCollections.Order
edDict{String, Float64}}, DifferentiableFlatten.var"#unflatten_to_Tuple#11"
{Tuple{Int64}, Tuple{Int64}, Tuple{DifferentiableFlatten.var"#unflatten_to_
Dict#16"{typeof(identity), OrderedCollections.OrderedDict{String, Float64}}
}}}}, Float64}, NonconvexCore.IneqConstraint{NonconvexCore.var"#80#82"{Nonc
onvexCore.FunctionWrapper{Main.var"##WeaveSandBox#225".var"#218#260"{Main.v
ar"##WeaveSandBox#225".var"#const_voltage_angle_difference_ub#252"{Dict{Int
64, Float64}}, Int64, Int64, Int64}}, DifferentiableFlatten.Unflatten{Tuple
{OrderedCollections.OrderedDict{String, Float64}}, DifferentiableFlatten.va
r"#unflatten_to_Tuple#11"{Tuple{Int64}, Tuple{Int64}, Tuple{DifferentiableF
latten.var"#unflatten_to_Dict#16"{typeof(identity), OrderedCollections.Orde
redDict{String, Float64}}}}}}, Float64}, NonconvexCore.IneqConstraint{Nonco
nvexCore.var"#80#82"{NonconvexCore.FunctionWrapper{Main.var"##WeaveSandBox#
225".var"#215#257"{Main.var"##WeaveSandBox#225".var"#const_thermal_limit#25
0"{Dict{Int64, Float64}}, Int64, Int64, Int64}}, DifferentiableFlatten.Unfl
atten{Tuple{OrderedCollections.OrderedDict{String, Float64}}, Differentiabl
eFlatten.var"#unflatten_to_Tuple#11"{Tuple{Int64}, Tuple{Int64}, Tuple{Diff
erentiableFlatten.var"#unflatten_to_Dict#16"{typeof(identity), OrderedColle
ctions.OrderedDict{String, Float64}}}}}}, Float64}, NonconvexCore.IneqConst
raint{NonconvexCore.var"#80#82"{NonconvexCore.FunctionWrapper{Main.var"##We
aveSandBox#225".var"#216#258"{Main.var"##WeaveSandBox#225".var"#const_therm
al_limit#250"{Dict{Int64, Float64}}, Int64, Int64, Int64}}, DifferentiableF
latten.Unflatten{Tuple{OrderedCollections.OrderedDict{String, Float64}}, Di
fferentiableFlatten.var"#unflatten_to_Tuple#11"{Tuple{Int64}, Tuple{Int64},
Tuple{DifferentiableFlatten.var"#unflatten_to_Dict#16"{typeof(identity), O
rderedCollections.OrderedDict{String, Float64}}}}}}, Float64}, NonconvexCor
e.IneqConstraint{NonconvexCore.var"#80#82"{NonconvexCore.FunctionWrapper{Ma
in.var"##WeaveSandBox#225".var"#217#259"{Main.var"##WeaveSandBox#225".var"#
const_voltage_angle_difference_lb#251"{Dict{Int64, Float64}}, Int64, Int64,
Int64}}, DifferentiableFlatten.Unflatten{Tuple{OrderedCollections.OrderedD
ict{String, Float64}}, DifferentiableFlatten.var"#unflatten_to_Tuple#11"{Tu
ple{Int64}, Tuple{Int64}, Tuple{DifferentiableFlatten.var"#unflatten_to_Dic
t#16"{typeof(identity), OrderedCollections.OrderedDict{String, Float64}}}}}
}, Float64}, NonconvexCore.IneqConstraint{NonconvexCore.var"#80#82"{Nonconv
exCore.FunctionWrapper{Main.var"##WeaveSandBox#225".var"#218#260"{Main.var"
##WeaveSandBox#225".var"#const_voltage_angle_difference_ub#252"{Dict{Int64,
Float64}}, Int64, Int64, Int64}}, DifferentiableFlatten.Unflatten{Tuple{Or
deredCollections.OrderedDict{String, Float64}}, DifferentiableFlatten.var"#
unflatten_to_Tuple#11"{Tuple{Int64}, Tuple{Int64}, Tuple{DifferentiableFlat
ten.var"#unflatten_to_Dict#16"{typeof(identity), OrderedCollections.Ordered
Dict{String, Float64}}}}}}, Float64}}}}}}, Float64}, Float64, 1})
Closest candidates are:
(::Type{T})(::Real, !Matched::RoundingMode) where T<:AbstractFloat
@ Base rounding.jl:207
(::Type{T})(::T) where T<:Number
@ Core boot.jl:792
Float64(!Matched::Int16)
@ Base float.jl:159
...model, res = solve_opf_optim(dataset);
resDict{String, Any} with 8 entries:
"cost" => 77.9548
"variables" => 44
"constraints" => 53
"case" => "../../benchmarks/OptimizationFrameworks/opf_
data…
"time_build" => 0.000490626
"time_solve_compilation" => 21.6869
"time_solve" => 17.3996
"feasible" => falsefile_name = "../../benchmarks/OptimizationFrameworks/opf_data/pglib_opf_case3_lmbd.m"
dataset = load_and_setup_data(file_name);Error: IOError: stream is closed or unusablemodel, res = solve_opf_optimization(dataset);
resDict{String, Any} with 8 entries:
"cost" => 17551.9
"variables" => 44
"constraints" => 53
"case" => "../../benchmarks/OptimizationFrameworks/opf_
data…
"time_build" => 0.178576
"time_solve_compilation" => 0.0861387
"time_solve" => 0.084448
"feasible" => truemodel, res = solve_opf_jump(dataset);
resDict{String, Any} with 8 entries:
"cost" => 17551.9
"variables" => 44
"constraints" => 53
"case" => "../../benchmarks/OptimizationFrameworks/opf_
data…
"time_build" => 0.00204671
"time_solve_compilation" => 0.0165238
"time_solve" => 0.0159231
"feasible" => truemodel, res = solve_opf_nlpmodels(dataset);
resDict{String, Any} with 8 entries:
"cost" => 17551.9
"variables" => 44
"constraints" => 53
"case" => "../../benchmarks/OptimizationFrameworks/opf_
data…
"time_build" => 0.0252522
"time_solve_compilation" => 0.0391716
"time_solve" => 0.0384374
"feasible" => truemodel, res = solve_opf_nonconvex(dataset);
resError: MethodError: no method matching Float64(::ForwardDiff.Dual{ForwardDi
ff.Tag{NonconvexUtils.var"#101#108"{NonconvexUtils.var"#100#107"{NonconvexU
tils.var"#97#104"{NonconvexCore.VectorOfFunctions{Tuple{NonconvexCore.IneqC
onstraint{NonconvexCore.var"#80#82"{NonconvexCore.FunctionWrapper{Main.var"
##WeaveSandBox#225".var"#215#257"{Main.var"##WeaveSandBox#225".var"#const_t
hermal_limit#250"{Dict{Int64, Float64}}, Int64, Int64, Int64}}, Differentia
bleFlatten.Unflatten{Tuple{OrderedCollections.OrderedDict{String, Float64}}
, DifferentiableFlatten.var"#unflatten_to_Tuple#11"{Tuple{Int64}, Tuple{Int
64}, Tuple{DifferentiableFlatten.var"#unflatten_to_Dict#16"{typeof(identity
), OrderedCollections.OrderedDict{String, Float64}}}}}}, Float64}, Nonconve
xCore.IneqConstraint{NonconvexCore.var"#80#82"{NonconvexCore.FunctionWrappe
r{Main.var"##WeaveSandBox#225".var"#216#258"{Main.var"##WeaveSandBox#225".v
ar"#const_thermal_limit#250"{Dict{Int64, Float64}}, Int64, Int64, Int64}},
DifferentiableFlatten.Unflatten{Tuple{OrderedCollections.OrderedDict{String
, Float64}}, DifferentiableFlatten.var"#unflatten_to_Tuple#11"{Tuple{Int64}
, Tuple{Int64}, Tuple{DifferentiableFlatten.var"#unflatten_to_Dict#16"{type
of(identity), OrderedCollections.OrderedDict{String, Float64}}}}}}, Float64
}, NonconvexCore.IneqConstraint{NonconvexCore.var"#80#82"{NonconvexCore.Fun
ctionWrapper{Main.var"##WeaveSandBox#225".var"#217#259"{Main.var"##WeaveSan
dBox#225".var"#const_voltage_angle_difference_lb#251"{Dict{Int64, Float64}}
, Int64, Int64, Int64}}, DifferentiableFlatten.Unflatten{Tuple{OrderedColle
ctions.OrderedDict{String, Float64}}, DifferentiableFlatten.var"#unflatten_
to_Tuple#11"{Tuple{Int64}, Tuple{Int64}, Tuple{DifferentiableFlatten.var"#u
nflatten_to_Dict#16"{typeof(identity), OrderedCollections.OrderedDict{Strin
g, Float64}}}}}}, Float64}, NonconvexCore.IneqConstraint{NonconvexCore.var"
#80#82"{NonconvexCore.FunctionWrapper{Main.var"##WeaveSandBox#225".var"#218
#260"{Main.var"##WeaveSandBox#225".var"#const_voltage_angle_difference_ub#2
52"{Dict{Int64, Float64}}, Int64, Int64, Int64}}, DifferentiableFlatten.Unf
latten{Tuple{OrderedCollections.OrderedDict{String, Float64}}, Differentiab
leFlatten.var"#unflatten_to_Tuple#11"{Tuple{Int64}, Tuple{Int64}, Tuple{Dif
ferentiableFlatten.var"#unflatten_to_Dict#16"{typeof(identity), OrderedColl
ections.OrderedDict{String, Float64}}}}}}, Float64}, NonconvexCore.IneqCons
traint{NonconvexCore.var"#80#82"{NonconvexCore.FunctionWrapper{Main.var"##W
eaveSandBox#225".var"#215#257"{Main.var"##WeaveSandBox#225".var"#const_ther
mal_limit#250"{Dict{Int64, Float64}}, Int64, Int64, Int64}}, Differentiable
Flatten.Unflatten{Tuple{OrderedCollections.OrderedDict{String, Float64}}, D
ifferentiableFlatten.var"#unflatten_to_Tuple#11"{Tuple{Int64}, Tuple{Int64}
, Tuple{DifferentiableFlatten.var"#unflatten_to_Dict#16"{typeof(identity),
OrderedCollections.OrderedDict{String, Float64}}}}}}, Float64}, NonconvexCo
re.IneqConstraint{NonconvexCore.var"#80#82"{NonconvexCore.FunctionWrapper{M
ain.var"##WeaveSandBox#225".var"#216#258"{Main.var"##WeaveSandBox#225".var"
#const_thermal_limit#250"{Dict{Int64, Float64}}, Int64, Int64, Int64}}, Dif
ferentiableFlatten.Unflatten{Tuple{OrderedCollections.OrderedDict{String, F
loat64}}, DifferentiableFlatten.var"#unflatten_to_Tuple#11"{Tuple{Int64}, T
uple{Int64}, Tuple{DifferentiableFlatten.var"#unflatten_to_Dict#16"{typeof(
identity), OrderedCollections.OrderedDict{String, Float64}}}}}}, Float64},
NonconvexCore.IneqConstraint{NonconvexCore.var"#80#82"{NonconvexCore.Functi
onWrapper{Main.var"##WeaveSandBox#225".var"#217#259"{Main.var"##WeaveSandBo
x#225".var"#const_voltage_angle_difference_lb#251"{Dict{Int64, Float64}}, I
nt64, Int64, Int64}}, DifferentiableFlatten.Unflatten{Tuple{OrderedCollecti
ons.OrderedDict{String, Float64}}, DifferentiableFlatten.var"#unflatten_to_
Tuple#11"{Tuple{Int64}, Tuple{Int64}, Tuple{DifferentiableFlatten.var"#unfl
atten_to_Dict#16"{typeof(identity), OrderedCollections.OrderedDict{String,
Float64}}}}}}, Float64}, NonconvexCore.IneqConstraint{NonconvexCore.var"#80
#82"{NonconvexCore.FunctionWrapper{Main.var"##WeaveSandBox#225".var"#218#26
0"{Main.var"##WeaveSandBox#225".var"#const_voltage_angle_difference_ub#252"
{Dict{Int64, Float64}}, Int64, Int64, Int64}}, DifferentiableFlatten.Unflat
ten{Tuple{OrderedCollections.OrderedDict{String, Float64}}, DifferentiableF
latten.var"#unflatten_to_Tuple#11"{Tuple{Int64}, Tuple{Int64}, Tuple{Differ
entiableFlatten.var"#unflatten_to_Dict#16"{typeof(identity), OrderedCollect
ions.OrderedDict{String, Float64}}}}}}, Float64}, NonconvexCore.IneqConstra
int{NonconvexCore.var"#80#82"{NonconvexCore.FunctionWrapper{Main.var"##Weav
eSandBox#225".var"#215#257"{Main.var"##WeaveSandBox#225".var"#const_thermal
_limit#250"{Dict{Int64, Float64}}, Int64, Int64, Int64}}, DifferentiableFla
tten.Unflatten{Tuple{OrderedCollections.OrderedDict{String, Float64}}, Diff
erentiableFlatten.var"#unflatten_to_Tuple#11"{Tuple{Int64}, Tuple{Int64}, T
uple{DifferentiableFlatten.var"#unflatten_to_Dict#16"{typeof(identity), Ord
eredCollections.OrderedDict{String, Float64}}}}}}, Float64}, NonconvexCore.
IneqConstraint{NonconvexCore.var"#80#82"{NonconvexCore.FunctionWrapper{Main
.var"##WeaveSandBox#225".var"#216#258"{Main.var"##WeaveSandBox#225".var"#co
nst_thermal_limit#250"{Dict{Int64, Float64}}, Int64, Int64, Int64}}, Differ
entiableFlatten.Unflatten{Tuple{OrderedCollections.OrderedDict{String, Floa
t64}}, DifferentiableFlatten.var"#unflatten_to_Tuple#11"{Tuple{Int64}, Tupl
e{Int64}, Tuple{DifferentiableFlatten.var"#unflatten_to_Dict#16"{typeof(ide
ntity), OrderedCollections.OrderedDict{String, Float64}}}}}}, Float64}, Non
convexCore.IneqConstraint{NonconvexCore.var"#80#82"{NonconvexCore.FunctionW
rapper{Main.var"##WeaveSandBox#225".var"#217#259"{Main.var"##WeaveSandBox#2
25".var"#const_voltage_angle_difference_lb#251"{Dict{Int64, Float64}}, Int6
4, Int64, Int64}}, DifferentiableFlatten.Unflatten{Tuple{OrderedCollections
.OrderedDict{String, Float64}}, DifferentiableFlatten.var"#unflatten_to_Tup
le#11"{Tuple{Int64}, Tuple{Int64}, Tuple{DifferentiableFlatten.var"#unflatt
en_to_Dict#16"{typeof(identity), OrderedCollections.OrderedDict{String, Flo
at64}}}}}}, Float64}, NonconvexCore.IneqConstraint{NonconvexCore.var"#80#82
"{NonconvexCore.FunctionWrapper{Main.var"##WeaveSandBox#225".var"#218#260"{
Main.var"##WeaveSandBox#225".var"#const_voltage_angle_difference_ub#252"{Di
ct{Int64, Float64}}, Int64, Int64, Int64}}, DifferentiableFlatten.Unflatten
{Tuple{OrderedCollections.OrderedDict{String, Float64}}, DifferentiableFlat
ten.var"#unflatten_to_Tuple#11"{Tuple{Int64}, Tuple{Int64}, Tuple{Different
iableFlatten.var"#unflatten_to_Dict#16"{typeof(identity), OrderedCollection
s.OrderedDict{String, Float64}}}}}}, Float64}, NonconvexCore.IneqConstraint
{NonconvexCore.var"#80#82"{NonconvexCore.FunctionWrapper{Main.var"##WeaveSa
ndBox#225".var"#215#257"{Main.var"##WeaveSandBox#225".var"#const_thermal_li
mit#250"{Dict{Int64, Float64}}, Int64, Int64, Int64}}, DifferentiableFlatte
n.Unflatten{Tuple{OrderedCollections.OrderedDict{String, Float64}}, Differe
ntiableFlatten.var"#unflatten_to_Tuple#11"{Tuple{Int64}, Tuple{Int64}, Tupl
e{DifferentiableFlatten.var"#unflatten_to_Dict#16"{typeof(identity), Ordere
dCollections.OrderedDict{String, Float64}}}}}}, Float64}, NonconvexCore.Ine
qConstraint{NonconvexCore.var"#80#82"{NonconvexCore.FunctionWrapper{Main.va
r"##WeaveSandBox#225".var"#216#258"{Main.var"##WeaveSandBox#225".var"#const
_thermal_limit#250"{Dict{Int64, Float64}}, Int64, Int64, Int64}}, Different
iableFlatten.Unflatten{Tuple{OrderedCollections.OrderedDict{String, Float64
}}, DifferentiableFlatten.var"#unflatten_to_Tuple#11"{Tuple{Int64}, Tuple{I
nt64}, Tuple{DifferentiableFlatten.var"#unflatten_to_Dict#16"{typeof(identi
ty), OrderedCollections.OrderedDict{String, Float64}}}}}}, Float64}, Noncon
vexCore.IneqConstraint{NonconvexCore.var"#80#82"{NonconvexCore.FunctionWrap
per{Main.var"##WeaveSandBox#225".var"#217#259"{Main.var"##WeaveSandBox#225"
.var"#const_voltage_angle_difference_lb#251"{Dict{Int64, Float64}}, Int64,
Int64, Int64}}, DifferentiableFlatten.Unflatten{Tuple{OrderedCollections.Or
deredDict{String, Float64}}, DifferentiableFlatten.var"#unflatten_to_Tuple#
11"{Tuple{Int64}, Tuple{Int64}, Tuple{DifferentiableFlatten.var"#unflatten_
to_Dict#16"{typeof(identity), OrderedCollections.OrderedDict{String, Float6
4}}}}}}, Float64}, NonconvexCore.IneqConstraint{NonconvexCore.var"#80#82"{N
onconvexCore.FunctionWrapper{Main.var"##WeaveSandBox#225".var"#218#260"{Mai
n.var"##WeaveSandBox#225".var"#const_voltage_angle_difference_ub#252"{Dict{
Int64, Float64}}, Int64, Int64, Int64}}, DifferentiableFlatten.Unflatten{Tu
ple{OrderedCollections.OrderedDict{String, Float64}}, DifferentiableFlatten
.var"#unflatten_to_Tuple#11"{Tuple{Int64}, Tuple{Int64}, Tuple{Differentiab
leFlatten.var"#unflatten_to_Dict#16"{typeof(identity), OrderedCollections.O
rderedDict{String, Float64}}}}}}, Float64}, NonconvexCore.IneqConstraint{No
nconvexCore.var"#80#82"{NonconvexCore.FunctionWrapper{Main.var"##WeaveSandB
ox#225".var"#215#257"{Main.var"##WeaveSandBox#225".var"#const_thermal_limit
#250"{Dict{Int64, Float64}}, Int64, Int64, Int64}}, DifferentiableFlatten.U
nflatten{Tuple{OrderedCollections.OrderedDict{String, Float64}}, Differenti
ableFlatten.var"#unflatten_to_Tuple#11"{Tuple{Int64}, Tuple{Int64}, Tuple{D
ifferentiableFlatten.var"#unflatten_to_Dict#16"{typeof(identity), OrderedCo
llections.OrderedDict{String, Float64}}}}}}, Float64}, NonconvexCore.IneqCo
nstraint{NonconvexCore.var"#80#82"{NonconvexCore.FunctionWrapper{Main.var"#
#WeaveSandBox#225".var"#216#258"{Main.var"##WeaveSandBox#225".var"#const_th
ermal_limit#250"{Dict{Int64, Float64}}, Int64, Int64, Int64}}, Differentiab
leFlatten.Unflatten{Tuple{OrderedCollections.OrderedDict{String, Float64}},
DifferentiableFlatten.var"#unflatten_to_Tuple#11"{Tuple{Int64}, Tuple{Int6
4}, Tuple{DifferentiableFlatten.var"#unflatten_to_Dict#16"{typeof(identity)
, OrderedCollections.OrderedDict{String, Float64}}}}}}, Float64}, Nonconvex
Core.IneqConstraint{NonconvexCore.var"#80#82"{NonconvexCore.FunctionWrapper
{Main.var"##WeaveSandBox#225".var"#217#259"{Main.var"##WeaveSandBox#225".va
r"#const_voltage_angle_difference_lb#251"{Dict{Int64, Float64}}, Int64, Int
64, Int64}}, DifferentiableFlatten.Unflatten{Tuple{OrderedCollections.Order
edDict{String, Float64}}, DifferentiableFlatten.var"#unflatten_to_Tuple#11"
{Tuple{Int64}, Tuple{Int64}, Tuple{DifferentiableFlatten.var"#unflatten_to_
Dict#16"{typeof(identity), OrderedCollections.OrderedDict{String, Float64}}
}}}}, Float64}, NonconvexCore.IneqConstraint{NonconvexCore.var"#80#82"{Nonc
onvexCore.FunctionWrapper{Main.var"##WeaveSandBox#225".var"#218#260"{Main.v
ar"##WeaveSandBox#225".var"#const_voltage_angle_difference_ub#252"{Dict{Int
64, Float64}}, Int64, Int64, Int64}}, DifferentiableFlatten.Unflatten{Tuple
{OrderedCollections.OrderedDict{String, Float64}}, DifferentiableFlatten.va
r"#unflatten_to_Tuple#11"{Tuple{Int64}, Tuple{Int64}, Tuple{DifferentiableF
latten.var"#unflatten_to_Dict#16"{typeof(identity), OrderedCollections.Orde
redDict{String, Float64}}}}}}, Float64}, NonconvexCore.IneqConstraint{Nonco
nvexCore.var"#80#82"{NonconvexCore.FunctionWrapper{Main.var"##WeaveSandBox#
225".var"#215#257"{Main.var"##WeaveSandBox#225".var"#const_thermal_limit#25
0"{Dict{Int64, Float64}}, Int64, Int64, Int64}}, DifferentiableFlatten.Unfl
atten{Tuple{OrderedCollections.OrderedDict{String, Float64}}, Differentiabl
eFlatten.var"#unflatten_to_Tuple#11"{Tuple{Int64}, Tuple{Int64}, Tuple{Diff
erentiableFlatten.var"#unflatten_to_Dict#16"{typeof(identity), OrderedColle
ctions.OrderedDict{String, Float64}}}}}}, Float64}, NonconvexCore.IneqConst
raint{NonconvexCore.var"#80#82"{NonconvexCore.FunctionWrapper{Main.var"##We
aveSandBox#225".var"#216#258"{Main.var"##WeaveSandBox#225".var"#const_therm
al_limit#250"{Dict{Int64, Float64}}, Int64, Int64, Int64}}, DifferentiableF
latten.Unflatten{Tuple{OrderedCollections.OrderedDict{String, Float64}}, Di
fferentiableFlatten.var"#unflatten_to_Tuple#11"{Tuple{Int64}, Tuple{Int64},
Tuple{DifferentiableFlatten.var"#unflatten_to_Dict#16"{typeof(identity), O
rderedCollections.OrderedDict{String, Float64}}}}}}, Float64}, NonconvexCor
e.IneqConstraint{NonconvexCore.var"#80#82"{NonconvexCore.FunctionWrapper{Ma
in.var"##WeaveSandBox#225".var"#217#259"{Main.var"##WeaveSandBox#225".var"#
const_voltage_angle_difference_lb#251"{Dict{Int64, Float64}}, Int64, Int64,
Int64}}, DifferentiableFlatten.Unflatten{Tuple{OrderedCollections.OrderedD
ict{String, Float64}}, DifferentiableFlatten.var"#unflatten_to_Tuple#11"{Tu
ple{Int64}, Tuple{Int64}, Tuple{DifferentiableFlatten.var"#unflatten_to_Dic
t#16"{typeof(identity), OrderedCollections.OrderedDict{String, Float64}}}}}
}, Float64}, NonconvexCore.IneqConstraint{NonconvexCore.var"#80#82"{Nonconv
exCore.FunctionWrapper{Main.var"##WeaveSandBox#225".var"#218#260"{Main.var"
##WeaveSandBox#225".var"#const_voltage_angle_difference_ub#252"{Dict{Int64,
Float64}}, Int64, Int64, Int64}}, DifferentiableFlatten.Unflatten{Tuple{Or
deredCollections.OrderedDict{String, Float64}}, DifferentiableFlatten.var"#
unflatten_to_Tuple#11"{Tuple{Int64}, Tuple{Int64}, Tuple{DifferentiableFlat
ten.var"#unflatten_to_Dict#16"{typeof(identity), OrderedCollections.Ordered
Dict{String, Float64}}}}}}, Float64}}}}}}, Float64}, Float64, 1})
Closest candidates are:
(::Type{T})(::Real, !Matched::RoundingMode) where T<:AbstractFloat
@ Base rounding.jl:207
(::Type{T})(::T) where T<:Number
@ Core boot.jl:792
Float64(!Matched::Int16)
@ Base float.jl:159
...model, res = solve_opf_optim(dataset);
resDict{String, Any} with 8 entries:
"cost" => 77.9548
"variables" => 44
"constraints" => 53
"case" => "../../benchmarks/OptimizationFrameworks/opf_
data…
"time_build" => 0.000575455
"time_solve_compilation" => 17.7061
"time_solve" => 17.6302
"feasible" => falseusing DataFrames, PrettyTables
function multidata_multisolver_benchmark(dataset_files; sizelimit = SIZE_LIMIT)
cases = String[]
vars = Int[]
cons = Int[]
optimization_time = Float64[]
mtk_time = Float64[]
jump_time = Float64[]
nlpmodels_time = Float64[]
nonconvex_time = Float64[]
optim_time = Float64[]
optimization_time_modelbuild = Float64[]
mtk_time_modelbuild = Float64[]
jump_time_modelbuild = Float64[]
nlpmodels_time_modelbuild = Float64[]
nonconvex_time_modelbuild = Float64[]
optim_time_modelbuild = Float64[]
optimization_time_compilation = Float64[]
mtk_time_compilation = Float64[]
jump_time_compilation = Float64[]
nlpmodels_time_compilation = Float64[]
nonconvex_time_compilation = Float64[]
optim_time_compilation = Float64[]
optimization_cost = Float64[]
mtk_cost = Float64[]
jump_cost = Float64[]
nlpmodels_cost = Float64[]
nonconvex_cost = Float64[]
optim_cost = Float64[]
for file in dataset_files
@show file
dataset = load_and_setup_data(file)
prob = build_opf_optimization_prob(dataset)
@info "Number of Variables: $(length(prob.u0))"
@info "Number of Constraints: $(length(prob.lcons))"
if length(prob.u0) > sizelimit
@info "Variable size over global limit. Skipping for now"
continue
end
@info "Running Optimization.jl"
model, res = solve_opf_optimization(dataset)
push!(cases, split(file,"/")[end])
push!(vars, res["variables"])
push!(cons, res["constraints"])
push!(optimization_time, res["time_solve"])
push!(optimization_time_modelbuild, res["time_build"])
push!(optimization_time_compilation, res["time_solve_compilation"])
push!(optimization_cost, res["cost"])
@info "Running ModelingToolkit.jl"
model, res = solve_opf_mtk(dataset)
push!(mtk_time, res["time_solve"])
push!(mtk_time_modelbuild, res["time_build"])
push!(mtk_time_compilation, res["time_solve_compilation"])
push!(mtk_cost, res["cost"])
@info "Running JuMP.jl"
model, res = solve_opf_jump(dataset)
push!(jump_time, res["time_solve"])
push!(jump_time_modelbuild, res["time_build"])
push!(jump_time_compilation, res["time_solve_compilation"])
push!(jump_cost, res["cost"])
@info "Running NLPModels.jl"
model, res = solve_opf_nlpmodels(dataset)
push!(nlpmodels_time, res["time_solve"])
push!(nlpmodels_time_modelbuild, res["time_build"])
push!(nlpmodels_time_compilation, res["time_solve_compilation"])
push!(nlpmodels_cost, res["cost"])
#=
@info "Running Nonconvex.jl"
model, res = solve_opf_nonconvex(dataset)
push!(nonconvex_time, res["time_solve"])
push!(nonconvex_time_modelbuild, res["time_build"])
push!(nonconvex_time_compilation, res["time_solve_compilation"])
push!(nonconvex_cost, res["cost"])
=#
if length(prob.u0) > 400
@info "Running Optim.jl"
model, res = solve_opf_optim(dataset)
push!(optim_time, NaN)
push!(optim_time_modelbuild, NaN)
push!(optim_time_compilation, NaN)
push!(optim_cost, NaN)
else
@info "Running Optim.jl"
model, res = solve_opf_optim(dataset)
push!(optim_time, res["time_solve"])
push!(optim_time_modelbuild, res["time_build"])
push!(optim_time_compilation, res["time_solve_compilation"])
push!(optim_cost, res["cost"])
end
end
DataFrame(:case => cases, :vars => vars, :cons => cons,
:optimization => optimization_time, :optimization_modelbuild => optimization_time_modelbuild, :optimization_wcompilation => optimization_time_compilation, :optimization_cost => optimization_cost,
:mtk => mtk_time, :mtk_time_modelbuild => mtk_time_modelbuild, :mtk_time_wcompilation => mtk_time_compilation, :mtk_cost => mtk_cost,
:jump => jump_time, :jump_modelbuild => jump_time_modelbuild, :jump_wcompilation => jump_time_compilation, :jump_cost => jump_cost,
:nlpmodels => nlpmodels_time, :nlpmodels_modelbuild => nlpmodels_time_modelbuild, :nlpmodels_wcompilation => nlpmodels_time_compilation, :nlpmodels_cost => nlpmodels_cost,
#:nonconvex => nonconvex_time, :nonconvex_modelbuild => nonconvex_time_modelbuild, :nonconvex_wcompilation => nonconvex_time_compilation, :nonconvex_cost => nonconvex_cost,
:optim => optim_time, :optim_modelbuild => optim_time_modelbuild, :optim_wcompilation => optim_time_compilation, :optim_cost => optim_cost)
end
test_datasets = [
"../../benchmarks/OptimizationFrameworks/opf_data/pglib_opf_case3_lmbd.m",
"../../benchmarks/OptimizationFrameworks/opf_data/pglib_opf_case5_pjm.m"
]2-element Vector{String}:
"../../benchmarks/OptimizationFrameworks/opf_data/pglib_opf_case3_lmbd.m"
"../../benchmarks/OptimizationFrameworks/opf_data/pglib_opf_case5_pjm.m"timing_data = multidata_multisolver_benchmark(test_datasets)file = "../../benchmarks/OptimizationFrameworks/opf_data/pglib_opf_case3_lmbd.m"
Error: IOError: stream is closed or unusableio = IOBuffer()
println(io, "```@raw html")
pretty_table(io, timing_data; backend = Val(:html))
# show(io, "text/html", pretty_table(timing_data; backend = Val(:html)))
println(io, "```")
Text(String(take!(io)))Error: UndefVarError: timing_data not defined
Run the Full Benchmark
using LibGit2
tmpdir = Base.Filesystem.mktempdir()
LibGit2.clone("https://github.com/power-grid-lib/pglib-opf", tmpdir)
benchmarkfiles = readdir(tmpdir)
benchmarkfiles = benchmarkfiles[endswith(".m").(benchmarkfiles)]
benchmark_datasets = joinpath.((tmpdir,),benchmarkfiles)66-element Vector{String}:
"/tmp/jl_BZO4gN/pglib_opf_case10000_goc.m"
"/tmp/jl_BZO4gN/pglib_opf_case10192_epigrids.m"
"/tmp/jl_BZO4gN/pglib_opf_case10480_goc.m"
"/tmp/jl_BZO4gN/pglib_opf_case118_ieee.m"
"/tmp/jl_BZO4gN/pglib_opf_case1354_pegase.m"
"/tmp/jl_BZO4gN/pglib_opf_case13659_pegase.m"
"/tmp/jl_BZO4gN/pglib_opf_case14_ieee.m"
"/tmp/jl_BZO4gN/pglib_opf_case162_ieee_dtc.m"
"/tmp/jl_BZO4gN/pglib_opf_case179_goc.m"
"/tmp/jl_BZO4gN/pglib_opf_case1803_snem.m"
⋮
"/tmp/jl_BZO4gN/pglib_opf_case6515_rte.m"
"/tmp/jl_BZO4gN/pglib_opf_case7336_epigrids.m"
"/tmp/jl_BZO4gN/pglib_opf_case73_ieee_rts.m"
"/tmp/jl_BZO4gN/pglib_opf_case78484_epigrids.m"
"/tmp/jl_BZO4gN/pglib_opf_case793_goc.m"
"/tmp/jl_BZO4gN/pglib_opf_case8387_pegase.m"
"/tmp/jl_BZO4gN/pglib_opf_case89_pegase.m"
"/tmp/jl_BZO4gN/pglib_opf_case9241_pegase.m"
"/tmp/jl_BZO4gN/pglib_opf_case9591_goc.m"timing_data = multidata_multisolver_benchmark(benchmark_datasets)file = "/tmp/jl_BZO4gN/pglib_opf_case10000_goc.m"
Error: IOError: stream is closed or unusableio = IOBuffer()
println(io, "```@raw html")
pretty_table(io, timing_data; backend = Val(:html))
# show(io, "text/html", pretty_table(timing_data; backend = Val(:html)))
println(io, "```")
Text(String(take!(io)))Error: UndefVarError: timing_data not defined
Appendix
Error: ArgumentError: Package SciMLBenchmarks not found in current path, ma
ybe you meant `import/using ..SciMLBenchmarks`.
- Otherwise, run `import Pkg; Pkg.add("SciMLBenchmarks")` to install the Sc
iMLBenchmarks package.