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