linearizebmi & solvebmi

Introduction, Download, Manuals, Examples, About Overbounding Approximation.      [ English / Japanese ]

Sample Programs

Some sample files are attached under demo folder of the distributed zip-file. For numerical examples, the coefficient matrices of the space state representations are selected from COMPleib (Problem 'HE1'), which is referred below. The plant data can be described by following program. In addition, the SDP solver SeDuMi is required to run the example programs. 
      %%% Define problem
      % Define plant data
      a=[-0.0366,  0.0271,  0.0188, -0.4555;...
          0.0482, -1.01,    0.0024, -4.0208;...
          0.1002,  0.3681, -0.7070,  1.42;...
          0,       0,       1,       0];
      b1=[4.678e-2, 0; 4.572e-2, 9.88e-3; 4.369e-2, 1.11e-3; -2.179e-2, 0];
      b2=[0.4422 0.1761;3.5446 -7.5922;-5.52 4.49;0 0];
      c1=(1/sqrt(2))*[2,0,0,0;0,1,0,0];
      c2=[0 1 0 0];
      nx=size(a,1);
      nw=size(b1,2);
      nu=size(b2,2);
      nz=size(c1,1);
      ny=size(c2,1);
      d11=zeros(nz,nw);
      d12=(1/sqrt(2))*[1, 0; 0, 1];
      d21=zeros(ny,nw);

linearizebmi

The following example program shows the way to use linearizebmi to convert BMI into LMI, and implement the ovebounding approximation method.
H∞ control problem for static output feedback are used as a BMI problem.

H∞ control problem

    %%% SDP solver settings (YALMIP)
    opts=sdpsettings;
    opts.solver='sedumi';   % SDP solver
    opts.verbose=0;

    %%% Define problem
    % Define decision matrices
    p=sdpvar(nx,nx,'symmetric');    % Lyapunov matrix
    k=sdpvar(nu,ny,'full');         % Controller (static gain)
    g=sdpvar(1,1);                  % H-inf norm


    % Define feasible solutions for 
    % each decision matrix in a bilinear term.
    p0=sdpvar(nx,nx,'symmetric');
    k0=sdpvar(nu,ny,'full');


    % Define decomposition matrix
    G =sdpvar(size(k,1));   % decision matrix
    G0=eye(size(G));        % feasible solution


    % BMI as a string
    Fstr = "[p*(a+b2*k*c2)+(p*(a+b2*k*c2))',p*(b1+b2*k*d21),(c1+d12*k*c2)';"   +...
            "(p*(b1+b2*k*d21))',            -g*eye(nw,nw),     (d11+d12*k*d21)';" +...
            "c1+d12*k*c2,                   d11+d12*k*d21,  -g*eye(nz)]";

        
    %%% Use linearizebmi() to convert BMI into dilated LMI
    %%% Select one of the dilation types.

    % (case 1) G is a constant matrix
    [LMIauto,LMIstr]=linearizebmi(Fstr,{'p','k'},{'p0','k0'});
    
    % (case 2) G is a decision matrix
    % [LMIauto,LMIstr]=linearizebmi(Fstr,{'p','k','G'},{'p0','k0','G0'});

    % Define constraints
    LMI=[LMIauto<=0];


    %%% Overbounding Aprroximation Method
    % Initial feasible solutions
    p0init=zeros(size(p));
    k0init=zeros(size(k));
    G0init=eye(size(G));

    % sequentially optimization
    lcmax=200;	% maximum step numbers
    ggall=[];
    for lc=1:lcmax
      extLMI=LMI;
      extLMI=replace(extLMI,p0,p0init);
      extLMI=replace(extLMI,k0,k0init);
      extLMI=replace(extLMI,k0,k0init);

      optimize(extLMI,g,opts);

      p0init=double(p);
      k0init=double(k);

      % Print out achieved optimal value
      gg=double(g);
      ggall=[ggall,gg];
      fprintf('Loop#%03d: %9.4f\n',lc,gg);

    end

    % Figure achieved value g
    figure;
    plot(ggall);
    figure;
    semilogy(ggall);
[Back to Top]

solvebmi

The following example programs solve BMI problems using solvebmi. H∞ and H2 control problems for static output feedback are used as BMI problems.

H∞ control problem

      %%% Define decision matrices
      p=sdpvar(nx,nx,'symmetric');    % Lyapunov matrix
      k=sdpvar(nu,ny,'full');         % Controller (static gain)
      g=sdpvar(1,1);                  % H-inf norm

      % Assign initial values to sdpvar
      assign(p,eye(nx,nx))
      assign(k,zeros(nu,ny))
      assign(g,0)


      %%% YALMIP optimize settings
      yalmipopts=sdpsettings;
      yalmipopts.solver='sedumi'; % SDP solver
      yalmipopts.verbose=0;       % No details

      %%% solvebmi options:
      opts = solvebmiOptions;
      opts.yalmip = yalmipopts;   % yalmip options
      opts.lcmax = 200;           % Maximum step times
      opts.regterm = 1;           % Use reguralation terms


      %%% BMI as a string
      Fstr = "[p*(a+b2*k*c2)+(p*(a+b2*k*c2))',p*(b1+b2*k*d21),(c1+d12*k*c2)';"   +...
              "(p*(b1+b2*k*d21))',            -g*eye(nw,nw),     (d11+d12*k*d21)';" +...
              "c1+d12*k*c2,                   d11+d12*k*d21,  -g*eye(nz)]";

          
      % All constraints (describe matrix variables)
      % In solvebmi, "Fstr" can be define [Fstr<=0]
      Flist = {Fstr, "-p"};


      %%% Use solvebmi() to execute ovebounding aproximation method 
      %%% to solve the BMI proble
      %%% Use 2 dilation type to solve BMI

      % (Type 1) Decomposition matrix G is a constant matrix
      % :Solve speed is fast
      opts.method = 0;
      [gg,vars,output] = solvebmi(Flist,{'p','k'},g,opts);

      % (Type 2) Decomposition matrix G is a decision matrix
      % :Converges quickly
      opts.method = 1;
      [gg2,vars2,output2] = solvebmi(Flist,{'p','k'},g,opts);


      %%% Output
      % Final acheived optimal value
      gg
      gg2

      % Designed controller
      K1 = vars.k
      K2 = vars2.k
[Back to Top]

H2 control problem

      %%% Define decision matrices
      p=sdpvar(nx,nx,'symmetric');    % Lyapunov matrix
      k=sdpvar(nu,ny,'full');         % Controller (static gain)
      g=sdpvar(1,1);                  % H-inf norm
      R=sdpvar(nz,nz,'symmetric');    % H2norm > sqrt(trace(R))

      % Assign initial values to sdpvar
      assign(p,eye(nx,nx))
      assign(k,zeros(nu,ny))
      assign(g,0)
      assign(R,eye(nz,nz))


      %%% YALMIP optimize settings
      yalmipopts=sdpsettings;
      yalmipopts.solver='sedumi'; % SDP solver
      yalmipopts.verbose=0;       % No details

      %%% solvebmi options:
      opts = solvebmiOptions;     % default option settings
      opts.yalmip = yalmipopts;   % yalmip options
      opts.lcmax = 200;           % Maximum step times


      %%% BMI as a string
      Fstr1 = "[p*(a+b2*k*c2)+(p*(a+b2*k*c2))',p*(b1+b2*k*d21);"+...
              "(p*(b1+b2*k*d21))',            -eye(nw,nw)];";

      % others constraints
      Fstr2 = "trace(R)-g";
      Fstr3 = "[-R              -(c1+d12*k*c2);"+...
              "-(c1+d12*k*c2)'  -p];";
          
      % All constraints (describe matrix variables)
      % In solvebmi, "Fstr" can be define [Fstr<=0]
      Flist = {Fstr1, Fstr2, Fstr3};


      %%% Use solvebmi() to execute ovebounding aproximation method 
      %%% to solve the BMI proble
      %%% Use 2 dilation type to solve BMI

      % (Type 1) Decomposition matrix G is a constant matrix
      % :Solve speed is fast
      opts.dilate = 0;
      [gg,vars,output] = solvebmi(Flist,{'p','k'},g,opts);

      % (Type 2) Decomposition matrix G is a decision matrix
      % :Converges quickly
      opts.dilate = 1;
      [gg2,vars2,output2] = solvebmi(Flist,{'p','k'},g,opts);


      %%% output
      % Final acheived optimal value
      gg
      gg2

      % Designed controller
      K1 = vars.k
      K2 = vars2.k
[Back to Top]
Return to Sebe's home page
sebe[a]ics.kyutech.ac.jp
Last modified: Mon Mar 7 12:42:01 JST 2022