## 2 Materials and Methods
In this research, the Hessian matrix \(H(J)=\frac{\partial^{2}I}{\partial\overline{\theta}^{2}}+\lambda\cdot J\) was calculated for this aim, where \(\lambda\) is the regularization coefficient, and the computational technique from [25; 26; 27; 28] is additionally applied. As is known, a positive definite matrix, \(H(J)\), at the point \(\overline{\theta}^{*}\) indicates that the minimum for criterion (6) has been found:  
\[H(J)=\left[ \begin{array}{cc}\frac{\partial^{2}I}{\partial\overline{\theta} _{1}^{2}}+\lambda\cdot J&\frac{\partial^{2}I}{\partial\overline{\theta}_{1} \partial\overline{\theta}_{2}^{2}}+\lambda\cdot J\\ \frac{\partial^{2}I}{\partial\overline{\theta}_{2}\partial\overline{\theta}_ {1}}+\lambda\cdot J&\frac{\partial^{2}I}{\partial\overline{\theta}_{2}^{2}}+ \lambda\cdot J\end{array} \right]=\left[ \begin{array}{cc}\frac{\partial^{2} J}{\partial\overline{\theta}_{i}\partial\overline{\theta}_{j}}+\lambda\cdot J \end{array} \right],\ (i,j=\overline{1,2}). \tag{15}\]  
[PERSON]'s criterion for a positive definite matrix (15) is presented as follows:  
\[\frac{\partial^{2}J}{\partial\overline{\theta}_{1}^{2}}+\lambda\cdot J>0,\ \frac{ \partial^{2}J}{\partial\overline{\theta}_{2}^{2}}+\lambda\cdot J>0,\ \det(H(J))>0. \tag{16}\]  
The definition of the elements in matrix (15) is based on the second-order sensitivity functions; that is,  
\[ \begin{bmatrix}\frac{\partial\overline{\theta}_{1}^{2}(t)}{\partial\overline{ \theta}_{1}}&\frac{\partial\overline{\theta}_{1}^{2}(t)}{\partial\overline{ \theta}_{2}}\\ \frac{\partial\overline{\theta}_{2}^{2}(t)}{\partial\overline{\theta}_{1}}& \frac{\partial\overline{\theta}_{2}^{2}(t)}{\partial\overline{\theta}_{2}} \end{bmatrix} =\big{[}\xi_{ij}(t)\big{]},\ (i,j=\overline{1,2}). \tag{17}\]  
The second-order sensitivity function characterizes a change's effect on the change rate of another sensitivity function, allowing one to take into account the system parameters' interactions under nonlinear dynamics [16; 17]. This approach provides a more in-depth analysis of the parameters' influence on the system's behaviour than first-order methods and allows one to more accurately adjust the controller's settings. This is especially important for systems with significant delays and complex nonlinearities, where traditional optimization methods may not be effective enough [16]. Empirical research has confirmed that the second-order sensitivity functions used lead to accelerated convergence of the optimization algorithm and increased system stability [18]. The analytical search for sensitivity functions (15) is a labour-intensive task [8; 12; 14]. Therefore, in this research, a more straightforward method based on the sensitivity model [16; 17; 18; 25; 26] was used. Within this research framework, similarly to [25], a second-order sensitivity function model was obtained for system (1) with controller (5):  
\[\xi_{ij}(t)=G_{p}(s)\cdot\frac{\partial^{2}u(t)}{\partial t\partial \partial t},\ (i,j=\overline{1,2}),\] \[\frac{\partial^{2}u(t)}{\partial t\partial\partial t}=\frac{ \partial^{2}G_{c}(s\overline{s})}{\partial t\partial\partial t}\cdot\varepsilon \big{(}t,\overline{\theta}\big{)}-\frac{\partial G_{c}(s\overline{\theta})}{ \partial t\partial\partial t}\cdot\xi_{j}(t)-\frac{\partial G_{c}(s\overline {\theta})}{\partial t\partial\partial}\cdot\xi_{i}(t)-G_{c}\big{(}s,\overline {\theta}\big{)}\cdot\xi_{ij}(t),\] \[\frac{\partial^{2}u(t)}{\partial t\partial\partial\theta}=\left\{ \begin{array}{l}-2\cdot\xi_{1}(t)-\frac{1}{s}\cdot\xi_{11}(t),g\big{(}t, \varepsilon(t,\overline{\theta}),\dot{\varepsilon}(t,\overline{\theta})\big{)} =1,\\ -2\cdot\xi_{1}(t)-\frac{1}{s}\cdot 0(t),g\big{(}t,\varepsilon(t,\overline{\theta}), \dot{\varepsilon}(t,\overline{\theta})\big{)}=0,\end{array} \right. \tag{18}\] \[\frac{\partial^{2}u(t)}{\partial t\partial\theta_{2}}=\left\{