头脑风暴算法（BSO），原理公式详解，附matlab代码

头脑风暴算法（Brain Storming Optimization Algorithm, BSO）是一种新型群智能优化算法，该算法灵感来源于头脑风暴法，当一群人围绕一个特定的兴趣领域产生新观点的时候，这种情境就叫做头脑风暴。头脑风暴法是一种充分开发人类创造性思维解决问题的方法，具有进化能力强、搜索速度快、寻优能力强的特点。该成果于2011年发表在知名SCI期刊Advances in Swarm Intelligence上，目前在谷歌学术上被引773次。

在一个头脑风暴的过程中，通常会涉及到三种人，即主持人、问题负责人和一组头脑风暴人员。主持人主要是为了引导头脑风暴过程顺利进行而不带偏见，这意味着一个有主持经验但不了解问题的主持人将是首选；头脑风暴小组是指在主持人的引导下，通过头脑风暴过程实际产生想法的一群人；问题所有者负责在每一轮创意生成中挑选他们认为更好的创意。一般来说，每一个头脑风暴的过程都要经过三轮的想法产生过程，因为经过三轮的想法产生过程，人就会筋疲力尽，效率也就不高了。

1、算法原理

（1）头脑风暴阶段

在BSO算法中，一个人代表一个想法，而这个想法又代表了待解决问题的潜在解决方案。类似于其他群体智能算法，均匀随机初始化是首选，特别是当没有先验知识是已知的问题要解决的。在原始的BSO中，通过根据以下公式将高斯随机值添加到选定的个体来生成新的个体：

其中，xi new表示新生成的个体的第i个维度，可以看作是产生的新想法；xi selected表示基于当前个体群体的所选个体的第i个维度，可以看作是原来的想法；random()是一个随机函数，它是原BSO中的高斯随机函数。

在原始BSO中，系数ξ(t)是加权高斯随机值对新生成值的贡献的步长，ξ(t)值由下式求得

其中，logsig()为对数sigmoid传递函数，T为最大迭代次数，t为当前迭代次数，k是控制logsig( )斜率的系数，rand()返回(0,1)范围内的随机值。

其中，原来的想法xi selected可以根据以下公式获得

其中，xi selected为xold1和xold2的加权和，rand为随机产生的0到1之间的系数，xold1和xold2是分别从两个集群中选取的想法，通过这个随机值实现两个个体想法的交互和融合。

（2）算法伪代码：

---

1.初始化种群规模I和维度dim、最大迭代次数1

2.初始化参数和种群 Xi(i = 1,2,..., )

3. While t < T do

4. 利用 k-means 聚类算法将 n 个个体聚类为 m个聚类

5. 计算每个个体的适应度值并排名，将最佳个体记录为每个聚类中的聚类中心

6. If rand < 0.2 then

7. 随机生成一个个体来替换随机选取的一个聚类的中心

8. End If

9. For i = 1 to N do

10. If rand < 0.8 then

11. 选取一个聚类的中心或随机个体

12. Else

13. 选取两个聚类的中心或随机个体

14. End lf

15. 通过式(1)产生新的个体，计算适应度值

16. 将新生成的个体与现有个体进行比较，将保留较好的个体

17. End For

18. End While

19. 输出最终解

---

 MATLAB核心代码

function best_fitness = bso(fun,n_p,n_d,n_c,rang_l,rang_r,max_iteration)
% fun = fitness_function
% n_p; population size
% n_d; number of dimension
% n_c: number of clusters
% rang_l; left boundary of the dynamic range
% rang_r; right boundary of the dynamic range


prob_one_cluster = 0.8; % probability for select one cluster to form new individual; 


stepSize = ones(1,n_d); % effecting the step size of generating new individuals by adding random values


popu = rang_l + (rang_r - rang_l) * rand(n_p,n_d); % initialize the population of individuals


popu_sorted  = rang_l + (rang_r - rang_l) * rand(n_p,n_d); % initialize the  population of individuals sorted according to clusters


n_iteration = 0; % current iteration number


% initialize cluster probability to be zeros
prob = zeros(n_c,1);


best = zeros(n_c,1);  % index of best individual in each cluster


centers = rang_l + (rang_r - rang_l) * rand(n_c,n_d);  % initialize best individual in each cluster
centers_copy = rang_l + (rang_r - rang_l) * rand(n_c,n_d);  % initialize best individual-COPY in each cluster FOR the purpose of introduce random best


best_fitness = 1000000*ones(max_iteration,1);
fitness_popu = 1000000*ones(n_p,1);  % store fitness value for each individual
fitness_popu_sorted = 1000000*ones(n_p,1);  % store  fitness value for each sorted individual


indi_temp = zeros(1,n_d);  % store temperary individual


% calculate fitness for each individual in the initialized population
for idx = 1:n_p
    fitness_popu(idx,1) = fun(popu(idx,:));
end


while n_iteration < max_iteration
    
    cluster = kmeans(popu, n_c,'Distance','cityblock','Start',centers,'EmptyAction','singleton'); % k-mean cluster


    % clustering    
    fit_values = inf*ones(n_c,1);  % assign a initial big fitness value  as best fitness for each cluster in minimization problems


    number_in_cluster = zeros(n_c,1);  % initialize 0 individual in each cluster
           
    for idx = 1:n_p
        number_in_cluster(cluster(idx,1),1)= number_in_cluster(cluster(idx,1),1) + 1;
               
        % find the best individual in each cluster
        if fit_values(cluster(idx,1),1) > fitness_popu(idx,1)  % minimization
            fit_values(cluster(idx,1),1) = fitness_popu(idx,1);
            best(cluster(idx,1),1) = idx;
        end
            
    end
    
    % form population sorted according to clusters
    counter_cluster = zeros(n_c,1);  % initialize cluster counter to be 0 
    
    acculate_num_cluster = zeros(n_c,1);  % initialize accumulated number of individuals in previous clusters
    
    for idx =2:n_c
        acculate_num_cluster(idx,1) = acculate_num_cluster((idx-1),1) + number_in_cluster((idx-1),1);
    end
    
    
    %start form sorted population
    for idx = 1:n_p
        counter_cluster(cluster(idx,1),1) = counter_cluster(cluster(idx,1),1) + 1 ;
        temIdx = acculate_num_cluster(cluster(idx,1),1) +  counter_cluster(cluster(idx,1),1);
        popu_sorted(temIdx,:) = popu(idx,:);
        fitness_popu_sorted(temIdx,1) = fitness_popu(idx,1);
    end
       
    % record the best individual in each cluster
    for idx = 1:n_c
        centers(idx,:) = popu(best(idx,1),:);        
    end
    
    centers_copy = centers;  % make a copy
    
    if (rand() < 0.2) %  select one cluster center to be replaced by a randomly generated center
        cenIdx = ceil(rand()*n_c);
        centers(cenIdx,:) = rang_l + (rang_r - rang_l) * rand(1,n_d);
    end 
           
    % calculate cluster probabilities based on number of individuals in
    % each cluster
    for idx = 1:n_c
        prob(idx,1) = number_in_cluster(idx,1)/n_p;
        if idx > 1
            prob(idx,1) = prob(idx,1) + prob(idx-1,1);
        end
    end
    
    % generate n_p new individuals by adding Gaussian random values
                   
    for idx = 1:n_p
        r_1 = rand();  % probability for select one cluster to form new individual
        if r_1 < prob_one_cluster % select one cluster
            r = rand();
            for idj = 1:n_c
                if r < prob(idj,1)                      
                    if rand() < 0.4  % use the center
                       indi_temp(1,:) = centers(idj,:); 
                    else % use one randomly selected  cluster
                        indi_1 = acculate_num_cluster(idj,1) + ceil(rand() * number_in_cluster(idj,1));
                        indi_temp(1,:) = popu_sorted(indi_1,:);  
                    end
                    break
                end
            end
        else % select two clusters
            % pick two clusters 
            cluster_1 = ceil(rand() * n_c);
            indi_1 = acculate_num_cluster(cluster_1,1) + ceil(rand() * number_in_cluster(cluster_1,1));
            
            cluster_2 = ceil(rand() * n_c);
            indi_2 = acculate_num_cluster(cluster_2,1) + ceil(rand() * number_in_cluster(cluster_2,1));
            
            tem = rand();
            if rand() < 0.5 %use center
                indi_temp(1,:) = tem * centers(cluster_1,:) + (1-tem) * centers(cluster_2,:); 
            else   % use randomly selected individuals from each cluster            
                indi_temp(1,:) = tem * popu_sorted(indi_1,:) + (1-tem) * popu_sorted(indi_2,:); 
            end
        end         
        
        stepSize = logsig(((0.5*max_iteration - n_iteration)/20)) * rand(1,n_d);
        indi_temp(1,:) = indi_temp(1,:) + stepSize .* normrnd(0,1,1,n_d);
        % if better than the previous one, replace it
        fv = fun(indi_temp);
        if fv < fitness_popu(idx,1)  % better than the previous one, replace
            fitness_popu(idx,1) = fv;
            popu(idx,:) = indi_temp(1,:);
        end 
        
    end
        
    % keep the best for each cluster
    for idx = 1:n_c
        popu(best(idx,1),:) = centers_copy(idx,:);  
        fitness_popu(best(idx,1),1) = fit_values(idx,1);
    end
       
    n_iteration = n_iteration +1;
    
    % record the best fitness in each iteration
    best_fitness(n_iteration, 1) = min(fit_values);
end

参考文献

[1] Shi Y. Brain storm optimization algorithm[C]//Advances in Swarm Intelligence: Second International Conference, ICSI 2011, Chongqing, China, June 12-15, 2011, Proceedings, Part I 2. Springer Berlin Heidelberg, 2011: 303-309.

完整代码获取方式：后台回复关键字：

TGDM100