CN106368772A - Scr***尿素喷射控制方法 - Google Patents

Abstract

一种SCR***尿素喷射控制方法，属于柴油机控制技术领域。本发明的目的是从尿素SCR***化学反应机理出发，研究***的偏微分方程模型（PDE），选择合适的逼近方法简化模型，利用滚动优化控制处理时滞和不确定性优势的SCR***尿素喷射控制方法。本发明的步骤：a、尿素SCR***偏微分建模，b、偏微分***预测控制器设计。本发明基于实际的尿素SCR***控制需求，将分布参数***建模与预测控制器设计结合在一起，能够在精确描述***动态性能的同时，解决该***的约束条件下的非线性优化问题。

Description

SCR***尿素喷射控制方法

技术领域

本发明属于柴油机控制技术领域。

背景技术

与汽油机相比，柴油机具有更高的燃油经济性和动力输出。然而，由于其稀燃的特性，柴油机会比汽油机产生更多的NO_x有害气体。随着全世界范围内针对NO_x排放越来越严格的法规出台，多种为降低NO_x的排放后处理***面世了。这些技术就包括尿素选择性催化还原(SCR)***。尿素SCR***工作时则不需额外燃油，而且尿素消耗也相对较低，凭借这些优势，已经在汽车工业界占据一定优势。在我国，目前的实际国情是燃油中硫含量较高，而且许多种排放控制技术推广都受到限制。所以，凭借其对硫的敏感性较低的特性，尿素SCR排放后处理技术在我国的发展更具优势。

尿素SCR技术的基本原理是利用NOx与氨(NH3)之间的氧化还原反应，而所用的氨一般都来源于32.5％的尿素溶液(添蓝溶液)。虽然氨能够还原NOx，但其较高的排放也是对人体有害的，并且有着刺鼻的气味。为实现较高的NOx转化效率，要有充分的氨做为还原剂；但是，这一点反过来会增加氨的逃逸量，这一矛盾成为了尿素SCR***研究面临的主要挑战之一。目前更为普遍的共识是，通过改进尿素喷射控制技术达到上述目标，是一种较便捷且经济的方法。针对尿素SCR控制***的研究，当前较为先进的是一些基于模型的反馈、前馈控制方法。但是，这些模型大多是基于常微分方程(ODE)建立的。SCR***作为一个典型的多物质沿催化器轴向不均匀传播、并发生复杂化学反应的分布参数***，不但具有无限维特性，复杂的多物质轴向传播的环境也使***具有较大的滞后和不确定性，上述基于集中参数的ODE模型的控制方法很难达到理想的控制效果。此外，对于尿素SCR***参数来说，也存在较多的约束条件。例如，***输入(尿素喷射器喷射量)有最大值限制，***输出(NO_x与NH₃)受到排放法规的限制等。这些带约束的优化控制问题，适合于在模型预测控制的框架内解决。

发明内容

本发明的目的是从尿素SCR***化学反应机理出发，研究***的偏微分方程模型(PDE)，选择合适的逼近方法简化模型，利用滚动优化控制处理时滞和不确定性优势的SCR***尿素喷射控制方法。

本发明的步骤：

a、尿素SCR***偏微分建模

基于尿素SCR***各主要物质沿轴向传输特性，主要参考了动态传输反应模型：

$\frac{\∂C_k^g}{\∂t}=-\frac{\∂J_k^g}{\∂z}-K_m^k(C_k^g-C_k^w),---\left(1\right)$

磁通密度公式：

J_k＝θ_zC_k， (2)

以及一种物质浓度在催化器表面的动态及稳态假设得到的：

$-K_m^k(C_k^g-C_k^w)==\frac{(n_{k,in}^g-n_{k,out}^g)}{{\mathrm{\ϵV}}_c}+r_j,---\left(3\right)$

建立得到***的一阶双曲型PDE模型：

$\frac{\∂C_{{\mathrm{NO}}_x}^g}{\∂t}=-{\mathrm{\θ}}_z\frac{\∂C_{{\mathrm{NO}}_x}^g}{\∂z}+\frac{(n_{{\mathrm{NO}}_x,in}^g-n_{{\mathrm{NO}}_x,out}^g)}{{\mathrm{\ϵV}}_c}-r_{SCR},$

$\frac{{\mathrm{d\Θ}}_{{\mathrm{NO}}_x}}{dt}=r_{ads}-r_{des}-r_{SCR}-r_{ox},$

$\frac{\∂C_{{\mathrm{NH}}_3}^g}{\∂t}=-{\mathrm{\θ}}_z\frac{\∂C_{{\mathrm{NH}}_3}^g}{\∂z}+\frac{(n_{{\mathrm{NH}}_3,in}^g-n_{{\mathrm{NH}}_3,out}^g)}{{\mathrm{\ϵV}}_c}-r_{ads}+r_{des}.---\left(4\right)$

其中，

$r_{ads}=c_s{S}_c{\mathrm{\α}}_{prob}\sqrt{\frac{RT}{2{\mathrm{\πM}}_{{\mathrm{NH}}_3}}}{C}_{{\mathrm{NH}}_3}(1-{\mathrm{\Θ}}_{{\mathrm{NH}}_3}),$

$r_{des}=c_s{k}_{des}{e}^{-\frac{E_{a,des}}{RT}}{\mathrm{\Θ}}_{{\mathrm{NH}}_3},$

$r_{SCR}=c_s{\mathrm{RTk}}_{SCR}{e}^{\left(\frac{-E_{a,\mathrm{Re}d}}{RT}\right)}{C}_{{\mathrm{NO}}_x}{\mathrm{\Θ}}_{{\mathrm{NH}}_3},$

$r_{ox}=c_c{k}_{Ox}{e}^{\left(\frac{-E_{a,Ox}}{RT}\right)}{\mathrm{\Θ}}_{{\mathrm{NH}}_3}.---\left(5\right)$

进一步整理得到：

$\frac{\∂C_{{\mathrm{NO}}_x}^g}{\∂t}=-{\mathrm{\θ}}_z\frac{\∂C_{{\mathrm{NO}}_x}^g}{\∂z}+a_1{n}_{{\mathrm{NO}}_x,in}^g-C_{{\mathrm{NO}}_x}^g(a_0{a}_1{m}_{EG}^{*}T+a_4(T\left){\mathrm{\Θ}}_{{\mathrm{NH}}_3}\right),$

$c_s\frac{{\mathrm{d\Θ}}_{{\mathrm{NH}}_3}}{dt}=a_2\left(T\right)(1-{\mathrm{\Θ}}_{{\mathrm{NH}}_3})C_{{\mathrm{NH}}_3}^g-\[a_3\left(T\right)+a_4\left(T\right)C_{{\mathrm{NO}}_x}^g+a_5\left(T\right)\]{\mathrm{\Θ}}_{{\mathrm{NH}}_3},$

$\frac{\∂C_{{\mathrm{NH}}_3}^g}{\∂t}=-{\mathrm{\θ}}_z\frac{\∂C_{{\mathrm{NH}}_3}^g}{\∂z}+a_1{n}_{{\mathrm{NH}}_3,in}^g-C_{{\mathrm{NH}}_3}^g\[a_0{a}_1{m}_{EG}^{*}T+a_2\left(T\right)(1-{\mathrm{\Θ}}_{{\mathrm{NH}}_3})\]+a_3\left(T\right){\mathrm{\Θ}}_{{\mathrm{NH}}_3}.---\left(6\right)$

其中，

$a_0=\frac{R_{S,EG}}{P_{amb}},a_1=\frac{1}{{\mathrm{\ϵV}}_c},$

$a_2\left(T\right)=c_s{S}_c{\mathrm{\α}}_{\mathrm{Pr}ob}\sqrt{\frac{RT}{2{\mathrm{\πM}}_{{\mathrm{NH}}_3}}},$

$a_3\left(T\right)=c_s{k}_{Des}{e}^{\left(\frac{-E_{a,Des}}{RT}\right)},$

$a_4\left(T\right)=c_s{\mathrm{RTk}}_{SCR}{e}^{\left(\frac{-E_{a,SCR}}{RT}\right)},$

$a_5\left(T\right)=c_s{k}_{Ox}{e}^{\left(\frac{-E_{a,Ox}}{RT}\right)},$

$a_6=\frac{c_{p,EG}}{c_{p,c}{m}_c},a_7=\frac{{\mathrm{\ϵ}}_{rad,scr}{\mathrm{\σ}}_{sb}{A}_{rad,scr}}{c_{p,c}{m}_c}.---\left(7\right)$

表1和表2分别显示了模型中所有常量和变量的相关定义及参数名义参考值：

表1常量命名法

表2变量命名法

将T，C_NOx等看做可测量变量，考虑的状态量仅为和定义状态变量时变参数控制输入变量被控输出为得到***面向控制模型：

$\left\{\begin{array}{c}\stackrel{\·}{x_1}=f_{11}(x,p)+f_{12}(x,p)x_2,\\ \stackrel{\·}{x_2}=f_{21}(x,p)+a_{1}u-{\mathrm{\θ}}_z\frac{\∂x_2}{\∂z},\end{array}\right.---\left(8\right)$

其中，

$\left\{\begin{array}{c}{f}_{11}(x,p)=-\frac{1}{c_s}\[a_3\left(T\right)+a_4\left(T\right)C_{{\mathrm{NO}}_x}+a_5\left(T\right)\]x_1,\\ f_{12}(x,p)=\frac{1}{c_s}{a}_2\left(T\right)(1-x_1),\\ f_{21}(x,p)=-x_2\[a_0{a}_1{m}_{EG}^{*}T+a_2\left(T\right)(1-x_1)\]+a_3\left(T\right)x_1,\end{array}\right.---\left(9\right)$

进一步整理成状态空间形式：

$\begin{array}{c}\frac{\∂X(z,t)}{\∂t}=\left[\begin{array}{cc}0& 0\\ 0& -{\mathrm{\θ}}_z\end{array}\right]\frac{\∂X(z,t)}{\∂z}+\left[\begin{array}{c}{f}_{11}(x,p)+f_{12}(x,p)x_2\\ f_{21}(x,p)+a_{1}u\end{array}\right],\\ y(z,t)=\left[\begin{array}{cc}1& 0\end{array}\right]X(z,t),\end{array}---\left(10\right)$

上述一阶双曲型***表示为以下状态空间形式：

$\frac{\∂X(z,t)}{\∂t}=A\frac{\∂X(z,t)}{\∂z}+B\left(X\right(z,t),U(z,t\left)\right)---\left(11\right)$

其中，

$\begin{array}{c}A=\left[\begin{array}{cc}0& 0\\ 0& -{\mathrm{\θ}}_z\end{array}\right],\\ B\left(X\right(z,t),U(z,t\left)\right)=\left[\begin{array}{c}{f}_{11}(x,p)+f_{12}(x,p)x_2\\ f_{21}(x,p)+a_{1}u\end{array}\right],\\ C=\left[\begin{array}{cc}1& 0\end{array}\right].\end{array}---\left(12\right);$

b、偏微分***预测控制器设计：

选取Haar小波尺度函数φ(x)，将时空变量展开如下形式：

$X(z,t)\≈X_j(z,t)=\underset{n=0}{\overset{2^j}{\Σ}}{X}_j(z_n,t)\·{\mathrm{\φ}}_j(z-z_n)---\left(13\right)$

一阶空间偏导数表示为：

$\frac{\∂X(z,t)}{\∂z}\≈\frac{\∂X_j^{\left(1\right)}(z,t)}{\∂z}=\underset{n=0}{\overset{2^j}{\Σ}}{X}_j(z_n,t)\·{\mathrm{\φ}}_j^{\left(1\right)}(z-z_n)---\left(14\right)$

采用Eular前向差分法对时间偏导数进行离散：

$\frac{\∂X(z_n,t)}{\∂t}\≈\frac{\∂X_j^{\left(1\right)}(z_n,t)}{\∂z}\≈\frac{X_j(z_n,i+1)-X_j(z_n,i)}{\mathrm{\Δ}t}---\left(15\right)$

应用上述方法，同时将时空偏导数替换得到：

$\begin{array}{c}\frac{X_j(z_n,i+1|i)-X_j(z_n,i)}{\mathrm{\Δ}t}=A\·\underset{n=0}{\overset{2^j}{\Σ}}{X}_j(z_n,t)\·{\mathrm{\φ}}_j^{\left(1\right)}(z-z_n)\\ +B\left(\underset{n=0}{\overset{2^j}{\Σ}}{X}_j\right(z_n,t)\·{\mathrm{\φ}}_j(z-z_n),U(z,i))\end{array}---\left(16\right)$

对上式进行整理，得到状态变量(z_n,i+1|i)的预测值为：

$\begin{array}{c}{X}_j(z_n,i+1|i)=X_j(z_n,i)+\mathrm{\Δ}t\·(A\·\underset{n=0}{\overset{2^j}{\Σ}}{X}_j(z_n,t)\·{\mathrm{\φ}}_j^{\left(1\right)}(z-z_n)\\ +B\left(\underset{n=0}{\overset{2^j}{\Σ}}{X}_j\right(z_n,t)\·{\mathrm{\φ}}_j(z-z_n),U(z,i))\end{array}---\left(17\right)$

由y＝C·X得P步向前预测输出值：

$\begin{array}{c}{\hat{y}}_j(z_n,i+1|i)=y_j(z_n,i)+\mathrm{\Δ}t\·(A\·\underset{n=0}{\overset{2^j}{\Σ}}{y}_j(z_n,t)\·{\mathrm{\φ}}_j^{\left(1\right)}(z-z_n)\\ +B\left(\underset{n=0}{\overset{2^j}{\Σ}}{y}_j\right(z_n,i)\·{\mathrm{\φ}}_j(z-z_n),U(z,i))\\ \·\\ \·\\ \·\\ {\hat{y}}_j(z_n,i+P+1|i+P-1)={\hat{y}}_j(z_n,i+P-1|i+P-2)+\mathrm{\Δ}t\·\\ (A\·\underset{n=0}{\overset{2^j}{\Σ}}{\hat{y}}_j(z_n,i+P-I|i+P-2)\·{\mathrm{\φ}}_j^{\left(1\right)}(z-z_n)+\\ B\left(\underset{n=0}{\overset{2^j}{\Σ}}{\hat{y}}_j\right(z_n,i)\·{\mathrm{\φ}}_j(z-z_n),U(z,i)\left)\right)\end{array}---\left(18\right)$

为了快速而平滑地得到***设定值，定义如下参考轨迹：

y^*(z_n,k+i|k) (19)

预测控制输出为：

$\widehat{y_c}(z,k+i|k)=\hat{y}(z,k+i|k)+d(z,k|k)---\left(20\right)$

定义***误差函数为：

$e(z,k+i|k)=y^{*}(z,k+i|k)-{\hat{y}}_c(z,k+i|k)---\left(21\right).$

本发明基于实际的尿素SCR***控制需求，将分布参数***建模与预测控制器设计结合在一起，能够在精确描述***动态性能的同时，解决该***的约束条件下的非线性优化问题。其优点：

1、本发明提出的偏微分建模方法能够准确描述SCR***动态。

2、本发明基于实际的尾气排放处理***，将偏微分建模过程与模型预测控制器推导相结合。

3、本发明提出的模型预测控制器，能够处理排放约束的优化控制问题。

附图说明

图1是基于分布式参数建模及预测控制的SCR***尿素喷射控制器设计技术路线图。

具体实施方式

本发明主要技术路线包括：根据尿素SCR***控制需求与控制问题的提炼与描述，建立偏微分模型；参照化学反应机理以及数学逼近方法，推导出偏微分***预测控制器；对优化问题进行归纳，并给出优化问题的求解方式；最后，得到满足排放法规要求的控制优化算法。其的步骤是：

a、尿素SCR***偏微分建模

基于尿素SCR***各主要物质沿轴向传输特性，主要参考了动态传输反应模型：

$\frac{\∂C_k^g}{\∂t}=-\frac{\∂J_k^g}{\∂z}-K_m^k(C_k^g-C_k^w),---\left(1\right)$

磁通密度公式：

J_k＝θ_zC_k， (2)

以及一种物质浓度在催化器表面的动态及稳态假设得到的：

$-K_m^k(C_k^g-C_k^w)==\frac{(n_{k,in}^g-n_{k,out}^g)}{{\mathrm{\ϵV}}_c}+r_j,---\left(3\right)$

建立得到***的一阶双曲型PDE模型：

$\frac{\∂C_{{\mathrm{NO}}_x}^g}{\∂t}=-{\mathrm{\θ}}_z\frac{\∂C_{{\mathrm{NO}}_x}^g}{\∂z}+\frac{(n_{{\mathrm{NO}}_x,in}^g-n_{{\mathrm{NO}}_x,out}^g)}{{\mathrm{\ϵV}}_c}-r_{SCR},$

$\frac{{\mathrm{d\Θ}}_{{\mathrm{NO}}_x}}{dt}=r_{ads}-r_{des}-r_{SCR}-r_{ox},$

$\frac{\∂C_{{\mathrm{NH}}_3}^g}{\∂t}=-{\mathrm{\θ}}_z\frac{\∂C_{{\mathrm{NH}}_3}^g}{\∂z}+\frac{(n_{{\mathrm{NH}}_3,in}^g-n_{{\mathrm{NH}}_3,out}^g)}{{\mathrm{\ϵV}}_c}-r_{ads}+r_{des}.---\left(4\right)$

其中，

$r_{ads}=c_s{S}_c{\mathrm{\α}}_{prob}\sqrt{\frac{RT}{2{\mathrm{\πM}}_{{\mathrm{NH}}_3}}}{C}_{{\mathrm{NH}}_3}(1-{\mathrm{\Θ}}_{{\mathrm{NH}}_3}),$

$r_{des}=c_s{k}_{des}{e}^{-\frac{E_{a,des}}{RT}}{\mathrm{\Θ}}_{{\mathrm{NH}}_3},$

$r_{SCR}=c_s{\mathrm{RTk}}_{SCR}{e}^{\left(\frac{-E_{a,\mathrm{Re}d}}{RT}\right)}{C}_{{\mathrm{NO}}_x}{\mathrm{\Θ}}_{{\mathrm{NH}}_3},$

$r_{ox}=c_c{k}_{Ox}{e}^{\left(\frac{-E_{a,Ox}}{RT}\right)}{\mathrm{\Θ}}_{{\mathrm{NH}}_3}.---\left(5\right)$

进一步整理得到：

$\frac{\∂C_{{\mathrm{NO}}_x}^g}{\∂t}=-{\mathrm{\θ}}_z\frac{\∂C_{{\mathrm{NO}}_x}^g}{\∂z}+a_1{n}_{{\mathrm{NO}}_x,in}^g-C_{{\mathrm{NO}}_x}^g(a_0{a}_1{m}_{EG}^{*}T+a_4(T\left){\mathrm{\Θ}}_{{\mathrm{NH}}_3}\right),$

$c_s\frac{{\mathrm{d\Θ}}_{{\mathrm{NH}}_3}}{dt}=a_2\left(T\right)(1-{\mathrm{\Θ}}_{{\mathrm{NH}}_3})C_{{\mathrm{NH}}_3}^g-\[a_3\left(T\right)+a_4\left(T\right)C_{{\mathrm{NO}}_x}^g+a_5\left(T\right)\]{\mathrm{\Θ}}_{{\mathrm{NH}}_3},$

$\frac{\∂C_{{\mathrm{NH}}_3}^g}{\∂t}=-{\mathrm{\θ}}_z\frac{\∂C_{{\mathrm{NH}}_3}^g}{\∂z}+a_1{n}_{{\mathrm{NH}}_3,in}^g-C_{{\mathrm{NH}}_3}^g\[a_0{a}_1{m}_{EG}^{*}T+a_2\left(T\right)(1-{\mathrm{\Θ}}_{{\mathrm{NH}}_3})\]+a_3\left(T\right){\mathrm{\Θ}}_{{\mathrm{NH}}_3}.---\left(6\right)$

其中，

$a_0=\frac{R_{S,EG}}{P_{amb}},a_1=\frac{1}{{\mathrm{\ϵV}}_c},$

$a_2\left(T\right)=c_s{S}_c{\mathrm{\α}}_{\mathrm{Pr}ob}\sqrt{\frac{RT}{2{\mathrm{\πM}}_{{\mathrm{NH}}_3}}},$

$a_3\left(T\right)=c_s{k}_{Des}{e}^{\left(\frac{-E_{a,Des}}{RT}\right)},$

$a_4\left(T\right)=c_s{\mathrm{RTk}}_{SCR}{e}^{\left(\frac{-E_{a,SCR}}{RT}\right)},$

$a_5\left(T\right)=c_s{k}_{Ox}{e}^{\left(\frac{-E_{a,Ox}}{RT}\right)},$

$a_6=\frac{c_{p,EG}}{c_{p,c}{m}_c},a_7=\frac{{\mathrm{\ϵ}}_{rad,scr}{\mathrm{\σ}}_{sb}{A}_{rad,scr}}{c_{p,c}{m}_c}.---\left(7\right)$

表1和表2分别显示了模型中所有常量和变量的相关定义及参数名义参考值：

表1常量命名法

表2变量命名法

为了实现氨覆盖率控制，将T，C_NOx等看做可测量变量，考虑的状态量仅为和定义状态变量时变参数控制输入变量被控输出为得到***面向控制模型：

$\left\{\begin{array}{c}\stackrel{\·}{x_1}=f_{11}(x,p)+f_{12}(x,p)x_2,\\ \stackrel{\·}{x_2}=f_{21}(x,p)+a_{1}u-{\mathrm{\θ}}_z\frac{\∂x_2}{\∂z},\end{array}\right.---\left(8\right)$

其中，

$\left\{\begin{array}{c}{f}_{11}(x,p)=-\frac{1}{c_s}\[a_3\left(T\right)+a_4\left(T\right)C_{{\mathrm{NO}}_x}+a_5\left(T\right)\]x_1,\\ f_{12}(x,p)=\frac{1}{c_s}{a}_2\left(T\right)(1-x_1),\\ f_{21}(x,p)=-x_2\[a_0{a}_1{m}_{EG}^{*}T+a_2\left(T\right)(1-x_1)\]+a_3\left(T\right)x_1,\end{array}\right.---\left(9\right)$

进一步整理成状态空间形式：

$\begin{array}{c}\frac{\∂X(z,t)}{\∂t}=\left[\begin{array}{cc}0& 0\\ 0& -{\mathrm{\θ}}_z\end{array}\right]\frac{\∂X(z,t)}{\∂z}+\left[\begin{array}{c}{f}_{11}(x,p)+f_{12}(x,p)x_2\\ f_{21}(x,p)+a_{1}u\end{array}\right],\\ y(z,t)=\left[\begin{array}{cc}1& 0\end{array}\right]X(z,t),\end{array}---\left(10\right)$

上述一阶双曲型***表示为以下状态空间形式：

$\frac{\∂X(z,t)}{\∂t}=A\frac{\∂X(z,t)}{\∂z}+B\left(X\right(z,t),U(z,t\left)\right)---\left(11\right)$

其中，

$\begin{array}{c}A=\left[\begin{array}{cc}0& 0\\ 0& -{\mathrm{\θ}}_z\end{array}\right],\\ B\left(X\right(z,t),U(z,t\left)\right)=\left[\begin{array}{c}{f}_{11}(x,p)+f_{12}(x,p)x_2\\ f_{21}(x,p)+a_{1}u\end{array}\right],\\ C=\left[\begin{array}{cc}1& 0\end{array}\right].\end{array}---\left(12\right).$

b、偏微分***预测控制器设计：

应用时空分解理论，选择恰当的空间基函数进行逼近。为了降低计算量，本发明选取简单实用的Haar小波尺度函数φ(x)，将时空变量展开如下形式：

$X(z,t)\≈X_j(z,t)=\underset{n=0}{\overset{2^j}{\Σ}}{X}_j(z_n,t)\·{\mathrm{\φ}}_j(z-z_n)---\left(13\right)$

一阶空间偏导数表示为：

$\frac{\∂X(z,t)}{\∂z}\≈\frac{\∂X_j^{\left(1\right)}(z,t)}{\∂z}=\underset{n=0}{\overset{2^j}{\Σ}}{X}_j(z_n,t)\·{\mathrm{\φ}}_j^{\left(1\right)}(z-z_n)---\left(14\right)$

采用Eular前向差分法对时间偏导数进行离散：

$\frac{\∂X(z_n,t)}{\∂t}\≈\frac{\∂X_j^{\left(1\right)}(z_n,t)}{\∂z}\≈\frac{X_j(z_n,i+1)-X_j(z_n,i)}{\mathrm{\Δ}t}---\left(15\right)$

应用上述方法，同时将时空偏导数替换得到：

$\begin{array}{c}\frac{X_j(z_n,i+1|i)-X_j(z_n,i)}{\mathrm{\Δ}t}=A\·\underset{n=0}{\overset{2^j}{\Σ}}{X}_j(z_n,t)\·{\mathrm{\φ}}_j^{\left(1\right)}(z-z_n)\\ +B\left(\underset{n=0}{\overset{2^j}{\Σ}}{X}_j\right(z_n,t)\·{\mathrm{\φ}}_j(z-z_n),U(z,i))\end{array}---\left(16\right)$

对上式进行整理，得到状态变量(z_n,i+1|i)的预测值为：

$\begin{array}{c}{X}_j(z_n,i+1|i)=X_j(z_n,i)+\mathrm{\Δ}t\·(A\·\underset{n=0}{\overset{2^j}{\Σ}}{X}_j(z_n,t)\·{\mathrm{\φ}}_j^{\left(1\right)}(z-z_n)\\ +B\left(\underset{n=0}{\overset{2^j}{\Σ}}{X}_j\right(z_n,t)\·{\mathrm{\φ}}_j(z-z_n),U(z,i))\end{array}---\left(17\right)$

由y＝C·X得P步向前预测输出值：

$\begin{array}{c}{\hat{y}}_j(z_n,i+1|i)=y_j(z_n,i)+\mathrm{\Δ}t\·(A\·\underset{n=0}{\overset{2^j}{\Σ}}{y}_j(z_n,t)\·{\mathrm{\φ}}_j^{\left(1\right)}(z-z_n)\\ +B\left(\underset{n=0}{\overset{2^j}{\Σ}}{y}_j\right(z_n,i)\·{\mathrm{\φ}}_j(z-z_n),U(z,i))\\ \·\\ \·\\ \·\\ {\hat{y}}_j(z_n,i+P+1|i+P-1)={\hat{y}}_j(z_n,i+P-1|i+P-2)+\mathrm{\Δ}t\·\\ (A\·\underset{n=0}{\overset{2^j}{\Σ}}{\hat{y}}_j(z_n,i+P-I|i+P-2)\·{\mathrm{\φ}}_j^{\left(1\right)}(z-z_n)+\\ B\left(\underset{n=0}{\overset{2^j}{\Σ}}{\hat{y}}_j\right(z_n,i)\·{\mathrm{\φ}}_j(z-z_n),U(z,i)\left)\right)\end{array}---\left(18\right)$

为了快速而平滑地得到***设定值，定义如下参考轨迹：

y^*(z_n,k+i|k) (19)

考虑扰动，预测控制输出为：

$\widehat{y_c}(z,k+i|k)=\hat{y}(z,k+i|k)+d(z,k|k)---\left(20\right)$

定义***误差函数为：

$e(z,k+i|k)=y^{*}(z,k+i|k)-{\hat{y}}_c(z,k+i|k)---\left(21\right).$

上述预测控制求解，是典型的约束非线性二次规划问题，可以采用粒子群优化方法求解。

随着分布参数***理论与预测控制理论发展，已经有学者将两者相结合，并提出了一些理论方法，例如：基于特征线法的鲁棒预测控制算法、基于小波正交逼近的分布参数***预测控制方法等。有学者已经将一些能够处理真实分布参数***约束问题的预测控制算法，应用到如薄膜沉积过程、化工生产过程以及超深井温度场控制等方面。从尿素SCR***化学反应机理出发，研究***的偏微分方程模型(PDE)，选择合适的逼近方法简化模型，利用滚动优化控制处理时滞和不确定性的优势，是得到高精度和较强适应能力的控制策略的有效途径。

根据尿素SCR***控制需求与控制问题的提炼与描述，建立偏微分模型。依据该***化学反应特性，选定工业中的常用一阶双曲型***描述流体流动过程。选择Haar小波作为空间基函数，将分布参数***在配置点上转化为集中参数***。将一阶空间偏导数投影到拟Haar小波基上，采用前向Eular法对时间变量离散，得到***的自回归模型。在处理***边界问题时，采用对称延拓的方法抑制边界效应。选择标准二次优化性能指标，将***转化为一个有约束滚动优化问题，设计基于低阶模型的非线性预测控制器，有效处理***的输入-输出约束。最后，采用了粒子群优化方法求解该约束非线性二次规划问题。同时，复杂多变的运行工况和环境给控制带来了大范围的不确定性和时变因素，这就要求指定的控制算法具有自适应的功能，能够随着工况和环境参数变化而进行自调整。因此，选择了自适应滚动优化控制方法优化排放。

仿真在TESIS公司开发的enDYNA发动机精确建模软件平台上进行。选择瞬态标准测试循环ECE，并将氨覆盖率的理想值设定为连续的阶跃信号，验证本发明提出的偏微分模型预测控制器的效果。仿真结果如下图所示。初始阶段氨覆盖率的理想值为0.6；第60秒时刻，氨覆盖率的理想值阶跃到0.8；在第120秒时刻该值阶跃下降到0.7。尿素喷射器约束为0.003mol/s，执行情况在图的下部位置。结果显示，基于尾气排放SCR后处理***的实际情况，本发明设计的控制器能够完成较好的跟踪控制效果。

Claims (1)

1.一种SCR***尿素喷射控制方法，其特征在于：

a、尿素SCR***偏微分建模

基于尿素SCR***各主要物质沿轴向传输特性，主要参考了动态传输反应模型：

$\frac{\∂C_k^g}{\∂t}=-\frac{\∂J_k^g}{\∂z}-K_m^k(C_k^g-C_k^w),---\left(1\right)$

磁通密度公式：

J_k＝θ_zC_k， (2)

以及一种物质浓度在催化器表面的动态及稳态假设得到的：

$-K_m^k(C_k^g-C_k^w)==\frac{(n_{k,in}^g-n_{k,out}^g)}{{\mathrm{\ϵV}}_c}+r_j,---\left(3\right)$

建立得到***的一阶双曲型PDE模型：

$\frac{\∂C_{{\mathrm{NO}}_x}^g}{\∂t}=-{\mathrm{\θ}}_z\frac{\∂C_{{\mathrm{NO}}_x}^g}{\∂z}+\frac{\left(n_{{\mathrm{NO}}_x,in}^g-n_{{\mathrm{NO}}_x,out}^g\right)}{{\mathrm{\ϵV}}_c}-r_{SCR},$

$\frac{{\mathrm{d\Θ}}_{{\mathrm{NO}}_x}}{dt}=r_{ads}-r_{des}-r_{SCR}-r_{ox},$

$\frac{\∂C_{{\mathrm{NH}}_3}^g}{\∂t}=-{\mathrm{\θ}}_z\frac{\∂C_{{\mathrm{NH}}_3}^g}{\∂z}+\frac{(n_{{\mathrm{NH}}_3,in}^g-n_{{\mathrm{NH}}_3,out}^g)}{{\mathrm{\ϵV}}_c}-r_{ads}+r_{des}.---\left(4\right)$

其中，

$r_{ads}=c_s{S}_c{\mathrm{\α}}_{prob}\sqrt{\frac{RT}{2{\mathrm{\πM}}_{{\mathrm{NH}}_3}}}{C}_{{\mathrm{NH}}_3}(1-{\mathrm{\Θ}}_{{\mathrm{NH}}_3}),$

$r_{des}=c_s{k}_{des}{e}^{-\frac{E_{a,des}}{RT}}{\mathrm{\Θ}}_{{\mathrm{NH}}_3},$

$r_{SCR}=c_s{\mathrm{RTk}}_{SCR}{e}^{\left(\frac{-E_{a,\mathrm{Re}d}}{RT}\right)}{C}_{{\mathrm{NO}}_x}-{\mathrm{\Θ}}_{{\mathrm{NH}}_3},$

$r_{ox}=c_c{k}_{Ox}{e}^{\left(\frac{-E_{a,Ox}}{RT}\right)}{\mathrm{\Θ}}_{{\mathrm{NH}}_3}.---\left(5\right)$

进一步整理得到：

$\frac{\∂C_{{\mathrm{NO}}_x}^g}{\∂t}=-{\mathrm{\θ}}_z\frac{\∂C_{{\mathrm{NO}}_x}^g}{\∂z}+a_1{n}_{{\mathrm{NO}}_x,in}^g-C_{{\mathrm{NO}}_x}^g(a_0{a}_1{m}_{EG}^{*}T+a_4(T\left){\mathrm{\Θ}}_{{\mathrm{NH}}_3}\right),$

$c_s\frac{{\mathrm{d\Θ}}_{{\mathrm{NH}}_3}}{dt}=a_2\left(T\right)(1-{\mathrm{\Θ}}_{{\mathrm{NH}}_3})C_{{\mathrm{NH}}_3}^g-\[a_3\left(T\right)+a_4\left(T\right)C_{{\mathrm{NO}}_x}^g+a_5\left(T\right)\]{\mathrm{\Θ}}_{{\mathrm{NH}}_3},$

$\frac{\∂C_{{\mathrm{NH}}_3}^g}{\∂t}=-{\mathrm{\θ}}_z\frac{\∂C_{{\mathrm{NH}}_3}^g}{\∂z}+a_1{n}_{{\mathrm{NH}}_{3,in}}^g-C_{{\mathrm{NH}}_3}^g\[a_0{a}_1{m}_{EG}^{*}T+a_2\left(T\right)(1-{\mathrm{\Θ}}_{{\mathrm{NH}}_3})\]+a_3\left(T\right){\mathrm{\Θ}}_{{\mathrm{NH}}_3}.---\left(6\right)$

其中，

$a_0=\frac{R_{S,EG}}{P_{amb}},a_1=\frac{1}{{\mathrm{\ϵV}}_c},$

$a_2\left(T\right)=c_s{S}_c{\mathrm{\α}}_{\mathrm{Pr}ob}\sqrt{\frac{RT}{2{\mathrm{\πM}}_{{\mathrm{NH}}_3}}},$

$a_3\left(T\right)=c_s{k}_{Des}{e}^{\left(\frac{-E_{a,Des}}{RT}\right)},$

$a_4\left(T\right)=c_s{\mathrm{RTk}}_{SCR}{e}^{\left(\frac{-E_{a,SCR}}{RT}\right)},$

$a_5\left(T\right)=c_s{k}_{Ox}{e}^{\left(\frac{-E_{a,Ox}}{RT}\right)},$

$a_6=\frac{c_{p,EG}}{c_{p,c}{m}_c},a_7=\frac{{\mathrm{\ϵ}}_{rad,scr}{\mathrm{\σ}}_{sb}{A}_{rad,scr}}{c_{p,c}{m}_c}.---\left(7\right)$

表1和表2分别显示了模型中所有常量和变量的相关定义及参数名义参考值：

表1常量命名法

表2变量命名法

将T，C_NOx等看做可测量变量，考虑的状态量仅为和定义状态变量时变参数控制输入变量被控输出为得到***面向控制模型：

$\left\{\begin{array}{c}{\stackrel{\·}{x}}_1=f_{11}(x,p)+f_{12}(x,p)x_2,\\ {\stackrel{\·}{x}}_2=f_{21}(x,p)+a_{1}u-{\mathrm{\θ}}_z\frac{\∂x_2}{\∂z},\end{array}\right.---\left(8\right)$

其中，

$\left\{\begin{array}{c}{f}_{11}\left(x,p\right)=-\frac{1}{c_s}\[a_3\left(T\right)+a_4\left(T\right)C_{{\mathrm{NO}}_x}+a_5\left(T\right)\]x_1,\\ f_{12}\left(x,p\right)=\frac{1}{c_s}{a}_2\left(T\right)\left(1-x_1\right),\\ f_{21}\left(x,p\right)=-x_2\[a_0{a}_1{m}_{EG}^{*}T+a_2\left(T\right)\left(1-x_1\right)\]+a_3\left(T\right)x_1,\end{array}\right.---\left(9\right)$

进一步整理成状态空间形式：

$\begin{array}{c}\frac{\∂X(z,t)}{\∂t}=\left[\begin{array}{cc}0& 0\\ 0& -{\mathrm{\θ}}_z\end{array}\right]\frac{\∂X(z,t)}{\∂z}+\left[\begin{array}{c}{f}_{11}\left(x,p\right)+f_{12}\left(x,p\right)x_2\\ f_{21}\left(x,p\right)+a_{1}u\end{array}\right],\\ y(z,t)=\[\begin{array}{cc}1& 0\end{array}\]X(z,t),\end{array}---\left(10\right)$

上述一阶双曲型***表示为以下状态空间形式：

$\frac{\∂X(z,t)}{\∂t}=A\frac{\∂X(z,t)}{\∂z}+B\left(X\right(z,t),U(z,t\left)\right)---\left(11\right)$

其中，

$\begin{array}{c}A=\left[\begin{array}{cc}0& 0\\ 0& -{\mathrm{\θ}}_z\end{array}\right],\\ B\left(X\left(z,t\right),U\left(z,t\right)\right)=\left[\begin{array}{c}{f}_{11}\left(x,p\right)+f_{12}\left(x,p\right)x_2\\ f_{21}\left(x,p\right)+a_{1}u\end{array}\right],\\ C=\left[\begin{array}{cc}1& 0\end{array}\right].\end{array}---\left(12\right);$

b、偏微分***预测控制器设计：

选取Haar小波尺度函数φ(x)，将时空变量展开如下形式：

$X(z,t)\≈X_j(z,t)=\underset{n=0}{\overset{2^j}{\Σ}}{X}_j(z_n,t)\·{\mathrm{\φ}}_j(z-z_n)---\left(13\right)$

一阶空间偏导数表示为：

$\frac{\∂X(z,t)}{\∂z}\≈\frac{\∂X_j^{\left(1\right)}(z,t)}{\∂z}=\underset{n=0}{\overset{2^j}{\Σ}}{X}_j(z_n,t)\·{\mathrm{\φ}}_j^{\left(1\right)}(z-z_n)---\left(14\right)$

采用Eular前向差分法对时间偏导数进行离散：

$\frac{\∂X(z_n,t)}{\∂t}\≈\frac{\∂X_j^{\left(1\right)}(z_n,t)}{\∂t}\≈\frac{X_j(z_n,i+1)-X_j(z_n,i)}{\mathrm{\Δ}t}---\left(15\right)$

应用上述方法，同时将时空偏导数替换得到：

$\begin{array}{c}\frac{X_j(z_n,i+1|i)-X_j(z_n,i)}{\mathrm{\Δ}t}=A\·\underset{n=0}{\overset{2^j}{\Σ}}{X}_j(z_n,i)\·{\mathrm{\φ}}_j^{\left(1\right)}(z-z_n)\\ +B(\underset{n=0}{\overset{2^j}{\Σ}}{X}_j(z_n,i)\·{\mathrm{\φ}}_j(z-z_n),U(z,i\left)\right)\end{array}---\left(16\right)$

对上式进行整理，得到状态变量(z_n,i+1|i)的预测值为：

$\begin{array}{c}{X}_j(z_n,i+1|i)=X_j(z_n,i)+\mathrm{\Δ}t\·(A\·\underset{n=0}{\overset{2^j}{\Σ}}{X}_j(z_n,i)\·{\mathrm{\φ}}_j^{\left(1\right)}(z-z_n)\\ +B\left(\underset{n=0}{\overset{2^j}{\Σ}}{X}_j\right(z_n,i)\·{\mathrm{\φ}}_j(z-z_n),U(z,i)))\end{array}---\left(17\right)$

由y＝C·X得P步向前预测输出值：

$\begin{array}{c}{\hat{y}}_j(z_n,i+1|i)=y_j(z_n,i)+\mathrm{\Δ}t\·(A\·\underset{n=0}{\overset{2^j}{\Σ}}{y}_j(z_n,i)\·{\mathrm{\φ}}_j^{\left(1\right)}(z-z_n)\\ +B\left(\underset{n=0}{\overset{2^j}{\Σ}}{y}_j\right(z_n,i)\·{\mathrm{\φ}}_j(z-z_n),U(z,i)))\\ .\\ .\\ .\\ {\hat{y}}_j(z_n,i+P+1|i+P-1)={\hat{y}}_j(z_n,i+P-1|i+P-2)+\mathrm{\Δ}t\·\\ (A\·\underset{n=0}{\overset{2^j}{\Σ}}{\hat{y}}_j\left(z_n,i+P-1|i+P-2\right)\·{\mathrm{\φ}}_j^{\left(1\right)}\left(z-z_n\right)+\\ B\left(\underset{n=0}{\overset{2^j}{\Σ}}{\hat{y}}_j\right(z_n,i)\·{\mathrm{\φ}}_j(z-z_n),U(z,i\left)\right))\end{array}---\left(18\right)$

为了快速而平滑地得到***设定值，定义如下参考轨迹：

y^*(z_n,k+i|k) (19)

预测控制输出为：

${\hat{y}}_c(z,k+i|k)=\hat{y}(z,k+i|k)+d(z,k|k)---\left(20\right)$

定义***误差函数为：

$e(z,k+i|k)=y^{*}(z,k+i|k)-{\hat{y}}_c(z,k+i|k)---\left(21\right).$