I am trying to extend the program on the following link https://github.com/shellei503/Linear-Programming-Examples/blob/master/9.3%20(The%20Shortest%20Path%20Problem)/prob_1/9.3-1.ipynb for single-source and multiple destinations.

Some modifications to the existing network include that each edge in the network is bidirectional. I have added the following constraints to the existing code to achieve the objective:

m.add_constraint(m.sum(x[(i,j)] for i,j in last_arc if j == final_END_NODE)== 1, ctname='last_arc')
for dd in remDst:
   m.add_constraint(m.sum(x[(i,j)] for i,j in subEdges if j==dd)>= 1, ctname='last_arc')
m.add_constraint(m.sum(x[(i,j)] for i,j in subEdges if i==src)== 1, ctname='start_arc_F')
m.add_constraint(m.sum(x[(i,j)] for i,j in subEdges if j==src)== 0, ctname='start_arc_B')​

where final_END_NODE is computed based on the distance from the source node to the farthest destination, remDst includes the remaining destinations, and subEdges includes the edges in the network which are bidirectionally visitable.

The output shows discontinuous edges while computing a single shortest path from the source to all destinations. For example, O is the source node and A, C, and E are the destination nodes. The output shows the edges:

O-C

C-E

A-B

B-A

The solution misses the edge from E-B. What should be the constraint to address this problem. Any help in this regard is appreciated.

minimize sum (<i,j> in Edges) dist[<i,j>]*x[<i,j>];​

minimize sum (<i,j> in Edges) dist[<i,j>]*x[<i,j>]-max(j in Cities:<1,j> in Edges) dist[<1,j>]*x[<1,j>];​

using CP;

execute
{
  cp.param.timelimit=10;
}

{string} nodes={"O","A","B","C","D","E","T"};

tuple edge
{
  key string o;
  key string d;
  int time;
}

{edge} edges with o,d in nodes=
{
  <"O","A",40>,
  <"O","B",60>,
  <"O","C",50>,
  <"A","B",10>,
  <"C","B",20>,
  <"A","D",70>,
  <"B","D",55>,
  <"B","E",40>,
  <"D","E",10>,
  <"D","T",60>,
  <"E","T",80>
};

string origin="O";
{string} targets={"A","C"};
int bigvalue=1000;
int repeat=3;
range ranks=1..repeat;

{edge} edgeswithsym=edges union {<d,o,t> | <o,d,t> in edges};

{edge} transitions=edgeswithsym union {<o,d,bigvalue> | o,d in nodes: <o,d> not in edgeswithsym};

tuple triplet {int id1; int id2; int value;};
{triplet} M = {<ord(nodes,tr.o)+1,ord(nodes,tr.d)+1,tr.time> | tr in transitions};

assert card(transitions)==card(nodes)*(card(nodes));

// Interval for visiting a node
dvar interval itvs[nodes][ranks] optional;

// Sequence means visits
dvar sequence seq in all(n in nodes,r in ranks)itvs[n][r] types 
all(n in nodes,r in ranks)(ord(nodes,n)+1); 

// First we want to minimize the time for latest visit
// Second we want to minimize present intervals
minimize 
staticLex(max(t in targets) min(r in ranks) startOf(itvs[t,r],bigvalue),
sum(n in nodes,r in ranks) presenceOf(itvs[n][r]));
subject to
{
  // We visit origin at time 0
  startOf(itvs[origin,1])==0;
  
  // origin and destinations should be present
  presenceOf(itvs[origin,1])==1;
  forall(t in targets) presenceOf(itvs[t,1])==1;
  
  // noOverlap constraint will enforce time matrix
  noOverlap(seq,M,true);
  
  // break sym
  forall(n in nodes) forall(r in ranks:r!=1) 
    presenceOf(itvs[n,r-1])>=presenceOf(itvs[n,r]);
  
}

int nbActiveIntervals=sum(n in nodes,r in ranks) presenceOf(itvs[n][r]);
range R=1..nbActiveIntervals;

// Display solution
execute {
  

  writeln(origin);
  var s=seq.first(); 
  for(var i in R) 
  {
   if (i!=nbActiveIntervals)
      writeln(Opl.item(nodes,-1+Opl.typeOfNext(seq,s,bigvalue,0)));
   s=seq.next(s); 
   
  } 
   
}

/*

gives

O
A
B
C

*/
​