基于三次样条插值和单纯形法的加氢站选址优化

问题描述

已知定点交通流量，求解加氢站建设位置

求解方法

已知国道或省道定点交通流量若干，根据已知交通流量插值得到每3km对应的交通流量。如图所示。

将该问题转换为p-中值问题：

其中，
需求点位置集合=每3km一个需求点

在i点的客户人数=i点（每个需求点）的车流量

设施总数=需要建设的加氢站数量

i点与j点之间的距离由第4步构建

整体模型思路参考https://blog.youkuaiyun.com/m0_45008035/article/details/127017479，

将该问题转换到加氢站选址，最终实现整个山西省吕梁市的加氢站规划建设

1. 加氢站数量选定和数据集读取

loc_nums = 6
sheet_name = 'S223'
df_1 = pd.read_excel('路段平均日交通量(1).xlsx',sheet_name=sheet_name)

2. 累计公里数计算

def re_cal_km(df_car):
    if len(df_car)>1:
        df_car = df_car.drop([0, 1, 2, len(df_1)-1])
        df_car = df_car.reset_index(drop=True)
    df_car = df_car[['汽车（辆/日）','观测里程（公里）']]
    if len(df_car)>1:
        df_car.at[0,'累计公里'] = df_car.at[0,'观测里程（公里）']
        df_car.at[1,'累计公里'] = df_car.at[1,'观测里程（公里）'] + df_car.at[0,'观测里程（公里）']
        for i in range(2,len(df_car)):
            df_car.at[i,'累计公里'] = df_car.at[i,'观测里程（公里）'] + df_car.at[i-1,'累计公里']
        df_car = df_car.drop('观测里程（公里）',axis=1)
        df_car.at[0,'累计公里'] = 0
    else:
        df_car.rename(columns={'观测里程（公里）': '累计公里'}, inplace=True)
        total_km = df_car['累计公里'].values[0]
        medin_km = total_km/2
        medin_cars = df_car['汽车（辆/日）'].values[0]
        r_num = np.random.randint(0, 10)
        cars_start = medin_cars - medin_km*3*r_num
        cars_end = medin_cars + medin_km*3*r_num
        df_car.at[0,'累计公里'] = 0
        df_car.at[1,'累计公里'] = total_km
        df_car.at[0,'汽车（辆/日）'] = cars_start
        df_car.at[1,'汽车（辆/日）'] = cars_end
        df_car = df_car.sort_index()
    return df_car

3. 交通流量插值

from scipy.interpolate import interp1d
import numpy as np

def inter_traffic(df_car,loc_nums):
    # 使用累计公里作为X轴，汽车（辆/日）作为Y轴
    x = df_car['累计公里']
    y = df_car['汽车（辆/日）']

    # 插值函数
    if len(x) <= 2:
        interp_function = interp1d(x, y, fill_value=0, bounds_error=False)
    else:    
        interp_function = interp1d(x, y, kind='cubic', fill_value=0, bounds_error=False)

    max_km = df_car.tail(1)['累计公里'].values[0]
    dist_col = int(max_km/loc_nums)
    if (loc_nums * 6) > max_km:
        if dist_col == 1:
            dist_col = 1.5
        print('公里数过小',max_km,loc_nums,dist_col)
        new_x = np.arange(0, x.max(),dist_col)
    else:
        new_x = np.arange(0, x.max(), 6)
        dist_col = 6
    # 使用插值函数计算新的Y轴数据
    new_y = interp_function(new_x)

    # 创建新的DataFrame以查看结果
    new_df = pd.DataFrame({'累计公里': new_x, '汽车（辆/日）': new_y})
    return new_df,dist_col

4. 距离矩阵构建

def cal_dist(new_df):
    # 计算两两之间的距离
    num_points = len(new_df)
    dist_matrix = {}

    for i in range(num_points):
        for j in range(num_points):
            # 计算绝对差值
            dist = abs(new_df['累计公里'].iloc[i] - new_df['累计公里'].iloc[j])
            # 存储在字典中
            dist_matrix[(i, j)] = dist
    return dist_matrix

5. 优化选址

from gurobipy import *
import numpy as np
import random

def gen_result(new_df,loc_nums,dist_matrix,dist_col):
    # 参数
    num_points = len(new_df)
    random.seed(0)
    np.random.seed(0)

    # 假设有预定义的距离矩阵 dist_matrix,
    num_people = new_df['汽车（辆/日）'].values

    # 创建模型
    m = Model("p-median Problem")
    m.setParam('MIPGap', 0.01)  # 设置1%的最优间隙
    m.setParam('MIPFocus', 2)  # 优先提高解的质量

    # 创建变量
    select = m.addVars(num_points, vtype=GRB.BINARY, name='Select') # X
    assign = m.addVars(num_points, num_points, vtype=GRB.BINARY, name='Assign') # Y

    # 添加约束
    m.addConstr((select.sum() == loc_nums), name='Num_limit')
    m.addConstrs((assign[i, j] <= select[j] for i in range(num_points) for j in range(num_points)), name='Assign_before_locate')
    m.addConstrs((assign.sum(i, '*') == 1 for i in range(num_points)), name='Unique_assign')

    # 设置目标
    m.setObjective(sum(num_people[i] * dist_matrix[i, j] * assign[i, j] for i in range(num_points) for j in range(num_points)), GRB.MINIMIZE)

    # 优化模型
    m.optimize()

    # 打印结果
    selected = [i for i in range(num_points) if select[i].X > 0.5]
    assigned = [(i, j) for i in range(num_points) for j in range(num_points) if assign[i, j].X > 0.5]
    print("\n选址间距(km) = "+str(dist_col)+" \nSelected positions ="+ str(selected))
    # print("Assigned relationships =", assigned)
    # print('Objective value = %g' % m.objVal)
    return dist_col,selected

运行结果