8.4.1 リカッチ方程式と最適レギュレータ


 関数 "care" によるリカッチ方程式の求解と最適レギュレータ ...... M ファイル optimal_care.m (p.174)
clear
format compact

A = [  0  1
     -10 -1 ];
B = [ 0
      1 ];

Q = diag([ 300  60 ]);
R = 1;

P = care(A,B,Q,R)

K = - inv(R)*B'*P
>> optimal_care
P =
  170.0000   10.0000
   10.0000    8.0000
K =
  -10.0000   -8.0000

 関数 "lqr" による最適レギュレータの設計 ...... M ファイル optimal_lqr.m (p.175)
clear
format compact

A = [  0  1
     -10 -1 ];
B = [ 0
      1 ];

Q = diag([ 300  60 ]);
R = 1;

K = - lqr(A,B,Q,R)
>> optimal_lqr
K =
  -10.0000   -8.0000

 リカッチ方程式の求解(数式処理) ...... M ファイル symb_riccati.m (p.175)
clear
format compact

A = [  0  1
     -10 -1 ];
B = [ 0
      1 ];

Q = diag([ 300  60 ]);
R = 1;

syms p11 p12 p22

P = [ p11  p12
      p12  p22 ];
M = P*A + A'*P - P*B*inv(R)*B'*P + Q;

[p11 p12 p22] = solve(M(1,1),M(1,2),M(2,2), 'p11,p12,p22');

num = size(p11,1);
for i = 1:num
  P = [ p11(i)  p12(i)
        p12(i)  p22(i) ]
  lambda = double(eig(P))
  
  if real(lambda) > 0 & imag(lambda) == 0
    disp('   ====> P は正定行列');
    K = - inv(R)*B'*P;
    K = double(K)
  else
    disp('   ====> P は正定行列ではない');
  end
  disp(' ');
end
>> symb_riccati
P =
[ -30, -30]
[ -30,   0]
lambda =
  -48.5410
   18.5410
   ====> P は正定行列ではない
 
P =
[  10, -30]
[ -30,  -2]
lambda =
  -26.5941
   34.5941
   ====> P は正定行列ではない
 
P =
[ 170, 10]
[  10,  8]
lambda =
    7.3851
  170.6149
   ====> P は正定行列
K =
   -10    -8
 
P =
[ -190,  10]
[   10, -10]
lambda =
 -190.5539
   -9.4461
   ====> P は正定行列ではない
 

 Simulink を利用した最適レギュレータ制御のシミュレーションを行うためのファイル

 (1) 関数 "lqr" による最適レギュレータの設計 ...... M ファイル optimal_lqr.m (p.175)
clear
format compact

A = [  0  1
     -10 -1 ];
B = [ 0
      1 ];

Q = diag([ 300  60 ]);
R = 1;

K = - lqr(A,B,Q,R)

 (2) シミュレーション ...... Simulink モデル (p.90)
simulink_sfbk_R2010b.mdl(R2010b 以降のバージョン用)
simulink_sfbk_R2007b.mdl(R2007b 以降のバージョン用)
simulink_sfbk_R13SP1.mdl(R13SP1 以降のバージョン用)
 (3) シミュレーション結果の表示 ...... M ファイル plot_data_optimal.m (p.176)
figure(1);
plot(t,x(:,1));
grid;
xlabel('time [s]');
ylabel('x1 [m]');

figure(2);
plot(t,x(:,2));
grid;
xlabel('t [s]');
ylabel('x2 [m/s]');
 (3') シミュレーション結果の表示(グラフのカスタマイズ)
    ...... M ファイル plot_data_optimal2.m (本には記載していない)
% =======================================================
figure(1);
plot(t,x(:,1),'r','linewidth',2); 
grid;
 
% --- グラフのカスタマイズ --------------------------------
xlim([0 4]);
ylim([-0.5 1.5]);
 
set(gca,'xtick',[0:1:4]);
set(gca,'ytick',[-0.5:0.5:1.5]);
 
set(gca,'fontname','times','fontsize',20);
 
xlabel('{\it{t}} [s]','fontname','times','fontsize',22);
ylabel('{\it{x}}_{1}({\it{t}}) [m]','fontname','times','fontsize',22);
 
 
% =======================================================
figure(2);
plot(t,x(:,2),'r','linewidth',2); 
grid;
 
% --- グラフのカスタマイズ --------------------------------
xlim([0 4]);
ylim([-2.4 0.8]);
 
set(gca,'xtick',[0:1:4]);
set(gca,'ytick',[-2.4:0.8:0.8]);
 
set(gca,'fontname','times','fontsize',20);
 
xlabel('{\it{t}} [s]','fontname','times','fontsize',22);
ylabel('{\it{x}}_{2}({\it{t}}) [m/s]','fontname','times','fontsize',22);


 実行結果

 (1) 最適レギュレータの設計と初期状態 x(0) の設定
>> optimal_lqr
K =
  -10.0000   -8.0000
>> x0 = [ 1  0 ]';
 (2) シミュレーション開始

    もしくは
    もしくは
>> sim('simulink_sfbk_R2010b')
 (3) シミュレーション結果の表示
>> plot_data_optimal

 (3') シミュレーション結果の表示(グラフのカスタマイズ)
>> plot_data_optimal2