সিমুলেটেড অ্যানিলিং#

সিমুলেটেড অ্যানিলিং (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$ এর একটি প্রতিবেশী স্টেট $s_{next}$ বেছে নিই। যদি $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)#

সিস্টেমের তাপমাত্রা উচ্চ এনার্জি বিশিষ্ট স্টেট গ্রহণে অ্যালগরিদমের ইচ্ছুকতা পরিমাপ করে। ক্ষয় হলো একটি ধ্রুবক যা অ্যালগরিদমের “কুলিং রেট” পরিমাপ করে। ধীর কুলিং রেট (বৃহত্তর $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})$ সম্ভাবনায় মুভ করতে পারে, যা গিবস মেজার এর সাথে সঙ্গতিপূর্ণ।

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() {
        // Generate the initial state
    }
    state next() {
        state s_next;
        // Modify s_next to a random neighboring state
        return s_next;
    }
    double E() {
        // Implement the energy function here
    };
};


pair<double, state> simAnneal() {
    state s = state();
    state best = s;
    double T = 10000; // Initial temperature
    double u = 0.995; // decay rate
    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$ : ক্ষয়। কুলিংয়ের হার নির্ধারণ করে। ধীর কুলিং রেট (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 = {%raw%} {{0,0},{2,2},{0,2},{2,0},{0,1},{1,2},{2,1},{1,0}} {%endraw%};
    }
    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";
    }
}

অ্যালগরিদমের আরও পরিবর্তন:#

• TLE প্রতিরোধে while লুপে একটি সময়-ভিত্তিক প্রস্থান শর্ত যোগ করুন
• উপরে ইমপ্লিমেন্ট করা ক্ষয়টি একটি এক্সপোনেনশিয়াল ক্ষয়। আপনি এটিকে আপনার প্রয়োজন অনুযায়ী একটি ক্ষয় ফাংশন দ্বারা প্রতিস্থাপন করতে পারেন।
• উপরে দেওয়া প্রবাবিলিটি অ্যাক্সেপ্টেন্স ফাংশন, এক্সপোনেন্টের গণনায় $E_{next} - E$ ফ্যাক্টরের কারণে কম এনার্জির স্টেট গ্রহণ করতে পছন্দ করে। আপনি এই ফ্যাক্টরটি সরিয়ে PAF কে এনার্জির পার্থক্য থেকে স্বাধীন করতে পারেন।
• এনার্জির পার্থক্য $E_{next} - E$ এর PAF এর উপর প্রভাব নিচে দেখানো অনুযায়ী এক্সপোনেন্টের বেস বাড়িয়ে/কমিয়ে বাড়ানো/কমানো যায়:

bool P(double E, double E_next, double T, mt19937 rng) {
    double e = 2; // set e to any real number greater than 1
    double prob =  pow(e,-(E_next-E)/T);
    if (prob > 1)
        return true;
    else {
        bernoulli_distribution d(prob);
        return d(rng);
    }
}

সমস্যা#

• USACO Jan 2017 - Subsequence Reversal
• Deltix Summer 2021 - DIY Tree
• AtCoder Contest Scheduling