[bisosprob] Update s-polynomials from constraints

[renatoloureiro]

The current way of solving bilinear sum-of-squares problems implies that the s-polynomials need to be specified with the monomials of a certain degree, this may be problematic because if these are chosen incorrectly without any concern for the constraints and the polynomials in it, most likely the solver will indicate that the problem is unfeasible or the algorithm will take a lot more time. This issue was created to address this usability problem and enable the user to set an a priori check of the s-polynomials and update them in case the initial specified degree/set of monomials is incorrect.

[tcunis]

This issue assumes that there are s-polynomials in the bilinear problem statement, or other multipliers with a similar role (certificates for correctness rather than decision variables that describe, e.g., a region of attraction. One way to approach the issue, in my opinion, would be to have a set-type constraint as syntactic sugar, viz.

prob.addsetinclusion(p, q, ...)

with the meaning "{x | p(x) ≥ 0} is included in {x | q(x) ≥ 0} which automatically creates a new decision variable s and encodes the constraint prob.ge(q, s*p, ...). The user then may or may not specify monomials for s. Such a constraint should be able to handle intersections as well as level sets with equality, {x | h(x) = 0}.

Another option (not necessarily alternative) would be to have the user to declare a decision variable as "certificate" (or similar) in the sense that the toolbox should select the monomials for it.

[tcunis]

As for the degree of the s-polynomial, if you haven't done so, please have a look at Appendix A.1.1. and A.1.3 of Tan's thesis [1]. There also is a selection of monomials implemented in sosopt's pcontain:

% Set default multiplier
if isempty(z)
    % This is roughly based on guidelines from Appendix A.1.1 
    % of the Ph.D. Thesis by Weehong Tan
    zmind = ceil(p1.mindeg/2);
    zmaxd = ceil( (p1.maxdeg-p2.maxdeg) / 2 );    
    if zmaxd<zmind
        z=1;
    else
        z = monomials(x, zmind:zmaxd );
    end
end

[1]: Tan, W. (2006). Nonlinear Control Analysis and Synthesis using Sum-of-Squares Programming. PhD thesis, University of California, Berkeley, CA.

[renatoloureiro]

Unfortunately I already implemented my version for the degree of the s-polynomials, i just did a push of it to the newly created branch spolydeg. Right now it may be a bit nasty due to nested loops and so on but I will try to fix that in the next couple of hours and add more comments.

[renatoloureiro]

prob.addsetinclusion(p, q, ...)

That doesn't seem too complicated, I will work on it then.

[renatoloureiro]

prob.addsetinclusion(p, q, ...)

With the current developments I'm already able to get the degree of p and q, so the only thing i should do is create a new variable s# and add a new constraint le for example. The problem is that all the elements in p are bilinear for sure due to the additional s# but the elements of q might not be entirely linear. Question: Should I work on a function to determine which variables are linear and bilinear in a certain expression or for the moment the user is mandated to specify the elements in q that are bilinear?

[tcunis]

To clarify, the inputs to addsetinclusion should be polynomials that are linear in the decision variables. With the multipliers from the generalized S-procedure (aka sufficient condition of the P-Satz) the full constraint would be bilinear. Thus, the user would currently need to enumerate the decision variables when calling this new function; something like:

function [prob,s] = addsetinclusion(prob, p, q, ..., pvar, qvar)
% Adds constraint that {x | q(x) ≥ 0} is subset of {x | p(x) ≥ 0}, 
% where pvar and qvar are cell arrays of variable names in p and q, respectively.

    ...

    bvar = [qvar svar]; % svar is cell array of multiplier variable names 

    prob = ge(prob, p - s*r, 0, pvar, bvar);
    prob = ge(prob, s, 0, svar);
end

[renatoloureiro]

With the latest commit https://github.com/tcunis/bisosprob/commit/484448f590212321b3a3b5fb968ee714cc3a1d71 the function addsetinclusion is now working