最近在做论文时,涉及到最优化问题,而最优化里面很多时候涉及的是二次约束二次规划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:显示求解的答案