Fix bug in validation script; create large tests

.github/workflows/test.yml

@@ -9,8 +9,8 @@ jobs:
     runs-on: ${{ matrix.os }}
     strategy:
       matrix:
-        julia-version: ['1.4', '1.5', '1.6']
-        julia-arch: [x64]
+        julia-version: ['1.3', '1.4', '1.5', '1.6']
+        julia-arch: [x64, x86]
         os: [ubuntu-latest, windows-latest, macOS-latest]
         exclude:
           - os: macOS-latest

CHANGELOG.md

@@ -11,10 +11,9 @@ All notable changes to this project will be documented in this file.
 [semver]: https://semver.org/spec/v2.0.0.html
 [pkjjl]: https://pkgdocs.julialang.org/v1/compatibility/#compat-pre-1.0
 
-## [0.2.2] - 2021-07-21
+## [0.2.2] - Unreleased
 ### Fixed
 - Fix small bug in validation scripts related to startup costs
-- Fix duplicated startup constraints (@mtanneau, #12)
 
 ## [0.2.1] - 2021-06-02
 ### Added

docs/format.md

@@ -28,14 +28,13 @@ Each section is described in detail below. For a complete example, see [case14](
 
 ### Parameters
 
-This section describes system-wide parameters, such as power balance and reserve shortfall penalties, and optimization parameters, such as the length of the planning horizon and the time.
+This section describes system-wide parameters, such as power balance penalties,  optimization parameters, such as the length of the planning horizon and the time.
 
 | Key                            | Description                                       | Default  | Time series?
 | :----------------------------- | :------------------------------------------------ | :------: | :------------:
-| `Time horizon (h)` | Length of the planning horizon (in hours). | Required | N
+| `Time horizon (h)`                     | Length of the planning horizon (in hours). | Required | N
 | `Time step (min)` | Length of each time step (in minutes). Must be a divisor of 60 (e.g. 60, 30, 20, 15, etc). | `60` | N
 | `Power balance penalty ($/MW)` | Penalty for system-wide shortage or surplus in production (in $/MW). This is charged per time step. For example, if there is a shortage of 1 MW for three time steps, three times this amount will be charged. | `1000.0` | Y
-| `Reserve shortfall penalty ($/MW)` | Penalty for system-wide shortage in meeting reserve requirements (in $/MW). This is charged per time step. Negative value implies reserve constraints must always be satisfied. | `-1` | Y
 
 
 #### Example
@@ -43,8 +42,7 @@ This section describes system-wide parameters, such as power balance and reserve
 {
     "Parameters": {
         "Time horizon (h)": 4,
-        "Power balance penalty ($/MW)": 1000.0,
-        "Reserve shortfall penalty ($/MW)": -1.0
+        "Power balance penalty ($/MW)": 1000.0
     }
 }
 ```

docs/model.md

@@ -148,7 +148,7 @@ for g in instance.units
 end
 ```
 
-### Fixing variables, modifying objective function and adding constraints
+### Modifying the model
 
 Since we now have a direct reference to the JuMP decision variables, it is possible to fix variables, change the coefficients in the objective function, or even add new constraints to the model before solving it. The script below shows how can this be accomplished. For more information on modifying an existing model, [see the JuMP documentation](https://jump.dev/JuMP.jl/stable/manual/variables/).
 
@@ -190,54 +190,6 @@ JuMP.set_objective_coefficient(
 UnitCommitment.optimize!(model)
 ```
 
-### Adding new component to a bus
-
-The following snippet shows how to add a new grid component to a particular bus. For each time step, we create decision variables for the new grid component, add these variables to the objective function, then attach the component to a particular bus by modifying some existing model constraints. 
-
-```julia
-using Cbc
-using JuMP
-using UnitCommitment
-
-# Load instance and build base model
-instance = UnitCommitment.read_benchmark("matpower/case118/2017-02-01")
-model = UnitCommitment.build_model(
-    instance=instance,
-    optimizer=Cbc.Optimizer,
-)
-
-# Get the number of time steps in the original instance
-T = instance.time
-
-# Create decision variables for the new grid component.
-# In this example, we assume that the new component can
-# inject up to 10 MW of power at each time step, so we
-# create new continuous variables 0 ≤ x[t] ≤ 10.
-@variable(model, x[1:T], lower_bound=0.0, upper_bound=10.0)
-
-# For each time step
-for t in 1:T
-
-    # Add production costs to the objective function.
-    # In this example, we assume a cost of $5/MW.
-    set_objective_coefficient(model, x[t], 5.0)
-
-    # Attach the new component to bus b1, by modifying the
-    # constraint `eq_net_injection`.
-    set_normalized_coefficient(
-        model[:eq_net_injection]["b1", t],
-        x[t],
-        1.0,
-    )
-end
-
-# Solve the model
-UnitCommitment.optimize!(model)
-
-# Show optimal values for the x variables
-@show value.(x)
-```
-
 References
 ----------
 * [KnOsWa20] **Bernard Knueven, James Ostrowski and Jean-Paul Watson.** "On Mixed-Integer Programming Formulations for the Unit Commitment Problem". INFORMS Journal on Computing (2020). [DOI: 10.1287/ijoc.2019.0944](https://doi.org/10.1287/ijoc.2019.0944)

src/instance/read.jl

@@ -98,10 +98,6 @@ 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 = [-1.0 for t in 1:T],
-    )
 
     # Read buses
     for (bus_name, dict) in json["Buses"]
@@ -268,7 +264,6 @@ function _from_json(json; repair = true)
     instance = UnitCommitmentInstance(
         T,
         power_balance_penalty,
-        shortfall_penalty,
         units,
         buses,
         lines,

src/instance/structs.jl

@@ -72,8 +72,6 @@ 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}

src/model/formulations/GenMorRam2017/startstop.jl

@@ -1,96 +0,0 @@
-# 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 Knueven et al. (2020).
-
-Creates constraints `eq_startstop_limit`,  `eq_startup_limit`, and `eq_shutdown_limit`
-using variables `Gar1962.StatusVars`, `prod_above` from `Gar1962.ProdVars`, and `reserve`.
-
-Constraints
----
-* `eq_startstop_limit`
-* `eq_startup_limit`
-* `eq_shutdown_limit`
-"""
-function _add_startup_shutdown_limit_eqs!(
-    model::JuMP.Model,
-    g::Unit,
-    formulation_prod_vars::Gar1962.ProdVars,
-    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
-
-    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)
-    end
-
-    for t in 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 Knueven 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.0
-                )
-            )
-        else
-            ## Startup limits
-            # Equation (23a) in Knueven 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.0
-                )
-            )
-
-            ## Shutdown limits
-            if t < T
-                # Equation (23b) in Knueven 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.0
-                    ) -
-                    max(0, g.shutdown_limit - g.startup_limit) *
-                    switch_on[gi, t]
-                )
-            end
-        end
-    end
-end

src/model/formulations/MorLatRam2013/scosts.jl

@@ -12,14 +12,14 @@ function _add_startup_cost_eqs!(
     S = length(g.startup_categories)
     startup = model[:startup]
     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(
-            model,
-            model[:switch_on][g.name, t] ==
-            sum(startup[g.name, t, s] for s in 1:S)
-        )
-
         for s in 1:S
+            # If unit is switching on, we must choose a startup category
+            eq_startup_choose[g.name, t, s] = @constraint(
+                model,
+                model[:switch_on][g.name, t] ==
+                sum(startup[g.name, t, s] for s in 1:S)
+            )
+
             # If unit has not switched off in the last `delay` time periods, startup category is forbidden.
             # The last startup category is always allowed.
             if s < S

src/model/formulations/MorLatRam2013/startstop.jl

@@ -1,89 +0,0 @@
-# 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 Knueven et al. (2020).
-
-Uses variable `prod_above` from `Gar1962.ProdVars`, the variables in `Gar1962.StatusVars`, and `reserve`
-to generate constraints below.
-
-Constraints
----
-* :eq_startstop_limit
-* :eq_startup_limit
-* :eq_shutdown_limit
-"""
-function _add_startup_shutdown_limit_eqs!(
-    model::JuMP.Model,
-    g::Unit,
-    formulation_prod_vars::Gar1962.ProdVars,
-    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
-
-    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 in 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 Knueven 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 Knueven 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 Knueven 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.0) <=
-                    (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
-    end
-end

src/model/formulations/base/bus.jl

@@ -4,11 +4,15 @@
 
 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
         net_injection[b.name, t] = AffExpr(-b.load[t])
 
+        # Reserves
+        reserve[b.name, t] = AffExpr()
+
         # Load curtailment
         curtail[b.name, t] =
             @variable(model, lower_bound = 0, upper_bound = b.load[t])

src/model/formulations/base/system.jl

@@ -11,12 +11,12 @@ end
 function _add_net_injection_eqs!(model::JuMP.Model)::Nothing
     T = model[:instance].time
     net_injection = _init(model, :net_injection)
-    eq_net_injection = _init(model, :eq_net_injection)
+    eq_net_injection_def = _init(model, :eq_net_injection_def)
     eq_power_balance = _init(model, :eq_power_balance)
     for t in 1:T, b in model[:instance].buses
         n = net_injection[b.name, t] = @variable(model)
-        eq_net_injection[b.name, t] =
-            @constraint(model, -n + model[:expr_net_injection][b.name, t] == 0)
+        eq_net_injection_def[t, b.name] =
+            @constraint(model, n == model[:expr_net_injection][b.name, t])
     end
     for t in 1:T
         eq_power_balance[t] = @constraint(
@@ -29,28 +29,13 @@ end
 
 function _add_reserve_eqs!(model::JuMP.Model)::Nothing
     eq_min_reserve = _init(model, :eq_min_reserve)
-    instance = model[:instance]
-    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]
+    for t in 1:model[:instance].time
         eq_min_reserve[t] = @constraint(
             model,
-            sum(model[:reserve][g.name, t] for g in instance.units) +
-            (shortfall_penalty >= 0 ? model[:reserve_shortfall][t] : 0.0) >=
-            instance.reserves.spinning[t]
+            sum(
+                model[:expr_reserve][b.name, t] for b in model[:instance].buses
+            ) >= model[:instance].reserves.spinning[t]
         )
-
-        # Account for shortfall contribution to objective
-        if shortfall_penalty >= 0
-            add_to_expression!(
-                model[:obj],
-                shortfall_penalty,
-                model[:reserve_shortfall][t],
-            )
-        end
     end
     return
 end

src/model/formulations/base/unit.jl

@@ -35,12 +35,7 @@ function _add_unit!(model::JuMP.Model, g::Unit, formulation::Formulation)
         formulation.status_vars,
     )
     _add_startup_cost_eqs!(model, g, formulation.startup_costs)
-    _add_startup_shutdown_limit_eqs!(
-        model,
-        g,
-        formulation.prod_vars,
-        formulation.status_vars,
-    )
+    _add_startup_shutdown_limit_eqs!(model, g)
     _add_status_eqs!(model, g, formulation.status_vars)
     return
 end
@@ -49,16 +44,12 @@ _is_initially_on(g::Unit)::Float64 = (g.initial_status > 0 ? 1.0 : 0.0)
 
 function _add_reserve_vars!(model::JuMP.Model, g::Unit)::Nothing
     reserve = _init(model, :reserve)
-    reserve_shortfall = _init(model, :reserve_shortfall)
     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
         end
-        reserve_shortfall[t] =
-            (model[:instance].shortfall_penalty[t] >= 0) ?
-            @variable(model, lower_bound = 0) : 0.0
     end
     return
 end
@@ -81,22 +72,7 @@ function _add_startup_shutdown_vars!(model::JuMP.Model, g::Unit)::Nothing
     return
 end
 
-"""
-    _add_startup_shutdown_limit_eqs!(model::JuMP.Model, g::Unit)::Nothing
-
-Creates startup/shutdown limit constraints below based on variables `Gar1962.StatusVars`, `prod_above` from `Gar1962.ProdVars`, and `reserve`.
-
-Constraints
----
-* :eq_startup_limit
-* :eq_shutdown_limit
-"""
-function _add_startup_shutdown_limit_eqs!(
-    model::JuMP.Model,
-    g::Unit,
-    formulation_prod_vars::Gar1962.ProdVars,
-    formulation_status_vars::Gar1962.StatusVars,
-)::Nothing
+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)
     is_on = model[:is_on]
@@ -234,5 +210,11 @@ 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

src/solution/solution.jl

@@ -51,12 +51,6 @@ function solution(model::JuMP.Model)::OrderedDict
     sol["Switch on"] = timeseries(model[:switch_on], instance.units)
     sol["Switch off"] = timeseries(model[:switch_off], instance.units)
     sol["Reserve (MW)"] = timeseries(model[:reserve], instance.units)
-    sol["Reserve shortfall (MW)"] = OrderedDict(
-        t =>
-            (instance.shortfall_penalty[t] >= 0) ?
-            round(value(model[:reserve_shortfall][t]), digits = 5) : 0.0 for
-        t in 1:instance.time
-    )
     sol["Net injection (MW)"] =
         timeseries(model[:net_injection], instance.buses)
     sol["Load curtail (MW)"] = timeseries(model[:curtail], instance.buses)

src/validation/validate.jl

@@ -324,16 +324,11 @@ function _validate_reserve_and_demand(instance, solution, tol = 0.01)
         # Verify spinning reserves
         reserve =
             sum(solution["Reserve (MW)"][g.name][t] for g in instance.units)
-        reserve_shortfall =
-            (instance.shortfall_penalty[t] >= 0) ?
-            solution["Reserve shortfall (MW)"][t] : 0
-
-        if reserve + reserve_shortfall < instance.reserves.spinning[t] - tol
+        if reserve < instance.reserves.spinning[t] - tol
             @error @sprintf(
-                "Insufficient spinning reserves at time %d (%.2f + %.2f should be %.2f)",
+                "Insufficient spinning reserves at time %d (%.2f should be %.2f)",
                 t,
                 reserve,
-                reserve_shortfall,
                 instance.reserves.spinning[t],
             )
             err_count += 1

test/model/formulations_test.jl

@@ -20,14 +20,8 @@ if ENABLE_LARGE_TESTS
 end
 
 function _small_test(formulation::Formulation)::Nothing
-    instances = ["matpower/case118/2017-02-01", "test/case14"]
-    for instance in instances
-        # Should not crash
-        UnitCommitment.build_model(
-            instance = UnitCommitment.read_benchmark(instance),
-            formulation = formulation,
-        )
-    end
+    instance = UnitCommitment.read_benchmark("matpower/case118/2017-02-01")
+    UnitCommitment.build_model(instance = instance, formulation = formulation)  # should not crash
     return
 end