簡(jiǎn)介面向電力系統(tǒng)分析與潮流計(jì)算學(xué)習(xí)者的C工具包聚焦穩(wěn)態(tài)網(wǎng)絡(luò)電壓、電流與功率分布的求解適合課程設(shè)計(jì)、科研復(fù)現(xiàn)及工程入門(mén)者參考。壓縮包共46個(gè)文件、約2.23MB以18個(gè)txt算例與結(jié)果數(shù)據(jù)、3個(gè)xlsx表格、2個(gè)cpp源文件和2個(gè)h頭文件為主另含Visual Studio工程及編譯中間文件可直接編譯運(yùn)行并對(duì)照輸入輸出。程序采用牛頓法迭代求解非線性潮流方程并借助tinny2稀疏矩陣技術(shù)降低大規(guī)模節(jié)點(diǎn)網(wǎng)絡(luò)的存儲(chǔ)與計(jì)算開(kāi)銷(xiāo)提供IEEE14至300等標(biāo)準(zhǔn)系統(tǒng)算例便于驗(yàn)證收斂性與效率。隨附報(bào)告涵蓋設(shè)計(jì)思路、算法描述與結(jié)果分析已有840人學(xué)習(xí)。讀者可據(jù)此理解潮流計(jì)算從方程線性化、雅可比矩陣組裝到稀疏求解的完整實(shí)現(xiàn)路徑也能借助多規(guī)模算例開(kāi)展性能比較與二次開(kāi)發(fā)。1. 一個(gè)“收斂了但結(jié)果不對(duì)”的 30 節(jié)點(diǎn)算例去年幫一位做配網(wǎng)規(guī)劃的朋友看一個(gè) 30 節(jié)點(diǎn)算例程序收斂得干干凈凈可現(xiàn)場(chǎng)拿到的母線電壓與程序結(jié)果差了 0.08 p.u.查了一整天才發(fā)現(xiàn)變壓器變比填在了 to 端代碼卻按 from 端讀。潮流計(jì)算powerflow最麻煩的情況往往是“算出來(lái)了但你不知道能不能信”并非“算不出來(lái)”。這套電力系統(tǒng)潮流計(jì)算通用程序針對(duì)的就是這一環(huán)節(jié)節(jié)點(diǎn)、支路、變壓器、發(fā)電機(jī)拆成四張標(biāo)準(zhǔn)表內(nèi)置牛頓-拉夫遜、快速解耦、高斯-賽德?tīng)柸N解法支持 PV/PQ/平衡節(jié)點(diǎn)混排、無(wú)功越限自動(dòng)切換、批量算例掃描輸出母線電壓、支路潮流、網(wǎng)損三張報(bào)表。想快速驗(yàn)證網(wǎng)架方案的可以直接用想把求解器拆開(kāi)換算法、接自己數(shù)據(jù)格式的也能改。2. 架構(gòu)與數(shù)據(jù)模型四張表?yè)纹鹫麄€(gè)潮流計(jì)算這套程序最核心的設(shè)計(jì)取舍是把“數(shù)據(jù)”和“算法”徹底切開(kāi)。求解器不認(rèn)得 bus 表長(zhǎng)什么樣它只吃三樣?xùn)|西節(jié)點(diǎn)導(dǎo)納矩陣Ybus、注入功率向量Sbus、以及平衡節(jié)點(diǎn)/PV 節(jié)點(diǎn)/PQ 節(jié)點(diǎn)的索引數(shù)組。這條邊界劃清楚了換算法就是換一個(gè)函數(shù)換數(shù)據(jù)格式就是換一個(gè) reader兩邊互不干擾邊界沒(méi)劃清楚改一個(gè) CSV 列名就能讓整個(gè)牛頓迭代崩掉這類(lèi)返工我見(jiàn)過(guò)不止一次。2.1 目錄結(jié)構(gòu)io、model、solver、report 四層拿到源碼包先別急著python run_pf.py把目錄看一遍比跑十次都值。這套程序的目錄基本長(zhǎng)這樣powerflow/ ├── run_pf.py # 命令行入口參數(shù)解析 調(diào)度 ├── config/ │ └── default.yaml # 默認(rèn)容差、迭代上限、輸出開(kāi)關(guān) ├── pf/ │ ├── model.py # Bus/Branch/Gen/Load 四張表 字段校驗(yàn) │ ├── ybus.py # 節(jié)點(diǎn)導(dǎo)納矩陣組裝含變壓器、對(duì)地并聯(lián) │ ├── caseio.py # JSON/CSV 讀寫(xiě)標(biāo)幺化統(tǒng)一在這一層做 │ ├── qlimit.py # 無(wú)功越限檢測(cè)與 PV→PQ 切換 │ ├── solver/ │ │ ├── registry.py # 求解器注冊(cè)表二次開(kāi)發(fā)的掛載點(diǎn) │ │ ├── newton.py # 牛頓-拉夫遜 │ │ ├── fdpf.py # 快速解耦 │ │ └── gauss.py # 高斯-賽德?tīng)?│ └── report.py # 母線/支路/網(wǎng)損三張報(bào)表 ├── cases/ │ ├── case14.json │ └── case30.json ├── ref/ │ └── case14_voltage_ref.csv └── tests/ └── test_newton_case14.py四層的職責(zé)很明確。caseio負(fù)責(zé)把外部文件JSON、CSV、甚至別的工具導(dǎo)出的文本讀成內(nèi)存里的四張表并且一次性完成標(biāo)幺化——電壓除以 baseKV、功率除以 baseMVA、阻抗除以 Zbase。很多新手版本把標(biāo)幺化散落在各個(gè)模塊里最后一定會(huì)出現(xiàn)“線路阻抗是標(biāo)幺值、負(fù)荷功率是有名值”的混合狀態(tài)算出來(lái)的電壓能到 0.3 p.u.。model只做校驗(yàn)不做計(jì)算檢查有沒(méi)有孤立的節(jié)點(diǎn)編號(hào)、支路兩端是否都在 bus 表里、PV 節(jié)點(diǎn)的 Vm 是否越界、平衡節(jié)點(diǎn)是不是只有一個(gè)。solver是純算法層。report把求解器吐出來(lái)的復(fù)數(shù)電壓向量還原成人能看的 kV、MW、MVar。提示先看config/default.yaml里面決定了容差和迭代上限默認(rèn)容差如果被改成1e-4牛頓法通常 3 次就“收斂”但結(jié)果只能信到小數(shù)點(diǎn)后兩位。2.2 四張核心表bus、branch、gen、load 的字段設(shè)計(jì)數(shù)據(jù)模型是一套潮流程序的骨架。這套程序用的是電力系統(tǒng)里最通用的“四表結(jié)構(gòu)”字段命名也和主流工具對(duì)齊方便你把別處的算例直接導(dǎo)進(jìn)來(lái)。bus 表字段類(lèi)型含義說(shuō)明bus_idint節(jié)點(diǎn)編號(hào)從 1 開(kāi)始內(nèi)部索引統(tǒng)一減 1typeint節(jié)點(diǎn)類(lèi)型1PQ2PV3平衡節(jié)點(diǎn)Vmfloat電壓幅值p.u.初值一般給 1.0Vafloat電壓相角度求解器內(nèi)部轉(zhuǎn)弧度Pd/Qdfloat負(fù)荷p.u.Pg/Qgfloat發(fā)電機(jī)出力p.u.PQ 節(jié)點(diǎn)通常為 0Qmax/Qminfloat無(wú)功上下限越限切換用baseKVfloat基準(zhǔn)電壓報(bào)表還原用branch 表字段含義備注f, t首端、末端節(jié)點(diǎn)號(hào)變壓器變比默認(rèn)掛在 f 端r, x, b電阻、電抗、對(duì)地電納p.u.b 是全長(zhǎng)充電電納rateA長(zhǎng)期允許容量MVA報(bào)表里標(biāo)越限tap變比0 視為 1.0shift移相角度status投運(yùn)狀態(tài)0 表示退出組裝 Ybus 時(shí)跳過(guò)gen表和load表是 bus 表的補(bǔ)充一個(gè)節(jié)點(diǎn)可以掛多臺(tái)機(jī)組、多個(gè)負(fù)荷caseio在加載階段把它們按bus_id聚合回Pg/Qg/Pd/Qd。這種“明細(xì)表 聚合”的寫(xiě)法比直接在 bus 表里堆字段靈活尤其適合后續(xù)要接時(shí)序數(shù)據(jù)的場(chǎng)景。常見(jiàn)做法是在model.py里定義一個(gè)輕量的Case數(shù)據(jù)類(lèi)把四張表裝進(jìn)去順便緩存一次索引from dataclasses import dataclass import numpy as np dataclass class Case: bus: np.ndarray # (nb, ncol) 浮點(diǎn)/整型混合用列號(hào)訪問(wèn) branch: np.ndarray # (nl, ncol) gen: np.ndarray load: np.ndarray baseMVA: float 100.0 def __post_init__(self): # 節(jié)點(diǎn)類(lèi)型索引求解器只認(rèn)這三個(gè)數(shù)組不認(rèn) bus 表 types self.bus[:, 1].astype(int) self.ref np.where(types 3)[0] # 平衡節(jié)點(diǎn)必須恰好 1 個(gè) self.pv np.where(types 2)[0] self.pq np.where(types 1)[0] if len(self.ref) ! 1: raise ValueError(f平衡節(jié)點(diǎn)數(shù)量應(yīng)為 1實(shí)際 {len(self.ref)})__post_init__里做索引緩存是這套代碼里最劃算的一處優(yōu)化牛頓迭代每輪都要用pvpq去切雅可比矩陣如果每輪都重新where一遍30 節(jié)點(diǎn)看不出來(lái)118 節(jié)點(diǎn)、上千次批量掃描時(shí)的開(kāi)銷(xiāo)就很明顯了。ref必須恰好一個(gè)這條校驗(yàn)也放在這里因?yàn)槠胶夤?jié)點(diǎn)缺失是“雅可比矩陣奇異”最常見(jiàn)的誘因早報(bào)錯(cuò)比讓np.linalg.solve拋LinAlgError友好得多。2.3 Ybus 組裝稀疏矩陣、變壓器變比與對(duì)地并聯(lián)節(jié)點(diǎn)導(dǎo)納矩陣是整個(gè)程序的地基寫(xiě)錯(cuò)了后面全是白算。標(biāo)準(zhǔn)支路含變壓器的四項(xiàng)導(dǎo)納貢獻(xiàn)是這樣的import numpy as np from scipy.sparse import lil_matrix def build_ybus(case, include_shuntTrue): nb case.bus.shape[0] Y lil_matrix((nb, nb), dtypecomplex) for k in range(case.branch.shape[0]): br case.branch[k] if int(br[10]) 0: # status0停運(yùn)支路直接跳過(guò) continue f, t int(br[0]) - 1, int(br[1]) - 1 r, x, b br[2], br[3], br[4] tap, shift br[8], br[9] if tap 0: tap 1.0 ys 1.0 / complex(r, x) # 串聯(lián)導(dǎo)納 tap_c tap * np.exp(1j * np.deg2rad(shift)) # 復(fù)數(shù)變比 # 變壓器按 f 端建模f 側(cè)除以 |tap|^2t 側(cè)不變 Y[f, f] (ys 1j * b / 2) / (tap_c * np.conj(tap_c)) Y[t, t] (ys 1j * b / 2) Y[f, t] -ys / np.conj(tap_c) Y[t, f] -ys / tap_c if include_shunt: for i in range(nb): Y[i, i] complex(case.bus[i, 9], case.bus[i, 10]) # Gs jBs return Y.tocsr()這段代碼有三個(gè)地方值得反復(fù)看。第一是對(duì)地電納要除以 2b在 branch 表里存的是全長(zhǎng)充電電納等分到兩端所以(ys 1j*b/2)如果算例給的是半長(zhǎng)值這里就會(huì)多出一倍表現(xiàn)為輕載時(shí)電壓偏高、網(wǎng)損偏小屬于典型的“結(jié)果看起來(lái)挺合理”的錯(cuò)誤。第二是變比的共軛位置。Y[f,t] -ys / conj(tap_c)、Y[t,f] -ys / tap_c這兩個(gè)共軛一旦寫(xiě)反變壓器兩側(cè)的有功方向就整體反了。我見(jiàn)過(guò)有人在調(diào)試時(shí)把這兩行對(duì)調(diào)“試試看”結(jié)果潮流方向?qū)α?、網(wǎng)損不對(duì)因?yàn)閷?duì)調(diào)只是把變比挪到了 t 側(cè)而表里存的是 f 側(cè)。第三是稀疏結(jié)構(gòu)的選擇。lil_matrix適合逐元素累加組裝完立刻tocsr()轉(zhuǎn)成行壓縮格式后面Ybus V的矩陣向量乘法就快了。用稠密numpy數(shù)組在 30 節(jié)點(diǎn)沒(méi)問(wèn)題上萬(wàn)節(jié)點(diǎn)的配網(wǎng)會(huì)直接把內(nèi)存吃滿。注意status字段的位置第 10 列在不同算例文件里可能不一致導(dǎo)入前一定對(duì)著表頭核一遍。索引錯(cuò)一列程序不會(huì)報(bào)錯(cuò)只會(huì)安靜地算錯(cuò)。3. 核心算法拆解一次牛頓迭代到底發(fā)生了什么牛頓-拉夫遜法之所以是潮流計(jì)算的主力是因?yàn)樗谡＿\(yùn)行方式下二次收斂——誤差從 1e-2 掉到 1e-4再到 1e-8通常 4 到 5 次迭代就夠。但“二次收斂”是有前提的初值不能太離譜雅可比矩陣不能奇異不平衡量的符號(hào)不能搞反。這一章就沿著newton_pf這條鏈路把一次迭代拆成三段看。3.1 不平衡量從 Sbus 到 F 向量迭代開(kāi)始前要先拼出指定注入功率Sbus。規(guī)則很樸素PQ 節(jié)點(diǎn)取-(Pd jQd)PV 節(jié)點(diǎn)取(Pg - Pd) j(-Qd)平衡節(jié)點(diǎn)的有功無(wú)功都當(dāng)成未知量、不參與方程。def make_sbus(case): nb case.bus.shape[0] Sbus np.zeros(nb, dtypecomplex) for i in range(nb): Pd, Qd case.bus[i, 4], case.bus[i, 5] Pg, Qg case.bus[i, 6], case.bus[i, 7] if i in case.pq: Sbus[i] complex(-Pd, -Qd) # PQ注入 -負(fù)荷 elif i in case.pv: Sbus[i] complex(Pg - Pd, -Qd) # PVQ 是未知量置 0 else: Sbus[i] complex(Pg - Pd, Qg - Qd) # 平衡節(jié)點(diǎn)實(shí)際不參與方程 return SbusPV 節(jié)點(diǎn)的無(wú)功為什么要置 0因?yàn)樗臒o(wú)功出力是待求量方程里只保留有功平衡無(wú)功方程直接不列。等迭代收斂后再用Q imag(V * conj(Ybus V))反算出這臺(tái)機(jī)組的實(shí)際無(wú)功去跟Qmax/Qmin比。這一步在qlimit.py里做越限就把該節(jié)點(diǎn)從 PV 改成 PQ無(wú)功固定在上限或下限重新迭代。有了Sbus不平衡量就是一行S V * np.conj(Ybus V) mis np.r_[(Sbus.real - S.real)[pvpq], (Sbus.imag - S.imag)[pq]]注意這里是指定值減計(jì)算值。符號(hào)反過(guò)來(lái)迭代的方向就反了表現(xiàn)為電壓越算越離譜、兩三次之后直接發(fā)散。判斷依據(jù)很簡(jiǎn)單如果第一次迭代的不平衡量是正數(shù)、第二次變成更大的正數(shù)八成是符號(hào)搞反了。pvpq是np.r_[case.pv, case.pq]順序無(wú)所謂但整個(gè)程序里必須保持一致——雅可比矩陣的行列、狀態(tài)量的回寫(xiě)位置都依賴(lài)這個(gè)順序。中途有人把pvpq改成np.sort(...)結(jié)果收斂了但電壓對(duì)不上就是順序不一致導(dǎo)致的錯(cuò)位。3.2 雅可比矩陣的四個(gè)分塊雅可比矩陣是[dP/dVa, dP/dVm; dQ/dVa, dQ/dVm]的四個(gè)分塊規(guī)模是(npvpqnpq) × (npvpqnpq)。手推公式容易出錯(cuò)實(shí)踐中更穩(wěn)的寫(xiě)法是先算復(fù)數(shù)域的dS/dVa和dS/dVm再取實(shí)部虛部。def dSbus_dV(Ybus, V): 返回 dS/dVa 與 dS/dVm均為 nb×nb 復(fù)矩陣 Y Ybus.toarray() if hasattr(Ybus, toarray) else np.asarray(Ybus) Ibus Y V diagV, diagI np.diag(V), np.diag(Ibus) diagVn np.diag(V / np.abs(V)) # 單位相量 dS_dVm diagV np.conj(Y diagVn) np.conj(diagI) diagVn dS_dVa 1j * diagV np.conj(diagI - Y diagV) return dS_dVa, dS_dVm這兩行公式看著唬人來(lái)源其實(shí)很清楚S diag(V) · conj(Y·V)對(duì)Va求導(dǎo)時(shí)V的幅值不變、只轉(zhuǎn)角度所以多了個(gè)1j對(duì)Vm求導(dǎo)時(shí)V Vm · e^{jθ}幅值變化等價(jià)于乘單位相量。diagI - Y diagV這一項(xiàng)是鏈?zhǔn)椒▌t帶出來(lái)的漏掉它雅可比矩陣就退化成一個(gè)近似版本迭代次數(shù)從 4 次變成十幾次節(jié)點(diǎn)一多直接不收斂。分塊組裝dS_dVa, dS_dVm dSbus_dV(Ybus, V) J11 dS_dVa[np.ix_(pvpq, pvpq)].real J12 dS_dVm[np.ix_(pvpq, pq)].real J21 dS_dVa[np.ix_(pq, pvpq)].imag J22 dS_dVm[np.ix_(pq, pq)].imag J np.block([[J11, J12], [J21, J22]])np.ix_是這里的常用技巧行列同時(shí)切片拿到的是子矩陣如果寫(xiě)成dS_dVa[pvpq][:, pvpq]兩次索引會(huì)先復(fù)制再切節(jié)點(diǎn)上千時(shí)性能差一個(gè)量級(jí)而且中間結(jié)果占內(nèi)存。J12的列是pq只有 PQ 節(jié)點(diǎn)有電壓幅值方程J21的行是pq只有 PQ 節(jié)點(diǎn)有無(wú)功方程這兩個(gè)索引最容易被統(tǒng)一寫(xiě)成pvpq一寫(xiě)就錯(cuò)。提示調(diào)試時(shí)可以檢查J的對(duì)角元素。正常情況下J11的對(duì)角元是負(fù)數(shù)有功對(duì)相角J22的對(duì)角元也是負(fù)數(shù)無(wú)功對(duì)幅值如果出現(xiàn)大面積正值基本可以斷定Ybus的符號(hào)寫(xiě)反了。3.3 修正方程求解與狀態(tài)量回寫(xiě)把上面幾段拼起來(lái)就是完整的迭代循環(huán)def newton_pf(Ybus, Sbus, V0, case, tol1e-8, max_iter20): V V0.astype(complex).copy() pvpq np.r_[case.pv, case.pq] npvpq, npq len(pvpq), len(case.pq) history [] for it in range(max_iter): S V * np.conj(Ybus V) mis np.r_[(Sbus.real - S.real)[pvpq], (Sbus.imag - S.imag)[case.pq]] err np.max(np.abs(mis)) history.append(err) if err tol: return V, it, history # 收斂返回迭代次數(shù) dS_dVa, dS_dVm dSbus_dV(Ybus, V) J np.block([ [dS_dVa[np.ix_(pvpq, pvpq)].real, dS_dVm[np.ix_(pvpq, case.pq)].real], [dS_dVa[np.ix_(case.pq, pvpq)].imag, dS_dVm[np.ix_(case.pq, case.pq)].imag], ]) dx np.linalg.solve(J, mis) # F Sspec - ScalcF -J故 dx J?1·mis Va, Vm np.angle(V), np.abs(V) Va[pvpq] dx[:npvpq] Vm[case.pq] dx[npvpq:] V Vm * np.exp(1j * Va) raise RuntimeError(f牛頓法 {max_iter} 次未收斂末次不平衡量 {history[-1]:.3e})dx solve(J, mis)這個(gè)符號(hào)關(guān)系值得單獨(dú)說(shuō)。設(shè)F(x) Sspec - Scalc(x)我們要F0牛頓步是x_{k1} x_k - F/F而F -J因?yàn)镴是Scalc的導(dǎo)數(shù)所以dx J?1 · F直接加而不是減。很多實(shí)現(xiàn)里寫(xiě)成dx -solve(J, mis)然后不平衡量又定義成Scalc - Sspec兩次負(fù)號(hào)抵消看起來(lái)也能跑但只要有人改動(dòng)其中一個(gè)整套就崩?；貙?xiě)只更新Va[pvpq]和Vm[pq]平衡節(jié)點(diǎn)的幅值相角都固定PV 節(jié)點(diǎn)的幅值固定、只更新相角。寫(xiě)成Vm[:] ...是最常見(jiàn)的越界錯(cuò)誤PV 節(jié)點(diǎn)電壓被推著走表現(xiàn)為這臺(tái)機(jī)組的無(wú)功越算越大、最后報(bào)越限。max_iter默認(rèn) 20 已經(jīng)足夠?qū)捤烧Ｋ憷?5 次內(nèi)收斂。如果 20 次還不收斂別急著調(diào)大——先看history的衰減速度誤差不降反升是發(fā)散誤差在 1e-3 附近反復(fù)橫跳是震蕩兩種情況的原因完全不同前者多半是初值或數(shù)據(jù)問(wèn)題后者常見(jiàn)于重載或接近電壓崩潰點(diǎn)的算例。3.4 快速解耦與高斯-賽德?tīng)柺裁磮?chǎng)景該降級(jí)牛頓法精度高但每輪都要重新組裝并分解雅可比矩陣118 節(jié)點(diǎn)以下的算例感覺(jué)不到上千節(jié)點(diǎn)、批量掃描幾百個(gè)方式時(shí)就慢了。這套程序提供了兩種降級(jí)方案?？焖俳怦頟Q 解耦法利用高壓網(wǎng)里R X的特點(diǎn)把雅可比矩陣近似成兩個(gè)常數(shù)矩陣迭代中不再更新def fast_decoupled_pf(Ybus, Sbus, V0, case, tol1e-8, max_iter30): Y Ybus.toarray() if hasattr(Ybus, toarray) else Ybus V V0.astype(complex).copy() pvpq, pq np.r_[case.pv, case.pq], case.pq Bp -Y.imag[np.ix_(pvpq, pvpq)] # B有功-相角 Bpp -Y.imag[np.ix_(pq, pq)] # B無(wú)功-幅值 for it in range(max_iter): S V * np.conj(Y V) dP (Sbus.real - S.real)[pvpq] dQ (Sbus.imag - S.imag)[pq] if max(np.max(np.abs(dP)), np.max(np.abs(dQ))) tol: return V, it Va, Vm np.angle(V), np.abs(V) Va[pvpq] np.linalg.solve(Bp, dP / Vm[pvpq]) Vm[pq] np.linalg.solve(Bpp, dQ / Vm[pq]) V Vm * np.exp(1j * Va) raise RuntimeError(快速解耦未收斂)Bp、Bpp是常數(shù)矩陣可以在循環(huán)外做一次 LU 分解之后每輪只做前代回代這就是它在批量掃描里快的原因。代價(jià)是收斂變成線性迭代次數(shù)通常 8 到 15 次。使用時(shí)有一條邊界條件線路的 R/X 比值不能太大。城市配網(wǎng)電纜的 R/X 能到 1 以上解耦前提不成立這時(shí)快速解耦會(huì)震蕩甚至發(fā)散應(yīng)換回牛頓法。這也是為什么程序默認(rèn)--method nr而不默認(rèn) fdpf。高斯-賽德?tīng)柗ㄖ蛔鳛槎档缀徒虒W(xué)用收斂慢但每輪代價(jià)極低、內(nèi)存占用小適合幾千節(jié)點(diǎn)的稀疏配網(wǎng)先跑個(gè)粗略初值再把結(jié)果作為牛頓法的熱啟動(dòng)。這套代碼里它就是一個(gè)簡(jiǎn)單循環(huán)def gauss_seidel_pf(Ybus, Sbus, V0, case, tol1e-8, max_iter500): Y Ybus.toarray() V V0.astype(complex).copy() for it in range(max_iter): for i in range(len(V)): if i in case.ref or i in case.pv: continue # 把第 i 行的自導(dǎo)納項(xiàng)單獨(dú)解出來(lái) Ii Sbus[i] / np.conj(V[i]) V[i] (Ii - Y[i, :] V Y[i, i] * V[i]) / Y[i, i] if np.max(np.abs(V - V0)) tol: return V, it V0 V.copy() raise RuntimeError(高斯-賽德?tīng)栁词諗?三種解法的選擇經(jīng)驗(yàn)是小于 300 節(jié)點(diǎn)、要精確結(jié)果用nr高壓網(wǎng)、上千節(jié)點(diǎn)、要批量掃方式用fdpf數(shù)據(jù)質(zhì)量差、初值離譜、只是想先看看電壓大概形態(tài)用gs熱啟動(dòng)再切nr。程序里--method就是這三個(gè)字符串。4. 跑起來(lái)環(huán)境、命令行參數(shù)與結(jié)果校驗(yàn)算法講清楚之后得先讓它跑出第一組數(shù)否則后面全是紙上談兵。這一章按“裝環(huán)境 → 跑單個(gè)算例 → 批量掃 → 對(duì)參考解”的順序來(lái)每步都能獨(dú)立驗(yàn)證。4.1 依賴(lài)安裝與最小可跑路徑依賴(lài)很輕numpyscipy是核心pandas只用于報(bào)表導(dǎo)出pyyaml讀配置matplotlib可選畫(huà)收斂曲線。建議開(kāi)虛擬環(huán)境別裝到系統(tǒng) Python 里。# 1. 建環(huán)境Python 3.9 以上都行3.10/3.11 最穩(wěn) python -m venv .venv source .venv/bin/activate # Windows: .venv\Scripts\activate # 2. 裝依賴(lài) pip install numpy scipy pandas pyyaml # 或者 pip install -r requirements.txt # 3. 自檢跑測(cè)試用例先確認(rèn)算法層是好的 python -m pytest tests/ -q # 4. 跑第一個(gè)算例 python run_pf.py --case cases/case14.json --method nr --verbose第 3 步容易被跳過(guò)但它是排錯(cuò)時(shí)的分水嶺測(cè)試通過(guò)說(shuō)明算法和 Ybus 組裝沒(méi)問(wèn)題后面跑不通就是數(shù)據(jù)或參數(shù)的問(wèn)題測(cè)試不通過(guò)就先別碰算例文件。tests/test_newton_case14.py里通常會(huì)硬編碼兩三個(gè)節(jié)點(diǎn)的參考電壓跑一下就知道牛頓法有沒(méi)有被改壞。常見(jiàn)的新手卡點(diǎn)是ModuleNotFoundError: No module named pf原因是沒(méi)在項(xiàng)目根目錄執(zhí)行或者run_pf.py里的sys.path沒(méi)把當(dāng)前目錄加進(jìn)去。簡(jiǎn)單辦法是在根目錄執(zhí)行或者pip install -e .裝成可編輯包。4.2 三種入口單算例、批量掃描、配置驅(qū)動(dòng)run_pf.py的參數(shù)不多但每個(gè)都影響結(jié)果值得對(duì)著表過(guò)一遍參數(shù)默認(rèn)值說(shuō)明--case無(wú)單個(gè)算例文件路徑JSON/CSV--case-dir無(wú)批量模式目錄下所有算例逐個(gè)求解--methodnrnr/fdpf/gs--tol1e-8功率不平衡量收斂閾值p.u.--max-iter20迭代上限--sweep-load無(wú)負(fù)荷倍率掃描格式起:止:步長(zhǎng)--outresults輸出目錄--verbosefalse打印每輪最大不平衡量單算例是最常用的python run_pf.py --case cases/case30.json --method nr \ --tol 1e-8 --max-iter 20 \ --out results/case30 --verbose加--verbose后終端會(huì)打出類(lèi)似iter 1: 2.31e-01 / iter 2: 8.42e-03 / iter 3: 1.9e-05 / iter 4: 6.1e-10的序列這條曲線比任何日志都好用正常應(yīng)該是階梯式下降每次大約掉兩個(gè)數(shù)量級(jí)如果第 2 次就掉到 1e-9八成是算例太簡(jiǎn)單或者容差被改過(guò)如果一直卡在 1e-3去看第 5 章的越限問(wèn)題。負(fù)荷掃描是規(guī)劃場(chǎng)景里最實(shí)用的功能一行命令能出一族結(jié)果python run_pf.py --case cases/case14.json --sweep-load 0.8:1.2:0.1 \ --method fdpf --out results/sweep14它會(huì)把負(fù)荷倍率從 0.8 到 1.2 按 0.1 步長(zhǎng)各算一次每次結(jié)果單獨(dú)存目錄最后在summary.json里匯總各方式的收斂情況和最低電壓。做重載校核時(shí)這條命令比手工改十遍文件靠譜得多。批量模式適合回歸驗(yàn)證python run_pf.py --case-dir cases/ --method nr --out results/batch # 任何一個(gè)算例不收斂進(jìn)程返回非零退出碼方便掛到 CI 上4.3 輸出報(bào)表與跟參考解對(duì)齊的方法跑完之后results/case30/下會(huì)有三張表和一份匯總results/case30/ ├── bus_result.csv # bus_id, Vm, Va_deg, Pg, Qg, Pd, Qd ├── branch_result.csv # f, t, Pf, Qf, Pt, Qt, Ploss, Qloss, loading ├── summary.json # 迭代次數(shù)、末次誤差、總網(wǎng)損、最大負(fù)載率 └── converge.log # 每輪不平衡量便于畫(huà)圖branch_result.csv里的Ploss/Qloss是兩端功率之差loading是max(|Pf|,|Pt|)/rateA這條列直接決定你要不要回頭改網(wǎng)架。所有報(bào)表的功率單位是 MVA/MW電壓是 p.u. 加 kV 兩列省得自己乘baseKV。結(jié)果到底可不可信靠對(duì)參考解。ref/case14_voltage_ref.csv存了 14 節(jié)點(diǎn)測(cè)試系統(tǒng)的標(biāo)準(zhǔn)電壓解寫(xiě)一段校驗(yàn)?zāi)_本import json import numpy as np import pandas as pd ref pd.read_csv(ref/case14_voltage_ref.csv) # 列bus_id, vm_ref, va_ref cur pd.read_csv(results/case14/bus_result.csv) err_v np.max(np.abs(cur[Vm].values - ref[vm_ref].values)) err_a np.max(np.abs(cur[Va_deg].values - ref[va_ref].values)) print(f電壓幅值最大偏差 {err_v:.3e} p.u.相角最大偏差 {err_a:.3e} 度) assert err_v 1e-6, 電壓幅值對(duì)不上優(yōu)先查標(biāo)幺化和 Ybus assert err_a 1e-4, 相角對(duì)不上優(yōu)先查變壓器移相角和變比方向 # 功率平衡校驗(yàn)總注入 總負(fù)荷 總網(wǎng)損精度取決于 tol loss pd.read_csv(results/case14/branch_result.csv)[[Ploss, Qloss]].sum() summary json.load(open(results/case14/summary.json)) gap abs(summary[total_ploss] - loss[Ploss]) assert gap 1e-4, f網(wǎng)損匯總不一致差 {gap:.3e} MW這兩條斷言是我每次改完代碼都會(huì)跑的。第一條卡電壓第二條卡網(wǎng)損兩者同時(shí)通過(guò)基本可以認(rèn)定結(jié)果可信。相角容差放到1e-4度是因?yàn)椴煌瑢?shí)現(xiàn)的角度基準(zhǔn)可能差一點(diǎn)點(diǎn)幅值容差1e-6p.u. 則沒(méi)什么商量余地。注意參考解文件里的節(jié)點(diǎn)順序必須和算例一致。有些參考文件按bus_id排序而算例里是亂序的直接按行號(hào)比較會(huì)得到一堆莫名其妙的偏差。穩(wěn)妥做法是先merge到bus_id上再比。5. 避坑與排查五個(gè)讓我返工過(guò)的問(wèn)題潮流計(jì)算程序的坑有一個(gè)共同特征程序不報(bào)錯(cuò)結(jié)果看著也挺合理。下面五條都是我自己或者身邊的人真實(shí)踩過(guò)的按“現(xiàn)象 → 原因 → 解決”寫(xiě)遇到了可以對(duì)號(hào)入座。5.1 基準(zhǔn)值不一致電壓算出來(lái)是 0.4 而不是 1.0現(xiàn)象第一個(gè)算例跑完bus_result.csv里一半節(jié)點(diǎn)的Vm在 0.3 到 0.5 之間但收斂正常、迭代 4 次、網(wǎng)損也在合理范圍。原因標(biāo)幺化做了兩次。caseio里已經(jīng)按baseMVA100把功率折算成 p.u.而run_pf.py在傳給求解器本文還有配套的精品資源點(diǎn)擊獲取