Alternating optimization in YALMIP

[Richard Southern]

Hi there,

First, let me say a big thank you for developing this fantastic tool!

I've got a question about Yalmip usage for SOS problems. I'm trying to evaluate the stability region of a simple Van Der Pol oscillator based on the method in here:
http://arxiv.org/abs/1010.2241

This essentially requires implementing 4 SOS constraints,2 of which are Lagrange multipliers. The solution to this problem is primal infeasible (run the code in the attachment to plot the constraints), but the paper on page 5, 2nd column, around the middle states:
"This optimization is bilinear in the decision parameters, however a reasonable approach which has proven successful is to iterate between optimizing over multipliers and optimizing over V."

So assuming I was going to develop a loop of some sort which alternatively solves for multipliers and then optimizes for V, how do I tell Yalmip to do this? Do I need to completely reformulate the problem?

This may be a lot to digest: if there are any examples of using Yalmip with alternating solutions that would probably be sufficient!

I am not a mathematician, so if my question is daft, then apologies.
 - R

[Richard Southern]

A quick follow up: is the answer essentially "replace"?

[Johan Löfberg]

Several issues with the code

rho is an sdpvar, and you divide by this. Hence, you model is definitely not polynomial, if you want rho to be one of the variables (which it is now, as you have only declared cL and cM] as parametric decision variables). If you want it to be parametric, the resulting problem would not be a linear SDP. Additionally, you use an objective which is x_perp, which immediately tells YALMIP that x_perp is a parametric decision variable. From the construction of your other code, that does not appear reasonable.

Ffrom the looks of it, you want OneOverRho (not rho, as it introduces a division), cL and cM to be parametric variables

As you want find the maximum of x_perp'*x_perp, you have to introduce another parametric decision variable t, and minimize t s.t. sos(t-xperp'*xperp)

[Richard Southern]

So here's a working version of the alternative solution, in case someone else has a similar query. We've still not quite mastered the problem, but we're getting feasible solutions in most of the regions.

If there are any other solutions to this I'd love to hear them.

Thanks again for your help!

    DV = 1.9102*xp^2-274.3428*rho*xp^2+3.5527e-14*rho*xp+0.7161*xp^3-244.6211*rho*xp^3-74.3883*rho*xp^4-7.9622*rho*xp^5;
    V = 10*xp^2;
    d = 3.820474944+1.4322*xp;
    delta = 1.0e-6;
   
    % Create some lagrange multipliers
    pdeg = 4;
    [L, cL, vL] = polynomial(xp, pdeg, 0);
    [M, cM, vM] = polynomial(xp, pdeg, 0);  
   
    % Set up our constraints     
    C1 = -DV - L*(1-V);
    C2 = d - delta - M*(1-V);
   
    next_rho = 1.0;
    this_rho = 10.0;
   
    tol = 1.0e-5;
    i = 0;
    while (abs(next_rho - this_rho) > tol)
        this_rho = next_rho;
        i = i+1;
       
        disp(strcat('step ', num2str(i), ', rho=',num2str(this_rho)));
       
        % First step: optimize for Legrange multiplier coefficients
        sos1 = sos(replace(C1,rho,this_rho));
        sos2 = sos(replace(C2,rho,this_rho));
        sos3 = sos(L);
        sos4 = sos(M);    
        coefList = [cL;cM];
        consList = [sos1;sos2;sos3;sos4];
        [sol, q, Q, res] = solvesos(consList,[], options, coefList);
       
        if (sol.problem ~= 0)
            disp(strcat('Error solving for lagrange multiplier coefficients at [',num2str(x_n),'] "',yalmiperror(sol.problem),'"'));
            continue;
        end
       
        % Assuming this worked, fix the coefficients and solve for rho
        sos1 = sos(replace(C1,cL,value(cL)));
        sos2 = sos(replace(C2,cM,value(cM)));
        sos3 = sos(replace(L,cL,value(cL)));
        sos4 = sos(replace(M,cM,value(cM)));
        coefList = rho;
        consList = [sos1;sos2;sos3;sos4];
        [sol, q, Q, res] = solvesos(consList,[], options, coefList);
       
        if (sol.problem ~= 0)
            disp(strcat('Error solving for rho at [',num2str(x_n),'] "',yalmiperror(sol.problem),'"'));
            continue;
        end
       
        next_rho = value(rho);
    end;

[Johan Löfberg]

Why don't you minimize rho in the second sos problem. It appears to me as if you want to find the smallest possible rho, but now you just iteratively solve feasibility problems and hope that rho decreases in each solve

[Johan Löfberg]

This would be an implementation based on the optimizer object

P = [sos(C1),sos(C2),sos(L),sos(M)];
LagrangeProblem = optimizer(P,[],[],rho,{cL,cM});
RhoProblem =  optimizer(P,rho,[],{cL,cM},rho);

while (abs(next_rho - this_rho) > tol)
    this_rho = next_rho;
    i = i+1;
   
    disp(strcat('step ', num2str(i), ', rho=',num2str(this_rho)));

   
    [csol,problem] = LagrangeProblem{this_rho}
   
    if (problem)

        disp(strcat('Error solving for lagrange multiplier coefficients at [',num2str(x_n),'] "',yalmiperror(sol.problem),'"'));
        continue;
    end

   
    [rhosol,problem] = RhoProblem{csol};
   
    if (problem)

        disp(strcat('Error solving for rho at [',num2str(x_n),'] "',yalmiperror(sol.problem),'"'));
        continue;
    end

   
    next_rho = rhosol
end

[Richard Southern]

egads you're right! You, sir, are a genius.

Switching second solvesos with this:
[sol, q, Q, res] = solvesos(consList,rho, options);

We still have some instability issues, but that will be due to choices of polynomial degrees and delta's.

Thanks!