基于SOA海鸥优化算法的二阶时滞系统PID控制器最优参数计算matlab仿真

1.算法仿真效果
matlab2022a仿真结果如下：

2.算法涉及理论知识概要

   一般来说，海鸥是群居性的，它们用自己的智慧来寻找并攻击猎物。海鸥最重要的是它们的迁徙和攻击行为。迁徙被定义为海鸥从一个地方到另一个地方的季节性迁移，以寻找最丰富食物来源，以提供足够的能量。该行为可描述为：

    在迁徙过程中，它们成群结队地迁徙。海鸥的初始位置不同，以避免相互碰撞；
    在一个群体中，海鸥可以朝着最适合生存的海鸥的方向前进；
    其他海鸥可以根据最适合的海鸥更新它们的初始位置。
    海鸥在迁徙过程中，会经常攻击其他候鸟，在攻击过程中，它们可以做出螺旋形的自然运动，如图1所示。

   海鸥的迁徙和攻击行为。海鸥根据季节更替进行迁徙，迁徙飞行时海鸥会避免相互碰撞；海鸥会攻击猎物，攻击时呈螺旋形的运动形态；在一个群体中，海鸥朝着最佳位置的方向前进。

1.初始化参数；
2.初始化种群位置；
3.计算适应度值并保留全局最优位置；
4.迁徙，全局搜索：
抽象海鸥的迁徙行为主要有三步，第一要满足避免碰撞条件（这样可以确保种群多样性）；第二要计算最佳位置的方向；第三要根据该方向移动到新的位置。

SOA的迭代伪代码如下图6所示：

3.MATLAB核心程序

Umax=0.9500;%最大隶属度值
Umin=0.0111;%最小隶属度值
Wmax=0.9;%权重最大值
Wmin=0.1;%权重最小值
Dim = 3;            % 维数
SwarmSize =30;    % 粒子群规模
MaxIter = 100;      % 最大迭代次数
MinFit = 10;       % 最小适应值
Ub = [100 100 100];
Lb = [0 0 0];
%%种群初始化
Range = ones(SwarmSize,1)*(Ub-Lb);
Swarm = rand(SwarmSize,Dim).*Range + ones(SwarmSize,1)*Lb;    % 初始化粒子群
fSwarm = zeros(SwarmSize,1);
for i=1:SwarmSize
    fSwarm(i,:) = PID_SOA(Swarm(i,:));                       % 粒子群的适应值
end
%%个体极值和群体极值
[bestf bestindex]=min(fSwarm);
zbest=Swarm(bestindex,:);   % 全局最佳
gbest=Swarm;                % 个体最佳
fgbest=fSwarm;              % 个体最佳适应值
fzbest=bestf;               % 全局最佳适应值
%迭代寻优
Di=0*rand(SwarmSize,Dim);
Buchang=0*rand(SwarmSize,Dim);
C=0*rand(SwarmSize,Dim);
Diego=0*rand(SwarmSize,Dim);
Dialt=0*rand(SwarmSize,Dim);
Dipro=0*rand(SwarmSize,Dim);
iter = 0;
y_fitness = zeros(1,MaxIter);   % 预先产生4个空矩阵
K_p = zeros(1,MaxIter);         
K_i = zeros(1,MaxIter);
K_d = zeros(1,MaxIter);
while( (iter < MaxIter) && (fzbest > MinFit) )
    for i=1:SwarmSize
        W=Wmax-iter*(Wmax-Wmin)/MaxIter;
        Diego(i,:)=sign(gbest(i,:) -Swarm(i,:));%确定利己方向
        Dialt(i,:)=sign(zbest -Swarm(i,:));%确定利他方向
        if PID_SOA(gbest(i,:))>=PID_SOA(Swarm(i,:))%确定预动方向
            Dipro(i,:)=-Di(i,:);
        else
            Dipro(i,:)=Di(i,:);
        end
        Di(i,:)=sign(W* Dipro(i,:)+0.5*Diego(i,:)+0.5*Dialt(i,:));%确定经验梯度方向
        [Orderfgbest,Indexfgbest]=sort(fgbest,'descend');
        u=Umax-(SwarmSize-Indexfgbest(i))*(Umax-Umin)/(SwarmSize-1);
        U=u+(1-u)*rand;
        H=(MaxIter-iter)/MaxIter;%迭代过程中权重的变化
        C(i,:)=H*abs(zbest-10*rand(1,3));%确定高斯函数的参数
        T=sqrt(-log(U));
        Buchang(i,:)=C(i,:)*T;%确定搜索不长的大
        Buchang(i,find(Buchang(i,:)>3*max(C(i,:))))=3*max(C(i,:));
        %更新位置
        Swarm(i,:)=Swarm(i,:)+Di(i,:).*Buchang(i,:);
        Swarm(i,find(Swarm(i,:)>100))=100;
        Swarm(i,find(Swarm(i,:)<0))=0;
        % 适应值
        fSwarm(i,:) =PID_SOA(Swarm(i,:));
        % 个体最优更新     
        if fSwarm(i) < fgbest(i)
            gbest(i,:) = Swarm(i,:);
            fgbest(i) = fSwarm(i);
        end
        % 群体最优更新
        if fSwarm(i) < fzbest
            zbest = Swarm(i,:);
            fzbest = fSwarm(i);
        end
    end 
    iter = iter+1;                      % 迭代次数更新
    y_fitness(1,iter) = fzbest;         % 为绘图做准备
    K_p(1,iter) = zbest(1);
    K_i(1,iter) = zbest(2);
    K_d(1,iter) = zbest(3);
end
%%绘图
figure(1)      % 绘制性能指标ITAE的变化曲线
plot(y_fitness,'LineWidth',2)
title('最优个体适应值');
xlabel('迭代次数');ylabel('适应值');
grid on
 
figure(2)      % 绘制PID控制器参数变化曲线
plot(K_p,'LineWidth',2)
hold on
plot(K_i,'k','LineWidth',2)
hold on
plot(K_d,'--r','LineWidth',2)
title('Kp、Ki、Kd 优化曲线');
xlabel('迭代次数');
ylabel('参数值');
legend('Kp','Ki','Kd');
grid on