最小成本流 - R 中的网络优化
Minimum Cost Flow - network optimization in R
我正在尝试在 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)
这是一个具有方向边的网络，一个 "source node" (1) 和一个 "sink node" (9)。每条边都有一个 "capacity"（这里通常表示为 99 表示无限制）和一个 "cost"（一个单位流经该边的成本）。我想找到流的整数向量 (x, length = 9)，它可以在通过网络传输预定义流时最小化成本（假设 50 个单位，从节点 1 到节点 9）。
免责声明：this post 提出了类似的问题，但没有得到令人满意的答案，而且已经过时（2012 年）。
我一直在寻找一个功能，但没有成功。原始函数调用另一个：
res <- .Call("R_igraph_maxflow", graph, source - 1, target -
1, capacity, PACKAGE = "igraph")
而且我不知道怎么处理。
目前我反转成本路径值以便在相反方向使用相同的函数:
E2 <- E # create another table
E2[, 3] <- max(E2[, 3]) + 1 - E2[, 3] # invert values
E2
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 1
g2 <- graph_from_data_frame(as.data.frame(E2)) # create 2nd graph
# Get maximum flow
m1 <- 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 2
m2$partition2 # Route on minimal cost
+ 3/6 vertices, named:
[1] 3 4 2
我在纸上画了图，我的代码与手动解决方案一致
正如我在评论
中提到的，应该使用真实已知值测试此方法
如果有兴趣，这里是我最终解决这个问题的方法。我使用带有 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)
希望它有趣/有用 :)
我正在尝试在 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)
这是一个具有方向边的网络，一个 "source node" (1) 和一个 "sink node" (9)。每条边都有一个 "capacity"（这里通常表示为 99 表示无限制）和一个 "cost"（一个单位流经该边的成本）。我想找到流的整数向量 (x, length = 9)，它可以在通过网络传输预定义流时最小化成本（假设 50 个单位，从节点 1 到节点 9）。
免责声明：this post 提出了类似的问题，但没有得到令人满意的答案，而且已经过时（2012 年）。
我一直在寻找一个功能，但没有成功。原始函数调用另一个：
res <- .Call("R_igraph_maxflow", graph, source - 1, target -
1, capacity, PACKAGE = "igraph")
而且我不知道怎么处理。
目前我反转成本路径值以便在相反方向使用相同的函数:
E2 <- E # create another table
E2[, 3] <- max(E2[, 3]) + 1 - E2[, 3] # invert values
E2
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 1
g2 <- graph_from_data_frame(as.data.frame(E2)) # create 2nd graph
# Get maximum flow
m1 <- 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 2
m2$partition2 # Route on minimal cost
+ 3/6 vertices, named:
[1] 3 4 2
我在纸上画了图，我的代码与手动解决方案一致 正如我在评论
中提到的，应该使用真实已知值测试此方法如果有兴趣，这里是我最终解决这个问题的方法。我使用带有 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)
希望它有趣/有用 :)