自适应控制方法在连续搅拌反应器（CSTR)放热反应过程的应用仿真验证

月生言己

已于 2024-03-12 11:14:13 修改

阅读量25

点赞数

分类专栏：论文仿真文章标签：其他

于 2020-11-29 19:44:40 首次发布

本文链接：https://blog.csdn.net/cabbage_hhh/article/details/110292978

版权

论文仿真专栏收录该内容

1 篇文章 0 订阅

订阅专栏

前言

本文为《一种构造化工过程被控变量的方法》论文的仿真实验记录，供日后温习回顾，也希望能让大家对神经网络结合自适应控制的方法有所了解。我在文中跳过了一些概念，想了解的，建议直接找这篇论文。

一、自适应控制概念理解

自适应控制的目标是找到这样一些变量，它们在不确定扰动下最优点的值保持相对不变，扰动发生时不需要重新进行优化计算，而是在反馈控制器作用下将这些变量维持在原先的设定点，也能使得过程运行在最优点附近。

NCO条件是维持系统运行在最优点的必要条件。NCO条件包括积极约束和简约梯度，保持这两部分为0，可趋使过程接近最优点。
在满足某些约束的前提下最小化成本或是最大化经济效益，可用
$min_uJ(u,d,y)$
$s . t . G (u, d, y) < = 0$
来表示，其中J为目标函数，u,d,y分别是操纵变量、扰动变量和测量变量；G是约束条件。由于稳态时，y是u和d的函数，（测量变量在稳态时只和操纵变量和干扰量相关），所以式子又可以写成：
$min_uJ(u,d)$
$s . t . G (u, d) < = 0$
上式是一个带约束的多变量优化问题，假设u=u_opt应该满足KKT条件
在这里插入图片描述
由上式可得积极约束和简约梯度式子

如果能在过程运行时保证这两个等式成立或者尽量成立"就能最大限度逼近最优点"提高系统经济益。
观察上式可知，等式右边为常数0，因此可以选择Ga和简约梯度为被控变量，并将它们设定点都设为0。
在这里插入图片描述
一般地，工程中积极约束是涉及生产能力)安全生产或者产品质量的约束，较容易获取。而对于简约梯度，其形式含有微分部分，难以直接测量，并且依赖不可测扰动d，因此很难获得。所以，如何在未知d的情况下，获得简约梯度的值成了问题的关键。
这篇论文的重点在于不是像过往方法以测量变量的线性组合来近似简约梯度值，而是通过采集很多干扰情况下测量变量y与简约梯度的非线性关系，用神经网络的方法找出非线性函数Z(y)。

二、实例研究

某连续搅拌反应器（CSTR）放热反应过程，系统的过程模型为
在这里插入图片描述
其反应速率r为
其中停留时间为60s，操纵变量u、测量变量y和扰动变量d如下式描述
目标函数定义为：J=-2.009CB+(0.001657Ti)^2
各变量物理含义及名义工作点见下表

扰动变量可能的波动范围为0.7<=CAi<=1.3, 0<=CBi<=0.3,优化的操作目的为在CAi和CBi波动的情况下最小化J。存在最优的Ti，使得经济效益最大化。

此例的自由度为1，故不存在积极约束Ga,因此，简约梯度可简化为dJ/dTi,满足dJ/dTi=0即可。
故在0.7<=CAi<=1.3, 0<=CBi<=0.3, Ti-424.292<=10%的范围内各自将其十等分，分别计算这 $11^3$ 种情况的稳态y和dJ/dTi值，其中dJ/dTi通过摄动法获得。

在得到系统微分方程和各个变量标称值时，可以将标称值代入微分方程计算，并结合物理意义加以思考。这里的化学过程是放热反应，所以出口T一定会比入口Ti要大，而反应速率一定是正值，要想系统存在稳态，即等式左端 $d T / d t$ 可为0，那么后面一定是符号错误，所以可以大胆判断论文的
在这里插入图片描述
这一式子中，后面的-5r应该为+5r，（后来经证实，确实是论文发刊的时候没改过来），所以在看微分方程时不能仅仅看式子，而是要结合物理意义，进行思考。
想解出系统微分方程的稳态，可以将等式左端直接置0，然后用fsolve求解。

三、仿真过程

为评估控制效果，定义稳态损失函数L
$L=J(u,d)-J_{opt}(d)$
式中 $J (u, d)$ 为实际采用控制系统闭环后的稳态目标函数值， $J_{opt}(d)$ 为真正最优点下的目标函数值。L值越小表明越接近最优点。表2计算比较了6种不同的控制结构在100组随机扰动下的损失函数之和。
在这里插入图片描述

3.1最优值计算

稳态为控制变量为满足系统微分方程而稳定不变的情况，最优值计算则是在操纵变量可变的情况求L的最小值，求稳态情况一般用fsolve，求最优值一般用fmincon。
计算在100种干扰量下的最优值

d=load('data_system_T.txt');
J_raw=load('J_raw_T.txt');
data=zeros(100,2);
x=zeros(100,5);  f=zeros(100,1);  flag=zeros(100,1);
L=zeros(100,1);
L_sum=0;
for i=1:100
   data(i,1)=d(i,1);
   data(i,2)=d(i,2);
end
for i=1:100
    [x(i,:),f(i,1),flag(i,1)]=rctOpt(data(i,:));
    L(i,1)=J_raw(i,1)-f(i,1);
    L_sum=L_sum+L(i,1);
end
L_sum
function [x,f,flag]=rctOpt(d)
options=optimset('LargeScale','off','Display','iter','TolX',1e-6,'TolFun',1e-6,'MaxFunEvals',1e6,'MaxIter',1e6,'algorithm','sqp');
x0=[0.498;0.502;426.803;424;0];
flag=0;
    while ~flag
        [x,f,flag] = fmincon(@(x,d)(-2.009*x(2)+(0.001657*x(4))^2),x0,[],[],[],[],[],[],@rctMdl,options,d);
        x0=x;
    end
end

function [c,ceq,gc,gceq]=rctMdl(x,d)
C_A=x(1);C_B=x(2);T=x(3);T_i=x(4);r=x(5);
C_Ai=d(1);C_Bi=d(2);
tao=60;
eq1=(C_Ai-C_A)/tao-r;%material balance
eq2=(C_Bi-C_B)/tao+r;%material balance
eq3=(T_i-T)/tao-5*r;%energe balance
eq4=5000*exp(-10000/(1.987*T))*C_A-1000000*exp(-15000/(1.987*T))*C_B-r;
c=[];ceq=[eq1;eq2;eq3;eq4];
if nargout>2
    gc=[];
    gceq=[];
end
end

3.2稳态目标函数值计算

解微分方程稳态时，由于是稳态所以可以确定被控变量为设定值常值，比如把CA当作被控变量时就把它定义为常量0.498，把CB当成被控变量，就=0.502来解。但六个控制系统，求最优解时，这个设定值就可以不管了，而且是六个控制系统共用一个最优值，虽然说少了一个设定值的条件，但是求最优值用的是fmincon，有目标函数最小的附带限制条件，所以依然可以求解。

3.2.1 $T_i$ 为被控变量，设定值为424.292

求解100种干扰量下的采用控制系统闭环后的稳态目标函数值。
$T_i$ 设定值为424.292。因此在fsolve求解中，将Ti直接赋值424.292。

set_var=zeros(100,2);

for k=1:100
    set_var(k,1)=0.7+0.6*rand;  %CAI  CBI  TI
    set_var(k,2)=0.3*rand;
    %set_var(k,3)=424.192+0.2*rand;
end
tao=60;  %停留时间
J_raw=zeros(100,1);
fun=@fun_F;
X=zeros(100,4);
for i=1:100
        X(i,:)=fsolve(@(X)fun(X,set_var(i,1),set_var(i,2)),[0.498;0.502;426.803;0.008]);  % X  为CB  T  y=CA CB T Ti
        J_raw(i)=-2.009*X(i,2)+(0.001657*424.292)^2;
        dlmwrite('data_system_T.txt',set_var,'delimiter', '\t','precision',6) 
end
J_raw
dlmwrite('J_raw_T.txt',J_raw,'delimiter', '\t','precision',6) 
function F= fun_F (X,CAi,CBi)
tao=60;
C_A=X(1);C_B=X(2);T=X(3);Ti=424.292;r=X(4);
F(1)= ( CAi -C_A)/tao-r ;  
F(2)= (CBi -C_B)/tao+r ;
F(3)= (Ti -T )/tao+5*r;
F(4)=5000*exp(-10000/(1.987*T))*C_A-1000000*exp(-15000/(1.987*T))*C_B-r;
end

L_sum = 0.5181

3.2.2 $T$ 为被控变量，设定值为426.803

set_var=zeros(100,2);
for k=1:100
    set_var(k,1)=0.7+0.6*rand;  %CAI  CBI  TI
    set_var(k,2)=0.3*rand;
    %set_var(k,3)=424.192+0.2*rand;
end
tao=60;  %停留时间
J_raw=zeros(100,1);
fun=@fun_F;
X=zeros(100,4);
for i=1:100
        X(i,:)=fsolve(@(X)fun(X,set_var(i,1),set_var(i,2)),[0.498;0.502;424.292;0.008]);  % X  为CB  T  y=CA CB T Ti
        J_raw(i)=-2.009*X(i,2)+(0.001657*X(i,3))^2;
        dlmwrite('disturb_T.txt',set_var,'delimiter', '\t','precision',6) 
end
J_raw
dlmwrite('J_raw_T.txt',J_raw,'delimiter', '\t','precision',6) 

function F= fun_F (X,CAi,CBi)
tao=60;
%以后打开MATLAB，先把上一次打开的文件全部关闭，确保编辑的文件都是在我所在的那个文件里面。
%否则，编辑其他文件的同名文件，最后调用的并不是它
C_A=X(1);C_B=X(2);T=426.803;Ti=X(3);r=X(4);
F(1)= ( CAi -C_A)/tao-r ;  
F(2)= (CBi -C_B)/tao+r ;
F(3)= (Ti -T )/tao+5*r;
F(4)=5000*exp(-10000/(1.987*T))*C_A-1000000*exp(-15000/(1.987*T))*C_B-r;
end

L_sum =0.5525

3.2.3 $C_A$ 为被控变量，设定值为0.498

%CA为被控变量
set_var=zeros(100,2);
for k=1:100
    set_var(k,1)=0.7+0.6*rand;  %CAI  CBI  TI
    set_var(k,2)=0.3*rand;
end
tao=60;  %停留时间
J_raw=zeros(100,1);
fun=@fun_F;
X=zeros(100,4);
for i=1:100
        X(i,:)=fsolve(@(X)fun(X,set_var(i,1),set_var(i,2)),[0.502;426.803;424.292;0]);  % X 为CA CB Ti  
        J_raw(i)=-2.009*X(i,1)+(0.001657*X(i,3))^2;
        dlmwrite('disturb_CA.txt',set_var,'delimiter', '\t','precision',6) 
end
dlmwrite('J_raw_CA.txt',J_raw,'delimiter', '\t','precision',6) 

function F= fun_F (X,CAi,CBi)
tao=60;
% T 为被控变量
C_A=0.498;C_B=X(1);T=X(2);T_i=X(3);r=X(4);
F(1)= ( CAi -C_A)/tao-r ;   
F(2)= (CBi -C_B)/tao+r ;
F(3)= (T_i -T )/tao+5*r;
F(4)=5000*exp(-10000/(1.987*T))*C_A-1000000*exp(-15000/(1.987*T))*C_B-r;
end

L_sum = 20.0575，这个地方算的和原论文差别是最大的。原论文是6.158，

3.2.4 $C_B$ 为被控变量，设定值为0.502

set_var=zeros(100,2);
for k=1:100
    set_var(k,1)=0.7+0.6*rand;  %CAI  CBI  TI
    set_var(k,2)=0.3*rand;
    %set_var(k,3)=424.192+0.2*rand;
end
tao=60;  %停留时间
J_raw=zeros(100,1);
fun=@fun_F;
X=zeros(100,4);
for i=1:100
        X(i,:)=fsolve(@(X)fun(X,set_var(i,1),set_var(i,2)),[0.498;426.803;424.292;0.008]);  % X  为CB  T  y=CA CB T Ti
        J_raw(i)=-2.009*0.502+(0.001657*X(i,3))^2;
        dlmwrite('disturb_CB.txt',set_var,'delimiter', '\t','precision',6) 
end
J_raw
dlmwrite('J_raw_CB.txt',J_raw,'delimiter', '\t','precision',6) 

function F= fun_F (X,CAi,CBi)
tao=60;
C_A=X(1);C_B=0.502;T=X(2);Ti=X(3);r=X(4);
F(1)= ( CAi -C_A)/tao-r ;  
F(2)= (CBi -C_B)/tao+r ;
F(3)= (Ti -T )/tao+5*r;
F(4)=5000*exp(-10000/(1.987*T))*C_A-1000000*exp(-15000/(1.987*T))*C_B-r;
end

L_sum = 14.6505，与原文数据相差还在可接受范围内。

3.2.5 $0.761515C_A+0.648127C_B+0.000113T+0.005187Ti$ 为被控变量，设定值为2.1956

set_var=zeros(100,2);
for k=1:100
    set_var(k,1)=0.7+0.6*rand;  %CAI  CBI  TI
    set_var(k,2)=0.3*rand;
    %set_var(k,3)=424.192+0.2*rand;
end
tao=60;  %停留时间
J_raw=zeros(100,1);
fun=@fun_F;
X=zeros(100,5);
for i=1:100
        X(i,:)=fsolve(@(X)fun(X,set_var(i,1),set_var(i,2)),[0.498;0.502;426.803;424.292;0.008]);  % X  为CB  T  y=CA CB T Ti
        J_raw(i)=-2.009*X(i,2)+(0.001657*X(i,4))^2;
        dlmwrite('disturb_linear.txt',set_var,'delimiter', '\t','precision',6) 
end
J_raw
dlmwrite('J_raw_linear.txt',J_raw,'delimiter', '\t','precision',6) 

function F= fun_F (X,CAi,CBi)
tao=60;
C_A=X(1);C_B=X(2);T=X(3);Ti=X(4);r=X(5);
F(1)= ( CAi -C_A)/tao-r ;  
F(2)= (CBi -C_B)/tao+r ;
F(3)= (Ti -T )/tao+5*r;
F(4)=5000*exp(-10000/(1.987*T))*C_A-1000000*exp(-15000/(1.987*T))*C_B-r;
F(5)=-0.761515*C_A+0.648127*C_B+0.000113*T+0.005187*Ti-2.1956;
end

L_sum =0.0159，这里已经比上面的任意单一被控变量要好很多了，可见取测量变量的线性组合作为被控变量确实有可取之处。下面来介绍本文核心方法。

3.2.6 神经网络得出的 $Z (y)$ 为被控变量，设定值为0

$Z (y)$ 为以 $y$ 为自变量， $dJ/dT_i$ 为因变量的非线性函数。在0.7<=CAi<=1.3, 0<=CBi<=0.3, Ti-424.292<=10%的范围内各自将其十等分，分别计算这 $11^3$ 种情况的稳态y和dJ/dTi值，其中dJ/dTi通过摄动法获得。然后用神经网络工具箱得到 $Z (y)$ 。

clc;%计算11*11*11种干扰情况下，的稳态y和dJ/dTi值，得到以y为自变量，dJ/dTi为因变量的非线性函数Z(y)
clear;
tao=60;  %停留时间
CAi=linspace(0.7,1.3,11);
CBi=linspace(0,0.3,11);
Ti=linspace(381.8628,466.7212,11);
y=zeros(1331,4);index=0;
DJ=zeros(1331,1);
for i=1:11
   for j=1:11
       for k=1:11
       index=index+1;
       options=optimset('Display','iter','TolX',1e-10,'TolFun',1e-10,'MaxFunEvals',1e6,'MaxIter',1e6);
       flag=0;flag1=0;
       while ~flag
          [X,fval1,flag]=fsolve(@(X)fun_ZY(X,CAi(i),CBi(j),Ti(k)),[0.498;0.502;426.803],options);  
       end
        y(index,1)=X(1);   y(index,2)=X(2);     y(index,3)=X(3);      y(index,4)=Ti(k);
        J=-2.009*X(2)+(0.001657*Ti(k))^2;
        while ~flag1
       [X,fval2,flag1]=fsolve(@(X)fun_ZY(X,CAi(i),CBi(j),Ti(k)+0.01),[0.498;0.502;426.803],options); 
        end
       J_new=-2.009*X(2)+(0.001657*(Ti(k)+0.01))^2;
       DJ(index,1)=(J_new-J)/0.01;
       end
   end
end
dlmwrite('Zy_y.txt',y,'delimiter', '\t','precision',6) 
dlmwrite('Zy_DJ.txt',DJ,'delimiter', '\t','precision',6) 

function F= fun_ZY (X,CAi,CBi,Ti)
tao=60;
C_A=X(1);C_B=X(2);T=X(3);
r=5000*exp(-10000/(1.987*T))*C_A-1000000*exp(-15000/(1.987*T))*C_B;
F(1)= ( CAi -C_A)/tao-r ;  
F(2)= (CBi -C_B)/tao+r ;
F(3)= (Ti -T )/tao+5*r;
end

算出DJ, y, 使用神经网络工具箱求出非线性函数 $Z (y)$ 作为求解微分方程的约束条件之一。如果调试出了什么问题，不要过分去怀疑神经网络的训练集、测试集、验证集比例，而应该优先怀疑数据是否正确，数据好了，超参数即使没选好，拟合的效果也不会太差。
类似(T-424.292)<=10% 这种式子，下次就不要认为是Ti<=424.392了，，，原数据直接错了，这还神经网络个鬼。

clc %计算100组干扰下，满足微分方程的y=[X(i,1),X(i,2),X(i,3),X(i,4)]'，进而计算J_raw
clear
set_var=zeros(100,2);
for k=1:100
    set_var(k,1)=0.7+0.6*rand;  %CAI  CBI  TI
    set_var(k,2)=0.3*rand;
end

tao=60;  %停留时间
J_raw=zeros(100,1);
X=zeros(100,5);
BPOUTPUT=zeros(100,1);
fval_error_cv=zeros(100,1);
for i=1:100
        flag_cv=0;
        options=optimset('Display','iter','TolX',1e-6,'TolFun',1e-6,'MaxFunEvals',1e6,'MaxIter',1e6);
        while ~flag_cv
        [X(i,:),fval,flag_cv]=fsolve(@(X)fun_F(X,set_var(i,1),set_var(i,2)),[0.498;0.502;426.803;424.292;0.0057],options);  % 解微分方程的初值也是列矩阵，要么分号，要么加’,
     %  另外，只要求微分方程，那解出的大致就是达到最优值的y_opt接近的点，初值要么设标称点，要么设
        end
        fval_error_cv(i,1)=sum(sum(fval.*fval));
        y=[X(i,1),X(i,2),X(i,3),X(i,4)]';  %千万别忘了 输入数据的矩阵要转置
        J_raw(i,1)=-2.009*X(i,2)+(0.001657*X(i,4))^2;
        dlmwrite('disturb_ZY.txt',set_var,'delimiter', '\t','precision',6) 
end
xlswrite('fval_zy.xls', fval_error_cv)
dlmwrite('J_raw_ZY.txt',J_raw,'delimiter', '\t','precision',6) 

function F= fun_F (X,CAi,CBi)
tao=60;
C_A=X(1);C_B=X(2);T=X(3);Ti=X(4);r=X(5);
F(1)= ( CAi -C_A)/tao-r ;  
F(2)= (CBi -C_B)/tao+r ;
F(3)= (Ti -T )/tao+5*r;
F(4)=5000*exp(-10000/(1.987*T))*C_A-1e6*exp(-15000/(1.987*T))*C_B-r;
y=[C_A, C_B, T ,Ti];  %千万别忘了 输入数据的矩阵要转置
DJ=myNeuralNetworkFunction(y);
F(5)=DJ;
end

100组随机扰动下的损失函数之和 L_sum =0.0088
虽然比论文中还大了一个数量级，不过证明这个方法已经很优秀了。

四、实验结果

No.	CV	sum(L)	论文原文
1	CA	20.0575	6.158
2	CB	14.6505	20.582
3	T	0.5525	0.545
4	Ti	0.5181	0.499
5	-0.761515C_A+0.648127C_B+0.000113T+0.005187Ti	0.0159	0.0158
6	Z(y)	0.0088	0.00043

MATLAB备忘录：

以后打开MATLAB，先把上一次打开的文件全部关闭，确保编辑的文件都是在我所在的那个文件里面。否则，编辑其他文件的同名文件，最后调用的并不是它

月生言己

关注

0
点赞
踩
0

收藏

觉得还不错? 一键收藏
0
评论
自适应控制方法在连续搅拌反应器（CSTR)放热反应过程的应用仿真验证

前言本文为《一种构造化工过程被控变量的方法》论文的仿真实验记录，供日后温习回顾，也希望能让大家对神经网络结合自适应控制的方法有所了解。我在文中跳过了一些概念，想了解的，建议直接找这篇论文。一、自适应控制概念理解自适应控制的目标是找到这样一些变量，它们在不确定扰动下最优点的值保持相对不变，扰动发生时不需要重新进行优化计算，而是在反馈控制器作用下将这些变量维持在原先的设定点，也能使得过程运行在最优点附近。NCO条件是维持系统运行在最优点的必要条件。NCO条件包括积极约束和简约梯度，保持这两部分为0，可趋
复制链接

扫一扫