提高gurobi建模的数值稳定性

详细参考： https://www.gurobi.com/documentation/current/refman/guidelines_for_numerical_i.html

1 各种条件条件影响求解的难度和准确性：

（1）计算机舍入误差的存在
（2）求解器验证解的完整性（1.0e-5）， 可行性(1.0e-6)， 最优性(1.0e-6)的所使用的容忍误差的精度

2 构建一个健壮的模型是无比重要的。 首先gurobi的官方文档建议了如下方案：

2.1 约束矩阵右端项和目标的数量级

约束矩阵的右端项的大小应该 小于1.0e4的范围内
最优解的数量级在1.0e4的范围内

2.2 谨慎选择大M的值

大M方法在线性规划中广泛使用， 特别是当x >0时， 二进制变量y=1, 否则y=0. 看下面的反例：

当x=9.999 时， y可以等于0.000009999， 结合完整性的误差1.0e-5, 求解器会判定为0。 当然可以通过调整IntFeasTol 来解决这个问题。
比较好的做法是M取值为x的上界， 如果这两个方法都行不通， 则使用 SOS constraints。

2.3 约束矩阵的取值范围

约束的系数矩阵的值范围建议在**[1.0e-3, 1.0e6]** 之间，可通过变量的取值范围的放缩来调整其系数的取值范围。 举例：

这个例子中， 系数的取值范围为[1.0e-7, 1.0e4], 可以对x进行放缩例如 $x=1.0e^5{x}^{\prime }$， 那么约束将变为下面的形式：

此时，系数的取值范围为[1.0e-2, 1.0e2]

2.4 存在数值问题的debug策略

gurobi 求解问题时有个presolve 阶段， 在极端的情况下Aggregate 的开启会导致数值误差的积累。如果有数值问题， 可以通过设置参数 Aggregate 和Presolve的值来对比是否有改善。
首先， 读取模型并打印presolved模型的统计信息。

m = read('gurobi.rew')
p = m.presolve()
p.printStats()

如果数值范围看起来比原始的模型更糟糕， 那么尝试设置 Aggregate=0后， 再次检查统计信息。

m.reset()
m.Params.Aggregate = 0
p = m.presolve()
p.printStats()

如果结果仍旧不好， 那么通过设置Presolve=0干脆关掉presolve.

m.reset()
m.Params.Presolve = 0
p = m.presolve()
p.printStats()

如果设置Aggregate=0 or Presolve=0后，通过统计信息来判断有所改善，我们需要进一步测试这两个参数。对于LP问题， 可以直接测试， 对于MIP问题， 需要对比LP 松弛后的文件。 参考下面的脚本。

m = read('gurobi.rew')
r = m.relax()
r.write('gurobi.relax-nopre.rew')
p = m.presolve()
r = p.relax()
r.write('gurobi.relax-pre.rew')
m.reset()
m.Params.Aggregate = 0
p = m.presolve()
r = p.relax()
r.write('gurobi.relax-agg0.rew')

需要注意的是 Aggregate=0 会帮助改善数值问题， 但是会有性能损失。 尝试使用 AggFill=0 来替代。

2.5 算法的选择

Gurobi 内置了单纯性方法和内点法来求解LP问题或者松弛后的MIP问题。

对于大型， 困难的模型， 内点法通常快一点， 但是由于数值问题的存在， 随后的crossover 算法可能会停滞（stall).

单纯形方法是另一个选择， 对数值问题较不敏感， 可以通过设置Method为1或者0分别使用对偶单纯形或者原始单纯形。

在很多优化问题中， 并不存在数值问题。 gurobi 默认支持并发优化器， 可以分别同时使用多种方法来测试对于特定的问题哪种方法更好。

2.6 降低算法敏感性的参数

最好的办法是重构模型来解决数值问题， 如果实在没有办法来解决， 可以尝试这些参数， 并且请参考 gurobi文档。 参考链接
ScaleFlag, ObjScale (All models)
NumericFocus (All models)
NormAdjust (Simplex)
BarHomogeneous (Barrier)
CrossoverBasis (Barrier)
GomoryPasses (MIP)
Cuts (MIP)

继续研究~~~