This article deals with the queueing system with unlimited number of facilities. Arrival process is a Markov renewal process. An arrival customer is the customer of the first type with the probability p(1) and the customer of the second type with the probability p(2). Every customer comes into any vacant server where it is served during a stochastic time period distributed according to the exponential law with the parameter mu(1) for customers of the first type and with the parameter mu(2) for customers of the second type. The problem consists to study the process, which characterized by a number of occupied services. The system of Kolmogorov's differential equations is derived. Using characteristic function, the main equation of research is obtained: partial derivative H(u, w, z)/partial derivative z + partial derivative H(u, w, 0)/partial derivative z[p(1)e(ju)PD(z) + p(2)e(jw)PD(z) - I] - mu(1)j(e(-ju)-1)partial derivative H(u, w, z)/partial derivative u - mu(2)j(e(-jw)-1)partial derivative H(u, w, z)/partial derivative w = 0. In particular, the first and the second order initial moments of the number of busy servers of different types are obtained. Furthermore, we found the expression for the correlation coefficient between the number of the different types of busy servers. The resulting correlation coefficient indicates that the number of busy servers of different type in the system are dependent. The systems under considerations are studied using asymptotic analysis. Namely, the expressions for the asymptotic of the first and second orders are received for the characteristic functions of the busy servers of any type in the system MR/GI/8 with the heterogeneous service. Therefore, we make conclusion that the system under study (with two types of customers) cannot be reduced in two separate systems (each with one type of customers) The method of asymptotic analysis in a condition of equivalent growing service time is offered. The asymptotic characteristic functions of the first and second orders are derived. It is shown that the asymptotic characteristic function of the second order by two-dimensional probability distribution of the number of the occupied devices of each type in the system has the Gaussian distribution: h(2)(u(1), u(2)) approximate to exp{lambda j(p(1)u(1)/mu(1) + p(2) u(2)/mu(2)) + +j(2)[lambda/2(p(1) u(1)(2)/mu(1) + p(2) u(2)(2)/mu(2)) + f'(0) E((p(1)u(2))(2)/2 mu(1) + (p(2)u(2))(2)/2 mu(2)+ 2 p(1)p(2) u(1)u(2)/mu(1)+mu(2))]}. The numerical analysis of the convergence of the exact and asymptotic algorithms computing the main probabilistic characteristics of the system is carried out.