基于粒子群优化（PSO）的模糊控制器的MATLAB实现
%% 基于粒子群优化的模糊控制器设计
% 作者: DeepSeek
% 创建日期: 2023-10-18
%% 系统初始化
clear; clc; close all;
% 被控系统模型 (二阶系统)
sys = tf(1, [1, 2, 1]); % G(s) = 1/(s^2 + 2s + 1)
% 仿真参数
sim_time = 10; % 仿真时间 (秒)
Ts = 0.01; % 采样时间 (秒)
t = 0:Ts:sim_time;
r = ones(size(t)); % 参考输入 (阶跃信号)
% PSO参数
n_particles = 30; % 粒子数量
max_iter = 50; % 最大迭代次数
w = 0.7; % 惯性权重
c1 = 1.5; % 个体学习因子
c2 = 1.8; % 群体学习因子
v_max = 0.5; % 最大速度
%% 模糊控制器结构定义
% 输入1: 误差 (e)
e_range = [-2, 2]; % 误差范围
e_mf = 3; % 误差的隶属函数数量 (负, 零, 正)
% 输入2: 误差变化率 (de)
de_range = [-1, 1]; % 误差变化率范围
de_mf = 3; % 误差变化率的隶属函数数量 (负, 零, 正)
% 输出: 控制量 (u)
u_range = [-10, 10]; % 控制量范围
u_mf = 5; % 控制量的隶属函数数量
% 计算参数向量维度
% 每个隶属函数有2个参数(中心和宽度) + 规则表参数
n_params = (e_mf * 2) + (de_mf * 2) + (u_mf * 2) + (e_mf * de_mf);
fprintf('模糊控制器参数总数: %d\n', n_params);
%% 粒子群初始化
% 位置和速度初始化
particles = zeros(n_particles, n_params);
velocity = zeros(n_particles, n_params);
% 参数范围 (统一归一化到[0,1])
lb = zeros(1, n_params); % 下界
ub = ones(1, n_params); % 上界
% 随机初始化粒子位置和速度
for i = 1:n_particles
particles(i, :) = rand(1, n_params);
velocity(i, :) = (rand(1, n_params) * 2 - 1) * v_max;
end
% 初始化个体和全局最优
pbest = particles; % 个体最优位置
pbest_fitness = inf(1, n_particles); % 个体最优适应度
gbest = zeros(1, n_params); % 全局最优位置
gbest_fitness = inf; % 全局最优适应度
% 存储收敛历史
fitness_history = zeros(max_iter, 1);
control_history = cell(max_iter, 1);
%% PSO主循环
for iter = 1:max_iter
fprintf('迭代 %d/%d: ', iter, max_iter);
% 评估所有粒子
for i = 1:n_particles
% 从粒子位置解码模糊控制器参数
fis = decode_fis(particles(i, :), e_range, de_range, u_range, e_mf, de_mf, u_mf);
% 仿真系统并计算适应度
fitness = evaluate_fis(fis, sys, r, t, Ts);
% 更新个体最优
if fitness < pbest_fitness(i)
pbest_fitness(i) = fitness;
pbest(i, :) = particles(i, :);
% 更新全局最优
if fitness < gbest_fitness
gbest_fitness = fitness;
gbest = particles(i, :);
best_fis = fis; % 保存最佳模糊控制器
fprintf('新全局最优: ITAE = %.4f | ', gbest_fitness);
end
end
end
% 存储最佳适应度
fitness_history(iter) = gbest_fitness;
control_history{iter} = best_fis; % 保存当前最佳控制器
% 更新粒子速度和位置
for i = 1:n_particles
% 更新速度
r1 = rand(1, n_params);
r2 = rand(1, n_params);
velocity(i, :) = w * velocity(i, :) + ...
c1 * r1 .* (pbest(i, :) - particles(i, :)) + ...
c2 * r2 .* (gbest - particles(i, :));
% 速度限制
velocity(i, :) = min(max(velocity(i, :), -v_max), v_max);
% 更新位置
particles(i, :) = particles(i, :) + velocity(i, :);
% 位置限制
particles(i, :) = min(max(particles(i, :), lb), ub);
end
fprintf('全局最优适应度: %.4f\n', gbest_fitness);
end
%% 结果显示与分析
% 解码最佳模糊控制器
best_fis = decode_fis(gbest, e_range, de_range, u_range, e_mf, de_mf, u_mf);
% 绘制收敛曲线
figure('Name', 'PSO收敛曲线', 'Position', [100, 100, 800, 400]);
plot(1:max_iter, fitness_history, 'b-o', 'LineWidth', 1.5);
xlabel('迭代次数');
ylabel('ITAE性能指标');
title('PSO优化收敛曲线');
grid on;
% 显示最佳模糊控制器
figure('Name', '优化后的模糊控制器', 'Position', [100, 100, 1200, 800]);
subplot(3,1,1);
plotmf(best_fis, 'input', 1);
title('误差(e)隶属函数');
xlabel('误差');
ylabel('隶属度');
subplot(3,1,2);
plotmf(best_fis, 'input', 2);
title('误差变化率(de)隶属函数');
xlabel('误差变化率');
ylabel('隶属度');
subplot(3,1,3);
plotmf(best_fis, 'output', 1);
title('控制量(u)隶属函数');
xlabel('控制量');
ylabel('隶属度');
% 绘制控制曲面
figure('Name', '模糊控制曲面', 'Position', [100, 100, 800, 600]);
gensurf(best_fis);
title('优化后的模糊控制曲面');
xlabel('误差(e)');
ylabel('误差变化率(de)');
zlabel('控制量(u)');
%% 性能对比
% 仿真最佳控制器
[~, y_opt, u_opt] = evaluate_fis(best_fis, sys, r, t, Ts, true);
% 初始随机控制器（用于对比）
init_fis = decode_fis(particles(1, :), e_range, de_range, u_range, e_mf, de_mf, u_mf);
[~, y_init, u_init] = evaluate_fis(init_fis, sys, r, t, Ts, true);
% 传统PID控制器（用于对比）
Kp = 1.5; Ki = 0.8; Kd = 0.2;
[y_pid, u_pid] = pid_controller(sys, r, t, Ts, Kp, Ki, Kd);
% 绘制响应曲线
figure('Name', '系统响应对比', 'Position', [100, 100, 1000, 800]);
% 系统输出
subplot(2,1,1);
plot(t, r, 'k--', 'LineWidth', 1.5); hold on;
plot(t, y_init, 'b-', 'LineWidth', 1.5);
plot(t, y_pid, 'g-', 'LineWidth', 1.5);
plot(t, y_opt, 'r-', 'LineWidth', 1.5);
xlabel('时间 (秒)');
ylabel('系统输出');
title('系统响应对比');
legend('参考输入', '初始模糊控制', 'PID控制', '优化模糊控制', 'Location', 'southeast');
grid on;
% 控制信号
subplot(2,1,2);
plot(t, u_init, 'b-', 'LineWidth', 1.5); hold on;
plot(t, u_pid, 'g-', 'LineWidth', 1.5);
plot(t, u_opt, 'r-', 'LineWidth', 1.5);
xlabel('时间 (秒)');
ylabel('控制信号');
title('控制信号对比');
legend('初始模糊控制', 'PID控制', '优化模糊控制', 'Location', 'northeast');
grid on;
%% 性能指标计算
metrics = struct();
metrics.initial = calculate_metrics(y_init, r, t, u_init);
metrics.pid = calculate_metrics(y_pid, r, t, u_pid);
metrics.optimized = calculate_metrics(y_opt, r, t, u_opt);
% 显示性能指标
disp(' ');
disp('============== 性能指标对比 ==============');
fprintf('%-20s | %-8s | %-8s | %-8s | %-8s\n', '控制器', 'ITAE', 'ISE', 'IAE', '控制能量');
fprintf('%-20s | %-8.4f | %-8.4f | %-8.4f | %-8.4f\n', ...
'初始模糊控制', metrics.initial.itae, metrics.initial.ise, metrics.initial.iae, metrics.initial.energy);
fprintf('%-20s | %-8.4f | %-8.4f | %-8.4f | %-8.4f\n', ...
'PID控制', metrics.pid.itae, metrics.pid.ise, metrics.pid.iae, metrics.pid.energy);
fprintf('%-20s | %-8.4f | %-8.4f | %-8.4f | %-8.4f\n', ...
'优化模糊控制', metrics.optimized.itae, metrics.optimized.ise, metrics.optimized.iae, metrics.optimized.energy);
%% ================== 辅助函数 ==================
%% 从粒子位置解码模糊控制器
function fis = decode_fis(position, e_range, de_range, u_range, e_mf, de_mf, u_mf)
% 创建初始FIS
fis = newfis('pso_fis', 'mamdani', 'min', 'max', 'min', 'max', 'centroid');
% 添加输入变量: 误差(e)
fis = addvar(fis, 'input', 'e', e_range);
% 添加输入变量: 误差变化率(de)
fis = addvar(fis, 'input', 'de', de_range);
% 添加输出变量: 控制量(u)
fis = addvar(fis, 'output', 'u', u_range);
% 参数索引
idx = 1;
% 添加误差(e)的隶属函数
for i = 1:e_mf
center = position(idx) * diff(e_range) + e_range(1);
width = position(idx+1) * 0.5 * diff(e_range) + 0.05 * diff(e_range);
fis = addmf(fis, 'input', 1, sprintf('mf%d', i), 'gaussmf', [width, center]);
idx = idx + 2;
end
% 添加误差变化率(de)的隶属函数
for i = 1:de_mf
center = position(idx) * diff(de_range) + de_range(1);
width = position(idx+1) * 0.5 * diff(de_range) + 0.05 * diff(de_range);
fis = addmf(fis, 'input', 2, sprintf('mf%d', i), 'gaussmf', [width, center]);
idx = idx + 2;
end
% 添加控制量(u)的隶属函数
for i = 1:u_mf
center = position(idx) * diff(u_range) + u_range(1);
width = position(idx+1) * 0.5 * diff(u_range) + 0.05 * diff(u_range);
fis = addmf(fis, 'output', 1, sprintf('mf%d', i), 'gaussmf', [width, center]);
idx = idx + 2;
end
% 添加规则表 (e_mf × de_mf 条规则)
rules = zeros(e_mf * de_mf, 3);
for i = 1:e_mf
for j = 1:de_mf
% 规则输出索引 (1到u_mf)
output_mf = round(position(idx) * (u_mf - 1)) + 1;
rules((i-1)*de_mf + j, :) = [i, j, output_mf, 1, 1];
idx = idx + 1;
end
end
fis = addrule(fis, rules);
end
%% 评估模糊控制器性能
function [fitness, y, u] = evaluate_fis(fis, sys, r, t, Ts, return_output)
if nargin < 6, return_output = false; end
% 离散化系统
sys_d = c2d(sys, Ts);
% 初始化变量
y = zeros(size(t));
u = zeros(size(t));
e = zeros(size(t));
de = zeros(size(t));
% 初始条件
y(1) = 0;
e(1) = r(1) - y(1);
% 仿真循环
for k = 2:length(t)
% 计算误差变化率
de(k) = (e(k-1) - (r(k-1) - y(k-1))) / Ts;
% 模糊推理
u(k) = evalfis(fis, [e(k-1), de(k)]);
% 限制控制信号
u(k) = min(max(u(k), -10), 10);
% 系统响应
y(k) = sys_d.num{1}*[u(k); u(k-1)] - sys_d.den{1}(2)*y(k-1) - sys_d.den{1}(3)*y(k-2);
% 更新误差
e(k) = r(k) - y(k);
end
% 计算适应度 (ITAE = ∫t|e|dt)
itae = sum(t' .* abs(e)) * Ts;
if return_output
fitness = itae;
else
fitness = itae;
y = [];
u = [];
end
end
%% PID控制器仿真
function [y, u] = pid_controller(sys, r, t, Ts, Kp, Ki, Kd)
% 离散化系统
sys_d = c2d(sys, Ts);
% 初始化变量
y = zeros(size(t));
u = zeros(size(t));
e = zeros(size(t));
ie = 0; % 积分项
de = 0; % 微分项
e_prev = 0;
% 仿真循环
for k = 2:length(t)
% 计算误差
e(k) = r(k) - y(k-1);
% PID控制
ie = ie + e(k) * Ts;
de = (e(k) - e_prev) / Ts;
u(k) = Kp * e(k) + Ki * ie + Kd * de;
% 限制控制信号
u(k) = min(max(u(k), -10), 10);
% 系统响应
y(k) = sys_d.num{1}*[u(k); u(k-1)] - sys_d.den{1}(2)*y(k-1) - sys_d.den{1}(3)*y(k-2);
% 更新误差
e_prev = e(k);
end
end
%% 计算性能指标
function metrics = calculate_metrics(y, r, t, u)
% 误差计算
e = r - y;
% 性能指标
metrics.itae = sum(t' .* abs(e)) * (t(2)-t(1)); % ITAE
metrics.ise = sum(e.^2) * (t(2)-t(1)); % ISE
metrics.iae = sum(abs(e)) * (t(2)-t(1)); % IAE
metrics.energy = sum(u.^2) * (t(2)-t(1)); % 控制能量
end
算法原理与实现解析
1. 模糊控制器结构
本实现设计了一个二维模糊控制器：
- 输入1：误差 (e = 参考输入 - 系统输出)
- 输入2：误差变化率 (de = de/dt)
- 输出：控制量 (u)
% 模糊控制器结构参数
e_range = [-2, 2]; % 误差范围
e_mf = 3; % 误差的隶属函数数量 (负, 零, 正)
de_range = [-1, 1]; % 误差变化率范围
de_mf = 3; % 误差变化率的隶属函数数量
u_range = [-10, 10]; % 控制量范围
u_mf = 5; % 控制量的隶属函数数量
2. 粒子群优化参数编码
PSO算法优化以下模糊控制器参数：
- 输入/输出变量的隶属函数参数（中心位置和宽度）
- 模糊规则表的输出隶属函数索引
% 参数向量结构:
% 前(e_mf*2)个参数: 误差e的隶属函数参数 [中心1, 宽度1, 中心2, 宽度2, ...]
% 接着(de_mf*2)个参数: 误差变化率de的隶属函数参数
% 接着(u_mf*2)个参数: 控制量u的隶属函数参数
% 最后(e_mf*de_mf)个参数: 规则表输出索引
3. PSO优化过程
for iter = 1:max_iter
% 评估所有粒子
for i = 1:n_particles
% 解码模糊控制器
fis = decode_fis(particles(i, :), ...);
% 评估控制器性能
fitness = evaluate_fis(fis, sys, r, t, Ts);
% 更新个体最优和全局最优
if fitness < pbest_fitness(i)
pbest_fitness(i) = fitness;
pbest(i, :) = particles(i, :);
if fitness < gbest_fitness
gbest_fitness = fitness;
gbest = particles(i, :);
end
end
end
% 更新粒子速度和位置
for i = 1:n_particles
% 速度更新公式
velocity(i, :) = w * velocity(i, :) + ...
c1 * rand() * (pbest(i, :) - particles(i, :)) + ...
c2 * rand() * (gbest - particles(i, :));
% 位置更新
particles(i, :) = particles(i, :) + velocity(i, :);
end
end
4. 性能评估指标
使用ITAE（积分时间绝对误差）作为主要性能指标：
% ITAE = ∫t|e|dt
itae = sum(t' .* abs(e)) * Ts;
同时计算其他性能指标用于对比：
- ISE（积分平方误差）
- IAE（积分绝对误差）
- 控制能量消耗
参考代码 基于粒子群优化的模糊控制器 youwenfan.com/contentalb/46317.html
本实现提供了完整的PSO优化模糊控制器框架，用户可根据具体被控对象调整系统模型和控制器参数，实现高性能的自适应控制。
