Synthesizing SOS Controllers

[Joel]

Hi,

I would like to use the SOS controller example at:  http://www.control.auc.dk/~sboe02/public_html/tools/yalmip/htmldata/sos.htm#matrix
as a template for synthesizing SOS controllers.  If I define the controller "K" as a vector with arbitrary values then the code runs.  If I leave the controller "K" defined as a variable then the code does not run.  I have PENBMI via TOMLAB and have requested a copy from the developers as noted in the "Yalmip and Penbmi" thread.  Is this issue due to the solver?  Is there a way to solve BMI problems without PENBMI?

Thanks!
-Joel

The code is below:

clear all
yalmip('clear');

% States...
sdpvar x1 x2
x = [x1;x2];

% Non-quadratic Lyapunov z'Pz 
z = [x1;x2;x1^2];
P = sdpvar(3,3);
V = z'*P*z;

% Non-linear state feedback
v = [x1;x2;x1^2];
K = sdpvar(1,3);
u = K*v;

% System x' = f(x)+Bu
f = [1.5*x1^2-0.5*x1^3-x2; 3*x1-x2];
B = [0;1];

% Closed loop system, u = Kv
fc = f+B*K*v;

% Stability and performance constraint dVdt < -x'x-u'u
% NOTE : This polynomial is bilinear in P and K
F = set(sos(-jacobian(V,x)*fc-x'*x-u'*u));

% P is positive definite, bound P and K for numerical reasons
F = F + set(P>0)+set(25>P(:)>-25)+set(25>K>-25);

% Minimize trace(P)
% Parametric variables P and K automatically detected
% by YALMIP since they are both constrained
solvesos(F,trace(P));


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

double(K)

ans =

   NaN   NaN   NaN


Output with K = -0.0001*ones(1, 3) and F = F + set(P > 0)
-------------------------------------------------------------------------
YALMIP SOS module started...
-------------------------------------------------------------------------
Detected 6 parametric variables and 2 independent variables.
Detected 0 linear inequalities, 0 equality constraints and 1 LMIs.
Using kernel representation (options.sos.model=1).
Initially 7 monomials in R^2
Newton polytope (0 LPs).........Keeping 4 monomials (0.15625sec)
Finding symmetries..............Found no symmetries (0sec)
 
The coefficient matrix is not full row rank, numerical problems may occur.
SeDuMi 1.3 by AdvOL, 2005-2008 and Jos F. Sturm, 1998-2003.
Alg = 2: xz-corrector, theta = 0.250, beta = 0.500
eqs m = 11, order n = 8, dim = 26, blocks = 3
nnz(A) = 34 + 0, nnz(ADA) = 121, nnz(L) = 66
 it :     b*y       gap    delta  rate   t/tP*  t/tD*   feas cg cg  prec
  0 :            5.98E+000 0.000
  1 :  6.38E-001 1.81E+000 0.000 0.3032 0.9000 0.9000   1.52  1  1  1.7E+000
  2 :  4.45E+000 6.17E-001 0.000 0.3400 0.9000 0.9000  -0.12  1  1  1.4E+000
  3 :  1.23E+001 1.50E-001 0.000 0.2425 0.9000 0.9000  -0.28  1  1  8.8E-001
  4 :  3.60E+001 4.52E-002 0.000 0.3020 0.9000 0.9000  -0.54  1  1  8.1E-001
  5 :  3.20E+002 1.42E-002 0.000 0.3148 0.9000 0.9000  -1.26  1  1  2.6E+000
  6 :  7.02E+002 2.88E-003 0.000 0.2026 0.9000 0.9000  -3.11  1  1  1.1E+000
  7 :  7.07E+005 1.46E-005 0.000 0.0051 0.9990 0.9990  -1.15  1  1  5.4E+000
  8 :  1.07E+012 2.47E-012 0.000 0.0000 1.0000 1.0000  -1.00  1  1  
Primal infeasible, dual improving direction found.
iter seconds  |Ax|    [Ay]_+     |x|       |y|
  8      0.4  0.0e+000  0.0e+000  0.0e+000  7.4e+000

Detailed timing (sec)
   Pre          IPM          Post
1.570E-001    2.810E-001    3.100E-002    
Max-norms: ||b||=2.253167e-001, ||c|| = 1,
Cholesky |add|=0, |skip| = 1, ||L.L|| = 3.65754.
 
-> Solver reported unboundness of the dual problem.
-> Your SOS problem is probably infeasible.

[Johan Löfberg]

Yes, the reason is that you do not have a YALMIP interfaced BMI solver.

I gave up supporting TOMLABs PENBMI some years ago, when they decided to have their own slightly changed version of interface to PENBMI. Hence, YALMIP assumes you have the binaries directly from the developer of PENOPT instead. Their version is called penbmim.xxx, whereas TOMLABs is called penbmi.xxx. You could try to change the name on the mex to see if they are compatible now (doubt it)

Having said that, you should not bet on a method which requires a BMI solver to solver your problems. PENBMI happens to be lucky on this particular example and actually finds (*a*) solution to the problem, but this is not a generic property, it would more often fail I guess.

Finally, note that the code is old and slightly obsoleted. It would now be written as

F = sos(-jacobian(V,x)*fc-x'*x-u'*u);
F = [F, P>=0, 25>=P(:)>=-25, 25>=K>=-25];

[Joel]

Thank you for your help.  The .MEX file from TOMLAB appears to have the penbmi filename rather then penbmim. 
I'll keep trying to find a way to solve the problem though :-)

[Johan Löfberg]

Yes, that is the file I meant you should try to rename to penbmim instead. Probably doesn't work, but worth a try. 

I might take a look at the interface soon, to see if they converged enough to make it simple to interface TOMLABs interface (the PENOPT penbmim interface was horribly complex to work with, that is why I don't want to touch these things too much)

[Joel]

I renamed the file to penbmim.  The code ran but did not find a good solution (the controller is NaN):

-------------------------------------------------------------------------
YALMIP SOS module started...
-------------------------------------------------------------------------

Detected 6 parametric variables and 5 independent variables.

Detected 0 linear inequalities, 0 equality constraints and 1 LMIs.
Using kernel representation (options.sos.model=1).

Initially 32 monomials in R^5
Newton polytope (0 LPs).........Keeping 7 monomials (0.26563sec)

Finding symmetries..............Found no symmetries (0sec)
 
The coefficient matrix is not full row rank, numerical problems may occur.
SeDuMi 1.3 by AdvOL, 2005-2008 and Jos F. Sturm, 1998-2003.
Alg = 2: xz-corrector, theta = 0.250, beta = 0.500

eqs m = 29, order n = 11, dim = 59, blocks = 3
nnz(A) = 58 + 0, nnz(ADA) = 805, nnz(L) = 417

 it :     b*y       gap    delta  rate   t/tP*  t/tD*   feas cg cg  prec

  0 :            4.41E+000 0.000
  1 :  3.70E+000 1.31E+000 0.000 0.2981 0.9000 0.9000   0.71  1  1  2.5E+000
  2 :  5.32E+001 4.21E-001 0.000 0.3205 0.9000 0.9000  -1.14  1  1  7.3E+000
  3 :  7.05E+002 2.91E-002 0.000 0.0692 0.9900 0.9900  -2.12  1  1  5.2E+000
  4 :  4.74E+005 1.70E-005 0.000 0.0006 0.9999 0.9999  -1.09  1  1  2.0E+000
  5 :  4.74E+012 1.70E-012 0.000 0.0000 1.0000 1.0000  -1.00  1  1 

Primal infeasible, dual improving direction found.
iter seconds  |Ax|    [Ay]_+     |x|       |y|

  5      0.5  0.0e+000  0.0e+000  0.0e+000  2.0e+000

Detailed timing (sec)
   Pre          IPM          Post

4.530E-001    4.840E-001    7.800E-002   
Max-norms: ||b||=4.369320e-001, ||c|| = 1,
Cholesky |add|=0, |skip| = 1, ||L.L|| = 3.45311.

 
-> Solver reported unboundness of the dual problem.
-> Your SOS problem is probably infeasible.

[Joel]

I would like to get the YALMIP solver output, submit it to the PENBMI solver on the NEOS server (http://www.neos-server.org/neos/), and do any post processing in YALMIP.

The code I used is:
[model,recoverymodel] = export(F, trace(P), sdpsettings('solver','PENBMI'));

The output is:
Error using
lmi/categorizeproblem (line
267)
Report bug in problem
classification (polynomial
constraint)

Error in compileinterfacedata
(line 242)
[ProblemClass,integer_variables,binary_variables,parametric_variables,uncertain_variables,semicont_variables,quad_info]
=
categorizeproblem(F,logdetStruct,h,options.relax,parametric,evaluation_based,
Error in export (line 126)
    [interfacedata,recoverdata,solver,diagnostic,F]
    =
    compileinterfacedata(F,[],logdetStruct,h,options,findallsolvers);



I would appreciate any suggestions.

Thanks!

-Joel

[Johan Löfberg]

In order to export, you first have to tell YALMIP to flatten the SOS problem to an SDP (this should be don automatically, I will try to fix that...)

[F,h] = compilesos(F,h,ops,parameters)
model = export(F,h,ops)

Then, if you want to upload that to neos

pen = model.penstruct
save uploadtoneos pen

I tried that with the bilinear example you posted in the first post, which runs fin on my machine (PENBMI 2.1), but it crashed on the neos server.

As I said, don't bet your research on penbmi...

[Johan Löfberg]

FYI, the crash on NEOS is probably not a YALMIP issue. NEOS/PENBMI crashes already on the test-examples they have for penbmi. I will report this to the maintainers of NEOS.

[Johan Löfberg]

The NEOS maintainers have fixed the bug. I tested to upload the pen structure created above, and it solved the problem, computing the same solution as when calling PENBMI directly from YALMIP.

[Joel]

Great!  Thanks for letting  me know.  I tried it and it works for me too.  

You mentioned that the SOS controller synthesized with YALMIP/PENBMI may not be the best.  What do you recommend for work in this area?  I've been looking into Gloptipoly and SOSTOOLS.  Are these good for synthesizing SOS controllers and analyzing their robustness?

[Johan Löfberg]

I would say synthesizing controllers based on nonconvex SOS problems is a shaky idea to begin with. Gloptipoly and SOSTOOLS do not support nonconvex SOS and similar extension but only solve standard problems (corresponding to standard linearly parameterized problems in solvemoment and solvesos)

If you really really want to go down this avenue, you probably have to come up with your own methods to solve the resulting problems, or at least guide the search for a solver such as penbmi by supplying clever initial guesses and good cuts. A generic blind approach will probably not work.