Matlab粒子群算法(PSO)進行空間-變壓器局部放電定位

%%ParticleSwarmOptimization(PSO),粒子群優(yōu)化算法計算局放源位置%%%%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;%粒子最大迭代次數(shù)nvars=3;%目標函數(shù)的自變量個數(shù)particlesize=30;%粒子群規(guī)模c1=2;%每個粒子的個體學習因子，也稱為加速常數(shù)1c2=2;%每個粒子的社會學習因子，也稱為加速常數(shù)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定義適應度函數(shù)以便將子函數(shù)文件與主程序文件放在一起%目標函數(shù)為：y=1+(2.1*(1-x+2*x.^2).*exp(-x.^2/2))注意不同于原始目標函數(shù)%inline命令定義適應度函數(shù)如下%fitness=inline('1./(1+(2.1*(1-x+2*x.^2).*exp(-x.^2/2)))','x');%inline定義的適應度函數(shù)會使程序運行速度大大降低%%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指令可以實現(xiàn)字符的組合輸出%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;適應度函數(shù)：%%PDLocationPSO1適應度函數(shù)functionF1=Fitness1(Source,Sensor_xyz0,Sensor_xyz1,Sensor_xyz2,Sensor_xyz3,t10,t2