制御の基本について

PID制御、安定性評価の基本についてまとめていきます。

・PIDの基本

・ラプラス変換の基本

各種伝達関数の単位ステップ応答

・PID Tuning

・安定性評価

・基本シーケンス

PIDの基本

教科書的なPIDの形は以下の基本形になります。

設定値SP(t)と、測定値PV(t)の偏差E(t)=PV(t)-SP(t)とすると、制御出力MV(t)は、以下の式で表されます。

$$MV(t)=Kp\left[E(t)+\frac{1}{Ti}\int E(t) dt + Td \frac{dE(t)}{dt}\right]$$

上記の基本形:Ideal型の他にも以下Series型やPrallel型があります。

$$MV(t)=Kp\left[E(t)+\frac{1}{Ti}\int E(t) dt \right] \left[ 1+ Td \frac{dE(t)}{dt}\right]$$ $$MV(t)=KpE(t)+\frac{1}{Ti}\int E(t) dt + Td \frac{dE(t)}{dt}$$

上記のPIDの基本形は連続時間における表現となっておりますが、実際にはDCSのスキャンタイムは通常0.5sec程度で固定されていますので、実際の制御を模擬する場合には上式を離散化する必要があります。

スキャンタイムをΔtとすると、位置型は以下のようになります。

$$MV(t_n)=Kp\left[E(t_n)+\frac{1}{Ti} \displaystyle \sum_{i=0}^{n}E(t_n) \Delta t + Td \frac{\Delta E(t_n)}{\Delta t}\right]$$ $$ \Delta E(t_n) = E(t_n)-E(t_{n-1})$$

速度型は以下のようになります。

$$ \Delta MV(t_n)=Kp\left[ \Delta E(t_n)+\frac{1}{Ti} E(t_n) \Delta t + Td \frac{\Delta (\Delta E(t_n))}{\Delta t}\right]$$ $$ \Delta E(t_n) = E(t_n)-E(t_{n-1})$$ $$ \Delta (\Delta E(t_n)) = (E(t_n)-E(t_{n-1}))-(E(t_{n-1})-E(t_{n-2}))$$

速度型は、前回の出力値に対して、現在時間における修正分ΔMVを足し合わせて出力値を求めます。

このため、Manual⇒Autoへの切り替え時にも大きなバンプを抑えることができます。

逆に位置型は、PV値が流量計のようにノイズを含みやすい場合に、安定的に制御することができます。(特に、後述のI-PD型の場合)

※参考書籍

Instrument Engineers' Handbook, Volume One: Process Measurement and Analysis (English Edition)

前述の基本形(PID型)に加え、以下の2つの制御方法があります。(位置型で記載しています)

I-PD型(比例微分先行型)

$$MV(t_n)=Kp\left[PV(t_n)+\frac{1}{Ti} \displaystyle \sum_{i=0}^{n}E(t_n) \Delta t + Td \frac{\Delta PV(t_n)}{\Delta t}\right]$$ $$ \Delta PV(t_n) = PV(t_n)-PV(t_{n-1})$$

PI-D型(微分先行型)

$$MV(t_n)=Kp\left[E(t_n)+\frac{1}{Ti} \displaystyle \sum_{i=0}^{n}E(t_n) \Delta t + Td \frac{\Delta PV(t_n)}{\Delta t}\right]$$ $$ \Delta PV(t_n) = PV(t_n)-PV(t_{n-1})$$

比例微分先行型は、目標値変更によるキックを抑制することができ、オーバーシュートも小さくすることができますが、目標値変更に対する追従性が遅いという問題があります。

カスケード制御における2次ループのように、設定値変更に対する追従性を良くしたい場合は、微分先行型を採用した方が良いことになります。

後述の2自由度制御とも重なりますが、PIDは理想的には①「設定値の変更に対して素早く追従すること」に加え、②「外乱が入ってきた際にも素早く制定できること」が望ましいです

上記の2つの目的を1つの制御で実現することは難しいですが、特定の外乱に対して、制御すべき方向がある程度事前に予測できる場合は、フィードフォワードとフィードバックを組み合わせることで制御性能を向上させることができます。

こちらの本に記載されているFigure6.4[Feedback trim for flow-to-flow control of the hot water process]を参考に、以前AVEVA社のDYNSIMを使って作成したモデルを以下参考のために示します。

FFFB

上記の例は冷却水をスチームで加温し温水を作るプロセスにおいて、冷却水の流量や温度が変化した場合に必要となるスチーム量を計算し、フィードフォワードで先行して制御に加算する形としています。

トレンドグラフを見ていただくと、フィードバックのみ(FB only)の場合に比べ、フィードフォワード+フィードバック(FF + FB)の方が温水の温度変動を抑制できていることが確認できます。

なお、上図ではフィードバック制御のみの時と、フィードフォワード+フィードバック制御の切り替えをボタンで行っていますが、DYNSIMではこのようにボタンによって特定のイベントを起こすことが可能になっています。(シナリオという機能もあり、スクリプト形式でイベントを書き下すこともできますが、GUI上から操作できるという点でボタンを使った方が便利なケースもあります)

検出端と操作バルブが離れており、大きな無駄時間があるような系に対して、通常のPID制御だけでコントロールをかけてしまうと、PID出力の結果が反映されるまでの間の誤差も累積されてしまい、ハンチングが生じてしまいます。

こうしたケースでは、一定の周期で制御とサンプリングを繰り返すサンプルPID制御が有効です。

以下はAVEVA社の非定常シミュレーターDYNSIMを用いた簡易検証モデルです。強制的に長い無駄時間を設けた熱量調整モデルですが、サンプルPID制御を行わなっていないNorPIDと比べて、制御時間2秒、サンプリング時間20秒としたサンプルPIDの方が安定的に制御できていることが確認できます。

samplePID

PIDは、理想的には①「設定値の変更に対して素早く追従すること」に加え、②「外乱が入ってきた際にも素早く制定できること」が望ましいです。

①、②を1自由度で同時に満たすことはできないため、1自由度で制御する場合はどちらかを犠牲にする形になります。

2自由度制御

ラプラス変換の基本

関数f(t)のラプラス変換

$$\mathfrak{L}\left[ f( t) \right] = \int_0^\infty [f( t ){e^{ - st}}dt]$$

原関数(original function) f(t)のラプラス変換

$$\mathfrak{L}\left[ f( t) \right] = F( s )$$

像関数(image function) F(s)の逆ラプラス変換

$$\mathfrak{L}^-1\left[ F( s) \right] = f( t )$$

各種伝達関数の単位ステップ応答

各プログラミング言語ごとの応答の確認はそれぞれ以下をご参照ください。

VBA Scilab

1次遅れプロセスの例として、Perry Handbook8版の8章を参照に以下記載します。

A成分が一定流量F[m3/sec]、濃度\(C_{in}\)[kg-mol/m3]で反応器に入り、反応速度定数k[1/sec]でB成分に変化する反応を考えます。

ある時間における反応器内のA成分の濃度を$C_A$[kg-mol/m3]とすると、物質収支は以下となります。

$$V\frac{dC_A}{dt}=FC_{in}-FC_A-kVC_A$$

液相反応でボリュームは非圧縮性、反応速度定数は一定(F,V,k一定)とすると、両辺をF+kVで割り、

$$\frac{V}{F+kV}\frac{dC_A}{dt}=\frac{F}{F+kV}C_{in}-C_A$$

V/(F+kV)をτ、F/(F+kV)をKとすると、上式は以下となります。

$$\tau\frac{dC_A}{dt}=KC_{in}-C_A$$

上式でτが時定数、Kがプロセスゲインとなります。

時定数τ、プロセスゲインKの1次遅れプロセスの伝達関数は以下となります。

$$G(s)=\frac{K}{\tau s +1}$$

1次応答プロセスに対するインパルス応答、ステップ応答は以下となります。

インパルス応答 $$f(t)=\mathcal{L^{-1}}\left[G(s)U(s)\right]$$

U(s)=1なので、

$$f(t)=\frac{K}{\tau}\mathcal{L^{-1}}\left[\frac{1}{s+1/ \tau}\right]=\frac{K}{\tau}e^{-t/\tau}$$ ステップ応答 $$f(t)=\mathcal{L^{-1}}[G(s)U(s)]$$

U(s)=1/sなので、

$$f(t)=K\mathcal{L^{-1}}\left[\frac{1}{s \cdot (\tau s+1)}\right]=K\mathcal{L^{-1}}\left[\frac{1}{s}-\frac{1} {s +1/ \tau}\right]$$ $$f(t)=K(1-e^{-t/\tau})$$

逆ラプラス変換については、こちらのサイトが非常に丁寧で分かりやすいです。

PID Tuning

限界感度法は無駄時間L/時定数Tが1より小さい範囲で適用します。

無駄時間が支配的なプロセスや、セットポイントのトラッキングテスト等では有効な効果が得られません。

限界感度法参照テーブル

Kp Ti Td
P 0.5Kc - -
PI 0.45Kc Pu/1.2 -
PID 0.6Kc 0.5Pu Pu/8

過渡応答法参照テーブル

Kp Ti Td
P 1/RL - -
PI 0.9/RL L/0.3 -
PID 1.2/RL 2L 0.5L

IMC法参照テーブル(1次遅れ系)

Kp Ti Td
PI (Tp+L/2)/Kp(F+L/2) Tp+L/2 -
PID (Tp+L/2)/Kp(F+L/2) Tp+L/2 (Tp・L/2)/(Tp+L/2)

安定性評価について

Bode線図の基本をまとめていきます。

以下のシンプルな伝達関数を例として、Bode線図の作成方法をご紹介します。

$$G(s)=\frac{1}{5s+1}$$

上式に\(s=\omega j \)を代入します。

$$G(\omega j)=\frac{1}{5\omega j+1}$$

実数部:G(Re)と虚数部:G(Im)を求めます。

$$G(Re)=\frac{1}{25\omega^2+1}$$ $$G(Im)=\frac{-5\omega j}{25\omega^2+1}$$

ゲイン、位相を求めます。

$$Gain=20 \cdot log10 |G(j\omega)|$$ $$Phase=tan^{-1} \left(\frac{G(Im)}{G(Re)} \right)$$

Nyquist線図の基本をまとめていきます。

Matlab、Scilab、Excelを使ったBode線図、Nyquist線図の作成方法をまとめていきます。

・Bode線図作成ツール

Bode線図の作成ツールをsilverlightで作成したのでご紹介します。

伝達関数、計算条件を下図のように設定します。(2017/6/4現在、1次遅れのみの実装です)

ControlTool

Bode線図が確認できます。

ControlTool

各データの出力結果がテーブル形式で確認できます

ControlTool

閉ループのBode線図、Nyquist線図も実装しました。(2017/6/4現在Impulse応答のみ対応しています)

以下Matalbとの比較になります。

ControlTool ControlTool

・ExcelでのNyquist線図、Bode線図の作成方法

作成方法:準備中

例として以下の伝達関数の各線図を紹介します。

$$\frac{K}{(s+1)(s+2)(s+3)}$$

Bode線図

Nyquist線図

・ScilabでのNyquist線図(重ね書き)、Bode線図の作成方法

・MatlabでのBode線図,Nyquist線図の作成方法

transfer functionの作成 : tf(Numerator,Denominator)

ex)tf([0, 1],[5 1])

Bode線図の作成:bode(transfer function)

ex)bode(tf([0 1],[5 1]))

Nyquist線図の作成:nyquist(transfer function)

ex)nyquist(tf([0 1],[5 1]))

システム同定の基本

・予測誤差法(Prediction Error Method)

・ARXモデル(Auto-Regressive with extra input model)

基本シーケンスについて

AND, ORユニット LATCHユニット

LATCHユニットはAND/ORユニットと比べると直感的に把握しにくいかもしれませんが、Advanced Digital Design with the Verilog HDLが素晴らしく分かりやすいです。

当サイトに不具合、ご意見等ございましたらCEsolutionにお知らせください。

Page Top