基于遗传算法的Matlab 16阵元天线优化
1. 设计要求
题目一 ：阵列天线综合
利用Matlab编制遗传算法或粒子群算法程序，并实现对间距为半波长均匀直线阵综合，指标如下：
1、阵元数：16元
2、副瓣电平：<-30dB
3、增益：>11dB
4、要求撰写设计报告，内容包括：所采用的算法基本原理，目标函数的设计，各个参数的设置源代码，仿真结果（增益方向图），参考文献。
2. 遗传算法
2.1 遗传算法的生物学基础
遗传算法的生物学基础自然选择学说认为适者生存，生物要存活下去，就必须进行生存斗争。生存斗争包括种内斗争、种间斗争以及生物跟环境之间的斗争三个方面。在生存斗争中，具有有利变异的个体容易存活下来，并且有更多的机会将有利变异传给后代，具有不利变异的个体就容易被淘汰，产生后代的机会就少得多。因此，凡是在生存斗争中获胜的个体都是对环境适应性比较强的个体。
2.2 遗传算法介绍
遗传算法(GA)是模拟生物在自然环境下的遗传和进化过程而形成的一种自适应全局优化概率搜索方法。其采纳了自然进化模型，从代表问题可能潜在解集的一个种群开始，种群由经过基因编码的一定数目的个体组成。每个个体实际上是染色体带有特征的实体。初始种群产生后，按照适者生存和优胜劣汰的原理，逐代演化产生出越来越好的解。在每一代，概据问题域中个体的适应度大小挑选个体，并借助遗传算子进行组合交叉和主客观变异，产生出代表新的解集的种群。这一过程循环执行，直到满足优化准则为止。最后，末代个体经解码，生成近似最优解。基于种群进化机制的遗传算法如同自然界进化一样，后生代种群比前生代更加适应于环境，通过逐代进化，逼近最优解。
2.3 算法流程
通过随机方式产生若干由确定长度（长度与待求解问题的精度有关）编码的初始群体。通过适应度函数对每个个体进行评价，选择适应度值高的个体参与遗传操作，适应度低的个体被淘汰。经遗传操作（复制、交叉、变异）的个体集合形成新一代种群，直到满足停止准则（进化代数GEN>=?）。将后代中变现最好的个体作为遗传算法的执行结果。其中，GEN是当前代数；M是种群规模，i代表种群数量。
图2.1 遗传算法流程图
2.4 选择
指个体被选中并遗传到下一代群体中的概率与该个体的适应度大小成正比。
式中 \(\ p_{i}\) ——个体i被选中的概率。
\(\ f_{i}\)——个体i的适应度。
\(\sum f_{i}\)——群体的累加适应度。
- 显然，个体适应度愈高，被选中的概率愈大。但是，适应度小的个体也有可能被选中，以便增加下一代群体的多样性。执行概率选择的手段是轮盘选择。
- 个体被选中的概率取决于个体的相对适应度：从统计意义讲，适应度大的个体，其刻度长，被选中的可能性大；反之，适应度小的个体被选中的可能性小，但有时也会被“破格”选中。
2.5 交叉
Ⅰ. 对群体中的个体进行两两随机配对。
若群体大小为M，则共有M/2 对相互配对的个体组。
Ⅱ. 每一对相互配对的个体，随机设置某一基因座之后的位置为交叉点。
若染色体的长度为l ，则共有l-1个可能的交叉点位置。
Ⅲ. 对每一对相互配对的个体，依设定的交叉概率pc在其交叉点处相互交换两个个体的部分染色体，从而产生出两个新的个体。
2.6 变异
对于基本遗传算法中用二进制编码符号串所表示的个体，若需要进行变异操作的某一基因座上的原有基因值为0，则变异操作将该基因值变为1，反之，若原有基因值为1，则变异操作将其变为0。
基本位变异因子的具体执行过程是：
Ⅰ. 对个体的每一个基因座，依变异概率pm指定其为变异点。
Ⅱ. 对每一个指定的变异点，对其基因值做取反运算或用其它等位基因值来代替， 从而产生出一个新的个体。
3. 阵列天线原理
直线阵:由N个相同的单元天线等间距地排列在一条直线上构成。
均匀直线阵:若各单元上的馈电电流大小相同，而相位沿线均匀递增或递减，这样的直线阵称为均匀直线阵。
N元直线阵的阵因子函数：
以阵列中第一个单元天线作为相位基准：
第二个单元天线：辐射场的幅值差异不计，相位差为
则第N个单元天线的相位差为：
N个相似元阵元在p处的辐射场叠加，表示为
整理可得远区辐射场表达式为：
4. matlab 程序设计
% matlab程序设计
clear all;
close all;
clc;
NP=50; %个体数
Pc = 0.8; %选择概率
Pm = 0.05; %变异概率
f=30000000;%信号频率
c=30000000;
lamda = c/f;%波长
d = lamda/2;%间距
beta = 2*pi/lamda;%波数
seta0 = 0;%波束方向
G = 100; %最大遗传代数
L = 72; %编码数
NL = 16; %阵元数
NN = 1800; %划分刻度
E0 = 1; %电压电流
%%%生成初始种群
f = randi(NP,L); %产生一个L*NP的随机矩阵，正态分布
%初始种群的分布
maxE = 100;
for i = 1:NP
for j = 1:L
if f(i,j)==1
plot(i,j)
hold on
end
end
end
%%%%遗传算法循环
for k = 1:G
k
%解码，个体实际数
% 计算峰值旁瓣比，即适应度
for i = 1:NP
F(i) = -MSLL3(d,beta,NN,NL,seta0,f(i,:),L,maxE); % 取成正的，下面方便取最大值并且处理
end% 取得NP个个体的旁瓣电平
maxF = max(F)
minF = min(F)
rr = find(F == maxF);
fBest = f(rr(1),:);
Fin = (F-minF)/(maxF-minF); %归一化
%rr(1,1)
%fBest
%maxFit
%%%基于轮赌盘的选择操作
sum_Fin = sum(Fin);
fitvalue = Fin./sum_Fin;
fitvalue = cumsum(fitvalue);%cumsum(A) 返回一个和A同行同列的矩阵，
%矩阵中第m行第n列元素是A中第1行到第m
%行的所有第n列元素的累加和
ms = sort(rand(NP,1));
fiti = 1;
newi = 1;
while newi <= NP
if(ms(newi)<fitvalue(fiti))
nf(newi,:) = f(fiti,:);
newi = newi+1;
else
fiti = fiti+1;
end
end
%概率交叉
for i = 1:2:NP
P=rand;
if P<Pc
q = randi(1,L);%生成1*L随机矩阵，0或1
for j = 1:L
if q(j) == 1
temp = nf(i+1,j);
nf(i+1,j) = nf(i,j);
nf(i,j) = temp;
end
end
end
end
%%%概率变异
for m = 1:NP
for n = 1:L
r = rand(1,1);
if r < Pm
nf(m,n) = ~nf(m,n);
end
end
end
f = nf;
f(1,:) = fBest;
trace(k) = maxF;
end
figure
plot(trace)
xlabel('迭代次数')
ylabel('目标函数值')
title('适应度进化曲线')
save fBestn.mat fBest
save L.mat L
save E.mat E0
%%%计算最大旁瓣%%%
function MSLL3 = MSLL3(d,beta,NN,NL,seta0,f,L,maxE)
seta = linspace(-pi,pi,NN);
alfa = (2*pi*(f(1:L/2)*(2.^(0:L/2-1)'))/(2^(L/2)-1)); %解码
x = maxE/8*f(L/2+1:L)*(2.^(0:L/2-1)')/(2^(L/2)-1); %解码
I1 = linspace(maxE-(NL/2)*x,maxE,NL/2);
I2 = linspace(maxE,maxE-(NL/2*x),NL/2);
I = [I1,I2]-x;
for m = 1:NN
fai = beta*d*([0:NL-1]+1-(NL+1)/2)*cos(seta(m)-seta0)-[0:NL-1]*alfa;
Fit(m) = abs(sum(I*exp(sqrt(-1)*(fai)')));
end
%plot(alfa,fit)
%hold on
FdB = 20*log10(Fit/max(Fit));
%那如果需要找到其中满足一定条件的元素索引
%%%取旁瓣电平，即最大值之外的次大值
FdB = FdB-3;
for i = 1:2
maxFdB = max(FdB);
rr = find(FdB == maxFdB); %find()函数的功能是找到向量或者矩阵中不为0的元素，
mm = rr(1);
tu_up = 0;
while(FdB(mm+tu_up) >= FdB(mm+tu_up+1))
tu_up = tu_up+1;
if mm+tu_up+1>1800
mm = -tu_up+1;
end
end
mm = rr(1);
tu_down = 0;
while (FdB(mm-tu_down) >= FdB(mm-tu_down-1))
tu_down = tu_down+1;
if mm-tu_down-1<1
mm = 1800+tu_down;
end
end
FdB(mm-tu_down:mm+tu_up) = -50;
end
MSLL3 = max(FdB);
end
clear all;
close all;
clc;
load('fBestn.mat');
load('L.mat')
load('E.mat')
fbest = 0;
NP=50; %个体数
Pc = 0.8; %选择概率
Pm = 0.01; %变异概率
f=30000000;%信号频率
c=30000000;
lamda = c/f;%波长
d = lamda/2;%间距
beta = 2*pi/lamda;%波数
seta0 = 0;%波束方向
G = 200; %最大遗传代数
NL = 16; %阵元数
NN = 1800; %划分刻度
e = 1;
maxE = 8;
seta = linspace(-pi,pi,NN);
%for m = 1:NN
%alfa = sum(2*pi*fBest*(2.^(0:L-1)')/(2^L-1)); %解码
%fai = alfa+beta*d.*cos(seta(m));
%F1(m) = abs(sin(NL*(fai./2))./sin(fai./2));
%end
alfa = (2*pi*(fBest(1:L/2)*(2.^(0:L/2-1)'))/(2^(L/2)-1)); %解码
x = E0*fBest(L/2+1:L)*(2.^(0:L/2-1)')/(2^(L/2)-1); %解码
I1 = linspace(maxE-(NL/2)*x,maxE,NL/2);
I2 = linspace(maxE,maxE-(NL/2*x),NL/2);
I = [I1,I2]-x;
for m = 1:NN
fai = beta*d*([0:NL-1]+1-(NL+1)/2)*cos(seta(m)-seta0)-[0:NL-1]*alfa;
F1(m) = abs(sum(I*exp(sqrt(-1)*(fai)')));
end
F1 = abs(F1);
D1 = 0;
for m = 1:NL-1;
D1 = (NL-m)/(m*beta*d)*sin(m*beta*d)*cos(m*alfa);
end
D = 1/(1/NL+(2/NL^2)*D1)
%plot(alfa,fit)
%hold on
plot(seta,F1);
FdB = 20*log10(F1/max(F1));
% 归一化方向图
FdB = FdB-3;
figure
plot(seta*180/pi,FdB)
xlabel('\theta/(度)')
ylabel('阵列增益/dB')
axis([0 180 -50 0]);%可为x轴和y轴设置一个极限范围
grid on
%%%求最大副瓣电平
for i = 1:2
maxFdB = max(FdB);
rr = find(FdB == maxFdB); %find()函数的功能是找到向量或者矩阵中不为0的元素，
maxFdB = max(FdB);
rr = find(FdB == maxFdB); %find()函数的功能是找到向量或者矩阵中不为0的元素，
mm = rr(1);
tu_up = 0;
while(FdB(mm+tu_up) >= FdB(mm+tu_up+1))
tu_up = tu_up+1;
if mm+tu_up+1>1800
mm = -tu_up+1;
end
end
mm = rr(1);
tu_down = 0;
while (FdB(mm-tu_down) >= FdB(mm-tu_down-1))
tu_down = tu_down+1;
if mm-tu_down-1<1
mm = 1800+tu_down;
end
end
FdB(mm-tu_down:mm+tu_up) = -50;
end
figure(3)
plot(seta*180/pi,FdB);
axis([-180,180,-30,10])
MSLL = max(FdB) %那如果需要找到其中满足一定条件的元素索引 找到最大旁瓣电平
5. 仿真结果
在Matlab中运行程序可得迭代过程中的适应度变化曲线。由下图可知程序设计的较为合理，能较好的满足实验需求。
图5.1 适应度变化曲线
运行程序可得归一化后的方向图如下图所示:
运行程序找出最大副瓣结果如下所示，满足题意
MSLL = -30.0622
注：本文或许会有很多不足之处，希望取其精华去其糟粕，仅供参考，各位小伙伴若有其他更好的方案，欢迎批评指正，望共同进步。
