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}()