实验内容：
1、哈密顿环问题：
(a)实现基于树的基本搜索算法，主要实现广度优先搜索和深度优先搜索
(b)在树搜索中利用爬山法的思想，考虑在搜索过程中如何选择节点进行展开搜索，设计并实现搜索的“个性化”优化策略
2、最小哈密顿环问题：实现求解最小哈密顿环问题的分支界限算法
实验过程及结果(matlab):
1、N个节点的无向图随机生成
可以按照两种类型生成，第一种为0,1矩阵，1表示两点之间有连接
第二种为带权重矩阵，且为全连接图，0表示无连接，有连接的为正数，但是最大权重小于MaxWeight
MaxWeight=20;%最大权重
edge=zeros(N);%N节点0值方阵
for i=1:N
forj=i+1:N
r=rand;
if type==1
if r>=0.5
if type==1
edge(i,j)=1;
edge(j,i)=1;
end
end
elseif type==2
edge(i,j)=MaxWeight*rand;
edge(j,i)=edge(i,j);
end
end
end
2、DFS深度优先搜索求解哈密顿环
（1）实验步骤
a. 扫描当前节点的可用子节点，按照节点的编号由小到大选取下一个遍历的点
fori=knum(k)+1:n %从该回溯的节点开始
if G(route(k),i)>0&&~visited(i)
t=i;
knum(k)=t;%计算到t
flag=1;
break
else
flag=0;
end
end
b. 若有新节点被加入路径，判断当前路径长度，若长度等于节点数，则判断是否为哈密顿路径，若没有新节点加入，回溯到上一层
if flag %i如果找到新节点则加入到路径中
if k~=n
visited(t)=1;
k=k+1;
route=[route t];
else
break
end
else %如果没找到就回溯
visited(route(k))=0;
route=route(1:size(route,2)-1);
knum(k)=0;
k=k-1;
end
（2）实验结果
以十个节点的图作为例子，邻接矩阵：
0 0 1 0 0 1 0 1 1 0
0 0 0 1 0 0 0 0 1 1
1 0 0 1 0 1 0 1 0 1
0 1 1 0 0 1 1 1 1 0
0 0 0 0 0 0 1 0 1 1
1 0 1 1 0 0 1 1 1 1
0 0 0 1 1 1 0 0 1 1
1 0 1 1 0 1 0 0 0 0
1 1 0 1 1 1 1 0 0 1
0 1 1 0 1 1 1 0 1 0
哈密顿路径：
时间：0.0087s
路径 1 3 4 2 9 5 7 10 6 8 1
（3）性能曲线
用节点个数8-20的无向图进行测试
（4）结果分析
从图中可以看出DFS搜索的时间通常较短，且哈密顿路径的搜索时间并不严格按照节点数增加而增加，所以途中出现随机增大现象，这是由于输入图的边的分布情况不同会对深度优先搜索的运行时间造成严重的影响。
3、BFS广度优先搜索求解哈密顿环
（1）实验步骤
a. 维护一个待搜索节点的队列，每个点都记录到达该点的路径，检索队列中所有结点是否有新的可用子节点
for i=1:size(player,2)
t=[];
for j=1:n %找出能够到达的节点并加入到队列中
ifG(player(i).node,j)>0&&(~player(i).visited(j))
ct=ct+1;
t.node=j;
t.route=[player(i).route j];
t.visited=player(i).visited;
t.visited(j)=1;
clayer=[clayer t];
flag=1;
end
end
end
b. 若有可用新节点，则将可用节点的加入相应路径，并判断该路径是否为哈密顿路径。若没有可用新节点，说明已遍历所有可能的路径，找不到哈密顿回路。
If flag %如果找到新节点
for i=1:size(clayer,2)
if size(clayer(i).route,2)==n %是否有路径长度为N
if G(1,clayer(i).route(size(clayer(i).route,2)))>0%是否是哈密顿环
oroute=clayer(i).route;
return
end
end
end
player=clayer;
clayer=[];
else %没找到哈密顿环
return;
end
（2）实验结果
选取10个节点的无向图作为测试：
0 1 1 0 1 1 1 0 0 0
1 0 1 1 1 1 1 1 0 1
1 1 0 0 1 1 1 1 1 0
0 1 0 0 1 1 1 1 0 1
1 1 1 1 0 0 0 0 1 1
1 1 1 1 0 0 1 0 1 1
1 1 1 1 0 1 0 0 0 0
0 1 1 1 0 0 0 0 1 0
0 0 1 0 1 1 0 1 0 1
0 1 0 1 1 1 0 0 1 0
哈密顿路径：
运行时间： 11.6400
路径：1 3 2 4 6 9 5 10 7 8 1
（3）性能曲线
选取6-12个节点的无向图进行测试，由于13个以上节点的输入图所带来的时间开销极为巨大，因此在这里不列出来了。
（4）结果分析
从图中可以看出广度优先搜索哈密顿回路所花费的时间远远大于深度优先搜索，并且其时间开销受边的分布情况影响不大，主要受节点个数影响，大约为O(n!)，因为在bfs搜索出结果前，已经遍历了几乎全部的路径。
4、爬山发求解哈密顿环
（1）实验步骤
a. 维护一个记录访问节点的栈，若上一次循环并非回溯，则取栈顶元素，找到其所有可用的子节点并计算子节点可用的出度，将子节点按出度排序，压栈，并记录当前节点的可用子节点数。这里优先扩展出度较大的点，理由是出度较大的点通常会带来更多的可行路线，增大在当前路径找到哈密顿回路的可能性。
if flag2 % 如果回溯，则不计算
for i=2:n
count=0;
ifG(route(k),i)>0&&~visited(i) %计算当前节点的出度以及连接点
for j=2:n
if j~=route(k)&&G(i,j)>0&&~visited(j)
count=count+1; %out-degree
end
end
flag=1;
knum(k)=knum(k)+1;%表示该节点有几个子节点需要访问
cnode=cnode+1;%
a(cnode).node=i; %记录节点序号
a(cnode).count=count; %记录几点出度
end
end
end
b. 若a中找到可用的新子节点，则取栈顶节点，加入队列，判断是否为哈密顿回路。若a中没有找到新的子节点，则判断当前队列是否为哈密顿回路，若不为，则需要回溯。回溯时判断路径上上一点是否有其他可用的子节点，若有，取栈顶元素替换当前路径点，若没有，从路径中删除当前元素。
if flag %如果找到新节点
[f,seq]=sort([a.count]);%求出度最大的，爬山策略
for i=1:size(f,2)
t.node=a(seq(i)).node;
t.count=f(i);
s=[t s]; %按照从大到小的顺序压栈
end
end
ifflag %找到新节点
if k~=n %如果当前路径不是全路径，就把当前节点放入路径中
t=s(1).node;
s=s(2:size(s,2)); %弹栈
visited(t)=1;
knum(k)=knum(k)-1;%Kth节点栈中剩余的个数
k=k+1;
route=[route t];%记录路径;
elseif k==n&&G(1,route(k))>0 %如果当前节点的路径长度为n，表示已找到哈密顿环
oroute=route;
break;
end
else %如果没找到新节点
if k==n&&G(1,route(k))>0 %如果当前节点的路径长度为n，表示已找到哈密顿环
oroute=route;
break;
elseif k==1 %如果当前节点的路径长度为1，表示没有哈密顿环
oroute=[1];
break;
end
if knum(k-1)>0 %该层栈中仍有剩余
t=s(1).node;
s=s(2:size(s,2));
visited(t)=1;
visited(route(k))=0; %访问的节点变为未访问
route(k)=t; %记录当前节点已访问
knum(k-1)=knum(k-1)-1;%该层栈中个数减一
knum(k)=0; %访问的个数为零
else %该层栈中个数为零
visited(route(k))=0; %当前节点未访问
route=route(1:size(route,2)-1); %放弃当前节点
k=k-1; %路径长度减一
flag2=0; %开始回溯
end
end
（2）实验结果
以十个节点的图作为输入：
0 0 1 1 0 1 1 0 0 0
0 0 1 1 0 0 0 0 0 1
1 1 0 0 0 1 1 0 0 1
1 1 0 0 0 1 1 0 0 0
0 0 0 0 0 1 0 1 0 0
1 0 1 1 1 0 1 0 0 1
1 0 1 1 0 1 0 1 1 1
0 0 0 0 1 0 1 0 0 0
0 0 0 0 0 0 1 0 0 1
0 1 1 0 0 1 1 0 1 0
哈密顿路径：
运行时间：0.0893
路径：1 6 5 8 7 9 10 3 2 4 1（
3）性能曲线
选取8-20个点的图作为输入
（4）结果分析
由于是基于深度优先的爬山法，从图中可以看出其性能曲线与深度优先的性能曲线有相似处，即输入图的规模不能代表算法的运行时间，算法的实际运行时间受输入图边的分布影响明显。同时爬山法所花费的时间要比单纯的dfs略多，这是由于爬山法会花费一些时间在扩展节点的选择上，而dfs则不需要考虑这种问题，图中两点间生成一条边的概率为0.5，边相对较多，爬山法的选择策略无法较好的发挥。
5、分支限界法求解最小哈密顿环
（1）实验步骤
a. 采用BSF与分支限界法，因此需要维护一个最小堆，其中权重为起点到当前点的代价，堆有插入与提取最小元素两种操作。
%a simple insert to a minimum heap
heap=[oheap new];
n=size(heap,2);
c=n;
while c~=1&&heap(c).weight<heap(floor(c/2)).weight
t=heap(c);
heap(c)=heap(floor(c/2));
heap(floor(c/2))=t;
c=floor(c/2);
end
output=heap;
%simply return the min value from a minmium heap,and update rest oftheheap
n=size(oheap,2);
min=oheap(1);
if n==1
heap=[];
return
end
heap=[oheap(n) oheap(:,2:n-1)];
c=1;
while c*2<=n-1
left=heap(c*2).weight;
if (2*c+1<=n-1)
right=heap(c*2+1).weight;
ifleft<=right&&left>heap(c).weight
t=heap(c*2);
heap(c*2)=heap(c);
heap(c)=t;
c=c*2;
elseifright<left&&right>heap(c).weight
t=heap(c*2+1);
heap(c*2+1)=heap(c);
heap(c)=t;
c=c*2+1;
else
break;
end
else
if left<heap(c).weight
t=heap(c*2);
heap(c*2)=heap(c);
heap(c)=t;
c=c*2;
else
break;
end
end
end
b. 取堆中最小的元素，若当前路径的代价已大于历史最小回路，则停止扩展，进行下一次循环，若小于历史回路，则找出其所有可用子节点，分别计算到这些子节点的权重，加入堆中。
[current heap]=heapmin(heap); %获取对中的最大权重节点
flag=0; %新节点标志，1表示节点需要介入堆中
ifminw>0&¤t.weight>minw %路线权重超过阈值，则不进行拓展
continue;
else
for i=1:n
if G(current.node,i)>0&&(~current.visited(i)) %f寻找新的子节点
t.node=i;
t.route=[current.route i];
t.visited=current.visited;
t.visited(i)=1;
t.weight=current.weight+G(current.node,i);
heap=heapinsert(heap,t);%addnew nodes to heap
flag=1;
end
end
end
c. 如果没有新节点加入，则判断当前路径是否为哈密顿回路，若为哈密顿回路，则比较当前回路与历史最小回路的代价大小，若当前路径代价小于历史最小路径，则更新历史最小路径与代价。
if~flag
if size(current.route,2)==n; %路径长度为n
if G(1,current.route(n))>0 %如果目前是哈密顿环，则更新阈值
oroute=current.route;
minw=current.weight;
oweight=current.weight;
end
end
end
（2）实验结果
以十个节点的完全图作为例子：
最小哈密顿路径：
1 7 2 8 5 6 4 9 10 3 1
代价为：7.0531
（3）性能曲线
（4）结果分析
由于需要遍历所有可能的路径，在节点数大于12时求最小哈密顿路径所花费的时间远大于节点数小于等于12时花费的时间，因此不画入图中。从图中可以看出因为要遍历所有可能的路径，算法所花费的时间与点数量的成指数关系。
详细实现代码：
matlab版本：
1.随机生成图n*n
function [ edge ] = randgraph( N,type )
%UNTITLED5 Summary of this function goes here
% Detailed explanation goes here
MaxWeight=20;
edge=zeros(N);
for i=1:N
for j=i+1:N
r=rand;
if type==1
if r>=0.5
if type==1
edge(i,j)=1;
edge(j,i)=1;
end
end
elseif type==2
edge(i,j)=MaxWeight*rand;
edge(j,i)=edge(i,j);
end
end
end
end
2.深度优先搜索
启动项：
clear;
clc;
N=10;%number of node
G=randgraph(N,1);
tic
route=mydfs(G);
a=toc
if size(route,2)<=1
disp('no hamiltonian path');
else
disp([route route(1)]);
end
核心：
function [ oroute ] = mydfs( G )
%UNTITLED6 Summary of this function goes here
% Detailed explanation goes here
n=size(G,1);
route=[1];
visited=zeros(n,1);%not visited=0,visited=1
k=1;%how many nodes are in current route
visited(1)=1;
knum=zeros(n,1);%kth node of current route's viewed adjective
while k>0
t=0;flag=0;
if k==n&&G(1,route(k))>0 %k==n,judge current route is hamilton route or not
oroute=route;
break;
end
for i=knum(k)+1:n %find next adjective node
if G(route(k),i)>0&&~visited(i)
t=i;
knum(k)=t;
flag=1;
break
else
flag=0;
end
end
if flag %if new adjective is found, add to route
if k~=n
visited(t)=1;
k=k+1;
route=[route t];
else
break
end
else %if not found, regress to last round
visited(route(k))=0;
route=route(1:size(route,2)-1);
knum(k)=0;
k=k-1;
end
end
oroute=route;
end
3.广度优先
clear;
clc;
N=10;
G=randgraph(N,1);
tic
route=mybfs(G);
a=toc
if size(route,2)<=1
disp('no hamiltonian path');
else
disp([route route(1)]);
endfunction [ oroute ] = hc2( G )
%HC2 Summary of this function goes here
% Detailed explanation goes here
n=size(G,1);
route=[1];
visited=zeros(n,1);%not visited=0,visited=1
k=1;%how many nodes are in current route
visited(1)=1;
a.node=1;
a.count=1;
s=[];%stack records node that needs to be visited
knum=zeros(n,1);%kth node’s last adjective
flag2=1; %flag of regress, 0 means last loop was regression
while k>=1
t=0;flag=0;a=[];
cnode=0;
if flag2 % if last loop was regression, dont search for new nodes, if not, search
for i=2:n
count=0;
if G(route(k),i)>0&&~visited(i) %find all possible nodes of current node
for j=2:n
if j~=route(k)&&G(i,j)>0&&~visited(j)
count=count+1; %out-degree
end
end
flag=1;
knum(k)=knum(k)+1;
cnode=cnode+1;
a(cnode).node=i; %record it’s name
a(cnode).count=count; %recode out-degree
end
end
end
flag2=1;
if flag %if new nodes are found
[f,seq]=sort([a.count]);%use out-degree as measurement,the larger the better
for i=1:size(f,2)
t.node=a(seq(i)).node;
t.count=f(i);
s=[t s]; %push to stack
end
end
if flag %if new nodes are found
if k~=n %if current route length is not equal to n,pop stack and add to current route
t=s(1).node;
s=s(2:size(s,2));
visited(t)=1;
knum(k)=knum(k)-1;
k=k+1;
route=[route t];
elseif k==n&&G(1,route(k))>0 %if current route’s length is n and it’s a hamilton route
oroute=route;
break;
end
else %if no new nodes are found
if k==n&&G(1,route(k))>0 %if current route’s length is n and it’s a hamilton route
oroute=route;
break;
elseif k==1 %if current route’s length is 1, which means no hamilton are found
oroute=[1];
break;
end
if knum(k-1)>0 %any nodes left for current layer, pop and replace current node
t=s(1).node;
s=s(2:size(s,2));
visited(t)=1; %set current node unvisited
visited(route(k))=0; %set the new node visited
route(k)=t; %replace
knum(k-1)=knum(k-1)-1;%number of nodes avaliable for this layer is reduced by 1
knum(k)=0; %set number of current node’s avaliable sons to 0
else %no nodes left for this layer,regress to upper layer
visited(route(k))=0; %set current node unvisited
route=route(1:size(route,2)-1); %discard current node
k=k-1; %length-1
flag2=0; %set regression flag to 0
end
end
end
oroute=route;
end
4.爬山法
function [ oroute ] = hc2( G )
%HC2 Summary of this function goes here
% Detailed explanation goes here
n=size(G,1);
route=[1];
visited=zeros(n,1);%not visited=0,visited=1
k=1;%how many nodes are in current route
visited(1)=1;
a.node=1;
a.count=1;
s=[];%stack records node that needs to be visited
knum=zeros(n,1);%kth node's last adjective
flag2=1; %flag of regress, 0 means last loop was regression
while k>=1
t=0;flag=0;a=[];
cnode=0;
if flag2 % if last loop was regression, dont search for new nodes, if not, search
for i=2:n
count=0;
if G(route(k),i)>0&&~visited(i) %find all possible nodes of current node
for j=2:n
if j~=route(k)&&G(i,j)>0&&~visited(j)
count=count+1; %out-degree
end
end
flag=1;
knum(k)=knum(k)+1;
cnode=cnode+1;
a(cnode).node=i; %record it's name
a(cnode).count=count; %recode out-degree
end
end
end
flag2=1;
if flag %if new nodes are found
[f,seq]=sort([a.count]);%use out-degree as measurement,the larger the better
for i=1:size(f,2)
t.node=a(seq(i)).node;
t.count=f(i);
s=[t s]; %push to stack
end
end
if flag %if new nodes are found
if k~=n %if current route length is not equal to n,pop stack and add to current route
t=s(1).node;
s=s(2:size(s,2));
visited(t)=1;
knum(k)=knum(k)-1;
k=k+1;
route=[route t];
elseif k==n&&G(1,route(k))>0 %if current route's length is n and it's a hamilton route
oroute=route;
break;
end
else %if no new nodes are found
if k==n&&G(1,route(k))>0 %if current route's length is n and it's a hamilton route
oroute=route;
break;
elseif k==1 %if current route's length is 1, which means no hamilton are found
oroute=[1];
break;
end
if knum(k-1)>0 %any nodes left for current layer, pop and replace current node
t=s(1).node;
s=s(2:size(s,2));
visited(t)=1; %set current node unvisited
visited(route(k))=0; %set the new node visited
route(k)=t; %replace
knum(k-1)=knum(k-1)-1;%number of nodes avaliable for this layer is reduced by 1
knum(k)=0; %set number of current node's avaliable sons to 0
else %no nodes left for this layer,regress to upper layer
visited(route(k))=0; %set current node unvisited
route=route(1:size(route,2)-1); %discard current node
k=k-1; %length-1
flag2=0; %set regression flag to 0
end
end
end
oroute=route;
end
clear;
clc;
warning off all;
N=10;
G=randgraph(N,1);
tic
route=hc2(G);
a=toc
if size(route,2)<=1
disp('no hamilton route');
else
disp([route route(1)]);
end
5.分支界限
function [ edge ] = randgraph( N,type )
%UNTITLED5 Summary of this function goes here
% Detailed explanation goes here
MaxWeight=20;
edge=zeros(N);
for i=1:N
for j=i+1:N
r=rand;
if type==1
if r>=0.5
if type==1
edge(i,j)=1;
edge(j,i)=1;
end
end
elseif type==2
edge(i,j)=MaxWeight*rand;
edge(j,i)=edge(i,j);
end
end
end
end
function [ min,heap ] = heapmin( oheap )
%UNTITLED2 Summary of this function goes here
% Detailed explanation goes here
%simply return the min value from a minmium heap,and update rest of theheap
n=size(oheap,2);
min=oheap(1);
if n==1
heap=[];
return
end
heap=[oheap(n) oheap(:,2:n-1)];
c=1;
while c*2<=n-1
left=heap(c*2).weight;
if (2*c+1<=n-1)
right=heap(c*2+1).weight;
if left<=right&&left>heap(c).weight
t=heap(c*2);
heap(c*2)=heap(c);
heap(c)=t;
c=c*2;
elseif right<left&&right>heap(c).weight
t=heap(c*2+1);
heap(c*2+1)=heap(c);
heap(c)=t;
c=c*2+1;
else
break;
end
else
if left<heap(c).weight
t=heap(c*2);
heap(c*2)=heap(c);
heap(c)=t;
c=c*2;
else
break;
end
end
end
end
function [ output ] = heapinsert( oheap,new )
%UNTITLED Summary of this function goes here
% Detailed explanation goes here
%a simple insert to a minimum heap
heap=[oheap new];
n=size(heap,2);
c=n;
while c~=1&&heap(c).weight<heap(floor(c/2)).weight
t=heap(c);
heap(c)=heap(floor(c/2));
heap(floor(c/2))=t;
c=floor(c/2);
end
output=heap;
end
function [ oroute,oweight ] = BnB( G )
%UNTITLED Summary of this function goes here
% Detailed explanation goes here
oroute=[];
oweight=[];
n=size(G,1);
%queue records node need to be visited of the same layer
myqueue.node=1;
myqueue.route=1; %route to this node
myqueue.weight=0; %weight of this route
myqueue.visited=zeros(n,1);
myqueue.visited(1)=1;
heap=myqueue;%minimum heap
minw=-1;%minimum weight, -1 stands for +infinity
while size(heap,1)>0
[current heap]=heapmin(heap); %get the min node in heap
flag=0; %flag of new nodes, 1 means new nodes are added
if minw>0&¤t.weight>minw %if greater than minimum,cut
continue;
else
for i=1:n
if G(current.node,i)>0&&(~current.visited(i)) %find new son nodes
t.node=i;
t.route=[current.route i];
t.visited=current.visited;
t.visited(i)=1;
t.weight=current.weight+G(current.node,i);
heap=heapinsert(heap,t);%add new nodes to heap
flag=1;
end
end
end
if ~flag %if no new nodes are added
if size(current.route,2)==n; %if length of current route equals n
if G(1,current.route(n))>0 %if current route is hamilton route, update output and minimun weight
oroute=current.route;
minw=current.weight;
oweight=current.weight;
end
end
end
end
end
clear;
clc;
N=10;
G=randgraph(N,2);
tic
[route weight]=BnB(G);%branch& bounding used in best first search
a=toc
if size(route,2)<=1
disp('no hamiltonian path');
else
disp([route route(1)]);
disp(['weight ' num2str(weight)]);
end
