CN106483853A - 炼油加热炉炉膛压力的分数阶分布式动态矩阵控制方法 - Google Patents

Abstract

本发明公开了一种炼油加热炉炉膛压力的分数阶分布式动态矩阵控制方法。本发明首先通过采集多变量过程的输入输出数据建立该过程的分数阶模型，然后将其近似为整数阶高阶模型，并获取对象的阶跃响应模型向量，再将多变量过程的在线优化实施问题转化为各个小规模子系统的优化求解问题，把网络环境下的每个子系统看为一个智能体。通过在各智能体的性能指标中引入PID算子，设计各智能体的分数阶PID型动态矩阵控制器，再将当前时刻的即时控制律作用于每个智能体，并将时域滚动至下一时刻，重复上述优化过程，从而完成整个大规模系统的优化任务。本发明有效改善了系统的控制性能，并提高了控制器参数设计的自由度。

Description

炼油加热炉炉膛压力的分数阶分布式动态矩阵控制方法

技术领域

本发明属于自动化技术领域，涉及一种炼油加热炉炉膛压力的分数阶分布式动态矩阵控制方法。

背景技术

随着工业生产中对产品质量及安全操作要求的不断提高，实际过程中普遍存在的复杂高维大规模系统已无法用整数阶微分方程进行精确描述，而分数阶微分方程恰能更精确地描述此类系统的内在特性，因而研究一类分数阶多变量过程的预测控制在线实施问题便显得尤为重要。分布式动态矩阵控制(DDMC)作为预测控制在分布式控制结构中的典型应用，综合利用网络通信技术与控制理论，将一个复杂高维系统的在线优化问题分散到各个低维子系统中分布求解，有效降低了问题的规模与复杂性，能很好的控制存在多变量、强耦合的被控对象。然而针对一类分数阶的多变量过程，传统的整数阶DDMC方法往往不能获得理想的控制效果。如果能够对传统整数阶DDMC方法进行改进，并将其推广至分数阶系统中，将有效弥补整数阶DDMC方法在分数阶系统控制中的不足，并促进DDMC方法在分数阶系统中的应用。

发明内容

本发明目的是针对传统整数阶DDMC在控制分数阶的多变量过程中存在的不足之处，提出了一种炼油加热炉炉膛压力的分数阶分布式动态矩阵控制方法。该方法利用PID控制对常规DDMC方法进行改进，并将其推广至一类多入多出分数阶线性系统中，在保证系统控制品质的同时，增加了控制参数设计的灵活性。

本发明方法首先通过采集多变量过程的输入输出数据建立该过程的分数阶模型，然后采用Oustaloup近似方法将其近似为整数阶高阶模型，并基于高阶模型获取对象的阶跃响应模型向量，再将多变量过程的在线优化实施问题转化为各个小规模子系统的优化求解问题，把网络环境下的每个子系统看为一个智能体，同时各智能体间通过网络通信完成信息交换以保证系统整体性能。通过在各智能体的性能指标中引入PID算子，并依据纳什最优思想来设计各智能体的分数阶PID型动态矩阵控制器，再将当前时刻的即时控制律作用于每个智能体，并将时域滚动至下一时刻，重复上述优化过程，从而完成整个大规模系统的优化任务。

本发明的技术方案是通过数据采集、模型建立、预测机理、优化等手段，确立了一种炼油加热炉炉膛压力的分数阶分布式动态矩阵控制方法，利用该方法能很好的处理一类多入多出分数阶系统的控制问题，在改善系统控制性能的同时，有效提高了控制参数调节的自由度。

本发明方法的步骤包括：

步骤1.通过炼油加热炉炉膛压力对象的实时阶跃响应数据建立被控对象的分数阶动态矩阵控制的阶跃响应模型向量，具体方法是：

1.1根据分布式预测控制思想，将一个N输入N输出的大规模分数阶系统分散为N个智能体子系统。

1.2在稳态工况下，采集多变量过程对象各智能子系统的实时输入输出数据，利用该数据建立t时刻过程对象第j(1≤j≤N)个输入对第i(1≤i≤N)个输出的分数阶微分方程模型，形式如下：

其中，β_ij为第j个输入对第i个输出的分数阶微分阶次，为输出对应的系数，为输入对应的系数，τ_ij为第j个输入对第i个输出的滞后时间，y_i(t),u_j(t)分别为t时刻第i个智能体的输出和第j个智能体的输入。

1.3根据分数阶微积分定义，对步骤1.2中的模型进行拉氏变换，得到过程对象第j个输入对第i个输出的传递函数形式如下：

其中，s为拉普拉斯变换算子。

1.4根据Oustaloup近似方法得到微分算子s^β的近似表达式：

其中，β为分数阶微分阶次，且满足0<β<1，N₁为选定的近似阶次，w_b和w_h分别为选定的拟合频率的下限和上限。

1.5根据步骤1.4中的方法，将步骤1.3中的分数阶模型近似为整数阶高阶模型，给该高阶模型施加一个单位阶跃信号，并记录其阶跃响应曲线。

1.6将步骤1.5得到的阶跃响应曲线进行滤波处理，然后拟合成一条光滑曲线，记录光滑曲线上每个采样时刻对应的阶跃响应数据，第一个采样时刻为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用来表示系统的统一建模时域，则有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)，得到第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个智能体的阶跃响应数据建立的矩阵。

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

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

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

进一步得到第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)分别表示k时刻第i个智能体模型的修正值，h为误差补偿的权矩阵，α为误差校正系数。

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

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

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

2.3计算第i个智能体在k时刻加入各智能体的输入控制增量序列Δu_1,M(k),Δu_2,M(k),…,Δu_n,M(k)后的预测输出值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,M(k+1|k),y_i,M(k+2|k),…,y_i,M(k+P|k)分别表示第i个智能体在k时刻对k+1,k+2,…,k+P时刻的模型预测输出值,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)，形式如下：

minJ_i(k)＝(w_i(k)-y_i,PM(k))^TK_I ⁱ(w_i(k)-y_i,PM(k))+(Δw_i(k)-Δy_i,PM(k))^TKⁱ _p(Δw_i(k)-Δy_i,PM(k))+

(Δ²w_i(k)-Δ²y_i,PM(k))^TKⁱ _d(Δ²w_i(k)-Δ²y_i,PM(k))+Δu_i,M(k)^TR_iΔu_i,M(k)

其中，

w_i(k)＝[ω_i(k+1),ω_i(k+2),…,ω_i(k+P)]^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

Δ²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+ε)＝λ^εy_i(k)+(1-λ^ε)c(k)(ε＝1,2,…,P)

Δω_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个智能体给定期望输出的参考轨迹，y_i(k)为k时刻第i个智能体的过程实际输出，c(k)为k时刻第i个智能体的期望输出，λ为参考轨迹的柔化因子。

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

其中，

进一步得到

同理可得

其中，

引入矩阵

进而有

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

其中，

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

其中，D_ii＝(R_i+A_ii ^TQ_iA_ii)^-1A_ii ^TQ_i；

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将k时刻第i个智能体的纳什最优解的首项作为即时控制增量Δu_i(k)，进而得到实际控制量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算子，设计了一种分数阶系统的分布式PID型动态矩阵控制器，有效改善了系统的控制性能，并提高了控制器参数设计的自由度，同时进一步推广了DDMC方法在分数阶系统中的应用。

具体实施方式

以炼油加热炉炉膛压力过程控制为例：

炼油加热炉炉膛压力控制系统是一个典型的多变量过程对象，调节手段采用调节烟道挡板的阀门开度。

步骤1.通过炼油加热炉炉膛压力对象的实时阶跃响应数据建立被控对象的分数阶动态矩阵控制的阶跃响应模型向量，具体方法是：

1.1根据分布式预测控制思想，将一个N输入N输出的大规模分数阶系统分散为N个智能体子系统。

1.2在稳态工况下，采集炼油加热炉炉膛压力控制系统各炉膛的实时输入输出数据，利用该数据建立t时刻炉膛压力对象第j(1≤j≤N)个输入对第i(1≤i≤N)个输出的分数阶微分方程模型，形式如下：

其中，β_ij为第j个输入对第i个输出的分数阶微分阶次，为炉膛压力对应的系数，为阀门开度对应的系数，τ_ij为第j个输入对第i个输出的滞后时间，y_i(t),u_j(t)分别为t时刻第i个炉膛的压力和第j个炉膛的阀门开度。

1.3根据分数阶微积分定义，对步骤1.2中的模型进行拉氏变换，得到炉膛压力对象第j个输入对第i个输出的传递函数形式如下：

其中，s为拉普拉斯变换算子。

1.4根据Oustaloup近似方法得到微分算子s^β的近似表达式：

其中，β为分数阶微分阶次，且满足0<β<1，N₁为选定的近似阶次，w_b和w_h分别为选定的拟合频率的下限和上限。

1.5根据步骤1.4中的方法，将步骤1.3中的分数阶模型近似为整数阶高阶模型，给该高阶模型施加一个单位阶跃信号，并记录其阶跃响应曲线。

1.6将步骤1.5得到的阶跃响应曲线进行滤波处理，然后拟合成一条光滑曲线，记录光滑曲线上每个采样时刻对应的阶跃响应数据，第一个采样时刻为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用来表示系统的统一建模时域，则有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)，得到第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个炉膛的阶跃响应数据建立的矩阵。

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

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

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

进一步得到第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)分别表示k时刻第i个炉膛模型的修正值，h为误差补偿的权矩阵，α为误差校正系数。

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

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

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

2.3计算第i个炉膛在k时刻加入各炉膛的输入控制增量序列Δu_1,M(k),Δu_2,M(k),…,Δu_n,M(k)后的预测输出值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,M(k+1|k),y_i,M(k+2|k),…,y_i,M(k+P|k)分别表示第i个炉膛在k时刻对k+1,k+2,…,k+P时刻的模型预测输出值,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

Δ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+ε)＝λ^εy_i(k)+(1-λ^ε)c(k)(ε＝1,2,…,P)

Δω_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个炉膛给定期望输出的参考轨迹，y_i(k)为k时刻第i个炉膛的过程实际输出，c(k)为k时刻第i个炉膛的期望输出，λ为参考轨迹的柔化因子。

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

其中，

进一步得到

同理可得

其中，

引入矩阵

进而有

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

其中，

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

其中，D_ii＝(R_i+A_ii ^TQ_iA_ii)^-1A_ii ^TQ_i。

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将k时刻第i个炉膛的纳什最优解的首项作为即时控制增量Δu_i(k)，进而得到实际控制量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.炼油加热炉炉膛压力的分数阶分布式动态矩阵控制方法，其特征在于该方法包括以下步骤：

步骤1.通过炼油加热炉炉膛压力对象的实时阶跃响应数据建立被控对象的分数阶动态矩阵控制的阶跃响应模型向量，具体是：

1.1根据分布式预测控制思想，将一个N输入N输出的大规模分数阶系统分散为N个智能体子系统；

1.2在稳态工况下，采集多变量过程对象各智能子系统的实时输入输出数据，利用该数据建立t时刻过程对象第j(1≤j≤N)个输入对第i(1≤i≤N)个输出的分数阶微分方程模型，形式如下：

$b_1^{ij}{y_i}^{\left({\mathrm{\&beta;}}_{ij}\right)}\left(t\right)+b_0^{ij}{y}_i\left(t\right)=c_0^{ij}{u}_j(t-{\mathrm{\&tau;}}_{ij})$

其中，β_ij为第j个输入对第i个输出的分数阶微分阶次，为输出对应的系数，为输入对应的系数，τ_ij为第j个输入对第i个输出的滞后时间，y_i(t),u_j(t)分别为t时刻第i个智能体的输出和第j个智能体的输入；

1.3根据分数阶微积分定义，对步骤1.2中的模型进行拉氏变换，得到过程对象第j个输入对第i个输出的传递函数形式如下：

$G_{ij}\left(s\right)=\frac{c_0^{ij}{e}^{-{\mathrm{\&tau;}}_{ij}s}}{b_1^{ij}{s}^{{\mathrm{\&beta;}}_{ij}}+b_0^{ij}}$

其中，s为拉普拉斯变换算子；

1.4根据Oustaloup近似方法得到微分算子s^β的近似表达式：

$s^{\mathrm{\&beta;}}\&ap;K\underset{n_1=1}{\overset{N_1}{\mathrm{\&Pi;}}}\frac{s+w_{n_1}^{\&prime;}}{s+w_{n_1}}$

其中，β为分数阶微分阶次，0<β<1，N₁为选定的近似阶次，w_b和w_h分别为选定的拟合频率的下限和上限；

1.5根据步骤1.4中的Oustaloup近似方法，将步骤1.3中的分数阶模型近似为整数阶高阶模型，给该高阶模型施加一个单位阶跃信号，并记录其阶跃响应曲线；

1.6将步骤1.5得到的阶跃响应曲线进行滤波处理，然后拟合成一条光滑曲线，记录光滑曲线上每个采样时刻对应的阶跃响应数据，第一个采样时刻为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用来表示系统的统一建模时域，则有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)，

得到第i个智能体的模型预测值y_i,P(k-1)：

$y_{i,P}(k-1)=y_{i,0}(k-1)+A_{ii,0}{\mathrm{\&Delta;u}}_1(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个智能体的阶跃响应数据建立的矩阵；

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

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

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

进一步得到第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)分别表示k时刻第i个智能体模型的修正值，h为误差补偿的权矩阵，α为误差校正系数；

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

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

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

2.3计算第i个智能体在k时刻加入各智能体的输入控制增量序列△u_1,M(k),△u_2,M(k),…,△u_n,M(k)后的预测输出值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}{\mathrm{\&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,M(k+1|k),y_i,M(k+2|k),…,y_i,M(k+P|k)分别表示第i个智能体在k时刻对k+1,k+2,…,k+P时刻的模型预测输出值,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\left(k\right)-y_{i,PM}\left(k\right)\right)}^T{K}_I^i\left(w_i\left(k\right)-y_{i,PM}\left(k\right)\right)+{\left({\mathrm{\&Delta;w}}_i\left(k\right)-{\mathrm{\&Delta;y}}_{i,PM}\left(k\right)\right)}^T{K}_p^i\left({\mathrm{\&Delta;w}}_i\left(k\right)-{\mathrm{\&Delta;y}}_{i,PM}\left(k\right)\right)+\\ {\left({\mathrm{\&Delta;}}^2{w}_i\left(k\right)-{\mathrm{\&Delta;}}^2{y}_{i,PM}\left(k\right)\right)}^T{K}_d^i\left({\mathrm{\&Delta;}}^2{w}_i\left(k\right)-{\mathrm{\&Delta;}}^2{y}_{i,PM}\left(k\right)\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

△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+ε)＝λ^εy_i(k)+(1-λ^ε)c(k)(ε＝1,2,…,P)

△ω_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个智能体给定期望输出的参考轨迹，y_i(k)为k时刻第i个智能体的过程实际输出，c(k)为k时刻第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\left(k+1\right),{\mathrm{\&Delta;e}}_0^i\left(k+2\right),\mathrm{...},{\mathrm{\&Delta;e}}_0^i\left(k+P\right)\&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\left(k+1\right),{\mathrm{\&Delta;}}^2{e}_0^i\left(k+2\right),\mathrm{...},{\mathrm{\&Delta;}}^2{e}_0^i\left(k+P\right)\&rsqb;}^T$

进一步得到

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

同理可得

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

其中，

引入矩阵

进而有

$\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.6依据纳什最优的思想，以△u_i,M(k)为控制变量，由求解第i个智能体的最优控制律，形式如下：

${\mathrm{\&Delta;u}}_{i,M}^{*}\left(k\right)=D_{ii}\left(w_i\left(k\right)-y_{i,P0}\left(k\right)-\underset{j=1,j\&NotEqual;i}{\overset{n}{\mathrm{\&Sigma;}}}{A}_{ij}{u}_{j,M}^{*}\left(k\right)\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}{\mathrm{\&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(\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将k时刻第i个智能体的纳什最优解的首项作为即时控制增量△u_i(k)，进而得到实际控制量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)，并依次循环。