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¶

In [ ]:

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()
graph.add_nodes_from(Sources + Consumers)
for n, m in product(range(N), range(M)):
    graph.add_edges_from([(Sources[n], Consumers[m])], weight=cost_matrix[n, m])


# Plot the graph
plt.figure(figsize=(6, 5))
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.tight_layout()
plt.axis("off")
plt.show()

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() graph.add_nodes_from(Sources + Consumers) for n, m in product(range(N), range(M)): graph.add_edges_from([(Sources[n], Consumers[m])], weight=cost_matrix[n, m]) # Plot the graph plt.figure(figsize=(6, 5)) 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.tight_layout() plt.axis("off") plt.show()

building the pyomo model - classical combinatorial optimization problem¶

In [ ]:

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,
)

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

In [ ]:

# opt_model.pprint()

# opt_model.pprint()

Solving with Classiq¶

We take a specific example: the one outlined above.

First generating parameters for the qauntum circuit¶

In [ ]:

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.

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:

In [ ]:

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

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.

In [ ]:

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(transpilation_option="none", timeout_seconds=3000)

qmod = set_preferences(qmod, preferences)

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(transpilation_option="none", timeout_seconds=3000) qmod = set_preferences(qmod, preferences)

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

In [ ]:

with open("electric_grid_optimization.qmod", "w") as f:
    f.write(qmod)

with open("electric_grid_optimization.qmod", "w") as f: f.write(qmod)

In [ ]:

from classiq import show, synthesize

qprog = synthesize(qmod)
show(qprog)

from classiq import show, synthesize qprog = synthesize(qmod) show(qprog)

We now execute the problem iteratively using the generated circuit by using the execute method:

In [ ]:

from classiq import execute

res = execute(qprog).result()

from classiq import execute res = execute(qprog).result()

We can check the convergence of the run:

In [ ]:

from classiq.execution import VQESolverResult

vqe_result = res[1].value
vqe_result.convergence_graph

from classiq.execution import VQESolverResult vqe_result = res[1].value vqe_result.convergence_graph

Get best solution statistics¶

In [ ]:

import pandas as pd

optimization_result = pd.DataFrame.from_records(res[0].value)
optimization_result.sort_values(by="cost", ascending=True).head(5)

import pandas as pd optimization_result = pd.DataFrame.from_records(res[0].value) optimization_result.sort_values(by="cost", ascending=True).head(5)

In [ ]:

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)

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¶

In [ ]:

# 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()
    graph_sol.add_nodes_from(Sources + Consumers)
    for n, m in product(range(N), range(M)):
        if x_sol_to_mat[n, m] > 0:
            graph_sol.add_edges_from(
                [(Sources[n], Consumers[m])], weight=cost_matrix[n, m]
            )

    plt.figure(figsize=(6, 5))
    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.tight_layout()
    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)

# 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() graph_sol.add_nodes_from(Sources + Consumers) for n, m in product(range(N), range(M)): if x_sol_to_mat[n, m] > 0: graph_sol.add_edges_from( [(Sources[n], Consumers[m])], weight=cost_matrix[n, m] ) plt.figure(figsize=(6, 5)) 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.tight_layout() 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)

Compare to a classical solver¶

In [ ]:

from pyomo.opt import SolverFactory

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

# opt_model.display()

from pyomo.opt import SolverFactory solver = SolverFactory("couenne") solver.solve(opt_model) # opt_model.display()

In [ ]:

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,
)

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, )

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