最新のリリースでは、このページがまだ翻訳されていません。 このページの最新版は英語でご覧になれます。

化学反応器のゲイン スケジュール制御

この例では、化学反応器の低変換率から高変換率への遷移のためのゲイン スケジュール コントローラーの設計および調整方法を示します。詳細については、Seborg, D.E. et al., "Process Dynamics and Control", 2nd Ed., 2004, Wiley, pp. 34-36 を参照してください。

連続攪拌タンク反応器

ここで取り上げるプロセスは、低変換率から高変換率 (高残留濃度から低残留濃度) への遷移中の連続攪拌タンク反応器 (CSTR) です。化学反応は発熱する (熱を生成する) ため、熱暴走を防ぐために反応器の温度を制御しなければなりません。制御タスクは、プロセス ダイナミクスは非線形であり、変換率が高くなるにつれて安定から不安定へ、そして安定に戻るというように遷移するという事実によって複雑になります。反応器のダイナミクスは Simulink でモデル化されます。制御変数は残留濃度 Cr および反応器の温度 Tr であり、操作変数は反応器の冷却被覆を循環する冷却水の温度 Tc です。

open_system('rct_CSTR_OL')

8.57 kmol/m^3 の残留濃度から遷移して、最初に 2 kmol/m^3 まで下げます。プロセス ダイナミクスが残留濃度 Cr に伴って変化する様子を理解するには、8.57 から 2 の間の Cr の 5 つの値について平衡条件を特定し、各平衡のプロセス ダイナミクスを線形化します。各平衡点における反応器と冷却水の温度をログに記録します。

CrEQ = linspace(8.57,2,5)';  % concentrations
TrEQ = zeros(5,1);           % reactor temperatures
TcEQ = zeros(5,1);           % coolant temperatures

% Specify trim conditions
opspec = operspec('rct_CSTR_OL',5);
for k=1:5
   % Set desired residual concentration
   opspec(k).Outputs(1).y = CrEQ(k);
   opspec(k).Outputs(1).Known = true;
end

% Compute equilibrium condition and log corresponding temperatures
[op,report] = findop('rct_CSTR_OL',opspec,...
   findopOptions('DisplayReport','off'));
for k=1:5
   TrEQ(k) = report(k).Outputs(2).y;
   TcEQ(k) = op(k).Inputs.u;
end

% Linearize process dynamics at trim conditions
G = linearize('rct_CSTR_OL', 'rct_CSTR_OL/CSTR', op);
G.InputName = {'Cf','Tf','Tc'};
G.OutputName = {'Cr','Tr'};

各平衡点における濃度の関数として反応器と冷却水の温度をプロットします。

subplot(311), plot(CrEQ,'b-*'), grid, title('Residual concentration'), ylabel('CrEQ')
subplot(312), plot(TrEQ,'b-*'), grid, title('Reactor temperature'), ylabel('TrEQ')
subplot(313), plot(TcEQ,'b-*'), grid, title('Coolant temperature'), ylabel('TcEQ')

Cr=8.57 および Cr=2 の平衡の間を滑らかに遷移するために、開ループ制御戦略は、上記の冷却水温度のプロファイルに従うことから成ります。ただし、この戦略は、中間の範囲での反応は不安定であり、熱暴走を防ぐには適切に冷却しなければならないという事実によって失敗します。これは、上記で取り上げた 5 つの平衡点の線形化されたモデルの極を検査することで確認されます (5 つのモデルのうち 3 つが不安定)。

pole(G)
ans(:,:,1) =

  -0.5225 + 0.0000i
  -0.8952 + 0.0000i


ans(:,:,2) =

   0.1733 + 0.0000i
  -0.8866 + 0.0000i


ans(:,:,3) =

   0.5114 + 0.0000i
  -0.8229 + 0.0000i


ans(:,:,4) =

   0.0453 + 0.0000i
  -0.4991 + 0.0000i


ans(:,:,5) =

  -1.1077 + 1.0901i
  -1.1077 - 1.0901i

ボード線図では、Cr=8.57 から Cr=2 に遷移する間の、プラント ダイナミクスにおける大きな変動がさらに強調表示されます。

clf, bode(G(:,'Tc'),{0.01,10})

フィードバック制御戦略

残留濃度を低下させながら熱暴走を防ぐには、フィードバック制御を使用して残留濃度 Cr および反応器の温度 Tr の測定値に基づいて冷却水温度 Tc を調整します。このアプリケーションの場合、内側のループが反応器の温度を調整して外側のループが濃度の設定点を追従するカスケード制御のアーキテクチャを使用します。両方のフィードバック ループは 0.5 分のサンプリング周期のデジタルです。

open_system('rct_CSTR')

ターゲット濃度 Cref は t=10 で 8.57 kmol/m^3 から t=36 で 2 kmol/m^3 へ低下します (遷移は 26 分間続きます)。反応器の温度に対応するプロファイル Tref は、平衡化された解析から平衡値 TrEQ を内挿することで取得されます。コントローラーは、Cr=8.57 の場合の初期平衡値 TcEQ(1)=297.98 に対して、冷却水温度調整 dTc を計算します。最初は、"Concentration controller" ブロックの出力 TrSP が反応器の温度に一致して、調整 dTc はゼロであり、冷却水温度 Tc が平衡値 TcEQ(1) となるようにモデルが設定されることに注意してください。

clf
t = [0 10:36 45];
C = interp1([0 10 36 45],[8.57 8.57 2 2],t);
subplot(211), plot(t,C), grid, set(gca,'ylim',[0 10])
title('Target residual concentration'), ylabel('Cref')
subplot(212), plot(t,interp1(CrEQ,TrEQ,C))
title('Corresponding reactor temperature at equilibrium'), ylabel('Tref'), grid

制御目的

TuningGoal オブジェクトを使用して、設計要件を取得します。まず、Cr は応答時間が約 5 分の設定点 Cref に従う必要があります。

R1 = TuningGoal.Tracking('Cref','Cr',5);

内側のループ (温度) は、必要に足るダンピングと十分に速い減衰によって反応のダイナミクスを安定させる必要があります。

MinDecay = 0.2;
MinDamping = 0.5;
% Constrain closed-loop poles of inner loop with the outer loop open
R2 = TuningGoal.Poles('Tc',MinDecay,MinDamping);
R2.Openings = 'TrSP';

コントローラー出力での Rate Limit ブロックは、冷却水温度 Tc が 1 分間に 10 度を越える速さで変化できないことを指定します。これはコントローラーの影響力に対する深刻な制限であり、これが無視されると性能の低下や不安定を招きかねません。このレート制限を考慮するために、0.25 kmol/m^3/min で変化する Cref を観察します。Tc が 10 度/分を越える速さで変化しないことを確認するには、Cref から Tc へのゲインは 10/0.25=40 未満である必要があります。

R3 = TuningGoal.Gain('Cref','Tc',40);

最後に、プラント入力 Tc において少なくとも 7 dB のゲイン余裕と 45 度の位相余裕が必要です。

R4 = TuningGoal.Margins('Tc',7,45);

ゲイン スケジュール コントローラー

これらの要件を達成するには、PI コントローラーを外側のループで、進み補償器を内側のループで使用します。サンプリング レートが遅いため、中間の範囲の濃度 Cr = 5.28 kmol/m^3/min での化学反応を適切に安定化させるために進み補償器が必要です。反応のダイナミクスは濃度によって大幅に変化するため、コントローラー ゲインを濃度の関数としてさらにスケジュールします。図 1 および図 2 に示すように、これは Lookup Table ブロックを使用して Simulink でモデル化されます。

図 1: 濃度ループのゲイン スケジュール PI コントローラー

図 2: 温度ループのゲイン スケジュール進み補償器

このゲイン スケジュール コントローラーの調整は、濃度値の範囲でのルックアップ テーブル データの調整に相当します。個々のルックアップ テーブルのエントリを調整するのではなく、コントローラー ゲイン Kp,Ki,Kt,a,b を、たとえば Cr の 2 次多項式としてパラメーター化します。

調整する変数の数が減るだけでなく、この方法によって Cr の変動に従う滑らかなゲイン遷移が確保されます。systune を使用して、上で計算された 5 つの平衡点において要件 R1-R4 を満たすように、係数 を自動的に調整できます。これは、Cref 軌跡に沿った 5 つの設計点におけるゲイン スケジュール コントローラーの調整に相当します。tunableSurface オブジェクトを使用して、各ゲインを Cr の 2 次関数としてパラメーター化します。"調整グリッド" は 5 つの濃度 CrEQ に対して設定され、2 次パラメーター化の基底関数は です。大部分のゲインは完全に同じようにゼロに初期化されています。

TuningGrid = struct('Cr',CrEQ);
ShapeFcn = @(Cr) [Cr , Cr^2];

Kp = tunableSurface('Kp', 0, TuningGrid, ShapeFcn);
Ki = tunableSurface('Ki', -2, TuningGrid, ShapeFcn);
Kt = tunableSurface('Kt', 0, TuningGrid, ShapeFcn);
a = tunableSurface('a', 0, TuningGrid, ShapeFcn);
b = tunableSurface('b', 0, TuningGrid, ShapeFcn);

コントローラーの調整

ターゲットの帯域幅がナイキスト周波数の 1 decade 内に収まっているため、コントローラーを離散領域内で直接調整するのは簡単です。サンプル時間が 0.5 分の線形化プロセス ダイナミクスを離散化します。デジタル コントローラーが連続時間プラントと通信を行う方法を反映するために、ZOH メソッドを使用します。

Ts = 0.5;
Gd = c2d(G,Ts);

上記で導入された 2 次ゲイン スケジュールの調整のために slTuner インターフェイスを作成します。ブロック置換を使用して、非線形プラント モデルを設計点 CrEQ で取得された 5 つの離散化線形モデル Gd に置き換えます。setBlockParam を使用して調整可能なゲイン関数 KpKiKtab を、同じ名前の Lookup Table ブロックに関連付けます。

BlockSubs = struct('Name','rct_CSTR/CSTR','Value',Gd);
ST0 = slTuner('rct_CSTR',{'Kp','Ki','Kt','a','b'},BlockSubs);
ST0.Ts = Ts;  % sample time for tuning

% Register points of interest
ST0.addPoint({'Cref','Cr','Tr','TrSP','Tc'})

% Parameterize look-up table blocks
ST0.setBlockParam('Kp',Kp);
ST0.setBlockParam('Ki',Ki);
ST0.setBlockParam('Kt',Kt);
ST0.setBlockParam('a',a);
ST0.setBlockParam('b',b);

systune を使用して、要件 R1-R4 に対してコントローラー係数を調整します。安定余裕要件を厳密な制約にして、残りの要件を最適化します。

ST = systune(ST0,[R1 R2 R3],R4);
Final: Soft = 1.21, Hard = 0.99887, Iterations = 210

結果の設計は厳密な制約 (Hard<1) を満たし、残りの要件 (1 に近い Soft) をほぼ満たします。この設計を検証するには、Cref と同じ勾配をもつ濃度の低下に対する応答をシミュレートします。各プロットは、5 つの設計点 CrEQ における線形応答を示します。

t = 0:Ts:20;
uC = interp1([0 2 5 20],(-0.25)*[0 0 3 3],t);
subplot(211), lsim(getIOTransfer(ST,'Cref','Cr'),uC)
grid, set(gca,'ylim',[-1.5 0.5]), title('Residual concentration')
subplot(212), lsim(getIOTransfer(ST,'Cref','Tc'),uC)
grid, title('Coolant temperature variation')

冷却水温度の変化率は、物理的制約の範囲内に収まっていることに注意してください (1 分あたり 10 度またはサンプル周期あたり 5 度)。

コントローラーの検証

遷移中に各ゲインが Cr によってどのように変動するかを検査します。

% Access tuned gain schedules
TGS = getBlockParam(ST);

% Plot gain profiles
clf
subplot(321), viewSurf(TGS.Kp), ylabel('Kp')
subplot(322), viewSurf(TGS.Ki), ylabel('Ki')
subplot(323), viewSurf(TGS.Kt), ylabel('Kt')
subplot(324), viewSurf(TGS.a), ylabel('a')
subplot(325), viewSurf(TGS.b), ylabel('b')

Simulink でゲイン スケジュール コントローラーを検証するには、まず writeBlockValue を使用して、調整結果を Simulink モデルに適用します。各 Lookup Table ブロックについて、これはテーブルのブレークポイントの対応する 2 次ゲインの式を評価して、テーブル データを適宜更新します。

writeBlockValue(ST)

次に [再生] ボタンを押して、調整後のゲイン スケジュールによる応答をシミュレートします。シミュレーション結果を図 3 に示します。ゲイン スケジュール コントローラーは、適切な応答時間がありレート制限の飽和がない遷移からの反応を正常に促進します (コントローラー出力 du は有効温度変動 dTc に一致します)。反応器の温度がその平衡値 Tref の近くに留まることから、コントローラーが熱暴走を防ぐ一方で、平衡付近の反応を維持していることが示されます。

図 3: ゲイン スケジュール カスケード コントローラーによる遷移。

MATLAB でのコントローラーの調整

あるいは、slTuner インターフェイスを使用せずに MATLAB で直接ゲイン スケジュールを調整できます。最初に、上記で行ったようにゲインを Cr の 2 次関数としてパラメーター化します。

TuningGrid = struct('Cr',CrEQ);
ShapeFcn = @(Cr) [Cr , Cr^2];

Kp = tunableSurface('Kp', 0, TuningGrid, ShapeFcn);
Ki = tunableSurface('Ki', -2, TuningGrid, ShapeFcn);
Kt = tunableSurface('Kt', 0, TuningGrid, ShapeFcn);
a = tunableSurface('a', 0, TuningGrid, ShapeFcn);
b = tunableSurface('b', 0, TuningGrid, ShapeFcn);

これらのゲインを使用して PI コントローラーおよびリード コントローラーを作成します。

PI = pid(Kp,Ki,'Ts',Ts,'TimeUnit','min');
PI.u = 'ECr';   PI.y = 'TrSP';

LEAD = Kt * tf([1 -a],[1 -b],Ts,'TimeUnit','min');
LEAD.u = 'ETr';   LEAD.y = 'Tc';

connect を使用して、5 つの設計点において制御システム全体の閉ループ モデルを作成します。ループが開かれ、これらの位置で安定余裕が評価されるように、コントローラー出力 TrSP および Tc を "解析ポイント" としてマークします。閉ループ モデル T0 は、Kp,Ki,Kt,a,b の調整可能な係数により異なる線形モデルの 5 行 1 列の配列です。各モデルは離散しており、30 秒ごとにサンプリングされます。

Gd.TimeUnit = 'min';
S1 = sumblk('ECr = Cref - Cr');
S2 = sumblk('ETr = TrSP - Tr');
T0 = connect(Gd(:,'Tc'),LEAD,PI,S1,S2,'Cref','Cr',{'TrSP','Tc'});

最後に、systune を使用してゲイン スケジュール係数を調整します。

T = systune(T0,[R1 R2 R3],R4);
Final: Soft = 1.21, Hard = 0.999, Iterations = 204

結果は、上記で取得されたものと同様です。T の調整後の係数を使用して Cr の関数としてゲインをプロットすることで確認します。

clf
subplot(321), viewSurf(setBlockValue(Kp,T)), ylabel('Kp')
subplot(322), viewSurf(setBlockValue(Ki,T)), ylabel('Ki')
subplot(323), viewSurf(setBlockValue(Kt,T)), ylabel('Kt')
subplot(324), viewSurf(setBlockValue(a,T)), ylabel('a')
subplot(325), viewSurf(setBlockValue(b,T)), ylabel('b')

各設計点において線形応答をシミュレーションすることで設計をさらに検証できます。ただし、Simulink に戻り、ゲイン スケジュール コントローラーの非線形応答をシミュレーションする必要があります。

参考

| |

関連する例

詳細