## 2 Materials and Methods
### Model Predictive Control  
MPC uses predictions and estimates of the current state of the controlled plant to solve an optimal control cost function over a fixed horizon. The solution is an array of present and future control moves to push the plant towards optimal performance based on the adopted objective function. The first element of the array is used as the current control input, and the procedure is iteratively repeated. The most recent developments and findings in MPC are based on the state space formulation [22; 23], but the implementation results are similar across various formulations. Furthermore, the state space formulation with a velocity setup offers several advantages, such as ensuring offset-free control [23].  
Therefore, in this work, we adopt the discrete state space model with the augmented velocity model described in Equations (10) and (11), respectively.  
\[(k+1)=A_{m}x_{m}(k)+B_{m}u(k)y(k)=C_{m}x_{m}(k) \tag{10}\]  
\[y(k+1)=Ax(k)+B\Delta u(k)y(k)=Cx(k) \tag{11}\]  
\begin{table}
\begin{tabular}{c c c c} \hline \hline
**Process** & \(K_{\rm C}\) & \(\tau_{l}\) & \(\tau_{D}\) \\ \hline Underdamped & \(\frac{\tau_{l}}{K_{p}(2\lambda+\theta)}\) & \(2\xi\tau-\frac{2\lambda^{2}-\theta^{2}}{2(2\lambda+\theta)}\) & \(\tau_{l}-2\xi\tau+\frac{\tau^{2}-\frac{\theta^{2}}{2(2\lambda+\theta)}}{\tau_ {l}}\) \\ Critically damped & \(\frac{\tau_{l}}{K_{p}(2\lambda+\theta)}\) & \(2\tau-\frac{2\lambda^{2}-\theta^{2}}{2(2\lambda+\theta)}\) & \(\tau_{l}-2\tau+\frac{\tau^{2}-\frac{\theta^{2}}{2(2\lambda+\theta)}}{\tau_ {l}}\) \\ Overdamped & \(\frac{\tau_{l}}{K_{p}(2\lambda+\theta)}\) & \((\tau_{\rm P1}+\tau_{\rm P2})-\frac{2\lambda^{2}-\theta^{2}}{2(2\lambda+\theta)}\) & \(\tau_{l}-(\tau_{\rm P1}+\tau_{\rm P2})+\frac{\tau_{\rm P1}\tau_{\rm P2}-\frac {\theta^{2}}{2(2\lambda+\theta)}}{\tau_{l}}\) \\ \hline \hline \end{tabular}
\end{table}
Table 1: PID Parameters for IMC—Mac.  
The state space matrices, \(A\in\mathbb{R}^{n\times n}\), \(B\in\mathbb{R}^{n\times m}\) and \(C\in\mathbb{R}^{p\times n}\), are defined according to Equation (12):  
\[A=\begin{bmatrix}A_{m}&0_{p}^{\intercal}\\ C_{m}A_{m}&I_{p}\end{bmatrix},\ \ \begin{bmatrix}B_{m}\\ C_{m}B_{m}\end{bmatrix},\ \ C=\begin{bmatrix}0_{n\times p}&I_{p}\end{bmatrix} \tag{12}\]  
where the variables \(0\) and \(I\) represent matrices of zeros and ones or appropriate dimensions, respectively.  
The optimization cost function adopted for the MPC allows for the penalization of the tracking error and change in control moves, as shown in Equation (13). However, we fix the penalty on the controlled variable so that only the manipulated variable is penalized, leading to a more precise tuning.  
\[J=\sum_{i=1}^{p}\left\|r(k+i)-y(k+i)\right\|^{2}+\sum_{i=1}^{M}\left\|\Delta u( k+i)\right\|_{r_{w}}^{2} \tag{13}\]  
The variable \(r_{w}\) is the input rate weighting factor; \(r\) is the setpoint; \(P\) is the prediction horizon; and \(M\) is the control horizon.  
The prediction horizon, \(P\), is set to the process settling time, including dead time, according to Equation (14) to capture process dynamics fully. The control horizon, \(M\), is fixed to a value of \(M=3\) throughout this work. This choice is supported by the fact that \(M=3\) to \(5\) is sufficient for most systems [24], and setting the value higher than five offers no significant benefits. Pre-fixing the identified parameters reduces the tuning parameters from \(4\) to \(1\), making controller tuning transparent.  
\[P=\frac{\theta+5\tau}{T_{s}} \tag{14}\]  
where the variable \(T_{s}\) is the sampling period. Therefore, the weighting \(r_{w}\) dictates the aggressiveness of control, hence its robustness. The weight \(r_{w}\) was initially varied over the range \(0.001\leq r_{w}\leq 100\). This value is adjusted for different system models, as highlighted in the simulation setup and relevant sections in the results section.