In the previous tutorial, we used UnitCommitment.jl to solve benchmark and user-provided instances using a default mathematical formulation for the problem. In this tutorial, we will explore how to customize this formulation.
This tutorial is not required for using UnitCommitment.jl, unless you plan to make changes to the problem formulation. In this page, we assume familiarity with the JuMP modeling language. Please see JuMP's official documentation for resources on getting started with JuMP.
By default, UnitCommitment.build_model uses a formulation that combines modeling components from different publications, and that has been carefully tested, using our own benchmark scripts, to provide good performance across a wide variety of instances. This default formulation is expected to change over time, as new methods are proposed in the literature. You can, however, construct your own formulation, based on the modeling components that you choose, as shown in the next example.
We start by importing the necessary packages and reading a benchmark instance:
using HiGHS
using JuMP
using UnitCommitment
@@ -22,11 +22,9 @@ instance = UnitCommitment.read_benchmark("matpower/case14/2017-01-01")
[ Info: Computing line outage factors...
[ Info: Computed LODF in 0.00 seconds
[ Info: Applying PTDF and LODF cutoffs (0.00800, 0.00300)
-[ Info: Built model in 0.94 seconds
In the previous tutorial, we saw how to access the optimal solution through UnitCommitment.solution. While this approach works well for basic usage, it is also possible to get a direct reference to the JuMP decision variables and query their values, as the next example illustrates.
First, we load a benchmark instance and solve it, as before.
instance = UnitCommitment.read_benchmark("matpower/case14/2017-01-01");
-model = UnitCommitment.build_model(
- instance=instance,
- optimizer=HiGHS.Optimizer,
-);
+[ Info: Built model in 0.87 seconds
In the previous tutorial, we saw how to access the optimal solution through UnitCommitment.solution. While this approach works well for basic usage, it is also possible to get a direct reference to the JuMP decision variables and query their values, as the next example illustrates.
First, we load a benchmark instance and solve it, as before.
instance = UnitCommitment.read_benchmark("matpower/case14/2017-01-01");
+model =
+ UnitCommitment.build_model(instance = instance, optimizer = HiGHS.Optimizer);
UnitCommitment.optimize!(model)
[ Info: Building model...
[ Info: Building scenario s1 with probability 1.0
[ Info: Computing injection shift factors...
@@ -51,7 +49,7 @@ Solving MIP model with:
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.1s
+ 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
@@ -64,8 +62,8 @@ Solving report
0 (bound viol.)
0 (int. viol.)
0 (row viol.)
- Timing 0.10 (total)
- 0.06 (presolve)
+ Timing 0.07 (total)
+ 0.04 (presolve)
0.00 (postsolve)
Nodes 1
LP iterations 861 (total)
@@ -73,26 +71,17 @@ Solving report
0 (separation)
0 (heuristics)
[ Info: Verifying transmission limits...
-[ Info: Verified transmission limits in 0.57 seconds
-[ Info: No violations found
At this point, it is possible to obtain a reference to the decision variables by calling model[:varname][index]. For example, model[:is_on]["g1",1] returns a direct reference to the JuMP variable indicating whether generator named "g1" is on at time 1. For a complete list of decision variables available, and how are they indexed, see the problem definition.
@show JuMP.value(model[:is_on]["g1",1])
1.0
To access second-stage decisions, it is necessary to specify the scenario name. UnitCommitment.jl models deterministic instances as a particular case in which there is a single scenario named "s1", so we need to use this key.
@show JuMP.value(model[:prod_above]["s1", "g1", 1])
145.16900936037393
When testing variations of the unit commitment problem, it is often necessary to modify the objective function, variables and constraints of the formulation. UnitCommitment.jl makes this process relatively easy. The first step is to construct the standard model using UnitCommitment.build_model:
instance = UnitCommitment.read_benchmark("matpower/case14/2017-01-01");
-model = UnitCommitment.build_model(
- instance=instance,
- optimizer=HiGHS.Optimizer,
-);
[ Info: Building model...
+[ Info: Verified transmission limits in 0.54 seconds
+[ Info: No violations found
At this point, it is possible to obtain a reference to the decision variables by calling model[:varname][index]. For example, model[:is_on]["g1",1] returns a direct reference to the JuMP variable indicating whether generator named "g1" is on at time 1. For a complete list of decision variables available, and how are they indexed, see the problem definition.
@show JuMP.value(model[:is_on]["g1", 1])
1.0
To access second-stage decisions, it is necessary to specify the scenario name. UnitCommitment.jl models deterministic instances as a particular case in which there is a single scenario named "s1", so we need to use this key.
@show JuMP.value(model[:prod_above]["s1", "g1", 1])
145.16900936037393
When testing variations of the unit commitment problem, it is often necessary to modify the objective function, variables and constraints of the formulation. UnitCommitment.jl makes this process relatively easy. The first step is to construct the standard model using UnitCommitment.build_model:
instance = UnitCommitment.read_benchmark("matpower/case14/2017-01-01");
+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
Now, before calling UnitCommitment.optimize, we can make any desired changes to the formulation. In the previous section, we saw how to obtain a direct reference to the decision variables. It is possible to modify them by using standard JuMP methods. For example, to fix the commitment status of a particular generator, we can use JuMP.fix:
JuMP.fix(model[:is_on]["g1",1], 1.0, force=true)
To modify the cost coefficient of a particular variable, we can use JuMP.set_objective_coefficient:
JuMP.set_objective_coefficient(
- model,
- model[:switch_on]["g1",1],
- 1000.0,
-)
It is also possible to make changes to the set of constraints. For example, we can add a custom constraint, using the JuMP.@constraint macro:
@constraint(
- model,
- model[:is_on]["g3",1] + model[:is_on]["g4",1] <= 1,
-);
We can also remove an existing model constraint using JuMP.delete. See the problem definition for a list of constraint names and indices.
JuMP.delete(model, model[:eq_min_uptime]["g1",1])
After we are done with all changes, we can call UnitCommitment.optimize and extract the optimal solution:
UnitCommitment.optimize!(model)
+[ Info: Built model in 0.01 seconds
Now, before calling UnitCommitment.optimize, we can make any desired changes to the formulation. In the previous section, we saw how to obtain a direct reference to the decision variables. It is possible to modify them by using standard JuMP methods. For example, to fix the commitment status of a particular generator, we can use JuMP.fix:
JuMP.fix(model[:is_on]["g1", 1], 1.0, force = true)
To modify the cost coefficient of a particular variable, we can use JuMP.set_objective_coefficient:
JuMP.set_objective_coefficient(model, model[:switch_on]["g1", 1], 1000.0)
It is also possible to make changes to the set of constraints. For example, we can add a custom constraint, using the JuMP.@constraint macro:
@constraint(model, model[:is_on]["g3", 1] + model[:is_on]["g4", 1] <= 1,);
We can also remove an existing model constraint using JuMP.delete. See the problem definition for a list of constraint names and indices.
JuMP.delete(model, model[:eq_min_uptime]["g1", 1])
After we are done with all changes, we can call UnitCommitment.optimize and extract the optimal solution:
UnitCommitment.optimize!(model)
@show 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…
@@ -109,18 +98,16 @@ model = UnitCommitment.build_model(
"Up-flexiramp shortfall (MW)" => OrderedDict{Any, Any}()
"Down-flexiramp (MW)" => OrderedDict{Any, Any}()
"Down-flexiramp shortfall (MW)" => OrderedDict{Any, Any}()
In this section we demonstrate how to add a new grid component to a particular bus in the network. This is useful, for example, when developing formulations for a new type of generator, energy storage, or any other grid device. We start by reading the instance data and buliding a standard model:
instance = UnitCommitment.read_benchmark("matpower/case118/2017-02-01")
-model = UnitCommitment.build_model(
- instance=instance,
- optimizer=HiGHS.Optimizer,
-);
[ Info: Building 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.25 seconds
+[ Info: Computed ISF in 0.22 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.46 seconds
Next, we create decision variables for the new grid component. In this example, we assume that the new component can inject up to 10 MW of power at each time step, so we create new continuous variables $0 \leq x_t \leq 10$.
T = instance.time
-@variable(model, x[1:T], lower_bound=0.0, upper_bound=10.0);
Next, we add the production costs to the objective function. In this example, we assume a generation cost of $5/MW:
for t in 1:T
+[ Info: Built model in 0.42 seconds
Next, we create decision variables for the new grid component. In this example, we assume that the new component can inject up to 10 MW of power at each time step, so we create new continuous variables $0 \leq x_t \leq 10$.
T = instance.time
+@variable(model, x[1:T], lower_bound = 0.0, upper_bound = 10.0);
Next, we add the production costs to the objective function. In this example, we assume a generation cost of $5/MW:
for t in 1:T
set_objective_coefficient(model, x[t], 5.0)
end
We then attach the new component to bus b1 by modifying the net injection constraint (eq_net_injection):
for t in 1:T
set_normalized_coefficient(
@@ -145,7 +132,7 @@ Solving MIP model with:
Proc. InQueue | Leaves Expl. | BestBound BestSol Gap | Cuts InLp Confl. | LpIters Time
0 0 0 0.00% -501782.068184 inf inf 0 0 0 0 0.5s
- R 0 0 0 0.00% 5748240.007416 5748250.232723 0.00% 0 0 0 8320 1.0s
+ R 0 0 0 0.00% 5748240.007416 5748250.232723 0.00% 0 0 0 8320 0.9s
Solving report
Status Optimal
@@ -157,8 +144,8 @@ Solving report
0 (bound viol.)
2.22044604925e-15 (int. viol.)
0 (row viol.)
- Timing 1.04 (total)
- 0.39 (presolve)
+ Timing 0.91 (total)
+ 0.37 (presolve)
0.00 (postsolve)
Nodes 1
LP iterations 8320 (total)
@@ -187,4 +174,4 @@ Solving report
10.0
10.0
10.0
- 10.0
![](\"../assets/cost_curve.png\")
\n