Adv Syst Sci Appl 2023; 02:99–114 Published online at https://ijassa.ipu.ru. On Estimating the Characteristics of a Fork-Join Queueing System with Poisson Input and Exponential Service Times Anastasia V. Gorbunova1*, Alexey V. Lebedev2 1V.A. Trapeznikov Institute of Control Sciences of Russian Academy of Sciences, Moscow, Russia 2Lomonosov Moscow State University, Moscow, Russia Abstract: The paper studies the classical fork-join queueing system with M |M |1 subsystems. The analysis of this system is still relevant due to the lack of exact solutions for assessing its performance characteristics if the number of subsystems exceeds two. In addition, the fork-join system is a mathematical model of parallel or distributed computing systems that have become widespread as one of the most effective methods for processing Big Data. An approach based on graphical analysis, non-linear regression, and the use of the Nelder-Mead optimization method is proposed to estimate the mathematical expectation and dispersion of the response time of a fork-join system. As a result, the authors managed to modify the known approximations and significantly (many times) improve their approximation quality. The paper also examines the quality of the experimental data of simulation modeling used to estimate the approximation error of the obtained expressions. As a rule, this issue remains outside the scope of ongoing research in the field of this topic due to the complexity of such an analysis. And sometimes, it is due to the underestimation of the importance of this issue. The article proposes an approach to finding confidence intervals for simulation results. It provides an algorithm for their construction and also gives some recommendations. Keywords: fork-join queueing system, parallel service, parallel computing, average response time, response time variance, simulation. 1. INTRODUCTION The fork-join queueing system (QS) is a system with parallel service where the task is divided into a fixed number K of smaller components (subtasks). Then each of the subtasks enters the queue for service in one of the subsystems in the system with an infinite storage capacity and one server. After each subtasks, initially constituting one task, is serviced, the task will be considered fully serviced and will be able to leave the system. The described system is a mathematical model of many real-life and similarly functioning physical systems. The most striking example of it is data-intensive systems. It is because dividing a large task into subtasks with their subsequent parallel processing saves time [1,2]. The study of the characteristics of such systems is currently relevant due to the breadth of the so-called “Big Data” phenomenon and, accordingly, the emergence of an impressive number of services that provide services for processing and analyzing Big Data [3–5]. Studies of fork-join systems have been conducted for a relatively long time. However, the exact result for such an essential characteristic of any QS as its average response time was obtained only in the relatively simplest case. Namely, when the number of subsystems is K = 2, the input flow is Poissonian, and the service time on servers is exponential (in other ∗Corresponding author: avgorbunova@list.ru 100 A.V. GORBUNOVA, A.V. LEBEDEV words, subsytems are Markov queues of type M |M |1) [6]. If the number of subsystems is more than two, then only approximations of varying degrees of accuracy [6–9] are known. The reason for the complexity of the analysis of the described system is the presence of dependencies between the sojourn times of subtasks in subsystems. This is the main difference between the described system and the set K of parallel functioning QS. In this case, the average response time of the entire system will be determined by the longest sojourn time in the corresponding subsystem of one of its subtasks. Thus, the average response time of a fork-join QS is essentially the mathematical expectation of the maximum of K dependent random variables of sojourn times in subsystems, i.e., the K-th order statistic in the sequence of specified random variables. Therefore, elements of the theory of order statistics [10] were used along with the methods of queuing theory to approximate the average response time. However, even in this case, the obtained estimates are approximate since the theory of order statistics focuses more on the study of samples of independent random variables. In addition, if we are not talking about the Poisson input flow and exponential service times, it is quite difficult to determine the distribution of the sojourn time of the task in the QS. In more recent publications in the field of this topic, an approach is applied to the analysis of fork-join systems using machine learning methods, in particular, artificial neural networks or multiple non-linear regression [11, 12]. As a result, approximations of the mathematical expectation and higher order moments of the response time were obtained for rather complex fork-join QS structures, for example, with the Pareto distribution for service time. Details of the approach based on machine learning methods can be found in [13]. There are other approaches to deriving an estimate of the average response time of fork-join systems. Many of them are described in the review [14]. You can also highlight [1, 15–17] among recent works. This article proposes another approach to estimating the response time of a fork-join QS. This approach is based on well-known approximations for the mathematical expectation and dispersion of the sojourn time of the task in the QS. Using graphical analysis, as well as using the Nelder-Mead optimization method, the authors managed to modify the known approximations for fork-join QS with K subsystems M |M |1 and significantly (several times) improve their quality of the approximation. Analytical formulas were obtained using a limited set of experimental data; however, as shown by a numerical experiment, they are quite accurate for large values of K. In addition, the proposed approach can be used to estimate the response time of a fork-join QS not only with subsystems M |M |1 but also for more complex variants of distributions of input flow and service time. The vast majority of solutions for fork-join QS are approximate. Therefore, the question of the reliability of the experimental data used to assess the quality of the approximation of the obtained expressions comes to the fore. As a rule, numerical data can only be obtained through simulation. The values obtained using the simulation are considered reference or true. However, it is clear that in order for the simulation data to become the most accurate, a large number of realizations of the estimated random variable are required within one run of the simulation model. At the same time, determining a sufficient number of realizations and confidence intervals of the obtained estimates remains a rather difficult question due to the correlation of these data (more on this will be discussed in the corresponding section). The article proposes an approach to finding confidence intervals, provides an algorithm for their construction, as well some recommendations. The algorithm is considered in the example of a fork-join QS with subsystems of the form M |M |1; however, it can also be extended to the case with more complex QS as subsystems of the fork-join system. So, the article is organized as follows. The second section describes an approach to constructing estimates for the mathematical expectation and variance of the response time, presents the obtained estimates in the form of analytical expressions, as well as the results of a numerical experiment; the third section raises the question of the correct organization of simulation modeling, evaluates the correlation between realizations of the random value Copyright © 2023 ASSA. Adv Syst Sci Appl (2023) ASSA LATEX TEMPLATE 101 of the response time (in within the framework of the simulation), on the basis of which a confidence interval is built using the proposed algorithm. 2. EVALUATION OF THE MAIN PERFORMANCE CHARACTERISTICS OF FORK-JOIN QS In this section, we present new estimates for such important characteristics of a fork-join QS with K subsystems M |M |1 as its average response time and variance. The estimates proposed by the authors are based on the modification of known approximations obtained earlier. Experimental data, which we obtained by modeling, were needed to correct these estimates. In particular, using a simulation model written in the Python software environment, the values of the mathematical expectation and dispersion of the response time were calculated for a constant arrival rate λ = 1 and a changing service rate value µ. Namely, the load factor ρ = λ/µ took values from 0.1 to 0.9 inclusive with a step of 0.05, and the number of subsystems K varied from 3 to 20 inclusive (because for K = 2 the exact value of the average response time). The number of tests (r.v. implementations) within one run of the simulation model varied from 5 million in the case of a low system load to 10 million in the case of a high system load. Next, we describe the process of deriving new estimates based on the results of the simulation and what it leads to. 2.1. Average response time Let E[RK ] and √ Var[Rk] denote the mean and standard deviation of the response time in a fork-join QS with K subsystems of type Mλ|Mµ|1. According to the approximate Nelson- Tantawi formula for estimating the average response time [6], which gives the smallest relative error compared to other known [11] formulas, which does not exceed 5% for K ⩽ 32, we have E[RK ] ≈ E[RK ]NT = [ HK H2 + 4 11 ( 1− HK H2 ) ρ ] 12− ρ 8 1 µ− λ , (2.1) where HK = ∑K i=1 1/i is the partial sum of the harmonic series. Since the dependence of the average response time on the parameters ρ, (µ− λ) and (HK/H2 − 1) is observed in the (2.1) formula, it is natural to suggest that our correction depends on from the same settings. In particular, suppose that the improved estimate for the average response time is E[RK ] ≈ ( HK H2 − 1 ) ρ µ− λ · µ̃+ E[RK ]NT , (2.2) where we assume that the Nelson-Tantawi estimate is sharp for K = 2 and asymptotically sharp for ρ → 0, while the scale factor 1/(µ− λ) is preserved. Next, to refine µ̃, plot the modified mean response time expression as a function of ρ and (HK/H2 − 1). In this case, by the modified expression, we mean the relation µ̃ = (E[RK ]− E[RK ]NT )(µ− λ) (HK/H2 − 1)ρ , (2.3) the calculation of the right side of which is possible due to the results of simulation modeling for E[RK ] and calculations using the formula (2.1) for E[RK ]NT . In the figures 2.1 and 2.2 one can catch the linear dependence of µ̃ on ρ and HK/H2 − 1 (though not very strict). Thus, the analysis of graphs allows us to propose an expression of Copyright © 2023 ASSA. Adv Syst Sci Appl (2023) 102 A.V. GORBUNOVA, A.V. LEBEDEV 0.00 0.05 0.10 0.15 0.20 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 µ ρ Fig. 2.1. Dependence µ̃ from (2.3) on ρ. 0.00 0.05 0.10 0.15 0.20 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 µ HK/H2 – 1 Fig. 2.2. Dependence µ̃ from (2.3) on (HK/H2 − 1). the form µ̃ = µ̃(ρ,HK) ≈ C1 − C2 ( HK H2 − 1 ) + C3ρ, (2.4) i.e. E[RK ] ≈ ( HK H2 − 1 ) ρ µ− λ · ( C1 − C2 ( HK H2 − 1 ) + C3ρ ) + E[RK ]NT . (2.5) We use the Nelder-Mead optimization method [18, 19] to determine the values of the coefficients Ci, i = 1, 2, 3. The idea of using this method is to minimize the maximum value of the modulus of the relative error of approximation of the average response time by the Copyright © 2023 ASSA. Adv Syst Sci Appl (2023) ASSA LATEX TEMPLATE 103 expression (2.5) to its “true” values calculated using simulation max ∣∣∣∣ Ê[RK ]− E[RK ] E[RK ] · 100% ∣∣∣∣ → min . (2.6) As a result of minimizing the error (2.6) on the entire set of data obtained through simulation, the optimal values of the coefficients C1, C2 and C3 are determined, which best reflect the dependence (2.5) within the chosen optimization method. As a result, after applying the Nelder-Mead optimization method to (2.6) as a function of several Ci variables in the Python programming environment, we get C1 ≈ 0.087197, C2 ≈ 0.070236, C3 ≈ 0.09638. (2.7) In this case, for the above values λ = 1, ρ ∈ [0.1, 0.9] with step 0.05 and K = 3, ..., 20, the following holds MaxAPE ≈ 0.343207%, MinAPE ≈ 0.001150%, MAPE ≈ 0.142565%, while for Nelson-Tantawi formula (2.1) with the same values of λ, ρ and K the following is true MaxAPE ≈ 3.944700%, MinAPE ≈ 0.000628%, MAPE ≈ 1.322934%. Thus, our correction improves MaxAPE by a factor of 11.5 and MAPE by a factor of 9.3. The results of the approximation of the average response time using the (2.5) formula in comparison with the simulation results are shown in the figure. 0 5 10 15 20 25 0 5 10 15 20 25 F o rm u la (5 ) Simulation Fig. 2.3. Average response time. The efficiency of the expression (2.5) is not limited to the maximum value of K = 20. This formula is expected to be good for a larger number of subsystems due to the built-in logic of the proposed approach to approximating the average response time. In particular, for K = 100 and values ρ ∈ [0.1, 0.9] with a step of 0.05 we have the following approximation errors MaxAPE ≈ 2.666589%, MinAPE ≈ 0.024029%, MAPE ≈ 0.642155%. This is better than the Nelson-Tantawi formula, which gives MaxAPE ≈ 2.741596%, MinAPE ≈ 0.100878%, MAPE ≈ 1.084311%. Copyright © 2023 ASSA. Adv Syst Sci Appl (2023) 104 A.V. GORBUNOVA, A.V. LEBEDEV For the value of K = 1000, due to the saving of computational and time resources of simulation modeling, we limited ourselves to checking the quality of the approximation under low load ρ = {0.1, 0.2, 0.3, 0.4, 0.5}, since with increasing values it does not only K, but also ρ, it is necessary to significantly increase the number of tests (realizations of r.v.) within one run of the simulation model. As a result of the approximation error, the formulas (2.5) amounted to MaxAPE ≈ 1.485234%, MinAPE ≈ 0.010079%, MAPE ≈ 0.491061%, and by the Nelson-Tantawi formula, we have MaxAPE ≈ 3.389362%, MinAPE ≈ 1.090899%, MAPE ≈ 2.530636%. Note that if we add one more term to the (2.4) formula, i.e. µ̃ ≈ C1 − C2 ( HK H2 − 1 ) + C3ρ+ C4ρ ( HK H2 − 1 ) , (2.8) and recalculate the optimal coefficients after substituting (2.8) into (2.2), which turn out to be equal to C1 ≈ 0.150293, C2 ≈ 0.152088, C3 ≈ 0.012684, C4 ≈ 0.105121, then we obtain a further significant improvement in the estimates for K from 3 to 20, namely, MaxAPE ≈ 0.252453%, MinAPE ≈ 0.000495%, MAPE ≈ 0.105624%. For K = 100 we get MaxAPE ≈ 1.548824%, MinAPE ≈ 0.467368%, MAPE ≈ 1.185327%. An increase in MAPE may indicate the influence of a small systematic error, which becomes noticeable at large K. It is possible that adjusting the coefficients (while keeping the general formulas) considering the simulation for K > 20 can further improve the approximation. 2.2. Response Time Standard Deviation In [20] a formula was proposed to estimate the variance of the fork-join QS response time with K subsystems Mλ|Mµ|1 in the general case of different service rates µi, 1 ≤ i ≤ K. In the case of a single service rate µ, it takes the form Var[RK ] ≈ 2 (µ− λ)2 K∑ i=1 ( K i ) (−1)i−1 1 i2 − ( 1 µ− λ K∑ i=1 1 i )2 , (2.9) which in turn simplifies to Var[RK ] ≈ QK (µ− λ)2 , where QK = K∑ i=1 1 i2 . We note that the above formula, by its construction, proceeds from the independence of the sojourn times in subsystems, i. e., we are talking about the dispersion of the maximum of independent exponential random variables. Copyright © 2023 ASSA. Adv Syst Sci Appl (2023) ASSA LATEX TEMPLATE 105 Suppose that our estimate will depend, as well as the (2.9) expression itself, on the variance of the average response time in QS M |M |1, i.e. on 1/(µ− λ)2, and on the partial sum of the inverse square series. Let us put Var[RK ] ≈ QK − 1 (µ− λ)2 · σ̃ + 1 (µ− λ)2 , (2.10) where we assume that the estimate is sharp for K = 1. Next, to specify the expression for σ̃, we construct and analyze plots of the modified variance versus ρ, (QK − 1), and (HK − 1). Under the modified variance expression, we mean the following: σ̃ = Var[RK ](µ− λ)2 − 1 QK − 1 , (2.11) where the values of Var[RK ] were obtained by simulation. 1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2 0.35 0.4 0.45 0.5 0.55 0.6 σ QK – 1 Fig. 2.4. Dependence σ̃ from (2.11) on (QK − 1). 1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2 0.65 1.15 1.65 2.15 2.65 σ HK – 1 Fig. 2.5. Dependence σ̃ from (2.11) on (HK − 1). Copyright © 2023 ASSA. Adv Syst Sci Appl (2023) 106 A.V. GORBUNOVA, A.V. LEBEDEV 1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 σ ρ Fig. 2.6. Dependence σ̃ from (2.11) on ρ. In figures 2.4 and 2.5 at fixed values of ρ, σ̃ is observed to depend quadratically on (QK − 1) and on (HK − 1), respectively. However, the analysis of the experimental data for large K showed that the dependence on (HK − 1) is closer to the truth since the actual values of the variance grow without limit with the growth of K. Otherwise, the dispersion values would quickly stabilize since, as is known, QK has a finite limit at K → ∞ equal to π2/6. Further, in the figure 2.6 there is an explicit linear dependence on ρ with a slope depending on K, and the straight lines pass through the point (0, 1). So let us assume that σ̃ = σ̃(ρ,HK) = 1 + ρ ( C1 + C2(HK − 1) + C3(HK − 1)2 ) (2.12) and correspondingly, Var[RK ] ≈ QK − 1 (µ− λ)2 · ( 1 + ρ ( C1 + C2(HK − 1) + C3(HK − 1)2 )) + 1 (µ− λ)2 . (2.13) Then, similarly, using the Nelder-Mead method, we solve the optimization problem of minimizing the maximum approximation error of the formula (2.13) max ∣∣∣∣V̂ar[RK ]− Var[RK ] Var[RK ] · 100% ∣∣∣∣ → min . (2.14) As a result, we obtain the following values of the coefficients Ci C1 ≈ −0.113658, C2 ≈ 0.339780, C3 ≈ 0.053745. (2.15) The figure 2.7 shows the results obtained using the (2.13) formula and simulation modeling to illustrate the quality of the obtained approximation. For the standard deviation of the response time for λ = 1, ρ ∈ {0.1, 0.15, 0.20, ..., 0.90} and K = 3, ..., 20, we have MaxAPE ≈ 0.564247%, MinAPE ≈ 0.000755%, MAPE ≈ 0.188812%, while for the formula (2.9) for the same values λ, ρ and K is true MaxAPE ≈ 14.475134%, MinAPE ≈ 0.173239%, MAPE ≈ 6.055298%. Thus, our estimate improves MaxAPE by a factor of 25.7 and MAPE by a factor of 32.1. Copyright © 2023 ASSA. Adv Syst Sci Appl (2023) ASSA LATEX TEMPLATE 107 0 1 2 3 4 5 6 7 8 9 10 11 12 13 0 1 2 3 4 5 6 7 8 9 10 11 12 13 F o rm u la (1 3 ) Simulation Fig. 2.7. The standard deviation of the response time. As in the case of the mathematical expectation of the response time, the formula for the standard deviation of the response time turns out to be good not only for the maximum value of K = 20, but also for a much larger number of subsystems. Thus, for K = 100 and values ρ ∈ [0.1, 0.9] with a step of 0.05 we have the following approximation errors MaxAPE ≈ 1.211417%, MinAPE ≈ 0.094622%, MAPE ≈ 0.724898%. For the value K = 1000, for the same reasons as in the case of the formula for the mathematical expectation, we present the results of the approximation error only under conditions of low system load ρ = {0.1, 0.2, 0.3, 0.4, 0.5} MaxAPE ≈ 1.029874%, MinAPE ≈ 0.203567%, MAPE ≈ 0.592974%. As a result, a very good approximation takes place. 3. FEATURES OF FORK-JOIN QS SIMULATION The fork-join queuing system is one of the systems that is difficult to study. Exact analytical results are known only for a small number of particular cases. In most works on this topic, many approximate methods, including numerical algorithms, have been developed to analyze the performance characteristics of fork-join QS. The effectiveness of a particular method can be determined by various criteria. The most significant among them are economy and high accuracy. Efficiency is understood as the required amount of resources: the less computer memory is required, or the shorter the time of the algorithm itself, the better. Regarding the accuracy of the method, data for comparison are required here. The only possible tool for obtaining them in this context can only be simulation modeling. Thus, the assessment of the accuracy of new methods for studying fork-join QS directly depends on the accuracy of the values obtained as a result of simulation modeling. Simulation modeling is based on the Monte Carlo method, which consists of multiple reproductions of the random process of the system functioning and further statistical processing of the results obtained. As a result of the fork-join QS simulation, we can get a data set or a vector, which is a sequence of realizations of a random value of the time the request stays in the system. Further, using these data, it is necessary to construct a point estimate Copyright © 2023 ASSA. Adv Syst Sci Appl (2023) 108 A.V. GORBUNOVA, A.V. LEBEDEV of the average response time and the corresponding confidence interval. However, the main difficulty in constructing a confidence interval is that the data obtained are correlated. In particular, it is clear that the sojourn time of the n-th task in the system will depend on the sojourn time of the previous one, i.e. (n− 1)-th task, the sojourn time of which, in turn, depends on the sojourn time of the (n− 2)-th task, and so on. Therefore, the question of a sufficient number of realizations of the studied quantity during the duration of one run of the simulation model is quite relevant. Of course, in practice, it is possible to experimentally determine the required size of the sequence, gradually increasing the number of its elements with each run of the model and observing how the value of the sample average response time changes. If these changes become insignificant, the experiment is stopped, and the final result is chosen. Despite the fact that the described approach to determining the duration of the run of the simulation model is the most common, it requires a lot of time, but it does not guarantee the result due to the lack of knowledge about the confidence level of the obtained values. Therefore, this paper proposes an approach to constructing confidence intervals for the average response time of a fork-join QS under the conditions of correlation of simulation data. An estimate of the correlation between realizations of the random value of the fork-join response time of the Ri system during the duration of one run, which are at a distance of n, n = 1, 2, ... is determined by the expression r̂(n) = Ê[RiRi+n]− Ê[R]2 V̂ar[R] , (3.16) where Ê[RiRi+n] = 1 N − n N−n∑ i=1 Ri ·Ri+n, σ̂2 = V̂ar[R] = Ê[R2]− Ê[R]2, R̂ = Ê[R] = 1 N N∑ i=1 Ri, Ê[R2] = 1 N N∑ i=1 R2 i . Next, we calculate the variance of the response time estimate obtained using simulation modeling Var[R̂] = Var [ 1 N N∑ i=1 Ri ] = 1 N2 Var [ N∑ i=1 Ri ] = = 1 N2 ( N∑ i=1 Var[Ri] + 2 N−1∑ i=1 N∑ j=i+1 Cov(Ri, Rj) ) = = σ2 N ( 1 + 2 N N−1∑ i=1 N∑ j=i+1 Corr(Ri, Rj) ) = σ2 N ( 1 + 2∆N ) , where ∆N = 1 N N−1∑ i=1 N∑ j=i+1 Corr(Ri, Rj) = 1 N N−1∑ n=1 (N − n)r(n), r(n) is the correlation coefficient between a pair of elements (response times) separated from each other (according to the numbers of serviced tasks) by a distance of n. Copyright © 2023 ASSA. Adv Syst Sci Appl (2023) ASSA LATEX TEMPLATE 109 Then ∆N = 1 N N−1∑ n=1 (N − n)r(n) = N−1∑ n=1 ( 1− n N ) r(n) −−−→ N→∞ ∞∑ n=1 r(n) = ∆. The question arises of how to estimate the sum of an infinite series ∆, while in practice, we can estimate only a finite number of autocorrelations r(n). In the dissertation [21] it was proposed to use for this the exponential asymptotics r(n) ∼ ae−bn, n → ∞, (3.17) observed for the types of QS considered there. Our calculations also show good agreement with the exponential law of decreasing autocorrelations (for large ρ). For clarity, the figure 3.8 shows a graph of the dependence of sample autocorrelations calculated by the formula (3.16) on the step length n at ρ = 0.7. 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 20 40 60 80 100 r( n ) n K=2 K=10 K=20 Fig. 3.8. Correlation coefficient between the sojourn times of i-th and (i+ n)-th tasks in the system, ρ = 0.7, n = 1, 2, ...100. We can now estimate ∆ by calculating the values of r̂(n) using simulation for a finite number of steps, i.e. for n = 1, 2, ...,m, ∆ = +∞∑ n=1 r(n) = m∑ n=1 r(n) + +∞∑ n=m+1 r(n), (3.18) taking into account the fact that due to (3.17) it is true +∞∑ n=m+1 r(n) ∼ a e−b(m+1) 1− e−b , m → ∞. (3.19) Two situations are possible. We sum r̂(n) as long as it is true that r̂(n) ⩾ 2/ √ N , otherwise we calculate the values of r̂(n) statistically not significant. In the first case, the number Copyright © 2023 ASSA. Adv Syst Sci Appl (2023) 110 A.V. GORBUNOVA, A.V. LEBEDEV of obtained elements r̂(n), n = 1, 2, ...,m, may be sufficient for summation, i.e., r̂(m) < 2/ sqrtN . Otherwise, when r̂(m) ⩾ 2/ √ N , we will need to estimate summation residuals for n from (m+ 1) to +∞ from (3.19). This can be done by building a regression model. We add the resulting estimate for the tail of the sum to the sum over n from 1 to m. As a result of the described actions, we will calculate some estimate ∆̂ for ∆ from (3.18). Thus, for specific values of the parameters of the system under consideration, for example, for a certain value of the load factor ρ and the number of subsystems K: 1. choose some value m; 2. based on the simulation results for a large number of implementations N , we calculate the set of values r̂(n) using the formula (3.16), n = 1, 2, ...m; 3. sum the elements of r̂(n) until r̂(n) < 2/ √ N , n = 1, 2, ...m; 4. if r̂(m) ⩾ 2/ √ N 1) we build a regression model of the form (3.17), which is reduced to a linear one by taking the logarithm of both parts of the expression: ln r(n) = ln a− bn+ ε (ε — random error of the model); 2) estimate the parameters a and b of the constructed regression model using the least squares method; 3) calculate the estimate of the summation tail by the formula (3.19), substituting into it the values of the obtained estimates of the regression parameters a and b, as well as m; 4) calculate the estimate ∆ according to the formula (3.18), substituting the corresponding values into it ∆̂ = m∑ n=1 r̂(n) + â e−b̂(m+1) 1− e−b̂ ; 5. if r̂(n∗ + 1) < 2/ √ N , and n∗ ⩽ m, we get the estimate ∆̂ = n∗∑ n=1 r̂(n). Due to the asymptotic normality of the estimate of the mean, we can build asymptotic confidence intervals for arbitrary levels of reliability, for example, by applying the “three sigma” rule. Then with probability p = 0.997, it will be true E[R] = Ê[R]± 3σÊ[R] ≈ Ê[R]± 3σ√ N √ 1 + 2∆N ≈ Ê[R]± 3σ̂√ N √ 1 + 2∆̂. (3.20) As you can see, compared to the case of independent random variables, the number of trials N during one model run will be required √ 1 + 2∆N times more due to the existing correlation between the data. Next, we construct an estimate for ∆ for a fork-join QS with K subsystems of type M |M |1 in the form of some analytic expression depending on ρ and K. To do this, using simulation for various combinations of pairs of values K = 2, ..., 20 and ρ = 0.1, 0.2, ...0.9, according to the algorithm described above, we calculate the values ∆̂ = ∆̂(K, ρ). The values of m were chosen depending on ρ, namely m = 10, 20, 30, 40, 50, 50, 100, 250, 1000. Note that here we start the simulation with K = 2, since for this characteristic, even in this case, there are still no exact results, just as there were no approximate ones, i.e. research is completely new. Copyright © 2023 ASSA. Adv Syst Sci Appl (2023) ASSA LATEX TEMPLATE 111 -0.05 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.10 0.20 0.30 0.40 0.50 0.60 0.70 0.80 0.90 K=2 K=3 K=4 K=5 K=6 K=7 K=8 K=9 K=10 K=11 K=12 K=13 Fig. 3.9. Dependence [∆̂/(2ρ/(1− ρ2))− 1] on ρ. Now let us define a specific type of functional dependence by analyzing the results of constructing various graphs. In particular, the plot of ∆̂/(2ρ/(1− ρ2))− 1 as a function of ρ resembles a bunch of straight lines with different slopes passing through the point (0, 0), taking into account fluctuations caused by random errors (Fig. 3.9). Therefore, it remains to find a formula that describes these lines. Taking into account the approximations obtained in earlier publications for the average response time of a fork-join system in which there is a partial sum of the harmonic series, it is natural to assume that ∆ will similarly depend on K not directly but through HK . With respect to HK , a dependence close to linear was also observed. Therefore, suppose that the lines can be described by the formula C · ρ(HK − 1). Finally, we get an expression for the ∆ estimate of the following form: ∆(K, ρ) ≈ 2ρ (1− ρ)2 ( 1 + C · ρ(HK − 1) ) , (3.21) where C is some coefficient whose value is to be found. The search for the C coefficient is carried out using the Nelder-Mead optimization method. Thus, we get C ≈ 0.186722. To verify the quality of the obtained approximation, we construct the corresponding graph (Fig. 3.10) and analyze the approximation errors. Let APE = ∆̂−∆ ∆ · 100%, then we get MaxAPE ≈ 5.27449%, MinAPE ≈ 0.013941%, where the average value of the modules of relative approximation errors MAPE ≈ 0.835525%, Copyright © 2023 ASSA. Adv Syst Sci Appl (2023) 112 A.V. GORBUNOVA, A.V. LEBEDEV 0 50 100 150 200 250 0 50 100 150 200 250 F o rm u la (2 1 ) Simulation Fig. 3.10. ∆ coefficient: comparison of simulation results with the formula (3.21). which indicates the good quality of the resulting expression. In this case, we can assume that the formula (3.21) can be successfully used to construct confidence intervals not only in the region 2 ⩽ K ⩽ 20, but also for large values of K. We also note that the values of ∆ observed in the studied range of parameters reach approximately 252, which means, in this case, an increase in the required number of trials N by more than 22 times compared to the case of independent trials (to achieve the same accuracy of the estimate average response time), and this effect cannot be neglected. The algorithm presented above can be used to build confidence intervals of the form (3.20) to estimate the average response time not only in the case of subsystems of the M |M |1 type but also in more complex cases, for example, for subsystems of the type G|G|1 and others. As regards the method of constructing a functional dependence for the coefficient ∆, more laborious work is possible here on the selection of the type of dependence in each specific case, i.e., for specific distributions of input flow and service times. In particular, dependence on HK naturally arises for an exponential distribution, and for example, in the case of distributions with heavy (power-law) tails (for example, as in [12]), dependence on some degree of K will be natural. However, the general approach to constructing an estimate for ∆, described by the example of deriving the (3.21) expression, can be applied to other types of QS fork-join subsystems. 4. CONCLUSION The paper presents an approach to obtaining more accurate than known estimates for the expectation and response time variance of a fork-join system with K subsystems M |M |1. As shown by the numerical experiment, the obtained analytical formulas are valid for a wide range of values of the parameters included in them. The considered approach can be used to construct good quality estimates in more complex cases. For example, when the subsystems are queuing systems of the form G|G|1, which can be the subject of further research. In addition, the article describes the essential aspects of fork-join QS simulation and the evaluation of its results, which is an integral part of the organization of a numerical experiment to check the quality of approximation of any approximation. An algorithm for Copyright © 2023 ASSA. Adv Syst Sci Appl (2023) ASSA LATEX TEMPLATE 113 constructing confidence intervals for sample estimates obtained by simulating the functioning of a fork-join QS is presented, and some recommendations are also given. REFERENCES 1. Nguyen, M., Alesawi, S., Li, N., Che, H. & Jiang, H. (2020). A Black-Box Fork- Join Latency Prediction Model for Data-Intensive Applications, IEEE Transactions on Parallel and Distributed Systems, 31(9), 1983-2000. 2. Peng, P., Soljanin, E. & Whiting, P. (2022). Diversity/Parallelism Trade-Off in Distributed Systems With Redundancy IEEE Transactions on Information Theory, 68(2), 1279-1295. 3. Samyukta Sethuraman (2022). Analysis of Fork-Join Systems: Network of Queues with Precedence Constraints. Boca-Raton, FL: CRC Press. 4. de Oliveira, D. C. M., Liu, J. & Pacitti, E. (2019). Data-Intensive Workflow Management: For Clouds and Data-Intensive and Scalable Computing Environments. San Rafael, CA: Morgan & Claypool Publishers. 5. Khalid, M. & Yousaf, M. M. (2021). A comparative analysis of big data frameworks: An adoption perspective, Applied Sciences, 11(22), 11033. 6. Nelson, R. & Tantawi, A. N. (1988). Approximate analysis of fork/join synchronization in parallel queues, IEEE Transactions on Computers, 37(6), 739–743. 7. Varki, E., Merchant, A. & Chen, H. (2002). The M/M/1 fork-join queue with variable subtasks, [Online]. Available: http://www.cs.unh.edu/ varki/publication/2002- nov-open.pdf. 8. Varma, S. & Makowski, A. M. (1994). Interpolation Approximations for Symmetric Fork-Join Queues, Performance Evaluation, 20, 245–265. 9. Kemper, B. & Mandjes, M. (2012). Mean sojourn time in two-queue fork-join systems: Bounds and approximations, OR Spectrum, 34, 723–742. 10. David, H. A. & Nadaraja, H. N. (2003). Order Statistics. Hoboken, NJ: John Wiley & Sons. 11. Gorbunova, A. V. & Vishnevsky, V. M. (2020). Estimating the response time of a cloud computing system with the help of neural networks, Advances in Systems Science and Applications, 20(3), 105–112. 12. Gorbunova, A. V. & Lebedev, A. V. (2022). Response Time Estimate for a Fork-join System with Pareto Distributed Service Time as a Model of a Cloud Computing System Using Neural Networks, Communications in Computer and Information Science, 1552, 318–332. 13. Vishnevsky, V. M. & Gorbunova, A. V. (2022). Application of Machine Learning Methods to Solving Problems of Queuing Theory, Communications in Computer and Information Science, 1605, 304–316. 14. Thomasian, A. (2014). Analysis of fork/join and related queueing systems, ACM Computing Surveys (CSUR), 47(2), 17:1–17:71. 15. Qiu, Z., Perez, J. F. & Harrison, P. G. (2015). Beyond the mean in fork-join queues: Efficient approximation for response-time tails, Performance Evaluation, 91, 99–116. 16. Wang, W., Harchol-Balter, M., Jiang, H., Scheller-Wolf, A. & Srikant, R. (2019). Delay asymptotics and bounds for multitask parallel jobs, Queueing Syst., 91, 207–239. 17. Gorbunova, A. V. & Lebedev, A. V. (2020). Bivariate Distributions of Maximum Remaining Service Times in Fork-Join Infinite-Server Queues, Problems of Information Transmission, 56(1), 73–90. 18. Nelder, J. A. & Mead, R. (1965). A Simplex Method for Function Minimization, Computer Journal, 7, 308–313. 19. Kochenderfer, M. J. & Wheeler, T. A. (2019). Algorithms for Optimization. Cambridge, UK: MIT Press. Copyright © 2023 ASSA. Adv Syst Sci Appl (2023) 114 A.V. GORBUNOVA, A.V. LEBEDEV 20. Gorbunova, A., Zaryadov, I, Matushenko, S., Sopin, E. (2016). The Estimation of Probability Characteristics of Cloud Computing Systems with Splitting of Requests, Communications in Computer and Information Science, 678, 418–429. 21. Zaharenkova, E. (2019). Analitiko-statisticheskie metody rascheta i optimizacii sistem i setej massovogo obsluzhivanija so stepennymi hvostami raspredelenij [Analytical and statistical methods for calculating and optimizing queuing systems and networks with power tails of distributions]. Ph. D. thesis, Omsk, [in Russian]. Copyright © 2023 ASSA. Adv Syst Sci Appl (2023) Introduction Evaluation of the main performance characteristics of fork-join QS Average response time Response Time Standard Deviation Features of fork-join QS simulation Conclusion