【PID优化】基于粒子群算法优化PID控制附simulink仿真2

✅作者简介：热爱科研的Matlab仿真开发者，擅长数据处理、建模仿真、程序设计、完整代码获取、论文复现及科研仿真。

🍎 往期回顾关注个人主页：Matlab科研工作室

🍊个人信条：格物致知,完整Matlab代码获取及仿真咨询内容私信。

🔥 内容介绍

基于粒子群算法（PSO）优化 PID 控制，是针对传统 PID 控制器 “参数调试依赖经验、难以兼顾多性能指标” 的高效改进方案。核心逻辑是：将 PID 的比例系数（Kp）、积分系数（Ki）、微分系数（Kd）作为 PSO 的优化变量，通过粒子群的全局搜索能力找到使控制系统 “响应速度快、超调量小、稳态误差小” 的最优参数组合，解决传统试凑法、Ziegler-Nichols 法在复杂系统中调参效率低、控制效果差的问题。

一、核心基础：PID 控制与 PSO 的适配逻辑

⛳️ 运行结果

📣 部分代码

% --- Parámetros y condiciones iniciales (igual que tu código) ---

m = 1.0;          % Masa 

g = 9.81;         % Gravedad 

Ix = 0.1; Iy = 0.1; Iz = 0.2;  % Momentos de inercia 

x0 = [0; 0; 0; 0; 0; 0];       % Posición inicial [x, y, z, roll, pitch, yaw]

xdot0 = [0; 0; 0; 0; 0; 0];    % Velocidades iniciales

X0 = [x0; xdot0];              % Estado completo inicial

% Tiempo de simulación

tspan = [0 20];

% Deseados (mantener constantes para prueba ZN)

z_des = 1; phi_des = 0; theta_des = 0; psi_des = 0;

% --- PASO 1: Encontrar Ganancia Crítica Ku y Periodo de Oscilación Pu para cada canal ---

% Esto se hace incrementando Kp mientras Ki y Kd son 0, hasta que la salida oscile en amplitud constante.

% Función para encontrar Ku y Pu (búsqueda manual o automática)

% Aquí se muestra un ejemplo manual para canal 'z' (altitud). 

% Repite similar para phi, theta y psi.

% Valores iniciales para búsqueda

Ki_test = 0; 

Kd_test = 0;

% Rango de Kp a probar (incremental)

Kp_values = 0:1:50;  % Ajusta según tu sistema

Ku_z = NaN; Pu_z = NaN; % Inicializar

for Kp_test = Kp_values

    % Reiniciar variables globales de integral

    global integral_z integral_phi integral_theta integral_psi;

    integral_z = 0; integral_phi = 0; integral_theta = 0; integral_psi = 0;

    % Simular con Kp_test en z, Ki=0 y Kd=0

    [t, X] = ode45(@(t, X) quadrotor_dynamics(t, X, m, g, Ix, Iy, Iz,...

        Kp_test, Ki_test, Kd_test, 0,0,0, 0,0,0, 0,0,0,...

        z_des, phi_des, theta_des, psi_des), tspan, X0);

    % Extraer señal de salida (altitud)

    z = X(:,3);

    % Revisar si la salida presenta oscilaciones sostenidas

    % Puedes usar análisis simple de oscilaciones:

    % - buscar cruces por cero en la diferencia z-z_des

    z_error = z - z_des;

    zero_crossings = find(diff(sign(z_error)));

    if length(zero_crossings) > 5  % Si tiene más de 5 cruces, puede ser oscilación sostenida

        % Calcular periodo Pu como promedio de diferencias entre cruces

        Pu_z = mean(diff(t(zero_crossings)));

        Ku_z = Kp_test;

        fprintf('Altitud: Ku = %.2f, Pu = %.2f\n', Ku_z, Pu_z);

        break;

    end

end

% Si quieres, repite el mismo proceso para phi, theta y psi (control de actitud)

% --- PASO 2: Calcular parámetros PID con fórmulas Ziegler-Nichols ---

% Fórmulas para controlador PID clásico ZN:

% Kp = 0.6 * Ku

% Ki = 2 * Kp / Pu

% Kd = Kp * Pu / 8

% Para altitud (z)

Kp_z = 0.6 * Ku_z;

Ki_z = 2 * Kp_z / Pu_z;

Kd_z = Kp_z * Pu_z / 8;

% Imprimir resultados

fprintf('Parámetros ZN para altitud (z): Kp=%.2f, Ki=%.2f, Kd=%.2f\n', Kp_z, Ki_z, Kd_z);

% Para actitud (roll, pitch, yaw) se recomienda hacer lo mismo (en este ejemplo uso valores aproximados)

Kp_phi = 0.6 * Ku_z;    Ki_phi = 2 * Kp_phi / Pu_z;    Kd_phi = Kp_phi * Pu_z / 8;

Kp_theta = Kp_phi;      Ki_theta = Ki_phi;             Kd_theta = Kd_phi;

Kp_psi = Kp_phi;        Ki_psi = Ki_phi;               Kd_psi = Kd_phi;

% --- PASO 3: Simular con parámetros PID obtenidos ---

global integral_z integral_phi integral_theta integral_psi;

integral_z = 0; integral_phi = 0; integral_theta = 0; integral_psi = 0;

[t, X] = ode45(@(t, X) quadrotor_dynamics(t, X, m, g, Ix, Iy, Iz,...

    Kp_z, Ki_z, Kd_z, Kp_phi, Ki_phi, Kd_phi,...

    Kp_theta, Ki_theta, Kd_theta, Kp_psi, Ki_psi, Kd_psi,...

    z_des, phi_des, theta_des, psi_des), tspan, X0);

% --- PASO 4: Graficar resultados ---

figure;

subplot(2,2,1);

plot(t, X(:,3), 'b', t, z_des*ones(size(t)), 'r--');

title('Altitud con PID Ziegler-Nichols');

xlabel('Tiempo (s)'); ylabel('Altitud (m)');

legend('Salida','Deseado'); grid on;

subplot(2,2,2);

plot(t, X(:,4), 'b', t, phi_des*ones(size(t)), 'r--');

title('Roll con PID Ziegler-Nichols');

xlabel('Tiempo (s)'); ylabel('Roll (rad)');

legend('Salida','Deseado'); grid on;

subplot(2,2,3);

plot(t, X(:,5), 'b', t, theta_des*ones(size(t)), 'r--');

title('Pitch con PID Ziegler-Nichols');

xlabel('Tiempo (s)'); ylabel('Pitch (rad)');

legend('Salida','Deseado'); grid on;

subplot(2,2,4);

plot(t, X(:,6), 'b', t, psi_des*ones(size(t)), 'r--');

title('Yaw con PID Ziegler-Nichols');

xlabel('Tiempo (s)'); ylabel('Yaw (rad)');

legend('Salida','Deseado'); grid on;

% --- PASO 5: Métricas de desempeño para altitud ---

z = X(:,3);

error_z = z_des - z;

% Overshoot (%)

overshoot = (max(z) - z_des) / z_des * 100;

% MSE (Mean Squared Error)

mse = mean((z - z_des).^2);

% Tiempo de establecimiento (±5%)

tolerance = 0.05 * z_des;

idx_settle = find(abs(z - z_des) > tolerance);

if isempty(idx_settle)

    settling_time = 0;

else

    settling_time = t(find(abs(z - z_des) > tolerance, 1, 'last'));

end

% Tiempo de subida (10% a 90%)

rise_start = z_des * 0.1;

rise_end = z_des * 0.9;

try

    t_rise_start = t(find(z >= rise_start, 1));

    t_rise_end = t(find(z >= rise_end, 1));

    rise_time = t_rise_end - t_rise_start;

catch

    rise_time = NaN; % No se detectó subida completa

end

% Mostrar resultados

fprintf('\n--- MÉTRICAS PARA ALTITUD (z) ---\n');

fprintf('Kp = %.2f, Ki = %.2f, Kd = %.2f\n', Kp_z, Ki_z, Kd_z);

fprintf('Overshoot: %.2f %%\n', overshoot);

fprintf('MSE: %.5f\n', mse);

fprintf('Tiempo de establecimiento: %.2f s\n', settling_time);

fprintf('Tiempo de subida: %.2f s\n', rise_time);

function dXdt = quadrotor_dynamics(t, X, m, g, Ix, Iy, Iz,...

        Kp_z, Ki_z, Kd_z, Kp_phi, Ki_phi, Kd_phi,...

        Kp_theta, Ki_theta, Kd_theta, Kp_psi, Ki_psi, Kd_psi,...

        z_des, phi_des, theta_des, psi_des)

    global integral_z integral_phi integral_theta integral_psi;

    % Extraer estados

    pos = X(1:6);       % [x, y, z, ϕ, θ, ψ]

    vel = X(7:12);      % [dx, dy, dz, dϕ, dθ, dψ]

    % Cálculo de errores

    errores = [z_des - pos(3);    % Error altitud

              phi_des - pos(4);   % Error Roll

              theta_des - pos(5); % Error Pitch

              psi_des - pos(6)];  % Error Yaw

    % Actualizar integrales

    integral_z = integral_z + errores(1);

    integral_phi = integral_phi + errores(2);

    integral_theta = integral_theta + errores(3);

    integral_psi = integral_psi + errores(4);

    % Control PID

    U1 = Kp_z*errores(1) + Ki_z*integral_z + Kd_z*(-vel(3));

    U2 = Kp_phi*errores(2) + Ki_phi*integral_phi + Kd_phi*(-vel(4));

    U3 = Kp_theta*errores(3) + Ki_theta*integral_theta + Kd_theta*(-vel(5));

    U4 = Kp_psi*errores(4) + Ki_psi*integral_psi + Kd_psi*(-vel(6));

    % Dinámica traslacional

    acc_lin = [... 

        (cos(pos(4))*sin(pos(5))*cos(pos(6)) + sin(pos(4))*sin(pos(6)))*U1/m;

        (cos(pos(4))*sin(pos(5))*sin(pos(6)) - sin(pos(4))*cos(pos(6)))*U1/m;

        (cos(pos(4))*cos(pos(5))*U1/m) - g];

    % Dinámica rotacional

    acc_ang = [... 

        (U2 + (Iy - Iz)*vel(5)*vel(6))/Ix;

        (U3 + (Iz - Ix)*vel(4)*vel(6))/Iy;

        (U4 + (Ix - Iy)*vel(4)*vel(5))/Iz];

    % Vector derivadas

    dXdt = [vel; acc_lin; acc_ang];

end

🔗 参考文献

🎈 部分理论引用网络文献，若有侵权联系博主删除

 👇 关注我领取海量matlab电子书和数学建模资料 

🏆团队擅长辅导定制多种科研领域MATLAB仿真，助力科研梦：

🌟 各类智能优化算法改进及应用

生产调度、经济调度、装配线调度、充电优化、车间调度、发车优化、水库调度、三维装箱、物流选址、货位优化、公交排班优化、充电桩布局优化、车间布局优化、集装箱船配载优化、水泵组合优化、解医疗资源分配优化、设施布局优化、可视域基站和无人机选址优化、背包问题、 风电场布局、时隙分配优化、 最佳分布式发电单元分配、多阶段管道维修、 工厂-中心-需求点三级选址问题、 应急生活物质配送中心选址、 基站选址、 道路灯柱布置、 枢纽节点部署、 输电线路台风监测装置、 集装箱调度、 机组优化、 投资优化组合、云服务器组合优化、 天线线性阵列分布优化、CVRP问题、VRPPD问题、多中心VRP问题、多层网络的VRP问题、多中心多车型的VRP问题、 动态VRP问题、双层车辆路径规划（2E-VRP）、充电车辆路径规划（EVRP）、油电混合车辆路径规划、混合流水车间问题、 订单拆分调度问题、 公交车的调度排班优化问题、航班摆渡车辆调度问题、选址路径规划问题、港口调度、港口岸桥调度、停机位分配、机场航班调度、泄漏源定位、冷链、时间窗、多车场等、选址优化、港口岸桥调度优化、交通阻抗、重分配、停机位分配、机场航班调度、通信上传下载分配优化

🌟 机器学习和深度学习时序、回归、分类、聚类和降维

2.1 bp时序、回归预测和分类

2.2 ENS声神经网络时序、回归预测和分类

2.3 SVM/CNN-SVM/LSSVM/RVM支持向量机系列时序、回归预测和分类

2.4 CNN|TCN|GCN卷积神经网络系列时序、回归预测和分类

2.5 ELM/KELM/RELM/DELM极限学习机系列时序、回归预测和分类

2.6 GRU/Bi-GRU/CNN-GRU/CNN-BiGRU门控神经网络时序、回归预测和分类

2.7 ELMAN递归神经网络时序、回归\预测和分类

2.8 LSTM/BiLSTM/CNN-LSTM/CNN-BiLSTM/长短记忆神经网络系列时序、回归预测和分类

2.9 RBF径向基神经网络时序、回归预测和分类

2.10 DBN深度置信网络时序、回归预测和分类

2.11 FNN模糊神经网络时序、回归预测

2.12 RF随机森林时序、回归预测和分类

2.13 BLS宽度学习时序、回归预测和分类

2.14 PNN脉冲神经网络分类

2.15 模糊小波神经网络预测和分类

2.16 时序、回归预测和分类

2.17 时序、回归预测预测和分类

2.18 XGBOOST集成学习时序、回归预测预测和分类

2.19 Transform各类组合时序、回归预测预测和分类

方向涵盖风电预测、光伏预测、电池寿命预测、辐射源识别、交通流预测、负荷预测、股价预测、PM2.5浓度预测、电池健康状态预测、用电量预测、水体光学参数反演、NLOS信号识别、地铁停车精准预测、变压器故障诊断

🌟图像处理方面

图像识别、图像分割、图像检测、图像隐藏、图像配准、图像拼接、图像融合、图像增强、图像压缩感知

🌟 路径规划方面

旅行商问题（TSP）、车辆路径问题（VRP、MVRP、CVRP、VRPTW等）、无人机三维路径规划、无人机协同、无人机编队、机器人路径规划、栅格地图路径规划、多式联运运输问题、 充电车辆路径规划（EVRP）、 双层车辆路径规划（2E-VRP）、 油电混合车辆路径规划、 船舶航迹规划、 全路径规划规划、 仓储巡逻、公交车时间调度、水库调度优化、多式联运优化

🌟 无人机应用方面

无人机路径规划、无人机控制、无人机编队、无人机协同、无人机任务分配、无人机安全通信轨迹在线优化、车辆协同无人机路径规划、

🌟 通信方面

传感器部署优化、通信协议优化、路由优化、目标定位优化、Dv-Hop定位优化、Leach协议优化、WSN覆盖优化、组播优化、RSSI定位优化、水声通信、通信上传下载分配

🌟 信号处理方面

信号识别、信号加密、信号去噪、信号增强、雷达信号处理、信号水印嵌入提取、肌电信号、脑电信号、信号配时优化、心电信号、DOA估计、编码译码、变分模态分解、管道泄漏、滤波器、数字信号处理+传输+分析+去噪、数字信号调制、误码率、信号估计、DTMF、信号检测

🌟电力系统方面

微电网优化、无功优化、配电网重构、储能配置、有序充电、MPPT优化、家庭用电、电/冷/热负荷预测、电力设备故障诊断、电池管理系统（BMS）SOC/SOH估算（粒子滤波/卡尔曼滤波）、 多目标优化在电力系统调度中的应用、光伏MPPT控制算法改进（扰动观察法/电导增量法）、电动汽车充放电优化、微电网日前日内优化、储能优化、家庭用电优化、供应链优化

🌟 元胞自动机方面

交通流 人群疏散 病毒扩散 晶体生长 金属腐蚀

🌟 雷达方面

卡尔曼滤波跟踪、航迹关联、航迹融合、SOC估计、阵列优化、NLOS识别

🌟 车间调度

零等待流水车间调度问题NWFSP 、 置换流水车间调度问题PFSP、 混合流水车间调度问题HFSP 、零空闲流水车间调度问题NIFSP、分布式置换流水车间调度问题 DPFSP、阻塞流水车间调度问题BFSP

👇

5 往期回顾扫扫下方二维码