获得化学反应器的计划控制

此示例使用：

金宝app 金宝app

打开脚本

本示例说明了如何为从低转化率过渡到高转化率的化学反应器设计和调整增益计划控制器。有关背景信息，请参阅Seborg，D.E.等人，“过程动力学和控制”，第二版，2004年，Wiley，第34-36页。

连续搅拌槽式反应器

这里考虑的过程是一个连续搅拌釜式反应器（CSTR），在从低转化率到高转化率（高到低残余浓度）的过渡过程中。因为化学反应是放热的（产生热量），必须控制反应器温度以防止热失控。由于过程动力学是非线性的，并且随着转化率的增加，从稳定到不稳定再回到稳定，控制任务变得复杂。反应器动力学在Simulink中建模。控制变量为剩余浓度离子金宝appCr以及反应堆温度Tr，而被操纵的变量是温度Tc在反应堆冷却套中循环的冷却剂。

开放式系统(“rct_CSTR_OL”）

我们想从最初的残留浓度8.57 kmol/m^3降到2 kmol/m^3。了解工艺动态如何随着残留浓度的变化而变化Cr，求五个值的平衡条件Cr在8.57和2之间，并在每个平衡周围线性化过程动力学。记录反应堆和冷却剂在每个平衡点的温度。

CrEQ=linspace（8.57,2,5）'；%浓度TrEQ = 0(5、1);%反应堆温度TcEQ = 0(5、1);%冷却液温度%指定修剪条件opspec = operspec (“rct_CSTR_OL”5);为k=1:5%设置所需的残留浓度opspec (k) .Outputs(1)。y = CrEQ (k);opspec (k) .Outputs(1)。知道= true;结束计算平衡条件并记录相应温度(op、报告)= findop (“rct_CSTR_OL”opspec,．..findopOptions (“DisplayReport”,“关”));为(2).y =1;TcEQ (k) = op (k) .Inputs.u;结束%在微调条件下线性化过程动力学G=线性化(“rct_CSTR_OL”,“rct_CSTR_OL /装运箱”, op);G.InputName = {Cf的,“助教”,“Tc”};G.OutputName = {“Cr”,“Tr”};

绘制反应堆和冷却剂的平衡温度作为浓度的函数。

次要情节(311),情节(CrEQ“b - *”)、网格、标题(的残余浓度), ylabel (“CrEQ”)次要情节(312)、情节(TrEQ“b - *”)、网格、标题(“反应堆温度”), ylabel (“特雷克”)子地块（313）、地块（TcEQ、，“b - *”)、网格、标题(冷却液温度的), ylabel (“TcEQ”）

一种开环控制策略包括遵循冷却剂温度分布以上平滑过渡之间Cr= 8.57,Cr= 2平衡。然而，这一策略注定要失败，因为反应在中程是不稳定的，必须适当冷却以避免热失控。这是通过考察上述五个均衡点的线性化模型的极点来证实的(五个模型中有三个是不稳定的)。

极(G)

Ans (:，:，3) = 0.5114 + 0.00000 i -0.8229 + 0.00000 i Ans (:，:，4) = 0.0453 + 0.00000 i -0.4991 + 0.00000 i Ans (:，:，5) = -1.1077 + 1.0901i -1.1077 - 1.0901i

Bode图进一步突出了从Cr= 8.57Cr=2.

波德(G (clf:,“Tc”), {0.01, 10})

反馈控制策略

为了防止热失控而逐步降低残留浓度，使用反馈控制来调整冷却剂温度Tc基于残留浓度的测量Cr和反应堆温度Tr．对于这个应用，我们使用一个串级控制架构，其中内环调节反应堆温度，而外环跟踪浓度设置。两个反馈回路都是数字化的，采样周期为0.5分钟。

开放式系统(“rct_CSTR”）

目标浓度Cref从t=10时的8.57 kmol/m^3下降到t=36时的2 kmol/m^3（过渡持续26分钟）。相应的剖面图特雷夫反应器的温度是通过对平衡值的插值得到的TrEQ从分析。控制器计算冷却剂温度调节直接转矩相对于初始平衡值TcEQ（1）= 297.98Cr= 8.57。注意，模型是这样设置的，最初，输出TrSP“浓度控制器”块的温度与反应器温度匹配，调整直接转矩温度为零，冷却液温度为零Tc是否处于均衡值TcEQ（1）．

CLF t = [0 10:36 45];C = interp1([0 10 36 45]，[8.57 8.57 2 2]，t);子图(211)，图(t,C)，网格，集合(gca，“ylim”，[0 10]）标题(“目标剩余浓度”), ylabel (“Cref”) subplot(212)， plot(t,interp1(CrEQ,TrEQ,C)) title(“平衡时相应的反应堆温度”), ylabel (“不洁净的”)、网格

控制目标

使用TuningGoal对象来捕获设计需求。首先,Cr应该遵循定位点Cref响应时间约为5分钟。

R1 = TuningGoal。跟踪(“Cref”,“Cr”5);

内环(温度)应该以足够的阻尼和足够快的衰减来稳定反应动力学。

MinDecay = 0.2;MinDamping = 0.5;%约束内环的闭环极点，使外环打开R2 = TuningGoal。波兰人(“Tc”，MinDecay，MinDamping）；R2.开口=“TrSP”；

控制器输出的速率限制块指定冷却液温度Tc变化速度不能超过每分钟10度。这是对控制器权限的严重限制，如果忽略它，可能会导致性能低下或不稳定。要考虑这个利率限制，请注意这一点Cref以0.25 kmol/m^3/min的速率变化，以确保Tc变化速度不超过10度/分钟Cref来Tc应小于10/0.25=40。

R3 = TuningGoal。获得(“Cref”,“Tc”, 40);

最后，在设备输入端要求至少7dB的增益裕度和45度的相位裕度Tc．

R4 = TuningGoal。利润(“Tc”7、45);

增益调度控制器

为了实现这些要求，我们在外环采用PI控制器，在内环采用超前补偿器。由于取样速率缓慢，需要铅补偿器来充分稳定化学反应在中等浓度Cr= 5.28 kmol / m ^ 3 /分钟。因为反应动力学在很大程度上随浓度而变化，我们进一步将控制器增益设定为浓度的函数。这是在Simulink中使用查找表块建金宝app模的，如图1和2所示。

图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，可以自动调整系数 $K_{p0}， K_{p1}， K_{p2}， K_{i0}， \ldots$ 符合要求R1-R4在上面计算的五个平衡点。这相当于在沿线路的五个设计点调整增益计划控制器Cref轨迹，使用tunableSurface对象将每个增益参数化为的二次函数Cr。将“调谐网格”设置为五种浓度CrEQ二次参数化的基函数为 C_r美元,C_r ^ 2美元 .大多数增益初始化为相同的零。

TuningGrid =结构(“Cr”, CrEQ);ShapeFcn = @(Cr) [Cr, Cr^2];Kp = tunableSurface (“金伯利进程”， 0, TuningGrid, ShapeFcn);Ki = tunableSurface (“吻”，-2，调整网格，ShapeFcn）；Kt=可调表面(“Kt”， 0, TuningGrid, ShapeFcn);一个= tunableSurface (“一个”， 0, TuningGrid, ShapeFcn);b = tunableSurface (“b”， 0, TuningGrid, ShapeFcn);

控制器调优

由于目标带宽在奈奎斯特频率的10年内，因此更容易在离散域直接调整控制器。将线性化的过程动态离散，采样时间为0.5分钟。使用ZOH方法来反映数字控制器如何与连续时间工厂相互作用。

t = 0.5;Gd =汇集(G, Ts);

创建一个slTuner用于调整上面介绍的二次增益计划的接口。采用块代换法将非线性对象模型替换为5个离散线性模型Gd在设计点获得CrEQ．使用setBlockParam关联可调增益函数Kp,Ki,Kt,一个,b使用同名的查找表块。

BlockSubs =结构(“名字”,“rct_CSTR /装运箱”,“价值”、Gd);ST0 = slTuner (“rct_CSTR”, {“金伯利进程”,“吻”,“Kt”,“一个”,“b”}, BlockSubs);ST0。Ts = Ts;用于调优的样本时间%登记兴趣点ST0.addPoint ({“Cref”,“Cr”,“Tr”,“TrSP”,“Tc”})参数化查找表块ST0.setBlockParam (“金伯利进程”Kp);ST0.setBlockParam (“吻”、Ki);ST0.setBlockParam (“Kt”、Kt);ST0.setBlockParam (“一个”,);ST0.setBlockParam (“b”, b);

您现在可以使用systune根据要求调整控制器系数R1-R4.将稳定裕度要求作为硬约束条件，并优化其余要求。

ST = systune(ST0，[R1 R2 R3]，R4);

Final: Soft = 1.23, Hard = 0.99991, Iterations = 195

结果设计满足硬约束(硬< 1)，几乎满足其余的要求(软为了验证这一设计，模拟浓度相同的斜坡的响应Cref. 每个图显示了五个设计点的线性响应CrEQ．

t=0:Ts:20；uC=interp1（[02520]，（-0.25）*[033]，t）；子批次（211），lsim（getIOTransfer）（ST，“Cref”,“Cr”)，uC）网格，集合（gca，“ylim”[-1.5 - 0.5])、标题(的残余浓度)次要情节(212),lsim (getIOTransfer(圣,“Cref”,“Tc”)、加州大学)网格,标题(冷却剂温度变化的）

请注意，冷却液温度的变化率保持在物理极限范围内（每分钟10度或每个采样周期5度）。

控制器验证

检查每个增益随时间的变化Cr在过渡期间。

%访问调整后的增益计划TGS=getBlockParam（ST）；绘制增益曲线clf subplot(321)， viewSurf(TGS.Kp)， ylabel(“金伯利进程”)子地块（322）、viewSurf（TGS.Ki）、ylabel(“吻”)子地块（323）、viewSurf（TGS.Kt）、ylabel(“Kt”) subplot(324)， viewSurf(TGS.a)， ylabel(“一个”) subplot(325)， viewSurf(TGS.b)， ylabel(“b”）

为了在Simulink中验证增益调度控制器，首先使用金宝appwriteBlockValue将调优结果应用到Simulink模型中。金宝app对于每个Lookup Table块，将在表断点处计算相应的二次增益公式，并相应地更新表数据。

writeBlockValue（ST）

接下来，按下Play按钮，用调整后的增益计划模拟响应。仿真结果如图3所示。增益调度控制器成功地驱动反应通过过渡充分的响应时间和没有饱和的速率限制(控制器输出杜匹配有效温度变化直接转矩)．反应堆温度保持在它的平衡值附近特雷夫，表示控制器保持反应接近平衡，同时防止热失控。

图3：增益调度级联控制器的转换。

MATLAB中的控制器整定

或者，您可以在MATLAB中直接调整增益调度而不使用slTuner接口。首先将增益参数化为Cr正如上面做的。

TuningGrid =结构(“Cr”, CrEQ);ShapeFcn = @(Cr) [Cr, Cr^2];Kp = tunableSurface (“金伯利进程”， 0, TuningGrid, ShapeFcn);Ki = tunableSurface (“吻”，-2，调整网格，ShapeFcn）；Kt=可调表面(“Kt”， 0, TuningGrid, ShapeFcn);一个= tunableSurface (“一个”， 0, TuningGrid, ShapeFcn);b = tunableSurface (“b”， 0, TuningGrid, ShapeFcn);

使用这些增益来构建PI和先导控制器。

π= pid (Kp,吻,“t”Ts,“时间单位”,“敏”）;π。u =ECr的；π。y =“TrSP”；LEAD = Kt * tf([1 -a]，[1 -b]，Ts，“时间单位”,“敏”）;领先。u =“ETr”；领先。y =“Tc”；

使用连接在五个设计点上建立整体控制系统的闭环模型。标记控制器输出TrSP和Tc作为“分析点”，以便可以打开循环，并在这些位置评估稳定性裕度。闭环模型T0线性模型的5乘1数组是否取决于可调系数Kp, Ki, Kt, a, b．每个模型都是离散的，每半分钟采样一次。

时间单位=“敏”；S1 = sumblk (‘ECr=Cref-Cr’); S2=sumblk('ETr = TrSP - Tr'）;T0 =连接(Gd (:,“Tc”)、铅、π,S1, S2,“Cref”,“Cr”, {“TrSP”,“Tc”}）;

最后，使用systune调整增益计划系数。

T = systune(T0，[R1 R2 R3]，R4);

Final: Soft = 1.23, Hard = 0.99987, Iterations = 193

结果与上述结果相似。通过绘制增益的函数来确定Cr在中使用调谐系数T．

clf subplot(321)， viewSurf(setBlockValue(Kp,T))， ylabel(“金伯利进程”) subplot(322)， viewSurf(setBlockValue(Ki,T))， ylabel(“吻”)子地块（323）、viewSurf（setBlockValue（Kt，T））、ylabel(“Kt”) subplot(324)， viewSurf(setBlockValue(a,T))， ylabel(“一个”) subplot(325)， viewSurf(setBlockValue(b,T))， ylabel(“b”）