一种大规模电网潮流修正方程并行求解方法
【专利摘要】本发明提供一种大规模电网潮流修正方程并行求解方法，根据因子表求解线性程方程运算过程的特点及潮流修正方程系数矩阵的稀疏性，将可以同时进行规格化的行分组，并行完成规格化运算和消去运算，同时，利用求因子表过程中得到的行分组信息，完成并行前代运算，从而提高大规模电网潮流求解效率，满足电网调度员潮流、安全分析、AVC等在线分析应用对潮流计算的实时性要求，便于调度员更快发现电网运行的潜在风险，从而保障电网的安全稳定运行。
【专利说明】一种大规模电网潮流修正方程并行求解方法
【技术领域】
[0001] 本发明涉及一种求解方法，具体涉及一种大规模电网潮流修正方程并行求解方 法。
【背景技术】
[0002] 当前，在国家电网公司的领导下，智能电网调度控制系统在线分析软件实用化工 作已取得了阶段性成果，并在各级调度机构的调度技术支持系统中推广试用，为各调度中 心的日常分析和调度工作起到了非常重要的支撑作用。从在线分析软件计算速度上来看， 需要研究适合大电网实时在线分析软件的并行计算技术和手段，进一步提升在线分析软件 的计算速度，满足大电网全模型实时分析计算的实时性要求。而潮流计算是大电网实时在 线分析软件的重要功能，随着电力系统的发展和电网规模的扩大，在线潮流软件的计算规 模越来越大，同时智能电网精细化调度对大电网在线潮流计算的实时性要求越来越高，使 得传统串行潮流软件不能很好的满足未来智能电网调度需求。为提高大电网在线潮流软件 的计算速度，除了采用高性能服务器及先进算法外，还可以采用分布式计算技术和并行计 算技术来提升软件计算性能。目前各调控中心的调度自动化系统一般使用高性能服务器， CPU个数在4个以上，单个CPU -般为6核及以上，内存一般大于16G，具有良好的并行处理 能力。因此，通过对智能电网调度控制系统现有应用软件进行并行化改造，使得硬件资源得 到充分利用，可以提高系统计算性能。
[0003] 在共享内存模型中，一个并行程序由多个共享内存的并行任务组成，数据的交换 通过隐含地使用共享数据来完成。此编程模式一般仅需指定可以并行执行的循环，而不需 考虑计算与数据如何划分，以及如何进行任务间通信，编译器会自动完成上述功能。目前流 行的共享内存模型开发标准是OpenMP。OpenMP定义了一套编译指导语句，用于指定程序的 并行性、数据的共享/私有等信息，同时OpenMP提供一个支持函数库。目前OpenMP得到了 广泛的支持和使用。微软的VC，linux操作系统上广泛使用的GCC编译器均支持OpenMP编 程模型，此外，IBM、HP、Sun、SGI、Intel等硬件厂商均有支持OpenMP的编译器产品。OpenMP 是一种标准的并行模式fork/join式并行模式，其基本思想是，程序开始时只有一个主线 程，程序中的串行部分都由主线程执行，并行的部分是通过派生其他线程来执行，但是如果 并行部分没有结束时是不会执行串行部分的。利用编译指导语句可将程序分为多个线程并 行执行。
[0004] C++定义了一个内容丰富的抽象数据类型标准库。其中最重要的标准库类型是 string、vector和map，它们分别定义了大小可变的字符串、集合和键-值对的集合。string 和vector往往将迭代器用作配套类型，用于访问string中的字符或者vector中的元素。 map类型为关联数组，可使用键作为下标来获取一个值，正如内置数组类型一样。其关联的 本质在于元素的值与某个特定的键相关联，并非通过元素在数组中的位置来获取。map容器 提供begin和end运算，以生产用于遍历整个容器的迭代器。这些标准库类型是语言组成 部分中更基本的那些数据类型（如数组和指针）的抽象。
[0005] 现有技术中，线性方程组求解主要包括三个步骤：因子分解生成因子表、前代、回 代。线性方程Bx = z的求解过程如下所述：
[0006] (1)因子分解生成因子表
[0007] 因子分解过程包括规格化运算和消去运算两个主要步骤。
[0008] ①规格化运算：对于ηΧη阶矩阵，第i行的规格化只需要对上三角部分中的第i 行非零元进行计算，若第i行j列元素为非零元，则i行对角元bn = Ι/by i行非零对角 兀 bij - bijXbii。
[0009] ②消去运算：利用完成规格化运算的第i行，对i行上三角非零元列号对应的行进 行消去运算，即对第j行进行消去运算，则消去行对角元b"_ = + 。若bu为非 零元，则还需对元素 b#进行消去运算，= 。在规格化过程中会有注入元 素，需要通过节点优化排序进行优化。
[0010] ⑵前代运算
[0011] 将独立矢量Z赋值到工作矢量X中。前代计算从小号节点到大号节点依次进行， 利用因子表上三角参与消去运算，若上三角第i行，第k列元素为非零元，0?（i-Ι)行已完 成消去运算，则x k = xk-bikXi，前代运算的串行程序可以用以下伪代码表示：
[0012] for i=0 Ιο η//循环节点，从0到/卜1 iter-即],begin()//|/J 始化迭代器 iter++/7迭代器 Ijlll while^terHBpj.cndOH+z/判断迭代器是超出范Μ k一iter->first//将键值赋给变.? k x[k]= x[k]-iler->second*x[i】//根据公式；=.υ4 - .Φ-:新 if·算 x* iter++//迭代器 1
[0013] }
[0014] 规格化算法可用以下伪代码表示：
[0015] for i = Oto η// 循环节点，从 0 到 n-1
[0016] iter - B[i]· begin ()// 初始化迭代器
[0017] x[i] = iter->second*x[i]//对非零元Xi进行规格化计算
[0018] (3)回代运算。
[0019] 回代运算从大节点号到小节点号，若i?（η-I)行已完成消去运算，i-Ι行记为k， k行的第j列元素为非零元，则xk = xk_bkjXj，回代运算的串行程序可用以下伪代码表示：
[0020] Γοι· Mm-1 to (V/掀坏 T/点，从/i-l 到(） k-i-\"k=i-\ iter-B[/].begin〇//初始化迭代器 iier++//迭代器加1 whi ilcr! =B [/] ,ctid〇 ).| //判断作超出范 _ itcr->firse//将键值赋治变.Μ:/ x[A，]= x[Ar]-iter->sceond*x|)·]//根公式=? -\.ν; Φ:新?Η?:Λ? iter++//选代器加1 }
[0021] 其中：i，j，k e {〇, 1，2…η}。B为对称系数矩阵，b为对称矩阵非零元，X为解矢 量，z为独立矢量。
【发明内容】
[0022] 为了克服上述现有技术的不足，本发明提供一种大规模电网潮流修正方程并行求 解方法，根据因子表求解线程方程运算过程的特点及潮流修正方程系数矩阵的稀疏性，将 可以同时进行规格化的行分组，并行完成规格化运算和消去运算，同时，利用求因子表过程 中得到的行分组信息，完成并行前代运算，从而提高大规模电网潮流求解效率，满足电网调 度员潮流、安全分析、AVC等在线分析应用对潮流计算的实时性要求，便于调度员更快发现 电网运行的潜在风险，从而保障电网的安全稳定运行。
[0023] 为了实现上述发明目的，本发明采取如下技术方案：
[0024] 本发明提供一种大规模电网潮流修正方程并行求解方法，所述方法包括以下步 骤：
[0025] 步骤1 :稀疏矩阵存储；
[0026] 步骤2 :并行因子分解生成因子表；
[0027] 步骤3 :进行并行前代运算；
[0028] 步骤4:回代运算；
[0029] 步骤5 :潮流迭代运算。
[0030] 所述步骤1中，在P-Q分解法的电网潮流修正方程并行求解过程中完成节点导纳 矩阵、修正方程系数矩阵B'和B"，以及B'和B"因子分解后生成的因子表的存储；修正 方程系数矩阵B'和B"均为稀疏矩阵。
[0031] 采用P-Q分解法的电网潮流修正方程表示为：
[0032] ΔΡ/υ = -Β，ΙΙΔΘ (1)
[0033] AQ/U = _B" AU (2)
[0034] 其中，ΛΡ为计算节点有功修正量矢量，AQ为计算节点无功修正量矢量，Λ θ为 相角修正量矢量，AU为电压修正量矢量，U为节点电压幅值，Θ为节点电压相角，Ρ为节点 注入有功功率，Q为节点注入无功功率，稀疏矩阵元素为对应支路电抗倒数取负，稀疏 矩阵B"元素为对应节点导纳矩阵虚部，i、j为节点编号。
[0035] 所述步骤1具体包括以下步骤：
[0036] 步骤1-1 :利用C++抽象数据类型标准库提供的顺序容器vector和关联容器map 实现P-Q分解法的电网潮流修正方程系数矩阵的存储；
[0037] 定义对象vector〈map〈int, double?Bl，B2, Y分别保存稀疏矩阵B' , B"和节点 导纳矩阵；其中map〈int，double〉用于保存行元素，列号由int型的键索引表示，稀疏矩阵 非零元素由double型的关联值表示；
[0038] 步骤1-2 :节点优化排序；
[0039] 统计节点关联支路数，按节点支路数由少到多对节点进行静态排序；且按照PQ节 点、PV节点、平衡节点的顺序依次排序；在排序的同时，形成内部节点和计算母线编号的对 照表；
[0040] 步骤1-3 :计算节点导纳矩阵；
[0041] 根据节点导纳矩阵维数初始化vector对象Y维数；节点i自导纳即对角元元素等 于该节点直接连接的所有支路导纳的总和，节点i、j之间的互导纳等于连接节点i、j支路 导纳的负值；将计算完的节点导纳阵元素保存到对象Y中，完成节点导纳矩阵的稀疏矩阵 存储，且仅保存节点导纳矩阵上三角；
[0042] 步骤1-4 :稀疏矩阵B'和B "保存到对象B1和B2中。
[0043] 所述步骤2包括以下步骤：
[0044] 步骤2-1 :指定并行处理进程数；
[0045] 步骤2-2 :首先通过稀疏矩阵非零元素的遍历，统计每行元素被消去次数， k行元素被消去次数保存在整型数组变量count [k]中；同时在map〈int, double?B 中增加注入元数位置，将最大被消去次数记为max ;定义变量vector〈int>muster和 vector〈vector〈int>>vec_muster (muster, max)用于收集被消去次数相同的行；vec_ muster[m]中保存被消去次数为m次的所有行号；将被消去0次的行存入vec_master[0] 中；
[0046] 步骤2-3 :根据并行进程数将vec_muster[m]中的行号分成N组，被消去次数m = m+1 ;
[0047] 步骤 2-4 :利用编译指导语句 #pragma omp parallel for num_threads (N)将主 进程分成N个线进程，用多个线程进行并行规格化运算和消去运算，被消去行被消去次数 减1，即count[k]-l，收集被消去次数变成0的行号，并保存在vec_muster[m]中；
[0048] 步骤2-5 :判断VeC_mUster[m]数组是否为空，若为空则并行因子分解结束；否则 返回步骤2-3,继续并行规格化运算和消去运算。
[0049] 所述步骤3包括以下步骤：
[0050] 步骤3-1 :根据潮流有功残差方程和潮流无功残差方程计算可得有功修正量矢量 ΛΡ和无功修正量矢量AQ ;
[0051] 潮流有功残差方程和潮流无功残差方程分别表示为：
[0052]
【权利要求】
1. 一种大规模电网潮流修正方程并行求解方法，其特征在于：所述方法包括以下步 骤： 步骤1 :稀疏矩阵存储； 步骤2 :并行因子分解生成因子表； 步骤3 :进行并行前代运算； 步骤4:回代运算； 步骤5 :潮流迭代运算。
2. 根据权利要求1所述的大规模电网潮流修正方程并行求解方法，其特征在于：所述 步骤1中，在P-Q分解法的电网潮流修正方程并行求解过程中完成节点导纳矩阵、修正方程 系数矩阵B'和B"，以及B'和B"因子分解后生成的因子表的存储；修正方程系数矩阵 和B"均为稀疏矩阵。
3. 根据权利要求2所述的大规模电网潮流修正方程并行求解方法，其特征在于：采用 P-Q分解法的电网潮流修正方程表示为： ΔΡ/U = -B/ UA θ (1) AQ/U = _B" AU (2) 其中，ΛΡ为计算节点有功修正量矢量，AQ为计算节点无功修正量矢量，Λ θ为相角 修正量矢量，AU为电压修正量矢量，U为节点电压幅值，Θ为节点电压相角，Ρ为节点注入 有功功率，Q为节点注入无功功率，稀疏矩阵Β'元素为对应支路电抗倒数取负，稀疏矩阵 Β"元素为对应节点导纳矩阵虚部，i、j为节点编号。
4. 根据权利要求2所述的大规模电网潮流修正方程并行求解方法，其特征在于：所述 步骤1具体包括以下步骤： 步骤1-1 :利用C++抽象数据类型标准库提供的顺序容器vector和关联容器map实现 P-Q分解法的电网潮流修正方程系数矩阵的存储； 定义对象vector〈map〈int, double?Bl，B2, Y分别保存稀疏矩阵Β'，B"和节点导纳 矩阵；其中map〈int，double〉用于保存行元素，列号由int型的键索引表示，稀疏矩阵非零 元素由double型的关联值表示； 步骤1-2 :节点优化排序； 统计节点关联支路数，按节点支路数由少到多对节点进行静态排序；且按照PQ节点、 PV节点、平衡节点的顺序依次排序；在排序的同时，形成内部节点和计算母线编号的对照 表； 步骤1-3 :计算节点导纳矩阵； 根据节点导纳矩阵维数初始化vector对象Y维数；节点i自导纳即对角元元素等于该 节点直接连接的所有支路导纳的总和，节点i、j之间的互导纳等于连接节点i、j支路导纳 的负值；将计算完的节点导纳阵元素保存到对象Y中，完成节点导纳矩阵的稀疏矩阵存储， 且仅保存节点导纳矩阵上三角； 步骤1-4 :稀疏矩阵K和B"保存到对象B1和B2中。
5. 根据权利要求1所述的大规模电网潮流修正方程并行求解方法，其特征在于：所述 步骤2包括以下步骤： 步骤2-1 :指定并行处理进程数； 步骤2-2 :首先通过稀疏矩阵非零元素的遍历，统计每行元素被消去次数，k行 元素被消去次数保存在整型数组变量count [k]中；同时在map〈int, double?B中 增加注入元数位置，将最大被消去次数记为max ;定义变量vector〈int>muster和 vector〈vector〈int>>vec_muster (muster, max)用于收集被消去次数相同的行；vec_ muster[m]中保存被消去次数为m次的所有行号；将被消去0次的行存入vec_master[0] 中； 步骤2-3 :根据并行进程数将vec_muster[m]中的行号分成N组，被消去次数m = m+l ; 步骤2-4:利用编译指导语句#pragma omp parallel for num_threads(N)将主进程 分成N个线进程，用多个线程进行并行规格化运算和消去运算，被消去行被消去次数减1， 即count[k]-l，收集被消去次数变成0的行号，并保存在vec_muster[m]中； 步骤2-5 :判断vec_muster[m]数组是否为空，若为空则并行因子分解结束；否则返回 步骤2-3,继续并行规格化运算和消去运算。
6.根据权利要求1所述的大规模电网潮流修正方程并行求解方法，其特征在于：所述 步骤3包括以下步骤： 步骤3-1 :根据潮流有功残差方程和潮流无功残差方程计算可得有功修正量矢量Λ P 和无功修正量矢量Λ Q ; 潮流有功残差方程和潮流无功残差方程分别表示为：
(3) (4) 其中，Λ Pi为节点i有功修正量，Λ Qi为节点i无功修正量；Λ Θ为相角修正量矢量， Λ U为电压修正量矢量；％和％分别为节点i和j电压幅值；Θ u为节点i和j的相位差即 θ υ = Θ Θ yPi为节点i注入有功功率，Qi为节点i注入无功功率；为节点导纳矩阵 实部，节点导纳矩阵虚部； 步骤3-2 :将独立矢量z赋值到工作矢量X中，指定进程数N ; 独立矢量z表示方程组（1)、（2)的常数项Λ P/U和Λ Q/U，工作矢量X表示方程组（1)、 ⑵的待求项Λ Θ和AU ; 步骤3-3 :初始化被消去次数m = 0 ; 步骤3-4 :取数组vecjnuster [m]中的行，若行数大于进程数N则分成N组，否则不分 组，进行并行前代消去运算，被消去次数m = m+Ι ; 步骤3-5:前代计算从小号节点到大号节点依次进行，利用因子表上三角参与并 行前代消去运算，若上三角第i行，第k列元素为非零元素，0?（i-1)行已完成并行 前代消去运算，则xk = xk-bikXi ;前代的消去过程和生成因子表的消去过程一致，利用 vector〈vector〈int?vec_muster中的行分组信息进行并行前代消去运算； 步骤3-6 :被消去次数m = m+Ι，判断被消去次数是否等于最大消去次数，若相等则前代 消去运算结束，若小于则返回步骤3-4继续进行并行前代消去运算； 步骤3-7 :前代规格化运算； 工作矢量X的第i个元素Xi = Xi/bii，前代规格化运算用编译指导语句#pragma omp parallel for完成for循环的前代规格化运算，即根据CPU最大并行处理线程数前代规格 化运算。
7. 根据权利要求1所述的大规模电网潮流修正方程并行求解方法，其特征在于：所述 步骤4中，回代运算从大节点号到小节点号，若i?（n-1)行已完成回代消去运算，i-Ι行 记为k，k行的第j列元素为非零元素，则x k = xk-bkjXj ;按式xk = xk-bkjXj，k由大到小逐个 对工作矢量的元素进行求解，最终的到方程组（1)、（2)的解。
8. 根据权利要求1所述的大规模电网潮流修正方程并行求解方法，其特征在于：所述 步骤5包括以下步骤： 步骤5-1 :根据得到的相角修正量矢量Λ Θ和电压修正量矢量AU，对节点电压幅值U 和节点电压相角Θ进行修正，即θ = θ+ΛΘ，υ = υ+Λυ; 步骤5-2 :按式（3)、（4)计算各节点的有功修正量和无功修正量； 步骤5-3 :若有功修正量小于有功收敛精度，且无功修正量小于无功收敛精度，则计算 完成；否则，返回步骤3和4重新进行并行前代运算和回代运算，并求解式（1)、（2)。
【文档编号】H02J3/00GK104158182SQ201410405872
【公开日】2014年11月19日 申请日期:2014年8月18日 优先权日:2014年8月18日
【发明者】李静, 罗雅迪, 徐杰, 贾育培, 郎燕生, 郭子明, 邹昱, 施贵荣, 窦成龙, 张 浩, 陈利杰, 王磊, 王淼, 王少芳 申请人:国家电网公司, 中国电力科学研究院, 国网冀北电力有限公司