一种基于双重并行计算的在线Prony分析方法技术

本发明专利技术公开了一种基于双重并行计算的在线Prony分析方法，涉及电力系统调度自动化领域和计算机高性能计算领域。该方法针对大电网发生低频振荡情况下，振荡参数辨识采用传统基于串行Prony算法资源利用率低、计算速度慢的弊端，提出了多计算节点的分布式并行，有效地进行任务调度和负载均衡，极大地降低系统的响应时间。首次实现了Prony数学模型上的并行化设计，采用多线程并行计算技术实现了Prony的多线程并行计算。本发明专利技术能够同步实现对电网多支路、多个电气量的振荡幅值、频率、初相位和衰减因子等参数的在线辨识，有效提高了计算分析速度，能够更好的适应大电网同步在线计算的要求。

【技术实现步骤摘要】
一种基于双重并行计算的在线Prony分析方法
本专利技术涉及电力系统调度自动化领域和计算机高性能计算领域，更具体是提供一种基于双重并行计算的在线Prony分析方法，在电网发生低频振荡过程中，实现对电网多支路、多个电气量的同步在线参数辨识。
技术介绍
电力系统在扰动下会发生电机转子间的相对摇摆并在缺乏阻尼的情况下引起持续振荡，振荡频率范围在0.1～2.5Hz，故称为低频振荡。低频振荡问题属于小扰动稳定范畴，随着互联电力系统的规模扩大、远距离重负荷输电系统的投入运行、快速自动励磁调节器和快速励磁系统的应用，国内外不少电力系统出现了低频振荡问题，低频振荡是影响电力系统的安全稳定运行的重要因素之一。对于由外界持续周期性功率扰动引起的系统低频振荡即强迫功率振荡，最有效的措施是快速定位并切除扰动源，为此专利《一种区域互联电网强迫功率振荡扰动源位置判断方法》(专利号：ZL201110390520.4)提出了“能流方向因子”的物理概念，采用Prony分析方法对支路的有功功率、起始节点频率数据进行参数辨识，求解各支路的能流方向因子，通过能流方向因子确定支路势能的大小及流动方向，并据此判断扰动源位置。广域测量系统(WAMS)由基于全球定位系统(GPS)的同步相量测量单元(PMU)及其通信系统组成，能够在广域电力系统中同步、高速采集机组和运行设备的有功功率、无功功率、电压、电流、相角以及重要的开关信号，是一种能对电力系统动态过程进行监测和分析的工具。广域测量系统为电网的频振荡的监测、振荡事件分析和振荡预防及抑制等几个方面提供了新的技术手段。Prony算法是一种能够根据采样值直接估算出信号频率、衰减、幅值和初相位的分析方法。它针对等间距采样点，假设模型是由一系列的具有任意振幅、相位、频率和衰减因子的指数函数的线性组合。Prony算法的数学模型离散时间的函数形式：视为X(n)的近似，式中bi和zi假定为复数，k为迭代次数，N为采样点数即：式中：Ai为振幅，θi为相位，αi为衰减因子，fi为频率，Δt为采样间隔。Prony方法的关键是认识到的拟合是一常系数线性差分方程的齐次解。即：亦有：X(n)＝-[α1X(n-1)+α2X(n-2)+...+αkX(n-k)](4)其中：n＝k,k+1,......N-1计算流程根据上式中的数学思想，Prony计算的处理过程如下：1、最小二乘计算式(4)是一个差分方程，对参数αi进行最小二乘估计，使误差平均和为最小。为了简化书写，用xn＝x(n)来表示第n个采样值。为保证计算精度，取迭代次数k＝N/2,可得N-k个方程如下：简写为：H·α＝-h(6)由公式(5)可知总平方误差W1为：由W1分别对α1、α2、……、αk求导，从而得出Prony的法方程：其中式(8)中的各元素计算如下：2、求解线性方程组由式(8)(9)(10)的计算结果，利用高斯约旦消元法求解线性方程组。得到α＝R-1·r；α＝(α1,α2,...,αk)T；(11)3、解高次方程根式(4)转换为x(n)+α1x(n-1)+α2x(n-2)+...+αkx(n-k)＝0(12)将式(1)代入式(12)得到特征多项式如下：zk+α1zk-1+α2zk-2+...+αk＝0(12)解高次方程根的过程可以转化为求解如下矩阵的特征值过程：求解矩阵特征值得到z＝(z1、z2、……、zk)T。.4、矩阵乘法及求解方程组这部分是Prony算法中最耗时的计算过程，为了确定待定系数b1、b2、…、bk,将z代入式(1)，得到方程组如下：式(14)左右两侧矩阵和左侧矩阵的转置矩阵做乘法运算，然后进行高斯消元法求出b：b＝(b1、b2、...、bk)T；5、结果输出根据以上结果可以算出幅值Ai、相位θi、频率fi和衰减因子αi。Prony的并行化：本专利技术利用的单机并行化工具是OpenMP，运用OpenMP的两种并行方式：隐式任务并行和显式任务并行进行Prony的并行化改造。隐式任务并行：隐式任务是由OpenMP提供的for和Section等编译指导语句形成的隐式并行区域生成的任务，隐式任务并行是对隐式任务指定调度方式、任务划分方式和数据控制方式并并行执行。隐式任务并行只需指定需要并行的部分，即可由OpenMP自动完成并行化，并在每次计算完成后自动汇总结果。隐式任务并行主要适用于结构规则的程序。这一类程序中，无数据相关性且划分后子任务的大小差异不大，如规则的for循环等。对此类程序，隐式任务并行能达到很高的效率。显式任务并行：显式任务是使用OpenMP的task指令指定的任务。task指令定义了与任务及其数据环境关联的代码。显式任务并行中，只要线程遇到任务构造(指定为task的部分)，就会生成新任务，并将新任务放置到OpenMP的任务共享队列(work-sharing)中。work-sharing队列是OpenMP内部队列的实现方式，各个核心(线程)共享这个队列，并依次从任务队列中取任务来并行计算。显示并行下，需要在每次并行计算完成后显式汇总结果。显式任务并行主要适用于不规则程序的并行化。在不规则程序中，数据可能存在相关性，划分得到的子任务差异可能较大，此时适宜采用显式任务并行进行程序并行化设计。静态调度：OpenMP中使用schedule调度子句时，schedule(static,size)指定的调度方式，当编译指导语句没有带schedule调度子句时，大部分系统中默认采用static调度方式。这种调度方式非常简单，静态调度可以使用参数size也可以不使用参数size，有n次循环迭代，t个线程时，使用参数size分配给每个线程size次连续的迭代计算，不使用size分配给每个线程n/t次连续的迭代计算。动态调度：OpenMP中使用schedule调度子句时，schedule(dynamic,size)指定的调度方式，动态调度依赖于运行时的状态动态确定线程所执行的迭代，由于线程启动和执行完的时间不确定，所以迭代被分配到哪个线程是无法事先知道的。动态调度可以使用参数size也可以不使用参数size，不指定size参数时是将迭代逐个地分配到各个线程，使用size参数时，每次分配给线程的迭代次数为指定的size次。伪共享现象：伪共享是指几个在逻辑上并不包含在同一个内存单元内的数据，由于被CPU加载在同一个缓存行当中，当在多线程环境下，被不同的核心执行，导致缓存行失效而引起的大量的缓存命中率降低的现象。在专利《一种区域互联电网强迫功率振荡扰动源位置判断方法》(专利号：ZL201110390520.4)所述方法中，计算量最大的步骤在于通过Prony方法实现能流方向因子的参数辨识。由于能流方向因子计算需要对电网中的多条振荡支路的有功功率和频率分量同时进行Prony参数辨识计算，当前主要是基于串行Prony计算方法，随着电网规模的扩大采用传统的串行计算方法计算速度不能满足在线分析的要求，需要采用并行计算。其次，由于Prony算法的上下文关联性较高，无法在结构上将计算任务划分并行，属于不规则的并行结构，因此在并行计算中需要基于Prony算法的数学模型进一步优化其并行结构，提高算法的并行效率，为此提出以下专利技术方法。
技术实现思路
本专利技术提供了一种基于双重并行计算的在本文档来自技高网...

【技术保护点】
一种基于双重并行计算的在线Prony分析方法，其特征在于，该方法包含下列步骤：(a)在由I条支路和J个节点组成的交流互联电网中，各支路和所有节点均装设同步测量单元PMU；计算平台由一个客户端、一个管理节点和C个计算节点组成的机群构成，其中Ci表示第i个计算节点，Pi为Ci的核数，Cij表示Ci中的第j个核，i＝1，2，…，C，j＝1，2，3，…，Pi；(b)当电网中发生低频振荡时，将PMU采集的M条支路E个电气量包括：有功功率ΔP、频率Δω，取同一时间段采样数据为一组发送给计算平台进行Prony分析，则需要分析的任务数为N＝EM，M＝1，2，…，I；每一类电气量中包含Len个采样数据，500≤Len≤3000；(c)多机并行计算过程，首先，进行计算节点能力测试，即客户端将接收到的N个任务发送给管理节点，管理节点从N个任务中任意选取一个任务，同时发送给Ci节点上进行计算，i＝1，2，…，C，记录每个节点对该任务的执行时间，节点Ci的执行时间记为Ti，根据第i节点的执行时间与C个节点总执行时间的比例对各计算节点的计算能力进行判别，根据计算能力的强弱对各计算节点进行排序，作为任务分配的依据，然后，将N个计算任务分配到各计算节点开展多机并行计算，节点Ci分配到的任务数为Ni：T=Σi=1CTi]]>Ni=[TiTN],i=1,2,...,C]]>将Ni个任务发送到计算节点Ci上计算，未能整除余下的任务根据计算能力排序依次分配到各计算节点，由Ci进行单节点多线程并行计算；(d)单节点多核并行计算过程：单节点多线程的并行化采用OpenMP的隐式任务并行和显式任务并行技术实现；当节点Ci接收到的任务数为Ni，每个任务的数据长度为Len，第m个采样数据为Xm，m＝1，2，…，Len；Prony计算过程中取迭代次数k＝Len/2；此步骤包含以下计算过程：(d1)最小二乘计算元素求解Ru,v=Σm=km=N-1Xm-vXm-u,u,v=1,2,......,k]]>ru=-Σm=km=N-1XmXm-v]]>采用隐式任务并行方式对Ru,v、ru进行并行计算，采用静态调度方式进行子任务的分配与调度，将Ru,v、ru平均划分为Pi个子任务，每个子任务负责[k/Pi]×k个数据；由OpenMP将划分得到的Pi个子任务自动分配到核Cij上执行，i＝1,2,…，C；j＝1，2，3，…，Pi；并行计算后得到Ru,v和ru，u＝1,2,…，k；v＝1，2,…，k；(d2)采用高斯约旦消元法求解线性方程组：R1,1R1,2R1,3...R1,kR2,1R2,2R2,3...R2,kR3,1R3,2R3,3...R3,k...............Rk,1Rk,2Rk,3...Rk,k·α1α2α3...αk=r1r2r3...rk]]>得到α＝(α1、α2、...、αk)T采用隐式并行策略实现高斯约旦消元法求解线性方程组的并行化；高斯约旦消元法要求进行k次消元，每次消元过程需要对k个方程进行消元，将k个方程的消元过程分配到Pi个核上并行执行，每个核分配[k/Pi]个方程的消元；划分时，使任务划分的数据间隔为一个高速缓存行的整数倍大小，防止伪共享现象，并在调度上采用动态调度方式。划分后的子任务将由OpenMP自动调度到核Ci1、Ci2、…、CiPi上计算；计算后得到α1、α2、…、αk；(d3)将α1、α2、…、αk代入式X(n)＝‑[α1X(n‑1)+α2X(n‑2)+...+αkX(n‑k)]求解高次方程zk+α1zk‑1+α2zk‑2+...+αk＝0即求解如下矩阵的特征值：-α1-α2...-αk-1-αk10...00...............0010000010]]>采用显式任务并行方式实现矩阵特征值的并行求解；求解过程需要k次循环，每次循环过程分为两个部分的计算，一是与其他任务有数据相关的部分，即对矩阵元素更新部分，二是无数据相关的部分，即计算出特征值部分；主线程只负责执行有数据相关的部分，无数据相关的部分利用task指令显式构造并行任务，添加到任务队列中，在OpenMP的调度下由其他线程执行。核Ci1执行主线程，在其对元素更新之后，将其余计算部分显式构造成子任务添加到任务队列中，进而由空闲的核取子任务并执行；最后显式汇总求得：z1、z2、……、zk。(d4)求解线性方程组111...1z1z2z3...zkz12z22z32...zk2...............z1N-1z2N-1z3N-1...zkN-1·b1b2b3...bk=x1x2x3...xN-1]]>采用隐...

【技术特征摘要】
1.一种基于双重并行计算的在线Prony分析方法，其特征在于，该方法包含下列步骤：(a)在由I条支路和J个节点组成的交流互联电网中，各支路和所有节点均装设同步测量单元PMU；计算平台由一个客户端、一个管理节点和C个计算节点组成的机群构成，其中Ci表示第i个计算节点，Pi为Ci的核数，Cij表示Ci中的第j个核，i＝1，2，…，C，j＝1，2，3，…，Pi；(b)当电网中发生低频振荡时，PMU采集的M条支路E个电气量包括：有功功率ΔP、频率Δω，取同一时间段采样数据为一组发送给计算平台进行Prony分析，则需要分析的任务数为N＝EM，M＝1，2，…，I；每一类电气量中包含Len个采样数据，500≤Len≤3000；(c)多机并行计算过程，首先，进行计算节点能力测试，即客户端将接收到的N个任务发送给管理节点，管理节点从N个任务中任意选取一个任务，同时发送给Ci节点上进行计算，i＝1，2，…，C，记录每个节点对该任务的执行时间，节点Ci的执行时间记为Ti，根据第i节点的执行时间与C个节点总执行时间的比例对各计算节点的计算能力进行判别，根据计算能力的强弱对各计算节点进行排序，作为任务分配的依据，然后，将N个计算任务分配到各计算节点开展多机并行计算，节点Ci分配到的任务数为Ni：将Ni个任务发送到计算节点Ci上计算，未能整除余下的任务根据计算能力排序依次分配到各计算节点，由Ci进行单节点多线程并行计算；(d)单节点多核并行计算过程：单节点多线程的并行化采用OpenMP的隐式任务并行和显式任务并行技术实现；当节点Ci接收到的任务数为Ni，每个任务的数据长度为Len，第m个采样数据为Xm，m＝1，2，…，Len；Prony计算过程中取迭代次数k＝Len/2；此步骤包含以下计算过程：(d1)最小二乘计算元素求解采用隐式任务并行方式对Ru,v、rv进行并行计算，采用静态调度方式进行子任务的分配与调度，将Ru,v、rv平均划分为Pi个子任务，每个子任务负责[k/Pi]×k个数据；由OpenMP将划分得到的Pi个子任务自动分配到核Cij上执行，i＝1,2,…，C；j＝1，2，3，…，Pi；并行计算后得到Ru,v和rv，u＝1,2,…，k；v＝1，2,…，k；(d2)采用高斯约旦消元法求解线性方程组：1得到α＝(α1、α2、...、αk)T采用隐...

【专利技术属性】
技术研发人员：杨东俊，王多强，许汉平，徐敬友，黄家祺，
申请(专利权)人：国家电网公司，国网湖北省电力公司经济技术研究院，
类型：发明
国别省市：北京;11