function myTSP1
%城市數(shù)目
N=5;
%城市坐標(biāo)及城市間距離
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];
for i=1:1:N
    for j=1:1:N
        d(i,j)=sqrt((cityx(i)-cityx(j))^2+(cityy(i)-cityy(j))^2);
    end
end
%網(wǎng)絡(luò)參數(shù)
A=500;
B=500;
C=1000;
D=500;
u0=0.02;
tao=1;
lamda=0.0001;
%求得一個合法解
%統(tǒng)計每次求得一個合法解要經(jīng)過多少次非法解
total=0;
%結(jié)束標(biāo)志
toend=0;
time=clock;
display(['current time is  ',num2str(time(1,4:6))])
while toend==0
    total=total+1
    %換位陣及初始化
    V=rand(N,N);
    U=atanh(2*V-1)*u0;
    %狀態(tài)更新
    for renew=1:1:1000
        %同步更新
        for ux=1:1:N
            for ui=1:1:N
                m1=0;
                m2=0;
                m3=0;
                m4=0;
                %求導(dǎo)公式第一項(xiàng)
                for j=1:1:N
                    if j~=ui
                        m1=m1+V(ux,j);
                    end                   
                end
                m1=-A*m1;
                %求導(dǎo)公式第二項(xiàng)
                for y=1:1:N
                    if y~=ux
                        m2=m2+V(y,ui);
                    end
                end
                m2=-B*m2;
                %求導(dǎo)公式第三項(xiàng)
                for x=1:1:N
                    for j=1:1:N
                        m3=m3+V(x,j);
                    end
                end
                m3=-C*(m3-N);
                %求導(dǎo)公式第四項(xiàng)
                for y=1:1:N
                    if y~=ux
                        if ui==1
                            m4=m4+d(ux,y)*(V(y,ui+1)+V(y,N));
                        elseif ui==N
                            m4=m4+d(ux,y)*(V(y,ui-1)+V(y,1));
                        else
                            m4=m4+d(ux,y)*(V(y,ui+1)+V(y,ui-1));
                        end
                    end  
                end
                m4=-D*m4;
                Udao(ux,ui)=-U(ux,ui)+m1+m2+m3+m4;
            end
        end
        %導(dǎo)數(shù)及狀態(tài)更新                
        U=U+lamda*Udao;
        V=(1+tanh(U/u0))/2;
        for ux=1:1:N
            for ui=1:1:N
                if V(ux,ui)<0.3
                    V(ux,ui)=0;
                end
                if V(ux,ui)>0.7
                    V(ux,ui)=1;
                end 
            end
        end
    end
    V;
    %判斷是否為合法解
    %換位陣全局約束，要求總共有N個1
    test1=0;
    for ux=1:1:N
        for ui=1:1:N
            test1=test1+V(ux,ui);
        end
    end  
    %城市行約束，每行不多于一個1 
    test2=0;
    for x=1:1:N
        for i=1:1:N-1
            for j=i+1:1:N
                test2=test2+V(x,i)*V(x,j);
            end
        end
    end    
    %城市列約束，每列不多于一個1   
    test3=0;
    for i=1:1:N
        for x=1:1:N-1
            for y=x+1:1:N
                test3=test3+V(x,i)*V(y,i);
            end
        end
    end
    %當(dāng)為合法解時，跳出循環(huán)
    if test1==N && test2==0 && test3==0
        toend = 1;
    else
        toend=0;
    end
end
time=clock;
display(['end time is  ',num2str(time(1,4:6))])
V
total
for j=1:1:N
    for i=1:1:N
        if V(i,j)==1
            cityx_final(j)=cityx(i);
            cityy_final(j)=cityy(i);
        end
    end
end
cityx_final(N+1)=cityx_final(1);
cityy_final(N+1)=cityy_final(1);
cityx_final;
cityy_final;
%總路程長度
td=0;
for i=1:1:N-1
    td=td+sqrt((cityx_final(i)-cityx_final(i+1))^2+(cityy_final(i)-cityy_final(i+1))^2);
end
td=td+sqrt((cityx_final(N)-cityx_final(1))^2+(cityy_final(N)-cityy_final(1))^2);
td
plot(cityx_final,cityy_final,'o-');