【OR】YALMIP Sum of squares（SOS）规划

sos programming

YALMIP内置了求解sos问题的模块，在sos问题中，需要对多项式$p(x)$进行分解，使得$p(x)=h^T(x)h(x)$，或者分解为$p(x)=v^T(x)Qv(x)$，其中$Q$是半正定矩阵，$v(x)$是向量.

YALMIP also supports linearly and nonlinearly parameterized problems, decomposition of matrix valued polynomials.

手动实现sos的建模求解过程

%% 定义多项式
x = sdpvar(1, 1);
y = sdpvar(1, 1);
p = (1+x)^4+(1-y)^2;

% 分解
v = monolist([x y],degree(p)/2);
Q = sdpvar(length(v));
p_sos = v'*Q*v;

%  系数匹配
F = [coefficients(p-p_sos,[x y]) == 0, Q >= 0];
% optimize(F)

% 设置从对偶形式求解
optimize(F, [], sdpsettings('dualize', 1))

找到多项式下界

%% 找到多项式下界
sdpvar t
F = [coefficients((p-t)-p_sos,[x y]) == 0, Q >= 0];
optimize(F,-t)

整合优化流程如下

sdpvar x y
P = [1+x^2 -x+y+x^2;-x+y+x^2 2*x^2-2*x*y+y^2];
m = size(P,1);
v = monolist([x y],degree(P)/2);
Q = sdpvar(length(v)*m);
R = kron(eye(m),v)'*Q*kron(eye(m),v)-P;
s = coefficients(R(find(triu(R))),[x y]);
optimize([Q >= 0, s==0]);
sdisplay(clean(kron(eye(m),v)'*value(Q)*kron(eye(m),v),1e-6))

sos optimization

使用命令sos定义（SOS约束），使用命令solvesos求解，sosd得到sos-decomposition结果

x = sdpvar(1,1);y = sdpvar(1,1);
p = (1+x)^4 + (1-y)^2;
F = sos(p);
solvesos(F);
h = sosd(F);
sdisplay(h)

检测结果，如果分解结果正确，那么$p(x)-h^T(x)h(x)$的值应该接近于0，可以使用命令clean忽略系数很小的向量

clean(p-h'*h, 1e-6)

得到分解式$p(x)=v^T(x)Qv(x)$

x = sdpvar(1,1);y = sdpvar(1,1);
p = (1+x)^4 + (1-y)^2;
F = sos(p);
[sol,v,Q] = solvesos(F);
clean(p-v{1}'*Q{1}*v{1},1e-6)

可以使用check指令返回系数最大值或者使用solvesos的第四个参数等效

checkset(F)
e = checkset(F(is(F, 'sos')))
[sol,v,Q,res] = solvesos(F);
res

References

pre and post processing sum-of-squares program in practice
Sum-of-squares programming