交直流混合系統(tǒng)的潮流計(jì)算麻煩的從來不是交流那一半。交流部分的 PQ 分解、牛頓迭代、稀疏求解隨便找個(gè)開源庫(kù)都能跑。真正讓人抓頭的是直流側(cè)那幾個(gè)換流器它們既是交流系統(tǒng)的負(fù)荷又是電源還帶著觸發(fā)角、熄弧角、換相電抗這一堆相互糾纏的參數(shù)。你改一個(gè)直流電流交流側(cè)兩個(gè)節(jié)點(diǎn)的有功無(wú)功同時(shí)變直流電壓也跟著變最后連控制方式都可能從一個(gè)切到另一個(gè)。我最初做這套東西是為了給一個(gè)含兩端直流聯(lián)絡(luò)線的測(cè)試系統(tǒng)算穩(wěn)態(tài)工作點(diǎn)。一開始圖省事用了交替迭代交流算一遍、直流算一遍、來回倒騰。小擾動(dòng)工況下能收斂一旦把直流功率推到接近額定值或者在逆變側(cè)把熄弧角壓到 16 度附近迭代就開始來回振蕩殘差卡在 1e-3 上下死活下不去。后來老老實(shí)實(shí)改成統(tǒng)一迭代法把交流節(jié)點(diǎn)電壓、相角、直流電壓、直流電流、控制角全部塞進(jìn)同一個(gè)未知量向量用同一個(gè)牛頓步一起修正收斂性立刻不一樣了——典型工況 4 到 6 次迭代就到 1e-8。這篇就把這套基于 Matlab 的統(tǒng)一迭代法交直流潮流計(jì)算從頭到尾拆一遍模型怎么建、方程怎么配平、雅可比矩陣那四個(gè)塊怎么拼、代碼骨架長(zhǎng)什么樣、以及我在調(diào)試時(shí)踩過的那些坑。適合已經(jīng)會(huì)寫純交流牛頓法潮流、想往交直流方向走一步的讀者也適合手頭有個(gè)混合系統(tǒng)要算穩(wěn)態(tài)點(diǎn)、不想從零推公式的人直接抄結(jié)構(gòu)。1. 為什么交直流潮流要單獨(dú)搭一套求解器1.1 純交流潮流模型在換流站面前會(huì)直接失靈純交流潮流里節(jié)點(diǎn)只有三類PQ、PV、平衡節(jié)點(diǎn)。每個(gè)節(jié)點(diǎn)的注入功率是常數(shù)或者電壓幅值給定雅可比矩陣的元素都是電壓幅值和相角的函數(shù)形式干凈利落。換流站一進(jìn)來這個(gè)假設(shè)就崩了。換流站對(duì)交流系統(tǒng)來說有功和無(wú)功都不是定值。有功由直流側(cè)傳輸功率決定無(wú)功由換流器的功率因數(shù)決定而功率因數(shù)又跟觸發(fā)角、換相重疊角、直流電流掛鉤。這意味著你在交流節(jié)點(diǎn)上寫不出一個(gè)只含 V 和 θ 的功率方程——它里面必然夾著直流側(cè)的未知量。更直接點(diǎn)說換流節(jié)點(diǎn)是一個(gè)特性由另一套方程描述的負(fù)荷你沒法把它簡(jiǎn)化成恒功率或者恒阻抗。還有一層換流變壓器變比、換相電抗、每極換流器個(gè)數(shù)、直流接線方式單極、雙極、金屬回線這些結(jié)構(gòu)參數(shù)都會(huì)改變方程的形態(tài)。雙極系統(tǒng)里正負(fù)極的直流電流在理想對(duì)稱時(shí)相等但接地回路和金屬回線的處理方式不同方程數(shù)量也不同。我見過有人直接把單極代碼改個(gè)系數(shù)套到雙極上結(jié)果雅可比維度都配不平跑起來直接報(bào)維度不匹配。所以結(jié)論很明確交直流潮流不是在交流潮流上加個(gè)直流負(fù)荷而是要重新建立一套包含兩類變量、兩類方程的聯(lián)合模型。這一點(diǎn)想清楚了后面的工作才有方向。1.2 交替迭代和統(tǒng)一迭代的取舍邏輯工程上算交直流潮流有兩條路。第一條是交替迭代也叫順序求解先把直流側(cè)等效成交流節(jié)點(diǎn)的注入功率解一遍交流潮流再把交流母線電壓代回?fù)Q流器方程解一遍直流側(cè)兩個(gè)結(jié)果來回傳遞直到前后兩次的換流站功率差小于閾值。第二條是統(tǒng)一迭代把所有未知量放進(jìn)一個(gè)向量所有方程放進(jìn)一個(gè)殘差向量用牛頓法一次性求解。交替迭代的優(yōu)點(diǎn)是實(shí)現(xiàn)簡(jiǎn)單、代碼復(fù)用率高交流部分可以直接拿現(xiàn)成的牛頓法程序。但它的收斂性依賴一個(gè)隱含條件交流側(cè)和直流側(cè)之間的耦合不能太強(qiáng)。直流功率越大、換流站電氣距離越近、控制方式越接近切換邊界這個(gè)耦合就越強(qiáng)交替迭代的收斂速度會(huì)急劇下降甚至出現(xiàn)振蕩。我實(shí)測(cè)過一個(gè)直流功率 800 MW 的兩端系統(tǒng)交替迭代要 20 多次才勉強(qiáng)收斂把功率提到 950 MW 就直接發(fā)散了。統(tǒng)一迭代的代價(jià)是要自己推導(dǎo)并組裝擴(kuò)展雅可比矩陣工作量大概多出兩三天。換來的是牛頓法本身的二次收斂特性——只要初值不離譜殘差平方級(jí)下降。而且因?yàn)轳詈享?xiàng)顯式出現(xiàn)在雅可比里迭代過程對(duì)控制方式切換、直流功率大范圍變化的適應(yīng)性明顯更好。我現(xiàn)在這套代碼直流功率從 0 掃到 1.05 倍額定值除了極端工況需要調(diào)整初值基本都能收斂。選型上我的建議是如果只是做個(gè)粗略估算、直流容量占比很小交替迭代夠用但凡要做參數(shù)掃描、要做 N-1 分析、要研究控制方式切換直接上統(tǒng)一迭代省得后期反復(fù)返工。1.3 整套程序的結(jié)構(gòu)劃分代碼我按四層來組織這個(gè)分層方式用下來最省心。第一層是數(shù)據(jù)層把交流網(wǎng)絡(luò)參數(shù)支路、節(jié)點(diǎn)類型、負(fù)荷、直流網(wǎng)絡(luò)參數(shù)換流站位置、換流變變比、換相電抗、控制目標(biāo)值、直流電阻分開存成結(jié)構(gòu)體。第二層是模型層負(fù)責(zé)根據(jù)當(dāng)前狀態(tài)變量算出所有殘差和雅可比塊的元素。第三層是求解層做稀疏線性方程求解和狀態(tài)量更新。第四層是輸出層把結(jié)果整理成表格、算換流站的功率因數(shù)、校驗(yàn)功率平衡。這么分層的好處是調(diào)試時(shí)能單獨(dú)拎出某一層來驗(yàn)證。比如懷疑雅可比算錯(cuò)了就把模型層的數(shù)值雅可比和解析雅可比做一次差分對(duì)比誤差應(yīng)該在 1e-6 量級(jí)。這種自檢手段在后期定位問題時(shí)非常關(guān)鍵比盯著不收斂的殘差曲線瞎猜有效得多。整套代碼的核心文件我控制在五個(gè)左右主腳本、數(shù)據(jù)文件、殘差函數(shù)、雅可比函數(shù)、結(jié)果輸出函數(shù)。文件太多反而增加跳轉(zhuǎn)成本。2. 把物理系統(tǒng)拆成可求解的方程2.1 交流側(cè)節(jié)點(diǎn)功率平衡方程怎么寫交流側(cè)跟純交流潮流完全一樣每個(gè)節(jié)點(diǎn)兩組方程有功不平衡和無(wú)功不平衡。對(duì) PQ 節(jié)點(diǎn)兩個(gè)都寫對(duì) PV 節(jié)點(diǎn)只寫有功、無(wú)功由電壓幅值給定的條件替代平衡節(jié)點(diǎn)兩個(gè)都不寫。這部分我不贅述細(xì)節(jié)只強(qiáng)調(diào)一點(diǎn)換流節(jié)點(diǎn)必須按 PQ 節(jié)點(diǎn)處理因?yàn)樗挠泄o(wú)功都隨直流狀態(tài)變化不能指定電壓幅值。節(jié)點(diǎn)注入功率寫成標(biāo)準(zhǔn)形式P_i V_i * Σ_j V_j * (G_ij*cosθ_ij B_ij*sinθ_ij) Q_i V_i * Σ_j V_j * (G_ij*sinθ_ij - B_ij*cosθ_ij)換流節(jié)點(diǎn)上這個(gè) P_i 和 Q_i 里要額外加上直流注入項(xiàng)。注意符號(hào)約定我統(tǒng)一取交流系統(tǒng)向換流站注入為正那么整流站功率從交流流向直流的注入是負(fù)的逆變站是正的。符號(hào)約定如果一開始沒定死后面推導(dǎo)偏導(dǎo)時(shí)會(huì)非常痛苦很容易出現(xiàn)一整塊雅可比符號(hào)全反的情況。2.2 換流器穩(wěn)態(tài)方程與等值換相電阻這是整個(gè)模型的核心也是新手最容易栽跟頭的地方。理想換流器的穩(wěn)態(tài)方程忽略換相重疊時(shí)是U_d (3√2/π) * n_T * U_s * cosα考慮換相重疊后直流電壓被吃掉一塊用等值換相電阻表示U_d (3√2/π) * n_T * U_s * cosα - (3/π) * X_c * I_d這里 n_T 是換流變壓器變比折算到交流側(cè)電壓基準(zhǔn)X_c 是折算到同一基準(zhǔn)下的換相電抗I_d 是直流電流。把 (3/π) * X_c 記作 R_c就是所謂的等值換相電阻它不代表真實(shí)的有功損耗只是描述換相過程對(duì)直流電壓的壓降效應(yīng)。這個(gè)區(qū)分很重要——算功率平衡的時(shí)候R_c 上消耗的功率不是熱損耗而是換流器無(wú)功消耗的一部分體現(xiàn)。逆變側(cè)方程形式上要改符號(hào)因?yàn)橹绷麟娏鞯姆较蛳鄬?duì)逆變器是流出的U_d (3√2/π) * n_T * U_s * cosγ R_c * I_dγ 是熄弧角跟觸發(fā)角的關(guān)系是 γ π - α - μμ 是換相重疊角。實(shí)際工程里逆變側(cè)一般直接拿 γ 當(dāng)控制變量比用 α 更直觀因?yàn)?γ 直接對(duì)應(yīng)換相失敗的風(fēng)險(xiǎn)。寫代碼時(shí)我把常數(shù) K (3√2/π) * n_T 單獨(dú)算出來存著每次迭代只更新 cosα 和 I_d能省一點(diǎn)重復(fù)計(jì)算。別小看這點(diǎn)優(yōu)化幾百次迭代、上千個(gè)節(jié)點(diǎn)的時(shí)候累加起來的開銷是可觀的。2.3 直流網(wǎng)絡(luò)方程與控制方程的配對(duì)直流網(wǎng)絡(luò)本身很簡(jiǎn)單。兩端系統(tǒng)就是一條支路直流電阻 R_d兩端電壓 U_dr 和 U_di直流電流 I_dU_dr - U_di - R_d * I_d 0多端系統(tǒng)就寫成節(jié)點(diǎn)電流平衡的形式跟交流網(wǎng)絡(luò)的節(jié)點(diǎn)導(dǎo)納矩陣完全一個(gè)思路只是直流側(cè)只有電阻沒有電抗導(dǎo)納矩陣是實(shí)數(shù)對(duì)稱陣。我建議直流網(wǎng)絡(luò)也建成矩陣形式這樣擴(kuò)展到三端、四端的時(shí)候不用改代碼結(jié)構(gòu)只改數(shù)據(jù)。真正需要?jiǎng)幽X子的是控制方程。直流系統(tǒng)必須有兩個(gè)獨(dú)立的控制量才能確定工作點(diǎn)缺一個(gè)方程就欠定多一個(gè)就超定。常見的組合有這么幾種控制方式組合整流側(cè)方程逆變側(cè)方程適用場(chǎng)景定電流 定熄弧角I_d - I_set 0γ - γ_set 0傳統(tǒng)直流輸電最常見定功率 定直流電壓U_dr*I_d - P_set 0U_di - U_set 0背靠背、柔性直流常見定電流 定直流電壓I_d - I_set 0U_di - U_set 0一端定電流、一端穩(wěn)壓定觸發(fā)角 定熄弧角α - α_set 0γ - γ_set 0特殊研究工況這里有個(gè)關(guān)鍵細(xì)節(jié)我要專門說未知量的選擇直接決定方程形式。如果你把 I_d 當(dāng)成一個(gè)全局變量單極單回線只有一條直流電流那么定電流控制就是把 I_d 釘在給定值上方程直接是 I_d - I_set 0非常干凈。但如果兩端各有獨(dú)立的直流電流變量定電流方程就要分別寫變量數(shù)也多出來容易配不平。我的做法是每條直流支路只有一個(gè)電流變量多端系統(tǒng)里用節(jié)點(diǎn)電流平衡代替支路電流這樣變量數(shù)最少雅可比也最緊湊。2.4 未知量向量和方程向量的完整配平以雙端單極系統(tǒng)、n 個(gè)交流節(jié)點(diǎn)為例其中一個(gè)平衡節(jié)點(diǎn)、npv 個(gè) PV 節(jié)點(diǎn)未知量包括所有非平衡節(jié)點(diǎn)的電壓相角n-1 個(gè)所有 PQ 節(jié)點(diǎn)的電壓幅值n - 1 - npv 個(gè)直流側(cè)變量U_dr、U_di、I_d、α或整流側(cè)控制角、γ 共 5 個(gè)方程包括非平衡節(jié)點(diǎn)的有功不平衡方程n-1 個(gè)PQ 節(jié)點(diǎn)的無(wú)功不平衡方程n - 1 - npv 個(gè)直流回路方程 1 個(gè)整流側(cè)換流器方程 1 個(gè)逆變側(cè)換流器方程 1 個(gè)整流側(cè)控制方程 1 個(gè)逆變側(cè)控制方程 1 個(gè)數(shù)量剛好對(duì)上。我第一版代碼就是在這上面栽了把逆變側(cè)的換流器方程漏寫了維度差 1Matlab 報(bào)的是矩陣接近奇異排查了大半天才發(fā)現(xiàn)是方程少了一個(gè)。建議在寫代碼前先在紙上把變量清單和方程清單各列一遍一一對(duì)應(yīng)打勾。這一步花十分鐘能省后面十個(gè)小時(shí)。3. 統(tǒng)一迭代法雅可比矩陣的拼接方法3.1 擴(kuò)展雅可比的分塊結(jié)構(gòu)統(tǒng)一迭代的整體迭代格式還是牛頓法那一套J * Δx -F(x) x_new x Δx其中 F 是所有方程的殘差向量J 是 F 對(duì) x 的偏導(dǎo)矩陣。按交流變量和直流變量分塊J 長(zhǎng)這樣J [ J_AA J_AD ] [ J_DA J_DD ]J_AA 就是傳統(tǒng)交流潮流的雅可比這部分代碼可以直接復(fù)用。J_AD 是交流方程對(duì)直流變量的偏導(dǎo)J_DA 是直流方程對(duì)交流變量的偏導(dǎo)J_DD 是直流方程對(duì)直流變量的偏導(dǎo)。這兩個(gè)耦合塊是統(tǒng)一迭代的靈魂也是交替迭代法里根本沒有的東西。很多教程講到這里就停了只說把耦合項(xiàng)加上但不告訴你具體加什么。下面我把每一塊都攤開。3.2 交流側(cè)對(duì)直流變量的耦合項(xiàng)如何推導(dǎo)換流節(jié)點(diǎn)上直流注入交流的有功是P_dc U_d * I_d注意這是有功的絕對(duì)值符號(hào)靠前面約定的方向來定。對(duì)整流側(cè)注入是 -U_dI_d對(duì)逆變側(cè)是 U_dI_d。所以?P_dc/?U_d ±I_d ?P_dc/?I_d ±U_d ?P_dc/?α 0有功部分很干凈控制角不出現(xiàn)在有功表達(dá)式里。無(wú)功部分麻煩一些因?yàn)橐^功率因數(shù)這一關(guān)。換流器的功率因數(shù)有個(gè)漂亮的關(guān)系式可以從換流器方程反推出來cosφ U_d / (K * U_s)推導(dǎo)過程很簡(jiǎn)單cosα (U_d R_cI_d)/(KU_s)cos(αμ) (U_d - R_cI_d)/(KU_s)兩者相加除以二就是 cosφ。這個(gè)式子說明什么說明換流器的功率因數(shù)完全由直流電壓和交流母線電壓的比值決定跟觸發(fā)角沒有直接關(guān)系。那么無(wú)功功率Q_dc P_dc * tanφ U_d * I_d * sqrt(1 - (U_d/(K*U_s))2) / (U_d/(K*U_s)) I_d * sqrt((K*U_s)2 - U_d2)這個(gè)形式很好用偏導(dǎo)直接能算?Q_dc/?U_d -I_d * U_d / sqrt((K*U_s)2 - U_d2) ?Q_dc/?I_d sqrt((K*U_s)2 - U_d2) ?Q_dc/?U_s I_d * K2 * U_s / sqrt((K*U_s)2 - U_d2)最后一項(xiàng) ?Q_dc/?U_s 是容易被漏掉的關(guān)鍵項(xiàng)。它落的位置是 J_AA 塊里換流節(jié)點(diǎn)無(wú)功方程對(duì)自身電壓幅值的偏導(dǎo)。你如果只按傳統(tǒng)公式寫 J_AA這一項(xiàng)就丟了結(jié)果是雅可比不精確牛頓法退化成擬牛頓法收斂次數(shù)從 5 次漲到 15 次甚至不收斂。我第一版就是這個(gè)問題收斂慢得離譜后來做數(shù)值差分對(duì)比才發(fā)現(xiàn)漏了這塊。3.3 直流側(cè)對(duì)交流變量的耦合項(xiàng)怎么處理直流方程里顯式含有交流母線電壓的只有換流器方程。以整流側(cè)為例殘差寫成f_conv U_d - K * U_s * cosα R_c * I_d 0各偏導(dǎo)?f/?U_d 1 ?f/?U_s -K * cosα ?f/?α K * U_s * sinα ?f/?I_d R_c逆變側(cè)符號(hào)略有不同cosα 換成 cosγR_c 項(xiàng)符號(hào)反過來。直流回路方程對(duì)交流變量沒有偏導(dǎo)因?yàn)樗簧婕爸绷髁??？刂品匠桃彩嵌娏?、定熄弧角、定功率都不含交流電壓。所?J_DA 塊其實(shí)只有換流器方程那兩行是稀疏的非零其余全是零。這個(gè)稀疏性在組裝的時(shí)候要利用起來不然矩陣會(huì)白白變稠。3.4 標(biāo)幺化處理與稀疏矩陣組裝標(biāo)幺化這件事我必須單獨(dú)強(qiáng)調(diào)。交直流混合系統(tǒng)里交流側(cè)電壓是幾百千伏量級(jí)直流電流是千安量級(jí)功率是兆瓦量級(jí)如果不統(tǒng)一基準(zhǔn)雅可比矩陣的條件數(shù)會(huì)非常糟糕線性求解器給出的解精度很差。我的處理方式交流側(cè)沿用系統(tǒng)基準(zhǔn) S_base 和交流電壓基準(zhǔn) U_base阻抗基準(zhǔn) Z_base U_base2/S_base。直流側(cè)取直流電壓基準(zhǔn) U_dB 等于換流變交流側(cè)電壓基準(zhǔn)折算到直流側(cè)的值直流電流基準(zhǔn) I_dB S_base/U_dB直流電阻基準(zhǔn) R_dB U_dB/I_dB這正好等于 Z_base。這樣所有量都在 1 附近雅可比元素?cái)?shù)量級(jí)一致求解穩(wěn)定。代碼里稀疏組裝用三元組方式% 預(yù)分配按最大非零元數(shù)估算 max_nnz 4 * n 20; rows zeros(max_nnz, 1); cols zeros(max_nnz, 1); vals zeros(max_nnz, 1); idx 0; % 交流塊填充示例實(shí)際按行遍歷更省事 for i 1:n-1 % ... 逐行填充 J_AA end % 最后一次性組裝 J sparse(rows(1:idx), cols(1:idx), vals(1:idx), nvar, nvar);用 sparse 而不是 zeros 建滿矩陣10 個(gè)節(jié)點(diǎn)看不出區(qū)別100 個(gè)節(jié)點(diǎn)就能差出一個(gè)數(shù)量級(jí)的求解時(shí)間。而且 Matlab 的 \ 運(yùn)算符對(duì)稀疏矩陣會(huì)自動(dòng)選合適的求解器比手動(dòng)調(diào) lu 靠譜。4. Matlab 實(shí)現(xiàn)全流程與代碼骨架4.1 數(shù)據(jù)結(jié)構(gòu)怎么組織才不容易亂我用四個(gè)結(jié)構(gòu)體bus交流節(jié)點(diǎn)、branch交流支路、dc直流系統(tǒng)、opt求解選項(xiàng)。bus struct(num, [], ... % 節(jié)點(diǎn)編號(hào) type, [], ... % 1PQ, 2PV, 3平衡, 4換流站 V, [], ... % 電壓幅值初值標(biāo)幺 theta, [], ... % 相角初值弧度 Pg, [], Qg, [], Pl, [], Ql, []); dc struct(conv, [], ... % 換流站數(shù)組每站一個(gè)結(jié)構(gòu) line, []); % 直流線路參數(shù) dc.conv(k).bus_ac 5; % 接入的交流節(jié)點(diǎn)號(hào) dc.conv(k).nT 1.0; % 換流變變比 dc.conv(k).Xc 0.1; % 換相電抗標(biāo)幺 dc.conv(k).mode rect; % 整流還是逆變 dc.conv(k).ctrl I; % 控制變量類型 dc.conv(k).ctrl_val 1.0; % 控制目標(biāo)值用結(jié)構(gòu)體數(shù)組而不是元胞數(shù)組好處是可以直接用點(diǎn)號(hào)訪問代碼讀起來清楚。換流站數(shù)量變化時(shí)不用改任何邏輯程序里全部按 length(dc.conv) 循環(huán)。有個(gè)細(xì)節(jié)換流站的交流接入節(jié)點(diǎn)類型要標(biāo)成 4這樣在主循環(huán)里能快速定位哪些節(jié)點(diǎn)需要加直流注入項(xiàng)、哪些節(jié)點(diǎn)要跳過電壓幅值更新。用 type 字段區(qū)分比單獨(dú)維護(hù)一個(gè)索引數(shù)組更不容易出錯(cuò)。4.2 主迭代循環(huán)的完整骨架整個(gè)牛頓迭代循環(huán)結(jié)構(gòu)如下max_iter 20; tol 1e-8; x init_state(bus, dc); % 初始狀態(tài)向量 history zeros(max_iter, 1); for iter 1:max_iter [F, J] build_residual_jacobian(x, bus, branch, dc, Ybus, ...); err max(abs(F)); history(iter) err; if err tol fprintf(第 %d 次迭代收斂最大殘差 %.2e\n, iter, err); break; end if iter max_iter warning(達(dá)到最大迭代次數(shù)仍未收斂當(dāng)前殘差 %.2e, err); end dx -J \ F; % 稀疏左除自動(dòng)選求解器 x x dx; % 越限檢查和鉗位 x enforce_limits(x, dc); end幾個(gè)實(shí)操要點(diǎn)第一dx -J \ F這一行是稀疏求解的核心Matlab 會(huì)自動(dòng)判斷用 LU 還是 Cholesky不要手寫 inv(J)*F數(shù)值上差很多。第二enforce_limits這個(gè)函數(shù)判斷控制角的合理性。觸發(fā)角一般限制在 5 到 90 度熄弧角限制在 10 到 60 度。牛頓迭代偶爾會(huì)給出負(fù)的觸發(fā)角或超過 90 度的值這時(shí)候如果直接鉗位殘差會(huì)跳變迭代可能被帶偏我的做法是鉗位的同時(shí)把對(duì)應(yīng)的控制方程臨時(shí)改成角度不變相當(dāng)于人為把這個(gè)變量?jī)鼋Y(jié)一步讓其他變量先收斂。第三迭代日志一定要存。我把每次的殘差存進(jìn) history 數(shù)組收斂完成后畫個(gè)半對(duì)數(shù)圖能直接看出收斂類型——二次收斂是斜率越來越陡線性收斂是直線振蕩是鋸齒。這個(gè)圖對(duì)調(diào)試的價(jià)值極高。4.3 換流器殘差與雅可比的子函數(shù)實(shí)現(xiàn)把換流器的計(jì)算單獨(dú)封裝成函數(shù)輸入是直流狀態(tài)量、交流電壓輸出是殘差和對(duì)應(yīng)的雅可比元素function [f, Jrow] conv_equation(Ud, Id, Us, alpha, K, Rc, mode) if strcmp(mode, rect) f Ud - K * Us * cos(alpha) Rc * Id; Jrow struct(dUd, 1, ... dUs, -K * cos(alpha), ... dalpha, K * Us * sin(alpha), ... dId, Rc); else f Ud - K * Us * cos(alpha) - Rc * Id; Jrow struct(dUd, 1, ... dUs, -K * cos(alpha), ... dalpha, K * Us * sin(alpha), ... dId, -Rc); end end無(wú)功對(duì)電壓幅值的偏導(dǎo)那個(gè) sqrt((KU_s)2 - U_d2) 分母有開方當(dāng) U_d 接近 KU_s 時(shí)cosφ 接近 1換流器運(yùn)行在接近單位功率因數(shù)的極端狀態(tài)這個(gè)偏導(dǎo)會(huì)趨于零反過來 ?Q/?U_d 會(huì)趨于無(wú)窮。實(shí)際系統(tǒng)里 U_d/(K*U_s) 一般在 0.85 到 0.95 之間離奇異點(diǎn)還有距離但代碼里要加個(gè)保護(hù)ratio Ud / (K * Us); ratio min(max(ratio, 0.5), 0.999); % 防止開方為負(fù)或除零這個(gè)保護(hù)不是多余的。我遇到過一次逆變側(cè)電壓被迭代推到很低、ratio 超過 1sqrt 直接返回復(fù)數(shù)整個(gè)迭代崩掉。加上鉗位之后至少不會(huì)出復(fù)數(shù)雖然結(jié)果不物理但能讓你從殘差曲線看出是哪一步跑偏的。4.4 初值選取與收斂控制牛頓法對(duì)初值敏感交直流系統(tǒng)尤其如此。我的初值策略是分層給交流側(cè)全部給 1.0 幅值、0 相角這個(gè)標(biāo)準(zhǔn)做法基本夠用。直流側(cè)要稍微算一下直流電壓初值按額定值給直流電流初值由目標(biāo)功率反算 I_d P_set / U_d整流側(cè)觸發(fā)角給 15 度約 0.26 弧度逆變側(cè)熄弧角給 18 度。為什么給 15 到 18 度這個(gè)范圍因?yàn)檫@是工程上實(shí)際運(yùn)行點(diǎn)的典型值牛頓法從這附近出發(fā)基本一步就能進(jìn)收斂域。你要是給 0 度或者 80 度第一第二個(gè)迭代步會(huì)跑得很遠(yuǎn)有時(shí)候直接跑到非物理解那邊去了。收斂判據(jù)我用最大殘差而不是二范數(shù)。二范數(shù)在大系統(tǒng)里會(huì)被節(jié)點(diǎn)數(shù)量放大閾值不好定最大殘差直觀1e-8 就是每個(gè)方程都滿足到 8 位有效數(shù)字物理意義明確。還有一個(gè)技巧前兩次迭代可以用放寬的雅可比也就是用初值處的雅可比多迭代幾步省下重新組裝和分解的時(shí)間。但對(duì)交直流系統(tǒng)來說耦合項(xiàng)變化快這個(gè)技巧收益不大我最后沒采用老老實(shí)實(shí)每步重算。5. 實(shí)際調(diào)試中踩過的坑5.1 不收斂的幾種典型癥狀和定位路徑殘差曲線是最好的診斷工具。我總結(jié)了幾種形態(tài)和對(duì)應(yīng)的病因殘差表現(xiàn)可能原因排查手段首次迭代后殘差反而變大雅可比符號(hào)錯(cuò)誤或初值嚴(yán)重偏離數(shù)值差分校驗(yàn)雅可比殘差線性下降每次只降一個(gè)數(shù)量級(jí)耦合項(xiàng)漏寫退化成擬牛頓檢查 J_AD 和 J_DA 是否為零殘差在 1e-3 附近來回振蕩控制方式在切換邊界反復(fù)跳打印每步觸發(fā)角看是否觸發(fā)限值殘差突然變成 NaN 或復(fù)數(shù)開方項(xiàng)為負(fù)、除零、角度越界在各處加 assert 斷言前三次下降很快之后停滯雅可比接近奇異算 cond(J)檢查是否有零對(duì)角元數(shù)值差分校驗(yàn)雅可比是排查的第一步也是最有效的一步。做法很簡(jiǎn)單對(duì)每個(gè)狀態(tài)變量給一個(gè) 1e-6 的擾動(dòng)用前向差分算數(shù)值偏導(dǎo)跟解析偏導(dǎo)逐個(gè)比較相對(duì)誤差超過 1e-4 的就是有問題。我第二次遇到收斂慢的問題時(shí)就是這么定位的查出來是 ?Q/?U_s 那一項(xiàng)少了個(gè) K2 因子。5.2 收斂到非物理解的表現(xiàn)和規(guī)避牛頓法有個(gè)討厭的特性它只管把殘差推到零不管解的物理意義。交直流系統(tǒng)里很容易出現(xiàn)這種情況——迭代收斂了殘差 1e-10但打印出來發(fā)現(xiàn)整流側(cè)觸發(fā)角是負(fù)的 20 度或者逆變側(cè)熄弧角是 3 度這明顯是換相失敗工況不可能是穩(wěn)態(tài)解。原因通常有兩個(gè)。一是初值離物理解太遠(yuǎn)牛頓步跨過了物理可行域落到了數(shù)學(xué)解的另一支上。二是控制方式配錯(cuò)了比如整流側(cè)明明應(yīng)該是定電流你按定功率寫了方程系統(tǒng)會(huì)找到一個(gè)數(shù)學(xué)上滿足定功率但物理上不可行的點(diǎn)。規(guī)避辦法每次迭代后檢查所有控制角的物理范圍超出范圍就打印警告并記錄收斂后如果仍有越界就說明這次求解無(wú)效需要換初值重來。我在主程序里加了自動(dòng)重啟機(jī)制如果首次求解出現(xiàn)越界就把觸發(fā)角初值從 15 度調(diào)到 20 度重新算一次通常第二次就能落在合理區(qū)域。5.3 結(jié)果校驗(yàn)的三種手段算完之后必須校驗(yàn)不能收斂了就用。我固定做三件事第一全系統(tǒng)有功平衡校驗(yàn)。所有交流電源的有功加上直流注入的有功減去負(fù)荷有功和網(wǎng)絡(luò)損耗應(yīng)該接近零。誤差在 1e-6 量級(jí)說明有功方程沒問題。這里要注意等值換相電阻 R_c 上的功率不是損耗不能算進(jìn)網(wǎng)損這是新手常犯的錯(cuò)誤。第二換流站交直流功率對(duì)照。交流側(cè)算出來的注入有功跟直流側(cè) U_d * I_d 應(yīng)該相等。這兩條路徑算出來的值如果差得多說明換流器方程或者功率因數(shù)公式用錯(cuò)了。我一般要求相對(duì)誤差小于 1e-6。第三跟交替迭代的結(jié)果做交叉驗(yàn)證。同一組數(shù)據(jù)用兩套方法算結(jié)果應(yīng)該一致到小數(shù)點(diǎn)后四五位。這個(gè)交叉驗(yàn)證做過一次之后后面改代碼就有信心了。5.4 我的幾條實(shí)操心得換相電抗的單位是最容易出錯(cuò)的地方。數(shù)據(jù)手冊(cè)上給的往往是折算到換流變閥側(cè)的歐姆值你必須先折算到交流母線側(cè)再除以阻抗基準(zhǔn)變成標(biāo)幺值。少除一個(gè)變比平方結(jié)果能差出百分之幾而且不容易發(fā)現(xiàn)因?yàn)檎諛邮諗?。直流電阻的處理也一樣雙極系統(tǒng)的直流線路電阻是單極的一半兩極并聯(lián)金屬回線方式下回路電阻還要加上回線電阻。這些細(xì)節(jié)數(shù)據(jù)文件里寫錯(cuò)了結(jié)果就是錯(cuò)的程序本身不會(huì)報(bào)錯(cuò)。還有一點(diǎn)換流變壓器的變比定義有歧義。有的資料定義成閥側(cè)額定電壓/交流側(cè)額定電壓有的定義成交流側(cè)電壓/閥側(cè)電壓兩者互為倒數(shù)。用錯(cuò)了直流電壓直接差一個(gè)變比平方功率就差得更多。我的建議是在數(shù)據(jù)文件里明確注釋清楚定義并且用一個(gè)已知工況反算驗(yàn)證一遍。6. 這套代碼后續(xù)能怎么擴(kuò)展現(xiàn)在這套代碼基本框架是穩(wěn)的往幾個(gè)方向擴(kuò)展都不難。多端直流就是把直流網(wǎng)絡(luò)從單支路推廣到節(jié)點(diǎn)導(dǎo)納矩陣形式方程結(jié)構(gòu)不變只是變量和方程數(shù)量增加雅可比矩陣的稀疏性要重新組織一下。加入直流網(wǎng)絡(luò)損耗的精細(xì)模型也不難把線路電阻從常數(shù)改成隨溫度變化的函數(shù)或者加上線路的分布參數(shù)效應(yīng)。如果要做含新能源接入的混合系統(tǒng)難點(diǎn)不在潮流本身而在于新能源的出力不確定性——那就變成概率潮流或者場(chǎng)景分析了需要外層套蒙特卡洛或者點(diǎn)估計(jì)法。這時(shí)候統(tǒng)一迭代法的優(yōu)勢(shì)更明顯因?yàn)閱未吻蠼饪鞄浊Т螆?chǎng)景掃描才跑得動(dòng)。另外一個(gè)我覺得值得做的方向是把換流器的詳細(xì)模型引進(jìn)來考慮換相重疊角的動(dòng)態(tài)變化和閥的導(dǎo)通壓降。這樣算出來的穩(wěn)態(tài)點(diǎn)更接近實(shí)際尤其是低功率運(yùn)行的時(shí)候理想模型的誤差會(huì)比較大。代價(jià)是方程更復(fù)雜雅可比的解析推導(dǎo)工作量翻倍可能需要借助符號(hào)計(jì)算工具輔助推導(dǎo)。最后分享一個(gè)我一直在用的小習(xí)慣每次改動(dòng)代碼之后先用一個(gè) 3 節(jié)點(diǎn)的小系統(tǒng)跑一遍冒煙測(cè)試確認(rèn)能收斂、功率平衡滿足再去跑 30 節(jié)點(diǎn)以上的實(shí)際算例。小系統(tǒng)的迭代過程可以直接在命令行里 print 出來一眼就能看出哪一步不對(duì)大系統(tǒng)一旦出問題定位成本高十倍。這個(gè)習(xí)慣幫我省下的時(shí)間比我優(yōu)化代碼省下的時(shí)間多得多。