3. 饱和函数替代符号函数

任务目标

在FOC无感控制中，SMO（滑模观测器）使用符号函数sign(x)作为开关函数，导致反电动势估计中存在高频抖振。本案例将符号函数替换为饱和函数sat(x)，通过边界层厚度参数调节切换增益，在保持观测器收敛性的前提下抑制抖振，并对比两种方案的反电动势观测波形。

核心原理

饱和函数定义为：

sat(x) = x/ε, 当 |x| ≤ ε
sat(x) = sign(x), 当 |x| > ε

其中ε为边界层厚度。当ε→0时，饱和函数退化为符号函数。ε越大，抖振抑制越强，但观测误差也随之增大。实际应用中需折中选取ε。

完整实现代码

import numpy as np
import matplotlib.pyplot as plt

# ========== 电机与观测器参数 ==========
Rs = 2.0          # 定子电阻(Ω)
Ls = 0.01         # 定子电感(H)
Ke = 0.1          # 反电动势系数(V·s/rad)
J = 0.002         # 转动惯量(kg·m²)
B = 0.001         # 阻尼系数(N·m·s/rad)
pole_pairs = 4    # 极对数

# 滑模增益与饱和边界层
K_smo = 50.0      # 滑模切换增益
epsilon = 0.5     # 饱和函数边界层厚度

# 仿真参数
dt = 1e-5         # 步长(s)
t_end = 0.2       # 仿真时长(s)
t = np.arange(0, t_end, dt)
n = len(t)

# 电气时间常数远小于机械，采用双速率：电气步长1e-5，机械每100步更新
mech_step = 100

# ========== 参考转速与负载 ==========
speed_ref = 100.0 * 2 * np.pi / 60  # 100rpm → rad/s
TL = 0.05          # 负载转矩(N·m)

# ========== 状态初始化 ==========
# 真实状态
theta = 0.0
omega = 0.0
i_alpha, i_beta = 0.0, 0.0
e_alpha, e_beta = 0.0, 0.0

# 观测器状态
i_alpha_hat, i_beta_hat = 0.0, 0.0
e_alpha_hat, e_beta_hat = 0.0, 0.0

# 存储
results = {
    'sign': {'e_alpha': [], 'e_beta': [], 'i_err_alpha': []},
    'sat':  {'e_alpha': [], 'e_beta': [], 'i_err_alpha': []}
}

# ========== 符号函数与饱和函数 ==========
def sign_func(x):
    return np.sign(x)

def sat_func(x, eps):
    return np.clip(x / eps, -1.0, 1.0)

# ========== 仿真主循环 ==========
for method in ['sign', 'sat']:
    # 重置状态
    theta, omega = 0.0, 0.0
    i_alpha, i_beta = 0.0, 0.0
    e_alpha, e_beta = 0.0, 0.0
    i_alpha_hat, i_beta_hat = 0.0, 0.0
    e_alpha_hat, e_beta_hat = 0.0, 0.0
    
    for k in range(n):
        # --- 电流采样（真实电机模型） ---
        # 反电动势（真实值）
        e_alpha = -Ke * omega * np.sin(theta)
        e_beta = Ke * omega * np.cos(theta)
        
        # 定子电压（假设采用id=0控制，v_alpha/beta由电流环输出）
        # 简化：直接给定正弦电压驱动电机旋转
        v_alpha = 5.0 * np.cos(omega * k * dt)
        v_beta = 5.0 * np.sin(omega * k * dt)
        
        # 电流微分方程
        di_alpha = (v_alpha - Rs * i_alpha - e_alpha) / Ls
        di_beta = (v_beta - Rs * i_beta - e_beta) / Ls
        i_alpha += di_alpha * dt
        i_beta += di_beta * dt
        
        # --- 滑模观测器 ---
        # 电流误差
        err_alpha = i_alpha_hat - i_alpha
        err_beta = i_beta_hat - i_beta
        
        # 切换函数
        if method == 'sign':
            switch_alpha = sign_func(err_alpha)
            switch_beta = sign_func(err_beta)
        else:
            switch_alpha = sat_func(err_alpha, epsilon)
            switch_beta = sat_func(err_beta, epsilon)
        
        # 观测器电流方程
        di_alpha_hat = (v_alpha - Rs * i_alpha_hat - e_alpha_hat - K_smo * switch_alpha) / Ls
        di_beta_hat = (v_beta - Rs * i_beta_hat - e_beta_hat - K_smo * switch_beta) / Ls
        i_alpha_hat += di_alpha_hat * dt
        i_beta_hat += di_beta_hat * dt
        
        # 反电动势观测（等效控制输出，经低通滤波）
        lpf_coeff = 2 * np.pi * 500 * dt  # 500Hz低通
        e_alpha_hat += lpf_coeff * (K_smo * switch_alpha - e_alpha_hat)
        e_beta_hat += lpf_coeff * (K_smo * switch_beta - e_beta_hat)
        
        # --- 机械更新（每mech_step步） ---
        if k % mech_step == 0:
            # 电磁转矩（由真实电流与反电动势计算）
            Te = 1.5 * pole_pairs * (e_alpha * i_alpha + e_beta * i_beta) / omega if omega != 0 else 0
            domega = (Te - B * omega - TL) / J
            omega += domega * (dt * mech_step)
            theta += omega * (dt * mech_step)
        
        # 存储（每0.1ms存一次）
        if k % 10 == 0:
            results[method]['e_alpha'].append(e_alpha)
            results[method]['e_beta'].append(e_beta)
            results[method]['i_err_alpha'].append(err_alpha)

# ========== 绘图对比 ==========
fig, axes = plt.subplots(2, 2, figsize=(12, 8))
t_plot = t[::10] * 1000  # ms

# 反电动势α轴对比
axes[0,0].plot(t_plot, results['sign']['e_alpha'], 'r-', lw=0.8, label='真实值')
axes[0,0].plot(t_plot, results['sign']['e_alpha'], 'b--', lw=0.8, alpha=0.7, label='符号函数观测')
axes[0,0].set_title('符号函数 - 反电动势α轴')
axes[0,0].set_xlabel('时间(ms)')
axes[0,0].legend()
axes[0,0].grid(True)

axes[0,1].plot(t_plot, results['sat']['e_alpha'], 'r-', lw=0.8, label='真实值')
axes[0,1].plot(t_plot, results['sat']['e_alpha'], 'g--', lw=0.8, alpha=0.7, label='饱和函数观测')
axes[0,1].set_title(f'饱和函数(ε={epsilon}) - 反电动势α轴')
axes[0,1].set_xlabel('时间(ms)')
axes[0,1].legend()
axes[0,1].grid(True)

# 电流误差对比
axes[1,0].plot(t_plot, results['sign']['i_err_alpha'], 'b-', lw=0.6)
axes[1,0].set_title('符号函数 - 电流观测误差α轴')
axes[1,0].set_xlabel('时间(ms)')
axes[1,0].grid(True)

axes[1,1].plot(t_plot, results['sat']['i_err_alpha'], 'g-', lw=0.6)
axes[1,1].set_title(f'饱和函数(ε={epsilon}) - 电流观测误差α轴')
axes[1,1].set_xlabel('时间(ms)')
axes[1,1].grid(True)

plt.tight_layout()
plt.savefig('smo_sat_comparison.png', dpi=150)
plt.show()

# ========== 定量评估 ==========
# 稳态段（后50%数据）
steady_start = len(results['sign']['e_alpha']) // 2

def calc_ripple(data):
    seg = np.array(data[steady_start:])
    return np.std(seg)

print("=== 抖振抑制效果对比（稳态段） ===")
print(f"符号函数  反电动势α轴标准差: {calc_ripple(results['sign']['e_alpha']):.4f} V")
print(f"饱和函数  反电动势α轴标准差: {calc_ripple(results['sat']['e_alpha']):.4f} V")
print(f"符号函数  电流误差α轴标准差: {calc_ripple(results['sign']['i_err_alpha']):.6f} A")
print(f"饱和函数  电流误差α轴标准差: {calc_ripple(results['sat']['i_err_alpha']):.6f} A")

运行与验证

运行环境：Python 3.8+，需安装NumPy和Matplotlib。

运行方式：直接执行脚本，将生成对比图smo_sat_comparison.png并在控制台输出定量指标。

通过判定标准：

图中饱和函数方案的反电动势观测曲线毛刺明显少于符号函数方案；
控制台输出的稳态标准差中，饱和函数方案的反电动势和电流误差标准差均小于符号函数方案；
两种方案的反电动势观测均值均能跟踪真实值，无发散。

参数调节说明

增大epsilon可进一步抑制抖振，但会增大稳态观测误差；减小epsilon则反之。实际应用中建议从epsilon = 0.1~1.0范围开始调试，观察电流误差波形确定合适边界层厚度。

📚 课程目录