Multiple Solutions in Minimum Spanning Tree example

In [1], I discussed some LP and MIP formulations for the Minimum Spanning Tree (MST) problem. 

Minimum Spanning Tree visualized through Google Maps

Here, I focus on two formulations: a multicommodity network approach (this can be solved as a large LP) and a MIP formulation based on techniques we know from the Traveling Salesman Problem (TSP). The main issue I want to discuss is the presence of multiple optimal solutions.

Data set

The data looks like:

set   i cities /                                  c1   'Manchester, N.H.'           c2   'Montpelier, Vt.'            c3   'Detroit, Mich.'             c4   'Cleveland, Ohio'            c5   'Charleston, W.Va.'          c6   'Louisville, Ky.'            c7   'Indianapolis, Ind.'         c8   'Chicago, Ill.'              c9   'Milwaukee, Wis.'            c10  'Minneapolis, Minn.'         c11  'Pierre, S.D.'               c12  'Bismarck, N.D.'             c13  'Helena, Mont.'              c14  'Seattle, Wash.'             c15  'Portland, Ore.'             c16  'Boise, Idaho'               c17  'Salt Lake City, Utah'       c18  'Carson City, Nevada'        c19  'Los Angeles, Calif.'        c20  'Phoenix, Ariz.'             c21  'Santa Fe, N.M.'             c22  'Denver, Colo.'              c23  'Cheyenne, Wyo.'             c24  'Omaha, Neb.'                c25  'Des Moines, Iowa'           c26  'Kansas City, Mo.'           c27  'Topeka, Kans.'              c28  'Oklahoma City, Okla.'       c29  'Dallas, Tex.'               c30  'Little Rock, Ark.'          c31  'Memphis, Tenn.'             c32  'Jackson, Miss.'             c33  'New Orleans, La.'           c34  'Birmingham, Ala.'           c35  'Atlanta, Ga.'               c36  'Jacksonville, Fla.'         c37  'Columbia, S.C.'             c38  'Raleigh, N.C.'              c39  'Richmond, Va.'              c40  'Washington, D.C.'           c41  'Boston, Mass.'              c42  'Portland, Me.'         /;   alias (i,j,k);   table d(i,j) 'from dantzig42 data set in TSPLIB'      c1  c2  c3  c4  c5  c6  c7  c8  c9 c10 c11 c12 c13 c14 c15 c16 c17 c18 c19 c20 c21 c2    8 c3   39  45 c4   37  47   9 c5   50  49  21  15 c6   61  62  21  20  17 c7   58  60  16  17  18   6 c8   59  60  15  20  26  17  10 c9   62  66  20  25  31  22  15   5 c10  81  81  40  44  50  41  35  24  20 c11 103 107  62  67  72  63  57  46  41  23 c12 108 117  66  71  77  68  61  51  46  26  11 c13 145 149 104 108 114 106  99  88  84  63  49  40 c14 181 185 140 144 150 142 135 124 120  99  85  76  35 c15 187 191 146 150 156 142 137 130 125 105  90  81  41  10 c16 161 170 120 124 130 115 110 104 105  90  72  62  34  31  27 c17 142 146 101 104 111  97  91  85  86  75  51  59  29  53  48  21 c18 174 178 133 138 143 129 123 117 118 107  83  84  54  46  35  26  31 c19 185 186 142 143 140 130 126 124 128 118  93 101  72  69  58  58  43  26 c20 164 165 120 123 124 106 106 105 110 104  86  97  71  93  82  62  42  45  22 c21 137 139  94  96  94  80  78  77  84  77  56  64  65  90  87  58  36  68  50  30 c22 117 122  77  80  83  68  62  60  61  50  34  42  49  82  77  60  30  62  70  49  21 c23 114 118  73  78  84  69  63  57  59  48  28  36  43  77  72  45  27  59  69  55  27 c24  85  89  44  48  53  41  34  28  29  22  23  35  69 105 102  74  56  88  99  81  54 c25  77  80  36  40  46  34  27  19  21  14  29  40  77 114 111  84  64  96 107  87  60 c26  87  89  44  46  46  30  28  29  32  27  36  47  78 116 112  84  66  98  95  75  47 c27  91  93  48  50  48  34  32  33  36  30  34  45  77 115 110  83  63  97  91  72  44 c28 105 106  62  63  64  47  46  49  54  48  46  59  85 119 115  88  66  98  79  59  31 c29 111 113  69  71  66  51  53  56  61  57  59  71  96 130 126  98  75  98  85  62  38 c30  91  92  50  51  46  30  34  38  43  49  60  71 103 141 136 109  90 115  99  81  53 c31  83  85  42  43  38  22  26  32  36  51  63  75 106 142 140 112  93 126 108  88  60 c32  89  91  55  55  50  34  39  44  49  63  76  87 120 155 150 123 100 123 109  86  62 c33  95  97  64  63  56  42  49  56  60  75  86  97 126 160 155 128 104 128 113  90  67 c34  74  81  44  43  35  23  30  39  44  62  78  89 121 159 155 127 108 136 124 101  75 c35  67  69  42  41  31  25  32  41  46  64  83  90 130 164 160 133 114 146 134 111  85 c36  74  76  61  60  42  44  51  60  66  83 102 110 147 185 179 155 133 159 146 122  98 c37  57  59  46  41  25  30  36  47  52  71  93  98 136 172 172 148 126 158 147 124 121 c38  45  46  41  34  20  34  38  48  53  73  96  99 137 176 178 151 131 163 159 135 108 c39  35  37  35  26  18  34  36  46  51  70  93  97 134 171 176 151 129 161 163 139 118 c40  29  33  30  21  18  35  33  40  45  65  87  91 117 166 171 144 125 157 156 139 113 c41   3  11  41  37  47  57  55  58  63  83 105 109 147 186 188 164 144 176 182 161 134 c42   5  12  55  41  53  64  61  61  66  84 111 113 150 186 192 166 147 180 188 167 140   +     c22 c23 c24 c25 c26 c27 c28 c29 c30 c31 c32 c33 c34 c35 c36 c37 c38 c39 c40 c41 c23   5 c24  32  29 c25  40  37   8 c26  36  39  12  11 c27  32  36   9  15   3 c28  36  42  28  33  21  20 c29  47  53  39  42  29  30  12 c30  61  62  36  34  24  28  20  20 c31  64  66  39  36  27  31  28  28   8 c32  71  78  52  49  39  44  35  24  15  12 c33  76  82  62  59  49  53  40  29  25  23  11 c34  79  81  54  50  42  46  43  39  23  14  14  21 c35  84  86  59  52  47  51  53  49  32  24  24  30   9 c36 105 107  79  71  66  70  70  60  48  40  36  33  25  18 c37  97  99  71  65  59  63  67  62  46  38  37  43  23  13  17 c38 102 103  73  67  64  69  75  72  54  46  49  54  34  24  29  12 c39 102 101  71  65  65  70  84  78  58  50  56  62  41  32  38  21   9 c40  95  97  67  60  62  67  79  82  62  53  59  66  45  38  45  27  15   6 c41 119 116  86  78  84  88 101 108  88  80  86  92  71  64  71  54  41  32  25 c42 124 119  90  87  90  94 107 114  77  86  92  98  80  74  77  60  48  38  32   6 ;

 

Notice that the distance matrix contains integers, and that there are many duplicates.

Multiple optimal solutions in LPs

The (possible) presence of multiple solutions in an LP can be detected by looking at the marginals (duals). If there are marginals with EPS values, we have a model with a dual degenerate solution [2]. The EPS value indicates both: (1) the dual is zero and (2) the row (or variable) is non-basic. When we use the multicommodity flow formulation, solutions are automatically integer-valued, so we can use an LP solver. If we look at the solution, in particular the duals of constraint nodebal, we see:

We see a lot of EPS values. This indicates that there are multiple solutions. It is noted that some of these may just be different dual solutions. To investigate further, we use a MIP model to enumerate all optimal integer solutions.

Enumerating integer solutions in a MIP model

The previous LP model is not really suited for enumerating integer solutions. It is very big. Its main advantage is that it can be solved as an LP.  As we need explicit binary variables anyway for our enumeration scheme, it is better to use a different, smaller MIP model. From [1], we use the lifted MTZ model.

There are three approaches we can try to enumerate all optimal integer solutions:

• Standard no-good constraints that forbid a precious found solution
• A simpler no-good constraint that exploits some model structure
• The solution pool in the MIP solver 

We focus on the binary variables that describe the tree: \[\color{darkred}x_{i,j} = \begin{cases}1 & \text{if arc $i \rightarrow j$ is part of the tree} \\  0 & \text{otherwise}\end{cases} \] It is noted that the models use a directed network (I am more familiar with that).

The standard no-good constraint to forbid a previously found solution \(\color{darkblue}x^*_{i,j}\in \{0,1\}\) can look like: \[\sum_{i,j|\color{darkblue}x^*_{i,j}=1}\color{darkred}x_{i,j} - \sum_{i,j|\color{darkblue}x^*_{i,j}=0}\color{darkred}x_{i,j} \le \sum_{i,j|\color{darkblue}x^*_{i,j}=1}1 - 1 \] A different way to say the same thing is: 

\[\sum_{i,j} \color{darkblue}x^*_{i,j}\cdot\color{darkred}x_{i,j} - \sum_{i,j} \left(1-\color{darkblue}x^*_{i,j}\right)\cdot \color{darkred}x_{i,j} \le \sum_{i,j}\color{darkblue}x^*_{i,j} - 1 \] This forbids exactly the solution \(\color{darkred}x_{i,j}=\color{darkblue}x^*_{i,j}\) but nothing else.

As we know that for any tree \[\sum_{i,j} \color{darkred}x_{i,j} = \color{darkblue}n-1\] where \( \color{darkblue}n\) is the number of nodes, we can use a much simpler cut: \[ \sum_{i,j} \color{darkblue}x^*_{i,j}\cdot\color{darkred}x_{i,j} \le \color{darkblue}n-2\] This can be easily derived by plugging  \[\sum_{i,j} \color{darkred}x_{i,j} = \sum_{i,j} \color{darkblue}x^*_{i,j} = \color{darkblue}n-1\] into our no-good cut.

The algorithm is to repeat solving and adding cuts until the objective starts to deteriorate:

1. Number of cuts = 0
2. Solve problem
3. If the number of cuts > 0 and the objective starts to deteriorate: STOP
4. Add cut to the problem to prevent the current solution from being repeated
5. Go to step 2

Basically, we keep on finding the next best solution as long as the objective remains the same. This approach shows:

----    396 PARAMETER info  

first objective                     591.000,    last objective                      592.000
number of optimal integer solutions  24.000

So we have 24 solutions with objective 591. We stop when we see the 25-th result with an objective of 592.

The Cplex solution pool finds the 24 optimal integer solutions very fast (2.8 seconds). We see the same problem as in [3]: if we use SolnPoolAGap=0 and Threads=1, Cplex is missing a solution, and it reports 23 optimal solutions.

Here are all the optimal solutions.

The plot was produced using simple HTML and SVG (Scalable Vector Graphics).

Theory + Practice

There is a theorem in Graph Theory: if the lengths in the data set for a minimum spanning tree problem are unique (no duplicates), then the solution will be unique. The proof (by contradiction) can even be found on Wikipedia [4].

To test this (seeing is believing), I added a random number to the lengths drawn from the uniform distribution \(U(-0.1,0.1)\). We can numerically prove the optimal solution is unique by adding a no-good cut to the problem and observing that the objective function starts to deteriorate immediately.

----    494 PARAMETER deterioration  objective should be worse

obj solve 1 590.266,    obj solve 2 590.276

And this is exactly what we see.

Conclusions

This particular Minimum Spanning Tree (MST) problem has multiple optimal solutions. We detected this by looking at the solution of an LP formulation: it shows the presence of dual generacy. Using a MIP formulation, we can easily enumerate all 24 optimal solutions. When we make the lengths of the arcs unique, the optimal solution is also unique. We numerically verified this. 

References

1. Minimum Spanning Trees in Math Programming Models, https://yetanothermathprogrammingconsultant.blogspot.com/2021/03/minimum-spanning-trees-in-math.html
2. What is this EPS in a GAMS solution? https://yetanothermathprogrammingconsultant.blogspot.com/2017/09/what-is-this-eps-in-gams-solution.html
3. N-Queens and Solution Pool, https://yetanothermathprogrammingconsultant.blogspot.com/2024/09/n-queens-and-solution-pool.html
4. Minimum Spanning Tree, https://en.wikipedia.org/wiki/Minimum_spanning_tree
5. Donald E. Knuth, Art of Computer Programming, Volume 4, Fascicle 4: Generating All Trees--History of Combinatorial Generation, 2006. 

Appendix: GAMS Model

$ontext     Uniqueness of optimal solutions in Minimum Spanning Tree Models.     Distances from TSPLIB     Sections:      1.  LP model based on multi-commodity flow    2.  MIP model based on lifted MTZ formulation (TSP like model)    3.  Find all optimal integer solutions using No-good cuts    4.  Find all optimal integer solutions using Cplex solution pool    5.  Verify solution is unique when distances are unique    6.  Plot single solution as Google map    7.  Plot all optimal integer solutions   $offtext   *------------------------------------------------------------------------------ * macros *------------------------------------------------------------------------------   * google map of single solution * make_google_map=1 : make google map, make_google_map=0 : skip google map $set make_google_map  0 $set google_map_key   xxxxx $set r_script         plot_map.R $set r_gdx            map_data.gdx $set r_exe            'C:\Program Files\R\R-4.4.1\bin\Rscript.exe' $set png              google_map_plot.png     * solution pool gdx $set pool_gdx         solution_pool.gdx   * HTML/SVG plots of all optimal integer solutions $set svg              svg_plots.html     *------------------------------------------------------------------------------ * data *------------------------------------------------------------------------------   set   i cities /                                  c1   'Manchester, N.H.'           c2   'Montpelier, Vt.'            c3   'Detroit, Mich.'             c4   'Cleveland, Ohio'            c5   'Charleston, W.Va.'          c6   'Louisville, Ky.'            c7   'Indianapolis, Ind.'         c8   'Chicago, Ill.'              c9   'Milwaukee, Wis.'            c10  'Minneapolis, Minn.'         c11  'Pierre, S.D.'               c12  'Bismarck, N.D.'             c13  'Helena, Mont.'              c14  'Seattle, Wash.'             c15  'Portland, Ore.'             c16  'Boise, Idaho'               c17  'Salt Lake City, Utah'       c18  'Carson City, Nevada'        c19  'Los Angeles, Calif.'        c20  'Phoenix, Ariz.'             c21  'Santa Fe, N.M.'             c22  'Denver, Colo.'              c23  'Cheyenne, Wyo.'             c24  'Omaha, Neb.'                c25  'Des Moines, Iowa'           c26  'Kansas City, Mo.'           c27  'Topeka, Kans.'              c28  'Oklahoma City, Okla.'       c29  'Dallas, Tex.'               c30  'Little Rock, Ark.'          c31  'Memphis, Tenn.'             c32  'Jackson, Miss.'             c33  'New Orleans, La.'           c34  'Birmingham, Ala.'           c35  'Atlanta, Ga.'               c36  'Jacksonville, Fla.'         c37  'Columbia, S.C.'             c38  'Raleigh, N.C.'              c39  'Richmond, Va.'              c40  'Washington, D.C.'           c41  'Boston, Mass.'              c42  'Portland, Me.'         /;   alias (i,j,k);   set coord 'coordinates' /x,y/;   table xy(i,coord) coordinates          x    y   c1   170.0  85.0   c2   166.0  88.0   c3   133.0  73.0   c4   140.0  70.0   c5   142.0  55.0   c6   126.0  53.0   c7   125.0  60.0   c8   119.0  68.0   c9   117.0  74.0  c10    99.0  83.0  c11    73.0  79.0  c12    72.0  91.0  c13    37.0  94.0  c14     6.0 106.0  c15     3.0  97.0  c16    21.0  82.0  c17    33.0  67.0  c18     4.0  66.0  c19     3.0  42.0  c20    27.0  33.0  c21    52.0  41.0  c22    57.0  59.0  c23    58.0  66.0  c24    88.0  65.0  c25    99.0  67.0  c26    95.0  55.0  c27    89.0  55.0  c28    83.0  38.0  c29    85.0  25.0  c30   104.0  35.0  c31   112.0  37.0  c32   112.0  24.0  c33   113.0  13.0  c34   125.0  30.0  c35   135.0  32.0  c36   147.0  18.0  c37   147.5  36.0  c38   154.5  45.0  c39   157.0  54.0  c40   158.0  61.0  c41   172.0  82.0  c42   174.0  87.0 ;   table d(i,j) 'from dantzig42 data set in TSPLIB'      c1  c2  c3  c4  c5  c6  c7  c8  c9 c10 c11 c12 c13 c14 c15 c16 c17 c18 c19 c20 c21 c2    8 c3   39  45 c4   37  47   9 c5   50  49  21  15 c6   61  62  21  20  17 c7   58  60  16  17  18   6 c8   59  60  15  20  26  17  10 c9   62  66  20  25  31  22  15   5 c10  81  81  40  44  50  41  35  24  20 c11 103 107  62  67  72  63  57  46  41  23 c12 108 117  66  71  77  68  61  51  46  26  11 c13 145 149 104 108 114 106  99  88  84  63  49  40 c14 181 185 140 144 150 142 135 124 120  99  85  76  35 c15 187 191 146 150 156 142 137 130 125 105  90  81  41  10 c16 161 170 120 124 130 115 110 104 105  90  72  62  34  31  27 c17 142 146 101 104 111  97  91  85  86  75  51  59  29  53  48  21 c18 174 178 133 138 143 129 123 117 118 107  83  84  54  46  35  26  31 c19 185 186 142 143 140 130 126 124 128 118  93 101  72  69  58  58  43  26 c20 164 165 120 123 124 106 106 105 110 104  86  97  71  93  82  62  42  45  22 c21 137 139  94  96  94  80  78  77  84  77  56  64  65  90  87  58  36  68  50  30 c22 117 122  77  80  83  68  62  60  61  50  34  42  49  82  77  60  30  62  70  49  21 c23 114 118  73  78  84  69  63  57  59  48  28  36  43  77  72  45  27  59  69  55  27 c24  85  89  44  48  53  41  34  28  29  22  23  35  69 105 102  74  56  88  99  81  54 c25  77  80  36  40  46  34  27  19  21  14  29  40  77 114 111  84  64  96 107  87  60 c26  87  89  44  46  46  30  28  29  32  27  36  47  78 116 112  84  66  98  95  75  47 c27  91  93  48  50  48  34  32  33  36  30  34  45  77 115 110  83  63  97  91  72  44 c28 105 106  62  63  64  47  46  49  54  48  46  59  85 119 115  88  66  98  79  59  31 c29 111 113  69  71  66  51  53  56  61  57  59  71  96 130 126  98  75  98  85  62  38 c30  91  92  50  51  46  30  34  38  43  49  60  71 103 141 136 109  90 115  99  81  53 c31  83  85  42  43  38  22  26  32  36  51  63  75 106 142 140 112  93 126 108  88  60 c32  89  91  55  55  50  34  39  44  49  63  76  87 120 155 150 123 100 123 109  86  62 c33  95  97  64  63  56  42  49  56  60  75  86  97 126 160 155 128 104 128 113  90  67 c34  74  81  44  43  35  23  30  39  44  62  78  89 121 159 155 127 108 136 124 101  75 c35  67  69  42  41  31  25  32  41  46  64  83  90 130 164 160 133 114 146 134 111  85 c36  74  76  61  60  42  44  51  60  66  83 102 110 147 185 179 155 133 159 146 122  98 c37  57  59  46  41  25  30  36  47  52  71  93  98 136 172 172 148 126 158 147 124 121 c38  45  46  41  34  20  34  38  48  53  73  96  99 137 176 178 151 131 163 159 135 108 c39  35  37  35  26  18  34  36  46  51  70  93  97 134 171 176 151 129 161 163 139 118 c40  29  33  30  21  18  35  33  40  45  65  87  91 117 166 171 144 125 157 156 139 113 c41   3  11  41  37  47  57  55  58  63  83 105 109 147 186 188 164 144 176 182 161 134 c42   5  12  55  41  53  64  61  61  66  84 111 113 150 186 192 166 147 180 188 167 140   +     c22 c23 c24 c25 c26 c27 c28 c29 c30 c31 c32 c33 c34 c35 c36 c37 c38 c39 c40 c41 c23   5 c24  32  29 c25  40  37   8 c26  36  39  12  11 c27  32  36   9  15   3 c28  36  42  28  33  21  20 c29  47  53  39  42  29  30  12 c30  61  62  36  34  24  28  20  20 c31  64  66  39  36  27  31  28  28   8 c32  71  78  52  49  39  44  35  24  15  12 c33  76  82  62  59  49  53  40  29  25  23  11 c34  79  81  54  50  42  46  43  39  23  14  14  21 c35  84  86  59  52  47  51  53  49  32  24  24  30   9 c36 105 107  79  71  66  70  70  60  48  40  36  33  25  18 c37  97  99  71  65  59  63  67  62  46  38  37  43  23  13  17 c38 102 103  73  67  64  69  75  72  54  46  49  54  34  24  29  12 c39 102 101  71  65  65  70  84  78  58  50  56  62  41  32  38  21   9 c40  95  97  67  60  62  67  79  82  62  53  59  66  45  38  45  27  15   6 c41 119 116  86  78  84  88 101 108  88  80  86  92  71  64  71  54  41  32  25 c42 124 119  90  87  90  94 107 114  77  86  92  98  80  74  77  60  48  38  32   6 ;   *------------------------------------------------------------------------------ * extra data *------------------------------------------------------------------------------   * I needed this for the visualization parameter box(*,coord) 'min/max coordinates'; box('min',coord) = smin(i,xy(i,coord)); box('max',coord) = smax(i,xy(i,coord)); display box;     sets    s(i) 'source node (can be any node)' /c40/    t(i) 'terminal nodes' ; t(i)$(not s(i)) = yes;   parameter c(i,j) 'cost directed arcs'; c(i,j) = d(i,j) + d(j,i);   set a(i,j) 'existing arcs (exclude extries without c(i,j))'; a(i,j)$c(i,j) = yes;   parameter size(*) 'number of nodes/arcs'; size('network nodes') = card(i); size('directed arcs') = card(a); option size:0; display size;    *------------------------------------------------------------------------------ * Section 1. LP model based on multi-commodity flow * This is a large LP *------------------------------------------------------------------------------   parameter b(i,k) 'right-hand side (exogenous net inflow)'; b(s,t) = 1; b(t,t) = -1;   variables    z         'objective'    x(i,j)    '0-1 variable indicating tree'    y(i,j,k)  'flow variables (last index is terminal node)' ; binary variable x,y;   equations    objective      'minimize total cost'    nodebal(i,k)   'node balance'    compare(i,j,k) 'implication x=0 => y=0' ;   objective..    z =e= sum(a, c(a)*x(a));   nodebal(i,t)..    sum(a(i,j), y(a,t)) =e= sum(a(j,i), y(a,t)) + b(i,t);   compare(a,t)..    y(a,t) =l= x(a);   * solve as LP model m1 /objective,nodebal,compare/; solve m1 minimizing z using rmip; display    y.l, x.l, z.l,    "observe there are many EPS values in the marginals",    nodebal.m;    *------------------------------------------------------------------------------ * Section 2. MIP model based on lifted MTZ formulation (TSP like model) *------------------------------------------------------------------------------   scalar n 'number of nodes'; n = card(i);   positive variables   v(i)      'ordering of nodes (>=hops from source node)' ; v.fx(s) = 0; v.lo(t) = 1; v.up(t) = n-1;   equations    degree1        'out-degree of source node'    degree2(k)     'in-degree of all other nodes'    mst2(i,j)      'lifted MST constraint' ;   degree1..    sum(a(s,j), x(a)) =g= 1;   degree2(t)..    sum(a(j,t), x(a)) =e= 1;   mst2(a(i,j))$(t(i) and t(j))..    v(i) - v(j) + (n-1)*x(i,j) + (n-3)*x(j,i) =l= n-2;    model m4 /objective,degree1,degree2,mst2/; solve m4 minimizing z using mip; display x.l;   *------------------------------------------------------------------------------ * save solution for Google map *------------------------------------------------------------------------------ parameter tree(i,j,*) 'optimal solution for Google map'; tree(i,j,'ordi')$(x.l(i,j)>0.5) = ord(i); tree(i,j,'ordj')$(x.l(i,j)>0.5) = ord(j); display tree;   *------------------------------------------------------------------------------ * Section 3. Find all optimal integer solutions using No-good cuts * two versions: *   1. general formulation *   2. Simpler version (works because of structure of the model) *------------------------------------------------------------------------------   set   p 'number of solutions' /sol1*sol50/   pcut(p) 'active cuts'    psols(p) 'found solutions (saved for plotting)'  ;   parameters    xstore(p,i,j) 'stored values'    xsols(p,i,j)  'found solutions (saved for plotting)'    zsols(p)      'objective values' ;   equations    cut1(p)       'general cut'    cut2(p)       'specialized cut' ;   * general cuts (using dynamic set pcut) cut1(pcut).. sum(a, xstore(pcut,a)*x(a)) - sum(a, (1-xstore(pcut,a))*x(a)) =l= sum(a,xstore(pcut,a))-1;   * specialized cuts (using dynamic set pcut) * we know sum(a,x(a))=n-1 so we can do: cut2(pcut).. sum(a, xstore(pcut,a)*x(a))  =l= n-2;     Models    m4cut1 /m4,cut1/    m4cut2 /m4,cut2/ ;   * select one of the above models $set m4cut m4cut2   * reduce output %m4cut%.solprint = %solprint.silent%; * faster solves option solvelink = %solveLink.loadLibrary%;     scalar zbest 'best objective value'; pcut(p) = no; xstore(p,a) = 0; loop(p,    solve %m4cut% minimizing z using mip; * first solve gives the best obj      zbest$(ord(p)=1) = z.l; * stop if obj deteriorates (use a tolerance)    break$(z.l > zbest+0.01); * add solution to xstore/pcut      pcut(p) = yes;    xstore(p,a) = round(x.l(a));  );   parameter info(*); info('first objective') = zbest; info('last objective') = z.l; info('number of optimal integer solutions') = card(pcut); display info;   * save for plotting psols(p) = pcut(p); xsols(p,a) = xstore(p,a); zsols(p) = sum(a,c(a)*xsols(p,a));   *------------------------------------------------------------------------------ * Section 4. Find all optimal integer solutions using Cplex solution pool * objective is integer valued so we can safely use SolnPoolAGap=0.1 * There are 24 solutions * with SolnPoolAGap=0 and threads=1 we get 23 solutions *------------------------------------------------------------------------------   $onecho > cplex.opt SolnPoolAGap=0.1 solnpoolintensity=4 solnpoolpop=2 populatelim=100000 solnpoolmerge=%pool_gdx%   *to reproduce Cplex bug use: *SolnPoolAGap=0 *threads=1 $offecho   m4.optfile=1; solve m4 minimizing z using mip;   * we should see: * --- Dumping 24 solutions from the solution pool...   *-------------------------------------------------------------------- * load all solutions *--------------------------------------------------------------------   Sets    index0'register elements' /soln_m4_p1*soln_m4_p100/    index(index0)  ;   parameter allsols(index0,i,j);   execute_load "%pool_gdx%" index,allsols=x; display index,allsols;     *-------------------------------------------------------------------- * Section 5. Verify solution is unique when distances are unique * this is a bit messy *--------------------------------------------------------------------   alias (i,ii),(j,jj);   Set   gt(i,j) 'only i>j part (lower triangular part)'   lt(i,j) 'only i<j part (upper triangular part)'   gt2(i,j,ii,jj) '(ii,jj)>(i,j) meaning (ii=i and jj>j) or ii>i' ; gt(i,j) = ord(i) > ord(j); lt(i,j) = ord(i) < ord(j); gt2(i,j,ii,jj) = gt(i,j) and gt(ii,jj) and                  ((ord(ii)=ord(i) and ord(jj)>ord(j)) or ord(ii)>ord(i));       * count number of duplicates (there are a ton of them) scalar ndup 'number of duplicate c(i,j)'; set isdup(i,j,ii,jj) 'c(i,j) and c(ii,jj) are duplicates'; isdup(gt2(i,j,ii,jj)) = c(i,j)=c(ii,jj); ndup = card(isdup); option isdup:0:0:4; display "original c(i,j), i>j",ndup,isdup;   * make unique by perturbing c(gt) = c(gt) + uniform(-0.1,0.1); c(lt(i,j)) = c(j,i); display c;   * now the number of duplicate c(i,j) should be zero (except for symmetry) isdup(gt2(i,j,ii,jj)) = c(i,j)=c(ii,jj); ndup = card(isdup); display "perturbed c(i,j), i>j",ndup,isdup;   * do first solve without cuts pcut(p) = no; xstore(p,a) = 0; %m4cut%.solprint = %solprint.on%; solve %m4cut% minimizing z using mip;   * second solve with single cut should have an obj that is worse pcut('sol1') = yes; xstore('sol1',a) = round(x.l(a)); zbest = z.l; solve %m4cut% minimizing z using mip;   parameter deterioration(*) 'objective should be worse'; deterioration('obj solve 1') = zbest; deterioration('obj solve 2') = z.l; display deterioration; abort$(abs(z.l-zbest)<0.001) "expected deterioration of objective"; display "Optimal solution is unique!!"     *-------------------------------------------------------------------- * Section 6. Plot single solution as Google map * make sure your R version has installed the * following packages: *   tidyverse *   ggmap *   gamstransfer * output: %png% i.e. google_map_plot.png   * Note: you need a Google Map Key for this.   *--------------------------------------------------------------------   if(%make_google_map%=1,      execute_unload "%r_gdx%",i,tree;      execute 'echo ====== R script ====== & "%r_exe%" %r_script%';      executetool 'win32.ShellExecute "%png%"'; );     $onecho > %script%   # # This script is called from GAMS #   library(tidyverse) library(ggmap) library(gamstransfer)   gdx = Container$new("%r_gdx%")   gdx$describeSets() gdx$describeParameters()   cat("==== set i =====\n") df <- gdx["i"]$records df   cat("==== parameter tree =====\n") tree <- gdx["tree"]$records tree   cat("====pivot=====\n") df2 <- pivot_wider(tree,names_from=uni,values_from=value) df2    register_google(key = "%google_map_key%") map <- get_map(location="united states", zoom=4, maptype = "terrain",                source="google",color="color")   location <- geocode(df$element_text) df$lat <- location$lat df$lon <- location$lon df   df2$lat1 <- location$lat[df2$ordi] df2$lat2 <- location$lat[df2$ordj] df2$lon1 <- location$lon[df2$ordi] df2$lon2 <- location$lon[df2$ordj] df2   ggmap(map) +   geom_point(data=df, aes(x=lon, y=lat), color="navy", size=1) +   geom_segment(data=df2, aes(x=lon1,y=lat1,xend=lon2,yend=lat2),color="navy",linewidth=0.5)   ggsave("%png%") $offecho   *-------------------------------------------------------------------- * Section 7.  Plot all optimal integer solutions * make a plot of all optimal solutions (all 26 of them) * Output is an HTML file (svg_plots.html) *--------------------------------------------------------------------   scalar colnum 'counter' /1/;   file f /%svg%/; put f;   put '<style>table,th,td {border-collapse: collapse;}</style>'/; put '<h2>Multiple Optimal Solutions in Minimum Spanning Tree problem</h2>'/; put '<p>This problem has integer lengths with duplicates. This produces ',card(psols):0:0,' solutions</p>'/;   put '<table border=1>'/; put '<tr>'/;   loop(psols(p),    put '<td>'/;    put '<svg height="200" width="300" viewbox="0 10 180 115">'/;    loop(a(i,j)$(xsols(p,a)=1),       put '<line x1="',xy(i,'x'):0:0,'" y1="',xy(i,'y'):0:0,'" x2="',xy(j,'x'):0:0,'" y2="',xy(j,'y'):0:0,'" style="stroke:blue;stroke-width:1" />'/;    );    loop(i,       put '<circle cx="',xy(i,'x'):0:0,'" cy="',xy(i,'y'):0:0,'" r="1" fill="red" stroke="red" />'/;    );    put '</svg><br>'/;    put 'Solution:',ord(p):0:0,' of ',card(psols):0:0,' Length=',zsols(p):0:1/;    put '</td>'/;      if (colnum=6 and ord(p)<card(psols),       colnum = 0;       put '</tr>'/;       put '<tr>'/;      );    colnum=colnum+1; ); put '</tr>'/; put '</table>'/;   putclose;   executetool 'win32.ShellExecute "%svg%"';