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

スクリプトを開く

この例では、化学反応器の低変換率から高変換率への遷移のためのゲインスケジュールコントローラーの設計および調整方法を示します。詳細については、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
opt = findopOptions(DisplayReport="off");
[op,report] = findop("rct_CSTR_OL",opspec,opt);
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(3,1,1)
plot(CrEQ,"b-*")
grid
title("Residual concentration")
ylabel("CrEQ")
subplot(3,1,2)
plot(TrEQ,"b-*")
grid
title("Reactor temperature")
ylabel("TrEQ")
subplot(3,1,3)
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
bodeplot(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(2,1,1)
plot(t,C)
grid
set(gca,"ylim",[0 10])
title("Target residual concentration")
ylabel("Cref")
subplot(2,1,2)
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 の二次多項式としてパラメーター化します。

$K_p(C_r) = K_{p0} + K_{p1} C_r + K_{p2} C_r^2 .$

調整する変数の数が減るだけでなく、この方法によって Cr の変動に従う滑らかなゲイン遷移が確保されます。systune を使用して、上で計算された 5 つの平衡点において要件 R1-R4 を満たすように、係数 $K_{p0}, K_{p1}, K_{p2}, K_{i0}, \ldots$ を自動的に調整できます。これは、Cref の軌跡に沿った 5 つの設計点におけるゲインスケジュールコントローラーの調整に相当します。tunableSurface オブジェクトを使用して、各ゲインを Cr の二次関数としてパラメーター化します。"調整グリッド" は 5 つの濃度 CrEQ に対して設定され、二次パラメーター化の基底関数は $C_r, C_r^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);

上記で導入された二次ゲインスケジュールの調整のために、slTuner インターフェイスを作成します。ブロック置換を使用して、非線形プラントモデルを、設計点 CrEQ で取得した 5 つの離散化線形モデル Gd に置き換えます。setBlockParam を使用して、調整可能なゲイン関数 Kp、Ki、Kt、a、b を同名の 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.23, Hard = 0.99976, Iterations = 198

結果の設計は厳密な制約を満たし (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(2,1,1)
lsimplot(getIOTransfer(ST,"Cref","Cr"),uC)
grid
set(gca,"ylim",[-1.5 0.5])
title("Residual concentration")
subplot(2,1,2)
lsimplot(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(3,2,1)
viewSurf(TGS.Kp)
ylabel("Kp")
subplot(3,2,2)
viewSurf(TGS.Ki)
ylabel("Ki")
subplot(3,2,3)
viewSurf(TGS.Kt)
ylabel("Kt")
subplot(3,2,4)
viewSurf(TGS.a)
ylabel("a")
subplot(3,2,5)
viewSurf(TGS.b)
ylabel("b")

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

writeBlockValue(ST)

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

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

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

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

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.23, Hard = 0.99997, Iterations = 204

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

clf
subplot(3,2,1)
viewSurf(setBlockValue(Kp,T))
ylabel("Kp")
subplot(3,2,2)
viewSurf(setBlockValue(Ki,T))
ylabel("Ki")
subplot(3,2,3)
viewSurf(setBlockValue(Kt,T))
ylabel("Kt")
subplot(3,2,4)
viewSurf(setBlockValue(a,T))
ylabel("a")
subplot(3,2,5)
viewSurf(setBlockValue(b,T))
ylabel("b")

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

参考

slTuner | tunableSurface | setBlockParam

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