c语言利用模拟退火算法解决最短路径问题

模拟退火算法来源于固体退火原理，将固体加温至充分高，再让其徐徐冷却，加温时，固体内部粒子随温升变为无序状，内能增大，而徐徐冷却时粒子渐趋有序，在每个温度都达到平衡态，最后在常温时达到基态，内能减为最小。根据Metropolis准则，粒子在温度T时趋于平衡的概率为e(-ΔE/(kT))，其中E为温度T时的内能，ΔE为其改变量，k为Boltzmann常数。用固体退火模拟组合优化问题，将内能E模拟为目标函数值f，温度T演化成控制参数t，即得到解组合优化问题的模拟退火算法：由初始解i和控制参数初值t开始，对当前解重复“产生新解→计算目标函数差→接受或舍弃”的迭代，并逐步衰减t值，算法终止时的当前解即为所得近似最优解，这是基于蒙特卡罗迭代求解法的一种启发式随机搜索过程。退火过程由冷却进度表(Cooling Schedule)控制，包括控制参数的初值t及其衰减因子Δt、每个t值时的迭代次数L和停止条件S。

模拟退火算法的模型

1模拟退火算法可以分解为解空间、目标函数和初始解三部分。
2模拟退火的基本思想:
(1) 初始化：初始温度T(充分大)，初始解状态S(是算法迭代的起点)，每个T值的迭代次数L
(2) 对k=1, …, L做第(3)至第6步：
(3) 产生新解S′
(4) 计算增量ΔT=C(S′)-C(S)，其中C(S)为评价函数
(5) 若ΔT<0则接受S′作为新的当前解，否则以概率exp(-ΔT/T)接受S′作为新的当前解.
(6) 如果满足终止条件则输出当前解作为最优解，结束程序。
终止条件通常取为连续若干个新解都没有被接受时终止算法。
(7) T逐渐减少，且T->0，然后转第2步。

模拟退火算法的步骤

模拟退火算法新解的产生和接受可分为如下四个步骤：
第一步是由一个产生函数从当前解产生一个位于解空间的新解；为便于后续的计算和接受，减少算法耗时，通常选择由当前新解经过简单地变换即可产生新解的方法，如对构成新解的全部或部分元素进行置换、互换等，注意到产生新解的变换方法决定了当前新解的邻域结构，因而对冷却进度表的选取有一定的影响。
第二步是计算与新解所对应的目标函数差。因为目标函数差仅由变换部分产生，所以目标函数差的计算最好按增量计算。事实表明，对大多数应用而言，这是计算目标函数差的最快方法。
第三步是判断新解是否被接受,判断的依据是一个接受准则，最常用的接受准则是Metropolis准则: 若ΔT<0则接受S′作为新的当前解S，否则以概率exp(-ΔT/T)接受S′作为新的当前解S。
第四步是当新解被确定接受时，用新解代替当前解，这只需将当前解中对应于产生新解时的变换部分予以实现，同时修正目标函数值即可。此时，当前解实现了一次迭代。可在此基础上开始下一轮试验。而当新解被判定为舍弃时，则在原当前解的基础上继续下一轮试验。

模拟退火算法与初始值无关，算法求得的解与初始解状态S(是算法迭代的起点)无关；模拟退火算法具有渐近收敛性，已在理论上被证明是一种以概率l 收敛于全局最优解的全局优化算法；模拟退火算法具有并行性。

#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <cstring>
#include <ctime>

#define T_start 5000.0      //初始温度
#define T_end (1e-8)        //结束温度 //1e-8：1乘以十的负八次方
#define q 0.98              //退火系数
#define L 1000              //每个温度最大迭代次数
#define N 34                //城市个数
int city_result[N];         //城市列表的解空间
double city[N][2] = {{9932, 4439}, {10109, 4351}, {11552, 3472}, {10302, 3290}, {8776, 3333}, {7040, 4867}, {9252, 4278}, {9395, 4539}, {11101, 2540}, {9825, 5087}, {10047, 4879}, {10227, 4648}, {100027, 4229}, {9878, 4211}, {9087, 4065}, {10438, 4075}, {10382, 3865}, {11196, 3563}, {11075, 3543}, {11544, 3365}, {11915, 2900}, {11305, 3189}, {11073, 3137}, {10950, 3394}, {11576, 2575}, {12239, 2785}, {11529, 2226}, {9328, 4006}, {10012, 3811}, {9952, 3410}, {10612, 2954}, {10349, 2784}, {11747, 2469}, {11673, 2461}}; //中国34个城市实际距离坐标x坐标y坐标
//函数声明
//double city[N][2] = {{0,0},{3,0},{3,4}};//测试数据
double distance(double *city1, double *city2); //计算两城市间距离
double path(int city_result[N]);               //计算总路径
void init();                                   //初始化解空间
void creat();                                  //生成新的解空间

int main()
{
    time_t start, end;      //记录程序开始结束时间
    double time_sum;        //记录程序运行时间
    start = clock();        //程序开始时间
    int i, count = 0;       //降温计数器
    int city_copyresult[N]; //拷贝解空间
    double path1, path2;    //原有解空间，新解空间的总路径
    double dE;              //原有解空间与新解空间的差值
    double r;               //随机产生0~1的值，是否接受新的解
    double T;               //当前温度

    srand(time(nullptr));//设置一个随机种子，每次运行都能保证随机种子不同，rand()函数可以用来产生随机数，但是这不是真正意义上的随机数，是一个伪随机数，它是根据一个数，我们可以称它为种子，为基准以某个递推公式推算出来的一系数，但这不是真正的随机数，当计算机正常开机后，这个种版子的值是定了的，除非你破坏了系统，为了改变这个种子的值
    init(); //初始化解空间

    T = T_start;//初始温度赋值

    while (T > T_end) //当前温度大于结束温度
    {
        for (i = 0; i < L; i++)
        {
            memcpy(city_copyresult, city_result, N * sizeof(int));
            creat(); //产生新的解空间
            path1 = path(city_copyresult);
            path2 = path(city_result);
            dE = path2 - path1;
            if (dE > 0) //Metropolis准则，以一定概率接受其作为新的解
            {
                r = rand() / (RAND_MAX);
                if (exp(-dE / T) <= r) //exp函数：以e为底的 -dE / T 次方 //高温状态下，可以接受能量差值较大的新状态；低温状态下，则只能接受能量差值较小的新状态
                    memcpy(city_result, city_copyresult, N * sizeof(int)); //保留原来的解
                    //void *memcpy(void *dest, const void *src, size_t n) 它的功能是从src的开始位置拷贝n个字节的数据到dest。如果dest存在数据，将会被覆盖。memcpy函数的返回值是dest的指针。memcpy函数定义在string.h头文件里。
            }
        }
        T *= q; //降温
        count++;
    }

    end = clock();  //程序结束时间
    time_sum = (double)(end - start ) / (CLOCKS_PER_SEC);    //换算成秒
    printf("共降温：%d次\n", count);
    printf("经过模拟退火算法得出最优路径长度为：%f\n", path(city_result));
    for (i = 0; i < N; i++)
    {
        printf("%d->", city_result[i]);
    }
    printf("%d\n", city_result[0]);
    printf("程序共耗时%f秒.\n",time_sum);
    system("pause");
    return 0;
}

double distance(double *city1, double *city2) //计算两个城市之间的距离
{
    double x1, x2, y1, y2, dis;
    x1 = *city1;       //第一个城市的x坐标
    x2 = *city2;       //第二个城市的x坐标
    y1 = *(city1 + 1); //第一个城市的y坐标
    y2 = *(city2 + 1); //第二个城市的y坐标
    dis = sqrt(pow((x2 - x1), 2) + pow((y2 - y1), 2));
    return dis;
}

double path(int city_result[N]) //计算总路径
{
    int i, city1_num, city2_num; //解空间中的两个城市序号
    double sum = 0;              //路径总长度
    for (i = 0; i < N - 1; i++)  //解空间中首位到末位的总路径
    {
        city1_num = city_result[i];
        city2_num = city_result[i + 1];
        sum += distance(city[city1_num], city[city2_num]);
    }
    sum += distance(city[0], city[N - 1]); //加上解空间中末位到首位的路径
    return sum;
}

void init() //初始化解空间
{
    int i;
    for (i = 0; i < N; i++)
        city_result[i] = i; //顺序生成解空间
}

void creat() //生成新的解空间
{
    int point1, point2, temp;
    point1 = rand() % N;    //产生1~N的一个随机整数
    point2 = rand() % N;
    //交换city_result[point1]和city_result[point2]
    temp = city_result[point1];
    city_result[point1] = city_result[point2];
    city_result[point2] = temp;
}

本文参考：
https://baike.baidu.com/item/模拟退火算法/355508?fromtitle=%E9%80%80%E7%81%AB%E7%AE%97%E6%B3%95&fromid=6081525&fr=aladdin

https://blog.csdn.net/zhuzilong2013/article/details/82821348?utm_medium=distribute.pc_relevant_t0.none-task-blog-BlogCommendFromMachineLearnPai2-1.nonecase&depth_1-utm_source=distribute.pc_relevant_t0.none-task-blog-BlogCommendFromMachineLearnPai2-1.nonecase