m基于梯度优化的混沌PSO磁悬浮球系统模型优化的matlab仿真

1.算法描述

     基本确定融合后的优化控制算法，即混沌粒子群优化算法，将该算法应用于对PID的参数整定上，通过仿真验证算法的可行性。将混沌思想引入PSO算法，前期工作首先对混沌算法局部搜索能力差和算法可能需要花费长时间才能取得较好优化性能的不足进行了改进，提出梯度优化混沌算法的思想，基本实现了梯度算法与混沌算法的融合，实现了用梯度算法改善混沌算法的预期目的，达到了取长补短的效果，使得该算法即高效又不容易陷入局部极小，为后期将优化后的混沌算法引入PSO算法做好了准备。主要研究如何把PSO算法与得到的混沌算法有效的结合，实现预期目的，用混沌算法来改善PSO算法，使其能持续在全局范围搜索，避免算法早熟。研究实现如下思路的混沌粒子群优化算法：将混沌映射直接映入PSO算法的迭代过程，使算法能持续在全局范围搜索，避免算法早熟，但对于算法后期收敛精度不高和收敛速度下降的问题，此时考虑再将梯度算法引入，以粒子群当前获得的最优粒子所在位置为初始点进行梯度搜索，提高局部搜索效率。再将混沌PSO和梯度算法的搜索过程进行循环，得到新的混沌PSO优化算法。将该算法应用到PID参数整定上，与前期所设计的算法进行比较，得到更优的。

   我们需要构建磁悬浮球系统模型的控制对象模型。这里，我们通过传递函数的形式来实现磁悬浮球的表达式，具体的推导过程如下所示：

   上图中，Xo为磁悬浮球的平衡位置的间隙，x为钢球在Y方向上的偏离平衡位置的位移。Io和i分别为电磁铁线圈的偏置电流和控制电流，F为电磁铁对钢球所产生的电磁铁，mg为钢球所受重力。为了保持平衡，合理F-mg必须为0.假设平衡状态下，电流为i0，钢球和电磁铁之间的气隙为x0，那么这个时候有如下的表达式：

   其中K0和P，通过预先给定的参数，可以获得，其为两个固定的常数，这个可以根据实际的参数自己设置，我们将上述的Gs作为控制器系统的控制对象进行控制处理。

2.仿真效果预览
matlab2022a仿真结果如下：

3.MATLAB核心程序

global yd y timef
%定义磁悬浮球传递函数
tmac = 1000;
den  =[1,32,2,1];
 
wmax = 0.9;
wmin = 0.1;
c1   = 0.3;
c2   = 0.3;
%Kp的范围
amin = 0;
amax = 100;
%Ki的范围
bmin = 0;
bmax = 100;
%Kd的范围
cmin = 0;
cmax = 100;
%速度的范围
vmin =-0.1;
vmax = 0.1;
%粒子数目
Pop  = 5;
BsJ  = 0;
%迭代次数
T    = 200;
%粒子初始化
for i=1:Pop
    Kpid(i,1)   = rand(1)*(amax-amin)+amin;
    a_best(i)   = Kpid(i,1);
    Kpid(i,2)   = rand(1)*(bmax-bmin)+bmin;
    b_best(i)   = Kpid(i,2);
    Kpid(i,3)   = rand(1)*(cmax-cmin)+cmin;
    c_best(i)   = Kpid(i,3);
    va(i)       =(vmax-vmin)*rand(1)+vmin;
    vb(i)       =(vmax-vmin)*rand(1)+vmin;
    vc(i)       =(vmax-vmin)*rand(1)+vmin;
    Kpidi       = Kpid(i,:);  
    [BsJ,Y]     = func_pid_controller_fitness(Kpidi,tmac,den);   
    BsJi(i)     = BsJ;
end
a            = Kpid(:,1);
b            = Kpid(:,2);
c            = Kpid(:,3);
[minJi,index]= min(BsJi);
Ta_best      = a(index);
Tb_best      = b(index);
Tc_best      = c(index); 
 
for t=1:T
    
    time(t) = t;
    w       = wmax-t*(wmax-wmin)/T;
    for i=1:Pop
        va(i)     = w*va(i)+c1*rand(1)*(a_best(i)-Kpid(i,1))+c2*rand(1)*(Ta_best-Kpid(i,1));
        Kpid(i,1) = Kpid(i,1)+va(i);
        if Kpid(i,1)>=amax
           Kpid(i,1)=amax;
        end
        if Kpid(i,1)<=amin
           Kpid(i,1)=amin;
        end
        
        vb(i)     = w*vb(i)+c1*rand(1)*(b_best(i)-Kpid(i,2))+c2*rand(1)*(Tb_best-Kpid(i,2));
        Kpid(i,2) = Kpid(i,2)+vb(i);
        if Kpid(i,2)>=bmax
           Kpid(i,2)=bmax;
        end
        if Kpid(i,2)<=bmin
           Kpid(i,2)=bmin;
        end
        
        vc(i)     = w*vc(i)+c1*rand(1)*(c_best(i)-Kpid(i,3))+c2*rand(1)*(Tc_best-Kpid(i,3));
        Kpid(i,3) = Kpid(i,3)+vc(i);
        if Kpid(i,3)>=cmax
           Kpid(i,3)=cmax;
        end
        if Kpid(i,3)<=cmin
           Kpid(i,3)=cmin;
        end
    
        Kpidi       = Kpid(i,:);
        [BsJ,Y]     = func_pid_controller_fitness(Kpidi,tmac,den);
 
        if BsJ<BsJi(i)
           BsJi(i)   = BsJ;
           a_best(i) = Kpid(i,1);
           b_best(i) = Kpid(i,2);
           c_best(i) = Kpid(i,3);
        end
        if BsJi(i)<minJi
           minJi   = BsJi(i);
           Ta_best = Kpid(i,1);
           Tb_best = Kpid(i,2);
           Tc_best = Kpid(i,3);
        end
    end
    Jibest(t) = minJi;
    A_best(t) = Ta_best;
    B_best(t) = Tb_best;
    C_best(t) = Tc_best;
end
 
figure;
plot(Jibest,'r','linewidth',2);
xlabel('迭代次数');
ylabel('J');
grid on
 
 
save R1.mat Jibest timef yd y
 
figure;
subplot(311);
plot(time,A_best,'b','linewidth',2);
xlabel('times(s)');
ylabel('Kp');
grid on
subplot(312);
plot(time,B_best,'b','linewidth',2);
xlabel('times(s)');
ylabel('Ki');
grid on
subplot(313);
plot(time,C_best,'b','linewidth',2);
xlabel('times(s)');
ylabel('Kd');
grid on