数学建模（二） 启发式算法

启发式算法

在很长一段时间里，提起数学建模，我就能想到遗传算法，模拟退火，粒子群。提起遗传算法，我就能想到函数GA（）。我打开CSDN，输入这些关键字，就可以很轻松的看到这些算法的原理，过程清楚明白。然而，我始终不明白的是，这些算法跟实际问题有啥关系？
带着疑问，我试图寻找使用它们的例子，好的，我又一次次看到了上一篇中的优化问题，求最小，求最大。我的疑惑更深了，这些问题不是用一个函数就解出来了吗？或者实在不行，我也可以拽个GA解出来啊。首先是我不知道数学建模的题如何抽象成数学问题，也就是建模，其次我不知道数学问题如何用特定方法求解。直到我看了几篇论文的讲解视频，我才仿佛明白了一点，我又重新第二遍开始看这些算法的时候，我觉得我终于把这三个点联系到一起了（其实建模跟方法根本就是一个点吧）。
由于这些算法之间作用的相似性，这里就介绍两种，然而，真的要比赛的话，我觉得懂一种就可以了（并不值得借鉴）。

1.模拟退火

背景：统计力学表明材料中粒子的不同结构对应于粒子的不同能量水平。在高温条件下，粒子的能量较高，可以自由运动和重新排列。在低温条件下，粒子能量较低。如果从高温开始，非常缓慢地降温（这个过程被称为退火），粒子就可以在每个温度下达到热平衡。当系统完全被冷却时，最终形成处于低能状态的晶体。
如果用粒子的能量定义材料的状态， Metropolis 算法用一个简单的数学模型描述了
退火过程。假设材料在状态i 之下的能量为 E(i) ，那么材料在温度T 时从状态i 进入状态 j 就遵循如下规律：

p是状态Si转为Sj的概率，如果如第一行，就表示一定会变为状态Sj；相反如果Sj的目标函数小于Si，则以一定概率由Si变为Sj，KB为物理学中玻尔兹曼常数，T为材料温度。
过程：1）解空间：求初始解
2）目标函数（代价函数）
3）新解产生：变换等
4）代价函数差
5）接受准则（p）
6）降温准则：T=alpha*T
7）结束：终止温度
关于TSP旅行商问题的求解:
问题就是求从一个图的某点出发遍历一圈儿回到起点的最短路径！！！
首先随机生成一个初始解，目标函数就是总路径长度，设置降温步长是100，看代码吧：

temperature = 1000;                 % Initialize the temperature.
cooling_rate = 0.94;                % cooling rate
iterations = 1;                     % Initialize the iteration number.

% Initialize random number generator with "seed". 
rand('seed',0);                    

route = randperm(numberofcities);%路径格式：[1,2...n]
%目标函数：
previous_distance = totaldistance(route,dis);%根据路径顺序计算总距离

%100次迭代后冷却当前温度的标志
temperature_iterations = 1;

while 1.0 < temperature
    %生成随机临解：
    temp_route = perturb(route,'reverse');
    % 计算当前解的路径长
    current_distance = totaldistance(temp_route, dis);
    % 计算距离改变
    diff = current_distance - previous_distance;
    
    % Metropolis Algorithm
    if (diff < 0) || (rand < exp(-diff/(temperature)))
        route = temp_route;         %接受新路径
        previous_distance = current_distance;
        
        % 更新迭代
        temperature_iterations = temperature_iterations + 1;
        plot_iterations = plot_iterations + 1;
        iterations = iterations + 1;
    end
    
    % 每100次迭代降温一次
    if temperature_iterations >= 100
       temperature = cooling_rate*temperature;
       temperature_iterations = 0;
    end
end

2.遗传算法

遗传算法（ Genetic Algorithms，简称 GA）是一种基于自然选择原理和自然遗传机制的搜索（寻优）算法，它是模拟自然界中的生命进化机制，在人工系统中实现特定目标的优化。遗传算法的实质是通过群体搜索技术，根据适者生存的原则逐代进化，最终得到最优解或准最优解。它必须做以下操作：初始群体的产生、求每一个体的适应度、根据适者生存的原则选择优良个体、被选出的优良个体两两配对，通过随机交叉其染色体的基因并随机变异某些染色体的基因后生成下一代群体，按此方法使群体逐代进化，直到满足进化终止条件。其实现方法如下：
（ 1） 根据具体问题确定可行解域，确定一种编码方法，能用数值串或字符串表示
可行解域的每一解。（说的好抽象啊，其实就是你所解决问题的解，把你的解的形式变成一种其他特定的形式！！！）
（ 2） 对每一解应有一个度量好坏的依据，它用一函数表示，叫做适应度函数（类似于前面的能量E，目标函数/代价函数！！！）。
（ 3） 确定进化参数群体规模 M 、交叉概率pc 、变异概率 pm 、进化终止条件。
染色体为一个解，种群为一组染色体，也就是一组解，基因是染色体上的一个单元，即解的一个参数。
例子依然是TSP问题：
main.m

popSize = 100;                      % 种群规模
max_generation = 1000;              % 代数
probmutation = 0.16;                % 突变率
 
rand('seed',103)
% 随机初始化种群，可以看到即生成100个解
pop = zeros(popSize,numberofcities); 
for i=1:popSize
    pop(i,:)=randperm(numberofcities); 
    %这里其实没用什么特别的编码形式，就是34个城市编号的一个序列
end

for generation = 1:max_generation   % 代数循环
    
    % 评价: 对于种群中的每个个体计算适应度(1/总距离)，因为适应度一般应设计为越大越好，所以取倒数！！！ 
    popDist = totaldistance(pop,dis); %dis是已知的距离矩阵，返回1*100（种群个数）的矩阵
    fitness = 1./popDist;
   
    % 寻找最佳路径（然而这里并没有用适应度函数呵呵！！！）
    [mindist, bestID] = min(popDist); 
    bestPop = pop(bestID, :);       % best route
    
    %选择
    pop = select(pop, fitness, popSize,'competition');
    
    % 交叉
    pop = crossover(pop);
    
    % 变异
    pop = mutation(pop, probmutation);
   
    % 保存精英(best path) 直接放入下一代！！！
    pop = [bestPop; pop];
end

% 返回结果
[mindist, bestID]=min(popDist); 
bestPop = pop(bestID, :);

选择

运算的使用是对个体进行优胜劣汰: 从父代群体中选取一些适应度高的个体，遗传到下一代群体。
常用：轮盘赌，按比例选择算子
两两竞争，从父代中随机选择两个根据适应度比较
select.m

function popselected = select(pop, fitness, nselected, method)

popSize = size(pop,1);

switch method
    case 'roulette' %轮盘赌
        p=fitness/sum(fitness); % 选择概率
        cump=cumsum(p);         % 概率累加
        I = interp1([0 cump],1:(popSize+1),rand(1,nselected),'linear');
        % 线性插值
        I = floor(I); %向下取整
    case 'competition' %两两竞争
        % randomly generated two sets of population
        i1 = ceil( popSize*rand(1,nselected) );
        i2 = ceil( popSize*rand(1,nselected) );
        % compare the fitness and select the fitter
        I = i1.*( fitness(i1)>=fitness(i2) ) + ...
            i2.*( fitness(i1)< fitness(i2) );
        
end
popselected=pop(I,:);

交叉

交叉略微复杂一点，这张图是过程中的一步，大概就是首先要把父1整个换为子2，父2换为子1，然后把2个子代的中段从左到右依次和父代对应位置交换，每换一个子代中就会有一对儿重复的，保留中段里的，把两边那个重复的替换成中段没交换前的。不知道我讲没讲清，但我不太明白为啥要搞这么复杂。简单来说还可以选择
直接把中段都换掉，或者选择单点交叉，也就是把左上和右上换一下，左下和右上再换一下。
crossover.m

function children = crossover(parents)
% Mapped Crossover (PMX) example:     
%           _                          _                          _
%    [1 2 3|4 5 6 7|8 9]  |-> [4 2 3|1 5 6 7|8 9]  |-> [4 2 3|1 8 6 7|5 9]
%    [3 5 4|1 8 7 6|9 2]  |   [3 5 1|4 8 7 6|9 2]  |   [3 8 1|4 5 7 6|9 2]
%           |             |            |           |              |            
%           V             |            V           |              |  
%    [* 2 3|1 5 6 7|8 9] _|   [4 2 3|1 8 6 7|* 9] _|              V
%    [3 5 *|4 8 7 6|9 2]      [3 * 1|4 5 7 6|9 2]           ... ... ...
%

[popSize, numberofcities] = size(parents);    
children = parents; % childrens

for i = 1:2:popSize % pairs counting
    parent1 = parents(i+0,:);  child1 = parent1;
    parent2 = parents(i+1,:);  child2 = parent2;
    %随机选择两个交叉点
    InsertPoints = sort(ceil(numberofcities*rand(1,2)));
    for j = InsertPoints(1):InsertPoints(2)
        if parent1(j)~=parent2(j)
            child1(child1==parent2(j)) = child1(j);
            child1(j) = parent2(j);
            
            child2(child2==parent1(j)) = child2(j);
            child2(j) = parent1(j);
        end
    end
    % 两个子代
    children(i+0,:)=child1;     children(i+1,:)=child2;
end

变异我理解上要简单一点，大概就是对换几个点之类的吧

function children = mutation(parents, probmutation)
%     
% swap:    _         _    slide:    _ _________    flip:     ---------->
%     [1 2|3 4 5 6 7 8|9]      [1 2|3 4 5 6 7 8|9]      [1 2|3 4 5 6 7 8|9] 
%                                   _________ _              <----------
%     [1 2|8 4 5 6 7 3|9]      [1 2|4 5 6 7 8 3|9]      [1 2|8 7 6 5 4 3|9]
%

[popSize, numberofcities] = size(parents);
children = parents;
for k=1:popSize
    if rand < probmutation
       InsertPoints = ceil(numberofcities*rand(1,2));
       I = min(InsertPoints);  J = max(InsertPoints);
       switch ceil(rand*6) %这里用了三种
           case 1    % swap %交换
             children(k,[I J]) = parents(k,[J I]);
           case 2    % slide #滑动
             children(k,[I:J]) = parents(k,[I+1:J I]);
           otherwise % flip #翻转吧
             children(k,[I:J]) = parents(k,[J:-1:I]);
       end
    end
end

大概就是这样。