高性能计算

模拟退火算法

TSP问题

TSP（Traveling Salesman Problem）问题，即旅行商问题。
假设有一个旅行商人要拜访N个城市，他必须选择所要走的路径，路径的限制是每个城市只能拜访一次，而且最后要回到原来出发的城市。路径的选择目标是要求得的路径路程为所有路径之中的最小值。
它是图论中最著名的问题之一，即“已给一个n个点的完全图，每条边都有一个长度，求总长度最短的经过每个顶点正好一次的封闭回路。
TSP问题是一个组合优化问题，该问题被证明具有NPC计算复杂性。

什么是NPC计算复杂性：

NPC问题

退火思路

模拟退火算法是基于Monte-Carlo迭代求解策略的一种随机寻优算法。
其出发点是基于物理中固体物质的退火过程与一般组合优化问题之间的相似性。模拟退火算法从某一较高初温出发，伴随温度参数的不断下降，结合概率突跳特性在解空间中随机寻找目标函数的全局最优解，即在局部最优解能概率性地跳出并最终趋于全局最优。
模拟退火算法是一种通用的优化算法，理论上算法具有概率的全局优化性能。

模拟退火算法的基础是Metropolis采样：
Metropolis 采样的基本思路是，从一个已知的形式较为简单的分布中采样，并以一定的概率接受这个样本作为目标分布的近似样本。
根据Metropolis准则，粒子在温度T时趋于平衡的概率为e(-ΔE/(kT))，其中E为温度T时的内能，ΔE为其改变量，k为Boltzmann常数。
由初始解i和控制参数初值t开始，对当前解重复“产生新解→计算目标函数差→接受或舍弃”的迭代，并逐步衰减t值，算法终止时的当前解即为所得近似最优解。

模拟退火串行算法参考链接：

模拟退火算法（SA）求解TSP 问题（C语言）

本文以下代码根据上述链接进行改进（并行算法）。

主要函数

distance()函数用来计算两个结点（城市）间的距离，参数为两个结点（城市）的坐标。

double distance(double* city1, double* city2)
{
    double x1 = *city1;
    double y1 = *(city1 + 1);
    double x2 = *(city2);
    double y2 = *(city2 + 1);
    double dis = sqrt((x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2));
    return dis;
}

path_len()函数用来计算路径长度，参数为解数组（即一组城市的下标序列）。

double path_len(int* arr)
{
    double path = 0;  // 初始化路径长度
    int index = *arr; // 定位到第一个数字(城市序号)
    for (int i = 0; i < NCITY - 1; i++)
    {
        int index1 = *(arr + i);
        int index2 = *(arr + i + 1);
        double dis = distance(city_pos[index1 - 1], city_pos[index2 - 1]);
        path += dis;
    }
    int last_index = *(arr + NCITY - 1); // 最后一个城市序号
    int first_index = *arr;              // 第一个城市序号
    double last_dis = distance(city_pos[last_index - 1], city_pos[first_index - 1]);
    path = path + last_dis;
    // 返回总的路径长度
    return path; 
}

init()为初始化函数，用来初始化变量值。

void init()
{
    T = T0;//初始温度
    count = 0;//初始降温次数
    best_Val = 0;//初始路径长度
    for (int i = 0; i < NCITY; i++)
    {
        //初始化一个解
        best_Road[i] = i + 1;
        tempCity_list[i] = i + 1;
        subProcess_Road[i] = i + 1;
    }
}

create_new()函数用来产生一个新解———采用随机交叉两个位置的方式产生新的解

void create_new()
{
    double r1 = ((double)rand()) / (RAND_MAX + 1.0);
    double r2 = ((double)rand()) / (RAND_MAX + 1.0);
    int pos1 = (int)(NCITY * r1); //第一个交叉点的位置
    int pos2 = (int)(NCITY * r2); //第二个交叉点的位置
    // 交换两个点
    int temp = tempCity_list[pos1];
    tempCity_list[pos1] = tempCity_list[pos2];
    tempCity_list[pos2] = temp;
}

模拟退火算法部分片段

for (int i = 0; i < L; i++)
{
	//复制数组
	memcpy(tempCity_list, subProcess_Road, NCITY * sizeof(int));
	//产生新解
	create_new();
	f1 = path_len(subProcess_Road);
	f2 = path_len(tempCity_list);
	df = f2 - f1;
	//以下是Metropolis准则
	if (df >= 0)
	{
		subProcess_Val = f2;
		memcpy(subProcess_Road, tempCity_list, NCITY * sizeof(int));
		r = ((double)rand()) / (RAND_MAX);
		if (exp(-df / T) <= r) //保留原来的解
		{
			subProcess_Val = f1;
			memcpy(tempCity_list, subProcess_Road, NCITY * sizeof(int));
		}
		else
		{
			subProcess_Val = f2;
			memcpy(subProcess_Road, tempCity_list, NCITY * sizeof(int));
		}
	}
}

并行思路

思路One（本文）

使用两个节点的多个进程，主进程会将所有城市的坐标信息分发给每个次结点，之后每个结点都运行模拟退火算法获取当前进程的城市“最短”路径，每个节点同样可以使用OpenMPI在内部进行并行。
每隔一定迭代次数（初步设置为10000），主进程会设置“栅栏”收集每个结点当前的“最短”路径，并将收集到的“最短”路径信息分发到每个结点中继续进行迭代。
当温度达到终止条件，收集到最后获得的“最短”路径。

思路Two（Else）

城市坐标数据集分成两份，主节点将其中一份发给次节点，分别由两个节点进行计算各部分的相对最短路径，在单个节点中使用OpenMPI对退火过程中的循环进行并行。
两台机器都算完后，结果汇聚到主节点，将两个路径通过最短的边连接形成一个路径（主节点参与运算），该路径即为最短路径。
在次节点发送结果与主节点接收结果之前设置栅栏。

思路Two文章链接: 用MPI并行模拟退火解决TSP问题的程序

全部代码

#include<stdio.h>
#include<stdlib.h>
#include<string.h>
#include<time.h>
#include<math.h>
#include "mpi.h"

#define NCITY 31        //城市数量
#define T0 50000.0  	//初始温度
#define T_end (1e-8)	//终止温度
#define q  0.98     	//退火系数
#define L 1000      	//每个温度的迭代次数，即链长


double T;                   //当前温度
int count;                  //降温次数
int tempCity_list[NCITY];   //临时路径
time_t start, finish;       //程序用时

double best_Val;            //当前最优路径值
int best_Road[NCITY];       //当前最优路径

double subProcess_Val;      //子进程处最优路径值
int subProcess_Road[NCITY]; //子进程处最优路径

// 中国城市坐标
double city_pos[NCITY][2] =
{
    {1304,2312},{3639,1315},{4177,2244},{3712,1399},
    {3488,1535},{3326,1556},{3238,1229},{4196,1004},
    {4312,7901},{4386,5702},{3007,1970},{2562,1756},
    {2788,1491},{2381,1676},{1332,6953},{3715,1678},
    {3918,2179},{4061,2370},{3780,2212},{3676,2578},
    {4029,2838},{4263,2931},{3429,1908},{3507,2367},
    {3394,2643},{3439,3201},{2935,3240},{3140,3550},
    {2545,2357},{2778,2826},{2370,2975}
};

// 距离函数
double distance(double* city1, double* city2)
{
    double x1 = *city1;
    double y1 = *(city1 + 1);
    double x2 = *(city2);
    double y2 = *(city2 + 1);
    double dis = sqrt((x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2));
    return dis;
}

// 计算路径长度
double path_len(int* arr)
{
    double path = 0; // 初始化路径长度
    int index = *arr; // 定位到第一个数字(城市序号)
    for (int i = 0; i < NCITY - 1; i++)
    {
        int index1 = *(arr + i);
        int index2 = *(arr + i + 1);
        double dis = distance(city_pos[index1 - 1], city_pos[index2 - 1]);
        path += dis;
    }
    int last_index = *(arr + NCITY - 1); // 最后一个城市序号
    int first_index = *arr; // 第一个城市序号
    double last_dis = distance(city_pos[last_index - 1], city_pos[first_index - 1]);
    path = path + last_dis;
    return path; // 返回总的路径长度
}

// 初始化函数
void init()
{
    T = T0;//初始温度
    count = 0;//初始降温次数
    best_Val = 0;//初始路径长度
    for (int i = 0; i < NCITY; i++)
    {
        //初始化一个解
        subProcess_Road[i] = i + 1;
        best_Road[i] = i + 1;
        tempCity_list[i] = i + 1;
    }
}

// 产生一个新解--采用随机交叉两个位置的方式产生新的解
void create_new()
{
    double r1 = ((double)rand()) / (RAND_MAX + 1.0);
    double r2 = ((double)rand()) / (RAND_MAX + 1.0);
    int pos1 = (int)(NCITY * r1); //第一个交叉点的位置
    int pos2 = (int)(NCITY * r2);
    int temp = tempCity_list[pos1];
    tempCity_list[pos1] = tempCity_list[pos2];
    tempCity_list[pos2] = temp;   // 交换两个点
}

// 主函数
int main(int argc, char* argv[])
{
    //初始化随机数种子
    srand((unsigned)time(NULL));

    int myid, NPROCESS, namelen;
    char processor_name[MPI_MAX_PROCESSOR_NAME];
    MPI_Status status;
    MPI_Init(&argc, &argv);			         		 	//初始化MPI环境
    MPI_Comm_rank(MPI_COMM_WORLD, &myid);	         	//获取当前进程编号
    MPI_Comm_size(MPI_COMM_WORLD, &NPROCESS);	     	//获取进程的个数
    MPI_Get_processor_name(processor_name, &namelen);   //获得本进程的机器名
    //processor_name为返回的机器名字符串，namelen为返回的机器名长度

    //用来保存机器名
    char **nameofm = (char**)malloc(sizeof(char*)*NPROCESS);
    //从子进程处接收的最优路径值
    double *recover_Val = (double*)malloc(sizeof(double)*NPROCESS);
    //从子进程处接收的最优路径
    int **recover_Road = (int**)malloc(sizeof(int*)*NPROCESS);
    for(int i=0;i<NPROCESS;i++)
    {
    	nameofm[i] = (char*)malloc(sizeof(char) * 128);
    	recover_Road[i] = (int*)malloc(sizeof(int) * NCITY);
    }

    //Test--第一次通迅！
    if (myid == 0)
    {
        for (int i = 1; i < NPROCESS; i++)
        {
            MPI_Recv(nameofm[i], 128, MPI_CHAR, i, 99, MPI_COMM_WORLD, &status);
            printf("我是机器:%s的主进程%d，收到子进程%d。\n", nameofm[i], myid, i);
        }
    }
    else
    {
        printf("我是子进程：%d\n", myid);
        MPI_Send(processor_name, sizeof(processor_name), MPI_CHAR, 0, 99, MPI_COMM_WORLD);
    }

    //城市坐标初始化操作
    if (myid == 0)
    {
    	//程序开始计时
        start = MPI_Wtime();
        //主进程（0）将城市坐标送至各个子进程
        for (int i = 1; i < NPROCESS; i++)
        {
            MPI_Send(*city_pos, NCITY * 2, MPI_DOUBLE, i, 99, MPI_COMM_WORLD);
        }
        init();//初始温度和初始解
    }
    else
    {
        //各个子进程从主进程处（0）接受城市坐标
        MPI_Recv(*city_pos, NCITY * 2, MPI_DOUBLE, 0, 99, MPI_COMM_WORLD, &status);
        init();//初始温度和初始解
    }

    //f1为初始解目标函数值，f2为新解目标函数值，df为二者差值
    double f1, f2, df;
    //0-1之间的随机数，用来决定是否接受新解
    double r;
    while (T > T_end)//当温度低于结束温度时，退火结束
    {
        //主进程每10000代从子进程中获取当前最优路径并分发下去
        //子进程通过模拟退火算法获取当前进程内最优路径
        if (myid == 0)
        {
            if ((count % 10 == 0 || T * q < T_end) && count != 0)
            {
                for (int i = 0; i < NPROCESS - 1; i++)
                {
                    //从子进程处回收最佳路径值及路径，分别存放在recover_Val & recover_Road
                    MPI_Recv(&recover_Val[i], 1, MPI_DOUBLE, i + 1, 99, MPI_COMM_WORLD, &status);
                    MPI_Recv(recover_Road[i], NCITY, MPI_INT, i + 1, 99, MPI_COMM_WORLD, &status);
                }
                //选出子进程中最好的路径值及路径
                if (best_Val == 0)
                {
                    best_Val = recover_Val[0];
                }
                for (int i = 0; i < NPROCESS - 1; i++)
                {
                    if (recover_Val[i] <= best_Val)
                    {
                        best_Val = recover_Val[i];
                        memcpy(best_Road, recover_Road[i], NCITY * sizeof(int));
                    }
                }
                //将最佳路径值及路径分发给子进程
                for (int i = 0; i < NPROCESS - 1; i++)
                {
                    MPI_Send(&best_Val, 1, MPI_DOUBLE, i + 1, 99, MPI_COMM_WORLD);
                    MPI_Send(best_Road, NCITY, MPI_INT, i + 1, 99, MPI_COMM_WORLD);
                }
                //输出当前最优路径
                printf("每10000代收集的“最优”路径：");
                for (int i = 0; i < NCITY - 1; i++)
                {
                    printf("%d--->", best_Road[i]);
                }
                printf("%d\n\n", best_Road[NCITY - 1]);
            }
        }
        else
        {
            for (int i = 0; i < L; i++)
            {
                //复制数组
                memcpy(tempCity_list, subProcess_Road, NCITY * sizeof(int));
                //产生新解
                create_new();
                f1 = path_len(subProcess_Road);
                f2 = path_len(tempCity_list);
                df = f2 - f1;
                //以下是Metropolis准则
                if (df >= 0)
                {
                    subProcess_Val = f2;
                    memcpy(subProcess_Road, tempCity_list, NCITY * sizeof(int));
                    r = ((double)rand()) / (RAND_MAX);
                    if (exp(-df / T) <= r) //保留原来的解
                    {
                        subProcess_Val = f1;
                        memcpy(tempCity_list, subProcess_Road, NCITY * sizeof(int));
                    }
                }
                else
                {
                    subProcess_Val = f2;
                    memcpy(subProcess_Road, tempCity_list, NCITY * sizeof(int));
                }
            }
            //子进程每隔10000代，向主进程发送自己的最优路径值和最优路径，再接受所有子进程中最优的
            if ((count % 10 == 0 || T * q < T_end) && count != 0)
            {
                //发送
                MPI_Send(&subProcess_Val, 1, MPI_DOUBLE, 0, 99, MPI_COMM_WORLD);
                MPI_Send(subProcess_Road, NCITY, MPI_INT, 0, 99, MPI_COMM_WORLD);
                //接受
                MPI_Recv(&subProcess_Val, 1, MPI_DOUBLE, 0, 99, MPI_COMM_WORLD, &status);
                MPI_Recv(subProcess_Road, NCITY, MPI_INT, 0, 99, MPI_COMM_WORLD, &status);
            }
        }
        T *= q; //降温
        count++;
    }
    if (myid == 0)
    {
        //输出最优路径
        printf("最终主进程收集的最优路径：Best_Road：");
        for (int i = 0; i < NCITY - 1; i++)
        {
            printf("%d--->", best_Road[i]);
        }
        printf("%d\n", best_Road[NCITY - 1]);
        printf("最优路径长度为:%lf\n\n", best_Val);
        //计算时间
        finish = MPI_Wtime();
        double duration = ((double)(finish - start)) / CLOCKS_PER_SEC *1000;
        printf("模拟退火算法，初始温度T0=%.2f,降温系数q=%.2f,每个温度迭代%d次,共降温%d次.\n", T0, q, L, count);
        printf("程序运行耗时:%lf秒.\n", duration);
    }

    MPI_Finalize();//结束MPI程序的运行
    return 0;
}