State feedback controller gain design

[Ahmed El Ebiary]

Hello 

 I am trying to use YALMIP and Sedumi solver to design the gains of the primary controller used for DC DC converter as in the attached paper page 7  but this should be done by solving an optimization problem as in the attached image but sedumi always gives me "Dual infeasible, primal improving direction found." So i need to know where is the problem 

Also the first LMI constraint has arbitrary entries i need to know how to define these arbitrary entries

 Thanks 

Best Regards

[Johan Löfberg]

https://yalmip.github.io/debugginginfeasible

[Ahmed El Ebiary]

Thank you for your reply Prof.Johan

I tried the steps in the debugging infeasible page, but still give me infeasible solution, also i checked the constraints and found that one of the constraints has primal negative residual as in the first screenshot which is M1Y as in the code below and the gain of the controller obtained is too small as in the second screenshot and doesn't give the desired reference tracking. 

I have also attached the problem and the code .

Thanks 

Best Regards

%%The Code 

%Parameters

  Vs1 = 100; Vs2 = 100 ; 

  Vref1 = 48 ; Vref2 = 48 ; 

  R12 = 0.05 ;R21 = 0.05; % TL Resistances definition

  L12 = 2.10e-06;

  Ct1 =0.0022 ;Ct2 =  0.0022;   %Capacitors definitions

  Lt1 = 0.0018; Lt2 = 0.002; % Inductors definitions

  Rt1 = 0.2; Rt2 = 0.3 ;  % Resistance definitions

  

% Primary controller for DG1 system matrices

  T1 = -(1/(R12*Ct1));

  A11 = [T1  1/Ct1  0 ; -1/Lt1  -Rt1/Lt1 0 ; -1 0 0];

  B1 = [0  ; 1/Lt1   ;0  ];

  T2 = -(1/(R21*Ct2));

  A22 = [T2  1/Ct2  0 ; -1/Lt2  -Rt2/Lt2  0 ; -1 0 0];

  B2 = [0  ; 1/Lt2  ;0 ];

%% Primary Controller Design DGU1

% Define Decision variables
Y1 = sdpvar(3,3,'full');
G1Y =  sdpvar(1,3);
gama11 = sdpvar(1,1);
gama21 = sdpvar(1,1);
gama31 = sdpvar(1,1);
beta1 = sdpvar(1,1);
zeta1 = sdpvar(1,1);

% Define positive weights
alpha11 = 10^-6;
alpha21 = 10^-2;
alpha31 =10^-3;
alpha41 =10^-4;
alpha51 =10^-5;
% Define Constraints
Taw1 = [ gama11 0 0 ; 0 gama21 0; 0 0 gama31];
Y1 = [((10*Ct1)^-1) 0 0; 0 0.1 0.2 ; 0 0.4 0.2];
M11 = (Y1*transpose(A11)) + (A11*Y1)+(transpose(G1Y)*transpose(B1))+(B1*G1Y);
M1Y = [ M11 Y1 ; Y1 -Taw1];  % This is the violated Constraint 
M21 = [-beta1*eye(3,3)  transpose(G1Y)];
M22 = [G1Y -eye(1,1)];
M2Y =  [M21 ; M22]; 
M3Y = [Y1 eye(3); eye(3) zeta1*eye(3)];
const1 = [Y1>=0,M1Y<=0,M2Y<=0,M3Y>=0,gama11>=0,gama21>=0,gama31>=0,beta1>=0,zeta1>=0];
  % Define objective function       
 h1 = alpha11*gama11 + alpha21*gama21 + alpha31*gama31 + alpha41*beta1 + alpha51*zeta1;
 ops = sdpsettings('solver','sedumi');
 optimize(const1,h1,ops);
 G1Y = double(G1Y);
 Y1 = double(Y1);
 K1 = G1Y * inv(Y1);

Thanks

[Johan Löfberg]

You've missed transpose

>> M1Y

Linear matrix variable 6x6 (full, real, 6 variables)

Coeffiecient range: 0.1 to 826446.281

>> M1Y<=0

++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

|   ID|                     Constraint|   Coefficient range|

++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

|   #1|   Element-wise inequality 20x1|   0.1 to 25207.0707|

++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

and/or you probably intended Y1 to be symmetric too which it isn't

>> norm(Y1-Y1')

ans =

     2.000000000000000e-01

and the code generally doesn't make sense as Y1 is constant, but you add a constraint on it??

const1 = [Y1>=0,M1Y<=0

generally it looks like you don't display at the expressions and constraint you are creating

BTW, transpose(X) is a very cumbersome way of writing X'

[Johan Löfberg]

and the correct command is value, double is obsolete since a decade back

onsdag 16 juni 2021 kl. 10:50:15 UTC+2 skrev ahmedha...@gmail.com:

[Ahmed El Ebiary]

Y1 in the paper is the first constraint with bottom right matrix as an arbitrary entries so does this means that these arbitrary entries should be a sdpvar in the code .

Also i displayed the constraints and found that three out of  four constraints are displayed as element wise inequalities and they should be matrix inequalities as defined in the optimization problem 

|   #1|    Element-wise inequality 4x1|              1 to 1|           % Y1 after using sdpvar for arbitrary entries

|   #2|   Element-wise inequality 20x1|     1 to 25252.5253| % M1Y

|   #3|          Matrix inequality 4x4|              1 to 1|                    % M2Y

|   #4|    Element-wise inequality 7x1|              1 to 1|              % M3Y

[Johan Löfberg]

The weird notation is probably short for blkdiag(1/eta,sdpvar(2))

You are still using an unsymmetric Y1 (I guess) as M1Y has turned out unsymmetric 

[Johan Löfberg]

and M3Y too

[Ahmed El Ebiary]

Is there a condition that all variables should be symmetric as I am new to YALMIP, also Prof. You told me that ''   You've missed transpose ''  

>> M1Y

Linear matrix variable 6x6 (full, real, 6 variables)

Coeffiecient range: 0.1 to 826446.281

I didn't understand where i have missed it. 

Thanks 

[Johan Löfberg]

If you create Y1 symmetric,

>>Y1 = blkdiag(((10*Ct1)^-1),sdpvar(2))

Linear matrix variable 3x3 (symmetric, real, 3 variables)

Coeffiecient range: 1 to 45.4545

then M1Y will be symmetric (since M11 is symmetric)

K>> M1Y

Linear matrix variable 6x6 (symmetric, real, 9 variables)

Coeffiecient range: 1 to 826446.281

had you defined Y1 without symmetry

K>> Y1 = blkdiag(((10*Ct1)^-1),sdpvar(2,2,'full'))

Linear matrix variable 3x3 (full, real, 4 variables)

Coeffiecient range: 1 to 45.4545

then M1Y would not be symmetric

K>> M1Y = [ M11 Y1 ; Y1 -Taw1]

Linear matrix variable 6x6 (full, real, 12 variables)

Coeffiecient range: 1 to 826446.281

as it would need a transpose on Y1 in that case to give symmetry

K>> M1Y = [ M11 Y1' ; Y1 -Taw1]

Linear matrix variable 6x6 (symmetric, real, 12 variables)

Coeffiecient range: 1 to 826446.281

however an Y1 without symmetry is obviously not what you want as you have Y1 in a diagonal block in M3Y

[Ahmed El Ebiary]

Thank You Prof. I will try this solution Thanks a lot for your time and support.

[Ahmed El Ebiary]

  Good Morning Prof Johan

The solution works and becomes feasible but the second constraint gives -ve primal residual when i check the constraint as shown in capture 2, also I tried solving another optimization problem of the same type as the one in capture 1 and the same happened,  and down below is the code. 

And can you please till me if the dimension of LMI constraint M2Y is true to be 4x4 or it should be 6x6 knowing that G1Y is 1x3 vector 

Thanks Prof for your help and time.

%% the code

Line dependent approach for the primary controller design

  % TL Resistances definition

  R12 = 0.05 ;R21 = 0.05;

  R13 = 0.07;
  R14 = 0.03;

   
  L12 = 2.10e-06;
  Ct1 =0.0022 ;Ct2 =  0.0022;   %Capacitors definitions
  Lt1 = 0.0018; Lt2 = 0.002; % Inductors definitions

  Rt1 = 0.2; Rt2 = 0.3 ;  % Resistance defintions

 
% Primary controller for DG1 system matrices
  T1 = -(1/(R12*Ct1))  
  A11 = [T1  1/Ct1  0 ; -1/Lt1  -Rt1/Lt1 0 ; -1 0 0];
  B1 = [0  ; 1/Lt1   ;0  ];
  T2 = -(1/(R21*Ct2));
  A22 = [T2  1/Ct2  0 ; -1/Lt2  -Rt2/Lt2  0 ; -1 0 0];
  B2 = [0  ; 1/Lt2  ;0 ];

 
%% Primary Controller Design DGU1

% Define unknown variables
yalmip('clear');
Y1 = sdpvar(3,3);
G1Y =  sdpvar(1,3)
gama1 = sdpvar(1,1)
beta1 = sdpvar(1,1)
delta1 = sdpvar(1,1)

% Define positive weights
alpha11 = 10^-6;
alpha21 = 10^-2;
alpha31 =10^-3;

%Constraints definition
Y1 = blkdiag(inv(10^-4),sdpvar(2))
M11 = (Y1*(A11')) + (A11*Y1)+((G1Y')*(B1'))+(B1*G1Y)
M1Y = [ M11 Y1 ; Y1 -gama1*eye(3)]
M21 = [-beta1*eye(3,3) (G1Y')]
M22 = [G1Y -eye(1,1)]
M2Y =  [M21 ; M22] % error here most probably has dimension of 4x4
M3Y = [Y1 eye(3); eye(3) delta1*eye(3)]
const1 = [Y1>=0,M1Y<=0 ,M2Y<=0,M3Y>=0,gama1>=0,beta1>=0,delta1>=0]

           
%objective function

 h1 = alpha11*gama1 + alpha21*beta1 + alpha31*delta1
 ops = sdpsettings('solver', 'sedumi', 'sedumi.eps', 1e-8, ...
                'sedumi.cg.qprec', 1, 'sedumi.cg.maxiter', 49, ...
                'sedumi.stepdif', 2);
%ops = sdpsettings('solver','sedumi');
optimize(const1,h1,ops);
G1Y = value(G1Y)
Y1 = value(Y1)

K1 = G1Y * inv(Y1)

[Johan Löfberg]

both sedumi and mosek screams about numerics, and indeed your data looks nasty in one constraint

>> const1

+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

|   ID|                    Constraint|     Coefficient range|

+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

|   #1|         Matrix inequality 3x3|            1 to 10000|

|   #2|         Matrix inequality 6x6|   1 to 181818181.8182|

|   #3|         Matrix inequality 4x4|                1 to 1|

|   #4|         Matrix inequality 6x6|            1 to 10000|

|   #5|   Element-wise inequality 1x1|                1 to 1|

|   #6|   Element-wise inequality 1x1|                1 to 1|

|   #7|   Element-wise inequality 1x1|                1 to 1|

+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

[Ahmed El Ebiary]

Does that means that i should define the constraints in a better way and  there is something wrong in the code

[Johan Löfberg]

it just means you have large numbers (Y1 is on the order of 10^4 and so is A11, so you will get stuff in the order of 10^8). Bad theory/data which needs improvements

[Piotr Balik]

It appears that state transition matrices are indeed badly conditioned cond(A11)=inf, rank(A11)=2.

Could that be the cause? Here I just removed the states to test

% Primary controller for DG1 system matrices
T1 = -(1/(R12*Ct1));

A11 = [T1  1/Ct1  ; -1/Lt1  -Rt1/Lt1 ; -1 0 0];

B1 = [0  ; 1/Lt1   ;0  ];
T2 = -(1/(R21*Ct2));
A22 = [T2  1/Ct2  0 ; -1/Lt2  -Rt2/Lt2  0 ; -1 0 0];
B2 = [0  ; 1/Lt2  ;0 ];

N=2; %state removal
A11 = A11(1:N,1:N);
A22 = A22(1:N,1:N);
B1=B1(1:N);
B2=B2(1:N);

%% Primary Controller Design DGU1
% Define unknown variables
yalmip('clear');

Y1 = sdpvar(N,N);
G1Y =  sdpvar(1,N)

gama1 = sdpvar(1,1)
beta1 = sdpvar(1,1)
delta1 = sdpvar(1,1)


% Define positive weights
alpha11 = 10^-6;
alpha21 = 10^-2;
alpha31 = 10^-3;

%Constraints definition

Y1 = blkdiag(inv(10^-4),sdpvar(N-1))

M11 = (Y1*(A11')) + (A11*Y1)+((G1Y')*(B1'))+(B1*G1Y)

M1Y = [ M11 Y1 ; Y1 -gama1*eye(N)]
M21 = [-beta1*eye(N) (G1Y')]

M22 = [G1Y -eye(1,1)]
M2Y =  [M21 ; M22] % error here most probably has dimension of 4x4

M3Y = [Y1 eye(N); eye(N) delta1*eye(N)]

const1 = [Y1>=0,M1Y<=0 ,M2Y<=0,M3Y>=0,gama1>=0,beta1>=0,delta1>=0]

%objective function
h1 = alpha11*gama1 + alpha21*beta1 + alpha31*delta1
ops = sdpsettings('solver', 'sedumi', 'sedumi.eps', 1e-8, ...
    'sedumi.cg.qprec', 1, 'sedumi.cg.maxiter', 49, ...
    'sedumi.stepdif', 2);
%ops = sdpsettings('solver','sedumi');
optimize(const1,h1,ops);
G1Y = value(G1Y)
Y1 = value(Y1)
K1 = G1Y * inv(Y1)

check(const1)
value(h1)

[Johan Löfberg]

those properties are completely uninteresting. 0 would be bad then...

the problem is the scale on numbers

https://yalmip.github.io/inside/debuggingnumerics/

[Ahmed El Ebiary]

I checked the constraints Professor and their coefficients range and they are as in the attached screenshots so it is the second constraints M1Y that causes the problem as it ranges from 1 to 181818181.8182 and it is also the one with negative prime residual, so how can i solve this problem can i scale that constraint as i checked the debugging solutions and it is suggesting that i am using wrong units but this is not the case.

Thanks Prof.

Best Regards

[Johan Löfberg]

There is no magic trick. Starting with a balanced realization never hurts. Where does the magic number 1/10^-4 come from. does it have to be that large. why not make it a decision variables and solve some regularized problem

[Ahmed El Ebiary]

1- I tried that solution professor and it didn't give the desired performance, also in the original paper the authors defined eta in Y1 and Y2 to be 10^-4 as in the attached feasibility and infeasibility screenshot.  

2- Also, is defining the arbitrary entries in Y1 and Y2 as sdpvar the proper solution or they should be constants?

3- Also I would like to know to judge the feasibility and infeasibility using Sedumi solver and YALMIP as the problem is solved telling me that it is a numerical problem but saying nothing about feasibility. 

4- Can you please Prof. Check my code and see if there is something wrong compared with the attached optimization problem required to be solved 

I have also attached the original paper for more information about the primary state feedback controller design that is formulated in the form of this optimization problem.

%% The Code 

clc , clear 

warning('off','YALMIP:strict');  yalmip('clear')

format shortG

  Vs1 = 100; Vs2 = 100; 

  Vref1 = 48 ; Vref2 = 48; 

  % TL Resistances definition

  R12 = 0.05 ;

  R13 = 0.07;

  R14 = 0.03;

  L12 = 2.10e-06;

  Ct1 =0.0022 ;Ct2 =  0.0022; Ct3 =  0.0022 ; Ct4 =  0.0022;  %Capacitors definitions

  Lt1 = 0.0018; Lt2 = 0.002; Lt3 = 0.0022; Lt4 = 0.003; % Inductors definitions

  Rt1 = 0.2; Rt2 = 0.3 ;Rt3 = 0.1; Rt4 = 0.5;  % Resistance defintions

  

% Primary controller for DG1 system matrices

  T1  = -(1/(R12*Ct1));  % for connecting DG1 with DG2 only and gives feasible when used in A11

  T12 = -((1/(R12*Ct1))+(1/(R13*Ct1))); % for connecting DG1 with DG2 & DG3  gives feasible solution when used in A11 instead of T1

  T13 = -((1/(R12*Ct1))+(1/(R13*Ct1))+(1/(R14*Ct1))); % for connecting DG1 with DG2 & DG3 & DG4 must give infeasible solution as per paper when used in A11 instead of T1

  A11 = [T1  1/Ct1  0 ; -1/Lt1  -Rt1/Lt1 0 ; -1 0 0]; % for Dg1 and Dg2 use T1 here for Dg1&2&3 use T12 for Dg1&2&3&4 use T13

  B1  = [0  ; 1/Lt1   ;0  ];

  

  T2  = -(1/(R12*Ct2));

  A22 = [T2  1/Ct2  0 ; -1/Lt2  -Rt2/Lt2  0 ; -1 0 0];

  B2  = [0  ; 1/Lt2  ;0 ];

  

 

%% Primary Controller Design DGU1

% Define unknown variables

G1    = sdpvar(1,3);

gama1  = sdpvar(1) ;

beta1  = sdpvar(1);

delta1 = sdpvar(1);

% Define positive weights

alpha11 = 10^-6;

alpha21 = 10^-2;

alpha31 = 10^-3;

% Define Constraints

Y1  = blkdiag(inv(10^-4),sdpvar(2));

M1 = [(Y1*(A11'))+(Y1*(A11'))'+(B1*G1)'+(B1*G1)             Y1;

      Y1'                                       -gama1*eye(3)]

M2 = [-beta1*eye(3,3)       G1';

       G1            -eye(1,1)]

M3 = [Y1               eye(3);

      eye(3)   delta1*eye(3)]

% comment const1 = [Y1>=0,trace(Y1) == 1,M1<=0 ,M2<=0,trace(M2) == 1,M3>=0,trace(M3) == 1,gama1>=0,trace(gama1) == 1,beta1>=0,trace(beta1) == 

% comment 1,delta1>=0,trace(delta1) == 1]% there is a problem in M1Y constrain

const1 = [ Y1>=0 , M1<=0 , M2<=0 , M3>=0,gama1>=0,beta1>=0,delta1>=0]% there is a problem in M1 constrain

            

% Objective Function Definition

Obj = alpha11*gama1 + alpha21*beta1 + alpha31*delta1;

% Optimization 

ops = sdpsettings('solver','sedumi','verbose',0);

Sol = optimize(const1,Obj,ops)

Info.Objective = value(Obj);

if Sol.problem == 0

    Info.Solution = 'Successfully solved';

    elseif Sol.problem == 4

    Info.Solution = 'Numerical problems';

    elseif Sol.problem == 1

    Info.Solution = 'Infeasible';

    elseif Sol.problem == -4

    Info.Solution = 'Solver not applicable';

end

Info.Nu_LMI_variables = length(depends(const1)); % Number of LMI variables

Info.LMI_Constraints = length(const1); % Number of LMI constraints

% Finding Controller Gains

Info.G1 = value(G1);

Info.Y1 = value(Y1);

Info.P1 = inv(Info.Y1);

Info.K1 = value(G1) * inv(Info.Y1);

Info.L1 = eig(A11+B1*Info.K1);

Info

Thanks 

Best Regards

[Piotr Balik]

Balancing helps reducing the order of coefficients:

%need to be corrected

C1=[1 1 1];

C2=C1;

sys=ss(A11,B1,C1,[])
[sysb,g]=balreal(sys);
[~,remove] = min(g);

sysr=modred(sysb,remove,'del')
A11=sysr.A;
B1=sysr.B;

sys=ss(A22,B2,C2,[])
[sysb,g]=balreal(sys);
[~,remove] = min(g);
sysr=modred(sysb,remove,'del')
A22=sysr.A;
B2=sysr.B;

Where we get the following result

+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
|   ID|                    Constraint|   Coefficient range|
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
|   #1|         Matrix inequality 4x4|        0.2 to 10000|
|   #2|         Matrix inequality 3x3|              1 to 1|
|   #3|         Matrix inequality 4x4|          1 to 10000|
|   #4|   Element-wise inequality 1x1|              1 to 1|

|   #5|   Element-wise inequality 1x1|              1 to 1|
|   #6|   Element-wise inequality 1x1|              1 to 1|

instead of the 1818...
And then, you'd have to transform the controller for reduced state to full object

[Ahmed El Ebiary]

This reduces the matrix A11 and A22 to 2x2 matrix and also the state vector will be 2x1 vector which is not the case required for the controller as the state feedback controller gains should be 1x3 vector