基于免疫算法的认知无线电资源分配优化算法的matlab仿真

1.算法描述

   认知无线电（CR）的概念来自Joseph Mitolo博士1999年的开创性工作。它自适应地调整内部通信机制，通过学习，了解等实时变化特定的无线电操作参数（功率，载波调制和编码等），适应外部无线电环境，并独立空闲频谱它可以被搜索和使用。这有助于用户选择用于无线传输的最佳和最合适的服务，甚至允许基于现有或未来的无线电资源来延迟甚至主动地开始传输。

免疫算法基本步骤：

1、随机产生初始父代种群A1，根据先验知识抽取疫苗；

2、若当前群体中包含最佳个体，则算法停止运行并输出结果；否则，继续；

3、对当前第k代父本种群Ak进行交叉操作，得到种群Bk；

4、对Bk进行变异操作，得到种群Ck；

5、对Ck进行接种疫苗操作，得到种群Dk；

6、对Dk进行免疫选择操作，得到新一代父本Ak+1，转至第二步。

免疫算法进行子载波分配流程
免疫算法的主要操作如下：

抗体的编码和产生每个抗体由M位十进制整数组成，对应N个子载波，每个整数范围从1到M，表明子载体被分配给M个认知用户中的一个。抗体群是K，对应于K个解决方案。第一代抗体通常是随机产生的。
计算和排序亲和力 根据贪婪算法的基本思想，完成每种抗体的安全传输速率调整。安全传输速率调整完成后，获得每个抗体对应的加权安全吞吐量，然后排序
（3）存储单元的生成在（2）中排序后对应于较高加权吞吐量的抗体被保存为存储单元以避免后续操作中的丢失。

（4）蚂蚁的选择ibodies

浓度概率，亲和概率和选择概率有三个概率。

2.仿真效果预览
matlab2022a仿真结果如下：

3.MATLAB核心程序

    
    pl_cu_pu(1,m)=((loc_cu(m,1)-loc_pu(1,1))^2+(loc_cu(m,2)-loc_pu(1,2))^2)^0.5;
end
pl_cu_pu=10^(0.1*randn*Sigma)./((pl_cu_pu/3).^Eta);  %产生主用户和认知用户的路径损耗
pl_cu_cbs=10^(0.1*randn*Sigma)./((pl_cu_cbs/3).^Eta);  %产生用户和认知基站的路径损耗
for m=1:M
    for i=1:M
        pl_cu_cu(m,i)=((loc_cu(m,1)-loc_cu(i,1))^2+(loc_cu(m,2)-loc_cu(i,2))^2)^0.5;
    end
end
pl_cu_cu=10^(0.1*randn*Sigma)./((pl_cu_cu/3).^Eta);   %产生认知用户和认知用户之间的路径损耗
 
 
%对主用户干扰
to_pu_g=raylrnd(1,[M,N]);   %对主用户的信道干扰增益（第ｍ行是第ｍ个小ｆ中各个子载波的增益）
lamda=zeros(1,N);  %对主用户距离干扰增益
Distance=zeros(1,N);   %第ｎ个子载波和主用户的频带距离
for n=1:N
    Distance(1,n)=((N-(2*n-1)/2)*Delta_f+W/2);
end
F=@(f)((sin(pi*Ts*f)./(pi*Ts*f)).^2);      %积分函数
for n=1:N      %Ｌａｍｄａ计算
    lamda(1,n)=Ts*quad(F,(Distance(1,n)-W/2),(Distance(1,n)+W/2));
end
 
for m=1:M
    for n=1:N
        to_pu_g(m,n)=to_pu_g(m,n)*lamda(1,n);   %对主用户的总的干扰增益
    end
end
 
for m=1:M
    to_pu_g(m,:)=to_pu_g(m,:)*pl_cu_pu(1,m);
end
 
 
%同频相互干扰
co_ch_g=zeros(M,M,N);  %同频干扰矩阵，第ｎ页是第ｎ个子载波的同频干扰
 
for n=1:N
    for i=1:M
        for j=i:M              %若要不同的同频干扰，则这里改为１：Ｍ
            co_ch_g(i,j,n)=raylrnd(1);
            co_ch_g(j,i,n)=co_ch_g(i,j,n);
        end
    end
end
for n=1:N
    for m=1:M
        co_ch_g(m,m,n)=0;
    end
end
for n=1:N
    co_ch_g(:,:,n)=(co_ch_g(:,:,n).^2).*pl_cu_cu;
end
for n=1:N
    for i=1:M
        co_ch_g(i,i,n)=0;
    end
end
 
%—有关各个子载波的信道增益——
ch_g=raylrnd(1,[M,N]);   %每个子载波的信道增益
for m=1:M
    ch_g(m,:)=(ch_g(m,:).^2)*pl_cu_cbs(1,m);
end
%＝＝＝—有关运算结果显示＝＝＝＝＝二＝
rate_gen1=zeros(1,genMax);   %每代的最优速率
rate_gen2=zeros(1,genMax);      %每代的最优速率
%rate_gen3=zeros(1,generaMax);   %每代的最优小速率
ratio_gen1=zeros(1,genMax);     %每代的最优速率功率比
ratio_gen2=zeros(1,genMax);     %每代的最优速率功率比
%ratio_gen3=zeros(1,genMax);     %每代的最优速率功率比
Pth=10;
for Ith=[0.001 0.005 0.01]
    for qq=1:10
        %＝＝＝＝＝＝＝＝＝程序初始化＝＝＝＝＝＝＝＝
        temp_antibody(:)=Pth/N;       %初始抗体功率是平均分配
        temp_memorycell(:)=Pth/N;      %初始记忆细胞功率也是平均分配
        antibody(:)=temp_antibody(:);      %抗体变量装入初始抗体
        memorycell(:)=temp_memorycell(:);      %记忆细胞变量装入记忆细胞
        
        
        %＝＝＝＝＝＝＝＝＝主函数＝＝＝＝＝＝＝＝＝
        for gen=1:genMax
            gen
            
            %记忆细胞添加            
            antibody(:,:,(K-L+1):K)=memorycell;
            %antibody_falg((K-L+1):K)=1;                       %记录成型抗体
            
            
            
            %——抗体合格检测
            exceed_Ith=0;   %超过干扰门限的量
            exceed_Pth=0;      %超过功率门限的量
            decr_var_I=zeros(M,N);      %有关干扰的归一化的减少量
            decr_val_P=zeros(M,N);     %有关功率的归一化的减少量
            ratio_P=zeros(1,3);      %计算当前功率与公平性的比值
            standerd_site=0;      %作为基准的用户
            
            for k=1:K
                decr_val_P=antibody(:,:,k)/sum(sum(antibody(:,:,k)));   %根据现有功率计算归一化功率减少量
                if(sum(sum(antibody(:,:,k)))-Pth>eps)
                    exceed_Pth=sum(sum(antibody(:,:,k)))-Pth;   %计算超出多少
                    antibody(:,:,k)=antibody(:,:,k)-exceed_Pth*decr_val_P(:,:);
                end
                
                itfere_pu(:,:,k)=to_pu_g.*antibody(:,:,k);
                decr_val_I=itfere_pu(:,:,k)/sum(sum(itfere_pu(:,:,k)));      %根据现有干扰计算归一化功率减少量
                if(sum(sum(itfere_pu(:,:,k)))-Ith>eps)
                    exceed_Ith=sum(sum(itfere_pu(:,:,k)))-Ith;         %计算超出多少
                    
                    antibody(:,:,k)=antibody(:,:,k)-(exceed_Ith*decr_val_I)./(to_pu_g);
                end
                itfere_pu(:,:,k)=to_pu_g.*antibody(:,:,k);      %更新对主用户的干扰
            end
            
            
            
            
            
            
            %---速率计算
            h=zeros(M,N,K);   %各个子载波的信噪比
            i_ss=zeros(M,N,K);   %各个子载波受到的同频干扰
            for k=1:K
                for n=1:N
                    for m=1:M
                        for i=1:M
                            i_ss(m,n,k)=i_ss(m,n,k)+co_ch_g(i,m,n)*antibody(i,n,k);   %计算各个子载波所受到的同频干扰
                        end
                    end
                end
            end
            
            temp=sum(i_ss);
            
            for k=1:K
                for n=1:N
                    for m=1:M
                        i_ss(m,n,k)=temp(1,n,k)-i_ss(m,n,k);      %计算每个子载波所受到的同频干扰
                    end
                end
            end
            for k=1:K
                h(:,:,k)=ch_g./(N0*Delta_f+i_ss(:,:,k));           %计算信噪比
            end
            
            %rate=Delta_f*log2(1+antibody.*h/Tau);          %计算每个抗体每个小区的速率
            rate=log2(1+antibody.*h/Tau);         %计算每个抗体每个小区的速率
            
            
            %总的公平性调整
            for o=1:O
                %rate_m=zeros(1,3,K);      %各个小区的总速率
                %ratio_r=zeros(1,3);             %计算当前速率与公平性的比值
                standard_site=0;      %作为基准的用户
                decr_val_r=zeros(M,N);        %有关功率的归一化的减少量
                exceed_r=0;            %超过干扰门限的量
                rate_m=sum(rate,2);
                for k=1:K
                    ratio_r=rate_m(:,:,k)./Rate_m';
                    [~,standard_site]=min(ratio_r);     %选取基准功率
                    switch(standard_site)      %进行功率公平性调整
                        case 1   %第一个为基准
                            for m=[2 3]
                                exceed_r=rate_m(1,:,k)/2-rate_m(m,:,k);
                                decr_val_r(m,:)=(rate(m,:,k)/sum(rate(m,:,k)))*exceed_r;
                                if(exceed_r<0)
                                    for n=1:N
                                        antibody(m,n,k)=(2^decr_val_r(m,n)*(Tau+antibody(m,n,k)*h(m,n,k))-Tau)/h(m,n,k);
                                    end
                                end
                            end
                            
                        case 2            %第二个为基准
                            exceed_r=rate_m(2,:,k)*2-rate_m(1,:,k);
                            decr_val_r(1,:)=(rate(1,:,k)/sum(rate(1,:,k)))*exceed_r;
                            if(exceed_r<0)
                                for n=1:N
                                    antibody(1,n,k)=(2^decr_val_r(1,n)*(Tau+antibody(1,n,k)*h(1,n,k))-Tau)/h(1,n,k);
                                end
                            end
                            exceed_r=rate_m(2,:,k)-rate_m(3,:,k);
                            decr_cal_r(3,:)=(rate(3,:,k)/sum(rate(3,:,k)))*exceed_r;
                            if(exceed_r<0)
                                for n=1:N
                                    antibody(3,n,k)=(2^decr_cal_r(3,n)*(Tau+antibody(3,n,k)*h(3,n,k))-Tau)/h(3,n,k);
                                end
                            end
                            
                             
                        case 3              %第三个为基准
                            exceed_r=rate_m(3,:,k)*2-rate_m(1,:,k);
                            decr_val_r(1,:)=(rate(1,:,k)/sum(rate(1,:,k)))*exceed_r;
                            if(exceed_r<0)
                                for n=1:N
                                    antibody(1,n,k)=(2^decr_val_r(1,n)*(Tau+antibody(1,n,k)*h(1,n,k))-Tau)/h(1,n,k);
                                end
                            end
                            exceed_r=rate_m(3,:,k)-rate_m(2,:,k);
                            decr_val_r(2,:)=(rate(2,:,k)/sum(rate(2,:,k)))*exceed_r;
                            if(exceed_r<0)
                                for n=1:N
                                    antibody(2,n,k)=(2^decr_val_r(2,n)*(Tau+antibody(2,n,k)*h(2,n,k))-Tau)/h(2,n,k);
                                end
                            end
                    end
                end
                temp=antibody<=eps;
                antibody(temp)=eps;
                h=zeros(M,N,K);      %各个子载波的信噪比
                i_ss=zeros(M,N,K);          %各个子载波受到的同频干扰
                for k=1:K
                    for n=1:N
                        for m=1:M
                            for i=1:M
                                i_ss(m,n,k)=i_ss(m,n,k)+co_ch_g(i,m,n)*antibody(i,n,k);      %计算各个子载波所受到的同频干扰
                            end
                        end
                    end
                end
                temp=sum(i_ss);
               
                for k=1:K
                    for n=1:N
                        for m=1:M
                            i_ss(m,n,k)=temp(1,n,k)-i_ss(m,n,k);      %计算每个子载波所受到的同频干扰
                        end
                    end
                end
                
                for k=1:K
                    h(:,:,k)=ch_g./(N0*Delta_f+i_ss(:,:,k));             %计算信噪比
                end
                %rate=Delta_f*log2(1+antibody.*h/Tau);   %计算每个抗体每个小区的速率
                rate=log2(1+antibody.*h/Tau);   %计算每个抗体每个小区的速率
            end
            
            %速率功率比计算
            %if gen==genMax
            %rate=floor(rate);
            %end
            rate_antibody=sum(sum(rate));               %计算每个抗体速率
            ratio=rate_antibody./sum(sum(antibody));           %计算每个抗体的速率功率比
            Objective=ratio*Alpha+rate_antibody*(1-Alpha);               %计算优化目标的值
            
            
            
            %亲和力排序
            for i=1:K-1
                for j=i+1:K
                    if(Objective(i)<Objective(j))
                        temp=rate_antibody(i);         %每代速度交换
                        rate_antibody(i)=rate_antibody(j);
                        rate_antibody(j)=temp;
                        
                        temp=rate(:,:,i);               %速率交换
                        rate(:,:,i)=rate(:,:,j);
                        rate(:,:,j)=temp;
                        
                        
                        temp=antibody(:,:,i);            %功率交换
                        antibody(:,:,i)=antibody(:,:,j);
                        antibody(:,:,j)=temp;
                        
                        
                        temp=itfere_pu(:,:,i);            %对主用户干扰交换
                        itfere_pu(:,:,i)=itfere_pu(:,:,j);
                        itfere_pu(:,:,j)=temp;
                        
                        
                        temp=Objective(i);           %每代速度交换
                        Objective(i)=Objective(j);
                        Objective(j)=temp;
                        
                        
                        temp=ratio(i);                 %每代速度交换
                        ratio(i)=ratio(j);
                        ratio(j)=temp;
                        
%                         temp=antibody_flag(i);
%                         antibody_flag(i)=antibody_flag(j);
%                         antibody_flag(j)=temp;
                        
                    end
                end
            end
            
            
            
            %提取记忆细胞
            memorycell=antibody(:,:,1:L);                %前Ｌ个个体作为记忆细胞提取
            
            
            %算法结果等待
            rate_gen1(gen)=rate_antibody(1);   %记录每代最优的最小速度
            ratio_gen1(gen)=Objective(1);      %每代最优速率功率比
            
            %抗体选择
            %pd=zeros(1,K);             %浓度概率
            pf=zeros(1,K);                 %适应度概率
            %p=zeros(1,K);                %选择概率
            resemble=zeros(1,K);         %抗体的相似度
            res=zeros(M,N,K);              %各个子载波的距离
            index=zeros(1,K);               %轮盘赌抗体索引
            
            for k=1:K
                for i=1:K
                    res(:,:,i)=abs(antibody(:,:,k)-antibody(:,:,i));               %计算各个子载波的距离
                end
                temp=(res<=Gama*Pth);             %距离过近的子载波个数
                temp=sum(sum(temp));                %计算与每个抗体的距离
                for i=1:K
                    resemble(k)=resemble(k)+temp(i);
                end
            end
            pd=1-resemble/(M*N*K);
            
            for k=1:K
                pf(k)=(Objective(k)-Objective(K)+eps)/(Objective(1)-Objective(K)+eps);
            end
            p=Rho*pf+(1-Rho)*pd;
            
            for i=1:K-1                   %抗体按照选择概率排序
                for j=i+1:K
                    if(p(i)<p(j))
                        temp=p(i);
                        p(i)=p(j);
                        p(j)=temp;
                        temp=antibody(i);
                        antibody(i)=antibody(j);
                        antibody(j)=temp;
                    end
                end
            end
            
            
            for k=1:K                        %选出K个抗体
                temp=rand*sum(p);
                for i=1:K
                    temp=temp-p(i);
                    if(temp<=0)
                        index(k)=i;             %赌出第ｋ个索引值
                        break;
                    end
                end
            end
            temp=antibody;
            antibody=temp(:,:,index);
            
            
            %两种交叉
            [antibody]=dou_atb_excha(N,K,Pc,antibody);               %调用交叉互换函数ｉ现ｉ叉互换
            %两种突变
            [antibody]=dou_atb_mutat(K,N,M,Pm,antibody);            %调用突变函数实现抗体突变
            
        end
        obj(qq,:)=ratio_gen1;
        obj11(qq,:)=2*rate_gen1;
    end
    hold on;
    oobj=mean(obj);
    oobj11=mean(obj11);
    plot(oobj11);
    xlabel('代数');ylabel('亲和力');title('收敛性能');
    legend('Ith=0.001','Ith=0.005','Ith=0.01','Location','NorthEast');
    
 
end