The adjoint equations are index 2 DAEs ([1]}), so the Runge-Kutta scheme can be used to numerically integrate them ([2]}). The Runge-Kutta scheme for adjoint equations (REF )-() is
\(&&\lambda _{fi}(k)= \nonumber \\ && \lambda _f(k+1) - \bar{h}(k+1) \sum _{j=1}^{\bar{s}} \bar{a}_{ij} \left[ -f_x^T\left(x\left( \bar{t}_j(k) \right) , z\left( \bar{t}_j(k) \right), \bar{u}(k+1)\right) \lambda _{fj}(k) \right.\nonumber \\&&\qquad \qquad \qquad \qquad \qquad \qquad \ \ \left. -g_x^T\left(x\left( \bar{t}_j(k) \right) \right) \lambda _{gj}(k) \right] =  \\&& \lambda _f(k+1) + \bar{h}(k+1) \sum _{j=1}^{\bar{s}} \bar{a}_{ij} \left[ f_x^T\left(x\left( \bar{t}_j(k) \right) , z\left( \bar{t}_j(k) \right), \bar{u}(k+1)\right) \lambda _{fj}(k) \right.\nonumber \\&&\qquad \qquad \qquad \qquad \qquad \qquad \ \ \left. +g_x^T\left(x\left( \bar{t}_j(k) \right) \right) \lambda _{gj}(k) \right], \nonumber \\&& 0 = -f^T_z\left( x\left( \bar{t}_i(k) \right) \right) \lambda _{fi}(k)\) 
