You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.

109 lines
3.3 KiB

2 months ago
  1. /*
  2. * Author: Heinrich Schuchardt <xypron.glpk@gmx.de>
  3. *
  4. * This model solves a clustering problem:
  5. *
  6. * Out of 50 towns select 3 to be cluster centers and assign the other
  7. * towns to the clusters such that the sum of the population weighted
  8. * euclidian distances between towns and centers is minimized.
  9. *
  10. * The solution is saved as a scalable vector graphic file with a
  11. * pseudo-random file name.
  12. */
  13. # Output file
  14. param fn, symbolic := "00000" & 100000 * Uniform01();
  15. param f, symbolic := "ct" & substr(fn, length(fn) - 4) & ".svg";
  16. # Centers
  17. param nc := 3;
  18. set C := {1 .. nc};
  19. # Towns
  20. param nt := 50;
  21. set T := {1 .. nt};
  22. param xt{T} := Uniform01();
  23. param yt{T} := Uniform01();
  24. param pt{T} := ceil(1000 * Uniform01());
  25. # Image size
  26. param scale := 1000;
  27. # Colors
  28. # saturation [0, 255]
  29. param sat := 192;
  30. param hue{c in C} := 6 * (c - 1) / nc;
  31. param red{c in C} :=
  32. if hue[c] <= 1 or hue[c] >= 5 then 255
  33. else (if hue[c] >=2 and hue[c] <= 4 then 255 - sat
  34. else (if hue[c] <=2 then 255 - sat + sat * (2-hue[c])
  35. else 255 - sat + sat * (hue[c]-4) ));
  36. param green{c in C} :=
  37. if hue[c] >= 1 and hue[c] <= 3 then 255
  38. else (if hue[c] >= 4 then 255 - sat
  39. else (if hue[c] <=1 then 255 - sat + sat * hue[c]
  40. else 255 - sat + sat * (4-hue[c]) ));
  41. param blue{c in C} :=
  42. if hue[c] >= 3 and hue[c] <= 5 then 255
  43. else (if hue[c] <=2 then 255 - sat
  44. else (if hue[c] <=3 then 255 - sat + sat * (hue[c]-2)
  45. else 255 - sat + sat * (6-hue[c]) ));
  46. var x{T};
  47. var y{T,T}, binary;
  48. minimize obj : sum{c in T, t in T} y[c,t] * pt[t]
  49. * sqrt((xt[c] - xt[t])^2 + (yt[c] - yt[t])^2);
  50. s.t. sumx : sum{c in T} x[c] = nc;
  51. s.t. cxy{c in T, t in T} : y[c,t] <= x[c];
  52. s.t. sumy{t in T} : sum{c in T} y[c,t] = 1;
  53. solve;
  54. for {c in T : x[c] > .5} {
  55. printf "Center %5.4f %5.4f\n", xt[c], yt[c];
  56. for {t in T : y[c,t] > .5} {
  57. printf " Town %5.4f %5.4f (%5.0f)\n", xt[t], yt[t], pt[t];
  58. }
  59. }
  60. # Output the solution as scalable vector graphic
  61. # header
  62. printf "<?xml version=""1.0"" standalone=""no""?>\n" > f;
  63. printf "<!DOCTYPE svg PUBLIC ""-//W3C//DTD SVG 1.1//EN"" \n" >> f;
  64. printf """http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd"">\n" >> f;
  65. printf "<svg width=""%d"" height=""%d"" version=""1.0"" \n",
  66. 1.2 * scale, 1.2 * scale >> f;
  67. printf "xmlns=""http://www.w3.org/2000/svg"">\n" >> f;
  68. # background
  69. printf "<rect x=""0"" y=""0"" width=""%d"" height=""%d""" &
  70. " stroke=""none"" fill=""white""/>\n",
  71. 1.2 * scale, 1.2 * scale>> f;
  72. # border
  73. printf "<rect x=""%d"" y=""%d"" width=""%d"" height=""%d""" &
  74. " stroke=""black"" stroke-width="".5"" fill=""white""/>\n",
  75. .1 * scale, .1 * scale, scale, scale >> f;
  76. # circles for towns
  77. for {t in T}
  78. printf {s in T, c in C : y[s,t] > .5
  79. && c = floor( .5 + sum{u in T : u <= s} x[u])}
  80. "<circle cx=""%f"" cy=""%f"" r=""%f"" stroke=""black"" " &
  81. "stroke-width=""1"" fill=""rgb(%d,%d,%d)""/>\n",
  82. (.1 + xt[t]) * scale, (.1 + yt[t]) * scale, .001 * sqrt(pt[t]) * scale,
  83. red[c], green[c] , blue[c] >> f;
  84. # lines from towns to assigned centers
  85. for {t in T, c in T : y[c,t] > .5}
  86. printf "<line x1=""%f"" y1=""%f"" x2=""%f"" y2=""%f""" &
  87. " style=""stroke:black;stroke-width:.5""/>\n",
  88. (.1 + xt[c]) * scale, (.1 + yt[c]) * scale,
  89. (.1 + xt[t]) * scale, (.1 + yt[t]) * scale >> f;
  90. printf "</svg>\n" >> f;
  91. end;