一種用于矩陣指數(shù)的暫態(tài)仿真多時(shí)間尺度輸出方法
【專利摘要】一種用于矩陣指數(shù)的暫態(tài)仿真多時(shí)間尺度輸出方法，首先確定待仿真電力系統(tǒng)的結(jié)構(gòu)與參數(shù)，得到微分方程形式的暫態(tài)模型，初始化并啟動(dòng)仿真程序；在每一個(gè)仿真步長內(nèi)采用矩陣指數(shù)方法形成狀態(tài)矩陣，并確定本步長可用的連續(xù)平方環(huán)節(jié)個(gè)數(shù)；計(jì)算矩陣指數(shù)函數(shù)的Padé逼近矩陣，并利用這一矩陣的連續(xù)平方運(yùn)算和矩陣向量乘法運(yùn)算得到離散時(shí)間節(jié)點(diǎn)之間多個(gè)中間節(jié)點(diǎn)上的狀態(tài)向量的近似值，仿真推進(jìn)一個(gè)步長；依此迭代進(jìn)行，直到仿真結(jié)束。本發(fā)明實(shí)現(xiàn)了矩陣指數(shù)積分方法的仿真計(jì)算步長與結(jié)果輸出步長之間的解耦，使矩陣指數(shù)積分方法良好的數(shù)值精度和剛性處理能力能得到充分的發(fā)揮，解決了矩陣指數(shù)積分方法在多時(shí)間尺度特性的電力系統(tǒng)暫態(tài)仿真領(lǐng)域的適用性問題。
【專利說明】一種用于矩陣指數(shù)的暫態(tài)仿真多時(shí)間尺度輸出方法

【技術(shù)領(lǐng)域】
[0001] 本發(fā)明涉及一種暫態(tài)仿真的稠密輸出方法。特別是涉及一種用于矩陣指數(shù)的暫態(tài) 仿真多時(shí)間尺度輸出方法。

【背景技術(shù)】
[0002] 電力系統(tǒng)是由發(fā)電、輸電、配電、用電的大量設(shè)備及其相關(guān)輔助系統(tǒng)組成的復(fù)雜系 統(tǒng)，涉及電磁暫態(tài)、機(jī)械暫態(tài)、電化學(xué)過程、熱力學(xué)過程等多種動(dòng)態(tài)過程，有顯著的多時(shí)間尺 度特征。
[0003] 在傳統(tǒng)電力系統(tǒng)數(shù)字仿真領(lǐng)域，針對電力系統(tǒng)動(dòng)態(tài)過程的不同時(shí)間尺度分別發(fā)展 出了電磁暫態(tài)仿真和機(jī)電暫態(tài)仿真兩種建模仿真方法，二者從數(shù)學(xué)模型到仿真算法都有不 同的特征。其中，電磁暫態(tài)仿真主要反映系統(tǒng)中電場與磁場相互耦合影響產(chǎn)生的電氣量變 化過程，采用詳細(xì)的動(dòng)態(tài)建模和微秒級(jí)仿真步長，得到從工頻到幾十kHz頻譜范圍內(nèi)的三 相電壓電流瞬時(shí)值波形；機(jī)電暫態(tài)仿真則主要關(guān)注故障、切機(jī)切負(fù)荷等大擾動(dòng)發(fā)生后，發(fā)電 機(jī)轉(zhuǎn)子功角搖擺等機(jī)械變化過程，對電力網(wǎng)絡(luò)及與之相連的機(jī)組繞組采用交流穩(wěn)態(tài)建模， 模型較電磁暫態(tài)仿真進(jìn)行了簡化，通常采用毫秒級(jí)仿真步長，能夠捕捉工頻以下的系統(tǒng)動(dòng) 態(tài)過程。
[0004] 近年來，我國電力系統(tǒng)實(shí)際運(yùn)行出現(xiàn)的一些仿真場景，為傳統(tǒng)電力系統(tǒng)數(shù)字仿真 提出了新的問題。輸電系統(tǒng)方面，隨著特高壓技術(shù)和直流輸電技術(shù)的發(fā)展與應(yīng)用，交直流互 聯(lián)電網(wǎng)規(guī)模逐漸擴(kuò)大；配電系統(tǒng)方面，可再生電源和電力電子裝置的大量接入，也使得傳統(tǒng) 的無源配網(wǎng)成為有源配網(wǎng)。這些運(yùn)行工況的復(fù)雜暫態(tài)特性，客觀上需要使用建模更詳細(xì)的 電磁暫態(tài)仿真程序進(jìn)行研究分析，然而，較大的系統(tǒng)規(guī)模和較長的仿真時(shí)間無疑給電磁暫 態(tài)仿真提出了新的挑戰(zhàn)，需要結(jié)合問題特性從仿真算法層面提出針對性的改進(jìn)。
[0005] 電力系統(tǒng)電磁暫態(tài)仿真本質(zhì)上可歸結(jié)為對動(dòng)力學(xué)系統(tǒng)時(shí)域響應(yīng)的求取，它包括系 統(tǒng)本身的數(shù)學(xué)模型和與之相適應(yīng)的數(shù)值算法。
[0006] 當(dāng)前，電力系統(tǒng)電磁暫態(tài)仿真基本框架可分為兩類，包括節(jié)點(diǎn)分析法（Nodal Analysis)和狀態(tài)變量分析法（State-Variable Analysis)?；诠?jié)點(diǎn)分析框架的電磁暫 態(tài)仿真方法可概括為先采用某種數(shù)值積分方法（通常為梯形積分法）將系統(tǒng)中動(dòng)態(tài)元件的 特性方程差分化，得到等效的計(jì)算電導(dǎo)與歷史項(xiàng)電流源并聯(lián)形式的諾頓等效電路，對其求 解即可得到系統(tǒng)中各節(jié)點(diǎn)電壓的瞬時(shí)值。節(jié)點(diǎn)分析法廣泛應(yīng)用于EMTP、PSCAD/EMTDC等專 業(yè)的電力系統(tǒng)電磁暫態(tài)仿真程序中，工程上也稱基于節(jié)點(diǎn)分析框架的電磁暫態(tài)仿真工具為 EMTP類程序。節(jié)點(diǎn)分析法的主要優(yōu)勢體現(xiàn)在程序?qū)崿F(xiàn)難度和仿真計(jì)算效率方面，但由于節(jié) 點(diǎn)電導(dǎo)方程本身已將數(shù)值積分方法與系統(tǒng)模型融為一體，導(dǎo)致EMTP類程序在求解算法選 擇方面缺乏靈活性與開放性，同時(shí)節(jié)點(diǎn)電導(dǎo)方程也不能給出系統(tǒng)本身的特征信息。
[0007] 狀態(tài)變量分析法屬于一般性建模方法（general purpose modeling)，不僅適于 電路與電力系統(tǒng)仿真，同樣也適于其它工程領(lǐng)域的動(dòng)力學(xué)系統(tǒng)的建模與仿真。Matlab/ SimPowerSystems軟件是狀態(tài)變量分析框架下暫態(tài)仿真程序的典型代表。與節(jié)點(diǎn)分析框架 相比，狀態(tài)方程在模型的計(jì)算求解方面具有高度的開放性和靈活性，可方便地選擇與問題 相適應(yīng)的數(shù)值積分方法，同時(shí)能夠提供關(guān)于系統(tǒng)各種特征的豐富信息（如系統(tǒng)的特征值）， 進(jìn)而能夠從全局角度了解系統(tǒng)的動(dòng)態(tài)特性，為各種快速、準(zhǔn)確、高效的仿真算法的開發(fā)與測 試工作提供了便利條件。
[0008] 應(yīng)用狀態(tài)變量分析的基礎(chǔ)是形成狀態(tài)方程描述的電力系統(tǒng)暫態(tài)仿真模型。改進(jìn)節(jié) 點(diǎn)法Modified Nodal Analysis(MNA)通過KCUKVL等約束關(guān)系以及元件伏安特性構(gòu)造得到 MNA模型，再經(jīng)過一定的正規(guī)化處理（regularization)即可得到狀態(tài)方程模型；也可采用 一般支路類方法，如Automated State Model Generator (ASMG)方法直接構(gòu)造得到?；?這些方法得到的電力系統(tǒng)模型能夠很容易的與本發(fā)明所提出的算法進(jìn)行接口。
[0009] 在數(shù)值算法方面，傳統(tǒng)數(shù)值積分方法可分為顯式和隱式兩類，不同積分方法所具 有的數(shù)值穩(wěn)定性和數(shù)值精度各不相同。一般來說，隱式方法處理仿真模型中剛性特征的能 力較強(qiáng)。電力系統(tǒng)由于動(dòng)態(tài)過程時(shí)間尺度差異較大，系統(tǒng)模型表現(xiàn)出一定剛性，這使得主流 電磁暫態(tài)軟件EMTP類程序采用隱式方法以保證數(shù)值穩(wěn)定性。從計(jì)算開銷方面來看，隱式方 法在每一時(shí)步內(nèi)需求解線性方程組，極大限制了其在大規(guī)模系統(tǒng)的應(yīng)用能力。與之相對的， 傳統(tǒng)顯式方法無需迭代，在每一時(shí)步內(nèi)的運(yùn)算量較小，但其有限的數(shù)值穩(wěn)定域使得仿真步 長受到約束，綜合來看對剛性系統(tǒng)的仿真性能不佳。
[0010] 矩陣指數(shù)積分方法（Exponential Integrator)是近年來從應(yīng)用數(shù)學(xué)領(lǐng)域發(fā)端的 一種數(shù)值積分方法。它使用矩陣指數(shù)算子eM精確描述動(dòng)態(tài)系統(tǒng)的線性變化規(guī)律，可以準(zhǔn)確 求解形如
[0011] X'（t) = Ax (t)+Bu ⑴
[0012] 的線性動(dòng)力學(xué)系統(tǒng)，并具有計(jì)算效率高、剛性處理能力強(qiáng)等特點(diǎn)。矩陣指數(shù)積分方 法已經(jīng)在諸如應(yīng)用物理、化學(xué)工程等領(lǐng)域得到一定應(yīng)用。然而，對于隱式梯形法主導(dǎo)的電力 系統(tǒng)仿真領(lǐng)域，矩陣指數(shù)積分方法的單步計(jì)算量較大，在相同步長條件下，盡管矩陣指數(shù)積 分方法能提供更為精確的仿真結(jié)果，但仿真軟件使用者更為關(guān)心的計(jì)算時(shí)長問題反而更為 惡化；如果按照得到與梯形法相同誤差水平的仿真結(jié)果為準(zhǔn)則設(shè)置仿真步長，計(jì)算時(shí)間問 題雖然得到解決，但此時(shí)仿真步長較大，往往無法提供足夠的采樣頻率，可能遺漏系統(tǒng)動(dòng)態(tài) 過程中的快動(dòng)態(tài)細(xì)節(jié)。這一兩難的處境極大限制了矩陣指數(shù)積分方法在電力系統(tǒng)數(shù)字仿真 領(lǐng)域的應(yīng)用，是矩陣指數(shù)暫態(tài)仿真急需解決的問題。


【發(fā)明內(nèi)容】

[0013] 本發(fā)明所要解決的技術(shù)問題是，提供一種用于矩陣指數(shù)的暫態(tài)仿真多時(shí)間尺度輸 出方法，能夠在每一個(gè)時(shí)步內(nèi)，獲得具有良好精度的依用戶設(shè)置的若干稠密輸出點(diǎn)，使得仿 真可以在采用較大仿真步長計(jì)算的同時(shí)，得到較小步長的輸出結(jié)果，從而覆蓋待研究電力 系統(tǒng)在多個(gè)時(shí)間尺度上的動(dòng)態(tài)過程。
[0014] 本發(fā)明所采用的技術(shù)方案是：一種用于矩陣指數(shù)的暫態(tài)仿真多時(shí)間尺度輸出方 法，包括如下步驟：
[0015] 1)根據(jù)仿真算例輸入文件，確定待仿真電力系統(tǒng)的結(jié)構(gòu)與參數(shù)，得到微分方程形 式的暫態(tài)模型x' = f (X)，X為狀態(tài)變量；
[0016] 2)設(shè)置仿真結(jié)束時(shí)間T，設(shè)置仿真計(jì)算時(shí)刻tn = 0,進(jìn)行仿真初始化計(jì)算；
[0017] 3)記tn時(shí)刻的仿真步長為h，狀態(tài)向量為xn，設(shè)定需要進(jìn)行稠密輸出的平方環(huán)節(jié) 的個(gè)數(shù)Ns，Ns是正整數(shù)；
[0018] 4)采用矩陣指數(shù)方法形成時(shí)步[tn，tn+h]所需的狀態(tài)矩陣A n ;
[0019] 5)確定[tn，tn+h]時(shí)步內(nèi)，矩陣指數(shù)計(jì)算所用的連續(xù)平方環(huán)節(jié)個(gè)數(shù)s，若s〈N s，臨時(shí) 將s的值賦值給Ns，否則Ns維持不變，記本時(shí)步的輸出為X，X是2Ns個(gè)列向量組成的矩陣；
[0020] 6)計(jì)算Pad6逼近矩陣F = rkm(h/2sAn)，'是矩陣指數(shù)函數(shù)的Pad6逼近函數(shù)；
[0021] 7)進(jìn)行（S-Ns)次F的連續(xù)平方，并將結(jié)果矩陣賦值給F，將F和Xn的矩陣向量乘 積賦值給X的第一列，作為t n+h/2Ns時(shí)刻狀態(tài)向量的值，如果Ns = 1，跳至步驟11)，否則設(shè) 定循環(huán)變量i = 1，并進(jìn)入下一步驟；
[0022] 8) F = F2,將F和Xn的矩陣向量乘積賦值給X的第21列，作為tn+h/2 Ns_i時(shí)亥Ij狀態(tài) 向量的值；
[0023] 9)記X的1至2^1列組成的矩陣為Xi，將F和Xi的矩陣乘積賦值給X的第2kl 至 2i+1-l 列，作為{^+0+1)/2?, tn+(2i+2)/2Nsh，…，tn+(2i+1-l)/2 Nsh}這 2Ll 個(gè)時(shí)刻的狀 態(tài)向量的值；
[0024] 10)更新i = i+1，如果i = Ns，結(jié)束本時(shí)步的仿真計(jì)算和輸出，將X的值寫入輸出 文件后進(jìn)入步驟11)，否則返回步驟8);
[0025] 11)更新tn = tn+h，仿真向前推進(jìn)一個(gè)步長，如果tn小于仿真結(jié)束時(shí)間T，返回步 驟3)繼續(xù)進(jìn)行仿真計(jì)算，否則仿真結(jié)束。
[0026] 步驟4)所述的采用矩陣指數(shù)方法形成時(shí)步[tn，tn+h]所需的狀態(tài)矩陣A n的具體 形式是：
[0027]

【權(quán)利要求】
1. 一種用于矩陣指數(shù)的暫態(tài)仿真多時(shí)間尺度輸出方法，其特征在于，包括如下步驟： 1) 根據(jù)仿真算例輸入文件，確定待仿真電力系統(tǒng)的結(jié)構(gòu)與參數(shù)，得到微分方程形式的 暫態(tài)模型X' = f (X)，X為狀態(tài)變量； 2) 設(shè)置仿真結(jié)束時(shí)間T，設(shè)置仿真計(jì)算時(shí)刻tn = 0,進(jìn)行仿真初始化計(jì)算； 3) 記tn時(shí)刻的仿真步長為h，狀態(tài)向量為Xn，設(shè)定需要進(jìn)行稠密輸出的平方環(huán)節(jié)的個(gè) 數(shù)Ns，Ns是正整數(shù)； 4) 采用矩陣指數(shù)方法形成時(shí)步[tn，tn+h]所需的狀態(tài)矩陣An; 5) 確定[tn，tn+h]時(shí)步內(nèi)，矩陣指數(shù)計(jì)算所用的連續(xù)平方環(huán)節(jié)個(gè)數(shù)s，若s〈Ns，臨時(shí)將s 的值賦值給Ns，否則Ns維持不變，記本時(shí)步的輸出為X，X是2Ns個(gè)列向量組成的矩陣； 6) 計(jì)算Pad6逼近矩陣F = rkm(h/2sAn)，rkm是矩陣指數(shù)函數(shù)的Pad6逼近函數(shù)； 7) 進(jìn)行（S-Ns)次F的連續(xù)平方，并將結(jié)果矩陣賦值給F，將？和\的矩陣向量乘積賦 值給X的第一列，作為t n+h/2Ns時(shí)刻狀態(tài)向量的值，如果Ns = 1，跳至步驟11)，否則設(shè)定循 環(huán)變量i = 1，并進(jìn)入下一步驟； 8. F = F2,將F和Xn的矩陣向量乘積賦值給X的第21列，作為tn+h/2 Ns4時(shí)刻狀態(tài)向量 的值； 9) 記X的1至2^1列組成的矩陣為Xi，將F和Xi的矩陣乘積賦值給X的第2^1至 2 i+1_l 列，作為{tn+(2^1)/2?, tn+(2^2)/2?,…，tn+(2i+1-l)/2Nsh}這 2^1 個(gè)時(shí)刻的狀態(tài) 向量的值； 10) 更新i = i+Ι，如果i =Ns，結(jié)束本時(shí)步的仿真計(jì)算和輸出，將X的值寫入輸出文件 后進(jìn)入步驟11)，否則返回步驟8); 11) 更新tn = tn+h，仿真向前推進(jìn)一個(gè)步長，如果tn小于仿真結(jié)束時(shí)間T，返回步驟3) 繼續(xù)進(jìn)行仿真計(jì)算，否則仿真結(jié)束。
2. 權(quán)利要求1所述的一種用于矩陣指數(shù)的暫態(tài)仿真多時(shí)間尺度輸出方法，其特征在 于，步驟4)所述的采用矩陣指數(shù)方法形成時(shí)步[t n，tn+h]所需的狀態(tài)矩陣An的具體形式 是：
其中，是非線性函數(shù)f (X)在Xn取值時(shí)的Jacobian矩陣。
3. 權(quán)利要求1所述的一種用于矩陣指數(shù)的暫態(tài)仿真多時(shí)間尺度輸出方法，其特征在 于，步驟5)所述的確定[tn，t n+h]時(shí)步內(nèi)矩陣指數(shù)計(jì)算所用的連續(xù)平方環(huán)節(jié)個(gè)數(shù)s的具體 方法為： (1) 計(jì)算仿真步長h與狀態(tài)矩陣An的乘積矩陣的列范數(shù)，記為Nm ; (2) 如果NmS Θ，其中Θ是通過數(shù)值分析得到的閾值，則本時(shí)步?jīng)]有能夠進(jìn)行稠密輸 出的平方環(huán)節(jié)，取s = 0,否則進(jìn)入下一步驟； ⑶選取小數(shù)M和整數(shù)E，使得0·5彡M〈l，MX2E = Nm/θ; (4)如果 M = 0· 5,取 s = E-1，否則取 s = E。
【文檔編號(hào)】G06F17/50GK104376158SQ201410619247
【公開日】2015年2月25日 申請日期:2014年11月5日 優(yōu)先權(quán)日:2014年11月5日
【發(fā)明者】王成山, 富曉鵬, 李鵬, 丁承第, 于浩, 于力, 郭曉斌, 許愛東, 董旭柱, 吳爭榮, 田兵, 王建邦, 魏文瀟 申請人:天津大學(xué), 中國南方電網(wǎng)有限責(zé)任公司電網(wǎng)技術(shù)研究中心, 南方電網(wǎng)科學(xué)研究院有限責(zé)任公司