آنیلینگ شبیه‌سازی شده¶
آنیلینگ شبیه‌سازی شده (Simulated Annealing - SA) یک الگوریتم تصادفی است که بهینه سراسری یک تابع را تقریب می‌زند. به آن الگوریتم تصادفی گفته می‌شود زیرا در جستجوی خود از مقداری تصادف استفاده می‌کند و بنابراین خروجی آن برای ورودی یکسان می‌تواند متفاوت باشد.
مسئله¶
تابعی به ما داده شده است، $E(s)$، که انرژی حالت $s$ را محاسبه می‌کند. وظیفه ما یافتن حالت $s_{best}$ است که در آن $E(s)$ کمینه می‌شود. SA برای مسائلی مناسب است که حالت‌ها گسسته بوده و $E(s)$ چندین کمینه محلی داشته باشد. ما مثال مسئله فروشنده دوره‌گرد (TSP) را بررسی خواهیم کرد.
مسئله فروشنده دوره‌گرد (TSP)¶
مجموعه‌ای از گره‌ها در فضای دو بعدی به شما داده می‌شود. هر گره با مختصات $x$ و $y$ خود مشخص می‌شود. وظیفه شما یافتن ترتیبی از گره‌ها است که مسافت طی شده هنگام بازدید از این گره‌ها به آن ترتیب را کمینه کند.
انگیزه¶
آنیلینگ (بازپخت) یک فرآیند متالورژیکی است که در آن یک ماده گرم شده و سپس به آن اجازه داده می‌شود تا خنک شود تا اتم‌های داخل آن بتوانند خود را در آرایشی با حداقل انرژی داخلی بازآرایی کنند، که به نوبه خود باعث می‌شود ماده خواص متفاوتی پیدا کند. «حالت» همان آرایش اتم‌ها و «انرژی داخلی» تابعی است که در حال کمینه شدن است. می‌توانیم حالت اولیه اتم‌ها را به عنوان یک کمینه محلی برای انرژی داخلی آن‌ها در نظر بگیریم. برای اینکه ماده اتم‌های خود را بازآرایی کند، باید آن را تشویق کنیم تا از منطقه‌ای عبور کند که انرژی داخلی آن کمینه نیست تا به کمینه سراسری برسد. این انگیزه با گرم کردن ماده تا دمای بالاتر فراهم می‌شود.
آنیلینگ شبیه‌سازی شده، به معنای واقعی کلمه، این فرآیند را شبیه‌سازی می‌کند. ما با یک حالت تصادفی (ماده) شروع می‌کنیم و دمای بالایی را تنظیم می‌کنیم (آن را گرم می‌کنیم). اکنون، الگوریتم آماده پذیرش حالت‌هایی است که انرژی بالاتری نسبت به حالت فعلی دارند، زیرا توسط دمای بالا انگیزه پیدا کرده است. این کار از گیر افتادن الگوریتم در کمینه‌های محلی جلوگیری کرده و آن را به سمت کمینه سراسری سوق می‌دهد. با گذشت زمان، الگوریتم خنک می‌شود و حالت‌های با انرژی بالاتر را رد می‌کند و به سمت نزدیک‌ترین کمینه‌ای که پیدا کرده است حرکت می‌کند.
تابع انرژی E(s)¶
$E(s)$ تابعی است که باید کمینه (یا بیشینه) شود. این تابع هر حالت را به یک عدد حقیقی نگاشت می‌کند. در مورد TSP، تابع $E(s)$ مسافت طی کردن یک دور کامل به ترتیب گره‌های موجود در حالت را برمی‌گرداند.
حالت¶
فضای حالت، دامنه تابع انرژی $E(s)$ است و یک حالت، هر عنصری است که به فضای حالت تعلق دارد. در مورد TSP، تمام مسیرهای ممکن که می‌توانیم برای بازدید از همه گره‌ها طی کنیم، فضای حالت است و هر یک از این مسیرها را می‌توان یک حالت در نظر گرفت.
حالت همسایه¶
حالت همسایه، حالتی در فضای حالت است که به حالت قبلی نزدیک است. این معمولاً به این معنی است که می‌توانیم حالت همسایه را از حالت اصلی با استفاده از یک تبدیل ساده به دست آوریم. در مورد مسئله فروشنده دوره‌گرد، یک حالت همسایه با انتخاب تصادفی ۲ گره و جابجایی موقعیت آن‌ها در حالت فعلی به دست می‌آید.
الگوریتم¶
ما با یک حالت تصادفی $s$ شروع می‌کنیم. در هر مرحله، یک حالت همسایه $s_{next}$ از حالت فعلی $s$ را انتخاب می‌کنیم. اگر $E(s_{next}) < E(s)$، آنگاه $s = s_{next}$ را به‌روزرسانی می‌کنیم. در غیر این صورت، از یک تابع پذیرش احتمالاتی $P(E(s),E(s_{next}),T)$ استفاده می‌کنیم که تصمیم می‌گیرد آیا باید به $s_{next}$ برویم یا در $s$ بمانیم. در اینجا T دما است، که در ابتدا روی مقدار بالایی تنظیم شده و با هر مرحله به آرامی کاهش می‌یابد. هرچه دما بالاتر باشد، احتمال حرکت به $s_{next}$ بیشتر است. همزمان، ما بهترین حالت $s_{best}$ را در طول تمام تکرارها نیز پیگیری می‌کنیم. این فرآیند تا رسیدن به همگرایی یا تمام شدن زمان ادامه می‌یابد.
نمایش بصری آنیلینگ شبیه‌سازی شده، در حال جستجو برای بیشینه این تابع با چندین بیشینه محلی.
دما (T) و نرخ کاهش (u)¶
دمای سیستم، میزان تمایل الگوریتم برای پذیرش حالتی با انرژی بالاتر را مشخص می‌کند. نرخ کاهش (Decay) یک ثابت است که "نرخ خنک شدن" الگوریتم را تعیین می‌کند. مشخص شده است که نرخ خنک شدن آهسته (u بزرگتر) نتایج بهتری می‌دهد.
تابع پذیرش احتمالاتی (PAF)¶
$P(E,E_{next},T) = \begin{cases} \text{True} &\quad\text{if } \mathcal{U}_{[0,1]} \le \exp(-\frac{E_{next}-E}{T}) \\ \text{False} &\quad\text{otherwise}\\ \end{cases}$
در اینجا، $\mathcal{U}_{[0,1]}$ یک مقدار تصادفی پیوسته و یکنواخت در بازه $[0,1]$ است. این تابع، حالت فعلی، حالت بعدی و دما را دریافت کرده و یک مقدار بولین برمی‌گرداند که به جستجوی ما می‌گوید آیا باید به $s_{next}$ برود یا در $s$ بماند. توجه داشته باشید که برای $E_{next} < E$، این تابع همیشه True برمی‌گرداند، در غیر این صورت همچنان می‌تواند با احتمال $\exp(-\frac{E_{next}-E}{T})$ حرکت کند که با اندازه گیبس (Gibbs measure) مطابقت دارد.
bool P(double E,double E_next,double T,mt19937 rng){
double prob = exp(-(E_next-E)/T);
if(prob > 1) return true;
else{
bernoulli_distribution d(prob);
return d(rng);
}
}
قالب کد¶
class state {
public:
state() {
// حالت اولیه را تولید کنید
}
state next() {
state s_next;
// s_next را به یک حالت همسایه تصادفی تغییر دهید
return s_next;
}
double E() {
// تابع انرژی را اینجا پیاده‌سازی کنید
};
};
pair<double, state> simAnneal() {
state s = state();
state best = s;
double T = 10000; // دمای اولیه
double u = 0.995; // نرخ کاهش
double E = s.E();
double E_next;
double E_best = E;
mt19937 rng(chrono::steady_clock::now().time_since_epoch().count());
while (T > 1) {
state next = s.next();
E_next = next.E();
if (P(E, E_next, T, rng)) {
s = next;
if (E_next < E_best) {
best = s;
E_best = E_next;
}
E = E_next;
}
T *= u;
}
return {E_best, best};
}
نحوه استفاده:¶
توابع کلاس state را به شکل مناسب پر کنید. اگر به دنبال یافتن یک بیشینه سراسری و نه کمینه هستید، اطمینان حاصل کنید که تابع $E()$ منفی تابعی را که می‌خواهید بیشینه کنید برمی‌گرداند و در انتها $-E_{best}$ را چاپ کنید. پارامترهای زیر را بر اساس نیاز خود تنظیم کنید.
پارامترها¶
- $T$: دمای اولیه. اگر می‌خواهید جستجو برای مدت طولانی‌تری اجرا شود، آن را روی مقدار بالاتری تنظیم کنید.
- $u$: نرخ کاهش (Decay). نرخ خنک شدن را تعیین می‌کند. نرخ خنک شدن آهسته‌تر (مقدار بزرگتر u) معمولاً نتایج بهتری می‌دهد، اما به قیمت اجرای طولانی‌تر. اطمینان حاصل کنید که $u < 1$.
تعداد تکرارهای حلقه از عبارت زیر به دست می‌آید:
$N = \lceil -\log_{u}{T} \rceil$
نکاتی برای انتخاب $T$ و $u$: اگر کمینه‌های محلی زیادی و فضای حالت وسیعی وجود دارد، $u = 0.999$ را برای نرخ خنک شدن آهسته تنظیم کنید، که به الگوریتم اجازه می‌دهد تا امکانات بیشتری را کاوش کند. از طرف دیگر، اگر فضای حالت محدودتر است، $u = 0.99$ کافی خواهد بود. اگر مطمئن نیستید، با تنظیم $u = 0.998$ یا بالاتر، جانب احتیاط را رعایت کنید. پیچیدگی زمانی یک تکرار الگوریتم را محاسبه کنید و از آن برای تخمین مقداری برای $N$ استفاده کنید که از TLE جلوگیری کند، سپس از فرمول زیر برای به دست آوردن $T$ استفاده کنید.
$T = u^{-N}$
پیاده‌سازی نمونه برای TSP¶
class state {
public:
vector<pair<int, int>> points;
std::mt19937 mt{ static_cast<std::mt19937::result_type>(
std::chrono::steady_clock::now().time_since_epoch().count()
) };
state() {
points = {{0,0},{2,2},{0,2},{2,0},{0,1},{1,2},{2,1},{1,0}} ;
}
state next() {
state s_next;
s_next.points = points;
uniform_int_distribution<> choose(0, points.size()-1);
int a = choose(mt);
int b = choose(mt);
s_next.points[a].swap(s_next.points[b]);
return s_next;
}
double euclidean(pair<int, int> a, pair<int, int> b) {
return hypot(a.first - b.first, a.second - b.second);
}
double E() {
double dist = 0;
int n = points.size();
for (int i = 0;i < n; i++)
dist += euclidean(points[i], points[(i+1)%n]);
return dist;
};
};
int main() {
pair<double, state> res;
res = simAnneal();
double E_best = res.first;
state best = res.second;
cout << "Lenght of shortest path found : " << E_best << "\n";
cout << "Order of points in shortest path : \n";
for(auto x: best.points) {
cout << x.first << " " << x.second << "\n";
}
}
اصلاحات بیشتر در الگوریتم:¶
- یک شرط خروج مبتنی بر زمان به حلقه while اضافه کنید تا از TLE جلوگیری شود.
- کاهش دمای پیاده‌سازی شده در بالا، یک کاهش نمایی است. شما همیشه می‌توانید این را با یک تابع کاهش دیگر متناسب با نیاز خود جایگزین کنید.
- تابع پذیرش احتمالاتی ارائه شده در بالا، به دلیل وجود عامل $E_{next} - E$ در صورت کسر توان، پذیرش حالت‌هایی با انرژی کمتر را ترجیح می‌دهد. شما می‌توانید به سادگی این عامل را حذف کنید تا PAF مستقل از تفاوت انرژی‌ها شود.
- تأثیر تفاوت انرژی‌ها، $E_{next} - E$، بر روی PAF را می‌توان با افزایش/کاهش پایه توان، همانطور که در زیر نشان داده شده، افزایش/کاهش داد:
bool P(double E, double E_next, double T, mt19937 rng) { double e = 2; // e را هر عدد حقیقی بزرگتر از 1 قرار دهید double prob = pow(e,-(E_next-E)/T); if (prob > 1) return true; else { bernoulli_distribution d(prob); return d(rng); } }