Optimising capacity allocation in networks of stochastic loss systems: A functional-form approach

TL;DR

This paper presents a hybrid analytical-simulation method for efficiently optimizing capacity allocation in complex networks of stochastic loss systems, balancing congestion mitigation and cost.

Contribution

It introduces a novel iterative algorithm that combines analytical approximations with minimal simulation, significantly reducing computational effort while achieving near-optimal solutions.

Findings

Achieves near-optimal solutions comparable to advanced optimization methods.

Requires substantially less computational time than traditional simulation-based approaches.

Effective for complex networks with high computational efficiency.

Abstract

Motivated by a wide variety of applications, this paper introduces a general class of networks of stochastic loss systems in which congestion renders lost revenue due to customers or jobs being permanently removed from the system. We seek to balance the trade-off between mitigating congestion by increasing service capacity and maintaining low costs for the service capacity provided. Given the lack of analytical results and the computational burden of simulation-based methods, we propose a hybrid functional-form approach for finding the optimal resource allocation in general networks of stochastic loss systems that combines the speed of an analytical approach with the accuracy of simulation-based optimisation. The key insight is a core iterative algorithm that replaces the computationally expensive gradient estimation in simulation optimisation with a closed-form analytical approximation that is calibrated using a simple simulation run. Extensive computational experiments on complex networks show that our approach renders near-optimal solutions with objective function values that are comparable to those obtained using stochastic approximation, surrogate optimisation and Bayesian optimisation methods while requiring significantly less computational effort.

Tables

Table 1. Table 1 : Criss-cross network Model I experiment 1 capacity costs.

$l$	1	2	3	4	5	6	7	8	9	10
${\theta }_l$	0.20195	0.15437	0.16738	0.14339	0.1553	0.16866	0.27243	0.13134	0.12818	0.25142

Table 2. Table 2 : Criss-cross network Model II experiment 1 capacity costs.

$l$	1	2	3	4	5	6	7	8	9	10
${\theta }_l$	0.29064	0.12747	0.21388	0.29513	0.20067	0.23353	0.10684	0.19122	0.13117	0.19521
${\omega }_{1,l}$	9.8628	5.0975	3.3692	2.7648	8.6216	1.4089	7.3632	3.4109	2.171	1.2418
${\omega }_{2,l}$	8.7339	2.7184	2.1695	3.4332	1.567	7.5708	3.5566	6.668	5.6143	3.1974

Equations

$$f(\boldsymbol{c})=-\langle\boldsymbol{c},\boldsymbol{{\theta}}\rangle+\sum_{r\in\mathcal{R}}\lambda_{r}\,{\omega}_{r}\prod_{i=1}^{N_{r}}\big{(}1-p_{r,i}(\boldsymbol{c})\big{)}\,.$$

$$f(\boldsymbol{c})=-\langle\boldsymbol{c},\boldsymbol{{\theta}}\rangle+\sum_{r\in\mathcal{R}}\lambda_{r}\sum_{i=1}^{N_{r}}{\omega}_{r,i}\big{(}1-p_{r,i}(\boldsymbol{c})\big{)}\prod_{j=1}^{i-1}p_{r,j}(\boldsymbol{c})\,.$$

$${\mathbf{c}}^{\star }:=\underset{\mathbf{c}\in {\mathbb{Z}}_{+}^L}{\mathrm{argmax}}f(\mathbf{c}),$$

$$f(\boldsymbol{c})=-{\theta}_{1}\,c_{1}-{\theta}_{2}\,c_{2}+\lambda_{1}\,{\omega}_{1}\,\big{(}1-p_{1,1}(c_{1})\big{)}\big{(}1-p_{1,2}(\boldsymbol{c})\big{)}\,.$$

$$f(\boldsymbol{c})=-{\theta}_{1}\,c_{1}-{\theta}_{2}\,c_{2}+\lambda_{1}\,{\omega}_{1,1}\,\big{(}1-p_{1,1}(c_{1})\big{)}+\lambda_{1}\,{\omega}_{1,2}\,\big{(}1-p_{1,2}(\boldsymbol{c})\big{)}p_{1,1}(c_{1})\,.$$

$${\mathbf{c}}^{\star }:=\underset{\mathbf{c}\in \mathcal{C}}{\mathrm{argmax}}f(\mathbf{c})$$

$${\tilde{\mathbf{c}}}^{\star }=\underset{\mathbf{c}\in \mathcal{C}}{\mathrm{argmax}}\tilde{f}(\mathbf{c},\mathbf{\tau })$$

$${\tilde{f}}^{(n)}(\cdot ):=\tilde{f}(\cdot ,{\mathbf{\tau }}^{(n)}),n=1,2,\dots ,$$

$${\mathbf{c}}^{(n)}:=\underset{\mathbf{c}\in \mathcal{C}}{\mathrm{argmax}}{\tilde{f}}^{(n)}(\mathbf{c}),n=1,2,\dots ,$$

$$\hat{f}({\mathbf{c}}^{(n-1)}):=\tilde{f}({\mathbf{c}}^{(n-1)},{\mathbf{\tau }}^{(n)}),n=1,2,\dots .$$

$$\tilde{p}(c,{\tau }_{r,i})=\exp \left(-(c{\tau }_{r,i}{)}^2\right).$$

$${\hat{p}}_{r,i}({\mathbf{c}}^{(0)})=\tilde{p}\left(c_{l(r,i)}^{(0)},{\tau }_{r,i}^{(1)}\right).$$

$$\tau_{r,i}^{(1)}=\frac{1}{c_{l(r,i)}^{(0)}}\Big{(}-\log\big{(}\widehat{p}_{r,i}\left(\boldsymbol{c}^{(0)}\right)\big{)}\Big{)}^{1/2}.$$

$$\underset{\boldsymbol{c}^{(1)}\in\mathbb{[}0,\infty)^{L}}{\operatorname*{arg\,max}}\quad-\langle\boldsymbol{c}^{(1)},\boldsymbol{{\theta}}\rangle+\sum_{r\in\mathcal{R}}\lambda_{r}\,{\omega}_{r}\,\prod_{j=1}^{i}\left(1-\widetilde{p}\big{(}c_{l(r,i)}^{(1)},\,\tau_{r,i}^{(1)}\big{)}\right)\,,$$

$$\underset{\boldsymbol{c}^{(1)}\in\mathbb{[}0,\infty)^{L}}{\operatorname*{arg\,max}}\quad-\langle\boldsymbol{c}^{(1)},\boldsymbol{{\theta}}\rangle+\sum_{r\in\mathcal{R}}\lambda_{r}\sum_{i=1}^{N_{r}}{\omega}_{r,i}\Big{(}1-\widetilde{p}\big{(}c_{l(r,i)}^{(1)},\,\tau_{r,i}^{(1)}\big{)}\Big{)}\prod_{j=1}^{i-1}\widetilde{p}\big{(}c_{l(r,j)}^{(1)},\,\tau_{r,j}^{(1)}\big{)}\,,$$

$$\boldsymbol{c}^{(n)}=\underset{\boldsymbol{c}\in\mathbb{[}0,\infty)^{L}}{\operatorname*{arg\,max}}\quad-\langle\boldsymbol{c},\boldsymbol{{\theta}}\rangle+\sum_{r\in\mathcal{R}}\lambda_{r}\,{\omega}_{r}\,\prod_{i=1}^{N_{r}}\Big{(}1-\exp\left(-\big{(}c_{l(r,i)}\tau_{r,i}^{(n)}\big{)}^{2}\right)\Big{)}\,,$$

$$\boldsymbol{c}^{(n)}=\underset{\boldsymbol{c}\in\mathbb{[}0,\infty)^{L}}{\operatorname*{arg\,max}}\quad-\langle\boldsymbol{c},\boldsymbol{{\theta}}\rangle+\sum_{r\in\mathcal{R}}\lambda_{r}\sum_{i=1}^{N_{r}}{\omega}_{r,i}\Big{(}1-\exp\left(-\big{(}c_{l(r,j)}\tau_{r,i}^{(n)}\big{)}^{2}\right)\Big{)}\prod_{j=1}^{i-1}\exp\left(-\big{(}c_{l(r,j)}\tau_{r,j}^{(n)}\big{)}^{2}\right)\,.$$

$$\frac{{\lambda }^c/c!}{{\sum }_{i=0}^c{\lambda }^i/i!}.$$

$$\frac{1.96\sqrt{1.96^2-4A_{r,i}^{(n)}({\hat{p}}_{r,i}(\mathbf{c})-1){\hat{p}}_{r,i}(\mathbf{c})}}{1.96^2+A_{r,i}^{(n)}},$$

$${\kappa }_1:={\lambda }_1{\omega }_1(1-p_{1,1}(c_1))(1-p_{1,2}(\mathbf{c})),{\kappa }_{1,1}:={\lambda }_1{\omega }_{1,2}(1-p_{1,1}(\mathbf{c})),{\kappa }_{1,2}:={\lambda }_1{\omega }_{1,2}(1-p_{1,2}(\mathbf{c}))p_{1,1}(c_1),$$

$$\mathbb{E}{N}_1(t)$$

$$\mathcal{D}={\mathcal{D}}_0+{\mathcal{D}}_1$$

$${\mathcal{D}}_0=\text{diag}([{\mathcal{Q}}_0,\dots ,{\mathcal{Q}}_{c_2}])+\text{diag}([0,\dots ,1])\otimes {\mathcal{M}}_1+I_{c_2+1}\otimes \mathcal{L}+{\mathcal{M}}_2\otimes I_{c_1+1},{\mathcal{D}}_1={\overline{I}}_{c_2+1}\otimes {\mathcal{M}}_1,$$

$$\mathbb{E}{N}_1(t)=\tilde{\mathbf{\pi }}{\tilde{\mathcal{D}}}_1{1}_{(c_1+1)(c_2+1)}t+\text{o}(t),\mathbb{E}{N}_2(t)=\tilde{\mathbf{\pi }}{\tilde{\mathcal{D}}}_2{1}_{(c_1+1)(c_2+1)}t+\text{o}(t),$$

$$\tilde{\mathcal{D}}={\tilde{\mathcal{D}}}_0+{\tilde{\mathcal{D}}}_1+{\tilde{\mathcal{D}}}_2$$

$${\tilde{\mathcal{D}}}_0=\text{diag}([{\mathcal{Q}}_0,\dots ,{\mathcal{Q}}_{c_2}])+{\mathcal{M}}_2\otimes I_{c_1+1}+I_{c_2+1}\otimes {\mathcal{M}}_1-{\mathring{I}}_{c_2+1}\otimes \Lambda ,{\tilde{\mathcal{D}}}_1=I_{c_2+1}\otimes \mathcal{L},{\tilde{\mathcal{D}}}_2={\overline{I}}_{c_2+1}\otimes \Lambda .$$

$$\sum _{n=1}^{\infty }{\alpha }^{(n)}=\infty ,\sum _{n=1}^{\infty }{\alpha }^{(n)}{\delta }^{(n)}<\infty ,\sum _{n=1}^{\infty }{{\alpha }^{(n)}}^2{{\delta }^{(n)}}^{-2}<\infty .$$

Full text

1 Introduction

Stochastic networks encountered in practice often involve losses due to various forms of congestion. A key design decision for such networks concerns the allocation of (expensive) resources to limit losses while maintaining low costs. For example, in cellular wireless networks, the calls of moving users may be dropped if the receiving cell tower has insufficient capacity available for a successful handover (see, e.g., [1, 27, 40]).

Given the ubiquity of wireless networks, the allocation of resources in these networks plays a crucial role in many real-world applications experienced throughout daily life.

A similar situation may arise in mobile cloud computing, where computational jobs are served at small “cloudlets” close to the user, rather than at a remote cloud computing facility (see, e.g., [31, 20, 28]). While the use of cloudlets reduces delays [20], each cloudlet has limited capacity and migrating jobs between cloudlets due to user mobility may lead to job losses and performance issues.

Given the growth of such cloudlet-based mobile computing, resource allocation in these environments plays a critical role in various real-world applications including augmented reality, crowd-sourcing video processing, and facial recognition [31].

These processes based on wireless networks and cloudlet mobile computing can be formulated as part of a general class of networks of stochastic loss systems (SLSs).

Another example within the general class of networks of SLSs concerns supply chain networks designed for the manufacturing and maintaining of computer and network systems with uncertain lifetimes. Any downtime of these systems can be expensive and disruptive, and should be minimised to the extent possible. When such a computer or network system breaks down, the necessary spare parts are quickly dispatched from the nearest stock point. If a spare part is unavailable at this location, losses occur and the parts must be sent from a more distant central warehouse, thus increasing downtime (see, e.g., [25, 34, 39]) and in turn making resource allocation crucially important. An analogous problem arises in the context of emergency services, where ambulances and fire trucks are geographically dispersed to ensure a fast response to emergencies (see, e.g., [7, 42, 43]). Requests for these vehicles (i.e., a nearby emergency) may be rejected because they are already in service, resulting in the risk of loss of life and requiring a rerouting of the request to other vehicles/resources with associated delays. The allocation of ambulance and fire fighting resources, possibly through computer-automated support, obviously plays a critically important role in emergency services. These processes based on supply chain and emergency service networks can also be formulated as part of a general class of networks of SLSs.

In this paper we consider such a general class of networks of SLSs, with a particular focus on both of the above stochastic network model instances within the general class. Given the importance of resource allocation in this class of networks of SLSs, we observe that there exist essentially two categories of capacity management approaches in stochastic networks in general: analytical and simulation-based. Analytical approaches typically approximate the complex stochastic network with a much simpler network and solve the corresponding optimisation problem. These approximations may be based on product-form networks (see, e.g., [5, 15, 22]), independence assumptions among network components (see, e.g., [4, 36]), aggregation of network components (see, e.g., [41]), fluid-limit scalings (see, e.g., [16]), and large-scale systems (see, e.g.,  [6, 17]). These closed-form approximations can then be used to solve the capacity management problem, either analytically (see, e.g., [24, 44]), numerically, or using a combination of both (see, e.g., [38]).

Simulation-based approaches typically evaluate and update capacity allocation using simulation of the stochastic network. One canonical simulation-based optimisation approach is stochastic approximation (SA) (see, e.g., [3, 19, 23, 37] and Appendix B), where gradients are estimated using simulation in order to apply gradient descent in a stochastic setting. Similar adaptations of Newton and quasi-Newton methods, based on higher-order derivatives, are currently a very active area of research (see, e.g., [10]). Another well-known approach is Bayesian optimisation (BO) (see, e.g., [9, 29]), where a prior is placed over the unknown objective function which is then updated via noisy observations. Surrogate optimisation (SO) is yet another popular approach (see, e.g., [14, 21, 35]), which optimises over a response surface that is generated by radial basis function interpolation between function evaluations.

Both analytical and simulation approaches have limitations that prevent application to the complex stochastic loss networks of interest. In analytical approaches, it is often impossible to find an approximation for the complex network that is sufficiently accurate and captures all relevant features; as a result, any capacity management decisions based on the approximate model may not, in general, work well in the original model. Simulation-based approaches achieve better accuracy, but can suffer from large computational costs which often make their application to high-dimensional complex stochastic loss networks, as well as their application at finer time scales, infeasible. These large computational costs may be driven by the need to estimate a (high-dimensional) gradient, which requires many simulation runs to obtain estimates with sufficient accuracy.

A recent paper [13] introduced a general hybrid approach for a class of stochastic networks without losses that exploits its theoretical properties to significantly reduce the computational costs of simulation-based optimisation. The proposed approach consists of two phases where simulation is used in the first phase to fit a closed-form functional-form approximation to the performance metric at the current resource allocation level. This approximation is then used to analytically determine a new resource allocation, rather than using gradient descent, and the process is repeated iteratively. A general version of the corresponding first-phase iterative process is illustrated in Figure 1. It is shown in [13] that the resource allocation obtained in the first phase is close to optimal while only requiring a fraction of the computational effort of existing methods, due to reducing the number of iterations and due to not having to estimate gradients from simulation in each iteration. The approach is similar to response surface or surrogate optimisation methods (see, e.g., [2, Section 3.2]), but differs due to the significantly higher level of detail that is analytically garnered about system structure and incorporated into the functional-form approximation model. Once the first-phase iterative process converges, the resulting near-optimal solution is then used as a starting point for the SA method as part of the proposed two-phase approach.

The first-phase approximation approach of Dieker et al. [13] does not carry over to the general class of networks of SLSs studied in this paper, since the addition of losses fundamentally changes the dynamics and other characteristics of the stochastic network. In particular, the objective changes from a function of the queue lengths in [13] to a function of the probabilities of losses; we are also interested in discrete capacities for the number of servers at each station of the network rather than continuous per-station service rates. Moreover, because the problem in [13] considers each queue to have infinite buffer capacity and thus no jobs are lost, the arrival rate into every station is independent of the capacity allocation subject to stability constraints. Then under the first-phase functional-form approximation in [13], the problem can be essentially decoupled into smaller sub-problems for each station. Since losses remove jobs from the system and losses occur more frequently when there is reduced capacity, we can no longer exploit in our setting the property that the effective long-term arrival rate into a station is independent of the capacity allocation of the upstream station(s). This added level of interaction among the stations substantially increases the problem complexity.

In this paper our focus is on the first-phase approximation approach of Dieker et al. [13] and we propose a general iterative approach for optimising capacity allocation in networks of SLSs that extends this first-phase approach in [13] to address the important complexities of these networks. Our core iterative algorithm explicitly preserves the dependence among stations by replacing all appearances of a loss probability in the objective function with a unique functional-form approximation—the capacity vector is then jointly optimised using numerical methods. The efficiency of our core algorithm is driven by the known relationships between losses at each station for each route, which are explicitly incorporated into the optimisation method. In order to encapsulate the full range of applications listed earlier, we detail our core algorithm for two model variants within the overarching general class of networks of SLSs. The primary difference between the two models is the objective of the jobs in the network. For the first model (Model I), each job aims to visit a subset of the stations and is lost as soon it encounters a station which is at capacity; whereas for the second model (Model II), each job aims to visit one station out of some subset of the stations and is only lost if all of the stations in the subset are at capacity. The core iterative algorithm of our first-phase approximation approach is not a generic simulation optimisation algorithm; it has been tailor-made specifically for the general class of networks of SLSs and should be used as such. At the same time, by detailing two models within this general class we demonstrate the general principles of how to use our first-phase approach within the class.

Our investigation into the performance of the core iterative algorithm of our first-phase approximation approach consists of extensive computational experiments on four example topologies of general networks of SLSs: Markovian tandem networks, criss-cross networks, ring networks, and a complex class of canonical networks. We consider all four topologies in the context of both of the model variants introduced above. Many of our examples feature renewal arrival processes and/or non-exponential service times and are thus represented by non-Markovian stochastic processes. The first of our examples is (just) simple enough that a benchmark for the accuracy of our core iterative algorithm against true optimal solutions is possible for both model variants. In addition to comparing the accuracy of our core algorithm with the exact solution for this simplest network topology example, we gauge the efficiency of our approach by benchmarking it against the SA, BO and SO algorithms introduced earlier. We also use these benchmark algorithms to gauge the accuracy and efficiency of our core iterative algorithm for the three more complex example network toplogies, SLSs of interest, with our focus on extending the first-phase approach in [13] to support the general class of networks of SLSs of interest.

For the Markovian tandem network example, we observe that our core iterative algorithm directly achieves near-optimal solutions with very low computational effort. By comparison, the other benchmark algorithms are occasionally slightly more accurate, but even for this relatively simple model these benchmark algorithms by themselves already exhibit much higher computational costs. For the remaining example topologies, we observe that our core algorithm makes tremendous strides within a few computationally inexpensive iterations to greatly improve the performance of the system from the initial starting point of random capacity allocation. We show that our core iterative algorithm compares quite favourably against SA both in terms of objective value found and in terms of computational efficiency. We further provide evidence that our core algorithm is frequently faster than SO and usually much faster than BO. These massive benefits in computational efficiency are at worse realized with only a very minor reduction in the objective values found, and indeed, for many of the scenarios we consider, our core algorithm outputs better solutions than the benchmark methods. Hence, the solution of our core iterative algorithm can be either used at finer time scales to quickly obtain a high-quality solution or used at coarser time scales to quickly provide a high-quality initial starting point for a preferred benchmark algorithm to efficiently obtain a refined solution, following the two-phase approach in [13].

In summary, our contributions are as follows:

We introduce a new general class of networks of SLSs with a wide range of applications across many areas, which has been previously ignored because of its intractability; 2. 2.

By extending the first-phase approach in [13], we devise a novel core iterative algorithm of our general first-phase approximation approach that supports this new class of stochastic models; 3. 3.

We apply our first-phase approach to two model variants within the general class of networks of SLSs, and show that it performs remarkably well.

The remainder of this paper is organised as follows. Section 2 presents our general class of networks of SLSs, including two canonical instances of this class of stochastic network models together with their corresponding objective functions. Section 3 provides an overview of our iterative functional-form optimisation approach for first-phase solutions to our general class of networks of SLSs. Section 4 demonstrates that our core algorithm works well on realistic network examples under a variety of scenarios, outperforming several benchmark approaches. This is followed by concluding remarks in Section 5. In addition, Appendix A shows how to compute the objective function for a small Markovian network instance, which is used as a benchmark in one of our computational experiments.

2 Mathematical models

We now define a general class of networks of SLSs and the corresponding optimisation problem formulation, followed by a small example for illustration and elucidation. Consider a set of stations $\mathcal{L}=[L]:=\{1,\dots,L\}$; e.g., each station could represent a base station in a cellular wireless network or a stock point in a computer spare-parts network. Station $l$ has $c_{l}$ servers that can each serve one customer at any given time; we refer to $c_{l}$ as the capacity of station $l$ and define $\boldsymbol{c}:={(c_{l})}_{l\in[L]}$. Define $\mathcal{R}:=[R]$ to be the set of customer classes. Customers of class $r\in\mathcal{R}$ arrive to the network according to some general exogenous arrival process with rate $\lambda_{r}$. No further assumptions are made on the arrival processes, allowing as examples renewal type with independent and identically distributed interarrival times (see, e.g., [8, 32]) and Markovian arrival processes with non-identically distributed interarrival times (see, e.g., [18, 30]). Associated with each class $r$ is a path comprising $N_{r}$ network stations defined by $\psi_{r}:={(l(r,i))}_{i\in[N_{r}]}$, where $l(r,i)\in\mathcal{L}$ is the $i$-th station in the class-$r$ path. We assume no station appears more than once in any path $\psi_{r}$, i.e., $l(r,i)\neq l(r,j)$, $\forall~{}i\neq j$. A class-$r$ customer arrives to the network at station $l(r,1)\in\mathcal{L}$.

As pointed out in the introduction, two variants of this class of stochastic networks are to be distinguished. The key differences between these variants are the goal of the customers, the effect that this has on their behaviour in response to congestion, and the implications for the revenue stream of the network operator. Our general functional-form optimisation approach then exploits underlying theoretical properties of the objective function (i.e., existing domain knowledge) to augment information from simulation and speed up the iteration process of finding a good approximation to the optimal solution. The two variants are described as follows.

• •

In the first variant, referred to as Model I, the objective of customers is to be served, without interruption, by all of the stations along their path. For this variant, when a customer arrives at a station to find no available server, it will immediately depart the network and we say the customer has been lost. However, if a station to which a customer arrives has an available server, it is accepted for service and remains at the server of the station for a generally distributed amount of time with finite expectation—we impose no other restrictions on the service distributions. Upon completion of service at station $l(r,i)$, for $i\in[N_{r}-1]$, the customer transitions to station $l(r,i+1)$ and either receives service or is lost, depending on whether a server is available at the station. The customer leaves the network after receiving service at station $l(r,N_{r})$. Customers who successfully traverse all the stations of their path $\psi_{r}$ pay an amount ${\omega}_{r}$ to the network operator. Customers who are lost while traversing the network do not pay the network operator. Model I is suitable, for example, as a representation of a cellular wireless network where calls of moving users may be dropped if the receiving cell tower has insufficient available capacity for a successful handover. A similar application is cloudlet-based mobile computing, where the jobs of mobile users may be lost and/or encounter performance issues if the limited capacity of cloudlets is not sufficient to properly handle the migration of these jobs along the path of the user.

• •

In the second variant, referred to as Model II, the objective of customers is to be served at one and only one of the stations along their path. For this variant, when a customer arrives at a station to find no available server, it will immediately attempt service at the next station on its path. If a customer is denied service from all of the stations in its path, then it is lost. If a station to which a customer arrives has an available server, it is accepted for service and remains at the server of the station for a generally distributed amount of time with finite expectation—we impose no other restrictions on the service distributions. Upon acceptance for service at station $l(r,i)$ a customer pays ${\omega}_{r,i}$ to the network operator. Typically we would expect that ${\omega}_{r,i+1}<{\omega}_{r,i}$, reflecting a preference for stations that are earlier in the path (though not a requirement of our approach). Upon completion of service at station $l(r,i)$ the customer departs the network. Model II is suitable, for example, as a representation of a supply chain network where customers who cannot obtain a spare part for a computer system at their preferred depot will seek to obtain the part at other depots. Another application is emergency services, where if all of the units for a particular district are occupied when requested, then units from nearby districts can be redirected at a higher cost (human and financial).

Objective functions. Our study addresses the fundamental trade-off between maintaining low capacity costs and collecting high rewards (by ensuring few losses). Consistent with the dynamics described above, the objective function takes on a different form depending on which variant of the model is being considered. For both variants, each server at station $l$ costs ${\theta}_{l}>0$ per unit time, and we define $\boldsymbol{{\theta}}:={({\theta}_{l})}_{l\in[L]}$. This means that the network operator pays $\langle\boldsymbol{c},\boldsymbol{{\theta}}\rangle$ per unit time, where $\langle\cdot\rangle$ denotes the usual Euclidean inner product. Let $p_{r,i}(\boldsymbol{c})$ denote the fraction of class-$r$ jobs arriving at the $i$-th station on their path $\psi_{r}$ that are then lost at the $i$-th station. We express the performance of the above system in terms of the expected net rate of reward generation by the system operating in equilibrium. For both variants, increased losses result in reduced revenue for the network operator. The difference between the variants is in how rewards are collected by the network operator, which is described as follows.

• •

For Model I, once a customer of class $r$ has successfully entered all of the stations in its path $\psi_{r}$, the network operator receives revenue ${\omega}_{r}$. This means that in equilibrium with capacity $\boldsymbol{c}$ the network operator receives payments of value ${\omega}_{r}$ at rate $\lambda_{r}$ thinned by \prod_{i=1}^{N_{r}}\big{(}1-p_{r,i}(\boldsymbol{c})\big{)}. The latter product represents the proportion of class-$r$ customers that are not lost at any of the stations in $\psi_{r}$. Summing these revenue streams across all of the classes $r\in\mathcal{R}$ together with the costs of providing servers $\langle\boldsymbol{c},\boldsymbol{{\theta}}\rangle$ results in the objective function

[TABLE]

This expression serves as our Model I objective function, which is inspired by the capacity value function investigated in [11, 12] within the context of individual Erlang-B loss systems and later generalised in [33] within the context of provisioning cloud computing platforms.

• •

For Model II, when a customer of class $r$ successfully enters station $l(r,i)$ along its path $\psi_{r}$, the network operator receives revenue ${\omega}_{r,i}$. In order for this revenue to be obtained, the customer must be lost at stations $l(r,1),\ldots,l(r,i-1)$ before being accepted at station $l(r,i)$. Hence, in equilibrium with capacity $\boldsymbol{c}$, the network operator receives payments of value ${\omega}_{r,i}$ at rate $\lambda_{r}$ thinned by \big{(}1-p_{r,i}(\boldsymbol{c})\big{)}\prod_{j=1}^{i-1}p_{r,j}(\boldsymbol{c}). This product represents the proportion of class-$r$ customers that are lost at stations $l(r,1),\ldots,l(r,i-1)$ and then accepted at station $l(r,i)$. Combining these revenue streams across all of the classes $r\in\mathcal{R}$ together with the costs of providing servers $\langle\boldsymbol{c},\boldsymbol{{\theta}}\rangle$ results in the objective function

[TABLE]

This expression serves as our Model II objective function.

For both model variants, our optimisation formulation seeks to determine

[TABLE]

where ${\mathbb{Z}}_{+}$ denotes the non-negative integers.

Example. In order to elucidate the two variants of our general class of networks of SLSs, we now briefly discuss a simple illustrative example. Consider a system with two stations, a single customer class, and capacities $c_{1}$, $c_{2}$;. Figure 2(a) illustrates the Model I variant of this system and Figure 2(b) illustrates the Model II variant of this system. Customers arrive to station $1$ according to a Poisson process with rate $\lambda_{1}$; upon successful service at station $1$ under Model I, the customer attempts service at station 2, whereas upon a loss at station 1 under Model II, the customer attempts service at station $2$. Hence, $\mathcal{L}=\{1,2\}$, $\mathcal{R}=\{1\}$, $\psi_{1}=(1,2)$, $l(1,1)=1$, $l(1,2)=2$. Observe that, while $p_{1,1}(c_{1})$ is independent of $c_{2}$, $p_{1,2}(\boldsymbol{c})$ depends on both $c_{1}$ and $c_{2}$ because larger $c_{1}$ renders a larger (smaller) number of arrivals to station $2$ under Model I (Model II). Capacity costs are ${\theta}_{1}$ and ${\theta}_{2}$ per unit time. In the case of Model I, rewards ${\omega}_{1}$ are obtained by the network operator per successful entry at station $2$, hence the Model I objective function (1) is given by

[TABLE]

In the case of Model II, rewards ${\omega}_{1,1}$ and ${\omega}_{1,2}$ are obtained by the network operator per successful entry to stations $1$ and $2$ respectively, hence the Model II objective function (2) is given by

[TABLE]

In both variants the capacity $\boldsymbol{c}$ is to be maximized over a set of possible network parameters $\mathcal{C}={\mathbb{Z}}_{+}^{2}$.

3 Functional-form approach

We now present our general approximation approach that seeks to extend the first-phase approach in [13] to our setting of networks of SLSs. A detailed description of our general approach and the underlying ideas is provided first, and then we turn our attention to applying this approach in the specific setting of interest.

3.1 Overview of general approach

The objective of our study is to determine parameters $\boldsymbol{c}$ that maximise a function $f:\mathcal{C}\to\mathbb{R}$ representing a performance metric of interest. In this paper we consider $\boldsymbol{c}$ to be the number of servers in a network of SLSs and $f(\boldsymbol{c})$ to be a combination of the weighted ergodic fraction of lost customers and the server provisioning costs. We note that Dieker et al. [13] consider $f(\boldsymbol{c})$ to be the weighted expected number of customers in a stochastic network (without losses) in equilibrium.

A primary difficulty in determining the value

[TABLE]

comes from the fact that the function $f$ is typically unknown, owing to the complexity of stochastic networks (especially in non-Markovian settings), and therefore we assume it can only be evaluated as a random variable (r.v.) $F$ such that $\mathbb{E}F(\boldsymbol{c})=f(\boldsymbol{c})$. With time-consuming simulation, we assume it is possible to obtain independent samples of $F$, denoted by $\widehat{f}$, such that the information contained in these samples can be used to find an approximation of $\boldsymbol{c}^{\star}$. The objective of this paper is to develop a method for closely approximating $\boldsymbol{c}^{\star}$ using as few evaluations of $F$ as possible, for each of the model instances detailed in Section 2.

Our iterative functional-form optimisation approach exploits underlying theoretical properties of $f$ (i.e., existing domain knowledge) to augment information from simulation and speed up the iteration process of finding a good approximation to $\boldsymbol{c}^{\star}$. This structural information is expressed as a closed-form function $\widetilde{f}(\boldsymbol{c},\boldsymbol{\tau})$ of both the network parameters $\boldsymbol{c}$ and some additional (potentially vector-valued) parameters $\boldsymbol{\tau}$. The parameters $\boldsymbol{\tau}$ are then used to tune $\widetilde{f}$ so that it fits $f$ well locally. Depending on the complexity of $\widetilde{f}$ we then solve, either analytically or numerically, for

[TABLE]

to approximate $\boldsymbol{c}^{\star}$. The function $\widetilde{f}$ is selected to ensure that (7) can be solved in closed form or using a fast numerical procedure, which is in stark contrast to solving for (6).

The quality of the approximation (7) relies in large part on choosing an $\widetilde{f}$ that properly represents fundamental behaviours of the stochastic network. We assume $\boldsymbol{\tau}$ is selected from an appropriate set such that the set of functions \big{\{}\widetilde{f}(\cdot,\boldsymbol{\tau})\big{\}}_{\boldsymbol{\tau}} has elements that approximate $f$ well and that can be reliably identified in terms of $\boldsymbol{\tau}$. Let us initially assume a good functional form is known (which is addressed below in Section 3.2) and provide an iterative procedure for choosing $\boldsymbol{\tau}$, where the sequence of samples $\widehat{f}(\boldsymbol{c}^{(1)}),\ldots,\widehat{f}(\boldsymbol{c}^{(N)})$ is used to guide the selection of $\boldsymbol{\tau}$ over a maximum number of iterations $N$.

Starting with an initial value of $\boldsymbol{c}=\boldsymbol{c}^{(0)}$, which can be chosen at random or via a crude analytical approximation, we first evaluate $\widehat{f}(\boldsymbol{c}^{(0)})$ using simulation. Here $\widehat{f}(\boldsymbol{c}^{(0)})$ represents the objective value at $\boldsymbol{c}=\boldsymbol{c}^{(0)}$ obtained through simulation. The next step is to set $\widehat{f}(\boldsymbol{c}^{(0)})=\widetilde{f}(\boldsymbol{c}^{(0)},\boldsymbol{\tau}^{(1)})$ and solve for $\boldsymbol{\tau}^{(1)}$. This renders our iteration $1$ approximation function $\widetilde{f}^{(1)}(\cdot):=\widetilde{f}(\cdot,\boldsymbol{\tau}^{(1)})$, as depicted in Figure 3(a)111Although Figure 3 depicts the domain of $f$ to be one-dimensional for ease of illustration, we are generally interested in high-dimensional stochastic networks, which represent a primary cause of complexity for the problems under consideration.. Then, from (7), we solve $\boldsymbol{c}^{(1)}:=\operatorname*{arg\,max}_{\boldsymbol{c}}\,\widetilde{f}^{(1)}(\boldsymbol{c})$ to obtain our iteration $1$ approximation for $\boldsymbol{c}^{\star}$.

Continuing in this manner, we evaluate $\widehat{f}(\boldsymbol{c}^{(1)})$ via simulation and solve for $\boldsymbol{\tau}^{(2)}$ in $\widehat{f}(\boldsymbol{c}^{(1)})=\widetilde{f}(\boldsymbol{c}^{(1)},\boldsymbol{\tau}^{(2)})$ to obtain our iteration $2$ approximation $\widetilde{f}^{(2)}(\cdot):=\widetilde{f}(\cdot,\boldsymbol{\tau}^{(2)})$ that intersects $f$ at $\boldsymbol{c}=\boldsymbol{c}^{(1)}$ in expectation; see Figure 3(b). Observe that, while selection of the functional form $\widetilde{f}$ requires fundamental insights into the stochastic network of interest, our approach does not rely solely on inaccurate queueing formulas, in contrast to purely analytic approximations, because simulation is used to evaluate the stochastic network.

In general, our iteration process consists of the sequence of approximation functions

[TABLE]

with corresponding maximisers

[TABLE]

where the tuning parameter $\boldsymbol{\tau}^{(n)}$ in each iteration is obtained by solving

[TABLE]

We then iterate (8) – (10) until the difference $||\boldsymbol{c}^{(n)}-\boldsymbol{c}^{(n-1)}||$ is sufficiently small or a maximum number of iterations $N$ has been reached.

Our general approach relies on the fact that $\widetilde{f}(\cdot,\boldsymbol{\tau}^{(n)})$ provides a good approximation for $f$ around $\boldsymbol{c}=\boldsymbol{c}^{(n-1)}$ and thus the algorithm is likely to move in a proper direction at each iteration. This also provides significant computational improvements over purely simulation-based methods: Instead of, for instance, running many expensive simulations to estimate the Jacobian (or Hessian) as in SA, our iterative functional-form approach essentially uses $\widetilde{f}$ (a function $\mathcal{C}\to\mathbb{R}$) to provide $\nabla\widetilde{f}(\boldsymbol{c})$ as an approximation for the gradients $\nabla f(\boldsymbol{c})$, which requires only a single evaluation of $\widehat{f}$ per iteration.

Following the two-phase approach in [13], when the final solution $\boldsymbol{c}^{(n)}$ from our our general first-phase approximation approach needs to be further refined, then $\boldsymbol{c}^{(n)}$ can be used as a high-quality starting point for a simulation-based optimisation approach with guaranteed convergence properties, such as SA. Such a second-phase solution may provide improved accuracy over our general first-phase approach above, but it also incurs additional computational costs. However, by first obtaining a near-optimal solution using our fast iterative functional-form approach, the more expensive purely simulation-based approach requires far fewer iterations to find the optimal solution, thus significantly reducing the overall computational costs over exclusively using simulation-based optimisation [13]. The solution of our first-phase approximation approach therefore can be used directly at finer time scales or indirectly at coarser time scales as a high-quality initial starting point for a preferred simulation-based algorithm.

3.2 Application to networks of stochastic loss systems

We aim to solve the resource optimisation problem (3) by determining the server capacity $\boldsymbol{c}\in{\mathbb{Z}}_{+}^{L}$ that maximises objective function (1) or (2). In most cases of interest, we are unable to do so directly because these objective functions must be evaluated via simulation. The costs $\boldsymbol{{\theta}}$, rewards $\boldsymbol{{\omega}}$ and external arrival rates $(\lambda_{r})_{r\in\mathcal{R}}$ are assumed to be known constant vectors, with each $\lambda_{r}$ independent of the resource allocation $\boldsymbol{c}$. Recall that $p_{r,i}(\boldsymbol{c})$ denotes the proportion of class-$r$ customers arriving at station $l(r,i)$ (i.e., $i$-th station of their path $\psi_{r}$) that are then lost at $l(r,i)$ as a function of $\boldsymbol{c}$. The function $p_{r,i}(\boldsymbol{c})$ depends on the long-term behaviour of the underlying stochastic process and is not explicitly known in general.

The complexity of $f$ in (1) or (2) is contained in the set of functions $\{p_{r,i}\,:\,r\in\mathcal{R},\,i\in\mathcal{L}\}$, for which no closed-form expression is known. We therefore propose to replace each function $p_{r,i}$ with an approximation $\widetilde{p}(\cdot,\tau_{r,i})$ that follows a functional form which incorporates features that can be reasonably believed to be present in $p_{r,i}$. Here $\tau_{r,i}$ is an unknown parameter of the functional form indexed by class $r$ and station $i$. Define $p_{r,0}(\boldsymbol{c})\equiv 0$. For networks of SLSs with objective function (1) and $p_{r,i-1}(\boldsymbol{c})<1$, it is reasonable to expect that: (i) $p_{r,i}(\boldsymbol{c})$ is strictly decreasing in $c_{l(r,i)}$; (ii) $p_{r,i}(\boldsymbol{c})=1$ when $c_{l(r,i)}=0$; and (iii) $p_{r,i}(\boldsymbol{c})\to 0$ as $c_{l(r,i)}\to\infty$. Expectations (i)–(iii) similarly hold for networks with objective function (2) and $p_{r,i-1}(\boldsymbol{c})>0$. By definition, the case of $p_{r,i-1}(\boldsymbol{c})=1$ for objective (1) implies $p_{r,i}(\boldsymbol{c})=0$ for all $c_{l(r,i)}$ (given no class $r$ arrivals to station $l(r,i)$); similarly, by definition, the case of $p_{r,i-1}(\boldsymbol{c})=0$ for objective (2) implies $p_{r,i}(\boldsymbol{c})=0$ for all $c_{l(r,i)}$. In these instances, changes to $c_{l(r,i)}$ do not affect the objective function via $p_{r,i}(\boldsymbol{c})$. We therefore focus on cases where $c_{l(r,i)}$ influences the objective function by finding an approximation $\widetilde{p}(\cdot,\tau_{r,i})$ for $p_{r,i}$ that satisfies properties (i)–(iii). Since these properties depend only on the resource allocation $c_{l(r,i)}$ at the relevant station $l(r,i)$, we take $\widetilde{p}(\cdot,\tau_{r,i})$ to be a univariate function of $c_{l(r,i)}$ (disregarding other components of $\boldsymbol{c}$). We observe that cumulative probability distribution functions for r.v.s with non-negative support form a class of functions with exactly these properties. Many such candidate functions were tested, including but not limited to (in the one-dimensional case) $\exp(-(c\tau)^{k})$ and $1/(1+(c\tau)^{k})$ for a wide range of $k$, as well as $(\tau^{c}/c!)\left(\sum_{i=0}^{c}\tau^{i}/i!\right)^{-1}$. Due to space constraints, we do not present results on all of these forms but instead focus on the functional form found to work best for the networks of SLSs of interest:

[TABLE]

As we will see in Section 4, the above choice indeed performs very well. Although the capacity values in (3) are discrete, the application of our optimisation approach (described below) using (11) requires continuous functions of $\boldsymbol{c}$. We therefore apply the following assumption during the optimisation procedure.

Assumption 1**.**

When a customer arrives to a station with (non-integer) capacity $c$ and finds $\lfloor c\rfloor+1$ customers, the customer is lost. If an arrival encounters $\lfloor c\rfloor$ customers, the newly arriving customer is accepted with probability $c-\lfloor c\rfloor$, and lost otherwise. When an arrival finds less than or equal to $\lfloor c\rfloor-1$ customers, the arrival is accepted.

This assumption results in a continuous relaxation of the objective functions that matches (1) or (2) at integer points. Whenever our (approximate) optimal solution is non-integer, it is rounded to the closest integer. We note that it is not uncommon to use a continuous relaxation of a discrete problem in this way (see, e.g., [2, Section 3.1]).

Our algorithm is formulated following (8) – (10) and formally presented in Algorithm 1 below. Let $\widehat{p}_{r,i}(\boldsymbol{c})$ denote an unbiased estimator of $p_{r,i}(\boldsymbol{c})$ obtained from simulation for an allocation $\boldsymbol{c}$. Given an initial guess (possibly at random) for the optimal capacity allocation $\boldsymbol{c}^{(0)}$, we simulate the system to obtain the loss probabilities $\widehat{p}_{r,i}(\boldsymbol{c}^{(0)})$ for each $(r,i)$ pair of class and station indices. Using this we can determine the parameters $\tau_{r,i}^{(1)}$ by solving the implicit equations

[TABLE]

Namely, we want to calibrate our approximation $\widetilde{p}(\cdot,\tau_{r,i}^{(1)})$ by choosing $\tau_{r,i}^{(1)}$ such that the predicted loss probability at $\boldsymbol{c}^{(0)}$ matches the observed loss probability. After some algebraic manipulation we have

[TABLE]

Once these $\tau_{r,i}^{(1)}$ are determined, for Model I we can obtain a new approximation to the optimal resource allocation problem by evaluating the corresponding version of (1) and (3):

[TABLE]

with $\widetilde{p}(\cdot,\tau_{r,i})$ given by (11). Solving (13), which can be done efficiently using standard optimisation techniques, renders a new resource allocation, which is used as the basis for the next step. Similarly, for Model II, once the collection of $\tau_{r,i}^{(1)}$ are determined we can obtain a new approximation to the optimal resource allocation problem by evaluating the corresponding version of (2) and (3):

[TABLE]

with $\widetilde{p}(\cdot,\tau_{r,i})$ given by (11). Solving (14) yields a new resource allocation, which is used as the basis for the next step. Iteratively applying these steps under either Model I or Model II for a desired level of accuracy $\epsilon>0$ or for a maximum number of steps $N\in{\mathbb{Z}}^{+}$ leads to Algorithm 1. This algorithm has (at least) two advantages over other optimisation techniques: (i) it explicitly incorporates the functional form of the objective function (1) or (2); and (ii) the functional-form approximation $\widetilde{p}(c,\,\tau_{r,i})=\exp\left(-(c\,\tau_{r,i})^{2}\right)$ given in (11) explicitly utilises known fundamental features of the intractable quantities $p_{r,i}(\boldsymbol{c})$.

Assuming a maximum number of servers available for allocation, (13) and (14) define a mapping for which the existence of a fixed point can be established using Brouwer’s theorem. Showing that this is also a contraction mapping (and therefore possesses a unique fixed point) can be very difficult because doing so relies on using properties of $p_{r,i}(\boldsymbol{c})$ which are not known in general. Even in the simple case of a single station with one class, Poisson arrivals and exponential service times, showing that the mapping is a contraction involves handling highly complicated expressions in terms of the derivatives (with respect to $c$) of the function

[TABLE]

For these reasons, we leave for future work the problem of showing that Algorithm 1 provably converges, noting that this currently remains an (highly challenging) open problem. By the inherent nature of the problem at hand, the efficiency of any implementation of Algorithm 1 will be dominated by the speed with which the simulation samples $\widehat{p}(\boldsymbol{c})$ can be generated. Our algorithm takes as input a vector of size $L$ and then in each iteration, after generating the set of estimates $\widehat{p}_{r,i}(\boldsymbol{c})$, the quantities $\tau_{r,i}$ must be computed for each class-station pair. Relative to obtaining the estimates $\widehat{p}_{r,i}$, computation of $\tau_{r,i}$ is very efficient, as is solving the optimisation problem in Step 6.

4 Computational experiments

In this section we perform extensive computational experiments to investigate and compare the performance of the core iterative algorithm of our first-phase approximation approach, reporting and discussing experimental results for four example networks of SLSs: (a) a Markovian tandem network; (b) a tandem network with cross-routes; (c) a ring network; and (d) a canonical network. These network topologies are illustrated in Figure 4. The goal of the experiments is to gauge and compare the accuracy and efficiency of our core algorithm, presented in the previous section, against benchmark algorithms and results. To our knowledge there are no examples of other algorithms tailored specifically to optimise the networks of SLSs of interest herein. For the Markovian tandem network example, we are able to derive explicit matrix-analytic expressions for the objective function (see Appendix A) that provide a benchmark for accuracy. Such expressions are not available for the other three example topologies of general networks of SLSs. For this reason, we compare the accuracy and efficiency of our core iterative algorithm against three benchmark algorithms: (i) a fine-tuned SA implementation (see Appendix B); (ii) SO (as implemented in Matlab 2020b with the surrogateopt function); and (iii) BO (as implemented in Matlab 2020b with the bayesopt function).

In each of the experiments, all the algorithms are based on mean-value Monte Carlo simulations of the models. Given a capacity allocation $\boldsymbol{c}$, the system is simulated from an initial empty state to obtain $\widehat{p}_{r,i}(\boldsymbol{c})$. Time is then batched into prespecified periods of length $\widetilde{T}$. At the end of each time period, approximate $95$% confidence interval widths for all $\widehat{p}_{r,i}(\boldsymbol{c})$ are computed using the score method for Bernoulli-type data (see e.g., [26, Example 6.16]) as follows:

[TABLE]

where $A^{(n)}_{r,i}$ is the number of arrivals to station-class pair $(r,i)$ in time period $n$. These approximate confidence intervals are appropriate since each arrival to a station can be viewed as a Bernoulli trial with parameter representing the probability that the arrival is accepted for processing. When all of these confidence interval widths fall below a specified level of $0.01$ or a maximum simulation clock is reached, the simulation ends and a value is returned. In our experiments, we follow this procedure $10$ times and then return the average of these values as our estimates $\widehat{p}_{r,i}(\boldsymbol{c})$. After extensive tuning, we found that setting the maximum simulation clock to $100$ provided a good balance between speed and accuracy (so that our numerous experiments could be completed within a reasonable time-frame). Some capacity allocations result in the true loss proportion being extremely close to [math], which can result in an inaccurate estimate of [math] being returned with high probability. Since our model does not allow for loss proportions with value [math] unless capacity is also [math], we set the estimate $\widehat{p}_{r,i}(\boldsymbol{c})$ to $10^{-6}$ whenever a loss proportion of [math] is estimated with positive capacity. There are a variety of other ways that these estimates could be obtained, which could also be used within our optimisation framework, but this method was chosen for its simplicity.

The computational experiments are performed in Matlab 2020b and the numerical optimisation package to implement Step 6 of Algorithm 1 is the interior point algorithm built into this platform. The experiments were conducted on an Intel Core i7-3770 3.4GHz workstation running Windows 10 with 16GB of RAM. Code for these experiments is available online at https://github.com/bpatch/networks-of-stochastic-loss-systems, where we also provide data files containing the results of completed experiments. All experiments are conducted with a fixed seed (included in the code) so that they can be easily replicated. Due to the large number of parameters needed to generate all of the scenarios we consider, it is impractical to provide all of these parameters in table form here and therefore the interested reader is referred to the data files provided online. In all of the box plots in this paper, outliers are defined as points that are greater than $q_{75}+w(q_{75}-q_{25})$ or less than $q_{25}-w(q_{75}-q_{25})$, where $q_{n}$ is the $n$-th percentile and $w$ is approximately $2.7$ times the sample standard deviation (which corresponds to 99.3% coverage, if the data are normally distributed).

4.1 Markovian tandem networks

Tandem Experiment 1 (illustrative). We start with an experiment that illustrates how our algorithm and benchmarks evolve as they search for the optimal solution in a specific problem instance. Consider the example from Section 2, displayed in Figure 2, with arrival rate $\lambda_{1}=16$, capacity costs ${\theta}_{1}=0.2$ and ${\theta}_{2}=0.3$, and suppose that the service requirements of jobs are exponentially distributed with rate $0.8$ at the first station and $0.6$ at the second station. We consider both variants of our model: For Model I, let rewards for successful service completions be ${\omega}_{1}=1.9$; and for Model II, let ${\omega}_{1,1}=1$ and ${\omega}_{1,2}=0.9$.

For this simple network, it is possible to explicitly evaluate the objection functions (4) and (5) using matrix-analytic methods (MAMs) and then search for the optimal solution in a brute-force manner. Although this approach does not extend to more complex settings, we present these derivations in Appendix A to enable an exact evaluation of the accuracy of Algorithm 1 compared to the SA, SO and BO benchmarks. Computations from the MAM approach for Model I provide the true maximum $f(\boldsymbol{c}^{\star})=13.4975$ at $\boldsymbol{c}^{\star}=(26,32)$. On the other hand, for Model II, the true maximum is $f(\boldsymbol{c}^{\star})=10.2049$ at $\boldsymbol{c}^{\star}=(26,0)$. Figure 5 displays trajectories of the values returned by Algorithm 1, SA, SO and BO as these algorithms search for the optimal function values $13.4975$ and $10.2049$, where the true maximum is displayed as a dashed line and the mean of the trajectories is displayed as a solid blue line, together with a light blue shading that encompasses all trajectories; and each of the trajectories displayed in grey. The displayed trajectories are from noisy evaluations of the objective function (and thus it is possible that they exceed the true maximum). Each trajectory starts from a random initial condition (which is held constant across approaches). Our key observations for the results in this figure, which will be remarked upon soon, are as follows.

Within the allocated number of iterations, all four methods are able to greatly improve upon the objective function value from that of the initial random capacity allocation. Specifically, sample paths of Algorithm 1 for both models achieve on average 95% of the potential improvement that is possible from the given initial conditions; one sample path under Model II exhibits an improvement of only 56%, suggesting this this trajectory has become stuck in a local solution. The objective function values returned by SO and BO consistently achieve 100% of the potential improvement. The objective function values returned by SA were also highly competitive, respectively achieving 96.5% and 99.7% of the potential improvements for Model 1 and Model II.

• -

Trajectories of Algorithm 1 typically make an immediate large jump towards the optimal function value (for Model I, eight of the trajectories make a gain of more than $88\%$ by $n{=}2$; and for Model II, eight of the trajectories make a gain of more than $99\%$ by $n{=}2$), much more quickly than the other methods, and then they experience minor oscillations (less than $1$ and $0.5$ after $n{=}3$ for Model I and Model II, respectively).

• -

Trajectories of SA and BO make consistent but small movements towards the optimal function value, whereas SO experiences mid-sized jumps towards the optimal function value separated by long periods of no movement (sometimes for more than $10$ iterations).

We follow up on our observations of the solution quality from the four algorithms with an investigation into the computing power needed to obtain these solutions. Figure 6 displays the cumulative CPU time required to generate the objective values displayed in Figure 5 for each trajectory (grey) and on average (blue). Our key observations for the results in this figure are as follows.

On average, to reach an objective function value within $1\%$ of the optimal value returned, the CPU times for SA (Model I: 3.6060 seconds; Model II: 6.8583 seconds) and for BO (Model I: 10.2423 seconds; Model II: 9.0879 seconds) are substantially greater than the CPU times for Algorithm 1 (Model I: 0.2321 seconds; Model II: 0.1493 seconds) and for SO (Model I: 1.29 seconds; Model II: 1.132 seconds), with those for SO significantly greater than those for Algorithm 1. Specifically, in comparison with Algorithm 1 on Model I and Model II, SA requires a multplicative factor of $15.5$ and $45.9$ more CPU time, BO requires a multplicative factor of $44.1$ and $60.9$ more CPU time, and SO requires a multplicative factor of $5.6$ and $7.6$ more CPU time.

• -

The best solution for Model I returned by Algorithm 1, SA, SO and BO in these experiments were $(23,30)$, $(25,30)$, $(24,30)$ and $(26,30)$, respectively. For Model II the best solution for Algorithm 1, SA and BO was $(26,0)$, while the best solution returned by SO was $(27,0)$.

One important takeaway from this experiment is that all four algorithms are capable of reliably generating solutions which are close or equal to the optimal solution with very small optimality gaps. The optimal function values achieved by SO and BO were the best, but were followed closely by Algorithm 1 and SA. Our algorithm was able to do so in the most efficient manner, and at speeds which were more than $5$-$7$ times faster than the next competitive method, namely SO. This experiment also suggests that SA and BO are vastly less efficient. For Model II, we observed that our algorithm can become trapped in a local solution, something that did not occur for the benchmarks in this experiment. At this point, based on the initial simple experiment, it appears that our algorithm may offer vast speed improvements but may suffer slightly in terms of accuracy, suggesting that SO may be best when accuracy is most important and Algorithm 1 is best when efficiency is most important. On the other hand, as previously discussed, the solution from our algorithm can be used as a high-quality starting point for SO and then benefit from both accuracy and efficiency.

Tandem Experiment 2 (scenarios). To investigate the accuracy of our algorithms more generally, we now consider a broader range of randomly chosen parameters for each model variant of the Markovian tandem network displayed in Figure 2. For Model I, we study $20$ different randomly sampled scenarios each corresponding to different parameter values, where the rewards, service rates, arrival rate and capacity costs are all sampled uniformly at random from $(0,5)$, $(0,3)$, $(0,20)$ and $(0.1,0.3)$, respectively. Similarly, for Model II, we sampled $20$ additional scenarios according to the same scheme, except with the rewards and capacity costs sampled uniformly at random from $(0,3)$ and $(0,0.4)$, respectively. To ensure the MAM-based computations for finding the exact solution (Appendix A) did not become too burdensome, we only experimented with scenarios where the optimal value as determined by these exact methods is an element of $[1,60)^{2}$. We conducted $10$ computational experiments for each scenario where each experiment was performed with a different randomly selected initial condition (which is held constant across methods).

Focusing first on accuracy, panels (a) and (b) of Figure 8 display the differences between the function values returned by each method and the values returned by Algorithm 1 (including the function values corresponding to the random initial conditions). Values greater (less) than 0 exhibit better (worse) performance compared to Algorithm 1. For these experiments, the reported objective function values are exact evaluations of the solutions outputted by the algorithms (using the methodology in Appendix A). The values given by the MAM method are the true optimal values. Since the experiment corresponding to each scenario is replicated 10 times with different initial conditions, each box plot summarises 200 data points. Our key observations for the results in these panels are as follows.

For Model I, Algorithm 1 improved upon the the objective function values associated with the random initial conditions by a median amount of $13.2222$, an improvement which is approximately (within 1.5 units) matched by the other methods (and exactly equal for SA). This improvement falls short of the best possible improvement in the median (given by MAM as $14.7344$) by $1.5122$ units. For SO and BO this shortfall is smaller, at only $0.0527$ and $0.5357$, respectively. This means that the values returned by MAM, SO and BO are better than the values returned by Algorithm 1 in more than $50\%$ of the experiments (for this specific simple instance of Model I), sometimes by up to 6.5855 units for SO and $4.3250$ units for BO (excluding outliers).

• -

These differences in values (panel (a)) are reflected by the Euclidean distances from the optimal solution, as displayed in panel (b). All of the algorithms move the solution much closer to the optimal solution, relative to the random initial condition, with median distances from the optimal solution of $3.8643$, $15.0830$, $1$ and $1.4142$ for Algorithm 1, SA, SO and BO, respectively. For Algorithm 1, SO and BO, these distances are also highly concentrated with interquartile ranges of $8.7097$, $2.2361$ and $3.8028$, respectively. The distances from the optimal solution are more spread out for SA, with an interquartile range of $39.7744$.

• -

For Model II, Algorithm 1 improved upon the the objective function values associated with the random initial conditions by a median amount of $19.6475$, an improvement which is approximately (within 0.3 units) matched by the other methods. This is close to optimal, lying within $0.1158$ units from the maximal median improvement (as given by MAM). Furthermore, since the interquartile ranges of the difference between MAM, SO, and BO with Algorithm 1 are all less than $0.35$, such results are achieved with very litte variability for these algorithms. The interquartile range of the difference in objective values is higher for SA at $3.0077$.

• -

These differences in values (panel (d)) are again reflected in the Euclidean distances from the optimal solution, as displayed in panel (e). The median distance from the optimal solution for the random initial conditions of $75.4307$ is reduced down to $2$ by Algorithm 1 and 1 for SO and BO. The reduction is not as great for SA, achieving a median of difference of $11.6110$. These reductions are also achieved reliably for Algorithm 1, SO and BO, with interquartile ranges of $2.2839$, $1$ and $1.2071$, respectively. For SA, the reductions are less reliably achieved with an interquartile range of $31.6012$.

We follow up on our observations of the solution accuracy from the four algorithms with an investigation into the computing power needed to obtain these solutions, where panels (c) and (f) of Figure 8 report the CPU times needed to achieve values that are within 1% of the values reported in panels (a) and (d) of the figure. The reason we allow a 1% gap is that SA, SO and BO may spend large amounts of CPU time making minor improvements in the objective value. Allowing this gap ensures that these algorithms are not penalised for such minimal refinements, thus providing a favourable case for these algorithms. Our key observations for the results in these panels are as follows.

Algorithm 1 is significantly more efficient than the benchmarks for both models. The median CPU times (seconds) for Algorithm 1 were $0.0978$ and $0.0783$ for Model 1 and Model 2, respectively. This is much more efficient than the corresponding times for SA and SO, which across both models were in the range $0.5756$ to $1.2012$, and vastly more efficient than BO, which had median CPU times of $21.275$ and $16.1226$ for Model 1 and Model 2, respectively.

Based on the observations from the results for the simple network in this subsection, we remark that there is evidence our algorithm and the benchmarks are all capable of substantially improving the performance of networks of SLSs. SA appears to be the weakest algorithm in terms of both accuracy and efficiency (although it may be improved with further tuning, we do not consider this since the other approaches do not benefit from tuning and it is not always possible to provide a tuning which works well for all perturbations of the model parameters and initial conditions). There is strong evidence that SO and BO reliably output a solution which is very close to the optimal solution. Our algorithm also outputs a solution which is close to the optimal solution, and it manages to do so much more efficiently than SA, SO or BO. Indeed, there is evidence that Algorithm 1 is typically (in terms of median performance) $10$ times faster than SA and SO, and $100$ times faster than BO. At this stage of our experiments, the results thus far suggest that our algorithm should be used when one needs a near optimal solution in a very efficient manner. When a greater computational budget is available, our algorithm can be used to provide a high-quality starting point for SO, BO or SA following [13] with the combination of algorithms being run to achieve asymptotic optimality. We next conduct further experiments to examine whether these findings continue to hold true for more complex networks.

4.2 Criss-cross networks with Markov-modulated-Poisson-process arrivals

To extend our investigation of the performance of our algorithm and the benchmarks into a setting with many stations and multiple customer classes, we consider a criss-cross network as displayed in Figure 4(b). This linear network consists of ten stations $\mathcal{L}=[10]$ and two classes where the customers of one class follow the path $\psi_{1}=(1,2,\dots,10)$ and the customers of the other class follow the path $\psi_{2}=(10,9,\dots,1)$. Customers arrive to the network according to a Markov modulated Poisson process (MMPP) with three background states; arrivals occur at a rate of $10$, $20$ or $30$ depending on the state of the background process. The background process jumps between states at a constant rate, and upon jumping moves to one of the other two states with equal probability. Arriving customers are equally likely to be of class $1$ or class $2$. The service requirements at each station are exponentially distributed. Since no analytical results are known for this network, we must rely on our benchmarks to test the accuracy of Algorithm 1.

Criss-cross Experiment 1 (illustrative). We again begin with an experiment that illustrates the trajectories of our algorithm and benchmarks in terms of accuracy and efficiency. Consider the system above where the background process jumps between states at unit rate and the service rate at each station is $1.2$. In the case of Model I, the reward for successfully processing a class 1 customer is ${\omega}_{1}=9.8628$ and the reward for successfully processing a class 2 customer is ${\omega}_{2}=8.7339$ (where these parameter values were chosen, prior to conducting the experiment, uniformly at random from $(0,20)$). The capacity costs were sampled uniformly at random from $(0.1,0.3)$, as given in Table 1. In the case of Model II, rewards for successfully processing customers are uniformly sampled (per class, per station) from $(0,10)$ and the costs ${\boldsymbol{{\theta}}}$ of capacity are uniformly sampled (per station) from $(0.1,\,0.3)$, as reported in Table 2.

For each model of the networks just described, the initial capacity allocations were sampled uniformly at random from $(1,2,\dots,100)^{10}$. Each of these capacity allocations was used to initialise our algorithm and the benchmarks. In Figure 9, the improvement in the objective value from running the $3$ methods is displayed. Beyond explicitly displaying trajectories of the algorithms, the figure also displays the average of the sample paths and a shaded region containing all of the sample paths. Our key observations for the results in this figure are as follows.

For Model I, Algorithm 1, SO and BO output comparable best objective values on average ($113.4697$, $102.4320$ and $104.4699$, respectively), which are substantially better than the average initial value (1.5387), with Algorithm 1 demonstrating some advantage over the other two. The spread of the output values closely surround these average values (with the range of the final output values being $29.7253$, $32.3170$, $43.5180$ and $21.2739$ for Algorithm 1, SA, SO and BO, respectively).

• -

For Model II, Algorithm 1 ($161.5652$) and SO ($147.4500$) output somewhat similar best objective values on average, which are substantially better than SA ($108.6893$), BO ($111.0690$) and the random initial function value ($31.4809$). Algorithm 1 again appears to have a sample path that has become trapped in a local solution.

• -

For both models, SA improves the objective value relative to the random initial solution, but also appears to require many more iterations before it may be able to achieve the performance of our algorithm or SO (which in turn requires vastly more CPU time than SO is allocated, as will be seen next).

We follow up on our above observations on the solution quality obtained by our algorithm and the benchmarks with an investigation into the computing power needed to obtain these solutions. In Figure 10, the cumulative CPU time needed to generate the objective values displayed in Figure 9 are displayed for each sample path and on average. Our key observations for the results in this figure are as follows.

BO uses vastly more CPU time and SA uses substantially more CPU time to generate iterates than does Algorithm 1 or SO (note the different scale of the axes).

• -

For Model I, the amount of CPU time needed by Algorithm 1 and SO is comparable. However, after an initial substantial improvement in value, SO endures a costly period of slow increase towards a best solution that is comparable to that given by Algorithm 1.

• -

Similarly, for Model II, Algorithm 1 very quickly (in terms of iterations and CPU usage) reaches a best solution comparable to the best solution of SO (as discussed above). SO experiences a slower steady increase towards its best solution, requiring substantial CPU time.

This set of experiments has illustrated that our algorithm can more efficiently obtain solutions comparable (and for this network topology they appear to be slightly better) to our benchmarks. Further evidence has been provided that SO is the superior benchmark. The potential for Algorithm 1 to become trapped in a local solution has appeared again. We now investigate whether these findings hold for a broader range of parameter values.

Criss-cross Experiment 2 (scenarios). For each model we sampled 20 scenarios. In each scenario for both models, the arrival rates corresponding to the three background states are sampled uniformly from $(0,20)$, $(0,40)$, and $(0,60)$. Rewards are sampled from $(0,10$), service rates (per station) from $(0,2.4)$, capacity costs from $(0.1,0.3)$, and the jump rate of the background chain from $(0,2)$, all uniformly at random.

To investigate the accuracy of our algorithms, as before, panels (a) and (d) of Figure 12 report the differences between the objective function values returned by each method and the values returned by Algorithm 1. All algorithms again improve substantially upon the random initial solutions. In contrast to the experiments for tandem networks however (as displayed in Figure 8), for both models in this case Algorithm 1 clearly demonstrates higher accuracy relative to SA and BO (which have medians that are $45.3431$ and $17.7394$ lower than the median of Algorithm 1). For both models, Algorithm 1 and SO have essentially the same median objective function values (for Model I the difference is in fact 0 and for Model II SO is lower by only $0.3598$). The superiority of Algorithm 1 and SO in terms of accuracy is also evident in panels (b) and (e), where it is shown that they return solutions which are typically much closer to the solution associated with the best observed function value across all experiments for any particular scenario. Finally, panels (c) and (f) of Figure 12 report the CPU times needed to achieve values which are within 1% of the values underlying panels (a) and (d). For both models, it can be seen that SA and BO are substantially less efficient than Algorithm 1 and SO. For Model I, in many cases SO and Algorithm 1 appear to have comparable (near 0 CPU time) efficiency, but there are a clear set of cases where SO is much less efficient (taking up to $46.8400$ second and having a 75$$th percentile of $23.6700$ seconds). For Model II, it is clear that Algorithm 1 is more efficient, except possibly in a small handful of cases. The key general observations made for tandem networks continue to hold true for criss-cross networks. Specifically, our algorithm and SO return the best solutions, and our algorithm typically does so much more efficiently than SO. The SA and BO benchmarks certainly improve the objective value but are much less efficient than our approach or SO.

4.3 Ring networks with MMPP arrivals

In this subsection we consider a network consisting of a set of $10$ stations arranged in a ring, as displayed in Figure 4(c). Associated with arrivals at each station $i$ is a route $\psi_{i}=(i,i+1,\dots,1,2,\dots,i-1)$ utilising all other stations in the network sequentially. This means that the system consists of $1000$ types of job when the jobs are classified by the combination of their route and the station at which they are being processed. We consider a parameterisation of the model where parameter values are randomly chosen following the same scheme as for criss-cross networks in the previous subsection (e.g., rewards are sampled from $(0,10)$).

For the Model I and Model II variants of the ring network just described, $10$ initial capacity allocations were sampled uniformly at random from $(1,2,\dots,100)^{10}$ for each scenario. Each of these $10$ capacity allocations was used to initialise our algorithm and the benchmarks for the ring network replications of the experiments previously performed for the tandem (Figure 8) and criss-cross (Figure 12) networks. Since we are optimising over a $10$ dimensional set, the results from the SA algorithm were not competitive in this case and therefore we only compare our algorithm to the SO and BO benchmarks. In Figure 14, we see that our earlier findings generalise to these more complex systems; namely our algorithm returns comparable objective values to the best performing benchmark, which is SO. For Model I, SO has a median performance that is only $0.3239$ greater than Algorithm 1, but it can also be up to $40.9869$ units better or up to $26.9052$ units worse. Similarly, for Model II, SO has a median performance that is $11.0364$ units lower than Algorithm 1, but it can also be up to $22.6228$ units better or up to $109.5804$ units worse. For Model I, the solutions returned by SO and Algorithm 1 have highly similar performance (in panel (b) the medians and interquartile ranges for the distance to the best observed solution are within 3.7929 units and 9.2617 units of each other, respectively). In the case of Model II, however, Algorithm 1 returns solutions that are much closer to the best observed solution (Algorithm 1: $5.3852$, SO: $57.4540$, BO: $68.0955$). Importantly, our algorithm achieves these results much more efficiently than SO and BO with median CPU usage for Model I at $3.6807$, $65.9950$ and $279.0879$ for Algorithm 1, SO and BO, respectively; and median CPU usage for Model II at $5.4399$, $93.0550$ and $143.2500$ for Algorithm 1, SO and BO, respectively. The BO benchmark also returns good objective values, but these are noticeably inferior to those provided by Algorithm 1 or SO. In summary, the overwhelming superiority of Algorithm 1 in terms of efficiency is again evident—our approach is approximately at least $10$ times faster than SO or BO. At the same time, Algorithm 1 returns excellent objective function values that are competitive with the benchmarks for Model I and are usually superior to the benchmarks for Model II.

4.4 Canonical networks with renewal arrivals and two-state Coxian distributed service requirements

In this subsection we evaluate the performance of Algorithm 1 on a six-station, twelve-class system, as displayed in Figure 4(c). Here $\boldsymbol{c}=(c_{1},\dots,c_{6})$, $\mathcal{L}=[6]$, $\mathcal{R}=[12]$, and the customer classes follow the network paths $\psi_{1}=(1)$, $\psi_{2}=(1,3)$, $\psi_{3}=(1,3,5)$, $\psi_{4}=(1,4)$, $\psi_{5}=(1,4,5)$, $\psi_{6}=(1,4,6)$, $\psi_{7}=(2)$, $\psi_{8}=(2,4)$, $\psi_{9}=(2,4,6)$, $\psi_{10}=(2,3)$, $\psi_{11}=(2,3,6)$, and $\psi_{12}=(2,3,5)$. Customers of each class arrive to the first per-path station according to a per-class renewal process; the service time at each station is the same for all classes, while the service times differ across the various stations.

In all of the experiments, the interarrival times and service requirements are generated from two-stage Coxian distributions with parameters set and scaled to match a collection of coefficients of variation (CoVs). We consider $20$ randomly chosen (unique) scenarios where the interarrival times have CoVs randomly sampled from $\{0.75,2,3.25\}$, service times at stations $1$ and $2$ have CoVs randomly sampled from $\{2,3.75,5\}$, service times at stations $3$ and $4$ have CoVs randomly sampled from $\{1.5,3,4.5\}$, and service times at stations $5$ and $6$ have CoVs equal to $1.5$. The mean service requirements were randomly sampled (per station) from $(0,0.2)$. The interarrival times for customers of class $i$ are sampled as \frac{1}{2\lambda_{i}}\big{(}\,E_{1}+\frac{1}{q_{i}}E_{2}\,B_{i}\big{)}, where $E_{1}$ and $E_{2}$ are independent unit-mean exponentially distributed r.v.s, $q_{i}\in(0,1)$ is a real number, and $B_{i}$ is a Bernoulli r.v. with parameter $q_{i}$; the parameter $\lambda_{i}$ is sampled from $(0,0.2)$. Then, the parameter $q_{i}$ is set such that the CoV of the interarrival time matches the desired value for the scenario being studied. The service times are specified on a per-station basis, with the service times of all customers at station $i$ sampled as \frac{1}{2\mu_{i}}\big{(}\,E_{1}+\frac{1}{q_{i}}E_{2}\,B\big{)}, where $E_{1}$, $E_{2}$, $B_{i}$, and $q_{i}$ are as before (chosen such that the desired CoV is achieved). The capacity costs ${\theta}_{i}$ are sampled from $(0.1,0.3)$ and the rewards (${\omega}_{r}$, or ${\omega}_{r,i}$) are sampled from $(0,50)$.

For the canonical network just described, an initial capacity allocation was sampled from $(1,2,\dots,100)^{6}$ for each scenario. These capacity allocations were used to initialise trajectories of our algorithm and the SO and BO benchmarks for both model variants. The results from the SA algorithm were again not competitive; hence we only compare our algorithm to the SO and BO benchmarks. We performed similar experiments to those which resulted in Figure 8, Figure 12, and Figure 14. The results from these experiments, now in the context of the canonical networks of this section, are displayed in Figure 16. As with the other network types, we see for both models that our algorithm is able to find its best approximate solution much faster than SO and BO. For Model I, we have that Algorithm 1, SO and BO respectively exhibit median CPU times (displayed in panel (c)) of $3.6967$, $52.5600$ and $104.0930$; and for Model II, we have that Algorithm 1, SO and BO respectively exhibit median CPU times (displayed in panel (f)) of $3.7164$, $98.3150$ and $114.8809$. All three algorithms obtain massive value improvements over the initial random capacity allocation (median more than $8.2910$ for Model I and median more than $56.3184$ for Model II). As shown in panels (a) and (d), the objective function values obtained by our algorithm are very similar to those obtained by SO (median $0.7953$ units lower for Model I, and median $0.6316$ units lower for Model II) and typically slightly better than those obtained by BO (median $2.1895$ units higher for Model I, and median $4.3081$ units higher for Model II).

Algorithm 1 again demonstrates that it is able to obtain high-quality solutions much more efficiently than the benchmark approaches.

5 Conclusion

In this paper we introduced a general class of networks of SLSs, motivated by applications arising in a wide variety of areas including cellular wireless networks, cloud edge computing, computer and network spare-part supply chains, and emergency services—just to name a few. Since revenues are lost due to congestion and loss of customers and since costs are incurred due to allocating servers to network stations, determining the optimal allocation is critically important and requires a careful balance between network performance and service costs. We therefore investigated the optimization of resource allocation in networks of SLSs with the goal of balancing the trade-off between mitigating congestion by increasing service capacity and maintaining low costs for the service capacity provided. Two model variants of such networks are considered as being of particular interest. The first model is applicable to situations where customers traverse a sequence of stations and may be lost at any one of them, e.g., cellular wireless networks or cloud edge computing, whereas the second model is applicable to situations where customers are redirected to new stations until they find one with sufficient capacity, e.g., computer spare-part supply chains or emergency services.

Given the lack of analytical results and the computational burden of simulation-based methods, we proposed a novel hybrid functional-form approach for the optimal allocation of resources SLSs that combines the speed of an analytical approach with the accuracy of simulation-based optimisation. Our core iterative algorithm replaces the computationally expensive gradient estimation in simulation optimisation with a closed-form analytical approximation that is calibrated using a single simple simulation run. Extensive computational experiments on complex networks within the context of four relevant example network topologies—two-station Markovian networks, criss-cross networks, ring networks, canonical networks—show that our approach consistently and with low variability yields near-optimal solutions for both network model variants under these four network topologies, with objective function values that are comparable (and sometimes even superior) to those obtained using stochastic approximation, surrogate optimisation and Bayesian optimisation methods while requiring dramatically less computational effort. Hence, our general solution approach has several advantages: use at finer time scales to obtain near-optimal solutions, including situations that require dynamic reallocation of capacities over frequently time-varying periods (e.g., emergency services in a widespread disaster); use at coarser time scales as a high-quality initial starting point for a simulation-based optimization approach; use to address very large networks (e.g., $1000$s of stations or more) for which simulation-based optimization approaches can be infeasible due to their computational costs.

Our fundamental approach can be extended to more general network settings (e.g., incorporating re-trialling customers, priority classes or rerouting mechanisms). In order to do so, new functional-form approximations need to be formulated and evaluated. As we already considered two different model types, our study provides a solid foundation for these extensions and demonstrates how our general framework can be translated between models.

Acknowledgements. An earlier version of this work was completed while BP was supported by the Australian Research Council (ARC) through the ARC Centre of Excellence for the Mathematical and Statistical Frontiers (ACEMS) under grant number CE140100049, the Netherlands Organisation for Scientific Research (NWO) through the Gravitation project Networks under grant number 024.002.003, and an Australian Government Research Training Program (RTP) scholarship. This version was completed while BP was supported by NWO grant 639.033.413 and ARC grant DP200101281. This work was initiated while BP was visiting IBM Research in Yorktown Heights.

Appendix A Matrix derivations for Markovian tandem network

While evaluating the objective functions (1) and (2) in Section 2 for a single-station, single-class Markovian system (i.e., an Erlang-B system) is quite straightforward on a modern computer, this is not the case for more complex networks of SLSs. The power of the approach described in [33], based on matrix-analytic methods (MAMs), is that more general systems can be considered. In this section we use these methods to explicitly compute the objective functions for the relatively simple example provided in Section 2, so that it can be used to provide an exact evaluation of the accuracy of our general approach as well as provide comparisons with the various benchmark methods. More specifically, we explain how to use matrix-analytic methods to obtain explicit expressions for the reward terms in (4) and (5), i.e., for

[TABLE]

where $\kappa_{1}$ is the reward term in (4) and $\kappa_{1,1}$ and $\kappa_{1,2}$ are the reward terms in (5).

A.1 Model I

First consider Model I and let \big{(}N_{1}(t),\,t\in\mathbb{R}_{0}\big{)} count the number of customers that successfully enter the second station. The rate at which $N_{1}$ increases corresponds to $\kappa_{1}$; this is the quantity we wish to identify. Let $X_{1}(t)$ and $X_{2}(t)$ be the number of customers that are processed by the system at the first and second stations, respectively. Upon an arrival to the first station, as long as $X_{1}<c_{1}$, then $X_{1}$ jumps up by $1$. Upon a successful service completion at the first station, if $X_{2}<c_{2}$, then $X_{2}$ jumps up by $1$ and so too does $N_{1}$. It can be seen that the bivariate process $(X_{1},X_{2})$ is a background process for an encompassing MAP that experiences arrivals when customers enter the second station (see, e.g., [18, Section 2.3] for more details on MAPs). In order to find $\kappa_{1}$ explicitly, we need to carefully construct the structure of the Markov chain $(X_{1},X_{2})$, use its relationship to $N_{1}$, and then apply [30, Theorem 1]. Figure 17(a) displays the transition diagram for the bivariate Markov chain $(X_{1},X_{2})$. The transitions in blue (dotted) correspond to service completions at the first station that result in a customer entering the second station, meaning $N_{1}$ jumps by $1$ for each occurrence of these transitions; whereas the transitions of $(X_{1},X_{2})$ in black (solid) correspond to those that do not result in a jump in $N_{1}$.

In order to obtain an expression for $\kappa_{1}$, we require some additional notation. Let

• •

$\boldsymbol{1}_{k}$ be a $k$-tuple containing all unit entries;

• •

$I_{k}$ be a $k\times k$ identity matrix;

• •

$\overline{I}_{k}$ be a $k\times k$ matrix with upper diagonal containing unit entries and otherwise 0;

• •

$\underline{I}_{k}$ be a $k\times k$ matrix with lower diagonal containing unit entries and otherwise 0;

• •

$\>\otimes\>$ denote the usual Kronecker product for matrices;

• •

$\text{o}(t)\to 0$ denote a function $h(t)$ such that $h(t)/t\to 0$ as $t\to\infty$.

Finally, using [30, Theorem 1] (see also, e.g., [18, Equation (2.110)]), the expected number of customers successfully entering the second station is given by

[TABLE]

where $\boldsymbol{\pi}$ is the stationary distribution of the Markov chain generated by the infinitesimal generator $\mathscr{D}$. Namely, $\boldsymbol{\pi}\,\mathscr{D}=0$ and $\boldsymbol{\pi}\boldsymbol{1}_{(c_{1}+1)(c_{2}+1)}=1$) with

[TABLE]

where

[TABLE]

utilizing the $(c_{1}+1){\times}(c_{1}+1)$ matrices $\mathscr{L}$, $\mathscr{M}_{1}$, and $(\mathscr{Q}_{j},\;j=0,\dots,c_{2})$, with non-zero entries $(\mathscr{L})_{i,i+1}=\lambda_{1}$, $(\mathscr{M}_{1})_{i+1,i}=i\mu_{1}$ and $(\mathscr{Q}_{j})_{i,i}=-(\lambda_{1}+i\,\mu_{1}+j\,\mu_{2})$ for $i\in[c_{1}]$, and utilizing the $(c_{2}+1){\times}(c_{2}+1)$ matrix $\mathscr{M}_{2}$ with non-zero entries $(\mathscr{M}_{2})_{i+1,i}=i\mu_{2}$ for $i\in[c_{2}]$.

Finally, the quantity of interest $\kappa_{1}$ follows immediately from differentiation of this expression and multiplication by the appropriate cost term ${\theta}_{1}$. Based on this we can explicitly compute (4) for a tandem network of loss systems of the Model I variant.

A.2 Model II

Next consider Model II and let \big{(}N_{1}(t),\,t\in\mathbb{R}_{0}\big{)} and \big{(}N_{2}(t),\,t\in\mathbb{R}_{0}\big{)} count the number of customers that successfully enter the first and second stations, respectively. The rates at which $N_{1}$ and $N_{2}$ increase correspond to $\kappa_{1,1}$ and $\kappa_{1,2}$, respectively, which are the quantities we wish to identify. As above, let $X_{1}(t)$ and $X_{2}(t)$ be the number of customers that are processed by the system at the first and second stations, respectively. When an arrival occurs, as long as $X_{1}<c_{1}$, then $X_{1}$ jumps by 1 and, for this variant of the model, so does $N_{1}$. In the case that $X_{1}=c_{1}$ at the time of an arrival, if $X_{2}<c_{2}$, then $X_{2}$ and $N_{2}$ jump up by 1. Alternatively, if $X_{1}=c_{1}$ and $X_{2}=c_{2}$ at the time of an arrival, then $X_{1}$, $X_{2}$, $N_{1}$, and $N_{2}$ remain unchanged—the arrival is simply ignored. It can be seen in this case that the bivariate process $(X_{1},X_{2})$ is a background process for the encompassing Marked MAP that experiences arrivals of different types when customers enter the first or second station (see, e.g., [18, Section 2.5] for more details on Marked MAPs). In order to explicitly find $\kappa_{1,1}$ and $\kappa_{1,2}$, we need to carefully devise the structure of the Markov chain $(X_{1},X_{2})$, use its relationship to $N_{1}$ and $N_{2}$, and then apply [30, Theorem 1]. Figure 17(b) displays the transition diagram for the bivariate Markov chain $(X_{1},X_{2})$. The transitions in blue (dotted) correspond to accepted arrivals at the first station, meaning $N_{1}$ jumps by $1$ for each occurrence of these transitions; the transitions in red (dashed) correspond to accepted arrivals at the second station, meaning $N_{2}$ jumps by $1$ for each occurrence of these transitions; and the transitions of $(X_{1},X_{2})$ in black (solid) correspond to those which do not result in an increase of either $N_{1}$ or $N_{2}$.

Then, using [30, Theorem 1] (see also, e.g., [18, Equation (2.110)]), the expected number of customers successfully entering the first and second stations are given by

[TABLE]

respectively, where $\widetilde{\boldsymbol{\pi}}$ is the stationary distribution of the Markov chain generated by the infinitesimal generator $\widetilde{\mathscr{D}}$. Namely, $\widetilde{\boldsymbol{\pi}}\,\widetilde{\mathscr{D}}=0$ and $\widetilde{\boldsymbol{\pi}}\boldsymbol{1}_{(c_{1}+1)(c_{2}+1)}=1$) with

[TABLE]

where

[TABLE]

Here we utilize the matrices defined previously for Model I with $\Lambda$ being a $(c_{1}+1){\times}(c_{1}+1)$ matrix of zeros except for $(\Lambda)_{c_{1}+1,c_{1}+1}=\lambda_{1}$, and $\mathring{I}_{c_{1}+1}$ a $(c_{2}+1){\times}(c_{2}+1)$ matrix with unit diagonal entries except for $(\mathring{I}_{c_{1}+1})_{c_{2}+1,c_{2}+1}=0$.

Appendix B Stochastic approximation implementation

For $\mathcal{C}\subset\mathbb{R}^{d}$, let $\Pi_{\mathcal{C}}$ be a function $\mathbb{R}^{d}\to\mathcal{C}$ that, for $\boldsymbol{x}\in\mathbb{R}^{d}$, returns $\boldsymbol{c}\in\mathcal{C}$ which minimises $||\boldsymbol{c}-\boldsymbol{x}||$. Let $\mathbf{e}_{l}$ be a vector with all components [math] except for component $l$, which is unitary, and let $\{\delta^{(n)}\}_{n\in{\mathbb{Z}}^{+}}$ and $\{\alpha^{(n)}\}_{n\in{\mathbb{Z}}^{+}}$ be sequences satisfying

[TABLE]

Then, if $\widehat{f}$ are uniformly bounded r.v.s, Algorithm 2 is a convergent (in probability) stochastic approximation algorithm [23].

In our implementation for each iteration of this algorithm, the step size $\alpha^{(n)}$ is determined using the following backtracking line-search method. First choose $\beta>0$, $\varrho_{1}\in(0,1)$, $\varrho_{2}\in(0,1)$, $D>0$. Set $\alpha^{(n)}=\beta^{(n)}$ and $d=1$, and then:

Determine $F=\widehat{f}\left(\Pi_{\mathcal{C}}\left(\boldsymbol{c}^{(n-1)}+\alpha^{(n)}\widehat{\nabla_{\boldsymbol{c}}f}(\boldsymbol{c}^{(n-1)}\right)\right)$; 2. 2.

If $F\geq\widehat{f}(\boldsymbol{c}^{(n-1)})+\varrho_{2}\alpha^{(n)}\widehat{\nabla_{\boldsymbol{c}}f}(\boldsymbol{c}^{(n-1)})\widehat{\nabla_{\boldsymbol{c}}f}(\boldsymbol{c}^{(n-1)})^{\sf T}$ or $d\geq D$: output $\alpha^{(n)}$;

Else: set $\alpha^{(n)}=\varrho_{1}\alpha^{(n)}$ and $d=d+1$, and return to (i).

Choosing $\beta^{(n)}=\beta n^{-1/3}$ and $\delta^{(n)}=\delta n^{-1/6}$ where $\beta,\delta\in\mathbb{R}_{+}$ are positive real numbers ensures that (18) holds and the algorithm therefore asymptotically converges to the true optimiser with probability 1. In all of our experiments we used $\beta=150$, $\delta=5$, $\varrho_{1}=0.8$, $\varrho_{2}=0.5$ and $D=20$.