一次遅れ＋むだ時間系のモデルと数値シミュレーション（外乱応答編）
前回の続きで、「一次遅れ＋むだ時間」系のモデルに外乱が加わった際のステップ応答シミュレーションについて書きたいと思います。外乱はプロセス変数（温度、圧力、流量 etc.）を乱す要因です。プロセス制御では外乱を相手にどれだけ一定の値で制御（定値制御）できるかが重要だったりするので、外乱を考慮した制御系設計は身に付けてしかるべきです。一方で、外乱を考慮したシミュレーションは、実装の形（モデルの表し方）や解釈の仕方に一癖も二癖もあるので、今回は一番直観的なものを考えていきます。
キーワード：モデル、差分方程式、一次遅れ＋むだ時間、外乱応答、プログラミング
１．外乱の生成
外乱といってもその時間応答波形は様々です。ここでは、白色ノイズ（平均：0、分散：\(s^2\)）と成形フィルタによって生成される有色ノイズを考えてみます。適当な外乱データがあればよかったのですが、白色と有色の外乱の面白い関係を見るために、あえて人工的な信号を作ります。
１．１．白色ノイズ
MATLAB／Octaveで、以下のように\({\tt randn()}\)関数を使うと任意の白色ノイズ、いわゆる雑音が作れます。図１の上真ん中が\({\tt xi }\)です。
N = 260; % データサイズ m_xi = 0.0; % 平均 v_xi = 1.0; % 分散 xi = m_xi + sqrt(v_xi)*randn(N);
１．２．有色ノイズ
有色ノイズは、無作為に変動する白色ノイズに比べ何かしらの傾向が見える外乱です。ここで、白色ノイズは周波数に対してパワースペクトルが平坦な信号のことですが、言い換えると、あらゆる周波数をもったsin波の重ね合わせにより信号が無造作に（見えるように）変動するとしたら、パワースペクトルを偏らせるとどうなるでしょうか。今回は、これをローパスフィルタ（成形フィルタ）で検証してみます。白色ノイズを低周波帯域から高周波帯域に向かってカットした信号を生成してみます。以下図１で、下左は400、下真ん中は40、下右は4のフィルタ時定数（カットオフ周波数の逆数）からなる一次遅れフィルタを先ほどの白色ノイズに適用した信号です。
図１：白色ノイズ（\(m_{\xi}=0.0\)、\( s_{\xi}^2=1.0 \)）と成形フィルタによって生成された有色ノイズ
若干の作為的な操作をしましたが、理想的には、パワースペクトルが低周波帯域に偏るほど傾向をもった波形（有色ノイズ）に近づきます。直感的に言えば、左のローパス特性が強いフィルタを適用すれば、遅い周波数のsin波が残り、傾向が支配的な波形に、逆に、右のローパス特性が弱いフィルタを適用すれば、速い周波数のsin波も依然として残っているので白色ノイズに近いままの波形になります。以上が、白色と有色の外乱の面白い関係でした。傾向をもった波形は偏ったパワースペクトルをもっているのです。そして、有色の外乱がどうやって発生するのかが”遅れをもったシステム（ローバス特性）”の観点から見えてきそうな気がします。
今回、ローパスフィルタの適用に際して実装方法（ソースコード）は省略しましたが、「フィルタ」をお題にして別途記事を書きたいと思います。本題として、外乱を含んだシミュレーションは今回生成したものを使って次の章で解説します。
２．外乱応答シミュレーション
紹介する外乱応答シミュレーションは直感的なものだと先述しました。ブロック線図と伝達関数で表すと、図２や（\ref{eq:eq8-1}）式のような出力端に外乱が加わる状況を想定します。
図2：出力端に外乱があるブロック線図
\begin{align} Y(s) &= \frac{K}{1+Ts}e^{-Ls}U(s) + D(s) \label{eq:eq8-1} \end{align}
出力\(Y(s)\)に外乱がそのまま重畳する形であり、（\ref{eq:eq8-1}）式からも微分方程式の解に外乱項が加算されたわりかしイメージ通り（？）な形となっています。それでは、（\ref{eq:eq8-1}）式を離散化（後進差分近似）した上で整理した差分方程式を以下に示します。登場する係数も計算の過程も前回の記事と全く同じなので、最終形だけにします。
\begin{align} y(k) = -a_1y(k - 1) + b_0u(k - t_d - 1) + d(k) + a_1d(k - 1) \label{eq:eq8-2}\end{align}
ただし、\(d(k)\)は離散時間領域における外乱項です。（\ref{eq:eq8-2}）式を実装したソースコードは末尾に載せました。また、１章で生成した４つの外乱を加えた時間応答シミュレーションは図３に示します。上から白色ノイズ、\(T_c=400, 40, 4\)の成形フィルタで生成した各有色外乱です。なお、信号の強度がステップ応答に比べ弱かったことから、それぞれの外乱を少しばかし増幅して実行しました。
図3：外乱を考慮した時間応答シミュレーション
各外乱の特徴がそのままステップ応答波形に現れており、制御量を乱しています。温度がどんどん上昇、ハンチング気味、ノイズのような微小振動．．．
ちなみに、（\ref{eq:eq8-2}）式の離散時間領域における伝達関数（パルス伝達関数）も以下のような出力端に外乱項が加わる形できちんと表せます。
\begin{align} y(k) = \frac{z^{-(t_d + 1)}b_0}{1 + a_1z^{-1}}u(k) + d(k) \label{eq:eq8-3} \end{align}
この形のモデルはOEM(Output Error Model)と呼ばれるそうです。主にシミュレーション向きのモデルとして使われます。
３．まとめ
今回、「出力端」に関する外乱の時間応答シミュレーションの実装方法を示しました。これで自由自在に本物のデータさえ引っ張ってくれば、ほぼほぼリアルな？プロセスのシミュレーションができます。しかしながら、外乱自体の話はこれだけでは終われません。外乱は様々なところから作用してきます。操作量に外乱が印可し、プロセスを介して出力に現れる「入力端」外乱もあります。そもそも、外乱がどの段階でどの変数に影響を与えるのかはよくよくプロセスや設備特性を考えてからモデルを構築・実装しなければなりません。方や、制御系設計で設計者がどちらの外乱で考える方が都合が良いみたいな設計寄りの目的で選ばれるときもあります。ちなみに、今回の直感的なシミュレーションモデル（OEM）は、よくセンサノイズが外乱として系に加わる例えとしてよく見られます。
番外編として、記事を書いていて気づいたのが、まさかのLaTeXのσが使えない！自動リンク機能の文字列に\sigmaが当てはまるみたいです。実は他にも-1のマイナスと1をくっつけて書くと自動リンクでMathJaxがエラー．．．
FILENAME = 'D.txt';
DV = load(FILENAME);
% 設計パラメータ
Ts = 1; % サンプリング時間
N = 260; % サンプル時刻の終点（シミュレーション回数）
K = 2; % システムゲイン
T = 40; % 時定数
L = 10; % むだ時間
a1 = -T/(Ts+T);
b0 = K*Ts/(Ts+T);
td = floor(L/Ts);
offset_time = max([1, 0, td+1]);
y = zeros(N + offset_time, 1);
u = [zeros(offset_time, 1); ones(N, 1) * 5];
d = [zeros(offset_time, 1); D(:, 2)];
for k = (1 : N) + offset_time
y(k) = - a1*y(k - 1) + b0*u(k - td - 1) + d(k) + a1*d(k - 1);
end
y = y((1 : N) + offset_time, 1);
u = u((1 : N) + offset_time, 1);
d = d((1 : N) + offset_time, 1);
ソースコードを前回のものから若干変えました。MATLABで因果関係を実装するのはやりづらいです。センスがないので、もっと効率的なコードがあるかもしれません。シミュレーションの開始・現在時刻を\(k=1\)から始めているので、MATLABで\({\tt y(k - 1)}\)などのインデックスが0や負になるような過去データを参照するとエラーが出てしまうので、対策として、\({\tt for}\)文中で参照する際には\({\tt offset_time}\)でずらしています。\({\tt offset_time}\)は基本的に\({\tt max}\)(伝達関数の分母の次数, 伝達関数の分子の次数, むだ時間+1)で選ぶか、適当な長さを指定します。
