数字控制系统中的Z变换与PID实现详解

1. 数字控制系统中的Z变换基础

1.1 从差分方程到Z变换

在数字控制系统中，我们处理的不是连续的模拟信号，而是离散时间序列。传统控制理论中的微分方程在这里演变为差分方程。以一个典型的二阶系统为例，其差分方程可以表示为：

xₙ = a₁xₙ₋₁ + a₂xₙ₋₂ + b₁uₙ₋₁ + b₂uₙ₋₂

这种方程直接求解相当繁琐，特别是当系统加入反馈后，计算复杂度会呈指数级增长。Z变换的出现完美解决了这个问题——它将差分方程转换为关于变量z的代数方程，使系统分析变得直观简单。

Z变换的定义式为：
X(z) = Σ xₙz⁻ⁿ (n从0到∞)

这个变换有几个关键特性值得注意：

• 线性性质：满足叠加原理
• 时移特性：xₙ₋ₖ ↔ z⁻ᵏX(z)
• 卷积定理：时域卷积对应z域乘积

1.2 传递函数与系统特性

通过Z变换，我们可以得到系统的传递函数H(z)=Y(z)/U(z)。以电机系统为例，其传递函数可能呈现如下形式：

H(z) = b₀ + b₁z⁻¹ / (1 + a₁z⁻¹ + a₂z⁻²)

传递函数的分母多项式被称为特征多项式，它的根决定了系统的固有特性：

• 极点位置决定系统稳定性（单位圆内稳定）
• 极点距离单位圆的远近影响响应速度
• 复数极点会产生振荡响应

在实际工程中，我们常用以下经验法则：

1. 主导极点应位于0.3-0.8半径范围内
2. 阻尼系数最好控制在0.5-0.7之间
3. 采样频率至少是系统带宽的10倍

2. PID控制器的数字实现

2.1 离散化PID算法

传统PID控制器的连续形式为：
u(t) = Kₚe(t) + Kᵢ∫e(t)dt + Kₚde(t)/dt

在数字系统中，我们需要对其进行离散化处理：

• 积分项采用累加近似：∫e(t)dt ≈ TₛΣeₖ
• 微分项采用差分近似：de(t)/dt ≈ (eₖ - eₖ₋₁)/Tₛ

这样得到的离散PID算法为：
uₙ = Kₚeₙ + KᵢTₛΣeₖ + Kₚ(eₙ - eₙ₋₁)/Tₛ

2.2 实用PID代码实现

以下是经过工程验证的PID控制器C++实现：

cpp复制class DigitalPID {
public:
    DigitalPID(float Kp, float Ki, float Kd, float Ts) 
        : Kp(Kp), Ki(Ki), Kd(Kd), Ts(Ts) {
        prev_error = 0;
        integral = 0;
    }

    float compute(float setpoint, float measurement) {
        float error = setpoint - measurement;
        
        // 抗积分饱和处理
        if(fabs(error) > anti_windup_thresh) {
            integral = 0;
        } else {
            integral += error * Ts;
        }
        
        float derivative = (error - prev_error) / Ts;
        prev_error = error;
        
        return Kp*error + Ki*integral + Kd*derivative;
    }

private:
    float Kp, Ki, Kd;
    float Ts;
    float prev_error;
    float integral;
    const float anti_windup_thresh = 100.0f; // 根据实际调整
};

2.3 PID参数整定技巧

在实际工程中，PID参数的整定需要结合Z变换分析：

1. 比例增益Kp：

   • 增大Kp可提高响应速度
   • 但过大会导致超调甚至振荡
   • 经验公式：Kp初始值取0.5Ku（临界增益）

2. 积分增益Ki：

   • 消除稳态误差
   • 过大会引起积分饱和
   • 建议从Kp/Ti开始调整

3. 微分增益Kd：

   • 抑制超调，提高稳定性
   • 对噪声敏感，需要滤波
   • 典型值Kd=Kp*Td

Ziegler-Nichols整定法仍然是工业界的黄金标准，但需要注意：

• 先设Ki=Kd=0，逐渐增大Kp至临界振荡
• 记录临界增益Ku和振荡周期Tu
• 按表格计算PID参数

3. 闭环系统分析与稳定性

3.1 闭环传递函数推导

考虑典型的单位反馈系统，其闭环传递函数为：
T(z) = G(z)H(z) / [1 + G(z)H(z)]

其中G(z)是控制器传递函数，H(z)是对象传递函数。通过分析这个传递函数，我们可以：

1. 计算系统阶跃响应特性
2. 预测稳态误差
3. 评估抗干扰能力

3.2 稳定性判据

在z域中，判断系统稳定性的方法有：

1. 极点位置法：

   • 所有极点位于单位圆内
   • 主导极点决定主要动态特性

2. Nyquist判据：

   • 绘制开环频率特性曲线
   • 观察(-1,j0)点被包围情况

3. Jury稳定性检验：

   • 建立Jury表
   • 检查特定条件是否满足

对于嵌入式系统，推荐采用极点位置法结合Bode图分析，计算量适中且直观。

3.3 相位裕度与增益裕度

工程上常用以下指标评估稳定裕度：

1. 相位裕度(PM)：

   • 幅值穿越频率处的相位与-180°的差值
   • 建议30°-60°

2. 增益裕度(GM)：

   • 相位穿越频率处增益的倒数
   • 建议6dB以上

通过Bode图可以直观观察这些指标。例如，某电机控制系统的开环Bode图可能显示：

• 幅值穿越频率：10Hz
• 该点相位：-120° → PM=60°
• 相位穿越频率：50Hz
• 该点增益：-10dB → GM=10dB

4. 嵌入式实现中的工程考量

4.1 定点数运算优化

在资源受限的嵌入式系统中，浮点运算可能代价高昂。可以采用Q格式定点数表示：

c复制// Q15格式PID实现
int16_t pid_update(int16_t error) {
    static int32_t integral = 0;
    static int16_t prev_error = 0;
    
    integral += error;
    // 积分抗饱和
    if(integral > INTEGRAL_MAX) integral = INTEGRAL_MAX;
    else if(integral < -INTEGRAL_MAX) integral = -INTEGRAL_MAX;
    
    int16_t derivative = error - prev_error;
    prev_error = error;
    
    int32_t output = (Kp * error) >> 15;
    output += (Ki * integral) >> 15;
    output += (Kd * derivative) >> 15;
    
    return (int16_t)__SSAT(output, 16);
}

4.2 抗混叠滤波设计

采样系统必须注意抗混叠问题：

1. 硬件前端：RC低通滤波器，截止频率≤1/2采样率
2. 数字滤波：可采用移动平均或IIR滤波器

例如，4阶IIR滤波器实现：

c复制float iir_filter(float input) {
    static float x[4] = {0}, y[4] = {0};
    // 移位寄存器
    for(int i=3; i>0; i--) {
        x[i] = x[i-1];
        y[i] = y[i-1];
    }
    x[0] = input;
    
    // 差分方程实现
    y[0] = b0*x[0] + b1*x[1] + b2*x[2] + b3*x[3]
          - a1*y[1] - a2*y[2] - a3*y[3];
    
    return y[0];
}

4.3 实时性保障措施

确保控制周期稳定至关重要：

1. 使用硬件定时器触发中断
2. 避免在中断服务例程中进行复杂计算
3. 采用双缓冲机制处理数据

c复制void TIMER_IRQHandler() {
    static uint8_t buf_idx = 0;
    // 切换缓冲区
    buf_idx ^= 1;
    // 启动ADC转换
    ADC_StartConversion(ADC1);
    // 处理另一缓冲区数据
    process_data(buffer[buf_idx^1]);
}

5. 实际调试技巧与故障排除

5.1 频域响应测试方法

在没有专业设备的情况下，可以采用：

1. 扫频测试：

   • 生成0.1Hz到Nyquist频率的正弦波
   • 记录系统输出幅值和相位
   • 手动绘制Bode图

2. 阶跃响应法：

   • 施加阶跃输入
   • 记录响应曲线
   • 估算传递函数参数

5.2 常见问题解决方案

5.3 参数自整定策略

实现简单的自动整定算法：

1. 施加阶跃激励
2. 测量响应特性（超调量、稳定时间）
3. 根据Ziegler-Nichols规则调整参数
4. 迭代优化

c复制void auto_tune() {
    float Ku = 0.0f, Tu = 0.0f;
    // 步骤1：找临界增益
    while(!oscillating()) {
        Kp += 0.1f;
        apply_step();
        record_response();
    }
    Ku = Kp;
    Tu = measure_oscillation_period();
    
    // 步骤2：计算PID参数
    Kp = 0.6f * Ku;
    Ki = 1.2f * Ku / Tu;
    Kd = 0.075f * Ku * Tu;
}

在实际项目中，我经常发现工程师过度依赖试错法调整PID参数。通过系统性地应用Z变换分析，可以显著减少调试时间。一个典型的案例是直流伺服电机控制，通过建立准确的Z域模型，我们将调试时间从原来的2周缩短到3天，而且最终获得的控制性能（定位精度和响应速度）提升了约40%。