長いデッド タイムをもつプロセス制御: スミス予測器

スクリプトを開く

この例では、長いデッド タイムをもつプロセスの PI 制御における制約を示し、「スミス予測器」と呼ばれる制御方法の利点について説明します。

この例は次の文献に基づいて作成されたものです。

A. Ingimundarson and T. Hagglund, "Robust Tuning Procedures of
Dead-Time Compensating Controllers," Control Engineering Practice,
9, 2001, pp. 1195-1208.

プロセス モデル

プロセスの開ループ応答は、40.2 秒の時定数と 93.9 秒の時間遅れをもつデッド タイムを 1 次プロセスに付加したものとしてモデル化されます。

s = tf('s');
P = exp(-93.9*s) * 5.6/(40.2*s+1);
P.InputName = 'u';
P.OutputName = 'y';
P

P =
 
  From input "u" to output "y":
                    5.6
  exp(-93.9*s) * ----------
                 40.2 s + 1
 
Continuous-time transfer function.

遅れが時定数の 2 倍以上になっていることに注目してください。これは、多くの化学的プロセスの典型的なモデルです。このモデルのステップ応答は次のとおりです。

step(P), grid on

PI 制御

比例-積分 (PI) 制御は、プロセス制御においてよく使用されるテクニックです。これに対応する制御アーキテクチャは次に示すとおりです。

補償器 C は、比例ゲイン Kp と積分時間 Ti の 2 つの調整パラメーターをもつ標準型 PI コントローラーです。PIDTUNE コマンドを使用して、0.006 rad/s の開ループ帯域幅をもつ PI コントローラーを設計します。

Cpi = pidtune(P,pidstd(1,1),0.006);
Cpi

Cpi =
 
             1      1 
  Kp * (1 + ---- * ---)
             Ti     s 

  with Kp = 0.0501, Ti = 47.3
 
Continuous-time PI controller in standard form

PI コントローラーの性能を評価するには、フィードバック ループを閉じて、基準信号 ysp と出力外乱信号 d におけるステップ変化への応答をシミュレートします。フィードバック パスに遅れがあるため、SS コマンドを使用して、P または Cpi を状態空間表現に変換する必要があります。

Tpi = feedback([P*Cpi,1],1,1,1);  % closed-loop model [ysp;d]->y
Tpi.InputName = {'ysp' 'd'};

step(Tpi), grid on

閉ループ応答では、オーバーシュートは許容可能ですが、多少の遅れがあります (約 600 秒の整定時間)。比例ゲイン Kp を大きくすると、応答は高速化しますが、オーバーシュートの大幅な増加も招くため、すぐに不安定性につながります。

Kp3 = [0.06;0.08;0.1];      % try three increasing values of Kp
Ti3 = repmat(Cpi.Ti,3,1);   % Ti remains the same
C3 = pidstd(Kp3,Ti3);       % corresponding three PI controllers
T3 = feedback(P*C3,1);
T3.InputName = 'ysp';

step(T3)
title('Loss of stability when increasing Kp')

デッド タイムが長いと、PI 制御の性能が厳しく制限されます。これは、PI コントローラーがデッド タイムを認識せず、実際の出力 y が目的の設定点 ysp と一致しないときに "性急" に反応するためです。これと似た状況は、シャワーを浴びるときに水温調整に時間がかかった場合に誰でも経験したことがあるはずです。この場合、性急に反応したことが原因で、水温を上げすぎて熱湯を浴び、水温を下げすぎて冷水を浴びるといったことを繰り返すことがよくあります。制御方法をより効果的なものとするには、温度設定に加えた変更が有効になるまで待ってから、さらなる調整を加えることが必要です。また、どのノブ設定で好ましい温度が得られるかがわかったら、シャワーが反応するのにかかる時間と同時に正しい温度を実現できるようになります。次に説明するスミス予測器は、この "最適" な制御方法を基本的な考えとしています。

スミス予測器

スミス予測器とは、次の図に示す制御構造のことです。

スミス予測器は、内部モデル Gp を使用して、プロセスの遅延なし応答 yp を予測します (特定のノブ設定から水温がどうなるかなど)。次に、この予測値 yp をターゲット設定点 ysp と比較して、どのような調整が必要であるか (制御量 u) を決定します。ドリフトを防ぎ外乱を抑圧するため、スミス予測器は、デッド タイムを考慮に入れた予測値 y1 と実際のプロセス出力の比較も行います。ギャップ dy=y-y1 はフィルター F を通過してフィードバックされ、全体的な誤差信号 e に寄与します。dy は、シャワーが反応するまで十分な時間を "待ってから"、温度の知覚差異と等しくなります。

スミス予測器スキームを展開するには、以下が必要となります。

  • プロセス ダイナミクスのモデル Gp、およびプロセス デッド タイムの推定値 tau

  • 補償器とフィルター ダイナミクスの適切な設定 (C と F)

プロセス モデルに基づいて、次の式を使用します。

$$G_p(s) = {5.6 \over 1 + 40.2 s } , \tau = 93.9 $$

F には、低周波数外乱を取り込むための 20 秒の時定数をもつ 1 次フィルターを使用します。

F = 1/(20*s+1);
F.InputName = 'dy';
F.OutputName = 'dp';

C については、PI コントローラーによって認識されるプラント全体を使用して、PI コントローラーを再設計します。これには、P、Gp、F、デッド タイムからのダイナミクスが含まれます。スミス予測器制御構造を使用すれば、開ループ帯域幅を上げて高速応答を実現したり、位相余裕を大きくしてオーバーシュートを削減したりすることができます。

% Process
P = exp(-93.9*s) * 5.6/(40.2*s+1);
P.InputName = 'u';
P.OutputName = 'y0';

% Prediction model
Gp = 5.6/(40.2*s+1);
Gp.InputName = 'u';
Gp.OutputName = 'yp';

Dp = exp(-93.9*s);
Dp.InputName = 'yp'; Dp.OutputName = 'y1';

% Overall plant
S1 = sumblk('ym = yp + dp');
S2 = sumblk('dy = y0 - y1');
Plant = connect(P,Gp,Dp,F,S1,S2,'u','ym');

% Design PI controller with 0.08 rad/s bandwidth and 90 degrees phase margin
Options = pidtuneOptions('PhaseMargin',90);
C = pidtune(Plant,pidstd(1,1),0.08,Options);
C.InputName = 'e';
C.OutputName = 'u';
C

C =
 
             1      1 
  Kp * (1 + ---- * ---)
             Ti     s 

  with Kp = 0.574, Ti = 40.2
 
Continuous-time PI controller in standard form

PI コントローラーとスミス予測器の比較

2 つの設計の性能を比較するには、まず、ysp,d から y への閉ループ伝達関数をスミス予測器アーキテクチャ用に導出します。関連するブロックすべての連結作業を容易にするため、入力チャネルと出力チャネルのすべてに名前を付け、CONNECT を使用して配線を自動的に行います。

% Assemble closed-loop model from [y_sp,d] to y
Sum1 = sumblk('e = ysp - yp - dp');
Sum2 = sumblk('y = y0 + d');
Sum3 = sumblk('dy = y - y1');
T = connect(P,Gp,Dp,C,F,Sum1,Sum2,Sum3,{'ysp','d'},'y');

STEP を使用して、スミス予測器 (青) を PI コントローラー (赤) と比較します。

step(T,'b',Tpi,'r--')
grid on
legend('Smith Predictor','PI Controller')

スミス予測器の方が応答が格段に速く、オーバーシュートが生じません。また、ysp から y への閉ループ ボード線図をプロットすると、周波数領域の違いも明らかになります。スミス予測器の方が帯域幅が高いことに注目してください。

bode(T(1,1),'b',Tpi(1,1),'r--',{1e-3,1})
grid on
legend('Smith Predictor','PI Controller')

モデル不一致のロバスト性

前の節で示した解析では、内部モデル

$$ G_p(s) e^{-\tau s} $$

は、プロセス モデル P と完全に一致しました。実際の状況では、内部モデルは真のプロセス ダイナミクスの近似にすぎないため、プロセス ダイナミクスとデッド タイムの不確かさに対してスミス予測器がどの程度ロバストであるかを把握しておくことが重要です。

摂動の影響を受けるプラント モデルが 2 つあり、プロセス パラメーターに対する次の範囲の不確かさをそれぞれ表しているとします。

P1 = exp(-90*s) * 5/(38*s+1);
P2 = exp(-100*s) * 6/(42*s+1);

bode(P,P1,P2), grid on
title('Nominal and Perturbed Process Models')

ロバスト性を解析するには、ノミナル モデルと摂動モデルをプロセス モデルの配列に収集して、PI 設計およびスミス予測器設計の閉ループ伝達関数をリビルドし、閉ループ応答をシミュレートします。

Plants = stack(1,P,P1,P2);  % array of process models
T1 = connect(Plants,Gp,Dp,C,F,Sum1,Sum2,Sum3,{'ysp','d'},'y'); % Smith
Tpi = feedback([Plants*Cpi,1],1,1,1);   % PI

step(T1,'b',Tpi,'r--')
grid on
legend('Smith Predictor 1','PI Controller')

どちらの設計も、閉ループ ボード線図で確認できるとおり、モデルの不一致の影響を受けます。

bode(T1(1,1),'b',Tpi(1,1),'r--')
grid on
legend('Smith Predictor 1','PI Controller')

ロバスト性の改善

モデル化誤差に対するスミス予測器の感度を下げるには、内側のループと外側のループの安定余裕をチェックします。内側のループ C には開ループ伝達関数 C*Gp があるため、安定余裕は次によって取得されます。

margin(C * Gp)
title('Stability Margins for the Inner Loop (C)')

内側のループには快適なゲイン余裕と位相余裕があるため、次に外側のループに注目します。CONNECT を使用して、内側のループを閉じた状態で ysp から dp への開ループ伝達関数 L を導出します。

Sum1o = sumblk('e = ysp - yp');  % open the loop at dp
L = connect(P,Gp,Dp,C,F,Sum1o,Sum2,Sum3,{'ysp','d'},'dp');

bodemag(L(1,1))

この伝達関数は本質的にゼロです。これは、プロセス モデルと予測モデルが完全に一致する場合に予想される値です。外側のループの安定余裕を詳しく把握するには、摂動プロセス モデル (P1 など) のいずれかを使用する必要があります。

H = connect(Plants(:,:,2),Gp,Dp,C,Sum1o,Sum2,Sum3,{'ysp','d'},'dy');
H = H(1,1);  % open-loop transfer ysp -> dy
L = F * H;

margin(L)
title('Stability Margins for the Outer Loop (F)')
grid on;
xlim([1e-2 1]);

このゲイン曲線には 0.04 rad/sec の近傍に凸部があり、これが原因でゲイン余裕が下がり、閉ループ ステップ応答の隆起が大きくなっています。この問題を解決するには、早期により速くロールオフするフィルター F を選択します。

F = (1+10*s)/(1+100*s);
F.InputName = 'dy';
F.OutputName = 'dp';

0.04 rad/sec の位相交差の近傍でゲイン余裕が改善されたことを確認します。

L = F * H;
margin(L)
title('Stability Margins for the Outer Loop with Modified F')
grid on;
xlim([1e-2 1]);

最後に、修正されたフィルターを使用して、閉ループ応答をシミュレートします。

T2 = connect(Plants,Gp,Dp,C,F,Sum1,Sum2,Sum3,{'ysp','d'},'y');

step(T2,'b',Tpi,'r--')
grid on
legend('Smith Predictor 2','PI Controller')

変更後の設計では、ノミナル応答は多少遅くなりますが、より一貫性の高い性能が実現します。

外乱の抑制の改善

d から y への閉ループ伝達関数の式では、F の最適な選択肢が次のように示されています。

$$ F(s) = e^{\tau s} $$

ここで、tau は内部モデルのデッド タイムです。この選択肢では、P と Gp の不一致にかかわらず、外乱を完全に抑制できます。ただし、これらの "負の遅れ" は問題の原因ではないため、実装できません。次の文献

  Huang, H.-P., et al., "A Modified Smith Predictor with an Approximate
  Inverse of Dead Time," AiChE Journal, 36 (1990), pp. 1025-1031

の中で著者は、位相進み近似を使用することを提案しています。

$$ e^{\tau s} \approx { 1 + B(s) \over 1 + B(s) e^{-\tau s} }$$

ここで、B は内部モデル Gp と同じ時定数をもつローパス フィルターです。このスキームは次のようにテストできます。

B(s) と F(s) を定義します。

B = 0.05/(40*s+1);
tau = totaldelay(Dp);
F = (1+B)/(1+B*exp(-tau*s));
F.InputName = 'dy';
F.OutputName = 'dp';

帯域幅が狭い PI コントローラーを再設計します。

Plant = connect(P,Gp,Dp,F,S1,S2,'u','ym');
C = pidtune(Plant,pidstd(1,1),0.02,pidtuneOptions('PhaseMargin',90));
C.InputName = 'e';
C.OutputName = 'u';
C

C =
 
             1      1 
  Kp * (1 + ---- * ---)
             Ti     s 

  with Kp = 0.144, Ti = 40.1
 
Continuous-time PI controller in standard form

閉ループ モデル T3 を計算します。

T3 = connect(Plants,Gp,Dp,C,F,Sum1,Sum2,Sum3,{'ysp','d'},'y');

T3 と T2 および Tpi と比較します。

step(T2,'b',T3,'g',Tpi,'r--')
grid on
legend('Smith Predictor 2','Smith Predictor 3','PI Controller')

この比較が示すとおり、最後の設計では、設定点の追従が低速化する代わりに外乱の抑制を高速化できます。