Extra material: Energy dispatch problem#

To meet the energy demand, power plants run day and night across the country to produce electricity from a variety of sources such as fossil fuels and renewable energy. On the short-time scale, the best operating levels for electric power plants are derived every 15 minutes by solving the so-called Optimal Power Flow (OPF) model. The OPF model is an optimization problem with the objective of minimizing the total energy dispatching cost, while ensuring that the generation meets the total energy demand. Furthermore, the model takes into account many constraints, among which operational and physical constraints.

# install dependencies and select solver
%pip install -q amplpy matplotlib networkx numpy pandas

SOLVER = "highs"

from amplpy import AMPL, ampl_notebook

ampl = ampl_notebook(
    modules=["highs"],  # modules to install
    license_uuid="default",  # license to use
)  # instantiate AMPL object and register magics
Using default Community Edition License for Colab. Get yours at: https://ampl.com/ce
Licensed to AMPL Community Edition License for the AMPL Model Colaboratory (https://colab.ampl.com).
import time
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import networkx as nx

Background: Power networks and power flow physics#

We model the nation-wide transmission power network as a directed graph \(G=(V, E)\), where \(V\) represents the set of nodes (e.g., cities, industrial districts, power plants) and \(E\) denotes the set of directed edges (e.g., physical transmission lines).

Each node \(i \in V\) has a power injection \(p_i\) and demand \(d_i\). The set of nodes are separated into generator and load nodes. The set of generators \(\mathcal{G} \subseteq V\) corresponds to the nodes \(i \in V\) for which \(p_i \geq 0\) and \(d_i = 0\). Each generator \(i \in \mathcal{G}\) has a minimum \(p_i^{\min}\) and maximum \(p_i^{\max}\) production capacity. The set of load nodes \(\mathcal{D} \subseteq V\) corresponds to the nodes for which \(p_i = 0\) and \(d_i \geq 0\). The load nodes thus correspond to the places where electricity is being consumed, e.g., cities and industrial districts. We say that supply and demand is matched if \(\sum_{i \in V} p_i - d_i = 0\). Since we cannot store electricity in large volumes, supply must meet demand at all times, hence adjusting to it every 15 minutes by solving the OPF.

Each edge \((i, j) \in E\) carries a power flow \(f_{ij} \in R\) and has a capacity \(f_{ij}^{\max} \geq 0\), i.e., the maximum power flow that it may carry. Note that our choice to model a directed graph is to make the modeling of the network easier. In particular, a directed edge \((i,j)\) may carry a ā€˜negativeā€™ flow \(f_{ij} < 0\), which implies that there is flow going from \(j\) to \(i\) where \(f_{ji} = -f_{ij}\). The capacity does not depend on the direction of the flow, implying that the flow capacity constraints are given by \(|f_{ij}| = |f_{ji}| \leq f_{ij}^{\max}\).

One crucial difference of power flow problems compared to typical network flow problems is that the power flows cannot be controlled directly. Instead, as you might recall from high-school physics, the power flows are determined by the laws of electricity, which we will now present as the power flow equations. Ignore for a moment the flow capacity constraints. Let \(\theta_{i} \in \mathbb{R}\) denote the phase angle of node \(i\). For each edge \((i,j)\), let \(b_{ij} > 0\) denote the line susceptance. Assuming that supply and demand is matched, i.e., \(\sum_{i=1}^{n} p_i - d_i = 0\), the power flows \(\mathbf{f} \in \mathbb{R}^{m}\) and phase angles \(\mathbf{\theta} \in \mathbb{R}^{n}\) are obtained by solving the following linear system of equations:

\[\begin{split}\begin{align} p_i - d_i &= \sum_{j: (i, j) \in E} f_{ij} - \sum_{j: (j, i) \in E} f_{ji}, & \forall \, i \in V,\\ f_{ij} &= b_{ij}(\theta_i - \theta_j), & \forall \, (i, j) \in E. \end{align}\end{split}\]

The first set of constraints ensures flow conservation and the second set of constrations captures the flow dependency on susceptances and angle differences. The DC power flow equations admit a unique power flow solution \(\mathbf{f}\) given matched power injections \(\mathbf{p}\) and demand \(\mathbf{d}\).

For a given matched supply and demand vector \(\mathbf{p}\) and \(\mathbf{d}\), we can compute the power flows on the network by solving the linear equations as described above. There are exactly \(|V|\) and \(|E|\) equations for the \(\theta_i\) variables and \(f_{ij}\) variables, meaning that this system of equations admit a solution.

Optimal Power Flow#

We assumed above that the power injections \(\mathbf{p}\) were given. However, in practice, the power injections need to be determined for each generator in the power network, where some types of generators may be cheaper than others. Moreover, we need to take into account operational constraints, such as the maximum flow and generator limits.

On the short-time scale, the power injections are calculated for each generator by solving the so-called Optimal Power Flow (OPF) problem. The goal of the OPF problem is to determine a solution \((\mathbf{p}, \mathbf{f}, \mathbf{\theta})\) with minimal costs such that:

  • Supply meets demand

  • Line capacity constraints are met

  • Generator capacity constraints are met

Let \(c_i > 0\) be the cost associated with the production of a unit energy by generator \(i\). Then, the OPF problem can be formulated as

\[\begin{split}\begin{align} \begin{array}{llll} \min & \sum_{i \in V} c_i p_i \\ \text{s.t.} & \sum_{j: (i, j) \in E} f_{ij} - \sum_{j: (j, i) \in E} f_{ji} = p_i - d_i & \forall \, i \in V,\\ & f_{ij} = b_{ij}(\theta_i - \theta_j), & \forall \, (i, j) \in E, \\ & |f_{ij}| \leq f_{ij}^{\max} & \forall (i, j) \in E,\\ & p_{i}^{\min } \leq p_{i} \leq p_{i}^{\max } & \forall i \in V, \\ & p_i \in \mathbb{R}_{\geq 0} & \forall i \in V, \\ & \theta_i \in \mathbb{R} & \forall i \in V, \\ & f_{ij} \in \mathbb{R} & \forall (i, j) \in E \\ \end{array} \end{align}\end{split}\]

For simplicity, you may assume that all load nodes do not produce energy, i.e., \(p_i = p_i^{\min} = p_i^{\max} = 0\) for all \(i \in \mathcal{D}\). You may therefore model \(p_i\) as decision variables for all nodes (both generator and load nodes). Similarly, you may assume that all generator nodes have no demand, i.e., \(d_i = 0\) for all \(i \in \mathcal{G}\).

To summarize, the decision variables in the OPF problem are:

  • \(p_i\) power injections

  • \(\theta_i\) phase angles

  • \(f_{ij}\) power flows

All the other quantities are instance dependent parameters.

Data#

You will solve the OPF problem on the IEEE-118 power network, which is a standard test network consisting of 118 nodes and 179 edges. In the following, we will load the data into the notebook and provide a description of the data.

Data import#

# Download the data
base_url = (
    "https://raw.githubusercontent.com/ampl/mo-book.ampl.com/dev/notebooks/04/"
)
nodes_df = pd.read_csv(base_url + "nodes.csv").set_index(["node_id", "instance"])
edges_df = pd.read_csv(base_url + "edges.csv").set_index(["node_id1", "node_id2"])

# Replace 'na' by an empty string
nodes_df.fillna("", inplace=True)

# Split the initial nodes data frame by instance
nodes_by_instance = [
    y.reset_index(level=1).drop("instance", axis=1)
    for x, y in nodes_df.groupby("instance")
]
I = [{"nodes": n, "edges": edges_df} for n in nodes_by_instance]

# Initialize a network for demonstration purposes
network = I[0]

Network data#

def visualize_network(network, edge_flows=None, ax=None):
    """Visualize a network instance, highlighting the generators in orange and the load buses in green."""
    plt.figure(figsize=[12, 10])
    g = nx.DiGraph(list(network["edges"].index))
    pos = nx.layout.kamada_kawai_layout(g, weight=None)

    color_mapping = {
        "solar": "#ffcb36",
        "wind": "white",
        "hydro": "#a5efff",
        "coal": "#686868",
        "gas": "#00ab4e",
        "": "#b6b6b6",
    }

    vertex2color = {
        k: color_mapping[v] for k, v in network["nodes"]["energy_type"].items()
    }
    v2c_list = [vertex2color[i] for i in g.nodes]  # Order based on networkx

    nodes = nx.draw_networkx_nodes(
        g,
        pos,
        node_size=250,
        node_color=v2c_list,
        linewidths=2,
    )
    edges = nx.draw_networkx_edges(
        g,
        pos,
        width=2,
        edge_color="#595959",
    )

    # Gives node colors
    ax = plt.gca()
    ax.collections[0].set_edgecolor("#595959")
    ax.set_axis_off()
visualize_network(network)
../../_images/aac15cc3a1ef019afb1a95a5356db92be86601eb12ab8c0155842cc4407e34c2.png

The IEEE-118 network has 18 generators, of which:

  • 2 hydro plants (cyan)

  • 4 coal plants (dark gray)

  • 4 solarfarms (yellow)

  • 2 windmills (white)

  • 6 gas plants (green)

Moreover, the network has 100 load nodes (gray). Each node has the following parameters:

  • node_id

  • d: demand

  • p_min: minimum production capacity

  • p_max: maximum production capacity

  • c_var: variable production costs

  • is_generator: boolean value to indicate whether the node is a generator or not

  • energy_type: the type of energy used for electricity production

All generator nodes can filtered by the is_generator parameter. All generators have a zero demand \(d_i=0\). For the renewable energy sources (i.e., hydro, solar and wind) there are no variable costs c_var. For solar and wind, the production is fixed, i.e., p_min = p_max, meaning that all available solar and wind energy must be produced.

example_nodes = network["nodes"]
example_nodes[example_nodes.is_generator]
d p_min p_max c_var is_generator energy_type
node_id
9 0.0 0.000000 400.000000 0.000000 True hydro
11 0.0 0.000000 200.000000 0.000000 True hydro
24 0.0 0.000000 397.800000 28.948321 True coal
25 0.0 0.000000 873.000000 22.220980 True coal
30 0.0 0.000000 612.000000 25.993982 True coal
45 0.0 0.000000 720.000000 24.202306 True coal
48 0.0 0.000000 0.000000 0.000000 True solar
53 0.0 0.000000 0.000000 0.000000 True solar
58 0.0 0.000000 0.000000 0.000000 True solar
60 0.0 0.000000 0.000000 0.000000 True solar
64 0.0 38.484861 38.484861 0.000000 True wind
65 0.0 180.602933 180.602933 0.000000 True wind
79 0.0 0.000000 916.200000 50.330000 True gas
86 0.0 0.000000 360.000000 50.330000 True gas
88 0.0 0.000000 1146.600000 50.330000 True gas
99 0.0 0.000000 1175.400000 50.330000 True gas
102 0.0 0.000000 194.400000 50.330000 True gas
110 0.0 0.000000 142.200000 50.330000 True gas

For the load nodes, the only important parameter is the demand \(d_i \geq 0\). All other parameters are either zero, False, or an empty string ''.

example_nodes[~example_nodes.is_generator]
d p_min p_max c_var is_generator energy_type
node_id
0 28.186145 0.0 0.0 0.0 False
1 10.921628 0.0 0.0 0.0 False
2 22.884795 0.0 0.0 0.0 False
3 21.961682 0.0 0.0 0.0 False
4 0.000000 0.0 0.0 0.0 False
... ... ... ... ... ... ...
113 4.396000 0.0 0.0 0.0 False
114 11.864042 0.0 0.0 0.0 False
115 108.311328 0.0 0.0 0.0 False
116 10.703998 0.0 0.0 0.0 False
117 18.242759 0.0 0.0 0.0 False

100 rows Ɨ 6 columns

Edges#

The IEEE-118 network has 179 edges. Each edge has the following parameters:

  • node_id1: first node of a given edge

  • node_id2: second node of a given edge

  • b: the line susceptance, and

  • f_max: the maximum edge capacity.

edges_df
b f_max
node_id1 node_id2
0 1 10.0100 271
2 23.5849 271
3 4 125.3133 316
2 4 9.2593 315
4 5 18.5185 316
... ... ... ...
64 65 28.9059 1427
67 68 28.9059 1427
80 79 28.9059 1427
86 85 4.8216 253
67 115 246.9136 12992

179 rows Ɨ 2 columns

Instances#

Since the OPF problem is solved every 15 minutes in practice, you will be given 24 * 4 = 96 network instances that need to be solved, hence solving an entire dayā€™s worth of OPF problems. The first instance thus corresponds to the power generation within time window [00:00, 00:15], the second instance corresponds to [00:15, 00:30], and so on. The data takes into account a realistic demand and (renewable energy) generation pattern. We assume that the decisions in each time window are independent of the previous and subsequent time windows, so every 15 minutes a new OPF instance is solved independently.

The network instances are stored in the variable I as a list and can be accessed using indexing, i.e., I[0] gives the first network instance and I[95] gives the 96th network instance.


Solving OPF#

Observe that the stated OPF problem contains absolute decision values. We first rewrite the OPF problem into a linear optimization problem.

\[\begin{split}\begin{align} \begin{array}{llll} \min & \sum_{i \in V} c_i p_i \\ \text{s.t.} & \sum_{j: (i, j) \in E} f^+_{ij} - f^-_{ij} - \sum_{j: (j, i) \in E} f^+_{ji} - f_{ji}^- = p_i - d_i & \forall \, i \in V,\\ & f^+_{ij} - f^-_{ij} = b_{ij}(\theta_i - \theta_j), & \forall \, (i, j) \in E, \\ & f^+_{ij} + f^-_{ij} \leq f_{ij}^{\max} & \forall (i, j) \in E,\\ & p_{i}^{\min } \leq p_{i} \leq p_{i}^{\max } & \forall i \in V, \\ & p_i \in \mathbb{R}_{\geq 0} & \forall i \in V, \\ & \theta_i \in \mathbb{R} & \forall i \in V, \\ & f_{ij}^+, f_{ij}^- \in \mathbb{R} & \forall (i, j) \in E \\ \end{array} \end{align}\end{split}\]

We then implement the model using AMPL and solve it for all instances I[0] to I[95], reporting the average objective value across all the 96 instances.

%%writefile opf1.mod

# Indexing sets
set NODES;
set EDGES within {NODES,NODES};

# Parameters
param d{NODES};
param p_min{NODES};
param p_max{NODES};
param c_var{NODES};
param is_generator{NODES};
param energy_type{NODES} symbolic;

param b{EDGES};
param f_max{EDGES};

# Declare decision variables
var p{i in NODES} >= p_min[i], <= p_max[i];
var theta{NODES};
var fabs{(i,j) in EDGES} >= 0, <= f_max[i, j];
var fp{EDGES} >= 0;
var fm{EDGES} >= 0;

# Declare objective value
minimize energy_cost: sum{i in NODES: is_generator[i] == 1} c_var[i] * p[i];

# Auxiliary variables for incoming and outgoing flows
var incoming_flow{i in NODES} = sum{j in NODES: (j,i) in EDGES} (fp[j, i] - fm[j, i]);
var outgoing_flow{i in NODES} = sum{j in NODES: (i,j) in EDGES} (fp[i, j] - fm[i, j]);

# Constraints
s.t. flow_conservation{i in NODES}:
    outgoing_flow[i] - incoming_flow[i] == p[i] - d[i];

s.t. susceptance{(i,j) in EDGES}:
    fp[i, j] - fm[i, j] == b[i, j] * (theta[i] - theta[j]);

s.t. abs_flow{(i,j) in EDGES}:
    fabs[i, j] == fp[i, j] + fm[i, j];
Overwriting opf1.mod
def OPF1(network):
    """
    Input:
    - network: a dictionary containing:
      - a dictionary of nodes with a dictionary of attributes
      - a dictionary of edges with a dictionary of attributes

    Output:
    - total energy dispatching cost
    """

    nodes = network["nodes"]
    edges = network["edges"]

    m = AMPL()
    m.read("opf1.mod")

    m.set_data(nodes, "NODES")
    m.set_data(edges, "EDGES")

    m.option["solver"] = SOLVER
    m.solve()

    return m.obj["energy_cost"].value()
OPF1_results = [OPF1(instance) for instance in I]
print(f"The average objective value over all instances is: {np.mean(OPF1_results):.2f}")
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 25556.38586
527 simplex iterations
0 barrier iterations
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 25291.14543
524 simplex iterations
0 barrier iterations
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 25014.46578
517 simplex iterations
0 barrier iterations
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 24033.92447
547 simplex iterations
0 barrier iterations
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 24607.51391
526 simplex iterations
0 barrier iterations
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 23894.80106
523 simplex iterations
0 barrier iterations
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 23539.85465
522 simplex iterations
0 barrier iterations
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 23344.3297
532 simplex iterations
0 barrier iterations
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 22863.0752
516 simplex iterations
0 barrier iterations
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 21810.52787
521 simplex iterations
0 barrier iterations
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 21900.80263
532 simplex iterations
0 barrier iterations
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 21547.3871
526 simplex iterations
0 barrier iterations
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 20462.85553
530 simplex iterations
0 barrier iterations
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 20203.29637
525 simplex iterations
0 barrier iterations
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 19631.87938
528 simplex iterations
0 barrier iterations
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 19053.61942
514 simplex iterations
0 barrier iterations
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 19052.65246
529 simplex iterations
0 barrier iterations
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 19361.24172
530 simplex iterations
0 barrier iterations
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 19838.31763
527 simplex iterations
0 barrier iterations
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 20490.77796
527 simplex iterations
0 barrier iterations
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 20961.18848
525 simplex iterations
0 barrier iterations
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 21810.71275
538 simplex iterations
0 barrier iterations
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 23204.42688
520 simplex iterations
0 barrier iterations
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 23961.68452
560 simplex iterations
0 barrier iterations
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 25710.80611
523 simplex iterations
0 barrier iterations
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 28528.05224
519 simplex iterations
0 barrier iterations
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 30606.15924
529 simplex iterations
0 barrier iterations
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 32931.97675
535 simplex iterations
0 barrier iterations
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 35933.96987
524 simplex iterations
0 barrier iterations
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 40878.89847
543 simplex iterations
0 barrier iterations
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 41272.13149
553 simplex iterations
0 barrier iterations
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 41625.78447
548 simplex iterations
0 barrier iterations
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 43657.80972
524 simplex iterations
0 barrier iterations
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 44523.08766
539 simplex iterations
0 barrier iterations
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 46653.79048
528 simplex iterations
0 barrier iterations
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 48123.99084
542 simplex iterations
0 barrier iterations
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 49986.24726
526 simplex iterations
0 barrier iterations
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 51594.08932
539 simplex iterations
0 barrier iterations
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 48767.34223
536 simplex iterations
0 barrier iterations
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 46412.44208
555 simplex iterations
0 barrier iterations
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 45234.51683
563 simplex iterations
0 barrier iterations
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 47192.70567
533 simplex iterations
0 barrier iterations
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 44287.62719
536 simplex iterations
0 barrier iterations
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 44509.58397
566 simplex iterations
0 barrier iterations
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 42934.44805
523 simplex iterations
0 barrier iterations
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 42869.82759
535 simplex iterations
0 barrier iterations
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 39555.8754
540 simplex iterations
0 barrier iterations
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 38670.71153
538 simplex iterations
0 barrier iterations
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 37154.52425
531 simplex iterations
0 barrier iterations
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 36860.74842
537 simplex iterations
0 barrier iterations
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 38740.0701
543 simplex iterations
0 barrier iterations
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 36354.1365
525 simplex iterations
0 barrier iterations
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 36912.55781
551 simplex iterations
0 barrier iterations
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 37668.23898
545 simplex iterations
0 barrier iterations
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 37040.66645
546 simplex iterations
0 barrier iterations
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 34794.47952
527 simplex iterations
0 barrier iterations
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 36325.10891
532 simplex iterations
0 barrier iterations
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 36452.65775
537 simplex iterations
0 barrier iterations
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 35163.06581
536 simplex iterations
0 barrier iterations
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 35116.35087
537 simplex iterations
0 barrier iterations
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 36112.04806
535 simplex iterations
0 barrier iterations
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 34697.75853
529 simplex iterations
0 barrier iterations
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 35180.60733
529 simplex iterations
0 barrier iterations
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 35021.77992
522 simplex iterations
0 barrier iterations
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 36530.09028
531 simplex iterations
0 barrier iterations
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 37554.629
525 simplex iterations
0 barrier iterations
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 41679.08392
542 simplex iterations
0 barrier iterations
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 46836.06918
543 simplex iterations
0 barrier iterations
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 54206.38573
541 simplex iterations
0 barrier iterations
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 57125.77445
522 simplex iterations
0 barrier iterations
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 59781.78181
554 simplex iterations
0 barrier iterations
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 62049.83588
546 simplex iterations
0 barrier iterations
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 61861.88082
547 simplex iterations
0 barrier iterations
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 62773.3659
569 simplex iterations
0 barrier iterations
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 63568.70959
548 simplex iterations
0 barrier iterations
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 63537.75177
544 simplex iterations
0 barrier iterations
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 65871.0868
554 simplex iterations
0 barrier iterations
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 64744.56999
562 simplex iterations
0 barrier iterations
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 64044.66726
549 simplex iterations
0 barrier iterations
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 62795.39073
552 simplex iterations
0 barrier iterations
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 60679.41848
543 simplex iterations
0 barrier iterations
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 60587.38654
578 simplex iterations
0 barrier iterations
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 56395.11375
541 simplex iterations
0 barrier iterations
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 54159.88066
545 simplex iterations
0 barrier iterations
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 51612.43852
553 simplex iterations
0 barrier iterations
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 47266.24936
544 simplex iterations
0 barrier iterations
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 43900.49913
540 simplex iterations
0 barrier iterations
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 40067.91939
535 simplex iterations
0 barrier iterations
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 36097.07305
551 simplex iterations
0 barrier iterations
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 34979.65852
519 simplex iterations
0 barrier iterations
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 34133.27621
516 simplex iterations
0 barrier iterations
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 33746.40968
511 simplex iterations
0 barrier iterations
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 32602.68668
514 simplex iterations
0 barrier iterations
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 31745.40199
520 simplex iterations
0 barrier iterations
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 30839.27015
528 simplex iterations
0 barrier iterations
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 29427.14541
529 simplex iterations
0 barrier iterations
 
The average objective value over all instances is: 38507.23

Strict fossil fuel policy pt.1#

Gas and coal plants emit CO2, while renewable energy sources are carbon neutral. For this reason, the Dutch government has decided to constrain the number of active fossil-fuel-based power plants for the generation of electricity. More specifically, a maximum of 2 gas plants and 1 coal plant may be active during a single OPF instance. Any plant that is set inactive for a specific instance cannot produce any electricity.

We first write down the new model. To this end, we introduce new decision variables \(x_i, i\in V\) to the model, which indicate whether a generator \(i\) is active or inactive.

\[\begin{split}\begin{align} \begin{array}{llll} \min & \sum_{i \in V} c_i p_i \\ \text{s.t.} & \sum_{j: (i, j) \in E} f^+_{ij} - f^-_{ij} - \sum_{j: (j, i) \in E} f^+_{ji} - f_{ji}^- = p_i - d_i & \forall \, i \in V,\\ & f^+_{ij} - f^-_{ij} = b_{ij}(\theta_i - \theta_j), & \forall \, (i, j) \in E, \\ & f^{abs}_{ij} = f^+_{ij} + f^-_{ij}, & \forall \, (i, j) \in E, \\ & f_{ij}^{abs} \leq f_{ij}^{\max} & \forall (i, j) \in E,\\ & p_{i}^{\min } x_i \leq p_{i} \leq p_{i}^{\max } x_i & \forall i \in V, \\ & \sum_{i \in \mathcal{G}_{\text{gas}}} x_i \leq 2 & \\ & \sum_{i \in \mathcal{G}_{\text{coal}}} x_i \leq 1 & \\ & p_i \in \mathbb{R}_{\geq 0} & \forall i \in V, \\ & \theta_i \in \mathbb{R} & \forall i \in V, \\ & f_{ij} \in \mathbb{R} & \forall (i, j) \in E \\ & x_i \in \{0, 1\} & \forall i \in V\\ \end{array} \end{align}\end{split}\]

We then implement the new model using AMPL and solve it for all instances I[0] to I[95], reporting the average objective value across the instances.

%%writefile opf2.mod

# Indexing sets
set NODES;
set EDGES within {NODES,NODES};

# Parameters
param d{NODES};
param p_min{NODES};
param p_max{NODES};
param c_var{NODES};
param is_generator{NODES};
param energy_type{NODES} symbolic;

param b{EDGES};
param f_max{EDGES};

param max_gas_plants;
param max_coal_plants;

# Declare decision variables
var p{i in NODES} >= 0;
var theta{NODES};
var x{NODES} binary;
var fabs{(i,j) in EDGES} >= 0, <= f_max[i, j];
var fp{EDGES} >= 0;
var fm{EDGES} >= 0;

# Declare objective value
minimize energy_cost: sum{i in NODES: is_generator[i] == 1} c_var[i] * p[i];

# Auxiliary variables for incoming and outgoing flows
var incoming_flow{i in NODES} = sum{j in NODES: (j,i) in EDGES} (fp[j, i] - fm[j, i]);
var outgoing_flow{i in NODES} = sum{j in NODES: (i,j) in EDGES} (fp[i, j] - fm[i, j]);

# Constraints
s.t. flow_conservation{i in NODES}:
    outgoing_flow[i] - incoming_flow[i] == p[i] - d[i];

s.t. susceptance{(i,j) in EDGES}:
    fp[i, j] - fm[i, j] == b[i, j] * (theta[i] - theta[j]);

s.t. abs_flow{(i,j) in EDGES}:
    fabs[i, j] == fp[i, j] + fm[i, j];

s.t. generation_upper_bound{i in NODES}: p[i] <= p_max[i] * x[i];
s.t. generation_lower_bound{i in NODES}: p[i] >= p_min[i] * x[i];

s.t. gas_plants_limit: sum{i in NODES: energy_type[i] == 'gas'} x[i] <= max_gas_plants;
s.t. coal_plants_limit: sum{i in NODES: energy_type[i] == 'coal'} x[i] <= max_coal_plants;
Overwriting opf2.mod
def OPF2(network, max_gas_plants, max_coal_plants):
    """
    Input:
    - network: a dictionary containing:
      - a dictionary of nodes with a dictionary of attributes
      - a dictionary of edges with a dictionary of attributes

    Output:
    - total energy dispatching cost
    """

    nodes = network["nodes"]
    edges = network["edges"]

    m = AMPL()
    m.read("opf2.mod")

    m.set_data(nodes, "NODES")
    m.set_data(edges, "EDGES")

    m.param["max_gas_plants"] = max_gas_plants
    m.param["max_coal_plants"] = max_coal_plants

    m.option["solver"] = SOLVER
    m.solve()

    return m.obj["energy_cost"].value()
max_gas_plants = 2
max_coal_plants = 1
OPF2_results = [OPF2(instance, max_gas_plants, max_coal_plants) for instance in I]
print(f"The average objective value over all instances is: {np.mean(OPF2_results):.2f}")
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 32014.85844
578 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 31513.92689
576 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 30885.43208
569 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 28881.74194
558 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 30037.14784
553 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 28659.49442
730 simplex iterations
3 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 28098.31029
706 simplex iterations
3 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 27643.13369
737 simplex iterations
3 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 26906.38333
666 simplex iterations
3 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 25000.0886
788 simplex iterations
3 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 25172.05815
774 simplex iterations
3 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 24597.59628
773 simplex iterations
3 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 22704.40086
810 simplex iterations
3 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 22272.96622
802 simplex iterations
3 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 21282.46505
794 simplex iterations
3 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 20149.63104
768 simplex iterations
3 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 20299.75691
747 simplex iterations
3 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 20796.63521
766 simplex iterations
3 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 21606.25852
820 simplex iterations
3 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 22744.74502
831 simplex iterations
3 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 23499.3076
835 simplex iterations
3 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 25010.50854
770 simplex iterations
3 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 27655.6242
725 simplex iterations
3 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 28739.85416
674 simplex iterations
3 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 32306.05117
570 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 38259.59372
592 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 42596.56579
572 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 47474.46445
587 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 52182.27728
549 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 57946.94606
549 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 58983.93317
556 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 59611.09337
562 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 62318.4829
610 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 63377.649
553 simplex iterations
1 branching nodes
absmipgap=7.27596e-12, relmipgap=1.14803e-16
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 66355.25327
648 simplex iterations
1 branching nodes
absmipgap=2.91038e-11, relmipgap=4.38606e-16
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 67718.3282
559 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 70199.13353
601 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 72074.24181
574 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 69357.83301
587 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 66541.23238
554 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 65359.14645
580 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 67603.72994
551 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 64424.45835
543 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 64578.61294
550 simplex iterations
1 branching nodes
absmipgap=8.73115e-11, relmipgap=1.35202e-15
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 62683.24322
544 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 62712.94063
537 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 58810.9962
539 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 57847.89312
533 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 54468.61779
534 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 54196.02103
534 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 56244.40675
538 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 52907.81025
515 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 53769.82573
558 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 55263.93141
549 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 54378.86264
545 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 50680.88639
535 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 53087.4615
541 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 52929.63938
539 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 51552.75812
524 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 51058.5326
534 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 52840.59651
545 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 50645.90154
522 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 51995.97075
554 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 51721.08012
537 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 54156.67505
542 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 55970.40312
544 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 61124.14254
537 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 67438.51907
557 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 75846.37597
606 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 79266.53736
570 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 82505.30565
580 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 85139.05154
577 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 84584.98267
590 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 85269.31644
592 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 86348.42938
573 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 86331.61954
568 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 88598.84186
630 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 87298.4441
595 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 86180.15407
560 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 84762.76321
569 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 82124.08662
549 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 82071.14498
575 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 77315.23968
564 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 74464.39753
582 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 71600.93405
609 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 66396.59052
621 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 62377.56686
564 simplex iterations
1 branching nodes
absmipgap=1.16415e-10, relmipgap=1.8663e-15
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 57665.26521
562 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 53248.29462
553 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 51707.5513
547 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 50039.81283
581 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 49235.30303
584 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 46856.8694
595 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 45074.09967
569 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 43164.18792
577 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 40195.99971
560 simplex iterations
1 branching nodes
 
The average objective value over all instances is: 53120.81

Strict fossil fuel policy pt.2#

The restriction on the number of gas and coal plants may pose a threat to the availability of electricity production when renewable energy sources fail to deliver. For this reason, the grid operators have decided to slightly change the constraint that was introduced for OPF2. If the total production of energy from renewable energy sources (i.e., solar, wind and hydro) is above 1000, then the number of gas and coal plants is restricted to 2 and 1, respectively. Otherwise, the restriction on the number of gas and coal plants is lifted. These constraints can be modeled as either-or constraints.

We first write down the new model, using big-\(M\) constraints to model the either-or constraint.

\[\begin{split}\begin{align} \begin{array}{llll} \min & \sum_{i \in V} c_i p_i \\ \text{s.t.} & \sum_{j: (i, j) \in E} f^+_{ij} - f^-_{ij} - \sum_{j: (j, i) \in E} f^+_{ji} - f_{ji}^- = p_i - d_i & \forall \, i \in V,\\ & f^+_{ij} - f^-_{ij} = b_{ij}(\theta_i - \theta_j), & \forall \, (i, j) \in E, \\ & f^{abs}_{ij} = f^+_{ij} + f^-_{ij}, & \forall \, (i, j) \in E, \\ & f_{ij}^{abs} \leq f_{ij}^{\max} & \forall (i, j) \in E,\\ & p_{i}^{\min } x_i \leq p_{i} \leq p_{i}^{\max } x_i & \forall i \in V, \\ & \sum_{i \in \mathcal{G}_{\text{gas}}} x_i \leq 2 + (1-y)M_0 & \\ & \sum_{i \in \mathcal{G}_{\text{coal}}} x_i \leq 1 + (1-y)M_1 & \\ & \sum_{i \in \mathcal{G}_\text{solar}} p_i + \sum_{i \in \mathcal{G}_\text{wind}} p_i + \sum_{i \in \mathcal{G}_\text{hydro}} p_i \leq 1000 + yM_2&\\ & p_i \in \mathbb{R}_{\geq 0} & \forall i \in V, \\ & \theta_i \in \mathbb{R} & \forall i \in V, \\ & f_{ij} \in \mathbb{R} & \forall (i, j) \in E \\ & x_i \in \{0, 1\} & \forall i \in V\\ & y \in \{0, 1\} &\\ \end{array} \end{align}\end{split}\]

We now implement the new model using AMPL and solve it for all instances I[0] to I[95], reporting the average objective value across the instances.

%%writefile opf3.mod

# Indexing sets
set NODES;
set EDGES within {NODES,NODES};

# Parameters
param d{NODES};
param p_min{NODES};
param p_max{NODES};
param c_var{NODES};
param is_generator{NODES};
param energy_type{NODES} symbolic;

param b{EDGES};
param f_max{EDGES};

param max_gas_plants;
param max_coal_plants;

# Big-Ms
param M0;
param M1;
param M2;

# Declare decision variables
var p{i in NODES} >= 0;
var theta{NODES};
var x{NODES} binary;
var fabs{(i,j) in EDGES} >= 0, <= f_max[i, j];
var fp{EDGES} >= 0;
var fm{EDGES} >= 0;
var y binary;

# Declare objective value
minimize energy_cost: sum{i in NODES: is_generator[i] == 1} c_var[i] * p[i];

# Auxiliary variables for incoming and outgoing flows
var incoming_flow{i in NODES} = sum{j in NODES: (j,i) in EDGES} (fp[j, i] - fm[j, i]);
var outgoing_flow{i in NODES} = sum{j in NODES: (i,j) in EDGES} (fp[i, j] - fm[i, j]);

# Constraints
s.t. flow_conservation{i in NODES}:
    outgoing_flow[i] - incoming_flow[i] == p[i] - d[i];

s.t. susceptance{(i,j) in EDGES}:
    fp[i, j] - fm[i, j] == b[i, j] * (theta[i] - theta[j]);

s.t. abs_flow{(i,j) in EDGES}:
    fabs[i, j] == fp[i, j] + fm[i, j];

s.t. generation_upper_bound{i in NODES}: p[i] <= p_max[i] * x[i];
s.t. generation_lower_bound{i in NODES}: p[i] >= p_min[i] * x[i];

s.t. gas_plants_limit: sum{i in NODES: energy_type[i] == 'gas'} x[i] <= max_gas_plants + (1 - y) * M0;

s.t. coal_plants_limit: sum{i in NODES: energy_type[i] == 'coal'} x[i] <= max_coal_plants + (1 - y) * M1;

s.t. renewable_energy_production:
    sum{i in NODES: energy_type[i] in {'solar', 'wind', 'hydro'}} p[i] <= 1000 + y * M2;
Overwriting opf3.mod
def OPF3(network, max_gas_plants, max_coal_plants):
    """
    Input:
    - network: a dictionary containing:
      - a dictionary of nodes with a dictionary of attributes
      - a dictionary of edges with a dictionary of attributes

    Output:
    - total energy dispatching cost
    """

    nodes = network["nodes"]
    edges = network["edges"]

    m = AMPL()
    m.read("opf3.mod")

    m.set_data(nodes, "NODES")
    m.set_data(edges, "EDGES")

    m.param["max_gas_plants"] = max_gas_plants
    m.param["max_coal_plants"] = max_coal_plants

    m.param["M0"] = 4
    m.param["M1"] = 3
    m.param["M2"] = nodes["p_max"].sum()

    m.option["solver"] = SOLVER
    m.solve()

    return m.obj["energy_cost"].value()
OPF3_results = [OPF3(instance, max_gas_plants, max_coal_plants) for instance in I]
print(f"The average objective value over all instances is: {np.mean(OPF3_results):.2f}")
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 25556.38586
542 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 25291.14543
541 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 25014.46578
532 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 24033.92447
532 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 24607.51391
527 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 23894.80106
548 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 23539.85465
531 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 23344.3297
555 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 22863.0752
538 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 21810.52787
535 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 21900.80263
531 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 21547.3871
524 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 20462.85553
543 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 20203.29637
534 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 19631.87938
534 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 19053.61942
535 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 19052.65246
539 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 19361.24172
535 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 19838.31763
539 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 20490.77796
531 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 20961.18848
531 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 21810.71275
536 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 23204.42688
535 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 23961.68452
533 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 25710.80611
521 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 28528.05224
533 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 30606.15924
540 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 32931.97675
518 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 35933.96987
563 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 40878.89847
540 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 41272.13149
551 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 41625.78447
530 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 43657.80972
556 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 44523.08766
522 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 46653.79048
531 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 48517.02462
610 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 51469.76795
793 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 53735.57811
667 simplex iterations
3 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 52478.34218
792 simplex iterations
3 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 51082.21976
789 simplex iterations
3 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 50820.70479
650 simplex iterations
3 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 52284.09621
607 simplex iterations
3 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 51770.03428
722 simplex iterations
3 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 52556.32959
1132 simplex iterations
3 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 52063.08667
934 simplex iterations
9 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 50973.40711
968 simplex iterations
3 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 49101.46232
1292 simplex iterations
7 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 47859.57387
1175 simplex iterations
3 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 47812.22458
1084 simplex iterations
10 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 47735.03354
1065 simplex iterations
7 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 49249.56069
1073 simplex iterations
10 branching nodes
absmipgap=2.11363, relmipgap=4.29168e-05
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 47079.99596
1159 simplex iterations
11 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 47106.46001
1103 simplex iterations
11 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 47698.24825
1418 simplex iterations
8 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 47740.93595
1170 simplex iterations
10 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 45854.74172
1287 simplex iterations
14 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 46738.37993
1057 simplex iterations
11 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 46817.47832
1160 simplex iterations
14 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 45851.54466
1257 simplex iterations
14 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 46189.99144
1222 simplex iterations
14 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 46339.28074
1165 simplex iterations
14 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 45723.75827
1106 simplex iterations
14 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 44508.52721
1083 simplex iterations
12 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 44574.67924
903 simplex iterations
12 branching nodes
absmipgap=0.500375, relmipgap=1.12255e-05
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 46015.12654
1007 simplex iterations
3 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 46672.0612
793 simplex iterations
3 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 49744.16264
1167 simplex iterations
7 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 54261.30821
631 simplex iterations
3 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 59898.76421
605 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 62696.04055
587 simplex iterations
1 branching nodes
absmipgap=1.00408e-09, relmipgap=1.60151e-14
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 63259.92706
747 simplex iterations
3 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 65614.15585
689 simplex iterations
3 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 64099.97548
730 simplex iterations
3 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 63826.39248
584 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 63624.76669
580 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 63537.75177
570 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 65871.0868
554 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 64744.56999
564 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 64044.66726
549 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 62795.39073
551 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 60679.41848
543 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 60587.38654
535 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 56395.11375
546 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 54159.88066
567 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 51612.43852
552 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 47266.24936
536 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 43900.49913
535 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 40067.91939
529 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 36097.07305
527 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 34979.65852
522 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 34133.27621
524 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 33746.40968
544 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 32602.68668
521 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 31745.40199
547 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 30839.27015
533 simplex iterations
1 branching nodes
 
HiGHS 1.5.3: HiGHS 1.5.3: optimal solution; objective 29427.14541
535 simplex iterations
1 branching nodes
 
The average objective value over all instances is: 41608.73

Comparing the three models#

For all three implemented models, we plot the objective values for all instances and explain the differences between the objective values in view of the different feasible regions.

objs = [OPF1_results, OPF2_results, OPF3_results]
plt.plot(objs[0], color="blue", label="OPF1")
plt.plot(objs[1], color="green", label="OPF2")
plt.plot(objs[2], color="orange", label="OPF3")
plt.title("Optimal objective values over all instances and models")
plt.xlabel("Instance number")
plt.ylabel("Optimal objective value")
plt.legend()
plt.show()
../../_images/fcc62a85fb344ffb1d63498812a9a316a78f32078f9f9dd5a6da978859b4d74b.png

The goal of the OPF problem is to minimize the total costs of dispatching energy. OPF1 is the original OPF formulation, whereas OPF2 and OPF3 are restricted versions of OPF1. This means that the feasible region of OPF1 is at least as large as OPF2 and OPF3, where we may assume that \(x_i = 1\) for all the generators \(i\).

If we let \(F_1, F_2, F_3\) denote the feasible region of OPF1, OPF2 and OPF3, then we observe that which explains that the optimal objective value of OPF1 always remains below the optimal objectives values of OPF2 and OPF3.

The most important observation to make is that based on the problem descriptions, we would expect OPF2 to take on the objective values of OPF1 or OPF3, depending on which of the either-or-constraints are activated.

  • OPF3 uses expensive generators, and possibly not all the renewable energy

  • OPF2 does only uses 1000 renewable energy power, because then it may keep using all the gas and coal generators.

  • OPF1 uses all renewable energy at all times, because it has the flexibility to use all generators in order to mitigate operational restrictions due to line flows.