the choice of basis in sos

[Yalmiper]

What is the choice of polynomial basis in sos programming? It seems to be very unstable for my problem. When I change solver the solution can change a lot. Is there a way to change the choice of polynomial basis in Yalmip?

[Johan Löfberg]

No such options, so you would have to set it up manually as explained in the SOS tutorial on the Wiki

Care to share your model to see if anything can be done quickly to imrove things?

[Yalmiper]

My problem is a Lyapunov function construction problem, so I am dealing with parametrized polynomials. I can parametrized the candidate Lyapunov function with some predefined basis, but I think when I use sos(p), the decomposition z^TQz will still use standard monomials as z? Will my choice of basis on the left hand side make numerical differences? (I do not want do sos programming by hand because I want to use the exploiting sparsity feature in Yalmip sos function).

The code might be a little too complicated:

clear;
x = sdpvar(2,1);
%Linear State-Space Transformation: Numerical Scaling for Polynomial
%Optimizations
T = [1,0;0,100];
y = T*x;
recRuls = [1 0;-1 0;0 1;0 -1];
ktc = 0.01;
kta = 0.01;
kdx = 5e-5;
kdm = 0.02;
recRate = [ktc;kdm*y(1);kta*y(1);kdx*y(2)];

%define Lyaopunv function
mono = monolist(x,3);
monoq = replace(mono,x,x.^2);
P = sdpvar(length(mono));
V = mono'*P*mono;
rV = monoq'*P*monoq;
C = [P];
C = C + sos(rV);

%calculate generator
qV=0;
for r = 1:size(recRuls,1)
    ds = y + recRuls(r,:)';
    ds = T\ds;
    qV = qV + recRate(r)*(replace(V,x,ds) - V);
end

[coq,moq] = coefficients(qV,x);
moqq = replace(moq,x,x.^2);
qqV = coq'*moqq;

%define the second constraint
b = 0.01;
s = sdpvar(1); %scaling parameter
lhs = b*s-qqV;
C = C + [s>=0] + sos(lhs);

% %define the third constraint
% % %box boundary
B = [0,0;10,400];
B = B/T;
a = B(1,:)';
c = B(2,:)';
% g = [flipud(a)-flipud(x),c-x];
g = c-x.^2;

% %the controlled specious
ind = 2;
% %define approximation function for indicator function
pd = 12;
pmo = monolist(x(ind),pd/2);
AM = sdpvar(length(pmo));
pf = pmo'*AM*pmo;
pfq = replace(pf,x,x.^2);
h = [a(ind)-x(ind);x(ind)-c(ind)];

l = [];
for i=1:length(h)
    Li = sdpvar(length(pmo));
    C = C + [Li];
    l = [l, pmo'*Li*pmo];
end
C = C + sos(l) + sos(pf+l*h); %ensure p is positive on chosen box B

for i=1:length(g)
    lambda = sdpvar(1,length(g(i,:)));
    C = C + [lambda>=0] + sos(-(qqV+1)+lambda*g(i,:)');
end
C = C + sos(pf-1-s-qV);

% formulate object function (i.e. the VOLUME!)
[coe,monomal] = coefficients(pf,x);
intco = zeros(length(monomal),1);
for i = 1:length(monomal)
    monog = monomal(i);
    for j = 1:length(ind)
        k = ind(j);
        monog = int(monog,x(k),a(k),c(k));
    end
    intco(i) = monog;
end
obj = intco'*coe;
%setting and optimize
ops = sdpsettings('solver','sdpt3');
optimize(C,obj,ops);


%%evaluate
eP = value(P);
es = value(s);

syms x1 x2;
X = [x1;x2];
state = T*X;
monos = [1;x1;x2;x1^2;x1*x2;x2^2;x1^3;x1^2*x2;x1*x2^2;x2^3];
% monos = [1;x1;x2;x3;x1^2;x1*x2;x2^2;x1*x3;x2*x3;x3^2];
% monos = [1;state];
eP(1,1) = 0;

eV = (1/es)*(monos'*eP*monos);
recrate = [ktc;kdm*state(1);kta*state(1);kdx*state(2)];

eqV=0;
for r = 1:size(recRuls,1)
    ds = state + recRuls(r,:)';
    ds = T\ds;
    eqV = eqV + recrate(r)*(subs(eV,X,ds) - eV);
end
eqV = subs(eqV,X,T\X);

%plot contour to check
qv = matlabFunction(eqV);

[X1,X2] = meshgrid(0:19,0:399);
av = zeros(size(X1));
for i = 1:size(X1,1)
    for j = 1:size(X1,2)
        av(i,j) = qv(X1(i,j),X2(i,j));
    end
end
b = max(max(av))+1;
mv = zeros(size(X1));
for i = 1:size(X1,1)
    for j = 1:size(X1,2)
        if av(i,j)>-1
            theta = av(i,j)+1;
            mv(i,j) = (1-theta)/(b-theta);
            if mv(i,j)<0
                mv(i,j) = 0;
            end
        else
            mv(i,j) = 1;
        end
    end
end
figure;
contour(X1,X2,av,linspace(-1,b,10));
colorbar;
figure;
contour(X1,X2,mv,[0.01,0.1,0.3,0.5,0.8,0.9]);
colorbar;

The objective function is trying to approximately minimize the volume of the contour qV>=-1, while restrict the maximum value of qV. When I use "SDPT3" it performs very well, but "mosek" or "sedumi" gives a polynomial that have a contour with much larger volume.

[Johan Löfberg]

Your issue is most likely the fact that the polynomials you study have  coefficients spanning 12 orders of magnitude

plot(semilogy(sort(abs(getbase(lhs)))))

[Yalmiper]

All right. I tried to scale it better but apparently I failed. Thank you for the advices. So if the coefficients are very badly scaled then it tends to have such problems?

[Yalmiper]

Actually I am optimizing over parametrized polynomials with variables x, you are plotting numerical coefficients of all the variables. Can you explain more what lead to a better numerical behaviour? I should try to scale these numerical coefficients to a smaller range of magnitude so the sdp program of these decision variables (i.e. the parameters of polynomals) will have a better numerical behaviour?

[Johan Löfberg]

with such numerics you can expect anything and nothing

[Johan Löfberg]

Good scaling of polynomials is a research problem on its own. Sum-of-squares is simply ill-conditioned in many cases

[Yalmiper]

Dear Professor Loefberg,

I don't quite understand it. So I am requiring the parametrized polynomials with variables x to be sos. Thus it generate a sdp program with decision variables of the parameters of the polynomial. So in general I need to keep the numerical coefficients of the decision variables of the sdp program in a smaller region of magitude? What if the optimal solutions really has values span a large range of magnitudes?

Thank you!

[Johan Löfberg]

Then you have problems. Sum-of-squares doesn't magically just always work. On the contrary actually...

[Yalmiper]

That's just very sad. Anyway, thanks for the useful advices.

[Yalmiper]

Dear Professor Loefberg,

Just to make to sure. I was asking how to change the choice of basis when using sos in Yalmip. If I want to use sos function in Yalmip, it means the right hand side of sos program (p = z^tQz), I am using standard monomials as z. So the only freedom I have is try different ways of parameterize p. Is that correct? (I probably would not construct sos program by hand as I want to use YALMIP's sos function to exploit sparsity, which to my knowledge is the only parser that implemented both Newton Polytope and Symmetry).

Thank you!

[Johan Löfberg]

standard monomial basis only possibility in solvesos