Partial Diffusion Markov Model of Heterogeneous TCP Link: Optimization with Incomplete Information

Abstract

The paper presents a new mathematical model of TCP (Transmission Control Protocol) link functioning in a heterogeneous (wired/wireless) channel. It represents a controllable, partially observable stochastic dynamic system. The system state describes the status of the modeled TCP link and expresses it via an unobservable controllable MJP (Markov jump process) with finite-state space. Observations are formed by low-frequency counting processes of packet losses and timeouts and a high-frequency compound Poisson process of packet acknowledgments. The information transmission through the TCP-equipped channel is considered a stochastic control problem with incomplete information. The main idea to solve it is to impose the separation principle on the problem. The paper proposes a mathematical framework and algorithmic support to implement the solution. It includes a solution to the stochastic control problem with complete information, a diffusion approximation of the high-frequency observations, a solution to the MJP state filtering problem given the observations with multiplicative noises, and a numerical scheme of the filtering algorithm. The paper also contains the results of a comparative study of the proposed state-based congestion control algorithm with the contemporary TCP versions: Illinois, CUBIC, Compound, and BBR (Bottleneck Bandwidth and RTT).

1. Introduction

Despite its age of almost 50 years, the Transmission Control Protocol (TCP) [1] is still an object of permanent modernization and improvement, and this evolution represents a natural perpetual process. The root of this development lies in incessant challenges caused by a wide variety of computer networks, impetuous progress in the communication devices design, and strengthening of requirements to the information transmission [2,3,4]. Meanwhile, guaranteeing data transfer independent of the hardware platform is the key task of the TCP algorithm; both the stable functioning and effective use of the available channel bandwidth are also the performance characteristics of each specific version of TCP. The congestion control algorithms are responsible for the implementation of all these functions. They use two characteristics as the control actions. The basic one is the congestion window size (cwnd), i.e., the number of packets sent without acknowledgment. A less influential one is the retransmission timeout, i.e., some waiting time for the acknowledgment of the successful packet reception, which excess is treated by the congestion control algorithm as a packet loss.

When most channels were wire channels and had a relatively small capacity and queue waiting time “Additive Increase–Multiple Decrease” (AIMD) congestion control rule demonstrated good performance. This presumed a linear growth of the cwnd between two successive packet losses when the cwnd abruptly decreased in a jump-like manner. The effectiveness of this strategy for such channels was transparent. First, the small channel capacity gave a chance to reach a bandwidth limit linearly without losses for a rather short time. Second, wired hops were so reliable that the fact of a sudden packet loss presumed congestion at some “bottleneck” almost surely. Therefore, the loss indicated the necessity to reduce the sending rate. This simple reason was a base to develop such loss-based versions of TCP as Tahoe, New Reno, etc. [5].

In the case of the “long fat” channels (ones with huge capacity and long queue waiting times), AIMD-based versions of TCP turned out to be ineffective: they underused the channel bandwidth significantly. In the case of the channels with high capacity, the linear growth does not allow for the congestion window to quickly achieve values close to the available bandwidth. Plus, a loss of at least one packet decreases the data transferring speed even more. In addition, if a channel includes a wireless hop, facts of single packet losses are not an explicit congestion indicator. The round-trip time (RTT) parameter starts to play a remarkable role in the congestion control algorithm, and this brings to the variety of the TCP versions: delay-sensitive, hybrid loss-delay, bandwidth estimation-based, etc. [2]. All the modifications make the congestion control algorithm more tolerant to packet losses: after each loss, it decreases cwnd not multiplicative but more sparingly. At the same time, the cwnd growth speed is more aggressive to reach the channel bandwidth faster. The bandwidth value is unknown but estimated given all past statistics of the channel functioning. The algorithm probes more or less gentle cwnd enlargement to give a chance to use all channel resources. Hence, the typical cwnd curve between two packet losses demonstrates a concave [6] or mixed concave-convex character [7].

The ubiquitous application of wireless technologies in computer networks is a challenge to TCP protocol performance and claims its subsequent enhancement. Jitter and periodical signal fading in the wireless channel hops are extra sources of uncertainty of the channel real throughput. These physical phenomena affect both the new mathematical models of the channel functioning and the congestion control algorithms.

Mathematical models of computer network traffic are also developed intensively. With no goal to present a comprehensive overview of these models, we only mention their major classes

• Markov and hidden Markov models [8,9,10,11],
• queuing systems [12,13],
• models, based on the fluid or diffusion approximation of jump processes [14,15,16],
• network calculus models [17,18,19],
• models involving selfsimilar processes [20,21,22],
• concurrent models and games [23,24,25], etc.

Generally speaking, a prospective mathematical model of a channel should satisfy the conditions below.

• A model should describe the data transferring process adequately.
• A model should represent a trade-off between a complicated object with many parameters, their uncertainty along with the uncertainty introduced by the external disturbances, and simplicity.
• A model should operate with the same collection of statistical information as the one available in the real channel.
• A model should provide a possibility to simulate the collection of recent “concurrent” versions of TCP.
• The chosen model presumes the presence of the developed mathematical framework for the solution to the complex of all the analysis, estimation/identification and optimization/control problems. Availability of both the theoretical solution to the problems above and their efficient numerical realization is strongly encouraged.

The aim of the paper is two-fold. First, this is a presentation of a new mathematical model of the TCP link functioning based on the heterogeneous (wired/wireless) channel. It represents a controllable, partially observable stochastic dynamic system. The system state describes the status of the modeled TCP link and expresses it via the controllable Markov jump process (MJP) with a finite-state space. This space can be chosen arbitrarily depending on the desired detailing of the link description. Below in this paper, we consider four possible channel states:

• $e_1$: the channel is idle,
• $e_2$: the channel is loaded moderately,
• $e_3$: congestion in the wired segment,
• $e_4$: signal fading in the wireless hop.

Looking rather simple, this model admits successful description of such a problematic link phenomenon as congestion in a channel “bottleneck” and the carrier radio signal fading.

The observations included into the model correspond to those available to a TCP control algorithm on the sending side. Two observable processes describe the flow of packet losses and the flow of timeouts. They are represented by controllable Cox processes with intensity that depends both on the control and unobserved link state. The third observation is a flow of the acknowledgments concerning the successful packet reception on the receiving node. The flow is expressed in terms of a compound Poisson processes (CPP). Its first component represents a counting process of acknowledgment reception moments, and the second one registers corresponding individual values of the Round-Trip Time (RTT).

In the paper, we control the TCP varying the cwnd value only; however, the proposed model allows other control parameters, e.g., RTO (retransmission timeout). We also demonstrate how the proposed mathematical model can describe various contemporary versions of the TCP: Illinois, CUBIC, BBR, and Compound.

The second aim of the paper is presentation of a new TCP prototype version. Its mathematical background is both the solution to the optimal MJP state control under complete information, and the solution to the optimal MJP state filtering given the diffusion and counting observations. The performance of the proposed prototype is demonstrated on the complex of the numerical experiments.

The paper is organized as follows. Section 2 contains a detailed description of the TCP link mathematical model in terms of the controllable stochastic observation system, along with the optimization problem of data transmission through this link.

One can enhance the use of the channel resources in terms of the optimal stochastic control with incomplete information. However, this approach promises complications during its realization: starting from the proof of the optimal solution existence and concluding by bulky numerical algorithms of its realization. Hence, we propose a rather simple suboptimal solution to the problem along with its effective numerical implementation.

To develop the TCP prototype, we need a substantial mathematical framework, which is introduced in Section 3:

• Section 3.1 contains the solution to the optimal MJP control problem with instant geometric control constraints and complete information [26],
• Section 3.2 introduces a diffusion approximation for the high-frequency CPP describing the packet acknowledgment flow [27],
• Section 3.3 presents a solution to the optimal MJP state filtering problem given both counting and diffusion observations with state-dependent noise [28],
• Section 3.4 contains a numerical algorithm for the optimal filtering realization [28].

In general, the articles [26,27,28] represent a formal, detailed mathematical background of all applied inferences presented in this paper. We use it in Section 4 to develop a new congestion control algorithm as follows. At the first stage, we calculate a high-precision channel state estimate based on the available observations discretized by time. At the second stage, we apply a separation principle: the obtained filtering estimate replaces the actual MJP state during the process of the optimal control synthesis with the complete information.

The aim of Section 5 is two-fold. First, it demonstrates the potential of the proposed mathematical model to describe various versions of the TCP: classic AIMD congestion control scheme and TCP Illinois (Section 5.1), TCP CUBIC (Section 5.2), TCP Compound (Section 5.3), TCP BBR (Section 5.4).

Second, the section contains the comparison of the proposed state-based TCP with versions mentioned above: Section 5.5 highlights some details of the numerical realization of the proposed TCP version, and Section 5.6 represents the summary of the performed numerical experiments. Section 6 contains concluding remarks.

2. Problem of Optimal Data Transmission through TCP Channel

On the canonical Wiener-Poisson space with filtration $(\Omega ,\mathcal{F},\mathcal{P},\left\{{\mathcal{F}}_t\right\})$ [29,30] we consider the following controllable stochastic system, describing the TCP link functioning

$$X_t=X_0+{\int }_0^{t}A\left(u_s\right)X_{s}ds+{\alpha }_t,$$

(1)

$$Y_t={\int }_0^{t}B\left(u_s\right)X_{s}ds+{\beta }_t,$$

(2)

$$Z_t={\int }_0^{t}C\left(u_s\right)X_{s}ds+{\gamma }_t,$$

(3)

$${\left\{({\tau }_n,V_n)\right\}}_{n\in \mathbb{N}}.$$

(4)

Here the TCP link state $X_t$ is a controllable finite-state MJP with values in the set ${\mathbb{S}}^N\triangleq \{e_1,\dots ,e_n\}$ formed by unit coordinate vectors of the Euclidean space ${\mathbb{R}}^N$. The initial value $X_0$ has a known distribution $\pi$, $A\left(u\right)=\parallel A^{ij}\left(u\right){\parallel }_{i,j=\overline{1,N}}$ is a controllable transition intensity matrix and ${\alpha }_t$ is a ${\mathcal{F}}_t$-adapted martingale with the quadratic characteristic [31]

$${\langle \alpha ,\alpha \rangle }_t={\int }_0^t\left(diag\left(A\left(u_s\right)X_s\right)-A\left(u_s\right)diag\left(X_s\right)-diag\left(X_s\right)A^{\top }\left(u_s\right)\right)ds.$$

The link state is unobservable, and the complex of observations $(Y_t,Z_t,\left\{({\tau }_n,V_t)\right\})$ includes three components.

• $Y_t$ is a counting process (flow) of packet losses described by its martingale representation (2): ${\beta }_t$ is an ${\mathcal{F}}_t$-adapted martingale with the quadratic characteristic

  $${\langle \beta ,\beta \rangle }_t={\int }_0^{t}B\left(u_s\right)X_{s}ds,$$

  $B\left(u\right)\triangleq row(B^1\left(u\right),\dots ,B^N\left(u\right))$ represents the collection of the loss intensities of the flow given the conditions $X_t=e_n$, $n=\overline{1,N}$.
• $Z_t$ is a counting process (flow) of packet timeouts described by its martingale representation (3): ${\gamma }_t$ is an ${\mathcal{F}}_t$-adapted martingale with the quadratic characteristic

  $${\langle \gamma ,\gamma \rangle }_t={\int }_0^{t}C\left(u_s\right)X_{s}ds,$$

  $C\left(u\right)\triangleq row(C^1\left(u\right),\dots ,C^N\left(u\right))$ represents the collection of the timeout intensities of the flow given the conditions $X_t=e_n$, $n=\overline{1,N}$.
• $\left\{({\tau }_n,V_t)\right\}$ is a flow of successful packet acknowledgments: here ${\tau }_n$ stands for the time instant of the n-th acknowledgment arrival and $V_t$ does for the specific RTT of the n-th acknowledgment. It represents controllable compound Poisson process (CPP) with the intensity driven by the Markov state $X_t$: the predictable measure generated by $\left\{({\tau }_n,V_t)\right\}$ conditioned by the MJP state X takes the form

  $${\mu }_p(\omega ,dt,dv)=\lambda \left(u_t\right)diag\left(X_{t-}\right)\Lambda (u_t,v)dtdv.$$

  Here $\lambda \left(u_t\right)\triangleq row({\lambda }^1\left(u_t\right),\dots ,{\lambda }^N\left(u_t\right))$ is a vector-valued function with continuous positive components, its nth component represent conditional intensity of acknowledgment arrivals given $X_t=e_n$; $\Lambda (u_t,v)\triangleq col({\Lambda }^1(u_t,v),\dots ,{\Lambda }^N(u_t,v))$ is a vector-valued function with continuous components, its nth component represent conditional probability density function (pdf) with respect to v given $X_t=e_n$ for each fixed $u_t$.

All martingale terms in the processes X, Y, Z and $(\tau ,V)$ are strongly orthogonal.

The control $u_t$ represents a current size of the congestion window, i.e., portion of packets which can be instantly transmitted. The set of admissible control contains all ${\mathsf{O}}_t$-predictable processes (${\mathsf{O}}_t\triangleq \sigma \{Y_s,Z_s,({\tau }_n,V_n):s,{\tau }_n\in [0,t]\}$ stands for a natural filtration induced by all observations available up to the moment t) with the geometric constraint:

$$u_s\in \mathsf{U}\triangleq [\underline{u},\overline{u}]\subset {\mathbb{R}}_{+}\mathcal{P}-a.s.\mathrm{for} \mathrm{all}s⩾0.$$

(5)

The intensity of acknowledgment arrivals is much more than all the state transition, packet loss and timeout ones:

$$\underset{n,u}{min}{\lambda }^n\left(u\right)\gg \underset{n,u}{max}\left(\right|A^n\left(u\right)|,B^n\left(u\right),C^n\left(u\right)).$$

The performance criterion

$${\displaystyle J\left(U\right)\triangleq {\mathsf{E}}_{}\left\{\psi X_T+{\int }_0^T\left(\varphi \left(u_s\right)-u_s\xi \right)X_{s}ds\right\}\to \underset{U}{max}}$$

(6)

represents an average profit for the transmitted information, which should be maximized. Here

• $\psi \triangleq row({\psi }^1,\dots ,{\psi }^N)$ is a vector of conditional gains given the terminal state $X_T$,
• $\varphi \left(u_s\right)\triangleq row({\varphi }^1\left(u_s\right),\dots ,{\varphi }^N\left(u_s\right))$ includes strictly concave components, which represent conditional instant gains for the transmitted information given the current link state $X_s$,
• $\xi \triangleq row({\xi }^1,\dots ,{\xi }^N)$ is a vector of specific transmission expenses per information unit in each link state.

The problem under consideration is challenging. First, in general, optimal control problems of stochastic jump processes with incomplete information are rather complicated [31,32,33,34]. Their proper statement and solution depends on the answer to several auxiliary questions/problems: the martingale one [35], the one of strong solution existence and uniqueness and the one of measurable control selection (see [36] and references within). Without positive answers to the questions, we cannot use the martingale theory [35,37] to express optimal control in terms of either variation inequalities (dynamic programming equation as the preferable outcome) or stochastic maximum principle. Please note that negative answers presumes only impossibility to use the mathematical tools mentioned above. Apparently, the control problem can be modified slightly to provide its solution existence which can be found involving other still undiscovered frameworks.

Second, both the dynamic programming equation and stochastic maximum principle have forward-backward form which complicates synthesis of the optimal control in the explicit form. The authors of [36] have solved the analogous problem of the MJP state (1) control observing the flow of packet losses (2) only. The theoretical optimal solution has been characterized both via the dynamic programming equation and the maximum principle. At the same time, the authors have presented a numerical realization of the obtained result only for the case when the transition intensity matrix of the MJP is independent of the control (i.e., the state is uncontrollable), and control affects the intensity of the losses only. Despite the restrictive conditions the obtained practical results have looked rather prospective: the optimal policy has demonstrated piecewise concave nature similar to the modern versions of TCP: Illinois [6], CUBIC without probe phase [38,39], Compound [40,41] etc.

Third, the essential weak points of the optimal control implementation are its poor robustness relating to the imprecise knowledge of the control system characteristics and small perturbations of the synthesized control to its performance. This means that either control system parameters slightly misspecified towards its unknown nominal, or “instrumental errors” in control caused by imperfection of its numerical realization could nullify gain of the sophisticated optimal control in comparison with a stable suboptimal algorithm.

Fourth, the flow of packet acknowledgments has high intensity and hence leads to a high-frequency control, which is resource intensive.

Keeping in mind all arguments above we avoid the direct solution to the optimal stochastic control problem (6) of the MJP (1) state given the observations (2), (3) and (4) including the martingale problem and the ones of the solution existence and uniqueness. Instead of this we use solutions to a complex of adjacent problems and propose a suboptimal control algorithm of high performance.

3. Mathematical Background

As a basis of the proposed suboptimal control algorithm, we use the following arguments and mathematical results. We derive the algorithm basing on the following mathematical results and reasons.

• The solution to the optimal stochastic control of the MJP (1) state with the complete information does exist and can be defined as a solution to the equation of dynamic programming [26].
• The high frequency allows us to approximate the observable controlled CPP (4) by a drifting Brownian motion [42] with the parameters modulated by the MJP state [27]. We can describe the distribution of the diffusion approximation via some moment characteristics only, and this fact leads to robustness of the subsequent state filtering algorithm towards the imprecise knowledge of the specific distribution of compound Poisson process jumps.
• The conversion of high-frequency acknowledgment flow to a diffusion process gives a possibility to use the solution to the optimal MJP (1) state filtering problem given the “diffusion” and counting observations [43]. This is extension of the Wonham filter [44] to the case of the diffusion observations with state-dependent noises. Under rather mild identifiability conditions the optimal filtering estimate coincides with the exact MJP state.
• The dynamic programming equation corresponding to the control problem with complete information mentioned at item 1, represents the system of ordinary differential equations with well-developed methods of numerical solution. By contrast, the equations of the generalized Wonham filter [43] require design of special numerical procedures similar to [28].
• To complete the control synthesis, we postulate a separation principle. This means we put the state filtering estimate mentioned at items 3, 4 into the control strategy defined at item 1.

3.1. Optimal Control Strategy with Complete Information

Let us consider the controllable MJP (1) which should be optimized with respect to the optimality criterion (6) where the set $\mathcal{U}$ of all admissible controls U includes all ${\mathcal{O}}_t$-predictable processes with the geometric constraint (5).

Let us define the Bellman function $V(t,x):[0,T]\times {\mathbb{S}}^N\to \mathbb{R}$:

$$\mathbf{B}(t,x)\triangleq \underset{U\in \mathcal{U}}{sup}{\mathsf{E}}_{}\left\{\psi X_T+{\int }_t^T\left(\varphi \left(u_s\right)-u_s\xi \left(s\right)\right)X_{s}ds|X_t=x\right\}.$$

(7)

Obviously, the function $\mathbf{B}(t,x)$ can be presented in the form $\mathbf{B}(t,x)={\eta }^{\top }\left(t\right)x$, where $\eta \left(t\right)\triangleq col({\eta }^1\left(t\right),\dots ,{\eta }^N\left(t\right))=col(\mathbf{B}(t,e_1),\dots ,\mathbf{B}(t,e_N))$ is a vector-valued function.

Theorem 1. 

The assertions below are true [26].

1.

The function $\eta \left(t\right)$ is the unique solution to the Cauchy problem

$$\left\{\begin{array}{ccc}{\displaystyle {\dot{\eta }}^n\left(t\right)=\underset{u\in \mathsf{U}}{max}\left[\sum _{j=1}^N{A}^{jn}\left(u\right){\eta }^j\left(t\right)+{\xi }^n\left(u\right)\right],}\hfill & n=\overline{1,N},\hfill & 0⩽t<T,\hfill \\ {\displaystyle {\eta }^n\left(T\right)={\psi }^n,}\hfill & n=\overline{1,N}.\hfill \end{array}\right.$$

(8)

2.

There exists a Borel function ${\widehat{u}}_t\left(x\right):[0,T]\times {\mathbb{S}}^N\to \mathsf{U}$, such that

$${\displaystyle {\widehat{u}}_t\left(x\right)\in \underset{u\in \mathsf{U}}{\mathrm{Argmax}}\left[\sum _{n=1}^N\left(\sum _{j=1}^N{A}^{jn}\left(u\right){\eta }^j\left(t\right)+{\xi }^n\left(u\right)\right)x^n\right]}$$

(9)

for any $(t,s)\in [0,T]\times {\mathbb{S}}^N$.

3.

The random process ${\widehat{U}}_t\triangleq {\widehat{u}}_t\left(X_{t-}\right)$ is an optimal control strategy for the problem (1), (6).

4.

The optimal value of criterion (6) has the form ${max}_{U\in \mathsf{U}}J(U)=J(\widehat{U})={\eta }^{\top }\left(0\right)\pi$; moreover, supremum in (7) is attained for any $(t,x)\in [0,T]\times {\mathbb{S}}^N$ at the strategy $\{{\widehat{U}}_s,s\in [t,T]\}$.

The theorem establishes the base of the practical control realization. Indeed, all variants of possible optimal controls (9) can be calculated and stored in advance via solution to (8), before the control synthesis. The synthesis itself represents the selection of suitable control from the set of possible ones using the “current” MJP state $X_{t-}$.

3.2. Diffusion Approximation of High-Frequency Counting Observations

Use of the “genuine” acknowledgments flow (4) to synthesize the control leads to discontinuous one with high frequency. Its calculation may be resource intensive: each new-coming acknowledgment triggers the control recalculation algorithm. The contemporary TCP versions are exactly like this, but they are relatively simple, so not too “costly”.

Once we consider (4) discretized by time with some appropriate time increment, we can see the probability distribution of the observation increments look like mixtures of some Gaussians due to the central limit theorem for renewal-reward processes (CLTRRP). In this subsection we answer two questions. First, we determine characteristics of these mixtures. Second, we form recommendations how to choose time increment value to provide appropriate closeness of the real discretized observation distribution to the theoretical mixture above.

First, to perceive the nature of diffusion approximation, we investigate the CPPs with a fixed control $u\in \mathsf{U}$. We consider a collection of the CPPs ${\{({\tau }_n^j,V_n^j)\}}_{\genfrac{}{}{0.0pt}{}{n\in \mathbb{N},j=\overline{1,N},}{u\in \mathsf{U}}}$ with the predictable measures ${\{{\mu }_p^j(dt,dv)\}}_{\genfrac{}{}{0.0pt}{}{s>0,}{u\in \mathsf{U}}}$:

$${\mu }_p^j(dt,dv)\triangleq {\lambda }^j\left(u\right){\Lambda }^j(u,v)dsdv.$$

Probabilistically they correspond to initial CPP $\left\{({\tau }_n,V_n)\right\}$ staying in the “single mode”: $X_t\equiv e_j$ and a fixed control value $u_t\equiv u$. Each CPP generates a stochastic measure

$${\mu }^j(\omega ,dt,dv)\triangleq \sum _{n\in \mathbb{N}}{\delta }_{({\tau }_n^j\left(\omega \right),V_n^j\left(\omega \right))}(dt,dv).$$

Keeping in mind the specific form of the predictable measures ${\mu }_p^j$, we can compute the moment characteristics for one jump of the CPPs:

$$m_{\tau }^j\triangleq {\mathsf{E}}_{}\left\{{\tau }_1^j\right\}=\frac{1}{{\lambda }^j\left(u\right)},m_V^j\triangleq {\mathsf{E}}_{}\left\{V_1^j\right\}={\int }_{\mathbb{R}}v{\Lambda }^j(u,v)dv,$$

(10)

$${\sigma }_{\tau }^j\triangleq \sqrt{var\left({\tau }_1^j\right)}=\frac{1}{{\lambda }^j\left(u\right)},{\sigma }_V^j\triangleq \sqrt{var\left(V_1^j\right)}=\sqrt{{\int }_{\mathbb{R}}{v}^2{\Lambda }^j(u,v)dv-{\left(m_V^j\right)}^2},$$

$${\kappa }^j\triangleq cov({\tau }_1^j,V_1^j)=0.$$

We investigate the asymptotic behavior of the distribution of the two-dimensional random process

$${\Theta }_t^j\triangleq \left[\begin{array}{c}{\displaystyle {\int }_{[0,t]\times \mathbb{R}}{\mu }^j(ds,dv)}\\ {\displaystyle {\int }_{[0,t]\times \mathbb{R}}v{\mu }^j(ds,dv)}\end{array}\right]=\left[\begin{array}{c}{\displaystyle \sum _{n\in \mathbb{N}}\mathbf{I}(t-{\tau }_n^j)}\\ {\displaystyle \sum _{n\in \mathbb{N}}{V}_n\mathbf{I}(t-{\tau }_n^j)}\end{array}\right]$$

(11)

when $t\to \infty$. The first component represents the total number of acknowledgments received at the sender over the time interval $[0,t]$, the second component, in turn, stands for the corresponding cumulative RTT value. The author of [42] proved a version of CLTRRP:

$$\frac{1}{\sqrt{{\lambda }^j\left(u\right)t}}\left({\Theta }_t^j-\left[\begin{array}{c}{\displaystyle {\lambda }^j\left(u\right)t}\\ {\displaystyle m_V^j{\lambda }^j\left(u\right)t}\end{array}\right]\right)\stackrel{Law}{⟶}\mathcal{N}\left(\left[\begin{array}{c}0\\ 0\end{array}\right],\left[\begin{array}{cc}1& m_V^j\\ m_V^j& {(m_V^j)}^2+{({\sigma }_V^j)}^2\end{array}\right]\right)$$

(12)

as $t\to \infty$. In other words, for rather huge t

$$\frac{1}{\sqrt{t}}{\Theta }_t^j\simeq \mathcal{N}\left(\left[\begin{array}{c}{\displaystyle {\left({\lambda }^j\left(u\right)\right)}^{\frac{3}{2}}\sqrt{t}}\\ {\displaystyle m_V^j{\left({\lambda }^j\left(u\right)\right)}^{\frac{3}{2}}\sqrt{t}}\end{array}\right],\left[\begin{array}{cc}{\lambda }^j\left(u\right)& {\lambda }^j\left(u\right)m_V^j\\ {\lambda }^j\left(u\right)m_V^j& {\lambda }^j\left(u\right)\left[{(m_V^j)}^2+{({\sigma }_V^j)}^2\right]\end{array}\right]\right).$$

Let us complicate the model, mixing the CPPs ${\{({\tau }_n^j,V_n^j)\}}_{\genfrac{}{}{0.0pt}{}{n\in \mathbb{N},j=\overline{1,N},}{u\in \mathsf{U}}}$ above with probabilities $\pi =col({\pi }^1,\dots ,{\pi }^N)$

$$\left[\begin{array}{c}{\overline{\tau }}_n\\ {\overline{V}}_n\end{array}\right]=\sum _{j=1}^N{X}_0^j\left[\begin{array}{c}{\tau }_n^j\\ V_n^j\end{array}\right].$$

(13)

Here $X_0\triangleq col(X_0^1,\dots ,X_0^N)\in {\mathbb{S}}^N$ is an ${\mathcal{F}}_0$-measurable random vector, independent of ${\{({\tau }_n^j,V_n^j)\}}_{\genfrac{}{}{0.0pt}{}{n\in \mathbb{N},j=\overline{1,N},}{u\in \mathsf{U}}}$; $X_0\sim {\pi }_0$. It is easy to verify that the predictable measure generated by $\{({\overline{\tau }}_n,{\overline{V}}_n)\}$, conditioned by $X_0$, takes the form

$${\overline{\mu }}_p(\omega ,dt,dv)=\lambda \left(u_t\right)diag\left(X_0\right)\Lambda (u_t,v)dtdv.$$

Please note that the mixed CPP (13) represents a specific case of the observations (4) with “single mode” MJP X: $A\left(u\right)\equiv 0$, $X_0\sim \pi$.

Making inferences as above we can conclude that for rather huge t

$$\begin{array}{c}\frac{1}{\sqrt{t}}\left[\begin{array}{c}{\displaystyle \sum _{k\in \mathbb{N}}\mathbf{I}(t-{\overline{\tau }}_k)}\\ {\displaystyle \sum _{k\in \mathbb{N}}{\overline{V}}_k\mathbf{I}(t-{\overline{V}}_k)}\end{array}\right]\hfill \\ \simeq \sum _{j=1}^N{\pi }^j\mathcal{N}\left(\left[\begin{array}{c}{\displaystyle {\left({\lambda }^j\left(u\right)\right)}^{\frac{3}{2}}\sqrt{t}}\\ {\displaystyle m_V^j{\left({\lambda }^j\left(u\right)\right)}^{\frac{3}{2}}\sqrt{t}}\end{array}\right],\left[\begin{array}{cc}{\lambda }^j\left(u\right)& {\lambda }^j\left(u\right)m_V^j\\ {\lambda }^j\left(u\right)m_V^j& {\lambda }^j\left(u\right)\left[{(m_V^j)}^2+{({\sigma }_V^j)}^2\right]\end{array}\right]\right).\hfill \end{array}$$

(14)

Therefore, given some MJP state $X_s$ distribution (conditional or unconditional) at the time instant s and a constant control $u_q\equiv u\in \mathsf{U},q\in [s,s+h)$ we assume that the cumulative observation increment over the interval $[s,s+h)$ is distributed approximately in the following way

$$\begin{array}{c}\frac{1}{\sqrt{h}}\left[\begin{array}{c}{\displaystyle \sum _{k\in \mathbb{N}}\mathbf{I}(t-{\tau }_k)\mathbf{I}({\tau }_k-s)}\\ {\displaystyle \sum _{k\in \mathbb{N}}{V}_k\mathbf{I}(t-{\tau }_k)\mathbf{I}({\tau }_k-s)}\end{array}\right]\hfill \\ \simeq \sum _{j=1}^N{\widehat{X}}_s^j\mathcal{N}\left(\left[\begin{array}{c}{\displaystyle {\left({\lambda }^j\left(u\right)\right)}^{\frac{3}{2}}\sqrt{t}}\\ {\displaystyle m_V^j{\left({\lambda }^j\left(u\right)\right)}^{\frac{3}{2}}\sqrt{t}}\end{array}\right],\left[\begin{array}{cc}{\lambda }^j\left(u\right)& {\lambda }^j\left(u\right)m_V^j\\ {\lambda }^j\left(u\right)m_V^j& {\lambda }^j\left(u\right)\left[{(m_V^j)}^2+{({\sigma }_V^j)}^2\right]\end{array}\right]\right).\hfill \end{array}$$

(15)

By analogy with (15) for the cumulative process, corresponding to the acknowledgment flow (4)

$$Q_t\triangleq \left[\begin{array}{c}{\displaystyle \sum _{k\in \mathbb{N}}\mathbf{I}(t-{\tau }_k)}\\ {\displaystyle \sum _{k\in \mathbb{N}}{V}_k\mathbf{I}(t-{\tau }_k)}\end{array}\right]$$

(16)

we propose the following approximate diffusion model

$$Q_t={\int }_0^{t}D\left(u_s\right)X_{s}ds+{\int }_0^t\sum _{n=1}^N{e}_n^{\top }{X}_s{E}_n^{\frac{1}{2}}\left(u_s\right)d{W}_s,$$

(17)

where

$$D\left(u\right)\triangleq \left[\begin{array}{cccc}{\left({\lambda }^1\left(u\right)\right)}^{\frac{3}{2}}& {\left({\lambda }^2\left(u\right)\right)}^{\frac{3}{2}}& \dots & {\left({\lambda }^N\left(u\right)\right)}^{\frac{3}{2}}\\ m_V^1{\left({\lambda }^1\left(u\right)\right)}^{\frac{3}{2}}& m_V^2{\left({\lambda }^2\left(u\right)\right)}^{\frac{3}{2}}& \dots & m_V^N{\left({\lambda }^N\left(u\right)\right)}^{\frac{3}{2}}\end{array}\right],$$

$$E_n\left(u\right)\triangleq \left[\begin{array}{cc}{\lambda }^n\left(u\right)& {\lambda }^n\left(u\right)m_V^{u,n}\\ {\lambda }^n\left(u\right)m_V^{u,n}& {\lambda }^n\left(u\right)\left({\left(m_V^{u,n}\right)}^2+{({\sigma }_V^{u,n})}^2\right)\end{array}\right]$$

Model (17) gives a chance both to solve the MJP state filtering problem given the diffusion and counting observations and develop corresponding algorithms of the numerical solution to the filtering problem.

By contrast with weak convergence in (12), any convergence in (15) is absent. First, the right-hand side (RHS) of (15) contains the mathematical expectation which is increasing function of t. Second, we determine (15) under hypothesis that the MJP state X remains unchanged over the discretization interval: $X_q\equiv X_s,q\in [s,s+t)$. In the general case, the probability of MJP state transition increases to 1 when the interval length t increases infinitely.

Use of the time-discretized observations (4) at the first stage of the control synthesis–MJP state filtering–presumes calculation of likelihood ratios for the single Gaussian modes and their mixtures. Therefore, the filtering performance depends on both the “theoretical” pdf (15) and the closeness of real distribution of the observation increments to (15).

We form recommendations for appropriate choice of the time interval for discretization of (4). On the one hand, the length should provide the appropriate performance of the diffusion approximation (15), when there is no MJP state transitions over the time interval. On the other hand, the interval length should be small enough to guarantee small probability of those state transitions.

In the CLT the closeness of the limit distribution and the pre-limit one is described by the Berry–Esseen inequality in terms of either the uniform metric or the total variation one [45,46,47]. By contrast, we are interested in closeness of the corresponding PDFs, and the appropriated results are valid for the case of the “classic” CLT, not for CLTRRP.

We propose some heuristic technique choose the discretization interval length, basing on a performance criterion of the distribution approximation.

We refer to the “single mode” processes ${\Theta }_t^j$ and construct the processes

$${\overline{\Theta }}_h^j\triangleq {\left(\sqrt{{\Theta }_h^{j,1}}\right)}^{+}\frac{1}{{\sigma }_V^j}\left({\Theta }_h^{j,2}-m_V^j{\Theta }_h^{j,1}\right).$$

(18)

From the definition one can conclude that ${\overline{\Theta }}_h^j$ represents the normalized sum of the random number of independent equally distributed normalized random summands. We investigate closedness of its distribution to the standard Gaussian one depending on time h.

Below in the filtering algorithm we operate with various likelihood ratios calculated via the pdfs, hence we need to characterize a distance between the pre-limit pdf and its limit one. The precise distance is difficult to calculate, and we must turn to some upper bound of this quantity.

Let $\mu \left(dx\right)$ be some positive measure on $(\mathbb{R},\mathcal{B}(\mathbb{R}\left)\right)$, and there exist both the pdf $\frac{d{P}_a}{d\mu }$ of the pre-limit distribution and the limit one $\frac{d{P}_{\ell }}{d\mu }$. Then the relative approximation error takes the form

$$\Delta \left(x\right)\triangleq \frac{\left|\frac{d{P}_a}{d\mu }-\frac{d{P}_{\ell }}{d\mu }\right|}{\frac{d{P}_{\ell }}{d\mu }}\left(x\right),$$

and its average

$$\int \Delta \left(x\right)\frac{d{P}_{\ell }}{d\mu }\left(x\right)\mu \left(dx\right)=Var(P_a,P_{\ell })$$

coincides with the total variation distance (TVD) between $P_a$ and $P_{\ell }$.

We use the notation $P^j(x,h)\triangleq \mathcal{P}\left\{{\overline{\Theta }}_h^j⩽x\right\}$ for the pre-limit distribution function, $P_n^j\left(x\right)$ stands for the distribution function of the normalized sum of n independent equally distributed normalized random summands with the pdf ${\Lambda }^j\left(u\right)$, and $\mathsf{\Phi }\left(x\right)\triangleq {\int }_{-\infty }^x\frac{1}{\sqrt{2\pi }}{e}^{-\frac{z^2}{2}}dz$ does for the distribution function of the standard Gaussian random value. From the total probability formula, it follows that

$$P^j(x,h)=e^{-{\lambda }^j\left(u\right)h}\left(\mathbf{I}\left(x\right)+\sum _{n\in \mathbb{N}}\frac{{\left({\lambda }^j\left(u\right)h\right)}^n}{n!}{P}_n^j\left(x\right)\right),$$

(19)

where $\mathbf{I}\left(x\right)$ is the Heaviside function.

Proposition 1.

For ${\lambda }^j\left(u\right)h⩾\frac{3+\sqrt{13}}{2}$ an approximate upper bound of $Var(P^j,\mathsf{\Phi })$ can be written as

$$J^j\left(h\right)=e^{-{\lambda }^j\left(u\right)h}\left(2+C_1\left(2\mathsf{\Phi }(-3)+\left(\frac{1}{\sqrt{1-\frac{3}{\sqrt{{\lambda }^j\left(u\right)h}}}}+\frac{1}{\sqrt{1+\frac{3}{\sqrt{{\lambda }^j\left(u\right)h}}}}\right)\right)\right),$$

(20)

where $C_1=C_1\left({\Lambda }^j(u,·)\right)$ is some parameter.

Proof. 

From (19) and the results of [48] (Theorem 1.1) and [49] (Theorem 2.6) the following inequalities are true

$$Var(P^j,\mathsf{\Phi })⩽e^{-{\lambda }^j\left(u\right)h}\left(2+{\sum }_{n\in \mathbb{N}}\frac{{\left({\lambda }^j\left(u\right)h\right)}^2}{n!}Var(P_n^j,\mathsf{\Phi })\right)⩽e^{-{\lambda }^j\left(u\right)h}\left(2+C_1{\sum }_{n\in \mathbb{N}}\frac{{\left({\lambda }^j\left(u\right)h\right)}^2}{\sqrt{n}n!}\right),$$

(21)

where $C_1=C_1\left({\Lambda }^j(u,·)\right)$ is some parameter (see [48,49] for details).

Under the Proposition conditions the approximation of the Poisson distribution by the Gaussian one is valid

$$\begin{array}{cc}\hfill {\displaystyle \sum _{n\in \mathbb{N}}\frac{{\left({\lambda }^j\left(u\right)h\right)}^2}{\sqrt{n}n!}}& {\displaystyle \approx {\int }_1^{\infty }\frac{1}{\sqrt{x}}\frac{1}{\sqrt{2\pi {\lambda }^j\left(u\right)h}}{e}^{\frac{{(x-{\lambda }^j\left(u\right)h)}^2}{2{\lambda }^j\left(u\right)h}}dx}\hfill \\ & {\displaystyle ⩽2\mathsf{\Phi }(-3)+\underset{{\lambda }^j\left(u\right)h-3\sqrt{{\lambda }^j\left(u\right)h}}{\overset{{\lambda }^j\left(u\right)h+3\sqrt{{\lambda }^j\left(u\right)h}}{\int }}(ax+b)\frac{1}{\sqrt{x}}\frac{1}{\sqrt{2\pi {\lambda }^j\left(u\right)h}}{e}^{\frac{{(x-{\lambda }^j\left(u\right)h)}^2}{2{\lambda }^j\left(u\right)h}}dx,}\hfill \end{array}$$

(22)

where

$$\left\{\begin{array}{c}a\triangleq \frac{\sqrt{{\lambda }^j\left(u\right)h-3\sqrt{{\lambda }^j\left(u\right)h}}-\sqrt{{\lambda }^j\left(u\right)h+3\sqrt{{\lambda }^j\left(u\right)h}}}{6\sqrt{{\lambda }^j\left(u\right)h}\sqrt{{\left({\lambda }^j\left(u\right)h\right)}^2-9{\lambda }^j\left(u\right)h}},\hfill \\ b\triangleq \frac{1}{\sqrt{{\lambda }^j\left(u\right)h-3\sqrt{{\lambda }^j\left(u\right)h}}}-\frac{{\lambda }^j\left(u\right)h-3\sqrt{{\lambda }^j\left(u\right)h}-\sqrt{{\lambda }^j\left(u\right)h-9}}{6\sqrt{{\lambda }^j\left(u\right)h+3\sqrt{{\lambda }^j\left(u\right)h}}}.\hfill \end{array}\right.$$

(23)

Coefficients a and b above correspond to a piecewise linear majorant for $y\left(x\right)=\frac{1}{\sqrt{x}}$ over the interval $[1,+\infty )$ (see Figure 1).

We can calculate the last integral analytically

$$\begin{array}{c}\sum _{n\in \mathbb{N}}\frac{{\left({\lambda }^j\left(u\right)h\right)}^2}{\sqrt{n}n!}\lesssim 2\mathsf{\Phi }(-3)+(1-2\mathsf{\Phi }(-3))(a{\lambda }^j\left(u\right)h+b)⩽2\mathsf{\Phi }(-3)+a{\lambda }^j\left(u\right)h+b\hfill \\ \hfill =2\mathsf{\Phi }(-3)+\frac{1}{2\sqrt{{\lambda }^j\left(u\right)h}}\left(\frac{1}{\sqrt{1-\frac{3}{\sqrt{{\lambda }^j\left(u\right)h}}}}+\frac{1}{\sqrt{1+\frac{3}{\sqrt{{\lambda }^j\left(u\right)h}}}}\right).\end{array}$$

(24)

Using the RHS of (24) in (21) we obtain the approximate upper bound (20). This ends the sketch of the proof of the Proposition. □

To characterize the distance between the $Q_t$ (16) increment distribution and its diffusion approximation (17) we should take into account the chance of the MJP transition during the discretization interval. Let us suppose $X_t^u=e_j$, then, taking into account (20), the upper bound of $Var(P^u,\mathsf{\Phi }|X_t=e_j)$ can be obtained by the total probability formula:

$$\begin{array}{c}Var(P,\mathsf{\Phi }|X_t=e_j)⩽{\mathcal{J}}^j(u,h)\triangleq \hfill \\ \hfill \triangleq e^{(A^{jj}\left(u\right)-{\lambda }^j\left(u\right))h}\left(2+C_1\left(2\mathsf{\Phi }(-3)+\left(\frac{1}{\sqrt{1-\frac{3}{\sqrt{{\lambda }^j\left(u\right)h}}}}+\frac{1}{\sqrt{1+\frac{3}{\sqrt{{\lambda }^j\left(u\right)h}}}}\right)\right)\right)+2\left(1-e^{A^{jj}\left(u\right)h}\right).\end{array}$$

(25)

The second summand in (25) answers the chance the MJP can leave the state $e_j$ during the time interval with probability $1-e^{A^{jj}\left(u\right)h}$, and the multiplier 2 is the upper bound of the TVD for any distributions.

To take into account the statistical uncertainty of the current state $X_t^u$, we must consider the following averaged criterion:

$$\mathbf{J}(u,p_1,\dots ,p_N,h)\triangleq \sum _{j=1}^N{p}_j{\mathcal{J}}^j(u,h),$$

(26)

which describes the guaranteeing estimate of distribution distance for the case of the fixed control $u\in \mathsf{U}$ and $X_t^u\sim col(p_1,\dots ,p_N)$.

From the practical point of view, the “rational” value of the time increment h can be chosen following to the one of policies:

• Numerical analysis of the values ${\mathcal{J}}^j(u,h)$ for various $(j,u,h)$ for the choice of an appropriate value for h.
• Solution to the individual minimax problems

  $${\mathcal{J}}^j(u,h)\to \underset{h:{\lambda }^j\left(u\right)h⩾\frac{3+\sqrt{13}}{2},u\in \mathsf{U}}{min}\underset{u\in \mathsf{U}}{max},j=\overline{1,N}$$

  with subsequent choice of the maximal h from the set of the individual solutions.
• Solution to the general minimax problem

  $$\mathbf{J}(u,p_1,\dots ,p_N,h))\to \underset{h:{\lambda }^j\left(u\right)h⩾\frac{3+\sqrt{13}}{2},u\in \mathsf{U},j=\overline{1,N}}{min}\underset{\genfrac{}{}{0.0pt}{}{u\in \mathsf{U},}{(p_1,\dots ,p_N)\in \Pi }}{max}.$$

In this paper, we use the first policy as the most economical one.

3.3. Optimal Filtering of MJP State Given Counting and Diffusion Observations

In this section, we investigate MJP state (1) filtering problem given counting (2), (3) and diffusion observations (17). Without loss of generality to simplify the presentation and subsequent analysis of the solution to the MJP filtering problem we must introduce below the additional assumptions.

• The control $u_t$ represents an observable nonrandom cádlág-process.
• The noises in $Q_t$ are uniformly nondegenerate [50], i.e., ${\displaystyle \underset{\genfrac{}{}{0.0pt}{}{1⩽n⩽N,}{u\in \mathsf{U}}}{min}{E}_n\left(u\right)>\alpha I}$ for some $\alpha >0$.
• The processes $K_{ij}\left(u_t\right)\triangleq {\mathbf{I}}_{\left\{\mathbf{0}\right\}}(E_i\left(u_t\right)-E_j\left(u_t\right))$, $i,j=\overline{1,N}$ has a finite local variation (here and below $\mathbf{0}$ stands for a zero matrix of appropriate dimensionality); $K\left(u_t\right)\triangleq {\parallel K_{ij}\left(u_t\right)\parallel }_{i,j=\overline{1,N}}$ is the corresponding $N\times N$-dimensional matrix-valued function.
  The optimal filtering problem is to find a Conditional Mathematical Expectation (CME) ${\widehat{X}}_t^u\triangleq {\mathsf{E}}_{}\left\{X_t^u|O_{t+}\right\}$, where $O_t\triangleq \sigma \{Y_s,Z_s,Q_s,s\in [0,t]\}$ is a natural flow of $\sigma$-algebras generated by the observations (2), (3) and (17).

The noise intensity in the observations (17) depends on the estimated state X, and this fact prevents to apply the known results of the optimal nonlinear filtering [37]. To overcome this obstacle, we use a special transformation of available diffusion observations [28]. Here we present a sketch of this transformation.

The Ito rule gives a possibility to obtain the observable quadratic characteristics of Q:

$${\langle Q,Q\rangle }_t={\int }_0^t\sum _{n=1}^N{e}_n^{\top }{X}_s{E}_n\left(u_s\right)ds.$$

(27)

We use the normalized diffusion observations

$${\overline{Q}}_t\triangleq {\int }_0^t{\left(\frac{d{\langle Q,Q\rangle }_s}{ds}\right)}^{-\frac{1}{2}}d{Q}_s.$$

(28)

as the first block component of the transformed observations. The model of this process is the following

$${\overline{Q}}_t={\int }_0^t\overline{D}\left(u_s\right)X_{s}ds+{\overline{W}}_s,$$

(29)

where $\overline{D}\left(u\right)\triangleq {\sum }_{n=1}^N{E}_n^{-\frac{1}{2}}\left(u\right)D\left(u\right)diag{e}_n$, and ${\overline{W}}_t$ is a standard Wiener process of appropriate dimensionality.

The quadratic characteristics $\langle Q,Q\rangle$ contains essential statistical information which should be included in the estimation algorithm. This process is a linear transformation of the estimated MJP state.

It is easy to verify that

$$F(u_t,X_t)\triangleq \frac{d{\langle Q,Q\rangle }_s}{ds}{|}_{s=t+}=\sum _{n=1}^N{e}_n^{\top }{X}_t{E}_n\left(u_t\right),$$

however, result of the direct derivation is a matrix-valued function with the excess dimensionality. All its statistical information is included in the complete preimage of F:

$$F=F(u,x)\stackrel{F^{-1}}{\to }\{e_n\in {\mathbb{S}}^N:E_n\left(u\right)=F\}.$$

In [28] we explain in detail how to reduce the “rough” process F to the N-dimensional “compressed” process $H_t$, which has the model

$$H_t=L\left(u_t\right)X_t,$$

(30)

where $L\left(u_t\right)$ is an $N\times N$-dimensional matrix-valued function with cádlág components; its rows are orthogonal and contains 0 or 1 only.

One can rewrite the process $H_t$ as a cumulative sum of the jumps occurred at some nonrandom (or $O_t$-predictable) moments $\tau$ (the term $H_t^D$) and one, which accumulates jumps at the random (totally inaccessible) moments (the term $H_t^R$):

$$H_t=\underset{\triangleq H_t^D}{\underbrace{L\left(u_0\right)X_0+\sum _{\tau ⩽t}\Delta L\left(u_{\tau }\right)X_{\tau }}}+\underset{\triangleq H_t^R}{\underbrace{{\int }_0^{t}L\left(u_s\right)d{X}_s}}.$$

The process $H_t^D$ represents the second block component of the transformed diffusion observations. To obtain the third component we must express $H_t^R$ through the equivalent complex of the counting processes $G_t=col(G_t^1,\dots ,G_t^N)$:

$$G_t\triangleq {\int }_0^t(I-diag{H}_{s-})d{H}_s-H_t^D.$$

The components of the process have the following properties.

• Each component $G_t^n$ has the martingale representation

  $$G_t^n={\int }_0^t\mathbf{1}{\Gamma }_n\left(u_s\right)X_{s}ds+{\int }_0^t(1-L_n\left(u_s\right)X_{s-})L_n\left(u_s\right)d{\alpha }_s^u,$$

  (31)

  where ${\alpha }_t^u$ is the martingale from the state representation (1), $L_n\left(u\right)\triangleq e_n^{\top }L\left(u\right)$ and

  $${\Gamma }^n\left(u\right)\triangleq diag\left(L_n\left(u\right)\right){\Lambda }^{\top }\left(u\right)\left(I-diag\left(L_n\left(u\right)\right)\right).$$
• ${[G^n,G^m]}_t\equiv 0$ for any $n\ne m$, and ${\langle G^n,G^n\rangle }_t={\int }_0^t\mathbf{1}{\Gamma }_n\left(u_s\right)X_{s}ds$.

Below we present a stochastic system for the CME ${\widehat{X}}_t$ along with its properties.

Proposition 2. 

The following assertions are true.

1.

The CME $\widehat{X}$ is the unique strong solution to the stochastic system

$$\begin{array}{c}{\widehat{X}}_t={\left({\left(H_0^D\right)}^{\top }L\left(u_0\right){\pi }_0\right)}^{+}diag\left(H_0^D\right)L\left(u_0\right){\pi }_0+{\int }_0^t{\Lambda }^{\top }\left(u_s\right){\widehat{X}}_{s}ds+{\int }_0^t{\widehat{k}}_s{\overline{D}}^{\top }\left(u_s\right)d{\omega }_s\hfill \\ +\sum _{n=1}^N{\int }_0^t\left({\Gamma }_n\left(u_s\right)-\mathbf{1}{\Gamma }_n\left(u_s\right){\widehat{X}}_{s-}I\right){\widehat{X}}_{s-}{\left(\mathbf{1}{\Gamma }_n\left(u_s\right){\widehat{X}}_{s-}\right)}^{+}d{\nu }_s^n\\ +{\int }_0^t{\widehat{k}}_s{B}^{\top }\left(u_s\right){(B\left(u_s\right){\widehat{X}}_{s-})}^{+}d{\widehat{\beta }}_s+{\int }_0^t{\widehat{k}}_s{C}^{\top }\left(u_s\right){(C\left(u_s\right){\widehat{X}}_{s-})}^{+}d{\widehat{\gamma }}_s\\ \hfill +\sum _{\tau ⩽t}\left({\left({(\Delta H_{\tau }^D)}^{\top }\Delta L\left(u_{\tau }\right){\widehat{X}}_{\tau -}\right)}^{+}diag(\Delta H_{\tau }^D)L\left(u_{\tau }\right)-I\right){\widehat{X}}_{\tau -},\end{array}$$

(32)

where

$${\displaystyle \begin{array}{c}{\widehat{k}}_t\triangleq diag{\widehat{X}}_t-{\widehat{X}}_t{\left({\widehat{X}}_t\right)}^{\top }={\mathsf{E}}_{}\left\{({\widehat{X}}_t-X_t){({\widehat{X}}_t-X_t)}^{\top }|O_{t+}\right\},\hfill \\ {\displaystyle {\omega }_t\triangleq {\int }_0^t(d{\overline{Q}}_s-\overline{D}\left(u_s\right){\widehat{X}}_{s}ds),}\hfill \\ {\displaystyle {\nu }_t^n\triangleq {\int }_0^t(d{G}_s^n-\mathbf{1}{\Gamma }_n\left(u_s\right){\widehat{X}}_{s-}ds),n=\overline{1,N},}\hfill \\ {\displaystyle {\widehat{\beta }}_t\triangleq {\int }_0^t(d{Y}_s-B\left(u_s\right){\widehat{X}}_{s-}ds),}\hfill \\ {\displaystyle {\widehat{\gamma }}_t\triangleq {\int }_0^t(d{Z}_s-C\left(u_s\right){\widehat{X}}_{s-}ds).}\hfill \end{array}}$$

2.

The estimate of the maximum a posteriori probability (MAP) ${\tilde{X}}_t=e_n$: ${\displaystyle n\in \underset{1⩽m⩽N}{\mathrm{Argmax}}{e}_m^{\top }{\widehat{X}}_t}$ minimizes the ${\mathcal{L}}_1$-criterion, i.e., ${\displaystyle {\tilde{X}}_t\in \underset{{\overline{X}}_t}{\mathrm{Argmin}}{\mathsf{E}}_{}\left\{\parallel {\overline{X}}_t-X_t{\parallel }_1\right\}}$.

3.

If $E_n\left(u\right)\ne E_m\left(u\right)$ for any $n\ne m$ almost everywhere on $[0,t]$, then ${\widehat{X}}_t=X_t$ $\mathcal{P}$a.s.

The validity of items 1 and 3 in Proposition 2 can be proved by complete analogy with [28] (Theorem 1, Corollary 1), meanwhile the one of item 2 is proved in [51].

The theoretical assertions above are also meaningful from the practical point of view for subsequent design of the suboptimal control of MJP state under incomplete information. First, the CME ${\widehat{X}}_t$ represents a solution to some closed finite-dimensional stochastic system, by contrast with the general case of the optimal filtering problem [37]. Second, the paths of the CME ${\widehat{X}}_t$ usually are piecewise continuous functions with values in $\Pi$, meanwhile the MJP X state trajectories are $\mathcal{P}$-a.s. piecewise constant functions with values in ${\mathbb{S}}^N$. Therefore, we cannot directly substitute the state X by its estimate $\widehat{X}$, imposing the separation principle to this control problem. The CME $\widehat{X}$ can be easily transformed into the MAP estimate $\tilde{X}$ with the paths with the same properties as the ones of X. Assertion 2 of Proposition indicates that the proposed MAP estimate is also ${\mathcal{L}}_1$-optimal. Third, if the observation system satisfies the identifiability conditions (see Assertion 3 of Proposition) then the MJP state can be restored exactly given the indirect noisy observations. This crucial property gives a chance to reduce the initial control problem with incomplete information to the one with complete information. Obviously, any numerical realization of the filtering estimate leads to some approximation errors, nevertheless Assertion 3 allows one to hope that the small filtering errors cause acceptable control performance.

At the same time, results of Proposition 2 are difficult for the direct application. First, due to the approximation of the acknowledgment flow (4) by the diffusion model (17), the former one is valid and can be effectively applied only for the observation increments over the time interval of significant length (see Section 3.2). Second, the process $H_t$, playing the key role in the estimation, is not observable directly, and represents a result of some stochastic limit passage since it is based on the quadratic characteristic $\langle Q,Q\rangle$. Due to the boundedness from below of the diffusion observation time increment, direct calculation of $H_t$ looks impossible. In the next subsection, basing on the time-discretized diffusion observations we present a special numerical algorithm of the nonlinear filtering together with its performance characteristics.

3.4. Numerical Realization of Filtering Algorithm

To construct the numerical algorithm of the MJP state filtering given the combination of both the diffusion and counting observations we consider a time-invariant version of the observation system (1), (3), (2), (17) given the observations discretized by time with the time increment $h>0$ ($t_r\triangleq rh,r\in \mathbb{N}$):

$$X_t=X_0+{\int }_0^{t}A{X}_{s}ds+{\alpha }_t,$$

(33)

$${\mathcal{Y}}_r={\int }_{t_{r-1}}^{t_r}B{X}_{s}ds+({\beta }_{t_{r-1}}-{\beta }_{t_r}),$$

(34)

$${\mathcal{Z}}_r={\int }_{t_{r-1}}^{t_r}C{X}_{s}ds+({\gamma }_{t_{r-1}}-{\gamma }_{t_{r-1}}),$$

(35)

$${\mathcal{Q}}_r={\int }_{t_{r-1}}^{t_r}D{X}_{s}ds+{\int }_{t_{r-1}}^{t_r}\sum _{n=1}^N{e}_n^{\top }{X}_s{E}_n^{\frac{1}{2}}d{W}_s,$$

(36)

and ${\mathcal{O}}_r\triangleq \sigma \{{\mathcal{Y}}_n,{\mathcal{Z}}_n,{\mathcal{Q}}_n,n⩽r\}$ is a natural filtration generated by the discretized observations.

An assumption that coefficients $A,B,C,D$ and E are constant, is not too restrictive in practice because below we will construct the MJP control which will be constant during the time discretization intervals. Please note that the discretized observations ${\mathcal{Y}}_r,{\mathcal{Z}}_r$ and ${\mathcal{Q}}_r$ are conditionally independent given ${\mathcal{F}}_{t_r}^X\vee {\mathcal{O}}_{r-1}$ due to the properties of the Wiener-Poisson canonical space and the result of [50] (Lemma 7.5). Specifically, the distribution of ${\mathcal{Y}}_r,{\mathcal{Z}}_r$ and ${\mathcal{Q}}_r$ depends on the random vector ${\eta }_r=col({\eta }_r^1,\dots ,{\eta }_r^N)={\int }_{t_{r-1}}^{t_r}{X}_{s}ds$ is a random vector composed of the occupation times of the state X in each state $e_n$ during the interval $[t_{r-1},t_r]$. Then

• conditional distribution of ${\mathcal{Y}}_r$ given ${\mathcal{F}}_{t_r}^X\vee {\mathcal{O}}_{r-1}$ is the Poisson one with the parameter $B{\eta }_r$,
• conditional distribution of ${\mathcal{Z}}_r$ given ${\mathcal{F}}_{t_r}^X\vee {\mathcal{O}}_{r-1}$ is the Poisson one with the parameter $C{\eta }_r$,
• conditional distribution of ${\mathcal{Q}}_r$ given ${\mathcal{F}}_{t_r}^X\vee {\mathcal{O}}_{r-1}$ is the Gaussian one with the mean $D{\eta }_r$ and covariance matrix ${\sum }_{n=1}^N{\eta }_r^n{E}_n$.

Below in the presentation we use the following notations:

• $\overline{A}\triangleq {max}_{n=\overline{1,N}}\left|A_{nn}\right|$;
• $\mathcal{D}\triangleq \{u=col(u^1,\dots ,u^N):u^n⩾0,{\sum }_{n=1}^N{u}^n=h\}$ is an $(N-1)$-dimensional simplex in the space ${\mathbb{R}}^M$; $\mathcal{D}$ is a distribution support of the vector ${\upsilon }_r$;
• $\Pi \triangleq \{\pi =col({\pi }^1,\dots ,{\pi }^N):{\pi }^n⩾0,{\sum }_{n=1}^N{\pi }^n=1\}$ is a “probabilistic simplex” formed by the possible values of $\pi$;
• $N_r^X$ is a random number of the state $X_t$ transitions, occurred on the interval $[t_{r-1},t_r]$,
• ${\rho }^{k,\ell ,q}\left(du\right)$ is a conditional distribution of the vector $X_{t_r}^{\ell }{\mathbf{I}}_{\left\{q\right\}}\left(N_r^X\right){\upsilon }_r$ given $X_{t_{r-1}}=e_k$, i.e., for any $\mathcal{G}\in \mathcal{B}\left({\mathbb{R}}^M\right)$ the following equality is true:

  $${\mathsf{E}}_{}\left\{{\mathbf{I}}_{\mathcal{G}}\left({\upsilon }_r\right){\mathbf{I}}_{\left\{q\right\}}\left(N_r^X\right)X_{t_r}^{\ell }|X_{t_{r-1}}=e_k\right\}={\int }_{\mathcal{G}}{\rho }^{k,\ell ,q}\left(du\right);$$
• ${\parallel \alpha \parallel }_K^2\triangleq {\alpha }^{\top }K\alpha$, ${\langle \alpha ,\beta \rangle }_K\triangleq {\alpha }^{\top }K\beta$;
• ${\displaystyle \mathcal{N}(q,m,K)\triangleq {\left(2\pi \right)}^{-M/2}{\det }^{-1/2}Kexp\left\{-\frac{1}{2}{\parallel y-m)\parallel }_{K^{-1}}^2\right\}}$ is an M-dimensional Gaussian probability density function (pdf) with the expectation m and nondegenerate covariance matrix K;
• ${\displaystyle \mathcal{P}(n,a)\triangleq e^{-a}\frac{a^n}{n!}}$ is a Poisson distribution with the parameter a;
• ${\displaystyle {\mathsf{{\rm Y}}}^{k,j,s}(y,z,q)\triangleq {\int }_{\mathcal{D}}\mathcal{P}(y,Bv)\mathcal{P}(z,Cv)\mathcal{N}(q,Dv,\sum _{i=1}^N{v}^i{E}_i){\rho }^{k,j,s}\left(dv\right)}$.

Below is an assertion introducing the calculation algorithm of the MJP state given the discretized observations ${\widehat{\mathcal{X}}}_r\triangleq {\mathsf{E}}_{}\left\{X_{t_r}|{\mathcal{O}}_r\right\}$.

Proposition 3.

The filtering estimate ${\widehat{\mathcal{X}}}_r$ can be calculated be the following recursive algorithm

$${\widehat{\mathcal{X}}}_r^j=\frac{{\displaystyle \sum _{n_1=1}^N{\widehat{\mathcal{X}}}_r^{n_1}\sum _{s_1=0}^{\infty }{\mathsf{{\rm Y}}}^{n_1,j,s_1}({\mathcal{Y}}_r,{\mathcal{Z}}_r,{\mathcal{Q}}_r)}}{{\displaystyle \sum _{n_2,j_2=1}^N{\widehat{\mathcal{X}}}_r^{n_2}\sum _{s_2=0}^{\infty }{\mathsf{{\rm Y}}}^{n_2,j_2,s_2}({\mathcal{Y}}_r,{\mathcal{Z}}_r,{\mathcal{Q}}_r)}},j=\overline{1,N},$$

(37)

and initial condition

$${\widehat{\mathcal{X}}}_0={\pi }_0.$$

(38)

Proof of Proposition 3 can be performed similarly to [28] (Lemma 3).

To construct a numerically realizable algorithm we must restrict the sums both in the numerator and denominator of (37)

$${\overline{\mathcal{X}}}_r^j\left(S\right)=\frac{{\displaystyle \sum _{n_1=1}^N{\overline{\mathcal{X}}}_r^{n_1}\sum _{s_1=0}^S{\mathsf{{\rm Y}}}^{n_1,j,s_1}({\mathcal{Y}}_r,{\mathcal{Z}}_r,{\mathcal{Q}}_r)}}{{\displaystyle \sum _{n_2,j_2=1}^N{\overline{\mathcal{X}}}_r^{n_2}\sum _{s_2=0}^S{\mathsf{{\rm Y}}}^{n_2,j_2,s_2}({\mathcal{Y}}_r,{\mathcal{Z}}_r,{\mathcal{Q}}_r)}},j=\overline{1,N},$$

(39)

and obtain the analytical approximation of the Sth order.

We present some summands $\mathsf{{\rm Y}}$ of the low order s:

$${\displaystyle {\mathsf{{\rm Y}}}^{k,j,0}(y,z,q)={\delta }_{kj}{e}^{A_{kk}h}\mathcal{P}(y,B^{k}h)\mathcal{P}(y,C^{k}h)\mathcal{N}(q,h{D}^k,h{E}_k),}$$

$$\begin{array}{c}{\displaystyle {\mathsf{{\rm Y}}}^{k,j,1}(y,z,q)=(1-{\delta }_{kj})A_{jk}{e}^{A_{jj}h}}\hfill \\ \times {\int }_0^h{e}^{(A_{kk}-A_{jj})v}\mathcal{P}(y,B^{k}v+B^j(h-v))\mathcal{P}(z,C^{k}v+C^j(h-v))\\ \hfill \times \mathcal{N}(q,v{D}^k+(h-v)D^j,v{E}_k+(h-v)E_j)dv,\end{array}$$

$$\begin{array}{c}{\displaystyle {\mathsf{{\rm Y}}}^{k,j,2}(y,z,q)}\hfill \\ =\sum _{i:i\ne k,i\ne j}{A}_{ik}{A}_{ji}{e}^{A_{jj}h}{\int }_0^h{\int }_0^{h-v^k}{e}^{(A_{kk}-A_{ii})v^k+(A_{ii}-A_{jj})v^j}\mathcal{P}(y,B^k{v}^i+B^i{v}^i+B^j(h-v^k-v^j))\\ \times \mathcal{P}(z,C^k{v}^i+C^i{v}^i+C^j(h-v^k-v^j))\\ \hfill \times \mathcal{N}(q,v^k{D}^k+v^i{D}^i+(h-v^k-v^i)D^j,v^k{E}_k+v^i{E}_i+(h-v^k-v^i)E_j)d{v}^{i}d{v}^k,\end{array}$$

where $D^k$ is the kth column of the matrix D. Other summands are also determined by the total probability formula and have complicated form. Obviously, the integrals above cannot be calculated analytically, and we approximate them by some integral sums

$${\widetilde{\mathsf{{\rm Y}}}}^{k,j,s}(y,z,q)\triangleq \sum _{\ell }^L\mathcal{P}(y,B{v}_{\ell })\mathcal{P}(z,C{v}_{\ell })\mathcal{N}(q,D{v}_{\ell },\sum _{i=1}^N{v}_{\ell }^i{E}_i){\varrho }_{\ell }^{kj},$$

(40)

where ${\left\{v_{\ell }\right\}}_{\ell =\overline{1,L}}\subset \mathcal{D}$ is a collection of points, and ${\left\{{\varrho }_{\ell }^{kj}\right\}}_{\ell =\overline{1,L}}$ are corresponding weights, such that ${\sum }_{j=1}^N{\sum }_{\ell =1}^L{\varrho }_{\ell }^{kj}⩽1$. Therefore, we calculate the filtering estimate by the recursion

$${\tilde{\mathcal{X}}}_r^j\left(S\right)=\frac{{\displaystyle \sum _{n_1=1}^N{\tilde{\mathcal{X}}}_r^{n_1}\sum _{s_1=0}^S{\mathsf{{\rm Y}}}^{n_1,j,s_1}({\mathcal{Y}}_r,{\mathcal{Z}}_r,{\mathcal{Q}}_r)}}{{\displaystyle \sum _{n_2,j_2=1}^N{\tilde{\mathcal{X}}}_r^{n_2}\sum _{s_2=0}^S{\mathsf{{\rm Y}}}^{n_2,j_2,s_2}({\mathcal{Y}}_r,{\mathcal{Z}}_r,{\mathcal{Q}}_r)}},j=\overline{1,N},$$

(41)

and refer it as the numerical approximation of the Sth order, corresponding to a chosen numerical integration scheme.

Let us fix a time instant t, and consider the asymptotic performance of approximation (41) as $h\to 0$. The performance index is ${sup}_{\pi \in \Pi }{\mathsf{E}}_{\pi }\left\{\parallel {\tilde{X}}_r-{\widehat{X}}_r{\parallel }_1\right\}$, i.e., an average of the ${\mathcal{L}}_1$-norm of the filtering error calculated at the step r for the worst initial distribution of the MJP.

Proposition 4.

If the condition

$$\underset{k=\overline{1,N},(y,z)\in {\mathbb{Z}}_{+}^2}{max}\sum _{j=1}^N{\int }_{{\mathbb{R}}^2}\left|\sum _{s=0}^S{\widetilde{\mathsf{{\rm Y}}}}^{n_1,j,s_1}(y,z,q)-{\mathsf{{\rm Y}}}^{n_1,j,s_1}(y,z,q)\right|dq<\delta ,$$

holds, then for small enough h

$$\underset{\pi \in \Pi }{sup}{\mathsf{E}}_{\pi }\left\{\parallel {\tilde{X}}_{t/h}-{\widehat{X}}_{t/h}{\parallel }_1\right\}⩽2t\left(2\overline{A}\frac{{\left(\overline{A}h\right)}^S}{(S+1)!}+\frac{\delta }{h}\right).$$

(42)

Proof of Proposition 4 can be performed similarly to [28] (Lemma 4, Theorem 2). The first term in (42) characterizes the error of the analytical approximation: formula (39) takes into account at most S possible state transitions occurred during the time discretization interval $[t_{r-1},t_r]$. The second term in (42) describes an impact of numerical integration error to the overall performance of the filtering approximation. We can deduce that the effective choice of the integration scheme should provide the equal contribution of both summands in (42).

For the numerical study we choose the analytical approximation of the 1st order realized by the middle-point scheme:

$${\widetilde{\mathsf{{\rm Y}}}}^{k,j,0}(y,z,q)={\mathsf{{\rm Y}}}^{k,j,0}(y,z,q),$$

$${\widetilde{\mathsf{{\rm Y}}}}^{k,j,1}(y,z,q)=(1-{\delta }_{kj})A_{jk}{e}^{\frac{h}{2}(A_{jj}+A_{jj})}h\mathcal{P}(y,\frac{h}{2}(B^k+B^j))\mathcal{P}(z,\frac{h}{2}(C^k+C^j))\mathcal{N}(q,\frac{h}{2}(D^k+D^j),\frac{h}{2}(E_k+E_j)).$$

4. State-Based Modification of TCP

In this section, we describe a TCP channel mathematical model we later use for simulation of some modern TCP versions and their comparison with the state-based optimal control policy. The model we use here is in general following the one of [52]. The main distinctive characteristic of this model is the channel state allocation: we use three states to describe the wire channel condition and add one extra state to cover the issues of the wireless connection. This allocation presents a reasonable trade-off between a comprehensive connection state model taking into account all possible features (including the data flows from every channel user, the current packet distribution in all the channel hops and buffers’ queues, and signal quality in the wireless channel segment) and the feasibility of the mathematical modeling.

Thus, we suppose that the link state from a sender to a receiver is described by a controllable MJP $X_t$ (1) with four possible states:

• $e_1$ is assigned for low channel load,
• $e_2$ is for moderate load,
• $e_3$ is for wired segment congestion,
• $e_4$ is for signal fading in the wireless segment.

The intensity matrix $A\left(u\right)=\parallel A^{ij}\left(u\right){\parallel }_{i,j=\overline{1,4}}$ is defined based on the following assumptions: the link has a single bottleneck device, which remains the same during the whole transmission, this bottleneck device uses Random Early Detection (RED) queuing discipline [53], its buffer capacity is Q, and the RED threshold of guaranteed packet rejection is $W^{″}$ ($W^{″}⩽Q$).

We also assume that the wireless connection quality does not dependent on the data flow, hence the intensities $A^{·4}$ and $A^{4·}$ corresponding to the transitions from/into the state $e_4$ are independent of the control $u_s$. Furthermore, the direct transitions between the $e_1$ and $e_3$ without passing through the $e_2$ are assumed impossible, i.e., $A^{13}=A^{31}\equiv 0$.

The controllable components of $A\left(u_t\right)$ have the form

$$\begin{array}{c}{A}^{21}\left(u_t\right)=\left\{\begin{array}{cc}{A}_0^{21}+\frac{{\displaystyle C^{21}}}{{\displaystyle U_{bdp}-u_t}},\hfill & \mathrm{if}{u}_t<U_{bdp},\hfill \\ \overline{A},\hfill & \mathrm{otherwise};\hfill \end{array}\right.\hfill \\ A^{12}\left(u_t\right)=A_0^{12}+C^{12}max(U_{bdp}-u_t,0);\hfill \\ A^{32}\left(u_t\right)=\left\{\begin{array}{cc}{A}_0^{32}+\frac{{\displaystyle C^{32}}}{{\displaystyle W^{″}-u_t}},\hfill & \mathrm{if}{u}_t<W^{″},\hfill \\ \overline{A},\hfill & \mathrm{otherwise};\hfill \end{array}\right.\hfill \\ A^{23}\left(u_t\right)=A_0^{23}+C^{23}max(W^{″}-u_t,0),\hfill \end{array}$$

where $U_{bdp}$ is the control, which corresponds to the bandwidth-delay product (BDP), in other words—the maximum window size yielding throughput equal to channel bandwidth. The constant $\overline{A}$ is a level of intensity which guarantees the state transition during the forthcoming RTT.

The dependence of $A^{ji}\left(u_t\right)$ on control $u_t$ is straightforward. In the state $e_1$, the number of packets in the link is less than $U_{bdp}$; and in the state $e_2$ the “bottleneck” buffer begins to fill. The inverse proportionality of $A^{21}\left(u_t\right)$ on $u_t$ and guaranteeing intensity $\overline{A}$ provides the increasing probability of $e_1\to e_2$ transition as $u_t$ approaches to $U_{bdp}$ and guarantees the transition when the threshold $U_{bdp}$ is reached. The constant additive term $A_0^{21}$ stands for a chance of the $e_1\to e_2$ transition under low control values $u<U_{bdp}$, which are probable due to the external flows. When $u_t$ decreases to levels less than $U_{bdp}$, the probability of backward transition $e_2\to e_1$ increases linearly due to the constant flow processing rate. The transition intensities $e_2\leftrightarrows e_3$ act the same way, but with a different threshold, namely $W^{″}$.

The conditional intensities of the acknowledgment arrivals ${\lambda }^j\left(u\right)$ depend on the control u and, according to (10), are inversely proportional to the average time between the acknowledgment arrivals:

$${\lambda }^j\left(u\right)=\frac{{\displaystyle 1}}{{\displaystyle m_{\tau }^j\left(u\right)}}.$$

We assume that if no packets are lost, then during each RTT cycle, the sender receives back the acknowledgments for all the packets currently being sent into the network; hence we assume that the following relation is valid:

$$m_{\tau }^j\left(u\right)=\frac{{\displaystyle m_V^j\left(u\right)}}{{\displaystyle u}}.$$

The average RTT for each state $m_V^j\left(u\right)$ is assumed to be a sum of the following components:

• constant propagation delay, ${\delta }_0$,
• average queuing delay caused by external data, flows $m_{V,ext}^j$,
• average queuing delay caused by the data flow under control, $u·m_{V,self}^j$.

Summing up the assumptions, we have the following relation for the conditional intensity of the acknowledgment arrivals:

$${\lambda }^j\left(u\right)=\frac{{\displaystyle u}}{{\displaystyle {\delta }_0+m_{V,ext}^j+u{m}_{V,self}^j}}.$$

(43)

The counting processes for loss (2) and timeouts (3) can now be defined as thinned versions of the acknowledgment flow with following conditional intensities:

$$\begin{array}{c}{B}^j\left(u\right)=B_0^j+{\lambda }^j\left(u\right)P_l\left(u\right),\hfill \\ C^j\left(u\right)={\lambda }^j\left(u\right)P_{to}^j.\hfill \end{array}$$

(44)

Here $P_{to}^j$ denotes the conditional probabilities of a timeout in the corresponding states. For the states $e_{1,2,3}$, which are related to the wired part of the link, we assume that the only cause for a timeout is a temporary communication hardware fault; and hence the probabilities for these states are constant and equal to each other: $P_{to}^1=P_{to}^2=P_{to}^3$. In the state $e_4$, the timeouts follow the wireless carrier signal fading; hence the probability of a timeout $P_{to}^4$ is different but still independent of the control u.

The packet loss conditional probabilities, on the contrary, are the functions of the control u. If the control value is less than the RED threshold $u<W^{″}$, then

$$P_l^1\left(u\right)=P_0,P_l^2\left(u\right)=P_0+max\left(\frac{{\displaystyle U_t-W^{\prime }}}{{\displaystyle W^{″}-W^{\prime }}}(P_1-P_0),0\right),P_l^3\left(u\right)=1,P_l^4\left(U_t\right)=P_l^4,$$

where $P_0$ is the probability of a packet loss in the wired segment during its propagation through the media, $W^{\prime }$ is the lower RED threshold ($W^{\prime }<W^{″}$). If the threshold of guaranteed packet loss is exceeded, then the loss is inevitable, thus $P^j\left(u\right)=1$ for any j, if $u\ge W^{″}$.

To conclude the definition of the loss and timeout intensities, it remains to mention that the additive terms $B_0^j$ in the loss intensity $B\left(u\right)$ stand for the losses caused by the external flows.

5. Comparative Study with Modern Versions of TCP

We have completely described the observation system (1)–(4) and its parameters’ dependence on the control u. Let ${\mathcal{O}}_t\triangleq \{Y_s,Z_s,Q_s,0⩽s⩽t\}$ be the natural filtration generated by the observations available up to the moment t. Generally speaking, any ${\mathcal{O}}_t$-predictable nonnegative control $U_t$ is admissible to (1)–(4).

In this section, we present the control processes, which describe the modern versions of TCP in terms of the presented model of channel state and observations. We also present here a state-based TCP control modification, which is based on the optimal state filtering and optimal control strategy. The section will be concluded by a comparative analysis of the TCP versions’ performance.

In what follows we will assume that the constant values $U_{bdp}$, $W^{″}$, ${\delta }_0$ and $m_{V,self}^j$ are selected so as to comply with the link of $C=100$ Mbps capacity, propagation delay of ${\delta }_0=0.1$ s, bottleneck queue limit of $Q=100$ packets, and $MSS=1000$ bytes:

$$U_{bdp}=\frac{{\displaystyle {10}^{6}C{\delta }_0}}{{\displaystyle 8MSS}}=1250,W^{″}=U_{bdp}+Q=1350,m_{V,self}^j=\frac{{\displaystyle 8MSS}}{{\displaystyle {10}^{6}C}}=8·{10}^{-5}.$$

5.1. AIMD Scheme and TCP Illinois

In [52] we presented an AIMD type control $u_t$ policy, which remains the same for the present channel model:

$$\left\{\begin{array}{c}{\displaystyle u_t=u_0+{\int }_0^t{\mathbf{I}}_{[\underline{W},W_t^{th})}\left(u_{s-}\right)\frac{{\displaystyle u_{s-}}}{{\displaystyle r_{s-}}}ds+{\int }_0^t{\mathbf{I}}_{[W_t^{th},+\infty )}\left(u_{s-}\right)\frac{{\displaystyle {\alpha }_s}}{{\displaystyle r_{s-}}}ds-{\int }_0^t{\beta }_s{u}_{s-}d{Y}_s+{\int }_0^t(\underline{W}-u_{s-})d{Z}_s,}\hfill \\ {\displaystyle W_t^{th}=W_0^{th}+{\int }_0^t\left(\frac{{\displaystyle 1}}{{\displaystyle 2}}{u}_{s-}-W_{s-}^{th}\right)d{Y}_s,}\hfill \end{array}\right.$$

(45)

where

• ${\mathbf{I}}_S\left(u\right)$ is an indicator function equal to one, if $u\in S$, and zero otherwise,
• $\underline{W}$ is the minimal window size,
• $W_t^{th}$ is a threshold actuating congestion avoidance phase,
• $r_t$ is the exponential smoothing estimate of RTT,
• ${\alpha }_t$ and ${\beta }_t$ are ${\mathcal{O}}_t$-predictable coefficients of additive increase and multiplicative decrease.

The first term in (45) describes the slow start mode, the second and the third stand for the linear increase and the multiplicative decrease in the congestion avoidance phase, and the fourth provides the window rollback to the minimal value $\underline{W}$ and return to the slow start mode when a timeout event occurs.

In the case ${\alpha }_t\equiv 1$ and ${\beta }_t\equiv 0.5$ Equation (45) represents the New Reno algorithm. The Illinois concave control policy is defined by convex ${\alpha }_t$ and increasing linear ${\beta }_t$ functions of the average queuing delay $d_a=\sum _{j=1}^4(m_{V,ext}^j+u{m}_{V,self}^j)e_j^T{X}_t$:

$$\begin{array}{c}{\alpha }_t\left(d_a\right)=\left\{\begin{array}{cc}{\alpha }_{max}\hfill & \mathrm{if}{d}_a⩽d_1\hfill \\ \frac{{\displaystyle {\kappa }_1}}{{\displaystyle {\kappa }_2+d_a}}\hfill & \mathrm{otherwise},\hfill \end{array}\right.\hfill \\ {\beta }_t\left(d_a\right)=\left\{\begin{array}{cc}{\beta }_{min}\hfill & \mathrm{if}{d}_a⩽d_2\hfill \\ {\kappa }_3+{\kappa }_4{d}_a\hfill & if{d}_2<d_a<d_3\hfill \\ {\beta }_{max}\hfill & \mathrm{otherwise},\hfill \end{array}\right.\hfill \end{array}$$

(46)

The parameters ${\kappa }_i$ and $d_i$ and other details of the Illinois control scheme can be found in [6]. It should be noted that the most important parameters are the maximum and minimum additive increase and multiplicative decrease coefficients, which for the standard implementation are set to $[{\alpha }_{min},{\alpha }_{max}]=[0.3,10]$, $[{\beta }_{min},{\beta }_{max}]=[0.125,0.5]$. In Figure 2, we present the simulation results for the Illinois TCP control policy for these standard parameters. The upper plot presents the channel parameters’ dynamics, including RTT (in red), losses (black triangles), and timeouts (red crosses). The filling color indicates the channel states: white for idle, green for moderate load, red for congestion in the wired segment, and grey for the wireless segment signal fading. The lower plot shows the control dynamics and the critical thresholds: $U_{bdp}$, which corresponds to the channel bandwidth-delay product and buffer overflow low bound $U_{bdp}+W^{″}$.

One can notice that by processing only the RTT information, the algorithm succeeds in the determination of the $U_{bdp}$ and becomes much more prudent once the bottleneck buffer starts to fill. This results in long periods of relatively high transmission rates without buffer overflows and rare losses. Nevertheless, during the intervals, when the channel is idle, the control values growth speed is insufficient, which results in underuse of the channel resources and, in the end, in lower average transmission rate.

5.2. TCP CUBIC

In contrast with TCP Illinois, this version of TCP does not rely on RTT observations most of the time. Instead, it considers the control value, at which a loss occurred last time,

$$W_t^{max}=\underset{0}{\overset{t}{\int }}(u_{s-}-W_{s-}^{max})d{Y}_s,$$

as the highest network use control and tends to form a plateau in the close region to this point. To that end, it keeps counting the time since the last loss or timeout,

$$T_t^{loss}=t-\underset{0}{\overset{t}{\int }}td{Y}_s-\underset{0}{\overset{t}{\int }}td{Z}_s,$$

and sets the control according to a cubic function of $T_t^{loss}$ forming two regions: a concave region to reach the last maximum control value of $W_t^{max}$, and then a convex region of network probing, where the control growth speed becomes higher as the time without loss increases. Upon the loss event, the control is reduced according to a constant multiplicative decrease coefficient $\beta$, and when a timeout occurs, the control is reset to a minimal window size $\underline{W}$. Summing up, the TCP CUBIC control can be represented as follows:

$$u\left(t\right)=W_t^{max}+C{\left(T_t^{loss}-{\left(\frac{{\displaystyle W_t^{max}(1-\beta )}}{{\displaystyle C}}\right)}^{\frac{1}{3}}\right)}^3-\underset{0}{\overset{t}{\int }}\beta u_{s-}d{Y}_s+{\int }_0^t(\underline{W}-u_{s-})d{Z}_s,$$

(47)

where C is a constant fixed to determine the aggressiveness of control growth: with higher C values (for example, $C=4.0$), CUBIC tends to be more aggressive, which can be quite useful in high BDP networks.

In Figure 3, we present the simulation results for the TCP CUBIC control with multiplicative decrease coefficient $\beta =0.9$ and scale constant $C=4.0$. It should be noted that this simulation is based on a more precise model of the protocol described in [38] and takes into account such details as TCP-friendly region and fast convergence heuristics. These details were not reflected in Equation (47) to avoid unnecessary complications. As in the previous Figure, the upper plot presents the channel dynamics (RTT, losses, timeouts, and state), and the lower plot shows the dynamics of the control.

One can see that TCP CUBIC manages to keep the control close to the desired $U_{bdp}$ value, allowing fast recovery after losses. At the same time, the probing phase, which is symmetrical to the recovery phase, is too aggressive, and the average throughput would benefit from longer “plateau” periods. Another advantage, which must be mentioned, is the ability to adjust to dramatic changes in the media: in contrast with TCP Illinois, the CUBIC protocol keeps the control at low values throughout the whole period of wireless signal degradation, which results in fewer losses.

5.3. TCP Compound

The TCP Compound algorithm tries to benefit both from the loss-based and congestion-based approach. To that end, the authors enhance the standard AIMD congestion avoidance scheme with an additional component, which allows faster growth on an idle channel when standard AIMD control underuses the resources [40]. When the congestion is detected, the window is adjusted to avoid packet losses. To estimate the congestion, the TCP Compound scheme compares the estimated number of backlogged packets (bottleneck queue size) $d_t$ with a known threshold value $\gamma$. The estimate of the queue size is computed as follows:

$$d_t=u_t\left(1-\frac{{\displaystyle V_t^{min}}}{{\displaystyle V_t}}\right),$$

where $V_t$ is current, and $V_t^{min}$ is a minimum registered RTT value.

The entire TCP Compound control scheme can be represented by the following expression:

$$\begin{array}{c}{\displaystyle u_t=u_0+{\int }_0^t{\mathbf{I}}_{[\underline{W},W_t^{th})}\left(u_{s-}\right)\frac{{\displaystyle u_{s-}}}{{\displaystyle r_{s-}}}ds}\hfill \\ {\displaystyle +{\int }_0^t{\mathbf{I}}_{[W_t^{th},+\infty )}\left(u_{s-}\right)\left({\mathbf{I}}_{[0,\gamma )}\left(d_{s-}\right)u_t^{\kappa }\frac{{\displaystyle \alpha }}{{\displaystyle r_{s-}}}-{\mathbf{I}}_{[\gamma ,+\infty )}\left(d_{s-}\right)\zeta d_{s-}\right)ds}\hfill \\ {\displaystyle -{\int }_0^t\beta u_{s-}d{Y}_s+{\int }_0^t(\underline{W}-u_{s-})d{Z}_s,}\hfill \end{array}$$

(48)

where

• ${\mathbf{I}}_{[\underline{W},W_t^{th})}\left(u_{s-}\right)$ is the slow start indicator,
• ${\mathbf{I}}_{[0,\gamma )}\left(d_{s-}\right)$ is the congestion indicators,
• $\alpha$, $\beta$, $\kappa$, $\zeta$ are tunable protocol parameters.

In (48), the first term describes the slow start mode, the second term reflects the growth phase and correction upon congestion detection, the third stands for the multiplicative decrease, and the fourth provides the window rollback and return to the slow start mode when a timeout event occurs.

In Figure 4, we present the simulation results for the TCP Compound protocol with standard parameter values: $\alpha =\beta =0.125$, $\kappa =0.75$, $\zeta =1.0$. The backlog estimate threshold value for congestion indication is set to $\gamma =80$. The upper plot presents the channel dynamics (RTT, losses, timeouts, and state), and the lower plot shows the dynamics of the control (in black) and the estimated backlog size $d_t$ (in blue). The figure illustrates the correction of the control when the backlog size estimate reaches the threshold and high control values when the bottleneck buffer queue is assumed empty. It should be noted that TCP Compound, such as the Illinois version, fails to quickly adapt to the wireless signal degradation, demonstrating high instability and a big number of losses during this channel state.

5.4. TCP BBR

The TCP BBR algorithm is purely delay-based [54]. It is designed with the idea of maintaining the total data in the channel equal to the BDP. At this load, a connection runs with the highest throughput and lowest delay. The BDP value is estimated as a product of $RTprop$—round-trip propagation time and $BtlBw$—bottleneck bandwidth or delivery rate. An estimate for the propagation time is the minimum registered RTT over a long time:

$$RTpro{p}_t=min\left\{RT{T}_s\right\},s\in [t-W_R,t],$$

where $W_R$ typically varies from tens of seconds to minutes. To estimate the delivery rate, BBR calculates the ratio of the portion of data delivered to the time elapsed from the delivery start. Since this ratio is calculated for every acknowledgment received, it is natural to take the data “inflight” at the moment the packet was sent as a portion and the RTT of this acknowledgment as the time elapsed from the delivery start. The estimated delivery rate then is a maximum of such ratios taken over a period $W_B$ equal to 6–10 RTTs:

$$BtlB{w}_t=max\left\{\frac{{\displaystyle u(t-RT{T}_t)}}{{\displaystyle RT{T}_t}}\right\},s\in [t-W_B,t].$$

The main problem of this approach is that the propagation time and the delivery rate cannot be observed at the same time. Indeed, the bottleneck buffer must be empty to observe RTT values close to the propagation time and, to observe the capacity of the channel, it must be overfilled. This problem is solved by two modes of the steady-state regime: ProbeBW and ProbeRTT. In ProbeBW, the algorithm cycles through eight phases with the following pacing gain values: $p_t=(5/4,3/4,1,1,1,1,1,1)$. The length of each phase is equal to the current estimate of the propagation time $RTpro{p}_t$. Thus, the capacity of the channel is achieved by a periodical increase of the sending rate followed by a rollback for the queue drain. ProbRTT is turned on when the value of $RTpro{p}_t$ is not updated for a long time. In this mode, the transmission barely stops for a short time to fully drain the queue. Simulation experiments show that in the present model, the last mode is redundant since BBR manages to maintain a very precise estimate of the propagation delay spending the whole time in ProbBW mode. Plus, we excluded from consideration the Startup and Drain modes since they are usually very short.

Thus, finally, the BBR control is defined as follows:

$$u_t=RTpro{p}_t·BtlB{w}_t·p_t^{T}e[(t/RTpro{p}_t)\%8+1],$$

(49)

where $e\left[k\right]\in {\mathbb{R}}^8$ is a vector with unity on k-th place and zeros on all others, and % is the modulo operator.

In Figure 5, we present the simulation results for the TCP BBR protocol. The upper plot presents the channel dynamics (RTT, losses, timeouts, and state), and the lower plot shows the dynamics of the control (in black) and the estimate of the BDP control equal to $RTpro{p}_t·BtlB{w}_t$ (in blue). One can notice that this estimate is quite precise, nevertheless, the channel is congested almost the whole time. This means that the BBR algorithm is too aggressive for the channel at hand parameters: the bottleneck buffer size is not enough to accommodate the periodical $25\%$ sending rate increase.

5.5. State-Based TCP

To obtain the state-based TCP control strategy, the optimization problem (6) needs to be solved for some predefined gains (instantaneous and terminal) and transmission expenses.

It is natural to bind the transmission expense function $\xi ={({\xi }^1,\dots ,{\xi }^4)}^T$ with the intensity of losses, which we aim to minimize, hence set

$${\xi }^j\left(u\right)=k^j{B}^j\left(u\right),$$

(50)

where $k^1,\dots ,k^4$ are coefficients, which reflect the gravity of losses in particular channel states.

We take the same instantaneous gain, as in [36]:

$${\varphi }^j\left(u\right)=-\frac{{\displaystyle a^j}}{{\displaystyle m_V^j\left(u\right)u}}=-\frac{{\displaystyle a^j}}{{\displaystyle {\delta }_0+m_{V,ext}^j+u{m}_{V,self}^j}},$$

(51)

where $a^1,\dots ,a^4$ are coefficients, which define the utility of the traffic, depending on the channel state.

Analyzing the behavior of the TCP versions described earlier in the present paper, we may conclude that the most beneficial in terms of the throughput and losses is the state $e_2$ (moderate load). Hence, it is natural to design the state-based version with the goal of spending most of the time in this state. Terminal gains ${\psi }^j$, satisfying the condition $max\left\{{\psi }^j\right\}={\psi }^2$, would reflect this idea.

In Figure 6 (left), we present a solution to the problem (6) with transmission expenses and instantaneous gains given by (50)–(51) with $k={({10}^{-4},10,{10}^2,1)}^T$ and $a={(100,100,1,100)}^T$. The terminal gains are $\psi =-{10}^6·{(2,1,2,4)}^T$, and the right bound of the observation interval is set to a rather small value of the propagation delay $T={\delta }_0=0.1$ so that the impact of the terminal gains on the criterion would be more valuable. The controls for the states $e_1$ (idle), $e_2$ (moderate load), $e_3$ (congestion), $e_4$ (wireless signal fading) are given in grey, green, red, and black colors, respectively.

One can observe that the optimal control we obtained is almost constant. This is a very useful property in terms of the scalability of the results. Indeed, the control strategy equal to the mean of the optimal controls

$${\overline{u}}_t=\frac{1}{T}\underset{0}{\overset{T}{\int }}{u}_{s}ds,$$

(52)

does not depend on the interval, where the original optimization problem (6) was defined.

In Figure 6 (right), we present three plots, which illustrate the behavior of state occupation probabilities of the channel $X_t$ with constant controls (52) given three different initial states: $X_0=e_1$, $X_0=e_2$, $X_0=e_3$. The color scheme is the same: grey, green, red, and black lines show the occupation probabilities for respectively $e_1$, $e_2$, $e_3$, $e_4$ states. With solid lines, we show the probabilities obtained as a result of the Kolmogorov equation solution, and with dotted lines, we show the same probabilities obtained through the Monte-Carlo sampling (with 1000 trajectories). One can see that even on a bigger time interval ($T=5$ s), the goal of the state-based control is achieved: from any given initial condition, the channel manages to revert to (or maintain) the most favorable state $e_2$.

In Figure 7, we present the simulation results for the state-based control policy. The upper plot presents the channel dynamics (RTT, losses, timeouts, and state), and the lower plot shows the optimal channel estimate ${\widehat{X}}_t$ in the form of a stack plot: the height of the white/green/red/grey area at a certain point of time corresponds to the conditional probability of state idle/moderate load/congestion/wireless signal fading. This plot demonstrates that the quality of the estimates is good and that the hidden channel state may be adequately revealed based on the available information. In the lower plot of Figure 7, we also show the dynamic of the control

$${\widehat{u}}_t={\overline{u}}_t^T{\widehat{X}}_t,$$

where ${\overline{u}}_t$ is given by (52). One can see that even on a larger interval, the main property of the proposed control strategy remains: the channel spends most of the time in the state $e_2$, which results in better throughput and fewer losses.

5.6. Comparison

To compare the performance of the TCP control schemes discussed above, we use statistical modeling. The performance metrics, namely the average throughput (a measure of bandwidth usage effectiveness) and the loss percentage (a measure of predisposition to congestions, which affect other users), are calculated on samples long enough to make the variance negligible. This way is preferable in comparison with taking the average on a bunch of short-term samples since it diminishes the effect of transient phases: initial probing for available channel characteristics, which is implemented differently but is an essential part of all TCP protocol versions.

On samples of ${10}^6$ seconds, we compare the state-based control with TCP Illinois, CUBIC, Compound, and BBR versions. To make the comparison fairer, we variate, where available, the parameters of TCP control algorithms to achieve better performance. For TCP Cubic, we take three values of multiplicative decrease coefficient $\beta \in \{0.7,0.8,0.9\}$; for TCP Compound, we consider nine values of the backlog estimate threshold $\gamma \in \{10,20,30,40,50,60,70,80,90\}$. Other parameters of the protocol are the same as they were defined in Section 5.2 and Section 5.3 since they have little or negative effect on the performance.

For the state-based version described in Section 5.5, one can tune the protocol behavior by choosing different optimization criteria (6). Nevertheless, since, in our case, the optimal control is constant, instead of the variation of the coefficients of the transmission expenses (50) and instantaneous gain (51), we can directly manipulate these constant values assigned for the channel states. The experiments show that changing controls for states $e_1$, $e_2$, $e_3$, which correspond to the wired part of the transmission channel, makes the performance worse. At the same time, the variation of the control for the state $e_4$ (wireless signal fading) can bring value; hence we consider four cases: $u_t^4\in \{20,50,100,200\}$.

The simulation results are summarized in Figure 8, where we present the average throughput and loss percentage and are detailed in

Table 1. Performance metrics.

Protocol	Parameter	Throughput	% loss	${\mathit{e}}_1$	${\mathit{e}}_2$	${\mathit{e}}_3$	${\mathit{e}}_4$
Illinois		63.97	0.011	15.3%	37.7%	24.8%	22.2%
CUBIC	$\beta =0.7$	59.85	0.005	25.7%	31.2%	20.9%	22.2%
CUBIC	$\beta =0.8$	63.99	0.006	17.7%	30.8%	29.3%	22.2%
CUBIC	$\beta =0.9$	68.74	0.007	8.9%	24.2%	44.7%	22.2%
Compound	$\gamma =10$	61.81	0.021	14.3%	42.6%	20.9%	22.2%
Compound	$\gamma =20$	65.63	0.019	10.7%	43.4%	23.7%	22.2%
Compound	$\gamma =30$	66.61	0.019	9.6%	39.9%	28.3%	22.2%
Compound	$\gamma =40$	67.02	0.021	9.0%	31.1%	37.7%	22.2%
Compound	$\gamma =50$	68.08	0.022	8.5%	26.8%	42.5%	22.2%
Compound	$\gamma =60$	68.63	0.022	8.3%	24.1%	45.4%	22.2%
Compound	$\gamma =70$	68.68	0.023	8.3%	22.8%	46.7%	22.2%
Compound	$\gamma =80$	68.76	0.024	8.3%	22.9%	46.6%	22.2%
Compound	$\gamma =90$	68.77	0.024	8.3%	22.8%	46.7%	22.2%
BBR		88.65	1.219	0.7%	8.2%	68.9%	22.2%
State-based	$u_t^4=20$	76.15	0.004	1.9%	74.2%	1.7%	22.2%
State-based	$u_t^4=50$	76.68	0.007	1.8%	74.3%	1.7%	22.2%
State-based	$u_t^4=100$	77.64	0.012	1.7%	74.4%	1.7%	22.2%
State-based	$u_t^4=200$	79.29	0.022	1.6%	74.5%	1.7%	22.2%

, where one can also find the control algorithm parameters and state occupation times.

One can immediately observe the same occupation time value for the state $e_4$, which is an indirect indicator of the sufficiency of the chosen simulation sample length: since the transition to and from the state of wireless signal fading does not depend on the control values, the limit probability for the corresponding state should be the same.

The highest occupation time for the state $e_2$ of moderate channel load is demonstrated by the state-based control. In addition, it can be confirmed that this allows this control algorithm to demonstrate better performance: for the case of $u_t^4=20$, the losses are minimal, and the average throughput is second best. It should be noted that the best throughput value demonstrated by the BBR protocol is only possible at the cost of huge losses. This is a characteristic feature of this control algorithm on shallow buffers [55]: it is too aggressive for a channel with chosen characteristics, and a small buffer cannot accommodate frequent 25% speed jumps.

The last thing, which is worth mentioning, is the ability of the state-based protocol to be tuned specifically for the cases of wireless channel issues. Depending on the application, it may try to maintain the maximal possible transmission rate at a cost of huge losses, or, vice versa, drop the speed and wait for the connection to restore to the full speed.

6. Conclusions

The class of controllable Markov jump processes equipped by the stochastic analysis framework represents an effective tool for the description of a TCP governed communication connection. The hidden channel state is described by a Markov jump process with a finite-state space, characterizing both the current channel load and physical “health status”. The state equation admits both to include various types of existing congestion control algorithms (Illinois, CUBIC, Compound, BBR, etc.) and to incorporate some novelties.

The available observations represent the Markov jump processes, namely the Cox processes of the packet losses and timeouts and compound Poisson processes of the packet reception acknowledgments.

The available mathematical framework admits designing the complete technological chain of the TCP congestion control optimization, namely:

• to describe properly the congestion control problem as the stochastic control one,
• to solve the problem above in the case of complete information under the admissible controls with geometric constraints,
• to simplify the mathematical model of available observations, replacing the high-frequency packet acknowledgments flow by its diffusion limit,
• to solve the connection state filtering by the available observations and obtain high-precision state estimates,
• to design effective numerical algorithms for the filtering and control problems solution,
• to apply the separation principle and the loop of congestion control synthesis, using the connection state estimates instead of their exact values.

The result of this optimization represents the proposed state-based version of TCP. The paper contains a comparative analysis of the proposed algorithm against the other contemporary TCP versions and demonstrates its advantages.

The potential of the controllable Markov jump processes for the description of the transport and applied layer communication protocols is far from being exhausted. In perspective, one can use it both for the enhancement of the existing protocols (see, e.g., multi-path TCP [56]) and for the development of new ones (see, e.g., “TCP-free” protocols such as QUIC [57]).

In conclusion, we should also note that the mathematical potential of Markov chains/ Markov jump processes allows designing complete technological chains “mathematical model-properly formulated mathematical problem-theoretical solution-efficient numerical algorithm” to solve many applied problems of the analysis, estimation, and control in such areas as biology [58,59,60], epidemiology [61,62,63], inventory control [64], mathematical finance [65], insurance [66,67], etc.