Optimal Design of Queuing Systems via Compositional Stochastic Programming

TL;DR

This paper introduces a stochastic optimization framework for designing queuing systems with arbitrary, unknown process distributions, using a compositional gradient descent algorithm that handles complex expectation-based objectives and constraints.

Contribution

It develops a novel constrained stochastic compositional gradient descent method for queuing system design, accommodating arbitrary process distributions and providing theoretical performance guarantees.

Findings

The proposed algorithm effectively handles complex queuing design problems.

Theoretical analysis characterizes the iteration complexity of the method.

Numerical tests validate the algorithm's efficacy and theoretical results.

Abstract

Well-designed queuing systems form the backbone of modern communications, distributed computing, and content delivery architectures. Designs balancing infrastructure costs and user experience indices require tools from teletraffic theory and operations research. A standard approach to designing such systems involves formulating optimization problems that strive to maximize the pertinent utility functions while adhering to quality-of-service and other physical constraints. In many cases, formulating such problems necessitates making simplistic assumptions on arrival and departure processes to keep the problem simple. This work puts forth a stochastic optimization framework for designing queuing systems where the exogenous processes may have arbitrary and unknown distributions. We show that many such queuing design problems can generally be formulated as stochastic optimization problems where the objective and constraint are non-linear functions of expectations. The compositional structure obviates the use of classical stochastic approximation approaches where the stochastic gradients are often required to be unbiased. To this end, a constrained stochastic compositional gradient descent algorithm is proposed that utilizes a tracking step for the expected value functions. The non-asymptotic performance of the proposed algorithm is characterized via its iteration complexity. Numerical tests allow us to validate the theoretical results and demonstrate the efficacy of the proposed algorithm.

Tables

Table 1. TABLE I : Summery of bounds on optimality gap and constraint violation

Choice of constants $a,b,c$	Optimality gap	Constraint violation	$\gamma$
$a=0.9167$, $b=0.5$, $c=0.75$	$𝒪\left(T^{-1/12}\right)$	0	$𝒪\left(T^{-1/12}\right)$
$a=0.9167$, $b=0.5$, $c=0.75$	$𝒪\left(T^{-1/12}\right)$	$𝒪\left(T^{-1/12}\right)$	0
$a=0.875$, $b=0.5$, $c=0.75$	$𝒪\left(T^{-1/8}\right)$	$𝒪\left(T^{-1/16}\right)$	0
$a=0.75$, $b=0.5$, $c=0.75$	$𝒪\left(T^{-1/4}\right)$	$𝒪\left(1\right)$	0

Equations

$${\mathbf{x}}^{\star }=\arg \underset{\mathbf{x}\in \mathcal{X}}{\min }f(\mathbb{E}\left[\mathbf{g}(\mathbf{x},\mathbf{\zeta })\right])\text{s.t. }\mathbf{q}(\mathbb{E}\left[\mathbf{h}(\mathbf{x},\mathbf{\zeta })\right])\le 0$$

$$U(\mu )=r\frac{\lambda }{\mu }-h\frac{\lambda /\mu }{\mu -\lambda }$$

$$[\mathbf{g}(\mathbf{\lambda },\mathbf{L}){]}_i=\{\begin{array}{ll}{\lambda }_i{L}_i& 1\le i\le N\\ {\lambda }_{i-N}{L}_{i-N}^2& N+1\le i\le 2N\end{array},f(\mathbf{y})=\sum _{i=1}^N\left[\phi \left(\frac{y_{N+i}}{2C_i(C_i-y_i)}\right)-\psi (y_i)\right]$$

$$\underset{1\le i\le N}{\max }\frac{{\lambda }_i\mathbb{E}{L}_i^2}{2C_i(C_i-{\lambda }_i\mathbb{E}{L}_i)}\le D^{\max }$$

$$\underset{\mathbf{\lambda }\in \mathbf{\Lambda }}{\min }U(\mathbf{\lambda })\text{s.t. }\underset{1\le i\le N}{\max }\frac{{\lambda }_i\mathbb{E}{L}_i^2}{2C_i(C_i-{\lambda }_i\mathbb{E}{L}_i)}\le D^{\max },$$

$$U(\mathbf{\lambda },\mathbf{p})=\sum _{i=1}^N\left[\psi ({\lambda }_i)-{\phi }_i\left(\frac{{\lambda }_i\mathbb{E}\left[\frac{1}{b_i^2(p_i,{\zeta }_i)}\right]}{2\left(1-{\lambda }_i\mathbb{E}\left[\frac{1}{b_i(p_i,{\zeta }_i)}\right]\right)}\right)\right].$$

$$[\mathbf{g}(\mathbf{\lambda },\mathbf{p},\mathbf{\zeta }){]}_i=\{\begin{array}{ll}{\lambda }_i& 1\le i\le N\\ \frac{{\lambda }_{i-N}}{b_{i-N}(p_{i-N},{\zeta }_{i-N})}& N+1\le i\le 2N\\ \frac{{\lambda }_{i-2N}}{b_{i-2N}^2(p_{i-2N},{\zeta }_{i-2N})}& 2N+1\le i\le 3N\end{array},f(\mathbf{y})=\sum _{i=1}^N\left[{\phi }_i\left(\frac{y_{i+2N}}{2(1-y_{i+N})}\right)-\psi (y_i)\right],$$

$$r_i(p_i,{\zeta }_i)=\{\begin{array}{ll}1& R_i<b_i(p_i,{\zeta }_i)\\ 0& R_i\ge b_i(p_i,{\zeta }_i).\end{array}$$

$$[\mathbf{g}(\mathbf{p},\mathbf{h}){]}_i$$

$$f(\mathbf{y})$$

$${\alpha }_i({\theta }_i)=\underset{t\to \infty }{\lim }\frac{1}{t{u}_i}\log \mathbb{E}\left[e^{u_i{A}_i(t)}\right]{|}_{u_i={\theta }_i}=\underset{t\to \infty }{\lim }-\frac{1}{t{u}_i}\log \mathbb{E}\left[e^{-u_i{C}_i(t)}\right]{|}_{u_i={\theta }_i},$$

$$m_i^a+\frac{u_i}{2}{\left({\sigma }_i^a\right)}^2{|}_{u_i={\theta }_i}=m_i^c(p_i)-\frac{u_i}{2}{\left({\sigma }_i^c(p_i)\right)}^2{|}_{u_i={\theta }_i},$$

$${\alpha }_i(p_i)$$

$$[\mathbf{g}(\mathbf{p},\mathbf{\zeta }){]}_i=\{\begin{array}{ll}{b}_i(p_i,{\zeta }_i)& 1\le i\le N\\ b_i^2(p_i,{\zeta }_i)& N+1\le i\le 2N\end{array},$$

$$r_i+l(p_{i+1}-p_i)\le r_{i+1}\le r_i+u(p_{i+1}-p_i).$$

$$a_i(\mathbf{r},C,\text{boldmathζ})$$

$$[\mathbf{g}(\mathbf{r},C,\text{boldmathζ}){]}_i$$

$$\nabla F(\mathbf{x}):=\mathbb{E}\tilde{\nabla }\mathbf{g}(\mathbf{x},\mathbf{\zeta })\nabla f\left(\mathbb{E}\mathbf{g}(\mathbf{x},\mathbf{\zeta })\right),\nabla Q_j(\mathbf{x}):=\mathbb{E}\tilde{\nabla }\mathbf{h}(\mathbf{x},\mathbf{\zeta })\nabla q_j\left(\mathbb{E}\mathbf{h}(\mathbf{x},\mathbf{\zeta })\right)$$

$$\mathbb{E}\left[{\Vert \tilde{\nabla }\mathbf{g}(\mathbf{x},\mathbf{\zeta })\Vert }^2\right]$$

$$\mathbb{E}\left[{\Vert \tilde{\nabla }\mathbf{h}(\mathbf{x},\mathbf{\zeta })\Vert }^2\right]$$

$${\Vert \nabla f(\mathbf{y})\Vert }^2$$

$${\Vert \nabla \mathbf{q}(\mathbf{z})\Vert }^2$$

$${\ell }_j(x):=\{\begin{array}{ll}\frac{1}{2}{x}^2& 0\le x\le C_{\ell }\\ C_{\ell }x-\frac{C_{\ell }^2}{2}& x>C_{\ell }\\ 0& x<0\end{array}$$

$${\mathbf{y}}_{t+1}$$

$${\mathbf{x}}_{t+1}$$

$${\mathbf{x}}^{\gamma }={\mathrm{argmin}}_{\mathbf{x}\in \mathcal{X}}F(\mathbf{x}),\text{s.t. }\mathbf{Q}(\mathbf{x})+\gamma 1\le 0$$

$$\sum _{t=T/2}^T\mathbb{E}\left[F({\mathbf{x}}_t)+\frac{{\delta }_t}{{\alpha }_t}L({\mathbf{x}}_t)-F({\mathbf{x}}^{\gamma })\right]\le \mathcal{O}\left(\frac{1}{{\alpha }_T}+\sum _{t=T/2}^T\left[\frac{{\delta }_t^2}{{\alpha }_t}+\frac{{\beta }_t^2}{{\alpha }_t}+\frac{{\delta }_t^2}{{\alpha }_t{\beta }_t}\right]\right).$$

$$\omega =2\frac{D_1}{T{\alpha }_T}+2\frac{D_2}{T}\sum _{t=T/2}^T\left[\frac{{\delta }_t^2}{{\alpha }_t}+\frac{{\beta }_t^2}{{\alpha }_t}+\frac{{\delta }_t^2}{{\alpha }_t{\beta }_t}\right]$$

$$\mathbb{E}\left[F(\hat{\mathbf{x}})\right]-F^{\star }\le \omega +\frac{2\sqrt{C_g{C}_f}{D}_x}{\sigma }\gamma ,\underset{j}{\max }{Q}_j(\hat{\mathbf{x}})\le \sqrt{\frac{J{T}^{c-a}\left(\omega +\sqrt{C_f{C}_g}{D}_x\right)}{2^{c-a}}}-\gamma$$

$$\sum _{t=T/2}^T\mathbb{E}\left[F({\mathbf{x}}_t)+\frac{{\delta }_t}{{\alpha }_t}L({\mathbf{x}}_t)-F({\mathbf{x}}^{\gamma })\right]\le \omega \frac{T}{2}.$$

$$\frac{2}{T}\sum _{t=T/2}^T\mathbb{E}\left[F({\mathbf{x}}_t)-F({\mathbf{x}}^{\star })\right]\le \omega +\frac{2}{T}\sum _{t=T/2}^T\mathbb{E}\left[F({\mathbf{x}}^{\gamma })-F({\mathbf{x}}^{\star })\right]\le \omega +\frac{2\sqrt{C_g{C}_f}{D}_x}{\sigma }\gamma .$$

$$\frac{2}{T}\sum _{t=T/2}^T\frac{{\delta }_t}{{\alpha }_t}\mathbb{E}\left[L({\mathbf{x}}_t)\right]\le \omega +\sqrt{C_f{C}_g}{D}_x.$$

Full text

This work puts forth a stochastic optimization framework for designing queuing systems where the exogenous processes may have arbitrary and unknown distributions. We show that many such queuing design problems can generally be formulated as stochastic optimization problems where the objective and constraint are non-linear functions of expectations. The compositional structure obviates the use of classical stochastic approximation approaches where the stochastic gradients are often required to be unbiased. To this end, a constrained stochastic compositional gradient descent algorithm is proposed that utilizes a tracking step for the expected value functions. The non-asymptotic performance of the proposed algorithm is characterized via its iteration complexity. Numerical tests allow us to validate the theoretical results and demonstrate the efficacy of the proposed algorithm.

I Introduction

As the Internet traffic continues to grow, meeting service-level agreement (SLA) targets necessitates proactive capacity planning and resource provisioning [1, 2, 3]. Sustaining the users’ quality of experience (QoE), for instance, requires the backbone and access networks to be designed with sufficiently high available capacity. Virtual machine resources, such as CPU, memory, and storage, must likewise be provisioned to tackle the demand variability in data centers. On the other hand, economic considerations generally dictate that resources be utilized efficiently and over-provisioning be minimized. Well-designed systems must strive to achieve a balance between resource utilization and excess capacity.

Teletraffic theory, utilizing advances in queuing theory and operations research, has been widely used for the design of communication networks [4]. Within this framework, the random exogenous factors such as resource demands and channel gains are studied in isolation and statistically modeled a priori (see Chapter 3 of [5]). Simple models are generally preferred, so as to yield closed-form expressions for various performance metrics such as waiting times, throughput, delay, waiting costs, provisioning costs, and quality-of-service (QoS). Subsequently, the underlying decision variables, such as the arrival rates or link capacities, are tuned to maximize the performance or minimize the costs, while adhering to the various QoS constraints. The latter step may be accomplished by solving a constrained optimization problem; see [4] and references therein.

The classical approach is however not suited to systems where the distribution of the exogenous variables is either unknown or too complicated to allow closed-form expressions. As a simple example, the average ergodic capacity of a single antenna wireless fading channel is not expressible in closed-form [6]. Consequently, existing approaches for wireless system design either ignore the delay or rely on approximations even for simplistic fading models [7, 8, 9, 10, 11]. More generally, the performance metrics of almost all such queuing systems can be expressed in terms of expectations involving the random variables. This work considers the design of queuing systems from a stochastic optimization perspective. The network design problem is formulated as a constrained optimization problem where the objective and constraint functions involve expectations with respect to these random variables. We consider handling such problems within the stochastic approximation rubric where a few sequentially available independent realizations of the random process are sufficient to yield near-optimal designs. Of particular interest are low-complexity first-order algorithms that are capable of running in real-time. A supervised learning approach is adopted where the parameters are learned during a training phase that is separate from the significantly longer operational phase. The proposed framework is general, and finds applications in communication networks, traffic systems, data centers, and supply chain.

From a theoretical vantage point, the objective function and constraints in these design problems usually depend on quantities such as queuing delay and throughput, that are themselves non-linear functions of expectations. The compositional structure of these problems obviates the use of standard stochastic gradient descent (SGD) methods where the gradient estimates are required to be unbiased or at least strongly consistent (see e.g., [12, 13, 14]). The more general case involving non-linear functions of expectations has only been recently considered in [15], where the stochastic compositional gradient descent (SCGD) algorithm was proposed. However, compositional problems involving expectation constraints have never been considered in the literature, and the proposed SCGD method for constrained problems is the first such algorithm. The non-asymptotic optimality gap and constraint violation incurred by the proposed algorithm is characterized in terms of the number of samples, and constitutes the key theoretical contribution of the present work. The sample error complexity results developed here are of independent interest and can be applied to any setting involving constrained stochastic optimization problems.

The rest of the paper is organized as follows. The next subsection briefly reviews the related work on optimization of queuing systems and on algorithms for solving constrained stochastic problems. Sec. II formulates the general constrained problem and discusses several examples from queuing literature. Sec. III details the proposed algorithm and provides the relevant theoretical guarantees. The simulation results are presented in Sec. IV and finally, Sec. V concludes the paper.

I-A Related Work

General compositional stochastic problems, where the goal is to minimize compositions of expected-value functions, were first discussed in [15] and solved via the SCGD algorithm. The corresponding finite-sum variant of the problem have subsequently been considered in [16] and solved via the variance-reduced SCGD. Other variants include the accelerated SCGD [17], those utilizing Markov samples [18], and generalized versions involving multiple nested compositions [19]. A functional variant of SCGD has recently been proposed in [20]. Compositional stochastic problems involving expectation constraints were mentioned in [21, Sec. 6.10] but have not been thoroughly studied.

Within the classical stochastic optimization framework, where the compositional form does not arise, problems with stochastic constraints have been considered. Relevant algorithms include the stochastic dual-descent algorithm [22, 23, 24, 25] and the stochastic variant of the Arrow–Hurwicz saddle point method [26]. The other problem formulations with infinite time horizon objective and constraint functions are discussed in [27]. A comprehensive article on constrained stochastic problems of large data can be found in [28]. It is remarked that all these classes of algorithms apply to problems with a specific structure, and sample complexity results for the general versions have not been developed. Further, these algorithms cannot be directly applied to problems involving compositional forms. Interestingly, the compositional stochastic optimization problem considered here subsumes the setting with only expectation constraints.

Optimization approaches for the design of queuing systems have a long history; see [4] and references therein. Formulated as deterministic optimization problems, these network designs consider a wide variety of objectives and systems. As stated earlier, the existing optimization-based approaches do not suffice when the problem involves random variables with unknown distributions. Optimization of queuing systems where the virtual network functions to be executed in an order are in queue are discussed in [29]. The objective function is formulated as a time average processing cost.

The system design problems considered here are static in nature, and the optimal decision variable obtained at the end of the training phase is utilized during the operational phase. Such a setting is in contrast to that of network control, where the goal is to learn optimal policies that depend on the observed environmental state or queue lengths. Utility-optimal control algorithms such as the backpressure algorithm have been extensively studied [30, 31, 32, 33, 34, 35, 36]. In network control, the optimal policy is updated at every time instant, and the algorithm continues to run forever. However, the focus here is on determining the optimal provisioning or dimensioning parameters that must be set a priori. Accordingly, the performance of the proposed approach is studied in terms of its sample complexity bounds.

II Problem Formulation and Examples

This section formulates the constrained stochastic optimization problem at hand and provides motivating examples within the context of queuing theory. Of particular interest are problems arising in queuing system design, where the distribution of the service times are general and a priori unknown.

II-A Problem Formulation

We begin with describing a general stochastic optimization problem where the formulation involves expectations with respect to a collection of $k$ random variables (collected into the vector ${\bm{\zeta}}\in{\mathbb{R}}^{k}$) with unknown distributions. Going beyond the classical stochastic optimization setting where the goal is to minimize the expectation of a function, we consider the more general case where the expectation appears as an argument of an outer non-linear function. In general, such a composite structure makes it difficult to obtain unbiased samples of the objective gradient, rendering the SGD and its variants inapplicable.

Consider the constrained optimization problem

[TABLE]

where the expectation is taken with respect to $\mathbf{\zeta}$. Here, $f:{\mathbb{R}}^{m}\rightarrow{\mathbb{R}}$, ${\mathbf{g}}:{\mathbb{R}}^{n}\times{\mathbb{R}}^{k}\rightarrow{\mathbb{R}}^{m}$, ${\mathbf{h}}:{\mathbb{R}}^{n}\times{\mathbb{R}}^{k}\rightarrow{\mathbb{R}}^{d}$, and ${\mathbf{q}}:{\mathbb{R}}^{d}\rightarrow{\mathbb{R}}^{J}$ are continuous functions. The optimization variable ${\mathbf{x}}$ belongs to the convex set ${\mathcal{X}}\subset{\mathbb{R}}^{n}$, described by simple constraints that are easy to project onto; examples include box or norm-ball constraints. Specific assumptions regarding convexity or smoothness of various functions will be explicitly stated in Sec III. In order to avoid trivial scenarios however, we assume that the problem in ($\mathcal{P}$) is feasible and has a finite solution.

When the distribution of ${\bm{\zeta}}$ is not known a priori, the expectations appearing in ($\mathcal{P}$) cannot be evaluated in closed-form. Motivated by classical stochastic approximation methods, the goal is to solve ($\mathcal{P}$) in an online fashion using only independent realizations ${\bm{\zeta}}_{1}$, ${\bm{\zeta}}_{2}$, $\ldots$ that are revealed sequentially. For many problems, such samples yield unbiased (or at least strongly consistent) estimates of the objective gradient, paving the way for the application of the popular SGD algorithm and its variants [37, 38]. In the present case however, unbiased estimates of the gradient ${\mathbb{E}}[{\bm{\nabla}}{\mathbf{g}}({\mathbf{x}},{\bm{\zeta}})]\nabla f({\mathbb{E}}{\mathbf{g}}({\mathbf{x}},{\bm{\zeta}}))$ cannot be obtained from a single sample ${\bm{\zeta}}$ due to the presence of the non-linear outer function $f$. Moreover, projected SGD and its variants cannot be used since the constraint also involves an expectation that cannot be evaluated in closed form.

The set-constrained version of ($\mathcal{P}$) was first considered in [15], where a quasi-gradient approach was proposed. The main idea was to run two iterations: one for performing quasi-gradient steps for estimating ${\mathbf{x}}^{\star}$ and another for keeping track of the vector {\mathbb{E}}\left[{{{\mathbf{g}}({\mathbf{x}}^{\star},\text{\boldmath\mathbf{\zeta}})}}\right]. The approach in [15] is however still not applicable to ($\mathcal{P}$) due to the presence of stochastic constraints.

Stochastic constraints can be dealt with algorithms operating in the dual domain. For instance, dual descent algorithms proposed in [22, 23] can handle expectation constraints of a specific form. Likewise, primal-dual or saddle point algorithms have also been proposed for problems with expectation constraints [39, 40, 41]. However, none of these approaches can handle non-linear functions of expectations either in the objective or constraint functions. In summary, existing algorithms do not cater to problems of the form in ($\mathcal{P}$). Nevertheless, such problems are commonplace in queuing systems where the delay is generally a non-linear function of expectation with unknown distributions.

II-B *Design of queuing systems *

In this section, we describe various queuing design examples that adhere to the formulation in ($\mathcal{P}$). As a simple example, consider a first-in-first-out (FIFO) M/M/1 queue with a fixed arrival rate $\lambda$ and service rate $\mu$, which is a design variable. As in many packet networks, the system should be designed to balance the utilization or throughput against the delay. Towards this end, consider an objective function that takes the form

[TABLE]

where $r$ and $h$ are positive scaling factors. The first term in (1) is proportional to the system throughput while the second term depends on the average delay experienced by a packet while waiting to be served. Maximizing $U(\mu)$ with respect to $\mu$ thus allows the designer to trade off throughput with queuing delay. In the simplest case, it can be seen that the optimal solution is given by $\mu^{\star}=\lambda+u+\sqrt{u(u+\lambda)}$ where $u=h/r$. In practice, constraints may be imposed so as to restrict $\mu$ to lie within a physically plausible region and may take the form of upper bounds on $\mu$, delay, queue length, and delay variance. Such an optimization framework is general and is applicable to a large class of queuing systems that admit closed-form expressions for utility and delay. See [4] and [42] for more examples where deterministic optimization algorithms have been applied to optimize the performance of a queuing systems.

In this section, several queuing system design examples are discussed and the resulting stochastic problems are shown to be of the form in ($\mathcal{P}$). For ease of exposition, we consider the standard utility-minus-penalty objective similar to (1) subject to various functional constraints. Other metrics can also be used, provided that they satisfy the assumptions in Sec. III. Alternatively, all requirements may be expressed as constraints only, so that ($\mathcal{P}$) becomes a feasibility problem.

II-B1 Example 1: Parallel M/G/1 queues

To begin with, consider a set of $N$ parallel $M/G/1$ queues as shown in Fig. 1(a). The system consists of a generic source, mimicking a router or an encoder, and capable of supplying at most $\lambda^{\max}$ packets per second on an average. Each packet generated by the source carries a destination address and is accordingly routed to one of the queues. The length of a packet sent to the $i$-th queue is measured in bits and denoted by the random variable $L_{i}$ and collected into the random vector ${\mathbf{L}}$. The lengths are exponentially distributed, but the means ${\mathbb{E}}L_{i}$ are unknown and different for each queue. The difference in mean packet lengths may arise due to the heterogeneity in the traffic stream arriving at each queue. The departures from the $i$-th queue are limited by the link capacity $C_{i}$ bits per second and the packet service times $X_{i}:=L_{i}/C_{i}$.

The optimization variables $\{\lambda_{i}\}$ are collected into the vector ${\bm{\lambda}}\in{\mathbb{R}}^{N}_{++}$ and are required to satisfy box constraints $\lambda_{i}\in[\lambda^{\min},\lambda^{\max}]$ that may arise from various physical considerations. As in (1), the goal is to maximize the utility function that depends on the throughput, and the delay. While throughput is still given by $\lambda_{i}{\mathbb{E}}L_{i}$, the queuing delay must be calculated using the Pollaczek-Khinchin (P-K) formula, and the objective function becomes $U(\bm{\lambda})=\sum_{i=1}^{N}\left[\psi_{i}(\lambda_{i}{\mathbb{E}}L_{i})-\varphi\left(\lambda_{i}{\mathbb{E}}L^{2}_{i}/\left({2C_{i}(C_{i}-\lambda_{i}{\mathbb{E}}L_{i})}\right)\right)\right]$, where $\psi_{i}:{\mathbb{R}}\rightarrow{\mathbb{R}}$ and $\varphi_{i}:{\mathbb{R}}\rightarrow{\mathbb{R}}$ are pre-specified functions. Examples include linear utility and cost functions, fairness-inducing log utility $\psi(x)=\bar{\psi}_{i}\log(x)$, and squared delay penalty $\varphi_{i}(x)=\bar{\varphi}_{i}x^{2}$, where $\bar{\psi}_{i}$ and $\bar{\varphi}_{i}$ are positive constants. It can be seen that the objective function depends on a non-linear function of the expectation and can be written in the form required in ($\mathcal{P}$) as

[TABLE]

where $m=2N$ and $n=k=N$. The waiting delay of a packet should not exceed the maximum tolerable value, and hence the following constraint is imposed:

[TABLE]

where $D^{\max}$ is the maximum tolerable delay. Comparing with the formulation in ($\mathcal{P}$), it can be seen that ${\mathbf{h}}({\bm{\lambda}},{\mathbf{L}})={\mathbf{g}}({\bm{\lambda}},{\mathbf{L}})$ and $q({\mathbf{z}})=\max_{1\leq i\leq N}z_{i+N}/\left({2C_{i}(C_{i}-z_{i})}\right)-D^{\max}$ so that $d=2N$. The overall optimization problem can be written as

[TABLE]

where ${\bm{\Lambda}}:=\{{\bm{\lambda}}|\lambda^{\min}\leq\lambda_{i}\leq\lambda^{\max}\forall i=1,\ldots,N,\sum_{i=1}^{N}\lambda_{i}\leq\lambda^{lim}\}$.

II-B2 Example 2: Parallel M/G/1 queues in wireless networks (ergodic capacity)

Consider a system of $N$ parallel M/G/1 queues operating within a common base station (BS) and communicating with $N$ users over $N$ orthogonal channels; see Fig. 1(b). The arrivals are still Poisson and the arrival rate vector ${\bm{\lambda}}$ is an optimization variable. Different from the earlier examples however, the packet lengths are constant and the service time instead depends on the channel capacity and the power allocation. Each queue is allocated with an average power of $p_{i}$ and the queues adhere to the total power budget $P^{\max}$, that is, $\sum_{i=1}^{N}p_{i}\leq P^{\max}.$ The allocated powers are collected into the vector ${\mathbf{p}}\in{\mathbb{R}}^{N}_{+}$. Together, the arrival rates and power allocations belong to the sets ${\bm{\Lambda}}:=\{{\bm{\lambda}}|\lambda^{\min}\leq\lambda_{i}\leq\lambda^{\max},i=1,\ldots,N,\sum_{i}\lambda_{i}\leq\lambda^{\max}\}$ and ${\bm{\Pi}}:=\{{\mathbf{p}}\mid p_{i}\geq P^{\min},\sum_{i}p_{i}\leq P^{\max}\}$.

The channel gain from the BS to the $i$-th user is denoted by $\zeta_{i}$, collected in the vector ${\bm{\zeta}}\in{\mathbb{R}}^{N}$, and remains constant over the transmission of a single packet. Assuming sufficiently long coherence times, the transmission rate is close to the ergodic capacity of the channel, given by $b_{i}(p_{i},\zeta_{i})=B_{i}\log\left(1+\zeta_{i}p_{i}\right)$ $\forall i$, where $B_{i}$ is the capacity of the $i$-th channel and is known a priori. If the capacity is measured in packets per second, the service time is the random variable $X_{i}=1/b_{i}(p_{i},\zeta_{i})$ measured in seconds. The design goal is to determine the optimal values of ${\bm{\lambda}}$ and ${\mathbf{p}}$ so as to maximize the utility function. Applying the PK formula for queuing delay, the utility function can be written as

[TABLE]

It can be seen that the objective function in (4) is of the compositional form in ($\mathcal{P}$) with

[TABLE]

and $m=3N$, $n=2N$, and $k=N$. In order to ensure quality of service to every user, it is also required that the physical layer rate for any user does not drop below a specified level, i.e., ${\mathbb{E}}[\min_{i}b_{i}(p_{i},\zeta_{i})]>R^{\min}$. In other words, on average, the minimum physical layer rate is not below $R^{\min}$. The constraint is again in the form required in ($\mathcal{P}$), with $q(z)=z$, and $h({\bm{\lambda}},{\mathbf{p}},{\bm{\zeta}})=-\min_{i}b_{i}(p_{i},\zeta_{i})$ so that $d=1$. Unlike in Example 1 however, it can be seen that empirical estimates of the moments of ${\bm{\zeta}}$ no longer suffice. Instead, the proposed algorithm will entail tracking ${\mathbb{E}}\left[{{{\mathbf{g}}({\bm{\lambda}}^{\star},{\mathbf{p}}^{\star},{\bm{\zeta}})}}\right]$ and ${\mathbb{E}}\left[{{{\mathbf{h}}({\bm{\lambda}}^{\star},{\mathbf{p}}^{\star},{\bm{\zeta}})}}\right]$ explicitly.

II-B3 Example 3: Parallel M/G/1 queues in wireless networks (outage capacity)

Ergodic capacity is generally difficult to achieve in practical systems with smaller coherence intervals. We consider the joint design of rate and power for systems where both outage and retransmission events are common. This example builds upon the queuing system described in [43].

Within the outage-limited setting, a packet in the $i$-th queue is transmitted at a fixed rate $R_{i}$ regardless of the instantaneous channel conditions. The transmission is deemed successful if $R_{i}<b_{i}(p_{i},\zeta_{i})$ and an outage is declared otherwise. In the case of outage, the server continues to attempt to retransmit a packet until it is successful. Define the indicator function for a successful transmission as

[TABLE]

Then the service times depend on the number of retransmission attempts, which in turn depends on the probability of successful transmission $\rho_{i}={\mathbb{E}}[r_{i}(p_{i},\zeta_{i})]$. In particular, it can be shown that for this case, the mean waiting time for a packet is given by $w_{i}(\lambda_{i},p_{i})=\lambda_{i}(1+\rho_{i})/\left({2R_{i}(1-\rho_{i})(R_{i}(1-\rho_{i})-\lambda_{i})}\right)$ while the throughput is given by $R_{i}(1-\rho_{i})$. The overall utility function can be written as $U({\bm{\lambda}},{\mathbf{p}})=\sum_{i=1}^{N}\left[\psi(R_{i}(1-\rho_{i}))-\varphi_{i}\left(w_{i}(\lambda_{i},p_{i})\right)\right]$ which is also of the form required in ($\mathcal{P}$) with

[TABLE]

Within the context of M/G/1 queues, we briefly discuss some more examples that serve to illustrate the flexibility and generality of the problem considered here.

Delay variance captures the jitter and may be an important QoS parameter in streaming and gaming applications. The delay variance can be minimized by adding terms of the form $({\mathbb{E}}[w_{i}(\lambda_{i},p_{i})])^{2}-{\mathbb{E}}[w_{i}(\lambda_{i},p_{i})^{2}]$ in Example 3. Likewise, various moments of delay or throughput may also be incorporated. 2. 2.

More generally, tail probabilities of the delay or throughput may be considered. For instance in Example 3, where the service time $X_{i}$ is a random variable, the tail probability can be written as $-\lim_{T\rightarrow\infty}\left({1/T}\right)\mathbb{P}\left(\left({1/T}\right)\sum_{t=1}^{T}X_{i,t}>x\right)=\max_{\theta}\left[\theta x-\log{\mathbb{E}}[\exp(\theta X_{i})]\right]$, and may therefore be incoporated within the objective function. 3. 3.

Worst-case delay can be considered by using the term $\max_{1\leq i\leq N}w_{i}(\lambda_{i},p_{i})$ in the utility function. Likewise, other non-linear combinations of the delays and throughput may be used. 4. 4.

Expectation constraints can also be incorporated via penalty methods. For instance, instead of imposing the constraint ${\mathbb{E}}[b_{i}(p_{i},\zeta_{i})]\geq\epsilon$ in Example II-B2, it is possible to use a log-barrier penalty of the form $\log\left({\mathbb{E}}[b_{i}(p_{i},\zeta_{i})]-\epsilon\right)$ can be used in the objective function. Likewise, instead of imposing the constraint of maximum delay $w_{i}(\lambda_{i},p_{i})\leq w^{\max}$, a log-barrier penalty of the form $\log\left(w^{\max}-w_{i}(\lambda_{i},p_{i})\right)$ can be used. Log penalties have been widely used in the context of interior point methods. 5. 5.

Probabilistic constraints can be incorporated much in the same way as in Example II-B3. For instance, a constraint of the form $\mathbb{P}(b_{i}(p_{i},\zeta_{i})>R_{i})>0.9$ can be converted into an expectation constraint ${\mathbb{E}}[r_{i}(p_{i},\zeta_{i})]>0.9$ and dealt with using an appropriate penalty function.

II-B4 Example 4: Parallel G/G/1 queues in wireless networks (effective capacity)

This example considers the more complicated G/G/1 queue where the expressions for delay cannot be characterized in closed form. The effective capacity, defined as the minimum rate required to guarantee a certain delay (say $W$), is an approximation that allows handling queues with arbitrary arrival and service distributions. Effective capacity has been widely studied for fading channels in systems where a maximum delay of $W$ units may be tolerated [7, 9, 10].

Let us consider a similar model as in Fig. 1(b) but with $N$ G/G/1 queues. The packet lengths are constant and the departures depend on the random channel gains and allocated powers $\{p_{i}\in{\mathbb{R}}_{+}\}_{i=1}^{N}$. As in Examples 2, 3, the power vector ${\mathbf{p}}$ is required to lie in the compact set ${\bm{\Pi}}$. Different from the earlier examples and for the sake of simplicity, it is assumed that the traffic arrival process at the $i$-th queue has mean $m_{i}^{a}$ and variance $(\sigma_{i}^{a})^{2}$, both of which are known a priori.

Before formulating the problem, let us briefly review the theory of effective capacity. At any time slot $t$, let $a_{i}(t)$ denote the instantaneous arrival rate at the $i$-th queue and $c_{i}(t)$ denote the instantaneous channel capacity of the $i$-th channel. The corresponding accumulated rates are denoted by $A_{i}(t):=\sum_{\tau=1}^{t}a_{i}(t)$ and $C_{i}(t)=\sum_{\tau=1}^{t}c_{i}(t)$, respectively. The effective capacity of the $i$-th user is given by [9, Eq. (7)]

[TABLE]

where $\theta_{i}$ is the QoS exponent and is the solution to the above equation [10, Eq. (3)]. The steady-state approximation is often used to simplify these expressions. Specifically, as $t\rightarrow\infty$, the application of the central limit theorem allows us to assume the averages $A_{i}(t)/t$ and $C_{i}(t)/t$ as Gaussian random variables with finite means and variances.

In the present case, recall that the instantaneous channel capacity is given by $b_{i}(p_{i},\zeta_{i})$ where $\zeta_{i}$ is the instantaneous channel. Therefore the mean and variance of $C_{i}(t)/t$ can be approximated as $m^{c}_{i}(p_{i})={\mathbb{E}}\left[{{b_{i}(p_{i},\zeta_{i})}}\right]$, and $\sigma^{c}_{i}(p_{i})=\sqrt{{\mathbb{E}}\left[{{b_{i}^{2}(p_{i},\zeta_{i})}}\right]-\left({{\mathbb{E}}\left[{{b_{i}(p_{i},\zeta_{i})}}\right]}\right)^{2}}$ respectively. As stated earlier, we assume that the means $\{m_{i}^{a}\}$ and variances $\{(\sigma_{i}^{a})^{2}\}$ of the arrival processes are known. Hence (9) can be rewritten as [10, Eq. (6)]

[TABLE]

implying that the QoS quotient $\theta_{i}$ is given by $\theta_{i}(p_{i})=\left({m^{c}_{i}(p_{i})-m^{a}_{i}}\right)/\left({\left({\sigma^{a}_{i}}\right)^{2}+\left({\sigma^{c}_{i}(p_{i})}\right)^{2}}\right).$ where we have explicated the dependence on the optimization variables $\{p_{i}\}$. Substituting into (9), we obtain the effective capacity as

[TABLE]

Proceeding along similar lines, the approximation for the probability of delay exceeding the pre-specified value of $W$ time slots is given by ${\mathbb{P}}\left\{{W_{i}(t)|_{t\rightarrow\infty}>W}\right\}\approx\eta e^{-\theta_{i}(p_{i})\alpha_{i}(p_{i})W}$ where $\eta$ is a normalization constant. So the utility function for this case can be defined as $U({\mathbf{p}})=\sum_{i=1}^{N}\left[\psi_{i}\left(\alpha_{i}(p_{i})\right)-\varphi_{i}\left(e^{-\theta_{i}(p_{i})\alpha_{i}(p_{i})W}\right)\right]$. The utility is therefore in the compositional form with the inner function

[TABLE]

and outer function is given by $f({\mathbf{y}})=\sum_{i=1}^{N}\left[{\psi_{i}\left(\alpha_{i}(y_{i},y_{i+N})\right)+\varphi_{i}\left(e^{-\theta_{i}(y_{i},y_{i+N})\alpha_{i}(y_{i},y_{i+N})W}\right)}\right]$, where $\theta_{i}(u,v)=\frac{u-m_{i}^{a}}{(\sigma_{i}^{a})^{2}+v-u^{2}}$ and $\alpha_{i}(u,v)=m_{i}^{a}+\frac{\theta_{i}(u,v)}{2}(\sigma_{i}^{a})^{2}$. In conclusion, it can be seen that the proposed formulation in ($\mathcal{P}$) is applicable to generalized capacity definitions that also involve non-linear functions of expectations.

II-B5 Example 5: Resource provisioning in cloud services

Cloud computing service companies seek to ensure that resources such as power, CPU, and storage are made available to the users as per their demands. Given the uncertainties in the user demands however, the task of distributing the resources to the users is challenging. Ideally, the resources must be provisioned so as to minimize the probability of SLA violation. This example formulates the resource provisioning problem from [2, 3] for the case when demand distributions are unknown.

Consider a service provider that offers $N$ different classes of service and offers a total of $C$ resource units. A customer subscribes to the $i$-th class by paying a service charge of $p_{i}$ units in return for $r_{i}$ units of resources. The service provider also incurs a maintenance charge of $\chi$ units per unit time, and per unit resource. An incoming customer may attempt to subscribe to class $i$ of its choice, but may be blocked if the available resources are less than $r_{i}$ units. The goal is to determine the optimal number of resource units $r_{i}$ corresponding to each resource class and the total number of resource units to be provisioned.

Let $\zeta_{i}$ denote the load for class $i$, which represents the number of customers willing to subscribe to class $i$ if the service provider was not resource-constrained. Then the blocking probability for the $i$-th class is given by [3, Eq. 1] $P^{B}_{i}({\mathbf{r}},C)={\mathbb{P}}(C-r_{i}<\sum_{j=1}^{N}\zeta_{j}r_{j}\leq C)/{\mathbb{P}}(\sum_{j=1}^{N}\zeta_{j}r_{j}\leq C).$ Roughly, the net revenue obtained by the service provider is given by $R=\sum_{i}p_{i}(1-P^{B}_{i}({\mathbf{r}},C))$. On the other hand, the maintenance costs incurred by the service provider amount to $\chi C$ units per unit time. Consequently, the provider’s average profit is given by $U({\mathbf{r}},C)=\sum_{i=1}^{N}p_{i}N_{i}(1-P^{B}_{i}({\mathbf{r}},C))-\chi C$, where $N_{i}$ is the average number of customers actually subscribing to class $i$ per unit time. Finally, the revenue offered for the $i$-th class should be a function of the price $p_{i}$ for that class. Without loss of generality, let $p_{1}<p_{2}<\ldots<p_{N}$. Then, the resource units offered for class $i+1$ should be sufficiently different from that offered for class $i$, and may be restricted to

[TABLE]

where $l$ and $u$ are constants. The goal is to maximize the utility $U({\mathbf{r}},C)$ subject to constraints in (11). The problem can be written in the compositional form by defining

[TABLE]

so that the inner function is given by

[TABLE]

and the outer function is $f({\mathbf{y}})=\sum_{i=1}^{N}p_{i}N_{i}\left({1-\frac{y_{i}}{y_{N+i}}}\right)-\chi y_{2N+1}$. This example demonstrates the applicability of the compositional form to problems involving probabilistic functions such as $P^{B}_{i}({\mathbf{r}},C)$.

Before concluding, it is remarked that the objective function and constraints in the examples are specifically selected for simplifying the exposition. More general utility functions and constraints are also possible, as long as they adhere to the compositional structure of ($\mathcal{P}$).

III Stochastic Compositional Gradient Descent for Constrained Problems

This section details the proposed algorithm for solving ($\mathcal{P}$) and provides the corresponding convergence rates. For the sake of brevity, we define $F({\mathbf{x}}):=f({\mathbb{E}}\left[{{{\mathbf{g}}({\mathbf{x}},{\bm{\zeta}})}}\right])$ and ${\mathbf{Q}}({\mathbf{x}}):={\mathbf{q}}({\mathbb{E}}\left[{{{\mathbf{h}}({\mathbf{x}},{\bm{\zeta}})}}\right])$. In the unconstrained case, ($\mathcal{P}$) can be solved using the SCGD algorithm proposed in [15] for convex non-smooth or smooth non-convex $F$. This section focuses on the constrained case when both $F$ and $Q$ are convex but possibly non-smooth.

III-A Assumptions

We begin with discussing the necessary assumptions on the functions $f$, ${\mathbf{g}}$, ${\mathbf{q}}$ and ${\mathbf{h}}$.

A1*.*

The outer functions $f$, ${\mathbf{q}}$ are continuously differentiable and the inner functions ${\mathbf{g}}$, ${\mathbf{h}}$ are continuous. Consequently, the sub-gradients of the objective and constraint functions are well-defined, with

[TABLE]

where $[{\widetilde{\nabla}}{\mathbf{g}}({\mathbf{x}},{\bm{\zeta}})]_{ij}:=\partial g_{j}({\mathbf{x}},{\bm{\zeta}})/\partial x_{i}$ and $[{\widetilde{\nabla}}{\mathbf{h}}({\mathbf{x}},{\bm{\zeta}})]_{ij}:=\partial h_{j}({\mathbf{x}},{\bm{\zeta}})/\partial x_{i}$ with $\partial$ denoting the generalized gradient operation.

A2*.*

The problem ($\mathcal{P}$) is a convex optimization problem and the constraints can be satisfied strictly, i.e., there exists $\widetilde{{\mathbf{x}}}$ such that $\max_{j}Q_{j}(\widetilde{{\mathbf{x}}})+\sigma\leq 0$ for some constant $\sigma>0$. The set ${\mathcal{X}}$ is closed and compact, i.e., $\sup_{{\mathbf{x}},{\mathbf{x}}^{\prime}\in{\mathcal{X}}}\left\|{\mathbf{x}}-{\mathbf{x}}^{\prime}\right\|^{2}\leq D_{x}<\infty$.

A3*.*

The random variables ${\bm{\zeta}}_{1},{\bm{\zeta}}_{2},...$ are independent and identically distributed.

A4*.*

The functions ${\mathbf{g}}$, ${\mathbf{h}}$ are Lipschitz continuous in expectation and have bounded second order moments, i.e.,

[TABLE]

A5*.*

The functions $f$, $q$ are smooth and have bounded gradients, i.e.,

[TABLE]

These assumptions are standard and hold for most problems at hand. In particular, we do not require ${\mathbf{g}}$ and ${\mathbf{h}}$ to be smooth or convex. Only the outer functions $f$, $q$ are required to be smooth and only the composite functions $F$ and ${\mathbf{Q}}$ are required to be convex. As in [15], the matrices ${\widetilde{\nabla}}{\mathbf{g}}({\mathbf{x}},{\bm{\zeta}})$ and ${\widetilde{\nabla}}{\mathbf{h}}({\mathbf{x}},{\bm{\zeta}})$ may consist of gradients, sub-gradients, or other directional gradients satisfying (14). Further, none of the functions are required to be strongly convex. Note that the Lipschitz continous function ${\mathbf{Q}}$ is defined over a compact domain, and hence is bounded.

III-B Proposed Algorithm

In order to describe the proposed algorithm, let us also define the smooth and convex function $\ell({\mathbf{w}})=\sum_{j=1}^{J}\ell_{j}(w_{j}+\gamma)$ where $\ell_{j}$ is defined as

[TABLE]

and $\gamma$ and $C_{\ell}$ are algorithm parameters. It can be seen that ${\nabla}\ell_{j}(x)=\max\{x,C_{\ell}\}$ for $x>0$ and zero otherwise. Such a function allows us to define a convex penalty function $L({\mathbf{x}}):=\ell({\mathbf{q}}({\mathbb{E}}\left[{{{\mathbf{h}}({\mathbf{x}},{\bm{\zeta}})}}\right]))$ that quantifies the constraint violation and allows us to control it through the choice of the parameter $\gamma$. The penalty function defined in (19) is quadratic and imposes square penalty on the constraint violation. Intuitively, the update takes the iterate towards the descent direction of optimal function $F$ as well as towards the feasible region $\{{\mathbf{x}}|{\mathbf{Q}}({\mathbf{x}})\leq 0\}$.

Recall that the SCGD algorithm for the unconstrained version of $\mathcal{P}$ entails carrying out updates of the form ${\mathbf{x}}_{t+1}={\mathbf{x}}_{t}-\alpha_{t}\hat{{\nabla}}F({\mathbf{x}}_{t})$ where $\hat{{\nabla}}F({\mathbf{x}}_{t})$ is an approximation of the actual gradient ${\nabla}F({\mathbf{x}}_{t})$ and $\alpha_{t}$ is the step size. To handle the constrained case, we consider the iterations ${\mathbf{x}}_{t+1}={\mathbf{x}}_{t}-\alpha_{t}\hat{{\nabla}}F({\mathbf{x}}_{t})-\delta_{t}\hat{{\nabla}}L({\mathbf{x}}_{t})$ where $\delta_{t}$ is also a step size parameter and $\hat{{\nabla}}L({\mathbf{x}}_{t})$ approximates the true gradient ${\nabla}L({\mathbf{x}}_{t})$. Finally, the approximate gradient $\hat{{\nabla}}L({\mathbf{x}}_{t})$, and $\hat{{\nabla}}F({\mathbf{x}}_{t})$ are calculated by introducing new variables ${\mathbf{z}}_{t}$, and ${\mathbf{y}}_{t}$ that track the quantities ${\mathbb{E}}\left[{{{\mathbf{h}}({\mathbf{x}}^{\star},{\bm{\zeta}})}}\right]$, and ${\mathbb{E}}\left[{{{\mathbf{g}}({\mathbf{x}}^{\star},{\bm{\zeta}})}}\right]$ are respectively. The step sizes $\alpha_{t}$, and $\delta_{t}$ are chosen such that the updates converge to a point that is not only optimal but also feasible. The full updates are summarized in Algorithm 1. The operation $\Pi_{{\mathcal{X}}}$ denotes projection onto the set ${\mathcal{X}}$. The computational complexity of the proposed algorithm is only slightly higher than that of the classical stochastic gradient descent method. At each time instant, only gradients of ${\mathbf{g}}$ and ${\mathbf{h}}$ are needed in order to carry out the updates. As compared to the SCGD algorithm for unconstrained setting [15], the proposed algorithm incurs only a single additional update. Finally, the iteration complexity results specify the precise number of data points needed in order to achieve a specified performance, making the algorithm readily applicable to any real-world application.

III-C Performance Analysis

This section provides the performance analysis of the proposed algorithm. We begin with describing the surrogate problem

[TABLE]

for a given value of $\gamma$. The subsequent analysis utilizes the following key steps (a) bounding the optimality gap $F({\hat{{\mathbf{x}}}})-F({\mathbf{x}}^{\gamma})$ and the penalty $L({\hat{{\mathbf{x}}}})$; (b) relating the solutions ${\mathbf{x}}^{\star}$ and ${\mathbf{x}}^{\gamma}$; and (c) obtaining the corresponding bounds on the optimality gap $F({\hat{{\mathbf{x}}}})-F({\mathbf{x}}^{\star})$ and the constraint violation $\sum_{j}\max\{Q_{j}({\hat{{\mathbf{x}}}}),0\}$.

For the first step, we establish the following Lemma, whose proof is provided in Appendix A.

Lemma 1*.*

Under Assumptions A1-A5, the iterates generated by Algorithm 1 satsify

[TABLE]

The constants within ${\mathcal{O}}(\cdot)$ in (23) depend only on the $D_{x},C_{q},C_{h},C_{g},J,C_{\ell},L_{f},L_{q},V_{g},$ and $V_{h}$. The proof of Lemma 1 follows along the lines of the SCGD proof in [15], with appropriate modifications incorporated to handle the two variables ${\mathbf{y}}_{t}$ and ${\mathbf{z}}_{t}$, and the required changes in obtaining the bound. Specifically, the proof requires us to re-define appropriate sequences that allow the telescopic sums to be taken.

For the next step, we relate the solutions of (22) and ($\mathcal{P}$). The proof is provided in Appendix B.

Lemma 2*.*

Under Assumption (A2) and for $0<\gamma<\sigma/2$, it holds that $F({\mathbf{x}}^{\gamma})-F({\mathbf{x}}^{\star})\leq 2\sqrt{C_{g}C_{f}}D_{x}\gamma/\sigma$.

Lemma 2 follows from the strong duality of (22) and from standard results in duality theory. An implication of Lemma 2 is that the gap between the solutions of ($\mathcal{P}$) and (22) cannot be large if $\gamma$ is small. Having stated the preliminary results, we now state the main result of the paper as the following theorem. For the sake of brevity, let us define

[TABLE]

where $D_{1}$ and $D_{2}$ are constants as defined in Appendix A.

Theorem 1*.*

Let the step-sizes be chosen as

either diminishing with $\alpha_{t}=t^{-a}$, $\beta_{t}=t^{-b}$, $\delta_{t}=t^{-c}$ for $1\leq t\leq T$; 2. 2.

or constant with $\alpha_{t}=T^{-a}$, $\beta_{t}=T^{-b}$, $\delta_{t}=T^{-c}$ for all $1\leq t\leq T$,

with $1>a\geq c\geq b>0$. Under (A1)-(A5) and for $C_{\ell}$ chosen such that $\max_{{\mathbf{x}}\in{\mathcal{X}}}Q_{j}({\mathbf{x}})\leq C_{\ell}-\gamma$, we have the following bounds

[TABLE]

where $\omega={\mathcal{O}}\left({T^{a-1}+T^{a-2c}+T^{a-2b}+T^{a+b-2c}}\right)$.

Proof:

The proof is written in two parts. We first begin with analyzing the optimality gap and subsequently bound the constraint violation.

Since the right-hand side of (23) is equal to $\omega T/2$ (cf. (68)), we have that

[TABLE]

Since ${\mathbb{E}}\left[{{L({\mathbf{x}}_{t})}}\right]\geq 0$, it follows that $\sum_{t=T/2}^{T}{\mathbb{E}}\left[{{F({\mathbf{x}}_{t})-F({\mathbf{x}}^{\gamma})}}\right]\leq\omega\frac{T}{2}$. Therefore the optimality gap can be bounded as

[TABLE]

Since $F$ is convex, the bound in (26) implies that ${\mathbb{E}}\left[{{F({\hat{{\mathbf{x}}}})}}\right]-F^{\star}\leq\mathcal{O}(\omega+\gamma)$ where both $\omega$ and $\gamma$ may depend on $T$.

Next observe that since $F$ is Lipschitz with parameter $\sqrt{C_{f}C_{g}}$ and the set ${\mathcal{X}}$ is compact, it follows that $F({\mathbf{x}}_{t})-F({\mathbf{x}}^{\gamma})\geq-\sqrt{C_{f}C_{g}}\left\|{\mathbf{x}}_{t}-{\mathbf{x}}^{\gamma}\right\|\geq-\sqrt{C_{f}C_{g}}D_{x}$. Substituting into (25) and rearranging yields

[TABLE]

Since the sequence $\{\frac{\alpha_{t}}{\delta_{t}}\}$ is either diminishing or constant, it holds that $\frac{\alpha_{t}}{\delta_{t}}\leq\frac{\alpha_{T/2}}{\delta_{T/2}}$ for all $T/2\leq t\leq T$. Recalling that ${\hat{{\mathbf{x}}}}=\frac{2}{T}\sum_{t=T/2}^{T}{\mathbf{x}}_{t}$ and that the function $L$ is convex, the left-hand side of (27) can be lower bounded as

[TABLE]

Recall that $Q$ is bounded and therefore there exists $C_{\ell}$ such that $\max_{j}Q_{j}({\mathbf{x}})\leq C_{\ell}-\gamma$ for $\gamma\in(0,1)$ or equivalently, $L({\hat{{\mathbf{x}}}})=\sum_{j}[Q_{j}({\hat{{\mathbf{x}}}})+\gamma]_{+}^{2}$. Therefore, the lower bound in (28) becomes

[TABLE]

Substituting (29) into (27) and taking square root, we obtain

[TABLE]

Finally the required order bound on $\omega$ can be obtained simply by substituting the constant/diminishing step-sizes and ignoring all constant terms. ∎

Theorem 1 provides a general bound on the optimality gap and the constraint violation that depends on the choice of constants $a$, $b$, $c$, and $\gamma$. For instance, it can be seen that if $\gamma=\sqrt{\frac{JT^{c-a}\left({\omega+\sqrt{C_{f}C_{g}}D_{x}}\right)}{2^{c-a}}}$ the proposed algorithm incurs zero constraint violation. By carefully choosing the parameters $a$, $b$, and $c$, it is then possible to minimize the bounds on the optimality gap and the constraint violation, as shown in Table I. As expected, the optimality gap decays slowly when stricter requirements are imposed on the constraint violation. It is remarked that the SCGD algorithm proposed in [15] for unconstrained problems also incurs an optimality gap of ${\mathcal{O}}(T^{-1/4})$, implying that the bounds obtained are tight in the limiting case.

IV Simulation Results

This section provides numerical tests for all examples discussed in Sec. II-B. For each example, Algorithm 1 is implemented over a network of $N=3$ parallel queues for $\psi_{i}(x)=\bar{\psi}_{i}\log(x)$, and $\varphi_{i}(x)=\bar{\varphi}_{i}x$ for $1\leq i\leq N$. We use the values $\bar{\psi}_{1}=1$, $\bar{\psi}_{2}=1.5$, $\bar{\psi}_{3}=2$, $\bar{\varphi}_{1}=10$, $\bar{\varphi}_{2}=15$, and $\bar{\varphi}_{3}=20$ in all the examples. The step sizes are chosen as constant with the selection $a=0.9167$, $b=0.5$, and $c=0.75$. For all cases, the evolution of the optimality gap $F({\hat{{\mathbf{x}}}})-F^{\star}$ as well as the constraint violation $\max_{j}Q_{j}({\hat{{\mathbf{x}}}})$ with iterations is shown.

IV-A Example 1: Parallel M/G/1 Queues

For the problem considered in Example 1, it is clear that outer functions $f$ and $q$ are differentiable (cf. (2) and (3)). Likewise, the inner functions ${\mathbf{g}}$ and ${\mathbf{h}}$ are continuous, as required in (A1). Further, observe that both the objective function as well as the constraint functions are convex in ${\bm{\lambda}}$. In its current form however, the gradient of the outer functions $f$ and $q$ is unbounded. To this end, let $L^{\max}$ denote the maximum length of the packet and let $\lambda^{\max}$ be such that $\lambda^{\max}\leq\epsilon C/L^{\max}$ for some given constant $\epsilon\in(0,1)$. Such a requirement ensures that arrival rate is not arbitrarily close to the system capacity.

Given $0<\sigma<D^{\max}$, it can be seen that any $\lambda_{i}<2C_{i}^{2}(D^{\max}-\sigma)/({\mathbb{E}}L^{2}_{i}+2C_{i}(D^{\max}-\sigma){\mathbb{E}}L_{i})$ is a strictly feasible solution to (3), as required in (A2). Next, (A3) holds if the packet lengths are i.i.d. and (A4) holds with $C_{g}=V_{g}=C_{h}=V_{h}=\sum_{i=1}^{N}\lambda^{\max}\left({{\mathbb{E}}L^{2}_{i}+{\mathbb{E}}L^{4}_{i}}\right)$ where the right-hand side is bounded since $L^{\max}<\infty$. Finally, (A5) holds with constants

[TABLE]

It can be seen that the constants depend on the specified value of $\epsilon$. In other words, attempt to operate the system very close to the boundary may lead to larger constants and slower rates of convergence.

Fig. 2 considers an example system with $D^{\max}=50ms$, $C=[100~{}200~{}500]$ kbps, $L^{\max}=[20~{}30~{}60]$ kb, $\lambda^{\min}=0.1$, $\lambda=15$, $\lambda^{\max}=[5~{}7~{}9]$, and mean packet lengths $\bar{L}=[15~{}20~{}35]$ bits. The simulation is run for different packet length distributions, all truncated to respect the maximum packet length requirements. It is remarked that Algorithm 1 needs one sample of $L_{i}$ at each time instant, in addition to the memory required to store the different iterates. As can be seen from Fig. 2, $100$ samples are enough for the algorithm to converge. Note that the problem is simpler if the distribution of packet lengths is known and $F^{\star}$ can be directly calculated as a function of the first and second moments of $L_{i}$. As expected, the optimality gap and the constraint violations decay with iterations.

IV-B Example 2: Parallel M/G/1 queues in wireless networks (ergodic capacity)

It can be seen that the inner functions ${\mathbf{g}}$ and $h$ are continuous as required in (A1). For ease of analysis, let us assume that the support of the random variable $\zeta_{i}$ is bounded away from zero, i.e., $\zeta_{i}\geq G$. Coupled with the fact that $p_{i}\geq P^{\min}$, it follows that $b_{i}(p_{i},\zeta_{i})>B_{i}\log(1+P^{\min}G)$. Therefore, it is always possible to find $p_{i}\in(P^{\min},P^{\max})$ such that ${\mathbb{E}}\left[{{\min_{i}b_{i}(p_{i},\zeta_{i})}}\right]>R^{\min}-\sigma$ as required in (A2). Further, since ${\mathbf{g}}$ is locally Lipschitz on compact sets ${\bm{\Lambda}}$ and ${\bm{\Pi}}$, it follows that both ${\mathbf{g}}({\bm{\lambda}},{\mathbf{p}},{\bm{\zeta}})$ and $\nabla{\mathbf{g}}({\bm{\lambda}},{\mathbf{p}},{\bm{\zeta}})$ are bounded for all ${\mathbf{p}}\in{\bm{\Pi}}$ and ${\bm{\lambda}}\in{\bm{\Lambda}}$. Indeed, as required in (A4), we have that

[TABLE]

$C_{h}=V_{h}=\frac{\max_{i}{B_{i}}}{\left({1+P^{\min}G}\right)^{2}}$ where $B_{i}$ is the capacity of the $i$-th channel.

Observe that the outer function $f$ as defined in Example 2 is unbounded and has a restricted domain. Instead, we redefine the outer functions as $\tilde{f}$ that approximates $f$ while also satisfying the regularity conditions. Let $\tilde{f}({\mathbf{y}})=\sum_{i}\varphi_{i}(\varsigma_{i}(y_{i+2N},y_{i+N}))-\psi(y_{i})$ where

[TABLE]

Note that since $b_{i}(p_{i},\zeta_{i})>B_{i}\log(1+P^{\min}G)$, the iterate ${\mathbf{y}}_{t}$ takes values in a compact domain, and consequently, it suffices to ensure that $\tilde{f}({\mathbf{y}})$ is locally Lipschitz, implying that $\nabla\tilde{f}({\mathbf{y}})$ is bounded. As required in (A5), we have that

[TABLE]

As $q(z)=z$, it is easy to see that $L_{q}=0$, and $C_{q}=1$. Finally, (A3) holds if the channel gains are i.i.d.

Verifying the convexity of the redefined objective function is challenging, as the expression for the Hessian is quite complicated. Instead, we resort to a numerical verification by observing that the objective function is smooth in ${\bm{\lambda}}\in{\bm{\Lambda}}$ and ${\mathbf{p}}\in{\bm{\Pi}}$. The Hessian for a summand of $\tilde{f}$ is evaluated over a two dimensional grid in ${\bm{\Lambda}}\times{\bm{\Pi}}$ and is found to be always strictly positive definite for the parameters considered in the simulations. Numerical verification is deferred to Appendix C.

For experiments we consider a setting with $B_{i}=10$ kbps for $i=1,\ldots,N$, $\lambda=37$, $\epsilon=0.95$, $G=0.25$, $R^{\min}=35$ kbps, $P^{\min}=14$ units, $P^{\max}=100$ units, $\lambda^{\min}=0.1$, $\lambda^{\max}=15$. Rayleigh fading is assumed so that $\zeta_{i}$ are chi-squared distributed with support truncated to $G$. Two sets of simulations are carried out with $5$, and $10$ antennas at the BS, and the performance is shown in Fig. 3. It can be seen that the optimality gap decays to zero as expected. Interestingly, observe that the constraints are never violated when the number of antennas are large, likely due to the diversity gains available to the system. Different from Example 1 however, $F^{\star}$ cannot be obtained in closed form but requires brute force search, since the moments of the inner function are themselves functions of the powers. From the Fig. 3, it is evident that for this example, 300 samples are required to obtain a feasible solution that is close to optimum.

IV-C Example 3: Parallel M/G/1 queues in wireless networks (outage capacity)

The objective function in Example 3 is not convex, and therefore the results in Theorem III do not apply. Nevertheless, the results from [15] are still applicable under (A1), (A3), (A4), and (A5). First, note that ${\bm{\Lambda}}$ and ${\bm{\Pi}}$ are defined as in Example 2. The channel is also assumed to have a support bounded away from zero, i.e., $\zeta_{i}\geq G$. Since the domain of $f$ is restricted, it is possible to utilize an approximate version as in Example 2. Alternatively, if the problem parameters satisfy $\epsilon R_{i}\left({1-r_{i}(P^{\min},G)}\right)\geq\lambda^{\max}$ for some $\epsilon$, then the outer function’s gradient is bounded. In the current form, inner function $g$ (see. (7)) is also not continuous as required in (A1). To this end we approximate the indicator function in (6) with the sigmoid function $r_{i}(p_{i},\zeta_{i})=1/\left({1+e^{-\eta\left({R_{i}-b_{i}(p_{i},\zeta_{i})}\right)}}\right)$ where $\eta\geq 1$ is a parameter that controls the tightness of the approximation. Outer function $f$ in (8) is clearly continuously differentiable (as required in (A1)) and (A3) holds if the channel gains are i.i.d. Further (A4) holds with

[TABLE]

Since the domain of $({\mathbf{p}},{\bm{\lambda}})$ is compact and $f$ is locally Lipschitz, appropriate constants $C_{f}$ and $L_{f}$ can likewise be found for (A5) as well.

For experiments we consider a setting with $B_{i}=100$ kbps, $\lambda=45$, $\epsilon=0.95$, $G=0.25$, $R=[30,35,40]$ kbps, $P^{\min}=10$ units, $P^{\max}=100$ units, $\lambda^{\min}=0.1$, and $\lambda_{i}^{\max}=25$. Channel gains are exponentially distributed with mean 1 and the base station has a single antenna. The SCGD algorithm [15] is implemented for this example and the simulations are shown in Fig. 4. For this case, the algorithm is only guaranteed to converge to a stationary point. However, in order to also evaluate the quality of the solution, we minimize the objective using the trust-region reflective method as implemented in the fmincon function of MATLAB. It is remarked that solving the problem using such an algorithm is generally impractical since it requires carrying out simulations in order to calculate the expected gradient and Hessian at every iteration. As evident from Fig. 4, it can be seen that the objective function values converge to a local minimum value $\tilde{F}$ obtained from MATLAB, while incurring significantly lower computational cost.

IV-D Example 4: Parallel G/G/1 queues in wireless networks (effective capacity)

Similar to the previous example, the objective function here is not convex, and the results from [15] apply under assumptions (A1), (A3), (A4), and (A5). Channel gains at consecutive update intervals are considered to be i.i.d., and therefore (A3) holds. As defined in Example 2, 3, $\Pi(\cdot)$ acts as the ball constraint, and because of which the domain of $f$ is restricted. As required in (A5), the outer function’s gradient is bounded as long as ${\mathbb{E}}\left[{{b_{i}(P^{\max},\zeta_{i})}}\right]$ and ${\mathbb{E}}\left[{{\left({b_{i}(P^{\max},\zeta_{i})}\right)^{2}}}\right]$ are finite. Additionally, outer function $f$ is differentiable as required in (A1). Further (A4) holds with

[TABLE]

For experiments, we consider a setting with $B_{i}=100$ kbps, $P^{\min}=0.1$ units, $P^{\max}=0.9$ units, and $W=500$ ms. Channel gains are exponentially distributed with means $[0.8,0.9,1]$, and the base station has a single antenna. Simulations for this example are shown in Fig. 5. Since the problem is non convex, the algorithm [15] is only guaranteed to converge to a stationary point. The local optimum objective $\tilde{F}$ is calculated by using the approximations in [10]. Since the objective function is of composition form, the inner functions ${\mathbb{E}}\left[{{b_{i}(p_{i},\zeta_{i})}}\right]$, ${\mathbb{E}}\left[{{\left({b_{i}(p_{i},\zeta_{i})}\right)^{2}}}\right]$ are approximated to a deterministic functions of $p_{i}$, and ${\mathbb{E}}\left[{{\zeta}}\right]$, for Rayleigh channels. In other words, composition form can be avoided if the channel gain $\zeta_{i}$ is of exponentially distributed. However the approximations are of some complicated functions, and valid only for Rayleigh channels. In contrast, the SCGD algorithm is applicable to any channel with arbitrary distributions. As evident from Fig. 5, the objective function value converges to a local minimum $\tilde{F}$.

V Conclusion

The problem of queuing system design is formulated as a constrained stochastic optimization problem. As the objective and the constraints of the problem are non-linear functions of expectations, classical stochastic approximation methods cannot be used. Instead a constrained variant of the stochastic compositional gradient descent method is proposed where the expected value functions are tracked explicitly. The sample complexity of the proposed algorithm is characterized and it is shown that after $T$ iterations, the average optimality gap and constraint violation are at most ${\mathcal{O}}\left({T^{-1/12}}\right)$. Alternative choices of the parameters also yield zero constraint violation but a larger optimality gap of ${\mathcal{O}}\left({T^{-1/12}}\right)$ or an ${\mathcal{O}}(1)$ constraint violation and the best possible optimality gap of ${\mathcal{O}}\left({T^{-1/4}}\right)$. The proposed algorithm is tested over the different queuing design examples and it is shown that the optimality gap and the constraint violation decay as expected.

Appendix A Proof of Lemma 1

The overall proof is split into four lemmas. For the sake of brevity, let us define the function $p({\mathbf{z}}):=\ell\left({{\mathbf{q}}({\mathbf{z}})}\right)$. The following preliminary lemma establishes the necessary properties of $p$.

Lemma 3*.*

The function $p({\mathbf{z}})$ has Lipschitz continuous and bounded gradients,

[TABLE]

Proof:

From the chain rule, we have that ${\nabla}p({\mathbf{z}})={\nabla}{\mathbf{q}}({\mathbf{z}}){\nabla}\ell({\mathbf{q}}({\mathbf{z}}))$. Using triangle inequality, Assumption (A5), and the definition of $\ell$, we obtain $\left\|{\nabla}p({\mathbf{z}})\right\|^{2}\leq JC_{q}C_{\ell}^{2}$ for all ${\mathbf{z}}$. In a similar vein, it can be seen that $p$ is smooth since both $\ell$ and ${\mathbf{q}}$ are smooth functions. In particular, we have that

[TABLE]

Taking norm, and applying the Cauchy-Schwartz inequality, we obtain

[TABLE]

where we have used (A5) and the implication that ${\mathbf{q}}$ is Lipschitz continuous with parameter $\sqrt{C_{q}}$. ∎

Building upon Lemma 3, we now characterize the consecutive iterate difference and show that it is ${\mathcal{O}}(\alpha^{2}_{t}+\delta_{t}^{2})$.

Lemma 4*.*

The updates in (21) yields

[TABLE]

for all $t\geq 1$.

Proof:

From Algorithm 1, since ${\mathbf{x}}_{t}\in{\mathcal{X}}$, the non-expansiveness property of the projection operation implies that

[TABLE]

where (35) and (36) follow from the Peter-Paul and Cauchy-Schwartz inequalities, respectively. The last inequality follows from Lemma 3 and (A5). Taking expectations and using (A4), we obtain

[TABLE]

which is the required result. ∎

Next, we establish a key result that legitimizes the tracking steps in (20). To this end, define functions ${\bar{{\mathbf{g}}}}({\mathbf{x}}):={\mathbb{E}}_{{\bm{\zeta}}}[{\mathbf{g}}({\mathbf{x}},{\bm{\zeta}})]$ and ${\bar{{\mathbf{h}}}}({\mathbf{x}}):={\mathbb{E}}_{{\bm{\zeta}}}[{\mathbf{h}}({\mathbf{x}},{\bm{\zeta}})]$ for all ${\mathbf{x}}\in{\mathcal{X}}$. Also, let ${\mathcal{F}}_{t}$ be the sigma algebra generated by the random variables $\{{\bm{\zeta}}_{1},\ldots,{\bm{\zeta}}_{t-1}\}$ so that ${\mathbb{E}}\left[{{{\mathbf{g}}({\mathbf{x}}_{t},{\bm{\zeta}}_{t})\mid{\mathcal{F}}_{t}}}\right]={\bar{{\mathbf{g}}}}({\mathbf{x}}_{t})$ from (A3).

Lemma 5*.*

For the choice $\alpha_{t}\leq\delta_{t}$ and $\delta_{t}\leq\beta_{t}$, the updates (20) in Algorithm 1 yield

[TABLE]

where $D_{y}$ and $D_{z}$ are initialization dependent constants.

Proof:

We provide the proof of (40) while observing that (41) follows along similar lines. Let ${\mathbf{e}}_{t}:=(1-\beta_{t})\left({\bar{{\mathbf{g}}}}({\mathbf{x}}_{t})-{\bar{{\mathbf{g}}}}({\mathbf{x}}_{t-1})\right)$ be the difference between evaluations of ${\bar{{\mathbf{g}}}}$ at consecutive iterates, scaled by $(1-\beta_{t})$. From the Peter-Paul inequality, we have that

[TABLE]

Since ${\bar{{\mathbf{g}}}}$ is Lipschitz (cf. (A4)), the second term of (42) can be bounded as

[TABLE]

For the first term of (42), we have that

[TABLE]

Taking squared norm, applying conditional expectation given ${\mathcal{F}}_{t}$, and using the definition of ${\bar{{\mathbf{g}}}}$, we obtain

[TABLE]

where the last inequality follows from (A4). Taking full expectation in (45) and plugging into (42), we obtain the required result as

[TABLE]

where we have also used the fact that $\beta_{t}<1$. To establish the upper bound, take $D_{y}:={\mathbb{E}}\left[{{\left\|{\mathbf{y}}_{2}-{\bar{{\mathbf{g}}}}({\mathbf{x}}_{1})\right\|^{2}}}\right]+2V_{g}+2C_{g}\left({C_{f}C_{g}+JC_{\ell}^{2}C_{q}C_{h}}\right)$ and observe that ${\mathbb{E}}\left[{{\left\|{\mathbf{y}}_{3}-{\bar{{\mathbf{g}}}}({\mathbf{x}}_{2})\right\|^{2}}}\right]\leq(1-\beta_{t})D_{y}+2V_{g}\beta_{t}^{2}+2\frac{\alpha_{t}^{2}}{\beta_{t}}C_{f}C_{g}^{2}+2\frac{\delta_{t}^{2}}{\beta_{t}}JC_{\ell}^{2}C_{g}C_{q}C_{h}\leq D_{y}$ since $\alpha_{t}/\beta_{t}\leq 1$ and $\delta_{t}/\beta_{t}\leq 1$. The argument can inductively be applied to all values of $t$.

The derivation of (41) follows along similar lines. For the upper bound, we set

$D_{z}={\mathbb{E}}\left[{{\left\|{\mathbf{z}}_{2}-{\bar{{\mathbf{h}}}}({\mathbf{x}}_{1})\right\|^{2}}}\right]+2V_{h}+2C_{h}\left({C_{f}C_{g}+JC_{\ell}^{2}C_{q}C_{h}}\right)$ and proceed inductively as before. ∎

Next, we analyze the distance from the optimal set $\left\|{\mathbf{x}}_{t}-{\mathbf{x}}^{\gamma}\right\|^{2}$, and establish its relation with the optimality gap and the constraint violation.

Lemma 6*.*

For convex functions $F$ and ${\mathbf{Q}}$, we have the following bound.

[TABLE]

Proof:

Since ${\mathbf{x}}^{\gamma}\in{\mathcal{X}}$, the non-expansiveness of the $\Pi_{{\mathcal{X}}}$ operator implies that

[TABLE]

where $u_{t}=2\alpha_{t}\left({{\mathbf{x}}_{t}-{\mathbf{x}}^{\gamma}}\right)^{T}{\widetilde{\nabla}}{\mathbf{g}}({\mathbf{x}}_{t},{\bm{\zeta}}_{t})\left({{\nabla}f({\bar{{\mathbf{g}}}}({\mathbf{x}}_{t}))-{\nabla}f({\mathbf{y}}_{t+1})}\right)$, and

$v_{t}=2\delta_{t}\left({{\mathbf{x}}_{t}-{\mathbf{x}}^{\gamma}}\right)^{T}{\widetilde{\nabla}}{\mathbf{h}}({\mathbf{x}}_{t},{\bm{\zeta}}_{t})\left({{\nabla}p({\bar{{\mathbf{h}}}}({\mathbf{x}}_{t}))-{\nabla}p({\mathbf{z}}_{t+1})}\right)$. Taking conditional expectation given ${\mathcal{F}}_{t}$ in (55), we obtain

[TABLE]

where (57) follows from Lemma 4 and the convexity of functions $F$ and $L$. Next, noting that $L({\mathbf{x}}^{\gamma})=0$ and taking full expectation in (57), we obtain

[TABLE]

In order to bound $u_{t}$, we utilize the Cauchy-Schwartz and Peter-Paul inequalities as

[TABLE]

where (62) follows from (A2). Taking expectation on both sides of (62), we obtain $E{u_{t}}\leq\beta_{t}{\mathbb{E}}\left[{{\left\|{\mathbf{y}}_{t+1}-{\bar{{\mathbf{g}}}}({\mathbf{x}}_{t})\right\|^{2}}}\right]+\frac{\alpha_{t}^{2}}{\beta_{t}}L_{f}D_{x}C_{g}$, and we can also obtain ${\mathbb{E}}\left[{{v_{t}}}\right]\leq\beta_{t}{\mathbb{E}}\left[{{\left\|{\mathbf{z}}_{t+1}-{\bar{{\mathbf{h}}}}({\mathbf{x}}_{t})\right\|^{2}}}\right]+\frac{\delta_{t}^{2}}{\beta_{t}}\left({\sqrt{J}L_{q}C_{\ell}+C_{q}}\right)D_{x}C_{h}$. The required result follows from substituting into (58). ∎

Having established the required recursive relationships, we proceed with applying telescopic sums in order to obtain the required bounds. Let

[TABLE]

and observe that $I_{t}\leq D_{x}+D_{y}+D_{z}=:D_{1}$. Multiplying (40) and (41) with $\left({1+\beta_{t}}\right)$ and substituting in (49), we obtain

[TABLE]

where $D_{2}$ is the maximum of $4\left({C_{g}+C_{h}}\right)C_{f}C_{g}+4\left({C_{g}+C_{h}}\right)JC_{\ell}^{2}C_{q}C_{h}\\ +\left({L_{f}C_{g}+\sqrt{J}L_{q}C_{h}C_{\ell}+C_{h}C_{q}}\right)D_{x}$, $2\left({C_{f}C_{g}+C_{\ell}^{2}C_{q}C_{h}}\right)$, and $4\left({V_{g}+V_{h}}\right)$. Summing over

$t=T/2$ to $T$, noting that $\alpha_{t}\leq\alpha_{T/2}$, and rearranging, we obtain the required result as

[TABLE]

where we have dropped the negative terms from the right-hand side. which is the required result.

Appendix B Proof of Lemma 2

Let us associate dual variable ${\bm{\lambda}}$ with the constraint in (22) and let ${\bm{\lambda}}^{\gamma}\geq 0$ be the dual optimal solution. Since $\gamma<\sigma/2<\sigma$, there exists a strictly feasible $\widetilde{{\mathbf{x}}}$ such that $Q_{j}(\widetilde{{\mathbf{x}}})+\gamma<Q_{j}(\widetilde{{\mathbf{x}}})+\sigma\leq 0$, and consequently strong duality holds for (22). Therefore, we have that

[TABLE]

where we have used the fact that ${\mathbf{Q}}({\mathbf{x}}^{\star})\leq 0$. Let $\widetilde{{\mathbf{x}}}$ be a strictly feasible point specified in (A2). Then from (69), it also follows that

[TABLE]

or equivalently

[TABLE]

where the last inequality follows since $F$ is Lipschitz continuous with parameter $\sqrt{C_{g}C_{g}}$ and since ${\mathcal{X}}$ is compact. The required bound follows from substituting (71) into (69).

Appendix C Numerical verification of convexity of the formulation in Example 2

Lemma 7*.*

The following function is convex with respect to $\lambda$, and $p$ for the box constraints $0.1\leq\lambda\leq 15$, and $14\leq p\leq 100$

[TABLE]

where $\zeta$ is a Chi-squared distributed random variable with $2K$ degrees of freedom.

Proof:

Since the expression in (72) is of non linear expectation function, it is more difficult to get an analytical proof. Instead, we try to derive the Hessian of (72) numerically, and show that the minimum eigenvalue of Hessian is nonnegative. To do that, we take help of two MATLAB functions namely ‘numerical gradient’, and ‘numerical integration’. The results in Fig. 6 confirm that the Hessian of the expression in (72) is positive semidefinite matrix.

∎