{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "
24 April, 2017
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Exercise 2\n", "\n", "### by Anne Peter (anne.peter@uni-weimar.de)\n", "---" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Implementing the Travelling Salesman Problem" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 1. Travelling Salesman Problem" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Given a list of cities and the distances between each pair of cities, what is the shortest possible route that visits each city exactly once and returns to the origin city?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The TSP is an **NP-Hard Problem**. That does not necessarily mean any one instance of the problem will be hard to solve, it just means that we do not currently have an algorithm that can give us the guaranteed best solution for all problems in “polynomial time”. We can’t make predictions about how long it might take to find the best solution to the TSP from just looking at the data. We have no way of knowing how long a problem that is twice as large as one that took 2 minutes to solve will take.\n", "\n", "Although we might not be able to make predications about finding the best solution, we often only want a _good solution_ to the TSP. We aren’t always so worried if we find a route amongst 1000 cities that is only a few miles longer than the best solution – particularly if it would take an inordinate amount of computing time to get from the good solution we already have to the best solution." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " \n", "
\n", "
\n", " TSP can be modelled as an undirected weighted graph:\n", " - cities are the graph's vertices\n", " - paths are the graph's edges\n", " - a path's distance is the edge's weight" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Adjacency matrix** for the graph above:\n", "\n", "- undirected graphs are always symetric\n", "\n", "| | A | B | C | D |\n", "| :----: | :---: | :---: | :---: | :---: |\n", "| **A** | 0 | 20 | 42 | 35 |\n", "| **B** | 20 | 0 | 30 | 34 |\n", "| **C** | 42 | 30 | 0 | 12 |\n", "| **D** | 35 | 34 | 12 | 0 |" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 1.1 TSP in Python" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from itertools import permutations\n", "\n", "\n", "def distance(point1, point2):\n", " \"\"\"\n", " Returns the Euclidean distance of two points in the Cartesian Plane.\n", "\n", " >>> distance([3,4],[0,0])\n", " 5.0\n", " \n", " \"\"\"\n", " return ((point1 - point2)**2 + (point1 - point2)**2) ** 0.5" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The Euclidean distance between P1 = (x1, y1) and P2 = (x2, y2) is:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$d_{P_{1},P_{2}} = \\sqrt{(x_{1} - x_{2})^2 + (y_{1} - y_{2})^2}$" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "9.0\n" ] } ], "source": [ "print(distance([3,6],[7,6]) + distance([7,6],[12,6]))" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def total_distance(points):\n", " \"\"\"\n", " Returns the length of the path passing throught\n", " all the points in the given order.\n", "\n", " >>> total_distance([[1,2],[4,6]])\n", " 5.0\n", " >>> total_distance([[3,6],[7,6],[12,6]])\n", " 9.0\n", " \"\"\"\n", " return sum([distance(point, points[index + 1]) for index, point in enumerate(points[:-1])])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "- [:-1] means \"all elements of the sequence but the last\"\n", "- _enumerate_ enumerates all elements of the given sequence (here: points), start is 0" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "[(0, 'Spring'), (1, 'Summer'), (2, 'Fall'), (3, 'Winter')]" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# example for enumerate\n", "seasons = ['Spring', 'Summer', 'Fall', 'Winter']\n", "list(enumerate(seasons))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "- the for loop goes over every _index_ of the new enumarated sequence and every _point_ of the sequence\n", "- the method _distance_ computes the distance between each point and its successor, added together by _sum_" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def travelling_salesman(points, start = None):\n", " \"\"\"\n", " Finds the shortest route to visit all the cities by bruteforce.\n", " Time complexity is O(N!), so never use on long lists.\n", "\n", " >>> travelling_salesman([[0,0],[10,0],[6,0]])\n", " ([0, 0], [6, 0], [10, 0])\n", " >>> travelling_salesman([[0,0],[6,0],[2,3],[3,7],[0.5,9],[3,5],[9,1]])\n", " ([0, 0], [6, 0], [9, 1], [2, 3], [3, 5], [3, 7], [0.5, 9])\n", " \"\"\"\n", " if start is None:\n", " start = points\n", " return min([perm for perm in permutations(points) if perm == start], key = total_distance)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "- _permutations_ returns tuples with all possible orderings of _points_ and no repeated elements\n", "- returns the minimum of all possible tuples by the help of *total\\_distance*" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The minimum distance to visit all the following points: [[0, 0], [1, 5.7], [2, 3], [3, 7], [0.5, 9], [3, 5], [9, 1], [10, 5]]\n", " starting at [0, 0] is 25.90302275027582 and takes this route: ([0, 0], [2, 3], [3, 5], [1, 5.7], [0.5, 9], [3, 7], [10, 5], [9, 1]).\n" ] } ], "source": [ "def main():\n", " points = [[0, 0], [1, 5.7], [2, 3], [3, 7],\n", " [0.5, 9], [3, 5], [9, 1], [10, 5]]\n", " result = travelling_salesman(points)\n", " distance_result = total_distance(result)\n", " print(\"\"\"The minimum distance to visit all the following points: {}\n", " starting at {} is {} and takes this route: {}.\"\"\".format(\n", " points,\n", " points,\n", " distance_result,\n", " result,))\n", "\n", "\n", "if __name__ == \"__main__\":\n", " main()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 1.2 Solving TSP using Hill Climbing" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Recap: Hill Climbing" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Idea:\n", "- use only your local solution and evaluate your \n", "neighbors to find a better one\n", "- repeat this step until no better neighbor exists" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Pros:\n", "- requires few resources (current state and neighbors)\n", "- finds local optimum (global is possible)\n", "- useful if the search space is huge (even unlimited)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Cons:\n", "\n", "- is prone to get stuck at the top of local maximum and on plateaus\n", "- strongly depends on “good” initialization" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Hill Climbing Algorithm" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We will use the standard Python lists to represent a route/tour through a collection of cities. Each city will simply be assigned a number from 0 to N-1 (where N is the number of cities) and therefore our list of cities will be a list of unique numbers between 0 and N-1." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Adjacency Matrix " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We also need to specify a “distance matrix” that we can use to find out the distance between two cities. To generate a distance matrix for a set of x,y coordinates the following will do nicely:" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def cartesian_matrix(coords):\n", " '''\n", " Creates a distance matrix for the city coords using straight line distances\n", " computed by the Euclidean distance of two points in the Cartesian Plane.\n", " '''\n", " matrix = {}\n", " for i,(x1,y1) in enumerate(coords):\n", " for j,(x2,y2) in enumerate(coords):\n", " dx,dy = x1-x2,y1-y2\n", " dist = (dx*dx + dy*dy) ** 0.5\n", " matrix[i,j] = dist\n", " return matrix" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "- cartesian_matrix() takes a Python list of (x,y) tuples and outputs a Python dictionary that contains the distance between the distances between any two cities:" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "{(0, 0): 0.0, (0, 1): 1.0, (0, 2): 1.4142135623730951, (1, 0): 1.0, (1, 1): 0.0, (1, 2): 1.0, (2, 0): 1.4142135623730951, (2, 1): 1.0, (2, 2): 0.0}\n", "\n", "1.0\n" ] } ], "source": [ "test_matrix = cartesian_matrix([(0,0),(1,0),(1,1)])\n", "\n", "print(test_matrix)\n", "print()\n", "print(test_matrix[1,2])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Where matrix[1,2] gives the distance between city number 1 and city number 2. In our case this is the same as matrix[2,1], but for some TSP’s it may not be (for example if there is a one way street between cities that means we have to take an indirect route)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Read City Coodinates from File" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In addition to generating the distance matrix we will probably also want to read the city coordinates from a text file (one x,y per line):" ] }, { "cell_type": "code", "execution_count": 36, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def read_coords(coord_file):\n", " coords=[]\n", " for line in coord_file:\n", " x,y = line.strip().split(',')\n", " coords.append((float(x),float(y)))\n", " return coords\n", "\n", "\n", "coord_file = open('city100.txt', 'r') \n", "coords = read_coords(coord_file)\n", "\n", "matrix = cartesian_matrix(coords)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "That should be sufficient for generating distance matrices for now. On real world problems generating a distance matrix may be more complicated – you might need to take map data and calculate what the actual distance by road between any two cities is.\n", "This process can be done offline, before we start optimising our routes and is a subject for another time.\n", "\n", "Ok, so now we can read in a list of cities from a file and generate our distance matrix. What next? Well it would be good if we knew how long a route was:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Compute Total Distance" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# matrix is a distance matrix and tour is a list of cities (as integers)\n", "def tour_length(matrix, tour):\n", " '''Sum up the total length of the tour based on the distance matrix'''\n", " total = 0\n", " num_cities = len(list(tour))\n", " for i in list(range(num_cities)):\n", " j = (i+1) % num_cities\n", " city_i = tour[i]\n", " city_j = tour[j]\n", " total += matrix[city_i, city_j]\n", " return total" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Implementing the Tweak Operators\n", "\n", "We are going to implement the two tweak operators as generator functions, that will return (in a random order) all of the possible versions of a route that can be made in one step of the operator.
\n", "Generators are iterators, but you can only iterate over them once. It's because they do not store all the values in memory, they generate the values on the fly.
\n", "By using a generator function we can assess each different possibility and perhaps decide to not generate any more variations – which saves us the overhead of generating all of the combinations in one go.\n", "\n", "Both tweak operators rely on the following generator function:" ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import random\n", "\n", "def all_pairs(size, shuffle = random.shuffle):\n", " r1 = list(range(size))\n", " r2 = list(range(size))\n", " if shuffle:\n", " shuffle(r1)\n", " shuffle(r2)\n", " for i in r1:\n", " for j in r2:\n", " yield (i,j)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Which will generate all pairings of the numbers from 0 to size as (i,j) tuples in a random order (needs the random module to work).\n", "\n", "So each tweak operator then looks like:" ] }, { "cell_type": "code", "execution_count": 38, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[2, 1, 3]\n", "[3, 2, 1]\n", "[1, 3, 2]\n", "\n", "[3, 2, 1]\n", "[2, 1, 3]\n", "[2, 3, 1]\n", "[3, 1, 2]\n", "[1, 3, 2]\n" ] } ], "source": [ "from copy import deepcopy\n", "\n", "def swapped_cities(tour):\n", " '''Generator to create all possible variations where two cities have been swapped.'''\n", " for i,j in all_pairs(len(tour)):\n", " if i < j:\n", " copy = deepcopy(tour)\n", " copy[i],copy[j] = tour[j],tour[i]\n", " yield copy\n", "\n", "def reversed_sections(tour):\n", " '''\n", " Generator to return all possible variations where the section between two cities are swapped.\n", " It preserves entire sections of a route, yet still affects the ordering of multiple cities in one go.\n", " '''\n", " for i,j in all_pairs(len(tour)):\n", " if i != j:\n", " copy = deepcopy(tour)\n", " if i < j:\n", " copy[i:j+1] = reversed(tour[i:j+1])\n", " else:\n", " copy[i+1:] = reversed(tour[:j])\n", " copy[:j] = reversed(tour[i+1:])\n", " if copy != tour: # no point returning the same tour\n", " yield copy\n", "\n", "# usage\n", "for tour in swapped_cities([1,2,3]):\n", " print(tour)\n", "\n", "print()\n", " \n", "for tour in reversed_sections([1,2,3]):\n", " print(tour)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Getting Started" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To begin with the hill-climbing code we need two functions:\n", "\n", "- an initialisation function – that will return a random solution\n", "- an objective function – that will tell us how “good” a solution is\n", "\n", "For the TSP the initialisation function will just return a tour of the correct length that has the cities arranged in a random order.\n", "\n", "The objective function will return the negated length of a tour/solution. We do this because we want to maximise the objective function, while at the same time minimise the tour length.\n", "\n", "As the hill-climbing code won’t know specifically about the TSP we need to ensure that the initialisation function takes no arguments and returns a tour of the correct length and the objective function takes one argument (the solution tour) and returns the negated length.\n", "\n", "So assuming we have our city coordinates in a variable _coords_ and our distance matrix in _matrix_ we can define the objective function and initialisation functions using the function *init\\_random\\_tour*:" ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def init_random_tour(tour_length):\n", " tour = list(range(tour_length))\n", " random.shuffle(list(tour))\n", " return tour\n", "\n", "init_function = lambda: init_random_tour(len(coords))\n", "objective_function = lambda tour: -tour_length(matrix, tour) # note the negation" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Short Explanation of Lambda Functions" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "- creation of anonymous functions" ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "25\n", "25\n" ] } ], "source": [ "# normal function definition\n", "def f(x): return x**2\n", "\n", "# lambda funtion definition\n", "g = lambda x: x**2\n", "\n", "print(f(5))\n", "print(g(5))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "- lambda definition does not include a \"return\" statement\n", "- it always contains an expression which is returned\n", "- you can put a lambda definition anywhere a function is expected\n", "- you don't have to assign it to a variable" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Basic Hill-Climb" ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def hillclimb(init_function, move_operator, objective_function, max_evaluations):\n", " '''Hillclimb until either max_evaluations is reached or we are at a local optima.'''\n", " best = init_function()\n", " best_score = objective_function(best)\n", " \n", " num_evaluations = 1\n", " \n", " while num_evaluations < max_evaluations:\n", " # examine moves around our current position\n", " move_made = False\n", " for next in move_operator(best):\n", " if num_evaluations >= max_evaluations:\n", " break\n", " \n", " next_score = objective_function(next)\n", " num_evaluations += 1\n", " # see if this move is better than the current\n", " if next_score > best_score:\n", " best = next\n", " best_score = next_score\n", " move_made = True\n", " break # depth first search\n", " \n", " if not move_made:\n", " break # we couldn't find a better move (must be at a local maximum)\n", " \n", " return (num_evaluations,best_score,best)" ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "collapsed": true }, "outputs": [], "source": [ "from PIL import Image, ImageDraw, ImageFont\n", "\n", "def write_tour_to_img(coords, tour, title, img_file):\n", " padding = 20\n", " # shift all coords in a bit\n", " coords = [(x+padding,y+padding) for (x,y) in coords]\n", " maxx, maxy = 0,0\n", " for x,y in coords:\n", " maxx = max(x,maxx)\n", " maxy = max(y,maxy)\n", " maxx += padding\n", " maxy += padding\n", " img = Image.new(\"RGB\",(int(maxx), int(maxy)), color=(255,255,255))\n", " \n", " font=ImageFont.load_default()\n", " d=ImageDraw.Draw(img);\n", " num_cities = len(tour)\n", " for i in range(num_cities):\n", " j = (i+1) % num_cities\n", " city_i = tour[i]\n", " city_j = tour[j]\n", " x1,y1 = coords[city_i]\n", " x2,y2 = coords[city_j]\n", " d.line((int(x1), int(y1), int(x2), int(y2)), fill=(0,0,0))\n", " d.text((int(x1)+7, int(y1)-5), str(i), font=font, fill=(32,32,32))\n", " \n", " \n", " for x,y in coords:\n", " x,y = int(x), int(y)\n", " d.ellipse((x-5, y-5, x+5, y+5), outline=(0,0,0), fill=(196,196,196))\n", " \n", " d.text((1,1), title, font=font, fill=(0,0,0))\n", " \n", " del d\n", " img.save(img_file, \"PNG\")" ] }, { "cell_type": "code", "execution_count": 18, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "-16035.517246233385\n", "[71, 98, 86, 30, 68, 44, 0, 81, 17, 29, 82, 70, 85, 77, 21, 20, 5, 66, 69, 99, 9, 14, 34, 78, 43, 87, 2, 1, 41, 96, 16, 54, 94, 47, 39, 31, 8, 57, 89, 92, 28, 25, 75, 52, 49, 64, 63, 80, 11, 12, 32, 61, 18, 22, 67, 23, 10, 56, 59, 74, 62, 90, 3, 97, 15, 26, 51, 48, 65, 4, 83, 36, 24, 27, 91, 46, 50, 35, 73, 60, 37, 6, 7, 33, 79, 88, 38, 19, 84, 45, 76, 58, 72, 42, 40, 93, 95, 55, 53, 13]\n" ] } ], "source": [ "# move_operator = swapped_cities\n", "move_operator = reversed_sections # better\n", "\n", "max_evaluations = 500\n", "\n", "num_evaluations, best_score, best = hillclimb(init_function, move_operator, objective_function, max_evaluations)\n", "\n", "print(best_score)\n", "print(best)" ] }, { "cell_type": "code", "execution_count": 19, "metadata": { "collapsed": false }, "outputs": [], "source": [ "filename = \"test500.PNG\"\n", "\n", "write_tour_to_img(coords, best, filename, open(filename, 'ab'))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " " ] }, { "cell_type": "code", "execution_count": 20, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "-5746.654404651131\n", "[69, 60, 87, 33, 24, 70, 27, 91, 40, 25, 2, 93, 12, 64, 45, 37, 6, 75, 4, 83, 62, 53, 71, 32, 15, 94, 55, 76, 95, 72, 58, 13, 52, 65, 48, 84, 54, 41, 16, 61, 36, 18, 34, 21, 10, 88, 67, 22, 90, 1, 63, 23, 96, 19, 31, 20, 79, 8, 14, 38, 43, 74, 49, 78, 26, 3, 51, 97, 47, 11, 80, 85, 28, 7, 46, 89, 29, 17, 98, 42, 81, 57, 59, 82, 99, 92, 5, 68, 0, 35, 56, 39, 86, 50, 9, 30, 73, 44, 66, 77]\n" ] } ], "source": [ "max_evaluations = 5000\n", "\n", "num_evaluations, best_score, best = hillclimb(init_function, move_operator, objective_function, max_evaluations)\n", "\n", "print(best_score)\n", "print(best)" ] }, { "cell_type": "code", "execution_count": 21, "metadata": { "collapsed": true }, "outputs": [], "source": [ "filename = \"test5000.PNG\"\n", "\n", "write_tour_to_img(coords, best, filename, open(filename, 'ab'))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " " ] }, { "cell_type": "code", "execution_count": 22, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "-4160.9456573136\n", "[68, 0, 35, 39, 56, 86, 73, 21, 34, 18, 22, 67, 88, 10, 96, 38, 23, 63, 1, 41, 61, 36, 16, 26, 54, 84, 3, 78, 90, 97, 51, 48, 65, 47, 52, 71, 53, 15, 32, 94, 55, 76, 95, 72, 58, 37, 45, 6, 93, 12, 2, 40, 91, 27, 17, 29, 89, 24, 70, 46, 98, 25, 42, 4, 28, 85, 80, 75, 62, 64, 13, 11, 49, 83, 74, 43, 8, 14, 19, 31, 30, 9, 50, 99, 20, 79, 59, 57, 7, 81, 33, 87, 82, 60, 69, 77, 44, 66, 92, 5]\n" ] } ], "source": [ "max_evaluations = 50000\n", "\n", "num_evaluations, best_score, best = hillclimb(init_function, move_operator, objective_function, max_evaluations)\n", "\n", "print(best_score)\n", "print(best)" ] }, { "cell_type": "code", "execution_count": 23, "metadata": { "collapsed": true }, "outputs": [], "source": [ "filename = \"test50000.PNG\"\n", "\n", "write_tour_to_img(coords, best, filename, open(filename, 'ab'))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 1.3 Solving TSP using Simulated Annealing" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As we discussed before, it is possible for an algorithm to find a solution that is “locally optimal”, but not necessarily “globally optimal”. That is to say we may find ourselves with a solution that is the best thing nearby, but it might not be the best thing. This happens with hill-climbing, because when we are offered the choice between two solutions we always take the best solution. The algorithm is greedy and short sighted." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "So instead we could try occasionally choosing something that’s worse. By doing that the algorithm can go “downhill” sometimes and hopefully reach new areas of the solution landscape." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Recap: Simulated Annealing" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Simulated annealing is essentially hill-climbing, but with the ability to go downhill (sometimes)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "It introduces a temperature variable. When the temperature is high a worse solution will have a higher chance of being chosen. It work’s like this:\n", "\n", "1. pick an initial solution\n", "2. set an initial temperature\n", "3. choose the next neighbour of the current solution:\n", " - if the neighbour is better, make that neighbour the current solution\n", " - if the neighbour is worse, probabilistically make this neighbour the current solution, based on the current\n", " temperature and how much worse the neighbour is\n", "4. decrease the temperature slightly\n", "5. go to 3.\n", "\n", "By slowly cooling the temperature we become less likely to choose worse solutions over time. Initially we are able to make some pretty big jumps around the solution landscape. By the end of a run we’ll be jumping around less and less. In fact if we lower the temperature enough we end up with plain old hill-climbing." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Simulated Annealing Algorithm" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Cooling Down" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Temperature is a key part of simulated annealing. How we lower the temperature over time is therefore very important. There are a couple of possible approaches, but I’ll show the one outlined by Kirkpatrick et al:" ] }, { "cell_type": "code", "execution_count": 24, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# you should think of a better cooling function\n", "def kirkpatrick_cooling(start_temp, alpha):\n", " T = start_temp\n", " while True:\n", " yield T\n", " T = alpha * T" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The keyword _yield_ is used like return, except the function will return a generator.\n", "Generators are iterators, but you can only iterate over them once. It's because they do not store all the values in memory, they generate the values on the fly" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This generator function takes an initial start temperature (start_temp) and returns a series of temperatures that are alpha times the size, where alpha is less than one. So we end up with a temperature that drops off quite quickly and then slowly decreases to practically nothing." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Probabilistically Choosing a Neighbor" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Below is the Python code to decide if what probability we will assign to moving from a solution with a score of prev_score to a solution with a value of next_score at the current temperature." ] }, { "cell_type": "code", "execution_count": 25, "metadata": { "collapsed": true }, "outputs": [], "source": [ "import math\n", "\n", "def P(prev_score, next_score, temperature):\n", " if next_score > prev_score:\n", " return 1.0\n", " elif temperature == 0:\n", " return 0.0\n", " else:\n", " return math.exp(-abs(next_score - prev_score)/temperature)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To keep later logic simpler it's returning 1.0 if next_score is better – so we’ll always choose better solutions.\n", "\n", "When the prev_score is worse we create a probability based on the difference between prev_score and next_score scaled by the current temperature. If we chart the probabilities versus the difference in scores we get (with a temperature of 1.0):" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As can be seen, for small differences (relative to the current temperature) we will have a high probability. This then tails off very quickly, so solutions that are much worse are increasingly less likely to be chosen.\n", "\n", "The net-effect being that solutions that are only a little bit worse are still fairly likely to be chosen. Much worse solutions may still be chosen, but it’s much less likely." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Remember the Best Solution" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "One other minor, but key, implementation detail is saving the best solution we find during the annealing process.\n", "\n", "During hill-climbing the current solution was always the best solution found, but simulated annealing will deliberately accept worse solutions at times. So we need to make sure we don’t just throw away the best we see. To avoid complicating the algorithm itself with extra checks of scores etc.\n", "\n", "I am going to use a class to wrap the objective function. I’ll override the \\_\\_call\\_\\_ method of the class, so that I can use the instance of the class like a function – in place of the normal objective function:" ] }, { "cell_type": "code", "execution_count": 26, "metadata": { "collapsed": true }, "outputs": [], "source": [ "class ObjectiveFunction:\n", " '''\n", " Class to wrap an objective function and \n", " keep track of the best solution evaluated.\n", " '''\n", " def __init__(self, objective_function):\n", " self.objective_function = objective_function\n", " self.best = None\n", " self.best_score = None\n", " \n", " def __call__(self, solution):\n", " score = self.objective_function(solution)\n", " if self.best is None or score > self.best_score:\n", " self.best_score = score\n", " self.best = solution\n", " return score" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can then access then best and best_score fields when we have finished our annealing." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Simulated Annealing itself" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The code below represents the simulated annealing algorithm. In many respects it is pretty similar to hill-climbing, but we are also concerned with a current temperature and we have introduced a probabilistic element to choosing the next solution." ] }, { "cell_type": "code", "execution_count": 27, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def anneal(init_function, move_operator, objective_function, max_evaluations, start_temp,alpha):\n", " \n", " # wrap the objective function (so we record the best)\n", " objective_function = ObjectiveFunction(objective_function)\n", " \n", " current = init_function()\n", " current_score = objective_function(current)\n", " num_evaluations = 1\n", " \n", " cooling_schedule = kirkpatrick_cooling(start_temp, alpha)\n", " \n", " for temperature in cooling_schedule:\n", " done = False\n", " # examine moves around our current position\n", " for next in move_operator(current):\n", " if num_evaluations >= max_evaluations:\n", " done = True\n", " break\n", " \n", " next_score = objective_function(next)\n", " num_evaluations += 1\n", " \n", " # probablistically accept this solution\n", " # always accepting better solutions\n", " p = P(current_score, next_score, temperature)\n", " if random.random() < p:\n", " current = next\n", " current_score = next_score\n", " break\n", " # see if completely finished\n", " if done: break\n", " \n", " best_score = objective_function.best_score\n", " best = objective_function.best\n", " \n", " return (num_evaluations, best_score, best)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The parameters are much the same as hill-climbing, but there are two extra specific to simulated annealing:\n", "\n", "- init_function - the function used to create our initial solution\n", "- move_operator - the function we use to iterate over all possible \"moves\" for a given solution\n", "- objective_function - used to assign a numerical score to a solution - how \"good\" the solution is\n", "- max_evaluations - used to limit how much search we will perform (how many times we'll call the objective_function)\n", "- start_temp - the initial starting temperature for annealing\n", "- alpha - should be less than one. controls how quickly the temperature reduces\n", "\n", "We are also only reducing the temperature after either accepting a new solution or evaluating all neighbours without choosing any of them. This is done so that temperature will only decrease as we start accepting moves. As that will be less frequent than just evaluating moves we cooling will happen at a slower pace. If we are accepting lots of moves then this will drop the temperature quite quickly. If we are not accepting many moves the temperature will stay steadier - maintaining the likelihood of accepting other \"worse\" moves. That latter point is useful, as if we are starting to get stuck on a local maximum the temperature won't decrease - hopefully helping us get unstuck." ] }, { "cell_type": "code", "execution_count": 28, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "-15428.808243629654\n", "[99, 6, 15, 16, 84, 41, 38, 90, 53, 80, 58, 89, 87, 14, 8, 64, 12, 76, 81, 91, 27, 33, 70, 17, 59, 46, 54, 51, 71, 37, 55, 66, 77, 82, 60, 30, 7, 92, 29, 40, 98, 24, 50, 67, 42, 47, 48, 72, 45, 2, 49, 36, 23, 22, 21, 9, 56, 39, 35, 44, 68, 18, 19, 88, 97, 63, 1, 26, 78, 25, 93, 95, 94, 65, 3, 4, 85, 83, 75, 62, 13, 10, 74, 52, 11, 32, 57, 43, 96, 34, 79, 20, 31, 28, 86, 61, 73, 0, 5, 69]\n" ] } ], "source": [ "# move_operator = swapped_cities\n", "move_operator = reversed_sections # better\n", "\n", "max_iterations = 500\n", "start_temp = 10\n", "alpha = 0.9999\n", "\n", "iterations, best_score, best = anneal(init_function, move_operator, objective_function, max_iterations, start_temp, alpha)\n", "\n", "print(best_score)\n", "print(best)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Comparison to Hill Climbing" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
EvaluationsAlgorithmAverages.d.WorstBest
50000Hill Climbing-4228.50126.45-4627.07-3942.03
50000SA (start_temp=10, alpha=0.9999)-4145.6996.56-4422.04-3924.34
100000Hill Climbing-4154.2590.60-4513.11-3946.65
100000SA (start_temp=10, alpha=0.99995)-4077.4071.72-4294.97-3907.19