51基于布谷鸟算法的导航波束形成问题

Halcom

基于布谷鸟算法的导航波束形成问题：
百度网盘链接：
链接：http://pan.baidu.com/s/1hsmpdJ2

具体链接在halcom.cn论坛，联系人QQ：3283892722
该论坛是一个学习交流平台，我会逐一的和大家分享学习。
欢迎大家录制视频，你可在论坛进行打赏分享。
视频专用播放器：http://halcom.cn/forum.php?mod=viewthread&tid=258&extra=page%3D1\

具体的代码如下：

1. % 布谷鸟函数优化分析
   
2. % cuckoo_search -- CS算法
   
3. %% 清空环境
   
4. clc          % 清屏
   
5. clear all;   % 删除workplace变量
   
6. close all;   % 关掉显示图形窗口
   
7. warning off; % 消除警告
   
8. tic
   
9. %% 参数初始化
   
10. maxiter = 20;    % 迭代次数
    
11. sizepop = 20;    % 种群规模
    
12. pa=0.25;         % 发现新巢的概率Discovery rate of alien eggs/solutions
    
13. %        M   K  phi
    
14. popmin = [2, 0, 331, 31, 91, 151, 211, 271];   % x 最小值
    
15. popmax = [6, 1, 390, 90, 150, 210, 270, 330];  % x 最大值
    
16. nvar = 8;  % 未知量个数
    
17. 
    
18. % 导航-波束形成问题
    
19. theta_signal=50;    % 俯仰角
    
20. phi_signal=200;     % 方位角
    
21. 
    
22. % 产生初始粒子
    
23. for i=1:sizepop
    
24.     nest(i,:)= popmin+(popmax-popmin).*rand(size(popmin)) ;
    
25.     fitness(i)= fun( [nest(i,1), nest(i,2)*50000, nest(i,3:nvar)], theta_signal,phi_signal );  % 适应度函数值--目标函数值--求解全局的极大值
    
26. end
    
27. % 找当前最好的个体
    
28. [fmax,bestindex]=max(fitness);
    
29. bestnest=nest(bestindex,:);   % 全局最佳
    
30. %% 寻优
    
31. for i=1:maxiter
    
32.     disp(['Iteration ' num2str(i)]);
    
33.     % 粒子位置更新
    
34.     new_nest=get_cuckoos(nest,bestnest,popmin,popmax); % Levy flights 粒子位置更新
    
35.     % 比较更新
    
36.     [fnew,best,nest,fitness]=get_best_nest(nest,new_nest,fitness, theta_signal,phi_signal);
    
37.     % 发现新巢、随机粒子的位置
    
38.     new_nest=empty_nests(nest,popmin,popmax,pa) ;
    
39.     % 继续评估新巢的适应度值
    
40.     [fnew,best,nest,fitness]=get_best_nest(nest,new_nest,fitness, theta_signal,phi_signal);
    
41.     % 寻找全局最优解  
    
42.     if fnew>fmax,
    
43.         fmax=fnew;
    
44.         bestnest=best;
    
45.     end
    
46.    
    
47.     CS_yy(i)=fmax;   % 保存极大值
    
48. end  % 结束迭代循环
    
49. CS_fitness_iter = CS_yy;
    
50. time = toc
    
51. save CS_fitness_iter.mat
    
52. %% 结果分析
    
53. plot(CS_yy,'bo-','Linewidth',2)
    
54. title(['适应度曲线  ' '终止代数＝' num2str(maxiter)]);
    
55. grid on
    
56. xlabel('进化代数');ylabel('适应度');
    
57. % 结果输出
    
58. disp('最优解')
    
59. zbest_result = fix( [ bestnest(1),bestnest(2)*50000, bestnest(3:nvar)] );
    
60. disp(zbest_result)   % 最佳个体值
    
61. fprintf('\n')
    
62. 
    
63. %% 结果
    
64. tic
    
65. [fitness_zbest, F]= fun(zbest_result, theta_signal, phi_signal);
    
66. time1 = toc
    
67. figure,
    
68. plot(0:359, F(:,50)); grid on   %截面图：
    
69. xlabel('方位角');
    
70. ylabel('增益(dB)');
    
71. title('俯仰角50度时截面')

复制代码

寻找新的解：

1. function new_nest=empty_nests(nest,Lb,Ub,pa)
   
2. % A fraction of worse nests are discovered with a probability pa
   
3. n=size(nest,1);
   
4. % Discovered or not -- a status vector
   
5. K=rand(size(nest))>pa;
   
6. %% New solution by biased/selective random walks
   
7. stepsize=rand*(nest(randperm(n),:)-nest(randperm(n),:));
   
8. new_nest= nest+stepsize.*K;
   
9. for j=1:size(new_nest,1)
   
10. 
    
11.     for k=1:length(nest(j,:))
    
12.        if new_nest(j,k)>Ub(k)
    
13.            new_nest(j,k)=Ub(k);
    
14.        end
    
15.        if new_nest(j,k)<Lb(k)
    
16.            new_nest(j,k)=Lb(k);
    
17.        end
    
18.     end
    
19.   
    
20. end

复制代码

寻找当前最优的解：

1. function [fmax,best,nest,fitness]=get_best_nest(nest,newnest,fitness, theta_signal,phi_signal)
   
2. % Evaluating all new solutions
   
3. for j=1:size(nest,1),
   
4. %     fnew=Ackley_Fun(newnest(j,:));  % 适应度值
   
5.     fnew= fun( [newnest(j,1), newnest(j,2)*50000,newnest(j, 3:end) ], theta_signal,phi_signal );  % 适应度函数值--目标函数值--求解全局的极大值
   
6.    
   
7.     if fnew>=fitness(j),
   
8.        fitness(j)=fnew;
   
9.        nest(j,:)=newnest(j,:);
   
10.     end
    
11. end
    
12. % Find the current best
    
13. [fmax,K]=max(fitness) ;
    
14. best=nest(K,:);

复制代码

位置更新：

1. % 位置更新
   
2. function nest=get_cuckoos(nest,best,Lb,Ub)
   
3. % Levy flights
   
4. n=size(nest,1);
   
5. % Levy exponent and coefficient
   
6. beta=3/2;
   
7. sigma=(gamma(1+beta)*sin(pi*beta/2)/(gamma((1+beta)/2)*beta*2^((beta-1)/2)))^(1/beta);
   
8. 
   
9. for j=1:n,
   
10.     s=nest(j,:);
    
11.     % implementing Levy flights
    
12.     % For standard random walks, use step=1;
    
13.     %% Levy flights by Mantegna's algorithm
    
14.     u=randn(size(s))*sigma;
    
15.     v=randn(size(s));
    
16.     step=u./abs(v).^(1/beta);  
    
17.     stepsize=0.01*step.*(s-best);
    
18.     nest(j,:) = s+stepsize.*randn(size(s));
    
19.    % Apply simple bounds/limits
    
20.    for k=1:length(nest(j,:))
    
21.        if nest(j,k)>Ub(k)
    
22.            nest(j,k)=Ub(k);
    
23.        end
    
24.        if nest(j,k)<Lb(k)
    
25.            nest(j,k)=Lb(k);
    
26.        end
    
27.    end
    
28.    
    
29. end
    

复制代码

适应度函数：

1. function [fitness, F]= fun(x, theta_signal, phi_signal)
   
2. % 输入：
   
3. %      theta_signal=50;    % 俯仰角
   
4. %      phi_signal=200;     % 方位角
   
5. %      x:
   
6. %         N=7;      % 阵元数目
   
7. %         M=6;      % 延迟数目
   
8. %         K=5000;   % 快拍数
   
9. % 输出：
   
10. %　　　fitness：适应度值，主瓣-旁瓣的差值最大
    
11. %      F：输出的beampattern
    
12. 
    
13. % x=fix(x); % 取整
    
14. N = 7;                   % 阵元数目
    
15. M = fix(x(1));                   % 延迟数目
    
16. K = fix(x(2));                   % 快拍数
    
17. m=[0:M-1]';               
    
18. 
    
19. c=3*10^8;                   % 光速
    
20. f0=1.26852*10^9;            % 信号中心频率
    
21. fs=5*10^9;                  % 采样频率
    
22. Ts=1/fs;                    % 单位时延
    
23. t = 0:Ts:(K*Ts-Ts);   
    
24. wl=c/f0;                    % 信号波长
    
25. r=wl/2;                     % 圆阵半径
    
26. % phi1=[0, 60, 120, 180, 240, 300];  % 圆周上阵元角度位置 6围1圆阵
    
27. % phi1 = 0:360/(N-1):360-360/(N-1);
    
28. phi1 = fix( x(3:end) );
    
29. P=[zeros(2,1), r*[cos(phi1*pi/180);sin(phi1*pi/180)]]; % 位置
    
30. 
    
31. % theta_signal=50;    % 俯仰角
    
32. % phi_signal=200;     % 方位角
    
33. 
    
34. theta_inter=40;     % 干扰1
    
35. phi_inter=50;
    
36. 
    
37. theta_inter2=60;    % 干扰2
    
38. phi_inter2=250;
    
39. 
    
40. theta_inter3=50;    % 干扰3
    
41. phi_inter3=60;
    
42. 
    
43. theta_inter4=30;    % 干扰4
    
44. phi_inter4=120;
    
45. 
    
46. theta_inter5=30;    % 干扰5
    
47. phi_inter5=80;
    
48. 
    
49. waven_signal=[sin(theta_signal*pi/180)*cos(phi_signal*pi/180), sin(theta_signal*pi/180)*sin(phi_signal*pi/180) ];  % 信号波数（wavenumber）
    
50. 
    
51. waven_inter=[sin(theta_inter*pi/180)*cos(phi_inter*pi/180), sin(theta_inter*pi/180)*sin(phi_inter*pi/180)];
    
52. waven_inter2=[sin(theta_inter2*pi/180)*cos(phi_inter2*pi/180), sin(theta_inter2*pi/180)*sin(phi_inter2*pi/180)];
    
53. waven_inter3=[sin(theta_inter3*pi/180)*cos(phi_inter3*pi/180), sin(theta_inter3*pi/180)*sin(phi_inter3*pi/180)];
    
54. waven_inter4=[sin(theta_inter4*pi/180)*cos(phi_inter4*pi/180), sin(theta_inter4*pi/180)*sin(phi_inter4*pi/180)];
    
55. waven_inter5=[sin(theta_inter5*pi/180)*cos(phi_inter5*pi/180), sin(theta_inter5*pi/180)*sin(phi_inter5*pi/180)];
    
56. 
    
57. snr1=-20;                   % 信噪比
    
58. 
    
59. inr1=30;                    % 干噪比
    
60. inr2=30;                    
    
61. inr3=30;
    
62. inr4=30;
    
63. inr5=30;
    
64. 
    
65. s_a0=exp(j*2*pi/wl*(waven_signal*P)');    % 期望信号1
    
66. s_a1=exp(j*2*pi/wl*(waven_inter*P)');     % 干扰1
    
67. s_a2=exp(j*2*pi/wl*(waven_inter2*P)');    % 干扰2
    
68. s_a4=exp(j*2*pi/wl*(waven_inter3*P)');
    
69. s_a5=exp(j*2*pi/wl*(waven_inter4*P)');
    
70. s_a6=exp(j*2*pi/wl*(waven_inter5*P)');
    
71. 
    
72. t_a=exp(j*2*pi*f0*Ts*m);                          
    
73. 
    
74. steer_signal=kron(t_a,s_a0);           % Kronecker积，期望及干扰导向矢量
    
75. 
    
76. steer_jamming1=kron(t_a,s_a1);
    
77. steer_jamming2=kron(t_a,s_a2);
    
78. steer_jamming3=kron(t_a,s_a4);
    
79. steer_jamming4=kron(t_a,s_a5);
    
80. steer_jamming5=kron(t_a,s_a6);
    
81. 
    
82. C=[steer_signal];      % 构造线性约束  C'W=G
    
83. g=[1]; % 两个目标的无失真约束
    
84. 
    
85. amp_signal=sqrt(2*10^(snr1/10))*cos(2*pi*f0*t+randn);     % 信号加上个随机相位
    
86. 
    
87. %白噪声干扰，同时也属于宽带干扰
    
88. amp_jamming1=sqrt(10^(inr1/10))*randn(1,K);            
    
89. amp_jamming2=sqrt(10^(inr2/10))*randn(1,K);
    
90. amp_jamming3=sqrt(10^(inr3/10))*randn(1,K);
    
91. amp_jamming4=sqrt(10^(inr4/10))*randn(1,K);           
    
92. amp_jamming5=sqrt(10^(inr5/10))*randn(1,K);
    
93. %amp_jamming1=sqrt(10^(inr1/10))*wgn(1,K,0,'complex');
    
94. %amp_jamming2=sqrt(10^(inr2/10))*wgn(1,K,0,'complex');
    
95. %amp_jamming3=sqrt(10^(inr3/10))*wgn(1,K,0,'complex');
    
96. %amp_jamming4=sqrt(10^(inr4/10))*wgn(1,K,0,'complex');
    
97. %amp_jamming5=sqrt(10^(inr5/10))*wgn(1,K,0,'complex');
    
98. 
    
99. % 单音干扰
    
100. %amp_jamming1= sqrt(2*10^(inr1/10))*sin(2*pi*f0*t+randn);
     
101. %amp_jamming2= sqrt(2*10^(inr2/10))*sin(2*pi*f0*t+randn);
     
102. %amp_jamming3= sqrt(2*10^(inr3/10))*sin(2*pi*f0*t+randn);
     
103. %amp_jamming4= sqrt(2*10^(inr4/10))*sin(2*pi*f0*t+randn);
     
104. %amp_jamming5= sqrt(2*10^(inr5/10))*sin(2*pi*f0*t+randn);
     
105. 
     
106. noise=(randn(N*M,K));
     
107. 
     
108. jam=steer_jamming1*amp_jamming1+steer_jamming2*amp_jamming2+steer_jamming3*amp_jamming3...
     
109.    +steer_jamming4*amp_jamming4+steer_jamming5*amp_jamming5;
     
110. %jam=steer_jamming2*amp_jamming2;
     
111. SS=steer_signal*amp_signal;
     
112. X=SS+jam+noise;                             % 接收信号（信号+干扰+噪声）
     
113. 
     
114. OO=jam+noise;
     
115. Rs = 1/K*SS*SS';
     
116. Rni=1/K*OO*OO';
     
117. Rn=1/K*noise*noise';
     
118. Ri=1/K*jam*jam';
     
119. Rr=Rn+Ri;
     
120. Rx = 1/K*X*X';
     
121. invR=inv(Rx);
     
122. 
     
123. W_LCMV=g'*inv(C'*invR*C)*C'*invR;
     
124. 
     
125. theta=(0:1:89)*pi/180;         
     
126. phi =(0:1:359)*pi/180;
     
127. beam = zeros(360, 90);
     
128. for ii=1:length(phi)
     
129.     tempphi=phi(ii);
     
130.     waven_sweep=[(sin(theta).*cos(tempphi))', (sin(theta)*sin(tempphi))'];
     
131.     svsweep=exp(j*2*pi/wl*(waven_sweep*P)');
     
132.     steer = kron(t_a, svsweep);
     
133.     beam(ii,:)=W_LCMV*steer;
     
134. end
     
135.    
     
136. F=20*log10(abs(beam)/max((max(abs(beam)))));    % beampattern
     
137. 
     
138. % ref_point = 200;  % 该坐标位置值越大越好
     
139. % fitness1 = fun_F(F, 70, ref_point);
     
140. % fitness2 = fun_F(F, 60, ref_point);
     
141. % fitness3 = fun_F(F, 50, ref_point);
     
142. % fitness4 = fun_F(F, 40, ref_point);
     
143. % fitness5 = fun_F(F, 30, ref_point);
     
144. % fitness = fitness1 + fitness2 + fitness3 + fitness4 + fitness5;
     
145. 
     
146. ref_point = phi_signal;  % 该坐标位置值越大越好
     
147. fitness = fun_F(F, theta_signal, ref_point);
     
148. 
     
149. clear invR jam C c beam amp_jamming1 amp_jamming2 amp_jamming3 amp_jamming4 amp_jamming5 SS
     
150. clear steer_jamming1 steer_jamming2 steer_jamming3 steer_jamming4 steer_jamming5

复制代码

funF极值求解函数：

游客，如果您要查看本帖隐藏内容请回复

参考文献：【1】Symmetric Linear Antenna Array Geometry Synthe
【2】THINNING OF PLANAR CIRCULAR ARRAY ANTENNAS USING FIREFLY ALGORITHM

wal1985

太牛了，会这么多算法

mzch

非常好的代码，谢谢

huwenping

基于布谷鸟算法的导航波束形成问题

德帝

太棒了，我正找这个

fangh1993

学习了，谢谢楼主分享~

短离别

基于布谷鸟算法的导航波束形成问题，万分感谢

dlutww

学习了，谢谢楼主分享~

zhangsen1264

感谢 研究研究 ！！

木子900

谢谢分享，学习一下！