跳转至

模拟退火算法

概述

模拟退火（Simulated Annealing, SA）是一种随机搜索算法，模拟金属退火过程，通过温度下降控制搜索过程，逐渐收敛到全局最优解。

物理背景

退火是将金属加热到高温，然后缓慢冷却的过程。高温时原子剧烈运动，可能跳出局部势阱；冷却后原子趋于稳定，达到最低能量状态。

模拟退火特点

  • 随机搜索：具有随机性，可能跳出局部最优
  • 温度控制：高温时接受劣解概率高，低温时趋于贪婪
  • 全局优化：理论上可收敛到全局最优
  • 通用性强：适用于各种组合优化问题

核心思想

模拟退火通过Metropolis准则接受新解：

  • 若新解更优，一定接受
  • 若新解较差，以概率 exp(-ΔE/T) 接受

其中ΔE为新旧解的差值，T为当前温度。

算法框架

C
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
#include <stdlib.h>
#include <math.h>
#include <time.h>

#define INITIAL_TEMP 10000.0
#define FINAL_TEMP 0.001
#define ALPHA 0.99
#define ITERATIONS 1000

double randomDouble() {
    return (double)rand() / RAND_MAX;
}

void simulatedAnnealing(
    void *state,
    double (*energy)(void*),
    void (*neighbor)(void*, void*),
    double initialTemp,
    double finalTemp,
    double alpha,
    int iterationsPerTemp
) {
    double currentTemp = initialTemp;
    double currentEnergy = energy(state);
    void *newState = malloc(sizeof(state));
    
    while (currentTemp > finalTemp) {
        for (int i = 0; i < iterationsPerTemp; i++) {
            neighbor(state, newState);
            double newEnergy = energy(newState);
            double delta = newEnergy - currentEnergy;
            
            if (delta < 0 || randomDouble() < exp(-delta / currentTemp)) {
                memcpy(state, newState, sizeof(state));
                currentEnergy = newEnergy;
            }
        }
        
        currentTemp *= alpha;
    }
    
    free(newState);
}

旅行商问题（TSP）

C
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
#define MAX_CITIES 100

typedef struct {
    int cities[MAX_CITIES];
    int n;
    double dist[MAX_CITIES][MAX_CITIES];
} TSPState;

double tspEnergy(TSPState *state) {
    double total = 0;
    for (int i = 0; i < state->n; i++) {
        int from = state->cities[i];
        int to = state->cities[(i + 1) % state->n];
        total += state->dist[from][to];
    }
    return total;
}

void tspNeighbor(TSPState *current, TSPState *neighbor) {
    *neighbor = *current;
    
    int i = rand() % neighbor->n;
    int j = rand() % neighbor->n;
    
    while (i == j) {
        j = rand() % neighbor->n;
    }
    
    int temp = neighbor->cities[i];
    neighbor->cities[i] = neighbor->cities[j];
    neighbor->cities[j] = temp;
}

void tsp2OptNeighbor(TSPState *current, TSPState *neighbor) {
    *neighbor = *current;
    
    int i = rand() % neighbor->n;
    int j = rand() % neighbor->n;
    
    if (i > j) {
        int temp = i;
        i = j;
        j = temp;
    }
    
    while (i < j) {
        int temp = neighbor->cities[i];
        neighbor->cities[i] = neighbor->cities[j];
        neighbor->cities[j] = temp;
        i++;
        j--;
    }
}

double solveTSP_SA(TSPState *state) {
    double currentTemp = INITIAL_TEMP;
    double currentEnergy = tspEnergy(state);
    TSPState neighbor;
    
    double bestEnergy = currentEnergy;
    int bestCities[MAX_CITIES];
    for (int i = 0; i < state->n; i++) {
        bestCities[i] = state->cities[i];
    }
    
    while (currentTemp > FINAL_TEMP) {
        for (int i = 0; i < ITERATIONS; i++) {
            tsp2OptNeighbor(state, &neighbor);
            double newEnergy = tspEnergy(&neighbor);
            double delta = newEnergy - currentEnergy;
            
            if (delta < 0 || randomDouble() < exp(-delta / currentTemp)) {
                *state = neighbor;
                currentEnergy = newEnergy;
                
                if (currentEnergy < bestEnergy) {
                    bestEnergy = currentEnergy;
                    for (int j = 0; j < state->n; j++) {
                        bestCities[j] = state->cities[j];
                    }
                }
            }
        }
        
        currentTemp *= ALPHA;
    }
    
    for (int i = 0; i < state->n; i++) {
        state->cities[i] = bestCities[i];
    }
    
    return bestEnergy;
}

函数优化问题

C
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
typedef struct {
    double x;
    double y;
} Point;

double rastrigin(Point *p) {
    double A = 10.0;
    double value = 2 * A;
    value += p->x * p->x - A * cos(2 * M_PI * p->x);
    value += p->y * p->y - A * cos(2 * M_PI * p->y);
    return value;
}

void pointNeighbor(Point *current, Point *neighbor, double stepSize) {
    double angle = 2 * M_PI * randomDouble();
    double r = stepSize * randomDouble();
    
    neighbor->x = current->x + r * cos(angle);
    neighbor->y = current->y + r * sin(angle);
}

Point optimizeFunction_SA() {
    Point current = {rand() % 10 - 5, rand() % 10 - 5};
    double currentEnergy = rastrigin(&current);
    
    double currentTemp = INITIAL_TEMP;
    Point neighbor;
    
    while (currentTemp > FINAL_TEMP) {
        for (int i = 0; i < ITERATIONS; i++) {
            pointNeighbor(&current, &neighbor, currentTemp / 100);
            double newEnergy = rastrigin(&neighbor);
            double delta = newEnergy - currentEnergy;
            
            if (delta < 0 || randomDouble() < exp(-delta / currentTemp)) {
                current = neighbor;
                currentEnergy = newEnergy;
            }
        }
        
        currentTemp *= ALPHA;
    }
    
    return current;
}

0-1背包问题

C
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
#define MAX_ITEMS 100

typedef struct {
    int selected[MAX_ITEMS];
    int n;
    int weight[MAX_ITEMS];
    int value[MAX_ITEMS];
    int capacity;
} KnapsackState;

double knapsackEnergy(KnapsackState *state) {
    int totalWeight = 0;
    int totalValue = 0;
    
    for (int i = 0; i < state->n; i++) {
        if (state->selected[i]) {
            totalWeight += state->weight[i];
            totalValue += state->value[i];
        }
    }
    
    if (totalWeight > state->capacity) {
        return -totalWeight * 1000;
    }
    
    return -totalValue;
}

void knapsackNeighbor(KnapsackState *current, KnapsackState *neighbor) {
    *neighbor = *current;
    
    int idx = rand() % neighbor->n;
    neighbor->selected[idx] = !neighbor->selected[idx];
}

C++ 实现

C++
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
template<typename State>
class SimulatedAnnealing {
private:
    function<double(const State&)> energy;
    function<State(const State&, double)> neighbor;
    
    double initialTemp;
    double finalTemp;
    double coolingRate;
    int iterationsPerTemp;
    
public:
    SimulatedAnnealing(
        function<double(const State&)> e,
        function<State(const State&, double)> n,
        double initT = 10000.0,
        double finalT = 0.001,
        double alpha = 0.99,
        int iter = 1000
    ) : energy(e), neighbor(n), initialTemp(initT), 
        finalTemp(finalT), coolingRate(alpha), 
        iterationsPerTemp(iter) {}
    
    State solve(State initial, int maxIterations = 100000) {
        State current = initial;
        double currentEnergy = energy(current);
        
        State best = current;
        double bestEnergy = currentEnergy;
        
        double temp = initialTemp;
        int iteration = 0;
        
        while (temp > finalTemp && iteration < maxIterations) {
            for (int i = 0; i < iterationsPerTemp; i++) {
                State next = neighbor(current, temp);
                double nextEnergy = energy(next);
                double delta = nextEnergy - currentEnergy;
                
                if (delta < 0 || 
                    (double)rand() / RAND_MAX < exp(-delta / temp)) {
                    current = next;
                    currentEnergy = nextEnergy;
                    
                    if (currentEnergy < bestEnergy) {
                        best = current;
                        bestEnergy = currentEnergy;
                    }
                }
                
                iteration++;
            }
            
            temp *= coolingRate;
        }
        
        return best;
    }
};

class TSP_SA {
public:
    vector<int> solve(const vector<vector<double>>& dist) {
        int n = dist.size();
        
        vector<int> current(n);
        iota(current.begin(), current.end(), 0);
        random_shuffle(current.begin(), current.end());
        
        auto energy = [&](const vector<int>& path) {
            double total = 0;
            for (int i = 0; i < n; i++) {
                total += dist[path[i]][path[(i + 1) % n]];
            }
            return total;
        };
        
        auto neighbor = [&](const vector<int>& path, double temp) {
            vector<int> newPath = path;
            int i = rand() % n;
            int j = rand() % n;
            if (i > j) swap(i, j);
            reverse(newPath.begin() + i, newPath.begin() + j + 1);
            return newPath;
        };
        
        SimulatedAnnealing<vector<int>> sa(energy, neighbor);
        return sa.solve(current);
    }
};

参数选择

冷却方案

方案 公式 特点
指数冷却 T(k+1) = α × T(k) 最常用，简单有效
线性冷却 T(k+1) = T(k) - β 收敛较快
对数冷却 T(k) = T₀ / ln(k+1) 理论最优，收敛慢

参数建议

C
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
typedef struct {
    double initialTemp;
    double finalTemp;
    double coolingRate;
    int iterationsPerTemp;
} SAParams;

SAParams getRecommendedParams(int problemSize) {
    SAParams params;
    
    params.initialTemp = 10000.0;
    params.finalTemp = 0.001;
    params.coolingRate = 0.99;
    params.iterationsPerTemp = problemSize * 100;
    
    return params;
}

时间复杂度

  • 时间复杂度：O(iterations × iterationsPerTemp)
  • 实际运行时间取决于参数设置
  • 通常远少于精确算法

收敛性

理论保证：若温度下降足够慢（T_k = T₀/ln(k)），模拟退火以概率1收敛到全局最优。

实际应用：使用更快的冷却方案，虽然不保证全局最优，但通常能得到较好的解。

应用场景

  1. 组合优化：TSP、背包、调度
  2. 函数优化：多峰函数求极值
  3. 机器学习：神经网络训练
  4. 电路设计：VLSI布局布线
  5. 图像处理：图像重建、分割

优缺点

优点

  1. 全局搜索能力：可跳出局部最优
  2. 易于实现：框架简单通用
  3. 鲁棒性强：对初始解不敏感
  4. 可并行化：可多次运行取最优

缺点

  1. 参数敏感：需要调参
  2. 收敛速度慢：需要足够迭代
  3. 无最优性保证：实际应用中可能陷入局部最优
  4. 效率依赖问题：邻域结构影响搜索效率

改进策略

  1. 自适应退火：根据接受率调整温度
  2. 多次运行：从不同初始解出发
  3. 混合策略：结合局部搜索
  4. 并行模拟退火：多个进程同时搜索

参考资料