一迭代法建模與實現(xiàn))
1. 交直流潮流為什么最后落在了統(tǒng)一迭代上第一次把交直流混合系統(tǒng)往 Matlab 里塞的時候我的做法和大多數(shù)人一樣交流側(cè)用一套牛頓法算完把換流站交流側(cè)的有功當已知量丟給直流網(wǎng)絡(luò)直流網(wǎng)絡(luò)解出來再回代更新?lián)Q流站功率來回往復(fù)直到兩次迭代的換流站功率差小于閾值。這套交替迭代法也有人叫它順序求解法跑純主從控制的簡單算例時確實舒服代碼短、調(diào)試直觀、交流子程序和直流子程序各自獨立出問題一眼就知道是哪半邊。但只要換流站數(shù)量一多、控制方式一混搭、再加上換流器損耗隨電壓波動這種細枝末節(jié)交替迭代就開始露怯迭代次數(shù)從七八次漲到三四十次偶爾還會在兩個解之間來回振蕩你把收斂閾值放寬一點它就收斂了檢查結(jié)果發(fā)現(xiàn)直流電壓偏離參考值上千分之幾功率平衡差得離譜?；?Matlab 的交直流潮流計算統(tǒng)一迭代法建模與實現(xiàn)本質(zhì)上就是在講一件事把交流節(jié)點方程和直流節(jié)點方程塞進同一個牛頓迭代框架里用一個大雅可比矩陣同時求解用額外的矩陣維度換取收斂的魯棒性和解的唯一性。它不算什么高深算法真正難的地方全在建模細節(jié)上——換流站的穩(wěn)態(tài)方程怎么擺、標幺基準怎么統(tǒng)一、耦合偏導(dǎo)要不要留、直流松弛節(jié)點怎么處理。這些細節(jié)一旦搞錯程序照樣能跑出數(shù)來但結(jié)果是錯的而且錯得很隱蔽。這篇東西適合三類人看第一類是在做交直流混合系統(tǒng)相關(guān)課題、需要自己寫潮流程序而不是只會點軟件按鈕的第二類是想搞清楚為什么我的交替迭代不收斂的第三類是有一定 Matlab 基礎(chǔ)、能看懂矩陣運算但沒接觸過電力系統(tǒng)潮流內(nèi)部結(jié)構(gòu)的。我會把方程怎么來、雅可比怎么裝、代碼怎么寫、坑在哪里踩一條線講下來。默認你會基礎(chǔ)的交流牛頓-拉夫遜潮流如果連 Ybus 怎么組、H 矩陣怎么來的都不清楚建議先把純交流版本手推一遍再往下看不然讀到雅可比那節(jié)會很痛苦。需要提前說明一點下面給出的算例參數(shù)和迭代記錄是我為了講清楚流程自己搭的一個小系統(tǒng)不是任何實際工程數(shù)據(jù)你在復(fù)現(xiàn)時數(shù)值對不上很正常重點看結(jié)構(gòu)和邏輯。2. 方案選型交替迭代與統(tǒng)一迭代到底差在哪2.1 交替迭代法在工程里卡在什么地方交替迭代的思路非常工程師把大問題切成兩個能各自求解的小問題輪流解用邊界量互相傳遞。它的收斂性依賴一個隱含假設(shè)——交流側(cè)和直流側(cè)的耦合不能太強。主從控制下定直流電壓的那個換流站幾乎把兩側(cè)釘死了耦合確實弱交替兩三輪就穩(wěn)。但一旦上 droop 控制直流電壓本身參與功率分配換流站有功不再由交流側(cè)給定而是由直流電壓決定此時交流方程里出現(xiàn)未知量 Vdc直流方程里出現(xiàn)未知量 Vm兩個子系統(tǒng)互相牽著走交替迭代的收斂速度會肉眼可見地掉下來。我在一個五端直流的算例上做過對比全主從控制時交替迭代 6 次收斂閾值 1e-8把其中兩端改成下垂控制后交替迭代需要 31 次而且下垂系數(shù)取大一點就發(fā)散。這不是算法實現(xiàn)有 bug而是交替迭代本質(zhì)上是一個塊高斯-賽德爾迭代耦合項的譜半徑隨耦合強度增大而逼近 1收斂自然就慢。統(tǒng)一迭代相當于把它換成了整體牛頓二次收斂的性質(zhì)量級上不一樣。還有一個更現(xiàn)實的問題交替迭代的收斂判據(jù)不好定。你按換流站功率差判斷可能出現(xiàn)功率差已經(jīng)很小但電壓還在漂的情況你按電壓差判斷又可能在功率嚴重不平衡時誤判收斂。統(tǒng)一迭代沒有這個煩惱不平衡量向量統(tǒng)一范數(shù)一個判據(jù)管到底。2.2 統(tǒng)一迭代付出的代價值不值代價很直接雅可比維度變大裝配邏輯變復(fù)雜調(diào)試難度上升。一個 nb 個交流節(jié)點、nd 個直流節(jié)點的系統(tǒng)純交流雅可比是 (2nb-1) 階左右加上直流后大概多出 nd 階而且多出來的耦合塊元素推導(dǎo)容易出錯。更要命的是一旦某個耦合偏導(dǎo)符號寫反程序不一定發(fā)散而是收斂到一個物理上不成立的解你拿萬用表都量不出來。但從我自己的經(jīng)驗看這筆賬是劃算的。第一牛頓法對初值敏感但一旦進入收斂域就是二次收斂交直流系統(tǒng)初值可以很方便地用平啟動 直流電壓全置 1.0給出實際很少不收斂第二統(tǒng)一迭代天然支持各種控制方式混搭主從、下垂、定功率、定電壓只需要改對應(yīng)行的方程框架不動第三擴展性好以后要加換流器限流、交流側(cè)無功限值、直流線路損耗都只是在已有方程上加項不像交替迭代那樣要重新設(shè)計外層循環(huán)。代價還有一條容易被忽略統(tǒng)一迭代對模型的完整性要求高。交替迭代可以容忍某一側(cè)用簡化模型比如直流側(cè)只當恒功率負載統(tǒng)一迭代不行你寫進雅可比的每一項都必須和方程嚴格對應(yīng)否則收斂性立刻崩。這其實是好事逼著我把換流器方程老老實實推了一遍。2.3 什么情況下我會勸你別用統(tǒng)一迭代不是所有場景都值得上統(tǒng)一迭代。如果系統(tǒng)里只有一個背靠背換流站兩側(cè)交流網(wǎng)絡(luò)規(guī)模都很小交替迭代兩三次就收斂代碼量還只有統(tǒng)一迭代的三分之一那就沒必要折騰。如果直流側(cè)是完全理想的恒功率負載模型損耗忽略、控制理想交直流之間嚴格解耦統(tǒng)一迭代的雅可比是分塊對角的等于把兩個獨立問題拼在一起解白白浪費計算量。真正值得上統(tǒng)一迭代的典型場景是多端直流、下垂控制、換流站數(shù)量大于等于三、需要嚴格校驗功率平衡、后續(xù)要做 N-1 掃描或連續(xù)潮流。這種情況下交替迭代的迭代次數(shù)不穩(wěn)定批量計算時總有幾個算例卡住統(tǒng)一迭代的穩(wěn)定收斂反而是省時間的。對比維度交替迭代法統(tǒng)一迭代法雅可比規(guī)模兩個小矩陣分別求解一個大矩陣一次求解收斂階數(shù)線性塊高斯-賽德爾二次牛頓控制方式混搭需要改外層調(diào)度邏輯只改對應(yīng)行方程下垂控制支持收斂慢易振蕩與主從控制同框架調(diào)試難度低兩側(cè)獨立高耦合項易錯代碼量我的實現(xiàn)約 300 行約 600 行擴展新設(shè)備改動大加行加列即可3. 模型拆解把換流站寫成能進雅可比的方程3.1 節(jié)點分類與未知量編號統(tǒng)一迭代第一個要想清楚的問題不是方程而是誰是未知量。交流側(cè)沿用經(jīng)典分類PQ 節(jié)點未知量是電壓幅值 Vm 和相角 VaPV 節(jié)點未知量只有相角幅值由控制給定平衡節(jié)點兩個都給定。直流側(cè)稍微特殊一點直流網(wǎng)絡(luò)沒有相角概念只有電壓幅值 Vdc所以每個直流節(jié)點的未知量就是 Vdc 本身但定直流電壓的那個節(jié)點例外——它的 Vdc 已知功率未知等價于交流側(cè)的平衡節(jié)點。這個例外是很多人第一次寫統(tǒng)一迭代時栽跟頭的地方。直覺上你會覺得直流節(jié)點都應(yīng)該有功率方程于是給定了電壓的節(jié)點也列了一條 ΔPdc 方程結(jié)果雅可比對應(yīng)行是零行或者接近零行矩陣直接奇異。正確做法是定電壓節(jié)點的 Vdc 從未知量集合里剔除對應(yīng)的 ΔPdc 方程也不列它的功率在收斂后由 Pdc Vdc * Σ G Vdc 反算出來。我用的是未知量索引映射的辦法先掃一遍表把所有未知量編號再用映射數(shù)組把物理節(jié)點號翻譯成未知量序號。這樣后面裝配雅可比時不用到處判斷類型索引一查就行。具體做法在第四節(jié)給代碼。3.2 換流站的穩(wěn)態(tài)方程怎么擺換流站建模有兩條路平均模型和開關(guān)模型。潮流計算用平均模型忽略開關(guān)過程和 PWM 諧波只保留基波功率平衡。我采用的是最常見的功率平衡寫法交流側(cè)注入電網(wǎng)的有功 Pc直流側(cè)從直流網(wǎng)絡(luò)取用的有功 Pdc。兩者之差是換流器損耗Pdc Pc Ploss。交流側(cè)無功 Qc 由控制方式?jīng)Q定定 Q 控制時直接給定。損耗模型取 Ploss a b·I c·I2其中 I 是換流站交流側(cè)電流的標幺值I sqrt(Pc2 Qc2) / Vm。選這個損耗模型不是因為它精確而是因為它同時包含了常數(shù)項、線性項和二次項能覆蓋導(dǎo)通損耗和開關(guān)損耗的主要成分而且對 Vm 的偏導(dǎo)好推——這一點在裝配耦合塊時非常關(guān)鍵。有些文獻用恒定效率模型Pdc Pc / η簡單但那樣 Ploss 對 Vm 的偏導(dǎo)恒為零交直流之間就沒有耦合項了統(tǒng)一迭代退化成兩個獨立問題失去了演示意義。注意損耗模型的選擇直接決定耦合塊是否存在。如果你的算例里所有換流站都用恒定效率、所有控制都是定 P 定 Q那統(tǒng)一迭代的雅可比是分塊對角的收斂性和交替迭代一樣這時候沒必要上統(tǒng)一迭代用它只是自己給自己找麻煩。3.3 標幺基準怎么統(tǒng)一交直流混合系統(tǒng)最容易出錯的地方是基準值。交流側(cè)一般取三相功率基準 Sb 和線電壓基準 Vb直流側(cè)如果直接沿用交流的 Sb直流電壓基準要單獨取——我習慣取 Vdb 使得直流側(cè)基準阻抗與交流側(cè)一致也就是 Vdb2 / Sb Vb2 / Sb這樣 Vdb Vb這里的 Vb 是交流相電壓基準若用線電壓基準則要乘系數(shù)。如果不統(tǒng)一直流網(wǎng)絡(luò)的電導(dǎo)矩陣 G 數(shù)值會偏離好幾個數(shù)量級雅可比條件數(shù)爆炸。我吃過一次虧直流側(cè)電導(dǎo)一開始用的是實際歐姆值沒轉(zhuǎn)標幺結(jié)果雅可比對角元量級從 1e1 跳到 1e4稀疏求解器雖然沒報錯但迭代過程出現(xiàn)了明顯的數(shù)值抖動收斂判據(jù)在 1e-6 附近停住不動。改成標幺后一次就收斂了。所以我的建議是建模階段就把所有量統(tǒng)一到同一套基準下代碼里只出現(xiàn)標幺值物理量只在輸入輸出環(huán)節(jié)轉(zhuǎn)換。量基準轉(zhuǎn)換關(guān)系常見錯誤交流功率SbP_pu P / Sb單相/三相混淆交流電壓Vb線電壓Vm_pu V / Vb用相電壓基準直流電壓Vdb VbVdc_pu Vdc / Vdb直接用電網(wǎng)額定值直流電阻Zdb Vdb2/SbR_pu R / Zdb忘了轉(zhuǎn)標幺換流器電流Ib Sb/(√3·Vb)I_pu I / Ib與功率混淆損耗系數(shù)按 Sb 歸一a_pu a / Sb系數(shù)單位不統(tǒng)一3.4 完整不平衡量向量的構(gòu)成把所有方程按順序排好不平衡量向量 ΔF 分成三塊第一塊是交流有功不平衡 ΔP對所有非平衡交流節(jié)點列共 (nb - 1) 條。對純負荷節(jié)點ΔP_i P_i^calc - P_i^spec對掛換流站的節(jié)點注入功率里要減去換流站交流側(cè)取用的 Pc。第二塊是交流無功不平衡 ΔQ對所有 PQ 節(jié)點列。掛換流站的節(jié)點如果換流站定 Q 控制那 Qc 作為給定值進 P_spec 的位置。第三塊是直流功率不平衡 ΔPdc對所有非定電壓直流節(jié)點列ΔPdc_k Vdc_k·Σ_m G_km·Vdc_m - Pdc_k^spec。這里有個符號約定的坑必須講清楚。換流站整流交流到直流時 Pc 為正、Pdc 為負還是反過來不同文獻寫法不一樣。我統(tǒng)一約定Pc 為換流站從交流電網(wǎng)吸收的有功正表示吸收Pdc 為換流站注入直流網(wǎng)絡(luò)的有功正表示注入那么整流時 Pc 0、Pdc 0 的絕對值小于 Pc關(guān)系是 Pdc -(Pc - Ploss)。逆變時反過來。約定一旦定下來全套代碼都要遵守中間改一次符號所有耦合項都得重推。4. 雅可比矩陣的結(jié)構(gòu)與裝配細節(jié)4.1 分塊結(jié)構(gòu)與稀疏索引統(tǒng)一迭代的雅可比是一個三行三列的分塊矩陣J [H N Ndc; M L Ldc; Md Hd D ]H、N、M、L 是大家熟悉的交流子塊分別是 ΔP 對 Va、ΔP 對 Vm、ΔQ 對 Va、ΔQ 對 Vm 的偏導(dǎo)。D 是直流子塊ΔPdc 對 Vdc 的偏導(dǎo)。剩下四塊是耦合塊Ndc 是 ΔP 對 VdcLdc 是 ΔQ 對 VdcMd 是 ΔPdc 對 VaHd 是 ΔPdc 對 Vm。在我采用的換流站模型下Ndc 和 Ldc 恒為零——因為交流側(cè)注入功率 Pc、Qc 是控制給定的不直接依賴直流電壓。Md 也是零——直流功率不依賴交流相角。只有 Hd 非零它來自損耗項對交流電壓幅值的依賴。這個結(jié)論很重要它說明統(tǒng)一迭代的耦合其實只有一個塊在起作用如果你把這個塊也簡化掉損耗恒定矩陣就是分塊下三角可以用兩次回代求解本質(zhì)上又退回了交替迭代。所以我的原則是耦合塊哪怕數(shù)值很小也必須保留它是統(tǒng)一迭代區(qū)別于交替迭代的物理來源。實際算例里 Hd 的元素量級大概在 1e-3 到 1e-2確實小但正是這個小項在迭代初期提供了方向信息減少了來回振蕩。4.2 交流子塊的經(jīng)典寫法回顧交流部分沿用標準極坐標牛頓法公式我不重復(fù)推導(dǎo)只強調(diào)兩個實現(xiàn)細節(jié)。第一Ybus 要用稀疏矩陣Matlab 里 sparse 存儲裝配功率時用 Ybus 乘以電壓向量別用循環(huán)。第二對角元和非對角元的偏導(dǎo)公式不同裝配時按行列分別處理不要圖省事用一個循環(huán)套所有元素。H 的對角元H_ii -Q_i - B_ii·Vm_i2 H 的非對角元H_ij Vm_i·Vm_j·(G_ij·sin θ_ij - B_ij·cos θ_ij)N 的對角元N_ii P_i/Vm_i G_ii·Vm_i N 的非對角元N_ij Vm_i·(G_ij·cos θ_ij B_ij·sin θ_ij)M、L 類似符號別搞反。這部分是教科書內(nèi)容但我在實際 debug 中發(fā)現(xiàn)很多人栽在 θ_ij 的定義上——θ_ij θ_i - θ_j 還是 θ_j - θ_i兩種寫法都有人用混用會導(dǎo)致 H 的非對角元整體符號反了程序仍然收斂但結(jié)果完全錯誤。我的做法是全程用 θ_ij θ_i - θ_j并且寫一個單元測試兩節(jié)點系統(tǒng)手算一遍對拍。4.3 直流子塊與耦合子塊的偏導(dǎo)推導(dǎo)直流注入功率 Pdc_k Vdc_k · Σ_m G_km · Vdc_m。對 Vdc_m 求偏導(dǎo)m ≠ k 時?Pdc_k/?Vdc_m Vdc_k · G_kmm k 時?Pdc_k/?Vdc_k Σ_m G_km·Vdc_m Vdc_k·G_kk這個對角元很容易寫錯成 2·Vdc_k·G_kk或者漏掉求和項里的 k 項。正確形式是該節(jié)點注入電流 自身電導(dǎo)乘電壓物理上就是電流對電壓的靈敏度。耦合塊 Hd 的推導(dǎo)ΔPdc_k 里包含 Pdc_k^spec而 Pdc_k^spec Pc_k Ploss_k按前面的符號約定取絕對值Ploss_k a b·I c·I2I sqrt(Pc2 Qc2)/Vm_iVm_i 是換流站所接交流節(jié)點的電壓幅值。?Ploss/?Vm (b 2c·I) · ?I/?Vm (b 2c·I) · (-I/Vm)所以 ?ΔPdc_k/?Vm_i -?Ploss/?Vm_i (b 2c·I)·I/Vm_i。符號為正物理意義是電壓幅值升高、電流減小、損耗減小、直流側(cè)所需注入功率增加合理。提示這個偏導(dǎo)的量級很小如果你在調(diào)試時發(fā)現(xiàn)它一直是零先檢查是不是損耗系數(shù) b、c 都設(shè)成了零或者 Vm 初值取 1.0 導(dǎo)致靈敏度看起來沒變化?？梢杂脭?shù)值差分驗證一遍方法在下一節(jié)。4.4 用數(shù)值差分給雅可比做體檢手推偏導(dǎo)最怕符號錯。我的習慣是寫完解析雅可比后立刻用中心差分對拍一次誤差超過 1e-4 就說明有地方錯了。代碼不到二十行但能省掉幾個晚上的調(diào)試function Jnum jac_check(x, Ffun, nx, h) % x : 當前狀態(tài)向量 % Ffun: 不平衡量函數(shù)句柄 F Ffun(x) % nx : 未知量個數(shù) % h : 差分步長取 1e-6 左右 F0 Ffun(x); Jnum zeros(numel(F0), nx); for k 1:nx xp x; xp(k) xp(k) h; xm x; xm(k) xm(k) - h; Jnum(:, k) (Ffun(xp) - Ffun(xm)) / (2*h); end end用法是先在第一次迭代時調(diào)用把解析雅可比 J 和 Jnum 相減取絕對值最大元素打印出來。如果最大值出現(xiàn)在耦合塊對應(yīng)的行說明耦合偏導(dǎo)有問題如果出現(xiàn)在交流子塊說明 Ybus 或者 θ 約定有誤。差分步長不宜太小1e-6 到 1e-8 之間試太大會有截斷誤差太小會撞上雙精度舍入誤差。5. Matlab 實現(xiàn)從數(shù)據(jù)結(jié)構(gòu)到求解循環(huán)5.1 三張表的數(shù)據(jù)結(jié)構(gòu)我不用面向?qū)ο笠膊挥?struct 數(shù)組就用三個二維矩陣列定義寫在注釋里。好處是索引快、易于向量化、保存成 mat 文件方便。缺點是列含義靠記憶所以我把列定義固化成常量。% bus : 交流節(jié)點表每行一個節(jié)點 % 1 id | 2 type(1PQ,2PV,3平衡) | 3 Vm | 4 Va(rad) % 5 Pg | 6 Qg | 7 Pd | 8 Qd | 9 conv_id(0無換流站) % % dc : 直流節(jié)點表 % 1 id | 2 Vdc | 3 type(1定電壓,2定功率,3下垂) | 4 Pdc_spec % 5 kdroop | 6 Vdc_ref | 7 conv_id % % conv : 換流站表 % 1 id | 2 ac_bus | 3 dc_node | 4 mode(1定P定Q,2定Vdc) % 5 Pset | 6 Qset | 7 a | 8 b | 9 ctype字段我堅持用整數(shù)而不是字符串因為后面要大量做邏輯判斷和索引字符串比較慢且容易拼錯。conv_id這個外鍵放在 bus 和 dc 表里是雙向索引的關(guān)鍵裝配耦合塊時靠它從交流節(jié)點找到對應(yīng)的直流節(jié)點反之亦然。5.2 未知量映射與主循環(huán)骨架主循環(huán)的結(jié)構(gòu)固定算不平衡量 → 判斷收斂 → 裝雅可比 → 解修正方程 → 更新狀態(tài)。真正需要小心的是狀態(tài)更新的邊界處理以及修正量對應(yīng)的物理量。function res acdc_nr(bus, dc, conv, opt) nb size(bus,1); nd size(dc,1); % 交流未知量索引 iVa find(bus(:,2) ~ 3); % 非平衡節(jié)點的相角 iVm find(bus(:,2) 1); % PQ 節(jié)點的幅值 % 直流未知量索引剔除定電壓節(jié)點 iVd find(dc(:,3) ~ 1); nx numel(iVa) numel(iVm) numel(iVd); % 狀態(tài)向量 x 與物理表的映射 mapVa zeros(nb,1); mapVa(iVa) 1:numel(iVa); mapVm zeros(nb,1); mapVm(iVm) numel(iVa) (1:numel(iVm)); mapVd zeros(nd,1); mapVd(iVd) numel(iVa)numel(iVm) (1:numel(iVd)); Vm bus(:,3); Va bus(:,4); Vdc dc(:,2); Y make_ybus(bus); G make_gdc(dc); for it 1:opt.maxit [dP, dQ] ac_mismatch(Y, Vm, Va, bus, conv, Vdc); dPdc dc_mismatch(G, Vdc, dc, conv, Vm); F [dP(iVa); dQ(iVm); dPdc(iVd)]; if norm(F, inf) opt.tol break; end J build_jac(Y, Vm, Va, G, Vdc, bus, dc, conv, ... iVa, iVm, iVd, mapVa, mapVm, mapVd); dx J \ (-F); % 稀疏左除 % 回代更新 Va(iVa) Va(iVa) dx(mapVa(iVa)); Vm(iVm) Vm(iVm) dx(mapVm(iVm)); Vdc(iVd) Vdc(iVd) dx(mapVd(iVd)); end res pack_result(bus, dc, conv, Vm, Va, Vdc, it, norm(F,inf)); end這段骨架里有個細節(jié)值得說J \ (-F)用的是左除Matlab 會自動識別 J 是稀疏矩陣并調(diào)用合適的求解器。如果你寫inv(J)*(-F)一是慢二是精度差三是在大系統(tǒng)上會直接內(nèi)存溢出。另外norm(F, inf)取無窮范數(shù)比二范數(shù)更嚴格因為它盯的是單個最大的不平衡量不會被平均掉。5.3 交流不平衡量的計算function [dP, dQ] ac_mismatch(Y, Vm, Va, bus, conv, Vdc) V Vm .* exp(1j*Va); S V .* conj(Y * V); P real(S); Q imag(S); Ps (bus(:,5) - bus(:,7)) / 1.0; % 給定注入已標幺 Qs (bus(:,6) - bus(:,8)) / 1.0; % 換流站交流側(cè)取用功率計入負荷 for k 1:size(conv,1) i conv(k,2); c k; [Pc, Qc, ~] conv_ac_power(conv(c,:), Vm(i), Vdc(conv(c,3))); Ps(i) Ps(i) - Pc; Qs(i) Qs(i) - Qc; end dP P - Ps; dQ Q - Qs; end這里conv_ac_power是核心它要根據(jù)控制方式返回換流站交流側(cè)實際取用的有功和無功。定 P 定 Q 模式下 Pc Pset、Qc Qset直接返回定直流電壓模式下Pc 是由直流側(cè)功率反推的未知量不能在這里給死需要靠方程間的關(guān)系解決——這也是統(tǒng)一迭代相比交替迭代省事的地方未知量 Pc 不出現(xiàn)在狀態(tài)向量里但它通過直流方程隱含求解。嚴格來說定電壓站的 Pc 在每次迭代時應(yīng)該由當前 Vdc 和損耗反算Pc Pdc - Ploss而 Pdc Vdc·ΣG·Vdc。這就是所謂的內(nèi)層更新不影響雅可比結(jié)構(gòu)但影響收斂速度。我一開始沒做這個內(nèi)層更新直接用上一次迭代的 Pc 值結(jié)果迭代次數(shù)多了快一倍。5.4 直流不平衡量與稀疏裝配function dPdc dc_mismatch(G, Vdc, dc, conv, Vm) Idc G * Vdc; Pdc Vdc .* Idc; Pspec zeros(size(Vdc)); for k 1:size(conv,1) j conv(k,3); [~, ~, Pdc_k] conv_ac_power(conv(k,:), Vm(conv(k,2)), Vdc(j)); Pspec(j) Pdc_k; end dPdc Pdc - Pspec; end雅可比裝配我采用三元組預(yù)分配的方式先用spalloc估一個非零元上界用三個長向量 irow、jcol、jval 累積最后一次性sparse。這種方式比逐元素賦值快一個數(shù)量級在節(jié)點數(shù)上百之后差別非常明顯。nnz_est 4*(nb nd) 16*size(conv,1); irow zeros(nnz_est,1); jcol irow; jval irow; cnt 0; add (r,c,v) deal_helper(); % 這里用局部函數(shù)累加略實際寫的時候我把累加封裝成局部函數(shù)避免匿名函數(shù)帶來的性能損失function push(r, c, v) cnt cnt 1; irow(cnt) r; jcol(cnt) c; jval(cnt) v; end這種寫法在嵌套函數(shù)里能直接改外部變量比傳結(jié)構(gòu)體清爽。最后J sparse(irow(1:cnt), jcol(1:cnt), jval(1:cnt), nx, nx)。5.5 收斂判據(jù)與迭代控制判據(jù)我一般用兩層主判據(jù)是無窮范數(shù)小于 1e-8標幺輔助判據(jù)是修正量范數(shù)小于 1e-6兩個都滿足才認為收斂。只判不平衡量的話在病態(tài)系統(tǒng)上可能出現(xiàn)不平衡量很小但電壓還在緩慢漂移的情況。迭代次數(shù)上限設(shè) 30超過就認為不收斂打印最后一次的不平衡量分布看是哪個節(jié)點貢獻最大。這個診斷信息非常有用——如果最大不平衡量集中在某個換流站節(jié)點八成是控制方式設(shè)置和方程不匹配如果分散在所有節(jié)點八成是初值或者標幺化的問題。還有一個實用技巧記錄每次迭代的不平衡量范數(shù)畫成半對數(shù)圖。正常牛頓法應(yīng)該是一條陡峭的直線下降如果出現(xiàn)平臺期、抖動或者先降后升說明雅可比有問題或者系統(tǒng)本身無解。這個圖我?guī)缀趺看握{(diào)試新模型都會畫比看數(shù)字直觀得多。6. 一個可復(fù)現(xiàn)的小算例與結(jié)果校驗6.1 算例配置我搭了一個五節(jié)點交流 三節(jié)點直流的系統(tǒng)規(guī)模小到能手算校驗又能體現(xiàn)出交直流耦合的特征。交流側(cè)節(jié)點 1 是平衡節(jié)點節(jié)點 2、3 是純負荷節(jié)點 4、5 各接一個換流站。直流側(cè)節(jié)點 1 是定電壓站對應(yīng)交流節(jié)點 4 的換流站節(jié)點 2、3 是定功率站對應(yīng)交流節(jié)點 5 的換流站和一個純直流負荷。參數(shù)數(shù)值說明交流基準功率100 MVA全系統(tǒng)統(tǒng)一交流基準電壓10 kV線電壓交流線路阻抗0.01 j0.08 pu環(huán)形結(jié)構(gòu)負荷節(jié)點 20.6 j0.2 pu恒功率負荷節(jié)點 30.4 j0.15 pu恒功率換流站 1 控制定 Vdc 1.0主從控制的從站換流站 2 控制定 P -0.3 pu逆變運行直流線路電阻0.02 pu三節(jié)點三角形損耗系數(shù) a0.002 pu常數(shù)項損耗系數(shù) b0.01線性項損耗系數(shù) c0.02二次項需要注意的是這里的定功率是站在直流側(cè)說的換流站 2 從直流網(wǎng)絡(luò)取 0.3 pu 的功率送到交流側(cè)同時因為損耗它從直流網(wǎng)絡(luò)實際取用的功率要略大于 0.3。6.2 迭代過程記錄平啟動所有電壓置 1.0相角置 0。前三次迭代的不平衡量范數(shù)大致是迭代次數(shù)不平衡量范數(shù)說明06.2e-1初始不平衡主要是負荷節(jié)點18.7e-2大幅下降牛頓方向正確24.1e-3進入二次收斂區(qū)39.6e-6繼續(xù)下降45.2e-11收斂這個下降模式是典型的牛頓法特征。如果你跑出來第一次迭代范數(shù)反而升了多半是雅可比符號錯了或者修正量沒有乘上正確的方向。我特意對比了一下忽略損耗b c 0的情況迭代次數(shù)變成 3 次收斂但收斂后的直流電壓和考慮損耗時差了約 0.008 pu。這個差別在電壓等級上就是 80 V 左右對于某些對電壓精度要求高的場景是不能忽略的。6.3 結(jié)果校驗的三種手段第一種是功率平衡校驗。全系統(tǒng)的總有功注入應(yīng)該等于總負荷加總損耗交直流兩側(cè)分別算誤差應(yīng)該在 1e-6 以內(nèi)。這個方法能抓出大部分符號錯誤和單位錯誤。第二種是對拍。我用一個最簡單的兩節(jié)點交直流系統(tǒng)一個換流站、一條直流線路把方程手工展開成二元方程組用 fsolve 解一遍和自己的牛頓法結(jié)果對比。fsolve 不依賴雅可比結(jié)果可信度高。這個方法有點笨但對新模型非常有效。第三種是擾動測試。把某個負荷增加 1%看直流電壓的變化方向是否符合物理直覺。逆變側(cè)負荷增加直流電壓應(yīng)該下降因為要從直流網(wǎng)絡(luò)多取功率線路壓降增大。如果方向反了說明某個偏導(dǎo)符號錯了哪怕收斂結(jié)果看起來挺合理。提示三種方法里我覺得擾動測試最有價值。功率平衡校驗只能抓失衡抓不出符號錯誤對拍要寫兩套代碼成本高擾動測試幾行代碼就能做而且直接檢驗物理正確性。7. 常見問題與排查實錄7.1 不收斂的六類原因速查現(xiàn)象可能原因排查方法處理第一次迭代范數(shù)就變大雅可比符號反了數(shù)值差分對拍逐個塊檢查范數(shù)卡在 1e-3 不動耦合塊漏了或為零打印耦合塊元素補上損耗偏導(dǎo)迭代振蕩不收斂下垂系數(shù)過大減小 kdroop 試或改用主從控制矩陣奇異報錯定電壓節(jié)點未剔除檢查未知量索引剔除對應(yīng)行列收斂但功率不平衡標幺基準不統(tǒng)一檢查 G 矩陣量級統(tǒng)一轉(zhuǎn)標幺收斂到明顯錯誤解換流站符號約定混亂檢查 Pc/Pdc 正方向全程統(tǒng)一約定這六類我基本都踩過。印象最深的是第三類在一個五端直流算例上下垂系數(shù)取 20 時怎么都不收斂取 5 就正常。當時的直覺是系數(shù)大響應(yīng)快應(yīng)該更好收斂實際上下垂系數(shù)大意味著直流電壓對功率變化更敏感雅可比對角占優(yōu)性變差牛頓法反而更容易跳出收斂域。這也是為什么實際工程里下垂系數(shù)都要留裕度。7.2 雅可比奇異與病態(tài)除了定電壓節(jié)點不剔除導(dǎo)致的六異還有一種更隱蔽的情況兩個換流站接在同一個交流節(jié)點上。這時交流側(cè)注入功率是兩者之和雅可比對應(yīng)行的偏導(dǎo)也是兩者疊加如果你的 conv_id 映射寫錯把兩個站當成一個處理矩陣會出現(xiàn)重復(fù)行奇異但不報錯Matlab 會給警告但很多人不看警告。病態(tài)則多半來自量級差異。交流側(cè)的偏導(dǎo)量級在 1e0 到 1e1直流側(cè)的電導(dǎo)偏導(dǎo)在 1e1 到 1e2耦合塊在 1e-3。這種跨四個數(shù)量級的矩陣條件數(shù)很容易上 1e6雖然稀疏求解器能處理但精度會下降。我的做法是迭代前對雅可比做一次行列均衡Matlab 的equilibrate或者手動縮放實測能把條件數(shù)降一到兩個數(shù)量級。7.3 結(jié)果看著對但不對的隱性錯誤最麻煩的不是不收斂而是收斂到一個錯誤解。我遇到過兩次一次是 θ_ij 的定義搞反交流潮流的相角全部反號但因為系統(tǒng)對稱功率分布看起來完全正常直到我檢查某條線路的功率方向才發(fā)現(xiàn)不對。另一次是直流網(wǎng)絡(luò)的 G 矩陣沒有把對角元置為負的非對角元之和導(dǎo)致直流注入功率算出來偏大一倍。這個問題在純直流測試時沒暴露因為測試用例只有一條線路一接上交直流系統(tǒng)收斂后的直流電壓全部偏低我才回頭查 G 矩陣的組裝。預(yù)防這類問題的辦法只有一個每個模塊寫完就單獨測別等整體跑通再測。Ybus 組完先跑個純交流潮流G 矩陣組完先跑個純直流網(wǎng)絡(luò)都對上了再拼起來。拼起來之后再用手算小算例對拍一次。7.4 性能與工程化的幾點建議節(jié)點數(shù)上到幾百之后性能瓶頸會從算法轉(zhuǎn)移到實現(xiàn)。幾條經(jīng)驗Ybus 和 G 矩陣必須稀疏裝配用三元組一次性構(gòu)造別用循環(huán)逐元素賦值。不平衡量計算里避免 for 循環(huán)遍歷換流站改成向量化索引換流站數(shù)量少的時候差別不大上百之后能差幾倍。雅可比裝配是熱點能復(fù)用的部分比如交流子塊的結(jié)構(gòu)在迭代間不變可以只更新數(shù)值不重建稀疏結(jié)構(gòu)。如果要批量掃算例比如 N-1 分析把外層循環(huán)改成 parfor注意每個 worker 內(nèi)部別共享大矩陣內(nèi)存會爆。收斂判據(jù)不要用相對誤差標幺制下絕對誤差更直觀1e-8 對標幺量來說已經(jīng)足夠嚴格。8. 我個人踩過之后想說的幾句這套代碼我前后重構(gòu)過三次。第一次是照著交替迭代的思路硬改成統(tǒng)一迭代結(jié)果雅可比裝錯了收斂是收斂了但結(jié)果和交替迭代對不上查了兩天才發(fā)現(xiàn)耦合塊的符號反了。第二次是把換流站方程重推了一遍用數(shù)值差分逐塊驗證這次對了但代碼寫得極其啰嗦每次加新控制方式都要改五六個地方。第三次才想明白控制方式的變化本質(zhì)上是哪條方程換成哪個形式于是把方程抽象成可配置的模板主循環(huán)一行沒動擴展新控制方式只加一個分支。如果讓我給剛開始做這塊的人一條建議那就是先把兩節(jié)點交直流系統(tǒng)做到絕對正確再往上加節(jié)點。兩節(jié)點系統(tǒng)里每一個偏導(dǎo)都能手算每一個符號錯誤都會立刻表現(xiàn)為不收斂或者明顯偏離這時候調(diào)試成本最低。等系統(tǒng)大了錯誤會被網(wǎng)絡(luò)的平均效應(yīng)掩蓋你看到的只是收斂慢或者結(jié)果差一點點根本定位不到問題。另外別迷信商用軟件的結(jié)果。我用某商業(yè)仿真軟件對拍過一個算例兩邊差了 0.3%查了半天發(fā)現(xiàn)是損耗模型不同——它默認用恒定效率我用的是二次損耗。這種差異不是 bug但如果你不知道就會以為是自己的程序錯了然后去改一個本來正確的公式。模型對齊比結(jié)果對齊重要。最后說個小事Matlab 的稀疏左除在某些版本上對近奇異矩陣會給 warning 但不報錯如果你的循環(huán)里沒有檢查lastwarn很可能跑完一整批算例都不知道其中有幾個是錯的。我現(xiàn)在的做法是每次求解后強制檢查rank(full(J))小系統(tǒng)或者條件數(shù)估計超閾值就直接標記該算例為可疑不讓它污染統(tǒng)計結(jié)果。這個習慣幫我省了很多返工。 基于Matlab的交直流潮流統(tǒng)一迭代法建模與實現(xiàn)
1. 交直流潮流為什么最后落在了統(tǒng)一迭代上第一次把交直流混合系統(tǒng)往 Matlab 里塞的時候我的做法和大多數(shù)人一樣交流側(cè)用一套牛頓法算完把換流站交流側(cè)的有功當已知量丟給直流網(wǎng)絡(luò)直流網(wǎng)絡(luò)解出來再回代更新?lián)Q流站功率來回往復(fù)直到兩次迭代的換流站功率差小于閾值。這套交替迭代法也有人叫它順序求解法跑純主從控制的簡單算例時確實舒服代碼短、調(diào)試直觀、交流子程序和直流子程序各自獨立出問題一眼就知道是哪半邊。但只要換流站數(shù)量一多、控制方式一混搭、再加上換流器損耗隨電壓波動這種細枝末節(jié)交替迭代就開始露怯迭代次數(shù)從七八次漲到三四十次偶爾還會在兩個解之間來回振蕩你把收斂閾值放寬一點它就收斂了檢查結(jié)果發(fā)現(xiàn)直流電壓偏離參考值上千分之幾功率平衡差得離譜?；?Matlab 的交直流潮流計算統(tǒng)一迭代法建模與實現(xiàn)本質(zhì)上就是在講一件事把交流節(jié)點方程和直流節(jié)點方程塞進同一個牛頓迭代框架里用一個大雅可比矩陣同時求解用額外的矩陣維度換取收斂的魯棒性和解的唯一性。它不算什么高深算法真正難的地方全在建模細節(jié)上——換流站的穩(wěn)態(tài)方程怎么擺、標幺基準怎么統(tǒng)一、耦合偏導(dǎo)要不要留、直流松弛節(jié)點怎么處理。這些細節(jié)一旦搞錯程序照樣能跑出數(shù)來但結(jié)果是錯的而且錯得很隱蔽。這篇東西適合三類人看第一類是在做交直流混合系統(tǒng)相關(guān)課題、需要自己寫潮流程序而不是只會點軟件按鈕的第二類是想搞清楚為什么我的交替迭代不收斂的第三類是有一定 Matlab 基礎(chǔ)、能看懂矩陣運算但沒接觸過電力系統(tǒng)潮流內(nèi)部結(jié)構(gòu)的。我會把方程怎么來、雅可比怎么裝、代碼怎么寫、坑在哪里踩一條線講下來。默認你會基礎(chǔ)的交流牛頓-拉夫遜潮流如果連 Ybus 怎么組、H 矩陣怎么來的都不清楚建議先把純交流版本手推一遍再往下看不然讀到雅可比那節(jié)會很痛苦。需要提前說明一點下面給出的算例參數(shù)和迭代記錄是我為了講清楚流程自己搭的一個小系統(tǒng)不是任何實際工程數(shù)據(jù)你在復(fù)現(xiàn)時數(shù)值對不上很正常重點看結(jié)構(gòu)和邏輯。2. 方案選型交替迭代與統(tǒng)一迭代到底差在哪2.1 交替迭代法在工程里卡在什么地方交替迭代的思路非常工程師把大問題切成兩個能各自求解的小問題輪流解用邊界量互相傳遞。它的收斂性依賴一個隱含假設(shè)——交流側(cè)和直流側(cè)的耦合不能太強。主從控制下定直流電壓的那個換流站幾乎把兩側(cè)釘死了耦合確實弱交替兩三輪就穩(wěn)。但一旦上 droop 控制直流電壓本身參與功率分配換流站有功不再由交流側(cè)給定而是由直流電壓決定此時交流方程里出現(xiàn)未知量 Vdc直流方程里出現(xiàn)未知量 Vm兩個子系統(tǒng)互相牽著走交替迭代的收斂速度會肉眼可見地掉下來。我在一個五端直流的算例上做過對比全主從控制時交替迭代 6 次收斂閾值 1e-8把其中兩端改成下垂控制后交替迭代需要 31 次而且下垂系數(shù)取大一點就發(fā)散。這不是算法實現(xiàn)有 bug而是交替迭代本質(zhì)上是一個塊高斯-賽德爾迭代耦合項的譜半徑隨耦合強度增大而逼近 1收斂自然就慢。統(tǒng)一迭代相當于把它換成了整體牛頓二次收斂的性質(zhì)量級上不一樣。還有一個更現(xiàn)實的問題交替迭代的收斂判據(jù)不好定。你按換流站功率差判斷可能出現(xiàn)功率差已經(jīng)很小但電壓還在漂的情況你按電壓差判斷又可能在功率嚴重不平衡時誤判收斂。統(tǒng)一迭代沒有這個煩惱不平衡量向量統(tǒng)一范數(shù)一個判據(jù)管到底。2.2 統(tǒng)一迭代付出的代價值不值代價很直接雅可比維度變大裝配邏輯變復(fù)雜調(diào)試難度上升。一個 nb 個交流節(jié)點、nd 個直流節(jié)點的系統(tǒng)純交流雅可比是 (2nb-1) 階左右加上直流后大概多出 nd 階而且多出來的耦合塊元素推導(dǎo)容易出錯。更要命的是一旦某個耦合偏導(dǎo)符號寫反程序不一定發(fā)散而是收斂到一個物理上不成立的解你拿萬用表都量不出來。但從我自己的經(jīng)驗看這筆賬是劃算的。第一牛頓法對初值敏感但一旦進入收斂域就是二次收斂交直流系統(tǒng)初值可以很方便地用平啟動 直流電壓全置 1.0給出實際很少不收斂第二統(tǒng)一迭代天然支持各種控制方式混搭主從、下垂、定功率、定電壓只需要改對應(yīng)行的方程框架不動第三擴展性好以后要加換流器限流、交流側(cè)無功限值、直流線路損耗都只是在已有方程上加項不像交替迭代那樣要重新設(shè)計外層循環(huán)。代價還有一條容易被忽略統(tǒng)一迭代對模型的完整性要求高。交替迭代可以容忍某一側(cè)用簡化模型比如直流側(cè)只當恒功率負載統(tǒng)一迭代不行你寫進雅可比的每一項都必須和方程嚴格對應(yīng)否則收斂性立刻崩。這其實是好事逼著我把換流器方程老老實實推了一遍。2.3 什么情況下我會勸你別用統(tǒng)一迭代不是所有場景都值得上統(tǒng)一迭代。如果系統(tǒng)里只有一個背靠背換流站兩側(cè)交流網(wǎng)絡(luò)規(guī)模都很小交替迭代兩三次就收斂代碼量還只有統(tǒng)一迭代的三分之一那就沒必要折騰。如果直流側(cè)是完全理想的恒功率負載模型損耗忽略、控制理想交直流之間嚴格解耦統(tǒng)一迭代的雅可比是分塊對角的等于把兩個獨立問題拼在一起解白白浪費計算量。真正值得上統(tǒng)一迭代的典型場景是多端直流、下垂控制、換流站數(shù)量大于等于三、需要嚴格校驗功率平衡、后續(xù)要做 N-1 掃描或連續(xù)潮流。這種情況下交替迭代的迭代次數(shù)不穩(wěn)定批量計算時總有幾個算例卡住統(tǒng)一迭代的穩(wěn)定收斂反而是省時間的。對比維度交替迭代法統(tǒng)一迭代法雅可比規(guī)模兩個小矩陣分別求解一個大矩陣一次求解收斂階數(shù)線性塊高斯-賽德爾二次牛頓控制方式混搭需要改外層調(diào)度邏輯只改對應(yīng)行方程下垂控制支持收斂慢易振蕩與主從控制同框架調(diào)試難度低兩側(cè)獨立高耦合項易錯代碼量我的實現(xiàn)約 300 行約 600 行擴展新設(shè)備改動大加行加列即可3. 模型拆解把換流站寫成能進雅可比的方程3.1 節(jié)點分類與未知量編號統(tǒng)一迭代第一個要想清楚的問題不是方程而是誰是未知量。交流側(cè)沿用經(jīng)典分類PQ 節(jié)點未知量是電壓幅值 Vm 和相角 VaPV 節(jié)點未知量只有相角幅值由控制給定平衡節(jié)點兩個都給定。直流側(cè)稍微特殊一點直流網(wǎng)絡(luò)沒有相角概念只有電壓幅值 Vdc所以每個直流節(jié)點的未知量就是 Vdc 本身但定直流電壓的那個節(jié)點例外——它的 Vdc 已知功率未知等價于交流側(cè)的平衡節(jié)點。這個例外是很多人第一次寫統(tǒng)一迭代時栽跟頭的地方。直覺上你會覺得直流節(jié)點都應(yīng)該有功率方程于是給定了電壓的節(jié)點也列了一條 ΔPdc 方程結(jié)果雅可比對應(yīng)行是零行或者接近零行矩陣直接奇異。正確做法是定電壓節(jié)點的 Vdc 從未知量集合里剔除對應(yīng)的 ΔPdc 方程也不列它的功率在收斂后由 Pdc Vdc * Σ G Vdc 反算出來。我用的是未知量索引映射的辦法先掃一遍表把所有未知量編號再用映射數(shù)組把物理節(jié)點號翻譯成未知量序號。這樣后面裝配雅可比時不用到處判斷類型索引一查就行。具體做法在第四節(jié)給代碼。3.2 換流站的穩(wěn)態(tài)方程怎么擺換流站建模有兩條路平均模型和開關(guān)模型。潮流計算用平均模型忽略開關(guān)過程和 PWM 諧波只保留基波功率平衡。我采用的是最常見的功率平衡寫法交流側(cè)注入電網(wǎng)的有功 Pc直流側(cè)從直流網(wǎng)絡(luò)取用的有功 Pdc。兩者之差是換流器損耗Pdc Pc Ploss。交流側(cè)無功 Qc 由控制方式?jīng)Q定定 Q 控制時直接給定。損耗模型取 Ploss a b·I c·I2其中 I 是換流站交流側(cè)電流的標幺值I sqrt(Pc2 Qc2) / Vm。選這個損耗模型不是因為它精確而是因為它同時包含了常數(shù)項、線性項和二次項能覆蓋導(dǎo)通損耗和開關(guān)損耗的主要成分而且對 Vm 的偏導(dǎo)好推——這一點在裝配耦合塊時非常關(guān)鍵。有些文獻用恒定效率模型Pdc Pc / η簡單但那樣 Ploss 對 Vm 的偏導(dǎo)恒為零交直流之間就沒有耦合項了統(tǒng)一迭代退化成兩個獨立問題失去了演示意義。注意損耗模型的選擇直接決定耦合塊是否存在。如果你的算例里所有換流站都用恒定效率、所有控制都是定 P 定 Q那統(tǒng)一迭代的雅可比是分塊對角的收斂性和交替迭代一樣這時候沒必要上統(tǒng)一迭代用它只是自己給自己找麻煩。3.3 標幺基準怎么統(tǒng)一交直流混合系統(tǒng)最容易出錯的地方是基準值。交流側(cè)一般取三相功率基準 Sb 和線電壓基準 Vb直流側(cè)如果直接沿用交流的 Sb直流電壓基準要單獨取——我習慣取 Vdb 使得直流側(cè)基準阻抗與交流側(cè)一致也就是 Vdb2 / Sb Vb2 / Sb這樣 Vdb Vb這里的 Vb 是交流相電壓基準若用線電壓基準則要乘系數(shù)。如果不統(tǒng)一直流網(wǎng)絡(luò)的電導(dǎo)矩陣 G 數(shù)值會偏離好幾個數(shù)量級雅可比條件數(shù)爆炸。我吃過一次虧直流側(cè)電導(dǎo)一開始用的是實際歐姆值沒轉(zhuǎn)標幺結(jié)果雅可比對角元量級從 1e1 跳到 1e4稀疏求解器雖然沒報錯但迭代過程出現(xiàn)了明顯的數(shù)值抖動收斂判據(jù)在 1e-6 附近停住不動。改成標幺后一次就收斂了。所以我的建議是建模階段就把所有量統(tǒng)一到同一套基準下代碼里只出現(xiàn)標幺值物理量只在輸入輸出環(huán)節(jié)轉(zhuǎn)換。量基準轉(zhuǎn)換關(guān)系常見錯誤交流功率SbP_pu P / Sb單相/三相混淆交流電壓Vb線電壓Vm_pu V / Vb用相電壓基準直流電壓Vdb VbVdc_pu Vdc / Vdb直接用電網(wǎng)額定值直流電阻Zdb Vdb2/SbR_pu R / Zdb忘了轉(zhuǎn)標幺換流器電流Ib Sb/(√3·Vb)I_pu I / Ib與功率混淆損耗系數(shù)按 Sb 歸一a_pu a / Sb系數(shù)單位不統(tǒng)一3.4 完整不平衡量向量的構(gòu)成把所有方程按順序排好不平衡量向量 ΔF 分成三塊第一塊是交流有功不平衡 ΔP對所有非平衡交流節(jié)點列共 (nb - 1) 條。對純負荷節(jié)點ΔP_i P_i^calc - P_i^spec對掛換流站的節(jié)點注入功率里要減去換流站交流側(cè)取用的 Pc。第二塊是交流無功不平衡 ΔQ對所有 PQ 節(jié)點列。掛換流站的節(jié)點如果換流站定 Q 控制那 Qc 作為給定值進 P_spec 的位置。第三塊是直流功率不平衡 ΔPdc對所有非定電壓直流節(jié)點列ΔPdc_k Vdc_k·Σ_m G_km·Vdc_m - Pdc_k^spec。這里有個符號約定的坑必須講清楚。換流站整流交流到直流時 Pc 為正、Pdc 為負還是反過來不同文獻寫法不一樣。我統(tǒng)一約定Pc 為換流站從交流電網(wǎng)吸收的有功正表示吸收Pdc 為換流站注入直流網(wǎng)絡(luò)的有功正表示注入那么整流時 Pc 0、Pdc 0 的絕對值小于 Pc關(guān)系是 Pdc -(Pc - Ploss)。逆變時反過來。約定一旦定下來全套代碼都要遵守中間改一次符號所有耦合項都得重推。4. 雅可比矩陣的結(jié)構(gòu)與裝配細節(jié)4.1 分塊結(jié)構(gòu)與稀疏索引統(tǒng)一迭代的雅可比是一個三行三列的分塊矩陣J [H N Ndc; M L Ldc; Md Hd D ]H、N、M、L 是大家熟悉的交流子塊分別是 ΔP 對 Va、ΔP 對 Vm、ΔQ 對 Va、ΔQ 對 Vm 的偏導(dǎo)。D 是直流子塊ΔPdc 對 Vdc 的偏導(dǎo)。剩下四塊是耦合塊Ndc 是 ΔP 對 VdcLdc 是 ΔQ 對 VdcMd 是 ΔPdc 對 VaHd 是 ΔPdc 對 Vm。在我采用的換流站模型下Ndc 和 Ldc 恒為零——因為交流側(cè)注入功率 Pc、Qc 是控制給定的不直接依賴直流電壓。Md 也是零——直流功率不依賴交流相角。只有 Hd 非零它來自損耗項對交流電壓幅值的依賴。這個結(jié)論很重要它說明統(tǒng)一迭代的耦合其實只有一個塊在起作用如果你把這個塊也簡化掉損耗恒定矩陣就是分塊下三角可以用兩次回代求解本質(zhì)上又退回了交替迭代。所以我的原則是耦合塊哪怕數(shù)值很小也必須保留它是統(tǒng)一迭代區(qū)別于交替迭代的物理來源。實際算例里 Hd 的元素量級大概在 1e-3 到 1e-2確實小但正是這個小項在迭代初期提供了方向信息減少了來回振蕩。4.2 交流子塊的經(jīng)典寫法回顧交流部分沿用標準極坐標牛頓法公式我不重復(fù)推導(dǎo)只強調(diào)兩個實現(xiàn)細節(jié)。第一Ybus 要用稀疏矩陣Matlab 里 sparse 存儲裝配功率時用 Ybus 乘以電壓向量別用循環(huán)。第二對角元和非對角元的偏導(dǎo)公式不同裝配時按行列分別處理不要圖省事用一個循環(huán)套所有元素。H 的對角元H_ii -Q_i - B_ii·Vm_i2 H 的非對角元H_ij Vm_i·Vm_j·(G_ij·sin θ_ij - B_ij·cos θ_ij)N 的對角元N_ii P_i/Vm_i G_ii·Vm_i N 的非對角元N_ij Vm_i·(G_ij·cos θ_ij B_ij·sin θ_ij)M、L 類似符號別搞反。這部分是教科書內(nèi)容但我在實際 debug 中發(fā)現(xiàn)很多人栽在 θ_ij 的定義上——θ_ij θ_i - θ_j 還是 θ_j - θ_i兩種寫法都有人用混用會導(dǎo)致 H 的非對角元整體符號反了程序仍然收斂但結(jié)果完全錯誤。我的做法是全程用 θ_ij θ_i - θ_j并且寫一個單元測試兩節(jié)點系統(tǒng)手算一遍對拍。4.3 直流子塊與耦合子塊的偏導(dǎo)推導(dǎo)直流注入功率 Pdc_k Vdc_k · Σ_m G_km · Vdc_m。對 Vdc_m 求偏導(dǎo)m ≠ k 時?Pdc_k/?Vdc_m Vdc_k · G_kmm k 時?Pdc_k/?Vdc_k Σ_m G_km·Vdc_m Vdc_k·G_kk這個對角元很容易寫錯成 2·Vdc_k·G_kk或者漏掉求和項里的 k 項。正確形式是該節(jié)點注入電流 自身電導(dǎo)乘電壓物理上就是電流對電壓的靈敏度。耦合塊 Hd 的推導(dǎo)ΔPdc_k 里包含 Pdc_k^spec而 Pdc_k^spec Pc_k Ploss_k按前面的符號約定取絕對值Ploss_k a b·I c·I2I sqrt(Pc2 Qc2)/Vm_iVm_i 是換流站所接交流節(jié)點的電壓幅值。?Ploss/?Vm (b 2c·I) · ?I/?Vm (b 2c·I) · (-I/Vm)所以 ?ΔPdc_k/?Vm_i -?Ploss/?Vm_i (b 2c·I)·I/Vm_i。符號為正物理意義是電壓幅值升高、電流減小、損耗減小、直流側(cè)所需注入功率增加合理。提示這個偏導(dǎo)的量級很小如果你在調(diào)試時發(fā)現(xiàn)它一直是零先檢查是不是損耗系數(shù) b、c 都設(shè)成了零或者 Vm 初值取 1.0 導(dǎo)致靈敏度看起來沒變化?？梢杂脭?shù)值差分驗證一遍方法在下一節(jié)。4.4 用數(shù)值差分給雅可比做體檢手推偏導(dǎo)最怕符號錯。我的習慣是寫完解析雅可比后立刻用中心差分對拍一次誤差超過 1e-4 就說明有地方錯了。代碼不到二十行但能省掉幾個晚上的調(diào)試function Jnum jac_check(x, Ffun, nx, h) % x : 當前狀態(tài)向量 % Ffun: 不平衡量函數(shù)句柄 F Ffun(x) % nx : 未知量個數(shù) % h : 差分步長取 1e-6 左右 F0 Ffun(x); Jnum zeros(numel(F0), nx); for k 1:nx xp x; xp(k) xp(k) h; xm x; xm(k) xm(k) - h; Jnum(:, k) (Ffun(xp) - Ffun(xm)) / (2*h); end end用法是先在第一次迭代時調(diào)用把解析雅可比 J 和 Jnum 相減取絕對值最大元素打印出來。如果最大值出現(xiàn)在耦合塊對應(yīng)的行說明耦合偏導(dǎo)有問題如果出現(xiàn)在交流子塊說明 Ybus 或者 θ 約定有誤。差分步長不宜太小1e-6 到 1e-8 之間試太大會有截斷誤差太小會撞上雙精度舍入誤差。5. Matlab 實現(xiàn)從數(shù)據(jù)結(jié)構(gòu)到求解循環(huán)5.1 三張表的數(shù)據(jù)結(jié)構(gòu)我不用面向?qū)ο笠膊挥?struct 數(shù)組就用三個二維矩陣列定義寫在注釋里。好處是索引快、易于向量化、保存成 mat 文件方便。缺點是列含義靠記憶所以我把列定義固化成常量。% bus : 交流節(jié)點表每行一個節(jié)點 % 1 id | 2 type(1PQ,2PV,3平衡) | 3 Vm | 4 Va(rad) % 5 Pg | 6 Qg | 7 Pd | 8 Qd | 9 conv_id(0無換流站) % % dc : 直流節(jié)點表 % 1 id | 2 Vdc | 3 type(1定電壓,2定功率,3下垂) | 4 Pdc_spec % 5 kdroop | 6 Vdc_ref | 7 conv_id % % conv : 換流站表 % 1 id | 2 ac_bus | 3 dc_node | 4 mode(1定P定Q,2定Vdc) % 5 Pset | 6 Qset | 7 a | 8 b | 9 ctype字段我堅持用整數(shù)而不是字符串因為后面要大量做邏輯判斷和索引字符串比較慢且容易拼錯。conv_id這個外鍵放在 bus 和 dc 表里是雙向索引的關(guān)鍵裝配耦合塊時靠它從交流節(jié)點找到對應(yīng)的直流節(jié)點反之亦然。5.2 未知量映射與主循環(huán)骨架主循環(huán)的結(jié)構(gòu)固定算不平衡量 → 判斷收斂 → 裝雅可比 → 解修正方程 → 更新狀態(tài)。真正需要小心的是狀態(tài)更新的邊界處理以及修正量對應(yīng)的物理量。function res acdc_nr(bus, dc, conv, opt) nb size(bus,1); nd size(dc,1); % 交流未知量索引 iVa find(bus(:,2) ~ 3); % 非平衡節(jié)點的相角 iVm find(bus(:,2) 1); % PQ 節(jié)點的幅值 % 直流未知量索引剔除定電壓節(jié)點 iVd find(dc(:,3) ~ 1); nx numel(iVa) numel(iVm) numel(iVd); % 狀態(tài)向量 x 與物理表的映射 mapVa zeros(nb,1); mapVa(iVa) 1:numel(iVa); mapVm zeros(nb,1); mapVm(iVm) numel(iVa) (1:numel(iVm)); mapVd zeros(nd,1); mapVd(iVd) numel(iVa)numel(iVm) (1:numel(iVd)); Vm bus(:,3); Va bus(:,4); Vdc dc(:,2); Y make_ybus(bus); G make_gdc(dc); for it 1:opt.maxit [dP, dQ] ac_mismatch(Y, Vm, Va, bus, conv, Vdc); dPdc dc_mismatch(G, Vdc, dc, conv, Vm); F [dP(iVa); dQ(iVm); dPdc(iVd)]; if norm(F, inf) opt.tol break; end J build_jac(Y, Vm, Va, G, Vdc, bus, dc, conv, ... iVa, iVm, iVd, mapVa, mapVm, mapVd); dx J \ (-F); % 稀疏左除 % 回代更新 Va(iVa) Va(iVa) dx(mapVa(iVa)); Vm(iVm) Vm(iVm) dx(mapVm(iVm)); Vdc(iVd) Vdc(iVd) dx(mapVd(iVd)); end res pack_result(bus, dc, conv, Vm, Va, Vdc, it, norm(F,inf)); end這段骨架里有個細節(jié)值得說J \ (-F)用的是左除Matlab 會自動識別 J 是稀疏矩陣并調(diào)用合適的求解器。如果你寫inv(J)*(-F)一是慢二是精度差三是在大系統(tǒng)上會直接內(nèi)存溢出。另外norm(F, inf)取無窮范數(shù)比二范數(shù)更嚴格因為它盯的是單個最大的不平衡量不會被平均掉。5.3 交流不平衡量的計算function [dP, dQ] ac_mismatch(Y, Vm, Va, bus, conv, Vdc) V Vm .* exp(1j*Va); S V .* conj(Y * V); P real(S); Q imag(S); Ps (bus(:,5) - bus(:,7)) / 1.0; % 給定注入已標幺 Qs (bus(:,6) - bus(:,8)) / 1.0; % 換流站交流側(cè)取用功率計入負荷 for k 1:size(conv,1) i conv(k,2); c k; [Pc, Qc, ~] conv_ac_power(conv(c,:), Vm(i), Vdc(conv(c,3))); Ps(i) Ps(i) - Pc; Qs(i) Qs(i) - Qc; end dP P - Ps; dQ Q - Qs; end這里conv_ac_power是核心它要根據(jù)控制方式返回換流站交流側(cè)實際取用的有功和無功。定 P 定 Q 模式下 Pc Pset、Qc Qset直接返回定直流電壓模式下Pc 是由直流側(cè)功率反推的未知量不能在這里給死需要靠方程間的關(guān)系解決——這也是統(tǒng)一迭代相比交替迭代省事的地方未知量 Pc 不出現(xiàn)在狀態(tài)向量里但它通過直流方程隱含求解。嚴格來說定電壓站的 Pc 在每次迭代時應(yīng)該由當前 Vdc 和損耗反算Pc Pdc - Ploss而 Pdc Vdc·ΣG·Vdc。這就是所謂的內(nèi)層更新不影響雅可比結(jié)構(gòu)但影響收斂速度。我一開始沒做這個內(nèi)層更新直接用上一次迭代的 Pc 值結(jié)果迭代次數(shù)多了快一倍。5.4 直流不平衡量與稀疏裝配function dPdc dc_mismatch(G, Vdc, dc, conv, Vm) Idc G * Vdc; Pdc Vdc .* Idc; Pspec zeros(size(Vdc)); for k 1:size(conv,1) j conv(k,3); [~, ~, Pdc_k] conv_ac_power(conv(k,:), Vm(conv(k,2)), Vdc(j)); Pspec(j) Pdc_k; end dPdc Pdc - Pspec; end雅可比裝配我采用三元組預(yù)分配的方式先用spalloc估一個非零元上界用三個長向量 irow、jcol、jval 累積最后一次性sparse。這種方式比逐元素賦值快一個數(shù)量級在節(jié)點數(shù)上百之后差別非常明顯。nnz_est 4*(nb nd) 16*size(conv,1); irow zeros(nnz_est,1); jcol irow; jval irow; cnt 0; add (r,c,v) deal_helper(); % 這里用局部函數(shù)累加略實際寫的時候我把累加封裝成局部函數(shù)避免匿名函數(shù)帶來的性能損失function push(r, c, v) cnt cnt 1; irow(cnt) r; jcol(cnt) c; jval(cnt) v; end這種寫法在嵌套函數(shù)里能直接改外部變量比傳結(jié)構(gòu)體清爽。最后J sparse(irow(1:cnt), jcol(1:cnt), jval(1:cnt), nx, nx)。5.5 收斂判據(jù)與迭代控制判據(jù)我一般用兩層主判據(jù)是無窮范數(shù)小于 1e-8標幺輔助判據(jù)是修正量范數(shù)小于 1e-6兩個都滿足才認為收斂。只判不平衡量的話在病態(tài)系統(tǒng)上可能出現(xiàn)不平衡量很小但電壓還在緩慢漂移的情況。迭代次數(shù)上限設(shè) 30超過就認為不收斂打印最后一次的不平衡量分布看是哪個節(jié)點貢獻最大。這個診斷信息非常有用——如果最大不平衡量集中在某個換流站節(jié)點八成是控制方式設(shè)置和方程不匹配如果分散在所有節(jié)點八成是初值或者標幺化的問題。還有一個實用技巧記錄每次迭代的不平衡量范數(shù)畫成半對數(shù)圖。正常牛頓法應(yīng)該是一條陡峭的直線下降如果出現(xiàn)平臺期、抖動或者先降后升說明雅可比有問題或者系統(tǒng)本身無解。這個圖我?guī)缀趺看握{(diào)試新模型都會畫比看數(shù)字直觀得多。6. 一個可復(fù)現(xiàn)的小算例與結(jié)果校驗6.1 算例配置我搭了一個五節(jié)點交流 三節(jié)點直流的系統(tǒng)規(guī)模小到能手算校驗又能體現(xiàn)出交直流耦合的特征。交流側(cè)節(jié)點 1 是平衡節(jié)點節(jié)點 2、3 是純負荷節(jié)點 4、5 各接一個換流站。直流側(cè)節(jié)點 1 是定電壓站對應(yīng)交流節(jié)點 4 的換流站節(jié)點 2、3 是定功率站對應(yīng)交流節(jié)點 5 的換流站和一個純直流負荷。參數(shù)數(shù)值說明交流基準功率100 MVA全系統(tǒng)統(tǒng)一交流基準電壓10 kV線電壓交流線路阻抗0.01 j0.08 pu環(huán)形結(jié)構(gòu)負荷節(jié)點 20.6 j0.2 pu恒功率負荷節(jié)點 30.4 j0.15 pu恒功率換流站 1 控制定 Vdc 1.0主從控制的從站換流站 2 控制定 P -0.3 pu逆變運行直流線路電阻0.02 pu三節(jié)點三角形損耗系數(shù) a0.002 pu常數(shù)項損耗系數(shù) b0.01線性項損耗系數(shù) c0.02二次項需要注意的是這里的定功率是站在直流側(cè)說的換流站 2 從直流網(wǎng)絡(luò)取 0.3 pu 的功率送到交流側(cè)同時因為損耗它從直流網(wǎng)絡(luò)實際取用的功率要略大于 0.3。6.2 迭代過程記錄平啟動所有電壓置 1.0相角置 0。前三次迭代的不平衡量范數(shù)大致是迭代次數(shù)不平衡量范數(shù)說明06.2e-1初始不平衡主要是負荷節(jié)點18.7e-2大幅下降牛頓方向正確24.1e-3進入二次收斂區(qū)39.6e-6繼續(xù)下降45.2e-11收斂這個下降模式是典型的牛頓法特征。如果你跑出來第一次迭代范數(shù)反而升了多半是雅可比符號錯了或者修正量沒有乘上正確的方向。我特意對比了一下忽略損耗b c 0的情況迭代次數(shù)變成 3 次收斂但收斂后的直流電壓和考慮損耗時差了約 0.008 pu。這個差別在電壓等級上就是 80 V 左右對于某些對電壓精度要求高的場景是不能忽略的。6.3 結(jié)果校驗的三種手段第一種是功率平衡校驗。全系統(tǒng)的總有功注入應(yīng)該等于總負荷加總損耗交直流兩側(cè)分別算誤差應(yīng)該在 1e-6 以內(nèi)。這個方法能抓出大部分符號錯誤和單位錯誤。第二種是對拍。我用一個最簡單的兩節(jié)點交直流系統(tǒng)一個換流站、一條直流線路把方程手工展開成二元方程組用 fsolve 解一遍和自己的牛頓法結(jié)果對比。fsolve 不依賴雅可比結(jié)果可信度高。這個方法有點笨但對新模型非常有效。第三種是擾動測試。把某個負荷增加 1%看直流電壓的變化方向是否符合物理直覺。逆變側(cè)負荷增加直流電壓應(yīng)該下降因為要從直流網(wǎng)絡(luò)多取功率線路壓降增大。如果方向反了說明某個偏導(dǎo)符號錯了哪怕收斂結(jié)果看起來挺合理。提示三種方法里我覺得擾動測試最有價值。功率平衡校驗只能抓失衡抓不出符號錯誤對拍要寫兩套代碼成本高擾動測試幾行代碼就能做而且直接檢驗物理正確性。7. 常見問題與排查實錄7.1 不收斂的六類原因速查現(xiàn)象可能原因排查方法處理第一次迭代范數(shù)就變大雅可比符號反了數(shù)值差分對拍逐個塊檢查范數(shù)卡在 1e-3 不動耦合塊漏了或為零打印耦合塊元素補上損耗偏導(dǎo)迭代振蕩不收斂下垂系數(shù)過大減小 kdroop 試或改用主從控制矩陣奇異報錯定電壓節(jié)點未剔除檢查未知量索引剔除對應(yīng)行列收斂但功率不平衡標幺基準不統(tǒng)一檢查 G 矩陣量級統(tǒng)一轉(zhuǎn)標幺收斂到明顯錯誤解換流站符號約定混亂檢查 Pc/Pdc 正方向全程統(tǒng)一約定這六類我基本都踩過。印象最深的是第三類在一個五端直流算例上下垂系數(shù)取 20 時怎么都不收斂取 5 就正常。當時的直覺是系數(shù)大響應(yīng)快應(yīng)該更好收斂實際上下垂系數(shù)大意味著直流電壓對功率變化更敏感雅可比對角占優(yōu)性變差牛頓法反而更容易跳出收斂域。這也是為什么實際工程里下垂系數(shù)都要留裕度。7.2 雅可比奇異與病態(tài)除了定電壓節(jié)點不剔除導(dǎo)致的六異還有一種更隱蔽的情況兩個換流站接在同一個交流節(jié)點上。這時交流側(cè)注入功率是兩者之和雅可比對應(yīng)行的偏導(dǎo)也是兩者疊加如果你的 conv_id 映射寫錯把兩個站當成一個處理矩陣會出現(xiàn)重復(fù)行奇異但不報錯Matlab 會給警告但很多人不看警告。病態(tài)則多半來自量級差異。交流側(cè)的偏導(dǎo)量級在 1e0 到 1e1直流側(cè)的電導(dǎo)偏導(dǎo)在 1e1 到 1e2耦合塊在 1e-3。這種跨四個數(shù)量級的矩陣條件數(shù)很容易上 1e6雖然稀疏求解器能處理但精度會下降。我的做法是迭代前對雅可比做一次行列均衡Matlab 的equilibrate或者手動縮放實測能把條件數(shù)降一到兩個數(shù)量級。7.3 結(jié)果看著對但不對的隱性錯誤最麻煩的不是不收斂而是收斂到一個錯誤解。我遇到過兩次一次是 θ_ij 的定義搞反交流潮流的相角全部反號但因為系統(tǒng)對稱功率分布看起來完全正常直到我檢查某條線路的功率方向才發(fā)現(xiàn)不對。另一次是直流網(wǎng)絡(luò)的 G 矩陣沒有把對角元置為負的非對角元之和導(dǎo)致直流注入功率算出來偏大一倍。這個問題在純直流測試時沒暴露因為測試用例只有一條線路一接上交直流系統(tǒng)收斂后的直流電壓全部偏低我才回頭查 G 矩陣的組裝。預(yù)防這類問題的辦法只有一個每個模塊寫完就單獨測別等整體跑通再測。Ybus 組完先跑個純交流潮流G 矩陣組完先跑個純直流網(wǎng)絡(luò)都對上了再拼起來。拼起來之后再用手算小算例對拍一次。7.4 性能與工程化的幾點建議節(jié)點數(shù)上到幾百之后性能瓶頸會從算法轉(zhuǎn)移到實現(xiàn)。幾條經(jīng)驗Ybus 和 G 矩陣必須稀疏裝配用三元組一次性構(gòu)造別用循環(huán)逐元素賦值。不平衡量計算里避免 for 循環(huán)遍歷換流站改成向量化索引換流站數(shù)量少的時候差別不大上百之后能差幾倍。雅可比裝配是熱點能復(fù)用的部分比如交流子塊的結(jié)構(gòu)在迭代間不變可以只更新數(shù)值不重建稀疏結(jié)構(gòu)。如果要批量掃算例比如 N-1 分析把外層循環(huán)改成 parfor注意每個 worker 內(nèi)部別共享大矩陣內(nèi)存會爆。收斂判據(jù)不要用相對誤差標幺制下絕對誤差更直觀1e-8 對標幺量來說已經(jīng)足夠嚴格。8. 我個人踩過之后想說的幾句這套代碼我前后重構(gòu)過三次。第一次是照著交替迭代的思路硬改成統(tǒng)一迭代結(jié)果雅可比裝錯了收斂是收斂了但結(jié)果和交替迭代對不上查了兩天才發(fā)現(xiàn)耦合塊的符號反了。第二次是把換流站方程重推了一遍用數(shù)值差分逐塊驗證這次對了但代碼寫得極其啰嗦每次加新控制方式都要改五六個地方。第三次才想明白控制方式的變化本質(zhì)上是哪條方程換成哪個形式于是把方程抽象成可配置的模板主循環(huán)一行沒動擴展新控制方式只加一個分支。如果讓我給剛開始做這塊的人一條建議那就是先把兩節(jié)點交直流系統(tǒng)做到絕對正確再往上加節(jié)點。兩節(jié)點系統(tǒng)里每一個偏導(dǎo)都能手算每一個符號錯誤都會立刻表現(xiàn)為不收斂或者明顯偏離這時候調(diào)試成本最低。等系統(tǒng)大了錯誤會被網(wǎng)絡(luò)的平均效應(yīng)掩蓋你看到的只是收斂慢或者結(jié)果差一點點根本定位不到問題。另外別迷信商用軟件的結(jié)果。我用某商業(yè)仿真軟件對拍過一個算例兩邊差了 0.3%查了半天發(fā)現(xiàn)是損耗模型不同——它默認用恒定效率我用的是二次損耗。這種差異不是 bug但如果你不知道就會以為是自己的程序錯了然后去改一個本來正確的公式。模型對齊比結(jié)果對齊重要。最后說個小事Matlab 的稀疏左除在某些版本上對近奇異矩陣會給 warning 但不報錯如果你的循環(huán)里沒有檢查lastwarn很可能跑完一整批算例都不知道其中有幾個是錯的。我現(xiàn)在的做法是每次求解后強制檢查rank(full(J))小系統(tǒng)或者條件數(shù)估計超閾值就直接標記該算例為可疑不讓它污染統(tǒng)計結(jié)果。這個習慣幫我省了很多返工。