ディジタルPID制御の時間応答シミュレーション - とりあえずフィードバック制御

　ディジタルPID制御のシミュレーションをやっていこうかと思います。ソースコードを実装するための言語はMATLAB／Octaveになりますが、可能な限り言語に依存しない書き方で作ってみました。

キーワード：MATLAB、Octave、シミュレーション、PID制御、ディジタル制御

はじめに

　シミュレーションに必要な前知識は以下の通りです。

■ 制御対象の離散化

■ PIDコントローラの離散化

PID制御の実装－位置型と速度型－ - とりあえずフィードバック制御

■ ソースコードと解説

　この記事以外でも今回のソースコードをベースに書いていきます。

＜MATLAB／Octaveコード＞

% ------------------------------------------------
% 
%                   PID Simulator
% 
% ・プラントは一次遅れ＋むだ時間系のモデル
% ・PIDコントローラはPID制御則
% ------------------------------------------------

clear;
close   all;
format  compact;

% ------------- 初期化 --------------

% プラントモデルと外乱
K  = 2.0;               % @システムゲイン
T  = 40.0;              % @時定数
L  = 1.0;               % @むだ時間
DV = load('DV.txt');    % @外乱

% PIDコントローラ
kp = 1.5;               % @比例ゲイン
Ti = 50.0;              % @積分時間
Td = 0.0;               % @微分時間

% シミュレーションに必要な設定
Ts = 1.0;               % @サンプリング時間（制御周期）
N  = length(DV(:, 1));  % シミュレーションのデータサイズ

% 離散化したモデルのパラメータ
a1 = -T/(Ts+T);
b0 = K*Ts/(Ts+T);
td = floor(L/Ts);

% プロセス変数
k_offset = max([1, 0+td+1, 2]);          % 因果性のためのオフセット[分母の次数、分子の次数、コントローラの次数]
y = zeros(N + k_offset, 1);              % 制御量
r = [zeros(k_offset, 1); ones(N, 1)*5];  % @目標値
u = zeros(N + k_offset, 1);              % 操作量
d = [zeros(k_offset, 1); DV(:, 1)] * 0;  % @外乱



% ---------- PID制御シミュレーション -----------

for k = (1:N)+k_offset
    
    % ①プラントシミュレーションと制御量の計測
    y(k) = -a1*y(k-1)+b0*u(k-td-1)+d(k)+a1*d(k-1);
    
    % ②制御偏差の計算
    e(k) = r(k)-y(k);
    
    % ③PID制御コントローラによる操作量の計算（PID制御則）
    du   = kp*((e(k)-e(k-1))+Ts/Ti*e(k)+Td/Ts*(e(k)-2*e(k-1)+e(k-2)));
    u(k) = u(k-1)+du;
    
end

自由に設定ができる部分は、コメントに「@」付きで示したものです。
「@外乱」は、末尾の0 or 1で外乱の印可を切替。「DV(:, 1)」の列要素は1~4までの中から一つを指定して、前回の記事で生成した4つの外乱を選択できます。外乱のサンプルデータは、こちら（まだ適当なアップローダを用意していないので．．．）。
「@目標値」も*5を書き換えることで変更可能。
「k_offset」は、時刻が0以前のような過去データを変数が表すとき、そのままMATLAB／Octaveのベクトルで感覚的にy(- 1)みたいに指定するとエラーが出るので、参照する過去データの最長の値で指定する要素に下駄をはかせています。
「for」文は、時間（正しくはサンプル時刻）の経過を意味し、中身の処理は、①制御量の計算（計測）⇒②制御偏差の計算⇒③操作量の計算（出力）がワンセットとなっており、一周で1サンプル時刻が進みます。図で表すと以下の通り。

f:id:hi-ctrl:20180119234356p:plain

PID制御のシミュレーション

例題１

　凡例は間違っているので無視してください。では、前節のソースコードでそのまま実行した結果のグラフ（青色）を示します。図1の左は全体、右は拡大図です。例えると、目標値が 0から5[℃]に変更された温度追従制御。操作量であるヒータの出力が一度オーバーシュートするような自動操作で、制御量である温度がおおよそ100[sample]後に目標値（マゼンタ色）に追従しています。更に、ステップ状の手動操作（赤色）に比べてPID制御の方が速く追従することがわかります。拡大図を見てみると、MATLAB／Octaveに実装されているtf関数（連続時間の伝達関数で実行、緑色）との比較がわかります。まずまずです。

f:id:hi-ctrl:20180120003201p:plain

図1：例題１（目標値変更に対する追従制御の時間応答シミュレーション）

例題２

　凡例は間違っているので無視してください。続いて、r=0の変更以外のパラメータはそのままで、今度は外乱を付加してみます。外乱は2番目の緩い傾向をもった外乱で、手動制御（u(k)=一定、赤色）と比較します。平均が0で例えにくいですが、手動制御では温度を一定値に制御できず、あさっての方向に温度が向かっています。一方で、PID制御（青色）だと、目標値0[℃]周りで、制御できていることがわかります。

f:id:hi-ctrl:20180120001612p:plain