YALMIP的简单说明

时间:2022-10-22 14:37:27

最近在做论文时,涉及到最优化问题,而最优化里面很多时候涉及的是二次约束二次规划QCQP这样的非凸问题,一般地,这样的非凸问题是得不到全局精确的最优解的,需要另辟蹊径。常用的有半定松弛SDR。将非线性松弛为线性,以致可以SDP,利用CVX这样的专有凸规划工具箱解决;但是缺点也是很明显的:就是灵活性和可靠性不高,当约束条件里面还有正定性未知的矩阵时,SDR力不从心,这时引入了QCQP下的二阶锥规格SOCP,利用YALMIP可以很好解决这样最优化问题的近似解。所以有必要说一下YALMIP的简单应用:

YALMIP总共有四个流程:

1:设置变量 ;

2::设定目标函数f;

3:设置限定条件;

4:利用工具箱求。

1:设置变量

变量设置常用的有三种:sdqvar() 设置实型,intvar()设置整形,binvar()设置二进制形。以sdqvar()为例。Yalmip中最重要的指令就是sdpvar。这个指令是用来定义决策变量。用n行和m列定义一个矩阵(或标量)P,我们写做:

P=sdqvar(n,m),这是一个默认对称的实数方阵,类似的还有简单形式sdqvar(n)以及详细形式sdqvar(n,m,'symmetric');

第三个参数可以用来获得许多预定义类型的变量,如Toeplitz, Hankel, diagonal(对角线), symmetric(对称)和skew-symmetric(反对称)矩阵,有关信息请参阅sdpvar帮助文档;

当然为了获得一个完全参数化(不一定是对称)的方阵,第三个参数是必要的:P = sdpvar(n,n,'full');

一个复数形式的完全参数化(不一定是对称)的方阵:P = sdpvar(n,n,'full','complex');

同时,也可以补上第三个和第四个参数:P = sdpvar(n,n,'sy','co')

2:sdpsettings

sdpsettings是YALMIP和solver之间通信的桥梁,一般用作最优化函数optimize, solvesos, solvemoment 和solvemp的第三个入口参数(本文以solvesdp):

options = sdpsettings('field',value,'field',value,...)
optimize(Constraints, Objective, options)

比如:ops = sdpsettings('solver','sdpa','sdpa.maxIteration',100);

可以从ops = sdpsettings;指令中看到其所有选择的默认参数:

ops = sdpsettings;

ops
               solver: ''
              verbose: 1
              warning: 1
         cachesolvers: 0
                debug: 1
        beeponproblem: [-5 -4 -3 -2 -1]
         showprogress: 0
            saveduals: 1
     removeequalities: 0
     savesolveroutput: 0
      savesolverinput: 0
    convertconvexquad: 1
               radius: Inf
                relax: 0
                usex0: 0             
            savedebug: 0
                  sos: [1x1 struct]
               moment: [1x1 struct]
                  bnb: [1x1 struct]
                bpmpd: [1x1 struct]
               bmibnb: [1x1 struct]
               cutsdp: [1x1 struct]
               global: [1x1 struct]
                  cdd: [1x1 struct]
                  clp: [1x1 struct]
                cplex: [1x1 struct]
                 csdp: [1x1 struct]
                 dsdp: [1x1 struct]
                 glpk: [1x1 struct]
                 kypd: [1x1 struct]
               lmilab: [1x1 struct]
              lmirank: [1x1 struct]
              lpsolve: [1x1 struct]
               maxdet: [1x1 struct]
                  nag: [1x1 struct]
               penbmi: [1x1 struct]
               pennlp: [1x1 struct]
               pensdp: [1x1 struct]
                 sdpa: [1x1 struct]
                sdplr: [1x1 struct]
                sdpt3: [1x1 struct]
               sedumi: [1x1 struct]
                qsopt: [1x1 struct]
               xpress: [1x1 struct]
             quadprog: [1x1 struct]
              linprog: [1x1 struct]
             bintprog: [1x1 struct]
              fmincon: [1x1 struct]
           fminsearch: [1x1 struct]

在ops里面有semudi选项,然后ops.semudi,查看semudi工具箱参数例如semudi.sps=1e-12:

ops.sedumi

ans = 

          alg: 2
         beta: 0.5000
        theta: 0.2500
         free: 1
          sdp: 0
      stepdif: 0
            w: [1 1]
           mu: 1
          eps: 1.0000e-09
       bigeps: 1.0000e-03
      maxiter: 150
        vplot: 0
       stopat: -1
         denq: 0.7500
         denf: 10
       numtol: 5.0000e-07
    bignumtol: 0.9000
      numlvlv: 0
         chol: [1x1 struct]
           cg: [1x1 struct]
    maxradius: Inf

3:optimize(最优化函数本文以solvesdp为例)

optimize是解决最优化问题的常见函数,本文里solvesdp就是其中一种。格式:diagnostics = optimize(Constraints,Objective,options)

举个简单例子:我们面对一个线性规划LP:{min cTx subject to Ax<= b},用本文方法的话:

x = sdpvar(length(c),1);
F = [A*x<=b];
h = c'*x;
optimize(F,h);
solution = value(x);
如果我们仅仅考虑灵活性那么可以忽略目标函数h,得到: diagnostics =  optimize ( F ) ;返回的diag可以用在检验灵活性:

if diagnostics.problem == 0
 disp('Feasible')
elseif diagnostics.problem == 1
 disp('Infeasible')
else
 disp('Something else happened')
end

4:double

通过上面函数获得最优化结果后,需要显示,这时就需要value(X),double是其中一种,即从决策向量里抽取数值。

整体效果可以再以一个例子说明:

已知非线性整数规划为:
Max z=x1^2+x2^2+3*x3^2+4*x4^2+2*x5^2-8*x1-2*x2-3*x3-x4-2*x5
s.t.
0<=xi<=99(i=1,2,...,5)
x1+x2+x3+x4+x5<=400
x1+2*x2+2*x3+x4+6*x5<=800
2*x1+x2+6*x3<=800
x3+x4+5*x5<=200
在matlab中输入 x=intvar(1,5);
f=[1 1 3 4 2]*(x'.^2)-[8 2 3 1 2]*x';F=set(0<=x<=99);
F=F+set([1 1 1 1 1]*x'<=400)+set([1 2 2 1 6]*x'<=800)+set(2*x(1)+x(2)+6*x(3)<=800);
F=F+set(x(3)+x(4)+5*x(5)<=200);solvesdp(F,-f)
double(f)    80199
double(x)    53    99    99    99     0
intvar(m,n):生成整数型变量;
sdpvar(m,n):生产变量;
solvesdp(F,f):求解最优解(最小值),其中F为约束条件(用set连接),f为目标函数
double:显示求解的答案