Electric Grid Optimization using QAOA

We have a set of N power plants (sources), and M consumers. The goal is to supply power to all consumers, meeting the constraints of the power plants, and minimize the total cost of supplying power. The presented model has small variation from [1].

Mathematical model, minimize the objective function:

$z = \sum_{i=1}^{n} \sum_{j=1}^{m} Z_{ij}x_{ij}$

where $$x_{ij}$$ is the required values of the transmitted power from source $$A_i$$ to consumer $$B_j$$.

The unit cost of transmitting power from node $$A_i$$ to node $$B_j$$ is $$Z_{ij}$$.

Constraint to:

The sum of powers flowing from power plant transmission lines to all customer nodes must be up to the power of the source $$A_i$$

$\sum_{j=1}^{M} x_{ij} \leq A_{i} \quad i=1,2,...,N$

Each consumer recieve power $$B_{j}$$:

$\sum_{i=1}^{N} x_{ij} = B_{j} \quad j=1,2,...,M$

We take $$B_{j} = 1$$ and $$A_{i} = 2$$ in the following example.

Notice that we use 2 kinds of constraints, equality and inequality.

Building the problem

from itertools import product

import matplotlib.pyplot as plt
import networkx as nx  # noqa
import numpy as np
import pandas as pd

# building data matrix, it doesn't need to be a symmetric matrix.
cost_matrix = np.array(
[[0.5, 1.0, 1.0, 2.1], [1.0, 0.6, 1.4, 1.0], [1.0, 1.4, 0.4, 2.3]]
)

Sources = ["A1", "A2", "A3"]
Consumers = ["B1", "B2", "B3", "B4"]

# number of sources
N = len(Sources)
# number of consumers
M = len(Consumers)

graph = nx.DiGraph()
for n, m in product(range(N), range(M)):

# Plot the graph
plt.figure(figsize=(10, 6))
left = nx.bipartite.sets(graph)[0]
pos = nx.bipartite_layout(graph, left)

nx.draw_networkx(graph, pos=pos, nodelist=Consumers, font_size=22, font_color="None")
nx.draw_networkx_nodes(
graph, pos, nodelist=Consumers, node_color="#119DA4", node_size=500
)
for fa in Sources:
x, y = pos[fa]
plt.text(
x,
y,
s=fa,
bbox=dict(facecolor="#F43764", alpha=1),
horizontalalignment="center",
fontsize=15,
)

nx.draw_networkx_edges(graph, pos, width=2)
labels = nx.get_edge_attributes(graph, "weight")
nx.draw_networkx_edge_labels(graph, pos, edge_labels=labels, font_size=12)
nx.draw_networkx_labels(
graph,
pos,
labels={co: co for co in Consumers},
font_size=15,
font_color="#F4F9E9",
)

plt.axis("off")
plt.show()


building the pyomo model - classical combinatorial optimization problem

from typing import List, Tuple, cast  # noqa

import pyomo.environ as pyo
from IPython.display import Markdown, display

opt_model = pyo.ConcreteModel()

sources_lst = range(N)
consumers_lst = range(M)

opt_model.x = pyo.Var(sources_lst, consumers_lst, domain=pyo.Binary)

@opt_model.Constraint(sources_lst)
def source_supply_rule(model, n):  # constraint (1)
return sum(model.x[n, m] for m in consumers_lst) <= 2

@opt_model.Constraint(consumers_lst)
def each_consumer_is_supplied_rule(model, m):  # constraint (2)
return sum(model.x[n, m] for n in sources_lst) == 1

opt_model.cost = pyo.Objective(
expr=sum(
cost_matrix[n, m] * opt_model.x[n, m]
for n in sources_lst
for m in consumers_lst
),
sense=pyo.minimize,
)


Printing the classical optimization problem

# opt_model.pprint()


Solving with Classiq

We take a specific example: the one outlined above.

First generating parameters for the qauntum circuit

from classiq import (
Preferences,
construct_combinatorial_optimization_model,
set_preferences,
)
from classiq.applications.combinatorial_optimization import OptimizerConfig, QAOAConfig

qaoa_config = QAOAConfig(num_layers=5, penalty_energy=3.0)
# the penalty_energy influence the convergence to the right solution.


For the classical optimization part of the QAOA algorithm (OptimizerConfig) we define the maximum number of iterations of the hybrid algorithm(max_iteration) and the $$\alpha$$-parameter (alpha_cvar) for running CVaR-QAOA, an improved variation of the QAOA algorithm:

optimizer_config = OptimizerConfig(
max_iteration=100,
alpha_cvar=1,
)


Combining all the parts together using construct_combinatorial_optimization_model to get a SeralizedModel which include pyomo.ConcreteModel, QAOAConfig, and OptimizerConfig.

qmod = construct_combinatorial_optimization_model(
pyo_model=opt_model,
qaoa_config=qaoa_config,
optimizer_config=optimizer_config,
)

# defining cosntraint such as computer and parameters for a quicker and more optimized circuit.
preferences = Preferences(timeout_seconds=3000)

qmod = set_preferences(qmod, preferences)


Now we can create a quantum circuit using the synthesize command and show it

from classiq import write_qmod

write_qmod(qmod, "electric_grid_optimization")

from classiq import show, synthesize

qprog = synthesize(qmod)
show(qprog)

Opening: https://platform.classiq.io/circuit/f21658ea-a358-4aad-a5e8-096d4c16bed7?version=0.41.0.dev39%2B79c8fd0855


We now execute the problem iteratively by calling the execute function on the quantum program we have generated:

from classiq import execute

res = execute(qprog).result()


We can check the convergence of the run:

from classiq.execution import VQESolverResult

vqe_result = res[0].value
vqe_result.convergence_graph


Get best solution statistics

import pandas as pd

from classiq.applications.combinatorial_optimization import (
get_optimization_solution_from_pyo,
)

solution = get_optimization_solution_from_pyo(
opt_model, vqe_result=vqe_result, penalty_energy=qaoa_config.penalty_energy
)

optimization_result = pd.DataFrame.from_records(solution)

probability cost solution count
160 0.001 3.6 [1, 0, 0, 1, 0, 1, 0, 0, 0, 0, 1, 0] 1
0 0.003 4.1 [0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 1, 0] 3
435 0.001 4.1 [0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 1, 0] 1
5 0.002 4.1 [0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 1, 0] 2
802 0.001 4.3 [0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 1, 1] 1
optimization_result["cost"].plot(
kind="hist", bins=30, edgecolor="black", weights=optimization_result["probability"]
)
plt.ylabel("Probability", fontsize=16)
plt.xlabel("cost", fontsize=16)
plt.tick_params(axis="both", labelsize=14)


Best solution visulalization

# This function plots the solution in a table and a graph

def plotting_sol(x_sol, cost, is_classic: bool):
x_sol_to_mat = np.reshape(np.array(x_sol), [N, M])  # vector to matrix
# opened facilities will be marked in red
opened_fac_dict = {}
for fa in range(N):
if sum(x_sol_to_mat[fa, m] for m in range(M)) > 0:
opened_fac_dict.update({Sources[fa]: "background-color: #F43764"})

# classical or quantum
if is_classic == True:
display(Markdown("**CLASSICAL SOLUTION**"))
print("total cost= ", cost)
else:
display(Markdown("**QAOA SOLUTION**"))
print("total cost= ", cost)

# plotting in a table
df = pd.DataFrame(x_sol_to_mat)
df.columns = Consumers
df.index = Sources
plotable = df.style.apply(lambda x: x.index.map(opened_fac_dict))
display(plotable)

# plotting in a graph
graph_sol = nx.DiGraph()
for n, m in product(range(N), range(M)):
if x_sol_to_mat[n, m] > 0:
[(Sources[n], Consumers[m])], weight=cost_matrix[n, m]
)

plt.figure(figsize=(10, 6))
left = nx.bipartite.sets(graph_sol, top_nodes=Sources)[0]
pos = nx.bipartite_layout(graph_sol, left)

nx.draw_networkx(
graph_sol, pos=pos, nodelist=Consumers, font_size=22, font_color="None"
)
nx.draw_networkx_nodes(
graph_sol, pos, nodelist=Consumers, node_color="#119DA4", node_size=500
)
for fa in Sources:
x, y = pos[fa]
if fa in opened_fac_dict.keys():
plt.text(
x,
y,
s=fa,
bbox=dict(facecolor="#F43764", alpha=1),
horizontalalignment="center",
fontsize=15,
)
else:
plt.text(
x,
y,
s=fa,
bbox=dict(facecolor="#F4F9E9", alpha=1),
horizontalalignment="center",
fontsize=15,
)

nx.draw_networkx_edges(graph_sol, pos, width=2)
labels = nx.get_edge_attributes(graph_sol, "weight")
nx.draw_networkx_edge_labels(graph, pos, edge_labels=labels, font_size=12)
nx.draw_networkx_labels(
graph_sol,
pos,
labels={co: co for co in Consumers},
font_size=15,
font_color="#F4F9E9",
)

plt.axis("off")
plt.show()

best_solution = optimization_result.loc[optimization_result.cost.idxmin()]

plotting_sol(best_solution.solution, best_solution.cost, is_classic=False)


QAOA SOLUTION

total cost=  3.599999999999998

B1 B2 B3 B4
A1 1 0 0 1
A2 0 1 0 0
A3 0 0 1 0

Compare to a classical solver

from pyomo.opt import SolverFactory

solver = SolverFactory("couenne")
solver.solve(opt_model)

# opt_model.display()

{'Problem': [{'Lower bound': -inf, 'Upper bound': inf, 'Number of objectives': 1, 'Number of constraints': 0, 'Number of variables': 12, 'Sense': 'unknown'}], 'Solver': [{'Status': 'ok', 'Message': 'couenne\\x3a Optimal', 'Termination condition': 'optimal', 'Id': 3, 'Error rc': 0, 'Time': 0.016887187957763672}], 'Solution': [OrderedDict([('number of solutions', 0), ('number of solutions displayed', 0)])]}

best_classical_solution = np.array(
[pyo.value(opt_model.x[idx]) for idx in np.ndindex(cost_matrix.shape)]
).reshape(cost_matrix.shape)

plotting_sol(
np.round([pyo.value(opt_model.x[idx]) for idx in np.ndindex(cost_matrix.shape)]),
pyo.value(opt_model.cost),
is_classic=True,
)


CLASSICAL SOLUTION

total cost=  2.4999999996418167

B1 B2 B3 B4
A1 1.000000 0.000000 0.000000 0.000000
A2 0.000000 1.000000 0.000000 1.000000
A3 0.000000 -0.000000 1.000000 0.000000

[1]: Solving optimization problems when designing power supply circuits. https://www.e3s-conferences.org/articles/e3sconf/pdf/2019/50/e3sconf_ses18_04011.pdf