Traveling Salesman Evolutionary Algorithm Bugs

I'm writing an evolutionary algorithm to solve the traveling salesman problem (i.e. salesman must go to n cities and back to starting city and wants to find route with least cost) but I'm encountering a few bugs. In my algorithm I start with some random candidate routes and I use mutation and ox crossover operator to get new ones. After each iteration I get rid of the worst performing half of the candidate routes and replace them with mutation of each route in the upper quarter and the children from the crossover of the second quarter.

For the mutation of the first quarter, I just pick two random cities in each route and swap them.

For the ox crossover of the second quarter, consider two parents with cut points shown by "|".

p_1=(123|4567|89) p_2=(452|1876|93)

First the segments between cut points are copied into offspring

o_1=(xxx|4567|xx) o_2=(xxx|1876|xx)

then starting from the second cut point of one parent, the cities from the other parent are copied in the same order omitting symbols already present. Reaching the end of the string, we continue from the first place of the string. Therefore the offspring are

o_1=(218|4567|93) o_2=(345|1876|92)

I have shown my code below where there are 17 cities and I iterate my algorithm 250 times. It works OK but sometimes it outputs an optimal route that is an invalid route. i.e. a route that visits one city twice and leaves out some cities. This only happens in about 2 out of 10 times but it is very frusterating. If anyone could tell me what I should change to fix this problem, it would be greatly appreciated

    %let n be the number of cities that must be visited
  n=17;

%let numCS be the initial number of candidate solutions
numCS=128;

%Will first initialise the population of random candidate solutions
for i=1:numCS
    tempCandSol(i,1:n)=randperm(n);
end

%But a tour must end at the city it started from so that city must be added
candSol=[tempCandSol,tempCandSol(1:numCS,1)];

%Now we must evaluate the respective fitnesses of the candidate solutions

%must load the data set that tells us the traveling costs
load tspadata1.txt
%Will call this the costMatrix
costMatrix=tspadata1;

%Will now evaluate each candidate solution in terms of its total cost.
%Will create fitness numCSx1 matrix where each row is the fitness of each
%candiate solution.
sum=0;
for j=1:numCS
    for k=1:n
    fitness(j,1)=costMatrix(candSol(j,k),candSol(j,k+1))+sum;
    sum=fitness(j,1);
    end
    sum=0;
end

%will now sort the candidate solutions wrt their fitness in acsending order
%order

%put fitness in candSol matrix to make candSolFitness matrix
candSolFitness=[candSol,fitness];

%now sort in ascending order wrt last fitness column
[sortedCol,sorter]=sort(candSolFitness(:,n+2),'ascend');

%Will call the sorted matrix rankedSol
rankedSol=candSolFitness(sorter,:);

%now we will choose to top performing half to breed and mutate and will let
%the bottom performing half die off.  Will mutate the top performing quater
%but keep the the original solutions that were mutated too.  Will use crossover
%on the bottom quarter and keep both parents and offspring so that in the
%end we will have kept the population constant.

for e=1:250
%Will now mutate the first quarter using reciprocal exchange/swap operator and put the mutated versions of the
%first quarter in the thrid quarter while keeping the originals in the
%first quarter.

%Will first copy top quarter of rankedSol into third quarter (i.e.overwritting it)
for r=1:numCS/4
    rankedSol((numCS/2)+r,:)=rankedSol(r,:);
end

%now mutate the third quarter
 for m=numCS/2+1:(3/4)*numCS
     %pick two random cities to swap
     city1Pos=randi(17);
     city2Pos=randi(17);
     %cannot pick first city or else invalid route
     if city1Pos==1
          city1Pos=2;
     end
     if city2Pos==1
          city2Pos=2;
     end
    %swap city1Pos entry with city2Pos entry
    temp_mutated1=rankedSol(m,city1Pos);
    temp_mutated2=rankedSol(m,city2Pos);
    rankedSol(m,city2Pos)=temp_mutated1;
    rankedSol(m,city1Pos)=temp_mutated2; 
 end

   %will now use ox crossover operator on second quarter and put the offspring in the fourth quarter thus overwritting the original fourth quarter
   %must first define cut points
   %using n-1 instead of n because I don't want first city to be changed or
   %else would be invalide route.
   cutPoint1=floor((1/3)*(n-1));
   cutPoint2=ceil((2/3)*(n-1));

   %Initialize offspring as zero vectors
   %offspring1=zeros(1,n);
   %offspring2=zeros(1,n);

   for y=(numCS/4)+1:((1/2)*((numCS/2)-(numCS/4))+(numCS/4))
       %must preserve first city and last city to ensure valid route
       offspring1(y-(numCS/4),1)=rankedSol(y,1);
       offspring1(y-(numCS/4),n+1)=rankedSol(y,n+1);
       %offspring2(y-(numCS/4),1)=rankedSol(y+1,1);
       %offspring2(y-(numCS/4),n+1)=rankedSol(y+1,n+1);
       %copy cities between cutpoints to offspring
       offspring1(y-(numCS/4),cutPoint1:cutPoint2)=rankedSol(y,cutPoint1:cutPoint2);
       %offspring2(1,cutPoint1:cutPoint2)=rankedSol(y+1,cutPoint1:cutPoint2);
       %Now construct vector of cities outside cutpoints for both parents
   end

   for y=((1/2)*((numCS/2)-(numCS/4))+(numCS/4))+1:numCS/2

       %must preserve first city and last city to ensure valid route
       offspring2(y-(((1/2)*((numCS/2)-(numCS/4))+(numCS/4))),1)=rankedSol(y,1);
       offspring2(y-(((1/2)*((numCS/2)-(numCS/4))+(numCS/4))),n+1)=rankedSol(y,n+1);
       offspring2(y-(((1/2)*((numCS/2)-(numCS/4))+(numCS/4))),cutPoint1:cutPoint2)=rankedSol(y,cutPoint1:cutPoint2);
   end

     for y=(numCS/4)+1:((1/2)*((numCS/2)-(numCS/4))+(numCS/4))
       for s=1:n-(cutPoint2+1)+1
          newVect1(y-(numCS/4),s)=rankedSol(y,cutPoint2+s);
       end
       %newVect1(1,1:n-(cutPoint2+1))=rankedSol(cutPoint2+1:n);
       for w=n-(cutPoint2+1)+2:n-(cutPoint2-cutPoint1)-1
          newVect1(y-(numCS/4),w)=rankedSol(y,w-(cutPoint2-cutPoint1));
       end
       %newVect1(1,(n-(cutPoint2+1))+1:((n-(cutPoint2+1))+1)+((cutPoint1-1)-2))=rankedSol(y,2:cutPoint1-1);
     end

     for y=((1/2)*((numCS/2)-(numCS/4))+(numCS/4))+1:numCS/2
       for a=1:n-(cutPoint2+1)+1
          newVect2(y-(((1/2)*((numCS/2)-(numCS/4))+(numCS/4))),a)=rankedSol(y,cutPoint2+a);
       end
       %newVect2(1,1:(n-(cutPoint2+1)))=rankedSol(y+1,cutPoint2+1:n);
         for b=n-(cutPoint2+1)+2:n-(cutPoint2-cutPoint1)-1
          newVect2(y-(((1/2)*((numCS/2)-(numCS/4))+(numCS/4))),b)=rankedSol(y,b-(cutPoint2-cutPoint1));
       end
     end

newNewVect1=[newVect1,offspring1(:,cutPoint1:cutPoint2)];
newNewVect2=[newVect2,offspring2(:,cutPoint1:cutPoint2)];

%start with first offspring

       for y=1:numCS/8
          for h=cutPoint2+1:n
              %make sure that no city is repeated
              if isRepeated(newNewVect2(y,h-cutPoint2),offspring1,y)==0
                  offspring1(y,h)=newNewVect2(y,h-cutPoint2);
              else
                  for u=h-cutPoint2:n
                      if isRepeated(newNewVect2(y,u),offspring1,y)==0
                          break
                      end
                  end
                  offspring1(y,h)=newNewVect2(y,u);
                  %offspring1(y,h)=newNewVect2(y,nextNewIndex(newNewVect2,h-cutPoint2,offspring1,y));
             end

          end
       end

       for y=1:numCS/8
       for t=2:cutPoint1-1
           if isRepeated(newNewVect2(y,(n-cutPoint2)+t-1),offspring1,y)==0
               offspring1(y,t)=newNewVect2(y,(n-cutPoint2)+t-1);
           else
               for u=(n-cutPoint2)+t-1:n
                      if isRepeated(newNewVect2(y,u),offspring1,y)==0
                          break
                      end
                  end
                  offspring1(y,t)=newNewVect2(y,u);
               %offspring1(y,t)=newNewVect2(y,nextNewIndex(newNewVect2,(n-cutPoint2)+t-1,offspring1,y));
           end

       end
       end
      %Now for second offspring

      for y=1:numCS/8

       for q=cutPoint2+1:n
       if isRepeated(newNewVect1(y,q-cutPoint2),offspring2,y)==0
           offspring2(y,q)=newNewVect1(y,q-cutPoint2);
       else
           for u=q-cutPoint2:n
                      if isRepeated(newNewVect1(y,u),offspring2,y)==0
                          break
                      end
                  end
                  offspring2(y,q)=newNewVect1(y,u);
           %offspring2(y,q)=newNewVect1(y,nextNewIndex(newNewVect1,q-cutPoint2,offspring2,y));
       end

       end
      end

      for y=1:numCS/8
       for c=2:cutPoint1-1
           if isRepeated(newNewVect1(y,(n-cutPoint2)+c-1),offspring2,y)==0
               offspring2(y,c)=newNewVect1(y,(n-cutPoint2)+c-1);
           else
               for u=(n-cutPoint2)+c-1:n
                      if isRepeated(newNewVect1(y,u),offspring2,y)==0
                          break
                      end
                  end
                  offspring2(y,c)=newNewVect1(y,u);
               %offspring2(y,c)=newNewVect1(y,nextNewIndex(newNewVect1,(n-cutPoint2)+c-1,offspring2,y));
           end

       end
      end
      %now that we have the offspring we can overwrite them into the fourth
      %quarter.

      rankedSol((3/4)*numCS+1:(7/8)*numCS,1:n+1)=offspring1;
      rankedSol((7/8)*numCS+1:numCS,1:n+1)=offspring2;

      %fitness column is out-of-date so must compute fitness of the new
      %solutions
      %take away old fitnesses
      unrankedSol=rankedSol(:,1:n+1);

      %evaluate new actual fitness
      sum2=0;
      for j=1:numCS
          for k=1:n
              fitness(j,1)=costMatrix(unrankedSol(j,k),unrankedSol(j,k+1))+sum2;
              sum2=fitness(j,1);
          end
          sum2=0;
      end

      %will now sort the candidate solutions wrt their fitness in acsending order
      %order

      %put fitness in candSol matrix to make candSolFitness matrix
      newCandSolFitness=[unrankedSol,fitness];

      %now sort in ascending order wrt last fitness column
      [sortedCol,sorter]=sort(newCandSolFitness(:,n+2),'ascend');

      %Will call the sorted matrix rankedSol
      newRankedSol=newCandSolFitness(sorter,:);

      rankedSol=newRankedSol;

end

disp('the best route that this algorithm was able to find is:');
disp(rankedSol(1,1:n+1));
disp('the cost associated with this route is:')
disp(rankedSol(1,n+2));

and it uses the following function to determine if a city has been repeated

function [repeated]=isRepeated(element,vectToCheck,row)
[N,M]=size(vectToCheck);

d=0;
for v=1:M
    if element==vectToCheck(row,v);
        d=1;
    end
end

if d==0
    repeated=0;
end

if d==1
    repeated=1;
end

Again, if anyone can tell me why some routes have repeated cities, it would be greatly appreciated

Answer 1

[Vito]

"The task of the direct-sales representative" is among combinatorial tasks. The finite decision represents a cycle. It is the task also belongs to the class of the transport task or the task about optimal resource allocation. But it it is impossible to solve a simplex a method.

The algorithm of its decision consists in search of all cycles.

Error in what not all cycles are considered.

The code more low visually shows the decision of the classical transport task. Such task can be solved and a combinatorial method and a simplex.

In the task we have three points of departure and four assignments. The task consists in a choice of a route with the least expenses. GUI the interface visually shows a choice of cycles in the table. Demonstration is presented in the form of slow animation.

On the termination of iterations it is possible to look at the graph of ways.

Probably, it will help to solve a problem.

Comments in Russian. If there will be questions, will be glad to help. ---------------

Create a file - linearAlgebraGIn and copy there the following code. ------------------

 function  linearAlgebraGIn(arg)
  global  Mtemp CC

C = {['ba'] [5] [10] [20] [15]; [10] [8 Inf] [3 Inf] [5 Inf] [2 Inf]; [15] [4 Inf] [1 Inf] [6 Inf] [7 Inf]; [25] [1 Inf] [9 Inf] [4 Inf] [3 Inf]};

% матрица результатов

Mtemp =C;

if nargin == 0;
   initialize
elseif arg == 0
  action(Mtemp)
elseif arg == 1
  ResultGraph(CC)  
%elseif arg < 0
%   setmode(arg)
else
   initialize(arg);
end

end

function  initialize % инициализация исходных параметров

global  Mtemp

%---------------------------
% графический интерфейс
figNumber=figure( ...
        'Visible','on', ...
        'NumberTitle','of', ...
        'Name','Linear - programming, metod of optimization');
    axes( ...
        'Units','normalized', ...
        'Position',[-0.5 0.4 1.7 0.57])
colormap(Cool);
%set(figNumber, 'Renderer', 'OpenGL');
%set(figNumber,'DoubleBuffer','on')
 % Set up the MiniCommand Window
    top=0.30;
    left=0.05;
    right=0.70;
    bottom=0.05;
    labelHt=0.05;
    spacing=0.005;
    % First, the MiniCommand Window frame
    frmBorder=0.02;
    frmPos=[left-frmBorder bottom-frmBorder ...
        (right-left)+2*frmBorder (top-bottom)+2*frmBorder];
    uicontrol( ...
        'Style','frame', ...
        'Units','normalized', ...
        'Position',frmPos, ...
        'BackgroundColor',[0.50 0.50 0.50]);
  % Then the text label
    labelPos=[left top-labelHt (right-left) labelHt];
    uicontrol( ...
        'Style','text', ...
        'Units','normalized', ...
        'Position',labelPos, ...
        'BackgroundColor',[0.50 0.50 0.50], ...
        'ForegroundColor',[1 1 1], ...
        'String','Info Window');
 % Then the editable text field
   fval=SimplexResult(Mtemp);

     strcat('Simplex result =',num2str(fval))

      mcwPos=[left bottom (right-left) top-bottom-labelHt-spacing];
      mcwHndl=uicontrol( ...
          'Style','edit', ...
          'HorizontalAlignment','left', ...
          'Units','normalized', ...
          'Max',20, ...
          'BackgroundColor',[1 1 1], ...
          'Position',mcwPos, ...
          'String',strcat('Simplex result =',num2str(fval)));    
   %====================================
      % Information for all buttons
      labelColor=[0.8 0.8 0.8];
      top=0.95;
      bottom=0.05;
      left=0.75;
      yInitLabelPos=0.90;
      left=0.75;
      labelWid=0.20;
      labelHt=0.05;
      btnWid=0.20;
      btnHt=0.05;
      % Spacing between the label and the button for the same command
      btnOffset=0.003;
      % Spacing between the button and the next command's label
      spacing=0.05;
   %====================================
      % The CONSOLE frame
      frmBorder=0.02;
      yPos=0.05-frmBorder;
      frmPos=[left-frmBorder yPos btnWid+2*frmBorder 0.9+2*frmBorder];
      h=uicontrol( ...
          'Style','frame', ...
          'Units','normalized', ...
          'Position',frmPos, ...
          'BackgroundColor',[0.50 0.50 0.50]);   

   % The PLOT TYPE command popup button
      btnNumber=1;
      yLabelPos=top-(btnNumber-1)*(btnHt+labelHt+spacing);
      labelStr=' Optimization type';
      labelList=' Combinatoric classic| Combinatoric optimization';
      cmdList=str2mat( ...
          ' 1 ',' 2 ');

    % Generic label information
      labelPos=[left yLabelPos-labelHt labelWid labelHt];
      uicontrol( ...
          'Style','text', ...
          'Units','normalized', ...
          'Position',labelPos, ...
          'BackgroundColor',labelColor, ...
          'HorizontalAlignment','left', ...
          'String',labelStr);  

    % Generic popup button information
      btnPos=[left yLabelPos-labelHt-btnHt-btnOffset btnWid btnHt];
      hndl1=uicontrol( ...
          'Style','popup', ...
          'Units','normalized', ...
          'Position',btnPos, ...
          'String',labelList, ...
          'UserData',cmdList);

     %==================================== 
   uicontrol('Style', 'text','BackgroundColor',[1 1 0],  'String', 'Number Iteration', 'Units','normalized',...
      'Position', [left bottom+4*btnHt+spacing btnWid 0.5*btnHt]);

h_angle = uicontrol('Style', 'edit','FontWeight','bold','FontSize',[10],'String', '0', 'Units','normalized',...
    'Position', [left bottom+2*btnHt+spacing btnWid 2*btnHt], 'Tag','PlotText');

callbackStr='set(gca,''Userdata'',-1)';
% The close button.
%callbackStr='ResultGraph @Mtemp';
stopHndl= uicontrol( ...
      'Style','pushbutton', ...
        'Units','normalized', ...
        'Position',[left 7*bottom btnWid 2*btnHt], ...
        'String','Show resylt graph', ...
        'Callback', 'linearAlgebraGIn(1)');   
%set(stopHndl,'Enable','off');

%---------------------------
% исходная матрица
%C = {['ba'] [5] [10] [20] [15]; [10] [8 0] [3 0] [5 0] [2 0];[15] [4 0] [1 0] [6 0] [7 0];[25] [1 0] [9 0] [4 0] [3 0]};

%global Mtemp

C=Mtemp;

for i = 2 : 4  % цикл по пунктам отправления 
     for j =2 : 5  % цикл по пунктам назначения

           C {i, j}=num2str(C {i, j});

       end
  end

cellplot(C)
hold on; % транспонирует график
%[x,y] = ginput(1)

%-------------------
tempMin=0;
temp=0;
tempSum=0;
%action(Mtemp,C)

callbackStr='action(''Mtemp'',''C'')';

startHndl=uicontrol(...
   'style','pushbutton', ...
   'units','normalized', ...
   'position',[left bottom btnWid 2*btnHt], ...
   'string','Start', ...
   'callback', 'linearAlgebraGIn(0)');  

% Uncover the figure
        hndlList=[mcwHndl h_angle stopHndl startHndl hndl1];
set(gcf,'Visible','on', ...
            'UserData',hndlList);

end % function

function Mtemp= action(Mtemp)

CG=Mtemp; % сохраняем таблицу перед преобразованием 
% конвертируем в строку для вывода на экран

[m,n]=size (Mtemp)

for i = 2 : m  % цикл по пунктам отправления 
     for j =2 : n  % цикл по пунктам назначения

           Mtemp {i, j}=num2str(Mtemp {i, j});
           %C {i, j}=num2str(C {i, j});

       end
  end

Mtemp=CG;

drawnow
cla

cellplot(Mtemp)
hold on; % транспонирует график

hndlList=get(gcf,'UserData');

set(hndlList(1),'String','Wait...');

[m,n]=size (Mtemp)

%m = 4; % количество пунктов оправления
%n = 5; % количество пунктов назначения

Jtemp = 0; % временный индекс

for i = 2 : m  % цикл по пунктам отправления (по строкам матрицы) 

   for k = 2:n % получаем значения элементов c строки матрицы 
      V{k-1} = Mtemp{i,k}(1);
  end        

[tempMin,J] = min( [V{1:4}]); % 1. В первой строке матрицы С=(сij) отыскиваем наименьший элемент. Пусть это будет элемент с(1,j1).        

temp=min (Mtemp {i, 1},Mtemp {1, J+1} ); % минимальное значение для пункта отправки и потребления         

if Mtemp {1, J+1} > SumCol(J+1,i,Mtemp) % проверим, есть ли потребность в данном пункте?
Mtemp{i,J+1}(2)= temp; % 2. Полагаем x(1,j1)= min(a1, bj1).
 end %if

if Mtemp {i, 1} > Mtemp {1, J+1} 
    %  3. Если a1 > bj1 отыскиваем в этой же строке второй наименьший элемент с(1,j2) 
    %и полагаем x(1,j2)= min(a1 - x(1, j1), bj2).  
 while Mtemp {i, 1}~= SumRow(i,n,Mtemp)% 4. Эти шаги продолжаем, пока не удовлетворим первое условие a1 = x(1,1j) + x(1,2j)+...+ x(1, nj). 

   V{J}= Inf;
   Jtemp =J;
   [tempMin,J] = min( [V{1:4}]); % отыскиваем в этой же строке второй наименьший элемент с(1,j2)   

   if  (Mtemp {i, 1}- SumRow(i,n,Mtemp))== Mtemp {1, J+1} 

    % 5. Если на некотором шаге этого процесса окажется, что остаток от a1
    % точно равен соответствующему b(jk), то    
    % положив x(1,jk)равным этому остатку x(1, jk) = b(jk)    
    Mtemp{i,J+1}(2)= (Mtemp {i, 1}- SumRow(i,n,Mtemp));
    % 6. Полагаем следующее по величине значение переменной x(1, jk+1)=0 (выбранный нуль - дополняет набор до ациклического ) 
    V{J}= Inf;
   Jtemp =J;
   [tempMin,J] = min( [V{1:4}]); 
   % 7. В столбцы соответствующие уравнения которых (потребности пунктов потребления) b1 =  x(1i,1) + x(2i,1)+...+ x(mi, 1)
  % полностью удовлетворены нули не записываются. 
   if  SumCol(J+1,i,Mtemp) ~= Mtemp {1, J+1} % 
   Mtemp{i,J+1}(2)=0;
    end %if

    elseif Mtemp {1, J+1} > SumCol(J+1,i,Mtemp) && Mtemp{i,Jtemp+1}(2) ~= Inf % только если есть потребность в данном пункте, и была в предыдущем 

   Mtemp{i,J+1}(2)=min (Mtemp {i, 1} - Mtemp{i,Jtemp+1}(2),Mtemp {1, J+1} );  %и полагаем x(1,j2)= min(a1 - x(1, j1), bj2).  

   elseif Mtemp {1, J+1} > SumCol(J+1,i,Mtemp) && Mtemp{i,Jtemp+1}(2) == Inf  % только если есть потребность в данном пункте и не было в предыдущем.

       Mtemp{i,J+1}(2)=min (Mtemp {i, 1},Mtemp {1, J+1}- SumCol(J+1,i,Mtemp));%и полагаем x(1,j2)= min(a1 , b(i,j2) - SumCol (b(i,j2))). 

   end %if

   end % while

end %if

end % конец цикла по пунктам отправления

CG=Mtemp; % сохраняем таблицу перед преобразованием 
% конвертируем в строку для вывода на экран
for i = 2 : m  % цикл по пунктам отправления 
     for j =2 : n  % цикл по пунктам назначения

           Mtemp {i, j}=num2str(Mtemp {i, j});
           %C {i, j}=num2str(C {i, j});

       end
  end

cellplot(Mtemp)
hold on; % транспонирует график

Mtemp=CG;
%if tempMin==Inf
%text(0,0,'Error! No solve!','FontSize',18,'Color','r') 
%end
%cellplot(C)

% выводим на экран граф

%G=[
%0 0 0 0 0 0 0 0 0 0 0 0;
%0 0 0 0 0 0 0 0 0 0 0 0;
%0 0 0 0 0 0 0 0 0 0 0 0;
%0 0 0 0 0 0 1 0 0 0 0 0;
%0 0 0 0 0 1 0 0 0 0 0 0;
%0 0 0 0 1 0 1 0 0 0 0 0;
%0 0 0 1 0 1 0 0 0 0 1 0;
%0 0 0 0 0 0 0 0 0 0 0 0;
%0 0 0 0 0 0 0 0 0 0 0 0;
%0 0 0 0 0 0 0 0 0 0 0 0;
%0 0 0 0 0 0 1 0 0 0 0 1;
%0 0 0 0 0 0 0 0 0 0 1 0 ;  
%];% матрица смежности

%V=[1 4; 2 4; 3 4; 4 4;
%   1 3; 2 3; 3 3; 4 3;
%   1 2; 2 2; 3 2; 4 2;
%];% координаты вершин

[m,n] = size(CG); % размер исходного массива
% Вначале преобразуем матрицу, удалив из нее информационные строки 

CT= cell(m-1,n-1);

for i = 1:m-1
    for j = 1:n-1
    CT{i,j}=CG{i+1,j+1};    
    end
end

[m,n] = size(CT); % размер исходного массива
% Теперь вектор (аналог всов вершин)

%vwieght =cell(1,m*n); % можно не объявлять

%vwieght =[1:m*n];% можно и так

c=1; % счетчик

for i = 1:m
    for j = 1:n
    vwieght{c}=CT{i,j};
    vwieghtX{c}=CT{i,j}(2); 
    c=c+1;   
    end
end

% Теперь создание матрицы смежности (adjacency matrix)

[m,n] = size(vwieght); % размер исходного вектора

[J] = find( [vwieghtX{1:n}]==0); % находим нулевой элемент

for i = m:n % обнуляем массив
   for j = m:n
   AM(i,j)=0;
   end
end

[J] = find( [vwieghtX{1:n}]~=Inf); % находим индексы x выбранных яччеек (узлов)
for i = m:n
  if  vwieght{i}(2)~= Inf 
   for j = i+1:n
 if  vwieght{j}(2)~= Inf &&  vwieght{j}(2)==0
 AM(i,j)=1; 
 AM(j,i)=1;    
   break; 
 elseif vwieght{j}(2)~= Inf &&  vwieght{j}(2)~=0 && i> J(1) % нуль элемент 
     %обычно служит для дополнения и вернее всего(на данный момент уверенности 
     % конечно нет) следует после первого в цепи выбранного элемента    
 AM(i,j)=1; 
 AM(j,i)=1;
    break;
 end % if
 end
     end % if
end

V=[1.5 2-0.2; 2.5 2-0.2; 3.5 2-0.2; 4.5 2-0.2; % из -за включения режима hold переворачиваем граф
   1.5 3-0.2; 2.5 3-0.2; 3.5 3-0.2; 4.5 3-0.2;
   1.5 4-0.2; 2.5 4-0.2; 3.5 4-0.2; 4.5 4-0.2;
];% координаты вершин
n= 1:12; % количество вершин

gplot(AM(n,n),V,'-bo');
axis off
for j = 1:12, text(V(j,1),V(j,2), int2str(j),'FontSize',15); end

pause(1.9);
%text(1,4.3,['Result=', int2str(20)], 'HorizontalAlignment','center', 'BackgroundColor',[.7 .9 .7],'FontSize',15);
%text(4,4.3,['Simplex result=', num2str(SimplexResult(Mtemp))], 'HorizontalAlignment','center', 'BackgroundColor',[.7 .9 .7],'FontSize',15);

% форматируем матрицу для вывода в информационное окно
Cfor=Mtemp
[m,n]=size (Cfor)
 stc=[]
 stv=[]
for i = 1 : m  % цикл по пунктам отправления 
     for j =1 : n  % цикл по пунктам назначения
       if i==1 && j>1  
         Cfor {i, j}=num2str(Cfor {i, j});
       stc= strcat(stc,'[',Cfor{i,j},'     ]')

         else

          Cfor {i, j}=num2str(Cfor {i, j});
         stc= strcat(stc,'[',Cfor{i,j},']')    
         end    

       end
       stv=strvcat(stv, stc)
       stc=[]
  end

set(hndlList(1),'ForegroundColor','b')
set(hndlList(1),'FontWeight',' bold')

stv=strvcat(stv, strcat('Simplex result =',num2str(SimplexResult(Mtemp))))
set(hndlList(1),'String',stv);

%strcat('Simplex result =',num2str(fval))
%set(hndlList(1),'ForegroundColor','k')
%set(hndlList(1),'FontWeight','normal')

%---------------------
% оптимизируем по классической схеме

hndlList=get(gcf,'UserData'); % берем данные из массива 
OptimizationType= get(hndlList(5),'Value');% тип оптимизации 

if OptimizationType==1
combClassic(Mtemp)
elseif OptimizationType==2
combOptimization(Mtemp)
end    

end % function

function tempSum= SumCol(k,i,Mtemp)

tempSum=0;

for b = 2: i % суммируем данные в столбце 

     if Mtemp{b,k}(2) ~= Inf  

    tempSum=tempSum+Mtemp{b,k}(2);

     end 
     end 
     %------------------------------ 

end %function

function tempSum= SumRow(k,n,Mtemp)

tempSum=0;

for b = 2: n % суммируем данные в строке 

     if Mtemp{k, b}(2) ~= Inf  

    tempSum=tempSum+Mtemp{k, b}(2);

     end 
     end 
     %------------------------------ 

end %function

Create file - ResultGraph and copy there the following code.

function ResultGraph(arg)
%arg = {['ba'] [5] [10] [20] [15];
%     [10] [8 Inf] [3 Inf] [5 Inf] [2 10];
%     [15] [4 5] [1 10] [6 0] [7 Inf];
%     [25] [1 Inf] [9 Inf] [4 20] [3 5]};
figNumber=figure( ...
        'Visible','on', ...
        'NumberTitle','of', ...
        'Name','Linear - programming, result graph',...
        'Color',[1.0 1.0 1.0],...
        'Toolbar', 'none');
    axes( ...
        'Units','normalized', ...
        'Position',[0.1 0.3 0.8 0.59])

[m,n] = size(arg); % размер исходного массива
% Вначале преобразуем матрицу, удалив из нее информационные строки

CT= cell(m-1,n-1);

for i = 1:m-1
    for j = 1:n-1
    CT{i,j}=arg{i+1,j+1};    
    end
end

[m,n] = size(CT); % размер исходного массива

AM=zeros(m+n,m+n);

for i = 1:m
   AM(:)=0
    for j = 1:n
 if  CT{i,j}(2)~= Inf %&&  vwieght{j}(2)==0
 AM(i,j+n-1)=1; 
 AM(j+n-1,i)=1;    
 TT{i}=AM 

   end % if
   end

end

% из -за включения режима hold переворачиваем граф

V=[2 2-0.5; 2 3-0.5; 2 4-0.5; 
   3 1.5; 3 2.5; 3 3.5; 3 4.5; ];% координаты вершин
n= 1:m+n; % количество вершин
hold on
%subplot(2,1,1);
%for i = 1:m
gplot(TT{1}(n,n),V,'-og');
gplot(TT{2}(n,n),V,'-ob');
gplot(TT{3}(n,n),V,'-om');

%end
%gplot(AM(n,n),V,'-og');
axis off
label=['a1'; 'a2'; 'a3'; 'b1'; 'b2'; 'b3'; 'b4']

for j = n, text(V(j,1),V(j,2), label(j,1:end),'FontSize',15); end

Push button "Start", watch of iterations, after the optimization termination push the button "Show reesult graph". From the falling out menu select algorithm of optimization - "combinatoric classic | combinatoric opntimization"

Answer 2

[Alex]

I think I figured it out. I had to clear the offspring after every iteration. If I didn't clear them then the ox crossover operator would create offspring that were invalid routes. It seems to work now.