function saEightQueensExample
close all; clc
%  Copyright (c) 2010, The MathWorks, Inc.
%  All rights reserved.
plot_board;

q = 1:8;
plot_board;
plot_queens(q, true);

global prev_q
prev_q = q;
title({'Place 8 queens in non attacking positions. Press any key...'})
pause

disp('Calling Simulated Annealing:')
saopts = saoptimset('DataType', 'custom', ...
    'AnnealingFcn', @my_anneal_fcn, ...
    'OutputFcns', {@my_output_fcn}, ...
    'ObjectiveLimit', 0, ...
    'InitialTemperature', 10);
[q,fval,exitFlag,output] = ...
    simulannealbnd(@cost_function, q, [], [], saopts); %#ok<*NASGU>

end


% Cost function: Compute number of queens that can attack each other
function f = cost_function(q)

% Create a matrix with 1's where queens are present
A = zeros(8);
for i = 1:numel(q); A(i,q(i)) = 1; end

f = 0;    % Initialize cost function

% The way the 'x' vector is set up, it is guaranteed that any
% row or column will have exactly one queen.
% We also want to ensure there is only one queen per diagonal

% Go along regular diags and count the number of queens that share a diag
for i = -6:6

    % Count number of queens on this diag
    x = sum(diag(A,i));

    % If more than one queen on this diag, increment cost function
    if x > 1
        f = f + x;
    end

end

% Go along rotated diags i.e. anti-diagonals, and do the same thing
B = rot90(A);
for i = -6:6

    % Count number of queens on this anti-diagonal
    x = sum(diag(B,i));

    % If more than one queen on this anti-diagonal, increment cost function
    if x > 1
        f = f + x;
    end

end

% Thats it! We have computed the cost function

end

% The annealing function swaps out two random queens
% Higher the temperature, more pairs are swapped
function q = my_anneal_fcn(optimValues, problem) %#ok<*INUSD>

q = optimValues.x;
N = numel(q);
T = fix(mean(optimValues.temperature)) + 1;

% Swap two random queen positions on board 'T' times
for i = 1:T

    % Make sure we have two distinct queens to swap
    while 1
        idx = randi(N, 2, 1);
        if idx(1) ~= idx(2); break; end
    end

    % Swap the two queens picked out randomly
    temp = q(idx(1));
    q(idx(1)) = q(idx(2));
    q(idx(2)) = temp;

end

end


function [stop, options, optchanged] = ...
    my_output_fcn(options, optimValues, flag)
global prev_q
stop = false; optchanged = false;
fprintf('Iteration = %4d, Cost = %d, Temperature = %d\n', ...
    optimValues.iteration, optimValues.fval, ...
    fix(mean(optimValues.temperature))+1);
end_q = optimValues.x;

% Move queens from positions 'prev_q' to 'end_q'
q = prev_q;
if all(isequal(q, end_q))
    title(sprintf('Number of attacking queens = %d', optimValues.fval))
    plot_queens(q);
else
    while 1
        for i = 1:numel(q)
            if q(i) < end_q(i); q(i) = q(i) + 1; end
            if q(i) > end_q(i); q(i) = q(i) - 1; end
        end
        title(sprintf('Number of attacking queens = %d', optimValues.fval))
        plot_queens(q);
        pause(.05)
        if all(isequal(q, end_q)); break; end
    end
end
prev_q = end_q;
end

% Function to plot a blank board with no pieces on it
function plot_board
persistent b

if isempty(b)
    b = zeros(800, 800);
    for i = 1:8
        for j = 1:8
            b((i-1)*100+1:i*100, (j-1)*100+1:j*100) = mod(i+j+1,2)+1;
        end
    end
end
image(b);
colormap([25 48 92; 255 255 255; ] ./ 255)
axis image; axis([1 800 1 800]);
set(gca, 'XTick', [50:100:750])
set(gca, 'YTick', [50:100:750])
set(gca, 'XTickLabel', char('a'+(0:7)'))
set(gca, 'YTickLabel', char('0'+(8:-1:1)'))
hold on

end

% Function to position queens on the board according to configuration
function plot_queens(pos, init)
persistent queens

if nargin > 1
    row = 1; col = pos(row);
    X = [(col-1)*100+1 col*100];
    Y = [(row-1)*100+1 row*100];
    queens(1) = image(X, Y, imread('queen.png'));
    c = get(queens(1), 'CData'); sc = sum(c,3);
    A = ones(size(sc)); A( sc == 255*3 ) = 0;
    set(queens(1), 'AlphaData', A);
    for k = 2:8
        queens(k) = copyobj(queens(1), gca);
    end
end

% Position the "queen" images according to "pos"
for row = 1:numel(pos)
    col = pos(row);
    X = [(col-1)*100+1 col*100];
    Y = [(row-1)*100+1 row*100];
    set(queens(row), 'XData', X, 'YData', Y);
end

end
Calling Simulated Annealing:
Iteration =    0, Cost = 8, Temperature = 11
Iteration =    1, Cost = 6, Temperature = 10
Iteration =    2, Cost = 4, Temperature = 9
Iteration =    3, Cost = 4, Temperature = 9
Iteration =    4, Cost = 4, Temperature = 8
Iteration =    5, Cost = 4, Temperature = 8
Iteration =    6, Cost = 4, Temperature = 7
Iteration =    7, Cost = 8, Temperature = 7
Iteration =    8, Cost = 8, Temperature = 7
Iteration =    9, Cost = 8, Temperature = 6
Iteration =   10, Cost = 2, Temperature = 6
Iteration =   11, Cost = 2, Temperature = 6
Iteration =   12, Cost = 2, Temperature = 6
Iteration =   13, Cost = 2, Temperature = 5
Iteration =   14, Cost = 2, Temperature = 5
Iteration =   15, Cost = 8, Temperature = 5
Iteration =   16, Cost = 6, Temperature = 5
Iteration =   17, Cost = 6, Temperature = 4
Iteration =   18, Cost = 6, Temperature = 4
Iteration =   19, Cost = 6, Temperature = 4
Iteration =   20, Cost = 6, Temperature = 4
Iteration =   21, Cost = 6, Temperature = 4
Iteration =   22, Cost = 6, Temperature = 4
Iteration =   23, Cost = 6, Temperature = 3
Iteration =   24, Cost = 6, Temperature = 3
Iteration =   25, Cost = 8, Temperature = 3
Iteration =   26, Cost = 6, Temperature = 3
Iteration =   27, Cost = 6, Temperature = 3
Iteration =   28, Cost = 6, Temperature = 3
Iteration =   29, Cost = 6, Temperature = 3
Iteration =   30, Cost = 6, Temperature = 3
Iteration =   31, Cost = 6, Temperature = 2
Iteration =   32, Cost = 6, Temperature = 2
Iteration =   33, Cost = 6, Temperature = 2
Iteration =   34, Cost = 6, Temperature = 2
Iteration =   35, Cost = 6, Temperature = 2
Iteration =   36, Cost = 6, Temperature = 2
Iteration =   37, Cost = 5, Temperature = 2
Iteration =   38, Cost = 5, Temperature = 2
Iteration =   39, Cost = 4, Temperature = 2
Iteration =   40, Cost = 4, Temperature = 2
Iteration =   41, Cost = 4, Temperature = 2
Iteration =   42, Cost = 4, Temperature = 2
Iteration =   43, Cost = 5, Temperature = 2
Iteration =   44, Cost = 5, Temperature = 1
Iteration =   45, Cost = 5, Temperature = 1
Iteration =   46, Cost = 8, Temperature = 1
Iteration =   47, Cost = 8, Temperature = 1
Iteration =   48, Cost = 7, Temperature = 1
Iteration =   49, Cost = 7, Temperature = 1
Iteration =   50, Cost = 8, Temperature = 1
Iteration =   51, Cost = 6, Temperature = 1
Iteration =   52, Cost = 6, Temperature = 1
Iteration =   53, Cost = 5, Temperature = 1
Iteration =   54, Cost = 5, Temperature = 1
Iteration =   55, Cost = 5, Temperature = 1
Iteration =   56, Cost = 5, Temperature = 1
Iteration =   57, Cost = 2, Temperature = 1
Iteration =   58, Cost = 2, Temperature = 1
Iteration =   59, Cost = 2, Temperature = 1
Iteration =   60, Cost = 2, Temperature = 1
Iteration =   61, Cost = 2, Temperature = 1
Iteration =   62, Cost = 2, Temperature = 1
Iteration =   63, Cost = 2, Temperature = 1
Iteration =   64, Cost = 2, Temperature = 1
Iteration =   65, Cost = 2, Temperature = 1
Iteration =   66, Cost = 2, Temperature = 1
Iteration =   67, Cost = 2, Temperature = 1
Iteration =   68, Cost = 2, Temperature = 1
Iteration =   69, Cost = 2, Temperature = 1
Iteration =   70, Cost = 2, Temperature = 1
Iteration =   71, Cost = 2, Temperature = 1
Iteration =   72, Cost = 0, Temperature = 1
Iteration =   72, Cost = 0, Temperature = 1
Optimization terminated: best function value reached options.ObjectiveLimit.