Matlab粒子群算法(PSO)进行空间-变压器局部放电定位

%%ParticleSwarmOptimization(PSO),粒子群优化算法计算局放源位置%%%%ClearScreenclc;closeall;clearall;%%InitializationLx=4;Ly=4;Lz=4;%变压器长宽高Sensor_xyz0=[1,0,1]';Sensor_xyz1=[0,1.5,0.5]';Sensor_xyz2=[2,3,0.5]';Sensor_xyz3=[3,1.5,0.7]';%四个传感器位置Source=[2,2,2]';%PD位置t=linspace(0,50e-9,1001);tc=10e-9;c=3e8;delta_t0=sqrt(sum((Source-Sensor_xyz0).^2))/c;delta_t1=sqrt(sum((Source-Sensor_xyz1).^2))/c;delta_t2=sqrt(sum((Source-Sensor_xyz2).^2))/c;delta_t3=sqrt(sum((Source-Sensor_xyz3).^2))/c;y0=-(exp(-(t-tc-delta_t0)/0.2e-9)-exp(-(t-tc-delta_t0)/3e-9)).*stepfun(t,tc+delta_t0).*sin(2*pi*1e9*(t-tc-delta_t0));%传感器0接收到PD信号y1=-(exp(-(t-tc-delta_t1)/0.2e-9)-exp(-(t-tc-delta_t1)/3e-9)).*stepfun(t,tc+delta_t1).*sin(2*pi*1e9*(t-tc-delta_t1));%传感器1接收到PD信号y2=-(exp(-(t-tc-delta_t2)/0.2e-9)-exp(-(t-tc-delta_t2)/3e-9)).*stepfun(t,tc+delta_t2).*sin(2*pi*1e9*(t-tc-delta_t2));%传感器2接收到PD信号y3=-(exp(-(t-tc-delta_t3)/0.2e-9)-exp(-(t-tc-delta_t3)/3e-9)).*stepfun(t,tc+delta_t3).*sin(2*pi*1e9*(t-tc-delta_t3));%传感器3接收到PD信号%plot(y0)%holdon%plot(y1)%plot(y2)%plot(y3)tic;%程序运行计时%[~,R]=gccphat(y1',y0');[~,YI]=max(R);t10=(YI-length(t))*diff(t(1:2));%GCCPHAT计算时延1%[~,R]=gccphat(y2',y0');[~,YI]=max(R);t20=(YI-length(t))*diff(t(1:2));%GCCPHAT计算时延2%[~,R]=gccphat(y3',y0');[~,YI]=max(R);t30=(YI-length(t))*diff(t(1:2));%GCCPHAT计算时延3[~,YI]=max(xcorr(y1,y0));t10=(YI-length(t))*diff(t(1:2));%GCCPHAT计算时延1[~,YI]=max(xcorr(y2,y0));t20=(YI-length(t))*diff(t(1:2));%GCCPHAT计算时延2[~,YI]=max(xcorr(y3,y0));t30=(YI-length(t))*diff(t(1:2));%GCCPHAT计算时延3t10=delta_t1-delta_t0;t20=delta_t2-delta_t0;t30=delta_t3-delta_t0;%Lx=1;Ly=1;Lz=0.8;%变压器长宽高%Sensor_xyz0=[0,0.2,0.4]';Sensor_xyz1=[0,0.6,0.3]';Sensor_xyz2=[0.3,1,0.2]';Sensor_xyz3=[0.6,1,0.5]';%四个传感器位置%Source=[0.61,0.32,0.35]';%PD位置%t10=0.036;t20=0.1;t30=0.053;%PSO初始化error=1e-6;%允许误差MaxNum=50;%粒子最大迭代次数nvars=3;%目标函数的自变量个数particlesize=30;%粒子群规模c1=2;%每个粒子的个体学习因子，也称为加速常数1c2=2;%每个粒子的社会学习因子，也称为加速常数2w=0.6;%惯性因子a=1;%约束因子，控制速度的权重vmax=0.8;%粒子的最大飞翔速度%x=3*rand(particlesize,nvars);%粒子的所在位置x(:,1)=Lx*rand(particlesize,1);%初始化x(:,2)=Ly*rand(particlesize,1);x(:,3)=Lz*rand(particlesize,1);%x=2*ones(particlesize,nvars);v=2*rand(particlesize,nvars);%粒子的飞翔速度%%ObjectiveFunction%用inline定义适应度函数以便将子函数文件与主程序文件放在一起%目标函数为：y=1+(2.1*(1-x+2*x.^2).*exp(-x.^2/2))注意不同于原始目标函数%inline命令定义适应度函数如下%fitness=inline('1./(1+(2.1*(1-x+2*x.^2).*exp(-x.^2/2)))','x');%inline定义的适应度函数会使程序运行速度大大降低%%Particleswarmoptimizationf=zeros(particlesize,1);fori=1:particlesize%forj=1:nvars%f(i)=fitness(x(i,j));%endf(i)=Fitness1(x(i,:),Sensor_xyz0',Sensor_xyz1',Sensor_xyz2',Sensor_xyz3',t10,t20,t30,c);endpersonalbest_x=x;personalbest_faval=f;[globalbest_faval,YI]=min(personalbest_faval);globalbest_x=personalbest_x(YI,:);k=1;Trace_error=[];Trace_pos=[];whilek<=MaxNumfori=1:particlesize%forj=1:nvars%f(i)=fitness(x(i,j));%endf(i)=Fitness1(x(i,:),Sensor_xyz0',Sensor_xyz1',Sensor_xyz2',Sensor_xyz3',t10,t20,t30,c);iff(i)<personalbest_faval(i)%判断当前位置是否是历史上的最佳位置personalbest_faval(i)=f(i);personalbest_x(i,:)=x(i,:);endend[globalbest_faval,i]=min(personalbest_faval);globalbest_x=personalbest_x(i,:);fori=1:particlesize%更新粒子群里的每个个体的最新位置v(i,:)=w*v(i,:)+c1*rand*(personalbest_x(i,:)-x(i,:))+c2*rand*(globalbest_x-x(i,:));forj=1:nvars%判断每个粒子的飞翔速度是否超过了最大飞翔速度ifv(i,j)>vmaxv(i,j)=vmax;elseifv(i,j)<-vmaxv(i,j)=-vmax;endendx(i,:)=x(i,:)+a*v(i,:);ifx(i,1)>Lxx(i,1)=Lx*rand;endifx(i,2)>Lyx(i,2)=Ly*rand;endifx(i,3)>Lzx(i,3)=Lz*rand;endendTrace_error=[Trace_error,globalbest_faval];Trace_pos=[Trace_pos;globalbest_x];ifabs(globalbest_faval)<errorbreak;endk=k+1;end%%PostProcessingsubplot(2,1,1)plot(Trace_error,'k--','LineWidth',2)gridonsubplot(2,1,2)plot3(Trace_pos(:,1),Trace_pos(:,2),Trace_pos(:,3),'--','LineWidth',2)holdonplot3(Trace_pos(end,1),Trace_pos(end,2),Trace_pos(end,3),'r*')xlim([0,Lx])ylim([0,Ly])zlim([0,Lz])plot3(Source(1),Source(2),Source(3),'ok')gridonValue2=globalbest_x;Value2=num2str(Value2);disp(strcat('Thecorrespondingcoordinate=',Value2))%Value1=1/globalbest_faval-1;%Value1=num2str(Value1);%%PostProcessing%%strcat指令可以实现字符的组合输出%disp(strcat('Themaximumvalue=','',Value1))%%输出最大值所在的横坐标位置%Value2=globalbest_x;%Value2=num2str(Value2);%disp(strcat('Thecorrespondingcoordinate=',Value2))%x=-5:0.01:5;%y=2.1*(1-x+2*x.^2).*exp(-x.^2/2);%plot(x,y,'m-','LineWidth',3)%holdon%plot(globalbest_x,1/globalbest_faval-1,'kp','LineWidth',4)%legend('ObjectiveFunction','Themaximumvalue')%xlabel('x')%ylabel('y')%gridontoc;适应度函数：%%PDLocationPSO1适应度函数functionF1=Fitness1(Source,Sensor_xyz0,Sensor_xyz1,Sensor_xyz2,Sensor_xyz3,t10,t2