CN106444388A - 一种焦炭炉炉膛压力的分布式pid型动态矩阵控制方法 - Google Patents

Abstract

本发明公开了一种焦炭炉炉膛压力的分布式PID型动态矩阵控制方法。本发明首先通过采集实时阶跃响应数据建立被控对象的动态矩阵模型向量，再将大规模***的在线优化实现问题转化为各个小规模子***的优化求解问题，并把网络环境下的每个子***看作一个智能体，同时各智能体之间通过网络通信完成信息交换。然后通过引入PID算子建立一种改进的性能指标，并依据纳什最优思想来设计各智能体的PID型动态矩阵控制器，再将当前时刻的即时控制律作用于每个智能体，并将时域滚动到下一时刻，重复上述优化过程，从而完成整个大规模***的优化任务。本发明在保证***整体控制品质的同时，有效弥补了传统 DDMC方法的不足，并提高了控制参数设计的自由度。

Description

一种焦炭炉炉膛压力的分布式PID型动态矩阵控制方法

技术领域

本发明属于自动化技术领域，涉及一种焦炭炉炉膛压力的分布式PID型动态矩阵控制方法。

背景技术

实际过程中往往存在着大量复杂高维的大规模***，其变量和约束的数目通常有几十个甚至上百个之多，因而针对这样一类大规模预测控制的在线实施问题便显得尤为重要。分布式动态矩阵控制(DDMC)作为分布式预测控制(DMPC)的一个主要分支，综合利用网络通信技术和控制理论，把一个复杂大规模***的在线求解问题分散到各个子***中分布实现，有效降低了问题的规模和复杂性，能很好的控制存在多变量、强耦合、不确定的被控对象，改善了***的控制性能。然而传统DDMC方法通常在线计算量较大，采样周期也不宜过小，从而不能使扰动受到快速抑制。因而对于模型阶次，环境扰动、非线性等方面存在不确定性较大的控制***，常规的DDMC方法通常难以达到预期的控制效果。传统的PID控制由于其控制结构简单、操作方便、鲁棒性强等优点，一直被广泛的应用于实际工业过程中。如果能够在实际过程中将PID控制同DDMC方法相结合，在保证***良好控制性能的同时，又提高了控制参数设计的自由度。

发明内容

本发明目的是针对传统DDMC方法应用于模型阶次，环境扰动、非线性等方面存在不确定性较大的控制***的不足之处，提出了一种焦炭炉炉膛压力的分布式PID型动态矩阵控制方法。该方法将PID控制同常规的DDMC方法相结合，弥补了传统DDMC方法的不足，并保证了良好的控制性能。本发明方法首先通过采集实时阶跃响应数据建立被控对象的动态矩阵模型向量，再将大规模***的在线优化实现问题转化为各个小规模子***的优化求解问题，并把网络环境下的每个子***看作一个智能体，同时各智能体之间通过网络通信完成信息交换，以保证***整体的控制品质。然后通过引入PID算子建立一种改进的性能指标，并依据纳什最优思想来设计各智能体的PID型动态矩阵控制器，再将当前时刻的即时控制律作用于每个智能体，并将时域滚动到下一时刻，重复上述优化过程，从而完成整个大规模***的优化任务。

本发明的技术方案是通过数据采集、模型建立、预测机理、优化等手段，确立了一种焦炭炉炉膛压力的分布式PID型动态矩阵控制方法，利用该方法能很好的处理***模型阶次，环境扰动、非线性等方面存在不确定性较大的控制问题，并在保证良好控制品质的同时，提高控制参数设计的自由度。

本发明方法的步骤包括：

步骤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＝L_ijT_s后趋于平稳，当a_ij(t)(t>L_ij)与a_ij(L_ij)的误差和测量误差有相同的数量级时，即可认为a_ij(L_ij)近似等于阶跃响应的稳态值。建立第j个输入对第i个输出的阶跃响应模型向量a_ij：

a_ij＝[a_ij(1),a_ij(2),…,a_ij(L_ij)]^T

其中T为矩阵的转置符号,a_ij(k)为t＝kT_s时刻第j个输入对第i个输出的阶跃响应采样值,L_ij为第j个输入对第i个输出的建模时域。

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

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

其中A_ij为第j个智能体输入对第i个智能体输出的P×M阶动态矩阵，P为动态矩阵控制算法的优化时域，M为动态矩阵控制算法的控制时域，且L_ij＝L(1≤i≤3,1≤j≤3),M<P<L，L为***的统一建模时域，N＝3为输入输出个数；

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,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)表示第i个智能体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)

其中，

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

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

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

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

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

2.3计算第i个智能体在M个连续的输入控制增量Δu_i(k),Δu_i(k+1),…,Δu_i(k+M-1)下的模型预测输出y_i,PM，具体方法是：

其中，

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)为第i个智能体k时刻对k+1,k+2,…,k+P时刻的初始预测输出值；

2.4建立第i个智能体的性能指标J_i(k)，形式如下：

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

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

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

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

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

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

Δω_i(k+ε)＝ω_i(k+ε)-ω_i(k+ε-1)

Δy_i,M(k+ε|k)＝y_i,M(k+ε|k)-y_i,M(k+ε-1|k)

Δ²ω_i(k+ε)＝Δω_i(k+ε)-Δω_i(k+ε-1)

Δ²y_i,M(k+ε|k)＝Δy_i,M(k+ε|k)-Δy_i,M(k+ε-1|k)

其中ω_i(k+ε)为第i个智能体给定期望输出的参考轨迹，y_i(k)为k时刻第i个智能体的过程实际输出，c(k)为k时刻第i个智能体的期望输出，λ为参考轨迹的柔化因子，分别为第i个智能体的比例系数矩阵、积分系数矩阵、微分系数矩阵；为第i个智能体的控制加权系数矩阵。

2.5将步骤2.4中的性能指标变换为如下形式：

进一步得到

同理可得

其中

引入矩阵

进而有

进一步可将性能指标转化为

其中，

2.6依据纳什最优的思想，以Δu_i,M(k)为控制变量，由求解最优控制律，形式如下：

其中，

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

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

其中：

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

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

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

本发明提出了一种焦炭炉炉膛压力的分布式PID型动态矩阵控制方法。该方法通过采集实时阶跃响应数据建立被控对象的动态矩阵模型向量，并依据纳什优化的思想设计了一种新型分布式动态矩阵控制器，在保证***整体控制品质的同时，有效弥补了传统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＝L_ijT_s后趋于平稳，当a_ij(t)(t>L_ij)与a_ij(L_ij)的误差和测量误差有相同的数量级时，即可认为a_ij(L_ij)近似等于阶跃响应的稳态值。建立第j个输入对第i个输出的阶跃响应模型向量a_ij：

a_ij＝[a_ij(1),a_ij(2),…,a_ij(L_ij)]^T

其中T为矩阵的转置符号,a_ij(k)为t＝kT_s时刻第j个输入对第i个输出的阶跃响应采样值,L_ij为第j个输入对第i个输出的建模时域。

步骤2.设计第i个炉膛的PID型动态矩阵控制器，具体方法如下：

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

其中A_ij为第j个炉膛输入对第i个炉膛输出的P×M阶动态矩阵，P为动态矩阵控制算法的优化时域，M为动态矩阵控制算法的控制时域，且L_ij＝L(1≤i≤3,1≤j≤3),M<P<L，L为***的统一建模时域，N＝3为输入输出个数；

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,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)

其中，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 y_i,cor(k|k),y_i,cor(k+1|k),…,y_i,cor(k+L-1|k)分别表示第i个炉膛在k时刻模型的修正值，h为误差补偿的权矩阵，α为误差校正系数；

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

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

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

2.3计算第i个炉膛在M个连续的阀门开度增量Δu_i(k),Δu_i(k+1),…,Δu_i(k+M-1)下的预测输出值y_i,PM，具体方法是：

其中，

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)为第i个炉膛k时刻对k+1,k+2,…,k+P时刻的初始预测输出值；

2.4建立第i个炉膛的性能指标J_i(k)，形式如下：

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

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

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

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

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

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

Δω_i(k+ε)＝ω_i(k+ε)-ω_i(k+ε-1)

Δy_i,M(k+ε|k)＝y_i,M(k+ε|k)-y_i,M(k+ε-1|k)

Δ²ω_i(k+ε)＝Δω_i(k+ε)-Δω_i(k+ε-1)

Δ²y_i,M(k+ε|k)＝Δy_i,M(k+ε|k)-Δy_i,M(k+ε-1|k)

其中ω_i(k+ε)为第i个炉膛给定期望输出的参考轨迹，y_i(k)为k时刻第i个炉膛的过程实际输出，c(k)为k时刻第i个炉膛的期望输出，λ为参考轨迹的柔化因子，分别为第i个炉膛的比例系数矩阵、积分系数矩阵、微分系数矩阵；为第i个炉膛的控制加权系数矩阵。

2.5将步骤2.4中的性能指标变换为如下形式：

进一步得到

同理可得

其中

引入矩阵

进而有

进一步可将性能指标转化为

其中，

2.6依据纳什最优的思想，以Δu_i,M(k)为控制变量，由求解最优控制律，形式如下：

其中，

2.7由步骤2.2到步骤2.6，可以得到k时刻第i个炉膛的新一轮迭代最优解为：

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

其中：

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

2.8将第i个炉膛k时刻的纳什最优解首项作为即时控制增量Δu_i(k)，得到第i个炉膛的实际控制量u_i(k)＝u_i(k-1)+Δu_i(k)作用于第i个炉膛；

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

Claims (1)

1.一种焦炭炉炉膛压力的分布式PID型动态矩阵控制方法，其特征在于该方法包括以下步骤：

步骤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＝L_ijT_s后趋于平稳，当a_ij(t)(t>L_ij)与a_ij(L_ij)的误差和测量误差有相同的数量级时，即可认为a_ij(L_ij)近似等于阶跃响应的稳态值；建立第j个输入对第i个输出的阶跃响应模型向量a_ij：

a_ij＝[a_ij(1),a_ij(2),…,a_ij(L_ij)]^T

其中T为矩阵的转置符号,a_ij(k)为t＝kT_s时刻第j个输入对第i个输出的阶跃响应采样值,L_ij为第j个输入对第i个输出的建模时域；

步骤2.设计第i个智能体的PID型动态矩阵控制器，具体如下：

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

其中A_ij为第j个智能体输入对第i个智能体输出的P×M阶动态矩阵，P为动态矩阵控制算法的优化时域，M为动态矩阵控制算法的控制时域，且L_ij＝L(1≤i≤3,1≤j≤3),M<P<L，L为***的统一建模时域，N＝3为输入输出个数；

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\&NotEqual;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)表示第i个智能体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)

其中，

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

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

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

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

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

2.3计算第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\&NotEqual;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)为第i个智能体k时刻对k+1,k+2,…,k+P时刻的初始预测输出值；

2.4建立第i个智能体的性能指标J_i(k)，形式如下：

$\begin{array}{c}\min J_i\left(k\right)={\left(w_i\right(k)-y_{i,PM}(k\left)\right)}^T{K}_I^i\left(w_i\right(k)-y_{i,PM}(k\left)\right)+{\left({\mathrm{\&Delta;w}}_i\right(k)-{\mathrm{\&Delta;y}}_{i,PM}(k\left)\right)}^T{K}_p^i\left({\mathrm{\&Delta;w}}_i\right(k)-{\mathrm{\&Delta;y}}_{i,PM}(k\left)\right)\\ {\left({\mathrm{\&Delta;}}^2{w}_i\right(k)-{\mathrm{\&Delta;}}^2{y}_{i,PM}(k\left)\right)}^T{K}_d^i\left({\mathrm{\&Delta;}}^2{w}_i\right(k)-{\mathrm{\&Delta;}}^2{y}_{i,PM}(k\left)\right)+{\mathrm{\&Delta;u}}_{i,M}{\left(k\right)}^T{R}_i{\mathrm{\&Delta;u}}_{i,M}\left(k\right)\end{array}$

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

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

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

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

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

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

Δω_i(k+ε)＝ω_i(k+ε)-ω_i(k+ε-1)

Δy_i,M(k+ε|k)＝y_i,M(k+ε|k)-y_i,M(k+ε-1|k)

Δ²ω_i(k+ε)＝Δω_i(k+ε)-Δω_i(k+ε-1)

Δ²y_i,M(k+ε|k)＝Δy_i,M(k+ε|k)-Δy_i,M(k+ε-1|k)

其中ω_i(k+ε)为第i个智能体给定期望输出的参考轨迹，y_i(k)为k时刻第i个智能体的过程实际输出，c(k)为k时刻第i个智能体的期望输出，λ为参考轨迹的柔化因子， 分别为第i个智能体的比例系数矩阵、积分系数矩阵、微分系数矩阵；为第i个智能体的控制加权系数矩阵；

2.5将步骤2.4中的性能指标变换为如下形式：

$\min J_i\left(k\right)=E_0^i{\left(k\right)}^T{K}_I^i{E}_0^i\left(k\right)+{\mathrm{\&Delta;E}}_0^i{\left(k\right)}^T{K}_p^i{\mathrm{\&Delta;E}}_0^i\left(k\right)+{\mathrm{\&Delta;}}^2{E}_0^i{\left(k\right)}^T{K}_d^i{\mathrm{\&Delta;}}^2{E}_0^i\left(k\right)+{\mathrm{\&Delta;u}}_{i,M}{\left(k\right)}^T{R}_i{\mathrm{\&Delta;u}}_{i,M}\left(k\right)$

$E_0^i\left(k\right)=w_i\left(k\right)-y_{i,PM}\left(k\right)={\&lsqb;e_0^i(k+1),e_0^i(k+2),...,e_0^i(k+P)\&rsqb;}^T$

${\mathrm{\&Delta;E}}_0^i\left(k\right)={\mathrm{\&Delta;w}}_i\left(k\right)-{\mathrm{\&Delta;y}}_{i,PM}\left(k\right)={\&lsqb;{\mathrm{\&Delta;e}}_0^i(k+1),{\mathrm{\&Delta;e}}_0^i(k+2),...,{\mathrm{\&Delta;e}}_0^i(k+P)\&rsqb;}^T$

${\mathrm{\&Delta;}}^2{E}_0^i\left(k\right)={\mathrm{\&Delta;}}^2{w}_i\left(k\right)-{\mathrm{\&Delta;}}^2{y}_{i,PM}\left(k\right)={\&lsqb;{\mathrm{\&Delta;}}^2{e}_0^i(k+1),{\mathrm{\&Delta;}}^2{e}_0^i(k+2),...,{\mathrm{\&Delta;}}^2{e}_0^i(k+P)\mathrm{\&rsqb;}}^T$

进一步得到

$\begin{array}{c}{\mathrm{\&Delta;e}}_0^i(k+\mathrm{\&epsiv;})={\mathrm{\&Delta;\&omega;}}_i(k+\mathrm{\&epsiv;})-{\mathrm{\&Delta;y}}_{i,M}(k+\mathrm{\&epsiv;}|k)\\ ={\mathrm{\&omega;}}_i(k+\mathrm{\&epsiv;})-y_{i,M}(k+\mathrm{\&epsiv;}|k)-\left({\mathrm{\&omega;}}_i\right(k+\mathrm{\&epsiv;}-1)-y_{i,M}(k+\mathrm{\&epsiv;}-1|k\left)\right)\\ =e_0^i(k+\mathrm{\&epsiv;})-e_0^i(k+\mathrm{\&epsiv;}-1)\end{array}$

同理可得

${\mathrm{\&Delta;}}^2{e}_0^i(k+\mathrm{\&epsiv;})={\mathrm{\&Delta;e}}_0^i(k+\mathrm{\&epsiv;})-{\mathrm{\&Delta;e}}_0^i(k+\mathrm{\&epsiv;}-1)$

其中

引入矩阵

进而有

$\left\{\begin{array}{c}{\mathrm{\&Delta;E}}_0^i\left(k\right)=S_1{E}_0^i\left(k\right)\\ {\mathrm{\&Delta;}}^2{E}_0^i=S_1{\mathrm{\&Delta;E}}_0^i\left(k\right)={S_1}^2{E}_0^i\left(k\right)\end{array}\right.$

进一步可将性能指标转化为

$\begin{array}{c}\min J_i\left(k\right)=E_0^i{\left(k\right)}^T{K}_I^i{E}_0^i\left(k\right)+E_0^i{\left(k\right)}^T{S_1}^T{K}_p^i{S}_1{E}_0^i\left(k\right)\\ +E_0^i{\left(k\right)}^T{\left({S_1}^2\right)}^T{K}_d^i\left({S_1}^2\right)E_0^i\left(k\right)+{\mathrm{\&Delta;u}}_{i,M}{\left(k\right)}^T{R}_i{\mathrm{\&Delta;u}}_{i,M}\left(k\right)\\ =E_0^i{\left(k\right)}^T{Q}_i{E}_0^i\left(k\right)+{\mathrm{\&Delta;u}}_{i,M}{\left(k\right)}^T{R}_i{\mathrm{\&Delta;u}}_{i,M}\left(k\right)\end{array}$

其中，

2.6依据纳什最优的思想，以Δu_i,M(k)为控制变量，由求解最优控制律，形式如下：

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

其中，

2.7由步骤2.2到步骤2.6，得到k时刻第i个智能体的新一轮迭代最优解为：

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

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

${\mathrm{\&Delta;u}}_M^{l+1}\left(k\right)=D_1\left(\mathrm{\&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)={\&lsqb;{\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)\&rsqb;}^T$

${\mathrm{\&Delta;u}}_M^l\left(k\right)={\&lsqb;{\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)\&rsqb;}^T$

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

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

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