Simulated Annealing

1
#include <bits/stdc++.h>
2
using namespace std;
3
4
// Rastrigin function in 2D: global minimum at (0,0) with value 0
5
static inline double rastrigin(double x, double y) {
6
    const double A = 10.0;
7
    return 2.0 * A + (x * x - A * cos(2.0 * M_PI * x)) + (y * y - A * cos(2.0 * M_PI * y));
8
}
9
10
int main() {
11
    ios::sync_with_stdio(false);
12
    cin.tie(nullptr);
13
14
    mt19937_64 rng(chrono::steady_clock::now().time_since_epoch().count());
15
    normal_distribution<double> gauss(0.0, 1.0); // standard normal for perturbations
16
    uniform_real_distribution<double> uni01(0.0, 1.0);
17
18
    // Initial solution near the origin
19
    double x = 5.0, y = 5.0;
20
    double fx = rastrigin(x, y);
21
22
    // Best found
23
    double bx = x, by = y, bfx = fx;
24
25
    // Annealing parameters
26
    double T = 1e2;           // initial temperature
27
    const double T_end = 1e-6; // stopping temperature
28
    const double beta = 0.9995; // exponential cooling multiplier
29
    const int iters_per_T = 500; // proposals per temperature level
30
31
    // Typical step size; scaled by temperature for exploration
32
    const double base_step = 1.0;
33
34
    while (T > T_end) {
35
        for (int it = 0; it < iters_per_T; ++it) {
36
            // Propose a neighbor using Gaussian perturbation scaled by sqrt(T)
37
            double nx = x + gauss(rng) * base_step * sqrt(T);
38
            double ny = y + gauss(rng) * base_step * sqrt(T);
39
40
            // Optional: clamp to search bounds if any (here we allow any real value)
41
            double fn = rastrigin(nx, ny);
42
            double delta = fn - fx; // we minimize
43
44
            if (delta <= 0.0) {
45
                // Always accept improvements
46
                x = nx; y = ny; fx = fn;
47
            } else {
48
                // Accept uphill move with probability exp(-delta / T)
49
                double p = exp(-delta / T);
50
                if (uni01(rng) < p) {
51
                    x = nx; y = ny; fx = fn;
52
                }
53
            }
54
55
            // Track best
56
            if (fx < bfx) {
57
                bfx = fx; bx = x; by = y;
58
            }
59
        }
60
        // Cool down
61
        T *= beta;
62
    }
63
64
    // Output best solution found: x y value
65
    cout << fixed << setprecision(6) << bx << ' ' << by << ' ' << bfx << '\n';
66
    return 0;
67
}
68

1
#include <bits/stdc++.h>
2
using namespace std;
3
4
struct Pt { double x, y; };
5
6
static inline double dist(const Pt& a, const Pt& b) {
7
    double dx = a.x - b.x, dy = a.y - b.y;
8
    return sqrt(dx*dx + dy*dy);
9
}
10
11
int main() {
12
    ios::sync_with_stdio(false);
13
    cin.tie(nullptr);
14
15
    mt19937_64 rng(chrono::steady_clock::now().time_since_epoch().count());
16
    uniform_real_distribution<double> uni01(0.0, 1.0);
17
    uniform_int_distribution<int> uN;
18
19
    // Generate synthetic instance
20
    int n = 300; // number of cities
21
    vector<Pt> p(n);
22
    for (int i = 0; i < n; ++i) {
23
        p[i].x = (double) (rng() % 100000) / 100.0; // 0..1000
24
        p[i].y = (double) (rng() % 100000) / 100.0;
25
    }
26
27
    // Initial tour: 0..n-1
28
    vector<int> tour(n);
29
    iota(tour.begin(), tour.end(), 0);
30
31
    auto tour_len = [&](const vector<int>& t){
32
        double L = 0.0;
33
        for (int i = 0; i < n; ++i) {
34
            int a = t[i], b = t[(i+1)%n];
35
            L += dist(p[a], p[b]);
36
        }
37
        return L;
38
    };
39
40
    double cur = tour_len(tour);
41
    double best = cur;
42
    vector<int> best_tour = tour;
43
44
    // Annealing parameters
45
    double T = 1000.0;
46
    const double T_end = 1e-3;
47
    const double beta = 0.9995; // exponential cooling
48
    const int iters_per_T = 2000;
49
    const int KMAX = 50; // cap segment length to keep reversal cost bounded
50
51
    while (T > T_end) {
52
        for (int iter = 0; iter < iters_per_T; ++iter) {
53
            // Choose i < j with limited segment length
54
            int i = (int)(rng() % n);
55
            int j = i + 1 + (int)(rng() % min(KMAX, n-1));
56
            if (j >= n) j = j - n; // wrap around
57
58
            // Ensure i < j in circular sense by mapping to linear indices
59
            int a = i, b = j;
60
            if (a > b) swap(a, b);
61
62
            // Identify endpoints: (a-1,a) and (b,b+1) will be replaced by (a-1,b) and (a,b+1)
63
            int a_prev = (a - 1 + n) % n;
64
            int b_next = (b + 1) % n;
65
            int A = tour[a_prev], B = tour[a], C = tour[b], D = tour[b_next];
66
67
            double before = dist(p[A], p[B]) + dist(p[C], p[D]);
68
            double after  = dist(p[A], p[C]) + dist(p[B], p[D]);
69
            double delta = after - before; // change in tour length if we reverse [a..b]
70
71
            bool accept = false;
72
            if (delta <= 0.0) accept = true;
73
            else {
74
                double prob = exp(-delta / T);
75
                if (uni01(rng) < prob) accept = true;
76
            }
77
78
            if (accept) {
79
                // Reverse segment [a..b]
80
                reverse(tour.begin() + a, tour.begin() + b + 1);
81
                cur += delta;
82
                if (cur < best) {
83
                    best = cur;
84
                    best_tour = tour;
85
                }
86
            }
87
        }
88
        T *= beta;
89
    }
90
91
    // Output best length and first few cities
92
    cout << fixed << setprecision(6) << best << '\n';
93
    for (int i = 0; i < min(10, n); ++i) cout << best_tour[i] << (i+1<min(10,n)? ' ' : '\n');
94
    return 0;
95
}
96