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