Implementation

The following Mosel model roadworks.mos implements the robust optimization model from the previous section. Notice the use of cardinality on the uncertains delay to limit the total number of occurrences of roadworks to at most NWORK (the parameter for N) with a default value of 2. The uncertains that are used in such a cardinality constraint need to have explicit lower and upper bounds set. The solution reporting displays the selected arcs and the values of the uncertains for which a nonzero delay is chosen by the opponent.

model "road network"
 uses "mmxprs", "mmrobust"

 parameters
   DATAFILE="roads_9.dat"
   NWORK = 2                                     ! Number of roadworks
 end-parameters

 declarations
   Nodes: range                                  ! Set of nodes
   ARC: array(Arcs:range,1..2) of integer        ! Arc origins/destinations
   LEN,MAXDELAY: array(Arcs) of real             ! Length and max delay per arc
   Source,Sink: integer                          ! Source and sink node numbers

   use: array(Arcs) of mpvar                     ! 1 iff arc is used
   TotalLength: robctr                           ! Objective function
   delay: array(Arcs) of uncertain               ! Uncertain delay
 end-declarations

 !**** Input datafile ****
 initializations from DATAFILE
   Nodes  Arcs  Source  Sink  ARC
   [LEN,MAXDELAY] as "ArcData"
 end-initializations

!**** Robust problem formulation ****
 forall(a in Arcs) use(a) is_binary

! Sink and source of flow
 sum(a in Arcs | ARC(a,1)=Source) use(a)=1
 sum(a in Arcs | ARC(a,2)=Sink) use(a)=1
 sum(a in Arcs | ARC(a,2)=Source) use(a)=0
 sum(a in Arcs | ARC(a,1)=Sink) use(a)=0

! Flow balance in intermediate nodes
 forall(n in Nodes-{Source,Sink})
   sum(a in Arcs | ARC(a,2)=n) use(a) = sum(a in Arcs | ARC(a,1)=n) use(a)

! Random construction work on NWORK arcs
 forall(a in Arcs) delay(a) <= MAXDELAY(a)
 forall(a in Arcs) delay(a) >= 0
 cardinality(union(a in Arcs) {delay(a)}, NWORK)

! Shortest path length
 TotalLength:= sum(a in Arcs) (LEN(a)+delay(a))*use(a)

!**** Solving ****
 minimize(TotalLength)

! Solution reporting
 if getprobstat=XPRS_OPT then	
   writeln("Robust shortest path: ", getobjval)
   writeln("Delay: ", sum(a in Arcs | use(a).sol>0) delay(a).sol)	
   forall(a in Arcs | use(a).sol>0)
     writeln(if(ARC(a,1)=Source, "Source", string(ARC(a,1))), " -> ",
             if(ARC(a,2)=Sink, "Sink", string(ARC(a,2))), ";  ")
   forall(a in Arcs | delay(a).sol>0)
     writeln("Delay on (", ARC(a,1), ",", ARC(a,2), "): ", delay(a).sol)
 else
   writeln("No solution")
 end-if    	

end-model