## 3 Control Method
### Solving the Predictive Control Problem  
The predicted state for next time step \(\hat{\mathbf{x}}_{k+1}\) was defined based on the current estimated state \(\hat{\mathbf{x}}_{k}\) filtered by Kalman filter and using the given predicted values \(\hat{\mathbf{u}}_{k}\) and \(\hat{\mathbf{w}}_{k}\) as follows:  
\[\hat{\mathbf{x}}_{k+1}=\mathbf{A}\hat{\mathbf{x}}_{k}+\mathbf{B}\hat{\mathbf{u}}_{k}+\mathbf{B}_{e}\hat {\mathbf{w}}_{k},\ \ \hat{\mathbf{y}}_{k}=\mathbf{C}\hat{\mathbf{x}}_{k}. \tag{18}\]  
The predicted control input \(\hat{\mathbf{u}}_{k}\) can be expressed as the sum of the predicted input difference \(\Delta\hat{\mathbf{u}}_{k}\) and control input at the previous step \(\mathbf{u}_{k-1}\) as follows:  
\[\hat{\mathbf{u}}_{k} = \Delta\hat{\mathbf{u}}_{k}+\mathbf{u}_{k-1}. \tag{19}\]  
\(\hat{\mathbf{w}}_{k}\) can be estimated using the ESN in the proposed ESN-MPC. Conversely, a conventional EDT-MPC follows the method that is described in Section 3.5. Assume that the variable \(\Delta\mathbf{U}_{k}=[\Delta\hat{\mathbf{u}}_{k},\cdots,\Delta\hat{\mathbf{u}}_{k+H_{u}-1} ]^{\intercal}\) represents the optimal control sequence at the time step \(k\). \(\Delta\hat{\mathbf{u}}_{k}\) in Equation (19) is derived by Equation (22) as the control sequence \(\Delta\mathbf{U}_{k}\), which is described later. The predicted output \(\mathbf{Y}_{k}\) can be expressed as follows [19]:  
\[\mathbf{Y}_{k}=\mathbf{\
abla}\hat{\mathbf{x}}_{k}+\mathbf{\gamma}\mathbf{u}_{k-1}+\mathbf{\Theta} \Delta\mathbf{U}_{k}+\mathbf{\Xi}\mathbf{W}_{mk} \tag{20}\]  
where  
\[\mathbf{Y}_{k}= \begin{bmatrix}\hat{\mathbf{y}}_{k+1}\\ \vdots\\ \hat{\mathbf{y}}_{k+H_{u}}\\ \hat{\mathbf{y}}_{k+H_{u}+1}\\ \vdots\\ \hat{\mathbf{y}}_{k+H_{p}}\end{bmatrix} ,\ \ \ \mathbf{W}_{mk}= \begin{bmatrix}\hat{\mathbf{w}}_{k}\\ \vdots\\ \hat{\mathbf{w}}_{k+H_{u}-1}\\ \hat{\mathbf{w}}_{k+H_{u}}\\ \vdots\\ \hat{\mathbf{w}}_{k+H_{p}-1}\end{bmatrix} ,\ \ \ \mathbf{\Psi}= \begin{bmatrix}\mathbf{C}\mathbf{A}\\ \vdots\\ \mathbf{C}\mathbf{A}^{H_{u}}\\ \mathbf{C}\mathbf{A}^{H_{u}+1}\\ \vdots\\ \mathbf{C}\mathbf{A}^{H_{p}}\end{bmatrix} ,\ \ \ \mathbf{\gamma}= \begin{bmatrix}\mathbf{C}\mathbf{B}\\ \vdots\\ \sum_{i=0}^{H_{u}-1}\mathbf{C}\mathbf{A}^{i}\mathbf{B}\\ \sum_{i=0}^{H_{u}}\mathbf{C}\mathbf{A}^{i}\mathbf{B}\\ \vdots\\ \sum_{i=0}^{H_{p}-1}\mathbf{C}\mathbf{A}^{i}\mathbf{B}\end{bmatrix} ,\]\[\mathbf{\Theta}=\left[ \begin{array}{cccc}\mathbf{CB}&0&\cdots&0\\ \mathbf{CAB}+\mathbf{CB}&\mathbf{CB}&\cdots&0\\ \vdots&\cdots&\ddots&\vdots\\ \sum_{i=0}^{H_{u}-1}\mathbf{CA}^{i}\mathbf{B}&\cdots&\cdots&\mathbf{CB}\\ \sum_{i=0}^{H_{u}}\mathbf{CA}^{i}\mathbf{B}&\cdots&\cdots&\mathbf{CAB}+\mathbf{CB}\\ \vdots&\cdots&\ddots&\vdots\\ \sum_{i=0}^{H_{p}-1}\mathbf{CA}^{i}\mathbf{B}&\sum_{i=0}^{H_{p}-2}\mathbf{CA}^{i}\mathbf{B}& \cdots&\sum_{i=0}^{H_{p}-H_{u}}\mathbf{CA}^{i}\mathbf{B}\end{array} \right],\quad\mathbf{ \Xi}=\left[ \begin{array}{cccc}\mathbf{CB}_{e}&\cdots 0&\cdots&0\\ \mathbf{CAB}_{e}&\mathbf{CB}_{e}&\cdots&0\\ \vdots&\vdots&\ddots&\vdots\\ \mathbf{CA}^{H_{p}-1}\mathbf{B}_{e}&\mathbf{CA}^{H_{p}-2}\mathbf{B}_{e}&\cdots&\mathbf{CB}_{e} \end{array} \right]\]  
where \(H_{u}\) and \(H_{p}\) represent the control horizon and prediction horizon, respectively, and \(H_{u}\leq H_{p}\).