CPLEX Studio中VRP鲁棒优化不确定性集实现问题咨询

现有代码核心错误 #

• 鲁棒建模逻辑完全错误：你写的execute块是脚本执行块，不能用于定义决策变量约束，当前的u变量赋值逻辑不会进入优化模型，本质是固定了偏差边，没有实现「最多Γ条边取最大延误」的最坏场景覆盖，只能算单场景解。
• 基础约束缺失/硬编码错误：
  • NB5约束硬编码排除节点6，通用场景下要删除该判断
  • 注释掉的出发时间基准约束s[0,k]=0、同节点自环禁止约束x[i,i,k]=0是模型必须约束，否则会出现无效解
  • 需求鲁棒约束NB12逻辑错误：多余嵌套了j的循环，且错误乘了旅行时间参数t[i,j]，不符合Bertsimas预算鲁棒的对偶形式
  • 旅行时间鲁棒约束NB13逻辑错误：当前写法没有通过对偶变量刻画最坏场景下的最大延误总和，只是固定加了偏差值，没有受不确定性预算G的约束。

鲁棒建模核心逻辑（对应你的Γ延误场景） #

你用的是Bertsimas提出的预算不确定性集：对任意车辆行驶路径，最多有Γ（即你代码里的参数G）条边会触发最大30分钟延误，其余边取标称旅行时间，鲁棒模型不需要枚举所有场景，通过强对偶定理可以转化为线性整数规划，核心是引入对应不确定性预算的对偶变量：

• 为每辆车k引入非负对偶变量pi_t[k]：对应旅行时间不确定性预算G的对偶变量
• 为每条边(i,j)、每辆车k引入非负对偶变量nu_t[i,j,k]：对应单条边是否触发最大延误的对偶变量
• 需求侧的鲁棒约束同理，修正你现有NB11、NB12的错误写法即可。

修正后可运行OPL代码 #

// parameters & sets
int n = ...; // 客户点数量
{int} N = asSet(1..n);  // 客户点集合
{int} N0 = {0} union N union {n+1}; // 含起点0、终点n+1的全节点集合
int m = ...; // 车辆数量
{int} M = asSet(1..m); // 车辆集合
float Q = ...; // 车辆额定载重
int r[N]= ...; // 客户点标称需求
int b[N]= ...;// 客户点时间窗截止时间
tuple Edge {int i ; int j;};
{Edge} A  = {<i,j> | i in N0, j in N0, i!=j}; // 排除自环的边集合
int K = 99999;   // 足够大的M常数，注意不要太小导致约束失效
float delay_edge[N0,N0] = ...; // 单条边最大延误值，你的场景下固定为30即可
float demand_dev[N] = ...; // 单个客户需求最大偏差
float G = ...;  // 旅行时间不确定性预算，即你说的Γ
float L = ...;  // 需求不确定性预算
float t[N0, N0] = ...; // 标称旅行时间矩阵
int c[N0, N0] = ...; // 边行驶成本矩阵

// 决策变量
dvar boolean x[N0, N0, M]; // x[i,j,k]=1表示车辆k走边(i,j)
dvar float+ s[N0, M]; // s[i,k]表示车辆k到达节点i的时间
dvar float+ p_demand[N,M]; // 需求侧鲁棒对偶变量
dvar float+ z_demand[M]; // 需求侧预算对偶变量
dvar float+ pi_t[M]; // 旅行时间侧预算对偶变量
dvar float+ nu_t[N0,N0,M]; // 旅行时间侧边对偶变量

// 目标函数：最小化总行驶成本
minimize sum(k in M, <i,j> in A) c[i,j]*x[i,j,k];

// 约束
subject to {
  // 1. 基础VRP流约束
  // 每个客户点恰好被访问一次
  forall (i in N) {
    sum (k in M, j in N0: <i,j> in A) x[i,j,k] == 1;
  };
  // 每辆车从起点0出发
  forall (k in M) {
    sum (j in N0: <0,j> in A) x[0,j,k] == 1;
  };
  // 流平衡：进入节点的流量等于流出节点的流量
  forall (i in N, k in M) {
    sum (j in N0: <i,j> in A) x[i,j,k] - sum (j in N0: <j,i> in A) x[j,i,k] == 0;
  };   
  // 每辆车最终回到终点n+1
  forall (k in M) {
     sum (i in N0: <i,n+1> in A) x[i,n+1,k] == 1;
  };
  // 车辆从起点出发时间为0
  forall (k in M) {
    s[0,k] == 0;
  };

  // 2. 需求侧鲁棒约束（预算L）
  forall (k in M) {
    sum (i in N) r[i] * sum (j in N0: <i,j> in A) x[i,j,k] + L*z_demand[k] + sum (i in N) p_demand[i,k] <= Q ;
  };    
  forall (k in M, i in N) {
    z_demand[k] + p_demand[i,k] >= demand_dev[i] * sum (j in N0: <i,j> in A) x[i,j,k];
  };    

  // 3. 旅行时间侧鲁棒约束（预算G，即最多G条边触发最大延误）
  forall (<i,j> in A, k in M) {
    s[j,k] >= s[i,k] + t[i,j] * x[i,j,k] + delay_edge[i,j] * x[i,j,k] - K*(1-x[i,j,k]) - pi_t[k] - nu_t[i,j,k];
  };
  forall (<i,j> in A, k in M) {
    pi_t[k] + nu_t[i,j,k] >= delay_edge[i,j] * x[i,j,k];
  };
  // 到达客户点时间不超过时间窗截止时间
  forall (i in N, k in M) {
    0 <= s[i,k] <= b[i];  
  };
};

其他实现思路参考 #

• Python实现：直接用docplex库（IBM官方CPLEX Python API），建模逻辑和上述OPL完全一致，不需要做逻辑改动，适合需要和其他数据处理脚本联动的场景；如果用Gurobi求解器，鲁棒对偶转化逻辑相同。
• 不建议用Excel实现：VRP属于NP难整数规划问题，Excel的规划求解插件只能处理10个节点以内的极小算例，完全无法支撑鲁棒VRP的求解需求。

参数填写说明 #

你的示例场景下，delay_edge[i,j]统一设为30即可，G设为你定义的单场景最多延误边数Γ，模型求解后得到的路径方案即可满足：任意路径上即使出现最多Γ条边30分钟延误，所有客户点的时间窗依然不会被打破，同时车辆载重不会超限。

内容的提问来源于stack exchange，提问作者user19390232