Skip to content

Facility Location Problem (P-median)

We consider the following optimization problem: we have a set of \(M\) customers and a set of \(N\) potential locations for opening a facility. Given transportation costs between facilities and customers and how many facilities we would like to open (\(P\)), we need to determine which facilities to open such that the total transportation cost between facilities and customers is minimal, under the constraint that each customer is allocated to only one facility.

Possible extensions:

\(\bullet\) Add different demand for each customer.

\(\bullet\) Add cost for opening a facility at each location.

\(\bullet\) Different categories of facilities, and customer can have multiple allocations for different facilities.

Mathematical Modeling

The input of the model is a set of \(M\) customers \(\{1,\dots,M\}\), a set of \(N\) potential locations for facilities \(\{1,\dots,N\}\), an \(N\times M\) matrix \(d\) where \(d_{nm}\) is the cost of customer \(m\) buying from facility \(n\), and the total number of facilities we want to open \(P\).

We define a binary variable for the optimization problem: an \(N\times M\) matrix \(x\) such that \(x_{nm}=1\) if customer \(m\) is allocated to facility \(n\).

The objective function to minimize is the total cost function: $$ \min_{x} \sum_{n,m} d_{nm}x_{nm} $$ We have 2 constraints:

(1) Each customer is supplied: \(\forall m\in[0,M] \,\,\, \sum_n x_{nm}=1\).

(2) Total number of opened facilities is \(P\): \(\sum_n\Pi_m (1-x_{nm})=N-P\) (The inner product is zero if the \(n-\)th facility is not open).

Alternative modeling

We note that there is an alternative modeling, in which one adds another variable to the model: a binary vector \(y\) of size \(N\) which indicates which facilities are open. In this formulation the second constraint can be written as \(\sum_n y_n=P\) together with an inequality constraint \(\forall n,m:\, x_{nm} \leq y_n\). A model which combines equality and inequality constraints will be available in the future.

Note that in this alternative modeling we have a QUBO problem (compared to the formulation above where constraint (2) is polynomial of degree \(m\)). However, the alternative modelling has more variables to minimize on, and thus refers to more qubits.

An Example

Let us say we can open facilities in Japan, USA, and France, and we have 4 customers whose costs for buying from these three locations are given. We would like to open in total \(P=2\) facilities. The optimization problem is to find where to open the facilities and which customer is allocated to which facility.

We now draw this specific example on a graph. We have \(N=3\) locations and \(M=4\) customers, where the weights of the edges between them signify the costs:

Suggestions: - Give general problem description and then givens, or provide givens after each mention. - Standardize writing numbers appearing in sentences out when less than 10. "4" vs "four".

# Import relevant packages

from itertools import product

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

# Declare givens from problem statement

Facilities = ["Japan", "USA", "France"]
Customers = ["A", "B", "C", "D"]

N = len(Facilities)  # potential facility count
M = len(Customers)  # customer count
P = 2  # allocated facility count

# costs of customers using facilities
d = np.array(
    [[0.02, 0.14, 0.62, 0.11], [0.99, 0.22, 0.91, 0.09], [0.4, 0.76, 0.95, 0.61]]
)


graph = nx.DiGraph()
graph.add_nodes_from(Facilities + Customers)
for n, m in product(range(N), range(M)):
    graph.add_edges_from([(Facilities[n], Customers[m])], weight=d[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=Customers, font_size=22, font_color="None")
nx.draw_networkx_nodes(
    graph, pos, nodelist=Customers, node_color="#119DA4", node_size=500
)
for fa in Facilities:
    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 Customers},
    font_size=22,
    font_color="#F4F9E9",
)

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

png

Building the Pyomo model from a matrix of distances input

from typing import List, Tuple, cast  # noqa

import pyomo.environ as pyo
from IPython.display import Markdown, display
## We define a function which gets a matrix of costs between customers and potential facilities
## and the number of facilities to open.


def pmedian(cost_mat: np.ndarray, P: int) -> pyo.ConcreteModel:
    model = pyo.ConcreteModel("pmedian")

    N = cost_mat.shape[0]  # potential facility amount
    M = cost_mat.shape[1]  # customer count
    Locations = range(N)
    Customers = range(M)

    model.x = pyo.Var(Locations, Customers, domain=pyo.Binary)

    @model.Constraint(Customers)
    def each_customer_is_supplied_rule(model, m):  # constraint (1)
        return sum(model.x[n, m] for n in Locations) == 1

    def is_location_alocated(n):  # constraint (2)
        return np.prod([(1 - model.x[n, m]) for m in Customers])

    model.num_facilities = pyo.Constraint(
        expr=sum(is_location_alocated(n) for n in Locations) == N - P
    )

    model.cost = pyo.Objective(
        expr=sum(cost_mat[n, m] * model.x[n, m] for n in Locations for m in Customers),
        sense=pyo.minimize,
    )

    return model

Solving with Classiq

We take a specific example: the one outlined above

pmedian_model = pmedian(d, P)
from classiq import construct_combinatorial_optimization_model
from classiq.applications.combinatorial_optimization import OptimizerConfig, QAOAConfig

qaoa_config = QAOAConfig(num_layers=5, penalty_energy=10)

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

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

Lastly, we define the CombinatorialOptimization instance which we can use to solve the problem:

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

We also set the quantum backend we want to execute on:

from classiq import set_execution_preferences
from classiq.execution import ClassiqBackendPreferences, ExecutionPreferences

backend_preferences = ExecutionPreferences(
    backend_preferences=ClassiqBackendPreferences(backend_name="aer_simulator")
)

qmod = set_execution_preferences(qmod, backend_preferences)
from classiq import write_qmod

write_qmod(qmod, "facility_location")

We define a function that plot solutions

# 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({Facilities[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 = Customers
    df.index = Facilities
    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(Facilities + Customers)
    for n, m in product(range(N), range(M)):
        if x_sol_to_mat[n, m] > 0:
            graph_sol.add_edges_from([(Facilities[n], Customers[m])], weight=d[n, m])

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

    nx.draw_networkx(
        graph_sol, pos=pos, nodelist=Customers, font_size=22, font_color="None"
    )
    nx.draw_networkx_nodes(
        graph_sol, pos, nodelist=Customers, node_color="#119DA4", node_size=500
    )
    for fa in Facilities:
        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 Customers},
        font_size=22,
        font_color="#F4F9E9",
    )

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

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

from classiq import show, synthesize

qprog = synthesize(qmod)
show(qprog)
Opening: https://platform.classiq.io/circuit/d804bf78-fab4-4ab0-b897-5fddca1b51f2?version=0.41.0.dev39%2B79c8fd0855

We now solve the problem 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

png

Optimizer statistics

import pandas as pd

from classiq.applications.combinatorial_optimization import (
    get_optimization_solution_from_pyo,
)

solution = get_optimization_solution_from_pyo(
    pmedian_model, vqe_result=vqe_result, penalty_energy=qaoa_config.penalty_energy
)
optimization_result = pd.DataFrame.from_records(solution)
optimization_result.sort_values(by="cost", ascending=True).head(5)
probability cost solution count
535 0.001 0.95 [1, 0, 1, 0, 0, 1, 0, 1, 0, 0, 0, 0] 1
504 0.001 0.97 [1, 0, 1, 1, 0, 1, 0, 0, 0, 0, 0, 0] 1
257 0.001 1.27 [0, 1, 1, 1, 0, 0, 0, 0, 1, 0, 0, 0] 1
536 0.001 1.39 [1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1] 1
35 0.004 1.51 [1, 0, 1, 1, 0, 0, 0, 0, 0, 1, 0, 0] 4
optimization_result.hist("cost", weights=optimization_result["probability"])
array([[<Axes: title={'center': 'cost'}>]], dtype=object)

png

Best solution

We plot the quantum result only if we get the right solution (to avoid problems with printing table and graph).

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

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

QAOA SOLUTION

total cost=  0.9499999999999922
  A B C D
Japan 1 0 1 0
USA 0 1 0 1
France 0 0 0 0

png

Compare to a classical solver

from pyomo.opt import SolverFactory

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

pmedian_model.display()
Model pmedian

  Variables:
    x : Size=12, Index=x_index
        Key    : Lower : Value                  : Upper : Fixed : Stale : Domain
        (0, 0) :     0 :                    1.0 :     1 : False : False : Binary
        (0, 1) :     0 :                    1.0 :     1 : False : False : Binary
        (0, 2) :     0 :                    1.0 :     1 : False : False : Binary
        (0, 3) :     0 :                    0.0 :     1 : False : False : Binary
        (1, 0) :     0 :                    0.0 :     1 : False : False : Binary
        (1, 1) :     0 :  6.846989283059948e-10 :     1 : False : False : Binary
        (1, 2) :     0 :                    0.0 :     1 : False : False : Binary
        (1, 3) :     0 :                    1.0 :     1 : False : False : Binary
        (2, 0) :     0 :                    0.0 :     1 : False : False : Binary
        (2, 1) :     0 : -6.846989283059948e-10 :     1 : False : False : Binary
        (2, 2) :     0 :                    0.0 :     1 : False : False : Binary
        (2, 3) :     0 :                    0.0 :     1 : False : False : Binary

  Objectives:
    cost : Size=1, Index=None, Active=True
        Key  : Active : Value
        None :   True : 0.8699999996302625

  Constraints:
    each_customer_is_supplied_rule : Size=4
        Key : Lower : Body : Upper
          0 :   1.0 :  1.0 :   1.0
          1 :   1.0 :  1.0 :   1.0
          2 :   1.0 :  1.0 :   1.0
          3 :   1.0 :  1.0 :   1.0
    num_facilities : Size=1
        Key  : Lower : Body              : Upper
        None :   1.0 : 1.000000000684699 :   1.0
best_classical_solution = np.array(
    [pyo.value(pmedian_model.x[idx]) for idx in np.ndindex(d.shape)]
).reshape(d.shape)

plotting_sol(best_classical_solution, pyo.value(pmedian_model.cost), is_classic=True)

CLASSICAL SOLUTION

total cost=  0.8699999996302625
  A B C D
Japan 1.000000 1.000000 1.000000 0.000000
USA 0.000000 0.000000 0.000000 1.000000
France 0.000000 -0.000000 0.000000 0.000000

png