diff --git a/src/loop.jl b/src/loop.jl index 7bf05e27..7503196f 100644 --- a/src/loop.jl +++ b/src/loop.jl @@ -42,12 +42,11 @@ function interval_loop(factory_like, model_kwargs::Dict, storage = model_kwargs["storage"] demand_flexibility = model_kwargs["demand_flexibility"] sets = _make_sets(case, storage) - storage_enabled = (sets.num_storage > 0) unused_load_shed_intervals_turnoff = 14 # Start looping for i in 1:n_interval # These must be declared global so that they persist through the loop. - global m, voi, pg0, storage_e0, intervals_without_loadshed + global m, pg0, storage_e0, intervals_without_loadshed @show ("load_shed_enabled" in keys(model_kwargs)) @show ("BarHomogeneous" in keys(solver_kwargs)) interval_start = start_index + (i - 1) * interval @@ -60,22 +59,22 @@ function interval_loop(factory_like, model_kwargs::Dict, end if i == 1 # Build a model with no initial ramp constraint - if storage_enabled + if storage.enabled model_kwargs["storage_e0"] = storage.sd_table.InitialStorage end m = new_model(factory_like) JuMP.set_optimizer_attributes(m, pairs(solver_kwargs)...) - m, voi = _build_model(m; symbolize(model_kwargs)...) + m = _build_model(m; symbolize(model_kwargs)...) elseif i == 2 # Build a model with an initial ramp constraint model_kwargs["initial_ramp_enabled"] = true model_kwargs["initial_ramp_g0"] = pg0 - if storage_enabled + if storage.enabled model_kwargs["storage_e0"] = storage_e0 end m = new_model(factory_like) JuMP.set_optimizer_attributes(m, pairs(solver_kwargs)...) - m, voi = _build_model(m; symbolize(model_kwargs)...) + m = _build_model(m; symbolize(model_kwargs)...) else # Reassign right-hand-side of constraints to match profiles bus_demand = _make_bus_demand(case, interval_start, interval_end) @@ -87,48 +86,48 @@ function interval_loop(factory_like, model_kwargs::Dict, case.wind[interval_start:interval_end, 2:end])) for t in 1:interval, b in sets.load_bus_idx JuMP.set_normalized_rhs( - voi.powerbalance[b, t], bus_demand[b, t]) + m[:powerbalance][b, t], bus_demand[b, t]) end if (("load_shed_enabled" in keys(model_kwargs)) && (model_kwargs["load_shed_enabled"] == true)) for t in 1:interval, i in 1:length(sets.load_bus_idx) JuMP.set_normalized_rhs( - voi.load_shed_ub[i, t], bus_demand[sets.load_bus_idx[i], t] + m[:load_shed_ub][i, t], bus_demand[sets.load_bus_idx[i], t] ) end end for t in 1:interval, g in 1:sets.num_hydro JuMP.set_normalized_rhs( - voi.hydro_fixed[g, t], simulation_hydro[g, t]) + m[:hydro_fixed][g, t], simulation_hydro[g, t]) end for t in 1:interval, g in 1:sets.num_solar JuMP.set_normalized_rhs( - voi.solar_max[g, t], simulation_solar[g, t]) + m[:solar_max][g, t], simulation_solar[g, t]) end for t in 1:interval, g in 1:sets.num_wind JuMP.set_normalized_rhs( - voi.wind_max[g, t], simulation_wind[g, t]) + m[:wind_max][g, t], simulation_wind[g, t]) end # Re-assign right-hand-side for initial conditions noninf_ramp_idx = findall(case.gen_ramp30 .!= Inf) for g in noninf_ramp_idx rhs = case.gen_ramp30[g] * 2 + pg0[g] - JuMP.set_normalized_rhs(voi.initial_rampup[g], rhs) + JuMP.set_normalized_rhs(m[:initial_rampup][g], rhs) rhs = case.gen_ramp30[g] * 2 - pg0[g] - JuMP.set_normalized_rhs(voi.initial_rampdown[g], rhs) + JuMP.set_normalized_rhs(m[:initial_rampdown][g], rhs) end - if storage_enabled + if storage.enabled for s in 1:sets.num_storage - JuMP.set_normalized_rhs(voi.initial_soc[s], storage_e0[s]) + JuMP.set_normalized_rhs(m[:initial_soc][s], storage_e0[s]) end end if demand_flexibility.enabled for t in 1:interval, i in 1:length(sets.load_bus_idx) JuMP.set_upper_bound( - voi.load_shift_up[i, t], bus_flex_amt[sets.load_bus_idx[i], t] + m[:load_shift_up][i, t], bus_flex_amt[sets.load_bus_idx[i], t] ) JuMP.set_upper_bound( - voi.load_shift_dn[i, t], bus_flex_amt[sets.load_bus_idx[i], t] + m[:load_shift_dn][i, t], bus_flex_amt[sets.load_bus_idx[i], t] ) end end @@ -143,13 +142,13 @@ function interval_loop(factory_like, model_kwargs::Dict, status = JuMP.termination_status(m) if status == JuMP.MOI.OPTIMAL f = JuMP.objective_value(m) - results = get_results(f, voi, model_kwargs["case"]) + results = get_results(f, model_kwargs["case"]) break elseif ((status == JuMP.MOI.LOCALLY_SOLVED) & ("load_shed_enabled" in keys(model_kwargs))) # if load shedding is enabled, we'll accept 'suboptimal' f = JuMP.objective_value(m) - results = get_results(f, voi, model_kwargs["case"]) + results = get_results(f, model_kwargs["case"]) break elseif ((status in numeric_statuses) & (JuMP.solver_name(m) == "Gurobi") @@ -165,7 +164,7 @@ function interval_loop(factory_like, model_kwargs::Dict, println("rebuild with load shed") m = new_model(factory_like) JuMP.set_optimizer_attributes(m, pairs(solver_kwargs)...) - m, voi = _build_model(m; symbolize(model_kwargs)...) + m = _build_model(m; symbolize(model_kwargs)...) intervals_without_loadshed = 0 elseif ((JuMP.solver_name(m) == "Gurobi") & !("BarHomogeneous" in keys(solver_kwargs))) @@ -178,7 +177,7 @@ function interval_loop(factory_like, model_kwargs::Dict, println("rebuild with load shed") m = new_model(factory_like) JuMP.set_optimizer_attributes(m, pairs(solver_kwargs)...) - m, voi = _build_model(m; symbolize(model_kwargs)...) + m = _build_model(m; symbolize(model_kwargs)...) intervals_without_loadshed = 0 else # Something has gone very wrong @@ -189,7 +188,7 @@ function interval_loop(factory_like, model_kwargs::Dict, if (("load_shed_enabled" in keys(model_kwargs)) && (model_kwargs["load_shed_enabled"] == true)) # Display where load shedding is occurring - load_shed_values = JuMP.value.(voi.load_shed) + load_shed_values = JuMP.value.(m[:load_shed]) load_shed_indices = findall(load_shed_values .> 1e-6) if length(load_shed_indices) > 0 @show load_shed_indices @@ -203,7 +202,7 @@ function interval_loop(factory_like, model_kwargs::Dict, # Save initial conditions for next interval pg0 = results.pg[:,end] - if storage_enabled + if storage.enabled storage_e0 = results.storage_e[:,end] end @@ -228,7 +227,7 @@ function interval_loop(factory_like, model_kwargs::Dict, delete!(model_kwargs, "load_shed_enabled") m = new_model(factory_like) JuMP.set_optimizer_attributes(m, pairs(solver_kwargs)...) - m, voi = _build_model(m; symbolize(model_kwargs)...) + m = _build_model(m; symbolize(model_kwargs)...) end end end diff --git a/src/model.jl b/src/model.jl index e4c799ae..c1f4edb0 100644 --- a/src/model.jl +++ b/src/model.jl @@ -170,7 +170,7 @@ function _make_sets(case::Case, storage::Union{Storage,Nothing})::Sets num_segments = convert(Int, maximum(case.gencost[:, NCOST])) - 1 segment_idx = 1:num_segments # Storage - storage_enabled = isa(storage, Storage) && (size(storage.gen, 1) > 0) + storage_enabled = isa(storage, Storage) && storage.enabled num_storage = storage_enabled ? size(storage.gen, 1) : 0 storage_idx = storage_enabled ? (1:num_storage) : nothing @@ -190,6 +190,357 @@ function _make_sets(case::Case, storage::Union{Storage,Nothing})::Sets end +function _add_constraint_power_balance!( + m::JuMP.Model, + case::Case, + sets::Sets, + storage::Storage, + demand_flexibility::DemandFlexibility, + bus_demand::Matrix, + load_shed_enabled::Bool, +) + # Generator topology matrix + gen_map = _make_gen_map(case) + # Branch connectivity matrix + branch_map = _make_branch_map(case) + + gen_injections = JuMP.@expression(m, gen_map * m[:pg]) + line_injections = JuMP.@expression(m, branch_map * m[:pf]) + injections = JuMP.@expression(m, gen_injections + line_injections) + if load_shed_enabled + injections = JuMP.@expression(m, injections + sets.load_bus_map * m[:load_shed]) + end + withdrawals = JuMP.@expression(m, bus_demand) + if storage.enabled + storage_bus_idx = [sets.bus_id2idx[b] for b in storage.gen[:, 1]] + storage_map = sparse( + storage_bus_idx, sets.storage_idx, 1, sets.num_bus, sets.num_storage + )::SparseMatrixCSC + injections = JuMP.@expression(m, injections + storage_map * m[:storage_dis]) + withdrawals = JuMP.@expression(m, withdrawals + storage_map * m[:storage_chg]) + end + if demand_flexibility.enabled + injections = JuMP.@expression( + m, injections + sets.load_bus_map * m[:load_shift_dn] + ) + withdrawals = JuMP.@expression( + m, withdrawals + sets.load_bus_map * m[:load_shift_up] + ) + end + JuMP.@constraint(m, powerbalance, (injections .== withdrawals)) + println("powerbalance, setting names: ", Dates.now()) + interval_length = size(bus_demand)[2] + for i in sets.bus_idx, j in 1:interval_length + JuMP.set_name( + powerbalance[i, j], "powerbalance[" * string(i) * "," * string(j) * "]" + ) + end +end + + +function _add_constraint_load_shed!( + m::JuMP.Model, + case::Case, + sets::Sets, + demand_flexibility::DemandFlexibility, + bus_demand::Matrix, +) + demand_for_load_shed = JuMP.@expression( + m, + [i=1:sets.num_load_bus, j=1:interval_length], + bus_demand[sets.load_bus_idx[i], j], + ) + if demand_flexibility.enabled + demand_for_load_shed = JuMP.@expression( + m, + [i=1:sets.num_load_bus, j=1:interval_length], + demand_for_load_shed[i, j] + m[:load_shift_up][i, j] - m[:load_shift_dn][i, j], + ) + end + JuMP.@constraint( + m, + load_shed_ub[i in 1:sets.num_load_bus, j in 1:interval_length], + m[:load_shed][i, j] <= demand_for_load_shed[i, j], + ) +end + + +function _add_constraints_storage_operation!( + m::JuMP.Model, + case::Case, + sets::Sets, + storage::Storage, + interval_length::Int, + storage_e0::Array{Float64,1}, +) + num_hour = interval_length + + println("storage soc_tracking: ", Dates.now()) + JuMP.@constraint( + m, + soc_tracking[i in 1:sets.num_storage, h in 1:(num_hour-1)], + m[:storage_soc][i, h+1] == ( + m[:storage_soc][i, h] * (1 - storage.sd_table.LossFactor[i]) + + storage.sd_table.InEff[i] * m[:storage_chg][i, h+1] + - (1 / storage.sd_table.OutEff[i]) * m[:storage_dis][i, h+1] + ), + container=Array, + ) + println("storage initial_soc: ", Dates.now()) + JuMP.@constraint( + m, + initial_soc[i in 1:sets.num_storage], + m[:storage_soc][i, 1] == ( + storage_e0[i] + + storage.sd_table.InEff[i] * m[:storage_chg][i, 1] + - (1 / storage.sd_table.OutEff[i]) * m[:storage_dis][i, 1] + ), + container=Array, + ) + println("storage final_soc_min: ", Dates.now()) + JuMP.@constraint( + m, + soc_terminal_min[i in 1:sets.num_storage], + m[:storage_soc][i, num_hour] >= storage.sd_table.ExpectedTerminalStorageMin[i], + container=Array, + ) + println("storage final_soc_max: ", Dates.now()) + JuMP.@constraint( + m, + soc_terminal_max[i in 1:sets.num_storage], + m[:storage_soc][i, num_hour] <= storage.sd_table.ExpectedTerminalStorageMax[i], + container=Array, + ) +end + + +function _add_constraints_demand_flexibility!( + m::JuMP.Model, + case::Case, + sets::Sets, + demand_flexibility::DemandFlexibility, + interval_length::Int, +) + if demand_flexibility.rolling_balance && ( + demand_flexibility.duration < interval_length + ) + println("rolling load balance: ", Dates.now()) + JuMP.@constraint( + m, + rolling_load_balance[ + i in 1:sets.num_load_bus, + k in 1:(interval_length - demand_flexibility.duration) + ], + sum( + m[:load_shift_up][i, j] - m[:load_shift_dn][i, j] + for j in k:(k + demand_flexibility.duration) + ) >= 0, + ) + end + if demand_flexibility.interval_balance + println("interval load balance: ", Dates.now()) + JuMP.@constraint( + m, + interval_load_balance[i in 1:sets.num_load_bus], + sum( + m[:load_shift_up][i, j] - m[:load_shift_dn][i, j] for j in 1:interval_length + ) >= 0, + ) + end +end + + +function _add_constraints_initial_ramping!( + m::JuMP.Model, case::Case, sets::Sets, initial_ramp_g0 +) + noninf_ramp_idx = findall(case.gen_ramp30 .!= Inf) + println("initial rampup: ", Dates.now()) + JuMP.@constraint( + m, + initial_rampup[i in noninf_ramp_idx], + m[:pg][i, 1] - initial_ramp_g0[i] <= case.gen_ramp30[i] * 2, + ) + println("initial rampdown: ", Dates.now()) + JuMP.@constraint( + m, + initial_rampdown[i in noninf_ramp_idx], + case.gen_ramp30[i] * -2 <= m[:pg][i, 1] - initial_ramp_g0[i], + ) +end + + +function _add_constraints_ramping!( + m::JuMP.Model, case::Case, sets::Sets, interval_length::Int +) + noninf_ramp_idx = findall(case.gen_ramp30 .!= Inf) + println("rampup: ", Dates.now()) + JuMP.@constraint( + m, + rampup[i in noninf_ramp_idx, h in 1:(interval_length-1)], + m[:pg][i, h+1] - m[:pg][i, h] <= case.gen_ramp30[i] * 2, + ) + println("rampdown: ", Dates.now()) + JuMP.@constraint( + m, + rampdown[i in noninf_ramp_idx, h in 1:(interval_length-1)], + case.gen_ramp30[i] * -2 <= m[:pg][i, h+1] - m[:pg][i, h], + ) +end + + +function _add_constraints_generator_segments!( + m::JuMP.Model, case::Case, sets::Sets, hour_idx, +) + segment_width = (case.gen_pmax - case.gen_pmin) ./ sets.num_segments + println("segment_max: ", Dates.now()) + JuMP.@constraint( + m, + segment_max[i in sets.noninf_pmax, s in sets.segment_idx, h in hour_idx], + m[:pg_seg][i, s, h] <= segment_width[i], + ) + println("segment_add: ", Dates.now()) + JuMP.@constraint( + m, + segment_add[i in sets.noninf_pmax, h in hour_idx], + m[:pg][i, h] == case.gen_pmin[i] + sum(m[:pg_seg][i, :, h]), + ) +end + + +function _add_constraints_branch_flow_limits!( + m::JuMP.Model, case::Case, sets::Sets, trans_viol_enabled, hour_idx +) + branch_pmin = vcat(-1 * case.branch_rating, case.dcline_pmin) + branch_pmax = vcat(case.branch_rating, case.dcline_pmax) + if trans_viol_enabled + JuMP.@expression(m, branch_limit_pmin, branch_pmin - m[:trans_viol]) + JuMP.@expression(m, branch_limit_pmax, branch_pmax + m[:trans_viol]) + else + JuMP.@expression( + m, + branch_limit_pmin[br in sets.branch_idx, h in hour_idx], + branch_pmin[br], + ) + JuMP.@expression( + m, + branch_limit_pmax[br in sets.branch_idx, h in hour_idx], + branch_pmax[br], + ) + end + println("branch_min, branch_max: ", Dates.now()) + JuMP.@constraint( + m, + branch_min[br in sets.noninf_branch_idx, h in hour_idx], + branch_limit_pmin[br, h] <= m[:pf][br, h], + ) + println("branch_max: ", Dates.now()) + JuMP.@constraint( + m, + branch_max[br in sets.noninf_branch_idx, h in hour_idx], + m[:pf][br, h] <= branch_limit_pmax[br, h], + ) +end + + +function _add_branch_angle_constraints!(m::JuMP.Model, case::Case, sets::Sets, hour_idx) + # Explicit numbering here so that we constrain AC branches but not DC + JuMP.@constraint( + m, + branch_angle[br in 1:sets.num_branch_ac, h in hour_idx], + ( + case.branch_reactance[br] * m[:pf][br, h] + == ( + m[:theta][sets.branch_to_idx[br], h] + - m[:theta][sets.branch_from_idx[br], h] + ) + ) + ) +end + + +function _add_profile_generator_limits!( + m::JuMP.Model, case::Case, sets::Sets, hour_idx, start_index, interval_length +) + end_index = start_index + interval_length - 1 + # Generation segments + simulation_hydro = Matrix(case.hydro[start_index:end_index, 2:end]) + simulation_solar = Matrix(case.solar[start_index:end_index, 2:end]) + simulation_wind = Matrix(case.wind[start_index:end_index, 2:end]) + println("hydro_fixed: ", Dates.now()) + JuMP.@constraint( + m, + hydro_fixed[i in 1:sets.num_hydro, h in hour_idx], + m[:pg][sets.gen_hydro_idx[i], h] == simulation_hydro[h, i], + ) + println("solar_max: ", Dates.now()) + JuMP.@constraint( + m, + solar_max[i in 1:sets.num_solar, h in hour_idx], + m[:pg][sets.gen_solar_idx[i], h] <= simulation_solar[h, i], + ) + println("wind_max: ", Dates.now()) + JuMP.@constraint( + m, + wind_max[i in 1:sets.num_wind, h in hour_idx], + m[:pg][sets.gen_wind_idx[i], h] <= simulation_wind[h, i], + ) +end + + +function _add_objective_function!( + m::JuMP.Model, + case::Case, + sets::Sets, + storage::Storage, + interval_length::Int, + load_shed_enabled::Bool, + load_shed_penalty::Number, + trans_viol_enabled::Bool, + storage_e0::Array{Float64,1}, +) + # Positional indices from mpc.gencost + COST = 5 + fixed_cost = case.gencost[:, COST+1] + segment_width = (case.gen_pmax - case.gen_pmin) ./ sets.num_segments + segment_slope = _build_segment_slope(case, sets.segment_idx, segment_width) + + # Start with generator variable O & M, piecewise + obj = JuMP.@expression( + m, + sum(segment_slope[sets.noninf_pmax, :] .* m[:pg_seg][sets.noninf_pmax, :, :]), + ) + # Add fixed costs + JuMP.add_to_expression!( + obj, JuMP.@expression(m, interval_length * sum(fixed_cost)) + ) + # Add load shed penalty (if necessary) + if load_shed_enabled + JuMP.add_to_expression!( + obj, JuMP.@expression(m, load_shed_penalty * sum(m[:load_shed])) + ) + end + # Add transmission violation penalty (if necessary) + if trans_viol_enabled + JuMP.add_to_expression!( + obj, JuMP.@expression(m, trans_viol_penalty * sum(m[:trans_viol])) + ) + end + # Pay for ending with less storage energy than initial + if storage.enabled + storage_penalty = JuMP.@expression( + m, + sum( + (storage_e0 - m[:storage_soc][:, end]) + .* storage.sd_table.TerminalStoragePrice + ), + ) + JuMP.add_to_expression!(obj, storage_penalty) + end + # Finally, set as objective of model + JuMP.@objective(m, Min, obj) +end + + """ _build_model(m; case=case, storage=storage, start_index=x, interval_length=y[, kwargs...]) @@ -212,16 +563,13 @@ function _build_model( initial_ramp_enabled::Bool=false, initial_ramp_g0::Array{Float64,1}=Float64[], storage_e0::Array{Float64,1}=Float64[] -)::Tuple{JuMP.Model, VariablesOfInterest} - # Positional indices from mpc.gencost - COST = 5 +)::JuMP.Model # Positional indices from mpc.gen PMAX = 9 PMIN = 10 println("building sets: ", Dates.now()) # Sets - time periods - num_hour = interval_length hour_idx = 1:interval_length end_index = start_index + interval_length - 1 # Sets - static @@ -229,33 +577,7 @@ function _build_model( println("parameters: ", Dates.now()) # Parameters - # Generator topology matrix - gen_map = _make_gen_map(case) - # Generation segments - segment_width = (case.gen_pmax - case.gen_pmin) ./ sets.num_segments - fixed_cost = case.gencost[:, COST+1] - segment_slope = _build_segment_slope(case, sets.segment_idx, segment_width) - # Branch connectivity matrix - branch_map = _make_branch_map(case) - branch_pmin = vcat(-1 * case.branch_rating, case.dcline_pmin) - branch_pmax = vcat(case.branch_rating, case.dcline_pmax) - # Demand by bus - bus_demand = _make_bus_demand(case, start_index, end_index) - bus_demand *= demand_scaling - simulation_hydro = Matrix(case.hydro[start_index:end_index, 2:end]) - simulation_solar = Matrix(case.solar[start_index:end_index, 2:end]) - simulation_wind = Matrix(case.wind[start_index:end_index, 2:end]) - # Storage parameters (if present) - storage_enabled = (sets.num_storage > 0) - if storage_enabled - storage_max_dis = storage.gen[:, PMAX] - storage_max_chg = -1 * storage.gen[:, PMIN] - storage_min_energy = storage.sd_table.MinStorageLevel - storage_max_energy = storage.sd_table.MaxStorageLevel - storage_bus_idx = [sets.bus_id2idx[b] for b in storage.gen[:, 1]] - storage_map = sparse(storage_bus_idx, sets.storage_idx, 1, - sets.num_bus, sets.num_storage)::SparseMatrixCSC - end + bus_demand = _make_bus_demand(case, start_index, end_index) * demand_scaling # Demand flexibility parameters (if present) if demand_flexibility.enabled bus_demand_flex_amt = _make_bus_demand_flexibility_amount( @@ -284,7 +606,11 @@ function _build_model( 0 <= trans_viol[i in sets.branch_idx, j in hour_idx], container=Array) end - if storage_enabled + if storage.enabled + storage_max_dis = storage.gen[:, PMAX] + storage_max_chg = -1 * storage.gen[:, PMIN] + storage_min_energy = storage.sd_table.MinStorageLevel + storage_max_energy = storage.sd_table.MaxStorageLevel JuMP.@variables(m, begin (0 <= storage_chg[i in sets.storage_idx, j in hour_idx] <= storage_max_chg[i]), (container=Array) @@ -315,231 +641,56 @@ function _build_model( # Constraints println("powerbalance: ", Dates.now()) - gen_injections = JuMP.@expression(m, gen_map * pg) - line_injections = JuMP.@expression(m, branch_map * pf) - injections = JuMP.@expression(m, gen_injections + line_injections) - if load_shed_enabled - injections = JuMP.@expression(m, injections + sets.load_bus_map * load_shed) - end - withdrawals = JuMP.@expression(m, bus_demand) - if storage_enabled - injections = JuMP.@expression(m, injections + storage_map * storage_dis) - withdrawals = JuMP.@expression(m, withdrawals + storage_map * storage_chg) - end - if demand_flexibility.enabled - injections = JuMP.@expression(m, injections + sets.load_bus_map * load_shift_dn) - withdrawals = JuMP.@expression( - m, withdrawals + sets.load_bus_map * load_shift_up - ) - end - JuMP.@constraint(m, powerbalance, (injections .== withdrawals)) - println("powerbalance, setting names: ", Dates.now()) - for i in sets.bus_idx, j in hour_idx - JuMP.set_name(powerbalance[i, j], - "powerbalance[" * string(i) * "," * string(j) * "]") - end + _add_constraint_power_balance!( + m, case, sets, storage, demand_flexibility, bus_demand, load_shed_enabled, + ) if load_shed_enabled - demand_for_load_shed = JuMP.@expression( - m, - [i=1:sets.num_load_bus, j=1:interval_length], - bus_demand[sets.load_bus_idx[i], j] - ) - if demand_flexibility.enabled - demand_for_load_shed = JuMP.@expression( - m, - [i=1:sets.num_load_bus, j=1:interval_length], - demand_for_load_shed[i, j] + load_shift_up[i, j] - load_shift_dn[i, j] - ) - end - JuMP.@constraint( - m, - load_shed_ub[i in 1:sets.num_load_bus, j in 1:interval_length], - load_shed[i, j] <= demand_for_load_shed[i, j] - ) + _add_constraint_load_shed!(m, case, sets, demand_flexibility, bus_demand) end - if storage_enabled - println("storage soc_tracking: ", Dates.now()) - JuMP.@constraint(m, - soc_tracking[i in 1:sets.num_storage, h in 1:(num_hour-1)], - storage_soc[i, h+1] == ( - storage_soc[i, h] * (1 - storage.sd_table.LossFactor[i]) - + storage.sd_table.InEff[i] * storage_chg[i, h+1] - - (1 / storage.sd_table.OutEff[i]) * storage_dis[i, h+1]), - container=Array) - println("storage initial_soc: ", Dates.now()) - JuMP.@constraint(m, - initial_soc[i in 1:sets.num_storage], - storage_soc[i, 1] == ( - storage_e0[i] - + storage.sd_table.InEff[i] * storage_chg[i, 1] - - (1 / storage.sd_table.OutEff[i]) * storage_dis[i, 1]), - container=Array) - println("storage final_soc_min: ", Dates.now()) - JuMP.@constraint(m, - soc_terminal_min[i in 1:sets.num_storage], - storage_soc[i, num_hour] >= storage.sd_table.ExpectedTerminalStorageMin[i], - container=Array) - println("storage final_soc_max: ", Dates.now()) - JuMP.@constraint(m, - soc_terminal_max[i in 1:sets.num_storage], - storage_soc[i, num_hour] <= storage.sd_table.ExpectedTerminalStorageMax[i], - container=Array) + if storage.enabled + _add_constraints_storage_operation!( + m, case, sets, storage, interval_length, storage_e0 + ) end if demand_flexibility.enabled - if demand_flexibility.rolling_balance && ( - demand_flexibility.duration < interval_length + _add_constraints_demand_flexibility!( + m, case, sets, demand_flexibility, interval_length ) - println("rolling load balance: ", Dates.now()) - JuMP.@constraint( - m, - rolling_load_balance[ - i in 1:sets.num_load_bus, - k in 1:(interval_length - demand_flexibility.duration) - ], - sum( - load_shift_up[i, j] - load_shift_dn[i, j] - for j in k:(k + demand_flexibility.duration) - ) >= 0 - ) - end - if demand_flexibility.interval_balance - println("interval load balance: ", Dates.now()) - JuMP.@constraint( - m, - interval_load_balance[i in 1:sets.num_load_bus], - sum( - load_shift_up[i, j] - load_shift_dn[i, j] for j in 1:interval_length - ) >= 0 - ) - end end - noninf_ramp_idx = findall(case.gen_ramp30 .!= Inf) if initial_ramp_enabled - println("initial rampup: ", Dates.now()) - JuMP.@constraint(m, - initial_rampup[i in noninf_ramp_idx], - pg[i, 1] - initial_ramp_g0[i] <= case.gen_ramp30[i] * 2) - println("initial rampdown: ", Dates.now()) - JuMP.@constraint(m, - initial_rampdown[i in noninf_ramp_idx], - case.gen_ramp30[i] * -2 <= pg[i, 1] - initial_ramp_g0[i]) + _add_constraints_initial_ramping!(m, case, sets, initial_ramp_g0) end if length(hour_idx) > 1 - println("rampup: ", Dates.now()) - JuMP.@constraint(m, - rampup[i in noninf_ramp_idx, h in 1:(num_hour-1)], - pg[i, h+1] - pg[i, h] <= case.gen_ramp30[i] * 2) - println("rampdown: ", Dates.now()) - JuMP.@constraint(m, - rampdown[i in noninf_ramp_idx, h in 1:(num_hour-1)], - case.gen_ramp30[i] * -2 <= pg[i, h+1] - pg[i, h]) + _add_constraints_ramping!(m, case, sets, interval_length) end - println("segment_max: ", Dates.now()) - JuMP.@constraint(m, - segment_max[ - i in sets.noninf_pmax, s in sets.segment_idx, h in hour_idx], - pg_seg[i, s, h] <= segment_width[i]) - println("segment_add: ", Dates.now()) - JuMP.@constraint(m, - segment_add[i in sets.noninf_pmax, h in hour_idx], - pg[i, h] == case.gen_pmin[i] + sum(pg_seg[i, :, h])) + _add_constraints_generator_segments!(m, case, sets, hour_idx) - if trans_viol_enabled - JuMP.@expression(m, branch_limit_pmin, branch_pmin - trans_viol) - JuMP.@expression(m, branch_limit_pmax, branch_pmax + trans_viol) - else - JuMP.@expression(m, - branch_limit_pmin[br in sets.branch_idx, h in hour_idx], - branch_pmin[br]) - JuMP.@expression(m, - branch_limit_pmax[br in sets.branch_idx, h in hour_idx], - branch_pmax[br]) - end - println("branch_min, branch_max: ", Dates.now()) - JuMP.@constraint(m, - branch_min[br in sets.noninf_branch_idx, h in hour_idx], - branch_limit_pmin[br, h] <= pf[br, h]) - println("branch_max: ", Dates.now()) - JuMP.@constraint(m, - branch_max[br in sets.noninf_branch_idx, h in hour_idx], - pf[br, h] <= branch_limit_pmax[br, h]) + _add_constraints_branch_flow_limits!(m, case, sets, trans_viol_enabled, hour_idx) println("branch_angle: ", Dates.now()) - # Explicit numbering here so that we constrain AC branches but not DC - JuMP.@constraint(m, - branch_angle[br in 1:sets.num_branch_ac, h in hour_idx], - (case.branch_reactance[br] * pf[br, h] - == (theta[sets.branch_to_idx[br], h] - - theta[sets.branch_from_idx[br], h]))) + _add_branch_angle_constraints!(m, case, sets, hour_idx) # Constrain variable generators based on profiles - println("hydro_fixed: ", Dates.now()) - JuMP.@constraint(m, - hydro_fixed[i in 1:sets.num_hydro, h in hour_idx], - pg[sets.gen_hydro_idx[i], h] == simulation_hydro[h, i]) - println("solar_max: ", Dates.now()) - JuMP.@constraint(m, - solar_max[i in 1:sets.num_solar, h in hour_idx], - pg[sets.gen_solar_idx[i], h] <= simulation_solar[h, i]) - println("wind_max: ", Dates.now()) - JuMP.@constraint(m, - wind_max[i in 1:sets.num_wind, h in hour_idx], - pg[sets.gen_wind_idx[i], h] <= simulation_wind[h, i]) + _add_profile_generator_limits!(m, case, sets, hour_idx, start_index, interval_length) println("objective: ", Dates.now()) - # Start with generator variable O & M, piecewise - obj = JuMP.@expression(m, - sum(segment_slope[sets.noninf_pmax, :] - .* pg_seg[sets.noninf_pmax, :, :])) - # Add fixed costs - JuMP.add_to_expression!( - obj, JuMP.@expression(m, num_hour * sum(fixed_cost))) - # Add load shed penalty (if necessary) - if load_shed_enabled - JuMP.add_to_expression!( - obj, JuMP.@expression(m, load_shed_penalty * sum(load_shed))) - end - # Add transmission violation penalty (if necessary) - if trans_viol_enabled - JuMP.add_to_expression!( - obj, JuMP.@expression(m, trans_viol_penalty * sum(trans_viol))) - end - # Pay for ending with less storage energy than initial - if storage_enabled - storage_penalty = JuMP.@expression(m, - sum((storage_e0 - storage_soc[:, end]) - .* storage.sd_table.TerminalStoragePrice)) - JuMP.add_to_expression!(obj, storage_penalty) - end - # Finally, set as objective of model - JuMP.@objective(m, Min, obj) + _add_objective_function!( + m, + case, + sets, + storage, + interval_length, + load_shed_enabled, + load_shed_penalty, + trans_viol_enabled, + storage_e0, + ) println(Dates.now()) - # For non-existent variables/constraints, define as `nothing` - load_shed = load_shed_enabled ? load_shed : nothing - load_shift_up = demand_flexibility.enabled ? load_shift_up : nothing - load_shift_dn = demand_flexibility.enabled ? load_shift_dn : nothing - storage_dis = storage_enabled ? storage_dis : nothing - storage_chg = storage_enabled ? storage_chg : nothing - storage_soc = storage_enabled ? storage_soc : nothing - initial_soc = storage_enabled ? initial_soc : nothing - initial_rampup = initial_ramp_enabled ? initial_rampup : nothing - initial_rampdown = initial_ramp_enabled ? initial_rampdown : nothing - load_shed_ub = load_shed_enabled ? load_shed_ub : nothing - voi = VariablesOfInterest(; - # Variables - pg=pg, pf=pf, - load_shed=load_shed, load_shift_up=load_shift_up, load_shift_dn=load_shift_dn, - storage_soc=storage_soc, storage_dis=storage_dis, storage_chg=storage_chg, - # Constraints - branch_min=branch_min, branch_max=branch_max, powerbalance=powerbalance, - initial_soc=initial_soc, load_shed_ub=load_shed_ub, - initial_rampup=initial_rampup, initial_rampdown=initial_rampdown, - hydro_fixed=hydro_fixed, solar_max=solar_max, wind_max=wind_max) - return (m, voi) + return m end diff --git a/src/query.jl b/src/query.jl index 3476516f..f0ce024d 100644 --- a/src/query.jl +++ b/src/query.jl @@ -1,17 +1,17 @@ """ - get_results(model) + get_results(model, case) Extract the results of a simulation, store in a struct. """ -function get_results(f::Float64, voi::VariablesOfInterest, case::Case)::Results +function get_results(f::Float64, case::Case)::Results status = "OPTIMAL" sets = _make_sets(case) # These variables will always be in the results - pg = JuMP.value.(voi.pg) - pf = JuMP.value.(voi.pf) - lmp = -1 * JuMP.shadow_price.(voi.powerbalance) - congl_temp = -1 * JuMP.shadow_price.(voi.branch_min) - congu_temp = -1 * JuMP.shadow_price.(voi.branch_max) + pg = JuMP.value.(m[:pg]) + pf = JuMP.value.(m[:pf]) + lmp = -1 * JuMP.shadow_price.(m[:powerbalance]) + congl_temp = -1 * JuMP.shadow_price.(m[:branch_min]) + congu_temp = -1 * JuMP.shadow_price.(m[:branch_max]) # If DC lines are present, separate their results # Initialize with empty arrays, to be discarded later if they stay empty pf_dcline = zeros(0, 0) @@ -34,13 +34,13 @@ function get_results(f::Float64, voi::VariablesOfInterest, case::Case)::Results storage_pg = zeros(0, 0) storage_e = zeros(0, 0) try - storage_dis = JuMP.value.(voi.storage_dis) - storage_chg = JuMP.value.(voi.storage_chg) - storage_e = JuMP.value.(voi.storage_soc) + storage_dis = JuMP.value.(m[:storage_dis]) + storage_chg = JuMP.value.(m[:storage_chg]) + storage_e = JuMP.value.(m[:storage_soc]) storage_pg = storage_dis - storage_chg catch e - if isa(e, MethodError) - # Thrown when storage variables are `nothing` + if isa(e, KeyError) + # Thrown when storage variables are not defined in the model else # Unknown error, rethrow it rethrow(e) @@ -51,11 +51,11 @@ function get_results(f::Float64, voi::VariablesOfInterest, case::Case)::Results # Initialize with empty arrays, to be discarded later if they stay empty load_shed = zeros(0, 0) try - load_shed_temp = JuMP.value.(voi.load_shed) + load_shed_temp = JuMP.value.(m[:load_shed]) load_shed = sets.load_bus_map * load_shed_temp catch e - if isa(e, MethodError) - # Thrown when load_shed is `nothing` + if isa(e, KeyError) + # Thrown when load_shed is not defined in the model else # Unknown error, rethrow it rethrow(e) @@ -67,12 +67,12 @@ function get_results(f::Float64, voi::VariablesOfInterest, case::Case)::Results load_shift_up = zeros(0, 0) load_shift_dn = zeros(0, 0) try - load_shift_up_temp = JuMP.value.(voi.load_shift_up) - load_shift_dn_temp = JuMP.value.(voi.load_shift_dn) + load_shift_up_temp = JuMP.value.(m[:load_shift_up]) + load_shift_dn_temp = JuMP.value.(m[:load_shift_dn]) load_shift_up = sets.load_bus_map * load_shift_up_temp load_shift_dn = sets.load_bus_map * load_shift_dn_temp catch e - if isa(e, MethodError) + if isa(e, KeyError) # Thrown when load shift variables are `nothing` else # Unknown error, rethrow it diff --git a/src/read.jl b/src/read.jl index 80ecbcf9..cac4e871 100644 --- a/src/read.jl +++ b/src/read.jl @@ -75,7 +75,9 @@ end """Read input matfile (if present), return parsed data in a Storage struct.""" function read_storage(filepath)::Storage # Fallback dataframe, in case there's no case_storage.mat file - storage = Dict("gen" => zeros(0, 21), "sd_table" => DataFrames.DataFrame()) + storage = Dict( + "enabled" => false, "gen" => zeros(0, 21), "sd_table" => DataFrames.DataFrame() + ) try case_storage_file = MAT.matopen(joinpath(filepath, "case_storage.mat")) storage_mat_data = read(case_storage_file, "storage") @@ -83,6 +85,7 @@ function read_storage(filepath)::Storage # Convert N x 1 array of strings into 1D array of Symbols (length N) column_symbols = Symbol.(vec(storage_mat_data["sd_table"]["colnames"])) storage = Dict( + "enabled" => true, "gen" => storage_mat_data["gen"], "sd_table" => DataFrames.DataFrame( storage_mat_data["sd_table"]["data"], column_symbols diff --git a/src/types.jl b/src/types.jl index a934f798..31ff587d 100644 --- a/src/types.jl +++ b/src/types.jl @@ -37,6 +37,7 @@ end Base.@kwdef struct Storage + enabled::Bool gen::Array{Float64,2} sd_table::DataFrames.DataFrame end @@ -69,28 +70,6 @@ Base.@kwdef struct Results end -Base.@kwdef struct VariablesOfInterest - pg::Array{JuMP.VariableRef,2} - pf::Array{JuMP.VariableRef,2} - load_shed::Union{Array{JuMP.VariableRef,2},Nothing} - load_shed_ub::Union{JuMP.Containers.DenseAxisArray,Nothing} - load_shift_up::Union{Array{JuMP.VariableRef,2},Nothing} - load_shift_dn::Union{Array{JuMP.VariableRef,2},Nothing} - powerbalance::Array{JuMP.ConstraintRef,2} - branch_min::JuMP.Containers.DenseAxisArray - branch_max::JuMP.Containers.DenseAxisArray - initial_rampup::Union{JuMP.Containers.DenseAxisArray,Nothing} - initial_rampdown::Union{JuMP.Containers.DenseAxisArray,Nothing} - hydro_fixed::JuMP.Containers.DenseAxisArray - solar_max::JuMP.Containers.DenseAxisArray - wind_max::JuMP.Containers.DenseAxisArray - storage_dis::Union{Array{JuMP.VariableRef,2},Nothing} - storage_chg::Union{Array{JuMP.VariableRef,2},Nothing} - storage_soc::Union{Array{JuMP.VariableRef,2},Nothing} - initial_soc::Union{Array{JuMP.ConstraintRef,1},Nothing} -end - - Base.@kwdef struct Sets # Branch and branch subsets num_branch::Int64