佘 東，嚴(yán) 正，王 晗，蒙 飛，李宏強(qiáng)，徐瀟源，戴玉臣

（1. 電力傳輸與功率變換控制教育部重點(diǎn)實(shí)驗(yàn)室（上海交通大學(xué)），上海市 200240；2. 國(guó)網(wǎng)寧夏電力有限公司，寧夏回族自治區(qū)銀川市 750010；3. 南瑞集團(tuán)有限公司（國(guó)網(wǎng)電力科學(xué)研究院有限公司），江蘇省南京市 211106）

0 引言

電力系統(tǒng)暫態(tài)分析主要研究系統(tǒng)受大干擾后的同步運(yùn)行能力與隨時(shí)間變化的動(dòng)態(tài)行為［1-2］。目前，時(shí)域仿真法（又稱(chēng)逐步積分法）是應(yīng)用最為廣泛的電力系統(tǒng)暫態(tài)分析方法，其通過(guò)建立并求解描述系統(tǒng)動(dòng)態(tài)行為的微分-代數(shù)方程組，以獲得系統(tǒng)狀態(tài)變量和代數(shù)變量隨時(shí)間變化的動(dòng)態(tài)曲線，進(jìn)而根據(jù)發(fā)電機(jī)轉(zhuǎn)子間的相對(duì)搖擺曲線判別系統(tǒng)穩(wěn)定性［3］。隨著電力系統(tǒng)中不確定性因素的增加，傳統(tǒng)確定性時(shí)域仿真計(jì)算逐漸向時(shí)域仿真不確定性分析轉(zhuǎn)變，以滿足強(qiáng)隨機(jī)性、波動(dòng)性場(chǎng)景下電力系統(tǒng)暫態(tài)分析的需要。計(jì)及不確定性因素的影響，如何快速、準(zhǔn)確地實(shí)現(xiàn)電力系統(tǒng)時(shí)域仿真也成為現(xiàn)階段研究的重點(diǎn)。

電力系統(tǒng)時(shí)域仿真不確定性分析具有時(shí)變性、非線性、隨機(jī)性的特點(diǎn)，針對(duì)該問(wèn)題的主要研究方法有蒙特卡洛模擬（Monte Carlo simulation，MCS）法［4-5］、隨 機(jī) 響 應(yīng) 面 法（stochastic response surface method，SRSM）［6］、多 項(xiàng) 式 混 沌 展 開(kāi)（polynomial chaos expansion，PCE）方法［7-9］、區(qū)間算法［10］、概率配點(diǎn)法［11-12］、軌跡靈敏度法［13-14］，隨機(jī)方程法［15-16］等。MCS 法旨在通過(guò)隨機(jī)變量的統(tǒng)計(jì)試驗(yàn)和隨機(jī)模擬求解不確定性因素影響下的隨機(jī)問(wèn)題，其原理簡(jiǎn)單且可保證計(jì)算結(jié)果的準(zhǔn)確性，但計(jì)算代價(jià)大。SRSM 和PCE 方法的基本原理一致，均是利用正交多項(xiàng)式的混沌展開(kāi)之和來(lái)建立代理模型，以擬合原模型的概率輸出響應(yīng)結(jié)果；相較于MCS 法，SRSM和PCE 方法避免了原模型中反復(fù)的數(shù)值積分計(jì)算，在盡可能保證計(jì)算精度的同時(shí)，提高了電力系統(tǒng)時(shí)域仿真的計(jì)算效率。然而，SRSM 和PCE 方法應(yīng)用于電力系統(tǒng)時(shí)域仿真不確定性分析時(shí)，無(wú)法考慮輸出變量在時(shí)間維度上的連續(xù)性，構(gòu)建代理模型的過(guò)程需依賴仿真時(shí)間內(nèi)多個(gè)獨(dú)立時(shí)間斷面的輸入-輸出樣本，樣本獲取過(guò)程復(fù)雜。同時(shí)，隨著不確定性因素的增加，SRSM 和PCE 方法會(huì)遇到維數(shù)災(zāi)問(wèn)題，難以應(yīng)用于含高維隨機(jī)變量的問(wèn)題［17］。區(qū)間算法利用區(qū)間數(shù)描述模型參數(shù)的不確定性，并采用區(qū)間泰勒級(jí)數(shù)法求解區(qū)間微分方程組，用以分析模型參數(shù)不確定性對(duì)輸出結(jié)果的影響，但采用該方法所得結(jié)果往往存在著一定的保守性［10］。概率配點(diǎn)法基于廣義混沌多項(xiàng)式原理進(jìn)行采樣［18］，其求解精度受配點(diǎn)方式的影響較大。軌跡靈敏度法通過(guò)在額定運(yùn)行點(diǎn)處建立表征系統(tǒng)輸入-輸出關(guān)系的一階線性近似表達(dá)式，分析參數(shù)變化對(duì)系統(tǒng)狀態(tài)變量的影響，但該方法難以應(yīng)用于非線性較強(qiáng)的系統(tǒng)［19］。隨機(jī)方程法針對(duì)系統(tǒng)中的隨機(jī)擾動(dòng)，采用隨機(jī)過(guò)程對(duì)系統(tǒng)動(dòng)態(tài)過(guò)程建模，并引入Euler 或Milstein 等方法實(shí)現(xiàn)隨機(jī)微分方程的求解，但其原理較為復(fù)雜，計(jì)算也較為煩瑣。

近年來(lái)，文獻(xiàn)［19-22］提出了基于Galerkin 法的多項(xiàng)式逼近方法，并將其應(yīng)用于電力系統(tǒng)時(shí)域仿真不確定性分析中，相較于SRSM 和PCE 方法，基于Galerkin 法的多項(xiàng)式逼近方法引入了時(shí)序項(xiàng)，能夠通過(guò)數(shù)值積分實(shí)現(xiàn)系統(tǒng)動(dòng)態(tài)仿真軌跡的遞推計(jì)算，避免了在代理模型構(gòu)建過(guò)程中需要大量多時(shí)間斷面采樣的問(wèn)題。但該方法同樣面臨維數(shù)災(zāi)問(wèn)題，需要借助其他方法予以解決，如靈敏度分析方法等。

針對(duì)現(xiàn)有方法中存在的不足，本文提出了一種基于含時(shí)序項(xiàng)代理模型的電力系統(tǒng)時(shí)域仿真不確定性分析方法。該方法利用非線性有源自回歸多項(xiàng)式混沌展開(kāi)（nonlinear autoregressive with exogenous input-polynomial chaos expansion，NARX-PCE）方法構(gòu)建了含時(shí)序項(xiàng)的代理模型，兼顧了不確定性因素的影響和時(shí)域仿真過(guò)程的動(dòng)態(tài)特性。其中，時(shí)序項(xiàng)是指由系統(tǒng)不同時(shí)刻點(diǎn)的輸出變量所組成的代理模型基函數(shù)項(xiàng)，以此反映系統(tǒng)輸出隨時(shí)間動(dòng)態(tài)變化的情況。此外，在代理模型的構(gòu)建中，本文引入了最小角回歸（least angle regression，LAR）策略，以克服高維隨機(jī)變量所引起的維數(shù)災(zāi)問(wèn)題。

1 考慮不確定性因素的電力系統(tǒng)時(shí)域仿真方法

1.1 電力系統(tǒng)隨機(jī)暫態(tài)模型

考慮不確定性因素的影響，建立刻畫(huà)電力系統(tǒng)動(dòng)態(tài)過(guò)程的隨機(jī)微分-代數(shù)方程，如式（1）所示［18］。

式中:x為表征系統(tǒng)動(dòng)態(tài)特性的狀態(tài)變量，如發(fā)電機(jī)轉(zhuǎn)子角等；y為表征系統(tǒng)運(yùn)行狀態(tài)的代數(shù)變量，如電壓幅值等；ξ為輸入隨機(jī)變量，包含發(fā)電機(jī)端的隨機(jī)出力、系統(tǒng)負(fù)荷水平等；f（?）和g（?）分別為描述系統(tǒng)運(yùn)行特性的微分方程組與代數(shù)方程組。

在確定性場(chǎng)景下，基于數(shù)值積分算法的時(shí)域仿真方法是一種求解微分-代數(shù)方程組的有效方法。針對(duì)式（1）所示的隨機(jī)微分-代數(shù)方程，在時(shí)域仿真方法的基礎(chǔ)上，可結(jié)合不確定性分析手段，獲得不確定性因素影響下的電力系統(tǒng)狀態(tài)變量和代數(shù)變量隨時(shí)間變化的情況，以評(píng)估不確定性場(chǎng)景下系統(tǒng)運(yùn)行的穩(wěn)定性與安全性。

1.2 電力系統(tǒng)時(shí)域仿真不確定性分析方法

目前，常用的電力系統(tǒng)時(shí)域仿真不確定性分析方法主要包括MCS 法和代理模型方法。

1.2.1 MCS 法

MCS 法是一種常用的不確定性分析方法，其通過(guò)采集大量的輸入隨機(jī)變量樣本，將隨機(jī)問(wèn)題轉(zhuǎn)化為一系列的確定性問(wèn)題求解，最終獲得系統(tǒng)輸出響應(yīng)的分布結(jié)果［23］。MCS 法原理簡(jiǎn)單且易于實(shí)現(xiàn)，但大量的重復(fù)計(jì)算在提高計(jì)算精度的同時(shí)也極大降低了計(jì)算效率?；贛CS 法的電力系統(tǒng)時(shí)域仿真不確定分析方法在獲得系統(tǒng)輸出變量動(dòng)態(tài)軌跡分布的同時(shí)，將耗費(fèi)大量計(jì)算成本用于反復(fù)的數(shù)值積分計(jì)算，難以實(shí)現(xiàn)電力系統(tǒng)動(dòng)態(tài)過(guò)程的快速分析。因此，基于MCS 法的時(shí)域仿真不確定性分析方法常作為驗(yàn)證其他方法準(zhǔn)確性的基準(zhǔn)方法使用［17］。

1.2.2 代理模型方法

基于代理模型的不確定性分析方法是指通過(guò)構(gòu)建表示原系統(tǒng)輸入-輸出關(guān)系的代理模型進(jìn)行計(jì)算，在盡可能保證計(jì)算精度的條件下，降低原系統(tǒng)模型不確定性分析過(guò)程的計(jì)算負(fù)擔(dān)。常用于電力系統(tǒng)時(shí)域仿真不確定性分析的代理模型方法包括PCE 方法和基于Galerkin 法的多項(xiàng)式逼近方法。

1）PCE 方法

PCE 方法利用含獨(dú)立隨機(jī)變量的正交多項(xiàng)式之和建立原系統(tǒng)的代理模型，并通過(guò)配點(diǎn)法或最小二乘法求解多項(xiàng)式系數(shù)，以獲得表征原系統(tǒng)輸入-輸出關(guān)系的表達(dá)式。基于PCE 方法，系統(tǒng)輸出變量Y的代理模型可表示為:

式中:?i（ξ）為含隨機(jī)變量ξ的正交多項(xiàng)式；αi為對(duì)應(yīng)于第i個(gè)多項(xiàng)式的系數(shù)。

當(dāng)應(yīng)用于電力系統(tǒng)時(shí)域仿真不確定性分析時(shí)，PCE 方法以系統(tǒng)輸出變量動(dòng)態(tài)過(guò)程的各個(gè)時(shí)間斷面為采樣點(diǎn)獲取輸入-輸出樣本，對(duì)各個(gè)時(shí)間斷面建立代理模型。在給定新的系統(tǒng)輸入隨機(jī)變量樣本時(shí)，基于PCE 方法的代理模型可計(jì)算得到各時(shí)間斷面輸出變量的結(jié)果，并最終形成不確定性場(chǎng)景下系統(tǒng)輸出變量的動(dòng)態(tài)響應(yīng)軌跡，以作為原系統(tǒng)數(shù)值積分計(jì)算過(guò)程的近似結(jié)果。需要說(shuō)明的是，常用的SRSM［24］本質(zhì)上是PCE 方法采用標(biāo)準(zhǔn)正態(tài)分布隨機(jī)變量為輸入變量、Hermite 正交多項(xiàng)式為基函數(shù)時(shí)的一種特殊形式。

2）基于Galerkin 法的多項(xiàng)式逼近方法

基于Galerkin 法的多項(xiàng)式逼近方法通過(guò)對(duì)式（1）中的狀態(tài)變量與代數(shù)變量進(jìn)行多項(xiàng)式展開(kāi)，得到對(duì)應(yīng)的代理模型［19-21］:

式中:xj(t)和yj(t)分別為t時(shí)刻x和y的第j個(gè)元素；αi，j，x(t)和αi，j，y(t)分 別 為t時(shí) 刻xj(t)和yj(t)進(jìn)行多項(xiàng)式展開(kāi)后對(duì)應(yīng)第i個(gè)多項(xiàng)式的系數(shù)。

通過(guò)構(gòu)造Galerkin 投影方程可得到含時(shí)間變量的確定性方程，進(jìn)而利用數(shù)值積分算法計(jì)算獲得不確定性場(chǎng)景下系統(tǒng)輸出變量的動(dòng)態(tài)響應(yīng)軌跡。

表1 對(duì)比了上述3 種電力系統(tǒng)時(shí)域仿真不確定性分析方法的特點(diǎn)。相較于MCS 法，PCE 方法和基于Galerkin 法的多項(xiàng)式逼近方法可通過(guò)構(gòu)建代理模型提升電力系統(tǒng)時(shí)域仿真不確定性分析過(guò)程的效率。相較于依賴多斷面樣本構(gòu)建代理模型的PCE方法，基于Galerkin 法的多項(xiàng)式逼近方法含有時(shí)間變量，可在確定多項(xiàng)式形式后進(jìn)行數(shù)值積分計(jì)算，無(wú)須進(jìn)行多斷面采樣。但PCE 方法和基于Galerkin 法的多項(xiàng)式逼近方法均存在維數(shù)災(zāi)問(wèn)題，隨著輸入隨機(jī)變量維度的增加，代理模型構(gòu)建過(guò)程的計(jì)算負(fù)擔(dān)將迅速增大，并嚴(yán)重影響代理模型方法在時(shí)域仿真不確定性分析中的應(yīng)用。

表1 電力系統(tǒng)時(shí)域仿真不確定性分析方法對(duì)比Table 1 Comparison of uncertainty analysis methods for time-domain simulation of power system

2 含時(shí)序項(xiàng)代理模型的構(gòu)建與求解

2.1 基于NARX-PCE 方法的含時(shí)序項(xiàng)代理模型

為了解決現(xiàn)有代理模型方法所遇到的維數(shù)災(zāi)問(wèn)題，并盡可能降低代理模型構(gòu)建對(duì)電力系統(tǒng)動(dòng)態(tài)過(guò)程各時(shí)間斷面采樣的依賴，本文提出了基于NARXPCE 方法的含時(shí)序項(xiàng)代理模型，并將其用于電力系統(tǒng)時(shí)域仿真不確定性分析中。

2.1.1 考慮不確定性因素的NARX 模型構(gòu)建方法

NARX 方法常用于表示系統(tǒng)某一時(shí)刻的輸出與系統(tǒng)歷史時(shí)刻輸出以及外部激勵(lì)信號(hào)之間的關(guān)系［25-26］，其所構(gòu)建的模型可表示如下:

式中:Y（t）為系統(tǒng)t時(shí)刻的輸出，對(duì)于電力系統(tǒng)而言，可為相應(yīng)狀態(tài)變量以及代數(shù)變量等；u（t）=［X（t），X（t-1），…，X（t-tX），Y（t-1），Y（t-2），…，Y（t-tY）］，其中X（t），X（t-1），…，X（t-tX）分別為t，t-1，…，t-tX時(shí)刻的外部輸入激勵(lì)，Y（t-1），Y（t-2），…，Y（t-tY）分別為系統(tǒng)在t-1，t-2，…，t-tY時(shí)刻（t時(shí)刻之前的歷史時(shí)刻）的輸出，tX和tY分別為所考慮的外部輸入激勵(lì)、系統(tǒng)歷史時(shí)刻輸出的最大時(shí)延；θi為NARX 模型中的待求系數(shù)；gi(u(t))為時(shí)序項(xiàng)基函數(shù)；Ng為模型中基函數(shù)的項(xiàng)數(shù)；εt為NARX 模型的殘差，一般假設(shè)其服從均值為0、方差為σ（t）的正態(tài)分布。

上述模型從系統(tǒng)動(dòng)態(tài)特性的角度表征了系統(tǒng)輸出隨時(shí)間變化的情況。進(jìn)一步考慮不確定性因素對(duì)系統(tǒng)輸出的影響，式（4）可擴(kuò)展為:

式中:Y(t，ξ)為輸入隨機(jī)變量ξ影響下系統(tǒng)t時(shí)刻的輸出；θi(ξ)為含輸入隨機(jī)變量ξ的系數(shù)；εt(t，ξ)為不確定性因素影響下NARX 模型在t時(shí)刻的殘差。

2.1.2 NARX-PCE 模型

由式（5）可知，NARX 模型中的系數(shù)θi(ξ)含有輸入隨機(jī)變量ξ，可引入PCE 模型表征輸入隨機(jī)變量對(duì)待求系數(shù)的影響:

式 中:αi，j，θ為PCE 模 型 中 的 待 求 系 數(shù)；i=1，2，…，Ng；φj(ξ)為關(guān)于輸入隨機(jī)變量ξ的多元正交多項(xiàng)式的基函數(shù)；Nφ為正交多項(xiàng)式的項(xiàng)數(shù)；εi為PCE 模型的截?cái)嗾`差。

結(jié)合式（5）和式（6），可得到基于NARX-PCE方法的含時(shí)序項(xiàng)代理模型，具體展開(kāi)式如下:

式中:ε(t，ξ)為NARX-PCE 方法所構(gòu)建代理模型的總誤差。

上述基于NARX-PCE 方法的含時(shí)序項(xiàng)代理模型具有以下2 個(gè)特點(diǎn):1）所含時(shí)序項(xiàng)可表征系統(tǒng)隨時(shí)間變化的動(dòng)態(tài)特性；2）所含隨機(jī)參數(shù)可表征系統(tǒng)受不確定性因素影響的隨機(jī)特性。以上特點(diǎn)使得基于NARX-PCE 方法的含時(shí)序項(xiàng)代理模型可用于不確定性環(huán)境下電力系統(tǒng)的時(shí)域仿真分析。

2.2 含時(shí)序項(xiàng)代理模型基函數(shù)確定與待求系數(shù)計(jì)算

結(jié)合式（5）—式（7）可知，構(gòu)建基于NARX-PCE方法的含時(shí)序項(xiàng)代理模型需要確定NARX 模型中的時(shí)序項(xiàng)基函數(shù)gi(u(t))、PCE 模型中正交多項(xiàng)式基函數(shù)φj(ξ)以及相應(yīng)的待求系數(shù)αi，j，θ。

2.2.1 基于LAR 策略的基函數(shù)選取方法

LAR 策 略 由Efron［27］于2004 年 提 出，用 于 從 待選基函數(shù)集合中依次篩選出與目標(biāo)向量相關(guān)性最大的基函數(shù)，以降低目標(biāo)向量量化表征過(guò)程中所需基函數(shù)的數(shù)量，克服該過(guò)程中可能遇到的過(guò)擬合和維數(shù)災(zāi)問(wèn)題。本文在確定NARX 模型中的時(shí)序項(xiàng)基函數(shù)gi(u(t))和PCE 模型中正交多項(xiàng)式基函數(shù)φj(ξ)時(shí)均利用了LAR 策略，具體方法如下。

1）NARX 模型中時(shí)序項(xiàng)基函數(shù)gi(u(t))的選取

首先，給定NARX 模型中待選時(shí)序項(xiàng)基函數(shù)的最大階數(shù)m、最大時(shí)延tY，生成NARX 模型時(shí)序項(xiàng)待選基函數(shù)集合。本文暫不考慮電力系統(tǒng)時(shí)域仿真過(guò)程中系統(tǒng)外部激勵(lì)的影響，因此，時(shí)序項(xiàng)待選基函數(shù)的構(gòu)建將不包含X（t）?；谳敵鲎兞縔（t）的歷史數(shù)據(jù)，待選基函數(shù)個(gè)數(shù)與最大階數(shù)m、最大時(shí)延數(shù)tY的關(guān)系如附錄A 表A1 所示。在實(shí)際應(yīng)用中，待選基函數(shù)的階數(shù)m一般不超過(guò)3，以降低階數(shù)過(guò)高導(dǎo)致的過(guò)擬合現(xiàn)象［26］。以m=3、tY=3 為例，待選基函數(shù)共包含19 項(xiàng)，如表2 所示。

表2 m=3、tY=3 時(shí)待選時(shí)序項(xiàng)基函數(shù)集合Table 2 Set of time-series-term basis functions to be selected with m=3，tY=3

其次，生成N個(gè)輸入隨機(jī)變量樣本ξS={ξ(1)，ξ(2)，…，ξ(N)}。對(duì) 于t=1，2，…，TS，其 中TS為采樣區(qū)間的時(shí)間斷面數(shù)，通過(guò)電力系統(tǒng)時(shí)域仿真計(jì)算獲得第k個(gè)輸入隨機(jī)變量樣本ξ(k)所對(duì)應(yīng)的輸出變量樣本Y(t，ξ)（TS＞tY）。

再次，在給定樣本ξ(k)下，將每個(gè)時(shí)間斷面tS（TS≥tS＞tY）的樣本點(diǎn)代入式（5）展開(kāi)，獲得如式（8）所示的矩陣形式。

進(jìn)一步，將式（8）寫(xiě)成緊湊形式為:

式中:Yk、ψk、θ(ξ(k))、ε(ξ(k))分別為代入第k個(gè)輸入隨機(jī)變量樣本ξ(k)后計(jì)算獲得的輸出變量向量、基函數(shù)信息矩陣、待求系數(shù)向量、殘差向量。

最后，給定需要選取的基函數(shù)數(shù)量Ng或擬合誤差εNARX，利用LAR 策略確定滿足數(shù)量要求或誤差要求的NARX 模型時(shí)序項(xiàng)基函數(shù)，具體過(guò)程如附錄B所示。需要說(shuō)明的是，基于LAR 策略選擇基函數(shù)的過(guò)程采用了LOO（leave-one-out）誤差進(jìn)行交叉驗(yàn)證，以提高代理模型的泛化能力［28］。其中，LOO 誤差的計(jì)算公式可見(jiàn)附錄C。

2）PCE 模型正交多項(xiàng)式基函數(shù)φj（ξ）的選取

在確定NARX 模型中時(shí)序項(xiàng)基函數(shù)后，可采用最小二乘法（具體算法將在2.2.2 節(jié)給出）獲得N個(gè)輸入隨機(jī)變量樣本ξS對(duì)應(yīng)的待求系數(shù)樣本矩陣。

以上述待求系數(shù)樣本矩陣的第i列θi(ξ(N))為輸出變量樣本，在給定正交多項(xiàng)式基函數(shù)選取的數(shù)量Nφ或擬合誤差εPCE的條件下，利用LAR 策略選取如式（6）所示的PCE 模型中的正交多項(xiàng)式基函數(shù)，其中，待選基函數(shù)類(lèi)型包括了Hermite 基函數(shù)、Legendre 基函 數(shù)、Laguerre 基函數(shù)等［29］。PCE 模型已在電力系統(tǒng)中有較多應(yīng)用，更為詳細(xì)的基函數(shù)選取流程可參見(jiàn)文獻(xiàn)［30-31］。

2.2.2 待求系數(shù)計(jì)算方法

1）NARX 模型待求系數(shù)θi(ξ)的計(jì)算

基于2.2.1 節(jié)中選取的NARX 模型的時(shí)序項(xiàng)基函數(shù)，依次將輸入隨機(jī)變量樣本ξ(k)（k=1，2，…，N）代入計(jì)算，獲得基函數(shù)信息矩陣ψk和輸出向量Yk；進(jìn)而，利用最小二乘法計(jì)算待求系數(shù)向量θ(ξ(k))=[θ1(ξ(k))，θ2(ξ(k))，…，θNg(ξ(k))]T。

2）PCE 模型待求系數(shù)αi，j，θ的計(jì)算

基于2.2.1 節(jié)中選取的PCE 模型正交多項(xiàng)式基函數(shù)，代入輸入隨機(jī)變量樣本ξ(k)（k=1，2，…，N），利用最小二乘法計(jì)算對(duì)應(yīng)θi（ξ）的待求系數(shù)向量

式中:φ為代入輸入隨機(jī)變量樣本ξ(k)（k=1，2，…，N）后正交多項(xiàng)式基函數(shù)的信息矩陣；θi(ξ)為式（10）中待求系數(shù)樣本矩陣的第i列，其中，式（10）中的各項(xiàng)元素已通過(guò)式（11）計(jì)算得到。

2.3 含時(shí)序項(xiàng)代理模型構(gòu)建步驟

一方面，含時(shí)序項(xiàng)代理模型利用NARX 模型表征了原模型輸出與歷史輸出之間的時(shí)間相關(guān)性，能夠刻畫(huà)系統(tǒng)的動(dòng)態(tài)特性；另一方面，采用PCE 模型表征了不確定性因素對(duì)原模型輸出的影響。

在含時(shí)序項(xiàng)代理模型構(gòu)建過(guò)程中，LAR 策略能夠有效應(yīng)對(duì)高階基函數(shù)和高維輸入隨機(jī)變量所帶來(lái)的過(guò)擬合和維數(shù)災(zāi)問(wèn)題。

在給定新的不確定性場(chǎng)景ξ(N+1)下，首先，利用傳統(tǒng)PCE 方法獲得時(shí)延范圍內(nèi)各時(shí)間斷面所對(duì)應(yīng)的輸出變量樣本Y(t，ξ(N+1))（t=1，2，…，tY），然后，利用本文所提出的含時(shí)序項(xiàng)代理模型逐步隨時(shí)間（t=tY+1，tY+2，…）遞推計(jì)算輸出變量，從而實(shí)現(xiàn)任意給定不確定性場(chǎng)景下輸出響應(yīng)的模擬和預(yù)測(cè)。綜上，利用所提含時(shí)序項(xiàng)代理模型對(duì)系統(tǒng)時(shí)域仿真中輸出響應(yīng)進(jìn)行不確定性分析時(shí)，主要計(jì)算步驟如圖1 所示。

圖1 基于含時(shí)序項(xiàng)代理模型的電力系統(tǒng)時(shí)域仿真不確定性分析流程圖Fig.1 Flow chart of uncertainty analysis of time-domain simulation of power system based on surrogate model containing time series terms

3 算例分析

本文采用IEEE 9 節(jié)點(diǎn)和118 節(jié)點(diǎn)系統(tǒng)進(jìn)行分析，并將所提基于NARX-PCE 方法的含時(shí)序項(xiàng)代理模型計(jì)算結(jié)果與傳統(tǒng)MCS 法、PCE 方法求解結(jié)果進(jìn)行對(duì)比，以驗(yàn)證所提方法的優(yōu)越性。本文測(cè)試環(huán)境為Intel Core i5-9300H 四核CPU、8 GB 內(nèi)存，仿真程序的編譯與測(cè)試在軟件PSCAD X4（4.6）、PSS@E 及其相應(yīng)的Python API 和MATLAB R2019a 上完成。

3.1 IEEE 9 節(jié)點(diǎn)系統(tǒng)仿真結(jié)果

IEEE 9 節(jié)點(diǎn)系統(tǒng)接線圖如附錄D 圖D1 所示，發(fā)電機(jī)（G1、G2 和G3）采用PSCAD 中的高階模型，部分模型參數(shù)可參見(jiàn)文獻(xiàn)［32］。在線路7-8 之間設(shè)置三相接地短路故障，故障開(kāi)始時(shí)間為t=0.5 s，故障持續(xù)時(shí)間為0.08 s，即t=0.58 s 時(shí)故障清除。系統(tǒng)穩(wěn)態(tài)運(yùn)行下，假定發(fā)電機(jī)G1（節(jié)點(diǎn)1）和發(fā)電機(jī)G3（節(jié)點(diǎn)3）的有功功率P1和P3為輸入隨機(jī)變量，其與發(fā)電機(jī)出力基準(zhǔn)值PG1和PG3之間的關(guān)系為:P1=（1+ξ1）PG1；P3=（1+ξ2）PG3，其中，ξ1和ξ2為均值為0.1、方差為0.02 的正態(tài)分布隨機(jī)變量。

3.1.1 算法對(duì)比

基于上述仿真數(shù)據(jù)，選取發(fā)電機(jī)G1 和G2 之間的相對(duì)功角差δ1-2為輸出變量，分別生成暫態(tài)功角動(dòng)態(tài)特性曲線的訓(xùn)練樣本集（2 000 組）和測(cè)試樣本集（100 組）。

分 別 采 用MCS 法、PCE 方 法 和NARX-PCE 方法開(kāi)展電力系統(tǒng)時(shí)域仿真不確定性分析，以驗(yàn)證本文所提方法的準(zhǔn)確性和高效性。其中，MCS 法的計(jì)算結(jié)果作為基準(zhǔn)值；PCE 方法可利用Uqlab 軟件［33］實(shí)現(xiàn)，展開(kāi)式最高階數(shù)設(shè)為5，輸出動(dòng)態(tài)響應(yīng)軌跡的時(shí)間斷面共取300 個(gè)（即需要構(gòu)建300 個(gè)PCE 代理模型）；基于NARX-PCE 方法的含時(shí)序項(xiàng)代理模型中最大延時(shí)參數(shù)tY=10，最大階數(shù)m=3，采樣區(qū)間內(nèi)的時(shí)間斷面數(shù)TS=60，待選時(shí)序項(xiàng)基函數(shù)為285 個(gè)。給定2 000 個(gè)輸入隨機(jī)變量樣本，利用PCE方法和NARX-PCE 方法所構(gòu)建的代理模型進(jìn)行電力系統(tǒng)時(shí)域仿真計(jì)算，所得相對(duì)功角差δ1-2的均值和方差隨時(shí)間動(dòng)態(tài)變化的軌跡分別如圖2（a）和圖2（b）所示。

圖2 相對(duì)功角差均值和方差的動(dòng)態(tài)軌跡Fig.2 Dynamic trajectories of mean value and variance of relative power angle difference

由仿真結(jié)果可知，基于NARX-PCE 方法所構(gòu)建代理模型的時(shí)域仿真結(jié)果與MCS 法、PCE 方法的計(jì)算結(jié)果基本一致，驗(yàn)證了所提方法的準(zhǔn)確性。其中，利用LAR 策略、基于NARX-PCE 方法所構(gòu)建的含時(shí)序項(xiàng)代理模型最終確定了6 個(gè)時(shí)序項(xiàng)基函數(shù):Y（t-1）、Y（t-10）、Y（t-1）Y（t-1）、Y（t-10）·Y（t-10）、Y（t-1）Y（t-1）Y（t-1）、Y（t-9）Y（t-10）Y（t-10），其LOO 誤差值為1.732 3×10-4。需要說(shuō)明的是，為方便對(duì)比不同方法所得時(shí)域仿真不確定性分析結(jié)果，本文3.1 節(jié)中將故障清除時(shí)刻作為0 時(shí)刻（時(shí)間坐標(biāo)軸原點(diǎn)）繪制時(shí)域仿真結(jié)果圖。

進(jìn)一步，圖3 對(duì)比了3 種方法在隨機(jī)選擇15 組測(cè)試樣本下所得的相對(duì)功角差曲線。由各條相對(duì)功角差動(dòng)態(tài)軌跡對(duì)比可知，在給定不確定性場(chǎng)景下，本文所提基于NARX-PCE 方法的含時(shí)序項(xiàng)代理模型能夠獲得與MCS 法基本相同的結(jié)果。同時(shí)，相比于基于PCE 方法的代理模型，本文所提代理模型的構(gòu)建不需要采集仿真時(shí)間內(nèi)所有時(shí)間斷面的樣本數(shù)據(jù)，僅需要仿真初始階段的時(shí)間斷面樣本數(shù)據(jù)（與建模過(guò)程中設(shè)置的TS有關(guān)）。因此，在后續(xù)無(wú)樣本數(shù)據(jù)為構(gòu)建代理模型提供訓(xùn)練的條件下（如圖3 中2.0～3.5 s 的區(qū)間內(nèi)），本文所提方法依然可以通過(guò)時(shí)序項(xiàng)遞推計(jì)算獲得功角差曲線變化的軌跡，而PCE 方法則無(wú)法給出該區(qū)間內(nèi)相對(duì)功角差的動(dòng)態(tài)軌跡。

圖3 不同測(cè)試樣本下相對(duì)功角差的動(dòng)態(tài)軌跡Fig.3 Dynamic trajectories of relative power angle difference with different test samples

考慮3.5 s 的時(shí)域仿真時(shí)間，利用MCS 法、PCE方法和NARX-PCE 方法獲得2 100 組（2 000 組訓(xùn)練集、100 組測(cè)試集）暫態(tài)相對(duì)功角差曲線的仿真計(jì)算總時(shí)間如表3 所示。

表3 不同方法計(jì)算時(shí)間對(duì)比Table 3 Comparison of computation time with different methods

由表3 可知，在用于電力系統(tǒng)時(shí)域仿真不確定性分析計(jì)算時(shí)，本文所提基于NARX-PCE 方法的含時(shí)序項(xiàng)代理模型相較于另外2 種方法具有較高的計(jì)算效率。盡管在代理模型構(gòu)建和給定新測(cè)試樣本計(jì)算過(guò)程中，NARX-PCE 方法相較于PCE 方法耗時(shí)較多，但由于NARX-PCE 方法不需要獲取所有時(shí)間斷面的樣本，在樣本集獲取過(guò)程中能夠節(jié)省大量仿真時(shí)間，從而算法總的計(jì)算時(shí)間顯著減少。值得說(shuō)明的是，隨著輸入隨機(jī)變量維度的增加，PCE方法所需樣本的數(shù)量會(huì)急劇增大，這將進(jìn)一步使得樣本獲取階段的耗時(shí)增加；而由于所提NARX-PCE方法利用了LAR 策略選擇代理模型中的關(guān)鍵基函數(shù)，降低了代理模型待求系數(shù)求解過(guò)程對(duì)樣本數(shù)量的依賴，可有效應(yīng)對(duì)高維輸入隨機(jī)變量引起的建模過(guò)程耗時(shí)增加的問(wèn)題，從而保證電力系統(tǒng)時(shí)域仿真不確定性分析過(guò)程的高效性。

3.1.2 時(shí)域仿真不確定性分析

如圖4（a）和圖4（b）所示，利用本文所提基于NARX-PCE 方法的含時(shí)序項(xiàng)代理模型分別獲得10 000 個(gè)隨機(jī)輸入場(chǎng)景下相對(duì)功角差δ1-2和節(jié)點(diǎn)3電壓V3的動(dòng)態(tài)軌跡。

圖4 基于NARX-PCE 方法的含時(shí)序項(xiàng)代理模型計(jì)算得到的電壓和功角差的動(dòng)態(tài)軌跡Fig.4 Dynamic trajectories of voltage and power angle difference obtained by calculation with NARX-PCE method based surrogate model containing time series terms

考慮G1 和G3 注入功率的隨機(jī)波動(dòng)，10 000 個(gè)隨機(jī)輸入場(chǎng)景下δ1-2保持在-10°～8°的區(qū)間內(nèi)波動(dòng)，并隨著仿真時(shí)間的推移逐漸趨于穩(wěn)定。同時(shí)，暫態(tài)電壓V3超過(guò)1.05 p.u.的概率為15.4%。通過(guò)上述計(jì)算結(jié)果可知，本文所提方法構(gòu)建的含時(shí)序項(xiàng)代理模型能夠代替原時(shí)域仿真計(jì)算模型，避免了大量隨機(jī)場(chǎng)景下重復(fù)的數(shù)值積分計(jì)算，顯著提升了計(jì)算效率，可用于實(shí)現(xiàn)電力系統(tǒng)時(shí)域仿真不確定性快速、準(zhǔn)確的分析。

3.2 IEEE 118 節(jié)點(diǎn)系統(tǒng)仿真結(jié)果

IEEE 118 節(jié)點(diǎn)系統(tǒng)接線圖如附錄D 圖D2 所示，發(fā)電機(jī)模型采用PSS@E 中的高階模型，部分模型參數(shù)可參見(jiàn)文獻(xiàn)［32］。通過(guò)不同故障場(chǎng)景的設(shè)置，給出了暫態(tài)穩(wěn)定、暫態(tài)不穩(wěn)定和暫態(tài)臨界穩(wěn)定3 種場(chǎng)景下的分析結(jié)果。

3.2.1 場(chǎng)景1:暫態(tài)穩(wěn)定

在母線54 與母線49 之間設(shè)置三相接地故障，故障開(kāi)始時(shí)間為t=0.2 s，故障持續(xù)時(shí)間為0.1 s，即t=0.3 s 故障清除。系統(tǒng)穩(wěn)態(tài)運(yùn)行下，以母線10、12、25、26 處的發(fā)電機(jī)有功功率P10、P12、P25、P26以及母線11、45、74、103 處的負(fù)荷L11、L45、L74、L103作為輸入隨機(jī)變量，且輸入隨機(jī)變量滿足正態(tài)分布，其均值即為發(fā)電機(jī)出力或負(fù)荷的基準(zhǔn)值，方差為各輸入隨機(jī)變量均值的1%。需要說(shuō)明的是，本文3.2 節(jié)中時(shí)刻原點(diǎn)即為仿真初始時(shí)刻。

選取母線54 處發(fā)電機(jī)的轉(zhuǎn)子角δ54作為輸出變量，利用2 000 組樣本構(gòu)建基于NARX-PCE 方法的含時(shí)序項(xiàng)代理模型，并進(jìn)行不確定性場(chǎng)景下的電力系統(tǒng)時(shí)域仿真計(jì)算，所得轉(zhuǎn)子角δ54的均值和方差隨時(shí)間動(dòng)態(tài)變化的軌跡分別如附錄E 圖E1（a）和圖E1（b）所示。模型構(gòu)建過(guò)程中最大時(shí)延參數(shù)tY=10，LOO 誤差值為7.206 3×10-6。

由附錄E 圖E1 中的仿真結(jié)果可知，在當(dāng)前故障場(chǎng)景下，考慮發(fā)電機(jī)隨機(jī)出力和負(fù)荷隨機(jī)波動(dòng)，轉(zhuǎn)子角δ54的均值逐步趨于穩(wěn)定，而通過(guò)功角差的計(jì)算可知系統(tǒng)處于穩(wěn)定狀態(tài)。同時(shí)，本文所提方法與MCS法的結(jié)果基本一致，表明本文所提方法能夠適用于含多維隨機(jī)變量（發(fā)電機(jī)隨機(jī)出力和負(fù)荷波動(dòng)）的電力系統(tǒng)時(shí)域仿真不確定性分析。

3.2.2 場(chǎng)景2:暫態(tài)不穩(wěn)定

在母線54 與母線49 之間設(shè)置三相接地故障，故障開(kāi)始時(shí)間為t=0.2 s，且故障一直存在。輸入隨機(jī)變量的選取與3.3.1 節(jié)保持一致。以發(fā)電機(jī)54 的轉(zhuǎn)子角δ54作為輸出變量，利用2 000 組樣本構(gòu)建基于NARX-PCE 方法的含時(shí)序項(xiàng)代理模型，并進(jìn)行不確定性場(chǎng)景下的電力系統(tǒng)時(shí)域仿真計(jì)算，所得轉(zhuǎn)子角δ54的均值和方差隨時(shí)間動(dòng)態(tài)變化的軌跡分別如附錄E 圖E2（a）和圖E2（b）所示。模型構(gòu)建過(guò)程中最大時(shí)延參數(shù)tY=10，LOO 誤差值為8.315 1×10-5。由圖中結(jié)果可知，在當(dāng)前故障場(chǎng)景下，轉(zhuǎn)子角δ54的均值呈現(xiàn)逐步增大的趨勢(shì)，通過(guò)計(jì)算發(fā)電機(jī)功角差可知系統(tǒng)處于不穩(wěn)定狀態(tài)，此時(shí)本文所提基于NARX-PCE 方法的含時(shí)序項(xiàng)代理模型也能夠得到與MCS 法基本一致的結(jié)果（方差結(jié)果的相對(duì)誤差在0.5%以內(nèi)），進(jìn)一步驗(yàn)證了本文所提方法的有效性。

3.2.3 場(chǎng)景3:暫態(tài)臨界穩(wěn)定

在母線100 與母線103 之間設(shè)置三相接地故障，故障開(kāi)始時(shí)間為t=0.2 s，故障持續(xù)時(shí)間為0.1 s，即t=0.3 s 故障清除。以母線10、12、25、26 處的發(fā)電機(jī)有功功率P10、P12、P25、P26為輸入隨機(jī)變量，且輸入隨機(jī)變量滿足正態(tài)分布，其均值即為發(fā)電機(jī)出力的基準(zhǔn)值，方差為各輸入隨機(jī)變量均值的10%。以發(fā)電機(jī)111 的轉(zhuǎn)子角δ111作為輸出變量，利用2 000 組樣本構(gòu)建基于NARX-PCE 方法的含時(shí)序項(xiàng)代理模型，并進(jìn)行不確定性場(chǎng)景下的電力系統(tǒng)時(shí)域仿真計(jì)算。隨機(jī)選取計(jì)算結(jié)果中的20 條轉(zhuǎn)子角δ111動(dòng)態(tài)軌跡，所得結(jié)果如圖5 所示。在當(dāng)前故障場(chǎng)景下，系統(tǒng)受發(fā)電機(jī)出力影響會(huì)同時(shí)存在暫態(tài)穩(wěn)定和不穩(wěn)定的情況，利用本文所提出的基于NARX-PCE 方法的含時(shí)序項(xiàng)代理模型仍舊可獲得不確定性場(chǎng)景下準(zhǔn)確的時(shí)域仿真結(jié)果（與MCS 方法結(jié)果基本一致），表明本文所提方法能夠應(yīng)用于暫態(tài)穩(wěn)定與不穩(wěn)定情況同時(shí)存在的復(fù)雜不確定性場(chǎng)景。

圖5 不同測(cè)試樣本下轉(zhuǎn)子角的動(dòng)態(tài)軌跡Fig.5 Dynamic trajectories of rotor angle with different test samples

4 結(jié)語(yǔ)

計(jì)及電力系統(tǒng)內(nèi)不確定性因素的影響，本文提出了一種基于含時(shí)序項(xiàng)代理模型的電力系統(tǒng)時(shí)域仿真不確定性分析方法，并將其應(yīng)用于不確定性環(huán)境下IEEE 9 節(jié)點(diǎn)和118 節(jié)點(diǎn)系統(tǒng)的暫態(tài)分析，所得結(jié)論如下:

1）利用NARX-PCE 方法所構(gòu)建的含時(shí)序項(xiàng)代理模型兼顧了系統(tǒng)隨時(shí)間變化的動(dòng)態(tài)特性和不確定性因素影響下的隨機(jī)特性，能夠適用于不確定性環(huán)境下電力系統(tǒng)時(shí)域仿真不確定性分析問(wèn)題；

2）相較于傳統(tǒng)MCS 法和PCE 方法，本文所提出的基于含時(shí)序項(xiàng)代理模型的電力系統(tǒng)時(shí)域仿真不確定性分析方法在保證計(jì)算精度的同時(shí)，降低了對(duì)大量時(shí)間斷面采樣過(guò)程的依賴，提高了不確定性環(huán)境下電力系統(tǒng)暫態(tài)分析的效率。

下一步工作將研究NARX-PCE 方法在電力系統(tǒng)時(shí)域仿真靈敏度分析中的應(yīng)用，進(jìn)一步辨識(shí)影響電力系統(tǒng)暫態(tài)穩(wěn)定的關(guān)鍵不確定性因素。

附錄見(jiàn)本刊網(wǎng)絡(luò)版（http：//www.aeps-info.com/aeps/ch/index.aspx），掃英文摘要后二維碼可以閱讀網(wǎng)絡(luò)全文。