$title Multicast Problem, with log(x+1) cost function replaced by concave piecewise linear approximation

$ontext

We consider a network in which the N nodes are connected in a straight
line. Each node has a fixed demand for content, which emanates from one of
M servers which are located at nodes to be determined. Nonserver nodes
receive their content through the network from nearby server nodes. The
cost of sending content along an arc of the network is a concave function
of the amount of traffic x on the arc; namely, log(x+1).

The aim is to choose the locations of the M servers in such a way that the
total cost of the network traffic (obtained by summing the costs over all
the arcs) is minimized. The concave function log(x+1) is approximated by a
piecewise linear function, implemented using SOS2 variables.

(Thanks to Jussara Almeida and Mary Vernon, UW-Madison)

$offtext

$setglobal N 10
$setglobal M 3
$setglobal DEMAND 10

option limrow=0, limcol=0, solprint=off;

* limit on the number of servers
scalar numServers/%M%/;

* define node set and its aliases
set nodes /1*%N%/;
alias (i,j,nodes);

* define the arcs
set arcs(i,j);

arcs(i,j) = no;
arcs(i,i+1)=yes;

* define client nodes and the set of possible server sites
* (allow all nodes to be possible servers)
set clients(nodes) /1*%N%/;
set possibleServers(nodes) /1*%N%/;


* define demands at the nodes - uniform in this example
parameter demand(nodes);
demand(nodes) = %DEMAND%;

* "big M" limit: make it larger than the acticipated total flow on
* any arc by setting it to the sum of demands
scalar bigM;
bigM = sum(nodes, demand(nodes));

* abscissae of the piecewise-linear bw cost function
* pieces: 0->10, 10->100, 100-100000
set pieces /1*4/;

* define abscissae and ordinates of piecewise linear function
parameter Xflow(pieces), Xbw(pieces);
Xflow('1') = 0;
Xbw('1') = 0;
Xflow('2') = 10;
Xbw('2') = log(10+1);
Xflow('3') = 100;
Xbw('3') = log(100+1);
Xflow('4') = 100000;
Xbw('4') = log(100000+1);

* for each arc (i,j), gamma(i,j,.) contains the SOS2 variables
* that define the piecewise linear approximation to the log(x+1)
* cost function
sos2 variables gamma(i,j,pieces)
;

binary variable servers(nodes) indicates which server locations are used 
;
variable flow(i,j) flow from i to j (negative indicates flow in reverse dir)
;
positive variable
    absflow(i,j)   absolute value of the flow
    arcPWcost(i,j) cost of flow along each arc
;

variable
    PWcost        piecewise linear objective
;

equations
    objectivePW      sum of piecewise objectives over all arcs
    balance(nodes)   balance constraints at nodes that aren't server candidates
    balanceLo(nodes)
    balanceUp(nodes) balance constraints at server candidate nodes
    EQabsflowNeg, EQabsflowPos  define absflow
    serverLimit      limit on number of servers
    EQpiecewise      constraints on SOS2 variables
    EQxflow	     figure flow from pieces
    EQarcPWcost      figure cost of flow from piecewise linear function
;

* restriction on the total number of servers allowed
serverLimit..
    sum(possibleServers, servers(possibleServers)) =l= numServers;

objectivePW..
    PWcost =e= sum((i,j)$arcs(i,j), arcPWcost(i,j));

* balance constraints at client nodes
balance(nodes)$(not possibleServers(nodes))..
    sum(i$arcs(i,nodes), flow(i,nodes))
   -sum(j$arcs(nodes,j), flow(nodes,j))
    =e= demand(nodes);

* balance constraints at possible servers: lower bounds. Essentially
* this constraint is "turned off" if the node is chosen as a server.
balanceLo(nodes)$possibleServers(nodes)..
    sum(i$arcs(i,nodes), flow(i,nodes))
   -sum(j$arcs(nodes,j), flow(nodes,j))
      =g= demand(nodes) - servers(nodes)*bigM;

* balance constraints at possible servers: upper bounds. Similarly,
* this constraint is "turned off" if the node is chosen as a server.
balanceUp(nodes)$possibleServers(nodes)..
    sum(i$arcs(i,nodes), flow(i,nodes))
   -sum(j$arcs(nodes,j), flow(nodes,j))
      =l= demand(nodes) + servers(nodes)*bigM;

* define the absolute flow
EQabsflowNeg(i,j)$arcs(i,j)..   absflow(i,j) =g= -flow(i,j);
EQabsflowPos(i,j)$arcs(i,j)..   absflow(i,j) =g=  flow(i,j);

* constraints on SOS2 variables 
EQpiecewise(i,j)$arcs(i,j)..
sum(pieces,gamma(i,j,pieces))=e= 1;

* figure flow from i to j from pieces
EQxflow(i,j)$arcs(i,j)..
  absflow(i,j) =e= sum(pieces,gamma(i,j,pieces)*Xflow(pieces));

EQarcPWcost(i,j)$arcs(i,j)..
  arcPWcost(i,j) =e= sum(pieces,gamma(i,j,pieces)*Xbw(pieces));

model toynetPW/objectivePW, balance, balanceLo, balanceUp, EQabsflowNeg, EQabsflowPos, serverLimit, EQpiecewise, EQxflow, EQarcPWcost/;

* define parameters for the mip solver
toynetPW.optca = 0;
toynetPW.optcr = 0.000001;
toynetPW.reslim = 1000000;
toynetPW.iterlim = 10000000;
toynetPW.workspace = 200;

solve toynetPW using mip minimizing PWcost;
display PWcost.l, flow.l, servers.l;