Note that low-cost for-profit solvers are much faster and necessary for large problems.
(Inspired by Example 2.1.8 in C. Floudas Textbook)

library(lpSolveAPI, quietly = TRUE)
library(knitr, quietly = TRUE)

CASE A: Linear Cost Structure

Load Problem Data

Minimize total cost to satisfy demand at Demand Points given supply capacity as Supply Points
Cost is a linear function of quantity transported from Supply Point to Demand Point

df = read.csv("transprob.csv", row.names=1)

Label Demand and Supply Points

demand.points = paste("D",1:4,sep="")
supply.points = paste("S",1:6,sep="")

Ensure total supply equals total demand.
Not necessary for linear formulation, but helpful later.

stopifnot(sum(df[supply.points,'LS']) == sum(df['LD',demand.points]))

kable(df,caption="Cost Structure for Linear Problem and Supply/Demand Limits")
Cost Structure for Linear Problem and Supply/Demand Limits
D1 D2 D3 D4 LS
S1 300 270 460 800 8
S2 740 600 540 380 24
S3 300 490 380 760 20
S4 430 250 390 600 24
S5 210 830 470 680 16
S6 360 290 400 310 12
LD 29 41 13 21 104

Ingest raw data to R objects

cost.data = data.matrix(df[supply.points,demand.points])
n = length(supply.points)
m = length(demand.points)

Build Linear Program

lps.model = make.lp(0,n*m)
name.lp(lps.model, "Transportation Problem")

Model Constraints

r = 0

Supply Constraints

tmpl = seq(0,n*m-1,n)
for(i in 1:n){
  r = r + 1
  add.constraint(lps.model, rep(1,m), "<=", df[i,"LS"], tmpl+i)
}

Demand Constraints

tmpl = seq(1,n)
for(j in 1:m){
  r = r + 1
  add.constraint(lps.model, rep(1,n), ">=", df["LD",j], tmpl+n*(j-1))
}
dimnames(lps.model)[[1]] = c(supply.points, demand.points)
dimnames(lps.model)[[2]] = c(outer(supply.points,demand.points,paste,sep='_'))

Model Objective

set.objfn(lps.model, c(cost.data))

Optimal Solution

Call solve function from lpSolveAPI

print(solve(lps.model))
## [1] 0

If status code is 0, optimal solution was found.
Optimal Objective Value

print(get.objective(lps.model))
## [1] 31490

Optimal Flow from S to D Points

v = get.variables(lps.model)
constr = get.constraints(lps.model)

soln = rbind(
  cbind(
    matrix(round(v,3),n,m, dimnames=list(supply.points,demand.points)),
    TS=constr[1:n],
    LS=df[1:n,"LS"]
  ),
  TD=c(constr[n + (1:m)],0,0),
  LD=c(df["LD",1:m],0,0)
)
kable(soln,caption="Flow from S point to D Point, T: Total, L: Limit")
Flow from S point to D Point, T: Total, L: Limit
D1 D2 D3 D4 TS LS
S1 3 5 0 0 8 8
S2 0 0 3 21 24 24
S3 10 0 10 0 20 20
S4 0 24 0 0 24 24
S5 16 0 0 0 16 16
S6 0 12 0 0 12 12
TD 29 41 13 21 0 0
LD 29 41 13 21 0 0

Save LP

write.lp(lps.model,"transprob.lp")

Get Dual Solution

x  = get.dual.solution(lps.model)

Dual Prices for Supply

print(data.frame(dp=x[1+(1:n)],row.names=supply.points))
##      dp
## S1 -160
## S2    0
## S3 -160
## S4 -180
## S5 -250
## S6 -140

Dual Prices for Demand

print(data.frame(dp=x[1+n+(1:m)],row.names=demand.points))
##     dp
## D1 460
## D2 430
## D3 540
## D4 380

Cost Differential to Change Optimal Flow

print(matrix(x[-(1:(1+n+m))],n,m,dimnames=list(supply.points,demand.points)))
##     D1  D2  D3  D4
## S1   0   0  80 580
## S2 280 170   0   0
## S3   0 220   0 540
## S4 150   0  30 400
## S5   0 650 180 550
## S6  40   0   0  70

CASE B: Quadratic Costs - Negative Externality (e.g. congestion costs)

Introducing a Positive Definite Term in the Cost Function

Q = data.frame(
  D1 = c(7,12,13,7,4,17),
  D2 = c(4,9,12,9,10,9),
  D3 = c(6,14,8,16,21,8),
  D4 = c(8,7,4,8,13,4),
  row.names = supply.points
)
kable(Q)
D1 D2 D3 D4
S1 7 4 6 8
S2 12 9 14 7
S3 13 12 8 4
S4 7 9 16 8
S5 4 10 21 13
S6 17 9 8 4

Ingest data in objects agreeing with quadprog library

D-Matrix

Dmat = 2 * diag(c(data.matrix(Q)))

d-vector

dvec = c(cost.data)

A-Matrix

Amat = matrix(0,n+m,n*m)
r = 0

Supply Constraints

tmpl = seq(0,n*m-1,n)
for(i in 1:n){
  r = r + 1
  Amat[r,tmpl+i] = 1
}

Demand Constraints

tmpl = seq(1,n)
for(j in 1:m){
  r = r + 1
  Amat[r,tmpl+n*(j-1)] = 1
}
Amat = rbind(Amat,diag(n*m))

b-vector

bvec = c(unlist(df[1:n,"LS"]), unlist(df["LD",1:m]), rep(0, n*m))

Minimize Cost

qp = quadprog:::solve.QP(Dmat,-dvec,t(Amat),bvec)

kable(matrix(round(qp$solution,3),n,m, dimnames=list(supply.points,demand.points)),
      caption="Solution for Quadratic Cost Structure (Case B)")
Solution for Quadratic Cost Structure (Case B)
D1 D2 D3 D4
S1 0.000 8.000 0.000 0.000
S2 0.000 4.933 2.561 16.506
S3 8.464 3.791 7.744 0.000
S4 4.455 16.850 2.694 0.000
S5 16.000 0.000 0.000 0.000
S6 0.080 7.426 0.000 4.494

Compare to Solution Arising When Cost Structure is Linear (Case A)
Now, more supply destination options are utilized and optimal flow is not integral.

kable(matrix(v,n,m, dimnames=list(supply.points,demand.points)),
      caption="Linear Cost Solution (Case A)")
Linear Cost Solution (Case A)
D1 D2 D3 D4
S1 3 5 0 0
S2 0 0 3 21
S3 10 0 10 0
S4 0 24 0 0
S5 16 0 0 0
S6 0 12 0 0

CASE C: Quadratic Costs - Cost Discount (e.g. economies of scale)

Quadratic Term is not Positive Definite, so Quadratic Programming is not an option.

Quick Solution Strategy: Use Global Optimization.

  • This is not computationally efficient.
  • There are a better ways to do this, but using a generic solver is tempting when problem size is small.

Set Up Functions to help the DEoptim formulation

Couldn’t find an element-wise min function, so here is a home-made one

elementwise.min = function(x,y) ifelse(x>y,y,x)

Assign upper flow limits depending on the corresponding demand and supply constraints

ul = outer(df[supply.points,'LS'],
           as.numeric(df['LD',demand.points]),
           elementwise.min)

Make sure flow is not large enough to make any S/D transporation cost negative

ul = mapply(elementwise.min, as.data.frame(cost.data / Q), as.data.frame(ul))

Cost Evaluation Function

cost.eval = function(x){
  sum(x * (cost.data - x*Q) )
}

Decision Variables

Initial Solution: Chosen arbitrarily and leaving room for improvement

x.init = matrix(0,n,m)
x.init[1,1] = 2
x.init[1,2] = 6
x.init[2,1] = 3
x.init[2,4] = 21
x.init[3,1] = 19
x.init[3,3] = 1
x.init[4,1] = 2
x.init[4,2] = 22
x.init[5,1] = 3
x.init[5,2] = 1
x.init[5,3] = 12
x.init[6,2] = 12

Consider a perturbation flow, which if added to the initial solution flow can give any desired feasible flow.
The dimension of the perturbation is (n-1) X (m-1)
Here is a function to give a new solution, given a perturbation and the initial solution

perb2full = function(u,x0) {
  try(if(!is.matrix(x0)) stop("pmat is of class matrix"))
  d = dim(x0)
  out = matrix(0,d[1],d[2])
  out[-d[1],-d[2]] = matrix(u,d[1]-1,d[2]-1)
  out[-d[1], d[2]]= -apply(out[-d[1],-d[2]],1,sum)
  out[d[1],] = -apply(out[-d[1],],2,sum)
  x0+out
}

Problem Cost (adds penalty for negative flows)

probcost = function(u, penalty=cost.data+1000) {
  x = perb2full(u,x.init)
  cost.eval(x) - sum( penalty[x<0] * x[x<0] )
}

Minimize Cost

Cost of Initial Solution

probcost(rep(0,(n-1)*(m-1)))
## [1] 18622

Compute Optimum (Be Patient!)

set.seed(2020)
print(Sys.time())
## [1] "2024-12-11 15:57:58 EST"
start.time = Sys.time()
soln = DEoptim:::DEoptim(
  probcost,
  -c(x.init[-n,-m]),
  c(ul[-n,-m]-x.init[-n,-m]),
  control=DEoptim:::DEoptim.control(NP=150,itermax=5000,trace=FALSE,strategy=5,c=0.4),
  fnMap = round # requires integrality
)
print(Sys.time()-start.time)
## Time difference of 16.50448 mins

Min Cost Flow Obtained Using the Global Optimization Heuristic

Near-Optimal Cost

print(soln$optim$bestval)
## [1] 15639

Near-Optimal Flow

round(perb2full(soln$optim$bestmem,x.init),2)
##      [,1] [,2] [,3] [,4]
## [1,]    6    2    0    0
## [2,]    0    3    0   21
## [3,]   20    0    0    0
## [4,]    0   24    0    0
## [5,]    3    0   13    0
## [6,]    0   12    0    0