## 2 Methodology of Parallel System Modeling
### System of Heterogeneous Parallel Servers  
Next, we considered a multiclass or heterogeneous server system with an arrival rate \(\lambda\) with a Poisson distribution. We considered a three-parallel-server system with service rates \(\mu_{1}\), \(\mu_{2}\), and \(\mu_{3}\), also with Poisson distributions. Parts or customers were serviced in entry-ordered discipline. The state diagram of the system, namely, model A, is given in Figure 1.  
Each server can only allocate one unit and the three-tuple indicates which server is in service (1) or empty (0), read from right to left, like a binary number. Thus, (110) denotes  
Figure 1: Model A transition rate diagram of 3 parallel heterogeneous entry-ordered servers. In the case where \(\mu_{1}=\mu_{2}=\mu_{3}\), it can be directly compared with the Erlang first model for \(k=3\) using (2). The throughput was calculated for model A by applying the PASTA property, and thus, \(\lambda_{\textit{eff}}=\lambda(1-p_{111})=\lambda(1-p_{f})\), where \(p_{f}\) is the probability of the maximum size or full system. There were significant differences between the experimentally simulated throughput and that calculated by model A, which mostly overestimated the throughput when the arrival rates were high (saturation). Note that in this conventional transition diagram that enumerates the states of the servers, the probability of the full system with all busy servers estimated the probability of loss. The transition to the full system (111) happened at the rate \(\lambda\) (birth) from the system of size 2--with states (110), (011), or (01)--and the system leaves (death) the full state from size 3 to 2 at a rate \(\Sigma\mu\). We also considered the homogeneous system M/M/m/m with \(m=3\), and thus, we could apply the Erlang-B Formula (2) for the calculation of the effective throughput \(\lambda_{\textit{eff}}=\lambda(1-p_{111})=\lambda(1-B(3,r))\); see Table 1. For the Erlang model, the results were completely in agreement with those from the simulation, and thus, the identical set of servers of M/M/m/m behaved with no difference from an ordered entry system, but since the servers were indistinguishable, the order became irrelevant and there would not be any difference when serving according to ordered entry, random entry, or another service discipline [31]. Next, we considered a new formulation of the transition diagram, namely, model B. Different from model A, when the system is full, the lost load is redirected to a virtual server with a queue of infinite capacity, herein called the loss server, where arrivals enter when the real servers of the system are occupied (Figure 2). When arrivals are lost, the busy system enters the state >\(111=b\). The transition matrix is given by (4), and the mean number of customers L in the system or the utilization U was calculated using (5). Note that when the system falls into the loss server, state >\((111)=b\), the three-server subsystem is occupied and the calculation of L in the servers by (5) is taken into account. In the transition diagram, the transition (birth) from the busy system (111) into the loss state happened at the rate \(\lambda\); therefore, the arrival goes into the loss server when servers are busy and the system leaves this loss state from size 3 to 2 at the rate \(\Sigma\mu\) as soon any server is available to serve, as represented. The results of the simulation, as shown in Table 1, were obtained with Arena software by Rockwell Inc. using \(\lambda=\mu_{i}=1\), a simulation run of 30 replications of 1000 h, 1000 h of warm-up for the stationary regime, and where the throughput is expressed as dimensionless by dividing by \(\Sigma\mu_{i}\). The effective throughput rate was calculated using \(\lambda_{\textit{eff}}=\lambda(1-p_{\textit{busy}})\), where \(p_{\textit{busy}}\) is the percentage of time the three-server system is busy.