diff --git a/.gitignore b/.gitignore index 29126e4..7cd2812 100644 --- a/.gitignore +++ b/.gitignore @@ -22,3 +22,5 @@ docs/site/ # committed for packages, but should be committed for applications that require a static # environment. Manifest.toml + +data/*.jld2 diff --git a/Project.toml b/Project.toml index e9a854d..7e83bc9 100644 --- a/Project.toml +++ b/Project.toml @@ -3,35 +3,37 @@ uuid = "2fff4b78-0b6c-428d-bac8-85ccea8c4bdf" version = "0.1.0" [deps] +CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab" -PGLib = "07a8691f-3d11-4330-951b-3c50f98338be" -PowerModels = "c36e90e8-916a-50a6-bd94-075b64ef4655" +Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6" ExaModels = "1037b233-b668-4ce9-9b63-f9f681f55dd2" +ExaPowerIO = "14903efe-9500-4d7f-a589-7ab7e15da6de" +GOC3Benchmark = "3a45b339-860d-44d0-a64b-5f943cdd120b" JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819" JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6" -GOC3Benchmark = "3a45b339-860d-44d0-a64b-5f943cdd120b" -Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6" + +[sources] +ExaPowerIO = {rev = "003970f", url = "https://github.com/MadNLP/ExaPowerIO.jl"} +GOC3Benchmark = {rev = "588f356", url = "https://github.com/lanl-ansi/GOC3Benchmark.jl"} [compat] -ExaModels = "~0.9" -PowerModels = "~0.21" -JLD2 = "~0.5" CUDA = "5 - 5.8.2" +ExaModels = "0.9.1" +ExaPowerIO = "0.2.0" +JLD2 = "~0.5" MadNLP = "~0.8" MadNLPGPU = "~0.7" [extras] CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" +Ipopt = "b6b21f68-93f8-5de0-b562-5493be1d77c9" +JuMP = "4076af6c-e467-56ae-b986-b466b2749572" KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c" MadNLP = "2621e9c9-9eb4-46b1-8089-e8c72242dfb6" MadNLPGPU = "d72a61cc-809d-412f-99be-fd81f4b8a598" -Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" -Ipopt = "b6b21f68-93f8-5de0-b562-5493be1d77c9" -JuMP = "4076af6c-e467-56ae-b986-b466b2749572" NLPModelsJuMP = "792afdf1-32c1-5681-94e0-d7bf7a5df49e" +PowerModels = "c36e90e8-916a-50a6-bd94-075b64ef4655" +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["Test", "MadNLP", "MadNLPGPU", "KernelAbstractions", "CUDA", "Ipopt", "JuMP", "NLPModelsJuMP"] - -[sources] -GOC3Benchmark = {url = "https://github.com/lanl-ansi/GOC3Benchmark.jl", rev = "588f356"} +test = ["Test", "MadNLP", "MadNLPGPU", "KernelAbstractions", "CUDA", "Ipopt", "JuMP", "NLPModelsJuMP", "PowerModels"] diff --git a/docs/Project.toml b/docs/Project.toml index 2d90f1a..95e6a79 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -10,6 +10,7 @@ ExaModelsPower = "2fff4b78-0b6c-428d-bac8-85ccea8c4bdf" [compat] CUDA = "5 - 5.8.2" +ExaModels = "0.9.1" [sources] ExaModelsPower = {path = ".."} diff --git a/src/ExaModelsPower.jl b/src/ExaModelsPower.jl index eaecd71..a072f8a 100644 --- a/src/ExaModelsPower.jl +++ b/src/ExaModelsPower.jl @@ -3,8 +3,6 @@ module ExaModelsPower import JLD2 import Downloads import ExaModels: ExaCore, variable, constraint, ExaModel, objective, constraint!, convert_array, solution -import PGLib -import PowerModels include("parser.jl") @@ -32,7 +30,6 @@ function __init__() global TMPDIR = joinpath(@__DIR__,"..","data") mkpath(TMPDIR) end - PowerModels.silence() end end # module ExaModelsExamples diff --git a/src/constraint.jl b/src/constraint.jl index 4cae701..a90b6c6 100644 --- a/src/constraint.jl +++ b/src/constraint.jl @@ -1,5 +1,6 @@ +# cost1 function gen_cost(g, pg) - return g.cost1 * pg^2 + g.cost2 * pg + g.cost3 + return g.c[1] * pg^2 + g.c[2] * pg + g.c[3] end function c_ref_angle_polar(va) @@ -46,7 +47,7 @@ end #no coordinates specified function c_thermal_limit(b, p,q) - return p^2 + q^2 - b.rate_a_sq + return p^2 + q^2 - b.rate_a^2 end #only for mp @@ -125,15 +126,15 @@ function c_ohms_rect(pst, qst, vr, vim, I2) end function c_stor_state(s, E0, E1, pstc, pstd) - return E0 - E1 - (s.etac*pstc - pstd/s.etad) + return E0 - E1 - (s.charge_efficiency*pstc - pstd/s.discharge_efficiency) end function c_storage_state_smooth(s, E0, E1, discharge_func::Function, pstd) - return E0 - E1 + discharge_func(pstd, s.Srating) + return E0 - E1 + discharge_func(pstd, s.thermal_rating) end function c_transfer_lim(s, pst, qst) - return pst^2 + qst^2 - s.Srating^2 + return pst^2 + qst^2 - s.thermal_rating^2 end #used for charge and discharge limits @@ -147,4 +148,4 @@ end function c_comp(pstd, pstc) return pstd*pstc -end \ No newline at end of file +end diff --git a/src/mpopf.jl b/src/mpopf.jl index f525469..4675804 100644 --- a/src/mpopf.jl +++ b/src/mpopf.jl @@ -2,7 +2,7 @@ using DelimitedFiles function parse_mp_power_data(filename, N, corrective_action_ratio) - data, dicts = parse_ac_power_data(filename) + data = parse_ac_power_data(filename) nbus = length(data.bus) @@ -12,49 +12,67 @@ function parse_mp_power_data(filename, N, corrective_action_ratio) ; data..., refarray = [(i,t) for i in data.ref_buses, t in 1:N], - barray = [(;b..., t = t) for b in data.branch, t in 1:N ], - busarray = [(;b..., t = t) for b in data.bus, t in 1:N ], - arcarray = [(;a..., t = t) for a in data.arc, t in 1:N ], - genarray = [(;g..., t = t) for g in data.gen, t in 1:N ], - storarray = isempty(data.storage) ? empty_data = empty_stor : [(;s..., t = t) for s in data.storage, t in 1:N], + barray = [(;b, t = t) for b in data.branch, t in 1:N ], + busarray = [(;b, t = t) for b in data.bus, t in 1:N ], + arcarray = [(;a, t = t) for a in data.arc, t in 1:N ], + genarray = [(;g, t = t) for g in data.gen, t in 1:N ], + storarray = isempty(data.storage) ? empty_data = empty_stor : [(;s, t = t) for s in data.storage, t in 1:N], Δp = corrective_action_ratio .* (data.pmax .- data.pmin) ) - - return data, dicts + + return data end -#Curve as input function update_load_data(busarray, curve) for t in eachindex(curve) for x in 1:size(busarray, 1) b = busarray[x, t] busarray[x, t] = ( - i = b.i, - pd = b.pd*curve[t], - gs = b.gs, - qd = b.qd*curve[t], - bs = b.bs, - bus_type = b.bus_type, - t = t + b=ExaPowerIO.BusData( + b.b.i, + b.b.bus_i, + b.b.type, + b.b.pd*curve[t], + b.b.qd*curve[t], + b.b.gs*curve[t], + b.b.bs*curve[t], + b.b.area, + b.b.vm, + b.b.va, + b.b.baseKV, + b.b.zone, + b.b.vmax, + b.b.vmin, + ), t=t ) end end end #Pd, Qd as input -function update_load_data(busarray, pd, qd, baseMVA, busdict) +function update_load_data(busarray, pd, qd, baseMVA) for (idx ,pd_t) in pairs(pd) - b = busarray[busdict[idx[1]], idx[2]] - busarray[busdict[idx[1]], idx[2]] = ( - i = b.i, - pd = pd_t/ baseMVA, - gs = b.gs, - qd = qd[idx[1], idx[2]] / baseMVA, - bs = b.bs, - bus_type = b.bus_type, - t = idx[2] - ) + b = busarray[idx[1], idx[2]] + busarray[idx[1], idx[2]] = ( + b=ExaPowerIO.BusData( + b.b.i, + b.b.bus_i, + b.b.type, + pd_t / baseMVA, + qd[idx[1], idx[2]] / baseMVA, + b.b.gs, + b.b.bs, + b.b.area, + b.b.vm, + b.b.va, + b.b.baseKV, + b.b.zone, + b.b.vmax, + b.b.vmin, + ), + t=idx[2], + ) end end @@ -63,30 +81,30 @@ function build_base_mpopf(core, data, N) #active, reactive power generated pg = variable(core, size(data.gen, 1), N; lvar = repeat(data.pmin, 1, N), uvar = repeat(data.pmax, 1, N)) - qg = variable(core, size(data.gen, 1), N; lvar = repeat(data.qmin, 1, N), uvar = repeat(data.qmax, 1, N)) + qg = variable(core, size(data.gen, 1), N; lvar = repeat(data.qmin, 1, N), uvar = repeat(data.qmax, 1, N)) #active, reactive power at each arc p = variable(core, size(data.arc, 1), N; lvar = repeat(-data.rate_a, 1, N), uvar = repeat(data.rate_a, 1, N)) q = variable(core, size(data.arc, 1), N; lvar = repeat(-data.rate_a, 1, N), uvar = repeat(data.rate_a, 1, N)) - o = objective(core, gen_cost(g, pg[g.i,g.t]) for g in data.genarray) + o = objective(core, gen_cost(g, pg[g.i, t]) for (g, t) in data.genarray) c_from_thermal_limit = constraint( core, - c_thermal_limit(b, p[b.f_idx, b.t], q[b.f_idx, b.t]) for b in data.barray; + c_thermal_limit(b, p[b.f_idx, t], q[b.f_idx, t]) for (b, t) in data.barray; lcon = fill(-Inf, size(data.barray)) ) c_to_thermal_limit = constraint( core, - c_thermal_limit(b, p[b.t_idx, b.t], q[b.t_idx, b.t]) for b in data.barray; + c_thermal_limit(b, p[b.t_idx, t], q[b.t_idx, t]) for (b, t) in data.barray; lcon = fill(-Inf, size(data.barray)) ) c_ramp_rate = constraint( core, - c_ramp(pg[g.i, g.t-1], pg[g.i, g.t]) for g in data.genarray[:, 2:N]; + c_ramp(pg[g.i, t-1], pg[g.i, t]) for (g, t) in data.genarray[:, 2:N]; lcon = repeat(-data.Δp, 1, N-1), ucon = repeat( data.Δp, 1, N-1) ) @@ -100,7 +118,7 @@ function build_base_mpopf(core, data, N) vars = ( pg = pg, qg = qg, - p = p, + p = p, q = q, ) @@ -120,31 +138,26 @@ function add_mpopf_cons(core, data, N, Nbus, vars, cons, form) uvar = repeat(data.vmax, 1, N), ) - c_ref_angle = constraint(core, c_ref_angle_polar(va[i,t]) for (i,t) in data.refarray) - - c_to_active_power_flow = constraint(core, c_to_active_power_flow_polar(b, p[b.f_idx, b.t], - vm[b.f_bus, b.t], vm[b.t_bus, b.t], va[b.f_bus, b.t], va[b.t_bus, b.t]) for b in data.barray) - - c_to_reactive_power_flow = constraint(core, c_to_reactive_power_flow_polar(b, q[b.f_idx, b.t], - vm[b.f_bus, b.t], vm[b.t_bus, b.t], va[b.f_bus, b.t], va[b.t_bus, b.t]) for b in data.barray) - - c_from_active_power_flow = constraint(core, c_from_active_power_flow_polar(b, p[b.t_idx, b.t], - vm[b.f_bus, b.t], vm[b.t_bus, b.t], va[b.f_bus, b.t], va[b.t_bus, b.t]) for b in data.barray) - - c_from_reactive_power_flow = constraint(core, c_from_reactive_power_flow_polar(b, q[b.t_idx, b.t], - vm[b.f_bus, b.t], vm[b.t_bus, b.t], va[b.f_bus, b.t], va[b.t_bus, b.t]) for b in data.barray) - + c_ref_angle = constraint(core, c_ref_angle_polar(va[i, t]) for (i, t) in data.refarray) + + c_to_active_power_flow = constraint(core, c_to_active_power_flow_polar(b, p[b.f_idx, t], vm[b.f_bus, t], vm[b.t_bus, t], va[b.f_bus, t], va[b.t_bus, t]) for (b, t) in data.barray) + + c_to_reactive_power_flow = constraint(core, c_to_reactive_power_flow_polar(b, q[b.f_idx, t], vm[b.f_bus, t], vm[b.t_bus, t], va[b.f_bus, t], va[b.t_bus, t]) for (b, t) in data.barray) + + c_from_active_power_flow = constraint(core, c_from_active_power_flow_polar(b, p[b.t_idx, t], vm[b.f_bus, t], vm[b.t_bus, t], va[b.f_bus, t], va[b.t_bus, t]) for (b, t) in data.barray) + + c_from_reactive_power_flow = constraint(core, c_from_reactive_power_flow_polar(b, q[b.t_idx, t], vm[b.f_bus, t], vm[b.t_bus, t], va[b.f_bus, t], va[b.t_bus, t]) for (b, t) in data.barray) + c_phase_angle_diff = constraint( - core, - c_phase_angle_diff_polar(b, va[b.f_bus, b.t], va[b.t_bus, b.t]) for b in data.barray; - lcon = repeat(data.angmin, 1, N), - ucon = repeat(data.angmax, 1, N), + core, + c_phase_angle_diff_polar(b, va[b.f_bus, t], va[b.t_bus, t]) for (b, t) in data.barray; + lcon = repeat(data.angmin, 1, N), + ucon = repeat(data.angmax, 1, N), ) - - c_active_power_balance = constraint(core, c_active_power_balance_demand_polar(b, vm[b.i, b.t]) for b in data.busarray) - c_reactive_power_balance = constraint(core, c_reactive_power_balance_demand_polar(b, vm[b.i, b.t]) for b in data.busarray) - - + + c_active_power_balance = constraint(core, c_active_power_balance_demand_polar(b, vm[b.i, t]) for (b, t) in data.busarray) + c_reactive_power_balance = constraint(core, c_reactive_power_balance_demand_polar(b, vm[b.i, t]) for (b, t) in data.busarray) + cons = (;cons..., c_ref_angle = c_ref_angle, c_to_active_power_flow = c_to_active_power_flow, @@ -153,7 +166,8 @@ function add_mpopf_cons(core, data, N, Nbus, vars, cons, form) c_from_reactive_power_flow = c_from_reactive_power_flow, c_phase_angle_diff = c_phase_angle_diff, c_active_power_balance = c_active_power_balance, - c_reactive_power_balance = c_reactive_power_balance) + c_reactive_power_balance = c_reactive_power_balance + ) vars = (;vars..., va = va, vm = vm) elseif form == :rect @@ -162,32 +176,28 @@ function add_mpopf_cons(core, data, N, Nbus, vars, cons, form) vim = variable(core, Nbus, N;) c_ref_angle = constraint(core, c_ref_angle_rect(vr[i, t], vim[i, t]) for (i, t) in data.refarray) - - c_to_active_power_flow = constraint(core, c_to_active_power_flow_rect(b, p[b.f_idx, b.t], - vr[b.f_bus, b.t], vr[b.t_bus, b.t], vim[b.f_bus, b.t], vim[b.t_bus, b.t]) for b in data.barray) - - c_to_reactive_power_flow = constraint(core, c_to_reactive_power_flow_rect(b, q[b.f_idx, b.t], - vr[b.f_bus, b.t], vr[b.t_bus, b.t], vim[b.f_bus, b.t], vim[b.t_bus, b.t]) for b in data.barray) - - c_from_active_power_flow = constraint(core, c_from_active_power_flow_rect(b, p[b.t_idx, b.t], - vr[b.f_bus, b.t], vr[b.t_bus, b.t], vim[b.f_bus, b.t], vim[b.t_bus, b.t]) for b in data.barray) - - c_from_reactive_power_flow = constraint(core, c_from_reactive_power_flow_rect(b, q[b.t_idx, b.t], - vr[b.f_bus, b.t], vr[b.t_bus, b.t], vim[b.f_bus, b.t], vim[b.t_bus, b.t]) for b in data.barray) - + + c_to_active_power_flow = constraint(core, c_to_active_power_flow_rect(b, p[b.f_idx, t], vr[b.f_bus, t], vr[b.t_bus, t], vim[b.f_bus, t], vim[b.t_bus, t]) for (b, t) in data.barray) + + c_to_reactive_power_flow = constraint(core, c_to_reactive_power_flow_rect(b, q[b.f_idx, t], vr[b.f_bus, t], vr[b.t_bus, t], vim[b.f_bus, t], vim[b.t_bus, t]) for (b, t) in data.barray) + + c_from_active_power_flow = constraint(core, c_from_active_power_flow_rect(b, p[b.t_idx, t], vr[b.f_bus, t], vr[b.t_bus, t], vim[b.f_bus, t], vim[b.t_bus, t]) for (b, t) in data.barray) + + c_from_reactive_power_flow = constraint(core, c_from_reactive_power_flow_rect(b, q[b.t_idx, t], vr[b.f_bus, t], vr[b.t_bus, t], vim[b.f_bus, t], vim[b.t_bus, t]) for (b, t) in data.barray) + c_phase_angle_diff = constraint( core, - c_phase_angle_diff_rect(b, vr[b.f_bus, b.t], vr[b.t_bus, b.t], vim[b.f_bus, b.t], vim[b.t_bus, b.t]) for b in data.barray; + c_phase_angle_diff_rect(b, vr[b.f_bus, t], vr[b.t_bus, t], vim[b.f_bus, t], vim[b.t_bus, t]) for (b, t) in data.barray; lcon = repeat(data.angmin, 1, N), ucon = repeat(data.angmax, 1, N), ) - - c_active_power_balance = constraint(core, c_active_power_balance_demand_rect(b, vr[b.i, b.t], vim[b.i, b.t]) for b in data.busarray) - c_reactive_power_balance = constraint(core, c_reactive_power_balance_demand_rect(b, vr[b.i, b.t], vim[b.i, b.t]) for b in data.busarray) + + c_active_power_balance = constraint(core, c_active_power_balance_demand_rect(b, vr[b.i, t], vim[b.i, t]) for (b, t) in data.busarray) + c_reactive_power_balance = constraint(core, c_reactive_power_balance_demand_rect(b, vr[b.i, t], vim[b.i, t]) for (b, t) in data.busarray) c_voltage_magnitude = constraint( - core, c_voltage_magnitude_rect(vr[b.i, b.t], vim[b.i, b.t]) - for b in data.busarray; + core, c_voltage_magnitude_rect(vr[b.i, t], vim[b.i, t]) + for (b, t) in data.busarray; lcon = repeat(data.vmin, 1, N).^2, ucon = repeat(data.vmax, 1, N).^2 ) @@ -200,16 +210,17 @@ function add_mpopf_cons(core, data, N, Nbus, vars, cons, form) c_phase_angle_diff = c_phase_angle_diff, c_active_power_balance = c_active_power_balance, c_reactive_power_balance = c_reactive_power_balance, - c_voltage_magnitude = c_voltage_magnitude) + c_voltage_magnitude = c_voltage_magnitude + ) vars = (;vars..., vr = vr, vim = vim) - + end - c_active_power_balance_arcs = constraint!(core, c_active_power_balance, a.bus + Nbus*(a.t-1) => p[a.i, a.t] for a in data.arcarray) - c_reactive_power_balance_arcs = constraint!(core, c_reactive_power_balance, a.bus + Nbus*(a.t-1) => q[a.i, a.t] for a in data.arcarray) + c_active_power_balance_arcs = constraint!(core, c_active_power_balance, a.bus + Nbus*(t-1) => p[a.i, t] for (a, t) in data.arcarray) + c_reactive_power_balance_arcs = constraint!(core, c_reactive_power_balance, a.bus + Nbus*(t-1) => q[a.i, t] for (a, t) in data.arcarray) - c_active_power_balance_gen = constraint!(core, c_active_power_balance, g.bus + Nbus*(g.t-1) => -pg[g.i, g.t] for g in data.genarray) - c_reactive_power_balance_gen = constraint!(core, c_reactive_power_balance, g.bus + Nbus*(g.t-1) => -qg[g.i, g.t] for g in data.genarray) + c_active_power_balance_gen = constraint!(core, c_active_power_balance, g.bus + Nbus*(t-1) => -pg[g.i, t] for (g, t) in data.genarray) + c_reactive_power_balance_gen = constraint!(core, c_reactive_power_balance, g.bus + Nbus*(t-1) => -qg[g.i, t] for (g, t) in data.genarray) return vars, cons end @@ -224,7 +235,7 @@ function build_mpopf(data, Nbus, N, form; backend = nothing, T = Float64, storag vars, cons = build_mpopf_stor_main(core, data, N, Nbus, vars, cons, form) vars, cons = add_piecewise_cons(core, data, N, vars, cons, storage_complementarity_constraint) end - + model = ExaModel(core; kwargs...) return model, vars, cons end @@ -269,20 +280,20 @@ function build_mpopf_stor_main(core, data, N, Nbus, vars, cons, form) c_active_power_balance = cons.c_active_power_balance c_reactive_power_balance = cons.c_reactive_power_balance - c_active_power_balance_stor = constraint!(core, c_active_power_balance, s.bus + Nbus*(s.t-1) => pst[s.c, s.t] for s in data.storarray) - c_reactive_power_balance_stor = constraint!(core, c_reactive_power_balance, s.bus + Nbus*(s.t-1) => qst[s.c, s.t] for s in data.storarray) + c_active_power_balance_stor = constraint!(core, c_active_power_balance, s.storage_bus + Nbus*(t-1) => pst[s.i, t] for (s, t) in data.storarray) + c_reactive_power_balance_stor = constraint!(core, c_reactive_power_balance, s.storage_bus + Nbus*(t-1) => qst[s.i, t] for (s, t) in data.storarray) - c_reactive_storage_power = constraint(core, c_reactive_stor_power(s, qst[s.c, s.t], qint[s.c, s.t], I2[s.c, s.t]) for s in data.storarray) + c_reactive_storage_power = constraint(core, c_reactive_stor_power(s, qst[s.i, t], qint[s.i, t], I2[s.i, t]) for (s, t) in data.storarray) - c_storage_transfer_thermal_limit = constraint(core, c_transfer_lim(s, pst[s.c, s.t], qst[s.c, s.t]) for s in data.storarray; lcon = fill(-Inf, size(data.storarray))) + c_storage_transfer_thermal_limit = constraint(core, c_transfer_lim(s, pst[s.i, t], qst[s.i, t]) for (s, t) in data.storarray; lcon = fill(-Inf, size(data.storarray))) if form == :polar vm = vars.vm - c_ohms = constraint(core, c_ohms_polar(pst[s.c, s.t], qst[s.c, s.t], vm[s.bus, s.t], I2[s.c, s.t]) for s in data.storarray) + c_ohms = constraint(core, c_ohms_polar(pst[s.i, t], qst[s.i, t], vm[s.storage_bus, t], I2[s.i, t]) for (s, t) in data.storarray) elseif form == :rect vr = vars.vr vim = vars.vim - c_ohms = constraint(core, c_ohms_rect(pst[s.c, s.t], qst[s.c, s.t], vr[s.bus, s.t], vim[s.bus, s.t], I2[s.c, s.t]) for s in data.storarray) + c_ohms = constraint(core, c_ohms_rect(pst[s.i, t], qst[s.i, t], vr[s.storage_bus, t], vim[s.storage_bus, t], I2[s.i, t]) for (s, t) in data.storarray) end cons = (;cons..., c_reactive_storage_power = c_reactive_storage_power, c_storage_transfer_thermal_limit = c_storage_transfer_thermal_limit, c_ohms=c_ohms) @@ -299,23 +310,23 @@ function add_piecewise_cons(core, data, N, vars, cons, storage_complementarity_c I2 = vars.I2 E = vars.E - c_active_storage_power = constraint(core, c_active_stor_power(s, pst[s.c, s.t], pstd[s.c, s.t], pstc[s.c, s.t], I2[s.c, s.t]) for s in data.storarray) + c_active_storage_power = constraint(core, c_active_stor_power(s, pst[s.i, t], pstd[s.i, t], pstc[s.i, t], I2[s.i, t]) for (s, t) in data.storarray) - c_storage_state = constraint(core, c_stor_state(s, E[s.c, s.t], E[s.c, s.t - 1], pstc[s.c, s.t], pstd[s.c, s.t]) for s in data.storarray[:, 2:N]) + c_storage_state = constraint(core, c_stor_state(s, E[s.i, t], E[s.i, t - 1], pstc[s.i, t], pstd[s.i, t]) for (s, t) in data.storarray[:, 2:N]) - c_storage_state_init = constraint(core, c_stor_state(s, E[s.c, s.t], s.Einit, pstc[s.c, s.t], pstd[s.c, s.t]) for s in data.storarray[:, 1]) + c_storage_state_init = constraint(core, c_stor_state(s, E[s.i, t], s.energy, pstc[s.i, t], pstd[s.i, t]) for (s, t) in data.storarray[:, 1]) - c_discharge_thermal_limit = constraint(core, c_discharge_lim(pstd[s.c, s.t], pstc[s.c, s.t]) for s in data.storarray; lcon = -repeat(data.srating, 1, N), ucon = repeat(data.srating, 1, N)) + c_discharge_thermal_limit = constraint(core, c_discharge_lim(pstd[s.i, t], pstc[s.i, t]) for (s, t) in data.storarray; lcon = -repeat(data.srating, 1, N), ucon = repeat(data.srating, 1, N)) - c_discharge_positivity = constraint(core, pstd[s.c, s.t] for s in data.storarray; ucon = fill(Inf, size(data.storarray))) + c_discharge_positivity = constraint(core, pstd[s.i, t] for (s, t) in data.storarray; ucon = fill(Inf, size(data.storarray))) #Complimentarity constraint if storage_complementarity_constraint - c_complementarity = constraint(core, c_comp(pstc[s.c, s.t], pstd[s.c, s.t]) for s in data.storarray) + c_complementarity = constraint(core, c_comp(pstc[s.i, t], pstd[s.i, t]) for (s, t) in data.storarray) cons = (;cons..., c_complementarity = c_complementarity) end - cons = (;cons..., + cons = (;cons..., c_active_storage_power = c_active_storage_power, c_storage_state = c_storage_state, c_storage_state_init = c_storage_state_init, @@ -331,15 +342,15 @@ function add_smooth_cons(core, data, N, vars, cons, discharge_func) I2 = vars.I2 E = vars.E - c_active_storage_power = constraint(core, c_active_storage_power_smooth(s, pst[s.c, s.t], pstd[s.c, s.t], I2[s.c, s.t]) for s in data.storarray) + c_active_storage_power = constraint(core, c_active_storage_power_smooth(s, pst[s.i, t], pstd[s.i, t], I2[s.i, t]) for (s, t) in data.storarray) - c_storage_state = constraint(core, c_storage_state_smooth(s, E[s.c, s.t], E[s.c, s.t - 1], discharge_func, pstd[s.c, s.t]) for s in data.storarray[:, 2:N]) + c_storage_state = constraint(core, c_storage_state_smooth(s, E[s.i, t], E[s.i, t - 1], discharge_func, pstd[s.i, t]) for (s, t) in data.storarray[:, 2:N]) - c_storage_state_init = constraint(core, c_storage_state_smooth(s, E[s.c, s.t], s.Einit, discharge_func, pstd[s.c, s.t]) for s in data.storarray[:, 1]) + c_storage_state_init = constraint(core, c_storage_state_smooth(s, E[s.i, t], s.energy, discharge_func, pstd[s.i, t]) for (s, t) in data.storarray[:, 1]) - c_discharge_thermal_limit = constraint(core, c_discharge_limit_smooth(pstd[s.c, s.t]) for s in data.storarray; lcon = -repeat(data.srating, 1, N), ucon = repeat(data.srating, 1, N)) + c_discharge_thermal_limit = constraint(core, c_discharge_limit_smooth(pstd[s.i, t]) for (s, t) in data.storarray; lcon = -repeat(data.srating, 1, N), ucon = repeat(data.srating, 1, N)) - cons = (;cons..., + cons = (;cons..., c_active_storage_power = c_active_storage_power, c_storage_state = c_storage_state, c_storage_state_init = c_storage_state_init, @@ -399,7 +410,7 @@ function mpopf_model( ) @assert length(curve) > 0 - data, dicts = parse_mp_power_data(filename, N, corrective_action_ratio) + data = parse_mp_power_data(filename, N, corrective_action_ratio) update_load_data(data.busarray, curve) data = convert_data(data,backend) Nbus = size(data.bus, 1) @@ -423,9 +434,9 @@ function mpopf_model( storage_complementarity_constraint = false, kwargs..., ) - - data, dicts = parse_mp_power_data(filename, N, corrective_action_ratio) - update_load_data(data.busarray, pd, qd, data.baseMVA[], dicts.bus) + + data = parse_mp_power_data(filename, N, corrective_action_ratio) + update_load_data(data.busarray, pd, qd, data.baseMVA[]) data = convert_data(data,backend) Nbus = size(data.bus, 1) @assert Nbus == size(pd, 1) @@ -449,7 +460,7 @@ function mpopf_model( ) @assert length(curve) > 0 - data, dicts = parse_mp_power_data(filename, N, corrective_action_ratio) + data = parse_mp_power_data(filename, N, corrective_action_ratio) update_load_data(data.busarray, curve) data = convert_data(data,backend) Nbus = size(data.bus, 1) @@ -473,10 +484,10 @@ function mpopf_model( storage_complementarity_constraint = false, kwargs..., ) - - - data, dicts = parse_mp_power_data(filename, N, corrective_action_ratio) - update_load_data(data.busarray, pd, qd, data.baseMVA[], dicts.bus) + + + data = parse_mp_power_data(filename, N, corrective_action_ratio) + update_load_data(data.busarray, pd, qd, data.baseMVA[]) data = convert_data(data,backend) Nbus = size(data.bus, 1) @assert Nbus == size(pd, 1) @@ -486,3 +497,4 @@ function mpopf_model( end return build_mpopf(data, Nbus, N, discharge_func, form, backend = backend, T = T, kwargs...) end + diff --git a/src/opf.jl b/src/opf.jl index ced6140..c2db6e8 100644 --- a/src/opf.jl +++ b/src/opf.jl @@ -206,7 +206,7 @@ function opf_model( kwargs..., ) - data, _ = parse_ac_power_data(filename) + data = parse_ac_power_data(filename) data = convert_data(data, backend) if form == :polar @@ -216,4 +216,4 @@ function opf_model( else error("Invalid coordinate symbol - valid options are :polar or :rect") end -end \ No newline at end of file +end diff --git a/src/parser.jl b/src/parser.jl index 28586fa..16c6a61 100644 --- a/src/parser.jl +++ b/src/parser.jl @@ -1,145 +1,51 @@ +using ExaPowerIO + convert_data(data::N, backend) where {names,N<:NamedTuple{names}} = NamedTuple{names}(convert_array(d, backend) for d in data) function parse_ac_power_data(filename) - d, f = splitdir(filename) - name, ext = splitext(f) + _, f = splitdir(filename) + name, _ = splitext(f) if isfile(joinpath(TMPDIR, name) * ".jld2") @info "Loading cached JLD2 file" loaded = JLD2.load(joinpath(TMPDIR, name) * ".jld2") - return loaded["data"], loaded["dicts"] - else - ff = if isfile(filename) - filename - elseif isfile(joinpath(TMPDIR, name) * ".m") - joinpath(TMPDIR, name) * ".m" - else - joinpath(PGLib.PGLib_opf, name * ".m") - end - @info "Loading MATPOWER file" - return process_ac_power_data(ff) + return loaded["data"] end -end - -function process_ac_power_data(filename) - data = PowerModels.parse_file(filename) - PowerModels.standardize_cost_terms!(data, order = 2) - PowerModels.calc_thermal_limits!(data) - - ref = PowerModels.build_ref(data)[:it][:pm][:nw][0] - - dicts = ( - arc = Dict(a => k for (k, a) in enumerate(ref[:arcs])), - bus = Dict(k => i for (i, (k, v)) in enumerate(ref[:bus])), - gen = Dict(k => i for (i, (k, v)) in enumerate(ref[:gen])), - branch = Dict(k => i for (i, (k, v)) in enumerate(ref[:branch])), + @info "Loading matpower file" + + library = isfile(filename) ? nothing : :pglib + data = ExaPowerIO.parse_matpower(filename; library) + + data = ( + baseMVA = [data.baseMVA], + bus = data.bus, + gen = data.gen, + arc = data.arc, + branch = data.branch, + storage = isempty(data.storage) ? empty_data = Vector{NamedTuple{(:i,), Tuple{Int64}}}() : data.storage, + ref_buses = [i for i in 1:length(data.bus) if data.bus[i].type == 3], + vmax = [bu.vmax for bu in data.bus], + vmin = [bu.vmin for bu in data.bus], + pmax = [g.pmax for g in data.gen], + pmin = [g.pmin for g in data.gen], + qmax = [g.qmax for g in data.gen], + qmin = [g.qmin for g in data.gen], + angmax = [br.angmax for br in data.branch], + angmin = [br.angmin for br in data.branch], + rate_a = [a.rate_a for a in data.arc], + vm0 = [b.vm for b in data.bus], + va0 = [b.va for b in data.bus], + pg0 = [g.pg for g in data.gen], + qg0 = [g.qg for g in data.gen], + pdmax = isempty(data.storage) ? Vector{NamedTuple{(:i,), Tuple{Int64}}}() : [s.charge_rating for s in data.storage], + pcmax = isempty(data.storage) ? Vector{NamedTuple{(:i,), Tuple{Int64}}}() : [s.discharge_rating for s in data.storage], + srating = isempty(data.storage) ? Vector{NamedTuple{(:i,), Tuple{Int64}}}() : [s.thermal_rating for s in data.storage], + emax = isempty(data.storage) ? Vector{NamedTuple{(:i,), Tuple{Int64}}}() : [s.energy_rating for s in data.storage], ) - data = ( - baseMVA = [ref[:baseMVA]], - bus = [ - begin - bus_loads = [ref[:load][l] for l in ref[:bus_loads][k]] - bus_shunts = [ref[:shunt][s] for s in ref[:bus_shunts][k]] - pd = sum(load["pd"] for load in bus_loads; init = 0.0) - gs = sum(shunt["gs"] for shunt in bus_shunts; init = 0.0) - qd = sum(load["qd"] for load in bus_loads; init = 0.0) - bs = sum(shunt["bs"] for shunt in bus_shunts; init = 0.0) - (i = dicts.bus[k], pd = pd, gs = gs, qd = qd, bs = bs, bus_type = v["bus_type"]) - end for (k, v) in ref[:bus] - ], - gen = [ - ( - i = dicts.gen[k], - cost1 = v["cost"][1], - cost2 = v["cost"][2], - cost3 = v["cost"][3], - bus = dicts.bus[v["gen_bus"]], - ) for (k, v) in ref[:gen] - ], - arc = [ - (i = k, rate_a = ref[:branch][l]["rate_a"], bus = dicts.bus[i]) for - (k, (l, i, j)) in enumerate(ref[:arcs]) - ], - branch = [ - begin - f_idx = dicts.arc[i, branch["f_bus"], branch["t_bus"]] - t_idx = dicts.arc[i, branch["t_bus"], branch["f_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"] - c1 = (-g * tr - b * ti) / ttm - c2 = (-b * tr + g * ti) / ttm - c3 = (-g * tr + b * ti) / ttm - c4 = (-b * tr - g * ti) / ttm - c5 = (g + g_fr) / ttm - c6 = (b + b_fr) / ttm - c7 = (g + g_to) - c8 = (b + b_to) - ( - i = dicts.branch[i], - j = 1, - f_idx = f_idx, - t_idx = t_idx, - f_bus = dicts.bus[branch["f_bus"]], - t_bus = dicts.bus[branch["t_bus"]], - c1 = c1, - c2 = c2, - c3 = c3, - c4 = c4, - c5 = c5, - c6 = c6, - c7 = c7, - c8 = c8, - rate_a_sq = branch["rate_a"]^2, - ) - end for (i, branch) in ref[:branch] - ], - storage = isempty(ref[:storage]) ? empty_data = Vector{NamedTuple{(:i,), Tuple{Int64}}}() : [ - begin - (c = i, - Einit = stor["energy"], - etac = stor["charge_efficiency"], - etad = stor["discharge_efficiency"], - Srating = stor["thermal_rating"], - Zr = stor["r"], - Zim = stor["x"], - Pexts = stor["ps"], - Qexts = stor["qs"], - bus = dicts.bus[stor["storage_bus"]]) - end for (i, stor) in ref[:storage] - ], - ref_buses = [dicts.bus[i] for (i, k) in ref[:ref_buses]], - vmax = [v["vmax"] for (k, v) in ref[:bus]], - vmin = [v["vmin"] for (k, v) in ref[:bus]], - pmax = [v["pmax"] for (k, v) in ref[:gen]], - pmin = [v["pmin"] for (k, v) in ref[:gen]], - qmax = [v["qmax"] for (k, v) in ref[:gen]], - qmin = [v["qmin"] for (k, v) in ref[:gen]], - rate_a = [ref[:branch][l]["rate_a"] for (k, (l, i, j)) in enumerate(ref[:arcs])], - angmax = [b["angmax"] for (i, b) in ref[:branch]], - angmin = [b["angmin"] for (i, b) in ref[:branch]], - vm0 = [v["vm"] for (k, v) in ref[:bus]], - va0 = [v["va"] for (k, v) in ref[:bus]], - pg0 = [v["pg"] for (k, v) in ref[:gen]], - qg0 = [v["qg"] for (k, v) in ref[:gen]], - pdmax = isempty(ref[:storage]) ? Vector{NamedTuple{(:i,), Tuple{Int64}}}() : [s["charge_rating"] for (i, s) in ref[:storage]], - pcmax = isempty(ref[:storage]) ? Vector{NamedTuple{(:i,), Tuple{Int64}}}() : [s["discharge_rating"] for (i, s) in ref[:storage]], - srating = isempty(ref[:storage]) ? Vector{NamedTuple{(:i,), Tuple{Int64}}}() : [s["thermal_rating"] for (i, s) in ref[:storage]], - emax = isempty(ref[:storage]) ? Vector{NamedTuple{(:i,), Tuple{Int64}}}() : [s["energy_rating"] for (i, s) in ref[:storage]], - ) - - @info "Saving JLD2 cache file" - d, f = splitdir(filename) - name,ext = splitext(f) - JLD2.save(joinpath(TMPDIR, name * ".jld2"), "data", data, "dicts", dicts) + JLD2.save(joinpath(TMPDIR, name * ".jld2"), "data", data) - return data, dicts + return data end - diff --git a/src/sc_parser.jl b/src/sc_parser.jl index 838e2d4..c60b895 100644 --- a/src/sc_parser.jl +++ b/src/sc_parser.jl @@ -1,5 +1,24 @@ +is_pr(uid::Int, L_J_pr::Int, L_J_cs::Int, producers_first::Bool)::Bool = producers_first ? uid < L_J_pr : uid >= L_J_cs +function is_pr(uid_str::String, L_J_pr::Int, L_J_cs::Int, producers_first::Bool)::Bool + uid = parse(Int, match(r"\d+", uid_str).match) + is_pr(uid, L_J_pr, L_J_cs, producers_first) +end + +function get_j_prcs(uid_str::String, L_J_pr::Int, L_J_cs::Int, producers_first::Bool) + parse(Int, match(r"\d+", uid_str).match) + 1 +end + +function get_j_pr(uid_str::String, L_J_pr::Int, L_J_cs::Int, producers_first::Bool) + offset = Int64(producers_first ? 0 : (-L_J_cs)) + parse(Int, match(r"\d+", uid_str).match) + 1 + offset +end -function parse_sc_data_static(data) +function get_j_cs(uid_str::String, L_J_pr::Int, L_J_cs::Int, producers_first::Bool) + offset = Int64(producers_first ? (-L_J_pr) : 0) + parse(Int, match(r"\d+", uid_str).match) + 1 + offset +end + +function parse_sc_data_static(data, producers_first) L_J_xf = length(data.twt_lookup) L_J_ln = length(data.ac_line_lookup) L_J_ac = L_J_ln + L_J_xf @@ -25,12 +44,12 @@ function parse_sc_data_static(data) begin ts_val = data.sdd_ts_lookup[key] j = parse(Int, match(r"\d+", val["uid"]).match) + L_J_br + 1 - j_prcs = parse(Int, match(r"\d+", val["uid"]).match) + 1 - j_pr = parse(Int, match(r"\d+", val["uid"]).match) + 1 + j_prcs = get_j_prcs(val["uid"], L_J_pr, L_J_cs, producers_first) + j_pr = get_j_pr(val["uid"], L_J_pr, L_J_cs, producers_first) bus = parse(Int, match(r"\d+", val["bus"]).match) + 1 uid = val["uid"] cost = ts_val["cost"] - (j = j, j_prcs = j_prcs, j_pr = j_pr, bus=bus, uid = uid, cost=cost) + (j = j, bus=bus, uid = uid, cost=cost, j_pr=j_pr, j_prcs=j_prcs) end for (key, val) in data.sdd_lookup if val["device_type"] == "producer" ], by = x -> x.j) @@ -41,12 +60,12 @@ function parse_sc_data_static(data) begin ts_val = data.sdd_ts_lookup[key] j = parse(Int, match(r"\d+", val["uid"]).match) + L_J_br + 1 - j_prcs = parse(Int, match(r"\d+", val["uid"]).match) + 1 - j_cs = parse(Int, match(r"\d+", val["uid"]).match) - L_J_pr + 1 + j_prcs = get_j_prcs(val["uid"], L_J_pr, L_J_cs, producers_first) + j_cs = get_j_cs(val["uid"], L_J_pr, L_J_cs, producers_first) bus = parse(Int, match(r"\d+", val["bus"]).match) + 1 uid = val["uid"] cost = ts_val["cost"] - (j = j, j_prcs = j_prcs, j_cs = j_cs, bus=bus, uid = uid, cost=cost) + (j = j, bus=bus, uid = uid, cost=cost, j_cs=j_cs, j_prcs=j_prcs) end for (key, val) in data.sdd_lookup if val["device_type"] == "consumer" ], by = x -> x.j) @@ -84,24 +103,24 @@ function parse_sc_data_static(data) uid = val["uid"] to_bus = parse(Int, match(r"\d+", val["to_bus"]).match)+1 fr_bus = parse(Int, match(r"\d+", val["fr_bus"]).match)+1 - c_su = val["connection_cost"] - c_sd = val["disconnection_cost"] - s_max = val["mva_ub_nom"] - r = val["r"] - x = val["x"] - g_sr = r / (x^2 + r^2) - b_sr = -x / (x^2 + r^2) - b_ch = val["b"] + c_su = Float64(val["connection_cost"]) + c_sd = Float64(val["disconnection_cost"]) + s_max = Float64(val["mva_ub_nom"]) + r = Float64(val["r"]) + x = Float64(val["x"]) + g_sr = Float64(r / (x^2 + r^2)) + b_sr = Float64(-x / (x^2 + r^2)) + b_ch = Float64(val["b"]) if val["additional_shunt"] == 1 - g_fr = val["g_fr"] - g_to = val["g_to"] - b_fr = val["b_fr"] - b_to = val["b_to"] + g_fr = Float64(val["g_fr"]) + g_to = Float64(val["g_to"]) + b_fr = Float64(val["b_fr"]) + b_to = Float64(val["b_to"]) else - g_fr = 0 - g_to = 0 - b_fr = 0 - b_to = 0 + g_fr = 0.0 + g_to = 0.0 + b_fr = 0.0 + b_to = 0.0 end (j = j, j_ac = j_ac, j_ln = j_ln, uid = uid, to_bus = to_bus, fr_bus = fr_bus, c_su = c_su, c_sd = c_sd, s_max = s_max, g_sr = g_sr, b_sr = b_sr, b_ch = b_ch, g_fr = g_fr, g_to = g_to, b_fr = b_fr, b_to = b_to) @@ -208,18 +227,18 @@ function parse_sc_data_static(data) [ begin ts_val = data.sdd_ts_lookup[key] - j = parse(Int, match(r"\d+", val["uid"]).match) + L_J_br + 1 - j_prcs = parse(Int, match(r"\d+", val["uid"]).match) + 1 - j_pr = parse(Int, match(r"\d+", val["uid"]).match) + 1 - bus = parse(Int, match(r"\d+", val["bus"]).match) + 1 - uid = val["uid"] - c_on = val["on_cost"] - c_su = val["startup_cost"] - c_sd = val["shutdown_cost"] - p_ru = val["p_ramp_up_ub"] - p_rd = val["p_ramp_down_ub"] - p_ru_su = val["p_startup_ramp_ub"] - p_rd_sd = val["p_shutdown_ramp_ub"] + j = Int(parse(Int, match(r"\d+", val["uid"]).match) + L_J_br + 1) + j_pr = get_j_pr(val["uid"], L_J_pr, L_J_cs, producers_first) + j_prcs = get_j_prcs(val["uid"], L_J_pr, L_J_cs, producers_first) + bus = Int(parse(Int, match(r"\d+", val["bus"]).match) + 1) + uid = String(val["uid"]) + c_on = Float64(val["on_cost"]) + c_su = Float64(val["startup_cost"]) + c_sd = Float64(val["shutdown_cost"]) + p_ru = Float64(val["p_ramp_up_ub"]) + p_rd = Float64(val["p_ramp_down_ub"]) + p_ru_su = Float64(val["p_startup_ramp_ub"]) + p_rd_sd = Float64(val["p_shutdown_ramp_ub"]) c_rgu = convert(Vector{Float64}, ts_val["p_reg_res_up_cost"]) #these need checks to see if empty c_rgd = convert(Vector{Float64}, ts_val["p_reg_res_down_cost"]) c_scr = convert(Vector{Float64}, ts_val["p_syn_res_cost"]) @@ -230,47 +249,46 @@ function parse_sc_data_static(data) c_rrd_off = convert(Vector{Float64}, ts_val["p_ramp_res_down_offline_cost"]) c_qru = convert(Vector{Float64}, ts_val["q_res_up_cost"]) c_qrd = convert(Vector{Float64}, ts_val["q_res_down_cost"]) - p_rgu_max = val["p_reg_res_up_ub"] - p_rgd_max = val["p_reg_res_down_ub"] - p_scr_max = val["p_syn_res_ub"] - p_nsc_max = val["p_nsyn_res_ub"] - p_rru_on_max = val["p_ramp_res_up_online_ub"] - p_rru_off_max = val["p_ramp_res_up_offline_ub"] - p_rrd_on_max = val["p_ramp_res_down_online_ub"] - p_rrd_off_max = val["p_ramp_res_down_offline_ub"] - p_0 = val["initial_status"]["p"] - q_0 = val["initial_status"]["q"] + p_rgu_max = Float64(val["p_reg_res_up_ub"]) + p_rgd_max = Float64(val["p_reg_res_down_ub"]) + p_scr_max = Float64(val["p_syn_res_ub"]) + p_nsc_max = Float64(val["p_nsyn_res_ub"]) + p_rru_on_max = Float64(val["p_ramp_res_up_online_ub"]) + p_rru_off_max = Float64(val["p_ramp_res_up_offline_ub"]) + p_rrd_on_max = Float64(val["p_ramp_res_down_online_ub"]) + p_rrd_off_max = Float64(val["p_ramp_res_down_offline_ub"]) + p_0 = Float64(val["initial_status"]["p"]) + q_0 = Float64(val["initial_status"]["q"]) p_max = convert(Vector{Float64}, ts_val["p_ub"]) p_min = convert(Vector{Float64}, ts_val["p_lb"]) q_max = convert(Vector{Float64}, ts_val["q_ub"]) q_min = convert(Vector{Float64}, ts_val["q_lb"]) - sus = val["startup_states"] + sus = convert(Vector{Vector{Float64}}, val["startup_states"]) - (j = j, j_prcs = j_prcs, j_pr = j_pr, bus=bus, uid = uid, c_on = c_on, c_su = c_su, c_sd = c_sd, p_ru = p_ru, p_rd = p_rd, p_ru_su = p_ru_su, p_rd_sd = p_rd_sd, + (j = j, j_pr=j_pr, j_prcs=j_prcs, bus=bus, uid = uid, c_on = c_on, c_su = c_su, c_sd = c_sd, p_ru = p_ru, p_rd = p_rd, p_ru_su = p_ru_su, p_rd_sd = p_rd_sd, c_rgu = c_rgu, c_rgd = c_rgd, c_scr = c_scr, c_nsc = c_nsc, c_rru_on = c_rru_on, c_rru_off = c_rru_off, c_rrd_on = c_rrd_on, c_rrd_off = c_rrd_off, c_qru = c_qru, c_qrd = c_qrd, p_rgu_max = p_rgu_max, p_rgd_max = p_rgd_max, p_scr_max = p_scr_max, p_nsc_max = p_nsc_max, p_rru_on_max = p_rru_on_max, p_rru_off_max=p_rru_off_max, p_rrd_on_max=p_rrd_on_max, p_rrd_off_max=p_rrd_off_max, p_0=p_0, q_0=q_0, p_max=p_max, p_min=p_min, q_max=q_max, q_min=q_min, sus=sus) end for (key, val) in data.sdd_lookup if val["device_type"] == "producer" ], by = x -> x.j), - #Consumers cons = sort( [ begin ts_val = data.sdd_ts_lookup[key] - j = parse(Int, match(r"\d+", val["uid"]).match) + L_J_br + 1 - j_prcs = parse(Int, match(r"\d+", val["uid"]).match) + 1 - j_cs = parse(Int, match(r"\d+", val["uid"]).match) - L_J_pr + 1 - bus = parse(Int, match(r"\d+", val["bus"]).match) + 1 - uid = val["uid"] - c_on = val["on_cost"] - c_su = val["startup_cost"] - c_sd = val["shutdown_cost"] - p_ru = val["p_ramp_up_ub"] - p_rd = val["p_ramp_down_ub"] - p_ru_su = val["p_startup_ramp_ub"] - p_rd_sd = val["p_shutdown_ramp_ub"] + j = Int(parse(Int, match(r"\d+", val["uid"]).match) + L_J_br + 1) + j_prcs = get_j_prcs(val["uid"], L_J_pr, L_J_cs, producers_first) + j_cs = get_j_cs(val["uid"], L_J_pr, L_J_cs, producers_first) + bus = Int(parse(Int, match(r"\d+", val["bus"]).match) + 1) + uid = String(val["uid"]) + c_on = Float64(val["on_cost"]) + c_su = Float64(val["startup_cost"]) + c_sd = Float64(val["shutdown_cost"]) + p_ru = Float64(val["p_ramp_up_ub"]) + p_rd = Float64(val["p_ramp_down_ub"]) + p_ru_su = Float64(val["p_startup_ramp_ub"]) + p_rd_sd = Float64(val["p_shutdown_ramp_ub"]) c_rgu = convert(Vector{Float64}, ts_val["p_reg_res_up_cost"]) c_rgd = convert(Vector{Float64}, ts_val["p_reg_res_down_cost"]) c_scr = convert(Vector{Float64}, ts_val["p_syn_res_cost"]) @@ -281,26 +299,26 @@ function parse_sc_data_static(data) c_rrd_off = convert(Vector{Float64}, ts_val["p_ramp_res_down_offline_cost"]) c_qru = convert(Vector{Float64}, ts_val["q_res_up_cost"]) c_qrd = convert(Vector{Float64}, ts_val["q_res_down_cost"]) - p_rgu_max = val["p_reg_res_up_ub"] - p_rgd_max = val["p_reg_res_down_ub"] - p_scr_max = val["p_syn_res_ub"] - p_nsc_max = val["p_nsyn_res_ub"] - p_rru_on_max = val["p_ramp_res_up_online_ub"] - p_rru_off_max = val["p_ramp_res_up_offline_ub"] - p_rrd_on_max = val["p_ramp_res_down_online_ub"] - p_rrd_off_max = val["p_ramp_res_down_offline_ub"] - p_0 = val["initial_status"]["p"] - q_0 = val["initial_status"]["q"] + p_rgu_max = Float64(val["p_reg_res_up_ub"]) + p_rgd_max = Float64(val["p_reg_res_down_ub"]) + p_scr_max = Float64(val["p_syn_res_ub"]) + p_nsc_max = Float64(val["p_nsyn_res_ub"]) + p_rru_on_max = Float64(val["p_ramp_res_up_online_ub"]) + p_rru_off_max = Float64(val["p_ramp_res_up_offline_ub"]) + p_rrd_on_max = Float64(val["p_ramp_res_down_online_ub"]) + p_rrd_off_max = Float64(val["p_ramp_res_down_offline_ub"]) + p_0 = Float64(val["initial_status"]["p"]) + q_0 = Float64(val["initial_status"]["q"]) p_max = convert(Vector{Float64}, ts_val["p_ub"]) p_min = convert(Vector{Float64}, ts_val["p_lb"]) q_max = convert(Vector{Float64}, ts_val["q_ub"]) q_min = convert(Vector{Float64}, ts_val["q_lb"]) - sus = val["startup_states"] + sus = convert(Vector{Vector{Float64}}, val["startup_states"]) - (j = j, j_prcs = j_prcs, j_cs = j_cs, bus=bus, uid = uid, c_on = c_on, c_su = c_su, c_sd = c_sd, p_ru = p_ru, p_rd = p_rd, p_ru_su = p_ru_su, p_rd_sd = p_rd_sd, + (j = j, j_cs=j_cs, j_prcs=j_prcs, bus=bus, uid = uid, c_on = c_on, c_su = c_su, c_sd = c_sd, p_ru = p_ru, p_rd = p_rd, p_ru_su = p_ru_su, p_rd_sd = p_rd_sd, c_rgu = c_rgu, c_rgd = c_rgd, c_scr = c_scr, c_nsc = c_nsc, c_rru_on = c_rru_on, c_rru_off = c_rru_off, c_rrd_on = c_rrd_on, c_rrd_off = c_rrd_off, c_qru = c_qru, c_qrd = c_qrd, p_rgu_max = p_rgu_max, p_rgd_max = p_rgd_max, p_scr_max = p_scr_max, p_nsc_max = p_nsc_max, p_rru_on_max = p_rru_on_max, - p_rru_off_max=p_rru_off_max, p_rrd_on_max=p_rrd_on_max, p_rrd_off_max=p_rrd_off_max, p_0=p_0, q_0=q_0, p_max=p_max, p_min=p_min, q_max=q_max, q_min=q_min, sus=sus) + p_rru_off_max=p_rru_off_max, p_rrd_on_max=p_rrd_on_max, p_rrd_off_max=p_rrd_off_max, p_0=p_0, q_0=q_0, p_max=p_max, p_min=p_min, q_max=q_max, q_min=q_min, sus=sus) end for (key, val) in data.sdd_lookup if val["device_type"] == "consumer" ] , by = x -> x.j), @@ -341,43 +359,60 @@ function parse_sc_data_static(data) ], by = x -> x.n), active_reserve_set_pr = [ - (i = parse(Int, match(r"\d+", bus["uid"]).match) + 1, j = parse(Int, match(r"\d+", device["uid"]).match) + L_J_br + 1, n = parse(Int, match(r"\d+", uid).match) + 1, n_p = parse(Int, match(r"\d+", uid).match) + 1, - j_prcs = parse(Int, match(r"\d+", device["uid"]).match) + 1, j_pr = parse(Int, match(r"\d+", device["uid"]).match) + 1) + (i = parse(Int, match(r"\d+", bus["uid"]).match) + 1, + j = parse(Int, match(r"\d+", device["uid"]).match) + L_J_br + 1, + n = parse(Int, match(r"\d+", uid).match) + 1, + n_p = parse(Int, match(r"\d+", uid).match) + 1, + j_pr=get_j_pr(device["uid"], L_J_pr, L_J_cs, producers_first), + j_prcs=get_j_prcs(device["uid"], L_J_pr, L_J_cs, producers_first), + ) for uid in data.azr_ids for bus in values(data.bus_lookup) if uid in bus["active_reserve_uids"] for device in values(data.sdd_lookup) - if device["bus"] == bus["uid"] && parse(Int, match(r"\d+", device["uid"]).match) < L_J_pr + if device["bus"] == bus["uid"] && is_pr(device["uid"], L_J_pr, L_J_cs, producers_first) ], active_reserve_set_cs = [ - (i = parse(Int, match(r"\d+", bus["uid"]).match) + 1, j = parse(Int, match(r"\d+", device["uid"]).match) + L_J_br + 1, n = parse(Int, match(r"\d+", uid).match) + 1, n_p = parse(Int, match(r"\d+", uid).match) + 1, - j_prcs = parse(Int, match(r"\d+", device["uid"]).match) + 1, j_cs = parse(Int, match(r"\d+", device["uid"]).match) + 1 - L_J_pr) + (i = parse(Int, match(r"\d+", bus["uid"]).match) + 1, + j = parse(Int, match(r"\d+", device["uid"]).match) + L_J_br + 1, + n = parse(Int, match(r"\d+", uid).match) + 1, + n_p = parse(Int, match(r"\d+", uid).match) + 1, + j_cs=get_j_cs(device["uid"], L_J_pr, L_J_cs, producers_first), + j_prcs=get_j_prcs(device["uid"], L_J_pr, L_J_cs, producers_first),) for uid in data.azr_ids for bus in values(data.bus_lookup) if uid in bus["active_reserve_uids"] for device in values(data.sdd_lookup) - if device["bus"] == bus["uid"] && parse(Int, match(r"\d+", device["uid"]).match) >= L_J_pr + if device["bus"] == bus["uid"] && !is_pr(device["uid"], L_J_pr, L_J_cs, producers_first) ], reactive_reserve_set_pr = [ - (i = parse(Int, match(r"\d+", bus["uid"]).match) + 1, j = parse(Int, match(r"\d+", device["uid"]).match) + L_J_br + 1, n = parse(Int, match(r"\d+", uid).match) + L_N_p + 1, n_q = parse(Int, match(r"\d+", uid).match) + 1, - j_prcs = parse(Int, match(r"\d+", device["uid"]).match) + 1, j_pr = parse(Int, match(r"\d+", device["uid"]).match) + 1) + (i = parse(Int, match(r"\d+", bus["uid"]).match) + 1, + j = parse(Int, match(r"\d+", device["uid"]).match) + L_J_br + 1, + n = parse(Int, match(r"\d+", uid).match) + L_N_p + 1, + n_q = parse(Int, match(r"\d+", uid).match) + 1, + j_pr=get_j_pr(device["uid"], L_J_pr, L_J_cs, producers_first), + j_prcs=get_j_prcs(device["uid"], L_J_pr, L_J_cs, producers_first),) for uid in data.rzr_ids for bus in values(data.bus_lookup) if uid in bus["reactive_reserve_uids"] for device in values(data.sdd_lookup) - if device["bus"] == bus["uid"] && parse(Int, match(r"\d+", device["uid"]).match) < L_J_pr + if device["bus"] == bus["uid"] && is_pr(device["uid"], L_J_pr, L_J_cs, producers_first) ], reactive_reserve_set_cs = [ - (i = parse(Int, match(r"\d+", bus["uid"]).match) + 1, j = parse(Int, match(r"\d+", device["uid"]).match) + L_J_br + 1, n = parse(Int, match(r"\d+", uid).match) + L_N_p + 1, n_q = parse(Int, match(r"\d+", uid).match) + 1, - j_prcs = parse(Int, match(r"\d+", device["uid"]).match) + 1, j_cs = parse(Int, match(r"\d+", device["uid"]).match) + 1 - L_J_pr) + (i = parse(Int, match(r"\d+", bus["uid"]).match) + 1, + j = parse(Int, match(r"\d+", device["uid"]).match) + L_J_br + 1, + n = parse(Int, match(r"\d+", uid).match) + L_N_p + 1, + n_q = parse(Int, match(r"\d+", uid).match) + 1, + j_cs=get_j_cs(device["uid"], L_J_pr, L_J_cs, producers_first), + j_prcs=get_j_prcs(device["uid"], L_J_pr, L_J_cs, producers_first),) for uid in data.rzr_ids for bus in values(data.bus_lookup) if uid in bus["reactive_reserve_uids"] for device in values(data.sdd_lookup) - if device["bus"] == bus["uid"] && parse(Int, match(r"\d+", device["uid"]).match) >= L_J_pr + if device["bus"] == bus["uid"] && !is_pr(device["uid"], L_J_pr, L_J_cs, producers_first) ], ) @@ -426,7 +461,9 @@ function get_as(dt, t) end function parse_sc_data(data, uc_data, data_json) - sc_data, lengths, cost_vector_pr, cost_vector_cs = parse_sc_data_static(data) + producers_first = data.sdd_lookup[minimum(keys(data.sdd_lookup))]["device_type"] == "producer" + sc_data, lengths, cost_vector_pr, cost_vector_cs = parse_sc_data_static(data, producers_first) + @info typeof(sc_data) (L_J_xf, L_J_ln, L_J_ac, L_J_dc, L_J_br, L_J_cs, L_J_pr, L_J_cspr, L_J_sh, I, L_T, L_N_p, L_N_q) = lengths @@ -438,18 +475,17 @@ function parse_sc_data(data, uc_data, data_json) add_status_flags(uc_data["time_series_output"]["simple_dispatchable_device"], data.sdd_lookup) ε_time = 1e-6 - T_supc_pr = [ (j = parse(Int, match(r"\d+", val["uid"]).match) + L_J_br + 1, - j_prcs = parse(Int, match(r"\d+", val["uid"]).match) + 1, - j_pr = parse(Int, match(r"\d+", val["uid"]).match) + 1, + j_pr=get_j_pr(val["uid"], L_J_pr, L_J_cs, producers_first), + j_prcs=get_j_prcs(val["uid"], L_J_pr, L_J_cs, producers_first), t = t, t_prime = t_prime, p_supc = data.sdd_ts_lookup[key]["p_lb"][t_prime] - val["p_startup_ramp_ub"]*(get_as(dt, t_prime)[3] - get_as(dt, t)[3]), u_su = uc["su_status"][t_prime] ) for (key, val) in data.sdd_lookup - if parse(Int, match(r"\d+", val["uid"]).match) < L_J_pr + if is_pr(val["uid"], L_J_pr, L_J_cs, producers_first) for t in periods for t_prime in periods if t_prime > t && data.sdd_ts_lookup[key]["p_lb"][t_prime] - val["p_startup_ramp_ub"]*(get_as(dt, t_prime)[3] - get_as(dt, t)[3]) > 0 @@ -469,15 +505,15 @@ function parse_sc_data(data, uc_data, data_json) T_supc_cs = [ (j = parse(Int, match(r"\d+", val["uid"]).match) + L_J_br + 1, - j_prcs = parse(Int, match(r"\d+", val["uid"]).match) + 1, - j_cs = parse(Int, match(r"\d+", val["uid"]).match) + 1 - L_J_pr, + j_cs=get_j_cs(val["uid"], L_J_pr, L_J_cs, producers_first), + j_prcs=get_j_prcs(val["uid"], L_J_pr, L_J_cs, producers_first), t = t, t_prime = t_prime, p_supc = data.sdd_ts_lookup[key]["p_lb"][t_prime] - val["p_startup_ramp_ub"]*(get_as(dt, t_prime)[3] - get_as(dt, t)[3]), u_su = uc["su_status"][t_prime] ) for (key, val) in data.sdd_lookup - if parse(Int, match(r"\d+", val["uid"]).match) >= L_J_pr + if !is_pr(val["uid"], L_J_pr, L_J_cs, producers_first) for t in periods for t_prime in periods if t_prime > t && data.sdd_ts_lookup[key]["p_lb"][t_prime] - val["p_startup_ramp_ub"]*(get_as(dt, t_prime)[3] - get_as(dt, t)[3]) > 0 @@ -502,9 +538,9 @@ function parse_sc_data(data, uc_data, data_json) for t_prime in periods for (key, val) in data.sdd_lookup if t_prime == 1 && t >= t_prime - p_sdpc[parse(Int, match(r"\d+", val["uid"]).match)+1, t, t_prime] = val["initial_status"]["p"] - val["p_shutdown_ramp_ub"]*(get_as(dt, t)[3] - get_as(dt, t_prime)[1]) + p_sdpc[get_j_prcs(val["uid"], L_J_pr, L_J_cs, producers_first), t, t_prime] = val["initial_status"]["p"] - val["p_shutdown_ramp_ub"]*(get_as(dt, t)[3] - get_as(dt, t_prime)[1]) elseif t >= t_prime - p_sdpc[parse(Int, match(r"\d+", val["uid"]).match)+1, t, t_prime] = data.sdd_ts_lookup[key]["p_lb"][t_prime-1]- val["p_shutdown_ramp_ub"]*(get_as(dt, t)[3] - get_as(dt, t_prime)[1]) + p_sdpc[get_j_prcs(val["uid"], L_J_pr, L_J_cs, producers_first), t, t_prime] = data.sdd_ts_lookup[key]["p_lb"][t_prime-1]- val["p_shutdown_ramp_ub"]*(get_as(dt, t)[3] - get_as(dt, t_prime)[1]) end end end @@ -512,15 +548,15 @@ function parse_sc_data(data, uc_data, data_json) T_sdpc_pr = [ (j = parse(Int, match(r"\d+", val["uid"]).match) + L_J_br + 1, - j_prcs = parse(Int, match(r"\d+", val["uid"]).match) + 1, - j_pr = parse(Int, match(r"\d+", val["uid"]).match) + 1, + j_pr=get_j_pr(val["uid"], L_J_pr, L_J_cs, producers_first), + j_prcs=get_j_prcs(val["uid"], L_J_pr, L_J_cs, producers_first), t = t, t_prime = t_prime, p_sdpc = p_sdpc[parse(Int, match(r"\d+", val["uid"]).match)+1, t, t_prime], u_sd = uc["sd_status"][t_prime] ) for (key, val) in data.sdd_lookup - if parse(Int, match(r"\d+", val["uid"]).match) < L_J_pr + if is_pr(val["uid"], L_J_pr, L_J_cs, producers_first) for t in periods for t_prime in periods if t_prime <= t && p_sdpc[parse(Int, match(r"\d+", val["uid"]).match)+1, t, t_prime] > 0 @@ -539,15 +575,15 @@ function parse_sc_data(data, uc_data, data_json) T_sdpc_cs = [ (j = parse(Int, match(r"\d+", val["uid"]).match) + L_J_br + 1, - j_prcs = parse(Int, match(r"\d+", val["uid"]).match) + 1, - j_cs = parse(Int, match(r"\d+", val["uid"]).match) + 1 - L_J_pr, + j_cs=get_j_cs(val["uid"], L_J_pr, L_J_cs, producers_first), + j_prcs=get_j_prcs(val["uid"], L_J_pr, L_J_cs, producers_first), t = t, t_prime = t_prime, p_sdpc = p_sdpc[parse(Int, match(r"\d+", val["uid"]).match)+1, t, t_prime], u_sd = uc["sd_status"][t_prime] ) for (key, val) in data.sdd_lookup - if parse(Int, match(r"\d+", val["uid"]).match) >= L_J_pr + if !is_pr(val["uid"], L_J_pr, L_J_cs, producers_first) for t in periods for t_prime in periods if t_prime <= t && p_sdpc[parse(Int, match(r"\d+", val["uid"]).match)+1, t, t_prime] > 0 @@ -567,10 +603,17 @@ function parse_sc_data(data, uc_data, data_json) W_en_max_pr = Vector{@NamedTuple{w_en_max_pr_ind::Int, j::Int, j_prcs::Int, j_pr::Int, a_en_max_start::Float64, a_en_max_end::Float64, e_max::Float64}}() w_en_max_pr_ind = 1 for val in values(data.sdd_lookup) - if parse(Int, match(r"\d+", val["uid"]).match) < L_J_pr + if is_pr(val["uid"], L_J_pr, L_J_cs, producers_first) for w in val["energy_req_ub"] - push!(W_en_max_pr, (w_en_max_pr_ind=w_en_max_pr_ind, j=parse(Int, match(r"\d+", val["uid"]).match) + L_J_br + 1, j_prcs = parse(Int, match(r"\d+", val["uid"]).match) + 1, - j_pr = parse(Int, match(r"\d+", val["uid"]).match) + 1, a_en_max_start = w[1], a_en_max_end = w[2], e_max = w[3])) + push!(W_en_max_pr, ( + w_en_max_pr_ind=w_en_max_pr_ind, + j=parse(Int, match(r"\d+", val["uid"]).match) + L_J_br + 1, + j_prcs=get_j_prcs(val["uid"], L_J_pr, L_J_cs, producers_first), + j_pr=get_j_pr(val["uid"], L_J_pr, L_J_cs, producers_first), + a_en_max_start = w[1], + a_en_max_end = w[2], + e_max = w[3]) + ) w_en_max_pr_ind +=1 end end @@ -581,10 +624,17 @@ function parse_sc_data(data, uc_data, data_json) W_en_max_cs = Vector{@NamedTuple{w_en_max_cs_ind::Int, j::Int, j_prcs::Int, j_cs::Int, a_en_max_start::Float64, a_en_max_end::Float64, e_max::Float64}}() w_en_max_cs_ind = 1 for val in values(data.sdd_lookup) - if parse(Int, match(r"\d+", val["uid"]).match) >= L_J_pr + if !is_pr(val["uid"], L_J_pr, L_J_cs, producers_first) for w in val["energy_req_ub"] - push!(W_en_max_cs, (w_en_max_cs_ind=w_en_max_cs_ind, j=parse(Int, match(r"\d+", val["uid"]).match) + L_J_br + 1, j_prcs = parse(Int, match(r"\d+", val["uid"]).match) + 1, - j_cs = parse(Int, match(r"\d+", val["uid"]).match) + 1 - L_J_pr, a_en_max_start = w[1], a_en_max_end = w[2], e_max = w[3])) + push!(W_en_max_cs, ( + w_en_max_cs_ind=w_en_max_cs_ind, + j=parse(Int, match(r"\d+", val["uid"]).match) + L_J_br + 1, + j_prcs=get_j_prcs(val["uid"], L_J_pr, L_J_cs, producers_first), + j_cs=get_j_cs(val["uid"], L_J_pr, L_J_cs, producers_first), + a_en_max_start = w[1], + a_en_max_end = w[2], + e_max = w[3]) + ) w_en_max_cs_ind +=1 end end @@ -595,12 +645,12 @@ function parse_sc_data(data, uc_data, data_json) W_en_min_pr = Vector{@NamedTuple{w_en_min_pr_ind::Int, j::Int, j_prcs::Int, j_pr::Int, a_en_min_start::Float64, a_en_min_end::Float64, e_min::Float64}}() w_en_min_pr_ind = 1 for val in values(data.sdd_lookup) - if parse(Int, match(r"\d+", val["uid"]).match) < L_J_pr + if is_pr(val["uid"], L_J_pr, L_J_cs, producers_first) for w in val["energy_req_lb"] push!(W_en_min_pr, (w_en_min_pr_ind = w_en_min_pr_ind, j = parse(Int, match(r"\d+", val["uid"]).match) + L_J_br + 1, - j_prcs = parse(Int, match(r"\d+", val["uid"]).match) + 1, - j_pr = parse(Int, match(r"\d+", val["uid"]).match) + 1, + j_prcs=get_j_prcs(val["uid"], L_J_pr, L_J_cs, producers_first), + j_pr=get_j_pr(val["uid"], L_J_pr, L_J_cs, producers_first), a_en_min_start = w[1], a_en_min_end = w[2], e_min = w[3])) @@ -614,12 +664,12 @@ function parse_sc_data(data, uc_data, data_json) W_en_min_cs = Vector{@NamedTuple{w_en_min_cs_ind::Int, j::Int, j_prcs::Int, j_cs::Int, a_en_min_start::Float64, a_en_min_end::Float64, e_min::Float64}}() w_en_min_cs_ind = 1 for val in values(data.sdd_lookup) - if parse(Int, match(r"\d+", val["uid"]).match) >= L_J_pr + if !is_pr(val["uid"], L_J_pr, L_J_cs, producers_first) for w in val["energy_req_lb"] push!(W_en_min_cs, (w_en_min_cs_ind = w_en_min_cs_ind, j = parse(Int, match(r"\d+", val["uid"]).match) + L_J_br + 1, - j_prcs = parse(Int, match(r"\d+", val["uid"]).match) + 1, - j_cs = parse(Int, match(r"\d+", val["uid"]).match) + 1 - L_J_pr, + j_prcs=get_j_prcs(val["uid"], L_J_pr, L_J_cs, producers_first), + j_cs=get_j_cs(val["uid"], L_J_pr, L_J_cs, producers_first), a_en_min_start = w[1], a_en_min_end = w[2], e_min = w[3])) @@ -633,15 +683,15 @@ function parse_sc_data(data, uc_data, data_json) T_w_en_max_pr = Vector{@NamedTuple{w_en_max_pr_ind::Int, j::Int, j_prcs::Int, j_pr::Int, t::Int, dt::Float64}}() w_en_max_pr_ind = 0 for val in values(data.sdd_lookup) - if parse(Int, match(r"\d+", val["uid"]).match) < L_J_pr + if is_pr(val["uid"], L_J_pr, L_J_cs, producers_first) for w in val["energy_req_ub"] w_en_max_pr_ind += 1 for t in periods if w[1] + ε_time < get_as(dt, t)[2] && get_as(dt, t)[2] <= w[2] + ε_time push!(T_w_en_max_pr, (w_en_max_pr_ind = w_en_max_pr_ind, j = parse(Int, match(r"\d+", val["uid"]).match) + L_J_br + 1, - j_prcs = parse(Int, match(r"\d+", val["uid"]).match) + 1, - j_pr = parse(Int, match(r"\d+", val["uid"]).match) + 1, + j_prcs=get_j_prcs(val["uid"], L_J_pr, L_J_cs, producers_first), + j_pr=get_j_pr(val["uid"], L_J_pr, L_J_cs, producers_first), t = t, dt=dt[t])) end @@ -653,15 +703,15 @@ function parse_sc_data(data, uc_data, data_json) T_w_en_max_cs = Vector{@NamedTuple{w_en_max_cs_ind::Int, j::Int, j_prcs::Int, j_cs::Int, t::Int, dt::Float64}}() w_en_max_cs_ind = 0 for val in values(data.sdd_lookup) - if parse(Int, match(r"\d+", val["uid"]).match) >= L_J_pr + if !is_pr(val["uid"], L_J_pr, L_J_cs, producers_first) for w in val["energy_req_ub"] w_en_max_cs_ind += 1 for t in periods if w[1] + ε_time < get_as(dt, t)[2] && get_as(dt, t)[2] <= w[2] + ε_time push!(T_w_en_max_cs, (w_en_max_cs_ind = w_en_max_cs_ind, j = parse(Int, match(r"\d+", val["uid"]).match) + L_J_br + 1, - j_prcs = parse(Int, match(r"\d+", val["uid"]).match) + 1, - j_cs = parse(Int, match(r"\d+", val["uid"]).match) + 1 - L_J_pr, + j_prcs=get_j_prcs(val["uid"], L_J_pr, L_J_cs, producers_first), + j_cs=get_j_cs(val["uid"], L_J_pr, L_J_cs, producers_first), t = t, dt=dt[t])) end @@ -673,14 +723,14 @@ function parse_sc_data(data, uc_data, data_json) T_w_en_min_pr = Vector{@NamedTuple{w_en_min_pr_ind::Int, j::Int, j_prcs::Int, j_pr::Int, t::Int, dt::Float64}}() w_en_min_pr_ind = 0 for val in values(data.sdd_lookup) - if parse(Int, match(r"\d+", val["uid"]).match) < L_J_pr + if is_pr(val["uid"], L_J_pr, L_J_cs, producers_first) for w in val["energy_req_lb"] w_en_min_pr_ind += 1 for t in periods if w[1] + ε_time < get_as(dt, t)[2] && get_as(dt, t)[2] <= w[2] + ε_time push!(T_w_en_min_pr, (w_en_min_pr_ind=w_en_min_pr_ind, j = parse(Int, match(r"\d+", val["uid"]).match) + L_J_br + 1, - j_prcs = parse(Int, match(r"\d+", val["uid"]).match) + 1, - j_pr = parse(Int, match(r"\d+", val["uid"]).match) + 1, + j_prcs=get_j_prcs(val["uid"], L_J_pr, L_J_cs, producers_first), + j_pr=get_j_pr(val["uid"], L_J_pr, L_J_cs, producers_first), t = t, dt = dt[t])) end @@ -692,14 +742,14 @@ function parse_sc_data(data, uc_data, data_json) T_w_en_min_cs = Vector{@NamedTuple{w_en_min_cs_ind::Int, j::Int, j_prcs::Int, j_cs::Int, t::Int, dt::Float64}}() w_en_min_cs_ind = 0 for val in values(data.sdd_lookup) - if parse(Int, match(r"\d+", val["uid"]).match) >= L_J_pr + if !is_pr(val["uid"], L_J_pr, L_J_cs, producers_first) for w in val["energy_req_lb"] w_en_min_cs_ind += 1 for t in periods if w[1] + ε_time < get_as(dt, t)[2] && get_as(dt, t)[2] <= w[2] + ε_time push!(T_w_en_min_cs, (w_en_min_cs_ind=w_en_min_cs_ind, j = parse(Int, match(r"\d+", val["uid"]).match) + L_J_br + 1, - j_prcs = parse(Int, match(r"\d+", val["uid"]).match) + 1, - j_cs = parse(Int, match(r"\d+", val["uid"]).match) + 1 - L_J_pr, + j_prcs=get_j_prcs(val["uid"], L_J_pr, L_J_cs, producers_first), + j_cs=get_j_cs(val["uid"], L_J_pr, L_J_cs, producers_first), t = t, dt = dt[t])) end @@ -717,7 +767,7 @@ function parse_sc_data(data, uc_data, data_json) for (t, cost_t) in enumerate(pc.cost) for (m, cost_tm) in enumerate(cost_t) c_en, p_max = cost_tm - push!(p_jtm_flattened_pr, (flat_k=flat_k, j=j, j_prcs=j_prcs, j_pr=j_pr, t=t, m=m, c_en=c_en, p_max=p_max)) + push!(p_jtm_flattened_pr, (flat_k=flat_k, j=j, j_prcs=j_prcs, j_pr=pc.j_pr, t=t, m=m, c_en=c_en, p_max=p_max)) flat_k+=1 end end @@ -732,7 +782,7 @@ function parse_sc_data(data, uc_data, data_json) for (t, cost_t) in enumerate(pc.cost) for (m, cost_tm) in enumerate(cost_t) c_en, p_max = cost_tm - push!(p_jtm_flattened_cs, (flat_k=flat_k, j=j, j_prcs=j_prcs, j_cs=j_cs, t=t, m=m, c_en=c_en, p_max=p_max)) + push!(p_jtm_flattened_cs, (flat_k=flat_k, j=j, j_prcs=j_prcs, j_cs=pc.j_cs, t=t, m=m, c_en=c_en, p_max=p_max)) flat_k+=1 end end @@ -840,11 +890,11 @@ function parse_sc_data(data, uc_data, data_json) for uc in uc_data["time_series_output"]["simple_dispatchable_device"] if p.uid == uc["uid"]], - prarray_pqbounds = isempty(val for val in values(data.sdd_lookup) if parse(Int, match(r"\d+", val["uid"]).match) < L_J_pr && val["q_bound_cap"]==1) ? + prarray_pqbounds = isempty(val for val in values(data.sdd_lookup) if is_pr(val["uid"], L_J_pr, L_J_cs, producers_first) && val["q_bound_cap"]==1) ? empty_data = Vector{NamedTuple{(:j, :jprcs, :j_pr, :u_on, :sum2_T_supc_pr_jt, :sum2_T_sdpc_pr_jt, :beta_max, :beta_min, :q_max_p0, :q_min_p0, :t), Tuple{Int64, Int64, Int64, Int64, Int64, Int64, Float64, Float64, Float64, Float64, Int64}}}() : [ (;j = parse(Int, match(r"\d+", val["uid"]).match) + L_J_br + 1, - j_prcs = parse(Int, match(r"\d+", val["uid"]).match) + 1, - j_pr = parse(Int, match(r"\d+", val["uid"]).match) + 1, + j_prcs=get_j_prcs(val["uid"], L_J_pr, L_J_cs, producers_first), + j_pr=get_j_pr(val["uid"], L_J_pr, L_J_cs, producers_first), u_on = uc["on_status"][t], sum2_T_supc_pr_jt=sum2_T_supc_pr[p.j_pr, t], sum2_T_sdpc_pr_jt=sum2_T_sdpc_pr[p.j_pr, t], @@ -855,15 +905,15 @@ function parse_sc_data(data, uc_data, data_json) t=t ) for val in values(data.sdd_lookup), t in periods - if parse(Int, match(r"\d+", val["uid"]).match) < L_J_pr && val["q_bound_cap"]==1 + if is_pr(val["uid"], L_J_pr, L_J_cs, producers_first) && val["q_bound_cap"]==1 for uc in uc_data["time_series_output"]["simple_dispatchable_device"] if p.uid == uc["uid"]], - prarray_pqe = isempty(val for val in values(data.sdd_lookup) if parse(Int, match(r"\d+", val["uid"]).match) < L_J_pr && val["q_bound_cap"]==1) ? + prarray_pqe = isempty(val for val in values(data.sdd_lookup) if is_pr(val["uid"], L_J_pr, L_J_cs, producers_first) && val["q_bound_cap"]==1) ? empty_data = Vector{NamedTuple{(:j, :jprcs, :j_pr, :u_on, :sum2_T_supc_pr_jt, :sum2_T_sdpc_pr_jt, :beta, :q_p0, :t), Tuple{Int64, Int64, Int64, Int64, Int64, Int64, Float64, Float64, Int64}}}() : [ (;j = parse(Int, match(r"\d+", val["uid"]).match) + L_J_br + 1, - j_prcs = parse(Int, match(r"\d+", val["uid"]).match) + 1, - j_pr = parse(Int, match(r"\d+", val["uid"]).match) + 1, + j_prcs=get_j_prcs(val["uid"], L_J_pr, L_J_cs, producers_first), + j_pr=get_j_pr(val["uid"], L_J_pr, L_J_cs, producers_first), u_on = uc["on_status"][t], sum2_T_supc_pr_jt=sum2_T_supc_pr[p.j_pr, t], sum2_T_sdpc_pr_jt=sum2_T_sdpc_pr[p.j_pr, t], @@ -872,7 +922,7 @@ function parse_sc_data(data, uc_data, data_json) t=t ) for val in values(data.sdd_lookup), t in periods - if parse(Int, match(r"\d+", val["uid"]).match) < L_J_pr && val["q_linear_cap"]==1 + if is_pr(val["uid"], L_J_pr, L_J_cs, producers_first) && val["q_linear_cap"]==1 for uc in uc_data["time_series_output"]["simple_dispatchable_device"] if p.uid == uc["uid"]], @@ -889,11 +939,11 @@ function parse_sc_data(data, uc_data, data_json) for uc in uc_data["time_series_output"]["simple_dispatchable_device"] if p.uid == uc["uid"]], - csarray_pqbounds = isempty(val for val in values(data.sdd_lookup) if parse(Int, match(r"\d+", val["uid"]).match) < L_J_pr && val["q_bound_cap"]==1) ? + csarray_pqbounds = isempty(val for val in values(data.sdd_lookup) if !is_pr(val["uid"], L_J_pr, L_J_cs, producers_first) && val["q_bound_cap"]==1) ? empty_data = Vector{NamedTuple{(:j, :jprcs, :j_cs, :u_on, :sum2_T_supc_cs_jt, :sum2_T_sdpc_cs_jt, :beta_max, :beta_min, :q_max_p0, :q_min_p0, :t), Tuple{Int64, Int64, Int64, Int64, Int64, Int64, Float64, Float64, Float64, Float64, Int64}}}() : [ (;j = parse(Int, match(r"\d+", val["uid"]).match) + L_J_br + 1, - j_prcs = parse(Int, match(r"\d+", val["uid"]).match) + 1, - j_cs = parse(Int, match(r"\d+", val["uid"]).match) + 1 - L_J_pr, + j_cs=get_j_cs(val["uid"], L_J_pr, L_J_cs, producers_first), + j_prcs=get_j_prcs(val["uid"], L_J_pr, L_J_cs, producers_first), u_on = uc["on_status"][t], sum2_T_supc_cs_jt=sum2_T_supc_cs[p.j_cs, t], sum2_T_sdpc_cs_jt=sum2_T_sdpc_cs[p.j_cs, t], @@ -904,15 +954,15 @@ function parse_sc_data(data, uc_data, data_json) t=t ) for val in values(data.sdd_lookup), t in periods - if parse(Int, match(r"\d+", val["uid"]).match) >= L_J_pr && val["q_bound_cap"]==1 + if !is_pr(val["uid"], L_J_pr, L_J_cs, producers_first) && val["q_bound_cap"]==1 for uc in uc_data["time_series_output"]["simple_dispatchable_device"] if p.uid == uc["uid"]], - csarray_pqe = isempty(val for val in values(data.sdd_lookup) if parse(Int, match(r"\d+", val["uid"]).match) < L_J_pr && val["q_bound_cap"]==1) ? + csarray_pqe = isempty(val for val in values(data.sdd_lookup) if !is_pr(val["uid"], L_J_pr, L_J_cs, producers_first) && val["q_bound_cap"]==1) ? empty_data = Vector{NamedTuple{(:j, :jprcs, :j_cs, :u_on, :sum2_T_supc_cs_jt, :sum2_T_sdpc_cs_jt, :beta, :q_p0, :t), Tuple{Int64, Int64, Int64, Int64, Int64, Int64, Float64, Float64, Int64}}}() : [ (;j = parse(Int, match(r"\d+", val["uid"]).match) + L_J_br + 1, - j_prcs = parse(Int, match(r"\d+", val["uid"]).match) + 1, - j_cs = parse(Int, match(r"\d+", val["uid"]).match) + 1 - L_J_pr, + j_cs=get_j_cs(val["uid"], L_J_pr, L_J_cs, producers_first), + j_prcs=get_j_prcs(val["uid"], L_J_pr, L_J_cs, producers_first), u_on = uc["on_status"][t], sum2_T_supc_cs_jt=sum2_T_supc_cs[p.j_cs, t], sum2_T_sdpc_cs_jt=sum2_T_sdpc_cs[p.j_cs, t], @@ -921,7 +971,7 @@ function parse_sc_data(data, uc_data, data_json) t=t ) for val in values(data.sdd_lookup), t in periods - if parse(Int, match(r"\d+", val["uid"]).match) >= L_J_pr && val["q_linear_cap"]==1 + if !is_pr(val["uid"], L_J_pr, L_J_cs, producers_first) && val["q_linear_cap"]==1 for uc in uc_data["time_series_output"]["simple_dispatchable_device"] if p.uid == uc["uid"]], @@ -984,19 +1034,22 @@ function parse_sc_data(data, uc_data, data_json) lengths = (L_J_xf, L_J_ln, L_J_ac, L_J_dc, L_J_br, L_J_cs, L_J_pr, L_J_cspr, L_J_sh, I, L_T, L_N_p, L_N_q, L_W_en_min_pr, L_W_en_min_cs, L_W_en_max_pr, L_W_en_max_cs, K) - return sc_time_data, lengths + return sc_time_data, lengths, producers_first end -function save_go3_solution(uc_filename, solution_name, result, vars, lengths) +function save_go3_solution(uc_filename, solution_name, result, vars, lengths, producers_first) uc_data = JSON.parsefile(uc_filename) (L_J_xf, L_J_ln, L_J_ac, L_J_dc, L_J_br, L_J_cs, L_J_pr, L_J_cspr, L_J_sh, I, L_T, L_N_p, L_N_q, L_W_en_min_pr, L_W_en_min_cs, L_W_en_max_pr, L_W_en_max_cs, K) = lengths #Update simple dispatchable devices for line in uc_data["time_series_output"]["simple_dispatchable_device"] - solution_index = parse(Int, match(r"\d+", line["uid"]).match) + 1 #This corresponds to j_prcs - if solution_index > L_J_pr + raw_uid = parse(Int, match(r"\d+", line["uid"]).match) + solution_index = raw_uid + 1 #This corresponds to j_prcs + if !is_pr(raw_uid, L_J_pr, L_J_cs, producers_first) #This section corresponds to consuming devices - solution_index -= L_J_pr + if producers_first + solution_index -= L_J_pr + end line["p_syn_res"] = Array(solution(result, vars.p_jt_scr_cs))[solution_index,:] line["p_ramp_res_up_online"] = Array(solution(result, vars.p_jt_rru_on_cs))[solution_index,:] line["p_nsyn_res"] = zeros(L_T) @@ -1011,6 +1064,9 @@ function save_go3_solution(uc_filename, solution_name, result, vars, lengths) line["p_ramp_res_down_offline"] = Array(solution(result, vars.p_jt_rrd_off_cs))[solution_index,:] else #Producing devices + if !producers_first + solution_index -= L_J_cs + end line["p_syn_res"] = Array(solution(result, vars.p_jt_scr_pr))[solution_index,:] line["p_ramp_res_up_online"] = Array(solution(result, vars.p_jt_rru_on_pr))[solution_index,:] line["p_nsyn_res"] = Array(solution(result, vars.p_jt_nsc_pr))[solution_index,:] diff --git a/src/scopf.jl b/src/scopf.jl index 06db084..9f5fb89 100644 --- a/src/scopf.jl +++ b/src/scopf.jl @@ -13,7 +13,8 @@ function scopf_model( uc_data = JSON.parsefile(uc_filename) data = GOC3Benchmark.get_data_from_file(filename) data_json = JSON.parsefile(filename) - sc_data, lengths = parse_sc_data(data, uc_data, data_json) + sc_data, lengths, producers_first = parse_sc_data(data, uc_data, data_json) + @info "parsed data" (L_J_xf, L_J_ln, L_J_ac, L_J_dc, L_J_br, L_J_cs, L_J_pr, L_J_cspr, L_J_sh, I, L_T, L_N_p, L_N_q, L_W_en_min_pr, @@ -674,6 +675,7 @@ function scopf_model( end + @info "built model" return model, sc_data_array, vars, lengths end diff --git a/test/opf_tests.jl b/test/opf_tests.jl index a163055..024db46 100644 --- a/test/opf_tests.jl +++ b/test/opf_tests.jl @@ -1,48 +1,49 @@ -function test_case3(result, result_pm, result_nlp_pm, dicts, pg, qg, p, q) - test_static_case(result, result_pm, result_nlp_pm, dicts, pg, qg) +function test_case3(result, result_pm, result_nlp_pm, pg, qg, p, q) + test_static_case(result, result_pm, result_nlp_pm, pg, qg) #Branches are encoded differently in solutions, so matches are hard coded vars_dict = Dict("p" => p, "q" => q) for st_var in ["p", "q"] var = vars_dict[st_var] - @test isapprox(Array(solution(result, var))[1], result_pm["solution"]["branch"]["2"][string(st_var, "f")], atol = result.options.tol*100) - @test isapprox(Array(solution(result, var))[2], result_pm["solution"]["branch"]["3"][string(st_var, "f")], atol = result.options.tol*100) - @test isapprox(Array(solution(result, var))[3], result_pm["solution"]["branch"]["1"][string(st_var, "f")], atol = result.options.tol*100) - @test isapprox(Array(solution(result, var))[4], result_pm["solution"]["branch"]["2"][string(st_var, "t")], atol = result.options.tol*100) - @test isapprox(Array(solution(result, var))[5], result_pm["solution"]["branch"]["3"][string(st_var, "t")], atol = result.options.tol*100) - @test isapprox(Array(solution(result, var))[6], result_pm["solution"]["branch"]["1"][string(st_var, "t")], atol = result.options.tol*100) + for i in 1:length(result_pm["solution"]["branch"]) + @test isapprox(Array(solution(result, var))[i], result_pm["solution"]["branch"][string(i)][string(st_var, "f")], atol = result.options.tol*100) + end end end -function test_static_case(result, result_pm, result_nlp_pm, dicts, pg, qg) +function test_static_case(result, result_pm, result_nlp_pm, pg, qg) @test result.status == result_nlp_pm.status @test isapprox(result.objective, result_pm["objective"], rtol = result.options.tol*100) - for key in keys(dicts.gen) - @test isapprox(Array(solution(result, pg))[dicts.gen[key]], result_pm["solution"]["gen"][string(key)]["pg"], atol = result.options.tol*1000) - @test isapprox(Array(solution(result, qg))[dicts.gen[key]], result_pm["solution"]["gen"][string(key)]["qg"], atol = result.options.tol*1000) + for i in 1:length(result_pm["solution"]["gen"]) + @test isapprox(Array(solution(result, pg))[i], result_pm["solution"]["gen"][string(i)]["pg"], atol = result.options.tol*1000) + @test isapprox(Array(solution(result, qg))[i], result_pm["solution"]["gen"][string(i)]["qg"], atol = result.options.tol*1000) end end -function test_polar_voltage(result, result_pm, dicts, va, vm) - for key in keys(dicts.bus) - @test isapprox(Array(solution(result, va))[dicts.bus[key]], result_pm["solution"]["bus"][string(key)]["va"], atol = result.options.tol*100) - @test isapprox(Array(solution(result, vm))[dicts.bus[key]], result_pm["solution"]["bus"][string(key)]["vm"], rtol = result.options.tol*100) +function test_polar_voltage(result, result_pm, va, vm) + @info "ppolar" + @info result_pm + for i in 1:length(result_pm["solution"]["bus"]) + @test isapprox(Array(solution(result, va))[i], result_pm["solution"]["bus"][string(i)]["va"], atol = result.options.tol*100) + @test isapprox(Array(solution(result, vm))[i], result_pm["solution"]["bus"][string(i)]["vm"], rtol = result.options.tol*100) end end -function test_rect_voltage(result, result_pm, dicts, vr, vim) - for key in keys(dicts.bus) - @test isapprox(Array(solution(result, vr))[dicts.bus[key]], result_pm["solution"]["bus"][string(key)]["vr"], rtol = result.options.tol*100) - @test isapprox(Array(solution(result, vim))[dicts.bus[key]], result_pm["solution"]["bus"][string(key)]["vi"], atol = result.options.tol*100) +function test_rect_voltage(result, result_pm, vr, vim) + @info "rrect" + @info result_pm + for i in 1:length(result_pm["solution"]["bus"]) + @test isapprox(Array(solution(result, vr))[i], result_pm["solution"]["bus"][string(i)]["vr"], rtol = result.options.tol*100) + @test isapprox(Array(solution(result, vim))[i], result_pm["solution"]["bus"][string(i)]["vi"], atol = result.options.tol*100) end end -function test_case5(result, result_pm, result_nlp_pm, dicts, pg, qg, p, q) - test_static_case(result, result_pm, result_nlp_pm, dicts, pg, qg) +function test_case5(result, result_pm, result_nlp_pm, pg, qg, p, q) + test_static_case(result, result_pm, result_nlp_pm, pg, qg) end -function test_case14(result, result_pm, result_nlp_pm, dicts, pg, qg, p, q) - test_static_case(result, result_pm, result_nlp_pm, dicts, pg, qg) +function test_case14(result, result_pm, result_nlp_pm, pg, qg, p, q) + test_static_case(result, result_pm, result_nlp_pm, pg, qg) end function test_float32(m, m64, result, backend) @@ -63,4 +64,4 @@ end function test_mp_case(result, true_sol) @test result.status == MadNLP.SOLVE_SUCCEEDED || result.status == MadNLP.SOLVED_TO_ACCEPTABLE_LEVEL @test isapprox(result.objective, true_sol, rtol = result.options.tol*100) -end \ No newline at end of file +end diff --git a/test/runtests.jl b/test/runtests.jl index 699fa81..e06001e 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -3,10 +3,10 @@ using Test, ExaModelsPower, MadNLP, MadNLPGPU, KernelAbstractions, CUDA, PowerMo include("opf_tests.jl") const CONFIGS = [ - (Float32, nothing), (Float64, nothing), - (Float32, CPU()), (Float64, CPU()), + (Float32, nothing), + (Float32, CPU()), ] if CUDA.has_cuda_gpu() @@ -48,180 +48,149 @@ mp_stor_test_cases = [("../data/pglib_opf_case3_lmbd_mod.m", "case3", "../data/c ("../data/pglib_opf_case5_pjm_mod.m", "case5", "../data/case5_5split.Pd", "../data/case5_5split.Qd", true_sol_case5_curve_stor, true_sol_case5_curve_stor_func, true_sol_case5_pregen_stor, true_sol_case5_pregen_stor_func)] +static_forms = [("rect", :rect, ACRPowerModel, test_rect_voltage), + ("polar", :polar, ACPPowerModel, test_polar_voltage)] + function example_func(d, srating) return d + 20/srating*d^2 end -function runtests() - @testset "ExaModelsPower test" begin - - for (T, backend) in CONFIGS +function sc_tests(filename) + uc_filename = "$filename.pop_solution.json" + model, sc_data_array, vars, lengths = ExaModelsPower.scopf_model(filename, uc_filename; backend=CUDABackend()) + @info "built model" + result = madnlp(model; print_level = MadNLP.ERROR, tol=8e-3, linear_solver=MadNLPGPU.CUDSSSolver) + JLD2.save("result.jld2", "solution", result, "vars", vars, "lens", lens) + ExaModelsPower.save_go3_solution(uc_filename, "solution_go3", result, vars, lengths) +end - for (filename, case, test_function) in test_cases - #Test static opf - data, dicts = ExaModelsPower.parse_ac_power_data(filename) - data_pm = PowerModels.parse_file(filename) +PowerModels.silence() - #Polar tests - m, v, c = opf_model(filename; T=T, backend = backend) - result = madnlp(m; print_level = MadNLP.ERROR) - va, vm, pg, qg, p, q = v +function parse_pm(filename) + data = PowerModels.parse_file(filename) + PowerModels.standardize_cost_terms!(data, order = 2) + PowerModels.calc_thermal_limits!(data) - nlp_solver = JuMP.optimizer_with_attributes(Ipopt.Optimizer, "tol"=>Float64(result.options.tol), "print_level"=>0) - result_pm = solve_opf(filename,ACPPowerModel, nlp_solver) + return data +end - - m_pm = JuMP.Model() - pm = instantiate_model(data_pm, ACPPowerModel, PowerModels.build_opf, jump_model = m_pm) - nlp_pm = MathOptNLPModel(m_pm) - result_nlp_pm = madnlp(nlp_pm; print_level = MadNLP.ERROR) - - - @testset "$(case), static, $(T), $(backend), polar" begin - if T == Float32 - m64, v64, c64 = opf_model(filename; T=Float64, backend = backend) - result = madnlp(m64; print_level = MadNLP.ERROR) - test_float32(m, m64, result, backend) - else - eval(test_function)(result, result_pm, result_nlp_pm, dicts, pg, qg, p, q) - test_polar_voltage(result, result_pm, dicts, va, vm) - end - end +function runtests() + @testset "ExaModelsPower test" begin - #Rectangular tests - m, v, c = eval(opf_model)(filename; T=T, backend = backend, form = :rect) - result = madnlp(m; print_level = MadNLP.ERROR) - vr, vim, pg, qg, p, q = v - - nlp_solver = JuMP.optimizer_with_attributes(Ipopt.Optimizer, "tol"=>Float64(result.options.tol), "print_level"=>0) - result_pm = solve_opf(filename, ACRPowerModel, nlp_solver) - - m_pm = JuMP.Model() - pm = instantiate_model(data_pm, ACRPowerModel, PowerModels.build_opf, jump_model = m_pm) - nlp_pm = MathOptNLPModel(m_pm) - result_nlp_pm = madnlp(nlp_pm; print_level = MadNLP.ERROR) - - @testset "$(case), static, $(T), $(backend), rect" begin - if T == Float32 - m64, v64, c64 = opf_model(filename; T=Float64, backend = backend, form = :rect) - result = madnlp(m64; print_level = MadNLP.ERROR) - test_float32(m, m64, result, backend) - else - eval(test_function)(result, result_pm, result_nlp_pm, dicts, pg, qg, p, q) - test_rect_voltage(result, result_pm, dicts, vr, vim) + for (T, backend) in CONFIGS + for (filename, case, test_function) in test_cases + data_pm = parse_pm(filename) + for (form_str, form, power_model, test_voltage) in static_forms + m32, v32, c32 = opf_model(filename; T=Float32, backend = backend, form=form) + result32 = madnlp(m32; print_level = MadNLP.ERROR) + va32, vm32, pg32, qg32, p32, q32 = v32 + + m64, v64, c64 = opf_model(filename; T=Float64, backend = backend, form=form) + result64 = madnlp(m64; print_level = MadNLP.ERROR) + va64, vm64, pg64, qg64, p64, q64 = v64 + + nlp_solver = JuMP.optimizer_with_attributes(Ipopt.Optimizer, "tol"=>Float64(result64.options.tol), "print_level"=>0) + result_pm = solve_opf(filename, power_model, nlp_solver) + + m_pm = JuMP.Model() + pm = instantiate_model(data_pm, power_model, PowerModels.build_opf, jump_model = m_pm) + nlp_pm = MathOptNLPModel(m_pm) + result_nlp_pm = madnlp(nlp_pm; print_level = MadNLP.ERROR) + + @info form_str + @testset "$case, static, $backend, $form_str" begin + test_float32(m32, m64, result64, backend) + eval(test_function)(result64, result_pm, result_nlp_pm, pg64, qg64, p64, q64) + test_voltage(result64, result_pm, va64, vm64) end end end #Test MP - for (form_str, symbol) in [("polar", :polar), ("rect", :rect)] + for (form_str, symbol) in [("rect", :rect), ("polar", :polar)] for (filename, case, Pd_pregen, Qd_pregen, true_sol_curve, true_sol_pregen) in mp_test_cases #Curve = [1, .9, .8, .95, 1] - m, v, c = eval(mpopf_model)(filename, [1, .9, .8, .95, 1]; T = T, backend = backend, form = symbol) - result = madnlp(m; print_level = MadNLP.ERROR) + m32, v32, c32 = eval(mpopf_model)(filename, [1, .9, .8, .95, 1]; T = Float32, backend = backend, form = symbol) + result32 = madnlp(m32; print_level = MadNLP.ERROR) + m64, v64, c64 = eval(mpopf_model)(filename, [1, .9, .8, .95, 1]; T = Float64, backend = backend, form = symbol) + result64 = madnlp(m64; print_level = MadNLP.ERROR) @testset "$(case), MP, $(T), $(backend), curve, $(form_str)" begin - if T == Float32 - m64, v64, c64 = eval(mpopf_model)(filename, [1, .9, .8, .95, 1]; T=Float64, backend = backend, form = symbol) - result = madnlp(m64; print_level = MadNLP.ERROR) - test_float32(m, m64, result, backend) - else - test_mp_case(result, true_sol_curve) - end + test_float32(m32, m64, result64, backend) + test_mp_case(result64, true_sol_curve) end #w function - m, v, c = eval(mpopf_model)(filename, [1, .9, .8, .95, 1], example_func; T = T, backend = backend, form = symbol) - result = madnlp(m; print_level = MadNLP.ERROR) + m32, v32, c32 = eval(mpopf_model)(filename, [1, .9, .8, .95, 1], example_func; T = Float32, backend = backend, form = symbol) + result32 = madnlp(m32; print_level = MadNLP.ERROR) + m64, v64, c64 = eval(mpopf_model)(filename, [1, .9, .8, .95, 1], example_func; T = Float64, backend = backend, form = symbol) + result64 = madnlp(m64; print_level = MadNLP.ERROR) @testset "$(case), MP, $(T), $(backend), curve, $(form_str), func" begin - if T == Float32 - m64, v64, c64 = eval(mpopf_model)(filename, [1, .9, .8, .95, 1], example_func; T=Float64, backend = backend, form = symbol) - result = madnlp(m64; print_level = MadNLP.ERROR) - test_float32(m, m64, result, backend) - else - test_mp_case(result, true_sol_curve) - end + test_float32(m32, m64, result64, backend) + test_mp_case(result64, true_sol_curve) end #Pregenerated Pd and Qd - m, v, c = eval(mpopf_model)(filename, Pd_pregen, Qd_pregen; T = T, backend = backend, form = symbol) - result = madnlp(m; print_level = MadNLP.ERROR) + m32, v32, c32 = eval(mpopf_model)(filename, Pd_pregen, Qd_pregen; T = Float32, backend = backend, form = symbol) + result32 = madnlp(m32; print_level = MadNLP.ERROR) + m64, v64, c64 = eval(mpopf_model)(filename, Pd_pregen, Qd_pregen; T = Float64, backend = backend, form = symbol) + result64 = madnlp(m64; print_level = MadNLP.ERROR) @testset "$(case), MP, $(T), $(backend), pregen, $(form_str)" begin - if T == Float32 - m64, v64, c64 = eval(mpopf_model)(filename, Pd_pregen, Qd_pregen; T=Float64, backend = backend, form = symbol) - result = madnlp(m64; print_level = MadNLP.ERROR) - test_float32(m, m64, result, backend) - else - test_mp_case(result, true_sol_pregen) - end + test_float32(m32, m64, result64, backend) + test_mp_case(result64, true_sol_pregen) end #w function - m, v, c = eval(mpopf_model)(filename, Pd_pregen, Qd_pregen, example_func; T = T, backend = backend, form = symbol) - result = madnlp(m; print_level = MadNLP.ERROR) + m32, v32, c32 = eval(mpopf_model)(filename, Pd_pregen, Qd_pregen, example_func; T = Float32, backend = backend, form = symbol) + result32 = madnlp(m32; print_level = MadNLP.ERROR) + m64, v64, c64 = eval(mpopf_model)(filename, Pd_pregen, Qd_pregen, example_func; T = Float64, backend = backend, form = symbol) + result64 = madnlp(m64; print_level = MadNLP.ERROR) @testset "$(case), MP, $(T), $(backend), pregen, $(form_str), func" begin - if T == Float32 - m64, v64, c64 = eval(mpopf_model)(filename, Pd_pregen, Qd_pregen, example_func; T=Float64, backend = backend, form = symbol) - result = madnlp(m64; print_level = MadNLP.ERROR) - test_float32(m, m64, result, backend) - else - test_mp_case(result, true_sol_pregen) - end + test_float32(m32, m64, result64, backend) + test_mp_case(result64, true_sol_pregen) end end - #Test MP w storage + # Test MP w storage for (filename, case, Pd_pregen, Qd_pregen, true_sol_curve_stor, true_sol_curve_stor_func, true_sol_pregen_stor, true_sol_pregen_stor_func) in mp_stor_test_cases - m, v, c = eval(mpopf_model)(filename, [1, .9, .8, .95, 1]; T = T, backend = backend, form = symbol) - result = madnlp(m; print_level = MadNLP.ERROR) + m32, v32, c32 = eval(mpopf_model)(filename, [1, .9, .8, .95, 1]; T = Float32, backend = backend, form = symbol) + result32 = madnlp(m32; print_level = MadNLP.ERROR) + m64, v64, c64 = eval(mpopf_model)(filename, [1, .9, .8, .95, 1]; T = Float64, backend = backend, form = symbol) + result64 = madnlp(m64; print_level = MadNLP.ERROR) @testset "MP w storage, $(case), $(T), $(backend), curve, $(form_str)" begin - if T == Float32 - m64, v64, c64 = eval(mpopf_model)(filename, [1, .9, .8, .95, 1]; T=Float64, backend = backend, form = symbol) - result = madnlp(m64; print_level = MadNLP.ERROR) - test_float32(m, m64, result, backend) - else - test_mp_case(result, true_sol_curve_stor) - end + test_float32(m32, m64, result64, backend) + test_mp_case(result64, true_sol_curve_stor) end #With function - m, v, c = eval(mpopf_model)(filename, [1, .9, .8, .95, 1], example_func; T = T, backend = backend, form = symbol) - result = madnlp(m; print_level = MadNLP.ERROR) + m32, v32, c32 = eval(mpopf_model)(filename, [1, .9, .8, .95, 1], example_func; T = Float32, backend = backend, form = symbol) + result32 = madnlp(m32; print_level = MadNLP.ERROR) + m64, v64, c64 = eval(mpopf_model)(filename, [1, .9, .8, .95, 1], example_func; T = Float64, backend = backend, form = symbol) + result64 = madnlp(m64; print_level = MadNLP.ERROR) @testset "MP w storage, $(case), $(T), $(backend), curve, $(form_str), func" begin - if T == Float32 - m64, v64, c64 = eval(mpopf_model)(filename, [1, .9, .8, .95, 1], example_func; T=Float64, backend = backend, form = symbol) - result = madnlp(m64; print_level = MadNLP.ERROR) - test_float32(m, m64, result, backend) - else - test_mp_case(result, true_sol_curve_stor_func) - end + test_float32(m32, m64, result64, backend) + test_mp_case(result64, true_sol_curve_stor_func) end #Pregenerated Pd and Qd - m, v, c = eval(mpopf_model)(filename, Pd_pregen, Qd_pregen; T = T, backend = backend, form = symbol) - result = madnlp(m; print_level = MadNLP.ERROR) + m32, v32, c32 = eval(mpopf_model)(filename, Pd_pregen, Qd_pregen; T = Float32, backend = backend, form = symbol) + result32 = madnlp(m32; print_level = MadNLP.ERROR) + m64, v64, c64 = eval(mpopf_model)(filename, Pd_pregen, Qd_pregen; T = Float64, backend = backend, form = symbol) + result64 = madnlp(m64; print_level = MadNLP.ERROR) @testset "MP w storage, $(case), $(T), $(backend), pregen, $(form_str)" begin - if T == Float32 - m64, v64, c64 = eval(mpopf_model)(filename, Pd_pregen, Qd_pregen; T=Float64, backend = backend, form = symbol) - result = madnlp(m64; print_level = MadNLP.ERROR) - test_float32(m, m64, result, backend) - else - test_mp_case(result, true_sol_pregen_stor) - end + test_float32(m32, m64, result64, backend) + test_mp_case(result64, true_sol_pregen_stor) end #With function - m, v, c = eval(mpopf_model)(filename, Pd_pregen, Qd_pregen, example_func; T = T, backend = backend, form = symbol) - result = madnlp(m; print_level = MadNLP.ERROR) + m32, v32, c32 = eval(mpopf_model)(filename, Pd_pregen, Qd_pregen, example_func; T = Float32, backend = backend, form = symbol) + result32 = madnlp(m32; print_level = MadNLP.ERROR) + m64, v64, c64 = eval(mpopf_model)(filename, Pd_pregen, Qd_pregen, example_func; T = Float64, backend = backend, form = symbol) + result64 = madnlp(m64; print_level = MadNLP.ERROR) @testset "MP w storage, $(case), $(T), $(backend), pregen, $(form_str), func" begin - if T == Float32 - m64, v64, c64 = eval(mpopf_model)(filename, Pd_pregen, Qd_pregen, example_func; T=Float64, backend = backend, form = symbol) - result = madnlp(m64; print_level = MadNLP.ERROR) - test_float32(m, m64, result, backend) - else - test_mp_case(result, true_sol_pregen_stor_func) - end + test_float32(m32, m64, result64, backend) + test_mp_case(result64, true_sol_pregen_stor_func) end end end