基于BPA潮流的電力系統中長期頻率電壓穩定仿真方法

夏冰清，周 華，何哲楠，華 文，吳 浩

（1.浙江大學電氣工程學院，杭州 310027；2.國網浙江省電力公司，杭州 310007；3.國網浙江省電力公司電力科學研究院，杭州 310014）

在電力需求不斷攀升和電力系統規?？焖侔l展的今天，為了保障電力系統的穩定運行，電網應用了大量的自動控制系統和保護裝置，使得系統的動態行為變得非常復雜。由于不同控制系統或裝置間的惡性交互可能危及較長時間后的系統穩定性，故在系統暫態穩定得以緩解的條件下，中長期動態過程對穩定性的影響逐步受到重視[1?3]。頻率穩定和電壓穩定是電力系統穩定的兩個重要方面，隨著大規模直流輸電、高比例可再生能源接入技術的廣泛應用，受端側電網中本地常規機組開機日益減少，使得電網調控能力下降，在直流故障導致的受電有功功率下降和無功功率需求上升的情況下，系統可能發生頻率穩定或電壓穩定問題，甚至兩種穩定問題交織出現[4?5]，因此，在研究實際系統中長期動態過程時，需要同時分析頻率穩定和電壓穩定，也即研究電力系統中長期頻率電壓穩定問題。

考慮到頻率穩定和電壓穩定的失穩機理不同，現有研究大多沿襲電力系統穩定性分類方法將兩種穩定問題獨立研究。文獻[6?7]針對中長期電壓穩定問題建立了研究模型，分析了負荷恢復、發電機過勵磁限制器、有載調壓變壓器分接頭調整等因素對中長期電壓穩定的影響，但忽略了中長期頻率動態的相關模型；文獻[8]研究了交直流混聯系統的中長期頻率穩定問題，分析了機組不等率、調頻死區等因素對系統一次調頻性能的影響，但忽略了中長期電壓變化對系統有功平衡的影響。盡管傳統中通常獨立研究頻率穩定性和電壓穩定性，但頻率和電壓的變化卻存在一定的相互影響，如系統電壓的下降通常會減小負荷的有功需求，一定程度上有利于系統的頻率穩定性；系統頻率的下降卻通常會增大系統的無功需求，一定程度上不利于系統的電壓穩定性；而低壓或低頻減載不僅有針對性地緩解系統低壓或低頻這一側面的威脅，而且會影響到頻率或電壓的另一側面。因此，為了更準確地分析系統的中長期穩定性，有必要綜合考慮影響系統頻率穩定和電壓穩定的因素。

時域仿真法是各類穩定問題的基本的研究工具[9?11]。文獻[12?13]針對電力系統中長期動態過程的分析需要，提出了準穩態QSS（quasi?steady?state）近似的思想，即在暫態穩定的前提下，將暫態過程用準穩態平衡方程代替，從而設計了準穩態仿真方法，簡化電力系統中長期動態仿真，提高仿真效率，獲得了廣泛關注[14?16]。顯然，該思想及其方法同樣適用于中長期頻率電壓穩定問題，能夠快速仿真各類中長期動態設備的動態響應，從而分析不同因素對系統穩定產生的影響。

本文針對電力系統中長期頻率電壓穩定問題，建立了綜合考慮頻率和電壓穩定影響因素的準穩態仿真模型，提出了一種基于既有潮流計算程序的準穩態時域仿真方法，同時，在計算效率方面，提出了變步長仿真方法和基于多時步預測的初始值設定方法。最后基于BPA軟件的潮流計算功能和Mathematica軟件[17]的開發平臺，實現了所提出的中長期頻率電壓準穩態仿真方法，并通過Nordic32算例系統和某大區電網算例系統，研究了所提方法的有效性和對大規模系統的適應性。

1 中長期頻率電壓動態過程的準穩態模型

1.1 考慮頻率電壓動態的準穩態模型

對電力系統動態，若主要關注中長期過程，則在準穩態近似思想下，可對短期暫態過程進行簡化，即認為系統始終處于暫態平衡點上，將描述暫態過程的微分方程替換為代數方程，從而僅保留描述中長期動態的方程，形成電力系統動態的準穩態模型（簡稱QSS模型）。

為了準確反映電力系統的中長期頻率電壓動態過程，在QSS模型中，需要考慮影響中長期過程的各類因素，包括負荷恢復特性、有載調壓變壓器OLTC（on?load tap changer）、發電機過勵限制器OEL（over excitation limiter）、自動發電控制 AGC（auto?maticgenerationcontrol）等。一般地，電力系統中長期頻率電壓動態過程的QSS模型可抽象地寫為

式中：x為暫態過程的狀態變量，代表發電機轉速、勵磁大小、調速器信號量等；y為代數變量，代表節點電壓、相角、注入功率等；zc和zd分別為中長期動態的連續狀態變量和離散狀態變量，代表節點負荷恢復值、發電機節點AGC信號、OLTC變比等；k為離散動態的仿真時刻計數。式（1）和式（2）分別為網絡方程的潮流平衡部分和節點功率頻率電壓特性修正部分；式（3）為系統暫態過程所對應微分方程的平衡形式；式（4）和式（5）分別為描述系統中長期動態過程的連續和離散方程。各式的具體內容分述如下。

1.2 中長期動態過程方程

中長期動態過程方程式（4）和式（5）反映了負荷恢復、OEL、OLTC、AGC等元件的慢動態過程，這些元件是引起中長期動態過程的主要元件，因此在QSS模型中，需保留其動態方程。由于篇幅限制，本文所涉及模型的詳細內容見文獻[18?22]。

1.3 暫態過程的平衡方程

暫態過程的平衡方程式（3）反映了發電機組調速器、自動電壓調節器、勵磁系統等暫態元件的穩態平衡形式。其QSS模型推導的基本思想為：在中長期時間尺度下，這類元件的暫態過程已經結束，從而可用其暫態模型對應的穩態平衡形式描述其穩態特性。顯然，暫態模型所對應的目標穩態不同時，其穩態平衡形式也將不同。經典模型[6，23]下，暫態元件的QSS模型的具體形式如下。

1）發電機組調速器QSS模型

QSS模型中可認為系統各節點頻率始終保持一致，即系統具有統一頻率f，從而系統的頻率偏差Δf為

式中，f0為額定頻率。

發電機組調速器的作用在QSS模型中可表示為頻率的一次調整，則發電機i的有功功率PGi為

式中：PGi0為額定有功功率；KGi為單位調節功率。

2）自動電壓調節器QSS模型

發電機自動電壓調節器AVR（automatic volt?age regulator）的作用在QSS模型中可表示為：在OEL未動作時，發電機i機端電壓保持恒定，即

式中：Ui為發電機機端電壓；Ui0為參考電壓。

3）勵磁系統QSS模型

在QSS模型中，發電機勵磁電壓Efd和功角δ與發電機有功出力PG和無功出力QG的關系[6]可表示為

式中：Ui和θi為發電機i的機端電壓幅值和相角；rai、xdi和xqi分別為發電機i的定子電阻、d軸電抗和q軸電抗。需說明的是，由于系統處于準穩態過程，故可認為發電機功角δ與節點相角θ的差不受頻率f的影響。

4）自動發電控制AGC

自動發電控制作用是彌補調速器一次調頻后的頻率偏差。AGC的典型控制方法一般以區域控制偏差ACE（area control error）信號作為控制標準，根據控制目的不同ACE可主要分為3種模式：恒定頻率控制、恒定交換功率控制以及聯絡線和頻率偏差控制，如恒定頻率控制下的ACE=Δf。

由于在QSS模型中已經忽略調速器的暫態過程，所以認為發電機出力響應可以跟隨AGC控制信號變化，但是其變化率依然受限于調頻機組的爬坡速度（v）。對于每個發電機的從沒有系統偏差的t0時刻至tk時刻，AGC信號表示為

式中：正負號根據頻率的調整方向決定；B為偏差調節系數，一般設為1；PG0為發電機額定功率；v為機組爬坡速度上限，典型值為30 MW/min，Δt為仿真步長。

5）有載調壓變壓器OLTC

有載調壓變壓器一般用于調節負荷側電壓大小，起到穩壓作用，是典型的電壓調節離散型設備。典型模型中需要設置控制節點及其控制范圍以及動作時延，OLTC的動作時延為8～12 s，動作邏輯可表示為

式中：nk表示在k時刻OLTC設置的檔位；Δn為單次調節檔位數，一般為1；U2為控制電壓（變壓器二次側電壓）；為二次側目標電壓；d為控制死區；nmax為單向最大調節檔位數。

1.4 考慮頻率電壓特性的網絡方程

網絡方程的潮流平衡部分式（1）可用節點功率平衡方程描述，即

式中：Pi和Qi分別為節點i的注入有功功率和無功功率；Ui和Uj分別為節點i和j的電壓幅值；θij、Gij和Bij分別為節點i和j間的相角差、電導和電納。同樣地，由于系統具有統一頻率f，故θij不受頻率f的影響。

顯然，考慮系統的一次調頻和AGC特性、發電機的自動電壓調節器和OEL特性以及負荷的頻率電壓靜態特性和恢復特性等后，式（10）的節點注入功率需要修正。方便起見，本文將相關修正方程匯總為式（2），具體分析如下。

（1）對發電機的有功功率PGi，在一次調頻特性式（7）的基礎上，考慮AGC后，可修正為

式中，ΔPGi為AGC指令出力調整值。

（2）對發電機的無功功率QGi，在發電機OEL動作的情況下可取固定最大值Qlim，即

式中：ra、xd、xq、δ、θ、U、Efd分別為發電機的定子電阻、d軸電抗、q軸電抗、功角、機端電壓幅值、機端電壓相角和勵磁電勢，Efd,max一般設為額定值的180%。

在發電機OEL未動作的情況下，可認為AVR將Vi控制為參考值，即式（8），其相應的QGi應保證節點無功功率平衡方程成立。

（3）對應負荷的有功功率PLi和無功功率QLi，在考慮負荷頻率電壓特性后，可描述為

式中：Pd和Qd為節點的功率需求；PL、QL為實際有功功率與無功功率；P0、Q0為額定電壓及頻率下的有功功率與無功功率；αP、βP、γP分別為有功負荷恒功率、恒電流、恒阻抗部分的比例，αP+βP+γP=1，TP和TQ分別為有功功率和無功功率的恢復時間常數，典型值為30s；U為實際電壓標幺值；KP、KQ分別為有功功率與無功功率的頻率因子，典型值為4%/Hz和0.5%/Hz；Δf為負荷母線頻率偏差標幺值。

綜上，考慮發電機和負荷的頻率電壓特性及中長期動態特性后，式（12）中的Pi和Qi分別為

此即節點功率的頻率電壓特性修正方程式（2）。

表1匯總了本文QSS模型涉及的變量，其中，n、m和k分別為節點數、發電機數和負荷數。

表1 QSS模型涉及的變量Tab.1 Variables associated with QSS model

需說明的是，式（1）對應于式（12），其方程數為2n；式（2）對應于式（17），包含發電機和負荷的有功功率和無功功率修正方程，其方程數為2（m+k）；由于AVR等暫態過程平衡方程的影響已在式（2）中考慮，故式（3）僅對應于式（6）和式（9），相關方程數為2m+1。因此，式（1）～式（3）的方程數為2n+4m+2k+1，其數與表1所涉及的變量數相等。

2 基于潮流計算的準穩態仿真方法

2.1 基于交替求解法的準穩態仿真

與暫態仿真類似，準穩態仿真（簡稱QSS仿真）需解取QSS模型式（1）～式（5）的聯立解，其關鍵在于代數方程式（1）～式（3）和微分?差分方程式（4）～式（5）的交接處理[23]。相較于聯立式（1）～式（5）的直接求解方法，交替求解式（1）～式（3）和式（4）～式（5）的方法具有更好的靈活性，可分別為兩者選取合適的算法，且各自修改或增減時，互不影響，更有利于編程實現。因此，本文采用交替求解法求解式（1）～式（5）的聯立解。

為方便敘述，將全部變量(x,y,zc,zd)記為X。由于改進歐拉法較好地兼顧了計算效率和精度，故本文基于該方法實現交替求解，其具體步驟如下：

步驟1設仿真步長Δt和仿真時間T，置t=t0；

步驟2根據式（4）和式（5），計算t+Δt時刻的zc和zd，有，并據此修改QSS模型參數和元件狀態；

步驟3根據式（1）～式（3），由求解代數方程組，得到t+Δt時刻x和y的取值；

步驟4根據式（4）～式（5），更新t+Δt時刻的zc和zd分別為：，并據此更新QSS模型參數和元件狀態；

步驟5根據式（1）～式（3），由zc(t+Δt)和zd(t+Δt)再次求解代數方程組，得到t+Δt時刻x和y的值x(t+Δt)、y(t+Δt)；

步驟6令t=t+Δt，若t≥T，則停止仿真；否則，返回步驟2。

2.2 基于潮流計算的代數方程求解

上述交替求解法步驟3和步驟5中，需求解代數方程式（1）～式（3）。該方程組可作為一般的非線性方程組，用牛頓法直接求解，但需要進行大量的編程開發。

從代數方程的具體形式來看，式（1）相對應的式（12）是代數方程的主體，在實際系統中占比通常超過50%，且去除某一發電機節點功率方程、并考慮發電機電壓/無功控制方程式（8）和式（12）后，與普通的潮流方程完全相同；式（2）對應的式（15）和式（16）在給定U時，P、Q與f之間為簡單的線性關系，而除去簡單的式（6），式（3）包含的式（9）在給定P、Q、U和θ求取Efd和δ時為一系列獨立的二元非線性方程組。針對上述特征，本文基于現有的潮流計算程序，求解代數方程式（1）～式（3），以提高程序的開發效率。

基于潮流計算的代數方程具體求解步驟如下：

步驟1由發電機OEL狀態及AVR參考電壓，確定發電機節點作為PV節點時的U或PQ節點時的Q，將無電壓控制能力的負荷等節點作為PQ節點，并選取某參與AGC調節且儲備較大的發電機節點作為平衡節點，設其相角為參考相角；

步驟2設相關待求變量的初始值（以上標0表示），置迭代計數器k=1；

步驟3將除平衡節點外的Pk?1、PQ節點的Qk?1以及PV節點和平衡節點的Uk?1代入除平衡節點功率方程的式（12），通過潮流計算，求取相關節點的Uk和θk，進而計算平衡節點的功率、PQ節點的無功功率以及系統的有功網損；

步驟4將Uk代入式（15），并假設系統有功網損不變，根據系統總發電與網損和總用電平衡的恒等關系，計算fk，進而由式（15）和式（16）更新Pk、Qk；

步驟5若變量P、Q、U、θ、f的第k-1次和第k次迭代值之差的最大絕對值大于允許誤差，且k小于最大迭代次數，則k=k+1，返回步驟3；

步驟6由式（9）計算各發電機的Efd和δ，停止迭代。

需說明的是，由于發電機能否維系機端電壓恒定取決于QSS模型式（5）確定的發電機OEL狀態，故上述潮流迭代計算過程中，不考慮發電機節點在PV節點和PQ節點之間的類型轉換。

圖1為本文所提基于潮流計算的代數方程求解方法，其中潮流計算由BPA程序完成，其余計算在Mathematica軟件中實現。

圖1 基于潮流計算的代數方程求解示意Fig.1 Schematic of solving algebraic equation based on power flow calculation

2.3 提升QSS仿真效率的方法

QSS仿真的總計算時間大致正比于總仿真時步數Nstep與單個時步計算時間Tstep的乘積，因此，為了提升QSS仿真的效率，可通過增加步長的方法減少Nstep，或通過改善迭代初始值的方法減少Tstep。

2.3.1 考慮離散動作特性的自適應變步長策略

QSS仿真步長的增加可減小Nstep，但會帶來一定的誤差。誤差的主要來源是數值積分方法的舍入誤差、截斷誤差和交替求解中的交接誤差。舍入誤差隨步長的增加而減小[24]，而QSS仿真步長通?？蛇_ 0.2～1.0 s[25?26]，故可忽略不計。交接誤差受限于BPA潮流計算輸出結果的精度，故本文不考慮。

截斷誤差與數值積分方法相關。采用改進歐拉積分法對y(t)積分時，截斷誤差ε與步長的三次方和y(t)三階導數y…(t)的絕對值成正比[24]，即

由于系統中長期動態主要受離散設備動作的影響，故在離散設備動作后，系統的連續動態逐步趨于平穩，從而積分函數三階導數y…(t)的絕對值也逐步減小。因此，為了平衡仿真效率和截斷誤差，可在故障或離散設備動作后，以小步長仿真，并隨仿真的推進，逐步增大仿真步長。據此，2次離散動作之間的變步長策略為

式中：hi和hi+1分別為第i步和第i+1步的步長；αi為變步長系數。

變步長系數αi的選取應使第i+1步的截斷誤差e小于等于最大允許截斷誤差εmax。為了確保該約束始終滿足，引入裕度系數0.9，從而有

由此，變步長系數αi最大為

在實際計算中y(t)的一階導數即式（4）已知，三階導數…y(t)可由前3個時刻的一階導數通過差分計算近似獲得。若前序時刻不足3，則αi設為1。由于QSS模型的zc為向量，對應于多個αi，故需取最小的αi作為QSS仿真的變步長系數。

將2次離散動作時刻之間的變步長策略式（19）與離散動作時刻的仿真步長要求相結合，可得如下考慮離散動作特性的自適應變步長策略：

（1）設最小和最大仿真步長分別為hmin和hmax，并以仿真步長hmin開始仿真；

（2）第i步仿真完成后，計算變步長系數αi；

（3）由離散設備動作條件、時延設定及計時器，求取下一離散動作發生時刻與當前時刻的差；

（4）第i+1時步的步長為；

（5）離散動作發生后，重置hi+1=hmin。

2.3.2 基于多時步預測的初始值設定方法

由于代數方程式（1）～式（3）的迭代求解時間主導了QSS仿真單時步的計算時間Tstep，故減少Tstep的關鍵是減少代數方程的迭代求解次數。顯然，改善代數方程迭代求解的初始值可達到這一目標，同時，這也意味著減少調用潮流計算程序的次數，從而節約了較為耗時的外部程序調用時間。

考慮到在2次離散動作之間的仿真時段內狀態變量的變化主要受AGC、負荷恢復等影響，在有限的時間內，中長期動態變化通常不大，故可用合適階數的多項式擬合先前時步的結果，再外推預測下一時步的變量值。具體地，在離散動作發生后，本文采用如下基于多時步預測的初始值設定策略：

（1）第1時步，以離散動作前的值作為初始值；

（2）第2時步，以前1時步的計算結果作為初始值；

（3）第3時步，以前2時步的線性預測作為初始值；

（4）第4時步及以后，以前3時步的二次預測作為初始值。

2.4 算法總流程

圖2為本文所提基于潮流計算的中長期頻率電壓穩定仿真算法總流程。

圖2 基于潮流計算的中長期仿真方法算法流程Fig.2 Flow chart of mid-long term simulation method based on power flow calculation

3 仿真算例

3.1 Nordic32測試系統

該系統為中長期電壓穩定研究常用的74節點算例系統，其單線圖如圖3所示[19]。本算例中，假設所有發電機都參與AGC調節，發電機G7和G17裝有過勵限制器，并考慮負荷3、4、41的恢復特性及其對應變壓器的OLTC控制。本算例仿真條件為CPU=i7?6700@3.4GHz，RAM=4 GB。

圖3 Nordic32系統單線圖Fig.3 Single-line diagram of Nordic32 System

設系統擾動為負荷功率突然增加5%，分別采用全時域仿真方法和3種QSS仿真方法（定步長、變步長以及變步長加多時步預測）（分別簡稱為方法1～4）對系統頻率變化以及節點41和節點3的電壓變化進行仿真分析，仿真結果如圖4所示。其中，定步長方法的步長為0.5 s，變步長方法的hmin=0.5 s、hmax=4.0 s。

圖4 不同方法的仿真結果對比Fig.4 Comparison of simulation results among different methods

由圖4（a）可見，系統頻率在擾動后由50 Hz迅速降低至49.89 Hz，隨后在AGC的作用下，調頻機組增加出力，頻率逐步恢復，約30 s達到額定值。由圖4（b）和4（c）可見，發電機G17和G7分別在t=16 s和29 s時觸發過勵限制，使得負荷節點3和節點41的電壓低于對應OLTC的控制死區下限0.95，從而節點3對應OLTC在t=39 s時增加一個檔位，節點41對應OLTC在t=41 s和69 s時增加一個檔位，使得相關負荷節點電壓恢復到0.95以上。

比較圖4中的各條仿真曲線可見，全時域仿真結果較好地反映了系統的頻率和電壓暫態，而3種QSS方法的仿真結果非常接近，雖然不能反映系統的短期暫態，但均較準確地反映了系統的中長期動態，驗證了本文所提中長期頻率電壓動態過程QSS仿真方法的準確性。

表2比較了4種仿真方法的計算時間、仿真時步數以及3種QSS方法的總迭代次數和每時步平均迭代次數。由表可見：與全時域仿真方法相比，QSS方法大大減小了計算時間和仿真步數；與定步長方法相比，變步長QSS方法的計算效率可提高近3倍；與變步長方法相比，變步長多時步預測QSS方法的仿真時步數不變，但平均迭代次數減小，從而減小了總迭代次數、進而提高了計算效率。

表2 不同方法的仿真性能Tab.2 Simulation performance of different methods

綜上，與通常的定步長QSS仿真方法相比，本文所提的變步長?多時步預測QSS仿真方法具有更好的仿真效率。

為了研究hmax對仿真性能的影響，在其他條件相同的情況下，分別以hmax=2、4、8 s進行仿真測試，其對比結果如圖5所示。由圖可見，3種hmax的變步長仿真結果均能夠正確地反映中長期動態變化，但在hmax=8 s時，連續變化時間段內的部分曲線細節有所損失，即仿真準確度下降。

圖5 hmax=2,4,8 s時的QSS仿真結果對比Fig.5 Comparison of QSS simulation results athmaxof 2，4，and 8 s

表3比較了hmax=2、4、8 s下QSS仿真的性能。由表可見，當hmax從2 s增大到4 s以及從4 s增大到8 s時，計算時間分別減少2.78 s和1.13 s，仿真時步數分別減少26次和7次。顯然，后者提升的仿真效率不如前者，其原因在于，本文的變步長策略受限于離散設備的動作時間，hmax的倍增并不意味著仿真時步數的倍減，從而計算效率顯著提升。

表3 hmax=2,4,8 s的QSS仿真結果對比Tab.3 Comparison of QSS simulation results athmaxof 2，4，and 8 s

綜上，在本算例中，hmax=4 s時較好地平衡了仿真準確度和仿真效率。

3.2 某大區電網等值3762節點系統

假設系統中對應于某省的401個負荷均有恢復特性，75臺發電機均裝有過勵限制器并參與AGC調節，并設其112臺變壓器為OLTC，動作時延為8～12 s。圖6為系統在t=1 s時直流單回閉鎖（損失3 000 MW功率）條件下系統頻率以及部分節點電壓、OLTC檔位、發電機無功功率的中長期動態變化曲線。算例由變步長+多時步預測QSS仿真方法完成，總仿真時步數162，計算時間約75 s。本算例仿真條件為CPU=i7?6700@3.4GHz，RAM=4 GB。

圖6 某大區電網系統中長期仿真結果Fig.6 Mid-long term simulation results of one large-scale power grid system

由圖6（a）和6（b）可見，故障發生后，系統頻率驟降至49.56 Hz，部分節點電壓亦下降0.02～0.08 p.u.不等，且由于負荷恢復特性的影響，電壓仍有進一步下降的趨勢；由圖6（c）可見，OLTC由于目標節點電壓持續低于其控制死區下限而連續動作，直至達到控制極限；由圖6（d）可見，受擾動影響，部分發電機持續增大無功出力，在觸發OEL后，失去電壓支撐，導致如圖6（b）在50 s的節點電壓驟降。

4 結語

本文針對電力系統中長期頻率電壓穩定動態仿真問題，提出了基于潮流計算的中長期準穩態時域仿真方法，降低了仿真實現的開發難度；實現了自適應變步長方法和多時步預測的初始值設定方法，提高了仿真效率。對Nordic32系統和某大區電網系統的仿真結果表明，該方法能滿足大規模系統的分析需求，具有計算準確、實現簡便的優點。

隨著大規模直流輸電、可再生能源接入、需求側響應等技術的廣泛應用，有必要在電力系統中長期頻率電壓穩定仿真分析中考慮這些因素。建立其準穩態仿真模型并研究其高效仿真方法是本文后續進一步的工作。