## 4 Gaussian Approximation of the Probability Distribution of the Total Resource Amounts
\begin{array}{c}\mathbf{g}\mathbf{Q}=\mathbf{r}(\lambda\mathbf{I}-\mathbf{ \Lambda}),\\ \mathbf{g}\mathbf{e}=0;\end{array} \right.\]
* _For the renewal process:_ \[\lambda=\left(\int\limits_{0}^{\infty}(1-A(x))dx\right)^{-1},\quad\kappa= \lambda^{3}\left(\sigma^{2}-a^{2}\right),\] (15) _where_ \(a\) _and_ \(\sigma^{2}\) _are the expectation and variance of a random variable given by the distribution function_ \(A(z),\,\mathbf{I}= \begin{pmatrix}1&0&\ldots&0\\ 0&1&\ldots&0\\ \vdots&\vdots&\ddots&\vdots\\ 0&0&\ldots&1\end{pmatrix} _{K\times K}.\)__  
Proof of Theorem 1.: Step 1. For a system with MMPP. Let us perform substitutions  
\[\varepsilon=\frac{1}{N},\,u_{i}=\varepsilon x_{i},\,i=1,\ldots,n,\,\mathbf{h} (u_{1},\ldots,u_{n},t)=\mathbf{f}_{1}(x_{1},\ldots,x_{n},t,\varepsilon). \tag{16}\]  
Then problem (4) takes the form  
\[ \begin{split} e\frac{\partial\mathbf{f}_{1}(x_{1},\ldots,x_{n},t,\varepsilon)}{\partial t}=\mathbf{f}_{1}(x_{1},\ldots,x_{n},t,\varepsilon) \Bigg{[}\mathbf{Q}+\sum_{i=1}^{n}S_{i}(t)\{G^{*}(\varepsilon x_{i})-1\}\mathbf{ \Lambda}\Bigg{]},\\ \mathbf{f}_{1}(x_{1},\ldots,x_{n},t_{0},\varepsilon)=\mathbf{r}. \end{split}  \tag{17}\]Let \(\varepsilon\to 0\), then problem (17) has the form  
\[ \begin{split}&\mathbf{f}_{1}(x_{1},\dots,x_{n},t)\mathbf{Q}=0,\\ &\mathbf{f}_{1}(x_{1},\dots,x_{n},t_{0})=\mathbf{r},\end{split} \]  
and the asymptotic solution \(\mathbf{f}_{1}(x_{1},\dots,x_{n},t)\) of (17) should be searched in the following form:  
\[ \begin{split}&\mathbf{f}_{1}(x_{1},\dots,x_{n},t)=\mathbf{r} \Phi_{1}(x_{1},\dots,x_{n},t),\\ &\Phi_{1}(x_{1},\dots,x_{n},t_{0})=1,\end{split}  \tag{18}\]  
where \(\Phi_{1}(x_{1},\dots,x_{n},t)\) is some scalar function. We then substitute (18) into (17) and multiply the resulting equation by \(\mathbf{e}\). Then, given that \(\mathbf{r}\Lambda\mathbf{e}=\lambda\), we get  
\[\varepsilon\frac{\partial\Phi_{1}(x_{1},\dots,x_{n},t)}{\partial t}=\lambda \Phi_{1}(x_{1},\dots,x_{n},t)\sum\limits_{i=1}^{n}S_{i}(t)\{G^{*}(\varepsilon x _{i})-1\}. \tag{19}\]  
We use the expansion of the exponent in the Taylor series for \(G^{*}(\varepsilon x_{i})\):  
\[ \begin{split} G^{*}(\varepsilon x_{i})=\int\limits_{0}^{ \infty}e^{j\varepsilon x_{i}y}& dG_{i}(y)=\int\limits_{0}^{ \infty}\Bigl{(}1+j\varepsilon x_{i}y+O(\varepsilon^{2})\Bigr{)}dG_{i}(y)=\\ &=1+j\varepsilon x_{i}a_{i}+O(\varepsilon^{2}),\end{split}  \tag{20}\]  
where \(a_{i}=\int\limits_{0}^{\infty}ydG_{i}(y)\). Substituting this into (19) and then divide the resulting equation by \(\varepsilon\left(\varepsilon\to 0\right)\), we derive  
\[\frac{\partial\Phi_{1}(x_{1},\dots,x_{n},t)}{\partial t}=j\lambda\Phi_{1}(x_{ 1},\dots,x_{n},t)\sum\limits_{i=1}^{n}S_{i}(t)x_{i}a_{i}. \tag{21}\]  
The solution of the differential Equation (21) is  
\[\Phi_{1}(x_{1},\dots,x_{n},t)=\exp\bigg{\{}j\lambda\sum\limits_{i=1}^{n}x_{i}a _{i}\int\limits_{0}^{t}S_{i}(\tau)d\tau\bigg{\}}. \tag{22}\]  
By performing back-substitutions and putting \(t_{0}=-\infty\), \(T=t=0\), we can write the first-order asymptotic characteristic function \(h_{1}(u_{1},\dots,u_{n})\) in the form  
\[ \begin{split}& h_{1}(u_{1},\dots,u_{n})=\mathbf{h}_{1}(u_{1}, \dots,u_{n},0)\mathbf{e}\approx\mathbf{f}_{1}(x_{1},\dots,x_{n},0)\mathbf{e}= \\ &=\mathbf{r}\exp\bigg{\{}j\lambda\sum\limits_{i=1}^{n}x_{i}a_{i} \int\limits_{-\infty}^{0}S_{i}(\tau)d\tau\bigg{\}}\mathbf{e}=\exp\bigg{\{}jN \lambda\sum\limits_{i=1}^{n}u_{i}a_{i}b_{i}\bigg{\}}.\end{split}  \tag{23}\]  
To construct the second-order asymptotic characteristic function \(\mathbf{h}_{2}(u_{1},\dots,u_{n},t)\), we represent the function \(\mathbf{h}(u_{1},\dots,u_{n},t)\) in the following form:  
\[\mathbf{h}(u_{1},\dots,u_{n},t)=\mathbf{h}_{2}(u_{1},\dots,u_{n},t)\exp\bigg{\{} jN\lambda\sum\limits_{i=1}^{n}u_{i}a_{i}b_{i}\bigg{\}}. \tag{24}\]  
Substituting this expression into (4), we obtain the equation regarding \(\mathbf{h}_{2}(u_{1},\dots,u_{n},t)\)  
\[