Consider the following discrete-time system
$$x_{i+1} = f(x_i, u_i)$$
with initial state $x_0 \in \mathbb R^n$. Let the cost function be
$$ C(u) := \sum_{i = 0}^{n-1} g_i(x_i, u_i) + h(x_n) $$
It is known that by defining the adjoint state as follows
$$p_n = J_h(x_n), \qquad p_i = p_{i+1}J_{f, x}(x_i, u_i) + J_{g_i, x}(x_i, u_i),$$
where $J_{f,x}$ denotes the Jacobian of $f$ with respect to $x$, we can get the gradient of $C$ with respect to $u$ from the following formula
$$ \mathrm dC(u) = \sum_{i = 0}^{n-1} (J_{g_i, u} + p_{i+1} J_{f, u})\, \mathrm du_i.$$
Are there any similar methods to get the Hessian of $C$?
References discussing this problem will also be appreciated.