量调节方式下区域热电系统的联合最优潮流

韩 赫，张沛超，杜 炜，刘学智，徐博强，姚 程

（1.电力传输与功率变换控制教育部重点实验室（上海交通大学）,上海市200240；2.国电南瑞科技股份有限公司,江苏省南京市211106；3.国网天津市电力公司电力科学研究院,天津市300384）

0 引言

区域供热网（district heating network,DHN）在现代社会中具有重要作用。DHN 可以与其他能源形式构成综合能源系统［1-2］,例如通过热电联供（combined heat and power,CHP）机组构成区域热电系统（district heat-electric system,DHES）,从而实现多能耦合,有效提高能源利用效率［3］。

DHN 有质调节和量调节这2 种基本调节方式。前者采用“变温恒流”,而后者采用“恒温变流”［4］。由于质调节的水力工况稳定、易于控制,现有研究多集中于质调节方式。例如,文献［5］假设管道热水流量恒定,通过改变CHP 机组供热蒸汽流量来调节管网供热温度；文献［6］在质调节模式下提出考虑机组调峰主动性的电热协调调度方法；文献［7-8］基于广义电路分析理论,建立多能源系统的支路及网络模型,描述热网的动态特性。然而,质调节存在热传输时延且单纯使用时易造成能源浪费［9］。对此,文献［10］提出应在供暖期不同阶段采用不同的调节方式。例如,在供暖初、末期负荷较小时,系统以最小设计流量按质调节运行；在供暖中期,当供热温度达到设计值后按量调节运行。可见,实际运行中量调节与质调节方式应构成有效互补。但是,针对质调节的研究以恒流作为基本假设,难以应用于量调节。为此,本文专门研究量调节方式。

量调节方式下由于热力工况稳定,无须考虑热传输时延问题。但因流量可变,其水力模型难以线性化处理。同时,其热力模型也为非线性,并与水力模型存在耦合。故采用量调节的文献均使用非线性、非凸模型［11-12］。但一般的非凸优化问题缺乏稳定、高效的解法,且难以得到全局最优解。为此,在保证足够精度的前提下,本文研究如何将上述一般的非凸优化问题转换为易于求解的问题。本文的主要贡献如下。

1）提出水力支路压降的二阶锥松弛方法。水路管道压降与流量呈平方关系,文献［13］假设水力工况为层流,从而实现压降方程线性化。但DHN 为避免流量过小引发水力失调,实际工况常处于紊流阻力平方区［14］,线性近似误差较大。本文提出基于凸包络的压降方程松弛方法,并证明该方法为精确松弛。

2）提出热力支路的等效热损方法。针对DHN水力、热力耦合的非线性因素,文献［13,15］先采用定流量求解热力支路模型,再迭代更新流量直至算法收敛。本文根据量调节特点,提出热网等效热损模型,具有较高精度且无须交替求解水、热模型。

3）在上述基础上,建立DHES 联合最优潮流的凸凹规划（convex-concave programming,CCP）模型,并将其转换为二阶锥规划（second-order cone programming,SOCP）问题进行序贯求解,显著提高了求解效率。

1 DHN 模型及其凸化方法

DHN 由空间同构的供、回水网组成［11］,如图1所示。本文DHN 采用量调节方式,热源温度Tsrc和热负荷出口温度To恒定。Ts1,Ts2,Ts3分别表示节点1,2,3 的供热温度,热水流经热负荷后,温度从供热温度降为To。Tr1,Tr2,Tr3分别表示回水网中节点1,2,3 的回热温度,由于回水网温度低、热损小,可假设各节点回热温度近似一致且等于To［11,13］,这样研究中仅需关注供水网。本章研究如何对DHN 中的非线性因素做松弛和近似变换。

图1 DHN 结构Fig.1 Structure of DHN

1.1 支路模型

1.1.1 水力支路模型松弛方法

管道压降可用Darcy-Weisbach 公式描述［13］,即

式中：h和mb分别为管道压降和流量（即质量流率）；L和D分别为管道长度和内径；ρ为热水密度；g为重力加速度；λ为Darcy 摩擦系数；Λ为阻力系数。

由式（1）可知,h与mb符号一致。在紊流阻力平方区下,λ由Prandtl-Karman 公式确定,为只取决于管道自身特性、不随流量变化的常数［14］,即

式中：ε为管道绝对粗糙度。

式（1）为非仿射等式约束。为保证水力工况的相对稳定,本文假设管道流向在一个调度周期内保持不变,且可以在日前获知［16］。设该流向为正方向,然后将式（1）松弛为二阶锥不等式约束,即

上述松弛扩大了问题的可行域,对此将在1.2 节解决。

1.1.2 热力支路模型近似方法

管道温度下降方程描述了热水与外界环境的热交换过程,表现为热水流经管道后温度下降,即

式中：Tstart,Tend,Ta分别为管道始端温度、末端温度及环境温度；β为管道热传导系数；σ为热水比热容。

进一步可得管道热损Φloss为：

流量和温度的乘积项表现了DHN 中水、热网的非线性耦合关系。本文提出热网等效热损模型,方法如下。

首先,通常情况下有βL≪σmb,因此可将式（4）中指数项做一阶泰勒展开［16-18］,即

将其代入式（5）,可得：

其次,由于量调节方式下热网热力工况稳定,且随着DHN 向第4 代发展［19］,热损较低,节点供热温度Ts一般仅比热源温度Tsrc低1～2 ℃［4,15］,远小于供给热负荷的温降。因此,本文将各节点供热温度均近似为Tsrc,从而式（7）可改写为式（8）,式（8）实际上代表管道对环境的稳定热传导。

最后,将上述热损等效叠加至管道末端的热负荷上,如图2 所示。经过以上处理,在量调节方式下可实现DHN 中水、热的解耦分析。需要指出的是,上述热力支路近似模型也会给管道流量带来一定的计算误差,详细分析见附录A。

图2 热力支路等效热损模型Fig.2 Equivalent heat loss model of thermal branch

1.2 拓扑约束模型

1.2.1 水力拓扑约束

类似基尔霍夫第一定律,各节点应满足如下流量平衡方程。

式 中：mni为 节 点i的 净 注 出 流 量；mbk为 管 道k的 流量；aik为节点关联矩阵A中的元素,定义为

类似基尔霍夫第二定律,环路压降平衡方程要求任意环路的管道压降代数和应等于0,即

式中：hk为管道k的压降；bk为环路关联矩阵B中的元素,定义为

为收紧式（3）造成的松弛间隙,本文在目标函数中添加惩罚项fh,通过最小化fh可实现零间隙,即

式中：Ω为惩罚系数；Λk为支路k的阻力系数。

收紧管道压降方程松弛如图3 所示,由图3 可见,式（3）的实质是用凸包络（图中蓝色区域）代替原可行域以达到凸化目的；而惩罚项又重新收紧了可行域,实现了精确松弛。详细证明见附录B。

图3 收紧管道压降方程松弛Fig.3 Tightening relaxation of pipe pressuredrop equation

1.2.2 热力拓扑约束

在管道的连结节点处,根据焓守恒原理可以导出如下方程。

式中：min和mout分别为流入、流出混合节点的管道流量；Tin为流入节点管道的末端温度；Tout为流出节点管道的始端温度。

但在量调节方式下,由于各节点供热温度近似一致,式（14）可以用式（9）的流量平衡方程代替。

2 DHES 联合潮流模型

2.1 配电网模型

交流潮流模型为典型的非线性方程组。节点i的功率注入模型为：

式中：Pni和Qni分别为节点i的净注入有功和无功功率；Vi,Vj,θij分别为节点i和j的节点电压和相角差；Gij和Bij分别为节点导纳阵第i行、第j列元素的实部和虚部。

对于辐射型配电网,可将节点注入模型式（15）—式（16）改写成DistFlow 支路潮流模型［20］。同时,可仿照热网,对线损做相同等效处理,可得到如下潮流方程。

式中：Pbk和Qbk分别为支路k的有功和无功功率,其中支路k由节点i流向节点j；J为节点j的流出支路k'的集合；Rk和Xk分别为支路k的电阻和电抗；lk为支路k的电流平方值；ui为节点i的电压平方值。

将式（20）松弛为具有二阶锥形式的不等式约束式（21）,从而建立电网潮流方程的凸优化模型［21］。

文献［21］证明了对于辐射型配电网,上述支路潮流松弛为精确松弛。

2.2 DHES 联合潮流模型的矩阵形式

至此,可将DHES 联合潮流模型统一表达为矩阵形式。定义DHES 的联合节点关联矩阵为：

式中：AE和AH分别为配电网和DHN 的节点关联矩阵。可见,矩阵A仍保留了强稀疏性。

定义Apos为矩阵A的非负部分,表征节点与流入节点的支路之间的对应关系,其元素apos,ik为：

DHES 联合潮流模型统一将支路的功率损失等效到支路末端的负荷上,这样,DHES 的系统功率平衡约束可以统一写成：

式中：Pb,Pn,Ploss,Psrc,Pload分别为有功功率列向量P的支路潮流、节点潮流、支路损耗、节点电源、节点负荷列向量；无功功率列向量Q及热功率列向量Φ以此类推。

各类支路损耗分别为：

式中：R和X分别为支路电阻和电抗的列向量；l为支路电流的平方列向量；L和β分别为支路管道长度和管道热传导系数列向量。向量中无该类参数及变量的元素取为0；“∘”运算符表示2 个向量按元素相乘。

DHN 环路压降平衡约束为：

式中：h,Λ,mb分别为支路管道压降、支路阻力系数、支路管道流量的向量形式。

其余电、热潮流约束也写为矩阵形式,即

式中：u为节点电压平方值列向量。

最后,各变量还需满足如下运行约束。

式中：Pb,max,Qb,max,Imax分别为支路有功潮流、无功潮流、电 流 上 限 列 向 量；Vmax,Vmin和mb,max,mb,min分 别为节点电压和支路流量的上、下限列向量。

3 DHES 联合最优潮流

3.1 CCP 模型

DHES 优化的目标函数为：

式中：fP,t,fΦ,t,fchp,t分别为t时段的购电成本、购热成本和CHP 成本；fh,t为式（13）中t时 段的压降惩罚成本。

购电与购热成本皆为线性函数,即

式中：CP,t和Pbuy,t分别为t时段的电价和购电量；CΦ,t和Φbuy,t分别为t时段的热价和购热量。

CHP 运行成本包括发电、产热和耦合成本项［16］,可表示为：

式中：α0至α5为成本系数；Pchp,t和Φchp,t分别为t时段CHP 机组电、热出力。

考虑如下约束条件。

1）热电联合潮流约束

热电联合潮流约束见式（24）—式（35）。

2）设备运行约束

CHP 设备运行约束见文献［22］。储能（包括电储能和热储能）运行约束见文献［23］,针对储能的非线性充放互斥约束,该文献提出了一种线性松弛方法。

3）购能约束

式中：Pbuy,max和Φbuy,max分别为购电和购热上限。

分析上述优化问题（记为问题P1）可知,其约束条件皆为仿射约束或二阶锥约束；在目标函数式（36）中,惩罚项fh,t为凹函数,其余为凸函数。这样,基于第1 章和第2 章的工作,本文将DHES 联合最优潮流转换为一类特殊优化问题——CCP 问题［24］。虽其仍为非凸问题,但却有可靠的求解方法,并已在电力系统中获得应用［20］。

3.2 序贯SOCP 求解方法

本文将上述CCP 问题转换为可序贯求解的凸优化子问题,具体步骤如下。

步骤1：首先忽略式（36）中的惩罚项fh,t,得到凸优化子问题（记为问题P2）,求解问题P2并获得各支路流量的初值。问题P2保留了问题P1的全部约束条件,则所求得初值必为问题P1的可行解,这能确保问题P1在后续迭代过程中收敛［24］。

步骤2：考虑惩罚项fh,t,但将其在当前流量值Mb处做一阶泰勒展开,得到凸优化子问题并求解。该凸优化子问题的目标函数（记为问题P3）如式（42）所示。

式中：hk,t和mbk,t分别为t时段管道k的压降和流量；Mbk,t为Mb中 元 素。

在这一步中,仅对惩罚项做了线性处理,尽量保留了问题P1的信息,这有助于快速收敛。

步骤3：根据优化结果更新流量值,重复步骤2,直至目标函数值改变量小于容差（本文设为10-3）。至此得到问题P1的最优解,该最优解已重新收紧水力支路压降的可行域。

但在凸优化子问题P2和P3中,目标函数中仍包含CHP 的二次成本函数项。本文引入辅助变量δt,并增加二阶锥不等式约束,将最小化fchp,t转变为如下等价问题。

这样,子问题P2和P3被进一步变换为SOCP 问题,变换过程详见附录C。最终,问题P1被转换为序贯SOCP 问题。相比一般凸优化问题,SOCP 问题具有非常高效的求解方法。

4 算例验证

4.1 算例设置

DHES 算例的网络拓扑如附录D 图D1 所示,包含IEEE 33 节 点 配 电 网 和 巴 厘 岛32 节 点DHN［11］。热电耦合设备为背压式CHP 机组和抽凝式CHP 机组,其热电可行域参数见文献［22］,电热出力范围均为0.2～1 MW,成本系数见附录D 表D1。

电网参数为：节点1 购电功率上限为3 MW；节点25 连接风电场；节点6 配置电储能,参数见附录D表D2；节点3,8,13,28 均配置无功补偿装置,最大容量均为0.6 Mvar；线路潮流有功和无功功率的上限分别为3 MW 和2 Mvar；设置节点1 为平衡节点,基准电压为12.66 kV,其余节点电压标幺值范围为0.95～1.05；负荷及主网电价见图D2,风电预测曲线见图D3。

热网参数为：节点1 向一级热网购热功率的上限为1 MW；节点5 连接太阳能热厂；节点14 配置热储能,参数见附录D 表D2；管道流量下限取为2.2 倍内径,以保证雷诺系数不低于104（紊流区）,上限取为6 kg/s,其余参数见文献［11］；热源温度为设计温度70 ℃,热负荷出口温度为30 ℃；热负荷及主网热价见图D2,环境温度及光热预测曲线见图D3。

总调度周期为24 h,单个滚动优化周期包含6 个时段,每时段为1 h。

4.2 模型精度分析

本节验证本文第3 章所建立的CCP 模型的求解精度。对比模型为未做松弛与近似变换的非凸精确模型,上述2 类模型的区别见附录C 表C1。本节设置如下3 个算例：算例a 采用CCP 模型；算例b 采用非凸精确模型,取零初值；算例c 同样采用非凸精确模型,但以算例a 的结果作为初值。由于初值选取恰当,认为算例c 求得了全局最优解。在求解方面,算例a 采用本文3.2 节的序贯SOCP 方法,SOCP子 问 题 的 求 解 器 为SeDuMi［25］；而 算 例b 和c 均 使 用内点法求解。

1）水路压降松弛方法验证

CCP 模型中环路管道计算压降的代数和（取正值）如附录E 图E1 所示。计算压降特指由流量mb优化结果所计算的压降值Λm2b。由图E1 可见,各时段环路压降代数和均在10-8～10-6量级之间,近似为0,证明本文对水路压降的松弛方法能有效满足环路压降平衡方程。表E1 中还给出了CCP 模型与非凸精确模型关于流量与压降的详细对比结果。

2）热网等效热损模型验证

算例a 与c 的管网总热损对比如图4 所示。由图4 可知,由于CCP 模型将节点供热温度近似为热源温度,导致热损始终偏高。进一步对照附录D 图D2 中热负荷曲线及图D3 中环境温度趋势可以发现,外界环境温度越高,DHN 热损越小；热负荷越小,CCP 模型的热损相对于精确热损的误差越大,但最大误差不超过2%。考虑到总热损仅占总负荷的3%～5%,上述绝对误差很小。而附录E 表E1 已表明水力计算结果误差很小。综上可见,本文虽然对水力、热力做解耦建模并独立分析,但等效保留了管道热水向外散热的损耗影响,计算结果足以满足工程需要。

图4 各时段总热损对比Fig.4 Comparison of total heat loss in each hour

为了验证热网在量调节方式下热力工况稳定的特点,本文基于非凸精确模型,给出了调度期间各节点温度的实际变化与分布情况,见附录E 图E2。

3）优化结果验证

各时段总成本对比如图5 所示。比较算例a 和c可见,两者求得的系统总成本几乎一致,绝对误差均保持在10-2量级,最大和最小相对误差分别出现在时段5（0.046%）和时段20（0.004%）。这说明本文方法保留了很高的精度。设备出力与成本详细对比结果见附录E 表E2。

图5 各时段总成本对比Fig.5 Comparison of total cost in each hour

进一步对比算例b 和c,发现部分时段算例b 的总成本明显偏高。这说明虽然算例b 和c 均采用精确模型,但由于所建立的优化问题非凸,其结果受初值影响较大,容易陷入局部最优。

4.3 模型计算效率分析

本节基于算例a 和c 进一步验证模型的计算效率。统计各算例在每个调度周期的平均求解时间发现,算例c（非凸精确模型+内点法）的平均求解时间 为49.96 s,而 算 例a（CCP 模 型+序 贯SOCP 方法）的平均求解时间为6.21 s。可见,本文方法在保证精度的同时,还显著提高了计算效率,这非常有利于大规模DHES 的优化计算。

附录E 图E3 为部分时段序贯SOCP 的收敛过程。由收敛曲线可见,求解初值时（问题P2）因不考虑压降惩罚项,目标函数值误差很大,说明经松弛后压降求解结果远远偏离与流量的二次关系曲线；当进入序贯SOCP 迭代后（问题P3）,惩罚项迫使压降松弛间隙逐渐逼近0,目标函数值经不超过10 次的迭代后快速收敛。

5 结语

量调节是DHN 的一种基本调节方式,本文针对其特点提出了DHES 的联合最优潮流模型及其高效解法,主要结论如下。

1）量调节方式下难以对水力模型做线性化处理。可利用凸包络对热网的水力压降方程进行松弛,并通过设计惩罚函数来收紧松弛间隙。

2）量调节方式具有热力工况稳定的特点,据此可实现稳态水力、热力模型的解耦,结果表明该近似具有很高精度,并能等效计及管网热损。

3）DHES 的联合最优潮流模型可被表征为特殊的非凸优化问题,CCP 问题进而转换为序贯SCOP问题。结果表明该求解方法能显著提高计算效率,有利于大规模DHES 的运行优化。

本文和已有文献均假定在一个调度周期内仅有一种调节方式。下一步将研究量调节与质调节的最优协调与方式切换问题。

附录见本刊网络版（http：//www.aeps-info.com/aeps/ch/index.aspx），扫英文摘要后二维码可以阅读网络全文。