启发式算法-2012课件

1、启发式算法四川大学数学学院谭英谊Email: 一、组合优化问题二、启发式算法三、模拟退火算法四、遗传算法解决离散的优化问题运筹学分支。通过数学方法的研究去寻找离散事件的最优编排、分组、次序或筛选等，可以涉及信息技术、经济管理、工业工程、交通运输和通信网络等许多方面。一.组合优化问题1.1 组合优化问题的基本概念一般模型目标函数约束函数0-1背包问题（Knapsack Problem）加工调度问题（Scheduling Problem）旅行商问题（Travelling Salesman Problem-TSP）装箱问题（Bin Packing Problem）图着色问题（Graph Colori

2、ng Problem）经典的组合优化问题：0-1背包问题一个旅行者要从n种物品中选取b公斤重的物品，每个物品至多选一件。问这个旅行者应该怎样选取，是所选物品的总价值最大。旅行商问题给定n个城市和每两个城市间的距离。一个货郎自某一城市出发巡回售货，问这个货郎应该如何选择路线，使每个城市经过一次且仅一次，并且路径长度最短。距离矩阵，路径，目标函数例：TSP的近似算法：最近邻法、最近插入法、最小支撑树法、局部搜索法等算法（最近邻法）：(1)任选一个城市c1为出发点；(2)根据与出发城市最近的原则，在c1以外的n-1个城市中选取第二个经过城市；(3)再以c2为出发城市，在尚未经过的按最近邻原则选取c3

3、，如此重复，直至所有城市都被经过，最后到达的城市是cn；（4）从cn返回c1，构成一条路径。例：0-1背包问题：贪婪算法每一步只在可选输入中选取一个满足约束条件的输入来构造一个可行解，实现某个优化测度（目标函数或别的测度）下的最优。后继步对前趋步的结果进行扩充，直到满足约束条件的输入选完，不能再扩充为止。1.3 邻域结构与局部最优1.4 局部搜索算法改善局部搜索算法性能的途径：(1)对大量初始解执行算法，再从中选优(2)引入更复杂的邻域结构，使算法能对解空间的更大范围进行搜索(3)改变局部搜索算法只接受优化解迭代的准则，在一定限度接受恶化解。3.2 模拟退火算法流程：1）随机产生一个初始解 ，

4、令 ，并计算目标函数值2）设置初始温度 t=T(0)，终止温度 降温系数3）while tTmini)for j=1:kii)对当前最优解 按照某一邻域函数,产生一新的解 计算新的目标函数值 并计算目标函数值的增量iii) 如果 则iv）如果 则如果否则v）end for4）降温5）end while6）输出当前最优点，计算结束。3.3 模拟退火算法要素1）状态空间与状态生成函数（邻域函数）搜索空间（状态空间）：由经过编码的可行解的集合所组成状态生成函数（邻域函数）应尽可能保证产生的候选解遍布全部解空间。通常由两部分组成，即产生候选解的方式和候选解产生的概率分布候选解一般采用按照某一概率密度函

5、数对解空间进行随机采样来获得概率分布可以是均匀分布、正态分布、指数分布等等2）状态转移概率（接受概率）p状态转移概率是指从一个状态xold（一个可行解）向另一个状态xnew（另一个可行解）的转移概率通俗的理解是接受一个新解为当前解的概率它与当前的温度参数T有关，随温度下降而减小一般采用Metropolis准则3.4 模拟退火算法实验性能分析：1）模拟退火法与局部搜索算法的差异2）优点：高效、健壮、通用、灵活3）不足：返回一个高质量近似解的时间花费较多，当问题规模不可避免增大时，难于承受的运行时间将使算法丧失可行性。如何改进？3.5 模拟退火算法的matlab实现利用模拟退火算法求函数极小：无约

6、束或有界约束x = simulannealbnd(fun,x0)x = simulannealbnd(fun,x0,lb,ub)x = simulannealbnd(fun,x0,lb,ub,options)x,fval = simulannealbnd(.)x,fval,exitflag = simulannealbnd(.)x,fval,exitflag,output = simulannealbnd(.)dejong5fcnx0 = 0 0;x,fval = simulannealbnd(dejong5fcn,x0)Optimization terminated: change in b

7、est function value less than options.TolFun.x = 0.0216 -31.9955fval = 2.9821Optimization terminated: change in best function value less than options.TolFun.x = -15.9669 -31.9749fval = 1.9920exitflag = 1output = iterations: 1608 funccount: 1621 message: 1x80 char randstate: 625x1 uint32 randnstate: 2

8、x1 double problemtype: unconstrained temperature: 2x1 double totaltime: 0.8268x,fval,exitflag,output = simulannealbnd(dejong5fcn,0,0)x0 = 0 0;lb = -64 -64;ub = 64 64;x,fval = simulannealbnd(dejong5fcn,x0,lb,ub)Optimization terminated: change in best function value less than options.TolFun.x = -32.01

9、69 -31.9879fval = 0.9980fun = (x) 3*sin(x(1)+exp(x(2);x = simulannealbnd(fun,1;1,0 0)Optimization terminated: change in best function value less than options.TolFun.x = 896.9234 0.0000options = saoptimset(PlotFcns,saplotbestf,saplottemperature,saplotf,saplotstopping);simulannealbnd(dejong5fcn,x0,lb,

10、ub,options);options = saoptimset(InitialTemperature,300 50);options = saoptimset(options,TemperatureFcn,temperaturefast);options = saoptimset(options,ReannealInterval,50);options = saoptimset(options,Display,iter,DisplayInterval,400);options = saoptimset(options,TolFun,1e-5);inputcities = 1150.0 630

11、.0 40.0 750.0 750.0 1030.0 1650.0 1490.0 . 790.0 710.0 840.0 1170.0 970.0 510.0 750.0 1280.0 230.0 460.0 . 1040.0 590.0 830.0 490.0 1840.0 1260.0 1280.0 490.0 1460.0 . 1260.0 360.0; 1760.0 1660.0 2090.0 1100.0 2030.0 2070.0 650.0 . 1630.0 2260.0 1310.0 550.0 2300.0 1340.0 700.0 900.0 1200.0 . 590.0

12、860.0 950.0 1390.0 1770.0 500.0 1240.0 1500.0 790.0 . 2130.0 1420.0 1910.0 1980.0;2）计算总距离的函数该函数根据给定路径计算该路径对应总路程。function d = distance(inputcities)% d = distance(inputcities)计算TSP问题中n个城市间某路径的总路程% 输入参数为2*n数组，其中n为城市数目，列为对应城市坐标 d = 0; for n = 1 : length(inputcities) if n = length(inputcities) d = d + no

13、rm(inputcities(:,n) - inputcities(:,1); else d = d + norm(inputcities(:,n) - inputcities(:,n+1); end end3）绘制路线的函数function f = plotcities(inputcities)% plotcities(inputcities)从inputcities参数中绘制城市的位置%输入参数为2*n数组，其中n为城市数目，列为对应城市坐标temp_1 = plot(inputcities(1,:),inputcities(2,:),b*);set(temp_1,erasemode,no

14、ne);temp_2 = line(inputcities(1,:),inputcities(2,:),Marker,*);set(temp_2,color,r);x = inputcities(1,1) inputcities(1,length(inputcities);y = inputcities(2,1) inputcities(2,length(inputcities);x1 = 10*round(max(inputcities(1,:)/10);y1 = 10*round(max(inputcities(2,:)/10);if x1 = 0 x1 = 1;endif y1 = 0

15、y1 = 1;endaxis(0 x1 0 y1);temp_3 = line(x,y);set(temp_3,color,r);dist = distance(inputcities);distance_print = sprintf(.The roundtrip length for %d cities is % 4.6f units. ,length(inputcities),dist);text(x1/15,1.05*y1,distance_print,fontweight,bold);drawnow;4）交换城市的函数% s = swapcities(inputcities,n) n

16、个城市随机的交换返回一组m个城市s = inputcities;for i = 1 : n city_1 = round(length(inputcities)*rand(1); if city_1 1 city_1 = 1; end city_2 = round(length(inputcities)*rand(1); if city_2 1 city_2 = 1; end temp = s(:,city_1); s(:,city_1) = s(:,city_2); s(:,city_2) = temp;end5）执行模拟退火的函数function s = simulatedannealin

17、g(inputcities,initial_temperature,.cooling_rate,threshold,numberofcitiestoswap)% s = simulatedannealing (inputcities,initial_temperature,cooling_rate,)通过% n 个城市的TSP问题的最优解返回一个新的城市构型。global iterations;temperature = initial_temperature;initial_cities_to_swap = numberofcitiestoswap;iterations = 1;comple

18、te_temperature_iterations = 0;while iterations threshold previous_distance = distance(inputcities); temp_cities = swapcities(inputcities,numberofcitiestoswap); current_distance = distance(temp_cities); diff = abs(current_distance - previous_distance);if current_distance = 10 temperature = cooling_ra

19、te*temperature; complete_temperature_iterations = 0; endnumberofcitiestoswap = round(numberofcitiestoswap*exp(-diff/(iterations*temperature); if numberofcitiestoswap = 0 numberofcitiestoswap = 1; end iterations = iterations + 1; complete_temperature_iterations = complete_temperature_iterations + 1;

20、else if rand(1) exp(-diff/(temperature) inputcities = temp_cities; plotcities(inputcities); numberofcitiestoswap = round(numberofcitiestoswap*exp(-diff/(iterations*temperature); if numberofcitiestoswap = 0 numberofcitiestoswap = 1; end complete_temperature_iterations = complete_temperature_iteration

21、s + 1; iterations = iterations + 1; end end clc fprintf(tttIterations = %dn,iterations); fprintf(tttTemperature = %3.8fn,temperature); fprintf(tttIn current temperature for %d timesn,complete_temperature_iterations);endprevious_distancesimulatedannealing(inputcities,2000,0.97,2000,2)四. 遗传算法4.1 遗传算法简

22、介19世纪60年代 Holland一族通过模拟自然进化过程搜索最优解的方法遗传算法是一种概率搜索算法，它是利用某种编码技术作用于称为染色体的二进制数串，其基本思想是模拟由这些串组成的群体的进化过程。遗传算法通过有组织地然而是随机的信息交换来重新结合那些适应性好的串，在每一代中，利用上一代串结构中适应性好的位和段来生成一个新的串的群体；作为额外添加，偶尔也要在串结构中尝试用新的位和段来替代原来的部分。与自然界相似，遗传算法对求解问题的本身一无所知，它所需要的仅是对算法所产生的每个染色体进行评价，并基于适应值来选择染色体，是适应性好的染色体比适应性差的染色体有更多的繁殖机会。遗传算法利用简单的编码

23、技术和繁殖机制来表现复杂的现象，从而解决非常困难的问题。特别是它不受搜索空间的限制性假设的约束，不必要求诸如连续性、导数存在、单峰等假设，以及其固有的并行性。2）初始化：设置进化代数计数器 t=0，设置最大进化代数T，随机生成M个个体作为初始种群P(0)5 2 1 4 3 52 4 3 1 5 23 4 2 1 5 31 4 5 2 3 1例如：1）编码：将解空间的解数据表示成遗传空间的基因型串结构数据；4.2 遗传算法流程3）个体评价：计算种群P(t)中各个个体的适应度4）选择运算：将选择算子作用于群体。选择的目的是把优化的个体直接遗传到下一代或通过配对交叉产生新的个体再遗传到下一代。选择操

24、作是建立在群体中个体的适应度评估基础上的。经典的选择算子：轮盘赌被选中的概率为：5）交叉运算：将交叉算子作用于群体。所谓交叉是指把两个父代个体的部分结构加以替换重组而生成新个体的操作。遗传算法中起核心作用的就是交叉算子5 2 | 1 4 | 3 52 4 | 3 1 | 5 23 4 | 2 1 | 5 31 4 | 5 2 | 3 1例：有序杂交法：将个体两两配对，并以交叉概率Pc把配对的父代个体加以替换重组而生成新个体5 2 | 3 1 | 3 52 4 | 1 4 | 5 2倘若交换基因后出现了重复，则将该城市取代字串B1与A1中的城市具有相同位置的新城市，直到与字串A1中的城市均不重复

25、为止。5 2 | 3 1 | 1 52 4 | 1 4 | 5 25 2 | 3 1 | 4 52 4 | 1 4 | 5 26）变异运算：将变异算子作用于群体。即是对群体中的个体串的某些基因座上的基因作变动。以变异概率Pm对群体中个体串某些基因位上的基因值作变动，若变异后子代的适应度值更加优异，则保留子代染色体，否则，仍保留父代染色体。例：5 2 | 3 1 | 4 55 2 | 1 3 | 4 5群体P(t)经过选择、交叉和变异运算之后得到下一代群体P(t+1)7）终止条件判断：若t=T，则以进化过程中所得到的具有最大适应度个体作为最优解输出，终止计算。4.3 遗传算法的matlab实现x

26、 = ga(fitnessfcn,nvars)x = ga(fitnessfcn,nvars,A,b)x = ga(fitnessfcn,nvars,A,b,Aeq,beq)x = ga(fitnessfcn,nvars,A,b,Aeq,beq,LB,UB)x = ga(fitnessfcn,nvars,A,b,Aeq,beq,LB,UB,nonlcon)x = ga(fitnessfcn,nvars,A,b,Aeq,beq,LB,UB,nonlcon,options)x, fval = ga(.)x, fval, exitflag, output, population, scores =

27、ga()A = 1 1; -1 2; 2 1;b = 2; 2; 3;lb = zeros(2,1);x, fval, exitflag = ga(lincontest6,2,A,b,lb)Optimization terminated: average change in the fitness value less than options.TolFun.x = 0.7904 1.2094fval = -8.0180exitflag = 1function y = simple_fitness(x) y = 100 * (x(1)2 - x(2) 2 + (1 - x(1)2;functi

28、on c, ceq = simple_constraint(x) c = 1.5 + x(1)*x(2) + x(1) - x(2); -x(1)*x(2) + 10; ceq = ;ObjectiveFunction = simple_fitness;nvars = 2; % Number of variablesLB = 0 0; % Lower boundUB = 1 13; % Upper boundConstraintFunction = simple_constraint;x,fval = ga(ObjectiveFunction,nvars,LB,UB, ConstraintFu

29、nction)Optimization terminated: average change in the fitness value less than options.TolFun and constraint violation is less than options.TolCon.x = 0.8122 12.3122fval = 1.3578e+004gatoolFitness function：ackleyfcnNumber of variables：10选中Plots中的Best fitnessX = gamultipoj(fitnessfcn,nvars,A,b,Aeq,beq

30、,lb,ub,options)可以处理线性不等式约束、线性等式约束和有界约束解多目标规划function y = simple_multiobjective(x) y(1) = (x+2)2 - 10; y(2) = (x-2)2 + 20;FitnessFunction = simple_multiobjective;numberOfVariables = 1;x,fval = gamultiobj(FitnessFunction,numberOfVariables);Optimization terminated: average change in the spread of Paret

31、o solutions less than options.TolFun.size(x)ans = 15 1size(fval)ans = 15 24.4 遗传算法对TSP的应用load(usborder.mat,x,y,xx,yy);plot(x,y,Color,red); hold on;cities = 40;locations = zeros(cities,2);n = 1;while (n = cities) xp = rand*1.5; yp = rand; if inpolygon(xp,yp,xx,yy) locations(n,1) = xp; locations(n,2) = yp; n = n+1; endendplot(locations(:,1),locations(:,2),bo);distances = zeros(cities);for count1=1:cities, for count2=1:count1, x1 = locations(count1,1); y1 = locations(count1,2); x2 = loca