基于OpenMP并行简约空间内点法的暂态稳定紧急控制

王 云，江全元

（浙江大学 电气工程学院，浙江 杭州 310027）

0 引言

暂态稳定紧急控制是维持电力系统安全稳定运行的重要手段。暂态稳定紧急控制问题通常建模为考虑暂态稳定约束和系统经济运行的最优化问题。

目前，求解暂态稳定紧急控制问题的算法主要有试凑法或启发式算法[1-3]、基于暂态能量函数的直接法和基于最优控制理论的优化算法。暂态稳定紧急控制问题的难点之一在于暂态稳定约束的描述，基于能量函数的方法在暂态稳定预防控制和暂态稳定紧急控制方面已有广泛的研究[4-11]。时域仿真法是解决暂态稳定紧急控制的另一种方法，基于最优控制理论，各种紧急决策计算均描述为以控制代价最小为目标、以维持系统安全稳定性为约束条件的最优控制问题，这类最优控制问题实质上就是最优参数选取问题。文献[12-13]将切负荷问题描述为最优控制问题，提出了一种系统化地决定最小切负荷量及其在各切负荷点分配的算法。文献[14]基于时域仿真得到的系统受扰轨迹给出暂态稳定紧急控制的非线性模型，然后采用近似规划法和拟贪婪法求解。

简约空间技术主要用来求解自由度较小的非线性规划问题，已在化工领域获得了广泛应用[15]。文献[16-17]用简约空间内点算法求解暂态稳定最优潮流问题，计算效率获得极大的提升。文献[18]应用简约空间二次规划法求解暂态电压稳定紧急控制问题，显著提高了计算效率。由于暂态稳定紧急控制问题中可切负荷/发电机节点相对于问题规模非常小，即该类问题自由度很低，本文运用简约空间内点算法进行求解。此外，算法采用C++编程实现，对于算法的关键耗时环节，充分利用多线程并行技术。多个算例测试结果表明，本文提出的基于OpenMP并行简约空间内点法在求解紧急控制问题中，计算时间明显缩短，占用内存减小，能够求解大规模系统的紧急控制问题。

1 暂态稳定紧急控制数学模型

电力系统最优切机切负荷控制问题是一个大规模动态优化问题，其一般形式为：

其中，x为优化变量；F为优化问题目标函数；H为含微分方程等式约束；G为含微分方程不等式约束；、分别为不等式约束上、下限。

采用隐式梯形积分方法对上述动态优化问题差分化，动态优化问题（1）可转化为如下非线性规划问题：

1.1 目标函数

暂态稳定紧急控制的目的是切除最少的发电机和负荷来保持系统暂态稳定，因此，本文将优化目标函数设为切机和切负荷量加权最小：

其中，u=[uG，uL]，为优化控制变量，uG、uL分别为各可切发电机节点切机比例组成的向量和可切负荷节点切负荷比例组成的向量；cG、cL分别为切机、切负荷控制代价的权值，用来表征可切发电机和负荷的重要性，本文中将权值均取为1，即所有可切发电机、可切负荷重要性相同；PG、PL分别为各可切发电机节点装机容量和各可切负荷节点总负荷。

1.2 等式约束

a.发电机动态方程。

电力系统经典模型对于系统第一摆稳定分析非常有效实用[19]。本文中发电机模型采用经典模型，负荷采用恒阻抗模型，因此，各发电机可由如下两阶微分方程描述：

其中，i=1，2，…，NG（NG为发电机节点数）；为发电机暂态电势，经典模型下为恒定值，E′i、x′di分别为发电机节点 i暂态电势、直轴暂态电抗，分别为发电机节点i在t时刻的功角、节点电压、节点电压实部、节点电压虚部；TJi、Pmi分别为发电机节点i惯性时间常数和有功出力；ωti为发电机节点i在t时刻的角速度；ωs为同步角速度有名值。

b.电网方程。

将发电机表示成电流源形式，并将发电机等值导纳和负荷等值导纳并入网络，电网方程可写成：

其中，Ixti、Ityi分别为t时刻发电机节点i的注入电流实部和虚部；Gij、Bij分别为节点i和节点j节点导纳实部和虚部；nB为系统节点数。

对于可切机的发电机节点，有：

对于不可切机的发电机节点，有：

对于可切负荷节点，有：

对于不可切负荷节点，有：

其中，G′ii、B′ii分别为发电机等值导纳和负荷等值导纳并入网络前节点自导纳实部和虚部；GLi、BLi分别为负荷等值导纳的实部和虚部；uGi、uLi分别为切机点切机比例和切负荷点切负荷比例。

1.3 不等式约束

a.控制变量上下界约束：b.暂态稳定约束。

暂态稳定判据采用中性惯量形式表示：

2 基于OpenMP并行简约空间内点法

2.1 简约空间内点算法

采用预测－校正内点法求解非线性规划问题（2），其一阶最优必需条件（KKT条件）求解包括预测和校正2步，其主要计算量在求解以下对称线性方程组[20]：

简约空间技术适用于自由度较低的非线性规划问题，用来求解对称线性方程组（12）。将变量x分成程空间DY和零空间DZ两部分，DY由m个非独立变量组成，DZ由n-m个独立变量组成。变量x的解可写成：

其中，pYєDY；pZєDZ；YєRn×m、ZєRn×（n-m）分别称为程空间和零空间基矩阵，[Y Z]非奇异。零空间基矩阵Z满足：

则 pY、pZ、Δy 分别如下求得：

详细的简约空间内点算法原理及实施细则可参考文献[16]。采用简约空间内点算法求解紧急控制问题时，方程组（15）—（17）的求解通过调用 KLU[21]计算，其对矩阵C的分解可复用，每次迭代只需做1次LU分解。算例分析表明，简约空间下紧急控制算法的主要计算时间集中在矩阵B和Zm的计算上，本文通过多线程技术提高这部分的计算效率。

2.2 OpenMP并行算法

2.2.1 OpenMP编程模型

为充分利用多核处理器计算平台的计算性能，针对算法可解耦部分利用多线程技术实现并行计算从而加快计算速度。OpenMP[22]提供了一套共享存储体系结构下的多线程编程模型，从而实现本文算法多线程并行技术。

如图1所示，OpenMP采用Fork-Join的并行执行模型：程序首先以单线程（主线程）开始，类似一个串行程序，当遇到并行结构，主线程会产生一组线程（Fork动作），每个线程都执行并行动态扩展域中的代码；并行结构执行完后，只有主线程继续执行，其他所有的线程结束执行（Join动作）。本文程序并行化部分，系统给每个线程分配一个CPU核心。

图1 Fork-Join模型Fig.1 Fork-Join model

2.2.2 并行简约空间内点算法

在简约空间内点算法框架下，多线程技术可通过下面2种方式实现。

a.矩阵Zm的多线程求解。

矩阵N是由n-m个m维列向量构成，在回代过程，矩阵N各个列向量回代过程相互独立。将N分解成p（p为线程数）个m维列向量组成的矩阵，即N=[N1N2… Np]。可将这p个矩阵分配给各计算线程实现并行回代，从而加快计算速度，如图2所示。KLU解法器是一个串行解法器，不支持并行计算。本文对KLU代码进行了修改，使其按图2所示支持多线程回代。

图2 Zm多线程计算Fig.2 Multi-thread calculation of Zm

b.矩阵B的多线程计算。

本文算法采用C++编程实现，算法实施过程中涉及到较多矩阵操作，尤其是B的计算占据较多时间。目前有许多高性能线程代数库可供调用，很多都支持多处理器平台下的多线程并行。在本文中，基础线程代数库（BLAS）被用来作一些矩阵与矩阵相乘运算和矩阵与向量相乘运算，在计算矩阵B时主要调用DCSRMM和DGEMM函数进行求解，这2个函数均支持OpenMP下多线程计算，程序实现时通过调用MKL_NUM_THREADS函数设置线程数，开启多线程从而加快计算速度。此外，线性代数解法器（LAPACK）被用来求解稠密线性方程组（18）。

2.3 算法实施流程

采用本文提出的算法求解暂态稳定紧急控制问题的主要步骤如下：

a.读入算例数据，包括故障信息；

b.对故障后系统进行暂态稳定仿真计算，若系统保持稳定，不需要切机切负荷，跳到步骤g；

c.按照第1节内容对暂态稳定紧急控制问题数学建模；

d.对系统变量及收敛条件初始化，k=0，算法优化开始；

制定智能化建筑用电节能设计时，一定要整体把握用电设施的布局分配、负荷容量及功率大小等信息，以此为基础，选择恰当、合理的节能供配电设备，保证用电设备正常工作的同时，最大限度地降低能耗的损失。

e.按照第2节提出的并行简约空间内点算法计算变量x更新步长；

f.判断是否达到收敛条件，若未收敛，k=k+1，跳到步骤e；

g.算法退出，输出结果。

3 算例分析

3.1 测试算例

本文采用基于OpenMP并行简约空间预测校正内点法对4个算例进行了计算，这4个算例是在IEEE标准算例的基础上修改得到，根据切机量和切负荷量相对功角稳定的灵敏度大小来确定可切发电机和负荷。表1给出4个测试算例的参数，其中nB、nL、NG、nG、Nload、nload分别表示系统节点数、系统线路数、系统发电机数、系统可切发电机数、系统负荷数、系统可切负荷数。对所有测试系统，在t=0时刻发生三相短路故障，t=0.2 s切除故障线路，t=0.3 s执行切机切负荷操作。tsim为算法仿真时间，即为切负荷时刻至末端时刻的时间差。采用隐式梯形积分法进行暂态稳定分析时，积分步长在本文中取为Δt。程序中各发电机相对于惯性中心的偏差不超过±110°作为暂态稳定判据。

表1 测试算例参数Table 1 Parameters for case tests

表2给出了4个测试算例采用隐式梯形积分差分化方法离散微分方程和网络方程后问题规模，其中n表示问题状态变量数，m表示问题等式约束数，δDOF表示问题自由度即为n与m之差，rDIM、NNZ分别表示原对偶系统式（12）的维数和非零元数量。可以看出，差分化后优化问题自由度很小，即满足，从而该问题非常适合简约空间下求解。

表2 测试算例问题规模Table 2 Scales of case tests

本文算法采用C++编程实现，测试环境为Intel Xeon 2 quad-core CPU，OpenMP被用来实现多线程并行。通过调用英特尔提供的高性能多线程线性代数库（BLAS和LAPACK）实现算法中的一些矩阵运算。KLU解法器被用来对算法中稀疏线性方程组进行求解。

3.2 测试结果

本文算法是基于OpenMP并行简约空间内点法。为表述简洁，本文称传统内点算法为全空间内点算法。为了比较本文算法与常规算法的性能优劣，对上面4个算例用隐式梯形差分化方法离散微分方程和网络方程，并在全空间、简约空间以及并行简约空间下求解，测试结果如表3所示。表3中“×”表示因时间或者内存不足导致问题不可解。

表3 测试结果Table 3 Results of case tests

从表3中可以看出，采用隐式梯形积分差分微分方程和网络方程，简约空间下计算速度明显比全空间下计算速度要快，并且可以求解更大规模系统紧急控制问题。因简约空间技术只是在求取式（12）时采取的策略，其并不改变迭代过程，简约空间和全空间下问题迭代进程是完全一样的，得出的优化结果也一样。因此，本文将把全空间和简约空间得出问题的迭代次数统一列出。

基于OpenMP并行简约空间算法是在简约空间算法的框架下对一些计算步骤利用多线程技术并行求解，它不改变算法的计算进程，因而两者有相同的迭代进程和相同的优化结果。从表3可以看出，在8核处理器计算平台下，并行简约空间内点算法计算效率进一步提高。

本文选取切机切负荷作为紧急控制措施来保证暂态稳定性。以CASE300为例，列出优化后切机、切负荷量如表4所示。

表4 切机切负荷优化结果Table 4 Results of generator trip and load shedding optimization

3.3 并行性能分析

首先对串行简约空间下算法各部分计算时间进行分析如表5所示。从表5可以看出，简约空间内点算法下暂态稳定紧急控制问题的计算时间主要集中在B计算时间和Zm回代时间。因而本文对这部分进行多线程并行求解，提高计算效率。

表5 串行简约空间内点算法计算效率分析Table 5 Calculation efficiency analysis of serial reduced-space interior point method s

图3 Zm回代过程加速比Fig.3 Acceleration ratio of Zmback substitution

图4 B计算过程加速比Fig.4 Acceleration ratio of B calculation

为了分析B计算和Zm回代的并行效率，线程数p取1到8分别进行测试，测试结果如图3、图4所示。从图3、图4可以看出，当线程数较少时，本文算法几乎可以达到线性加速比，随着线程数的增加，加速比增长放缓，在8核处理器下，加速比最大可达4.5。一般而言，测试算例规模越大，本文算法可获得越好的加速比。因此，本文算法对大规模系统暂态稳定紧急控制问题的求解有较好的潜力。

3.4 时域仿真验证

通过时域仿真来验证上面4个算例的第一摆稳定性，图5给出了各个测试算例发电机摇摆曲线，图中实线为系统发生故障后不执行紧急控制操作的功角曲线；虚线为执行切机切负荷操作的功角曲线。因全空间内点算法、简约空间内点算法和并行简约空间内点算法三者有相同的迭代进程，因此3种算法得出的优化曲线一致，本文只给出简约空间下发电机摇摆曲线。由图5可以看出，电力系统发生故障后不执行切机切负荷操作将失去稳定。采用本文提出的算法执行切机切负荷操作，系统能够保持暂态稳定。

图5 测试算例时域仿真结果Fig.5 Time-domain simulative results of case tests

4 结语

电力系统暂态稳定紧急控制问题是一类大规模非线性动态优化问题。本文提出了一种基于OpenMP的并行简约空间内点法求解该类问题。测试结果表明，相对传统内点算法，本文提出的并行简约空间内点算法计算效率明显提高，能计算更大规模系统，在多核处理器平台下，多线程技术进一步提高了算法计算速度。随着多核处理器的普及，本文算法实际效益更为明显。