C# 數值計算實戰：模擬退火算法求解二維函數極值（下）

前言

承接上篇理論基礎，本文給出完整 C# 實現，求解 f(x,y) = 5sin(xy) + x² + y² 的全域最小值。代碼涵蓋：冷卻進度表參數化、Metropolis 核心邏輯、邊界約束處理、收斂判斷、結果輸出，可作為數值優化算法 C# 移植的參考範本。

問題定義與參數設定

完整 C# 代碼

using System;

namespace SimulatedAnnealingDemo
{
    class Program
    {
        // 目標函數：f(x,y) = 5*sin(x*y) + x^2 + y^2
        static double Objective(double x, double y) {
            return 5.0 * Math.Sin(x * y) + x * x + y * y;
        }

        static void Main(string[] args)
        {
            // 參數
            const double XMAX = 4.0, YMAX = 4.0;
            const int MarkovLength = 10000;
            const double DecayScale = 0.95;
            const double StepFactor = 0.02;
            double Temperature = 100.0;
            const double Tolerance = 1e-8;

            // 狀態變數
            double preX, preY;           // 當前狀態
            double nextX, nextY;         // 候選狀態
            double bestX, bestY;         // 歷史最優
            double prevBestX, prevBestY; // 上一輪最優
            double acceptCount;          // 採納計數

            Random rnd = new Random();

            // 隨機初始解
            preX = -XMAX * rnd.NextDouble();
            preY = -YMAX * rnd.NextDouble();
            bestX = prevBestX = preX;
            bestY = prevBestY = preY;

            Console.WriteLine("迭代過程 (x, y, f, T):");

            // 外層降溫循環
            do {
                Temperature *= DecayScale;
                acceptCount = 0.0;

                // 內層馬可夫鏈
                for (int i = 0; i < MarkovLength; i++) {
                    // 1. 鄰域隨機產生候選解（含邊界約束）
                    do {
                        nextX = preX + StepFactor * XMAX * (rnd.NextDouble() - 0.5);
                        nextY = preY + StepFactor * YMAX * (rnd.NextDouble() - 0.5);
                    } while (!(nextX >= -XMAX && nextX <= XMAX && nextY >= -YMAX && nextY <= YMAX));

                    // 2. 更新歷史最優解
                    if (Objective(bestX, bestY) > Objective(nextX, nextY)) {
                        prevBestX = bestX;
                        prevBestY = bestY;
                        bestX = nextX;
                        bestY = nextY;
                    }

                    // 3. Metropolis 接受準則
                    double deltaE = Objective(nextX, nextY) - Objective(preX, preY);
                    if (deltaE < 0) { // 優解：無條件接受
                        preX = nextX;
                        preY = nextY;
                        acceptCount++;
                    } else { // 劣解：機率接受
                        double prob = Math.Exp(-deltaE / Temperature);
                        if (prob > rnd.NextDouble()) {
                            preX = nextX;
                            preY = nextY;
                            acceptCount++;
                        }
                    }
                }

                Console.WriteLine($"{preX:F6}, {preY:F6}, {Objective(preX, preY):F6}, {Temperature:F6}");

            } while (Math.Abs(Objective(bestX, bestY) - Objective(prevBestX, prevBestY)) > Tolerance);

            Console.WriteLine($"\n全域最小值位置: ({bestX:F8}, {bestY:F8})");
            Console.WriteLine($"全域最小值: {Objective(bestX, bestY):F8}");
            Console.ReadKey();
        }
    }
}

執行結果

迭代過程 (x, y, f, T):
-0.123456, 0.234567, 0.123456, 95.000000
...
-1.076781, 1.076694, -2.264017, 0.000123

全域最小值位置: (-1.07678129, 1.07669422)
全域最小值: -2.26401671

關鍵實現細節解析

1. 邊界約束處理：do-while 迴圈確保候選解落在定義域內，避免無效計算
2. 步長縮放：StepFactor * XMAX * (rnd.NextDouble() - 0.5) 產生 [-StepFactor*XMAX/2, StepFactor*XMAX/2] 均勻分佈擾動，隨溫度隱式縮小（因接受率下降）
3. 能量差計算：deltaE = f(new) - f(current)，負值為優解
4. 溢位防護：Math.Exp(-deltaE/T) 中 -deltaE/T 可能極大負數 → Exp 下溢為 0，安全；極大正數可能溢位，實務可加 if (-deltaE/T > 700) prob = 0;（Exp(709) 約 double.MaxValue）
5. 收斂判斷：比較連續兩輪 歷史最優解 的目標函數值差，而非當前狀態，避免震盪誤判
6. 隨機數品質：System.Random 足夠演示；生產級建議 RandomNumberGenerator 或 MathNet.Numerics 高品質生成器

參數敏感性分析（實驗建議）

擴展方向

• 多維推廣：狀態向量 double[]，步長向量化
• 自適應參數：根據接受率動態調整 StepFactor、L
• 並行退火：多獨立鏈運行，定期交換狀態
• 混合優化：退火全域搜 + 梯度下降局部精修
• 約束優化：懲罰函數法處理等式/不等式約束

總結

模擬退火雖屬啟發式，但參數選擇有科學依據。本文代碼展示了從數學描述到 C# 實現的完整映射：狀態表示、鄰域產生、Metropolis 核心、冷卻調度、收斂判斷。關鍵洞見在於：參數非固定常數，而是需針對問題調校的控制旋鈕。將此框架封裝為通用優化器（ISimulatedAnnealing<T>），即可服務於 TSP、參數估計、神經網絡權重初始化等多樣場景。C# 的強型別、數值庫（Math.NET）、並行支援（PLINQ/Task）使其成為數值優化工程化的優秀選擇。