一次遅れ＋むだ時間系のモデルと数値シミュレーション（ステップ応答編）
前回の記事では、基本的なPID制御則のアルゴリズムについて書きました。今回は、それを使ってPID制御のシミュレーションをするためにもう一つ必要な制御対象の「モデル」の考え方・実装方法について書きます。これで晴れて制御系設計のためのシミュレーションができるようになります。
※言葉足らずな文章でしたので、若干加筆しました＠2017.12.05
キーワード：モデル、差分方程式、一次遅れ＋むだ時間、ステップ応答、プログラミング
はじめに
シミュレーションを使った設計は、なにも最新の研究や自動車開発の十八番ではありません。製造業においてもその恩恵はいかんなく発揮されるはずです。それでは、必要となるプラント、あるいは制御対象の微分方程式／数式モデルをどのように準備するかですが、実験室スケールの化学反応や明快な運動方程式をいつも考えるわけではありません。ことプロセス制御にいては。製造業における制御対象は非線形・システム変動・高次遅れなどが幾重にも重なり複雑極まりない大きなシステムです。一つひとつ方程式を立てていく訳にはいかなく、ある程度理にかなったブラックボックスモデルを使います。
一次遅れ＋むだ時間系のモデルによる近似
そこで、製造業では制御対象を（\ref{eq:eq7_1}）式の「一次遅れ＋むだ時間」系でばっさり近似したモデルをよく使います。
\begin{align} Y(s) &= \frac{K}{1+Ts}e^{-Ls}U(s) + D(s) \label{eq:eq7_1} \end{align}
\(K\)はシステムゲイン、\(T\)は時定数、\(L\)はむだ時間を表したシステムパラメータで、\(Y(s)\)、\(U(s)\)、\(D(s)\)はそれぞれ制御量、操作量、外乱です（ラプラス変換）。詳細な原理原則は目をつむり、実用性を重視したモデリングで、一次遅れ＋むだ時間系は
- ステップ応答波形からシステムパラメータが目視で確認できる。
- 化学工学で重要な収支モデル（一階の微分方程式）の伝達関数である。
- プラントの時間応答によく現れる波形である。
- システムパラメータの同定がし易く、動特性を簡単かつ最小限に表したモデルである。
- CHR法やZN法からPIDパラメータを求める公式が開発されており、現場で使いやすい。
といった特徴があります。
後進差分近似によるシミュレーション
続いて、（\ref{eq:eq7_1}）式を実装する方法について説明し、例題を示したいと思います。MATLABのSimulinkやtf関数を使えば一発ですが、中身は謎！この方法がわかれば、別にMATLABやC言語じゃなくてもExcelなんかで簡単に制御対象の時間応答シミュレーションができます。
本題として、実際に実装するのは（\ref{eq:eq7_1}）式の伝達関数ではなく、その微分方程式になります。以下がその微分方程式です。
\begin{align} y(t) + T\frac{dy(t)}{dt} = Ku(t-L) \label{eq:eq7_2} \end{align}
（\ref{eq:eq7_2}）式では外乱項は記事のフットワークを軽くするために省略しました。詳細な外乱項を導入したモデル（CARMA／CARIMAモデル）は、別途記事にしたいと思います。それでは、プログラミングで実装するために（\ref{eq:eq7_2}）式を後進差分近似（離散化）します。
\begin{align} y(k) + T\frac{y(k) - y(k - 1)}{T_s} &= Ku(k - t_d - 1)\notag \\ y(k) &= - a_1y(k - 1) + b_0u(k - t_d - 1) \label{eq:eq7_3}\end{align}
\begin{align} a_1=- \frac{T}{T_s + T}, b_0 = \frac{KT_s}{T_s + T}, t_d={\tt floor}(\frac{L}{T_s}) \notag \end{align}
要するに、微分の定義に基づく近似微分です。なお、\(T_s\)は微小な時間、サンプリング時間、\(k\)は整数のサンプル時刻です。ここで、今回はゼロ次ホールドまでを考慮した離散化を行わなかったので、操作量\(u(k - t_d {\bf - 1})\)の項のサンプル時刻に勝手に\(- 1\)を追加しています。ゼロ次ホールドについても別途記事にしたいと思いますので、今回の一先ずの理屈としては、操作量の変化は最低でも\(1\)サンプル時刻遅れて制御量に現れる／直達項はない、として下さい。最後に、連続時間領域（アナログ）における現象の振舞いを微分方程式として表現しますが、コンピュータのように値が飛び飛びの離散時間領域（ディジタル）では（\ref{eq:eq7_3}）式のような差分方程式で表されます。
以上でモデルの考え方と実装方法についての説明を終わります。他にもあれやこれや書きたいことはありますが、一先ず、（\ref{eq:eq7_3}）式をガリガリ計算すれば制御対象の一次遅れ＋むだ時間系の時間応答シミュレーションができます。以下では、外乱なしの例題を示します。システムパラメータの詳細は図の下に。記事の末尾にはソースコードも記載します。
図１：一次遅れ＋むだ時間系のステップ応答シミュレーション（\(K=2\)、\(T=40\)、\(L=10\)、\(T_s=1\)）
ソースコードは、操作量\(u(k)\)にステップ状の信号を与えた際に、制御量\(y(k)\)のステップ応答波形をシミュレートするものになっています。たとえば、燃料を供給するフィーダをスッテプ状に操作したとき、温度が時間とともに変化する状況を想定しています。\(u(k)\)に与える信号の形やシステムパラメータを変えて、いろいろな時間応答シミュレーションを試してみてください。・・・2017.12.05追加
↓↓外乱編の記事です↓↓
まとめ
今回の一次遅れ＋むだ時間系のシミュレーションは、実装するアルゴリズムも比較的簡単なので、製造の現場でもExcel一つで応答波形をシミュレーションでき、PID制御のアルゴリズムも合わせれば、制御系設計のシミュレーション（次回）が実際にできます。外乱がないので、「こんなきれいな波形で温度は変化しない！」にも対応できるよう外乱を考慮した記事も追って書きたいと思います。（記事を書きました）シミュレーションに基づいた開発・改善・解析は、実機だけのトライアル＆エラーによる負担やコストを軽減し、安定性の確認にもつながります。最後に、制御対象がどうしてあのモデルになるのかという理屈やプロセスの仕組み・設備仕様をすっ飛ばしてブラックボックスモデリングをやりましたが、システムパラメータの値はどこから来たんだ？という大きな疑問が残っています。モデルの選定やシステムパラメータの推定をひっくるめてシステム同定といいますが、これはこれで非常に大切な問題で、これについては日を改めて書きたいと思います。
ところで、学部生の時に受けた制御工学の講義では、\(s\)をどう実装して時間応答のシミュレーションをするのかお察しな疑問を抱いていました（笑）。その講義ではシミュレーション・実験まではやっていなく、まともなディジタル制御の講義がでてきたのは院生からでした（汗）。結局のところ、情報工学系の講義で数値シミュレーション、信号処理工学で離散化の勉強をしてやっと話がリンクし、伝達関数をそのまま実装するんじゃないことに気づかされました。
% 設計パラメータ
Ts = 1; % サンプリング時間
N = 260; % サンプル時刻の終点（シミュレーション回数）
K = 2; % システムゲイン
T = 40; % 時定数
L = 10; % むだ時間
a1 = -T/(Ts+T);
b0 = K*Ts(Ts+T);
td = floor(L/Ts);
y = zeros(N, 1);
u = zeros(N, 1);
d = zeros(N, 1);
u(:) = 1;
for k = 1 + max([1, 0, td+1]) : N
y(k) = - a1*y(k - 1) + b0*u(k - td - 1);
end
