优化 | 考虑随机性的车辆路径规划问题详解（Stochastic VRP）：模型、理论推导及Python调用Gurobi实现)

作者：樵溪子，清华大学，清华大学深圳国际研究生院，清华-伯克利深圳学院，硕士在读

审校：刘兴禄，清华大学，清华大学深圳国际研究生院，清华-伯克利深圳学院，博士在读

0. 摘要

本文首先介绍了带时间窗的车辆路径规划问题(VRPTW)的问题描述、模型建立，并用python调用Gurobi进行求解。之后简要概述了考虑随机性的VRP问题(SVRPTW)的问题类型和建模方式。主要针对考虑需求随机性的VRPTW问题(VRPTWSD)，以机遇约束规划 (Chance-constrainted programming) 的方式建模，通过理论推导将不确定性约束转化为确定性约束，然后用python调用Gurobi进行求解。

本文的模型参考自文献：

[1] Desaulniers, Guy, Jacques Desrosiers, and Marius M. Solomon, eds. Column generation. Vol. 5. Springer Science & Business Media, 2006.
[2] Stewart Jr, William R., and Bruce L. Golden. “Stochastic vehicle routing: A comprehensive approach.” European Journal of Operational Research 14.4 (1983): 371-385.

1. 带时间窗的车辆路径规划问题(VRPTW)

1.1 问题描述

VRPTW问题定义在一队车辆$\mathcal{V}$、一组客户$\mathcal{C}$和一个有向图$\mathcal{G}$上。图$\mathcal{G}$由$|\mathcal{C}|+2$个点组成，其中客户由$1,2,...,n$表示，车站由$0$（始发站）和$n+1$（终点站）表示。所有节点的集合可表示为 N = { 0 , 1 , . . . , n + 1 } \mathcal{N}=\{0,1,...,n+1\} N={0,1,...,n+1}。所有边的集合$\mathcal{A}$表示了车站与客户之间，以及客户之间的有向连接。其中没有边从$0$点终止或者从$n+1$点开始。每一条弧$(i,j),i\ne j$都对应一个成本$c_{ij}$和时间$t_{ij}$，时间可以包括在弧$(i,j)$上的行驶时间和在$i$点的服务时间。所有车辆通常是同质化的，每辆车都存在容量上限$Q$，每一个客户$i$点都有需要被满足的需求$d_i$，并且需要在时间窗$[a_i,b_i]$之内被服务[1]。待优化的问题即为，如何决策车辆访问客户的路径，使得在满足一定约束的条件下，实现最小化总成本的目标。

1.2 模型构建

（1） 参数

• $\mathcal{V}$：车辆集合， V = { 0 , 1 , . . . , V } \mathcal{V}=\{0,1,...,V\} V={0,1,...,V}；
• $\mathcal{C}$：客户点的集合， C = { 1 , 2 , . . . , n } \mathcal{C}=\{1,2,...,n\} C={1,2,...,n}；
• $\mathcal{N}$：所有节点的集合， N = { 0 , 1 , . . . , n + 1 } \mathcal{N}=\{0,1,...,n+1\} N={0,1,...,n+1}；
• $c_{ij}$：从$i$到$j$的行驶距离；
• $d_i$：$i$点的客户需求量；
• $Q$：车辆容量；
• $t_{ij}$：从 $i$点到 $j$点的行驶时间，以及服务$i$点的时间之和；
• $[a_i,b_i]$：访问 $i$点的时间窗；
• $M_{ij}$：足够大的正数。

（2） 决策变量

• $s_{ik}$：车辆$k$服务$i$点的开始时刻；
• $x_{ijk}$：车辆路径决策变量
  $$x_{ijk}=\{\begin{array}{ll}1,& 车辆k从弧(i,j)上经过\\ 0,& 其他\end{array}$$

（3）混合整数规划模型
基于上述决策变量和参数，VRPTW可以被建模为下面的混合整数规划模型。

min ⁡ ∑ k ∈ V ∑ i ∈ N ∑ j ∈ N c i j x i j k s . t . ∑ i ∈ C d i ∑ j ∈ N x i j k ⩽ Q ∀ k ∈ V , ∑ k ∈ V ∑ j ∈ N x i j k = 1 ∀ i ∈ C , ∑ j ∈ N x 0 j k = 1 ∀ k ∈ V , ∑ i ∈ N x i h k − ∑ j ∈ N x h j k = 0 ∀ h ∈ C , ∀ k ∈ V , ∑ i ∈ N x i , n + 1 , k = 1 ∀ k ∈ V , s i k + t i j − M i j ( 1 − x i j k ) ⩽ s j k , ∀ i , j ∈ N , ∀ k ∈ V , a i ⩽ s i k ⩽ b i ∀ i ∈ N , ∀ k ∈ V , x i j k ∈ { 0 , 1 } ∀ i , j ∈ N , ∀ k ∈ V s i k ≥ 0 ∀ i ∈ N , ∀ k ∈ V \begin{align} \min \quad &\sum_{k \in \mathcal{V}} \sum_{i \in \mathcal{N}}\sum_{j \in \mathcal{N}} c_{ij}x_{ijk} && \\ s.t. \quad &\sum_{i \in \mathcal{C}}d_i \sum_{j \in \mathcal{N}} x_{ijk} \leqslant Q \quad && \forall k \in \mathcal{V}, \\ & \sum_{k \in \mathcal{V}} \sum_{j \in \mathcal{N}} x_{ijk} =1 \quad &&\forall i \in \mathcal{C}, \\ &\sum_{j \in \mathcal{N}}x_{0jk}=1 \quad && \forall k \in \mathcal{V},\\ & \sum_{i \in \mathcal{N}} x_{ihk}-\sum_{j \in \mathcal{N}}x_{hjk}=0 \quad && \forall h \in \mathcal{C}, \forall k \in \mathcal{V},\\ & \sum_{i \in \mathcal{N}}x_{i,n+1,k}=1 \quad \forall k \in \mathcal{V},\\ %&x_{ijk}(s_{ik}+t_{ij}-s_{jk}) \leqslant 0 \quad &&\forall i,j \in \mathcal{N}, \forall k \in \mathcal{V},\\ & s_{ik}+t_{ij}-M_{ij}(1-x_{ijk}) \leqslant s_{jk}, \quad&& \forall i,j \in \mathcal{N}, \forall k \in \mathcal{V}, \\ &a_i \leqslant s_{ik} \leqslant b_i \quad && \forall i \in \mathcal{N}, \forall k \in \mathcal{V},\\ & x_{ijk} \in \{0,1\} \quad && \forall i,j \in \mathcal{N}, \forall k \in \mathcal{V}\\ & s_{ik} \geq 0\quad && \forall i \in \mathcal{N}, \forall k \in \mathcal{V} \end{align} mins.t.​k∈V∑​i∈N∑​j∈N∑​cij​xijk​i∈C∑​di​j∈N∑​xijk​⩽Qk∈V∑​j∈N∑​xijk​=1j∈N∑​x0jk​=1i∈N∑​xihk​−j∈N∑​xhjk​=0i∈N∑​xi,n+1,k​=1∀k∈V,sik​+tij​−Mij​(1−xijk​)⩽sjk​,ai​⩽sik​⩽bi​xijk​∈{0,1}sik​≥0​​∀k∈V,∀i∈C,∀k∈V,∀h∈C,∀k∈V,∀i,j∈N,∀k∈V,∀i∈N,∀k∈V,∀i,j∈N,∀k∈V∀i∈N,∀k∈V​​

下面是各个约束的具体含义：

• 目标函数(1)最小化总体行驶距离；
• 约束(2)说明车辆的载重量不能超过其容量上限；
• 约束(3)保证了每个客户点都被访问了一次；
• 约束(4-6)分别保证了每辆车必须从始发站出发，到达并离开每个客户点，并最终回到终点站；
• 约束(7)规定了车辆$k$开始服务$j$点的时刻$s_{jk}$，不能早于开始服务$i$点的时刻$s_{ik}$，加上从$i$点到 $j$点的行驶时间以及服务$i$点的时间之和$t_{ij}$，其中$M_{ij}$的取值下界为$b_i+t_{ij}-a_j$；
• 约束(8)使得开始服务$i$点的时刻是在规定的时间窗范围内；
• 约束(9-10)是对于决策变量的约束。

1.3 数值实验

本文使用的算例为Solomn算例中的c101.txt ，利用前25个客户点作为数据集并调用Gurobi求解。

• 求解日志

Gurobi Optimizer version 9.5.2 build v9.5.2rc0 (mac64[arm])
Thread count: 8 physical cores, 8 logical processors, using up to 8 threads
Optimize a model with 5513 rows, 5103 columns and 33684 nonzeros
Model fingerprint: 0x0167fefc
Variable types: 189 continuous, 4914 integer (4914 binary)
Coefficient statistics:
  Matrix range     [1e+00, 1e+04]
  Objective range  [1e+00, 4e+01]
  Bounds range     [1e+00, 1e+03]
  RHS range        [1e+00, 1e+04]
Presolve removed 4256 rows and 1953 columns
Presolve time: 0.04s
Presolved: 1257 rows, 3150 columns, 20472 nonzeros
Variable types: 189 continuous, 2961 integer (2961 binary)

Root relaxation: objective 1.840273e+02, 468 iterations, 0.01 seconds (0.01 work units)

    Nodes    |    Current Node    |     Objective Bounds      |     Work
 Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time

     0     0  184.02728    0   60          -  184.02728      -     -    0s
     0     0  184.02728    0   36          -  184.02728      -     -    0s
     0     0  184.02728    0   36          -  184.02728      -     -    0s
     0     0  184.02728    0   31          -  184.02728      -     -    0s
     0     0  184.02728    0   31          -  184.02728      -     -    0s
H    0     0                     193.0414573  184.02728  4.67%     -    0s
     0     0  184.02728    0   28  193.04146  184.02728  4.67%     -    0s
     0     0  184.02728    0   28  193.04146  184.02728  4.67%     -    0s
     0     0  184.04305    0   31  193.04146  184.04305  4.66%     -    0s
     0     0  184.04305    0   15  193.04146  184.04305  4.66%     -    0s
H    0     0                     192.6417404  184.04305  4.46%     -    0s
     0     0  184.04305    0   20  192.64174  184.04305  4.46%     -    0s
     0     0  184.04305    0   15  192.64174  184.04305  4.46%     -    0s
H    0     0                     191.8136198  184.04305  4.05%     -    0s
     0     0     cutoff    0       191.81362  191.81362  0.00%     -    0s

Cutting planes:
  Learned: 1
  Gomory: 5
  Implied bound: 3
  Clique: 14
  MIR: 4
  RLT: 1
  Relax-and-lift: 1

Explored 1 nodes (1814 simplex iterations) in 0.26 seconds (0.35 work units)
Thread count was 8 (of 8 available processors)

Solution count 3: 191.814 192.642 193.041 

Optimal solution found (tolerance 1.00e-04)
Best objective 1.918136197787e+02, best bound 1.918136197787e+02, gap 0.0000%

• 结果输出

Model status   : 2    
Objective value: 191.81
Vehicle capacity: 200.00

Vehicle id | Route 
1          | [0, 5, 3, 7, 8, 10, 11, 9, 6, 4, 2, 1, 26] 
4          |[0, 13, 17, 18, 19, 15, 16, 14, 12, 26] 
6          | [0, 20, 24, 25, 23, 22, 21, 26] 


Vehicle id |  Route length   Travel distance    Loads           
1          |  13             59.49              160.00          
4          |  10             95.88              190.00       
6          |  8              36.44              110.00    

• 结果可视化
  

2. 考虑随机性的车辆路径规划问题简介 (SVPR)

2.1 问题分类

考虑随机性的车辆路径规划问题，从随机性的内容出发可以分为多种类型[3]，例如：

• （1）随机需求(VRPSD)：客户的需求量是随机变量；
• （2）随机客户(VRPSC)：客户以一定概率出现，但是有确定的需求量；
• （3）随机客户和需求(VRPSCD)：客户以一定概率出现，且需求量也随机；
• （4）随机行驶时间(VRPTWSTT)：在两点之间的行驶时间，以及服务时间是随机变量。

2.2 建模方式

SVPR问题的建模方法主要有以下几种：

• （1）机会约束规划，Chance constrained program (CCP)
• （2）带补偿的随机规划，Stochastic program with recourse (SPR)
• （3）动态规划，Dynamic programming (DP)
• （4）鲁棒优化，Robust optimization (RO)

本文主要针对考虑需求随机性的VRPTW问题，以机遇约束规划 (Chance-constrainted programming) 的方式建模，通过理论推导将不确定性约束转化为确定性约束，然后调用Gurobi求解。

3. 考虑需求不确定性的VRPTW：数学模型、模型转化和代码实现

上文提到了多种参数的不确定性，例如需求不确定、客户不确定等，但是不同的不确定因素的处理方法有一定的相似之处。本小节就以需求量不确定为例来介绍考虑随机因素的VRP的建模和模型重构方法。此外，为了帮助读者能加深理解，本文提供了完整的实现代码。

3.1 数学模型

Chance-constrainted programming最早由Charnes 和Cooper提出[4,5]，该模型将VRPTW中的车辆容量确定性约束(2) 转化为机会约束(13)。

$$\begin{array}{rlrl}\min & \sum _{k\in \mathcal{V}}\sum _{i\in \mathcal{N}}\sum _{j\in \mathcal{N}}{c}_{ij}{x}_{ijk}& & \\ s.t.& (3)-(10),& & \\ & P[\sum _{i,j}{d}_i{x}_{ijk}⩽Q]⩾1-\alpha ,& & \forall k\in \mathcal{V}.\end{array}$$

其中$P[{\sum }_{i\in \mathcal{C}}{d}_i{\sum }_{j\in \mathcal{N}}{x}_{ijk}⩽Q]$代表了车辆载重未超过容量上限的概率，$\alpha$代表了该约束被违背的最大接受概率。$d_i$此处是一组独立分布的随机变量，代表了每个客户点的随机需求量。

注意，不同的$d_i$也许会服从不同的分布，但是，一些模型重构的方法对$d_i$的分布一般会有一定的要求，所以这里需要具体问题具体分析。

3.2 理论推导：模型转化

在推导模型之前，我们先给出模型推导的大体思路。一般情况下，处理机会约束的方法，就是通过模型转化，将带有概率表达式的约束转化为确定性的不等式，从而使模型可以直接求解。下面，就来具体介绍本文考虑的随机约束的转化方法。

在上述模型中，约束(13)可以被转化为确定性约束。对于所有客户$\forall i\in \mathcal{C}$而言，若其需求$d_i$是一组独立分布的随机变量，其需求的均值为${\mu }_i$，标准差为${\sigma }_i$。基于此，对于车辆$k$而言，其行驶路径上访问的客户的总需求的均值和标准差可以用下面的表达式来计算：

总需求的均值
$$\begin{array}{r}{M}_k=\sum _{i,j}{\mu }_i{x}_{ijk}\end{array}$$

总需求的标准差
$$\begin{array}{r}{S}_k=\sqrt{\sum _{i,j}{\sigma }_i^2{x}_{ijk}^2}\end{array}$$

若存在常数$\tau$使得
$$\begin{array}{r}P\left[\right(\sum _{i,j}{d}_i{x}_{ijk}-M_k)/S_k⩽\tau ]=1-\alpha ,\end{array}$$
则约束(16)可以被表示为如下确定性约束
$$\begin{array}{r}{M}_k+\tau S_k⩽Q,\forall k\in \mathcal{V}.\end{array}$$
若写为展开形式，可知约束(17)是非线性约束：
$$\begin{array}{r}\sum _{i,j}{\mu }_i{x}_{ijk}+\tau \sqrt{\sum _{i,j}{\sigma }_i^2{x}_{ijk}^2}⩽Q,\forall k\in \mathcal{V}.\end{array}$$
当需求的分布满足一定条件时，非线性约束(16)可以转化为关于${\mu }_i$的线性约束[2]。

定理 \quad 如果需求$d_i$满足：

• （1）需求$d_i$的概率分布是独立的；
• （2）$({\sum }_{i,j}{d}_i{x}_{ijk}-M_k)/S_k$和$(d_i-{\mu }_i)/{\sigma }_i$拥有相同的分布，且对所有$i$成立；
• （3）${\sigma }_i^2$是${\mu }_i$的常数倍，即${\sigma }_i^2=\pi {\mu }_i$，
  则存在常数$\overline{Q}$使得线性约束
  $$\begin{array}{r}\sum _{i,j}{\mu }_i{x}_{ijk}⩽\overline{Q},\forall k\in \mathcal{V}.\end{array}$$
  和机会约束(13)和(14)，以及确定性约束(15)和(16)是等价的。

证明 \quad VRP和SVRP问题中的决策变量$x_{ijk}$是0-1决策变量，则有$(x_{ijk}{)}^2=x_{ijk}$。因此，
$$\begin{array}{r}{S}_k=\sqrt{\sum _{i,j}{\sigma }_i^2{x}_{ijk}^2}=\sqrt{\sum _{i,j}\pi {\mu }_i{x}_{ijk}}=\sqrt{\pi \sum _{i,j}{\mu }_i{x}_{ijk}}=\sqrt{\pi M_k}.\end{array}$$
将(20)带入(17)则有

$$\begin{array}{r}{M}_k+\tau \sqrt{\pi M_k}⩽Q.\end{array}$$
进而得到
$$\begin{array}{r}{M}_k⩽\frac{2Q+{\tau }^2\pi -\sqrt{{\tau }^4{\pi }^2+4Q{\tau }^2\pi }}{2}=\overline{Q}.\end{array}$$
即
$$\begin{array}{r}\sum _{i,j}{\mu }_i{x}_{ijk}⩽\overline{Q}.\end{array}$$
多种分布均满足上述定理的要求，例如泊松分布、正态分布、二项分布、伽马分布等。
因此，若需求的分布满足上述定理的要求，则转化后的模型是如下容量约束为$\overline{Q}$的VRPTW问题：
$$\begin{array}{rlrl}\min & \sum _{k\in \mathcal{V}}\sum _{i\in \mathcal{N}}\sum _{j\in \mathcal{N}}{c}_{ij}{x}_{ijk}& & \\ s.t.& (3)-(10),& & \\ & \sum _{i,j}{\mu }_i{x}_{ijk}⩽\overline{Q},& & \forall k\in \mathcal{V}.\end{array}$$

3.3 Python调用Gurobi求解

为了方便读者加深理解，本小节使用Python调用Gurobi求解上文介绍的数学规划模型，并给出求解结果的可视化。

我们使用Solomn VRPTW benchmark instance算例中的c101.txt，利用前25个客户点作为数据集。假设客户需求是一组独立的随机变量$d_0,d_1,...,d_n$，其中每个$d_i$都服从正态分布，均值${\mu }_i$等于原始数据中的需求量${\overline{d}}_i$，方差${\sigma }_i^2=\pi {\mu }_i$ ($\pi =0.5$)。

• 主体代码（完整代码获取方式见文末）

class ModelBuilder(object):
    def __init__(self, data):
        self.data = data
        self.x = {}
        self.s = {}
        self.v = 5         # velocity of each vehicle
        self.M = 10000     # a sufficiently large number
        self.alpha = 0.1   #  the maximum allowable probability of route failure 
        self.mu = []       # mean
        self.sigma = []    # standard deviation 
    
    def stochastic_demand(self, pi):
        """ demand at each customer: normal distribution  """

        # mean demand = demand from data®
        self.mu = copy.deepcopy(self.data.demand)

        # std = sqrt(pi * mu)
        for i in range(len(self.mu)):
            self.sigma.append(math.sqrt(pi * self.mu[i]))

        # update demand value = random variable value
        for i in range(len(self.data.demand)):
            self.data.demand[i] = np.random.normal(self.mu[i],self.sigma[i],1)[0]

        # chance constraint: capacity
        r = norm.ppf(1-self.alpha, 0, 1)
        Q = self.data.vehicle_capacity
        self.Q = (2 * Q + r**2 * pi - math.sqrt(r**4 * pi**2 + 4 * Q * r**2 * pi))/2


    def build_model(self, data=None, stochastic=False, pi=0, result_path=''):

        self.model = None
        self.model = Model('SVRPTW')

        """ create the decision variables """
        for i in self.data.node_id:
            for k in range(self.data.vehicle_number):
                self.s[i,k] = self.model.addVar(lb=0, ub=self.data.due_time[-1], vtype=GRB.CONTINUOUS, name='s_'+str(i)+'_'+str(k))
                for j in self.data.node_id:
                    if ( i!= j ):
                        self.x[i,j,k] = self.model.addVar(lb=0, ub=1, vtype=GRB.BINARY, name='x_'+str(i)+'_'+str(j)+'_'+str(k))

        self.model.update()

        """ (1) set the objective function """
        obj = LinExpr()
        for i in self.data.node_id:
            for j in self.data.node_id:
                for k in range(self.data.vehicle_number):
                    if ( i!= j ):
                        c = self.data.distance_matrix[i][j]
                        obj.addTerms(c, self.x[i,j,k])
        self.model.setObjective(obj, GRB.MINIMIZE)    

        
        """ (2) add constraints : vehicle capacity limit (stochastic or not)  """
        if not stochastic:
            self.Q = self.data.vehicle_capacity
            for k in range(self.data.vehicle_number):
                lhs = LinExpr()
                for i in self.data.customer_id:
                    d = self.data.demand[i]
                    for j in self.data.node_id:
                        if ( i!= j ):
                            lhs.addTerms(d,self.x[i,j,k])
                self.model.addConstr(lhs <= self.Q, name='vehicle_capacity_'+str(k))
        else:
            self.stochastic_demand(pi)    # including update self.data.demand
            for k in range(self.data.vehicle_number):
                lhs = LinExpr()
                for i in self.data.customer_id:
                    for j in self.data.node_id:
                        if ( i!= j ):
                            lhs.addTerms(self.mu[i],self.x[i,j,k]) 
                self.model.addConstr(lhs <= self.Q, name='chance_constarint_'+str(k))


        """ (3) add constraints : each customer is visited """
        for i in self.data.customer_id: 
            lhs = LinExpr()
            for j in self.data.node_id:
                if ( i!= j ):
                    for k in range(self.data.vehicle_number):
                        lhs.addTerms(1,self.x[i,j,k])
            self.model.addConstr(lhs == 1, name='visit_customer_'+str(i))

        """ (4) add constraints : vehicle leave starting depot """
        for k in range(self.data.vehicle_number):
            lhs = LinExpr()
            for j in self.data.node_id:
                if ( j!= 0 ):
                    lhs.addTerms(1,self.x[0,j,k])
            self.model.addConstr(lhs == 1 , name='leave_depot'+str(k))

        """ (5) add constraints : flow conservation - vehicle visit customer then leave"""
        for k in range(self.data.vehicle_number):
            for h in self.data.customer_id:
                lhs1 = LinExpr()
                for i in self.data.node_id:
                    if ( h!= i ):
                        lhs1.addTerms(1,self.x[i,h,k])
                lhs2 = LinExpr()
                for j in self.data.node_id:
                    if ( h!= j ):
                        lhs2.addTerms(1,self.x[h,j,k])
                self.model.addConstr(lhs1 == lhs2 , name='flow_conservation'+str(k)+'_'+str(h))

        """ (6) add constraints : vehicle return to returning depot """
        for k in range(self.data.vehicle_number):
            lhs = LinExpr()
            for i in self.data.node_id:
                if ( i!= self.data.node_id[-1] ):
                    lhs.addTerms(1,self.x[i,self.data.node_id[-1],k])
            self.model.addConstr(lhs == 1 , name='return_depot'+str(k))

        """ (7) add constraints : vehivle departure time relationship  """
        for k in range(self.data.vehicle_number):
            for i in self.data.node_id:
                for j in self.data.node_id:
                    if ( i!= j ):
                        t = self.data.distance_matrix[i][j] / self.v
                        self.model.addConstr(self.s[i,k] + t - self.M * (1-self.x[i,j,k]) <= self.s[j,k], name='departure_time'+str(i)+'_'+str(j)+'_'+str(k))

        """ (8) add constraints : node time window  """
        for k in range(self.data.vehicle_number):
            for i in self.data.node_id:
                self.model.addConstr(self.s[i,k] >= self.data.ready_time[i], name='ready_time'+str(i)+'_'+str(k))
                self.model.addConstr(self.s[i,k] <= self.data.due_time[i], name='due_time'+str(i)+'_'+str(k))

        self.model.update()

        """ get the optimal solution """
        self.model.optimize()
        model_result = self.model.ObjVal, self.x, self.model.status, self.Q

        """ store the model file """
        self.data.customer_number = len(self.data.customer_id)
        if stochastic:
            type = 'SVRPTW'
            self.model.write(result_path + type + '_' + str(self.data.customer_number) + '_pi=' + str(pi) + '.lp')
            self.model.write(result_path + type + '_' + str(self.data.customer_number) + '_pi=' + str(pi) + '.rlp')
        else:
            type = 'VRPTW'
            self.model.write(result_path +type+'_'+str(self.data.customer_number)+'.lp')
            self.model.write(result_path +type+'_'+str(self.data.customer_number)+'.rlp')
        return model_result

• 求解日志

Gurobi Optimizer version 9.5.2 build v9.5.2rc0 (mac64[arm])
Thread count: 8 physical cores, 8 logical processors, using up to 8 threads
Optimize a model with 5513 rows, 5103 columns and 33684 nonzeros
Model fingerprint: 0x511b6a13
Variable types: 189 continuous, 4914 integer (4914 binary)
Coefficient statistics:
  Matrix range     [1e+00, 1e+04]
  Objective range  [1e+00, 4e+01]
  Bounds range     [1e+00, 1e+03]
  RHS range        [1e+00, 1e+04]
Presolve removed 4256 rows and 1953 columns
Presolve time: 0.04s
Presolved: 1257 rows, 3150 columns, 20472 nonzeros
Variable types: 189 continuous, 2961 integer (2961 binary)

Root relaxation: objective 1.840273e+02, 464 iterations, 0.01 seconds (0.01 work units)

    Nodes    |    Current Node    |     Objective Bounds      |     Work
 Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time

     0     0  184.02728    0   60          -  184.02728      -     -    0s
     0     0  184.02728    0   55          -  184.02728      -     -    0s
     0     0  184.02728    0   55          -  184.02728      -     -    0s
H    0     0                     262.3027177  184.02728  29.8%     -    0s
H    0     0                     255.4779601  184.02728  28.0%     -    0s
     0     0  184.02728    0   60  255.47796  184.02728  28.0%     -    0s
     0     0  184.02728    0   60  255.47796  184.02728  28.0%     -    0s
     0     0  184.02728    0   57  255.47796  184.02728  28.0%     -    0s
     0     0  184.02728    0   59  255.47796  184.02728  28.0%     -    0s
     0     0  184.02728    0   57  255.47796  184.02728  28.0%     -    0s
H    0     0                     253.0177705  184.02728  27.3%     -    0s
     0     0  184.02728    0   57  253.01777  184.02728  27.3%     -    0s
     0     0  189.28589    0   61  253.01777  189.28589  25.2%     -    0s
H    0     0                     251.7899329  189.28589  24.8%     -    0s
     0     0  189.28589    0   61  251.78993  189.28589  24.8%     -    0s
     0     0  189.28589    0   61  251.78993  189.28589  24.8%     -    0s
     0     0  189.50286    0   61  251.78993  189.50286  24.7%     -    0s
     0     0  191.60105    0   56  251.78993  191.60105  23.9%     -    0s
     0     0  191.60105    0   56  251.78993  191.60105  23.9%     -    0s
     0     0  191.60105    0   56  251.78993  191.60105  23.9%     -    0s
H    0     2                     234.3355243  191.60105  18.2%     -    0s
     0     2  191.60105    0   56  234.33552  191.60105  18.2%     -    0s
*   79    62               9     231.1985824  191.81362  17.0%  54.0    0s
H   80    62                     229.6371690  191.81362  16.5%  53.3    0s
H   90    64                     228.1675465  191.81362  15.9%  53.8    0s
H  373   161                     228.1260451  191.81362  15.9%  34.6    0s
*  667   275              18     228.0484854  191.81362  15.9%  28.9    1s
*  688   275              44     227.1573287  191.81362  15.6%  28.3    1s

Cutting planes:
  Learned: 6
  Gomory: 5
  Cover: 52
  Implied bound: 71
  MIR: 15
  StrongCG: 2
  Inf proof: 1
  Mod-K: 46
  RLT: 1
  Relax-and-lift: 1

Explored 3923 nodes (104094 simplex iterations) in 2.21 seconds (4.20 work units)
Thread count was 8 (of 8 available processors)

Solution count 10: 227.157 228.048 228.126 ... 255.478

Optimal solution found (tolerance 1.00e-04)
Best objective 2.271573287284e+02, best bound 2.271573287284e+02, gap 0.0000%

• 结果输出

Model status   : 2    
Objective value: 227.16
Vehicle capacity: 187.59
 [0, 13, 17, 18, 19, 11, 9, 6, 4, 2, 1, 26] 
Vehicle id | Route 
1          | [0, 13, 17, 18, 19, 11, 9, 6, 4, 2, 1, 26] 
4          |[0, 5, 3, 7, 8, 10, 15, 16, 14, 12, 26] 
6          | [0, 20, 24, 25, 23, 22, 21, 26] 

Vehicle id |  Route length   Travel distance    Loads      
1         |  12             97.55              184.07                   
4         |  11             93.16              199.28           
6          |  8              36.44              112.95         

• 结果可视化
  

4. VRPTW 和SVRPTW的结果分析和对比

为了方便直观对比VRPTW和SVRPTW的解的不同，我们再次将两个模型的解放在一起进行观察：

接下来我们改变$\pi$的取值，来观察随机模型的解的变化。

对于本文设定的算例而言，由于随机模型的需求量均值${\mu }_i$设定为与确定模型中的需求$d_i$相等，因此随机模型和确定模型的区别即在于车辆容量约束的容量限制$\overline{Q}$。由式(22)可知，在给定车辆载重未超过容量上限的概率$\alpha$（由式(16)可知，即给定$\tau$）和车辆实际容量上限$Q$情况下，$\overline{Q}$由参数$\pi$决定。下表即为数值试验中不同$\pi$的设定下，模型的目标函数值obj、车辆容量约束的容量限制$\overline{Q}$以及每辆车的实际载重量load。

$\pi$	$\overline{Q}$	obj	load
0.0	200.00	191.81	[160.0, 110.0, 190.0]
0.1	194.35	191.81	[107.43, 189.53, 165.83]
0.2	192.06	191.81	[111.33, 198.33, 171.94]
0.3	190.32	191.81	[110.57, 187.32, 156.46]
0.4	188.86	227.16	[192.2, 178.14, 113.74]
0.5	187.59	227.16	[198.05, 178.32, 103.5]

由表中结果可得:

• 若$\pi =0$，则随机模型退化为确定性模型，SVRPTW和VRPTW有相同的解。
• 若$\pi$越大，对应的$\overline{Q}$越小，模型的目标值也更大，表现更差。

直观来看，由于${\sigma }_i^2=\pi {\mu }_i$，因此$\pi$越大，则表示需求的方差越大，则波动越大。$\overline{Q}$代表了考虑随机性的情况下对于车辆的容量约束，如果需求波动大，则需要将$\overline{Q}$降低以保证前期分配的运输量比较保险，不会在实际运输时出现容量超载的情况。

5. 小结与思考

1. 随机性模型和确定性模型的不同

随机性存在的原因是，决策者在做决策时，真实需求尚未可知，仅能根据历史数据等信息作出估计，而估计的需求与真实实现的需求往往存在一定差异。本文中介绍的机会约束就是一种将输入数据不确定性纳入考虑的建模方法。

2. 随机性模型的建模方法和处理思路。

• 建模方法：
  考虑随机性的VRP问题可以用机会约束规划、带补偿的随机规划、动态规划、鲁棒优化等方式建模。其中本文采取机会约束规划的建模方式，将确定性模型中的车辆容量确定性约束，转化为机会约束，即该约束成立的概率要大于阈值。
• 处理思路：
  处理机会约束的一般思路为：通过模型转化，将带有概率表达式的约束转化为确定性的不等式或者其他能够直接处理的形式，从而使模型可以直接求解。

3. 随机性造成的影响。

随机模型相比于确定性模型直接求解，预期的成本会增加，但是模型的鲁棒性有所提升，能够抵抗一定程度的风险。在参数波动时，随机性模型的解仍然能够保持较好的表现，但是，确定性模型的解，有可能随着参数的波动产生较大的变化，有时甚至会变得不可行。随机模型和确定模型之间的差异，也可以视为是为了抵抗随机性付出的代价。

完整代码获取方式

本文实现了python调用Gurobi版本的VRPTW和SVRPTW，完整代码可以通过朋友圈集赞获得。转发该推文至朋友圈，集齐30赞，然后将集赞截图发送至公众号后台，并发送姓名+学校+电话+邮箱至公众号后台，即可免费获得完整的代码。

参考文献

[1] Desaulniers, Guy, Jacques Desrosiers, and Marius M. Solomon, eds. Column generation. Vol. 5. Springer Science & Business Media, 2006.

[2] Stewart Jr, William R., and Bruce L. Golden. “Stochastic vehicle routing: A comprehensive approach.” European Journal of Operational Research 14.4 (1983): 371-385.

[3] Gendreau, Michel, Gilbert Laporte, and René Séguin. “Stochastic vehicle routing.” European Journal of Operational Research 88.1 (1996): 3-12.

[4] A. Chames and W. Cooper, Chance-constrained programming, Management Sci. 6 (1959) 73-79.

[5] A. Charnes and W. Cooper, Deterministic equivalents for optimizing and satisficing under chance constraints, Operations Res. 11(1963) 18-39.