一种非自衡对象的分布式动态矩阵控制方法

摘要

本发明公开了一种非自衡对象的分布式动态矩阵控制方法。本发明首先通过采集阶跃响应数据建立含非自衡对象的多变量过程的矩阵模型向量，再将多变量过程的在线优化实施问题转化成各个小规模子系统的优化实施问题。然后选取合适的性能指标，通过不断迭代获取各智能体的纳什最优解，进而得到各智能体动态矩阵控制器的参数，再对每个智能体实施该时刻的即时控制律，并将时域滚动到下一时刻，重复上述优化过程，从而完成整个系统的优化任务。本发明在保证较高控制精度和稳定性的前提下，能够有效弥补了传统DDMC方法在含非自衡对象的多变量过程控制中的不足，并满足实际工业过程的需求。

说明书

技术领域

本发明属于自动化技术领域，涉及一种非自衡对象的分布式动态矩阵控制(DDMC)方法。

背景技术

实际过程中广泛存在着大量复杂高维的大规模系统，采用集中式的整体求解对计算机的性能和处理速度等往往要求很高，与实际工业系统中必须考虑的经济性相违背。分布式动态矩阵控制(DDMC)作为分布式预测控制(DMPC)的一个主要分支，综合利用计算机通信技术和控制理论，把一个复杂大规模系统的在线求解问题分散到各个子系统中去分布实现，有效降低了问题的规模和复杂性，能很好的控制存在多变量、强耦合、不确定的被控对象，改善了系统控制性能。然而在实际工业过程中，存在许多含非自衡对象的多变量过程，如一部分储罐、锅炉汽包水位、精馏塔液位等等。由于非自衡对象的传递函数中含有典型的积分环节，进而导致被控对象在定值阶跃下的响应趋于无穷，这就使得传统DDMC算法无法直接应用。如果能够在实际过程中对传统DDMC方法进行改进，便能有效弥补传统DDMC方法在含非自衡对象的多变量过程控制中的不足，使得DMPC在实际应用中得到进一步延伸和发展。

发明内容

本发明目的是针对传统DDMC方法在含非自衡对象的多变量过程控制中的不足之处，提出了一种非自衡对象的DDMC方法。

该方法首先通过采集阶跃响应数据建立含非自衡对象的多变量过程的矩阵模型向量，挖掘出基本的对象特性，再将多变量过程的在线优化实施问题转化成各个小规模子系统的优化实施问题，并结合多智能体中的理论和思想，把网络环境下的每个子系统看作为一个智能体，各智能体之间通过网络进行物质、能量及信息通信。然后通过将一种针对非自衡对象改进转移矩阵的方法和一种新的误差校正方法相结合，并选取合适的性能指标，基于纳什优化的思想通过不断迭代获取各智能体的纳什最优解，进而得到各智能体动态矩阵控制器的参数，再对每个智能体实施该时刻的即时控制律，并将时域滚动到下一时刻，重复上述优化过程，从而完成整个系统的优化任务。

本发明的技术方案是通过数据采集、模型建立、预测机理、优化等手段，确立了一种非自衡对象的分布式动态矩阵控制方法，利用该方法在保证较高控制精度和稳定性的前提下，能够有效弥补了传统DDMC方法在含非自衡对象的多变量过程控制中的不足，并满足实际工业过程的需求。

本发明方法的步骤包括：

步骤1.通过非自衡对象的实时阶跃响应数据建立相应的动态矩阵模型向量，具体方法是：

1.1依据分布式预测控制思想，将一个N输入N输出非自衡对象的大规模系统分散为N个智能体子系统；

1.2在稳态工况下，以第j个智能体控制量为输入对第i个智能体输出量进行阶跃响应实验，分别记录第j(1≤j≤N)个输入对第i(1≤i≤N)个输出的阶跃响应曲线；

1.3将步骤1.2得到的阶跃响应曲线进行滤波处理，然后拟合成一条光滑曲线，记录光滑曲线上每个采样时刻对应的阶跃响应数据，第一个采样时刻为T_s，相邻两个采样时刻的间隔时间为T_s，采样时刻顺序为T_s、2T_s、3T_s……；被控对象的阶跃响应数据将在某一个时刻t_L＝I_ijT_s开始呈现定斜率上升，以该时刻的数据为起始点，之前的数据分别记做建立第j个输入对第i个输出之间的阶跃响应模型向量a_ij：

$a_{ij}={[a_1^{ij},a_2^{ij},\mathrm{...},a_{I_{ij}-1}^{ij},a_{I_{ij}}^{ij}+\delta ,a_{I_{ij}}^{ij}+2\delta ,...,a_{L_{ij}}^{ij}]}^T$

$a_{L_{ij}}^{ij}=a_{I_{ij}}^{ij}+(L_{ij}-I_{ij})\delta$

其中T为矩阵的转置符号,δ为阶跃响应数据呈恒定斜率上升之后相邻两个数据之间的恒定差值，L_ij为设定的第j个输入对第i个输出的模型长度，L_ij≥I_ij+1。

步骤2.设计第i个智能体的动态矩阵控制器，具体方法是：

2.1利用步骤1获得的模型向量a_ij建立被控对象的动态矩阵，其形式如下：

其中A_ij为第j个智能体输入对第i个智能体输出的P×M阶动态矩阵，a_ij(k)为第j个输入对第i个输出的阶跃响应数据，P为动态矩阵控制算法的优化时域，M为动态矩阵控制算法的控制时域，L_ij＝L(1≤i≤3,1≤j≤3),M<P<L，N为输入输出个数；

2.2获取第i个智能体当前k时刻的模型预测初始响应值y_i,0(k)

首先，在k-1时刻加入控制增量△u₁(k-1),△u₂(k-1),…,△u_n(k-1)，得到第i个智能体的模型预测值y_i,P(k-1)：

$y_{i,P}(k-1)=y_{i,0}(k-1)+A_{ii,0}{\mathrm{\Delta u}}_i(k-1)+\underset{j=1,j\ne i}{\overset{n}{\Sigma }}{A}_{ij,0}{\mathrm{\Delta u}}_j(k-1)$

其中，

y_i,P(k-1)＝[y_i,1(k|k-1),y_i,1(k+1|k-1),…,y_i,1(k+L-1|k-1)]^T

y_i,0(k-1)＝[y_i,0(k|k-1),y_i,0(k+1|k-1),…,y_i,0(k+L-1|k-1)]^T,

A_ii,0＝[a_ii(1),a_ii(2),…,a_ii(L)]^T,A_ij,0＝[a_ij(1),a_ij(2),…,a_ij(L)]^T

y_i,1(k|k-1),y_i,1(k+1|k-1),…,y_i,1(k+L-1|k-1)分别表示第i个智能体在k-1时刻对k,k+1,…,k+L-1时刻的模型预测值，y_i,0(k|k-1),y_i,0(k+1|k-1),…,y_i,0(k+L-1|k-1)表示k-1时刻对k,k+1,…,k+L-1时刻的初始预测值，A_ii,0,A_ij,0分别为第i个智能体和第j个智能体输入对第i个智能体输出的阶跃响应数据建立的矩阵，△u₁(k-1),△u₂(k-1),…,△u_n(k-1)为k-1时刻各智能体的输入控制量；

然后，可以得到k时刻第i个智能体的模型预测误差值e_i(k)：

e_i(k)＝y_i(k)-y_i,1(k|k-1)

其中y_i(k)表示k时刻测得的第i个智能体的实际输出值；

进一步得到k时刻修正后的模型输出值y_i,cor(k)：

y_i,cor(k)＝y_i,0(k-1)+h₁*e_i(k)+h₂*e_i(k)

其中，

y_i,cor(k)＝[y_i,cor(k|k),y_i,cor(k+1|k),…,y_i,cor(k+L-1|k)]^T,

h₁＝[1,α,…,α]^T,h₂＝[0,1,…,L-1]^T

y_i,cor(k|k),y_i,cor(k+1|k),…,y_i,cor(k+L-1|k)分别表示第i个智能体在k时刻模型的修正值，h₁和h₂为误差补偿的权矩阵，α为误差校正系数，0<α≤1；

最后得到第i个智能体k时刻的模型预测的初始响应值y_i,0(k)：

y_i,0(k)＝Sy_i,cor(k)

其中，S为L×L阶的新的状态转移矩阵，

2.3依据步骤2.1获取第i个智能体在M个连续的控制增量△u_i(k),△u_i(k+1),…,△u_i(k+M-1)下的预测输出值y_i,PM，具体方法是：

$y_{i,PM}\left(k\right)=y_{i,P0}\left(k\right)+A_{ii}{\mathrm{\Delta u}}_{i,M}\left(k\right)+\underset{j=1,j\ne i}{\overset{n}{\Sigma }}{A}_{ij}{\mathrm{\Delta u}}_{j,M}\left(k\right)$

其中，

y_i,PM(k)＝[y_i,M(k+1|k),y_i,M(k+2|k),…,y_i,M(k+P|k)]^T

y_i,P0(k)＝[y_i,0(k+1|k),y_i,0(k+2|k),…,y_i,0(k+P|k)]^T

△u_i,M(k)＝[△u_i(k),△u_i(k+1),…,△u_i(k+M-1)]^T

△u_j,M(k)＝[△u_j(k),△u_j(k+1),…,△u_j(k+M-1)]^T

y_i,P0(k)是y_i,0(k)的前P项，y_i,0(k+1|k),y_i,0(k+2|k),…,y_i,0(k+P|k)为k时刻对k+1,k+2,…,k+P时刻的模型预测输出值；

2.4建立非自衡对象第i个智能体动态矩阵控制器的性能指标J_i(k)和参考轨迹ω_i(k)，形式如下：

minJ_i(k)＝(ω_i(k)-y_i,PM(k))^TQ_i(ω_i(k)-y_i,PM(k))+△u_i,M(k)^TR_i△u_i,M(k)

ω_i(k)＝[ω_i(k+1),ω_i(k+2),…,ω_i(k+P)]^T

ω_i(k+ε)＝β^εy(k)+(1-β^ε)c(k)(ε＝1,2,…,P)

其中为误差加权矩阵，为控制加权矩阵，和分别为Q_i,R_i中的权重系数，ω_i(k)为第i个智能体的参考轨迹，β为参考轨迹的柔化系数；

2.5依据纳什最优的思想，由步骤2.4中性能指标得到第i个智能体当前k时刻的纳什最优解：

${\mathrm{\Delta u}}_{i,M}^{*}\left(k\right)=D_{ii}\left({\omega }_i\right(k)-y_{i,P0}(k)-\underset{j=1,j\ne i}{\overset{n}{\Sigma }}{A}_{ij}{\mathrm{\Delta u}}_{j,M}^{*}(k\left)\right)$

其中：

2.6由步骤2.2到2.5可以得到在k时刻智能体i的新一轮迭代最优解为：

${\mathrm{\Delta u}}_{i,M}^{l+1}\left(k\right)=D_{ii}\left({\omega }_i\right(k)-y_{i,P0}(k)-\underset{j=1,j\ne i}{\overset{n}{\Sigma }}{A}_{ij}{\mathrm{\Delta u}}_{j,M}^l(k\left)\right)$

进一步得到整个系统在k时刻的最优控制律：

${\mathrm{\Delta u}}_M^{l+1}\left(k\right)=D_1\left(\omega \right(k)-y_{P0}(k\left)\right)-D_0{\mathrm{\Delta u}}_M^l\left(k\right)$

其中：

${\mathrm{\Delta u}}_M^{l+1}\left(k\right)={[{\mathrm{\Delta u}}_{1,M}^{l+1}\left(k\right),{\mathrm{\Delta u}}_{2,M}^{l+1}\left(k\right),...,{\mathrm{\Delta u}}_{n,M}^{l+1}\left(k\right)]}^T$

${\mathrm{\Delta u}}_M^l\left(k\right)={[{\mathrm{\Delta u}}_{1,M}^l\left(k\right),{\mathrm{\Delta u}}_{2,M}^l\left(k\right),...,{\mathrm{\Delta u}}_{n,M}^l\left(k\right)]}^T$

ω(k)＝[ω₁(k),ω₂(k),…,ω_n(k)]^T，y_P0(k)＝[y_1,P0(k),y_2,P0(k),…,y_n,P0(k)]^T

2.7将第i个智能体k时刻的纳什最优解首项作为即时控制律△u_i(k)，得到智能体i的实际控制量u_i(k)＝u_i(k-1)+△u_i(k)作用于第i个智能体；

2.8在下一时刻，重复步骤2.2到2.7继续求解第i个智能体的即时控制律△u_i(k+1)，进而得到整个系统的最优解△u(k+1)，并依次循环。

本发明提出了一种非自衡对象的DDMC方法。该方法在传统DDMC方法的基础上，将一种针对非自衡对象改进转移矩阵的方法和一种新的误差校正方法结合，在保证较高控制精度和稳定性的前提下，有效弥补了传统DDMC方法在含非自衡对象的多变量过程控制中的不足，并满足实际工业过程的需求。

具体实施方式

以锅炉汽包水位控制为例：

锅炉汽包水位控制系统是一个典型的带积分环节的多变量非自衡对象，调节手段采用控制给水阀阀门开度。

步骤1.通过锅炉汽包水位对象的实时阶跃响应数据建立相应的动态矩阵模型向量，具体方法是：

1.1依据分布式预测控制思想，将一个3输入3输出锅炉汽包水位对象的大规模系统分散为3个子系统；

1.2在稳态工况下，以第j个锅炉给水阀阀门开度为输入对第i个锅炉汽包水位进行阶跃响应实验，分别记录第j(1≤j≤3)个输入对第i(1≤i≤3)个输出的阶跃响应曲线；

1.3将步骤1.2得到的阶跃响应曲线进行滤波处理，然后拟合成一条光滑曲线，记录光滑曲线上每个采样时刻对应的阶跃响应数据，第一个采样时刻为T_s，相邻两个采样时刻的间隔时间为T_s，采样时刻顺序为T_s、2T_s、3T_s……；锅炉汽包水位的阶跃响应数据将在某一个时刻t_L＝I_ijT_s开始呈现定斜率上升，以该时刻的数据为起始点，之前的数据分别记做建立第j个锅炉输入对第i个锅炉输出之间的阶跃响应模型向量a_ij：

$a_{ij}={[a_1^{ij},a_2^{ij},\mathrm{...},a_{I_{ij}-1}^{ij},a_{I_{ij}}^{ij}+\delta ,a_{I_{ij}}^{ij}+2\delta ,...,a_{L_{ij}}^{ij}]}^T$

$a_{L_{ij}}^{ij}=a_{I_{ij}}^{ij}+(L_{ij}-I_{ij})\delta$

其中T为矩阵的转置符号,δ为阶跃响应数据呈恒定斜率上升之后相邻两个数据之间的恒定差值，L_ij为设定的第j个输入对第i个输出的模型长度，L_ij≥I_ij+1。

步骤2.设计第i个锅炉的动态矩阵控制器，具体是：

2.1利用步骤1获得的模型向量a_ij建立锅炉汽包水位的动态矩阵，其形式如下：

其中A_ij为第j个锅炉输入对第i个锅炉输出的P×M阶动态矩阵，a_ij(k)为第j个锅炉输入对第i个锅炉输出的阶跃响应数据，P为动态矩阵控制算法的优化时域，M为动态矩阵控制算法的控制时域，L_ij＝L(1≤i≤3,1≤j≤3),M<P<L，N＝3为输入输出个数；

2.2获取第i个锅炉当前k时刻的模型预测初始响应值y_i,0(k)

首先，在k-1时刻加入控制增量△u₁(k-1),△u₂(k-1),…,△u_n(k-1)(n＝3)，得到第i个锅炉的模型预测值y_i,P(k-1)：

$y_{i,P}(k-1)=y_{i,0}(k-1)+A_{ii,0}{\mathrm{\Delta u}}_i(k-1)+\underset{j=1,j\ne i}{\overset{n}{\Sigma }}{A}_{ij,0}{\mathrm{\Delta u}}_j(k-1)$

其中，

y_i,P(k-1)＝[y_i,1(k|k-1),y_i,1(k+1|k-1),…,y_i,1(k+L-1|k-1)]^T

y_i,0(k-1)＝[y_i,0(k|k-1),y_i,0(k+1|k-1),…,y_i,0(k+L-1|k-1)]^T,

A_ii,0＝[a_ii(1),a_ii(2),…,a_ii(L)]^T,A_ij,0＝[a_ij(1),a_ij(2),…,a_ij(L)]^T

y_i,1(k|k-1),y_i,1(k+1|k-1),…,y_i,1(k+L-1|k-1)分别表示第i个锅炉在k-1时刻对k,k+1,…,k+L-1时刻的模型预测值，y_i,0(k|k-1),y_i,0(k+1|k-1),…,y_i,0(k+L-1|k-1)表示k-1时刻对k,k+1,…,k+L-1时刻的初始预测值，A_ii,0,A_ij,0分别为第i个锅炉和第j个锅炉输入对第i个锅炉输出的阶跃响应数据建立的矩阵，△u₁(k-1),△u₂(k-1),…,△u_n(k-1)为k-1时刻各锅炉的输入给水阀阀门开度增量；

然后，可以得到k时刻第i个锅炉的模型预测误差值e_i(k)：

e_i(k)＝y_i(k)-y_i,1(k|k-1)

其中y_i(k)表示k时刻测得的第i个锅炉的实际输出值；

进一步得到k时刻修正后的模型输出值y_i,cor(k)：

y_i,cor(k)＝y_i,0(k-1)+h₁*e_i(k)+h₂*e_i(k)

其中，

y_i,cor(k)＝[y_i,cor(k|k),y_i,cor(k+1|k),…,y_i,cor(k+L-1|k)]^T,

h₁＝[1,α,…,α]^T,h₂＝[0,1,…,L-1]^T

y_i,cor(k|k),y_i,cor(k+1|k),…,y_i,cor(k+L-1|k)分别表示第i个锅炉在k时刻模型的修正值，h₁和h₂为误差补偿的权矩阵，α为误差校正系数，0<α≤1；

最后得到第i个锅炉k时刻的模型预测的初始响应值y_i,0(k)：

y_i,0(k)＝Sy_i,cor(k)

其中，S为L×L阶的新的状态转移矩阵，

2.3依据步骤2.1获取第i个锅炉在M个连续的控制增量△u_i(k),△u_i(k+1),…,△u_i(k+M-1)下的预测输出值y_i,PM，具体方法是：

$y_{i,PM}\left(k\right)=y_{i,P0}\left(k\right)+A_{ii}{\mathrm{\Delta u}}_{i,M}\left(k\right)+\underset{j=1,j\ne i}{\overset{n}{\Sigma }}{A}_{ij}{\mathrm{\Delta u}}_{j,M}\left(k\right)$

其中，

y_i,PM(k)＝[y_i,M(k+1|k),y_i,M(k+2|k),…,y_i,M(k+P|k)]^T

y_i,P0(k)＝[y_i,0(k+1|k),y_i,0(k+2|k),…,y_i,0(k+P|k)]^T

△u_i,M(k)＝[△u_i(k),△u_i(k+1),…,△u_i(k+M-1)]^T

△u_j,M(k)＝[△u_j(k),△u_j(k+1),…,△u_j(k+M-1)]^T

y_i,P0(k)是y_i,0(k)的前P项，y_i,0(k+1|k),y_i,0(k+2|k),…,y_i,0(k+P|k)为k时刻对k+1,k+2,…,k+P时刻的模型预测输出值；

2.4建立锅炉汽包水位对象第i个锅炉动态矩阵控制器的性能指标J_i(k)和参考轨迹ω_i(k)，形式如下：

minJ_i(k)＝(ω_i(k)-y_i,PM(k))^TQ_i(ω_i(k)-y_i,PM(k))+△u_i,M(k)^TR_i△u_i,M(k)

ω_i(k)＝[ω_i(k+1),ω_i(k+2),…,ω_i(k+P)]^T

ω_i(k+ε)＝β^εy(k)+(1-β^ε)c(k)(ε＝1,2,…,P)

其中为误差加权矩阵，为控制加权矩阵，和分别为Q_i,R_i中的权重系数，ω_i(k)为第i个锅炉的参考轨迹，β为参考轨迹的柔化系数；

2.5依据纳什最优的思想，由步骤2.4中性能指标得到第i个锅炉当前k时刻的纳什最优解：

${\mathrm{\Delta u}}_{i,M}^{*}\left(k\right)=D_{ii}\left({\omega }_i\right(k)-y_{i,P0}(k)-\underset{j=1,j\ne i}{\overset{n}{\Sigma }}{A}_{ij}{\mathrm{\Delta u}}_{j,M}^{*}(k\left)\right)$

其中：

2.6由步骤2.2到2.5可以得到在k时刻锅炉i的新一轮迭代最优解为：

${\mathrm{\Delta u}}_{i,M}^{l+1}\left(k\right)=D_{ii}\left({\omega }_i\right(k)-y_{i,P0}(k)-\underset{j=1,j\ne i}{\overset{n}{\Sigma }}{A}_{ij}{\mathrm{\Delta u}}_{j,M}^l(k\left)\right)$

进一步得到整个系统在k时刻的最优控制律：

${\mathrm{\Delta u}}_M^{l+1}\left(k\right)=D_1\left(\omega \right(k)-y_{P0}(k\left)\right)-D_0{\mathrm{\Delta u}}_M^l\left(k\right)$

其中：

${\mathrm{\Delta u}}_M^{l+1}\left(k\right)={[{\mathrm{\Delta u}}_{1,M}^{l+1}\left(k\right),{\mathrm{\Delta u}}_{2,M}^{l+1}\left(k\right),...,{\mathrm{\Delta u}}_{n,M}^{l+1}\left(k\right)]}^T$

${\mathrm{\Delta u}}_M^l\left(k\right)={[{\mathrm{\Delta u}}_{1,M}^l\left(k\right),{\mathrm{\Delta u}}_{2,M}^l\left(k\right),...,{\mathrm{\Delta u}}_{n,M}^l\left(k\right)]}^T$

ω(k)＝[ω₁(k),ω₂(k),…,ω_n(k)]^T，y_P0(k)＝[y_1,P0(k),y_2,P0(k),…,y_n,P0(k)]^T

2.7将第i个锅炉k时刻的纳什最优解首项作为即时控制律△u_i(k)，得到锅炉i的实际给水阀阀门开度u_i(k)＝u_i(k-1)+△u_i(k)作用于第i个锅炉；

2.8在下一时刻，重复步骤2.2到2.7继续求解第i个锅炉的即时控制律△u_i(k+1)，进而得到整个系统的最优控制律△u(k+1)，并依次循环。