sostutorial.m

%% Manual model
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)

F = [coefficients(p-p_sos,[x y]) == 0, Q >= 0];
optimize(F,[],sdpsettings('dualize',1))

%% Manual model to compute lower bound
sdpvar t
F = [coefficients((p-t)-p_sos,[x y]) == 0, Q >= 0];
optimize(F,-t)

%% Manual model to compute matrix SOS
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 using solvesos
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)
clean(p-h'*h,1e-6)

%% check fidelity
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)
checkset(F)
e = checkset(F(is(F,'sos')))

[sol,v,Q,res] = solvesos(F);
res
  

%% Include parametric variables, here to compute lower bound
sdpvar x y lower
p = (1+x*y)^2-x*y+(1-y)^2;
F = sos(p-lower);
solvesos(F,-lower,[],lower);
value(lower)

solvesos(F,-lower);
value(lower)

%% Design polynomials
sdpvar x y t
p1 = t*(1+x*y)^2-x*y+(1-y)^2;
p2 = (1-x*y)^2+x*y+t*(1+y)^2;
F = [sos(p1), sos(p2)];
solvesos(F,t);
sdisplay(sosd(F(1)))
sdisplay(sosd(F(2)))

%% Constrained polynomial programming via multipliers
sdpvar x y lower
p = (1+x*y)^2-x*y+(1-y)^2;
g = [1-x;1+x;1-y;1+y]
[s1,c1] = polynomial([x y],2);
[s2,c2] = polynomial([x y],2);
[s3,c3] = polynomial([x y],2);
[s4,c4] = polynomial([x y],2);

F = [sos(p-lower-[s1 s2 s3 s4]*g), sos(s1), sos(s2), sos(s3), sos(s4)];
solvesos(F,-lower,[],[c1;c2;c3;c4;lower]);