##################################### # This file is part of the # # Xpress-R interface examples # # # # (c) 2022-2024 Fair Isaac Corporation # ##################################### #' --- #' title: "Water Supply Management" #' author: Y. Gu #' date: Jul. 2021 #' --- #' #' #' ## ----setup, include=FALSE----------------------------------------------------- knitr::opts_chunk$set(echo = TRUE) knitr::opts_chunk$set(results = "hold") knitr::opts_chunk$set(warning = FALSE, message = FALSE) #' #' #' ## Brief Introduction To The Problem #' #' [This is a conversion of the Mosel example 'Water Supply Management'](https://www.fico.com/fico-xpress-optimization/docs/latest/examples/mosel/ApplBook/J_Service/j1water.mos). #' A brief introduction to this problem is given below, and to see the full #' mathematical modeling of this problem you may refer to section 15.1, page 229 of the #' book 'Applications of optimization with Xpress'. #' #' This example is a classical maximum flow problem. There are 2 reservoirs, 3 cities #' and 5 pumping stations. The three cities are supplied from the two reservoirs. #' All the 10 nodes are connected by a network of pipes, and we want to determine the #' maximum flow in this network to see whether the demands can be satisfied entirely. #' #' As a typical maximum flow problem, firstly we create two artificial nodes SOURCE and #' SINK, where SOURCE is connected to the two reservoirs by two arcs with capacities #' corresponding to the availability of water from the two reservoirs, and SINK is the #' node to which the three cities are connected by three arcs with capacities corresponding #' to the cities' requirement for water. #' #' Then, we define variables 'flow' for each arc and the 'flow' of each arc should not #' exceed its capacity. Note that we set the variable type as continuous, because it is #' possible to let non-integer flows pass through the arcs. But since all capacities here #' are integer, the simplex algorithm will automatically find integer solution values in #' the optimal solution to the linear problem. #' #' Our objective is to maximize the total flow, which is the sum of flows into SINK or #' the total flows out of SOURCE. An important constraint needed to be satisfied in #' network flow problems is the so-called flow conservation for all nodes except for #' SINK and SOURCE. That is, for the intermediate nodes, the the total flow arriving at #' any node also has to leave this node. #' #' After solving this problem, we can see the optimum total flow is 52, which is 1 less #' that the total demand from the three cities. So, the existing network will not be able #' to satisfy the demands. The method used here may be applied to other liquids, and also #' to road or telecommunications networks. #' #' #' For this example, we need packages 'xpress', 'dplyr' and 'igraph'. Besides, we use the #' function `pretty_name` to give the variables and constraints concise names. #' ## ----Load The Packages And The Function To Give Names------------------------- library(xpress) library(dplyr) library(igraph) pretty_name <- function(prefix, y) { "%s_%s" %>% sprintf( prefix, paste(lapply(names(y), function(name) paste(name, y[name], sep="_")), collapse = "_") ) } #' #' #' Create a new empty problem, set the objective sense as maximization and give the #' problem a suitable name. #' ## ----Create The Problem------------------------------------------------------- # create a new problem prob <- createprob() # change this problem to a maximization problem chgobjsense(prob,objsense = xpress:::OBJ_MAXIMIZE) # set the problem name setprobname(prob, "WaterSupply") #' #' #' Add the values we need for this example. #' ## ----Data--------------------------------------------------------------------- # nodes, where nodes 1 and 2 represent 2 reservoirs, nodes 8, 9 and 10 represent 3 cities, # nodes 3-7 represent pumping stations and nodes 11 and 12 are SOURCE and SINK NODES <- 1:12 SOURCE <- 11; SINK <- 12 # arcs & capacities ARCS.df <- data.frame(index = 1:20, head = c(1, 1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 5, 6, 7, 7, 8, 9, 10, 11, 11), tail = c(3, 5, 6, 5, 6, 4, 5, 8, 9, 8, 9, 10, 7, 9, 10, 12, 12, 12, 1, 2), capacity = c(20, 15, 12, 6, 22, 15, 10, 7, 10, 10, 15, 15, 22, 10, 10, 18, 15, 20, 35, 25)) #' #' #' Create variables 'flow' for each arc and set objective as mentioned in introduction. #' ## ----Add Columns-------------------------------------------------------------- # variables 'flow' for each arc ARCS.df$flow <- apply(ARCS.df, 1, function(x) xprs_newcol( prob, lb = 0, ub = x["capacity"], coltype = "C", name = pretty_name("flow", x["index"]), objcoef = NULL )) # change objective coefficients according to the definition of total flow chgobj(prob, colind = (ARCS.df %>% filter(tail == SINK))$flow, objcoef = rep(1, sum(ARCS.df$tail == SINK))) #' #' #' Add the flow conservation constraints. #' ## ----Add Rows----------------------------------------------------------------- # the total flow arriving at any node also has to leave this node (with the exception of the source and sink nodes) for (v in NODES[-c(SINK, SOURCE)]) { xprs_addrow( prob, colind = c((ARCS.df %>% filter(head == v))$flow, (ARCS.df %>% filter(tail == v))$flow), rowcoef = c(rep(1, sum(ARCS.df$head == v)), rep(-1, sum(ARCS.df$tail == v))), rowtype = "E", rhs = 0, name = paste0("balance_", v) ) } #' #' #' Now we can solve the problem. #' ## ----Solve The Problem-------------------------------------------------------- setoutput(prob) summary(xprs_optimize(prob)) #' #' #' Display the solutions here. #' ## ----The Solutions------------------------------------------------------------ paste("The optimum total flow is: ", getdblattrib(prob, xpress:::LPOBJVAL)) ARCS.df$solution <- xprs_getsolution(prob) invisible(apply(ARCS.df, 1, function(x) cat(x["head"], "->", x["tail"], ":", x["solution"], "/", x["capacity"], fill = TRUE))) # visualize the solution in network graph elist <- data.frame(tails = ARCS.df$head, heads = ARCS.df$tail) graph <- graph_from_data_frame(d = elist, directed = T) vlabel <- c( "1 reservoir", "2 reservoir", 3:7, "8 Gotham City", "9 Metropolis", "10 Spider Ville", "11 SOURCE", "12 SINK" ) elabel <- apply(ARCS.df, 1, function(x) paste0(x["solution"], "/", x["capacity"])) set.seed(66) plot( graph, vertex.color = "gold", vertex.size = 20, vertex.frame.color = "gray", vertex.label = vlabel, vertex.label.color = "black", vertex.label.cex = 0.8, layout = layout.auto(graph), edge.color = "gray", edge.label = elabel, edge.label.color = "red", edge.arrow.size = 0.5, edge.label.cex = 0.8 ) #' #'