hopfield_TSP程序

人工神经网络实验用CHNN算法求解TSP问题一．问题描述利用连续型Hopfield反馈网络求解10城市的旅行商（TSP）问题。其中10个城市的坐标给定如下：1(0.4000,0.4439),2(0.2439,0.1463),3(0.1707,0.2293),4(0.2293,0.7610),5(0.5171,0.9414),6(0.8732,0.6536),7(0.6878,0.5219),8(0.8488,0.3609),9(0.6683,0.2536),10(0.6195,0.2634),基本网络参数为：二．算法实现1．CHNN算法应用CHNN网络解决优化问题一般需要以下步骤：(1.)对于特定的问题，要选择一种合适的表示方法，使得神经网络的输出与问题的解相对应。(2.)构造网络的能量函数，使其最小值对应于问题的最佳解。(3.)将能量函数与CHNN算法标准形式相比较，推出神经网络权值与偏流表达式。(4.)推出网络状态更新公式，并利用更新公式迭代求问题的最优解。2．TSP问题为使用CHNN网络进行TSP问题的求解，根据上述步骤，可将问题转化为：(1.)对N个城市的TSP问题，用一个N*N的换位阵描述旅行路线，换位阵中每行每列有且只有一个元素为1，其余全为0。为1的元素其横坐标x表示城市名，纵坐标i表示该城市在访问路线中的位臵。(2.)网络的能量函数由四部分组成，分别用来保证换位阵的合法性以及最终路线长度的最短。(3.)将能量函数与标准形式相比较，得到网络权值与偏流表达式为(4.)从而，网络更新公式为：3．程序设计根据上述推导在MATLAB中设计CHNN网络求解TSP问题的程序（程序代码见附页）。(1.)程序说明本程序中有以下两点需要说明。迭代结束条件：理论上来说，当网络的能量函数不再减小时网络达到最优状态，但在实际中如果用能量函数的变化来判断程序的结束存在潜在的问题（如计算能量函数的复杂性以及误差导致判断不准确等），因此，本实验利用迭代次数控制程序结束，当迭代至1000次时，一次运行结束。程序输出规则：由于不能保证每次迭代结束所到的解都是合法解，而当城市个数较多时人为检查合法性又非常的不方便，因此每次迭结束后在程序中检查该次解的合法性，若为合法解，则输出该解，程序结束；否则，再次求解。参数调整：在实验中发现，当网络参数取为最初给定的值A=B=D500,C=200,u0=0.02时，几乎得不到合法解，观察每次迭代结束后的解，发现大部分下只有8个每行每列有且只有一个1的情况，另外还有两列全部为0。这说明在能量函数中保证有N个1的合法性所占的比重相对较小，也就是参数C相对于A、B、D来说较小，因此，将基本参数C调整为1000，其余不变。(2.)程序流程i.初始化：城市个数、城市坐标、网络参数ii.用随机数初始化换位阵及状态阵iii.对状态阵及换位阵，进行1000步同步更新，得最终换位阵的解Viv.判断所得V的合法性，若为合法解，给出访问次序，旅行路线图及路线总长度，程序结束；否则，转到第ii步。三．实验结果1．基本结果在城市个数取为10，网络的基本参数取为A=B=D=500,C=1000,u0=0.02,lamda=0.0001时运行程序并统计实验结果，得：（图见下页）运行次数200合法解次数29最优解次数1最优解（路线总长度）2.6907次优解次数1次优解（路线总长度）2.7693较优解（路线总长度）2.7782较优解（路线总长度）2.8352平均一次运行所需时间(s)0.88134图1最优路线（2.6907）图2最优解换位阵图3次优路线(2.7693)图4次优换位阵图5较优路线(2.7782)图6较优换位阵52．参数影响(1.)运行时间估计在城市数目N及更新步长lamda固定的情况下，每求解一次V所用的时间是固定的，因此，比较每次出现合法解所用的时间可通过比较循环次数进行。在下面参数影响的讨论中，均通过循环次数比较相对时间长短。(2.)权系数A、B、C权矩阵A、B、C、D的相对大小反映了对解的要求。其中A、B、C是为了保证合法解的项的权系数，A是保证每行最多一个1的权系数；B是保证每列最多一个1的权系数；C是保证共有N个1；D是保证路线总长度最短的项的权系数。当C相对于A和B较小（A=B=500，C=200）时，实验很难出现合法解，多数解都有两列全为0，程序往往陷入死循环。这说明，解的合法性的第三项没有得到足够的重视。因此，逐渐加大C并观察实验结果，当C为500时，上述情况仍没有明显改善；当C取为1000时，合法解出现频率明显提高（200次实验中，平均每6.7次出现一次合法解），其中也出现了最优解（见1中的实验结果）；当C取为2000是，平均每6次出现一次全法解，其中同样出现一次最优解。总结，C较小不能保证解的合法性，C较大时出现合法解的频率明显提高，但同时C较大时路线最短项的权系数D相对较小，因此，出现最优解的频率将有所下降。(3.)权系数D权系数D反映了路线长度在能量函数中所占的比重。当D取为200时，平均每1.5次出现一次合法解，但路线长度非常大，一般在4.0左右，几乎不能出现最优解；当D取为500时，平均每6.7次出现一次合法解，其中也出现了最优解（见1中的实验结果）；当D取为600时，出现合法解的频率有所下降；当D取为700时，151次运行中出现一次较优解；当D取为1000时，6程序几乎陷于死循环，出现合法解的几率极低。总结，D较小时，相对更强调解的合法性，因此出现合法解的频率较大，但路线长度很大；D较大时，出现合法解的频率有所降低，但路线长度明显变小，出现最优解的可能性相对增加；而当D过大时，由于过度强调路线长度，很难出现合法解，因此程序易冻结。(4.)步长lamda当lamda为0.0001时运行结果如1中所述；当lamda取为0.001时，平均每2.5次出现一次合法解。可见，lamda较大时，状态矩阵变化较大，会提高出现合法解的频率；但lamda过大时状态矩阵会由于变化剧烈而难以出现合法解；lamda较小时会导致更新速度过慢甚至冻结(5.)初值u0取为0.02时运行结果如1中所述；当取为0.005时，平均每3.6次出现一次最优解；取为0.3时，平均每41次出现一次全法解。可见，较小，激励函数趋近于离散值，缩短出现寻优时间，但不易出现最优解；较大，激励函数过于平坦，不利于收敛。3．改变城市数目下面分别给出城市数目为5和11，网络参数不变时的实验结果。由实验结果可知，城市数目下降，寻优时间缩短，得到合法解和最优解的频率明显增加。另外，由于网络参数对实验的影响同前面类似，此处不再赘述。(1.)城市数目为11（第11个城市的坐标为(0.9125,0.9568)）平均每7.5次出现一次合法解，其中一次较优解路线长度为3.1382，图形如下：7图7十一城市TSP问题较优解(2.)城市数目为5（取前5个城市的坐标）平均每5次出现一次合法解，35次实验中出现3次最优解。最优解为1.8324，其中还多次出现较优解1.8904，图形如下：图8五城市TSP问题的最优解图9五城市TSP问题较优解8四．附页（程序代码）functionmyTSP1%城市数目N=10;%5%11%城市坐标及城市间距离cityx=[0.4,0.2439,0.1707,0.2293,0.5171,0.8732,0.6878,0.8488,0.6683,0.6195,0.9125];cityy=[0.4439,0.1463,0.2293,0.761,0.9414,0.6536,0.5219,0.3609,0.2536,0.2634,0.9568];fori=1:1:Nforj=1:1:Nd(i,j)=sqrt((cityx(i)-cityx(j))^2+(cityy(i)-cityy(j))^2);endend%网络参数A=500;B=500;C=1000;D=500;u0=0.02;tao=1;lamda=0.0001;%求得一个合法解%统计每次求得一个合法解要经过多少次非法解total=0;%结束标志toend=0;time=clock;display(['currenttimeis',num2str(time(1,4:6))])whiletoend==0total=total+1%换位阵及初始化V=rand(N,N);U=atanh(2*V-1)*u0;%状态更新forrenew=1:1:1000%同步更新forux=1:1:Nforui=1:1:Nm1=0;m2=0;m3=0;m4=0;%求导公式第一项forj=1:1:Nifj~=uim1=m1+V(ux,j);endendm1=-A*m1;%求导公式第二项fory=1:1:Nify~=uxm2=m2+V(y,ui);endendm2=-B*m2;%求导公式第三项forx=1:1:Nforj=1:1:Nm3=m3+V(x,j);endendm3=-C*(m3-N);%求导公式第四项fory=1:1:Nify~=uxifui==1m4=m4+d(ux,y)*(V(y,ui+1)+V(y,N));elseifui==Nm4=m4+d(ux,y)*(V(y,ui-1)+V(y,1));elsem4=m4+d(ux,y)*(V(y,ui+1)+V(y,ui-1));endendendm4=-D*m4;Udao(ux,ui)=-U(ux,ui)+m1+m2+m3+m4;endend%导数及状态更新U=U+lamda*Udao;V=(1+tanh(U/u0))/2;forux=1:1:Nforui=1:1:NifV(ux,ui)0.3V(ux,ui)=0;endifV(ux,ui)0.7V(ux,ui)=1;endendendendV;%判断是否为合法解%换位阵全局约束，要求总共有N个1test1=0;forux=1:1:Nforui=1:1:Ntest1=test1+V(ux,ui);endend%城市行约束，每行不多于一个1test2=0;forx=1:1:Nfori=1:1:N-1forj=i+1:1:Ntest2=test2+V(x,i)*V(x,j);endendend%城市列约束，每列不多于一个1test3=0;fori=1:1:Nforx=1:1:N-1fory=x+1:1:Ntest3=test3+V(x,i)*V(y,i);endendend%当为合法解时，跳出循环iftest1==N&&test2==0&&test3==0toend=1;elsetoend=0;endendtime=clock;display(['endtimeis',num2str(time(1,4:6))])Vtotal%按结果重新排列城市坐标forj=1:1:Nfori=1:1:NifV(i,j)==1cityx_final(j)=cityx(i);cityy_final(j)=cityy(i);endendendcityx_final(N+1)=cityx_final(1);cityy_final(N+1)=cityy_final(1);cityx_finalcityy_final%总路程长度td=0;fori=1:1:N-1td=td+sqrt((cityx_final(i)-cityx_final(i+1))^2+(cityy_final(i)-cityy_final(i+1))^2);endtd=td+sqrt((cityx_final(N)-cityx_final(1))^2+(cityy_final(N)-cityy_final(1))^2);tdplot(cityx_final,cityy_final,'o-');