basicpso.m

粒子群算法

%函数说明：该函数只是最基本的PSO算法的版本，且由于本人学习matlab时间不长，所以质量不高，望谅解。
%进一步的工作可以有：发挥matlab图形显示的长处，三维显示整个优化过程，加入各种改进算法，如协同PSO,扰动等等。
%发现本程序运行结果不及VC，不知道是不是随机数的问题，不过基本程序框架是这样的
%任何问题请联系：缘来是你~（38171545），谢谢！

function exepso()
     [fitness,dimens,gbest]=MyPSO;%调用PSO函数
      %test function1
     %x=[-5.12:0.05:5.12];
     %y=zeros(1,205);
     %for i=1:204
	%	y(i)=y(i)+(x(i)*x(i)-10*cos(2*pi*x(i))+10);
    % end
    %plot(x,y,'--r*',x,fitness,'-b');
    %test function1
function [XX,fitness]=Calculate(XX,ddimens)%计算函数的适应度值，也就是把自己的优化函数放进去
%........................
    %fitness=0;
     %test function1
	%for i=1:ddimens
		%if(XX(i)>5.12||XX(i)<-5.12)
			%XX(i)=10.24*rand-5.12;
            %    end
            %    end
%	for i=1:ddimens
	%	fitness=fitness+(XX(i)*XX(i)-10*cos(2*pi*XX(i))+10);
    %  end
     %test function1
     for i=1:ddimens%检查是不越界，是则纠正
         if XX(i)>200||XX(i)<0
             XX(i)=200*rand;
         end
     end

	if ((0.0193*XX(1)-0.0625*XX(3))<=0)&&((0.00954*XX(1)-0.0625*XX(4))<=0)&&((750*1728-pi*XX(1)*XX(1)*XX(2)-4*pi*XX(1)*XX(1)*XX(1)/3)<=0)&&((XX(2)-240)<=0)
       	fitness=0.6224*0.0625*XX(3)*XX(1)*XX(2)+1.7781*0.0625*XX(4)*XX(1)*XX(1)+3.1661*0.0625*XX(3)*0.0625*XX(3)*XX(2)+19.84*0.0625*XX(3)*0.0625*XX(3)*XX(1);
    else
        fitness=10000000;
    end

function [s,dimens,Pp]=MyPSO()%PSO主函数
pso_size=30;%种群大小
%test function1
%dimens=30;
dimens=4;%待优化问题的维数
Pb=10000000;%适应度值
w_max=0.9;%最大权重
w_min=0.4;%最小权重
w=w_max;
c1=2;c2=2;%学习因子取值
run_max=5000;%迭代次数上限
for i=1:pso_size%粒子位置、速度初始化
    Yb(i)=0;
    for j=1:dimens
        Xb(i,j)=200*rand;%Xb记录个体最优位置
        X(i,j)=Xb(i,j);
        V(i,j)=150*rand-75;
    end
end
for i=1:pso_size
    [X(i,:),Yb(i)]=Calculate(X(i,:),dimens);%计算适应度值
    if Yb(i)<Pb
        Pb=Yb(i);%Pb为当前群体最优适应度值
        for j=1:dimens
            Pp(j)=X(i,j);%Pp记录群体最优位置
        end
    end
end
    
for count=1:run_max
    for k=1:pso_size
        for j=1:dimens
            V(k,j)=w*V(k,j)+c1*rand*(Xb(k,j)-X(k,j))+c2*rand*(Pp(j)-X(k,j));%速度更新公式
            if abs(V(k,j))>150%限制速度
                if V(k,j)>0
                    V(k,j)=150;
                else
                    V(k,j)=-150;
                end
            end
        end
            for i=1:dimens
                X(k,i)=X(k,i)+V(k,i);%位置更新公式
            end
               [X(i,:),temp]=Calculate(X(k,:),dimens);%计算适应度值
                if temp<=Yb(k)
                    Yb(k)=temp;
                    for i=1:dimens
                        Xb(k,i)=X(k,i);
                    end
                end
                    if Yb(k)<=Pb
                        Pb=Yb(k);
                        for i=1:dimens
                            Pp(i)=X(k,i);
                        end
                    end 
                        
            end
            w=w_max-count*(w_max-w_min)/run_max;%基本惯性权重的计算公式       
    end
    s=Pb