Constrained nonlinear network flow optimization problem
In this post, I am going to reuse the previously outlined network in order to solve an optimization problem whereby the cost function is based upon end-to-end delay of the traffic flow we intend to allocate. Such end-to-end delay will be computed based on M/M/1 queue modeling. In other words I will assume the inter-arrival and service time driven by a Poisson process. Queue theory tells us that the M/M/1 queue sojourn time is equal to:
\({M/M/1\:sojourn\:time} = \dfrac{1}{\mu - \lambda}\)
Such sojourn time takes into account the queueing time and the service time at nodes. The transportation delay introduced by links is modeled as fixed delay implied by the link lengths. The node delay is modeled by taking advantage of the M/M/1 queue delay computation whereby the service speed is determined by the outgoing link bandwidth.
Our sample network is the same of previous post, a network made of four nodes connected by five links.
.
Comments about link labels are given later.
Known variables are:
\(\delta_{i,j}\) the fixed delay through link connecting nodes i and j whereby a link between node i and j exists.
Unknown variables are:
\(x_{1,2}\) the traffic amount delivered through link connecting nodes 1 and 2
\(x_{1,3}\) the traffic amount delivered through link connecting nodes 1 and 3
\(x_{2,3}\) the traffic amount delivered through link connecting nodes 2 and 3
\(x_{3,2}\) the traffic amount delivered through link connecting nodes 3 and 2
\(x_{2,4}\) the traffic amount delivered through link connecting nodes 2 and 4
\(x_{3,4}\) the traffic amount delivered through link connecting nodes 3 and 4
In the order, above listed variables are the unknowns of our problem.
The delay implied in flowing through a specific link connecting node i with node j can be modeled as (per each delay units):
\(S(x_{i,j}) = H(x_{i,j})\)
where H is the Heaveside function. It means that if a flow \(x_{i,j}\) is travelling through link (i,j), some fixed delay is implied. Heaveside R code is herein below shown.
Heaveside <- function(x) { as.numeric(x > 0) }If the optimization problem solver routine would require to compute gradients (numerically, in the case) the Heaveside function can be replaced by a sigmoid function, such as for example:
\(S(x_{i,j}) = \dfrac{2}{1+e^{-10*x}} - 1\)
Sigmoid <- function(x) { 2/(1+exp(-10*x)) - 1}Our objective function to minimize is:
\(min_{\textbf{x}} F(\textbf{x}) = min_{\textbf{x}} (\delta_{1,2}*S(x_{1,2}) + \dfrac{1}{c_{1,2}-x_{1,2}} + \delta_{1,3}*S(x_{1,3}) + \dfrac{1}{c_{1,3}-x_{1,3}} + \delta_{2,3}*S(x_{2,3}) + \dfrac{1}{c_{2,3}-x_{2,3}} + \\ + \delta_{3,2}*S(x_{3,2}) + \dfrac{1}{c_{3,2}-x_{3,2}} + \delta_{2,4}*S(x_{2,4}) + \dfrac{1}{c_{2,4}-x_{2,4}} + \delta_{3,4}*S(x_{3,4}) + \dfrac{1}{c_{3,4}-x_{3,4}})\)
where:
\(\textbf{x} = (x_{1,2}, x_{1,3}, x_{2,3}, x_{3,2}, x_{2,4}, x_{3,4})\)
is the vector of unknowns. As you can see, the cost function is a sum of fixed and M/M/1 queue delay terms which takes into account all potential paths that can bring traffic from node 1 to node 4.
A first set of constraints is determined by the traffic demand at the source and destination nodes and the absence of traffic losses throughout the network.
\(x_{1,2} + x_{1,3} = D\)
\(x_{1,2} - x_{2,3} - x_{2,4} + x_{3,2} = 0\)
\(x_{1,3} + x_{2,3} - x_{3,2} + x_{3,4} = 0\)
\(x_{2,4} + x_{3,4} = D\)
where D is the traffic demand.
Further constraints are determined by the links capacity limit and the non-negative traffic flows amount value constraints, as herein shown.
\(x_{1,2} <= c_{1,2}\), \(x_{1,3} <= c_{1,3}\), \(x_{2,3} <= c_{2,3}\), \(x_{2,4} <= c_{2,4}\), \(x_{3,2} <= c_{2,3}\), \(x_{3,4} <= c_{3,4}\)
\(x_{1,2} >= 0\), \(x_{1,3} >= 0\), \(x_{2,3} >= 0\), \(x_{2,4} >= 0\), \(x_{3,2} >= 0\), \(x_{3,4} >= 0\)
To solve numerically this problem, I will take advantages of two R packages, igraph and Rsolnp.
The igraph package basically allows to create graphs of various types, associate attributes to vertexes and edges, plot them and to compute properties as defined the mathematical graph theory.
The Rsolnp package basically allows to define and solve non-linear constrained optimization problems, with both equality and inequality constraints.
The igraph object which represents the network as above displayes is herein defined.
suppressPackageStartupMessages(library(igraph))
g <- graph.formula(1-2, 1-3, 2-3, 2-4, 3-4)
E(g)$capacity <- c(10, 9, 8, 10, 10)
E(g)$delay <- c(10, 4, 6, 15, 12).
I defined the links capacities and delays and assigned to the igraph object as edge custom properties. Link delay values are used as link label in the picture. Links capacities will be involved in constraints definition, while the link delays will be involved in the functional cost definition.
As further input to our problem, the flow request is defined as below.
request <- list("node.source" = 1, "node.destination" = 4, "demand" = 15)The above request asks for determining a path able to carry traffic from node 1 to node 4. The traffic demand may also be forwarded by splitting it in two or more sub-flows crossing the network throughout different paths starting from node 1 and reaching node 4. That will be for sure implied by any demand whose value is strictly greater than the links capacities at some stage of the forwarding process. Considering the above demand value of 15 and the links capacities which are not strictly greater than 10, it is the case.
The objective function results:
\(min_{\textbf{x}} F(\textbf{x}) = min_{\textbf{x}} (10*S(x_{1,2}) + \dfrac{1}{c_{1,2}-x_{1,2}} + 4*S(x_{1,3}) + \dfrac{1}{c_{1,3}-x_{1,3}} + 6*S(x_{2,3}) + \dfrac{1}{c_{2,3}-x_{2,3}} + \\ 6*S(x_{3,2}) + \dfrac{1}{c_{3,2}-x_{3,2}} + 15*S(x_{2,4}) + \dfrac{1}{c_{2,4}-x_{2,4}} + 12*S(x_{3,4}) + \dfrac{1}{c_{3,4}-x_{3,4}})\)
as \(\delta_{3,2}\) is equal to \(\delta_{2,3}\)
The equality constraints are:
\(x_{1,2} + x_{1,3} = 15\)
\(x_{1,2} - x_{2,3} - x_{2,4} + x_{3,2} = 0\)
\(x_{1,3} + x_{2,3} - x_{3,2} + x_{3,4} = 0\)
\(x_{2,4} + x_{3,4} = 15\)
The inequality constraints are:
\(x_{1,2} <= 10\), \(x_{1,3} <= 9\), \(x_{2,3} <= 8\), \(x_{2,4} <= 10\), \(x_{3,2} <= 8\), \(x_{3,4} <= 10\)
\(x_{1,2} >= 0\), \(x_{1,3} >= 0\), \(x_{2,3} >= 0\), \(x_{2,4} >= 0\), \(x_{3,2} >= 0\), \(x_{3,4} >= 0\)
Following, R code to find the solution where basically:
the Rsolnp package is loaded
the Heaveside function is choosen for modeling fixed delays
the M/M/1 queue delay function is defined
the number of unknows is determined by the graph edges set
the function to be optimized, i.e. the source-to-node flow delay function is defined (fn_delay)
equality constraints function is defined (eqn_flow_graph)
inequality constraintes function is defined (eqn_flow_lim)
initial unknown variable state is defined (x0)
equality constraints right term is defined (eqB_rt)
inequality constraints left term is defined (ineqLB_lt)
inequality constraints right term is defined (ineqLB_rt)
suppressPackageStartupMessages(library(Rsolnp))
# modeling fixed delays using Heaveside function
S <- Heaveside
# M/M/1 queue delay function, c = link capacity, x = incoming flow bandwidth
MM1 <- function(c, x) { 1/(c-x) }
# unknowns of out optimization problem is equal to unknown flows
n_unknowns <- length(E(g))+1
# the following mapping between flow unknowns and x[] variables
# x12 = x[1]
# x13 = x[2]
# x23 = x[3]
# x32 = x[4]
# x24 = x[5]
# x34 = x[6]
#
# cost function
fn_delay <- function(x) {
E(g)$delay[1]*S(x[1]) +
E(g)$delay[2]*S(x[2]) +
E(g)$delay[3]*S(x[3]) +
E(g)$delay[3]*S(x[4]) +
E(g)$delay[4]*S(x[5]) +
E(g)$delay[5]*S(x[6]) +
MM1(E(g)$capacity[1], x[1]) +
MM1(E(g)$capacity[2], x[2]) +
MM1(E(g)$capacity[3], x[3]) +
MM1(E(g)$capacity[3], x[4]) +
MM1(E(g)$capacity[4], x[5]) +
MM1(E(g)$capacity[5], x[6])
}
# the h[4] constraint would be reported by solnp as redundant
# because it can be obtained by summing up constraints #2 and #3 while
# taking advantage of constraint#1
#
# equality constraint function
eqn_flow_graph <- function(x) {
h = rep(NA, 3)
h[1] = x[1]+x[2]-request$demand
h[2] = x[1]-x[3]-x[5]+x[4]
h[3] = x[2]+x[3]-x[4]-x[6]
# h[4] = x[5]+x[6]-request$demand ### redundant constraint
h
}
# inequality constraints are applied to all unknowns
eqn_flow_lim <- function(x) {
h = rep(NA, n_unknowns)
h[1] = x[1]
h[2] = x[2]
h[3] = x[3]
h[4] = x[4]
h[5] = x[5]
h[6] = x[6]
h
}
# initial state
x0 = rep(0, n_unknowns)
# equality constraint equation right term
eqB_rt <- rep(0, 3)
# inequality constraint left term
ineqLB_lt <- rep(0, n_unknowns)
# inequality constraint right term
ineqUB_rt <- c(E(g)$capacity[1], E(g)$capacity[2], E(g)$capacity[3],
E(g)$capacity[3], E(g)$capacity[4], E(g)$capacity[5])
# optimal solution
opt_delay <- solnp(pars = x0, fun = fn_delay,
eqfun = eqn_flow_graph, eqB = eqB_rt,
ineqfun = eqn_flow_lim, ineqLB = ineqLB_lt,
ineqUB = ineqUB_rt
)##
## Iter: 1 fn: 55.0574 Pars: 7.965424303 7.034575705 0.397922846 0.000004274 7.567505731 7.432494277
## Iter: 2 fn: 55.0574 Pars: 7.9654145682 7.0345854318 0.3979205355 0.0000006069 7.5674946396 7.4325053604
## solnp--> Completed in 2 iterationsIt takes only two interactions to find the optimal solution. The 0- discontinuity of the Heaveside function seems not to affect the solnp routine. Following, the optimal flow values are collected in a data frame.
unknowns_label <- c("X12", "X13", "X23", "X32", "X24", "X34")
opt_delay$pars## [1] 7.965415e+00 7.034585e+00 3.979205e-01 6.068992e-07 7.567495e+00
## [6] 7.432505e+00result <- data.frame("flows" = unknowns_label, "values" = opt_delay$pars)
result## flows values
## 1 X12 7.965415e+00
## 2 X13 7.034585e+00
## 3 X23 3.979205e-01
## 4 X32 6.068992e-07
## 5 X24 7.567495e+00
## 6 X34 7.432505e+00Further, the optimal delay value is:
opt_delay$values[length(opt_delay$values)]## [1] 55.05743Two questions arise:
why does the optimal solution prescribe flows to be exchanged both directions between nodes 2 and 3 ?
would it make any difference using the sigmoid function to model fixed delays ?
First issue appears to be some sub-optimal choice made by the solver and the second question suggests a comparison with the sigmoid function used for modeling fixed delays.
# modeling fixed delays using Heaveside function
S <- Sigmoid
# computing again the optimal solution with all the rest as unchanged
opt_delay <- solnp(pars = x0, fun = fn_delay,
eqfun = eqn_flow_graph, eqB = eqB_rt,
ineqfun = eqn_flow_lim, ineqLB = ineqLB_lt,
ineqUB = ineqUB_rt
)##
## Iter: 1 fn: 43.0715 Pars: 7.828435801256 7.171564206640 0.000000004262 0.000000004250 7.828435801242 7.171564206653
## Iter: 2 fn: 43.0715 Pars: 7.828435493636 7.171564506364 0.000000001112 0.000000001085 7.828435493609 7.171564506391
## solnp--> Completed in 2 iterations# reporting all new results
unknowns_label <- c("X12", "X13", "X23", "X32", "X24", "X34")
opt_delay$pars## [1] 7.828435e+00 7.171565e+00 1.111780e-09 1.085081e-09 7.828435e+00
## [6] 7.171565e+00result <- data.frame("flows" = unknowns_label, "values" = opt_delay$pars)
result## flows values
## 1 X12 7.828435e+00
## 2 X13 7.171565e+00
## 3 X23 1.111780e-09
## 4 X32 1.085081e-09
## 5 X24 7.828435e+00
## 6 X34 7.171565e+00opt_delay$values[length(opt_delay$values)]## [1] 43.07146The optimal solution so found makes more sense, and flows between nodes 2 and 3 can be assumed as equal to zero. Also notice that the cost function at optimal point has lower value than previous solution has. You may improve precision by adding to solnp() call the relative tolerance on feasibility and optimality tol parameter within the control list.
In this post, I have shown the optimal solution to a non-linear optimization problem with both equality and non-equality constraints that arises from a single commodity network problem based on flow delay cost functional.
