r - 最小成本流 - R 中的网络优化

Question

我正在尝试在 R 中实现“Minimum Cost Network Flow”运输问题解决方案。

我知道这可以使用类似 lpSolve 的东西从头开始实现。但是，我看到“Maximum Flow”有一个方便的 igraph 实现。这样一个预先存在的解决方案会方便得多，但我找不到 Minimum Cost 的等效函数。

有没有 igraph 函数可以计算最小成本网络流量解决方案，或者有没有办法将 igraph::max_flow 函数应用于最小成本问题？

igraph 网络示例:

library(tidyverse)
library(igraph)

edgelist <- data.frame(
  from = c(1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 5, 5, 6, 6, 7, 8),
  to = c(2, 3, 4, 5, 6, 4, 5, 6, 7, 8, 7, 8, 7, 8, 9, 9),
  capacity = c(20, 30, 99, 99, 99, 99, 99, 99, 99, 99, 99, 99, 99, 99, 99, 99),
  cost = c(0, 0, 1, 2, 3, 4, 3, 2, 3, 2, 3, 4, 5, 6, 0, 0))

g <- graph_from_edgelist(as.matrix(edgelist[,c('from','to')]))

E(g)$capacity <- edgelist$capacity
E(g)$cost <- edgelist$cost

plot(g, edge.label = E(g)$capacity)
plot(g, edge.label = E(g)$cost)

这是一个具有定向边、“源节点”(1) 和“汇节点”(9) 的网络。每条边都有一个“容量”(这里通常表示为 99 表示无限制)和一个“成本”(一个单位流经该边的成本)。我想找到流的整数向量 (x, length = 9)，它可以在通过网络传输预定义流时将成本降至最低(假设 50 个单位，从节点 1 到节点 9)。

免责声明:this post问了一个类似的问题，但没有得到令人满意的答案，而且已经过时了(2012 年)。

Answer 1

如有兴趣，以下是我最终解决此问题的方法。我使用了一个边缘数据框，其中包含 to 节点、from 节点、cost 属性和 capacity 属性，用于创建每个边约束矩阵。随后，我使用 lpSolve 包将其输入到线性优化中。下面逐步介绍。

---

从上面示例中的边缘列表数据框开始

library(magrittr)

# Add edge ID
edgelist$ID <- seq(1, nrow(edgelist))

glimpse(edgelist)

看起来像这样:

Observations: 16
Variables: 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

---

创建约束矩阵

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)

结果应该是这样的

$lhs
1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16
1       1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
2       0  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0
3       0  0  1  0  0  0  0  0  0  0  0  0  0  0  0  0
4       0  0  0  1  0  0  0  0  0  0  0  0  0  0  0  0
5       0  0  0  0  1  0  0  0  0  0  0  0  0  0  0  0
6       0  0  0  0  0  1  0  0  0  0  0  0  0  0  0  0
7       0  0  0  0  0  0  1  0  0  0  0  0  0  0  0  0
8       0  0  0  0  0  0  0  1  0  0  0  0  0  0  0  0
9       0  0  0  0  0  0  0  0  1  0  0  0  0  0  0  0
10      0  0  0  0  0  0  0  0  0  1  0  0  0  0  0  0
11      0  0  0  0  0  0  0  0  0  0  1  0  0  0  0  0
12      0  0  0  0  0  0  0  0  0  0  0  1  0  0  0  0
13      0  0  0  0  0  0  0  0  0  0  0  0  1  0  0  0
14      0  0  0  0  0  0  0  0  0  0  0  0  0  1  0  0
15      0  0  0  0  0  0  0  0  0  0  0  0  0  0  1  0
16      0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  1
2       1  0 -1 -1 -1  0  0  0  0  0  0  0  0  0  0  0
3       0  1  0  0  0 -1 -1 -1  0  0  0  0  0  0  0  0
4       0  0  1  0  0  1  0  0 -1 -1  0  0  0  0  0  0
5       0  0  0  1  0  0  1  0  0  0 -1 -1  0  0  0  0
6       0  0  0  0  1  0  0  1  0  0  0  0 -1 -1  0  0
7       0  0  0  0  0  0  0  0  1  0  1  0  1  0 -1  0
8       0  0  0  0  0  0  0  0  0  1  0  1  0  1  0 -1
source -1 -1  0  0  0  0  0  0  0  0  0  0  0  0  0  0
target  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

---

将约束条件馈送到 lpSolve 以获得理想的解决方案

library(lpSolve)

# Run lpSolve to find best solution
solution <- lp(
  direction = 'min',
  objective.in = edgelist$cost,
  const.mat = constraintsMatrix$lhs,
  const.dir = constraintsMatrix$dir,
  const.rhs = constraintsMatrix$rhs)

# Print vector of flow by edge
solution$solution

# Include solution in edge dataframe
edgelist$flow <- solution$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 cost
E(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 colour
plot(g, edge.width = E(g)$flow)

希望它有趣/有用:)