CN112817282B - 一种基于数字孪生系统的高炉炉顶压力控制回路参数优化方法 - Google Patents

Abstract

本发明公开了一种基于数字孪生系统的高炉炉顶压力控制回路参数优化方法。该方法步骤：首先利用三维重建技术构建高炉设备三维模型；随后根据高炉设备技术要求，设计相关动作引擎，在特定操作参数情况下，完成与显示实体相同的动作；设计Web页面，用于对相关控制器参数进行设置，并提供数据可视化功能。使用傅里叶变换提取振动数据的主要特征，进而将其与非振动数据重新整合以保障数据维数一致；采用Local‑CVA算法建立系统状态模型；设定Harris控制性能指标为优化目标，通过粒子群算法调整控制器参数数值，直到计算得到最优的控制性能指标；最终，计算得到的所有数据返回前端，三维模型进行相应动作，并且Web页面的曲线图及柱状图也根据返回数据进行实时显示。

Description

一种基于数字孪生系统的高炉炉顶压力控制回路参数优化 方法

技术领域

本发明涉及数字孪生与过程控制技术领域，具体涉及一种基于数字孪生系统的高炉炉顶压力控制回路参数优化方法。

背景技术

随着工业4.0的推进，数字化的概念已经进入了业界的每一个企业。所有企业都在就如何实现自身数字化，如何开发与设计更加具有数字化的产品不断做出尝试。

目前，对于企业的数字化转型，其中最炙手可热技术便是数字孪生。最近，由工信部中国电子技术标准化研究院牵头编写的2020年《数字孪生应用白皮书》在中国国际高新技术成果交易会中发布。其中，着重关注了智慧城市、智慧交通、智慧能源、智慧建筑、智能制造、智慧健康6大领域。其中，智能制造作为其中重要的组成部分收到了业界广泛关注。

现在，很多企业已经将数字孪生技术和自身企业发展相结合，做出了初步尝试。河钢和西门子合作，在钢材制造中引入数字化进行顶层设计。在传统钢铁工厂中，各个不同车间中不同的自动化和信息化系统，往往造成各自为政现象的产生，互相之间无法进行有效的信息共享与集成，进而形成了一个个信息孤岛。为此，西门子利用数字孪生技术实现的虚拟场景的工业流程验证和布局优化。中冶京诚工程技术有限公司也利用数字孪生技术实现了全流程虚拟钢铁工厂，实现了数字化的运维和管控。

然而，目前的大多数的数字孪生实践中只注重对实时数据的一个显示，这导致了我们仅仅实现了对数据的感知，而未对数据进行进一步处理。这不仅浪费了大量的工业数据资源，对耗时制作的三维模型的利用率也是极低的。

发明内容

针对上述现有技术中存在的缺陷，本发明提供了一种基于数字孪生系统的高炉炉顶压力控制回路参数优化方法。该方法将数字孪生技术与高炉炉顶压力控制，专注于高炉炼铁流程中各个工业流程的仿真控制和数字孪生体的构造。

为了实现上述发明目的，本发明所采用的技术方案为：

一种基于数字孪生系统的高炉炉顶压力控制回路参数优化方法，包括以下步骤：

步骤1：根据现实场景中的物理实体，利用Unity和3Dmax软件与三维重建技术，构建高炉设备三维模型；

步骤2：根据高炉设备技术要求，设计相关动作引擎，以完成在特定操作参数情况下，完成与高炉设备一致的动作；

步骤3；设计Web页面，使得现场操作人员可以进行相关控制器参数的设置，并提供了曲线图、柱状图等数据可视化的功能；

步骤4；使用傅里叶变换提取振动数据的主要特征，进而将其与非振动数据重新整合以保障数据维数一致；所述主要特征包括X、Y、Z轴有效值及其1倍频和2倍频的幅值；

步骤5：采用了Local-CVA算法建立系统状态模型，并计算Harris控制性能指标；Local-CVA算法采用局部线性嵌入方法对数据的非线性特性进行了进一步分析；

步骤6：设定Harris控制性能指标为优化目标，通过粒子群算法调整控制器参数数值，直到计算得到最优的控制性能指标；

步骤7：计算得到的所有数据返回前端，三维模型进行相应动作，并且Web页面的曲线图和柱状图也根据返回数据进行实时显示。

上述技术方案中，进一步地，所述步骤4中，使用傅里叶变换提取振动数据的主要特征，具体为：

步骤1)：获取训练数据：选取设备在正常运行状态下的震动数据作为训练数据，将传感器安装于设备表面，用于采集设备中需要进行监测的振动信号数据，每一次采样可以得到1×m的测量向量x_k,其中下标k为采样时刻，m为传感器数量，采样q次后得到采样矩阵

采样矩阵X包括X轴振动信号及其有效值、Y轴振动信号及其有效值、Z轴振动信号及其有效值等；

步骤2)：利用傅里叶变换提取X、Y、Z轴振动信号中的主要频域特征：将采集到的X轴振动信号、Y轴振动信号、Z轴振动信号根据时段分别使用傅里叶变换，得到各个时段的频率的幅值分布。随后，采用各个时间段中的1倍频和2倍频以表征各个时间段中的频率特征。

进一步地，所述步骤5具体包括如下步骤：

步骤1)：为了分析数据的动态特性，本发明首先通过数据的时序关系构建过去与将来数据矩阵：先对所有时刻的输入数据u和输出数据y进行均一化处理，对于给定k时刻，分别构建过去数据向量

和将来数据向量

最后采用所有的过去数据向量和将来数据向量分别组成过去数据矩阵

和将来数据矩阵

其中N＝n-f-p+1，p和f分别为过去与将来时滞参数，n为输入数据u的样本数。

步骤2)：为了进一步分析步骤1处理后数据的非线性特性，本发明使用局部线性嵌入方法对数据进行线性重构，其核心方法为将一段非线性数据分割为多个线性数据段的组合，为了最小化原始数据与线性重构后数据的误差，建立如下优化目标：

其中，l为邻近样本个数，

为线性重构权重，当所有样本的线性权重均计算完毕后，得到线性权重矩阵

进而，得到局部结构矩阵L：

L＝I-Ψ (公式2)

这里I为合适尺度的单位矩阵。基于以上方法，可以得到过去数据矩阵X_p的局部结构矩阵L_p和将来数据矩阵X_f的局部结构矩阵L_f。

步骤3)：分别计算过去和将来样本协方差：

和

过去与将来的交叉样本协方差：

为。其后，构建Hankel矩阵

步骤4)：对上述组成的Hankel矩阵Υ进行奇异值分解(SVD)：

Υ＝ГΛΔ^T (公式3)

其中，Г和Δ分别为左右奇异向量，Λ为奇异值矩阵。并根据奇异值矩阵中的相关系数将右奇异向量进行划分，将前c个奇异值所对应的右奇异向量记作Δ_c，则其状态变量子空间可以表示为：

步骤5)：通过如上的计算结构，可计算系统的状态空间方程：

x(t+1)＝Ax(t)+Bu(t)+e(t) (公式5)

y(t)＝Cx(t)+Du(t)+v(t) (公式6)

其中，A、B、C、D为系统矩阵，e、v为系统噪声，并且此处x(t)表示W中的第t列，同理u(t)、y(t)表示为U、Y的第t列，U和Y分别为输入数据u的矩阵和输出数据y的矩阵。

系统矩阵的估计

为：

其中Σ[.,.]表示样本协方差。并且，系统噪声的估计值

和

可计算得到：

随后可得到系统噪声方差

步骤6：采集监测时间段的输出值与设定值，计算实时输出值与设定值的偏差。随后计算该偏差序列的方差

最后根据下列公式求得该监测时段的Harris控制性能指标η：

进一步地，所述步骤6中，粒子群算法用一个D维的粒子来表示待优化问题的一个可能的解。每个粒子的状态改变需要根据四个信息进行，其中两个是每个粒子自身的基本属性：由各维数值构成的位置信息和用于改变位置的速度信息。其中每个粒子所在位置的好坏，交由一个待优化问题决定的适应度fitness函数进行评价。除此之外还有两个共享信息，一个是个体向种群的共享信息：自身到目前为止发现的位置中适应度最优的点，称为个体最优pbest，这个位置代表粒子个体的寻优经验。另一个是种群向个体的共享信息：整个种群到目前为止所有粒子个体经验中最优的位置，称为全局最优gbest，这个位置代表整个种群的寻优经验。这两个最优经验的共享机制，即体现了粒子之间的组织协同性。在获得这四个信息后，各粒子在每次迭代中将根据当前速度更新自己的位置，同时根据个体最优位置和种群最优位置信息更新下一次飞行的速度。

粒子群算法具体的执行过程采以伪代码方式表述为：

具体的数学表达和操作为，一个n个D维粒子组成的种群P＝{p₁,p₂,…,p_n}，每个粒子的位置为x_i＝(x_i,1,x_i,2,…,x_i,D)^T，速度为v_i＝(v_i,1,v_i,2,…,v_i,D)^T。粒子p_i的适应度为finessp_i，每个粒子的适应度即为Harris控制性能指标(，即把当前粒子下的Harris控制性能指标作为适应度)，自身的历史最优适应度为p_ibest，对应位置pbest＝(pbest_i,1,pbest_i,2,…,pbest_i,D)^T。通过对所有粒子的p_ibest进行比较排名统计出的种群最优适应度为gbest，对应的位置为(gbest_i,1,gbest_i,2,…,gbest_i,D)^T。从而可确定在最优Harris控制性能指标下的粒子位置，粒子位置即为最优情况下的控制器参数数值；粒子p_i在每次迭代中更新自身位置与速度的公式：

其中，rand₁、rand₂是介乎(0，1)之间的随机数。式中

被称为动量部分，表示新的搜索方向一定程度上延续了原有搜索方向，也可以视为粒子的飞行惯性。每个粒子的动量部分各不相同，体现了每个粒子搜索方向的随机性与差异性。

为粒子当前位置与自身个体最优位置之间方位矢量，因此可表示粒子的搜索过程中对自身经验的参考。

为粒子当前位置与全局历史最优位置之间方位矢量，表示粒子的搜索过程中对群体经验的参考，体现所有粒子搜索工作之间的协作性。系数C₁、C₂称为学习因子。C₁＝C₂表示个体经验和群体经验同等重要。C₁＝0表示不考虑个体经验，只考虑群体经验。这样一开始迭代，所有粒子就快速聚集在当前最优个体附近，结果容易陷入局部极值。C₂＝0表示所有粒子各行其是，仅按照自身个体经验改变搜索方向。此时的种群内部失去协同性，收敛速度会极慢。因此通常C₁、C₂二者相等。同时过低的C₁、C₂会使粒子惯性作用过大，难以转向。而过高的C₁、C₂可能会导致粒子飞过目标点，或在目标点附近徘徊。

粒子的速度更新公式中，第一项惯性部分影响算法的全局搜索能力，第二第三项经验部分影响局部搜索能力。为了平衡全局和局部所搜能力，提高收敛速度同时避免陷入局部极值，通常在第一项前面乘上一个惯性权重因子

的数值一般随迭代次数线性递减：

其中k_max为最大迭代世代数，

为惯性权重最大、最小值，k为当次迭代次数。取值上一般

在迭代的初期，

数值较高，粒子搜索的补偿较大，各粒子可以较快的找到解所在的大致区域。随着迭代世代数的增大，

数值渐渐下降，所有粒子的速度放慢，有利于随小区域进行精细的搜索。

进一步地，所述步骤1中，构建高炉设备的三维模型具体为，均采用利用Unity和3Dmax软件完成，所有模型均来源于实际现场设备，在外观上做到和现场设备一致，并拥有透视功能，实现对设备的全面解析；此外，对于模型间的管道，内部流动的气体通过数值控制流速，并用不同颜色来区分热气与冷气；此外，对于高炉内部的燃烧过程，通过不同颜色反应燃烧温度变化，并用小液滴表示反应生成物。

本发明的有益效果在于：

本发明将数字孪生技术与高炉顶压控制方法相结合，专注凸显数字孪生技术在用户交互方面的优越性，可以更为清晰地显示设备中各个变量的变化。对于工业过程中高采样率的振动数据，本发明使用了傅里叶变换提取振动数据的主要特征，并将主要特征与非振动数据对齐，为后续的数据建模奠定了基础。在系统建模中，本发明全面的分析了数据的动态性和非线性特性，采用了Local-CVA方法提高系统的模型精度。最后以Harris控制性能指标为目标，粒子群方法作为核心的参数寻优算法来解决现场控制器参数不佳的普遍问题。本发明采用的Local-CVA在传统CVA算法上，采用了局部线性嵌入的方法进一步分析了数据中的非线性特性，并将其应用于系统的状态空间方程构建中，能够有效提高所构建系统精度。

此外，本发明在人机交互方面，还可以从前端页面输入当前控制器参数，从而将输入的参数和后端寻优算法进行了关联。并且，能根据搜寻得到的最优参数进行模拟仿真，进一步将仿真数据传输到前端，在二维图表中实现对各个变量数值的显示，同时充分利用了数字孪生中所建立的大量三维模型，实现了在不同的数值情况下对三维模型的精准控制。

附图说明

图1为本发明总体方案实现流程图。

具体实施方式

以下结合附图和具体实施例对本发明做进一步阐述。

本实施例公开一种基于数字孪生系统的高炉炉顶压力控制回路参数优化方法，主要是面向各个炼钢厂中使用。在本方案中，需要收集高炉炼铁过程中高炉的情况及其内部运行过程中的主要产物外观和状态进行三维模型构建，需要根据实际情况进行状态空间模型的构建、控制性能指标的计算和搜寻最优控制器参数的算法构造，并且需要符合现场使用习惯的前端页面等。

总体方案实现流程图，如图1所示。

一种基于数字孪生系统的高炉炉顶压力控制回路参数优化方法，包括以下步骤：

步骤1：根据现实场景中的物理实体，利用Unity和3Dmax软件与三维重建技术，构建高炉设备三维模型；

步骤2：根据高炉设备技术要求，设计相关动作引擎，以完成在特定操作参数情况下，完成与高炉设备一致的动作；

步骤3；Web页面也被设计，使得现场操作人员可以进行相关控制器参数的设置，并提供了曲线图、柱状图等数据可视化的功能；

步骤4；使用傅里叶变换提取振动数据的主要特征，进而将其与非振动数据重新整合以保障数据维数一致；所述主要特征包括X、Y、Z轴有效值及其1倍频和2倍频的幅值；

步骤5：采用了Local-CVA算法建立系统状态模型，并计算Harris控制性能指标；Local-CVA算法采用局部线性嵌入方法对数据的非线性特性进行了进一步分析；

步骤6：设定Harris控制性能指标为优化目标，通过粒子群算法调整控制器参数数值，直到计算得到最优的控制性能指标；

步骤7：计算得到的所有数据返回前端，三维模型进行相应动作，并且Web页面的曲线图和柱状图也根据返回数据进行实时显示。

上述技术方案中，进一步地，所述步骤4中，使用傅里叶变换提取振动数据的主要特征，具体为：

步骤1)：获取训练数据：选取设备在正常运行状态下的震动数据作为训练数据，将传感器安装于设备表面，用于采集设备中需要进行监测的振动信号数据，每一次采样可以得到1×m的测量向量x_k,其中下标k为采样时刻，m为传感器数量，采样q次后得到采样矩阵

采样矩阵X包括X轴振动信号及其有效值、Y轴振动信号及其有效值、Z轴振动信号及其有效值等；

步骤2)：利用傅里叶变换提取X、Y、Z轴振动信号中的主要频域特征：将采集到的X轴振动信号、Y轴振动信号、Z轴振动信号根据时段分别使用傅里叶变换，得到各个时段的频率的幅值分布。随后，采用各个时间段中的1倍频和2倍频以表征各个时间段中的频率特征。

进一步地，所述步骤5具体包括如下步骤：

步骤1)：为了分析数据的动态特性，本发明首先通过数据的时序关系构建过去与将来数据矩阵：先对所有时刻的输入数据u和输出数据y进行均一化处理，对于给定k时刻，分别构建过去数据向量

和将来数据向量

最后采用所有的过去数据向量和将来数据向量分别组成过去数据矩阵

和将来数据矩阵

其中N＝n-f-p+1，p和f分别为过去与将来时滞参数，n为输入数据u的样本数。

步骤2)：为了进一步分析步骤1处理后数据的非线性特性，本发明使用局部线性嵌入方法对数据进行线性重构，其核心方法为将一段非线性数据分割为多个线性数据段的组合，为了最小化原始数据与线性重构后数据的误差，建立如下优化目标：

其中，l为邻近样本个数，

为线性重构权重，当所有样本的线性权重均计算完毕后，得到线性权重矩阵

进而，得到局部结构矩阵L：

L＝I-Ψ (公式2)

这里I为合适尺度的单位矩阵。基于以上方法，可以得到过去数据矩阵X_p的局部结构矩阵L_p和将来数据矩阵X_f的局部结构矩阵L_f。

步骤3)：分别计算过去和将来样本协方差：

和

过去与将来的交叉样本协方差：

为。其后，构建Hankel矩阵

步骤4)：对上述组成的Hankel矩阵Υ进行奇异值分解(SVD)：

Υ＝ГΛΔ^T (公式3)

其中，Γ和Δ分别为左右奇异向量，Λ为奇异值矩阵。并根据奇异值矩阵中的相关系数将右奇异向量进行划分，将前c个奇异值所对应的右奇异向量记作Δ_c，则其状态变量子空间可以表示为：

步骤5)：通过如上的计算结构，可计算系统的状态空间方程：

x(t+1)＝Ax(t)+Bu(t)+e(t) (公式5)

y(t)＝Cx(t)+Du(t)+v(t) (公式6)

其中，A、B、C、D为系统矩阵，e、v为系统噪声，并且此处x(t)表示W中的第t列，同理u(t)、y(t)表示为U、Y的第t列，U和Y分别为输入数据u的矩阵和输出数据y的矩阵。

系统矩阵的估计

为：

其中Σ[.,.]表示样本协方差。并且，系统噪声的估计值

和

可计算得到：

随后可得到系统噪声方差

步骤6：采集监测时间段的输出值与设定值，计算实时输出值与设定值的偏差。随后计算该偏差序列的方差

最后根据下列公式求得该监测时段的Harris控制性能指标η：

进一步地，所述步骤6中，粒子群算法用一个D维的粒子来表示待优化问题的一个可能的解。每个粒子的状态改变需要根据四个信息进行，其中两个是每个粒子自身的基本属性：由各维数值构成的位置信息和用于改变位置的速度信息。其中每个粒子所在位置的好坏，交由一个待优化问题决定的适应度fitness函数进行评价。除此之外还有两个共享信息，一个是个体向种群的共享信息：自身到目前为止发现的位置中适应度最优的点，称为个体最优pbest，这个位置代表粒子个体的寻优经验。另一个是种群向个体的共享信息：整个种群到目前为止所有粒子个体经验中最优的位置，称为全局最优gbest，这个位置代表整个种群的寻优经验。这两个最优经验的共享机制，即体现了粒子之间的组织协同性。在获得这四个信息后，各粒子在每次迭代中将根据当前速度更新自己的位置，同时根据个体最优位置和种群最优位置信息更新下一次飞行的速度。

粒子群算法具体的执行过程采以伪代码方式表述为：

具体的数学表达和操作为，一个n个D维粒子组成的种群P＝{p₁,p₂,…,p_n}，每个粒子的位置为x_i＝(x_i,1,x_i,2,…,x_i,D)^T，速度为v_i＝(v_i,1,v_i,2,…,v_i,D)^T。粒子p_i的适应度为finessp_i，每个粒子的适应度即为Harris控制性能指标(，即把当前粒子下的Harris控制性能指标作为适应度)；自身的历史最优适应度为p_ibest，对应位置pbest＝(pbest_i,1,pbest_i,2,…,pbest_i,D)^T。通过对所有粒子的p_ibest进行比较排名统计出的种群最优适应度为gbest，对应的位置为(gbest_i,1,gbest_i,2,…,gbest_i,D)^T。从而可确定在最优Harris控制性能指标下的粒子位置，粒子位置即为最优情况下的控制器参数数值；粒子p_i在每次迭代中更新自身位置与速度的公式：

其中，rand₁、rand₂是介乎(0，1)之间的随机数。式中

被称为动量部分，表示新的搜索方向一定程度上延续了原有搜索方向，也可以视为粒子的飞行惯性。每个粒子的动量部分各不相同，体现了每个粒子搜索方向的随机性与差异性。

为粒子当前位置与自身个体最优位置之间方位矢量，因此可表示粒子的搜索过程中对自身经验的参考。

为粒子当前位置与全局历史最优位置之间方位矢量，表示粒子的搜索过程中对群体经验的参考，体现所有粒子搜索工作之间的协作性。系数C₁、C₂称为学习因子。C₁＝C₂表示个体经验和群体经验同等重要。C₁＝0表示不考虑个体经验，只考虑群体经验。这样一开始迭代，所有粒子就快速聚集在当前最优个体附近，结果容易陷入局部极值。C₂＝0表示所有粒子各行其是，仅按照自身个体经验改变搜索方向。此时的种群内部失去协同性，收敛速度会极慢。因此通常C₁、C₂二者相等。同时过低的C₁、C₂会使粒子惯性作用过大，难以转向。而过高的C₁、C₂可能会导致粒子飞过目标点，或在目标点附近徘徊。

粒子的速度更新公式中，第一项惯性部分影响算法的全局搜索能力，第二第三项经验部分影响局部搜索能力。为了平衡全局和局部所搜能力，提高收敛速度同时避免陷入局部极值，通常在第一项前面乘上一个惯性权重因子

的数值一般随迭代次数线性递减：

其中k_max为最大迭代世代数，

为惯性权重最大、最小值，k为当次迭代次数。取值上一般

在迭代的初期，

数值较高，粒子搜索的补偿较大，各粒子可以较快的找到解所在的大致区域。随着迭代世代数的增大，

数值渐渐下降，所有粒子的速度放慢，有利于随小区域进行精细的搜索。

所述步骤1中，构建高炉设备的三维模型具体为，采用利用Unity和3Dmax软件完成，所有模型均来源于实际现场设备，在外观上做到和现场设备一致，并拥有透视功能，实现对设备的全面解析；此外，对于模型间的管道，内部流动的气体通过数值控制流速，并用不同颜色来区分热气与冷气；此外，对于高炉内部的燃烧过程，通过不同颜色反应燃烧温度变化，并用小液滴表示反应生成物。

以上实施例仅用以说明本发明的技术方案而非限制，本领域普通技术人员对本发明的技术方案所做的其他修改或者等同替换，只要不脱离本发明技术方案的精神和范围，均应涵盖在本发明的权利要求范围中。

Claims (4)

1.一种基于数字孪生系统的高炉炉顶压力控制回路参数优化方法，其特征在于，包括以下步骤：

步骤1：利用Unity和3Dmax软件与三维重建技术，构建高炉设备三维模型；

步骤2：根据高炉设备技术要求，设计相关动作引擎，在特定操作参数情况下，完成与高炉设备一致的动作；

步骤3；设计Web页面，用于对相关控制器参数进行设置，并提供曲线图、柱状图数据可视化的功能；

步骤4；使用傅里叶变换提取振动数据的主要特征，进而将其与非振动数据重新整合以保障数据维数一致；所述主要特征包括X、Y、Z轴有效值及其1倍频和2倍频的幅值；

步骤5：采用了Local-CVA算法建立系统状态模型，并计算Harris控制性能指标；Local-CVA算法采用局部线性嵌入方法对数据的非线性特性进行了进一步分析；

步骤6：设定Harris控制性能指标为优化目标，通过粒子群算法调整控制器参数数值，直到计算得到最优的控制性能指标；

步骤7：计算得到的所有数据返回前端，三维模型进行相应动作，并且Web页面的曲线图和柱状图也根据返回数据进行实时显示；

所述步骤5具体为：

步骤1)：分析整合后的数据的动态特性，首先通过数据的时序关系构建过去与将来数据矩阵：先对所有时刻的输入数据u和输出数据y进行均一化处理，对于给定k时刻，分别构建过去数据向量

和将来数据向量

最后采用所有的过去数据向量和将来数据向量分别组成过去数据矩阵

和将来数据矩阵

其中N＝n-f-p+1，p和f分别为过去与将来的时滞参数，n为输入数据u的样本数；

步骤2)：进一步分析步骤1处理后数据的非线性特性，使用局部线性嵌入方法对数据进行线性重构，具体为将一段非线性数据分割为多个线性数据段的组合，为了最小化原始数据与线性重构后数据的误差，建立如下优化目标：

其中，l为邻近样本个数，

为线性重构权重，当所有样本的线性权重均计算完毕后，得到线性权重矩阵

进而，得到局部结构矩阵L：

L＝I-Ψ (公式2)

其中I为合适尺度的单位矩阵；基于以上方法，可以得到过去数据矩阵P的局部结构矩阵L_p和将来数据矩阵F的局部结构矩阵L_f；

步骤3)：分别计算过去和将来样本协方差：

和

过去与将来的交叉样本协方差：

其后，构建Hankel矩阵

步骤4)：对Hankel矩阵γ进行奇异值分解：

γ＝ΓΛΔ^T (公式3)

其中，Г和Δ分别为左右奇异向量，Λ为奇异值矩阵；并根据奇异值矩阵中的相关系数将右奇异向量进行划分，将前c个奇异值所对应的右奇异向量记作Δ_c，则其状态变量子空间可以表示为：

步骤5)：计算系统的状态空间方程：

x(t+1)＝Ax(t)+Bu(t)+e(t) (公式5)

y(t)＝Cx(t)+Du(t)+v(t) (公式6)

其中，A、B、C、D为系统矩阵，e、v为系统噪声，并且此处x(t)表示W中的第t列，同理u(t)、y(t)表示为U、Y的第t列，U和Y分别为输入数据u的矩阵和输出数据y的矩阵；

系统矩阵的估计

为：

其中Σ[.，.]表示样本协方差；系统噪声的估计值

和

为：

随后可得到系统噪声方差

步骤6)：采集监测时间段的输出值与设定值，计算实时输出值与设定值的偏差；随后计算该偏差序列的方差

最后根据下列公式求得该监测时段的Harris控制性能指标η：

2.根据权利要求1所述的基于数字孪生系统的高炉炉顶压力控制回路参数优化方法，所述步骤4中，使用傅里叶变换提取振动数据的主要特征，具体为：

步骤1)：获取训练数据：选取设备在正常运行状态下的振动数据作为训练数据，将传感器安装于设备表面，用于采集设备中需要进行监测的振动信号数据，每一次采样得到1×m的测量向量x_k，其中下标k为采样时刻，m为传感器数量，采样q次后得到采样矩阵

采样矩阵X包括X轴振动信号及其有效值、Y轴振动信号及其有效值、Z轴振动信号及其有效值；

步骤2)：利用傅里叶变换提取X、Y、Z轴振动信号中的主要频域特征：将采集到的X轴振动信号、Y轴振动信号、Z轴振动信号根据时段分别使用傅里叶变换，得到各个时段的频率的幅值分布；随后，采用各个时间段中的1倍频和2倍频以表征各个时间段中的频率特征。

3.根据权利要求1所述的基于数字孪生系统的高炉炉顶压力控制回路参数优化方法，其特征在于：所述的步骤6具体为，

一个由n个D维粒子组成的种群P＝{p₁，p₂，…，p_n}，每个粒子的位置为x_i＝(x_i，1，x_i，2，…，x_i，D)^T，速度为v_i＝(v_i，1，v_i，2，…，v_i，D)^T；粒子p_i的适应度为finessp_i，每个粒子的适应度即为Harris控制性能指标，自身的历史最优适应度为p_ibest，对应位置pbest＝(pbest_i，1，pbest_i，2，…，pbest_i，D)^T；通过对所有粒子的p_ibest进行比较排名统计出的种群最优适应度为gbest，对应的位置为(gbest_i，1，gbest_i，2，…，gbest_i，D)^T；从而可确定在最优Harris控制性能指标下的粒子位置，粒子位置即为最优情况下的控制器参数数值；

粒子p_i在每次迭代中更新自身位置与速度的公式为：

其中，rand₁、rand₂是介乎(0，1)之间的随机数；

为惯性权重因子，

被称为动量部分，表示新的搜索方向一定程度上延续了原有搜索方向；

为粒子当前位置与自身个体最优位置之间的方位矢量，

为粒子当前位置与全局历史最优位置之间方位矢量，系数C₁、C₂称为学习因子，C₁＝C₂表示个体经验和群体经验同等重要，

的取值随迭代次数线性递减：

其中k_max为最大迭代世代数，

为惯性权重最大、最小值，k为当次迭代次数，

4.根据权利要求1所述的基于数字孪生系统的高炉炉顶压力控制回路参数优化方法，其特征在于：

所述步骤1中，构建高炉设备的三维模型具体为，采用利用Unity和3Dmax软件完成，所有模型均来源于实际现场设备，在外观上做到和现场设备一致，并拥有透视功能，实现对设备的全面解析；此外，对于模型间的管道，内部流动的气体通过数值控制流速，并用不同颜色来区分热气与冷气；此外，对于高炉内部的燃烧过程，通过不同颜色反应燃烧温度变化，并用小液滴表示反应生成物。

Images