Minimum Cost Flow - network optimization in R

In case of interest, here is how I ended up solving this problem. I used an edge dataframe with to nodes, from nodes, cost property and capacity property for each edge to create a constraint matrix. Subsequently, I fed this into a linear optimisation using the lpSolve package. Step-by-step below.

Start with the edgelist dataframe from my example above

library(magrittr)# Add edge IDedgelist$ID <- seq(1, nrow(edgelist))glimpse(edgelist)

Looks like this:

Observations: 16Variables: 4$ from     <dbl> 1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 5, 5, 6, 6, 7, 8$ to       <dbl> 2, 3, 4, 5, 6, 4, 5, 6, 7, 8, 7, 8, 7, 8, 9, 9$ capacity <dbl> 20, 30, 99, 99, 99, 99, 99, 99, 99, 99, 99, 99, 99, 99, 99, 99$ cost     <dbl> 0, 0, 1, 2, 3, 4, 3, 2, 3, 2, 3, 4, 5, 6, 0, 0$ ID       <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16

Create constraints matrix

createConstraintsMatrix <- function(edges, total_flow) {  # Edge IDs to be used as names  names_edges <- edges$ID  # Number of edges  numberof_edges <- length(names_edges)  # Node IDs to be used as names  names_nodes <- c(edges$from, edges$to) %>% unique  # Number of nodes  numberof_nodes <- length(names_nodes)  # Build constraints matrix  constraints <- list(    lhs = NA,    dir = NA,    rhs = NA)  #' Build capacity constraints ------------------------------------------------  #' Flow through each edge should not be larger than capacity.  #' We create one constraint for each edge. All coefficients zero  #' except the ones of the edge in question as one, with a constraint  #' that the result is smaller than or equal to capacity of that edge.  # Flow through individual edges  constraints$lhs <- edges$ID %>%    length %>%    diag %>%    set_colnames(edges$ID) %>%    set_rownames(edges$ID)  # should be smaller than or equal to  constraints$dir <- rep('<=', times = nrow(edges))  # than capacity  constraints$rhs <- edges$capacity  #' Build node flow constraints -----------------------------------------------  #' For each node, find all edges that go to that node  #' and all edges that go from that node. The sum of all inputs  #' and all outputs should be zero. So we set inbound edge coefficients as 1  #' and outbound coefficients as -1. In any viable solution the result should  #' be equal to zero.  nodeflow <- matrix(0,                     nrow = numberof_nodes,                     ncol = numberof_edges,                     dimnames = list(names_nodes, names_edges))  for (i in names_nodes) {    # input arcs    edges_in <- edges %>%      filter(to == i) %>%      select(ID) %>%      unlist    # output arcs    edges_out <- edges %>%      filter(from == i) %>%      select(ID) %>%      unlist    # set input coefficients to 1    nodeflow[      rownames(nodeflow) == i,      colnames(nodeflow) %in% edges_in] <- 1    # set output coefficients to -1    nodeflow[      rownames(nodeflow) == i,      colnames(nodeflow) %in% edges_out] <- -1  }  # But exclude source and target edges  # as the zero-sum flow constraint does not apply to these!  # Source node is assumed to be the one with the minimum ID number  # Sink node is assumed to be the one with the maximum ID number  sourcenode_id <- min(edges$from)  targetnode_id <- max(edges$to)  # Keep node flow values for separate step below  nodeflow_source <- nodeflow[rownames(nodeflow) == sourcenode_id,]  nodeflow_target <- nodeflow[rownames(nodeflow) == targetnode_id,]  # Exclude them from node flow here  nodeflow <- nodeflow[!rownames(nodeflow) %in% c(sourcenode_id, targetnode_id),]  # Add nodeflow to the constraints list  constraints$lhs <- rbind(constraints$lhs, nodeflow)  constraints$dir <- c(constraints$dir, rep('==', times = nrow(nodeflow)))  constraints$rhs <- c(constraints$rhs, rep(0, times = nrow(nodeflow)))  #' Build initialisation constraints ------------------------------------------  #' For the source and the target node, we want all outbound nodes and  #' all inbound nodes to be equal to the sum of flow through the network  #' respectively  # Add initialisation to the constraints list  constraints$lhs <- rbind(constraints$lhs,                           source = nodeflow_source,                           target = nodeflow_target)  constraints$dir <- c(constraints$dir, rep('==', times = 2))  # Flow should be negative for source, and positive for target  constraints$rhs <- c(constraints$rhs, total_flow * -1, total_flow)  return(constraints)}constraintsMatrix <- createConstraintsMatrix(edges, 30)

Should result in something like this

$lhs1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 161       1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  02       0  1  0  0  0  0  0  0  0  0  0  0  0  0  0  03       0  0  1  0  0  0  0  0  0  0  0  0  0  0  0  04       0  0  0  1  0  0  0  0  0  0  0  0  0  0  0  05       0  0  0  0  1  0  0  0  0  0  0  0  0  0  0  06       0  0  0  0  0  1  0  0  0  0  0  0  0  0  0  07       0  0  0  0  0  0  1  0  0  0  0  0  0  0  0  08       0  0  0  0  0  0  0  1  0  0  0  0  0  0  0  09       0  0  0  0  0  0  0  0  1  0  0  0  0  0  0  010      0  0  0  0  0  0  0  0  0  1  0  0  0  0  0  011      0  0  0  0  0  0  0  0  0  0  1  0  0  0  0  012      0  0  0  0  0  0  0  0  0  0  0  1  0  0  0  013      0  0  0  0  0  0  0  0  0  0  0  0  1  0  0  014      0  0  0  0  0  0  0  0  0  0  0  0  0  1  0  015      0  0  0  0  0  0  0  0  0  0  0  0  0  0  1  016      0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  12       1  0 -1 -1 -1  0  0  0  0  0  0  0  0  0  0  03       0  1  0  0  0 -1 -1 -1  0  0  0  0  0  0  0  04       0  0  1  0  0  1  0  0 -1 -1  0  0  0  0  0  05       0  0  0  1  0  0  1  0  0  0 -1 -1  0  0  0  06       0  0  0  0  1  0  0  1  0  0  0  0 -1 -1  0  07       0  0  0  0  0  0  0  0  1  0  1  0  1  0 -1  08       0  0  0  0  0  0  0  0  0  1  0  1  0  1  0 -1source -1 -1  0  0  0  0  0  0  0  0  0  0  0  0  0  0target  0  0  0  0  0  0  0  0  0  0  0  0  0  0  1  1$dir[1] "<=" "<=" "<=" "<=" "<=" "<=" "<=" "<=" "<=" "<=" "<=" "<=" "<=" "<="[15] "<=" "<=" "==" "==" "==" "==" "==" "==" "==" "==" "=="$rhs[1]  20  30  99  99  99  99  99  99  99  99  99  99  99  99  99  99   0[18]   0   0   0   0   0   0 -30  30

Feed constraints to lpSolve for the ideal solution

library(lpSolve)# Run lpSolve to find best solutionsolution <- lp(  direction = 'min',  objective.in = edgelist$cost,  const.mat = constraintsMatrix$lhs,  const.dir = constraintsMatrix$dir,  const.rhs = constraintsMatrix$rhs)# Print vector of flow by edgesolution$solution# Include solution in edge dataframeedgelist$flow <- solution$solution

Now we can convert our edges to a graph object and plot the solution

library(igraph)g <- edgelist %>%  # igraph needs "from" and "to" fields in the first two colums  select(from, to, ID, capacity, cost, flow) %>%  # Make into graph object  graph_from_data_frame()# Get some colours in to visualise costE(g)$color[E(g)$cost == 0] <- 'royalblue'E(g)$color[E(g)$cost == 1] <- 'yellowgreen'E(g)$color[E(g)$cost == 2] <- 'gold'E(g)$color[E(g)$cost >= 3] <- 'firebrick'# Flow as edge size,# cost as colourplot(g, edge.width = E(g)$flow)

Hope it's interesting / useful :)

I was looking for a function but didn't success. The original function calls another: res <- .Call("R_igraph_maxflow", graph, source - 1, target - 1, capacity, PACKAGE = "igraph")

And I don't know how to deal with it.

For the moment I inverted the cost path values in order to use the same function in oposite direction:

E2 <- E # create another tableE2[, 3] <- max(E2[, 3]) + 1 - E2[, 3] # invert valuesE2     from to capacity[1,]    1  3        8[2,]    3  4       10[3,]    4  2        9[4,]    1  5       10[5,]    5  6        9[6,]    6  2        1g2 <- graph_from_data_frame(as.data.frame(E2)) # create 2nd graph# Get maximum flowm1 <- max_flow(g1, source=V(g1)["1"], target=V(g1)["2"])m2 <- max_flow(g2, source=V(g2)["1"], target=V(g2)["2"])m1$partition2 # Route on maximal cost+ 4/6 vertices, named:[1] 4 5 6 2m2$partition2 # Route on minimal cost+ 3/6 vertices, named:[1] 3 4 2

I draw in paper the graph and my code agree with manual solutionThis method should be tested with real know values, as I mentioned in the comment