Smart home gateway system over Bluetooth low energy with wireless energy transfer capability

Abstract

As billions of sensors and smart meters connect to the Internet of Things (IoT), current wireless technologies are taking decisive steps to ensure their sustainable operation. One popular IoT scenario features a smart home service gateway, which becomes the central point of user’s home environment facilitating a multitude of tasks. Given that most IoT devices connected to residential gateway are small-scale and battery-powered, the key challenge is to extend their lifetime without recharging/replacing batteries. To this end, a novel radio technology named Bluetooth low energy (BLE) has recently been completed to enable energy-efficient data transfer. Another inspiring innovation is the capability of sensors to harvest wireless energy in their local environment. In this work, we envision a scenario where many in-home sensors are communicating with a smart gateway over the BLE protocol, while at the same time harvesting RF energy transmitted from the gateway wirelessly via a dedicated radio interface. We thoroughly investigate performance limitations of such wireless energy transfer interface (WETI) with dynamic analytical model and with important practical considerations. Our methodology delivers the upper bound on WETI operation coupled with BLE-based communication, which characterizes ultimate system performance over the class of practical radio and energy resource management algorithms.

1 Introduction and background

Industry experts predict that on the order of 50 billion unattended devices be connected to the Internet by the year 2020, thus bringing revolutionary changes to how people, business, and society interact [1]. Such massive connectivity enables a plethora of emerging applications, ranging from utilities and vehicular telematics to healthcare and consumer electronics. Consequently, wireless technology has recently taken decisive steps to ensure sustainable operation of diverse Internet of Things (IoT) applications. As a response, mobile network operators are beginning to deploy novel IoT solutions allowing their customers additional degrees of control over already familiar environments.

One very attractive customer scenario is a smart home system, where residential gateways have already been used to provide broadband connectivity, quality of service (QoS), and software applications to the end users. As these systems evolve, operators foresee that smart home service gateway will be increasingly employed to control a wide variety of additional user functions, from home energy management and automation to social media interaction, connected storage, as well as multi-device printing and media streaming [2]. Hence, smart home gateway (HG) is becoming the central point of user’s home environment, and it is deemed increasingly important to ensure its flexible and stable operation.

Given that the majority of sensors, actuators, and smart meters connected to the smart HG are small-scale and battery-powered, the paramount concern in delivering sustainable home networking is energy efficiency. Indeed, the incentive to achieve all-time wireless connectivity should not, ideally, compromise the need to preserve the battery life of small-scale wireless sensors and maximize their operation time without recharging. With this constraint, intricate performance trade-offs arise between an individual device and the entire network. Ultimately, the goal of respective research is to not just extend device battery lifetime in order to provide a better user experience but also to be efficient with energy as a resource within a broader system and environmental context.

To this end, several exciting innovations have recently emerged, which become very relevant in the context of smart HG scenario. One of them is the introduction of novel energy-efficient technology for wireless sensor networks, named Bluetooth low energy (BLE). The BLE has become a result of Bluetooth Core Specification 4.0 [3] in 2010. The major purpose of developing the protocol was to enable transceivers with lower power consumption, reduced complexity, and more competitive price than the ones available with the classic Bluetooth [3, 4]. Although BLE operates in the same 2.4 GHz band and has inherited some features of the classical Bluetooth, the two technologies are not compatible [5, 6].

The initial studies (e.g., [5–7]) have shown that BLE technology can provide 1.5–2 times higher throughput, is more energy efficient, and has lower cost compared to alternative market solutions (e.g., IEEE 802.15.4). Another significant benefit of BLE, which is likely to boost its utilization even further, is its compatibility with most of mobile phones, tablets, and computers available on the market today. As the result, in the recent years, the majority of large consumer electronics manufacturers has started to replace in their product radio transceivers supporting only Bluetooth Classic with the ones compatible with both Bluetooth Classic and BLE (see e.g., [8] for the list of existing devices).

Another crucial innovation, which has captured the interest of academia and industry alike, is the energy harvesting capability of modern wireless equipment, driven by the evolution of ultra low-power electronics. With energy harvesting, devices may scavenge additional energy from the ambient environment, including solar power, electromagnetic waves, thermal energy, wind energy, salinity gradients, and kinetic energy [9]. The challenges and constraints faced by energy harvesting devices differ from those faced by devices powered by conventional energy sources and thus bring a fundamental shift in design principles compared to traditional systems with battery-operated devices [10].

Recent literature has already captured a rich set of unique energy harvesting features: the intermittent nature of harvested energy [11], the limited capacity and leakage of energy storage devices [12], as well as the constraints on device complexity [13]. Existing publications differ in their target scenarios (single, multiple, and cooperative links) [14, 15], various optimization targets and approaches (e.g., on-line vs. off-line optimization) [16, 17], key practical constraints (battery capacity and leakage, circuit power consumption, etc.) [18], as well as energy arrival considerations (continuous vs. discrete flows) [19]. Most recently, radio engineers have shown the possibility for a mobile device to use the excess availability of RF energy present in its local environment. Some call this the true future of energy harvesting technology, when a user does not need to charge the mobile device and only replaces a battery whenever it quits recharging [20].

The overview of the recent advances related to far-field wireless power transfer is presented in [21], where the authors focus on architectures of wireless energy receivers and compare the energy conversion efficiencies of the state-of-the-art hardware implementations. Moreover, [22] provides a comprehensive survey of the various energy transfer techniques and overviews a wide variety of design issues related to wireless power transfer and its practical applications. The paper discusses in detail the cases of single- and multi-hop energy transfer, presents the results, and discusses the capabilities of using multi-antenna wireless energy receivers and transmitters, as well as analyzes the perspectives of constructing RF-powered cognitive networks. The problem of optimizing the waveform for improving efficiency of wireless power transfer has been explored in [23].

It has been shown that the signals featuring high peak-to-average power ratio (PAPR) are able to improve RF-DC conversion efficiency in rectifier circuits. Several real-life RF-powered applications for different practical scenarios have recently been implemented and described in [24–26]. Particularly, in [24], the authors describe a biomedical sensor/actuator implanted in human body which is remotely powered by a sink located on the body surface. The solution provided by [25] scavenges ambient energy from AM band and buffers it for further use by structural health monitoring wireless sensor. The authors claim that their scavenger could operate having the input radio signal as weak as −39 dBm. Finally, in [26], the authors present a battery-less WSN transceiver featuring simultaneous energy and information transfer, where the transceiver receives energy and data requests in 915 MHz downlink channel and replies in 2.45 GHz uplink channel.

Inspired by the above developments, we envision that wireless energy harvesting would be increasingly employed in future smart HG deployments to prolong their operation time without recharging/replacing batteries while communicating over BLE technology. Consequently, in this work, we study a system, which comprises two major components: the energy-limited sensor nodes and the smart HG (see Fig. 1). The gateway is equipped with multiple BLE transceivers that allow it to establish several active data sessions at a time. Additionally, it features a dedicated non-BLE radio interface which could transfer energy wirelessly to the sensors and thus replenish their energy supply. To this end, we are interested in an upper bound on wireless energy harvesting capability by the sensors, which reveals the limits of any practical energy distribution mechanism. Particularly, we aim at evaluation of the key system parameters, such as overall data delay, as well as the probability of battery drain (power outage), and the average level of remaining battery power.

Fig. 1

Envisioned network architecture with a home gateway and sensors

Fundamentally, there exist two possibilities for the analysis of our system, namely, discrete and hybrid continuous/discrete approaches, which differ in the properties of energy circulation process. Discrete techniques have been widely studied in the past literature, while hybrid techniques may be preferred when dealing with small granularity (which would otherwise result in prohibitive number of discrete states) and require solving a system of differential equations. A comprehensive overview of hybrid techniques (fluid models) may be found in [27]. Further, the work in [28] offers a general overview of fluid systems and numerical approaches to analyze them, while [29] describes in detail many special case problems for more general systems. Given the peculiarities of the considered smart home deployment, we focus on offering a hybrid fluid model in the course of this paper.

In summary, the primary contributions of this paper are:

• 1. A novel concept design of a smart home system, featuring BLE-based communication and wireless energy transmission from a residential gateway to the battery-powered sensors. The proposed architecture is novel and may be attractive for future commercial deployments.

• 2. A dynamic flow-level 3D model for spatial sensor deployment and communication, which is more adequate for evaluating the data and energy transfer capabilities. Such model, when sensors are assumed to be located within a sphere, is a step ahead by contrast to past 2D models, full-buffer traffic assumptions, and symmetric systems where channel characteristics remain constant.

• 3. An oracle-assisted upper bound corresponding to the best-case energy acquisition and consumption, which allows to quantify the optimistic regions of the system operation, hence characterizing ultimate system performance over the class of practical radio and energy resource management algorithms.

• 4. A detailed technical discussion on the state-of-the-art communication and energy transfer capabilities with contemporary wireless technology. In particular, it reveals the practical limitations of wireless energy transfer, major impediments to sustainable energy operation, as well as indicates the way forward.

The rest of the text is organized as follows. Section 2 introduces the primary assumptions of the proposed system model, while Section 3 details our mathematical approach to characterize the system under study. Further, Section 4 summarizes the technology-related discussions on practical limitations of information and energy transfer, followed by some numerical results and conclusions.

2 Proposed dynamic system model

In this section, we introduce our integrated mathematical model, which captures system dynamics, spatial distribution of the sensing nodes, as well as wireless energy transfer/harvesting capabilities. The proposed system model is a powerful abstraction able to shed light on the operation of a realistic wireless network with rechargeable transmitters. We summarize the core assumptions of the model in what follows.

Our proposed model investigates the flow-level dynamics of a centralized wireless sensor network observed in stationary conditions. This network comprises a fixed number M of wireless devices performing the sensing functionality, that is, collecting temperature, motion, utility usage, or any other data, depending on the particular smart home application. The sensors are spatially distributed around the smart HG (home gateway) as shown in Fig. 1.

Assumption 1.

To characterize the averaged system performance spatially, we employ a stochastic spatial model in contrast to imposing fixed sensor locations. To this end, we assume that the locations of the sensors are driven from a Poisson point process (PPP) onto $\mathbb {R}^{3}$.

Assumption 2.

Each sensor is equipped with two radio interfaces (as per our proposal in Fig. 2 b). The first radio interface (BLE transceiver) enables the sensor to communicate its data to the gateway at the fixed transmission rate r (according to the BLE protocol, see Section 4). Secondly, each sensor node hosts a separate wireless energy transfer interface (WETI), which harvests energy of the radio signal from the gateway and stores it in the rechargeable battery.

Fig. 2

Proposed structures of smart home gateway with energy transfer (a) and wireless sensor with energy harvesting (b) capabilities

Generally, BLE communication enables both uplink (from sensors to smart HG) and downlink (reverse) data transmissions, but for the purpose of exposition, we omit the consideration of downlink direction thus focusing solely on uplink traffic. The downlink communication may easily be taken into account by means of considering the necessary BLE timings (as detailed in Section 4).

Assumption 3.

The smart HG is equipped with multiple BLE transceivers that enable the gateway to handle up to K simultaneous data links at a time (see Fig. 2 a). We further assume that the BLE transceivers at the gateway are synchronized and do not interfere with each other. The gateway additionally hosts a WETI radio that can share/transfer electromagnetic energy to the sensors and thus refill their power supply (see Section 4 for related practical considerations).

For the sake of notation, we term the sensor active, if it has some sensed data to transmit, or inactive otherwise.

Assumption 4.

Once a sensor acquires data to transmit, it enters the active state and remains there until these data are transmitted successfully (i.e., served by the gateway). We assume that an inactive sensor may generate new data of random size with the mean rate λ. The size of data is exponentially distributed with the average of μ ⁻¹.

In practice, once sensed data have arrived, the reporting sensor begins to send advertisement packets in the BLE advertisement channels (see Section 4 for detailed description) to notify the HG, which is subject to some delay. The gateway then constructs a particular schedule to serve these data by employing the available (unoccupied) data channels. The corresponding scheduling algorithm is typically implementation-specific and rather complex, attempting to ensure that the system operates with the highest possible efficiency.

To abstract away the intricate features of practical scheduling and provide good first-order understanding of the BLE operation with WETI capability, we seek to offer an optimistic estimate of the system operation. Correspondingly, we consider that the gateway is capable of constructing the optimal (in some sense) transmission schedule for the active sensors and that the use of advertising channels to do so does not impose overhead on the system. We thus arrive at an upper bound on any practical radio resource management algorithm.

Following such optimal schedule, the gateway establishes the connections with different sensors using the available channels. If there is no available channel for a newly activated sensor, it should be waiting for service according to the FIFO discipline. Once all the sensed data are forwarded to the HG, the sensor returns to inactive mode until the next activation. We further emphasize that the data transmission (service) time is strongly dependent on their size, as well as on the employed transmission rate.

Assumption 5.

In certain conditions, the connection between a sensor and the gateway may be closed during service (e.g., due to transmission link errors, etc.) after some time, which we assume random and exponentially distributed with the average of ν ⁻¹.

Further, we consider a rate-power function, which characterizes the data rate achievable at some transmit power level by a particular modulation and coding scheme.

Assumption 6.

The transmit power p of a sensor and its corresponding data rate r are coupled by the Shannon’s formula:

$$ r = w \log\left(1+\text{SNR}\right) =w \log\left(1+ \frac{p \gamma}{N_{0}} \right), $$

((1))

where p is the output power of the RF power amplifier, w is the channel bandwidth, and N ₀ is the noise power.

The parameter γ stands for the channel gain, and it is assumed to be proportional to some power function of the distance d between the transmitter and the receiver. We intentionally omit the consideration of fading effects, focusing on the long-term averages, as dictated by our problem formulation.

Assumption 7.

To avoid the unbounded growth of value γ, we impose a respective upper bound, such that the signal propagation is defined by:

$$ \gamma =\min \left(G d^{-k}, \gamma_{max} \right), $$

((2))

where G is the propagation constant, k is the propagation exponent. We note that γ _max characterizes the maximum level of SNR at the distance $d_{0} = \left (\frac {G}{\gamma _{\textit {max}}} \right)^{\frac {1}{k}}$, so that its further increase does not impact the achievable data rate.

Moreover, to make our model more realistic, we introduce practical restrictions on the transmit power, that is, the minimum and the maximum allowed power levels (p _min and p _max, respectively). Interestingly, the value of p _max defines the strict border on the distance between the sensors and the gateway. Equivalently, all sensors able to maintain the required data rate r should be positioned within a sphere of radius R:

$$ R = \left[ \frac{p_{\text{max}} G }{N_{0} \left(e^{\frac{r}{w}} - 1 \right)} \right]^{\frac{1}{k}}. $$

((3))

Additionally, from Assumption 1 on the PPP follows the uniform distribution in a certain area. Therefore, we may further consider that all of the sensors are uniformly distributed within the sphere of radius R.

By contrast to much conventional research, our model also enables a more detailed consideration of the power consumption based on a particular device profile.

Assumption 8.

Therefore, we assume the following levels of power consumption at the sensor:

(i) transmit power consumption P _tx =p+P _c , as long as the sensor is sending data, where p is the radiated power and P _c is the constant circuit power;

(ii) active power consumption P _a , when the sensor has some data and awaits service;

(iii) inactive power consumption, when the sensor has no data to transmit, which is assumed to be zero without loss of generality.

While the stored energy is spent according to the levels P _tx and P _a , we also employ the following model for energy replenishment (see Fig. 3). Generally, the WETI of a sensor converts the received RF energy to electrical energy, which is then stored in a rechargeable battery. Since the corresponding wireless technology is not widely available at the moment, below we assume a generic operation of the respective energy harvesting mechanism to provide a first-order evaluation. It can be further tailored to any practical operation, when the details of a particular energy transfer technology are known.

Fig. 3

Considered oracle-assisted model of continuous energy replenishment

Assumption 9.

We further assume that all the sensors may be split into L≥1 classes by the volume of the harvested RF energy h ₁,…,h _L , which are dependent on the particular WETI technology, as well as on the distance between the gateway and the sensor. For the energy replenishment process, we assume similar propagation parameters differentiating them with the index H (e.g., k _H , G _H , etc.).

The introduction of L classes reflects our intuition on that the volume of energy harvested by WETI drops with the increasing distance (if the energy transfer parameters at the gateway are fixed). However, the actual energy degradation parameters, as well as the exact practical number of such classes may be difficult to define before corresponding systems are implemented as prototypes (even though we share some of our thinking on this in Section 4). Regardless of the final properties shown by WETI, we propose that the system would realize priority service for the sphere edge sensors to balance the harvested/transmitted energy in the system.

To continue with the optimistic estimate for the power consumption/replenishment analysis, we consider that the system is managed by an oracle, which collects all the harvested energy by the sensors energy into a common rechargeable battery (buffer) of capacity B and then distributes it between the active sensors (see Fig. 3). Even though in a practical system such rechargeable batteries naturally belong to individual sensors, the assumption of a shared energy storage allows us to construct an upper bound on the feasible energy management mechanisms.

Generally, the total energy inflow into the oracle’s energy buffer should be modeled by an appropriate random process for the given type of energy source (e.g., WETI type). Here, to provide an initial evaluation, we assume that such replenishment process is continuous with the rate of h _l for the sensor of class l. Finally, the harvested energy is instantaneously stored in the common buffer and may be spent immediately after.

3 Our mathematical approach

The proposed analytical approach is constructed as follows. First, we address the spatial system features and deliver the stationary distribution of distances between the sensors and the gateway. This distribution translates into the respective distributions of transmit powers, as well as levels of harvested energy, depending on the distance to the gateway. Further, we provide an underlying solution, characterizing system behavior without energy replenishment accounted for, which essentially is the baseline Markov process. In some special cases (e.g., when all sensors belong to the same energy replenishment class, L=1), this process may be reformulated as one-dimensional birth-death process, or as an L-dimensional process for L>1. Basing on such process, we introduce a Markov process-driven fluid queue model of energy level evolution, which captures data transmission together with energy expenditure/replenishment processes.

3.1 Spatial characteristics

Given that the sensors are distributed in the 3D space around the gateway (see our system model in Section 2), we begin with calculating their average power consumption in the system. To this end, the distribution of the individual transmit powers may be given as follows.

Proposition 1.

The distribution of the individual transmit powers of the sensors is delivered by:

$$ f_{p}(p) = \left\{ \begin{array}{ll} p^{\frac{3}{k}-1} \cdot \frac{3}{k} (p_{\max})^{-\frac{3}{k}}, \quad p_{\min} \leq p < p_{\max},\\ F_{d}(p_{\min}) = \left[\frac{p_{\min}}{p_{\max}}\right]^{\frac{3}{k}}. \end{array} \right. $$

((4))

Proof.

We recall that the spatial locations of the sensors are distributed uniformly within a sphere of a known radius R centered at the gateway. Hence, we obtain the distribution of distances between a sensor and the gateway, that is, the distance between a random point within the sphere and its center. Due to the uniform distribution, the probability density function of sensor coordinates f(x,y) is defined as:

$$ f(x,y,z) = \frac{1}{\frac{4}{3}\pi R^{3}}, \quad x^{2}+y^{2}+z^{2} \leq R^{2}. $$

In terms of spherical coordinates, we establish:

$${} \begin{aligned} f_{r, \theta,\phi}(r,\theta,\phi) &=f_{x,y,z}(r\sin \theta \cos \phi, r \sin \theta \sin \phi,r \cos \theta) \cdot J \\&= \frac{3r^{2} \sin \theta}{4\pi R^{3}}, \quad 0< r < R. \end{aligned} $$

Therefore, we may conclude, that

$${} f_{\phi}(\phi) =\frac{1}{2 \pi}, 0 \leq \phi \leq 2\pi, \quad f_{\theta}(\theta) =\frac{1}{2} \sin \theta, 0 \leq \theta \leq \pi, $$

and, more importantly, the distribution of the distances to the center of the sphere is $f_{r}(r) = \frac {3r^{2}}{R^{3}}, 0< r < R$. It implies, in turn, that the distribution of the distances between an arbitrary sensor and the gateway is given by:

$$ f_{d}(d) = \frac{3d^{2}}{R^{3}}, 0< d < R, \quad F_{d}(d) = \frac{d^{3}}{R^{3}}, 0< d < R. $$

We also remind that the transmit power consumption is lower-limited by p _min. Without loss of generality, we specify that d ₀ (or γ _max-limited SNR, see above) corresponds to p _min (otherwise, it can be forced by updating, e.g., d ₀ to the maximum value).

Consequently, if the transmit power is defined as:

$${} {\fontsize{8.8pt}{9.6pt}\selectfont{\begin{aligned} p &= \max \left[\frac{N_{0}}{\gamma}\left(e^{\frac{r}{w}} -1\right), p_{\text{min}}\right] =\max \left[d^{k} \frac{N_{0}}{G }\left(e^{\frac{r}{w} }-1\right), p_{\text{min}}\right],\\ d &= \left[\frac{p G }{N_{0} \left(e^{\frac{r}{w}} - 1 \right)} \right]^{\frac{1}{k}}, \quad d'_{p} = \frac{1}{k} p^{\frac{1}{k}-1} \left[ \frac{G }{N_{0} \left(e^{\frac{r}{w}} - 1 \right)} \right]^{\frac{1}{k}}, \end{aligned}}} $$

and the individual power distribution may be given by (4).

Therefore, the average power may be delivered by:

$$ \begin{array}{l} \!\!\! E\left[p\right]=\!\!\! \int \limits_{p_{\min}}^{p_{\max}}\!\!\! pf_{p}(p) dp =\frac{kp_{\min}}{3+k} \left[\frac{p_{\min}}{p_{\max}} \right]^{\frac{3}{k}} + \frac{3p_{\max}}{3+k}, \end{array} $$

((5))

while the variance may be obtained via the second moment:

$$ \begin{array}{l} \!\!\! E\left[p^{2}\right]=\!\!\! \int \limits_{p_{\min}}^{p_{\max}}\!\!\! p^{2} f_{p}(p) dp =\frac{2kp^{2}_{\min}}{3+2k} \left[ \frac{p_{\min}}{p_{\max}} \right]^{\frac{3}{k}} +\frac{3p^{2}_{\max}}{3+2k}. \end{array} $$

((6))

We emphasize that in case when i independent sensors are active, the total consumed power is determined by the random value y _i , where $y_{i} =\sum _{n=1}^{\min (i,K)} p_{n}$ is the sum of random variables, the distributions of which we already know. The straightforward convolution of K distributions at hand is an intricate technical problem, and therefore, we suggest an alternative approach. Correspondingly, we calculate the convolution of two random variables and then investigate the distribution of y _i+1=y _i +p, given that y _i ,i≥2 may be approximated by the normal distribution. Using the expectation and the variance from (5)–(6) only, we arrive at a tight approximation for the distribution $F_{y_{i}}(y_{i})$ and the corresponding expectation E[y _i ], which we do not detail in the paper due to the simplicity of its calculation.

Generalization 1.

Importantly, by similar reasoning, we may characterize the arbitrary distribution of sensors around the gateway, if required in practice.

Further, given the estimated average power consumption of i sensors in service (as we were not focusing on those actively waiting, but on those transmitting), we find the random power consumption P _i for i sensors. Then, we include into consideration both transmit and circuit power levels of the random transmitting sensors, as well as the active power of the waiting sensors:

$$ P_{i} = \sum\limits_{n=1}^{i_{0}} \left(p_{n} +P_{c} - P_{a}\right)\frac{1}{K} +P_{a} i, \quad 0 \leq i \leq M, $$

where i ₀= min(i,K) is the number of transmitting sensors and p _n is the transmit power of a random sensor n located at a particular point. We remind that during the data transmission by the tagged sensor, the rest (i−1) sensors remain active according to the BLE operation and consume P _a power units.

Proposition 2.

The spatially averaged total system power consumption $\bar {p}_{i}$ may be determined by averaging the sum of point-dependent levels p _n ,0≤n≤i of individual sensor power consumptions:

$${} {\fontsize{8.9pt}{9.6pt}\selectfont{\begin{aligned} \bar{p}_{i} = \frac{1}{K} E\left[y_{\min(i,K)}\right] + \frac{\min(i,K)}{K} \left(P_{c} - P_{a}\right) +P_{a} i, \quad 0 \leq i \leq M, \end{aligned}}} $$

where $y_{i} \,=\,\sum _{n=1}^{i}\! p_{n}$ is the sum of random transmit powers.

Generalization 2.

For the case of L>1, there are simultaneously i=i ₁+…+i _L sensors of different energy replenishment classes in the system. Then, the spatially averaged power distribution would be dependent on the discipline of sensors’ service. For example, assuming priority in serving the sphere edge sensors, where the energy class corresponds to the priority class (i.e., worse energy conditions yield higher service priority), we arrive at the following formulation:

$$ \begin{aligned} \bar{p}_{i} &= \frac{1}{K} \left (E_{L}[y_{i_{L}}] + \ldots + E_{s}[y_{i_{s}}] +E_{s-1}[y_{K-i_{L}-\ldots-i_{s}}] \right) \\ &\quad+ \frac{\min(i,K)}{K} \left(P_{c} - P_{a}\right) +P_{a} i, \quad 0 \leq i \leq M, \end{aligned} $$

where $E_{s}[y_{i}] \!=E\left [\sum _{n=1}^{i}\! p_{n}\right ]$ is the expectation of the sum of random transmit powers across a given energy class. The solution for the system under this or any similar condition constitutes a separate technical problem and remains out of scope of this paper.

3.2 Underlying Markov process

In what follows, we incorporate a dynamic component into the static spatial model detailed above and apply the theory of Markov processes to analyze our system in its stationary mode. Hence, we consider the system at the embedded points, corresponding to the data arrivals and departures. According to the system model in Section 2, the data service time has an exponentially distributed random boundary with the rate ν, which gives a combination of service rates.

Proposition 3.

For one class of energy replenishment, the number of sensors in service N(t) constitutes the one-dimensional birth-death process (BDP, see “horizontal” process in Fig. 4) with the following transitions:

$$\begin{array}{*{20}l} \lambda_{i} &= \lambda(M-i), \\ \mu_{i}&= \left\{ \begin{array}{ll} i(r\mu+\nu),\quad \text{if}\ i < K, \\ K(r\mu+\nu),\quad\text{if}\ i\geq K. \end{array} \right. \end{array} $$

Fig. 4

Baseline hybrid continuous-discrete process

The infinitesimal generator Q(x)∈R ^(M+1)×(M+1), conditioning on the particular level of energy x>0, does not depend on x and is given as:

$$ {Q(x)}= \left(\begin{array}{cccccc} -\lambda_{0} & \lambda_{0}& 0&\ldots& 0& 0\\ \mu_{1} & - \left(\lambda_{1} + \mu_{1}\right)&\lambda_{1} & \ldots& 0& 0\\ & \ldots&& & & \\ 0 & 0& & \ldots& \mu_{M}& -\mu_{M} \end{array} \right). $$

Interestingly, in the particular case of unconstrained energy supply, the considered system corresponds to the M/M/K/M/M model. In the range 0≤n≤K<M, we obtain:

$$ \pi_{n} = \pi_{0} \left(\frac{\lambda}{r\mu} \right)^{n} {M\choose n}. $$

In the range K≤n≤M, the stationary probabilities are:

$$ \begin{array}{rl} \pi_{n} = \pi_{0} \left(\frac{\lambda}{r\mu} \right)^{n} {M\choose n} \frac{n!}{K!}K^{K-n}, \end{array} $$

where

$$ \pi_{0} =\left(1+ \sum\limits_{n=1}^{M} \pi_{n} \right)^{-1}. $$

Generalization 3.

Accounting for several classes of energy replenishment, we transform the BDP into an L-dimensional Markov process. Therefore, such finite Markov process has (M+1)^L states and the corresponding infinitesimal generator Q has to be updated accordingly.

For instance, the transitions of this process from Generalization 2 (when L>1) may be defined as:

$${} {\fontsize{8.4pt}{9.6pt}\selectfont{\begin{aligned} &\lambda_{(n_{1},\ldots,n_{L}) \rightarrow (n_{1},\ldots,n_{j}+1,..n_{L}) }= \lambda\left(M-\sum\limits_{i=1}^{L} n_{i}\right)q_{j}, \quad \sum\limits_{i=1}^{L} n_{i} \leq M, \\ &\mu_{(n_{1},\ldots,n_{L}) \rightarrow (n_{1},\ldots,n_{j}-1,..n_{L}) }= \min\!\left(\!n_{j},\max\left(\!K\,-\,\sum\limits_{i=j+1}^{L}n_{i},0\!\right)\!\right)(r\mu+\nu), \end{aligned}}} $$

((7))

where q _j is the probability that a sensor belongs to class j.

Even for the case of two classes, it is not straightforward (however, possible) to analyze the corresponding system in the closed form. The following approach allows numerical derivation of the system parameters for the arbitrary number of classes L.

3.3 Fluid queue model

Here, we concentrate on characterizing the energy-related properties of our system. We remind that every sensor is equipped with a rechargeable battery of volume B, which may provide energy for communication, as well as recharge with the energy harvested from the gateway by WETI. That, all together, may be represented as a fluid queue driven by the above BDP with the state N(t). Hence, we describe the fluid model by the process S(t)={N(t),X(t);t≥0} describing the state of the fluid queue with a single (common) energy buffer of the oracle, where X(t)∈[0,B] is the fluid level of the energy buffer at time t and B is the maximum fluid level (see example in Fig. 4). We note that the system driven by the one-dimensional process may be extended for the case of L>1.

The fluid level distribution is continuous across different fluid levels, but it also has probability masses at the border values. Following the approach in [29] and [28], we define the stationary fluid density π(x) and the fluid mass function c(x) as:

$$ \begin{array}{l} \pi_{i}(x) ={\lim}_{t \rightarrow \infty}{\lim}_{\Delta \rightarrow 0} \frac{\Pr\left\{ N(t) = j, x \leq X(t) <x+\Delta \right\}}{\Delta}, \\ c_{i}(x) ={\lim}_{t \rightarrow \infty} \Pr\left\{ N(t) = j, X(t)=x \right\}. \end{array} $$

Over the continuous interval (0,B), the following system of the ordinary differential equations (ODE) should be satisfied:

$$ \begin{aligned} r_{0} \frac{d}{dx} \pi_{0}(x) &= {\sum\nolimits}_{j=0}^{M}\pi_{j}(x) q_{j0}, \\[-5pt] &\quad\ldots\\[-5pt] r_{i}\frac{d}{dx} \pi_{i}(x) &= {\sum\nolimits}_{j=0}^{M}\pi_{j}(x) q_{ji},\\[-5pt] &\quad\ldots\\[-5pt] r_{M} \frac{d}{dx} \pi_{M}(x) &= {\sum\nolimits}_{j=0}^{M}\pi_{j}(x) q_{jM}, \end{aligned} $$

((8))

or, in the matrix form:

$$ \frac{d}{dx} \left(\pi(x) R \right)= \pi(x) Q, $$

((9))

where q _ij is an element of the infinitesimal generator matrix Q, π(x) ∈ R ^1×(M+1) is the transposed vector of stationary fluid density at the corresponding state, R∈R ^(M+1)×(M+1) is the diagonal matrix consisting of elements r _i , and r _i is the fluid rate at the state i. The fluid rate equals to:

$$ r_{i} = h - \bar{p}_{i}, $$

((10))

where h=M h ₁ is the average total harvested energy in case of L=1 energy replenishment classes. For the case of L>1, r _i is recalculated accordingly.

Importantly, in case of unconstrained energy supply we may set $\frac {d}{dx} \pi (x) = 0$. Therefore, the equation (9) transforms (for the stationary mode) into a system of linear equations, which follows from Chapman-Kolmogorov equations of Markov process described above and may be solved easily.

Proposition 4.

We assume here that the spatial distribution and the service process are independent, which may not seem practical, but in the longer-term and averaged across space yields similar results. Therefore, in our analysis we abstract away the exact history of sensor arrivals and departures. The latter effectively translates into the assumption that the distribution of powers by the transmitting sensors is time-independent and we consider all the transmitting sensors at the state i to be distributed according to (4). The spatially-averaged power consumption is calculated according to Proposition 2.

Furthermore, we note that the rate r _i depends on the number of active sensors and does not depend on the fluid level and, hence, continuous. Otherwise, it is possible to address a number of particular cases according to [29].

Generalization 4.

Assuming that all the sensors are divided into energy replenishment classes by a number of distance-based thresholds d ₁,…,d _L−1, where d _L =R, with the harvested power h ₁,..,h _L , we may further calculate the average harvested energy for the system E[h] by taking into account the distribution of sensors from different energy classes. By derivations similar to those in Proposition 2, we may estimate the average transmit power consumed by i sensors of a given class l basing on:

$${} {\fontsize{9.2pt}{9.6pt}\selectfont{\begin{aligned} E\left[p\right]= \int \limits_{p(d_{l-1})}^{p(d_{l})} pf_{p}(p) dp =\frac{k_{H}p(d_{l-1})}{3+k_{H}} \left[ \frac{p(d_{l-1})}{p(d_{l})} \right]^{\frac{3}{k_{H}}} + \frac{3p(d_{l})}{3+k_{H}}. \end{aligned}}} $$

The corresponding variance characterizing a sensor of class l may be found by using:

$$ E\left[p^{2}\right]= \frac{2k_{H}p^{2}(d_{l-1})}{3+2k_{H}} \left[ \frac{p(d_{l-1})}{p(d_{l})} \right]^{\frac{3}{k_{H}}} + \frac{3p^{2}(d_{l})}{3+2k_{H}}. $$

Further, we continue with the boundary conditions at the fluid levels 0 and B, which may be obtained similarly to [28]:

$$ \begin{aligned} -\pi_{j}(0) r_{j} + {\sum\nolimits}_{k} c_{k}(0) q_{kj}=0,\quad 0\leq j \leq M, \\ \pi_{j}(B) r_{j} + {\sum\nolimits}_{k} c_{k}(B) q_{ij}=0,\quad 0\leq j \leq M. \end{aligned} $$

((11))

The boundary conditions in (11) constitute an intricate part of our solution. In Appendix, we detail the proposed approach to establish the boundary conditions, solve a corresponding system of differential equations minding the specific operation of our system, and finally obtain the stationary probability distribution π(x) (27) and the probability mass functions c(0),c(B) (26).

Given the stationary distribution ${\pi }(x) = \sum _{i=0}^{M} \alpha _{i} e^{\lambda _{i} x}\phi _{i}$, we may derive the stationary distribution of the number of active sensors as:

$$\begin{array}{@{}rcl@{}}{} \pi_{i}^{(n)} &=& {\lim}_{t \rightarrow \infty}{\Pr \left \{ N(t) = j \right\}} = c_{i}(0) + c_{i}(B) + \int\limits_{0}^{B} \pi_{i}(x)dx \\[-5pt] &=& c_{i}(0) + c_{i}(B) + \sum\limits_{i=1}^{M} \frac{\alpha_{i}}{\lambda_{i}}\phi_ i \left(e^{\lambda_{i} B} -1 \right). \end{array} $$

((12))

Here, we omit a separate consideration of λ=0 (all the corresponding integrals may be easily obtained for the term α _i ϕ _i instead of $\phantom {\dot {i}\!}{\alpha _{i}}\phi _{i} e^{\lambda _{i}x}$).

In summary, we calculate several important stationary state metrics in our system. Basing on the stationary distribution π(x), we produce the average number of active sensors $\bar {N}$ and the average data delay $\bar {T}$ as:

$$ \bar{N} =\sum\limits_{i=1}^{M}i \pi_{i}^{(n)}\quad \text{and}\quad\bar{T}=\frac{\bar{N}}{\lambda\left(M-\bar{N}\right)}. $$

((13))

Similarly, we may calculate the waiting time and the time interval between the arrival and the beginning of service. The probability of a full or empty energy buffer follows immediately from the vector x (26):

$$ P_{\text{empty}} =\sum\limits_{i=0}^{M} c_{i}(0)\quad \text{and}\quad P_{\text{full}}=\sum\limits_{i=0}^{M} c_{i}(B). $$

((14))

Finally, the average level of remaining battery charge may be established from the following:

$$ \bar{X}= B P_{\text{full}} + \sum\limits_{i=0}^{M} {\int\limits_{0}^{B}}x \pi_{i}(x)dx. $$

((15))

4 Summary of technology-based considerations

This section summarizes our considerations on the state-of-the-art BLE technology and how it translates to our system model in Section 2, as well as on the potential performance of wireless energy transfer. In particular, we investigate the technology-related limitations of WETI together with the most prominent solutions to improve its performance.

4.1 Communication over BLE technology

The below text describes the operation of the BLE technology in more detail, so that to substantiate the system model built in Section 2. On the physical layer (PHY), BLE operates in the unlicensed industrial, scientific, and medical (ISM) 2.4 GHz band and uses the Gaussian frequency shift keying (GFSK) with the bandwidth-time product equal to 0.5 and the symbol rate of 1 mega-symbol per second. This leads to a fixed transmission rate r in Assumption 2.

For a BLE transmitter, according to the current specification [3], the power of output radio signal should range from −20 to +10 dBm and the sensitivity level of the BLE receiver should be below −70 dBm, which translates to the limits p _min and p _max on the allowed transmit power. In order to simplify the design of the transceivers, the maximum length of PHY BLE packets with all the headers is set at 47 bytes [6].

The available frequency band is divided into 40 2-MHz wide channels. Three of the channels are assigned specifically for advertising and discovery of the services and are called advertising channels. The remaining 37 data channels might be used for peer-to-peer data transfers. We envision that in practice, fewer channels may actually be available for smart home operation (as some of the channels may be reserved by the human user) and thus arrive at the variable number of channels K in Assumption 3. The data transfer between BLE devices is bound to time units known as advertising and connection events.

The advertising events are used by the BLE transceivers to broadcast small blocks of data or to agree on the parameters of the connection to be established in the data channels. As detailed in Fig. 5, at the beginning of an advertising event, the advertiser, i.e., the device that has some data to transmit, sends an advertising frame. Using the different advertising frame types, the advertiser either encapsulates up to 31 bytes of data directly in the advertisement frame or announces its willingness to establish a connection in data channels.

Fig. 5

Illustration of BLE advertising events and related timing

In the latter case, having transmitted the frame, the advertiser switches to receive mode and waits for possible connection establishment requests. If the connection request from a device (which is referred to as an initiator) is received, the two devices start peer-to-peer connection in the data channels. Depending on the specifics of the implementation and the application, a BLE advertiser might either send its advertisements in a single advertising channel or sequentially switch between different advertising channels. The minimum time period between the starts of two consecutive advertising events is defined as:

$$T_{\text{advEvent}}= \text{advInterval} + \text{advDelay}, $$

where advInterval is an integer multiple of 0.625 ms in the range of 20 ms to 10.24 s and advDelay is a pseudo-random value with a range of 0 to 10 ms generated anew for each advertising event. The above procedure is used in practice to construct the transmission schedule at the gateway, and we choose to omit the consideration of the scheduling operation in this paper for the sake of simplification. Above, we provide a bound on the operation of any practical radio resource management mechanisms employing the advertising channels.

Once a connection over a data channel is established, the initiator becomes the master and the advertiser becomes the slave. As shown in Fig. 6, at the beginning of each connection event (referred to as the connection event anchor point), the used radio channel is changed following the pre-agreed sequence. The communication in each connection event is started by the master, which sends a frame to the slave. The master and the slave alternate sending the frames on the same data channel while at least one of the devices has data to transmit or until the end of the current connection event. This yields our assumption on non-interfering and synchronous data channels in Assumption 3.

Fig. 6

Illustration of BLE data transmission in data channels and related timing

In the case if either the master or the slave receives two consecutive frames with CRC errors, the connection event is closed. The same happens if either device misses a radio packet. Once the connection event is closed, both master and slave might switch to low-power mode up to the start of the next connection event. The parameters of the connection (e.g., the connection event interval – connInterval or the list of used data channels) might be updated dynamically. The connection is closed either by the devices once the link is not required any more or automatically due to the connection timeout (ranges from 100 ms to 32 s). All of the above motivates our Assumption 5 on the connection closing interval.

In more detail, the timing of connection events is determined using two parameters, namely, the connInterval, and the slave latency (connSlaveLatency). The connInterval is a multiple of 1.25 ms and has values ranging from 7.5 ms to 4.0 s. The connSlaveLatency (connSlaveLatency≤500) defines the maximum number of consecutive connection events in which a slave device is not required to listen to the master and enables energy saving. The period between the frames on the same data channel is equal to the Interframe Space period (IFS) set at 150 μ s. The maximum link layer (LL) payload of a BLE data frame is just 27 bytes [6, 30]. These timings help us produce numerical results in the following section.

The illustration of the entire procedure including advertising, connection establishment, and data transfer is provided in Fig. 7.

Fig. 7

Example advertising, connection establishment, and data transfer operation of BLE

4.2 Technical limitations of wireless energy transfer

The state-of-the-art on the radio energy powered radios today is well illustrated by Figure 1 in [31]. Accordingly, to enable wireless energy transfer, the power of the received radio signal should be above −30 dBm. The reason for this is that less powerful radio signals are unable to provide the minimum voltage required to commutate the stages of diode rectifiers, which are used for retrieving the energy of the received radio signal.

Another major challenge related to wireless energy harvesting/transfer is that the efficiency of RF-DC conversion strongly depends on the power of the received radio signal (see, e.g., Figures 17 and 19 in [32] or Figure 3 in [33]). Whereas the latter problem has been partially solved in [34] by employing rectifiers with the reconfigurable number of stages, the efficiency of corresponding systems still does not exceed 60 %.

Finally, in practical cases, the efficiency of energy transfer from the rectifier’s output to the energy buffer is also below 100 %, and it varies for different rectifier output voltages (e.g., the work in [35] indicates the best-case efficiency of around 90 %).

In light of the above technological limitations and accounting for the restrictions by the frequency regulation agencies (i.e., Finnish Communications Regulatory Authority, FICORA) on the use of the unlicensed ISM channels [36], we decided to investigate the following important questions: a) which frequency bands would be most attractive for WETI (wireless energy transfer interface), b) what is the maximum feasible distance of wireless energy transfer, and c) what is the actual amount of energy one can potentially transfer from the home gateway to the sensors.

The input data used for our calculations is summarized in Table 1 and the results are detailed in Table 2.

Table 1 Input parameters for WETI efficiency calculations

Table 2 Expected efficiency of practical WETI operation

Our results indicate that for high frequencies, due to lower maximum permitted transmit power and higher path losses, it is practically impossible to enable sustainable WETI operation for the distances above 1 m. Further, although low frequencies have sufficient maximum permitted transmit power and significantly lower path losses, the major challenge in employing those for energy transfer is the design of the antennas. In our calculations, we assumed the halfwave dipole antennas. As revealed in Table 2, the length of such an antenna, e.g., for 40 MHz band, is about 3.6 m. Hardly, the antenna of such length is feasible for most of contemporary smart home applications.

Meanwhile, the use of smaller antennas would significantly reduce the energy of received radio signal and might, in turn, introduce some undesired directivity issues as well. However, an attractive compromise between the maximum permitted transmit power, path losses, and antenna efficiency is achieved in 869.5 MHz band (see Table 2). Nonetheless, even there, the gateway transmitting with the full permitted power might currently provide only 0.03 % of the energy required for continuous operation of a BLE transceiver located at a distance of 5 m away. Note that in our calculations, we used the free-space propagation model and the parameters corresponding to the best reported research prototypes, and thus the results for a real-life system might be even less optimistic.

Therefore, while enablement of RF energy transfer promises many inspiring benefits, there is much further work required for making such systems a reality. Among the primary technical challenges is: how to produce the signal strong enough to enable energy harvesting at the output of the sensor’s antenna. This problem might be approached from different angles. First, if we are to go for the sub-100 MHz band with lower path losses, there is a pressing need for the high-efficient small-size antennas, which are not available today. The second opportunity is to stay at higher frequencies, but to increase the power of the energy carrying radio signal. This might be done by updating the conservative limitations in the currently available bands, as well as by introducing novel sub-band(s), specifically dedicated for wireless energy transfer.

Both paths require the involvement of the frequency regulation authorities and demand further studies on the possible effects of the radio waves on the environment and human health in order to determine the permitted power. Alternatively, the WETI efficiency problem might be solved by outlining the way to reduce the minimum required level of the radio signal to enable energy harvesting (i.e., well below −30 dBm). This can be done either via improving parameters of the diodes used in current rectifier designs or by suggesting new rectifier architectures. Finally, the most realistic option considering the current state of the technology seems to employ highly directional antennas. This, in turn, would require simple and precise node localization mechanisms, as well as technology to enable cheap and simple antennas with a changeable radiation pattern. However, we are hopeful that the outlined technical issues are solved in the very near future.

5 Selected numerical results and conclusions

In this section, we introduce selected numerical results detailing the performance of the smart home system investigated above. A realistic operation scenario is considered with the parameters following from Section 4 (see Table 3) and analyzed by means of the proposed oracle-based approach. We remind that the essence of the oracle is to combine and centrally manage all the energy flows in our system, and the proposed model thus delivers an optimistic performance estimate to a practical distributed sensor network. Additionally, our model employs realistic considerations on spatial deployment of sensors and their data dynamics. The main purpose of this study is thus to illustrate the potential of our system behavior under the various conditions, rather than conduct exact performance evaluation for a fixed set of parameters.

Table 3 Employed numerical parameters

In the following results, we outline the unstable region of system operation, which corresponds to the situation when the power outage probability P _empty or, literally, the battery drain probability (when the energy level X hits its minimum of 0) exceeds some threshold ε for particular data and energy arrival rates. In particular, if P _empty>ε, then the respective operation region corresponds to the system-wide energy outage, which in practice yields frequent events when there is no energy at the sensor to continue communication. Even though sensors might be equipped with pre-charged batteries (which would translate into initial charge level in our model), system instability threatens sustainable sensor operation and may drain its battery power quickly.

We continue by illustrating the unstable region of operation in Fig. 8 for several system profiles (corresponding to different numbers of sensors and available data channels). In these plots, we show the region of interest bounded by the harvested power per sensor, as well as by the data arrival rate per sensor. By the power outage probability calculated as (26), the power rate region splits into two parts, one of which corresponds to the region of unstable system operation with a threat of a system-wide energy outage. Naturally, the region of instability expands with either increase in the number of users or decrease in the number of channels, but our approach allows us to quantify it exactly for the given system parameters.

Fig. 8

System-wide energy outage region for ε=0.01. a M = 100, K = 17. b M = 100, K = 37. c M = 50, K = 37

Fig. 9 continues our evaluations and demonstrates an example of actual relation between the estimated power outage, the estimated power charge level (15) normalized to the maximum battery capacity, and the power/rate parameters. The constructed surface clearly indicates that there is a set of outage probability levels, including the case when the buffer empties in the stationary mode. We note that computer simulation of such process may be non-trivial and definitely is time-consuming due to longer transition times. On the contrary, our proposed analysis may easily shed light on the system operation in the stationary mode.

Fig. 9

An illustration of the estimated power outage probability and power charge level (normalized by the maximum capacity) evolution

In summary, we conclude that our mathematical approach may be very useful to provide initial understanding of the smart home system operation with energy transfer/harvesting functionality for a wide range of practical parameters. Its applicability goes far beyond simple delay and load estimations, allowing to investigate and predict more complex system events, such as a large-scale power outage. We expect the proposed framework to be applicable at the design time of the future smart home gateway deployments enhanced with wireless energy transfer capability.

In the future, we plan to prototype the envisioned smart home system in hardware and test it in practice. This will enable us to study the validity of our current assumptions, as well as provide data for calibrating our model. The first prototype will be based on the currently available commercial radio energy transfer solutions for 868 MHz band and off-the-shelf BLE systems-on-chip. The conceptual functionality of the prototype system is as follows.

Once acquiring some sensed data for the gateway, a sensor enables its BLE transceiver and starts advertising. The event interval is set based on the energy arrival rate. Using the information on the sensor data at each node, its importance for the over-the-top applications, and accounting for the available sensors energy, the gateway builds the schedule for its BLE transceivers. Following the defined schedule, the transceivers establish the connections with the sensors over the data channels and receive the data. Once all the data from a sensor are transferred, the node switches to the low-power mode until a new portion of data arrives.

Note that the gateway may be made aware of the coexisting wireless networks (e.g., WiFi, ZigBee, etc.) and might thus ban the use of radio channels utilized by those for BLE. Given that the BLE transceivers on the gateway are properly synchronized, only the casual interference might affect communication. In our experiments, we plan to study the operation of the system both for the nodes with periodic (e.g., temperature/humidity sensors) and event-based (e.g., movement sensors) traffic generation.

6 Appendix

This appendix describes the process of obtaining solution to the differential equations under the given boundary conditions. Importantly, depending on the system parameters for some state j ₀, it may happen that $r_{j_{0}}(0) = 0$. In existing literature, this situation is referred to as zero state and requires special treatment.

Following the specific character of our problem, we may assume that there are three levels of r _j : negative, positive, and zero. We introduce the corresponding sets of indexes: (.)⁺,(.)⁰,(.)⁻ with respect to the sign of r _j and denote π ⁽⁺⁾,π ⁽⁰⁾,π ⁽⁻⁾ as the probability density of the states j with the respective level of r _j . Therefore, we rewrite the system of equations (9) as:

$${} {\fontsize{9.4pt}{9.6pt}\selectfont{\begin{aligned} \frac{d}{dx} \left(\pi^{+} (x) R^{+} \right) &= \pi^{+} (x) Q^{(++)} \! + \pi_{j_{0}} (x) Q^{(0+)} + \pi^{-} (x) Q^{(-+) }, \\ 0&= \pi^{+}(x) Q^{(+0)} + \pi_{j_{0}}(x) q_{j_{0},j_{0}} + \pi^{-}(x) Q^{(-0)}, \\ \frac{d}{dx} \left(\pi^{-} (x) R^{-} \right) &= \pi^{+} (x) Q^{(+-)} \! + \pi_{j_{0}} (x) Q^{(0-)} + \pi^{-} (x) Q^{(--)}, \end{aligned}}} $$

((16))

where Q ^(..),R ^(.) are the sub-matrices of Q and R, corresponding to the states with negative, zero, and positive rates.

From the second matrix equation, we obtain:

$$ \pi_{j_{0}}(x) = - \frac{1}{q_{j_{0},j_{0}}}\pi^{+}(x) Q^{(+0)} - \frac{1}{q_{j_{0},j_{0}}} \pi^{-}(x) Q^{(-0)}, $$

((17))

where $q_{j_{0},j_{0}} = - \left (\lambda _{j_{0}} + \mu _{j_{0}}\right)$.

The expression (17) may be substituted into the system (16), such that:

$${} \begin{aligned} \frac{d}{dx}\left(\pi^{+}(x) R^{+}\right)= \pi^{+}(x) \! \left[ Q^{(++)}-\frac{1}{q_{j_{0},j_{0}}}Q^{(+0)}Q^{(0+)} \right] +\\ + \pi^{-}(x)\! \left[ Q^{(-+) }- \frac{1}{q_{j_{0},j_{0}}}Q^{(-0)}Q^{(0+)} \right],\\ \frac{d}{dx}\left(\pi^{-}(x) R^{-}\right)= \pi^{+}(x)\! \left[ Q^{(+-)} - \frac{1}{q_{j_{0},j_{0}}}Q^{(+0)}Q^{(0-)} \right] + \\ + \pi^{-}(x)\! \left[ Q^{(--)}- \frac{1}{q_{j_{0},j_{0}}} Q^{(-0)} Q^{(0-)} \right], \end{aligned} $$

which is equivalent to solving a system:

$$ \frac{d}{dx} \left(\tilde{\pi}(x) \tilde {R}\right)= \tilde{\pi} \tilde{Q} \quad\text{or}\quad \frac{d}{dx} \left(\tilde{\pi}(x) \right)= \tilde{\pi} \tilde{Q} \tilde{R}^{-1}, $$

((18))

where $\tilde {Q} \in \mathbf {R}^{M\times M}$ is the corresponding matrix:

$${} \left(\begin{array}{rr} Q^{(++)}-\frac{1}{q_{j_{0},j_{0}}}Q^{(+0)}Q^{(0+)} &Q^{(+-)} - \frac{1}{q_{j_{0},j_{0}}}Q^{(+0)}Q^{(0-)} \\ Q^{(-+) }- \frac{1}{q_{j_{0},j_{0}}}Q^{(-0)}Q^{(0+)} &Q^{(--)}- \frac{1}{q_{j_{0},j_{0}}} Q^{(-0)} Q^{(0-)} \end{array} \right), $$

and $\tilde {R}$, $\tilde {\pi }$ are obtained by excluding the components of zero state.

6.1 Boundary conditions transform

The boundary conditions at the fluid levels 0 and B may be rewritten as:

$$ \begin{aligned} \pi^{+}(0) R^{+} &= c^{+}(0) Q^{(++)} + c_{j_{0}}(0) Q^{(0+)} + c^{-}(0) Q^{(-+) }, \\ 0 &= c^{+}(0) Q^{(+0)} + c_{j_{0}}(0) Q^{(00)} + c^{-}(0) Q^{(-0) }, \\ \pi^{-}(0) R^{-} &= c^{+}(0) Q^{(+-)} + c_{j_{0}}(0)Q^{(0-)} + c^{-}(0) Q^{(--) }, \end{aligned} $$

((19))

and

$${} \begin{aligned} -\pi^{+} (B) R^{+}&= c^{+}(B) Q^{(++)} + c_{j_{0}}(B) Q^{(0+)} + c^{-}(B) Q^{(-+) }, \\ 0&=c^{+}(B) Q^{(+0)} + c_{j_{0}}(B) Q^{(00)} + c^{-}(B) Q^{(-0) }, \\ -\pi^{-}(B) R^{-}&= c^{+}(B) Q^{(+-)} + c_{j_{0}}(B) Q^{(0-)} + c^{-}(B) Q^{(--) }, \end{aligned} $$

((20))

so that c ^.(0) and c ^.(B) are the parts of the transposed vectors of probabilities (corresponding to empty and full buffer c(0) and c(B)∈R ^(M+1)×1, respectively).

Further, the probability that the energy buffer is empty or full is zero in cases when rate increases (positive indexes) or decreases (negative indexes), correspondingly:

$$ c^{+}(0) = 0~ \text{and }~ c^{-}(B) = 0. $$

Therefore, transforming the boundary conditions (19), (20), and expressing $c_{j_{0}}(0)$ and $c_{j_{0}}(B)$, we obtain:

$$ \begin{aligned} \pi^{+}(0) R^{+}&= c^{-}(0) \left[ - \frac{1}{q_{j_{0},j_{0}}} Q_{0}^{(-0) }Q^{(0+)} + Q^{(-+) }\right], \\ \pi^{-}(0) R^{-}&= c^{-}(0) \left[- \frac{1}{q_{j_{0},j_{0}}} Q^{(-0) }Q^{(0-)} + Q^{(--) }\right],\\ c_{j_{0}}(0) &= - \frac{1}{q_{j_{0},j_{0}}}c^{-}(0) Q^{(-0) }, \end{aligned} $$

and

$$ \begin{aligned} \pi^{+}(B) R^{+}&= -c^{+}(B) \left[ Q^{(++)} - \frac{1}{q_{j_{0},j_{0}}} Q^{(+0)} Q^{(0+)}\right], \\ \pi^{-}(B) R^{-}&= -c^{+}(B) \left[ Q^{(+-)} - \frac{1}{q_{j_{0},j_{0}}} Q^{(+0)} Q^{(0-)} \right], \\ c_{j_{0}}(B) &= - \frac{1}{q_{j_{0},j_{0}}} c^{+}(B) Q^{(+0)}, \end{aligned} $$

which may be rewritten more compactly as:

$$ \begin{array}{l} \!\! \tilde{\pi}(0)\tilde{R}= c^{-}\!(0) \left[Q_{1}, Q_{2}\right],\quad\tilde{\pi} (B) \tilde{R}= c^{+}(B) \left[Q_{3}, Q_{4}\right ]\!, \end{array} $$

((21))

where

$$ \begin{aligned} Q_{1} &= - \frac{1}{q_{j_{0},j_{0}}(0)} Q^{(-0) }(0)Q^{(0+)}(0) + Q^{(-+) }(0), \\ Q_{2} &=- \frac{1}{q_{j_{0},j_{0}}(0)} Q^{(-0) }(0) Q^{(0-)}(0) + Q^{(--) }(0),\\ Q_{3} &= -Q^{(++)} + \frac{1}{q_{j_{0},j_{0}}} Q^{(+0)} Q^{(0+)},\\ Q_{4} &= -Q^{(+-)} + \frac{1}{q_{j_{0},j_{0}}} Q^{(+0)} Q^{(0-)}. \end{aligned} $$

6.2 Solving the system of differential equations

The solution to the system of first-order differential equations (18) may be obtained in the form $\tilde {\pi }(x) = e^{\lambda x} \phi$, ϕ∈R ^1×M [28], such that after substituting it back into our system we arrive at the characteristic equation:

$$ \phi \left(\lambda I - \tilde{Q} \tilde{R}^{-1} \right) = 0, $$

((22))

where I is the identity matrix, λ is the eigenvalue of the matrix $\tilde {Q} \tilde {R}^{-1}$, and ϕ is the respective transposed eigenvector.

Consequently, we represent our solution in the form:

$$ \tilde{\pi}(x) = \sum\limits_{i=1}^{M} \alpha_{i} e^{\lambda_{i} x} \phi_{i}, $$

((23))

where λ _i and ϕ _i are the eigenvalues and the corresponding transposed eigenvectors of the matrix (Q R ⁻¹), α _i are auxiliary coefficients, which may be obtained from the initial conditions.

Hence, in order to find the coefficients α _i , we substitute (23) into the boundary conditions (21):

$$ \begin{aligned} &{\sum\nolimits}_{i=1}^{M} \alpha_{i} \phi_{i} \tilde{R}= c^{-}\!(0) \left[Q_{1}, Q_{2} \right],\\ &{\sum\nolimits}_{i=1}^{M} \alpha_{i} e^{\lambda_{i} B} \phi_{i} \tilde{R}= c^{+}(B) \left[ Q_{3}, Q_{4} \right]. \end{aligned} $$

((24))

In the linear notation, the latter expression modifies as follows:

$$ \begin{aligned} &{\sum\nolimits}_{i=1}^{M} \alpha_{i} \phi_{1i} \tilde{r}_{1}= {\sum\nolimits}_{k=1}^{M}c^{*}_{k} g_{k1},\\ &{\sum\nolimits}_{i=1}^{M} \alpha_{i} \phi_{ji} \tilde{r}_{j}= {\sum\nolimits}_{k=1}^{M}c^{*}_{k} g_{kj},\\[-5pt] &\ldots,\\ &{\sum\nolimits}_{i=1}^{M} \alpha_{i} \phi_{1i} \tilde{r}_{1} e^{\lambda_{i} B}= {\sum\nolimits}_{k=1}^{M}c^{*}_{k} g_{k(M+1)},\\ &{\sum\nolimits}_{i=1}^{M} \alpha_{i} \phi_{ji} \tilde{r}_{j} e^{\lambda_{i} B}= {\sum\nolimits}_{k=1}^{M}c^{*}_{k} g_{k(M+j)},\\[-5pt] &\ldots \end{aligned} $$

((25))

where $c^{*}_{k}$ are the elements of the unknown row c ^∗=(c ⁻ (0),c ⁺ (B))∈R ^1×M, ϕ _ij are the elements of the matrix of eigenvectors Φ∈R ^M×M, and g _ij are the elements of the matrix G:

$$ {G}= \left(\begin{array}{cccc} Q_{1} &Q_{2}&0&0\\ 0&0&Q_{3}&Q_{4} \end{array} \right). $$

The linear system (25) contains exactly 2M equations together with 2M-dimensional vector of unknown variables:

$$ x = \left(\alpha_{1},\ldots,\alpha_{L},c^{*}_{1}\ldots,c^{*}_{M} \right). $$

((26))

Therefore, in order to derive all the unknown variables, we solve the system A x=0, where:

$$ A = \left(\begin{array}{ccccccc} \phi_{11} \tilde{r}_{1} &\ldots&\phi_{1L} \tilde{r}_{1}\\ &\ldots&\\ \phi_{M1} \tilde{r}_{M} &\ldots&\phi_{ML} \tilde{r}_{M}\\ \phi_{11}e^{\lambda_{1} B} \tilde{r}_{1} &\ldots&\phi_{1L} e^{\lambda_{L}B}\tilde{r}_{1}\\ &\ldots&\\ \phi_{M1}e^{\lambda_{1}B} \tilde{r}_{M} &\ldots &\phi_{ML} e^{\lambda_{L}B}\tilde{r}_{M}\\ \end{array} \left| \begin{array}{cc} -{Q_{1}^{T}}&0\\ -{Q_{2}^{T}}&0\\ 0&-{Q_{3}^{T}}\\ 0&-{Q_{4}^{T}}\\ \end{array} \right. \right). $$

In addition, we employ the auxiliary equation following from the normalization conditions:

$$ \sum\limits_{i=1}^{M} c^{*}_{i} +\sum\limits_{i=0}^{M} {\int\limits_{0}^{B}} \pi_{i}(x)dx +c_{j_{0}}^{-}(0) +c_{j_{0}}^{+}(B) = 1. $$

The solution $\tilde {\pi }(x)$ to the system above completely determines the solution π(x):

$$ \tilde{\pi}(x) = \sum\limits_{i=1}^{M} \alpha_{i} e^{\lambda_{i} x} \phi_{i}, $$

((27))

$$ \pi_{j_{0}}(x) = - \frac{1}{q_{j_{0},j_{0}}}\pi^{+}(x) Q^{(+0)} - \frac{1}{q_{j_{0},j_{0}}} \pi^{-}(x) Q^{(-0)}. $$

((28))

We note that if there is no zero state, we exploit the same way of solving the differential equation, by assuming Q ^(..)=0 corresponding to zero state. Another technical issue that might cause numerical instability is the low rank of the matrix A, when a particular element becomes large enough (e.g., for positive eigenvalues). This may be avoided by setting the corresponding α=0 and reformulating A.

By means of the above calculations, we establish the stationary distribution π(x) from (27) and the mass probability functions c(0),c(B) in (26), such that all the stationary parameters in our system may finally be derived.