一、MATLAB全局优化概况

MATLAB中有个全局优化工具箱，该工具箱集成了几个主流的全局优化算法，包含全局搜索、多初始点、模式搜索、遗传算法、多目标遗传算法、模拟退火求解器和粒子群求解器， 如图所示。对于目标函数或约束函数连续、不连续、随机、导数不存在以及包含仿真或黑箱函数的优化问题，都可使用这些求解器来求解。 另外，还可通过设置选项和自定义创建、更新函数来改进求解器效率。可以使用自定义数据类型，配合遗传算法和模拟退火求解器，来描绘采用标准数据类型不容易表达的问题。利用混合函数选项，可在第一个求解器之后应用第二个求解器来改进解算。

二、遗传算法

1.基本过程

2.遗传算法步骤

(1)初始参数

(2)染色体编码

由种群个体的表现型集合所组成的空间称为问题空间，由种群基因型个体所组成的空间称为编码空间。由问题空间向编码空间的映射称作编码，而由编码空间向问题空间的映射称为解码。常用的编码方式有：二进制编码和浮点数（实数）编码。

(3)适应度函数

适应度函数是用来衡量个体优劣，度量个体适应度的函数。适应度函数值越大的个体越好，反之，适应值越小的个体越差。在遗传算法中根据适应值对个体进行选择，以保证适应性能好的个体有更多的机会繁殖后代，使优良特性得以遗传。

(4)约束条件的处理

在遗传算法中必须对约束条件进行处理，但目前尚无处理各种约束条件的一般方法。根据具体问题，可选择下列三种方法：罚函数法、搜索空间限定法、可行解变换法。

(5)遗传算子

在遗传算法中必须对约束条件进行处理，但目前尚无处理各种约束条件的一般方法。根据具体问题，可选择下列三种方法：罚函数法、搜索空间限定法、可行解变换法。 遗传算法中包含了3个模拟生物基因遗传操作的遗传算子：选择（复制）、交叉（重组）和变异（突变）。遗传算法利用遗传算子产生新一代群体来实现群体进化，算子的设计是遗传策略的主要组成部分，也是调整和控制进化过程的基本工具。

(6)搜索终止条件

遗传算法的终止条件有以下两个，满足任何一个条件，搜索就结束。 （1） 遗传操作中连续多次前后两代群体中最优个体的适应度相差在某个任意小的正数x所确定的范围内。 （2）达到遗传操作的最大进化代数t。

(7)、遗传算法的实例

例题1 解题步骤： （1）用MATLAB编写一个命名为simple_fitness.m的函数。 （2）对于约束条件，同样可以创建一个命名为simple_constraint.m的函数来表示。 （3）直接调用ga函数来实现用遗传算法对以上优化问题的求解。 结果： 例题2(最短路径之旅行商问题) function main a=zeros(6); a(1,2)=56;a(1,3)=35;a(1,4)=21;a(1,5)=51;a(1,6)=60; a(2,3)=21;a(2,4)=57;a(2,5)=78;a(2,6)=70; a(3,4)=36;a(3,5)=68;a(3,6)=68; a(4,5)=51;a(4,6)=61; a(5,6)=13; a=a+a’; L=size(a,1); c=[5 1:4 6 5]; %选取初始圈 [circle,long]=modifycircle(a,L,c) %调用下面修改圈的子函数 %******************************************* %以下为修改圈的子函数 %******************************************* function [circle,long]=modifycircle(a,L,c); for k=1:L flag=0; %退出标志 for m=1:L-2 %m为算法中的i for n=m+2:L %n为算法中的j if a(c(m),c(n))+a(c(m+1),c(n+1))<a(c(m),c(m+1))+a(c(n),c(n+1)) c(m+1:n)=c(n: -1:m+1); flag=flag+1; %修改一次，标志加1 end end end if flag==0 %一条边也没有修改,就返回 long=0; %圈长的初始值 for i=1:L long=long+a(c(i),c(i+1)); %求改良圈的长度 end circle=c; %返回修改圈 return end end 总结： 1.善于用对称矩阵读取数据

三、模拟退火算法

1.原理

模拟退火算法（Simulated Annealing，SA）是一种通用概率算法，用来在一个大的搜寻空间内寻找问题的最优解。 模拟退火算法的思想源于固体的退火过程：将固体加热至足够高的温度，再缓慢冷却；升温时，固体内部粒子随温度升高变为无序状，内能增大，而缓慢冷却使粒子又逐渐趋于有序。从理论上讲，如果冷却过程足够缓慢，那么冷却中任一温度时，固体都能达到热平衡，而冷却到低温时，将达到这一低温下的内能最小状态。物理退火过程和模拟退火算法的类比关系如图所示。

2.步骤

(1)符号说明 (2)算法基本步骤

(3)算法说明 1)状态表达。 状态表达即指：实际问题的解（即状态）如何以一种合适的数学形式被表达出来，它应当适用于SA的求解，又能充分表达实际问题，这需要仔细地设计。 2)新解的产生 新解产生机制的基本要求是能够尽量遍及解空间的各个区域，这样，在某一恒定温度不断产生新解时，就可能跳出当前区域以搜索其他区域。 3)收敛的一般性条件 收敛到全局最优的一般性条件是：初始温度足够高；热平衡时间足够长；终止温度足够低；降温过程足够缓慢。 4)参数的选择。控制参数 T的初值T0；控制参数T的衰减函数；Markov链长度。 5)算法停止准则。 一般Tf应设为一个足够的小的正数，比如0.01～5，但这只是一个粗糙的经验，更精细的设置及其他的终止准则需要根据具体问题做进一步的研究后再设定。

3.实例

例题1 算法设计步骤： (1) TSP问题的解空间和初始解

(2) 目标函数 (3) 新解的产生 新解可通过分别或者交替使用以下两种方法来产生：二变换法、三变换法。 (4) 目标函数差 (5) Metropolis接受准则 例题2 TSPLIB是一组各类TSP问题的实例集合，这里以TSPLIB的berlin52为例进行求解，berlin52有52座城市。 答案： clear clc a = 0.99; % 温度衰减函数的参数 t0 = 97; tf = 3; t = t0; Markov_length = 10000; % Markov链长度 coordinates = [ 1 565.0 575.0; 2 25.0 185.0; 3 345.0 750.0; 4 945.0 685.0; 5 845.0 655.0; 6 880.0 660.0; 7 25.0 230.0; 8 525.0 1000.0; 9 580.0 1175.0; 10 650.0 1130.0; 11 1605.0 620.0; 12 1220.0 580.0; 13 1465.0 200.0; 14 1530.0 5.0; 15 845.0 680.0; 16 725.0 370.0; 17 145.0 665.0; 18 415.0 635.0; 19 510.0 875.0; 20 560.0 365.0; 21 300.0 465.0; 22 520.0 585.0; 23 480.0 415.0; 24 835.0 625.0; 25 975.0 580.0; 26 1215.0 245.0; 27 1320.0 315.0; 28 1250.0 400.0; 29 660.0 180.0; 30 410.0 250.0; 31 420.0 555.0; 32 575.0 665.0; 33 1150.0 1160.0; 34 700.0 580.0; 35 685.0 595.0; 36 685.0 610.0; 37 770.0 610.0; 38 795.0 645.0; 39 720.0 635.0; 40 760.0 650.0; 41 475.0 960.0; 42 95.0 260.0; 43 875.0 920.0; 44 700.0 500.0; 45 555.0 815.0; 46 830.0 485.0; 47 1170.0 65.0; 48 830.0 610.0; 49 605.0 625.0; 50 595.0 360.0; 51 1340.0 725.0; 52 1740.0 245.0; ]; coordinates(:,1) = []; % 除去第一列的城市编号 amount = size(coordinates,1); % 城市的数目 % 通过向量化的方法计算距离矩阵 dist_matrix = zeros(amount, amount); coor_x_tmp1 = coordinates(:,1) * ones(1,amount); coor_x_tmp2 = coor_x_tmp1’; coor_y_tmp1 = coordinates(:,2) * ones(1,amount); coor_y_tmp2 = coor_y_tmp1’; dist_matrix = sqrt((coor_x_tmp1-coor_x_tmp2).^2 + … (coor_y_tmp1-coor_y_tmp2).^2);

sol_new = 1:amount;         % 产生初始解

% sol_new是每次产生的新解；sol_current是当前解；sol_best是冷却中的最好解； E_current = inf;E_best = inf; % E_current是当前解对应的回路距离； % E_new是新解的回路距离； % E_best是最优解的 sol_current = sol_new; sol_best = sol_new; p = 1;

while t>=tf
    for r=1:Markov_length       % Markov链长度
        % 产生随机扰动
        if (rand < 0.5) % 随机决定是进行两交换还是三交换
            % 两交换
            ind1 = 0; ind2 = 0;
            while (ind1 == ind2)
                ind1 = ceil(rand.*amount);
                ind2 = ceil(rand.*amount);
            end
            tmp1 = sol_new(ind1);
            sol_new(ind1) = sol_new(ind2);
            sol_new(ind2) = tmp1;
        else
            % 三交换
            ind1 = 0; ind2 = 0; ind3 = 0;
            while (ind1 == ind2) || (ind1 == ind3) ...
                || (ind2 == ind3) || (abs(ind1-ind2) == 1)
                ind1 = ceil(rand.*amount);
                ind2 = ceil(rand.*amount);
                ind3 = ceil(rand.*amount);
            end
            tmp1 = ind1;tmp2 = ind2;tmp3 = ind3;
            % 确保ind1 < ind2 < ind3
            if (ind1 < ind2) && (ind2 < ind3)
                ;
            elseif (ind1 < ind3) && (ind3 < ind2)
                ind2 = tmp3;ind3 = tmp2;
            elseif (ind2 < ind1) && (ind1 < ind3)
                ind1 = tmp2;ind2 = tmp1;
            elseif (ind2 < ind3) && (ind3 < ind1) 
                ind1 = tmp2;ind2 = tmp3; ind3 = tmp1;
            elseif (ind3 < ind1) && (ind1 < ind2)
                ind1 = tmp3;ind2 = tmp1; ind3 = tmp2;
            elseif (ind3 < ind2) && (ind2 < ind1)
                ind1 = tmp3;ind2 = tmp2; ind3 = tmp1;
            end

            tmplist1 = sol_new((ind1+1):(ind2-1));
            sol_new((ind1+1):(ind1+ind3-ind2+1)) = ...
                sol_new((ind2):(ind3));
            sol_new((ind1+ind3-ind2+2):ind3) = ...
                tmplist1;
        end

        %检查是否满足约束

        % 计算目标函数值（即内能）
        E_new = 0;
        for i = 1 : (amount-1)
            E_new = E_new + ...
                dist_matrix(sol_new(i),sol_new(i+1));
        end
        % 再算上从最后一个城市到第一个城市的距离
        E_new = E_new + ...
            dist_matrix(sol_new(amount),sol_new(1));

        if E_new < E_current
            E_current = E_new;
            sol_current = sol_new;
            if E_new < E_best

% 把冷却过程中最好的解保存下来 E_best = E_new; sol_best = sol_new; end else % 若新解的目标函数值小于当前解的， % 则仅以一定概率接受新解 if rand < exp(-(E_new-E_current)./t) E_current = E_new; sol_current = sol_new; else sol_new = sol_current; end end end t=t.*a; % 控制参数t（温度）减少为原来的a倍 end

disp(最优解为：)
disp(sol_best)
disp(最短距离：)
disp(E_best)

多次执行几次上面的脚本文件，以减少其中的随机数可能带来的影响，得到的最好结果如下： 例题3 用SA算法求解经典的山峰问题 (1)定义优化问题 (2)用常规最优算法求解 [x,f] = fmincon(problem) (3)用SA算法寻找全局最小值 problem.solver = ‘simulannealbnd’; problem.objective = @(x) peaks(x(1),x(2)) + (x(1)^2 + x(2)^2 - 9); problem.options = saoptimset(‘OutputFcn’,@peaksPlotIterates,… ‘Display’,‘iter’,… ‘InitialTemperature’,10,… ‘MaxIter’,300)

[x,f] = simulannealbnd(problem) 运行结果：

四、全局优化求解器汇总