0 写在前面的

本文是对2022年华为软件精英挑战赛（普朗克计划）区域初赛的一个解读。首先说明的是本文的算法无法直接拿来参赛的，因为区域初赛的要求是不能调用其它的算法包，python的话只能用numpy，其它的第三方包一概不能使用。而本文的方法是基于数学规划模型，需要调用求解器来求解（例如：Gurobi, Cplex等）。因此本文是作为一个运筹优化爱好者的角度来看看如何解决一个实际的优化问题，如果能给对比赛感兴趣的同学一点启发的话那就更好了。关于区域初赛的说明文档，测试数据和本文算法的代码，我已经全部上传在Github，感兴趣的同学可以自行下载：

https://github.com/WenYuZhi/HuaWeiSoftwareCompetition2022

1 赛题背景

在视频直播场景中，网络成本是影响服务成本的关键因素之一，不同的流量调度方案 会产生不同的网络使用成本。本赛题以华为云视频直播服务流量调度问题为基础，并进行一定的抽象、调整和简化。参赛选手需要设计高效的调度算法，在满足客户要求 的前提下，通过对流量的合理调度，最小化网络使用成本。

详细的题目介绍就不赘述了，在参赛文件里边都有。下面我们直接进入模型部分。

2 优化模型

在参赛文件里边已经给出了基本的数学模型，我们这里是在参数文件的基础上进一步规范化一下模型以便于我们后边的求解，其本质和参赛文档中的数学模型基本一致：
目标函数：
$$J=\underset{x}{\min }\sum _{j=1}^N{s}_j(1)$$
Qos约束和非负约束：
$$\begin{gathered} X_{ij}^t=0,if{Y}_{ij}\ge Q,\forall t\in T,1\le i\le M,1\le j\le N(2) \\ {X}_{ij}^t\in Z^{+},\forall t\in T,1\le i\le M,1\le j\le N(3) \end{gathered}$$
需求约束：
$$\sum _{j=1}^N{X}_{ij}^t=D_i^t,\forall t\in T,1\le i\le M(4)$$
边缘节点带宽约束：
$$\sum _{i=1}^M{X}_{ij}^t\le C_j,\forall t\in T,1\le j\le N(5)$$
中间变量约束：
$$w_j^t=\sum _{i=1}^M{X}_{ij}^t,\forall t\in T,1\le j\le N(6)$$
95分位数约束：
s j = Quantile 95 { w j , t ∣ ∀ t ∈ T } ,    1 ≤ j ≤ N      ( 7 ) s_{j}^{}=\text{Quantile}95 \left\{ w_{j,t}|\forall t\in T \right\} ,\,\,1\le j\le N\;\;(7) sj​=Quantile95{wj,t​∣∀t∈T},1≤j≤N(7)
其中Quantile95表示95%分位数。

3 优化模型转化

从上面的优化问题中可以看出，目标函数和约束（1-6）都是线性的，因此相对来说是比较简单的约束。本题的关键就在于约束(7) 95分位数的约束，因为95分位数是一个比较难处理的约束。显然，在数学规划中我们很难直接用一个close form去表示出95分位数的表达式，那么在我们不能直接处理这个优化问题的时候往往会采用转化的方式。可以说在平常的数学规划模型中我们还是较少遇到类似95分位数的约束，这确实也是造成这个问题比较困难的主要原因。接下来我们就集中在如何处理这个95分位数约束上：

转化思路1：用最大值代替95分位数

根据95分位数的定义，它是一个比95%的数都要大的一个数。那么我们可以用最大值来近似一下95分位数，即将式（7）中 95分位数约束改为

s ^ j = max ⁡ { w j , t ∣ ∀ t ∈ T } ,    1 ≤ j ≤ N      ( 8 ) \hat{s}_{j}^{}=\max \left\{ w_{j,t}|\forall t\in T \right\} ,\,\,1\le j\le N\;\;(8) s^j​=max{wj,t​∣∀t∈T},1≤j≤N(8)
进而改写成新的目标函数：
$$J_1=\underset{x}{\min }\sum _{j=1}^N{\hat{s}}_j(9)$$
约束（2-6）继续保留， 约束（7）替换为约束（8）
因为最大值肯定是大于等于95分位数的，因此新的优化问题的目标函数大于原优化问题目标函数：$J_1\ge J$

转化思路2：采样95%的数据的最大值来代替95分位数

在上面的转化思路1中，虽然我们很容易通过最大值得到了95分位数的近似，但是这个近似还是比较简单粗暴的，我们实际上还没有利用95%的信息来做文章。因此接下来我们希望利用这个95%的信息，来得到比转化思路1的目标函数更接近原目标函数的新目标函数。
s ~ j = max ⁡ { w j , t ∣ ∀ t ∈ Φ } ,    1 ≤ j ≤ N      ( 10 ) Φ ⊆ T ,   ∣ Φ ∣ = ceil ( 0.95 ∣ T ∣ )      ( 11 ) \tilde{s}_j=\max \left\{ w_{j,t}|\forall t\in \varPhi \right\} ,\,\,1\le j\le N\;\;(10) \\ \varPhi \subseteq T,\ \left| \varPhi \right|=\text{ceil}\left( 0.95\left| T \right| \right) \;\;(11) s~j​=max{wj,t​∣∀t∈Φ},1≤j≤N(10)Φ⊆T, ∣Φ∣=ceil(0.95∣T∣)(11)
式（11）即表示集合 $\Phi$ 是集合 $T$的子集，并且集合 $\Phi$ 要拥有95%的 集合 $T$ 的元素（我们都默认集合 $\Phi$ 和集合 $T$都不存在重复的元素。其中 $\text{ceil}$ 表示上取整符号。

进而利用新定义出的 ${\tilde{s}}_j$ 给出新的目标函数：
$$J_2=\underset{x}{\min }\sum _{j=1}^N{\tilde{s}}_j(12)$$
约束（2-6）继续保留， 约束（7）替换为约束（10-11）

显然由于$J_2$是任选的95%的数据最大值，因此 $J\le J_2$ 。又因为 $J_1$ 是全部数据的最大值，而 $J_2$ 只是95%的数据的最大值，因此 $J_2\le J_1$ 。综上所述：
$$J\le J_2\le J_1(13)$$
至此，我们发现 $J_2$ 确实要比 $J_1$要更接近原目标函数。这也是因为我们在 $J_2$ 中利用到了95%的信息。恰恰是因为在约束（10-11）中我们只要求95%的数据最大值即可，这样使得目标函数 $J_2$ 比转化思路1中更加紧一些。

显然采样次数越多，则 $J_2$ 越接近 $J$ ，当采样次数趋于充分大的时候，必然有$J_2=J$。当然在实际操作层面上来讲，由于计算代价的问题我们不可能让采样次数太大。

4 基于 Gurobi 实现求解

本章我们就采用转化思路2中的优化问题：
目标函数
$$J_2=\underset{x}{\min }\sum _{j=1}^N{\tilde{s}}_j(14)$$
约束（2-6）和约束（10-11）
如下是在python中的主干代码：

from data_set import Data_Set
from model import Model

data_set = Data_Set()
data_set.load() # 读入数据
data_set.handle_data() # 数据预处理

flow_schedule = Model(data_set=data_set, qos_ub=400) # qos_ub表示qos的上限
flow_schedule.sample_time(quantile=0.95) # 生成采样集合
flow_schedule.set_objective()  # 设置目标函数
flow_schedule.set_constrs()    # 设置约束条件
flow_schedule.solve(time_limit = 300)  # 调用求解器求解
flow_schedule.get_solution()  # 获得最优解
flow_schedule.get_quantile()  # 获得原目标函数值
flow_schedule.save_solution()  # 保存解到本地文件
flow_schedule.write_model()   # 保存优化模型

具体优化模型实现的代码：

import gurobipy as gp
from gurobipy import GRB
import random
import math
from matplotlib.pyplot import axis
import numpy as np
import pandas as pd
FILE_PATH = './'+ "results//"

class Model():
    def __init__(self, data_set, qos_ub) -> None:
        self.qos, self.demand, self.bandwidth, self.site_name = data_set.qos, data_set.demand_values, data_set.bandwidth, data_set.site_name
        self.n_site, self.n_cust = data_set.n_site, data_set.n_cust # n_site 基站数量; n_cust:用户数
        self.qos_ub, self.n_time = qos_ub,  data_set.n_time # qos的上限; n_time:时刻数
        self.m = gp.Model('flow scheduled model') # 建立优化模型
        self.x = self.m.addVars(self.n_cust, self.n_site, self.n_time, lb=0, ub=GRB.INFINITY,vtype=GRB.INTEGER, name ='x') 
        self.w = self.m.addVars(self.n_site, self.n_time, lb=0, ub=GRB.INFINITY,vtype=GRB.INTEGER, name ='w')
        self.w_max = self.m.addVars(self.n_site, lb=0, ub=GRB.INFINITY,vtype=GRB.INTEGER, name ='w_max') # w_max:代表模型中的 s_hat
        
    def sample_time(self, quantile):
        self.quantile_point = math.ceil(quantile*self.n_time)
        self.sample_t = random.sample([i for i in range(self.n_time)], self.quantile_point) 
    
    def set_objective(self):
        self.m.setObjective(gp.quicksum(self.w_max[j] for j in range(self.n_site))) #目标函数 对应式(14)
    
    def set_qos_constrs(self):
        self.m.addConstrs((self.x[i,j,t] == 0 for i in range(self.n_cust) for j in range(self.n_site) for t in range(self.n_time) \
        if self.qos[i,j] >= self.qos_ub), name ='qos') # qos约束 对应式(2)
    
    def set_demand_constrs(self):
        self.m.addConstrs((gp.quicksum(self.x[i,j,t] for j in range(self.n_site)) == self.demand[i,t] for i in range(self.n_cust) \
        for t in range(self.n_time)), name = 'demand') # 需求约束 对应式(4)
    
    def set_bandwidth_constrs(self):
        self.m.addConstrs((self.w[j,t] <= self.bandwidth[j] for j in range(self.n_site) for t in range(self.n_time)), name = 'bandwidth') # 边缘节点带宽约束 对应式(5)
    
    def set_link_constrs(self):
        self.m.addConstrs((gp.quicksum(self.x[i,j,t] for i in range(self.n_cust)) == self.w[j,t] for j in range(self.n_site) \
        for t in range(self.n_time)), name ='link') # 中间变量约束 对应式(6)
    
    def set_quantile95_constrs(self):
        self.m.addConstrs((self.w_max[j] >= self.w[j,t] for j in range(self.n_site) for t in range(self.n_time) if t in self.sample_t), name = 'max') # 95分位数约束 对应式(10-11)

    def set_constrs(self):
        self.set_qos_constrs()
        self.set_demand_constrs()
        self.set_link_constrs()
        self.set_bandwidth_constrs()
        self.set_quantile95_constrs()
    
    def solve(self, time_limit):
        self.m.Params.TimeLimit = time_limit 
        self.m.optimize()
    
    def write_model(self):
        self.m.write('model.lp')
    
    def get_solution(self):
        vars = self.m.getVars() #获取求解结果
        self.x_dim = self.n_cust*self.n_site*self.n_time
        self.x_sol = np.array([v.x for v in vars[0:self.x_dim]]).reshape((self.n_cust, self.n_site, self.n_time))
        self.w_dim = self.n_site*self.n_time
        self.w_sol = np.array([v.x for v in vars[self.x_dim:self.x_dim+self.w_dim]]).reshape((self.n_site, self.n_time)).T
        self.w_max_dim = self.n_site
        self.w_max_sol = np.array([v.x for v in vars[self.x_dim+self.w_dim:self.x_dim+self.w_dim+self.w_max_dim]]).reshape((self.n_site, 1))
    
    def get_quantile(self):
        self.w_sorted = np.sort(self.w_sol, axis=0)
        self.w_q = self.w_sorted[self.quantile_point-1,:].reshape((self.n_site, 1))
        print("the 95 percent objective function: {}".format(np.sum(self.w_q)))

    def save_solution(self):
        pd.DataFrame(self.w_sol, columns=self.site_name).to_csv(FILE_PATH + "w_sol.csv")
        pd.DataFrame(np.hstack((self.w_max_sol, self.w_q, self.w_max_sol-self.w_q)), index=self.site_name, \
        columns=['max','95','diff']).to_csv(FILE_PATH + "w_max_sol.csv")
        pd.DataFrame(self.w_sorted, columns=self.site_name).to_csv(FILE_PATH + "w_sorted.csv")

参数设置：我们设置Qos的上限一律为400（这个参数可以在代码中修改qos_ub变量来实现）。
由于算法存在一定的随机性，这是因为我们每次会随机生成不同的采样集合 $\Phi$。由于集合$\Phi$在每次运行的时候都不一样，那么计算结果也可能会不一样。

我们运行了10次，结果如下所示：

上面运行结果最好的是182939，当然我们前面已经说了运行次数越多结果会越好。

5 总结

以上仅仅是我自己的想法，当然因为无法在线提交，所以也没有办法在线的测试算法的效果，姑且作为一个求解思路分享出来吧。