Análisis numérico de flujo a costo mínimo en redes con estaciones de espera: el caso M/M/1 Numerical analysis of minimum cost network flow with queuing stations: the M/M/1 case

In a network the nodes represent stations, warehouses, distribution centers and customers and not just materials but also information circulate, so the minimum cost flow model, that only takes transport costs into account is one of the tools used to support the decision-making. In reality, the nodes provide a service which requires a service time, the servicing follows a discipline and also a queue is formed, generating the respective service and queuing costs. A modified version of the minimum cost flow model is proposed in this paper for the optimization of the flow in a queuing network. A variety of cases were solved. An acceptable level of accuracy was observed in the calculation of the time cycle and work in progress. The results indicate that the optimum solution Hernández González, Salvador et al. No 18, Vol. 9 (1), 2017. ISSN 2007 – 0705, pp.: 257 289 258 tends to balance the workload along the network. This paper is of interest to administrators and people in charge of supply chains and is useful for decision-making in the medium or long time.


Introduction
In industry and services, supply systems have a high degree of complexity, owing to the existence of the multiple relations that each element has with others within the same system. In consequence, it is often observed that work should be done with a network of stations/plants, warehouses/suppliers, containers/distribution centers and customers. Throughout such a network not only materials but also information circulates (Bhaskar and Lallement, 2010).
Some examples of the decisions that the administrators of these systems have to take are: how much product should be manufactured or how much should be sent from each workplace to the following link in the chain, with the total cost of dispatch being the common performance measurement ( Figure 1). Consider also a cross-docking system where goods and materials are transferred without actually entering the warehouse or being put away into storage (Kulwiec, 2004); a manager must determine the service time or the number of servers that reduces the residence time within the facilities, in other words the cycle time in the system. Optimal network design plays an important role in the supply chain operation, as good logistic distribution network can save transportation costs as well as improve service levels (Marmolejo, Rodríguez, Cruz-Mejía & Saucedo, 2016). The minimum cost network flow model is a useful tool for analyzing these types of systems.
Nº 18, Vol. 9 (1), 2017. ISSN 2007 -0705, pp.: 257 -289 -260 -It is worth mentioning that the nodes in a network, such as the one mentioned above are usually stations that provide a service, in other words, we are talking about queuing networks.
Analytical queuing models are used for determining properties such as work in progress, time in the system (cycle time), congestion of the system, measuring the quality of the service as well as quantifying the operating cost of said systems.
Networks may be open (the customers can enter and leave through any station on the network) or closed (the customer circulates around the network). Queuing network models are employed to represent not only production or manufacturing lines or, but also service and supply chain systems (Bolch, Greiner, de Meer and Trivedi, 2006).

Problem
The minimum cost network flow model only takes the transport costs into account; however every node on the network is often a station that provides a service (a factory, a warehouse). This operation generates a service cost (the products have to be manufactured) and a cost for the time taken to be served (the orders follow a service discipline and a queue is formed). In this context, managers and practitioners also need information about cycle time, work-in-process as well as the availability of resources at plants and warehouses. For a given topology of plants, warehouses and demand centers, it is evident that the solution obtained by a minimum cost flow model will change if each node is considered as a queueing system but, how much is the difference?; how sensitive is the solution if practitioners include transportation cost, waiting cost and service cost; what a practitioner and a manager can expect if they wish to measure supply chain performance as queueing networks.

Background
Most papers on queuing networks focus on the determination of the properties of the networks and use both analytical means and simulation. These models can be used for finding properties, control, planning and optimization, which is the case with this paper.
In the paper of Pourbabai, Blanc and van der Duyn Schouten (1996) Kerbache and MacGregor Smith (2000) propose a mathematical programming model that determines the optimum path that maximizes the throughput, the properties are calculated by applying the Generalized Expansion Method. Given that this is a mixed integer nonlinear model, an algorithm is developed to find good solutions, Yenisey (2006) applies concepts from Pourbabai et al (1996) to construct a mathematical programming model for MRP systems.
MacGregor Smith (2011) considers M/G/c, M/G/c/c systems and proposes a procedure for finding the optimum flow paths in a queuing network: In this case, the objective function is the system's throughput and is applied to the modeling of passenger flow in underground transport systems. Morabito, de Souza and Vázquez (2014) propose an analytical method for the shortest path problem in G/G/c queuing networks as well as the application of equations to improve the accuracy of the calculations of the system's properties, while, when applicable, the objective function also maximizes the throughput of the system. Van Woensel T., Cruz F.R.B. (2014) extended the model to G/G/c/K networks to determine the shortest route in systems with multiple

servers.
Some examples of applied cases we can mention are Bhaskar and Lallement's paper (2010) for the analysis of the supply chain of an Indian textile company, Srivathsan and Kamath (2012); they propose a model for the analysis of supply chains, but raise the idea of using the model to plan and find the optimum values of decision variables. He and Hu (2014) propose a model for decision-making in disaster prevention systems for the Shanghai area, while Beard and Chamberlain (2013) employ the principles proposed in Pourbabai. et al (1996) for the analysis of the flow of information in computer network systems.
Other analytical techniques for the modelling of supply chains are Petri Nets, Series parallel graphs, Markov chains and system dynamics (Srinivasa &Viswanadham, 2001).
This paper's contributions are: 1. This paper is a variant on the proposal of Pourbabai, et al (1996), but in this case the queuing and node operating costs are included in the objective function.
2. It differs from Morabito, Souza and Vázquez's paper (2014) in that different versions of the cost model are employed and assessed; the analysis is, for the moment, confined to the fact that the service time and demand follow an exponential distribution. 4. We are contributing with an approach for the process of modeling and analysis of a supply chain as a network of queues.
5. We obtain an optimal solution and calculate the work in process and cycle time of proeuction orders.
6. Numerical examples are solved and the results are compared with those obtained from a conventional minimum cost flow model. Simulation is used to validate analytical results of cycle time and work-in-process.
This paper is of interest to administrators and/or people in charge of supply chains, purchasing schedulers or production schedulers and is useful for decision-making in the medium or long time. With this model we get information about the cycle time, the work in process and the congestion at each station as well as the quantity of a product that must be sent to each customer in the network, minimizing the total operating cost of the system, this approach takes ino account the uncertainty of demand and service time of plants and distributors.

The minimum cost network flow problem
Consider a directed network that consists of a finite set of nodes N and a set of directed arcs A that join pairs of N nodes. Arc (i, j) is said to be incident to nodes i and j and is directed from node i to node j. Every node i is assigned a number fi that represents the available supply of an article (f > 0) or the demand for some article (f < 0). Nodes with supply are known as source nodes and the demand nodes are known as destination nodes. If fi= 0 then it is a transfer node, given that there is no supply or demand at said node.
Each arc is associated with a value xij that is the flow of material from node i to node j and a constant cij that is the cost of circulating through said arc. We shall that the total supply and the total demand are equal, in other words . The minimum cost network flow problem can be posed as follows: Shipping the available supply through the network in order to satisfy the demand at a minimum cost. Mathematically, this problem is expressed as follows: Minimize: (3) The constraints (2) are called network flow conservation equations. In conservation equations, the term represents the total flow that leaves node i towards any node j, while the term is the flow that enters node i from any other node k.
The equations require the net flow that leaves node i to be equal to fi. Finally, it must be pointed out that there can be flow capacity through the arcs.
The minimum cost flow model can be used to represent a logistics network where people and/or materials move through various points, the movement of locomotives, communications systems, oil pipelines, tanker schedules, etc (Bazaraa, Jarvis andSherali, 1974, Wagner, 1975).
The minimum cost network flow problem can be posed as a linear programming problem and the simplex method can be applied to find the optimum solution for a particular instance.
One of the assumptions in the aforementioned flow model is that no operation is done in the nodes and that the material flows without undergoing any kind of delay. However the people in charge of the administration and planning of the operations must take into account the time customers and materials spend in a queue at the station while waiting for the corresponding service; the time required to process a product, order or job; and the cost of the service, while they must also determine and analyze the system's properties such as the amount of work in process (WIP) and the cycle time in the system (from the time when the order is received until it is delivered to the end customer).
The aforementioned properties are performance measurements for queuing networks. A summary of the method for calculating said properties is given below.

Queuing networks
A queuing network is one where there are at least two stations connected to each other ( Figure 2).
Each node on the network represents a resource (for example a station or a factory). In this paper we use the following premises: 1. The stations have infinite capacity in the line.
2. Customers arrive at the network at the rate of λ customers per unit of time and travel from one station to another within the network.
3. Each station has si servers with a capacity of µi customers per unit of time.

Properties of queuing networks
Assume a network of n stations, the transit of the customers from one station to another is defined by the transit matrix P where each element pij is the probability that the customer goes from station i to station j. P is a square matrix. Assume that, through node i, γ0i customers per unit of time can enter the network from the exterior. Therefore the total flow to any node i on the network satisfies the following expression (Buzacott and Shanthikumar, 1993, Curry and Feldman, 2011, Bolch et al, 2006: Expressing this result in a matrix form: Where γ and λ are the flow column vectors from the exterior to each node and the entry flow to the node, respectively and P T is the transpose of the probability matrix. The system can be solved to obtain the entry flows to each station λ: Where I is the identity matrix. As can be observed, given γ and P, the equation calculates the entry flows to each station on the network.

Quadratic coefficients of the entry flows
After the equation (4), the entry flow to each station is made up by two elements: the entry from outside to the system through node i and the flow from any station within the network.
Assume that the external flow to a node j has a coefficient of variation that is denoted as , likewise, the coefficient of variation for the flow from a node k to node j is denoted as .
If the station has s servers, and each customer requires units of time to be served, with a coefficient of variation , also if each station that makes up the network is stable, in other words it holds that (the used capacity or congestion of the work station) therefore the quadratic coefficients of variation in the entry flows to each station are obtained with the following approximation: Once again this approximation can be expressed in matrix form and is solved as a system of linear equations, whose solution shall be: Where each element of the matrix Q is calculated as follows: And each element of the vector β is obtained with: Once the flows and their respective coefficient of variation have been obtained, the next step will consist of assessing the properties of each one of the network's station applying the decomposition principle. The expressions for the cycle time of station i and variability of departure from station i are: A summary of the procedure is given in figure 3 (Buzacott andShanthikumar, 1993, Curry andFeldman, 2011).
1. Read the vector of service times tS, the vector of coefficient of variation Cts 2 , the vector of entry flows to each station γ and the vector of coefficients of variation Ca 2 at the entry.

Read probability matrix P
3. Calculate the vector of entry flows λ at each station 4. Assess ρ for each station, if ρ > 1 for any station, stop, the

Calculate the vector of quadratic coefficients of variation
Ca 2 at the entry.
6. Calculate the properties of each station: Work in process, cycle time.

Cost function in a queuing network
The operating cost of a queuing network is expressed by the following equation: Conventionally, only the service cost and queuing cost are included in the operating cost (OC) in a queuing network. Likewise, only the network flow cost is included in a minimum cost flow problem.
Equation (13) includes the three terms: the cost arising from the operation of the servers ( ), the cost of having inventory in process in the station ( ) and finally the cost of transiting from station i to station j ( ).
As can be appreciated, flow xij between pairs of nodes is calculated using (4). Assuming that the external flows only appear at the supply nodes, then the term γ0i is zero, and the flow for any node other than the supply nodes is: -268 -

Optimization model
Assume that we have a network with a set of nodes and set of arcs . There is a set of nodes that represent the customers and generate a demand of units per unit of time; said demand follows a probability distribution G.
The material receives a transformation at each node that, on average, requires units of time and follows a probability distribution G. Therefore this is a network with G/G/s type stations. If and , then the station is Markovian or M/M/s, in accordance with the notation of Kendall.
There is no external flow to the demand nodes. Moreover, reflows are only considered for representing reprocesses, whereas exits from the system through any node other than the demand nodes are not considered.
The flow between the network's stations can be represented as a probability matrix P.
There is a service cost, Service Costi, for serving or producing a product unit, and there is also a cost per unit of work in process, Work in process Costi. Finally there is also a cost for transiting from one station to another cij. The proposed optimization model is as follows: Minimize: Subject to: , The quantity of material that should circulate between pairs of nodes or stations should be determined in such a way as to minimize the value of operation cost (16). It is worth mentioning that in the model, the decision variable is the probability or fraction of material that travels through the arc, with which the flow through the arc is calculated afterwards.
The flow is constrained to the entrythroughput balance at each node (17) being satisfied, as well as to the flow in the system being such that the congestion or used capacity at each station ρ (used capacity or congestion) are fulfilled or, if not, to ρ < 1 for all the stations (18).
Likewise, the way in which the flow is directed through the network is constrained to the range given in (19). In the model, the constraint (20) ensures that the material does not leave the network through any node that is not a demand node. In any event this constraint can be modified if necessary. The flows are calculated using (21).
The problem of routing customers, packets and other entity flows in finite queueing networks is an integer nonlinear programming. Our model is a variant of the Optimal Routing Problem in queueing networks which is known to be NP-hard Daskalaki, 1988, Kerbache andMacGregor Smith, 2000). The proposed model can incorporate constraints for the work in process or the cycle time. It is also feasible that the service time of one or more stations and the number of servers are also decision variables.
It is worth mentioning that, like any model, it is necessary to be careful to take the necessary factors into account and thus avoid situations such as getting a very complex model that does not have a feasible solution.
This model can be applied to the analysis and solution of problems involved with the planning of supply chain operations, using a networks approach and properties, such as the cycle time, work in process, the load at the stations or production or processing centers.

Implementation on a spreadsheet
It was necessary, for the study, to construct a spreadsheet template using the steps of the The template automatically calculates the total entry flows of each node using the equations (4) -(6); the flows through each arc with their respective cost; the costs of the stations; cycle time and amount variability of departures of each station; as well as the total operating cost of the system using the equations (7) - (12).
Work in process of each station is obtained with Little's Law: The template allows us to model systems that have up to 10 stations with a network arrangement or else with a pure series arrangement. Further details will be encountered in Appendix 1.

Numerical examples
A numerical analysis was carried out of several cases with a proposed minimum cost flow model, estimating properties (Work in process and cycle time), while the results were validated using simulation.
All cases were built taking into account that the maximum number of stages that are considered in mathematical models of supply chains are four (Olhager, Pashaei and Sternberg, 2015). Case A is a conventional minimum cost flow model without operation at the nodes. In series B, the minimum cost flow model with operation at the nodes was studied with and without costs (Table 1).
In case B1, decision variables are the probabilities pij of the customer going from node i to node j.
In case B2, the decision variables are pij as well as the service time at each station or node. In case B3, which is the model with service cost and queuing cost, the decision variable is the fraction pij that travels between pairs of nodes. .

B2 Yes Yes No
Fraction pij that is sent from node i to node j and service time at each node. The flow is calculated using (13).

B3 Yes Yes Yes
Fraction pij that is sent from node i to node j. The flow is calculated using (13).

B4 Yes Yes Yes
Fraction pij that is sent from node i to node j and service time at each node. The flow is calculated using (13

Series A
Consider a system such as that of figure 5. Nodes 1 and 2 are the sources of the supply, nodes 3 and 4 are manufacturing nodes, 5, 6 and 7 are transfer nodes and nodes 8 and 9 are the demand nodes. The costs per unit transported between pairs of nodes are given in table 3. The demand at nodes 8 and 9 is 30 units, while the supply available at nodes 1 and 2 is 30 units.  On solving the minimum cost flow model, using Solver and the data, we get the solution given in figure 6: 30 units should be sent from node 1 to node 3, and 30 units from node 2 to node 4. 30 units shall be sent from node 3 to node 5, 30 units shall be sent from node 4 to node 7.

Figure 6. Solution to case A.
Finally node 5 shall send 30 units to demand node 9 and node 7 shall supply 30 units to demand node 8 ( Table 4). The operating cost is $3570.

Series B. Minimum cost flow model with operation at the nodes
The analysis of the system is presented below, incorporating the data about the operation at each node, as well as the constraints on flow and use of the capacity.

Case B1
Consider the above minimum cost flow problem, taking into account that the nodes are stations where the customers receive a service, the data corresponding to the service time at each station are given in table 5. We will assume, in this case, that the service times are exponentially distributed, there are only transport costs between stations and that there is only one server at each node.
The decision variable is fraction pij which must be sent from station i to station j, this value is used to calculate amount of the flow of material between each station (Table 6).  The constraints for the value of pij are given in table 6. For this example in particular, the used capacity is constrained to the (25% -90%) range.
In this case we assume that the total demand of nodes 8 and 9 is shared out equally between nodes 1 and 2 (this could be the result of a decision taken beforehand by the person in charge in accordance with some prior criterion) in consequence, the entry flow to said stations is 30 each. The same table 6 shows the flow constraints for the supply nodes 1 and 2, as well as the demand nodes 8 and 9. Of transfer nodes 3 -7, only the equations corresponding to nodes 3 and 4 are shown.
The flows between stations were obtained by running the Solver complement and are given in table 7. The first significant point is the fact that the flow between stations is not an integer value, unlike a conventional minimum cost flow model. In the new solution, unlike Case A, the flow is The flow between stations requires the sending of material from station 1 to 4 and from station 2 to 3. The same happens with the flow from 3 to 5, 3 to 6, 3 to 7, 4 to 6, 4 to 7, 5 to 8, 6 to 8 and 7 to 9.
The operating cost obtained, being $3,776.25, is higher, which represents a difference of $206.25 or 2.8% in respect of case A. It is worth mentioning that this is normal as, owing to the conditions imposed on each node and the flow between stations, arcs were employed that had not been considered, thus raising said cost.
Additional information is obtained with model B1 for the administrators and decisionmakers: We observe, for example, that station 1 is only used at 30% of its used capacity ( or congestion); stations 2 and 3 are used at 84 and 81% of their capacity, respectively; while stations 6 and 7 are used at 67.5% of their capacity (Table 8). One relevant factor is the identification of the bottleneck. One accepted definition establishes that the bottleneck is the station where the most work in process is to be observed (Lawrence and Buss, 1995). In this example, stations 4 and 5 can both be considered to set the limits for the production capacity, as they are used at 90% of their capacity.
The cycle time for the process is the total time that elapses from the moment when the production order is received until it arrives at its destination with the customer: In the example, the total cycle time is 0.535 days and if we assume a basis of 1 month= 30 days, then the cycle time will be 16.05 days.
The amount of work (customers or production orders queuing in front of the station waiting for service) in the overall process is 32.096.
As regards the properties obtained by simulation, the work in process obtained by the analytical method differs by -2.46% in respect of results of the simulation, while the cycle time differs by-1.43%.

Case B2
The flows and cost that were found by analytical methods are given in table 9.In this case we observe that on being compared with case B1, the solution obtained presents several changes in both the flows and the properties of the system. In the flows we observe, for example, that the arc 2 -3 is not used. The total operating cost, which in this case is only the flow between pairs of nodes, is less in respect of case B1, presenting a reduction in cost of 0.84%.
As regards the performance measurements, case B2 presents a bottleneck at stations 4 and 5, where both are employed at 90% of their capacity, also there is a drop in the work in process and time cycle (Table 10).
A saving of 5.7% is obtained in the cycle time and there is a reduction of 5.67% in the work in process with respect to case B1.The analytical results for the work in process and the cycle time differ by -4.24% and -4.30%, respectively.  Figure 8 shows the comparison in cost of case A with cases B1 and B2. The cost of model B2 gives a better solution than case B1, which means it is worth considering the service time as a variable in order to improve the performance of the system.
As we have already mentioned, case A gives a lower cost, however this does not imply that it is a better solution, on the contrary, this is because additional elements are being incorporated into the model, to with: constraint on the flow, constraints on each node; all the above results in the feasible region shrinking.

Case B3
This model incorporates the service costs and queuing cost at each station. The values, which are given in table 11, are arbitrary as the only criterion for assigning them was for the queuing cost to be more expensive than the service cost. The rest of the model does not have any alterations. The comparisons shall only be made with model B4, which is described further on.  With the flows given in table 12, the properties of the system are shown below: used capacity, cycle time and work in process (Table 13). Here it is worth mentioning that a desirable goal for any administrator is to shorten the customer's queuing time and decrease the inventory in process. In both cases the customer would perceive an improvement in the service they receive: less of a wait to receive their product or else less of a wait to be served.
We can also see that there is only one bottleneck in this example and this is at station 5, which is at 90% of its capacity, having on average 9 units in a queue, which is followed, in terms of holdups, by stations 2, 3 and 4.
The analytical results for the work in process now show a difference of 6.21% and in the case of the cycle time, of 5.3% in respect of the results obtained by simulation.
As for the costs, having analyzed table 14 in detail, we observe that station 3 represents 26.7% of the total costs and what stands out is that this has the highest service cost, in other words, service is expensive at this station. Furthermore, when we look at the queuing cost column, we can see that station 5 has the highest congestion cost, therefore, from a monetary perspective, queuing at this station is the most expensive. Table 15 shows the costs of transport, the operating cost of the stations and the total cost of the system. In this example in particular, the cost of the transport represents approximately 60% of the total cost. should not be ruled out.
In a situation like the one in the example being analyzed, the decision-maker has several possible directions where they could propose improvements: they can take individual (service, congestion) costs into account; consider the overall cost; focus on suggesting improvements either to temper the effect of the production bottleneck or for the cost of transport which, in this case, represents 60% of the system's total operating cost.

Case B4
In this variant we must now determine the service times for stations that minimize the cost of flow in the network and take the values pij between pairs of nodes as fixed, so there is no change in the flows. The values are taken from the model B3 solution and Table 16 shows the flows that are obtained by analytical methods.  Station 3 contributes with 31.9% of the total costs of the system, followed by station 4 (Table   18). It is worth pointing out that both the transport cost and the operating cost of the stations are slightly lower than the ones obtained in case B4 (Table 19).  Figure 9. We can again see that it is a good idea to treat the service time as a variable, as this makes it possible to get a better solution: The savings in costs we obtain through finding the optimum service time and the optimum flow between stations is 6.58%, while the time cycle decreases by 10.83%. The comparison of B3 and B4 with cases A, B1 and B2 is not considered valid because said models represent other systems as they do not include either service costs or queuing costs.

Conclusions
The incorporation of new elements into a known model enriches the information that can be extracted on obtaining its solution.
In the case of the networks flow, the conventional minimum cost flow model can be expanded by assuming that each node is a station that provides a service j has a mean of Ts. In this case the addition of elements from the queuing networks makes it possible to obtain additional information about how the resources are being employed as well as, for example, getting metrics such as the time cycle and quantity of work in process. Equations of queeing systems take into account the stochastic component of a supply chain.
In this paper we tested a model for a single product and a single server, and can be applied to the analysis of a supply chain, however it can be broadened to work with a variety of products, Experiments showed that the optimal solution obtained with the minimum cost flow model changes if the nodes are considered as queueing systems. Another fact observed is that, when the nodes are queuing systems, the optimal solution obtained with the model tends to balance the workload in the nodes and also tends to distribute the flows through the network.
The equations for the properties (WIP and cycle time) are analytic approximations. We observe in the tests that were done that there is a difference between the analytical and simulated results, with the biggest difference being observed in case B3 and the smallest difference being observed in case B4. Even so, we can see that the analytical results have an acceptable level of accuracy.
Under the assumption of Poisson demand, service time exponentianlly distributed and stability of the system, managers and practitionares can expect that the proposed analytical model can be considered as a reliable tool for the analysis of a supply chain.
The conventional minimum cost flow model can be expanded to incorporate elements from queuing networks such as variability of demand, variability of service time, and variability due to failures or set-ups. -288 -destination costs, flow between nodes, cost of the flow between nodes, Operating cost of the stations, transport cost and total cost of the system (figure A2).
The shaded cells are the data that should be entered, for example: External flows (C16:C25) and their standard deviation of the arrivals (D16:D25), service times (O3:O12) with their standard deviation (P3:P12). Other data are service costs (O16:O25), queuing costs (P16:P25) and the costs of transfer between nodes(X3:AG12). Figure A2. Partial view of the spreadsheet template: Costs of transport and flows.

Flow in the network
First it is necessary evaluate the arrival flow to each station with eq. (6). Once the user has entered all data, the spreadsheet evaluates the transpose with the following function: After, the spreadsheet calculates the difference P t -I and then gets its inverse:

Variability of the arrivals
To solve eq. (8) the procedure is as follows: Each element of Q is calculated with eq. (9). Function IF is used to avoid division by zero.