王鑫, 徐祖华
, 赵均, 邵之江
浙江大学 控制科学与工程学院, 流程生产质量优化与控制国际联合研究中心, 杭州 310027
收稿日期：2018-10-10
基金项目：国家重点研发计划（2017YFB0603703）；国家自然科学基金-浙江两化融合联合基金项目（U1509209）；国家自然科学基金项目（61773340）；中央高校基本科研业务费专项（2018QNA5011）
作者简介：王鑫(1993-), 男, 硕士研究生
通信作者：徐祖华, 副教授, E-mail:
zhxu@zju.edu.cn 摘要：在分布曲线控制中，传统的积分平方误差指标仅考虑输出曲线与目标曲线的误差面积，忽略了分布曲线的内在特征结构。该文从曲线相似度度量出发，提出了一种基于梯度特征的分布曲线模型预测控制算法。该算法采用B样条模型描述分布曲线对象，基于梯度特征度量曲线相似度，综合曲线数值信息与梯度信息构造优化命题，并通过复合梯形方法对优化命题进行离散化，求解得到最优控制策略。仿真结果表明：该算法可提高曲线切换过程中输出曲线与目标曲线的相似度，实现了曲线形状的自然过渡。
关键词：模型预测控制分布曲线曲线相似度梯度特征
Gradient feature-based model predictive controlalgorithm of distribution processes
WANG Xin
, XU Zuhua
, ZHAO Jun
, SHAO Zhijiang
National Center for International Research on Quality-Targeted Process Optimization and Control, College of Control Science and Engineering, Zhejiang University, Hangzhou 310027, China
Abstract: In the control of distribution processes, the traditional integral square error performance index only considers the area between the output curve and the target curve, which ignores the structural features of the distribution curve. A gradient feature-based model predictive control algorithm that takes into account the curve similarities is developed for distribution processes. The algorithm first models the distribution process curve with B-splines. Then, the algorithm quantifies the similarity between the curves based on gradient features and optimizes the design by combining numerical and gradient information. The composite trapezoidal rule is then used to discretize the optimization proposition. Finally, the optimization proposition is solved to get the optimal solution. Simulations show that this algorithm improves the similarity between the output curve and the target curve during curve switching with natural transitions of the curve shape.
Key words: model predictive controldistribution processcurve similaritygradient feature
分布曲线控制问题在很多工业应用场合中都有涉及，如温度场分布、纸张密度分布等。以聚合产品的生产过程为例，在传统工业控制中，常采用宏观指标作为被控目标，如平均相对分子质量或熔融指数等。但是，宏观指标不足以包含足够信息，如涂料和纸张涂层需要准确指定相对分子质量分布
[1-2]。当分布曲线形状具有双峰或尾部较高时，使用平均相对分子质量可能会引起误导。链长分布不同的聚合物可能具有相同的分散性指数，但最终性质明显不同
[3]。因此，有必要研究以完整分布曲线为目标的控制算法。
近十多年来，国内外****对分布曲线控制问题进行了大量的研究。王宏
[4]提出了随机分布控制系统模型，利用B样条函数逼近随机系统输出的概率密度函数(probability density function, PDF)，并先后设计了鲁棒控制
[5]、自适应控制
[6]、迭代学习控制
[7]等一系列控制算法。基于有理平方根B样条近似方法，Yue等
[8]通过以权重动态模型中的实际权重代替伪权重，降低了建模时的计算负荷，提出了利用测量数据建立动态输出PDF模型的方法和针对完整分布曲线的预测控制算法。为同时控制聚合物组成和相对分子质量分布，Vicente等
[9]采用迭代动态规划算法进行开环控制，取得了良好的控制效果。Ali等
[3]提出了一种非线性模型预测控制算法, 用于在线控制聚合物相对分子质量分布，并通过Kalman滤波器引入相对分子质量离线测量值。针对聚合反应相对分子质量分布对象，申珊华等
[10]利用Legendre正交多项式组合神经网络建立灰箱模型，并提出了基于分布函数矩的预测控制算法。针对管式反应器内的温度分布控制问题，曹柳林等
[11]将线性递归神经网络和B样条神经网络相结合，建立了温度分布模型，并提出了扩展积分方差控制指标，设计了自适应控制策略。周靖林等
[12]以每个时刻的分布跟踪误差满足某一确定的上限作为控制目标，设计了对系统不确定性不敏感的鲁棒跟踪控制器。Buehler等
[13]针对具有无界随机不确定性的非线性系统，提出了一种基于Lyapunov的随机非线性模型预测控制算法，用Fokker-Planck方程来描述概率分布函数的动态特性，用Hellinger距离量化PDF间的相似性。康岳群等
[14]基于积分平方误差(integral square error, ISE)和稳态可行性分析，将传统的预测控制算法推广到分布曲线对象，提出一种分布曲线无偏预测控制算法。
已有分布曲线控制方法大多根据输出曲线和目标曲线的积分平方误差来完成分布曲线的最优切换控制，仅考虑曲线间的误差面积，忽略了分布曲线的内在几何特征结构，如波峰、波谷等，曲线切换过程中易出现形状过渡不自然的现象，难以得到令人满意的控制品质。在实现分布曲线无偏跟踪的基础上，优化曲线的切换过程，使曲线几何特征变换过程单调平滑，将有助于提高控制系统性能。模型预测控制(model predictive control, MPC)算法已成功应用在石油、化工等过程工业中，能有效处理复杂过程的多变量约束控制问题。本文提出的分布曲线模型预测控制算法采用B样条函数建立对象模型，基于曲线梯度特征度量输出曲线与目标曲线的相似度，综合曲线数值信息与梯度信息构造优化命题，最后通过求解优化命题，得到最优控制策略，从而实现曲线切换过程中几何特征的自然过渡。
1 分布曲线对象模型本文采用的分布曲线对象模型，包括静态非线性和动态线性两部分。其中：静态非线性利用B样条基函数描述，动态线性部分的输出为B样条基函数的权值。模型表达式如下：
| ${\gamma _k}(y) = \mathit{\boldsymbol{D}}{(y)^{\rm{T}}}\mathit{\boldsymbol{v}}(k) + {e_k}(y),y \in [a,b].$ | (1) |
| $\begin{array}{*{20}{c}}{{\mathit{\boldsymbol{x}}_{\rm{m}}}\left( k \right) = {\mathit{\boldsymbol{A}}_{\rm{m}}}{\mathit{\boldsymbol{x}}_{\rm{m}}}\left( {k - 1} \right) + {\mathit{\boldsymbol{B}}_{\rm{m}}}\mathit{\boldsymbol{u}}\left( {k - 1} \right),}\\{\mathit{\boldsymbol{v}}\left( k \right) = {\mathit{\boldsymbol{C}}_{\rm{m}}}{\mathit{\boldsymbol{x}}_{\rm{m}}}\left( k \right).}\end{array}$ | (2) |
其中：
D(
y)=[
D1(
y),
D2(
y), …,
Dnv(
y)]
T是B样条基函数向量，
v(
k)=[
v1(
k),
v2(
k), …,
vnv(
k)]
T为权值向量，
xm(
k)为状态向量，
u(
k)为输入向量，
γk(
y)和
ek(
y)分别为
k时刻系统输出和不可测输出扰动的分布函数。
令Δ
xm(
k)=
xm(
k)-
xm(
k-1)，
η(
k)=
v(
k)，
Δu(
k)=
u(
k)-
u(
k-1)，将方程(2)转化为式(3)的增量表示形式：
| $\begin{array}{*{20}{c}}{\left[ {\begin{array}{*{20}{c}}{\Delta {\mathit{\boldsymbol{x}}_{\rm{m}}}(k + 1)}\\{\mathit{\boldsymbol{\eta }}(k + 1)}\end{array}} \right] = \left[ {\begin{array}{*{20}{c}}{{\mathit{\boldsymbol{A}}_{\rm{m}}}}&{\bf{0}}\\{{\mathit{\boldsymbol{C}}_{\rm{m}}}{\mathit{\boldsymbol{A}}_{\rm{m}}}}&I\end{array}} \right]\left[ {\begin{array}{*{20}{c}}{\Delta {\mathit{\boldsymbol{x}}_{\rm{m}}}(k)}\\{\mathit{\boldsymbol{\eta }}(k)}\end{array}} \right] + }\\{\left[ {\begin{array}{*{20}{c}}{{\mathit{\boldsymbol{B}}_{\rm{m}}}}\\{{\mathit{\boldsymbol{C}}_{\rm{m}}}{\mathit{\boldsymbol{B}}_{\rm{m}}}}\end{array}} \right]\Delta \mathit{\boldsymbol{u}}(k),}\\{\mathit{\boldsymbol{v}}(\mathit{\boldsymbol{k}}) = \left[ {\begin{array}{*{20}{l}}{\bf{0}}&\mathit{\boldsymbol{I}}\end{array}} \right]\left[ {\begin{array}{*{20}{c}}{\Delta {\mathit{\boldsymbol{x}}_{\rm{m}}}(k)}\\{\mathit{\boldsymbol{\eta }}(k)}\end{array}} \right].}\end{array}$ | (3) |
令
$x(k) = \left[ {\begin{array}{*{20}{c}}{\Delta {x_{\rm{m}}}(k)}\\{\mathit{\boldsymbol{\eta }}(k)}\end{array}} \right], \mathit{\boldsymbol{A}} = \left[ {\begin{array}{*{20}{c}}{{\mathit{\boldsymbol{A}}_{\rm{m}}}}&\mathit{\boldsymbol{0}}\\{{\mathit{\boldsymbol{C}}_{\rm{m}}}{\mathit{\boldsymbol{A}}_{\rm{m}}}}&\mathit{\boldsymbol{I}}\end{array}} \right], \;\mathit{\boldsymbol{B}} = \left[ {\begin{array}{*{20}{c}}{{\mathit{\boldsymbol{B}}_{\rm{m}}}}\\{{\mathit{\boldsymbol{C}}_{\rm{m}}}{\mathit{\boldsymbol{B}}_{\rm{m}}}}\end{array}} \right], \;\mathit{\boldsymbol{C}} = \left[ {\begin{array}{*{20}{l}}\mathit{\boldsymbol{0}}&\mathit{\boldsymbol{I}}\end{array}} \right]$，则分布曲线对象模型可以表示为：
| ${\gamma _k}(y) = \mathit{\boldsymbol{D}}{(y)^{\rm{T}}}\mathit{\boldsymbol{v}}(k) + {e_k}(y),y \in [a,b].$ | (4) |
| $\begin{array}{*{20}{c}}{\mathit{\boldsymbol{x}}(k + 1) = \mathit{\boldsymbol{Ax}}(k) + \mathit{\boldsymbol{B}}\Delta \mathit{\boldsymbol{u}}(k),}\\{\mathit{\boldsymbol{v}}(k) = \mathit{\boldsymbol{Cx}}(k).}\end{array}$ | (5) |
2 曲线相似度度量积分平方误差准则是一种常用的量化输出值与设定值偏差大小的指标，取值范围为[0, +∞)，其值越小表明系统输出越接近设定值。对于两条分布范围为[
a,
b]的曲线
x(
l)和
y(
l)，其ISE表达式为
| $\mathrm{ISE}=\int_{a}^{b}[x(l)-y(l)]^{2} \mathrm{d} l.$ | (6) |
ISE准则是基于Euclid距离计算的，而Euclid距离作为典型的距离度量方法，主要关注数值差异的绝对值，对序列中的突变或异常点极其敏感，但不能很好地反映曲线几何特征。此外，ISE指标的数值大小实质为输出曲线与目标曲线间的误差面积，忽视了曲线上点与点之间的关联，ISE指标大小受曲线均值影响较大。这一特点导致当曲线均值相近、形状不同时，ISE指标难以有效反映曲线间的相似度。
与距离度量方法不同，Pearson相似度作为一种典型的相似度度量方法，更关注序列间的趋势性和方向性差异，两个变量的位置和尺度变化不会引起该系数的改变，取值范围为[-1，+1]。其中：0代表无相关性，负值表示负相关，正值为正相关，1代表输出曲线与目标曲线一致。对于长度相等的序列
X和
Y，两者之间的Pearson相似度为
| $\rho (X,Y) = \frac{{{\rm{E}}\left[ {\left( {X - {\mu _X}} \right)\left( {Y - {\mu _Y}} \right)} \right]}}{{{\sigma _X}{\sigma _Y}}}.$ | (7) |
其中：
μX，
μY为
X,
Y的均值；
σX，
σY为
X,
Y的标准差；E表示期望值。
以
图 1中
A、
B两条曲线为例，比较两条曲线与目标曲线的相似度。观察曲线整体形状可以发现，曲线
B具有明显的波峰、波谷特征，与目标曲线更为接近。
表 1为
图 1中曲线根据式(6)和(7)计算得到的两个相似度。根据
表 1中结果，采用ISE指标难以区分曲线
A、
B与目标曲线的相似度，而Pearson相似度显示出曲线
B与目标曲线更为相似。
表 1 两条曲线的相似度
| 相似度指标 | 曲线A | 曲线B |
| ISE | 0.072 | 0.075 |
| Pearson相似度 | 0.927 | 0.994 |
与ISE指标相比，Pearson相似度计算的是相对距离，对绝对数值不敏感，是从方向上区分差异的。与此相似，梯度值的计算与绝对数值无关，且曲线的局部梯度信息反映了该局部区间内曲线的变化方向与速率。因此，可利用完整的曲线梯度信息描述曲线沿分布方向的变化趋势，并利用输出曲线与目标曲线间梯度信息的偏差量化曲线间的趋势相似度。本文以曲线斜率作为梯度信息，则曲线
X、
Y的梯度相似度
SG(
X,
Y)的表达式为
| $S_{G}(X, Y)=\int_{a}^{b}[\dot{x}(l)-\dot{y}(l)]^{2} \mathrm{d} l.$ | (8) |
特别地，针对模型预测控制算法设计问题，为实现曲线的无偏控制和切换过程中几何特征自然过渡，滚动优化性能指标函数应包含曲线数值信息和梯度信息。本文所提MPC算法采用ISE指标与梯度相似度的加权和量化输出曲线与目标曲线的偏差大小，取值范围为[0, +∞)，值越小表明曲线越相似，0代表输出曲线与目标曲线一致。
| $\begin{aligned} S_{\mathrm{C}}(X, Y)=& \int_{a}^{b}\left\{m_{1}[x(l)-y(l)]^{2}+\right.\\ m_{2} &[\dot{x}(l)-\dot{y}(l)]^{2} \} \mathrm{d} l \end{aligned}.$ | (9) |
式(9)中：
m1和
m2为数值信息和梯度信息的权重系数；
a和
b为曲线分布范围的最小值和最大值。
3 基于梯度特征的模型预测控制3.1 预测模型在采样时刻
k，根据对象状态的估计值
x(
k|
k)和
k至
k+
M-1时刻的控制增量Δ
U(
k)，可计算
k至
k+
P时刻的权值预测值
V(
k+1),
| $\mathit{\boldsymbol{V}}(k + 1) = \mathit{\boldsymbol{Fx}}(k|k) + \mathit{\boldsymbol{ \boldsymbol{\varPhi} }}\Delta \mathit{\boldsymbol{U}}(k).$ | (10) |
其中各变量定义如下所示：
| $\begin{array}{*{20}{c}}{\mathit{\boldsymbol{V}}(k + 1) = \left[ {\mathit{\boldsymbol{v}}{{(k + 1|k)}^{\rm{T}}}} \right.,}\\{\mathit{\boldsymbol{v}}{{(k + 2|k)}^{\rm{T}}}, \cdots ,\mathit{\boldsymbol{v}}{{(k + P|k)}^{\rm{T}}}{]^{\rm{T}}},}\end{array}$ | (11) |
| $\begin{array}{*{20}{c}}{\Delta \mathit{\boldsymbol{U}}(k) = \left[ {\Delta \mathit{\boldsymbol{u}}{{(k)}^{\rm{T}}},\Delta \mathit{\boldsymbol{u}}{{(k + 1)}^{\rm{T}}}, \cdots ,} \right.}\\{\Delta \mathit{\boldsymbol{u}}{{(k + M - 1)}^{\rm{T}}}{]^{\rm{T}}},}\end{array}$ | (12) |
| $\boldsymbol{F}=\left[\boldsymbol{C A}, \boldsymbol{C A}^{2}, \cdots, \boldsymbol{C A}^{P}\right]^{\mathrm{T}},$ | (13) |
| $\mathit{\boldsymbol{ \boldsymbol{\varPhi} }} = \left[ {\begin{array}{*{20}{c}}{\mathit{\boldsymbol{CB}}}&{}&{}&{}\\{\mathit{\boldsymbol{CAB}}}&{\mathit{\boldsymbol{CB}}}&{}&{}\\ \vdots & \vdots &{}&{}\\{\mathit{\boldsymbol{C}}{\mathit{\boldsymbol{A}}^{P - 1}}\mathit{\boldsymbol{B}}}&{\mathit{\boldsymbol{C}}{\mathit{\boldsymbol{A}}^{P - 2}}\mathit{\boldsymbol{B}}}& \cdots &{\mathit{\boldsymbol{C}}{\mathit{\boldsymbol{A}}^{P - M}}\mathit{\boldsymbol{B}}}\end{array}} \right].$ | (14) |
3.2 滚动优化根据式(9)，综合曲线数值信息和梯度信息，构造滚动优化性能指标函数，即式(15)。在采样时刻
k，求解相应优化命题可确定最优控制序列。
| $\begin{array}{*{20}{c}}{J = \sum\limits_{i = 1}^P {{q_i}} \int_a^b {{{\left[ {{\gamma _{k + i}}(y) - {\gamma _{{\rm{tar}}}}(y)} \right]}^2}} {\rm{d}}y + \sum\limits_{i = 1}^P {{w_i}} \cdot }\\{\int_a^b {{{\left[ {{{\dot \gamma }_{k + i}}(y) - {{\dot \gamma }_{{\rm{tar}}}}(y)} \right]}^2}} {\rm{d}}y + \sum\limits_{i = 0}^{M - 1} {{r_i}} \Delta {{\bf{u}}^2}(k + i).}\end{array}$ | (15) |
将性能指标函数
J分为以下3个部分：
曲线数值信息,
| ${J_1} = \sum\limits_{i = 1}^P {{q_i}} \int_a^b {{{\left[ {{\gamma _{k + i}}(y) - {\gamma _{{\mathop{\rm tar}\nolimits} }}(y)} \right]}^2}} {\rm{d}}y.$ | (16) |
曲线梯度信息,
| ${J_2} = \sum\limits_{i = 1}^P {{w_i}} \int_a^b {{{\left[ {{{\dot \gamma }_{k + i}}(y) - {{\dot \gamma }_{{\rm{ tar }}}}(y)} \right]}^2}} {\rm{d}}y$ | (17) |
控制作用软约束,
| ${J_3} = \sum\limits_{i = 0}^{M - 1} {{r_i}} \Delta {\mathit{\boldsymbol{u}}^2}(k + i).$ | (18) |
其中：
γtar(
y)、
${\dot \gamma _{{\rm{ tar }}}}\left( y \right)$为目标曲线的数值和斜率；
γk+i(
y)、
${\dot \gamma _{{\rm{ tar }}}}\left( y \right)$为
k+
i时刻预测输出曲线的数值和斜率；
qi、
wi、
ri分别为数值特征误差、梯度特征误差和控制增量在
k+
i时刻的权重系数。
把曲线分布范围[
a,
b]等分为
n段，根据复合梯形方法分别对积分项
J1和
J2进行离散化得到：
| $\begin{array}{*{20}{c}}{{J_1} = \frac{h}{2}\sum\limits_{i = 1}^P {{q_i}} \left\{ {{{\left[ {{\gamma _{k + i}}(a) - {\gamma _{{\rm{ tar }}}}(a)} \right]}^2} + } \right.}\\{2\sum\limits_{j = 1}^{n - 1} {{{\left[ {{\gamma _{k + i}}\left( {{y_j}} \right) - {\gamma _{{\rm{tar}}}}\left( {{y_j}} \right)} \right]}^2}} + }\\{{{\left[ {{\gamma _{k + i}}(b) - {\gamma _{{\mathop{\rm tar}\nolimits} }}(b)} \right]}^2}\} ,}\end{array}$ | (19) |
| $\begin{array}{*{20}{c}}{{J_2} = \frac{h}{2}\sum\limits_{i = 1}^P {{w_i}} \left\{ {{{\left[ {{{\dot \gamma }_{k + i}}(a) - {{\dot \gamma }_{{\rm{tar}}}}(a)} \right]}^2} + } \right.}\\{2\sum\limits_{j = 1}^{n - 1} {{{\left[ {{{\dot \gamma }_{k + i}}\left( {{y_j}} \right) - {{\dot \gamma }_{{\rm{tar}}}}\left( {{y_j}} \right)} \right]}^2}} + }\\{{{\left[ {{{\dot \gamma }_{k + i}}(b) - {{\dot \gamma }_{{\mathop{\rm tar}\nolimits} }}(b)} \right]}^2}\} .}\end{array}$ | (20) |
其中
h=(
b-
a)/
n。曲线梯度信息用相邻两点间的斜率表示：
| $\dot{\gamma}_{k+i}\left(y_{j}\right)=\frac{\gamma_{k+i}\left(y_{j}\right)-\gamma_{k+i}\left(y_{j-1}\right)}{h},$ | (21) |
| $\dot{\gamma}_{\mathrm{tar}}\left(y_{j}\right)=\frac{\gamma_{\mathrm{tar}}\left(y_{j}\right)-\gamma_{\mathrm{tar}}\left(y_{j-1}\right)}{h}.$ | (22) |
将式(19)、(20)代入式(15)，即可得到性能指标函数的离散化表达式。代入预测方程(10)，可得到性能指标函数，
| $\begin{array}{*{20}{c}}{J = \left\| {\begin{array}{*{20}{c}}{{\mathit{\boldsymbol{D}}_{\rm{I}}}\left( a \right)\left[ {\mathit{\boldsymbol{Fx}}\left( k \right) + \mathit{\boldsymbol{ \boldsymbol{\varPhi} }}\Delta \mathit{\boldsymbol{U}}} \right] + {\mathit{\boldsymbol{e}}_{\rm{I}}}\left( a \right) - {\mathit{\boldsymbol{\gamma }}_{{\rm{tar}},{\rm{I}}}}\left( a \right)}\\{{\mathit{\boldsymbol{D}}_{\rm{I}}}\left( {{y_1}} \right)\left[ {\mathit{\boldsymbol{Fx}}\left( k \right) + \mathit{\boldsymbol{ \boldsymbol{\varPhi} }}\Delta \mathit{\boldsymbol{U}}} \right] + {\mathit{\boldsymbol{e}}_{\rm{I}}}\left( {{y_1}} \right) - {\mathit{\boldsymbol{\gamma }}_{{\rm{tar}},{\rm{I}}}}\left( {{y_1}} \right)}\\ \vdots \\{{\mathit{\boldsymbol{D}}_{\rm{I}}}\left( b \right)\left[ {\mathit{\boldsymbol{Fx}}\left( k \right) + \mathit{\boldsymbol{ \boldsymbol{\varPhi} }}\Delta \mathit{\boldsymbol{U}}} \right] + {\mathit{\boldsymbol{e}}_{\rm{I}}}\left( b \right) - {\mathit{\boldsymbol{\gamma }}_{{\rm{tar}},{\rm{I}}}}\left( b \right)}\end{array}} \right\|_\mathit{\boldsymbol{Q}}^2 + }\\{\left\| {\begin{array}{*{20}{c}}{{\mathit{\boldsymbol{D}}_{\rm{G}}}\left( a \right)\left[ {\mathit{\boldsymbol{Fx}}\left( k \right) + \mathit{\boldsymbol{ \boldsymbol{\varPhi} }}\Delta \mathit{\boldsymbol{U}}} \right] + {\mathit{\boldsymbol{e}}_{\rm{G}}}\left( a \right) - {\mathit{\boldsymbol{\gamma }}_{{\rm{tar}},{\rm{G}}}}\left( a \right)}\\{{\mathit{\boldsymbol{D}}_{\rm{G}}}\left( {{y_{\rm{I}}}} \right)\left[ {\mathit{\boldsymbol{Fx}}\left( k \right) + \mathit{\boldsymbol{ \boldsymbol{\varPhi} }}\Delta \mathit{\boldsymbol{U}}} \right] + {\mathit{\boldsymbol{e}}_{\rm{G}}}\left( {{y_{\rm{1}}}} \right) - {\mathit{\boldsymbol{\gamma }}_{{\rm{tar}},{\rm{G}}}}\left( {{y_1}} \right)}\\ \vdots \\{{\mathit{\boldsymbol{D}}_{\rm{G}}}\left( b \right)\left[ {\mathit{\boldsymbol{Fx}}\left( k \right) + \mathit{\boldsymbol{ \boldsymbol{\varPhi} }}\Delta \mathit{\boldsymbol{U}}} \right] + {\mathit{\boldsymbol{e}}_{\rm{G}}}\left( b \right) - {\mathit{\boldsymbol{\gamma }}_{{\rm{tar}},{\rm{G}}}}\left( b \right)}\end{array}} \right\|_\mathit{\boldsymbol{W}}^2 + }\\{\left\| {\Delta \mathit{\boldsymbol{U}}} \right\|_\mathit{\boldsymbol{R}}^2.}\end{array}$ | (23) |
其中
R为控制增量权重矩阵。其他变量表达式为：
| ${\mathit{\boldsymbol{Q}}_P} = h\left[ {{q_1},{q_2}, \cdots ,{q_P}} \right],$ | (24) |
| $\mathit{\boldsymbol{Q}} = {\rm{block - diag}}\left( {\frac{1}{2}{\mathit{\boldsymbol{Q}}_P},{\mathit{\boldsymbol{Q}}_P}, \cdots ,{\mathit{\boldsymbol{Q}}_P},\frac{1}{2}{\mathit{\boldsymbol{Q}}_P}} \right),$ | (25) |
| ${\mathit{\boldsymbol{D}}_{\rm{I}}}\left( {{y_j}} \right) = {\rm{block - diag}}\left( {\mathit{\boldsymbol{D}}{{\left( {{y_j}} \right)}^{\rm{T}}}, \cdots ,\mathit{\boldsymbol{D}}{{\left( {{y_j}} \right)}^{\rm{T}}}} \right),$ | (26) |
| ${\mathit{\gamma }_{{\rm{tar,I}}}}\left( {{y_j}} \right) = {\left[ {{\gamma _{{\rm{tar}}}}\left( {{y_j}} \right), \cdots ,{\gamma _{{\rm{tar}}}}\left( {{y_j}} \right)} \right]^{\rm{T}}},$ | (27) |
| ${\mathit{\boldsymbol{e}}_{\rm{I}}}\left( {{y_j}} \right) = {\left[ {{e_k}\left( {{y_j}} \right), \cdots ,{e_k}\left( {{y_j}} \right)} \right]^{\rm{T}}},$ | (28) |
| ${\mathit{\boldsymbol{W}}_P} = \frac{1}{h}\left[ {{w_1},{w_2}, \cdots ,{w_P}} \right],$ | (29) |
| $\mathit{\boldsymbol{W}} = {\rm{block - diag}}\left( {\frac{1}{2}{\mathit{\boldsymbol{W}}_P},{\mathit{\boldsymbol{W}}_P}, \cdots ,{\mathit{\boldsymbol{W}}_P},\frac{1}{2}{\mathit{\boldsymbol{W}}_P}} \right),$ | (30) |
| ${\mathit{\boldsymbol{\gamma }}_{{\rm{tar}},{\rm{G}}}}\left( {{y_j}} \right) = \left[ {\begin{array}{*{20}{c}}{{\gamma _{{\rm{tar}}}}\left( {{y_j}} \right) - {\gamma _{{\rm{tar}}}}\left( {{y_{j - 1}}} \right)}\\ \vdots \\{{\gamma _{{\rm{tar}}}}\left( {{y_j}} \right) - {\gamma _{{\rm{tar}}}}\left( {{y_{j - 1}}} \right)}\end{array}} \right],$ | (31) |
| $\begin{array}{*{20}{c}}{{\mathit{\boldsymbol{D}}_{\rm{G}}}\left( {{y_j}} \right) = }\\{\left[ {\begin{array}{*{20}{c}}{\mathit{\boldsymbol{D}}{{\left( {{y_j}} \right)}^{\rm{T}}} - \mathit{\boldsymbol{D}}{{\left( {{y_{j - 1}}} \right)}^{\rm{T}}}}&{}&{}\\{}& \ddots &{}\\{}&{}&{\mathit{\boldsymbol{D}}{{\left( {{y_j}} \right)}^{\rm{T}}} - \mathit{\boldsymbol{D}}{{\left( {{y_{j - 1}}} \right)}^{\rm{T}}}}\end{array}} \right],}\end{array}$ | (32) |
| ${\mathit{\boldsymbol{e}}_{\rm{G}}}\left( {{y_j}} \right) = \left[ {\begin{array}{*{20}{c}}{{e_k}\left( {{y_j}} \right) - {e_k}\left( {{y_{j - 1}}} \right)}\\ \vdots \\{{e_k}\left( {{y_j}} \right) - {e_k}\left( {{y_{j - 1}}} \right)}\end{array}} \right].$ | (33) |
在输入输出无约束条件的情况下，由极值必要条件?
J/?
U=0求得最优解表达式(34)，
| $\begin{array}{*{20}{c}}{\Delta \mathit{\boldsymbol{U}} = - {{\left( {\frac{1}{2}\mathit{\boldsymbol{f}}(a) + \sum\limits_{j = 1}^{n - 1} \mathit{\boldsymbol{f}} \left( {{y_j}} \right) + \frac{1}{2}\mathit{\boldsymbol{f}}(b) + \mathit{\boldsymbol{R}}} \right)}^{ - 1}} \cdot }\\{\left[ {\frac{1}{2}\mathit{\boldsymbol{g}}(a) + \sum\limits_{j = 1}^{n - 1} \mathit{\boldsymbol{g}} \left( {{y_j}} \right) + \frac{1}{2}\mathit{\boldsymbol{g}}(b)} \right].}\end{array}$ | (34) |
其中：
| $\begin{array}{*{20}{c}}{\mathit{\boldsymbol{f}}\left( {{y_j}} \right) = {\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}^{\rm{T}}}{\mathit{\boldsymbol{D}}_{\rm{I}}}{{\left( {{y_j}} \right)}^{\rm{T}}}{\mathit{\boldsymbol{Q}}_P}{\mathit{\boldsymbol{D}}_{\rm{I}}}\left( {{y_j}} \right)\mathit{\boldsymbol{ \boldsymbol{\varPhi} }} + }\\{{\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}^{\rm{T}}}{\mathit{\boldsymbol{D}}_{\rm{G}}}{{\left( {{y_j}} \right)}^{\rm{T}}}{\mathit{\boldsymbol{W}}_P}{\mathit{\boldsymbol{D}}_{\rm{G}}}\left( {{y_j}} \right)\mathit{\boldsymbol{ \boldsymbol{\varPhi} }},}\end{array}$ | (35) |
| $\begin{array}{*{20}{c}}{\mathit{\boldsymbol{g}}\left( {{y_j}} \right) = {\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}^{\rm{T}}}{\mathit{\boldsymbol{D}}_{\rm{I}}}{{\left( {{y_j}} \right)}^{\rm{T}}}{\mathit{\boldsymbol{Q}}_P}\left( {{\mathit{\boldsymbol{D}}_{\rm{I}}}\left( {{y_j}} \right)\mathit{\boldsymbol{Fx}}(k) + } \right.}\\{\left. {{\mathit{\boldsymbol{e}}_{\rm{I}}}\left( {{y_j}} \right) - {\mathit{\boldsymbol{\gamma }}_{{\rm{tar,I}}}}\left( {{y_j}} \right)} \right) + {\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}^{\rm{T}}}{\mathit{\boldsymbol{D}}_{\rm{G}}}{{\left( {{y_j}} \right)}^{\rm{T}}} \cdot }\\{{\mathit{\boldsymbol{W}}_P}\left( {{\mathit{\boldsymbol{D}}_{\rm{G}}}\left( {{y_j}} \right)\mathit{\boldsymbol{Fx}}\left( k \right) + {\mathit{\boldsymbol{e}}_{\rm{G}}}\left( {{y_j}} \right) - {\mathit{\boldsymbol{\gamma }}_{{\rm{tar}},{\rm{G}}}}\left( {{y_j}} \right)} \right).}\end{array}$ | (36) |
3.3 反馈校正由于存在模型失配、不可测干扰等因素，模型预测输出与对象真实输出可能存在偏差。利用测量数据可进行反馈校正，具体步骤如下：
1) 利用测量数据得到k时刻对象的实测分布曲线
γk(
y),
y∈
a,
b。
2) 根据
k-1时刻的预测输出曲线
γk|k-1(
y)与实测分布曲线
γk(
y)，计算预测误差曲线
ek(
y)。
3) 将预测误差曲线
ek(
y)代入性能指标函数，用于下一步的滚动优化计算。
4 仿真结果本文针对如下的分布曲线对象设计控制器：
| ${\gamma _k}(y) = \mathit{\boldsymbol{D}}{(y)^{\rm{T}}}\mathit{\boldsymbol{v}}(k) + {e_k}(y),y \in [a,b].$ | (37) |
| $\begin{array}{*{20}{c}}{{\mathit{\boldsymbol{x}}_{\rm{m}}}(k) = {\mathit{\boldsymbol{A}}_{\rm{m}}}{\mathit{\boldsymbol{x}}_{\rm{m}}}(k - 1) + {\mathit{\boldsymbol{B}}_{\rm{m}}}\mathit{\boldsymbol{u}}(k - 1),}\\{\mathit{\boldsymbol{v}}(k) = {\mathit{\boldsymbol{C}}_{\rm{m}}}{\mathit{\boldsymbol{x}}_{\rm{m}}}(k).}\end{array}$ | (38) |
| ${\mathit{\boldsymbol{A}}_{\rm{m}}} = \left[ {\begin{array}{*{20}{c}}{0.5}&0&0\\0&{0.8}&0\\0&0&{0.6}\end{array}} \right],$ | (39) |
| ${\mathit{\boldsymbol{B}}_{\rm{m}}} = \left[ {\begin{array}{*{20}{c}}{0.45}&{0.1}\\{0.1}&{0.65}\\{0.4}&{0.4}\end{array}} \right],$ | (40) |
| ${\mathit{\boldsymbol{C}}_{\rm{m}}} = \left[ {\begin{array}{*{20}{c}}1&{0.2}&{0.2}\\{0.2}&1&{0.2}\\{0.2}&{0.2}&1\end{array}} \right].$ | (41) |
其中
D(
y)是3阶B样条基函数，按式(42)和(43)递推得到：
| ${D_{j,0}} = \left\{ {\begin{array}{*{20}{l}}{1,}&{{x_j} < x < {x_{j + 1}};}\\{0,}&{其他}.\end{array}} \right.$ | (42) |
| $\begin{array}{*{20}{c}}{{D_{j,q}}(x) = \frac{{x - {x_j}}}{{{x_{j + q}} - {x_j}}}{D_{j,q - 1}}(x) + }\\{\frac{{{x_{j + q + 1}} - x}}{{{x_{j + q + 1}} - {x_{j + 1}}}}{D_{j + 1,q - 1}}(x).}\end{array}$ | (43) |
其中：
x表示变量，
q表示B样条基函数的阶次，
j表示第
j个基函数。
图 2为B样条基函数曲线与目标曲线示意图。仿真目标曲线表达式为
| ${\gamma _{{\rm{tar}}}}(y) = 1.26{D_1}(y) + 0.3{D_2}(y) + 1.22{D_3}(y).$ | (44) |
作为对比，除本文所提算法(Gradient-MPC)外，对基于ISE指标的MPC算法(ISE-MPC)进行仿真，其滚动优化性能指标为
| $\begin{array}{*{20}{c}}{J = \sum\limits_{i = 1}^P {{q_i}} \int_a^b {{{\left[ {{\gamma _{k + i}}\left( y \right) - {\gamma _{{\rm{tar}}}}\left( y \right)} \right]}^2}{\rm{d}}y} + }\\{\sum\limits_{i = 0}^{M - 1} {{r_i}\Delta {\mathit{\boldsymbol{u}}^2}\left( {k + i} \right)} .}\end{array}$ | (45) |
预测控制中，通常根据过程的上升时间选择优化时域
P，控制系统的快速性和稳定性受控制时域
M影响。ISE-MPC中，选择优化时域
P=6，控制时域
M=3。为简单起见，输出偏差权重
qi和控制增量权重系数
ri取值均为1。Gradient-MPC的参数中，优化时域
P=6，控制时域
M=3，控制增量权重系数
ri=1。
qi和
wi(
i=1, 2, …,
P)为数值特征误差和梯度特征误差的加权系数，为保持对控制增量软约束效果与ISE-MPC相近，取
qi=0.5,
wi=0.008。
定义
k时刻控制器控制增量的能量值为
Eu(
k)=||Δ
Uk||
R2，
图 3中实线、虚线分别为ISE-MPC和Gradient-MPC仿真中
Eu的变化曲线。可见，仿真过程中两控制器对控制增量的软约束效果相近。
图 4为控制输入信号的响应曲线, 两条曲线的变化幅度基本相似，而Gradient-MPC的输入信号更快、更平稳地达到稳态值。
图 5和
6分别为Gradient-MPC和ISE-MPC的仿真结果。在相同的采样时刻，
图 5中的输出曲线具有更明显的双峰特征，整体形状更接近目标曲线，曲线切换过程中曲线几何特征过渡更自然。经过17个控制周期，Gradient-MPC输出曲线已跟踪上目标曲线，而ISE-MPC输出曲线与目标曲线仍存在明显偏差。
图 7为两算法仿真输出曲线与目标曲线的相似度(ISE、Pearson相似度)的变化曲线。两算法仿真结果的ISE指标变化曲线基本一致，而Gradient-MPC对应仿真结果的Pearson相似度明显优于ISE-MPC。
5 总结本文采用B样条基函数与状态空间模型来描述分布曲线对象。针对ISE指标无法有效反映曲线形状特征的问题，在控制器设计中考虑曲线特征的控制问题，利用曲线数值信息和梯度信息，重新构造性能指标函数。最后，通过求解优化命题，得到最优控制策略。
本文基于曲线梯度特征设计MPC算法，实现了曲线切换过程中形状自然过渡。与基于ISE指标的MPC算法相比，切换过程中本文算法的输出曲线与目标曲线的Pearson相似度更高，证明了本文算法的有效性。
参考文献 | [1] | CONGALIDIS J P, RICHARDS J R. Process control of polymerization reactors:An industrial perspective[J]. Polymer Reaction Engineering, 1998, 6(2): 71-111. DOI:10.1080/10543414.1998.10744484
|
| [2] | SAYER C, ARZAMENDI G, ASUA J M, et al. Dynamic optimization of semicontinuous emulsion copolymerization reactions:Composition and molecular weight distribution[J]. Computers & Chemical Engineering, 2001, 25(4-6): 839-849.
|
| [3] | ALI M A, AJBAR E A A H, ALHUMAIZI K. Control of molecular weight distribution of polyethylene in gas-phase fluidized bed reactors[J]. Korean Journal of Chemical Engineering, 2010, 27(1): 364-372.
|
| [4] | WANG H. Bounded dynamic stochastic distributions:Modelling and control[M]. London, UK: Springer-Verlag, 2000.
|
| [5] | WANG H. Robust control of the output probability density functions for multivariable stochastic systems[C]//Proceedings of the 37th IEEE Conference on Decision and Control. Piscataway, USA: IEEE Press, 1998: 1305-1310.
|
| [6] | WANG H. Model reference adaptive control of the output stochastic distributions for unknown linear stochastic systems[J]. International Journal of Systems Science, 1999, 30(7): 707-715. DOI:10.1080/002077299292029
|
| [7] | YUE H, WANG H, ZHANG J F. Shaping of molecular weight distribution by iterative learning probability density function control strategies[J]. Proceedings of the Institution of Mechanical Engineers, Part I:Journal of Systems and Control Engineering, 2008, 222(7): 639-653. DOI:10.1243/09596518JSCE584
|
| [8] | ZHANG J F, YUE H, ZHOU J L. Predictive PDF control in shaping of molecular weight distribution based on a new modeling algorithm[J]. Journal of Process Control, 2015, 30: 80-89. DOI:10.1016/j.jprocont.2014.12.009
|
| [9] | VICENTE M, SAYER C, LEIZA J R, et al. Dynamic optimization of non-linear emulsion copolymerization systems:Open-loop control of composition and molecular weight distribution[J]. Chemical Engineering Journal, 2002, 85(2-3): 339-349. DOI:10.1016/S1385-8947(01)00180-2
|
| [10] | 申珊华, 曹柳林, 王晶. 基于分布函数矩的聚合物分子量分布预测控制[J]. 化工学报, 2013, 64(12): 4379-4384.
SHEN S H, CAO L L, WANG J. Predictive control of molecular weight distribution in polymerization reaction based on moment of MWD[J]. Journal of Chemical Industry and Engineering (China), 2013, 64(12): 4379-4384. (in Chinese)
|
| [11] | CAO L L, LI D Z, ZHANG C Y, et al. Control and modeling of temperature distribution in a tubular polymerization process[J]. Computers & Chemical Engineering, 2007, 31(11): 1516-1524.
|
| [12] | ZHOU J L, LI G, WANG H. Robust tracking controller design for non-Gaussian singular uncertainty stochastic distribution systems[J]. Automatica, 2014, 50(4): 1296-1303. DOI:10.1016/j.automatica.2014.02.032
|
| [13] | BUEHLER E A, PAULSON J A, MESBAH A. Lyapunov-based stochastic nonlinear model predictive control: Shaping the state probability distribution functions[C]//2016 American Control Conference. Philadelphia, USA: American Automatic Control Council, 2016: 5389-5394.
|
| [14] | 康岳群, 徐祖华, 赵均, 等. 分布曲线对象的无偏模型预测控制算法[J]. 化工学报, 2016, 67(3): 701-706.
KANG Y Q, XU Z H, ZHAO J, et al. Offset-free model-predictive control algorithm of distribution process[J]. Journal of Chemical Industry and Engineering (China), 2016, 67(3): 701-706. (in Chinese)
|
