Genetic algorithm (n-queens problem) in MATLAB is stuck

Question

I need to solve the n-queens problem in MATLAB and I have tried with the code attached below. The problem that I have is that, at some point of the simulation the "checks" of the populations get stuck in a certain value (close to zero but not zero as I need). This happens 90% of the time, the other 10% gives a solution very fast (in a very few iterations). Maybe some parts are not well optimized but that is my best try. Feel free to ask anything about it, sorry if the comments are in spanish, thank you.

This is the code I used, if you run it multiple times, just a few of them will give you a solution.

% No siempre da soluciones, se suele quedar estancado con muy pocos jaques
% por población.

clc, clear
% Parámetros y vectores iniciales
N = 10000; % Número de iteraciones
n = 8; % Número de reinas
l = 8; % Longitud de cada lado del tablero
% Condición para que no haya más reinas que lados del tablero
while n > l
    n = n-1;
end
pob = 5; % Número de población
f = @(x) n*l - x; % Función de aptitud, entre mayor cantidad de jaques, menor será la puntuación.
tableros = cell(1,pob+2); % Celda en donde se guardarán los vectores de cada población
posiciones = cell(1,pob+2);
hijo1 = zeros(n,2);
hijo2 = zeros(n,2);
hijo1(1:n,1) = 0:n-1;
hijo2 = hijo1;

[posiciones_iniciales, tableros_iniciales, cont_jaques, puntuacion] = pob_inicial(pob, n, l, f);
posiciones = posiciones_iniciales;

% Padres
[valoresMayores, padres] = mink(cont_jaques,2);


cont = 1;
resultado = 0;
while cont ~= N
        pc = randi([1,n]); % Punto de cruzamiento
        hijo1 = posiciones{padres(1)};
        hijo2 = posiciones{padres(2)};
        hijo1(pc+1:n,2) = posiciones{padres(2)}(pc+1:n,2);
        hijo2(pc+1:n,2) = posiciones{padres(1)}(pc+1:n,2);
    if rand > 0.2 % Probabilidad de mutación del 80% del hijo 1
        m =  randi([1,n]);
        hijo1(m, 2) = randi([0,n-1]);     
    end
    if rand > 0.2 % Probabilidad de mutación del 80% del hijo 1
        m =  randi([1,n]);
        hijo2(m, 2) = randi([0,n-1]);     
    end
    posiciones{pob + 1} = hijo1;
    posiciones{pob + 2} = hijo2;

    % Posiciones de los tableros (para el cálculo de los jaques)
    for ii = 1:pob+2
    tableros{ii} = zeros(l,l); % Tablero (si en la casilla hay un 1, significa que en esa posición hay una reina)
    for k = 1:n
        tableros{ii}(posiciones{ii}(k,1)+1, posiciones{ii}(k,2)+1) = 1;
    end

    % Checar los jaques de los tableros iniciales
    n_jaques = jaques(posiciones{ii}(:,2)', l, ii, tableros);
    cont_jaques(ii) = n_jaques;
    
    % Prueba de aptitud
    puntuacion(ii) = f(cont_jaques(ii)); 
    end

    % Eliminados (los de menor aptitud)
    [valoresMenores, eliminados] = maxk(cont_jaques, 2);
    posiciones(:, eliminados) = [];
    cont_jaques(eliminados) = [];
    puntuacion(eliminados) = [];

    % Padres 
    [valoresMayores, padres] = mink(cont_jaques, 2);

    % Checar si hay soluciones
    for jj = 1:pob
        if cont_jaques(jj) == 0
            resultado = 1;
            break
        end
    end

    if resultado == 1
        break
    else
        cont = cont + 1;
    end
    
end

if cont == N
    fprintf('No se ha encontrado solución en %d iteraciones. \n', N)
else
    fprintf('Se encontró la solución despues de %d iteraciones. \n', cont)
end

for ii = 1:pob
    if cont_jaques(ii) == 0
        subplot(2, 3, 2)
        tablero(l)
        title('Tablero inicial', ii)
        reinas(posiciones_iniciales{ii}(:,1),posiciones_iniciales{ii}(:,2), n, 'red')
        subplot(2, 3, 5)
        tablero(l)
        title('Tablero final', ii)
        reinas(posiciones{ii}(:,1),posiciones{ii}(:,2), n, 'green')
    end
end

%% Función que genera las poblaciones iniciales aleatoriamente
function [posiciones_iniciales, tableros_iniciales, cont_jaques, puntuacion] = pob_inicial(pob, n, l, f)
% pob: Número de población
% n: Número de reinas
% l: Lado del tablero
% f: Función de aptitud

tableros_iniciales = cell(1,pob); % Celda en donde se guardarán los vectores de cada población
posiciones_iniciales = cell(1,pob); % Celda en donde se guardarán los vectores de cada población
cont_jaques = zeros(1,pob); % Vector de contador de jaques de cada tablero
puntuacion = zeros(1,pob); % Vector que guarda las puntuaciones de aptitud
for ii = 1:pob
    % Poblaciones iniciales aleatorias
    r = randi([0,l-1],1,n);
    for jj = 0:n-1
        posiciones_iniciales{ii}(jj+1,1) = jj;
        posiciones_iniciales{ii}(jj+1,2) = r(jj+1); 
    end
    
    % Posiciones del tablero
    tableros_iniciales{ii} = zeros(l,l); % Tablero (si en la casilla hay un 1, significa que en esa posición hay una reina)
    for k = 1:n
        tableros_iniciales{ii}(posiciones_iniciales{ii}(k,1)+1, posiciones_iniciales{ii}(k,2)+1) = 1;
    end

    % Checar los jaques de los tableros iniciales
    n_jaques = jaques(posiciones_iniciales{ii}(:,2), l, ii, tableros_iniciales);
    cont_jaques(ii) = n_jaques;
    
    % Prueba de aptitud
    puntuacion(ii) = f(cont_jaques(ii)); 
end
end


%% Función de comprobación de jaques
function n_jaques = jaques(ry, l, n_tablero, celda_tableros)
% Función que comprueba cuántas reinas hay en jaque dado un tablero
n_jaques = 0;
% Verticales
By = unique(ry);
n_repetidos_y = histc(ry, By); 
c_unicos_y = size(n_repetidos_y);
for jj = 1:c_unicos_y(2)
    if n_repetidos_y(jj) > 1
        combinatoria_y = factorial(n_repetidos_y(jj))/(factorial(n_repetidos_y(jj) - 2));
        n_jaques = n_jaques + combinatoria_y; % Se suman la cantidad de jaques posibles de la columna
    end
end

% Diagonales con pendiente negativa
for k = 0:l
% Diagonales con k positiva
diagonal1 = diag(celda_tableros{n_tablero},k);
n_reinas_d1 = numel(find(diagonal1==1));
if n_reinas_d1 == 2
   n_jaques = n_jaques + 1;
elseif n_reinas_d1 > 2
    combinatoria_d1 = factorial(n_reinas_d1)/(factorial(2)*(factorial(n_reinas_d1 - 2)));
    n_jaques = n_jaques + combinatoria_d1;
end
% Diagonales con k negativa
end
for k = -1:-1:-l
diagonal2 = diag(celda_tableros{n_tablero},k);
n_reinas_d2 = numel(find(diagonal2==1));
if n_reinas_d2 == 2
   n_jaques = n_jaques + 1;
elseif n_reinas_d2 > 2
    combinatoria_d2 = factorial(n_reinas_d2)/(factorial(2)*(factorial(n_reinas_d2 - 2)));
    n_jaques = n_jaques + combinatoria_d2; 
end
end

% Diagonales con pendiente positiva
for k = 0:l
% Diagonales con k positiva
diagonal1 = diag(fliplr(celda_tableros{n_tablero}),k);
n_reinas_d1 = numel(find(diagonal1==1));
if n_reinas_d1 == 2
   n_jaques = n_jaques + 1;
elseif n_reinas_d1 > 2
    combinatoria_d1 = factorial(n_reinas_d1)/(factorial(2)*(factorial(n_reinas_d1 - 2)));
    n_jaques = n_jaques + combinatoria_d1;
end
% Diagonales con k negativa
end
for k = -1:-1:-l
diagonal2 = diag(fliplr(celda_tableros{n_tablero}),k);
n_reinas_d2 = numel(find(diagonal2==1));
if n_reinas_d2 == 2
   n_jaques = n_jaques + 1;
elseif n_reinas_d2 > 2
    combinatoria_d2 = factorial(n_reinas_d2)/(factorial(2)*(factorial(n_reinas_d2 - 2)));
    n_jaques = n_jaques + combinatoria_d2; 
end

end
end

%% Función de gráfica del tablero
function tablero(l)
%   Función que crea un tablero de lado l
x1=0; x2=l;
y1=0; y2=l;
x = [x1, x2, x2, x1, x1];
y = [y1, y1, y2, y2, y1];
plot(x, y, 'k-', 'LineWidth', 1.5);
hold on
xlim([-1, l+1]);
ylim([-1, l+1]);
% Grid
for i = 1:l
plot([0,l],[i,i], 'k-', 'LineWidth', 1.5)
plot([i,i],[0,l],'k-', 'LineWidth', 1.5)
end
end

%% Función de gráfica de las reinas
function reinas(rx,ry,n,color)
% Función que grafica n cantidad de reinas en las posiciones en x y y dadas por los vectores rx y ry.
for n = 1:n
    g_reinas = plot(rx(n)+0.5,ry(n)+0.5, 'or');
    g_reinas.MarkerFaceColor = [color];
    g_reinas.MarkerSize = 6;
    g_reinas.MarkerEdgeColor = [color];
    % legend(reinas,'Reinas')
end
hold off
end

Answer 1

Perhaps using a random function to find all populations does not allow you to find all solutions quickly.

You could search through all possible permutations n! of the column indexes, example of the different arrangements of the column indexes are (1. 5. 8. 6. 3. 7. 2. 4) And (1. 6. 8. 3. 7. 4. 2. 5) which in the first case correspond to positions of the queens (1,1);(2,5);(3,8);(4,6);(5,3);(6,7);(7,2) ;(8,4).

Or use an optimized solution like this one that doesn't evaluate all permutations but jumps are made if certain conditions are met.

This code has been tested with SCILAB but should also work in MATLAB:

n=8;
count=0;
Index_col=ones(1,n);
i = 2;
while (1 == 1)
    ck_var = 1;
    while(ck_var == 1);
        ck_var = 0;
        for j=1:i-1
            if ((Index_col(1,i) == Index_col(1,j)) || (abs(Index_col(1,i) - Index_col(1,j)) == i - j))
                if (Index_col(1,i) < n)
                    Index_col(1,i)=Index_col(1,i)+1;
                    ck_var = 1;
                else
                    ck_var = 2;
                end
                break;
            end
        end
    end
    if (i == n || ck_var == 2)
        if (ck_var == 0)
            //           print board
            for jc=1:n
                Board(jc,:)=zeros(1,n);
                Board(jc,Index_col(1,jc))=1;
            end
            Board
            //
            count=count+1;
        end
        while (Index_col(1,i) == n && i > 1)
            i=i-1;
        end
        if (Index_col(1,1) == n && i == 1)
            break;  
        else                
            Index_col(1,i)=Index_col(1,i)+1;
            Index_col(1,i+1:n) = 1;
            i = 2;
        end
    else
        i=i+1;
    end
end
count