Missed function definition

.github/workflows/test.yml

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

CHANGELOG.md

@@ -11,9 +11,10 @@ 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] - Unreleased
+## [0.2.2] - 2021-07-21
 ### 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,13 +28,14 @@ Each section is described in detail below. For a complete example, see [case14](
 
 ### Parameters
 
-This section describes system-wide parameters, such as power balance penalties,  optimization parameters, such as the length of the planning horizon and the time.
+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.
 
 | Key                            | Description                                       | Default  | Time series?
 | :----------------------------- | :------------------------------------------------ | :------: | :------------:
 | `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
@@ -42,7 +43,8 @@ This section describes system-wide parameters, such as power balance penalties,
 {
     "Parameters": {
         "Time horizon (h)": 4,
-        "Power balance penalty ($/MW)": 1000.0
+        "Power balance penalty ($/MW)": 1000.0,
+        "Reserve shortfall penalty ($/MW)": -1.0
     }
 }
 ```

src/instance/read.jl

@@ -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 = [-1.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,

src/instance/structs.jl

@@ -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}

src/model/formulations/Gar1962/status.jl

@@ -2,6 +2,12 @@
 # Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
 # Released under the modified BSD license. See COPYING.md for more details.
 
+"""
+    _add_status_vars!
+
+Adds symbols identified by `Gar1962.StatusVars` to `model`.
+Fix variables if a certain generator _must_ run or based on initial conditions.
+"""
 function _add_status_vars!(
     model::JuMP.Model,
     g::Unit,
@@ -10,15 +16,93 @@ function _add_status_vars!(
     is_on = _init(model, :is_on)
     switch_on = _init(model, :switch_on)
     switch_off = _init(model, :switch_off)
+    FIX_VARS = !formulation_status_vars.fix_vars_via_constraint
+    is_initially_on = _is_initially_on(g) > 0
     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
         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)
+
+        # Use initial conditions and whether a unit must run to fix variables
+        if FIX_VARS
+            # Fix variables using fix function
+            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)
+            elseif t == 1
+                if is_initially_on
+                    # 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
+        else
+            # Add explicit constraint if !FIX_VARS
+            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
+            elseif t == 1
+                if is_initially_on
+                    switch_on[g.name, t] = 0.0
+                else
+                    switch_off[g.name, t] = 0.0
+                end
+            end
+        end
+
+        # Use initial conditions and whether a unit must run to fix variables
+        if FIX_VARS
+            # Fix variables using fix function
+            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)
+            elseif t == 1
+                if is_initially_on
+                    # 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
+        else
+            # Add explicit constraint if !FIX_VARS
+            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
+            elseif t == 1
+                if is_initially_on
+                    switch_on[g.name, t] = 0.0
+                else
+                    switch_off[g.name, t] = 0.0
+                end
+            end
         end
     end
     return

src/model/formulations/Gar1962/structs.jl

@@ -17,8 +17,53 @@ import ..PiecewiseLinearCostsFormulation
 import ..ProductionVarsFormulation
 import ..StatusVarsFormulation
 
+"""
+Variables
+---
+* `prod_above`:
+        [gen, t];
+        *production above minimum required level*;
+        lb: 0, ub: Inf.
+        KnuOstWat2020: `p'_g(t)`
+* `segprod`:
+        [gen, segment, t];
+        *how much generator produces on cost segment in time t*;
+        lb: 0, ub: Inf.
+        KnuOstWat2020: `p_g^l(t)`
+"""
 struct ProdVars <: ProductionVarsFormulation end
+
 struct PwlCosts <: PiecewiseLinearCostsFormulation end
-struct StatusVars <: StatusVarsFormulation end
+
+"""
+Variables
+---
+* `is_on`:
+        [gen, t];
+        *is generator on at time t?*
+        lb: 0, ub: 1, binary.
+        KnuOstWat2020: `u_g(t)`
+* `switch_on`:
+        [gen, t];
+        *indicator that generator will be turned on at t*;
+        lb: 0, ub: 1, binary.
+        KnuOstWat2020: `v_g(t)`
+* `switch_off`: binary;
+        [gen, t];
+        *indicator that generator will be turned off at t*;
+        lb: 0, ub: 1, binary.
+        KnuOstWat2020: `w_g(t)`
+
+Arguments
+---
+* `fix_vars_via_constraint`:
+        indicator for whether to set vars to a constant using `fix` or by adding an explicit constraint
+        (particulary useful for debugging purposes).
+"""
+struct StatusVars <: StatusVarsFormulation
+    fix_vars_via_constraint::Bool
+
+    StatusVars() = new(false)
+end
 
 end

src/model/formulations/base/bus.jl

@@ -4,15 +4,11 @@
 
 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

@@ -29,13 +29,28 @@ end
 
 function _add_reserve_eqs!(model::JuMP.Model)::Nothing
     eq_min_reserve = _init(model, :eq_min_reserve)
-    for t in 1:model[:instance].time
+    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]
         eq_min_reserve[t] = @constraint(
             model,
-            sum(
+            sum(model[:reserve][g.name, t] for g in instance.units) +
-                model[:expr_reserve][b.name, t] for b in model[:instance].buses
+            (shortfall_penalty >= 0 ? model[:reserve_shortfall][t] : 0.0) >=
-            ) >= model[:instance].reserves.spinning[t]
+            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

@@ -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,7 +44,12 @@ 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)
+    _add_startup_shutdown_limit_eqs!(
+        model,
+        g,
+        formulation.status_vars,
+        formulation.prod_vars,
+    )
     _add_status_eqs!(model, g, formulation.status_vars)
     return
 end
@@ -44,12 +58,16 @@ _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
@@ -72,7 +90,22 @@ function _add_startup_shutdown_vars!(model::JuMP.Model, g::Unit)::Nothing
     return
 end
 
-function _add_startup_shutdown_limit_eqs!(model::JuMP.Model, g::Unit)::Nothing
+"""
+    _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_status_vars::Gar1962.StatusVars,
+    formulation_prod_vars::Gar1962.ProdVars,
+)::Nothing
     eq_shutdown_limit = _init(model, :eq_shutdown_limit)
     eq_startup_limit = _init(model, :eq_startup_limit)
     is_on = model[:is_on]
@@ -91,8 +124,15 @@ function _add_startup_shutdown_limit_eqs!(model::JuMP.Model, g::Unit)::Nothing
         )
         # Shutdown limit
         if g.initial_power > g.shutdown_limit
+            # 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?)
+            if formulation_status_vars.fix_vars_via_constraint
                 eq_shutdown_limit[g.name, 0] =
-                @constraint(model, switch_off[g.name, 1] <= 0)
+                    @constraint(model, model[:switch_off][g.name, 1] <= 0.0)
+            else
+                fix(model[:switch_off][g.name, 1], 0.0; force = true)
+            end
         end
         if t < T
             eq_shutdown_limit[g.name, t] = @constraint(
@@ -210,11 +250,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

src/solution/solution.jl

@@ -51,6 +51,12 @@ 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,11 +324,16 @@ 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)
-        if reserve < instance.reserves.spinning[t] - tol
+        reserve_shortfall =
+            (instance.shortfall_penalty[t] >= 0) ?
+            solution["Reserve shortfall (MW)"][t] : 0
+
+        if reserve + reserve_shortfall < instance.reserves.spinning[t] - tol
             @error @sprintf(
-                "Insufficient spinning reserves at time %d (%.2f should be %.2f)",
+                "Insufficient spinning reserves at time %d (%.2f + %.2f should be %.2f)",
                 t,
                 reserve,
+                reserve_shortfall,
                 instance.reserves.spinning[t],
             )
             err_count += 1

test/model/formulations_test.jl

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