function saGalacticTravelingSalesmanExample
% Ref: http://www.essex1.com/people/speer/starmodel.html

%  Copyright (c) 2010, The MathWorks, Inc.
%  All rights reserved.
clc; close all;
global xc yc zc N
global bestfval_log

% load star positions from 'stars.txt'
stars = importdata('stars.txt');
N = size(stars.data, 1);
xc = stars.data(:,1); yc = stars.data(:,2); zc = stars.data(:,3);
names = stars.rowheaders;
plot_stars(xc, yc, zc, names);
bestfval_log = [];

% set up schedule, and plot it
schedule = 1:N;
disp('Initial Schedule:')
disp(schedule);
plot_schedule(schedule);
title(sprintf('Initial Schedule: %s \n Distance = %f light year', ...
    num2str(schedule), cost_function(schedule)))
rotate3d
disp('Initial Total Distance:')
disp(cost_function(schedule));
pause

% set up options and call simulated annealing
rotate3d;
saopts = saoptimset('DataType', 'custom', ...
    'AnnealingFcn', @sa_anneal_fcn, ...
    'OutputFcns', @sa_output_fcn, ...
    'StallIterLimit', 100, ...
    'ReannealInterval', Inf, ...
    'InitialTemperature', 20);
[schedule,fval,exitFlag,output] = ...
    simulannealbnd(@cost_function, schedule, [], [], saopts); %#ok<*NASGU>

figure; plot(bestfval_log, '.-');
title 'Convergence for Galactic Traveling Salesman Problem';
xlabel 'Iteration'; ylabel 'Total Distance Traveled (light year)';

end

function d = cost_function(schedule)
global xc yc zc N
d = 0;
for i = 1:N-1
    d = d + ( ( xc(schedule(i)) - xc(schedule(i+1)) )^2 ...
        + ( yc(schedule(i)) - yc(schedule(i+1)) )^2 ...
        + ( zc(schedule(i)) - zc(schedule(i+1)) )^2) ^.5 ;
end
i = N;
d = d + ( ( xc(schedule(i)) - xc(schedule(1)) )^2 ...
    + ( yc(schedule(i)) - yc(schedule(1)) )^2 ...
    + ( zc(schedule(i)) - zc(schedule(1)) )^2) ^.5 ;
end

function schedule = sa_anneal_fcn(optimValues,problem) %#ok<*INUSD>
global N
schedule = optimValues.x;
T = fix(mean(optimValues.temperature)) + 1;
for i = 1:T
    idx = randi(N, 2, 1);
    while idx(1) == idx(2)
        idx = randi(N, 2, 1);
    end
    temp = schedule(idx(1));
    schedule(idx(1)) = schedule(idx(2));
    schedule(idx(2)) = temp;
end
end


function plot_stars(xc, yc, zc, names)
N = numel(xc);
f = figure('NumberTitle', 'off', ...
    'Name', 'Galactic Traveling Salesman Problem');
ax = axes('Color', [0 0 0]);
hold on;
textargs = {'VerticalAlignment', 'top', ...
    'HorizontalAlignment', 'center', ...
    'FontName', 'Gill Sans MT'};
plot3(xc(1), yc(1), zc(1), 'y*');
text(xc(1), yc(1), zc(1)-.5, names{1}, textargs{:}, 'Color', 'y');
for i = 2:N
    plot3(xc(i), yc(i), zc(i), 'w*');
    text(xc(i), yc(i), zc(i)-.5, names{i}, textargs{:}, 'Color', 'w');
end
% line([0 0], [0 0], [-15 +15], 'Color', 'g')
% line([-15 +15], [0 0], [0 0], 'Color', 'g')
% line([0 0], [-15 +15], [0 0], 'Color', 'g')
view(3); rotate3d;
end

function plot_schedule(schedule)
global xc yc zc
persistent l
if isempty(l) || ~ishghandle(l); l = plot3(0,0,0,'b-'); rotate3d; end
set(l, ...
    'XData', [xc(schedule); xc(schedule(1))], ...
    'YData', [yc(schedule); yc(schedule(1))], ...
    'ZData', [zc(schedule); zc(schedule(1))] );
title(sprintf('Best schedule found so far: %s \n Distance = %f light year', ...
    num2str(schedule), cost_function(schedule)))
end

function [stop,options,optchanged] = ...
    sa_output_fcn(options,optimValues,flag)
global bestfval_log
stop = false; optchanged = false;
if strcmpi(flag, 'init'); figure(1); end
fprintf('Iter=%d, Fval=%8.6f, Steps = %d\n', ...
    optimValues.iteration, optimValues.fval, ...
    fix(mean(optimValues.temperature))+1);
bestfval_log(end+1) = optimValues.bestfval;
plot_schedule(optimValues.x);
camorbit(1,0);
pause(.05);
end
Initial Schedule:
  Columns 1 through 14

     1     2     3     4     5     6     7     8     9    10    11    12    13    14

  Columns 15 through 16

    15    16

Initial Total Distance:
  209.2406

Iter=0, Fval=209.240605, Steps = 21
Iter=1, Fval=185.807661, Steps = 19
Iter=2, Fval=185.807661, Steps = 18
Iter=3, Fval=198.240108, Steps = 17
Iter=4, Fval=198.240108, Steps = 16
Iter=5, Fval=195.845394, Steps = 15
Iter=6, Fval=225.972711, Steps = 14
Iter=7, Fval=210.939784, Steps = 14
Iter=8, Fval=210.939784, Steps = 13
Iter=9, Fval=196.012875, Steps = 12
Iter=10, Fval=197.121722, Steps = 12
Iter=11, Fval=197.121722, Steps = 11
Iter=12, Fval=200.420900, Steps = 11
Iter=13, Fval=200.420900, Steps = 10
Iter=14, Fval=200.420900, Steps = 10
Iter=15, Fval=198.545301, Steps = 9
Iter=16, Fval=198.545301, Steps = 9
Iter=17, Fval=198.545301, Steps = 8
Iter=18, Fval=195.840213, Steps = 8
Iter=19, Fval=195.840213, Steps = 8
Iter=20, Fval=187.229283, Steps = 7
Iter=21, Fval=197.208381, Steps = 7
Iter=22, Fval=189.479454, Steps = 7
Iter=23, Fval=189.479454, Steps = 6
Iter=24, Fval=189.479454, Steps = 6
Iter=25, Fval=189.479454, Steps = 6
Iter=26, Fval=189.479454, Steps = 6
Iter=27, Fval=189.479454, Steps = 5
Iter=28, Fval=188.259237, Steps = 5
Iter=29, Fval=188.259237, Steps = 5
Iter=30, Fval=188.259237, Steps = 5
Iter=31, Fval=188.259237, Steps = 4
Iter=32, Fval=175.534203, Steps = 4
Iter=33, Fval=175.534203, Steps = 4
Iter=34, Fval=162.941756, Steps = 4
Iter=35, Fval=162.941756, Steps = 4
Iter=36, Fval=162.941756, Steps = 3
Iter=37, Fval=162.941756, Steps = 3
Iter=38, Fval=162.941756, Steps = 3
Iter=39, Fval=162.941756, Steps = 3
Iter=40, Fval=162.941756, Steps = 3
Iter=41, Fval=157.565950, Steps = 3
Iter=42, Fval=157.565950, Steps = 3
Iter=43, Fval=157.565950, Steps = 3
Iter=44, Fval=157.565950, Steps = 2
Iter=45, Fval=157.565950, Steps = 2
Iter=46, Fval=157.565950, Steps = 2
Iter=47, Fval=157.565950, Steps = 2
Iter=48, Fval=157.565950, Steps = 2
Iter=49, Fval=157.565950, Steps = 2
Iter=50, Fval=157.565950, Steps = 2
Iter=51, Fval=157.565950, Steps = 2
Iter=52, Fval=157.565950, Steps = 2
Iter=53, Fval=157.565950, Steps = 2
Iter=54, Fval=157.565950, Steps = 2
Iter=55, Fval=157.565950, Steps = 2
Iter=56, Fval=157.565950, Steps = 2
Iter=57, Fval=152.431211, Steps = 2
Iter=58, Fval=152.431211, Steps = 1
Iter=59, Fval=151.483126, Steps = 1
Iter=60, Fval=151.483126, Steps = 1
Iter=61, Fval=151.483126, Steps = 1
Iter=62, Fval=151.483126, Steps = 1
Iter=63, Fval=151.483126, Steps = 1
Iter=64, Fval=144.872901, Steps = 1
Iter=65, Fval=144.872901, Steps = 1
Iter=66, Fval=144.872901, Steps = 1
Iter=67, Fval=144.872901, Steps = 1
Iter=68, Fval=144.872901, Steps = 1
Iter=69, Fval=144.849413, Steps = 1
Iter=70, Fval=144.849413, Steps = 1
Iter=71, Fval=144.849413, Steps = 1
Iter=72, Fval=144.849413, Steps = 1
Iter=73, Fval=144.589925, Steps = 1
Iter=74, Fval=144.589925, Steps = 1
Iter=75, Fval=144.589925, Steps = 1
Iter=76, Fval=144.589925, Steps = 1
Iter=77, Fval=144.589925, Steps = 1
Iter=78, Fval=144.589925, Steps = 1
Iter=79, Fval=144.589925, Steps = 1
Iter=80, Fval=144.589925, Steps = 1
Iter=81, Fval=144.589925, Steps = 1
Iter=82, Fval=144.589925, Steps = 1
Iter=83, Fval=144.589925, Steps = 1
Iter=84, Fval=144.589925, Steps = 1
Iter=85, Fval=144.589925, Steps = 1
Iter=86, Fval=144.589925, Steps = 1
Iter=87, Fval=144.589925, Steps = 1
Iter=88, Fval=144.589925, Steps = 1
Iter=89, Fval=142.892770, Steps = 1
Iter=90, Fval=142.892770, Steps = 1
Iter=91, Fval=141.368615, Steps = 1
Iter=92, Fval=141.368615, Steps = 1
Iter=93, Fval=141.368615, Steps = 1
Iter=94, Fval=141.368615, Steps = 1
Iter=95, Fval=141.368615, Steps = 1
Iter=96, Fval=141.368615, Steps = 1
Iter=97, Fval=141.368615, Steps = 1
Iter=98, Fval=141.368615, Steps = 1
Iter=99, Fval=141.368615, Steps = 1
Iter=100, Fval=141.368615, Steps = 1
Iter=101, Fval=141.368615, Steps = 1
Iter=102, Fval=141.368615, Steps = 1
Iter=103, Fval=141.368615, Steps = 1
Iter=104, Fval=141.368615, Steps = 1
Iter=105, Fval=141.368615, Steps = 1
Iter=106, Fval=141.368615, Steps = 1
Iter=107, Fval=141.368615, Steps = 1
Iter=108, Fval=141.368615, Steps = 1
Iter=109, Fval=140.768329, Steps = 1
Iter=110, Fval=140.768329, Steps = 1
Iter=111, Fval=140.768329, Steps = 1
Iter=112, Fval=140.369359, Steps = 1
Iter=113, Fval=140.369359, Steps = 1
Iter=114, Fval=140.369359, Steps = 1
Iter=115, Fval=140.369359, Steps = 1
Iter=116, Fval=140.369359, Steps = 1
Iter=117, Fval=140.369359, Steps = 1
Iter=118, Fval=140.369359, Steps = 1
Iter=119, Fval=140.369359, Steps = 1
Iter=120, Fval=140.369359, Steps = 1
Iter=121, Fval=140.369359, Steps = 1
Iter=122, Fval=140.369359, Steps = 1
Iter=123, Fval=140.369359, Steps = 1
Iter=124, Fval=140.369359, Steps = 1
Iter=125, Fval=140.369359, Steps = 1
Iter=126, Fval=140.369359, Steps = 1
Iter=127, Fval=140.369359, Steps = 1
Iter=128, Fval=139.556289, Steps = 1
Iter=129, Fval=139.556289, Steps = 1
Iter=130, Fval=139.556289, Steps = 1
Iter=131, Fval=139.556289, Steps = 1
Iter=132, Fval=139.556289, Steps = 1
Iter=133, Fval=139.556289, Steps = 1
Iter=134, Fval=139.556289, Steps = 1
Iter=135, Fval=139.556289, Steps = 1
Iter=136, Fval=139.556289, Steps = 1
Iter=137, Fval=139.556289, Steps = 1
Iter=138, Fval=139.556289, Steps = 1
Iter=139, Fval=139.556289, Steps = 1
Iter=140, Fval=139.556289, Steps = 1
Iter=141, Fval=139.556289, Steps = 1
Iter=142, Fval=139.556289, Steps = 1
Iter=143, Fval=139.556289, Steps = 1
Iter=144, Fval=139.556289, Steps = 1
Iter=145, Fval=139.556289, Steps = 1
Iter=146, Fval=139.556289, Steps = 1
Iter=147, Fval=139.556289, Steps = 1
Iter=148, Fval=139.556289, Steps = 1
Iter=149, Fval=139.556289, Steps = 1
Iter=150, Fval=138.579097, Steps = 1
Iter=151, Fval=138.579097, Steps = 1
Iter=152, Fval=138.579097, Steps = 1
Iter=153, Fval=138.579097, Steps = 1
Iter=154, Fval=138.579097, Steps = 1
Iter=155, Fval=138.579097, Steps = 1
Iter=156, Fval=138.579097, Steps = 1
Iter=157, Fval=138.579097, Steps = 1
Iter=158, Fval=138.579097, Steps = 1
Iter=159, Fval=138.579097, Steps = 1
Iter=160, Fval=138.579097, Steps = 1
Iter=161, Fval=138.579097, Steps = 1
Iter=162, Fval=138.579097, Steps = 1
Iter=163, Fval=138.579097, Steps = 1
Iter=164, Fval=138.579097, Steps = 1
Iter=165, Fval=138.579097, Steps = 1
Iter=166, Fval=132.467485, Steps = 1
Iter=167, Fval=127.570022, Steps = 1
Iter=168, Fval=127.570022, Steps = 1
Iter=169, Fval=127.570022, Steps = 1
Iter=170, Fval=127.570022, Steps = 1
Iter=171, Fval=127.570022, Steps = 1
Iter=172, Fval=127.570022, Steps = 1
Iter=173, Fval=127.570022, Steps = 1
Iter=174, Fval=127.570022, Steps = 1
Iter=175, Fval=127.570022, Steps = 1
Iter=176, Fval=127.570022, Steps = 1
Iter=177, Fval=127.570022, Steps = 1
Iter=178, Fval=127.570022, Steps = 1
Iter=179, Fval=127.570022, Steps = 1
Iter=180, Fval=127.570022, Steps = 1
Iter=181, Fval=127.570022, Steps = 1
Iter=182, Fval=127.570022, Steps = 1
Iter=183, Fval=127.570022, Steps = 1
Iter=184, Fval=127.570022, Steps = 1
Iter=185, Fval=127.570022, Steps = 1
Iter=186, Fval=127.570022, Steps = 1
Iter=187, Fval=127.570022, Steps = 1
Iter=188, Fval=127.570022, Steps = 1
Iter=189, Fval=127.570022, Steps = 1
Iter=190, Fval=127.570022, Steps = 1
Iter=191, Fval=127.570022, Steps = 1
Iter=192, Fval=127.570022, Steps = 1
Iter=193, Fval=127.570022, Steps = 1
Iter=194, Fval=127.570022, Steps = 1
Iter=195, Fval=127.570022, Steps = 1
Iter=196, Fval=127.570022, Steps = 1
Iter=197, Fval=127.570022, Steps = 1
Iter=198, Fval=127.570022, Steps = 1
Iter=199, Fval=127.570022, Steps = 1
Iter=200, Fval=127.570022, Steps = 1
Iter=201, Fval=127.570022, Steps = 1
Iter=202, Fval=127.570022, Steps = 1
Iter=203, Fval=127.570022, Steps = 1
Iter=204, Fval=127.570022, Steps = 1
Iter=205, Fval=127.570022, Steps = 1
Iter=206, Fval=127.570022, Steps = 1
Iter=207, Fval=127.570022, Steps = 1
Iter=208, Fval=127.570022, Steps = 1
Iter=209, Fval=127.570022, Steps = 1
Iter=210, Fval=127.570022, Steps = 1
Iter=211, Fval=127.570022, Steps = 1
Iter=212, Fval=127.570022, Steps = 1
Iter=213, Fval=127.570022, Steps = 1
Iter=214, Fval=127.570022, Steps = 1
Iter=215, Fval=127.570022, Steps = 1
Iter=216, Fval=120.914359, Steps = 1
Iter=217, Fval=120.914359, Steps = 1
Iter=218, Fval=120.914359, Steps = 1
Iter=219, Fval=120.914359, Steps = 1
Iter=220, Fval=120.914359, Steps = 1
Iter=221, Fval=120.914359, Steps = 1
Iter=222, Fval=120.914359, Steps = 1
Iter=223, Fval=120.914359, Steps = 1
Iter=224, Fval=120.914359, Steps = 1
Iter=225, Fval=120.914359, Steps = 1
Iter=226, Fval=120.914359, Steps = 1
Iter=227, Fval=120.914359, Steps = 1
Iter=228, Fval=120.914359, Steps = 1
Iter=229, Fval=119.048937, Steps = 1
Iter=230, Fval=119.048937, Steps = 1
Iter=231, Fval=119.048937, Steps = 1
Iter=232, Fval=119.048937, Steps = 1
Iter=233, Fval=119.048937, Steps = 1
Iter=234, Fval=119.048937, Steps = 1
Iter=235, Fval=119.048937, Steps = 1
Iter=236, Fval=119.048937, Steps = 1
Iter=237, Fval=119.048937, Steps = 1
Iter=238, Fval=119.048937, Steps = 1
Iter=239, Fval=119.048937, Steps = 1
Iter=240, Fval=119.048937, Steps = 1
Iter=241, Fval=119.048937, Steps = 1
Iter=242, Fval=119.048937, Steps = 1
Iter=243, Fval=119.048937, Steps = 1
Iter=244, Fval=119.048937, Steps = 1
Iter=245, Fval=119.048937, Steps = 1
Iter=246, Fval=119.048937, Steps = 1
Iter=247, Fval=119.048937, Steps = 1
Iter=248, Fval=119.048937, Steps = 1
Iter=249, Fval=119.048937, Steps = 1
Iter=250, Fval=119.048937, Steps = 1
Iter=251, Fval=119.048937, Steps = 1
Iter=252, Fval=119.048937, Steps = 1
Iter=253, Fval=119.048937, Steps = 1
Iter=254, Fval=119.048937, Steps = 1
Iter=255, Fval=119.048937, Steps = 1
Iter=256, Fval=119.048937, Steps = 1
Iter=257, Fval=119.048937, Steps = 1
Iter=258, Fval=119.048937, Steps = 1
Iter=259, Fval=119.048937, Steps = 1
Iter=260, Fval=119.048937, Steps = 1
Iter=261, Fval=119.048937, Steps = 1
Iter=262, Fval=119.048937, Steps = 1
Iter=263, Fval=119.048937, Steps = 1
Iter=264, Fval=119.048937, Steps = 1
Iter=265, Fval=119.048937, Steps = 1
Iter=266, Fval=119.048937, Steps = 1
Iter=267, Fval=119.048937, Steps = 1
Iter=268, Fval=119.048937, Steps = 1
Iter=269, Fval=119.048937, Steps = 1
Iter=270, Fval=119.048937, Steps = 1
Iter=271, Fval=119.048937, Steps = 1
Iter=272, Fval=119.048937, Steps = 1
Iter=273, Fval=119.048937, Steps = 1
Iter=274, Fval=119.048937, Steps = 1
Iter=275, Fval=119.048937, Steps = 1
Iter=276, Fval=119.048937, Steps = 1
Iter=277, Fval=119.048937, Steps = 1
Iter=278, Fval=119.048937, Steps = 1
Iter=279, Fval=119.048937, Steps = 1
Iter=280, Fval=119.048937, Steps = 1
Iter=281, Fval=119.048937, Steps = 1
Iter=282, Fval=119.048937, Steps = 1
Iter=283, Fval=119.048937, Steps = 1
Iter=284, Fval=119.048937, Steps = 1
Iter=285, Fval=119.048937, Steps = 1
Iter=286, Fval=119.048937, Steps = 1
Iter=287, Fval=119.048937, Steps = 1
Iter=288, Fval=119.048937, Steps = 1
Iter=289, Fval=119.048937, Steps = 1
Iter=290, Fval=119.048937, Steps = 1
Iter=291, Fval=119.048937, Steps = 1
Iter=292, Fval=119.048937, Steps = 1
Iter=293, Fval=119.048937, Steps = 1
Iter=294, Fval=119.048937, Steps = 1
Iter=295, Fval=119.048937, Steps = 1
Iter=296, Fval=119.048937, Steps = 1
Iter=297, Fval=119.048937, Steps = 1
Iter=298, Fval=119.048937, Steps = 1
Iter=299, Fval=119.048937, Steps = 1
Iter=300, Fval=119.048937, Steps = 1
Iter=301, Fval=119.048937, Steps = 1
Iter=302, Fval=119.048937, Steps = 1
Iter=303, Fval=119.048937, Steps = 1
Iter=304, Fval=119.048937, Steps = 1
Iter=305, Fval=119.048937, Steps = 1
Iter=306, Fval=119.048937, Steps = 1
Iter=307, Fval=119.048937, Steps = 1
Iter=308, Fval=119.048937, Steps = 1
Iter=309, Fval=119.048937, Steps = 1
Iter=310, Fval=119.048937, Steps = 1
Iter=311, Fval=119.048937, Steps = 1
Iter=312, Fval=119.048937, Steps = 1
Iter=313, Fval=119.048937, Steps = 1
Iter=314, Fval=119.048937, Steps = 1
Iter=315, Fval=119.048937, Steps = 1
Iter=316, Fval=119.048937, Steps = 1
Iter=317, Fval=119.048937, Steps = 1
Iter=318, Fval=119.048937, Steps = 1
Iter=319, Fval=119.048937, Steps = 1
Iter=320, Fval=119.048937, Steps = 1
Iter=321, Fval=119.048937, Steps = 1
Iter=322, Fval=119.048937, Steps = 1
Iter=323, Fval=119.048937, Steps = 1
Iter=324, Fval=119.048937, Steps = 1
Iter=325, Fval=119.048937, Steps = 1
Iter=326, Fval=119.048937, Steps = 1
Iter=327, Fval=119.048937, Steps = 1
Iter=328, Fval=119.048937, Steps = 1
Iter=329, Fval=119.048937, Steps = 1
Iter=329, Fval=119.048937, Steps = 1
Optimization terminated: change in best function value less than options.TolFun.