狀態(tài)估計實現(xiàn))
電力系統(tǒng)狀態(tài)估計這個話題老早以前是SCADA的天下調(diào)度員靠RTU傳來的有功、無功和幅值再用非線性加權(quán)最小二乘去迭代一套下來動不動幾十次迭代碰上壞數(shù)據(jù)還得來回排查。這幾年P(guān)MU相量測量單元鋪開之后局面變了不少——它直接給出帶時標的電壓相量幅值有了相角也有了狀態(tài)估計從非線性一下子被拉回到線性模型算起來清爽太多。今天想聊的這個項目就是用MATLAB搭一套基于PMU量測的電力系統(tǒng)狀態(tài)估計標題叫《基于matlab PMU相量測量單元電力系統(tǒng)狀態(tài)估計》帶源碼編號14925期。我按自己做同類項目的經(jīng)驗把整體設(shè)計思路、數(shù)學模型、MATLAB實現(xiàn)細節(jié)和踩過的坑全部拆開講一遍給正在做電力系統(tǒng)課設(shè)、畢設(shè)或者入門廣域監(jiān)測系統(tǒng)的朋友一個可以直接落地的參考。這套東西能解決什么問題簡單說在系統(tǒng)里裝了一批PMU之后如何在已知網(wǎng)絡(luò)拓撲、線路參數(shù)和部分量測數(shù)據(jù)的情況下把全網(wǎng)各母線電壓的幅值和相角給估出來。有了這個東西后續(xù)的潮流追蹤、擾動定位、靜態(tài)穩(wěn)定分析才有基礎(chǔ)數(shù)據(jù)。適合誰看電氣工程相關(guān)專業(yè)的學生、剛接觸WAMS廣域測量系統(tǒng)的工程師以及想把MATLAB數(shù)值計算能力和電力系統(tǒng)分析結(jié)合起來練手的開發(fā)者。下面按項目從設(shè)計到代碼再到排錯的順序從頭捋一遍。1. 項目整體設(shè)計與思路拆解1.1 核心需求為什么有了PMU狀態(tài)估計變簡單了傳統(tǒng)狀態(tài)估計用的是SCADA量測量測類型主要是節(jié)點注入有功、無功、支路潮流和電壓幅值。注意這里面沒有相角量測因為SCADA根本沒有統(tǒng)一時標去測相角。于是狀態(tài)量和量測量之間是非線性關(guān)系比如支路有功潮流公式是P_ij V_i2G_ij - V_iV_j(G_ij cosθ_ij B_ij sinθ_ij)求解必須靠高斯-牛頓迭代每次迭代都要重新算雅可比矩陣計算量大而且在初值離真值遠的時候還可能發(fā)散。PMU把這個痛點直接戳掉了。它能以GPS/北斗授時同步采樣輸出的就是帶絕對時標的電壓相量幅值V_i和相角θ_i是直接測出來的。因此狀態(tài)量和量測量之間變成了線性關(guān)系z Hx e。這種情況下狀態(tài)估計本質(zhì)上就是一個線性加權(quán)最小二乘問題不用迭代一次矩陣運算就出結(jié)果穩(wěn)定性和速度都有質(zhì)的提升。打個比方SCADA就像你只知道一個人的身高卻要推斷他的站姿PMU則是直接給你一張帶角度的照片身高、傾斜角度全都標好了。不是不需要估計但估計的工作量小了很多。1.2 方案選型為什么用線性WLS而不是卡爾曼濾波這個項目走的是加權(quán)最小二乘WLS路線而且是線性版本。選它有幾個理由第一靜態(tài)狀態(tài)估計是基礎(chǔ)。WLS是經(jīng)典框架理論成熟、代碼易懂、結(jié)果可以直接跟真值做誤差對比。對課設(shè)和畢設(shè)來說這個路線做出來有理有據(jù)答辯也好講。第二線性WLS不需要迭代矩陣的反復(fù)分解。傳統(tǒng)WLS每輪迭代都要做一次LU分解而線性模型只需要一次求解對MATLAB這種解釋型語言來說性能友好很多。第三卡爾曼濾波雖然能處理動態(tài)過程但它需要模型噪聲和量測噪聲的先驗協(xié)方差調(diào)參難度大。這里先做靜態(tài)版本后續(xù)要擴展成動態(tài)估計可以在此框架上加狀態(tài)轉(zhuǎn)移方程那屬于進階玩法。結(jié)論以線性WLS為骨架把PMU量測方程、權(quán)重矩陣、可觀測性分析、壞數(shù)據(jù)檢測這幾個模塊串起來是性價比最高的方案。1.3 數(shù)據(jù)流和模塊劃分整個項目的核心流程可以切分成五個模塊數(shù)據(jù)輸入系統(tǒng)節(jié)點數(shù)、支路參數(shù)、PMU安裝位置及量測值?？捎^測性分析檢查量測是否足夠讓H矩陣列滿秩。核心估計構(gòu)建z、H、R求解WLS狀態(tài)估計。壞數(shù)據(jù)檢測計算殘差做檢測與剔除。結(jié)果輸出展示估計值、誤差指標可視化對比。這個劃分也直接映射到MATLAB的代碼結(jié)構(gòu)上每個功能對應(yīng)一個函數(shù)main腳本負責串聯(lián)。后面第三章會給出各模塊的代碼骨架和關(guān)鍵參數(shù)選取邏輯。2. 核心模型與算法拆解2.1 PMU量測方程與狀態(tài)向量的數(shù)學形式先定義狀態(tài)向量。對于節(jié)點數(shù)為n的電力系統(tǒng)若所有節(jié)點都裝了PMU每個節(jié)點的狀態(tài)是電壓幅值V_i和相角θ_i那么理論上的狀態(tài)量總數(shù)是2n。但因為全網(wǎng)相角需要一個參考基準通常把某個參考節(jié)點的相角錨定為0實際上待估計的狀態(tài)數(shù)是2n-1。PMU的量測方程可以寫成z Hx e其中x是狀態(tài)向量e是量測噪聲向量。假設(shè)在母線i上裝了一臺PMU它能同時量測該母線的電壓幅值和相角還能量測與該母線相連的支路電流相量。以電壓量測為例對應(yīng)H矩陣的行是第i個幅值狀態(tài)位為1第i個相角狀態(tài)位為0對應(yīng)相角量測的行是第i個相角狀態(tài)位為1幅值狀態(tài)位為0。如果PMU還提供支路電流相量那需要基于線路的π型等值電路把電流相量表達式轉(zhuǎn)換為狀態(tài)量的線性組合這部分關(guān)鍵是要精確填H矩陣的系數(shù)稍有差錯估計結(jié)果就會偏掉。這里面有個細節(jié)很多人第一次做容易漏PMU給的是絕對相角以UTC為基準的全球同步相角而狀態(tài)估計里的狀態(tài)量是相對參考母線的相角。處理方式是把所有相角量測減去參考母線的相角量測再進入估計。否則H矩陣會多一列全零或者出現(xiàn)秩虧。2.2 加權(quán)矩陣R的選取為什么要按量測類型分權(quán)重WLS的目標函數(shù)是J (z-Hx)?R?1(z-Hx)R是量測誤差協(xié)方差矩陣。PMU的相量測量單元精度比傳統(tǒng)RTU高很多幅值誤差通常在0.1%量級相角誤差在0.01°~0.02°量級。但不同PMU通道、不同幅值和相角的方差差異還是存在所以R不能簡單設(shè)成單位陣。實踐中常用做法是查PMU的精度指標把幅值標準差和相角標準差換算成方差電壓幅值量測標準差約0.001~0.002 p.u.對應(yīng)R對角線元素約為1e-6~4e-6。電壓相角量測標準差約0.0002~0.0004 rad對應(yīng)R對角線元素約4e-8~1.6e-7。電流幅值和相角量測取決于CT/PT的精度等級通常比重會比電壓量測的方差大一些。權(quán)重的本質(zhì)是量測的信任度方差越小權(quán)重越大在求解時對結(jié)果的貢獻越大。這個道理和加權(quán)平均是一樣的實際項目里我習慣先用均勻權(quán)重跑一遍看殘差分布再用殘差方差反推R做一次迭代定權(quán)。這個小技巧能明顯改善結(jié)果。2.3 可觀測性分析為什么H矩陣必須滿秩線性狀態(tài)估計能解出唯一解的前提是量測方程個數(shù)不小于狀態(tài)數(shù)而且H矩陣列滿秩。列滿秩意味著每個狀態(tài)量都被足夠的獨立量測覆蓋不存在某個母線電壓怎么測都測不到的情況。在MATLAB里判斷很簡單計算秩rank(H)列滿秩的判據(jù)是rank(H)等于狀態(tài)數(shù)2n-1。同時還可以看條件數(shù)cond(H)條件數(shù)太大說明H矩陣近似病態(tài)哪怕滿秩數(shù)值上也可能解出漫天亂跳的結(jié)果。條件數(shù)控制在1e6以內(nèi)比較好超過這個量級就要警惕。這就引出一個部署問題PMU數(shù)量不夠怎么辦工程上常見的是PMU只裝在部分關(guān)鍵節(jié)點剩下的節(jié)點靠SCADA量測補齊。這種混合量測場景下模型重新變成非線性得用傳統(tǒng)WLS迭代。如果非要保持線性模型可以假設(shè)SCADA區(qū)域的狀態(tài)初始值已知那其實就不叫狀態(tài)估計了屬于擾動分析邏輯上要分清。2.4 壞數(shù)據(jù)檢測標準化殘差怎么用PMU數(shù)據(jù)也不是百分百干凈通信丟包、相量計算異常、GPS失步都會產(chǎn)生壞數(shù)據(jù)。線性模型下殘差r z - Hx_est理論上服從零均值高斯分布。采用的是基于標準化殘差的檢測r_i_normalized r_i / sqrt(R_ii * (I - H(H?R?1H)?1H?R?1)_ii)分子是第i個量測的殘差分母是殘差方差的開方。標準化之后r_i_normalized近似服從標準正態(tài)分布用閾值λ一般取3.0對應(yīng)99.7%置信度去卡超過閾值就判壞數(shù)據(jù)。有一個容易踩的坑多個壞數(shù)據(jù)同時存在時逐次剔除比一次性剔除更穩(wěn)。因為壞數(shù)據(jù)可能互相掩蓋殘差會被拉平單次殘差檢驗會漏掉。每次只剔除標準化殘差最大的那個量測重新做一遍估計再檢直到?jīng)]有超閾值的點為止。后面第四章會展開講這個問題的具體表現(xiàn)。3. MATLAB實操實現(xiàn)與核心環(huán)節(jié)3.1 數(shù)據(jù)準備我用IEEE 9節(jié)點系統(tǒng)作為測試床我復(fù)現(xiàn)這個項目時用的算例是IEEE 9節(jié)點系統(tǒng)節(jié)點數(shù)少、拓撲清楚、Matpower里有現(xiàn)成數(shù)據(jù)適合驗證算法。數(shù)據(jù)準備階段要明確幾樣?xùn)|西節(jié)點表9個節(jié)點包括基準電壓、類型PQ/PV/平衡。支路表每條支路的電阻、電抗、對地電納以及變壓器變比。PMU位置試驗時我假定節(jié)點1、3、6、9裝了PMU量測覆蓋這四點的電壓相量以及相連支路的電流相量。在MATLAB里我推薦用struct組織這些數(shù)據(jù)別用一堆散變量data.n 9; data.branch [ 1 4 0.0000 0.0576 0.0000 0; 4 5 0.0170 0.0920 0.1580 0; 5 6 0.0390 0.1700 0.3580 0; ... ]; data.pmu [1; 3; 6; 9]; data.z_meas YOUR_MEAS_VECTOR;如果手上沒有實測PMU數(shù)據(jù)可以用潮流計算結(jié)果作為真值再疊加高斯噪聲生成量測值。這個做法對驗證代碼正確性特別有用——因為你知道真值就能算誤差。3.2 核心求解函數(shù)從測量向量到狀態(tài)量這一段是代碼的樞紐。給出一個核心函數(shù)框架讀者可以直接改成自己的數(shù)據(jù)規(guī)模。function x_est pmu_wls_se(z, H, R) % z: 量測向量 m x 1 % H: 量測矩陣 m x (2n-1) % R: 量測誤差協(xié)方差矩陣 m x m G H * (R \ H); % 信息矩陣 b H * (R \ z); % 右端項 x_est G \ b; % 最小二乘解 end實際項目里H不是手工填的而是根據(jù)PMU位置和網(wǎng)絡(luò)拓撲動態(tài)生成。構(gòu)建H矩陣的邏輯分三步第一步建立狀態(tài)索引。每個節(jié)點分配兩個索引幅值索引idxV_i 2*(i-1)1相角索引idxTheta_i 2*(i-1)2。參考節(jié)點的相角索引要特殊處理要么在H中刪掉該列要么在x中固定為0并同步調(diào)整量測方程。第二步填電壓量測行。對于母線i的PMU幅值量測行在idxV_i位置填1相角量測行在idxTheta_i位置填1。相角量測如果用的是全局相角記得統(tǒng)一減去參考節(jié)點的全局相角后再進估計。第三步填支路電流量測行。以π型等值電路為準先算線路導(dǎo)納Y_ij G jB再根據(jù)電流相量I_ij Y_ij(V_i - V_j) jBsh/2 * V_i把實部和虛部對狀態(tài)量的偏導(dǎo)算出來填到對應(yīng)位置。這一步最容易寫錯強烈建議先用一個簡單兩節(jié)點系統(tǒng)驗證H矩陣的正確性。3.3 結(jié)果驗證誤差分析和殘差分析狀態(tài)估計算完不能直接交差得驗證。我的做法是把估計值跟潮流真值做對比計算每個節(jié)點的幅值誤差和相角誤差。通常用RMSE均方根誤差來評價整體精度公式是RMSE sqrt(mean((x_est - x_true).^2))從我的測試結(jié)果看在R矩陣設(shè)置合理、量測噪聲標準差符合PMU實際水平的前提下9節(jié)點系統(tǒng)的幅值估計誤差大約在1e-4 p.u.量級相角估計誤差大約在1e-3 rad量級。如果誤差偏大一個數(shù)量級以上優(yōu)先檢查H矩陣有沒有填錯其次是R矩陣是否給了不合理的權(quán)重。另外要畫殘差分布圖。把標準化殘差畫成條形圖直觀能看出有沒有異常量測。正常情況下殘差密布在±3之間且沒有明顯單點突出如果有突出點先別急著刪確認是數(shù)據(jù)問題還是H矩陣問題。3.4 整體腳本結(jié)構(gòu)與運行流程main腳本的結(jié)構(gòu)可以這樣安排% 第1步加載系統(tǒng)數(shù)據(jù)和PMU量測 system_data load_pmu_system(ieee9); % 第2步構(gòu)建H矩陣和R矩陣 [H, R, z, idx] build_linear_model(system_data); % 第3步可觀測性檢查 assert(rank(H) size(H,2), H矩陣秩虧系統(tǒng)不可觀測); % 第4步狀態(tài)估計 x_est pmu_wls_se(z, H, R); % 第5步壞數(shù)據(jù)檢測 [r_norm, bad_idx] bad_data_detect(z, H, R, x_est); % 第6步結(jié)果輸出 plot_result(x_est, system_data);這里每步調(diào)用一個函數(shù)函數(shù)內(nèi)部再細分可讀性和可維護性都比寫一個兩三百行的主腳本強得多。有一個小建議運行前用tic/toc記錄時間線性模型下9節(jié)點系統(tǒng)的計算時間應(yīng)該在毫秒級G矩陣是17x17MATLAB分解起來非?？臁Ｈ绻艹鰜硪脦酌牖究梢源_定有冗余循環(huán)或者H矩陣構(gòu)建邏輯低效需要排查。4. 常見問題與排查技巧實錄4.1 H矩陣奇異或條件數(shù)過大怎么辦這是我在項目里遇到最多的一個問題幾乎每個初做PMU狀態(tài)估計的人都會撞上一次。表現(xiàn)就是rank(H) size(H,2)或者cond(H)在1e12以上估計值亂跳?？赡茉蛴腥齻€。一是PMU覆蓋不足某些母線完全沒有量測覆蓋對應(yīng)H行全是零狀態(tài)不可觀。二是參考相角沒有處理干凈H矩陣里含有一列全零或者兩列線性相關(guān)。三是線路參數(shù)填錯導(dǎo)致支路電流量測對應(yīng)行與電壓量測行產(chǎn)生線性相關(guān)關(guān)系。排查順序第一步打印H矩陣的稀疏模式用spy(H)看哪列全零哪兩列是成比例關(guān)系。第二步逐個PMU檢查量測方程個數(shù)和類型確認覆蓋范圍。第三步用一個只有兩臺PMU的兩節(jié)點系統(tǒng)做單元測試H矩陣規(guī)模小一眼能看出問題。提示如果是可觀測性不足不要試圖在代碼層面打補丁。老老實實增加PMU量測或者把部分SCADA量測補進模型用混合量測的非線性WLS。強行求解只會得到數(shù)值上看起來正常、實際上完全沒意義的解。4.2 相角參考基準沖突PMU量測給出的是全球同步的絕對相角不同廠家的PMU在接入同一系統(tǒng)時由于GPS信號處理延遲的差異可能會有微小的角度偏移。如果直接把不同PMU的相角拿來拼成同一個z向量容易出現(xiàn)系統(tǒng)性偏差。我的處理辦法是在數(shù)據(jù)預(yù)處理階段先把所有PMU相角量測減去同一個參考PMU的相角量測得到相對相角序列再進入估計。這樣即使PMU本身有固定延遲誤差只要延遲在短時間內(nèi)穩(wěn)定相減后誤差會被抵消掉。這個操作在代碼里就是一行z_theta z_theta - z_theta_ref。4.3 壞數(shù)據(jù)檢測的誤檢與漏檢誤檢通常是因為R矩陣的方差設(shè)得太小量測噪聲本來沒那么高精度標準化之后殘差就會偏大超過閾值。漏檢則常見于兩個壞數(shù)據(jù)互相抵消的場景。舉個例子某條支路兩端的電流量測同時壞掉它們的殘差可能方向相反平均下來標準化殘差不大就會被漏掉。實操中我的建議是別只依賴一次殘差檢驗。做一個循環(huán)剔除法——每次只刪標準化殘差最大的那個量測重新估計后再檢。同時把檢驗閾值從3.0放寬到2.8多捕獲一些邊緣可疑點寧可多剔除一個可疑量測也不能放壞數(shù)據(jù)進門。當然剔除的量測數(shù)不能太多一般超過總量測數(shù)的5%就要回頭檢查是不是數(shù)據(jù)質(zhì)量整體不行。4.4 量測噪聲設(shè)置與實際不符測試階段想要模擬真實場景可以在潮流真值上加高斯白噪聲但噪聲標準差的選擇要有依據(jù)。我見過有人直接用randn加噪聲標準差設(shè)成0.1結(jié)果估計誤差大得離譜還以為是算法有問題實際上是噪聲水平跟PMU的真實指標差了百倍。PMU的幅值測量精度典型值在0.1%相角測量精度在0.01°左右對應(yīng)弧度約1.7e-4 rad。按這個量級設(shè)置噪聲估計結(jié)果才能反映算法本身的性能。做完之后可以統(tǒng)計殘差的標準差和設(shè)置的噪聲標準差做對比如果差太多說明H矩陣或者加權(quán)有問題。4.5 計算效率與內(nèi)存小技巧PMU量測點數(shù)一旦多了H矩陣和R矩陣的維度會漲得很快。比如IEEE 118節(jié)點系統(tǒng)全裝PMU狀態(tài)量就是235個量測可能有上千行。這時候如果H還是稠密矩陣求逆操作會越來越慢。兩個改進方向一是用稀疏矩陣存儲H和RMATLAB里直接sparse(H)二是用信息矩陣G HRH然后對G做Cholesky分解而不是直接求H的偽逆。這兩步能把計算時間下降兩個數(shù)量級。我在118節(jié)點系統(tǒng)上試過從幾十秒降到幾百毫秒效果非常明顯。另外一個細節(jié)R矩陣是對角陣R \ H這一步不要寫成inv(R) * H直接用左除效率更高數(shù)值也穩(wěn)定得多。MATLAB里對稀疏對角陣的除法有專門優(yōu)化一定要利用上。最后再分享一個擴展思路。這個靜態(tài)估計框架跑通之后如果還想往深做可以嘗試兩個方向一是把PMU量測數(shù)據(jù)按時間序列連續(xù)輸入加入狀態(tài)轉(zhuǎn)移模型升級成動態(tài)狀態(tài)估計這時就該上卡爾曼濾波了二是在現(xiàn)有框架中把H矩陣的構(gòu)建推廣到三相不平衡系統(tǒng)就能用來處理配電網(wǎng)狀態(tài)估計。我個人做下來最大的體會是這套模型的代碼骨架一旦搭好往各種方向擴展都很快關(guān)鍵是前期的數(shù)據(jù)結(jié)構(gòu)和H矩陣構(gòu)建邏輯要寫干凈別為省幾行代碼把后續(xù)的擴展性毀了。 基于MATLAB PMU相量測量單元的電力系統(tǒng)狀態(tài)估計實現(xiàn)
電力系統(tǒng)狀態(tài)估計這個話題老早以前是SCADA的天下調(diào)度員靠RTU傳來的有功、無功和幅值再用非線性加權(quán)最小二乘去迭代一套下來動不動幾十次迭代碰上壞數(shù)據(jù)還得來回排查。這幾年P(guān)MU相量測量單元鋪開之后局面變了不少——它直接給出帶時標的電壓相量幅值有了相角也有了狀態(tài)估計從非線性一下子被拉回到線性模型算起來清爽太多。今天想聊的這個項目就是用MATLAB搭一套基于PMU量測的電力系統(tǒng)狀態(tài)估計標題叫《基于matlab PMU相量測量單元電力系統(tǒng)狀態(tài)估計》帶源碼編號14925期。我按自己做同類項目的經(jīng)驗把整體設(shè)計思路、數(shù)學模型、MATLAB實現(xiàn)細節(jié)和踩過的坑全部拆開講一遍給正在做電力系統(tǒng)課設(shè)、畢設(shè)或者入門廣域監(jiān)測系統(tǒng)的朋友一個可以直接落地的參考。這套東西能解決什么問題簡單說在系統(tǒng)里裝了一批PMU之后如何在已知網(wǎng)絡(luò)拓撲、線路參數(shù)和部分量測數(shù)據(jù)的情況下把全網(wǎng)各母線電壓的幅值和相角給估出來。有了這個東西后續(xù)的潮流追蹤、擾動定位、靜態(tài)穩(wěn)定分析才有基礎(chǔ)數(shù)據(jù)。適合誰看電氣工程相關(guān)專業(yè)的學生、剛接觸WAMS廣域測量系統(tǒng)的工程師以及想把MATLAB數(shù)值計算能力和電力系統(tǒng)分析結(jié)合起來練手的開發(fā)者。下面按項目從設(shè)計到代碼再到排錯的順序從頭捋一遍。1. 項目整體設(shè)計與思路拆解1.1 核心需求為什么有了PMU狀態(tài)估計變簡單了傳統(tǒng)狀態(tài)估計用的是SCADA量測量測類型主要是節(jié)點注入有功、無功、支路潮流和電壓幅值。注意這里面沒有相角量測因為SCADA根本沒有統(tǒng)一時標去測相角。于是狀態(tài)量和量測量之間是非線性關(guān)系比如支路有功潮流公式是P_ij V_i2G_ij - V_iV_j(G_ij cosθ_ij B_ij sinθ_ij)求解必須靠高斯-牛頓迭代每次迭代都要重新算雅可比矩陣計算量大而且在初值離真值遠的時候還可能發(fā)散。PMU把這個痛點直接戳掉了。它能以GPS/北斗授時同步采樣輸出的就是帶絕對時標的電壓相量幅值V_i和相角θ_i是直接測出來的。因此狀態(tài)量和量測量之間變成了線性關(guān)系z Hx e。這種情況下狀態(tài)估計本質(zhì)上就是一個線性加權(quán)最小二乘問題不用迭代一次矩陣運算就出結(jié)果穩(wěn)定性和速度都有質(zhì)的提升。打個比方SCADA就像你只知道一個人的身高卻要推斷他的站姿PMU則是直接給你一張帶角度的照片身高、傾斜角度全都標好了。不是不需要估計但估計的工作量小了很多。1.2 方案選型為什么用線性WLS而不是卡爾曼濾波這個項目走的是加權(quán)最小二乘WLS路線而且是線性版本。選它有幾個理由第一靜態(tài)狀態(tài)估計是基礎(chǔ)。WLS是經(jīng)典框架理論成熟、代碼易懂、結(jié)果可以直接跟真值做誤差對比。對課設(shè)和畢設(shè)來說這個路線做出來有理有據(jù)答辯也好講。第二線性WLS不需要迭代矩陣的反復(fù)分解。傳統(tǒng)WLS每輪迭代都要做一次LU分解而線性模型只需要一次求解對MATLAB這種解釋型語言來說性能友好很多。第三卡爾曼濾波雖然能處理動態(tài)過程但它需要模型噪聲和量測噪聲的先驗協(xié)方差調(diào)參難度大。這里先做靜態(tài)版本后續(xù)要擴展成動態(tài)估計可以在此框架上加狀態(tài)轉(zhuǎn)移方程那屬于進階玩法。結(jié)論以線性WLS為骨架把PMU量測方程、權(quán)重矩陣、可觀測性分析、壞數(shù)據(jù)檢測這幾個模塊串起來是性價比最高的方案。1.3 數(shù)據(jù)流和模塊劃分整個項目的核心流程可以切分成五個模塊數(shù)據(jù)輸入系統(tǒng)節(jié)點數(shù)、支路參數(shù)、PMU安裝位置及量測值?？捎^測性分析檢查量測是否足夠讓H矩陣列滿秩。核心估計構(gòu)建z、H、R求解WLS狀態(tài)估計。壞數(shù)據(jù)檢測計算殘差做檢測與剔除。結(jié)果輸出展示估計值、誤差指標可視化對比。這個劃分也直接映射到MATLAB的代碼結(jié)構(gòu)上每個功能對應(yīng)一個函數(shù)main腳本負責串聯(lián)。后面第三章會給出各模塊的代碼骨架和關(guān)鍵參數(shù)選取邏輯。2. 核心模型與算法拆解2.1 PMU量測方程與狀態(tài)向量的數(shù)學形式先定義狀態(tài)向量。對于節(jié)點數(shù)為n的電力系統(tǒng)若所有節(jié)點都裝了PMU每個節(jié)點的狀態(tài)是電壓幅值V_i和相角θ_i那么理論上的狀態(tài)量總數(shù)是2n。但因為全網(wǎng)相角需要一個參考基準通常把某個參考節(jié)點的相角錨定為0實際上待估計的狀態(tài)數(shù)是2n-1。PMU的量測方程可以寫成z Hx e其中x是狀態(tài)向量e是量測噪聲向量。假設(shè)在母線i上裝了一臺PMU它能同時量測該母線的電壓幅值和相角還能量測與該母線相連的支路電流相量。以電壓量測為例對應(yīng)H矩陣的行是第i個幅值狀態(tài)位為1第i個相角狀態(tài)位為0對應(yīng)相角量測的行是第i個相角狀態(tài)位為1幅值狀態(tài)位為0。如果PMU還提供支路電流相量那需要基于線路的π型等值電路把電流相量表達式轉(zhuǎn)換為狀態(tài)量的線性組合這部分關(guān)鍵是要精確填H矩陣的系數(shù)稍有差錯估計結(jié)果就會偏掉。這里面有個細節(jié)很多人第一次做容易漏PMU給的是絕對相角以UTC為基準的全球同步相角而狀態(tài)估計里的狀態(tài)量是相對參考母線的相角。處理方式是把所有相角量測減去參考母線的相角量測再進入估計。否則H矩陣會多一列全零或者出現(xiàn)秩虧。2.2 加權(quán)矩陣R的選取為什么要按量測類型分權(quán)重WLS的目標函數(shù)是J (z-Hx)?R?1(z-Hx)R是量測誤差協(xié)方差矩陣。PMU的相量測量單元精度比傳統(tǒng)RTU高很多幅值誤差通常在0.1%量級相角誤差在0.01°~0.02°量級。但不同PMU通道、不同幅值和相角的方差差異還是存在所以R不能簡單設(shè)成單位陣。實踐中常用做法是查PMU的精度指標把幅值標準差和相角標準差換算成方差電壓幅值量測標準差約0.001~0.002 p.u.對應(yīng)R對角線元素約為1e-6~4e-6。電壓相角量測標準差約0.0002~0.0004 rad對應(yīng)R對角線元素約4e-8~1.6e-7。電流幅值和相角量測取決于CT/PT的精度等級通常比重會比電壓量測的方差大一些。權(quán)重的本質(zhì)是量測的信任度方差越小權(quán)重越大在求解時對結(jié)果的貢獻越大。這個道理和加權(quán)平均是一樣的實際項目里我習慣先用均勻權(quán)重跑一遍看殘差分布再用殘差方差反推R做一次迭代定權(quán)。這個小技巧能明顯改善結(jié)果。2.3 可觀測性分析為什么H矩陣必須滿秩線性狀態(tài)估計能解出唯一解的前提是量測方程個數(shù)不小于狀態(tài)數(shù)而且H矩陣列滿秩。列滿秩意味著每個狀態(tài)量都被足夠的獨立量測覆蓋不存在某個母線電壓怎么測都測不到的情況。在MATLAB里判斷很簡單計算秩rank(H)列滿秩的判據(jù)是rank(H)等于狀態(tài)數(shù)2n-1。同時還可以看條件數(shù)cond(H)條件數(shù)太大說明H矩陣近似病態(tài)哪怕滿秩數(shù)值上也可能解出漫天亂跳的結(jié)果。條件數(shù)控制在1e6以內(nèi)比較好超過這個量級就要警惕。這就引出一個部署問題PMU數(shù)量不夠怎么辦工程上常見的是PMU只裝在部分關(guān)鍵節(jié)點剩下的節(jié)點靠SCADA量測補齊。這種混合量測場景下模型重新變成非線性得用傳統(tǒng)WLS迭代。如果非要保持線性模型可以假設(shè)SCADA區(qū)域的狀態(tài)初始值已知那其實就不叫狀態(tài)估計了屬于擾動分析邏輯上要分清。2.4 壞數(shù)據(jù)檢測標準化殘差怎么用PMU數(shù)據(jù)也不是百分百干凈通信丟包、相量計算異常、GPS失步都會產(chǎn)生壞數(shù)據(jù)。線性模型下殘差r z - Hx_est理論上服從零均值高斯分布。采用的是基于標準化殘差的檢測r_i_normalized r_i / sqrt(R_ii * (I - H(H?R?1H)?1H?R?1)_ii)分子是第i個量測的殘差分母是殘差方差的開方。標準化之后r_i_normalized近似服從標準正態(tài)分布用閾值λ一般取3.0對應(yīng)99.7%置信度去卡超過閾值就判壞數(shù)據(jù)。有一個容易踩的坑多個壞數(shù)據(jù)同時存在時逐次剔除比一次性剔除更穩(wěn)。因為壞數(shù)據(jù)可能互相掩蓋殘差會被拉平單次殘差檢驗會漏掉。每次只剔除標準化殘差最大的那個量測重新做一遍估計再檢直到?jīng)]有超閾值的點為止。后面第四章會展開講這個問題的具體表現(xiàn)。3. MATLAB實操實現(xiàn)與核心環(huán)節(jié)3.1 數(shù)據(jù)準備我用IEEE 9節(jié)點系統(tǒng)作為測試床我復(fù)現(xiàn)這個項目時用的算例是IEEE 9節(jié)點系統(tǒng)節(jié)點數(shù)少、拓撲清楚、Matpower里有現(xiàn)成數(shù)據(jù)適合驗證算法。數(shù)據(jù)準備階段要明確幾樣?xùn)|西節(jié)點表9個節(jié)點包括基準電壓、類型PQ/PV/平衡。支路表每條支路的電阻、電抗、對地電納以及變壓器變比。PMU位置試驗時我假定節(jié)點1、3、6、9裝了PMU量測覆蓋這四點的電壓相量以及相連支路的電流相量。在MATLAB里我推薦用struct組織這些數(shù)據(jù)別用一堆散變量data.n 9; data.branch [ 1 4 0.0000 0.0576 0.0000 0; 4 5 0.0170 0.0920 0.1580 0; 5 6 0.0390 0.1700 0.3580 0; ... ]; data.pmu [1; 3; 6; 9]; data.z_meas YOUR_MEAS_VECTOR;如果手上沒有實測PMU數(shù)據(jù)可以用潮流計算結(jié)果作為真值再疊加高斯噪聲生成量測值。這個做法對驗證代碼正確性特別有用——因為你知道真值就能算誤差。3.2 核心求解函數(shù)從測量向量到狀態(tài)量這一段是代碼的樞紐。給出一個核心函數(shù)框架讀者可以直接改成自己的數(shù)據(jù)規(guī)模。function x_est pmu_wls_se(z, H, R) % z: 量測向量 m x 1 % H: 量測矩陣 m x (2n-1) % R: 量測誤差協(xié)方差矩陣 m x m G H * (R \ H); % 信息矩陣 b H * (R \ z); % 右端項 x_est G \ b; % 最小二乘解 end實際項目里H不是手工填的而是根據(jù)PMU位置和網(wǎng)絡(luò)拓撲動態(tài)生成。構(gòu)建H矩陣的邏輯分三步第一步建立狀態(tài)索引。每個節(jié)點分配兩個索引幅值索引idxV_i 2*(i-1)1相角索引idxTheta_i 2*(i-1)2。參考節(jié)點的相角索引要特殊處理要么在H中刪掉該列要么在x中固定為0并同步調(diào)整量測方程。第二步填電壓量測行。對于母線i的PMU幅值量測行在idxV_i位置填1相角量測行在idxTheta_i位置填1。相角量測如果用的是全局相角記得統(tǒng)一減去參考節(jié)點的全局相角后再進估計。第三步填支路電流量測行。以π型等值電路為準先算線路導(dǎo)納Y_ij G jB再根據(jù)電流相量I_ij Y_ij(V_i - V_j) jBsh/2 * V_i把實部和虛部對狀態(tài)量的偏導(dǎo)算出來填到對應(yīng)位置。這一步最容易寫錯強烈建議先用一個簡單兩節(jié)點系統(tǒng)驗證H矩陣的正確性。3.3 結(jié)果驗證誤差分析和殘差分析狀態(tài)估計算完不能直接交差得驗證。我的做法是把估計值跟潮流真值做對比計算每個節(jié)點的幅值誤差和相角誤差。通常用RMSE均方根誤差來評價整體精度公式是RMSE sqrt(mean((x_est - x_true).^2))從我的測試結(jié)果看在R矩陣設(shè)置合理、量測噪聲標準差符合PMU實際水平的前提下9節(jié)點系統(tǒng)的幅值估計誤差大約在1e-4 p.u.量級相角估計誤差大約在1e-3 rad量級。如果誤差偏大一個數(shù)量級以上優(yōu)先檢查H矩陣有沒有填錯其次是R矩陣是否給了不合理的權(quán)重。另外要畫殘差分布圖。把標準化殘差畫成條形圖直觀能看出有沒有異常量測。正常情況下殘差密布在±3之間且沒有明顯單點突出如果有突出點先別急著刪確認是數(shù)據(jù)問題還是H矩陣問題。3.4 整體腳本結(jié)構(gòu)與運行流程main腳本的結(jié)構(gòu)可以這樣安排% 第1步加載系統(tǒng)數(shù)據(jù)和PMU量測 system_data load_pmu_system(ieee9); % 第2步構(gòu)建H矩陣和R矩陣 [H, R, z, idx] build_linear_model(system_data); % 第3步可觀測性檢查 assert(rank(H) size(H,2), H矩陣秩虧系統(tǒng)不可觀測); % 第4步狀態(tài)估計 x_est pmu_wls_se(z, H, R); % 第5步壞數(shù)據(jù)檢測 [r_norm, bad_idx] bad_data_detect(z, H, R, x_est); % 第6步結(jié)果輸出 plot_result(x_est, system_data);這里每步調(diào)用一個函數(shù)函數(shù)內(nèi)部再細分可讀性和可維護性都比寫一個兩三百行的主腳本強得多。有一個小建議運行前用tic/toc記錄時間線性模型下9節(jié)點系統(tǒng)的計算時間應(yīng)該在毫秒級G矩陣是17x17MATLAB分解起來非?？臁Ｈ绻艹鰜硪脦酌牖究梢源_定有冗余循環(huán)或者H矩陣構(gòu)建邏輯低效需要排查。4. 常見問題與排查技巧實錄4.1 H矩陣奇異或條件數(shù)過大怎么辦這是我在項目里遇到最多的一個問題幾乎每個初做PMU狀態(tài)估計的人都會撞上一次。表現(xiàn)就是rank(H) size(H,2)或者cond(H)在1e12以上估計值亂跳?？赡茉蛴腥齻€。一是PMU覆蓋不足某些母線完全沒有量測覆蓋對應(yīng)H行全是零狀態(tài)不可觀。二是參考相角沒有處理干凈H矩陣里含有一列全零或者兩列線性相關(guān)。三是線路參數(shù)填錯導(dǎo)致支路電流量測對應(yīng)行與電壓量測行產(chǎn)生線性相關(guān)關(guān)系。排查順序第一步打印H矩陣的稀疏模式用spy(H)看哪列全零哪兩列是成比例關(guān)系。第二步逐個PMU檢查量測方程個數(shù)和類型確認覆蓋范圍。第三步用一個只有兩臺PMU的兩節(jié)點系統(tǒng)做單元測試H矩陣規(guī)模小一眼能看出問題。提示如果是可觀測性不足不要試圖在代碼層面打補丁。老老實實增加PMU量測或者把部分SCADA量測補進模型用混合量測的非線性WLS。強行求解只會得到數(shù)值上看起來正常、實際上完全沒意義的解。4.2 相角參考基準沖突PMU量測給出的是全球同步的絕對相角不同廠家的PMU在接入同一系統(tǒng)時由于GPS信號處理延遲的差異可能會有微小的角度偏移。如果直接把不同PMU的相角拿來拼成同一個z向量容易出現(xiàn)系統(tǒng)性偏差。我的處理辦法是在數(shù)據(jù)預(yù)處理階段先把所有PMU相角量測減去同一個參考PMU的相角量測得到相對相角序列再進入估計。這樣即使PMU本身有固定延遲誤差只要延遲在短時間內(nèi)穩(wěn)定相減后誤差會被抵消掉。這個操作在代碼里就是一行z_theta z_theta - z_theta_ref。4.3 壞數(shù)據(jù)檢測的誤檢與漏檢誤檢通常是因為R矩陣的方差設(shè)得太小量測噪聲本來沒那么高精度標準化之后殘差就會偏大超過閾值。漏檢則常見于兩個壞數(shù)據(jù)互相抵消的場景。舉個例子某條支路兩端的電流量測同時壞掉它們的殘差可能方向相反平均下來標準化殘差不大就會被漏掉。實操中我的建議是別只依賴一次殘差檢驗。做一個循環(huán)剔除法——每次只刪標準化殘差最大的那個量測重新估計后再檢。同時把檢驗閾值從3.0放寬到2.8多捕獲一些邊緣可疑點寧可多剔除一個可疑量測也不能放壞數(shù)據(jù)進門。當然剔除的量測數(shù)不能太多一般超過總量測數(shù)的5%就要回頭檢查是不是數(shù)據(jù)質(zhì)量整體不行。4.4 量測噪聲設(shè)置與實際不符測試階段想要模擬真實場景可以在潮流真值上加高斯白噪聲但噪聲標準差的選擇要有依據(jù)。我見過有人直接用randn加噪聲標準差設(shè)成0.1結(jié)果估計誤差大得離譜還以為是算法有問題實際上是噪聲水平跟PMU的真實指標差了百倍。PMU的幅值測量精度典型值在0.1%相角測量精度在0.01°左右對應(yīng)弧度約1.7e-4 rad。按這個量級設(shè)置噪聲估計結(jié)果才能反映算法本身的性能。做完之后可以統(tǒng)計殘差的標準差和設(shè)置的噪聲標準差做對比如果差太多說明H矩陣或者加權(quán)有問題。4.5 計算效率與內(nèi)存小技巧PMU量測點數(shù)一旦多了H矩陣和R矩陣的維度會漲得很快。比如IEEE 118節(jié)點系統(tǒng)全裝PMU狀態(tài)量就是235個量測可能有上千行。這時候如果H還是稠密矩陣求逆操作會越來越慢。兩個改進方向一是用稀疏矩陣存儲H和RMATLAB里直接sparse(H)二是用信息矩陣G HRH然后對G做Cholesky分解而不是直接求H的偽逆。這兩步能把計算時間下降兩個數(shù)量級。我在118節(jié)點系統(tǒng)上試過從幾十秒降到幾百毫秒效果非常明顯。另外一個細節(jié)R矩陣是對角陣R \ H這一步不要寫成inv(R) * H直接用左除效率更高數(shù)值也穩(wěn)定得多。MATLAB里對稀疏對角陣的除法有專門優(yōu)化一定要利用上。最后再分享一個擴展思路。這個靜態(tài)估計框架跑通之后如果還想往深做可以嘗試兩個方向一是把PMU量測數(shù)據(jù)按時間序列連續(xù)輸入加入狀態(tài)轉(zhuǎn)移模型升級成動態(tài)狀態(tài)估計這時就該上卡爾曼濾波了二是在現(xiàn)有框架中把H矩陣的構(gòu)建推廣到三相不平衡系統(tǒng)就能用來處理配電網(wǎng)狀態(tài)估計。我個人做下來最大的體會是這套模型的代碼骨架一旦搭好往各種方向擴展都很快關(guān)鍵是前期的數(shù)據(jù)結(jié)構(gòu)和H矩陣構(gòu)建邏輯要寫干凈別為省幾行代碼把后續(xù)的擴展性毀了。