diff --git a/.github/workflows/TestOnPRs.yml b/.github/workflows/TestOnPRs.yml index 19a1f5a..932e8de 100644 --- a/.github/workflows/TestOnPRs.yml +++ b/.github/workflows/TestOnPRs.yml @@ -18,12 +18,18 @@ concurrency: jobs: test: + strategy: + fail-fast: false + matrix: + version: + - "1.9" + - "1" uses: ./.github/workflows/ReusableTest.yml with: os: ubuntu-latest - version: "1" + version: ${{ matrix.version }} arch: x64 allow_failure: false - run_codecov: true + run_codecov: ${{ matrix.version == '1' }} secrets: codecov_token: ${{ secrets.CODECOV_TOKEN }} diff --git a/CHANGELOG.md b/CHANGELOG.md index f1c2cbb..d0e0d88 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,6 +7,11 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [Unreleased] +### Added + +- Lightweight in-memory DART SCUC, SCED, generator N-1, rolling, and settlement + functionality. + ### Changed - Add `clean_energy_policy = 2` for local annual RPS compliance without diff --git a/Project.toml b/Project.toml index 237a169..1f129ed 100644 --- a/Project.toml +++ b/Project.toml @@ -14,7 +14,6 @@ DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" Dates = "ade2ca70-3891-5945-98fb-dc099432e06a" Distances = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7" GLPK = "60bf3e95-4087-53dc-ae20-288a0d20c6a6" -Gurobi = "2e9cd046-0924-5485-92f1-d5272153d98b" HiGHS = "87dc4568-4c63-4d18-b0c0-bb2238e4078b" JuMP = "4076af6c-e467-56ae-b986-b466b2749572" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" @@ -24,6 +23,7 @@ YAML = "ddb6d928-2868-570f-bddf-ab3f9cf99eb6" [weakdeps] CPLEX = "a076750e-1247-5638-91d2-ce28b192dca0" +Gurobi = "2e9cd046-0924-5485-92f1-d5272153d98b" SCIP = "82193955-e24f-5292-bf16-6f2c5261a85f" [extensions] diff --git a/README.md b/README.md index 788e99b..f49aa92 100644 --- a/README.md +++ b/README.md @@ -19,7 +19,7 @@ The **Holistic Optimization Program for Electricity (HOPE)** model is a transpar 1. `GTEP` mode: a generation & transmission expansion planning model 2. `PCM` mode: a production cost model -3. `DART` mode: a SCUC/SCED market model for simulating day-ahead and real-time markets +3. `DART` module: an individual-generator SCUC/SCED model for day-ahead and real-time markets 4. `OPF` mode: (under development): an optimal power flow model 5. `HOPE-AI` mode: an AI agent helps connect all HOPE modules and enables complex modeling workflows. The current HOPE-AI framework is powered by [PowerAgent](https://github.com/Power-Agent), while more specialized agents are under development. diff --git a/docs/make.jl b/docs/make.jl index 110af6c..c381994 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -34,6 +34,7 @@ pages = OrderedDict( "Notation" => "notation.md", "GTEP" => "GTEP.md", "PCM" => "PCM.md", + "DART" => "DART.md", ], "Input Data Explanation" => ["GTEP Inputs" => "GTEP_inputs.md", "PCM Inputs" => "PCM_inputs.md"], diff --git a/docs/src/95-reference.md b/docs/src/95-reference.md index fa36bc0..c62a02c 100644 --- a/docs/src/95-reference.md +++ b/docs/src/95-reference.md @@ -18,3 +18,9 @@ Pages = ["95-reference.md"] ```@autodocs Modules = [HOPE] ``` + +## DART module + +```@autodocs +Modules = [HOPE.DART] +``` diff --git a/docs/src/DART.md b/docs/src/DART.md new file mode 100644 index 0000000..3240261 --- /dev/null +++ b/docs/src/DART.md @@ -0,0 +1,106 @@ +```@meta +CurrentModule = HOPE.DART +``` + +# DART + +DART is HOPE's lightweight day-ahead and real-time operations module. It is an +in-memory API: callers construct typed Julia inputs and receive typed Julia +results. It does not introduce another case-file schema. + +The V1 model contains: + +- hourly individual-generator day-ahead SCUC; +- rolling real-time SCED, normally at five-minute resolution; +- lossless nodal PTDF transmission constraints; +- regulation, spinning, and quick-start non-spinning reserves; +- strict generator N-1 corrective redispatch by default; +- chronological storage with mutually exclusive charge/discharge operation; and +- day-ahead/real-time energy, reserve, and uplift settlements. + +The implementation is contained in `src/DART.jl`. The main entry points are +`solve_dart_scuc`, `solve_dart_sced`, `run_dart_rolling`, and +`calculate_dart_settlements`. + +## Formulation mapping + +| Formulation group | Implementation | +|:--|:--| +| SCUC-OBJ, UC-1--UC-4 | Commitment, startup/shutdown, minimum up/down time, and commitment costs | +| GEN-1--GEN-7, RAMP-1--RAMP-2 | Generator bounds, transition limits, hourly ramps, and reserve headroom | +| REN-1--REN-2 | Time-varying availability; zero-cost curtailment is eliminated algebraically | +| RES-1--RES-5 | Product capability, response time, requirements, and quick-start eligibility | +| STO-1--STO-6 | Exclusive charge/discharge modes, SOC chronology, and reserve-energy deliverability | +| NET-1--NET-4 | Nodal injection, balance, PTDF flow, and line limits | +| SEC-1--SEC-5 | One corrective scenario per eligible generator | +| RT-1--RT-3 | Fixed DA commitment and energy-only real-time ramping | +| COUP-1 | Binding RT dispatch and SOC are carried into the next solve | + +V1 deliberately uses system-wide reserve requirements, nodal load shedding, +generator contingencies, and the lossless PTDF network option. Line outages, +zonal transport, demand response, and policy constraints remain outside this +small operational core. + +Contingency load shedding is disabled by default. Set +`allow_emergency_contingency_shed = true` in `DARTConfig` only when a penalized +soft-security solve is preferable to infeasibility. + +Quick-start reserve is available only while a unit is offline, in service, +within its response-time qualification, beyond any remaining minimum-down +obligation, and within its availability- and forced-outage-derated capacity. + +Storage charge and discharge are separated by a binary operating-mode variable +in both SCUC and SCED. A horizon-end SOC can be imposed by setting +`terminal_storage_soc_mwh` in `DARTForecast`. It is optional because forcing a +terminal target at the end of every short rolling RT look-ahead can distort the +binding dispatch. During rolling simulation, a supplied target applies only to +a slice that reaches the end of the corresponding input forecast. + +## Operational safeguards + +`DARTDispatchResult` records the solver termination status, primal status, +solve time, and relative optimality gap. An optimal result is always accepted. +By default, a time-limited award is also accepted when it has a feasible primal +solution and a relative gap no greater than one percent. Configure this with +`accept_feasible_time_limit` and `maximum_relative_gap`. The fixed-integer LP +used to produce settlement prices must still solve to optimality. + +Explicit generator-contingency variables scale with the number of generators, +intervals, and monitored outages. DART estimates that block before creating it +and rejects a solve above `maximum_security_variables` (two million by +default). Reduce the look-ahead or the set of generators marked +`contingency_eligible`, or deliberately raise the limit after confirming the +available memory and solver capacity. + +Inputs are checked before model construction for dimensions, finite values, +physical limits, known reserve products, valid state values, and configuration +consistency. + +## Settlements + +Energy settles at fixed-integer nodal prices. Because contingency deliverability +can make reserve value resource-specific, V1 pays reserve awards at each +resource's submitted reserve offer instead of reporting a misleading uniform +security-reserve price. The result field +`reserve_requirement_shadow_price_per_mw_hour` is the informational dual of the +system reserve-requirement constraint; settlements do not use it. The +settlement output identifies this rule as `:pay_as_bid`. + +Settlement input must contain one binding column per result, a valid +chronological RT-to-DA mapping, and exactly one DA interval of RT duration for +each DA award. Reserve and uplift charges are allocated in proportion to +real-time served energy. If there is no served energy, the amounts are reported +explicitly as `unallocated_reserve_charge` and +`unallocated_uplift_charge`; they are included in the settlement balance rather +than silently disappearing. + +## Minimal example + +The ModelCases repository contains `DART_two_bus_example.jl`, a programmatic +example with no generated or committed data files. + +From the HOPE repository root, run it with: + +```bash +julia --project=. ModelCases/DART_two_bus_example.jl +``` diff --git a/docs/src/model_introduction.md b/docs/src/model_introduction.md index 392fb76..e7e4f31 100644 --- a/docs/src/model_introduction.md +++ b/docs/src/model_introduction.md @@ -60,3 +60,9 @@ Both `GTEP` and `PCM` support: - `2`: REG + SPIN + NSPIN - `summary_table`: - `1`: write summary analytics to `output/Analysis/Summary_*.csv` + +## DART Module + +The lightweight in-memory DART API provides individual-generator day-ahead +SCUC, rolling real-time SCED, generator N-1 security, and two-settlement +accounting. See [DART](@ref). diff --git a/src/DART.jl b/src/DART.jl new file mode 100644 index 0000000..0bfd0a1 --- /dev/null +++ b/src/DART.jl @@ -0,0 +1,1824 @@ +""" +Day-ahead security-constrained unit commitment, rolling real-time dispatch, +generator N-1 security, and two-settlement accounting. + +The implementation is intentionally contained in one file. The JuMP +constraints follow the notation groups in `HOPE_SCUC_SCED.pdf` directly. +""" +module DART + +using HiGHS +using JuMP + +const MOI = JuMP.MOI + +# ----------------------------------------------------------------------------- +# Data +# ----------------------------------------------------------------------------- + +Base.@kwdef struct DARTReserveProduct + name::Symbol + direction::Symbol + response_minutes::Int + requires_online::Bool = true +end + +Base.@kwdef struct DARTGenerator + name::String + bus::String + pmax_mw::Float64 + pmin_mw::Float64 = 0.0 + variable_cost_per_mwh::Float64 = 0.0 + no_load_cost_per_hour::Float64 = 0.0 + startup_cost::Float64 = 0.0 + shutdown_cost::Float64 = 0.0 + ramp_up_mw_per_hour::Float64 = Inf + ramp_down_mw_per_hour::Float64 = Inf + real_time_ramp_up_mw_per_hour::Float64 = ramp_up_mw_per_hour + real_time_ramp_down_mw_per_hour::Float64 = ramp_down_mw_per_hour + startup_limit_mw::Float64 = pmax_mw + shutdown_limit_mw::Float64 = pmax_mw + min_up_hours::Int = 0 + min_down_hours::Int = 0 + commitment_required::Bool = true + contingency_eligible::Bool = false + forced_outage_rate::Float64 = 0.0 + quick_start_eligible::Bool = false + quick_start_time_minutes::Int = typemax(Int) + reserve_capability_mw::Dict{Symbol,Float64} = Dict{Symbol,Float64}() + reserve_cost_per_mw_hour::Dict{Symbol,Float64} = Dict{Symbol,Float64}() +end + +Base.@kwdef struct DARTStorage + name::String + bus::String + energy_capacity_mwh::Float64 + charge_capacity_mw::Float64 + discharge_capacity_mw::Float64 + charge_efficiency::Float64 = 1.0 + discharge_efficiency::Float64 = 1.0 + variable_cost_per_mwh::Float64 = 0.0 + reserve_capability_mw::Dict{Symbol,Float64} = Dict{Symbol,Float64}() + reserve_cost_per_mw_hour::Dict{Symbol,Float64} = Dict{Symbol,Float64}() +end + +Base.@kwdef struct DARTNetwork + bus_names::Vector{String} + line_names::Vector{String} = String[] + ptdf::Matrix{Float64} = zeros(0, length(bus_names)) + line_limit_mw::Vector{Float64} = Float64[] + emergency_line_limit_mw::Vector{Float64} = copy(line_limit_mw) +end + +Base.@kwdef struct DARTSystemData + generators::Vector{DARTGenerator} + network::DARTNetwork + storage::Vector{DARTStorage} = DARTStorage[] + reserve_products::Vector{DARTReserveProduct} = default_dart_reserve_products() +end + +""" +Time-varying inputs for one optimization horizon. + +Rows follow generators or buses; columns are chronological intervals. +Interchange is positive for an injection and negative for a withdrawal. +""" +Base.@kwdef struct DARTForecast + interval_hours::Float64 + load_mw::Matrix{Float64} + availability::Matrix{Float64} + generator_in_service::Matrix{Bool} = trues(size(availability)) + interchange_mw::Matrix{Float64} = zeros(size(load_mw)) + reserve_requirement_mw::Dict{Symbol,Vector{Float64}} = Dict{Symbol,Vector{Float64}}() + terminal_storage_soc_mwh::Union{Nothing,Vector{Float64}} = nothing +end + +"""Physical state immediately before an optimization horizon.""" +Base.@kwdef struct DARTState + commitment::Vector{Int} + generation_mw::Vector{Float64} + on_duration_hours::Vector{Float64} + off_duration_hours::Vector{Float64} + storage_soc_mwh::Vector{Float64} = Float64[] + generator_in_service::Vector{Bool} = trues(length(commitment)) +end + +Base.@kwdef struct DARTConfig + value_of_lost_load_per_mwh::Float64 = 100_000.0 + emergency_load_shed_cost_per_mwh::Float64 = 50_000.0 + allow_emergency_contingency_shed::Bool = false + apply_for_derating::Bool = true + security_up_products::Vector{Symbol} = [:reg_up, :spin, :nspin] + security_down_products::Vector{Symbol} = [:reg_down] + accept_feasible_time_limit::Bool = true + maximum_relative_gap::Float64 = 0.01 + maximum_security_variables::Int = 2_000_000 + silent::Bool = true +end + +Base.@kwdef struct DARTDispatchResult + stage::Symbol + objective_value::Float64 + termination_status::MOI.TerminationStatusCode + primal_status::MOI.ResultStatusCode + solve_time_seconds::Float64 + relative_gap::Float64 + interval_hours::Float64 + load_mw::Matrix{Float64} + interchange_mw::Matrix{Float64} + generator_in_service::Matrix{Bool} + generation_mw::Matrix{Float64} + commitment::Matrix{Int} + startup::Matrix{Int} + shutdown::Matrix{Int} + generator_reserve_mw::Dict{Symbol,Matrix{Float64}} + storage_charge_mw::Matrix{Float64} + storage_discharge_mw::Matrix{Float64} + storage_soc_mwh::Matrix{Float64} + storage_reserve_mw::Dict{Symbol,Matrix{Float64}} + load_shed_mw::Matrix{Float64} + line_flow_mw::Matrix{Float64} + lmp_per_mwh::Matrix{Float64} + reserve_requirement_shadow_price_per_mw_hour::Dict{Symbol,Vector{Float64}} + emergency_load_shed_mw::Array{Float64,3} + contingency_generator_indices::Vector{Int} +end + +Base.@kwdef struct DARTRollingResult + day_ahead_results::Vector{DARTDispatchResult} + real_time_results::Vector{DARTDispatchResult} + rt_to_da_index::Vector{Int} + final_state::DARTState +end + +Base.@kwdef struct DARTSettlementResult + reserve_settlement_rule::Symbol + generator_day_ahead_energy::Vector{Float64} + generator_real_time_deviation::Vector{Float64} + generator_reserve_credit::Vector{Float64} + generator_uplift::Vector{Float64} + storage_day_ahead_energy::Vector{Float64} + storage_real_time_deviation::Vector{Float64} + storage_reserve_credit::Vector{Float64} + storage_uplift::Vector{Float64} + load_energy_payment::Vector{Float64} + load_reserve_charge::Vector{Float64} + load_uplift_charge::Vector{Float64} + unallocated_reserve_charge::Float64 + unallocated_uplift_charge::Float64 + interchange_credit::Vector{Float64} + merchandising_surplus::Float64 + settlement_balance::Float64 +end + +function default_dart_reserve_products() + return [ + DARTReserveProduct(name = :reg_up, direction = :up, response_minutes = 5), + DARTReserveProduct(name = :reg_down, direction = :down, response_minutes = 5), + DARTReserveProduct(name = :spin, direction = :up, response_minutes = 10), + DARTReserveProduct( + name = :nspin, + direction = :up, + response_minutes = 30, + requires_online = false, + ), + ] +end + +function default_dart_state(data::DARTSystemData) + commitment = [generator.commitment_required ? 0 : 1 for generator in data.generators] + return DARTState( + commitment = commitment, + generation_mw = zeros(length(data.generators)), + on_duration_hours = [status == 1 ? Inf : 0.0 for status in commitment], + off_duration_hours = [status == 0 ? Inf : 0.0 for status in commitment], + storage_soc_mwh = [0.5 * resource.energy_capacity_mwh for resource in data.storage], + ) +end + +# ----------------------------------------------------------------------------- +# Validation and small helpers +# ----------------------------------------------------------------------------- + +function _bus_indices(data::DARTSystemData) + index = Dict(name => n for (n, name) in pairs(data.network.bus_names)) + generator_bus = [index[generator.bus] for generator in data.generators] + storage_bus = [index[resource.bus] for resource in data.storage] + return generator_bus, storage_bus +end + +_is_nonnegative_or_infinite(value) = !isnan(value) && value >= 0 + +function _validate_inputs( + data::DARTSystemData, + forecast::DARTForecast, + state::DARTState, + config::DARTConfig, +) + G = length(data.generators) + S = length(data.storage) + N = length(data.network.bus_names) + L = length(data.network.line_names) + T = size(forecast.load_mw, 2) + + G > 0 || throw(ArgumentError("DART requires at least one generator.")) + N > 0 || throw(ArgumentError("DART requires at least one bus.")) + T > 0 || throw(ArgumentError("DART requires at least one interval.")) + length(unique(data.network.bus_names)) == N || + throw(ArgumentError("DART bus names must be unique.")) + length(unique(data.network.line_names)) == L || + throw(ArgumentError("DART line names must be unique.")) + length(unique(getfield.(data.generators, :name))) == G || + throw(ArgumentError("DART generator names must be unique.")) + length(unique(getfield.(data.storage, :name))) == S || + throw(ArgumentError("DART storage names must be unique.")) + all(generator.bus in data.network.bus_names for generator in data.generators) || + throw(ArgumentError("Every generator must reference a DART bus.")) + all(resource.bus in data.network.bus_names for resource in data.storage) || + throw(ArgumentError("Every storage resource must reference a DART bus.")) + size(forecast.load_mw) == (N, T) || + throw(ArgumentError("Load must be bus count by interval count.")) + size(forecast.interchange_mw) == (N, T) || + throw(ArgumentError("Interchange must be bus count by interval count.")) + size(forecast.availability) == (G, T) || + throw(ArgumentError("Availability must be generator count by interval count.")) + size(forecast.generator_in_service) == (G, T) || + throw(ArgumentError("Service status must be generator count by interval count.")) + size(data.network.ptdf) == (L, N) || + throw(ArgumentError("PTDF dimensions must be line count by bus count.")) + length(data.network.line_limit_mw) == L || + throw(ArgumentError("Every line requires a normal limit.")) + length(data.network.emergency_line_limit_mw) == L || + throw(ArgumentError("Every line requires an emergency limit.")) + all(isfinite, data.network.ptdf) || throw(ArgumentError("PTDF entries must be finite.")) + all(value -> isfinite(value) && value >= 0, data.network.line_limit_mw) || + throw(ArgumentError("Normal line limits must be finite and nonnegative.")) + all(value -> isfinite(value) && value >= 0, data.network.emergency_line_limit_mw) || + throw(ArgumentError("Emergency line limits must be finite and nonnegative.")) + all(data.network.emergency_line_limit_mw .>= data.network.line_limit_mw) || + throw(ArgumentError("Emergency line limits must be at least normal limits.")) + isfinite(forecast.interval_hours) && forecast.interval_hours > 0 || + throw(ArgumentError("Forecast interval length must be finite and positive.")) + all(isfinite, forecast.load_mw) || throw(ArgumentError("Load must be finite.")) + all(isfinite, forecast.interchange_mw) || + throw(ArgumentError("Interchange must be finite.")) + all(forecast.load_mw .>= 0) || throw(ArgumentError("Load cannot be negative.")) + all(value -> 0 <= value <= 1, forecast.availability) || + throw(ArgumentError("Availability must be between zero and one.")) + if forecast.terminal_storage_soc_mwh !== nothing + target = forecast.terminal_storage_soc_mwh + length(target) == S || + throw(ArgumentError("Terminal storage SOC must match the storage count.")) + all( + isfinite(target[s]) && 0 <= target[s] <= data.storage[s].energy_capacity_mwh for + s = 1:S + ) || throw(ArgumentError("Terminal storage SOC is outside its energy bounds.")) + end + + product_names = getfield.(data.reserve_products, :name) + length(unique(product_names)) == length(product_names) || + throw(ArgumentError("Reserve product names must be unique.")) + for product in data.reserve_products + product.direction in (:up, :down) || + throw(ArgumentError("Reserve direction must be :up or :down.")) + product.response_minutes > 0 || + throw(ArgumentError("Reserve response time must be positive.")) + product.requires_online || + product.direction == :up || + throw(ArgumentError("Offline reserve must be an upward product.")) + end + for (name, requirement) in forecast.reserve_requirement_mw + name in product_names || + throw(ArgumentError("Unknown reserve requirement $(name).")) + length(requirement) == T || + throw(ArgumentError("Reserve requirements must match the horizon.")) + all(value -> isfinite(value) && value >= 0, requirement) || + throw(ArgumentError("Reserve requirements must be finite and nonnegative.")) + end + product_by_name = Dict(product.name => product for product in data.reserve_products) + + for generator in data.generators + all( + isfinite, + ( + generator.pmax_mw, + generator.pmin_mw, + generator.variable_cost_per_mwh, + generator.no_load_cost_per_hour, + generator.startup_cost, + generator.shutdown_cost, + generator.startup_limit_mw, + generator.shutdown_limit_mw, + generator.forced_outage_rate, + ), + ) || throw(ArgumentError("Generator inputs must be finite for $(generator.name).")) + 0 <= generator.pmin_mw <= generator.pmax_mw || + throw(ArgumentError("Generator limits are invalid for $(generator.name).")) + 0 <= generator.forced_outage_rate < 1 || + throw(ArgumentError("FOR must be in [0, 1) for $(generator.name).")) + generator.min_up_hours >= 0 && generator.min_down_hours >= 0 || + throw(ArgumentError("Minimum up/down times cannot be negative.")) + 0 <= generator.startup_limit_mw <= generator.pmax_mw || + throw(ArgumentError("Startup limit is invalid for $(generator.name).")) + 0 <= generator.shutdown_limit_mw <= generator.pmax_mw || + throw(ArgumentError("Shutdown limit is invalid for $(generator.name).")) + all( + _is_nonnegative_or_infinite, + ( + generator.ramp_up_mw_per_hour, + generator.ramp_down_mw_per_hour, + generator.real_time_ramp_up_mw_per_hour, + generator.real_time_ramp_down_mw_per_hour, + ), + ) || throw(ArgumentError("Ramp limits are invalid for $(generator.name).")) + generator.quick_start_time_minutes >= 0 || throw( + ArgumentError("Quick-start time cannot be negative for $(generator.name)."), + ) + all(name in product_names for name in keys(generator.reserve_capability_mw)) || + throw(ArgumentError("Generator reserve capability uses an unknown product.")) + all(name in product_names for name in keys(generator.reserve_cost_per_mw_hour)) || + throw(ArgumentError("Generator reserve cost uses an unknown product.")) + all( + value -> isfinite(value) && 0 <= value <= generator.pmax_mw, + values(generator.reserve_capability_mw), + ) || throw(ArgumentError("Reserve capability is invalid for $(generator.name).")) + all(isfinite, values(generator.reserve_cost_per_mw_hour)) || + throw(ArgumentError("Reserve costs must be finite for $(generator.name).")) + end + for resource in data.storage + all( + isfinite, + ( + resource.energy_capacity_mwh, + resource.charge_capacity_mw, + resource.discharge_capacity_mw, + resource.charge_efficiency, + resource.discharge_efficiency, + resource.variable_cost_per_mwh, + ), + ) || throw(ArgumentError("Storage inputs must be finite for $(resource.name).")) + resource.energy_capacity_mwh >= 0 && + resource.charge_capacity_mw >= 0 && + resource.discharge_capacity_mw >= 0 || + throw(ArgumentError("Storage limits cannot be negative.")) + 0 < resource.charge_efficiency <= 1 && 0 < resource.discharge_efficiency <= 1 || + throw(ArgumentError("Storage efficiencies must be in (0, 1].")) + all(name in product_names for name in keys(resource.reserve_capability_mw)) || + throw(ArgumentError("Storage reserve capability uses an unknown product.")) + all(name in product_names for name in keys(resource.reserve_cost_per_mw_hour)) || + throw(ArgumentError("Storage reserve cost uses an unknown product.")) + for (name, capability) in resource.reserve_capability_mw + limit = + product_by_name[name].direction == :up ? resource.discharge_capacity_mw : + resource.charge_capacity_mw + isfinite(capability) && 0 <= capability <= limit || + throw(ArgumentError("Reserve capability is invalid for $(resource.name).")) + end + all(isfinite, values(resource.reserve_cost_per_mw_hour)) || + throw(ArgumentError("Reserve costs must be finite for $(resource.name).")) + end + + length(state.commitment) == G && + length(state.generation_mw) == G && + length(state.on_duration_hours) == G && + length(state.off_duration_hours) == G && + length(state.generator_in_service) == G || + throw(ArgumentError("Generator state vectors have the wrong length.")) + length(state.storage_soc_mwh) == S || + throw(ArgumentError("Storage state vector has the wrong length.")) + all(status in (0, 1) for status in state.commitment) || + throw(ArgumentError("Commitment state must contain only zero or one.")) + all( + isfinite(state.generation_mw[g]) && + -1e-8 <= state.generation_mw[g] <= data.generators[g].pmax_mw + 1e-8 for g = 1:G + ) || throw(ArgumentError("Initial generation is outside its physical bounds.")) + all(_is_nonnegative_or_infinite, state.on_duration_hours) && + all(_is_nonnegative_or_infinite, state.off_duration_hours) || + throw(ArgumentError("Initial on/off durations must be nonnegative.")) + all( + isfinite(state.storage_soc_mwh[s]) && + -1e-8 <= state.storage_soc_mwh[s] <= data.storage[s].energy_capacity_mwh + 1e-8 + for s = 1:S + ) || throw(ArgumentError("Initial storage SOC is outside its energy bounds.")) + isfinite(config.value_of_lost_load_per_mwh) && config.value_of_lost_load_per_mwh > 0 || + throw(ArgumentError("Value of lost load must be finite and positive.")) + isfinite(config.emergency_load_shed_cost_per_mwh) && + config.emergency_load_shed_cost_per_mwh > 0 || + throw(ArgumentError("Emergency load-shed cost must be finite and positive.")) + length(unique(config.security_up_products)) == length(config.security_up_products) || + throw(ArgumentError("Security-up products cannot contain duplicates.")) + length(unique(config.security_down_products)) == + length(config.security_down_products) || + throw(ArgumentError("Security-down products cannot contain duplicates.")) + all(name in product_names for name in config.security_up_products) || + throw(ArgumentError("A security-up product is not defined.")) + all(name in product_names for name in config.security_down_products) || + throw(ArgumentError("A security-down product is not defined.")) + all(product_by_name[name].direction == :up for name in config.security_up_products) || + throw(ArgumentError("Security-up products must provide upward reserve.")) + all( + product_by_name[name].direction == :down for name in config.security_down_products + ) || throw(ArgumentError("Security-down products must provide downward reserve.")) + isfinite(config.maximum_relative_gap) && config.maximum_relative_gap >= 0 || + throw(ArgumentError("Maximum relative gap must be finite and nonnegative.")) + config.maximum_security_variables > 0 || + throw(ArgumentError("Maximum security variables must be positive.")) + return nothing +end + +_for_factor(generator, config) = + config.apply_for_derating ? 1.0 - generator.forced_outage_rate : 1.0 + +function _available_pmax(data, forecast, config, g, t) + generator = data.generators[g] + service = forecast.generator_in_service[g, t] ? 1.0 : 0.0 + return service * + forecast.availability[g, t] * + generator.pmax_mw * + _for_factor(generator, config) +end + +function _effective_commitment(forecast, commitment, g, t) + return forecast.generator_in_service[g, t] ? commitment[g, t] : 0.0 +end + +function _new_model(optimizer, config) + model = Model(optimizer) + config.silent && set_silent(model) + return model +end + +function _model_relative_gap(model) + try + return relative_gap(model) + catch + return NaN + end +end + +function _model_solve_time(model) + try + return solve_time(model) + catch + return NaN + end +end + +function _require_solution(model, stage, config; require_optimal::Bool = false) + termination = termination_status(model) + primal = primal_status(model) + gap = _model_relative_gap(model) + termination == MOI.OPTIMAL && return nothing + + acceptable_time_limit = + !require_optimal && + config.accept_feasible_time_limit && + termination == MOI.TIME_LIMIT && + primal == MOI.FEASIBLE_POINT && + isfinite(gap) && + gap <= config.maximum_relative_gap + acceptable_time_limit && return nothing + + error( + "DART $(stage) solve failed with termination status $(termination), " * + "primal status $(primal), and relative gap $(gap).", + ) +end + +function _has_discrete_variables(model) + return any( + variable -> is_binary(variable) || is_integer(variable), + all_variables(model), + ) +end + +function _security_variable_count(data, forecast) + G = length(data.generators) + N = length(data.network.bus_names) + L = length(data.network.line_names) + T = size(forecast.load_mw, 2) + C = count(generator.contingency_eligible for generator in data.generators) + return Int128(2G + 2N + L) * T * C + T +end + +function _check_security_model_size(data, forecast, config) + count = _security_variable_count(data, forecast) + count <= config.maximum_security_variables || throw( + ArgumentError( + "Generator N-1 security would create approximately $(count) scenario " * + "variables, above maximum_security_variables=$(config.maximum_security_variables). " * + "Reduce the horizon or contingency set, or explicitly raise the safeguard.", + ), + ) + return nothing +end + +# ----------------------------------------------------------------------------- +# Equations shared by SCUC and SCED +# ----------------------------------------------------------------------------- + +function _add_dispatch!( + model, + data, + forecast, + state, + commitment, + config; + real_time::Bool = false, +) + G = length(data.generators) + S = length(data.storage) + N = length(data.network.bus_names) + L = length(data.network.line_names) + K = length(data.reserve_products) + T = size(forecast.load_mw, 2) + dt = forecast.interval_hours + generator_bus, storage_bus = _bus_indices(data) + + @variable(model, generation[1:G, 1:T] >= 0) + @variable(model, generator_reserve[1:G, 1:T, 1:K] >= 0) + @variable(model, storage_charge[1:S, 1:T] >= 0) + @variable(model, storage_discharge[1:S, 1:T] >= 0) + @variable(model, storage_soc[1:S, 1:T] >= 0) + @variable(model, storage_reserve[1:S, 1:T, 1:K] >= 0) + if S > 0 + @variable(model, storage_charging_mode[1:S, 1:T], Bin) + end + @variable(model, load_shed[1:N, 1:T] >= 0) + @variable(model, injection[1:N, 1:T]) + @variable(model, line_flow[1:L, 1:T]) + + online_up = [ + k for k = 1:K if data.reserve_products[k].direction == :up && + data.reserve_products[k].requires_online + ] + online_down = [ + k for k = 1:K if data.reserve_products[k].direction == :down && + data.reserve_products[k].requires_online + ] + offline_up = [ + k for k = 1:K if data.reserve_products[k].direction == :up && + !data.reserve_products[k].requires_online + ] + all_up = [k for k = 1:K if data.reserve_products[k].direction == :up] + all_down = [k for k = 1:K if data.reserve_products[k].direction == :down] + + # GEN-1 through GEN-5 and RES-1, RES-2, RES-5. + for g = 1:G, t = 1:T + generator = data.generators[g] + effective_commitment = _effective_commitment(forecast, commitment, g, t) + pmax = _available_pmax(data, forecast, config, g, t) + pmin = + generator.commitment_required ? + generator.pmin_mw * _for_factor(generator, config) * effective_commitment : 0.0 + ramp_up = + real_time ? generator.real_time_ramp_up_mw_per_hour : + generator.ramp_up_mw_per_hour + ramp_down = + real_time ? generator.real_time_ramp_down_mw_per_hour : + generator.ramp_down_mw_per_hour + + @constraint( + model, + generation[g, t] + + sum((generator_reserve[g, t, k] for k in online_up); init = 0.0) <= + pmax * commitment[g, t] + ) + @constraint( + model, + generation[g, t] - + sum((generator_reserve[g, t, k] for k in online_down); init = 0.0) >= pmin + ) + + for k = 1:K + product = data.reserve_products[k] + reserve = generator_reserve[g, t, k] + capability = get(generator.reserve_capability_mw, product.name, 0.0) + if capability == 0 + fix(reserve, 0.0; force = true) + elseif product.requires_online + @constraint(model, reserve <= capability * effective_commitment) + response_hours = product.response_minutes / 60 + if product.direction == :up && isfinite(ramp_up) + @constraint( + model, + reserve <= + ramp_up * + response_hours * + _for_factor(generator, config) * + effective_commitment + ) + elseif product.direction == :down && isfinite(ramp_down) + @constraint( + model, + reserve <= + ramp_down * + response_hours * + _for_factor(generator, config) * + effective_commitment + ) + end + elseif generator.quick_start_eligible && + product.direction == :up && + generator.quick_start_time_minutes <= product.response_minutes + available_quick_start = + min(capability, _available_pmax(data, forecast, config, g, t)) + @constraint( + model, + reserve <= available_quick_start * (1 - commitment[g, t]) + ) + else + fix(reserve, 0.0; force = true) + end + end + @constraint( + model, + sum((generator_reserve[g, t, k] for k in offline_up); init = 0.0) <= + pmax * (1 - commitment[g, t]) + ) + end + + # STO-1 through STO-6. The mode binary prevents simultaneous operation. + for s = 1:S, t = 1:T + resource = data.storage[s] + @constraint( + model, + storage_charge[s, t] <= + resource.charge_capacity_mw * storage_charging_mode[s, t] + ) + @constraint( + model, + storage_discharge[s, t] <= + resource.discharge_capacity_mw * (1 - storage_charging_mode[s, t]) + ) + @constraint(model, storage_soc[s, t] <= resource.energy_capacity_mwh) + previous_soc = t == 1 ? state.storage_soc_mwh[s] : storage_soc[s, t-1] + @constraint( + model, + storage_soc[s, t] == + previous_soc + dt * resource.charge_efficiency * storage_charge[s, t] - + dt / resource.discharge_efficiency * storage_discharge[s, t] + ) + @constraint( + model, + storage_discharge[s, t] + + sum((storage_reserve[s, t, k] for k in all_up); init = 0.0) <= + resource.discharge_capacity_mw + ) + @constraint( + model, + storage_charge[s, t] + + sum((storage_reserve[s, t, k] for k in all_down); init = 0.0) <= + resource.charge_capacity_mw + ) + for k = 1:K + capability = + get(resource.reserve_capability_mw, data.reserve_products[k].name, 0.0) + if capability == 0 + fix(storage_reserve[s, t, k], 0.0; force = true) + else + @constraint(model, storage_reserve[s, t, k] <= capability) + end + end + @constraint( + model, + sum( + ( + data.reserve_products[k].response_minutes / 60 * + storage_reserve[s, t, k] for k in all_up + ); + init = 0.0, + ) <= previous_soc + ) + @constraint( + model, + sum( + ( + data.reserve_products[k].response_minutes / 60 * + storage_reserve[s, t, k] for k in all_down + ); + init = 0.0, + ) <= resource.energy_capacity_mwh - previous_soc + ) + end + if forecast.terminal_storage_soc_mwh !== nothing + @constraint( + model, + storage_terminal_soc[s = 1:S], + storage_soc[s, T] == forecast.terminal_storage_soc_mwh[s] + ) + end + + # NET-1 through NET-4. + generators_at_bus = [[g for g = 1:G if generator_bus[g] == n] for n = 1:N] + storage_at_bus = [[s for s = 1:S if storage_bus[s] == n] for n = 1:N] + @constraint( + model, + load_shed_limit[n = 1:N, t = 1:T], + load_shed[n, t] <= forecast.load_mw[n, t] + ) + @constraint( + model, + injection_definition[n = 1:N, t = 1:T], + injection[n, t] == + sum((generation[g, t] for g in generators_at_bus[n]); init = 0.0) + + sum( + (storage_discharge[s, t] - storage_charge[s, t] for s in storage_at_bus[n]); + init = 0.0, + ) + + forecast.interchange_mw[n, t] - forecast.load_mw[n, t] + load_shed[n, t] + ) + @constraint(model, system_balance[t = 1:T], sum(injection[n, t] for n = 1:N) == 0) + @constraint( + model, + ptdf_flow[l = 1:L, t = 1:T], + line_flow[l, t] == sum(data.network.ptdf[l, n] * injection[n, t] for n = 1:N) + ) + @constraint( + model, + line_upper[l = 1:L, t = 1:T], + line_flow[l, t] <= data.network.line_limit_mw[l] + ) + @constraint( + model, + line_lower[l = 1:L, t = 1:T], + line_flow[l, t] >= -data.network.line_limit_mw[l] + ) + + # RES-3 and RES-4. + @constraint( + model, + reserve_requirement[k = 1:K, t = 1:T], + sum(generator_reserve[g, t, k] for g = 1:G) + + sum((storage_reserve[s, t, k] for s = 1:S); init = 0.0) >= + get(forecast.reserve_requirement_mw, data.reserve_products[k].name, zeros(T))[t] + ) + + @expression( + model, + dispatch_cost, + dt * sum( + data.generators[g].variable_cost_per_mwh * generation[g, t] for g = 1:G, t = 1:T + ) + + dt * sum( + get( + data.generators[g].reserve_cost_per_mw_hour, + data.reserve_products[k].name, + 0.0, + ) * generator_reserve[g, t, k] for g = 1:G, t = 1:T, k = 1:K + ) + + dt * sum( + data.storage[s].variable_cost_per_mwh * + (storage_charge[s, t] + storage_discharge[s, t]) for s = 1:S, t = 1:T + ) + + dt * sum( + get( + data.storage[s].reserve_cost_per_mw_hour, + data.reserve_products[k].name, + 0.0, + ) * storage_reserve[s, t, k] for s = 1:S, t = 1:T, k = 1:K + ) + + dt * config.value_of_lost_load_per_mwh * sum(load_shed[n, t] for n = 1:N, t = 1:T) + ) + return nothing +end + +# SEC-1 through SEC-5: one scenario for each eligible generator. +function _add_generator_security!(model, data, forecast, commitment, config) + G = length(data.generators) + N = length(data.network.bus_names) + L = length(data.network.line_names) + T = size(forecast.load_mw, 2) + products = Dict(product.name => k for (k, product) in pairs(data.reserve_products)) + up_products = [products[name] for name in config.security_up_products] + down_products = [products[name] for name in config.security_down_products] + offline_up = [k for k in up_products if !data.reserve_products[k].requires_online] + contingencies = findall(generator.contingency_eligible for generator in data.generators) + C = length(contingencies) + model[:contingency_generator_indices] = contingencies + C == 0 && return @expression(model, emergency_cost, 0.0) + _check_security_model_size(data, forecast, config) + + generator_bus, _ = _bus_indices(data) + generation = model[:generation] + reserve = model[:generator_reserve] + base_injection = model[:injection] + base_shed = model[:load_shed] + + @variable(model, contingency_up[1:G, 1:T, 1:C] >= 0) + @variable(model, contingency_down[1:G, 1:T, 1:C] >= 0) + @variable(model, emergency_load_shed[1:N, 1:T, 1:C] >= 0) + @variable(model, contingency_injection[1:N, 1:T, 1:C]) + @variable(model, contingency_line_flow[1:L, 1:T, 1:C]) + @variable(model, worst_contingency_shed[1:T] >= 0) + + if !config.allow_emergency_contingency_shed + for emergency_shed in emergency_load_shed + fix(emergency_shed, 0.0; force = true) + end + end + + for c = 1:C + failed = contingencies[c] + for t = 1:T + for g = 1:G + if g == failed + fix(contingency_up[g, t, c], 0.0; force = true) + fix(contingency_down[g, t, c], 0.0; force = true) + else + @constraint( + model, + contingency_up[g, t, c] <= + sum((reserve[g, t, k] for k in up_products); init = 0.0) + ) + @constraint( + model, + contingency_down[g, t, c] <= + sum((reserve[g, t, k] for k in down_products); init = 0.0) + ) + effective_commitment = _effective_commitment(forecast, commitment, g, t) + quick_start_reserve = + sum((reserve[g, t, k] for k in offline_up); init = 0.0) + @constraint( + model, + generation[g, t] + contingency_up[g, t, c] - + contingency_down[g, t, c] <= + _available_pmax(data, forecast, config, g, t) * + effective_commitment + quick_start_reserve + ) + generator = data.generators[g] + pmin = + generator.commitment_required ? + generator.pmin_mw * + _for_factor(generator, config) * + effective_commitment : 0.0 + @constraint( + model, + generation[g, t] + contingency_up[g, t, c] - + contingency_down[g, t, c] >= pmin + ) + end + end + + for n = 1:N + @constraint( + model, + emergency_load_shed[n, t, c] <= + forecast.load_mw[n, t] - base_shed[n, t] + ) + failed_output = generator_bus[failed] == n ? generation[failed, t] : 0.0 + response = sum( + ( + contingency_up[g, t, c] - contingency_down[g, t, c] for + g = 1:G if generator_bus[g] == n + ); + init = 0.0, + ) + @constraint( + model, + contingency_injection[n, t, c] == + base_injection[n, t] - failed_output + + response + + emergency_load_shed[n, t, c] + ) + end + @constraint(model, sum(contingency_injection[n, t, c] for n = 1:N) == 0) + @constraint( + model, + worst_contingency_shed[t] >= sum(emergency_load_shed[n, t, c] for n = 1:N) + ) + for l = 1:L + @constraint( + model, + contingency_line_flow[l, t, c] == sum( + data.network.ptdf[l, n] * contingency_injection[n, t, c] for n = 1:N + ) + ) + @constraint( + model, + -data.network.emergency_line_limit_mw[l] <= + contingency_line_flow[l, t, c] <= + data.network.emergency_line_limit_mw[l] + ) + end + end + end + + return @expression( + model, + emergency_cost, + forecast.interval_hours * + config.emergency_load_shed_cost_per_mwh * + sum(worst_contingency_shed[t] for t = 1:T) + ) +end + +# ----------------------------------------------------------------------------- +# Day-ahead SCUC and real-time SCED +# ----------------------------------------------------------------------------- + +function build_dart_scuc_model( + data::DARTSystemData, + forecast::DARTForecast, + state::DARTState; + config::DARTConfig = DARTConfig(), + optimizer = HiGHS.Optimizer, +) + _validate_inputs(data, forecast, state, config) + isapprox(forecast.interval_hours, 1.0; atol = 1e-9) || + throw(ArgumentError("Day-ahead SCUC requires hourly intervals.")) + G = length(data.generators) + T = size(forecast.load_mw, 2) + dt = forecast.interval_hours + model = _new_model(optimizer, config) + + @variable(model, commitment[1:G, 1:T], Bin) + @variable(model, startup[1:G, 1:T], Bin) + @variable(model, shutdown[1:G, 1:T], Bin) + + for g = 1:G + generator = data.generators[g] + if !generator.commitment_required + for t = 1:T + fix(commitment[g, t], 1.0; force = true) + fix(startup[g, t], 0.0; force = true) + fix(shutdown[g, t], 0.0; force = true) + end + continue + end + for t = 1:T + previous = t == 1 ? state.commitment[g] : commitment[g, t-1] + @constraint( + model, + commitment[g, t] - previous == startup[g, t] - shutdown[g, t] + ) + @constraint(model, startup[g, t] + shutdown[g, t] <= 1) + + first_up = max(1, t - ceil(Int, generator.min_up_hours / dt) + 1) + first_down = max(1, t - ceil(Int, generator.min_down_hours / dt) + 1) + @constraint( + model, + sum(startup[g, tau] for tau = first_up:t) <= commitment[g, t] + ) + @constraint( + model, + sum(shutdown[g, tau] for tau = first_down:t) <= 1 - commitment[g, t] + ) + end + + remaining_up = + state.commitment[g] == 1 ? + max(generator.min_up_hours - state.on_duration_hours[g], 0.0) : 0.0 + remaining_down = + state.commitment[g] == 0 ? + max(generator.min_down_hours - state.off_duration_hours[g], 0.0) : 0.0 + for t = 1:min(T, ceil(Int, remaining_up / dt - 1e-9)) + fix(commitment[g, t], 1.0; force = true) + end + for t = 1:min(T, ceil(Int, remaining_down / dt - 1e-9)) + fix(commitment[g, t], 0.0; force = true) + end + end + + _add_dispatch!(model, data, forecast, state, commitment, config) + generation = model[:generation] + reserve = model[:generator_reserve] + online_up = [ + k for k in eachindex(data.reserve_products) if + data.reserve_products[k].requires_online && + data.reserve_products[k].direction == :up + ] + online_down = [ + k for k in eachindex(data.reserve_products) if + data.reserve_products[k].requires_online && + data.reserve_products[k].direction == :down + ] + + # GEN-6/7 and RAMP-1/2. + for g = 1:G, t = 1:T + generator = data.generators[g] + previous_generation = t == 1 ? state.generation_mw[g] : generation[g, t-1] + previous_commitment = t == 1 ? state.commitment[g] : commitment[g, t-1] + previous_service = + t == 1 ? state.generator_in_service[g] : forecast.generator_in_service[g, t-1] + current_service = forecast.generator_in_service[g, t] + previous_effective_commitment = previous_service ? previous_commitment : 0.0 + current_effective_commitment = current_service ? commitment[g, t] : 0.0 + physical_startup = + current_service && !previous_service ? commitment[g, t] : startup[g, t] + + @constraint( + model, + generation[g, t] <= + generator.startup_limit_mw * physical_startup + + _available_pmax(data, forecast, config, g, t) * + (current_effective_commitment - physical_startup) + ) + if t < T + @constraint( + model, + generation[g, t] <= + generator.shutdown_limit_mw * shutdown[g, t+1] + + _available_pmax(data, forecast, config, g, t) * + (current_effective_commitment - shutdown[g, t+1]) + ) + end + + upward_reserve = sum((reserve[g, t, k] for k in online_up); init = 0.0) + downward_reserve = sum((reserve[g, t, k] for k in online_down); init = 0.0) + if isfinite(generator.ramp_up_mw_per_hour) + @constraint( + model, + generation[g, t] + upward_reserve - previous_generation <= + generator.ramp_up_mw_per_hour * + _for_factor(generator, config) * + dt * + previous_effective_commitment + + generator.startup_limit_mw * physical_startup + ) + end + observed_trip = previous_service && !current_service + if isfinite(generator.ramp_down_mw_per_hour) && !observed_trip + @constraint( + model, + previous_generation - generation[g, t] + downward_reserve <= + generator.ramp_down_mw_per_hour * + _for_factor(generator, config) * + dt * + current_effective_commitment + generator.shutdown_limit_mw * shutdown[g, t] + ) + end + end + + # An offline quick-start unit must first complete its minimum-down obligation. + for g = 1:G + generator = data.generators[g] + minimum_down_intervals = ceil(Int, generator.min_down_hours / dt - 1e-9) + remaining_down = + state.commitment[g] == 0 ? + ceil( + Int, + max(generator.min_down_hours - state.off_duration_hours[g], 0.0) / dt - + 1e-9, + ) : 0 + for (k, product) in pairs(data.reserve_products) + product.requires_online && continue + for t = 1:min(T, remaining_down) + fix(reserve[g, t, k], 0.0; force = true) + end + capability = get(generator.reserve_capability_mw, product.name, 0.0) + capability == 0 && continue + for t = 1:T + first_shutdown = max(1, t - minimum_down_intervals + 1) + @constraint( + model, + reserve[g, t, k] <= + capability * (1 - sum(shutdown[g, tau] for tau = first_shutdown:t)) + ) + end + end + end + + emergency_cost = _add_generator_security!(model, data, forecast, commitment, config) + @expression( + model, + commitment_cost, + dt * sum( + data.generators[g].no_load_cost_per_hour * commitment[g, t] for g = 1:G, t = 1:T + ) + sum( + data.generators[g].startup_cost * startup[g, t] + + data.generators[g].shutdown_cost * shutdown[g, t] for g = 1:G, t = 1:T + ) + ) + @objective(model, Min, model[:dispatch_cost] + commitment_cost + emergency_cost) + return model +end + +function _commitment_matrix(data, commitment, T) + G = length(data.generators) + matrix = + commitment isa AbstractVector ? repeat(reshape(Int.(commitment), G, 1), 1, T) : + Int.(commitment) + size(matrix) == (G, T) || + throw(ArgumentError("Fixed commitment has the wrong dimensions.")) + all(status in (0, 1) for status in matrix) || + throw(ArgumentError("Fixed commitment must contain only zero or one.")) + for g = 1:G + if !data.generators[g].commitment_required + matrix[g, :] .= 1 + end + end + return matrix +end + +function build_dart_sced_model( + data::DARTSystemData, + forecast::DARTForecast, + state::DARTState, + fixed_commitment; + config::DARTConfig = DARTConfig(), + optimizer = HiGHS.Optimizer, +) + _validate_inputs(data, forecast, state, config) + G = length(data.generators) + T = size(forecast.load_mw, 2) + dt = forecast.interval_hours + commitment_values = _commitment_matrix(data, fixed_commitment, T) + model = _new_model(optimizer, config) + @variable(model, 0 <= commitment[1:G, 1:T] <= 1) + for g = 1:G, t = 1:T + fix(commitment[g, t], commitment_values[g, t]; force = true) + end + + _add_dispatch!(model, data, forecast, state, commitment, config; real_time = true) + generation = model[:generation] + reserve = model[:generator_reserve] + + # RT-3: ordinary intervals use energy-only RT ramp limits. Commitment and + # equipment-status transitions use their physical boundary limits. + for g = 1:G, t = 1:T + generator = data.generators[g] + previous_generation = t == 1 ? state.generation_mw[g] : generation[g, t-1] + previous_commitment = + generator.commitment_required ? + (t == 1 ? state.commitment[g] : commitment_values[g, t-1]) : 1 + current_commitment = commitment_values[g, t] + previous_service = + t == 1 ? state.generator_in_service[g] : forecast.generator_in_service[g, t-1] + current_service = forecast.generator_in_service[g, t] + startup = current_commitment > previous_commitment + shutdown = current_commitment < previous_commitment + return_to_service = current_service && !previous_service + observed_trip = previous_service && !current_service + + if startup || return_to_service + @constraint( + model, + generation[g, t] - previous_generation <= generator.startup_limit_mw + ) + elseif isfinite(generator.real_time_ramp_up_mw_per_hour) + @constraint( + model, + generation[g, t] - previous_generation <= + generator.real_time_ramp_up_mw_per_hour * dt + ) + end + if !observed_trip + if shutdown + @constraint( + model, + previous_generation - generation[g, t] <= generator.shutdown_limit_mw + ) + elseif isfinite(generator.real_time_ramp_down_mw_per_hour) + @constraint( + model, + previous_generation - generation[g, t] <= + generator.real_time_ramp_down_mw_per_hour * dt + ) + end + end + end + + # Offline reserve remains unavailable until minimum down time is complete. + for g = 1:G + generator = data.generators[g] + off_hours = state.off_duration_hours[g] + previous_commitment = state.commitment[g] + for t = 1:T + if commitment_values[g, t] == 0 && off_hours + 1e-9 < generator.min_down_hours + for (k, product) in pairs(data.reserve_products) + product.requires_online || fix(reserve[g, t, k], 0.0; force = true) + end + end + off_hours = + commitment_values[g, t] == 0 ? + (previous_commitment == 0 ? off_hours + dt : dt) : 0.0 + previous_commitment = commitment_values[g, t] + end + end + + emergency_cost = _add_generator_security!(model, data, forecast, commitment, config) + @objective(model, Min, model[:dispatch_cost] + emergency_cost) + return model +end + +# ----------------------------------------------------------------------------- +# Solving and prices +# ----------------------------------------------------------------------------- + +function _fix_discrete_decisions!(pricing_model, award_model) + for award_variable in all_variables(award_model) + (is_binary(award_variable) || is_integer(award_variable)) || continue + pricing_variable = variable_by_name(pricing_model, name(award_variable)) + pricing_variable === nothing && + error("Pricing model is missing $(name(award_variable)).") + is_binary(pricing_variable) && unset_binary(pricing_variable) + is_integer(pricing_variable) && unset_integer(pricing_variable) + fix(pricing_variable, round(value(award_variable)); force = true) + end + return nothing +end + +function _result(award_model, pricing_model, data, forecast, stage) + G = length(data.generators) + S = length(data.storage) + N = length(data.network.bus_names) + T = size(forecast.load_mw, 2) + K = length(data.reserve_products) + dt = forecast.interval_hours + generator_reserve = Array(value.(award_model[:generator_reserve])) + storage_reserve = Array(value.(award_model[:storage_reserve])) + contingencies = + get(object_dictionary(award_model), :contingency_generator_indices, Int[]) + emergency = + isempty(contingencies) ? zeros(N, T, 0) : + Array(value.(award_model[:emergency_load_shed])) + + return DARTDispatchResult( + stage = stage, + objective_value = objective_value(award_model), + termination_status = termination_status(award_model), + primal_status = primal_status(award_model), + solve_time_seconds = _model_solve_time(award_model), + relative_gap = _model_relative_gap(award_model), + interval_hours = dt, + load_mw = copy(forecast.load_mw), + interchange_mw = copy(forecast.interchange_mw), + generator_in_service = copy(forecast.generator_in_service), + generation_mw = Array(value.(award_model[:generation])), + commitment = round.(Int, Array(value.(award_model[:commitment]))), + startup = haskey(object_dictionary(award_model), :startup) ? + round.(Int, Array(value.(award_model[:startup]))) : zeros(Int, G, T), + shutdown = haskey(object_dictionary(award_model), :shutdown) ? + round.(Int, Array(value.(award_model[:shutdown]))) : zeros(Int, G, T), + generator_reserve_mw = Dict( + data.reserve_products[k].name => generator_reserve[:, :, k] for k = 1:K + ), + storage_charge_mw = Array(value.(award_model[:storage_charge])), + storage_discharge_mw = Array(value.(award_model[:storage_discharge])), + storage_soc_mwh = Array(value.(award_model[:storage_soc])), + storage_reserve_mw = Dict( + data.reserve_products[k].name => storage_reserve[:, :, k] for k = 1:K + ), + load_shed_mw = Array(value.(award_model[:load_shed])), + line_flow_mw = Array(value.(award_model[:line_flow])), + lmp_per_mwh = [ + -dual(pricing_model[:injection_definition][n, t]) / dt for n = 1:N, t = 1:T + ], + reserve_requirement_shadow_price_per_mw_hour = Dict( + data.reserve_products[k].name => + [dual(pricing_model[:reserve_requirement][k, t]) / dt for t = 1:T] for + k = 1:K + ), + emergency_load_shed_mw = emergency, + contingency_generator_indices = contingencies, + ) +end + +function solve_dart_scuc( + data::DARTSystemData, + forecast::DARTForecast, + state::DARTState; + config::DARTConfig = DARTConfig(), + optimizer = HiGHS.Optimizer, +) + award = + build_dart_scuc_model(data, forecast, state; config = config, optimizer = optimizer) + optimize!(award) + _require_solution(award, :day_ahead, config) + + pricing = + build_dart_scuc_model(data, forecast, state; config = config, optimizer = optimizer) + _fix_discrete_decisions!(pricing, award) + optimize!(pricing) + _require_solution(pricing, :day_ahead_pricing, config; require_optimal = true) + return _result(award, pricing, data, forecast, :day_ahead) +end + +function solve_dart_sced( + data::DARTSystemData, + forecast::DARTForecast, + state::DARTState, + fixed_commitment; + config::DARTConfig = DARTConfig(), + optimizer = HiGHS.Optimizer, +) + award = build_dart_sced_model( + data, + forecast, + state, + fixed_commitment; + config = config, + optimizer = optimizer, + ) + optimize!(award) + _require_solution(award, :real_time, config) + + pricing = award + if _has_discrete_variables(award) + pricing = build_dart_sced_model( + data, + forecast, + state, + fixed_commitment; + config = config, + optimizer = optimizer, + ) + _fix_discrete_decisions!(pricing, award) + optimize!(pricing) + _require_solution(pricing, :real_time_pricing, config; require_optimal = true) + end + return _result(award, pricing, data, forecast, :real_time) +end + +# ----------------------------------------------------------------------------- +# Rolling DA/RT simulation +# ----------------------------------------------------------------------------- + +function _slice_forecast(forecast::DARTForecast, first_interval, last_interval) + columns = first_interval:last_interval + return DARTForecast( + interval_hours = forecast.interval_hours, + load_mw = forecast.load_mw[:, columns], + availability = forecast.availability[:, columns], + generator_in_service = forecast.generator_in_service[:, columns], + interchange_mw = forecast.interchange_mw[:, columns], + reserve_requirement_mw = Dict( + name => requirement[columns] for + (name, requirement) in forecast.reserve_requirement_mw + ), + terminal_storage_soc_mwh = last_interval == size(forecast.load_mw, 2) ? + forecast.terminal_storage_soc_mwh : nothing, + ) +end + +function _binding_result(result::DARTDispatchResult) + return DARTDispatchResult( + stage = result.stage, + objective_value = result.objective_value, + termination_status = result.termination_status, + primal_status = result.primal_status, + solve_time_seconds = result.solve_time_seconds, + relative_gap = result.relative_gap, + interval_hours = result.interval_hours, + load_mw = result.load_mw[:, 1:1], + interchange_mw = result.interchange_mw[:, 1:1], + generator_in_service = result.generator_in_service[:, 1:1], + generation_mw = result.generation_mw[:, 1:1], + commitment = result.commitment[:, 1:1], + startup = result.startup[:, 1:1], + shutdown = result.shutdown[:, 1:1], + generator_reserve_mw = Dict( + name => values[:, 1:1] for (name, values) in result.generator_reserve_mw + ), + storage_charge_mw = result.storage_charge_mw[:, 1:1], + storage_discharge_mw = result.storage_discharge_mw[:, 1:1], + storage_soc_mwh = result.storage_soc_mwh[:, 1:1], + storage_reserve_mw = Dict( + name => values[:, 1:1] for (name, values) in result.storage_reserve_mw + ), + load_shed_mw = result.load_shed_mw[:, 1:1], + line_flow_mw = result.line_flow_mw[:, 1:1], + lmp_per_mwh = result.lmp_per_mwh[:, 1:1], + reserve_requirement_shadow_price_per_mw_hour = Dict( + name => values[1:1] for + (name, values) in result.reserve_requirement_shadow_price_per_mw_hour + ), + emergency_load_shed_mw = result.emergency_load_shed_mw[:, 1:1, :], + contingency_generator_indices = result.contingency_generator_indices, + ) +end + +function advance_dart_state(data, state, result::DARTDispatchResult) + dt = result.interval_hours + commitment = vec(result.commitment[:, 1]) + on_hours = copy(state.on_duration_hours) + off_hours = copy(state.off_duration_hours) + for g in eachindex(commitment) + if commitment[g] == 1 + on_hours[g] = state.commitment[g] == 1 ? on_hours[g] + dt : dt + off_hours[g] = 0.0 + else + off_hours[g] = state.commitment[g] == 0 ? off_hours[g] + dt : dt + on_hours[g] = 0.0 + end + end + + generation = max.(vec(result.generation_mw[:, 1]), 0.0) + storage_soc = vec(result.storage_soc_mwh[:, 1]) + for s in eachindex(data.storage) + storage_soc[s] = clamp(storage_soc[s], 0.0, data.storage[s].energy_capacity_mwh) + end + return DARTState( + commitment = commitment, + generation_mw = generation, + on_duration_hours = on_hours, + off_duration_hours = off_hours, + storage_soc_mwh = storage_soc, + generator_in_service = vec(result.generator_in_service[:, 1]), + ) +end + +""" +Run hourly DA SCUC and rolling RT SCED from in-memory forecasts. + +`day_ahead_forecast` must be hourly. `real_time_forecast` is normally +five-minute data. One DA hour and one RT interval bind at each step. +""" +function run_dart_rolling( + data::DARTSystemData, + day_ahead_forecast::DARTForecast, + real_time_forecast::DARTForecast, + initial_state::DARTState = default_dart_state(data); + hours_to_run::Int = min( + size(day_ahead_forecast.load_mw, 2), + floor( + Int, + size(real_time_forecast.load_mw, 2) * real_time_forecast.interval_hours + 1e-9, + ), + ), + day_ahead_lookahead_hours::Int = 24, + real_time_lookahead_intervals::Int = 12, + config::DARTConfig = DARTConfig(), + optimizer = HiGHS.Optimizer, +) + isfinite(day_ahead_forecast.interval_hours) && + isapprox(day_ahead_forecast.interval_hours, 1.0; atol = 1e-9) || + throw(ArgumentError("Day-ahead rolling input must be hourly.")) + isfinite(real_time_forecast.interval_hours) && real_time_forecast.interval_hours > 0 || + throw(ArgumentError("RT interval length must be finite and positive.")) + intervals_per_hour = round(Int, 1 / real_time_forecast.interval_hours) + isapprox(intervals_per_hour * real_time_forecast.interval_hours, 1.0; atol = 1e-9) || + throw(ArgumentError("RT interval length must divide one hour.")) + hours_to_run > 0 || throw(ArgumentError("hours_to_run must be positive.")) + day_ahead_lookahead_hours > 0 || + throw(ArgumentError("day_ahead_lookahead_hours must be positive.")) + real_time_lookahead_intervals > 0 || + throw(ArgumentError("real_time_lookahead_intervals must be positive.")) + hours_to_run <= size(day_ahead_forecast.load_mw, 2) || + throw(ArgumentError("The DA forecast is shorter than hours_to_run.")) + hours_to_run * intervals_per_hour <= size(real_time_forecast.load_mw, 2) || + throw(ArgumentError("The RT forecast is shorter than hours_to_run.")) + + state = initial_state + day_ahead_results = DARTDispatchResult[] + real_time_results = DARTDispatchResult[] + rt_to_da_index = Int[] + + for hour = 1:hours_to_run + da_last = + min(size(day_ahead_forecast.load_mw, 2), hour + day_ahead_lookahead_hours - 1) + da_forecast = _slice_forecast(day_ahead_forecast, hour, da_last) + da = solve_dart_scuc( + data, + da_forecast, + state; + config = config, + optimizer = optimizer, + ) + push!(day_ahead_results, _binding_result(da)) + + for interval_in_hour = 1:intervals_per_hour + rt_first = (hour - 1) * intervals_per_hour + interval_in_hour + rt_last = min( + size(real_time_forecast.load_mw, 2), + rt_first + real_time_lookahead_intervals - 1, + ) + rt_forecast = _slice_forecast(real_time_forecast, rt_first, rt_last) + rt_commitment = + zeros(Int, length(data.generators), size(rt_forecast.load_mw, 2)) + for t in axes(rt_commitment, 2) + da_interval = 1 + div(interval_in_hour + t - 2, intervals_per_hour) + da_interval = min(da_interval, size(da.commitment, 2)) + rt_commitment[:, t] .= da.commitment[:, da_interval] + end + rt = solve_dart_sced( + data, + rt_forecast, + state, + rt_commitment; + config = config, + optimizer = optimizer, + ) + binding_rt = _binding_result(rt) + push!(real_time_results, binding_rt) + push!(rt_to_da_index, hour) + state = advance_dart_state(data, state, binding_rt) + end + end + + return DARTRollingResult( + day_ahead_results = day_ahead_results, + real_time_results = real_time_results, + rt_to_da_index = rt_to_da_index, + final_state = state, + ) +end + +# ----------------------------------------------------------------------------- +# Two-settlement accounting +# ----------------------------------------------------------------------------- + +function _allocate_to_load(total, served_energy) + sum(served_energy) > 0 || return zeros(length(served_energy)) + return total .* served_energy ./ sum(served_energy) +end + +function _validate_binding_result(data, result, expected_stage) + G = length(data.generators) + S = length(data.storage) + N = length(data.network.bus_names) + L = length(data.network.line_names) + product_names = Set(getfield.(data.reserve_products, :name)) + + result.stage == expected_stage || throw( + ArgumentError( + "Settlement result has stage $(result.stage), expected $(expected_stage).", + ), + ) + isfinite(result.interval_hours) && result.interval_hours > 0 || + throw(ArgumentError("Settlement intervals must be finite and positive.")) + isfinite(result.objective_value) || + throw(ArgumentError("Settlement results must have a finite objective value.")) + + expected_dimensions = ( + (:load_mw, result.load_mw, (N, 1)), + (:interchange_mw, result.interchange_mw, (N, 1)), + (:generator_in_service, result.generator_in_service, (G, 1)), + (:generation_mw, result.generation_mw, (G, 1)), + (:commitment, result.commitment, (G, 1)), + (:startup, result.startup, (G, 1)), + (:shutdown, result.shutdown, (G, 1)), + (:storage_charge_mw, result.storage_charge_mw, (S, 1)), + (:storage_discharge_mw, result.storage_discharge_mw, (S, 1)), + (:storage_soc_mwh, result.storage_soc_mwh, (S, 1)), + (:load_shed_mw, result.load_shed_mw, (N, 1)), + (:line_flow_mw, result.line_flow_mw, (L, 1)), + (:lmp_per_mwh, result.lmp_per_mwh, (N, 1)), + ) + for (name, values, dimensions) in expected_dimensions + size(values) == dimensions || throw( + ArgumentError("Settlement field $(name) must have dimensions $(dimensions)."), + ) + all(isfinite, values) || + throw(ArgumentError("Settlement field $(name) must contain finite values.")) + end + + Set(keys(result.generator_reserve_mw)) == product_names || throw( + ArgumentError("Generator reserve settlement products do not match system data."), + ) + Set(keys(result.storage_reserve_mw)) == product_names || throw( + ArgumentError("Storage reserve settlement products do not match system data."), + ) + Set(keys(result.reserve_requirement_shadow_price_per_mw_hour)) == product_names || + throw(ArgumentError("Reserve shadow-price products do not match system data.")) + for name in product_names + generator_reserve = result.generator_reserve_mw[name] + storage_reserve = result.storage_reserve_mw[name] + shadow_price = result.reserve_requirement_shadow_price_per_mw_hour[name] + size(generator_reserve) == (G, 1) || throw( + ArgumentError( + "Generator reserve result $(name) must be generator count by one.", + ), + ) + size(storage_reserve) == (S, 1) || throw( + ArgumentError("Storage reserve result $(name) must be storage count by one."), + ) + length(shadow_price) == 1 || + throw(ArgumentError("Reserve shadow price $(name) must contain one interval.")) + all(isfinite, generator_reserve) && + all(isfinite, storage_reserve) && + all(isfinite, shadow_price) || + throw(ArgumentError("Reserve settlement values must be finite.")) + end + + contingencies = result.contingency_generator_indices + length(unique(contingencies)) == length(contingencies) && + all(index in 1:G for index in contingencies) || + throw(ArgumentError("Contingency generator indices are invalid.")) + size(result.emergency_load_shed_mw) == (N, 1, length(contingencies)) || + throw(ArgumentError("Emergency load-shed settlement dimensions are invalid.")) + all(isfinite, result.emergency_load_shed_mw) || + throw(ArgumentError("Emergency load-shed values must be finite.")) + return nothing +end + +function _validate_settlement_inputs(data, rolling) + da_results = rolling.day_ahead_results + rt_results = rolling.real_time_results + isempty(da_results) && + throw(ArgumentError("Settlement requires at least one day-ahead result.")) + isempty(rt_results) && + throw(ArgumentError("Settlement requires at least one real-time result.")) + length(rolling.rt_to_da_index) == length(rt_results) || + throw(ArgumentError("Every RT result must map to a DA result.")) + issorted(rolling.rt_to_da_index) || + throw(ArgumentError("RT-to-DA settlement mappings must be chronological.")) + all(index in eachindex(da_results) for index in rolling.rt_to_da_index) || + throw(ArgumentError("An RT result maps to a nonexistent DA result.")) + + for result in da_results + _validate_binding_result(data, result, :day_ahead) + isapprox(result.interval_hours, 1.0; atol = 1e-9) || + throw(ArgumentError("Day-ahead settlement results must be hourly.")) + end + for result in rt_results + _validate_binding_result(data, result, :real_time) + end + + covered_hours = zeros(length(da_results)) + for (rt_index, da_index) in pairs(rolling.rt_to_da_index) + covered_hours[da_index] += rt_results[rt_index].interval_hours + end + for da_index in eachindex(da_results) + isapprox( + covered_hours[da_index], + da_results[da_index].interval_hours; + atol = 1e-8, + rtol = 1e-8, + ) || throw( + ArgumentError( + "RT results mapped to DA interval $(da_index) cover " * + "$(covered_hours[da_index]) hours instead of " * + "$(da_results[da_index].interval_hours).", + ), + ) + end + return nothing +end + +function calculate_dart_settlements(data::DARTSystemData, rolling::DARTRollingResult) + da_results = rolling.day_ahead_results + rt_results = rolling.real_time_results + _validate_settlement_inputs(data, rolling) + G = length(data.generators) + S = length(data.storage) + N = length(data.network.bus_names) + generator_bus, storage_bus = _bus_indices(data) + + generator_da = zeros(G) + generator_rt = zeros(G) + generator_reserve = zeros(G) + generator_cost = zeros(G) + storage_da = zeros(S) + storage_rt = zeros(S) + storage_reserve = zeros(S) + storage_cost = zeros(S) + load_energy = zeros(N) + interchange_credit = zeros(N) + served_energy = zeros(N) + + # Day-ahead energy settles at nodal prices. V1 reserve awards settle at + # resource offers because security deliverability can make their dual value + # resource-specific. + for da in da_results + dt = da.interval_hours + for g = 1:G + generator = data.generators[g] + generator_da[g] += + da.generation_mw[g, 1] * da.lmp_per_mwh[generator_bus[g], 1] * dt + generator_cost[g] += + generator.no_load_cost_per_hour * da.commitment[g, 1] * dt + + generator.startup_cost * da.startup[g, 1] + + generator.shutdown_cost * da.shutdown[g, 1] + for product in data.reserve_products + offer = get(generator.reserve_cost_per_mw_hour, product.name, 0.0) + generator_reserve[g] += + da.generator_reserve_mw[product.name][g, 1] * offer * dt + end + end + for s = 1:S + net = da.storage_discharge_mw[s, 1] - da.storage_charge_mw[s, 1] + storage_da[s] += net * da.lmp_per_mwh[storage_bus[s], 1] * dt + for product in data.reserve_products + offer = get(data.storage[s].reserve_cost_per_mw_hour, product.name, 0.0) + storage_reserve[s] += da.storage_reserve_mw[product.name][s, 1] * offer * dt + end + end + for n = 1:N + served = da.load_mw[n, 1] - da.load_shed_mw[n, 1] + load_energy[n] += served * da.lmp_per_mwh[n, 1] * dt + interchange_credit[n] += da.interchange_mw[n, 1] * da.lmp_per_mwh[n, 1] * dt + end + end + + # Real-time energy deviations settle at real-time nodal prices. Reserve + # deviations use the same resource offers as the day-ahead awards. + for (rt_index, rt) in pairs(rt_results) + da = da_results[rolling.rt_to_da_index[rt_index]] + dt = rt.interval_hours + for g = 1:G + generator = data.generators[g] + deviation = rt.generation_mw[g, 1] - da.generation_mw[g, 1] + generator_rt[g] += deviation * rt.lmp_per_mwh[generator_bus[g], 1] * dt + generator_cost[g] += + generator.variable_cost_per_mwh * rt.generation_mw[g, 1] * dt + for product in data.reserve_products + offer = get(generator.reserve_cost_per_mw_hour, product.name, 0.0) + reserve_deviation = + rt.generator_reserve_mw[product.name][g, 1] - + da.generator_reserve_mw[product.name][g, 1] + generator_reserve[g] += reserve_deviation * offer * dt + generator_cost[g] += + offer * rt.generator_reserve_mw[product.name][g, 1] * dt + end + end + for s = 1:S + resource = data.storage[s] + da_net = da.storage_discharge_mw[s, 1] - da.storage_charge_mw[s, 1] + rt_net = rt.storage_discharge_mw[s, 1] - rt.storage_charge_mw[s, 1] + storage_rt[s] += (rt_net - da_net) * rt.lmp_per_mwh[storage_bus[s], 1] * dt + storage_cost[s] += + resource.variable_cost_per_mwh * + (rt.storage_charge_mw[s, 1] + rt.storage_discharge_mw[s, 1]) * + dt + for product in data.reserve_products + reserve_deviation = + rt.storage_reserve_mw[product.name][s, 1] - + da.storage_reserve_mw[product.name][s, 1] + offer = get(resource.reserve_cost_per_mw_hour, product.name, 0.0) + storage_reserve[s] += reserve_deviation * offer * dt + storage_cost[s] += offer * rt.storage_reserve_mw[product.name][s, 1] * dt + end + end + for n = 1:N + da_served = da.load_mw[n, 1] - da.load_shed_mw[n, 1] + rt_served = rt.load_mw[n, 1] - rt.load_shed_mw[n, 1] + load_energy[n] += (rt_served - da_served) * rt.lmp_per_mwh[n, 1] * dt + served_energy[n] += rt_served * dt + interchange_credit[n] += + (rt.interchange_mw[n, 1] - da.interchange_mw[n, 1]) * + rt.lmp_per_mwh[n, 1] * + dt + end + end + + generator_market_credit = generator_da + generator_rt + generator_reserve + storage_market_credit = storage_da + storage_rt + storage_reserve + generator_uplift = max.(generator_cost - generator_market_credit, 0.0) + storage_uplift = max.(storage_cost - storage_market_credit, 0.0) + total_reserve = sum(generator_reserve) + sum(storage_reserve) + total_uplift = sum(generator_uplift) + sum(storage_uplift) + load_reserve = _allocate_to_load(total_reserve, served_energy) + load_uplift = _allocate_to_load(total_uplift, served_energy) + unallocated_reserve = total_reserve - sum(load_reserve) + unallocated_uplift = total_uplift - sum(load_uplift) + + energy_supplier_credit = + sum(generator_da) + + sum(generator_rt) + + sum(storage_da) + + sum(storage_rt) + + sum(interchange_credit) + merchandising_surplus = sum(load_energy) - energy_supplier_credit + settlement_balance = + sum(load_energy) + + sum(load_reserve) + + sum(load_uplift) + + unallocated_reserve + + unallocated_uplift - energy_supplier_credit - total_reserve - total_uplift - + merchandising_surplus + + return DARTSettlementResult( + reserve_settlement_rule = :pay_as_bid, + generator_day_ahead_energy = generator_da, + generator_real_time_deviation = generator_rt, + generator_reserve_credit = generator_reserve, + generator_uplift = generator_uplift, + storage_day_ahead_energy = storage_da, + storage_real_time_deviation = storage_rt, + storage_reserve_credit = storage_reserve, + storage_uplift = storage_uplift, + load_energy_payment = load_energy, + load_reserve_charge = load_reserve, + load_uplift_charge = load_uplift, + unallocated_reserve_charge = unallocated_reserve, + unallocated_uplift_charge = unallocated_uplift, + interchange_credit = interchange_credit, + merchandising_surplus = merchandising_surplus, + settlement_balance = settlement_balance, + ) +end + +export advance_dart_state +export build_dart_sced_model +export build_dart_scuc_model +export calculate_dart_settlements +export DARTConfig +export DARTDispatchResult +export DARTForecast +export DARTGenerator +export DARTNetwork +export DARTReserveProduct +export DARTRollingResult +export DARTSettlementResult +export DARTState +export DARTStorage +export DARTSystemData +export default_dart_reserve_products +export default_dart_state +export run_dart_rolling +export solve_dart_sced +export solve_dart_scuc + +end diff --git a/src/HOPE.jl b/src/HOPE.jl index 4a4b7c8..b84ca17 100644 --- a/src/HOPE.jl +++ b/src/HOPE.jl @@ -19,9 +19,17 @@ using Statistics # extensions so they remain optional for users who do not install them. using Cbc using HiGHS -using Clp using GLPK +# Clp.jl cannot initialize on Apple Silicon. Keep the rest of HOPE loadable on +# that platform; requesting Clp still produces a clear error in initiate_solver. +@static if Sys.isapple() && Sys.ARCH == :aarch64 + const CLP_AVAILABLE = false +else + using Clp + const CLP_AVAILABLE = true +end + #include HOPE module scripts include("constants.jl") #shared constants and configuration include("utils.jl") #utility functions @@ -33,6 +41,8 @@ include("GTEP.jl")#capacity expansion model include("PCM.jl")#production cost model include("write_output.jl")#write output module include("solver_config.jl")#setting solver parameters +include("DART.jl")#day-ahead and real-time operations +using .DART include("solve.jl") #solve model function include("run.jl") #run module include("erec.jl") #EREC postprocessing module @@ -44,8 +54,26 @@ export aggregate_gendata_gtep export aggregate_gendata_pcm export create_GTEP_model export create_PCM_model +export advance_dart_state +export build_dart_sced_model +export build_dart_scuc_model +export calculate_dart_settlements export calculate_erec export calculate_erec_from_output +export DART +export DARTConfig +export DARTDispatchResult +export DARTForecast +export DARTGenerator +export DARTNetwork +export DARTReserveProduct +export DARTRollingResult +export DARTSettlementResult +export DARTState +export DARTStorage +export DARTSystemData +export default_dart_reserve_products +export default_dart_state export default_aggregation_settings export default_erec_settings export default_rep_day_settings @@ -57,10 +85,13 @@ export load_erec_settings export load_rep_day_settings export load_postprocess_snapshot export resolve_rep_day_time_periods +export run_dart_rolling export run_hope export run_hope_holistic export run_hope_holistic_fresh export solve_model +export solve_dart_sced +export solve_dart_scuc export write_output export marginal_load_price_from_dual diff --git a/src/rep_day.jl b/src/rep_day.jl index e33c394..6620e46 100644 --- a/src/rep_day.jl +++ b/src/rep_day.jl @@ -1167,7 +1167,15 @@ function build_storage_rep_linkage(day_assignments::DataFrame) for (prev, cnt) in incoming weight = total_incoming > 0 ? cnt / total_incoming : 0.0 predecessor_weight[(prev, rep)] = weight - push!(transition_rows, (prev, rep, cnt, weight)) + push!( + transition_rows, + ( + PredecessorRepresentativePeriod = prev, + RepresentativePeriod = rep, + Count = cnt, + Weight = weight, + ), + ) end end @@ -1192,7 +1200,15 @@ function build_storage_rep_linkage(day_assignments::DataFrame) lengths = get(run_lengths, rep, Int[]) avg_len = isempty(lengths) ? 0.0 : mean(lengths) max_len = isempty(lengths) ? 0 : maximum(lengths) - push!(run_rows, (rep, length(lengths), avg_len, max_len)) + push!( + run_rows, + ( + RepresentativePeriod = rep, + NumRuns = length(lengths), + AverageRunLength = avg_len, + MaxRunLength = max_len, + ), + ) end return Dict( diff --git a/src/solver_config.jl b/src/solver_config.jl index f66c935..644a1a3 100644 --- a/src/solver_config.jl +++ b/src/solver_config.jl @@ -75,6 +75,39 @@ instantiate_jump_model(optimizer) = Base.invokelatest(Model, optimizer) instantiate_jump_direct_model(optimizer) = Base.invokelatest(direct_model, optimizer) +""" + initiate_solver(config::AbstractDict) + +Create the solver named by `config["solver"]` without case-specific settings. +""" +function initiate_solver(config::AbstractDict) + haskey(config, "solver") || + throw(ArgumentError("Solver configuration is missing the 'solver' key.")) + solver = lowercase(string(config["solver"])) + if solver == "cbc" + return optimizer_with_attributes(Cbc.Optimizer) + elseif solver == "clp" + if CLP_AVAILABLE + return optimizer_with_attributes(Clp.Optimizer) + end + @warn "Clp is unavailable on Apple Silicon; using HiGHS instead." + return optimizer_with_attributes(HiGHS.Optimizer) + elseif solver == "glpk" + return optimizer_with_attributes(GLPK.Optimizer) + elseif solver == "highs" + return optimizer_with_attributes(HiGHS.Optimizer) + elseif haskey(OPTIONAL_SOLVER_PACKAGES, solver) + _ensure_optional_solver_loaded(solver) + builder = + solver == "gurobi" ? _gurobi_optimizer : + solver == "scip" ? _scip_optimizer : _cplex_optimizer + return _worldsafe_optimizer_with_attributes( + Base.invokelatest(builder, Dict{String,Any}()), + ) + end + throw(ArgumentError("Unknown solver '$(solver)'.")) +end + function initiate_solver(case::AbstractString, solver::AbstractString) solver_settings_path = joinpath(case, "Settings", solver * "_settings.yml") solver_settings = open(solver_settings_path) do io @@ -124,6 +157,8 @@ function initiate_solver(case::AbstractString, solver::AbstractString) ) end if solver == "clp" + CLP_AVAILABLE || + error("solver='clp' is unavailable on Apple Silicon. Use solver='highs'.") # Optional solver parameters ############################################ Myfeasib_Tol = 1e-7 if (haskey(solver_settings, "Feasib_Tol")) diff --git a/test/runtests.jl b/test/runtests.jl index b2ab31b..0563597 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -15,4 +15,5 @@ using Test include("test-holistic.jl") include("test-lmp-sign-regression.jl") include("test-pcm-constraints.jl") + include("test-dart.jl") end diff --git a/test/test-dart.jl b/test/test-dart.jl new file mode 100644 index 0000000..fa8fbcd --- /dev/null +++ b/test/test-dart.jl @@ -0,0 +1,710 @@ +using HOPE +using JuMP +using Test + +single_bus() = DARTNetwork(bus_names = ["bus"]) + +function state(commitment, generation; service = trues(length(commitment)), soc = Float64[]) + return DARTState( + commitment = commitment, + generation_mw = generation, + on_duration_hours = [status == 1 ? 10.0 : 0.0 for status in commitment], + off_duration_hours = [status == 0 ? 10.0 : 0.0 for status in commitment], + storage_soc_mwh = soc, + generator_in_service = service, + ) +end + +@testset "DART generator N-1 and reserve response" begin + generators = [ + DARTGenerator( + name = "energy", + bus = "bus", + pmax_mw = 80.0, + variable_cost_per_mwh = 10.0, + contingency_eligible = true, + ), + DARTGenerator( + name = "response", + bus = "bus", + pmax_mw = 100.0, + variable_cost_per_mwh = 30.0, + ramp_up_mw_per_hour = 600.0, + reserve_capability_mw = Dict(:spin => 100.0), + ), + ] + data = DARTSystemData(generators = generators, network = single_bus()) + forecast = DARTForecast( + interval_hours = 1.0, + load_mw = fill(80.0, 1, 1), + availability = ones(2, 1), + ) + result = solve_dart_scuc(data, forecast, state([1, 1], [80.0, 0.0])) + + @test result.generation_mw[:, 1] ≈ [80.0, 0.0] atol = 1e-6 + @test result.generator_reserve_mw[:spin][2, 1] ≈ 80.0 atol = 1e-6 + @test sum(result.emergency_load_shed_mw) <= 1e-6 + @test result.contingency_generator_indices == [1] + @test result.termination_status == JuMP.MOI.OPTIMAL + @test result.primal_status == JuMP.MOI.FEASIBLE_POINT + @test result.solve_time_seconds >= 0 + @test result.relative_gap <= 1e-8 +end + +@testset "DART strict and soft N-1 security" begin + generators = [ + DARTGenerator( + name = "contingency", + bus = "bus", + pmax_mw = 10.0, + contingency_eligible = true, + ), + DARTGenerator( + name = "expensive_response", + bus = "bus", + pmax_mw = 10.0, + variable_cost_per_mwh = 60_000.0, + ramp_up_mw_per_hour = 600.0, + reserve_capability_mw = Dict(:spin => 10.0), + reserve_cost_per_mw_hour = Dict(:spin => 100_000.0), + ), + ] + data = DARTSystemData(generators = generators, network = single_bus()) + forecast = DARTForecast( + interval_hours = 1.0, + load_mw = fill(10.0, 1, 1), + availability = ones(2, 1), + ) + initial_state = state([1, 1], [10.0, 0.0]) + + strict = solve_dart_scuc(data, forecast, initial_state) + soft = solve_dart_scuc( + data, + forecast, + initial_state; + config = DARTConfig(allow_emergency_contingency_shed = true), + ) + + @test sum(strict.emergency_load_shed_mw) <= 1e-6 + @test strict.generation_mw[2, 1] ≈ 10.0 atol = 1e-6 + @test sum(soft.emergency_load_shed_mw) ≈ 10.0 atol = 1e-6 + @test soft.generation_mw[1, 1] ≈ 10.0 atol = 1e-6 +end + +@testset "DART product response limits and quick start" begin + products = [ + DARTReserveProduct(name = :reg_up, direction = :up, response_minutes = 5), + DARTReserveProduct(name = :spin, direction = :up, response_minutes = 10), + DARTReserveProduct( + name = :nspin, + direction = :up, + response_minutes = 30, + requires_online = false, + ), + ] + generators = [ + DARTGenerator( + name = "outage", + bus = "bus", + pmax_mw = 15.0, + variable_cost_per_mwh = 10.0, + contingency_eligible = true, + ), + DARTGenerator( + name = "online_response", + bus = "bus", + pmax_mw = 100.0, + variable_cost_per_mwh = 30.0, + ramp_up_mw_per_hour = 60.0, + reserve_capability_mw = Dict(:reg_up => 100.0, :spin => 100.0), + ), + ] + data = DARTSystemData( + generators = generators, + network = single_bus(), + reserve_products = products, + ) + forecast = DARTForecast( + interval_hours = 1.0, + load_mw = fill(15.0, 1, 1), + availability = ones(2, 1), + ) + config = DARTConfig( + security_up_products = [:reg_up, :spin, :nspin], + security_down_products = Symbol[], + ) + result = solve_dart_scuc(data, forecast, state([1, 1], [15.0, 0.0]); config = config) + + @test result.generator_reserve_mw[:reg_up][2, 1] ≈ 5.0 atol = 1e-6 + @test result.generator_reserve_mw[:spin][2, 1] ≈ 10.0 atol = 1e-6 + @test sum(result.emergency_load_shed_mw) <= 1e-6 + + quick_start = DARTGenerator( + name = "quick_start", + bus = "bus", + pmax_mw = 20.0, + variable_cost_per_mwh = 50.0, + quick_start_eligible = true, + quick_start_time_minutes = 20, + reserve_capability_mw = Dict(:nspin => 20.0), + ) + quick_data = DARTSystemData( + generators = [generators[1], quick_start], + network = single_bus(), + reserve_products = products, + ) + quick_result = + solve_dart_scuc(quick_data, forecast, state([1, 0], [15.0, 0.0]); config = config) + @test quick_result.commitment[2, 1] == 0 + @test quick_result.generator_reserve_mw[:nspin][2, 1] >= 15.0 - 1e-6 + @test quick_result.generator_reserve_mw[:nspin][2, 1] <= 20.0 + 1e-6 + @test sum(quick_result.emergency_load_shed_mw) <= 1e-6 + + unavailable_forecast = DARTForecast( + interval_hours = 1.0, + load_mw = fill(15.0, 1, 1), + availability = reshape([1.0, 0.0], 2, 1), + ) + unavailable = solve_dart_scuc( + quick_data, + unavailable_forecast, + state([1, 0], [15.0, 0.0]); + config = DARTConfig( + allow_emergency_contingency_shed = true, + security_up_products = [:reg_up, :spin, :nspin], + security_down_products = Symbol[], + ), + ) + @test unavailable.generator_reserve_mw[:nspin][2, 1] <= 1e-8 + @test sum(unavailable.emergency_load_shed_mw) ≈ 15.0 atol = 1e-6 + + minimum_down_generator = DARTGenerator( + name = "minimum_down_quick_start", + bus = "bus", + pmax_mw = 10.0, + min_down_hours = 2, + quick_start_eligible = true, + quick_start_time_minutes = 20, + reserve_capability_mw = Dict(:nspin => 10.0), + ) + minimum_down_data = DARTSystemData( + generators = [minimum_down_generator], + network = single_bus(), + reserve_products = products, + ) + minimum_down_forecast = DARTForecast( + interval_hours = 1.0, + load_mw = zeros(1, 2), + availability = ones(1, 2), + reserve_requirement_mw = Dict(:nspin => [0.0, 10.0]), + ) + @test_throws ErrorException solve_dart_scuc( + minimum_down_data, + minimum_down_forecast, + state([1], [0.0]); + config = config, + ) + + offline_products = [ + DARTReserveProduct( + name = :nspin_10, + direction = :up, + response_minutes = 10, + requires_online = false, + ), + DARTReserveProduct( + name = :nspin_30, + direction = :up, + response_minutes = 30, + requires_online = false, + ), + ] + aggregate_quick_start = DARTGenerator( + name = "aggregate_quick_start", + bus = "bus", + pmax_mw = 10.0, + quick_start_eligible = true, + quick_start_time_minutes = 10, + reserve_capability_mw = Dict(:nspin_10 => 10.0, :nspin_30 => 10.0), + ) + aggregate_data = DARTSystemData( + generators = [aggregate_quick_start], + network = single_bus(), + reserve_products = offline_products, + ) + aggregate_forecast = DARTForecast( + interval_hours = 1.0, + load_mw = zeros(1, 1), + availability = ones(1, 1), + reserve_requirement_mw = Dict(:nspin_10 => [10.0], :nspin_30 => [10.0]), + ) + @test_throws ErrorException solve_dart_scuc( + aggregate_data, + aggregate_forecast, + state([0], [0.0]); + config = DARTConfig( + security_up_products = Symbol[], + security_down_products = Symbol[], + ), + ) +end + +@testset "DART real-time ramp and equipment transitions" begin + generator = DARTGenerator( + name = "unit", + bus = "bus", + pmax_mw = 100.0, + startup_limit_mw = 5.0, + shutdown_limit_mw = 30.0, + real_time_ramp_up_mw_per_hour = 60.0, + real_time_ramp_down_mw_per_hour = 60.0, + ) + data = DARTSystemData(generators = [generator], network = single_bus()) + config = DARTConfig( + security_up_products = [:reg_up, :spin, :nspin], + security_down_products = [:reg_down], + ) + + forecast = DARTForecast( + interval_hours = 1 / 12, + load_mw = fill(15.0, 1, 1), + availability = ones(1, 1), + ) + ramp = solve_dart_sced(data, forecast, state([1], [10.0]), [1]; config = config) + @test ramp.generation_mw[1, 1] ≈ 15.0 atol = 1e-6 + + startup = solve_dart_sced(data, forecast, state([0], [0.0]), [1]; config = config) + @test startup.generation_mw[1, 1] ≈ 5.0 atol = 1e-6 + @test startup.load_shed_mw[1, 1] ≈ 10.0 atol = 1e-6 + + trip_forecast = DARTForecast( + interval_hours = 1 / 12, + load_mw = zeros(1, 1), + availability = ones(1, 1), + generator_in_service = falses(1, 1), + ) + trip = solve_dart_sced( + data, + trip_forecast, + state([1], [30.0]; service = [true]), + [1]; + config = config, + ) + @test trip.generation_mw[1, 1] <= 1e-8 + + return_to_service = solve_dart_sced( + data, + forecast, + state([1], [0.0]; service = [false]), + [1]; + config = config, + ) + @test return_to_service.generation_mw[1, 1] ≈ 5.0 atol = 1e-6 + + noncommit_data = DARTSystemData( + generators = [ + DARTGenerator( + name = "noncommit", + bus = "bus", + pmax_mw = 100.0, + commitment_required = false, + ), + ], + network = single_bus(), + ) + noncommit_result = solve_dart_sced( + noncommit_data, + forecast, + default_dart_state(noncommit_data), + [0]; + config = config, + ) + @test noncommit_result.commitment[1, 1] == 1 + @test noncommit_result.generation_mw[1, 1] ≈ 15.0 atol = 1e-6 + @test noncommit_result.load_shed_mw[1, 1] <= 1e-6 +end + +@testset "DART nodal PTDF and storage chronology" begin + network = DARTNetwork( + bus_names = ["one", "two"], + line_names = ["one_to_two"], + ptdf = reshape([0.0, -1.0], 1, 2), + line_limit_mw = [20.0], + emergency_line_limit_mw = [25.0], + ) + generators = [ + DARTGenerator( + name = "cheap", + bus = "one", + pmax_mw = 100.0, + variable_cost_per_mwh = 10.0, + commitment_required = false, + ), + DARTGenerator( + name = "local", + bus = "two", + pmax_mw = 100.0, + variable_cost_per_mwh = 30.0, + commitment_required = false, + ), + ] + data = DARTSystemData(generators = generators, network = network) + forecast = DARTForecast( + interval_hours = 1.0, + load_mw = reshape([0.0, 40.0], 2, 1), + availability = ones(2, 1), + ) + result = solve_dart_scuc(data, forecast, default_dart_state(data)) + @test result.generation_mw[:, 1] ≈ [20.0, 20.0] atol = 1e-6 + @test result.line_flow_mw[1, 1] ≈ 20.0 atol = 1e-6 + @test result.lmp_per_mwh[:, 1] ≈ [10.0, 30.0] atol = 1e-6 + + storage = DARTStorage( + name = "battery", + bus = "bus", + energy_capacity_mwh = 10.0, + charge_capacity_mw = 10.0, + discharge_capacity_mw = 10.0, + ) + storage_generator = DARTGenerator( + name = "expensive", + bus = "bus", + pmax_mw = 100.0, + variable_cost_per_mwh = 100.0, + commitment_required = false, + ) + storage_data = DARTSystemData( + generators = [storage_generator], + storage = [storage], + network = single_bus(), + ) + storage_forecast = DARTForecast( + interval_hours = 1.0, + load_mw = fill(10.0, 1, 1), + availability = ones(1, 1), + ) + storage_result = + solve_dart_scuc(storage_data, storage_forecast, state([1], [0.0]; soc = [10.0])) + @test storage_result.storage_discharge_mw[1, 1] ≈ 10.0 atol = 1e-6 + @test storage_result.storage_soc_mwh[1, 1] <= 1e-6 + @test storage_result.storage_charge_mw[1, 1] * + storage_result.storage_discharge_mw[1, 1] <= 1e-8 + + storage_rt = solve_dart_sced( + storage_data, + storage_forecast, + state([1], [0.0]; soc = [10.0]), + [1], + ) + @test storage_rt.storage_discharge_mw[1, 1] ≈ 10.0 atol = 1e-6 + @test isfinite(storage_rt.lmp_per_mwh[1, 1]) + + terminal_forecast = DARTForecast( + interval_hours = 1.0, + load_mw = reshape([0.0, 8.0], 1, 2), + availability = reshape([1.0, 0.0], 1, 2), + terminal_storage_soc_mwh = [0.0], + ) + terminal_result = + solve_dart_scuc(storage_data, terminal_forecast, state([1], [0.0]; soc = [0.0])) + @test terminal_result.storage_charge_mw[1, 1] ≈ 8.0 atol = 1e-6 + @test terminal_result.storage_discharge_mw[1, 2] ≈ 8.0 atol = 1e-6 + @test terminal_result.storage_soc_mwh[1, 2] ≈ 0.0 atol = 1e-6 + + must_run = DARTGenerator(name = "must_run", bus = "bus", pmax_mw = 20.0, pmin_mw = 20.0) + inefficient_storage = DARTStorage( + name = "inefficient", + bus = "bus", + energy_capacity_mwh = 10.0, + charge_capacity_mw = 100.0, + discharge_capacity_mw = 100.0, + charge_efficiency = 0.8, + discharge_efficiency = 0.8, + ) + physical_data = DARTSystemData( + generators = [must_run], + storage = [inefficient_storage], + network = single_bus(), + ) + zero_load = + DARTForecast(interval_hours = 1.0, load_mw = zeros(1, 1), availability = ones(1, 1)) + @test_throws ErrorException solve_dart_sced( + physical_data, + zero_load, + state([1], [20.0]; soc = [10.0]), + [1], + ) +end + +@testset "DART rolling simulation and settlements" begin + generators = [ + DARTGenerator( + name = "energy", + bus = "bus", + pmax_mw = 80.0, + variable_cost_per_mwh = 10.0, + contingency_eligible = true, + ), + DARTGenerator( + name = "response", + bus = "bus", + pmax_mw = 100.0, + variable_cost_per_mwh = 30.0, + ramp_up_mw_per_hour = 600.0, + real_time_ramp_up_mw_per_hour = 600.0, + real_time_ramp_down_mw_per_hour = 600.0, + reserve_capability_mw = Dict(:spin => 100.0), + reserve_cost_per_mw_hour = Dict(:spin => 5.0), + ), + ] + data = DARTSystemData(generators = generators, network = single_bus()) + day_ahead = DARTForecast( + interval_hours = 1.0, + load_mw = fill(70.0, 1, 1), + availability = ones(2, 1), + ) + real_time = DARTForecast( + interval_hours = 1 / 12, + load_mw = fill(70.0, 1, 12), + availability = ones(2, 12), + ) + rolling = run_dart_rolling( + data, + day_ahead, + real_time, + state([1, 1], [70.0, 0.0]); + hours_to_run = 1, + day_ahead_lookahead_hours = 1, + real_time_lookahead_intervals = 3, + ) + settlement = calculate_dart_settlements(data, rolling) + + @test length(rolling.day_ahead_results) == 1 + @test length(rolling.real_time_results) == 12 + @test all(result.generation_mw[1, 1] ≈ 70.0 for result in rolling.real_time_results) + @test rolling.final_state.generation_mw ≈ [70.0, 0.0] atol = 1e-6 + @test all( + result.reserve_requirement_shadow_price_per_mw_hour[:spin][1] ≈ 0.0 for + result in rolling.real_time_results + ) + @test settlement.generator_reserve_credit[2] ≈ 350.0 atol = 1e-6 + @test settlement.settlement_balance ≈ 0.0 atol = 1e-6 + + incomplete = DARTRollingResult( + day_ahead_results = rolling.day_ahead_results, + real_time_results = rolling.real_time_results[1:11], + rt_to_da_index = rolling.rt_to_da_index[1:11], + final_state = rolling.final_state, + ) + @test_throws ArgumentError calculate_dart_settlements(data, incomplete) +end + +@testset "DART numerical two-settlement accounting" begin + generator = DARTGenerator( + name = "supplier", + bus = "bus", + pmax_mw = 200.0, + variable_cost_per_mwh = 10.0, + no_load_cost_per_hour = 20.0, + commitment_required = false, + reserve_capability_mw = Dict(:spin => 50.0), + reserve_cost_per_mw_hour = Dict(:spin => 5.0), + ) + data = DARTSystemData(generators = [generator], network = single_bus()) + day_ahead = DARTForecast( + interval_hours = 1.0, + load_mw = fill(100.0, 1, 1), + availability = ones(1, 1), + interchange_mw = fill(10.0, 1, 1), + reserve_requirement_mw = Dict(:spin => [10.0]), + ) + real_time = DARTForecast( + interval_hours = 1.0, + load_mw = fill(110.0, 1, 1), + availability = ones(1, 1), + interchange_mw = fill(10.0, 1, 1), + reserve_requirement_mw = Dict(:spin => [12.0]), + ) + rolling = run_dart_rolling( + data, + day_ahead, + real_time, + state([1], [90.0]); + hours_to_run = 1, + day_ahead_lookahead_hours = 1, + real_time_lookahead_intervals = 1, + ) + settlement = calculate_dart_settlements(data, rolling) + da = only(rolling.day_ahead_results) + rt = only(rolling.real_time_results) + + @test da.generation_mw[1, 1] ≈ 90.0 atol = 1e-6 + @test rt.generation_mw[1, 1] ≈ 100.0 atol = 1e-6 + @test da.generator_reserve_mw[:spin][1, 1] ≈ 10.0 atol = 1e-6 + @test rt.generator_reserve_mw[:spin][1, 1] ≈ 12.0 atol = 1e-6 + @test da.reserve_requirement_shadow_price_per_mw_hour[:spin][1] ≈ 5.0 atol = 1e-6 + @test rt.reserve_requirement_shadow_price_per_mw_hour[:spin][1] ≈ 5.0 atol = 1e-6 + + @test settlement.generator_day_ahead_energy ≈ [900.0] atol = 1e-6 + @test settlement.generator_real_time_deviation ≈ [100.0] atol = 1e-6 + @test settlement.generator_reserve_credit ≈ [60.0] atol = 1e-6 + @test settlement.generator_uplift ≈ [20.0] atol = 1e-6 + @test settlement.load_energy_payment ≈ [1_100.0] atol = 1e-6 + @test settlement.load_reserve_charge ≈ [60.0] atol = 1e-6 + @test settlement.load_uplift_charge ≈ [20.0] atol = 1e-6 + @test settlement.interchange_credit ≈ [100.0] atol = 1e-6 + @test settlement.merchandising_surplus ≈ 0.0 atol = 1e-6 + @test settlement.settlement_balance ≈ 0.0 atol = 1e-6 +end + +@testset "DART zero-load settlement residual" begin + generator = DARTGenerator( + name = "reserve_supplier", + bus = "bus", + pmax_mw = 100.0, + commitment_required = false, + reserve_capability_mw = Dict(:spin => 10.0), + reserve_cost_per_mw_hour = Dict(:spin => 5.0), + ) + data = DARTSystemData(generators = [generator], network = single_bus()) + forecast = DARTForecast( + interval_hours = 1.0, + load_mw = zeros(1, 1), + availability = ones(1, 1), + reserve_requirement_mw = Dict(:spin => [10.0]), + ) + rolling = run_dart_rolling( + data, + forecast, + forecast; + hours_to_run = 1, + day_ahead_lookahead_hours = 1, + real_time_lookahead_intervals = 1, + ) + settlement = calculate_dart_settlements(data, rolling) + + @test settlement.reserve_settlement_rule == :pay_as_bid + @test settlement.generator_reserve_credit ≈ [50.0] atol = 1e-6 + @test settlement.load_reserve_charge ≈ [0.0] atol = 1e-6 + @test settlement.unallocated_reserve_charge ≈ 50.0 atol = 1e-6 + @test settlement.unallocated_uplift_charge ≈ 0.0 atol = 1e-6 + @test settlement.settlement_balance ≈ 0.0 atol = 1e-6 + + invalid_mapping = DARTRollingResult( + day_ahead_results = rolling.day_ahead_results, + real_time_results = rolling.real_time_results, + rt_to_da_index = [2], + final_state = rolling.final_state, + ) + @test_throws ArgumentError calculate_dart_settlements(data, invalid_mapping) +end + +@testset "DART 24-hour generator N-1 case and size safeguard" begin + generators = [ + DARTGenerator( + name = "unit_$(g)", + bus = "bus", + pmax_mw = 100.0, + variable_cost_per_mwh = Float64(g), + commitment_required = false, + contingency_eligible = true, + reserve_capability_mw = Dict(:spin => 100.0), + ) for g = 1:12 + ] + data = DARTSystemData(generators = generators, network = single_bus()) + forecast = DARTForecast( + interval_hours = 1.0, + load_mw = fill(600.0, 1, 24), + availability = ones(12, 24), + ) + result = solve_dart_scuc(data, forecast, default_dart_state(data)) + + @test result.termination_status == JuMP.MOI.OPTIMAL + @test size(result.emergency_load_shed_mw) == (1, 24, 12) + @test sum(result.emergency_load_shed_mw) <= 1e-6 + @test all( + isapprox(sum(result.generation_mw[:, t]), 600.0; atol = 1e-6) for + t in axes(result.generation_mw, 2) + ) + + @test_throws ArgumentError build_dart_scuc_model( + data, + forecast, + default_dart_state(data); + config = DARTConfig(maximum_security_variables = 100), + ) +end + +@testset "DART production input validation" begin + generator = DARTGenerator( + name = "unit", + bus = "bus", + pmax_mw = 10.0, + commitment_required = false, + ) + data = DARTSystemData(generators = [generator], network = single_bus()) + forecast = + DARTForecast(interval_hours = 1.0, load_mw = zeros(1, 1), availability = ones(1, 1)) + + invalid_availability = DARTForecast( + interval_hours = 1.0, + load_mw = zeros(1, 1), + availability = fill(NaN, 1, 1), + ) + @test_throws ArgumentError build_dart_scuc_model( + data, + invalid_availability, + default_dart_state(data), + ) + + invalid_terminal = DARTForecast( + interval_hours = 1.0, + load_mw = zeros(1, 1), + availability = ones(1, 1), + terminal_storage_soc_mwh = [0.0], + ) + @test_throws ArgumentError build_dart_scuc_model( + data, + invalid_terminal, + default_dart_state(data), + ) + + invalid_generator = + DARTGenerator(name = "unit", bus = "bus", pmax_mw = 10.0, startup_limit_mw = -1.0) + invalid_generator_data = + DARTSystemData(generators = [invalid_generator], network = single_bus()) + @test_throws ArgumentError build_dart_scuc_model( + invalid_generator_data, + forecast, + default_dart_state(invalid_generator_data), + ) + + invalid_network = DARTNetwork( + bus_names = ["bus"], + line_names = ["line"], + ptdf = zeros(1, 1), + line_limit_mw = [-1.0], + emergency_line_limit_mw = [0.0], + ) + invalid_network_data = + DARTSystemData(generators = [generator], network = invalid_network) + @test_throws ArgumentError build_dart_scuc_model( + invalid_network_data, + forecast, + default_dart_state(invalid_network_data), + ) + + @test_throws ArgumentError build_dart_scuc_model( + data, + forecast, + default_dart_state(data); + config = DARTConfig(maximum_security_variables = 0), + ) + @test_throws ArgumentError run_dart_rolling( + data, + forecast, + forecast; + hours_to_run = 1, + day_ahead_lookahead_hours = 0, + ) +end