شبیه سازی تبرید (Simulated Annealing) – به زبان ساده

۴۱۵۸ بازدید
آخرین به‌روزرسانی: ۱۲ مهر ۱۴۰۲
زمان مطالعه: ۱۴ دقیقه
شبیه سازی تبرید (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')
بر اساس رای ۱۷ نفر
آیا این مطلب برای شما مفید بود؟
اگر بازخوردی درباره این مطلب دارید یا پرسشی دارید که بدون پاسخ مانده است، آن را از طریق بخش نظرات مطرح کنید.
نظر شما چیست؟

نشانی ایمیل شما منتشر نخواهد شد. بخش‌های موردنیاز علامت‌گذاری شده‌اند *