C:\JUAN\POST\MARKDOWN
Linear Programming Assignment 4 - PFCM
2020-01-21
Propose a maximum-flow problem involving 11 nodes, defining the network with its corresponding arcs and capacities.
In the year 2230 AD, complex life is found for the first time on planets not very far from Earth. From that moment on, the profession of astrobiologist gains momentum, since those planets are rich in biodiversity with great potential for exploitation as a resource, so young researchers venture onto the wild exoplanets in search of a discovery that will make them rich and famous when they return home.
On one of these research and adventure trips, a group of astrobiologists funded by the USC discover a living being similar to a terrestrial plant that captures molecules dispersed in the air through an “antenna” and carries out a complex digestive process through nine “gastric nodules”, which ultimately send extremely nutritious and purified complex juices to the roots.
The USC researchers have experimentally verified the nutritional value of these juices and estimate that they could become rich by selling food supplements if they were able to optimise the process, since this being seems to use some of the fluids endogenously to produce other structures, such as defensive spines, petals or seeds, which would not be of commercial interest.
After a long deliberation, one of those biologists - who had taken the USC course in Linear and Integer Programming - very shrewdly proposed that an artificial-selection process could be carried out to improve the species (provisionally named Linae programatta) by selecting only the specimens that send the greatest amount of this juice directly to the root, where it will be easy to collect with a simple vial. It would be very useful to predict which branches are chiefly responsible for the net flow reaching the root and artificially force the less important ones to disappear. This would save hundreds of failed trials, time and money in the crossbreeding process.
“In linear-programming terms…” thought the astro-biostatistician, “…this is still a Minimum-Cost Network Flow problem. If we obtain the graph of this being's vascularisation from the data we have on its vascular system, we will know how to optimise the flows, which edges are unnecessary, and we will be able to select specimens that come closer to our desired characteristics until, after a few iterations, the production of juice that can be extracted directly through the root reaches its optimum.”
“We are going to be rich, and all thanks to Linear Programming!” they all shouted.
“And me most of all! I am the only one who knows how to program this problem in R!” shouted the astro-biostatistician, glad to have completed his Master's in Statistics in his younger days.
And this is how he solved the problem of optimising the flow of nutritious gastric juices towards the root in Linae programatta:
[Note: capacity measurements are expressed in millilitres of gastric juice per hour (ml/h)]
Corresponding graph: defining the nodes, edges and their associated bounds and costs in R.
library('igraph')##
## Attaching package: 'igraph'## The following objects are masked from 'package:stats':
##
## decompose, spectrum## The following object is masked from 'package:base':
##
## unionnodes <- c("A", "1", "2", "3", "4", "5", "6", "7", "8", "9", "R")
# A=Antena Filtradora
# R=Raiz
# 1-9: Nódulos Gástricos
# coordenadas en X, Y y color de los nodos
x <- c(3,1,5,2,4,5,3,1,2,3,3)
y <- c(1,2,2,3,3,3,4,4,5,6,7)
color <- c("red", rep("lightblue", 9), "gold")
nodos <- data.frame(nodes, x ,y, color)
# direccion aristas
from <- c('A','A','A','A','1','2','3','4','3','4','6','6','8', '8','9')
to <- c('1','3','4','2','3','4','4','5','6','6','9','8','7', '9','R')
capacity <- c( 5, 1, 1, 5, 5, 5, 3, 6, 5, 10, 2, 6 , 2, 15, 20)
aristas <- data.frame(from, to, capacity)
grafo <- graph_from_data_frame(vertices = nodos, d = aristas, directed = T)
E(grafo)$label <- c(
"(0,5,0)",
"(0,1,0)",
"(0,1,0)",
"(0,5,0)",
"(0,5,0)",
"(0,5,0)",
"(0,3,0)",
"(0,6,0)",
"(0,5,0)",
"(0,10,0)",
"(0,2,0)",
"(0,6,0)",
"(0,2,0)",
"(0,15,0)",
"(0,20,0)")
tkplot(grafo)## [1] 1# se usa tkPlot grafo para poder hacer edición de la apariencia a posteriori,
# se incluye a continuación la imagen corregida
“Graph of the gastric system of Linae programatta Sp.. The antenna (A) is shown in red, the root (R) in gold and the gastric nodules (1-9) in blue”
- With the following associated arcs and flow capacities:
aristas## from to capacity
## 1 A 1 5
## 2 A 3 1
## 3 A 4 1
## 4 A 2 5
## 5 1 3 5
## 6 2 4 5
## 7 3 4 3
## 8 4 5 6
## 9 3 6 5
## 10 4 6 10
## 11 6 9 2
## 12 6 8 6
## 13 8 7 2
## 14 8 9 15
## 15 9 R 20Write the linear-programming model for the proposed problem, explaining the meaning of the variables, the objective function and the constraints.
In mathematical language, the proposed linear-programming model can be defined as follows:
Let \(G=(N,M)\) be a directed graph with nodes \(N=({A,1,...,9,R})\) and arcs \(M=(a_{1},...,a_{15})\). This graph can be represented as an incidence matrix \(\overline{B}_{n* m}\) whose coordinates are assigned a \(1\) si \(i\) is the initial node of arc \(a_{k}\), \(-1\) si \(i\) is the terminal node of \(a_{k}\) and \(0\) otherwise. This matrix is created as an R object in the following section; nevertheless, I include it here as a matrix for visualisation purposes, with the edges on the X axis and the nodes on the Y axis:
In addition, each arc \(k\) has three associated parameters: a lower bound (minimum flow through arc \(k\)), an upper bound (maximum flow through arc \(k\)) and coste (positivo para costes de flujo y negativo para beneficios).
\[ \begin{pmatrix} Cota Inferior & l_{k} \geq 0 \\ Cota Superior & u_{k} \geq 0\\ Coste & c_{k}\\ \end{pmatrix} \]
Con la incidence matrix podemos plantear el problema correctamente como una función a minimizar (maximizar en nuestro caso pero será cuestión de hacer la inversa del resultado y listo) y una serie de restricciones.
- First, define the variables, which will be the different edges over which we calculate the flow:
\[ \begin{pmatrix} Variable & Arista \\ i_{1} & 1ª \\ i_{2} & 2ª \\ (...) & (...) \\ i_{15} & 15ª \\ i_{A} & Auxiliar \\ \end{pmatrix} \]
- And now we define the objective function and its constraints:
\(\quad\)\(\quad\) Min \(\quad\)\(\quad\) \(-i_{A}\)
Subject to:
\(\quad\)\(\quad\)\(\quad\)\(\quad\) \(i_{1}\)\(+i_{2}\)\(+i_{4}\)\(+i_{5}-i_{16}\)\(=\)\(0\)
\(\quad\)\(\quad\)\(\quad\)\(\quad\) \(-i_{1}\)\(+i_{3}\)\(=\)\(0\)
\(\quad\)\(\quad\)\(\quad\)\(\quad\) \(-i_{2}\)\(+i_{6}\)\(=\)\(0\)
\(\quad\)\(\quad\)\(\quad\)\(\quad\) \(-i_{3}\)\(-i_{4}\)\(+i_{7}\)\(+i_{8}\)\(=\)\(0\)
\(\quad\)\(\quad\)\(\quad\)\(\quad\) \(-i_{5}\)\(-i_{6}\)\(-i_{7}\)\(+i_{9}+i_{10}\)\(=\)\(0\)
\(\quad\)\(\quad\)\(\quad\)\(\quad\) \(-i_{10}\)\(=\)\(0\)
\(\quad\)\(\quad\)\(\quad\)\(\quad\) \(-i_{8}\)\(-i_{9}\)\(+i_{11}\)\(+i_{12}\)\(=\)\(0\)
\(\quad\)\(\quad\)\(\quad\)\(\quad\) \(-i_{13}\)\(=\)\(0\)
\(\quad\)\(\quad\)\(\quad\)\(\quad\) \(-i_{12}\)\(+i_{13}\)\(+i_{14}\)\(=\)\(0\)
\(\quad\)\(\quad\)\(\quad\)\(\quad\) \(-i_{11}\)\(-i_{14}\)\(=\)\(0\)
\(\quad\)\(\quad\)\(\quad\)\(\quad\) \(-i_{15}\)\(+i_{16}\)\(=\)\(0\)
\(\\\)
Solve the problem with R.
- Ya tenemos el problema de optimización PFCM diseñado gráficamente (como se puede visualizar en la Figura 1) y listo para ser resuelto con R. Procedemos entonces a resolver el problema mediante el uso de lpSolverAPI definiendo antes la incidence matrix, función objetivo a maximizar, restricciones… y demás.
library(lpSolveAPI)
nodos <- 11
arcos <- 16 # 15 más el arco auxiliar
# incidence matrix
B <- matrix(0, nrow = nodos, ncol = arcos)
# nodos iniciales
B[1, 1] <- 1
B[1, 2] <- 1
B[1, 4] <- 1
B[1, 5] <- 1
B[2, 3] <- 1
B[3, 6] <- 1
B[4, 7] <- 1
B[4, 8] <- 1
B[5, 9] <- 1
B[5, 10] <- 1
B[7, 11] <- 1
B[7, 12] <- 1
B[9, 13] <- 1
B[9, 14] <- 1
B[10, 15] <- 1
B[11, 16] <- 1
# nodos terminales
B[2, 1] <- -1
B[3, 2] <- -1
B[4, 3] <- -1
B[4, 4] <- -1
B[5, 5:7] <- -1
B[6, 10] <- -1
B[7, 8:9] <- -1
B[9, 12] <- -1
B[8, 13] <- -1
B[10, 11] <- -1
B[10, 14] <- -1
B[11, 15] <- -1
B[1,16] <- -1
B## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13]
## [1,] 1 1 0 1 1 0 0 0 0 0 0 0 0
## [2,] -1 0 1 0 0 0 0 0 0 0 0 0 0
## [3,] 0 -1 0 0 0 1 0 0 0 0 0 0 0
## [4,] 0 0 -1 -1 0 0 1 1 0 0 0 0 0
## [5,] 0 0 0 0 -1 -1 -1 0 1 1 0 0 0
## [6,] 0 0 0 0 0 0 0 0 0 -1 0 0 0
## [7,] 0 0 0 0 0 0 0 -1 -1 0 1 1 0
## [8,] 0 0 0 0 0 0 0 0 0 0 0 0 -1
## [9,] 0 0 0 0 0 0 0 0 0 0 0 -1 1
## [10,] 0 0 0 0 0 0 0 0 0 0 -1 0 0
## [11,] 0 0 0 0 0 0 0 0 0 0 0 0 0
## [,14] [,15] [,16]
## [1,] 0 0 -1
## [2,] 0 0 0
## [3,] 0 0 0
## [4,] 0 0 0
## [5,] 0 0 0
## [6,] 0 0 0
## [7,] 0 0 0
## [8,] 0 0 0
## [9,] 1 0 0
## [10,] -1 1 0
## [11,] 0 -1 1write.table(B, "./B.txt", sep=";")
# costes
costes <- rep(0, arcos)
costes[16] <- c(-1)
costes## [1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1# cotas inferiores
L <- rep(0, arcos)
# cotas superiores
U <- rep(1e30, arcos)
U[1:(arcos-1)] <- c(5,5,5,1,1,5,3,5,10,6,2,6,2,15,20)
U## [1] 5.0e+00 5.0e+00 5.0e+00 1.0e+00 1.0e+00 5.0e+00 3.0e+00 5.0e+00 1.0e+01
## [10] 6.0e+00 2.0e+00 6.0e+00 2.0e+00 1.5e+01 2.0e+01 1.0e+30# resolvemos el problema asociado
n_restricciones <- nrow(B)
n_variables <- ncol(B)
PFCM <- make.lp(nrow = n_restricciones, ncol = n_variables)
# matriz de restricciones
for (i in 1:n_restricciones){
set.row(PFCM, i, B[i, ])
}
set.objfn(PFCM, costes)
b <- rep(0, n_restricciones)
set.rhs(PFCM, b)
tipores <- rep("=", n_restricciones)
set.constr.type(PFCM, tipores)
set.bounds(PFCM, upper = U)
set.bounds(PFCM, lower = L)
PFCM## Model name:
## a linear program with 16 decision variables and 11 constraintsval <- solve(PFCM)
val## [1] 0obj<-get.objective(PFCM)
obj## [1] -8vars<-get.variables(PFCM)
vars## [1] 4 2 4 1 1 2 0 5 3 0 2 6 0 6 8 8The lpSolverAPI output tells us that there is an optimal solution to the stated problem.
The optimal solution entails a flow of 8 ml/h of gastric juice through the root, as indicated by the value of the function being optimised.
Finally, lpSolveAPI tells us which edges (connections between nodules) would have zero flow at the optimum. These are 7th edge (nodos 3-4), 10ª (4-5) and 13ª (8-7).
# cambiamos el color de los vértices (azul) que maximizan el flujo
vcol <- rep("gray", 14)
vars.n <- vars[-15]
vcol[vars.n != 0] <- "green"
E(grafo)$color = vcol
tkplot(grafo)## [1] 2
“Zoo-graph to be optimised through artificial selection of Linae programatta Sp. propuesto por R-lpSolveAPI. Con las aristas a conservar para optimizar flujo máximo en verde and las deshechables en gris.”
The resulting graph indicates that the edges with no flow towards the root - and which are therefore targets for the astrobotanists' artificial-selection efforts - are:
The 7th edge, which joined two intermediate nodes and probably had a structural function, will not be needed in the greenhouse, making it a perfect candidate for flow optimisation.
The 10th and 13th edges supported defensive spines that it will not need in a predator-free environment either, making them two more perfect candidates for flow optimisation.
The resulting example graph that astrobotanists should follow in their crops to optimise juice extraction through artificial selection can be seen in Figure 2.
Solve the problem with AMPL.
The astrobiologists have an internet connection for a couple of hours each month courtesy of the Xunta, so they use the time not only to program the PFCM problem in R, but also to do so in AMPL and solve it in the cloud with Gurobi. To do this, they submit the following text documents:
- flujo.dat with the data:
set INTER := A 1 2 3 4 5 6 7 8 9 R ;
param entr := A ;
param sali := R ;
param: CAMI: capa :=
A 1 5, A 2 5, A 3 1, A 4 1,
1 3 5,
2 4 5,
3 4 3, 3 6 5,
4 5 6, 4 6 10,
6 9 2, 6 8 6,
8 7 2, 8 9 15,
9 R 20;
- flujo.mod with the proposed model:
set INTER; # intersecciones
param entr symbolic in INTER; # entrada a la red
param sali symbolic in INTER, <> entr; # salida de la red
set CAMI within (INTER diff {sali}) cross (INTER diff {entr});
param capa {CAMI} >= 0; # capacidades de caminos
var Traf {(i,j) in CAMI} integer, <= capa[i,j], >=0; # tráfico
maximize Entrada_Traf: sum {(entr,j) in CAMI} Traf[entr,j];
subject to Equilibrio {k in INTER diff {entr,sali}}: sum {(i,k) in CAMI} Traf[i,k] = sum {(k,j) in CAMI} Traf[k,j];
- flujo.run with the commands and desired outputs:
#Ejemplo del problema del flujo máximo
solve;
display Traf, Entrada_Traf;
\(\\\)
- After waiting a few hours, they obtain the following flow-optimisation result:
Gurobi 8.1.0: optimal solution; objective 8 Traf := 1 3 2 2 4 5 3 4 2 3 6 0 4 5 0 4 6 8 6 8 6 6 9 2 8 7 0 8 9 6 9 R 8 A 1 2 A 2 5 A 3 0 A 4 1 ;
Entrada_Traf = 8
– The optimal-flow result is the same; it should reach \(8ml/h\), on which both models agree.
– Nevertheless, AMPL proposes certain new edges that could potentially be removed during the optimisation process:
\((3 - 6)\): Esta arista no había sido propuesta previamente and aunque parece ser central en el desarrollo de Linae programatta given its anatomical position, it also appears to be a bottleneck in the flow of gastric juices, so an attempt should be made to select varieties with greater capacity in this edge or, if that is not possible, perhaps eliminate it.
\((4 - 5)\): This edge had already been proposed. Since it is a spine, it is merely a useless sink for gastric juices. It must be optimised.
\((8 - 7)\): This edge had already been proposed. Since it is a spine, it is merely a useless sink for gastric juices. It must be optimised.
\((A - 3)\): This edge is also new and, although it appears central to the anatomy of Linae programatta, closer analysis indicates that it has very little flow, probably only to support petals and sepals and form seeds. This edge could also be optimised through crossbreeding.
“The graph of Linae programatta in an optimal-flow format according to AMPL would look something like the following, with the edges that optimise flow shown in green.”
CONCLUSION:
Thanks to formulating this as a Linear Programming problem for optimising a flow, the astrobiologists now have a much clearer understanding of the vascular anatomy and physiology of Linae programatta. Thanks to this, they will be able to carry out crosses and artificially select the offspring that come closest to meeting these standards established with R and/or AMPL (see Figure 4), until the species is stabilised with beneficial characteristics that allow it to be exploited as a commercial resource.
(…)
Meses después, los invernaderos de campaña y los entrecruzamientos han dado su fruto and los astro-biólogos -liderados por el astro-bio-estadista - vuelven a casa con variedades seleccionadas de Linae programatta ready to scale up production and, if they become as rich as they hope, establish a chair in spatial optimisation at the USC Faculty of Mathematics.
“Graph of the vascular system of the selected variety of Linae programatta. The edges on which both R and AMPL agreed as potential optimisation targets have been removed, along with intermediate edges with little capacity.”