CN104504257A - 一种基于双重并行计算的在线Prony分析方法 - Google Patents

Abstract

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

Description

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

技术领域

本发明涉及电力系统调度自动化领域和计算机高性能计算领域，更具体是提供一种基于双重并行计算的在线Prony分析方法，在电网发生低频振荡过程中，实现对电网多支路、多个电气量的同步在线参数辨识。

背景技术

电力系统在扰动下会发生电机转子间的相对摇摆并在缺乏阻尼的情况下引起持续振荡，振荡频率范围在0.1～2.5Hz，故称为低频振荡。低频振荡问题属于小扰动稳定范畴，随着互联电力系统的规模扩大、远距离重负荷输电系统的投入运行、快速自动励磁调节器和快速励磁系统的应用，国内外不少电力系统出现了低频振荡问题，低频振荡是影响电力系统的安全稳定运行的重要因素之一。

对于由外界持续周期性功率扰动引起的系统低频振荡即强迫功率振荡，最有效的措施是快速定位并切除扰动源，为此专利《一种区域互联电网强迫功率振荡扰动源位置判断方法》(专利号：ZL201110390520.4)提出了“能流方向因子”的物理概念，采用Prony分析方法对支路的有功功率、起始节点频率数据进行参数辨识，求解各支路的能流方向因子，通过能流方向因子确定支路势能的大小及流动方向，并据此判断扰动源位置。

广域测量系统(WAMS)由基于全球定位系统(GPS)的同步相量测量单元(PMU)及其通信系统组成，能够在广域电力系统中同步、高速采集机组和运行设备的有功功率、无功功率、电压、电流、相角以及重要的开关信号，是一种能对电力系统动态过程进行监测和分析的工具。广域测量系统为电网的频振荡的监测、振荡事件分析和振荡预防及抑制等几个方面提供了新的技术手段。

Prony算法是一种能够根据采样值直接估算出信号频率、衰减、幅值和初相位的分析方法。它针对等间距采样点，假设模型是由一系列的具有任意振幅、相位、频率和衰减因子的指数函数的线性组合。

Prony算法的数学模型

离散时间的函数形式：

$\hat{X}\left(n\right)=\underset{i=1}{\overset{k}{\mathrm{\Σ}}}{b}_i{z}_i^n,n=\mathrm{0,1,2}...,N-1---\left(1\right)$

视(n)为X(n)的近似，式中b_i和z_i假定为复数，k为迭代次数，N为采样点数即：

$\left\{\begin{array}{c}{b}_i=A_i\exp (j,{\mathrm{\θ}}_i)\\ z_i=\exp \left[({\mathrm{\α}}_i+j2\mathrm{\π}{f}_i)\mathrm{\Δt}\right]\end{array}\right.---\left(2\right)$

式中：A_i为振幅，θ_i为相位，α_i为衰减因子，f_i为频率，Δt为采样间隔。

Prony方法的关键是认识到(n)的拟合是一常系数线性差分方程的齐次解。即：

$\hat{X}\left(n\right)=-\underset{i=0}{\overset{k}{\mathrm{\Σ}}}{\mathrm{\α}}_i\hat{X}(n-i),n=k,k+1,......N-1;---\left(3\right)$

亦有：

X(n)＝-[α₁X(n-1)+α₂X(n-2)+...+α_kX(n-k)]    (4)

其中：n＝k,k+1,......N-1

计算流程

根据上式中的数学思想，Prony计算的处理过程如下：

1、最小二乘计算

式(4)是一个差分方程，对参数α_i进行最小二乘估计，使误差平均和为最小。为了简化书写，用x_n＝x(n)来表示第n个采样值。为保证计算精度，取迭代次数k＝N/2,可得

N-k个方程如下：

$\left[\begin{array}{ccccc}{X}_{k-1}& X_{k-2}& ...& X_1& X_0\\ X_k& X_{k-1}& ...& X_2& X_1\\ ...& ...& ...& ...& \\ X_{N-2}& X_{N-3}& ...& X_{N-k}& X_{N-1-k}\end{array}\right]\·\left[\begin{array}{c}{\mathrm{\α}}_1\\ {\mathrm{\α}}_2\\ ...\\ {\mathrm{\α}}_k\end{array}\right]=-\left[\begin{array}{c}{X}_k\\ X_{k+1}\\ ...\\ X_{N-1}\end{array}\right]---\left(5\right)$

简写为：H·α＝-h    (6)

由公式(5)可知总平方误差W₁为：

$W_1=\underset{m=k}{\overset{m=N-1}{\mathrm{\Σ}}}{(x_m+a_1{x}_{m-1}+...+a_i{x}_{m-i}+......a_k{x}_{m-k})}^2---\left(7\right)$

由W₁分别对α₁、α₂、……、α_k求导，从而得出Prony的法方程：

$\left[\begin{array}{ccccc}{R}_{\mathrm{1,1}}& R_{\mathrm{1,2}}& R_{\mathrm{1,3}}& ...& R_{1,k}\\ R_{\mathrm{2,1}}& R_{\mathrm{2,2}}& R_{\mathrm{2,3}}& ...& R_{2,k}\\ R_{\mathrm{3,1}}& R_{\mathrm{3,2}}& R_{\mathrm{3,3}}& ...& R_{3,k}\\ ...& ...& ...& ...& ...\\ R_{k,1}& R_{k,2}& R_{k,3}& ...& R_{k,k}\end{array}\right]\·\left[\begin{array}{c}{\mathrm{\α}}_1\\ {\mathrm{\α}}_2\\ {\mathrm{\α}}_3\\ ...\\ {\mathrm{\α}}_k\end{array}\right]=\left[\begin{array}{c}{r}_1\\ r_2\\ r_3\\ ...\\ r_k\end{array}\right]---\left(8\right)$

其中式(8)中的各元素计算如下：

$R_{i,j}=\underset{m=k}{\overset{m=N-1}{\mathrm{\Σ}}}{x}_{m-j}{x}_{m-i};i,j=\mathrm{1,2},......,k;---\left(9\right)$

$r_i=-\underset{m=k}{\overset{m=N-1}{\mathrm{\Σ}}}{x}_m{x}_{m-i};---\left(10\right)$

2、求解线性方程组

由式(8)(9)(10)的计算结果，利用高斯约旦消元法求解线性方程组。得到

α＝R^-1·r；α＝(α₁,α₂,...,α_k)^T；    (11)

3、解高次方程根

式(4)转换为

x(n)+α₁x(n-1)+α₂x(n-2)+...+α_kx(n-k)＝0    (12)

将式(1)代入式(12)得到特征多项式如下：

z^k+α₁z^k-1+α₂z^k-2+...+α_k＝0    (12)

解高次方程根的过程可以转化为求解如下矩阵的特征值过程：

$\left[\begin{array}{ccccc}-{\mathrm{\α}}_1& -{\mathrm{\α}}_2& ...& -{\mathrm{\α}}_{k-1}& -{\mathrm{\α}}_k\\ 1& 0& ...& 0& 0\\ ...& ...& ...& ...& ...\\ 0& 0& 1& 0& 0\\ 0& 0& 0& 1& 0\end{array}\right]---\left(13\right)$

求解矩阵特征值得到z＝(z₁、z₂、……、z_k)^T。.

4、矩阵乘法及求解方程组

这部分是Prony算法中最耗时的计算过程，为了确定待定系数b₁、b₂、…、b_k,将z代入式(1)，得到方程组如下：

$\left[\begin{array}{ccccc}1& 1& 1& ...& 1\\ z_1& z_2& z_3& ...& z_k\\ z_1^2& z_2^2& z_3^2& ...& z_k^2\\ ...& ...& ...& ...& ...\\ z_1^{N-1}& z_2^{N-1}& z_3^{N-1}& ...& z_k^{N-1}\end{array}\right]\·\left[\begin{array}{c}{b}_1\\ b_2\\ b_3\\ ...\\ b_k\end{array}\right]=\left[\begin{array}{c}{x}_1\\ x_2\\ x_3\\ ...\\ x_{N-1}\end{array}\right]---\left(14\right)$

式(14)左右两侧矩阵和左侧矩阵的转置矩阵做乘法运算，然后进行高斯消元法求出b：

b＝(b₁、b₂、...、b_k)^T；

5、结果输出

根据以上结果可以算出幅值A_i、相位θ_i、频率f_i和衰减因子α_i。

$\left\{\begin{array}{c}{A}_i=\left|b_i\right|\\ {\mathrm{\θ}}_i=\arctan (\mathrm{lm}\left(b_i\right)/\mathrm{Re}\left(b_i\right))/2\mathrm{\π\Δt}\\ {\mathrm{\α}}_i=\mathrm{lm}\left(z_i\right)/\mathrm{\Δt}\\ f_i=\arctan (\mathrm{lm}\left(z_i\right)/\mathrm{Re}\left(z_i\right))/2\mathrm{\π\Δt}\end{array}\right.i=\mathrm{1,2},...,k;---\left(15\right)$

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分析方法，该方法针对大电网发生低频振荡情况下，振荡参数辨识采用传统基于串行Prony算法资源利用率低、计算速度慢的弊端，提出了Prony的双重并行的计算方法，该方法能够实现多机上分布式并行和单机多核并行计算。双重并行计算提高了Prony方法的计算效率，能够同步实现对电网多支路、多个电气量的振荡幅值、频率、初相位和衰减因子等参数的在线辨识，为大电网支路能流方向因子计算提供了有效的参数在线辨识方法。

为了实现上述目的，本发明首先构建了由多机构成的分布式计算平台。该分布式计算平台中的节点类型分为两种：管理节点Master和计算节点Slave。Master节点包括三项功能：计算时间检测，即在任务分配前对当前分布式计算环境中各节点的计算时间进行校验；节点状态监测，即对Slave的当前状态进行监控，防止节点因意外宕机；负载均衡，Master动态调整节点任务量，使每个节点的任务量大致相等。Slave节点从Master节点接收任务并完成计算，将结果回送给Master节点并等待下一个任务的到来。

其次，本发明针对多线程并行计算模式完成了Prony数学模型的并行化，实现了单节点内多核计算平台上的Prony并行计算。本发明将Prony算法的最小二乘计算、求解线性方程组、解高次方程根、矩阵乘法及求解方程组等各步骤分别进行并行化处理，使核间任务达到负载均衡，提高计算节点资源利用率和计算效率。

基于上述设计，最终实现了一种基于双重并行计算的在线Prony分析方法。

一种基于双重并行计算的在线Prony分析方法，该方法包含下列步骤：

(a)在由I条支路和J个节点组成的交流互联电网中，各支路和所有节点均装设同步测量单元PMU；计算平台由一个客户端、一个管理节点和C个计算节点组成的机群构成，其中C_i表示第i个计算节点，P_i为C_i的核数，C_ij表示C_i中的第j个核，i＝1，2，…，C，j＝1，2，3，…，P_i；

(b)当电网中发生低频振荡时，将PMU采集的M条支路E个电气量包括：有功功率ΔP、频率Δω，取同一时间段采样数据为一组发送给计算平台进行Prony分析，则需要分析的任务数为N＝EM，M＝1，2，…，I；每一类电气量中包含Len个采样数据，500≤Len≤3000；

(c)多机并行计算过程，首先，进行计算节点能力测试，即客户端将接收到的N个任务发送给管理节点，管理节点从N个任务中任意选取一个任务，同时发送给C_i节点上进行计算，i＝1，2，…，C，记录每个节点对该任务的执行时间，节点C_i的执行时间记为T_i，根据第i节点的执行时间与C个节点总执行时间的比例对各计算节点的计算能力进行判别，根据计算能力的强弱对各计算节点进行排序，作为任务分配的依据。然后，将N个计算任务分配到各计算节点开展多机并行计算，节点C_i分配到的任务数为N_i：

$T=\underset{i=1}{\overset{c}{\mathrm{\Σ}}}{T}_{\text{i}}$

$N_i=\left[\frac{T_i}{T}N\right],i=\mathrm{1,2},...,C$

将N_i个任务发送到计算节点C_i上计算，未能整除余下的任务根据计算能力排序依次分配到各计算节点，由C_i进行单节点多线程并行计算；

(d)单节点多核并行计算过程：单节点多线程的并行化采用OpenMP的隐式任务并行和显式任务并行技术实现；当节点C_i接收到的任务数为N_i，每个任务的数据长度为Len，第m个采样数据为X_m，m＝1，2，…，Len；Prony计算过程中取迭代次数k＝Len/2；

此步骤包含以下计算过程：

(d1)最小二乘计算元素求解

$R_{u,v}=\underset{m=k}{\overset{m=N-1}{\mathrm{\Σ}}}{X}_{m-v}{X}_{m-u},u,v=\mathrm{1,2},......,k$

$r_u=-\underset{m=k}{\overset{m=N-1}{\mathrm{\Σ}}}{X}_m{X}_{m-v}$

采用隐式任务并行方式对R_u,v、r_u进行并行计算，采用静态调度方式进行子任务的分配与调度，将R_u,v、r_u平均划分为P_i个子任务，每个子任务负责[k/P_i]×k个数据；由OpenMP将划分得到的P_i个子任务自动分配到核C_ij上执行，i＝1,2,…，C；j＝1，2，3，…，P_i；并行计算后得到R_u,v和r_u，u＝1,2,…，k；v＝1，2,…，k；

(d2)采用高斯约旦消元法求解线性方程组：

$\left[\begin{array}{ccccc}{R}_{\mathrm{1,1}}& R_{\mathrm{1,2}}& R_{\mathrm{1,3}}& ...& R_{1,k}\\ R_{\mathrm{2,1}}& R_{\mathrm{2,2}}& R_{\mathrm{2,3}}& ...& R_{2,k}\\ R_{\mathrm{3,1}}& R_{\mathrm{3,2}}& R_{\mathrm{3,3}}& ...& R_{3,k}\\ ...& ...& ...& ...& ...\\ R_{k,1}& R_{k,2}& R_{k,3}& ...& R_{k,k}\end{array}\right]\·\left[\begin{array}{c}{\mathrm{\α}}_1\\ {\mathrm{\α}}_2\\ {\mathrm{\α}}_3\\ ...\\ {\mathrm{\α}}_k\end{array}\right]=\left[\begin{array}{c}{r}_1\\ r_2\\ r_3\\ ...\\ r_k\end{array}\right]$

得到

α＝(α₁、α₂、...、α_k)^T

采用隐式并行策略实现高斯约旦消元法求解线性方程组的并行化；高斯约旦消元法要求进行k次消元，每次消元过程需要对k个方程进行消元，将k个方程的消元过程分配到P_i个核上并行执行，每个核分配[k/P_i]个方程的消元；划分时，使任务划分的数据间隔为一个高速缓存行的整数倍大小，防止伪共享现象，并在调度上采用动态调度方式。划分后的子任务将由OpenMP自动调度到核C_i1、C_i2、…、C_iPi上计算；计算后得到α₁、α₂、…、α_k；(d3)将α₁、α₂、…、α_k代入式

X(n)＝-[α₁X(n-1)+α₂X(n-2)+...+α_kX(n-k)]

求解高次方程

z^k+α₁z^k-1+α₂z^k-2+...+α_k＝0

即求解如下矩阵的特征值：

$\left[\begin{array}{ccccc}-{\mathrm{\α}}_1& -{\mathrm{\α}}_2& ...& -{\mathrm{\α}}_{k-1}& -{\mathrm{\α}}_k\\ 1& 0& ...& 0& 0\\ ...& ...& ...& ...& ...\\ 0& 0& 1& 0& 0\\ 0& 0& 0& 1& 0\end{array}\right]$

采用显式任务并行方式实现矩阵特征值的并行求解；求解过程需要k次循环，每次循环过程分为两个部分的计算，一是与其他任务有数据相关的部分，即对矩阵元素更新部分，二是无数据相关的部分，即计算出特征值部分；主线程只负责执行有数据相关的部分，无数据相关的部分利用task指令显式构造并行任务，添加到任务队列中，在OpenMP的调度下由其他线程执行。核Ci1执行主线程，在其对元素更新之后，将其余计算部分显式构造成子任务添加到任务队列中，进而由空闲的核C_i2、...、C_iPi取子任务并执行；最后显式汇总求得：

z₁、z₂、……、z_k。

(d4)求解线性方程组

$\left[\begin{array}{ccccc}1& 1& 1& ...& 1\\ z_1& z_2& z_3& ...& z_k\\ z_1^2& z_2^2& z_3^2& ...& z_k^2\\ ...& ...& ...& ...& ...\\ z_1^{N-1}& z_2^{N-1}& z_3^{N-1}& ...& z_k^{N-1}\end{array}\right]\·\left[\begin{array}{c}{b}_1\\ b_2\\ b_3\\ ...\\ b_k\end{array}\right]=\left[\begin{array}{c}{x}_1\\ x_2\\ x_3\\ ...\\ x_{N-1}\end{array}\right]$

采用隐式任务并行方式完成并行计算，将转置矩阵分为P_i个较小的矩阵块，每个矩阵块大小为(k/P_i)×N，这样就产生P_i个矩阵相乘的子任务，指定调度方式为静态调度，P_i个子任务由OpenMP自动调度到核C_i1、C_i2、...、C_iPi上完成并行计算；并行计算后得到k个方程，用高斯约旦消元法和隐式任务并行方式进行线性方程组的并行求解，指定调度方式为动态调度；并行计算得到b₁、b₂、…、b_k；

(d5)将以上求出的z_u和b_u代入下式：

$\left\{\begin{array}{c}{A}_u=\left|b_u\right|\\ {\mathrm{\θ}}_u=\arctan (\mathrm{lm}\left(b_u\right)/\mathrm{Re}\left(b_u\right))/2\mathrm{\π\Δt}\\ {\mathrm{\α}}_u=\mathrm{lm}\left(z_u\right)/\mathrm{\Δt}\\ f_u=\arctan (\mathrm{lm}\left(z_u\right)/\mathrm{Re}\left(z_u\right))/2\mathrm{\π\Δt}\end{array}\right.uu=\mathrm{1,2},...,k$

计算出每个任务的幅值A_u、相位θ_u、频率f_u和衰减因子α_u，u＝1，2，…，k；

对分配到的N_i个计算任务重复步骤(d1)至(d5)，完成节点C_i上的N_i个计算任务，得到N_i个任务对应电气量的振幅、频率、初相位、衰减因子的值；

(e)在每个节点完成其负责的N_i个计算任务后，将各电气量的振幅、频率、初相位和衰减因子的值发送到管理节点；由管理节点汇总C个计算节点的结果，将结果返回给客户端。

本发明与现有技术相比，具有以下优点：

1.采用并行处理技术对Prony计算过程进行了全面的并行化改造。首次实现了Prony数学模型上的并行化设计，利用多核CPU计算平台，采用多线程并行计算技术实现了Prony的多核多线程并行计算，提高了Prony算法的计算速度。并行Prony计算过程在数据规模较大时能获得更高的加速比；同时随着核数的增多，并行Prony计算能达到更高的加速比，实施效果如图4所示。

2.针对大电网低频振荡情况下多支路同步在线计算需求，本发明设计了分布式计算平台，建立了两重并行结构的处理框架——计算节点间的分布式并行与计算节点内部的多线程并行。在分布式计算中，有效地进行任务调度和负载均衡，极大地降低系统的响应时间，提高了系统的处理速度和效率，并具有良好的可扩展性。

3.本发明设计的双重并行计算方式，以简单的接口对外提供服务，内部采用多种计算方法，但对用户层程序是完全透明的，因此使用方便、简单易用。同时，本发明所建立的Prony并行计算系统，结构简单、层次清晰，利用本发明定义的并行计算规则建立，克服了传统并行及分布式计算系统繁琐的约束，实现了更好的计算效率。

附图说明

图1为双重并行计算平台体系结构示意图。

图2为Prony计算的多机并行过程示意图。

图3为Prony计算的多核并行过程示意图。

图4为Prony计算的多核并行与串行比较的加速比曲线。

具体实施方式

本发明提供了一个基于双重并行计算的Prony分析方法，下面结合案例对本发明的多机和多核并行计算过程和实施效果做进一步说明。

实施例一

一种基于双重并行计算的在线Prony分析方法，该方法包含下列步骤：

(a)在由I条支路和J个节点组成的交流互联电网中，各支路和所有节点均装设同步测量单元PMU；计算平台结构如图1所示，由一个客户端、一个管理节点和C个计算节点组成的机群构成，其中C_i表示第i个计算节点，P_i为C_i的核数，C_ij表示C_i中的第j个核，i＝1，2，…，C，j＝1，2，3，…，P_i；

(b)当电网中发生低频振荡时，将PMU采集的M条支路E个电气量包括：有功功率ΔP、频率Δω，取同一时间段采样数据为一组发送给计算平台进行Prony分析，则需要分析的任务数为N＝EM，M＝1，2，…，I；每一类电气量中包含Len个采样数据，500≤Len≤3000；

(c)多机并行计算过程如图2所示。首先，进行计算节点能力测试，即客户端将接收到的N个任务发送给管理节点，管理节点从N个任务中任意选取一个任务，同时发送给C_i节点上进行计算，i＝1，2，…，C，记录每个节点对该任务的执行时间，节点C_i的执行时间记为T_i，根据第i节点的执行时间与C个节点总执行时间的比例对各计算节点的计算能力进行判别，根据计算能力的强弱对各计算节点进行排序，作为任务分配的依据。然后，将N个计算任务分配到各计算节点开展多机并行计算，节点C_i分配到的任务数为N_i：

$T=\underset{i=1}{\overset{c}{\mathrm{\Σ}}}{T}_{\text{i}}$

$N_i=\left[\frac{T_i}{T}N\right],i=\mathrm{1,2},...,C$

将N_i个任务发送到计算节点C_i上计算，未能整除余下的任务根据计算能力排序依次分配到各计算节点，由C_i进行单节点多线程并行计算；

(d)单节点多核并行计算过程如图3所示。单节点多线程的并行化采用OpenMP的隐式任务并行和显式任务并行技术实现；当节点C_i接收到的任务数为N_i，每个任务的数据长度为Len，第m个采样数据为X_m，m＝1，2，…，Len；Prony计算过程中取迭代次数k＝Len/2；

此步骤包含以下计算过程：

(d1)最小二乘计算元素求解

$R_{u,v}=\underset{m=k}{\overset{m=N-1}{\mathrm{\Σ}}}{X}_{m-v}{X}_{m-u},u,v=\mathrm{1,2},......,k$

$r_u=-\underset{m=k}{\overset{m=N-1}{\mathrm{\Σ}}}{X}_m{X}_{m-v}$

采用隐式任务并行方式对R_u,v、r_u进行并行计算，采用静态调度方式进行子任务的分配与调度，将R_u,v、r_u平均划分为P_i个子任务，每个子任务负责[k/P_i]×k个数据；由OpenMP将划分得到的P_i个子任务自动分配到核C_ij上执行，i＝1,2,…，C；j＝1，2，3，…，P_i；并行计算后得到R_u,v和r_u，u＝1,2,…，k；v＝1，2,…，k；

(d2)采用高斯约旦消元法求解线性方程组：

$\left[\begin{array}{ccccc}{R}_{\mathrm{1,1}}& R_{\mathrm{1,2}}& R_{\mathrm{1,3}}& ...& R_{1,k}\\ R_{\mathrm{2,1}}& R_{\mathrm{2,2}}& R_{\mathrm{2,3}}& ...& R_{2,k}\\ R_{\mathrm{3,1}}& R_{\mathrm{3,2}}& R_{\mathrm{3,3}}& ...& R_{3,k}\\ ...& ...& ...& ...& ...\\ R_{k,1}& R_{k,2}& R_{k,3}& ...& R_{k,k}\end{array}\right]\·\left[\begin{array}{c}{\mathrm{\α}}_1\\ {\mathrm{\α}}_2\\ {\mathrm{\α}}_3\\ ...\\ {\mathrm{\α}}_k\end{array}\right]=\left[\begin{array}{c}{r}_1\\ r_2\\ r_3\\ ...\\ r_k\end{array}\right]$

得到

α＝(α₁、α₂、...、α_k)^T

采用隐式并行策略实现高斯约旦消元法求解线性方程组的并行化；高斯约旦消元法要求进行k次消元，每次消元过程需要对k个方程进行消元，将k个方程的消元过程分配到P_i个核上并行执行，每个核分配[k/P_i]个方程的消元；划分时，使任务划分的数据间隔为一个高速缓存行的整数倍大小，防止伪共享现象，并在调度上采用动态调度方式。划分后的子任务将由OpenMP自动调度到核C_i1、C_i2、…、C_iPi上计算；计算后得到α₁、α₂、…、α_k；

(d3)将α₁、α₂、…、α_k代入式

X(n)＝-[α₁X(n-1)+α₂X(n-2)+...+α_kX(n-k)]

求解高次方程

z^k+α₁z^k-1+α₂z^k-2+...+α_k＝0

即求解如下矩阵的特征值：

$\left[\begin{array}{ccccc}-{\mathrm{\α}}_1& -{\mathrm{\α}}_2& ...& -{\mathrm{\α}}_{k-1}& -{\mathrm{\α}}_k\\ 1& 0& ...& 0& 0\\ ...& ...& ...& ...& ...\\ 0& 0& 1& 0& 0\\ 0& 0& 0& 1& 0\end{array}\right]$

采用显式任务并行方式实现矩阵特征值的并行求解；求解过程需要k次循环，每次循环过程分为两个部分的计算，一是与其他任务有数据相关的部分，即对矩阵元素更新部分，二是无数据相关的部分，即计算出特征值部分；主线程只负责执行有数据相关的部分，无数据相关的部分利用task指令显式构造并行任务，添加到任务队列中，在OpenMP的调度下由其他线程执行。核Ci1执行主线程，在其对元素更新之后，将其余计算部分显式构造成子任务添加到任务队列中，进而由空闲的核取子任务并执行；最后显式汇总求得：

z₁、z₂、……、z_k。

(d4)求解线性方程组

$\left[\begin{array}{ccccc}1& 1& 1& ...& 1\\ z_1& z_2& z_3& ...& z_k\\ z_1^2& z_2^2& z_3^2& ...& z_k^2\\ ...& ...& ...& ...& ...\\ z_1^{N-1}& z_2^{N-1}& z_3^{N-1}& ...& z_k^{N-1}\end{array}\right]\·\left[\begin{array}{c}{b}_1\\ b_2\\ b_3\\ ...\\ b_k\end{array}\right]=\left[\begin{array}{c}{x}_1\\ x_2\\ x_3\\ ...\\ x_{N-1}\end{array}\right]$

采用隐式任务并行方式完成并行计算，将转置矩阵分为P_i个较小的矩阵块，每个矩阵块大小为(k/P_i)×N，这样就产生P_i个矩阵相乘的子任务，指定调度方式为静态调度，P_i个子任务由OpenMP自动调度到核C_i1、C_i2、...、C_iPi上完成并行计算；并行计算后得到k个方程，用高斯约旦消元法和隐式任务并行方式进行线性方程组的并行求解，指定调度方式为动态调度；并行计算得到b₁、b₂、…、b_k；

(d5)将以上求出的z_u和b_u代入下式：

$\left\{\begin{array}{c}{A}_u=\left|b_u\right|\\ {\mathrm{\θ}}_u=\arctan (\mathrm{lm}\left(b_u\right)/\mathrm{Re}\left(b_u\right))/2\mathrm{\π\Δt}\\ {\mathrm{\α}}_u=\mathrm{lm}\left(z_u\right)/\mathrm{\Δt}\\ f_u=\arctan (\mathrm{lm}\left(z_u\right)/\mathrm{Re}\left(z_u\right))/2\mathrm{\π\Δt}\end{array}\right.uu=\mathrm{1,2},...,k$

计算出每个任务的幅值A_u、相位θ_u、频率f_u和衰减因子α_u，u＝1，2，…，k；

对分配到的N_i个计算任务重复步骤(d1)至(d5)，完成节点C_i上的N_i个计算任务，得到N_i个任务对应电气量的振幅、频率、初相位、衰减因子的值；

(e)在每个节点完成其负责的N_i个计算任务后，将各电气量的振幅、频率、初相位和衰减因子的值发送到管理节点；由管理节点汇总C个计算节点的结果，将结果返回给客户端。

实施例二

(1)广域测量系统对50条联络线的有功功率和频率数据采用Prony分析方法进行参数辨识，相当于100个Prony计算任务。

(2)设在有1个管理节点和2个8核计算节点C₁、C₂组成的机群中，当管理节点接收到客户端的任务N＝100之后，进行Prony并行计算。

(3)管理节点从N个任务中任意选择一个任务，将这个任务同时发送到C₁和C₂上进行计算，记录每个节点对该任务的执行时间。本例中C₁的执行时间为T₁＝0.037s，C₂执行时间为T₂＝0.035s。根据T₁/(T₁+T₂)和T₂/(T₁+T₂)将N₁＝51个任务发送到C₁上计算，N₂＝49个任务发送到C₂上计算。由C₁和C₂分别进行单机多线程并行计算。

(4)C₁接收到的任务数为N₁，C₂接收到的任务数为N₂。指定C₁用于并行的核数为P₁≤8，C₂用于并行的核数为P₂≤8。将N₁个任务依次用P₁个核并行处理，N₂个任务依次用P₂个核并行处理。

表1单个任务不同核数不同采样值计算时间和加速比

(5)C₁将计算得到的51个任务的幅值A、频率f、初相位θ和衰减因子α发送回管理节点，C₂将计算得到的49个任务的幅值A、频率f、初相位θ和衰减因子α发送回管理节点，管理节点将结果转发给客户端。

(6)客户端将C₁和C₂的结果汇总排序，得到50条联络线功率振荡幅值、频率振荡幅值、支路功率初相位、频率振荡初相位、功率振荡幅值和频率振荡幅值。将多个任务的串行执行时间和并行执行时间做对比，结果表明采用本发明方法有效提高了Prony算法的计算速度，结果比较如表2所示：

表2多个任务串行和并行计算时间和加速比

Claims (1)

1.一种基于双重并行计算的在线Prony分析方法，其特征在于，该方法包含下列步骤：

(a)在由I条支路和J个节点组成的交流互联电网中，各支路和所有节点均装设同步测量单元PMU；计算平台由一个客户端、一个管理节点和C个计算节点组成的机群构成，其中C_i表示第i个计算节点，P_i为C_i的核数，C_ij表示C_i中的第j个核，i＝1，2，…，C，j＝1，2，3，…，P_i；

(b)当电网中发生低频振荡时，将PMU采集的M条支路E个电气量包括：有功功率ΔP、频率Δω，取同一时间段采样数据为一组发送给计算平台进行Prony分析，则需要分析的任务数为N＝EM，M＝1，2，…，I；每一类电气量中包含Len个采样数据，500≤Len≤3000；

(c)多机并行计算过程，首先，进行计算节点能力测试，即客户端将接收到的N个任务发送给管理节点，管理节点从N个任务中任意选取一个任务，同时发送给C_i节点上进行计算，i＝1，2，…，C，记录每个节点对该任务的执行时间，节点C_i的执行时间记为T_i，根据第i节点的执行时间与C个节点总执行时间的比例对各计算节点的计算能力进行判别，根据计算能力的强弱对各计算节点进行排序，作为任务分配的依据，然后，将N个计算任务分配到各计算节点开展多机并行计算，节点C_i分配到的任务数为N_i：

$T=\underset{i=1}{\overset{C}{\mathrm{\Σ}}}{T}_i$

$N_i=\left[\frac{T_i}{T}N\right],i=\mathrm{1,2},...,C$

将N_i个任务发送到计算节点C_i上计算，未能整除余下的任务根据计算能力排序依次分配到各计算节点，由C_i进行单节点多线程并行计算；

(d)单节点多核并行计算过程：单节点多线程的并行化采用OpenMP的隐式任务并行和显式任务并行技术实现；当节点C_i接收到的任务数为N_i，每个任务的数据长度为Len，第m个采样数据为X_m，m＝1，2，…，Len；Prony计算过程中取迭代次数k＝Len/2；

此步骤包含以下计算过程：

(d1)最小二乘计算元素求解

$R_{u,v}=\underset{m=k}{\overset{m=N-1}{\mathrm{\Σ}}}{X}_{m-v}{X}_{m-u},u,v=\mathrm{1,2},......,k$

$r_u=-\underset{m=k}{\overset{m=N-1}{\mathrm{\Σ}}}{X}_m{X}_{m-v}$

采用隐式任务并行方式对R_u,v、r_u进行并行计算，采用静态调度方式进行子任务的分配与调度，将R_u,v、r_u平均划分为P_i个子任务，每个子任务负责[k/P_i]×k个数据；由OpenMP将划分得到的P_i个子任务自动分配到核C_ij上执行，i＝1,2,…，C；j＝1，2，3，…，P_i；并行计算后得到R_u,v和r_u，u＝1,2,…，k；v＝1，2,…，k；

(d2)采用高斯约旦消元法求解线性方程组：

$\left[\begin{array}{ccccc}{R}_{\mathrm{1,1}}& R_{\mathrm{1,2}}& R_{\mathrm{1,3}}& ...& R_{1,k}\\ R_{\mathrm{2,1}}& R_{\mathrm{2,2}}& R_{\mathrm{2,3}}& ...& R_{2,k}\\ R_{\mathrm{3,1}}& R_{\mathrm{3,2}}& R_{\mathrm{3,3}}& ...& R_{3,k}\\ ...& ...& ...& ...& ...\\ R_{k,1}& R_{k,2}& R_{k,3}& ...& R_{k,k}\end{array}\right]\·\left[\begin{array}{c}{\mathrm{\α}}_1\\ {\mathrm{\α}}_2\\ {\mathrm{\α}}_3\\ ...\\ {\mathrm{\α}}_k\end{array}\right]=\left[\begin{array}{c}{r}_1\\ r_2\\ r_3\\ ...\\ r_k\end{array}\right]$

得到

α＝(α₁、α₂、...、α_k)^T

采用隐式并行策略实现高斯约旦消元法求解线性方程组的并行化；高斯约旦消元法要求进行k次消元，每次消元过程需要对k个方程进行消元，将k个方程的消元过程分配到P_i个核上并行执行，每个核分配[k/P_i]个方程的消元；划分时，使任务划分的数据间隔为一个高速缓存行的整数倍大小，防止伪共享现象，并在调度上采用动态调度方式。划分后的子任务将由OpenMP自动调度到核C_i1、C_i2、…、C_iPi上计算；计算后得到α₁、α₂、…、α_k；

(d3)将α₁、α₂、…、α_k代入式

X(n)＝-[α₁X(n-1)+α₂X(n-2)+...+α_kX(n-k)]

求解高次方程

z^k+α₁z^k-1+α₂z^k-2+...+α_k＝0

即求解如下矩阵的特征值：

$\left[\begin{array}{ccccc}-{\mathrm{\α}}_1& -{\mathrm{\α}}_2& ...& -{\mathrm{\α}}_{k-1}& -{\mathrm{\α}}_k\\ 1& 0& ...& 0& 0\\ ...& ...& ...& ...& ...\\ 0& 0& 1& 0& 0\\ 0& 0& 0& 1& 0\end{array}\right]$

采用显式任务并行方式实现矩阵特征值的并行求解；求解过程需要k次循环，每次循环过程分为两个部分的计算，一是与其他任务有数据相关的部分，即对矩阵元素更新部分，二是无数据相关的部分，即计算出特征值部分；主线程只负责执行有数据相关的部分，无数据相关的部分利用task指令显式构造并行任务，添加到任务队列中，在OpenMP的调度下由其他线程执行。核Ci1执行主线程，在其对元素更新之后，将其余计算部分显式构造成子任务添加到任务队列中，进而由空闲的核取子任务并执行；最后显式汇总求得：

z₁、z₂、……、z_k。

(d4)求解线性方程组

$\left[\begin{array}{ccccc}1& 1& 1& ...& 1\\ z_1& z_2& z_3& ...& z_k\\ z_1^2& z_2^2& z_3^2& ...& z_k^2\\ ...& ...& ...& ...& ...\\ z_1^{N-1}& z_2^{N-1}& z_3^{N-1}& ...& z_k^{N-1}\end{array}\right]\·\left[\begin{array}{c}{b}_1\\ b_2\\ b_3\\ ...\\ b_k\end{array}\right]=\left[\begin{array}{c}{x}_1\\ x_2\\ x_3\\ ...\\ x_{N-1}\end{array}\right]$

采用隐式任务并行方式完成并行计算，将转置矩阵分为P_i个较小的矩阵块，每个矩阵块大小为(k/P_i)×N，这样就产生P_i个矩阵相乘的子任务，指定调度方式为静态调度，P_i个子任务由OpenMP自动调度到核上完成并行计算；并行计算后得到k个方程，用高斯约旦消元法和隐式任务并行方式进行线性方程组的并行求解，指定调度方式为动态调度；并行计算得到b₁、b₂、…、b_k；

(d5)将以上求出的z_u和b_u代入下式

$\left\{\begin{array}{c}{A}_u=\left|b_u\right|\\ {\mathrm{\θ}}_u=\arctan (\mathrm{lm}\left(b_u\right)/\mathrm{Re}\left(b_u\right))/2\mathrm{\π\Δt}\\ {\mathrm{\α}}_u=\mathrm{lm}\left(z_u\right)/\mathrm{\Δt}\\ f_u=\arctan (\mathrm{lm}\left(z_u\right)/\mathrm{Re}\left(z_u\right))/2\mathrm{\π\Δt}\end{array}\right.,u=\mathrm{1,2},...k$

计算出每个任务的幅值A_u、相位θ_u、频率f_u和衰减因子α_u，u＝1，2，…，k；

对分配到的N_i个计算任务重复步骤(d1)至(d5)，完成节点C_i上的N_i个计算任务，得到N_i个任务对应电气量的振幅、频率、初相位、衰减因子的值；

(e)在每个节点完成其负责的N_i个计算任务后，将各电气量的振幅、频率、初相位和衰减因子的值发送到管理节点；由管理节点汇总C个计算节点的结果，将结果返回给客户端。