知识分享|非凸问题求解方法及代码示例【分类迭代】【大M法】

主要内容   

之前发布了非线性问题线性化的几种方法，如知识分享|分段函数线性化及matlab测试，学习园地 | yalmip实用操作-线性化，非线性优化 | 非线性问题matlab+yalmip求解案例，但是在实际建模及编程过程中，会遇到各种各样的问题，比如下面的模型就出现了非凸的报错问题，主要原因在于目标函数为两个变量乘积。

该问题直接写成matlab代码如下：

clear;clc;
% ====== 变量定义 ======%
x=sdpvar(1,2,'full');
r=sdpvar(1,2,'full');
b1=binvar(2,1,'full');
b2=binvar(2,1,'full');
% ====== 约束表达 ======%
P=[30,70];
C=[];
C=[C,x(1)+x(2)==100];
C=[C,0.5*P(1)<=x(1)<=3*P(1),0.5*P(2)<=x(2)<=3*P(2)];
C=[C,sum(b1)==1,
    implies(b1(1),[x(1)>=0.8*P(1),r(1)==0.0200]);
    implies(b1(2),[x(1)<=0.8*P(1),r(1)==0.0162])];
C=[C,sum(b2)==1,
    implies(b2(1),[x(2)>=0.9*P(2),r(2)==0.0195]);
    implies(b2(2),[x(2)<=0.9*P(2),r(2)==0.0130])];
C=[C,x(1)>=0,x(2)>=0];
% ===== 目标及求解 ======%
F=(-1)*(r(1)*x(1)+r(2)*x(2));
ops = sdpsettings('solver', 'cplex');  % gurobi mosek
sol=optimize(C,F,ops);

求解之后会出现报错，主要是由于模型非凸所致，为了更好的求解该问题，需要对模型采用其他方式进行求解。

本次讨论两种解决办法，一种是分类迭代求解，另外一种是采用BigM法进行线性化处理。

1 分类迭代求解

通过观察能够看出来，r的非线性分四种情况，因此，可以分别求解四种情况下的最优解，然后通过比较确定最终最佳方案，代码如下：

clear all;clc;
P=[30,70];
Fc=200;
for i=1:4
    % ====== 变量定义 ======%
    x=sdpvar(1,2,'full');
%     r=sdpvar(2);
    % ====== 约束表达 ======%
    C=[];
    if i==1
        r(1)=0.0200;
        r(2)=0.0195;
        C=[C,x(1)>=0.8*P(1),x(2)>=0.9*P(2)];
    elseif i==2
        r(1)=0.0200;
        r(2)=0.0130;
        C=[C,x(1)>=0.8*P(1),x(2)<=0.9*P(2)];
    elseif i==3
        r(1)=0.0162;
        r(2)=0.0195;
        C=[C,x(1)<=0.8*P(1),x(2)>=0.9*P(2)];
    else
        r(1)=0.0162;
        r(2)=0.0130;
        C=[C,x(1)<=0.8*P(1),x(2)<=0.9*P(2)];
    end
C=[C,x(1)+x(2)==100];
C=[C,0.5*P(1)<=x(1)<=3*P(1),0.5*P(2)<=x(2)<=3*P(2)];
C=[C,x(1)>=0,x(2)>=0];
% ===== 目标及求解 ======%
F=(-1)*(r(1)*x(1)+r(2)*x(2));
ops = sdpsettings('solver', 'cplex');  % gurobi mosek
sol=optimize(C,F,ops);
% ===== 判别最佳结果 ======%
F1=value(F);
x1=value(x);
r1=value(r);
if Fc>F1
    Fc=F1;
    xc=x1;
    rc=r1;
end
yalmip('clear');
end

最终得到的F最优值为-1.9685，x为37和63。

2.BigM法非线性处理

因为该问题的目标是两个变量乘积项导致的非线性问题，可以通过线性化方式进行转化：

令

则关于v1的约束可以转化为如下形式：

同理，关于v2的约束也可转换为上述形式。

具体代码如下：

clear;clc;
% ====== 变量定义 ======%
x=sdpvar(1,2,'full');
r=sdpvar(1,2,'full');
beta=binvar(1,2,'full');
v=sdpvar(1,2,'full');
b1=binvar(2,1,'full');
b2=binvar(2,1,'full');
% ====== 约束表达 ======%
P=[30,70];
M=100;
C=[];
C=[C,x(1)+x(2)==100];
C=[C,0.5*P(1)<=x(1)<=3*P(1),0.5*P(2)<=x(2)<=3*P(2)];
C=[C,-beta(1)*M+0.02*x(1)<=v(1)<=beta(1)*M+0.02*x(1)];
C=[C,-(1-beta(1))*M+0.0162*x(1)<=v(1)<=(1-beta(1))*M+0.0162*x(1)];
C=[C,x(1)>=0.8*P(1)*(1-beta(1))-beta(1)*M];
C=[C,x(1)<=0.8*P(1)*beta(1)+(1-beta(1))*M];
C=[C,-beta(2)*M+0.0195*x(2)<=v(2)<=beta(2)*M+0.0195*x(2)];
C=[C,-(1-beta(2))*M+0.013*x(2)<=v(2)<=(1-beta(2))*M+0.013*x(2)];
C=[C,x(2)>=0.9*P(2)*(1-beta(2))-beta(2)*M];
C=[C,x(2)<=0.9*P(2)*beta(2)+(1-beta(2))*M];
C=[C,x(1)>=0,x(2)>=0];
% ===== 目标及求解 ======%
F=-sum(v);
ops = sdpsettings('solver', 'cplex');  % gurobi mosek
sol=optimize(C,F,ops);

最终得到的F最优值为-1.9685，x为37和63，和分类迭代得到的结果是一致的。