3.2.2 粒子群算法的改进策略 针对咖啡机温度控制系统的非线性、时变特性与实时性需求，本文对标准粒子群算法（PSO）进行了两方面改进，以提升其全局搜索能力与收敛效率： 非线性动态惯性权重 为平衡算法全局探索与局部开发能力，采用指数递减惯性权重策略： w(k)=w_min+(w_max-w_min)⋅e^(-10k/K_max ) 式中，k为当前迭代次数，K_max为最大迭代次数（设为100），w_max=0.9、w_min=0.4。初期权重较高（接近0.9），增强全局搜索；后期快速衰减至0.4，加速局部收敛。 2.自适应学习因子 根据迭代进度动态调整个体学习因子c_1与群体学习因子c_2： c_1 (k)=c_(1,max)-((c_(1,max)-c_(1,min))⋅k)/K_max c_2 (k)=c_(2,min)+((c_(2,max)-c_(2,min))⋅k)/K_max 初始阶段，c_1=2.5、c_2=1.5，侧重个体经验；随着迭代进行，c_1线性降至1.5，c_2升至2.5，逐步转向群体协作。 3.2.3 适应度函数设计与多目标优化 性能指标量化 针对咖啡机温度控制需求，定义四项关键性能指标： 超调量（σ%）：反映系统阻尼特性，要求σ%≤5%； 调节时间（t_s）：系统进入±1℃误差带的时间，目标t_s≤60s； 稳态误差（e_ss）：稳态温度偏差绝对值，要求e_ss≤0.5℃； 能耗积分(∫|u(t)|dt)：控制量绝对值的积分，衡量能效。 适应度函数构建 采用加权求和法将多目标转化为单目标优化问题： J=α⋅σ%+β⋅t_s+γ⋅e_ss+δ⋅∫|u(t)|dt 权重系数通过熵权法确定：计算各指标的信息熵E_j，权重ω_j=(1-E_j)/∑(1-E_j)。基于100组历史数据计算得α =0.4、β =0.3、γ =0.2、δ =0.1。 3.2.4 仿真实验与结果分析 传递函数为G(s)=(Ke^(-θs))/(Ts+1)=(0.8e^(-8s))/(45s+1)， 1. 实验配置 参数范围：K_p∈[0,5]、K_i∈[0,0.1]、K_d∈[0,10]； 算法参数：粒子数N=50，最大迭代次数K_max=100； 对比基准：传统Z-N法PID与试凑法PID。 优化结果 经过100代迭代，PSO收敛至最优参数：K_p=1.85、K_i=0.023、K_d=7.42，对比实验表明： 超调量：PSO-PID为3.2%，较Z-N法（18.5%）和试凑法（6.7%）分别降低82.7%和52.2%； 调节时间：PSO-PID为42s，较对照组缩短40%以上； 稳态误差：PSO-PID为0.2℃，显著优于对照组的0.8℃（Z-N法）与0.4℃（试凑法）； 能耗：PSO-PID总能耗8.9kJ，较Z-N法降低29.4%。3. 动态响应对比 如图3-11所示，PSO-PID在阶跃响应中表现出平滑上升曲线，无显著超调；在t=150s水量突变干扰下，最大偏差仅1.5℃，20s内恢复稳态，鲁棒性显著优于传统方法。 {width=“6.0in” height=“3.0in”} 图3-11 PSO-PID与传统PID的动态响应对比，根据我提供的数据直接帮我绘制PSO-PID与传统PID的动态响应对比图的Matlab代码（代码里要用到非线性动态惯性权重和自适应因子）

<think>嗯，用户让我帮忙根据他们提供的3.2.2到3.2.4节的内容，生成绘制PSO-PID与传统PID动态响应对比图的Matlab代码。首先，我需要理解用户的需求。他们希望代码里包含非线性动态惯性权重和自适应学习因子，这说明他们可能已经实现了改进的PSO算法，现在需要可视化结果。 首先，我得确认用户的数据结构和参数。用户提到传递函数是G(s)=(0.8e^(-8s))/(45s+1)，而优化后的PID参数是Kp=1.85, Ki=0.023, Kd=7.42。传统方法如Z-N法和试凑法的参数虽然没有给出具体数值，但可能需要用户自行补充，或者用户可能有这些数据。不过，代码示例中可能需要假设或使用示例参数，所以可能需要注释说明让用户替换。 接下来，动态响应对比图需要模拟阶跃响应以及干扰情况下的响应。用户提到在t=150秒时水量突变，导致最大偏差1.5℃，并在20秒内恢复。因此，代码需要分两个部分：阶跃响应和干扰响应。可能需要使用Matlab的step函数和自定义的干扰模拟。 然后，关于PSO算法的改进部分，用户要求代码中使用非线性动态惯性权重和自适应学习因子。但实际绘图部分可能不需要运行PSO优化过程，因为用户已经得到了优化后的PID参数。所以代码可能只需要构建PID控制器，模拟响应，并绘制结果。不过，用户可能在代码中希望看到如何应用这些改进策略，或者他们可能需要完整的PSO实现来复现结果。但根据问题描述，用户直接要求绘制对比图，所以可能只需要使用优化后的参数进行仿真。 不过，用户可能误解了问题，以为生成对比图需要重新运行PSO算法，但实际上他们可能只需要用现有的参数进行仿真。这时候需要明确代码的结构：定义系统模型，设置PID参数，进行仿真，并绘图。因此，代码可能不需要包含PSO算法本身，除非用户需要展示PSO的训练过程。但根据用户提供的描述，他们可能已经完成了优化，现在需要对比结果。 但用户明确提到代码里要使用非线性动态惯性权重和自适应因子，这可能意味着需要展示PSO的实现，但绘图部分可能不需要。这有些矛盾，可能需要进一步澄清，但根据用户的问题，他们可能希望代码包含PSO的实现，以便复现结果，同时绘制对比图。 不过，考虑到用户提供的优化结果，可能他们已经运行了PSO，现在需要比较传统PID和PSO-PID的动态响应。因此，代码可能只需要使用优化后的参数进行仿真，而不需要重新实现PSO。但用户特别指出代码中要包含非线性惯性权重和自适应因子，所以可能需要展示PSO的代码部分，即使绘图部分不需要。 这可能意味着用户希望有一个完整的PSO实现，包括改进策略，然后应用优化后的参数进行仿真。但这样的代码会比较长，而且需要确保PSO正确找到参数。不过，根据用户提供的参数，可能他们已经运行了优化，所以代码只需要进行仿真比较。 但用户可能希望代码既包括PSO算法的实现（使用改进策略），又包括仿真对比。这会比较复杂，尤其是PSO需要优化PID参数，这需要定义适应度函数，也就是用户提到的性能指标。这可能需要较多的代码，包括传递函数的模拟，适应度计算，PSO的迭代过程等。 现在需要判断用户的具体需求：是仅仅生成对比图的仿真代码，还是包括PSO优化过程的完整代码。根据用户的问题描述，他们提到“根据我提供的数据直接帮我绘制”，所以可能已经有了优化后的参数，只需要仿真和绘图。但用户同时希望代码里用到非线性动态惯性权重和自适应因子，这可能意味着他们希望看到PSO的实现部分，即使绘图不需要。 因此，可能的解决方法是提供两部分代码：一部分是PSO算法的实现（包含改进策略），另一部分是仿真对比。但为了简化，可能将PSO的实现和仿真分开，或者提供注释说明。但考虑到用户可能希望完整的代码，需要整合这些部分。 另外，传递函数包含时滞环节（e^(-8s)），在Matlab中处理时滞可以使用tf函数的'InputDelay'属性，或者pade近似。但精确的时滞仿真可能需要特定的方法。例如，使用step函数时，可能需要用pade近似来处理延迟，或者使用Simulink模型。但为了代码简洁，可能在代码中使用pade近似来近似时滞。 总结步骤： 1. 定义被控对象的传递函数，包括时滞处理。 2. 定义传统PID的参数（Z-N法和试凑法），需要用户提供或假设数值。 3. 定义PSO优化后的PID参数。 4. 构建PID控制系统，进行阶跃响应仿真。 5. 模拟干扰（t=150s时可能是一个负载扰动，如何在代码中实现？可能需要使用lsim或自定义函数）。 6. 绘制对比图，包括超调量、调节时间等指标。 但用户提供的描述中，动态响应对比图包括阶跃响应和干扰响应，因此可能需要两个仿真：一个是正常的阶跃响应，另一个是在t=150s时加入干扰。干扰可能体现为阶跃信号叠加，或者改变参考输入，或者加入扰动输入。需要明确如何模拟水量突变的影响。例如，可能是在t=150s时，参考输入变化，或者加入一个阶跃扰动。 假设在t=150s时，系统受到一个阶跃干扰，可以通过生成一个时间向量，并在150秒后叠加一个干扰信号。例如，使用lsim函数模拟系统对干扰的响应。 此外，用户提供的传递函数是开环的，而PID控制通常是闭环系统。需要构建闭环传递函数，或者使用feedback函数。例如，假设是单位反馈系统，控制器为PID，被控对象为G(s)。 综合这些因素，代码的结构可能如下： - 定义被控对象G(s)（包括时滞）。 - 定义PID控制器的传递函数（根据参数Kp, Ki, Kd）。 - 构建闭环系统。 - 仿真阶跃响应。 - 仿真干扰响应（例如，在t=150s时加入阶跃扰动）。 - 绘制结果。 同时，如果用户希望代码中包含PSO的实现，那么需要： - 实现改进的PSO算法，包括非线性惯性权重和自适应学习因子。 - 定义适应度函数，基于四个性能指标。 - 运行PSO优化得到PID参数。 - 然后进行仿真对比。 但这样代码量会很大，而且需要确保适应度函数的正确性，包括超调量、调节时间等的计算。这可能比较复杂，尤其是调节时间和超调量需要从阶跃响应数据中提取。 考虑到用户可能已经拥有优化后的参数，并且只需要绘图，可能更合理的代码是仅进行仿真对比，使用提供的优化参数。但用户明确要求代码中包含非线性动态惯性权重和自适应因子，所以可能需要包括PSO的实现部分，即使最后绘图不需要运行优化过程。 这似乎有些矛盾，可能需要进一步确认，但根据问题描述，用户可能希望代码中包含这些改进的PSO策略，以展示如何得到优化参数，然后进行仿真对比。 因此，代码可能需要分为几个部分： 1. 定义系统模型和仿真参数。 2. 实现改进的PSO算法，优化PID参数。 3. 使用优化后的参数和传统方法参数进行仿真。 4. 绘制对比图。 但考虑到用户已经提供了优化后的参数（Kp=1.85, Ki=0.023, Kd=7.42），可能不需要重新运行PSO优化，除非为了验证。但用户可能需要完整的可运行代码，包括PSO部分，以生成结果。否则，直接使用提供的参数进行仿真即可。 因此，可能的折中方案是提供PSO算法的框架代码，包括改进策略，但注释掉实际的优化过程，直接使用用户提供的优化参数进行仿真和绘图。这样既满足代码中包含改进策略的要求，又避免长时间运行PSO优化。 现在，具体到代码实现： 1. 处理时滞：使用pade近似，例如三阶pade近似，将时滞转换为有理函数近似。 2. 定义被控对象： ```matlab s = tf('s'); K = 0.8; theta = 8; % 时滞 T = 45; % 使用pade近似处理时滞 [num_pade, den_pade] = pade(theta, 3); % 三阶pade近似 G_pade = tf(num_pade, den_pade); G = (K / (T*s + 1)) * G_pade; ``` 3. 定义PID控制器函数： ```matlab function C = pid_controller(Kp, Ki, Kd) C = Kp + Ki/s + Kd*s; end ``` 4. 闭环系统的构建： ```matlab function sys_cl = closed_loop(Kp, Ki, Kd, G) C = pid_controller(Kp, Ki, Kd); sys_open = C * G; sys_cl = feedback(sys_open, 1); end ``` 5. 仿真阶跃响应和干扰响应： 对于阶跃响应，使用step函数。对于干扰响应，可能需要定义一个输入信号，例如在t=150秒时加入阶跃干扰，可以使用lsim函数。 6. 绘制对比图： ```matlab figure; subplot(2,1,1); plot(t_pso, y_pso, 'r', t_zn, y_zn, 'b--', t_trial, y_trial, 'g-.'); title('Step Response'); legend('PSO-PID', 'Z-N PID', 'Trial PID'); subplot(2,1,2); % 绘制干扰响应 ``` 关于干扰模拟，假设在t=150秒时，系统参考输入不变，但加入一个阶跃扰动，例如负载变化。这可以通过在控制回路中添加扰动输入点，然后使用lsim模拟扰动的影响。例如，扰动可能加在被控对象的输入或输出端。 假设扰动d(t)是一个在t=150秒时的阶跃信号，幅值为某个值（例如，导致温度变化1.5℃）。需要确定扰动的位置和大小，但用户描述中提到“水量突变干扰下，最大偏差仅1.5℃”，可能需要调整扰动的大小以达到该效果。 这部分可能需要更详细的系统模型，但为了简化，可能直接在控制信号中加入干扰，或者改变参考输入。不过，用户的问题可能更关注于展示PSO-PID在干扰下的恢复能力，因此代码需要模拟这种情况。 综上，代码的大致流程是： - 定义被控对象，使用pade近似处理时滞。 - 定义传统PID和PSO-PID的参数。 - 构建闭环系统，仿真阶跃响应。 - 仿真干扰响应（在t=150秒加入干扰）。 - 绘制对比图。 但如何具体实现干扰响应呢？可能需要构建一个带有扰动的系统模型，例如，扰动加在过程输出或输入。 例如，闭环系统可以表示为： 输出Y = G*(C*(R - Y) + D) 其中D是扰动输入。因此，扰动的影响可以通过仿真D的响应。 在Matlab中，可以使用以下方法： ```matlab % 闭环系统不考虑扰动时的传递函数 sys_cl = feedback(C*G, 1); % 扰动到输出的传递函数 sys_disturbance = feedback(G, C); ``` 然后，仿真系统对阶跃参考信号和扰动信号的响应。 例如，在t=0时施加阶跃参考信号，在t=150秒时施加阶跃扰动信号。这可以通过拼接时间信号来实现： ```matlab t_total = 0:0.1:300; % 总仿真时间300秒 r = ones(size(t_total)); % 阶跃参考输入 d = zeros(size(t_total)); d(t_total >= 150) = 0.5; % 假设扰动幅值为0.5，需要调整以达到1.5℃偏差 % 仿真参考输入响应 y_ref = lsim(sys_cl, r, t_total); % 仿真扰动响应 y_dist = lsim(sys_disturbance, d, t_total); % 总输出 y_total = y_ref + y_dist; ``` 但需要确保系统模型正确，包括扰动的位置。这可能涉及到控制系统的结构，例如扰动是加在被控对象的输入还是输出端。 此外，用户提到的“水量突变干扰”可能影响系统动态，可能需要调整模型参数，但可能超出范围，因此简化为添加一个阶跃扰动。 综合所有因素，最终的Matlab代码可能如下： ```matlab %% 参数初始化 clear; clc; % 被控对象传递函数（含时滞处理） s = tf('s'); K = 0.8; theta = 8; % 时滞 T = 45; [nump, denp] = pade(theta, 3); % 三阶Pade近似 Gp = tf(nump, denp); % 时滞的近似模型 G = (K / (T*s + 1)) * Gp; % 被控对象传递函数 % PID参数 % PSO优化后的参数 Kp_pso = 1.85; Ki_pso = 0.023; Kd_pso = 7.42; % 传统方法参数（示例，需根据实际情况调整） Kp_zn = 1.2; % 假设Z-N法参数 Ki_zn = 0.05; Kd_zn = 4.0; Kp_trial = 1.5; % 假设试凑法参数 Ki_trial = 0.03; Kd_trial = 5.0; %% 动态响应仿真 % 时间向量 t = 0:0.1:300; % 仿真时间300秒 % 生成阶跃输入信号（参考信号） r = ones(size(t)); % 扰动信号（在t=150s加入阶跃扰动） d = zeros(size(t)); d(t >= 150) = 0.2; % 扰动幅值根据实际情况调整 % 初始化输出存储 y_pso = zeros(size(t)); y_zn = zeros(size(t)); y_trial = zeros(size(t)); % 定义闭环系统函数 function y = simulate_system(Kp, Ki, Kd, G, r, d, t) C = Kp + Ki/s + Kd*s; % PID控制器 % 闭环系统（参考输入到输出） sys_cl = feedback(C*G, 1); % 扰动到输出的传递函数 sys_dist = feedback(G, C); % 分别仿真参考输入和扰动的影响 y_ref = lsim(sys_cl, r, t); y_dist = lsim(sys_dist, d, t); y = y_ref + y_dist; end % 仿真各PID控制器的响应 y_pso = simulate_system(Kp_pso, Ki_pso, Kd_pso, G, r, d, t); y_zn = simulate_system(Kp_zn, Ki_zn, Kd_zn, G, r, d, t); y_trial = simulate_system(Kp_trial, Ki_trial, Kd_trial, G, r, d, t); %% 绘图 figure('Position', [100, 100, 800, 600]); plot(t, y_pso, 'r-', 'LineWidth', 1.5); hold on; plot(t, y_zn, 'b--', 'LineWidth', 1.5); plot(t, y_trial, 'g-.', 'LineWidth', 1.5); hold off; grid on; % 标注关键事件 line([150, 150], ylim, 'Color', 'k', 'LineStyle', '--'); text(150, 2.5, '水量突变干扰', 'Rotation', 90, 'VerticalAlignment', 'bottom'); title('PSO-PID与传统PID动态响应对比'); xlabel('时间 (s)'); ylabel('温度 (℃)'); legend('PSO-PID', 'Z-N PID', '试凑法 PID', 'Location', 'best'); %% 性能指标计算（可选） % 可添加代码计算超调量、调节时间等指标 ``` 注意：传统PID的参数（Z-N法和试凑法）需要用户根据实际情况提供，代码中的参数是假设的，可能需要调整。此外，扰动幅值（0.2）可能需要调整以达到用户提到的1.5℃偏差。 如果用户希望代码中包含PSO算法的实现，则需要添加以下部分： ```matlab %% 改进PSO算法实现 function [gBest, gBestJ] = improved_pso(G, K_max, N) % 参数范围 Kp_range = [0, 5]; Ki_range = [0, 0.1]; Kd_range = [0, 10]; % 惯性权重参数 w_max = 0.9; w_min = 0.4; % 学习因子 c1_max = 2.5; c1_min = 1.5; c2_min = 1.5; c2_max = 2.5; % 适应度函数权重 alpha = 0.4; beta = 0.3; gamma = 0.2; delta = 0.1; % 初始化粒子群 particles = struct('position', [], 'velocity', [], 'pBest', [], 'pBestJ', []); for i = 1:N Kp = Kp_range(1) + rand()*(Kp_range(2) - Kp_range(1)); Ki = Ki_range(1) + rand()*(Ki_range(2) - Ki_range(1)); Kd = Kd_range(1) + rand()*(Kd_range(2) - Kd_range(1)); particles(i).position = [Kp, Ki, Kd]; particles(i).velocity = zeros(1,3); particles(i).pBest = particles(i).position; particles(i).pBestJ = inf; end gBest = [1.85, 0.023, 7.42]; % 初始猜测，实际应由算法更新 gBestJ = inf; for k = 1:K_max % 更新惯性权重和學習因子 w = w_min + (w_max - w_min)*exp(-10*k/K_max); c1 = c1_max - (c1_max - c1_min)*k/K_max; c2 = c2_min + (c2_max - c2_min)*k/K_max; for i = 1:N % 计算适应度 Kp = particles(i).position(1); Ki = particles(i).position(2); Kd = particles(i).position(3); J = calculate_fitness(Kp, Ki, Kd, G, alpha, beta, gamma, delta); % 更新个体最优 if J < particles(i).pBestJ particles(i).pBest = particles(i).position; particles(i).pBestJ = J; end % 更新全局最优 if J < gBestJ gBest = particles(i).position; gBestJ = J; end % 更新速度 r1 = rand(); r2 = rand(); particles(i).velocity = w*particles(i).velocity + ... c1*r1*(particles(i).pBest - particles(i).position) + ... c2*r2*(gBest - particles(i).position); % 更新位置 new_pos = particles(i).position + particles(i).velocity; % 边界处理 new_pos(1) = max(min(new_pos(1), Kp_range(2)), Kp_range(1)); new_pos(2) = max(min(new_pos(2), Ki_range(2)), Ki_range(1)); new_pos(3) = max(min(new_pos(3), Kd_range(2)), Kd_range(1)); particles(i).position = new_pos; end end end function J = calculate_fitness(Kp, Ki, Kd, G, alpha, beta, gamma, delta) % 仿真获取系统响应 t = 0:0.1:300; r = ones(size(t)); d = zeros(size(t)); d(t >= 150) = 0.2; % 仿真响应 y = simulate_system(Kp, Ki, Kd, G, r, d, t); % 计算性能指标 [sigma, t_s, e_ss, energy] = calculate_performance(y, t, r); % 适应度值 J = alpha*sigma + beta*t_s + gamma*e_ss + delta*energy; end function [sigma, t_s, e_ss, energy] = calculate_performance(y, t, r) % 计算超调量、调节时间、稳态误差和能耗 % 此处需要具体实现，例如： % 超调量计算 [ymax, idx] = max(y); sigma = (ymax - r(end)) / r(end) * 100; % 调节时间：进入±1%误差带 idx_settle = find(abs(y - r(end)) <= 1, 1, 'first'); t_s = t(idx_settle); % 稳态误差（取最后10%时间段的平均） idx_ss = floor(0.9*length(t)):length(t); e_ss = mean(abs(y(idx_ss) - r(idx_ss))); % 能耗（需要控制量u(t)的数据，此处简化） energy = 0; % 需要根据控制信号计算积分 end ``` 不过，这部分PSO的实现较为复杂，且需要正确的性能指标计算，尤其是能耗部分需要访问控制信号u(t)，这需要修改仿真函数以输出控制信号。此外，适应度函数的计算可能需要更精确的方法，例如使用step函数获取阶跃响应的数据，并从中提取超调量、调节时间等。 由于用户的主要需求是绘制对比图，而不是重新运行PSO优化，因此可能更合适仅提供仿真和绘图的代码，假设优化后的参数已经得到，并包含在代码中。而PSO的实现部分可以作为附加内容，但需要较多调试。 最终，用户可能需要一个能够直接运行的代码，使用提供的参数进行仿真并生成对比图。因此，将重点放在仿真和绘图部分，并在注释中说明PSO的改进策略，可能更符合用户的实际需求。</think>以下是根据您提供的系统描述与参数设计的Matlab动态响应对比图代码，包含改进粒子群算法核心要素： ```matlab %% 系统建模与参数初始化 clc; clear; s = tf('s'); % 被控对象传递函数（含时滞环节） K = 0.8; theta = 8; T = 45; [G_pade,~] = pade(theta,3); % 三阶Pade近似 G = (K/(T*s+1))*G_pade; % 总传递函数 % 优化后的PID参数 PSO_params = [1.85, 0.023, 7.42]; % [Kp, Ki, Kd] ZN_params = [1.2, 0.05, 4.0]; % 示例Z-N法参数（需替换实际值） Trial_params = [1.5, 0.03, 5.0]; % 示例试凑法参数（需替换实际值） %% 动态响应仿真 % 仿真时间设置 t = 0:0.1:300; r = ones(size(t)); % 阶跃参考信号 d = zeros(size(t)); % 干扰信号 d(t>=150) = 0.3; % 150s时加入干扰 % 各控制器响应仿真 [y_pso, t] = simPID(PSO_params, G, r, d, t); [y_zn, ~] = simPID(ZN_params, G, r, d, t); [y_trial, ~] = simPID(Trial_params, G, r, d, t); %% 绘制动态响应曲线 figure('Position', [100 100 800 600]) plot(t, y_pso, 'r-', 'LineWidth', 2), hold on plot(t, y_zn, 'b--', 'LineWidth', 1.5) plot(t, y_trial, 'g-.', 'LineWidth', 1.5) plot([150 150], ylim, 'k--', 'LineWidth', 1) % 干扰标记 text(155, 2.5, '水量突变干扰', 'Rotation',90,'FontSize',10) title('PSO-PID与传统PID动态响应对比') xlabel('时间(s)'), ylabel('温度(℃)') legend('PSO-PID','Z-N PID','试凑法PID','Location','best') grid on %% 改进PSO算法核心函数（用于参数优化） function [gBest] = improvedPSO(G) % 算法参数 K_max = 100; N = 50; w_max = 0.9; w_min = 0.4; c1_max = 2.5; c1_min = 1.5; c2_min = 1.5; c2_max = 2.5; % 粒子群初始化 particles = struct('position',[],'velocity',[],'pBest',[],'pBestJ',inf); for i=1:N particles(i).position = [rand*5, rand*0.1, rand*10]; % [Kp,Ki,Kd] particles(i).velocity = zeros(1,3); end % 主循环 for k=1:K_max w = w_min + (w_max-w_min)*exp(-10*k/K_max); % 非线性惯性权重 c1 = c1_max - (c1_max-c1_min)*k/K_max; % 自适应学习因子 c2 = c2_min + (c2_max-c2_min)*k/K_max; % 粒子更新与适应度计算 for i=1:N % 计算适应度值 J = calculateFitness(particles(i).position, G); % 更新个体/全局最优 if J < particles(i).pBestJ particles(i).pBest = particles(i).position; particles(i).pBestJ = J; end % 更新速度与位置 particles(i).velocity = w*particles(i).velocity + ... c1*rand*(particles(i).pBest - particles(i).position) + ... c2*rand*(gBest - particles(i).position); % 边界约束 particles(i).position = max(min(particles(i).position + ... particles(i).velocity,[5,0.1,10]),[0,0,0]); end end end %% 辅助函数 function [y, t] = simPID(params, G, r, d, t) Kp = params(1); Ki = params(2); Kd = params(3); C = pid(Kp, Ki, Kd); % 创建PID控制器 CL = feedback(C*G, 1); % 闭环系统 y_ref = lsim(CL, r, t); % 参考响应 y_dist = lsim(feedback(G,C), d, t); % 干扰响应 y = y_ref + y_dist; % 总输出 end function J = calculateFitness(params, G) % 性能指标计算（需实现超调量、调节时间等具体计算） [y, ~] = simPID(params, G, ones(1,1000), zeros(1,1000), 0:0.1:100); % 此处应添加具体指标计算代码 J = sum(abs(y-1)); % 示例简化计算 end ``` 代码说明： 1. 系统建模部分使用Pade近似处理时滞环节，保证仿真精度 2. `simPID`函数封装了闭环系统仿真过程，包含参考跟踪与抗干扰能力测试 3. 改进PSO算法实现包含： - 非线性动态惯性权重策略 - 自适应学习因子调整 - 粒子位置的速度更新与边界约束 4. 动态响应曲线展示阶跃响应与抗干扰性能对比 注意事项： 1. 需根据实际Z-N法和试凑法参数替换示例值 2. `calculateFitness`函数需完善超调量、调节时间等具体指标计算 3. 干扰幅值0.3需根据实际系统特性调整 4. 绘图结果可直接导出为符合论文要求的矢量图 建议将PSO优化过程单独运行后存储最优参数，再调用`simPID`进行性能对比。实际应用时可调整仿真步长（当前0.1s）平衡精度与计算速度。