智能调度，高效赋能：揭秘算力网络资源优化分配之道#

算力网络，简单来说，就是将计算能力像水电一样进行输送的网络。想象一下，我们平时使用的手机、电脑等设备，都需要处理器来计算和运行各种程序。而算力网络就像是一个超级大的"处理器仓库"，它将分散在各地的计算资源整合起来，通过互联网进行优化分配，以满足不同用户和场景的计算需求。

举个例子，当我们观看高清视频、玩大型游戏或者进行复杂的数据分析时，这些操作都需要大量的计算能力。算力网络就像是一个智能的"电力系统"，能在短时间内为我们提供所需的计算资源，确保我们的使用体验流畅、高效。这样一来，我们就可以随时随地享受到强大的计算能力，而不用担心设备性能不足的问题。

算力网络面临的挑战#

  • 资源分配效率: 算力网络需要高效地分配计算资源以满足不同用户和应用的多样化需求。这要求网络能够实时监测资源使用情况，并动态调整资源分配策略，以避免资源闲置或过载。

  • 网络延迟: 计算任务的快速响应需要低延迟的网络环境。算力网络必须解决数据在网络中传输的延迟问题，特别是在处理实时性要求高的应用时，如自动驾驶等。

  • 成本要低: 搭建和维护算力网络需要投入大量资金，所以需要在保证性能的同时，尽量降低成本。

算力网络的组成部分#

算力网络主要由三个部分组成：

  • 端侧设备: 也就是各种智能设备，它们负责收集数据。

  • 边缘服务器: 位于设备附近，负责处理部分数据，并减轻云端服务器的负担。

  • 云端服务器: 位于数据中心，拥有强大的计算能力，负责处理更复杂的数据。

赛题内容#

考虑一个特定区域的算力网络布局优化问题。将区域划分成若干相邻正方形网格，算力需求分布数据提供了各个网格内的算力需求。数据中的坐标值为所在网格的中心坐标，为简化问题，将各网格内的算力需求点统一到所在网格中心点（即可认为是一个网格对应一个需求点）。

网格内的算力需求由端侧设备产生，端侧设备是连接到网络的终端设备，例如传感器、智能手机、工业机器人等。算力网络中的算力需求由边缘服务器和云端服务器满足。边缘服务器位于网络的"边缘"，通常靠近终端用户或设备。它们的任务是在离用户更近的地方处理数据，以提高响应速度和效率。边缘服务器可以更快地处理请求，因为它们距离用户更近。边缘服务器可以减轻核心云基础设施的负担，提高整体操作效率。云端服务器位于远离用户的数据中心，具有强大的计算和存储能力。当边缘服务器容量不足时，云端服务器可以作为补充。边缘服务器和云端服务器之间的协同作用可以优化整个系统的性能和可靠性。

这是一个图片的描述

问题1: 在问题1中我们只考虑由边缘服务器满足计算区域内的需求。 假设要在算力需求分布的网格区域内开设2个边缘服务器，每个边缘服务器的覆盖半径为1。 要求建立QUBO模型，确定在哪些点放置边缘服务器可以覆盖最多的算力需求。

问题2: 当边缘服务器无法满足算力需求时，将由上游云端服务器提供计算服务。现在网格区域外部增设云端服务器，端侧和边缘服务器可以选择连接到云端节点，当边缘服务器接受的计算需求超过容量限制时，边缘服务器多余的需求将直接被分配到云端服务器。 每个端侧节点的算力需求都必须被满足，且只能由一个服务器服务，这个计算节点可以是云端节点，也可以是边缘服务器。由于边缘服务器存在资源容量上限约束，假设每个边缘服务器的可用计算资源容量为12，云端服务器的可用计算资源容量为无穷大， 即不考虑云端服务器的可用计算资源限制。服务器存在一定的覆盖半径，假设边缘服务器的覆盖半径为3，云端服务器的覆盖半径为无穷大，即不考虑云端服务器的覆盖范围。

开设边缘服务器通常需要成本，该成本由固定成本，计算成本和传输成本组成。其中固定成本与是否开设以及开设位置有关；计算成本与请求的计算资源量成正比，计算方法为：单位计算成本乘以计算量。云端服务器的单位计算量成本为1，边缘服务器的单位计算量成本为2。 同时，从端侧到边侧、从边侧到云侧、从端侧到云侧之间的传输也存在传输成本，传输成本计算方式为：传输的算力需求量乘以传输距离乘以单位传输成本，其中传输距离的计算使用欧式距离保留两位小数，按照单边距离计算（不考虑往返传输）。从端侧到边侧以及从边侧到云侧的单位传输成本为1，从端侧到云侧的单位传输成本为2。

如果要满足区域内全部的端侧算力需求，建立QUBO模型，求解给出整体成本最小的算力网络布局，即边缘服务器的位置、数量，以及端侧到边侧、边侧到云侧和端侧到云侧节点之间的连接关系。

问题本质#

在算力网络布局优化问题中，表面上看，我们可能会联想到经典的资源调度问题或网络流优化问题，并尝试通过传统的贪心算法或动态规划方法来求解"如何在不同区域间分配资源以最小化成本"。这些方法或许能够帮助我们找到某一时刻的局部最优解，但这个问题的复杂性远超表面。我们不仅需要考虑每个计算节点的资源分配，还要全局性地优化整个网络的布局，保证每个资源都得到合理利用，并且在多个区域之间保持高效的计算和通信。

深入分析#

算力网络布局优化问题的核心是如何在整个网络中找到一种资源分配和节点布局的最优方案，类似于在图G上寻找最优的计算节点布局和通信路径组合。这实际上是一个组合优化问题，需要从所有可能的网络布局中挑选出最优解，而这些可能的布局数目是巨大的。例如，如果我们有n个潜在的边缘服务器部署位置和m个计算任务点，那么这些位置和任务的组合数量将以指数级别增长，即n^m种可能的布局方式。

显然，简单的穷举法是不现实的，因为即使是中等规模的网络，其组合数量也是天文数字，远远超出了经典计算机的处理能力。就像旅行商问题中的路径选择一样，我们在这里面临着一个排列组合的难题，网络中的每个节点和资源配置都可能影响整体的优化结果。

换句话说，随着问题规模的扩大，算法的时间复杂度呈现非多项式增长，难以在合理时间内求解。由此可见，问题的关键并不在于如何局部优化某一个资源配置，而在于如何全局性地优化整个算力网络的布局，找到那条能够最小化整体计算成本的最优布局路径。这个问题的复杂性不仅仅是计算量的增加，更在于需要在多维度上同时进行优化，寻找全局最优解。

第一问参考模型#

第一问重述#

在一个4x4的网格中，每个网格代表一定的计算需求区域。我们的目标是确定在哪些网格中部署边缘计算节点，以便最大化覆盖的用户需求。

符号定义#

  • R=1R=1: 边缘计算节点的覆盖半径。

  • P=2P=2: 计划部署的边缘计算节点个数。

  • Lu\mathcal{L}_{u}: 用户位置集合。

  • Le\mathcal{L}_{e}: 候选边缘节点位置集合。

  • demandl(l∈Lu)\text{demand}_{l} (l\in \mathcal{L}_{u}): 位于网格 ll 处的算力需求量。

  • dijd_{ij}: 网格 ii 和 jj 之间的距离。

  • αij=1dij≤R\alpha_{ij}=\mathbb{1}_{d_{ij}\leq R}: 网格 ii 和网格 jj 之间的距离是否不超过边缘节点覆盖半径 RR。

决策变量#

  • xix_{i}: 一个二进制变量，表示是否在网格 ii 部署边缘计算节点。

  • ziz_{i}: 一个二进制变量，表示网格 ii 是否被覆盖。

目标函数#

最大化总覆盖的计算需求量：

max⁡∑i∈Luzi⋅demandi\max \sum_{i\in \mathcal{L}_{u}} {z_{i}\cdot \text{demand}_{i}}

约束条件#

  • 每个用户网格的覆盖状态不超过其周围边缘节点的覆盖状态：

    zi≤∑j∈Leαij⋅xj,∀i∈Luz_{i}\leq \sum_{j \in \mathcal{L}_{e}}{\alpha_{ij}\cdot x_{j}}, \forall i\in \mathcal{L}_{u}

  • 部署的边缘计算节点总数等于 PP：

    ∑j∈Lexj=P\sum_{j\in \mathcal{L}_{e}} x_{j} = P

简化模型（容斥原理）#

当 P=2P=2 时，我们可以使用容斥原理来简化模型：

max⁡∑i∈Lu(∑j∈Leαijxj−∑j,k∈Leαijαikxjxk)demandis.t.∑j∈Lexj=2\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}

第二问参考模型#

问题概述#

  • 6×6的网格，每个网格内都有一定的计算需求(只有11个网格中的计算需求为非零值)，需要满足所有计算需求

  • 5个开设边缘服务器的候选位置，边缘服务器具有计算容量限制; 当边缘服务器接收的计算请求超过容量时，将多出的计算请求发送给云端服务器

  • 1个位置开设了云端服务器，云端服务器无容量限制

  • 用户可以连接边缘服务器或者云端服务器(只能连一个); 边缘服务器容量不够可以连接云端服务器

  • 最小化成本：固定成本 + 可变成本(计算成本) + 传输成本(=单位成本*距离*传输量)

示意图

符号定义#

  • CedgeC_{\text{edge}}: 边侧服务器容量。

  • c-fix-edgej\text{c-fix-edge}_{j}: 边缘计算节点开设在网格 jj 的固定成本。

  • c-var-cloud\text{c-var-cloud}: 云侧节点的单位计算量成本。

  • c-var-edge\text{c-var-edge}: 边侧节点的单位计算量成本。

  • c-tranu,e\text{c-tran}_{u,e}: 从端侧到边侧的单位传输成本。

  • c-tranu,c\text{c-tran}_{u,c}: 从端侧到云侧的单位传输成本。

  • c-trane,c\text{c-tran}_{e,c}: 从边侧到云侧的单位传输成本。

  • dijd_{ij}: 网格 ii 和 jj 之间的距离。

  • did_{i}: 网格 ii 和云端服务器之间的距离。

  • αij=1dij≤R\alpha_{ij}=\mathbb{1}_{d_{ij} \leq R}: 网格 ii 和网格 jj 之间的距离是否不超过边缘节点覆盖半径 RR。

中间变量#

  • yiju,ey^{u,e}_{ij}: 网格 ii 的需求是否由位于网格 jj 的边侧服务器服务。

  • yje,cy^{e,c}_{j}: 网格 jj 的边侧服务器是否连接云端服务器。

  • yiu,cy^{u,c}_{i}: 网格 ii 的需求是否由云端服务器服务。

  • uj∈Nu_{j} \in \mathbb{N}: 超过边侧服务器 jj 的容量的计算需求中，由云端服务器服务的数量。

决策变量#

  • xedge,jx_{\text{edge}, j}: 是否在位置 jj 处开设边侧服务器。

数学模型#

目标函数:

min⁡Costfix+Costvar+Costtran\min \text{Cost}_{\text{fix}} + \text{Cost}_{\text{var}} + \text{Cost}_{\text{tran}}

其中

Costfix=∑j∈Lec-fix-edgej xedge,j\begin{aligned} \text{Cost}_{\text{fix}} = \sum_{j \in \mathcal{L}_{e}} \text{c-fix-edge}_{j}\, x_{\text{edge},j} \end{aligned}
Costvar=  c-var-cloud∑i∈Ludemandi yiu,c+c-var-edge∑j∈Le∑i∈Ludemandi yiju,e+(c-var-cloud−c-var-edge)∑j∈Leuj\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}
Costtran=  c-tranu,e∑i∈Lu∑j∈Leyiju,e demandi dij+c-tranu,c∑i∈Luyiu,c demandi di+c-trane,c∑j∈Ledj uj\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}

本模型中，我们将 uju_{j} 表示为

uj=yje,c⋅(∑i∈Ludemandi⋅yiju,e−Cedge)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)

同时我们添加约束条件，使得当且仅当边侧服务器 jj 收到的计算需求超过其容量限制时，yjku,ey^{u,e}_{jk} 才会取值为 11，否则取值为 00。虽然 uju_{j} 的表达式中出现了二次项，但是本模型中云端服务器没有容量限制，uju_{j} 只出现在目标函数中，而不会出现在约束中，从而保证了转化成QUBO模型时不会出现高阶项。

约束条件:

  • 算力需求点被分配到边侧或云侧，且只被一个服务

    ∑j∈Leyiju,e+yiu,c=1,∀i∈Lu\sum_{j \in \mathcal{L}_{e}}{y^{u,e}_{ij}} + y^{u,c}_{i} = 1, \forall i \in \mathcal{L}_{u}

  • 覆盖关系，且只有开了边侧才可以从需求点连接

    yiju,e≤αij⋅xedge,j,∀i∈Lu,j∈Ley^{u,e}_{ij} \leq \alpha_{ij} \cdot x_{\text{edge}, j}, \forall i \in \mathcal{L}_{u}, j \in \mathcal{L}_{e}

  • 边侧连接云侧和开启边侧的关系

    yje,c≤xedge,j,∀j∈Ley^{e,c}_{j} \leq x_{\text{edge}, j}, \forall j \in \mathcal{L}_{e}

  • 当且仅当边侧服务器 jj 收到的计算需求超过其容量限制时，yjku,ey^{u,e}_{jk} 取值为 11

    yje,c⋅Cedge≤∑i∈Ludemandi⋅αij⋅yiju,e(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)}
    ∑i∈Ludemandi⋅αij⋅yiju,e−Cedge≤(max-uj[j]−Cedge)⋅yje,c(2)\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)}

    其中 max-uj[j]\text{max-uj}[j] 为边侧服务器 jj 接收的计算需求的上界，这里我们令

    max-uj[j]=∑i∈Ludemandi⋅αij\text{max-uj}[j] = \sum_{i \in \mathcal{L}_{u}}{\text{demand}_{i} \cdot \alpha_{ij}}

预处理#

对于边侧服务器候选位置 jj，如果其接收的计算需求始终小于边侧服务器容量限制，例如

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

成立时，可以令 yje,c=0y^{e,c}_{j} = 0，并且对于该候选位置，可以不考虑不等式约束 (1) 和 (2)，减少比特数（松弛变量数量）。

第一问代码#

  1import math
  2import numpy as np
  3import kaiwu as kw
  4
  5
  6class EdgeCoveragePlanner:
  7    """Solver for Edge Coverage Planning problem."""
  8
  9    def __init__(
 10        self,
 11        num_nodes: int = 2,
 12        coverage_range: float = 2.0,
 13        penalty: float = 100.0,
 14        num_slack_bins: int = 1,
 15    ) -> None:
 16        """Initialize EdgeCoveragePlanner parameters.
 17
 18        Args:
 19            num_nodes: Number of edge nodes to be selected.
 20
 21            coverage_range: Node coverage threshold.
 22
 23            penalty: penalty coefficient.
 24
 25            num_slack_bins: Slack binary digits, used to override constraints.
 26        """
 27        self.num_nodes = num_nodes
 28        self.coverage_range = coverage_range
 29        self.penalty = penalty
 30        self.num_slack_bins = num_slack_bins
 31
 32        # Data containers
 33        self.demand = {}
 34        self.locations = []
 35        self.distances = {}
 36        self.coverage = {}
 37        self.dimension = 0
 38
 39        # QUBO model
 40        self.model = None
 41
 42    def prepare_data(self, demand_data):
 43        """prepare demand, distance, and coverage matrix data.
 44
 45        Args:
 46            demand_data: dict, the key is the coordinate string 'i, j',
 47                        and the value corresponds to the calculation requirement.
 48        """
 49        self.demand = demand_data
 50        # Create a list containing all location coordinates
 51        self.locations = list(demand_data.keys())
 52        # Calculate the number of location coordinates
 53        self.dimension = int(math.sqrt(len(self.demand)))
 54        # Iterate through all location coordinates to calculate distances between pairs
 55        self.distances = {
 56            (i, j): np.linalg.norm(
 57                np.array([int(coord) for coord in i.split(",")])
 58                - np.array([int(coord) for coord in j.split(",")])
 59            )
 60            for i in self.locations
 61            for j in self.locations
 62        }
 63        # Initialize a dictionary to store the coverage relationship between locations
 64        self.coverage = {
 65            (i, j): int(dist <= self.coverage_range)
 66            for (i, j), dist in self.distances.items()
 67        }
 68
 69    def prepare_model(self):
 70        """Building a Qubo Model"""
 71        # Create binary variable arrays x and z to represent edge computing node locations
 72        # and demand coverage conditions
 73        var_x = kw.core.ndarray((self.dimension, self.dimension), "x", kw.core.Binary)
 74        var_z = kw.core.ndarray((self.dimension, self.dimension), "z", kw.core.Binary)
 75        limit_z = kw.core.zeros((self.dimension, self.dimension))
 76        # Initialize slack variables
 77        slack = kw.core.ndarray(
 78            (self.dimension, self.dimension),
 79            "slack",
 80            kw.core.Integer,
 81            (0, (2**self.num_slack_bins) - 1),
 82        )
 83        # Define the objective function to minimize the total computing power demand
 84        obj = kw.core.quicksum(
 85            self.demand[f"{i + 1},{j + 1}"] * var_z[i, j]
 86            for i in range(self.dimension)
 87            for j in range(self.dimension)
 88        )
 89        # Iterate through all location coordinates to build the second constraint
 90        for i in range(self.dimension):
 91            for j in range(self.dimension):
 92                # For each location, the demand coverage z[i] should be less than or equal to
 93                # the total coverage provided by all edge computing nodes
 94                limit_z[i, j] = kw.core.quicksum(
 95                    self.coverage[f"{i + 1},{j + 1}", f"{i1 + 1},{j1 + 1}"]
 96                    * var_x[i1, j1]
 97                    for i1 in range(self.dimension)
 98                    for j1 in range(self.dimension)
 99                )
100
101        self.model = kw.core.QuboModel()
102        self.model.set_objective(-obj)
103        self.model.add_constraint(
104            var_x.sum() == self.num_nodes, "c1", penalty=self.penalty
105        )
106        self.model.add_constraint(
107            var_z <= limit_z, "c2", penalty=self.penalty, slack_var_expr=slack
108        )
109
110    def solve(self):
111        """Solving the QUBO model.
112
113        Returns:
114            tuple: Result dictionary and Result dictionary.
115
116            - dict: Result dictionary. The key is the variable name, and the value is the corresponding spin value.
117
118            - float: qubo value.
119        """
120        # Perform the Simulated Annealing algorithm
121        _solver = kw.classical.SimulatedAnnealingOptimizer(
122            initial_temperature=100000,
123            alpha=0.99,
124            cutoff_temperature=0.0001,
125            iterations_per_t=100,
126            rand_seed=10,
127        )
128        _sol_dict, _qubo_value = _solver.solve_qubo(self.model)
129        return _sol_dict, float(_qubo_value)
130
131    def recovery(self, sol_dict):
132        """Verify whether the solution is feasible"""
133        return self.model.verify_constraint(sol_dict)
134
135
136def assert_documented_results():
137    """Check that the documented model-building path stays valid."""
138    demand_data = {
139        "1,1": 38,
140        "1,2": 22,
141        "1,3": 65,
142        "1,4": 56,
143        "2,1": 53,
144        "2,2": 48,
145        "2,3": 76,
146        "2,4": 46,
147        "3,1": 56,
148        "3,2": 36,
149        "3,3": 7,
150        "3,4": 29,
151        "4,1": 50,
152        "4,2": 37,
153        "4,3": 48,
154        "4,4": 40,
155    }
156    solver = EdgeCoveragePlanner()
157    solver.prepare_data(demand_data)
158    solver.prepare_model()
159    assert solver.model is not None
160
161
162if __name__ == "__main__":
163    # Store computing power demand data in the DEM dictionary,
164    # where the keys are location coordinates and the values are computing power demands
165    # fmt: off
166    demand_data = {
167        '1,1': 38, '1,2': 22, '1,3': 65, '1,4': 56,
168        '2,1': 53, '2,2': 48, '2,3': 76, '2,4': 46,
169        '3,1': 56, '3,2': 36, '3,3': 7,  '3,4': 29,
170        '4,1': 50, '4,2': 37, '4,3': 48, '4,4': 40
171    }
172    # fmt: on
173    # Create an instance of the SPQCSolver class
174    solver = EdgeCoveragePlanner()
175
176    # Prepare data
177    solver.prepare_data(demand_data)
178
179    # Prepare the QUBO model with a specified penalty coefficient lambda
180    solver.prepare_model()
181
182    # Use the Simulated Annealing algorithm to find the optimal solution
183    best_sol_dict, qubo_value = solver.solve()
184
185    # Recover the original problem solution from the QUBO solution and check its feasibility
186    unsatisfied_count, result_dict = solver.recovery(best_sol_dict)
187    if unsatisfied_count == 0:
188        print("Find a feasible solution")
189        print("Objective value:", -qubo_value)
190    else:
191        print("No feasible solution")

第二问代码#

  1import math
  2from typing import Tuple
  3import numpy as np
  4import kaiwu as kw
  5from kaiwu.core import quicksum
  6
  7
  8class CloudEdgeUserSolver:
  9    """Solver for cloud-edge-user cost-minimization QUBO model."""
 10
 11    def __init__(self, edge_capacity: int = 12, edge_radius: float = 3.0) -> None:
 12        """Initialize core parameters and placeholders."""
 13        # capacities and radii
 14        self.edge_capacity = edge_capacity
 15        self.edge_radius = edge_radius
 16
 17        # Service node configurations
 18        self.user_locations = [
 19            "1,1",
 20            "1,4",
 21            "2,6",
 22            "3,3",
 23            "3,5",
 24            "4,4",
 25            "5,5",
 26            "6,3",
 27            "6,6",
 28        ]
 29        self.edge_locations = ["6,1", "2,3", "4,5", "6,5"]
 30        self.cloud_locations = ["4,0"]
 31
 32        # edge server costs
 33        self.fix_cost_edge = {"6,1": 60, "2,3": 42, "4,5": 46, "6,5": 54}
 34
 35        # variable costs
 36        self.var_cost_cloud = 1
 37        self.var_cost_edge = 2
 38
 39        # Unit transmission cost
 40        self.trans_cost_user_edge = 1
 41        self.trans_cost_user_cloud = 2
 42        self.trans_cost_edge_cloud = 1
 43
 44        # Demand data
 45        self.demand = {
 46            "1,1": 7,
 47            "1,4": 4,
 48            "2,6": 5,
 49            "3,3": 9,
 50            "3,5": 8,
 51            "4,4": 11,
 52            "5,5": 1,
 53            "6,3": 7,
 54            "6,6": 5,
 55        }
 56
 57        # Precomputed data
 58        self.num_user = 0
 59        self.num_edge = 0
 60        self.num_cloud = 0
 61        self.distances = {}
 62        self.coverage = {}
 63        self.qubo_model = None
 64
 65    def prepare_data(self):
 66        """Prepare data for QUBO model."""
 67        # initialize user demands
 68        # fmt: off
 69        loc_set_inner_full = ['1,1', '1,2', '1,3', '1,4', '1,5', '1,6',
 70                              '2,1', '2,2', '2,3', '2,4', '2,5', '2,6',
 71                              '3,1', '3,2', '3,3', '3,4', '3,5',
 72                              '4,1', '4,2', '4,3', '4,4', '4,5', '4,6',
 73                              '5,1', '5,2', '5,3', '5,4', '5,5', '5,6',
 74                              '6,1', '6,2', '6,3', '6,4', '6,5', '6,6']
 75        # fmt: on
 76        # all nodes
 77        loc_set = loc_set_inner_full + self.cloud_locations
 78
 79        self.num_cloud = len(self.cloud_locations)
 80        self.num_edge = len(self.edge_locations)
 81        self.num_user = len(self.user_locations)
 82
 83        # compute distances
 84        for i in loc_set:
 85            for j in loc_set:
 86                self.distances[(i, j)] = round(
 87                    np.sqrt(
 88                        (
 89                            int(i.split(",", maxsplit=1)[0])
 90                            - int(j.split(",", maxsplit=1)[0])
 91                        )
 92                        ** 2
 93                        + (
 94                            int(i.split(",", maxsplit=1)[1])
 95                            - int(j.split(",", maxsplit=1)[1])
 96                        )
 97                        ** 2
 98                    ),
 99                    2,
100                )
101
102        # compute coverage relationships
103        for i in loc_set:
104            for j in loc_set:
105                if self.distances[(i, j)] <= self.edge_radius:
106                    self.coverage[(i, j)] = 1
107                else:
108                    self.coverage[(i, j)] = 0
109
110    def prepare_model(
111        self,
112        eq_penalties: Tuple[float, float, float] = (1e4, 1e4, 1e4),
113        ineq_penalties: Tuple[float, float] = (1e4, 1e4),
114    ) -> None:
115        """
116        Build QUBO with 3 equality and 2 inequality constraints.
117        """
118        # compute the maximum demand for each edge
119        max_ujk = []
120        for j in self.edge_locations:
121            max_ujk.append(
122                sum(
123                    self.demand[i] * self.coverage[i, j] for i in self.user_locations
124                )
125            )
126
127        # initialize decision variables
128        x_edge = kw.core.ndarray(self.num_edge, "x_edge", kw.core.Binary)
129        y_ij = kw.core.ndarray((self.num_user, self.num_edge), "yij", kw.core.Binary)
130        y_jk = kw.core.ndarray((self.num_edge, self.num_cloud), "yjk", kw.core.Binary)
131        y_ik = kw.core.ndarray((self.num_user, self.num_cloud), "yik", kw.core.Binary)
132
133        # for each edge, if the maximum demand does not exceed the edge's capacity, set the corresponding yjk to 0
134        for j in range(self.num_edge):
135            if max_ujk[j] <= self.edge_capacity:
136                y_jk[j][0] = 0
137
138        # initialize ujk related variables
139        u_jk = np.zeros(
140            shape=(self.num_edge, self.num_cloud), dtype=kw.core.BinaryExpression
141        )
142        ujk_residual = np.zeros(shape=self.num_edge, dtype=kw.core.BinaryExpression)
143        for j in range(self.num_edge):
144            ujk_residual[j] = (
145                quicksum(
146                    self.demand[self.user_locations[i]] * y_ij[i][j]
147                    for i in range(self.num_user)
148                )
149                - self.edge_capacity
150            )
151            for k in range(self.num_cloud):
152                u_jk[j][k] = y_jk[j][k] * ujk_residual[j]
153        # build objective
154        obj = self._build_objective_components(u_jk, x_edge, y_ij, y_ik)
155        # build equality constraint
156        constraint1, constraint2, constraint3 = self._build_eq_constraint(
157            x_edge, y_ij, y_ik, y_jk
158        )
159        # build inequality constraint
160        ineq_qubo1, ineq_qubo2 = self._build_ineq_constraint(
161            max_ujk, ujk_residual, y_ij, y_jk
162        )
163
164        # building the final model
165        self.qubo_model = kw.core.QuboModel()
166        self.qubo_model.set_objective(obj)
167        self.qubo_model.add_constraint(constraint1, name="c1", penalty=eq_penalties[0])
168        self.qubo_model.add_constraint(constraint2, name="c2", penalty=eq_penalties[1])
169        self.qubo_model.add_constraint(constraint3, name="c3", penalty=eq_penalties[2])
170        self.qubo_model.add_constraint(ineq_qubo1, name="c4", penalty=ineq_penalties[0])
171        self.qubo_model.add_constraint(ineq_qubo2, name="c5", penalty=ineq_penalties[1])
172
173    def _build_ineq_constraint(self, max_ujk, ujk_residual, y_ij, y_jk):
174        # inequality constraint 1: after subtracting the edge's maximum capacity from demand,
175        # the yjk constraint should hold
176        ineq_constraint1 = []
177        ineq_qubo1 = kw.core.BinaryExpression(coefficient={}, offset=0)
178        len_slack1 = math.ceil(math.log2(max(max_ujk) + 1))
179        slack1 = kw.core.ndarray(
180            (self.num_edge, self.num_cloud, len_slack1), "slack1", kw.core.Binary
181        )
182        for j in range(self.num_edge):
183            ineq_constraint1.append([])
184            for k in range(self.num_cloud):
185                if y_jk[j][k] == 0:
186                    ineq_constraint1[j].append(0)
187                else:
188                    ineq_constraint1[j].append(
189                        ujk_residual[j] - (max_ujk[j] - self.edge_capacity) * y_jk[j][k]
190                    )
191                    ineq_qubo1 += (
192                        ineq_constraint1[j][k]
193                        + quicksum(slack1[j][k][_] * (2**_) for _ in range(len_slack1))
194                    ) ** 2
195        # inequality constraint 2: the capacity of an edge should be greater than or equal to the demand
196        ineq_qubo2 = kw.core.BinaryExpression(coefficient={}, offset=0)
197        ineq_constraint2 = []
198        len_slack2 = math.ceil(math.log2(max(max_ujk) + 1))
199        slack2 = kw.core.ndarray(
200            (self.num_edge, self.num_cloud, len_slack2), "slack2", kw.core.Binary
201        )
202        for j in range(self.num_edge):
203            ineq_constraint2.append([])
204            for k in range(self.num_cloud):
205                if y_jk[j][k] == 0:
206                    ineq_constraint2[j].append(0)
207                else:
208                    ineq_constraint2[j].append(
209                        y_jk[j][k] * self.edge_capacity
210                        - quicksum(
211                            self.demand[self.user_locations[i]] * y_ij[i][j]
212                            for i in range(self.num_user)
213                        )
214                    )
215                    ineq_qubo2 += (
216                        ineq_constraint2[j][k]
217                        + quicksum(slack2[j][k][_] * (2**_) for _ in range(len_slack2))
218                    ) ** 2
219        return ineq_qubo1, ineq_qubo2
220
221    def _build_eq_constraint(self, x_edge, y_ij, y_ik, y_jk):
222        # constraint 1: ensure that each user's service demand is assigned to only one location (either edge or cloud)
223        constraint1 = 0
224        for i in range(self.num_user):
225            constraint1 += (
226                quicksum(y_ij[i][j] for j in range(self.num_edge))
227                + quicksum(y_ik[i][k] for k in range(self.num_cloud))
228                - 1
229            ) ** 2
230        # constraint 2: initialize constraint expression
231        constraint2 = 0
232        for i in range(self.num_user):
233            for j in range(self.num_edge):
234                if self.coverage[(self.user_locations[i], self.edge_locations[j])] == 0:
235                    y_ij[i][j] = 0
236                else:
237                    constraint2 += y_ij[i][j] * (1 - x_edge[j])
238        # constraint 3: ensure the relationship between yjk and x_edge
239        constraint3 = 0
240        for j in range(self.num_edge):
241            for k in range(self.num_cloud):
242                constraint3 += y_jk[j][k] * (1 - x_edge[j])
243        return constraint1, constraint2, constraint3
244
245    def _build_objective_components(self, u_jk, x_edge, y_ij, y_ik):
246        # objective function
247        c_fix = quicksum(
248            self.fix_cost_edge[self.edge_locations[j]] * x_edge[j]
249            for j in range(self.num_edge)
250        )
251        c_var = self.var_cost_cloud * quicksum(
252            quicksum(
253                self.demand[self.user_locations[i]] * y_ik[i][k]
254                for i in range(self.num_user)
255            )
256            for k in range(self.num_cloud)
257        )
258        c_var += self.var_cost_edge * quicksum(
259            quicksum(
260                self.demand[self.user_locations[i]] * y_ij[i][j]
261                for i in range(self.num_user)
262            )
263            for j in range(self.num_edge)
264        )
265        c_var += (self.var_cost_cloud - self.var_cost_edge) * quicksum(
266            quicksum(u_jk[j][k] for j in range(self.num_edge))
267            for k in range(self.num_cloud)
268        )
269        c_tran = self.trans_cost_user_edge * quicksum(
270            quicksum(
271                self.demand[self.user_locations[i]]
272                * self.distances[(self.user_locations[i], self.edge_locations[j])]
273                * y_ij[i][j]
274                for i in range(self.num_user)
275            )
276            for j in range(self.num_edge)
277        )
278        c_tran += self.trans_cost_user_cloud * quicksum(
279            quicksum(
280                self.demand[self.user_locations[i]]
281                * self.distances[(self.user_locations[i], self.cloud_locations[k])]
282                * y_ik[i][k]
283                for i in range(self.num_user)
284            )
285            for k in range(self.num_cloud)
286        )
287        c_tran += self.trans_cost_edge_cloud * quicksum(
288            quicksum(
289                self.distances[(self.edge_locations[j], self.cloud_locations[k])]
290                * u_jk[j][k]
291                for j in range(self.num_edge)
292            )
293            for k in range(self.num_cloud)
294        )
295        return c_fix + c_tran + c_var
296
297    def solve(
298        self,
299        max_iterations: int = 10,
300        init_temp: float = 1e3,
301        decay: float = 0.99,
302        min_temp: float = 1e-4,
303        iter_per_temp: int = 10,
304    ):
305        """
306        Perform simulated annealing.
307        """
308        best = math.inf
309        for num in range(max_iterations):
310            solver = kw.classical.SimulatedAnnealingOptimizer(
311                initial_temperature=init_temp,
312                alpha=decay,
313                cutoff_temperature=min_temp,
314                iterations_per_t=iter_per_temp,
315            )
316            sol, val = solver.solve_qubo(self.qubo_model)
317            feasible = self.recovery(sol)
318            if feasible and val < best:
319                best = val
320                print(f"Iter {num}: best={val}")
321        print(f"Optimal={best}")
322
323    def recovery(self, sol_dict: dict) -> bool:
324        """
325        Validate solution via kw.core.get_val interface.
326        """
327        unsatisfied_count, _ = self.qubo_model.verify_constraint(sol_dict)
328        return not unsatisfied_count
329
330
331def assert_documented_results():
332    """Check that the documented model-building path stays valid."""
333    solver = CloudEdgeUserSolver()
334    solver.prepare_data()
335    solver.prepare_model()
336    assert solver.qubo_model is not None
337
338
339if __name__ == "__main__":
340    # create an instance of the solver class
341    solver = CloudEdgeUserSolver()
342
343    # prepare the data
344    solver.prepare_data()
345
346    # prepare the qubo model, set the penalty coefficient lambda
347    solver.prepare_model()
348
349    # use simulated annealing to find the optimal solution
350    solver.solve()
On this page