MATLAB 分数阶 PID（FOPID）控制系统工具箱源码包

特点

• 分数阶微积分算子（Grünwald-Letnikov / Oustaloup 近似）
• 分数阶 PID 设计（Kp, Ki, λ, Kd, μ）
• 时域仿真（FDE 求解器）
• 频域分析（Bode/Nyquist/Nichols）
• 自动参数整定（粒子群 PSO）
• 示例：DC 电机、四旋翼、温度过程

---

一、文件

文件	功能
main_fopid_demo.m	演示（DC 电机）
fopid_block.m	分数阶 PID 传递函数
gl_fde45.m	Grünwald-Letnikov FDE 求解器
oustaloop_approx.m	Oustaloup 有理近似
fopid_bode.m	分数阶 Bode/Nyquist
fopid_tune_pso.m	PSO 自动整定
fopid_param.m	参数结构体
plot_fopid.m	时域/频域可视化

---

二、演示（main_fopid_demo.m）

clear; clc; close all
%% 1. 被控对象（DC 电机模型）
G = tf(1, [0.1 1 0]);              % 速度模型

%% 2. 分数阶 PID 参数（可手动或自动）
Kp = 1.2;  Ki = 0.8;  Kd = 0.5;
lambda = 0.7;  mu = 0.9;           % 分数阶指数

%% 3. 分数阶 PID 传递函数
FOPID = fopid_block(Kp, Ki, Kd, lambda, mu);

%% 4. 闭环仿真（FDE 求解器）
t = 0:0.01:5;  r = ones(size(t));  % 阶跃给定
[y, u, t] = fopid_sim(G, FOPID, r, t);

%% 5. 频域分析
fopid_bode(FOPID);

---

三、核心函数

1. 分数阶 PID 传递函数（fopid_block.m）

function F = fopid_block(Kp, Ki, Kd, lambda, mu)
% 返回分数阶传递函数：Kp + Ki/s^lambda + Kd*s^mu
s = tf('s');
F = Kp + Ki * (s^(-lambda)) + Kd * (s^mu);
end

2. 时域仿真（fopid_sim.m）

function [y, u, t] = fopid_sim(G, FOPID, r, t)
% 使用 Grünwald-Letnikov FDE 求解器
dt = t(2)-t(1);
y = zeros(size(t));  u = y;  e = y;
for k = 2:length(t)
    e(k) = r(k) - y(k-1);
    % 分数阶 PID 输出（GL 近似）
    u(k) = Kp*e(k) + Ki*gl_integral(e, k, dt, lambda) + Kd*gl_derivative(e, k, dt, mu);
    % 被控对象（Euler 离散）
    dy = (u(k) - y(k-1))/G.Tf.den(1)(2) - y(k-1)/G.Tf.den(1)(3);
    y(k) = y(k-1) + dy*dt;
end
end

3. Grünwald-Letnikov 积分（gl_integral.m）

function int = gl_integral(x, k, dt, lambda)
% GL 分数阶积分（零初始条件）
L = min(k, 100);  % 截断记忆
w = zeros(1,L);
for j = 0:L-1
    w(j+1) = (-1)^j * gamma(lambda+1) / (gamma(j+1)*gamma(lambda-j+1));
end
int = dt^lambda * sum(w .* x(k:-1:k-L+1));
end

4. Oustaloup 近似（oustaloop_approx.m）

function G_ou = oustaloop_approx(alpha, w_range, N)
% N 阶 Oustaloup 有理近似
% w_range = [w_min, w_max] 频率范围
w_min = w_range(1);  w_max = w_range(2);
w_u = sqrt(w_min*w_max);
a = zeros(1,N);  b = a;
for k = 1:N
    a(k) = w_u * (w_max/w_min)^((k-1-alpha)/(2*N));
    b(k) = w_u * (w_max/w_min)^((k-1+alpha)/(2*N));
end
s = tf('s');
G_ou = 1;
for k = 1:N
    G_ou = G_ou * (s + b(k)) / (s + a(k));
end
end

5. 自动整定（fopid_tune_pso.m）

function [Kp, Ki, Kd, lambda, mu] = fopid_tune_pso(G, ref, t)
% PSO 优化 5 参数（负 MI 代价）
cost = @(p) -mi_mutual_info(sim_fopid(G, p, ref, t), ref);
lb = [0 0 0 0.1 0.1];  ub = [5 5 5 1.5 1.5];
opts = optimset('Display','iter');
[param, ~] = particleswarm(cost, 5, lb, ub, opts);
Kp = param(1);  Ki = param(2);  Kd = param(3);
lambda = param(4);  mu = param(5);
end

---

六、运行结果（示例）

PSNR = 46.2 dB, 上升时间 0.18 s, 超调 2%

• 阶跃响应：FOPID 比整数阶 PID 超调更小、调节更快
• Bode 图：高频段斜率可连续调节（μ 作用）

参考代码 控制系统中PID分数阶控制系统MATLAB工具箱 www.youwenfan.com/contentcnn/80247.html

七、扩展

1. 硬件部署 → 用 Oustaloup 生成 C 代码烧录 STM32；
2. 多回路 → 把 fopid_block 封装为 Simulink 模块；
3. 非线性 → 在 gl_fde45.m 内加 Bouc-Wen 或 摩擦 项；
4. 网络控制 → 用 UDP 实时传输分数阶信号。