Economic dispatch¶
Economic dispatch is a classic optimization problem in power systems. It is used to determine the optimal power output of a set of generators to meet the demand at the lowest cost. The problem can be formulated as a linear programming problem or a quadratic programming problem depending on the formulation used to model the economic cost.
In this example, we will show how to use PyOptInterface to solve an economic dispatch problem using the quadratic programming formulation.
Problem Formulation¶
We consider a system with \(N_G\) generators and \(N_D\) demands. The decision variables are the power output of the generators \(P_{i,t}\) for \(i=1,\ldots,N_G\) and \(t=1,\ldots,T\). The objective function is the total cost of the generators, which is the sum of the quadratic cost, the linear cost, and the constant cost. The first constraint ensures that the total power output of the generators meets the demand. The second and third constraints are the power output limits of the generators. The last constraint is the ramping limits of the generators.
Implementation¶
Firstly, we need to create a model object. We will use the HiGHS solver in this example because it is an excellent open-source solver and supports quadratic programming problems. You need to read the installation guide to install the HiGHS solver manually and set the path to shared library of HiGHS as described in the guide.
import pyoptinterface as poi
from pyoptinterface import highs
model = highs.Model()
Then, we need to create new variables in the model to represent the power output of the generators. We will use the add_variables method to add multidimensional variables to the model.
T = 24
N_G = 3
P_min = [50, 50, 50]
P_max = [100, 100, 100]
a_cost = [0.01, 0.015, 0.02]
b_cost = [1, 1.5, 2]
c_cost = [0, 0, 0]
Rup = [20, 20, 20]
Rdn = [20, 20, 20]
P = model.add_variables(range(N_G), range(T), name="P")
Running HiGHS 1.8.1 (git hash: 4a7f24a): Copyright (c) 2024 HiGHS under MIT licence terms
For the lower and upper bound of variables, we can set the corresponding variable attribute to set them.
for i in range(N_G):
for t in range(T):
model.set_variable_attribute(P[i, t], poi.VariableAttribute.LowerBound, P_min[i])
model.set_variable_attribute(P[i, t], poi.VariableAttribute.UpperBound, P_max[i])
Next, we need to add the power balance constraint to the model. We will use the add_linear_constraint method to add a linear constraint to the model. The multidimensional constraint is managed by tupledict and we use the make_tupledict method to create a multidimensional constraint.
The total demand is assumed to be a constant in this example.
total_demand = [220 for _ in range(T)]
def con(t):
lhs = poi.quicksum(P[i, t] for i in range(N_G))
rhs = total_demand[t]
return model.add_linear_constraint(lhs, poi.Eq, rhs, name=f"powerbalance_{t}")
powebalance_constraints = poi.make_tupledict(range(T), rule=con)
def rampup_con(i, t):
if t == 0:
return None
lhs = P[i, t] - P[i, t-1]
rhs = Rup[i]
return model.add_linear_constraint(lhs, poi.Leq, rhs, name=f"rampup_{i}_{t}")
rampup_constraints = poi.make_tupledict(range(N_G), range(T), rule=rampup_con)
def rampdown_con(i, t):
if t == 0:
return None
lhs = P[i, t] - P[i, t-1]
rhs = -Rdn[i]
return model.add_linear_constraint(lhs, poi.Geq, rhs, name=f"rampdown_{i}_{t}")
rampdown_constraints = poi.make_tupledict(range(N_G), range(T), rule=rampdown_con)
Then, we need to add the quadratic objective function to the model. We will use the set_objective method to set the objective function of the model.
obj = poi.ExprBuilder()
for t in range(T):
for i in range(N_G):
obj += a_cost[i] * P[i, t] * P[i, t] + b_cost[i] * P[i, t] + c_cost[i]
model.set_objective(obj)
Finally, we can solve the model and query the solution.
model.optimize()
print(model.get_model_attribute(poi.ModelAttribute.TerminationStatus))
print("Objective value: ", model.get_value(obj))
Coefficient ranges:
Matrix [1e+00, 1e+00]
Cost [1e+00, 2e+00]
Bound [5e+01, 1e+02]
RHS [2e+01, 2e+02]
Iteration Objective NullspaceDim
0 14328.52 0 0.00s
100 12722.021 1 0.00s
108 12684.021 0 0.00s
Model status : Optimal
Simplex iterations: 63
QP ASM iterations: 108
Objective value : 1.2684000000e+04
HiGHS run time : 0.00
TerminationStatusCode.OPTIMAL
Objective value: 12684.0
The optimal value of decision variables can be queried via get_value function.
import numpy as np
P_value = np.fromiter(
(model.get_value(P[i, t]) for i in range(N_G) for t in range(T)), float
).reshape(N_G, T)
P_value
array([[100., 100., 100., 100., 100., 100., 100., 100., 100., 100., 100.,
100., 100., 100., 100., 100., 100., 100., 100., 100., 100., 100.,
100., 100.],
[ 70., 70., 70., 70., 70., 70., 70., 70., 70., 70., 70.,
70., 70., 70., 70., 70., 70., 70., 70., 70., 70., 70.,
70., 70.],
[ 50., 50., 50., 50., 50., 50., 50., 50., 50., 50., 50.,
50., 50., 50., 50., 50., 50., 50., 50., 50., 50., 50.,
50., 50.]])
Change the load and solve the model again¶
We can change the load and solve the model again without creating a new model from scratch by modifying the right-hand side of the power balance constraint.
For example, we increase the load and solve the model again.
total_demand = [220 + 1.0 * t for t in range(T)]
for t in range(T):
model.set_normalized_rhs(powebalance_constraints[t], total_demand[t])
model.optimize()
print(model.get_model_attribute(poi.ModelAttribute.TerminationStatus))
print("Objective value: ", model.get_value(obj))
P_value = np.fromiter(
(model.get_value(P[i, t]) for i in range(N_G) for t in range(T)), float
).reshape(N_G, T)
print(P_value)
Coefficient ranges:
Matrix [1e+00, 1e+00]
Cost [1e+00, 2e+00]
Bound [5e+01, 1e+02]
RHS [2e+01, 2e+02]
Iteration Objective NullspaceDim
0 15430.248 0 0.00s
100 14091.031 3 0.01s
169 13740.237 10 0.01s
Model status : Optimal
Simplex iterations: 86
QP ASM iterations: 169
Objective value : 1.3740213571e+04
HiGHS run time : 0.01
TerminationStatusCode.OPTIMAL
Objective value: 13740.213571429404
[[100. 100. 100. 100. 100.
100. 100. 100. 100. 100.
100. 100. 100. 100. 100.
100. 100. 100. 100. 100.
100. 100. 100. 100. ]
[ 70. 71. 72. 73. 74.
75. 76. 77. 78. 79.
80. 81. 82. 83. 83.71423796
84.28566633 84.85709469 85.42852306 85.99995143 86.5713798
87.14280816 87.71423653 88.2856649 88.85709327]
[ 50. 50. 50. 50. 50.
50. 50. 50. 50. 50.
50. 50. 50. 50. 50.28576204
50.71433367 51.14290531 51.57147694 52.00004857 52.4286202
52.85719184 53.28576347 53.7143351 54.14290673]]