Optimal Power Procurement for Green Cellular Wireless Networks Under Uncertainty and Chance Constraints

TL;DR

This paper proposes a new method for optimizing energy procurement in cellular networks using renewable energy while ensuring service quality.

Contribution

A novel time-continuous Lagrangian relaxation approach is introduced to handle probabilistic QoS constraints in energy procurement.

Findings

The proposed method efficiently solves the stochastic optimal control problem with probabilistic constraints.

Numerical results show the approach is computationally efficient and provides near-optimal solutions in practical timeframes.

Abstract

Given the increasing global emphasis on sustainable energy usage and the rising energy demands of cellular wireless networks, this work seeks an optimal short-term, continuous-time power-procurement schedule to minimize operating expenditure and the carbon footprint of cellular wireless networks equipped with energy-storage capacity, and hybrid energy systems comprising uncertain renewable energy sources. Despite the stochastic nature of wireless fading channels, the network operator must ensure a certain quality-of-service (QoS) constraint with high probability. This probabilistic constraint prevents using the dynamic programming principle to solve the stochastic optimal control problem. This work introduces a novel time-continuous Lagrangian relaxation approach tailored for real-time, near-optimal energy procurement in cellular networks, overcoming tractability problems associated with the probabilistic QoS constraint. The numerical solution procedure includes an efficient upwind finite-difference solver for the Hamilton–Jacobi–Bellman equation corresponding to the relaxed problem, and an effective combination of the limited memory bundle method (LMBM) for handling nonsmooth optimization and the stochastic subgradient method (SSM) to navigate the stochasticity of the dual problem. Numerical results, based on the German power system and daily cellular traffic data, demonstrate the computational efficiency of the proposed numerical approach, providing a near-optimal policy in a practical timeframe.

Figures

Full text

1. Introduction

Since 2021, worldwide mobile broadband traffic has increased by an average of 19.6% annually, reaching 1 zettabyte (ZB) in 2023 [1]. The 4G network coverage increased from 41% of the world population in 2015 to 92% of the world population in 2024 [1]. The estimated transmission network electricity consumption in 2022 was 240 to 340 TWh, an increase of 64% from 2015 [2]. Next-generation transmission technology (e.g., 5G) requires even more power to ensure sufficient mobile network coverage [3]. A consequence of the high energy demand is increased greenhouse gas emissions. Data-transmission networks accounted for around 330 Mt  $[eqn]$ -equivalent in 2020, around 0.9% of all energy-related emissions [2]. About 83% of the energy consumed by cellular networks came from fossil fuels in 2021, whereas only 9% came from renewables [4]. This problem motivates the coupling of cellular wireless networks with renewable energy to reduce the costs and carbon footprint of the telecommunications industry. Hybrid energy systems provide a consistent power supply to wireless networks by combining renewable energy sources (e.g., solar and wind) [5]. However, modeling these systems and deriving an optimal power-procurement policy can be challenging, due to the stochasticity of solar/wind power and wireless channels [6].

Previous studies have considered the optimal power-procurement problem in a discrete-time setting for isolated energy systems under an uncertain energy demand and incoming renewable energy [7,8,9]. The corresponding optimization problem for wireless networks has additional difficulties due to the uncertainty in wireless channels and probabilistic constraints due to the network quality of service (QoS) [10,11,12,13]. The solution method for this problem in the current literature can be generally categorized into robust optimization and stochastic optimization. The robust optimization approach constructs a bounded uncertainty set encompassing all realizations of the stochastic variables and optimizes the worst-case objective function value [14,15,16,17]. Robust optimization, while ensuring feasibility under the worst-case scenario, often leads to overly conservative policies that could underutilize renewable energy resources. The stochastic optimization approach assumes a known probability distribution for the random variables and optimizes the expected value of the objective function over the probability distribution of the uncertain parameters [18,19].

Chance-constrained programming is a subclass of stochastic optimization methods, first introduced in [20,21]. This approach allows constraints to be violated with a certain prescribed probability or risk level. The standard approach to solving discrete chance-constrained problems is to replace the chance constraint with a conservative but tractable approximation using the Bernstein inequality [22,23]. This approach is well-suited for wireless networks, since the associated QoS constraint is commonly defined probabilistically, allowing the operator to balance between conservatism and optimality. However, the computational complexity of the method in the discrete-time setting increases exponentially for large datasets and complex problems. For example, the authors in [24,25,26] solved the resulting discrete optimization problem using heuristic techniques, yielding suboptimal solutions.

This work proposes a novel mathematical modeling and numerical framework for the optimal management of hybrid energy systems to power cellular wireless networks under uncertainty and chance constraints in continuous time for a short-term planning horizon, where the system comprises locally installed renewable energy sources, fossil fuel power stations, and a storage capacity modeled by a single battery. We employ stochastic differential equations (SDEs) to model the incoming instantaneous renewable power and wireless fading channels. Data-driven SDEs have been employed to model instantaneous wind [27,28,29,30] and solar power [31,32,33]. These SDE models follow a mean reversion to deterministic wind speed/solar irradiance forecasts. The parameters of the SDE characterizing its uncertainty are calibrated from historic discrepancies (data) between the actual and forecasted quantities. Moreover, we employ SDEs to model stationary wireless fading channels. The authors in [34,35,36,37] constructed and configured SDEs to reproduce characteristics of various fading channel models, primarily concerning their stationary distribution.

This work formulates and solves a time-continuous stochastic optimal control problem to derive a near-optimal policy minimizing the expected operating expenditure and carbon footprint of cellular base stations. We optimize for both cost minimization and environmental impact, introducing a trade-off parameter  $[eqn]$  to balance these objectives. The optimal control is subject to demand constraints, a probabilistic and time-pointwise QoS constraint, and the dynamics and capacity limits of energy systems and wireless channels. The probabilistic constraint prevents the direct application of the dynamic programming principle because the constraint controls the joint distribution of the random variables. The optimal policy must be a function of the current state and time for the standard dynamic programming lemma to hold, not a function of its distribution [38]. To overcome this problem, we pose the corresponding dual problem by introducing a time-continuous Lagrangian relaxation to penalize violations of the probabilistic constraint [39]. This dual formulation yields a standard time-continuous stochastic optimization problem solved iteratively by optimizing the dual function. Each iteration involves solving the Hamilton–Jacobi–Bellman (HJB) partial differential equation (PDE) [40,41] using an upwind finite-difference scheme [42] to compute the dual objective value and subgradient. With access to only noisy subgradients, we solve the convex, nonsmooth dual problem using a novel combination of the limited memory bundle method (LMBM) [43,44,45] and the stochastic subgradient method [46,47]. Moreover, we design a refinement strategy for the Lagrangian multipliers to control the pointwise violation in the chance constraint. The numerical experiments and results based on the German power system and daily cellular traffic profiles validate the efficiency of the proposed models and approach.

The contributions can be summarized as follows:

1. We propose a novel time-continuous optimization framework for optimal power procurement for green wireless cellular networks, subject to SDE dynamics and chance constraints. Compared to the discrete-time formulation in previous studies [24,25,26], the proposed approach decouples the model development from the numerical approximation, enhancing the model fidelity (see Remark 3). This formulation also yields a continuous control curve over time, allowing its application for any time discretization scheme and eliminating the need for ad hoc interpolations [48].
2. We calibrate the data-driven SDE model developed for instantaneous wind power in [30] using German wind power data from the year 2023. The calibrated SDE is a driving dynamic for the stochastic optimal control problem.
3. We derive an SDE to model the instantaneous Nakagami wireless fading channel [49] with a shifted-gamma invariant distribution, using the Pearson class of diffusions [50]. This SDE is a driving dynamic for the stochastic optimal control problem.
4. We apply Lagrangian relaxation to the probabilistic QoS constraint, transforming the problem into a standard time-continuous stochastic optimization problem. Studies have explored numerical methods for time-continuous stochastic optimization with final-time chance constraints using Lagrangian relaxation [51,52,53] or reformulation as a stochastic target problem [54,55]. However, the proposed approach is novel in addressing a chance constraint that must be satisfied at every time point. Moreover, we implement this within the context of cellular wireless networks.
5. We develop an iterative algorithm to optimize the dual function within a finite-dimensional function class numerically. Each iteration involves solving the HJB PDE to compute the dual function value and its noisy subgradient. The proposed approach extends the work in [48] on the time-continuous deterministic optimization of coupled hydrothermal power systems to the stochastic setting.
6. With access only to noisy realizations of the subgradient, we combine the LMBM method [56] with the stochastic subgradient method [47] to solve the dual optimization problem. We iteratively refine the Lagrange multiplier to ensure compliance with the chance constraint at every time point.

The outline of the paper is as follows. Section 2 introduces the complete system model, detailing the costs, power generation, wireless network, and state dynamics. Next, Section 3 formulates the primal stochastic optimal control problem and the associated dual problem by Lagrangian relaxation of the chance constraint. Then, Section 4 introduces the numerical approach for solving the dual problem and presents all components of the final algorithm in detail. This section discusses the numerical solution of the HJB PDE, subgradient estimation, dual-problem optimization, and numerical errors arising from the proposed approach. Finally, Section 5 demonstrates the numerical results of applying the above approach when solving the stochastic optimization problem for a cellular base station based on the German power system and daily cellular traffic profiles.

2. System Model

This work considers a cellular network of multiple noninteracting base stations, each powered by renewable energy sources (solar panels or wind turbines) and battery storage. Figure 1 in [25] depicts the schematic network. Figure 1 illustrates the electrical power flow in each base station, each with a battery that stores the power generated from renewables ( $[eqn]$ ). The power transmitted by the base station to serve cellular users consists of power procured from the battery ( $[eqn]$ ) and bought from the grid ( $[eqn]$ ). Any extra energy stored in the battery could be sold back to the grid ( $[eqn]$ ) for revenue. We aim to optimize the energy procurement at each base station for a one-day operation cycle.

2.1. Base Station Model

The instantaneous power procured by the mobile operator from the traditional grid to power the base station is denoted by  $[eqn]$ , where 0 and T denote the time at the beginning and end of the day, respectively. The corresponding unit price is denoted by  $[eqn]$ . This price can vary during the day depending on strategies followed by the stakeholders in the energy market. Moreover, the base station has its own internal source of power (i.e., renewable power generator) whose generated power is denoted by  $[eqn]$ . This power is assumed to be free of charge. The internal power source is assumed to have a maximum power output capacity denoted by  $[eqn]$ . The base station has a battery, allowing the mobile operator to store the incoming renewable power. As the base station is interconnected with the traditional grid, the mobile operator chooses to use the stored renewable power to run the base station or sell the extra energy back to the grid for revenue. The instantaneous power drawn from the battery to run the base station is denoted by  $[eqn]$ . The instantaneous power from the battery sold back to the grid is denoted by  $[eqn]$ , and the corresponding unit price is denoted by  $[eqn]$ . We reasonably assume  $[eqn]$  for all  $[eqn]$ . The instantaneous power transmitted by the base station to serve its users is denoted by  $[eqn]$ . We assume the total transmitted power is equally divided among each average mobile user connected to the network. The instantaneous power transmitted by the base station per average user is denoted by  $[eqn]$  and  $[eqn]$  for all  $[eqn]$ , where  $[eqn]$  is the number of connected mobile users at time t. The instantaneous power balance equation for the base station is

$$[eqn]$$

where  $[eqn]$  represents a factor that scales with the transmitted power due to the amplifier and feeder losses in the network, and  $[eqn]$  models the offset power required to operate the base station irrespective of the transmitted power [25]. We assume that the base station can transmit a maximum power of  $[eqn]$ .

$$[eqn]$$

The term  $[eqn]$  in (1) is determined by the choice of the daily mobile user traffic model, heavily influencing the power demand in the base station. One such model is introduced in Section 5.1.

2.2. Cellular Network Model

According to the power law, the power transmitted by the base station decays with increased distance from the base station. Moreover, the transmitted power is uncertain due to atmospheric conditions, natural and human-made environmental obstacles, and interference from other signals [11]. Hence, the instantaneous power received by a user u at position  $[eqn]$  from a base station at  $[eqn]$  is given as  $[eqn]$ , where  $[eqn]$  denotes the path loss constant and exponent, respectively. Moreover,  $[eqn]$  denotes the Euclidean norm.  $[eqn]$  represents a stochastic process, that is greater than  $[eqn]$  almost surely (a.s.), modeling the wireless fading channel.

This work assumes the channel follows Nakagami fading [49]. The Nakagami fading model is a good data-fit model for cellular systems in urban and suburban areas [57]. In this model, the instantaneous power is gamma distributed, a Pearson type III distribution [58]. Hence, Pearson diffusion [50] is employed to model the wireless fading channel. Thus, this work constructs an SDE for the process  $[eqn]$  (3) with a mean-reverting linear drift and a squared diffusion, which is a second-order polynomial of the state, and whose invariant distribution is the shifted-gamma distribution [50] with shape parameter  $[eqn]$  and scale parameter  $[eqn]$ . We consider the probability space  $[eqn]$ , where  $[eqn]$  is the filtration of the standard Wiener process  $[eqn]$ :

$$[eqn]$$

where  $[eqn]$  is the invariant distribution of the process  $[eqn]$ .

Remark 1 (Computational efficiency of Pearson diffusions).  Pearson diffusion naturally models the continuous-time wireless fading channel, where process  $[eqn]$  lies in  $[eqn]$  for all  $[eqn]$ . This model also has the advantage of analytical tractability, providing closed-form solutions for the invariant distribution of the process  $[eqn]$  and offering better insight into its behavior. This approach makes the model choice computationally more efficient than other Markov-based models (e.g., discrete-state Markov chains), whose complexity increases quadratically as the number of possible states increases.

In this work, the network performance is characterized by the offered QoS, determined by the signal-to-noise ratio (SNR) threshold  $[eqn]$ . The SNR, in decibels (dB), at user u at  $[eqn]$  served by the base station can be expressed as follows:

$$[eqn]$$

The ambient noise level is denoted by  $[eqn]$ . The QoS of the base station is determined by the proportion of users whose SNR is above  $[eqn]$ . A user u is defined to be in outage at time t if  $[eqn]$ . The proportion of users in outage is denoted by  $[eqn]$  and can be expressed as follows:

$$[eqn]$$

where  $[eqn]$  can be interpreted as an empirical average over a user distribution. The following approximation assumes that each user is identical and that sufficient number of mobile users exist at all times:

$$[eqn]$$

where  $[eqn]$  denotes the distribution of mobile users at time t. In (6),  $[eqn]$  is the complementary cumulative distribution function (ccdf) of the user distribution. This quantity can be computed (even analytically) for many well-known distributions.

The objective of the network operator is to ensure  $[eqn]$  for all  $[eqn]$ , where  $[eqn]$  represents a threshold ratio. Due to the uncertainty associated with channel fading, achieving the objective in an a.s. sense is infeasible. Instead, we introduce the following chance constraint for the network QoS.

$$[eqn]$$

where  $[eqn]$ , for  $[eqn]$ , is the confidence level of satisfying the constraint in (7). The network operator also generates revenue from servicing the users throughout the day:

$$[eqn]$$

where  $[eqn]$  denotes the unit price paid by users to connect to the network, and  $[eqn]$  represents the number of users that are connected to the network (not in outage).

2.3. Renewable Power Model

This work also incorporates the uncertainty in renewable power generation. We apply a data-driven parametric SDE, whose solution defines a stochastic process that models the error in a given forecast for renewable power [30]. Let the process  $[eqn]$  denote the normalized generated renewable power.

$$[eqn]$$

where  $[eqn]$  is a deterministic forecast for the normalized power provided by an official source and  $[eqn]$  is its time derivative. In addition,  $[eqn]$  denotes the mean-reversion time-varying parameter ensuring that the process R is the unique strong solution to (9) for all  $[eqn]$  with range  $[eqn]$ , a.s. Moreover,  $[eqn]$  are parameters inferred from discrepancies between historical real-world observations and their forecasts. In addition,  $[eqn]$  is a standard Wiener process. The drift function in (9) is designed such that the process reverts to its mean  $[eqn]$ , with a time-varying speed  $[eqn]$  proportional to the deviation of the process  $[eqn]$  from its mean, and it tracks the time derivative  $[eqn]$ . The diffusion function in (9) is designed so that the process  $[eqn]$  avoids exiting from the range  $[eqn]$ . The net renewable power generated is then given by  $[eqn]$  for all  $[eqn]$ .

2.4. Battery Model

This work uses a simple model for the battery based on a circuit model of the state of the charge of a lithium-ion battery [59]. We assume that it has no energy loss when used and that no associated operational costs or aging effects exist. Process  $[eqn]$  denotes the normalized charge in the battery:

$$[eqn]$$

where  $[eqn]$  represents the maximum battery charge capacity,  $[eqn]$  denotes the maximum power the battery can absorb, and  $[eqn]$  indicates the maximum power that the battery can supply. These quantities are characteristic of the battery. Generally, the battery can supply more power than it can absorb; hence,  $[eqn]$  for all  $[eqn]$ . In addition,  $[eqn]$  represents the normalized charge in the battery at the start of the day. The requirement that  $[eqn]$  for all  $[eqn]$  is enforced via constraints on charging and discharging at the extremes:

$$[eqn]$$

The network operator generates revenue by selling some energy stored in the battery throughout the day:

$$[eqn]$$

Remark 2 (Battery model improvements).  Although the battery model (10) is straightforward and idealistic, realistic battery efficiency losses in the model can be incorporated by adding simple energy loss terms in the objective function of the optimal control problem. For example, this adjustment could include Ohmic efficiency losses [60] as a polynomial function of the battery discharge ( $[eqn]$ ), and state-of-charge losses [61] as an exponential function of the normalized battery charge A, which is outside the current scope of the work.

2.5. Grid Power Model

The network operator incurs the following cost of buying power from the traditional grid throughout the day:

$$[eqn]$$

The power from the traditional grid is primarily generated by fossil fuel stations, and has a highly negative environmental effect. This work does not impose an upper limit on how much power could be bought from the grid. To restrict the green cellular network from buying substantial amounts of grid power, we also impose an additional "environmental" cost of buying fossil fuel-based power [62].

$$[eqn]$$

where  $[eqn]$  and  $[eqn]$  denote the pollutant emission coefficients of grid power.

2.6. Running Horizon Framework

Although this work aims to formulate an optimal power-procurement problem for a one-day operation cycle, such a short-term approach could lead to spurious results at the end of the day. For example, the resulting optimal policy might have operators selling all stored energy at the end of the day to generate revenue. This scenario is unrealistic because the operator would prefer to keep some charge in the battery for the next day. This work devises a running horizon framework to avoid producing such impractical optimal policies and still work in a short-term finite horizon. The framework consists of two facets. First, we add fictitious revenue from storing battery charge at the end of the day. More specifically, we introduce a third revenue source, as follows:

$$[eqn]$$

where  $[eqn]$  represents the fictitious cost per unit of stored battery charge. For example, publicly available deterministic forecasts of normalized renewable power and day-ahead energy price forecasts can be applied to determine a suitable value of  $[eqn]$ .

Next, we formulate the optimal control problem for a two-day operation cycle (instead of one) but apply the optimal power-procurement policy only for one day. This approach requires two-day-ahead forecasts of renewable power, energy prices, and cellular user demand. The optimal power-procurement policy is updated daily, as and when the operators receive the two-day-ahead forecasts of these quantities. Regarding the problem formulation, the only change is that the time horizon is denoted by  $[eqn]$ , where 0 denotes beginning of the day, and T denotes the end of the second day (instead of at the end of the current day). Henceforth, this notation is applied while formulating the optimization problem.

2.7. Model Summary

The state and control vectors are denoted by  $[eqn]$  and  $[eqn]$ , respectively:

$$[eqn]$$ $$[eqn]$$

The state vector  $[eqn]$  is defined on a probability space  $[eqn]$  and each of its components satisfy the controlled dynamics in (3), (9), and (10), respectively, in  $[eqn]$ . The initial conditions are also random. The control  $[eqn]$  is an  $[eqn]$ -adapted Markovian control that is measurable, where  $[eqn]$  represents the filtration on the probability space.

Constraints 1 (A.s. instantaneous control constraints).  The controls must satisfy the following constraints for all  $[eqn]$ :

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$

The controls must also follow the following state-dependent constraint:

$$[eqn]$$

Furthermore, for a given deterministic forecast of the number of cellular users connected to the network,  $[eqn]$ , the controls are constrained to satisfy the following:

$$[eqn]$$

Constraints 2 (Probabilistic instantaneous control constraints).  The controls must satisfy the following probabilistic constraint for all  $[eqn]$ :

$$[eqn]$$

The shorthand notation  $[eqn]$  is used here. The constraints in (18) to (24) define the set of admissible Markov controls  $[eqn]$  at each time  $[eqn]$ . Then, we define the set of admissible policies as follows:

$$[eqn]$$

3. Stochastic Optimal Control Formulation

Problem 1 (Primal problem with the probabilistic constraint). Given the initial data  $[eqn]$ , deterministic forecasts for daily cellular user traffic profile  $[eqn]$ , energy spot prices  $[eqn]$  and  $[eqn]$ , renewable energy forecast  $[eqn]$ , and mobile network usage price  $[eqn]$  for  $[eqn]$ , we solve the following:

$$[eqn]$$

where the minimization is done over all controls satisfying Constraints 1 and 2 with dynamics (3), (9), and (10). We define  $[eqn]$  as the net utility of the base station. In addition,  $[eqn]$  is a Pareto parameter to weigh the two objectives: financial cost  $[eqn]$  and environmental cost ( $[eqn]$ ), where w can be chosen by the base station operators, depending on how environmentally conscious they are. The dependence of the costs on the controls  $[eqn]$  in (26) is suppressed for ease of notation.

Remark 3 (Discrete vs. continuous-time formulation).  The optimal solution to a discrete-time formulation leads to a discrete set of control decisions at each time step, dependent on the selected discretization. Applying the dynamic programming principle to a discrete-time formulation leads to a stochastic optimization problem at each time step. Finer grids substantially increase the computational work, making the problem computationally infeasible without guaranteeing that the optimal policies would converge. In a continuous-time formulation, the limit of the dynamic programming principle exists as the time step becomes infinitesimally small, leading to the derivation of the associated HJB PDE. Numerical discretization methods are applied just once to solve this PDE and obtain a solution at all time points. This approach facilitates discretization error control, ensuring convergence to the true optimal solution, and enables the flexible application of adaptive and higher-order numerical approximation schemes.

The probabilistic constraint in 2 does not allow classical dynamic programming to solve Problem 1 because the constraint controls the joint distribution of the Markovian dynamics, rather than each realization. The standard dynamic programming lemma [38] cannot be applied if the Markovian control function depended additionally on the state’s distribution. We overcome this problem using Langrangian relaxation of the probabilistic constraint for all  $[eqn]$ , yielding a continuous-time Lagrangian relaxation. For this approach, we define the time-continuous deterministic Lagrange multiplier function  $[eqn]$ , associated with the relaxed constraint in 2. We let  $[eqn]$  denote the admissible space of Lagrange multiplier functions. The new set of admissible controls is only defined by the constraints in (18)–(23) and is denoted by  $[eqn]$  at each time  $[eqn]$ . Then, we define the new set of admissible policies as follows:

$$[eqn]$$

We define the Lagrangian using terminology from the classical Lagrangian relaxation technique for constrained optimization problems [63] as follows:

$$[eqn]$$

Problem 2 (Relaxed problem). Given the initial data  $[eqn]$ , Lagrange multiplier function  $[eqn]$ , deterministic forecasts for daily cellular user traffic profile  $[eqn]$ , energy spot prices  $[eqn]$  and  $[eqn]$ , renewable energy forecast  $[eqn]$ , and mobile network usage price  $[eqn]$  for  $[eqn]$ , we solve the following:

$$[eqn]$$

where the minimization is done over all controls satisfying the constraint in 1 with dynamics (3), (9), and (10).

With the relaxation of the probabilistic constraint in 2, the dynamic programming principle can be applied to solve Problem 2. The dual problem associated with the primal Problem 1 is formulated below.

Problem 3 (Dual problem).  We determine

$$[eqn]$$

where the dual function  $[eqn]$  is given by

$$[eqn]$$

and  $[eqn]$  solves Problem 2 for a given  $[eqn]$ .

The admissible set  $[eqn]$  of Lagrange multiplier functions depends on the problem structure. For the considered optimization problem,  $[eqn]$  denotes a positive real-valued function of time t due to the deterministic nature of the relaxed constraint (24).

3.1. HJB Equation Related to Problem 2

As a first step to solving the dual problem (Problem 3), Problem 2 is solved to obtain  $[eqn]$  for a given  $[eqn]$ . Using terminology from standard continuous-time stochastic optimal control [38,40], we first define the cost-to-go function  $[eqn]$  associated with Problem 2 as follows:

$$[eqn]$$

We define the components of  $[eqn]$ , where a, r, and  $[eqn]$  denote the variables corresponding to components of the state vector  $[eqn]$ ,  $[eqn]$ , and  $[eqn]$ , respectively. In addition,  $[eqn]$  denotes the domain of  $[eqn]$ , in this case  $[eqn]$ , because  $[eqn]$  and  $[eqn]$  are  $[eqn]$ -valued random variables and  $[eqn]$  is a  $[eqn]$ -valued random variable for all  $[eqn]$  (see Equations (3), (9) and (10)). Using (32), the dual function is defined as  $[eqn]$ , where u, for given  $[eqn]$ , solves the following second-order nonlinear HJB final value PDE.

$$[eqn]$$

where  $[eqn]$  denotes the gradient vector, and  $[eqn]$  represents the Hessian matrix of a scalar-valued function. The Hamiltonian  $[eqn]$  associated with Problem 2 in (33) is defined as follows:

$$[eqn]$$

We obtain an approximation  $[eqn]$  of the dual function  $[eqn]$  by numerically approximating the solution to (33). Section 4 discusses this approach in detail.

3.2. Finite-Dimensional Approximation of Problem 3

Problem 3 is an infinite-dimensional optimization problem to satisfy the constraint in 2 pointwise in time. We discretize  $[eqn]$  into a piecewise constant function to approximate the infinite-dimensional dual problem, allowing for practical numerical optimization. This approach implies satisfying the constraint for a finite number of subintervals in  $[eqn]$ . That is, the integral of the pointwise violation of the constraint in 2 over each subinterval must be 0. We consider the time discretization  $[eqn]$  of the time domain  $[eqn]$ . The Lagrange multiplier is approximated as follows:

$$[eqn]$$

where  $[eqn]$  denotes components of the vector of Lagrange multipliers  $[eqn]$ . For the given time discretization, we define

$$[eqn]$$

With this approximation, the finite-dimensional approximation of Problem 3 is formulated.

Problem 4 (Finite-dimensional dual problem).  For a given ℓ, we determine the following:

$$[eqn]$$

The following Lagrangian is written for a given  $[eqn]$ :

$$[eqn]$$

with

$$[eqn]$$

where  $[eqn]$  denotes the  $[eqn]$  component of an ℓ-dimensional subgradient vector  $[eqn]$  of the dual function  $[eqn]$  evaluated at  $[eqn]$ . To solve Problem 4, the optimal ℓ-dimensional vector  $[eqn]$  that maximizes the dual function  $[eqn]$  is selected. Hence, Problem 4 is a convex, ℓ-dimensional optimization problem but is not generally smooth. Therefore, we apply subgradient methods to solve Problem 4. Moreover, the choice of ℓ and the corresponding time discretization grid is set to control the pointwise violation in Constraint 2. Section 4 details this approach.

4. Numerical Approach

4.1. Numerically Solving the HJB PDE

We employ the explicit Euler upwind finite-difference scheme to numerically approximate the solution to the HJB equation in (33). This method provides a convergent numerical scheme under certain conditions that can be proven for the considered problem [42].

We consider the following time and space discretization of  $[eqn]$ , where  $[eqn]$ , and the following finite-difference grid  $[eqn]$ :

$$[eqn]$$

For the dimension  $[eqn]$ , we approximate the domain  $[eqn]$  by  $[eqn]$  for an appropriately chosen  $[eqn]$ , which is scaled to  $[eqn]$ . We employ uniform grids in all dimensions; that is,

$$[eqn]$$

We combine the scaled drift and diffusion terms of dynamics (3), (9), and (10) into vector-valued functions F and G, respectively.

$$[eqn]$$

For convenience, we define the following function representing the running cost in (34).

$$[eqn]$$

For a function f, we define  $[eqn]$  and  $[eqn]$ . For a function f, we define  $[eqn]$ . Let  $[eqn]$  denote the numerical approximation of the cost-to-go function  $[eqn]$  at the point  $[eqn]$  on the grid  $[eqn]$ . Then, the explicit update rule for the upwind finite-difference scheme is given as follows:

$$[eqn]$$ $$[eqn]$$

In (41), we apply the notation  $[eqn]$  and  $[eqn]$  to denote the  $[eqn]$  component of the vector-valued functions F and G, respectively. Optimal controls at grid point  $[eqn]$ , denoted by  $[eqn]$ , are computed using the known derivatives (up to the second order) of the numerical cost-to-go-function as follows:

$$[eqn]$$

where  $[eqn]$  and  $[eqn]$  represent the upwind derivatives given by

$$[eqn]$$

In general, boundary conditions must be imposed on the HJB PDE (33) to obtain a unique numerical solution. However, boundary conditions at the boundaries  $[eqn]$ ,  $[eqn]$ , and  $[eqn]$  are unnecessary in this setting because the drift terms  $[eqn]$  naturally have a sign at the boundaries ensuring the use of an interior point to approximate first-order derivatives. Similarly, the diffusion terms  $[eqn]$  are naturally 0 at these boundaries, ensuring no exterior points are necessary to approximate second-order derivatives. At the boundary  $[eqn]$ , we impose a nonreflective boundary condition to indicate that the actual domain is not bounded at  $[eqn]$ . We impose the Courant–Friedrichs–Lewy (CFL) condition at each grid point in  $[eqn]$  to ensure the stability of the numerical scheme

$$[eqn]$$

Although the HJB PDE (33) is nonlinear, the numerical scheme given by Equations (41)–(43) is explicit because the cost-to-go function at time  $[eqn]$  is completely determined by its value and derivatives at time  $[eqn]$ . However, this relationship is nonlinear and cannot be vectorized. Moreover, the four-dimensional (4D) constrained optimization problem in (43) must be solved at every grid point in  $[eqn]$ . Depending on the structure of  $[eqn]$ , (43) can be solved analytically or numerically. Linear interpolation is applied to extend the numerical solution  $[eqn]$  from the grid  $[eqn]$  to the entire domain  $[eqn]$ . Next, the dual function  $[eqn]$  is approximated by  $[eqn]$  for a given  $[eqn]$ . Algorithm A1 in Appendix A presents a pseudo-algorithm of the numerical solver.

Remark 4 (Discretization error of the PDE solver).  The upwind finite-difference scheme in Section 4.1 is a first-order explicit scheme with a discretization error of  $[eqn]$ . When solving the overall optimization problem up to a relative tolerance of  $[eqn]$  with respect to the dual function value  $[eqn]$ , the PDE discretization error must be  $[eqn]$ . Thus,  $[eqn]$  must be  $[eqn]$ , implying the computational work of  $[eqn]$  to solve the HJB PDE once.

4.2. Estimating Subgradient DΘ

Upon solving the HJB PDE (33) numerically and obtaining an approximate cost-to-go function in the domain  $[eqn]$ , we can numerically simulate optimally controlled paths of the state variable  $[eqn]$  forward in time, from the initial condition  $[eqn]$  using dynamics (3), (9), and (10). We use the Euler–Maruyama time discretization of the SDEs in (3), (9), and (10). We consider the discretization  $[eqn]$  of the time domain  $[eqn]$  with  $[eqn]$  uniform time steps. It follows that  $[eqn]$ , and  $[eqn]$ . If  $[eqn]$  is the time-discretized version of the stochastic process  $[eqn]$ , then the Euler–Maruyama time discretization of the SDEs is expressed as follows:

$$[eqn]$$

where  $[eqn]$  for  $[eqn]$ , denotes independent and identically distributed (i.i.d.) three-dimensional standard normal random variables, and  $[eqn]$  denotes elementwise multiplication of two vectors. The optimal controls  $[eqn]$  in (46) are obtained by minimizing (43) at the current point of the state  $[eqn]$ , which may be a point outside the grid  $[eqn]$ . The values of the numerical derivatives  $[eqn]$  of the cost-to-go function at those points are linearly interpolated from the computed derivative values at the nearest points in grid  $[eqn]$ . The subgradient of  $[eqn]$  at  $[eqn]$  is computed according to (39), approximating the expected value using Monte Carlo sampling with  $[eqn]$  i.i.d. sample paths of the optimally controlled process  $[eqn]$ . Moreover, the integral in (39) is approximated using a forward Euler approximation. For each  $[eqn]$ , we consider the discretization  $[eqn]$  of the time domain  $[eqn]$ , corresponding to the  $[eqn]$  component of the subgradient, with  $[eqn]$  uniform time steps. It follows that  $[eqn]$ , and  $[eqn]$ . The numerical approximation of the  $[eqn]$  component of the ℓ-dimensional subgradient is denoted by  $[eqn]$  and given by

$$[eqn]$$

where  $[eqn]$  denotes the  $[eqn]$  i.i.d. realization of the random variables required to generate sample paths of  $[eqn]$ . The time discretization for subgradient computation and the Lagrange multiplier function  $[eqn]$  generally need not coincide. In this case, Brownian bridge interpolation is necessary to evaluate the optimally controlled paths  $[eqn]$  at the time discretization points  $[eqn]$  for each  $[eqn]$ . The same  $[eqn]$  realisations of the optimally controlled path  $[eqn]$  are used to estimate all ℓ components of  $[eqn]$ . Algorithm A2 in Appendix A presents the corresponding pseudo-algorithm.

Moreover, (47) is a Monte Carlo approximation of (39). It is subject to noisy statistical error controlled in a probabilistic sense by the choice of  $[eqn]$ . Thus, we only have access to noisy subgradients of the dual function. Hence, we apply the stochastic subgradient method (SSM) to solve the convex, nonsmooth optimization Problem 4.

Remark 5 (Subgradient estimation error).  The error in the subgradient estimation consists of two parts: a statistical Monte Carlo error from approximating the expected value in (39) using a sample average with  $[eqn]$  samples and a discretization error from approximating the integral in the expectation in (39) using a forward Euler summation with  $[eqn]$  steps. From the central limit theorem [64], the Monte Carlo approximation error is  $[eqn]$ . The discretization error for the first-order forward Euler summation is  $[eqn]$ , yielding a total error in the subgradient estimation of  $[eqn]$ . When solving the overall optimization problem up to a relative tolerance of  $[eqn]$  with respect to the subgradient, then we must estimate the subgradients up to an absolute tolerance of  $[eqn]$ . This requirement implies that  $[eqn]$  and  $[eqn]$ , yielding a total computational work of  $[eqn]$  for each subgradient estimation.

4.3. Dual Problem Optimization

After approximating the dual function  $[eqn]$  and its subgradient  $[eqn]$  for a given Lagrange multiplier function  $[eqn]$ , we must solve the nonsmooth convex dual optimization problem in Problem 4. We solve this in two stages: (i) by constructing an approximation  $[eqn]$ , of  $[eqn]$  (as in (35)), and (ii) by numerically optimizing the amplitudes  $[eqn]$  of the approximated function  $[eqn]$ . Algorithm 1 illustrates the overall procedure. Figure A1 depicts Algorithm 1 visually for further clarity.

Algorithm 1: Numerical dual optimization procedure Input:  $[eqn]$ ,  $[eqn]$ , max-iter,  $[eqn]$ ,  $[eqn]$ ,  $[eqn]$ ,             $[eqn]$ Output: Optimal controls  $[eqn]$ , optimal Lagrange multiplier function  $[eqn]$ Construct  $[eqn]$  with  $[eqn]$  using (35); Obtain  $[eqn]$  using initialization Algorithm A3 with inputs  $[eqn]$ , $[eqn]$ ,    $[eqn]$ ; Construct  $[eqn]$  with  $[eqn]$  using (35); Obtain  $[eqn]$  with the LMBM routine [65] with starting point  $[eqn]$ , number of iterations    $[eqn]$  and parameters specified in Table A3; Construct  $[eqn]$  with  $[eqn]$  using (35);

4.3.1. Lagrange Multiplier Refinement

This work constructs a finite-dimensional piecewise constant approximation,  $[eqn]$ , of the Lagrange multiplier function,  $[eqn]$ , as detailed in Section 3.2. Increasing the refinement level ℓ enhances the approximation quality of  $[eqn]$  by  $[eqn]$ , but also increases the optimization dimension (also ℓ). This work determines a sufficient ℓ that controls the violation of the relaxed constraint. We start with  $[eqn]$  and keep uniformly doubling it, until the violation of Constraint 2 is sufficiently controlled pointwise in time.

4.3.2. Numerical Optimization

This work solves the optimization of  $[eqn]$  for  $[eqn]$  with respect to the amplitudes  $[eqn]$  as discussed in Section 3.2. For a given evaluation point  $[eqn]$ , we have access to the dual function value (approximated in Section 4.1) and a noisy subgradient (approximated in (47)). Hence, the SSM [47] is employed to solve the nonsmooth optimization problem. The SSM is a subgradient-based search algorithm using a dual function evaluation and noisy evaluation of an arbitrary subgradient at each evaluation point. At each refinement level ℓ, the SSM runs until the subgradient norm is below a prescribed relative tolerance  $[eqn]$  or the number of iterations exceeds a threshold max-iter. The algorithm output at each ℓ is the optimal amplitude vector  $[eqn]$ . However, this optimal value is not necessarily reached at the final iteration of the SSM because the objective function to be maximized does not increase at every step of the subgradient method [47]. Hence, we track the highest value of the dual function and store its corresponding  $[eqn]$  as the optimal value.

The choices for the starting point and step-size are crucial for quick convergence of the SSM, especially in higher dimensions. This work uses a nonsummable diminishing step-size of  $[eqn]$ , where k denotes the iteration number. The convergence of the SSM with step-size is proven in [47]. The associated constant, denoted by  $[eqn]$  in Algorithm 1, is tuned once such that the step-size taken is  $[eqn]$ .  $[eqn]$  need not be tuned again for a different set of parameters with the current problem structure. For enhanced performance, the step-size at each iteration is adjusted by dividing the current subgradient by its norm. Furthermore, the SSM runs for the first  $[eqn]$  iterations using a constant step-size before reverting to a diminishing step size. The Lagrange multiplier function  $[eqn]$  constructed using the computed optimal amplitudes  $[eqn]$  is used as the starting point of the SSM at level ℓ. However, we must still choose a good starting point at level  $[eqn]$ . This work devises an initialization algorithm and complements it with a nonsmooth deterministic optimizer.

Remark 6 (SSM convergence rate).  A convergence rate of  $[eqn]$  for the SSM has been proven for a class of convex functions with suitable step-size choice [66]. When solving the optimization problem up to a relative tolerance of  $[eqn]$  with respect to the dual function value  $[eqn]$ , the SSM requires an expected number of  $[eqn]$  iterations to converge under the assumptions in [66].

Remark 7 (Primal feasible solution).  Note that Algorithm 1 uses the subgradient as a stopping criterion, instead of a duality gap. This is because the output of Algorithm 1 is already a primal feasible solution up to given relative tolerance  $[eqn]$ . That is, the optimal controls produce a solution that minimizes the primal cost while only violating Constraint 2 lesser than a small value of  $[eqn]$ .

4.3.3. LMBM-Boosted Initialization

Level  $[eqn]$  implies that the Lagrange multiplier function  $[eqn]$  is constant in time, with the constant denoted by  $[eqn]$  (see (35)). First, this work devises an initialization algorithm starting with the arbitrary point  $[eqn]$  and estimating the dual function value and a subgradient at that point. Then, if the subgradient at  $[eqn]$  is positive,  $[eqn]$  is continually increased by a factor  $[eqn]$  until the subgradient becomes negative or its norm reaches a prescribed relative tolerance  $[eqn]$ . Conversely, if the subgradient at  $[eqn]$  is negative,  $[eqn]$  is continually decreased by a factor  $[eqn]$ , until the subgradient becomes positive or its norm reaches relative tolerance  $[eqn]$ . The final obtained point is stored as  $[eqn]$ . Algorithm A3 in Appendix A presents the pseudo-algorithm of the devised initialization procedure. Figure A2 depicts Algorithm A3 visually for further clarity.

To obtain a better starting point,  $[eqn]$  is used as a starting point for a deterministic nonsmooth optimization routine called the LMBM [45,56]. Bundle methods are more robust than a simple subgradient-based method because they approximate the entire subdifferential of the objective function, enhancing its convergence speed [44]. However, the LMBM requires access to deterministic estimates of the objective function value and its subgradient at each evaluation point. Although this work only has access to noisy Monte Carlo estimates of subgradients (47), we can still employ the LMBM routine to reach a good starting point for the SSM, especially in the regime where changes in the objective function value considerably outweigh the noise in subgradient estimates. This work runs the LMBM for a fixed number of iterations  $[eqn]$ . Table A3 in Appendix B.4 specifies the parameters required to run the LMBM routine [65]. The output of this routine  $[eqn]$  is applied as a starting point for the SSM algorithm described in Section 4.3.2.

Remark 8 (*Scale of  $[eqn]$ *). With respect to the parameter  $[eqn]$  that determines the confidence level of satisfying Constraint 2, the optimal primal cost  $[eqn]$  associated with Problem 4 is an $[eqn]$ term, while the optimal subgradient would be  $[eqn]$  (ideally 0). Equation (38) then implies that the optimal Lagrange multiplier  $[eqn]$  is  $[eqn]$ . This yields a good guess for an initial point  $[eqn]$  for the LMBM-boosted initialization procedure.

5. Numerical Experiments and Results

This section presents a model example of a cellular base station powered by the German power grid. This section describes the system and all dynamics driving the operation of the base station, and provides the results of applying the proposed numerical approach to solve the optimal power-procurement problem for the system.

5.1. Description of Model Cellular Base Station System

Figure 1 schematically illustrates the considered base station model. Table A1 provides the descriptions and numerical values of all coefficients used to describe the model.

The daily mobile user traffic profile primarily drives the power demand of a cellular base station. The description of the daily traffic profile consists of two facets: (i) the number of people connected to the network  $[eqn]$ , and (ii) the physical distribution of users around the base station  $[eqn]$ . The numerical experiment applies the following sinusoidal profile to model  $[eqn]$ :

$$[eqn]$$

This model and the stochastic versions of it have been empirically demonstrated to approximate practical user patterns closely for calibrated values of the smoothness parameter  $[eqn]$  [67,68]. The parameter  $[eqn]$  determines the rate of increase or decrease in  $[eqn]$  during the day. A higher value of  $[eqn]$  indicates a steeper increase or decrease rate, and  $[eqn]$  in the numerical experiment. This model also considers two peak hours (about 09.00 to 10.00 and 20.00 to 21.00) and few connected users at night (00.00 to 05.00). Moreover,  $[eqn]$  and  $[eqn]$  denote the minimum and maximum number of users, respectively, that are connected to the base station at any time. This work sets  $[eqn]$  and  $[eqn]$ . These values were selected using the data on the average population densities of German cities [69] and the geographical area served by an average cellular base station in Germany [70]. Figure 2 depicts the model for  $[eqn]$ .

The user distribution  $[eqn]$  crucially determines the estimation of the outage proportion  $[eqn]$  at given time t (Section 2.1). Some simple 2D distributions have been employed to model the mobile user distribution in urban zones efficiently. For example, the uniform distribution can be employed to model the mobile user distribution in suburban areas very well, and the Gaussian distribution is a good model for mobile user distributions in industrial zones [71]. Figure 3 visualizes this model. Analytical expressions can be derived for the corresponding  $[eqn]$  for these distributions. Appendix B.2 provides the closed-form expressions. For the numerical experiments, we apply the symmetric 2D Gaussian distribution centered on the site of the base station  $[eqn]$  with a diagonal covariance matrix with variance  $[eqn]$  along both directions.

We apply the SDE in (3) to model the instantaneous Nakagami fading channel, whose invariant distribution is a shifted-gamma distribution with the shape parameter  $[eqn]$  and scale parameter  $[eqn]$ . We set  $[eqn]$  and  $[eqn]$  so that a quick, smooth convergence occurs to the invariant distribution and  $[eqn]$  is bounded for all  $[eqn]$ . We set the shift  $[eqn]$ , corresponding to the value of  $[eqn]$  that can still achieve an outage proportion of  $[eqn]$  in the worst-case scenario when  $[eqn]$ . Moreover, we set the initial distribution  $[eqn]$  the same as the invariant distribution of process  $[eqn]$ . Figure 4 visualizes the samples paths and the invariant distribution of the SDE in (3).

We set  $[eqn]$  EUR/h per person for  $[eqn]$ . The price was taken from the Vodafone 5G prices in Germany [72]. The constant price assumption is generally valid because regulators set these prices and do not move within a day.

We apply the SDE in (9) to model the instantaneous normalized generated renewable power. This work uses the 2023 German wind power forecast and production data with a 15-min frequency from the operator 50Hertz [73] to calibrate the SDE in (9). The forecast data are employed to construct  $[eqn]$  and  $[eqn]$ , and the discrepancy between the forecast and true production data is applied to calibrate the parameters  $[eqn]$  and  $[eqn]$  in (9). Apart from these parameters, we also calibrate an additional parameter  $[eqn]$ . The discrepancy between the forecast and production at time  $[eqn]$  is not typically zero in real-world applications.  $[eqn]$  is defined as the time before  $[eqn]$ , for which the forecast error can be set to zero. This approach ensures that uncertainty also exists in renewable power production at  $[eqn]$ . We follow the calibration procedure detailed in [30] (see [30] for details). The mean-reversion parameter  $[eqn]$  in (9) is defined as follows:

$$[eqn]$$

The calibrated values for the above data are  $[eqn]$ ,  $[eqn]$ , and  $[eqn]$ . Given the deterministic forecasts of wind power from the same operator in 2024, we can generate sample paths of the SDE in (9) with parameters  $[eqn]$  calibrated from 2023 data. Figure 5 illustrates the results for two sample days in 2024.

This work applies (10) to model the instantaneous normalized battery charge. As given in (10), the maximum charging capacity  $[eqn]$  and the maximum discharge capacity  $[eqn]$  are defined as functions of the stored battery charge. This function is typically called the characteristic curve of the battery. The numerical experiment employs a simple characteristic curve illustrated in Figure 6. The characteristic curve imposes (11), implying that the battery cannot discharge power when it has zero charge and can no longer charge itself when fully charged.

Forecasts of buying and selling prices,  $[eqn]$  and  $[eqn]$ , are constructed from publicly available day-ahead spot prices of grid power. We set  $[eqn]$ , implying that the buying and selling prices of power are the same. This assumption is valid because the market commonly ensures that this holds in real life. The numerical experiment employs the German day-ahead spot-price data from 2024 available in [74]. However, these data are generally available only at a 1-h frequency. Hence, we smoothed and extended the data to domain  $[eqn]$  by fitting a polynomial curve to them and using the polynomial function as the input  $[eqn]$ . Figure 7 illustrates an example of the spot-price data.

5.2. Numerical Results

We present the results of solving the stochastic optimal control problem for the model cellular base station (Section 5.1) using the numerical approach in Section 4.

5.2.1. Cost-to-Go Function of the Dual Problem

The numerical solution of the HJB PDE in (33) for a given Lagrange multiplier function  $[eqn]$  (35) is obtained using the scheme in Section 4.1. We apply a uniform rectangular spatial grid with  $[eqn]$ , where  $[eqn]$  is set as the lowest value that satisfies the CFL condition (45). In this case, we obtain  $[eqn]$ . The presence of second-order derivatives in the HJB PDE in (33) yields  $[eqn]$  for a numerically stable solution, explaining the high value obtained for  $[eqn]$ . These discretization parameter choices ensure that the relative numerical PDE discretization error is below  $[eqn]$ . We set  $[eqn]$  as the upper  $[eqn]$  quantile of the invariant distribution of process  $[eqn]$ .

The optimization problem in (43) must be solved at each point in the grid  $[eqn]$ . For the simple cases of uniform or Gaussian cellular user distribution, we obtain closed-form expressions for  $[eqn]$  (Appendix B.2). The resulting  $[eqn]$  ensures that (43) is a convex optimization problem and can be solved analytically. Appendix B.3 presents additional details. Figure 8 presents three slices of the numerical cost-to-go function  $[eqn]$  with the optimal Lagrange multiplier function  $[eqn]$ , each with two fixed state variables. The numerical dual function approximation  $[eqn]$  is obtained by evaluating  $[eqn]$ .

5.2.2. Obtaining Reference Costs for Comparison

This work compares the values for the dual cost above with reference values obtained from two related, easy-to-solve stochastic optimal control problems. This comparison was conducted to ensure that the output of the numerical algorithm makes sense with respect to these reference values. First, we solve the same optimization problem as in Problem 1, but with the probabilistic constraint replaced with a deterministic a.s. constraint, solving the following problem.

Problem 5 (Primal problem with the a.s. QoS constraint). Given the initial data  $[eqn]$ , and deterministic forecasts for the daily cellular user traffic profile  $[eqn]$ , energy spot prices  $[eqn]$ , renewable energy forecast  $[eqn]$ , and mobile network usage price  $[eqn]$  for  $[eqn]$ , we solve the following:

$$[eqn]$$

where the minimization is subject to dynamics (3), (9), and (10), and is done over all controls that satisfy the constraint in 1 and the following a.s. constraint (instead of the constraint in 2).,

$$[eqn]$$

where (50) provides a solution for the control  $[eqn]$  for all  $[eqn]$ .

Hence, the control vector for this problem has only three components,  $[eqn]$ . Problem 5 is a standard continuous-time stochastic optimal control problem; hence the dynamic programming principle can be used to derive and numerically solve the associated HJB PDE, analogous to Section 3.1. The cost-to-go function for this problem is expected to be higher than the dual cost of Problem 4 because the network operator would need to buy more power to keep  $[eqn]$  constant at all times, whereas the operator is allowed to relax the constraint probabilistically in Problem 4 as described in Constraint 2. Next, we solve the same optimization problem as Problem 1, but without Constraint 2.

Problem 6 (Primal problem with no QoS constraint). Given the initial data  $[eqn]$ , and deterministic forecasts for the daily cellular user traffic profile  $[eqn]$ , energy spot prices  $[eqn]$ , renewable energy forecast  $[eqn]$ , and mobile network usage price  $[eqn]$  for  $[eqn]$ , we solve the following:

$$[eqn]$$

where the minimization is taken over all controls satisfying the constraint in 1 with dynamics (3), (9), and (10).

Problem 6 is also a standard continuous-time stochastic optimal control problem; hence the dynamic programming principle can be applied to derive and numerically solve the associated HJB PDE. Without the Lagrangian relaxation, the cost-to-go function associated with this problem is the same as in (32). Thus, this problem can be solved using the numerical scheme in Section 4.1, setting  $[eqn]$ . The cost-to-go function for this problem is expected to be less than the dual cost of Problem 4 because the network operator only maximizes revenue without ensuring that QoS is achieved. The HJB PDEs corresponding to Problems 5 and 6 must be numerically solved only once because no associated dual problem needs to be solved.

5.2.3. Dual Optimization Results

This work runs Algorithm 1 with parameters set to values in Appendix B.4. Notably, we set  $[eqn]$  and  $[eqn]$ , implying that the network operator must ensure that the proportion of users in outage  $[eqn]$  should be less than  $[eqn]$  with a probability of  $[eqn]$  for all  $[eqn]$ . Figure 9 presents the results of the initialization Algorithm A3 and the LMBM routine. Figure 9a reveals the dual cost versus  $[eqn]$  at each iteration of Algorithm A3. The dual cost increases with the number of iterations, and it is within the two reference values attained by solving Problems 5 and 6. Figure 9b plots the corresponding estimated subgradient values. The initialization algorithm reaches the prescribed tolerance in three iterations. Figure 9 reveals that the LMBM offers a more optimal point than the output of Algorithm A3, providing an excellent starting point for the dual optimization in Algorithm 1.

Figure 10 illustrates the evolution of the dual function values and its corresponding subgradient at each iteration of the SSM for a fixed ℓ. Figure 10a verifies that not every step of the SSM is taken towards maximization. This behavior illustrates the difficulty of solving the nonsmooth optimization in Problem 4. Due to the satisfactory choice of starting values and step sizes (Section 4.3), Algorithm 1 converges for each ℓ within max-iter iterations with respect to the subgradient.

Figure 11 plots the evolving optimal dual function  $[eqn]$  and its corresponding subgradient  $[eqn]$  constructed using the optimal amplitudes  $[eqn]$  at each level ℓ of Algorithm 1. Table 1 lists the optimal dual cost at each stage of Algorithm 1. Table 1 presents minor increases in the dual function evaluations as the refinement level ℓ increases, implying that we successively obtain more optimal solutions as the Lagrange multiplier function  $[eqn]$  becomes finer. We normalize the subgradient vector norm by the number of dimensions for a fair comparison between the norm of vectors of varying dimensions.

Table 1 also highlights the advantage of applying a probabilistic constraint to achieve the QoS in cellular networks, instead of doing it a.s. for all  $[eqn]$  (as in Problem 5). The optimal dual cost achieved by Algorithm 1 is much closer to the expected optimal cost of the unconstrained Problem 6 and around four times less expensive than satisfying (50), a.s. (Problem 5). In rare occasions when an unreasonable surge in demand ( $[eqn]$ ) occurs, or a very low wireless fading channel ( $[eqn]$ ) exists due to environmental conditions or low availability of renewable power ( $[eqn]$ ), the numerical approach allows the network operator to relax the stringent QoS constraint to save substantial expenditure while still serving as many customers as possible. This approach contrasts the solution to Problem 5, which enforces (50) even in such periods, significantly increasing operating expenditure.

Figure 12 presents the evolution of the optimal Lagrange multiplier function as the refinement level ℓ increases in Algorithm 1. For the given example, Algorithm 1 reaches prescribed  $[eqn]$  at  $[eqn]$ , and does not refine the Lagrange multiplier function further. High-demand periods (Figure 2) correspond to high Lagrange multiplier function values. This behaviour is expected because a high  $[eqn]$  forces the mobile operator to satisfy Constraint 2 in times of high demand, although it significantly increases costs. The oscillation frequency of the base station power demand dictates the oscillations in the Lagrange multiplier function.

5.2.4. Optimal Power-Procurement Policy

With the optimal Lagrange multiplier function  $[eqn]$  from Algorithm 1, the dual cost is EUR  $[eqn]$  (Table 1). Figure 13 illustrates the optimal power-procurement policy for the model cellular base station for the next 24 h using the optimal solution, including the (a) net power consumption of the base station  $[eqn]$  (refer (1)), (b) power drawn from the battery to run the base station  $[eqn]$ , (c) power bought from the grid  $[eqn]$ , (d) power sold back to the grid  $[eqn]$ , and (e) optimal charge stored in battery  $[eqn]$ . Figure 13a indicates that the base station consumes high power during periods of high demand (refer Figure 2). Figure 13b–d reveal that the battery provides most of the power to run the base station throughout the day. The clean, renewable energy stored in the battery powers a large proportion of base station operations so that the mobile operator does not need to buy expensive grid power. This approach reduces costs and the carbon footprint of the base station. The operator must only buy some power from the grid to satisfy high demand during peak hours (Figure 13c). The base station sells a negligible amount of power back to the grid (Figure 13d). The optimal policy also advises the operator to store incoming renewable power and charge the battery during low-demand hours before using that power to serve customers during high-demand hours (Figure 13e). Figure 13e also indicates that the battery has some stored energy at the end of the day that could be used the next day. This observation establishes the effect of implementing the running horizon framework (see Section 2.6), although we solve the problem for a short-term horizon.

We estimate the time-pointwise violation of Constraint 2 using the Monte Carlo method and compare it with the final subgradient vector  $[eqn]$  from Algorithm 1 to validate the optimality of the solution further. We run approximate sample numerical paths of the optimally controlled dynamics, according to (46), with the optimal policy above. Figure 14b plots this information for a discrete set of points in  $[eqn]$ . We also compare this with the optimal policy attained from solving the unconstrained problem in Problem 6 in Figure 14a. Figure 14a reveals that the user outage proportion  $[eqn]$  for the optimal policy with no constraints on the QoS is almost consistently greater than the threshold  $[eqn]$ . In contrast, the proposed optimal solution controls the probabilistic constraint in an integral (averaged) sense (refer (39)), even though there are some oscillations around zero in the time-pointwise violation of Constraint 2. Figure 14b illustrates this behavior, with the time-averaged violation in the probabilistic Constraint 2 close to zero up to the prescribed relative tolerance  $[eqn]$ .

5.2.5. Sensitivity Analysis

This work investigates the performance of the proposed algorithm under certain scenarios that could be of interest to network operators, which are listed below.

Scenario A: No incoming renewable power for a day ( $[eqn]$  for all  $[eqn]$ ).Scenario B: The wireless fading channel is substantially low due to extreme weather ( $[eqn]$  is halved).Scenario C: High weight is assigned to minimizing operating expenditure ( $[eqn]$ ).Scenario D: The price users pay to connect to the network substantially reduces ( $[eqn]$  EUR/h per person).Scenario E: The probabilistic QoS constraint 2 must be satisfied with low confidence ( $[eqn]$ ).

All parameters and dynamics are set as in Table A1 and Section 5.1, except for the prescribed parameter change to simulate each scenario. For each scenario, we solve three problems: (i) dual optimization (Problem 4), (ii) Problem 5 with a.s. QoS constraint, and (iii) Problem 6 with no QoS constraint. Table 2 reports the optimal costs in each case for all five scenarios.

In all scenarios, the optimal dual cost lies between the primal costs of the corresponding unconstrained and a.s. QoS constrained problems. Moreover, the optimal dual cost is closer to the primal cost of the unconstrained problem in all cases because it is costly to satisfy the QoS constraint a.s. at all times. The optimal dual cost in Scenario A is 33% higher than the optimal dual cost in Table 1 because there is no “free” incoming renewable energy for the operator to run the base station, increasing the expenditure. The optimal dual cost of Scenario C is 102% higher than the optimal dual cost in Table 1, directly implying the high weight assigned to financial profits, with no regard for its carbon footprint. The optimal cost of Scenario E confirms the intuition that it is less expensive to satisfy the QoS with low confidence (80%), compared to 90% confidence. To quantify differences in the optimal procurement policy for each scenario, we produced plots analogous to Figure 13, and approximated the area below the mean curve for each plot to obtain the expected energy balance. Table 3 lists these values. The base scenario in Table 3 refers to the original problem described in Section 5.1.

In Scenario A, the operator buys six times more energy than the base scenario. The lack of incoming renewable power forces the network operator to buy more power from the grid to satisfy the QoS constraint. In Scenario C, the operator buys almost half the energy consumed by the base station from the grid and sells a considerable amount of energy stored in the battery. The optimal policy is to make as much financial profit as possible because a high weight is assigned to minimizing financial costs. The base station consumes 9% more energy in Scenario B than the base scenario because the energy demand of satisfying QoS during periods of low fading channel (extreme weather) is higher. An 8% decrease occurs in the energy consumed by the base station in Scenario E compared with the base scenario. Less energy is required to satisfy the QoS constraint with lower confidence.

This work also performs scenario simulations over probability distributions of both model and algorithmic parameters in Table A1 and Table A2, respectively. Table A4 lists the various probability distributions used to randomize the simulations. Algorithm 1 is run independently for 50 i.i.d. samples of the various parameters with no additional tuning between runs. Figure illustrates the optimal dual cost and subgradient achieved in each case. Figure 15b reveals that Algorithm 1 achieves the prescribed subgradient tolerance in all 50 cases, demonstrating the robustness of the proposed numerical procedure. These results emphasize the practicality of applying the proposed numerical approach as a decision-making tool for stakeholders in the wireless communications industry.

6. Conclusions

This work formulates a short-term, continuous-time stochastic optimal control problem that minimizes the operating expenditure and carbon footprint of a cellular wireless network powered by a hybrid energy system and equipped with an energy-storage capacity under uncertainty in renewable energy and wireless fading channels. We constructed an SDE with the shifted-gamma invariant distribution to model the stochastic Nakagami wireless fading channel. We employed an SDE model [30] to model the uncertain renewable power, with parameters calibrated using publicly available German energy market data from 2023. The QoS of the wireless network was posed as a probabilistic constraint to address uncertainties efficiently in the system. We developed a novel Lagrangian relaxation procedure in continuous time to relax the above probabilistic constraint. We developed a robust numerical procedure to solve the corresponding dual problem based on a finite-dimensional approximation of the Lagrange multiplier function, using an SSM algorithm driven by the numerical solution of the associated HJB PDE. An LMBM-boosted initialization procedure and an adaptive refinement of the Lagrange multiplier function enhanced algorithm performance. A running horizon framework resolved the spurious results at the end of the day due to the short-term horizon problem formulation. The practical application of the numerical procedure was demonstrated on a model cellular base station driven by the German power system and daily cellular traffic demand. The results illustrate the viability and efficiency of the proposed approach for real-world problems.

A potential future research direction is incorporating stochasticity in other dynamics, such as energy prices  $[eqn]$ , and the cellular user traffic profile  $[eqn]$  [68]. This research would build a more realistic algorithm that is robust to random fluctuations in these dynamics but would increase the dimension of the associated HJB PDE. Another limitation of the current work is that small values of  $[eqn]$  ( $[eqn]$ ) cannot be achieved, which modern-day cellular networks aim to achieve, because the naive Monte Carlo method is employed to estimate the subgradient. For a small  $[eqn]$ , the Fokker–Planck PDE for the joint distribution of the optimally controlled dynamics could be solved to get a deterministic estimate of the subgradient. Alternatively, importance sampling could be applied to reduce the relative statistical error of the Monte Carlo subgradient estimates. A further improvement is devising a better strategy for refining the Lagrange multiplier function  $[eqn]$ . The computed time-pointwise violation in constraint values could be employed to construct a clever strategy to refine  $[eqn]$  instead of equally spaced refinements.