ANL-CEEESA/UnitCommitment.jl
1
0
Fork 0
mirror of https://github.com/ANL-CEEESA/UnitCommitment.jl.git synced 2025-12-06 00:08:52 -06:00
Code Issues Packages Projects Releases Wiki Activity

Added comments on formulations, added start/stop constraints into MorLatRam and GenMorRam, added ability to add shortfall penalty.

Browse Source
This commit is contained in:
Aleksandr Kazachkov
2021-07-21 14:52:58 -04:00
parent 7a1b6f0f55
commit 483c679c4e
24 changed files with 919 additions and 86 deletions
Download Patch File Download Diff File Expand all files Collapse all files

5
src/instance/read.jl
View File

@@ -98,6 +98,10 @@ function _from_json(json; repair = true)
json["Parameters"]["Power balance penalty (\$/MW)"],
default = [1000.0 for t in 1:T],
)
shortfall_penalty = timeseries(
json["Parameters"]["Reserve shortfall penalty (\$/MW)"],
default=[0. for t in 1:T]
)
# Read buses
for (bus_name, dict) in json["Buses"]
@@ -264,6 +268,7 @@ function _from_json(json; repair = true)
instance = UnitCommitmentInstance(
T,
power_balance_penalty,
shortfall_penalty,
units,
buses,
lines,

2
src/instance/structs.jl
View File

@@ -72,6 +72,8 @@ end
mutable struct UnitCommitmentInstance
time::Int
power_balance_penalty::Vector{Float64}
"Penalty for failing to meet reserve requirement."
shortfall_penalty::Vector{Float64}
units::Vector{Unit}
buses::Vector{Bus}
lines::Vector{TransmissionLine}

2
src/model/build.jl
View File

@@ -21,6 +21,8 @@ Arguments
- `optimizer`:
the optimizer factory that should be attached to this model (e.g. Cbc.Optimizer).
If not provided, no optimizer will be attached.
- `formulation`:
the details of which constraints, variables, etc. to use
- `variable_names`:
If true, set variable and constraint names. Important if the model is going
to be exported to an MPS file. For large models, this can take significant

25
src/model/formulations/ArrCon2000/ramp.jl
View File

@@ -2,6 +2,26 @@
# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details.
"""
_add_ramp_eqs!
Ensure constraints on ramping are met.
Based on Arroyo and Conejo (2000).
Eqns. (24), (25) in Kneuven et al. (2020).
Variables
---
* :is_on
* :switch_off
* :switch_on
* :prod_above
* :reserve
Constraints
---
* :eq_ramp_up
* :eq_ramp_down
"""
function _add_ramp_eqs!(
model::JuMP.Model,
g::Unit,
@@ -22,7 +42,6 @@ function _add_ramp_eqs!(
reserve = model[:reserve]
eq_ramp_down = _init(model, :eq_ramp_down)
eq_ramp_up = _init(model, :eq_ramp_up)
is_initially_on = (g.initial_status > 0)
# Gar1962.ProdVars
prod_above = model[:prod_above]
@@ -35,7 +54,7 @@ function _add_ramp_eqs!(
for t in 1:model[:instance].time
# Ramp up limit
if t == 1
if is_initially_on
if _is_initially_on(g)
# min power is _not_ multiplied by is_on because if !is_on, then ramp up is irrelevant
eq_ramp_up[gn, t] = @constraint(
model,
@@ -66,7 +85,7 @@ function _add_ramp_eqs!(
# Ramp down limit
if t == 1
if is_initially_on
if _is_initially_on(g)
# TODO If RD < SD, or more specifically if
# min_power + RD < initial_power < SD
# then the generator should be able to shut down at time t = 1,

22
src/model/formulations/CarArr2006/pwlcosts.jl
View File

@@ -2,6 +2,28 @@
# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details.
"""
_add_production_piecewise_linear_eqs!
Ensure respect of production limits along each segment.
Based on Garver (1962) and Carrión and Arryo (2006),
which replaces (42) in Kneuven et al. (2020) with a weaker version missing the on/off variable.
Equations (45), (43), (44) in Kneuven et al. (2020).
NB: when reading instance, UnitCommitment.jl already calculates difference between max power for segments k and k-1
so the value of cost_segments[k].mw[t] is the max production *for that segment*.
===
Variables
* :segprod
* :is_on
* :prod_above
===
Constraints
* :eq_prod_above_def
* :eq_segprod_limit
"""
function _add_production_piecewise_linear_eqs!(
model::JuMP.Model,
g::Unit,

36
src/model/formulations/DamKucRajAta2016/ramp.jl
View File

@@ -2,6 +2,26 @@
# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details.
"""
_add_ramp_eqs!
Ensure constraints on ramping are met.
Based on Damcı-Kurt et al. (2016).
Eqns. (35), (36) in Kneuven et al. (2020).
Variables
---
* :prod_above
* :reserve
* :is_on
* :switch_on
* :switch_off],
Constraints
---
* :eq_str_ramp_up
* :eq_str_ramp_down
"""
function _add_ramp_eqs!(
model::JuMP.Model,
g::Unit,
@@ -14,12 +34,14 @@ function _add_ramp_eqs!(
RESERVES_WHEN_RAMP_UP = true
RESERVES_WHEN_RAMP_DOWN = true
RESERVES_WHEN_SHUT_DOWN = true
known_initial_conditions = true
is_initially_on = (g.initial_status > 0)
SU = g.startup_limit
SD = g.shutdown_limit
RU = g.ramp_up_limit
RD = g.ramp_down_limit
is_initially_on = _is_initially_on(g)
# The following are the same for generator g across all time periods
SU = g.startup_limit # startup rate
SD = g.shutdown_limit # shutdown rate
RU = g.ramp_up_limit # ramp up rate
RD = g.ramp_down_limit # ramp down rate
gn = g.name
eq_str_ramp_down = _init(model, :eq_str_ramp_down)
eq_str_ramp_up = _init(model, :eq_str_ramp_up)
@@ -94,7 +116,7 @@ function _add_ramp_eqs!(
on_last_period = 0.0
if t > 1
on_last_period = is_on[gn, t-1]
elseif (known_initial_conditions && g.initial_status > 0)
elseif is_initially_on
on_last_period = 1.0
end

22
src/model/formulations/Gar1962/prod.jl
View File

@@ -2,6 +2,11 @@
# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details.
"""
_add_production_vars!(model, unit, formulation_prod_vars)
Creates variables `:prod_above` and `:segprod`.
"""
function _add_production_vars!(
model::JuMP.Model,
g::Unit,
@@ -18,6 +23,23 @@ function _add_production_vars!(
return
end
"""
_add_production_limit_eqs!(model, unit, formulation_prod_vars)
Ensure production limit constraints are met.
Based on Garver (1962) and Morales-España et al. (2013).
Eqns. (18), part of (69) in Kneuven et al. (2020).
===
Variables
* :is_on
* :prod_above
* :reserve
===
Constraints
* :eq_prod_limit
"""
function _add_production_limit_eqs!(
model::JuMP.Model,
g::Unit,

22
src/model/formulations/Gar1962/pwlcosts.jl
View File

@@ -2,6 +2,28 @@
# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details.
"""
_add_production_piecewise_linear_eqs!
Ensure respect of production limits along each segment.
Based on Garver (1962).
Equations (42), (43), (44) in Kneuven et al. (2020).
NB: when reading instance, UnitCommitment.jl already calculates difference between max power for segments k and k-1,
so the value of cost_segments[k].mw[t] is the max production *for that segment*.
Added to `model`: `:eq_prod_above_def` and `:eq_segprod_limit`.
===
Variables
* :segprod
* :is_on
* :prod_above
===
Constraints
* :eq_prod_above_def
* :eq_segprod_limit
"""
function _add_production_piecewise_linear_eqs!(
model::JuMP.Model,
g::Unit,

106
src/model/formulations/Gar1962/status.jl
View File

@@ -2,28 +2,83 @@
# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details.
"""
_add_status_vars!
Create `is_on`, `switch_on`, and `switch_off` variables.
Fix variables if a certain generator _must_ run or based on initial conditions.
"""
function _add_status_vars!(
model::JuMP.Model,
g::Unit,
formulation_status_vars::Gar1962.StatusVars,
ALWAYS_CREATE_VARS = false
)::Nothing
is_on = _init(model, :is_on)
switch_on = _init(model, :switch_on)
switch_off = _init(model, :switch_off)
for t in 1:model[:instance].time
if g.must_run[t]
is_on[g.name, t] = 1.0
switch_on[g.name, t] = (t == 1 ? 1.0 - _is_initially_on(g) : 0.0)
switch_off[g.name, t] = 0.0
else
if ALWAYS_CREATE_VARS || !g.must_run[t]
is_on[g.name, t] = @variable(model, binary = true)
switch_on[g.name, t] = @variable(model, binary = true)
switch_off[g.name, t] = @variable(model, binary = true)
end
if ALWAYS_CREATE_VARS
# If variables are created, use initial conditions to fix some values
if t == 1
if _is_initially_on(g)
# Generator was on (for g.initial_status time periods),
# so cannot be more switched on until the period after the first time it can be turned off
fix(switch_on[g.name, 1], 0.0; force = true)
else
# Generator is initially off (for -g.initial_status time periods)
# Cannot be switched off more
fix(switch_off[g.name, 1], 0.0; force = true)
end
end
if g.must_run[t]
# If the generator _must_ run, then it is obviously on and cannot be switched off
# In the first time period, force unit to switch on if was off before
# Otherwise, unit is on, and will never turn off, so will never need to turn on
fix(is_on[g.name, t], 1.0; force = true)
fix(switch_on[g.name, t], (t == 1 ? 1.0 - _is_initially_on(g) : 0.0); force = true)
fix(switch_off[g.name, t], 0.0; force = true)
end
else
# If vars are not created, then replace them by a constant
if t == 1
if _is_initially_on(g)
switch_on[g.name, t] = 0.0
else
switch_off[g.name, t] = 0.0
end
end
if g.must_run[t]
is_on[g.name, t] = 1.0
switch_on[g.name, t] = (t == 1 ? 1.0 - _is_initially_on(g) : 0.0)
switch_off[g.name, t] = 0.0
end
end # check if ALWAYS_CREATE_VARS
end
return
end
"""
_add_status_eqs!
Variables
---
* is_on
* switch_off
* switch_on
Constraints
---
* eq_binary_link
* eq_switch_on_off
"""
function _add_status_eqs!(
model::JuMP.Model,
g::Unit,
@@ -35,27 +90,32 @@ function _add_status_eqs!(
switch_off = model[:switch_off]
switch_on = model[:switch_on]
for t in 1:model[:instance].time
if !g.must_run[t]
# Link binary variables
if t == 1
eq_binary_link[g.name, t] = @constraint(
model,
is_on[g.name, t] - _is_initially_on(g) ==
switch_on[g.name, t] - switch_off[g.name, t]
)
else
eq_binary_link[g.name, t] = @constraint(
model,
is_on[g.name, t] - is_on[g.name, t-1] ==
switch_on[g.name, t] - switch_off[g.name, t]
)
end
# Cannot switch on and off at the same time
eq_switch_on_off[g.name, t] = @constraint(
if g.must_run[t]
continue
end
# Link binary variables
# Equation (2) in Kneuven et al. (2020), originally from Garver (1962)
if t == 1
eq_binary_link[g.name, t] = @constraint(
model,
switch_on[g.name, t] + switch_off[g.name, t] <= 1
is_on[g.name, t] - _is_initially_on(g) ==
switch_on[g.name, t] - switch_off[g.name, t]
)
else
eq_binary_link[g.name, t] = @constraint(
model,
is_on[g.name, t] - is_on[g.name, t-1] ==
switch_on[g.name, t] - switch_off[g.name, t]
)
end
# Cannot switch on and off at the same time
# amk: I am not sure this is in Kneuven et al. (2020)
eq_switch_on_off[g.name, t] = @constraint(
model,
switch_on[g.name, t] + switch_off[g.name, t] <= 1
)
end
return
end

86
src/model/formulations/GenMorRam2017/startstop.jl Normal file
View File

@@ -0,0 +1,86 @@
# UnitCommitment.jl: Optimization Package for Security-Constrained Unit Commitment
# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details.
"""
_add_startup_shutdown_limit_eqs!(model::JuMP.Model, g::Unit)::Nothing
Startup and shutdown limits from Gentile et al. (2017).
Eqns. (20), (23a), and (23b) in Kneuven et al. (2020).
Variables
---
* :is_on
* :prod_above
* :reserve
* :switch_on
* :switch_off
Constraints
---
* :eq_startstop_limit
* :eq_startup_limit
* :eq_shutdown_limit
"""
function _add_startup_shutdown_limit_eqs!(model::JuMP.Model, g::Unit)::Nothing
# TODO: Move upper case constants to model[:instance]
RESERVES_WHEN_START_UP = true
RESERVES_WHEN_RAMP_UP = true
RESERVES_WHEN_RAMP_DOWN = true
RESERVES_WHEN_SHUT_DOWN = true
eq_startstop_limit = _init(model, :eq_startstop_limit)
eq_shutdown_limit = _init(model, :eq_shutdown_limit)
eq_startup_limit = _init(model, :eq_startup_limit)
is_on = model[:is_on]
prod_above = model[:prod_above]
reserve = model[:reserve]
switch_off = model[:switch_off]
switch_on = model[:switch_on]
T = model[:instance].time
gi = g.name
if g.initial_power > g.shutdown_limit
#eqs.shutdown_limit[gi, 0] = @constraint(mip, vars.switch_off[gi, 1] <= 0)
fix(switch_off[gi, 1], 0.; force = true)
end
for t = 1:T
## 2020-10-09 amk: added eqn (20) and check of g.min_uptime
# Not present in (23) in Kneueven et al.
if g.min_uptime > 1
# Equation (20) in Kneuven et al. (2020)
eqs.startstop_limit[gi,t] =
@constraint(model,
prod_above[gi, t] + reserve[gi, t]
<= (g.max_power[t] - g.min_power[t]) * is_on[gi, t]
- max(0, g.max_power[t] - g.startup_limit) * switch_on[gi, t]
- (t < T ? max(0, g.max_power[t] - g.shutdown_limit) * switch_off[gi, t+1] : 0.)
)
else
## Startup limits
# Equation (23a) in Kneuven et al. (2020)
eqs.startup_limit[gi, t] =
@constraint(model,
prod_above[gi, t] + reserve[gi, t]
<= (g.max_power[t] - g.min_power[t]) * is_on[gi, t]
- max(0, g.max_power[t] - g.startup_limit) * switch_on[gi, t]
- (t < T ? max(0, g.startup_limit - g.shutdown_limit) * switch_off[gi, t+1] : 0.)
)
## Shutdown limits
if t < T
# Equation (23b) in Kneuven et al. (2020)
eqs.shutdown_limit[gi, t] =
@constraint(model,
prod_above[gi, t] + reserve[gi, t]
<= (g.max_power[t] - g.min_power[t]) * xis_on[gi, t]
- (t < T ? max(0, g.max_power[t] - g.shutdown_limit) * switch_off[gi, t+1] : 0.)
- max(0, g.shutdown_limit - g.startup_limit) * switch_on[gi, t])
end
end # check if g.min_uptime > 1
end # loop over time
end # _add_startup_shutdown_limit_eqs!

24
src/model/formulations/KnuOstWat2018/pwlcosts.jl
View File

@@ -2,6 +2,30 @@
# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details.
"""
_add_production_piecewise_linear_eqs!
Ensure respect of production limits along each segment.
Based on Kneuven et al. (2018b).
Eqns. (43), (44), (46), (48) in Kneuven et al. (2020).
NB: when reading instance, UnitCommitment.jl already calculates difference between max power for segments k and k-1
so the value of cost_segments[k].mw[t] is the max production *for that segment*.
===
Variables
* :segprod
* :is_on
* :switch_on
* :switch_off
* :prod_above
===
Constraints
* :eq_prod_above_def
* :eq_segprod_limit_a
* :eq_segprod_limit_b
* :eq_segprod_limit_c
"""
function _add_production_piecewise_linear_eqs!(
model::JuMP.Model,
g::Unit,

116
src/model/formulations/KnuOstWat2018/scosts.jl Normal file
View File

@@ -0,0 +1,116 @@
# UnitCommitment.jl: Optimization Package for Security-Constrained Unit Commitment
# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details.
"""
_add_startup_cost_eqs!
Extended formulation of startup costs using indicator variables
based on Kneuven, Ostrowski, and Watson, 2020
--- equations (59), (60), (61).
Variables
---
* switch_on
* switch_off
* downtime_arc
Constraints
---
* eq_startup_at_t
* eq_shutdown_at_t
"""
function _add_startup_cost_eqs!(
model::JuMP.Model,
g::Unit,
formulation::MorLatRam2013.StartupCosts,
)::Nothing
S = length(g.startup_categories)
if S == 0
return
end
gn = g.name
_init(model, eq_startup_at_t)
_init(model, eq_shutdown_at_t)
switch_on = model[:switch_on]
switch_off = model[:switch_off]
downtime_arc = model[:downtime_arc]
DT = g.min_downtime # minimum time offline
TC = g.startup_categories[S].delay # time offline until totally cold
# If initial_status < 0, then this is the amount of time the generator has been off
initial_time_shutdown = (g.initial_status < 0 ? -g.initial_status : 0)
for t in 1:model[:instance].time
# Fix to zero values of downtime_arc outside the feasible time pairs
# Specifically, x(t,t') = 0 if t' does not belong to 𝒢 = [t+DT, t+TC-1]
# This is because DT is the minimum downtime, so there is no way x(t,t')=1 for t'<t+DT
# and TC is the "time until cold" => if the generator starts afterwards, always has max cost
#start_time = min(t + DT, T)
#end_time = min(t + TC - 1, T)
#for tmp_t in t+1:start_time
# fix(vars.downtime_arc[gn, t, tmp_t], 0.; force = true)
#end
#for tmp_t in end_time+1:T
# fix(vars.downtime_arc[gn, t, tmp_t], 0.; force = true)
#end
# Equation (59) in Kneuven et al. (2020)
# Relate downtime_arc with switch_on
# "switch_on[g,t] >= x_g(t',t) for all t' \in [t-TC+1, t-DT]"
eq_startup_at_t[gn, t] =
@constraint(model,
switch_on[gn, t]
>= sum(downtime_arc[gn,tmp_t,t]
for tmp_t in t-TC+1:t-DT if tmp_t >= 1)
)
# Equation (60) in Kneuven et al. (2020)
# "switch_off[g,t] >= x_g(t,t') for all t' \in [t+DT, t+TC-1]"
eqs.shutdown_at_t[gn, t] =
@constraint(model,
switch_off[gn, t]
>= sum(downtime_arc[gn,t,tmp_t]
for tmp_t in t+DT:t+TC-1 if tmp_t <= T)
)
# Objective function terms for start-up costs
# Equation (61) in Kneuven et al. (2020)
default_category = S
if initial_time_shutdown > 0 && t + initial_time_shutdown - 1 < TC
for s in 1:S-1
# If off for x periods before, then belongs to category s
# if -x+1 in [t-delay[s+1]+1,t-delay[s]]
# or, equivalently, if total time off in [delay[s], delay[s+1]-1]
# where total time off = t - 1 + initial_time_shutdown
# (the -1 because not off for current time period)
if t + initial_time_shutdown - 1 < g.startup_categories[s+1].delay
default_category = s
break # does not go into next category
end
end
end
add_to_expression!(model[:obj],
switch_on[gn, t],
g.startup_categories[default_category].cost)
for s in 1:S-1
# Objective function terms for start-up costs
# Equation (61) in Kneuven et al. (2020)
# Says to replace the cost of last category with cost of category s
start_range = max((t - g.startup_categories[s + 1].delay + 1),1)
end_range = min((t - g.startup_categories[s].delay),T-1)
for tmp_t in start_range:end_range
if (t < tmp_t + DT) || (t >= tmp_t + TC) # the second clause should never be true for s < S
continue
end
add_to_expression!(model[:obj],
downtime_arc[gn,tmp_t,t],
g.startup_categories[s].cost - g.startup_categories[S].cost)
end
end # iterate over startup categories
end # iterate over time
end # add_startup_costs_KneOstWat20

29
src/model/formulations/MorLatRam2013/ramp.jl
View File

@@ -2,6 +2,27 @@
# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details.
"""
_add_ramp_eqs!
Ensure constraints on ramping are met.
Needs to be used in combination with shutdown rate constraints, e.g., (21b) in Kneuven et al. (2020).
Based on Morales-España, Latorre, and Ramos, 2013.
Eqns. (26)+(27) [replaced by (24)+(25) if time-varying min demand] in Kneuven et al. (2020).
Variables
---
* :is_on
* :switch_off
* :switch_on
* :prod_above
* :reserve
Constraints
---
* :eq_ramp_up
* :eq_ramp_down
"""
function _add_ramp_eqs!(
model::JuMP.Model,
g::Unit,
@@ -15,10 +36,10 @@ function _add_ramp_eqs!(
RESERVES_WHEN_RAMP_DOWN = true
RESERVES_WHEN_SHUT_DOWN = true
is_initially_on = (g.initial_status > 0)
SU = g.startup_limit
SD = g.shutdown_limit
RU = g.ramp_up_limit
RD = g.ramp_down_limit
SU = g.startup_limit # startup rate
SD = g.shutdown_limit # shutdown rate
RU = g.ramp_up_limit # ramp up rate
RD = g.ramp_down_limit # ramp down rate
gn = g.name
eq_ramp_down = _init(model, :eq_ramp_down)
eq_ramp_up = _init(model, :eq_str_ramp_up)

64
src/model/formulations/MorLatRam2013/scosts.jl
View File

@@ -2,21 +2,57 @@
# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details.
"""
_add_startup_cost_eqs!
Extended formulation of startup costs using indicator variables
based on Muckstadt and Wilson, 1968;
this version by Morales-España, Latorre, and Ramos, 2013.
Eqns. (54), (55), and (56) in Kneuven et al. (2020).
Note that the last 'constraint' is actually setting the objective.
\tstartup[gi,s,t] ≤ sum_{i=s.delay}^{(s+1).delay-1} switch_off[gi,t-i]
\tswitch_on[gi,t] = sum_{s=1}^{length(startup_categories)} startup[gi,s,t]
\tstartup_cost[gi,t] = sum_{s=1}^{length(startup_categories)} cost_segments[s].cost * startup[gi,s,t]
Variables
---
* startup
* switch_on
* switch_off
Constraints
---
* eq_startup_choose
* eq_startup_restrict
"""
function _add_startup_cost_eqs!(
model::JuMP.Model,
g::Unit,
formulation::MorLatRam2013.StartupCosts,
)::Nothing
S = length(g.startup_categories)
if S == 0
return
end
# Constraints created
eq_startup_choose = _init(model, :eq_startup_choose)
eq_startup_restrict = _init(model, :eq_startup_restrict)
S = length(g.startup_categories)
# Variables needed
startup = model[:startup]
switch_on = model[:switch_on]
switch_off = model[:switch_off]
gn = g.name
for t in 1:model[:instance].time
# If unit is switching on, we must choose a startup category
eq_startup_choose[g.name, t] = @constraint(
# Equation (55) in Kneuven et al. (2020)
eq_startup_choose[gn, t] = @constraint(
model,
model[:switch_on][g.name, t] ==
sum(startup[g.name, t, s] for s in 1:S)
switch_on[gn, t] ==
sum(startup[gn, t, s] for s in 1:S)
)
for s in 1:S
@@ -26,25 +62,27 @@ function _add_startup_cost_eqs!(
range_start = t - g.startup_categories[s+1].delay + 1
range_end = t - g.startup_categories[s].delay
range = (range_start:range_end)
# If initial_status < 0, then this is the amount of time the generator has been off
initial_sum = (
g.initial_status < 0 && (g.initial_status + 1 in range) ? 1.0 : 0.0
)
eq_startup_restrict[g.name, t, s] = @constraint(
# Change of index version of equation (54) in Kneuven et al. (2020):
# startup[gi,s,t] ≤ sum_{i=s.delay}^{(s+1).delay-1} switch_off[gi,t-i]
eq_startup_restrict[gn, t, s] = @constraint(
model,
startup[g.name, t, s] <=
initial_sum + sum(
model[:switch_off][g.name, i] for i in range if i >= 1
)
startup[gn, t, s] <=
initial_sum + sum(switch_off[gn, i] for i in range if i >= 1)
)
end
end # if s < S (not the last category)
# Objective function terms for start-up costs
# Equation (56) in Kneuven et al. (2020)
add_to_expression!(
model[:obj],
startup[g.name, t, s],
startup[gn, t, s],
g.startup_categories[s].cost,
)
end
end
end # iterate over startup categories
end # iterate over time
return
end

86
src/model/formulations/MorLatRam2013/startstop.jl Normal file
View File

@@ -0,0 +1,86 @@
# UnitCommitment.jl: Optimization Package for Security-Constrained Unit Commitment
# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details.
"""
_add_startup_shutdown_limit_eqs!(model::JuMP.Model, g::Unit)::Nothing
Startup and shutdown limits from Morales-España et al. (2013a).
Eqns. (20), (21a), and (21b) in Kneuven et al. (2020).
Variables
---
* :is_on
* :prod_above
* :reserve
* :switch_on
* :switch_off
Constraints
---
* :eq_startstop_limit
* :eq_startup_limit
* :eq_shutdown_limit
"""
function _add_startup_shutdown_limit_eqs!(model::JuMP.Model, g::Unit)::Nothing
# TODO: Move upper case constants to model[:instance]
RESERVES_WHEN_START_UP = true
RESERVES_WHEN_RAMP_UP = true
RESERVES_WHEN_RAMP_DOWN = true
RESERVES_WHEN_SHUT_DOWN = true
eq_startstop_limit = _init(model, :eq_startstop_limit)
eq_shutdown_limit = _init(model, :eq_shutdown_limit)
eq_startup_limit = _init(model, :eq_startup_limit)
is_on = model[:is_on]
prod_above = model[:prod_above]
reserve = model[:reserve]
switch_off = model[:switch_off]
switch_on = model[:switch_on]
T = model[:instance].time
gi = g.name
for t = 1:T
## 2020-10-09 amk: added eqn (20) and check of g.min_uptime
if g.min_uptime > 1 && t < T
# Equation (20) in Kneuven et al. (2020)
# UT > 1 required, to guarantee that vars.switch_on[gi, t] and vars.switch_off[gi, t+1] are not both = 1 at the same time
eq_startstop_limit[gi,t] =
@constraint(model,
prod_above[gi, t] + reserve[gi, t]
<= (g.max_power[t] - g.min_power[t]) * is_on[gi, t]
- max(0, g.max_power[t] - g.startup_limit) * switch_on[gi, t]
- max(0, g.max_power[t] - g.shutdown_limit) * switch_off[gi, t+1])
else
## Startup limits
# Equation (21a) in Kneuven et al. (2020)
# Proposed by Morales-España et al. (2013a)
eqs_startup_limit[gi, t] =
@constraint(model,
prod_above[gi, t] + reserve[gi, t]
<= (g.max_power[t] - g.min_power[t]) * is_on[gi, t]
- max(0, g.max_power[t] - g.startup_limit) * switch_on[gi, t])
## Shutdown limits
if t < T
# Equation (21b) in Kneuven et al. (2020)
# TODO different from what was in previous model, due to reserve variable
# ax: ideally should have reserve_up and reserve_down variables
# i.e., the generator should be able to increase/decrease production as specified
# (this is a heuristic for a "robust" solution,
# in case there is an outage or a surge, and flow has to be redirected)
# amk: if shutdown_limit is the max prod of generator in time period before shutting down,
# then it makes sense to count reserves, because otherwise, if reserves ≠ 0,
# then the generator will actually produce more than the limit
eqs.shutdown_limit[gi, t] =
@constraint(model,
prod_above[gi, t]
+ (RESERVES_WHEN_SHUT_DOWN ? reserve[gi, t] : 0.) # amk added
<= (g.max_power[t] - g.min_power[t]) * is_on[gi, t]
- max(0, g.max_power[t] - g.shutdown_limit) * switch_off[gi, t+1])
end
end # check if g.min_uptime > 1
end # loop over time
end # _add_startup_shutdown_limit_eqs!

20
src/model/formulations/NowRom2000/scosts.jl Normal file
View File

@@ -0,0 +1,20 @@
# UnitCommitment.jl: Optimization Package for Security-Constrained Unit Commitment
# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details.
"""
_add_startup_cost_eqs!
Based on Nowak and Römisch, 2000.
Introduces auxiliary startup cost variable, c_g^SU(t) for each time period,
and uses startup status variable, u_g(t);
there are exponentially many facets in this space,
but there is a linear-time separation algorithm (Brandenburg et al., 2017).
"""
function _add_startup_cost_eqs!(
model::JuMP.Model,
g::Unit,
formulation::MorLatRam2013.StartupCosts,
)::Nothing
error("Not implemented.")
end

83
src/model/formulations/OstAnjVan2012/ramp.jl Normal file
View File

@@ -0,0 +1,83 @@
# UnitCommitment.jl: Optimization Package for Security-Constrained Unit Commitment
# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details.
"""
_add_ramp_eqs!
Ensure constraints on ramping are met.
Based on Ostrowski, Anjos, Vannelli (2012).
Eqn (37) in Kneuven et al. (2020).
Variables
---
* :is_on
* :prod_above
* :reserve
Constraints
---
* :eq_str_prod_limit
"""
function _add_ramp_eqs!(
model::JuMP.Model,
g::Unit,
formulation_prod_vars::Gar1962.ProdVars,
formulation_ramping::MorLatRam2013.Ramping,
formulation_status_vars::Gar1962.StatusVars,
)::Nothing
# TODO: Move upper case constants to model[:instance]
RESERVES_WHEN_START_UP = true
RESERVES_WHEN_RAMP_UP = true
RESERVES_WHEN_RAMP_DOWN = true
RESERVES_WHEN_SHUT_DOWN = true
is_initially_on = _is_initially_on(g)
gn = g.name
eq_str_prod_limit = _init(model, :eq_str_prod_limit)
# Variables that we need
reserve = model[:reserve]
# Gar1962.ProdVars
prod_above = model[:prod_above]
# Gar1962.StatusVars
is_on = model[:is_on]
switch_off = model[:switch_off]
# The following are the same for generator g across all time periods
UT = g.min_uptime
SU = g.startup_limit # startup rate
SD = g.shutdown_limit # shutdown rate
RU = g.ramp_up_limit # ramp up rate
RD = g.ramp_down_limit # ramp down rate
# TODO check initial conditions, but maybe okay as long as (35) and (36) are also used
for t in 1:model[:instance].time
Pbar = g.max_power[t]
#TRD = floor((Pbar - SU)/RD)
# TODO check amk changed TRD wrt Kneuven et al.
TRD = ceil((Pbar - SD) / RD) # ramp down time
if Pbar < 1e-7
# Skip this time period if max power = 0
continue
end
if UT >= 1
# Equation (37) in Kneuven et al. (2020)
KSD = min( TRD, UT-1, T-t-1 )
eq_str_prod_limit[gn, t] =
@constraint(model,
prod_above[gn, t] + g.min_power[t] * is_on[gn, t]
+ (RESERVES_WHEN_RAMP_DOWN ? reserve[gn, t] : 0.) # amk added; TODO: should this be RESERVES_WHEN_RAMP_DOWN or RESERVES_WHEN_SHUT_DOWN?
<= Pbar * is_on[gi, t]
- sum((Pbar - (SD + i * RD)) * switch_off[gi, t+1+i]
for i in 0:KSD)
)
end # check UT >= 1
end # loop over time
end

22
src/model/formulations/PanGua2016/ramp.jl
View File

@@ -2,6 +2,28 @@
# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details.
"""
_add_ramp_eqs!
Add tighter upper bounds on production based on ramp-down trajectory.
Based on (28) in Pan and Guan (2016).
But there is an extra time period covered using (40) of Kneuven et al. (2020).
Eqns. (38), (40), (41) in Kneuven et al. (2020).
Variables
---
* :prod_above
* :reserve
* :is_on
* :switch_on
* :switch_off
Constraints
---
* :str_prod_limit
* :prod_limit_ramp_up_extra_period
* :prod_limit_shutdown_trajectory
"""
function _add_ramp_eqs!(
model::JuMP.Model,
g::Unit,

1
src/model/formulations/base/bus.jl
View File

@@ -4,7 +4,6 @@
function _add_bus!(model::JuMP.Model, b::Bus)::Nothing
net_injection = _init(model, :expr_net_injection)
reserve = _init(model, :expr_reserve)
curtail = _init(model, :curtail)
for t in 1:model[:instance].time
# Fixed load

3
src/model/formulations/base/psload.jl
View File

@@ -2,6 +2,9 @@
# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details.
"""
_add_price_sensitive_load!(model::JuMP.Model, ps::PriceSensitiveLoad)
"""
function _add_price_sensitive_load!(
model::JuMP.Model,
ps::PriceSensitiveLoad,

2
src/model/formulations/base/sensitivity.jl
View File

@@ -18,7 +18,7 @@ function _injection_shift_factors(;
lines::Array{TransmissionLine},
)
susceptance = _susceptance_matrix(lines)
incidence = _reduced_incidence_matrix(lines = lines, buses = buses)
incidence = _reduced_incidence_matrix(buses = buses, lines = lines)
laplacian = transpose(incidence) * susceptance * incidence
isf = susceptance * incidence * inv(Array(laplacian))
return isf

29
src/model/formulations/base/structs.jl
View File

@@ -2,28 +2,37 @@
# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details.
abstract type TransmissionFormulation end
abstract type RampingFormulation end
abstract type AbstractTransmissionFormulation end
abstract type AbstractRampingFormulation end
abstract type PiecewiseLinearCostsFormulation end
abstract type StartupCostsFormulation end
abstract type AbstractStartupCostsFormulation end
abstract type StatusVarsFormulation end
abstract type ProductionVarsFormulation end
"""
struct Formulation
Every formulation has to specify components for setting production variables and limits,
startup costs and piecewise-linear costs, ramping variables and constraints,
status variables (on/off, starting up, shutting down), and transmission constraints.
Some of these components are allowed to be empty, as long as overall validity of the formulation is maintained.
"""
struct Formulation
prod_vars::ProductionVarsFormulation
pwl_costs::PiecewiseLinearCostsFormulation
ramping::RampingFormulation
startup_costs::StartupCostsFormulation
ramping::AbstractRampingFormulation
startup_costs::AbstractStartupCostsFormulation
status_vars::StatusVarsFormulation
transmission::TransmissionFormulation
transmission::AbstractTransmissionFormulation
function Formulation(;
prod_vars::ProductionVarsFormulation = Gar1962.ProdVars(),
pwl_costs::PiecewiseLinearCostsFormulation = KnuOstWat2018.PwlCosts(),
ramping::RampingFormulation = MorLatRam2013.Ramping(),
startup_costs::StartupCostsFormulation = MorLatRam2013.StartupCosts(),
ramping::AbstractRampingFormulation = MorLatRam2013.Ramping(),
startup_costs::AbstractStartupCostsFormulation = MorLatRam2013.StartupCosts(),
status_vars::StatusVarsFormulation = Gar1962.StatusVars(),
transmission::TransmissionFormulation = ShiftFactorsFormulation(),
transmission::AbstractTransmissionFormulation = ShiftFactorsFormulation(),
)
return new(
prod_vars,
@@ -61,7 +70,7 @@ Arguments
the cutoff that should be applied to the LODF matrix. Entries with magnitude
smaller than this value will be set to zero.
"""
struct ShiftFactorsFormulation <: TransmissionFormulation
struct ShiftFactorsFormulation <: AbstractTransmissionFormulation
isf_cutoff::Float64
lodf_cutoff::Float64
precomputed_isf::Union{Nothing,Matrix{Float64}}

64
src/model/formulations/base/system.jl
View File

@@ -2,12 +2,32 @@
# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details.
"""
_add_system_wide_eqs!(model::JuMP.Model)::Nothing
Calls `_add_net_injection_eqs!` and `add_reserve_eqs!`.
"""
function _add_system_wide_eqs!(model::JuMP.Model)::Nothing
_add_net_injection_eqs!(model)
_add_reserve_eqs!(model)
return
end
"""
_add_net_injection_eqs!(model::JuMP.Model)::Nothing
Adds `net_injection`, `eq_net_injection_def`, and `eq_power_balance` identifiers into `model`.
Variables
---
* expr_net_injection
* net_injection
Constraints
---
* eq_net_injection_def
* eq_power_balance
"""
function _add_net_injection_eqs!(model::JuMP.Model)::Nothing
T = model[:instance].time
net_injection = _init(model, :net_injection)
@@ -27,15 +47,49 @@ function _add_net_injection_eqs!(model::JuMP.Model)::Nothing
return
end
"""
_add_reserve_eqs!(model::JuMP.Model)::Nothing
Ensure constraints on reserves are met.
Based on Morales-España et al. (2013a).
Eqn. (68) from Kneuven et al. (2020).
Adds `eq_min_reserve` identifier to `model`, and corresponding constraint.
Variables
---
* reserve
* reserve_shortfall
Constraints
---
* eq_min_reserve
"""
function _add_reserve_eqs!(model::JuMP.Model)::Nothing
instance = model[:instance]
eq_min_reserve = _init(model, :eq_min_reserve)
for t in 1:model[:instance].time
for t in 1:instance.time
# Equation (68) in Kneuven et al. (2020)
# As in Morales-España et al. (2013a)
# Akin to the alternative formulation with max_power_avail
# from Carrión and Arroyo (2006) and Ostrowski et al. (2012)
shortfall_penalty = instance.shortfall_penalty[t]
eq_min_reserve[t] = @constraint(
model,
sum(
model[:expr_reserve][b.name, t] for b in model[:instance].buses
) >= model[:instance].reserves.spinning[t]
sum(model[:reserve][g.name, t] for g in instance.units)
+ (shortfall_penalty > 1e-7 ? model[:reserve_shortfall][t] : 0.)
>= instance.reserves.spinning[t]
)
end
# Account for shortfall contribution to objective
if shortfall_penalty > 1e-7
add_to_expression!(model.obj,
shortfall_penalty,
model[:reserve_shortfall][t])
else
# Not added to the model at all
#fix(model.vars.reserve_shortfall[t], 0.; force=true)
end
end # loop over time
return
end

134
src/model/formulations/base/unit.jl
View File

@@ -2,6 +2,15 @@
# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details.
"""
_add_unit!(model::JuMP.Model, g::Unit, formulation::Formulation)
Add production, reserve, startup, shutdown, and status variables,
and constraints for min uptime/downtime, net injection, production, ramping, startup, shutdown, and status.
Fix variables if a certain generator _must_ run or if a generator provides spinning reserves.
Also, add overflow penalty to objective for each transmission line.
"""
function _add_unit!(model::JuMP.Model, g::Unit, formulation::Formulation)
if !all(g.must_run) && any(g.must_run)
error("Partially must-run units are not currently supported")
@@ -35,6 +44,7 @@ function _add_unit!(model::JuMP.Model, g::Unit, formulation::Formulation)
formulation.status_vars,
)
_add_startup_cost_eqs!(model, g, formulation.startup_costs)
_add_shutdown_cost_eqs!(model, g)
_add_startup_shutdown_limit_eqs!(model, g)
_add_status_eqs!(model, g, formulation.status_vars)
return
@@ -42,26 +52,42 @@ end
_is_initially_on(g::Unit)::Float64 = (g.initial_status > 0 ? 1.0 : 0.0)
function _add_reserve_vars!(model::JuMP.Model, g::Unit)::Nothing
"""
_add_reserve_vars!(model::JuMP.Model, g::Unit)::Nothing
Add `:reserve` variable to `model`, fixed to zero if no spinning reserves specified.
"""
function _add_reserve_vars!(model::JuMP.Model, g::Unit, ALWAYS_CREATE_VARS = false)::Nothing
reserve = _init(model, :reserve)
reserve_shortfall = _init(model, :reserve_shortfall) # for accounting for shortfall penalty in the objective
for t in 1:model[:instance].time
if g.provides_spinning_reserves[t]
reserve[g.name, t] = @variable(model, lower_bound = 0)
else
reserve[g.name, t] = 0.0
if ALWAYS_CREATE_VARS
reserve[g.name, t] = @variable(model, lower_bound = 0)
fix(reserve[g.name, t], 0.0; force = true)
else
reserve[g.name, t] = 0.0
end
end
end
return
end
"""
_add_reserve_eqs!(model::JuMP.Model, g::Unit)::Nothing
"""
function _add_reserve_eqs!(model::JuMP.Model, g::Unit)::Nothing
reserve = model[:reserve]
for t in 1:model[:instance].time
add_to_expression!(expr_reserve[g.bus.name, t], reserve[g.name, t], 1.0)
end
# nothing to do here
return
end
"""
_add_startup_shutdown_vars!(model::JuMP.Model, g::Unit)::Nothing
Add `startup` to model.
"""
function _add_startup_shutdown_vars!(model::JuMP.Model, g::Unit)::Nothing
startup = _init(model, :startup)
for t in 1:model[:instance].time
@@ -72,6 +98,22 @@ function _add_startup_shutdown_vars!(model::JuMP.Model, g::Unit)::Nothing
return
end
"""
_add_startup_shutdown_limit_eqs!(model::JuMP.Model, g::Unit)::Nothing
Variables
---
* :is_on
* :prod_above
* :reserve
* :switch_on
* :switch_off
Constraints
---
* :eq_startup_limit
* :eq_shutdown_limit
"""
function _add_startup_shutdown_limit_eqs!(model::JuMP.Model, g::Unit)::Nothing
eq_shutdown_limit = _init(model, :eq_shutdown_limit)
eq_startup_limit = _init(model, :eq_startup_limit)
@@ -91,8 +133,13 @@ function _add_startup_shutdown_limit_eqs!(model::JuMP.Model, g::Unit)::Nothing
)
# Shutdown limit
if g.initial_power > g.shutdown_limit
eq_shutdown_limit[g.name, 0] =
@constraint(model, switch_off[g.name, 1] <= 0)
# TODO check what happens with these variables when exporting the model
# Generator producing too much to be turned off in the first time period
# (can a binary variable have bounds x = 0?)
#eqs.shutdown_limit[gi, 0] = @constraint(mip, vars.switch_off[gi, 1] <= 0)
fix(model.vars.switch_off[gi, 1], 0.; force = true)
#eq_shutdown_limit[g.name, 0] =
#@constraint(model, switch_off[g.name, 1] <= 0)
end
if t < T
eq_shutdown_limit[g.name, t] = @constraint(
@@ -107,10 +154,36 @@ function _add_startup_shutdown_limit_eqs!(model::JuMP.Model, g::Unit)::Nothing
return
end
"""
_add_shutdown_cost_eqs!
Variables
---
* :switch_off
"""
function _add_shutdown_cost_eqs!(model::JuMP.Modle, g::Unit)::Nothing
T = model[:instance].time
gi = g.name
for t = 1:T
shutdown_cost = 0.
if shutdown_cost > 1e-7
# Equation (62) in Kneuven et al. (2020)
add_to_expression!(model[:obj],
model[:switch_off][gi, t],
shutdown_cost)
end
end # loop over time
end # _add_shutdown_cost_eqs!
end
"""
_add_ramp_eqs!(model, unit, formulation)
"""
function _add_ramp_eqs!(
model::JuMP.Model,
g::Unit,
formulation::RampingFormulation,
formulation::AbstractRampingFormulation,
)::Nothing
prod_above = model[:prod_above]
reserve = model[:reserve]
@@ -153,6 +226,26 @@ function _add_ramp_eqs!(
end
end
"""
_add_min_uptime_downtime_eqs!(model::JuMP.Model, g::Unit)::Nothing
Ensure constraints on up/down time are met.
Based on Garver (1962), Malkin (2003), and Rajan and Takritti (2005).
Eqns. (3), (4), (5) in Kneuven et al. (2020).
Variables
---
* :is_on
* :switch_off
* :switch_on
Constraints
---
* :eq_min_uptime
* :eq_min_downtime
"""
function _add_min_uptime_downtime_eqs!(model::JuMP.Model, g::Unit)::Nothing
is_on = model[:is_on]
switch_off = model[:switch_off]
@@ -162,18 +255,24 @@ function _add_min_uptime_downtime_eqs!(model::JuMP.Model, g::Unit)::Nothing
T = model[:instance].time
for t in 1:T
# Minimum up-time
# Equation (4) in Kneuven et al. (2020)
eq_min_uptime[g.name, t] = @constraint(
model,
sum(switch_on[g.name, i] for i in (t-g.min_uptime+1):t if i >= 1) <= is_on[g.name, t]
sum(switch_on[g.name, i] for i in (t-g.min_uptime+1):t if i >= 1)
<= is_on[g.name, t]
)
# Minimum down-time
# Equation (5) in Kneuven et al. (2020)
eq_min_downtime[g.name, t] = @constraint(
model,
sum(
switch_off[g.name, i] for i in (t-g.min_downtime+1):t if i >= 1
) <= 1 - is_on[g.name, t]
sum(switch_off[g.name, i] for i in (t-g.min_downtime+1):t if i >= 1)
<= 1 - is_on[g.name, t]
)
# Minimum up/down-time for initial periods
# Equations (3a) and (3b) in Kneuven et al. (2020)
# (using :switch_on and :switch_off instead of :is_on)
if t == 1
if g.initial_status > 0
eq_min_uptime[g.name, 0] = @constraint(
@@ -196,6 +295,9 @@ function _add_min_uptime_downtime_eqs!(model::JuMP.Model, g::Unit)::Nothing
end
end
"""
_add_net_injection_eqs!(model::JuMP.Model, g::Unit)::Nothing
"""
function _add_net_injection_eqs!(model::JuMP.Model, g::Unit)::Nothing
expr_net_injection = model[:expr_net_injection]
for t in 1:model[:instance].time
@@ -210,11 +312,5 @@ function _add_net_injection_eqs!(model::JuMP.Model, g::Unit)::Nothing
model[:is_on][g.name, t],
g.min_power[t],
)
# Add to reserves expression
add_to_expression!(
model[:expr_reserve][g.bus.name, t],
model[:reserve][g.name, t],
1.0,
)
end
end