暂态稳定约束

暂态稳定约束（共7篇）

暂态稳定约束 篇1

摘要：提出了一种暂态稳定的发电机有功调度新算法。通过求取暂态稳定约束下的发电机有功输出极限,将暂态稳定约束转化为最优潮流中的发电机有功出力约束,从而将暂态稳定约束下的发电机有功调度问题转化为常规最优潮流问题进行求解。介绍了电力系统角半径的概念,通过仿真计算轨迹灵敏度,在此基础上可得到系统角半径对各发电机有功输出的灵敏度,该灵敏度信息可用于计算发电机有功输出的调整量。算法实现简单,可处理多个预想故障。在10机新英格兰典型电力系统上的算例分析验证了新方法的有效性和合理性。

关键词：电力系统,暂态稳定,最优潮流,轨迹灵敏度

0 引言

考虑各种约束的发电调度问题一直是研究人员关心的重点,最优潮流方法为该问题提供了很好的解决手段。电力系统最优潮流(Optimal Power Flow OPF)[1,2]是在满足各种运行约束条件下,通过调整系统中可利用的控制手段实现预定目标最优的非线性规划问题。学术界对OPF问题进行了大量的研究。随着电网规模的扩大和负荷需求的迅猛增长,电力系统面临的暂态稳定性问题变得日益突出。近年来,考虑暂态稳定约束的最优潮流(OPF with Transient Stability constraints,OTS)[3,4,5,6,7,8,9]引起了各国学者的广泛关注。

OTS的暂态稳定约束条件中包括一系列的微分方程,难以直接求解。文献[3-5]将暂态稳定约束中的微分方程差分化为一系列代数方程,使OTS问题可用常规的OPF方法求解,具有模型精确、结果可靠等优点,但会导致优化问题的变量和方程数显著增加,可能存在计算量大或收敛性不好等问题。文献[6]将势能边界面法与最小二乘法相结合,建立系统暂态稳定裕度与系统运行参数变化之间非线性关系的函数表达式,从而将暂态稳定约束条件引入到常规最优潮流模型中。文献[7]将暂态稳定裕度约束加入最优潮流模型,以此得到OTS的解。文献[8]提出了一种考虑联络线暂态稳定约束的最优潮流计算方法。文献[9]把暂态稳定约束转化为控制变量的不等式约束,与OPF问题相比只增加了少量约束条件,有效地简化了计算,是一种较实用的方法。

本文提出一种考虑暂态稳定约束的发电机有功再调度新方法,在迭代过程中OPF运行点上求取相关机组在暂态稳定约束下的有功输出极限,把微分方程约束条件转化为发电机有功出力的不等式约束,然后回到OPF主程序进行求解。当有危害故障都被排除后,即得到最终解。本文首先介绍了电力系统角半径的概念,系统角半径对发电机有功出力的灵敏度通过轨迹灵敏度计算得到,该灵敏度用于计算暂态稳定约束下的发电机有功输出极限。上述方法将暂态稳定约束与OPF有机地结合在一起,只需在OPF的基础上增加求取发电机有功输出极限的模块,易于用传统的OPF方法实现。最后在10机新英格兰典型电力系统上验证了本文算法的有效性和合理性。

1 轨迹灵敏度分析

1.1 电力系统角半径

对包含n台发电机的电力系统,定义δi和θi=δi-δCOI为第i台发电机转子相对于系统同步转速和系统惯性中心(COI)的角度,其中是系统惯性中心的角度,Mi为第i台发电机的惯性时间常数,发电机相对系统惯性中心的功角构成一个以惯性中心为坐标原点的n维角度空间。电力系统动态中,系统状态轨迹可映射为各发电机功角(θ1,θ2,,θn)对应的相点在该n维空间中的运动轨迹。令

式中,R称为电力系统角半径,在数学上,系统角半径是某时刻所有发电机相对惯性中心的功角θi的欧几里得范数,相当于n维空间中相点与原点(惯性中心)的距离,可作为发电机功角相对惯性中心摆开程度的一种度量。

R随时间的变化反映了电力系统在遭受大扰动后能否保持功角稳定运行[10]。若系统中一台或多台发电机失去同步,则至少有一台发电机的功角满足θi→∞|t→∞,由式(1)可知此时R→∞|t→∞;反之,若系统是功角稳定的,R将在某一固定范围内变化。上述特性在对不同规模的电力系统的仿真分析中得到了验证[10]。

1.2 轨迹灵敏度数学模型[11]

轨迹灵敏度分析是针对微分方程模型研究状态变量变化对参数依赖性的方法。

多机电力系统的动态行为可用如下的微分代数方程组描述

式中:x和y分别表示系统状态变量和代数变量;x0和y0为其对应的初值;α是可控参数,本文中为发电机有功出力。

式(2)的解记作x(α,t)和y(α,t),对其按泰勒级数展开并忽略高阶项可得

式中:xα(t)和yα(t)分别为x和y对α的灵敏度,称作轨迹灵敏度,反映了系统轨迹与可控参数增量的关系。

式(2)对α求偏导数可得

xα(t)的初值[12]为

假设gy非奇异,即gy可逆,由式(5)可得yα(t)的初值如下

式(5)~式(7)构成电力系统轨迹灵敏度数学模型,是系统式(2)的伴随系统。若系统仿真采用隐式梯形积分法,则只需在每一仿真时步增加一次线性方程求解的工作即可得到轨迹灵敏度系统式(5)~式(7)的解。

令式(1)对α求偏导数可得到系统角半径R对可控参数的灵敏度。

其中

式中可通过求解轨迹灵敏度系统式(5)~式(7)获得。

2 OPF模型

本文研究的最优潮流问题以发电成本最小为目标函数,写作

式中:Pgi是第i台发电机的有功出力,ai、bi和ci为第i台发电机的费用系数,SG是可调发电机的集合。

约束条件包括

式(11)~式(15)分别是功率平衡方程、电源有功出力上下限约束、无功电源出力上下限约束、节点电压上下限约束和线路热稳定约束。式中:Pg和Qr分别是母线有功和无功注入功率向量;PL和QL分别为母线有功和无功负荷向量;V和ϕ分别是母线电压的幅值和相角;Sij是流过线路(i,j)的视在功率;SR是可调无功电源集合;SB是系统所有节点集合;Sl为所有支路集合。

3 算法步骤

3.1 暂态稳定约束下的有功输出极限

式(10)~式(15)是OPF的数学模型,应在此基础上增加暂态稳定约束条件,通常采用式(16)所示的启发式判据作为暂态稳定约束。

式中:θmax是研究时间内发电机相对于惯性中心的角度的最大值;ρ为指定的功角门槛值。

直接对含暂态稳定约束的有功再调度问题进行求解是困难的,本文首先计算发电机有功输出极限,再用其增补约束条件式(12),从而将暂态稳定约束转化为发电机有功出力约束,最后求解常规OPF问题最终解。

对有危害故障,计算相关发电机组的有功调整量,由于电力系统是强非线性系统,因此还需要与迭代方法相结合,最终得到暂态稳定约束下的发电机组有功输出极限。对有危害故障,令

式中,tρ是θmax超过ρ的时刻,为节省计算时间,tf的取值由以下方法[13]确定:令ζ=(2.0~3.0)ρ,当θmax≥ζ时,结束仿真,以该时刻作为tf;若在仿真时间内始终有ρ<θmax<ζ,则令tf为发电机功角最大值对应的时刻。

本文采用C≤0作为最优调度的目标,即将系统角半径的变化控制在一定的范围内。上文推导了系统角半径R对α的灵敏度,由C的定义可知C对α的灵敏度S=∂C∂α直接可求。

以发电机有功出力P作为控制参数α,则控制目标可近似写为如下形式

式中:是C(P)对各发电机有功出力的灵敏度向量,∆P是各发电机有功出力调整量的向量。

调整机组出力后,系统发电成本增量为

式中,∆Pgi是第i台发电机的有功调整量。

将系统发电成本增量最小作为目标函数,其数学模型可表示为

式中发电机有功输出向量P为常数,待优化变量是发电机有功输出调整量向量∆P。本文的目的在于求取发电机在暂态稳定约束下的有功输出极限(上限),因此计算时可作以下近似处理:只考虑Si为正,即应减少有功输出的发电机组,令SG+表示Si为正的可调发电机的集合,则式(20)可简化为以下形式:

求解过程中,有功出力调整导致的有功不平衡量由Si为负的可调发电机组补偿。上述策略可以显著降低最优化问题的维数,提高计算效率。式(21)用于在迭代过程中计算有危害故障的发电机有功输出调整量,可以用一般的优化方法求解。由于电力系统是强非线性系统,上述方法需要与迭代计算相结合,以获得精确解。

有功输出极限的计算流程如图1所示。每一迭代时步中调整有功输出后可能出现两种情况:(1)系统仍不稳定,此时重新进行轨迹灵敏度计算并求解最优化问题式(21),得到∆Pk;(2)系统稳定,判断收敛条件是否满足(ε是收敛允许误差),若收敛,发电机有功输出极限Pcr=Pk,计算结束;否则调整修正量,继续进行下一步迭代计算,其中系数β∈(0,1)用于修正∆Pk,防止过度调整。

3.2 算法流程

电力系统的可能发生的故障很多,但大多数故障并不对系统的稳定运行构成威胁。因此,首先通过故障扫描滤除无危害故障,并得到有危害故障的临界切除时间(CCT)或其他稳定指标,将有危害故障按其失稳机组是否相同进行分组,只对每组中的最严重故障进行分析。

以常规OPF的解作为初值,假设通过故障扫描已得到有危害故障及其稳定指标,求解步骤如下:

1)OPF计算。

2)采用3.1节介绍的方法计算最严重故障对应的发电机有功输出极限,在约束条件式(12)中增加相应发电机的输出极限约束。

3)根据修改后的约束条件重新进行OPF计算,分析有危害故障集中的其他故障,若无暂态稳定问题,转步骤4),否则转步骤2)。

4)结束。

4 算例分析

4.1 计算结果

采用如图2所示的10机新英格兰典型电力系统验证本文算法的有效性,系统基本参数见文献[14],其中1号发电机作为平衡机,系统中所有发电机采用经典模型,负荷采用恒阻抗模型,仿真时间为4 s,暂态稳定功角门槛值ρ=180°,ζ=540°,发电机有功输出极限收敛精度为ε=0.1 MW,系数β=0.5。各发电机的有功输出上下限、费用系数和初始OPF解如表1所示。最优潮流计算程序采用MATPOWER 3.0软件包[15]。由表1数据计算得到未考虑暂态稳定约束的系统最小发电成本为61 745.4$($为某一货币单位)。

令故障切除时间tcl=0.12 s,对所有线路首末端三相短路故障进行扫描,滤出有危害故障集合,并计算其CCT所在区间。表2中第一列是故障线路,“29-28”表示29号母线在0时刻发生三相短路故障,故障后0.12 s切除该线路。第二列是调整前的CCT,第三列是考虑暂态稳定后重新计算的CCT,可见有危害故障都已被排除。

表3中第二列是初始状态下,对有危害故障“29-28”计算得到的灵敏度S,只有对9号发电机的灵敏度为正值,即只需计算9号发电机的有功输出极限,原问题简化为一维最优化问题,大大减轻了计算负担。通过计算,得到在暂态稳定约束下9号发电机的有功输出极限为739.42 MW,计算结果如表3所示。

4.2 算法效率分析

本文算法通过计算暂态稳定约束下的发电机有功输出极限,将暂态稳定约束转化为发电机有功出力约束,该方法不增加优化问题的维数,复杂度与常规OPF相同,适用于大规模电力系统的分析计算。

对失稳发电机组相同的有危害故障,只选取其最严重的进行处理后,通常情况下其他故障也转为稳定。例如在本文算例中,有危害故障的失稳发电机相同,均为9号发电机,线路“29-28”首端三相短路故障为最严重故障,对其计算有功输出极限并调整有功出力后,其他故障也转为稳定。该策略可提高分析效率。

系统角半径R包含了系统中所有发电机的功角信息,因此该方法适用于多机失稳的情况。在计算有功输出极限的迭代过程中,只考虑灵敏度为正的发电机,以及动态选择tf等措施,都有效地简化了计算。

算例中对故障“29-28”求取相关机组的有功输出极限进行了8次迭代,用时约16 s;调整后其他故障已无稳定问题,对这些故障进行仿真校验的时间为12 s左右;共使用MATPOWER进行了两次OPF计算,用时约8 s。因此,本文算例中计算总用时约为36 s。

5 结论

本文提出了一种考虑暂态稳定约束的发电再调度新方法。该方法利用轨迹灵敏度信息,结合迭代计算的方法,求解暂态稳定约束下的发电机有功输出极限,将暂态稳定约束转化为发电机有功出力约束。该方法不受元件模型限制,可处理多个预想故障,对新英格兰典型电力系统的分析结果验证了本文算法的有效性。

暂态稳定约束 篇2

随着我国电网规模的进一步扩大以及大容量发电机的不断投入,重要的500 k V枢纽变电站高压侧三相短路电流已经面临着超过断路器遮断电流的威胁,严重降低了系统的抗风险能力和调度的灵活性。采用串联电抗器限制500 k V母线三相短路电流是一种有效措施[1,2,3],但是串联电抗器价格较为昂贵,而且会增长电气距离,降低受端母线电压,并对电网稳定性带来不利影响,所以在系统中不宜加装过多。因此,对串抗器进行优化配置、减少加装串抗器的数量和阻抗显得尤为重要。

传统配置方法主要依靠经验和反复的试验,效率低,而且无法把握全局效果。文献[4]提出一种电力系统限流措施的优化方法,将限制短路电流问题转化为一个常规混合整数规划问题进行求解。文献[5]提出一种基于粒子群算法的串联电抗器优化方法,在测试系统及实际电网中取得良好效果。但是目前所有串联电抗器优化算法并没有在优化过程中计及暂态稳定约束,而是通过对优化结果的稳定性校验判断配置方案是否满足稳定性的要求,若不满足稳定性要求,则需要修改串联电抗器容量范围,重新进行优化。这一过程往往需要依赖实际工程经验,而且增加了工作量。所以这种不考虑稳定约束的串联电抗器优化配置算法并不全面,在优化过程中计及暂态稳定性是十分必要的。

在优化中处理暂态稳定约束的方法有2类:时域仿真法和暂态能量函数法。基于时域仿真的方法有差分化方法和基于约束转换技术的方法[6]。前者将系统动态方程差分为代数方程从而建立静态优化模型,采用优化方法求解。后者利用约束转换技术处理包含代数微分方程组的附加约束,将函数空间的优化问题转换为常规的静态优化问题求解。在实际应用时,差分化方法可能会出现维数灾问题,而约束转化技术虽然降低了系统的规模,但是每次计算负担过重。而暂态能量函数法主要是利用能量函数法计算系统稳定裕度及其灵敏度,进而判断系统是否满足约束条件,常用的方法有单机等值法[7]、故障模式法[8]、基于稳定域边界的主导不稳定平衡点法(BCU法)[9]等。相比而言,能量函数法易于求得系统稳定裕度,并能快速计算线路故障临界切除时间(CCT)。

本文在文献[5]的基础上利用BCU法计算暂态稳定裕度,提出考虑暂态稳定约束的串联电抗器优化算法。根据串联电抗器对系统稳定性影响的特点,确定暂态稳定约束所需考虑的线路故障范围,减少暂态稳定分析计算量。在优化过程中对这些线路故障的稳定裕度进行计算,从而保证优化结果满足稳定性的要求。通过测试系统的仿真及与常规优化配置算法结果的比较,验证了本文方法的有效性和可行性。

1 常规串联电抗器优化配置模型

1.1 串联电抗器数学模型

常用串联电抗器配置方式有串接于分段母线联络线方式和串接于线路方式[10]。前者需要母线可分段运行,而且要求分段母线间留有一定空间用于安装串联电抗器。考虑到配置的通用性,本文一律采用串接于线路方式配置串联电抗器。

本文采用一种简化串联电抗器模型[5],假设与超标母线相连的线路上都加装了串联电抗器,其等值电路如图1所示。

当Δzij=0时,表示线路上未配置串联电抗器;当Δzij>0时,表示线路上配置了串联电抗器,串联电抗器阻抗为Δzij。

需要注意的是对于双回线和三回线而言,若需加装串联电抗器,则每一回线都需要加装串联电抗器,且阻抗值必须相同。在计算导纳矩阵时,采用文献[4]的方法,避免反复重构导纳矩阵,节约计算时间。

1.2 常规串联电抗器优化配置模型

常规串联电抗器优化配置模型可以表示为

其中,F(x)表示目标函数;G(x)表示等式约束;H(x)表示不等式约束。

1.2.1 目标函数

本文从经济性角度,通过配置串联电抗器的总成本来评价配置方案的优劣。投入总成本包括串联电抗器的设备成本和安装成本。其设备成本可以近似认为与其阻抗值成正比,而其安装成本可认为与其配置数量成正比,则目标函数表示为

其中,g1表示单台串联电抗器安装成本;g2表示单位Ω串联电抗器的制造成本;xi表示第i条线路安装的串联电抗器阻抗值;NCLI表示安装串联电抗器的台数。

考虑到实际配置成本会受到原材料成本及劳动力成本的影响,而单台串联电抗器安装成本与单位Ω串联电抗器的制造成本的比值(即g1:g2的值)变化很小,可以采用成本系数进行计算。实际成本就在成本系数和的基础上乘以当前单位成本系数的花费。

1.2.2 等式约束条件

等式约束条件主要是指系统状态量满足式(2)所示潮流方程:

其中,PG.i和QG.i分别表示i节点发电机输出的有功功率和无功功率;PL.i和QL.i分别表示i节点负荷有功功率和无功功率;N表示系统节点数。

1.2.3 不等式约束条件

常规优化配置串联电抗器的不等式约束条件包含两方面,一方面是短路电流约束,另一方面是母线电压约束,如式(3)所示。

其中,Ire.i和Isc.i分别表示i母线的限流目标值和当前短路电流值;Ui.min和Ui.max分别表示电网正常运行时i母线的电压允许最小值和最大值。

1.2.4 优化算法

通过式(1)—(3)将串联电抗器优化配置问题转换为非线性规划问题。由于很难求得短路电流对线路阻抗灵敏度的解析解,所以在应用诸如内点法等确定性优化算法时存在困难。一般采用现代优化算法进行求解,如粒子群算法、遗传算法等。对于本文所采用的模型,推荐采用粒子群算法,因为粒子群算法具有更好的边界搜索能力。

2 基于BCU法的暂态稳定裕度分析

直接法是一种快速高效的暂态稳定分析方法,同时能够对系统的稳定性进行定量分析。目前成熟应用的方法主要有EEAC法[11]和BCU法[12]。本文选择BCU法作为稳定裕度计算的主要方法。

2.1 BCU法简述

BCU法即基于稳定域边界的主导不稳定平衡点法,是一种快速有效的暂态能量裕度的计算方法。该方法是PEBS法和概念性主导不稳定平衡点法的结合,主要根据稳定域边界理论和原始系统与相应梯度系统平衡点相关性理论,其特点是对于主导不稳定平衡点的计算具有较好的收敛性[12,13,14,15]。

BCU法利用故障后的主导不稳定平衡点(CUEP)的能量作为临界能量Wcr,而故障切除时的系统能量为Wcl,此时系统暂态稳定裕度为

若稳定裕度ΔW≥0,则表示故障切除时系统稳定;若稳定裕度ΔW<0,则表示故障切除时系统稳定裕度不足,即失稳。

2.2 能量函数

在COI坐标系下,发电机采用经典二阶模型,负荷采用恒阻抗模型,若不计机械阻尼,发电机转子运动方程如式(5)所示。

其中,i,j=1,2,…,NG。

采用首次积分法构造暂态能量函数[15]。在积分耗散项时,为简化计算过程提高计算速度,采用线性路径加以近似。采用近似后虽然可能得到偏于冒进的错误结果,但是在实际应用中误差满足工程应用的要求[16]。此时能量函数如式(6)所示。

其中,θs表示故障后系统稳定平衡点(SEP)。

2.3 BCU法的实现步骤

BCU法的实现主要包括时域仿真、故障后SEP和故障后CUEP的计算。本文程序主要是在基于Matlab的PST软件包中进行开发,时域仿真采用PST自带程序,步长设为0.002 s。在计算故障后的SEP时,以故障前的SEP为初值代入方程(7)中,采用牛顿-拉夫逊法即可求得。

求解故障后CUEP是实现BCU法的关键,也是BCU法的精髓。首先根据在平衡点处原系统的特点,假设一直为0,构造原系统的收缩系统,如式(8)所示。

再通过求解收缩系统的CUEP,从而得到原系统的CUEP。求解收缩系统CUEP的步骤如下:

a.对故障后系统进行时域仿真,检测故障后轨迹出口点θEP;

b.以θEP为初值,积分收缩系统方程,计算最小梯度点θMGP;

c.以θMGP为初值,采用Newton法求解式(7),即可得到收缩系统的CUEP为θCUEP,此时的(θCUEP,0)即为原系统的CUEP。

详细计算步骤可以参考文献[15]。

3 考虑暂态稳定约束的串联电抗器配置优化算法

3.1 算法模型

不失一般性,选择系统预想故障集为任一线路首端或末端发生三相接地短路。则考虑暂态稳定约束的串联电抗器优化算法模型可以表示为

其中,目标函数F(x)、等式约束G(x)和不等式约束H(x)与常规优化配置算法模型相同,分别采用式(1)—(3);系统暂态稳定约束表示为ΔUk#k-j≥0,ΔUk#k-j表示k-j线路k侧母线发生三相短路故障,切除线路k-j后的系统暂态稳定裕度。

3.2 算法设计

直接对模型式(11)采用粒子群算法进行优化,由于需要对全网进行稳定性分析,计算量大,并不可行。考虑到配置串联电抗器虽然会对系统的稳定性带来影响,但是对不同线路稳定性的影响程度不尽相同,一般会对装设串联电抗器的线路及其周围线路影响较大。利用这一特点可以缩小稳定性分析的范围,设计思想是:首先采用不考虑系统暂态稳定性的优化配置串联电抗器算法,计算得到最优配置方案;然后利用BCU法计算系统所有线路的临界切除时间tcc。若所有线路满足tcc要求,则所得配置方案就是最终配置方案;若有线路不满足tcc的要求,则采用考虑系统暂态稳定约束的优化配置串联电抗器算法,并根据如下条件确定稳定约束所考虑线路及故障范围:

a.故障时tcc>0.3 s的线路不纳入稳定约束考虑范围,因为线路的稳定裕度已经足够充足,合理地加装串联电抗器对该线路稳定性的影响有限,不至于使其稳定裕度不足;

b.故障时稳定裕度不足的线路及此故障母线所连接的其他线路,不得与a冲突;

c.加装串联电抗器的线路。

采用这种设计思想可以减少稳定约束所需考虑线路及故障的范围,将稳定性分析的对象从全网所有线路缩小到所需关注的部分线路,减少了计算量,而且不会影响结果的准确性。为保证结果的准确性,可以在得到考虑稳定性配置方案后,再次校验全网稳定性,若出现新的线路稳定裕度不足,则将该线路纳入稳定裕度考虑范围,再次进行优化,以保证结果满足稳定性的要求。

3.3 算法流程

具体算法流程如图2所示。

4 算例分析

4.1 测试系统

为更好地测试所提方法限制短路电流的效果,在New England 10机39节点标准测试系统[17]的基础上对其线路电抗值进行调整。调整后系统结构如图3所示,相应参数调整线路列于表1。

设定断路器最大遮断电流为30 k A,测试系统中三相短路电流超标母线及电流值列于表2。

kA

4.2 优化配置

算例限流目标设为30 k A,目标函数中设置g1=2,g2=1。根据规定,系统故障临界切除时间需满足tcc≥0.1 s,考虑到BCU算法可能出现的误差以及tcc需留有一定裕度,设定tcc≥0.14 s时系统满足稳定约束要求。采用文献[5]中提出的串联电抗器优化配置算法限制表2中超标母线短路电流。所得配置方案(即方案1)列于表3。

采用BCU法计算根据方案1加装串联电抗器后系统的暂态稳定裕度。结果表明线路15-16和16-21在近母线16侧发生三相短路故障时,tcc<0.14 s,不满足先前设定的稳定裕度要求。采用本文算法优化配置串抗器限制短路电流,所得结果(方案2)列于表4。表5和表6分别列出根据方案2加装串联电抗器后系统的短路电流和暂态稳定裕度,此时系统满足正常运行的要求。

方案2在满足限制短路电流的前提下,使所有线路稳定裕度都满足设定的要求。从配置方案角度分析,方案2需要在3条线路中加装串联电抗器,总加装电抗值为5.25Ω,总投入成本系数为11.25。方案1则需要在2条线路中加装串联电抗器,总加装电抗值为5.5Ω,总投入成本系数为9.5。相比于方案1,方案2总加装电抗值减少了4.5%,但由于多配置了一台串联电抗器造成配套设备成本及施工成本增加,使得总投入成本增加了18.4%。通过分析可知,本文算法所得方案既满足限制短路电流要求也满足暂态稳定裕度要求,虽然可能会造成总投入成本的增加,但是能够保证系统稳定运行,具有应用价值。

5 结语

暂态稳定约束 篇3

电力市场环境下,仅考虑静态安全约束的最优潮流(optimal power flow,OPF)使得系统运行点更容易接近稳定极限,已无法满足现代电力系统的要求,考虑暂态稳定约束的OPF(transient stability constrained OPF,TSOPF)成为新的研究方向。

OPF是一个大规模的非线性优化问题,内点法及其改进算法由于良好的计算效率和收敛特性在该问题求解中得到了广泛应用[1,2,3,4]。考虑暂态稳定约束后,增加了大量的微分代数方程约束,是一类特殊的无穷维优化问题[5]。目前关于TSOPF的研究集中于暂态稳定约束条件的处理。文献[6]将发电机暂态微分代数方程差分后作为等式约束加入到OPF;文献[7]则考虑了多个预想故障,并引入既约导纳矩阵以消去网络方程;文献[8]采用大步长对预想故障集进行故障预选和同调识别,对于小步长精确计算仅考虑“起作用”故障;文献[9,10]采用约束转换技术,通过欧氏空间变换将函数空间的TSOPF转换为相同规模的静态优化问题;文献[11]将该问题分解为最优潮流和最优控制2个子问题,通过交替求解这2个子问题得到原问题的解;文献[12]将暂态能量裕度作为附加约束引入到OPF模型中;文献[13]利用暂态稳定裕度指标及其灵敏度组成的不等式代替暂态稳定约束。逐次线性规划[6]、内点法[7]以及智能优化算法[14]均被成功应用于该问题的求解。

常规方式下,基于时域仿真差分思想的TSOPF将发电机转子运动方程离散化,并作为OPF问题的等式约束。由于任何差分方法都具有截断误差,将其作为等式约束过于苛刻。内点法求解过程中,这种处理方式使得修正方程组的阶数随着系统规模和故障个数急剧上升,导致计算耗时长、内存消耗大甚至不可解的“维数灾”。针对该问题,本文将描述系统暂态过程的差分方程转换为不等式约束,其上下限与所采用的差分方法精度相关,从而大大降低了内点法计算过程中修正方程组的阶数。对3个不同规模的系统进行测试,结果显示与常规方式相比,该方法计算时间明显缩短,占用内存减少,能够求解更大规模系统的TSOPF问题。

1 TSOPF问题描述

考虑多个预想故障的TSOPF问题可表示为:

$\{\begin{array}{l}\underset{\mathit{x}}{{\displaystyle \min }}f(\mathit{x})\\ \text{s}.\text{t}.\mathit{h}(\mathit{x})=0\\ \underset{¯}{\mathit{g}}\le \mathit{g}(\mathit{x})\le \overline{\mathit{g}}\end{array}(1)$

式中:f(x)为该问题的目标函数,为不失一般性,本文采用发电费用最小[7]为目标函数;h(x)为等式约束向量,包括潮流方程、发电机初值约束和暂态微分代数方程离散化后得到的差分方程;g(x)为不等式约束向量,包括发电机有功和无功出力、变压器变比、节点电压和线路有功潮流上下限等静态安全约束以及暂态稳定约束;$\overline{\mathit{g}}$和$\underset{¯}{\mathit{g}}$分别为不等式约束上下限。详细的问题模型可参阅文献[7,8]。

TSOPF问题的等式约束包含大量暂态差分方程。常规方式下,若发电机采用经典模型,负荷用恒定阻抗表示,由隐式梯形法可得如下差分方程[8]:

$\{\begin{array}{l}\text{Δ}\delta {}_i^t(k)=\delta {}_i^t(k)-\delta {}_i^{t-1}(k)-\frac{\omega {}_0\text{Δ}t}{2}(\omega {}_i^t(k)+\\ \omega {}_i^{t-1}(k))=0\\ \text{Δ}\omega {}_i^t(k)=\left(\frac{D{}_i\text{Δ}t}{2{\rm M}{}_i}+1\right)\omega {}_i^t(k)+\left(\frac{D{}_i\text{Δ}t}{2{\rm M}{}_i}-1\right)\cdot \\ \omega {}_i^{t-1}(k)-\frac{\text{Δ}t}{2{\rm M}{}_i}(2{\rm P}{}_{\text{G}i}-{\rm P}{}_{\text{e}i}^t(k)-\\ {\rm P}{}_{\text{e}i}^{t-1}(k))=0\end{array}(2)$

式中:i∈nG,k∈nC,t∈nt,nG,nC,nt分别为发电机、预想故障和积分时段数目;ω0为同步角速度有名值;δti(k),ωti(k)分别为发电机i第k个故障第t时段的功角和角速度相对同步转速的偏差;Δt为积分步长;Di,Mi,PGi分别为发电机i的阻尼系数、惯性常数和有功出力;Ptei(k)为发电机i的电磁功率,可采用既约导纳矩阵表示以消去网络方程[7]。

2 降阶内点算法

基于内点法的TSOPF计算其主要任务是求解对称非正定线性方程组:

$\left[\begin{array}{cc}{\mathit{{\rm H}}}^{\prime }& \nabla {}_{\mathit{x}}\mathit{h}(\mathit{x})\\ \nabla {}_{\mathit{x}}^{\text{Τ}}\mathit{h}(\mathit{x})& 0\end{array}\right]\left[\begin{array}{c}\text{Δ}\mathit{x}\\ \text{Δ}\mathit{y}\end{array}\right]=\left[\begin{array}{c}\mathit{L}{}_{{\mathit{x}}^{\prime }}\\ -\mathit{L}{}_{\mathit{y}}\end{array}\right](3)$

式中:x为控制和状态变量;y为等式约束乘子;H′为障碍矩阵和目标函数、等式及不等式约束对应海森矩阵的线性组合[7];Lx′和Ly为库恩—塔克(KKT)系统的右端残差向量,其定义可参看文献[15]。

式(3)的求解占用了绝大部分计算时间[7],其阶数将直接影响到内点法的计算耗时。由式(3)可知,系数矩阵的阶数取决于等式约束和变量的数目,与不等式约束个数无关。如果设法减少等式约束个数,则可降低式(3)的阶数。常规方式下,差分TSOPF等式约束中绝大部分为暂态差分方程。事实上,由于任何数值算法都存在误差,将这些差分方程作为严格等式约束并无必要。对于初值问题:

${\mathit{v}}^{\prime }=\text{ϕ}(t,\mathit{v})t\ge t{}_0,\mathit{v}(t{}_0)=\mathit{v}{}_0(4)$

采用隐式梯形法差分化,第j步到第j+1步产生的截断误差为[16]:

v(tj+1)-$\mathit{v}(t{}_j)+\frac{\text{Δ}t}{2}(\text{ϕ}(t{}_j,\mathit{v}(t{}_j))+\text{ϕ}(t{}_{j+1},$

v(tj+1)))$=-\frac{(\text{Δ}t){}^3}{12}{\mathit{v}}^{‴}(t{}_j)+{\rm O}((\text{Δ}t){}^4)(5)$

即隐式梯形法是二阶精度的差分方法,存在3阶局部截断误差。实际上,任何一种差分方法都存在误差,将差分方程作为等式约束是非常苛刻的条件,直接导致式(3)阶数增加,带来计算时间和内存消耗过大等问题。在大规模系统、多预想故障的TSOPF求解中,这些问题更加突出。若将这些差分方程转换为不等式,则可以大大减少等式约束个数,即式(2)转换为:

$\{\begin{array}{l}-\beta (\text{Δ}t){}^n\le \text{Δ}\delta {}_i^t(k)\le \beta (\text{Δ}t){}^n\\ -\beta (\text{Δ}t){}^n\le \text{Δ}\omega {}_i^t(k)\le \beta (\text{Δ}t){}^n\end{array}(6)$

式中:Δt为积分步长;n为所用差分方法的精度阶数。

差分不等式约束的上下限设置为±β(Δt)n,其中β为约束系数,β∈[0,1]。在内点法求解过程中,将1个差分方程由等式约束转换为上下限不等式约束,可减少1个等式约束及其乘子(y),同时增加2个松弛变量(l,u)以及2个不等式约束乘子(z,w)。由于松弛变量和乘子通过直接回代求解,计算耗时很少,故这些变量的增加对优化过程的总耗时影响甚小,测试结果也将证明这一点。

如图1所示,假设t0时刻系统中某处发生故障,tcr时刻故障被切除,tmax时刻暂态仿真过程结束。以tb表示暂态过程中等式约束和不等式约束的分隔时刻,即t0～tb时段的暂态差分方程视为等式约束,tb～tmax时段的暂态差分方程视为不等式约束。以隐式梯形法为例,说明本文采用的差分方程不等式化降阶过程。

步骤1:发电机采用经典模型,负荷视为恒定阻抗模型,采用隐式梯形法将发电机转子运动方程差分化,得到差分方程组。

步骤2:以tb为界限,将差分方程组分别转换为等式约束集S1和不等式约束集S2:

$\begin{array}{l}S{}_1=\{\begin{array}{l}\text{Δ}\delta {}_i^t(k)=0\\ \text{Δ}\omega {}_i^t(k)=0\end{array}t{}_0\le t\le t{}_{\text{b}}(7)\\ S{}_2=\{\begin{array}{l}-\beta (\text{Δ}t){}^n\le \text{Δ}\delta {}_i^t(k)\le \beta (\text{Δ}t){}^n\\ -\beta (\text{Δ}t){}^n\le \text{Δ}\omega {}_i^t(k)\le \beta (\text{Δ}t){}^n\end{array}t{}_{\text{b}}<t\le t{}_{\max }(8)\end{array}$

式中:tb∈[t0,tmax]。

如考虑到故障时段对整个暂态过程的重要影响,可令tb=tcr。若取tb=t0,即将所有差分方程处理为不等式约束;若取tb=tmax,即常规方式。对于隐式梯形法,精度阶数n为2。

步骤3:将S1和S2分别作为TSOPF问题的等式约束和不等式约束,采用内点法求解降阶后的TSOPF问题。

由于式(5)中v(tj)的计算比较困难,故差分不等式上下限中引入约束系数β。TSOPF问题中,各积分时段的δt, ωt都与前一时刻的δt-1,ωt-1有关,因此β的取值关系到降阶内点算法的精度。

3 算例分析

3.1 测试算例

本文采用预测—校正内点法对3个算例进行了计算。表1给出了3个测试算例的参数,其中nB,nL分别为系统节点和线路数目,t0=0 s。对于所有算例,均采用第2节提出的处理方式,分别取tb=tcr和tb=t0进行计算,并与常规处理方式所得结果进行对比。本文中也就β和tb的取值对计算耗时和目标函数值的影响进行了分析。程序采用各发电机功角相对于惯性中心(center of inertia,COI)的偏差不超过±100°作为暂态稳定判据。为简化计算,假定变压器变比恒定。计算中最大迭代次数设定为50。采用隐式梯形法进行暂态稳定分析时,积分步长可取0.01 s～ 0.02 s,甚至更长[15],本文采用0.02 s的积分步长,式(8)中精度阶数n为2。

采用32位MATLAB R2009a编程求解TSOPF问题,并用MATLAB集成的常微分方程解法器ODE45对程序求得的摇摆曲线进行验证。仿真环境为Intel Core 2 Duo E8400, 4 GB内存。程序已充分考虑矩阵稀疏性,且计算过程中已自动对稀疏矩阵进行优化编号以减少新增非零注入元。

3.2 测试结果

为叙述方便,称tb=tmax的常规方式为TSO1,tb=tcr和tb=t0的降阶处理方式分别为TSO2和TSO3。表2给出了上述3个算例在各处理方式下TSOPF问题的规模和计算结果。

表2中,re为等式约束数,Nd为线性方程组(3)的阶数,Pr为阶数百分比,Ni为内点法迭代求解次数,Tc为求解式(3)的总耗时,Ta为各次迭代求解式(3)的平均耗时,Tb为各次迭代回代求解及更新松弛变量和乘子的平均耗时,Vobj为目标函数值,Dobj′为当前目标函数值Vobj相对于常规方式下目标函数值VTSO1obj偏差的百分比,即

$D{}_{\text{o}\text{b}\text{j}}´=\frac{|V{}_{\text{o}\text{b}\text{j}}-V{}_{\text{o}\text{b}\text{j}}^{\text{Τ}\text{S}\text{Ο}1}|}{V{}_{\text{o}\text{b}\text{j}}^{\text{Τ}\text{S}\text{Ο}1}}\times 100％(9)$

附录A图A1～图A3给出了上述3个算例在各处理方式下所得的摇摆曲线,可知,在给定的预想故障下,时域仿真得到的摇摆曲线与TSOPF计算所得摇摆曲线基本一致,说明该降阶方法所得最优潮流结果满足预先设定的暂态稳定约束条件。表2中数据Dobj′说明降阶后所得目标函数值与常规方式的偏差在可接受的范围内。

对于所有算例,比较表2中各处理方式下的测试结果Nd,Pr,Tc和Ta可知,降阶后各算例中修正方程组的阶数降低至接近原来的一半,迭代求解该线性方程组的总耗时和平均耗时均大大减少,说明该降阶内点算法对所考虑的算例均具有明显的省时效果。其中对于规模较小的系统(IEEE39),由于常规方式的计算耗时本来就不多,该算法的省时效果并不明显;对于较大规模系统(CASE162),该算法较常规方式节省的时间在10倍左右;对于更大规模的系统(CASE300),常规方式下由于内存限制该问题已不可求解,而经过差分不等式降阶处理后也能够在较短时间内给出优化结果。这些都充分说明了该降阶内点算法在计算耗时和内存耗用量方面的优势。

表2的测试结果显示,每次迭代过程中回代求解松弛变量和乘子的时间Tb相对于求解线性方程组(3)的时间Ta非常少。对于所有算例,TSO2和TSO3对应的Tb与TSO1相比增加很少,说明降阶后由于松弛变量和不等式约束乘子增加带来的计算耗时增长有限,从而印证了第2节的结论。

3.3 β取值分析

图2给出了各算例取tb=tcr时,β取不同值对计算耗时和目标函数值的影响。其中,T1为β取不同值时各算例中求解式(3)的总耗时;T2为表2中各算例TSO2方式下的总耗时;Dobj为当前目标函数值Vobj相对于各算例TSO2方式下目标函数值VTSO2obj的偏差百分比,其定义为:

$D{}_{\text{o}\text{b}\text{j}}=\frac{|V{}_{\text{o}\text{b}\text{j}}-V{}_{\text{o}\text{b}\text{j}}^{\text{Τ}\text{S}\text{Ο}2}|}{V{}_{\text{o}\text{b}\text{j}}^{\text{Τ}\text{S}\text{Ο}2}}\times 100％(10)$

从图2(a)可知,β取值小使得总的计算时间增加(图中β=10-5时,程序达到迭代次数上限,各算例均存在收敛性问题)。图2(b)说明β过大将导致暂态稳定约束过于松弛,目标值偏差较大,优化结果不理想。结合图2(a)和图2(b)可知,对于所考虑的算例,采用隐式梯形法差分化后,β取0.001～0.020可在计算耗时和暂态仿真及目标函数值精度之间取得较好的折中。

表3给出了在所建议的β取值范围内,各算例通过ODE45进行暂态仿真所得发电机的最大相对功角δmax和目标值偏差百分比Dobj′(如式(9)所示)。

从表3结果看出,差分不等式降阶处理所得的暂态仿真和目标函数值精度均在可接受的范围内,说明这种将差分方程部分或全部作为不等式约束的方法保证了足够的计算结果精度,且节省了大量计算时间,因此是可行且有效的。

3.4 tb取值分析

图3给出了不同tb时各算例计算总耗时的变化曲线,其中β按表2所示各算例TSO2方式下取恒定值,T1为tb取不同值时各算例中每次迭代求解式(3)的平均耗时,T2为表2中各算例对应TSO2方式下的平均耗时。在所考虑的算例中,各故障的切除时间tcr均为0.2 s。

从图3看出,随着tb的增大,由于修正方程组阶数不断增大,计算耗时总体上是不断增加的,对于大规模系统则更加明显。这说明对于同一个算例,修正方程组的阶数是影响计算耗时的一个重要因素,从而进一步证明了本文所提出的降阶内点算法的合理性。当tb≤tcr时,计算耗时不会明显增加。注意对于算例CASE300,tb>0.2 s后该问题的求解由于内存不足无法完成。

4 结语

本文提出了一种多预想故障TSOPF问题的降阶内点算法。该方法考虑了差分方法的截断误差,将暂态差分方程转换为不等式,其上下限与所用差分方法的精度相关,并采用内点法求解降阶后的TSOPF问题。对3个不同规模系统的测试结果证明该降阶算法能有效减少计算时间,且测试结果显示出以下特点:

1)降阶内点算法中,线性方程组阶数比常规方式明显减小,其阶数接近常规方式的50%。对于不同规模的系统均可节省大量计算时间,甚至能够求解常规方式下不能求解的问题。

2)测试证明,β取值过小可能导致总的计算时间增加;β取值过大则导致优化结果不理想。若采用隐式梯形法差分化,对于所有算例β,取0.001～0.020时能够取得满意的结果。

3)测试证明,tb取值越大,修正方程组阶数越高,其求解耗时总体上增加,对于大规模系统尤其如此。对于本文的算例,tb≤tcr是一个较合理的取值。

附录见本刊网络版(http://www.aeps-info.com/aeps/ch/index.aspx)。

暂态稳定约束 篇4

中国将在目前大区和省独立电网基础上形成以超高压、直流输电为基础的多区互联的全国统一电网,以实现资源互补和优化。物理上将出现长距离、高电压、重负荷、多区域弱联系的交直流混合输电系统,而且在系统运行方面还将出现潮流多变、重负荷运行、靠近稳定运行极限等新的问题。对于这样的复杂系统,为了能够分析及仿真系统在正常工况、预想事故集下的动态行为,找出全系统或某一个区域的安全裕度,并将此信息快速、准确、直观地呈现给调度人员,对暂态稳定约束下的区域电网之间最大输电能力(TTC)的计算开展深入、全面的理论和应用研究是非常必要和迫切的。

针对预想故障集,暂态稳定约束的TTC评估本身计算强度大、时间长,依靠串行计算只能实现离线分析,无法满足在线分析的要求。近几年,随着高性能计算机软、硬件技术的高速发展,以及价格的不断下降,为研究开发新的分布式并行算法来实现大规模电力系统暂态稳定约束TTC评估提供了很好的机会。目前,国内外有关暂态稳定约束TTC评估的研究文献[1,2,3,4,5,6,7,8]主要是针对处理暂态稳定约束方面;研究范围大多局限于单个预想故障下、挖掘暂态稳定计算内在的并行性的暂态稳定约束TTC评估。值得一提的是,无论采取哪种并行算法都存在加速比饱和的问题[9]。

本文针对预想故障集下求解大规模互联电力系统的暂态稳定约束TTC问题,介绍了暂态稳定约束TTC的计算流程;提出了一种新颖的基于优化的工作站群计算方法的设计思想,并对其进行了实现。

1 暂态稳定约束TTC计算流程

暂态稳定约束TTC评估是在紧急故障集或预想故障(集)下,针对一种给定的发电—负荷增长模式(可定义为通过改变发电、负荷的量来满足特定的系统功率调度准则)来确定所研究系统的稳定裕度。

暂态稳定约束TTC评估的计算流程如附录A所示。在计算过程中,暂态稳定约束是基于时域仿真的方法,并利用一个基于功角的稳定判别指标[8]来判断系统是否暂态失稳;并由给定的发电—负荷增长模式计算出系统的稳定裕度。

从流程图可以看出,计算一次暂态稳定约束TTC需要多次迭代,而每次迭代都需要进行一次时域仿真,采用串行方式计算每个故障集下的暂态稳定约束TTC难以满足在线的要求。由此,如采用分布式并行计算模式,将多台计算机的计算能力加以综合,提高整体计算速度,从而达到在线应用具有重要意义。

2 优化的工作站群的计算方法设计与实现

2.1 分布式并行计算方法分析

电力系统中,传统意义上的并行计算一般分为单机并行计算及分布式并行计算。单机并行计算一般采取开启多线程并行计算方式,计算平台一般为单机多CPU系统。这种并行计算方法能够充分发挥高性能计算机的性能,但不能把多台计算机资源综合利用,扩展性较差。分布式并行计算方法一般采用多机分摊的任务分配策略[9],或根据电力系统特有的分层分区特性,采用基于区域故障集的任务分配策略[10]。计算平台一般为同构的PC机集群,这2种计算模式克服了单机并行计算模式扩展性差的缺点,综合利用了多台计算机的计算能力,在一定程度上可以满足在线计算的要求。然而,实际应用中,可利用的计算资源大多为异构的PC机,这些PC机的性能是不一样的,甚至区别很大。若利用异构的PC机构成一个分布式计算系统,不管采用多机分摊还是采用基于区域故障集的任务分配策略,系统总的计算耗时取决于系统中最慢计算节点对故障的处理时间。性能强的计算机势必得不到充分利用,浪费了本来可利用的计算资源。此外,同构的PC机集群系统构成的分布式计算平台,当遇到分布式任务中的子任务计算执行时间事先难以确定时,采用多机分摊的任务分配策略同样会有计算资源得不到充分利用的问题。

2.2 计算方法的设计

针对上述问题,本文提出了一种新颖的分布式并行计算方法。该方法是基于异构分布式并行计算系统,系统中的各台计算机拥有不同的计算性能;通过采用一种“能者多劳”的任务调度策略,最大限度地挖掘每个计算节点的计算能力。

该分布式计算模式包含2种类型的节点:①计算节点:装有计算引擎,实现子任务的计算并把计算临时结果发送给计算管理器节点;②计算管理器节点:作为一个中央站点提供任务的划分、任务调度和汇总等功能。优化工作站群的计算方法设计框图如图1所示。

根据图1,所设计的计算模式基本执行流程如下:

1)一个可利用的计算节点发送计算请求到计算管理器;

2)计算管理器响应请求并创建一个工作线程与其交互;

3)计算管理器分派任务给工作线程;

4)工作线程发送分配的任务给各计算节点;

5)各计算节点完成计算任务,发送计算结果到计算管理器;

6)各计算节点通知计算管理器并等待工作线程分派下一个任务。

从基本执行流程可知,计算节点处理的任务数与其计算能力是成正比的。在特定情况下,一些高性能的计算节点(比如多CPU硬件配置),可以同时开启多个计算引擎。从以上分析可知,每台计算机各尽所能,计算能力得到充分挖掘,从而有效达到了系统负载动态均衡。

此外,需要考虑一种特殊情况:如某些计算节点因为自身计算能力的不足,导致一直未能完成所分配的任务(或出于网络故障的原因),与此同时,其他计算节点处于空闲状态(此时所有的任务已分配完)。在这种情形下,提出了一种称之为条件资源抢占的计算策略,在一定条件下,如计算任务超时,空闲的计算节点会抢占未完成的任务并进行计算,从而使整个工作站群的计算效率得到进一步提高。

2.3 计算方法的实现

由2.2节分析可知,整个计算方法的核心在于计算管理器;在实现计算管理器过程中,本文采用了多进程与多线程编程技术[11]。根据计算管理器的职能,将计算管理器分解成多个进程对象,每个进程对象独立运行,负责各自的职能,某些进程对象会创建多个线程对象,它们之间可通过全局变量、信号、共享内存等方式进行通信。

计算管理器的逻辑结构由初始化进程、协调线程、监听线程、数据存储线程、工作线程链表以及计算进程(TTC评估)等模块组成。附录B给出了计算管理器逻辑结构流程图。

该计算管理器中的各个模块的特点与功能描述如下:

1)初始化进程:负责系统的初始化工作;

2)协调线程:负责接收来自中心管理器的潮流数据、故障扫描及生成计算任务数组并通知所有空闲线程进行计算;

3)监听线程:负责监听来自计算节点的连接,为每一个连接创建一个工作线程与连接对应的计算节点进程通信,监视计算节点;

4)工作线程:与计算节点进行通信,发送计算任务、接收计算结果,监视计算节点状态;

5)数据存储线程:将计算结果存储到数据库(当计算结果不需要通过计算管理器直接存储到数据库时,不需要安装);

6)计算进程:为安装在计算管理器上的计算引擎(当计算管理器不参与分布式并行计算时,不需要安装)。

工作线程链表、数据存储线程以及计算进程在计算管理器中是可选的,也就是说,可以根据实际的需要或根据服务器的性能进行组合。

3 计算方法的扩展性与容错性

3.1 扩展性

分布式并行计算平台的计算节点数越多,可以同时计算的子任务数就越多,整个平台的计算能力就越强。为了高效地利用分布在Internet上的大规模计算资源,分布式并行计算平台的扩展性是必须要考虑的。扩展性通常分为静态扩展和动态扩展。静态扩展指系统增删节点时,需要先停止当前运行的系统,然后根据节点的增删情况进行重新配置,最后重新启动整个系统。动态扩展也称为在线扩展,系统增删节点时能够自动适应这种变化,自动完成资源的重新配置,在无人干涉的情况下进行自我管理与自我维护。良好的静态扩展性是分布式计算方法具备的一般特征,但并不是所有的分布式计算方法都具备良好的动态扩展性。

本文提出的基于优化的工作站群的分布式并行计算方法把计算管理器作为服务端,计算节点作为客户端,采用了一个服务器对应多个客户端的架构,并且在进行暂态稳定约束TTC评估中,子任务之间不需要进行交互,当一个计算节点离开系统时,计算管理器能够自动把该计算节点还没有计算完的子任务置为未处理状态并转移给其他计算节点进行处理,从而使整个系统的分布式任务结果不会受到影响。

此外,作为客户端的计算节点在系统运行时可以动态地加入与退出,不会中断分布式任务的执行。因此,本文所设计的分布式并行计算方法除了具有良好的静态扩展性一般特征之外,还具有良好的动态扩展性。

3.2 容错性

为了保证计算平台实际运行的可靠性,本文在系统实现过程中使用了TCP/IP编程模型、FTP(文件传送协议),TCP/IP及FTP中的重传机制确保了数据传输的可靠性。

基于优化的工作站群分布式计算方法采用了子任务级的容错机制,这样做的好处是在不降低系统性能的情况下将容错过程简化,降低容错过程中的系统开销。该容错机制假定计算管理器是比较稳定的节点,而计算节点是不稳定的(可能随时会加入或退出系统),所以容错的处理主要是针对计算节点。当计算节点在计算子任务的过程中出现故障,无法继续工作时,计算管理器会及时发现此情况,同时置该计算节点处理的子任务为未处理状态并转移给空闲的计算节点,在一定程度上避免了因为计算节点出现故障而导致系统无限等待的情形。

在提高计算平台的安全性方面,为了防止外来系统对计算管理器的非法访问,本文采用了客户端认证机制,每个计算节点加入到系统中均需要在计算管理器进行合法性检查后方能进入到分布式系统。

4 算例分析

4.1 算例系统

应用本文提出的分布式并行计算方法对中国某实际电网进行了测试。该实际系统中有1 101个母线、1 485条交流线路、6条高压直流输电(HVDC)线路、97台发电机、387个负荷、5个区域。针对系统的50个单线路故障(500 kV节点之间线路)、49个切机故障、19个双回线故障(同杆架设),共118个故障进行暂态稳定约束TTC计算。其中,仿真步长为0.01 s,仿真时间为5.0 s。

4.2 计算环境

实验室中已建有由5台计算机构成的100 Mbit/s以太网络架构;该网络中各计算资源配置见附录C。

4.3 实验与性能测试

本文采用了3种分布式并行计算方案进行测试。

方案1:仅1台Xeon机参与计算,同时该计算节点开启2个计算线程进行并行计算。

方案2:采用3台Xeon机,同时每机分别开启2个计算线程,即共有6个计算线程参与计算。

方案3:实验室100 Mbit/s以太网络中的5台计算机(即3台Xeon机,2台Pentium机)均参与计算。

3种方案中均采用1台Pentium机作为计算管理器,需要指出的是,在方案3中这台Pentium机还同时作为计算节点参与计算。

为了测试所设计计算平台的计算稳定性,独立运行各方案10次,相应的计算时间如图2所示。

图2中各曲线波动性很小,说明该系统的稳定性非常好,此外,随着计算资源的不断加入,计算时间逐渐缩短。各方案平均计算时间见附录C。

3台同构计算机与1台同构计算机相比较,取得了很高的加速比;而2台Pentium机加入到3台Xeon组成的同构分布式计算环境中,形成由5台异构计算机组成的分布式计算环境,2台性能较差的Pentium机没有影响到系统的整体性能,而是分担了部分计算任务,进一步缩短了整体计算时间。由此可见,本文提出的基于优化的工作站群的分布式并行计算方法具有较高的并行效率,且能充分利用异构计算资源。方案3中各台计算机分摊的任务数如图3所示。图中X1,X2,X3 表示Xeon类型机的第1、第2和第3台机;P1和P2表示Pentium类型机的第1和第2台机。

从图3中的统计结果可知,由于Xeon的性能优于Pentium机,因此该类型的计算节点所承担的任务数量明显多于Pentium类型的计算节点,从而体现了“能者多劳”的任务调度策略。同类型计算节点承担的任务数存在较小的波动是因为故障类型的区别导致所需计算时间不同,如Pentium机节点计算1个单线路或切机故障的TTC大约需要40 s,而计算1个双回线路故障的TTC只需要3 s,因为该N-2线路故障只经过一次迭代计算就导致系统失稳,计算终止。

表1给出了各计算节点的计算时间与通信开销。根据表1,可以看出:Xeon机计算节点开启2个线程进行并行计算,线程之间相互抢占网络带宽,因此,Xeon机计算线程的通信开销较Pentium机要大,但是由于Xeon机计算节点采用多线程并行计算,实际上更为充分地利用了计算资源,提高了计算效率。由于计算线程与计算管理器同时运行于P1机之上,所以P1计算节点的通信开销接近0。总之,各计算线程的通信开销占总的处理时间比例相当小,实际运行中可以忽略不计。

5 结语

在异构计算资源广泛存在的背景下,本文采用“能者多劳”的任务调度策略,提出了基于优化的工作站群的分布式并行计算方法用于预想故障集下的电力系统暂态稳定约束TTC的评估。详细介绍了该计算方法的设计思想与实现方法。实验结果表明,该计算方法运行于同构分布式系统中时,能够得到很高的并行加速比;运行于异构分布式系统中时,能够充分利用异构计算资源,真正体现了“能者多劳”的任务调度策略,系统的整体计算性能得到了优化。所设计的分布式并行计算方法能方便地集成到现有的电力系统动态安全评估系统[11]中实现在线应用。

附录见本刊网络版(http://www.aeps-info.com/aeps/ch/index.aspx)。

参考文献

[1]刘明波,夏岩,吴捷.计及暂态稳定约束的可用传输容量计算.中国电机工程学报,2003,23(9):28-33.LI U Mingbo,XI A Yan,WU Jie.Calculation of available transfer capability with transient stability constraints.Proceedings of the CSEE,2003,23(9):28-33.

[2]李国庆,沈杰,申艳杰.考虑暂态稳定约束的可用功率交换能力计算的研究.电网技术,2004,28(15):67-71.LI Guoqing,SHENJie,SHEN Yanjie.Study on calculation of available transfer capability considering transient stability constraints.Power System Technology,2004,28(15):67-71.

[3]樊纪超,余贻鑫.考虑暂态稳定约束的临界割集功率传输和发电经济调度.电力系统自动化,2003,27(17):4-7.FANJichao,YU Yixin.Transient stability-constrained power transfer limit of a critical cutest and generation economic dispatching.Automation of Electric Power Systems,2003,27(17):4-7.

[4]杨新林,孙元章,王海风.考虑暂态稳定约束极限传输容量的计算方法.电力系统自动化,2004,28(10):29-33.YANG Xinlin,SUN Yuanzhang,WANG Haifeng.An optimization approach to assess stability-constrained total transfer capability.Automation of Electric Power Systems,2004,28(10):29-33.

[5]ZHANG X,SONG Y H,LU Q,et al.Dynamic available transfer capability(ATC)evaluation by dynamic constrained optimization.IEEE Trans on Power Systems,2004,19(2):1240-1242.

[6]OU Y,SI NGH C.Assessment of available transfer capability and margins.IEEE Trans on Power Systems,2002,17(2):463-468.

[7]薛禹胜,徐泰山,刘兵,等.暂态电压稳定性及电压跌落可接受性.电力系统自动化,1999,23(14):4-8.XUE Yusheng,XU Taishan,LI U Bing,et al.Quantitative assessments for transient voltage security.Automation ofElectric Power Systems,1999,23(14):4-8.

[8]AJJARAPU V,CHRISTY C.The continuation power flow:a tool for steady state voltage stability analysis.IEEE Trans on Power Systems,1992,7(1):416-423.

[9]郭琦,张伯明,王守相,等.基于计算机集群的电力系统暂态风险评估.电网技术,2005,29(15):13-18.GUO Qi,ZHANG Boming,WANG Shouxiang,et al.Personal computer cluster based power systemtransient risk assessment.Power System Technology,2005,29(15):13-18.

[10]薛巍,舒继武,王心丰,等.电力系统暂态稳定仿真并行算法的研究进展.系统仿真学报,2002,14(2):177-182.XUE Wei,SHU Ji wu,WANG Xinfeng,et al.Advance of parallel algorithm for power system transient stability simulation.Journal of System Simulation,2002,14(2):177-182.

暂态稳定约束 篇5

智能电网的目标是实现电网运行的可靠、安全、经济、高效、环境友好和使用安全。电力系统是一个动态系统,其安全问题的本质必然是动态的[1]。常规的最优潮流(Optimal Power Flow,OPF)难以保证电网运行的动态安全性,不能满足智能电网的要求。暂态稳定约束最优潮流(Transient Stability Constrained Optimal Power Flow,TSCOPF)是一种兼顾安全性与经济性的暂态稳定预防控制实现方法,得到了研究人员的广泛关注[2,3,4,5,6,7,8,9,10,11,12,13,14]。

TSCOPF是在常规OPF中引入一组等式和不等式约束:等式为系统的转子运动方程以及初值方程;不等式为功角稳定性约束。实际上,TSCOPF的约束条件为一组微分代数方程。对微分方程的处理目前主要有两种方法:第一种,将微分方程离散化为一系列的代数方程,直接加到OPF中[2,3,4,5,6,7];第二种,采用约束转换技术将TSCOPF问题化为同规模的OPF问题[8,9,10]。第二种方法虽然可以避免第一种方法中出现规模庞大的问题,但得不到原问题的最优解。现代内点法解优化问题主要的计算量是解修正方程。许多大型的电力系统优化问题都利用稀疏技术及节点优化编号来减少注入元以提高求解速度[15],提升的空间有限。文献[6]在第一种方法的基础上,采用预测校正内点法与简约空间技术相结合的算法对TSCOPF问题进行求解,缩减了求解修正方程所耗用的时间,有效地缓解了上述问题。

简约空间技术是一种高效的线性方程求解器,最初被用于化学工程及流程优化等领域[12,13]。简约空间技术按照一定的分解策略对修正方程进行分解,通过缩小修正方程的规模,减少解修正方程所耗用的时间,从而达到提高计算效率、降低内存损耗的目的。简约空间技术适合解自由度小的问题,其计算复杂度主要取决于问题自由度的大小,即控制变量与总原始变量的比值的大小。离散化后的TSCOPF问题虽然规模庞大,但其控制变量的个数相对很小,且不随故障重数、仿真时长和仿真步长的改变而改变,因此适于用简约空间技术求解。

转子运动方程包含了电磁功率方程,大部分文献中的电磁功率方程采用混合极坐标[2,3,4,5,6,7,9,11,12]。与直角坐标及混合极坐标形式相比,极坐标形式的模型最为简洁,便于记忆和编程。本文在文献[6]的基础上,建立了新的基于极坐标形式的TSCOPF模型,即潮流方程及电磁功率方程均采用极坐标。统一的坐标系,降低了程序的编制复杂度。采用了简约空间内点法对该模型进行求解。测试结果表明,该模型简洁可靠,所提算法计算效率高,占用内存少,适于求解大规模的TSCOPF问题。

1 极坐标形式的TSCOPF模型

OPF模型在数学上可归纳为如式(1)所示。

其中:x∈Rn为变量;f(x):Rn→R为目标函数;h(x):Rn→Rm为等式约束;g(x):Rn→Rc为不等式约束;分别为不等式约束的上下限。TSCOPF模型的差异主要体现在电磁功率方程及潮流方程的坐标形式的选取上。

1.1 三种坐标形式比较

以有功功率方程为例,目前混合极坐标、极坐标及直角坐标形式的潮流方程都已被采用,如式(2)~式(4)所示,但电磁功率方程还只是局限于混合极坐标,如式(5)所示。

其中:i∈SN,SN为所有节点集合;SG为发电机节点的集合;PGi、PDi为节点i处有功源输出功率、有功负荷;Ui、θi为节点i电压幅值、相角;Gij、Bij为系统导纳矩阵的实部、虚部;Yij、αij为节点导纳矩阵的幅值、相角;ei、fi为节点i电压值的实部、虚部;Ptei、δit为发电机i在t时刻的电磁功率、转子角度;Ei为发电机i的暂态电抗后电势;Geij、Beij为仅保留发电机节点的既约导纳阵实部、虚部。

文献[6]采用式(4)与式(5)组合模型一,文献[3]采用式(2)与式(5)组合模型二,文献[4]采用式(3)与式(5)组合模型三。组合模型一、模型三中出现两种坐标系,使得模型复杂,后续的公式推导及编程的工作量成倍数增加,且出错率增大。比较式(2)、式(3)与式(4)可以发现,混合极坐标形式的复杂度约为极坐标形式的两倍,其计算量远大于极坐标形式。虽然直角坐标模型的海森矩阵为常数阵,但不如极坐标模型简洁。基于此,本文提出基于完全极坐标形式的TSCOPF模型,即潮流方程及电磁功率方程均采用极坐标。极坐标形式的电磁功率方程如式(6)所示。

其中:Yeij、αeij为仅保留发电机节点的既约导纳矩阵的幅值、相角;其他变量定义与前相同。

1.2 极坐标形式的数学模型

本文考虑经典模型:发电机用暂态电抗Xd′和暂态电抗后电势E′相串联表示,且设E′恒定;发电机输入的机械功率Pm恒定;负荷用恒定阻抗表示。取系统网损最小为目标函数。极坐标形式的TSCOPF模型描述如下:

其中:,ST为积分时段的集合;SR为无功源点的集合;QRi、QDi为节点i处无功源输出功率、无功负荷;为发电机i在t时刻的角速度;Δt为积分步长;iT为发电机i的转动惯量;Xd′i为发电机i的暂态电抗;分别为有功源出力上下限、无功源出力上下限和节点电压幅值上下限;为发电机转子相对摇摆角度上下限,取为±100°;δ0COI为稳态时惯性中心角度;δtCOI为预想故障下各时段的惯性中心角度;ωN为角速度基准值;其他变量定义与前相同。

比较各模型可见,本文模型最为简洁,便于记忆及编程。统一的坐标系,使后续的梯度矩阵、海森矩阵的公式推导及编程实现的工作量成倍数减少,有效地降低了出错率,提高了计算效率。

2 简约空间内点法求解步骤

本节首先介绍简约空间技术的推导过程并分析其计算复杂性,然后阐述简约空间内点法解TSCOPF问题的具体步骤。

2.1 简约空间技术公式推导

线性系统问题的标准格式为

其中:D∈Rn×n;A∈Rm×n;Lx∈Rn;Ly∈Rm。

将原变量x分解为两个子空间DY和DZ,令

由式(17)即可得Δy。

2.2 计算复杂性分析

内点法既约修正方程即为如式(7)所示的线性系统。用符号W(·)表示式(7)中系数矩阵,则W(·)∈R(n+m)×(n+m)。W(·)为高度稀疏对称阵,常规的处理方法为采用LDLT分解直接求解。该部分的计算量消耗了整个程序的绝大部分计算时间,其阶数直接影响到内点法的计算效率。

简约空间技术的主要计算量集中于求解式(10)、式(14)、式(15)和式(17),其计算复杂度主要取决于问题自由度的大小。当系统自由度很小,即n-m/n<<1时,系数矩阵C∈Rm×m,B∈R(n-m)×(n-m),修正方程的阶数被大大缩小,能有效地减少解修正方程所耗用的时间。且观察式(10)、式(15)、式(17)的系数矩阵发现,对该三式只须进行一次系数矩阵的分解,可采用前代、回代对其他两式进行求解。

2.3 简约空间内点法求解步骤

原始-对偶内点算法[18]在电力系统中得到了很好的应用。文献[3]亦证明了原始-对偶内点法求解TSCOPF问题的有效性。基于简约空间内点法的TSCOPF问题的求解步骤如下。

步骤1:读入系统数据及故障信息,置kmax=50,中心参数取σ=0.1,计算精度为ε=10-5。

步骤2:常规潮流计算,得到电压初值。

步骤3:根据常规潮流结果计算各负荷等值阻抗,结合预想故障信息,消去网络节点,形成故障中及故障切除后仅含发电机节点的既约导纳阵。

步骤4:设定k=0。选择初值l>0,u>0,z>0,w<0,y=0。

步骤5:判断k≤kmax是否成立,如果不成立,转步骤13。

步骤6:计算互补间隙,如果CGap<ε,则输出最优结果,停机。

步骤7:计算扰动因子

步骤8:计算KKT条件中的Ly、Lw、Lz、Ll、Lu,目标函数、等式及不等式约束的雅克比矩阵和海森矩阵∇f(x)、∇h(x)、∇g(x)、∇2f(x)、∇2h(x)、确良∇2g(x),构造既约修正方程中的H′、Lx。既约修正方程为

其中:

步骤9:简约空间技术解修正方程。

步骤10:计算

步骤11:求最大原始和对偶步长。

步骤12:更新变量

k=k+1,转步骤5。

步骤13:提示“计算不收敛!”,停机。

3 测试结果与讨论

本节以表1所示的五个测试系统为例,验证上述简约空间内点法解极坐标形式的TSCOPF问题的有效性。从表1可以看出,各系统TSCOPF问题的自由度都很小。各系统预想故障信息如表2所示。目前只考虑单重故障。故障节点的故障类型为三相短路。采用32位Matlab R2010b。计算环境为Dual-Core AMD 2.2 GHz,2 GB内存。

互补间隙是衡量系统最优性的重要指标。各系统的互补间隙随迭代次数的变化曲线如图1所示。从图1可以看出,各测试系统的对偶间隙单调收敛至零,该算法具有良好的收敛性。

以CEPRI-22、NE-39和IEEE-118系统为例,发电机转子相对摇摆曲线如图2所示。所有发电机的转子相对角度均在-100°~100°范围内,满足了在仿真时间内各发电机不失稳的要求。

常规原始-对偶内点法与简约空间内点法的测试结果比较如表3所示。从表3可以看出,简约空间内点法的CPU总计算时间及解修正方程所耗用的时间均远小于常规原始-对偶内点法。

若定义简约空间内点法的CPU计算时间与常规原始-对偶内点法的CPU计算时间的比值为加速比,则测试系统的加速比如图3所示。可见,加速比曲线随着系统规模的增大呈单调递增趋势,系统规模越大,简约空间内点法加速越明显。

图4给出了常规原始-对偶内点法和简约空间内点法解修正方程的时间占CPU总计算时间的百分比。两条百分比曲线均随着系统规模的增大呈单调递增趋势,可见,系统规模越大,内点法中修正方程的求解问题越严重。由简约空间内点法的百分比曲线低于常规原始-对偶内点法的百分比曲线可以看出,简约空间内点法有效地缩减了解修正方程所耗用的时间,提高了计算效率。

4 结论

本文提出了极坐标形式的TSCOPF模型,该模型将潮流方程及转子运动方程统一到极坐标形式下,比其他形式的模型更简洁明了,便于记忆及编程。采用了简约空间内点法进行求解,这一方法通过缩小修正方程的规模,减少解修正方程所耗用的时间,达到了提高计算效率、降低内存损耗的目的。

暂态稳定约束 篇6

大停电发生后的系统恢复过程通常分为3个阶段, 即黑启动、网架重构及负荷恢复[1,2]。负荷恢复作为系统恢复的最后阶段, 其目标是在较短的时间内尽可能多地恢复负荷[1]。负荷完全恢复一般需要经历几个小时[3], 常通过分批多次投入来完成。单次负荷恢复可能是1条甚至多条线路的同时投入, 其时间窗口应为秒级, 加之恢复初期系统规模较小和大停电后的冷负荷启动特性, 若单次投入的负荷量过大, 可能引起电压、频率等问题;如果单次投入的负荷量过小, 则会增加操作次数, 延缓恢复进程。因此, 确定合适的负荷恢复步长对加快系统恢复进程具有重要意义[4]。

文献[5]在扩展潮流计算的基础上, 求出运行点的灵敏度, 仅考虑频率约束得出负荷恢复量, 通过最优潮流和松弛负荷量的方法优化该值。文献[6]考虑了冷负荷恢复的典型特性, 建立了调速器及出力爬坡速度的模型, 应用粒子群算法来寻找最优值。文献[7]研究了负荷恢复和系统并网问题, 考虑了网络N-1潮流不可行和电压、热稳定的约束限制, 且频率约束通过动态仿真来实现。文献[8]根据各种原动机的简化方框图, 近似求得其频率的响应率, 由此得出允许频率限制下的最大负荷恢复量。

目前的负荷恢复研究较为详细地考虑了频率的限制, 关于电压问题的考虑主要从稳态网络潮流的约束出发, 而对暂态电压变化过程分析较少。系统恢复初期, 如果单次负荷恢复量过大或功率传输的电气距离过长, 负荷投入时的初始冲击可能会导致暂态电压过低, 引起电动机堵转, 甚至电压失稳。

本文结合实际恢复方案制定中遇到的情况, 在确定变电站最大负荷恢复量时引入了暂态电压安全的约束, 构建了分析暂态电压安全的微分代数方程组, 将描述暂态电压安全问题的二元表作为暂态电压约束满足情况的判据。提出用改进的二分法对模型进行求解, 加快了搜索寻优时的收敛速度。

1 暂态电压安全约束

暂态电压安全包括暂态电压稳定 (TVS) 和暂态电压可接受性 (TVDA) 两方面。本文主要基于TVDA来考虑暂态电压安全问题。TVDA是从电力用户的角度对暂态电压安全性概念的补充, 即故障后系统能够保证负荷节点电压低于某个给定值的持续时间不超过预定时段, 否则认为电压是暂态不安全的[9]。

由于系统网架的未完全恢复, 负荷投入时的初始冲击可能会引起严重的暂态电压下降;另外, 负荷中电动机的启动影响较大, 其暂态电压响应速度一般要快于频率的响应过程。上述两方面原因造成在某些情况下, 暂态电压下降严重而频率仍满足要求。因此, 有必要在确定变电站单次最大负荷恢复量时引入暂态电压安全的约束。

对电力系统的暂态电压安全性进行分析时, 需要一种能够定量评估暂态电压安全与否的实用化判据。文献[10]对不同的暂态电压降落标准进行调查, 列出了各种标准下的暂态低电压水平和所能持续的时间限值 (见附录A表A1) 。文献[9]提出用一组包含电压跌落门槛值Vcr和可接受最大持续时间Tcr的二元表 (Vcr, Tcr) 来描述对TVDA的要求。文献[11]在该二元表的基础上, 给出了山东电网对扰动后各节点TVDA的要求, 并确定了低压切负荷装置的低电压定值 (标幺值) 和延时定值 (0.75, 0.2 s) 。

本文基于暂态电压降落标准和二元表的定义, 保守地确定负荷恢复过程中描述暂态电压安全的二元表为 (0.8, 0.5 s) 。

2 最大负荷恢复模型

系统恢复过程中, 根据电网的运行状态确定待恢复变电站的单次最大负荷恢复量时, 需要在考虑频率和稳态电压约束的基础上, 加入反映暂态电压安全的动态变化约束。

2.1 电压约束

电压约束包括稳态电压约束和暂态电压约束两方面。稳态电压约束通过求解给定负荷下的潮流断面即可检验。暂态电压约束需要在考虑适当的发电机、负荷和网络模型的基础上, 通过求解微分—代数方程组来分析考察。本文在分析暂态电压变化过程时, 对系统中各部分模型进行了适当简化, 在保证精度的同时提高了计算速度。

负荷模型为恒阻抗和电动机相结合的动态模型。其中电动机采用了考虑转子绕组暂态的机电暂态模型, 该3阶模型相对于机械暂态模型在暂态过程中具有更好的仿真精度, 其数学表达式为:

$\{\begin{array}{l}\dot{U}={\dot{E}}^{\prime }+(r{}_{\text{s}}+\text{j}{X}^{\prime })\dot{{\rm I}}\\ \frac{\text{d}{\dot{E}}^{\prime }}{\text{d}t}=-\text{j}s{\dot{E}}^{\prime }-\frac{\dot{E}-\text{j}(X-X^{\prime })\dot{{\rm I}}}{{\rm T}{}_{d0}´}\\ {\rm T}{}_{\text{j}}\frac{\text{d}s}{\text{d}t}={\rm T}{}_{\text{m}}-{\rm T}{}_{\text{e}}\end{array}(1)$

式中:$\dot{U}$为电动机端电压;${\dot{E}}^{\prime }$为转子侧暂态电势;rs为定子电阻;X′为暂态电抗;X为同步电抗;Td0′为定子开路暂态时间常数;Tj为惯性时间常数;Tm和Te分别为机械负载转矩和电磁转矩。

为了考虑电压的动态变化过程, 应加入励磁控制系统的影响, 此时发电机模型宜采用忽略阻尼绕组的Eq′变化模型[12]。

网络方程式用节点电压法可表示为:

$\mathit{Y}\mathit{U}=\mathit{{\rm I}}(2)$

式中:Y为导纳矩阵;U为节点电压向量;I为节点注入电流向量。

分析暂态电压变化时, 发电机和负荷通过修改导纳矩阵中相应节点的导纳值和注入电流来参与网络运算。考虑发电机凸极效应时, 式 (2) 中对应发电机节点i的方程须写成实部与虚部分开的形式:

$\begin{array}{l}\left[\begin{array}{c}{\rm I}{}_{xi}\\ {\rm I}{}_{yi}\end{array}\right]=\left[\begin{array}{cc}G{}_{ii}& -B{}_{ii}\\ B{}_{ii}& G{}_{ii}\end{array}\right]\left[\begin{array}{c}U{}_{xi}\\ U{}_{yi}\end{array}\right]+\\ {\displaystyle \sum _{\begin{array}{l}j=1\\ j\ne i\end{array}}^{{\rm N}}\left(\left[\begin{array}{cc}G{}_{ij}& -B{}_{ij}\\ B{}_{ij}& G{}_{ij}\end{array}\right]\left[\begin{array}{c}U{}_{xj}\\ U{}_{yj}\end{array}\right]\right)}(3)\end{array}$

式中:Gii和Bii分别为导纳矩阵中对角元素Yii的实部和虚部;Gij和Bij (j≠i) 分别为非对角元素Yij的实部和虚部;N为网络节点数目。

动态仿真计算时, 网络方程中发电机节点的相应部分是变化的。当发电机采用Eq′变化的模型时, 设其接在网络中的节点i, 则注入电流的表达式为:

$[\begin{array}{l}{\rm I}{}_{xi}\\ {\rm I}{}_{yi}\end{array}〗=[\begin{array}{l}b{}_{xi}\\ g{}_{yi}\\ \end{array}〗E{}_{qi}´-\left[\begin{array}{cc}G{}_{xi}& B{}_{xi}\\ B{}_{yi}& G{}_{yi}\end{array}\right][\begin{array}{l}U{}_{xi}\\ U{}_{yi}\end{array}〗(4)$

式中:bxi, gyi, Gxi, Bxi, Gyi, Byi为相应系数, 与发电机自身参数和其相对于系统的转子角度有关。

电动机启动时, 暂态电压响应速度一般快于频率的响应速度, 因此在建立网络模型时, 可以忽略发电机端功角的变化, 使得式 (4) 中与功角相关的各系数变为恒定量。由此, 修正后的系数导纳矩阵可简化为恒定矩阵, 迭代求解时, 可以一次求出矩阵的因子表, 在迭代过程中保持恒定不变, 提高了计算速度。通过数值方法计算上述模型即可求得电压的暂态变化曲线, 根据确定的二元表判断是否满足暂态电压安全的要求。

2.2 频率约束

负荷恢复时, 影响频率偏差大小的主要因素为功率缺额占整个系统容量的比例以及备用在不同类型原动机中的分布情况[1]。不同类型的原动机对于频率的响应率也大不相同, 本文在文献[1,8]的基础上, 结合实际的工程原则[13], 给出了一些典型原动机的频率响应率 (见附录B表B1) 。假设恢复初期机组备用量相对于待恢复负荷量足够大, 由此估算出给定负荷恢复量下的频率偏移为:

$\text{Δ}f=\frac{L{}_{\text{f}\text{i}\text{x}}}{\frac{S{}_{\text{S}\text{E}}}{r{}_{\text{S}\text{E}}}+\frac{S{}_{\text{C}\text{Τ}}}{r{}_{\text{C}\text{Τ}}}+\frac{S{}_{\text{Η}\text{E}}}{r{}_{\text{Η}\text{E}}}}(5)$

式中:Δf为系统频率偏差;Lfix为负荷恢复量;S和r分别为已恢复机组的容量及其频率响应率;下标表示原动机类型。

原动机的频率响应率还与初始负荷量的大小有关, 初始负荷越小, 频率偏移越大。附录B表B1的典型值均为最小负荷量下各原动机的频率响应率, 这是一种保守的估计。

2.3 数学模型

确定变电站单次最大负荷恢复量时, 考虑暂态电压、频率和稳态电压约束后, 其数学模型为:

$\{\begin{array}{l}\max L\\ \text{s}.\text{t}.\dot{\mathit{x}}=\mathit{f}(\mathit{x},\mathit{y},L)\\ \mathit{g}(\mathit{x},\mathit{y})=0\\ {\rm T}(V{}_i\le V{}_{\text{s}\text{e}\text{t}})\le {\rm T}{}_{\text{s}\text{e}\text{t}}(i=1,2,\cdots ,{\rm N})\\ f{}_{\min }\le f(L)\le f{}_{\max }\\ V{}_{\min }\le V{}_{i,\text{s}\text{t}}\le V{}_{\max }(i=1,2,\cdots ,{\rm N})\end{array}(6)$

式中:L为目标变电站的负荷恢复量;Vi, st为稳态电压, 通过潮流计算求出;Vmax和Vmin分别为稳态电压的上、下限。

微分方程表示系统元件的动态过程, 代数方程主要包括网络方程和发电机、电动机定子端的电压方程;求解微分—代数方程组可以分析系统暂态过程中节点电压Vi低于电压限定值Vset的时间T, 通过T≤Tset来反映暂态电压约束;通过式 (5) 估算出给定负荷恢复量的频率偏移, 根据实际工程原则中的上下限来约束频率。

3 模型求解

3.1 改进的二分法

由式 (6) 所示模型可以看出, 求解变电站单次最大负荷恢复量为一维非线性多约束问题。该问题可以通过二分法进行求解, 即根据当前约束条件的满足情况更新负荷量取值的上下限, 取该范围的中点为新的负荷量继续迭代, 直至满足收敛条件。

为了加快收敛速度, 本文提出一种改进的二分法对模型进行求解。该方法基于传统二分法的搜索思路, 根据每次迭代时负荷量取值的上下限及其相应的约束裕度, 通过线性拟合求出新的负荷量位置, 即寻优过程中考虑了每次迭代的约束裕度。通过该方法确定新的负荷量时, 其表达式为:

$L{}_{\text{n}\text{e}\text{w}}=\frac{|\eta {}_{\max }|L{}_{\min }+|\eta {}_{\min }|L{}_{\max }}{|\eta {}_{\max }|+|\eta {}_{\min }|}(7)$

式中:Lnew为新的负荷量;Lmax和Lmin分别为当前迭代过程中更新后的负荷量取值上下限;η为相应负荷量的约束裕度, 根据满足情况可正可负。

式 (7) 中的约束裕度η由暂态电压约束、频率或稳态电压约束来确定, 其中频率和稳态电压约束为1维不等式约束问题, 其裕度容易获得且与相关参数对应的曲线具有较好的光滑性。暂态电压约束为2维不等式问题, 其裕度—参数曲线具有高度的非线性和不光滑性, 以此修正负荷量时可能造成迭代次数的增加。为了改善这种情况, 本文采用文献[9]中的方法, 通过曲线拟合技术求取暂态电压约束裕度, 将2维不等式约束转化为1维的相应问题。

3.2 计算流程

采用改进的二分法对式 (6) 所示模型进行求解时, 其计算流程见附录C图C1。每一次迭代过程中, 根据当前负荷量依次考察暂态电压、频率和稳态电压约束的满足情况, 计算相应的约束裕度。根据约束裕度的正负更新负荷量的上下限, 然后通过式 (7) 确定新的负荷量进行下一次迭代, 直到满足收敛条件, 即所有约束均未越限, 并且存在约束裕度的绝对值小于给定正值ε的情况, 此时的负荷量即为最优负荷量。

4 算例分析

4.1 IEEE14节点算例

以IEEE 14节点测试系统为例, 其结构图和参数见附录D。假定已恢复的系统负荷为150 MW, 对于不同的负荷模型, 节点14单次最大负荷恢复量的计算结果如表1所示, 与其相对应的暂态电压变化曲线如图1所示。

表1和图1的结果表明, 系统在该运行状态下对节点14进行负荷恢复时, 频率和稳态电压并不是主要的约束因素 (频率偏差要求在±0.5 Hz内, 电压在0.9～1.1范围内) 。最大负荷恢复量主要受暂态电压约束的影响, 且其大小与负荷中电动机的比例成反比。

假定负荷模型采用表1中的模型1, 则对于不同的系统负荷水平, 节点14的最大负荷恢复量见附录E表E1。其结果同样表明, 3种负荷水平下节点14的最大负荷恢复量由暂态电压约束条件决定。且随系统负荷水平的增高, 最大负荷恢复量有所减小, 说明负荷水平的增加导致暂态电压下降幅度和持续时间都有所增加。

4.2 实际系统

以山东电网为例, 发生大规模停电事故后, 泰安抽水蓄能电厂 (泰抽电厂) 具有黑启动能力, 其1号机组启动后经泰山、天平、桃园、高余、石横乙站送石横乙5号机组厂用电。石横乙5号机组并网后, 与泰抽电厂1号机组形成稳定小系统, 需要恢复该地区部分重要负荷, 如图2所示。

泰抽电厂1号机组容量为295 MVA, 假定其运行在最小稳定功率附近。石横乙5号机组容量为350 MVA。假定系统当前负荷水平为175 MW, 对于不同的负荷模型, 南郊站能够恢复的最大负荷量如表2所示。

假定采用表2的负荷模型A, 对于不同的系统负荷水平, 南郊站最大负荷恢复量见附录E表E2。

表2和附录E表E2的结果同样表明在上述恢复状态下, 暂态电压约束是寻优过程中的主要制约因素, 且最大负荷恢复量与负荷模型以及系统负荷水平有很大关系。

根据表2所给的3种负荷模型, 求取最大负荷量时分别采用改进二分法和传统二分法, 其对比情况如表3所示。结果表明, 2种方法下得到的最大负荷恢复量差别不大, 但改进二分法具有更快的收敛速度。该方法在Borland C++Builder软件开发环境下实现, 当计算机硬件条件为Pentium Dual E2140 1.6 GHz、内存1 GB时, 表3中通过改进二分法得到的结果均在1.5 s～2.0 s内获得。

5 结语

负荷恢复的目标是在满足安全约束的条件下尽快地恢复负荷, 以减少社会影响和停电损失。本文引入暂态电压约束, 构建了变电站最大负荷恢复量模型, 并对系统各部分模型进行了适当简化, 提出用改进的二分法进行求解。IEEE 14节点算例和山东电网的仿真结果表明:

1) 在某些情况下, 负荷恢复会造成严重的暂态电压跌落, 暂态电压约束成为制约变电站单次最大负荷恢复量的主要因素。

2) 负荷构成和系统负荷水平等因素对负荷恢复量有重要的影响。

3) 所建立的数学模型能够有效地确定变电站单次最大负荷恢复量, 且模型的适当简化和改进二分法的使用加快了计算速度。

摘要：电力系统恢复过程中, 需要确定变电站的单次最大负荷恢复量。通常主要考虑系统频率和稳态电压的约束, 文中在此基础上加入了暂态电压安全的约束限制, 建立了确定变电站最大负荷恢复量的数学模型。基于描述暂态电压安全问题的二元表, 构建了分析暂态电压安全的微分代数方程组, 对系统中各部分组成进行适当简化, 在保证精度的前提下提高了计算速度。提出一种改进的二分法对模型进行求解, 该方法考虑了每次迭代的约束裕度, 加快了搜索寻优时的收敛速度。IEEE14节点算例和山东电网仿真结果表明了模型的有效性和方法的快速性。

暂态稳定约束 篇7

在半导体制造中, 由于晶圆直径越来越大, 加工工艺复杂, 组合设备被广泛的用来加工晶圆。组合设备是一类集成制造设备。在这样的一台设备中, 有几个加工模块 (Processing Module, PM) , 一个机械手以及真空锁 (Loadlock) 。机械手分为单臂机械手和双臂机械手, 如图1所示为单臂机械手, 因此叫单臂组合设备。此系统的工作过程大致如下:一批装在一个卡盒中待加工的晶圆首先放入真空锁中, 然后由机械手从卡盒中抓取一枚晶圆放到AL中, 经AL的视觉系统进行校准后, 由机械手按照工艺要求将该枚晶圆依次放入相应的PM中加工, 晶圆在每一个加工模块内完成一道工序, 当完成最后一道工序后, 晶圆被送入CL中冷却, 避免热的晶圆影响到真空锁内其他的晶圆, 最后晶圆被送回到真空锁。整个过程由程序自动控制, 当该批晶圆加工完成后, 真空锁中卸载, 并装入下一批晶圆。两个真空锁通过轮流使用实现对晶圆的装卸作业, 以保证组合装置加工的连续性, 使得系统绝大部分时间都是处在稳定运行的状态。这样的组合设备可以为半导体制造提供一个灵活的、可重置的和有效的环境, 减小生产节拍, 提高空间的利用率, 降低成本。一个晶圆通过真空所进入系统, 然后按照已知的加工路径在加工模块中加工, 当所有的加工工序完成后, 由机械手送回真空所。

为了调度组合设备, 许多学者在建模和对系统执行评估做了很多工作[1,2]。基于这些研究发现, 在稳态下, 当机械手有空闲时间时, 加工模块的加工时间决定系统的生产周期, 当机械手一直忙的时候, 系统的生产节拍由机械手的活动时间决定。在一台组合设备中, 仅仅当一个晶圆载入到加工模块中, 加工模块才开始加工晶圆。因此, 加工模块的活动跟随着机械手的任务[3]。这样, 调度组合设备的关键在于调度机械手的任务。众所周知, 组合设备中机械手从一个加工模块移动到另一个加工模块的时间可以看作为常数, 并且比晶圆的加工时间短很多[4]。这样, 对于单臂组合设备拉式调度是最优的[5,6]。然而, 以上的研究结果都是基于假设晶圆加工完成后可以在加工模块中停留无限长的时间。

对于一些晶圆加工工艺, 例如低压化学沉积工艺, 它对晶圆在加工模块中的滞留时间有严格的约束。这种约束叫晶圆驻留时间约束 (residency time constraint, RTC) [4]。RTC指一个晶圆在加工模块中加工完成后, 必须在规定的有限的时间内卸载出来, 否则此晶圆将成为废品。由于在加工模块之间没有缓冲区, 使得很难调度组合设备去满足晶圆在每个加工模块中的驻留时间约束。这样, 为了能够调度带驻留时间约束的组合设备, 文献[4]提出了一些方法为双臂组合设备找到最优的周期性调度。之后, 为了提高计算有效性, 伍乃骐等人分别对于单臂和双臂组合设备提出了充分必要的可行性调度条件[2,7]。如果可调度, 可通过由解析表达式组成的调度算法找到最优的周期性调度。然而在以上的研究中, 并没有给出怎么样进入稳态调度更好。也就是说, 不同文献给出进入稳态的调度方法不一样。这就意味着, 不同进入稳态的调度方法导致不同的稳态调度方案。这样, 本文基于earliest调度策略, 通过对系统暂态过程的分析, 给出了系统进入稳态后, 机械手的等待时间分布。

接下来文章主要介绍Petri网建模。基于此Petri网模型, 第二部分分析了系统进入稳态后机械手的等待时间分布。第三部分给出了实例验证了本文得出的结论。第四部分对文章进行了总结。

1 Petri网建模

本文用Petri网 (Petri net, PN) 对系统进行建模。Petri网的概念来自于文献[8-9], 对于没有重入加工工艺的组合设备, 晶圆的加工模式可以定义为 (m1, m2, …, mn) , 其中n是晶圆加工步骤数目, mi是晶圆加工步骤i中所包含的加工模块数目。

本文展示的Petri网模型中, 用库所pi表示步骤i的加工模块并有K (pi) =mi, iϵ{1, 2, …, n}。将真空锁看作一个加工步骤并称之为步骤0, 用p0表示。由于真空锁可以容纳设备中的所有晶圆, 因而有K (p0) =∞。这样, 整个系统就有n+1步加工步骤。让Ω={0, 1, …, n}表示包含步骤0的步骤集, 让Nn={1, …, n}表示可加工晶圆的步骤集。机械手由库所r表示, K (r) =1意味着机械手只有一只手臂。当库所r中有令牌时表示机械手是空闲状态。

对于带驻留时间约束的单臂组合设备, 机械手的等待时间用来平衡加工步骤之间的工作负载。用库所qi2和qi1表示机械手在缷载或安装晶圆到加工步骤i之前的等待时间。模型中变迁, si1, iϵNn, 和s01分别表示机械手安装一个晶圆到加工模块pi和真空锁p0中。变迁si2, iϵNn-1, 表示机械手从加工模块pi卸载一个晶圆然后移动到加工模块pi+1。变迁s02表示机械手从真空锁p0下载一个晶圆然后移动到加工模块p1。变迁sn2表示机械手从加工模块pn中下载一个晶圆然后移动到真空锁。yi表示机械手从第 (i+2) 道工序旋转至第i道工序, yn-1表示机械手从真空所旋转至第 (n-1) 道工序, yn表示机械手从第1道工序旋转至第n道工序, y0表示机械手从第2道工序旋转至第0道工序。图2展示了单臂组合设备具有n个加工步骤的Petri网模型。

假设组合设备按照给定的晶圆加工流模式 (m1, m2, …, mn) 对晶圆进行加工, 稳态时系统达到最大产能时有∑i0=1mi个晶圆同时在系统中加工。所以, 让初始状态M0为M0 (pi) =mi, iϵNn, M0 (r) =1来表示机械手为空, 让M0 (p0) =n表示系统中总有晶圆可被载入到系统中。为了避免图2中的Petri网模型死锁, 给出让Petri网避免死锁的控制策略。

定义1[2]:在标识M下, 如果M (pi+1) =mi+1-1, 则变迁yi, iϵNn-1∪{0}, 是可触发的;如果M (pi) =mi, iϵNn, 则变迁yn是可触发的。

通过定义1, 强迫Petri基于拉式策略运行。下面我们展示采用拉式策略的单臂组合设备的运行过程。稳态时一定有状态为M (pi) =mi, iϵNn, M (r) =1, 接着系统的运行过程如:, 这样系统走完一个周期。基于拉式策略, 系统不断的周期性重复以上的运行过程。

为了解决该系统的调度问题, 需将系统的活动时间与系统的Petri网模型联系起来。令机械手每次的运动时间为μ, 机械手从ρ0卸载一枚晶圆的时间为λ0, 而在其他步骤中卸载或安装晶圆的时间均为λ。这样, 有激活变迁yi, iϵΩ, 需要μ个单位时间, 激活变迁si1, iϵΩ, 需要λ个单位时间, 激活si2, iϵNn, 需要λ+μ个单位时间, 激活s02需要λ0+μ个单位时间。用ai表示晶圆在加工模块中的加工时间。由于考虑驻留时间约束, 这样用δi表示这样的约束, 即:晶圆在第i步加工完后, 需要在δi个单位时间内移走, 否则将变为废品。用ωi1和ωi2分别表示机械手在库所qi1和qi2中的等待时间。这样就完成了对单臂组合设备的Petri网建模。

2 基于Earliest调度策略, 机械手等待时间的分布规律

Earliest调度策略是指: (1) 机械手到达加工模块卸载晶圆时, 如果加工模块中的晶圆已经加工好, 那么机械手立即将其卸载下来; (2) 机械手到达加工模块卸载晶圆时, 如果加工模块中的晶圆还没加工好, 那么机械手在此等待, 一旦晶圆加工好, 机械手立即将其卸载下来; (3) 机械手拿着待加工的晶圆来到加工模块时, 如果该加工模块是空的, 立即将其载入加工模块中。暂态是指组合设备所有加工模块中从没有加工晶圆开始, 机械手将晶圆一个个的载入到组合设备中, 到所有的加工模块都有加工晶圆并最终组合设备进入稳态运行的过程。

根据文献[7], 设τi表示第i道工序的晶圆在加工模块中的停留时间, 如果晶圆在第i道工序满足驻留时间约束, 那么τiϵ[ai, ai, +δi]。根据图2中的Petri网模型, 第i步完成一个晶圆需要经过以下过程:触发si2 (λ+μ) →触发si+1 (λ) →触发yi-1 (μ) →机械手在第i-1步等待 (ω (i-1) 2) →触发s (i-1) 2 (λ+μ) →触发sil (λ) →晶圆在第i步加工 (τi) →触发si2 (λ+μ) 。则该过程所需时间为ρi=τi+4λ+3μ+ω (i-1) 2。这样, 考虑并行模块, 第i道工序完成一个晶圆允许的最短的时间为:

第i道工序完成一个晶圆允许的最长时间为:

第一道工序完成一个晶圆允许的最短时间为:

第一道工序完成一个晶圆允许的最长时间为:

在图2的Petri网模型中, 分析机械手走过一个周期的活动路径及时间:触发yn (μ) →机械手在qn步等待 (ωn) →触发sn2 (λ+μ) →触发s01 (λ) →触发yn-1 (μ) →机械手在qn-1步等待 (ωn-1) →触发s (n-1) 2 (λ+μ) →触发sn1 (λ) →…→触发yi (μ) →机械手在qi步等待ωi→触发Si2 (λ+μ) →触发s (i+1) 1 (λ) →触发yi-1 (μ) →机械手在q (i-1) 2步等待 (ω (i-1) 2) 触发s (i-1) 2 (λ+μ) →触发si1 (λ) →…→触发y1 (μ) →机械手在q12步等待 (ω1) →触发s12 (λ+μ) →触发s21 (λ) →触发y0 (μ) →机械手在q02步等待 (ω02) →触发s02 (λ0+μ) →触发s11 (λ) →再次触发yn (μ) , 这样一个周期完成, 若令机械手的节拍为ψ, 可得到:

式 (5) 中可以发现2 (n+1) μ+ (2n+1) λ+λ0是一个常数, 可令其为ψ1, 而是不确定的, 可令其为ψ2, 所以

根据文献[7], 晶圆在加工模块的停留时间如下式:

令表示第i道工序的生产节拍, 根据文献[7], 当所有的加工模块都有晶圆在加工时, 有下面的结论。

引理1:对于带驻留时间约束的单臂组合设备, 当所有的加工模块都有晶圆在加工时, 有θi=ψ, iϵNn。

让Πi L和Πi U分别表示第i道工序允许生产节拍的上限和下限, 那么有:

文献[7]针对带驻留时间约束的单臂组合设备的调度问题, 提出了使系统可调度的充分必要条件。在此基础上, 给出稳态下机械手的等待时间分布规律。假设第k道工序为瓶颈工序, 1≤k≤n。那么有下面的定理。

定理1:对于带驻留时间约束的单臂组合设备, 基于earliest调度策略, 当系统到达稳态后, 等待时间ωi1=0, iϵΩ和ωi=0, i=0, 1, ……, k-1, 成立。

证明:根据拉式调度策略, 当机械手拿着待加工的晶圆来到第i步时, 第i步一定为空。这样根据earliest调度策略定义, 机械手可以立即将此晶圆载入第i步中, 即ωi1=0。因为真空锁总有待加工的晶圆, 这样当机械手来到真空锁可立即将其中的晶圆卸载出一个, 所以有ω02=0。根据拉式调度策略, 当机械手来到第一步时, 其一定已经在瓶颈工序第k步完成一系列的机械手活动, 这样此时该系统的周期由瓶颈工序决定, 即为Πk L。这样由引理1和式子 (8) 可得τ1=m1×Πk L- (3λ+λ0+3μ+ω0) =m1×Πk L- (3λ+λ0+3μ) ≥m1×Π1L- (3λ+λ0+3μ) =a1, 所以当机械手来到第1步卸载晶圆时, 第1步中的晶圆已经加工好了。因此机械手可以立刻将其中的晶圆卸载出来, 同时不需要等待, 即ω1=0。当机械手来到第2步时, 同理由引理1和式子 (7) τ2=m2×Πk L- (4λ+3μ+ω1) =m2×Πk L- (4λ+3μ≥m2×Π2L- (4λ+3μ) =a2, 所以当机械手来到第2步卸载晶圆时, 第2步中的晶圆已经加工好了。因此机械手可以立刻将其中的晶圆卸载出来, 同时不需要等待, 即ω2=0。相似的, 可推出ω3=ω4=…=ωk-1=0。因此定理1成立。

根据文献[10], 可知机械手在qi2中的等待时间是有上界的, 否则此等待时间将会影响系统的最优性。让γi表示可分配到库所qi2最大的等待时间, 这样有γi-1=mi×Πk L- (ai+4λ+3μ) , 1≤i≤n和ωnϵ[0, ∞]成立[10]。

证明:根据定理1可知, ωi1=0, iϵΩ和ωi2=0, i=0, 1, …, k-1, 成立。假设机械手在第k步的等待时间ωk2>γk2, 根据引理1和式子 (7) 有τk+1=mk+1×Πk L- (4λ+3μ+ωk) ak+1, 根据earliest调度策略, 这就暗示着不断减小同时ωk2不断增加并且。直到当ωk2=γk2时, 有τk+1=mk+1×Πk L- (4λ+3μ+ωk2) =ak+1, 这时根据earliest调度策略, ωk2此时不再增加。同理可以得到, 如果, 而机械手来到第d步时, 由于ψ=Πk L, 这样, 而且此时有成立, 这就意味着ωi2=0, d+1in, 成立。如果, 那么ωi2=γi2, i

到现在为止, 分析了基于earliest调度策略, 稳态时机械手等待时间的分布情况。同时, 机械手的等待时间分布可通过解析表达式计算出来, 因此非常有效。

3 实例验证

一台组合设备有四道加工工序, 其中第三道工序有两个并行加工模块, 而其他工序只有一个加工模块。每道工序的加工时间分别为:a1=92 s, a2=98 s, a3=210 s, a4=95 s。机械手的活动时间分别为:λ=10s, λ0=15s, μ=2s。

基于工艺参数可以求得:ψ1=125s, Π1L=143s, Π2L=144s, Π3L=128s和Π4L=141s, 其中第二道工序为瓶颈工序。基于定理1和定理2得出, 稳态时, 机械手的等待时间分布为:ω02=ω12=0, ω22=19和ω32=ω42=0。通过基于e M-plant仿真, 采用earliest调度策略, 系统稳态时机械手等待时间的分布为ω02=ω12=0, ω22=19和ω32=ω42=0。此仿真结果与基于定理1和定理2得出的结果一致, 因此验证了本文的结论。

4 结束语

半导体加工是极复杂的制造系统, 它的调度问题非常复杂。本文对基于earliest调度策略的带驻留时间约束的单臂组合设备的暂态过程进行了分析, 运用Petri网模型得出了系统达到稳态后机械手的等待时间分布情况, 并且此分布可通过解析表达式计算得出, 因此非常有效。基于以上成果, 今后将着重研究基于earliest调度策略多组合设备到达稳态时机械手的等待时间分布情况。

参考文献

[1]J.Yi, S.Ding, M.Zhang.Steady-state throughput and scheduling analysis of multi-cluster tools:A decomposi tion approach[J].IEEE Trans.Automation Sci.Eng., Apr, 2008, 5 (2) :321-336.

[2]N.Q.Wu, M.C.Zhou.A closed-form solution for schedu lability and optimal scheduling of dual-arm cluster tools with wafer residency time constraint based on steady schedule analysis[J].IEEE Transactions on Automa tion Science and Engineering, 2010, 7 (2) :303-315.

[3]Y.-H.Shin, T.-E.Lee, J.-H.Kim, et al.Modeling and implementing a real-time scheduler for dual-armed cluster tools[J].Computers In Industry, 2001, 45 (1) :13-27.

[4]J.-H.Kim, T.-E.Lee, H.-Y.Lee, et al.Scheduling analysis of timed-constrained dual-armed cluster tools[J].IEEE Transactions on Semiconductor Manufactur ing, 2003, 16 (3) :521-534.

[5]T.-E.Lee, H.-Y.Lee, Y.-H.Shin.Workload balanc ing and scheduling of a single-armed cluster tool[A].Proceedings of the 5th APIEMS Conference[C], Gold Coast, Australia.2004:1-15.

[6]M.-J.Lopez, S.-C.Wood.Systems of multiple cluster tools-configuration, reliability, and performance[J].IEEE Transactions on Semiconductor Manufactur ing, 2003, 16 (2) :170-178.

[7]N.Q.Wu, C.B.Chu, F.Chu, et al.A Petri net method for schedulability and scheduling problems in single-arm cluster tools with wafer residency time constraints[J].IEEE Transactions on Semiconductor Manufacturing, 2008, 21 (2) :224-237.

[8]M.C.Zhou, K.Venkatesh.Modeling, simulation and control of flexible manufacturing systems[A].Petri net approach, World Scientific, Singapore, 1998.

[9]T.Murata.Petri nets:Properties, analysis and applica tions[A].Proceedings of the IEEE[C].1989, 77 (4) :541-580.