4. 基于反正切函数的转子位置提取方法
任务目标
本案例完成从滑模观测器（SMO）估计的反电动势（EMF）中，利用反正切函数提取转子位置 θ̂ 的完整实现，并定量分析提取误差的来源。核心目标是建立从 αβ 轴反电动势到转子电角度的最小可复现链路，并验证误差随转速和反电动势幅值变化的关系。
实现原理与误差来源
在表贴式PMSM的αβ静止坐标系下，SMO估计的反电动势分量 ê_α、ê_β 满足：
ê_α = -ψ_f·ω_e·sin(θ_e)，ê_β = ψ_f·ω_e·cos(θ_e)
转子位置可由四象限反正切直接提取：θ̂_e = atan2(-ê_α, ê_β)。误差来源主要有三类：
- 反电动势幅值缩放：
atan2对幅值不敏感，但若ê_α、ê_β存在非对称增益误差，会引入2倍频位置脉动； - 低通滤波器相位滞后：SMO输出经低通滤波后，反电动势相位滞后
Δφ = arctan(ω_e / ω_c)（ω_c为截止频率），该滞后直接叠加到位置估计上； - 滑模抖振与数字量化：符号函数切换带来的高频分量经滤波后残余纹波，以及ADC量化误差，表现为位置估计的高频噪声。
关键完整代码
import numpy as np
import matplotlib.pyplot as plt
# ---------- 参数定义 ----------
fs = 20000.0 # 采样频率 20kHz
Ts = 1.0 / fs
f_e = 100.0 # 电频率 100Hz（对应1500rpm@4极对）
omega_e = 2 * np.pi * f_e
psi_f = 0.05 # 永磁磁链 Wb
R = 0.5 # 定子电阻 Ω（用于SMO重构，此处仅作参数传递）
L = 0.008 # 电感 H
# 低通滤波器（一阶IIR，截止频率 fc = 500Hz）
fc = 500.0
alpha_lpf = 1.0 - np.exp(-2 * np.pi * fc * Ts)
# ---------- 真实反电动势（理想正弦） ----------
N = 2000 # 仿真0.1s
t = np.arange(N) * Ts
theta_e = omega_e * t
e_alpha_true = -psi_f * omega_e * np.sin(theta_e)
e_beta_true = psi_f * omega_e * np.cos(theta_e)
# ---------- 模拟SMO输出：真实EMF + 高频抖振噪声 ----------
# 抖振幅值取基波幅值的5%，频率为开关频率（10kHz）
chatter_amp = 0.05 * psi_f * omega_e
e_alpha_smo = e_alpha_true + chatter_amp * np.sign(np.sin(2 * np.pi * 10000 * t))
e_beta_smo = e_beta_true + chatter_amp * np.sign(np.cos(2 * np.pi * 10000 * t))
# ---------- 低通滤波（模拟SMO反电动势提取环节） ----------
e_alpha_filt = np.zeros(N)
e_beta_filt = np.zeros(N)
for k in range(1, N):
e_alpha_filt[k] = e_alpha_filt[k-1] + alpha_lpf * (e_alpha_smo[k] - e_alpha_filt[k-1])
e_beta_filt[k] = e_beta_filt[k-1] + alpha_lpf * (e_beta_smo[k] - e_beta_filt[k-1])
# ---------- 反正切位置提取 ----------
theta_hat = np.arctan2(-e_alpha_filt, e_beta_filt)
# 相位展开（消除atan2的±π跳变）
theta_hat_unwrapped = np.unwrap(theta_hat)
# ---------- 误差分析 ----------
# 理论相位滞后：Δφ = arctan(ω_e / ω_c)
omega_c = 2 * np.pi * fc
phase_lag = np.arctan(omega_e / omega_c) # 理论滞后量（rad）
# 实际误差（稳态段，取后50%数据）
steady_start = N // 2
theta_error = theta_hat_unwrapped[steady_start:] - theta_e[steady_start:]
# 去除平均偏差（直流偏置由相位滞后引起）
mean_error = np.mean(theta_error)
ripple_error = np.std(theta_error - mean_error)
print(f"理论相位滞后: {phase_lag:.4f} rad ({np.degrees(phase_lag):.2f}°)")
print(f"实际平均误差: {mean_error:.4f} rad ({np.degrees(mean_error):.2f}°)")
print(f"误差脉动(RMS): {ripple_error:.4f} rad ({np.degrees(ripple_error):.2f}°)")
# ---------- 可视化 ----------
plt.figure(figsize=(10, 6))
plt.subplot(2, 1, 1)
plt.plot(t[:500], theta_e[:500], 'b-', label='真实位置 θ_e', linewidth=2)
plt.plot(t[:500], theta_hat_unwrapped[:500], 'r--', label='提取位置 θ̂_e', linewidth=1.5)
plt.xlabel('时间 (s)'); plt.ylabel('电角度 (rad)'); plt.legend(); plt.grid(True)
plt.title('转子位置提取结果（前25ms）')
plt.subplot(2, 1, 2)
plt.plot(t[steady_start:], np.degrees(theta_error), 'g-', linewidth=1)
plt.xlabel('时间 (s)'); plt.ylabel('位置误差 (°)'); plt.grid(True)
plt.title(f'位置误差（稳态段，平均误差={np.degrees(mean_error):.2f}°，脉动RMS={np.degrees(ripple_error):.2f}°）')
plt.tight_layout()
plt.savefig('position_extraction_result.png', dpi=100)
plt.show()运行与验证
环境依赖：Python 3.8+，需安装 numpy 和 matplotlib。
运行方式：将上述代码保存为 smo_position_extraction.py，在终端执行：
python smo_position_extraction.py通过条件：
- 控制台输出的实际平均误差与理论相位滞后之差小于0.01 rad（即滤波器相位滞后模型准确）；
- 误差脉动RMS小于0.02 rad（约1.15°），证明抖振抑制有效；
- 图中蓝色真实位置与红色提取位置曲线在稳态段基本重合，误差曲线呈现以平均误差为中心的窄带波动。
误差来源定量验证：将 fc 从500改为200Hz重新运行，观察平均误差增大（因 ω_c 减小导致 Δφ 增大）；将 chatter_amp 系数从0.05改为0.15重新运行，观察脉动RMS增大。这两组对比实验可分别验证滤波器滞后和抖振对位置提取精度的影响。