差分变换法与双线性变换法的基本原理和代码实现

1.差分变换法

差分变换法就是把微分方程中的导数用有限差分来近似等效，得到一个与原微分方程逼近的差分方程。

差分变换法包括后向差分与前向差分。

1.1 后向差分法

差分变换如下：
$$\frac{de(t)}{dt}=\frac{e(kT)-e(kT-T)}{T}$$
对两边进行z变换，可得：
$$s=\frac{1-z^{-1}}{T}$$
后向差分变换法又称后向矩阵积分法，即用后向矩形面积来代替。因此当采样周期较大时，这种变换方法精度变差。

由变换过程可知：

• 若D(s)稳定，则D(z)一定稳定；对于不稳定的D(s)，D(z)也有可能稳定。
• 由于后向差分变换不再满足z变换定义z=e^(sT)，因此s平面与z平面的映射关系发生了改变，数字控制器D(z)的频率响应产生了较大的畸变。

严重的频率映射畸变导致变换精度下降，使后向差分变换的应用受到了一定限制，但当系统性能要求不是很高时，后向差分也起到一定作用。

例：
$$D(s)=\frac{20s+80}{s+10},T=0.015$$
采用后向差分，代入下式：
$$s=\frac{1-z^{-1}}{T}$$
可得：
$$D(z)=\frac{U(z)}{E(z)}=\frac{21.2-20z^{-1}}{1.15-z^{-1}}$$
进行z反变换，可得：
$$u(k)=0.87u(k-1)+18.43e(k)-17.39e(k-1)$$
进一步探究采样频率对系统映射的影响，matlab代码如下：

numerator = [20, 80];     
denominator = [1, 10]; 
G = tf(numerator, denominator);

G_discrete1 = c2d(G, 0.01, 'foh'); 
G_discrete2 = c2d(G, 0.001, 'foh');
G_discrete3 = c2d(G, 0.0001, 'foh');


figure;
hold on; 

bode(G, G_discrete1, G_discrete2, G_discrete3);

所得bode图为：

由bode图可知，当采样频率为0.01时，离散系统有较大影响；当采样频率为0.001甚至更高时，离散系统基本能反应连续系统的特性。

通过后向差分设计一个PID控制器，C语言伪代码实现如下：

// PID参数
#define KP 1.0f    // 比例增益
#define KI 0.1f    // 积分增益
#define KD 0.05f   // 微分增益

// 控制变量
float setpoint = 100.0f;     // 目标速度
float currentSpeed = 0.0f;   // 实际速度
float previousSpeed = 0.0f;  // 上一速度
float integral = 0.0f;       // 积分值
float lastTime = 0.0f;       // 上次计算时间
float deltaTime = 0.1f;      // 计算周期，假设为0.1秒

void PID_Controller() {

    float error = setpoint - currentSpeed;

    integral += error * deltaTime;

    float derivative = (currentSpeed - previousSpeed) / deltaTime;

    float output = (KP * error) + (KI * integral) - (KD * derivative);

    // 更新电机控制信号（假设控制信号范围为0-255）
    if (output > 255) {
        output = 255;
    } else if (output < 0) {
        output = 0;
    }

    setMotorSpeed((uint8_t)output);

    previousSpeed = currentSpeed;
}


void readCurrentSpeed() {
    currentSpeed = getMotorSpeedFromEncoder(); 
}

int main(void) {
    while (1) {
        readCurrentSpeed(); 
        PID_Controller();   
        Delay(100);     
    }
}

verilog实现如下，输出限幅等功能未实现：

module PID_Controller (
    input wire clk,               		// 时钟信号
    input wire rst,             	    // 复位信号
    input wire [15:0] setpoint,  	    // 目标值
    input wire [15:0] current_value,    // 当前值
    output reg [15:0] control_output    // 控制输出
);

parameter [15:0] KP = 16'd2; // 比例增益
parameter [15:0] KI = 16'd1; // 积分增益
parameter [15:0] KD = 16'd1; // 微分增益

reg [15:0] error;            // 误差
reg [15:0] previous_value;   // 上一时刻的当前值
reg [31:0] integral;         // 积分值
reg [15:0] derivative;       // 微分值

always @(posedge clk or posedge rst) begin
    if (rst) begin
        control_output <= 16'd0;
        previous_value <= 16'd0; 
        integral <= 32'd0;        
    end else begin
        error <= setpoint - current_value;

        integral <= integral + error;

        derivative <= current_value - previous_value;

        control_output <= (KP * error) + (KI * integral[15:0]) - (KD * derivative);

        previous_value <= current_value;
    end
end
endmodule

1.2 前向差分法

前向差分变换原理：
$$\frac{de(t)}{dt}=\frac{e(kT+T)-e(kT)}{T}$$
进行z变换可得：
$$s=\frac{z-1}{T}$$
前向差分性质如下：

• 若D(s)稳定，则D(z)不一定稳定。
• 数字控制器D(z)的频率响应会产生较大畸变。

因此，前向差分很少被使用。

2.双线性变换法

双线性变化法又被称为Tustin变换法，推导过程如下：
$$z=e^{Ts}=e^{\frac{Ts/2}{-Ts/2}}$$
对e^(Ts/2)部分进行泰勒级数展开，并取前两项近似：
$$z=\frac{2/T+s}{2/T-s}$$
得到双线性变换法的计算公式为：
$$s=\frac{2}{T}\frac{z-1}{z+1}$$
双线性变换的特点为：

• 若D(s)稳定，则D(z)也稳定。
• 变换前后的频率响应发生畸变。
• 不存在频率混叠现象。

对matlab函数c2d进行简要说明;

c2d：将模型从连续时间转换为离散时间

sysd = c2d(sysc,Ts,method) 指定离散化方法

离散化方法，指定为以下值之一：

• 'zoh' - 零阶保持（默认值）。假定控制输入在采样时间周期 Ts 内是分段常量。
• 'foh' - 三角形逼近（修正一阶保持）。假定控制输入在采样时间周期 Ts 内是分段线性值。
• 'impulse' - 冲激不变离散化
• 'tustin' - 双线性（突斯汀）方法。要使用频率预修正指定此方法（以前称为 'prewarp' 方法），请使用 c2dOptions 的 PrewarpFrequency 选项。
• 'matched' - 零极点匹配法
• 'least-squares' - 最小二乘法
• 'damped' - 仅基于稀疏模型的 TRBDF2 公式的阻尼突斯汀逼近。

首先对比相同采样频率下使用后向差分与双线性变换的性能影响：

numerator = [20, 80];     
denominator = [1, 10]; 
G = tf(numerator, denominator);

G_discrete1 = c2d(G, 0.015, 'foh'); 
G_discrete2 = c2d(G, 0.015, 'tustin');


figure;
hold on; 

bode(G, G_discrete1, G_discrete2);

由bode图可知，在该采样频率下，双线性变换与差分变化作用效果近似。

将采样频率降低到0.15：

numerator = [20, 80];     
denominator = [1, 10]; 
G = tf(numerator, denominator);

G_discrete1 = c2d(G, 0.15, 'foh'); 
G_discrete2 = c2d(G, 0.15, 'tustin');


figure;
hold on; 

bode(G, G_discrete1, G_discrete2);

观察bode图，此时使用双线性变化效果优于后向差分：

使用chatGPT对比这两种方法：

后向差分（Backward Difference）和双线性差分（Bilinear Difference）是数值分析中常用的差分方法，它们在不同的场景下各有优势。下面分别分析这两种方法的适用场景及优劣。

• 后向差分的优势场景

1. 稳定性要求高的系统：
   • 后向差分法是稳定的，尤其适合处理刚性系统。刚性系统的特点是某些动态特性变化非常快，这使得普通的显式方法难以稳定计算。后向差分能有效避免不稳定现象。
2. 计算时间延迟的情况：
   • 在一些实时系统中，后向差分法可以使用最新的已知数据点来估计当前状态，这在需要延迟计算的系统中非常有用，如控制系统中的实时反馈。
3. 简单实现：
   • 对于线性时间不变系统，后向差分的实现较为简单，特别是在不需要复杂的求解算法时。
4. 误差控制：
   • 后向差分通常能更好地控制截断误差，尤其在时间步长较大的情况下，避免了由大时间步带来的数值误差累积。

• 双线性差分的优势场景

1. 高频动态响应：
   • 双线性差分法能够更好地处理高频动态系统，尤其在需要高精度模拟频率响应特性的情况下，能够较好地保留系统的动态特性。
2. 非线性系统的处理：
   • 在处理非线性系统时，双线性差分法的表现往往优于后向差分。因为它能更有效地捕捉系统的变化，提供更准确的结果。
3. 数值稳定性与精度的折中：
   • 双线性差分法在保持数值稳定性的同时，能够提供更高的精度。特别是在需要优化精度的应用场景下，如电子电路模拟、信号处理等，双线性差分具有明显优势。
4. 适应多种时间步长：
   • 双线性差分能够更好地适应不同的时间步长，能够灵活应对多种动态情况，因此在自适应算法中较为常用。

• 总结

• 后向差分适合于需要稳定性和较简单实现的场合，尤其是在刚性系统和实时反馈控制中。
• 双线性差分适合于高频动态响应、非线性系统以及需要高精度的应用场合。

简而言之，后向差分已经适用于大部分场景，针对特殊情况，再考虑使用双线性差分。

C语言代码实现例子：

#include <stdio.h>

double Kp = 1.0; // 比例增益
double Ki = 0.5; // 积分增益
double Kd = 0.1; // 微分增益

// 采样周期
double Ts = 0.01;

// PID控制器状态
double integral = 0.0; // 积分项
double prev_error = 0.0; // 上一时刻的误差


double PID_Controller(double setpoint, double measured_value) {

    double error = setpoint - measured_value;
    integral += error * Ts;
    double derivative = (2 * (error - prev_error)) / Ts; // 双线性差分法计算微分
    double output = Kp * error + Ki * integral + Kd * derivative;
    
    prev_error = error;

    return output;
}

verilog代码：

module PID_Controller (
    input wire clk,
    input wire reset,
    input wire [15:0] setpoint,        // 目标值
    input wire [15:0] measured_value,  // 测量值
    output reg [15:0] control_signal   // 控制信号输出
);

parameter Kp = 16'd100; // 比例增益
parameter Ki = 16'd50;  // 积分增益
parameter Kd = 16'd10;  // 微分增益

reg [15:0] integral; 
reg [15:0] prev_error; 
reg [15:0] error; 
reg [15:0] derivative; 

always @(posedge clk or posedge reset) begin
    if (reset) begin
        integral <= 0;
        prev_error <= 0;
        control_signal <= 0;
    end else begin
        error <= setpoint - measured_value;

        integral <= integral + error;

        derivative <= (error - prev_error) * 2; 

        control_signal <= (Kp * error + Ki * integral + Kd * derivative) >> 8; 

        prev_error <= error;
    end
end
endmodule