考虑灾害响应的应急设施预定位选址

考虑灾害响应的应急设施预定位选址

真实背景：专业要求完成此课程设计并通过答辩，否则就会被强迫安排去实习~~~/(ㄒoㄒ)/~~

1.问题背景

近年来，随着全球气候异常多变，导致灾难频繁发生。因此，有必要根据历史灾难救援数据，进行灾害响应的应急设施预定位。在灾害发生前，应急设施预定位可以根据潜在的灾害风险和历史灾害的情况，结合地理信息系统、环境监测数据等科技手段，制定出最优的预案（应急储备仓库选址方案），并预留相应的资源和设备，以便在灾害发生后能够迅速响应。在灾害发生后，应急设施预定位可以通过实时监测灾情和通信系统，及时获取灾区信息和需求，快速响应并调配资源和人员，提供救援和支持。同时，应急设施预定位还可以在灾后进行复杂情况的应对和恢复工作，促进灾区重建和发展。
基于上述背景，相关决策人员为了在发生紧急情况时，能够快速有效地响应并提供紧急救援和支持，计划在全区域预设8个候选地点，作为建立应急设施储备库的可选方案。由此提高应急响应的效率和准确性，减少人员和财产的损失。

应急网络设计

目前，针对某地区潜在的48个灾难发生地的应急物资需求，相关决策人员选择了8个地点作为侯选位置，以备未来建设（或开放）物资储备库。令G(V, E)代表一个应急物资调配网络。V表示网络中节点的集合，其中包括48个潜在灾难点和8个候选应急设施库地点。E表示网络中的边的集合。总共考虑10个可能发生的应急场景，每个场景均以一定的概率发生，其详细信息如表1所示。另外，在每个场景中，考虑了12个救援周期。在每个周期内，每个节点需求不同。对于物资储备库，每个周期的容量也会因应急场景的不同可能存在差异。网络的总成本包括开设应急设施的建设成本，以及对应每个场景下的物资调配成本。

为了尽可能地符合实际的灾难场景，对于每个应急设施的物资储存过程，以及物资从设施库调往需求点的过程，存在一个物资损坏率（而且物资损坏率在每个场景下每个周期内可能各不相同）。在此基础上，对于每个开放的应急设施，在扣除物资损坏后，调配到各个灾难点的物资总和不能超过该设施每个周期的最大储存量。此外，每个灾难点仅由一个应急设施进行服务。

2 问题分析

已知48个灾难点在可能的10种灾难场景下的物质需求，在8个候选储备库中选着几个进行建立，构建物流网络，考虑运输成本和建设成本，选出储备库向灾难点运输物资的最小成本方案–（类似于指派问题）。

2.1目标

已知48个灾难点在可能的10种灾难场景下的物质需求，在8个候选储备库中选着几个进行建立，构建物流网络，考虑运输成本和建设成本，选出储备库向灾难点运输物资的最小成本方案–（类似于指派问题）。

2.2数据分析

物资需求矩阵

每个灾难点都可能发生10种灾难场景中的任意一种，不同的灾难场景，对应的物资需求量是不同的。灾难发生后，不是储备库向灾难点运输一次物资就ok的，要考虑12个救援周期，每个救援周期的物资需求量是不同的。这样的话，针对每个灾难点，在每个灾难场景，每个救援周期都有一个对应的物资需求，我们可以根据题目所给的数据构建一个物资需求矩阵，维度为（10，12，48）~在灾难场景k，周期l的情况下，灾难场景j的物资需求。

单位运价矩阵

在物资运输过程中，最重要的是运输成本，这也是优化目标之一。根据数据我们可以先提炼出一个距离矩阵（8，48）_{8个储备库候选点到48个灾难点的距离矩阵。题目中提到，12个周期，每个周期的的运输成本是不一样的（单位距离、单位重量）}（12，）。根据这一点，我们将两个矩阵进行合并，得到维度为（12，8，48）的单位重量运价矩阵。

物资完备率矩阵

储备库在不同灾难场景下的是存在物资损耗率的，因为有10个场景嘛，题目就给了10张表，意思是：针对每一个储备库，在不同场景、不同周期下的物资完备率。构建矩阵（10，8，12）~在场景情况下，储备库i在周期l的物资完备率。

2.3已知常量和决策变量

决策变量

• $X_{ij}^{kl}$ ~ 在场景k，周期l的情况下，储备库i向灾难点j运输的物资量
• $Y_i$ ~ i储备库是否开放

已知常量

• I ∈ { 1 , 2 , 3 ， . . . . . . 8 } I \in \{{1,2,3，......8}\} I∈{1,2,3，......8}
• J ∈ { 1 , 2 , 3 ， . . . . . . 48 } J \in \{{1,2,3，......48}\} J∈{1,2,3，......48}
• K ∈ { 1 , 2 , 3 ， . . . . . . 10 } K \in \{{1,2,3，......10}\} K∈{1,2,3，......10}
• L ∈ { 1 , 2 , 3 ， . . . . . . 12 } L \in \{{1,2,3，......12}\} L∈{1,2,3，......12}
• $K_{index}$ – 每个场景出现的概率~$[10]$
• $D_{ij}$ – 单位重量单位距离的距离矩阵 -维度~$[8,48]$
• $C_{ij}^l$ – 单位重量的运价矩阵 ~$[12,8,48]$
• $T_j^l$ – 储备库的物资储备量 ~$[8,12]$
• $J_i^{kl}$ – 物资完备率 ~$[10,8,12]$
• $A_j^{kl}$ – 物资需求量~$[10,8,12]$
• $P$ – 惩罚成本
• $f$ – 建设成本

辅助变量

• $W_{ij}^{kl}$ – 与$X_{ij}^{kl}$同步，为0-1变量，当$X_{ij}^{kl}$大于0时为1，反之亦然。
• ${\delta }_j^{kl}$ – 在场景$k$，周期$l$，灾难点$j$未被满足的物资量
• $Z_{ij}^{kl}$–在场景$k$，周期$l$，储备库能向灾难点$j$运输的物资量

指示变量

• y - 0-1变量

2.4 建立数学模型

目标函数

$$min\sum _{i=1}^8\sum _{j=1}^{48}\sum _{k=1}^{10}\sum _{l=1}^{12}{C}_{ij}^l\cdot X_{ij}^{kl}+\sum _{j=1}^{48}\sum _{k=1}^{10}\sum _{l=1}^{12}P\cdot {\delta }_j^{kl}+\sum _{i=1}^8{Y}_i\cdot f$$

约束条件

约束1：$$\sum _{i=1}^8{W}_{ij}^{kl}=1\forall k\in K\forall l\in L\forall j\in J$$
约束2：$$W_{ij}^{kl}\le Y_i\forall k\in K\forall l\in L\forall i\in I\forall j\in J$$
约束3：$$X_{ij}^{kl}<M\cdot W_{ij}^{kl}\forall k\in K\forall l\in L\forall i\in I\forall j\in J$$
约束4：$$\sum _{j=1}^{48}{X}_{ij}^{kl}\le T_i^l\cdot J_i^{kl}\forall k\in K\forall l\in L\forall i\in I$$
约束5：$$Z_{ij}^{kl}\le T_i^l\cdot J_i^{kl}-\sum _{!=j}^{48}{X}_{ij}^{kl}\forall k\in K\forall l\in L\forall i\in I\forall j\in J$$
约束6：$$\sum _{i=1}^8{Z}_{ij}^{kl}\cdot W_{ij}^{kl}\le A_j^{kl}+M(1-y)\forall k\in K\forall l\in L\forall j\in J$$
约束7：$$A_j^{kl}\le \sum _{i=1}^8{Z}_{ij}^{kl}\cdot W_{ij}^{kl}+M\cdot y\forall k\in K\forall l\in L\forall j\in J$$
约束8：$$\sum _{i=1}^8{X}_{ij}^{kl}=(1-y)A_j^{kl}+y\cdot \sum _{i=1}^8{Z}_{ij}^{kl}\cdot W_{ij}^{kl}\forall k\in K\forall l\in L\forall j\in J$$
约束9：$${\delta }_j^{kl}+\sum _{i=1}^8{X}_{ij}^{kl}=A_j^{kl}\forall k\in K\forall l\in L\forall j\in J$$

模型解释：

• 前三个约束联立起来看，$W_{ij}^{kl}与X_{ij}^{kl}$是一对的，当$X_{ij}^{kl}$是一个大于0的值时，$W_{ij}^{kl}$为1，反之亦然。三个约束联立表示-每个灾难点有且仅有一个储备库进行救援，只有被建立的候选储备库才能被使用
• 第四个约束表示$i$储备库向外输送的物资不能大于其本身的储备量
• 5-8约束联立起来看，第五个约束计算出在场景$k$，周期$l$,储备库$i$向灾难点$j$运输物资时，其当前的剩余物资量，第6~8约束表示为将取出（储备库当前剩余物资量，灾难点物资需求量）的较小一个值，赋值于$X_{ij}^{kl}$
• 约束9为计算出在场景k、周期l、灾难点j无法被满足的物资量。
• 为什么需要加入约束5-8呢，是为了计算出在储备库无法供应的时候的惩罚成本，如果之加入约束9没有约束5-8的话，就需要对惩罚系数进行精心设计，因为当惩罚系数较小的时候，模型可能就不会向灾难点运输物资了，直接赔钱给它得了，加入约束5-8后，模型就会尽全力进行物资供给，在物资不足的时候才考虑惩罚成本。

============================ 手动分割线 ===================================
python调用gurobi对模型进行求解，貌似约束不够紧致，导致搜索空间太大了，求解了好久都没求解出来

import numpy as np
import pandas as pd
from time import time
from scipy.spatial import distance
import torch
import gurobipy as gp
from collections import deque
import random
import matplotlib.pyplot as plt


class SelectAddress(object):
    def __init__(self):
        self.C2_1 = None
        self.C2_2 = None
        self.C2_3 = None
        self.C3 = None
        self.C4 = None
        self.C5 = None
        self.C6 = None
        # random.seed(3)

    def Freight_matrix(self, distance_matrix):
        """

        :param distance_matrix: (8, 48) 距离矩阵
        :return: (12, 8, 48) 特定周期内的价格单位运价矩阵
        """
        Freight = torch.tensor(np.array(self.C4).astype(np.float16)).reshape(12, 1, 1).repeat(1, 8, 48).numpy()
        result = Freight * distance_matrix
        return result

    def Calculate_distance(self):
        """

        :return: (8, 48) 距离矩阵
        """
        X_Y_disaster = self.C2_1.loc[:, ['Latitude', 'Longitude']]
        X_Y_rescue = self.C2_2.loc[:, ['Latitude', 'Longitude']]
        dis = distance.cdist(X_Y_rescue, X_Y_disaster)
        return dis

    # 不做期望
    def first_stage(self):
        """

        :return None
        """
        # 读取数据
        self.read_data()
        # 初始化常量
        i_ = 8
        j_ = 48
        k_ = 10
        l_ = 12
        k_probability = np.array(self.C5.drop(labels='W', axis=1))
        K_delta = np.array(self.C2_3)
        D_i_j = self.Calculate_distance()  # 距离矩阵(8, 48)
        C_l_i_j = self.Freight_matrix(D_i_j)  # 运价 (12, 8, 48)
        T_l_i = np.array(self.C3).T  # 每个周期 每个储备库的库存
        J_k_i_l = np.array(self.C6)  # 每个场景 每个周期 储备库物质的完备率
        A_j_l = np.array(self.C2_1.drop(labels=['State', 'Latitude', 'Longitude'], axis=1))  # 每个周期 每个灾难点的物资需求量
        k_j_l = torch.tensor(np.array(self.C2_3)).reshape(10, 1, 1).repeat(1, 48, 12).numpy()
        A_k_j_l = A_j_l * k_j_l  # 每个场景 每个周期 每个灾难点的物资需求量
        k_index_j = [self.get_random_index() for i in range(48)]
        k_index_i = [self.get_random_index() for i in range(10)]
        f_i = 2e9
        P = np.mean(C_l_i_j)
        M = 10000000

        # 初始化模型
        model = gp.Model('物流网络')

        # 初始化决策变量
        Y = {}
        X = {}
        Z = {}
        D = {}
        W = {}
        y = {}

        for i in range(i_):
            name = f'Y_{i}'
            Y[i] = model.addVar(lb=0, ub=1, vtype=gp.GRB.BINARY, name=name)

        for k in range(k_):
            for l in range(l_):
                for i in range(i_):
                    for j in range(j_):
                        name_X = f'X_{k}_{l}_{i}_{j}'
                        name_Z = f'Z_{k}_{l}_{i}_{j}'
                        name_W = f'W_{k}_{l}_{i}_{j}'
                        X[k, l, i, j] = model.addVar(lb=0, ub=1000000, vtype=gp.GRB.CONTINUOUS, name=name_X)
                        Z[k, l, i, j] = model.addVar(lb=0, ub=1000000, vtype=gp.GRB.CONTINUOUS, name=name_Z)
                        W[k, l, i, j] = model.addVar(lb=0, ub=1, vtype=gp.GRB.BINARY, name=name_W)

        for k in range(k_):
            for l in range(l_):
                for j in range(j_):
                    name_D = f'D_{k}_{l}_{j}'
                    name_y = f'y_{k}_{l}_{j}'
                    D[k, l, j] = model.addVar(lb=0, ub=10000, vtype=gp.GRB.CONTINUOUS, name=name_D)
                    y[k, l, j] = model.addVar(lb=0, ub=1, vtype=gp.GRB.BINARY, name=name_y)

        # 初始化目标函数
        obj = gp.LinExpr(0)
        obj += gp.quicksum(
            (C_l_i_j[l, i, j] * W[k, l, i, j] * A_k_j_l[k, j, l] * k_probability[k]) for k in range(k_) \
            for l in range(l_) for i in range(i_) for j in range(j_))
        obj += gp.quicksum((P * D[k, l, j] * k_probability[k]) for k in range(k_) for l in range(l_) for j in range(j_))
        obj += gp.quicksum((Y[i] * f_i) for i in range(i_))
        model.setObjective(obj, gp.GRB.MINIMIZE)

        # 约束1
        for k in range(k_):
            for l in range(l_):
                for j in range(j_):
                    expr = gp.quicksum(W[k, l, i, j] for i in range(i_))
                    model.addConstr(expr == 1, name=f'Const1_{k}_{l}_{j}')

        # 约束2
        for k in range(k_):
            for l in range(l_):
                for i in range(i_):
                    for j in range(j_):
                        model.addConstr(W[k, l, i, j] <= Y[i], name=f'Const2_{k}_{l}_{i}_{j}')

        # 约束3
        for k in range(k_):
            for l in range(l_):
                for i in range(i_):
                    for j in range(j_):
                        model.addConstr(X[k, l, i, j] <= M * W[k, l, i, j], name=f'Const3_{k}_{l}_{i}_{j}')

        # 约束4
        for k in range(k_):
            for l in range(l_):
                for i in range(i_):
                    expr = gp.quicksum(X[k, l, i, j] for j in range(j_))
                    model.addConstr(expr <= T_l_i[l, i] * J_k_i_l[k, i, l], name=f'Const4_{k}_{l}_{i}')

        # 约束5
        for k in range(k_):
            for l in range(l_):
                for i in range(i_):
                    for j in range(j_):
                        x_delta = gp.quicksum(X[k, l, i, h] for h in range(j_) if h != j)
                        model.addConstr(Z[k, l, i, j] == T_l_i[l, i] * J_k_i_l[k, i, l] - x_delta,
                                        name=f'Const5_{k}_{l}_{i}_{j}')

        # # 约束6 7 8 9
        for k in range(k_):
            for l in range(l_):
                for j in range(j_):
                    X_sum = gp.quicksum(X[k, l, i, j] for i in range(i_))
                    # model.addConstr(X_sum == A_k_j_l[k, j, l] - D[k, l, j], name=f'Const8_{k}_{l}_{j}')
                    delta = gp.quicksum((Z[k, l, i, j] * W[k, l, i, j]) for i in range(i_))
                    model.addConstr(delta <= A_k_j_l[k, j, l] + M * (1 - y[k, l, j]), name=f'Const6_{k}_{l}_{j}')
                    model.addConstr(A_k_j_l[k, j, l] - delta <= M * y[k, l, j], name=f'Const7_{k}_{l}_{j}')
                    if y:
                        model.addConstr(X_sum == delta, name=f'Const8_{k}_{l}_{j}')
                    else:
                        model.addConstr(X_sum == A_k_j_l[k, j, l])
                    model.addConstr(D[k, l, j] + X_sum == A_k_j_l[k, j, l], name=f'Const9_{k}_{l}_{j}')

        # 优化模型
        model.optimize()

        for i in range(i_):
            print(f'{Y[i].VarName} = {Y[i].X}')
        # 将变量的取值写入txt文件
        with open('solution.txt', 'w') as f:
            for v in model.getVars():
                if str(v.varName)[0] == 'W':
                    f.write(f'{v.varName}: {v.x}\n')

        self.Visual_map(0, 0)

    def Stroe_result(self):
        result = deque()
        with open('solution.txt', 'r') as f:
            for line in f:
                _, value = line.strip('\n').split(':')
                result.append(value)
        result_array = np.zeros((10, 12, 8, 48))
        result_array -= 1
        k, l, i, j = result_array.shape
        for k_ in range(k):
            for l_ in range(l):
                for i_ in range(i):
                    for j_ in range(j):
                        result_array[k_, l_, i_, j_] = result.popleft()
        return result_array

    def Visual_map(self, k=None, l=None):
        """

        :return img
        """
        edges = self.Stroe_result()[k, l]  # (8, 48)

        # 储备库 灾难点 位置信息
        X_Y_disaster = np.array(self.C2_1.loc[:, ['Latitude', 'Longitude']])  # (48, 2)
        X_Y_rescue = np.array(self.C2_2.loc[:, ['Latitude', 'Longitude']])  # (8, 2)

        plt.style.use('bmh')
        # 绘制图形
        fig, ax = plt.subplots(figsize=(15, 15))
        for i, (x, y) in enumerate(X_Y_disaster):
            # 绘制灾难点
            ax.plot(x, y, 'o', color=f'C{i % 10}', markersize=6)
            ax.text(x - 0.5, y - 0.5, f'{str(i + 1)}', fontsize=6, color='black', ha='center', va='center')

        for i, (x, y) in enumerate(X_Y_rescue):
            # 绘制救援点
            ax.plot(x, y, '*', color=f'C{i + 10}', markersize=12)
            ax.text(x - 0.5, y - 0.5, f'D{str(i + 1)}', fontsize=6, color='black', ha='center', va='center')

            # 绘制与当前救援点相连的灾难点

            xs = X_Y_disaster[edges[i] > 0, 0]
            ys = X_Y_disaster[edges[i] > 0, 1]
            ax.plot(xs, ys, 'o', color=f'C{i + 10}', markersize=6)

            # 绘制救援点与灾难点之间的连线
            for (xd, yd) in zip(xs, ys):
                ax.plot([x, xd], [y, yd], '--', color=f'C{i + 10}', linewidth=1)

        ax.set(title='Disaster and rescue locations')

        plt.show()

    def read_data(self):
        """

        :return None
        """
        print('=======================<开始读取数据>=======================')
        T1 = time()
        # 获取48个灾难点的信息
        self.C2_1 = pd.read_csv('data/C2_1.csv')
        # 获取8个候选点的信息
        self.C2_2 = pd.read_csv('data/C2_2.csv')
        # self.C2_2.loc[0, ['Latitude', 'Longitude']] = [40, -100]
        # 获取场景系数
        self.C2_3 = pd.read_csv('data/C2_3.csv').drop(labels='场景', axis=1)

        # C3 每个候选点每个时期的资源
        self.C3 = pd.read_excel('data/教指委案例大赛-Data.xlsx', sheet_name='C3')[2:].drop(labels='Unnamed: 0', axis=1)
        self.C3.columns = [i for i in range(12)]
        # self.C3 = self.C3.rename(columns={0: 'State'})

        # C4 Transportation cost from emergency facilities to demand points in period t  ($/(km*kg) )
        self.C4 = pd.read_excel('data/教指委案例大赛-Data.xlsx', sheet_name='C4')[1:]
        self.C4.columns = [i for i in range(12)]

        # C5 Table C5. 每个场景的概率
        self.C5 = pd.read_excel('data/教指委案例大赛-Data.xlsx', sheet_name='C5')[2:].loc[:, ['Unnamed: 0', 'Unnamed: 1']]
        self.C5.columns = ['W', 'P']

        # 获取C6~  物资完好率
        self.C6 = np.zeros((10, 8, 12))
        for i in range(1, 11):
            if i == 1:
                temp = pd.read_excel('data/教指委案例大赛-Data.xlsx', sheet_name=f'C6-Scenario {i}')[49:57].drop(
                    labels='Table C11. 物资完好率', axis=1)
            else:
                temp = pd.read_excel('data/教指委案例大赛-Data.xlsx', sheet_name=f'C6-Scenario {i}')[48:57].drop(
                    labels='Unnamed: 0', axis=1)
            self.C6[i - 1, :, :] = temp
        T2 = time()
        print(f'数据读取时间为： {(T2 - T1):04f}')
        print('=======================<数据读取完毕>=======================')

    def get_random_index(self):
        # 根据给定的选项和概率权重进行随机选择
        options = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]  # 定义选项
        weights = [0.1, 0.122, 0.026, 0.011, 0.09, 0.12, 0.098, 0.17, 0.234, 0.029]  # 定义概率权重
        random_index = random.choices(options, weights=weights)[0]
        return random_index


if __name__ == '__main__':
    Model = SelectAddress()
    Model.first_stage()

结果可视化

求解的gap始终无法到达0，且当前最优解几乎不会下降了，就把程序停了下了，出现了两种结果

• 开6库（星星标号为储备库的位置，圆圈为在难点位置）>开7库（星星标号为储备库的位置，圆圈为在难点位置）
  

在此之前我建立过一个比较简单的模型求解这个问题，但是为了应付答辩特意把模型改复杂了，结果貌似没有之前的好了，后续有时间再找一些智能算法求解一下吧，摆烂了摆烂了/(ㄒoㄒ)/~~