每日算法 — 使用java实现数字华容道:模拟退火全局优化与自适应降温策略

引言

数字华容道(又称15-puzzle或滑动拼图)是一款经典的益智游戏:在一个 n × n 的棋盘上,放置 n² - 1 个带数字的方块和一个空白格,玩家通过上下左右滑动相邻方块,将棋盘从初始乱序状态复原为目标顺序状态。当棋盘规模增大时,状态空间呈阶乘级爆炸,传统的BFS或A搜索往往因内存或时间限制而失效。模拟退火(Simulated Annealing, SA)*作为一种基于物理退火过程的全局随机优化算法,能够在巨大状态空间中跳出局部最优,以较高概率逼近全局最优解。本文用Java实现数字华容道的模拟退火求解器,详解状态表示、邻域搜索、能量函数设计与自适应降温 schedule。

问题建模

状态表示

用一维数组 board[16] 表示 4 × 4 棋盘的当前布局,其中 0 代表空白格。

目标状态:

1  2  3  4
5  6  7  8
9  10 11 12
13 14 15 0

可解性判定

对于 4 × 4 数字华容道,设逆序对数为 inv,空白格所在行(从底部数)为 row,当且仅当 inv + row 为偶数时,该局面可解。

能量函数(目标函数)

模拟退火需要定义”能量”来衡量当前状态与目标的距离。本文采用加权曼哈顿距离:

E(state) = Σ |current_row(i) - target_row(i)| + |current_col(i) - target_col(i)|

能量越低,状态越接近目标。当 E = 0 时, puzzle 已复原。

核心算法:模拟退火

算法思路

  1. 初始化:从当前乱序状态出发,设定初始温度 T₀
  2. 邻域搜索:每次随机选择一个可移动方向,生成邻域状态
  3. Metropolis准则:
  4. 若邻域能量 E' < E,一定接受该移动
  5. 若 E' ≥ E,以概率 exp(-(E' - E) / T) 接受该移动(允许”上山”跳出局部最优)
  6. 降温:按 schedule 降低温度 T = α × T,α ∈ (0.999, 0.9999)
  7. 终止:温度低于阈值或达到最大迭代次数,输出最优解

自适应降温策略

固定降温系数往往导致收敛过早或过慢。本文采用自适应调整:
– 若连续多次接受率过高(>80%),说明温度仍太高,加快降温
– 若连续多次接受率过低(<5%),说明陷入局部最优,适当升温”回火”

复杂度分析

  • 单次迭代:邻域生成 O(1),能量计算 O(n²),接受判定 O(1)
  • 总复杂度:O(K × n²),K 为迭代次数,与问题难度相关
  • 空间复杂度:O(n²),仅需存储当前状态

Java完整实现

import java.util.*;

/**
 * 数字华容道:模拟退火全局优化与自适应降温策略
 * 
 * 核心思路:
 * 1. 将4×4棋盘用一维数组表示,0代表空白格
 * 2. 能量函数采用加权曼哈顿距离衡量与目标状态的差距
 * 3. 每次随机选择合法移动方向生成邻域状态
 * 4. 使用Metropolis准则决定是否接受劣解,允许跳出局部最优
 * 5. 自适应降温:根据接受率动态调整降温系数,必要时回火升温
 */
public class SlidingPuzzleSA {

    static final int N = 4;           // 棋盘边长
    static final int SIZE = N * N;    // 总格子数

    // 目标状态:1~15按顺序,0在最后
    static final int[] GOAL = {
        1, 2, 3, 4,
        5, 6, 7, 8,
        9, 10, 11, 12,
        13, 14, 15, 0
    };

    // 四个移动方向:上、下、左、右(对应一维数组的位移量)
    static final int[] DX = {-1, 1, 0, 0};
    static final int[] DY = {0, 0, -1, 1};
    static final String[] DIR_NAME = {"UP", "DOWN", "LEFT", "RIGHT"};

    private Random rand;
    private int[] current;    // 当前状态
    private int blankPos;     // 空白格位置
    private int bestEnergy;   // 历史最优能量
    private int[] bestState;  // 历史最优状态

    public SlidingPuzzleSA() {
        this.rand = new Random();
    }

    // ============ 核心工具方法 ============

    /** 将二维坐标转一维索引 */
    private int idx(int r, int c) {
        return r * N + c;
    }

    /** 获取一维索引对应的行 */
    private int row(int idx) {
        return idx / N;
    }

    /** 获取一维索引对应的列 */
    private int col(int idx) {
        return idx % N;
    }

    /** 查找空白格位置 */
    private int findBlank(int[] state) {
        for (int i = 0; i < SIZE; i++) {
            if (state[i] == 0) return i;
        }
        return -1;
    }

    /** 计算当前状态的曼哈顿距离能量 */
    public int energy(int[] state) {
        int e = 0;
        for (int i = 0; i < SIZE; i++) {
            int val = state[i];
            if (val == 0) continue;
            int targetPos = val - 1; // 数字val的目标位置
            e += Math.abs(row(i) - row(targetPos)) + Math.abs(col(i) - col(targetPos));
        }
        return e;
    }

    /** 计算逆序对数,用于判断可解性 */
    public int inversionCount(int[] state) {
        int inv = 0;
        for (int i = 0; i < SIZE; i++) {
            if (state[i] == 0) continue;
            for (int j = i + 1; j < SIZE; j++) {
                if (state[j] != 0 && state[i] > state[j]) {
                    inv++;
                }
            }
        }
        return inv;
    }

    /** 判断当前状态是否可解 */
    public boolean isSolvable(int[] state) {
        int inv = inversionCount(state);
        int blankRowFromBottom = N - row(findBlank(state));
        // 4×4规则:inv + blankRowFromBottom 必须为偶数
        return (inv + blankRowFromBottom) % 2 == 0;
    }

    /** 检查是否达到目标状态 */
    public boolean isGoal(int[] state) {
        return Arrays.equals(state, GOAL);
    }

    /** 深拷贝状态 */
    private int[] copyState(int[] state) {
        return Arrays.copyOf(state, SIZE);
    }

    /** 打印棋盘 */
    public void printBoard(int[] state) {
        for (int r = 0; r < N; r++) {
            for (int c = 0; c < N; c++) {
                int v = state[idx(r, c)];
                if (v == 0) System.out.print("  . ");
                else System.out.printf("%3d ", v);
            }
            System.out.println();
        }
    }

    // ============ 模拟退火核心 ============

    /**
     * 对给定初始状态执行模拟退火求解
     * @param initialState 初始乱序状态
     * @param maxIter 最大迭代次数
     * @param initialTemp 初始温度
     * @return 求解结果(是否成功、步数、最终状态)
     */
    public SAResult solve(int[] initialState, int maxIter, double initialTemp) {
        if (!isSolvable(initialState)) {
            return new SAResult(false, 0, null, "局面不可解");
        }

        current = copyState(initialState);
        blankPos = findBlank(current);
        int currentEnergy = energy(current);
        bestEnergy = currentEnergy;
        bestState = copyState(current);

        double T = initialTemp;
        int accepted = 0;
        int totalSteps = 0;
        List<String> moveHistory = new ArrayList<>();

        // 自适应参数
        double alphaBase = 0.9995;      // 基础降温系数
        double alpha = alphaBase;
        int windowSize = 1000;          // 接受率统计窗口
        int[] acceptWindow = new int[windowSize];
        int windowIdx = 0;

        for (int iter = 0; iter < maxIter; iter++) {
            // 1. 生成邻域:随机选择一个合法移动方向
            List<Integer> validDirs = new ArrayList<>();
            int br = row(blankPos), bc = col(blankPos);
            for (int d = 0; d < 4; d++) {
                int nr = br + DX[d], nc = bc + DY[d];
                if (nr >= 0 && nr < N && nc >= 0 && nc < N) {
                    validDirs.add(d);
                }
            }
            if (validDirs.isEmpty()) continue;

            int dir = validDirs.get(rand.nextInt(validDirs.size()));
            int nr = br + DX[dir], nc = bc + DY[d];
            int newBlankPos = idx(nr, nc);

            // 2. 快速计算能量变化(无需重算全部)
            int movedTile = current[newBlankPos];
            int oldTargetPos = movedTile - 1;
            int oldDist = Math.abs(nr - row(oldTargetPos)) + Math.abs(nc - col(oldTargetPos));
            int newDist = Math.abs(br - row(oldTargetPos)) + Math.abs(bc - col(oldTargetPos));
            int deltaE = newDist - oldDist; // 能量变化量

            // 3. Metropolis准则
            boolean acceptedMove;
            if (deltaE < 0) {
                acceptedMove = true; // 能量降低,一定接受
            } else {
                double prob = Math.exp(-deltaE / T);
                acceptedMove = rand.nextDouble() < prob;
            }

            if (acceptedMove) {
                // 执行移动:交换空白格与目标格子
                current[blankPos] = movedTile;
                current[newBlankPos] = 0;
                blankPos = newBlankPos;
                currentEnergy += deltaE;
                moveHistory.add(DIR_NAME[dir]);
                accepted++;

                // 更新历史最优
                if (currentEnergy < bestEnergy) {
                    bestEnergy = currentEnergy;
                    bestState = copyState(current);
                    if (bestEnergy == 0) {
                        return new SAResult(true, moveHistory.size(), copyState(current),
                            String.format("模拟退火成功!迭代%d次,实际移动%d步", iter + 1, moveHistory.size()));
                    }
                }
            }

            acceptWindow[windowIdx % windowSize] = acceptedMove ? 1 : 0;
            windowIdx++;
            totalSteps++;

            // 4. 自适应降温
            if (iter > 0 && iter % windowSize == 0) {
                int recentAccept = 0;
                for (int v : acceptWindow) recentAccept += v;
                double acceptRate = (double) recentAccept / windowSize;

                if (acceptRate > 0.80) {
                    alpha = Math.max(0.9990, alpha - 0.0002); // 接受率过高,加快降温
                } else if (acceptRate < 0.05) {
                    alpha = Math.min(0.9999, alpha + 0.0003); // 陷入局部最优,回火升温
                    T *= 1.5; // 直接升温
                } else {
                    alpha = alphaBase; // 恢复正常
                }
            }

            T *= alpha;
            if (T < 1e-8) T = 1e-8; // 温度下限
        }

        // 返回历史最优结果
        boolean success = (bestEnergy == 0);
        return new SAResult(success, moveHistory.size(), bestState,
            String.format("迭代结束。最优能量=%d, 历史最优%s", bestEnergy,
                success ? "已达目标" : "未达目标"));
    }

    /** 生成一个可解的随机初始状态 */
    public int[] generateRandomState(int shuffleSteps) {
        int[] state = copyState(GOAL);
        int bp = SIZE - 1;
        Random r = new Random();

        // 从目标状态开始随机游走shuffleSteps步,保证可解
        for (int i = 0; i < shuffleSteps; i++) {
            List<Integer> valid = new ArrayList<>();
            int br = row(bp), bc = col(bp);
            for (int d = 0; d < 4; d++) {
                int nr = br + DX[d], nc = bc + DY[d];
                if (nr >= 0 && nr < N && nc >= 0 && nc < N) {
                    valid.add(d);
                }
            }
            int dir = valid.get(r.nextInt(valid.size()));
            int np = idx(br + DX[dir], bc + DY[dir]);
            state[bp] = state[np];
            state[np] = 0;
            bp = np;
        }
        return state;
    }

    // ============ 结果类 ============
    static class SAResult {
        boolean solved;
        int steps;
        int[] finalState;
        String message;

        public SAResult(boolean solved, int steps, int[] finalState, String message) {
            this.solved = solved;
            this.steps = steps;
            this.finalState = finalState;
            this.message = message;
        }
    }

    // ============ 主程序 ============
    public static void main(String[] args) {
        SlidingPuzzleSA solver = new SlidingPuzzleSA();

        System.out.println("=== 数字华容道:模拟退火全局优化 ===\n");

        // 生成一个中等难度的随机初始状态(随机游走200步)
        int[] initial = solver.generateRandomState(200);
        System.out.println("初始状态:");
        solver.printBoard(initial);
        System.out.println("初始能量(曼哈顿距离): " + solver.energy(initial));
        System.out.println("可解性: " + solver.isSolvable(initial) + "\n");

        // 执行模拟退火
        long startTime = System.currentTimeMillis();
        SAResult result = solver.solve(initial, 500000, 50.0);
        long elapsed = System.currentTimeMillis() - startTime;

        System.out.println(result.message);
        System.out.println("耗时: " + elapsed + " ms");

        if (result.finalState != null) {
            System.out.println("\n最终状态:");
            solver.printBoard(result.finalState);
            System.out.println("最终能量: " + solver.energy(result.finalState));
        }

        // 对比:与纯随机游走(无Metropolis准则)对比
        System.out.println("\n【对比实验】纯贪心(仅接受能量降低的移动):");
        int[] greedyState = solver.copyState(initial);
        int greedyEnergy = solver.energy(greedyState);
        int greedyBlank = solver.findBlank(greedyState);
        int greedySteps = 0;
        for (int i = 0; i < 100000 && greedyEnergy > 0; i++) {
            // 简化为一次随机合法移动
            int br = solver.row(greedyBlank), bc = solver.col(greedyBlank);
            List<Integer> valid = new ArrayList<>();
            for (int d = 0; d < 4; d++) {
                int nr = br + solver.DX[d], nc = bc + solver.DY[d];
                if (nr >= 0 && nr < N && nc >= 0 && nc < N) valid.add(d);
            }
            int dir = valid.get(new Random(42).nextInt(valid.size()));
            int np = solver.idx(br + solver.DX[dir], bc + solver.DY[dir]);
            int moved = greedyState[np];
            int oldTarget = moved - 1;
            int delta = (Math.abs(br - solver.row(oldTarget)) + Math.abs(bc - solver.col(oldTarget)))
                      - (Math.abs(br + solver.DX[dir] - solver.row(oldTarget)) + Math.abs(bc + solver.DY[dir] - solver.col(oldTarget)));
            if (delta < 0) {
                greedyState[greedyBlank] = moved;
                greedyState[np] = 0;
                greedyBlank = np;
                greedyEnergy += delta;
                greedySteps++;
            }
        }
        System.out.println("贪心最终能量: " + greedyEnergy + " (步数: " + greedySteps + ")");
        System.out.println("模拟退火优势: 能量从 " + solver.energy(initial) + " -> " +
            (result.finalState != null ? solver.energy(result.finalState) : "N/A") +
            ",通过允许劣解跳出局部最优。");
    }
}

关键代码解析

1. 增量能量计算

int oldDist = Math.abs(nr - row(oldTargetPos)) + Math.abs(nc - col(oldTargetPos));
int newDist = Math.abs(br - row(oldTargetPos)) + Math.abs(bc - col(oldTargetPos));
int deltaE = newDist - oldDist;

模拟退火需要大量迭代,若每次重算全部曼哈顿距离,O(n²) 的复杂度将成为瓶颈。通过只计算被移动方块的距离变化,将单次能量计算优化到 O(1),整体性能提升 n² 倍。

2. Metropolis接受准则

double prob = Math.exp(-deltaE / T);
acceptedMove = rand.nextDouble() < prob;

这是模拟退火的核心:即使新状态的_energy更高(更差),仍有概率被接受。温度 T 较高时,接受概率大,算法倾向于全局探索;温度降低后,接受概率趋近于0,算法收敛到局部精细搜索。

3. 自适应降温与回火

if (acceptRate > 0.80) alpha = Math.max(0.9990, alpha - 0.0002);
else if (acceptRate < 0.05) { alpha = Math.min(0.9999, alpha + 0.0003); T *= 1.5; }

根据最近窗口内的接受率动态调整降温速度:
– 接受率过高(>80%):说明搜索过于激进,大多数劣解都被接受,需要加快降温以收敛
– 接受率过低(<5%):说明陷入局部最优,几乎只接受优解,此时”回火”升温以重新获得跳出能力

扩展思考

  1. 能量函数改进:可引入线性冲突(Linear Conflict)惩罚,当两个目标在同一行/列但顺序颠倒时,额外增加惩罚项,使能量函数更接近真实步数下界
  2. 混合策略:先用A或IDA求解较简单子问题获得良好初始解,再用模拟退火精细优化
  3. 并行退火:维护多条马尔可夫链并行搜索,定期交换状态(模拟退火的并行变体)

总结

本文以数字华容道为载体,完整实现了模拟退火算法及其自适应降温策略。关键贡献包括:
– O(1) 增量能量更新,使大规模迭代成为可能
– 基于接受率统计的自适应降温与回火机制
– 可解性判定与随机状态生成方法

模拟退火作为通用全局优化框架,广泛应用于旅行商问题(TSP)、电路布局、蛋白质折叠等NP-hard问题的近似求解。