## 2 Background
### Arbitrary-Order Differentiation Using HYPAD
The mathematics of multidual numbers dictate that at every operation on the real part of the number, the derivatives are formed and stored in the dual imaginary parts of the number. Therefore,after every operation, the dual imaginary parts of the number contain the derivatives with respect to the inputs [71], and the real part remains unaffected. The expression to calculate the \(m^{th}\) order partial derivatives of a function of \(j\) variables \(f(\theta)\), with \(\theta=\{\theta_{1},\ldots,\theta_{j}\}\), using multidual numbers resultant from obtaining the hypercomplex Taylor series expansion of a function is given by [72]:  
\[\frac{\partial^{m}f(\theta_{1},\ldots,\theta_{j})}{\partial\theta_{1}^{n_{1}} \ldots\partial\theta_{j}^{n_{j}}}=Im_{1\ldots m}\bigg{[}f\bigg{(}\theta_{1}+ \sum_{g=1}^{b_{1}}e_{g},\ldots,\theta_{j}+\sum_{g=q}^{q+b_{j}}e_{g}\bigg{)} \bigg{]} \tag{7}\]  
where \(=\sum_{g=1}^{j-1}b_{g}\), \(\sum_{g=1}^{j}b_{g}=m\), and \(j\) is the number of input model parameters. A use aspect of the methodology is that for computing sensitivities of a specific order, lower-order sensitivities are obtained simultaneously because that information is always available by taking different dual imaginary axes. For instance, the third-order partial derivative, \(\partial^{3}f(\theta_{1},\theta_{2},\theta_{3})/\partial\theta_{1}\theta_{2} \theta_{3}\), can be obtained from Equation (7) using \(m=3\), and \(b_{1}=b_{2}=b_{3}=1\), which gives  
\[\frac{\partial^{3}f(\theta_{1},\theta_{2},\theta_{3})}{\partial\theta_{1} \partial\theta_{2}\partial\theta_{3}}=Im_{123}[f(\theta_{1}+e_{1},\theta_{2}+ e_{2},\theta_{3}+e_{3})] \tag{8}\]  
Similarly, first- and second-order partial derivatives can be found from the same function evaluation as:  
\[\frac{\partial^{2}f(\theta_{1},\theta_{2},\theta_{3})}{\partial\theta_{1} \partial\theta_{2}}=Im_{12}[f(\theta_{1}+e_{1},\theta_{2}+e_{2},\theta_{3}+ e_{3})] \tag{9a}\]  
\[\frac{\partial^{2}f(\theta_{1},\theta_{2},\theta_{3})}{\partial\theta_{1} \partial\theta_{3}}=Im_{13}[f(\theta_{1}+e_{1},\theta_{2}+e_{2},\theta_{3}+ e_{3})] \tag{9b}\]  
\[\frac{\partial^{2}f(\theta_{1},\theta_{2},\theta_{3})}{\partial\theta_{2} \partial\theta_{3}}=Im_{23}[f(\theta_{1}+e_{1},\theta_{2}+e_{2},\theta_{3}+ e_{3})] \tag{9c}\]  
\[\frac{\partial f(\theta_{1},\theta_{2},\theta_{3})}{\partial\theta_{1}}=Im_{1 }[f(\theta_{1}+e_{1},\theta_{2}+e_{2},\theta_{3}+e_{3})] \tag{9d}\]  
\[\frac{\partial f(\theta_{1},\theta_{2},\theta_{3})}{\partial\theta_{2}}=Im_{2 }[f(\theta_{1}+e_{1},\theta_{2}+e_{2},\theta_{3}+e_{3})] \tag{9e}\]  
\[\frac{\partial f(\theta_{1},\theta_{2},\theta_{3})}{\partial\theta_{3}}=Im_{3 }[f(\theta_{1}+e_{1},\theta_{2}+e_{2},\theta_{3}+e_{3})] \tag{9f}\]  
For more background about HYPAD and multidual numbers, the interested reader is directed to [48; 49; 69]. The use of HYPAD to compute sensitivities requires the use of multidual numbers. However, most numerical packages do not have hypercomplex algebra capabilities and external libraries are required. As a result, a hypercomplex library MultiZ [73] was developed that defines hypercomplex variable types, functions, and operations using operator overloading for both the Fortran and Python languages. MultiZ can handle arbitrary-order multidual and multicomplex numbers. In this work, MultiZ was used as the main tool to carry out operations involving multidual mathematics. As an alternative, multidual numbers can be represented using real only matrices using a matrix notation called the Cauchy Riemann (CR) notation. In the CR notation, the multidual operations are replaced by matrix operations [26; 73; 74; 75]; however, the computational cost is increased for every order of derivative and the programming effort is considerable higher compared to the operator overloading approach.