- 1 -
中国科技论文在线
基于最小二乘理论的供热系统阻力特性辨
识方法研究#
周志刚*
基金项目：高等学校博士学科点专项科研基金（20112302120046）
作者简介：周志刚（1978 ），男，副教授，城市集中供热系统参数辨识与故障诊断，智能热网
（哈尔滨工业大学市政环境工程学院，哈尔滨，150090） 5
摘要：基于最小二乘原理建立了供热系统阻力特性辨识的数学模型，将参数辨识问题转化非
线性优化求解。在改进遗传算法及适用于连续空间优化的蚂蚁算法的基础上，提出了一种新
型的混合遗传蚂蚁算法，计算性能获得了较大的提升。基于反映管网阻力特性系数变化提出
了流量与压力监测点优化布置方法。结合算例热网，说明了基于反映阻力特性变化的流量、10
压力监测点优化布置方法的实现过程，验证了了算法的优越性。
关键词：供热、供燃气、通风及空调工程；最小二乘；阻力特性；辨识
中图分类号：
Identification Resistance of District Heating System Based 15
on The Least Square Theory
Zhou Zhigang
(School of Municipal and Environmental Engineering,Harbin Institute of Technology,
Harbin,150090)
Abstract: The identification mathematic model of resistance coefficient of heat-supply networks is 20
established based on the least square theory. The identification problem is changed into non-linear
optimal problem by this model. Based on the improved genetic algorithm and the ant algorithm
that fits continuity space optimization, this new mixing algorithm improves computing efficiency.
The optimal arrangement method of flow measuring points and pressure measuring points, which
reflect the change of resistance coefficient of networks, is put forward based on the fuzzy 25
clustering theory. An example is shown to introduce the optimal method in detail, which aims at
identification of resistance coefficient of networks and reflects the flow and pressures that
character the resistance coefficient of networks. The superiority of the new algorithm is proved.
Key words: Heat,gas supply,ventilation and air conditioning engineering;The Least Square Theory;
Resistance; Identification 30
0 引言
城市供热系统管道阻力特性受管网结构与实际施工及使用情况的影响处于动态变化过
程中，准确获知其实际情况对于热网改扩建设计及运行调节具有重要的意义。常规通过理论
公式计算存在较大误差，逐一测量受管网规模及实施条件限制亦不现实，因此根据有限的管35
网监测数据，采用优化算法实现系统参数辨识是一种可行的方法。本文提出一种基于最小二
乘原理的供热系统阻力特性辨识方法，并对优化算法与测点布置进行了研究。通过算例对算
法的性能进行了验证。
1 供热系统阻力特性辨识数学模型
基于最小二乘原理，结合热网特点选择节点压力与管段流量作为监测数据，建立阻力特40
- 2 -
中国科技论文在线
性辨识数学模型。
2 2
, , 0 , , , 0 ,
1 1 1
,
( )
1
, ,
min ( ) ( ) ( )
0 1, 2
0 ( 1,2 )
.
- ( 1,2 )
k
yL x
pi t di t d i t gj t dj t d j t
t i j
m ij di
B
k
i k
i
ij ij m ij m ij ij
F s W p p W g g
g Q i N
p p k L
p s g g dh i N
S
min max ( 1,2 ) ijs S i N
(1)
式中 L ——用于阻力特性辨识的热网运行工况数；
s ——管网各管段的阻力特性系数， 2Pa/(t/h) ；
x ——热网中压力测点的数量； 45
y ——热网中流量测点的数量；
,di tp ——压力测点 i的辨识计算压力，Pa；
0 ,d i tp ——压力测点 i的实测压力，Pa；
,dj tg ——流量测点 j的辨识计算流量， t/h；
0 ,d j tg ——流量测点 j的实测流量， t/h； 50
,pi tW ——压力测点 i的权重因子；
,g j tW ——流量测点 j的权重因子；
,m ijg ——与节点 i相关联的管段质量流量代数和，t/h；
diQ ——节点 i的出流质量流量，t/h；
N ——节点个数； 55
j ——与节点 i关联的节点；
i ——节点下标；
( )
1
kB
k
i
i
p
——属于基本环路 k的管段阻力损失代数和，Pa；
kp ——环路平差精度；
L ——基本环路个数； 60
kB ——属于环路 k的管段数；
ijp ——与节点 i相关联管段的阻力损失，Pa；
ijs ——与节点 i相关联管段的阻力特性系数，
2Pa/(t/h) ；
,m ijg ——与节点 i相关联管段的质量流量，t/h；
ijdh ——与节点 i相关联管段的水泵扬程，Pa； 65
- 3 -
中国科技论文在线
minS ——管网管段阻力特性系数 S的下限向量；
maxS ——管网管段阻力特性系数 S的上限向量。
式(1)所示模型的物理意义为：在满足热网水力条件的约束下，在允许的范围内，通过
辨识管网管段阻力特性系数使管网中测压点实测压力与辨识计算压力的偏差以及测流管段
实测流量与辨识计算流量的偏差降至最小。理论上，所有因素均考虑的情况下，可调整到70
min ( ) 0F s ，但由于量测误差等原因，实际工程应用中即无必要也无可能达到这一要求，
因此最常用的方法是将误差控制在某一标准之下，即使min ( )F s 。这样，就将热网阻
力特性辨识问题就转化为非线性函数优化问题。
2 优化算法研究 75
供热系统阻力特性辨识的关键在于优化计算的时间与准确度。研究辨识算法的改进，减
少计算时间、提高计算精度对于该辨识方法能否在实际工程中应用具有十分重要的意义。本
文提出一种改进的优化算法，混合遗传蚁群算法。
遗传算法改进
为提高算法的计算效率，本文对标准遗传算法进行了必要的改进。引入父子竞争机制。80
在父子两代中选择最优的个体进入下一代，子代总是优于或者等于它们的父代，使进化总是
朝着最优的方向进行[1]。对交叉与变异概率进行自适应选择[2]。根据适应度自动改变交叉概
率 cP 和变异概率 mP ；如果适应度值高于平均适应度值，说明该个体性能优良，根据其适应
度值取相应的交叉率和变异率。
蚁群算法改进 85
热网阻力特性辨识这类连续空间的一般函数优化问题无法直接应用蚁群算法，通过改进
得到一种可用于连续函数优化的蚂蚁算法。
设连续空间优化函数为：
1 2
min max 2
min ( ), ( , , , )
. [ , ]
n
n
Z f X X x x x
X X X
(2)
首先确定蚁群初始位置，设蚁群规模为M ，将这M 个蚂蚁随机放置在优化空间 Z 等90
分后的各个区域，作为每个单蚁进行搜索的起点。其中单蚁初始所在各维子区间的长度 sal 为
max min
( ) ( )
( ) , 1,2, ,sa
X j X j
l j j n
M
(3)
完成一次循环后，蚂蚁 ( 1,2, , , )iX i M i best 将选择向找到最优解的蚂蚁 bestX 进
行全局转移，或是选择在原有位置的临域范围内进行随机搜索。蚂蚁 iX 向 bestX 转移的概率
及转移步长与两者之间的相对位置和信息素强度有关，其转移概率计算如下： 95
=e /eibest bestibestp
(4)
= -ibest best i (5)
式中 i ——蚂蚁在 i处的信息素强度，且 i best 。
- 4 -
中国科技论文在线
当蚂蚁向信息素强度高的地方移动时，可能会在移动过程中找到更优的解，此时可确定
第 i只蚂蚁向最好位置处转移时的步长为： 100
0( )
(-1,1)
i best i ibest
i
i sa
X X X p p
X
X rand l 否则
l
(6
式中
0p、 均为[0,1]区间的常数。所有的蚂蚁在完成全局搜索与局部搜索后，将进行
信息素强度更新，更新规则如下：
( , ) ( , )new i old i i (7)
式中
( , )n e w i ——更新后的信息素强度； 105
( , )old i ——原有的信息素强度；
——为信息素挥发系数，0 1 ；
i ——信息素增量。
其中 i 可按照下式计算：
sgn( ( )) ( )i i if X f X
(8) 110
式中 ——调节因子，一般可取0< <1 ；
( )if X ——更新前后目标函数差值， ( , ) ( , )( ) ( ) ( )i new i old if X f X f X 。
与经典搜索方法从一个孤立的初始点出发进行寻优的过程相比，该算法具有明显的优越
性和稳定性。
混合遗传蚂蚁算法 115
对比遗传算法与蚁群算法的特点可以发现，蚂蚁算法具有全局收敛，求解精度高的优点，
但是在初期由于信息素的缺乏，求解速度较慢[3]，而遗传算法虽然在算法后期容易发生冗余
迭代现象[4]，但是其算法搜索初期速度较快，两种算法具有较为明显的互补性。另外，对于
热网阻力特性辨识问题，初始解空间范围对于辨识效果具有至关重要的意义，由于初始解空
间范围过小有可能无法获取全局最优解，因此我们会将初始解空间范围定义的较大，而解空120
间范围过大又会导致计算耗时过长，甚至出现无法收敛，这就要求用于辨识的算法应该在计
算初期具备较高的收敛速度，能够快速的缩减解空间范围。基于以上原因，本文将遗传算法
与蚂蚁算法有机的融合，提出了一种新型算法――混合遗传蚁群算法。
125
130
135
- 5 -
中国科技论文在线
140
145
150
155
160
图 1 GASCAA 算法流程图
Fig. 1 Flow chart of program with GASCAA
算法的基本思路是在算法的前过程充分利用遗传算法的快速性、随机性、全局收敛性，
快速缩减解空间范围，生成有关问题的初始信息素分布，而在算法后过程充分利用蚁群算法165
并行性、正反馈性及求解效率高等的特点，求解问题的最优解。其具体的计算流程如图1所
示。
3 监测点优化布置研究
对于供热管网，各管段阻力特性变化是引起管段流量与节点压力变化的一个重要因素，
当管网某管段的阻力特性系数发生变化时，必然会引起全网各管段流量与节点压力发生变170
化，但其对于每一管段和节点的影响程度是不一样的，而由此引起变化最大的管段或节点极
可能就是布置监测点的最佳位置。基于这一思想，本文结提出了反映管网阻力特性变化的热
网监测点优化布置方法。
在实际阻力特性系数未知的情况下，选择管段设计阻力特性系数
d
S 作为测点优化布置
时的参考阻力特性系数
r
S 。管道设计阻力特性系数可根据下式直接求出： 175
3
10
dp zh
h
K
S l
d
(9)
式中
dp
S ——供、回水管网管段的设计阻力特性系数， 2Pa/(t/h) ；
K ——管壁的当量绝对粗糙度，m；
d ——管道的内径，m；
- 6 -
中国科技论文在线
zh
l ——管段的折算长度，m； 180
h
——水的密度， 3kg/m 。
流量监测点优化布置
设热网共有 B个管段，在 i管段处其参考阻力特性系数产生
i
S 变化值，则整个管网的
管段流量都会受到不同程度的影响，设其中被考察管段 j的流量变化值为
j
G 。则可用
/
j i
G S 表示 i管段处单位阻力特性系数变化在 j管段处产生的流量变化率，它反映了 j管185
段流量受其他管段阻力特性系数变化影响的大小。其物理意义表示由于 i管段参考阻力特性
系数的变化，引起管段 j流量波动的程度。
设 ( , ) /
G j i
X j i G S ，分别对
j
G 与
i
S 取偏微分有
( , ) / ( , 1, 2, , )
G j i
X j i G S j i B (10)
式中 ( , )
G
X j i 即为 i 管段阻力特性系数对 j 管段流量的影响度。 ( , )
G
X j i 可用矩阵190
,G g ji B Bx X 表示，称之为阻力特性系数对管段流量影响矩阵。
结合供热管网的特点，采用求解析解的方式计算
G
X 阵。即通过管段流量列向量G 对管
段阻力特性系数列向量
m
S 求偏导计算
G
X 。
G
/
m
X G S (11)
其中计算的关键为构造一转换矩阵 E 将管段阻力特性系数对角阵
1 2
( , )
B
diag s s s S195
转换为列向量 T
1 2
( , )
m B
s s s S
1 2 1ij ij ijk ijB BE e e e e (12)
其中元素 ijke ( 1,2, , ; 1,2 , ;i B j B 1,2, , )k B 是一个矩阵向量，其具体定义为：
1
0
ijk
i j k
e
当 时
其它
通过推导得到影响矩阵公式为 200
1 1 1 1 1
G S S S S
[ ( ) ]( )
T T
k k k k
X M A A M A A M M EG G (13)
即已知某一确定工况下的管段的阻力特性系数
m
S 与流量G ，则
G
X 可求。选择设计工
况作为基准工况，根据设计工况下的热网各管段的阻力特性系数
m
S 与流量计算求得
G
X 。
对于
G
X 阵中的第 j行第 i列元素 ,g jix ，其值越大表示 i管段阻力特性系数变化对 j管段流量
的影响也越大。 205
选择模糊聚类方法进行流量监测点优化。 根据模糊类聚理论，研究样本间的关系，需
1
m
1
0 0 0 0 00 0 0 0 0 0 0 0 0
0 0 0 00 0 0 0 0 0 0 0
= = + + + + 0 0 0 0 00 0 0 0 0 0 0 0 0
0 0 0 00 0 0 0 0 0 0 0
0 0 0 00 0 0 0 0 0 0 0 0 0
B
ijk k k
k
BB B B B B B
s
e s s
s
ES
- 7 -
中国科技论文在线
选择一个能反映研究对象之间亲疏关系的合适的统计量，即反映样本间相似程度的统计量，
根据这个量的大小形成分类系统。在此选择相似度作为分类的指标。管网两两管段之间的相
似程度关系用相似度可表示为：
, ,
2 1
, ,
1 2 2
, ,
1 1
1
(1 ( ) )/ 2 ( , 1,2, )
B
g ik g jkB
k
ij g ik g jk
B B
k
g ik g jk
k k
x x
x x i j B
B
x x
(14) 210
式中
,g ikx ， ,g jkx —— GX 阵中第 i j、 行各元素值。
采用编网法进行聚类。所谓编网法就是根据实际应用情况，选定某一阈值 [0,1] ，
作 γ 阵的截矩阵 γ 。 γ 阵的元素满足：
1
0
ij
ij
ij
(15)
确定截矩阵 γ 各元素后，将矩阵 γ 改为编网图，由对称性只取主对角线下面的一半元215
素，将对角线上的“1”改写成对应元素的编号，在对角线下方用“*”号代替“1”，“0”则
变成空格，其中画“*”号的位置称为结点，由各结点向矩阵的主对角线引竖直的经线和水
平的纬线，构成网络，经、纬线路沟通在一起的元素编成一片网，相应的元素聚为一类，同
一类元素间必有线路沟通,不同类元素间必无线路沟通，从而完成分类。
管网各管段通过模糊聚类完成分类后，在每一类组中选择一个最具代表性的点作为流量220
测点。设某一组别含m个管段，每个管段都与其余 1m 个管段存在欧氏距离。根据欧氏距
离的概念，与其余 1m 个管段平均欧氏距离最小的管段是最具代表性的，则选此管段作为
流量测点所在位置。
管段间的欧氏距离按下式计算：
2
, ,
1
1
( ) ( , 1,2, )
m
ij g ik g jk
k
r x x i j m
m
(16) 225
管段 i与其余 1M 个管段的平均欧氏距离为：
1
1
( 1,2, , )
1
m
i ij
j
i j
r r i m
m
(17)
最小平均欧氏距离按下式计算：
min min{ } ( 1,2, , )ir r i m (18)
压力监测点优化布置 230
设管网共有B个管段、N 个节点，在 i管段处其阻力特性系数产生 iS 变化值，则整个
管网除作为定压点的参考节点外的 1N 个节点的节点压力都将受到不同程度的影响，设其
中被考察节点 j的压力变化值为 jP 。则可用 /j iP S 表示 i管段处单位阻力特性系数变化
在 j节点处产生的压力变化率，其物理意义为 j节点压力随 i管段阻力特性系数的变化而发
生波动的程度。此时，可用 ( , )PX j i 表示 i管段阻力特性系数对 j节点压力的影响度，即： 235
- 8 -
中国科技论文在线
( , ) /P j iX j i P S (19)
分别对 jP 与 iS 取偏微分有
( , ) / ( 1,2, , 1; 1,2, , )P j iX j i P S j N i B (20)
各节点的阻力特性系数影响度 ( , )PX j i 可用矩阵 , ( 1)P p ji N Bx X 表示，称之为阻力
特性系数对节点压力的影响矩阵。其实质就是节点压力列向量 P 对管段阻力特性系数列向240
量
mS 的偏导数，即
P / m X P S (21)
通过矩阵变换得到
1 1 1
P S S( )
Tk k kX A M A A M EG G (22)
PX 阵中列元素值表示该列表示管段的阻力特性系数变化对所有节点压力的影响程度，245
行元素值则表示该行表示节点的压力受管段阻力特性系数变化的影响程度。阵中元素 ,p jix
越大，表示管段阻力特性系数的影响程度越大。
采用一个综合指标来反映管网中所有管段的阻力特性系数变化对某一节点的影响度，根
据数理统计的原理[5]，定义反映阻力特性系数变化的节点压力灵敏度为：
2
,
1
( 1,2, , 1; 1,2, , )
B
j p ji
i
x j N i B
(23) 250
式中 j —— j节点对管网所有管段阻力特性系数变化的灵敏度；
,p jix ——阻力特性系数影响矩阵 PX 中的元素。
某节点的灵敏度 越大，表示当管段阻力特性系数发生变化时，该节点的压力产生的
变化也越大。换言之，节点灵敏度高越高，就越能通过测量该节点的压力变化来了解管网阻
力特性系数的变化情况，该节点也就越可能成为反映阻力特性系数变化的压力监测点。从而255
为管网阻力特性系数的辨识提供更多的压力信息。
根据计算得到管网各节点的灵敏度，可通确定压力监测点的数目与位置，具体步骤为：
根据热网基础数据，计算管网各管段参考阻力特性系数；选择设计工况为基准工况，计算热
网各管段的阻力特性系数与流量，得到阻力特性系数对管段流量影响矩阵；计算反映阻力特
性系数变化的节点压力灵敏度；结合工程实际要求，根据灵敏度排序确定反映阻力特性系数260
变化的压力测点集。
4 算例分析
选择一多热源环状热网作为算例，其平面示意图如图 2 所示。
- 9 -
中国科技论文在线
图 2 算例热网平面图 265
Directed graph of plane sample heat-supply network
热网为一级管网，热源数 2s ，热力站数 13u 。供、回水平面网空间对称。热网的
空间结构如图 3 所示。
图 3 算例热网空间示意图
Sketch map of sample spatial heat-supply network 270
根据上文基于反映管网阻力特性系数变化对管段流量与节点压力影响程度的原则，最终
可计算确定算例热网的监测点布置方案。
流量测点所在管段名称：
{11-12, 4-5, 3 -16 ,14 -15 , 21-22, 16 -17 , 15-31,10 -27 u1, u7 ,u9, u11, u12, u13,s1, s2} ， ；
压力测点所在节点名称： 275
{ n26, n26 , n9, n8, n9 , n8 , n7, n10, n25 n7 , n27, n10 , n11, n25 , n27 ,n6 } ， 。
- 10 -
中国科技论文在线
图 4 热网监测设备分布示意图
Sketch map of monitor equipmenton the heat-supply network
采用文中提出的遗传算法与蚂蚁算法有机融合后的混合算法进行计算。算法在热网阻力280
特性辨识问题中应用的具体的流程与参数设置如下：
（1）遗传生成初始优化解：将所有管段按照管网热用户所在管段与管网管段、热源所
在管段划分为两类，辨识范围分别为0 30S 与0 50S ；由于遗传算法是为了快速
缩减解空间范围，因此选择较小的群体规模 popsize=60；锦标赛法选择操作，改进非一致交
叉，非重复一致变异，自适应交叉概率和变异概率参数 1 ＝ 、 2 ＝ 、 1 ＝ 、285
2 ＝ ；根据最大运行代数 maxruns=20 控制循环次数，生成一组优化解；
（2）蚁群算法求最优解：算法生成的优化解中两类管段的范围分别缩小为［0，10］与
［0，30］，因此设 5cR ，重新确定阻力特性系数范围为［0，15］与［0，35］；根据该范
围初始化信息素分布，取蚁群规模 M=40、信息素挥发系数 ＝；然后执行 CAA 相应
的操作；最后根据最大运行代数 maxruns=30 控制运算迭代次数，解码输出最优解。 290
混合优化算法的计算进程如图 5、图 6 所示，利用遗传算法快速收敛的特性，算法在最
初的 20 代快速收敛，利用蚁群算法良好的全局寻优特性，在第 31 代即找到最优解，最优解
对应的目标函数计算值为 。
混合优化算法作为一种混合型的算法，融合了遗传算法与蚂蚁算法的特点，算法的计算
环节较多，某些具体的设置与操作会对算法的计算效果产生较大的影响，因此对于不同的问295
题，算法的具体设置也可能不同。算法那在热网阻力特性辨识中应用时应注意以下几点：
300
- 11 -
中国科技论文在线
0 10 20 30 40 50
种
群
逐
代
最
优
个
体
目
标
函
数
值
GASCAA运算代数
0 10 20 30 40 50
种
群
逐
代
平
均
目
标
函
数
值
GASCAA运算代数
（1）算法群体规模与最大世代数的选取：由于遗传算法计算初期的收敛速度较快，为
了减少整体运算时间，算法初期采用遗传计算时其群体规模与最大世代数均不宜选得过大，
对于中小型热网，一般群体规模可取 80～100，最大世代数取 20 代即可，对于大型的热网305
可相应的增加。
（2）初始解空间的确定：对于热网阻力特性辨识问题，其初始解空间的确定对于计算
的效率与精度具有十分重要的影响。对于不同的热网，其阻力特性系数初始解空间的确定主
要参考设计参数以及具体的实际工程经验确定。一般可选择设计参数作为参考依据，并在此
基础上采用两种方法确定初始解空间。第一种方法是事先搜寻整个管网的最大设计阻力特性310
系数 maxS 与最小设计阻力特性系数 minS ，然后确定初始解空间。第二种方法是本文采用的
方法，对管网各管段按照设计阻力特性系数值进行划分归类，确定各类管段的阻力特性系数
变化区间。第一种方法管网所有管段阻力特性系数的求解区间统一，操作比较简单方便，并
且由于整个解空间区间比较大，减少了计算无法搜寻全局最优解现象出现的几率。但是这种
方法无形中增大了搜索空间，导致运算时间加长、计算速度降低，当热网规模很大时，有可315
能出现在指定的的计算条件下无法收敛的情况。与第一种方法相反，第二种方法需要对管网
各管段进行分类，根据不同的分类划分其求解区间，操作比较复杂。同时由于分类使得解空
间区间变小，也增加了无法搜寻全局最优解的可能性。但是，由于解空间区间变小，大大减
少了运算时间，因此对于大规模的热网阻力特性辨识还是很有实际意义的。
5 结论 320
基于最小二乘原理建立了供热系统阻力特性辨识数学模型，将热网阻力特性辨识问题转
化非线性优化问题。提出了基于反映管网阻力特性变化的热网流量监测点与压力监测点优化
布置方法，以模糊聚类分析理论为基础，在反映管段阻力特性系数变化对管网运行参数影响
的前提下，测点选择包含了热源与热用户信息，增强了该优化方法的针对性与实用性。对遗
传算法与蚂蚁算法进行改进构造的新型混合算法，收敛速度与计算精度获得了较大的提升。 325
图 5 混合优化算法平均辨识
目标函数值变化
The change of every generation
average objective function value of GASCAA
图 6 混合优化算法最优个体
辨识目标函数值变化
The change of every generation optimum objective
function value of GASCAA
- 12 -
中国科技论文在线
[参考文献] (References)
[1] 关旭，张春梅，王尚锦. 一种改进的自适应遗传算法[J]. 微机发展，2003，11：4l-44
[2] M. Srinivas, L. M. Patnaik. Adaptive Probabilities of Crossover and Mutations in GAs[J].. IEEE Trans. On
SMC, 1994, 24(4): 656-667
[3] 吴庆洪，张纪会，徐心和．具有变异特征的蚁群算法[J].．训算机研究与发展．1999, 36（10）：1240-1245 330
[4] 李敏强，徐博艺，寇纪淞．遗传算法与神经网络的结合[J].．系统工程理论与实践．1999, 19( 2) ：65- 69
[5] 何灿芝，喻胜华．应用统计[M]．湖南科技出版社，1997
