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

Initial public release

Browse Source
This commit is contained in:
Alinson S. Xavier
2020-11-03 15:49:12 -06:00
commit b2480ef356
162 changed files with 4397 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

15
src/UnitCommitment.jl Normal file
View File

@@ -0,0 +1,15 @@
# 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.
module UnitCommitment
include("log.jl")
include("dotdict.jl")
include("instance.jl")
include("screening.jl")
include("model.jl")
include("sensitivity.jl")
include("validate.jl")
include("convert.jl")
include("initcond.jl")
end

56
src/convert.jl Normal file
View File

@@ -0,0 +1,56 @@
# 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.
using DataStructures, JSON, GZip
function read_json(path::String)::OrderedDict
if endswith(path, ".gz")
file = GZip.gzopen(path)
else
file = open(path)
end
return JSON.parse(file, dicttype=()->DefaultOrderedDict(nothing))
end
function read_egret_solution(path::String)::OrderedDict
egret = read_json(path)
T = length(egret["system"]["time_keys"])
solution = OrderedDict()
is_on = solution["Is on"] = OrderedDict()
production = solution["Production (MW)"] = OrderedDict()
reserve = solution["Reserve (MW)"] = OrderedDict()
production_cost = solution["Production cost (\$)"] = OrderedDict()
startup_cost = solution["Startup cost (\$)"] = OrderedDict()
for (gen_name, gen_dict) in egret["elements"]["generator"]
if endswith(gen_name, "_T") || endswith(gen_name, "_R")
gen_name = gen_name[1:end-2]
end
if "commitment" in keys(gen_dict)
is_on[gen_name] = gen_dict["commitment"]["values"]
else
is_on[gen_name] = ones(T)
end
production[gen_name] = gen_dict["pg"]["values"]
if "rg" in keys(gen_dict)
reserve[gen_name] = gen_dict["rg"]["values"]
else
reserve[gen_name] = zeros(T)
end
startup_cost[gen_name] = zeros(T)
production_cost[gen_name] = zeros(T)
if "commitment_cost" in keys(gen_dict)
for t in 1:T
x = gen_dict["commitment"]["values"][t]
commitment_cost = gen_dict["commitment_cost"]["values"][t]
prod_above_cost = gen_dict["production_cost"]["values"][t]
prod_base_cost = gen_dict["p_cost"]["values"][1][2] * x
startup_cost[gen_name][t] = commitment_cost - prod_base_cost
production_cost[gen_name][t] = prod_above_cost + prod_base_cost
end
end
end
return solution
end

28
src/docs/css/custom.css Normal file
View File

@@ -0,0 +1,28 @@
.navbar-default {
border-bottom: 0px;
background-color: #fff;
box-shadow: 0px 0px 15px rgba(0, 0, 0, 0.2);
}
a, .navbar-default a {
color: #06a !important;
font-weight: normal;
}
.disabled > a {
color: #999 !important;
}
.navbar-default a:hover,
.navbar-default .active,
.active > a {
background-color: #f0f0f0 !important;
}
.icon-bar {
background-color: #666 !important;
}
.navbar-collapse {
border-color: #fff !important;
}

257
src/docs/format.md Normal file
View File

@@ -0,0 +1,257 @@
Data Format
===========
Instances are specified by JSON files containing the following main sections:
* [Parameters](#parameters)
* [Buses](#buses)
* [Generators](#generators)
* [Price-sensitive loads](#price-sensitive-loads)
* [Transmission lines](#transmission-lines)
* [Reserves](#reserves)
* [Contingencies](#contingencies)
Each section is described in detail below. For a complete example, see [case14.json](https://github.com/ANL-CEEESA/UnitCommitment.jl/blob/dev/instances/matpower-24h/case14.json).
## Parameters
This section describes system-wide parameters, such as power balance penalties, and optimization parameters, such as the length of the planning horizon.
| Key | Description | Default | Time series?
| :----------------------------- | :------------------------------------------------ | :------: | :------------:
| `Time (h)` | Length of the planning horizon (in hours) | Required | N
| `Power balance penalty ($/MW)` | Penalty for system-wide shortage or surplus in production (in $/MW). This is charged per time period. For example, if there is a shortage of 1 MW for three time periods, three times this amount will be charged. | `1000.0` | Y
### Example
```json
{
"Parameters": {
"Time (h)": 4,
"Power balance penalty ($/MW)": 1000.0
}
}
```
## Buses
This section describes the characteristics of each bus in the system.
| Key | Description | Default | Time series?
| :----------------- | :------------------------------------------------------------ | ------- | :-------------:
| `Load (MW)` | Fixed load connected to the bus (in MW). | Required | Y
### Example
```json
{
"Buses": {
"b1": {
"Load (MW)": 0.0
},
"b2": {
"Load (MW)": [
26.01527,
24.46212,
23.29725,
22.90897
]
}
}
}
```
## Generators
This section describes all generators in the system, including thermal units, renewable units and virtual units.
| Key | Description | Default | Time series?
| :------------------------ | :------------------------------------------------| ------- | :-----------:
| `Bus` | Identifier of the bus where this generator is located (string) | Required | N
| `Production cost curve (MW)` and `Production cost curve ($)` | Parameters describing the piecewise-linear production costs. See below for more details. | Required | Y
| `Startup costs ($)` and `Startup delays (h)` | Parameters describing how much it costs to start the generator after it has been shut down for a certain amount of time. If `Startup costs ($)` and `Startup delays (h)` are set to `[300.0, 400.0]` and `[1, 4]`, for example, and the generator is shut down at time `t`, then it costs 300 to start up the generator at times `t+1`, `t+2` or `t+3`, and 400 to start the generator at time `t+4` or any time after that. The number of startup cost points is unlimited, and may be different for each generator. Startup delays must be strictly increasing. | `[0.0]` and `[1]` | N
| `Minimum uptime (h)` | Minimum amount of time the generator must stay operational after starting up (in hours). For example, if the generator starts up at time 1 and `Minimum uptime (h)` is set to 4, then the generator can only shut down at time 5. | `1` | N
| `Minimum downtime (h)` | Minimum amount of time the generator must stay offline after shutting down (in hours). For example, if the generator shuts down at time 1 and `Minimum downtime (h)` is set to 4, then the generator can only start producing power again at time 5. | `1` | N
| `Ramp up limit (MW)` | Maximum increase in production from one time period to the next (in MW). For example, if the generator is producing 100 MW at time 1 and if this parameter is set to 40 MW, then the generator will produce at most 140 MW at time 2. | `+inf` | N
| `Ramp down limit (MW)` | Maximum decrease in production from one time period to the next (in MW). For example, if the generator is producing 100 MW at time 1 and this parameter is set to 40 MW, then the generator will produce at least 60 MW at time 2. | `+inf` | N
| `Startup limit (MW)` | Maximum amount of power a generator can produce immediately after starting up (in MW). | `+inf` | N
| `Shutdown limit (MW)` | Maximum amount of power a generator can produce immediately before shutting down (in MW). Specifically, the generator can only shut down at time `t+1` if its production at time `t` is below this limit. | `+inf` | N
| `Initial status (h)` | If set to a positive number, indicates the amount of time the generator has been on at the beginning of the simulation, and if set to a negative number, the amount of time the generator has been off. For example, if `Initial status (h)` is `-2`, this means that the generator was off at simulation time `-2` and `-1`. The simulation starts at time `0`. | Required | N
| `Initial power (MW)` | Amount of power the generator at time period `-1`, immediately before the planning horizon starts. | Required | N
| `Must run?` | If `true`, the generator should be committed, even that is not economical (Boolean). | `false` | Y
| `Provides spinning reserves?` | If `true`, this generator may provide spinning reserves (Boolean). | `true` | Y
### Production costs and limits
Production costs are represented as piecewise-linear curves. Figure 1 shows an example cost curve with three segments, where it costs 1400, 1600, 2200 and 2400 dollars to generate, respectively, 100, 110, 130 and 135 MW of power. To model this generator, `Production cost curve (MW)` should be set to `[100, 110, 130, 135]`, and `Production cost curve ($)` should be set to `[1400, 1600, 2200, 2400]`.
Note that this curve also specifies the production limits. Specifically, the first point identifies the minimum power output when the unit is operational, while the last point identifies the maximum power output.
<center>
<img src="../images/cost_curve.png" style="max-width: 500px"/>
<div><b>Figure 1.</b> Piecewise-linear production cost curve.</div>
<br/>
</center>
**Additional remarks:**
* For time-dependent production limits or time-dependent production costs, the usage of nested arrays is allowed. For example, if `Production cost curve (MW)` is set to `[5.0, [10.0, 12.0, 15.0, 20.0]]`, then the unit may generate at most 10, 12, 15 and 20 MW of power during time periods 1, 2, 3 and 4, respectively. The minimum output for all time periods is fixed to at 5 MW.
* There is no limit to the number of piecewise-linear segments, and different generators may have a different number of segments.
* If `Production cost curve (MW)` and `Production cost curve ($)` both contain a single element, then the generator must produce exactly that amount of power when operational. To specify that the generator may produce any amount of power up to a certain limit `P`, the parameter `Production cost curve (MW)` should be set to `[0, P]`.
* Production cost curves must be convex.
### Example
```json
{
"Generators": {
"gen1": {
"Bus": "b1",
"Production cost curve (MW)": [100.0, 110.0, 130.0, 135.0],
"Production cost curve ($)": [1400.0, 1600.0, 2200.0, 2400.0],
"Startup costs ($)": [300.0, 400.0],
"Startup delays (h)": [1, 4],
"Ramp up limit (MW)": 232.68,
"Ramp down limit (MW)": 232.68,
"Startup limit (MW)": 232.68,
"Shutdown limit (MW)": 232.68,
"Minimum downtime (h)": 4,
"Minimum uptime (h)": 4,
"Initial status (h)": 12,
"Must run?": false,
"Provides spinning reserves?": true,
},
"gen2": {
"Bus": "b5",
"Production cost curve (MW)": [0.0, [10.0, 8.0, 0.0, 3.0]],
"Production cost curve ($)": [0.0, 0.0],
"Provides spinning reserves?": true,
}
}
}
```
## Price-sensitive loads
This section describes components in the system which may increase or reduce their energy consumption according to the energy prices. Fixed loads (as described in the `buses` section) are always served, regardless of the price, unless there is significant congestion in the system or insufficient production capacity. Price-sensitive loads, on the other hand, are only served if it is economical to do so.
| Key | Description | Default | Time series?
| :---------------- | :------------------------------------------------ | :------: | :------------:
| `Bus` | Bus where the load is located. Multiple price-sensitive loads may be placed at the same bus. | Required | N
| `Revenue ($/MW)` | Revenue obtained for serving each MW of power to this load. | Required | Y
| `Demand (MW)` | Maximum amount of power required by this load. Any amount lower than this may be served. | Required | Y
### Example
```json
{
"Price-sensitive loads": {
"p1": {
"Bus": "b3",
"Revenue ($/MW)": 23.0,
"Demand (MW)": 50.0
}
}
}
```
## Transmission Lines
This section describes the characteristics of transmission system, such as its topology and the susceptance of each transmission line.
| Key | Description | Default | Time series?
| :--------------------- | :----------------------------------------------- | ------- | :------------:
| `Source bus` | Identifier of the bus where the transmission line originates. | Required | N
| `Target bus` | Identifier of the bus where the transmission line reaches. | Required | N
| `Reactance (ohms)` | Reactance of the transmission line (in ohms). | Required | N
| `Susceptance (S)` | Susceptance of the transmission line (in siemens). | Required | N
| `Normal flow limit (MW)` | Maximum amount of power (in MW) allowed to flow through the line when the system is in its regular, fully-operational state. May be `null` is there is no limit. | `+inf` | Y
| `Emergency flow limit (MW)` | Maximum amount of power (in MW) allowed to flow through the line when the system is in degraded state (for example, after the failure of another transmission line). | `+inf` | Y
| `Flow limit penalty ($/MW)` | Penalty for violating the flow limits of the transmission line (in $/MW). This is charged per time period. For example, if there is a thermal violation of 1 MW for three time periods, three times this amount will be charged. | `5000.0` | Y
### Example
```json
{
"Transmission lines": {
"l1": {
"Source bus": "b1",
"Target bus": "b2",
"Reactance (ohms)": 0.05917,
"Susceptance (S)": 29.49686,
"Normal flow limit (MW)": 15000.0,
"Emergency flow limit (MW)": 20000.0,
"Flow limit penalty ($/MW)": 5000.0
}
}
}
```
## Reserves
This section describes the hourly amount of operating reserves required.
| Key | Description | Default | Time series?
| :-------------------- | :------------------------------------------------- | --------- | :----:
| `Spinning (MW)` | Minimum amount of system-wide spinning reserves (in MW). Only generators which are online may provide this reserve. | `0.0` | Y
### Example
```json
{
"Reserves": {
"Spinning (MW)": [
57.30552,
53.88429,
51.31838,
50.46307
]
}
}
```
## Contingencies
This section describes credible contingency scenarios in the optimization, such as the loss of a transmission line or generator.
| Key | Description | Default
| :-------------------- | :----------------------------------------------- | ----------
| `Affected generators` | List of generators affected by this contingency. May be omitted if no generators are affected. | `[]`
| `Affected lines` | List of transmission lines affected by this contingency. May be omitted if no lines are affected. | `[]`
### Example
```json
{
"Contingencies": {
"c1": {
"Affected lines": ["l1", "l2", "l3"],
"Affected generators": ["g1"]
},
"c2": {
"Affected lines": ["l4"]
},
}
}
```
## Additional remarks
### Time series parameters
Many numerical properties in the JSON file can be specified either as a single floating point number if they are time-independent, or as an array containing exactly `T` elements, where `T` is the length of the planning horizon, if they are time-dependent. For example, both formats below are valid when `T=3`:
```json
{
"Load (MW)": 800.0,
"Load (MW)": [800.0, 850.0, 730.0]
}
```
### Current limitations
* All reserves are system-wide (no zonal reserves)
* Network topology remains the same for all time periods
* Only N-1 transmission contingencies are supported. Generator contingencies are not supported.

60
src/docs/illustrations.ipynb Normal file
View File

File diff suppressed because one or more lines are too long

BIN
src/docs/images/cost_curve.png Normal file
View File

Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 35 KiB

23
src/docs/index.md Normal file
View File

@@ -0,0 +1,23 @@
# UnitCommitment.jl
**UnitCommitment.jl** is an optimization package for the Security-Constrained Unit Commitment Problem (SCUC), a fundamental optimization problem in power systems which is used, for example, to clear the day-ahead electricity markets. The problem asks for the most cost-effective power generation schedule under a number of physical, operational and economic constraints.
### Package Components
* **Data Format:** The package proposes an extensible and fully-documented JSON-based data specification format for SCUC, developed in collaboration with Independent System Operators (ISOs), which describes the most important aspects of the problem.
* **Benchmark Instances:** The package provides a diverse collection of large-scale benchmark instances collected from the literature and extended to make them more challenging and realistic, based on publicly available data from ISOs.
* **Model Implementation**: The package provides a Julia/JuMP implementation of state-of-the-art formulations and solution methods for SCUC. Our goal is to keep this implementation up-to-date, as new methods are proposed in the literature.
* **Benchmark Tools:** The package provides automated benchmark scripts to accurately evaluate the performance impact of proposed code changes.
### Contents
* [Installation Guide](install.md)
* [Data Format Specification](format.md)
### Authors
* **Alinson Santos Xavier,** Argonne National Laboratory
* **Feng Qiu,** Argonne National Laboratory
### Collaborators
* **Yonghong Chen,** Midcontinent Independent System Operator
* **Feng Pan,** Pacific Northwest National Laboratory

21
src/docs/install.md Normal file
View File

@@ -0,0 +1,21 @@
# Installation Guide
This package was tested and developed with [Julia 1.4](https://julialang.org/). To install Julia, please
follow the [installation guide on their website](https://julialang.org/downloads/platform.html).
To install `UnitCommitment.jl`, run the Julia interpreter, type `]` to open the
package manager, then type:
```text
pkg> add https://github.com/ANL-CEEESA/UnitCommitment.jl.git
```
To test that the package has been correctly installed, run:
```text
pkg> test UnitCommitment
```
If all tests pass, the package should now be ready to be used by any Julia script on the machine. To try it out in the julia interpreter hit `backspace` to return to the regular interpreter, and type the following command:
```julia
using UnitCommitment
```

19
src/docs/isf.md Normal file
View File

@@ -0,0 +1,19 @@
Linear Sensitivity Factors
==========================
UnitCommitment.jl includes a number of functions to compute typical linear sensitivity
factors, such as [Injection Shift Factors](@ref) and [Line Outage Distribution Factors](@ref). These sensitivity factors can be used to quickly compute DC power flows in both base and N-1 contigency scenarios.
Injection Shift Factors
-----------------------
Given a network with `B` buses and `L` transmission lines, the Injection Shift Factors (ISF) matrix is an `L`-by-`B` matrix which indicates much power flows through a certain transmission line when 1 MW of power is injected at bus `b` and withdrawn from the slack bus. For example, `isf[:l7, :b5]` indicates the amount of power (in MW) that flows through line `l7` when 1 MW of power is injected at bus `b5` and withdrawn from the slack bus.
This matrix is computed based on the DC linearization of power flow equations and does not include losses.
To compute the ISF matrix, the function `injection_shift_factors` can be used. It is necessary to specify the set of lines, buses and the slack bus:
```julia
using UnitCommitment
instance = UnitCommitment.load("ieee_rts/case14")
isf = UnitCommitment.injection_shift_factors(lines = instance.lines,
buses = instance.buses,
slack = :b14)
@show isf[:l7, :b5]
```

8
src/docs/js/mathjax.js Normal file
View File

@@ -0,0 +1,8 @@
MathJax.Hub.Config({
"tex2jax": { inlineMath: [ [ '$', '$' ] ] }
});
MathJax.Hub.Config({
config: ["MMLorHTML.js"],
jax: ["input/TeX", "output/HTML-CSS", "output/NativeMML"],
extensions: ["MathMenu.js", "MathZoom.js"]
});

90
src/docs/model.md Normal file
View File

@@ -0,0 +1,90 @@
Benchmark Model
===============
UnitCommitment.jl includes a reference Mixed-Integer Linear Programming
(MILP), built with [JuMP](https://github.com/JuliaOpt/JuMP.jl), which can
either be used as-is to solve instances of the problem, or be extended to
build more complex formulations.
Building and Solving the Model
-------------------------------
Given an instance and a JuMP optimizer, the function `build_model` can be used to
build the reference MILP model. For example:
```julia
using UnitCommitment, JuMP, Cbc
instance = UnitCommitment.load("ieee_rts/case118")
model = build_model(instance, with_optimizer(Cbc.Optimizer))
```
The model enforces all unit constraints described in [Unit Commitment
Instances](@ref), including ramping, minimum-up and minimum-down times. Some
system-wide constraints, such as spinning reserves, are also enforced. The
model, however, does not enforce transmission or N-1 security constraints,
since these are typically generated on-the-fly.
A reference to the JuMP model is stored at `model.mip`. After constructed, the model can
be optimized as follows:
```julia
optimize!(model.mip)
```
Decision Variables
------------------
References to all decision variables are stored at `model.vars`.
A complete list of available decision variables is as follows:
| Variable | Description
| :---------------------------- | :---------------------------
| `model.vars.production[gi,t]` | Amount of power (in MW) produced by unit with index `gi` at time `t`.
| `model.vars.reserve[gi,t]` | Amount of spinning reserves (in MW) provided by unit with index `gi` at time `t`.
| `model.vars.is_on[gi,t]` | Binary variable indicating if unit with index `gi` is operational at time `t`.
| `model.vars.switch_on[gi,t]` | Binary variable indicating if unit with index `gi` was switched on at time `t`. That is, the unit was not operational at time `t-1`, but it is operational at time `t`.
| `model.vars.switch_off[gi,t]` | Binary variable indicating if unit with index `gi` was switched off at time `t`. That is, the unit was operational at time `t-1`, but it is no longer operational at time `t`.
| `model.vars.unit_cost[gi,t]` | The total cost to operate unit with index `gi` at time `t`. Includes start-up costs, no-load costs and any other production costs.
| `model.vars.cost[t]` | Total cost at time `t`.
| `model.vars.net_injection[bi,t]` | Total net injection (in MW) at bus with index `bi` and time `t`. Net injection is defined as the total power being produced by units located at the bus minus the bus load.
Accessing the Solution
----------------------
To access the value of a particular decision variable after the
optimization is completed, the function `JuMP.value(var)` can be used. The
following example prints the amount of power (in MW) produced by each unit at time 5:
```julia
for g in instance.units
@show value(model.vars.production[g.index, 5])
end
```
Modifying the Model
-------------------
Prior to being solved, the reference model can be modified by using the variable references
above and conventional JuMP macros. For example, the
following code can be used to ensure that at most 10 units are operational at time 4:
```julia
using UnitCommitment, JuMP, Cbc
instance = UnitCommitment.load("ieee_rts/case118")
model = build_model(instance, with_optimizer(Cbc.Optimizer))
@contraint(model.mip,
sum(model.vars.is_on[g.index, 4]
for g in instance.units) <= 10)
optimize!(model.mip)
```
It is not currently possible to modify the constraints included in the
reference model.
Reference
---------
```@docs
UnitCommitment.build_model
```

81
src/docs/solvers.md Normal file
View File

@@ -0,0 +1,81 @@
Benchmark Solver
================
Solving an instance of the Unit Commitment problem typically involves more
than simply building a Mixed-Integer Linear Programming and handing it over to the
solver. Since the number of transmission and N-1 security constraints can
easily exceed hundreds of millions for large instances of the problem, it is often
necessary to iterate between MILP optimization and contingency screening, so
that only necessary transmission constraints are added to the MILP.
`UnitCommitment.jl` includes a fast implementation of the contingency
screening method described in
[[1]](https://doi.org/10.1109/TPWRS.2019.2892620), which is able to
efficiently handle even ISO-scale instances of the problem. The method makes
use of Injection Shift Factors (ISFs) and Line Outage Distribution Factors
(LODFs) to model DC power flows and N-1 contingencies. If Julia is configured
to use multiple threads (through the environment variable `JULIA_NUM_THREADS`)
then multiple contingency scenarios are evaluated in parallel.
Usage
-----
To solve one of the benchmark instances using the included benchmark solver, use the method `UnitCommitment.solve`
as shown in the example below.
julia> UnitCommitment.solve("ieee_rts/case118")
[ Info: Loading instance: ieee_rts/case118
[ Info: 54 units
[ Info: 118 buses
[ Info: 186 lines
[ Info: Scaling problem (0.6 demands, 1.0 limits)...
[ Info: Using Cbc as MILP solver (0.001 gap, 4 threads)
[ Info: Computing sensitivity factors (0.001 ISF cutoff, 0.0001 LODF cutoff)...
[ Info: Building MILP model (24 hours, 0.01 reserve)...
[ Info: Optimizing...
[ Info: Optimal value: 4.033106e+06
[ Info: Solved in 8.73 seconds
With default settings, the solver does not consider any transmission or
security constraints, and the peak load is automatically set to 60% of the
installed capacity of the system. These, and many other settings, can be
configured using keyword arguments. See the reference section below for more
details. Sample usage:
julia> UnitCommitment.solve("ieee_rts/case118", demand_scale=0.7, security=true)
[ Info: Loading instance: ieee_rts/case118
[ Info: 54 units
[ Info: 118 buses
[ Info: 186 lines
[ Info: Scaling problem (0.7 demands, 1.0 limits)...
[ Info: Using Cbc as MILP solver (0.001 gap, 4 threads)
[ Info: Computing sensitivity factors (0.001 ISF cutoff, 0.0001 LODF cutoff)...
[ Info: Building MILP model (24 hours, 0.01 reserve)...
[ Info: Optimizing...
[ Info: Verifying flow constraints...
[ Info: Optimal value: 4.888740e+06
[ Info: Solved in 4.50 seconds
When transmission or N-1 security constraints are activated, the solver uses an
iterative method to lazily enforce them. See [this
paper](https://doi.org/10.1109/TPWRS.2019.2892620) for a detailed description
of the method. Injection Shift Factors (ISF) and Line Outage Distribution
Factors (LODF) are used for the computation of DC power flows.
!!! note
Many of the benchmark instances were not originally designed for N-1
security-constrained studies, and may become infeasible if these constraints
are enforced. To avoid infeasibilities, the transmission limits can be
increased through the keyword argument `limit_scale`.
By default, the MILP is solved using [Cbc, the COIN-OR Branch and Cut
solver](https://github.com/coin-or/Cbc). If `UnitCommitment` is loaded after
either [CPLEX](https://github.com/JuliaOpt/CPLEX.jl) or
[SCIP](https://github.com/SCIP-Interfaces/SCIP.jl), then these solvers will be
used instead. A detailed solver log can be displayed by setting `verbose=true`.
## Reference
```@docs
UnitCommitment.solve
```

68
src/dotdict.jl Normal file
View File

@@ -0,0 +1,68 @@
# 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.
struct DotDict
inner::Dict
end
DotDict() = DotDict(Dict())
function Base.setproperty!(d::DotDict, key::Symbol, value)
setindex!(getfield(d, :inner), value, key)
end
function Base.getproperty(d::DotDict, key::Symbol)
(key == :inner ? getfield(d, :inner) : d.inner[key])
end
function Base.getindex(d::DotDict, key::Int64)
d.inner[Symbol(key)]
end
function Base.getindex(d::DotDict, key::Symbol)
d.inner[key]
end
function Base.keys(d::DotDict)
keys(d.inner)
end
function Base.values(d::DotDict)
values(d.inner)
end
function Base.iterate(d::DotDict)
iterate(values(d.inner))
end
function Base.iterate(d::DotDict, v::Int64)
iterate(values(d.inner), v)
end
function Base.length(d::DotDict)
length(values(d.inner))
end
function Base.show(io::IO, d::DotDict)
print(io, "DotDict with $(length(keys(d.inner))) entries:\n")
count = 0
for k in keys(d.inner)
count += 1
if count > 10
print(io, " ...\n")
break
end
print(io, " :$(k) => $(d.inner[k])\n")
end
end
function recursive_to_dot_dict(el)
if typeof(el) == Dict{String, Any}
return DotDict(Dict(Symbol(k) => recursive_to_dot_dict(el[k]) for k in keys(el)))
else
return el
end
end
export recursive_to_dot_dict

76
src/initcond.jl Normal file
View File

@@ -0,0 +1,76 @@
# 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.
using JuMP
"""
generate_initial_conditions!(instance, optimizer)
Generates feasible initial conditions for the given instance, by constructing
and solving a single-period mixed-integer optimization problem, using the given
optimizer. The instance is modified in-place.
"""
function generate_initial_conditions!(instance::UnitCommitmentInstance,
optimizer)
G = instance.units
B = instance.buses
t = 1
mip = JuMP.Model(optimizer)
# Decision variables
@variable(mip, x[G], Bin)
@variable(mip, p[G] >= 0)
# Constraint: Minimum power
@constraint(mip,
min_power[g in G],
p[g] >= g.min_power[t] * x[g])
# Constraint: Maximum power
@constraint(mip,
max_power[g in G],
p[g] <= g.max_power[t] * x[g])
# Constraint: Production equals demand
@constraint(mip,
power_balance,
sum(b.load[t] for b in B) == sum(p[g] for g in G))
# Constraint: Must run
for g in G
if g.must_run[t]
@constraint(mip, x[g] == 1)
end
end
# Objective function
function cost_slope(g)
mw = g.min_power[t]
c = g.min_power_cost[t]
for k in g.cost_segments
mw += k.mw[t]
c += k.mw[t] * k.cost[t]
end
if mw < 1e-3
return 0.0
else
return c / mw
end
end
@objective(mip,
Min,
sum(p[g] * cost_slope(g) for g in G))
JuMP.optimize!(mip)
for g in G
if JuMP.value(x[g]) > 0
g.initial_power = JuMP.value(p[g])
g.initial_status = 24
else
g.initial_power = 0
g.initial_status = -24
end
end
end

349
src/instance.jl Normal file
View File

@@ -0,0 +1,349 @@
# 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.
using Printf
using JSON
using DataStructures
import Base: getindex, time
import GZip
mutable struct Bus
name::String
offset::Int
load::Array{Float64}
units::Array
price_sensitive_loads::Array
end
mutable struct CostSegment
mw::Array{Float64}
cost::Array{Float64}
end
mutable struct StartupCategory
delay::Int
cost::Float64
end
mutable struct Unit
name::String
bus::Bus
max_power::Array{Float64}
min_power::Array{Float64}
must_run::Array{Bool}
min_power_cost::Array{Float64}
cost_segments::Array{CostSegment}
min_uptime::Int
min_downtime::Int
ramp_up_limit::Float64
ramp_down_limit::Float64
startup_limit::Float64
shutdown_limit::Float64
initial_status::Union{Int,Nothing}
initial_power::Union{Float64,Nothing}
provides_spinning_reserves::Array{Bool}
startup_categories::Array{StartupCategory}
end
mutable struct TransmissionLine
name::String
offset::Int
source::Bus
target::Bus
reactance::Float64
susceptance::Float64
normal_flow_limit::Array{Float64}
emergency_flow_limit::Array{Float64}
flow_limit_penalty::Array{Float64}
end
mutable struct Reserves
spinning::Array{Float64}
end
mutable struct Contingency
name::String
lines::Array{TransmissionLine}
units::Array{Unit}
end
mutable struct PriceSensitiveLoad
name::String
bus::Bus
demand::Array{Float64}
revenue::Array{Float64}
end
mutable struct UnitCommitmentInstance
time::Int
power_balance_penalty::Array{Float64}
units::Array{Unit}
buses::Array{Bus}
lines::Array{TransmissionLine}
reserves::Reserves
contingencies::Array{Contingency}
price_sensitive_loads::Array{PriceSensitiveLoad}
end
function Base.show(io::IO, instance::UnitCommitmentInstance)
print(io, "UnitCommitmentInstance with ")
print(io, "$(length(instance.units)) units, ")
print(io, "$(length(instance.buses)) buses, ")
print(io, "$(length(instance.lines)) lines, ")
print(io, "$(length(instance.contingencies)) contingencies, ")
print(io, "$(length(instance.price_sensitive_loads)) price sensitive loads")
end
function read_benchmark(name::AbstractString) :: UnitCommitmentInstance
basedir = dirname(@__FILE__)
return UnitCommitment.read("$basedir/../instances/$name.json.gz")
end
function read(path::AbstractString)::UnitCommitmentInstance
if endswith(path, ".gz")
return read(GZip.gzopen(path))
else
return read(open(path))
end
end
function read(file::IO)::UnitCommitmentInstance
return from_json(JSON.parse(file, dicttype=()->DefaultOrderedDict(nothing)))
end
function from_json(json; fix=true)
units = Unit[]
buses = Bus[]
contingencies = Contingency[]
lines = TransmissionLine[]
loads = PriceSensitiveLoad[]
T = json["Parameters"]["Time (h)"]
name_to_bus = Dict{String, Bus}()
name_to_line = Dict{String, TransmissionLine}()
name_to_unit = Dict{String, Unit}()
function timeseries(x; default=nothing)
x != nothing || return default
x isa Array || return [x for t in 1:T]
return x
end
function scalar(x; default=nothing)
x != nothing || return default
x
end
# Read parameters
power_balance_penalty = timeseries(json["Parameters"]["Power balance penalty (\$/MW)"],
default=[1000.0 for t in 1:T])
# Read buses
for (bus_name, dict) in json["Buses"]
bus = Bus(bus_name,
length(buses),
timeseries(dict["Load (MW)"]),
Unit[],
PriceSensitiveLoad[])
name_to_bus[bus_name] = bus
push!(buses, bus)
end
# Read units
for (unit_name, dict) in json["Generators"]
bus = name_to_bus[dict["Bus"]]
# Read production cost curve
K = length(dict["Production cost curve (MW)"])
curve_mw = hcat([timeseries(dict["Production cost curve (MW)"][k]) for k in 1:K]...)
curve_cost = hcat([timeseries(dict["Production cost curve (\$)"][k]) for k in 1:K]...)
min_power = curve_mw[:, 1]
max_power = curve_mw[:, K]
min_power_cost = curve_cost[:, 1]
segments = CostSegment[]
for k in 2:K
amount = curve_mw[:, k] - curve_mw[:, k-1]
cost = (curve_cost[:, k] - curve_cost[:, k-1]) ./ amount
replace!(cost, NaN=>0.0)
push!(segments, CostSegment(amount, cost))
end
# Read startup costs
startup_delays = scalar(dict["Startup delays (h)"], default=[1])
startup_costs = scalar(dict["Startup costs (\$)"], default=[0.])
startup_categories = StartupCategory[]
for k in 1:length(startup_delays)
push!(startup_categories, StartupCategory(startup_delays[k],
startup_costs[k]))
end
# Read and validate initial conditions
initial_power = scalar(dict["Initial power (MW)"], default=nothing)
initial_status = scalar(dict["Initial status (h)"], default=nothing)
if initial_power == nothing
initial_status == nothing || error("unit $unit_name has initial status but no initial power")
else
initial_status != nothing || error("unit $unit_name has initial power but no initial status")
initial_status != 0 || error("unit $unit_name has invalid initial status")
if initial_status < 0 && initial_power > 1e-3
error("unit $unit_name has invalid initial power")
end
end
unit = Unit(unit_name,
bus,
max_power,
min_power,
timeseries(dict["Must run?"], default=[false for t in 1:T]),
min_power_cost,
segments,
scalar(dict["Minimum uptime (h)"], default=1),
scalar(dict["Minimum downtime (h)"], default=1),
scalar(dict["Ramp up limit (MW)"], default=1e6),
scalar(dict["Ramp down limit (MW)"], default=1e6),
scalar(dict["Startup limit (MW)"], default=1e6),
scalar(dict["Shutdown limit (MW)"], default=1e6),
initial_status,
initial_power,
timeseries(dict["Provides spinning reserves?"],
default=[true for t in 1:T]),
startup_categories)
push!(bus.units, unit)
name_to_unit[unit_name] = unit
push!(units, unit)
end
# Read reserves
reserves = Reserves(zeros(T))
if "Reserves" in keys(json)
reserves.spinning = timeseries(json["Reserves"]["Spinning (MW)"],
default=zeros(T))
end
# Read transmission lines
if "Transmission lines" in keys(json)
for (line_name, dict) in json["Transmission lines"]
line = TransmissionLine(line_name,
length(lines) + 1,
name_to_bus[dict["Source bus"]],
name_to_bus[dict["Target bus"]],
scalar(dict["Reactance (ohms)"]),
scalar(dict["Susceptance (S)"]),
timeseries(dict["Normal flow limit (MW)"],
default=[1e8 for t in 1:T]),
timeseries(dict["Emergency flow limit (MW)"],
default=[1e8 for t in 1:T]),
timeseries(dict["Flow limit penalty (\$/MW)"],
default=[5000.0 for t in 1:T]))
name_to_line[line_name] = line
push!(lines, line)
end
end
# Read contingencies
if "Contingencies" in keys(json)
for (cont_name, dict) in json["Contingencies"]
affected_units = Unit[]
affected_lines = TransmissionLine[]
if "Affected lines" in keys(dict)
affected_lines = [name_to_line[l] for l in dict["Affected lines"]]
end
if "Affected units" in keys(dict)
affected_units = [name_to_unit[u] for u in dict["Affected units"]]
end
cont = Contingency(cont_name, affected_lines, affected_units)
push!(contingencies, cont)
end
end
# Read price-sensitive loads
if "Price-sensitive loads" in keys(json)
for (load_name, dict) in json["Price-sensitive loads"]
bus = name_to_bus[dict["Bus"]]
load = PriceSensitiveLoad(load_name,
bus,
timeseries(dict["Demand (MW)"]),
timeseries(dict["Revenue (\$/MW)"]),
)
push!(bus.price_sensitive_loads, load)
push!(loads, load)
end
end
instance = UnitCommitmentInstance(T,
power_balance_penalty,
units,
buses,
lines,
reserves,
contingencies,
loads)
if fix
UnitCommitment.fix!(instance)
end
return instance
end
"""
slice(instance, range)
Creates a new instance, with only a subset of the time periods.
This function does not modify the provided instance. The initial
conditions are also not modified.
Example
-------
# Build a 2-hour UC instance
instance = UnitCommitment.read_benchmark("test/case14")
modified = UnitCommitment.slice(instance, 1:2)
"""
function slice(instance::UnitCommitmentInstance, range::UnitRange{Int})::UnitCommitmentInstance
modified = deepcopy(instance)
modified.time = length(range)
modified.power_balance_penalty = modified.power_balance_penalty[range]
modified.reserves.spinning = modified.reserves.spinning[range]
for u in modified.units
u.max_power = u.max_power[range]
u.min_power = u.min_power[range]
u.must_run = u.must_run[range]
u.min_power_cost = u.min_power_cost[range]
u.provides_spinning_reserves = u.provides_spinning_reserves[range]
for s in u.cost_segments
s.mw = s.mw[range]
s.cost = s.cost[range]
end
end
for b in modified.buses
b.load = b.load[range]
end
for l in modified.lines
l.normal_flow_limit = l.normal_flow_limit[range]
l.emergency_flow_limit = l.emergency_flow_limit[range]
l.flow_limit_penalty = l.flow_limit_penalty[range]
end
for ps in modified.price_sensitive_loads
ps.demand = ps.demand[range]
ps.revenue = ps.revenue[range]
end
return modified
end
export UnitCommitmentInstance

50
src/log.jl Normal file
View File

@@ -0,0 +1,50 @@
# 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.
import Logging: min_enabled_level, shouldlog, handle_message
using Base.CoreLogging, Logging, Printf
struct TimeLogger <: AbstractLogger
initial_time::Float64
file::Union{Nothing, IOStream}
screen_log_level
io_log_level
end
function TimeLogger(;
initial_time::Float64,
file::Union{Nothing, IOStream} = nothing,
screen_log_level = CoreLogging.Info,
io_log_level = CoreLogging.Info,
) :: TimeLogger
return TimeLogger(initial_time, file, screen_log_level, io_log_level)
end
min_enabled_level(logger::TimeLogger) = logger.io_log_level
shouldlog(logger::TimeLogger, level, _module, group, id) = true
function handle_message(logger::TimeLogger,
level,
message,
_module,
group,
id,
filepath,
line;
kwargs...)
elapsed_time = time() - logger.initial_time
time_string = @sprintf("[%12.3f] ", elapsed_time)
if level >= logger.screen_log_level
print(time_string)
println(message)
end
if logger.file != nothing && level >= logger.io_log_level
write(logger.file, time_string)
write(logger.file, message)
write(logger.file, "\n")
flush(logger.file)
end
end
export TimeLogger

646
src/model.jl Normal file
View File

@@ -0,0 +1,646 @@
# 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.
using JuMP, MathOptInterface, DataStructures
import JuMP: value, fix, set_name
# Extend some JuMP functions so that decision variables can be safely replaced by
# (constant) floating point numbers.
function value(x::Float64)
x
end
function fix(x::Float64, v::Float64; force)
abs(x - v) < 1e-6 || error("Value mismatch: $x != $v")
end
function set_name(x::Float64, n::String)
# nop
end
mutable struct UnitCommitmentModel
mip::JuMP.Model
vars::DotDict
eqs::DotDict
exprs::DotDict
instance::UnitCommitmentInstance
isf::Array{Float64, 2}
lodf::Array{Float64, 2}
obj::AffExpr
end
function build_model(;
filename::Union{String, Nothing}=nothing,
instance::Union{UnitCommitmentInstance, Nothing}=nothing,
isf::Union{Array{Float64,2}, Nothing}=nothing,
lodf::Union{Array{Float64,2}, Nothing}=nothing,
isf_cutoff::Float64=0.005,
lodf_cutoff::Float64=0.001,
optimizer=nothing,
model=nothing,
variable_names::Bool=false,
) :: UnitCommitmentModel
if (filename == nothing) && (instance == nothing)
error("Either filename or instance must be specified")
end
if filename != nothing
@info "Reading: $(filename)"
time_read = @elapsed begin
instance = UnitCommitment.read(filename)
end
@info @sprintf("Read problem in %.2f seconds", time_read)
end
if length(instance.buses) == 1
isf = zeros(0, 0)
lodf = zeros(0, 0)
else
if isf == nothing
@info "Computing injection shift factors..."
time_isf = @elapsed begin
isf = UnitCommitment.injection_shift_factors(lines=instance.lines,
buses=instance.buses)
end
@info @sprintf("Computed ISF in %.2f seconds", time_isf)
@info "Computing line outage factors..."
time_lodf = @elapsed begin
lodf = UnitCommitment.line_outage_factors(lines=instance.lines,
buses=instance.buses,
isf=isf)
end
@info @sprintf("Computed LODF in %.2f seconds", time_lodf)
@info @sprintf("Applying PTDF and LODF cutoffs (%.5f, %.5f)", isf_cutoff, lodf_cutoff)
isf[abs.(isf) .< isf_cutoff] .= 0
lodf[abs.(lodf) .< lodf_cutoff] .= 0
end
end
@info "Building model..."
time_model = @elapsed begin
if model == nothing
if optimizer == nothing
mip = Model()
else
mip = Model(optimizer)
end
else
mip = model
end
model = UnitCommitmentModel(mip,
DotDict(), # vars
DotDict(), # eqs
DotDict(), # exprs
instance,
isf,
lodf,
AffExpr(), # obj
)
for field in [:prod_above, :segprod, :reserve, :is_on, :switch_on, :switch_off,
:net_injection, :curtail, :overflow, :loads, :startup]
setproperty!(model.vars, field, OrderedDict())
end
for field in [:startup_choose, :startup_restrict, :segprod_limit, :prod_above_def,
:prod_limit, :binary_link, :switch_on_off, :ramp_up, :ramp_down,
:startup_limit, :shutdown_limit, :min_uptime, :min_downtime, :power_balance,
:net_injection_def, :min_reserve]
setproperty!(model.eqs, field, OrderedDict())
end
for field in [:inj, :reserve, :net_injection]
setproperty!(model.exprs, field, OrderedDict())
end
for lm in instance.lines
add_transmission_line!(model, lm)
end
for b in instance.buses
add_bus!(model, b)
end
for g in instance.units
add_unit!(model, g)
end
for ps in instance.price_sensitive_loads
add_price_sensitive_load!(model, ps)
end
build_net_injection_eqs!(model)
build_reserve_eqs!(model)
build_obj_function!(model)
end
@info @sprintf("Built model in %.2f seconds", time_model)
if variable_names
set_variable_names!(model)
end
return model
end
function add_transmission_line!(model, lm)
vars, obj, T = model.vars, model.obj, model.instance.time
for t in 1:T
overflow = vars.overflow[lm.name, t] = @variable(model.mip, lower_bound=0)
add_to_expression!(obj, overflow, lm.flow_limit_penalty[t])
end
end
function add_bus!(model::UnitCommitmentModel, b::Bus)
mip, vars, exprs = model.mip, model.vars, model.exprs
for t in 1:model.instance.time
# Fixed load
exprs.net_injection[b.name, t] = AffExpr(-b.load[t])
# Reserves
exprs.reserve[b.name, t] = AffExpr()
# Load curtailment
vars.curtail[b.name, t] = @variable(mip, lower_bound=0, upper_bound=b.load[t])
add_to_expression!(exprs.net_injection[b.name, t], vars.curtail[b.name, t], 1.0)
add_to_expression!(model.obj,
vars.curtail[b.name, t],
model.instance.power_balance_penalty[t])
end
end
function add_price_sensitive_load!(model::UnitCommitmentModel, ps::PriceSensitiveLoad)
mip, vars = model.mip, model.vars
for t in 1:model.instance.time
# Decision variable
vars.loads[ps.name, t] = @variable(mip, lower_bound=0, upper_bound=ps.demand[t])
# Objective function terms
add_to_expression!(model.obj, vars.loads[ps.name, t], -ps.revenue[t])
# Net injection
add_to_expression!(model.exprs.net_injection[ps.bus.name, t], vars.loads[ps.name, t], -1.0)
end
end
function add_unit!(model::UnitCommitmentModel, g::Unit)
mip, vars, eqs, exprs, T = model.mip, model.vars, model.eqs, model.exprs, model.instance.time
gi, K, S = g.name, length(g.cost_segments), length(g.startup_categories)
if !all(g.must_run) && any(g.must_run)
error("Partially must-run units are not currently supported")
end
if g.initial_power == nothing || g.initial_status == nothing
error("Initial conditions for $(g.name) must be provided")
end
is_initially_on = (g.initial_status > 0 ? 1.0 : 0.0)
# Decision variables
for t in 1:T
for k in 1:K
model.vars.segprod[gi, t, k] = @variable(model.mip, lower_bound=0)
end
model.vars.prod_above[gi, t] = @variable(model.mip, lower_bound=0)
if g.provides_spinning_reserves[t]
model.vars.reserve[gi, t] = @variable(model.mip, lower_bound=0)
else
model.vars.reserve[gi, t] = 0.0
end
for s in 1:S
model.vars.startup[gi, t, s] = @variable(model.mip, binary=true)
end
if g.must_run[t]
model.vars.is_on[gi, t] = 1.0
model.vars.switch_on[gi, t] = (t == 1 ? 1.0 - is_initially_on : 0.0)
model.vars.switch_off[gi, t] = 0.0
else
model.vars.is_on[gi, t] = @variable(model.mip, binary=true)
model.vars.switch_on[gi, t] = @variable(model.mip, binary=true)
model.vars.switch_off[gi, t] = @variable(model.mip, binary=true)
end
end
for t in 1:T
# Time-dependent start-up costs
for s in 1:S
# If unit is switching on, we must choose a startup category
eqs.startup_choose[gi, t, s] =
@constraint(mip, vars.switch_on[gi, t] == sum(vars.startup[gi, 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
range = (t - g.startup_categories[s + 1].delay + 1):(t - g.startup_categories[s].delay)
initial_sum = (g.initial_status < 0 && (g.initial_status + 1 in range) ? 1.0 : 0.0)
eqs.startup_restrict[gi, t, s] =
@constraint(mip, vars.startup[gi, t, s]
<= initial_sum + sum(vars.switch_off[gi, i] for i in range if i >= 1))
end
# Objective function terms for start-up costs
add_to_expression!(model.obj,
vars.startup[gi, t, s],
g.startup_categories[s].cost)
end
# Objective function terms for production costs
add_to_expression!(model.obj, vars.is_on[gi, t], g.min_power_cost[t])
for k in 1:K
add_to_expression!(model.obj, vars.segprod[gi, t, k], g.cost_segments[k].cost[t])
end
# Production limits (piecewise-linear segments)
for k in 1:K
eqs.segprod_limit[gi, t, k] =
@constraint(mip, vars.segprod[gi, t, k] <= g.cost_segments[k].mw[t] * vars.is_on[gi, t])
end
# Definition of production
eqs.prod_above_def[gi, t] =
@constraint(mip, vars.prod_above[gi, t] == sum(vars.segprod[gi, t, k] for k in 1:K))
# Production limit
eqs.prod_limit[gi, t] =
@constraint(mip,
vars.prod_above[gi, t] + vars.reserve[gi, t]
<= (g.max_power[t] - g.min_power[t]) * vars.is_on[gi, t])
# Binary variable equations for economic units
if !g.must_run[t]
# Link binary variables
if t == 1
eqs.binary_link[gi, t] =
@constraint(mip,
vars.is_on[gi, t] - is_initially_on ==
vars.switch_on[gi, t] - vars.switch_off[gi, t])
else
eqs.binary_link[gi, t] =
@constraint(mip,
vars.is_on[gi, t] - vars.is_on[gi, t-1] ==
vars.switch_on[gi, t] - vars.switch_off[gi, t])
end
# Cannot switch on and off at the same time
eqs.switch_on_off[gi, t] =
@constraint(mip, vars.switch_on[gi, t] + vars.switch_off[gi, t] <= 1)
end
# Ramp up limit
if t == 1
if is_initially_on == 1
eqs.ramp_up[gi, t] =
@constraint(mip,
vars.prod_above[gi, t] + vars.reserve[gi, t] <=
(g.initial_power - g.min_power[t]) + g.ramp_up_limit)
end
else
eqs.ramp_up[gi, t] =
@constraint(mip,
vars.prod_above[gi, t] + vars.reserve[gi, t] <=
vars.prod_above[gi, t-1] + g.ramp_up_limit)
end
# Ramp down limit
if t == 1
if is_initially_on == 1
eqs.ramp_down[gi, t] =
@constraint(mip,
vars.prod_above[gi, t] >=
(g.initial_power - g.min_power[t]) - g.ramp_down_limit)
end
else
eqs.ramp_down[gi, t] =
@constraint(mip,
vars.prod_above[gi, t] >=
vars.prod_above[gi, t-1] - g.ramp_down_limit)
end
# Startup limit
eqs.startup_limit[gi, t] =
@constraint(mip,
vars.prod_above[gi, t] + vars.reserve[gi, t] <=
(g.max_power[t] - g.min_power[t]) * vars.is_on[gi, t]
- max(0, g.max_power[t] - g.startup_limit) * vars.switch_on[gi, t])
# Shutdown limit
if g.initial_power > g.shutdown_limit
eqs.shutdown_limit[gi, 0] =
@constraint(mip, vars.switch_off[gi, 1] <= 0)
end
if t < T
eqs.shutdown_limit[gi, t] =
@constraint(mip,
vars.prod_above[gi, t] <=
(g.max_power[t] - g.min_power[t]) * vars.is_on[gi, t]
- max(0, g.max_power[t] - g.shutdown_limit) * vars.switch_off[gi, t+1])
end
# Minimum up-time
eqs.min_uptime[gi, t] =
@constraint(mip,
sum(vars.switch_on[gi, i]
for i in (t - g.min_uptime + 1):t if i >= 1
) <= vars.is_on[gi, t])
# # Minimum down-time
eqs.min_downtime[gi, t] =
@constraint(mip,
sum(vars.switch_off[gi, i]
for i in (t - g.min_downtime + 1):t if i >= 1
) <= 1 - vars.is_on[gi, t])
# Minimum up/down-time for initial periods
if t == 1
if g.initial_status > 0
eqs.min_uptime[gi, 0] =
@constraint(mip, sum(vars.switch_off[gi, i]
for i in 1:(g.min_uptime - g.initial_status) if i <= T) == 0)
else
eqs.min_downtime[gi, 0] =
@constraint(mip, sum(vars.switch_on[gi, i]
for i in 1:(g.min_downtime + g.initial_status) if i <= T) == 0)
end
end
# Add to net injection expression
add_to_expression!(exprs.net_injection[g.bus.name, t], vars.prod_above[g.name, t], 1.0)
add_to_expression!(exprs.net_injection[g.bus.name, t], vars.is_on[g.name, t], g.min_power[t])
# Add to reserves expression
add_to_expression!(exprs.reserve[g.bus.name, t], vars.reserve[gi, t], 1.0)
end
end
function build_obj_function!(model::UnitCommitmentModel)
@objective(model.mip, Min, model.obj)
end
function build_net_injection_eqs!(model::UnitCommitmentModel)
T = model.instance.time
for t in 1:T, b in model.instance.buses
net = model.vars.net_injection[b.name, t] = @variable(model.mip)
model.eqs.net_injection_def[t, b.name] =
@constraint(model.mip, net == model.exprs.net_injection[b.name, t])
end
for t in 1:T
model.eqs.power_balance[t] =
@constraint(model.mip, sum(model.vars.net_injection[b.name, t]
for b in model.instance.buses) == 0)
end
end
function build_reserve_eqs!(model::UnitCommitmentModel)
reserves = model.instance.reserves
for t in 1:model.instance.time
model.eqs.min_reserve[t] =
@constraint(model.mip, sum(model.exprs.reserve[b.name, t]
for b in model.instance.buses) >= reserves.spinning[t])
end
end
function enforce_transmission(;
model::UnitCommitmentModel,
violation::Violation,
isf::Array{Float64,2},
lodf::Array{Float64,2})::Nothing
instance, mip, vars = model.instance, model.mip, model.vars
limit::Float64 = 0.0
if violation.outage_line == nothing
limit = violation.monitored_line.normal_flow_limit[violation.time]
@info @sprintf(" %8.3f MW overflow in %-5s time %3d (pre-contingency)",
violation.amount,
violation.monitored_line.name,
violation.time)
else
limit = violation.monitored_line.emergency_flow_limit[violation.time]
@info @sprintf(" %8.3f MW overflow in %-5s time %3d (outage: line %s)",
violation.amount,
violation.monitored_line.name,
violation.time,
violation.outage_line.name)
end
fm = violation.monitored_line.name
t = violation.time
flow = @variable(mip, base_name="flow[$fm,$t]")
overflow = vars.overflow[violation.monitored_line.name, violation.time]
@constraint(mip, flow <= limit + overflow)
@constraint(mip, -flow <= limit + overflow)
if violation.outage_line == nothing
@constraint(mip, flow == sum(vars.net_injection[b.name, violation.time] *
isf[violation.monitored_line.offset, b.offset]
for b in instance.buses
if b.offset > 0))
else
@constraint(mip, flow == sum(vars.net_injection[b.name, violation.time] * (
isf[violation.monitored_line.offset, b.offset] + (
lodf[violation.monitored_line.offset, violation.outage_line.offset] *
isf[violation.outage_line.offset, b.offset]
)
)
for b in instance.buses
if b.offset > 0))
end
nothing
end
function set_variable_names!(model::UnitCommitmentModel)
@info "Setting variable and constraint names..."
time_varnames = @elapsed begin
set_jump_names!(model.vars)
set_jump_names!(model.eqs)
end
@info @sprintf("Set names in %.2f seconds", time_varnames)
end
function set_jump_names!(dict)
for name in keys(dict)
for idx in keys(dict[name])
idx_str = join(map(string, idx), ",")
set_name(dict[name][idx], "$name[$idx_str]")
end
end
end
function get_solution(model::UnitCommitmentModel)
instance, T = model.instance, model.instance.time
function timeseries(vars, collection)
return OrderedDict(b.name => [round(value(vars[b.name, t]), digits=5) for t in 1:T]
for b in collection)
end
function production_cost(g)
return [value(model.vars.is_on[g.name, t]) * g.min_power_cost[t] +
sum(Float64[value(model.vars.segprod[g.name, t, k]) * g.cost_segments[k].cost[t]
for k in 1:length(g.cost_segments)])
for t in 1:T]
end
function production(g)
return [value(model.vars.is_on[g.name, t]) * g.min_power[t] +
sum(Float64[value(model.vars.segprod[g.name, t, k])
for k in 1:length(g.cost_segments)])
for t in 1:T]
end
function startup_cost(g)
S = length(g.startup_categories)
return [sum(g.startup_categories[s].cost * value(model.vars.startup[g.name, t, s])
for s in 1:S)
for t in 1:T]
end
sol = OrderedDict()
sol["Production (MW)"] = OrderedDict(g.name => production(g) for g in instance.units)
sol["Production cost (\$)"] = OrderedDict(g.name => production_cost(g) for g in instance.units)
sol["Startup cost (\$)"] = OrderedDict(g.name => startup_cost(g) for g in instance.units)
sol["Is on"] = timeseries(model.vars.is_on, instance.units)
sol["Switch on"] = timeseries(model.vars.switch_on, instance.units)
sol["Switch off"] = timeseries(model.vars.switch_off, instance.units)
sol["Reserve (MW)"] = timeseries(model.vars.reserve, instance.units)
sol["Net injection (MW)"] = timeseries(model.vars.net_injection, instance.buses)
sol["Load curtail (MW)"] = timeseries(model.vars.curtail, instance.buses)
if !isempty(instance.lines)
sol["Line overflow (MW)"] = timeseries(model.vars.overflow, instance.lines)
end
if !isempty(instance.price_sensitive_loads)
sol["Price-sensitive loads (MW)"] = timeseries(model.vars.loads, instance.price_sensitive_loads)
end
return sol
end
function fix!(model::UnitCommitmentModel, solution)::Nothing
vars, instance, T = model.vars, model.instance, model.instance.time
for g in instance.units
for t in 1:T
is_on = round(solution["Is on"][g.name][t])
production = round(solution["Production (MW)"][g.name][t], digits=5)
reserve = round(solution["Reserve (MW)"][g.name][t], digits=5)
JuMP.fix(vars.is_on[g.name, t], is_on, force=true)
JuMP.fix(vars.prod_above[g.name, t], production - is_on * g.min_power[t], force=true)
JuMP.fix(vars.reserve[g.name, t], reserve, force=true)
end
end
end
function set_warm_start!(model::UnitCommitmentModel, solution)::Nothing
vars, instance, T = model.vars, model.instance, model.instance.time
for g in instance.units
for t in 1:T
JuMP.set_start_value(vars.is_on[g.name, t], solution["Is on"][g.name][t])
JuMP.set_start_value(vars.switch_on[g.name, t], solution["Switch on"][g.name][t])
JuMP.set_start_value(vars.switch_off[g.name, t], solution["Switch off"][g.name][t])
end
end
end
function optimize!(model::UnitCommitmentModel;
time_limit=3600,
gap_limit=1e-4,
two_phase_gap=true,
)::Nothing
function set_gap(gap)
try
JuMP.set_optimizer_attribute(model.mip, "MIPGap", gap)
@info @sprintf("MIP gap tolerance set to %f", gap)
catch
@warn "Could not change MIP gap tolerance"
end
end
instance = model.instance
initial_time = time()
large_gap = false
has_transmission = (length(model.isf) > 0)
if has_transmission && two_phase_gap
set_gap(1e-2)
large_gap = true
else
set_gap(gap_limit)
end
while true
time_elapsed = time() - initial_time
time_remaining = time_limit - time_elapsed
if time_remaining < 0
@info "Time limit exceeded"
break
end
@info @sprintf("Setting MILP time limit to %.2f seconds", time_remaining)
JuMP.set_time_limit_sec(model.mip, time_remaining)
@info "Solving MILP..."
JuMP.optimize!(model.mip)
has_transmission || break
violations = find_violations(model)
if isempty(violations)
@info "No violations found"
if large_gap
large_gap = false
set_gap(gap_limit)
else
break
end
else
enforce_transmission(model, violations)
end
end
nothing
end
function find_violations(model::UnitCommitmentModel)
instance, vars = model.instance, model.vars
length(instance.buses) > 1 || return []
violations = []
@info "Verifying transmission limits..."
time_screening = @elapsed begin
non_slack_buses = [b for b in instance.buses if b.offset > 0]
net_injections = [value(vars.net_injection[b.name, t])
for b in non_slack_buses, t in 1:instance.time]
overflow = [value(vars.overflow[lm.name, t])
for lm in instance.lines, t in 1:instance.time]
violations = UnitCommitment.find_violations(instance=instance,
net_injections=net_injections,
overflow=overflow,
isf=model.isf,
lodf=model.lodf)
end
@info @sprintf("Verified transmission limits in %.2f seconds", time_screening)
return violations
end
function enforce_transmission(model::UnitCommitmentModel, violations::Array{Violation, 1})
for v in violations
enforce_transmission(model=model,
violation=v,
isf=model.isf,
lodf=model.lodf)
end
end
export UnitCommitmentModel, build_model, get_solution, optimize!

198
src/screening.jl Normal file
View File

@@ -0,0 +1,198 @@
# 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.
# Copyright (C) 2019 Argonne National Laboratory
# Written by Alinson Santos Xavier <axavier@anl.gov>
using DataStructures
using Base.Threads
struct Violation
time::Int
monitored_line::TransmissionLine
outage_line::Union{TransmissionLine, Nothing}
amount::Float64 # Violation amount (in MW)
end
function Violation(;
time::Int,
monitored_line::TransmissionLine,
outage_line::Union{TransmissionLine, Nothing},
amount::Float64,
) :: Violation
return Violation(time, monitored_line, outage_line, amount)
end
mutable struct ViolationFilter
max_per_line::Int
max_total::Int
queues::Dict{Int, PriorityQueue{Violation, Float64}}
end
function ViolationFilter(;
max_per_line::Int=1,
max_total::Int=5,
)::ViolationFilter
return ViolationFilter(max_per_line, max_total, Dict())
end
function offer(filter::ViolationFilter, v::Violation)::Nothing
if v.monitored_line.offset ∉ keys(filter.queues)
filter.queues[v.monitored_line.offset] = PriorityQueue{Violation, Float64}()
end
q::PriorityQueue{Violation, Float64} = filter.queues[v.monitored_line.offset]
if length(q) < filter.max_per_line
enqueue!(q, v => v.amount)
else
if v.amount > peek(q)[1].amount
dequeue!(q)
enqueue!(q, v => v.amount)
end
end
nothing
end
function query(filter::ViolationFilter)::Array{Violation, 1}
violations = Array{Violation,1}()
time_queue = PriorityQueue{Violation, Float64}()
for l in keys(filter.queues)
line_queue = filter.queues[l]
while length(line_queue) > 0
v = dequeue!(line_queue)
if length(time_queue) < filter.max_total
enqueue!(time_queue, v => v.amount)
else
if v.amount > peek(time_queue)[1].amount
dequeue!(time_queue)
enqueue!(time_queue, v => v.amount)
end
end
end
end
while length(time_queue) > 0
violations = [violations; dequeue!(time_queue)]
end
return violations
end
"""
function find_violations(instance::UnitCommitmentInstance,
net_injections::Array{Float64, 2};
isf::Array{Float64,2},
lodf::Array{Float64,2},
max_per_line::Int = 1,
max_per_period::Int = 5,
) :: Array{Violation, 1}
Find transmission constraint violations (both pre-contingency, as well as post-contingency).
The argument `net_injection` should be a (B-1) x T matrix, where B is the number of buses
and T is the number of time periods. The arguments `isf` and `lodf` can be computed using
UnitCommitment.injection_shift_factors and UnitCommitment.line_outage_factors.
The argument `overflow` specifies how much flow above the transmission limits (in MW) is allowed.
It should be an L x T matrix, where L is the number of transmission lines.
"""
function find_violations(;
instance::UnitCommitmentInstance,
net_injections::Array{Float64, 2},
overflow::Array{Float64, 2},
isf::Array{Float64,2},
lodf::Array{Float64,2},
max_per_line::Int = 1,
max_per_period::Int = 5,
)::Array{Violation, 1}
B = length(instance.buses) - 1
L = length(instance.lines)
T = instance.time
K = nthreads()
size(net_injections) == (B, T) || error("net_injections has incorrect size")
size(isf) == (L, B) || error("isf has incorrect size")
size(lodf) == (L, L) || error("lodf has incorrect size")
filters = Dict(t => ViolationFilter(max_total=max_per_period,
max_per_line=max_per_line)
for t in 1:T)
pre_flow::Array{Float64} = zeros(L, K) # pre_flow[lm, thread]
post_flow::Array{Float64} = zeros(L, L, K) # post_flow[lm, lc, thread]
pre_v::Array{Float64} = zeros(L, K) # pre_v[lm, thread]
post_v::Array{Float64} = zeros(L, L, K) # post_v[lm, lc, thread]
normal_limits::Array{Float64,2} = [l.normal_flow_limit[t] + overflow[l.offset, t]
for l in instance.lines, t in 1:T]
emergency_limits::Array{Float64,2} = [l.emergency_flow_limit[t] + overflow[l.offset, t]
for l in instance.lines, t in 1:T]
is_vulnerable::Array{Bool} = zeros(Bool, L)
for c in instance.contingencies
is_vulnerable[c.lines[1].offset] = true
end
@threads for t in 1:T
k = threadid()
# Pre-contingency flows
pre_flow[:, k] = isf * net_injections[:, t]
# Post-contingency flows
for lc in 1:L, lm in 1:L
post_flow[lm, lc, k] = pre_flow[lm, k] + pre_flow[lc, k] * lodf[lm, lc]
end
# Pre-contingency violations
for lm in 1:L
pre_v[lm, k] = max(0.0,
pre_flow[lm, k] - normal_limits[lm, t],
- pre_flow[lm, k] - normal_limits[lm, t])
end
# Post-contingency violations
for lc in 1:L, lm in 1:L
post_v[lm, lc, k] = max(0.0,
post_flow[lm, lc, k] - emergency_limits[lm, t],
- post_flow[lm, lc, k] - emergency_limits[lm, t])
end
# Offer pre-contingency violations
for lm in 1:L
if pre_v[lm, k] > 1e-5
offer(filters[t], Violation(time=t,
monitored_line=instance.lines[lm],
outage_line=nothing,
amount=pre_v[lm, k]))
end
end
# Offer post-contingency violations
for lm in 1:L, lc in 1:L
if post_v[lm, lc, k] > 1e-5 && is_vulnerable[lc]
offer(filters[t], Violation(time=t,
monitored_line=instance.lines[lm],
outage_line=instance.lines[lc],
amount=post_v[lm, lc, k]))
end
end
end
violations = Violation[]
for t in 1:instance.time
append!(violations, query(filters[t]))
end
return violations
end
export Violation, ViolationFilter, offer, query, find_violations

80
src/sensitivity.jl Normal file
View File

@@ -0,0 +1,80 @@
# 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.
using SparseArrays, Base.Threads, LinearAlgebra, JuMP
"""
injection_shift_factors(; buses, lines)
Returns a (B-1)xL matrix M, where B is the number of buses and L is the number of transmission
lines. For a given bus b and transmission line l, the entry M[l.offset, b.offset] indicates
the amount of power (in MW) that flows through transmission line l when 1 MW of power is
injected at the slack bus (the bus that has offset zero) and withdrawn from b.
"""
function injection_shift_factors(; buses, lines)
susceptance = susceptance_matrix(lines)
incidence = reduced_incidence_matrix(lines = lines, buses = buses)
laplacian = transpose(incidence) * susceptance * incidence
isf = susceptance * incidence * inv(Array(laplacian))
return isf
end
"""
reduced_incidence_matrix(; buses::Array{Bus}, lines::Array{TransmissionLine})
Returns the incidence matrix for the network, with the column corresponding to the slack
bus is removed. More precisely, returns a (B-1) x L matrix, where B is the number of buses
and L is the number of lines. For each row, there is a 1 element and a -1 element, indicating
the source and target buses, respectively, for that line.
"""
function reduced_incidence_matrix(; buses::Array{Bus}, lines::Array{TransmissionLine})
matrix = spzeros(Float64, length(lines), length(buses) - 1)
for line in lines
if line.source.offset > 0
matrix[line.offset, line.source.offset] = 1
end
if line.target.offset > 0
matrix[line.offset, line.target.offset] = -1
end
end
matrix
end
"""
susceptance_matrix(lines::Array{TransmissionLine})
Returns a LxL diagonal matrix, where each diagonal entry is the susceptance of the
corresponding transmission line.
"""
function susceptance_matrix(lines::Array{TransmissionLine})
return Diagonal([l.susceptance for l in lines])
end
"""
line_outage_factors(; buses, lines, isf)
Returns a LxL matrix containing the Line Outage Distribution Factors (LODFs) for the
given network. This matrix how does the pre-contingency flow change when each individual
transmission line is removed.
"""
function line_outage_factors(;
buses::Array{Bus, 1},
lines::Array{TransmissionLine, 1},
isf::Array{Float64,2},
) :: Array{Float64,2}
n_lines, n_buses = size(isf)
incidence = Array(reduced_incidence_matrix(lines=lines,
buses=buses))
lodf::Array{Float64,2} = isf * transpose(incidence)
m, n = size(lodf)
for i in 1:n
lodf[:, i] *= 1.0 / (1.0 - lodf[i, i])
lodf[i, i] = -1
end
return lodf
end

26
src/sysimage.jl Normal file
View File

@@ -0,0 +1,26 @@
# 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.
using PackageCompiler
using DataStructures
using Documenter
using GLPK
using JSON
using JuMP
using MathOptInterface
using SparseArrays
using TimerOutputs
pkg = [:DataStructures,
:Documenter,
:GLPK,
:JSON,
:JuMP,
:MathOptInterface,
:SparseArrays,
:TimerOutputs]
@info "Building system image..."
create_sysimage(pkg, sysimage_path="build/sysimage.so")

334
src/validate.jl Normal file
View File

@@ -0,0 +1,334 @@
# 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.
using Printf
bin(x) = [xi > 0.5 for xi in x]
"""
fix!(instance)
Verifies that the given unit commitment instance is valid and automatically fixes
some validation errors if possible, issuing a warning for each error found.
If a validation error cannot be automatically fixed, issues an exception.
Returns the number of validation errors found.
"""
function fix!(instance::UnitCommitmentInstance)::Int
n_errors = 0
for g in instance.units
# Startup costs and delays must be increasing
for s in 2:length(g.startup_categories)
if g.startup_categories[s].delay <= g.startup_categories[s-1].delay
prev_value = g.startup_categories[s].delay
new_value = g.startup_categories[s-1].delay + 1
@warn "Generator $(g.name) has non-increasing startup delays (category $s). " *
"Changing delay: $prev_value → $new_value"
g.startup_categories[s].delay = new_value
n_errors += 1
end
if g.startup_categories[s].cost < g.startup_categories[s-1].cost
prev_value = g.startup_categories[s].cost
new_value = g.startup_categories[s-1].cost
@warn "Generator $(g.name) has decreasing startup cost (category $s). " *
"Changing cost: $prev_value → $new_value"
g.startup_categories[s].cost = new_value
n_errors += 1
end
end
for t in 1:instance.time
# Production cost curve should be convex
for k in 2:length(g.cost_segments)
cost = g.cost_segments[k].cost[t]
min_cost = g.cost_segments[k-1].cost[t]
if cost < min_cost - 1e-5
@warn "Generator $(g.name) has non-convex production cost curve " *
"(segment $k, time $t). Changing cost: $cost → $min_cost"
g.cost_segments[k].cost[t] = min_cost
n_errors += 1
end
end
# Startup limit must be greater than min_power
if g.startup_limit < g.min_power[t]
new_limit = g.min_power[t]
prev_limit = g.startup_limit
@warn "Generator $(g.name) has startup limit lower than minimum power. " *
"Changing startup limit: $prev_limit → $new_limit"
g.startup_limit = new_limit
n_errors += 1
end
end
end
return n_errors
end
function validate(instance_filename::String, solution_filename::String)
instance = UnitCommitment.read(instance_filename)
solution = JSON.parse(open(solution_filename))
return validate(instance, solution)
end
"""
validate(instance, solution)::Bool
Verifies that the given solution is feasible for the problem. If feasible,
silently returns true. In infeasible, returns false and prints the validation
errors to the screen.
This function is implemented independently from the optimization model in `model.jl`, and
therefore can be used to verify that the model is indeed producing valid solutions. It
can also be used to verify the solutions produced by other optimization packages.
"""
function validate(instance::UnitCommitmentInstance,
solution::Union{Dict,OrderedDict};
)::Bool
err_count = 0
err_count += validate_units(instance, solution)
err_count += validate_reserve_and_demand(instance, solution)
if err_count > 0
@error "Found $err_count validation errors"
return false
end
return true
end
function validate_units(instance, solution; tol=0.01)
err_count = 0
for unit in instance.units
production = solution["Production (MW)"][unit.name]
reserve = solution["Reserve (MW)"][unit.name]
actual_production_cost = solution["Production cost (\$)"][unit.name]
actual_startup_cost = solution["Startup cost (\$)"][unit.name]
is_on = bin(solution["Is on"][unit.name])
for t in 1:instance.time
# Auxiliary variables
if t == 1
is_starting_up = (unit.initial_status < 0) && is_on[t]
is_shutting_down = (unit.initial_status > 0) && !is_on[t]
ramp_up = max(0, production[t] + reserve[t] - unit.initial_power)
ramp_down = max(0, unit.initial_power - production[t])
else
is_starting_up = !is_on[t-1] && is_on[t]
is_shutting_down = is_on[t-1] && !is_on[t]
ramp_up = max(0, production[t] + reserve[t] - production[t-1])
ramp_down = max(0, production[t-1] - production[t])
end
# Compute production costs
production_cost, startup_cost = 0, 0
if is_on[t]
production_cost += unit.min_power_cost[t]
residual = max(0, production[t] - unit.min_power[t])
for s in unit.cost_segments
cleared = min(residual, s.mw[t])
production_cost += cleared * s.cost[t]
residual = max(0, residual - s.mw[t])
end
end
# Production should be non-negative
if production[t] < -tol
@error @sprintf("Unit %s produces negative amount of power at time %d (%.2f)",
unit.name, t, production[t])
err_count += 1
end
# Verify must-run
if !is_on[t] && unit.must_run[t]
@error @sprintf("Must-run unit %s is offline at time %d",
unit.name, t)
err_count += 1
end
# Verify reserve eligibility
if !unit.provides_spinning_reserves[t] && reserve[t] > tol
@error @sprintf("Unit %s is not eligible to provide spinning reserves at time %d",
unit.name, t)
err_count += 1
end
# If unit is on, must produce at least its minimum power
if is_on[t] && (production[t] < unit.min_power[t] - tol)
@error @sprintf("Unit %s produces below its minimum limit at time %d (%.2f < %.2f)",
unit.name, t, production[t], unit.min_power[t])
err_count += 1
end
# If unit is on, must produce at most its maximum power
if is_on[t] && (production[t] + reserve[t] > unit.max_power[t] + tol)
@error @sprintf("Unit %s produces above its maximum limit at time %d (%.2f + %.2f> %.2f)",
unit.name, t, production[t], reserve[t], unit.max_power[t])
err_count += 1
end
# If unit is off, must produce zero
if !is_on[t] && production[t] + reserve[t] > tol
@error @sprintf("Unit %s produces power at time %d while off",
unit.name, t)
err_count += 1
end
# Startup limit
if is_starting_up && (ramp_up > unit.startup_limit + tol)
@error @sprintf("Unit %s exceeds startup limit at time %d (%.2f > %.2f)",
unit.name, t, ramp_up, unit.startup_limit)
err_count += 1
end
# Shutdown limit
if is_shutting_down && (ramp_down > unit.shutdown_limit + tol)
@error @sprintf("Unit %s exceeds shutdown limit at time %d (%.2f > %.2f)",
unit.name, t, ramp_down, unit.shutdown_limit)
err_count += 1
end
# Ramp-up limit
if !is_starting_up && !is_shutting_down && (ramp_up > unit.ramp_up_limit + tol)
@error @sprintf("Unit %s exceeds ramp up limit at time %d (%.2f > %.2f)",
unit.name, t, ramp_up, unit.ramp_up_limit)
err_count += 1
end
# Ramp-down limit
if !is_starting_up && !is_shutting_down && (ramp_down > unit.ramp_down_limit + tol)
@error @sprintf("Unit %s exceeds ramp down limit at time %d (%.2f > %.2f)",
unit.name, t, ramp_down, unit.ramp_down_limit)
err_count += 1
end
# Verify startup costs & minimum downtime
if is_starting_up
# Calculate how much time the unit has been offline
time_down = 0
for k in 1:(t-1)
if !is_on[t - k]
time_down += 1
else
break
end
end
if t == time_down + 1
initial_down = unit.min_downtime
if unit.initial_status < 0
initial_down = -unit.initial_status
end
time_down += initial_down
end
# Calculate startup costs
for c in unit.startup_categories
if time_down >= c.delay
startup_cost = c.cost
end
end
# Check minimum downtime
if time_down < unit.min_downtime
@error @sprintf("Unit %s violates minimum downtime at time %d",
unit.name, t)
err_count += 1
end
end
# Verify minimum uptime
if is_shutting_down
# Calculate how much time the unit has been online
time_up = 0
for k in 1:(t-1)
if is_on[t - k]
time_up += 1
else
break
end
end
if t == time_up + 1
initial_up = unit.min_uptime
if unit.initial_status > 0
initial_up = unit.initial_status
end
time_up += initial_up
end
if (t == time_up + 1) && (unit.initial_status > 0)
time_up += unit.initial_status
end
# Check minimum uptime
if time_up < unit.min_uptime
@error @sprintf("Unit %s violates minimum uptime at time %d",
unit.name, t)
err_count += 1
end
end
# Verify production costs
if abs(actual_production_cost[t] - production_cost) > 1.00
@error @sprintf("Unit %s has unexpected production cost at time %d (%.2f should be %.2f)",
unit.name, t, actual_production_cost[t], production_cost)
err_count += 1
end
# Verify startup costs
if abs(actual_startup_cost[t] - startup_cost) > 1.00
@error @sprintf("Unit %s has unexpected startup cost at time %d (%.2f should be %.2f)",
unit.name, t, actual_startup_cost[t], startup_cost)
err_count += 1
end
end
end
return err_count
end
function validate_reserve_and_demand(instance, solution, tol=0.01)
err_count = 0
for t in 1:instance.time
load_curtail = 0
fixed_load = sum(b.load[t] for b in instance.buses)
production = sum(solution["Production (MW)"][g.name][t]
for g in instance.units)
if "Load curtail (MW)" in keys(solution)
load_curtail = sum(solution["Load curtail (MW)"][b.name][t]
for b in instance.buses)
end
balance = fixed_load - load_curtail - production
# Verify that production equals demand
if abs(balance) > tol
@error @sprintf("Non-zero power balance at time %d (%.2f - %.2f - %.2f != 0)",
t, fixed_load, load_curtail, production)
err_count += 1
end
# Verify spinning reserves
reserve = sum(solution["Reserve (MW)"][g.name][t] for g in instance.units)
if reserve < instance.reserves.spinning[t] - tol
@error @sprintf("Insufficient spinning reserves at time %d (%.2f should be %.2f)",
t, reserve, instance.reserves.spinning[t])
err_count += 1
end
end
return err_count
end