linearizebmi & solvebmi

はじめに, ダウンロード, マニュアル, サンプルプログラム, 逐次LMI化法. 　　[ English / Japanese ]

サンプルプログラム

いくつかのサンプルファイルが配布している zip-file の demo フォルダ以下に含まれている. ここでは
S. N. Singh and A. A. R. Coelho,
"Nonlinear control of mismatched uncertain linear systems and application to control of aircraft",
Journal of Dynamic Systems, Measurement and Control, 106, pp.203-210, 1984.
(これはベンチマーク問題集 COMPLieb 1.1 に "HE1" として収録されている数値例である.) の制御対象に対する以下の問題のサンプルプログラムを示す.

linearizebmi
- H_∞ control problem
solvebmi
- H_∞ 制御問題
- H₂ 制御問題

%%% SDP solver の設定
opts=sdpsettings;
opts.solver='sedumi';    % 利用する SDP solver
opts.verbose=0;

%%% 問題の定義
% 制御対象データの定義
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

定数出力フィードバックによる H_∞ 制御問題は, BMI 問題として定式化される. これを linearizebmi を用いて逐次 LMI 化法によって解くサンプルプログラムを以下に示す.

% 決定変数の定義
p=sdpvar(nx,nx,'symmetric');    % Lyapunov 行列
k=sdpvar(nu,ny,'full');         % 制御器(定数ゲイン)
g=sdpvar(1,1);                  % H∞ノルム

% 双線形項決定変数の暫定解の宣言 (dummy)
p0=sdpvar(nx,nx,'symmetric');
k0=sdpvar(nu,ny,'full');

% BMI 最適化問題の定義
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),      (d11+d12*k*d21)';" +...
       "c1+d12*k*c2,                   d11+d12*k*d21,   -g*eye(nz)]";

%%% 提案した linearizebmi() による BMI の十分条件化
[LMIauto,LMIstr]=linearizebmi(Fstr,{'p','k'},{'p0','k0'});
LMI=[LMIauto<=0];

%%% 逐次 LMI 化法
% 初期暫定解の宣言
p0init=zeros(size(p));
k0init=zeros(size(k));

lcmax=200;      % 繰り返し回数の設定
ggall=[];       % 達成値の更新過程保存用
for lc=1:lcmax
  extLMI=LMI;
  extLMI=replace(extLMI,p0,p0init);
  extLMI=replace(extLMI,k0,k0init);

  optimize(extLMI,g,opts);

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

  % 達成値表示のためのコード
  gg=double(g);
  ggall=[ggall,gg];
  fprintf('Loop#%03d: %9.4f\n',lc,gg);

end

% 達成値の更新過程の表示
figure;
plot(ggall);
figure;
semilogy(ggall);

上記サンプルプログラムのダウンロード

solvebmi

linearizebmiを利用した， BMI問題を自動的に解く solvebmiのサンプルプログラムを以下に示す．先の例題に対して定数フィードバックによる H_∞ 制御問題と H₂ 制御問題のサンプルプログラムです．

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]

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
      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]

瀬部昇のホームに戻る

Last modified: Mon Mar 7 12:41:49 JST 2022