下载文档原格式
运筹学算法matlab程序西北工业大学数学系2009级1.顺向Dijkstra 算法M=[ 0 5 9 Inf Inf Inf InfInf 0 Inf Inf 12 Inf InfInf 3 0 15 Inf 23 InfInf 6 Inf 0 Inf 8 7Inf 12 Inf 5 0 Inf 14Inf Inf Inf Inf Inf 0 10Inf Inf Inf Inf Inf Inf 0];first=1;last=7;[m,n]=size(M);L=zeros(1,m);symbol=zeros(1,m);direction=zeros(1,m);for i=1:mif(i~=first)L(i)=inf;enddirection(i)=first;endjudge=1;while judgefor i=1:mif(symbol(i)==0)min=L(i);temporary=i;breakendendfor i=1:mif(symbol(i)==0)if(L(i)<min)min=L(i);temporary=i;endendendk=temporary;for j=1:mif(symbol(1,j)==0)if(M(k,j)==inf)continue;elseif(L(k)+M(k,j)<L(j))L(j)=L(k)+M(k,j);direction(j)=k;endendendendsymbol(k)=1;num=0;for i=1:mif(symbol(i)==1)num=num+1;endendif(num==m)judge=0;endendp=last;arrow=zeros(1,m);arrow(1)=last;i=2;while p~=firstarrow(1,i)=direction(p);i=i+1;p=direction(p);enddistance=L(last);M=[ 0 5 9 Inf Inf Inf Inf Inf 0 Inf Inf 12 Inf InfInf 3 0 15 Inf 23 Inf Inf 6 Inf 0 Inf 8 7 Inf 12 Inf 5 0 Inf 14 Inf Inf Inf Inf Inf 0 10Inf Inf Inf Inf Inf Inf 0]; [m,n]=size(M);first=1;last=7;L=zeros(1,m);direction=zeros(1,m);symbol=zeros(1,m);for i=1:mdirection(i)=last;if(i~=last)L(i)=inf;endendjudge=1;while judgefor i=1:mif(symbol(i)==0)min=L(i);temporary=i;breakendendfor i=1:mif(symbol(i)==0)if(L(i)<min)min=L(i);temporary=i;endendendk=temporary;for i=1:mif(M(i,k)==inf)continueelseif(M(i,k)+L(k)<L(i))L(i)=L(k)+M(i,k);direction(i)=k;endendendsymbol(k)=1;sum=0;for i=1:mif(symbol(i)==1)sum=sum+1;endendif(sum==m)judge=0;endendp=first;i=2;arrow=zeros(1,m);arrow(1)=first;while p~=lastarrow(i)=direction(p);i=i+1;p=direction(p);endd=[0 7 5 12 inf infinf 0 inf 3 inf infinf inf 0 6 inf 1512 inf 6 0 inf 86 inf 13 inf 0 infinf 4 15 inf 9 0];[m,n]=size(d);p=zeros(m,n);for i=1:np(:,i)=i;endfor k=1:nfor i=1:mfor j=1:nif(d(i,k)+d(k,j)<d(i,j))d(i,j)=d(i,k)+d(k,j);p(i,j)=p(i,k);endendendend4.仿floyd 算法d=[inf 6 0 4 0 0 00 inf 0 0 5 0 04 7 inf 0 05 00 0 4 inf 0 3 00 0 2 0 inf 0 00 0 0 0 4 inf 50 0 0 0 6 0 inf];[m,n]=size(d);first=1;last=7;direction=zeros(m,m);for i=1:mdirection(:,i)=i;endfor i=1:mfor j=1:mfor k=1:msmall=min(d(i,k),d(k,j));if d(i,j)<smalld(i,j)=small;direction(i,j)=direction(i,k);endendendendarrow=zeros(1,m);arrow(1)=first;i=2;p=first;while p~=lastp=direction(p,last);arrow(i)=p;i=i+1;end—dijkstra算法d=[0 inf 3 5 inf10 0 14 inf 8inf inf 0 7 -6inf inf inf 0 infinf inf inf -1 0];[m,n]=size(d);first=2;last=4;L=zeros(1,n);z=zeros(m,n);symbol=zeros(1,n);direction=zeros(1,n);for i=1:nfor j=1:mif d(i,j)~=0if d(i,j)~=infz(i,j)=1;endendenddirection(i)=first;if i~=firstL(i)=inf;endendjudge=1;while judgemini=10;for j=1:nif symbol(j)==0sum=0;for i=1:mp=z(i,j)*(1-symbol(i));sum=sum+p;endif(sum==0)mini=j;breakendendendfor j=1:nif symbol(j)==0&&z(mini,j)==1if L(mini)+d(mini,j)<L(j)L(j)=L(mini)+d(mini,j);direction(j)=mini;endendendsymbol(mini)=1;num=0;for i=1:nif symbol(i)==1num=num+1;endendif num==m;judge=0;endendarrow=zeros(1,m);p=last;arrow(1)=last;i=2;while p~=firstp=direction(p);arrow(i)=p;i=i+1;end—dijkstra算法d=[0 inf 3 5 inf10 0 14 inf 8inf inf 0 7 -6inf inf inf 0 infinf inf inf -1 0];[m,n]=size(d);first=2;last=4;L=zeros(1,n);z=zeros(m,n);symbol=zeros(1,n);direction=zeros(1,n);for i=1:nfor j=1:mif d(i,j)~=0if d(i,j)~=infz(i,j)=1;endendenddirection(i)=last;if i~=lastL(i)=inf;endendjudge=1;while judgemini=10;for i=1:nif symbol(i)==0sum=0;for j=1:mp=z(i,j)*(1-symbol(j));sum=sum+p;endif(sum==0)mini=i;breakendendendfor i=1:nif symbol(i)==0&&z(i,mini)==1if L(mini)+d(i,mini)<L(i)L(i)=L(mini)+d(i,mini);direction(i)=mini;endendendsymbol(mini)=1;num=0;for i=1:nif symbol(i)==1num=num+1;endendif num==m;judge=0;endendarrow=zeros(1,m);p=first;arrow(1)=first;i=2;while p~=lastp=direction(p);arrow(i)=p;i=i+1;endM=[ 0 17 11 inf inf inf17 0 13 12 28 1511 13 0 inf 19 infinf 12 inf 0 inf 16inf 28 19 inf 0 10inf 15 inf 16 10 0];[m,n]=size(M);X=zeros(m,n);Y=zeros(m);Z=zeros(m);Y(1)=1;for i=2:mZ(i)=i;endjudge=1;while judgefor i=1:mif(Y(i)~=0)for j=1:mif(Z(j)~=0)min=M(i,j);a=i;b=j;endendendendfor i=1:mif(Y(i)~=0)for j=1:mif(Z(j)~=0)if(M(i,j)<min)min=M(i,j);a=i;b=j;endendendendendY(b)=b;Z(b)=0;X(a,b)=1;X(b,a)=1;c=0;for i=1:mif(Y(i)~=0)c=c+1;endendif(c==m)judge=0;endend网络最大流Ford—Fulkersen算法d=[inf 12 17 0 0 00 inf 0 8 0 00 6 inf 0 12 00 0 5 inf 0 150 0 0 4 inf 90 0 0 0 0 inf];[m,n]=size(d);X=zeros(m,n);first=1;last=6;recognize=1;while recognizeL=zeros(1,m);L(first)=inf;direction=ones(1,m);symbol=zeros(1,m);judge=1;while judgefor i=1:mif symbol(i)==0big=L(i);k=i;break;endendfor i=1:mif symbol(i)==0if L(i)>bigbig=L(i);k=i;endendendif k==nif L(n)==0breakendelsefor j=1:mif d(k,j)>0u=min(L(k),d(k,j)-X(k,j));if u>L(j)L(j)=u;direction(j)=k;endelseif d(j,k)>0u=min(L(k),X(j,k));if u>L(j)L(j)=u;direction(j)=k;endendendendendsymbol(k)=1;num=0;for i=1:mif symbol(i)==1num=num+1;endendif num==mjudge=0;endendafter=last;before=after;while before~=firstbefore=direction(after);if d(before,after)>0X(before,after)=X(before,after)+L(n); elseX(before,after)=X(before,after)-L(n); endafter=before;endif L(m)==0recognize=0;end end。
图论实验三个案例单源最短路径问题 1.1 Dijkstra 算法Dijkstra 算法是解单源最短路径问题的一个贪心算法。
其基本思想是，设置 一个顶点集合S 并不断地作贪心选择来扩充这个集合。
一个顶点属于集合S 当且 仅当从源到该顶点的最短路径长度已知。
设 v 是图中的一个顶点，记l(v)为顶点 v 到源点V 1的最短距离，V i,V jV ，若(V i,V j)E ，记“到百的权w 。
Dijkstra 算法：① S {V J I(V J 0 ； V V {可 1(V ) i i S V {V J ；J7JJJ7②S，停止，否则转③；l(v) min{ l(v) , d(V j ,v)}V j S④ 存在Vi 1，使l (V i l) min{l(V)}，V S ；⑤SSU{v i 1}S S {v i 1}i i 1实际上，Dijkstra 算法也是最优化原理的应用：如果V 1V 2LV n1Vn是从V1到Vn的最短路径，贝UV 1V 2L Vn1也必然是从V1到Vn 1的最优路径。
在下面的MATLA 实现代码中，我们用到了距离矩阵，矩阵第 i 行第j 行元 素表示顶点Vi到Vj的权Wj，若v 到V j无边，则W ijrealmax，其中realmax 是 MATLA 常量，表示最大的实数(1.7977e+308)function re=Dijkstra(ma)%用Dijkstra 算法求单源最短路径%俞入参量ma是距离矩阵%输出参量是一个三行n 列矩阵，每列表示顶点号及顶点到源的最短距离和前顶点n=size(ma,1);% 得到距离矩阵的维数s=ones(1,n);s(1)=0;% 标记集合S和S 的补r=zeros(3,n);r(1,:)=1:n;r(2,2:end)=realmax;% 初始化for i=2:n;% 控制循环次数mm=realmax;for j=find(s==0);% 集合S中的顶点for k=find(s==1);% 集合S补中的顶点if(r(2,j)+ma(j,k)<r(2,k))r(2,k)=r(2,j)+ma(j,k);r(3,k)=j;endif(mm>r(2,k))mm=r(2,k);t=k;endendends(1,t)=0;%找到最小的顶点加入集合Send re=r;1.2动态规划求解最短路径动态规划是美国数学家 Richard Bellman 在1951年提出来的分析一类多阶 段决策过程的最优化方法，在工程技术、工业生产、经济管理、军事及现代化控 制工程等方面均有着广泛的应用。
MATLAB 智能算法30个案例分析第 1章1、案例背景遗传算法（Genetic Algorithm，GA）是一种进化算法，其基本原理是仿效生物界中的“物竞天择、适者生存”的演化法则。
遗传算法的做法是把问题参数编码为染色体，再利用迭代的方式进行选择、交叉以及变异等运算来交换种群中染色体的信息，最终生成符合优化目标的染色体。
在遗传算法中，染色体对应的是数据或数组，通常是由一维的串结构数据来表示，串上各个位置对应基因的取值。
基因组成的串就是染色体，或者叫基因型个体( Individuals)。
一定数量的个体组成了群体（Population)。
群体中个体的数目称为群体大小（Population Size），也叫群体规模。
而各个个体对环境的适应程度叫做适应度( Fitness)。
2、案例目录：1.1理论基础1.1.1遗传算法概述1.编码2.初始群体的生成3.适应度评估4.选择5.交叉6.变异1.1.2设菲尔德遗传算法工具箱1.工具箱简介2.工具箱添加1.2案例背景1.2.1问题描述1.简单一元函数优化2.多元函数优化1.2.2解决思路及步骤1.3 MATLAB程序实现1.3.1工具箱结构1.3.2遗传算法中常用函数1.创建种群函数—crtbp2.适应度计算函数—ranking3.选择函数—select4.交叉算子函数—recombin5.变异算子函数—mut6.选择函数—reins7.实用函数—bs2rv8.实用函数—rep1.3.3遗传算法工具箱应用举例1.简单一元函数优化2.多元函数优化1.4延伸阅读1.5参考文献3、主程序：1.简单一元函数优化：clcclear allclose all%%画出函数图figure(1);hold on;lb=1;ub=2; %函数自变量范围【1,2】ezplot('sin(10*pi*X)/X',[lb,ub]); %画出函数曲线xlabel('自变量/X')ylabel('函数值/Y')%%定义遗传算法参数NIND=40; %个体数目MAXGEN=20; %最大遗传代数PRECI=20; %变量的二进制位数GGAP=0.95; %代沟px=0.7; %交叉概率pm=0.01; %变异概率trace=zeros(2,MAXGEN); %寻优结果的初始值FieldD=[PRECI;lb;ub;1;0;1;1]; %区域描述器Chrom=crtbp(NIND,PRECI); %初始种群%%优化gen=0; %代计数器X=bs2rv(Chrom,FieldD); %计算初始种群的十进制转换ObjV=sin(10*pi*X)./X; %计算目标函数值while gen<MAXGENFitnV=ranking(ObjV); %分配适应度值SelCh=select('sus',Chrom,FitnV,GGAP); %选择SelCh=recombin('xovsp',SelCh,px); %重组SelCh=mut(SelCh,pm); %变异X=bs2rv(SelCh,FieldD); %子代个体的十进制转换ObjVSel=sin(10*pi*X)./X; %计算子代的目标函数值[Chrom,ObjV]=reins(Chrom,SelCh,1,1,ObjV,ObjVSel); %重插入子代到父代，得到新种群X=bs2rv(Chrom,FieldD);gen=gen+1; %代计数器增加%获取每代的最优解及其序号，Y为最优解,I为个体的序号[Y,I]=min(ObjV);trace(1,gen)=X(I); %记下每代的最优值trace(2,gen)=Y; %记下每代的最优值endplot(trace(1,:),trace(2,:),'bo'); %画出每代的最优点grid on;plot(X,ObjV,'b*'); %画出最后一代的种群hold off%%画进化图figure(2);plot(1:MAXGEN,trace(2,:));grid onxlabel('遗传代数')ylabel('解的变化')title('进化过程')bestY=trace(2,end);bestX=trace(1,end);fprintf(['最优解:\nX=',num2str(bestX),'\nY=',num2str(bestY),'\n'])2.多元函数优化clcclear allclose all%%画出函数图figure(1);lbx=-2;ubx=2; %函数自变量 x范围【-2,2】lby=-2;uby=2; %函数自变量 y范围【-2,2】ezmesh('y*sin(2*pi*x)+x*cos(2*pi*y)',[lbx,ubx,lby,uby],50); %画出函数曲线hold on;%%定义遗传算法参数NIND=40; %个体数目MAXGEN=50; %最大遗传代数PRECI=20; %变量的二进制位数GGAP=0.95; %代沟px=0.7; %交叉概率pm=0.01; %变异概率trace=zeros(3,MAXGEN); %寻优结果的初始值FieldD=[PRECI PRECI;lbx lby;ubx uby;1 1;0 0;1 1;1 1]; %区域描述器Chrom=crtbp(NIND,PRECI*2); %初始种群%%优化gen=0; %代计数器XY=bs2rv(Chrom,FieldD); %计算初始种群的十进制转换X=XY(:,1);Y=XY(:,2);ObjV=Y.*sin(2*pi*X)+X.*cos(2*pi*Y); %计算目标函数值while gen<MAXGENFitnV=ranking(-ObjV); %分配适应度值SelCh=select('sus',Chrom,FitnV,GGAP); %选择SelCh=recombin('xovsp',SelCh,px); %重组SelCh=mut(SelCh,pm); %变异XY=bs2rv(SelCh,FieldD); %子代个体的十进制转换X=XY(:,1);Y=XY(:,2);ObjVSel=Y.*sin(2*pi*X)+X.*cos(2*pi*Y); %计算子代的目标函数值[Chrom,ObjV]=reins(Chrom,SelCh,1,1,ObjV,ObjVSel); %重插入子代到父代，得到新种群XY=bs2rv(Chrom,FieldD);gen=gen+1; %代计数器增加%获取每代的最优解及其序号，Y为最优解,I为个体的序号[Y,I]=max(ObjV);trace(1:2,gen)=XY(I,:); %记下每代的最优值trace(3,gen)=Y; %记下每代的最优值endplot3(trace(1,:),trace(2,:),trace(3,:),'bo'); %画出每代的最优点grid on;plot3(XY(:,1),XY(:,2),ObjV,'bo'); %画出最后一代的种群hold off%%画进化图figure(2);plot(1:MAXGEN,trace(3,:));grid onxlabel('遗传代数')ylabel('解的变化')title('进化过程')bestZ=trace(3,end);bestX=trace(1,end);bestY=trace(2,end);fprintf(['最优解:\nX=',num2str(bestX),'\nY=',num2str(bestY),'\nZ=',num2str(bestZ), '\n'])第 2章基于遗传算法和非线性规划的函数寻优算法1.1案例背景1.1.1非线性规划方法非线性规划是 20世纪 50年代才开始形成的一门新兴学科。
Dijkstra 算法：用矩阵n n a ⨯（n 为顶点个数）存放各边权的邻接矩阵，行向量pb 、1index 、2index、d 分别用来存放P 标号信息、标号顶点顺序、标号顶点索引、最短通路的值。
其中分量⎩⎨⎧=顶点未标号当第顶点已标号当第i i i pb 01)(；)(2i index存放始点到第i 点最短通路中第i 顶点前一顶点的序号；)(i d 存放由始点到第i 点最短通路的值。
求第一个城市到其它城市的最短路径的Matlab 程序如下： clear; clc; M=10000;a(1,:)=[0,50,M,40,25,10]; a(2,:)=[zeros(1,2),15,20,M,25]; a(3,:)=[zeros(1,3),10,20,M]; a(4,:)=[zeros(1,4),10,25]; a(5,:)=[zeros(1,5),55]; a(6,:)=zeros(1,6); a=a+a';pb(1:length(a))=0;pb(1)=1;index1=1;index2=ones(1,length(a)); d(1:length(a))=M;d(1)=0;temp=1; while sum(pb)<length(a) tb=find(pb==0);d(tb)=min(d(tb),d(temp)+a(temp,tb)); tmpb=find(d(tb)==min(d(tb))); temp=tb(tmpb(1)); pb(temp)=1;index1=[index1,temp];index=index1(find(d(index1)==d(temp)-a(temp,index1))); if length(index)>=2 index=index(1); endindex2(temp)=index; endd, index1, index2%dijkstra 最短路算法通用程序，用于求从起始点s 到其它各点的最短路%D 为赋权邻接矩阵，d 为s 到其它各点最短路径的长度，DD 记载了最短路径生成树 function [d,DD]=dijkstra_aiwa(D,s) [m,n]=size(D); d=inf.*ones(1,m); d(1,s)=0;dd=zeros(1,m);dd(1,s)=1;y=s;DD=zeros(m,m);DD(y,y)=1;counter=1;while length(find(dd==1))<mfor i=1:mif dd(i)==0d(i)=min(d(i),d(y)+D(y,i)); endendddd=inf;for i=1:mif dd(i)==0&&d(i)<dddddd=d(i);endendyy=find(d==ddd);counter=counter+1;DD(y,yy(1,1))=counter;DD(yy(1,1),y)=counter;y=yy(1,1);dd(1,y)=1;endFloyd算法：Matlab程序如下：clear;clc;M=10000;a(1,:)=[0,50,M,40,25,10];a(2,:)=[zeros(1,2),15,20,M,25];a(3,:)=[zeros(1,3),10,20,M];a(4,:)=[zeros(1,4),10,25];a(5,:)=[zeros(1,5),55];a(6,:)=zeros(1,6);b=a+a';path=zeros(length(b));for k=1:6for i=1:6for j=1:6if b(i,j)>b(i,k)+b(k,j)b(i,j)=b(i,k)+b(k,j);path(i,j)=k;end end end end b, pathprim 算法构造最小生成树：prim 算法如下：（i ）}{1v P =，Φ=Q ； （ii ）while V P =~},,min(P V v P p w pv pv -∈∈= }{v P P += }{pv Q Q += end用prim 算法求右图的最小生成树。
MATLAB 智能算法30个案例分析第1 章1、案例背景遗传算法（Genetic Algorithm，GA）是一种进化算法，其基本原理是仿效生物界中的“物竞天择、适者生存”的演化法则。
遗传算法的做法是把问题参数编码为染色体，再利用迭代的方式进行选择、交叉以及变异等运算来交换种群中染色体的信息，最终生成符合优化目标的染色体。
在遗传算法中，染色体对应的是数据或数组，通常是由一维的串结构数据来表示，串上各个位置对应基因的取值。
基因组成的串就是染色体，或者叫基因型个体( Individuals) 。
一定数量的个体组成了群体（Population)。
群体中个体的数目称为群体大小（Population Size），也叫群体规模。
而各个个体对环境的适应程度叫做适应度( Fitness) 。
2、案例目录：1.1 理论基础1.1.1 遗传算法概述1. 编码2. 初始群体的生成3. 适应度评估4. 选择5. 交叉6. 变异1.1.2 设菲尔德遗传算法工具箱1. 工具箱简介2. 工具箱添加1.2 案例背景1.2.1 问题描述1. 简单一元函数优化2. 多元函数优化1.2.2 解决思路及步骤1.3 MATLAB程序实现1.3.1 工具箱结构1.3.2 遗传算法中常用函数1. 创建种群函数—crtbp2. 适应度计算函数—ranking3. 选择函数—select4. 交叉算子函数—recombin5. 变异算子函数—mut6. 选择函数—reins7. 实用函数—bs2rv8. 实用函数—rep1.3.3 遗传算法工具箱应用举例1. 简单一元函数优化2. 多元函数优化1.4 延伸阅读1.5 参考文献3、主程序：1. 简单一元函数优化：clcclear allclose all%% 画出函数图figure(1);hold on;lb=1;ub=2; %函数自变量范围【1,2】ezplot('sin(10*pi*X)/X',[lb,ub]); %画出函数曲线xlabel('自变量/X')ylabel('函数值/Y')%% 定义遗传算法参数NIND=40; %个体数目MAXGEN=20; %最大遗传代数PRECI=20; %变量的二进制位数GGAP=0.95; %代沟px=0.7; %交叉概率pm=0.01; %变异概率trace=zeros(2,MAXGEN); %寻优结果的初始值FieldD=[PRECI;lb;ub;1;0;1;1]; %区域描述器Chrom=crtbp(NIND,PRECI); %初始种群%% 优化gen=0; %代计数器X=bs2rv(Chrom,FieldD); %计算初始种群的十进制转换ObjV=sin(10*pi*X)./X; %计算目标函数值while gen<MAXGENFitnV=ranking(ObjV); %分配适应度值SelCh=select('sus',Chrom,FitnV,GGAP); %选择SelCh=recombin('xovsp',SelCh,px); %重组SelCh=mut(SelCh,pm); %变异X=bs2rv(SelCh,FieldD); %子代个体的十进制转换ObjVSel=sin(10*pi*X)./X; %计算子代的目标函数值[Chrom,ObjV]=reins(Chrom,SelCh,1,1,ObjV,ObjVSel); %重插入子代到父代，得到新种群X=bs2rv(Chrom,FieldD);gen=gen+1; %代计数器增加%获取每代的最优解及其序号，Y为最优解,I为个体的序号[Y,I]=min(ObjV);trace(1,gen)=X(I); %记下每代的最优值trace(2,gen)=Y; %记下每代的最优值endplot(trace(1,:),trace(2,:),'bo'); %画出每代的最优点grid on;plot(X,ObjV,'b*'); %画出最后一代的种群hold off%% 画进化图figure(2);plot(1:MAXGEN,trace(2,:));grid onxlabel('遗传代数')ylabel('解的变化')title('进化过程')bestY=trace(2,end);bestX=trace(1,end);fprintf(['最优解:\nX=',num2str(bestX),'\nY=',num2str(bestY),'\n'])2. 多元函数优化clcclear allclose all%% 画出函数图figure(1);lbx=-2;ubx=2; %函数自变量x范围【-2,2】lby=-2;uby=2; %函数自变量y范围【-2,2】ezmesh('y*sin(2*pi*x)+x*cos(2*pi*y)',[lbx,ubx,lby,uby],50); %画出函数曲线hold on;%% 定义遗传算法参数NIND=40; %个体数目MAXGEN=50; %最大遗传代数PRECI=20; %变量的二进制位数GGAP=0.95; %代沟px=0.7; %交叉概率pm=0.01; %变异概率trace=zeros(3,MAXGEN); %寻优结果的初始值FieldD=[PRECI PRECI;lbx lby;ubx uby;1 1;0 0;1 1;1 1]; %区域描述器Chrom=crtbp(NIND,PRECI*2); %初始种群%% 优化gen=0; %代计数器XY=bs2rv(Chrom,FieldD); %计算初始种群的十进制转换X=XY(:,1);Y=XY(:,2);ObjV=Y.*sin(2*pi*X)+X.*cos(2*pi*Y); %计算目标函数值while gen<MAXGENFitnV=ranking(-ObjV); %分配适应度值SelCh=select('sus',Chrom,FitnV,GGAP); %选择SelCh=recombin('xovsp',SelCh,px); %重组SelCh=mut(SelCh,pm); %变异XY=bs2rv(SelCh,FieldD); %子代个体的十进制转换X=XY(:,1);Y=XY(:,2);ObjVSel=Y.*sin(2*pi*X)+X.*cos(2*pi*Y); %计算子代的目标函数值[Chrom,ObjV]=reins(Chrom,SelCh,1,1,ObjV,ObjVSel); %重插入子代到父代，得到新种群XY=bs2rv(Chrom,FieldD);gen=gen+1; %代计数器增加%获取每代的最优解及其序号，Y为最优解,I为个体的序号[Y,I]=max(ObjV);trace(1:2,gen)=XY(I,:); %记下每代的最优值trace(3,gen)=Y; %记下每代的最优值endplot3(trace(1,:),trace(2,:),trace(3,:),'bo'); %画出每代的最优点grid on;plot3(XY(:,1),XY(:,2),ObjV,'bo'); %画出最后一代的种群hold off%% 画进化图figure(2);plot(1:MAXGEN,trace(3,:));grid onxlabel('遗传代数')ylabel('解的变化')title('进化过程')bestZ=trace(3,end);bestX=trace(1,end);bestY=trace(2,end);fprintf(['最优解:\nX=',num2str(bestX),'\nY=',num2str(bestY),'\nZ=',num2str(bestZ), '\n']) 第2 章基于遗传算法和非线性规划的函数寻优算法1.1案例背景1.1.1 非线性规划方法非线性规划是20世纪50年代才开始形成的一门新兴学科。
超全的图论程序关注微信公众号“超级数学建模”，教你做有料、有趣的数模人程序一：可达矩阵算法function P=dgraf(A)n=size(A,1);P=A;for i=2:nP=P+A^i;endP(P~=0)=1;P;程序二：关联矩阵和邻接矩阵互换算法function W=incandadf(F,f)if f==0m=sum(sum(F))/2;n=size(F,1);W=zeros(n,m);k=1;for i=1:nfor j=i:nif F(i,j)~=0W(i,k)=1;W(j,k)=1;k=k+1;endendendelseif f==1m=size(F,2);n=size(F,1);W=zeros(n,n);for i=1:ma=find(F(:,i)~=0);W(a(1),a(2))=1;W(a(2),a(1))=1;endelsefprint('Please imput the right value of f');endW;程序三：有向图关联矩阵和邻接矩阵互换算法function W=mattransf(F,f)if f==0m=sum(sum(F));n=size(F,1);W=zeros(n,m);k=1;for i=1:nfor j=i:nif F(i,j)~=0W(i,k)=1;W(j,k)=-1;k=k+1;endendendelseif f==1m=size(F,2);n=size(F,1);W=zeros(n,n);for i=1:ma=find(F(:,i)~=0);if F(a(1),i)==1W(a(1),a(2))=1;elseW(a(2),a(1))=1;endendelsefprint('Please imput the right value of f'); endW;第二讲：最短路问题程序一：Dijkstra算法（计算两点间的最短路）function [l,z]=Dijkstra(W)n = size (W,1);for i = 1 :nl(i)=W(1,i);z(i)=0;endi=1;while i<=nfor j =1 :nif l(i)>l(j)+W(j,i)l(i)=l(j)+W(j,i);z(i)=j-1;if j<ii=j-1;endendendi=i+1;end程序二：floyd算法（计算任意两点间的最短距离）function [d,r]=floyd(a)n=size(a,1);d=a;for i=1:nfor j=1:nr(i,j)=j;endendr;for k=1:nfor i=1:nfor j=1:nif d(i,k)+d(k,j)<d(i,j)d(i,j)=d(i,k)+d(k,j); r(i,j)=r(i,k);endendendend程序三：n2short.m 计算指定两点间的最短距离function [P u]=n2short(W,k1,k2)n=length(W);U=W;m=1;while m<=nfor i=1:nfor j=1:nif U(i,j)>U(i,m)+U(m,j)U(i,j)=U(i,m)+U(m,j);endendendm=m+1;endu=U(k1,k2);P1=zeros(1,n);k=1;P1(k)=k2;V=ones(1,n)*inf;kk=k2;while kk~=k1for i=1:nV(1,i)=U(k1,kk)-W(i,kk);if V(1,i)==U(k1,i)P1(k+1)=i;kk=i;k=k+1;endendendk=1;wrow=find(P1~=0);for j=length(wrow):-1:1P(k)=P1(wrow(j));k=k+1;endP;程序四、n1short.m(计算某点到其它所有点的最短距离)function[Pm D]=n1short(W,k)n=size(W,1);D=zeros(1,n);for i=1:n[P d]=n2short(W,k,i);Pm{i}=P;D(i)=d;end程序五：pass2short.m(计算经过某两点的最短距离) function [P d]=pass2short(W,k1,k2,t1,t2)[p1 d1]=n2short(W,k1,t1);[p2 d2]=n2short(W,t1,t2);[p3 d3]=n2short(W,t2,k2);dt1=d1+d2+d3;[p4 d4]=n2short(W,k1,t2);[p5 d5]=n2short(W,t2,t1);[p6 d6]=n2short(W,t1,k2);dt2=d4+d5+d6;if dt1<dt2d=dt1;P=[p1 p2(2:length(p2)) p3(2:length(p3))]; elsed=dt1;p=[p4 p5(2:length(p5)) p6(2:length(p6))]; endP;d;第三讲：最小生成树程序一：最小生成树的Kruskal算法function [T c]=krusf(d,flag)if nargin==1n=size(d,2);m=sum(sum(d~=0))/2;b=zeros(3,m);k=1;for i=1:nfor j=(i+1):nif d(i,j)~=0b(1,k)=i;b(2,k)=j;b(3,k)=d(i,j);k=k+1;endendendelseb=d;endn=max(max(b(1:2,:)));m=size(b,2);[B,i]=sortrows(b',3);B=B';c=0;T=[];k=1;t=1:n;for i=1:mif t(B(1,i))~=t(B(2,i))T(1:2,k)=B(1:2,i);c=c+B(3,i);k=k+1;tmin=min(t(B(1,i)),t(B(2,i)));tmax=max(t(B(1,i)),t(B(2,i)));for j=1:nif t(j)==tmaxt(j)=tmin;endendendif k==nbreak;endendT;c;程序二：最小生成树的Prim算法function [T c]=Primf(a)l=length(a);a(a==0)=inf;k=1:l;listV(k)=0;listV(1)=1;e=1;while (e<l)min=inf;for i=1:lif listV(i)==1for j=1:lif listV(j)==0 & min>a(i,j)min=a(i,j);b=a(i,j);s=i;d=j;endendendendlistV(d)=1;distance(e)=b;source(e)=s;destination(e)=d;e=e+1;endT=[source;destination];for g=1:e-1c(g)=a(T(1,g),T(2,g));endc;另外两种程序最小生成树程序1（prim 算法构造最小生成树）a=[inf 50 60 inf inf inf inf;50 inf inf 65 40 inf inf;60 inf inf 52 inf inf 45;...inf 65 52 inf 50 30 42;inf 40 inf 50 inf 70 inf;inf inf inf 30 70 inf inf;...inf inf 45 42 inf inf inf];result=[];p=1;tb=2:length(a);while length(result)~=length(a)-1temp=a(p,tb);temp=temp(:);d=min(temp);[jb,kb]=find(a(p,tb)==d);j=p(jb(1));k=tb(kb(1));result=[result,[j;k;d]];p=[p,k];tb(find(tb==k))=[];endresult最小生成树程序2（Kruskal 算法构造最小生成树）clc;clear;a(1,2)=50; a(1,3)=60; a(2,4)=65; a(2,5)=40;a(3,4)=52;a(3,7)=45; a(4,5)=50; a(4,6)=30;a(4,7)=42; a(5,6)=70;[i,j,b]=find(a);data=[i';j';b'];index=data(1:2,:);loop=max(size(a))-1;result=[];while length(result)<looptemp=min(data(3,:));flag=find(data(3,:)==temp);flag=flag(1);v1=data(1,flag);v2=data(2,flag);if index(1,flag)~=index(2,flag)result=[result,data(:,flag)];endindex(find(index==v2))=v1;data(:,flag)=[];index(:,flag)=[];endresult第四讲：Euler图和Hamilton图程序一：Fleury算法（在一个Euler图中找出Euler环游）注：包括三个文件；fleuf1.m, edf.m, flecvexf.m function [T c]=fleuf1(d)%注：必须保证是Euler环游，否则输出T=0,c=0n=length(d);b=d;b(b==inf)=0;b(b~=0)=1;m=0;a=sum(b);eds=sum(a)/2;ed=zeros(2,eds);vexs=zeros(1,eds+1);matr=b;for i=1:nif mod(a(i),2)==1m=m+1;endendif m~=0fprintf('there is not exit Euler path.\n') T=0;c=0;endif m==0vet=1;flag=0;t1=find(matr(vet,:)==1);for ii=1:length(t1)ed(:,1)=[vet,t1(ii)];vexs(1,1)=vet;vexs(1,2)=t1(ii);matr(vexs(1,2),vexs(1,1))=0;flagg=1;tem=1;while flagg[flagg ed]=edf(matr,eds,vexs,ed,tem); tem=tem+1;if ed(1,eds)~=0 & ed(2,eds)~=0T=ed;T(2,eds)=1;c=0;for g=1:edsc=c+d(T(1,g),T(2,g));endflagg=0;break;endendendendfunction[flag ed]=edf(matr,eds,vexs,ed,tem)flag=1;for i=2:eds[dvex f]=flecvexf(matr,i,vexs,eds,ed,tem);if f==1flag=0;break;endif dvex~=0ed(:,i)=[vexs(1,i) dvex];vexs(1,i+1)=dvex;matr(vexs(1,i+1),vexs(1,i))=0;elsebreak;endendfunction [dvex f]=flecvexf(matr,i,vexs,eds,ed,temp) f=0;edd=find(matr(vexs(1,i),:)==1);dvex=0;dvex1=[];ded=[];if length(edd)==1dvex=edd;elsedd=1;dd1=0;kkk=0;for kk=1:length(edd)m1=find(vexs==edd(kk));if sum(m1)==0dvex1(dd)=edd(kk);dd=dd+1;dd1=1;elsekkk=kkk+1;endendif kkk==length(edd)tem=vexs(1,i)*ones(1,kkk);edd1=[tem;edd];for l1=1:kkklt=0;ddd=1;for l2=1:edsif edd1(1:2,l1)==ed(1:2,l2)lt=lt+1;endendif lt==0ded(ddd)=edd(l1);ddd=ddd+1;endendendif temp<=length(dvex1)dvex=dvex1(temp);elseif temp>length(dvex1) & temp<=length(ded)dvex=ded(temp);elsef=1;endend程序二：Hamilton改良圈算法（找出比较好的Hamilton路）function [C d1]= hamiltonglf(v)%d表示权值矩阵%C表示算法最终找到的Hamilton圈。
最短路法MATLAB 程序例1： 某公司在六个城市621,,,c c c 中有分公司，从i c 到j c的直接航程票价记在下述矩阵的),(j i 位置上。
（∞表示无直接航路），请帮助该公司设计一张城市1c 到其它城市间的票价最便宜的路线图。
⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎣⎡∞∞∞∞∞∞055252510550102025251001020402010015252015050102540500 用矩阵n n a ⨯（n 为顶点个数）存放各边权的邻接矩阵，行向量pb 、1index 、2index 、d 分别用来存放P 标号信息、标号顶点顺序、标号顶点索引、最短通路的值。
其中分量⎩⎨⎧=顶点未标号当第顶点已标号当第i i i pb 01)(； )(2i index 存放始点到第i 点最短通路中第i 顶点前一顶点的序号；)(i d 存放由始点到第i 点最短通路的值。
求第一个城市到其它城市的最短路径的Matlab 程序如下：程序一：clc,cleara=zeros(6); %邻接矩阵初始化a(1,2)=50;a(1,4)=40;a(1,5)=25;a(1,6)=10;a(2,3)=15;a(2,4)=20;a(2,6)=25;a(3,4)=10;a(3,5)=20;a(4,5)=10;a(4,6)=25;a(5,6)=55;a=a+a'; %将矩阵数据对称赋予矩阵a(find(a==0))=inf; %找到0值并赋值为infpb(1:length(a))=0;pb(1)=1;index1=1;index2=ones(1,length(a)); d(1:length(a))=inf;d(1)=0;temp=1; %最新的P标号的顶点while sum(pb)<length(a)tb=find(pb==0);d(tb)=min(d(tb),d(temp)+a(temp,tb)); %d指标号顶点索引、tmpb=find(d(tb)==min(d(tb)));temp=tb(tmpb(1)); %可能有多个点同时达到最小值，只取其中的一个pb(temp)=1;index1=[index1,temp]; %intex1指标号顶点顺序temp2=find(d(index1)==d(temp)-a(temp,index1));index2(temp)=index1(temp2(1)); %intex2指标号顶点索引endd, index1, index2程序二clear;clc;M=10000;a(1,:)=[0,50,M,40,25,10];a(2,:)=[zeros(1,2),15,20,M,25];a(3,:)=[zeros(1,3),10,20,M];a(4,:)=[zeros(1,4),10,25];a(5,:)=[zeros(1,5),55];a(6,:)=zeros(1,6);a=a+a';pb(1:length(a))=0;pb(1)=1;d(1:length(a))=M;d(1)=0;temp=1; while sum(pb)<length(a)tb=find(pb==0);d(tb)=min(d(tb),d(temp)+a(temp,tb));tmpb=find(d(tb)==min(d(tb)));temp=tb(tmpb(1)); pb(temp)=1;endd。
图论实验三个案例单源最短路径问题 1.1 Dijkstra 算法Dijkstra 算法是解单源最短路径问题的一个贪心算法。
其基本思想是，设置一个顶点集合S 并不断地作贪心选择来扩充这个集合。
一个顶点属于集合S 当且仅当从源到该顶点的最短路径长度已知。
设v 是图中的一个顶点，记()l v 为顶点v 到源点v 1的最短距离，,i j v v V∀∈，若(,)i j v v E∉，记i v到j v 的权ij w =∞。
Dijkstra 算法：① 1{}S v =,1()0l v =；1{}v V v ∀∉-,()l v =∞,1i =,1{}S V v =-； ② S φ=，停止，否则转③；③ ()min{(),(,)}j l v l v d v v =，j v S∈，v S ∀∈；④ 存在1i v +，使1()min{()}i l v l v +=，v S ∈；⑤1{}i S S v +=U ，1{}i S S v +=-，1i i =+，转②；实际上，Dijkstra 算法也是最优化原理的应用：如果121n nv v v v -L 是从1v 到nv 的最短路径，则121n v v v -L 也必然是从1v 到1n v -的最优路径。
在下面的MATLAB 实现代码中，我们用到了距离矩阵，矩阵第i 行第j 行元素表示顶点i v 到j v 的权ij w ，若i v到j v 无边，则realmax ij w =，其中realmax 是MATLAB 常量，表示最大的实数(1.7977e+308)。
function re=Dijkstra(ma)%用Dijkstra算法求单源最短路径%输入参量ma是距离矩阵%输出参量是一个三行n列矩阵，每列表示顶点号及顶点到源的最短距离和前顶点n=size(ma,1);%得到距离矩阵的维数s=ones(1,n);s(1)=0;%标记集合S和S的补r=zeros(3,n);r(1,:)=1:n;r(2,2:end)=realmax;%初始化for i=2:n;%控制循环次数mm=realmax;for j=find(s==0);%集合S中的顶点for k=find(s==1);%集合S补中的顶点if(r(2,j)+ma(j,k)<r(2,k))r(2,k)=r(2,j)+ma(j,k);r(3,k)=j;endif(mm>r(2,k))mm=r(2,k);t=k;endendends(1,t)=0;%找到最小的顶点加入集合S end re=r;1.2 动态规划求解最短路径动态规划是美国数学家Richard Bellman 在1951年提出来的分析一类多阶段决策过程的最优化方法，在工程技术、工业生产、经济管理、军事及现代化控制工程等方面均有着广泛的应用。
动态规划应用了最佳原理：假设为了解决某一优化问题，需要依次作出n 个决策12,,,nD D D L ，如若这个决策是最优的，对于任何一个整数k ，1<k <n ，不论前面k 个决策是怎样的，以后的最优决策只取决于由前面决策所确定的当前状态，即12,,,k k nD D D ++L 也是最优的。
如图1，从A 1点要铺设一条管道到A 16点，中间必须要经过5个中间站，第一站可以在{ A 2，A 3}中任选一个，第二、三、四、五站可供选择的地点分别是：{ A 4，A 5，A 6，A 7}，{ A 8，A 9，A 10}，{ A 11，A 12，A 13}，{ A 14，A 15}。
连接两地管道的距离用连线上的数字表示，要求选一条从A 1到A 16的铺管线路，使总距离最短。
图1 可选择的管道图解决此问题可以用穷举法，从A1到A16有48条路径，只须比较47次，就可得到最短路径为：A1→A2→A5→A8→A12→A15→A16，最短距离为18。
也可以使用Dijkstra算法。
这里，我们动态规划解决此问题。
注意到最短路径有这样一个特性，即如果最短路径的第k站通过P k，则这一最短路径在由P k出发到达终点的那一部分路径，对于始点为P k到终点的所有可能的路径来说，必定也是距离最短的。
根据最短路径这一特性，启发我们计算时从最后一段开始，从后向前逐步递推的方法，求出各点到A16的最短路径。
在算法中，我们用数组六元数组ss表示中间车站的个数(A1也作为中间车站)，用距离矩阵path表示该图。
为简便起见，把该图看作有向图，各边的方向均为从左到右，则path不是对称矩阵，如path(12,14)=5，而path(14,12)=0(用0表示不通道路)。
用3´16矩阵spath表示算法结果，第一行表示结点序号，第二行表示该结点到终点的最短距离，第三行表示该结点到终点的最短路径上的下一结点序号。
下面给出MATLAB实现算法。
function [scheme] = ShortestPath(path,ss)%利用动态规划求最短路径%path是距离矩阵,ss是车站个数n=size(path,1);%结点个数scheme=zeros(3,n);%构造结果矩阵scheme(1,:)=1:n;%设置结点序号scheme(2,1:n-1)=realmax;%预设距离值k=n-1;%记录第一阶段结点最大序号for i=size(ss,2):-1:1;%控制循环阶段数for j=k:-1:(k-ss(i)+1);%当前阶段结点循环for t=find(path(j,:)>0);%当前结点邻接结点if path(j,t)+scheme(2,t)<scheme(2,j)scheme(2,j)=path(j,t)+scheme(2,t);scheme(3,j)=t;endendendk=k-ss(i);移入下一阶段end先在MATLAB命令窗口中构造距离矩阵path，再输入：>> ShortestPath(path,ss)得到以下结果：1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 18 13 16 13 10 9 12 7 6 8 7 5 9 4 3 02 5 6 8 8 9 10 12 12 12 14 15 15 16 16 0将该结果表示为图，即为图1粗线所示。
棋盘覆盖问题 1.1 问题的提出在一个22k k⨯个方格组成的棋盘中，若恰有一个方格与其他方格不同，则称该方格为一特殊方格，且称该棋盘为一特殊的棋盘。
如图1就是当3k =时的特2所示4种不同形态的L 形骨牌覆盖一个1.2 问题的分析易知，用到的L 型骨牌个数恰为(41)/3k -。
利用分治策略，我们可以设计出解棋盘覆盖问题的一个简捷的算法。
图1 当k =3时的特殊棋盘图2 4种不同形态的L 型骨牌(a)(b)(c)(d)当k >0时，我们将22k k ⨯棋盘分割为4个1122k k --⨯子棋盘如图3两粗实线所示。
特殊方格必位于4个较小子棋盘之一中，其余3个子棋盘中无特殊方格。
为了将这3个无特殊方格的子棋盘转化为特殊棋盘，我们可以用一个L 型骨牌覆盖这3个较小棋盘的会合处，如图4中央L 型骨牌所示，这3个子棋盘上被L 型骨牌覆盖的方格就成为该棋盘上的特殊方格，从而将原问题转化为4个较小规模的棋盘覆盖问题。
递归地使用这种分割，直至棋盘简化为11⨯棋盘。
1.3 算法的MATLAB 实现首先特殊方格在棋盘中的位置可以用一个12⨯的数组sp 表示；对于图2所示的4种L 型骨牌，可用数字1,2,3,4表示；对于特殊棋盘的骨牌覆盖表示，只图3 棋盘分割 图4 关键结点1 2 3须注意到图4所示的关键点，对每个关键点，给定一种L 型骨牌，就能覆盖整个棋盘，所以对于22k k ⨯的特殊棋盘的骨牌覆盖，可用一个(21)(21)k k -⨯-的矩阵表示。
按照这种思想，图4的矩阵表示为：k =4，特殊方格位置为：[1,4]，覆盖矩阵为：1040102040002040302030003000104030204000304030403⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦下面是在MATLAB 中的棋盘覆盖实现程序。
function re = chesscover(k,sp) %解决棋盘的覆盖问题%棋盘为2^k*2^k ，sp 为特殊方格的棋盘位置 global covermatrixcovermatrix=zeros(2^k-1,2^k-1);even1=floor(sp(1,1)/2)*2==sp(1,1);%判断水平位置是否是偶数 even2=floor(sp(1,2)/2)*2==sp(1,2);%判断竖直位置是否是偶数 if even1==1&&even2==0%找出找出特殊方格相对关键结点的位置 i=4; elsei=even1+even2+1;endtempfun(1,1,k,[sp(1,1)-even1,sp(1,2)-even2,i]);re=covermatrix;function tempfun(top,left,k,tp)%子函数,tp为转换后特殊方格在棋盘网络的相对位置global covermatrixif k==1switch tp(1,3)case 1covermatrix(tp(1,1),tp(1,2))=3;case 2covermatrix(tp(1,1),tp(1,2))=4;case 3covermatrix(tp(1,1),tp(1,2))=1;case 4covermatrix(tp(1,1),tp(1,2))=2;endelsehalf=2^(k-1);i=top+half-1;j=left+half-1;if tp(1,1)<iif tp(1,2)<j%特殊方格在左上covermatrix(i,j)=3; %添加类型为3的L型骨牌tempfun(top,left,k-1,tp);tempfun(top,left+half,k-1,[i-1,j+1,4]);tempfun(top+half,left+half,k-1,[i+1,j+1,1] );tempfun(top+half,left,k-1,[i+1,j-1,2]);else %特殊方格在右上covermatrix(i,j)=4;%添加类型为4的L型骨牌tempfun(top,left,k-1,[i-1,j-1,3]);tempfun(top,left+half,k-1,tp);tempfun(top+half,left+half,k-1,[i+1,j+1,1] );tempfun(top+half,left,k-1,[i+1,j-1,2]);endelseif tp(1,2)>j%特殊方格在右下covermatrix(i,j)=1;%添加类型为3的L型骨牌tempfun(top,left,k-1,[i-1,j-1,3]);tempfun(top,left+half,k-1,[i-1,j+1,4]);tempfun(top+half,left+half,k-1,tp);tempfun(top+half,left,k-1,[i+1,j-1,2]);else %特殊方格在左下covermatrix(i,j)=2;%添加类型为4的L型骨牌tempfun(top,left,k-1,[i-1,j-1,3]);tempfun(top,left+half,k-1,[i-1,j+1,4]);tempfun(top+half,left+half,k-1,[i+1,j+1,1] );tempfun(top+half,left,k-1,tp);endendend在MATLAB命令窗口中输入指令chesscover(3,[1,4])将会得到如上面矩阵一样的结果。
