はじめに
実応用上優れた特性を有する方法として，Virtual reference feedback tuning(VRFT)と呼ばれる方法が提案されています。VRFTは一度実験するだけで所望の制御器を自動獲得可能なオートチューニング手法です。VRFTは様々な利点が存在しますが，国内ではあまり活用されていません。そういった経緯から，前回記事ではVRFT法を取り上げ、VRFTのどういった部分が優れているのか・制御器調整の原理・運用上での注意点等をいくつか解説しました。
前回は読者の皆様がVRFTを体感いただけるようにPythonを使ったツールのサンプルコードを載せて、簡単な数値例を通してPIDゲインを一度の入出力データのみで自動調整できることを解説しました。ありがたいことに記事自体は大変好評でした。ただ、一部の方から「ぜひMATLAB版のプログラムも配布してほしい！」というご意見をいただきました。
そういった経緯から，本記事ではMATLABユーザー様にもVRFTを体感いただけるようにMATLABを使ったVRFTのサンプルコードを解説していきたいと思います。なお，本記事のサンプルコードではMATLAB(R2021a)本体以外にControl system ToolBoxも必要なのでご注意ください。
VRFTの概要
まず，VRFTの概要を述べておきます。前回記事を既にご覧の方は飛ばしていただいても問題ありません。Fig.1のフィードバック制御系を考えます。オートチューニングの前提として，規範モデル$T_{d}$をあらかじめ与えます。$T_{d}$は閉ループ系の要求仕様を伝達関数として定量化したものだと考えてください。VRFTではFig.1の閉ループ伝達関数を規範モデル$T_{d}$に一致させる制御器$C$を導出することを目指します。
Fig.1 SISOの制御系
まず，開ループ/閉ループ実験のどちらかで制御対象の初期入出力データ(時系列データ)を取得します。次に，初期入出力データを用いて擬似誤差信号$e_{s}$と整形後制御入力$u_{s}$を
u_{s}=L u_{0}
e_{s}=L\left(\frac{1}{T_{d}} y_{0}-y_{0}\right)
と導出します。$u_{0}$と$y_{0}$は初期実験で取得した入出力応答の時系列データです。$L$はノイズ除去やノンプロパーな伝達関数計算を回避するためのプレフィルタです。上記の信号を用いて，VRFTの評価関数を
J_{VRFT}(C)=\left\|u_{s}-C e_{s}\right\|^{2}
と定義します。この評価関数は制御器$C$でパラメタライズされています(厳密な表記ではありませんが，わかりやすさ重視でこのような表記をしています)。何らかの数理最適化手法にて評価関数$J_{VRFT}$を最小化するような制御器$C$のパラメータを探索します。例えば，制御器$C$が評価関数に対して線形関数になれば最小二乗法にて評価関数を最小化できます。評価関数を十分小さくできれば，規範モデルと閉ループ伝達関数を一致させるような制御器$C$を導出できます。VRFTによる制御器更新手順は下記のようにまとめることができます。
1.プレフィルタ$L$と規範モデル$T_{d}$をあらかじめ与える。
2.開ループ/閉ループ実験のどちらかで入出力応答を測定する。
3.$J_{VRFT}$を最小化する制御器$C$のパラメータを数理最適化にて導出する。
VRFTは一組の実験データを測定するだけで未知システムに対する制御器をオートチューニングできます。制御器$C$がPID制御器であれば，最小二乗法にてPIDゲインを導出できます。つまり，最適化における計算コストが非常に低く，パラメータの最適性が保証されます。PIDゲインを最小二乗法で求解できる点は，VRFTの強い利点と言ってもいいでしょう。本記事はMATLABのサンプルコードを中心とした解説なので、VRFT自体の解説は概要レベルに留めることにいたします。詳しい理論や利点を知りたい方がいらっしゃれば，前回記事で解説しておりますので，そちらをご確認いただければと思います。
参考文献
Campi, M. C., Lecchini, A., & Savaresi, S. M. (2002). Virtual reference feedback tuning: a direct method for the design of feedback controllers. Automatica, 38(8), 1337-1346.
MATLABによるVRFTを用いたPID自動調整
それでは，本命であるMATLABによるPID制御器調整のデモをしてみましょう。実を言うと，VRFTのMATLABプログラム自体はVRFTの生みの親であるCampi先生のホームページで公開されています。
VRFT_toolboxというプログラムとそれの解説用国際会議論文が紹介されており，大変分かりやすく完成度の高い内容になっています。ただ、System identification Toolboxも必要ということもあって，人によっては少し敷居が高いかもしれません。そういった理由から，本記事ではCampi先生のプログラムはあえて利用せず，myVRFTという関数をMATLABで自作し，それを利用してPIDゲインを自動調整するデモを実施します。
例題
前回と同じく，簡単な例として一次遅れ系の制御対象に対するPID制御器調整問題を考えてみましょう。制御対象の伝達関数を
P=\frac{1}{s+1}
と設定します。また，制御器を一般的なPID制御器：
C=K_{P}+\frac{K_{I}}{s}+\frac{K_{D}s}{0.05s+1}
としました。本例題では，PIDゲインの初期値は$K_{P}=0.1$,$K_{I}=0.5$,$K_{D}=0.01$と設定しました。最後に，規範モデルを
T_{d}=\frac{1}{0.5s+1}
と与えました。以上の問題設定のもと，MATLABを使ってPIDゲインを調整することを目指しましょう。PIDゲインを自動調整するシミュレーションプログラムとVRFTの関数myVRFTを下記のように作成してみました。
%test_VRFT.m
close all
clear
%テストデータ生成処理(初期制御器による初期シミュレーション
kp0=0.1;%比例ゲイン(初期値)
ki0=0.5;%積分ゲイン(初期値)
kd0=0.01;%微分ゲイン(初期値)
P=tf([1],[1,1]);%制御対象
Cfb=kp0*tf([1],[1])+ki0*tf([1],[1,0])+kd0*tf([1,0],[0.05,1]);%初期制御器
t=(0:0.01:5)';%シミュレーション時間(サンプリング0.01s)
ref=ones(length(t),1);%目標値計算処理
ref(1)=0;%目標値計算処理(ステップ指令)
y=lsim(P*Cfb/(1+P*Cfb),ref,t);%閉ループ系の出力信号シミュレーション
u=lsim(Cfb/(1+P*Cfb),ref,t);%閉ループ系の制御入力シミュレーション
Ts=0.5;%規範モデルの時定数
Td=tf([1],[Ts,1]);%規範モデル
[kp,ki,kd]=myVRFT(y,u,t,Td);
Cfb=kp*tf([1],[1])+ki*tf([1],[1,0])+kd*tf([1,0],[0.05,1]);%制御器更新
ya=lsim(P*Cfb/(1+P*Cfb),ref,t);%閉ループ系の出力信号シミュレーション(調整後)
ua=lsim(Cfb/(1+P*Cfb),ref,t);%閉ループ系の制御入力シミュレーション(調整後)
% 初期シミュレーション結果のプロット
figure(1)
plot(t,y,'linewidth',8);
hold on;
plot(t,ref,'-.','linewidth',8);
hold on;
plot(t,lsim(Td,ref,t),':','linewidth',8);
hold on;
xlabel('Time [s]');
ylabel('Output');
legend('Output','reference','desired');
hold on;
set(gca,'FontSize',45);
xlim([0 5]);
% 調整後シミュレーション結果のプロット
figure(2)
plot(t,ya,'linewidth',8);
hold on;
plot(t,ref,'-.','linewidth',8);
hold on;
plot(t,lsim(Td,ref,t),':','linewidth',8);
hold on;
xlabel('Time [s]');
ylabel('Output');
legend('Output','reference','desired');
hold on;
set(gca,'FontSize',45);
xlim([0 5]);
function [kp,ki,kd]=myVRFT(y,u,t,Td)
L=Td;%プレフィルタ(カスタマイズした場合はここを変更)
pn=3;
A=ones(pn,length(y));
intg=tf([1],[1,0]); %積分
dif=tf([1,0],[0.05,1]);% 微分
A(1,:)=lsim(L*(1-Td)/Td,y,t);
A(2,:)=lsim(L*intg*(1-Td)/Td,y,t);
A(3,:)=lsim(L*dif*(1-Td)/Td,y,t);
A=A';
ul=lsim(L,u,t);
wp=pinv(A)*ul;
kp=wp(1);%調整後比例ゲイン
ki=wp(2);%調整後積分ゲイン
kd=wp(3);%調整後微分
このプログラムを実行すると，初期制御器による閉ループシミュレーションを実施し，その結果をmyVRFT関数に入力してPIDゲインを計算してくれます。最後に，初期シミュレーションの結果とVRFTで調整した後のシミュレーション結果をプロットしてくれます。実行してみると，初期のPIDゲインを使った閉ループシミュレーションの閉ループ応答がFig.2のようにプロットされます。出力が規範モデルの目標値応答と一致していないことがわかります。したがって，出力と規範モデルの応答が一致するようにPIDゲインを調整しないといけません。
Fig.2 初期データ(青：出力、黄色：規範モデルの目標値応答、赤：目標値)
この閉ループ応答データを用いて，myVRFTにてPIDゲインを調整してみます。その結果，PIDゲインの値をそれぞれ，$K_{P}= 2.00$, $K_{I}= 2.00$, $K_{D}= 0.0$と得ることができました(実際には計算されたゲインに対して，有効数字で切り捨て・切り上げしています)。制御器調整後の閉ループシミュレーションの結果が下記のFig.3になります。出力が規範モデルの応答と重なっており，望ましい制御器調整を実現できています。本結果から，MATLABを利用することで一度入出力データを測定するだけでPID制御器をオートチューニングできることがわかります。なお，プレフィルタに関しては前回記事で解説したように$L=T_{d}$と設定しています。プレフィルタ設定をカスタマイズすれば，さらに高性能なPIDゲイン調整もできるので，興味のある方は試してみてください。
Fig.3 調整後データ(青：出力、黄色：規範モデルの目標値応答、赤：目標値)
おわりに
本記事では，VRFTと呼ばれる一度実験するだけで所望の制御器を自動獲得可能なオートチューニング手法をMATLABで実装してみました。myVRFTの関数を利用すれば，実際の実験データを使ってPIDゲインをオートチューニングすることもできます。MATLABユーザーの方はぜひご自身の業務でも活用してみてください。