本发明属于能源互联网领域，涉及一种用于多能互补综合能源系统的动态潮流分析方法及动态潮流分析系统。
背景技术：
近年来，以能源互联网为特征的第三次工业革命已蓄势待发。能源互联网的特征是以电网为主干和平台，进行多种能源形式的耦合互补。多能耦合系统的网络分析是能源互联网领域的重要研究内容之一，是进行系统规划、运行调控、能源交易的计算依据与基础，这里不仅要考虑系统中供能侧、用能侧各自的能量转换互补，更要考虑多能系统的网络平衡。目前，在传统的电、热、气等各自领域，各系统分析方法相对成熟，例如，电力系统采用潮流计算；热力系统遵循流体和热力学定律，通过联立水力(或蒸汽)和热力两组方程计算；天然气系统遵循流体力学定律，用流体力学方程表征计算。此外，随着风机、光伏和chp/cchp等可再生能源/清洁能源的大量建设以及主动性配网、智能微网的兴起，也提出了较多针对单个设备和小型电力网络的建模及仿真方法，部分微网系统由于存在多种能源形式的综合利用，也提出了相应的数学模型。上述大部分模型主要从设备角度出发，考虑了能量、动量以及质量等“点”平衡，但均没有考虑网络平衡，并不适于系统级的分析计算。
相比而言，文献[1]和[2]分别研究了电网与天然气网络、热网的联合潮流分析，文献[3]、[4]使用网络流模型对美国的煤、天然气和电网进行了仿真分析，这些文献从系统角度出发，考虑了网络平衡以及系统的多能特性，其中的电网部分，考虑的主要是传统交流电网中的一些典型设备，计算方法也是以传统的牛顿-拉夫逊法为主，对于电网设备种类多样(如新能源电场、facts器件等)或采用交直流混联电网结构等复杂情况，上述方法未见说明，此外，现有技术也未提及用户可以自己开发模型并接入计算的方法。
另一方面，目前在电、热各自的领域都有专业的仿真软件进行辅助分析，例如电领域中的psasp、bpa、pscad等，热领域中的thermoflow、ansys、cycle-tempo等，然而这些软件并未考虑除自身领域以外的多能之间的耦合。
[1]martinez-maresa,fuerte-esquivelcr.aunifiedgasandpowerflowanalysisinnaturalgasandelectricitycouplednetworks[j].ieeetransactionsonpowersystems,2012,27(4):2156-2166.
[2]xliu,njenkins,jwu,etal.combinedanalysisofelectricityandheatnetworks.energyprocedia,2014,61:155-159.liux.combinedanalysisofelectricityandheatnetworks[d].cardiffuniversityinstituteofenergy,2013.
[3]quelhasa,gile,mccalleyjd,etal.amultiperiodgeneralizednetworkflowmodeloftheu.s.integratedenergysystem:parti—modeldescription[j].ieeetransactionsonpowersystems,2007,22(2):829-836.
[4]quelhasa,mccalleyjd.amultiperiodgeneralizednetworkflowmodeloftheu.s.integratedenergysystem:partii—simulationresults[j].ieeetransactionsonpowersystems,2007,22(2):837-844.
技术实现要素：
有鉴于此，本发明提供一种用于多能互补综合能源系统的动态潮流分析方法及动态潮流分析系统，通过构建多能互补综合能源系统的动态潮流模型及采用混合求解算法求解对多能耦合系统进行网络分析，能够有效支撑系统级的分析计算。
为达到上述目的，本发明的技术方案如下：
一种用于多能互补综合能源系统的动态潮流分析方法，包括如下步骤：
步骤一，基于图论将非电系统抽象成图进行网络拓扑分析，构建非电系统准稳态模型；
步骤二，构建电力系统潮流模型；
步骤三，用功率输出外特性稳态模型或一阶惯性动态模型表征各耦合设备，构建非电系统与电力系统之间进行耦合的耦合设备功率输出外特性稳态模型或一阶惯性动态模型；
步骤四，联立步骤一非电系统准稳态模型、步骤二电力系统潮流模型以及步骤三耦合设备功率输出外特性稳态模型或一阶惯性动态模型，构建多能互补综合能源系统的动态潮流模型；
步骤五，采用混合求解算法求解步骤四多能互补综合能源系统动态潮流模型，混合求解算法指非电系统潮流模型和电力系统潮流模型分别用不同迭代算法进行求解。
进一步的，步骤一中将非电系统中管道对应图的有向边，连接件对应图的顶点，每个管段定义流量正方向，用图矩阵描述非电系统拓扑结构，并将阀门作为管道的附属属性进行计算。
进一步的，非电系统准稳态模型由非电系统静态潮流模型和非电系统动态潮流模型构成，其中非电系统节点采用稳态方程描述，非电系统管道采用针对温度、热量状态量的动态方程
进一步的，非电系统静态潮流模型为：
其中，a为关联矩阵，au、ad分别是上关联矩阵与下关联矩阵，bf为回路矩阵。m*为b阶管道流量列向量，q*为n阶节点流入流量列向量，δh*为b阶管道压差列向量，z*为b阶管道首末端节点高度差列向量，te为b阶管道末端温度列向量(℃)，tn为n阶节点温度列向量(℃)，qj*为n阶节点热负荷列向量，ta为b阶环境温度列向量(℃)，e为b阶温度衰减系数对角阵，s为b阶管道阻力系数列向量，hp*为泵扬程列向量，hv*为阀门两侧压差列向量。
进一步的，非电系统动态潮流模型为：
其中，a为关联矩阵，
进一步的，非电系统与电力系统之间进行耦合的耦合设备功率输出外特性稳态模型为：
其中，pg*、qg*分别为联供机组电、热功率，pp*为泵电功率，η为泵效率，m*为泵内流体流量，hp*为泵扬程，qhp*、qeb*、qc*分别为热泵、电锅炉、制冷机的热功率，php*、peb*、pc*分别为热泵、电锅炉、制冷机的电功率。
进一步的，非电系统与电力系统之间进行耦合的耦合设备一阶惯性动态模型为：
其中，τ是时间常数，s为拉普拉斯算子，pg*、qg*分别为联供机组电、热功率，pp*为泵电功率，η为泵效率，m*为泵内流体流量，hp*为泵扬程，qhp*、qeb*、qc*分别为热泵、电锅炉、制冷机的热功率，php*、peb*、pc*分别为热泵、电锅炉、制冷机的电功率。
进一步的，多能互补综合能源系统的动态潮流模型为：
其中，fe＝0表示步骤二电力系统潮流模型，fh＝0表示步骤一非电系统非电系统准稳态模型，feh＝0表示步骤三耦合设备功率输出外特性稳态模型或一阶惯性动态模型，f＝0表示多能互补综合能源系统动态潮流模型。
进一步的，混合求解算法包含如下步骤：
(1)静态潮流计算过程：
进行非电系统静态潮流计算，计算收敛后，根据非电系统收敛计算结果通过耦合设备功率输出外特性稳态或一阶惯性动态模型修正耦合器件电功率；
进行电力系统潮流计算，收敛后，再根据电力系统收敛计算后耦合器件电功率通过步骤三耦合设备功率输出外特性稳态或一阶惯性动态模型修正电热耦合器件非电状态量；
若前后两次计算得到的耦合器件电功率误差不满足精度要求，则刷新电热耦合器件非电状态量后继续进行计算直至超出最大迭代次数，若满足，则停止迭代，计算结果将作为动态计算的0时刻初值；
(2)动态潮流计算过程：
检测该时刻是否存在扰动，若存在则将扰动量经过换算计入非电系统动态潮流模型中；
进行非电系统动态潮流计算，计算收敛后，根据非电系统收敛计算结果通过耦合设备功率输出外特性稳态或一阶惯性动态模型修正耦合器件电功率；
进行电力系统潮流计算，收敛后，再根据电力系统收敛计算后耦合器件电功率通过耦合设备功率输出外特性稳态或一阶惯性动态模型修正电热耦合器件非电状态量；
判断前后两次计算得到的耦合器件电功率误差是否满足精度要求或超出最大迭代次数，若否，则刷新电热耦合器件电功率，继续进行计算；若是，则进一步判断是否已经到达截止时间，若是，则计算结束，输出结果，否则进入下一时步计算。
具体的，本发明具体实施方式采用的混合求解算法包含如下步骤：
(1)对多能互补综合能源系统拓扑关系进行分析，获取关联矩阵；
(2)对模型计算数据进行数据标幺化处理，从而完成计算数据准备；
(3)初始化迭代次数k＝0，初始化耦合设备电力系统状态量xe(k)及热网系统状态量xh(k)，且初始值随机分配。
(4)进行热力管网静态潮流计算，然后判断静态潮流计算是否收敛，如果收敛跳入步骤(5)；否则计算结束；
(5)根据热力管网静态潮流计算结果修正耦合器件电功率pe(k)，记为pe(1)(k)，并利用耦合器件电功率pe(1)(k)进行电力系统潮流计算，然后判断是否达到收敛，若是，跳入步骤(6)；否则计算结束；
(6)记电力系统潮流计算收敛后耦合器件电功率pe(2)(k)，根据耦合器件电功率pe(2)(k)修正耦合器件热网状态量xh(k)；然后判断是否|pe(1)(k)-pe(2)(k)|＜ε或k＞kmax，其中ε，kmax取经验值；若k＞kmax，则计算结束；若|pe(1)(k)-pe(2)(k)|＜ε且k≤kmax，则跳入步骤(7)进行动态计算，且当前时刻t＝t+dt，其中dt表示一个时步的度量；否则，k＝k+1后返回步骤(4)；
(7)判断是否存在扰动，若存在，处理扰动信息，更新热网状态量xh(k)后，跳入步骤(8)；否则直接跳入步骤(8)；
(8)进行热力管网动态潮流计算，判断计算是否收敛，若是，跳入步骤(9)，否则计算结束；
(9)根据热力管网动态潮流计算结果修正耦合器件电功率pe(k)，记为pe(3)(k)，并利用耦合器件电功率pe(3)(k)进行电力系统潮流计算，然后判断电力系统潮流计算是否能够计算收敛，若是跳入步骤(10)，否则计算结束；
(10)记电力系统潮流计算收敛后耦合器件电功率pe(4)(k)，根据耦合器件电功率pe(4)(k)修正耦合器件热网状态量xh(k)；然后判断是否|pe(3)(k)-pe(4)(k)|＜ε或k＞kmax，其中ε，kmax取经验值；若是跳入步骤(11)；否则k＝k+1后返回步骤(8)；
(11)判断是否t>tend，其中tend表示截止时间，若是，计算结束，输出计算结果；否则，t＝t+dt跳入步骤(7)。
本发明还提供一种多能互补综合能源系统动态潮流分析系统，包括电力分析模块、非电力分析模块和耦合分析模块，非电力分析模块用于实现步骤一非电系统准稳态模型，电力分析模块用于实现步骤二电力模块潮流模型的潮流计算，耦合分析模块用于实现非电模块与电力模块之间耦合设备功率输出外特性稳态模型或一阶惯性动态模型的计算；电力分析模块、非电力分析模块和耦合分析模块数据交互实现步骤四和步骤五多能互补综合能源系统的动态潮流模型的求解。
进一步的，非电力分析模块包括管道及节点两个基类，在此基础上设计包含依赖于上述两基类的设备模型类，用于指代同一类设备，然后设计设备模型类的子类模型实例类，指代实际设备，设计模型实例类的子类模型实例集合类，实现各实际设备模型对象的集合，从而形成非电力分析模块。
进一步的，将非电力分析模块以及耦合分析模块的模型计算单独编译成为动态链接库文件，通过现有电力系统仿真分析软件通过其用户程序接口调用该文件以实现联合计算；现有电力系统仿真分析软件作为电力分析模块。
进一步的，将现有电力系统仿真分析软件作为电力分析模块单独编译成可执行程序，通过嵌入到非电力分析系统以及耦合分析系统的模型计算程序中，执行该程序以完成联合计算。
进一步的，用户通过图形化界面建立设备的自定义模型，系统自动生成模型文件，动态链接库程序解析模型文件并内置算法，并通过与系统主程序的接口实现数据的相互传递与协同计算。
进一步的，系统主程序接入用户自定义模型的接入程序计算流程包括：
(1)记时间t＝t，其中t为当前时刻，主程序完成自身动态潮流计算流程，并在判断已经达到收敛后，通过用户接口调用用户自定义模型，跳转至步骤(2)，若不能达到收敛，则程序结束；
(2)完成自定义模型的解析、赋初值及潮流计算流程，若潮流计算达到收敛，则跳转至步骤(3)，否则程序结束；
(3)判断收敛后是否满足程序结束约束条件，即将收敛后自定义模型的输出变量值和对应主程序接口变量值进行比较，若误差小于设定阈值则t＝t+dt，主程序继续进行下一时步动态潮流计算，否则程序结束。
本发明的有益效果是：
(1)本发明通过构建多能互补综合能源系统的动态潮流模型及采用混合求解算法求解对多能耦合系统进行网络分析，能够有效支撑系统级的分析计算；并考虑了综合能源系统时域上的状态变化因素，进行基于时域的动态潮流滚动计算，提高了计算的准确性。
(2)本发明针对不同系统采用不同迭代方法的混合求解算法，例如电力潮流采用最优乘子法，热力潮流采用牛顿-拉夫逊法，保证收敛并且不牺牲效率。
(3)本发明采用嵌入式程序开发方法，具有开放性和可扩展性，将其他能源形式网络以及耦合设备的模型计算单独编译成为动态链接库文件，电力系统仿真分析软件通过其用户程序接口调用该文件以实现联合计算，或将电力系统仿真分析软件单独编译成可执行程序，通过嵌入到其他能源形式网络以及耦合设备的模型计算程序中，执行完成联合计算。
(4)本发明采用调用动态链接库文件的方式外接用户自定义建模，具有开放性和可扩展性，用户通过图形化界面建立设备的自定义传递函数数学模型，程序自动生成模型文件，动态链接库程序解析模型文件并内置算法，并通过与主程序的接口实现数据的相互传递与协同计算；用户自定义建模采用面向对象的程序架构，本身也具有开放性和可扩展性。
(5)本发明中一方面利用现有成熟电力系统仿真软件计算功能，可以进行包括直流系统、新能源电场、facts器件等在内的多种电网设备仿真，减少了程序开发的工作量，同时可以采用更为先进的计算方法，保证了计算的收敛性和可靠性；另一方面，拓展了电力系统仿真软件功能，使其能够进行多能耦合系统的仿真计算；再一方面，可以依托电力系统仿真软件强大的模型库，在联合潮流计算基础上对多能耦合系统进行进一步的稳态或动态分析，而且具备开放性和可扩展性，本领域技术人员可以在此基础上进一步开发自定义模型参与仿真计算，极大地扩展了针对综合能源系统的仿真分析能力。
附图说明
通过以下参照附图对本发明实施例的描述，本发明的上述以及其它目的、特征和优点将更为清楚，在附图中：
图1为本发明提供的多能互补综合能源系统动态潮流计算流程(以电热耦合系统为例)；
图2为本发明提供的非电系统基于面向对象程序架构的类及关系；
图3为本发明提供的主程序接入自定义模型情况下的计算流程。
图4为本发明实施例提供的6母线电力系统拓扑结构图；
图5为本发明实施例提供的冷网系统拓扑结构图；
图6为本发明实施例提供的热网系统拓扑结构图；
图7为本发明实施例提供的混合迭代求解算法效果图；
图8为cchp2的节点温度变化趋势图；
图9为热网节点温度分布图；
图10为cchp2的输出电功率变化趋势图；
图11为cchp1、2的输出热功率变化趋势图。
具体实施方式
以下基于实施例对本发明进行描述，但是本发明并不仅仅限于这些实施例。
本发明提供了一种用于多能互补综合能源系统的动态潮流分析方法，通过对由热、气、冷等非电系统和电力系统构建的多能耦合系统的动态潮流计算，实现多能耦合系统的网络分析。具体包含如下步骤：
步骤一：构建非电系统准稳态模型
为便于描述，本发明动态潮流分析方法实施例仅以热网为例描述本发明提供的非电系统准稳态模型，气、冷网等非电系统与热网结构相似，本领域技术人员容易将本发明提供的非电系统准稳态模型扩展至气网、冷网等其他非电系统。
本发明所述的非电系统准稳态模型指网络系统中的节点采用稳态方程，而管道采用针对温度、热量状态量的动态方程。
热力管网主要由热力管道和连接件构成，本发明基于图论将热力系统抽象成图进行网络拓扑分析，其中热力管道对应图的有向边，连接件(热源、热负荷、管道连接件)对应图的顶点，将阀门作为管道的附属属性，每个管段定义流量正方向，例如取为热网设计时的流体流向，从而可用图的矩阵描述热力系统拓扑结构。设热力系统拓扑结构共包含n个节点及b条管道，热网抽象后，可得到热网关联矩阵，上关联矩阵，下关联矩阵和回路矩阵。
(1)静态潮流模型
分别对热网进行流体力学建模和热力工况建模，分别构建热力系统流体力学稳态方程(公式1第1-3表达式)和热力工况稳态方程(公式1第4-5表达式)，两者联立获得热网静态潮流模型，如下所示。
其中，a为热网关联矩阵，
其中表达式1、5用于计算非电系统节点，表达式2-4管道用于计算管道公式(2)动态潮流模型第0时刻起始值。
(2)动态潮流模型
针对温度、热量状态量的动态潮流模型如下所示：
a为热网关联矩阵，
步骤二：构建电力系统潮流模型
本发明采用的电力系统潮流模型为现有技术，当今电力系统潮流模型及计算比较成熟，并出现很多成熟的电力系统仿真软件，比如中国电科院开发的psasp软件、美国epri开发bpa软件等，因此在此不再赘述，仅列出常用电力系统潮流模型公式。
其中，pi、qi分别为电力系统网络节点i的注入有功、无功功率，gij、bij分别为节点i、j之间的互电导和互电纳，δij表示节点i、j之间的相角差，u为节点电压。以上均为标幺值。
步骤三：构建非电系统与电力系统之间进行耦合的设备功率输出外特性稳态或一阶惯性动态模型
非电系统与电力系统之间进行耦合的耦合设备模型构建是本发明进行潮流统一计算的前提和基础，以电热耦合为例，电热耦合器件通常包括联供机组、水泵、热泵、电锅炉、制冷机等。这些设备内部的工作机制各自不同且较为复杂，但对于电网和非电系统耦合潮流计算而言，只需用功率输出外特性稳态或一阶惯性动态模型表征各耦合设备，便可建立电网潮流与非电系统潮流间的关系，实现联合求解。模型如下公式所示，本模型涵盖了电与非电所有耦合设备的描述，在计算中根据实际耦合器件的选择进行相应的模型计算。
输出外特性稳态模型如公式(4)所示：
一阶惯性动态模型如公式(5)所示：
一阶惯性动态模型主要考虑了实际物理系统中，自变量发生变化到因变量响应的时延效应(过程)，是常用等值模拟。其中τ是时间常数，s为拉普拉斯算子，该模型是用传递函数表示的。
公式(4)和(5)中，pg*、qg*分别为联供机组电、热功率，pp*为泵电功率，η为泵效率，m*为泵内流体流量，hp*为泵扬程，qhp*、qeb*、qc*分别为热泵、电锅炉、制冷机的热功率，php*、peb*、pc*分别为热泵、电锅炉、制冷机的电功率。以上均为标幺值。
步骤四：构建多能互补综合能源系统的动态潮流模型
联立步骤一非电系统潮流模型、步骤二电力系统潮流模型以及步骤三耦合设备运行外特性模型，构建多能互补综合能源系统的动态潮流模型，如下所示。
其中，fe＝0表示上述电力系统潮流模型，fh＝0表示非电系统准稳态模型，feh＝0表示耦合设备运行外特性模型，f＝0表示多能互补综合能源系统动态潮流模型。
步骤五：采用混合求解算法求解上述多能互补综合能源系统动态潮流模型
电力系统与非电系统，二者在系统构成、运行机理、响应时间尺度、输出特性等方面有着显著的差别，实际运行过程中的系统控制相对独立。以电热耦合为例，电、热两部分通过电热耦合设备关联，二者是一种松耦合的关系。因此根据电力系统与非电系统耦合的这个特征，在具体计算过程中，不需要在一次迭代中对两部分进行联合求解，可以单独进行，仅在每次电力或非电潮流收敛后，根据耦合器件的输出特性方程进行部分状态量的调整，再转入另一种能源潮流迭代计算，如此反复，达到整体收敛。这样做可以避免将电力系统与非电系统作为整体考虑计算收敛性差、耗时长和占用较多计算机内存空间的问题，且使得计算机程序易于设计实现。
进一步的，本发明针对非电系统和电力系统分别采用不同迭代方法进行混合求解算法，比如电力系统潮流计算采用最优乘子法，非电系统潮流计算采用牛顿-拉夫逊法，能够实现更快速的收敛。
本发明提供的多能互补综合能源系统动态潮流计算流程如图1所示，以电热耦合系统为例，进行说明。具体包含如下步骤：
(1)对多能互补综合能源系统拓扑关系进行分析，获取关联矩阵(如步骤一热网关联矩阵a，上关联矩阵au，下关联矩阵ad和回路矩阵bf等)；
(2)对模型计算数据进行数据标幺化处理，从而完成计算数据准备；
(3)初始化迭代次数k＝0，初始化耦合设备电力系统状态量xe(k)＝[pt,qt,ut,δt]t(字符定义见步骤二)及热网系统状态量xh(k)＝[m*t,q*t,δh*t,tet,tnt,qj*t]t(字符定义见步骤一)为标幺值，且初始值随机分配。
(4)进行热力管网静态潮流计算，然后判断静态潮流计算是否收敛，如果收敛跳入步骤(5)；否则计算结束；
(5)根据热力管网静态潮流计算结果修正耦合器件电功率pe(k)＝[pg*t,pp*t,php*t,peb*t,pc*t]t(字符定义见步骤3)，记为pe(1)(k)，并利用耦合器件电功率pe(1)(k)进行电力系统潮流计算，然后判断是否达到收敛，若是，跳入步骤(6)；否则计算结束；
(6)记电力系统潮流计算收敛后耦合器件电功率pe(2)(k)，根据耦合器件电功率pe(2)(k)修正耦合器件热网状态量xh(k)；然后判断是否|pe(1)(k)-pe(2)(k)|＜ε或k＞kmax，其中ε，kmax取经验值。若k＞kmax，则计算结束；若|pe(1)(k)-pe(2)(k)|＜ε且k≤kmax，则跳入步骤(7)进行动态计算，且当前时刻t＝t+dt，其中dt表示一个时步的度量；否则，k＝k+1后返回步骤(4)；
(7)判断是否存在扰动，若存在，处理扰动信息，更新热网状态量xh(k)后，跳入步骤(8)；否则直接跳入步骤(8)；
(8)进行热力管网动态潮流计算，判断计算是否收敛，若是，跳入步骤(9)，否则计算结束；
(9)根据热力管网动态潮流计算结果修正耦合器件电功率pe(k)，记为pe(3)(k)，并利用耦合器件电功率pe(3)(k)进行电力系统潮流计算，然后判断电力系统潮流计算是否能够计算收敛，若是跳入步骤(10)，否则计算结束；
(10)记电力系统潮流计算收敛后耦合器件电功率pe(4)(k)，根据耦合器件电功率pe(4)(k)修正耦合器件热网状态量xh(k)；然后判断是否|pe(3)(k)-pe(4)(k)|＜ε或k＞kmax，其中ε，kmax取经验值；若是跳入步骤(11)；否则k＝k+1后返回步骤(8)；
(11)判断是否t>tend，其中tend表示截止时间，若是，计算结束，输出计算结果；否则，t＝t+dt跳入步骤(7)。
综上，本发明先进行静态潮流计算过程：首先进行步骤一非电系统静态潮流模型计算，计算收敛后，根据非电系统收敛计算结果通过步骤三耦合设备功率输出外特性稳态或一阶惯性动态模型修正耦合器件电功率，然后进行电力系统潮流计算，收敛后，再根据电力系统收敛计算后耦合器件电功率通过步骤三耦合设备功率输出外特性稳态或一阶惯性动态模型修正电热耦合器件非电状态量。若前后两次计算得到的耦合器件电功率误差不满足精度要求，则刷新电热耦合器件电功率，继续进行计算，则满足，则判断是否已经超出最大迭代次数。若超出最大迭代次数则计算结束，否则进入纯延迟动态计算。
纯延迟动态计算过程中，在消除扰动后，分别进行非电动态潮流计算和电力系统潮流计算，过程与上述静态潮流计算过程类似，并在耦合器件功率误差满足精度要求或达到最大迭代次数后，当时间到达截止时刻，计算结束，输出结果。
其中，电力系统潮流计算可由现有成熟的电力系统仿真软件，比如中国电科院开发的psasp软件、美国epri开发bpa软件等完成，因此在此不再赘述。
进一步的，本发明所述的非电系统可以为热\冷\气等非电能源的任一种能源形式，比如上述热网，实现电热耦合。
为便于数据交互，同时也为简化计算、提高收敛性及便于进行结果分析，本发明潮流计算如上述均采用标幺值，但本发明不仅限于采用标幺值，还包括本领域技术人员容易想到的其他单位制。
本发明还提供了一种用于多能互补综合能源系统的动态潮流分析系统，包括电力分析模块、非电力分析模块和耦合分析模块，其中非电力分析模块用于实现步骤一非电模块动态潮流模型的准稳态模型计算，电力分析模块用于实现步骤二电力模块潮流模型的潮流计算，耦合分析模块用于实现非电模块与电力模块之间耦合设备功率输出外特性稳态或一阶惯性动态模型的计算。
电力分析模块、非电力分析模块和耦合分析模块共同实现步骤四和步骤五多能互补综合能源系统的动态潮流模型的求解。具体的，从数据文件读取多能互补综合能源系统动态潮流模型参数，执行上述多能互补综合能源系统动态潮流计算流程。
本发明非电力分析模块基于c++面向对象程序架构的设计理念，如图2所示，包括管道(pipe)及节点(node)两个基类，在此基础上设计包含依赖于上述两基类的设备模型类，用于指代同一类设备，然后设计设备模型类的子类模型实例类，指代实际设备，设计模型实例类的子类模型实例集合类，实现各实际设备模型对象的集合，从而形成非电力分析模块。
本发明采用嵌入式程序开发方法，将非电力分析模块以及耦合分析模块的模型计算单独编译成为动态链接库文件，电力系统仿真分析软件通过其用户程序接口调用该文件以实现联合计算，此时可采用现有成熟仿真分析软件；也可以将现有电力系统仿真分析软件单独编译成可执行程序，通过嵌入到非电力分析系统以及耦合分析系统的模型计算程序中，执行该程序以完成联合计算。
通过上述采用调用动态链接库文件的方式，可以实现用户自定义建模外接，具有开放性和可扩展性，用户通过图形化界面建立设备的自定义传递函数数学模型，程序自动生成模型文件，动态链接库程序解析模型文件并内置算法，并通过与主程序的接口实现数据的相互传递与协同计算。
以电力系统仿真分析软件通过其用户程序接口调用动态链接库文件实现联合计算为例，比如psasp电力系统仿真分析，本发明可以利用psasp用户程序接口(upi)功能，实现psasp与用户自定义非电力分析系统以及耦合分析系统的模型计算程序模块兼容，使psasp变成一个开放的软件包。其中的潮流计算用户程序接口(lf/upi)即实现电力潮流计算模块和上述用户程序模块交替运行，共同完成一项基于潮流计算的新的任务。
如图3所示，实现主程序接入用户自定义模型的接入程序计算流程包括：
(1)记时间t＝t，其中t为当前时刻，主程序完成自身动态潮流计算流程，并在判断已经达到收敛后，通过用户接口调用用户自定义模型，跳转至步骤(2)，若不能达到收敛，则程序结束；
(2)完成自定义模型的解析、赋初值及潮流计算流程，若潮流计算达到收敛，则跳转至步骤(3)，否则程序结束；
(3)判断收敛后是否满足程序结束约束条件，即将收敛后自定义模型的输出变量值和对应主程序接口变量值进行比较，若误差小于设定阈值则t＝t+dt，主程序继续进行下一时步动态潮流计算，否则程序结束。
本发明的实施例：
以6母线电力系统接冷网、热网为例，如图4所示，母线1接外部大电网，母线2、6分别接两个冷热电三联供机组(以下简称cchp机组)，母线3、4、5接负荷；此外，母线3可接风电场，母线4可接光伏电站。电力系统网络拓扑参数如表1和表2所示：
表1电力系统线路参数
表2变压器参数
如图5和图6所示，冷、热网均为5节点系统，其中节点5、4分别接图4两个cchp机组供冷、热，节点1、2、3接冷、热负荷。管道参数如表3所示(回水管道参数与供水管道对称)：
表3冷/热网管道参数
cchp机组热电比：
设cchp1热电比满足：
cchp2热电比满足：
其中φ1，φ2表示cchp机组热功率，p1，p2表示cchp机组电功率。
冷热比取cop＝1.2。
设负荷1为商用负荷，负荷2为工业负荷，负荷3为民用负荷。取负荷的功率因数为0.95，风电、光伏的功率因数为0.8。
将e1设为平衡节点，e3、e4、e5(即负荷1、负荷2、负荷3)设为pq节点，e2、e6(即cchp1、cchp2)设为pv-pq节点，根据上述参数，列出步骤二所表述的电力系统静态潮流模型，其中n＝6，gij、bij分别为表1、表2中各线路/变压器的电阻、电抗的倒数，即线路、变压器的电导和电纳标幺值，也即步骤四中的方程fe＝0。
然后根据表3中的管道首末节点连接关系，利用图论得到冷/热网关联矩阵，即步骤一a为热网关联矩阵，
阵，bf为回路矩阵，假设冷/热网系统无泄露，则q*＝0，假设各管道在同一水平面内，则z*＝0，h0*、sp根据泵的铭牌参数换算，为给定值；ta即耦合设备所在的当地环境温度，也为给定值，m*,δh*,te*,tn*,qj*为待求的冷/热网状态量，因此可以列出步骤一所表述的冷/热网静态潮流模型和动态潮流模型，动态潮流模型中，密度、管长、管截面积、比热容、导热率在传输工质确定的情况下均为给定值，即公式步骤四中的方程fh＝0。
本实施例中的电热冷耦合设备只有cchp1、cchp2，因此可以只计算公式11，根据上述cchp机组热电比公式和冷热比公式，得到步骤三中的公式11，即方程p*＝f(qj*)，其中，p*为耦合设备电力系统状态量xe中的设备电磁功率标幺值，qj*为耦合设备热网/冷网系统状态量xh中的设备热功率(热负荷，对于cchp为负值)标幺值，也即步骤四中的方程feh＝0。由此可以得到多能互补综合能源耦合系统的动态潮流模型。
再根据步骤五的计算流程进行求解。其中，电力系统潮流计算交与现有电力系统仿真软件计算；冷/热网系统潮流计算及耦合设备潮流计算根据客户自定义生成对应模型的动态链接库，由电力系统仿真分析软件通过其用户程序接口调用动态链接库文件实现联合计算，计算结果为xh以及xe的最终收敛数值解。
如图7为电热或电冷耦合时，一个时步内分别采用统一迭代算法和不同迭代算法得到的结果。比如电力系统计算采用最优乘子法，热网、冷网计算采用牛顿-拉夫逊法，从图中可以看出采用不同迭代算法的混合求解算法能够实现更快速的收敛。
本实施例的扰动信号是cchp2的节点温度在10s时刻由80度缓慢变化为100度，如图8所示，计算时间域为24分钟。以电热耦合系统为例，cchp2节点温度变化，带来的对热力系统、电力系统后续设备状态的影响，由于管道的传热延迟，其他节点的温度变化根据与该节点距离的不同，相继发生温度的变化，冷热电联供机组出力变化根据其节点温度的变化而变化的，因为温度上升所有出力随之增加，但因为燃料供应没有变化，最终还是回到了初始水平附近。如图9表示热网节点温度，图10表示cchp2的输出电功率变化曲线，图11表示cchp1、2的输出热功率变化曲线图。
以上所述仅为本发明的优选实施例，并不用于限制本发明，另外，本发明可以有各种改动和变化。凡在本发明的精神和原理之内所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。