4. 动量轮卸载的经典控制律设计

任务定义

本案例针对卫星动量轮长期运行后因累积角动量接近饱和而丧失姿态控制能力的问题，设计经典卸载控制律。目标为：当飞轮转速超过设定阈值时，自动生成卸载指令，驱动磁力矩器产生卸载力矩，使飞轮转速回落至安全区间，同时保证卸载过程不干扰卫星正常姿态指向。

系统建模与问题描述

考虑卫星本体携带一个沿本体z轴安装的反作用飞轮，其角动量为 h_w = J_w · Ω，其中 J_w 为飞轮转动惯量，Ω 为飞轮转速。卫星受到环境干扰力矩 T_d（如气动、重力梯度、太阳光压等）持续作用，导致飞轮不断累积角动量。当 |Ω| ≥ Ω_sat（饱和阈值）时，需启动卸载。

卸载执行机构为沿卫星本体三轴安装的磁力矩器，产生磁控力矩 T_m = m × B，其中 m 为磁矩指令，B 为地磁场矢量在星体坐标系的投影。由于磁力矩方向始终垂直于磁场方向，无法直接产生任意方向的卸载力矩，因此采用经典交叉积卸载控制律：

m = -k · (h_w × B) / |B|²

该控制律使磁力矩方向始终指向减小飞轮角动量的方向，实现能量最优卸载。

控制律设计与实现

完整卸载控制律包含以下环节：

1. 卸载触发判断：当 |Ω| > Ω_dead（死区上限）时启动卸载，当 |Ω| < Ω_dead_low（死区下限）时停止卸载，防止频繁切换；
2. 卸载力矩计算：根据当前飞轮角动量与磁场测量值，计算期望卸载力矩 T_unload = -k_p · h_w（比例卸载）；
3. 磁矩指令生成：将期望力矩投影到磁场垂直平面，生成实际磁矩指令 m_cmd；
4. 指令限幅：限制磁矩幅值不超过执行机构最大输出 m_max。

Python仿真代码

import numpy as np
import matplotlib.pyplot as plt

# ===== 卫星与执行机构参数 =====
J_w = 0.02          # 飞轮转动惯量 (kg·m²)
Omega_sat = 3000    # 饱和转速 (rpm)
Omega_dead_hi = 2500  # 卸载启动阈值 (rpm)
Omega_dead_lo = 1500  # 卸载停止阈值 (rpm)
k_p = 0.01          # 卸载比例系数 (N·m·s/rad)
m_max = 0.5         # 磁力矩器最大磁矩 (A·m²)
dt = 0.1            # 仿真步长 (s)
T_total = 2000      # 总仿真时间 (s)

# ===== 环境模型 =====
# 地磁场模拟（轨道角速度约0.001 rad/s，磁场随轨道变化）
def magnetic_field(t):
    # 简化的近极轨道磁场模型，单位特斯拉
    B0 = 3e-5
    omega_orbit = 0.001
    Bx = B0 * np.sin(omega_orbit * t)
    By = B0 * np.cos(omega_orbit * t)
    Bz = B0 * 0.3 * np.sin(2 * omega_orbit * t)
    return np.array([Bx, By, Bz])

# 环境干扰力矩（常值+周期项，单位N·m）
def disturbance_torque(t):
    Td = np.array([
        2e-5 * np.sin(0.002 * t),
        1.5e-5 * np.cos(0.001 * t),
        1e-5 * np.sin(0.003 * t)
    ])
    return Td

# ===== 卸载控制律 =====
def unload_controller(Omega, B_vec, t):
    """
    经典交叉积卸载控制律
    输入：飞轮转速(rpm)，磁场矢量(T)，时间(s)
    输出：磁矩指令(A·m²)
    """
    # 飞轮角动量 (N·m·s)
    h_w = J_w * Omega * (2 * np.pi / 60)  # rpm -> rad/s
    
    # 卸载触发判断
    if abs(Omega) < Omega_dead_hi:
        return np.zeros(3)  # 未超过阈值，不卸载
    
    # 期望卸载力矩（比例反馈）
    T_unload = -k_p * h_w * np.array([0, 0, 1])  # 仅卸载z轴飞轮
    
    # 交叉积控制律：m = (T_unload × B) / |B|²
    B_mag_sq = np.dot(B_vec, B_vec)
    if B_mag_sq < 1e-12:
        return np.zeros(3)
    
    m_cmd = np.cross(T_unload, B_vec) / B_mag_sq
    
    # 磁矩限幅
    m_norm = np.linalg.norm(m_cmd)
    if m_norm > m_max:
        m_cmd = m_cmd * (m_max / m_norm)
    
    return m_cmd

# ===== 仿真主循环 =====
def simulate():
    # 初始状态：飞轮转速略低于饱和
    Omega = 2800.0  # rpm
    time_arr = []
    omega_arr = []
    m_cmd_arr = []
    unloading_active = False
    
    for i in range(int(T_total / dt)):
        t = i * dt
        B = magnetic_field(t)
        Td = disturbance_torque(t)
        
        # 控制律计算
        m_cmd = unload_controller(Omega, B, t)
        
        # 磁力矩实际输出
        T_m = np.cross(m_cmd, B)
        
        # 飞轮动力学：J_w * dΩ/dt = T_m_z + Td_z
        # 注意：磁力矩z分量 + 干扰力矩z分量共同作用于飞轮
        dOmega_dt = (T_m[2] + Td[2]) / J_w
        Omega += dOmega_dt * dt * (60 / (2 * np.pi))  # 转换为rpm
        
        # 记录数据
        time_arr.append(t)
        omega_arr.append(Omega)
        m_cmd_arr.append(np.linalg.norm(m_cmd))
    
    return np.array(time_arr), np.array(omega_arr), np.array(m_cmd_arr)

# ===== 运行仿真 =====
time, omega, m_cmd = simulate()

# ===== 结果可视化 =====
plt.figure(figsize=(10, 6))
plt.subplot(2, 1, 1)
plt.plot(time / 60, omega, 'b-', linewidth=1.5)
plt.axhline(y=Omega_sat, color='r', linestyle='--', label='饱和阈值')
plt.axhline(y=Omega_dead_hi, color='g', linestyle='--', label='卸载启动阈值')
plt.axhline(y=Omega_dead_lo, color='orange', linestyle='--', label='卸载停止阈值')
plt.xlabel('时间 (min)')
plt.ylabel('飞轮转速 (rpm)')
plt.legend()
plt.grid(True)
plt.title('动量轮卸载控制律仿真')

plt.subplot(2, 1, 2)
plt.plot(time / 60, m_cmd, 'r-', linewidth=1)
plt.xlabel('时间 (min)')
plt.ylabel('磁矩指令幅值 (A·m²)')
plt.grid(True)

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

验证方法与判定标准

运行方式：在Python环境中安装 numpy 和 matplotlib 后直接运行上述代码。

预期输出：仿真结果应显示飞轮转速从初始2800 rpm开始，在卸载控制作用下逐渐下降，最终稳定在1500 rpm（卸载停止阈值）附近波动。磁矩指令幅值在卸载期间非零，且不超过0.5 A·m²限幅值。

判定标准：

1. 飞轮转速在2000秒仿真时间内从2800 rpm降至1500 rpm以下，且未超过3000 rpm饱和阈值；
2. 卸载启动后磁矩指令持续输出，转速降至死区下限后磁矩归零；
3. 若增大干扰力矩幅值，飞轮转速仍能被控制在安全范围内，说明控制律具有鲁棒性。