基于遗传/粒子群算法的多类型变量&多优化目标下的优化仿真实践 - 以波束形成优化为例 博文对应的代码

%% 多变量、多优化目标下的波束形成优化 %% %多变量是指： 同时包含幅值、相位以及阵元相对位置关系三种变量 %多优化目标是指： 使波束指向特定的角度、波束的峰值旁瓣比越高越好 %分别使用粒子群算法和遗传算法进行仿真，甚至：对于不同的优化变量可以使用不同的优化算法？：幅值用粒子群去优化，相位用遗传算法去优化？% %本代码使用粒子群算法对全部的变量进行优化，并同时优化角度和峰值旁瓣比。 %author: https://blog.csdn.net/xhblair?spm=1010.2135.3001.5343 clear; close all; clc; fc = 77e9; c = 3e8; lambda = c/fc; %------阵列参数预设------% ElementNum = 8; %线阵，阵元个数 I_range = [0 10]; %各阵元馈电幅值的可选范围 phi_range = [-180 180]; %各阵元馈电相位的可选范围 arrayAperture = [0 10*lambda]; %阵列孔径：各阵元在此限定孔径下排布(为稀布阵) dmin = 0.5*lambda; %【允许的两两阵元之间的最短间隔】 Fov = (-90:1:90); %阵列水平向的方向图范围(我们优化的范围) %需求的波束方向 XitaNeed = 10; %我们希望波束指向该方向 DxitaMax = max(abs(Fov(1)-XitaNeed),abs(Fov(end) - XitaNeed)); %可能的最大的波束差，用这个值来将波束差变成(0,1)范围内 %------粒子群算法参数预设------% Num_partical = 20; %预设粒子群数量 length_partical = 3*ElementNum; %每个粒子的长度(1:8对应每个粒子的幅度设置、9:16对应相位、17:24对应位置) Vrange_I = [-1 1]; %预设粒子(各阵元)幅值变化速度的限定范围 Vrange_phi = [-20 20]; %预设粒子(各阵元)相位变化速度的限定范围 Vrange_site = [-1*lambda 1*lambda]; %预设粒子(各阵元)位置变化速度的限定范围 % w_I = [0.4 0.8]; %幅值惯性权重范围，本仿真中该值随迭代次数线性减小。 % C1_I = 0.5; %个体学习因子 % C2_I = 0.5; %群体学习因子 % w_phi = [0.4 0.8]; %相位惯性权重范围，本仿真中该值随迭代次数线性减小。 % C1_phi = 0.5; %个体学习因子 % C2_phi = 0.5; %群体学习因子 % w_site = [0.4 0.8]; %位置惯性权重范围，本仿真中该值随迭代次数线性减小。 % C1_site = 0.5; %个体学习因子 % C2_site = 0.5; %群体学习因子 w_I = [0.6 1.2]; %幅值惯性权重范围，本仿真中该值随迭代次数线性减小。 C1_I = 0.4; %个体学习因子 C2_I = 0.6; %群体学习因子 w_phi = [0.6 1.2]; %相位惯性权重范围，本仿真中该值随迭代次数线性减小。 C1_phi = 0.4; %个体学习因子 C2_phi = 0.6; %群体学习因子 w_site = [0.6 1.2]; %位置惯性权重范围，本仿真中该值随迭代次数线性减小。 C1_site = 0.4; %个体学习因子 C2_site = 0.6; %群体学习因子 Num_iteration = 200; %迭代次数(迭代终止条件) %------在预设参数下对各个粒子进行初始化------% PositionP = zeros(Num_partical,length_partical); VelocityP = zeros(Num_partical,length_partical); %初始化阵元的幅值及幅值速度 PositionP(:,1:ElementNum) = rand(Num_partical,ElementNum).*I_range(2); VelocityP(:,1:ElementNum) = rand(Num_partical,ElementNum).*(Vrange_I(2)- Vrange_I(1)) + Vrange_I(1); %初始化阵元的相位及相位速度 PositionP(:,ElementNum+1:2*ElementNum) = rand(Num_partical,ElementNum).*(phi_range(2)-phi_range(1)) + phi_range(1); VelocityP(:,ElementNum+1:2*ElementNum) = rand(Num_partical,ElementNum).*(Vrange_phi(2)- Vrange_phi(1)) + Vrange_phi(1); %初始化阵元的相对位置关系 Site_partical = sort( (rand(Num_partical,ElementNum).*(arrayAperture(2)-arrayAperture(1)) + arrayAperture(1)),2); %对粒子位置进行优化： 将全部阵元的位置往第一个阵元对齐：减去第一个阵元的位置值, 并按顺序判断和适当改变后续阵元的位置。 for ii = 1:Num_partical Site_partical(ii,:) = Site_partical(ii,:) - Site_partical(ii,1); %平移使得第一个阵元位置为0。 for jj = 1:ElementNum-1 if (Site_partical(ii,jj+1) - Site_partical(ii,jj)) < dmin Site_partical(ii,jj+1) = Site_partical(ii,jj) + dmin; %我们允许最后一个阵元略微超过阵列孔径，但是这种可能性很小。 end end if Site_partical(ii,end) < arrayAperture(2) %并尽量保证阵列的孔径为预设孔径（如果最后一个阵元的位置小于预设孔径，将它放在预设孔径处） Site_partical(ii,end) = arrayAperture(2); end end PositionP(:,2*ElementNum+1:3*ElementNum) = Site_partical; %初始化速度 VelocityP(:,2*ElementNum+1:3*ElementNum) = rand(Num_partical,ElementNum).*(Vrange_site(2)- Vrange_site(1)) + Vrange_site(1); %------预设一些用于可视化的变量，并定义适应度函数------% BestFit_partical_Fitness = zeros(Num_partical,1); %装载各粒子的历史最好适应度 BestFit_partical_partical = zeros(Num_partical,3*ElementNum); %装载各粒子最好适应度下的粒子 BestFit_partical_arrayPattern = zeros(Num_partical,length(Fov)); %装载各粒子最好适应度下对应的波束方向图 BestFit_global_Fitness = zeros(1,1); %装载粒子群的历史最好适应度 BestFit_global_partical = zeros(1,3*ElementNum); %装载对应的粒子 BestFit_global_arrayPattern = zeros(1,length(Fov)); %装载对应的波束方向图 BestFit_everyIteration_Fitness = zeros(Num_iteration+1,1); %装载历次迭代的最好适应度（第一行为初始化下的值） BestFit_everyIteration_partical = zeros(Num_iteration+1,3*ElementNum); %装载历次迭代的最好粒子 BestFit_everyIteration_arrayPattern = zeros(Num_iteration+1,length(Fov)); %装载历次迭代的最好波束方向图 %------进行迭代------% %可以看看历次迭代下的各粒子适应度的变化 figure(1); hold on; for ii = 1:(Num_iteration+1) %因为把初始化下的适应度等计算也放在这里面了。 %需要更新惯性权重 Wuse_I = w_I(2) - (w_I(2)-w_I(1))/Num_iteration*ii; Wuse_phi = w_phi(2) - (w_phi(2)-w_phi(1))/Num_iteration*ii; Wuse_site = w_site(2) - (w_site(2)-w_site(1))/Num_iteration*ii; %首先计算波束方向图 [arrayPattern] = CalArrayPattern(PositionP,Fov,Num_partical,ElementNum,lambda); %计算适应度值，1.评估主波束方向是否是需求的波束方向(越接近适应度越高) 2.评估峰值旁瓣比(越大适应度越高) %我们需要将这两个值归一化(如果不归一化很难直接加权值并相加)，并合并得到一个适应度值。 %计算各粒子与需求波束方向的差值(通过找到值最大处对应的方向与需求方向的做差得到) [BeamGap] = CalBeamGap(arrayPattern,Fov,XitaNeed,Num_partical); %计算各粒子下的峰值旁瓣比 [PSLR] = CalPSLR(arrayPattern,Num_partical); %计算适应度 [Fitness] = ClaFitness(BeamGap,PSLR,Num_partical,DxitaMax); %进行装载，并判断是否为Num_iteration+1 如果是，不需要进行后续处理，此时装载的就是第Num_iteration次迭代下的结果。 [BestFit_everyIteration_Fitness(ii,1),Index] = max(Fitness); BestFit_everyIteration_partical(ii,:) = PositionP(Index,:); BestF