☰
【电机滤波例程6】平方根扩展卡尔曼滤波(SR-EKF)原理与MATLAB例程:PMSM状态估计与QR协方差因子递推。附完整代码的下载链接
2026/10/10 2:45:17 网站建设 项目流程

`原创代码,请勿翻卖

文章目录

  • 程序简介
    • 算法原理
    • 实现流程
    • 输出说明
  • 运行结果
  • MATLAB源代码

程序简介

算法原理

采用表贴式永磁同步电机模型,满足Ld=Lq=p.L。通过已知αβ电压驱动机电状态预测,以αβ电流为量测修正电流、机械转速和电角度。 扩展卡尔曼滤波利用离散模型雅可比近似传播局部状态误差。 仿真采用非零初始转速和明确的初始估计误差,展示转速与电角度的估计收敛及变速跟踪过程。

实现流程

  1. 设置参数、固定随机种子,生成电压、真实电机状态和带噪电流。
  2. 初始化带误差的估计状态,核对解析雅可比。
  3. 逐步执行模型预测和当前案例的量测更新。
  4. 检查有限值与协方差正定性,统计指标并输出Figure。

输出说明

result.time为1×N时间;truth、estimate、error为状态数×N矩阵;observation和voltage为2×N矩阵。result.stateNames明确状态顺序:[i_alpha, i_beta, omega_m, theta_e],第5维为负载转矩T_L。

result.metrics保存统计指标;result.diagnostic保存新息、NIS及案例专属诊断。初始样本尚未量测更新,相关新息以NaN表示。

电压第k列驱动[t(k),t(k+1)]区间,电流第k列对应t(k);最后一列电压复制前一列,仅为输出长度一致,不用于额外传播。

各图窗的主题见“运行结果”。命令行仅输出该节列出的指标,其余误差和数值诊断从result.metrics与result.diagnostic读取。 转速和电角度主指标统计后80%仿真区间;电角度误差先按圆周环绕处理,再计算RMSE。

result.statistics.full保存全时段误差统计,post保存分区间统计,均含样本数、RMSE、均值、方差、标准差、平均绝对误差和最大绝对误差。方差采用N归一化;均值表示估计偏差。电角度先环绕再统计,参数使用物理量而非对数状态。current对比相同有效时刻的观测和滤波误差,缺失电流不计入观测统计。全时段包含初始样本;角度不参加信号算术均值/方差比较。

运行结果

转速与位置估计图:对比真实与估计机械转速,并放大展示电角度局部曲线;机械转速以r/min显示,电角度以rad显示。

估计误差图:显示机械转速与环绕处理后的电角度误差,用于观察收敛过程和变速响应。

电流观测与估计图:分别展示α、β电流的真值、带噪观测和滤波估计,观察电流平滑与跟踪情况。

负载转矩估计图:对比负载转矩真值和估计值,并给出负载估计误差,展示负载阶跃跟踪。

平方根数值诊断图:展示协方差最小特征值与平方根因子重构误差,核查因子递推的数值一致性。

命令行输出:对比全时段与分区间的RMSE、误差均值和方差,列出真值与估计值统计、电流滤波效果及数值诊断。

MATLAB源代码

部分代码如下:

% PMSM平方根EKF:独立MATLAB例程% 2026-09-26/Ver1,作者:matlabfilter,V同号,可获取代码定制、讲解等% 直接点击“运行”。只依赖基础 MATLAB;本文件末尾包含全部局部函数。% 固定随机种子0,直接运行并绘图;主输出为工作区中的 result。clear;clc;close all;rng(0);%% 集中参数:SI单位;wm为机械角速度,theta为电角度p.Rs=0.45;p.L=3e-3;p.psi=0.075;p.pairs=4;p.J=2e-3;p.B=5e-4;p.dt=1e-4;p.duration=0.65;p.ratedSpeed=150;t=0:p.dt:p.duration;count=numel(t);%% 真值与观测生成(理想有编码器驱动台架,独立于待验证估计器)% 台架的理想电流/速度闭环只用于生成物理一致的电压、电流数据。% 估计器不使用台架的真实角度、转速、速度给定或未知负载。% 本例验证无位置传感器“估计”,并不声称实现无传感器闭环启动。nx=5;plant=[0;0;80;0.3];truth=zeros(nx,count);voltage=zeros(2,count);currentObs=zeros(2,count);loadTrue=zeros(1,count);noiseStd=0.04*ones(2,count);outlierMask=false(1,count);speedReference=80+25*(1+tanh((t-0.16)/0.035))/2...-15*(1+tanh((t-0.44)/0.035))/2;loadTrue(:)=0.35;loadTrue(t>=0.25)=0.60;loadTrue(t>=0.43)=0.42;truth(1:4,1)=plant;ifnx==5,truth(5,:)=loadTrue;endcurrentObs(:,1)=plant(1:2)+noiseStd(:,1).*randn(2,1);speedIntegral=0;fork=2:count% 真实台架:有编码器FOC等效电压源;电压限幅为±60 V。wm=plant(3);angle=plant(4);rot=[cos(angle),sin(angle);-sin(angle),cos(angle)];idq=rot*plant(1:2);speedError=speedReference(k-1)-wm;speedIntegral=max(-8,min(8,speedIntegral+p.dt*speedError));iqRef=max(-10,min(10,0.28*speedError+3.0*speedIntegral+0.8));udq=[p.Rs*idq(1)-p.pairs*wm*p.L*idq(2);...p.Rs*idq(2)+p.pairs*wm*(p.L*idq(1)+p.psi)]...+2.8*([0;iqRef]-idq);u=max(-60,min(60,rot'*udq));voltage(:,k-1)=u;% 真值采用4个RK4子步,未知负载只传给真实电机。forsub=1:4h=p.dt/4;a=motorODE(plant,u,p,loadTrue(k-1));b=motorODE(plant+h*a/2,u,p,loadTrue(k-1));c=motorODE(plant+h*b/2,u,p,loadTrue(k-1));d=motorODE(plant+h*c,u,p,loadTrue(k-1));plant=plant+h*(a+2*b+2*c+d)/6;endplant(4)=wrapAngle(plant(4));truth(1:4,k)=plant;currentObs(:,k)=plant(1:2)+noiseStd(:,k).*randn(2,1);endvoltage(:,end)=voltage(:,end-1);available=true(1,count);currentObs(:,~available)=NaN;%% EKF初始化(故意设置初始误差,不从真值数组初始化)x=[0.15;-0.12;87;0.42];P=diag([0.3^2,0.3^2,12^2,0.18^2]);Q=diag([2e-6,2e-6,2e-3,2e-8]);ifnx==5x=[x;0.20];P=blkdiag(P,0.25^2);Q=blkdiag(Q,2e-5);endR=0.04^2*eye(2);H=[eye(2),zeros(2,nx-2)];estimate=zeros(nx,count);estimate(:,1)=x;innovation=nan(2,count);nis=nan(1,count);diagnostic=ones(2,count);rStdEstimate=0.04*ones(2,count);minEigenvalue=zeros(1,count);minEigenvalue(1)=min(eig(P));symmetryError=zeros(1,count);sqrtError=zeros(1,count);Lcov=chol(P,'lower');Lq=chol(Q,'lower');% 对离散预测函数作中心差分,核对链式解析雅可比。xcheck=x;xcheck(4)=0.6;jacobianError=checkJacobian(xcheck,[1;24],p,0.35);assert(jacobianError<1e-6,'离散模型雅可比检查失败');assert(abs(wrapAngle((-pi+0.01)-(pi-0.01))-0.02)<1e-12);%% 预测与量测更新:滤波器只使用电压、电流及已知参数%% 命令行统计:方差按N归一化,有限观测按共同有效时刻比较

完整代码:
https://download.csdn.net/download/callmeup/93509424

或前往专栏查看更多代码: https://blog.csdn.net/callmeup/category_13213136.html

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询