شبیه سازی تبرید (Simulated Annealing) – به زبان ساده
در ریاضیات و علوم کامپیوتر، «مساله بهینهسازی» (Optimization Problem)، در واقع مساله پیدا کردن «بهترین» راه حل، از میان کلیه راه حلهای «ممکن» برای مساله است. نوع خاصی از مسائل بهینهسازی وجود دارند که در آنها، به دلیل زیاد شدن تعداد «اشیا» (Objects) (منظور تعداد مشاهدات است)، مدیریت کردن آنها با استفاده از «روشهای ترکیبیاتی» (Combinatorial Methods) امکانپذیر نیست. یک مثال شناخته شده از چنین مسائلی، «مساله فروشنده دورهگرد» (Traveling Salesman Problem) است که از دسته مسائل «انپی کامل» (NP-Complete که در آن، NP سرنامی برای عبارت Non-Deterministic Polynomial است) محسوب میشود.
برای حل چنین مسائلی، یک الگوریتم بسیار کاربردی به نام «شبیه سازی تبرید» (Simulated Annealing) (به آن «تبرید شبیه سازی شده»، «بازپخت شبیه سازی شده» و یا به طور خلاصه، SA نیز گفته میشود) وجود دارد. از الگوریتم شبیه سازی تبرید اغلب برای تخمین «بهینه سراسری» (Global Optimization) در مسائل بهینهسازی که فضای جستجوی آنها بزرگ است، استفاده میشود.
همچنین، این الگوریتم معمولا هنگامی که فضای جستجو گسسته باشد، مورد استفاده قرار میگیرد (همه سفرهایی که از مجموعه مشخصی از شهرها طی آنها عبور میشود). در این مطلب، به الگوریتم شبیه سازی تبرید و کاربرد آن در حل مسائل بهینهسازی پرداخته میشود، اما پیش از پرداختن به الگوریتم، مفهوم فرایند بازپخت (انیلینگ | Annealing) در «متالورژی» (Metallurgy) که الگوریتم تبرید شبیه سازی شده در واقع تقلیدی از آن است، بیان خواهد شد. سپس، الگوریتم شبیه سازی تبرید شرح داده میشود و در نهایت به حل مساله فروشنده دورهگرد با آن پرداخته میشود. کد مربوط به پیادهسازی الگوریتم شبیه سازی تبرید در زبان برنامهنویسی «پایتون» (Python) و «جاوا» (Java) در انتخای این مطلب، ارائه شده است.
فرایند بازپخت یا انیلینگ
در متالورژی و علم مواد، بازپخت یا انیلینگ، به عملیات حرارتی گفته میشود که طی آن مشخصات فیزیکی و گاهی، مشخصات شیمیایی ماده تغییر میکند. این کار با هدف افزایش شکلپذیری و کاهش سختی ماده انجام میشود. طی این فرایند، ابتدا فلز گرم شده، سپس در یک درجه حرارت خاص نگه داشته و در نهایت، به تدریج سرد میشود.
با گرم کردن فلز، مولکولها آزادانه به هر سوی حرکت میکنند. با سرد کردن تدریجی ماده، این آزادی کاهش پیدا میکند. اگر فرایند سرد کردن به اندازه کافی کند باشد که بتوان اطمینان حاصل کرد فلز در هر مرحله در تعادل ترمودینامیکی قرار دارد، میتوان مطمئن شد که انرژی گرمایی به طور یکنواخت در جسم توزیع شده و بهترین ساختار بلوری در آن وجود دارد که متقارن و مقاوم است. در الگوریتم شبیهسازی تبرید، از فرایند مذکور الگوبرداری شده است.
مساله فروشنده دورهگرد
در مساله «فروشنده دورهگرد» (Travelling Salesman Problem | TSP)، پرسش زیر مطرح میشود:
لیستی از شهرها و فاصله بین هر جفت از شهرها موجود است. کوتاهترین مسیر ممکنی که با استفاده از آن بتوان از یک نقطه شروع کرد و پس از عبور از همه شهرها به نقطه اول بازگشت، کدام است؟
این مساله، یکی از مسائل «انپی-سخت» (NP-Hard) در بهینهسازی کامپیوتری محسوب میشود که در حوزه «تحقیق در عملیات» (Operational Research | OR) و «علوم نظری رایانه» (Theoretical Computer Science | TCS) حائز اهمیت است. در «نظریه پیچیدگی محاسباتی» (Theory of Computational Complexity)، نسخه تصمیم مساله فروشنده دورهگرد (که در آن یک طول L داده شده، وظیفه تصمیمگیری پیرامون آن است که آیا گراف دارای دوری کمتر از L است یا خیر) از جمله مسائل انپی-کامل محسوب میشود.
بنابراین، این امکان وجود دارد که بدترین زمان اجرا برای هر الگوریتم برای TSP، با افزایش تعداد شهرها به صورت «فوق چند جملهای» (Superpolynomially) افزایش پیدا کند (اما نه بیشتر از نمایی).
این مساله برای اولین بار در سال ۱۹۳۰ فرموله شد و یکی از مسائلی است که بیشترین مطالعات در حوزه بهینهسازی روی آن انجام شده است. از این مساله، به عنوان «بنچمارک» (Benchmark) برای بسیاری از روشهای بهینهسازی استفاده میشود. هرچند حل مساله فروشنده دورهگرد حتی به لحاظ کامپیوتری نیز دشوار است، تعداد زیادی از «روشهای تکاملی» (Evolutionary Algorithms) و «الگوریتمهای دقیق» (Exact Algorithm) برای حل این مساله شناخته شده هستند که میتوانند TSP را با تعداد هزاران و حتی میلیونها شهر حل کنند. مساله فروشنده دورهگرد را میتوان به صورت یک مساله گراف و یا «برنامهنویسی خطی صحیح» (Integer Linear Programming) فرموله کرد. تعداد کل پاسخهای ممکن برای مساله برابر با «$$\frac{1}{2}\left(n-1\right)!$$» برای n>2 است که در آن، n تعداد کل شهرها است. در واقع، این عدد تعداد دورهای همیلتونی در یک گراف کامل با n راس را نشان میدهد.
الگوریتم شبیه سازی تبرید
همانطور که پیشتر بیان شد، در الگوریتم شبیه سازی تبرید (یا تبرید شبیه سازی شده) از فرایند بازپخت که از مباحث رشته متالورژی و مواد محسوب میشود، الگو گرفته شده است. انتخاب نام شبیهسازی تبرید برای این الگوریتم، ریشه در فرایند دارد که از آن تقلید میکند. در بهینهسازی نیز مانند فرایند انیلینگ، آنچه در بخش پیشین پیرامون بازپخت مواد بیان شد، برای حل مسائل قابل انجام است.
یعنی در واقع، جوابهای یک مساله به خوبی گرم میشوند و با نوسانات زیادی تغییر میکنند؛ سپس، به تدریج دامنه تغییرات کم میشود و در واقع یک سری شیار به سمت جواب بهینه ساخته میشوند. الگوریتم شبیه سازی تبرید برای اولین بار در سال ۱۹۸۳، توسط «کریکپاتریک» (Kirkpatrick) و همکاران معرفی شد. شایان ذکر است، الگوریتم شبیه سازی تبرید از جمله الگوریتمهای فراابتکاری (فراتکاملی | فرااکتشافی | Metaheuristic) محسوب میشود. در الگوریتم شبیه سازی تبرید، از روش احتمالاتی برای حل مساله بهینهسازی استفاده میشود. علاقمندان به مشاهده ویدئوی آموزشی این الگوریتم به زبان فارسی، میتوانند از «آموزش شبیهسازی تبرید یا Simulated Annealing در متلب» بهرهمند شوند.
در الگوریتم شبیه سازی تبرید (SA)، نقطه s یک حالت از سیستم فیزیکی محسوب میشود و تابع (E(s مشابه با انرژی داخلی سیستم در حالت s است. هدف آن است که با شروع سیستم از یک حالت اولیه دلخواه (یک s0 دلخواه)، به حالتی رسیده شود (Sn) که تابع (E(s در آن کمینه است. در واقع، با شروع از یک حالت دلخواه از سیستم فیزیکی، به حالتی رسیده میشود که انرژی داخلی سیستم در آن حالت کمینه است (سیستم کمترین انرژی را در آن حالت خواهد داشت). برای انجام این کار، الگوریتم از یک نقطه دلخواه آغاز به کار و سپس، یک حالت همسایه را انتخاب میکند. پس از آن، به طور احتمالی تصمیم میگیرد که در حالت کنونی بماند و یا به حالت همسایه جا به جا شود. مجموع این جا به جاییهای احتمالی، سیستم را به سوی حالتی با انرژی داخلی کمتر هدایت میکند. این کار تا جایی انجام میشود که سیستم به یک حالت عقلانی برسد یا اینکه میزان محاسبات، از یک آستانه مشخص بیشتر شود. برای درک بهتر وجه تمایز الگوریتم شبیه سازی تبرید با برخی از دیگر روشهای ابتکاری مانند تپهنوری، مثال زیر قابل توجه است.
مساله فروشنده دورهگرد را میتوان به عنوان مثالی در نظر گرفت که الگوریتم شبیه سازی تبرید در حل آن کاربرد دارد. همانطور که پیشتر بیان شد، در این مساله، فروشنده باید از تعداد زیادی از شهرها در حالی عبور کند که مسافت کل پیموده شده کمینه باشد. اگر فروشنده از یک مسیر تصادفی حرکت خود را آغاز کند، بعدا میتواند با این امید که در هر تغییر شهر مسافت پیموده شده را کاهش دهد، به ترتیب از کلیه شهرها عبور کند.
چالش اصلی در حل مساله فروشنده دوره گرد با استفاده از روشی مانند الگوریتم «تپهنوردی» (Hill Climbing)، از آنجا نشات میگیرد که این الگوریتم با جا به جایی بین همسایهها معمولا به حالت کمینه میرسد و در همان نقطه متوقف میشود (در واقع در حالتی که نسبت به دیگر همسایگیهای بررسی شده کمینه باشد متوقف میشود؛ در صورتی که امکان دارد یک کمینه سراسری وجود داشته باشد). در واقع، الگوریتم تپهنوردی «کمینه محلی» (Local Minimum) را معمولا به سرعت پیدا میکند، اما ممکن است در همان جا متوقف شود و بنابراین از آنجا به «کمینه سراسری» (Global Minimum) نمیرسد. این در حالی است که شبیهسازی تبرید میتواند راه حل خیلی خوبی را حتی در حضور «نویز» (Noise) برای چنین مسالهای پیدا کند.
راهکار الگوریتم شبیه سازی تبرید برای غلبه بر این مشکل و بهبود استراتژی مذکور، استفاده از دو ترفند است. ترفند اول، از «الگوریتم متروپولیس» (Metropolis Algorithm) گرفته شده که در سال ۱۹۵۳ توسط متروپولیس و همکاران اون کشف شده است. در این الگوریتم، گاهی مسیرهای کوتاه پذیرفته نمیشوند، زیرا این عدم پذیرش منجر به آن میشود که الگوریتم «فضای راه حل ممکن» بزرگتری را اکتشاف کند. ترفند دوم، مجددا با تقلید از فرایند بازپخت فلز و کاهش دما به دمای پایینتر اتفاق میافتد. پس از آنکه بررسیها و جا به جاییهای زیادی انجام شد و مشاهده شد که تابع (E(s (این تابع در واقع همان تابع هزینه یا Cost Function است) به آرامی کاهش پیدا میکند (دما کاهش پیدا میکند)، اندازه انجام جا به جاییهای «بد» کنترل میشود. پس از چندین بار کاهش دما به مقدار کمتر، فرایند «فرونشانی» (Quench) با پذیرش جا به جاییهای خوب به منظور پیدا کردن کمینه محلی تابع هزینه اتفاق میافتد. «زمانبندیهای بازپخت» (Annealing Schedules) متعددی برای کاهش درجه حرارت وجود دارد، اما نتایج به طور کلی خیلی به جزئیات حساس نیستند.
اکنون، الگوریتم شبیه سازی تبرید به صورت ریاضیاتی و دقیق ارائه میشود. احتمال انتقال از یک حالت فعلی، مثلا s به حالت کاندید جدید مانند $$s\prime$$ با یک تابع احتمال پذیرش $$P(e,e\prime,T)$$ مشخص میشود که در آن، (e=E(s و $$e\prime=E(s\prime)$$ است (چنانچه پیشتر بیان شد، تابع E در فضای حالت نشانگر انرژی داخلی سیستم و T نشانگر دما است. دمای T بر اساس زمان تغییر میکند). بنابراین، با توجه به آنکه پیشتر گفته شد هدف آن است که انرژی سیستم کمینه باشد، پس حالتی که (منظور s است) در آن (E(s کمتر باشد، بهتر از حالتی است که در آن مقدار (E(s بیشتر باشد. نکته قابل توجه در الگوریتم تبرید شبیه سازی شده آن است که تابع احتمال P همواره و حتی در شرایطی که e کوچکتر از $$e\prime$$ است، باید مثبت باشد. این خصوصیت موجب میشود الگوریتم در بهینه محلی که نسبت به بهینه سراسری «بدتر» است، متوقف نشود. در واقع، در صورتی که $$s\prime$$ بهتر از s $$E(s)\ge E(s\prime)$$ باشد، $$s\prime$$ پذیرفته میشود و اگر $$s\prime$$ بدتر از s باشد و موجب $$E(s)< E(s\prime)$$ شود، $$s\prime$$ با یک احتمال پذیرفته میشود. تابع احتمال P به صورت زیر است:
$$p=e^{-\frac{E(s\prime)-E(s)}{T}}=e^{-\frac{\Delta f}{T}}$$
در صورت کاهش دمای T و میل کردن آن به سمت صفر، احتمال P نیز باید کاهش پیدا کند و یا به صفر و یا یک عدد مثبت میل کند. در این شرایط، P هنگامی به سمت صفر میل میکند که $$e<e\prime$$ باشد و در صورتی به یک عدد مثبت میل میکند که $$e\prime<e$$ باشد. همانطور که مشهود است، دما در کنترل تغییرات سیستم نقش اساسی دارد. همانطور که پیشتر هم بیان شد، دما در شبیهسازی به صورت تدریجی کاهش پیدا میکند. بنابراین، الگوریتم از $$T=\infty$$ و به عبارتی، از یک دمای بسیار بزرگ کار خود را آغاز میکند و در هر مرحله با توجه به یک زمانبندی تبرید از پیش مشخص شده، کاهش پیدا میکند. زمانبندی تبرید با توجه به این موضوع انجام میشود که اگر منابع مورد استفاده (مثلا میزان محاسبات) به پایان برسند، زمان انجام فرایند نیز به پایان میرسد. از همین رو، الگوریتم ابتدا در فضای بزرگی از راه حلها، صرف نظر از انرژی داخلی سیستم به دنبال پاسخ میگردد و به تدریج، به سمت مناطق دارای انرژی کمتر جا به جا میشود. این منطقه، به تدریج کوچکتر میشود و این کار، تا زمانی که بهینه سراسری یافته شود، ادامه پیدا میکند. برای یک مساله دارای فضای متناهی، افزایش زمان موجب میشود تا با احتمال یک، الگوریتم در نقطه بهینه سراسری متوقف شود؛ هر چند زمان لازم برای انجام این کار، بیشتر از زمان لازم برای جستجو در کل فضا است و بنابراین، این موضوع در عمل کاربردی ندارد.
فلوچارت الگوریتم شبیه سازی تبرید
حل مساله فروشنده دورهگرد با شبیه سازی تبرید
برنامهای که با بهرهگیری از آن میتوان مساله فروشنده دورهگرد را با الگوریتم شبیه سازی تبرید حل کرد، از این مسیر [+] در دسترس است. در ادامه، انیمیشن مربوط به حل مساله فروشنده دورهگرد با استفاده از این برنامه برای ۴۸ تا از مراکز استانهای ایالات متحده آمریکا، قابل مشاهده است.
شبه کد الگوریتم تبرید شبیه سازی شده
1Let s = s0
2For k = 0 through kmax (exclusive):
3T ← temperature( kmax/(k+1) )
4Pick a random neighbour, snew ← neighbour(s)
5If P(E(s), E(snew), T) ≥ random(0, 1):
6s ← snew
7Output: the final state s
کد الگوریتم تبرید شبیه سازی شده در جاوا
1// Java program to implement Simulated Annealing
2import java.util.*;
3
4public class SimulatedAnnealing {
5
6 // Initial and final temperature
7 public static double T = 1;
8
9 // Simulated Annealing parameters
10
11 // Temperature at which iteration terminates
12 static final double Tmin = .0001;
13
14 // Decrease in temperature
15 static final double alpha = 0.9;
16
17 // Number of iterations of annealing
18 // before decreasing temperature
19 static final int numIterations = 100;
20
21 // Locational parameters
22
23 // Target array is discretized as M*N grid
24 static final int M = 5, N = 5;
25
26 // Number of objects desired
27 static final int k = 5;
28
29
30 public static void main(String[] args) {
31
32 // Problem: place k objects in an MxN target
33 // plane yielding minimal cost according to
34 // defined objective function
35
36 // Set of all possible candidate locations
37 String[][] sourceArray = new String[M][N];
38
39 // Global minimum
40 Solution min = new Solution(Double.MAX_VALUE, null);
41
42 // Generates random initial candidate solution
43 // before annealing process
44 Solution currentSol = genRandSol();
45
46 // Continues annealing until reaching minimum
47 // temprature
48 while (T > Tmin) {
49 for (int i=0;i<numIterations;i++){
50
51 // Reassigns global minimum accordingly
52 if (currentSol.CVRMSE < min.CVRMSE){
53 min = currentSol;
54 }
55
56 Solution newSol = neighbor(currentSol);
57 double ap = Math.pow(Math.E,
58 (currentSol.CVRMSE - newSol.CVRMSE)/T);
59 if (ap > Math.random())
60 currentSol = newSol;
61 }
62
63 T *= alpha; // Decreases T, cooling phase
64 }
65
66 //Returns minimum value based on optimiation
67 System.out.println(min.CVRMSE+"\n\n");
68
69 for(String[] row:sourceArray) Arrays.fill(row, "X");
70
71 // Displays
72 for (int object:min.config) {
73 int[] coord = indexToPoints(object);
74 sourceArray[coord[0]][coord[1]] = "-";
75 }
76
77 // Displays optimal location
78 for (String[] row:sourceArray)
79 System.out.println(Arrays.toString(row));
80
81 }
82
83 // Given current configuration, returns "neighboring"
84 // configuration (i.e. very similar)
85 // integer of k points each in range [0, n)
86 /* Different neighbor selection strategies:
87 * Move all points 0 or 1 units in a random direction
88 * Shift input elements randomly
89 * Swap random elements in input sequence
90 * Permute input sequence
91 * Partition input sequence into a random number
92 of segments and permute segments */
93 public static Solution neighbor(Solution currentSol){
94
95 // Slight perturbation to the current solution
96 // to avoid getting stuck in local minimas
97
98 // Returning for the sake of compilation
99 return currentSol;
100
101 }
102
103 // Generates random solution via modified Fisher-Yates
104 // shuffle for first k elements
105 // Pseudorandomly selects k integers from the interval
106 // [0, n-1]
107 public static Solution genRandSol(){
108
109 // Instantiating for the sake of compilation
110 int[] a = {1, 2, 3, 4, 5};
111
112 // Returning for the sake of compilation
113 return new Solution(-1, a);
114 }
115
116
117 // Complexity is O(M*N*k), asymptotically tight
118 public static double cost(int[] inputConfiguration){
119
120 // Given specific configuration, return object
121 // solution with assigned cost
122 return -1; //Returning for the sake of compilation
123 }
124
125 // Mapping from [0, M*N] --> [0,M]x[0,N]
126 public static int[] indexToPoints(int index){
127 int[] points = {index%M, index/M};
128 return points;
129 }
130
131 // Class solution, bundling configuration with error
132 static class Solution {
133
134 // function value of instance of solution;
135 // using coefficient of variance root mean
136 // squared error
137 public double CVRMSE;
138
139 public int[] config; // Configuration array
140 public Solution(double CVRMSE, int[] configuration) {
141 this.CVRMSE = CVRMSE;
142 config = configuration;
143 }
144 }
145}
کد الگوریتم تبرید شبیه سازی شده در پایتون
1#!/usr/bin/env python
2# -*- coding: utf-8 -*-
3"""
4Module for implementing Simulated Annealing algorithm for single objective optimization.
5author := "Daniel Barker"
6"""
7
8import random
9import math
10import numpy as np
11import pandas as pd
12import plotly.plotly as py
13import plotly.graph_objs as go
14
15
16# Formatting variables
17HEADER = '++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++'
18SPACER = '===================================================================='
19FOOTER = '********************************************************************'
20
21
22class SimulatedAnnealing(object):
23 """ Class for performing SimulatedAnnealing optimization algorithm."""
24
25 def __init__(self, start_temp, stop_temp, alpha, num_sims):
26 """
27 Constructor Method.
28 """
29 self.start_temp = start_temp
30 self.stop_temp = stop_temp
31 self.alpha = alpha
32 self.num_sims = num_sims
33
34 def _random_int(self, lower_bound, upper_bound, step):
35 return random.randrange(lower_bound, upper_bound, step)
36
37 def _random_float(self, lower_bound, upper_bound, sig_dig):
38 return round(random.uniform(lower_bound, upper_bound), sig_dig)
39
40 def _get_neighbor_int(self, current_value, distance, step, lower_bound, upper_bound):
41 """ Return a neighboring random integer value from some current value.
42 :param distance: defines neighborhood as current value +- distance
43 :param step: step size range within neighborhood
44 :param lower_bound: hard cutoff lower bound that overrides distances lower than this threshold
45 :param upper_bound: hard cutoff upper bound that overrides distances higher than this threshold
46 :return: random value within neighborhood of current input value
47 """
48 upper_neighbor = current_value + distance
49 if upper_neighbor > upper_bound:
50 upper_neighbor = upper_bound
51 else:
52 upper_neighbor = upper_neighbor
53
54 lower_neighbor = current_value - distance
55 if lower_neighbor < lower_bound:
56 lower_neighbor = lower_bound
57 else:
58 lower_neighbor = lower_neighbor
59 neighbor = random.randrange(lower_neighbor, upper_neighbor, step)
60 return neighbor
61
62 def _get_neighbor_float(self, current_value, distance, lower_bound, upper_bound, sig_dig):
63 """ Return a neighboring random float from some current value.
64 :param distance: defines neighborhood as current value +- distance
65 :param sig_dig: number of significant digits in result
66 :param lower_bound: hard cutoff lower bound that overrides distances lower than this threshold
67 :param upper_bound: hard cutoff upper bound that overrides distances higher than this threshold
68 :return: random value within neighborhood of current input value
69 """
70 upper_neighbor = current_value + distance
71 if upper_neighbor > upper_bound:
72 upper_neighbor = upper_bound
73 else:
74 upper_neighbor = upper_neighbor
75
76 lower_neighbor = current_value - distance
77 if lower_neighbor < lower_bound:
78 lower_neighbor = lower_bound
79 else:
80 lower_neighbor = lower_neighbor
81 neighbor = random.uniform(lower_neighbor, upper_neighbor)
82 neighbor = round(neighbor, sig_dig)
83 return neighbor
84
85 def _roll_dice(self, delta, t_new):
86 """
87 Method to calculate the probability of acceptance
88 :param delta: difference between the new cost and the old cost
89 :param t_new: the current temperature of the annealing algorithm
90 :return: Accept the new design? (True/False)
91 """
92 pa = math.exp(delta / t_new)
93 print('Probability of Acceptance = ' + str(pa) + ' %')
94 r = random.random()
95 print('Rolling the dice... = ' + str(r) + '%')
96 decision = pa - r
97 if decision >= 0:
98 d = True
99 else:
100 d = False
101 return d
102
103 def run_annealing(self, xy_function, x_ub, x_lb, y_ub, y_lb, x_neighbor_distance, y_neighbor_distance):
104 """ Run the annealing algorithm.
105 INPUTS:
106 xy_function :: (function) 2-D function of x and y
107 x_ub :: (float) Upper X-boundary of search space
108 x_lb :: (float) Lower X-boundary of search space
109 y_ub :: (float) Upper Y-boundary of search space
110 y_lb :: (float) Lower Y-boundary of search space
111 x_neighbor_distance :: (float) Distance between neighbors for X parameter
112 y_neighbor_distance :: (float) Distance between neighbors for Y parameter
113 RETURNS:
114 df :: (pandas dataframe) The search path of the algorithm in x, y, z (cost) coordinates
115 """
116 print('\n' + HEADER)
117 print("Testing SimulatedAnnealing Module")
118 print(HEADER + '\n')
119
120 # 1 - 2. Initialize random solution and calculate cost
121 x_list = []
122 y_list = []
123 c_list = []
124 x_old = self._random_float(x_lb, x_ub, 12)
125 y_old = self._random_float(y_lb, y_ub, 12)
126 c_old = xy_function(x_old, y_old)
127 x_list.append(x_old)
128 y_list.append(y_old)
129 c_list.append(c_old)
130 print("Initial Cost = " + str(c_old))
131
132 # 3. Generate a new neighboring solution
133 t_new = self.start_temp
134 while t_new > self.stop_temp:
135 print(SPACER)
136 print('\n' + 'Running ' + str(self.num_sims) + ' simulations @ Temperature = ' + str(t_new) + '\n')
137 print(SPACER)
138 for i in range(self.num_sims): # TODO: Not very Pythonic...
139 print('\n' + 'Simulation Run #: ' + str(i + 1))
140 x_new = self._get_neighbor_float(current_value=x_old, distance=x_neighbor_distance, lower_bound=x_lb,
141 upper_bound=x_ub,
142 sig_dig=12)
143 y_new = self._get_neighbor_float(current_value=y_old, distance=y_neighbor_distance, lower_bound=y_lb,
144 upper_bound=y_ub,
145 sig_dig=12)
146 c_new = xy_function(x_new, y_new)
147
148 # Probability of acceptance
149 delta = c_new - c_old
150 print('delta = ' + str(delta))
151 if delta >= 0:
152 print('The new solution is better - moving to it!')
153 x_old = x_new
154 y_old = y_new
155 c_old = c_new
156 else:
157 print('The new solution is worse - rolling the dice to decide whether to accept it anyway...')
158 decision = self._roll_dice(delta=delta, t_new=t_new)
159 if decision == True:
160 print('The gods have spoken: selecting the new solution even though it is worse!')
161 x_old = x_new
162 y_old = y_new
163 c_old = c_new
164 else:
165 print(
166 'The gods have smiled on us: rejecting this worse solution and sticking with the old one!')
167 x_old = x_old
168 y_old = y_old
169 c_old = c_old
170
171 print('x_old = ' + str(x_old))
172 print('y_old = ' + str(y_old))
173 print('c_old = ' + str(c_old))
174 x_list.append(x_old)
175 y_list.append(y_old)
176 c_list.append(c_old)
177 print("New Cost = " + str(c_new))
178 t_new = t_new * self.alpha
179 print(FOOTER)
180 print('New Temperature = ' + str(t_new))
181 print(FOOTER)
182
183 dfx = pd.DataFrame(x_list)
184 dfx.columns = ['x']
185 dfy = pd.DataFrame(y_list)
186 dfy.columns = ['y']
187 dfc = pd.DataFrame(c_list)
188 dfc.columns = ['c']
189 results = pd.concat([dfx, dfy, dfc], axis=1)
190 return results
191
192
193def goldstein_price_function(x, y):
194 """ Method for testing algorithm using the Goldstein-Price Function as the problem to optimize.
195 See wiki page for optimization functions: https://en.wikipedia.org/wiki/Test_functions_for_optimization
196 """
197 cost = (1 + ((x + y + 1) ** 2) * (19 - 14 * x + 3 * (x ** 2) - 14 * y + 6 * x * y + 3 * (y ** 2))) * (
198 30 + ((2 * x - 3 * y) ** 2) * (18 - 32 * x + 12 * (x ** 2) + 48 * y - 36 * x * y + 27 * (y ** 2)))
199 return cost
200
201def ackleys_function(x, y):
202 """ Method for testing algorithm using the Goldstein-Price Function as the problem to optimize.
203 See wiki page for optimization functions: https://en.wikipedia.org/wiki/Test_functions_for_optimization
204 """
205 cost = -20 * math.exp((-0.2 * (math.sqrt(0.5 * (x ** 2 + y ** 2))))) - math.exp(
206 (0.5 * ((math.cos(2 * math.pi * x)) + math.cos(2 * math.pi * y)))) + math.exp(1) + 20
207 return cost
208
209
210def plot_function_points(xy_function):
211 """ Generates a dataframe of x, y, z coordinates to """
212 x = np.linspace(-2, 2, 50)
213 y = np.linspace(-2, 2, 50)
214 x_list = []
215 y_list = []
216 cost_list = []
217 for j in range(len(x)): # TODO: Again, lots of non-pythonic stuff
218 x1 = x[j]
219 for p in range(len(y)): # TODO: Again, lots of non-pythonic stuff
220 y1 = y[p]
221 cost = xy_function(x1, y1)
222 x_list.append(x1)
223 y_list.append(y1)
224 cost_list.append(cost)
225 dfx = pd.DataFrame(x_list)
226 dfx.columns = ['x']
227 dfy = pd.DataFrame(y_list)
228 dfy.columns = ['y']
229 dfz = pd.DataFrame(cost_list)
230 dfz.columns = ['z']
231 function_data = pd.concat([dfx, dfy, dfz], axis=1)
232 return function_data
233
234
235if __name__ == "__main__":
236
237 sa = SimulatedAnnealing(start_temp=1000, stop_temp=1, alpha=0.95, num_sims=250)
238 results = sa.run_annealing(xy_function=goldstein_price_function, x_ub=2.0001, x_lb=-2.0001, y_ub=2.0001,
239 y_lb=-2.0001, x_neighbor_distance=0.1000000, y_neighbor_distance=0.1000000)
240 function_data = plot_function_points(goldstein_price_function)
241 trace1 = go.Scatter3d(
242 x=results['x'],
243 y=results['y'],
244 z=results['c'],
245 mode='markers',
246 name='Optimization Path',
247 marker=dict(
248 size=9,
249 color=results['c'],
250 colorscale='Rainbow',
251 opacity=0.95
252 )
253 )
254 trace2 = go.Scatter3d(
255 x=function_data['x'],
256 y=function_data['y'],
257 z=function_data['z'],
258 mode='markers',
259 name='Function Surface',
260 marker=dict(
261 size=6,
262 color=function_data['z'],
263 colorscale='Blues',
264 opacity=0.65
265 )
266 )
267 data = [trace2, trace1]
268 layout = go.Layout(
269 autosize=True,
270 margin=dict(l=0, r=0, b=0, t=0),
271 showlegend=True,
272 title=' Simulated Annealing Algorithm ',
273 xaxis=dict(title='X'),
274 yaxis=dict(title='Y'),
275 font=dict(size=15)
276 )
277 fig = go.Figure(data=data, layout=layout)
278 py.iplot(fig, filename='3d-scatter-colorscale')