模拟退火算法实例（c++ 与 c# 实现）

此片文章主要参考CSDN博主里头的一篇文章， 将自己的理解写下来，以方便后期的查阅。

 一、C++ 实现

1. 已知平面上若干点坐标（xi, yi）, 求平面上一点p(x, y) , 到这些点的总距离最小。

思路： 取所有点的均值为目标点。计算全部点与目标点求差值的和，将目标点以一定系数朝着总和的方向移动，得到新的目标点。

// 求最小距离
// 限制条件： 1 <= n <= 100, 0<= xi, yi <= 1e4
#include<iostream>
#include<cstdio>
#include<cdtdlib>
#include<cmatch>

using namespace std;

struct Pt{
    double x, y;
}P[105]

double sqr(double x) {
    return x*x
}

// get distance between two point
double dist(Pt a, Pt b)
{
    return sqrt(sqr(a.x - b.x) + sqr(a.y - b.y));
}

double get_sum(Pt p0, int n)
{
    double ret = 0;
    for(int i = 0; i < n; ++i)
        ret += dist(p0, p[i]);


int main
{
    // 设置初始化点数据
    // p[n] = { ... ,..., ....}
    // 取所有点的平均位置，作为最近点的位置
    double x0=0, y0 =0;
    for(int i = 0; i < n; ++i)
    {
        x0 += p[i].x;
        y0 += p[i].y;
    }
    x0 /= n;
    y0 /= n;

    double ans = get_sum((pt){x0, y0}, n); // 当前 目标值
    double temp = le5;  // 初始化温度， 根据需要设定

    while(temp > 0.02)  // 0.02 为温度的下限， 若温度为 temp 达到下限， 则停止搜索
    {
        double x = 0, y = 0;
        for(int i = 0; i < n; ++i) { // 获取步长的规则根据要求设定
            x += (p[i].x - x0) / dist((Pt){x0, y0}, p[i]);
            y += (p[i].x - y0) / dist((Pt){x0, y0}, p[i]);
        }
        // 变化后，新的目标值
        // 此处变化的系数应该是逐渐减小的， temp 逐渐减小，符合要求
        double tmp = get_sum((Pt){x0 + x * temp, y0 + y * temp}, n); 

        // 进行判断
        if(temp < ans)
        {
            ans = temp;
            x0 += x * temp;
            y0 += y * temp;
        }
        // 退火算法的精髓 ： 当不满足移动条件时， 也按照一定的概率进行移动
        // 注意 ： 移动的概率应该逐渐减小  
        // e n次幂， n 应该小于0
        //  假设 random() 的作用 ： 产生 0- 1 之间的随机数
        else if(Math.exp((ans - temp) / temp) > random())
        {
            ans = temp;
            x0 += x * temp;
            y0 += y * temp;
        }
        temp *= 0.98; // 0.98 为降火速率（范围为0~1， 数字越大，得到的全局最优解概率越高，运行时间越长）
    }
    printf("The minimal dist is : ");
    printf("%.0f
", ans);
}

二、C# 实现代码

已知空间上若干点（xi, yi, zi）, 求空间上包含这些点的最小球半径 R, 以及球心坐标。

思路：球心与这些点的最大距离为半径， 球心与最大距离点生成向量，将球心朝着该向量方向移动若干距离，再计算半径的变化。

namespace Test_BST
{
    public class Program
    {
        static void Main(string[] args)
        {
            // 初始化输入点
            List<Point> originPoints = new List<Point>() { ............};
            double radius = AnnealAlgorithm(originPoints);
        }

        private struct Point
        {
            public double x;
            public double y;
            public double z;
        }

        // square of a number
        private static double Sqr(double x) { return x * x; }

        // 两点之间的距离
        private static double Dist(Point A, Point B)
        {
            return Math.Sqrt(Sqr(A.x - B.x) + Sqr(A.y - B.y) + Sqr(A.z - B.z));
        }

        // 求最大半径
        private static double GetMaxRadius(Point p0, List<Point> pts)
        {
            double maxRadius = 0;
            foreach (var point in pts)
            {
                double radius = Dist(p0, point);
                maxRadius = radius > maxRadius ? radius : maxRadius;
            }

            return maxRadius;
        }

        private static double AnnealAlgorithm(List<Point> originPts)
        {
            Point center = new Point();
            center.x = 0;
            center.y = 0;
            center.z = 0;

            // 将初始化中心点设置为所有点的代数平均位置
            foreach (var pt in originPts)
            {
                center.x += pt.x;
                center.y += pt.y;
                center.z += pt.z;
            }
            center.x /= originPts.Count;
            center.y /= originPts.Count;
            center.z /= originPts.Count;

            double temp = 1e3; // 初始温度
            double coolingFactor = 0.98; // 降温因子
            double ans = GetMaxRadius(center, originPts); // 当前最小半径
            var random = new Random();

            while (temp > 1e-5)
            {
                Point newCenter = new Point();
                double max_r = 0;
                // 找到与当前中心点距离最远的点，将中心向着改点移动
                for (int i = 0; i < originPts.Count; i++)
                {
                    double r = Dist(center, originPts[i]);
                    if (r > max_r)
                    {
                        newCenter.x = (originPts[i].x - center.x) / r;
                        newCenter.y = (originPts[i].y - center.y) / r;
                        newCenter.z = (originPts[i].z - center.z) / r;
                        max_r = r;
                    }
                }
                newCenter.x = center.x + newCenter.x * temp;
                newCenter.y = center.y + newCenter.y * temp;
                newCenter.z = center.z + newCenter.z * temp;

                // 移动后的最大半径
                double tmp = GetMaxRadius(newCenter, originPts);

                if (tmp < ans) 
                {
                    center.x += newCenter.x * temp;
                    center.y += newCenter.y * temp;
                    center.z += newCenter.z * temp;
                }
                else if (Math.Exp((ans -tmp)/temp) > random.NextDouble() )
                {
                    center.x += newCenter.x * temp;
                    center.y += newCenter.y * temp;
                    center.z += newCenter.z * temp;
                }

                temp *= coolingFactor;
            }
            double miniRadius = GetMaxRadius(center, originPts);
            Console.WriteLine("the cooridnate of the center is {0}, the radius value is {1}", center, miniRadius));

            return miniRadius;
        }
    }
}

 参考： http://blog.csdn.net/whai362/article/details/46980471#comments