簡介Power System ToolboxPST是一套運行于MATLAB環(huán)境的電力系統(tǒng)專業(yè)工具箱面向科研人員、工程師與學生用于解決復雜電力網(wǎng)絡的建模、穩(wěn)態(tài)/動態(tài)/暫態(tài)仿真及控制策略驗證問題。壓縮包共有211個文件其中209個m腳本/函數(shù)和2個pdf文檔整體僅1.1MBm文件涵蓋設備模型、仿真主程序、示例數(shù)據(jù)與繪圖工具pdf文檔提供使用與部署說明。核心功能支持發(fā)電機、變壓器、線路、負荷、斷路器、繼電保護裝置等模型詳細程度可從簡化靜態(tài)到復雜動態(tài)靈活選擇能夠模擬正常運行時的電壓、功率分布及故障后的系統(tǒng)穩(wěn)定性。工具箱還提供勵磁控制器、調速器、無功補償器等控制器的設計與測試環(huán)境配合波形圖、功率流圖、頻率響應曲線等可視化圖表并支持用戶自定義模型和Simulink集成便于新能源并網(wǎng)、保護配置優(yōu)化等應用擴展。目前已有761人學習下載借助內(nèi)置腳本和示例數(shù)據(jù)電力系統(tǒng)專業(yè)師生或從業(yè)者可快速建立仿真環(huán)境開展教學實驗與課題研究。1. Power System Toolbox 到底是干什么的三個問題幫你確認它值不值得裝Power System Toolbox 這類 MATLAB 仿真工具集最容易被誤解成“又一個畫潮流圖的插件”。實際上它是一套面向電力系統(tǒng)動態(tài)分析的命令行框架能算潮流、能搭動態(tài)模型、能跑暫態(tài)穩(wěn)定還能把你自己的控制策略寫成一段回調函數(shù)塞進仿真循環(huán)里。所以先別急著找圖形界面——如果你關心的是“我這個控制器掛到雙機系統(tǒng)上故障切除后轉子還能不能穩(wěn)住”那它就比一堆只做靜態(tài)潮流計算的工具箱對路得多。它適合三類人做機電暫態(tài)算法驗證的研究生做保護或控制策略研發(fā)的工程師以及要在論文里復現(xiàn)“某個典型系統(tǒng)”的仿真曲線、又不想從零手寫微分方程求解器的開發(fā)者。不適合只想要一個“點一下出一張單線圖”的可視化工具的人。反直覺的一點是這套工具真正值錢的地方不是它自帶的幾個測試系統(tǒng)而是那套“矩陣定義系統(tǒng)、函數(shù)驅動計算、回調注入控制”的松耦合結構。你的核心資產(chǎn)是算例數(shù)據(jù)和模型參數(shù)不是界面。2. 從數(shù)據(jù)格式入手母線、支路、發(fā)電機矩陣怎么排才不出錯2.1 bus 矩陣母線編號、電壓和功率注入怎么排這類工具箱最通用的數(shù)據(jù)組織方式是把整個電力系統(tǒng)拆成若干組固定列數(shù)的矩陣。第一組是母線矩陣每一行代表一條母線每一列對應一個物理量。列的順序一般沿襲一套歷史約定母線編號、電壓幅值初值、相角初值、有功負荷、無功負荷、并聯(lián)導納、所屬區(qū)域編號、節(jié)點類型以及電壓約束上下限。% 構造一個三節(jié)點測試系統(tǒng)的 bus 矩陣 % 列含義編號, V0(pu), theta0(rad), P_load(pu), Q_load(pu), Gsh(pu), Bsh(pu), 區(qū)域, Vmax, Vmin bus [ 1 1.00 0.00 1.20 0.40 0 0 1 1.10 0.90; % 節(jié)點1PQ負荷1.20.4j 2 1.00 0.00 0.50 0.20 0 0 1 1.10 0.90; % 節(jié)點2小負荷 3 1.00 0.00 0.30 0.10 0 0 1 1.10 0.90; % 節(jié)點3接發(fā)電機 ];這里每個數(shù)字都是標幺值pu基準容量一般取全系統(tǒng)統(tǒng)一的 MVA 值比如 100MVA。第 1、2 列的電壓和相角是潮流迭代的初值不是約束值真正的電壓約束在最后兩列。新手最容易踩的坑是把“初值”當成“給定值”導致發(fā)電機母線電壓永遠迭代不動。2.2 line 矩陣R/X/B/變比對應哪一列第二組是支路矩陣每條線路或每臺變壓器一行。列順序有很強的一致性支路編號、首端母線號、末端母線號、串聯(lián)電阻 R、串聯(lián)電抗 X、對地電納 B、變壓器變比以及相移。單位同樣是標幺值R/X 是基于系統(tǒng)基準值的串聯(lián)阻抗B 是線路充電電納的一半還是全量不同工具箱習慣不一但最常見的是全量電納除以 2 再填入。% 三節(jié)點系統(tǒng)的三條支路 % 列含義編號, from, to, R(pu), X(pu), B/2(pu), 變比k, 相移(度) line [ 1 1 2 0.020 0.060 0.030 0 0.0; % 線路1-2 2 2 3 0.030 0.090 0.040 0 0.0; % 線路2-3 3 3 1 0.040 0.120 0.050 1 0.0; % 變壓器3-1變比1.0 ];注意變壓器和線路混在同一張表里時變比列的填法有講究如果是標幺變比 1.0直接填 1如果是非標幺變比比如 10.5/121kV 的實際抽頭需要折算到基準電壓之后再填。X 為負不一定是你填錯了可能是串補電容但如果 R 和 X 填反潮流會直接給你一個“換相失敗”式的離譜相角。2.3 發(fā)電機與動態(tài)模型矩陣潮流模型和動態(tài)模型的銜接發(fā)電機矩陣是這兩類模型之間的橋梁。潮流計算只需要發(fā)電機母線號、有功出力、無功出力、有功上下限、無功上下限這幾個字段而暫態(tài)穩(wěn)定仿真還需要給出發(fā)電機的動態(tài)參數(shù)比如慣性常數(shù)、暫態(tài)電抗、阻尼系數(shù)以及勵磁系統(tǒng)和調速器對應的模型編號。常見做法是分兩張表一張管潮流注入一張管動態(tài)參數(shù)中間用母線號關聯(lián)。% 發(fā)電機潮流數(shù)據(jù)母線號, Pg(pu), Qg(pu), Pmax, Pmin, Qmax, Qmin gen [ 3 0.80 0.20 1.20 0.00 0.80 -0.40; % 節(jié)點3的平衡機 ]; % 發(fā)電機動態(tài)參數(shù)母線號, H(s), Xd(pu), Xq(pu), D, 模型階數(shù) mac [ 3 6.50 0.30 0.35 2.0 2; % 2階經(jīng)典模型 ];mac 表的模型階數(shù)直接決定后面積分器要不要解勵磁繞組暫態(tài)。2 階經(jīng)典模型只保留轉子運動方程適合看第一擺穩(wěn)定性想研究阻尼特性至少要 4 階以上模型把勵磁回路動態(tài)算進去。如果之后發(fā)現(xiàn)暫態(tài)波形衰減得過分平滑先回來檢查 D 是不是被填成了 0。3. 最小潮流計算跑通三組命令把結果拿到手3.1 準備算例文件把三組矩陣寫進工作區(qū)這類工具箱通常不要求你把數(shù)據(jù)存成特定格式的文件只要在 MATLAB 工作區(qū)里有bus、line、gen三組矩陣就能直接開工。最穩(wěn)妥的落地方式是把系統(tǒng)數(shù)據(jù)寫進一個獨立的腳本文件每次運行先執(zhí)行這個腳本再調求解器。這樣算例與求解分離換系統(tǒng)時只替換數(shù)據(jù)文件不碰主循環(huán)。% case3.m —— 三節(jié)點測試系統(tǒng)數(shù)據(jù)腳本 bus [ ... ]; % 粘貼 2.1 節(jié)的三行母線數(shù)據(jù) line [ ... ]; % 粘貼 2.2 節(jié)的三行支路數(shù)據(jù) gen [ ... ]; % 粘貼 2.3 節(jié)的發(fā)電機數(shù)據(jù)腳本文件沒有返回值運行之后矩陣直接進工作區(qū)后續(xù)函數(shù)按變量名取數(shù)據(jù)。變量名必須嚴格一致否則求解器會報“找不到匹配的 bus 變量”。養(yǎng)成一個習慣任何從外部導入的數(shù)據(jù)無論是文本還是 Excel都先轉成這三組矩陣再做后續(xù)計算不要試圖繞過這套數(shù)據(jù)結構。3.2 跑起來牛頓法還是高斯賽德爾怎么選以這套框架的常見命令風格為例潮流求解的主函數(shù)只需要傳入三組矩陣外加兩個關鍵選項迭代算法和收斂容差。牛頓法收斂快、適合病態(tài)系統(tǒng)但對初值敏感高斯賽德爾法實現(xiàn)簡單、初值要求低但大系統(tǒng)里收斂慢。實際工程里我一般直接用牛頓法容差設在 1e-6最大迭代次數(shù)給 30 次。% 運行潮流計算 [V, converged, iter] psflow(bus, line, gen, ... algorithm, newton, ... % 可選 newton 或 gauss tol, 1e-6, ... % 收斂判據(jù)電壓修正量范數(shù) maxiter, 30); % 迭代上限超過即判不收斂 % 檢查是否收斂 if ~converged error(潮流不收斂迭代 %d 次后仍未滿足精度, iter); end收斂判據(jù)這里的含義是相鄰兩次迭代的電壓修正向量范數(shù)小于 1e-6。設得越嚴計算越準但對病態(tài)電網(wǎng)越容易觸發(fā)上限退出。如果發(fā)現(xiàn)不收斂與其把容差調松到 1e-4不如先檢查負荷和發(fā)電機出力是否嚴重不平衡——多數(shù)不收斂是初值離解太遠不是算法不行。3.3 讀結果電壓、相角、線路潮流分別在哪潮流函數(shù)一般返回一個復數(shù)電壓向量V實部虛部分別對應直角坐標下的電壓分量或者用abs(V)和angle(V)轉成幅值與相角。要看線路有功和無功流向需要再用一組函數(shù)基于V和line矩陣回算線路功率。這里最容易踩坑的是V里元素的順序是按母線編號排的還是按 bus 矩陣的行號排的。兩者在多數(shù)算例里一致但如果你手動刪過 bus 某一行索引就會錯位。% 提取幅值和相角 Vmag abs(V); % 各母線電壓幅值單位 pu Vang angle(V) * 180/pi; % 各母線相角單位 度 % 線路功率回算返回首末端有功/無功 [Pf, Qf, Pt, Qt] lineflow(bus, line, V); % 打印節(jié)點1的電壓和節(jié)點1-2線路有功 fprintf(V1 %.4f pu, theta1 %.2f deg\n, Vmag(1), Vang(1)); fprintf(P12 %.4f pu\n, Pf(1));線路功率回算是基于節(jié)點電壓和線路阻抗支路的復功率公式方向定義為首端流向末端。做校核時有個快速粗查辦法所有線路首端有功之和加上所有發(fā)電機出力必須等于所有負荷有功加全系統(tǒng)網(wǎng)絡損耗誤差超過 0.5% 就說明某個矩陣填錯了。4. 把潮流結果接進暫態(tài)穩(wěn)定仿真三個必調參數(shù)4.1 仿真步長與總時長先定這兩個數(shù)暫態(tài)穩(wěn)定仿真的總時長和步長不是越大越好也不是越小越好。總時長要覆蓋“第一擺”到“第二擺”的完整過程經(jīng)典機電暫態(tài)取 5 到 10 秒就夠了研究弱阻尼或低頻振蕩要放到 20 秒。步長則取決于積分算法和你關心的最高頻率分量梯形隱式積分法在機電暫態(tài)尺度下一般取 0.01 秒極限試驗時可以到 0.005 秒再小就是純浪費時間。% 暫態(tài)穩(wěn)定求解常見命令風格示例 [t, x, Vt] psstab(bus, line, gen, mac, ... tend, 10.0, ... % 仿真總時長單位秒 dt, 0.01, ... % 固定步長單位秒 event, event_seq); % 故障事件序列見 4.2仿真總時長和事件時間有一個隱蔽關系如果事件發(fā)生時刻不在步長網(wǎng)格上求解器要么就近取整、要么插值兩種處理都會讓波形產(chǎn)生虛假抖動。我自己會先把所有事件時刻對齊到步長的整數(shù)倍比如 1.00、1.10、1.30而不要寫 1.05、1.12 這種步長整除不了的數(shù)。4.2 故障序列接地、斷線、重合閘的事件表事件序列是暫態(tài)仿真的核心輸入通常用一張表來定義每一行是一個事件列包含觸發(fā)時間、動作類型、作用母線。以常見實現(xiàn)為例三相短路接地、故障清除、線路切除、重合閘都是通過修改導納矩陣來模擬的動作的物理含義是“把某條母線對地導納改變到多少”。% 事件序列t1.0s 母線3三相短路1.1s清除1.3s重合閘 event_seq [ 1.00 3 3 -1 0.0 1; % t1.0母線3經(jīng)0.0歐接地三相短路 1.10 0 0 -1 0.0 0; % t1.1清除故障恢復正常拓撲 1.30 3 3 -1 0.0 1; % t1.3斷路器重合 ];事件表的時間列必須嚴格單調遞增同一個時刻只能有一個事件。想模擬單相接地需要把三相短路的“接地”改成“某相接地”并重新設置不對稱系數(shù)多數(shù)機電暫態(tài)框架默認只做正序網(wǎng)絡單相和三相短路的差異主要體現(xiàn)在故障期間的正序附加阻抗上不要過份追求不對稱細節(jié)。4.3 動態(tài)模型參數(shù)哪些是必調哪些可以默認發(fā)電機動態(tài)模型參數(shù)表中慣性常數(shù) H、暫態(tài)電抗 Xd、負序/次暫態(tài)參數(shù)、阻尼系數(shù) D 是最先要確認的四個。經(jīng)典 2 階模型只需要 H、Xd 和 D4 階模型要補勵磁繞組時間常數(shù) Td06 階模型還要再加交軸阻尼繞組參數(shù)。勵磁系統(tǒng)的模型參數(shù)更敏感——調壓器放大倍數(shù)設大了暫態(tài)過程會高頻振鈴設小了電壓恢復慢。參數(shù)典型范圍調大后的現(xiàn)象調小后的現(xiàn)象慣性常數(shù) H (s)3.0 ~ 8.0轉子擺動變慢第一擺幅度變小擺動變快更容易失穩(wěn)阻尼系數(shù) D1.0 ~ 5.0振蕩衰減加快穩(wěn)態(tài)偏差變大振蕩衰減慢持續(xù)搖擺暫態(tài)電抗 Xd (pu)0.2 ~ 0.4極限切除時間變短更難穩(wěn)定穩(wěn)定裕度增大勵磁放大倍數(shù) K_A20 ~ 200電壓響應快可能高頻振蕩電壓恢復慢動態(tài)壓降大這些參數(shù)在真實發(fā)電機銘牌和試驗報告里都有據(jù)可查做研究時建議從典型值起步不要一上來就調成極端值。一個很容易被忽視的細節(jié)是調速器模型參數(shù)中調差系數(shù)對潮流結果有影響但對暫態(tài)第一擺影響不大如果你只關心故障后 2 秒內(nèi)的行為調速器用默認參數(shù)即可把精力花在勵磁和發(fā)電機本體的參數(shù)上。5. 避坑指南暫態(tài)穩(wěn)定仿真最容易翻車的 4 個數(shù)據(jù)坑5.1 潮流不收斂但程序不報錯最后全輸出 1.0 pu現(xiàn)象潮流函數(shù)返回converged true但所有母線電壓都是 1.0相角全是 0和物理常識完全不符。原因典型的“初值等于解”的假收斂。當負荷和發(fā)電機出力都被填成 0或者系統(tǒng)基準 MVA 和發(fā)電機數(shù)據(jù)基準 MVA 不一致時比如發(fā)電機數(shù)據(jù)按 200MVA 基準填寫系統(tǒng)基準卻是 100MVA迭代方程在原點附近自然收斂但這個解沒有承載任何功率。解決先用功率平衡粗查。把全系統(tǒng) P_load 加起來減去所有 Pg差值應該在網(wǎng)損量級如果差值幾乎為 0大概率是基準值錯了。把發(fā)電機有功整體乘 0.5 再跑一次如果電壓和相角有明顯變化說明數(shù)據(jù)通路是通的只是基準需要統(tǒng)一。5.2 暫態(tài)波形在故障清除時刻出現(xiàn)鋸齒形突變現(xiàn)象搖擺曲線在 1.1 秒附近出現(xiàn)一個很窄的脈沖或凹坑看起來像數(shù)值振蕩但步長已經(jīng)很小了。原因事件發(fā)生時間和步長網(wǎng)格不重合。比如步長 0.01 秒事件寫的是 1.105 秒積分器會就近取整到 1.10 或 1.11再強行觸發(fā)代數(shù)變量跳變形成的波形就帶細鋸齒。這在重合閘時刻尤其明顯因為斷路器動作本身就會造成一個不連續(xù)的拓撲切換。解決把每個事件時刻都對齊到步長的整數(shù)倍。我通常先定步長再反推事件時間步長 0.01 秒事件就寫 1.00、1.10、1.30需要更精細的故障時刻就把步長改成 0.005 秒讓 1.105 變得合法。另一個辦法是檢查求解器有沒有事件插值開關開了反而會引入新的數(shù)值阻尼慎用。5.3 線路電抗為負功率方向全亂現(xiàn)象潮流算出來某條線路的有功方向是從負荷側流向電源側且該線路的 X 列是負數(shù)。原因兩種可能。一種是 X 列本來填的就是串補電容的容抗為負正常另一種更常見——R 和 X 填反了或者 bus/line 矩陣里的列順序和文檔對不上導致程序把某列當成了 X。這類錯誤最陰險因為潮流照樣收斂只是結果不物理。解決先把該線路單獨拿出來算一遍手算值已知兩端電壓幅值和相角差手算線路有功傳輸公式 P V1*V2/X * sin(theta1 - theta2)。如果手算是正方向而程序結果是負方向就是矩陣列錯位。我會把 bus 和 line 矩陣的表頭注釋永遠留在腳本第一段不刪除、不壓縮切換數(shù)據(jù)源時先拿一個已知算例驗證列順序。5.4 動態(tài)仿真初始功率和潮流結果對不上現(xiàn)象潮流計算顯示發(fā)電機出力是 0.8 pu但暫態(tài)仿真的第一幀結果里同一條母線的注入功率變成了 0.75 pu二者對不上。原因動態(tài)模型初始化沒有正確吃進潮流結果。發(fā)電機的內(nèi)電勢、勵磁電壓初值、調速器閥門開度都需要從潮流結果反推常見實現(xiàn)里這一步叫“初值計算”或“初始化”。如果動態(tài)參數(shù)表里的 Xd 與潮流計算用的 Xd 不一致反推出來的內(nèi)電勢就會帶誤差進而影響第一幀功率。解決檢查初始化函數(shù)是否基于同一份潮流結果。我的習慣是先跑潮流把 V 作為初始化函數(shù)的輸入不要初始化和潮流各算各的。另一個自查辦法把暫態(tài)仿真的前 0.1 秒結果取平均正常情況下它應該非常接近潮流結果如果差得明顯優(yōu)先懷疑動態(tài)電抗參數(shù)不一致而不是懷疑求解器。調試時把這個差值打印出來比盯著整條曲線找問題快得多。6. 進階用法做出自己的算例并用結果反向驗證6.1 把自己的數(shù)據(jù)寫成算例文件如果你拿到的是某電網(wǎng)或某設備廠商的導出數(shù)據(jù)常見格式是帶表頭的 CSV 或 Excel。不要手動復制進矩陣——用腳本自動拼。先讀表頭確認列名再按工具箱要求的列順序重新組織輸出成caseN.m腳本。關鍵步驟是給每列做一次單位折算發(fā)電機容量基準、電壓基準、阻抗基準三者必須一致。我會習慣在數(shù)據(jù)文件頭部用注釋寫明“基準容量 100MVA基準電壓按各區(qū)平均電壓取”防止三個月后回來看不懂。% 從CSV自動生成bus矩陣的示例邏輯 csv readtable(network.csv); bus zeros(height(csv), 10); bus(:,1) csv.bus_number; bus(:,3) (csv.load_kW / 100000) * exp(0); % kW - 標幺P bus(:,4) (csv.load_kvar / 100000); % kvar - 標幺Q % 列2電壓初值默認1.0列9/10為電壓限值默認1.1/0.96.2 三個自查手段功率平衡、電壓越限、曲線單調性做完仿真先別畫圖用三個腳本內(nèi)聯(lián)檢查代替肉眼全網(wǎng)功率平衡檢查、電壓越限掃描、搖擺曲線局部單調性檢查。功率平衡直接決定潮流數(shù)據(jù)可信度電壓越限掃描幫你發(fā)現(xiàn)會不會有接近失穩(wěn)的邊界曲線單調性則是數(shù)值振蕩的信號——真實機電暫態(tài)不會出現(xiàn)持續(xù)等幅鋸齒出現(xiàn)就說明步長或初值有問題。% 電壓越限掃描 Vfinal abs(Vt(end,:)); if any(Vfinal 1.1) || any(Vfinal 0.9) warning(存在電壓越限節(jié)點max%.3f pu, min%.3f pu, ... max(Vfinal), min(Vfinal)); end % 第一擺單調性轉子角前0.5s內(nèi)不應出現(xiàn)頻繁反向 theta1 angle(Vt(:,1)); % 取1號機轉子角隨時間序列 dtheta diff(theta1(1:round(0.5/t))); if sum(dtheta 0) length(dtheta) * 0.8 warning(轉子角前0.5s內(nèi)有劇烈反向疑似數(shù)值振蕩); end我最早做這類仿真的時候習慣拿到結果直接畫曲線眼睛看著差不多就完事。后來在一次故障時間參數(shù)微調中吃了虧——曲線整體變好看了但功率平衡檢查差了 3%追查半天才發(fā)現(xiàn)是某條線路的 R 填錯了一位小數(shù)。從那以后我把這三個檢查寫進每次仿真的收尾腳本里不改數(shù)據(jù)、不重跑結果不合格就拒絕看曲線。這算不上什么聰明技巧但確實讓我的仿真結果可信度上了一個臺階。希望這篇筆記能幫你少走這段彎路把精力花在真正要研究的控制策略和穩(wěn)定性問題上。本文還有配套的精品資源點擊獲取