做电力系统动态分析的同行应该都有过这种焦灼系统惯量这个指标拿不到后面一次调频、频率稳定性评估全都在踩软沙。传统做法无非两条路要么靠离线模型推导要么直接拿SCADA稳态数据硬凑可前者依赖模型算得准不准后者压根不知道惯量这种动态参数怎么从稳态数据里还原出来。实际上惯量是藏在“动态过程”里的想把它捞出来就得从扰动后的频率响应入手。我在这类项目里试过几种方案ARMAX模型是其中落地最快、效果也最稳的一个只需要一段带功率扰动的频率响应数据就能把系统等效惯量常数辨识出来核心函数不多MATLAB里几行就能跑通。这篇就把整个流程完整拆开附一份可直接复现的代码适合刚接触惯量评估、或者想换一条辨识路子的同学直接抄作业。先明确这篇笔记的定位。它不是教科书式的参数推导也不会堆一堆查不到背景的理论名词而是一条从真实需求出发的完整技术路线电力系统惯量到底是什么为什么ARMAX能把惯量捞出来数据到底要怎么预处理辨识结果才靠谱最后是一份我自己跑通的MATLAB代码以及几处调试时踩过的坑。文章里用到的数据、模型和参数全部可复现按步骤走一遍就能得到一张带噪声的仿真频率响应表以及一个和真实值非常接近的惯量辨识结果。这套方法做风机并网评估、孤网稳定性分析、微网惯量监测也同样适用不需要改动太多逻辑。1. 惯量评估在评估什么ARMAX为什么能上手1.1 惯量是藏在频率动态里的“惯性”电力系统频率动态的基础公式是摇摆方程形式不复杂但物理含义很重。在单机等值模型下发电机转子的动能会缓冲功率不平衡带来的转速变化写成标幺形式就是2H * (dΔf/dt) ΔP - D * Δf其中Δf是频率偏差标幺ΔP是机械功率减电磁功率的不平衡量标幺D是等效阻尼系数H就是我们要找的等效惯量常数单位是秒。H的物理意义可以理解成在额定功率下转子储存的动能还能维持发电机运行多少秒。常规火电、水电机组的H值通常在2到10秒之间新能源机组通过变流器接入后这一块等效惯量往往很小所以系统惯量水平逐年下降是个很现实的工程问题。关键就在扰动发生的瞬间比如某台机组突然跳闸ΔP瞬间阶跃上去但频率不可能瞬间掉下来。频率变化率RoCoF、功率缺额ΔP和惯量H之间存在一个近似关系dΔf/dt ≈ ΔP / (2H)。这个式子说明同样大小的扰动惯量越小频率跌得越快。所以惯量评估的本质就是从这一段频率动态过程里反推出H的大小。但难点在于测量到的频率信号是系统综合响应里面既有惯量信息也有阻尼信息还混着调速器动作和量测噪声。单纯取起点附近的斜率求RoCoF对噪声和扰动检测时延都很敏感离线建模虽然可以整定出H值但模型参数不确定性太大会让结果失去参考意义。这也是为什么数据驱动的辨识方法在惯量评估里越来越受关注。1.2 ARMAX模型在这个问题上的优势ARMAX全称是带外生输入的自回归滑动平均模型数学形式是A(q)y(t) B(q)u(t) C(q)e(t)。它做的事情很简单用过去的输出、过去的输入和过去噪声这三类信息去解释当前输出。在这里输入u是功率扰动输出y是频率偏差辨识出来的传递函数B(q)/A(q)就描述了功率扰动到频率响应之间的动态关系。ARMAX在处理惯量评估时有几个很实在的优势。第一它天然带噪声模型C(q)量测噪声和负荷随机波动都可以被它吸收掉这就比直接用最小二乘拟合一个传递函数要稳得多后者遇到有色噪声时参数估计会明显偏移。第二model函数可以直接给出A多项式和B多项式系数离散传递函数的极点和增益就藏在里面惯量时间常数可以精准提取。第三MATLAB系统辨识工具箱里的armax函数是成熟封装不需要自己写优化算法避免了一部分遗传算法或粒子群方案里反复调参、容易陷入局部最优的问题。我还对比过Prony分析、FFT谱估计这些方法。它们对稳态谐波和振荡模态分析有优势但做惯量评估时需要人为从频谱里挑主导模态操作起来繁琐而且低频段分辨率不够时非常容易出错。ARMAX是直接在时域里拟合带输入的动态过程思路更贴近物理模型工程上手也快。1.3 整体辨识思路整个辨识流程是线性的拆开看大概是五步。获取数据用仿真或实测记录一段“功率扰动频率响应”曲线注意扰动前后都要有足够的稳态段。数据预处理选择事件后数据段去基准偏差必要时重采样和低通滤波。模型辨识调用armax函数设置合适的模型阶次得到离散传递函数。参数提取从A多项式系数算离散极点从B多项式系数算静态增益再换算成等效时间常数和惯量。模型验证对比模型输出与实际数据、检查残差、观察不同阶次下参数是否稳定。这五步里最容易翻车的不是最后的公式代换而是第二步的数据预处理和第三步的阶次选择。很多同学上来就把原始信号直接喂给armax结果辨识出的A多项式完全不满足稳定条件或者惯量值飘得离谱然后把问题归咎于模型不行。其实绝大部分情况是数据没弄干净阶次也没想清楚。下面我把每一步的原理和代码都展开讲。2. 模型原理与关键细节2.1 离散时间常数与惯量的数学关系先做一个很标准的推导把连续域的惯性模型转换到离散域再看ARMAX系数和惯量之间怎么对应。对摇摆方程做拉普拉斯变换在无调速器参与、只保留惯量和阻尼的情况下G(s) Δf(s) / ΔP(s) 1 / (2Hs D)这是一个标准一阶惯性系统静态增益K 1/D时间常数τ 2H/D。用零阶保持器进行离散化采样时间Ts下这个连续系统的离散传递函数是G(z) K * (1 - α) / (z - α)其中 α exp(-Ts / τ)再看ARMAX辨识出来的离散传递函数。对于na1、nb1、nk1的一阶模型A(q) 1 a1 * q^-1B(q) b1 * q^-1传递函数就是G(q) b1 * q^-1 / (1 a1 * q^-1) b1 / (z a1)对比上面零阶保持离散化的形式就能得到两组关键关系α -a1也就是说 exp(-Ts/τ) -a1b1 K * (1 - α)静态增益K b1 / (1 - α)知道了α就可以反解时间常数τ -Ts / ln(α)同时K 1/D所以D 1/K。再结合τ 2H/D最终得到H τ * D / 2 τ / (2K)把K用b1和α带进去也可以直接写成H [-Ts / ln(-a1)] * (1 a1) / (2 * b1)这里有个细节很关键离散极点α通常非常接近1因为电力系统惯性时间常数有数秒而采样周期只有几十毫秒所以a1会非常接近-1。这种数值形态对噪声很敏感稍微拟合偏一点ln(-a1)就会明显变化导致H估计偏差。因此数据质量、模型阶次、辨识算法稳定性必须同时可控。下面这个表把符号关系整理了一遍方便对照写代码物理量连续域表示离散域ARMAX表示等效时间常数 τ2H/D-Ts / ln(-a1)静态增益 K1/Db1 / (1 a1)阻尼系数 D原值(1 a1) / b1惯量常数 HDτ/2[-Ts/ln(-a1)] * (1a1) / (2b1)2.2 噪声模型C(q)为什么重要ARMAX和ARX模型的区别就在C(q)。ARX相当于C(q)固定为1假设噪声是白噪声直接叠加在输出上。实际电力系统的频率量测噪声并不是纯白噪声它包含PMU量化误差、相位延迟、通信丢包带来的毛刺甚至还有负荷随机波动的低频成分。如果忽略这些噪声结构A多项式和B多项式的估计就会产生偏差惯量算出来自然不可信。ARMAX把噪声建模成MA过程等于给辨识过程加了一个“噪声吸震器”。在频响数据信噪比不高的情况下这个模型结构能明显提升参数估计的稳定性。我在仿真里对比过信噪比大约30dB时ARX辨识H的偏差可能在5%以上而ARMAX能做到1%左右。实际工程信号质量往往比实验室仿真差所以直接用ARMAX是更稳妥的选择。2.3 数据预处理决定成败的隐蔽层数据预处理是整套流程里最容易被忽视、又最能影响结果的一环。我见过不少人把半个小时长的频率曲线直接丢给armax结果A多项式拟合出来完全不满足一阶系统假设这是因为长期趋势、调频动作、相邻事件扰动都会干扰参数提取。第一步是截断数据段。一般取扰动前12秒和扰动后1030秒这段窗口既能覆盖主要动态过程又不会掺入太多无关低频分量。窗口太长系统频率回升段会引入调速器模型的影响让等效单机模型失真窗口太短频率还没进入稳态初始斜率信息不够τ辨识不出来。第二步是去基准。把扰动前的平均频率偏差作为基准值从整段信号里减掉这样可以消除稳态频率偏移的影响让辨识模型的稳态初始条件归零。这个操作对应到实际PMU数据里就是因为负荷水平不同导致的频率基线漂移。第三步是采样和滤波。如果原始数据采样率太高比如几百赫兹建议降采样到2050HzMatlab的resample函数可以直接处理。降采样前要过一遍低通滤波器截止频率取510Hz足够了这样既保留机电动态信息又滤掉高频噪声和可能的数值尖峰。去掉坏点是最后一步特别是PMU数据里的跳变毛刺要用中值滤波或人工检查的方式剔除否则一个异常点就足以把A多项式参数拉偏。3. 完整MATLAB代码实现3.1 第1段生成仿真数据搭一个可控的测试环境先搭一个可控的仿真数据环境。我们用一阶惯性系统模拟真实系统频率响应然后叠加量测噪声这样既知道真实H值又能检验后面辨识算法的精度。建议你也按这个流程做先把算法在仿真环境里打通再上实测数据否则问题定位会很痛苦。%% 1. 参数设置与仿真数据生成 clear; clc; rng(2026); H_true 5.0; % 真实惯量常数单位秒 D_true 2.0; % 等效阻尼系数标幺值 Ts 0.05; % 采样周期单位秒 Tfinal 30; % 仿真总时长单位秒 t (0:Ts:Tfinal); N length(t); % 功率扰动1% 阶跃模拟机组跳闸或负荷突增 u zeros(N, 1); u(t 2) 0.01; % 连续一阶系统G(s) (1/D) / ((2H/D)*s 1) sys tf(1/D_true, [2*H_true/D_true, 1]); % 理想频率响应 y_ideal lsim(sys, u, t); % 叠加测量噪声模拟PMU量测误差标准差0.0002 p.u.约为0.01Hz noise 0.0002 * randn(N, 1); y_meas y_ideal noise;这里H_true取5秒相当于一台常规火电机组或一组机组等值的惯量水平D_true取2.0是为了让系统的等效阻尼更接近工程实际也让稳态频率偏差不至于过大。1%的功率阶跃在工程上属于中等偏小的扰动在这种扰动幅度下把H辨识准方法才更有说服力。采样周期选0.05秒也就是20Hz跟多数PMU的50帧/秒相比略低但足够覆盖机电暂态的低频段而且能减轻后面辨识算法对高频噪声的敏感度。运行完这段y_meas就是带噪声的“实测频率偏差信号”它的物理单位是标幺值0.005对应的实际频率偏差是0.25Hz以50Hz为基准。仿真里真实惯量已知后面就知道辨识结果准不准。3.2 第2段事件段提取与数据预处理数据预处理我用一个独立小节来写因为这是我最想强调的部分。下面这段代码把仿真信号转成辨识用的数据对象。%% 2. 数据预处理 % 只保留扰动发生之后的数据段 ind t 2; t_seg t(ind); u_seg u(ind); y_seg y_meas(ind); % 计算扰动前基准值去基线 base_idx t 2; u_base mean(u(base_idx)); y_base mean(y_meas(base_idx)); u_seg u_seg - u_base; y_seg y_seg - y_base; % 构建系统辨识工具箱的数据对象 data iddata(y_seg, u_seg, Ts); data.Name Frequency response estimation data;注意去基准这段代码我用的是扰动前信号的平均值。在仿真数据里扰动前u和y都是0所以这段操作看起来像是“白做”但它对实测数据非常重要。实测PMU频率在扰动前往往在49.98Hz、50.02Hz附近波动如果不把这个偏移减掉等于给模型强加了一个非零的初始条件A多项式的辨识会明显受影响。到这一步已经可以直接调用armax了但我建议你先用iddata的plot命令快速看一眼数据形态确认事件检测准确、没有明显毛刺再往下走。数据问题是辨识问题里最磨人的一类越早发现越好。3.3 第3段ARMAX辨识与惯量提取核心代码来了。这一段完成模型辨识并从A、B多项式系数中提取惯量。%% 3. ARMAX模型辨识 na 1; nb 1; nc 1; nk 1; model armax(data, [na nb nc nk]); % 查看模型结构 disp(model); % 获取A、B多项式系数 A_coef model.A; % A(q) 1 a1*q^-1 B_coef model.B; % B(q) 0 b1*q^-1 或 b1*q^-1 a1 A_coef(2); b1 B_coef(end); % 取最后一个非零系数 %% 4. 从ARMAX系数提取惯量 alpha -a1; % 离散极点 alpha exp(-Ts/tau) if alpha 0 || alpha 1 error(辨识出的极点不在有效范围(0,1)内请检查数据或模型阶次。); end tau -Ts / log(alpha); % 等效时间常数 tau 2H/D K b1 / (1 - alpha); % 静态增益 K 1/D D_est 1 / K; H_est tau / (2 * K); % H D*tau/2 tau/(2K) fprintf(真实值: H %.3f s, D %.3f\n, H_true, D_true); fprintf(辨识值: H %.3f s, D %.3f\n, H_est, D_est); fprintf(等效时间常数 tau %.3f s\n, tau);这里我习惯写成nal、nb、nc、nk四个参数其中nk表示输入延迟的采样点数。对电力系统功率扰动到频率响应的物理过程来说延迟非常小一般取1就够也就是从当前采样时刻的输入开始影响到下一个采样时刻的输出。如果你手动设置更大的nk等于人为引入了一个滞后会浪费一部分动态信息不建议这么做。运行这段代码后你会看到类似下面的输出由于噪声是随机的每次结果会有细微差异真实值: H 5.000 s, D 2.000 辨识值: H 4.968 s, D 2.013 等效时间常数 tau 4.934 s辨识值和真实值相差不到1%说明整个链路是通的。tau对应2H/D真实值是5秒辨识值是4.93秒也符合预期。Maltab的armax函数返回的是idpoly对象model.A是多项式系数行向量从常数项开始排列所以A_coef(2)就是q^-1的系数a1。model.B的系数里可能会有前导零表示输入延迟所以取B_coef(end)是最稳的。如果你打印model对象看到的B(q)可能是“0.002 q^-1”这种形式对应B_coef就是[0, 0.002]。3.4 第4段结果验证与可视化参数算出来只是第一步还要验证模型靠不靠谱。这段代码把ARMAX模型的阶跃响应和原始测量数据叠加对比同时画出残差自相关快速判断拟合质量。%% 5. 模型验证与可视化 % 对比模型输出与实际数据 y_hat compare(data, model); figure; plot(t_seg, y_seg, k, LineWidth, 1.0); hold on; plot(t_seg, y_hat.OutputData, r--, LineWidth, 1.5); grid on; xlabel(时间 / s); ylabel(频率偏差 / p.u.); legend(测量数据, ARMAX模型输出, Location, best); title(ARMAX模型拟合效果对比); % 残差分析 figure; resid(model, data);compare函数会返回模型在同样输入下的仿真输出如果曲线和测量数据贴合良好说明辨识模型抓住了主导动态。resid命令会画出残差的自相关函数和互相关函数理想的残差应该接近白噪声也就是说自相关曲线应该在置信区间内小幅波动没有明显的周期成分。如果残差里还有明显相关结构说明模型阶次不够或者数据里还有没滤掉的干扰源。这里再补充一种更直观的验证方式把辨识出的离散模型转换成连续域模型然后画Bode图看低频增益和穿越频率是否与理论值吻合。代码就一行sys_est_id tf(model); bode(sys_est_id, sys);把辨识模型和真实连续模型画在同一张Bode图上如果两条曲线在机电动态频段0.011Hz基本重合那说明辨识系统不仅数值对动态特性也对。4. 常见问题与排查技巧4.1 辨识出的极点不在(0,1)范围怎么办我在实际调试里遇到过好几次alpha算出来是负值或者大于1的情况这时程序会直接报错。这种问题通常有三个来源。第一数据段取得太短频率动态还没走完A多项式拟合不到正确的极点位置。这种情况把时间窗口加长到10秒以上一般能解决。第二测量噪声太大一阶模型被噪声主导。此时需要加强低通滤波或者把采样率降下来减少高频噪声对极点估计的干扰。第三事件并非阶跃扰动比如是缓慢爬坡的功率变化这种输入激励对一阶系统的极点辨识激励度不够模型参数容易跑飞。还有一类特殊情况系统本身不是严格一阶可能是多机系统或含有调速器动态这时强迫用na1去拟合得到的a1是多种动态模式的折中结果不一定满足稳定域约束。解决办法是提高模型阶次让ARMAX去拟合更高阶动态然后从A多项式的根里挑主导极点来算惯量。4.2 模型阶次na、nb、nc到底怎么选模型阶次选择没有银弹但我给你一个足够稳妥的启动方案先从na1、nb1、nc1开始跑通后再逐步升阶。伺服看两个指标一个是残差是否接近白噪声另一个是辨识出的H在不同阶次下是否稳定在一个小范围内。下面是我用同一组仿真数据跑出来的结果对比你可以感受一下阶次的影响阶次组合 (na, nb, nc)辨识H (s)辨识D残差状态(1, 1, 1)4.972.01良好(2, 2, 2)5.032.02良好(3, 3, 3)5.112.05轻微过拟合(1, 1, 0) ARX5.211.94残差有相关性可以看到一阶ARMAX已经能给出不错的结果二阶模型也不会差太多但三阶以上反而开始过拟合H估计偏离真实值更多。ARX模型由于没有噪声模型残差里残留相关结构参数偏得也更快。所以我的经验是一阶模型作为基准二阶做交叉验证如果两个结果差异在5%以内那组参数基本可信。4.3 实测数据的几个坑实测数据和仿真数据之间的差距比很多初学者想象中大得多。第一个坑是触发时刻不精确。实测数据要确定功率扰动发生的准确时刻通常用频率变化率或者功率突变检测但检测算法本身有延迟数据段的起始点因此可能偏移几十毫秒。ARMAX对起点偏移有一定鲁棒性因为A多项式拟合的是整个数据段的动态但如果起点偏差超过一个时间常数结果就会明显变差。第二个坑是频率信号的质量。PMU的频率通道有内部滤波和相位延迟不同厂家的PMU延迟特性还不一样这相当于在输入输出通路里加了一个未知延迟。处理手法是把nk从小调到大比如从1测试到5观察H估计值是否稳定。如果nk增大后结果变化很大说明原始信号延迟问题比较严重可以考虑先对齐事件触发点。第三个坑是数据段内包含多次扰动。如果评估窗口内还有其他机组跳闸或者负荷突变ARMAX会把多次扰动当做一个叠加输入辨识结果必然失真。解决方法是先做事件检测把每个扰动事件单独切段每次评估一个事件然后把多次事件的H估计结果取平均或中位数。第四个坑是最隐蔽的多机系统的惯量能否用一个等效H来概括。ARMAX辨识出来的本质上是“从测量点到系统之间的等效动态”它反映的是这个测量点能观察到的系统惯量水平而不是全部机组惯量的简单求和。做系统级惯量监测时这个值很有意义但不要试图用单点测量去还原每台机的惯量分布。4.4 独家避坑技巧这里分享几个常规文档里不常写、但实操时非常管用的经验。做参数提取时我强烈建议不要把H和D分开单独看。因为从式子里看H和D是耦合的噪声扰动可能让H偏低、D偏高但两者的乘积对应的时间常数却相对稳定。如果多次辨识的τ值都稳定在某个区间而H和D波动较大那大概率是数据信噪比不足而不是模型错了这时候参考τ比参考H更有意义。数据段长度的选择有个经验规律扰动后窗口建议覆盖至少3倍时间常数。真实时间常数未知怎么办先粗拟合一次估计出τ再回头调整窗口长度重新辨识一遍。这种“迭代式窗口选择”能让结果稳定不少。另外用多次蒙特卡洛仿真来校验算法的精度是个好习惯。同一个真值系统换不同随机种子生成十组噪声数据辨识十次取均值和标准差。这样你能直观看到算法的重复性。如果标准差偏大说明该工况下辨识条件不理想直接看单次结果没有意义。最后再分享一个我觉得很实用的扩展思路。这套方法不只适用于单机等值系统在微网、风电场并网点、孤岛系统这些场景里只要你能获取并网点的功率扰动数据和频率响应数据就能用同样的流程辨识等效惯量。我做风电场并网评估时把这个流程封装成了一个脚本输入是CSV格式的功率和频率信号输出就是惯量H、阻尼D和置信区间。工程上拿来做趋势监测完全够用。希望这份带完整代码的笔记能帮你少走几个弯路。后面你如果遇到辨识结果不稳或者数据预处理上的问题欢迎在评论区把情况贴出来我们一起讨论。
用ARMAX模型辨识电力系统惯量：原理与MATLAB实现
做电力系统动态分析的同行应该都有过这种焦灼系统惯量这个指标拿不到后面一次调频、频率稳定性评估全都在踩软沙。传统做法无非两条路要么靠离线模型推导要么直接拿SCADA稳态数据硬凑可前者依赖模型算得准不准后者压根不知道惯量这种动态参数怎么从稳态数据里还原出来。实际上惯量是藏在“动态过程”里的想把它捞出来就得从扰动后的频率响应入手。我在这类项目里试过几种方案ARMAX模型是其中落地最快、效果也最稳的一个只需要一段带功率扰动的频率响应数据就能把系统等效惯量常数辨识出来核心函数不多MATLAB里几行就能跑通。这篇就把整个流程完整拆开附一份可直接复现的代码适合刚接触惯量评估、或者想换一条辨识路子的同学直接抄作业。先明确这篇笔记的定位。它不是教科书式的参数推导也不会堆一堆查不到背景的理论名词而是一条从真实需求出发的完整技术路线电力系统惯量到底是什么为什么ARMAX能把惯量捞出来数据到底要怎么预处理辨识结果才靠谱最后是一份我自己跑通的MATLAB代码以及几处调试时踩过的坑。文章里用到的数据、模型和参数全部可复现按步骤走一遍就能得到一张带噪声的仿真频率响应表以及一个和真实值非常接近的惯量辨识结果。这套方法做风机并网评估、孤网稳定性分析、微网惯量监测也同样适用不需要改动太多逻辑。1. 惯量评估在评估什么ARMAX为什么能上手1.1 惯量是藏在频率动态里的“惯性”电力系统频率动态的基础公式是摇摆方程形式不复杂但物理含义很重。在单机等值模型下发电机转子的动能会缓冲功率不平衡带来的转速变化写成标幺形式就是2H * (dΔf/dt) ΔP - D * Δf其中Δf是频率偏差标幺ΔP是机械功率减电磁功率的不平衡量标幺D是等效阻尼系数H就是我们要找的等效惯量常数单位是秒。H的物理意义可以理解成在额定功率下转子储存的动能还能维持发电机运行多少秒。常规火电、水电机组的H值通常在2到10秒之间新能源机组通过变流器接入后这一块等效惯量往往很小所以系统惯量水平逐年下降是个很现实的工程问题。关键就在扰动发生的瞬间比如某台机组突然跳闸ΔP瞬间阶跃上去但频率不可能瞬间掉下来。频率变化率RoCoF、功率缺额ΔP和惯量H之间存在一个近似关系dΔf/dt ≈ ΔP / (2H)。这个式子说明同样大小的扰动惯量越小频率跌得越快。所以惯量评估的本质就是从这一段频率动态过程里反推出H的大小。但难点在于测量到的频率信号是系统综合响应里面既有惯量信息也有阻尼信息还混着调速器动作和量测噪声。单纯取起点附近的斜率求RoCoF对噪声和扰动检测时延都很敏感离线建模虽然可以整定出H值但模型参数不确定性太大会让结果失去参考意义。这也是为什么数据驱动的辨识方法在惯量评估里越来越受关注。1.2 ARMAX模型在这个问题上的优势ARMAX全称是带外生输入的自回归滑动平均模型数学形式是A(q)y(t) B(q)u(t) C(q)e(t)。它做的事情很简单用过去的输出、过去的输入和过去噪声这三类信息去解释当前输出。在这里输入u是功率扰动输出y是频率偏差辨识出来的传递函数B(q)/A(q)就描述了功率扰动到频率响应之间的动态关系。ARMAX在处理惯量评估时有几个很实在的优势。第一它天然带噪声模型C(q)量测噪声和负荷随机波动都可以被它吸收掉这就比直接用最小二乘拟合一个传递函数要稳得多后者遇到有色噪声时参数估计会明显偏移。第二model函数可以直接给出A多项式和B多项式系数离散传递函数的极点和增益就藏在里面惯量时间常数可以精准提取。第三MATLAB系统辨识工具箱里的armax函数是成熟封装不需要自己写优化算法避免了一部分遗传算法或粒子群方案里反复调参、容易陷入局部最优的问题。我还对比过Prony分析、FFT谱估计这些方法。它们对稳态谐波和振荡模态分析有优势但做惯量评估时需要人为从频谱里挑主导模态操作起来繁琐而且低频段分辨率不够时非常容易出错。ARMAX是直接在时域里拟合带输入的动态过程思路更贴近物理模型工程上手也快。1.3 整体辨识思路整个辨识流程是线性的拆开看大概是五步。获取数据用仿真或实测记录一段“功率扰动频率响应”曲线注意扰动前后都要有足够的稳态段。数据预处理选择事件后数据段去基准偏差必要时重采样和低通滤波。模型辨识调用armax函数设置合适的模型阶次得到离散传递函数。参数提取从A多项式系数算离散极点从B多项式系数算静态增益再换算成等效时间常数和惯量。模型验证对比模型输出与实际数据、检查残差、观察不同阶次下参数是否稳定。这五步里最容易翻车的不是最后的公式代换而是第二步的数据预处理和第三步的阶次选择。很多同学上来就把原始信号直接喂给armax结果辨识出的A多项式完全不满足稳定条件或者惯量值飘得离谱然后把问题归咎于模型不行。其实绝大部分情况是数据没弄干净阶次也没想清楚。下面我把每一步的原理和代码都展开讲。2. 模型原理与关键细节2.1 离散时间常数与惯量的数学关系先做一个很标准的推导把连续域的惯性模型转换到离散域再看ARMAX系数和惯量之间怎么对应。对摇摆方程做拉普拉斯变换在无调速器参与、只保留惯量和阻尼的情况下G(s) Δf(s) / ΔP(s) 1 / (2Hs D)这是一个标准一阶惯性系统静态增益K 1/D时间常数τ 2H/D。用零阶保持器进行离散化采样时间Ts下这个连续系统的离散传递函数是G(z) K * (1 - α) / (z - α)其中 α exp(-Ts / τ)再看ARMAX辨识出来的离散传递函数。对于na1、nb1、nk1的一阶模型A(q) 1 a1 * q^-1B(q) b1 * q^-1传递函数就是G(q) b1 * q^-1 / (1 a1 * q^-1) b1 / (z a1)对比上面零阶保持离散化的形式就能得到两组关键关系α -a1也就是说 exp(-Ts/τ) -a1b1 K * (1 - α)静态增益K b1 / (1 - α)知道了α就可以反解时间常数τ -Ts / ln(α)同时K 1/D所以D 1/K。再结合τ 2H/D最终得到H τ * D / 2 τ / (2K)把K用b1和α带进去也可以直接写成H [-Ts / ln(-a1)] * (1 a1) / (2 * b1)这里有个细节很关键离散极点α通常非常接近1因为电力系统惯性时间常数有数秒而采样周期只有几十毫秒所以a1会非常接近-1。这种数值形态对噪声很敏感稍微拟合偏一点ln(-a1)就会明显变化导致H估计偏差。因此数据质量、模型阶次、辨识算法稳定性必须同时可控。下面这个表把符号关系整理了一遍方便对照写代码物理量连续域表示离散域ARMAX表示等效时间常数 τ2H/D-Ts / ln(-a1)静态增益 K1/Db1 / (1 a1)阻尼系数 D原值(1 a1) / b1惯量常数 HDτ/2[-Ts/ln(-a1)] * (1a1) / (2b1)2.2 噪声模型C(q)为什么重要ARMAX和ARX模型的区别就在C(q)。ARX相当于C(q)固定为1假设噪声是白噪声直接叠加在输出上。实际电力系统的频率量测噪声并不是纯白噪声它包含PMU量化误差、相位延迟、通信丢包带来的毛刺甚至还有负荷随机波动的低频成分。如果忽略这些噪声结构A多项式和B多项式的估计就会产生偏差惯量算出来自然不可信。ARMAX把噪声建模成MA过程等于给辨识过程加了一个“噪声吸震器”。在频响数据信噪比不高的情况下这个模型结构能明显提升参数估计的稳定性。我在仿真里对比过信噪比大约30dB时ARX辨识H的偏差可能在5%以上而ARMAX能做到1%左右。实际工程信号质量往往比实验室仿真差所以直接用ARMAX是更稳妥的选择。2.3 数据预处理决定成败的隐蔽层数据预处理是整套流程里最容易被忽视、又最能影响结果的一环。我见过不少人把半个小时长的频率曲线直接丢给armax结果A多项式拟合出来完全不满足一阶系统假设这是因为长期趋势、调频动作、相邻事件扰动都会干扰参数提取。第一步是截断数据段。一般取扰动前12秒和扰动后1030秒这段窗口既能覆盖主要动态过程又不会掺入太多无关低频分量。窗口太长系统频率回升段会引入调速器模型的影响让等效单机模型失真窗口太短频率还没进入稳态初始斜率信息不够τ辨识不出来。第二步是去基准。把扰动前的平均频率偏差作为基准值从整段信号里减掉这样可以消除稳态频率偏移的影响让辨识模型的稳态初始条件归零。这个操作对应到实际PMU数据里就是因为负荷水平不同导致的频率基线漂移。第三步是采样和滤波。如果原始数据采样率太高比如几百赫兹建议降采样到2050HzMatlab的resample函数可以直接处理。降采样前要过一遍低通滤波器截止频率取510Hz足够了这样既保留机电动态信息又滤掉高频噪声和可能的数值尖峰。去掉坏点是最后一步特别是PMU数据里的跳变毛刺要用中值滤波或人工检查的方式剔除否则一个异常点就足以把A多项式参数拉偏。3. 完整MATLAB代码实现3.1 第1段生成仿真数据搭一个可控的测试环境先搭一个可控的仿真数据环境。我们用一阶惯性系统模拟真实系统频率响应然后叠加量测噪声这样既知道真实H值又能检验后面辨识算法的精度。建议你也按这个流程做先把算法在仿真环境里打通再上实测数据否则问题定位会很痛苦。%% 1. 参数设置与仿真数据生成 clear; clc; rng(2026); H_true 5.0; % 真实惯量常数单位秒 D_true 2.0; % 等效阻尼系数标幺值 Ts 0.05; % 采样周期单位秒 Tfinal 30; % 仿真总时长单位秒 t (0:Ts:Tfinal); N length(t); % 功率扰动1% 阶跃模拟机组跳闸或负荷突增 u zeros(N, 1); u(t 2) 0.01; % 连续一阶系统G(s) (1/D) / ((2H/D)*s 1) sys tf(1/D_true, [2*H_true/D_true, 1]); % 理想频率响应 y_ideal lsim(sys, u, t); % 叠加测量噪声模拟PMU量测误差标准差0.0002 p.u.约为0.01Hz noise 0.0002 * randn(N, 1); y_meas y_ideal noise;这里H_true取5秒相当于一台常规火电机组或一组机组等值的惯量水平D_true取2.0是为了让系统的等效阻尼更接近工程实际也让稳态频率偏差不至于过大。1%的功率阶跃在工程上属于中等偏小的扰动在这种扰动幅度下把H辨识准方法才更有说服力。采样周期选0.05秒也就是20Hz跟多数PMU的50帧/秒相比略低但足够覆盖机电暂态的低频段而且能减轻后面辨识算法对高频噪声的敏感度。运行完这段y_meas就是带噪声的“实测频率偏差信号”它的物理单位是标幺值0.005对应的实际频率偏差是0.25Hz以50Hz为基准。仿真里真实惯量已知后面就知道辨识结果准不准。3.2 第2段事件段提取与数据预处理数据预处理我用一个独立小节来写因为这是我最想强调的部分。下面这段代码把仿真信号转成辨识用的数据对象。%% 2. 数据预处理 % 只保留扰动发生之后的数据段 ind t 2; t_seg t(ind); u_seg u(ind); y_seg y_meas(ind); % 计算扰动前基准值去基线 base_idx t 2; u_base mean(u(base_idx)); y_base mean(y_meas(base_idx)); u_seg u_seg - u_base; y_seg y_seg - y_base; % 构建系统辨识工具箱的数据对象 data iddata(y_seg, u_seg, Ts); data.Name Frequency response estimation data;注意去基准这段代码我用的是扰动前信号的平均值。在仿真数据里扰动前u和y都是0所以这段操作看起来像是“白做”但它对实测数据非常重要。实测PMU频率在扰动前往往在49.98Hz、50.02Hz附近波动如果不把这个偏移减掉等于给模型强加了一个非零的初始条件A多项式的辨识会明显受影响。到这一步已经可以直接调用armax了但我建议你先用iddata的plot命令快速看一眼数据形态确认事件检测准确、没有明显毛刺再往下走。数据问题是辨识问题里最磨人的一类越早发现越好。3.3 第3段ARMAX辨识与惯量提取核心代码来了。这一段完成模型辨识并从A、B多项式系数中提取惯量。%% 3. ARMAX模型辨识 na 1; nb 1; nc 1; nk 1; model armax(data, [na nb nc nk]); % 查看模型结构 disp(model); % 获取A、B多项式系数 A_coef model.A; % A(q) 1 a1*q^-1 B_coef model.B; % B(q) 0 b1*q^-1 或 b1*q^-1 a1 A_coef(2); b1 B_coef(end); % 取最后一个非零系数 %% 4. 从ARMAX系数提取惯量 alpha -a1; % 离散极点 alpha exp(-Ts/tau) if alpha 0 || alpha 1 error(辨识出的极点不在有效范围(0,1)内请检查数据或模型阶次。); end tau -Ts / log(alpha); % 等效时间常数 tau 2H/D K b1 / (1 - alpha); % 静态增益 K 1/D D_est 1 / K; H_est tau / (2 * K); % H D*tau/2 tau/(2K) fprintf(真实值: H %.3f s, D %.3f\n, H_true, D_true); fprintf(辨识值: H %.3f s, D %.3f\n, H_est, D_est); fprintf(等效时间常数 tau %.3f s\n, tau);这里我习惯写成nal、nb、nc、nk四个参数其中nk表示输入延迟的采样点数。对电力系统功率扰动到频率响应的物理过程来说延迟非常小一般取1就够也就是从当前采样时刻的输入开始影响到下一个采样时刻的输出。如果你手动设置更大的nk等于人为引入了一个滞后会浪费一部分动态信息不建议这么做。运行这段代码后你会看到类似下面的输出由于噪声是随机的每次结果会有细微差异真实值: H 5.000 s, D 2.000 辨识值: H 4.968 s, D 2.013 等效时间常数 tau 4.934 s辨识值和真实值相差不到1%说明整个链路是通的。tau对应2H/D真实值是5秒辨识值是4.93秒也符合预期。Maltab的armax函数返回的是idpoly对象model.A是多项式系数行向量从常数项开始排列所以A_coef(2)就是q^-1的系数a1。model.B的系数里可能会有前导零表示输入延迟所以取B_coef(end)是最稳的。如果你打印model对象看到的B(q)可能是“0.002 q^-1”这种形式对应B_coef就是[0, 0.002]。3.4 第4段结果验证与可视化参数算出来只是第一步还要验证模型靠不靠谱。这段代码把ARMAX模型的阶跃响应和原始测量数据叠加对比同时画出残差自相关快速判断拟合质量。%% 5. 模型验证与可视化 % 对比模型输出与实际数据 y_hat compare(data, model); figure; plot(t_seg, y_seg, k, LineWidth, 1.0); hold on; plot(t_seg, y_hat.OutputData, r--, LineWidth, 1.5); grid on; xlabel(时间 / s); ylabel(频率偏差 / p.u.); legend(测量数据, ARMAX模型输出, Location, best); title(ARMAX模型拟合效果对比); % 残差分析 figure; resid(model, data);compare函数会返回模型在同样输入下的仿真输出如果曲线和测量数据贴合良好说明辨识模型抓住了主导动态。resid命令会画出残差的自相关函数和互相关函数理想的残差应该接近白噪声也就是说自相关曲线应该在置信区间内小幅波动没有明显的周期成分。如果残差里还有明显相关结构说明模型阶次不够或者数据里还有没滤掉的干扰源。这里再补充一种更直观的验证方式把辨识出的离散模型转换成连续域模型然后画Bode图看低频增益和穿越频率是否与理论值吻合。代码就一行sys_est_id tf(model); bode(sys_est_id, sys);把辨识模型和真实连续模型画在同一张Bode图上如果两条曲线在机电动态频段0.011Hz基本重合那说明辨识系统不仅数值对动态特性也对。4. 常见问题与排查技巧4.1 辨识出的极点不在(0,1)范围怎么办我在实际调试里遇到过好几次alpha算出来是负值或者大于1的情况这时程序会直接报错。这种问题通常有三个来源。第一数据段取得太短频率动态还没走完A多项式拟合不到正确的极点位置。这种情况把时间窗口加长到10秒以上一般能解决。第二测量噪声太大一阶模型被噪声主导。此时需要加强低通滤波或者把采样率降下来减少高频噪声对极点估计的干扰。第三事件并非阶跃扰动比如是缓慢爬坡的功率变化这种输入激励对一阶系统的极点辨识激励度不够模型参数容易跑飞。还有一类特殊情况系统本身不是严格一阶可能是多机系统或含有调速器动态这时强迫用na1去拟合得到的a1是多种动态模式的折中结果不一定满足稳定域约束。解决办法是提高模型阶次让ARMAX去拟合更高阶动态然后从A多项式的根里挑主导极点来算惯量。4.2 模型阶次na、nb、nc到底怎么选模型阶次选择没有银弹但我给你一个足够稳妥的启动方案先从na1、nb1、nc1开始跑通后再逐步升阶。伺服看两个指标一个是残差是否接近白噪声另一个是辨识出的H在不同阶次下是否稳定在一个小范围内。下面是我用同一组仿真数据跑出来的结果对比你可以感受一下阶次的影响阶次组合 (na, nb, nc)辨识H (s)辨识D残差状态(1, 1, 1)4.972.01良好(2, 2, 2)5.032.02良好(3, 3, 3)5.112.05轻微过拟合(1, 1, 0) ARX5.211.94残差有相关性可以看到一阶ARMAX已经能给出不错的结果二阶模型也不会差太多但三阶以上反而开始过拟合H估计偏离真实值更多。ARX模型由于没有噪声模型残差里残留相关结构参数偏得也更快。所以我的经验是一阶模型作为基准二阶做交叉验证如果两个结果差异在5%以内那组参数基本可信。4.3 实测数据的几个坑实测数据和仿真数据之间的差距比很多初学者想象中大得多。第一个坑是触发时刻不精确。实测数据要确定功率扰动发生的准确时刻通常用频率变化率或者功率突变检测但检测算法本身有延迟数据段的起始点因此可能偏移几十毫秒。ARMAX对起点偏移有一定鲁棒性因为A多项式拟合的是整个数据段的动态但如果起点偏差超过一个时间常数结果就会明显变差。第二个坑是频率信号的质量。PMU的频率通道有内部滤波和相位延迟不同厂家的PMU延迟特性还不一样这相当于在输入输出通路里加了一个未知延迟。处理手法是把nk从小调到大比如从1测试到5观察H估计值是否稳定。如果nk增大后结果变化很大说明原始信号延迟问题比较严重可以考虑先对齐事件触发点。第三个坑是数据段内包含多次扰动。如果评估窗口内还有其他机组跳闸或者负荷突变ARMAX会把多次扰动当做一个叠加输入辨识结果必然失真。解决方法是先做事件检测把每个扰动事件单独切段每次评估一个事件然后把多次事件的H估计结果取平均或中位数。第四个坑是最隐蔽的多机系统的惯量能否用一个等效H来概括。ARMAX辨识出来的本质上是“从测量点到系统之间的等效动态”它反映的是这个测量点能观察到的系统惯量水平而不是全部机组惯量的简单求和。做系统级惯量监测时这个值很有意义但不要试图用单点测量去还原每台机的惯量分布。4.4 独家避坑技巧这里分享几个常规文档里不常写、但实操时非常管用的经验。做参数提取时我强烈建议不要把H和D分开单独看。因为从式子里看H和D是耦合的噪声扰动可能让H偏低、D偏高但两者的乘积对应的时间常数却相对稳定。如果多次辨识的τ值都稳定在某个区间而H和D波动较大那大概率是数据信噪比不足而不是模型错了这时候参考τ比参考H更有意义。数据段长度的选择有个经验规律扰动后窗口建议覆盖至少3倍时间常数。真实时间常数未知怎么办先粗拟合一次估计出τ再回头调整窗口长度重新辨识一遍。这种“迭代式窗口选择”能让结果稳定不少。另外用多次蒙特卡洛仿真来校验算法的精度是个好习惯。同一个真值系统换不同随机种子生成十组噪声数据辨识十次取均值和标准差。这样你能直观看到算法的重复性。如果标准差偏大说明该工况下辨识条件不理想直接看单次结果没有意义。最后再分享一个我觉得很实用的扩展思路。这套方法不只适用于单机等值系统在微网、风电场并网点、孤岛系统这些场景里只要你能获取并网点的功率扰动数据和频率响应数据就能用同样的流程辨识等效惯量。我做风电场并网评估时把这个流程封装成了一个脚本输入是CSV格式的功率和频率信号输出就是惯量H、阻尼D和置信区间。工程上拿来做趋势监测完全够用。希望这份带完整代码的笔记能帮你少走几个弯路。后面你如果遇到辨识结果不稳或者数据预处理上的问题欢迎在评论区把情况贴出来我们一起讨论。
