構造固定化された制御器に対するロバストチューニング

はじめに

前回の記事にてH∞制御を用いたロバスト制御系の設計をトライしました。

このような制御は新規に制御システムを設計するといった際には適用しやすいですが、例えば既存/レガシーの制御構造があって、そこから大きく変更することができない場合において、パラメーターチューニングのみで制御器をロバストに構成しないといけないといった課題に直面することがあると思います。
また仮に新規に制御システムを構築できる場合においても制御器の次元の問題によってコントローラー（マイコンなどのハード）に実装できないという問題に直面することもあります。

特に産業界では、プラントのプロセス制御や航空機の飛行制御系、自動車の運動制御などにおいてはPID制御が多用されており、信頼実績のある制御構造が決まっていて、運転に支障をきたすような新規性はユーザーに好まれないといったケースがあります。

この場合には、構造が固定化された制御器に対し、性能を改善するような施策が必要となります。
特に制御器については、ロバスト性が重要視されるようなケースがほとんどだと思います。

今回はこのようなケースを題材にどのように取り組んでいけるのか考えてみたいと思います。

航空機の縦方向の制御

題材としては、前回の記事と同様にB747の線形化された縦方向の運動モデルを扱っていきたいと思います。

上記において、エレベーター舵角のアクチュエーターダイナミクスは一次遅れ（$T_{\delta_e}$）＋むだ時間（$DT$）として考慮しています。

ここで、航空機の飛行制御系は歴史的に下図に示すようなSAS（Stability Augumentation System）とCAS(Control Augumentation System)という制御構造をベースに設計されていることが一般的です。

ここで、SASとは制御対象を安定化させるためのフィードバックループであり、通常はピッチレートなどの変化率信号にゲインを掛けてフィードバックすること減衰性を向上させる効果を狙います。

一方のCASはパイロットの操縦性を向上させるための補償機構であり、積分器などで角度指令値への良好な追従性を示すように設計されます。

今回はこの構造にならって、機体の縦方向の制御が以下のような状態フィードバック+積分器で構成される古典的なサーボ系であると考えます。

この中で設計者が自由に変えられるのは制御ゲイン$K_x$ならびに$K_i$のみです。
この2つのゲインだけを用いて、所定の制御性能を保証する制御系を設計しなければなりません。

不確かさの考慮

前回の記事と同じく、縦方向の線形システムのノミナルなパラメーター値については、ジャンボジェットの愛称で知られるBoeing747の巡行形態時（高度40000ft(12192m)、Mach 0.8）における値を使用します。

ノミナル値

Xu = -0.0028;
Xa = 9.1771;
Xq = -18.9201;
Xth = -9.7784;
Zu = -0.0003;
Za = -0.3191;
Zq = 0.9813;
Zth = -0.0034;
Mu = 0.0007;
Ma = -0.7841;
Mq = -0.4271;
Mth = 0.0003;
Xde = 0.4389;
Zde = -0.0233;
Mde = -1.1579;
Tde = 1;       %アクチュエーターの時定数
DT = 0.2;      %アクチュエーターのむだ時間
d = 1;         %外乱の振幅ゲイン

上記のノミナル値に対し、各パラメーラーの変動量を以下の様に見積もります。

不確かさの設定

%公称値に対する割合による表現
Xu = ureal('Xu',Xu,'percent',50);
Xa = ureal('Xa',Xa,'percent',50);
Xq = ureal('Xq',Xq,'percent',50);
Zu = ureal('Zu',Zu,'percent',50);
Za = ureal('Za',Za,'percent',50);
Zq = ureal('Zq',Zq,'percent',50);
Zth = ureal('Zth',Zth,'percent',50);
Mu = ureal('Mu',Mu,'percent',50);
Ma = ureal('Ma',Ma,'percent',50);
Mq = ureal('Mq',Mq,'percent',50);
Mth = ureal('Mth',Mth,'percent',50);
d = ureal('d',d,'percent',50);

%公称値を含む区間による表現
Xde = ureal('Xde',Xde,'Range',[0.3 0.5]);
Zde = ureal('Zde',Zde,'Range',[-0.05 -0.01]);

%公称値を基準とした増減による表現
Mde = ureal('Mde',Mde,'PlusMinus',[-0.5 0.5]);    
Tde = ureal('Tde',Tde,'PlusMinus',[-0.5 0.5]); 
DT = ureal('DT',DT,'PlusMinus',[-0.05 0.05]); 

以上の値を基にして制御対象のモデルを構築します。

制御対象

%機体の縦方向モデル
A = [Xu Xa Xq Xth;
     Zu Za Zq Zth;
     Mu Ma Mq Mth;
     0 0 1 0];
B = [Xde;Zde;Mde;0];
C = [0 0 0 1];
D = 0;
P = ss(A,B,C,D);
P.InputName = 'de';
P.OutputName = 'theta';
Pnom = P.NominalValue;

%アクチュエーターモデル
Pact = tf(1,[Tde 1]);
Pact.InputName = 'de_c';
Pact.OutputName = 'de';
Pdt = tf([-0.5*DT 1],[0.5*DT 1]);   %むだ時間は1次のパデ近似で表現

%全体モデル
Pall = P*Pact*Pdt;
Pall_nom = Pall.NominalValue;

ノミナル制御器の設計

まずベースラインとなるノミナルの制御器を設計します。
設計にあたっては、前節で設計したアクチュエーターのダイナミクスを含んだ全体モデルのノミナル値をベースにLQRで設計します。

今回はピッチ角$\theta$に関してサーボ系を構成するため、以下のように目標値$r$との偏差に関する積分器を用意して、それを元の制御対象と統合した拡大系システムにします。

LQRの設計

nx = size(Pall_nom.A,1); nu = size(Pall_nom.B,2);
Aa = [Pall_nom.A zeros(nx,1);
     -Pall_nom.C 0];
Ba = [Pall_nom.B;zeros(1,nu)];

% 拡大系のシステムを離散化
ts = 1e-3;                                     % Sampling time
M = expm(ts*[Aa Ba;zeros(1,size(Aa,1)+1)]);    % Discretized system
Aa = M(1:nx+1,1:nx+1);
Ba = M(1:nx+1,end);

% Design discrete LQR
Wx = diag([1 1 1 1 1 1 1e3]); Wu = .1;
K = dlqr(Aa,Ba,Wx,Wu);

LQRの設計に際して、拡大系を制御周期$ts$で離散化します。
離散化には行列指数のexpmを使用しています。

離散化したシステムを用いてLQRを設計する際にはdlqrコマンドを使用することが出来ます。

以上で求めたゲインを初期値として、これから不確かさに対しロバストなパラメーター設計を行っていきたいと思います。

手始めにノミナル値でチューニング

ロバストチューニングのステップに入る前に、手始めに制御対象がノミナルな場合における制御パラメーターのチューニングを行ってみて、前節で設計したLQRのゲインからどの程度応答性が改善するかを確認してみたいと思います。

チューニングに際してはSimulinkモデルをベースに、Simulink Control Designから提供されているチューニングタスク用のAPIコマンド $slTuner$ならびに$systune$を使っていきます。

設計のベースとなるSimulinikモデルは以下の通りです。

上図においてFCSというサブシステムがコントローラーモデル部であり、内部は積分制御と状態フィードバックで構成されています。なお、状態量の推定にはKalman Filterを使っています。

ノミナルチューニング

% SimulinkモデルにLQRで設計した値を反映
set_param([mdl '/FCS/Kx'],'Gain','K(1:end-1)')                              % mdlはモデル名
set_param([mdl '/FCS/Ki'],'Gain','K(end)')

% Simulinkモデルからチューニングする情報を抽出
blkpath = {'Ki','Kx'};                                                      % チューニング対象のブロック名
options = slTunerOptions('RateConversionMethod','tustin','SampleTime',ts);  % 離散化オプション    
ST0 = slTuner(mdl,blkpath,options);                                         % mdlはモデル名 
ST0.Ts = ts;                                                                % Specify sampling time

% 調整目標（Tuning Goal）の設定
Req1 = TuningGoal.Margins('de_c',6,40);                                     % ゲイン余裕：6dB、位相余裕：40[deg]
Req2 = TuningGoal.Transient('r','y',tf(1,[10 1]));                          % 閉ループの過渡応答を時定数10秒の規範モデルで規定
Req3 = TuningGoal.StepRejection('d','y',2*pi/180,10);                       % ステップ状外乱から出力に対する影響を振幅で2[deg]以内、10秒で抑え込む
Req4 = TuningGoal.Tracking('r','y',10,0.02,1);                              % 目標値への追従性を応答時間10秒、定常誤差2%、ピーク値の誤差比率1以内に抑える
soft = [Req2 Req3];                                                         % 柔軟な目標
hard = [Req1 Req4];                                                         % 厳密な目標

% 調整開始
ST = systune(ST0,soft,hard);

上記がslTuner及びsystuneコマンドを利用したSimulinkモデルに対するゲインのチューニングプログラムです。
まずslTunerコマンドを呼び出して、チューニング対象となるSimulinkモデルの指定ならびにブロック等の情報をST0と名付けたslTunerオブジェクトに抽出します。

続いて、設計における指標である調整目標（Tuning Goal）を指定します。
指定できる調整目標は時間領域、周波数領域、安定余裕など様々あり、それをSimulinkで指定した線形化ポイント上の信号ラインと紐づけていきます。

今回は、制御系の安定余裕の確保、目標値に対する良好な過渡応答特性ならびに追従性、外乱抑圧性を確保するために調整目標を4つ指定しています。

指定した調整目標をST0と共にsystuneコマンドに渡します。
なお、systuneについては与える調整目標を柔軟な目標、厳密な目標として区別して指定することが出来るようになっています。
systuneは非凸（ノンスムース）な問題を扱うことが可能な最適化ソルバーであり、与えられる調整目標を制約条件として所定の条件を満たす設計変数を求めようとします。

この際に、柔軟な目標、厳密な目標とは、ソフト制約、ハード制約にそれぞれ該当します。
systuneでの調整時にそれぞれの指標がスカラー値として返りますが、値が1以下であれば制約を満足していることを示します。

今回は、安定余裕と追従特性に関する調整目標を厳密な目標に、残りを柔軟な目標に指定して解いています。

それでは、得られた解について解析してみたいと思います。

ノミナルチューニングの解析

% 目標値から出力までの閉ループ伝達関数をモデルから抽出
Tini = getIOTransfer(ST0,'r','y');              % チューニング前
Tnom = getIOTransfer(ST,'r','y');               % チューニング後

% 外乱から出力までの閉ループ伝達関数をモデルから抽出
Tdini = getIOTransfer(ST0,'d','y');             % チューニング前
Tdnom = getIOTransfer(ST,'d','y');              % チューニング後

% 目標値に対するステップ応答の比較
clf,step(Tini,'m--',Tnom,'b-'),set(gca,'Ylim',[-1 5]),legend('Initial gain','Tuned gain')

% 外乱に対するステップ応答の比較
clf,step(Tdini,'m--',Tdnom,'b-'),set(gca,'Ylim',[-5 5]),legend('Initial gain','Tuned gain')

% 調整目標に対する性能評価
viewSpec([Req1 Req2 Req3 Req4],ST0)             % チューニング前
viewSpec([Req1 Req2 Req3 Req4],ST)              % チューニング後

応答解析に際して、getIOTransferコマンドを使用してモデルから指定したポイント間の伝達特性を抽出します。
抽出した伝達関数モデルはstepコマンド等で時間応答の解析が行えます。

上図は左が目標値、右が外乱に対するステップ応答を示しています。
図中、ピンクの破線がLQRの値、青の実線がsystuneによるチューニング後の結果を示しています。
結果を見ると分かる通り、LQRは制御が発散しています。
一方、systuneによる調整後の結果では応答が安定していることがわかります。

続いて、調整目標に対する結果も見ておきたいと思います。
まずは、LQRの場合です。

左上から右下に向かって、安定余裕、目標値に対する過渡応答特性、外乱抑圧性、目標値追従性を示しています。
時間領域のものもあれば周波数領域で示されているものもあります。
結果を見ると、どれも満足な結果とはなっていないことが分かります。

続いて、systuneの場合です。

安定余裕、目標値追従性の結果を見てみると、クリーム色でエッジされている領域に掛らなければ目標を十分に満たすことを示しますが、残念ながら一部が掛かってしまっていることが分かります。
しかし、LQRと比較すると十分な性能であると言えます。特に安定余裕が全周波数帯に亘って確保されました。
過渡特性、外乱抑圧性についてもまずまずの性能であると言えそうです。

それでは、この制御パラメーターを使用して、制御対象側のパラメーター変動を考慮した際にロバストであるかを確認してみます。

ロバスト性の評価

% Simulinikモデル上の該当するブロックに不確かさを含む線形化の定義を指定
blocksubs(1).Name = [mdl '/Longitudinal Plant/A'];
blocksubs(1).Value = A;
blocksubs(2).Name = [mdl '/Longitudinal Plant/B'];
blocksubs(2).Value = B;
blocksubs(3).Name = [mdl '/Actuator/Tde'];
blocksubs(3).Value = 1/Tde;
blocksubs(4).Name = [mdl '/Actuator/Delay'];
blocksubs(4).Value = c2d(usample(Pdt,5),ts,'tustin'); % むだ時間は1次のパデ近似モデルとして、ランダムに5サンプル抽出したモデルで置き換える
blocksubs(5).Name = [mdl '/dist'];
blocksubs(5).Value = d;

% 再度slTunerで情報を取得
UST0 = slTuner(mdl,{'Ki','Kx'},blocksubs,options);

% ノミナルチューニングで調整したゲイン値をUST0に反映し、応答を確認
setBlockValue(UST0,getBlockValue(ST));
Tnom = getIOTransfer(UST0,'r','y');
Tdnom = getIOTransfer(UST0,'d','y');
clf, step(Tnom,30), grid

変動を含んだ線形システムのA、B行列、アクチュエーターの時定数、むだ時間（1次のパデ近似）、外乱の振幅ゲインをそれぞれモデル上の各ブロックにslTunerを介して反映していきます。

変動を反映後に再度slTunerでUST0という名でオブジェクトを作成し、ノミナルチューニングで調整したゲイン結果をgetBlockValueコマンド、setBlockValueコマンドを使用してUST0に反映します。

getIOTransferコマンドを使用して、目標値、外乱に対する出力までの伝達関数を取得して、各ステップ応答を見てみましょう。

図の並びは先ほどと同じです。
ランダムにサンプリングされた代表的な変動ケースに対する応答が示されていますが、いずれも振動的かつ発散してしまうケースがあることが分かります。

以上からノミナルケースでチューニングした制御ゲインではロバスト性に問題があると結論付けられそうです。

いよいよロバストチューニング

ということでノミナルチューニングではロバスト性に難があるので変動を考慮したシステムに対してsystuneを掛けてみたいと思います。
すでに変動を考慮したslTunerオブジェクトUST0を作成済みであるため、あとはこれをsystuneに掛けるだけになります。

ロバストチューニング

Req1.Focus = [1e-3 10];    % 安定余裕の調整目標について、フォーカスする帯域を0.001~10[rad/s]に限定
soft = [Req2 Req3];
hard = [Req1 Req4];

opt = systuneOptions('Display','off','UseParallel',true);  % 並列化を有効
[UST, fsoft, ghard] = systune(UST0,soft,hard,opt);

なお、ロバストチューニングでは、調整目標が厳しすぎて解なしになることを防ぐため、安定余裕の目標については指定する周波数帯を限定（目標を緩和）しておくことにします。

またsystuneではParallel Computing Toolboxの機能を使って並列化のオプションを使用できるため、計算の高速化を狙って並列処理を有効にしておきます。

ロバストチューニングの解析

UST.Ts = ts;
Trob = getIOTransfer(UST,'r','y');
Tdrob = getIOTransfer(UST,'d','y');

clf, step(Tnom,Trob,100), grid
set(gca,'Ylim',[-2 2])
legend('Nominal tuning','Robust tuning')

clf, step(Tdnom,Tdrob,100), grid
legend('Nominal tuning','Robust tuning')

viewSpec([Req1 Req2 Req3 Req4],UST0)
viewSpec([Req1 Req2 Req3 Req4],UST)

以上の設定の下で、チューニングした結果が以下になります。
まずはステップ応答の確認です。

図中、赤の実線がロバストチューニング、青の実線が比較対象として載せたノミナルチューニングの結果です。
ランダムにサンプリングされた代表的な変動ケースに対し、ロバストチューニングの結果は良好であることが分かります。

続いて、調整目標に対する性能評価です。

上図はノミナルチューニングの結果ですが、安定余裕において2、6[rad/s]あたりで著しく減少していることが分かります。
また目標値に対する過渡応答特性及び外乱抑圧性の結果を見ると、振動的かつ不安定なケースがあることが分かります。

続いて、上図がロバストチューニングの結果となります。
安定余裕に関して目標を完全には達成できないものの、ノミナルチューニングの結果と比較して改善していることが分かります。
特に2[rad/s]における余裕が大きく改善しています。
（6[rad/s]におけるゲイン余裕が依然としてかなり厳しいですが、これは制御構造そのものを検討しなければならないかもしれません）

また目標値に対する過渡応答特性、外乱抑圧性もノミナルチューニングと比較して改善しています。

ただし、追従特性については低周波帯域で顕著に劣化していることが確認されます。
実際に先ほどのステップ目標に対する応答結果を見ても明らかですが、偏差のバラつきがノミナルチューニングの結果（発散ケースを除き）と比較して大きいです。
このことから複数の調整目標を一つの構造固定化されたコントローラーで完全に満足するということは難しく、一般にトレードオフが発生することが分かります。

以上から変動を含んだシステムに対し、ロバストチューニングの結果は概ね一定の改善効果を見込めるということが結論づけられました。

最後に元のSimulinkモデルにおいて、ロバストチューニングの結果を反映し、かつモンテカルロシミュレーションを実施してロバスト性の確認をしておきたいと思います。

モンテカルロシミュレーション

% parsminを用いたモンテカルロシミュレーション
rng('default')                  % 乱数のシード設定
sample = 100;                   % MCSの実行回数

% SimulationInputオブジェクトを利用して各実行シナリオを定義していく
for i = sample:-1:1
    in(i) = Simulink.SimulationInput(mdl);
    % setBlockParameterを利用してブロックパラメーターを生成したサンプルの値で置換していく
    in(i) = in(i).setBlockParameter([mdl '/Longitudinal Plant/A'],'Gain','usample(A,1)');
    in(i) = in(i).setBlockParameter([mdl '/Longitudinal Plant/B'],'Gain','usample(B,1)');
    in(i) = in(i).setBlockParameter([mdl '/Actuator/Tde'],'Gain','1/usample(Tde,1)');
    in(i) = in(i).setBlockParameter([mdl '/Actuator/Delay'],'DelayTime','usample(DT,1)');
end

% parsimで実行開始
out = parsim(in,'ShowProgress','on','TransferBaseWorkspaceVariables','on','ManageDependencies','on');

% 結果の描画
figure,clf
subplot(211)
for i = 1:sample
    plot(out(i).logsout.get('r').Values.Time,out(i).logsout.get('r').Values.Data*180/pi,'m--'),hold on
    plot(out(i).logsout.get('y').Values.Time,out(i).logsout.get('y').Values.Data*180/pi,'b-'),hold on
end
hold off
grid on
xlabel('Time [m]')
ylabel('Pitch angle [deg]')

subplot(212)
for i = 1:sample
    plot(out(i).logsout.get('de').Values.Time,out(i).logsout.get('de').Values.Data*180/pi,'m--'),hold on
end
hold off
grid on
xlabel('Time [m]')
ylabel('de [deg]')

上図はピッチ角の目標値として5[deg]のステップ信号を入れた際の応答結果（100回分）であり、上段がピッチ角、下段が制御入力であるエレベーター舵角のトレンドを示しています。
若干一部のケースで振動的な応答が見られますが、概ね安定した制御性能を示していることが確認できます。
（systune調整時はむだ時間を1次のパデ近似で線形化して評価していますが、元のSimulinkモデルは非線形ですのでその影響が出現しているかもしれません）

以下が代表ケースに対するUnreal Engineを用いたアニメーションによる可視化結果（3倍速）です。
特におかしな挙動はなさそうです。

最後に

今回は構造固定化された制御器に対する制御パラメーターのロバストなチューニングを試してみました。
制御系設計を行う際に、新たに制御系の構造そのものを設計できるケースというのは稀であり、通常は実績/信頼あるものを踏襲した形でものづくりが進められるケースの方が多いと思います。
そのような場合において、今回のようなアプローチが検討できるかもしれません。