Matlab调用Cplex的二三事

最近一个师弟来探讨为何问题feasible但是有一类等式约束明显不满足的问题，细说起来大概是前一天还是infeasible，不知动了什么突然有解了，不经大喜，冷静下来定睛一看发现解有问题这样一个过程，不经回想起曾经的我，把求解优化问题归为玄学。打了几年的交道，虽然还没有深入敌人内部，但也是积累了一些应对小技巧，借此机会总结分享一下~

在具体案例分析之前还是想客观的说一句，对于包括本人在内的一些不太了解求解器算法的外行来说，常常蜜汁自信，将无法理解的bug归为玄学，事实是永远不要轻易怀疑求解器，毕竟它已经是一个相对成熟的求解器了==绝大多数时候都是因为自己不经意间的一些错误或不规范操作导致的。

下面开始干货：

情形一：问题infeasible（exitflag=-2）

这个应该是第一次满怀信心把问题丢给求解器后碰到的最常见的情形。绝大多数情况都是由于模型本身存在冲突约束，或者自己敲的时候手误导致。

1、按照模型约束的物理意义进行初步推导，检查参数设置是否合理，是否可能导致可行域为空。

2、模型本身没有问题，手误导致可太正常了，这个时候先不要慌，首先要在保持模型变量与约束结构完整的前提下，把问题降维到最小。然后可以借助两种方法：

1）cplex.writeModel('model.lp')导出模型lp文件，大概长下面这个样子，截了一部分，小问题可以通过此方法逐一检查丢给求解器的模型是不是自己想的那样。

\ENCODING=ISO-8859-1
\Problem name: CPLEX

Minimize
 obj: x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9 + x10 + x11 + x12 + x13 + x14 + x15
      + x16 + x17 + x18 + x19 + x20 + x21 - x22 - x23 - x24 - x25 - x26 - x27
      - x28 - x29 - x30 - x31 - x32 - x33 - x34 - x35 - x36 - x37 - x38 - x39
      - x40 - x41
Subject To
 c1:   2 x42 - x108 <= 0
 c2:   4 x43 - x109 <= 0
 c3:   - x42 + x110 <= 0
 c4:   - 3 x43 + x111 <= 0 。。。
Bounds
 1 <= x1 <= 4
 0 <= x2 <= 1
 0 <= x3 <= 1
 0 <= x4 <= 1 。。。
Binaries
 x2  x3  x4  x5  x6  x7  x8  x9  x10  x11  x12  x13  x14  x15  x16  x17  x18 
 x19  x20  x21  x22  x23  x24  x25  x26  x27  x28  x29  x30  x31  x32  x33 
 x34  x35  x36  x37  x38  x39  x40  x41  x42  x43  x74  x75  x76  x77  x78 
 x79  x80  x81  x82  x83 
End

2）cplex.writeConflict('conflict.clp')运用冲突检测功能，类似的可以导出冲突部分的模型，它会告诉你哪条约束不满足，以及所涉及变量的其他约束、上下限定义等。

if cplex.Solution.status == 103
     cplex.refineConflict();
     cplex.writeConflict('conflict.clp');
 end

\ENCODING=ISO-8859-1
\Problem name: CPLEX_conflict

Minimize
 obj:
Subject To
 c9:   - x44 + x74 <= 0
 c14:  - x49 + x79 <= 0
 c15:  - x50 + x80 <= 0。。。
\Sum of equality rows in the conflict:
\ sum_eq: - x42 + x44 + x49 + x50 + x51 + x52 + x53 - x59 - x60 - x61 - x62
\         - x63 + 2 x64 + 2 x69 + 2 x70 + 2 x71 + 2 x72 + 2 x73 - 6 x74 - 5 x75
\         - 5 x76 - 5 x77 - 5 x78 - 5 x79 - 4 x80 - 3 x81 - 2 x82 - x83  = 0
Bounds
 0 <= x42 <= 1   
 0 <= x74 <= 5。。。
Binaries
 x42  x74  x75  x76  x77  x78  x79  x80  x81  x82  x83 
End

可以看出要想快速找出infeasible的原因， 原问题规模太大是很不友好的， 建议对于有明显结构特征的问题， 养成从小算例开始调试再扩展的好习惯~

情形2： 问题feasible， 但exitflag=5

这个是个进阶情形了，在含二次约束的较复杂问题中比较容易碰到。 好容易跑出了解， 先别高兴得太早看看标志位是不是1（ 毕竟即使是1解都有可能有问题==），5的意思是重要方向导数小于规定的容许范围并且约束违背小于options TolCom， 说人话我理解就是可能有一些约束不满足。如果所得的解在重要约束上是满足精度要求的那么roughly可以接受，否则这个解还是没有办法用。

1）诶那我直观地想我把tolerance调小一点行不行，但发现我碰到的问题这个不太好使。

opt = cplexoptimset('cplex');
opt.feasopt.tolerance=1e-9;（默认1e-6）

2）既然我改变不了求解器的内部算法逻辑，那只能通过改变我丢给它的这个玩意儿看看是否能躲过这个坑（即用魔法打败魔法）。对于等式约束，可以尝试在可接受范围内进行松弛，然后以罚函数形式加入目标函数中（与原目标函数值量级最好相差两个量级左右），罚函数同样也适用于不等式约束的松弛（即哪里不满足强调哪里）。

情形3：问题feasible，exitflag=1，欸但是有约束不满足！

这种时候是大家最容易相信玄学的时候……表面上看似乎很不科学，作为外行的我至今不太能解释这种情况出现的原因，因此应对措施也只能case by case见招拆招，比如我之前有个问题中把自变量的范围由[-1,1]调整为[0,1]，等式约束就满足了。似乎求解器在识别正负、等与不等存在着本质区别。另外还有可能是数据截断问题导致，说到这里就不得不强调下求解过程中需要注意模型的数值问题，在一个模型中尽量不要使得不同变量的数量级差太多（标幺大法保平安，不要轻易搞出很大或者很小的数，可行域会比较畸形），运用大M法进行松弛的时候，M的取值也需注意，一般高出2~3个数量级即可，无需任性的设置为1e100这样……

此外做一些别的小总结，常见报错：

1）CPLEX Error  3019: Failure to solve MIP subproblem.求解混合整数规划问题中遇到数值问题，注意下变量的量级

2）exitflag=-8这种诡异结果可能是内存问题……

3）CPLEX Error 5002: Q is not positive semi-definite.求解二次规划问题中二次项系数矩阵Q非半正定问题，首先检查下是否问题本身就不是凸的，不是的话出门左转不送（Cplex、Gurobi只能求解凸规划），Xpress和matlab自带的fmincon可以求解一些更一般的问题。另外将二阶锥约束sqrt(xTAx)<=y改为二次约束xTAx-y^2<=0，要对y加上约束y>=0，这样Q才不会报非半正定的错误（二阶锥是一种特殊的凸约束但Q的确不是半正定）

Cplex的调用方式：

1）通过中介公司yalmip，这个我一直没用过……据说比较人性化还是比较推荐。

2）暴力赋值系数矩阵，以MIQCP问题为例，有三种方法：

（1）

cplex = Cplex();
cplex.Model.sense = 'minimize';
cplex.Model.obj   = f;
cplex.Model.lb    = lb;
cplex.Model.ub    = ub;
cplex.Model.l    = l;
cplex.Model.Q    = Q;
cplex.Model.r    = r;
cplex.Model.ctype = ctype';
cplex.Model.A =  [Aineq;Aeq];
cplex.Model.lhs = [-inf*ones(rineq,1); 0*ones(req,1)];
cplex.Model.rhs = [bineq;beq];
cplex.solve();
if cplex.Solution.status == 103
    cplex.refineConflict();
    cplex.writeConflict('conflict.clp');
end
x=cplex.Solution.x;

注意这种写法如果等式约束常数项beq不为0，cplex会新增若干变量表示beq，但这些变量取值不是等于beq而是松弛到0<=x<=beq，导致解不正确，所以只能把这部分等式约束改为不等式约束。

（2）

problem.H=[];
problem.f=f;
problem.Aineq=Aineq;
problem.bineq=bineq;
problem.Aeq=Aeq;
problem.beq=beq;
problem.qc.a=l;
problem.qc.rhs=r;
problem.qc.Q=Q;
problem.sos=[];
problem.lb=lb;
problem.ub=ub;
problem.ctype=ctype';
problem.options = cplexoptimset('diagnostics','on','mip.tolerances.mipgap',0, 'mip.tolerances.absmipgap',0);%参数设置
[x,fval,exitflag,~]= cplexmiqcp(problem);

（3）

opt = cplexoptimset('cplex');%参数设置
opt.feasopt.tolerance=1e-9;

opt.ExportModel = 'model.lp';

[x,fval,exitflag]=cplexmiqcp([],f,Aineq,bineq,Aeq,beq,l,Q,r,[],[],[],lb,ub,ctype',[],opt);

贴一些Cplex的一些常用参数设置， 具体参考https://www.ibm.com/docs/zh/icos/12.10.0?topic=cplex-parameters，学会看官方文档， 其实很多问题都可以从中找到解答。

比如设置参数a=b

写法（1）

cplex.Param.a=b

写法（2）

problem.options = cplexoptimset('a','b'/b);

写法（3）

options = cplexoptimset('cplex');

options.a=b;

常用参数：

timelimit求解时间

mip.tolerances.mipgap整数规划问题的求解gap（默认0.0001，一般可以放大一些）

今天先写到这里，内容有错误或不严谨的地方还请各路大神多多批评指正， 本人也在不断学习优化中~