SOS problem: access to coefficients

[Denis Aßmann]

Dear all,

my goal is find a polynomial p(x) = sum_a p_a * x^a which is non-negative on some basic semi-algebraic subset B = {x | g_i(x) >= 0 for all i} but satisfies some additional constraints over its coefficients:  

min c^T p_a

p_a = 0   for some a \in A

p(x) = sum_a p_a * x^a  \in Q(g_i) = s_0 + g_1 * s_1 + ... + g_k * s_k

How can this kind of problem be modeled with YALMIP?

I am currently working with a small toy example where B = [-1, 1], p(x) = a + b x + c x^2 and c = [ 2.5 -0.625 1.45833] (which is actually the integral of p(x) between -1.5 and 1).

Furthermore, p(x) is bounded from above by 1, so we have

min Integral_[-1.5, 1] p(x)

st.

0 <= p(x) <= 1 for all x

deg(p) = 2 [i.e. p(x) = a + bx + c x^2]

My attempt to express this with YALMIP:

x = sdpvar(1, 1);

lb = 1 + x;  % x >= -1

ub = 1 - x;  % x <= 1  <=>  1 - x >= 0

% quadratic module: p0 + p1 * lb + p2 * ub

[p0, c0, v0] = polynomial(x, 2);

[p1, c1, v1] = polynomial(x, 0);

[p2, c2, v2] = polynomial(x, 0);

% unknown polynomial p

[p, c, v] = polynomial(symbols, 2);

F = [coefficients(p - p0 + p1 * lb + p2 * ub, [x]) == 0, sos(p0), sos(p1), sos(p2), sos(1-p)];

obj = [2.5 -0.6 1.45] * coefficients(p, x);

[sol,v,Q] = solvesos(F, obj, options, [c; c0; c1; c2]);

The model is dual infeasible and it is also missing the part where I want to set p_a = 0 for some subset of the coefficients. However, since I implemented this model with SOSTOOLS as well, I am pretty sure that it should have a feasible solution (so my yalmip attempt is wrong).

Ideally, I'd prefer to not introduce this ``p'' polynomial at all but rather access the coefficients of p0 + p1 * lb + p2 * ub directly...

To conclude, I am wondering how to access coefficients of unknown polynomials in a SOS-type problem with yalmip.

Any help is much appreciated!

Thank you very much and best regards,

Denis 

[Johan Löfberg]

Undefined function or variable 'symbols'.

 

[Denis Aßmann]

Sorry, I failed to test the code / clear matlab after simplifying it a bit. This should work:

x = sdpvar(1, 1);

lb = 1 + x;  % x >= -1

ub = 1 - x;  % x <= 1  <=>  1 - x >= 0

% quadratic module: p0 + p1 * lb + p2 * ub

[p0, c0, v0] = polynomial(x, 2);

[p1, c1, v1] = polynomial(x, 0);

[p2, c2, v2] = polynomial(x, 0);

% unknown polynomial p

[p, c, v] = polynomial(x, 2);

F = [coefficients(p - p0 + p1 * lb + p2 * ub, [x]) == 0, sos(p0), sos(p1), sos(p2), sos(1-p)];

obj = [2.5 -0.6 1.45] * coefficients(p, x);

[sol,v,Q] = solvesos(F, obj, [], [c; c0; c1; c2]); 

[Johan Löfberg]

The problem is feasible, but unbounded (simple solve the feasibility problem and you see that it is feasible).

...or by adding a limit on objective, and see that you always meet that limit

[sol,v,Q] = solvesos([F,obj >= -1e4], obj, [], [c; c0; c1; c2])
value(coefficients(p,x))

BTW, why not simply

F = [sos(p + p1 * lb + p2 * ub), sos(p1), sos(p2), sos(1-p)];

and