Error solving BMI using PENLAB

[Cindy Cindy]

Hi Johan,

I saw your post about BMI (https://yalmip.github.io/example/decayrate/)

And I want to do this:

We can ignore the polytopic part for now.

And here is my code

clear all

clc

fs = 10000.0;

Ts = 1/fs;

fg = 60;

wg = 2*pi*fg;

R = 0.5;

L = 1.7e-3;

A=[R/L wg;

   -wg R/L];

B=[1/L 0;

   0 1/L];

C=eye(2);

D=zeros(2);

   

[Ade, Bde] = c2d(A, B, Ts);

% Integral model

Ai = zeros(2,2);

Bi = eye(2);

[Aide, Bide] = c2d(Ai, Bi, Ts);

AA = [Ade, zeros(2,2); -Bide*C, Aide];

BB = [Bde; zeros(2,2)];

% declare variables

Q = sdpvar(4,4);         

Y = sdpvar(2,4);

Qi = sdpvar(4,4);  

alpha = sdpvar(1,1);

% declare LMIs

LMI1 = Qi;                   

LMI2 = [Qi (AA*Qi+BB*Y)';

            AA*Qi+BB*Y Q];

LMI = [LMI1 >= 0, LMI2 >= 0, Q<(alpha*Qi), 0<alpha, alpha<1];

options = sdpsettings('solver','penlab','verbose',1,'warning',1);

solution = optimize(LMI,alpha,options);

K = double(Y)*inv(double(Qi));

When I run the code, sometime I got message "PenLab didn't converge: unconstrained minimization failed" or the Matlab just stuck at some iteration (for example iteration number 8) and no result comes out even after I wait for a long time. And then I need to stop the running code because the Matlab won't continue the simulation.

Do you have any idea why?

Thank you.

[Johan Löfberg]

Penlab is known to have crashing bugs, so don't use it.

Solve it by bisection instead

[Cindy Cindy]

Thank you, Johan.

So I tried to implement bisection as below and I can get the result.

alpha_lower = 0.0;

alpha_upper = 1.0;

tol = 0.01;

alpha_works = alpha_lower

while (alpha_upper-alpha_lower)>tol

  alpha = (alpha_upper+alpha_lower)/2;

  disp([alpha_lower alpha_upper alpha])

  LMI1 = Qi;

  LMI2 = [Qi (AA*Qi+BB*Y)';

              AA*Qi+BB*Y Q];

  F = [LMI1>=0, LMI2>=0, Q<(alpha*Qi), alpha>alpha_lower, alpha<alpha_upper];

  ops = sdpsettings('verbose',0,'warning',0);

  sol = optimize(F,[],ops);

  if sol.problem==1

    alpha_upper = alpha;

  else

    alpha_lower = alpha;

    alpha_works = alpha;

    Pworks = value(Qi);

 end

end

I used lower and upper bound of alpha as given (0<alpha<1).

I can get the gain, however when I implemented the gain in simulink (state feedback control simulation), the result has very large overshoot in the starting point.

As I know, we cannot tune the parameter for LMI/BMI, unllike LQR.

So, how can I decrease the overshoot?

Thank you.

[Johan Löfberg]

The poor performance has nothing to do with the fact that you've solved a BMI. You simply have a bad choice of control design method. As far as I can tell, you are simply deriving a fastest possible decay controller. If you don't like the behavior of the resulting closed-loop system, you've obviously not solved a relevant problem. Why don't you solve the LQR problem then if you like that controller better

[Cindy Cindy]

Okay I understand.

I derived the decay rate like that because I'd like to reproduce an approach refer to a published paper.

And I use LMI because I'd like to deal with parameter uncertainties later on.

Do you have any reference books/papers recommendation related to decay controller?

Thank you for your help.

[Johan Löfberg]

Doing LQR with polytopic system is simple, it's basically this https://yalmip.github.io/example/lpvstatefeedback/ except that you don't parametrize the controller

[Cindy Cindy]

Thank you for your response. I really appreciate it.

I tried to derived the equations for discrete time system.

Here is the implementation in Matlab:

Q = eye(4);

R = 1*eye(2);

Y = sdpvar(4,4);

L = sdpvar(2,4,'full');

n = 4;

nu = 2;

LMI1 = Y;

LMI2 = [ Y,                         L',                    (AA*Y+BB*L)',   Y;

              L,                         inv(R),             zeros(nu,n) ,      zeros(nu,n);

              (AA*Y+BB*L),      zeros(n,nu),    Y,                       zeros(n,n);

              Y,                         zeros(n,nu),    zeros(n,n),         inv(Q)]; 

LMI = [LMI1>=0, LMI2>=0];

optimize(LMI,-trace(Y))

K = value(L)*inv(value(Y));

In this case, the result is really unstable. And I already tried to tune the weighting matrix Q and R (I just use the basic value for example in this code).

Should I tune it until I get satisfy result? What do you think?

Thank you.

[Johan Löfberg]

Works here. Can only assume you have a bad soler, or made an incorrect simulation

>> -dlqr(AA,BB,Q,R)

ans =

   -1.5671   -0.0293    0.8273   -0.4741

    0.0293   -1.5671    0.4741    0.8273

>> value(L)*inv(value(Y))

ans =

   -1.5671   -0.0293    0.8272   -0.4740

    0.0293   -1.5671    0.4740    0.8272

[Cindy Cindy]

Ah okay, so there is no problem with the code.

Actually I have different system from what I wrote here (because here is the simple version). And the gain I got is as below

K =

  Columns 1 through 4

    3.9890   -0.1118  -10.4656    0.1688

    0.1118    3.9890   -0.1688  -10.4656

  Columns 5 through 8

    0.2809   -0.0066    0.5442    0.2624

    0.0066    0.2809   -0.2624    0.5442

  Columns 9 through 12

    3.9619    3.8654    0.2586    0.2523

   -0.2586   -0.2523    3.9619    3.8654

  Columns 13 through 16

    1.9909    1.9720    0.1151    0.1140

   -0.1151   -0.1140    1.9909    1.9720

Even tho the eig are all inside the unit circle, my system is still unstable. And as for my simulation, I believe it's correct because I used it before to implement LQR and the result was as expected.

However now I want to try using LMI and as you suggested, I used LQR-LMI based now for implementation. That's why I'm wondering maybe I should tune the weighting matrix Q and R?

[Johan Löfberg]

impossible to say anything. if the solver works well, and you see that eigenvalues are stable, the simulation should be stable, unless you have horrible data causing numerical issues in the simulation