欢迎参加本次网络讲座（本字幕根据声音翻译，发现问题请联系help@gurobi.cn)
讲座的内容是关于数值计算问题
Gurobi 已经有十年的发展时间
Gurobi 团队在优化方面有超过20多年
丰富的数值计算和优化经验
实际上我们不断在帮助客户日常
使用Gurobi中积累更多的经验
客户的应用范围很广
从会计到程序编制
无所不包
今天的讲座得益于这些经验
让我们先看一下Gurobi 关于数值问题的建议
或者说发生问题的原因
更重要的是，如何避免
什么是数值问题
比较难确切定义
每当用户发现与预期不符的结果
或者同一个问题出现不一致的结果
就会想到可能是数值问题
但造成问题的原因有很多
最常见原因是这四个，可以容易处理
首先是四舍五入
然后是理论上的数值和计算机能够理解并且表征的数值之间的巨大差异
这个差异反映出我们容易忘记这些数值都有局限性
而这些局限性带来了精度预期的差异
最后我们讨论不恰当的条件设置
首先我们先理解这些问题的源头
做容易发生的是系数的四舍五入
可以用一个例子说明，这个问题可以很严重
这是一个简单的例子
我假设你输入了二个约束条件
看看这二个约束条件
你也许会说他们看一起来一样
数值非常接近, 1/3 和 0.333 一样
优化器解读你提供的数值是真实的数值
通过简单的代数替换
你会发现这二个约束条件意味着
x=0, y=1 是唯一的结果
你提供的系数是优化器解读的确切数值
换一种更高的近似精度
例如用这二个约束条件
描述了不同的故事
简单的代数替换产生这样的表达式
这样的表达式非常非常非常接近1
优化器认为等于1，因此二个约束条件冗余，删除其中一个
只要满足这个条件都是可行的
如果你采用之前的精度，那告诉优化器只能选择
x=1,y=0
但往往这不是你希望的结果
因此寻找问题首先考虑模型中的数值是不是全精度的
也就是16位小数
这是我采用的标准精度
一个好问题可以问自己
为什么 2*1E-16这个数值如此特殊
以致于Gurobi 会忽略它
这就说明这些浮点数据不是真正的实数
或者不是计算机处理的方式
有一个简单的验证方法
打开 Python 环境，输入这个表达式
会得到一个 True的结果，但看起来
不应该为 True
但这就是计算机理解的
也许你认为这可以是由于数值太小
近似为1不会有太大错误
如果数值很大的话就不会出现这个问题
但并非如此
做另外一个实验将 1E16 加 1 与 1E16 对比
Python 又一次说这是相等的
但奇怪这些都不是小数值
简而言之
这些和Python 无关
在 Excel,C,R,Java 等任何编程语言中也如此
只要使用标准的浮点表达式，结果都为True
这似乎与我们在学校中学习的基本代数规则相悖
再例如代数中的结合律规则
1 加上这个很小的数值之后再加上这个很小的数值
并不等同于这二个很小的数值相加之后再加1
人们的理解和计算机硬件存储数值的方式大相径庭
造成期望的结果和实际结果不一致
我相信有些变通的方法
例如 GNU-MP 和其他一些工具
但是这些工具依赖软件实现，速度会慢
你也许知道局限性的来源
在 Python 中，你输入 import system
然后打印出浮点信息
可以看到这里有个数值 epsilon， 非常小的数值
这个数值是可以加到 1 数值上，但仍然和 1 可以区分开的最小数值
比这个数值再小就无法存在计算机里，我们无法使用
这就是我说的避开不了的地方
也就是我说的实数不实的地方，至少在计算机里是这样
因此我们需要比较谨慎
在计算机里进行数学计算
稍微解释一下，这意味着
大多数系统提供了硬件实现数值近似的方法
这些硬件是取得好结果的关键
这些硬件处理可以比其他软件方法快10，100甚至1000倍
所以如果你在意性能
应该在限度内使用这些数值
这也是 Gurobi 采用的64位双精度数值原因
效果还不错
我可以演示给你看，如果合理使用
99.999% 的情况下结果有效
还有最后一点，正因为这些数据不是确切数值
将他们加在一起不是一个精确的运算
每当进行数值比较时
利用等式比较不一定准确，需要考虑到可能出现的精度误差
一个可行的方法是利用客户数值范围上的容差量 Tolerance
再强调一下，之前说过，Gurobi 认为
所有客户提供的数值都是准确的数值
如果有用户在模型中输入了1 加 1E-7 再乘以 x 小于或等于 1
x 是一个 0-1 变量
Presolve 也许会检测到这点，因为这个是不可见的
然后决定这意味着 x 只能等于0
从传给优化器的数值来看这个决定不错
我需要强调的是这种检测只是可能会发生，并不是一定发生
原因是在 Presolve 和 切平面时
我们需要限制花费在这种检测上的时间
不是任何时候都会做全部检测
但是当找到一个可行解，需要检测这个可行解是否真正可行
之前是说过，不能在约束条件上用等式检测
需要引入容差量 Tolerance 的概念
每当有一个可行解
至少需要检测三个方面
第一个检测是否是整数
我们定义了 IntInfTol
用于检测是否满足整数条件
同样，原始单纯形可行性检测引入了容差量 FeasibilityTol
每当检测不等式或者等式是否成立时
我们都留有一定空间用于判断
同样，用于检测对偶可行时引入了 OptimalityTol
这三个参数用户是可以调整的
但这些调整比较精细，需要谨慎说明
首先这些容差量都是绝对数值，不依赖输入的系数而改变
如果你改变了放大或缩小了约束和变量的系数倍数
你需要小心，这些容差量的含义也发生变化
同时也意味着一个可能满足之前容差量的可行解
可能会被 Presolve 或者切平面认为是不可行解而切掉
为什么这样会带来麻烦
因为这意味着对于数值比较敏感和具备挑战性的模型
最优解可能会随着算法的选择而变化
这意味着对于很难获得准确解的问题
存在一个灰色地带
在灰色地带里有时会被认为是可行解，有时则不可行
也就是说那些严格可行的仍然可行，那些在容差量之外的则判为不可行
每当问题存在一个较大的灰色地带时
往往容易出现数值问题
但到底这种情形会有多槽糕
目前为止听我说得似乎很糟糕
但使用优化器的人并不经常遇到这个问题
怎么会这样？应该有什么原因
的确是
首先Gurobi 默认的原始和对偶容差量是 1e-6
整数是1e-5
假设你定义的模型右端数值在 1e4 和 1e3 级别
意味着，有时可行有时不可行的灰色地带会很小
大概只有 1e-9 或者说十亿分之一
这个数值很小，对于大部分问题无足轻重
灰色地带不起什么作用
这就是我说的scaling（系数缩放）很关键
因为你来决定约束和变量的系数和变化范围
这里的容差量是绝对数值，而不是相对数值
如果你的约束变化范围很窄
优化器可能认为是等式
例如解在 1e-10 和 1e-9 之间移动
优化器可能就无法检测到
所以约束变化范围要明显一些，除非他们的确就是等式
对于变量而言，同样适用
变量范围不能太小
这些范围会限制原始单纯形约束的边界
同样变化范围也不要太大
特别对于整数变量
例如范围是100万
最好处理为连续变量
移动一个单位，对于百万级别而言效果微不足道
但优化需要付出的代价却是很高
对于目标函数而言，这些规则也适用
当你检测可行解的时候
也在检测对偶条件
确保线性最优性
我的建议是将问题处理成
最优目标取值范围在 -1e4 和 1e4 之间
否则在对偶可行性上需要增加很多计算负担
例如假设我们在处理百亿美金级别的问题
一个好的处理方法是改变目标函数中单位量级
比如以千或者万为单位
而不是元或者分为单位
否则目标函数将失控
造成优化困难
和之前所述约束变动范围不要太小的原因一样
尽量避免有意义的目标函数取值太小
我建议要么取值大于1或者小于-1
避免太接近0
这样的话在对偶部分容易处理一些
除非你的问题是约束满足问题，所有目标函数系数都为0
这样的话也不会牺牲性能
以上这些规则并不难实现
完全取决于你在行、列和目标函数中
如何调整系数缩放和单位
之前说过，Gurobi不会自动做这些工作
用户需要自己定义合适的范围
那么如何定义紧凑合适的约束和变量范围呢
首先需要利用模型中没有包含的信息
让约束和变量更紧凑
例如需要优化一个快递问题
理论上车辆的数量可以取任何数值
但如果我了解问题背景，知道只有三个地点可以存放车辆
每个地点不能容纳超过100辆车
就可以将这个信息输入给优化器
优化器就可以很好的利用这个信息缩小搜索
让优化快起来
如果没有这个信息
优化器不得不考虑各种可能情况，而有些情况现实中根本不会出现
因此首要事情是获得额外的信息
缩小变量的变化范围
然后设置合适的量纲
例如从以吨为单位到以2000吨为单位
以元为单位到以千元为单位
合适的量纲可以让约束和变量具有合适的变动范围
希望不会超过 1E4
最后还有一种情形，在目标函数中容易设置很大或者很小的系数
就是当用户处理多个层级目标时
将这些目标混在一起
最重要的目标给予非常大的系数
次重要的目标给予小的系数
会造成目标函数中差异很大的二个部分
现在 Gurobi 直接支持多目标优化
每个目标函数可以设置合理的数值
合理的系数
而且容易处理
我再强调一下为什么缩放很重要
可以快速做一个试验
这是一个模型，我将系数都缩放一个随机比例
理论上优化值不会变化
这里是修改目标函数的系数
这是缩放变量的上下界，如果是有限的话
运行这个代码，这部分是一个普通模型程序
这里是说明缩放比例是小于1的随机数
启动优化，会马上看到一个警告
提示最小系数和最大系数的差异过大
但Gurobi 可以很好求解
找到最优解，经过了1000次左右迭代
最优值是4.49等等这个数值
现在我放大系数到1000
矩阵系数差放大到 1e18 次方
仍然会看到警告
我们仍然可以得到同样最优解
但是，关键是，需要迭代的次数增加了
如果我放大更多，比如说100万
矩阵系数波动范围更疯狂
Gurobi 仍然可以获得最优值，但时间更长
同时，会出现更多警告信息
提示说尽管很努力处理非缩放时的约束违反
但在优化后无法找到更小的原始和对偶违反
对于这个例子这属于极限情况
最优目标值仍然正确，只是迭代次数又增加
无法找到一个让约束违反更小的结果
如果我们再走一步，缩放范围到1亿
这时候问题彻底崩溃
Gurobi 无法求解，只能显示不可行
Gurobi 信任模型中用户提供的数据
理论上缩放不会影响最优结果
但是数值精度的存储方式决定了缩放数值改变了这些数值的含义
面临着无法获得可行解的风险
同时请记住放大范围会让问题变得更难求解
第一次虽然Gurobi 可以求解，但花费时间更长了
所以再强调一下，合适地缩放模型
可以带来更快的优化速度
用户可以做的事情包括
对行和列做合适的缩放，设定合适的系数
我举例说明
请看这个问题
我可以用 1e4x' 替换掉 x
替换后的第一个约束条件如下
同样，可以简化第二个约束条件
这样可以从系数矩阵指数差10个级别降到4个级别
同时精度和容许值是绝对数值
所以即便问题一样，求解缩放后的问题更容易一些
因为缩放前的指数差更高
但我要求的容差量只有1e-6
而对于这个约束而言
1e-6 容差量就容易让优化器找到 x'和 z
不要骗自己。一个普适的原则就是
缩小行和列系数的变化范围
下次遇到这样的问题，很容易和自己说
我可以耍些小技巧，消除掉那些巨大数值
与其有 1e6 这样的数值，而我只有1和10
这样公式看起来更不错，不是吗
但是，你没有改变任何东西
同样不可行解的 epsilon 仍然存在于模型中
Presolve 会检测到这点，改写为原来的样子
所以这个小手段并没有改变问题
你可以采用的方式是
用 Y' 替换 Y
这是他们之间的表达式
系数矩阵的指数差降到3，但变量范围变大
这样好处是如果 Y'有epsilon 误差
则 X 的可能误差会很小，只有 1e-3，这很好
这会改善数值状况
而前面的简单替代对问题没有丝毫帮助
关于大M 约束
原则上如果你合适缩放约束和变量
就不应该有大M，因为变量限制在1e4范围内
计算处理一切良好
但出于某种原因避免不了
例如这里大的系数
而你无法重新缩放 x 的含义
你可以采取提升数值准确度的方法是
定义这些约束为 SOS 约束或者广义指标约束
这些 Gurobi 7.5 都支持
但是，虽然你会得到更可靠的结果
事实上这里意味着如果Y是0，则X也为0
这个条件也会满足
代价是需要花费额外的计算时间
所以保证数值稳定性的代价，是额外的计算时间
让我们做一个快速总结
变量取值一般来说不要超出 -1e4 到 1e4 范围
同样适用于约束条件
意味着约束条件在这个范围内移动
也适用于目标值，不应该超出这个范围太多
你也不应该比较二个绝对值差异极大数值
例如比较1e-10 和 1
会带来计算上的麻烦
尽量避免这些
如果你按此修改模型
你也许会惊喜发现性能有了提升
这是因为Gurobi每当发现数值问题
会努力解决，或者用更高的精度
或者内部调整算法
例如关闭一些原来让求解更快的代码而解决数值问题
所以选择合适的数值范围是改善性能的关键
一定要试试
接下来，如果你已经有了模型，怎么检查是否有数值问题
这是一个重要的问题
第一件事情是将模型和参数分开
设置这个参数 GURO_PAR_DUMP 等于1
在调用optimize()函数前设置
优化器会产生二个文件 REW 和 PRM
REW 是隐去所有变量和约束名称的模型文件
PRM是非默认取值的参数文件
这些参数很重要
用于指导优化器管理和选择何种算法，如何去优化
这二个文件可以让优化结果重现
获得这二个文件之后
可以在 Python 交互界面读入模型
输入 printStats 命令打印模型统计数据
会显示矩阵和目标系数范围等信息
你可以判断这些范围是否合理，是否有异常情况
如果看起来合理
可以读入PRM 参数文件
然后启动优化
这时候观察输出日志
Gurobi 输出日志提供了很多有用的信息
很多情况下提示了很多警告信息
常见的警告信息有这些
矩阵系数或者目标系数看起来不太对
或者你看到 Markowitz Tolerance tightened 到某个数值
这是一个明确信号，Gurobi 在努力寻找稳定解
这也提示模型的确有数值问题
如果看到 switch to quad precision
意味着Gurobi不但通过调整算法无法求解
还不得不扩展浮点精度
我稍后还会讲这个
当然如果你看到 Numeric error
或者 numerical trouble encountered
或者 restart crossover 这意味着
内点算法(障碍算法）遇到了困难
还有其他这些警告信息（请看PPT)
意味Gurobi无法满足容差设置
这些警告信息都是强烈信号，你应该看看模型哪里有问题
同时，你可以输入 printquality() 检测解的质量
检查边界违反的最糟糕情况
最糟糕的约束违反情况
最糟糕的整数条件违反情况
对于这个例子
边界情况不错，只有 1e-8
但是约束违反达到了 1e-4 接近 1e-3
这很糟糕
整数条件还不错
这说明得到的解不是特别的贴合
Gurobi 在努力找到一个解
用于满足容差数值
你也可以查看条件值 KappaExact
我稍后会讲
条件值如果很大，也是一个优化器遇到数值问题的信号
其他一些症状包括
修改参数，导致解发生变化，优化值发生变化
问题从不可行变成可行
这些都是强烈的信号说明你的模型有数值问题
绝大部分问题的性质都是数值方面的
意味着一般不是 Gurobi的 bugs
而是模型中的数值问题
可以通过一个小实验看到
就是缩小容差值的范围
例如从 1e-6 缩小为 1e-8
优化器可以稳定求解
这是一个强烈信号说明你的模型在缩放系数方面有欠缺
虽然减少容差值可以增强求解的稳定性
但代价是需要花费额外的数值计算时间
这点不要忘了
如果你已经改善了各种缩放比例
模型看起来不错
仍然出现很奇怪的问题
可以联系我们
以上是关于如何检测和评估模型数值问题
如果发现数值问题，你可以做些什么
Gurobi 也有一些参数可以帮助这些问题
但不是解决这些问题，而是变通化解
付出的代价就是额外的优化时间
首先说一下 Presolve(预优化)
Presolve可能会带来负面影响
检测的方法是导入Presolve 之后的模型
然后打印出模型统计参数
如果统计参数现实系数看起来很奇怪
那么Presolve作用是负面的
这种情况不经常发生，但有时会出现
如果出现这种情况
可以关闭掉 Presolve
或者将 Aggregate 设置为0
Aggregration(聚合)就是我们之前提到的那个自欺欺人的例子
通过多个约束条件，隐藏过大的系数
聚合的效果就是将这些多个约束整合在一起
改写为整合后的系数
你可以关闭这个参数
但是你需要记住，关闭参数没有消除模型本质上的过大系数
仍需谨慎
或者将 AggFill 关闭
Gurobi 有时会改写约束
也许会引入一些不合适的系数
如果是这样的话，关闭这个参数
看看是否有帮助
同时，考虑采用合适的算法
通常来说内点法（障碍法）速度最快，但不绝对
但它对数值问题比较敏感
所以你可以尝试原始、对偶和内点各种方法
看看是否有帮助，可以得到稳定一些的结果
可以采用这些参数选择根节点的算法
同时你也可以让算法对数值问题不太敏感
最简单的方法是采用参数 NumericFocus
默认值是-1，意味着Gurobi 自动内部处理
你可以设置为0，全部关闭
也可以设置为1,2,3
Gurobi 根据这些数值
决定数值计算的保守程度，包括
MarkowitzTol 处理单纯形算法
也决定什么时候转到 Quad 精度
它也会限制 GomoryPasses
减少可能会在一些具有挑战性的数值计算问题中带来问题
它也会限制CutPasses 和 CutAggregation
等等
所以这些是遇到数值计算问题时常用的方法
可以将 NumericFocus 设置为1到3
如果你希望改变单纯形法的定价方式，可以修改 NormAdjust
如果内点法（障碍法）遇到麻烦
可以采用 BarHomogeneous
也可以尝试 CrossoverBasis
或者关闭掉 Cuts
在一些具有数值挑战性的问题中，这些参数可能会导致数值问题
再重复一次，这些参数不是为了修复数值问题
他们只是面对数值问题的变通处理方法
真正修复数值问题的方法是回到模型中，建立一个合适的模型
这是本次讲座的最后一部分
关于条件数
到目前为止问题要么来源于四舍五入
没有提供足够多的精度
要么计算机数值表达的局限性
无法让一切严丝合缝
但是还有些问题
源自问题的几何涵义
这就是和条件数有关
也就是之前看到的 Kappa
优化问题的条件数
用来捕捉优化问题或者优化解
对于输入数据微小变化的敏感度
通常的定义，特别对于线性系统
称之为 Kappa，定义为
矩阵的范数和逆的乘积
我更喜欢几何意义上的解释
这和可行域圆形规整有关系
事实上你可以想象为
包含可行域的最小球
和被可行域包含的最大球
他们之间的比例关系
举二个例子
比如你可以通过很多很多边界线性不等式
建立这样的模型
条件数的几何意义就是外围淡灰色的圆
和内圈深灰色的圆之间的比值
这很接近1
这个条件数说明
如果比值很小，那么输入数据中的微小变化
导致解的变化也会很小
事实也如此，这是我的代码
我用1000个不等式定义了圆形区域
在系数上星星点点做了一点修改
然后再优化
看看原始单纯形解上有什么变化
这是500次迭代后得到的
系数的微小变化（1e-6) 导致原始解同等量级的变化
从几何意义上也可以看出
这也很合理
这些问题很容易解决
但如果我需要优化的可行域
是一个瘦长的三角形
可以想象如果我拉伸这个三角形
里面的圆形会越来越小
外面的圆形会越来越大
比值会越来越糟糕
我定义这个问题
优化它
然后在目标函数上做非常微小的变化
虽然只有 1e-6 的变化
但在原始单纯形解上有16次方的改变
这并非错误结果
其实后面这三列可以显示原始、对偶和其他边界违反的数值
几乎接近0
说明结果是正确的，但问题本身很糟糕
这说明哪怕在目标函数中改变很小的数值
都会带来原始解的巨大变化
这个和数值问题无关
这是和几何形态有关
这又回到之前提到的模型缩放问题
如果你的模型缩放合适
那么会看起来四四方方，很规整
那么一切就容易优化
如果你有一个很狭长的区域
那么你可能面临结果不稳定的情况
我们做一个总结
事实上数值稳定性源头不仅仅来自于数值
也可能来自于几何形态
Gurobi 不会对用户的数据进行猜测
用户提供什么精度，Gurobi 就采用什么精度
大部分情况下，处理数值问题最好的方式
是审视模型做恰当的缩放
或者改写模型公式
感谢你关注这个讲座（本字幕根据声音翻译，发现问题请联系help@gurobi.cn)