diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 987b69b..6804013 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -34,9 +34,10 @@ jobs: key: compiled-${{ runner.os }}-${{ hashFiles('**/Project.toml') }}-${{ hashFiles('**/Manifest.toml') }} restore-keys: compiled-${{ runner.os }}- + - name: Install package dependencies run: | - julia --project=. -e 'using Pkg; Pkg.instantiate(); Pkg.precompile()' + julia --project=. -e 'using Pkg; Pkg.add(PackageSpec(url="https://github.com/lanl-ansi/GOC3Benchmark.jl")); Pkg.instantiate(); Pkg.precompile()' - name: Check for CUDA availability id: check-cuda diff --git a/.github/workflows/docs.yml b/.github/workflows/docs.yml index b586b09..d2fc95a 100644 --- a/.github/workflows/docs.yml +++ b/.github/workflows/docs.yml @@ -17,7 +17,7 @@ jobs: with: version: '1' - name: Install dependencies - run: julia --project=docs/ -e 'using Pkg; Pkg.Registry.update(); Pkg.develop(PackageSpec(path=pwd())); Pkg.instantiate()' + run: julia --project=docs/ -e 'using Pkg; Pkg.add(PackageSpec(url="https://github.com/lanl-ansi/GOC3Benchmark.jl")); Pkg.Registry.update(); Pkg.develop(PackageSpec(path=pwd())); Pkg.instantiate(); Pkg.precompile()' - name: Build and deploy env: GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} # If authenticating with GitHub Actions token diff --git a/Project.toml b/Project.toml index a7d6947..0b83bf9 100644 --- a/Project.toml +++ b/Project.toml @@ -9,6 +9,9 @@ PGLib = "07a8691f-3d11-4330-951b-3c50f98338be" PowerModels = "c36e90e8-916a-50a6-bd94-075b64ef4655" ExaModels = "1037b233-b668-4ce9-9b63-f9f681f55dd2" JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819" +JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6" +GOC3Benchmark = "3a45b339-860d-44d0-a64b-5f943cdd120b" +Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6" [compat] ExaModels = "~0.8" diff --git a/README.md b/README.md index ed002b8..3a94602 100644 --- a/README.md +++ b/README.md @@ -9,7 +9,8 @@ ExaModelsPower.jl is an optimal power flow models using ExaModels.jl ## Usage ### Static optimal power flow ```julia -using ExaModelsPower, MadNLP, MadNLPGPU, CUDA +using ExaModelsPower, MadNLP, MadNLPGPU, CUDA, ExaModels, GOC3Benchmark, JSON + model, vars, cons = opf_model( "pglib_opf_case118_ieee.m"; @@ -21,11 +22,32 @@ result = madnlp(model; tol=1e-6) ### Security-constrained optimal power flow ```julia -model, vars, cons = scopf_model( - "pglib_opf_case118_ieee.m"; contingencies = [1,2], +#This model is based on the GOC3 formulation of the SCOPF problem +#https://www.pnnl.gov/publications/grid-optimization-competition-challenge-3-problem-formulation + +#The current implementation requires a UC solution to be provided, which is then parsed with +#the other input data to generate a structure of named tuples which can then interface with +#ExaModels to generate the full model. We do not make any relaxations or decompositions for this problem + +model, sc_data, vars, lengths = scopf_model( + "data/C3E4N00073D1_scenario_303.json", "data/C3E4N00073D1_scenario_303_solution.json"; backend = CUDABackend() ) -result = madnlp(model; tol=1e-6) # currently failing +result = madnlp(model; tol=1e-4) + +#Solution from GPU can be used to warm start a CPU solution or vice versa +model_cpu, sc_data, vars, lengths = scopf_model( + "data/C3E4N00073D1_scenario_303.json", "data/C3E4N00073D1_scenario_303_solution.json"; + result_set = [result, vars] +) +result_cpu = ipopt(model_cpu; tol=1e-8) + +#Additionally, the SC problem can be evaluated without contingencies +model, sc_data, vars, lengths = scopf_model( + "data/C3E4N00073D1_scenario_303.json", "data/C3E4N00073D1_scenario_303_solution.json"; + backend = CUDABackend(), include_ctg = false +) +result = madnlp(model; tol=1e-4) ``` ### Multi-period optimal power flow @@ -59,7 +81,7 @@ result = madnlp(model; tol=1e-6) #Alternatively, provide a smooth function for the charge/discharge efficiency to remove complementarity constraint function example_func(d, srating) - return d + .2/srating*d^2 + return -((s_rating/2)^d)+1 end model, vars, cons = mpopf_model( diff --git a/src/ExaModelsPower.jl b/src/ExaModelsPower.jl index 3aafdaa..eaecd71 100644 --- a/src/ExaModelsPower.jl +++ b/src/ExaModelsPower.jl @@ -1,15 +1,18 @@ module ExaModelsPower import JLD2 +import Downloads +import ExaModels: ExaCore, variable, constraint, ExaModel, objective, constraint!, convert_array, solution import PGLib -import ExaModels: ExaCore, variable, constraint, ExaModel, objective, constraint!, convert_array import PowerModels + include("parser.jl") include("opf.jl") include("scopf.jl") include("mpopf.jl") include("constraint.jl") +include("sc_parser.jl") const NAMES = filter(names(@__MODULE__; all = true)) do x str = string(x) diff --git a/src/mpopf.jl b/src/mpopf.jl index b5f7a16..f525469 100644 --- a/src/mpopf.jl +++ b/src/mpopf.jl @@ -58,17 +58,8 @@ function update_load_data(busarray, pd, qd, baseMVA, busdict) end end -#If no storage contraints, the "build_base_xxx_mpopf" returns the final version of the mpopf -function build_base_polar_mpopf(core, data, N, Nbus) - #voltage angle, voltage magnitude - va = variable(core, Nbus, N;) - vm = variable( - core, - Nbus, N; - start = ones(size(data.busarray)), - lvar = repeat(data.vmin, 1, N), - uvar = repeat(data.vmax, 1, N), - ) +#If no storage contraints, the "build_base_mpopf" returns the final version of the mpopf +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)) @@ -78,55 +69,8 @@ function build_base_polar_mpopf(core, data, N, Nbus) 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)) - #Storage specific variables - - #active/reactive power from bus into storage - pst = variable(core, size(data.storage, 1), N) - qst = variable(core, size(data.storage, 1), N) - - #current magnitude squared - I2 = variable(core, size(data.storage, 1), N; lvar = zeros(size(data.storarray))) - - #ability of converter to control generation/absorption of reactive power - qint = variable(core, size(data.storage, 1), N; lvar = -repeat(data.srating, 1, N), uvar = repeat(data.srating, 1, N)) - - #energy/ state of charge - E = variable(core, size(data.storage, 1), N; lvar = zeros(size(data.storarray)), uvar = repeat(data.emax, 1, N)) - o = objective(core, gen_cost(g, pg[g.i,g.t]) for g in data.genarray) - 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_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), - ) - - 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_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_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_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_from_thermal_limit = constraint( core, @@ -148,436 +92,260 @@ function build_base_polar_mpopf(core, data, N, Nbus) ) cons = ( - c_ref_angle = c_ref_angle, - c_to_active_power_flow = c_to_active_power_flow, - c_to_reactive_power_flow = c_to_reactive_power_flow, - c_from_active_power_flow = c_from_active_power_flow, - 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_from_thermal_limit = c_from_thermal_limit, c_to_thermal_limit = c_to_thermal_limit, c_ramp_rate = c_ramp_rate ) vars = ( - va = va, - vm = vm, pg = pg, qg = qg, p = p, - q = q, - pst = pst, - qst = qst, - I2 = I2, - qint = qint, - E = E + q = q, ) return vars, cons end -function build_base_rect_mpopf(core, data, N, Nbus) - #real, imaginary voltage - vr = variable(core, Nbus, N; start = ones(size(data.busarray))) - vim = variable(core, Nbus, 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)) - - #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)) - - #Storage specific variables - - #active/reactive power from bus into storage - pst = variable(core, size(data.storage, 1), N) - qst = variable(core, size(data.storage, 1), N) +function add_mpopf_cons(core, data, N, Nbus, vars, cons, form) + pg, qg, p, q = vars + if form == :polar + #voltage angle, voltage magnitude + va = variable(core, Nbus, N; lvar = -pi, uvar = pi) + vm = variable( + core, + Nbus, N; + start = ones(size(data.busarray)), + lvar = repeat(data.vmin, 1, N), + uvar = repeat(data.vmax, 1, N), + ) - #current magnitude squared - I2 = variable(core, size(data.storage, 1), N; lvar = zeros(size(data.storarray))) + 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) - #ability of converter to control generation/absorption of reactive power - qint = variable(core, size(data.storage, 1), N; lvar = -repeat(data.srating, 1, N), uvar = repeat(data.srating, 1, N)) + 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_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), + ) + + 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) + + + cons = (;cons..., + c_ref_angle = c_ref_angle, + c_to_active_power_flow = c_to_active_power_flow, + c_to_reactive_power_flow = c_to_reactive_power_flow, + c_from_active_power_flow = c_from_active_power_flow, + 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) + vars = (;vars..., va = va, vm = vm) - #energy/ state of charge - E = variable(core, size(data.storage, 1), N; lvar = zeros(size(data.storarray)), uvar = repeat(data.emax, 1, N)) + elseif form == :rect + #real, imaginary voltage + vr = variable(core, Nbus, N; start = ones(size(data.busarray))) + vim = variable(core, Nbus, N;) - o = objective(core, gen_cost(g, pg[g.i,g.t]) for g in data.genarray) + 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_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_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_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; + 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_voltage_magnitude = constraint( + core, c_voltage_magnitude_rect(vr[b.i, b.t], vim[b.i, b.t]) + for b in data.busarray; + lcon = repeat(data.vmin, 1, N).^2, + ucon = repeat(data.vmax, 1, N).^2 + ) + cons = (;cons..., + c_ref_angle = c_ref_angle, + c_to_active_power_flow = c_to_active_power_flow, + c_to_reactive_power_flow = c_to_reactive_power_flow, + c_from_active_power_flow = c_from_active_power_flow, + 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_voltage_magnitude = c_voltage_magnitude) + vars = (;vars..., vr = vr, vim = vim) + + end - 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_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; - 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_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_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_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_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; - 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; - lcon = fill(-Inf, size(data.barray)) - ) - - c_voltage_magnitude = constraint( - core, c_voltage_magnitude_rect(vr[b.i, b.t], vim[b.i, b.t]) - for b in data.busarray; - lcon = repeat(data.vmin, 1, N).^2, - ucon = repeat(data.vmax, 1, N).^2 - ) - - 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]; - lcon = repeat(-data.Δp, 1, N-1), - ucon = repeat( data.Δp, 1, N-1), - ) - - vars = ( - vr = vr, - vim = vim, - pg = pg, - qg = qg, - p = p, - q = q, - pst = pst, - qst = qst, - I2 = I2, - qint = qint, - E = E - ) - - cons = cons = ( - c_ref_angle = c_ref_angle, - c_to_active_power_flow = c_to_active_power_flow, - c_to_reactive_power_flow = c_to_reactive_power_flow, - c_from_active_power_flow = c_from_active_power_flow, - 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_from_thermal_limit = c_from_thermal_limit, - c_to_thermal_limit = c_to_thermal_limit, - c_voltage_magnitude = c_voltage_magnitude, - c_ramp_rate = c_ramp_rate - ) - return vars, cons end - -function build_polar_mpopf(data, Nbus, N; backend = nothing, T = Float64, storage_complementarity_constraint = false, kwargs...) +function build_mpopf(data, Nbus, N, form; backend = nothing, T = Float64, storage_complementarity_constraint = false, kwargs...) core = ExaCore(T; backend = backend) - vars, cons = build_base_polar_mpopf(core, data, N, Nbus) + vars, cons = build_base_mpopf(core, data, N) + vars, cons = add_mpopf_cons(core, data, N, Nbus, vars, cons, form) if length(data.storarray) > 0 - va, vm, pg, qg, p, q, pst, qst, I2, qint, E = vars - - (c_ref_angle, - c_to_active_power_flow, - c_to_reactive_power_flow, - c_from_active_power_flow, - c_from_reactive_power_flow, - c_phase_angle_diff, - c_active_power_balance, - c_reactive_power_balance, - c_from_thermal_limit, - c_to_thermal_limit, - c_ramp_rate) = cons - - #charge or discharge from battery to grid - pstc = variable(core, size(data.storage, 1), N; lvar = zeros(size(data.storarray)), uvar = repeat(data.pcmax, 1, N)) - pstd = variable(core, size(data.storage, 1), N; lvar = zeros(size(data.storarray)), uvar = repeat(data.pdmax, 1, N)) - - #adding storage constraints - 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_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_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_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_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_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 = lcon = fill(-Inf, size(data.storarray))) - - 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)) - - #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) - end - - - cons = ( - c_ref_angle = c_ref_angle, - c_to_active_power_flow = c_to_active_power_flow, - c_to_reactive_power_flow = c_to_reactive_power_flow, - c_from_active_power_flow = c_from_active_power_flow, - 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_from_thermal_limit = c_from_thermal_limit, - c_to_thermal_limit = c_to_thermal_limit, - c_ramp_rate = c_ramp_rate, - c_active_storage_power = c_active_storage_power, - c_reactive_storage_power = c_reactive_storage_power, - c_ohms = c_ohms, - c_storage_state = c_storage_state, - c_storage_state_init = c_storage_state_init, - c_storage_transfer_thermal_limit = c_storage_transfer_thermal_limit, - c_discharge_thermal_limit = c_discharge_thermal_limit - ) - - vars = va, vm, pg, qg, p, q, pst, qst, I2, qint, E, pstd, pstc + 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 -function build_rect_mpopf(data, Nbus, N; backend = nothing, T = Float64, storage_complementarity_constraint = false, kwargs...) +#different constraints used when a function is added to remove complementarity and make charge/discharge curve smooth +function build_mpopf(data, Nbus, N, discharge_func::Function, form; backend = nothing, T = Float64, kwargs...) core = ExaCore(T; backend = backend) - vars, cons = build_base_rect_mpopf(core, data, N, Nbus) + vars, cons = build_base_mpopf(core, data, N) + vars, cons = add_mpopf_cons(core, data, N, Nbus, vars, cons, form) if length(data.storarray) > 0 - vr, vim, pg, qg, p, q, pst, qst, I2, qint, E = vars - - (c_ref_angle, - c_to_active_power_flow, - c_to_reactive_power_flow, - c_from_active_power_flow, - c_from_reactive_power_flow, - c_phase_angle_diff, - c_active_power_balance, - c_reactive_power_balance, - c_from_thermal_limit, - c_to_thermal_limit, - c_voltage_magnitude, - c_ramp_rate) = cons - - #charge or discharge from battery to grid - pstc = variable(core, size(data.storage, 1), N; lvar = zeros(size(data.storarray)), uvar = repeat(data.pcmax, 1, N)) - pstd = variable(core, size(data.storage, 1), N; lvar = zeros(size(data.storarray)), uvar = repeat(data.pdmax, 1, N)) - - #adding storage constraints - 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) + vars, cons = build_mpopf_stor_main(core, data, N, Nbus, vars, cons, form) + vars, cons = add_smooth_cons(core, data, N, vars, cons, discharge_func) + end - 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) + model = ExaModel(core; kwargs...) + return model, vars, cons +end - 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) +function build_mpopf_stor_main(core, data, N, Nbus, vars, cons, form) - 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]) + #Storage specific variables - 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]) + #active/reactive power from bus into storage + pst = variable(core, size(data.storage, 1), N) + qst = variable(core, size(data.storage, 1), N) - 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 = lcon = fill(-Inf, size(data.storarray))) + #current magnitude squared + I2 = variable(core, size(data.storage, 1), N; lvar = zeros(size(data.storarray))) - 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)) + #ability of converter to control generation/absorption of reactive power + qint = variable(core, size(data.storage, 1), N; lvar = -repeat(data.srating, 1, N), uvar = repeat(data.srating, 1, N)) - #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) - end + #energy/ state of charge + E = variable(core, size(data.storage, 1), N; lvar = zeros(size(data.storarray)), uvar = repeat(data.emax, 1, N)) - cons = ( - c_ref_angle = c_ref_angle, - c_to_active_power_flow = c_to_active_power_flow, - c_to_reactive_power_flow = c_to_reactive_power_flow, - c_from_active_power_flow = c_from_active_power_flow, - 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_from_thermal_limit = c_from_thermal_limit, - c_to_thermal_limit = c_to_thermal_limit, - c_voltage_magnitude = c_voltage_magnitude, - c_ramp_rate = c_ramp_rate, - c_active_storage_power = c_active_storage_power, - c_reactive_storage_power = c_reactive_storage_power, - c_ohms = c_ohms, - c_storage_state = c_storage_state, - c_storage_state_init = c_storage_state_init, - c_storage_transfer_thermal_limit = c_storage_transfer_thermal_limit, - c_discharge_thermal_limit = c_discharge_thermal_limit - ) - vars = vr, vim, pg, qg, p, q, pst, qst, I2, qint, E, pstc, pstd - end + #discharge from battery to grid + pstd = variable(core, size(data.storage, 1), N; uvar = repeat(data.pdmax, 1, N)) + vars = (;vars..., pst=pst, qst=qst, I2=I2, qint=qint, E=E, pstd=pstd) - model = ExaModel(core; kwargs...) - return model, vars, cons -end + c_active_power_balance = cons.c_active_power_balance + c_reactive_power_balance = cons.c_reactive_power_balance -#different constraints used when a function is added to remove complementarity and make charge/discharge curve smooth -function build_polar_mpopf(data, Nbus, N, discharge_func::Function; backend = nothing, T = Float64, kwargs...) - core = ExaCore(T; backend = backend) + 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) - vars, cons = build_base_polar_mpopf(core, data, N, Nbus) + 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) - if length(data.storarray) > 0 - va, vm, pg, qg, p, q, pst, qst, I2, qint, E = vars + 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_ref_angle, - c_to_active_power_flow, - c_to_reactive_power_flow, - c_from_active_power_flow, - c_from_reactive_power_flow, - c_phase_angle_diff, - c_active_power_balance, - c_reactive_power_balance, - c_from_thermal_limit, - c_to_thermal_limit, - c_ramp_rate) = cons + 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) + 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) + 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) + return vars, cons +end - #discharge or charge from battery to grid (positive or negative) - pstd = variable(core, size(data.storage, 1), N; lvar = zeros(size(data.storarray)), uvar = repeat(data.pdmax, 1, N)) +function add_piecewise_cons(core, data, N, vars, cons, storage_complementarity_constraint) + #charge from battery to grid + pstc = variable(core, size(data.storage, 1), N; lvar = zeros(size(data.storarray)), uvar = repeat(data.pcmax, 1, N)) + vars = (;vars..., pstc=pstc) - 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) + pst = vars.pst + pstd = vars.pstd + I2 = vars.I2 + E = vars.E - 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_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_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_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_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_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_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 = lcon = fill(-Inf, size(data.storarray))) - - 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)) - - cons = ( - c_ref_angle = c_ref_angle, - c_to_active_power_flow = c_to_active_power_flow, - c_to_reactive_power_flow = c_to_reactive_power_flow, - c_from_active_power_flow = c_from_active_power_flow, - 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_from_thermal_limit = c_from_thermal_limit, - c_to_thermal_limit = c_to_thermal_limit, - c_ramp_rate = c_ramp_rate, - c_active_storage_power = c_active_storage_power, - c_reactive_storage_power = c_reactive_storage_power, - c_ohms = c_ohms, - c_storage_state = c_storage_state, - c_storage_state_init = c_storage_state_init, - c_storage_transfer_thermal_limit = c_storage_transfer_thermal_limit, - c_discharge_thermal_limit = c_discharge_thermal_limit - ) + 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_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_positivity = constraint(core, pstd[s.c, s.t] for s in data.storarray; ucon = fill(Inf, size(data.storarray))) - vars = va, vm, pg, qg, p, q, pst, qst, I2, qint, E, pstd + #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) + cons = (;cons..., c_complementarity = c_complementarity) end - model = ExaModel(core; kwargs...) - return model, vars, 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, + c_discharge_thermal_limit = c_discharge_thermal_limit) + + return vars, cons end -function build_rect_mpopf(data, Nbus, N, discharge_func::Function; backend = nothing, T = Float64, kwargs...) +function add_smooth_cons(core, data, N, vars, cons, discharge_func) - core = ExaCore(T; backend = backend) - vars, cons = build_base_rect_mpopf(core, data, N, Nbus) + pst = vars.pst + pstd = vars.pstd + I2 = vars.I2 + E = vars.E - if length(data.storarray) > 0 - vr, vim, pg, qg, p, q, pst, qst, I2, qint, E = vars - - (c_ref_angle, - c_to_active_power_flow, - c_to_reactive_power_flow, - c_from_active_power_flow, - c_from_reactive_power_flow, - c_phase_angle_diff, - c_active_power_balance, - c_reactive_power_balance, - c_from_thermal_limit, - c_to_thermal_limit, - c_voltage_magnitude, - c_ramp_rate) = cons - - - #discharge or charge from battery to grid (can be positive or negative) - pstd = variable(core, size(data.storage, 1), N; lvar = zeros(size(data.storarray)), uvar = repeat(data.pdmax, 1, N)) - - 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.c, s.t], pstd[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.c, s.t], qint[s.c, s.t], I2[s.c, s.t]) for s 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_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_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 = 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_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_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 = lcon = fill(-Inf, size(data.storarray))) - - 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)) - - cons = ( - c_ref_angle = c_ref_angle, - c_to_active_power_flow = c_to_active_power_flow, - c_to_reactive_power_flow = c_to_reactive_power_flow, - c_from_active_power_flow = c_from_active_power_flow, - 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_from_thermal_limit = c_from_thermal_limit, - c_to_thermal_limit = c_to_thermal_limit, - c_voltage_magnitude = c_voltage_magnitude, - c_ramp_rate = c_ramp_rate, - c_active_storage_power = c_active_storage_power, - c_reactive_storage_power = c_reactive_storage_power, - c_ohms = c_ohms, - c_storage_state = c_storage_state, - c_storage_state_init = c_storage_state_init, - c_storage_transfer_thermal_limit = c_storage_transfer_thermal_limit, - c_discharge_thermal_limit = c_discharge_thermal_limit - ) + 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)) - vars = vr, vim, pg, qg, p, q, pst, qst, I2, qint, E, pstd - end + cons = (;cons..., + c_active_storage_power = c_active_storage_power, + c_storage_state = c_storage_state, + c_storage_state_init = c_storage_state_init, + c_discharge_thermal_limit = c_discharge_thermal_limit) - model = ExaModel(core; kwargs...) - return model, vars, cons + return vars, cons end """ @@ -636,13 +404,11 @@ function mpopf_model( data = convert_data(data,backend) Nbus = size(data.bus, 1) - if form == :polar - return build_polar_mpopf(data, Nbus, N, backend = backend, T = T, storage_complementarity_constraint = storage_complementarity_constraint, kwargs...) - elseif form == :rect - return build_rect_mpopf(data, Nbus, N, backend = backend, T = T, storage_complementarity_constraint = storage_complementarity_constraint, kwargs...) - else + if form != :polar && form != :rect error("Invalid coordinate symbol - valid options are :polar or :rect") end + return build_mpopf(data, Nbus, N, form, backend = backend, T = T, storage_complementarity_constraint = storage_complementarity_constraint, kwargs...) + end function mpopf_model( @@ -664,13 +430,11 @@ function mpopf_model( Nbus = size(data.bus, 1) @assert Nbus == size(pd, 1) - if form == :polar - return build_polar_mpopf(data, Nbus, N, backend = backend, T = T, storage_complementarity_constraint = storage_complementarity_constraint, kwargs...) - elseif form == :rect - return build_rect_mpopf(data, Nbus, N, backend = backend, T = T, storage_complementarity_constraint = storage_complementarity_constraint, kwargs...) - else + if form != :polar && form != :rect error("Invalid coordinate symbol - valid options are :polar or :rect") end + return build_mpopf(data, Nbus, N, form, backend = backend, T = T, storage_complementarity_constraint = storage_complementarity_constraint, kwargs...) + end #Input to discharge_func should be discharge rate (or negative charge), output should be loss in battery level @@ -690,13 +454,11 @@ function mpopf_model( data = convert_data(data,backend) Nbus = size(data.bus, 1) - if form == :polar - return build_polar_mpopf(data, Nbus, N, discharge_func, backend = backend, T = T, kwargs...) - elseif form == :rect - return build_rect_mpopf(data, Nbus, N, discharge_func, backend = backend, T = T, kwargs...) - else + if form != :polar && form != :rect error("Invalid coordinate symbol - valid options are :polar or :rect") end + return build_mpopf(data, Nbus, N, discharge_func, form, backend = backend, T = T, kwargs...) + end function mpopf_model( @@ -719,11 +481,8 @@ function mpopf_model( Nbus = size(data.bus, 1) @assert Nbus == size(pd, 1) - if form == :polar - return build_polar_mpopf(data, Nbus, N, discharge_func, backend = backend, T = T, kwargs...) - elseif form == :rect - return build_rect_mpopf(data, Nbus, N, discharge_func, backend = backend, T = T, kwargs...) - else + if form != :polar && form != :rect error("Invalid coordinate symbol - valid options are :polar or :rect") end + return build_mpopf(data, Nbus, N, discharge_func, form, backend = backend, T = T, kwargs...) end diff --git a/src/parser.jl b/src/parser.jl index 03d2c84..28586fa 100644 --- a/src/parser.jl +++ b/src/parser.jl @@ -134,6 +134,7 @@ function process_ac_power_data(filename) 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) diff --git a/src/sc_parser.jl b/src/sc_parser.jl new file mode 100644 index 0000000..838e2d4 --- /dev/null +++ b/src/sc_parser.jl @@ -0,0 +1,1054 @@ + +function parse_sc_data_static(data) + L_J_xf = length(data.twt_lookup) + L_J_ln = length(data.ac_line_lookup) + L_J_ac = L_J_ln + L_J_xf + L_J_dc = length(data.dc_line_lookup) + L_J_br = L_J_ac + L_J_dc + L_J_cs = length(data.sdd_ids_consumer) + L_J_pr = length(data.sdd_ids_producer) + L_J_cspr = L_J_cs + L_J_pr + L_J_sh = length(data.shunt_lookup) + L_N_p = length(data.azr_lookup) + L_N_q = length(data.rzr_lookup) + I = length(data.bus_lookup) + L_T = length(data.dt) + + lengths = (L_J_xf=L_J_xf, L_J_ln=L_J_ln, L_J_ac=L_J_ac, L_J_dc=L_J_dc, L_J_br=L_J_br, L_J_cs=L_J_cs, + L_J_pr=L_J_pr, L_J_cspr = L_J_cspr, L_J_sh=L_J_sh, I=I, L_T=L_T, L_N_p, L_N_q) + + ε_time = 1e-6 + + cost_vector_pr = sort( + #producers + [ + 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"] + cost = ts_val["cost"] + (j = j, j_prcs = j_prcs, j_pr = j_pr, bus=bus, uid = uid, cost=cost) + end for (key, val) in data.sdd_lookup if val["device_type"] == "producer" + + ], by = x -> x.j) + + #Consumers + cost_vector_cs = 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"] + cost = ts_val["cost"] + (j = j, j_prcs = j_prcs, j_cs = j_cs, bus=bus, uid = uid, cost=cost) + end for (key, val) in data.sdd_lookup if val["device_type"] == "consumer" + ], by = x -> x.j) + + + sc_data = ( + bus = sort([ + begin + i = parse(Int, match(r"\d+", val["uid"]).match)+1 + uid = val["uid"] + v_min = val["vm_lb"] + v_max = val["vm_ub"] + (i = i, uid = uid, v_min = v_min, v_max = v_max) + end for val in values(data.bus_lookup) + ], by = x -> x.i), + + shunt = sort([ + begin + j = parse(Int, match(r"\d+", val["uid"]).match) + L_J_br + L_J_cspr+1 + j_sh = parse(Int, match(r"\d+", val["uid"]).match)+1 + uid = val["uid"] + bus = parse(Int, match(r"\d+", val["bus"]).match)+1 + g_sh = val["gs"] + b_sh = val["bs"] + (j = j, j_sh=j_sh, uid = uid, bus=bus, g_sh = g_sh, b_sh = b_sh) + end for val in values(data.shunt_lookup) + ], by = x -> x.j), + + acl_branch = sort( + # AC lines + [ + begin + j = parse(Int, match(r"\d+", val["uid"]).match)+1 + j_ac = j + j_ln = j + 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"] + 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"] + else + g_fr = 0 + g_to = 0 + b_fr = 0 + b_to = 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) + end for val in values(data.ac_line_lookup) + ],by = x -> x.j), + + # Transformers + acx_branch = sort( + [ + begin + j_xf = parse(Int, match(r"\d+", val["uid"]).match)+1 + j = j_xf + L_J_ln + j_ac = j + 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"] + 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"] + else + g_fr = 0 + g_to = 0 + b_fr = 0 + b_to = 0 + end + (j = j, j_ac = j_ac, j_xf=j_xf, 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) + end for val in values(data.twt_lookup) + ] + , by = x -> x.j), + #Variable phase difference + vpd = isempty(val for val in values(data.twt_lookup) if val["ta_lb"] < val["ta_ub"]) ? empty_data = Vector{NamedTuple{(:j,), Tuple{Int64}}}() : sort([ + begin + j_xf = parse(Int, match(r"\d+", val["uid"]).match)+1 + j = parse(Int, match(r"\d+", val["uid"]).match) + L_J_ln+1 + j_ac = j + phi_min = val["ta_lb"] + phi_max = val["ta_ub"] + (j = j, j_ac = j_ac, j_xf=j_xf, phi_min = phi_min, phi_max = phi_max) + end for val in values(data.twt_lookup) if val["ta_lb"] < val["ta_ub"] + ], by = x -> x.j), + #Fixed phase difference + fpd = isempty(val for val in values(data.twt_lookup) if val["ta_lb"] >= val["ta_ub"]) ? empty_data = Vector{NamedTuple{(:j,), Tuple{Int64}}}() : sort([ + begin + j = parse(Int, match(r"\d+", val["uid"]).match) + L_J_ln+1 + j_ac = j + j_xf = parse(Int, match(r"\d+", val["uid"]).match)+1 + phi_o = val["initial_status"]["ta"] + (j = j, j_ac = j_ac, j_xf=j_xf, phi_o = phi_o) + end for val in values(data.twt_lookup) if val["ta_lb"] >= val["ta_ub"] + ], by = x -> x.j), + #Variable winding ratio + vwr = isempty(val for val in values(data.twt_lookup) if val["tm_lb"] < val["tm_ub"]) ? empty_data = Vector{NamedTuple{(:j,), Tuple{Int64}}}() : sort([ + begin + j = parse(Int, match(r"\d+", val["uid"]).match) + L_J_ln+1 + j_ac = j + j_xf = parse(Int, match(r"\d+", val["uid"]).match)+1 + tau_min = val["tm_lb"] + tau_max = val["tm_ub"] + (j=j, j_ac=j_ac, j_xf=j_xf, tau_min=tau_min, tau_max=tau_max) + end for val in values(data.twt_lookup) if val["tm_lb"] < val["tm_ub"] + ], by = x -> x.j), + #Fixed winding ratio + fwr = isempty(val for val in values(data.twt_lookup) if val["tm_lb"] >= val["tm_ub"]) ? empty_data = Vector{NamedTuple{(:j,), Tuple{Int64}}}() : sort([ + begin + j = parse(Int, match(r"\d+", val["uid"]).match) + L_J_ln+1 + j_ac = j + j_xf = parse(Int, match(r"\d+", val["uid"]).match)+1 + tau_o = val["initial_status"]["tm"] + (j=j, j_ac=j_ac, j_xf=j_xf, tau_o=tau_o) + end for val in values(data.twt_lookup) if val["tm_lb"] >= val["tm_ub"] + + ], by = x -> x.j), + + dc_branch = sort([ + begin + j = parse(Int, match(r"\d+", val["uid"]).match) + L_J_ac+1 + j_dc = parse(Int, match(r"\d+", val["uid"]).match)+1 + uid = val["uid"] + pdc_max = val["pdc_ub"] + qdc_fr_min = val["qdc_fr_lb"] + qdc_to_min = val["qdc_to_lb"] + qdc_fr_max = val["qdc_fr_ub"] + qdc_to_max = val["qdc_to_ub"] + to_bus = parse(Int, match(r"\d+", val["to_bus"]).match)+1 + fr_bus = parse(Int, match(r"\d+", val["fr_bus"]).match)+1 + (j=j, j_dc = j_dc, uid=uid, pdc_max=pdc_max, qdc_fr_min=qdc_fr_min, qdc_to_min=qdc_to_min, qdc_fr_max=qdc_fr_max, qdc_to_max=qdc_to_max, to_bus=to_bus, fr_bus=fr_bus) + end for val in values(data.dc_line_lookup) + + ], by = x -> x.j), + + prod = sort( + # Producers + [ + 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"] + 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"]) + c_nsc = convert(Vector{Float64}, ts_val["p_nsyn_res_cost"]) + c_rru_on = convert(Vector{Float64}, ts_val["p_ramp_res_up_online_cost"]) + c_rru_off = convert(Vector{Float64}, ts_val["p_ramp_res_up_offline_cost"]) + c_rrd_on = convert(Vector{Float64}, ts_val["p_ramp_res_down_online_cost"]) + 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_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"] + + (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, + 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"] + 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"]) + c_nsc = convert(Vector{Float64}, ts_val["p_nsyn_res_cost"]) + c_rru_on = convert(Vector{Float64}, ts_val["p_ramp_res_up_online_cost"]) + c_rru_off = convert(Vector{Float64}, ts_val["p_ramp_res_up_offline_cost"]) + c_rrd_on = convert(Vector{Float64}, ts_val["p_ramp_res_down_online_cost"]) + 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_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"] + + (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, + 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"] == "consumer" + ] + , by = x -> x.j), + active_reserve = sort([ + begin + ts_val = data.azr_ts_lookup[key] + n = parse(Int, match(r"\d+", val["uid"]).match) + 1 + n_p = n + uid = val["uid"] + c_rgu = val["REG_UP_vio_cost"] + c_rgd = val["REG_DOWN_vio_cost"] + c_scr = val["SYN_vio_cost"] + c_nsc = val["NSYN_vio_cost"] + c_rru = val["RAMPING_RESERVE_UP_vio_cost"] + c_rrd = val["RAMPING_RESERVE_DOWN_vio_cost"] + σ_rgu = val["REG_UP"] + σ_rgd = val["REG_DOWN"] + σ_scr = val["SYN"] + σ_nsc = val["NSYN"] + p_rru_min = convert(Vector{Float64}, ts_val["RAMPING_RESERVE_UP"]) + p_rrd_min = convert(Vector{Float64}, ts_val["RAMPING_RESERVE_DOWN"]) + (n=n, n_p=n_p, uid=uid, c_rgu=c_rgu, c_rgd=c_rgd, c_scr=c_scr, c_nsc=c_nsc, c_rru=c_rru, c_rrd=c_rrd, σ_rgu=σ_rgu, σ_rgd=σ_rgd, σ_scr=σ_scr, + σ_nsc=σ_nsc, p_rru_min=p_rru_min, p_rrd_min=p_rrd_min) + end for (key, val) in data.azr_lookup + ], by = x -> x.n), + reactive_reserve = sort([ + begin + ts_val = data.rzr_ts_lookup[key] + n = parse(Int, match(r"\d+", val["uid"]).match) + L_N_p + 1 + n_q = parse(Int, match(r"\d+", val["uid"]).match) + 1 + uid = val["uid"] + c_qru = val["REACT_UP_vio_cost"] + c_qrd = val["REACT_DOWN_vio_cost"] + q_qru_min = convert(Vector{Float64}, ts_val["REACT_UP"]) + q_qrd_min = convert(Vector{Float64}, ts_val["REACT_DOWN"]) + (n=n, n_q=n_q, uid=uid, c_qru=c_qru, c_qrd=c_qrd, q_qru_min=q_qru_min, q_qrd_min=q_qrd_min) + end for (key, val) in data.rzr_lookup + ], 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) + 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 + ], + + 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) + 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 + ], + + 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) + 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 + ], + + 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) + 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 + ], + + ) + return sc_data, lengths, cost_vector_pr, cost_vector_cs +end + +function add_status_flags(uc_data, data) + for dict in uc_data + uid = dict["uid"] + on_status = dict["on_status"] + n = length(on_status) + u_on_o = data[uid]["initial_status"]["on_status"] + + su_status = zeros(Int, n) + sd_status = zeros(Int, n) + + # Handle first index using u_on_o if provided + + if u_on_o == 0 && on_status[1] == 1 + su_status[1] = 1 + elseif u_on_o == 1 && on_status[1] == 0 + sd_status[1] = 1 + end + + + # Iterate from 1 to n-1 for the rest + for i in 1:n-1 + if on_status[i] == 0 && on_status[i+1] == 1 + su_status[i+1] = 1 + end + if on_status[i] == 1 && on_status[i+1] == 0 + sd_status[i+1] = 1 + end + end + + dict["su_status"] = su_status + dict["sd_status"] = sd_status + end +end + +function get_as(dt, t) + a_end = sum(dt[1:t]) + a_start = a_end - dt[t] + a_mid = (a_start + a_end)/2 + return a_start, a_mid, a_end +end + +function parse_sc_data(data, uc_data, data_json) + sc_data, lengths, cost_vector_pr, cost_vector_cs = parse_sc_data_static(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 + periods = data.periods + dt = Float64.(data.dt) + K = length(data_json["reliability"]["contingency"]) + add_status_flags(uc_data["time_series_output"]["ac_line"], data.ac_line_lookup) + add_status_flags(uc_data["time_series_output"]["two_winding_transformer"], data.twt_lookup) + 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, + 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 + 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 + for uc in uc_data["time_series_output"]["simple_dispatchable_device"] + if val["uid"] == uc["uid"] + ] + + #This sum corresponds to constraint 69 (summing p_supc*u_su) + sum_T_supc_pr = zeros(L_J_pr, length(periods)) + #This sum corresponds to constraint 112/113 (summing u_su) + sum2_T_supc_pr = zeros(L_J_pr, length(periods)) + + for b in T_supc_pr + sum_T_supc_pr[b.j_pr, b.t] += b.p_supc*b.u_su + sum2_T_supc_pr[b.j_pr, b.t] += b.u_su + end + + 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, + 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 + 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 + for uc in uc_data["time_series_output"]["simple_dispatchable_device"] + if val["uid"] == uc["uid"] + ] + #This sum corresponds to constraint 69 (p_supc*u_su) + sum_T_supc_cs = zeros(L_J_cs, L_T) + #This sum corresponds to constraints 122-126 (u_su) + sum2_T_supc_cs = zeros(L_J_cs, L_T) + for b in T_supc_cs + sum_T_supc_cs[b.j_cs, b.t] += b.p_supc*b.u_su + sum2_T_supc_cs[b.j_cs, b.t] += b.u_su + end + + + #Build a p_sdpc set to be used for T_sdpc + #Index order is j_prcs, t, t_prime + p_sdpc = zeros(L_J_cspr, L_T, L_T) + + for t in periods + 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]) + 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]) + end + end + end + end + + 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, + 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 + 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 + for uc in uc_data["time_series_output"]["simple_dispatchable_device"] + if val["uid"] == uc["uid"] + ] + + #This sum corresponds to constraint 70 (summing p_sdpc*u_sd) + sum_T_sdpc_pr = zeros(L_J_pr, L_T) + #This sum corresponds to constraints 112 and 113 (summing u_sd) + sum2_T_sdpc_pr = zeros(L_J_pr, L_T) + for b in T_sdpc_pr + sum_T_sdpc_pr[b.j_pr, b.t] += b.p_sdpc*b.u_sd + sum2_T_sdpc_pr[b.j_pr, b.t] += b.u_sd + end + + 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, + 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 + 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 + for uc in uc_data["time_series_output"]["simple_dispatchable_device"] + if val["uid"] == uc["uid"] + ] + #This sum corresponds to constraint 70 (p_sdpc*u_sd) + sum_T_sdpc_cs = zeros(L_J_cs, L_T) + #This sum corresponds to constraint 122-126 (u_sd) + sum2_T_sdpc_cs = zeros(L_J_cs, L_T) + for b in T_sdpc_cs + sum_T_sdpc_cs[b.j_cs, b.t] += b.p_sdpc*b.u_sd + sum2_T_sdpc_cs[b.j_cs, b.t] += b.u_sd + end + + + 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 + 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])) + w_en_max_pr_ind +=1 + end + end + end + + L_W_en_max_pr = length(W_en_max_pr) + + 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 + 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])) + w_en_max_cs_ind +=1 + end + end + end + + L_W_en_max_cs = length(W_en_max_cs) + + 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 + 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, + a_en_min_start = w[1], + a_en_min_end = w[2], + e_min = w[3])) + w_en_min_pr_ind += 1 + end + end + end + + L_W_en_min_pr = length(W_en_min_pr) + + 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 + 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, + a_en_min_start = w[1], + a_en_min_end = w[2], + e_min = w[3])) + w_en_min_cs_ind += 1 + end + end + end + + L_W_en_min_cs = length(W_en_min_cs) + + 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 + 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, + t = t, + dt=dt[t])) + end + end + end + end + end + + 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 + 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, + t = t, + dt=dt[t])) + end + end + end + end + end + + 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 + 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, + t = t, + dt = dt[t])) + end + end + end + end + end + + 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 + 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, + t = t, + dt = dt[t])) + end + end + end + end + end + + p_jtm_flattened_pr = Vector{@NamedTuple{flat_k::Int, j::Int, j_prcs::Int, j_pr::Int, t::Int, m::Int, c_en::Float64, p_max::Float64}}() + flat_k=1 + for (pc_idx, pc) in pairs(cost_vector_pr) + j=pc.j + j_prcs=pc.j_prcs + j_pr = pc.j_pr + 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)) + flat_k+=1 + end + end + end + + p_jtm_flattened_cs = Vector{@NamedTuple{flat_k::Int, j::Int, j_prcs::Int, j_cs::Int, t::Int, m::Int, c_en::Float64, p_max::Float64}}() + flat_k=1 + for (pc_idx, pc) in pairs(cost_vector_cs) + j=pc.j + j_prcs=pc.j_prcs + j_cs = pc.j_cs + 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)) + flat_k+=1 + end + end + end + + jtk_ln_flattened = Vector{@NamedTuple{flat_jtk_ln::Int, ctg::Int, j::Int, j_ac::Int, j_ln::Int, to_bus::Int, fr_bus::Int, b_sr::Float64, s_max_ctg::Float64, u_on::Int, t::Int, dt::Float64}}() + flat_jtk_ln = 1 + for ctg in data_json["reliability"]["contingency"], t in periods + for component in ctg["components"] + for val in values(data.ac_line_lookup) + if val["uid"] != component + for uc in uc_data["time_series_output"]["ac_line"] + if val["uid"] == uc["uid"] + r = val["r"] + x = val["x"] + push!(jtk_ln_flattened, (flat_jtk_ln=flat_jtk_ln, ctg = parse(Int, match(r"\d+", ctg["uid"]).match)+1, j = parse(Int, match(r"\d+", val["uid"]).match)+1, j_ac = parse(Int, match(r"\d+", val["uid"]).match)+1, + j_ln = parse(Int, match(r"\d+", val["uid"]).match)+1, to_bus = parse(Int, match(r"\d+", val["to_bus"]).match)+1, + fr_bus = parse(Int, match(r"\d+", val["fr_bus"]).match)+1, b_sr = -x / (x^2 + r^2), s_max_ctg = val["mva_ub_em"], u_on=uc["on_status"][t], t=t, dt = dt[t])) + flat_jtk_ln += 1 + end + end + end + end + end + end + + + + jtk_xf_flattened = Vector{@NamedTuple{flat_jtk_xf::Int, ctg::Int, j::Int, j_ac::Int, j_xf::Int, to_bus::Int, fr_bus::Int, b_sr::Float64, s_max_ctg::Float64, u_on::Int, t::Int, dt::Float64}}() + flat_jtk_xf = 1 + for ctg in data_json["reliability"]["contingency"], t in periods + for component in ctg["components"] + for val in values(data.twt_lookup) + if val["uid"] != component + for uc in uc_data["time_series_output"]["two_winding_transformer"] + if val["uid"] == uc["uid"] + r = val["r"] + x = val["x"] + push!(jtk_xf_flattened, (flat_jtk_xf=flat_jtk_xf, ctg = parse(Int, match(r"\d+", ctg["uid"]).match)+1, j = parse(Int, match(r"\d+", val["uid"]).match)+1+ L_J_ln, j_ac = parse(Int, match(r"\d+", val["uid"]).match)+1+ L_J_ln, + j_xf = parse(Int, match(r"\d+", val["uid"]).match)+1, to_bus = parse(Int, match(r"\d+", val["to_bus"]).match)+1, + fr_bus = parse(Int, match(r"\d+", val["fr_bus"]).match)+1, b_sr = -x / (x^2 + r^2), s_max_ctg = val["mva_ub_em"], u_on=uc["on_status"][t], t=t, dt = dt[t])) + flat_jtk_xf += 1 + end + end + end + end + end + end + + jtk_dc_flattened = Vector{@NamedTuple{flat_jtk_dc::Int, ctg::Int, j::Int, j_dc::Int, to_bus::Int, fr_bus::Int, t::Int, dt::Float64}}() + flat_jtk_dc = 1 + for ctg in data_json["reliability"]["contingency"], t in periods + for component in ctg["components"] + for val in values(data.dc_line_lookup) + if val["uid"] != component + push!(jtk_dc_flattened, (flat_jtk_dc=flat_jtk_dc, ctg = parse(Int, match(r"\d+", ctg["uid"]).match)+1, j = parse(Int, match(r"\d+", val["uid"]).match)+1+ L_J_ac, + j_dc = parse(Int, match(r"\d+", val["uid"]).match)+1, to_bus = parse(Int, match(r"\d+", val["to_bus"]).match)+1, + fr_bus = parse(Int, match(r"\d+", val["fr_bus"]).match)+1, t=t, dt = dt[t])) + flat_jtk_dc += 1 + end + end + end + end + + empty_vpd = Vector{NamedTuple{(:j, :j_ac, :j_xf, :phi_min, :phi_max, :t), Tuple{Int64, Int64, Int64, Float64, Float64, Int64}}}() + empty_fpd = Vector{NamedTuple{(:j, :j_ac, :j_xf, :phi_o, :t), Tuple{Int64, Int64, Int64, Float64, Int64}}}() + empty_vwr = Vector{NamedTuple{(:j, :j_ac, :j_xf, :tau_min, :tau_max, :t), Tuple{Int64, Int64, Int64, Float64, Float64, Int64}}}() + empty_fwr = Vector{NamedTuple{(:j, :j_ac, :j_xf, :tau_o, :t), Tuple{Int64, Int64, Int64, Float64, Int64}}}() + + sc_time_data = ( + ; + periods = periods, + + c_p = [data.violation_cost["p_bus_vio_cost"]], + c_q = [data.violation_cost["q_bus_vio_cost"]], + c_s = [data.violation_cost["s_vio_cost"]], + c_e = [data.violation_cost["e_vio_cost"]], + busarray = [(;b..., t=t, dt=dt[t]) for b in sc_data.bus, t in periods], + k_busarray = [(;b..., t=t, k=k) for b in sc_data.bus, t in periods, k in 1:length(data_json["reliability"]["contingency"])], + shuntarray = [ + (;s..., t=t, u_sh = uc["step"][t]) + for s in sc_data.shunt, t in periods + for uc in uc_data["time_series_output"]["shunt"] + if s.uid == uc["uid"] + ], + k_shuntarray = [(;b..., t=t, k=k) for b in sc_data.shunt, t in periods, k in 1:length(data_json["reliability"]["contingency"])], + + + preservearray = [(;n=r.n, n_p=r.n_p, uid=r.uid, c_rgu=r.c_rgu, c_rgd=r.c_rgd, c_scr=r.c_scr, c_nsc=r.c_nsc, c_rru=r.c_rru, c_rrd=r.c_rrd, + σ_rgu=r.σ_rgu, σ_rgd=r.σ_rgd, σ_scr=r.σ_scr, σ_nsc=r.σ_nsc, p_rru_min=r.p_rru_min[t], p_rrd_min=r.p_rrd_min[t], + t=t, dt=dt[t]) for r in sc_data.active_reserve, t in periods], + + qreservearray = [(;n=q.n, n_q=q.n_q, uid=q.uid, c_qru=q.c_qru, c_qrd=q.c_qrd, q_qru_min=q.q_qru_min[t], q_qrd_min=q.q_qrd_min[t], t=t, dt = dt[t]) + for q in sc_data.reactive_reserve, t in periods], + + k_prarray = [(;j_prcs=b.j_prcs, j_pr=b.j_pr, bus = b.bus, uid=b.uid, t=t, k=k) for b in sc_data.prod, t in periods, k in 1:length(data_json["reliability"]["contingency"])], + prarray = [(;j=p.j, j_prcs=p.j_prcs, j_pr=p.j_pr, bus=p.bus, uid=p.uid, c_on=p.c_on, c_sd=p.c_sd, c_su = p.c_su, p_ru=p.p_ru, p_rd=p.p_rd, + p_ru_su=p.p_ru_su, p_rd_sd=p.p_rd_sd, c_rgu=p.c_rgu[t], c_rgd=p.c_rgd[t], c_scr=p.c_scr[t], c_nsc=p.c_nsc[t], c_rru_on=p.c_rru_on[t], + c_rru_off=p.c_rru_off[t], c_rrd_on=p.c_rrd_on[t], c_rrd_off=p.c_rrd_off[t], c_qru=p.c_qru[t], c_qrd=p.c_qrd[t], p_rgu_max=p.p_rgu_max, + p_rgd_max=p.p_rgd_max, p_scr_max=p.p_scr_max, p_nsc_max=p.p_nsc_max, p_rru_on_max=p.p_rru_on_max, p_rru_off_max=p.p_rru_off_max, + p_rrd_on_max=p.p_rrd_on_max, p_rrd_off_max=p.p_rrd_off_max, p_0=p.p_0, q_0=p.q_0, p_max=p.p_max[t], p_min=p.p_min[t], q_max=p.q_max[t], q_min=p.q_min[t], #sus = p.sus, + u_on = uc["on_status"][t], u_su = uc["su_status"][t], u_sd = uc["sd_status"][t], t=t, + sum_T_supc_pr_jt = sum_T_supc_pr[p.j_pr, t], sum_T_sdpc_pr_jt = sum_T_sdpc_pr[p.j_pr, 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], dt = dt[t]) + for p in sc_data.prod, t in periods + 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) ? + 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, + 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], + beta_max = val["beta_ub"], + beta_min = val["beta_lb"], + q_max_p0 = val["q_0_ub"], + q_min_p0 = val["q_0_lb"], + 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 + 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) ? + 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, + 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], + beta = val["beta"], + q_p0 = val["q_0"], + 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 + for uc in uc_data["time_series_output"]["simple_dispatchable_device"] + if p.uid == uc["uid"]], + + k_csarray = [(;j_prcs=b.j_prcs, j_cs=b.j_cs, bus = b.bus, uid=b.uid, t=t, k=k) for b in sc_data.cons, t in periods, k in 1:length(data_json["reliability"]["contingency"])], + + csarray = [(;j=p.j, j_prcs=p.j_prcs, j_cs=p.j_cs, bus=p.bus, uid=p.uid, c_on=p.c_on, c_sd=p.c_sd, c_su=p.c_su, p_ru=p.p_ru, p_rd=p.p_rd, + p_ru_su=p.p_ru_su, p_rd_sd=p.p_rd_sd, c_rgu=p.c_rgu[t], c_rgd=p.c_rgd[t], c_scr=p.c_scr[t], c_nsc=p.c_nsc[t], c_rru_on=p.c_rru_on[t], + c_rru_off=p.c_rru_off[t], c_rrd_on=p.c_rrd_on[t], c_rrd_off=p.c_rrd_off[t], c_qru=p.c_qru[t], c_qrd=p.c_qrd[t], p_rgu_max=p.p_rgu_max, + p_rgd_max=p.p_rgd_max, p_scr_max=p.p_scr_max, p_nsc_max=p.p_nsc_max, p_rru_on_max=p.p_rru_on_max, p_rru_off_max=p.p_rru_off_max, + p_rrd_on_max=p.p_rrd_on_max, p_rrd_off_max=p.p_rrd_off_max, p_0=p.p_0, q_0=p.q_0, p_max=p.p_max[t], p_min=p.p_min[t], q_max=p.q_max[t], q_min=p.q_min[t], #sus=p.sus, + u_on = uc["on_status"][t], u_su = uc["su_status"][t], u_sd = uc["sd_status"][t], t=t, + sum_T_supc_cs_jt = sum_T_supc_cs[p.j_cs, t], sum_T_sdpc_cs_jt = sum_T_sdpc_cs[p.j_cs, 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], dt = dt[t]) + for p in sc_data.cons, t in periods + 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) ? + 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, + 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], + beta_max = val["beta_ub"], + beta_min = val["beta_lb"], + q_max_p0 = val["q_0_ub"], + q_min_p0 = val["q_0_lb"], + 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 + 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) ? + 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, + 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], + beta = val["beta"], + q_p0 = val["q_0"], + 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 + for uc in uc_data["time_series_output"]["simple_dispatchable_device"] + if p.uid == uc["uid"]], + + acxbrancharray = [ + (;j=b.j, j_ac=b.j_ac, j_xf=b.j_xf, uid=b.uid, to_bus=b.to_bus, fr_bus=b.fr_bus, c_su=b.c_su, c_sd=b.c_sd, s_max=b.s_max, g_sr=b.g_sr, b_sr=b.b_sr, b_ch=b.b_ch, + g_fr=b.g_fr, g_to=b.g_to, b_fr=b.b_fr, b_to=b.b_to, u_on=uc["on_status"][t], u_su=uc["su_status"][t], u_sd=uc["sd_status"][t], t=t, dt = dt[t]) + for b in sc_data.acx_branch, t in periods + for uc in uc_data["time_series_output"]["two_winding_transformer"] + if b.uid == uc["uid"]], + + aclbrancharray = [ + (;j=b.j, j_ac=b.j_ac, j_ln=b.j_ln, uid=b.uid, to_bus=b.to_bus, fr_bus=b.fr_bus, c_su=b.c_su, c_sd=b.c_sd, s_max=b.s_max, g_sr=b.g_sr, b_sr=b.b_sr, b_ch=b.b_ch, + g_fr=b.g_fr, g_to=b.g_to, b_fr=b.b_fr, b_to=b.b_to, u_on=uc["on_status"][t], u_su=uc["su_status"][t], u_sd=uc["sd_status"][t], t=t, dt = dt[t]) + for b in sc_data.acl_branch, t in periods + for uc in uc_data["time_series_output"]["ac_line"] + if b.uid == uc["uid"]], + + fpdarray = isempty(sc_data.fpd) ? empty_data = empty_fpd : [(;b..., t=t) for b in sc_data.fpd, t in periods], + fwrarray = isempty(sc_data.fwr) ? empty_data = empty_fwr : [(;b..., t=t) for b in sc_data.fwr, t in periods], + vpdarray = isempty(sc_data.vpd) ? empty_data = empty_vpd : [(;b..., t=t) for b in sc_data.vpd, t in periods], + vwrarray = isempty(sc_data.vwr) ? empty_data = empty_vwr : [(;b..., t=t) for b in sc_data.vwr, t in periods], + dclinearray = [(;b..., t=t) for b in sc_data.dc_branch, t in periods], + + p_jt_fr_dc_max = [dc.pdc_max for dc in sc_data.dc_branch, t in periods], + p_jt_to_dc_max = [dc.pdc_max for dc in sc_data.dc_branch, t in periods], + q_jt_fr_dc_lvar = [dc.qdc_fr_min for dc in sc_data.dc_branch, t in periods], + q_jt_fr_dc_uvar = [dc.qdc_fr_max for dc in sc_data.dc_branch, t in periods], + q_jt_to_dc_lvar = [dc.qdc_to_min for dc in sc_data.dc_branch, t in periods], + q_jt_to_dc_uvar = [dc.qdc_to_max for dc in sc_data.dc_branch, t in periods], + + preservesetarray_pr = [(;b..., t=t) for b in sc_data.active_reserve_set_pr, t in periods], + preservesetarray_cs = [(;b..., t=t) for b in sc_data.active_reserve_set_cs, t in periods], + qreservesetarray_pr = [(;b..., t=t) for b in sc_data.reactive_reserve_set_pr, t in periods], + qreservesetarray_cs = [(;b..., t=t) for b in sc_data.reactive_reserve_set_cs, t in periods], + + v_lvar = repeat([b.v_min for b in sc_data.bus], 1, L_T), + v_uvar = repeat([b.v_max for b in sc_data.bus], 1, L_T), + + W_en_max_pr=W_en_max_pr, + W_en_max_cs=W_en_max_cs, + T_w_en_max_pr=T_w_en_max_pr, + T_w_en_max_cs=T_w_en_max_cs, + W_en_min_pr=W_en_min_pr, + W_en_min_cs=W_en_min_cs, + T_w_en_min_pr=T_w_en_min_pr, + T_w_en_min_cs=T_w_en_min_cs, + + tk_index = [(t=t, k=k) for t in periods, k in 1:length(data_json["reliability"]["contingency"])], + + p_jtm_flattened_pr=p_jtm_flattened_pr, + p_jtm_flattened_cs=p_jtm_flattened_cs, + p_jtm_pr_uvar = [b.p_max for b in p_jtm_flattened_pr], + p_jtm_cs_uvar = [b.p_max for b in p_jtm_flattened_cs], + jtk_ln_flattened=jtk_ln_flattened, + jtk_xf_flattened=jtk_xf_flattened, + jtk_dc_flattened=jtk_dc_flattened, + + ) + + 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 +end + +function save_go3_solution(uc_filename, solution_name, result, vars, lengths) + 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 + #This section corresponds to consuming devices + solution_index -= L_J_pr + 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) + line["p_reg_res_up"] = Array(solution(result, vars.p_jt_rgu_cs))[solution_index,:] + line["p_ramp_res_down_online"] = Array(solution(result, vars.p_jt_rrd_on_cs))[solution_index,:] + line["p_on"] = Array(solution(result, vars.p_jt_on_cs))[solution_index,:] + line["q"] = Array(solution(result, vars.q_jt_cs))[solution_index,:] + line["p_reg_res_down"] = Array(solution(result, vars.p_jt_rgd_cs))[solution_index,:] + line["p_ramp_res_up_offline"] = zeros(L_T) + line["q_res_down"] = Array(solution(result, vars.q_jt_qrd_cs))[solution_index,:] + line["q_res_up"] = Array(solution(result, vars.q_jt_qru_cs))[solution_index,:] + line["p_ramp_res_down_offline"] = Array(solution(result, vars.p_jt_rrd_off_cs))[solution_index,:] + else + #Producing devices + 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,:] + line["p_reg_res_up"] = Array(solution(result, vars.p_jt_rgu_pr))[solution_index,:] + line["p_ramp_res_down_online"] = Array(solution(result, vars.p_jt_rrd_on_pr))[solution_index,:] + line["p_on"] = Array(solution(result, vars.p_jt_on_pr))[solution_index,:] + line["q"] = Array(solution(result, vars.q_jt_pr))[solution_index,:] + line["p_reg_res_down"] = Array(solution(result, vars.p_jt_rgd_pr))[solution_index,:] + line["p_ramp_res_up_offline"] = Array(solution(result, vars.p_jt_rru_off_pr))[solution_index,:] + line["q_res_down"] = Array(solution(result, vars.q_jt_qrd_pr))[solution_index,:] + line["q_res_up"] = Array(solution(result, vars.q_jt_qru_pr))[solution_index,:] + line["p_ramp_res_down_offline"] = zeros(L_T) + end + end + #Update two winding transformers + for line in uc_data["time_series_output"]["two_winding_transformer"] + solution_index = parse(Int, match(r"\d+", line["uid"]).match) + 1 #This corresponds to j_xf + line["tm"] = Array(solution(result, vars.τ_jt_xf))[solution_index,:] + line["ta"] = Array(solution(result, vars.φ_jt_xf))[solution_index,:] + end + #ac line flows can be inferred from voltage and angle levels + + #Update DC lines + for line in uc_data["time_series_output"]["dc_line"] + solution_index = parse(Int, match(r"\d+", line["uid"]).match) + 1 #This corresponds to j_dc + line["qdc_fr"] = Array(solution(result, vars.q_jt_fr_dc))[solution_index,:] + line["pdc_fr"] = Array(solution(result, vars.p_jt_fr_dc))[solution_index,:] #pdc_to is implied from energy conservation + line["qdc_to"] = Array(solution(result, vars.q_jt_to_dc))[solution_index,:] + end + + #Update buses + for line in uc_data["time_series_output"]["bus"] + solution_index = parse(Int, match(r"\d+", line["uid"]).match) + 1 #This corresponds to i + line["vm"] = Array(solution(result, vars.v_it))[solution_index,:] + line["va"] = Array(solution(result, vars.θ_it))[solution_index,:] + end + + open(solution_name, "w") do io + JSON.print(io, uc_data, 4) + end +end diff --git a/src/scopf.jl b/src/scopf.jl index 717d9c4..06db084 100644 --- a/src/scopf.jl +++ b/src/scopf.jl @@ -1,214 +1,680 @@ -function parse_sc_power_data(data, contingencies, corrective_action_ratio, backend) - # Raw data - nbus = length(data.bus) - ngen = length(data.gen) - narc = length(data.arc) - # Contingencies - K = length(contingencies) + 1 # +1 accounts for base case - - # Build indexes - idx_ref = [(r, k) for r in data.ref_buses, k in 1:K] - idx_branch = [(b, 1) for b in data.branch] - idx_branch_down = similar(idx_branch, 0) - # Scan contingencies - for k in 2:K # k=1 is base case - line_down = contingencies[k-1] # ID of line down - for b in data.branch - if b.i != line_down - push!(idx_branch, (b, k)) - else - push!(idx_branch_down, (b, k)) - end - end - end - - idx_bus = [(b, k) for b in data.bus, k in 1:K] - idx_gen = [(b, k) for b in data.gen, k in 1:K] - idx_arc = [(b, k) for b in data.arc, k in 1:K] - - # Compute maximum injection in lines - max_inj = repeat(data.rate_a, K) - # Remove (pinj, qinj) for lines down by setting the bounds to 0 - # so the solver can convert it as a fixed variable. - for (b, k) in idx_branch_down - max_inj[b.f_idx + (k-1) * narc] = 0.0 - max_inj[b.t_idx + (k-1) * narc] = 0.0 - end - gen_buses = [b for b in data.bus if b[:bus_type] in [2, 3]] # Select buses with generators (=PV and REF buses) - - data = convert_data( - ( - data..., - idx_ref = idx_ref, - idx_branch = idx_branch, - idx_branch_down = idx_branch_down, - idx_bus = idx_bus, - idx_gen = idx_gen, - idx_arc = idx_arc, - max_inj = max_inj, - Δp = corrective_action_ratio .* (data.pmax .- data.pmin), - idx_gen_cont = [(b, k) for b in data.gen, k in 2:K], # discard base case - idx_bus_cont = [(b, k) for b in gen_buses, k in 2:K], - ), - backend) - return ( - data..., - K = K, - nbus = nbus, - ngen = ngen, - narc = narc, - ) -end +using GOC3Benchmark, JSON, ExaModels function scopf_model( - filename; + filename, uc_filename; backend = nothing, - data = parse_ac_power_data(filename), - contingencies = [b.i for b in data.branch], # full n-1 formulation by default - corrective_action_ratio = 0.1, T = Float64, + include_ctg = true, #contingencies will be specified in input data + result_set = [], kwargs... ) - data = parse_sc_power_data(data, contingencies, corrective_action_ratio, backend) + + 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) + + (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 + + sc_data_array = sc_data + sc_data = convert_data(sc_data, backend) + + core = ExaCore(T; backend =backend) + + if result_set != [] + initialize_vars = true + result = result_set[1] + vars = result_set[2] + else + initialize_vars = false + end + + #variables are indexed j,t,k or j,t (t always second if present) + b_jt_sh = variable(core, L_J_sh, L_T; start = initialize_vars ? Array(solution(result, vars.b_jt_sh)) : 1) + g_jt_sh = variable(core, L_J_sh, L_T; start = initialize_vars ? Array(solution(result, vars.g_jt_sh)) : 1) + #Split e_w_plus into separate sets for W_en_min and W_en_max and for pr, cs + #Bounds from 4.6.3 Maximum/minimum energy over multiple intervals (77) + e_w_plus_min_pr = variable(core, L_W_en_min_pr; lvar = 0, start = initialize_vars ? Array(solution(result, vars.e_w_plus_min_pr)) : 0) + e_w_plus_min_cs = variable(core, L_W_en_min_cs; lvar = 0, start = initialize_vars ? Array(solution(result, vars.e_w_plus_min_cs)) : 0) + e_w_plus_max_pr = variable(core, L_W_en_max_pr; lvar = 0, start = initialize_vars ? Array(solution(result, vars.e_w_plus_max_pr)) : 0) + e_w_plus_max_cs = variable(core, L_W_en_max_cs; lvar = 0, start = initialize_vars ? Array(solution(result, vars.e_w_plus_max_cs)) : 0) + p_it = variable(core, I, L_T; start = initialize_vars ? Array(solution(result, vars.p_it)) : 1) + p_it_plus = variable(core, I, L_T; start = initialize_vars ? Array(solution(result, vars.p_it_plus)) : 1) + #splitting p_jt and q_jt for shunts, producers, and consumers + p_jt_sh = variable(core, L_J_sh, L_T; start = initialize_vars ? Array(solution(result, vars.p_jt_sh)) : 1) + p_jt_pr = variable(core, L_J_pr, L_T; start = initialize_vars ? Array(solution(result, vars.p_jt_pr)) : 1) + p_jt_cs = variable(core, L_J_cs, L_T; start = initialize_vars ? Array(solution(result, vars.p_jt_cs)) : 1) + q_jt_sh = variable(core, L_J_sh, L_T; start = initialize_vars ? Array(solution(result, vars.q_jt_sh)) : 1) + q_jt_pr = variable(core, L_J_pr, L_T; start = initialize_vars ? Array(solution(result, vars.q_jt_pr)) : 1) + q_jt_cs = variable(core, L_J_cs, L_T; start = initialize_vars ? Array(solution(result, vars.q_jt_cs)) : 1) + #Splitting p on, sd , su into pr and cs + p_jt_on_pr = variable(core, L_J_pr, L_T; start = initialize_vars ? Array(solution(result, vars.p_jt_on_pr)) : 1) + p_jt_on_cs = variable(core, L_J_cs, L_T; start = initialize_vars ? Array(solution(result, vars.p_jt_on_cs)) : 1) + p_jt_su_pr = variable(core, L_J_pr, L_T; start = initialize_vars ? Array(solution(result, vars.p_jt_su_pr)) : 1) + p_jt_su_cs = variable(core, L_J_cs, L_T; start = initialize_vars ? Array(solution(result, vars.p_jt_su_cs)) : 1) + p_jt_sd_pr = variable(core, L_J_pr, L_T; start = initialize_vars ? Array(solution(result, vars.p_jt_sd_pr)) : 1) + p_jt_sd_cs = variable(core, L_J_cs, L_T; start = initialize_vars ? Array(solution(result, vars.p_jt_sd_cs)) : 1) + #p_jtm has been flattened and uses only one special index, k_flat + #Bounds from 4.6.9 Energy cost and value (129) + p_jtm_pr = variable(core, length(sc_data.p_jtm_flattened_pr); lvar = 0, uvar = sc_data.p_jtm_pr_uvar, start = initialize_vars ? Array(solution(result, vars.p_jtm_pr)) : 1) + p_jtm_cs = variable(core, length(sc_data.p_jtm_flattened_cs); lvar = 0, uvar = sc_data.p_jtm_cs_uvar, start = initialize_vars ? Array(solution(result, vars.p_jtm_cs)) : 1) + #to/from power split into ln, xf, and dc lines + #Bounds from 4.8.4 DC lines (152-155) + p_jt_fr_ln = variable(core, L_J_ln, L_T; start = initialize_vars ? Array(solution(result, vars.p_jt_fr_ln)) : 1) + p_jt_fr_xf = variable(core, L_J_xf, L_T; start = initialize_vars ? Array(solution(result, vars.p_jt_fr_xf)) : 1) + p_jt_fr_dc = variable(core, L_J_dc, L_T; lvar = -sc_data.p_jt_fr_dc_max, uvar = sc_data.p_jt_fr_dc_max, start = initialize_vars ? Array(solution(result, vars.p_jt_fr_dc)) : 1) + p_jt_to_ln = variable(core, L_J_ln, L_T; start = initialize_vars ? Array(solution(result, vars.p_jt_to_ln)) : 1) + p_jt_to_xf = variable(core, L_J_xf, L_T; start = initialize_vars ? Array(solution(result, vars.p_jt_to_xf)) : 1) + p_jt_to_dc = variable(core, L_J_dc, L_T; lvar = -sc_data.p_jt_to_dc_max, uvar = sc_data.p_jt_to_dc_max, start = initialize_vars ? Array(solution(result, vars.p_jt_to_dc)) : 1) + q_jt_fr_ln = variable(core, L_J_ln, L_T; start = initialize_vars ? Array(solution(result, vars.q_jt_fr_ln)) : 1) + q_jt_fr_xf = variable(core, L_J_xf, L_T; start = initialize_vars ? Array(solution(result, vars.q_jt_fr_xf)) : 1) + q_jt_fr_dc = variable(core, L_J_dc, L_T; lvar = sc_data.q_jt_fr_dc_lvar, uvar = sc_data.q_jt_fr_dc_uvar, start = initialize_vars ? Array(solution(result, vars.q_jt_fr_dc)) : 1) + q_jt_to_ln = variable(core, L_J_ln, L_T; start = initialize_vars ? Array(solution(result, vars.q_jt_to_ln)) : 1) + q_jt_to_xf = variable(core, L_J_xf, L_T; start = initialize_vars ? Array(solution(result, vars.q_jt_to_xf)) : 1) + q_jt_to_dc = variable(core, L_J_dc, L_T; lvar = sc_data.q_jt_to_dc_lvar, uvar = sc_data.q_jt_to_dc_uvar, start = initialize_vars ? Array(solution(result, vars.q_jt_to_dc)) : 1) + #p_jt rgu, rgd, scr, rru,on, rru,off, rrd,on, rrd,off and q_jt qru/qrd split into pr and cs + #bounds from 4.6.4 Device reserve variable domains (80-89) + p_jt_rgu_pr = variable(core, L_J_pr, L_T; lvar = 0, start = initialize_vars ? Array(solution(result, vars.p_jt_rgu_pr)) : 1) + p_jt_rgu_cs = variable(core, L_J_cs, L_T; lvar = 0, start = initialize_vars ? Array(solution(result, vars.p_jt_rgu_cs)) : 1) + p_jt_rgd_pr = variable(core, L_J_pr, L_T; lvar = 0, start = initialize_vars ? Array(solution(result, vars.p_jt_rgd_pr)) : 1) + p_jt_rgd_cs = variable(core, L_J_cs, L_T; lvar = 0, start = initialize_vars ? Array(solution(result, vars.p_jt_rgd_cs)) : 1) + p_jt_scr_pr = variable(core, L_J_pr, L_T; lvar = 0, start = initialize_vars ? Array(solution(result, vars.p_jt_scr_pr)) : 1) + p_jt_scr_cs = variable(core, L_J_cs, L_T; lvar = 0, start = initialize_vars ? Array(solution(result, vars.p_jt_scr_cs)) : 1) + p_jt_nsc_pr = variable(core, L_J_pr, L_T; lvar = 0, start = initialize_vars ? Array(solution(result, vars.p_jt_nsc_pr)) : 1) + p_jt_rru_on_pr = variable(core, L_J_pr, L_T; lvar = 0, start = initialize_vars ? Array(solution(result, vars.p_jt_rru_on_pr)) : 1) + p_jt_rru_on_cs = variable(core, L_J_cs, L_T; lvar = 0, start = initialize_vars ? Array(solution(result, vars.p_jt_rru_on_cs)) : 1) + p_jt_rru_off_pr = variable(core, L_J_pr, L_T; lvar = 0, start = initialize_vars ? Array(solution(result, vars.p_jt_rru_off_pr)) : 1) + p_jt_rrd_on_pr = variable(core, L_J_pr, L_T; lvar = 0, start = initialize_vars ? Array(solution(result, vars.p_jt_rrd_on_pr)) : 1) + p_jt_rrd_on_cs = variable(core, L_J_cs, L_T; lvar = 0, start = initialize_vars ? Array(solution(result, vars.p_jt_rrd_on_cs)) : 1) + p_jt_rrd_off_cs = variable(core, L_J_cs, L_T; lvar = 0, start = initialize_vars ? Array(solution(result, vars.p_jt_rrd_off_cs)) : 1) + q_jt_qru_pr = variable(core, L_J_pr, L_T; lvar = 0, start = initialize_vars ? Array(solution(result, vars.q_jt_qru_pr)) : 1) + q_jt_qru_cs = variable(core, L_J_cs, L_T; lvar = 0, start = initialize_vars ? Array(solution(result, vars.q_jt_qru_cs)) : 1) + q_jt_qrd_pr = variable(core, L_J_pr, L_T; lvar = 0, start = initialize_vars ? Array(solution(result, vars.q_jt_qrd_pr)) : 1) + q_jt_qrd_cs = variable(core, L_J_cs, L_T; lvar = 0, start = initialize_vars ? Array(solution(result, vars.q_jt_qrd_cs)) : 1) + + + p_nt_rgu_req = variable(core, L_N_p, L_T; start = initialize_vars ? Array(solution(result, vars.p_nt_rgu_req)) : 1) + p_nt_rgd_req = variable(core, L_N_p, L_T; start = initialize_vars ? Array(solution(result, vars.p_nt_rgd_req)) : 1) + p_nt_scr_req = variable(core, L_N_p, L_T; start = initialize_vars ? Array(solution(result, vars.p_nt_scr_req)) : 1) + p_nt_nsc_req = variable(core, L_N_p, L_T; start = initialize_vars ? Array(solution(result, vars.p_nt_nsc_req)) : 1) + + p_jt_pr_max = variable(core, L_N_p, L_T;) + #Bounds from 4.3.1 Reserve shortfall domains (20-27) + p_nt_rgu_plus = variable(core, L_N_p, L_T; lvar = 0, start = initialize_vars ? Array(solution(result, vars.p_nt_rgu_plus)) : 1) + p_nt_rgd_plus = variable(core, L_N_p, L_T; lvar = 0, start = initialize_vars ? Array(solution(result, vars.p_nt_rgd_plus)) : 1) + p_nt_scr_plus = variable(core, L_N_p, L_T; lvar = 0, start = initialize_vars ? Array(solution(result, vars.p_nt_scr_plus)) : 1) + p_nt_nsc_plus = variable(core, L_N_p, L_T; lvar = 0, start = initialize_vars ? Array(solution(result, vars.p_nt_nsc_plus)) : 1) + p_nt_rru_plus = variable(core, L_N_p, L_T; lvar = 0, start = initialize_vars ? Array(solution(result, vars.p_nt_rru_plus)) : 1) + p_nt_rrd_plus = variable(core, L_N_p, L_T; lvar = 0, start = initialize_vars ? Array(solution(result, vars.p_nt_rrd_plus)) : 1) + q_nt_qru_plus = variable(core, L_N_q, L_T; lvar = 0, start = initialize_vars ? Array(solution(result, vars.q_nt_qru_plus)) : 1) + q_nt_qrd_plus = variable(core, L_N_q, L_T; lvar = 0, start = initialize_vars ? Array(solution(result, vars.q_nt_qrd_plus)) : 1) + + q_it = variable(core, I, L_T; start = initialize_vars ? Array(solution(result, vars.q_it)) : 1) + q_it_plus = variable(core, I, L_T; start = initialize_vars ? Array(solution(result, vars.q_it_plus)) : 1) + + #s_jt_plus split on ln and xf + #Bounds from 4.8.1 AC branch flow limits and penalties (138) + s_jt_plus_ln = variable(core, L_J_ln, L_T; lvar = 0, start = initialize_vars ? Array(solution(result, vars.s_jt_plus_ln)) : 1) + s_jt_plus_xf = variable(core, L_J_xf, L_T; lvar = 0, start = initialize_vars ? Array(solution(result, vars.s_jt_plus_xf)) : 1) + + #Bounds form 4.2.4 Bus voltage (19) + + v_it = variable(core, I, L_T; lvar = sc_data.v_lvar, uvar = sc_data.v_uvar, start = initialize_vars ? Array(solution(result, vars.v_it)) : ones(I, L_T)) + + z_w_en_max_pr = variable(core, L_W_en_max_pr; start = initialize_vars ? Array(solution(result, vars.z_w_en_max_pr)) : 0) + z_w_en_max_cs = variable(core, L_W_en_max_cs; start = initialize_vars ? Array(solution(result, vars.z_w_en_max_cs)) : 0) + z_w_en_min_pr = variable(core, L_W_en_min_pr; start = initialize_vars ? Array(solution(result, vars.z_w_en_min_pr)) : 0) + z_w_en_min_cs = variable(core, L_W_en_min_cs; start = initialize_vars ? Array(solution(result, vars.z_w_en_min_cs)) : 0) + + #split z_jt_en and on into pr and cs + z_jt_en_pr = variable(core, L_J_pr, L_T; start = initialize_vars ? Array(solution(result, vars.z_jt_en_pr)) : 0) + z_jt_en_cs = variable(core, L_J_cs, L_T; start = initialize_vars ? Array(solution(result, vars.z_jt_en_cs)) : 0) + + z_it_p = variable(core, I, L_T; start = initialize_vars ? Array(solution(result, vars.z_it_p)) : 0) + z_it_q = variable(core, I, L_T; start = initialize_vars ? Array(solution(result, vars.z_it_q)) : 0) + + #z_jt_s split into ln and xf + z_jt_s_ln = variable(core, L_J_ln, L_T; start = initialize_vars ? Array(solution(result, vars.z_jt_s_ln)) : 0) + z_jt_s_xf = variable(core, L_J_xf, L_T; start = initialize_vars ? Array(solution(result, vars.z_jt_s_xf)) : 0) + #z_jt rgu, rgd, scr, nsc, rru, rrd, qru, qrd split into pr and cs + z_jt_rgu_pr = variable(core, L_J_pr, L_T; start = initialize_vars ? Array(solution(result, vars.z_jt_rgu_pr)) : 0) + z_jt_rgu_cs = variable(core, L_J_cs, L_T; start = initialize_vars ? Array(solution(result, vars.z_jt_rgu_cs)) : 0) + z_jt_rgd_pr = variable(core, L_J_pr, L_T; start = initialize_vars ? Array(solution(result, vars.z_jt_rgd_pr)) : 0) + z_jt_rgd_cs = variable(core, L_J_cs, L_T; start = initialize_vars ? Array(solution(result, vars.z_jt_rgd_cs)) : 0) + z_jt_scr_pr = variable(core, L_J_pr, L_T; start = initialize_vars ? Array(solution(result, vars.z_jt_scr_pr)) : 0) + z_jt_scr_cs = variable(core, L_J_cs, L_T; start = initialize_vars ? Array(solution(result, vars.z_jt_scr_cs)) : 0) + z_jt_nsc_pr = variable(core, L_J_pr, L_T; start = initialize_vars ? Array(solution(result, vars.z_jt_nsc_pr)) : 0) + z_jt_rru_pr = variable(core, L_J_pr, L_T; start = initialize_vars ? Array(solution(result, vars.z_jt_rru_pr)) : 0) + z_jt_rru_cs = variable(core, L_J_cs, L_T; start = initialize_vars ? Array(solution(result, vars.z_jt_rru_cs)) : 0) + z_jt_rrd_pr = variable(core, L_J_pr, L_T; start = initialize_vars ? Array(solution(result, vars.z_jt_rrd_pr)) : 0) + z_jt_rrd_cs = variable(core, L_J_cs, L_T; start = initialize_vars ? Array(solution(result, vars.z_jt_rrd_cs)) : 0) + z_jt_qru_pr = variable(core, L_J_pr, L_T; start = initialize_vars ? Array(solution(result, vars.z_jt_qru_pr)) : 0) + z_jt_qru_cs = variable(core, L_J_cs, L_T; start = initialize_vars ? Array(solution(result, vars.z_jt_qru_cs)) : 0) + z_jt_qrd_pr = variable(core, L_J_pr, L_T; start = initialize_vars ? Array(solution(result, vars.z_jt_qrd_pr)) : 0) + z_jt_qrd_cs = variable(core, L_J_cs, L_T; start = initialize_vars ? Array(solution(result, vars.z_jt_qrd_cs)) : 0) + + z_nt_rgu = variable(core, L_N_p, L_T; start = initialize_vars ? Array(solution(result, vars.z_nt_rgu)) : 0) + z_nt_rgd = variable(core, L_N_p, L_T; start = initialize_vars ? Array(solution(result, vars.z_nt_rgd)) : 0) + z_nt_scr = variable(core, L_N_p, L_T; start = initialize_vars ? Array(solution(result, vars.z_nt_scr)) : 0) + z_nt_nsc = variable(core, L_N_p, L_T; start = initialize_vars ? Array(solution(result, vars.z_nt_nsc)) : 0) + z_nt_rru = variable(core, L_N_p, L_T; start = initialize_vars ? Array(solution(result, vars.z_nt_rru)) : 0) + z_nt_rrd = variable(core, L_N_p, L_T; start = initialize_vars ? Array(solution(result, vars.z_nt_rrd)) : 0) + z_nt_qru = variable(core, L_N_q, L_T; start = initialize_vars ? Array(solution(result, vars.z_nt_qru)) : 0) + z_nt_qrd = variable(core, L_N_q, L_T; start = initialize_vars ? Array(solution(result, vars.z_nt_qru)) : 0) + + #Bounds added so that angle doesnt blow up and cause computational errors + θ_it = variable(core, I, L_T; lvar = 0, uvar = pi + 0.0001, start = initialize_vars ? Array(solution(result, vars.θ_it)) : 0) + + #split τjt and φjt into xf only, ln is fixed + τ_jt_xf = variable(core, L_J_xf, L_T; start = initialize_vars ? Array(solution(result, vars.τ_jt_xf)) : 1) + #Bounds added so that angle doesnt blow up and cause computational errors + φ_jt_xf = variable(core, L_J_xf, L_T; lvar = 0, uvar = pi + 0.0001, start = initialize_vars ? Array(solution(result, vars.φ_jt_xf)) : 0) + - core = ExaModels.ExaCore(T; backend = backend) - - # Voltage angle and magnitudes - va = ExaModels.variable(core, 1:data.nbus, 1:data.K; start=zeros(data.nbus, data.K)) - vm = ExaModels.variable( - core, - 1:data.nbus, 1:data.K; - start = ones(data.nbus, data.K), - lvar = repeat(data.vmin, data.K), - uvar = repeat(data.vmax, data.K), - ) - # Power generations - pg = ExaModels.variable( - core, - 1:data.ngen, 1:data.K; - lvar = repeat(data.pmin, data.K), - uvar = repeat(data.pmax, data.K), - ) - qg = ExaModels.variable( - core, - 1:data.ngen, 1:data.K; - lvar = repeat(data.qmin, data.K), - uvar = repeat(data.qmax, data.K), - ) - # Power injection - p = ExaModels.variable( - core, - 1:data.narc, 1:data.K; - lvar = -data.max_inj, - uvar = data.max_inj, - ) - q = ExaModels.variable( - core, - 1:data.narc, 1:data.K; - lvar = -data.max_inj, - uvar = data.max_inj, - ) - - # Objective (cost for base case only) - obj = ExaModels.objective( - core, - g.cost1 * pg[g.i, 1]^2 + g.cost2 * pg[g.i, 1] + g.cost3 for g in data.gen - ) - - # Constraints - # Voltage angle at reference buses is 0 - c1 = ExaModels.constraint(core, va[i, k] for (i, k) in data.idx_ref) - # Branch power injection - # from, w.r.t active power - c2 = ExaModels.constraint( - core, - p[b.f_idx, k] - b.c5 * vm[b.f_bus, k]^2 - - b.c3 * (vm[b.f_bus, k] * vm[b.t_bus, k] * cos(va[b.f_bus, k] - va[b.t_bus, k])) - - b.c4 * (vm[b.f_bus, k] * vm[b.t_bus, k] * sin(va[b.f_bus, k] - va[b.t_bus, k])) for - (b, k) in data.idx_branch - ) - # from, w.r.t reactive power - c3 = ExaModels.constraint( - core, - q[b.f_idx, k] + - b.c6 * vm[b.f_bus, k]^2 + - b.c4 * (vm[b.f_bus, k] * vm[b.t_bus, k] * cos(va[b.f_bus, k] - va[b.t_bus, k])) - - b.c3 * (vm[b.f_bus, k] * vm[b.t_bus, k] * sin(va[b.f_bus, k] - va[b.t_bus, k])) for - (b, k) in data.idx_branch - ) - # to, w.r.t active power - c4 = ExaModels.constraint( - core, - p[b.t_idx, k] - b.c7 * vm[b.t_bus, k]^2 - - b.c1 * (vm[b.t_bus, k] * vm[b.f_bus, k] * cos(va[b.t_bus, k] - va[b.f_bus, k])) - - b.c2 * (vm[b.t_bus, k] * vm[b.f_bus, k] * sin(va[b.t_bus, k] - va[b.f_bus, k])) for - (b, k) in data.idx_branch - ) - # to, w.r.t reactive power - c5 = ExaModels.constraint( - core, - q[b.t_idx, k] + - b.c8 * vm[b.t_bus, k]^2 + - b.c2 * (vm[b.t_bus, k] * vm[b.f_bus, k] * cos(va[b.t_bus, k] - va[b.f_bus, k])) - - b.c1 * (vm[b.t_bus, k] * vm[b.f_bus, k] * sin(va[b.t_bus, k] - va[b.f_bus, k])) for - (b, k) in data.idx_branch - ) - # Difference of angles - c6 = ExaModels.constraint( - core, - va[b.f_bus, k] - va[b.t_bus, k] for (b, k) in data.idx_branch; - lcon = repeat(data.angmin, data.K), - ucon = repeat(data.angmax, data.K), - ) - # Line-flow constraints - # from - lcon = fill!(similar(data.branch, Float64, data.K * length(data.branch)), -Inf) - c7 = ExaModels.constraint( - core, - p[b.f_idx, k]^2 + q[b.f_idx, k]^2 - b.rate_a_sq for (b, k) in data.idx_branch; - lcon = lcon, - ) - # to - c8 = ExaModels.constraint( - core, - p[b.t_idx, k]^2 + q[b.t_idx, k]^2 - b.rate_a_sq for (b, k) in data.idx_branch; - lcon = lcon, - ) - - # Power flow constraints - c9 = ExaModels.constraint(core, b.pd + b.gs * vm[b.i, k]^2 for (b, k) in data.idx_bus) - c10 = ExaModels.constraint(core, b.qd - b.bs * vm[b.i, k]^2 for (b, k) in data.idx_bus) - # TODO: shifts of index is hacky. - c11 = ExaModels.constraint!(core, c9, a.bus + (k-1)*data.nbus => p[a.i, k] for (a, k) in data.idx_arc) - c12 = ExaModels.constraint!(core, c10, a.bus + (k-1)*data.nbus => q[a.i, k] for (a, k) in data.idx_arc) - c13 = ExaModels.constraint!(core, c9, g.bus + (k-1)*data.nbus => -pg[g.i, k] for (g, k) in data.idx_gen) - c14 = ExaModels.constraint!(core, c10, g.bus + (k-1)*data.nbus=> -qg[g.i, k] for (g, k) in data.idx_gen) - - # Corrective OPF formulation - c_corrective_gen = ExaModels.constraint( - core, - pg[g.i, k] - pg[g.i, 1] for (g, k) in data.idx_gen_cont; - lcon = repeat(-data.Δp, data.K), - ucon = repeat(data.Δp, data.K), - ) - c_corrective_vmag = ExaModels.constraint( - core, - vm[b.i, k] - vm[b.i, 1] for (b, k) in data.idx_bus_cont; - ) - - model =ExaModels.ExaModel(core; kwargs...) + if include_ctg + #Bound from 4.9.1 Penalty on post-contingency AC branch overload (157) + s_jtk_plus_ln = variable(core, length(sc_data.jtk_ln_flattened); lvar = 0, start = initialize_vars ? Array(solution(result, vars.s_jtk_plus_ln)) : 0) + s_jtk_plus_xf = variable(core, length(sc_data.jtk_xf_flattened); lvar = 0, start = initialize_vars ? Array(solution(result, vars.s_jtk_plus_xf)) : 0) + z_jtk_s_ln = variable(core, length(sc_data.jtk_ln_flattened); start = initialize_vars ? Array(solution(result, vars.z_jtk_s_ln)) : 0) + z_jtk_s_xf = variable(core, length(sc_data.jtk_xf_flattened); start = initialize_vars ? Array(solution(result, vars.z_jtk_s_xf)) : 0) + p_jtk_ln = variable(core, length(sc_data.jtk_ln_flattened); start = initialize_vars ? Array(solution(result, vars.p_jtk_ln)) : 1) + p_jtk_xf = variable(core, length(sc_data.jtk_xf_flattened); start = initialize_vars ? Array(solution(result, vars.p_jtk_xf)) : 1) + θ_itk = variable(core, I, L_T, K; start = initialize_vars ? Array(solution(result, vars.θ_itk)) : 0) + p_t_sl = variable(core, L_T; start = initialize_vars ? Array(solution(result, vars.p_t_sl)) : 1) + z_tk_ctg = variable(core, L_T, K; start = initialize_vars ? Array(solution(result, vars.z_tk_ctg)) : 0) + z_t_ctg_min = variable(core, L_T; start = initialize_vars ? Array(solution(result, vars.z_t_ctg_min)) : 0) + z_t_ctg_avg = variable(core, L_T; start = initialize_vars ? Array(solution(result, vars.z_t_ctg_avg)) : 0) + end - vars = ( - va = va, - vm = vm, - pg = pg, - qg = qg, - p = p, - q = q - ) - - return model, vars + + #4.1 Market surplus objective + #constraint (6-9) + #All objectives are negative so that we can minimize + + #Removing all uc variables, which include z_on, z_su, z_sd, z_sus. Note that this means objective value from ExaModelsPower will not match the objective calculated for the full GOC3 model, even though we still are solving the same problem (with UC as constants) + + if include_ctg + o2 = objective(core, -z_t_ctg_min[t] for t in sc_data.periods) + o3 = objective(core, -z_t_ctg_avg[t] for t in sc_data.periods) + end + + o6_t_pr = objective(core, -(-z_jt_en_pr[pr.j_pr, pr.t] + - (z_jt_rgu_pr[pr.j_pr, pr.t] + z_jt_rgd_pr[pr.j_pr, pr.t] + z_jt_scr_pr[pr.j_pr, pr.t] + z_jt_nsc_pr[pr.j_pr, pr.t] + z_jt_rru_pr[pr.j_pr, pr.t] + + z_jt_rrd_pr[pr.j_pr, pr.t] + z_jt_qru_pr[pr.j_pr, pr.t] +z_jt_qrd_pr[pr.j_pr, pr.t])) for pr in sc_data.prarray) + o6_t_cs = objective(core, -(z_jt_en_cs[cs.j_cs, cs.t] + - (z_jt_rgu_cs[cs.j_cs, cs.t] + z_jt_rgd_cs[cs.j_cs, cs.t] + z_jt_scr_cs[cs.j_cs, cs.t] + z_jt_rru_cs[cs.j_cs, cs.t] + + z_jt_rrd_cs[cs.j_cs, cs.t] + z_jt_qru_cs[cs.j_cs, cs.t] +z_jt_qrd_cs[cs.j_cs, cs.t])) for cs in sc_data.csarray) + o6_t_ln = objective(core, -(- z_jt_s_ln[ln.j_ln, ln.t]) for ln in sc_data.aclbrancharray) + o6_t_xf = objective(core, -(- z_jt_s_xf[xf.j_xf, xf.t]) for xf in sc_data.acxbrancharray) + o6_t_i = objective(core, -(-(z_it_p[b.i, b.t] + z_it_q[b.i, b.t])) for b in sc_data.busarray) + o6_t_Np = objective(core, -(-(z_nt_rgu[n.n_p, n.t] + z_nt_rgd[n.n_p, n.t] + z_nt_scr[n.n_p, n.t] + z_nt_nsc[n.n_p, n.t] + z_nt_rru[n.n_p, n.t] + z_nt_rrd[n.n_p, n.t])) for n in sc_data.preservearray) + o6_t_Nq = objective(core, -(-(z_nt_qru[n.n_q, n.t] + z_nt_qrd[n.n_q, n.t])) for n in sc_data.qreservearray) + if L_W_en_max_pr > 0 + o6_en_max_pr = objective(core, z_w_en_max_pr[w] for w in 1:L_W_en_max_pr) + end + if L_W_en_max_cs > 0 + o6_en_max_cs = objective(core, z_w_en_max_cs[w] for w in 1:L_W_en_max_cs) + end + if L_W_en_min_pr > 0 + o6_en_min_pr = objective(core, z_w_en_min_pr[w] for w in 1:L_W_en_min_pr) + end + if L_W_en_min_cs > 0 + o6_en_min_cs = objective(core, z_w_en_min_cs[w] for w in 1:L_W_en_min_cs) + end + + + + if include_ctg + if K>0 + c4 = constraint(core, z_t_ctg_min[ind.t] - z_tk_ctg[ind.t, ind.k] for ind in sc_data.tk_index; lcon = fill(-Inf, size(sc_data.tk_index))) + c5 = constraint(core, z_t_ctg_avg[t]*K for t in sc_data.periods) + c5_a = constraint!(core, c5, ind.t => -z_tk_ctg[ind.t, ind.k] for ind in sc_data.tk_index) + else + c4 = constraint(core, z_t_ctg_min[t] for t in sc_data.periods) + c5 = constraint(core, z_t_ctg_avg[t] for t in sc_data.periods) + end + c10 = constraint(core, z_tk_ctg[ind.t, ind.k] for ind in sc_data.tk_index) + c10_ln = constraint!(core, c10, ln.t + L_T*(ln.ctg-1) => z_jtk_s_ln[ln.flat_jtk_ln] for ln in sc_data.jtk_ln_flattened) + c10_xf = constraint!(core, c10, xf.t + L_T*(xf.ctg-1) => z_jtk_s_xf[xf.flat_jtk_xf] for xf in sc_data.jtk_xf_flattened) + end + + #4.2.1 Bus power mismatch and penalized mismatch definitions + c_11 = constraint(core, p_it_plus[b.i, b.t] - p_it[b.i, b.t] for b in sc_data.busarray; ucon = fill(Inf, size(sc_data.busarray))) + c_12 = constraint(core, p_it_plus[b.i, b.t] + p_it[b.i, b.t] for b in sc_data.busarray; ucon = fill(Inf, size(sc_data.busarray))) + c_13 = constraint(core, q_it_plus[b.i, b.t] - q_it[b.i, b.t] for b in sc_data.busarray; ucon = fill(Inf, size(sc_data.busarray))) + c_14 = constraint(core, q_it_plus[b.i, b.t] + q_it[b.i, b.t] for b in sc_data.busarray; ucon = fill(Inf, size(sc_data.busarray))) + #4.2.2 Bus pwoer mismatch penalty + c15 = constraint(core, z_it_p[b.i, b.t] - b.dt*sum(sc_data.c_p)*p_it_plus[b.i, b.t] for b in sc_data.busarray) + c16 = constraint(core, z_it_q[b.i, b.t] - b.dt*sum(sc_data.c_q)*q_it_plus[b.i, b.t] for b in sc_data.busarray) + #4.2.3 Bus real and reactive power balance + c17 = constraint(core, p_it[b.i, b.t] for b in sc_data.busarray) + c17_pr = constraint!(core, c17, pr.bus + I*(pr.t-1) => p_jt_pr[pr.j_pr, pr.t] for pr in sc_data.prarray) + c17_cs = constraint!(core, c17, cs.bus + I*(cs.t-1) => -p_jt_cs[cs.j_cs, cs.t] for cs in sc_data.csarray) + c17_sh = constraint!(core, c17, sh.bus + I*(sh.t-1) => -p_jt_sh[sh.j_sh, sh.t] for sh in sc_data.shuntarray) + + #Reminder: fr and to split for ln, xf, and dc + c17_fr_ln = constraint!(core, c17, ln.fr_bus + I*(ln.t-1) => -p_jt_fr_ln[ln.j_ln, ln.t] for ln in sc_data.aclbrancharray) + c17_fr_xf = constraint!(core, c17, xf.fr_bus + I*(xf.t-1) => -p_jt_fr_xf[xf.j_xf, xf.t] for xf in sc_data.acxbrancharray) + c17_fr_dc = constraint!(core, c17, dc.fr_bus + I*(dc.t-1) => -p_jt_fr_dc[dc.j_dc, dc.t] for dc in sc_data.dclinearray) + c17_to_ln = constraint!(core, c17, ln.to_bus + I*(ln.t-1) => -p_jt_to_ln[ln.j_ln, ln.t] for ln in sc_data.aclbrancharray) + c17_to_xf = constraint!(core, c17, xf.to_bus + I*(xf.t-1) => -p_jt_to_xf[xf.j_xf, xf.t] for xf in sc_data.acxbrancharray) + c17_to_dc = constraint!(core, c17, dc.to_bus + I*(dc.t-1) => -p_jt_to_dc[dc.j_dc, dc.t] for dc in sc_data.dclinearray) + c18 = constraint(core, q_it[b.i, b.t] for b in sc_data.busarray) + c18_pr = constraint!(core, c18, pr.bus + I*(pr.t-1) => q_jt_pr[pr.j_pr, pr.t] for pr in sc_data.prarray) + c18_cs = constraint!(core, c18, cs.bus + I*(cs.t-1) => -q_jt_cs[cs.j_cs, cs.t] for cs in sc_data.csarray) + c18_sh = constraint!(core, c18, sh.bus + I*(sh.t-1) => -q_jt_sh[sh.j_sh, sh.t] for sh in sc_data.shuntarray) + #Reminder: fr and to split for ln, xf, and dc + c18_fr_ln = constraint!(core, c18, ln.fr_bus + I*(ln.t-1) => -q_jt_fr_ln[ln.j_ln, ln.t] for ln in sc_data.aclbrancharray) + c18_fr_xf = constraint!(core, c18, xf.fr_bus + I*(xf.t-1) => -q_jt_fr_xf[xf.j_xf, xf.t] for xf in sc_data.acxbrancharray) + c18_fr_dc = constraint!(core, c18, dc.fr_bus + I*(dc.t-1) => -q_jt_fr_dc[dc.j_dc, dc.t] for dc in sc_data.dclinearray) + c18_to_ln = constraint!(core, c18, ln.to_bus + I*(ln.t-1) => -q_jt_to_ln[ln.j_ln, ln.t] for ln in sc_data.aclbrancharray) + c18_to_xf = constraint!(core, c18, xf.to_bus + I*(xf.t-1) => -q_jt_to_xf[xf.j_xf, xf.t] for xf in sc_data.acxbrancharray) + c18_to_dc = constraint!(core, c18, dc.to_bus + I*(dc.t-1) => -q_jt_to_dc[dc.j_dc, dc.t] for dc in sc_data.dclinearray) + + + #4.3.2 Reserve shortfall penalties + c28 = constraint(core, z_nt_rgu[n.n_p, n.t] - n.dt*n.c_rgu*p_nt_rgu_plus[n.n_p, n.t] for n in sc_data.preservearray) + c29 = constraint(core, z_nt_rgd[n.n_p, n.t] - n.dt*n.c_rgd*p_nt_rgd_plus[n.n_p, n.t] for n in sc_data.preservearray) + c30 = constraint(core, z_nt_scr[n.n_p, n.t] - n.dt*n.c_scr*p_nt_scr_plus[n.n_p, n.t] for n in sc_data.preservearray) + c31 = constraint(core, z_nt_nsc[n.n_p, n.t] - n.dt*n.c_nsc*p_nt_nsc_plus[n.n_p, n.t] for n in sc_data.preservearray) + c32 = constraint(core, z_nt_rru[n.n_p, n.t] - n.dt*n.c_rru*p_nt_rru_plus[n.n_p, n.t] for n in sc_data.preservearray) + c33 = constraint(core, z_nt_rrd[n.n_p, n.t] - n.dt*n.c_rrd*p_nt_rrd_plus[n.n_p, n.t] for n in sc_data.preservearray) + c34 = constraint(core, z_nt_qru[n.n_q, n.t] - n.dt*n.c_qru*q_nt_qru_plus[n.n_q, n.t] for n in sc_data.qreservearray) + c35 = constraint(core, z_nt_qrd[n.n_q, n.t] - n.dt*n.c_qrd*q_nt_qrd_plus[n.n_q, n.t] for n in sc_data.qreservearray) + + + #4.3.3 Reserve requirements + c36 = constraint(core, p_nt_rgu_req[n.n_p, n.t]/n.σ_rgu for n in sc_data.preservearray) + c36_cs = constraint!(core, c36, cs.n + L_N_p*(cs.t-1) => -p_jt_cs[cs.j_cs, cs.t] for cs in sc_data.preservesetarray_cs) + c37 = constraint(core, p_nt_rgd_req[n.n_p, n.t]/n.σ_rgd for n in sc_data.preservearray) + c37_cs = constraint!(core, c37, cs.n + L_N_p*(cs.t-1) => -p_jt_cs[cs.j_cs, cs.t] for cs in sc_data.preservesetarray_cs) + #assuming c_scr and c_nsc are always positive + cmax38 = constraint(core, p_jt_pr_max[pr.n_p, pr.t] - p_jt_pr[pr.j_pr, pr.t] for pr in sc_data.preservesetarray_pr; ucon = fill(Inf, size(sc_data.prarray))) + c38 = constraint(core, p_nt_scr_req[n.n_p, n.t] - n.σ_scr*p_jt_pr_max[n.n_p, n.t] for n in sc_data.preservearray) + c39 = constraint(core, p_nt_nsc_req[n.n_p, n.t] - n.σ_nsc*p_jt_pr_max[n.n_p, n.t] for n in sc_data.preservearray) + + #4.3.4 Reserve balance + #Reminder, p and q sets have been split up for pr and cs + c40 = constraint(core, p_nt_rgu_plus[n.n_p, n.t] - p_nt_rgu_req[n.n_p, n.t] for n in sc_data.preservearray; ucon = fill(Inf, size(sc_data.preservearray))) + c40_pr = constraint!(core, c40, pr.n_p + L_N_p*(pr.t-1) => p_jt_rgu_pr[pr.j_pr, pr.t] for pr in sc_data.preservesetarray_pr) + c40_cs = constraint!(core, c40, cs.n_p + L_N_p*(cs.t-1) => p_jt_rgu_cs[cs.j_cs, cs.t] for cs in sc_data.preservesetarray_cs) + c41 = constraint(core, p_nt_rgd_plus[n.n_p, n.t] - p_nt_rgd_req[n.n_p, n.t] for n in sc_data.preservearray; ucon = fill(Inf, size(sc_data.preservearray))) + c41_pr = constraint!(core, c41, pr.n_p + L_N_p*(pr.t-1) => p_jt_rgd_pr[pr.j_pr, pr.t] for pr in sc_data.preservesetarray_pr) + c41_cs = constraint!(core, c41, cs.n_p + L_N_p*(cs.t-1) => p_jt_rgd_cs[cs.j_cs, cs.t] for cs in sc_data.preservesetarray_cs) + c42 = constraint(core, p_nt_scr_plus[n.n_p, n.t] - p_nt_rgu_req[n.n_p, n.t] - p_nt_scr_req[n.n_p, n.t] for n in sc_data.preservearray; ucon = fill(Inf, size(sc_data.preservearray))) + c42_pr = constraint!(core, c42, pr.n_p + L_N_p*(pr.t-1) => p_jt_rgu_pr[pr.j_pr, pr.t] + p_jt_scr_pr[pr.j_pr, pr.t] for pr in sc_data.preservesetarray_pr) + c42_cs = constraint!(core, c42, cs.n_p + L_N_p*(cs.t-1) => p_jt_rgu_cs[cs.j_cs, cs.t] + p_jt_scr_cs[cs.j_cs, cs.t] for cs in sc_data.preservesetarray_cs) + c43 = constraint(core, p_nt_nsc_plus[n.n_p, n.t] - p_nt_rgu_req[n.n_p, n.t] - p_nt_scr_req[n.n_p, n.t] - p_nt_nsc_req[n.n_p, n.t] for n in sc_data.preservearray; ucon = fill(Inf, size(sc_data.preservearray))) + c43_pr = constraint!(core, c43, pr.n_p + L_N_p*(pr.t-1) => p_jt_rgu_pr[pr.j_pr, pr.t] + p_jt_scr_pr[pr.j_pr, pr.t] + p_jt_nsc_pr[pr.j_pr, pr.t] for pr in sc_data.preservesetarray_pr) + c43_cs = constraint!(core, c43, cs.n_p + L_N_p*(cs.t-1) => p_jt_rgu_cs[cs.j_cs, cs.t] + p_jt_scr_cs[cs.j_cs, cs.t] for cs in sc_data.preservesetarray_cs) + c44 = constraint(core, p_nt_rru_plus[n.n_p, n.t] - n.p_rru_min for n in sc_data.preservearray; ucon = fill(Inf, size(sc_data.preservearray))) + c44_pr = constraint!(core, c44, pr.n_p + L_N_p*(pr.t-1) => p_jt_rru_on_pr[pr.j_pr, pr.t] + p_jt_rru_off_pr[pr.j_pr, pr.t] for pr in sc_data.preservesetarray_pr) + c44_cs = constraint!(core, c44, cs.n_p + L_N_p*(cs.t-1) => p_jt_rru_on_cs[cs.j_cs, cs.t] for cs in sc_data.preservesetarray_cs) + c45 = constraint(core, p_nt_rrd_plus[n.n_p, n.t] - n.p_rrd_min for n in sc_data.preservearray; ucon = fill(Inf, size(sc_data.preservearray))) + c45_pr = constraint!(core, c45, pr.n_p + L_N_p*(pr.t-1) => p_jt_rrd_on_pr[pr.j_pr, pr.t] for pr in sc_data.preservesetarray_pr) + c45_cs = constraint!(core, c45, cs.n_p + L_N_p*(cs.t-1) => p_jt_rrd_on_cs[cs.j_cs, cs.t] + p_jt_rrd_off_cs[cs.j_cs, cs.t] for cs in sc_data.preservesetarray_cs) + c46 = constraint(core, q_nt_qru_plus[n.n_q, n.t] - n.q_qru_min for n in sc_data.qreservearray; ucon = fill(Inf, size(sc_data.qreservearray))) + c46_pr = constraint!(core, c46, pr.n_q + L_N_q*(pr.t-1) => q_jt_qru_pr[pr.j_pr, pr.t] for pr in sc_data.qreservesetarray_pr) + c46_cs = constraint!(core, c46, cs.n_q + L_N_q*(cs.t-1) => q_jt_qru_cs[cs.j_cs, cs.t] for cs in sc_data.qreservesetarray_cs) + c47 = constraint(core, q_nt_qrd_plus[n.n_q, n.t] - n.q_qrd_min for n in sc_data.qreservearray; ucon = fill(Inf, size(sc_data.qreservearray))) + c47_pr = constraint!(core, c47, pr.n_q + L_N_q*(pr.t-1) => q_jt_qrd_pr[pr.j_pr, pr.t] for pr in sc_data.qreservesetarray_pr) + c47_cs = constraint!(core, c47, cs.n_q + L_N_q*(cs.t-1) => q_jt_qrd_cs[cs.j_cs, cs.t] for cs in sc_data.qreservesetarray_cs) + + #skipping constraints 48-67. Assume unit commitment is always satisfied + + #4.6.1 Producing and consuming device startup, shutdown, and dispatchable power + #p_jt variants are split on pr and cs + c68_pr = constraint(core, p_jt_on_pr[pr.j_pr, pr.t] + p_jt_su_pr[pr.j_pr, pr.t] + p_jt_sd_pr[pr.j_pr, pr.t] - p_jt_pr[pr.j_pr, pr.t] for pr in sc_data.prarray) + c68_cs = constraint(core, p_jt_on_cs[cs.j_cs, cs.t] + p_jt_su_cs[cs.j_cs, cs.t] + p_jt_sd_cs[cs.j_cs, cs.t] - p_jt_cs[cs.j_cs, cs.t] for cs in sc_data.csarray) + c69_pr = constraint(core, p_jt_su_pr[pr.j_pr, pr.t] - pr.sum_T_supc_pr_jt for pr in sc_data.prarray) + c69_cs = constraint(core, p_jt_su_cs[cs.j_cs, cs.t] - cs.sum_T_supc_cs_jt for cs in sc_data.csarray) + c70_pr = constraint(core, p_jt_sd_pr[pr.j_pr, pr.t] - pr.sum_T_sdpc_pr_jt for pr in sc_data.prarray) + c70_cs = constraint(core, p_jt_sd_cs[cs.j_cs, cs.t] - cs.sum_T_sdpc_cs_jt for cs in sc_data.csarray) + + #4.6.2 Ramping limits + #p split for pr and cs + c71_pr = constraint(core, p_jt_pr[pr.j_pr, pr.t] - pr.p_0 - pr.dt*(pr.p_ru*(pr.u_on-pr.u_su) + pr.p_ru_su*(pr.u_su+1-pr.u_on)) for pr in sc_data.prarray[1:L_J_pr]; + lcon = fill(-Inf, size(sc_data.prarray[1:L_J_pr]))) + c71_cs = constraint(core, p_jt_cs[cs.j_cs, cs.t] - cs.p_0 - cs.dt*(cs.p_ru*(cs.u_on-cs.u_su) + cs.p_ru_su*(cs.u_su+1-cs.u_on)) for cs in sc_data.csarray[1:L_J_cs]; + lcon = fill(-Inf, size(sc_data.csarray[1:L_J_cs]))) + c72_pr = constraint(core, p_jt_pr[pr.j_pr, pr.t] - p_jt_pr[pr.j_pr, pr.t-1] - pr.dt*(pr.p_ru*(pr.u_on-pr.u_su) + pr.p_ru_su*(pr.u_su+1-pr.u_on)) for pr in sc_data.prarray[L_J_pr+1:end]; + lcon = fill(-Inf, size(sc_data.prarray[L_J_pr+1:end]))) + c72_cs = constraint(core, p_jt_cs[cs.j_cs, cs.t] - p_jt_cs[cs.j_cs, cs.t-1] - cs.dt*(cs.p_ru*(cs.u_on-cs.u_su) + cs.p_ru_su*(cs.u_su+1-cs.u_on)) for cs in sc_data.csarray[L_J_cs+1:end]; + lcon = fill(-Inf, size(sc_data.csarray[L_J_cs+1:end]))) + c73_pr = constraint(core, p_jt_pr[pr.j_pr, pr.t] - pr.p_0 + pr.dt*(pr.p_rd*pr.u_on+pr.p_rd_sd*(1-pr.u_on)) for pr in sc_data.prarray[1:L_J_pr]; + ucon = fill(Inf, size(sc_data.prarray[1:L_J_pr]))) + c73_cs = constraint(core, p_jt_cs[cs.j_cs, cs.t] - cs.p_0 + cs.dt*(cs.p_rd*cs.u_on+cs.p_rd_sd*(1-cs.u_on)) for cs in sc_data.csarray[1:L_J_cs]; + ucon = fill(Inf, size(sc_data.csarray[1:L_J_cs]))) + c74_pr = constraint(core, p_jt_pr[pr.j_pr, pr.t] - p_jt_pr[pr.j_pr, pr.t-1] + pr.dt*(pr.p_rd*pr.u_on+pr.p_rd_sd*(1-pr.u_on)) for pr in sc_data.prarray[L_J_pr+1:end]; + ucon = fill(Inf, size(sc_data.prarray[L_J_pr+1:end]))) + c74_cs = constraint(core, p_jt_cs[cs.j_cs, cs.t] - p_jt_cs[cs.j_cs, cs.t-1] + cs.dt*(cs.p_rd*cs.u_on+cs.p_rd_sd*(1-cs.u_on)) for cs in sc_data.csarray[L_J_cs+1:end]; + ucon = fill(Inf, size(sc_data.csarray[L_J_cs+1:end]))) + + #4.6.3 Maximum/minimum energy over multiple intervals + #J_pr,cs has been split for pr and cs + c75_pr = constraint(core, e_w_plus_max_pr[w.w_en_max_pr_ind] + w.e_max for w in sc_data.W_en_max_pr; ucon = fill(Inf, size(sc_data.W_en_max_pr))) + c75_pr_a = constraint!(core, c75_pr, t.w_en_max_pr_ind => -t.dt*p_jt_pr[t.j_pr, t.t] for t in sc_data.T_w_en_max_pr) + c75_cs = constraint(core, e_w_plus_max_cs[w.w_en_max_cs_ind] + w.e_max for w in sc_data.W_en_max_cs; ucon = fill(Inf, size(sc_data.W_en_max_cs))) + c75_cs_a = constraint!(core, c75_cs, t.w_en_max_cs_ind => -t.dt*p_jt_cs[t.j_cs, t.t] for t in sc_data.T_w_en_max_cs) + c76_pr = constraint(core, e_w_plus_min_pr[w.w_en_min_pr_ind] - w.e_min for w in sc_data.W_en_min_pr; lcon = fill(-Inf, size(sc_data.W_en_min_pr))) + c76_pr_a = constraint!(core, c76_pr, t.w_en_min_pr_ind => -t.dt*p_jt_pr[t.j_pr, t.t] for t in sc_data.T_w_en_min_pr) + c76_cs = constraint(core, e_w_plus_min_cs[w.w_en_min_cs_ind] - w.e_min for w in sc_data.W_en_min_cs; lcon = fill(-Inf, size(sc_data.W_en_min_cs))) + c76_cs_a = constraint!(core, c76_cs, t.w_en_min_cs_ind => -t.dt*p_jt_cs[t.j_cs, t.t] for t in sc_data.T_w_en_min_cs) + c78_pr = constraint(core, z_w_en_max_pr[w.w_en_max_pr_ind] - sum(sc_data.c_e)*e_w_plus_max_pr[w.w_en_max_pr_ind] for w in sc_data.W_en_max_pr) + c78_cs = constraint(core, z_w_en_max_cs[w.w_en_max_cs_ind] - sum(sc_data.c_e)*e_w_plus_max_cs[w.w_en_max_cs_ind] for w in sc_data.W_en_max_cs) + c79_pr = constraint(core, z_w_en_min_pr[w.w_en_min_pr_ind] - sum(sc_data.c_e)*e_w_plus_min_pr[w.w_en_min_pr_ind] for w in sc_data.W_en_min_pr) + c79_cs = constraint(core, z_w_en_min_cs[w.w_en_min_cs_ind] - sum(sc_data.c_e)*e_w_plus_min_cs[w.w_en_min_cs_ind] for w in sc_data.W_en_min_cs) + + #4.6.5 Device reserve costs + #p_jt split into pr and cs + c90_pr = constraint(core, z_jt_rgu_pr[pr.j_pr, pr.t] - pr.dt*pr.c_rgu*p_jt_rgu_pr[pr.j_pr, pr.t] for pr in sc_data.prarray) + c90_cs = constraint(core, z_jt_rgu_cs[cs.j_cs, cs.t] - cs.dt*cs.c_rgu*p_jt_rgu_cs[cs.j_cs, cs.t] for cs in sc_data.csarray) + c91_pr = constraint(core, z_jt_rgd_pr[pr.j_pr, pr.t] - pr.dt*pr.c_rgd*p_jt_rgd_pr[pr.j_pr, pr.t] for pr in sc_data.prarray) + c91_cs = constraint(core, z_jt_rgd_cs[cs.j_cs, cs.t] - cs.dt*cs.c_rgd*p_jt_rgd_cs[cs.j_cs, cs.t] for cs in sc_data.csarray) + c92_pr = constraint(core, z_jt_scr_pr[pr.j_pr, pr.t] - pr.dt*pr.c_scr*p_jt_scr_pr[pr.j_pr, pr.t] for pr in sc_data.prarray) + c92_cs = constraint(core, z_jt_scr_cs[cs.j_cs, cs.t] - cs.dt*cs.c_scr*p_jt_scr_cs[cs.j_cs, cs.t] for cs in sc_data.csarray) + c93_pr = constraint(core, z_jt_nsc_pr[pr.j_pr, pr.t] - pr.dt*pr.c_nsc*p_jt_nsc_pr[pr.j_pr, pr.t] for pr in sc_data.prarray) + #due to c106, z_jt_nsc_cs = 0. This is realized in o6 + c94_pr = constraint(core, z_jt_rru_pr[pr.j_pr, pr.t] - pr.dt*(pr.c_rru_on*p_jt_rru_on_pr[pr.j_pr, pr.t] + pr.c_rru_off*p_jt_rru_off_pr[pr.j_pr, pr.t]) for pr in sc_data.prarray) + c94_cs = constraint(core, z_jt_rru_cs[cs.j_cs, cs.t] - cs.dt*(cs.c_rru_on*p_jt_rru_on_cs[cs.j_cs, cs.t]) for cs in sc_data.csarray) + c95_pr = constraint(core, z_jt_rrd_pr[pr.j_pr, pr.t] - pr.dt*(pr.c_rrd_on*p_jt_rrd_on_pr[pr.j_pr, pr.t]) for pr in sc_data.prarray) + c95_cs = constraint(core, z_jt_rrd_cs[cs.j_cs, cs.t] - cs.dt*(cs.c_rrd_on*p_jt_rrd_on_cs[cs.j_cs, cs.t] + cs.c_rrd_off*p_jt_rrd_off_cs[cs.j_cs, cs.t]) for cs in sc_data.csarray) + c96_pr = constraint(core, z_jt_qru_pr[pr.j_pr, pr.t] - pr.dt*pr.c_qru*q_jt_qru_pr[pr.j_pr, pr.t] for pr in sc_data.prarray) + c96_cs = constraint(core, z_jt_qru_cs[cs.j_cs, cs.t] - cs.dt*cs.c_qru*q_jt_qru_cs[cs.j_cs, cs.t] for cs in sc_data.csarray) + c97_pr = constraint(core, z_jt_qrd_pr[pr.j_pr, pr.t] - pr.dt*pr.c_qrd*q_jt_qrd_pr[pr.j_pr, pr.t] for pr in sc_data.prarray) + c97_cs = constraint(core, z_jt_qrd_cs[cs.j_cs, cs.t] - cs.dt*cs.c_qrd*q_jt_qrd_cs[cs.j_cs, cs.t] for cs in sc_data.csarray) + + #4.6.6 Absolute reserve limits, based on ramp rates + c98_pr = constraint(core, p_jt_rgu_pr[pr.j_pr, pr.t] - pr.p_rgu_max*pr.u_on for pr in sc_data.prarray; lcon = fill(-Inf, size(sc_data.prarray))) + c98_cs = constraint(core, p_jt_rgu_cs[cs.j_cs, cs.t] - cs.p_rgu_max*cs.u_on for cs in sc_data.csarray; lcon = fill(-Inf, size(sc_data.csarray))) + c99_pr = constraint(core, p_jt_rgd_pr[pr.j_pr, pr.t] - pr.p_rgd_max*pr.u_on for pr in sc_data.prarray; lcon = fill(-Inf, size(sc_data.prarray))) + c99_cs = constraint(core, p_jt_rgd_cs[cs.j_cs, cs.t] - cs.p_rgd_max*cs.u_on for cs in sc_data.csarray; lcon = fill(-Inf, size(sc_data.csarray))) + c100_pr = constraint(core, p_jt_rgu_pr[pr.j_pr, pr.t] + p_jt_scr_pr[pr.j_pr, pr.t] - pr.p_scr_max*pr.u_on for pr in sc_data.prarray; lcon = fill(-Inf, size(sc_data.prarray))) + c100_cs = constraint(core, p_jt_rgu_cs[cs.j_cs, cs.t] + p_jt_scr_cs[cs.j_cs, cs.t] - cs.p_scr_max*cs.u_on for cs in sc_data.csarray; lcon = fill(-Inf, size(sc_data.csarray))) + c101_pr = constraint(core, p_jt_nsc_pr[pr.j_pr, pr.t] - pr.p_nsc_max*(1-pr.u_on) for pr in sc_data.prarray; lcon = fill(-Inf, size(sc_data.prarray))) + #due to c107, c101_cs is not needed (assume unit commitment satisfies it) + c102_pr = constraint(core, p_jt_rgu_pr[pr.j_pr, pr.t] + p_jt_scr_pr[pr.j_pr, pr.t] + p_jt_rru_on_pr[pr.j_pr, pr.t] - pr.p_rru_on_max*pr.u_on for pr in sc_data.prarray; lcon = fill(-Inf, size(sc_data.prarray))) + c102_cs = constraint(core, p_jt_rgu_cs[cs.j_cs, cs.t] + p_jt_scr_cs[cs.j_cs, cs.t] + p_jt_rru_on_cs[cs.j_cs, cs.t] - cs.p_rru_on_max*cs.u_on for cs in sc_data.csarray; lcon = fill(-Inf, size(sc_data.csarray))) + c103_pr = constraint(core, p_jt_nsc_pr[pr.j_pr, pr.t] + p_jt_rru_off_pr[pr.j_pr, pr.t] - pr.p_rru_off_max*(1-pr.u_on) for pr in sc_data.prarray; lcon = fill(-Inf, size(sc_data.prarray))) + #due to c107, 108, c103_cs is not needed (assume unit commitment satisfies it) + c104_pr = constraint(core, p_jt_rgd_pr[pr.j_pr, pr.t] + p_jt_rrd_on_pr[pr.j_pr, pr.t] - pr.p_rrd_on_max*pr.u_on for pr in sc_data.prarray; lcon = fill(-Inf, size(sc_data.prarray))) + c104_cs = constraint(core, p_jt_rgd_cs[cs.j_cs, cs.t] + p_jt_rrd_on_cs[cs.j_cs, cs.t] - cs.p_rrd_on_max*cs.u_on for cs in sc_data.csarray; lcon = fill(-Inf, size(sc_data.csarray))) + #due to c106, c105_pr is not needed (assume unit commitment satisfies it) + c105_cs = constraint(core, p_jt_rrd_off_cs[cs.j_cs, cs.t] - cs.p_rrd_off_max*(1-cs.u_on) for cs in sc_data.csarray; lcon = fill(-Inf, size(sc_data.csarray))) + #c106 forces p_jt_rrd_off_pr = 0. This is realized in changes to c45, c95, and c105 + #c107 forces p_jt_nsc_cs = 0. This is realized in changes to c43, c93, c101, and c103 + #c108 forces p_jt_rru_off_cs = 0. This is realized in changes to c44, c94, and c103 + + #4.6.7 Relative reserve limits, based on headroom to max/min, producing devices + c109 = constraint(core, p_jt_on_pr[pr.j_pr, pr.t] + p_jt_rgu_pr[pr.j_pr, pr.t] + p_jt_scr_pr[pr.j_pr, pr.t] + p_jt_rru_on_pr[pr.j_pr, pr.t] - pr.p_max*pr.u_on for pr in sc_data.prarray; lcon = fill(-Inf, size(sc_data.prarray))) + c110 = constraint(core, p_jt_on_pr[pr.j_pr, pr.t] - p_jt_rgd_pr[pr.j_pr, pr.t] - p_jt_rrd_on_pr[pr.j_pr, pr.t] - pr.p_min*pr.u_on for pr in sc_data.prarray; ucon = fill(Inf, size(sc_data.prarray))) + c111 = constraint(core, p_jt_su_pr[pr.j_pr, pr.t] + p_jt_sd_pr[pr.j_pr, pr.t] + p_jt_nsc_pr[pr.j_pr, pr.t] + p_jt_rru_off_pr[pr.j_pr, pr.t] - pr.p_max*(1-pr.u_on) for pr in sc_data.prarray; lcon = fill(-Inf, size(sc_data.prarray))) + c112 = constraint(core, q_jt_pr[pr.j_pr, pr.t] + q_jt_qru_pr[pr.j_pr, pr.t] - pr.q_max*(pr.u_on + pr.sum2_T_supc_pr_jt + pr.sum2_T_sdpc_pr_jt) for pr in sc_data.prarray; lcon = fill(-Inf, size(sc_data.prarray))) + c113 = constraint(core, q_jt_pr[pr.j_pr, pr.t] - q_jt_qrd_pr[pr.j_pr, pr.t] - pr.q_min*(pr.u_on + pr.sum2_T_supc_pr_jt + pr.sum2_T_sdpc_pr_jt) for pr in sc_data.prarray; ucon = fill(Inf, size(sc_data.prarray))) + + c114 = constraint(core, q_jt_pr[pr.j_pr, pr.t] + q_jt_qru_pr[pr.j_pr, pr.t] - pr.q_max_p0*(pr.u_on + pr.sum2_T_supc_pr_jt + pr.sum2_T_sdpc_pr_jt) - pr.beta_max*p_jt_pr[pr.j_pr, pr.t] for pr in sc_data.prarray_pqbounds; + lcon = fill(-Inf, size(sc_data.prarray_pqbounds))) + c115 = constraint(core, q_jt_pr[pr.j_pr, pr.t] - q_jt_qrd_pr[pr.j_pr, pr.t] - pr.q_min_p0*(pr.u_on + pr.sum2_T_supc_pr_jt + pr.sum2_T_sdpc_pr_jt) - pr.beta_min*p_jt_pr[pr.j_pr, pr.t] for pr in sc_data.prarray_pqbounds; + ucon = fill(Inf, size(sc_data.prarray_pqbounds))) + c116 = constraint(core, q_jt_pr[pr.j_pr, pr.t] - pr.q_p0*(pr.u_on + pr.sum2_T_supc_pr_jt + pr.sum2_T_sdpc_pr_jt) - pr.beta*p_jt_pr[pr.j_pr, pr.t] for pr in sc_data.prarray_pqe) + #These constraints could be removed and the variables removed to simplify other constraints. However, they are kept for continuity + c117 = constraint(core, q_jt_qru_pr[pr.j_pr, pr.t] for pr in sc_data.prarray_pqe) + c118 = constraint(core, q_jt_qrd_pr[pr.j_pr, pr.t] for pr in sc_data.prarray_pqe) + + #4.6.8 Relative reserve limits, based on headroom to max/min, consuming devices + c119 = constraint(core, p_jt_on_cs[cs.j_cs, cs.t] + p_jt_rgd_cs[cs.j_cs, cs.t] + p_jt_rrd_on_cs[cs.j_cs, cs.t] - cs.p_max*cs.u_on for cs in sc_data.csarray; lcon = fill(-Inf, size(sc_data.csarray))) + c120 = constraint(core, p_jt_on_cs[cs.j_cs, cs.t] - p_jt_rgu_cs[cs.j_cs, cs.t] - p_jt_scr_cs[cs.j_cs, cs.t] - p_jt_rru_on_cs[cs.j_cs, cs.t] -cs.p_min*cs.u_on for cs in sc_data.csarray; ucon = fill(Inf, size(sc_data.csarray))) + c121 = constraint(core, p_jt_su_cs[cs.j_cs, cs.t] + p_jt_sd_cs[cs.j_cs, cs.t] + p_jt_rrd_off_cs[cs.j_cs, cs.t] - cs.p_max*(1-cs.u_on) for cs in sc_data.csarray; lcon = fill(-Inf, size(sc_data.csarray))) + c122 = constraint(core, q_jt_cs[cs.j_cs, cs.t] + q_jt_qrd_cs[cs.j_cs, cs.t] - cs.q_max*(cs.u_on + cs.sum2_T_supc_cs_jt + cs.sum2_T_sdpc_cs_jt) for cs in sc_data.csarray; lcon = fill(-Inf, size(sc_data.csarray))) + c123 = constraint(core, q_jt_cs[cs.j_cs, cs.t] - q_jt_qru_cs[cs.j_cs, cs.t] - cs.q_min*(cs.u_on + cs.sum2_T_supc_cs_jt + cs.sum2_T_sdpc_cs_jt) for cs in sc_data.csarray; ucon = fill(Inf, size(sc_data.csarray))) + c124 = constraint(core, q_jt_cs[cs.j_cs, cs.t] + q_jt_qrd_cs[cs.j_cs, cs.t] - cs.q_max_p0*(cs.u_on + cs.sum2_T_supc_cs_jt + cs.sum2_T_sdpc_cs_jt) - cs.beta_max*p_jt_cs[cs.j_cs, cs.t] for cs in sc_data.csarray_pqbounds; lcon = fill(-Inf, size(sc_data.csarray_pqbounds))) + c125 = constraint(core, q_jt_cs[cs.j_cs, cs.t] - q_jt_qru_cs[cs.j_cs, cs.t] - cs.q_min_p0*(cs.u_on + cs.sum2_T_supc_cs_jt + cs.sum2_T_sdpc_cs_jt) - cs.beta_min*p_jt_cs[cs.j_cs, cs.t] for cs in sc_data.csarray_pqbounds; ucon = fill(Inf, size(sc_data.csarray_pqbounds))) + c126 = constraint(core, q_jt_cs[cs.j_cs, cs.t] - cs.q_p0*(cs.u_on + cs.sum2_T_supc_cs_jt + cs.sum2_T_sdpc_cs_jt) - cs.beta*p_jt_cs[cs.j_cs, cs.t] for cs in sc_data.csarray_pqe) + c127 = constraint(core, q_jt_qru_cs[cs.j_cs, cs.t] for cs in sc_data.csarray_pqe) + c128 = constraint(core, q_jt_qrd_cs[cs.j_cs, cs.t] for cs in sc_data.csarray_pqe) + + #4.6.9 Energy cost and value + c130_pr = constraint(core, p_jt_pr[pr.j_pr, pr.t] for pr in sc_data.prarray) + c130_pr_a = constraint!(core, c130_pr, pr.j_pr + L_J_pr*(pr.t-1) => -p_jtm_pr[pr.flat_k] for pr in sc_data.p_jtm_flattened_pr) + c130_cs = constraint(core, p_jt_cs[cs.j_cs, cs.t] for cs in sc_data.csarray) + c130_cs_a = constraint!(core, c130_cs, cs.j_cs + L_J_cs*(cs.t-1) => -p_jtm_cs[cs.flat_k] for cs in sc_data.p_jtm_flattened_cs) + c131_pr = constraint(core, z_jt_en_pr[pr.j_pr, pr.t]/pr.dt for pr in sc_data.prarray) + c131_pr_a = constraint!(core, c131_pr, pr.j_pr + L_J_pr*(pr.t-1) => -pr.c_en*p_jtm_pr[pr.flat_k] for pr in sc_data.p_jtm_flattened_pr) + c131_cs = constraint(core, z_jt_en_cs[cs.j_cs, cs.t]/cs.dt for cs in sc_data.csarray) + c131_cs_a = constraint!(core, c131_cs, cs.j_cs + L_J_cs*(cs.t-1) => -cs.c_en*p_jtm_cs[cs.flat_k] for cs in sc_data.p_jtm_flattened_cs) + + #4.7 Shunt devices + c132 = constraint(core, p_jt_sh[sh.j_sh, sh.t] - g_jt_sh[sh.j_sh, sh.t]*v_it[sh.bus, sh.t]^2 for sh in sc_data.shuntarray) + c133 = constraint(core, q_jt_sh[sh.j_sh, sh.t] + b_jt_sh[sh.j_sh, sh.t]*v_it[sh.bus, sh.t]^2 for sh in sc_data.shuntarray) + c134 = constraint(core, g_jt_sh[sh.j_sh, sh.t] - sh.g_sh*sh.u_sh for sh in sc_data.shuntarray) + c135 = constraint(core, b_jt_sh[sh.j_sh, sh.t] - sh.b_sh*sh.u_sh for sh in sc_data.shuntarray) + #Assume (136-137) properly handled in uc solution + + #4.8.1 AC branch flow limits and penalties + #AC branches split into ln and xf + c139_ln = constraint(core, z_jt_s_ln[ln.j_ln, ln.t] - ln.dt*sum(sc_data.c_s)*s_jt_plus_ln[ln.j_ln, ln.t] for ln in sc_data.aclbrancharray) + c139_xf = constraint(core, z_jt_s_xf[xf.j_xf, xf.t] -xf.dt*sum(sc_data.c_s)*s_jt_plus_xf[xf.j_xf, xf.t] for xf in sc_data.acxbrancharray) + c140_ln = constraint(core, (p_jt_fr_ln[ln.j_ln, ln.t]^2 + q_jt_fr_ln[ln.j_ln, ln.t]^2)^.5 - ln.s_max - s_jt_plus_ln[ln.j_ln, ln.t] for ln in sc_data.aclbrancharray; + lcon = fill(-Inf, size(sc_data.aclbrancharray))) + c140_xf = constraint(core, (p_jt_fr_xf[xf.j_xf, xf.t]^2 + q_jt_fr_xf[xf.j_xf, xf.t]^2)^.5 - xf.s_max - s_jt_plus_xf[xf.j_xf, xf.t] for xf in sc_data.acxbrancharray; + lcon = fill(-Inf, size(sc_data.acxbrancharray))) + c141_ln = constraint(core, (p_jt_to_ln[ln.j_ln, ln.t]^2 + q_jt_to_ln[ln.j_ln, ln.t]^2)^.5 - ln.s_max - s_jt_plus_ln[ln.j_ln, ln.t] for ln in sc_data.aclbrancharray; + lcon = fill(-Inf, size(sc_data.aclbrancharray))) + c141_xf = constraint(core, (p_jt_to_xf[xf.j_xf, xf.t]^2 + q_jt_to_xf[xf.j_xf, xf.t]^2)^.5 - xf.s_max - s_jt_plus_xf[xf.j_xf, xf.t] for xf in sc_data.acxbrancharray; + lcon = fill(-Inf, size(sc_data.acxbrancharray))) + + #4.8.2 AC branch controls + #c142 forces φ_jt_ln = 0. This is realized in c148-151 and c161 + #c143 forces τ_jt_ln = 1. This is realized in c148-151 + c144 = constraint(core, φ_jt_xf[xf.j_xf, xf.t] - xf.phi_o for xf in sc_data.fpdarray) + c145 = constraint(core, τ_jt_xf[xf.j_xf, xf.t] - xf.tau_o for xf in sc_data.fwrarray) + c146 = constraint(core, φ_jt_xf[xf.j_xf, xf.t] for xf in sc_data.vpdarray; + lcon = [xf.phi_min for xf in sc_data.vpdarray], ucon = [xf.phi_max for xf in sc_data.vpdarray]) + c147 = constraint(core, τ_jt_xf[xf.j_xf, xf.t] for xf in sc_data.vwrarray; + lcon = [xf.tau_min for xf in sc_data.vwrarray], ucon = [xf.tau_max for xf in sc_data.vwrarray]) + + #4.8.3 AC branch flows + + + c148_ln = constraint(core, -p_jt_fr_ln[ln.j_ln, ln.t] + ln.u_on*((ln.g_sr+ln.g_fr)*v_it[ln.fr_bus, ln.t]^2 + (-ln.g_sr*cos(θ_it[ln.fr_bus, ln.t] - θ_it[ln.to_bus, ln.t]) + - ln.b_sr*sin(θ_it[ln.fr_bus, ln.t] - θ_it[ln.to_bus, ln.t]))*v_it[ln.fr_bus, ln.t]*v_it[ln.to_bus, ln.t]) for ln in sc_data.aclbrancharray) + c148_xf = constraint(core, -p_jt_fr_xf[xf.j_xf, xf.t] + xf.u_on*((xf.g_sr+xf.g_fr)*v_it[xf.fr_bus, xf.t]^2/(τ_jt_xf[xf.j_xf, xf.t]^2) + (-xf.g_sr*cos(θ_it[xf.fr_bus, xf.t] - θ_it[xf.to_bus, xf.t] - φ_jt_xf[xf.j_xf, xf.t]) + - xf.b_sr*sin(θ_it[xf.fr_bus, xf.t] - θ_it[xf.to_bus, xf.t] - φ_jt_xf[xf.j_xf, xf.t]))*v_it[xf.fr_bus, xf.t]*v_it[xf.to_bus, xf.t]/τ_jt_xf[xf.j_xf, xf.t]) for xf in sc_data.acxbrancharray) + c149_ln = constraint(core, -q_jt_fr_ln[ln.j_ln, ln.t] + ln.u_on*((-ln.b_sr - ln.b_fr - ln.b_ch/2)*v_it[ln.fr_bus, ln.t]^2 + (ln.b_sr*cos(θ_it[ln.fr_bus, ln.t] - θ_it[ln.to_bus, ln.t]) + - ln.g_sr*sin(θ_it[ln.fr_bus, ln.t] - θ_it[ln.to_bus, ln.t]))*v_it[ln.fr_bus, ln.t]*v_it[ln.to_bus, ln.t]) for ln in sc_data.aclbrancharray) + c149_xf = constraint(core, -q_jt_fr_xf[xf.j_xf, xf.t] + xf.u_on*((-xf.b_sr - xf.b_fr - xf.b_ch/2)*v_it[xf.fr_bus, xf.t]^2/(τ_jt_xf[xf.j_xf, xf.t]^2) + (xf.b_sr*cos(θ_it[xf.fr_bus, xf.t] - θ_it[xf.to_bus, xf.t] - φ_jt_xf[xf.j_xf, xf.t]) + - xf.g_sr*sin(θ_it[xf.fr_bus, xf.t] - θ_it[xf.to_bus, xf.t] - φ_jt_xf[xf.j_xf, xf.t]))*v_it[xf.fr_bus, xf.t]*v_it[xf.to_bus, xf.t]/τ_jt_xf[xf.j_xf, xf.t]) for xf in sc_data.acxbrancharray) + c150_ln = constraint(core, -p_jt_to_ln[ln.j_ln, ln.t] + ln.u_on*((ln.g_sr+ln.g_to)*v_it[ln.to_bus, ln.t]^2 + (-ln.g_sr*cos(θ_it[ln.fr_bus, ln.t] - θ_it[ln.to_bus, ln.t]) + + ln.b_sr*sin(θ_it[ln.fr_bus, ln.t] - θ_it[ln.to_bus, ln.t]))*v_it[ln.fr_bus, ln.t]*v_it[ln.to_bus, ln.t]) for ln in sc_data.aclbrancharray) + c150_xf = constraint(core, -p_jt_to_xf[xf.j_xf, xf.t] + xf.u_on*((xf.g_sr+xf.g_to)*v_it[xf.to_bus, xf.t]^2 + (-xf.g_sr*cos(θ_it[xf.fr_bus, xf.t] - θ_it[xf.to_bus, xf.t] - φ_jt_xf[xf.j_xf, xf.t]) + + xf.b_sr*sin(θ_it[xf.fr_bus, xf.t] - θ_it[xf.to_bus, xf.t] - φ_jt_xf[xf.j_xf, xf.t]))*v_it[xf.fr_bus, xf.t]*v_it[xf.to_bus, xf.t]/τ_jt_xf[xf.j_xf, xf.t]) for xf in sc_data.acxbrancharray) + c151_ln = constraint(core, -q_jt_to_ln[ln.j_ln, ln.t] + ln.u_on*((-ln.b_sr-ln.b_to-ln.b_ch/2)*v_it[ln.to_bus, ln.t]^2 + (ln.b_sr*cos(θ_it[ln.fr_bus, ln.t] - θ_it[ln.to_bus, ln.t]) + + ln.g_sr*sin(θ_it[ln.fr_bus, ln.t] - θ_it[ln.to_bus, ln.t]))*v_it[ln.fr_bus, ln.t]*v_it[ln.to_bus, ln.t]) for ln in sc_data.aclbrancharray) + c151_xf = constraint(core, -q_jt_to_xf[xf.j_xf, xf.t] + xf.u_on*((-xf.b_sr-xf.b_to-xf.b_ch/2)*v_it[xf.to_bus, xf.t]^2 + (xf.b_sr*cos(θ_it[xf.fr_bus, xf.t] - θ_it[xf.to_bus, xf.t] - φ_jt_xf[xf.j_xf, xf.t]) + + xf.g_sr*sin(θ_it[xf.fr_bus, xf.t] - θ_it[xf.to_bus, xf.t] - φ_jt_xf[xf.j_xf, xf.t]))*v_it[xf.fr_bus, xf.t]*v_it[xf.to_bus, xf.t]/τ_jt_xf[xf.j_xf, xf.t]) for xf in sc_data.acxbrancharray) + + + #4.8.4 DC lines + c156 = constraint(core, p_jt_fr_dc[dc.j_dc, dc.t] + p_jt_to_dc[dc.j_dc, dc.t] for dc in sc_data.dclinearray) + + #4.9.1 Penalty on post-contingency AC branch overload + if include_ctg + #ac split into ln and xf + c158_ln = constraint(core, ln.dt*sum(sc_data.c_s)*s_jtk_plus_ln[ln.flat_jtk_ln] - z_jtk_s_ln[ln.flat_jtk_ln] for ln in sc_data.jtk_ln_flattened) + c158_xf = constraint(core, xf.dt*sum(sc_data.c_s)*s_jtk_plus_xf[xf.flat_jtk_xf] - z_jtk_s_xf[xf.flat_jtk_xf] for xf in sc_data.jtk_xf_flattened) + + #4.9.2 Post-contingency AC power flow limits + c159_ln = constraint(core, (p_jtk_ln[ln.flat_jtk_ln]^2 + q_jt_fr_ln[ln.j_ln, ln.t]^2)^0.5 - ln.s_max_ctg - s_jtk_plus_ln[ln.flat_jtk_ln] for ln in sc_data.jtk_ln_flattened; + lcon = fill(-Inf, size(sc_data.jtk_ln_flattened))) + c159_xf = constraint(core, (p_jtk_xf[xf.flat_jtk_xf]^2 + q_jt_fr_xf[xf.j_xf, xf.t]^2)^0.5 - xf.s_max_ctg - s_jtk_plus_xf[xf.flat_jtk_xf] for xf in sc_data.jtk_xf_flattened; + lcon = fill(-Inf, size(sc_data.jtk_xf_flattened))) + c160_ln = constraint(core, (p_jtk_ln[ln.flat_jtk_ln]^2 + q_jt_to_ln[ln.j_ln, ln.t]^2)^0.5 - ln.s_max_ctg - s_jtk_plus_ln[ln.flat_jtk_ln] for ln in sc_data.jtk_ln_flattened; + lcon = fill(-Inf, size(sc_data.jtk_ln_flattened))) + c160_xf = constraint(core, (p_jtk_xf[xf.flat_jtk_xf]^2 + q_jt_to_xf[xf.j_xf, xf.t]^2)^0.5 - xf.s_max_ctg - s_jtk_plus_xf[xf.flat_jtk_xf] for xf in sc_data.jtk_xf_flattened; + lcon = fill(-Inf, size(sc_data.jtk_xf_flattened))) + + #4.9.3 Post-contingency AC branch real power flows + c161_ln = constraint(core, p_jtk_ln[ln.flat_jtk_ln] + ln.b_sr*ln.u_on*(θ_itk[ln.fr_bus, ln.t, ln.ctg] - θ_itk[ln.to_bus, ln.t, ln.ctg]) for ln in sc_data.jtk_ln_flattened) + c161_xf = constraint(core, p_jtk_xf[xf.flat_jtk_xf] + xf.b_sr*xf.u_on*(θ_itk[xf.fr_bus, xf.t, xf.ctg] - θ_itk[xf.to_bus, xf.t, xf.ctg] - φ_jt_xf[xf.j_xf, xf.t]) for xf in sc_data.jtk_xf_flattened) + + #4.9.4 Post-contingency real power balance + + c162 = constraint(core, p_t_sl[t] for t in 1:L_T) + c162_pr = constraint!(core, c162, pr.t => -p_jt_pr[pr.j_pr, pr.t] for pr in sc_data.prarray) + c162_cs = constraint!(core, c162, cs.t => p_jt_cs[cs.j_cs, cs.t] for cs in sc_data.csarray) + c162_sh = constraint!(core, c162, sh.t => p_jt_sh[sh.j_sh, sh.t] for sh in sc_data.shuntarray) + c163 = constraint(core, -p_t_sl[b.t]/I for b in sc_data.k_busarray) #alpha_i = 1/I as per the data format guidelines + c163_ln_fr = constraint!(core, c163, ln.fr_bus + I*(ln.t-1) + I*L_T*(ln.ctg-1) => -p_jtk_ln[ln.flat_jtk_ln] for ln in sc_data.jtk_ln_flattened) + c163_ln_to = constraint!(core, c163, ln.to_bus + I*(ln.t-1) + I*L_T*(ln.ctg-1) => p_jtk_ln[ln.flat_jtk_ln] for ln in sc_data.jtk_ln_flattened) + c163_xf_fr = constraint!(core, c163, xf.fr_bus + I*(xf.t-1) + I*L_T*(xf.ctg-1) => -p_jtk_xf[xf.flat_jtk_xf] for xf in sc_data.jtk_xf_flattened) + c163_xf_to = constraint!(core, c163, xf.to_bus + I*(xf.t-1) + I*L_T*(xf.ctg-1) => p_jtk_xf[xf.flat_jtk_xf] for xf in sc_data.jtk_xf_flattened) + c163_pr = constraint!(core, c163, pr.bus + I*(pr.t-1) + I*L_T*(pr.k-1) => p_jt_pr[pr.j_pr, pr.t] for pr in sc_data.k_prarray) + c163_cs = constraint!(core, c163, cs.bus + I*(cs.t-1) + I*L_T*(cs.k-1) => -p_jt_cs[cs.j_cs, cs.t] for cs in sc_data.k_csarray) + c163_sh = constraint!(core, c163, sh.bus + I*(sh.t-1) + I*L_T*(sh.k-1) => -p_jt_sh[sh.j_sh, sh.t] for sh in sc_data.k_shuntarray) + c163_dc_fr = constraint!(core, c163, dc.fr_bus + I*(dc.t-1) + I*L_T*(dc.ctg-1) => -p_jt_fr_dc[dc.j_dc, dc.t] for dc in sc_data.jtk_dc_flattened) + c163_dc_to = constraint!(core, c163, dc.to_bus + I*(dc.t-1) + I*L_T*(dc.ctg-1) => -p_jt_fr_dc[dc.j_dc, dc.t] for dc in sc_data.jtk_dc_flattened) + end + + + model = ExaModel(core; kwargs...) + + vars = (b_jt_sh = b_jt_sh, + g_jt_sh = g_jt_sh, + #Split e_w_plus into separate sets for W_en_min and W_en_max ad for pr, cs + e_w_plus_min_pr = e_w_plus_min_pr, + e_w_plus_min_cs = e_w_plus_min_cs, + e_w_plus_max_pr = e_w_plus_max_pr, + e_w_plus_max_cs = e_w_plus_max_cs, + p_it = p_it, + p_it_plus = p_it_plus, + #splitting p_jt and q_jt for shunts, producers, and consumers + p_jt_sh = p_jt_sh, + p_jt_pr = p_jt_pr, + p_jt_cs = p_jt_cs, + q_jt_sh = q_jt_sh, + q_jt_pr = q_jt_pr, + q_jt_cs = q_jt_cs, + #Splitting p on, sd , su into pr and cs + p_jt_on_pr = p_jt_on_pr, + p_jt_on_cs = p_jt_on_cs, + p_jt_su_pr = p_jt_su_pr, + p_jt_su_cs = p_jt_su_cs, + p_jt_sd_pr = p_jt_sd_pr, + p_jt_sd_cs = p_jt_sd_cs, + #p_jtm has been flattened and uses only one special index, k_flat + p_jtm_pr = p_jtm_pr, + p_jtm_cs = p_jtm_cs, + #to/from power split into ln, xf, and dc lines + p_jt_fr_ln = p_jt_fr_ln, + p_jt_fr_xf = p_jt_fr_xf, + p_jt_fr_dc = p_jt_fr_dc, + p_jt_to_ln = p_jt_to_ln, + p_jt_to_xf = p_jt_to_xf, + p_jt_to_dc = p_jt_to_dc, + q_jt_fr_ln = q_jt_fr_ln, + q_jt_fr_xf = q_jt_fr_xf, + q_jt_fr_dc = q_jt_fr_dc, + q_jt_to_ln = q_jt_to_ln, + q_jt_to_xf = q_jt_to_xf, + q_jt_to_dc = q_jt_to_dc, + #p_jt rgu, rgd, scr, rru,on, rru,off, rrd,on, rrd,off and q_jt qru/qrd split into pr and cs + p_jt_rgu_pr = p_jt_rgu_pr, + p_jt_rgu_cs = p_jt_rgu_cs, + p_jt_rgd_pr = p_jt_rgd_pr, + p_jt_rgd_cs = p_jt_rgd_cs, + p_jt_scr_pr = p_jt_scr_pr, + p_jt_scr_cs = p_jt_scr_cs, + p_jt_nsc_pr = p_jt_nsc_pr, + p_jt_rru_on_pr = p_jt_rru_on_pr, + p_jt_rru_on_cs = p_jt_rru_on_cs, + p_jt_rru_off_pr = p_jt_rru_off_pr, + p_jt_rrd_on_pr = p_jt_rrd_on_pr, + p_jt_rrd_on_cs = p_jt_rrd_on_cs, + p_jt_rrd_off_cs = p_jt_rrd_off_cs, + q_jt_qru_pr = q_jt_qru_pr, + q_jt_qru_cs = q_jt_qru_cs, + q_jt_qrd_pr = q_jt_qrd_pr, + q_jt_qrd_cs = q_jt_qrd_cs, + p_nt_rgu_req = p_nt_rgu_req, + p_nt_rgd_req = p_nt_rgd_req, + p_nt_scr_req = p_nt_scr_req, + p_nt_nsc_req = p_nt_nsc_req, + p_jt_pr_max = p_jt_pr_max, + p_nt_rgu_plus = p_nt_rgu_plus, + p_nt_rgd_plus = p_nt_rgd_plus, + p_nt_scr_plus = p_nt_scr_plus, + p_nt_nsc_plus = p_nt_nsc_plus, + p_nt_rru_plus = p_nt_rru_plus, + p_nt_rrd_plus = p_nt_rrd_plus, + q_nt_qru_plus = q_nt_qru_plus, + q_nt_qrd_plus = q_nt_qrd_plus, + q_it = q_it, + q_it_plus = q_it_plus, + #s_jt_plus split on ln and xf + s_jt_plus_ln = s_jt_plus_ln, + s_jt_plus_xf = s_jt_plus_xf, + v_it = v_it, + z_w_en_max_pr = z_w_en_max_pr, + z_w_en_max_cs = z_w_en_max_cs, + z_w_en_min_pr = z_w_en_min_pr, + z_w_en_min_cs = z_w_en_min_cs, + #split z_jt_en and on into pr and cs + z_jt_en_pr = z_jt_en_pr, + z_jt_en_cs = z_jt_en_cs, + z_it_p = z_it_p, + z_it_q = z_it_q, + #z_jt_s split into ln and xf + z_jt_s_ln = z_jt_s_ln, + z_jt_s_xf = z_jt_s_xf, + #z_jt rgu, rgd, scr, nsc, rru, rrd, qru, qrd split into pr and cs + z_jt_rgu_pr = z_jt_rgu_pr, + z_jt_rgu_cs = z_jt_rgu_cs, + z_jt_rgd_pr = z_jt_rgd_pr, + z_jt_rgd_cs = z_jt_rgd_cs, + z_jt_scr_pr = z_jt_scr_pr, + z_jt_scr_cs = z_jt_scr_cs, + z_jt_nsc_pr = z_jt_nsc_pr, + z_jt_rru_pr = z_jt_rru_pr, + z_jt_rru_cs = z_jt_rru_cs, + z_jt_rrd_pr = z_jt_rrd_pr, + z_jt_rrd_cs = z_jt_rrd_cs, + z_jt_qru_pr = z_jt_qru_pr, + z_jt_qru_cs = z_jt_qru_cs, + z_jt_qrd_pr = z_jt_qrd_pr, + z_jt_qrd_cs = z_jt_qrd_cs, + z_nt_rgu = z_nt_rgu, + z_nt_rgd = z_nt_rgd, + z_nt_scr = z_nt_scr, + z_nt_nsc = z_nt_nsc, + z_nt_rru = z_nt_rru, + z_nt_rrd = z_nt_rru, + z_nt_qru = z_nt_qru, + z_nt_qrd = z_nt_qrd, + θ_it = θ_it, + #split τjt and φjt into xf only, ln is fixed + τ_jt_xf = τ_jt_xf, + φ_jt_xf = φ_jt_xf) + + if include_ctg + vars = (;vars..., + s_jtk_plus_ln = s_jtk_plus_ln, + s_jtk_plus_xf = s_jtk_plus_xf, + z_jtk_s_ln = z_jtk_s_ln, + z_jtk_s_xf = z_jtk_s_xf, + p_jtk_ln = p_jtk_ln, + p_t_sl = p_t_sl, + p_jtk_xf = p_jtk_xf, + θ_itk = θ_itk, + z_tk_ctg = z_tk_ctg, + z_t_ctg_min = z_t_ctg_min, + z_t_ctg_avg = z_t_ctg_avg) + end + + + return model, sc_data_array, vars, lengths end diff --git a/test/runtests.jl b/test/runtests.jl index a0e66d4..699fa81 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -32,14 +32,14 @@ true_sol_case3_pregen = 29049.351564 true_sol_case5_curve = 78491.04247 true_sol_case5_pregen = 87816.396884 #W storage -true_sol_case3_curve_stor = 25358.827525 -true_sol_case3_curve_stor_func = 25354.331998 -true_sol_case3_pregen_stor = 29023.69118 -true_sol_case3_pregen_stor_func = 29019.172473 +true_sol_case3_curve_stor = 25358.8275 +true_sol_case3_curve_stor_func = 25352.57 +true_sol_case3_pregen_stor = 29023.691 +true_sol_case3_pregen_stor_func = 29019.32 true_sol_case5_curve_stor = 68782.0125 -true_sol_case5_curve_stor_func = 70235.91846 -true_sol_case5_pregen_stor = 79640.08493 -true_sol_case5_pregen_stor_func = 79630.14036 +true_sol_case5_curve_stor_func = 69271.9 +true_sol_case5_pregen_stor = 79640.085 +true_sol_case5_pregen_stor_func = 79630.4 mp_test_cases = [("../data/pglib_opf_case3_lmbd.m", "case3", "../data/case3_5split.Pd", "../data/case3_5split.Qd", true_sol_case3_curve, true_sol_case3_pregen), ("../data/pglib_opf_case5_pjm.m", "case5", "../data/case5_5split.Pd", "../data/case5_5split.Qd", true_sol_case5_curve, true_sol_case5_pregen)] @@ -49,7 +49,7 @@ mp_stor_test_cases = [("../data/pglib_opf_case3_lmbd_mod.m", "case3", "../data/c true_sol_case5_curve_stor, true_sol_case5_curve_stor_func, true_sol_case5_pregen_stor, true_sol_case5_pregen_stor_func)] function example_func(d, srating) - return d + .2/srating*d^2 + return d + 20/srating*d^2 end function runtests()