Intelligent Scheduling, Efficient Empowerment: Unveiling the Path to Optimal Resource Allocation in Computing Power Networks

A computing power network, in simple terms, is a network that delivers computational power just like electricity or water. Imagine that the phones, computers, and other devices we use daily all require processors to compute and run various programs. The computing power network is like a massive ‘processor warehouse’ that consolidates scattered computing resources and optimizes their distribution via the internet to meet the computational needs of different users and scenarios.

For example, when we watch high-definition videos, play large-scale games, or perform complex data analysis, these operations require substantial computing power. The computing power network acts like an intelligent ‘power system’ that can quickly provide the necessary computational resources, ensuring a smooth and efficient user experience. This way, we can enjoy powerful computing capabilities anytime and anywhere without worrying about insufficient device performance.

Challenges Faced by Computing Power Networks

• Resource Allocation Efficiency: A computing power network needs to efficiently allocate computational resources to meet the diverse needs of different users and applications. This requires the network to monitor resource usage in real-time and dynamically adjust resource allocation strategies to avoid resource idleness or overload.

• Network Latency: Quick response to computational tasks requires a low-latency network environment. A computing power network must address the latency issues in data transmission, especially when handling applications with high real-time requirements, such as autonomous driving.

• Low Cost: Building and maintaining a computing power network requires significant financial investment, so it is necessary to minimize costs while ensuring performance.

Components of a Computing Power Network

A computing power network mainly consists of three parts:

• Edge Devices: These are various smart devices that are responsible for collecting data.

• Edge Servers: Located near the devices, they are responsible for processing part of the data and alleviating the burden on cloud servers.

• Cloud Servers: Located in data centers, they have powerful computing capabilities and are responsible for processing more complex data.

Problem Description

Consider an optimization problem for the computing power network layout in a specific area. The area is divided into several adjacent square grids, and the computing power demand distribution data provides the computing power demand within each grid. The coordinate values in the data represent the center coordinates of the grid. To simplify the problem, the computing power demand points within each grid are unified to the grid center point (i.e., one grid corresponds to one demand point).

The computing power demand within the grid is generated by edge devices, which are terminal devices connected to the network, such as sensors, smartphones, and industrial robots. The computing power demand in the computing power network is met by edge servers and cloud servers. Edge servers are located at the ‘edge’ of the network, usually close to the end users or devices. Their task is to process data closer to the user to improve response speed and efficiency. Edge servers can process requests faster because they are closer to the users. Edge servers can also alleviate the burden on the core cloud infrastructure, improving overall operational efficiency. Cloud servers are located in data centers far from users and have powerful computing and storage capabilities. When the capacity of edge servers is insufficient, cloud servers can serve as a supplement. The collaboration between edge servers and cloud servers can optimize the performance and reliability of the entire system.

Problem 1: In Problem 1, we only consider meeting the demand within the computing region using edge servers. Suppose we want to set up 2 edge servers within the grid area with computing power demand distribution, with each edge server having a coverage radius of 1. A QUBO model is required to determine at which points to place the edge servers to cover the maximum computing power demand.

Question 2: When the edge servers cannot meet the computational demand, upstream cloud servers will provide the necessary computing services. Now, a cloud server is added outside the grid area, and both end nodes and edge servers can choose to connect to the cloud node. When the computational demand received by the edge server exceeds its capacity, the excess demand will be directly allocated to the cloud server. The computational demand of each end node must be met and can only be served by one server, which can be either a cloud node or an edge server. Since edge servers have a resource capacity limit, assume that the available computational resource capacity of each edge server is 12, and the available computational resource capacity of the cloud server is infinite, meaning that the cloud server’s available computational resource capacity is not considered. Servers have a certain coverage radius, with the edge server’s coverage radius assumed to be 3, and the cloud server’s coverage radius is infinite, meaning that the cloud server’s coverage range is not considered.

Establishing edge servers typically incurs costs, which consist of fixed costs, computation costs, and transmission costs. The fixed cost is related to whether and where the edge server is established. The computation cost is proportional to the amount of requested computational resources, calculated as the unit computation cost multiplied by the computation amount. The unit computation cost for cloud servers is 1, while for edge servers, it is 2. Additionally, there is a transmission cost between the end node and the edge, from the edge to the cloud, and from the end to the cloud. The transmission cost is calculated by multiplying the computational demand by the transmission distance and the unit transmission cost. The transmission distance is calculated using the Euclidean distance, rounded to two decimal places, and calculated as a one-way distance (round-trip transmission is not considered). The unit transmission cost from the end to the edge and from the edge to the cloud is 1, while the unit transmission cost from the end to the cloud is 2.

To meet all the end-side computational demands within the area, establish a QUBO model to solve for the network layout that minimizes the overall cost, including the location and number of edge servers, and the connections between end-to-edge, edge-to-cloud, and end-to-cloud nodes.

The essence of the problem

In the problem of computational network layout optimization, on the surface, we may think of classical resource scheduling or network flow optimization problems and try to solve ‘how to allocate resources among different areas to minimize costs’ using traditional greedy algorithms or dynamic programming methods. These methods might help us find a local optimal solution at a specific moment, but the complexity of this problem goes far beyond the surface. We not only need to consider the resource allocation of each computational node but also globally optimize the entire network layout, ensuring that each resource is reasonably utilized and that efficient computation and communication are maintained across multiple regions.

In-depth analysis

The core of the computational network layout optimization problem is how to find an optimal resource allocation and node layout scheme across the entire network, similar to finding the optimal computational node layout and communication path combination in graph G. This is essentially a combinatorial optimization problem, requiring the selection of the optimal solution from all possible network layouts, and the number of these possible layouts is enormous. For example, if we have n potential edge server deployment locations and m computational task points, the combination of these locations and tasks will grow exponentially, with n^m possible layout combinations.

Obviously, a simple brute-force method is unrealistic because even for a moderately sized network, the number of combinations is astronomical, far beyond the capabilities of classical computers. Just like the path selection problem in the traveling salesman problem, we face a combinatorial challenge here, where every node and resource configuration in the network could affect the overall optimization result.

In other words, as the problem size increases, the time complexity of the algorithm grows non-polynomially, making it difficult to solve within a reasonable time. Therefore, the key to the problem is not how to locally optimize a single resource allocation, but how to globally optimize the entire computational network layout, finding the optimal layout path that minimizes overall computation costs. The complexity of this problem lies not only in the increase in computation but also in the need for multi-dimensional optimization to find a global optimal solution.

Reference model for the first question

Restatement of the first question

In a 4x4 grid, each grid represents an area with computational demand. Our goal is to determine in which grids to deploy edge computing nodes to maximize the covered user demand.

Symbol Definitions

• $R=1$: Coverage radius of edge computing nodes

• $P=2$: Number of planned edge computing nodes to be deployed.

• $\mathcal{L}_{u}$: Set of user locations.

• $\mathcal{L}_{e}$: Set of candidate edge node locations.

• $\text{demand}_{l} (l\in \mathcal{L}_{u})$: Computational demand at grid $l$.

• $d_{ij}$: Distance between grid $i$ and grid $j$.

• $\alpha_{ij}=\mathbb{1}_{d_{ij}\leq R}$: Indicates whether the distance between grid $i$ and grid $j$ does not exceed the coverage radius $R$ of the edge node.

Decision Variables

• $x_{i}$: A binary variable indicating whether to deploy an edge computing node in grid $i$.

• $z_{i}$: A binary variable indicating whether grid $i$ is covered.

Objective Function

Maximize the total covered computational demand:

$$\max \sum_{i\in \mathcal{L}_{u}} {z_{i}\cdot \text{demand}_{i}}$$

Constraints

• The coverage status of each user grid does not exceed the coverage status of the surrounding edge nodes:

  $$z_{i}\leq \sum_{j \in \mathcal{L}_{e}}{\alpha_{ij}\cdot x_{j}}, \forall i\in \mathcal{L}_{u}$$

• The total number of deployed edge computing nodes is equal to $P$:

  $$\sum_{j\in \mathcal{L}_{e}} x_{j} = P$$

Simplified Model (Inclusion-Exclusion Principle)

When $P=2$, we can simplify the model using the inclusion-exclusion principle:

$$\begin{aligned} \max \quad & \sum_{i\in \mathcal{L}_{u}} \left( \sum_{j\in \mathcal{L}_{e}} \alpha_{ij} x_j - \sum_{j,k\in \mathcal{L}_{e}} \alpha_{ij}\alpha_{ik} x_j x_k \right) \mathrm{demand}_i \\[0.5ex] \text{s.t.} \quad & \sum_{j\in \mathcal{L}_{e}} x_j = 2 \end{aligned}$$

Reference Model for Question 2

Problem Overview

• A 6×6 grid where each grid has a certain computational demand (only 11 grids have non-zero demand), and all computational demands need to be met.

• There are 5 candidate locations for edge servers, each with limited computational capacity. When the computational requests received by an edge server exceed its capacity, the excess requests are sent to the cloud server.

• One location hosts a cloud server, which has no capacity limit.

• Users can connect to either the edge server or the cloud server (but only one); if the edge server’s capacity is insufficient, it can connect to the cloud server.

• Minimize cost: fixed cost + variable cost (computational cost) + transmission cost (= unit cost * distance * transmission volume)

Symbol Definitions

• $C_{\text{edge}}$: Edge server capacity.

• $\text{c-fix-edge}_{j}$: Fixed cost of establishing an edge computing node in grid $j$.

• $\text{c-var-cloud}$: Unit computational cost at the cloud node.

• $\text{c-var-edge}$: Unit computational cost at the edge node.

• $\text{c-tran}_{u,e}$: Unit transmission cost from user side to edge side.

• $\text{c-tran}_{u,c}$: Unit transmission cost from user side to cloud side.

• $\text{c-tran}_{e,c}$: Unit transmission cost from edge side to cloud side.

• $d_{ij}$: Distance between grid $i$ and grid $j$.

• $d_{i}$: The distance between grid $i$ and the cloud server.

• $\alpha_{ij}=\mathbb{1}_{d_{ij} \leq R}$: Whether the distance between grid $i$ and grid $j$ does not exceed the edge node coverage radius $R$

Intermediate Variables

• $y^{u,e}_{ij}$: Whether the demand of grid $i$ is served by the edge server located at grid $j$

• $y^{e,c}_{j}$: Whether the edge server at grid $j$ is connected to the cloud server.

• $y^{u,c}_{i}$: Whether the demand of grid $i$ is served by the cloud server.

• $u_{j} \in \mathbb{N}$: The number of computational demands exceeding the capacity of edge server $j$ that are served by the cloud server.

Decision Variables

• $x_{\text{edge}, j}$: Whether to open an edge server at location $j$.

Mathematical Model

Objective Function:

$$\min \text{Cost}_{\text{fix}} + \text{Cost}_{\text{var}} + \text{Cost}_{\text{tran}}$$

where

$$\begin{aligned} \text{Cost}_{\text{fix}} = \sum_{j \in \mathcal{L}_{e}} \text{c-fix-edge}_{j}\, x_{\text{edge},j} \end{aligned}$$

$$\begin{aligned} \text{Cost}_{\text{var}} =\;& \text{c-var-cloud} \sum_{i \in \mathcal{L}_{u}} \text{demand}_{i}\, y^{u,c}_{i} \\ &+ \text{c-var-edge} \sum_{j \in \mathcal{L}_{e}} \sum_{i \in \mathcal{L}_{u}} \text{demand}_{i}\, y^{u,e}_{ij} \\ &+ (\text{c-var-cloud}-\text{c-var-edge}) \sum_{j \in \mathcal{L}_{e}} u_j \end{aligned}$$

$$\begin{aligned} \text{Cost}_{\text{tran}} =\;& \text{c-tran}_{u,e} \sum_{i \in \mathcal{L}_{u}} \sum_{j \in \mathcal{L}_{e}} y^{u,e}_{ij}\, \text{demand}_{i}\, d_{ij} \\ &+ \text{c-tran}_{u,c} \sum_{i \in \mathcal{L}_{u}} y^{u,c}_{i}\, \text{demand}_{i}\, d_i \\ &+ \text{c-tran}_{e,c} \sum_{j \in \mathcal{L}_{e}} d_j\, u_j \end{aligned}$$

In this model, we express $u_{j}$ as

$$u_{j} = y^{e,c}_{j} \cdot \Big(\sum_{i \in \mathcal{L}_{u}} \text{demand}_{i} \cdot y^{u,e}_{ij} - C_{\text{edge}}\Big)$$

At the same time, we add constraints so that $y^{u,e}_{jk}$ will take the value of $1$ only when the computational demand received by edge server $j$ exceeds its capacity limit, otherwise it will be $0$. Although the expression of $u_{j}$ contains quadratic terms, in this model, the cloud server has no capacity limit, and $u_{j}$ only appears in the objective function, not in the constraints, ensuring that no higher-order terms appear when converted to the QUBO model.

Constraints:

• The computational demand points are allocated to either edge or cloud and are served by only one.

  $$\sum_{j \in \mathcal{L}_{e}}{y^{u,e}_{ij}} + y^{u,c}_{i} = 1, \forall i \in \mathcal{L}_{u}$$

• Coverage relationship, and only if the edge server is opened can it connect from the demand point.

  $$y^{u,e}_{ij} \leq \alpha_{ij} \cdot x_{\text{edge}, j}, \forall i \in \mathcal{L}_{u}, j \in \mathcal{L}_{e}$$

• The relationship between edge connection to the cloud and opening the edge server.

  $$y^{e,c}_{j} \leq x_{\text{edge}, j}, \forall j \in \mathcal{L}_{e}$$

• Only when the computational demand received by edge server $j$ exceeds its capacity limit, will $y^{u,e}_{jk}$ take the value of $1$.

  $$y^{e,c}_{j} \cdot C_{\text{edge}} \leq \sum_{i \in \mathcal{L}_{u}}{\text{demand}_{i} \cdot \alpha_{ij} \cdot y^{u,e}_{ij}} \quad \text{(1)}$$
  $$\sum_{i \in \mathcal{L}_{u}}{\text{demand}_{i} \cdot \alpha_{ij} \cdot y^{u,e}_{ij}} - C_{\text{edge}} \leq (\text{max-uj}[j] - C_{\text{edge}}) \cdot y^{e,c}_{j} \quad \text{(2)}$$

  where :math:text{max-uj}[j] is the upper bound of the computational demand received by edge server $j$, which we set as

  $$\text{max-uj}[j] = \sum_{i \in \mathcal{L}_{u}}{\text{demand}_{i} \cdot \alpha_{ij}}$$

Preprocessing

For candidate locations of edge servers $j$, if the computational demand they receive is always less than the edge server capacity limit, for example,

$$\text{max-uj}[j] \leq C_{\text{edge}}$$

When this is true, we can set $y^{e,c}_{j} = 0$, and for this candidate location, inequalities (1) and (2) can be ignored to reduce the number of bits (relaxation variables).

Code for Question 1

import math
import numpy as np
import kaiwu as kw


class EdgeCoveragePlanner:
    """Solver for Edge Coverage Planning problem."""

    def __init__(
        self,
        num_nodes: int = 2,
        coverage_range: float = 2.0,
        penalty: float = 100.0,
        num_slack_bins: int = 1,
    ) -> None:
        """Initialize EdgeCoveragePlanner parameters.

        Args:
            num_nodes: Number of edge nodes to be selected.

            coverage_range: Node coverage threshold.

            penalty: penalty coefficient.

            num_slack_bins: Slack binary digits, used to override constraints.
        """
        self.num_nodes = num_nodes
        self.coverage_range = coverage_range
        self.penalty = penalty
        self.num_slack_bins = num_slack_bins

        # Data containers
        self.demand = {}
        self.locations = []
        self.distances = {}
        self.coverage = {}
        self.dimension = 0

        # QUBO model
        self.model = None

    def prepare_data(self, demand_data):
        """prepare demand, distance, and coverage matrix data.

        Args:
            demand_data: dict, the key is the coordinate string 'i, j',
                        and the value corresponds to the calculation requirement.
        """
        self.demand = demand_data
        # Create a list containing all location coordinates
        self.locations = list(demand_data.keys())
        # Calculate the number of location coordinates
        self.dimension = int(math.sqrt(len(self.demand)))
        # Iterate through all location coordinates to calculate distances between pairs
        self.distances = {
            (i, j): np.linalg.norm(
                np.array([int(coord) for coord in i.split(",")])
                - np.array([int(coord) for coord in j.split(",")])
            )
            for i in self.locations
            for j in self.locations
        }
        # Initialize a dictionary to store the coverage relationship between locations
        self.coverage = {
            (i, j): int(dist <= self.coverage_range)
            for (i, j), dist in self.distances.items()
        }

    def prepare_model(self):
        """Building a Qubo Model"""
        # Create binary variable arrays x and z to represent edge computing node locations
        # and demand coverage conditions
        var_x = kw.core.ndarray((self.dimension, self.dimension), "x", kw.core.Binary)
        var_z = kw.core.ndarray((self.dimension, self.dimension), "z", kw.core.Binary)
        limit_z = kw.core.zeros((self.dimension, self.dimension))
        # Initialize slack variables
        slack = kw.core.ndarray(
            (self.dimension, self.dimension),
            "slack",
            kw.core.Integer,
            (0, (2**self.num_slack_bins) - 1),
        )
        # Define the objective function to minimize the total computing power demand
        obj = kw.core.quicksum(
            self.demand[f"{i + 1},{j + 1}"] * var_z[i, j]
            for i in range(self.dimension)
            for j in range(self.dimension)
        )
        # Iterate through all location coordinates to build the second constraint
        for i in range(self.dimension):
            for j in range(self.dimension):
                # For each location, the demand coverage z[i] should be less than or equal to
                # the total coverage provided by all edge computing nodes
                limit_z[i, j] = kw.core.quicksum(
                    self.coverage[f"{i + 1},{j + 1}", f"{i1 + 1},{j1 + 1}"]
                    * var_x[i1, j1]
                    for i1 in range(self.dimension)
                    for j1 in range(self.dimension)
                )

        self.model = kw.core.QuboModel()
        self.model.set_objective(-obj)
        self.model.add_constraint(
            var_x.sum() == self.num_nodes, "c1", penalty=self.penalty
        )
        self.model.add_constraint(
            var_z <= limit_z, "c2", penalty=self.penalty, slack_var_expr=slack
        )

    def solve(self):
        """Solving the QUBO model.

        Returns:
            tuple: Result dictionary and Result dictionary.

            - dict: Result dictionary. The key is the variable name, and the value is the corresponding spin value.

            - float: qubo value.
        """
        # Perform the Simulated Annealing algorithm
        _solver = kw.classical.SimulatedAnnealingOptimizer(
            initial_temperature=100000,
            alpha=0.99,
            cutoff_temperature=0.0001,
            iterations_per_t=100,
            rand_seed=10,
        )
        _sol_dict, _qubo_value = _solver.solve_qubo(self.model)
        return _sol_dict, float(_qubo_value)

    def recovery(self, sol_dict):
        """Verify whether the solution is feasible"""
        return self.model.verify_constraint(sol_dict)


def assert_documented_results():
    """Check that the documented model-building path stays valid."""
    demand_data = {
        "1,1": 38,
        "1,2": 22,
        "1,3": 65,
        "1,4": 56,
        "2,1": 53,
        "2,2": 48,
        "2,3": 76,
        "2,4": 46,
        "3,1": 56,
        "3,2": 36,
        "3,3": 7,
        "3,4": 29,
        "4,1": 50,
        "4,2": 37,
        "4,3": 48,
        "4,4": 40,
    }
    solver = EdgeCoveragePlanner()
    solver.prepare_data(demand_data)
    solver.prepare_model()
    assert solver.model is not None


if __name__ == "__main__":
    # Store computing power demand data in the DEM dictionary,
    # where the keys are location coordinates and the values are computing power demands
    # fmt: off
    demand_data = {
        '1,1': 38, '1,2': 22, '1,3': 65, '1,4': 56,
        '2,1': 53, '2,2': 48, '2,3': 76, '2,4': 46,
        '3,1': 56, '3,2': 36, '3,3': 7,  '3,4': 29,
        '4,1': 50, '4,2': 37, '4,3': 48, '4,4': 40
    }
    # fmt: on
    # Create an instance of the SPQCSolver class
    solver = EdgeCoveragePlanner()

    # Prepare data
    solver.prepare_data(demand_data)

    # Prepare the QUBO model with a specified penalty coefficient lambda
    solver.prepare_model()

    # Use the Simulated Annealing algorithm to find the optimal solution
    best_sol_dict, qubo_value = solver.solve()

    # Recover the original problem solution from the QUBO solution and check its feasibility
    unsatisfied_count, result_dict = solver.recovery(best_sol_dict)
    if unsatisfied_count == 0:
        print("Find a feasible solution")
        print("Objective value:", -qubo_value)
    else:
        print("No feasible solution")

Code for Question 2

import math
from typing import Tuple
import numpy as np
import kaiwu as kw
from kaiwu.core import quicksum


class CloudEdgeUserSolver:
    """Solver for cloud-edge-user cost-minimization QUBO model."""

    def __init__(self, edge_capacity: int = 12, edge_radius: float = 3.0) -> None:
        """Initialize core parameters and placeholders."""
        # capacities and radii
        self.edge_capacity = edge_capacity
        self.edge_radius = edge_radius

        # Service node configurations
        self.user_locations = [
            "1,1",
            "1,4",
            "2,6",
            "3,3",
            "3,5",
            "4,4",
            "5,5",
            "6,3",
            "6,6",
        ]
        self.edge_locations = ["6,1", "2,3", "4,5", "6,5"]
        self.cloud_locations = ["4,0"]

        # edge server costs
        self.fix_cost_edge = {"6,1": 60, "2,3": 42, "4,5": 46, "6,5": 54}

        # variable costs
        self.var_cost_cloud = 1
        self.var_cost_edge = 2

        # Unit transmission cost
        self.trans_cost_user_edge = 1
        self.trans_cost_user_cloud = 2
        self.trans_cost_edge_cloud = 1

        # Demand data
        self.demand = {
            "1,1": 7,
            "1,4": 4,
            "2,6": 5,
            "3,3": 9,
            "3,5": 8,
            "4,4": 11,
            "5,5": 1,
            "6,3": 7,
            "6,6": 5,
        }

        # Precomputed data
        self.num_user = 0
        self.num_edge = 0
        self.num_cloud = 0
        self.distances = {}
        self.coverage = {}
        self.qubo_model = None

    def prepare_data(self):
        """Prepare data for QUBO model."""
        # initialize user demands
        # fmt: off
        loc_set_inner_full = ['1,1', '1,2', '1,3', '1,4', '1,5', '1,6',
                              '2,1', '2,2', '2,3', '2,4', '2,5', '2,6',
                              '3,1', '3,2', '3,3', '3,4', '3,5',
                              '4,1', '4,2', '4,3', '4,4', '4,5', '4,6',
                              '5,1', '5,2', '5,3', '5,4', '5,5', '5,6',
                              '6,1', '6,2', '6,3', '6,4', '6,5', '6,6']
        # fmt: on
        # all nodes
        loc_set = loc_set_inner_full + self.cloud_locations

        self.num_cloud = len(self.cloud_locations)
        self.num_edge = len(self.edge_locations)
        self.num_user = len(self.user_locations)

        # compute distances
        for i in loc_set:
            for j in loc_set:
                self.distances[(i, j)] = round(
                    np.sqrt(
                        (
                            int(i.split(",", maxsplit=1)[0])
                            - int(j.split(",", maxsplit=1)[0])
                        )
                        ** 2
                        + (
                            int(i.split(",", maxsplit=1)[1])
                            - int(j.split(",", maxsplit=1)[1])
                        )
                        ** 2
                    ),
                    2,
                )

        # compute coverage relationships
        for i in loc_set:
            for j in loc_set:
                if self.distances[(i, j)] <= self.edge_radius:
                    self.coverage[(i, j)] = 1
                else:
                    self.coverage[(i, j)] = 0

    def prepare_model(
        self,
        eq_penalties: Tuple[float, float, float] = (1e4, 1e4, 1e4),
        ineq_penalties: Tuple[float, float] = (1e4, 1e4),
    ) -> None:
        """
        Build QUBO with 3 equality and 2 inequality constraints.
        """
        # compute the maximum demand for each edge
        max_ujk = []
        for j in self.edge_locations:
            max_ujk.append(
                sum(
                    self.demand[i] * self.coverage[i, j] for i in self.user_locations
                )
            )

        # initialize decision variables
        x_edge = kw.core.ndarray(self.num_edge, "x_edge", kw.core.Binary)
        y_ij = kw.core.ndarray((self.num_user, self.num_edge), "yij", kw.core.Binary)
        y_jk = kw.core.ndarray((self.num_edge, self.num_cloud), "yjk", kw.core.Binary)
        y_ik = kw.core.ndarray((self.num_user, self.num_cloud), "yik", kw.core.Binary)

        # for each edge, if the maximum demand does not exceed the edge's capacity, set the corresponding yjk to 0
        for j in range(self.num_edge):
            if max_ujk[j] <= self.edge_capacity:
                y_jk[j][0] = 0

        # initialize ujk related variables
        u_jk = np.zeros(
            shape=(self.num_edge, self.num_cloud), dtype=kw.core.BinaryExpression
        )
        ujk_residual = np.zeros(shape=self.num_edge, dtype=kw.core.BinaryExpression)
        for j in range(self.num_edge):
            ujk_residual[j] = (
                quicksum(
                    self.demand[self.user_locations[i]] * y_ij[i][j]
                    for i in range(self.num_user)
                )
                - self.edge_capacity
            )
            for k in range(self.num_cloud):
                u_jk[j][k] = y_jk[j][k] * ujk_residual[j]
        # build objective
        obj = self._build_objective_components(u_jk, x_edge, y_ij, y_ik)
        # build equality constraint
        constraint1, constraint2, constraint3 = self._build_eq_constraint(
            x_edge, y_ij, y_ik, y_jk
        )
        # build inequality constraint
        ineq_qubo1, ineq_qubo2 = self._build_ineq_constraint(
            max_ujk, ujk_residual, y_ij, y_jk
        )

        # building the final model
        self.qubo_model = kw.core.QuboModel()
        self.qubo_model.set_objective(obj)
        self.qubo_model.add_constraint(constraint1, name="c1", penalty=eq_penalties[0])
        self.qubo_model.add_constraint(constraint2, name="c2", penalty=eq_penalties[1])
        self.qubo_model.add_constraint(constraint3, name="c3", penalty=eq_penalties[2])
        self.qubo_model.add_constraint(ineq_qubo1, name="c4", penalty=ineq_penalties[0])
        self.qubo_model.add_constraint(ineq_qubo2, name="c5", penalty=ineq_penalties[1])

    def _build_ineq_constraint(self, max_ujk, ujk_residual, y_ij, y_jk):
        # inequality constraint 1: after subtracting the edge's maximum capacity from demand,
        # the yjk constraint should hold
        ineq_constraint1 = []
        ineq_qubo1 = kw.core.BinaryExpression(coefficient={}, offset=0)
        len_slack1 = math.ceil(math.log2(max(max_ujk) + 1))
        slack1 = kw.core.ndarray(
            (self.num_edge, self.num_cloud, len_slack1), "slack1", kw.core.Binary
        )
        for j in range(self.num_edge):
            ineq_constraint1.append([])
            for k in range(self.num_cloud):
                if y_jk[j][k] == 0:
                    ineq_constraint1[j].append(0)
                else:
                    ineq_constraint1[j].append(
                        ujk_residual[j] - (max_ujk[j] - self.edge_capacity) * y_jk[j][k]
                    )
                    ineq_qubo1 += (
                        ineq_constraint1[j][k]
                        + quicksum(slack1[j][k][_] * (2**_) for _ in range(len_slack1))
                    ) ** 2
        # inequality constraint 2: the capacity of an edge should be greater than or equal to the demand
        ineq_qubo2 = kw.core.BinaryExpression(coefficient={}, offset=0)
        ineq_constraint2 = []
        len_slack2 = math.ceil(math.log2(max(max_ujk) + 1))
        slack2 = kw.core.ndarray(
            (self.num_edge, self.num_cloud, len_slack2), "slack2", kw.core.Binary
        )
        for j in range(self.num_edge):
            ineq_constraint2.append([])
            for k in range(self.num_cloud):
                if y_jk[j][k] == 0:
                    ineq_constraint2[j].append(0)
                else:
                    ineq_constraint2[j].append(
                        y_jk[j][k] * self.edge_capacity
                        - quicksum(
                            self.demand[self.user_locations[i]] * y_ij[i][j]
                            for i in range(self.num_user)
                        )
                    )
                    ineq_qubo2 += (
                        ineq_constraint2[j][k]
                        + quicksum(slack2[j][k][_] * (2**_) for _ in range(len_slack2))
                    ) ** 2
        return ineq_qubo1, ineq_qubo2

    def _build_eq_constraint(self, x_edge, y_ij, y_ik, y_jk):
        # constraint 1: ensure that each user's service demand is assigned to only one location (either edge or cloud)
        constraint1 = 0
        for i in range(self.num_user):
            constraint1 += (
                quicksum(y_ij[i][j] for j in range(self.num_edge))
                + quicksum(y_ik[i][k] for k in range(self.num_cloud))
                - 1
            ) ** 2
        # constraint 2: initialize constraint expression
        constraint2 = 0
        for i in range(self.num_user):
            for j in range(self.num_edge):
                if self.coverage[(self.user_locations[i], self.edge_locations[j])] == 0:
                    y_ij[i][j] = 0
                else:
                    constraint2 += y_ij[i][j] * (1 - x_edge[j])
        # constraint 3: ensure the relationship between yjk and x_edge
        constraint3 = 0
        for j in range(self.num_edge):
            for k in range(self.num_cloud):
                constraint3 += y_jk[j][k] * (1 - x_edge[j])
        return constraint1, constraint2, constraint3

    def _build_objective_components(self, u_jk, x_edge, y_ij, y_ik):
        # objective function
        c_fix = quicksum(
            self.fix_cost_edge[self.edge_locations[j]] * x_edge[j]
            for j in range(self.num_edge)
        )
        c_var = self.var_cost_cloud * quicksum(
            quicksum(
                self.demand[self.user_locations[i]] * y_ik[i][k]
                for i in range(self.num_user)
            )
            for k in range(self.num_cloud)
        )
        c_var += self.var_cost_edge * quicksum(
            quicksum(
                self.demand[self.user_locations[i]] * y_ij[i][j]
                for i in range(self.num_user)
            )
            for j in range(self.num_edge)
        )
        c_var += (self.var_cost_cloud - self.var_cost_edge) * quicksum(
            quicksum(u_jk[j][k] for j in range(self.num_edge))
            for k in range(self.num_cloud)
        )
        c_tran = self.trans_cost_user_edge * quicksum(
            quicksum(
                self.demand[self.user_locations[i]]
                * self.distances[(self.user_locations[i], self.edge_locations[j])]
                * y_ij[i][j]
                for i in range(self.num_user)
            )
            for j in range(self.num_edge)
        )
        c_tran += self.trans_cost_user_cloud * quicksum(
            quicksum(
                self.demand[self.user_locations[i]]
                * self.distances[(self.user_locations[i], self.cloud_locations[k])]
                * y_ik[i][k]
                for i in range(self.num_user)
            )
            for k in range(self.num_cloud)
        )
        c_tran += self.trans_cost_edge_cloud * quicksum(
            quicksum(
                self.distances[(self.edge_locations[j], self.cloud_locations[k])]
                * u_jk[j][k]
                for j in range(self.num_edge)
            )
            for k in range(self.num_cloud)
        )
        return c_fix + c_tran + c_var

    def solve(
        self,
        max_iterations: int = 10,
        init_temp: float = 1e3,
        decay: float = 0.99,
        min_temp: float = 1e-4,
        iter_per_temp: int = 10,
    ):
        """
        Perform simulated annealing.
        """
        best = math.inf
        for num in range(max_iterations):
            solver = kw.classical.SimulatedAnnealingOptimizer(
                initial_temperature=init_temp,
                alpha=decay,
                cutoff_temperature=min_temp,
                iterations_per_t=iter_per_temp,
            )
            sol, val = solver.solve_qubo(self.qubo_model)
            feasible = self.recovery(sol)
            if feasible and val < best:
                best = val
                print(f"Iter {num}: best={val}")
        print(f"Optimal={best}")

    def recovery(self, sol_dict: dict) -> bool:
        """
        Validate solution via kw.core.get_val interface.
        """
        unsatisfied_count, _ = self.qubo_model.verify_constraint(sol_dict)
        return not unsatisfied_count


def assert_documented_results():
    """Check that the documented model-building path stays valid."""
    solver = CloudEdgeUserSolver()
    solver.prepare_data()
    solver.prepare_model()
    assert solver.qubo_model is not None


if __name__ == "__main__":
    # create an instance of the solver class
    solver = CloudEdgeUserSolver()

    # prepare the data
    solver.prepare_data()

    # prepare the qubo model, set the penalty coefficient lambda
    solver.prepare_model()

    # use simulated annealing to find the optimal solution
    solver.solve()