home *** CD-ROM | disk | FTP | other *** search
/ The C Users' Group Library 1994 August / wc-cdrom-cusersgrouplibrary-1994-08.iso / listings / v_10_07 / 1007075a < prev    next >
Text File  |  1992-05-13  |  6KB  |  236 lines
  1. /*
  2.   This is a simulated annealing solution to the
  3.   quadratic assignment problem (a.k.a. placement of
  4.   facilities).  The particular data sets were drawn from
  5.   [1] and solutions found were optimal (where optimal
  6.   solutions are known (5*5, 6*6, 7*7, 8*8, and 12*12
  7.   datasets) or better than or equal to previous known
  8.   results for the larger datasets (15*15, 20*20 and
  9.   30*30).  Time taken was surprisingly short.  Found a
  10.   better solution for 30*30 on second run, which lasted
  11.   only a few minutes on an 8MHz PC (compiled under
  12.   Turbo C++).
  13.  
  14.   [1] Nugent, C.E., Vollmann, T.E., and Ruml, J. --
  15.       "An Experimental Comparison of Techniques for
  16.       the Assignment of Facilities to Locations",
  17.       _Operations Research_ 16, pp.  150-173 January,
  18.       1968
  19.  
  20.   [2] Kirkpatrick, S., Gelatt, C.D. Jr., and Vecchi,
  21.       M.P. -- "Optimization by Simluated Annealing",
  22.       _Science_ 220:4598, pp. 671-680, May 1983.
  23. */
  24.  
  25. #include <stdlib.h>
  26. #include <stdio.h>
  27. #include <limits.h>
  28. #include <math.h>
  29. #include <time.h>
  30.  
  31. #ifdef _cplusplus
  32. #define INLINE inline
  33. #else
  34. #define INLINE
  35. #endif
  36.  
  37. #define MAX 30
  38.  
  39. /* prototypes */
  40. void init(void);
  41. int energy(void);
  42. int acceptable(int E, int E2);
  43. int metastable(int accepts, int rejects);
  44. int annealing(int E);
  45.  
  46. int n;                            /* data array size */
  47. double init_T;              /* initial "temperature" */
  48. double alpha;           /* cooling schedule variable */
  49. int accept_n;                /* metastable indicator */
  50. int reject_n;                /* metastable indicator */
  51. int freeze_n;                    /* system is frozen */
  52. double T;                   /* Current "temperature" */
  53.  
  54.  
  55. /* number of units to ship between plant i and plant j
  56.    is units[i][j] == units[j][i]
  57. */
  58. int units[MAX][MAX];
  59.  
  60. /* cost/unit shipped between site i and site j is
  61.    in cost[i][j] == cost[j][i]
  62. */
  63. int cost[MAX][MAX];
  64.  
  65. /* site[i] == p means plant p is at site i */
  66. int site[MAX];
  67.  
  68.  
  69. void init(void)
  70. {
  71.   int i, j, input_array[MAX][MAX];
  72.   char file[13];
  73.  
  74.   FILE *fp;
  75.  
  76.   scanf("%d", &n);
  77.   scanf("%12s", file);
  78.   scanf("%lf", &init_T);
  79.   scanf("%lf", &alpha);
  80.   scanf("%d", &accept_n);
  81.   scanf("%d", &reject_n);
  82.   scanf("%d", &freeze_n);
  83.  
  84.   printf("Matrix size           : %d\n", n);
  85.   printf("Data file             : %s\n", file);
  86.   printf("Initial temperature   : %lf\n", init_T);
  87.   printf("alpha (T = alpha * T) : %lf\n", alpha);
  88.   printf("accept max            : %d\n", accept_n);
  89.   printf("reject max            : %d\n", reject_n);
  90.   printf("freeze max            : %d\n", freeze_n);
  91.   printf("\n\n");
  92.  
  93.   fp = fopen(file, "rt");
  94.   if (!fp)
  95.   {
  96.     printf("cannot open file: %s\n", file);
  97.     exit(4);
  98.   }
  99.  
  100.   /* Get cost/unit array */
  101.   for (i = 0; i < n; i++)
  102.     for (j = 0; j < n; j++)
  103.       fscanf(fp, "%d", &input_array[i][j]);
  104.   fclose(fp);
  105.  
  106.  /* Unit deliveries are in upper triangular section
  107.     of the input array.  Transfer costs are in lower
  108.     triangular section of the input array.  Copy
  109.     each part into its own array, and reflect it
  110.     into the other triangular section so that
  111.     units[i][j] == units[j][i] and
  112.     cost[i][j] == cost[j][i].
  113.  */
  114.   for (i = 0; i < n; i++)
  115.     for (j = i + 1; j < n; j++)
  116.     {
  117.       units[i][j] = units[j][i] = input_array[i][j];
  118.       cost[i][j] = cost[j][i] = input_array[j][i];
  119.     }
  120.  
  121.   /* Generate initial plant-site layout */
  122.   for (i = 0; i < n; i++)
  123.     site[i] = i;
  124. }
  125.  
  126.  
  127. INLINE int energy(void)
  128. {
  129.   int plant_i, plant_j, E = 0;
  130.  
  131.   for (plant_i = 0; plant_i < n; plant_i++)
  132.     for (plant_j = plant_i + 1; plant_j < n; plant_j++)
  133.       E += units[plant_i][plant_j] *
  134.            cost[site[plant_i]][site[plant_j]];
  135.   return E;
  136. }
  137.  
  138.  
  139. INLINE int acceptable(int E, int E2)
  140. {
  141.   static double rand_max = (double) RAND_MAX;
  142.   double random_number;
  143.   double P;
  144.  
  145.   if (E2 <= E)
  146.     return 1;
  147.  
  148.   random_number = random(RAND_MAX) / rand_max;
  149.   P = exp((double) (E - E2) / T);
  150.  
  151.   return random_number < P;
  152. }
  153.  
  154.  
  155. INLINE int metastable(int accepts, int rejects)
  156. {
  157.   return accepts > accept_n || rejects > reject_n;
  158. }
  159.  
  160.  
  161. int annealing(int E)
  162. {
  163.   int accepts, rejects;
  164.   int swap1, swap2, temp;
  165.   int E2;
  166.  
  167.   accepts = 0;
  168.   rejects = 0;
  169.   do                              /* metastable loop */
  170.   {
  171.     do          /* randomly pick two distinct plants */
  172.     {
  173.  
  174.       swap1 = rand() % n;
  175.       swap2 = rand() % n;
  176.     } while (swap1 == swap2);
  177.     temp = site[swap1];          /* swap plant sites */
  178.     site[swap1] = site[swap2];
  179.     site[swap2] = temp;
  180.     E2 = energy();         /* energy of new solution */
  181.  
  182.     if (acceptable(E, E2))        /* Metropolis test */
  183.     {
  184.       accepts++;
  185.       E = E2;
  186.     }
  187.     else
  188.     {                          /* no good, throw out */
  189.       rejects++;
  190.       temp = site[swap1];          /* swap them back */
  191.       site[swap1] = site[swap2];
  192.       site[swap2] = temp;
  193.     }
  194.   } while (!metastable(accepts, rejects));
  195.  
  196.   return E;   /* energy of current, metastable state */
  197. }
  198.  
  199.  
  200. main()
  201. {
  202.   int i;
  203.   int E, old_E;
  204.   int freeze;
  205.  
  206.   randomize();
  207.   init();
  208.  
  209.   T = init_T;
  210.   old_E = E = energy();     /* E of inital solution */
  211.   freeze = 0;
  212.  
  213.   do                               /* freezing loop */
  214.   {
  215.     E = annealing(E);     /* anneal at current temp */
  216.  
  217.     if (E != old_E)             /* check for change */
  218.     {
  219.       freeze = 0;                 /* still changing */
  220.       old_E = E;
  221.     }
  222.     else
  223.       freeze++;           /* closing in on freezing */
  224.  
  225.     T *= alpha;                /* lower temperature */
  226.  
  227.   } while (freeze < freeze_n);
  228.  
  229.   /* display solution */
  230.   printf("Frozen at T = %f E = %d\n", T, E);
  231.   for (i = 0; i < n; i++)
  232.     printf("place plant %d at site %d\n", i, site[i]);
  233.  
  234.   return 0;
  235. }
  236.