CPLEX中MILP模型灵敏度分析实现及代码报错问题咨询

一、代码存在的问题 #

1. 原理性限制：CPLEX原生的.slack（松弛变量）、.dual（对偶值）属性仅支持纯线性规划（LP）问题，你的模型包含布尔型、整数型决策变量，属于混合整数线性规划（MILP），直接调用该属性必然触发不兼容报错。
2. 语法错误：你所有约束都定义在forall(t in 1..Time, j in J)循环内，属于带索引的约束数组，直接写c1.slack没有指定具体的t、j索引，哪怕是纯LP模型也无法正常运行。
3. 约束定义冗余：仓库平衡约束c2、库容约束c3、扩张/收缩约束c5-c7、产能约束c8-c9、初始值约束c10-c12均和产品索引j无关，你把它们写在遍历j的循环内，会导致同一约束被重复定义Prod次，属于逻辑错误。

二、MILP可采用的灵敏度分析方法 #

• 固定整数变量法（最常用）：先求解原MILP得到最优解，将所有整数、布尔变量的取值固定为最优值，原模型就退化为纯LP问题，此时即可正常调用松弛变量、对偶值做常规灵敏度分析，得到的结果对应当前最优整数解下的参数敏感度，参考价值最高。
• 参数扫描法：手动调整你需要分析的参数（比如需求、成本系数、库容上限等），多次运行模型记录最优目标值、决策变量的变化，直接观测参数波动对结果的影响，逻辑简单结果直观，适合待分析参数不多的场景。
• 鲁棒性验证法：调整CPLEX的最优性间隙阈值，测试参数在多大范围内变动时，当前得到的整数解仍然保持最优或满足精度要求，可用来判断解的稳定性。

三、修正代码参考 #

1. 约束块修正 #

把和j无关的约束从forall(t,j)循环中拆分出来：

subject to {
forall(t in 1..Time, j in J) {
  c1: i[t-1][j] + x[t][j] - d[t][j] == i[t][j];    //Inventory balance constraint
  c4:x[t][j] <= (sum(t1 in t..Time)d[t1][j])*y[t][j];  // dtT Cumulative demand from t to end T 
}

forall(t in 1..Time) {
  c2: w[t-1] + u[t] - v[t] == w[t];   //warehouse balance constraint
  c3:sum(j in J)i[t][j] <= w[t];   // inventory should be max. the warehouse level
  c5:u[t] <= (sum(t in T,j in J)d[t][j])*z1[t]; // Total demand from 1 to T * z1
  c6:v[t] <= (sum(t in T,j in J)d[t][j])*z2[t]; // Total demand from 1 to T *z2
  c7:w[t] <= 14000; // WH max available capacity 
  c8:x[t][1] <= 8000; // production capacity prod 1 in a period t 
  c9:x[t][2] <= 9000; // production capacity prod2 in a period t 
}  

c10:i[0][1]==100;   // initial inv prod 1
c11:i[0][2]==500;   // initial inv prod 2
c12:w[0]==1000;  //initial warehouse size 
}

2. 后处理及求解流程修正 #

新增主求解流程，先求解MILP再固定整数变量转LP，再循环输出带索引的松弛变量：

main {
  thisOplModel.generate();
  // 求解原MILP
  if (cplex.solve()) {
    writeln("原MILP最优成本：", cplex.getObjValue());
    // 固定所有整数/布尔变量取值，转为LP
    cplex.solveFixed();
    // 执行后处理输出松弛变量
    thisOplModel.postProcess();
  } else {
    writeln("原MILP无解");
  }
}

execute POSTPROCESS {
  for(var t=1;t<=Time;t++){
    for(var j=1;j<=Prod;j++){
      writeln("c1(t=",t,",j=",j,") 松弛值：",c1[t][j].slack);
      writeln("c4(t=",t,",j=",j,") 松弛值：",c4[t][j].slack);
    }
    writeln("c2(t=",t,") 松弛值：",c2[t].slack);
    writeln("c3(t=",t,") 松弛值：",c3[t].slack);
    writeln("c5(t=",t,") 松弛值：",c5[t].slack);
    writeln("c6(t=",t,") 松弛值：",c6[t].slack);
    writeln("c7(t=",t,") 松弛值：",c7[t].slack);
    writeln("c8(t=",t,") 松弛值：",c8[t].slack);
    writeln("c9(t=",t,") 松弛值：",c9[t].slack);
  }
  writeln("c10 松弛值：",c10.slack);
  writeln("c11 松弛值：",c11.slack);
  writeln("c12 松弛值：",c12.slack); 
}

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