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
注：本文或許會有很多不足之處，希望取其精華去其糟粕，僅供參考，各位小伙伴若有其他更好的方案，歡迎批評指正，望共同進步。