电机控制算法综述01-2-模糊PID控制器的Matlab仿真

文章目录

1. 1. 模糊PID控制的matlab脚本实现
   1. 1.1. 电机模型的搭建
   2. 1.2. PID控制器设计
   3. 1.3. 模糊PID控制器控制原则
   4. 1.4. 模糊化
   5. 1.5. 模糊推理
   6. 1.6. 解模糊
2. 2. 模糊PID控制的Simulink仿真模型搭建

这一专栏将介绍电机控制算法，为接下来的课题研究提供参考。
本文首先采用DeepSeek生成一个模糊PID控制电机转速的matlab脚本代码，结合代码介绍模糊PID控制原理。
然后介绍Simulink自带模糊控制工具箱(Fuzzy Logic Toolbox)的使用方法，采用Simulink搭建采用模糊PID控制电机转速的仿真模型。（未完待续）
最后对采用模糊PID前后的控制效果进行对比，另外对纯matlab脚本与Simulink仿真实现的控制效果进行对比。（未完待续）

模糊PID控制的matlab脚本实现

DeepSeek生成代码如下：

% 主程序：模糊PID控制的电机转速环控制
%clear; clc;

% 电机参数
J = 0.01;   % 转动惯量
B = 0.1;    % 阻尼系数
Kt = 0.15;  % 转矩常数
Ts = 0.001;  % 采样时间
T = 0:Ts:0.1; % 仿真时间

% 初始状态
speed = 0;
target_speed = 100; % 目标转速

% PID初始参数
Kp0 = 100;
Ki0 = 40;
Kd0 = 0.001;

% 初始化误差记录
e_prev = 0;
integral = 0;

% 存储结果
speed_history = zeros(size(T));
u_history = zeros(size(T));
kp_history = zeros(size(T));
ki_history = zeros(size(T));
kd_history = zeros(size(T));

% 初始化模糊控制器
fuzzy_params = init_fuzzy();

% 主控制循环
for k = 1:length(T)
    % 计算误差
    e = target_speed - speed;
    
    % 计算误差变化率
    if k == 1
        ec = 0;
    else
        ec = (e - e_prev) / Ts;
    end
    
    % 模糊推理调整参数
    [delta_Kp, delta_Ki, delta_Kd] = fuzzy_pid(e, ec, fuzzy_params);
    
    % 更新PID参数（并进行限幅）
    Kp = max(min(Kp0 + delta_Kp, 1000), 0);
    Ki = max(min(Ki0 + delta_Ki, 1000), 0);
    Kd = max(min(Kd0 + delta_Kd, 1000), 0);
    
    % 计算控制量
    integral = integral + e * Ts;
    derivative = ec;
    u = Kp * e + Ki * integral + Kd * derivative;
    
    % 执行限幅
    u = max(min(u, 220), -220);
    
    % 更新系统状态（简化电机模型）
    speed = (1 - B*Ts/J)*speed + (Kt*Ts/J)*u;
    
    % 保存数据
    speed_history(k) = speed;
    u_history(k) = u;
    kp_history(k) = Kp;
    ki_history(k) = Ki;
    kd_history(k) = Kd;
    e_prev = e;
end

% 绘制结果
figure(1);
subplot(2,1,1);
plot(T, speed_history, 'b', T, target_speed*ones(size(T)), 'r--');
title('转速响应'); xlabel('时间(s)'); ylabel('转速(rpm)'); legend('实际', '目标');
subplot(2,1,2);
plot(T, u_history, 'g');
title('控制量'); xlabel('时间(s)'); ylabel('控制电压(V)');
figure(2);
subplot(3,1,1);
plot(T, kp_history);
title('PID参数变化');
xlabel('时间(s)'); ylabel('Kp');
subplot(3,1,2);
plot(T, ki_history);
xlabel('时间(s)'); ylabel('Ki');
subplot(3,1,3);
plot(T, kd_history);
xlabel('时间(s)'); ylabel('Kd');

% 初始化模糊控制器函数
function params = init_fuzzy()
    % 输入变量隶属度函数参数
    params.e_mf = [       % 误差e的模糊集参数
        [-30, -30, -20]  % 1: NB
        [-30, -20, -10]    % 2: NM
        [-20, -10, 0]    % 3: NS
        [-10, 0, 10]     % 4: Z
        [0, 10, 20]      % 5: PS
        [10, 20, 30]      % 6: PM
        [20, 30, 30]     % 7: PB
    ];
    
    params.ec_mf = [      % 误差变化率ec的模糊集参数
        [-30, -30, -20]  % 1: NB
        [-30, -20, -10]    % 2: NM
        [-20, -10, 0]    % 3: NS
        [-10, 0, 10]     % 4: Z
        [0, 10, 20]      % 5: PS
        [10, 20, 30]      % 6: PM
        [20, 30, 30]     % 7: PB
    ];
    
    % 输出变量单点值
    params.dKp = [-0.3, -0.2, -0.1, 0, 0.1, 0.2, 0.3];    % ΔKp的模糊集值
    params.dKi = [-1.5, -1, -0.5, 0, 0.5, 1, 1.5];% ΔKi的模糊集值
    params.dKd = [-0.0003, -0.0002, -0.0001, 0, 0.0001, 0.0002, 0.0003];% ΔKd的模糊集值
    
    % 模糊规则表[E, EC, dKp, dKi, dKd]
    params.rules = [
    1 1 7 1 5;1 2 7 1 3;1 3 6 2 1;1 4 6 2 1;1 5 5 3 1;1 6 4 4 2;1 7 4 4 2;
    2 1 7 1 5;2 2 7 1 3;2 3 6 3 1;2 4 5 4 2;2 5 5 4 2;2 6 4 4 3;2 7 3 4 3;
    3 1 6 1 4;3 2 6 2 3;3 3 6 3 2;3 4 5 4 3;3 5 4 5 3;3 6 3 5 3;3 7 3 5 3;
    4 1 6 3 4;4 2 6 3 3;4 3 5 4 3;4 4 4 4 3;4 5 3 5 3;4 6 2 6 3;4 7 2 6 3;
    5 1 5 3 4;5 2 5 4 4;5 3 4 5 4;5 4 3 5 4;5 5 3 6 4;5 6 2 6 4;5 7 2 6 4;
    6 1 4 4 7;6 2 4 4 6;6 3 3 5 5;6 4 3 6 5;6 5 2 6 5;6 6 2 7 5;6 7 1 7 5;
    7 1 4 4 7;7 2 4 4 6;7 3 3 5 6;7 4 3 6 6;7 5 2 6 6;7 6 1 7 7;7 7 1 7 7;
    ];
end

% 模糊PID控制器函数
function [dKp, dKi, dKd] = fuzzy_pid(e, ec, params)
    % 模糊化
    E=-e;%误差e定义为r-y，但模糊控制表中E定义为y-r
    e_deg = calc_membership(E, params.e_mf);
    ec_deg = calc_membership(ec, params.ec_mf);    
    % 推理和解模糊
    [dKp, dKi, dKd] = defuzzify(e_deg, ec_deg, params);
end

% 隶属度计算函数
function degrees = calc_membership(x, mf_params)
    degrees = zeros(1, size(mf_params,1));
    for i = 1:size(mf_params,1)
        a = mf_params(i,1);
        b = mf_params(i,2);
        c = mf_params(i,3);
        
        if x <= a
            degrees(i) = 0;
        elseif x >= c
            degrees(i) = 0;
        elseif x >= a && x <= b
            if a == b
                degrees(i) = 1;
            else
                degrees(i) = (x - a)/(b - a);
            end
        else
            if b == c
                degrees(i) = 1;
            else
                degrees(i) = (c - x)/(c - b);
            end
        end
    end
end

% 解模糊函数
function [dKp, dKi, dKd] = defuzzify(e_deg, ec_deg, params)
    sum_w_kp = 0; sum_kp = 0;
    sum_w_ki = 0; sum_ki = 0;
    sum_w_kd = 0; sum_kd = 0;
    
    for i = 1:size(params.rules,1)
        rule = params.rules(i,:);
        e_idx = rule(1);
        ec_idx = rule(2);
        kp_idx = rule(3);
        ki_idx = rule(4);
        kd_idx = rule(5);
        
        % 规则激活强度
        strength = min(e_deg(e_idx), ec_deg(ec_idx));
        
        % 加权累加
        sum_kp = sum_kp + strength * params.dKp(kp_idx);
        sum_w_kp = sum_w_kp + strength;
        
        sum_ki = sum_ki + strength * params.dKi(ki_idx);
        sum_w_ki = sum_w_ki + strength;
        
        sum_kd = sum_kd + strength * params.dKd(kd_idx);
        sum_w_kd = sum_w_kd + strength;
    end
    
    % 计算解模糊值
    dKp = sum_kp / (sum_w_kp + eps);
    dKi = sum_ki / (sum_w_ki + eps);
    dKd = sum_kd / (sum_w_kd + eps);

end

电机模型的搭建

对于各种类型的旋转动力系统，其动力学方程式（机电能量转换方程）为：
\(T_e-T_L=J\frac{d^2 \theta}{dt^2}+B\frac{d\theta}{dt}=J\frac{d \omega}{dt}+B\omega\)。
式中转速\(\omega=\frac{d\theta}{dt}\)，\(T_e\)是电磁转矩，\(T_L\)是负载转矩，\(J\)为转动惯量，\(B\)为粘性摩擦系数（有的文献也用\(D\)表示）。
注：没错，这个公式不仅仅适用于电机，将电磁转矩\(T_e\)换成原动力转矩后，就可以适用于任何绕轴旋转的系统。对于电机而言原动力自然是电；如果这个转轴是由蒸汽驱动的，那么原动力就是蒸汽产生的热量；这个轴当然也可以是人力或畜力驱动的。
本文将该方程式转换为差分方程形式，其代码为：
speed = (1 - B*Ts/J)*speed + (Kt*Ts/J)*u;
将其整理成动力学方程式的样式如下：
\(K_t \cdot u=J\frac{speed(t+1)-speed(t)}{T_s}+B \cdot speed(t)\)
式中Ts为采样时间（步长），是仿真运行过程中两个时刻的间隔时间，单位为秒。
上文代码中原采样时间为0.01s，实际上这个采样时间太大了，为了使得仿真运行更精确，我将其修改为0.001s，也就是1ms，可以满足对算法可行性进行初步试验的需要。而对于实际的电机控制系统，采样时间通常取\(10^{-5}s\)(10μs)或\(10^{-6}s\)(1μs)，这样的仿真结果更为精确，可以满足毕业论文的需要。
u是模糊PID控制器输出量，对于电机控制系统这个量就是电压。
对于不同的电机，电压到转矩的转换过程各有不同，这取决于电机的动态特性（与电枢电阻、电感、磁链等相关）以及功率变换器件工作特性，但通常电压与转矩呈正比关系。
本文对这个转换过程进行简化，取\(T_e-T_L=K_t \cdot u\)，式中\(K_t\)是转矩常数，也就是电压转换到转矩的比例系数。
当前时刻转速为\(speed(t)\)，位于代码等号右侧，下一时刻转速为\(speed(t+1)\)，位于代码等号左侧。
由此得到的转速变化量\(\frac{d \omega}{dt}=\frac{speed(t+1)-speed(t)}{T_s}\)

PID控制器设计

PID控制器输出公式如下：
u = Kp * e + Ki * integral + Kd * derivative;
\(u=K_p \cdot e + K_i \cdot \int e + K_d \cdot ec\)
e是转速误差，即目标转速减去实际转速所得的差；
integral是误差累计和，程序中的初始值为0，其计算公式为\(integral(t+1) = integral(t) + e \cdot T_s\) derivative是误差变化率(ec)，其计算公式为\(ec=(e-e_{prev})/T_s\)。
式中\(e_{prev}\)是上一时刻的误差，该值在程序中的初始值为0。
但这里就会出现一个问题，若电机目标转速设定为100转，显然t=0时刻（程序中k=1时）实际转速为0，得到e=100，那么计算得到的ec值就会极大，这显然不符合实际情况，因此设定k=1时ec=0。

模糊PID控制器控制原则

模糊PID控制器的核心在于根据控制过程的不同阶段（启动阶段、动态调节阶段、稳态阶段）动态特性，通过误差（e）和误差变化率（ec）的模糊关系，动态调整PID参数（Kp,Ki,Kd），尽可能兼顾系统动态与稳态性能（响应速度、超调量、稳态误差等）要求。
对于不同的被控对象，由于其动态特性有所不同，因此得到的控制规则也会不同。因此需要根据控制系统实际情况修正控制规则，下面给出的规则仅做参考。

起动阶段：此时电机转速误差为较大正值(e=PB)，误差变化量为较大负值(ec=NB)。
此时需要增大\(Kp\)使其具有较高的响应速度，限制\(K_i\)避免积分饱和，适当增大\(K_d\)抑制超调。
（\(K_p\)要大，\(K_i\)要小，\(K_d\)较大）
调节阶段：此时电机转速误差与误差变化量绝对值不大不小(e=PM,NM;ec=NM,PM)。
此时需要减小\(K_p\)避免超调，逐步增大\(K_i\)以消除稳态误差，保持\(K_d\)不变实现稳定过渡。
（\(K_p\)要小，\(K_i\)适中，\(K_d\)不变）
稳定阶段：此时电机转速误差与误差变化量绝对值均很小(e=PS,ZO,NS;ec=NS,ZO,PS)。
此时需要维持较小\(K_p\)防止振荡，增大\(K_i\)消除残余误差，适当增大\(K_d\)增强抗干扰能力。
（\(K_p\)要小，\(K_i\)要大，\(K_d\)要大）
在另一篇文章中，要求\(K_d\)取较小值，减小被控过程的制动作用，进而补偿在调节过程初期由于\(K_d\)值较大所造成的调节过程的时间延长

模糊化

将精确值(e,ec)转化为模糊值(E,EC)的过程被称作模糊化。
论域：系统输入的误差e与误差变化率ec、输出的PID参数变化量\(\Delta K_p,\Delta K_i,\Delta K_d\)的取值范围被称作论域。
模糊集：模糊量E和EC的所有可能取值的集合是模糊集。本文代码中模糊集为{NB、NM、NS、ZO、PS、PM、PB}，依次代表负大(Negative Big)、负中(Negative Medium)、负小(Negative Small)、零(Zero)、正小(Positive Small)、正中(Positive Medium)、正大(Positive Big)。
量化因子：论域中的精确值转换到模糊集中的模糊值时的比例值为量化因子。
精确值x的论域是{x_L,x_H}，模糊集为{-n,-n+1,-n+2,...,-1,0,1,...,n-1,n,n+1}，则量化因子\(k_x=\frac{2n}{x_H-x_L}\)
文中代码误差e的论域为{-30,30}，模糊集{NB、NM、NS、ZO、PS、PM、PB}可以记作{-3,-2,-1,0,1,2,3}，n=3，求得量化因子\(k_e=\frac{1}{10}\)。
可以得到精确值到模糊值的转换公式\(y=k_x(x-\frac{x_H+x_L}{2})\)
比如若误差e的值为20，\(e_L=-30,e_H=30,k_e=10\)，则可计算出\(y=\frac{1}{10} \times (20-\frac{30+(-30)}{2})=2\)，其对应的模糊值为PM。
而若误差e的值为23，则计算得y=1.3，该值在1与2之间，所以对应模糊值为PM与PB。
比例因子：模糊值转换到输出量精确值的比例关系为比例因子。
\(k_u=\frac{2n}{u_H-u_L}\)，值得注意的是，这里的u是输出的PID参数变化量\(\Delta K_p,\Delta K_i,\Delta K_d\)，不是模糊PID控制器最后的输出量。
隶属度与隶属度函数：一般精确值不太可能对应一个模糊值，很可能对应多个模糊值，比如上文误差e=23时，对应的模糊值为PM与PB。那么这个精确值对应于PM的程度高还是对应于PB的程度更高，这就需要引入隶属度与隶属度函数的概念。
最常用的隶属度函数形式为三角形，这里就以三角形为例。
图1 三角形隶属度函数

隶属度\(\mu(x)=\left\{ \begin{gathered} 0 & x\le a 或 x\ge b \\ \frac{x-a}{b-a} & a \le x \le b \\ \frac{c-x}{c-b} & b \le x \le c \\ 1 & (a=b,x=a)或(b=c,x=c) \end{gathered}\right.\)
根据上述公式，可以计算出误差e=23时，PM的隶属度函数为0.7，PB的隶属度函数为0.3，NB、NM、NS、ZO、PS的隶属度函数值为0。

模糊推理

模糊推理是根据模糊控制表，将模糊值E与EC与PID参数变化量\(\Delta K_p,\Delta K_i,\Delta K_d\)的模糊值相对应。
本文程序中模糊控制表储存于params.rules中，其中第一行1 1 7 1 5表示当E=NB、EC=NB时，dKp=PB、dKi=NB、dKd=PS。
模糊控制表的建立基于长期积累所得的控制经验，本文所采用的模糊控制表是网络上流传最广的一种。
对于某个特定条件下的控制，需要结合被控对象的实际情况及电机控制需求对该规则表加以调整。
图2 ΔKp、ΔKi、ΔKd模糊规则表

上述模糊规则表目前存在一个比较引人注目的争议，该争议涉及模糊值E的定义问题。在绝大多数的文献与仿真模型中，误差e被定义为设定值减去实际值所得之差，即e=r-y，按此道理应当将r-y所得的e模糊化得到E。
但也有人认为，该模糊表中的E是由y-r所得的e模糊化而来，并且若定义e=y-r，则该模糊表更契合既定控制要求。
（可以设计实验来对比两种定义下的控制效果，实际控制效果比误差定义更重要）
参考文章：
CSDN-模糊PID控制的规则表一点理解
知乎-怎么理解模糊pid控制表？

解模糊

解模糊是模糊控制流程的最后一部，即将PID参数变化量\(\Delta K_p,\Delta K_i,\Delta K_d\)的模糊值转换为相应的精确值。
解模糊存在多种算法，其中包括重心法（Centroid Method，面积中心法）、最大值平均法（Mean of Maxima, MoM）、中心平均法（Center Average）等，更多细节可以参考这篇文章：
CSDN-推理过程的去模糊化（Defuzzification）方法详解
其中重心法是最为常用的解模糊算法，它取的是加权平均值，其公式如下：
\[\Delta K=\frac{\sum_{i=1}^n \mu_i c_i}{\sum_{i=1}^n \mu_i}\] 该公式的输出\(\Delta K\)为PID参数变化量\(\Delta K_p,\Delta K_i,\Delta K_d\)。
式中n是模糊规则的条数（上文代码中有49条规则），\(\mu_i\)为第i条规则的激活强度，\(c_i\)是对应的模糊集值（存储于params.dKp,params.dKi,params.dKd）。
激活强度\(\mu_i=\min(\mu_e,\mu_{ec})\)，也就是各个模糊集值的最小值。这体现了模糊逻辑中的最小化原则，该原则可避免高激活强度规则掩盖弱规则的影响，确保控制量调整更稳健。

最后，得到下一时刻的PID控制器参数值为：
\[K(n)=K(n-1)+\Delta K\]

模糊PID控制的Simulink仿真模型搭建

To Be Countinued未完待续......