## 3 Synthesis of the Algorithm for the Solution of the Zero-Sum Differential Game
### Gradient Calculation
 \tag{33}\]  
In (33) it is necessary to calculate \(\
abla_{\boldsymbol{\pi}}\hat{\mathbf{x}}(i)\), which is obtained based on the expression (27) as follows:  
\[ \begin{split}&\
abla_{\boldsymbol{\pi}}\hat{\mathbf{x}}(i)= \mathbf{a}\otimes\mathbf{B}(\mathbf{x}(i-1))+\
abla_{\mathbf{x}}[\mathbf{a} \otimes\mathbf{B}(\mathbf{x}(i-1))\boldsymbol{\pi}]\cdot\
abla_{\boldsymbol{ \pi}}\hat{\mathbf{x}}(i-1)\\ &+\
abla_{\mathbf{x}}\hat{\mathbf{f}}(\hat{\mathbf{x}}(i-1))\cdot \
abla_{\boldsymbol{\pi}}\hat{\mathbf{x}}(i-1)+\
abla_{\mathbf{x}}[\mathbf{a} \otimes\mathbf{D}(\mathbf{x}(i-1))\boldsymbol{\rho}]\cdot\
abla_{\boldsymbol{ \pi}}\hat{\mathbf{x}}(i-1),\end{split}  \tag{34}\]where  
\[\
abla_{\hat{\mathbf{x}}}[\mathbf{a}\otimes\mathbf{B}(\mathbf{x}(i-1)) \boldsymbol{\pi}] =\mathbf{a}\otimes\Big{(}\boldsymbol{\pi}^{\mathrm{T}}\otimes \mathbf{I}_{n}\Big{)}\cdot\
abla_{\hat{\mathbf{x}}}\mathrm{vec}(\mathbf{B}( \mathbf{x}(i-1))), \tag{35}\] \[\
abla_{\hat{\mathbf{x}}}[\mathbf{a}\otimes\mathbf{D}(\mathbf{x}( i-1))\boldsymbol{\rho}] =\mathbf{a}\otimes\Big{(}\boldsymbol{\rho}^{\mathrm{T}}\otimes \mathbf{I}_{n}\Big{)}\cdot\
abla_{\hat{\mathbf{x}}}\mathrm{vec}(\mathbf{D}( \mathbf{x}(i-1))). \tag{36}\]  
The expressions (34)-(36) represent recursive matrix relations for \(i=1,2,\ldots,N-1\) with initial conditions  
\[\
abla_{\boldsymbol{\pi}}\hat{\mathbf{x}}(0)=\mathbf{0}. \tag{37}\]  
Note that the functions, \(\hat{\mathbf{f}}(\hat{\mathbf{x}})\), \(\mathbf{h}(\mathbf{x})\), \(\mathbf{P}(\mathbf{x})\), \(\mathbf{R}(\mathbf{x})\), \(\mathbf{B}(\mathbf{x})\) and \(\mathbf{D}(\mathbf{x})\), which appear in previous expressions, are known from system dynamics and predetermined approximation functions basis and their derivatives can be easily calculated to machine precision by applying a complex-step method. The first-order derivative of \(\hat{\mathbf{f}}(\hat{\mathbf{x}})\) with respect to \(\hat{\mathbf{x}}\) using complex-step approximation is accomplished by approximating its component (let us say the \(k\)-th component) with a complex variable using a Taylor's series expansion [41; 42]  
\[\hat{f}_{k}(\hat{\mathbf{x}}+i\,h\,\mathbf{e}_{j})=\hat{f}_{k}(\hat{\mathbf{x }})+i\,h\,\
abla_{\hat{\mathbf{x}}}\hat{f}_{k}(\hat{\mathbf{x}})\cdot\mathbf{ e}_{j}+\frac{(i\,h)^{2}}{2!}\,\mathbf{e}_{j}^{\mathrm{T}}\cdot\
abla_{\hat{ \mathbf{x}}}^{2}\hat{f}_{k}(\hat{\mathbf{x}})\cdot\mathbf{e}_{j}+\ldots. \tag{38}\]  
where \(\mathbf{e}\) is unit vector. Taking only the imaginary parts gives  
\[\
abla_{\hat{\mathbf{x}}}\hat{f}_{k}(\hat{\mathbf{x}})\cdot\mathbf{e}_{j}= \frac{\mathrm{Im}\Big{\{}\hat{f}_{k}(\hat{\mathbf{x}}+i\,h\,\mathbf{e}_{j}) \Big{\}}}{h}+\mathcal{O}(h^{2}). \tag{39}\]  
In the above expressions, \(i\) represents an imaginary unit \(i^{2}=-1\), while previously with \(i\) we denoted the \(i\)-th time point. For this reason, the \(\hat{\mathbf{x}}(i)\) is omitted in these expressions, but this is implied. It can be seen that subtraction does not appear in the numerator of the expression (39), i.e., the subtractive cancellation problem has been removed. This means that the step \(h\) can be extremely small, so the derivative can be calculated to machine precision. The implementation of a complex-step method can be very simply achieved by using a high-level programming language that does not require a prior definition of the variable types, but it is possible to define complex variables automatically by applying built-in functions. The complex-step method for calculating \(\
abla_{\hat{\mathbf{x}}}\hat{\mathbf{f}}(\hat{\mathbf{x}})\) is described by the pseudocode given in Algorithm 4. The calculation of derivatives of other known functions can be obtained in an analogous manner.