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

　PID制御に限らず何かしらの計算アルゴリズムの実装を考えたとき、実績のあるパッケージなりを利用するのが鉄則だと思います。しかし、自分で実装してこそ制御則の理解が深まるときもあると思います。PID制御は半世紀前から現在に至るまで使われ続けており、その周辺技術もたくさんあります。多様なアルゴリズムをもつPID制御を使い
こなし、オリジナルの工夫をやるためにも、基本的なPID制御則の実装の話をします。　2018.03.17に記事をもう少し整理するために加筆・修正しました。

キーワード：PID制御、実装、位置型、速度型

ディジタルPIDコントローラ

　黒板や教科書での議論、あるいはアナログ回路による設計で使われるPIDコントローラの制御則は、以下のような連続時間領域の表現で与えられる。

\begin{align} u(t) &= k_c\biggl\{e(t) +\frac{1}{T_i}\int_0^{t}e(\mbox{τ})d\mbox{τ} + T_d\frac{de(t)}{dt}\biggr\} \label{eq:eq1_20161208} \end{align}

ここで、\(u(t)\)は現時刻\(t\)における操作量（MV: Manipulate Value）、\(e(t)\)は制御偏差（E: Error）で目標値\(r(k)\)（SV: Set Value）と制御量\(y(k)\)（PV: Process Value）の差\(r(t)-y(t)\)である。ところで、汎用計算機やマイコン等のディジタル回路にプログラムとして（\(\ref{eq:eq1_20161208}\)）式を実装しようとすると、積分動作や微分動作の実現方法を考えなければならない。ディジタル制御に限らないが、このようなとき離散化、つまり近似で実現する。では、プログラムに実装するために離散化したPIDコントローラの制御則を以下に示す。

\begin{align} u(k) &= k_c \biggl\{e(k) + \frac{1}{T_i}\sum_{\mbox{τ}_k=0}^{k} e(\mbox{τ}_k)T_s + T_d \frac{e(k) - e(k - 1)}{T_s} \biggr\} \label{eq:eq2_20161208} \end{align}

（\ref{eq:eq2_20161208}）式の第二項と第三項は、それぞれ離散化した積分動作と微分動作であるが、これは長方形近似と後進差分近似と呼ばれるものを使った。どちらも微小時間\(T_s\)を無限小に飛ばす前の積分と微分の定義式のことであり、ディジタル制御の場合は\(T_s\)をサンプリング時間あるは制御周期とする。また、時刻\(t\)に関して連続の関数となる比例動作や\(u(t)\)、\(e(t)\)については、\(t=kT_s\)のタイミングでサンプリングした離散値列\(u(k)\)、\(e(k)\)を用いる（\(T_s\)は表記を省略）。ちなみに、\(T_s\)は離散化の精度に直結するので、制御対象の応答性に対して相対的に大きな値となれば、おおざっぱで不正確な積分、微分、離散化列になってしまうので、制御性能が失われる原因の一つとなる。

位置型PID制御と速度型PID制御

　（\ref{eq:eq2_20161208}）式は過去や現在の制御偏差\(e(k)\)のみに基づいて操作量\(u(k)\)の絶対量を過去の\(u\)に関係なく計算する。このようなPID制御を位置型PID制御と呼ぶが、プログラムとして実装するに当たっては、主に以下のような扱いづらい点がある。

手動から自動へのモード変更やコントローラの切り替え時に不連続な変化（bump）が発生する。
\(u(k)\)のままだと諸所の演算処理（リセットワインドアップや変化量制限）が難しい。
半永久的な制御で\(\sum\)の計算のためにメモリが肥大化する（メモリが小さいマイコンなどで問題となる）。

そこで、\(u(k)\)を直接計算するのではなく、（\(\ref{eq:eq4_20161208}\)）式のような変化量\(\Delta u(k)\)をまずは計算する。

\begin{align} \Delta u(k) &= k_c \biggl\{e(k)-e(k - 1) + \frac{T_s}{T_i}e(k) + \frac{T_d}{T_s} \{e(k) - 2e(k - 1) + e(k - 2)\} \biggr\} \label{eq:eq4_20161208} \end{align}

そして、本来必要となる\(u(k)\)は、メモリに保存しておいた前回値\(u(k - 1)\)から

\begin{align} u(k) &= u(k - 1) + \Delta u(k) \label{eq:eq5_20161208} \end{align}

と計算するのである。このような前回値\(u(k - 1)\)を始点に、次に必要となる変化量\(\Delta u(k)\)分だけを計算・加算することで現在の\(u(k)\)を決定するPID制御を速度型PID制御と呼ぶ。主にプログラムとして実装する際は（\ref{eq:eq4_20161208}）と（\ref{eq:eq5_20161208}）を採用するのが良いと思われる。この計算の工夫により、バンプレスな切り替え（図1）、リセットワインドアップ、メモリ節約が可能となる。

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

図1：手動⇒自動の切り替えに関する位置型PIDと速度型PIDの比較
（k=0~25まではu(k)=2で手動制御し、それ以降は自動制御に切り替えている。位置型の場合、u(k)は不連続に2から4付近まで変化させて、y(k)がオーバーシュートしているが、速度型だとそのようなことが起きていない）

PIDコントローラの実装とその要領

　位置型PID制御（\ref{eq:eq2_20161208}）と速度型PID制御（\ref{eq:eq4_20161208}）（\ref{eq:eq5_20161208}）をMATLAB／Octaveで実装したソースコードを以下に示す。このソースコードのみではまだ十分に動かないが、主にPIDコントローラの部分でどのような書き方をするのかを参考程度にしてほしい（これが正解で効率的な書き方ではない）。位置型PID制御はないが実際に動くソースコードは、下記の記事を参照願う。

hi-ctrl.hatenablog.com

＜MATLAB／Octaveコード＞

% 設計パラメータ
Ts = 1;             % サンプリング時間
N  = XXX;           % サンプル時刻の終点（シミュレーション回数）

flg = 'v';  % p:position, v:verocity

offset_time = 2;

y  = zeros(N + offset_time, 1);
r  = [zeros(offset_time, 1); ones(N, 1) * 5];
e  = zeros(N + offset_time, 1);
u  = zeros(N + offset_time, 1);

for k = (1 : N) + offset_time
    
    % 制御量の読込み
    y(k) = read_PV();
    
    % 制御偏差の計算
    e(k) = r(k) - y(k);
    
    % PID制御則による操作量の計算
    if strcmp(flg, 'p')
        u(k) = kp * (e(k) + Ts/Ti*sum(e(1:k)) + Td/Ts*(e(k) - e(k-1)));
    else
        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
    
    % 制御量の書き込み
    write_MV(u(k));
    
end

PIDコントローラをMATLAB／Octaveで実装するに当たって、コーディングでは以下の点を留意してほしい。

時系列データをベクトルなどで実装したとき、各演算で過去データを参照する際は、例えば、u(- 1)やe(- 1)のような要素の指定がマイナスになるケースを回避する方法が必要である。上のコードの場合、offset_timeなど。

まとめ

　教科書なんかに載っている数式を実装し、シミュレーションを行いたい時は離散化などの近似が必要である場合がある。今回、PID制御則を離散化する方法を説明し、プログラムとして実装するための位置型PID制御（\ref{eq:eq2_20161208}）と速度型PID制御（\ref{eq:eq4_20161208}）（\ref{eq:eq5_20161208}）を示した。特に、後者の内容に関しては、

パッケージとしてすでに実装されているのは、殆どが速度型アルゴリズム（DCSなんかはそうだったり）
（\ref{eq:eq5_20161208}）式の変化量を計算する考え方はモデル予測制御や他の制御則でも使われるから意外と重要

などの点を補足として加える。