企业官网建设流程全解析

简介：Power System Toolbox（PST）是一套面向MATLAB环境开发的电力系统建模、分析与仿真工具包，适合电力专业研究人员、工程师及高校学生用于稳态/暂态计算、故障仿真、控制器设计与新能源并网分析等工作。资源共211个文件，以m脚本为主，含少量pdf说明文档，压缩包仅1.1MB，便于下载后快速部署。m文件覆盖发电机、变压器、线路等设备模型及多种仿真分析函数，可支撑从简单静态计算到复杂动态仿真的研究需求。目前已吸引762人学习，适合正在开展电力系统课程设计、稳定性评估或控制策略研究的用户参考。通过研读PST的源码结构与仿真流程，读者能掌握MATLAB下电力系统自定义建模、批量仿真与结果可视化的实现思路，进一步提升解决实际电力工程问题的效率。

1. Power System Toolbox：藏在MATLAB里的电力系统仿真利器

Power System Toolbox（PST）是一套基于MATLAB的电力系统动态仿真工具箱，很多人从PSASP或BPA入门电力系统仿真，等到要做自定义模型、批量参数扫描时才发现商业软件的黑匣子让人束手束脚。PST把潮流计算、暂态稳定时域仿真、小信号特征值分析全都以可读的.m源文件形式开放出来，适合研究生、开发者和算法验证工程师。这里以一次三相短路暂态仿真为主线，串起PST的数据格式、初始化流程、仿真参数与结果解读，最后把实际踩过的坑逐个拆开。

2. PST数据与初始化：母线、支路、发电机三类矩阵怎么填

2.1 为什么选PST：与商业仿真软件的核心差异

在决定用PST之前，先得搞清楚它的定位。商业软件的优势是模型库全、计算引擎成熟、撑得住大规模系统，但代价是自定义模型困难，很多内部迭代细节不开放，批量仿真还要依赖脚本接口，版本升级后模型文件兼容性焦虑总是存在。PST最大的价值在于透明：从一个数据文件到潮流结果，再到暂态曲线和特征值，每一步都能对应到源文件里的具体数学操作。

PST适合的典型场景有两类：一类是教学与科研场景，需要把暂态稳定、励磁系统、调速系统等背后的状态方程摊开来看；另一类是算法预研场景，新的控制策略要在仿真里反复改模型、改参数，商业软件里改一个环节要过层层封装，PST里直接改一个.m文件就行。代价是需要自己维护数据格式，对使用者有基础要求，至少得懂潮流计算和状态空间的基本概念，否则出了问题很难排查。

PST支持的系统规模通常在几十台机到几百台机的量级，和商业软件数千台机没法比，但用于方法验证和机理分析绰绰有余。从另一个角度看，PST与MATLAB生态无缝对接，后续做参数辨识、优化算法、结果绘图都能直接复用MATLAB的数据结构和算法库，这种便利性是外部仿真软件难以替代的。

2.2 数据文件结构：三类矩阵的行列语义

PST的数据文件本质上是一个MATLAB脚本，运行后在工作区生成一系列全局变量。最核心的三类是bus、branch、gen，分别描述母线、支路和发电机注入信息。这三类矩阵的列顺序有严格约定，填错一列后面全乱。

bus矩阵每行是一条母线，列含义依次是母线编号、母线类型（1摆动母线、2 PV母线、3 PQ母线）、电压幅值初始值（pu）、电压相角初始值（弧度）、有功负荷（MW）、无功负荷（Mvar）、发电机有功出力（MW）、发电机无功出力（Mvar）、无功补偿（Mvar）。母线类型决定了潮流计算里该母线的边界条件：摆动母线电压幅值和相角固定，PV母线固定电压幅值和有功出力，PQ母线固定负荷。

branch矩阵每行是一条支路，列依次是起始母线编号、终止母线编号、串联电阻R（pu）、串联电抗X（pu）、对地电纳B/2（pu）、变压器非标准变比、移相角等。支路参数的标幺值必须与数据文件基准容量一致，常见坑是原始数据从别的仿真格式转过来，基准容量不一致却直接照搬数值，导致潮流计算结果完全失真。

gen矩阵描述母线处发电机的无功上下限，列依次是母线编号、无功出力上限、无功出力下限。除了这三类，数据文件还可能包含mac_con（发电机动态参数）、exc_con（励磁系统）、pss_con（电力系统稳定器）等矩阵，这些只有在动态仿真初始化时才真正参与计算。所以如果只做稳态潮流分析，数据文件可以只写bus、branch、gen三块。

给一个最小示例，实际抄写数据时注意矩阵中每一列用空格分隔，负号前面留出空隙以避免MATLAB把负数当运算符处理。

% 最小数据文件示例（仅潮流计算所需） basmva = 100; % 基准容量 100MVA bus = [ 1 1 1.00 0.00 0.0 0.0 50.0 20.0 0.0 ; 2 2 1.00 0.00 20.0 10.0 80.0 30.0 0.0 ; 3 3 1.00 0.00 40.0 20.0 0.0 0.0 0.0 ; ]; branch = [ 1 2 0.010 0.10 0.02 0.0 0.0 ; 2 3 0.015 0.15 0.03 0.0 0.0 ; 1 3 0.020 0.20 0.04 0.0 0.0 ; ]; gen = [ 1 50.0 0.0 ; 2 80.0 0.0 ; ];

示例里的母线3是纯负荷节点，所以gen矩阵不包含母线3。basmva=100表示所有标幺值都以100MVA为基准，原始数据若是基准容量1000MVA的系统，需要先按容量比统一换算再填入矩阵。换算的思路很简单：新标幺值=旧标幺值×（旧基准容量/新基准容量），支路电阻、电抗、电纳全都按这条规则处理。每次迁移数据，我会把原始系统的总出力除以总负荷，检查比值是否接近1，这能快速暴露容量换算错误。

2.3 初始化流程：从数据文件到可仿真模型

拿到一份数据文件后，标准的三步初始化流程在PST里非常固定：加载数据、潮流计算、发电机模型初始化。

% 步骤1：加载数据文件 clear; clc; global basmva bus branch gen mac_con datafile = 'case9.m'; % 9节点测试系统数据文件 run(datafile); % 执行数据脚本，生成全局变量 % 步骤2：潮流计算 tol = 1e-6; % 收敛精度 itermax = 30; % 最大迭代次数 [bus, branch, gen] = pf(bus, branch, gen, tol, itermax); % 步骤3：初始化发电机动态模型 mac_con = mac_ini(bus, gen);

第一步中run(datafile)直接执行数据脚本，脚本内部会把全局变量填充好。注意必须先声明global变量名，否则数据文件里赋的值只在函数工作区生效，主函数根本看不到。第二步的pf是PST的潮流计算函数，按bus、branch、gen的顺序迭代求电压幅值和相角。tol和itermax分别是收敛精度和最大迭代次数，对9节点测试系统tol=1e-6、itermax=30通常一次就过。如果系统规模几百台机，itermax可以放到50到100，但调迭代次数是治标不治本，真正的问题多半是数据本身有矛盾。

第三步mac_ini把潮流结果和发电机的动态参数矩阵mac_con联立起来，计算初始转子角、初始暂态电动势等状态变量。这一步是时域仿真的前置条件，跳过它直接调用仿真函数，PST会报状态变量维度不一致的错误。值得注意的是，mac_ini要求工作区里已经存在mac_con矩阵，如果数据文件里没有定义这个矩阵，运行mac_ini会直接报变量未定义。

注意：mac_con矩阵必须在mac_ini调用之前就存在于工作区，且行数要与gen矩阵的行数一致。检查时先打开数据文件末尾，确认动态参数矩阵是否齐全。

完成这三步后，工作区里已经有一组自洽的稳态工作点。我自己的习惯是每次初始化完后，先把bus矩阵第3列的电压幅值打印出来看一眼，如果出现1.1以上的离谱电压，多半是数据换算出了问题，及时止损比后面排查快得多。

3. 暂态稳定仿真实战：三相短路故障下的完整流程

3.1 故障序列设置：sw_con矩阵的行列约定与时序设计

PST的时域仿真通过开关控制矩阵sw_con来描述扰动事件，每一行定义一次操作，包括操作时间、操作类型和操作对象。以一次典型的三相短路为例：

% 故障序列定义 % 第1列：操作时间（秒） % 第2列：操作类型，2表示三相短路接地，1表示清除故障 % 第3列：短路母线编号或线路起始母线 % 第4列：线路终止母线（单端接地故障时为0） sw_con = [ 0.05 2 3 0; % 0.05s 时在母线3设置三相短路 0.15 1 3 0; % 0.15s 时清除短路故障 ]; n_sw = size(sw_con, 1); % 开关操作次数

故障时刻选在0.05秒，让系统先运行几个周波，把稳态波形充分稳定下来。清除时刻选在0.15秒，对应的短路持续时间为0.1秒，即6个周波（60Hz系统）。这个持续时间在绝大多数测试系统上都足以考察暂态稳定行为：太短看不出振荡过程，太长系统大概率直接失稳。

操作类型的具体数值在不同PST版本间有差异，有些版本用2表示三相短路，有些则用其他编码。建议拿到代码后先翻一下源码里对sw_con第2列的注释，或者直接跑一个官方示例对照结果。第3列和第4列描述操作对象时需要注意：如果是母线接地故障，第4列填0；如果是线路中间故障，需要在线路上增设一个虚拟母线，PST本身不会自动折算线路任意点的短路。

故障母线位置对暂态稳定性影响很大。通常选择距离发电机最近的母线做三相短路，因为短路点越接近机端，故障期间发电机输出的电磁功率越低，转子加速面积越大，系统越容易被推向失稳。实际工程中除了验证最严重工况，还应该扫描不同母线位置的故障，找出系统的薄弱环节。

3.2 时域仿真执行：步长选择与仿真时长

仿真步长的选择直接决定稳定性和精度。PST的积分器是显式数值积分，步长过大会出现数值振荡甚至发散，步长过小则仿真耗时成倍增加。针对60Hz系统，我一般把步长设在0.002到0.01秒之间，最常用的值是0.005秒。如果系统里有小时间常数环节（比如快速励磁系统），建议把步长压到0.001秒，否则那些快状态会制造数值刚度问题。

仿真时长的设定取决于你关心的动态过程。对于暂态稳定而言，首摆稳定性和多摆稳定性分别对应不同时间尺度：看首摆，3到5秒足够；看后续阻尼特性，至少需要10秒以上。PST的仿真函数内部使用固定步长，这一点和MATLAB自带的变步长ODE求解器不同。固定步长意味着无法自动加密时间点到满足精度，所以步长选择必须由使用者拍板。

关于仿真时长的设置，我通常先用5秒快速筛一遍临界切除时间，筛选通过后再把仿真时长拉到20秒做精细的阻尼分析。时长拉长后数据量会变大，尤其当步长为0.005秒时，20秒仿真产生4000个时间点，每个时间点存几十个状态变量，对MATLAB来说不是压力，但绘图时记得只画需要的时间段，否则整张图密集得看不出振荡的衰减特性。

% 设置仿真参数并执行时域仿真 t0 = 0; % 起始时间（秒） tf = 5; % 仿真时长（秒） h = 0.005; % 积分步长（秒） % 执行仿真，生成时间变量t和状态变量矩阵 simu(t0, tf, h, n_sw, sw_con);

仿真结束后，工作区会多出若干时间序列矩阵。最常见的三个是：t（时间向量）、mac_var（发电机状态变量）、母线电压幅值序列等。mac_var的行对应时刻，列对应各台发电机堆叠的状态量，前几列固定是转子角（弧度）和转速偏差（pu）。

这里有个容易栽跟头的地方：不同版本里mac_var的列顺序可能不同，别凭记忆直接取第2列当功角。我每次在新环境里跑完仿真后，都会打开PST源码里的状态量拼接部分，确认状态变量的排列方式再动手绘图。这个过程花不了两分钟，但能省掉后面半天的返工。

PST仿真本质上是联立求解微分代数方程组，微分方程来自发电机和控制器动态，代数方程来自网络潮流约束。显式积分在每一步长内对网络代数方程做迭代求解，因此步长设置不宜超过系统主导振荡周期的1/20。如果系统中有串联补偿或HVDC这类响应速度极快的器件，时间常数可能小于1毫秒，这时候固定步长的思路已经不适用，要么把步长压到0.0005秒以下，要么在模型层面做简化处理。

3.3 结果解读：功角曲线、转速轨迹与稳定判据

仿真完成后，第一件要做的事情是画出所有发电机的转子角相对参考机的差值。如果所有相对功角都在有限范围内振荡并且收敛，系统暂态稳定；只要有任意一台机的相对功角持续单调增大，这台机已经与系统失去同步，属于暂态失稳。

% 绘制功角曲线 figure; plot(t, mac_var(:, 1) * 180/pi, 'LineWidth', 1.5); hold on; plot(t, mac_var(:, 2) * 180/pi, 'LineWidth', 1.5); plot(t, mac_var(:, 3) * 180/pi, 'LineWidth', 1.5); xlabel('时间 (s)', 'FontSize', 12); ylabel('转子角 (deg)', 'FontSize', 12); legend('G1', 'G2', 'G3'); grid on;

注意这段代码假设mac_var的前三列就是三台发电机的转子角。如果你的系统有更多发电机，需要根据实际列数调整绘图范围。把弧度乘以180/pi转成度展示更直观。

从曲线形状看首摆幅度：如果两台机之间的相对功角最大摆幅超过180度，在物理上意味着转子已经越过不稳定平衡点，即使后续角度回落也大概率是数值假象，判据上直接视为失稳。更严谨的做法是用仿真末段的阻尼比来判断：对功角曲线后1秒的数据做衰减拟合，提取振荡频率和阻尼比，阻尼比低于0.02的模式需要警惕。

转速轨迹也是重要佐证。失稳的发电机转速偏差会持续攀升（pu值远超0.01），而稳定机组的转速偏差会围绕零值振荡衰减。把功角曲线和转速曲线放在同一张图里看，能够更快判断失稳类型是首摆失稳还是多摆振荡失稳。

这里补充一个概念：临界切除时间（CCT）是暂态稳定分析中一个很实用的指标。它表示系统能承受的最大故障持续时间，超过这个时间就会失稳。故障切除时间与CCT的差就是稳定裕度，工程上通常要求这个裕度在0.1秒以上，这部分在第5章会给出批量扫描的实现代码。

3.4 小信号分析：特征值计算与弱阻尼模式定位

时域仿真只能告诉你"稳不稳"，特征值分析可以进一步告诉你"稳定裕度有多少"。PST提供small_sig函数，基于当前稳态工作点线性化整个系统，输出状态矩阵A，然后直接用MATLAB的eig函数求特征值。小信号分析的价值在于，它能在不做时域仿真的前提下快速判断工作点的局部稳定性，还能直接提取主导振荡模式的频率和阻尼比。

% 小信号稳定分析 A = small_sig(bus, branch, gen, mac_con); % 计算特征值 eig_val = eig(A); % 计算振荡频率和阻尼比 freq = imag(eig_val) / (2 * pi); damp = -real(eig_val) ./ abs(eig_val); % 只保留振荡模式（虚部频率足够高的共轭复数） idx = find(freq > 0.1); damp_osc = damp(idx); freq_osc = freq(idx);

对于9节点3机系统，A的阶数通常在十几阶左右，特征值中包含实根（非振荡模式）和共轭复根（振荡模式）。筛选条件freq > 0.1是为了滤掉数值上接近零频率的实特征值，避免把非振荡模式当作弱阻尼模式处理。阻尼比排序的核心逻辑是负实部除以模值，越接近零说明振荡越难衰减。

small_sig返回的A矩阵阶数等于系统的动态状态量总数，状态量的排列顺序与mac_var保持一致。如果后续要计算参与因子，需要读一下PST源码里状态量的堆叠顺序，再做特征向量的参与因子计算，这一步算完基本能把弱阻尼模式定位到具体机组和具体环节。

特征值分析还有个实用场景：快速对比不同运行方式下的阻尼变化。比如要评估某种负荷模型的接入对系统阻尼的影响，只需要修改负荷模型参数后重新计算潮流和小信号，比较主导模式的阻尼比变化趋势，不跑大量时域仿真，计算速度比时域快一个数量级。

4. PST避坑指南：五个高频故障与逐条排查路径

4.1 潮流不收敛：靠调迭代次数没用

现象：pf函数迭代几十次后报错，输出NaN，母线电压变成离谱值。

原因：潮流不收敛的根源绝大多数在数据一致性。最典型的有三类：一是bus矩阵第5列有功负荷为负值且绝对值大于该母线发电机出力，导致净注入功率为负；二是branch矩阵的R和X直接用有名值填写而不是标幺值，潮流的数值范围完全不对；三是bus矩阵第2列的母线类型设置有误，比如系统里没有摆动母线，或者同时设置了两条摆动母线。

解决：先把所有母线的净注入功率算一遍，确认每条母线的有功无功平衡关系。再看branch矩阵里最大的X和最小的X是否在三个数量级以内。最后检查bus矩阵第2列，确保有且只有一条母线类型为1。我每次遇到不收敛都会先把这几项过一遍，通常问题就出在这几个地方。如果上面都查不出问题，再检查基准容量basmva，母线功率总量和支路阻抗的标幺值基准不一致是另一个隐藏根源。

4.2 时域仿真中途发散：步长和模型刚度是主因

现象：仿真在故障清除后突然数值爆炸，功角曲线直接飞掉，甚至出现NaN。

原因：时域仿真发散有两种典型情况。第一种是步长过大，显式积分器的稳定性受最大步长约束，超过临界步长后误差指数放大。第二种是系统里存在时间常数很小的状态变量，例如励磁系统的快速调节器和PSS的高频滤波环节，时间常数可能只有几毫秒，对步长提出了苛刻要求。

解决：先把步长从0.005缩小到0.001再跑一遍，如果曲线恢复正常，说明是步长问题。如果缩小步长仍然发散，就需要检查自定义模型里是否有不合理的参数配置，比如阻尼系数D被填成负值。还有一种常见做法是改用刚性求解器，但PST内置积分器不支持自动切换，我通常直接缩短步长，简单有效。

顺便说一个亲历案例：某次在模拟项目X里给励磁系统加了PSS，仿真跑到0.4秒后功角曲线开始高频抖动，看起来像噪声。一开始以为是数值问题，缩短步长后抖动还在，后来检查PSS参数才发现增益填成了负值，导致PSS输出和励磁电压形成正反馈。这类问题靠调积分步长解决不了，必须回到模型参数本身排查。

4.3 数据文件格式翻车：指数写法与空格陷阱

现象：数据文件运行不报错，但MATLAB变量里某些值明显不对，比如一个本应填0.1的数值变成了1e-10。

原因：PST数据文件大多从Excel或别的仿真工具导出文本，粘贴进.m文件时容易出现格式误差。最常见的是指数符号含混：原始数据用1E-4表示，MATLAB里没问题，但如果在E和数字之间多了空格，就成了1E-4被拆成两个token。另外负数的负号与前面的空格连在一起时，MATLAB可能把-4误解成单独表达式。

解决：数据文件加载后先执行whos查看各变量的尺寸和数值范围，再抽查几个关键母线的电压和支路功率是否落在合理范围。我养成了一个习惯：每次写完数据文件，先跑一遍disp(max(abs(bus(:, 3))))，如果最大电压超过1.2或小于0.8，几乎可以确定数据录入有问题。

4.4 MATLAB版本差异导致API不兼容

现象：同一份PST代码，在某台机器的MATLAB里运行正常，换到新版MATLAB就报未定义函数，或者函数名冲突。

原因：PST是多年积累的开源项目，部分函数名在新版MATLAB中被占用。最典型的例子是pf这个函数名，在MATLAB部分工具箱里存在同名函数，执行路径顺序不同时可能调用到错误版本。另外有些老代码用i、j作为循环变量，在新版MATLAB中可能触发复数虚部的隐式覆盖。

解决：在运行PST之前，用which pf确认执行的是PST目录下的函数。把PST的src目录通过addpath添加到MATLAB路径并置于路径列表顶部，避免同名冲突。循环变量统一改用ii和jj，别跟虚数单位较劲。

4.5 自定义模型状态量顺序错位

现象：接入自定义励磁模型后，状态变量矩阵的列数对不上，或者仿真结果里某台机的励磁电压波形明显异常。

原因：PST里自定义模型的状态量是按固定顺序堆叠到状态向量中的，新模型的状态量必须和现有状态量按约定顺序拼接。如果没有准确登记新状态量在整体状态向量中的起始列和结束列，结果矩阵和绘图函数都会拿错列、错数据。

解决：新增模型时先在源文件里找到状态量拼接区，把新状态量的维度写清楚，并且在初始化阶段打印一次状态变量的维度检查语句。仿真后再做一次状态变量的逐列对比，画出励磁电压、机械功率等变量，确认每一列的确是预期中的物理量。

5. 技巧：批量扫描临界切除时间验证暂态稳定裕度

在做暂态稳定分析时，临界切除时间（CCT）是衡量系统稳定裕度的核心指标。CCT的物理含义很直观：故障发生的那一刻起，如果故障在CCT以内被切除，系统能恢复到稳定运行；如果晚于CCT，系统就回不来了。PST的CCT扫描思路也很直接：把故障切除时间当成一个可调参数，不断增大它，观察系统是否失稳，用二分搜索找到稳定与失稳的临界点。

5.1 二分搜索CCT的实现

% 二分搜索临界切除时间 t_low = 0.05; % 已知稳定 t_high = 0.30; % 已知失稳 tol = 0.001; % 搜索精度 for iter = 1:20 t_mid = (t_low + t_high) / 2; % 更新故障清除时间 sw_con(2, 1) = t_mid; % 重新初始化工作点并仿真 [bus, branch, gen] = pf(bus, branch, gen, tol, itermax); mac_con = mac_ini(bus, gen); simu(t0, tf, h, n_sw, sw_con); % 判断是否稳定（相对功角是否超过180度） delta = mac_var(:, 2) - mac_var(:, 1); if max(abs(delta)) * 180/pi < 180 t_low = t_mid; % 系统稳定，向上界收缩 else t_high = t_mid; % 系统失稳，向下界收缩 end end cct = (t_low + t_high) / 2; fprintf('临界切除时间 CCT = %.3f s\n', cct);

这个二分搜索收敛速度很快，20次迭代以内就能达到毫秒级精度。搜索前需要确认t_low对应的工况确实是稳定的，t_high对应的工况确实是失稳的，否则区间不闭合，结果没有意义。每次修改清除时间后都要重新执行潮流和mac_ini，保证状态变量在每次仿真前都被重新初始化。

稳定判据从单台功角改为最大相对功角更可靠，因为大系统中不同机组之间可能存在多模式振荡，单台机的功角可能碰巧落在小摆幅区间，而其他机组早已失稳。用max(abs(delta))能捕获最危险的那台机。

5.2 从CCT到批量裕度分析

如果有多条线路、多个故障点需要逐一验证，把上面的二分搜索放进一个外层循环，遍历所有故障场景即可。每次进入新场景时要重建sw_con矩阵，故障对象换一次，母线编号跟着换。批量扫描里要控制时间成本：单次仿真时长缩短到3秒足够判断首摆失稳，因为绝大多数暂态失稳在3秒内已经露出端倪。多机系统里可能出现多摆失稳的情况，这时候再单独把仿真时长拉长重跑即可。

从那以后，我每次接入新系统或修改参数后，都会强制跑一遍CCT扫描作为常规回归测试。哪怕只是改了一个阻尼系数或者一条线路参数，CCT的变化也能非常敏感地暴露问题：比如CCT突然从0.2秒掉到0.08秒，多半是新参数引入了不该有的负阻尼，或者模型状态量顺序又错位了。这种习惯帮我提前拦下过不少只有在深夜跑批量时才爆发的低级错误，希望帮到你。

本文还有配套的精品资源，点击获取