供热管网水力仿真模型辨识修正方法及系统、操作方法与流程

本发明属于供热领域，具体涉及一种供热管网水力仿真模型辨识修正方法及系统、操作方法。

背景技术：

建立符合实际供热管网运行特性的水力仿真模型是实现热网智慧化管理的重要基础。随着供热规模的日益扩大和网络结构的日趋复杂，建立严格精确的供热管网的水力仿真模型变得愈加困难。通常管网的连接结构、节点标高等参数较为容易得到，而管道阻力特性受多种物理参数影响，且随时间变化，通过严格机理计算获得的理论阻力特性往往与实际特性存在偏差，导致建立的管网水力仿真模型不能很好地模拟实际运行工况。为此，需要充分利用实际运行历史数据，采用辨识的方法修正管段阻力特性。目前大多数供热管网在热力站处有供、回水压力测点和流量测点，而其它管段和节点处并不存在测点，因而难以通过直接计算的方法获得整个管网管段的阻力系数。

技术实现要素：

本发明的目的是提供一种供热管网水力仿真模型辨识修正方法及系统，以解决供热管网水力仿真模型不能准确描述实际管网运行的技术问题。

为了解决上述技术问题，本发明提供了一种供热管网水力仿真模型辨识修正方法，包括如下步骤：

步骤S1，选取供热系统在稳态工况下的多组实测数据；

步骤S2，计算各管段的阻力系数理论值；以及

步骤S3，建立供热管网水力仿真模型，利用粒子群算法结合多组实测数据辨识计算各阻力系数修正量，并获得修正过的阻力系数。

进一步，所述步骤S1中选取供热系统在稳态工况下的多组实测数据，即

识别供热管网系统的多个稳态运行工况，选取多组稳态运行工况的实测数据，其方法包括：

步骤S11，工况数据预处理；

步骤S12，求各热力站的稳态工况时间区间；

步骤S13，求整个供热管网系统的稳态工况时间区间；以及

步骤S14，选取所述多组实测数据。

进一步，所述步骤S2中计算各管段的阻力系数理论值的方法包括：

通过如下机理公式计算出各管段的阻力系数理论值；即

δ_i为管段i的阻力系数(Pa/(t/h)²)；K为管壁的当量绝对粗糙度；d_i为管道内径，单位为m；l_i为管段长度，单位为m；l_di为管段局部阻力当量长度，单位为m；ρ为管内热水的平均密度，单位为kg/m³；D为管段数量。

进一步，所述步骤S3中建立供热管网水力仿真模型的方法包括：

将供热管网转换为由节点和区段构成的有向图模型，即管网图；其中

节点表示存在流量进出的点，用集合V表示，V＝{V₁,V₂,…,V_n}，式中，n为管网中的节点个数；

区段表示节点间的连接管段，用集合E表示，E＝{E₁,E₂,…,E_m}，式中D为区段数即管段数量；

转换为有向图，则表示为G＝<V,E>；

根据网络图论得到供热管网的关联矩阵A和基本回路矩阵B，其中A为n×D阶矩阵，B为s×D阶矩阵，s为基本回路个数s＝D-n+1；

构建供热管网水力仿真模型，即

AG^T＝0；

BΔH^T＝0；

式中，G为记录管网图中各管段内体积流量的行向量，即

G＝[G₁,G₂,…,G_D]；

ΔH为记录管网图中各管段总阻力损失的列向量，即

ΔH＝ζ*|G|*G+Z-H_b；

式中，ζ为管网中各管段的阻力系数矩阵；

ζ＝diag{δ₁+Δδ₁,δ₂+Δδ₂,…,δ_D+Δδ_D}；

其中δ₁、δ₂、……、δ_D为阻力系数理论值，Δδ₁、Δδ₂、……、Δδ_D为阻力系数修正量，即需要辨识的参数；D表示待辨识的阻力系数修正量的维数；

Z为管段两节点位能差的列向量Z＝[Z₁,Z₂,…,Z_D]^T；

H_b为热网中的水泵扬程列向量H_b＝[H_b1,H_b2,…,H_bD]^T；

在供热管网水力仿真模型应用于辨识时，将其抽象表示为如下形式：

Y＝Func(X,δ)

式中X表示输入的多组实测数据；Y表示输出向量；δ表示待辨识管段的阻力系数向量,δ＝[δ₁+Δδ₁,δ₂+Δδ₂,…,δ_D+Δδ_D]。

进一步，所述步骤S3中利用粒子群算法结合多组实测数据辨识计算各阻力系数修正量，即

采用粒子群算法在每一代产生不同的阻力系数修正量，并通过供热管网水力仿真模型中不同阻力系数修正量向量在多个工况下管网的流动状态与相应实测数据进行对比，辨识出最优的阻力系数修正量。

进一步，采用粒子群算法结合多组实测数据辨识计算各阻力系数修正量向量的方法包括如下步骤：

步骤S31，粒子群算法参数设置，即

微粒个数N、微粒维数D，以及决策变量上界U＝[Δδ_1u,Δδ_2u,…,Δδ_Du]和下界L＝[Δδ_1l,Δδ_2l,…,Δδ_Dl]；

步骤S32，初始化，即

读取阻力系数理论值，生成N个微粒，即N组阻力系数修正量向量；

步骤S33，建立用于阻力系数修正量辨识的目标函数，以计算各微粒的目标函数值；

所述目标函数如下：

当阻力系数修正量Δδ₁,Δδ₂,…,Δδ_D作为决策变量给定值后，通过供热管网水力仿真模型计算得到相应输出数据；

对每一个微粒k∈{1,2,…,N}，将其与阻力系数理论值向量[δ₁,δ₂,…,δ_D]相加，得到该微粒代表的阻力系数向量将δ_k输入供热管网水力仿真模型，并依次利用第1、2，……，第N_C组实测数据对其进行水力计算，得到相应工况下各热力站的压力计算值；

式中，

PC_ij：稳态工况j下由仿真模型计算得到的热力站i的压力计算值；

PO_ij：稳态工况j下热力站i对应的压力测量值；

N_P：压力测点数量；

步骤S34，微粒飞行，产生下一代微粒；

步骤S35，达到设定的运行代数，算法终止，否则回到步骤S33继续执行。

又一方面，本发明还提供了一种供热管网水力仿真模型辨识修正系统，包括：

粒子群算法单元，产生每一代的阻力系数修正量对应的微粒，并将微粒写入数据库，并且读取各微粒的目标函数值；

目标函数计算单元，建立用于阻力系数修正量辨识的目标函数，驱动多个供热管网水力仿真模型，即产生N个仿真模型进行并行的模拟计算，并从数据库获取微粒，通过仿真模型计算微粒的目标函数值并写回数据库。

数据库单元，存储阻力系数修正量向量的群体和目标函数值，以及目标函数的计算状态，并且还将稳态工况的实测数据以及每一代优化计算的结果也存入数据库中。

进一步，所述目标函数如下：

当阻力系数修正量Δδ₁,Δδ₂,…,Δδ_D作为决策变量给定值后，通过供热管网水力仿真模型计算得到相应输出数据；

对每一个微粒k∈{1,2,…,N}，将其与阻力系数理论值向量[δ₁,δ₂,…,δ_D]相加，得到该微粒代表的阻力系数向量将δ_k输入供热管网水力仿真模型，并依次利用第1、2，……，第N_C组实测数据对其进行水力计算，得到相应工况下各热力站的压力计算值；

式中，

PC_ij：稳态工况j下由仿真模型计算得到的热力站i的压力计算值；

PO_ij：稳态工况j下热力站i对应的压力测量值；

N_P：压力测点数量。

第三方面，本发明还提供了一种供热管网水力仿真模型辨识修正系统的操作方法，包括如下步骤：

步骤Sa，建立待辨识管网的拓扑结构模型，以生成管网各部件的唯一标识ID；

步骤Sb，获得工况实测数据和计算阻力系数理论值的准备工作；

步骤Sc，设置粒子群算法单元中的算法参数，启动该单元；

步骤Sd，启动目标函数计算单元。

进一步，所述步骤Sb中所述准备工作包括：将待辨识管段的ID和阻力系数理论值对应地整理到配置文件中，供粒子群算法单元读取；以及

将筛选的多组实测数据写入数据库中，供目标函数计算单元查询，每一组工况数据包括工况编号、热源或热力站ID、热源或热力站对应的流量和供回水压力。

本发明的有益效果是，本发明的供热管网水力仿真模型辨识修正方法及系统可显著提高供热管网水力仿真模型仿真计算的精度，其采用粒子群算法在管段的阻力系数理论值的基础上辨识阻力系数修正量，计算时利用了多组工况运行数据，在观测点不足且无法通过解析求解的情况下能够得到相对准确的结果；采用并行架构进行粒子群算法阻力系数修正量的辨识，计算效率高。供热管网阻力系数辨识不仅能够应用于建立高精度仿真模型，也为管网阻塞提供了一种诊断方法。

附图说明

下面结合附图和实施例对本发明进一步说明。

图1是本发明一种供热管网水力仿真模型辨识修正方法的流程图；

图2是本发明供热管网水力仿真模型辨识修正系统的结构图；

图3为本发明实施例的供热管网系统。

具体实施方式

现在结合附图对本发明作进一步详细的说明。这些附图均为简化的示意图，仅以示意方式说明本发明的基本结构，因此其仅显示与本发明有关的构成。

实施例1

如图1所示，本发明的一种供热管网水力仿真模型辨识修正方法，包括如下步骤：

步骤S1，选取供热系统在稳态工况下的多组实测数据；

步骤S2，计算各管段的阻力系数理论值；以及

步骤S3，建立供热管网水力仿真模型，利用粒子群算法结合多组实测数据辨识计算各阻力系数修正量，并获得修正过的阻力系数。

进一步，所述步骤S1选取供热系统在稳态工况下的多组实测数据根据热网历史运行工况数据选取多组实测数据，热网历史运行工况数据为一个供暖期的实测数据包括但不限于各热力站一次侧进、出口压力和流量，也可以采集热源的供、回水压力和流量等相关参数。

本实施例1将以热力站为例，对供热管网水力仿真模型辨识修正方法进行展开说明。

具体地，步骤S1展开包括如下步骤：

步骤S11，工况数据预处理，即作出各个热力站的时间——流量曲线图(折线图)，通过分析将认为明显不合理的数据予以剔除，如传感器或传输故障导致的过大或过小的数据，然后以被剔除数据前后数据的平均数填充。

步骤S12，求各热力站的稳态工况时间区间，即以流量稳定为工况稳定的条件，流量波动的次数即可近似为工况变化的次数，设各个热力站的时间——流量曲线图中相邻拐点流量的相对变化率例如但不限于大于3％表示流量波动。根据曲线图得出流量波动次数最多的热力站s_k，定义为工况基准站，并划分出各热力站的稳态工况时间区间其中s_i表示热力站i，n_i表示热力站i的流量历史数据可划分的稳态工况时间区间个数，……表示各时间区间。

步骤S13，求整个供热管网系统的稳态工况时间区间，即以流量波动最多的热力站s_k为基准，求其它各热力站的稳态工况时间区间与的交集即为整个热网系统的稳态工况时间区间，其中n_r为交集的个数，即稳态工况数，n_r≤n_k。

步骤S14，选取所述多组实测数据，即从每个稳态工况时间区间中选取一组完整的工况数据，同一稳态工况时间区间的多条数据采取求平均的方式，最终得到N_C组稳态工况数据，即为多组实测数据。

进一步，所述步骤S2计算各管段的阻力系数理论值通过如下机理公式计算出各管段的阻力系数理论值；即

δ_i为管段i的阻力系数(Pa/(t/h)²)；K为管壁的当量绝对粗糙度，一般K＝0.0005m；d_i为管道内径，单位为m；l_i为管段长度，单位为m；l_di为管段局部阻力当量长度，单位为m；ρ为管内热水的平均密度，单位为kg/m³；D为管段数量。

进一步，所述步骤S3中建立供热管网水力仿真模型的方法包括：

将供热管网转换为由节点和区段构成的有向图模型，即管网图；其中

节点表示存在流量进出的点，用集合V表示，V＝{V₁,V₂,…,V_n}，式中，n为管网中的节点个数；

区段表示节点间的连接管段，用集合E表示，E＝{E₁,E₂,…,E_m}，式中D为区段数即管段数量；

转换为有向图，则表示为G＝<V,E>；

根据网络图论得到供热管网的关联矩阵A和基本回路矩阵B，其中A为n×D阶矩阵，B为s×D阶矩阵，s为基本回路个数s＝D-n+1；

构建供热管网水力仿真模型，即

AG^T＝0；

BΔH^T＝0；

式中，G为记录管网图中各管段内体积流量的行向量，即

G＝[G₁,G₂,…,G_D]；

ΔH为记录管网图中各管段总阻力损失的列向量，即

ΔH＝ζ*|G|*G+Z-H_b；

式中，ζ为管网中各管段的阻力系数矩阵；

ζ＝diag{δ₁+Δδ₁,δ₂+Δδ₂,…,δ_D+Δδ_D}；

其中δ₁、δ₂、……、δ_D为阻力系数理论值，Δδ₁、Δδ₂、……、Δδ_D为阻力系数修正量，即需要辨识的参数；D表示待辨识的阻力系数修正量的维数；

Z为管段两节点位能差的列向量Z＝[Z₁,Z₂,…,Z_D]^T；

H_b为热网中的水泵扬程列向量H_b＝[H_b1,H_b2,…,H_bD]^T；

在供热管网水力仿真模型应用于辨识时，将其抽象表示为如下形式：

Y＝Func(X,δ)

式中X表示输入的多组实测数据，从选取的稳态工况实测数据中获取，包括但不限于总回水压力，各热力站一次侧流量和最不利站的压差(进出口压力之差)；热源的回水压力；Y表示输出向量，包括到达各热力站一次侧的进出口压力；δ表示待辨识管段的阻力系数向量,

δ＝[δ₁+Δδ₁,δ₂+Δδ₂,…,δ_D+Δδ_D]。

进一步，所述步骤S3中利用粒子群算法结合多组实测数据辨识计算各阻力系数修正量，即采用粒子群算法在每一代产生不同的阻力系数修正量向量，并通过供热管网水力仿真模型中不同阻力系数修正量向量在多个工况下管网的流动状态与相应实测数据进行对比，辨识出最优的阻力系数修正量。

具体地，采用粒子群算法结合多组实测数据辨识计算各阻力系数修正量的方法包括如下步骤：

步骤S31，粒子群算法参数设置，即

微粒个数N、微粒维数D，运行代数N_gen，以及决策变量上界U＝[Δδ_1u,Δδ_2u,…,Δδ_Du]和下界L＝[Δδ_1l,Δδ_2l,…,Δδ_Dl]；其中决策变量上下界一般设置为阻力系数理论值的+30％、-30％，可调整上述参数多次运行本算法以获得更好的辨识结果。

步骤S32，初始化，即

读取阻力系数理论值，在决策变量上下界之间随机生成N个微粒，即N组阻力系数修正量向量

步骤S33，建立用于阻力系数修正量辨识的目标函数，以计算各微粒的目标函数值；

所述目标函数如下：

当阻力系数修正量Δδ₁,Δδ₂,…,Δδ_D作为决策变量给定值后，通过供热管网水力仿真模型计算得到相应输出数据；

对每一个微粒k∈{1,2,…,N}，将其与阻力系数理论值向量[δ₁,δ₂,…,δ_D]相加，得到该微粒代表的阻力系数向量将δ_k输入供热管网水力仿真模型，并依次利用第1、2，……，第N_C组实测数据对其进行水力计算，得到相应工况下各热力站的压力计算值；

式中，

PC_ij：稳态工况j下由仿真模型计算得到的热力站i的压力计算值；

PO_ij：稳态工况j下热力站i对应的压力测量值；

N_P：压力测点数量；

步骤S34，微粒飞行，产生下一代微粒；

步骤S35，达到设定的运行代数N_gen，算法终止，否则回到步骤S33继续执行。

其中，运行代数N_gen可以根据需要设置相应数值，例如但不限于100、150、200、300等。

实施例2

在实施例1基础上，本实施例2提供了一种供热管网水力仿真模型辨识修正系统，包括：

粒子群算法单元，产生每一代的阻力系数修正量向量对应的微粒，并将微粒写入数据库，并且读取各微粒的目标函数值；

目标函数计算单元，建立用于阻力系数修正量辨识的目标函数，驱动多个供热管网水力仿真模型，即产生N个仿真模型进行并行的模拟计算，并从数据库获取微粒，通过仿真模型计算微粒的目标函数值并写回数据库。

数据库单元，存储阻力系数修正量向量的群体和目标函数值，以及目标函数的计算状态，并且还将稳态工况的实测数据以及每一代优化计算的结果也存入数据库中。

所述目标函数如下：

当阻力系数修正量Δδ₁,Δδ₂,…,Δδ_D作为决策变量给定值后，通过供热管网水力仿真模型计算得到相应输出数据；

对每一个微粒k∈{1,2,…,N}，将其与阻力系数理论值向量[δ₁,δ₂,…,δ_D]相加，得到该微粒代表的阻力系数向量将δ_k输入供热管网水力仿真模型，并依次利用第1、2，……，第N_C组实测数据对其进行水力计算，得到相应工况下各热力站的压力计算值；

式中，

PC_ij：稳态工况j下由仿真模型计算得到的热力站i的压力计算值；

PO_ij：稳态工况j下热力站i对应的压力测量值；

N_P：压力测点数量。

实施例3

在实施例1和2基础上，本实施例3提供了一种针对实施例2供热管网水力仿真模型辨识修正系统的操作方法，包括如下步骤：

步骤Sa，建立待辨识管网的拓扑结构模型，以生成管网各部件的唯一标识ID；

步骤Sb，获得工况实测数据和计算阻力系数理论值的准备工作；

步骤Sc，设置粒子群算法单元中的算法参数，启动该单元；

步骤Sd，启动目标函数计算单元。

所述步骤Sb中所述准备工作包括：将待辨识管段的ID和阻力系数理论值对应地整理到配置文件中，供粒子群算法单元读取；以及

将筛选的多组实测数据写入数据库中，供目标函数计算单元查询，每一组工况数据包括工况编号、热源或热力站ID、热源或热力站对应的流量和供回水压力。

实施例4

本实施例4为具体说明本发明采用粒子群算法计算模型中各管段的阻力系数修正量的方法。

如图3展示了一个供热管网系统，该系统包括1个热源，4个热力站，1个循环泵，供、回水管段共19根，分别用S1-S19表示。

表1为选取的5组稳态运行工况的实测数据。

表2为管网中管段的阻力系数理论值。

具体地，采用粒子群算法计算模型中各管段的阻力系数修正量的方法包括如下步骤：

步骤S31，粒子群算法参数设置：运行代数N_gen＝200，微粒个数N＝50、微粒维数即管段个数D＝19、决策变量上界U＝[2.913，2.913，0.5958，0.5958，1.3965，1.3965，2.7255，2.7255，2.6298，2.6298，0.1461，0.1461，3.0537，3.0537，2.2866，2.2866，0.4563，0.4563，2.9103]，

决策变量下界L＝-[2.913，2.913，0.5958，0.5958，1.3965，1.3965，2.7255，2.7255，2.6298，2.6298，0.1461，0.1461，3.0537，3.0537，2.2866，2.2866，0.4563，0.4563，2.9103]。

步骤S32，初始化。读取阻力系数理论值，在决策变量上界U和下界L之间随机生成N＝50个微粒，即50组供回水管道的阻力系数修正量向量。

步骤S33，计算微粒的目标函数值。

所建立的粒子群算法阻力系数修正量辨识的目标函数如下：

所述目标函数本身不是决策变量的显示表达式，当给定决策变量Δδ₁,Δδ₂,…,Δδ_D的值后，通过供热管网水力仿真模型计算得到压力数据。将生成的N＝50个微粒即50组阻力系数修正量向量与阻力系数理论值向量相加后分别输入管网水力仿真模型，然后依次利用第1，第2，……，第N_C＝5组工况数据对其进行水力计算，得到不同工况下热力站的压力计算值。

约束：通过调用供热管网水力仿真模型的水力计算模块可以自动满足质量、能量守恒约束。

式中，

Δδ₁,Δδ₂,…,Δδ_D：待辨识的阻力系数修正量，本实施例D＝19；

PC_ij：工况j下由模型计算得到的热力站i的压力计算值；

PO_ij：工况j下热力站i对应的压力测量值；

N_P：压力测点数量，本实施例4个热力站的供、回水压力测点数量为N_P＝4×2＝8；

步骤S34，微粒飞行，产生下一代微粒；

微粒飞行过程中微粒坐标不超出U和L限定的范围。

步骤S35，达到设定的运行代数N_gen＝200，算法终止，否则回到步骤S33继续执行。

通过辨识，得到优化的微粒即阻力系数修正量如表3所示，该微粒的目标函数值是0.14691。

表3为本发明实施例辨识计算的阻力系数修正量。

计算阻力系数理论值代表的微粒，即阻力系数修正量向量为[0,0,…,0]的目标函数值为0.15732，可见辨识后的阻力系数要优于阻力系数理论值，说明该组系数更能描述实际热网的运行状态。

以上述依据本发明的理想实施例为启示，通过上述的说明内容，相关工作人员完全可以在不偏离本项发明技术思想的范围内，进行多样的变更以及修改。本项发明的技术性范围并不局限于说明书上的内容，必须要根据权利要求范围来确定其技术性范围。