SOS with polynomial variable

[DirkPitt]

Hello, i have some huge problem understanding how to solve the following problem in SOS - yalmip (i read the wiki). Basically i have to insert in the problem some other polynomials variable, but i can't figure out how.

This is the problem i have to solve

V is a polynomial function V(x,u), umax is a bound over u. ( dim[x]=n, dim[u]=m ). P is given.

     min                                tr(Q)
V,s1,s2,s3                  
                              s.t.      V is sos(n+m)
                                         s_i is sos(n+m) , i=1,2,3
                                         - ( (1-V)*s1 + (umax - u'u)*s2+ (1-V)*(umax-u'*u)*s10 + [x' u'] P [x;u]-1) is sos(n+m)

Q should be the Gram matrix in V=z' Q z

I did something like this (without optimization):
x=sdpvar(n,1);
u=sdpvar(m,1);
umax %given
%define V, s_i
p=-((1-V)*s1 + (umax-u'*u)*s2 + (1-V)*(umax-u'*u)*s3 + [x' u']*P*[x;u]-1);
v=[sos(V), sos(s1), sos(s2), sos(s3), sos(p)];
[sol,v,Sol]=solvesos(v);

but i can't define the polynomial variable V and s_i with the function 'polynomial' (when i tried such as V=polynomial([x;u],n+m) i receive an error) . Can you suggest me how to do it?

Thanks in advance!

[Johan Löfberg]

Since you are using the Gramian in one of the decompositions, you must encode that SOS condition manually (since those things are created inside solvesos and thus not accsesible outside solvesos). Is it really necessary for you to minimize the trace of the Gramian of one of the conditions?

% Manual setup of sos(V)
d=ceil(degree(V)/2);
z = [x;u];
v = monolist(z,d);
Q = sdpvar(length(v));
Model = [Q >= 0, coefficients(v'*Q*v - V,z)==0];

% Quadratic multipliers
[s1,c1] = polynomial(z,2);
[s2,c2] = polynomial(z,2);
[s3,c3] = polynomial(z,2);

% Other SOS constraints where we don't touch the Gramians
p=-((1-V)*s1+(umax-u'*u)*s2+(1-V)*(umax-u'*u)*s3+[x' u']*P*[x;u]-1);
Model = [Model, sos(s1),sos(s2),sos(s3),sos(p)];
parameters = [Q(:);c1;c2;c3];

solvesos(Model,trace(Q),[],parameters)

[Johan Löfberg]

Now I saw that V is not given. Then you cannot formulate the problem as a convex SOS problem since you want to multiply V with s1, both parameterized.

You have to fix either V or s1. If you fix s1 to, e.g., 1,  you get

v = monolist(z,desireddegree);
Q = sdpvar(length(v));
V = v'*Q*v;
Model = [Q >= 0];

% Quadratic multipliers
s1 = 1;


[s2,c2] = polynomial(z,2);
[s3,c3] = polynomial(z,2);

% Other SOS constraints where we don't touch the Gramians
p=-((1-V)*s1+(umax-u'*u)*s2+(1-V)*(umax-u'*u)*s3+[x' u']*P*[x;u]-1);


Model = [Model,sos(s2),sos(s3),sos(p)];
parameters = [Q(:);c2;c3];

solvesos(Model,trace(Q),[],parameters)


and if you fix Q you solve

v = monolist(z,desireddegree);
V = v'*Q*v;



% Quadratic multipliers
[s1,c1] = polynomial(z,2);


[s2,c2] = polynomial(z,2);
[s3,c3] = polynomial(z,2);

% Other SOS constraints where we don't touch the Gramians
p=-((1-V)*s1+(umax-u'*u)*s2+(1-V)*(umax-u'*u)*s3+[x' u']*P*[x;u]-1);


Model = [sos(s1),sos(s2),sos(s3),sos(p)];
parameters = [c1;c2;c3];

solvesos(Model,[],[],parameters)


[DirkPitt]

Yes, i have to minimize a certain ellipsoid....

Maybe this will be a poor questions, but i need some further explanation: how can i define a-priori V (so that i can take the degree in the first line of your code) ? And how did you decide the degree of s_i woudl be 2?

Really thanks

[DirkPitt]

Il giorno martedì 21 maggio 2013 14:11:43 UTC+2, Johan Löfberg ha scritto:

Now I saw that V is not given. Then you cannot formulate the problem as a convex SOS problem since you want to multiply V with s1, both parameterized.

So it's maybe a bilinear problem?

[Johan Löfberg]

I assumed V was a given polynomial, hence I could apply the degree function.

The choice of quadratic multipliers was just an example.

[Johan Löfberg]

Yes, if V is unknown and parametrized by Q, and s1 is an unknown polynomial parametrized by c1, it is bilinear in Q and c1

[DirkPitt]

Il giorno martedì 21 maggio 2013 14:20:26 UTC+2, Johan Löfberg ha scritto:

I assumed V was a given polynomial, hence I could apply the degree function.

The choice of quadratic multipliers was just an example.

Ok me too tried with just an example but growing the degree of multipliers the program start running and...doesn't finish :(

So if I'm correct in my idea (V is not given) i can't use yalmip?

[Johan Löfberg]

Increasing the degree will lead to polynomially larger SDP, so it quickly grows to big

Bilinear SOS is basically not solvable. Optimistically, one would say just setup the bilinear SOS in YALMIP and tell YALMIP to solve the problem using the bilinear SDP solver penbmi. Does not work in practice though. You will have to come up with clever iterative heuristic approaches to attack the problem, and for that you can of course use YALMIP. The problem you are trying to solve (region of atraction analysis?) has been addressed many times using such approaches.

[DirkPitt]

I just re-checked the paper (IEEE) i am following. Even there it's written the problem is bilinear! But it doesn't suggest any heuristic.

For sake of clarity the problem is one step controllable set for non linear system

[Johan Löfberg]

Which paper is it.

[DirkPitt]

Mmm, i feel i will get some problem with my degree saying that :(

However you are cited on it:

"Condition (xx)-(xx) have a bilinear structure, therefore the BMI Optimization will be numerically solved by means of standard DSP packages, e.g. [27],[28]"

And [27] is: "J. Lofberg "A toolbox for modeling and optimization in matlab"

[DirkPitt]

SDP, not DSP of course!

[DirkPitt]

Why did you suggest penbmi? Is bmibnb not able in this case?

[Johan Löfberg]

Sure, the *global* bmisolver bmibnb is of course also applicable. It would however require the *local* bmi solver penbmi to solve node problems

A completely theoretical discussion of course, because in practice the problems would be too large for bmibnb, and already the local solver penbmi would have problems too (based on my experience)

[DirkPitt]

Ok, thank you so much for your help and explanations!

[Johan Löfberg]

A very optimistic statement. Or at least, it hides a massive amount of work to actually implement something which works reasonably well. Basically, do some heuristic approach, iterate of s1 and V, start local solver penbmi from good initial points, and add constraints based on things you know from simulating the system 

Related work has been parts of several papers and theses by, e.g., W. Tan, Ufuk Topcu, Andrew Packard, Jarvis-Wloszek etc etc, just to mention one cluster of researchers

http://www.cds.caltech.edu/~utopcu/papers/TPSW-acc07.pdf

[DirkPitt]

Yes, i knew something about it. But at this point i suppose have no other possibility from asking who give me the before mentioned paper.

[DirkPitt]

I tried to solve it just to see what happened, but now i can't undestand if the error was in my own code, or not. I tried with bmibnb and fmincon and this is the output

-------------------------------------------------------------------------
YALMIP SOS module started...
-------------------------------------------------------------------------
Detected 136 parametric variables and 7 independent variables.
Detected 0 linear inequalities, 0 equality constraints and 1 LMIs.
Using image representation (options.sos.model=2). Nonlinear parameterization found
Initially 8 monomials in R^7
Newton polytope (0 LPs).........Keeping 8 monomials (0sec)
Finding symmetries..............Found no symmetries (0sec)
Initially 115 monomials in R^7
Newton polytope (1 LPs).........Keeping 85 monomials (0.0468sec)
Finding symmetries..............Found no symmetries (0.0624sec)
 

ans =

    solvertime: 0
          info: 'No suitable solver'
       problem: -2
    yalmiptime: 107.4900

[Johan Löfberg]

To solve bilinear SDPs using the global solver, a nonlinear SDP solver is required. You do not have any (penbmi is the only possible choice)