CN106444362A - 一种废塑料裂解炉炉膛温度的分布式pid型预测函数控制方法 - Google Patents

Abstract

本发明公开了一种废塑料裂解炉炉膛温度的分布式PID型预测函数控制方法。本发明首先通过采集控制系统的实时阶跃响应数据建立过程对象的阶跃响应模型向量，再将大规模系统的在线优化问题转化为各个小规模子系统的优化求解问题，并把网络环境下的每个子系统看作一个智能体，同时各智能体之间通过网络通信实现信息共享。然后通过引入PID算子对DPFC方法的性能指标进行改进，并设计各智能体的PID型预测函数控制器。再将当前时刻所得即时控制律作用于每个智能体，并将时域滚动至下一时刻，最后通过依次迭代循环完成整个大规模系统的优化任务。本发明有效弥补了传统 DPFC方法的不足，并提高了控制参数设计的自由度。

Description

一种废塑料裂解炉炉膛温度的分布式PID型预测函数控制 方法

技术领域

本发明属于自动化技术领域，涉及一种废塑料裂解炉炉膛温度的分布式PID型预测函数控制方法。

背景技术

随着计算机网络技术的发展，控制作为信息处理手段的一种已不仅仅局限于集中式的控制结构，而是更多为分布式的控制结构所取代，给传统的控制问题带来了诸多挑战。分布式预测函数控制(DPFC)作为预测控制在分布式控制结构中的典型应用，在有效解决复杂高维大规模系统在线优化控制问题的同时，通常存在对模型要求高、稳态性能差、不能快速抑制未知干扰等不足。因此对于模型阶次，环境扰动、非线性等方面存在较大不确定性的控制系统，常规的DPFC方法往往难以达到预期的控制效果。而传统的PID控制由于其控制结构简单、操作方便、鲁棒性强等优点，至今仍被广泛的应用于实际工业过程中。如果能在实际过程中将PID控制与DPFC方法相结合，将进一步推进预测控制在分布式控制结构中的发展，在保证系统良好控制性能的同时，有效增加了控制参数设计的灵活性。

发明内容

本发明的目的是为了克服常规DPFC方法在应用于在模型阶次、环境扰动、非线性等方面存在不确定性较大的控制系统中的不足之处，提出了一种废塑料裂解炉炉膛温度的分布式PID型预测函数控制方法。该方法在常规的DPFC方法中引入传统PID控制，在保证系统良好控制性能的同时有效提高了控制参数设计的自由度。本发明方法首先通过采集控制系统的实时阶跃响应数据建立过程对象的阶跃响应模型向量，再将大规模系统的在线优化问题转化为各个小规模子系统的优化求解问题，并把网络环境下的每个子系统看作一个智能体，同时各智能体之间通过网络通信实现信息共享，保证了系统整体的控制品质。然后通过引入PID算子对DPFC方法的性能指标进行改进，并依据纳什最优思想设计各智能体的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

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

步骤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个智能体的控制量：

其中u_i(k+j)表示第i个智能体在k+j时刻的控制量，E表示为基函数的个数，μ_i,n表示基函数的线性加权系数，f_i,kn(n＝1,2,…,E)表示第i个智能体的基函数，f_i,kn(j)表示第i个智能体的基函数f_i,kn在t＝jT_s时刻的值，T_s表示采样周期。

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

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

△u_i(k+j)＝u_i(k+j)-u_i(k+j-1)

结合步骤2.2进一步推导可得

其中，

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

F_i,0E＝[f_i,k1(-1),f_i,k2(-1),…,f_i,kE(-1)],μ_i(k)＝[μ_i,1(-1),μ_i,2(-1),…,μ_i,E(-1)]^T

G_ij,0＝A_ij,0F_j,0E,u_i,0＝u_i(k-2)

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_i,0为第i个智能体k-2时刻的控制输入；

然后可以得到第i个智能体在k时刻的模型预测误差值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时刻对k,k+1,…,k+L-1时刻预测模型的修正值，h为误差补偿的权向量，α为误差校正系数；

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

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

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

2.4获取第i个智能体在M个连续的控制增量作用下的预测输出值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

A0_i,P0＝[A0_i1,P0,A0_i2,P0,…,A0_iN,P0]^T,A0_ij,P0＝[a_ij(1),a_ij(2),…,a_ij(M),…,a_ij(P)]^T

u_P0＝[u_1,P0,u_2,P0,…,u_N,P0]^T,μ_i(k)＝[μ_i,1(k),μ_i,2(k),…,μ_i,E(k)]^T

G_ij(k)＝A_ijF_j,E

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.5选取第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个智能体的比例系数矩阵、积分系数矩阵、微分系数矩阵，为第i个智能体的控制加权系数矩阵，ω_i(k+ε)为第i个智能体在k+ε时刻的参考轨迹，y_i(k)为k时刻第i个智能体的实际输出，c(k)为k时刻第i个智能体的期望输出，λ为参考轨迹的柔化因子。

2.6对步骤2.5中的性能指标进行转换，形式如下：

进一步得到

同理可得

其中

引入矩阵

进而有

进一步可将性能指标变换为

其中，

2.7依据纳什最优的概念，对性能指标求极值，可以得到形式如下的纳什最优解：

其中，

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

进一步得到k时刻整个系统的纳什最优解：

μ^l+1(k)＝D₁[w(k)-Y_P0(k)+A0_i,P0u_P0]+D₀μ^l(k)

其中：

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

2.9由第i个智能体k时刻的纳什最优解得到最优控制量u_i(k)，并将其作用于第i个智能体；

2.10在下一时刻，重复步骤2.2到2.9继续求解第i个智能体的纳什最优解进而得到整个大规模系统的最优解μ^*(k+1)，并依次循环。

本发明提出了一种废塑料裂解炉炉膛温度的分布式PID型预测函数控制方法。该方法通过采集系统的实时阶跃响应数据建立被控对象的阶跃响应模型向量，并依据纳什优化的思想设计了一种改进的分布式预测函数控制器，在保证系统整体控制品质的同时，有效弥补了传统DPFC方法的不足，并提高了控制参数设计的自由度。

具体实施方式

以废塑料裂解炉炉膛温度控制为例：

废塑料裂解炉炉膛温度控制系统是一个典型的多变量含滞后的耦合过程，调节手段采用燃烧火嘴开度。

步骤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

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

步骤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个炉膛的控制量：

其中u_i(k+j)表示第i个炉膛在k+j时刻的燃烧火嘴开度，E表示为基函数的个数，μ_i,n表示基函数的线性加权系数，f_i,kn(n＝1,2,…,E)表示第i个炉膛的基函数，f_i,kn(j)表示第i个炉膛的基函数f_i,kn在t＝jT_s时刻的值，T_s表示采样周期。

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

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

△u_i(k+j)＝u_i(k+j)-u_i(k+j-1)

结合步骤2.2进一步推导可得

其中，

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

F_i,0E＝[f_i,k1(-1),f_i,k2(-1),…,f_i,kE(-1)],μ_i(k)＝[μ_i,1(-1),μ_i,2(-1),…,μ_i,E(-1)]^T

G_ij,0＝A_ij,0F_j,0E,u_i,0＝u_i(k-2)

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_i,0为第i个炉膛k-2时刻的控制输入；

然后可以得到第i个炉膛在k时刻的模型预测误差值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时刻对k,k+1,…,k+L-1时刻预测模型的修正值，h为误差补偿的权向量，α为误差校正系数；

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

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

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

2.4获取第i个炉膛在M个连续的控制增量作用下的预测输出值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

A0_i,P0＝[A0_i1,P0,A0_i2,P0,…,A0_iN,P0]^T,A0_ij,P0＝[a_ij(1),a_ij(2),…,a_ij(M),…,a_ij(P)]^T

u_P0＝[u_1,P0,u_2,P0,…,u_N,P0]^T,μ_i(k)＝[μ_i,1(k),μ_i,2(k),…,μ_i,E(k)]^T

G_ij(k)＝A_ijF_j,E

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.5选取第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个炉膛的比例系数矩阵、积分系数矩阵、微分系数矩阵，为第i个炉膛的控制加权系数矩阵，ω_i(k+ε)为第i个炉膛在k+ε时刻的参考轨迹，y_i(k)为k时刻第i个炉膛的实际输出，c(k)为k时刻第i个炉膛的期望输出，λ为参考轨迹的柔化因子。

2.6对步骤2.5中的性能指标进行转换，形式如下：

进一步得到

同理可得

其中

引入矩阵

进而有

进一步可将性能指标变换为

其中，

2.7依据纳什最优的概念，对性能指标求极值，可以得到形式如下的纳什最优解：

其中，

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

进一步得到k时刻整个系统的纳什最优解：

μ^l+1(k)＝D₁[w(k)-Y_P0(k)+A0_i,P0u_P0]+D₀μ^l(k)

其中：

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

2.9由第i个炉膛k时刻的纳什最优解得到最优控制量u_i(k)，并将其作用于第i个炉膛；

2.10在下一时刻，重复步骤2.2到2.9继续求解第i个炉膛的纳什最优解进而得到整个分布式系统的最优解μ^*(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

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

步骤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个智能体的控制量：

$u_i(k+j)=\underset{n=1}{\overset{E}{\&Sigma;}}{\mathrm{\&mu;}}_{i,n}\left(k\right)\&times;f_{i,kn}\left(j\right),j=0,1,...,P-1$

其中u_i(k+j)表示第i个智能体在k+j时刻的控制量，E表示为基函数的个数，μ_i,n表示基函数的线性加权系数，f_i,kn(n＝1,2,…,E)表示第i个智能体的基函数，f_i,kn(j)表示第i个智能体的基函数f_i,kn在t＝jT_s时刻的值，T_s表示采样周期；

2.3获取第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)$

△u_i(k+j)＝u_i(k+j)-u_i(k+j-1)

结合步骤2.2进一步推导可得

$y_{i,P}(k-1)=y_{i,0}(k-1)-A_{ii,0}{u}_{i,0}+G_{ii,0}(k-1){\mathrm{\&mu;}}_i(k-1)+\underset{j=1,j\&NotEqual;i}{\overset{n}{\&Sigma;}}{G}_{ij,0}(k-1){\mathrm{\&mu;}}_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

F_i,0E＝[f_i,k1(-1),f_i,k2(-1),…,f_i,kE(-1)],μ_i(k)＝[μ_i,1(-1),μ_i,2(-1),…,μ_i,E(-1)]^T

G_ij,0＝A_ij,0F_j,0E,u_i,0＝u_i(k-2)

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_i,0为第i个智能体k-2时刻的控制输入；

然后得到第i个智能体在k时刻的模型预测误差值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时刻对k,k+1,…,k+L-1时刻预测模型的修正值，h为误差补偿的权向量，α为误差校正系数；

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

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

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

2.4获取第i个智能体在M个连续的控制增量作用下的预测输出值y_i,PM，具体是：

$y_{i,PM}\left(k\right)=y_{i,P0}\left(k\right)-A{0}_{i,P0}{u}_{P0}+G_{ii}\left(k\right){\mathrm{\&mu;}}_i\left(k\right)+\underset{j=1,j\&NotEqual;i}{\overset{n}{\&Sigma;}}{G}_{ij}\left(k\right){\mathrm{\&mu;}}_j\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

A0_i,P0＝[A0_i1,P0,A0_i2,P0,…,A0_iN,P0]^T,A0_ij,P0＝[a_ij(1),a_ij(2),…,a_ij(M),…,a_ij(P)]^T

u_P0＝[u_1,P0,u_2,P0,…,u_N,P0]^T,μ_i(k)＝[μ_i,1(k),μ_i,2(k),…,μ_i,E(k)]^T

G_ij(k)＝A_ijF_j,E

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.5选取第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个智能体的比例系数矩阵、积分系数矩阵、微分系数矩阵，为第i个智能体的控制加权系数矩阵，ω_i(k+ε)为第i个智能体在k+ε时刻的参考轨迹，y_i(k)为k时刻第i个智能体的实际输出，c(k)为k时刻第i个智能体的期望输出，λ为参考轨迹的柔化因子；

2.6对步骤2.5中的性能指标进行转换，形式如下：

$\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={\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)\&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\left(k\right)=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.7依据纳什最优的概念，对性能指标求极值，得到形式如下的纳什最优解：

${\mathrm{\&mu;}}_i^{*}\left(k\right)=D_{ii}\&lsqb;w_i\left(k\right)-y_{i,P0}\left(k\right)+A{0}_{i,P0}{u}_{P0}-\underset{j=1,j\&NotEqual;i}{\overset{n}{\&Sigma;}}{G}_{ij}{\mathrm{\&mu;}}_j^{*}\left(k\right)\&rsqb;$

其中，

2.8重复步骤2.2至步骤2.7，得到第i个智能体在k时刻的新一轮迭代最优解为：

${\mathrm{\&mu;}}_i^{l+1}\left(k\right)=D_{ii}\&lsqb;w_i\left(k\right)-y_{i,P0}\left(k\right)+A{0}_{i,P0}{u}_{P0}-\underset{j=1,j\&NotEqual;i}{\overset{n}{\&Sigma;}}{G}_{ij}{\mathrm{\&mu;}}_j^l\left(k\right)\&rsqb;$

进一步得到k时刻整个系统的纳什最优解：

μ^l+1(k)＝D₁[w(k)-Y_P0(k)+A0_i,P0u_P0]+D₀μ^l(k)

其中：

${\mathrm{\&mu;}}^{l+1}\left(k\right)={\&lsqb;{\mathrm{\&mu;}}_1^{l+1}\left(k\right),{\mathrm{\&mu;}}_2^{l+1}\left(k\right),...,{\mathrm{\&mu;}}_n^{l+1}\left(k\right)\&rsqb;}^T$

${\mathrm{\&mu;}}^l\left(k\right)={\&lsqb;{\mathrm{\&mu;}}_1^l\left(k\right),{\mathrm{\&mu;}}_2^l\left(k\right),...,{\mathrm{\&mu;}}_n^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.9由第i个智能体k时刻的纳什最优解得到最优控制量u_i(k)，并将其作用于第i个智能体；

2.10在下一时刻，重复步骤2.2到2.9继续求解第i个智能体的纳什最优解进而得到整个大规模系统的最优解μ^*(k+1)，并依次循环。