// A model for minimum spanning trees. To be used together with an OPL
// script. The idea is to iteratively add cuts as we find them in the
// current solution. We start by finding the cheapest n-1 edges. If they
// form a spanning tree we are done. If not we find (in the script)
// connected  components in the graph consisting of the edges chosen.
// We then add these as cuts in the model with the condition that inside
// a component of size m there can be at most m-1 edges and each component
// must have an edge out of it.
// All found during the execution are kept, meaning that we add more and
// more cuts.
 
// vertices

int n =...;
range vertices 1..n;

//edges

struct edge {int i; int j;};
{edge} Edges                 = {<i,j> | ordered i,j in vertices};
int cost[Edges]              = ...;

// Connected component data obtained from the script
// via the import command.

import Open componentset;
range  comprange [componentset.low..componentset.up];
int    compsize[k in comprange] = card(componentset[k]);

//Decision variables

var int x[Edges] in 0..1;
var int y[vertices,vertices] in 0..1;

// CPLEX settings
setting fracCuts = 2;

/****************************************************************************************
*
*  Model
*
*****************************************************************************************/

//objective

minimize sum (<i,j> in Edges) cost[<i,j>]*x[<i,j>]
subject to {
   forall (i,j in vertices) y[i,j]=y[j,i];
   forall (i in vertices ) y[i,i]=0;
   forall (<i,j> in Edges) y[i,j]=x[<i,j>];

// precisely n-1 edges used
sum (<i,j> in Edges) x[<i,j>] = n -1;

// connectivity constraints: at most size -1 edges inside a componentset and its complement

forall (k in [componentset.low..componentset.up]){
    sum (i,j in componentset[k] : i<j) x[<i,j>] <=compsize[k]-1;
    sum (i,j in vertices :(i not in componentset[k] & j not in componentset[k] & i<j ))x[<i,j>] <= n-compsize[k]-1;
};
// at least one edge out of each component

forall (k in [componentset.low..componentset.up])
    sum (i,j in vertices : (i in componentset[k] & j not in componentset[k])) y[i,j] >= 1;
};