## 2 Materials and Methods
### Theory Basics
#### 2.2.1 Model Predictive Control  
Consider the following system:  
\[\xi_{ra}(\theta+1)=M\xi_{ra}(\theta)+N\mu_{ra}(\theta)+E\omega_{ra}(\theta),a \in(1,2,...,n), \tag{5}\]where \(\xi_{ra}(\theta)\), \(\mu_{ra}(\theta)\), and \(\omega_{ra}(\theta)\) represent the state input, control input, and the disturbance term, respectively. As it was commented in the preliminaries, there is a disturbance term \(\omega_{ra}(\theta)\) to deal with. The disturbance would be considered constant over the entire prediction horizon. The augmented model results in  
\[\widetilde{\xi}_{ra}(\theta+1)=M_{d}\widetilde{\xi}_{ra}(\theta)+N_{d}\mu_{ra}( \theta), \tag{6}\]  
where  
\[\widetilde{\xi}_{ra}(\theta+1)=\begin{bmatrix}\xi_{ra}(\theta+1)\\ \omega_{ra}(\theta+1)\end{bmatrix};M_{d}=\begin{bmatrix}M&E\\ 0&I\end{bmatrix};N_{d}=\begin{bmatrix}N\\ 0\end{bmatrix}. \tag{7}\]  
The state-space controller is implemented using the incremental formula of the state-space model, which considers the control increment \(\Delta\mu_{ra}(\theta)\) instead of the control signal \(\mu_{ra}(\theta)\), where \(\Delta\mu_{ra}(\theta)=\mu_{ra}(\theta)-\mu_{ra}(\theta-1)\). Using this model, the disturbance term can be coupled to the state input and processed together. Therefore, defining a new state \(\chi_{ra}=\begin{bmatrix}\widetilde{\xi}_{ra}(\theta)&\mu_{rwi}(\theta-1) \end{bmatrix}=\begin{bmatrix}\xi_{ra}(\theta)&\omega_{ra}(\theta)&\mu_{ra}( \theta-1)\end{bmatrix}\) that depends on the previously enhanced state, the incremental model adopts the general form given by  
\[\chi_{ra}(\theta+1)=A_{r}\chi_{ra}(\theta)+B_{r}\Delta\mu_{ra}(\theta), \tag{8}\]  
where  
\[A_{r}=\begin{bmatrix}M_{rd}&N_{rd}\\ 0&I\end{bmatrix};B_{r}=\begin{bmatrix}N_{rd}\\ I\end{bmatrix}. \tag{9}\]  
Based on the state Equation (8), we formulate the following rolling optimization problem:  
**Problem 1**.: \[\begin{split}&\min_{\mu_{ra}(0;\theta_{\tau}),...,\mu_{ra}(N_{p}; \theta_{\tau})}\left.\right|_{d}(\chi_{ra},\mu_{ra})\\ =&\sum_{l=0}^{N_{p}}\left[\parallel\chi_{ra}(l;\theta_{\tau})-R_{ a,l}(\theta_{\tau})\parallel_{Q_{d}}+\parallel\mu_{ra}(l;\theta_{\tau}) \parallel_{R_{d}}\right]+\parallel\chi_{ra}(N_{p};\theta_{\tau})\parallel_{P_{ d}}\end{split}\] (10a) _subject to_ \[\chi_{ra}(l;\theta_{\tau})=A\chi_{ra}(l-1;\theta_{\tau})+B\mu_{ra} (l-1;\theta_{\tau}),\] (10b) \[\mu_{ra}(l;\theta_{\tau})\in\mathcal{U}_{a},\] (10c) \[\chi_{ra}(l;\theta_{\tau})\in\mathcal{X}_{a},\] (10d) \[\chi_{ra}(N_{p};\theta_{\tau})\in\mathcal{X}_{a}^{f}\in\mathcal{X} _{a},\] (10e) \[R_{ra,l}=\left\{\begin{array}{ll}r_{0},&\text{layer}1,\\ \frac{\sum_{l\in\mathcal{X}_{a}}\chi_{j}}{|\mathcal{X}_{a}|},&\text{layer}2,3, \ldots,n,\end{array}\right.\] (10f)  
_where \(N_{p}\) represents the prediction horizon, \(\mu_{ra}(l;\theta_{\tau})\) represents the variables to be optimized, and \(Q_{a}>0\), \(R_{a}>0\), and \(P_{a}>0\) are weight matrices._  
Equation (10b) represents the relationship satisfied by the predicted states before and after optimization, while 10c-e represent the corresponding constraints. By solving the cost function \(j_{a}(\chi_{ra},\mu_{ra},\chi_{a}^{a},\chi_{j}^{a})\), we can obtain the optimal control sequence \(U_{a}^{*}\). The first element is then used as the control input, i.e.,  
\[\mu_{ra}(\theta)=[1,0,...,0]U_{a}^{*}. \tag{11}\]