## 3 Control Method
### ESN with Gaussian Process Regression (GPR)  
ESN, initially developed by [PERSON], is a machine learning technique based on a recursive neural network [21; 22]. This architecture is structured to create a reservoir with a randomly connected recurrent network that encodes the underlying dynamics in the internal states. This configuration enables fast learning as a simple algorithm is used. Furthermore, the training procedure, where only the readout unit is trained, is relatively inexpensive. In the proposed method, ESN was applied to predict future disturbance. The time required for the optimization calculation was ensured to be within the time margin between the current and actual execution times (Figure 2). ESN involves three components, the input units, reservoir, and output units, as shown in Figure 3. The reservoir, which is entirely connected with \(K\) input and \(L\) output units, transforms the input data into nonlinear features in a high-dimensional space. To obtain the same output as the target value, the output weight matrix was updated in the output layer only. The activations of the input units, internal units, and output units at time step \(t\) are \(\mathbf{s}(t)\)\(=\)\((s_{1}(t),\cdots,s_{\mathrm{K}}(t))^{\mathrm{T}}\), \(\mathbf{r}(t)\)\(=\)\((r_{1}(t),\cdots,r_{N}(t))^{\mathrm{T}}\), and \(\mathbf{z}(t)\)\(=\)\((z_{1}(t),\cdots,z_{L}(t))^{\mathrm{T}}\), respectively. The internal activations were updated in the reservoir with leaky integrator neurons, as follows:  
\[\mathbf{r}(t) =(1-a)\mathbf{r}(t-1)+a\;f(\mathbf{W}^{in}\mathbf{s}(t)+\mathbf{W}\mathbf{r}(t-1)) \tag{33}\]  
where \(f\) is the activation function of the reservoir, and \(a\) is the leaking decay rate. The connection weights were randomly generated in an \(N\times K\) input weight matrix \(\mathbf{W}^{in}\) and an \(N\times N\) internal weight matrix \(\mathbf{W}\). Generally, these weights are derived using a normal or uniform distribution, which represents fixed values in the network. The output of the reservoir is computed as:  
\[\mathbf{z}(t) =f^{out}(\mathbf{W}^{out}(\mathbf{s}(t),\mathbf{r}(t),\mathbf{z}(t-1))) \tag{34}\]  
where \(f^{out}\) is the activation function of the output neurons. \(\mathbf{W}^{out}\) is an \(L\times(K+N+L)\) output weight matrix, which is trained using a linear regression algorithm. Furthermore, we assumed that the time-series data for learning the ESN included system noise. Additionally, a Gaussian process regression (GPR) [23] was introduced to estimate the model using the input data adaptively. GPR is a method that constructs a regression model of nonlinear time-series data by assuming the joint distribution of each value of time-series by following the Gaussian distribution. Generally, the regression process of time-series data involves estimating the regression model \(yi=f(t_{k})\) for arbitrary input \(t_{k}\) when the outputs \((y_{1},y_{2},...,y_{M})^{\text{T}}\) for time-series inputs \((t_{1},t_{2},...,t_{M})^{\text{T}}\) are given. Assuming that the joint distribution of the M-dimensional vector \(\mathbf{y}\) consisting of time-series data \(\mathbf{y}=(y_{1},y_{2},...,y_{m})^{\text{T}}\) follows a Gaussian distribution. The Gaussian process that provides the models for joint distribution can be described by the mean \(E[\mathbf{y}]\) and covariance \(Cov[\mathbf{y}]\). By estimating the mean and covariance, a probabilistic model can be generated for the given time-series data. By assuming the observed noise for \(\mathbf{y}\), the target regression model can be calculated by using the covariance matrix \(\mathbf{C}_{M}\), whose elements are the sum of the kernel functions \(k_{f}(t_{i},t_{j})\) that compose the Gram matrix and the observed noise \(\beta^{-1}\delta_{ij}\) of \(\mathbf{y}\), where \(\delta_{i,j}\) is a delta function that becomes 1 when \(i=j\). To estimate the output \(\hat{y}_{M+1}\) for the new input \(t_{M+1}\), we would have to consider the average of the Gaussian distribution at the input \(t_{M+1}\).