模擬退火演算法
阿新 • • 發佈：2018-12-04
模擬退火演算法是一種求函式最小值點的隨機演算法，最近工作中要用到優化演算法，因此研究了一下這個比較簡單的演算法。模擬退火最基本的思想來源於金屬退火過程，在退火過程中，熱運動的原子逐漸凍結在勢能最低的位置，從而保證了整塊金屬晶格結構的一致性，達到最佳效能。模擬退火演算法的基本步驟如下：
1. 確定搜尋起始點x0和E(x0)，其中E(x)為要求其最小值的函式。
2. 從x0出發在相空間中隨機走一步到x點（可固定步長，也可不斷減小步長加快收斂性）。
3. 比較E(x)和E(x0)的大小，如果E(x)<E(x0)，則令x0=x，重複2；否則取一個隨機數r，如果r < exp[-(E-E0)/(k*T)]，則令x0=x，重複2；如果r > exp[-(E-E0)/(k*T)]，則直接重複2。其中E=E(x), E0=E(x0)。
4. 2和3步迴圈n次後，令溫度T下降一點，繼續重複2-3步。如果溫度T小於某個值T0，則整個過程結束。
在以上過程中，始終記錄當前的最小值Emin和對應的xmin，在過程結束後給出Emin作為函式最小值，xmin作為函式取最小值的點作為整個問題的近似解。
在實際計算中，需要根據問題的特點，適當調節步長的取值、迴圈次數n和溫度下降的比例，從而得到比較好的效果。下面給出模擬退火演算法的程式碼，該程式碼改編自GSL中的Simulated Annealing程式碼，為了在windows下執行做了部分改動，並將隨機數生成部分的程式碼單獨提出來，保證整個程式碼的獨立性，不依賴於VS之外的其他庫。
siman.h
#ifndef __GSL_SIMAN_H__ #define __GSL_SIMAN_H__ #include <stdlib.h> #include <stdio.h> #ifdef THIS_IS_DLL #define PREFIXAPI __declspec(dllexport) #else #define PREFIXAPI __declspec(dllimport) #endif #undef __BEGIN_DECLS #undef __END_DECLS #ifdef __cplusplus # define __BEGIN_DECLS extern "C" { # define __END_DECLS } #else # define __BEGIN_DECLS /* empty */ # define __END_DECLS /* empty */ #endif __BEGIN_DECLS typedef struct { const char* name; unsigned long int max; unsigned long int min; size_t size; void (*set) (void *state, unsigned long int seed); unsigned long int(*get) (void *state); double(*get_double) (void *state); } gsl_rng_type; typedef struct { const gsl_rng_type *type; void *state; } gsl_rng; PREFIXAPI extern const gsl_rng_type *gsl_rng_default; unsigned long int gsl_rng_default_seed; PREFIXAPI const gsl_rng_type *gsl_rng_env_setup(const gsl_rng_type *r, unsigned long int s); PREFIXAPI gsl_rng *gsl_rng_alloc(const gsl_rng_type *T); PREFIXAPI void gsl_rng_free(gsl_rng *r); /// 獲得整型隨機數的函式 __inline unsigned long int gsl_rng_get(const gsl_rng *r) { return (r->type->get) (r->state); } /// 獲得實型隨機數的函式 __inline double gsl_rng_uniform(const gsl_rng *r) { return (r->type->get_double) (r->state); } /* types for the function pointers passed to gsl_siman_solve */ typedef double (*gsl_siman_Efunc_t) (void *xp); typedef void (*gsl_siman_step_t) (const gsl_rng *r, void *xp, double step_size); typedef double (*gsl_siman_metric_t) (void *xp, void *yp); typedef void (*gsl_siman_print_t) (void *xp); typedef void (*gsl_siman_copy_t) (void *source, void *dest); typedef void * (*gsl_siman_copy_construct_t) (void *xp); typedef void (*gsl_siman_destroy_t) (void *xp); /* this structure contains all the information needed to structure the search, beyond the energy function, the step function and the initial guess. */ typedef struct { int n_tries; /* how many points to try for each step */ int iters_fixed_T; /* how many iterations at each temperature? */ double step_size; /* max step size in the random walk */ /* the following parameters are for the Boltzmann distribution */ double k, t_initial, mu_t, t_min; } gsl_siman_params_t; /* prototype for the workhorse function */ PREFIXAPI void gsl_siman_solve(const gsl_rng * r, void *x0_p, gsl_siman_Efunc_t Ef, gsl_siman_step_t take_step, gsl_siman_metric_t distance, gsl_siman_print_t print_position, gsl_siman_copy_t copyfunc, gsl_siman_copy_construct_t copy_constructor, gsl_siman_destroy_t destructor, size_t element_size, gsl_siman_params_t params); PREFIXAPI void gsl_siman_solve_many (const gsl_rng * r, void *x0_p, gsl_siman_Efunc_t Ef, gsl_siman_step_t take_step, gsl_siman_metric_t distance, gsl_siman_print_t print_position, size_t element_size, gsl_siman_params_t params); __END_DECLS #endif /* __GSL_SIMAN_H__ */
siman.c
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <string.h>
#include <assert.h>
#define THIS_IS_DLL
#include "siman.h"
#define GSL_LOG_DBL_MIN (-7.0839641853226408e+02)
const gsl_rng_type *gsl_rng_default;
/// 隨機數初始化，實際上為底層隨機數生成器包了一層
const gsl_rng_type *gsl_rng_env_setup(const gsl_rng_type *r, unsigned long int s)
{
unsigned long int seed = s;
gsl_rng_default = r;
gsl_rng_default_seed = seed;
return gsl_rng_default;
}
/// 真正的隨機數初始化
void gsl_rng_set(const gsl_rng *r, unsigned long int seed)
{
(r->type->set) (r->state, seed);
}
gsl_rng *gsl_rng_alloc(const gsl_rng_type *T)
{
gsl_rng *r = (gsl_rng *)malloc(sizeof(gsl_rng));
if (r == 0) printf("failed to allocate space for rng struct");
r->state = calloc(1, T->size);
if (r->state == 0)
{
free(r);
printf("failed to allocate space for rng state");
}
r->type = T;
gsl_rng_set(r, gsl_rng_default_seed);
return r;
}
void gsl_rng_free(gsl_rng *r)
{
if (!r) return;
free(r->state);
free(r);
}
static __inline double
boltzmann(double E, double new_E, double T, gsl_siman_params_t *params)
{
double x = -(new_E - E) / (params->k * T);
/* avoid underflow errors for large uphill steps */
return (x < GSL_LOG_DBL_MIN) ? 0.0 : exp(x);
}
static __inline void
copy_state(void *src, void *dst, size_t size, gsl_siman_copy_t copyfunc)
{
if (copyfunc) {
copyfunc(src, dst);
} else {
memcpy(dst, src, size);
}
}
/* implementation of a basic simulated annealing algorithm */
void
gsl_siman_solve (const gsl_rng * r, void *x0_p, gsl_siman_Efunc_t Ef, // Ef是待求函式，x0_p是起始點，r是隨機數生成器
gsl_siman_step_t take_step, // 用於隨機走一步
gsl_siman_metric_t distance,
gsl_siman_print_t print_position, // 用於將中間結果打印出來
gsl_siman_copy_t copyfunc, // 用於實現變步長的模擬退火，將一個點的值複製給另一個點
gsl_siman_copy_construct_t copy_constructor, // 用於變步長的模擬退火，生成一個新座標點
gsl_siman_destroy_t destructor, // 回收座標點記憶體
size_t element_size,
gsl_siman_params_t params)
{
void *x, *new_x, *best_x;
double E, new_E, best_E;
int i;
double T, T_factor;
int n_evals = 1, n_iter = 0, n_accepts, n_rejects, n_eless;
/* this function requires that either the dynamic functions (copy,
copy_constructor and destrcutor) are passed, or that an element
size is given */
assert((copyfunc != NULL && copy_constructor != NULL && destructor != NULL)
|| (element_size != 0));
distance = 0 ; /* This parameter is not currently used */
E = Ef(x0_p);
if (copyfunc) {
x = copy_constructor(x0_p);
new_x = copy_constructor(x0_p);
best_x = copy_constructor(x0_p);
} else {
x = (void *) malloc (element_size);
memcpy (x, x0_p, element_size);
new_x = (void *) malloc (element_size);
best_x = (void *) malloc (element_size);
memcpy (best_x, x0_p, element_size);
}
best_E = E;
T = params.t_initial;
T_factor = 1.0 / params.mu_t;
if (print_position) {
printf ("#-iter #-evals temperature position energy\n");
}
while (1) {
n_accepts = 0;
n_rejects = 0;
n_eless = 0;
for (i = 0; i < params.iters_fixed_T; ++i) {
copy_state(x, new_x, element_size, copyfunc);
take_step (r, new_x, params.step_size);
new_E = Ef (new_x);
if(new_E <= best_E){
if (copyfunc) {
copyfunc(new_x,best_x);
} else {
memcpy (best_x, new_x, element_size);
}
best_E=new_E;
}
++n_evals; /* keep track of Ef() evaluations */
/* now take the crucial step: see if the new point is accepted
or not, as determined by the boltzmann probability */
if (new_E < E) {
if (new_E < best_E) {
copy_state(new_x, best_x, element_size, copyfunc);
best_E = new_E;
}
/* yay! take a step */
copy_state(new_x, x, element_size, copyfunc);
E = new_E;
++n_eless;
} else if (gsl_rng_uniform(r) < boltzmann(E, new_E, T, ¶ms)) {
/* yay! take a step */
copy_state(new_x, x, element_size, copyfunc);
E = new_E;
++n_accepts;
} else {
++n_rejects;
}
}
if (print_position) {
/* see if we need to print stuff as we go */
/* printf("%5d %12g %5d %3d %3d %3d", n_iter, T, n_evals, */
/* 100*n_eless/n_steps, 100*n_accepts/n_steps, */
/* 100*n_rejects/n_steps); */
printf ("%5d %7d %12g", n_iter, n_evals, T);
print_position (x);
printf (" %12g %12g\n", E, best_E);
}
/* apply the cooling schedule to the temperature */
/* FIXME: I should also introduce a cooling schedule for the iters */
T *= T_factor;
++n_iter;
if (T < params.t_min) {
break;
}
}
/* at the end, copy the result onto the initial point, so we pass it
back to the caller */
copy_state(best_x, x0_p, element_size, copyfunc);
if (copyfunc) {
destructor(x);
destructor(new_x);
destructor(best_x);
} else {
free (x);
free (new_x);
free (best_x);
}
}
/* implementation of a simulated annealing algorithm with many tries */
/// 模擬退火方法的一個變種
void
gsl_siman_solve_many (const gsl_rng * r, void *x0_p, gsl_siman_Efunc_t Ef,
gsl_siman_step_t take_step,
gsl_siman_metric_t distance,
gsl_siman_print_t print_position,
size_t element_size,
gsl_siman_params_t params)
{
/* the new set of trial points, and their energies and probabilities */
void *x, *new_x;
double *energies, *probs, *sum_probs;
double Ex; /* energy of the chosen point */
double T, T_factor; /* the temperature and a step multiplier */
int i;
double u; /* throw the die to choose a new "x" */
int n_iter;
if (print_position) {
printf ("#-iter temperature position");
printf (" delta_pos energy\n");
}
x = (void *) malloc (params.n_tries * element_size);
new_x = (void *) malloc (params.n_tries * element_size);
energies = (double *) malloc (params.n_tries * sizeof (double));
probs = (double *) malloc (params.n_tries * sizeof (double));
sum_probs = (double *) malloc (params.n_tries * sizeof (double));
T = params.t_initial;
T_factor = 1.0 / params.mu_t;
memcpy (x, x0_p, element_size);
n_iter = 0;
while (1)
{
Ex = Ef (x);
for (i = 0; i < params.n_tries - 1; ++i)
{ /* only go to N_TRIES-2 */
/* center the new_x[] around x, then pass it to take_step() */
sum_probs[i] = 0;
memcpy ((char *)new_x + i * element_size, x, element_size);
take_step (r, (char *)new_x + i * element_size, params.step_size);
energies[i] = Ef ((char *)new_x + i * element_size);
probs[i] = boltzmann(Ex, energies[i], T, ¶ms);
}
/* now add in the old value of "x", so it is a contendor */
memcpy ((char *)new_x + (params.n_tries - 1) * element_size, x, element_size);
energies[params.n_tries - 1] = Ex;
probs[params.n_tries - 1] = boltzmann(Ex, energies[i], T, ¶ms);
/* now throw biased die to see which new_x[i] we choose */
sum_probs[0] = probs[0];
for (i = 1; i < params.n_tries; ++i)
{
sum_probs[i] = sum_probs[i - 1] + probs[i];
}
u = gsl_rng_uniform (r) * sum_probs[params.n_tries - 1];
for (i = 0; i < params.n_tries; ++i)
{
if (u < sum_probs[i])
{
memcpy (x, (char *) new_x + i * element_size, element_size);
break;
}
}
if (print_position)
{
printf ("%5d\t%12g\t", n_iter, T);
print_position (x);
printf ("\t%12g\t%12g\n", distance (x, x0_p), Ex);
}
T *= T_factor;
++n_iter;
if (T < params.t_min)
{
break;
}
}
/* now return the value via x0_p */
memcpy (x0_p, x, element_size);
/* printf("the result is: %g (E=%g)\n", x, Ex); */
free (x);
free (new_x);
free (energies);
free (probs);
free (sum_probs);
}
下面是一個測試程式碼，其中重寫了部分GSL中的mt19937隨機數生成器
#include <math.h>
#include <stdlib.h>
#include <string.h>
#include "siman.h"
#include <stdio.h>
/****************************************************
以下為隨機數生成器部分
*****************************************************/
#define N 624
#define M 397
static const unsigned long UPPER_MASK = 0x80000000UL;
static const unsigned long LOWER_MASK = 0x7fffffffUL;
typedef struct
{
unsigned long mt[N];
int mti;
}
mt_state_t;
void mt_set(void *vstate, unsigned long int s)
{
mt_state_t *state = (mt_state_t *)vstate;
int i;
if (s == 0) s = 4357;
state->mt[0] = s & 0xffffffffUL;
for (i = 1; i < N; i++)
{
state->mt[i] = (1812433253UL * (state->mt[i - 1] ^ (state->mt[i - 1] >> 30)) + i);
state->mt[i] &= 0xffffffffUL;
}
state->mti = i;
}
inline unsigned long mt_get(void *vstate)
{
mt_state_t *state = (mt_state_t *)vstate;
unsigned long k;
unsigned long int *const mt = state->mt;
#define MAGIC(y) (((y)&0x1) ? 0x9908b0dfUL : 0)
if (state->mti >= N)
{
int kk;
for (kk = 0; kk < N - M; kk++)
{
unsigned long y = (mt[kk] & UPPER_MASK) | (mt[kk + 1] & LOWER_MASK);
mt[kk] = mt[kk + M] ^ (y >> 1) ^ MAGIC(y);
}
for (; kk < N - 1; kk++)
{
unsigned long y = (mt[kk] & UPPER_MASK) | (mt[kk + 1] & LOWER_MASK);
mt[kk] = mt[kk + (M - N)] ^ (y >> 1) ^ MAGIC(y);
}
unsigned long y = (mt[N - 1] & UPPER_MASK) | (mt[0] & LOWER_MASK);
mt[N - 1] = mt[M - 1] ^ (y >> 1) ^ MAGIC(y);
state->mti = 0;
}
k = mt[state->mti];
k ^= (k >> 11);
k ^= (k << 7) & 0x9d2c5680UL;
k ^= (k << 15) & 0xefc60000UL;
k ^= (k >> 18);
state->mti++;
return k;
}
double mt_get_double(void *vstate)
{
return mt_get(vstate) / 4294967296.0;
}
const gsl_rng_type mt_type =
{
"mt19937",
0xffffffffUL,
0,
sizeof(mt_state_t),
&mt_set,
&mt_get,
&mt_get_double
};
/******************************************************
以下為模擬退火演算法的測試部分
*******************************************************/
#define N_TRIES 500
#define ITERS_FIXED_T 2000
#define STEP_SIZE 0.02
#define K 1.0
#define T_INITIAL 1000
#define MU_T 1.006
#define T_MIN 1.0e-3
gsl_siman_params_t params
= {N_TRIES, ITERS_FIXED_T, STEP_SIZE,
K, T_INITIAL, MU_T, T_MIN};
double E1(void *xp)
{
double x = *((double *)xp);
return exp(-pow((x-1.0),2.0))*sin(8*x);
}
double E2(void *xp)
{
double x = *((double *)xp);
double y = *((double *)xp+1);
double z = 3*pow(1-x,2)*exp(-pow(x,2)-pow(y+1,2))
-10*(x/5-pow(x,3)-pow(y,5))*exp(-pow(x,2)-pow(y,2))
-exp(-pow(x+1,2)-pow(y,2))/3;
return -z;
}
double M1(void *xp, void *yp)
{
double x = *((double *)xp);
double y = *((double *)yp);
return fabs(x-y);
}
double M2(void *xp, void *yp)
{
double x0 = *((double *)xp);
double y0 = *((double *)xp+1);
double x1 = *((double *)yp);
double y1 = *((double *)yp+1);
return sqrt(pow(x1-x0,2)+pow(y1-y0,2));
}
void S1(const gsl_rng *r, void *xp, double step_size)
{
double old_x = *((double *)xp);
double new_x;
double u = gsl_rng_uniform(r);
new_x = u*2*step_size - step_size + old_x;
memcpy(xp, &new_x, sizeof(new_x));
}
void S2(const gsl_rng *r, void *xp, double step_size)
{
double *old_x = (double *)xp;
double *new_x = new double[2];
double u = gsl_rng_uniform(r);
double v = gsl_rng_uniform(r);
new_x[0] = u*2*step_size - step_size + old_x[0];
new_x[1] = v*2*step_size - step_size + old_x[1];
memcpy(xp, new_x, sizeof(double)*2);
}
void P1(void *xp)
{
printf("%12g", *((double *)xp));
}
void P2(void *xp)
{
printf("%12g,%12g", *((double *)xp), *((double *)xp+1));
}
int main(int argc, char *argv[])
{
const gsl_rng_type *T;
gsl_rng *r;
double x_initial1 = 15.5;
double x_initial2[2] = {1.5, 1.5};
gsl_rng_env_setup(&mt_type, 10);
T = gsl_rng_default;
r = gsl_rng_alloc(T);
/// 最小值應為 xmin = 1.36313, ymin = -0.872871
gsl_siman_solve(r, &x_initial1, E1, S1, M1, P1,
NULL, NULL, NULL,
sizeof(double), params);
getchar();
/// 最小值應為 xmin = (-0.0093, 1.5814), ymin = -8.1062
gsl_siman_solve_many(r, x_initial2, E2, S2, M2, P2,
sizeof(double)*2, params);
gsl_rng_free(r);
system("pause");
return 0;
}
下面是模擬退火演算法的經典問題，Travelling Salesman Problem的原始碼，此處改寫了隨機數生成器，採用c標準庫函式srand和rand來產生隨機數。
#include <math.h>
#include <string.h>
#include <stdio.h>
#include "siman.h"
#define M_PI 3.14159265358979323846264338328
void set(void *vstate, unsigned long int s)
{
srand(s);
}
__inline unsigned long get(void *vstate)
{
return rand();
}
double get_double(void *vstate)
{
return get(vstate) / ((double)RAND_MAX);
}
const gsl_rng_type mt_type =
{
"default",
RAND_MAX,
0,
0,
&set,
&get,
&get_double
};
/* set up parameters for this simulated annealing run */
#define N_TRIES 200 /* how many points do we try before stepping */
#define ITERS_FIXED_T 2000 /* how many iterations for each T? */
#define STEP_SIZE 1.0 /* max step size in random walk */
#define K 1.0 /* Boltzmann constant */
#define T_INITIAL 5000.0 /* initial temperature */
#define MU_T 1.002 /* damping factor for temperature */
#define T_MIN 5.0e-1
gsl_siman_params_t params = {N_TRIES, ITERS_FIXED_T, STEP_SIZE,
K, T_INITIAL, MU_T, T_MIN};
struct s_tsp_city {
const char * name;
double lat, longitude; /* coordinates */
};
typedef struct s_tsp_city Stsp_city;
void prepare_distance_matrix(void);
void exhaustive_search(void);
void print_distance_matrix(void);
double city_distance(Stsp_city c1, Stsp_city c2);
double Etsp(void *xp);
double Mtsp(void *xp, void *yp);
void Stsp(const gsl_rng * r, void *xp, double step_size);
void Ptsp(void *xp);
/* in this table, latitude and longitude are obtained from the US
Census Bureau, at http://www.census.gov/cgi-bin/gazetteer */
Stsp_city cities[] = {{"Santa Fe", 35.68, 105.95},
{"Phoenix", 33.54, 112.07},
{"Albuquerque", 35.12, 106.62},
{"Clovis", 34.41, 103.20},
{"Durango", 37.29, 107.87},
{"Dallas", 32.79, 96.77},
{"Tesuque", 35.77, 105.92},
{"Grants", 35.15, 107.84},
{"Los Alamos", 35.89, 106.28},
{"Las Cruces", 32.34, 106.76},
{"Cortez", 37.35, 108.58},
{"Gallup", 35.52, 108.74}};
#define N_CITIES (sizeof(cities)/sizeof(Stsp_city))
double distance_matrix[N_CITIES][N_CITIES];
/* distance between two cities */
double city_distance(Stsp_city c1, Stsp_city c2)
{
const double earth_radius = 6375.000; /* 6000KM approximately */
/* sin and cos of lat and long; must convert to radians */
double sla1 = sin(c1.lat*M_PI/180), cla1 = cos(c1.lat*M_PI/180),
slo1 = sin(c1.longitude*M_PI/180), clo1 = cos(c1.longitude*M_PI/180);
double sla2 = sin(c2.lat*M_PI/180), cla2 = cos(c2.lat*M_PI/180),
slo2 = sin(c2.longitude*M_PI/180), clo2 = cos(c2.longitude*M_PI/180);
double x1 = cla1*clo1;
double x2 = cla2*clo2;
double y1 = cla1*slo1;
double y2 = cla2*slo2;
double z1 = sla1;
double z2 = sla2;
double dot_product = x1*x2 + y1*y2 + z1*z2;
double angle = acos(dot_product);
/* distance is the angle (in radians) times the earth radius */
return angle*earth_radius;
}
/* energy for the travelling salesman problem */
double Etsp(void *xp)
{
/* an array of N_CITIES integers describing the order */
int *route = (int *) xp;
double E = 0;
unsigned int i;
for (i = 0; i < N_CITIES; ++i) {
/* use the distance_matrix to optimize this calculation; it had
better be allocated!! */
E += distance_matrix[route[i]][route[(i + 1) % N_CITIES]];
}
return E;
}
double Mtsp(void *xp, void *yp)
{
int *route1 = (int *) xp, *route2 = (int *) yp;
double distance = 0;
unsigned int i;
for (i = 0; i < N_CITIES; ++i) {
distance += ((route1[i] == route2[i]) ? 0 : 1);
}
return distance;
}
/* take a step through the TSP space */
void Stsp(const gsl_rng * r, void *xp, double step_size)
{
int x1, x2, dummy;
int *route = (int *) xp;
step_size = 0 ; /* prevent warnings about unused parameter */
/* pick the two cities to swap in the matrix; we leave the first
city fixed */
x1 = (gsl_rng_get (r) % (N_CITIES-1)) + 1;
do {
x2 = (gsl_rng_get (r) % (N_CITIES-1)) + 1;
} while (x2 == x1);
dummy = route[x1];
route[x1] = route[x2];
route[x2] = dummy;
}
void Ptsp(void *xp)
{
unsigned int i;
int *route = (int *) xp;
printf(" [");
for (i = 0; i < N_CITIES; ++i) {
printf(" %d ", route[i]);
}
printf("] ");
}
int main(void)
{
int x_initial[N_CITIES];
unsigned int i;
const gsl_rng * r = gsl_rng_alloc (gsl_rng_env_setup(&mt_type, 0)) ;
// gsl_ieee_env_setup ();
prepare_distance_matrix();
/* set up a trivial initial route */
printf("# initial order of cities:\n");
for (i = 0; i < N_CITIES; ++i) {
printf("# \"%s\"\n", cities[i].name);
x_initial[i] = i;
}
printf("# distance matrix is:\n");
print_distance_matrix();
printf("# initial coordinates of cities (longitude and latitude)\n");
/* this can be plotted with */
/* ./siman_tsp > hhh ; grep city_coord hhh | awk '{print $2 " " $3}' | xyplot -ps -d "xy" > c.eps */
for (i = 0; i < N_CITIES+1; ++i) {
printf("###initial_city_coord: %g %g \"%s\"\n",
-cities[x_initial[i % N_CITIES]].longitude,
cities[x_initial[i % N_CITIES]].lat,
cities[x_initial[i % N_CITIES]].name);
}
/* exhaustive_search(); */
gsl_siman_solve(r, x_initial, Etsp, Stsp, Mtsp, Ptsp, NULL, NULL, NULL,
N_CITIES*sizeof(int), params);
printf("# final order of cities:\n");
for (i = 0; i < N_CITIES; ++i) {
printf("# \"%s\"\n", cities[x_initial[i]].name);
}
printf("# final coordinates of cities (longitude and latitude)\n");
/* this can be plotted with */
/* ./siman_tsp > hhh ; grep city_coord hhh | awk '{print $2 " " $3}' | xyplot -ps -d "xy" > c.eps */
for (i = 0; i < N_CITIES+1; ++i) {
printf("###final_city_coord: %g %g %s\n",
-cities[x_initial[i % N_CITIES]].longitude,
cities[x_initial[i % N_CITIES]].lat,
cities[x_initial[i % N_CITIES]].name);
}
printf("# ");
fflush(stdout);
#if 0
system("date");
#endif /* 0 */
fflush(stdout);
system("pause");
return 0;
}
void prepare_distance_matrix()
{
unsigned int i, j;
double dist;
for (i = 0; i < N_CITIES; ++i) {
for (j = 0; j < N_CITIES; ++j) {
if (i == j) {
dist = 0;
} else {
dist = city_distance(cities[i], cities[j]);
}
distance_matrix[i][j] = dist;
}
}
}
void print_distance_matrix()
{
unsigned int i, j;
for (i = 0; i < N_CITIES; ++i) {
printf("# ");
for (j = 0; j < N_CITIES; ++j) {
printf("%15.8f ", distance_matrix[i][j]);
}
printf("\n");
}
}
/* [only works for 12] search the entire space for solutions */
static double best_E = 1.0e100, second_E = 1.0e100, third_E = 1.0e100;
static int best_route[N_CITIES];
static int second_route[N_CITIES];
static int third_route[N_CITIES];
static void do_all_perms(int *route, int n);
void exhaustive_search()
{
static int initial_route[N_CITIES] = {0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11};
printf("\n# ");
fflush(stdout);
#if 0
system("date");
#endif
fflush(stdout);
do_all_perms(initial_route, 1);
printf("\n# ");
fflush(stdout);
#if 0
system("date");
#endif /* 0 */
fflush(stdout);
printf("# exhaustive best route: ");
Ptsp(best_route);
printf("\n# its energy is: %g\n", best_E);
printf("# exhaustive second_best route: ");
Ptsp(second_route);
printf("\n# its energy is: %g\n", second_E);
printf("# exhaustive third_best route: ");
Ptsp(third_route);
printf("\n# its energy is: %g\n", third_E);
}
/* James Theiler's recursive algorithm for generating all routes */
static void do_all_perms(int *route, int n)
{
if (n == (N_CITIES-1)) {
/* do it! calculate the energy/cost for that route */
double E;
E = Etsp(route); /* TSP energy function */
/* now save the best 3 energies and routes */
if (E < best_E) {
third_E = second_E;
memcpy(third_route, second_route, N_CITIES*sizeof(*route));
second_E = best_E;
memcpy(second_route, best_route, N_CITIES*sizeof(*route));
best_E = E;
memcpy(best_route, route, N_CITIES*sizeof(*route));
} else if (E < second_E) {
third_E = second_E;
memcpy(third_route, second_route, N_CITIES*sizeof(*route));
second_E = E;
memcpy(second_route, route, N_CITIES*sizeof(*route));
} else if (E < third_E) {
third_E = E;
memcpy(route, third_route, N_CITIES*sizeof(*route));
}
} else {
int new_route[N_CITIES];
unsigned int j;
int swap_tmp;
memcpy(new_route, route, N_CITIES*sizeof(*route));
for (j = n; j < N_CITIES; ++j) {
swap_tmp = new_route[j];
new_route[j] = new_route[n];
new_route[n] = swap_tmp;
do_all_perms(new_route, n+1);
}
}
}
完整的VS2013工程可從如下網址下載：