分布估计算法(Estimation of Distribution Algorithm,EDA)介绍及MATLAB代码(求解背包问题KP)

EDA算法

EDA算法原理

      通过一个概率模型描述候选解在空间得分布,采用统计学习手段从群体宏观的角度建立一个描述解分布得概率模型,然后对概率模型随机采样产生新的种群,如此反复进行,实现种群得进化,直到终止条件。(—建模-采用-建模-采样----循环)

EDA的不同变体

      EDA有很多不同的变体,包括如下一些,想详细了解可以参照参考文献。
      变量无关:PBIL、UMDA、cGA算法
      双变量相关:MIMIC、BMDA算法
      多变量相关:ECGA、FDA、BOA算法

EDA算法流程

通用算法步骤为:
      Step1 随机生成M个个体作为初始种群 ;
      Step2 对第L代种群计算个体适应度,判断是否满足终止条件,若满足,终止循环;若不满足,继续进行。
      Step3 根据适应度数值从种群中选出前N个(N≤M)优势个体,组成第L+1代的优势子种群 ;
      Step4 根据优势子种群更新概率模型 ;
      Step5 对概率模型进行随机采样,生成新种群(规模M),返回Step2.

具体介绍两个变体UMDA(Univariate marginal distribution algorithm)和PBIL(Population based incremental learning),主要区别在于更新概率模型的公式不同

EDA算法变体UMDA

      由德国学者Muhlenbein在1996年提出,算法描述:
      Step1 随机产生M个个体作为初始种群;
      Step2 然后计算M个个体的适应值,如果符合终止条件,算法结束,否则继续进行;
      Step3 选择最优的N个个体用来更新概率向量p(x), N <= M
更新过程:
p_l (x)=p(x│D_l^S )=1/N ∑_(k=1)N▒x_lk
      Step4 由新的概率模型采样M次,得到新一代群体,返回Step2

EDA算法变体PBIL

     由美国卡耐基梅隆大学的Baluja在1994年提出,算法描述:
     Step1 随机产生M个个体作为初始种群;
     Step2 然后计算M个个体的适应值,如果符合终止条件,算法结束,否则继续进行;
     Step3 选择最优的N个个体用来更新概率向量p(x), N <= M
更新过程:
在这里插入图片描述
     Step4 由新的概率模型采样M次,得到新一代群体,返回Step2

测试算例:KP

     本测试算例为背包问题,来源自https://people.sc.fsu.edu/~jburkardt/datasets/knapsack_01/knapsack_01.html
算例如下表:
在这里插入图片描述

MATLAB代码

UMDA代码

function [bestFit]=binaryUMDA(weightMax,weight_Individual,value_Individual)
%EDA_UMDA Algorithm,二进制编程
%参数
iterations = 500;%迭代最大次数
populationSize = 200;%种群规模
dimensionality = size(weight_Individual,2);%维度
dominantNum = 20;%优势群体个数
% weightMax             背包最大容量/重量
% weight_Individual     每个变量的体积/重量
% value_Individual      每个变量的价值
  
%初始化种群
probability = 0.5*ones(1,dimensionality);%初始化概率0.5
Best_Individual=zeros(iterations,dimensionality+1);%每次迭代中的最优解,(:,1)存储适应度值

%循环迭代
for I=1:iterations
    %按照概率模型创建样本
    flag = 0;i=1;
    while i<=populationSize
        r=rand(1,dimensionality);
        Species(i,:)=1.*(r<probability);
        %判断是否超出容量范围
        weightSum = sum(Species(i,:).*weight_Individual,2);
        if flag>=20
            Species(i,:) = zeros(1,dimensionality);%多次仍未符合要求,舍弃
            flag = 0;
        elseif weightSum>weightMax
            i=i-1;flag=flag+1;
        else
            flag = 0;
        end
        i=i+1;
    end
    %计算适应度
    Fitness_Value=zeros(populationSize,1);%创建适应度计算
    for i=1:populationSize
         Fitness_Value(i,1)=sum(Species(i,:).*value_Individual,2);
    end
    %排序
    [Fitness,index] = sort(Fitness_Value);%=》大排序
    %每次迭代中的最优解
    Best_Individual(I,1) = I;
    Best_Individual(I,2) = Fitness_Value(index(populationSize));
    for i=3:dimensionality+2
        Best_Individual(I,i) = Species(index(populationSize),i-2);
    end
    %选取优势种群
    dominantSpecies=zeros(dominantNum,dimensionality);%创建选取的优势群体<populationSize,(:,1)存储适应度值
    for i=1:dominantNum
        dominantSpecies(i,:) = Species(index(populationSize-dominantNum+i),:);
    end
    %更新概率模型 划重点
    Ones_Number = sum(dominantSpecies);
    probability = Ones_Number/dominantNum;
end

%画图函数
bestFit =Best_Individual(iterations,[2:dimensionality+2]);
% 输出,画出fitness-iterations
plot(Best_Individual(:,1),Best_Individual(:,2));
disp(strcat('最优结果',':',num2str(Best_Individual(iterations,2))));

PBIL代码

function [bestFit]=binaryPBIL(weightMax,weight_Individual,value_Individual)
%EDA_PBIL Algorithm,二进制编程
%参数
iterations = 500;%迭代最大次数
populationSize = 200;%种群规模
dimensionality = size(weight_Individual,2);%维度
learningRate =0.3;% ((0.001*(dimensionality)^2)>1)*1+((0.001*(dimensionality)^2)<=1)*(0.001*(dimensionality)^2);
%学习效率太小会不收敛
dominantNum = 20;%优势群体个数
% weightMax             背包最大容量/重量
% weight_Individual     每个变量的体积/重量
% value_Individual      每个变量的价值
  
%初始化种群
probability = 0.5*ones(1,dimensionality);%初始化概率0.5
Best_Individual=zeros(iterations,dimensionality+1);%每次迭代中的最优解,(:,1)存储适应度值

%循环迭代
for I=1:iterations
    %按照概率模型创建样本
    flag = 0;i=1;
    while i<=populationSize
        r=rand(1,dimensionality);
        Species(i,:)=1.*(r<probability);
        %判断是否超出容量范围
        weightSum = sum(Species(i,:).*weight_Individual,2);
        if flag>=20
            Species(i,:) = zeros(1,dimensionality);%多次仍未符合要求,舍弃
            flag = 0;
        elseif weightSum>weightMax
            i=i-1;flag=flag+1;
        else
            flag = 0;
        end
        i=i+1;
    end
    %计算适应度
    Fitness_Value=zeros(populationSize,1);%创建适应度计算
    for i=1:populationSize
         Fitness_Value(i,1)=sum(Species(i,:).*value_Individual,2);
    end
    %排序
    [Fitness,index] = sort(Fitness_Value);%=》大
    %每次迭代中的最优解
    Best_Individual(I,1) = I;
    Best_Individual(I,2) = Fitness_Value(index(populationSize));
    for i=3:dimensionality+2
        Best_Individual(I,i) = Species(index(populationSize),i-2);
    end
    %选取优势种群
    dominantSpecies=zeros(dominantNum,dimensionality);%创建选取的优势群体<populationSize,(:,1)存储适应度值
    for i=1:dominantNum
        dominantSpecies(i,:) = Species(index(populationSize-dominantNum+i),:);
    end
    %更新概率模型  划重点
    Ones_Number = sum(dominantSpecies);
    probability = (1-learningRate)*probability+learningRate*Ones_Number/dominantNum;    
end

%画图函数
bestFit =Best_Individual(iterations,[2:dimensionality+2]);
% 输出,画出fitness-iterations
% plot(Best_Individual(:,1),Best_Individual(:,2));
disp(strcat('最优结果',':',num2str(Best_Individual(iterations,2))));

实验结果展示

UMDA结果

运行结果在如下参数的前提下得出
种群规模:200
优势群体个数:20
种群计算迭代次数:500
备注:算法运行结果为求50次运算结果的平均值,误差=(理论最优解-算法运行结果)/理论最优解*100%

在这里插入图片描述

PBIL结果

运行结果在如下参数的前提下得出
学习效率:0.3
种群规模:200
优势群体个数:20
种群计算迭代次数:500
备注:算法运行结果为求50次运算结果的平均值,误差=(理论最优解-算法运行结果)/理论最优解*100%

在这里插入图片描述

参考文献

[1] 周树德, 孙增圻. 分布估计算法综述[J]. 自动化学报, 2007, 33(2).

  • 11
    点赞
  • 93
    收藏
    觉得还不错? 一键收藏
  • 2
    评论
APES算法是一种高分辨率的DOA估计方法,通过对信号进行空间谱估计,可以得到信号的幅度和相位信息,从而实现目标的定位。以下是APES算法matlab代码实现: ``` %APES算法实现 clc; clear all; close all; %% 参数设置 N = 100; % 信号长度 M = 5; % 阵元数 d = 0.5; % 阵元间距 theta = [30 60]; % 信号方向 K = length(theta); % 信号数 SNR = 10; % 信噪比 fs = 1000; % 采样频率 f = [100 200]; % 信号频率 w = 2 * pi * f / fs; % 角频率 %% 生成信号和加噪声 s = zeros(K, N); for k = 1:K s(k,:) = exp(1i * w(k) * (0:N-1)) + SNR * randn(1,N) + SNR * 1i * randn(1,N); % 生成信号 end %% 矩阵构造 A = zeros(M,K); for m = 1:M for k = 1:K A(m,k) = exp(-1i * 2 * pi * d * (m-1) * sin(theta(k) * pi / 180)); % 构造矩阵 end end %% 噪声子空间信号幅度矩阵和相位矩阵 Rn = s * s' / N; % 噪声协方差矩阵 [V,D] = eig(Rn); % 噪声协方差矩阵特征分解 [~,index] = sort(diag(D),'ascend'); % 特征值升序排序 En = V(:,index(1:M-K)); % 噪声子空间 B = En * En'; % 噪声子空间信号幅度矩阵 P = angle(B); % 噪声子空间信号相位矩阵 %% 信号子空间谱估计 R = A * B * A'; % 信号协方差矩阵 [V,D] = eig(R); % 信号协方差矩阵特征分解 [~,index] = sort(diag(D),'descend'); % 特征值降序排序 Un = V(:,index(1:K)); % 信号子空间 %% DOA估计 theta_range = -90:0.1:90; % 方向搜索范围 Pmusic = zeros(1,length(theta_range)); for i = 1:length(theta_range) a = zeros(M,1); for m = 1:M a(m) = exp(-1i * 2 * pi * d * (m-1) * sin(theta_range(i) * pi / 180)); % 构造矩阵 end Pmusic(i) = 1 / (a' * inv(Rn) * a) / (a' * inv(R) * a); % 谱估计 end %% 绘图 figure; plot(theta_range,Pmusic); xlabel('DOA/°'); ylabel('谱估计'); title('APES算法DOA估计'); ``` 需要注意的是,APES算法的实现需要注意参数的设置和矩阵构造。在此代码中,我们使用了阵元数为5,阵元间距为0.5,信号方向为30度和60度,信噪比为10dB,采样频率为1000Hz,信号频率为100Hz和200Hz的信号进行了实验。在实现中,我们首先构造了信号矩阵和噪声子空间信号幅度矩阵和相位矩阵,然后通过信号子空间谱估计来实现DOA估计。最后,我们绘制了DOA估计的谱估计图像。通过调整参数,可以实现不同场景下的DOA估计

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 2
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值