Skip to main content

Documentation Index

Fetch the complete documentation index at: https://prod-mint.classiq.io/llms.txt

Use this file to discover all available pages before exploring further.

View on GitHub

Open this notebook in GitHub to run it yourself
import random

import numpy as np
import torch

random.seed(8)
np.random.seed(8)
torch.manual_seed(8)
Output:
<torch._C.

Generator at 0x7135940974b0>

​
Electric Grid Optimization using QAOA

For a set of N power plants (sources) and M consumers, the goal is to supply power to all consumers while meeting the constraints of the power plants and minimizing the total cost of supplying power. The model here is a minor variation of [1]. Mathematical model, minimizing the objective function: z=∑i=1n∑j=1mZijxijz = \sum_{i=1}^{n} \sum_{j=1}^{m} Z_{ij}x_{ij} where xijx_{ij} is the required values of the transmitted power from source AiA_i to consumer BjB_j. The unit cost of transmitting power from node AiA_i to node BjB_j is ZijZ_{ij}. Constraint: the sum of powers flowing from power plant transmission lines to all customer nodes must be up to the power of the source AiA_i: ∑j=1Mxij≤Aii=1,2,...,N\sum_{j=1}^{M} x_{ij} \leq A_{i} \quad i=1,2,...,N Each consumer receives power BjB_{j}: ∑i=1Nxij=Bjj=1,2,...,M\sum_{i=1}^{N} x_{ij} = B_{j} \quad j=1,2,...,M This example takes Bj=1B_{j} = 1 and Ai=2A_{i} = 2. Note the use of two 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()
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=(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()
output Build the Pyomo mjodel for a classical combinatorial optimization problem: 
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,
)
Print the classical optimization problem: 
opt_model.pprint()
Output:
1 Var Declarations
      x : Size=12, Index={0, 1, 2}*{0, 1, 2, 3}
          Key    : Lower : Value : Upper : Fixed : Stale : Domain
          (0, 0) :     0 :  None :     1 : False :  True : Binary
          (0, 1) :     0 :  None :     1 : False :  True : Binary
          (0, 2) :     0 :  None :     1 : False :  True : Binary
          (0, 3) :     0 :  None :     1 : False :  True : Binary
          (1, 0) :     0 :  None :     1 : False :  True : Binary
          (1, 1) :     0 :  None :     1 : False :  True : Binary
          (1, 2) :     0 :  None :     1 : False :  True : Binary
          (1, 3) :     0 :  None :     1 : False :  True : Binary
          (2, 0) :     0 :  None :     1 : False :  True : Binary
          (2, 1) :     0 :  None :     1 : False :  True : Binary
          (2, 2) :     0 :  None :     1 : False :  True : Binary
          (2, 3) :     0 :  None :     1 : False :  True : Binary

  1 Objective Declarations
      cost : Size=1, Index=None, Active=True
          Key  : Active : Sense    : Expression
          None :   True : minimize : 0.5*x[0,0] + x[0,1] + x[0,2] + 2.1*x[0,3] + x[1,0] + 0.6*x[1,1] + 1.4*x[1,2] + x[1,3] + x[2,0] + 1.4*x[2,1] + 0.4*x[2,2] + 2.3*x[2,3]

  2 Constraint Declarations
      each_consumer_is_supplied_rule : Size=4, Index={0, 1, 2, 3}, Active=True
          Key : Lower : Body                     : Upper : Active
            0 :   1.0 : x[0,0] + x[1,0] + x[2,0] :   1.0 :   True
            1 :   1.0 : x[0,1] + x[1,1] + x[2,1] :   1.0 :   True
            2 :   1.0 : x[0,2] + x[1,2] + x[2,2] :   1.0 :   True
            3 :   1.0 : x[0,3] + x[1,3] + x[2,3] :   1.0 :   True
      source_supply_rule : Size=3, Index={0, 1, 2}, Active=True
          Key : Lower : Body                              : Upper : Active
            0 :  -Inf : x[0,0] + x[0,1] + x[0,2] + x[0,3] :   2.0 :   True
            1 :  -Inf : x[1,0] + x[1,1] + x[1,2] + x[1,3] :   2.0 :   True
            2 :  -Inf : x[2,0] + x[2,1] + x[2,2] + x[2,3] :   2.0 :   True

  4 Declarations: x source_supply_rule each_consumer_is_supplied_rule cost

​
Solving with Classiq

Take the specific example outlined above.

​

​
Generating Parameters for the Quantum Circuit

from classiq import *
from classiq.applications.combinatorial_optimization import CombinatorialProblem

combi = CombinatorialProblem(pyo_model=opt_model, num_layers=4, penalty_factor=10)

qmod = combi.get_model()

​

​
Synthesizing the QAOA Circuit and Solving the Problem

Synthesize and view the QAOA circuit (ansatz) used to solve the optimization problem: 
qprog = combi.get_qprog()
show(qprog)
Output:

Quantum program link: https://platform.classiq.io/circuit/39Z3KwseLb8I9dfXZ3jjAPvBXrL
Output:
https://platform.classiq.io/circuit/39Z3KwseLb8I9dfXZ3jjAPvBXrL?login=True&version=17
Solve the problem by calling the optimize method of the CombinatorialProblem object. For the classical optimization part of QAOA, define the maximum number of classical iterations (maxiter) and the α\alpha-parameter (quantile) for running CVaR-QAOA, an improved variation of the QAOA algorithm [3]: 
optimized_params = combi.optimize(maxiter=100, quantile=1)
Check the convergence of the run: 
import matplotlib.pyplot as plt

plt.plot(combi.cost_trace)
plt.xlabel("Iterations")
plt.ylabel("Cost")
plt.title("Cost convergence")
Output:

Text(0.5, 1.0, 'Cost convergence')
output

​
Best Solution Statistics

optimization_result = combi.sample(optimized_params)
optimization_result.sort_values(by="cost").head(5)
solutionprobabilitycost
415{‘x’: [[0, 1, 0, 0], [1, 0, 0, 1], [0, 0, 1, 0…0.0004883.4
1022{‘x’: [[0, 1, 0, 1], [0, 0, 0, 0], [1, 0, 1, 0…0.0004884.5
348{‘x’: [[1, 0, 0, 0], [0, 0, 0, 1], [0, 0, 1, 0…0.00048821.9
169{‘x’: [[0, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 1…0.00097723.3
1215{‘x’: [[0, 0, 0, 1], [1, 0, 0, 0], [0, 0, 1, 0…0.00048823.5
Compare the optimized results to uniformly sampled results: 
uniform_result = combi.sample_uniform()
And compare the histograms: 
optimization_result["cost"].plot(
    kind="hist",
    bins=50,
    edgecolor="black",
    weights=optimization_result["probability"],
    alpha=0.6,
    label="optimized",
)
uniform_result["cost"].plot(
    kind="hist",
    bins=50,
    edgecolor="black",
    weights=uniform_result["probability"],
    alpha=0.6,
    label="uniform",
)
plt.legend()
plt.ylabel("Probability", fontsize=16)
plt.xlabel("cost", fontsize=16)
plt.tick_params(axis="both", labelsize=14)
output

​
Best Solution

# 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=(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["x"][i] for i in range(len(best_solution.solution["x"]))],
    best_solution.cost,
    is_classic=False,
)
Output:
<IPython.core.display.

Markdown object>
Output:
total cost=  3.4
B1B2B3B4
A10100
A21001
A30010
output

​
Comparing to a Classical Solver

from pyomo.opt import SolverFactory

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

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,
)
Output:
<IPython.core.display.

Markdown object>
Output:
total cost=  2.4999999996418167
B1B2B3B4
A11.0000000.0000000.0000000.000000
A20.0000001.0000000.0000001.000000
A30.000000-0.0000001.0000000.000000
output

​
Reference

[1] O. V. Shemelova, E. V. Yakovleva, T. G. Makuseva, I. I. Eremina, and O. N. Makusev. (2019). Solving optimization problems when designing power supply circuits. E3S Web of Conferences 124, 04011.