Getting started
Installing the package
UnitCommitment.jl was tested and developed with Julia 1.10. To install Julia, please follow the installation guide on the official Julia website. To install UnitCommitment.jl, run the Julia interpreter, type ]
to open the package manager, then type:
pkg> add UnitCommitment@0.4
To solve the optimization models, a mixed-integer linear programming (MILP) solver is also required. Please see the JuMP installation guide for more instructions on installing a solver. Typical open-source choices are HiGHS, Cbc and GLPK. In the instructions below, HiGHS will be used, but any other MILP solver should also be compatible.
Solving a benchmark instance
We start this tutorial by illustrating how to use UnitCommitment.jl to solve one of the provided benchmark instances. The package contains a large number of deterministic benchmark instances collected from the literature and converted into a common data format, which can be used to evaluate the performance of different solution methods. See Instances for more details. The first step is to import UnitCommitment
and HiGHS.
using HiGHS
using UnitCommitment
Next, we use the function UnitCommitment.read_benchmark
to read the instance.
instance = UnitCommitment.read_benchmark("matpower/case14/2017-01-01");
Now that we have the instance loaded in memory, we build the JuMP optimization model using UnitCommitment.build_model
:
model =
UnitCommitment.build_model(instance = instance, optimizer = HiGHS.Optimizer);
[ Info: Building model...
[ Info: Building scenario s1 with probability 1.0
[ Info: Computing injection shift factors...
[ Info: Computed ISF in 0.00 seconds
[ Info: Computing line outage factors...
[ Info: Computed LODF in 0.00 seconds
[ Info: Applying PTDF and LODF cutoffs (0.00500, 0.00100)
[ Info: Built model in 0.01 seconds
Next, we run the optimization process, with UnitCommitment.optimize!
:
UnitCommitment.optimize!(model)
[ Info: Setting MILP time limit to 86400.00 seconds
[ Info: Solving MILP...
Running HiGHS 1.6.0: Copyright (c) 2023 HiGHS under MIT licence terms
Presolving model
4382 rows, 2704 cols, 14776 nonzeros
3279 rows, 2195 cols, 11913 nonzeros
3160 rows, 2085 cols, 12983 nonzeros
Solving MIP model with:
3160 rows
2085 cols (865 binary, 0 integer, 0 implied int., 1220 continuous)
12983 nonzeros
Nodes | B&B Tree | Objective Bounds | Dynamic Constraints | Work
Proc. InQueue | Leaves Expl. | BestBound BestSol Gap | Cuts InLp Confl. | LpIters Time
0 0 0 0.00% 1.86264515e-09 inf inf 0 0 0 0 0.0s
R 0 0 0 0.00% 360642.328869 360642.544974 0.00% 0 0 0 861 0.1s
Solving report
Status Optimal
Primal bound 360642.544974
Dual bound 360642.328869
Gap 6e-05% (tolerance: 0.01%)
Solution status feasible
360642.544974 (objective)
0 (bound viol.)
0 (int. viol.)
0 (row viol.)
Timing 0.05 (total)
0.03 (presolve)
0.00 (postsolve)
Nodes 1
LP iterations 861 (total)
0 (strong br.)
0 (separation)
0 (heuristics)
[ Info: Verifying transmission limits...
[ Info: Verified transmission limits in 0.00 seconds
[ Info: No violations found
Finally, we extract the optimal solution from the model:
solution = UnitCommitment.solution(model)
OrderedCollections.OrderedDict{Any, Any} with 15 entries:
"Thermal production (MW)" => OrderedDict("g1"=>[181.92, 172.85, 166.8…
"Thermal production cost (\$)" => OrderedDict("g1"=>[7241.95, 6839.01, 657…
"Startup cost (\$)" => OrderedDict("g1"=>[0.0, 0.0, 0.0, 0.0, 0…
"Is on" => OrderedDict("g1"=>[1.0, 1.0, 1.0, 1.0, 1…
"Switch on" => OrderedDict("g1"=>[0.0, 0.0, 0.0, 0.0, 0…
"Switch off" => OrderedDict("g1"=>[0.0, 0.0, 0.0, 0.0, 0…
"Net injection (MW)" => OrderedDict("b1"=>[181.92, 172.85, 166.8…
"Load curtail (MW)" => OrderedDict("b1"=>[0.0, 0.0, 0.0, 0.0, 0…
"Line overflow (MW)" => OrderedDict("l1"=>[0.0, 0.0, 0.0, 0.0, 0…
"Spinning reserve (MW)" => OrderedDict("r1"=>OrderedDict("g1"=>[4.6…
"Spinning reserve shortfall (MW)" => OrderedDict("r1"=>[0.0, 0.0, 0.0, 0.0, 0…
"Up-flexiramp (MW)" => OrderedDict{Any, Any}()
"Up-flexiramp shortfall (MW)" => OrderedDict{Any, Any}()
"Down-flexiramp (MW)" => OrderedDict{Any, Any}()
"Down-flexiramp shortfall (MW)" => OrderedDict{Any, Any}()
We can then explore the solution using Julia:
@show solution["Thermal production (MW)"]["g1"]
36-element Vector{Float64}:
181.92024301829747
172.8503730182975
166.8067830182975
163.2384530182975
165.5149530182975
169.95052301829747
182.3719130182975
191.43277301829744
196.67234301829745
202.53524301829742
⋮
165.0181030182975
167.24040301829746
180.89039301829752
196.95238301829744
206.5552530182975
214.18877301829747
225.48092301829746
226.81792301829745
222.7256530182974
Or export the entire solution to a JSON file:
UnitCommitment.write("solution.json", solution)
Solving a custom deterministic instance
In the previous example, we solved a benchmark instance provided by the package. To solve a custom instance, the first step is to create an input file describing the list of elements (generators, loads and transmission lines) in the network. See Data Format for a complete description of the data format UC.jl expects. To keep this tutorial self-contained, we will create the input JSON file using Julia; however, this step can also be done with a simple text editor. First, we define the contents of the file:
json_contents = """
{
"Parameters": {
"Version": "0.4",
"Time horizon (h)": 4
},
"Buses": {
"b1": {
"Load (MW)": [100, 150, 200, 250]
}
},
"Generators": {
"g1": {
"Bus": "b1",
"Type": "Thermal",
"Production cost curve (MW)": [0, 200],
"Production cost curve (\$)": [0, 1000],
"Initial status (h)": -24,
"Initial power (MW)": 0
},
"g2": {
"Bus": "b1",
"Type": "Thermal",
"Production cost curve (MW)": [0, 300],
"Production cost curve (\$)": [0, 3000],
"Initial status (h)": -24,
"Initial power (MW)": 0
}
}
}
""";
Next, we write it to example.json
.
open("example.json", "w") do file
return write(file, json_contents)
end;
Now that we have the input file, we can proceed as before, but using UnitCommitment.read
instead of UnitCommitment.read_benchmark
:
instance = UnitCommitment.read("example.json");
model =
UnitCommitment.build_model(instance = instance, optimizer = HiGHS.Optimizer);
UnitCommitment.optimize!(model)
[ Info: Building model...
[ Info: Building scenario s1 with probability 1.0
[ Info: Built model in 0.00 seconds
[ Info: Setting MILP time limit to 86400.00 seconds
[ Info: Solving MILP...
Running HiGHS 1.6.0: Copyright (c) 2023 HiGHS under MIT licence terms
Presolving model
74 rows, 36 cols, 166 nonzeros
63 rows, 29 cols, 164 nonzeros
49 rows, 25 cols, 120 nonzeros
39 rows, 19 cols, 97 nonzeros
32 rows, 15 cols, 90 nonzeros
24 rows, 13 cols, 64 nonzeros
8 rows, 6 cols, 21 nonzeros
0 rows, 0 cols, 0 nonzeros
Presolve: Optimal
Solving report
Status Optimal
Primal bound 3750
Dual bound 3750
Gap 0% (tolerance: 0.01%)
Solution status feasible
3750 (objective)
0 (bound viol.)
0 (int. viol.)
0 (row viol.)
Timing 0.00 (total)
0.00 (presolve)
0.00 (postsolve)
Nodes 0
LP iterations 0 (total)
0 (strong br.)
0 (separation)
0 (heuristics)
Finally, we extract and display the solution:
solution = UnitCommitment.solution(model)
OrderedCollections.OrderedDict{Any, Any} with 14 entries:
"Thermal production (MW)" => OrderedDict("g1"=>[100.0, 150.0, 200.0, …
"Thermal production cost (\$)" => OrderedDict("g1"=>[500.0, 750.0, 1000.0,…
"Startup cost (\$)" => OrderedDict("g1"=>[0.0, 0.0, 0.0, 0.0], …
"Is on" => OrderedDict("g1"=>[1.0, 1.0, 1.0, 1.0], …
"Switch on" => OrderedDict("g1"=>[1.0, 0.0, 0.0, 0.0], …
"Switch off" => OrderedDict("g1"=>[0.0, 0.0, 0.0, 0.0], …
"Net injection (MW)" => OrderedDict("b1"=>[0.0, 0.0, 0.0, 0.0])
"Load curtail (MW)" => OrderedDict("b1"=>[0.0, 0.0, 0.0, 0.0])
"Spinning reserve (MW)" => OrderedDict{Any, Any}()
"Spinning reserve shortfall (MW)" => OrderedDict{Any, Any}()
"Up-flexiramp (MW)" => OrderedDict{Any, Any}()
"Up-flexiramp shortfall (MW)" => OrderedDict{Any, Any}()
"Down-flexiramp (MW)" => OrderedDict{Any, Any}()
"Down-flexiramp shortfall (MW)" => OrderedDict{Any, Any}()
@show solution["Thermal production (MW)"]["g1"]
4-element Vector{Float64}:
100.0
150.0
200.0
200.0
@show solution["Thermal production (MW)"]["g2"]
4-element Vector{Float64}:
0.0
0.0
0.0
50.0
Solving a custom stochastic instance
In addition to deterministic test cases, UnitCommitment.jl can also solve two-stage stochastic instances of the problem. In this section, we demonstrate the most simple form, which builds a single (extensive form) model containing information for all scenarios. See Decomposition for more advanced methods.
First, we need to create one JSON input file for each scenario. Parameters that are allowed to change across scenarios are marked as "uncertain" in the JSON data format page. It is also possible to specify the name and weight of each scenario, as shown below.
We start by creating example_s1.json
, the first scenario file:
json_contents_s1 = """
{
"Parameters": {
"Version": "0.4",
"Time horizon (h)": 4,
"Scenario name": "s1",
"Scenario weight": 3.0
},
"Buses": {
"b1": {
"Load (MW)": [100, 150, 200, 250]
}
},
"Generators": {
"g1": {
"Bus": "b1",
"Type": "Thermal",
"Production cost curve (MW)": [0, 200],
"Production cost curve (\$)": [0, 1000],
"Initial status (h)": -24,
"Initial power (MW)": 0
},
"g2": {
"Bus": "b1",
"Type": "Thermal",
"Production cost curve (MW)": [0, 300],
"Production cost curve (\$)": [0, 3000],
"Initial status (h)": -24,
"Initial power (MW)": 0
}
}
}
"""
open("example_s1.json", "w") do file
return write(file, json_contents_s1)
end;
Next, we create example_s2.json
, the second scenario file:
json_contents_s2 = """
{
"Parameters": {
"Version": "0.4",
"Time horizon (h)": 4,
"Scenario name": "s2",
"Scenario weight": 1.0
},
"Buses": {
"b1": {
"Load (MW)": [200, 300, 400, 500]
}
},
"Generators": {
"g1": {
"Bus": "b1",
"Type": "Thermal",
"Production cost curve (MW)": [0, 200],
"Production cost curve (\$)": [0, 1000],
"Initial status (h)": -24,
"Initial power (MW)": 0
},
"g2": {
"Bus": "b1",
"Type": "Thermal",
"Production cost curve (MW)": [0, 300],
"Production cost curve (\$)": [0, 3000],
"Initial status (h)": -24,
"Initial power (MW)": 0
}
}
}
""";
open("example_s2.json", "w") do file
return write(file, json_contents_s2)
end;
Now that we have our two scenario files, we can read them using UnitCommitment.read
. Note that, instead of a single file, we now provide a list.
instance = UnitCommitment.read(["example_s1.json", "example_s2.json"])
UnitCommitmentInstance(2 scenarios, 2 thermal units, 0 profiled units, 1 buses, 0 lines, 0 contingencies, 0 price sensitive loads, 4 time steps)
If we have a large number of scenario files, the Glob package can also be used to avoid having to list them individually:
using Glob
instance = UnitCommitment.read(glob("example_s*.json"))
UnitCommitmentInstance(2 scenarios, 2 thermal units, 0 profiled units, 1 buses, 0 lines, 0 contingencies, 0 price sensitive loads, 4 time steps)
Finally, we build the model and optimize as before:
model =
UnitCommitment.build_model(instance = instance, optimizer = HiGHS.Optimizer);
UnitCommitment.optimize!(model)
[ Info: Building model...
[ Info: Building scenario s1 with probability 0.75
[ Info: Building scenario s2 with probability 0.25
[ Info: Built model in 0.00 seconds
[ Info: Setting MILP time limit to 86400.00 seconds
[ Info: Solving MILP...
Running HiGHS 1.6.0: Copyright (c) 2023 HiGHS under MIT licence terms
Presolving model
115 rows, 47 cols, 251 nonzeros
47 rows, 43 cols, 113 nonzeros
33 rows, 33 cols, 93 nonzeros
23 rows, 25 cols, 61 nonzeros
7 rows, 9 cols, 21 nonzeros
5 rows, 7 cols, 13 nonzeros
Solving MIP model with:
5 rows
7 cols (4 binary, 0 integer, 0 implied int., 3 continuous)
13 nonzeros
Nodes | B&B Tree | Objective Bounds | Dynamic Constraints | Work
Proc. InQueue | Leaves Expl. | BestBound BestSol Gap | Cuts InLp Confl. | LpIters Time
0 0 0 0.00% 4562.5 inf inf 0 0 0 0 0.0s
T 0 0 0 0.00% 4562.5 5312.5 14.12% 0 0 0 2 0.0s
Solving report
Status Optimal
Primal bound 5312.5
Dual bound 5312.5
Gap 0% (tolerance: 0.01%)
Solution status feasible
5312.5 (objective)
0 (bound viol.)
0 (int. viol.)
0 (row viol.)
Timing 0.00 (total)
0.00 (presolve)
0.00 (postsolve)
Nodes 1
LP iterations 2 (total)
0 (strong br.)
0 (separation)
0 (heuristics)
The solution to stochastic instances follows a slightly different format, as shown below:
solution = UnitCommitment.solution(model)
OrderedCollections.OrderedDict{Any, Any} with 2 entries:
"s1" => OrderedDict{Any, Any}("Thermal production (MW)"=>OrderedDict("g1"=>[1…
"s2" => OrderedDict{Any, Any}("Thermal production (MW)"=>OrderedDict("g1"=>[2…
The solution for each scenario can be accessed through solution[scenario_name]
. For conveniance, this includes both first- and second-stage optimal decisions:
solution["s1"]
OrderedCollections.OrderedDict{Any, Any} with 14 entries:
"Thermal production (MW)" => OrderedDict("g1"=>[100.0, 150.0, 200.0, …
"Thermal production cost (\$)" => OrderedDict("g1"=>[500.0, 750.0, 1000.0,…
"Startup cost (\$)" => OrderedDict("g1"=>[0.0, 0.0, 0.0, 0.0], …
"Is on" => OrderedDict("g1"=>[1.0, 1.0, 1.0, 1.0], …
"Switch on" => OrderedDict("g1"=>[1.0, -0.0, 0.0, 0.0],…
"Switch off" => OrderedDict("g1"=>[0.0, 0.0, 0.0, 0.0], …
"Net injection (MW)" => OrderedDict("b1"=>[0.0, 0.0, 0.0, 0.0])
"Load curtail (MW)" => OrderedDict("b1"=>[0.0, 0.0, 0.0, 0.0])
"Spinning reserve (MW)" => OrderedDict{Any, Any}()
"Spinning reserve shortfall (MW)" => OrderedDict{Any, Any}()
"Up-flexiramp (MW)" => OrderedDict{Any, Any}()
"Up-flexiramp shortfall (MW)" => OrderedDict{Any, Any}()
"Down-flexiramp (MW)" => OrderedDict{Any, Any}()
"Down-flexiramp shortfall (MW)" => OrderedDict{Any, Any}()