QQ登录

只需要一步,快速开始

 注册地址  找回密码

tag 标签: 算法

相关帖子

版块 作者 回复/查看 最后发表
建模用到的算法(含代码) attachment 数模资源交流 yue1010126 2012-1-8 47 11520 YP20150401 2015-4-29 20:08
模拟退火算法 模拟退火法 xttataat 2012-1-13 18 7161 空木葬花 2014-3-7 21:25
传统 模拟退火算法 源代码(VB.net) 模拟退火法 xttataat 2012-1-13 15 9832 弘道 2014-7-29 12:19
PSO算法使用简介 attachment 蚁群、粒子群算法 xttataat 2012-1-13 18 9499 shiep5f 2015-9-16 10:59
调度算法,急~~~~ 最优化算法 赵佩娴 2012-1-13 2 3870 小杰大哥 2013-1-13 10:08
本人找到了一个计算机分解大合数算法。 数学基础 632158 2012-1-18 8 4596 lansatiankong 2012-2-12 12:40
美赛有哪些常用算法啊? 美国大学生数学建模竞赛(MCM/ICM) tang19930118 2012-1-18 26 4534 tang19930118 2012-2-9 15:12
图论算法及其MATLAB 程序代码 attachment MATLAB论坛 lansatiankong 2012-1-21 7 2114 asdjh 2014-8-2 18:42
算法大全(分章节) - [!price! 2 点体力] attachment 数模资源交流 飞逝 2012-1-27 15 5720 小石头150 2012-4-27 12:00
传些有用的资料关于遗传算法 attachment 遗传算法 shorloker 2012-2-2 18 5221 一盏清茶 2014-8-24 20:18
美赛资料,一种基于蜜蜂进化型遗传算法,欢迎大家下载 attachment 美国大学生数学建模竞赛(MCM/ICM) jinfeng@1 2012-2-3 24 3058 陈剑平 2012-2-7 15:33
BP算法 模拟退火法 qicheng 2012-2-3 3 771 alair008 2012-2-3 19:45
遗传算法一般求解什么问题比较好 遗传算法 whywby001 2012-2-4 2 2301 qzwqzw 2012-8-11 16:10
【2012美赛建模组老师压轴预测】——值得一看!!! attachment 美国大学生数学建模竞赛(MCM/ICM) Brambles.Q 2012-2-4 275 23160 jayla_s 2013-1-30 14:32
求韦伯斯特算法 最优化算法 jhonguofei 2012-2-5 3 2754 SkyWalker19 2013-7-26 14:24
精通matlab6.5 attachment Matlab 资料库 视频 教程 讲义 代码 忘却依恋 2012-2-6 2 1225 alair006 2012-2-8 09:39
总结的算法大全 数学建模用到的算法基本上都有用到~希望对大家有帮助 attachment 美国大学生数学建模竞赛(MCM/ICM) xiaomiyacai 2012-2-7 17 4782 MikeyShell 2013-9-8 11:18
matlab数学建模算法全收录(只求原文作者大名) attachment 美国大学生数学建模竞赛(MCM/ICM) 天之叶子 2012-2-8 53 11069 沙林欧阳 2013-1-22 10:43
用PSO算法优化控制图参数选择 蚁群、粒子群算法 cs871106 2012-7-2 2 4384 lsy2013 2013-10-13 17:57
数学建模十大算法漫谈 数模经验分享 异一 2012-8-1 35 10621 jiajinshan110 2013-8-7 10:18

相关日志

分享 建模10大算法汇总
书香宝儿 2013-1-18 14:42
1、蒙特卡罗算法(该算法又称随机性模拟算法,是通过计算机仿真来解决问题的算法,同时可以通过模拟来检验自己模型的正确性,是比赛时必用的方法) 2、数据拟合、参数估计、插值等数据处理算法(比赛中通常会遇到大量的数据需要处理,而处理数据的关键就在于这些算法,通常使用Matlab作为工具) 3、线性规划、整数规划、多元规划、二次规划等规划类问题(建模竞赛大多数问题属于最优化问题,很多时候这些问题可以用数学规划算法来描述,通常使用Lindo、Lingo软件实现) 4、图论算法(这类算法可以分为很多种,包括最短路、网络流、二分图等算法,涉及到图论的问题可以用这些方法解决,需要认真准备) 5、动态规划、回溯搜索、分支定界等计算机算法(这些算法是算法设计中比较常用的方法,很多场合可以用到竞赛中) 6、最优化理论的三大非经典算法:模拟退火法、神经网络、遗传算法(这些问题是用来解决一些较困难的最优化问题的算法,对于有些问题非常有帮助,但是算法的实现比较困难,需慎重使用) 7、网格算法和穷举法(网格算法和穷举法都是暴力搜索最优点的算法,在很多竞赛题中有应用,当重点讨论模型本身而轻视算法的时候,可以使用这种暴力方案,最好使用一些高级语言作为编程工具) 8、一些连续离散化方法(很多问题都是实际来的,数据可以是连续的,而计算机只认的是离散的数据,因此将其离散化后进行差分代替微分、求和代替积分等思想是非常重要的) 9、数值分析算法(如果在比赛中采用高级语言进行编程的话,那一些数值分析中常用的算法比如方程组求解、矩阵运算、函数积分等算法就需要额外编写库函数进行调用) 10、图象处理算法(赛题中有一类问题与图形有关,即使与图形无关,论文中也应该要不乏图片的,这些图形如何展示以及如何处理就是需要解决的问题,通常使用Matlab进行处理)
415 次阅读|0 个评论
分享 MCMC算法中的M-H抽样
liwenhui 2012-9-19 12:53
用MATLAB编写MCMC算法产生一条平稳分布为指定分布的马尔可夫链: function =mhs(f,p0,sigma,rt) k1=1;k2=0; gser=zeros(rt,1);gser(1)=p0; while(k1rt); p1=p0+sigma*randn; q1=feval(f,p0); q0=feval(f,p1); r=q0/q1; alpha=min(1,r); urand=rand(1); ifurandalpha; gser(k1+1)=p1;p0=p1;k1=k1+1; end k2=k2+1; end mcs=gser;efp=rt/k2; f可以用M文件定义,也可以用匿名函数定义。
个人分类: 在人间|1564 次阅读|0 个评论
分享 模拟退火算法
trxiao2011 2012-9-2 15:56
模拟退火算法   模拟退火算法来源于固体退火原理,将固体加温至充分高,再让其徐徐冷却,加温时,固体内部粒子随温升变为无序状, 内能增大,而徐徐冷却时粒子渐趋有序,在每个温度都达到平衡态,最后在常温时达到基态,内能减为最小。根据Metropolis 准则,粒子在温度T时趋于平衡的概率为e-ΔE/(kT),其中E为温度T时的内能,ΔE为其改变量,k为Boltzmann常数。用固体退 火模拟组合优化问题,将内能E模拟为目标函数值f,温度T演化成控制参数t,即得到解组合优化问题的模拟退火算法:由初始 解i和控制参数初值t开始,对当前解重复“产生新解→计算目标函数差→接受或舍弃”的迭代,并逐步衰减t值,算法终止时的 当前解即为所得近似最优解,这是基于蒙特卡罗迭代求解法的一种启发式随机搜索过程。退火过程由冷却进度表(Cooling Schedule)控制,包括控制参数的初值t及其衰减因子Δt、每个t值时的迭代次数L和停止条件S。 3.5.1 模拟退火算法的模型   模拟退火算法可以分解为解空间、目标函数和初始解三部分。  模拟退火的基本思想:   (1) 初始化:初始温度T(充分大),初始解状态S(是算法迭代的起点), 每个T值的迭代次数L   (2) 对k=1,……,L做第(3)至第6步:   (3) 产生新解S′   (4) 计算增量Δt′=C(S′)-C(S),其中C(S)为评价函数   (5) 若Δt′0则接受S′作为新的当前解,否则以概率exp(-Δt′/T)接受S′作为新的当前解.   (6) 如果满足终止条件则输出当前解作为最优解,结束程序。 终止条件通常取为连续若干个新解都没有被接受时终止算法。   (7) T逐渐减少,且T-0,然后转第2步。 算法对应动态演示图: 模拟退火算法新解的产生和接受可分为如下四个步骤:   第一步是由一个产生函数从当前解产生一个位于解空间的新解;为便于后续的计算和接受,减少算法耗时,通常选择由当 前新解经过简单地变换即可产生新解的方法,如对构成新解的全部或部分元素进行置换、互换等,注意到产生新解的变换方法 决定了当前新解的邻域结构,因而对冷却进度表的选取有一定的影响。   第二步是计算与新解所对应的目标函数差。因为目标函数差仅由变换部分产生,所以目标函数差的计算最好按增量计算。 事实表明,对大多数应用而言,这是计算目标函数差的最快方法。   第三步是判断新解是否被接受,判断的依据是一个接受准则,最常用的接受准则是Metropo1is准则: 若Δt′0则接受S′作 为新的当前解S,否则以概率exp(-Δt′/T)接受S′作为新的当前解S。   第四步是当新解被确定接受时,用新解代替当前解,这只需将当前解中对应于产生新解时的变换部分予以实现,同时修正 目标函数值即可。此时,当前解实现了一次迭代。可在此基础上开始下一轮试验。而当新解被判定为舍弃时,则在原当前解的 基础上继续下一轮试验。   模拟退火算法与初始值无关,算法求得的解与初始解状态S(是算法迭代的起点)无关;模拟退火算法具有渐近收敛性,已在 理论上被证明是一种以概率l 收敛于全局最优解的全局优化算法;模拟退火算法具有并行性。 3.5.2 模拟退火算法的简单应用   作为模拟退火算法应用,讨论货郎担问题(Travelling Salesman Problem,简记为TSP):设有n个城市,用数码1,…,n代表 。城市i和城市j之间的距离为d(i,j) i, j=1,…,n.TSP问题是要找遍访每个域市恰好一次的一条回路,且其路径总长度为最 短.。   求解TSP的模拟退火算法模型可描述如下:   解空间 解空间S是遍访每个城市恰好一次的所有回路,是{1,……,n}的所有循环排列的集合,S中的成员记为(w1,w2 ,… …,wn),并记wn+1= w1。初始解可选为(1,……,n)   目标函数 此时的目标函数即为访问所有城市的路径总长度或称为代价函数:   我们要求此代价函数的最小值。   新解的产生 随机产生1和n之间的两相异数k和m,若km,则将   (w1, w2 ,…,wk , wk+1 ,…,wm ,…,wn)   变为:   (w1, w2 ,…,wm , wm-1 ,…,wk+1 , wk ,…,wn).   如果是km,则将   (w1, w2 ,…,wk , wk+1 ,…,wm ,…,wn)   变为:   (wm, wm-1 ,…,w1 , wm+1 ,…,wk-1 ,wn , wn-1 ,…,wk).   上述变换方法可简单说成是“逆转中间或者逆转两端”。   也可以采用其他的变换方法,有些变换有独特的优越性,有时也将它们交替使用,得到一种更好方法。   代价函数差 设将(w1, w2 ,……,wn)变换为(u1, u2 ,……,un), 则代价函数差为: 根据上述分析,可写出用模拟退火算法求解TSP问题的伪程序: Procedure TSPSA:  begin   init-of-T; { T为初始温度}   S={1,……,n}; {S为初始值}   termination=false;   while termination=false    begin     for i=1 to L do       begin         generate(S′form S); { 从当前回路S产生新回路S′}         Δt:=f(S′))-f(S);{f(S)为路径总长}         IF(Δt0) OR (EXP(-Δt/T)Random-of- )         S=S′;         IF the-halt-condition-is-TRUE THEN         termination=true;       End;     T_lower;    End;  End   模拟退火算法的应用很广泛,可以较高的效率求解最大截问题(Max Cut Problem)、0-1背包问题(Zero One Knapsack Problem)、图着色问题(Graph Colouring Problem)、调度问题(Scheduling Problem)等等。 3.5.3 模拟退火算法的参数控制问题   模拟退火算法的应用很广泛,可以求解NP完全问题,但其参数难以控制,其主要问题有以下三点:   (1) 温度T的初始值设置问题。   温度T的初始值设置是影响模拟退火算法全局搜索性能的重要因素之一、初始温度高,则搜索到全局最优解的可能性大,但 因此要花费大量的计算时间;反之,则可节约计算时间,但全局搜索性能可能受到影响。实际应用过程中,初始温度一般需要 依据实验结果进行若干次调整。   (2) 退火速度问题。   模拟退火算法的全局搜索性能也与退火速度密切相关。一般来说,同一温度下的“充分”搜索(退火)是相当必要的,但这 需要计算时间。实际应用中,要针对具体问题的性质和特征设置合理的退火平衡条件。   (3) 温度管理问题。   温度管理问题也是模拟退火算法难以处理的问题之一。实际应用中,由于必须考虑计算复杂度的切实可行性等问题,常采 用如下所示的降温方式: T(t+1)=k×T(t) 式中k为正的略小于1.00的常数,t为降温的次数 使用SA解决TSP问题的Matlab程序: function out = tsp(loc) % TSP Traveling salesman problem (TSP) using SA (simulated annealing). % TSP by itself will generate 20 cities within a unit cube and % then use SA to slove this problem. % % TSP(LOC) solve the traveling salesman problem with cities' % coordinates given by LOC, which is an M by 2 matrix and M is % the number of cities. % % For example: % % loc = rand(50, 2); % tsp(loc); if nargin == 0, % The following data is from the post by Jennifer Myers ( jmyers@nwu . edu) edu) % to comp.ai.neural-nets. It's obtained from the figure in % Hopfield Tank's 1985 paper in Biological Cybernetics % (Vol 52, pp. 141-152). loc = ; end NumCity = length(loc); % Number of cities distance = zeros(NumCity); % Initialize a distance matrix % Fill the distance matrix for i = 1:NumCity, for j = 1:NumCity, distance(i, j) = norm(loc(i, - loc(j, ); distance(i, j) = norm(loc(i, - loc(j, ); end end % To generate energy (objective function) from path %path = randperm(NumCity); %energy = sum(distance((path-1)*NumCity + )); % Find typical values of dE count = 20; all_dE = zeros(count, 1); for i = 1:count path = randperm(NumCity); energy = sum(distance((path-1)*NumCity + )); new_path = path; index = round(rand(2,1)*NumCity+.5); inversion_index = (min(index):max(index)); new_path(inversion_index) = fliplr(path(inversion_index)); all_dE(i) = abs(energy - ... sum(sum(diff(loc( ,)'.^2))); end dE = max(all_dE); dE = max(all_dE); temp = 10*dE; % Choose the temperature to be large enough fprintf('Initial energy = %f\n\n',energy); % Initial plots out = ; plot(loc(out(, 1), loc(out(, 2),'r.', 'Markersize', 20); axis square; hold on h = plot(loc(out(, 1), loc(out(, 2)); hold off MaxTrialN = NumCity*100; % Max. # of trials at a temperature MaxAcceptN = NumCity*10; % Max. # of acceptances at a temperature StopTolerance = 0.005; % Stopping tolerance TempRatio = 0.5; % Temperature decrease ratio minE = inf; % Initial value for min. energy maxE = -1; % Initial value for max. energy % Major annealing loop while (maxE - minE)/maxE StopTolerance, minE = inf; minE = inf; maxE = 0; TrialN = 0; % Number of trial moves AcceptN = 0; % Number of actual moves while TrialN MaxTrialN AcceptN MaxAcceptN, new_path = path; index = round(rand(2,1)*NumCity+.5); inversion_index = (min(index):max(index)); new_path(inversion_index) = fliplr(path(inversion_index)); new_energy = sum(distance( ... (new_path-1)*NumCity+ )); if rand exp((energy - new_energy)/temp), % accept it! energy = new_energy; path = new_path; minE = min(minE, energy); maxE = max(maxE, energy); AcceptN = AcceptN + 1; end TrialN = TrialN + 1; end end % Update plot out = ; set(h, 'xdata', loc(out(, 1), 'ydata', loc(out(, 2)); drawnow; % Print information in command window fprintf('temp. = %f\n', temp); tmp = sprintf('%d ',path); fprintf('path = %s\n', tmp); fprintf('energy = %f\n', energy); fprintf(' = \n', minE, maxE); fprintf(' = \n\n', AcceptN, TrialN); % Lower the temperature temp = temp*TempRatio; end % Print sequential numbers in the graphic window for i = 1:NumCity, text(loc(path(i),1)+0.01, loc(path(i),2)+0.01, int2str(i), ... 'fontsize', 8); end 又一个相关的Matlab程序 function =zkp(w,c,M,t0,tf) =size(w); L=100*n; t=t0; clear m; x=zeros(1,n) xmax=x; fmax=0; while ttf for k=1:L xr=change(x) gx=g_0_1(w,x); gxr=g_0_1(w,xr); if gxr=M fr=f_0_1(xr,c); f=f_0_1(x,c); df=fr-f; if df0 x=xr; if frfmax fmax=fr; xmax=xr; end else p=rand; if pexp(df/t) x=xr; end end end end t=0.87*t end P=fmax; X=xmax; %下面的函数产生新解 function =exchange_2(pi0,d) =size(d); clear m; u=rand; u=u*(n-2); u=round(u); if u2 u=2; end if un-2 u=n-2; end v=rand; v=v*(n-u+1); v=round(v); if v1 v=1; end v=v+u; if vn v=n; end pi_1(u)=pi0(v); pi_1(u)=pi0(u); if u1 for k=1:(u-1) pi_1(k)=pi0(k); end end if v(u+1) for k=1:(v-u-1) pi_1(u+k)=pi0(v-k); end end if vn for k=(v+1):n pi_1(k)=pi0(k); end end d_f=0; if vn d_f=d(pi0(u-1),pi0(v))+d(pi0(u),pi0(v+1)); for k=(u+1):n d_f=d_f+d(pi0(k),pi0(k-1))-d(pi0(v),pi0(v+1)); end d_f=d_f-d(pi0(u-1),pi0(u)); for k=(u+1):n d_f=d_f-d(pi0(k-1),pi0(k)); end else d_f=d(pi0(u-1),pi0(v))+d(pi0(u),pi0(1))-d(pi0(u-1),pi0(u))-d(pi0(v),pi0(1)); for k=(u+1):n d_f=d_f-d(pi0(k),pi0(k-1)); end for k=(u+1):n d_f=d_f-d(pi0(k-1),pi0(k)); end end pi_r=pi_1; 遗传算法GA 遗传算法: 旅行商问题(traveling saleman problem,简称tsp): 已知n个城市之间的相互距离,现有一个推销员必须遍访这n个城市,并且每个城市只能访问一次,最后又必须返回出发城市。如何安排他对这些城市的访问次序,可使其旅行路线的总长度最短? 用图论的术语来说,假设有一个图 g=(v,e),其中v是顶点集,e是边集,设d=(dij)是由顶点i和顶点j之间的距离所组成的距离矩阵,旅行商问题就是求出一条通过所有顶点且每个顶点只通过一次的具有最短距离的回路。 这个问题可分为对称旅行商问题(dij=dji,,任意i,j=1,2,3,…,n)和非对称旅行商问题(dij≠dji,,任意i,j=1,2,3,…,n)。 若对于城市v={v1,v2,v3,…,vn}的一个访问顺序为t=(t1,t2,t3,…,ti,…,tn),其中ti∈v(i=1,2,3,…,n),且记tn+1= t1,则旅行商问题的数学模型为: min l=σd(t(i),t(i+1)) (i=1,…,n) 旅行商问题是一个典型的组合优化问题,并且是一个np难问题,其可能的路径数目与城市数目n是成指数型增长的,所以一般很难精确地求出其最优解,本文采用遗传算法求其近似解。 遗传算法: 初始化过程:用v1,v2,v3,…,vn代表所选n个城市。定义整数pop-size作为染色体的个数,并且随机产生pop-size个初始染色体,每个染色体为1到18的整数组成的随机序列。 适应度f的计算:对种群中的每个染色体vi,计算其适应度,f=σd(t(i),t(i+1)). 评价函数eval(vi):用来对种群中的每个染色体vi设定一个概率,以使该染色体被选中的可能性与其种群中其它染色体的适应性成比例,既通过轮盘赌,适应性强的染色体被选择产生后台的机会要大,设alpha∈(0,1),本文定义基于序的评价函数为eval(vi)=alpha*(1-alpha).^(i-1) 。 选择过程:选择过程是以旋转赌轮pop-size次为基础,每次旋转都为新的种群选择一个染色体。赌轮是按每个染色体的适应度进行选择染色体的。 step1 、对每个染色体vi,计算累计概率qi,q0=0;qi=σeval(vj) j=1,…,i;i=1,…pop-size. step2、从区间(0,pop-size)中产生一个随机数r; step3、若qi-1rqi,则选择第i个染色体 ; step4、重复step2和step3共pop-size次,这样可以得到pop-size个复制的染色体。 grefenstette编码:由于常规的交叉运算和变异运算会使种群中产生一些无实际意义的染色体,本文采用grefenstette编码《遗传算法原理及应用》可以避免这种情况的出现。所谓的grefenstette编码就是用所选队员在未选(不含淘汰)队员中的位置,如: 8 15 2 16 10 7 4 3 11 14 6 12 9 5 18 13 17 1 对应: 8 14 2 13 8 6 3 2 5 7 3 4 3 2 4 2 2 1。 交叉过程:本文采用常规单点交叉。为确定交叉操作的父代,从 到pop-size重复以下过程:从 中产生一个随机数r,如果rpc ,则选择vi作为一个父代。 将所选的父代两两组队,随机产生一个位置进行交叉,如: 8 14 2 13 8 6 3 2 5 7 3 4 3 2 4 2 2 1 6 12 3 5 6 8 5 6 3 1 8 5 6 3 3 2 1 1 交叉后为: 8 14 2 13 8 6 3 2 5 1 8 5 6 3 3 2 1 1 6 12 3 5 6 8 5 6 3 7 3 4 3 2 4 2 2 1 变异过程:本文采用均匀多点变异。类似交叉操作中选择父代的过程,在rpm 的标准下选择多个染色体vi作为父代。对每一个选择的父代,随机选择多个位置,使其在每位置按均匀变异(该变异点xk的取值范围为 ,产生一个 中随机数r,该点变异为x'k=ukmin+r(ukmax-ukmin))操作。如: 8 14 2 13 8 6 3 2 5 7 3 4 3 2 4 2 2 1 变异后: 8 14 2 13 10 6 3 2 2 7 3 4 5 2 4 1 2 1 反grefenstette编码:交叉和变异都是在grefenstette编码之后进行的,为了循环操作和返回最终结果,必须逆grefenstette编码过程,将编码恢复到自然编码。 循环操作:判断是否满足设定的带数xzome,否,则跳入适应度f的计算;是,结束遗传操作,跳出。 Matlab程序: function =ga(d,termops,num,pc,cxops,pm,alpha) % %———————————————————————— % =ga(d,termops,num,pc,cxops,pm,alpha) %d:距离矩阵 %termops:种群代数 %num:每代染色体的个数 %pc:交叉概率 %cxops:由于本程序采用单点交叉,交叉点的设置在本程序中没有很好的解决,所以本文了采用定点,即第cxops,可以随机产生。 %pm:变异概率 %alpha:评价函数eval(vi)=alpha*(1-alpha).^(i-1). %bestpop:返回的最优种群 %trace:进化轨迹 %------------------------------------------------ %####@@@##版权所有!欢迎广大网友改正,改进!##@@@#### %e-mail:tobysidney33@sohu.com %#################################################### % citynum=size(d,2); n=nargin; if n2 disp('缺少变量!!') disp('^_^开个玩笑^_^') end if n2 termops=500; num=50; pc=0.25; cxops=3; pm=0.30; alpha=0.10; end if n3 num=50; pc=0.25; cxops=3; pm=0.30; alpha=0.10; end if n4 pc=0.25; cxops=3; pm=0.30; alpha=0.10; end if n5 cxops=3; pm=0.30; alpha=0.10; end if n6 pm=0.30; alpha=0.10; end if n7 alpha=0.10; end if isempty(cxops) cxops=3; end =initializega(num,citynum); for i=1:termops =f(d,t); =find(l==max(l)); trace(i)=-l(y(1)); bestpop=t(y(1),:); =select(t,l,alpha); =grefenstette(t); =crossover(g,pc,cxops); =mutation(g1,pm); %均匀变异 =congrefenstette(g); end --------------------------------------------------------- function =initializega(num,citynum) for i=1:num t(i,:)=randperm(citynum); end ----------------------------------------------------------- function =f(d,t) =size(t); for k=1:m for i=1:n-1 l(k,i)=d(t(k,i),t(k,i+1)); end l(k,n)=d(t(k,n),t(k,1)); l(k)=-sum(l(k,:)); end ----------------------------------------------------------- function =select(t,l,alpha) =size(l); t1=t; =sort(l,2);%fsort from l to u for i=1:n aftersort(i)=aftersort1(n+1-i); %change end for k=1:n; t(k,:)=t1(aftersort(k),:); l1(k)=l(aftersort(k)); end t1=t; l=l1; for i=1:size(aftersort,2) evalv(i)=alpha*(1-alpha).^(i-1); end m=size(t,1); q=cumsum(evalv); qmax=max(q); for k=1:m r=qmax*rand(1); for j=1:m if j==1r=q(1) t(k,:)=t1(1,:); elseif j~=1rq(j-1)r=q(j) t(k,:)=t1(j,:); end end end -------------------------------------------------- function =grefenstette(t) =size(t); for k=1:m t0=1:n; for i=1:n for j=1:length(t0) if t(k,i)==t0(j) g(k,i)=j; t0(j)= =crossover(g,pc,cxops) =size(g); ran=rand(1,m); r=cxops; =find(ranpc); if ru=2 for k=1:2:length(ru)-1 g1(ru(k),:)= ),g(ru(k+1), )]; g(ru(k+1),:)= ),g(ru(k), )]; g(ru(k),:)=g1(ru(k),:); end end -------------------------------------------- function =mutation(g,pm) %均匀变异 =size(g); ran=rand(1,m); r=rand(1,3); %dai gai jin rr=floor(n*rand(1,3)+1); =find(ranpm); for k=1:length(mu) for i=1:length(r) umax(i)=n+1-rr(i); umin(i)=1; g(mu(k),rr(i))=umin(i)+floor((umax(i)-umin(i))*r(i)); end end --------------------------------------------------- function =congrefenstette(g) =size(g); for k=1:m t0=1:n; for i=1:n t(k,i)=t0(g(k,i)); t0(g(k,i))= =geneticTSP(D,n,C,m,alpha) =size(D); farm=zeros(n,N);%用于存储种群 for i=1:n farm(i,:)=randperm(N);%随机生成初始种群 end R=farm(1,:);%存储最优种群 len=zeros(n,1);%存储路径长度 fitness=zeros(n,1);%存储归一化适应值 counter=0; while counterC for i=1:n len(i,1)=myLength(D,farm(i,:));%计算路径长度 end maxlen=max(len); minlen=min(len); fitness=fit(len,m,maxlen,minlen);%计算归一化适应值 rr=find(len==minlen); R=farm(rr(1,1),:);%更新最短路径 FARM=farm;%优胜劣汰,nn记录了复制的个数 nn=0; for i=1:n if fitness(i,1)=alpha*rand nn=nn+1; FARM(nn,:)=farm(i,:); end end FARM=FARM(1:nn,:); =size(FARM);%交叉和变异 while aan if nn=2 nnper=randperm(2); else nnper=randperm(nn); end A=FARM(nnper(1),:); B=FARM(nnper(2),:); =intercross(A,B); FARM= ; =size(FARM); end if aan FARM=FARM(1:n,:);%保持种群规模为n end farm=FARM; clear FARM counter=counter+1 end Rlength=myLength(D,R); function =intercross(a,b) L=length(a); if L=10%确定交叉宽度 W=1; elseif ((L/10)-floor(L/10))=randL10 W=ceil(L/10); else W=floor(L/10); end p=unidrnd(L-W+1);%随机选择交叉范围,从p到p+W for i=1:W%交叉 x=find(a==b(1,p+i-1)); y=find(b==a(1,p+i-1)); =exchange(a(1,p+i-1),b(1,p+i-1)); =exchange(a(1,x),b(1,y)); end function =exchange(x,y) temp=x; x=y; y=temp; % 计算路径的子程序 function len=myLength(D,p) =size(D); len=D(p(1,N),p(1,1)); for i=1:(N-1) len=len+D(p(1,i),p(1,i+1)); end %计算归一化适应值子程序 function fitness=fit(len,m,maxlen,minlen) fitness=len; for i=1:length(len) fitness(i,1)=(1-((len(i,1)-minlen)/(maxlen-minlen+0.000001))).^m; end 一个C++的程序: //c++的程序 #includeiostream.h #includestdlib.h templateclass T class Graph { public: Graph(int vertices=10) { n=vertices; e=0; } ~Graph(){} virtual bool Add(int u,int v,const T w)=0; virtual bool Delete(int u,int v)=0; virtual bool Exist(int u,int v)const=0; int Vertices()const{return n;} int Edges()const{return e;} protected: int n; int e; }; templateclass T class MGraph:public GraphT { public: MGraph(int Vertices=10,T noEdge=0); ~MGraph(); bool Add(int u,int v,const T w); bool Delete(int u,int v); bool Exist(int u,int v)const; void Floyd(T** d,int** path); void print(int Vertices); private: T NoEdge; T** a; }; templateclass T MGraphT::MGraph(int Vertices,T noEdge) { n=Vertices; NoEdge=noEdge; a=new T* ; for(int i=0;in;i++){ a =new T ; a =0; for(int j=0;jn;j++)if(i!=j)a =NoEdge; } } templateclass T MGraphT::~MGraph() { for(int i=0;in;i++)delete ; delete ==NoEdge)return false; return true; } templateclass T bool MGraphT::Add(int u,int v,const T w) { if(u0||v0||un-1||vn-1||u==v||a !=NoEdge){ cerr"BadInput!"endl; return false; } a =w; e++; return true; } templateclass T bool MGraphT:delete(int u,int v) { if(u0||v0||un-1||vn-1||u==v||a ==NoEdge){ cerr"BadInput!"endl; return false; } a =NoEdge; e--; return true; } templateclass T void MGraphT::Floyd(T** d,int** path) { d=new T* ; path=new int* ; for(int i=0;in;i++){ d =new T ; path =new int ; for(int j=0;jn;j++){ d =a ; if(i!=ja NoEdge)path =i; else path =-1; } } for(int k=0;kn;k++){ for(i=0;in;i++) for(int j=0;jn;j++) if(d +d d ){ d =d +d ; path =path ; } } } templateclass T void MGraphT::print(int Vertices) { for(int i=0;iVertices;i++) for(int j=0;jVertices;j++) { couta ' ';if(j==Vertices-1)coutendl; } } #define noEdge 10000 #includeiostream.h void main() { cout"请输入该图的节点数:"endl; int vertices; cinvertices; MGraphfloat b(vertices,noEdge); cout"请输入u,v,w:"endl; int u,v; float w; cinuvw; while(w!=noEdge){ //u=u-1; b.Add(u-1,v-1,w); b.Add(v-1,u-1,w); cout"请输入u,v,w:"endl; cinuvw; } b.print(vertices); int** Path; int** path=Path; float** D; float** d=D; b.Floyd(d,path); for(int i=0;ivertices;i++){ for(int j=0;jvertices;j++){ coutPath ' '; if(j==vertices-1)coutendl; } } int *V; V=new int ; cout"请输入任意一个初始H-圈:"endl; for(int n=0;n=vertices;n++){ cinV ; } for(n=0;n55;n++){ for(i=0;in-1;i++){ for(int j=0;jn-1;j++) { if(i+10ji+1jn-1){ if(D ] ]+D ] ]D ] ]+D ] ]){ int l; l=V ;V =V ;V =l; } } } } } float total=0; cout"最小回路:"endl; for(i=0;i=vertices;i++){ coutV +1' '; } coutendl; for(i=0;ivertices;i++) total+=D ] ]; cout"最短路径长度:"endl; couttotal; } C语言程序 #includestdio.h #includestdlib.h #includemath.h #includealloc.h #includeconio.h #includefloat.h #includetime.h #includegraphics.h #includebios.h #define maxpop 100 #define maxstring 100 struct pp{unsigned char chrom ; float x,fitness; unsigned int parent1,parent2,xsite; }; struct pp *oldpop,*newpop,*p1; unsigned int popsize,lchrom,gem,maxgen,co_min,jrand; unsigned int nmutation,ncross,jcross,maxpp,minpp,maxxy; float pcross,pmutation,sumfitness,avg,max,min,seed,maxold,oldrand ; unsigned char x ,y ; float *dd,ff,maxdd,refpd,fm ; FILE *fp,*fp1; float objfunc(float); void statistics(); int select(); int flip(float); int crossover(); void generation(); void initialize(); void report(); float decode(); void crtinit(); void inversion(); float random1(); void randomize1(); main() {unsigned int gen,k,j,tt; char fname ; float ttt; clrscr(); co_min=0; if((oldpop=(struct pp *)farmalloc(maxpop*sizeof(struct pp)))==NULL) {printf("memory requst fail!\n");exit(0);} if((dd=(float *)farmalloc(maxstring*maxstring*sizeof(float)))==NULL) {printf("memory requst fail!\n");exit(0);} if((newpop=(struct pp *)farmalloc(maxpop*sizeof(struct pp)))==NULL) {printf("memory requst fail!\n");exit(0);} if((p1=(struct pp *)farmalloc(sizeof(struct pp)))==NULL) {printf("memory requst fail!\n");exit(0);} for(k=0;kmaxpop;k++) oldpop .chrom ='\0'; for(k=0;kmaxpop;k++) newpop .chrom ='\0'; printf("Enter Result Data Filename:"); gets(fname); if((fp=fopen(fname,"w+"))==NULL) {printf("cannot open file\n");exit(0);} gen=0; randomize(); initialize(); fputs("this is result of the TSP problem:",fp); fprintf(fp,"city: %2d psize: %3d Ref.TSP_path: %f\n",lchrom,popsize,refpd); fprintf(fp,"Pc: %f Pm: %f Seed: %f\n",pcross,pmutation,seed); fprintf(fp,"X site:\n"); for(k=0;klchrom;k++) {if((k%16)==0) fprintf(fp,"\n"); fprintf(fp,"%5d",x ); } fprintf(fp,"\n Y site:\n"); for(k=0;klchrom;k++) {if((k%16)==0) fprintf(fp,"\n"); fprintf(fp,"%5d",y ); } fprintf(fp,"\n"); crtinit(); statistics(oldpop); report(gen,oldpop); getch(); maxold=min; fm =100.0*oldpop .x/ff; do { gen=gen+1; generation(); statistics(oldpop); if(maxmaxold) {maxold=max; co_min=0; } fm =100.0*oldpop .x/ff; report(gen,oldpop); gotoxy(30,25); ttt=clock()/18.2; tt=ttt/60; printf("Run Clock: %2d: %2d: %4.2f",tt/60,tt%60,ttt-tt*60.0); printf("Min=%6.4f Nm:%d\n",min,co_min); }while((gen100)!bioskey(1)); printf("\n gen= %d",gen); do{ gen=gen+1; generation(); statistics(oldpop); if(maxmaxold) {maxold=max; co_min=0; } fm =100.0*oldpop .x/ff; report(gen,oldpop); if((gen%100)==0)report(gen,oldpop); gotoxy(30,25); ttt=clock()/18.2; tt=ttt/60; printf("Run Clock: %2d: %2d: %4.2f",tt/60,tt%60,ttt-tt*60.0); printf("Min=%6.4f Nm:%d\n",min,co_min); }while((genmaxgen)!bioskey(1)); getch(); for(k=0;klchrom;k++) {if((k%16)==0)fprintf(fp,"\n"); fprintf(fp,"%5d",oldpop .chrom ); } fprintf(fp,"\n"); fclose(fp); farfree(dd); farfree(p1); farfree(oldpop); farfree(newpop); restorecrtmode(); exit(0); } /*%%%%%%%%%%%%%%%%*/ float objfunc(float x1) {float y; y=100.0*ff/x1; return y; } /**/ void statistics(pop) struct pp *pop; {int j; sumfitness=pop .fitness; min=pop .fitness; max=pop .fitness; maxpp=0; minpp=0; for(j=1;jpopsize;j++) {sumfitness=sumfitness+pop .fitness; if(pop .fitnessmax) {max=pop .fitness; maxpp=j; } if(pop .fitnessmin) {min=pop .fitness; minpp=j; } } avg=sumfitness/(float)popsize; } /*%%%%%%%%%%%%%%%%%%%%*/ void generation() {unsigned int k,j,j1,j2,i1,i2,mate1,mate2; float f1,f2; j=0; do{ mate1=select(); pp:mate2=select(); if(mate1==mate2)goto pp; crossover(oldpop .chrom,oldpop .chrom,j); newpop .x=(float)decode(newpop .chrom); newpop .fitness=objfunc(newpop .x); newpop .parent1=mate1; newpop .parent2=mate2; newpop .xsite=jcross; newpop .x=(float)decode(newpop .chrom); newpop .fitness=objfunc(newpop .x); newpop .parent1=mate1; newpop .parent2=mate2; newpop .xsite=jcross; if(newpop .fitnessmin) {for(k=0;klchrom;k++) oldpop .chrom =newpop .chrom ; oldpop .x=newpop .x; oldpop .fitness=newpop .fitness; co_min++; return; } if(newpop .fitnessmin) {for(k=0;klchrom;k++) oldpop .chrom =newpop .chrom ; oldpop .x=newpop .x; oldpop .fitness=newpop .fitness; co_min++; return; } j=j+2; }while(jpopsize); } /*%%%%%%%%%%%%%%%%%*/ void initdata() {unsigned int ch,j; clrscr(); printf("-----------------------\n"); printf("A SGA\n"); printf("------------------------\n"); /*pause();*/clrscr(); printf("*******SGA DATA ENTRY AND INITILIZATION *******\n"); printf("\n"); printf("input pop size");scanf("%d",popsize); printf("input chrom length");scanf("%d",lchrom); printf("input max generations");scanf("%d",maxgen); printf("input crossover probability");scanf("%f",pcross); printf("input mutation prob");scanf("%f",pmutation); randomize1(); clrscr(); nmutation=0; ncross=0; } /*%%%%%%%%%%%%%%%%%%%%*/ void initreport() {int j,k; printf("pop size=%d\n",popsize); printf("chromosome length=%d\n",lchrom); printf("maxgen=%d\n",maxgen); printf("pmutation=%f\n",pmutation); printf("pcross=%f\n",pcross); printf("initial generation statistics\n"); printf("ini pop max fitness=%f\n",max); printf("ini pop avr fitness=%f\n",avg); printf("ini pop min fitness=%f\n",min); printf("ini pop sum fit=%f\n",sumfitness); } void initpop() {unsigned char j1; unsigned int k5,i1,i2,j,i,k,j2,j3,j4,p5 ; float f1,f2; j=0; for(k=0;klchrom;k++) oldpop .chrom =k; for(k=0;klchrom;k++) p5 =oldpop .chrom ; randomize(); for(;jpopsize;j++) {j2=random(lchrom); for(k=0;kj2+20;k++) {j3=random(lchrom); j4=random(lchrom); j1=p5 ; p5 =p5 ; p5 =j1; } for(k=0;klchrom;k++) oldpop .chrom =p5 ; } for(k=0;klchrom;k++) for(j=0;jlchrom;j++) dd =hypot(x -x ,y -y ); for(j=0;jpopsize;j++) {oldpop .x=(float)decode(oldpop .chrom); oldpop .fitness=objfunc(oldpop .x); oldpop .parent1=0; oldpop .parent2=0; oldpop .xsite=0; } } /**/ void initialize() {int k,j,minx,miny,maxx,maxy; initdata(); minx=0; miny=0; maxx=0;maxy=0; for(k=0;klchrom;k++) {x =rand(); if(x maxx)maxx=x ; if(x minx)minx=x ; y =rand(); if(y maxy)maxy=y ; if(y miny)miny=y ; } if((maxx-minx)(maxy-miny)) {maxxy=maxx-minx;} else {maxxy=maxy-miny;} maxdd=0.0; for(k=0;klchrom;k++) for(j=0;jlchrom;j++) {dd =hypot(x -x ,y -y ); if(maxdddd )maxdd=dd ; } refpd=dd ; for(k=0;klchrom;k++) refpd=refpd+dd ; for(j=0;jlchrom;j++) dd =4.0*maxdd; ff=(0.765*maxxy*pow(lchrom,0.5)); minpp=0; min=dd ; for(j=0;jlchrom-1;j++) {if(dd min) {min=dd ; minpp=j; } } initpop(); statistics(oldpop); initreport(); } /**/ void report(int l,struct pp *pop) {int k,ix,iy,jx,jy; unsigned int tt; float ttt; cleardevice(); gotoxy(1,1); printf("city:%4d para_size:%4d maxgen:%4d ref_tour:%f\n" ,lchrom,popsize,maxgen,refpd); printf("ncross:%4d Nmutation:%4d Rungen:%4d AVG=%8.4f MIN=%8.4f\n\n" ,ncross,nmutation,l,avg,min); printf("Ref.cominpath:%6.4f Minpath length:%10.4f Ref_co_tour:%f\n" ,pop .x/maxxy,pop .x,ff); printf("Co_minpath:%6.4f Maxfit:%10.8f" ,100.0*pop .x/ff,pop .fitness); ttt=clock()/18.2; tt=ttt/60; printf("Run clock:%2d:%2d:%4d.2f\n",tt/60,tt%60,ttt-tt*60.0); setcolor(1%15+1); for(k=0;klchrom-1;k++) {ix=x .chrom ]; iy=y .chrom ]+110; jx=x .chrom ]; jy=y .chrom ]+110; line(ix,iy,jx,jy); putpixel(ix,iy,RED); } ix=x .chrom ]; iy=y .chrom ]+110; jx=x .chrom ]; jy=y .chrom ]+110; line(ix,iy,jx,jy); putpixel(jx,jy,RED); setcolor(11); outtextxy(ix,iy,"*"); setcolor(12); for(k=0;k1%200;k++) {ix=k+280; iy=366-fm /3; jx=ix+1; jy=366-fm /3; line(ix,iy,jx,jy); putpixel(ix,iy,RED); } printf("GEN:%3d",l); printf("Minpath:%f Maxfit:%f",pop .x,pop .fitness); printf("Clock:%2d:%2d:%4.2f\n",tt/60,tt%60,ttt-tt*60.0); } /*###############*/ float decode(unsigned char *pp) {int j,k,l; float tt; tt=dd *lchrom+pp ]; for(j=0;jlchrom-1;j++) {tt=tt+dd *lchrom+pp ];} l=0; for(k=0;klchrom-1;k++) for(j=k+1;jlchrom;j++) {if(pp ==pp )l++;} return tt+4*l*maxdd; } /*%%%%%%%%%%%%%%%%%%*/ void crtinit() {int driver,mode; struct palettetype p; driver=DETECT; mode=0; initgraph(driver,mode,""); cleardevice(); } /*$$$$$$$$$$$$$$$$$$$$*/ int select() {double rand1,partsum; float r1; int j; partsum=0.0; j=0; rand1=random1()*sumfitness; do{ partsum=partsum+oldpop .fitness; j=j+1; }while((partsumrand1)(jpopsize)); return j-1; } /*$$$$$$$$$$$$$$$*/ int crossover(unsigned char *parent1,unsigned char *parent2,int k5) {int k,j,mutate,i1,i2,j5; int j1,j2,j3,s0,s1,s2; unsigned char jj,ts1 ,ts2 ; float f1,f2; s0=0;s1=0;s2=0; if(flip(pcross)) {jcross=random(lchrom-1); j5=random(lchrom-1); ncross=ncross+1; if(jcrossj5){k=jcross;jcross=j5;j5=k;} } else jcross=lchrom; if(jcross!=lchrom) {s0=1; k=0; for(j=jcross;jj5;j++) {ts1 =parent1 ; ts2 =parent2 ; k++; } j3=k; for(j=0;jlchrom;j++) {j2=0; while((parent2 !=ts1 )(j2k)){j2++;} if(j2==k) {ts1 =parent2 ; j3++; } } j3=k; for(j=0;jlchrom;j++) {j2=0; while((parent1 !=ts2 )(j2k)){j2++;} if(j2==k) {ts2 =parent1 ; j3++; } } for(j=0;jlchrom;j++) {newpop .chrom =ts1 ; newpop .chrom =ts2 ; } } else {for(j=0;jlchrom;j++) {newpop .chrom =parent1 ; newpop .chrom =parent2 ; } mutate=flip(pmutation); if(mutate) {s1=1; nmutation=nmutation+1; for(j3=0;j3200;j3++) {j1=random(lchrom); j=random(lchrom); jj=newpop .chrom ; newpop .chrom =newpop .chrom ; newpop .chrom =jj; } } mutate=flip(pmutation); if(mutate) {s2=1; nmutation=nmutation+1; for(j3=0;j3100;j3++) {j1=random(lchrom); j=random(lchrom); jj=newpop .chrom ; newpop .chrom =newpop .chrom ; newpop .chrom =jj; } } } j2=random(2*lchrom/3); for(j=j2;jj2+lchrom/3-1;j++) for(k=0;klchrom;k++) {if(k==j)continue; if(kj){i2=k;i1=j;} else{i1=k;i2=j;} f1=dd .chrom +newpop .chrom ]; f1=f1+dd .chrom + newpop .chrom ]; f2=dd .chrom + newpop .chrom ]; f2=f2+dd .chrom + newpop .chrom ]; if(f1f2){inversion(i1,i2,newpop .chrom);} } j2=random(2*lchrom/3); for(j=j2;jj2+lchrom/3-1;j++) for(k=0;klchrom;k++) {if(k==j)continue; if(kj){i2=k;i1=j;} else{i1=k;i2=j;} f1=dd .chrom +newpop .chrom ]; f1=f1+dd .chrom + newpop .chrom ]; f2=dd .chrom + newpop .chrom ]; f2=f2+dd .chrom + newpop .chrom ]; if(f1f2){inversion(i1,i2,newpop .chrom);} } return 1; } /*$$$$$$$$$$$$$$$*/ void inversion(unsigned int k,unsigned int j,unsigned char *ss) {unsigned int l1,i; unsigned char tt; l1=(j-k)/2; for(i=0;il1;i++) {tt=ss ; ss =ss ; ss =tt; } } /*%%%%%%%%%%%%%%%*/ void randomize1() {int i; randomize(); for(i=0;ilchrom;i++) oldrand =random(30001)/30000.0; jrand=0; } /*%%%%%%%%%%%*/ float random1() {jrand=jrand+1; if(jrand=lchrom) {jrand=0; randomize1(); } return oldrand ; } /*%%%%%%%%%%*/ int flip(float probability) {float ppp; ppp=random(20001)/20000.0; if(ppp=probability)return 1; return 0; } 改进后用来求解VBR问题的delphi程序 unit uEA; interface uses uUtilsEA, uIEA, uITSP, Classes, GaPara, windows, SysUtils, fEA_TSP; type TIndividual = class(TInterfacedObject, IIndividual) private // The internally stored fitness value fFitness: TFloat; fWeConstrain: integer; fBackConstrain: integer; fTimeConstrain: integer; procedure SetFitness(const Value: TFloat); function GetFitness: TFloat; function GetWeConstrain: integer; procedure SetWeConstrain(const Value: integer); procedure SetBackConstrain(const Value: integer); function GetBackConstrain: integer; function GetTimeConstrain: integer; procedure SetTimeConstrain(const Value: integer); public property Fitness : TFloat read GetFitness write SetFitness; property WeConstrain :integer read GetWeConstrain write SetWeConstrain; property BackConstrain :integer read GetBackConstrain write SetBackConstrain; property TimeConstrain :integer read GetTimeConstrain write SetTimeConstrain; end; TTSPIndividual = class(TIndividual, ITSPIndividual) private // The route we travel fRouteArray : ArrayInt; fWeConstrain: integer; fBackConstrain: integer; fTimeConstrain: integer; function GetRouteArray(I: Integer): Integer; procedure SetRouteArray(I: Integer; const Value: Integer); procedure SetSteps(const Value: Integer); function GetSteps: Integer; function GetWeConstrain: integer; procedure SetWeConstrain(const Value: integer); procedure SetBackConstrain(const Value: integer); procedure SetTimeConstrain(const Value: integer); function GetBackConstrain: integer; function GetTimeConstrain: integer; public // Constructor, called with initial route size constructor Create(Size : TInt); reintroduce; destructor Destroy; override; property RouteArray : Integer read GetRouteArray write SetRouteArray; // The number of steps on the route property Steps : Integer read GetSteps write SetSteps; property Fitness : TFloat read GetFitness write SetFitness; property WeConstrain :integer read GetWeConstrain write SetWeConstrain; property BackConstrain :integer read GetWeConstrain write SetBackConstrain; property TimeConstrain :integer read GetTimeConstrain write SetTimeConstrain; end; TTSPCreator = class(TInterfacedObject, ITSPCreator) private // The Control component we are associated with fController: ITSPController; function GetController: ITSPController; procedure SetController(const Value: ITSPController); public // Function to create a random individual function CreateIndividual : IIndividual; function CreateFeasibleIndividual: IIndividual; property Controller : ITSPController read GetController write SetController; end; TKillerPercentage = class(TInterfacedObject, IKillerPercentage) private fPer: TFloat; procedure SetPercentage(const Value: TFloat); function GetPercentage: TFloat; public function Kill(Pop : IPopulation): Integer; // Percentage of population to be killed property Percentage: TFloat read GetPercentage write SetPercentage; end; TParentSelectorTournament = class(TInterfacedObject, IParentSelector) public function SelectParent(Population: IPopulation): IIndividual; end; TTSPBreederCrossover = class(TInterfacedObject, IBreeder) public function BreedOffspring(PSelector: IParentSelector; Pop: IPopulation): IIndividual; end; TTSPMutator = class(TInterfacedObject, ITSPMutator) private fTrans: TFloat; fInv: TFloat; procedure SetInv(const Value: TFloat); procedure SetTrans(const Value: TFloat); function GetInv: TFloat; function GetTrans: TFloat; public procedure Mutate(Individual: IIndividual); published // Probability of doing a transposition property Transposition: TFloat read GetTrans write SetTrans; // Probability of doing an inversion property Inversion: TFloat read GetInv write SetInv; end; TTSPExaminer = class(TInterfacedObject, ITSPExaminer) private // The Control component we are associated with fController: ITSPController; function GetController: ITSPController; procedure SetController(const Value: ITSPController); public // Returns the fitness of an individual as a real number where 0 = best function GetFitness(Individual : IIndividual) : TFloat; property Controller : ITSPController read GetController write SetController; end; TPopulation = class(TInterfacedObject, IPopulation) private // The population fPop : TInterfaceList; // Worker for breeding fBreeder: IBreeder; // Worker for killing fKiller: IKiller; // Worker for parent selection fParentSelector: IParentSelector; // Worker for mutation fMutator: IMutator; // Worker for initial creation fCreator: ICreator; // Worker for fitness calculation fExaminer: IExaminer; // On Change event FOnChange: TNotifyEvent; procedure Change; // Getters and Setters function GetIndividual(I: Integer): IIndividual; function GetCount: Integer; function GetBreeder: IBreeder; function GetCreator: ICreator; function GetExaminer: IExaminer; function GetKiller: IKiller; function GetMutator: IMutator; function GetOnChange: TNotifyEvent; function GetParentSelector: IParentSelector; procedure SetBreeder(const Value: IBreeder); procedure SetCreator(const Value: ICreator); procedure SetExaminer(const Value: IExaminer); procedure SetKiller(const Value: IKiller); procedure SetMutator(const Value: IMutator); procedure SetOnChange(const Value: TNotifyEvent); procedure SetParentSelector(const Value: IParentSelector); // not interfaced procedure DanQuickSort(SortList: TInterfaceList; L, R: Integer; SCompare: TInterfaceCompare); procedure Sort(Compare: TInterfaceCompare); protected // Comparison function for Sort() function CompareIndividuals(I1, I2: IIndividual): Integer; // Sort the population procedure SortPopulation; public // The constructor constructor Create; // The destructor destructor Destroy; override; // Adds an individual to the population procedure Add(New : IIndividual); // Deletes an individual from the population procedure Delete(I : Integer); // Runs a single generation procedure Generation; // Initialise the population procedure Initialise(Size : Integer); // Clear ourselves out procedure Clear; // Get the fitness of an individual function FitnessOf(I : Integer) : TFloat; // Access to the population members property Pop : IIndividual read GetIndividual; default; // The size of the population property Count : Integer read GetCount; property ParentSelector : IParentSelector read GetParentSelector write SetParentSelector; property Breeder : IBreeder read GetBreeder write SetBreeder; property Killer : IKiller read GetKiller write SetKiller; property Mutator : IMutator read GetMutator write SetMutator; property Creator : ICreator read GetCreator write SetCreator; property Examiner : IExaminer read GetExaminer write SetExaminer; // An event property OnChange : TNotifyEvent read GetOnChange write SetOnChange; end; TTSPController = class(TInterfacedObject, ITSPController) private fXmin, fXmax, fYmin, fYmax: TFloat; { The array of 'cities' } fCities : array of TPoint2D; { The array of 'vehicles' } fVehicles : array of TVehicle; { The array of 'vehicle number' } fNoVehicles : ArrayInt;///////////////////// { The number of 'new cities' } fCityCount: Integer; { The number of 'old cities' } foldCityCount: Integer; { The number of 'travelers' } fTravelCount:Integer; /////////////////////// { The number of 'depots' } fDepotCount:Integer; /////////////////////// { Getters... } function GetCity(I: Integer): TPoint2D; function GetNoVehicle(I: Integer): TInt; function GetCityCount: Integer; function GetOldCityCount: Integer; function GetTravelCount:Integer; function GetDepotCount:Integer; function GetXmax: TFloat; function GetXmin: TFloat; function GetYmax: TFloat; function GetYmin: TFloat; { Setters... } procedure SetCityCount(const Value: Integer); procedure SetOldCityCount(const Value: Integer); procedure SetTravelCount(const Value: Integer); ///////////// procedure SetDepotCount(const Value: Integer); ///////////// procedure SetXmax(const Value: TFloat); procedure SetXmin(const Value: TFloat); procedure SetYmax(const Value: TFloat); procedure SetYmin(const Value: TFloat); function TimeCostBetween(C1, C2: Integer): TFloat; function GetTimeConstraint(Individual: IIndividual): TInt; function DateSpanToMin(d1, d2: TDateTime): integer; function GetVehicleInfo(routeInt: Tint): integer; procedure writeTimeArray; procedure writeCostArray; public { The constructor } constructor Create; { The destructor } destructor Destroy; override; { Get the distance between two cities } function DistanceBetween(C1, C2 : Integer) : TFloat; { Get the cost between two cities } function CostBetween(C1, C2: Integer): TFloat; function GetWeightConstraint( Individual: IIndividual): TInt; function GetBackConstraint( Individual: IIndividual): TInt; { Places the cities at random points } procedure RandomCities; { Area limits } property Xmin: TFloat read GetXmin write SetXmin; property Xmax: TFloat read GetXmax write SetXmax; property Ymin: TFloat read GetYmin write SetYmin; property Ymax: TFloat read GetYmax write SetYmax; { Properties... } property CityCount : Integer read GetCityCount write SetCityCount; property OldCityCount : Integer read GetOldCityCount write SetOldCityCount; property TravelCount : Integer read GetTravelCount write SetTravelCount; /////////// property DepotCount : Integer read GetDepotCount write SetDepotCount; /////////// { Access to the cities array } property Cities : TPoint2D read GetCity; property NoVehicles : TInt read GetNoVehicle; /////////////// end; implementation uses Math; { TIndividual } function TIndividual.GetFitness: TFloat; begin result := fFitness; end; function TIndividual.GetWeConstrain: integer; begin result := fWeConstrain; end; function TIndividual.GetBackConstrain: integer; begin result := fBackConstrain; end; function TIndividual.GetTimeConstrain: integer; begin result := fTimeConstrain; end; procedure TIndividual.SetBackConstrain(const Value: integer); begin fBackConstrain := Value; end; procedure TIndividual.SetFitness(const Value: TFloat); begin fFitness := Value; end; procedure TIndividual.SetWeConstrain(const Value: integer); begin fWeConstrain := Value; end; procedure TIndividual.SetTimeConstrain(const Value: integer); begin fTimeConstrain := Value; end; { TTSPIndividual } constructor TTSPIndividual.Create(Size: TInt); begin Inherited Create; SetLength(fRouteArray, Size); // fSteps := Size; end; destructor TTSPIndividual.Destroy; begin SetLength(fRouteArray, 0); inherited; end; function TTSPIndividual.GetRouteArray(I: Integer): Integer; begin result := fRouteArray ; end; function TTSPIndividual.GetSteps: Integer; begin result := Length(fRouteArray); end; procedure TTSPIndividual.SetSteps(const Value: Integer); begin SetLength(fRouteArray, Value); end; procedure TTSPIndividual.SetRouteArray(I: Integer; const Value: Integer); begin fRouteArray := Value; end; function TTSPIndividual.GetWeConstrain: integer; begin result := fWeConstrain; end; function TTSPIndividual.GetBackConstrain: integer; begin result := fBackConstrain; end; function TTSPIndividual.GetTimeConstrain: integer; begin result := fTimeConstrain; end; procedure TTSPIndividual.SetWeConstrain(const Value: integer); begin fWeConstrain := Value; end; procedure TTSPIndividual.SetBackConstrain(const Value: integer); begin fBackConstrain := Value; end; procedure TTSPIndividual.SetTimeConstrain(const Value: integer); begin fTimeConstrain := Value; end; { TTSPCreator } function TTSPCreator.CreateIndividual: IIndividual; var New: ITSPIndividual; i, j, Top, Temp : Integer; //trav:integer; begin // Get the number of cities Top := fController.CityCount; // Create the new individual New := TTSPIndividual.Create(Top); // Initialise it with a sequential route for i := 0 to Top - 1 do New.RouteArray := i; // Shuffle the route for i := Top - 1 downto 1 do begin j := Random(i); Temp := New.RouteArray ; New.RouteArray := New.RouteArray ; New.RouteArray := Temp; end; result := New; end; function TTSPCreator.CreateFeasibleIndividual: IIndividual; var New: ITSPIndividual; i, j, Top, Temp : Tint; Msg:TMsg; begin // Get the number of cities Top := fController.CityCount; // Create the new individual New := TTSPIndividual.Create(Top); // Initialise it with a sequential route repeat begin////////////////////////////////// for i := 0 to Top - 1 do New.RouteArray := i; // Shuffle the route for i := Top - 1 downto 1 do begin j := Random(i); Temp := New.RouteArray ; New.RouteArray := New.RouteArray ; New.RouteArray := Temp; end; //process message sequence////////// while PeekMessage(Msg,0,0,0,1) do/// begin /// if Msg.Message18 then /// begin /// TranslateMessage(Msg); /// DispatchMessage(Msg); /// end; /// end; /// //////////////////////////////////// end until (fController.GetWeightConstraint(New)=0)and(fController.GetBackConstraint(New)=0); result := New; end; function TTSPCreator.GetController: ITSPController; begin result := fController; end; procedure TTSPCreator.SetController(const Value: ITSPController); begin fController := Value; end; { TKillerPercentage } function TKillerPercentage.GetPercentage: TFloat; begin result := fPer; end; function TKillerPercentage.Kill(Pop: IPopulation): Integer; var KillCount, i : Integer; begin // Work out the number we have to kill KillCount := Floor(Pop.Count * (fPer / 100)); // Delete the worst individuals - assuming the population is sorted for i := 1 to KillCount do Pop.Delete(Pop.Count - 1); // Return the number killed Result := KillCount; end; procedure TKillerPercentage.SetPercentage(const Value: TFloat); begin fPer := Value; end; { TParentSelectorTournament } function TParentSelectorTournament.SelectParent( Population: IPopulation): IIndividual; var i1, i2 : Integer; begin // Select a random individual i1 := Random(Population.Count); // Select a *different* random individual repeat i2 := Random(Population.Count); until i1 i2; // Hold the tournament and return the fittest of the two if Population.FitnessOf(i1) Population.FitnessOf(i2) then Result := Population else Result := Population ; end; { TTSPBreederCrossover } function TTSPBreederCrossover.BreedOffspring(PSelector: IParentSelector; Pop: IPopulation): IIndividual; var Child, Mom, Dad, Parent1, Parent2 : ITSPIndividual; i, j, p : Integer; function AlreadyAssigned(City, x : Integer) : Boolean; var y : Integer; Found : Boolean; begin Found := False; for y := 0 to x - 1 do begin if Child.RouteArray = City then begin Found := True; Break; end; end; Result := Found; end; begin // Select a some parents... Mom := PSelector.SelectParent(Pop) as ITSPIndividual; Dad := PSelector.SelectParent(Pop) as ITSPIndividual; // Create a child Child := TTSPIndividual.Create(Mom.Steps); // Copy the route from parents to child for i := 0 to Child.Steps - 1 do begin // Choose a parent at random p := Random(2); if p = 0 then begin Parent1 := Mom; Parent2 := Dad; end else begin Parent1 := Dad; Parent2 := Mom; end; if not AlreadyAssigned(Parent1.RouteArray , i) then begin // Use city from Parent 1 unless used already Child.RouteArray := Parent1.RouteArray ; end else if not AlreadyAssigned(Parent2.RouteArray , i) then begin // Otherwise use city from Parent 2 unless used already Child.RouteArray := Parent2.RouteArray ; end else begin // If both assigned already then use a random city repeat j := Random(Child.Steps); until not AlreadyAssigned(j, i); Child.RouteArray := j; end; end; // Return the child Result := Child; end; { TTSPMutator } function TTSPMutator.GetInv: TFloat; begin result := fInv; end; function TTSPMutator.GetTrans: TFloat; begin result := fTrans; end; procedure TTSPMutator.Mutate(Individual: IIndividual); var P: Double; i, j, t : Integer; Start, Finish : Integer; begin with Individual as ITSPIndividual do begin // Should we do an inversion? P := Random * 100; if P FTrans then begin // Do an inversion (i.e. swap two cities at random) // Choose first city i := Random(Steps); // Choose a second city repeat j := Random(Steps); until i j; // Swap them over t := RouteArray ; RouteArray := RouteArray ; RouteArray := t; end; // Should we do a transposition? P := Random * 100; if P FInv then begin // Do a transposition (i.e. reverse a sub-route) // Choose random start and finish points Start := Random(Steps - 1); Finish := Start + Random(Steps - Start); // Reverse the sub-route for i := 0 to Floor((Finish - Start) / 2) do begin t := RouteArray ; RouteArray := RouteArray ; RouteArray := t; end; end; end; end; procedure TTSPMutator.SetInv(const Value: TFloat); begin fInv := Value; end; procedure TTSPMutator.SetTrans(const Value: TFloat); begin fTrans := Value; end; { TTSPExaminer } function TTSPExaminer.GetController: ITSPController; begin result := fController; end; function TTSPExaminer.GetFitness(Individual: IIndividual): TFloat; var i , weightConstraint, backConstraint : TInt; Distance , penaltyW, penaltyB : TFloat; Indi : ITSPIndividual; begin Indi := Individual as ITSPIndividual; Distance := 0; penaltyW:=FormGaPara.EditWeightConstrain.Value; penaltyB:=FormGaPara.EditBackConstrain.Value; for i := 0 to Indi.Steps - 2 do begin Distance := Distance + fController.DistanceBetween(Indi.RouteArray , Indi.RouteArray ); end; Distance := Distance + fController.DistanceBetween(Indi.RouteArray , Indi.RouteArray ); WeightConstraint:=fController.GetWeightConstraint(Indi); backConstraint:=fController.GetBackConstraint(Indi); Indi.WeConstrain:=WeightConstraint; Indi.BackConstrain:=backConstraint; Result := Distance+penaltyW*weightconstraint+penaltyB*backConstraint; end; procedure TTSPExaminer.SetController(const Value: ITSPController); begin fController := Value; end; { TPopulation } constructor TPopulation.Create; begin inherited; fPop := TInterfaceList.Create; end; destructor TPopulation.Destroy; begin fPop.Free; inherited; end; procedure TPopulation.Add(New: IIndividual); begin fPop.Add(New); end; procedure TPopulation.Clear; begin fPop.Clear; end; function TPopulation.CompareIndividuals(I1, I2: IIndividual): Integer; var A, B, D : TFloat; begin // Get the difference between the two individuals (real number) A := I1.Fitness; B := I2.Fitness; D := A - B; // Quickest way to convert that to an integer is... if D 0 then Result := 1 else if D 0 then Result := -1 else Result := 0; end; procedure TPopulation.Delete(I: Integer); begin fPop.Delete(I); end; function TPopulation.FitnessOf(I: Integer): TFloat; begin result := Pop .Fitness; end; procedure TPopulation.Change; begin if Assigned(fOnChange) then FOnChange(Self); end; procedure TPopulation.Generation; var Replace, i : Integer; New : IIndividual; begin // Kill some of the population Replace := fKiller.Kill(Self); for i := 1 to Replace do begin // Breed a new individual New := fBreeder.BreedOffspring(fParentSelector, Self); // Perform some mutation on the individual FMutator.Mutate(New); // Get the fitness of the new individual New.Fitness := fExaminer.GetFitness(New); // Add it to the population Add(New); end; // Sort the population into fitness order where first == best SortPopulation; Change; end; function TPopulation.GetBreeder: IBreeder; begin result := fBreeder; end; function TPopulation.GetCount: Integer; begin result := fPop.Count; end; function TPopulation.GetCreator: ICreator; begin result := fCreator; end; function TPopulation.GetExaminer: IExaminer; begin result := fExaminer; end; function TPopulation.GetIndividual(I: Integer): IIndividual; begin result := (fPop as IIndividual); end; function TPopulation.GetKiller: IKiller; begin result := fKiller; end; function TPopulation.GetMutator: IMutator; begin result := fMutator; end; function TPopulation.GetOnChange: TNotifyEvent; begin result := fOnChange; end; function TPopulation.GetParentSelector: IParentSelector; begin result := fParentSelector; end; procedure TPopulation.Initialise(Size: Integer); var i,feasibleCount: Integer; New: IIndividual; begin feasibleCount:=round(size*(FormGaPara.EditFeasible.Value)/100); //feasibleCount:=1; // Clear out the old stuff Clear; // Set the capacity first to save about 12 nanoseconds ;o) fPop.Capacity := Size; // Create the appropriate number of individuals for i := 1 to feasibleCount do begin // Create the individual New := fCreator.CreateFeasibleIndividual; // Get the fitness of the new individual New.Fitness := fExaminer.GetFitness(New); // Add to the population Add(New); end; for i := feasibleCount+1 to Size do begin // Create the individual New := fCreator.CreateIndividual; //////// // Get the fitness of the new individual New.Fitness := fExaminer.GetFitness(New); // Add to the population Add(New); end; SortPopulation; Change; end; procedure TPopulation.SetBreeder(const Value: IBreeder); begin fBreeder := Value; end; procedure TPopulation.SetCreator(const Value: ICreator); begin fCreator := Value; end; procedure TPopulation.SetExaminer(const Value: IExaminer); begin fExaminer := Value; end; procedure TPopulation.SetKiller(const Value: IKiller); begin fKiller := Value; end; procedure TPopulation.SetMutator(const Value: IMutator); begin fMutator := Value; end; procedure TPopulation.SetOnChange(const Value: TNotifyEvent); begin fOnChange := Value; end; procedure TPopulation.SetParentSelector(const Value: IParentSelector); begin fParentSelector := Value; end; procedure TPopulation.DanQuickSort(SortList: TInterfaceList; L, R: Integer; SCompare: TInterfaceCompare); var I, J: Integer; P: IIndividual; begin repeat I := L; J := R; P := SortList.Items as IIndividual; repeat while SCompare(SortList.Items as IIndividual, P) 0 do Inc(I); while SCompare(SortList.Items as IIndividual, P) 0 do Dec(J); if I = J then begin SortList.Exchange(I, J); Inc(I); Dec(J); end; until I J; if L J then DanQuickSort(SortList, L, J, SCompare); L := I; until I = R; end; procedure TPopulation.Sort(Compare: TInterfaceCompare); begin if Assigned(fPop) and (Count 0) then DanQuickSort(fPop, 0, Count - 1, Compare); end; procedure TPopulation.SortPopulation; begin Sort(self.CompareIndividuals); end; { TTSPController } constructor TTSPController.Create; begin inherited; end; destructor TTSPController.Destroy; begin SetLength(FCities, 0); SetLength(FNoVehicles, 0); SetLength(FVehicles, 0); inherited; end; { Standard euclidian distance between two 2D vectors... } function TTSPController.DistanceBetween( C1, C2: Integer): TFloat; var temp:TFloat; i,j,intTemp,intTemp2:TInt; begin intTemp:=0; temp:=FormGaPara.EditWrongConstrain.Value; {if (Cities .id=1)and(Cities .id=fOldCityCount)and(Cities .id0) then begin fCities .serviceDepot:= fCities .serviceDepot; end; //} //8 if (Cities .id=1)and(Cities .id=fOldCityCount)and(Cities .id=1)and(Cities .id=fOldCityCount) then begin temp:=CostArray ; end; //1 if Cities .id=0 then begin for i:=0 to fDepotCount-1 do begin intTemp:=intTemp+fNoVehicles ; if Cities .id =fOldCityCount + intTemp +1 then temp:=0; end; intTemp:=0; end; //2 if Cities .id=0 then begin for i:=1 to fDepotCount do begin intTemp:=intTemp+fNoVehicles ; if Cities .id =fOldCityCount + intTemp then temp:=0; end; intTemp:=0; end; //5 for i:=0 to fDepotCount-1 do begin intTemp:=intTemp+fNoVehicles ; { if (Cities .id=fOldCityCount + intTemp +1)and(Cities .id=Cities .id+1) then temp:=10; /////////////////////////// } if (Cities .id=fOldCityCount + intTemp +1)and(Cities .id=fOldCityCount + intTemp+fNoVehicles ) and(Cities .id=fOldCityCount + intTemp +1)and(Cities .id=fOldCityCount + intTemp+fNoVehicles ) then temp:=0;//} end; intTemp:=0; //7 if (Cities .id=Cities .id)and(Cities .id fOldCityCount) then begin temp:=0; end; //3 if (Cities .id fOldCityCount)and(Cities .id=1)and(Cities .id=fOldCityCount) then begin //temp := Sqrt(sqr(Cities .X - Cities .X)+sqr(Cities .Y - Cities .Y)); temp:=CostArray ; end; //4 if (Cities .id=fOldCityCount)and(Cities .id=1)and(Cities .id fOldCityCount) then begin //if Cities .serviceDepot=Cities .serviceDepot then //back to the start point temp:=CostArray ; end; //6 intTemp:=0; for i:=1 to fDepotCount do begin intTemp:=intTemp+fNoVehicles ; if Cities .id= fOldCityCount + intTemp then begin intTemp2:=0; for j:=0 to fDepotCount-1 do begin intTemp2:=intTemp2+fNoVehicles ; if Cities .id=fOldCityCount + intTemp2 +1 then if abs(Cities .id-Cities .id) fNoVehicles -1 then temp:=0; end; //} end; end; intTemp:=0; result:=temp; end; function TTSPController.CostBetween(C1, C2: Integer): TFloat; //matrix cij var distance:TFloat; begin distance := Sqrt(sqr(Cities .X - Cities .X)+ sqr(Cities .Y - Cities .Y)); //result:=distance+TimeCostBetween(C1,C2); result:=distance; end; function TTSPController.TimeCostBetween(C1, C2: Integer): TFloat; var cost:TFloat; i,j,penaltyW,penaltyD:TFloat; startTime:TDateTime; begin startTime:=strToDateTime(FormGa.EditStartTime.Text); penaltyW:=FormGaPara.EditWaitConstrain.Value; penaltyD:=FormGaPara.EditDelayConstrain.Value; if Cities .idfOldCityCount then fCities .totalTime:=0 else fCities .totalTime:=Cities .totalTime+Cities .serviceTime+timeArray ; fCities .waitTime:= max(0,DateSpanToMin(startTime,Cities .early)-Cities .totalTime); fCities .delayTime:=max(0,Cities .totalTime-DateSpanToMin(startTime,Cities .late)); if Cities .late0 then //consider time or not begin if Cities .early0 then //window or deadline cost:=penaltyW*fCities .waitTime +penaltyD*fCities .delayTime else cost:=penaltyD*fCities .delayTime; end else cost:=0; result:=cost; end; function TTSPController.DateSpanToMin(d1,d2:TDateTime):integer; var span:TDateTime; Year, Month, Day, Hour, Min, Sec, MSec: Word; begin span:=abs(d2-d1); DecodeDate(span, Year, Month, Day); DecodeTime(span, Hour, Min, Sec, MSec); result:=Min; end; //return the position in the vehicles array function TTSPController.GetVehicleInfo( routeInt:Tint):integer; begin result:=routeInt-fOldCityCount-1; end; function TTSPController.GetWeightConstraint( Individual: IIndividual): TInt; var Indi: ITSPIndividual; totalCapacity,maxCapacity: TFloat; i,j:TInt; tempArray:array of TInt; tempResult:TInt; begin Indi := Individual as ITSPIndividual; SetLength(tempArray, fCityCount+1); tempResult:=0; ///////////////////////////////////////////////////////// for i:=0 to fCityCount-1 do begin if Indi.RouteArray =fOldCityCount+1 then break; end; for j:=0 to fCityCount-i-1 do begin tempArray := Indi.RouteArray ; end; for j:=fCityCount-i to fCityCount-1 do begin tempArray := Indi.RouteArray ; end; tempArray := tempArray ; ////////////////////////////////////////////////////////// //totalCapacity:=fCities ].supply; //supply maxCapacity:=fVehicles )].volume; totalCapacity:=maxCapacity; for i:=0 to fCityCount do begin if (FCities ].id=fOldCityCount)and(FCities ].id0) then begin totalCapacity:=totalCapacity+FCities ].supply-FCities ].demand; if (totalCapacitymaxCapacity)or(totalCapacity0) then begin tempResult:=tempResult+1; //break; end; end; if FCities ].idfOldCityCount then begin //totalCapacity:=fCities ].supply; //supply maxCapacity:=fVehicles )].volume; totalCapacity:=maxCapacity; end; end; SetLength(tempArray,0); result:=tempResult; end; function TTSPController.GetBackConstraint( Individual: IIndividual): TInt; var Indi: ITSPIndividual; i,j:TInt; tempArray:array of TInt; tempResult:TInt; begin Indi := Individual as ITSPIndividual; SetLength(tempArray, fCityCount+1); tempResult:=0; for i:=0 to fCityCount-1 do begin if Indi.RouteArray =fOldCityCount+1 then break; end; for j:=0 to fCityCount-i-1 do begin tempArray := Indi.RouteArray ; end; for j:=fCityCount-i to fCityCount-1 do begin tempArray := Indi.RouteArray ; end; tempArray :=tempArray ; {tempArray :=11;tempArray :=5;tempArray :=8;tempArray :=7; tempArray :=9;tempArray :=6;tempArray :=12;tempArray :=10; tempArray :=2;tempArray :=4;tempArray :=3;tempArray :=1; tempArray :=0;tempArray :=11;tempArray :=3;tempArray :=1; tempArray :=4;tempArray :=11;//10,2,2} for i:=0 to fCityCount-1 do begin if (Cities ].id=fOldCityCount) then begin fCities ].serviceDepot:= fCities ].serviceDepot; end; if (Cities ].id=fOldCityCount)and(Cities ].id=1)and(Cities ].id fOldCityCount) then begin if Cities ].serviceDepotCities ].serviceDepot then //back to the start point begin tempResult:=tempResult+1; // break; end; end; end; SetLength(tempArray,0); result:=tempResult; end; function TTSPController.GetTimeConstraint( Individual: IIndividual): TInt; var Indi: ITSPIndividual; i,j:TInt; totalTimeCost:TFloat; tempArray:array of TInt; tempResult:TInt; begin Indi := Individual as ITSPIndividual; SetLength(tempArray, fCityCount+1); tempResult:=0; for i:=0 to fCityCount-1 do begin if Indi.RouteArray =fOldCityCount+1 then break; end; for j:=0 to fCityCount-i-1 do begin tempArray := Indi.RouteArray ; end; for j:=fCityCount-i to fCityCount-1 do begin tempArray := Indi.RouteArray ; end; tempArray :=tempArray ; totalTimeCost:=0; for i:=0 to fCityCount-1 do begin totalTimeCost:=totalTimeCost+timeCostBetween(tempArray ,tempArray ); end; if totalTimeCost0 then tempResult:=1; SetLength(tempArray,0); end; function TTSPController.GetCity(I: Integer): TPoint2D; begin result := fCities ; end; function TTSPController.GetNoVehicle(I: Integer): TInt; begin result := fNoVehicles ; end; function TTSPController.GetCityCount: Integer; begin result := fCityCount; end; function TTSPController.GetOldCityCount: Integer; begin result := fOldCityCount; end; function TTSPController.GetTravelCount: Integer; begin result := fTravelCount; end; function TTSPController.GetDepotCount: Integer; begin result := fDepotCount; end; function TTSPController.GetXmax: TFloat; begin result := fXmax; end; function TTSPController.GetXmin: TFloat; begin result := fXmin; end; function TTSPController.GetYmax: TFloat; begin result := fYmax; end; function TTSPController.GetYmin: TFloat; begin result := fYmin; end; procedure TTSPController.RandomCities; //from database var i,j,k,m,intTemp,totalVehicleCount: Integer; tempVehicle:TVehicle; begin ////////////////////////////////////////////////////////// fNoVehicles :=0; totalVehicleCount:=0; for i:=1 to fDepotCount do //from depots database begin fNoVehicles :=fTravelCount +1; totalVehicleCount:=totalVehicleCount+ fNoVehicles ; //real and virtual vehicles end; SetLength(fVehicles,totalVehicleCount); intTemp:=0; for i:=1 to fDepotCount do begin for j:=intTemp to intTemp+fNoVehicles -2 do begin fVehicles .index:=j+1; fVehicles .id:='real vehicle'; fVehicles .volume:=50; end; with fVehicles -1] do begin index:=intTemp+fNoVehicles ; id:='virtual vehicle'; volume:=0; end; intTemp:=intTemp+ fNoVehicles ; end; /////////////////////////////////////////////////////////// intTemp:=0; for i:=1 to fDepotCount do //depot 1--value begin intTemp:=intTemp + fNoVehicles ; end; for i := 0 to FOldCityCount do //from database begin FCities .id:= i; FCities .X := Xmin+0.5+(Xmax-Xmin-1.0)*Random; FCities .Y := Ymin+0.5+(Ymax-Ymin-1.0)*Random; FCities .early:=0; FCities .late:=0; //TDateTime FCities .serviceTime:=0; FCities .totalTime:=0; FCities .waitTime:=0; FCities .delayTime:=0; end; for i:=FOldCityCount+1 to FCityCount-1 do begin FCities .id:= i; if fDepotCount=1 then begin FCities .X := Xmin+0.5+(Xmax-Xmin-1.0)*RandomRange(2,4)/5; FCities .Y := Ymin+0.5+(Ymax-Ymin-1.0)*RandomRange(2,4)/5; end else begin FCities .X := Xmin+0.5+(Xmax-Xmin-1.0)*Random; FCities .Y := Ymin+0.5+(Ymax-Ymin-1.0)*Random; end; FCities .early:=0; FCities .late:=0; //TDateTime FCities .serviceTime:=0; FCities .totalTime:=0; FCities .waitTime:=0; FCities .delayTime:=0; end; for i := 0 to FOldCityCount do begin FCities .serviceDepot:=i; end; m:=FOldCityCount+1; for k:=1 to fDepotCount do begin for j:=0 to fNoVehicles -1 do begin FCities .serviceDepot:= fOldCityCount+k; m:=m+1; end; end; //supply and demand //////////////////////////from database FCities .demand:=0; FCities .supply:=0; for i:=1 to FOldCityCount do begin FCities .demand:=10; FCities .supply:=0; end; for i:=FOldCityCount+1 to FCityCount-1 do begin FCities .demand:=0; FCities .supply:=50; end; //////////////////////////////////////////////////////////// intTemp:=0; for i:=0 to fDepotCount-1 do begin intTemp:=intTemp+fNoVehicles ; for j:=2 to fNoVehicles do begin FCities .X :=FCities .X; FCities .Y :=FCities .Y; end; end; writeTimeArray; writeCostArray; end; procedure TTSPController.writeTimeArray; //database var i,j:integer; begin SetLength(timeArray,fCityCount,fCityCount); for i:=0 to fCityCount-1 do begin for j:=0 to fCityCount-1 do begin if i=j then timeArray :=0 else timeArray :=10; end; end; end; procedure TTSPController.writeCostArray; //database var i,j:integer; begin SetLength(costArray,fCityCount,fCityCount); for i:=0 to fCityCount-1 do begin for j:=0 to fCityCount-1 do begin if i=j then costArray :=0 else costArray :=costBetween(i,j); end; end; end; procedure TTSPController.SetCityCount(const Value: Integer); begin SetLength(fCities, Value); fCityCount := Value; RandomCities; end; procedure TTSPController.SetOldCityCount(const Value: Integer); begin fOldCityCount := Value; end; procedure TTSPController.SetTravelCount(const Value: Integer); /////////// begin fTravelCount := Value; end; procedure TTSPController.SetDepotCount(const Value: Integer); /////////// begin SetLength(fNoVehicles, Value+1); /////////////// fDepotCount := Value; end; procedure TTSPController.SetXmax(const Value: TFloat); begin fXmax := Value; end; procedure TTSPController.SetXmin(const Value: TFloat); begin fXmin := Value; end; procedure TTSPController.SetYmax(const Value: TFloat); begin fYmax := Value; end; procedure TTSPController.SetYmin(const Value: TFloat); begin fYmin := Value; end; end.
个人分类: 量子原胞自动机|500 次阅读|0 个评论
分享 欧拉回路算法
葱冲拌数学 2012-8-18 19:54
欧拉回路算法: 步骤1(起始通路) (a)令E是G的边集。 (b)选择一个顶点U,令通路C仅有U组成。 步骤2(扩充通路) while(E非空) 步骤2.1(为扩充选择一个起始点) (a)令V是C中的顶点,它与E中的某条边关联。 (b)令通路P恰好包含V。 步骤2.2(扩充P,使它成为一条从V到V的通路) (a)令W=V (b)while E中有邻接W的边e (a)从E中删除去e //从与顶点关联的边中任意选取一条边删去 (b)用邻接e的另一个顶点替换W (c)将边e和顶点W添加到通路P. endwhile 步骤2.3(扩充C) //将所有集合拼凑到一个大集合中 用通路P替换C中出现的任意的某一个V。 endwhile 步骤3(输出) 通路C是一条欧拉回路
298 次阅读|0 个评论
分享 具体些
zhangquan 2012-8-2 21:37
一、蒙特卡罗算法 1946年,美国拉斯阿莫斯国家实验室的三位科学家John von Neumann,Stan Ulam 和 Nick Metropolis 共同发明了,蒙特卡罗方法。 此算法被评为20世纪最伟大的十大算法之一,详情,请参见我的博文: http://blog.csdn.net/v_JULY_v/archive/2011/01/10/6127953.aspx 蒙特卡罗方法(Monte Carlo method),又称随机抽样或统计模拟方法,是一种以概率统计理论为指导 的一类非常重要的数值计算方法。此方法使用随机数(或更常见的伪随机数)来解决很多计算问题的方 法。 由于传统的经验方法由于不能逼近真实的物理过程,很难得到满意的结果,而蒙特卡罗方法由于能够真 实地模拟实际物理过程,故解决问题与实际非常符合,可以得到很圆满的结果。 蒙特卡罗方法的基本原理及思想如下: 当所求解问题是某种随机事件出现的概率,或者是某个随机变量的期望值时,通过某种“实验”的方法 ,以这种事件出现的频率估计这一随机事件的概率,或者得到这个随机变量的某些数字特征,并将其作 为问题的解。 有一个例子可以使你比较直观地了解蒙特卡洛方法: 假设我们要计算一个不规则图形的面积,那么图形的不规则程度和分析性计算(比如,积分)的复杂程 度是成正比的。蒙特卡洛方法是怎么计算的呢?假想你有一袋豆子,把豆子均匀地朝这个图形上撒,然 后数这个图形之中有多少颗豆子,这个豆子的数目就是图形的面积。当你的豆子越小,撒的越多的时候 ,结果就越精确。 在这里我们要假定豆子都在一个平面上,相互之间没有重叠。 蒙特卡罗方法通过抓住事物运动的几何数量和几何特征,利用数学方法来加以模拟,即进行一种数字模 拟实验。它是以一个概率模型为基础,按照这个模型所描绘的过程,通过模拟实验的结果,作为问题的 近似解。 蒙特卡罗方法与一般计算方法有很大区别,一般计算方法对于解决多维或因素复杂的问题非常困难,而 蒙特卡罗方法对于解决这方面的问题却比较简单。其特点如下: I、 直接追踪粒子,物理思路清晰,易于理解。 II、 采用随机抽样的方法,较真切的模拟粒子输运的过程,反映了统计涨落的规律。 III、不受系统多维、多因素等复杂性的限制,是解决复杂系统粒子输运问题的好方法。 等等。 此算法,日后还会在本BLOG 内详细阐述。 二、数据拟合、参数估计、插值等数据处理算法 我们通常会遇到大量的数据需要处理, 而处理数据的关键就在于这些算法,通常使用Matlab作为工具。 数据拟合在数学建模比赛中中有应用,与图形处理有关的问题很多与拟合有关系,一个例子就是98年数 学建模美国赛A题,生物组织切片的三维插值处理,94年A题逢山开路,山体海拔高度的插值计算,还有 吵的沸沸扬扬可能会考的“非典”问题也要用到数据拟合算法,观察数据的走向进行处理。 此类问题在 MATLAB 中有很多现成的函数可以调用,熟悉MATLAB,这些方法都能游刃有余的用好。 三、线性规划、整数规划、多元规划、二次规划等规划类问题 数学建模竞赛中很多问题都和数学规划有关,可以说不少的模型都可以归结为一组不等式作为约束条件 、几个函数表达式作为目标函数的问题,遇到这类问题,求解就是关键了,比如98年B题,用很多不等式 完全可以把问题刻画清楚,因此列举出规划后用 Lindo 、 Lingo 等软件来进行解决比较方便,所以还 需要熟悉这两个软件。 四、图论算法 这类问题算法有很多, 包括: Dijkstra 、 Floyd 、 Prim 、 Bellman-Ford ,最大流,二分匹配等问题。 关于此类图论算法,可参考Introduction to Algorithms--算法导论,关于图算法的第22章-第26章。 同时,本BLOG内经典算法研究系列,对Dijkstra算法有所简单描述, ----------- 经典算法研究系列:二、Dijkstra 算法初探 http://blog.csdn.net/v_JULY_v/archive/2010/12/24/6096981.aspx 更多,请关注本BLOG 日后更新的博文。 五、动态规划、回溯搜索、分治算法、分支定界等计算机算法 在数学建模竞赛中,如:92 年B题用分枝定界法, 97年B题是典型的动态规划问题, 此外 98 年 B 题体现了分治算法。 这方面问题和 ACM 程序设计竞赛中的问题类似, 推荐看一下算法导论,与《计算机算法设计与分析》(电子工业出版社)等与计算机算法有关的书。 六、最优化理论的三大经典算法:模拟退火法、神经网络、遗传算法 这十几年来最优化理论有了飞速发展,模拟退火法、神经网络、遗传算法这三类算法发展很快。 在数学建模竞赛中:比如97年A题的模拟退火算法,00年B题的神经网络分类算法,01年B题这种难题也可 以使用神经网络,还有美国竞赛89年A题也和 BP 算法有关系,当时是86年刚提出BP算法,89年就考了, 说明赛题可能是当今前沿科技的抽象体现。 03 年 B 题伽马刀问题也是目前研究的课题,目前算法最佳的是遗传算法。 另,本人对人工智能非常感兴趣,遗传算法已在本BLOG内有所阐述,敬请参见。 ---------- 经典算法研究系列:七、深入浅出遗传算法,透析GA本质 http://blog.csdn.net/v_JULY_v/archive/2011/01/12/6132775.aspx 其它俩大算法,模拟退火法,与神经网络,也定会在本BLOG内日后的博文更新中,详细阐述。 七、网格算法和穷举法 网格算法和穷举法一样,只是网格法是连续问题的穷举。 比如要求在 N 个变量情况下的最优化问题,那么对这些变量可取的空间进行采点, 比如在 区间内取 M +1 个点,就是 a; a +( b ? a ) =M; a +2 ¢ ( b ? a ) =M ; …;b 那么这样循环就需要进行 ( M + 1) N 次运算,所以计算量很大。 在数学建模竞赛中:比如 97 年 A 题、 99 年 B 题都可以用网格法搜索,这种方法最好在运算速度较 快的计算机中进行,还有要用高级语言来做,最好不要用 MATLAB 做网格,否则会算很久。 穷举法大家都熟悉,自不用多说了。 八、一些连续离散化方法 大部分物理问题的编程解决,都和这种方法有一定的联系。物理问题是反映我们生活在一个连续的世界 中,计算机只能处理离散的量,所以需要对连续量进行离散处理。 这种方法应用很广,而且和上面的很多算法有关。 事实上,网格算法、蒙特卡罗算法、模拟退火都用了这个思想。 九、数值分析算法 数值分析(numerical analysis),是数学的一个分支,主要研究连续数学(区别于离散数学)问题的 算法。 如果在比赛中采用高级语言进行编程的话,那一些数值分析中常用的算法比 如方程组求解、矩阵运算、 函数积分等算法就需要额外编写库函数进行调用。 这类算法是针对高级语言而专门设的,如果你用的是 MATLAB 、 Mathematica ,大可不必准备, 因为像数值分析中有很多函数一般的数学软件是具备的。 十、图象处理算法 在数学建模竞赛中:比如01 年 A 题中需要你会读 BMP 图象、美国赛 98 年 A 题需要你知道三维插值 计算, 03 年 B 题要求更高,不但需要编程计算还要进行处理,而数模论文中也有很多图片需要展示, 因此图象处理就是关键。做好这类问题,重要的是把MATLAB 学好,特别是图象处理的部分。
319 次阅读|0 个评论
分享 数学建模算法
遇到 2012-7-21 20:50
可以找我要
116 次阅读|0 个评论
分享 数模人的骄傲啊,遗传算法用在中国无人机计划上
sdccumcm 2012-7-9 10:56
英国新科学家网站报道,中国海军研究人员日前透露了他们如何利用舰射无人机猎杀潜艇的计划。该计划由中国海军大连舰艇学院制定,将利用遗传算法选择最佳无人机猎杀模式。遗传算法是一种通过模拟自然进化过程搜索最优解的方法。 通过该算法打造无人机,最终可达到使无人机充分利用燃料、应对空中与海上威胁以及联合水下声纳浮标共同作业的目的。 这种无人机的出现,或许能够在有关台湾冲突问题中派上用场。问题是,中国为何公布相关信息,特别是在当前其敌对方正发展潜艇反制措施之时? 众所周知,潜艇猎杀相关情报属机密信息。例如,二战时期,英国布莱切利公园成功破解了纳粹德国的 Enigma 密码,并利用破解情报成功定位德军潜艇。 这是大西洋战役获同盟国阵营赢得大西洋战役的关键,但英国在长时间内却对相关信息秘而不宣。 因此,中国公布该计划多少有些令人困惑。不过,这种情况也不是第一次出现。 2010 年,中国研究人员就曾发布过一篇有关如何大范围攻击、破坏美国电网的文章,令美方颇为愤怒。 附:遗传算法 (Genetic Algorithms) 简介 遗传算法( Genetic Algorithms )是基于生物进化理论的原理发展起来的一种广为应用的、高效的随机搜索与优化的方法。其主要特点是群体搜索策略和群体中个体之间的信息交换,搜索不依赖于梯度信息。 它是在 70 年代初期由美国密执根( Michigan )大学的霍兰( Holland )教授发展起来的。 1975 年霍兰教授发表了第一本比较系统论述遗传算法的专著《自然系统与人工系统中的适应性》(《 Adaptation in Natural and Artificial Systems 》)。 遗传算法最初被研究的出发点不是为专门解决最优化问题而设计的,它与进化策略、进化规划共同构成了进化算法的主要框架,都是为当时人工智能的发展服务的。迄今为止,遗传算法是进化算法中最广为人知的算法。 近几年来,遗传算法主要在复杂优化问题求解和工业工程领域应用,取得了一些令人信服的结果,所以引起了很多人的关注,而且在发展过程中,进化策略、进化规划和遗传算法之间差异越来越小。 遗传算法成功的应用包括:作业调度与排序、可靠性设计、车辆路径选择与调度、成组技术、设备布置与分配、交通问题等等。
183 次阅读|0 个评论
分享 数学建模中常用的方法与算法集锦(个人整理)
热度 1 liuzi0705 2012-5-7 11:04
在数学建模中常用的方法: 类比法、 二分法、 量纲分析法、 差分法、 变分法、 图论法、 层次分析法、 数据拟合法、 回归分析法、 数学规划(线性规划,非线性规划,整数规划,动态规划,目标规划)、 机理分析、 排队方法、 对策方法、 决策方法、 模糊评判方法、 时间序列方法、 灰色理论方法、 现代优化算法(禁忌搜索算法,模拟退火算法,遗传算法,神经网络)。 用这些方法可以解下列一些模型: 优化模型、 微分方程模型、 统计模型、 概率模型、 图论模型、 决策模型。 在数学建模中常用的算法: 1 :蒙特卡罗算法; 2 :数据拟合、参数估计、插值等数据处理算法(常用 matlab 实现); 3 :线性规划、整数规划、多元规划、二次规划 ( 用 lingo 、 lingdo 、 matlab 即可实现 ) ; 4 :图论算法(包括最短路、网络流、二分图); 5 :动态规划、回溯搜索、分治算法、分支界定; 6 :最优化理论的三大经典算法(模拟退火算法、神经网络算法、遗传算法); 7 :网格算法和穷举法; 8 :连续数据离散化; 9 :数值分析算法; 10 :图象处理算法(常用 matlab 来实现)。 拟合与插值方法 (给出一批数据点,确定满足特定要求的曲线或者曲面,从而反映对象整体的变化趋势): matlab 可以实现一元函数,包括多项式和非线性函数的拟合以及多元函数的拟合,即回归分析,从而确定函数; 同时也可以用 matlab 实现分段线性、多项式、样条以及多维插值。 在 优化方法 中,决策变量、目标函数(尽量简单、光滑)、约束条件、求解方法是四个关键因素。其中包括无约束规则(用 fminserch 、 fminbnd 实现)线性规则(用 linprog 实现)非线性规则、( 用 fmincon 实现)多目标规划(有目标加权、效用函数)动态规划(倒向和正向)整数规划。 回归分析 :对具有相关关系的现象,根据其关系形态,选择一个合适的数学模型,用来近似地表示变量间的平均变化关系的一种统计方法 (一元线性回归、多元线性回归、非线性回归),回归分析在一组数据的基础上研究这样几个问题:建立因变量与自变量之间的回归模型(经验公式);对回归模型的可信度进行检验;判断每个自变量对因变量的影响是否显著;判断回归模型是否适合这组数据;利用回归模型对进行预报或控制。相对应的有 线性回归、多元二项式回归、非线性回归。 逐步回归分析: 从一个自变量开始,视自变量作用的显著程度,从大到地依次逐个引入回归方程:当引入的自变量由于后面变量的引入而变得不显著时,要将其剔除掉;引入一个自变量或从回归方程中剔除一个自变量,为逐步回归的一步;对于每一步都要进行值检验,以确保每次引入新的显著性变量前回归方程中只包含对作用显著的变量;这个过程反复进行,直至既无不显著的变量从回归方程中剔除,又无显著变量可引入回归方程时为止。(主要用 SAS 来实现,也可以用 matlab 软件来实现)。 聚类分析 :所研究的样本或者变量之间存在程度不同的相似性,要求设法找出一些能够度量它们之间相似程度的统计量作为分类的依据,再利用这些量将样本或者变量进行分类。 系统聚类分析 — 将 n 个样本或者 n 个指标看成 n 类,一类包括一个样本或者指标,然后将性质最接近的两类合并成为一个新类,依此类推。最终可以按照需要来决定分多少类,每类有多少样本(指标)。 系统聚类方法步骤: 计算 n 个样本两两之间的距离 构成 n 个类,每类只包含一个样品 合并距离最近的两类为一个新类 计算新类与当前各类的距离(新类与当前类的距离等于当前类与组合类中包含的类的距离最小值),若类的个数等于 1 ,转 5 ,否则转 3 画聚类图 决定类的个数和类。 判别分析 :在已知研究对象分成若干类型,并已取得各种类型的一批已知样品的观测数据,在此基础上根据某些准则建立判别式,然后对未知类型的样品进行判别分类。 距离判别法 — 首先根据已知分类的数据,分别计算各类的重心,计算新个体到每类的距离,确定最短的距离(欧氏距离、马氏距离) Fisher 判别法 — 利用已知类别个体的指标构造判别式(同类差别较小、不同类差别较大),按照判别式的值判断新个体的类别 Bayes 判别法 — 计算新给样品属于各总体的条件概率,比较概率的大小,然后将新样品判归为来自概率最大的总体 模糊数学 :研究和处理模糊性现象的数学 (概念与其对立面之间没有一条明确的分界线)与模糊数学相关的问题:模糊分类问题 — 已知若干个相互之间不分明的模糊概念,需要判断某个确定事物用哪一个模糊概念来反映更合理准确;模糊相似选择 — 按某种性质对一组事物或对象排序是一类常见的问题,但是用来比较的性质具有边界不分明的模糊性;模糊聚类分析 — 根据研究对象本身的属性构造模糊矩阵,在此基础上根据一定的隶属度来确定其分类关系 ;模糊层次分析法 — 两两比较指标的确定;模糊综合评判 — 综合评判就是对受到多个因素制约的事物或对象作出一个总的评价,如产品质量评定、科技成果鉴定、某种作物种植适应性的评价等,都属于综合评判问题。由于从多方面对事物进行评价难免带有模糊性和主观性,采用模糊数学的方法进行综合评判将使结果尽量客观从而取得更好的实际效果 。 时间序列 是按时间顺序排列的、随时间变化且相互关联的数据序列 — 通过对预测目标自身时间序列的处理,来研究其变化趋势(长期趋势变动、季节变动、循环变动、不规则变动) 自回归模型: 一般自回归模型 AR(n)— 系统在时刻 t 的响应 X(t) 仅与其以前时刻的响应 X(t-1),… , X(t-n) 有关,而与其以前时刻进入系统的扰动无关 ;移动平均模型 MA(m)— 系统在时刻 t 的响应 X(t) ,与其以前任何时刻的响应无关,而与其以前时刻进入系统的扰动 a(t-1),…,a(t-m) 存在着一定的相关关系 ;自回归移动平均模型 ARMA(n,m)— 系统在时刻 t 的响应 X(t) ,不仅与其前 n 个时刻的自身值有关,而且还与其前 m 个时刻进入系统的扰动存在一定的依存关系 。 时间序列建模的基本步骤 数据的预处理:数据的剔取及提取趋势项 取 n=1 ,拟合 ARMA(2n,2n-1) (即 ARMA(2,1) )模型 n=n+1 ,拟合 ARMA(2n,2n-1) 模型 用 F 准则检验模型的适用性。若检验显著,则转入第 2 步。若检验不显著,转入第 5 步。 检查远端时刻的系数值的值是否很小,其置信区间是否包含零。若不是,则适用的模型就是 ARMA(2n,2n-1) 。若很小,且其置信区间包含零,则拟合 ARMA(2n-1,2n-2) 。 利用 F 准则检验模型 ARMA(2n,2n-1) 和 ARMA(2n-1,2n-2 ) ,若 F 值不显著,转入第 7 步;若 F 值显著,转入第 8 步。 舍弃小的 MA 参数,拟合 m2n-2 的模型 ARMA(2n-1,m) ,并用 F 准则进行检验。重复这一过程,直到得出具有最小参数的适用模型为止 舍弃小的 MA 参数,拟合 m2n-1 的模型 ARMA(2n,m) ,并用 F 准则进行检验。重复这一过程,直到得出具有最小参数的适用模型为止。 图论方法: 最短路问题:两个指定顶点之间的最短路径 — 给出了一个连接若干个城镇的铁路网络,在这个网络的两个指定城镇间,找一条最短铁路线 ( Dijkstra 算法 )每对顶点之间的最短路径 ( Dijkstra 算法、 Floyd 算法 )。 最小生成树问题:连线问题 — 欲修筑连接多个城市的铁路设计一个线路图,使总造价最低( prim 算法、 Kruskal 算法 )。 图的匹配问题:人员分派问题: n 个工作人员去做件 n 份工作,每人适合做其中一件或几件,问能否每人都有一份适合的工作?如果不能,最多几人可以有适合的工作? ( 匈牙利算法 ) 。 遍历性问题:中国邮递员问题 — 邮递员发送邮件时,要从邮局出发,经过他投递范围内的每条街道至少一次,然后返回邮局,但邮递员希望选择一条行程最短的路线 最大流问题。 运输问题: 最小费用最大流问题:在运输问题中,人们总是希望在完成运输任务的同时,寻求一个使总的运输费用最小的运输方案
个人分类: 数模经验|884 次阅读|1 个评论
分享 关于图论最短路算法的新设想
一枕清霜 2012-4-10 16:34
其实我看floyd算法是一种代试探性的穷举法,是以一种逐步的试探进行的,但我认为,其实光学基于折射定律,费马原理的光的折射模型,也是一种求最短路径,最少时间的算法,举个例子,正如我之前提到过的,最速降线涉及的时间最少的路径求解,是可以装换用折射模型求解的,只不过数学上的处理比较复杂,以为它涉及一个由分立到连续的转变 其实我认为,一个能用图论求解的问题,可以转化为以场论的方式求解,既是由分立转化为连续,整个模型的主要 问题也是将实际问题场化的方法探究,模型的求解是去解一个大的微分方程,最后的解的路线应该是连续曲线,还要把它优化成直线, 我认为这个模型的应用不会广泛,因为他的结果必与floyd算法是相同的,但是数学的难度上,所需的知识上,都要多余floyd算法,但他的求解却是一种确定性的,既是不需要以试探性性的方法进行,他可以给出新的路径,也就说通过他的结果,可以新建出更好的路径,而floyd算法基于的路径已经给出,折射模型会给出更优的路径,在折射模型中,路径可以是不先给出的! 仅是设想,供数学,物理爱好者探讨
342 次阅读|0 个评论
分享 数学建模十大算法
Uni—— 2012-3-26 09:43
数学建模十大算法机器代码实现,下载地址: http://download.csdn.net/detail/wanglihongwm/4152168
331 次阅读|0 个评论
分享 出道
guoqiangi1 2012-3-20 21:54
很高兴今天整整的认识了数学中国这个平台,在这里我可以学习我自己的想学的算法,happy,happy.....
个人分类: 杂谈|270 次阅读|0 个评论
分享 [转帖]TF/IDF算法
CQU_LUcifeR 2012-2-12 22:47
一直在看代码,发现切词和相关性中都用到了tf和idf,百度了一下,发现了一篇文章挺好的。转载一下。 TF/IDF TF/IDF(term frequency/inverse document frequency) 的概念被公认为信息检索中最重要的发明。 一。TF/IDF描述单个term与特定document的相关性 TF(Term Frequency): 表示一个term与某个document的相关性。 公式为这个term在document中出现的次数除以该document中所有term出现的总次数. IDF(Inverse Document Frequency)表示一个term表示document的主题的权重大小。主要是通过包含了该term的docuement的数量和docuement set的总数量来比较的。出现的次数越多,权重越小。 公式是log(D/Dt) D是docuemnt set的总数量, Dt是包含了该term的document的总数。 这样,根据关键字k1,k2,k3进行搜索结果的相关性就变成TF1*IDF1 + TF2*IDF2 + TF3*IDF3。比如document1的term总量为1000,k1,k2,k3在document1出现的次数是100,200,50。包含了k1, k2, k3的docuement总量分别是 1000, 10000,5000。document set的总量为10000。 TF1 = 100/1000 = 0.1 TF2 = 200/1000 = 0.2 TF3 = 50/1000 = 0.05 IDF1 = log(10000/1000) = log(10) = 2.3 IDF2 = log(10000/100000) = log(1) = 0; IDF3 = log(10000/5000) = log(2) = 0.69 这样关键字k1,k2,k3与docuement1的相关性= 0.1*2.3 + 0.2*0 + 0.05*0.69 = 0.2645 其中k1比k3的比重在document1要大,k2的比重是0. TF/IDF 的概念就是一个特定条件下、关键词的概率分布的交叉熵(Kullback-Leibler Divergence). 二。用TF/IDF来描述document的相似性。 假如document1和document2的term的TF/IDF分别是t11,t12,t13,...t1n和t21,t22,t23,...,t2n.他们之间的相似性可以用余弦定理来表示。则: cos(d1,d2) = d1和d2的内积/(d1的长度*d2的长度) = (t11*t21 + t12*t22 + t13*t23 + ... + t1n*t2n)/(|d1|*|d2|). d1 = sqrt(t11*t11 + t12*t12 + t13*t13 + ... + t1n*t1n); 夹角越大,相似性越大。为1则表示d1和d2一致。 在今日我们可以从网络上吸收大量资讯,有时候一堆文章看不完。如果我们想要吸收资讯,时间却又不够的时候,使用电脑帮我们过滤资讯,或是用电脑帮我们做个总整理,是个方法。如果今天手中有一篇文章,我们想要用电脑帮我们找出这篇文章最重要的关键字,要怎麽做呢?在资讯检索 (IR: Information Retrieval)领域里面,有个基础的方法,入门必学的方法,就是使用 TF 和 IDF (TF: Term Frequency, IDF: Inverse Document Frequency)。使用这两个估计值,可以让电脑具有计算重要关键字的能力,进而节省我们的时间。   接下来让我们看看,TF 和 IDF 个是甚麽东西呢?TF 全名是Term Frequency,也就是某个关键字出现的次数,譬如说某篇文章里面,「电脑」这个词出现很多次,或是「使用者需求」这个词出现很多次,那麽这些词句的出现频率,就会很高。一篇文章中出现很多次的词句,必定有其重要性。譬如说一篇论述「人工智慧」的文章,「人工智慧」这个词句再文章中出现的频率也一定很高。然而为甚麽除了 TF (Term Frequency) 以外,还要有 IDF (Inverse Document Frequency) 呢?   让我们先想想,如果单使用某个字词出现的频率,来判断一篇文章最重要的关键字,会有甚麽困难。首先,我们会遇到一些常用字词,出现的频率也很高,会和重要字词出现的频率一样高,让电脑因此无法分辨出,哪些是常用字词,那些是重要字词。如果就英文来说,有个规则是语言学家 (linguist) 归纳出来的规则,叫做 Zipf’s Law 引述中文维基百科的一段介绍如下:   从根本上讲, 齐夫定律 可以表述为, 在 自然语言 的 语素库 里, 一个单词出现的频率与它在频率表里的排名成 反比. 所以, 频率最高的单词出现的频率大约是出现频率第二位的单词的 2 倍,而出现频率第二位的单词则是出现频率第四位的单词的2倍。这个定律被作为任何与 power law probability distributions 有关的事物的参考。 这个 “定律” 是 Harvard linguist George Kingsley Zipf (IPA )发表的。 比如, 在 Brown 语库 , “the” 是最常见的单词,它在这个语库中出现了大约 7 %(10 万单词中出现 69971 次)。正如齐夫定律中所描述的一样,出现次数为第二位的单词 “of” 占了整个语库中的 3.5% (36411次), 之後的是”and” (28852次). 仅仅 135 但此项就占了 Brown 语库的一半。   所以我们现在知道问题在哪边了。如果只用词句出现的频率来判断某一篇文章里面最重要的关键字,我们可能会找到常用字,而不是最重要的字,像是英文里面的 “the”、”a”、”it”,都是常常出现的字,但是通常一篇文章里面最重要的字不是这些字,即使那些重要的字出现的频率也很高。   这个时候我们要怎麽办呢?IDF 在这个时候就帮上忙了。在了解 IDF 之前,我们先了解 DF 是甚麽。DF 就是Document Frequency,也就是说,如果今天我们手中有固定 N 篇文章,某个关键字的 Document Frquency (DF),就是说这个关键字在 N 篇文章里面出现了几次。Inverse Document Frequency (IDF) 则是把 DF 取倒数,如此一来,一个数字乘以 IDF,就等於是除以 DF 的意思。   有了 TF 和 IDF 以後,我们就可以计算 TF 乘上 IDF,对每一个关键字都算出一个分数。这个分数的高低,就代表了这个关键字在某篇文章中的重要程度。为甚麽我们说这样子可以找出重要的字,而不是常出现的字呢?因为 TF 会把某篇文章中,出现最多次的排在第一位,其次的排在第二位,以此类推。然而乘上 IDF 以後,也就是除以 DF,那些常常出现的字,像是英文中的 “the”、”a”、”it”,因为每一篇文章都会出现,所以 DF 就大。DF 大,取倒数之後的 IDF 就小,IDF 小,乘上 TF 以後,虽然”the”、”a”、”it”在某篇文章中出现的频率很高,但是因为 IDF 小,TF * IDF 一相乘,重要性就变低了,我们 (电脑程式) 就不会把这些常出现的字,误认为是重要的字了!   真正重要的字会得到甚麽样子的分数呢?如果这篇文章刚好在讲 AI,”AI” 出现很多次,因此 “AI” 在这篇文章里面的 TF 很高。然而我们电脑资料库里面的 N 篇文章,并不是每一篇都在讲 AI,也因此”AI”可能只有在 N 篇文章里面的某 3 篇文章出现,因此 DF 只有 3,IDF 变成 0.33,假设我们 N = 100 有 100 篇文章在资料库里面,其他常出现字像是 “the” 每一篇都出现,DF 就是 100,IDF 就是 0.01。所以 “AI” 的 IDF 会比 “the” 的 IDF 高,假设这篇文章中 “AI” 和 “the” 两个字出现的次数刚好一样,乘上 IDF 以後,”AI” 这个字的分数就比 “the” 这个字的分数来的高,电脑也就会判断 “AI” 是这篇文章重要的关键字,而 “the” 这个字并不是这篇文章的重要关键字。   所以经由 TF * IDF,我们可以计算某个关键字,在某篇文章里面的重要性。从这一个方向,我们可以计算一篇文章中重点的字有哪些,帮我们做一篇文章的总整理。从相反的方向,我们可以给定关键字,然後再每一篇文章里面为这个关键字计算一次 TF * IDF,然後比较哪一篇文章,这个关键字是最具重要性的,用这个方法找出和一个关键字最相关的文章。不管是从文章找出重点字词,或是由关键字找相关文章,TF * IDF 都是个基本且不错的方法。会写程式又还没?试过这个方法的读者,或许可以亲自试试看,不过可能要先自己准备文章资料库 (corpus),或是从网际网路上面用网页撷取器 (crawler) 存几篇有兴趣的网页,然後把 HTML 标签清理乾净,剩下纯文字,就可以用这个方法来小试身手罗!   我们也可以比较一下人类和电脑的不同。电脑做数学数字的计算,或是执行固定的步骤 ,非常擅长,速度也很快。人类可以了解一个字的意思,读完一篇文章以後,了解了意思,之後要找这篇文章最重要的关键字,是从「意义」开始,回忆出或做出结论,这篇文章重要的关键字是甚麽。   然而如果要电脑也遵照这个方向,先了解字的意义,再了解文章的意义,然後在做出结论,这篇文章的重要关键字,反而困难,因为要了解字的意义,电脑需要先有一个语意网路 (Semantic Network),或是知识的分类关系树 (Ontology),把字句依照语意分门别类,有如生物里面的「界门纲目科属种」一般的关系分类,才有办法了解一个字和其他字的关系。之後要了解一篇文章,又必须要了解一个句子,牵涉到自然语言处理 (NLP: Natural language Processing) 的问题,像是从句子里面找出主词、动词、和受词,以及补语,分辨出子句和主句,代名词的指称,以及前後文判断产生不同的剖析 (parsing)。了解完一句,才能了解整篇文章。 因此,TF * IDF 对於电脑来说,计算速度快,工程也不浩大,不用大型计算机就可以计算。这边也可以顺便提到 strong AI 和 weak AI 的关系。如果就工程的角度,TF * IDF 是个好方法,it works! 节省我们的时间,或是解决大问题中的一个小环节。然而 strong AI 在这边会提出「中文房间」(Chinese Room) 的论证,也就是说,电脑能够找出重要关键字,是否就代表电脑真的「知道」(understand) 关键字的意义呢?   中文房间 (Chinese Room) 简单地说,就是一个人关在房间里面,只留两个窗口,一个地方会送纸条出来,另一个地方会送纸条出去。房间里面有一本手册,里面写满对照表,记载者看到甚麽英文字,就应该输出甚麽中文字,以及一些指令的对照,譬如说窗口送一个指令说 COMBINE,就把两个中文字写在一起才送出去。接着我们在外面就开始送英文句子进去这个房间,另一个窗口就会有这句话的中文翻译跑出来。然而这个论证想要坦讨的就是,虽然这个房间看起来像是会把英文翻译成中文,但是在房间里面的那个操作人员并不懂中文,他指是按照指令,还有手册里面的对照表,机械式地动作,可是外面看起来像是这个房间会英翻中,因此这个房间应该懂得中文才对。   在这边我的看法是,也许就近程来看,我们只要有可以解决问题的解答就可以,不管电脑是否真的懂 (understand) 字的意义。然而长期来说,如果我们真的需要具有人类的智力的电脑出现,能够真的懂而不是行为上看起来懂,那麽就要仔细探讨中文房间这种论证。也许生物的方法,像是计算神经科学的方法,是一个方向。   我们可能又会问,神经元只有动作电位和静止两个状态,怎麽能了解意义?但是只有一个神经元,或许没办法了解意义,全部大脑的神经元交互作用,意义可能就因此被了解了!其中的奥妙,就是计算神经科学?试要解答的问题之一。有兴趣的读者也可以一起从人脑开始,解决 strong AI 的问题。或是有数学的高手,也许某一个数学理论,可以很漂亮地解决意义了解的问题也说不定,像是 manifolds,具有一个集合使用不同面向来观看的特性,同时具有 Global 和 Local 的性质,是个不错的候选选项。从这个方向去解决 strong AI 也是另一个可能性。总之,继续努力研究就是了! 来自: http://hi.baidu.com/hiei1125/blog/item/af9afeed9e86e34579f0552a.html
个人分类: 算法|267 次阅读|0 个评论
分享 大白话 模拟退火
hundunwei 2012-2-8 15:55
大白话解析模拟退火算法 Posted on 2010-12-20 17:01 heaad 阅读(10185) 评论(23) 编辑 收藏    优化算法入门系列文章目录(更新中):   1. 模拟退火算法   2. 遗传算法 一. 爬山算法 ( Hill Climbing ) 介绍模拟退火前,先介绍爬山算法。爬山算法是一种简单的贪心搜索算法,该算法每次从当前解的临近解空间中选择一个最优解作为当前解,直到达到一个局部最优解。 爬山算法实现很简单,其主要缺点是会陷入局部最优解,而不一定能搜索到全局最优解。如图1所示:假设C点为当前解,爬山算法搜索到A点这个局部最优解就会停止搜索,因为在A点无论向那个方向小幅度移动都不能得到更优的解。 图1 二. 模拟退火(SA,Simulated Annealing)思想 爬山法是完完全全的贪心法,每次都鼠目寸光的选择一个当前最优解,因此只能搜索到局部的最优值。模拟退火其实也是一种贪心算法,但是它的搜索过程引入了随机因素。模拟退火算法 以一定的概率 来接受一个比当前解要差的解,因此 有可能 会跳出这个局部的最优解,达到全局的最优解。以图1为例,模拟退火算法在搜索到局部最优解A后,会 以一定的概率 接受到E的移动。也许经过几次这样的不是局部最优的移动后会到达D点,于是就跳出了局部最大值A。 模拟退火算法描述: 若J( Y(i+1) )= J( Y(i) ) (即移动后得到更优解),则总是接受该移动 若J( Y(i+1) ) J( Y(i) ) (即移动后的解比当前解要差),则 以一定的概率接受移动,而且这个概率随着时间推移逐渐降低(逐渐降低才能趋向稳定)   这里的“一定的概率”的计算参考了金属冶炼的退火过程,这也是模拟退火算法名称的由来。   根据热力学的原理,在温度为T时,出现能量差为dE的降温的概率为P(dE),表示为:     P(dE) = exp( dE/(kT) )   其中k是一个常数,exp表示自然指数,且dE0。这条公式说白了就是:温度越高,出现一次能量差为dE的降温的概率就越大;温度越低,则出现降温的概率就越小。又由于dE总是小于0(否则就不叫退火了),因此dE/kT 0 ,所以P(dE)的函数取值范围是(0,1) 。   随着温度T的降低,P(dE)会逐渐降低。   我们将一次向较差解的移动看做一次温度跳变过程,我们以概率P(dE)来接受这样的移动。   关于爬山算法与模拟退火,有一个有趣的比喻:   爬山算法:兔子朝着比现在高的地方跳去。它找到了不远处的最高山峰。但是这座山不一定是珠穆朗玛峰。这就是爬山算法,它不能保证局部最优值就是全局最优值。   模拟退火:兔子喝醉了。它随机地跳了很长时间。这期间,它可能走向高处,也可能踏入平地。但是,它渐渐清醒了并朝最高方向跳去。这就是模拟退火。 下面给出模拟退火的伪代码表示。 三. 模拟退火算法伪代码 代码 !-- Code highlighting produced by Actipro CodeHighlighter (freeware) http://www.CodeHighlighter.com/ -- /* * J(y):在状态y时的评价函数值 * Y(i):表示当前状态 * Y(i+1):表示新的状态 * r: 用于控制降温的快慢 * T: 系统的温度,系统初始应该要处于一个高温的状态 * T_min :温度的下限,若温度T达到T_min,则停止搜索 */ while ( T T_min ) {   dE = J( Y(i + 1 ) ) - J( Y(i) ) ;   if ( dE = 0 ) // 表达移动后得到更优解,则总是接受移动 Y(i + 1 ) = Y(i) ; // 接受从Y(i)到Y(i+1)的移动   else   { // 函数exp( dE/T )的取值范围是(0,1) ,dE/T越大,则exp( dE/T )也 if ( exp( dE / T ) random( 0 , 1 ) ) Y(i + 1 ) = Y(i) ; // 接受从Y(i)到Y(i+1)的移动   }   T = r * T ; // 降温退火 ,0r1 。r越大,降温越慢;r越小,降温越快   /*   * 若r过大,则搜索到全局最优解的可能会较高,但搜索的过程也就较长。若r过小,则搜索的过程会很快,但最终可能会达到一个局部最优值   */   i ++ ; } 复制代码 四. 使用模拟退火算法解决旅行商问题   旅行商问题 ( TSP , Traveling Salesman Problem ) :有N个城市,要求从其中某个问题出发,唯一遍历所有城市,再回到出发的城市,求最短的路线。   旅行商问题属于所谓的NP完全问题,精确的解决TSP只能通过穷举所有的路径组合,其时间复杂度是O(N!) 。   使用模拟退火算法可以比较快的求出TSP的一条近似最优路径。(使用遗传算法也是可以的,我将在下一篇文章中介绍)模拟退火解决TSP的思路: 1. 产生一条新的遍历路径P(i+1),计算路径P(i+1)的长度L( P(i+1) ) 2. 若L(P(i+1)) L(P(i)),则接受P(i+1)为新的路径,否则以模拟退火的那个概率接受P(i+1) ,然后降温 3. 重复步骤1,2直到满足退出条件   产生新的遍历路径的方法有很多,下面列举其中3种: 1. 随机选择2个节点,交换路径中的这2个节点的顺序。 2. 随机选择2个节点,将路径中这2个节点间的节点顺序逆转。 3. 随机选择3个节点m,n,k,然后将节点m与n间的节点移位到节点k后面。 五. 算法评价 模拟退火算法是一种随机算法,并不一定能找到全局的最优解,可以比较快的找到问题的近似最优解。如果参数设置得当,模拟退火算法搜索效率比穷举法要高。  from here : http://www.cnblogs.com/heaad/ 转载请注明
254 次阅读|0 个评论
qq
收缩
  • 电话咨询

  • 04714969085

关于我们| 联系我们| 诚征英才| 对外合作| 产品服务| QQ

手机版|Archiver| |繁體中文 手机客户端  

蒙公网安备 15010502000194号

Powered by Discuz! X2.5   © 2001-2013 数学建模网-数学中国 ( 蒙ICP备14002410号-3 蒙BBS备-0002号 )     论坛法律顾问:王兆丰

GMT+8, 2024-5-24 13:55 , Processed in 0.900045 second(s), 36 queries .

回顶部