Stage-wise controls for IRK methods

[PierreMartinon]

A big thanks to @vmerc !

• stagewise control is now the default for GL2 and GL3 scheme, with the stepwise control still available through options gauss_legendre_2|3_constant_control
• [DONE] manual sparsity patterns are not available yet, but should not be too difficult to adapt from the stepwise control formulation
• solution currently outputs the 'averaged' control at each step, not the individual stage controls (which could be retrieved with a bit of work in the high-level generic getter). The averaged control uses the same weights b_j as the slopes in the discretized dynamics x_i+1 = x_i + h_i sum_j b_j f_ij. Recall that the f_ij are part of the NLP unknown, and subject to the equations f_ij = f(t_ij, x_ij, u_ij). The state values at the stages are treated as internal variables and computed via x_ij = x_i + sum_k a_jk f_ik. Alternately, it may be possible to use the x_ij as NLP unknowns and compute the f_ij instead with proper buffering. To be investigated.
• another related question would be the handling of path constraints, which are currently evaluated at time steps using the first stage controls, ie g(t_i, x_i, u_i1). We could use the averaged control instead of the first one, or conversely evaluate the constraint at the stage times, reusing the x_ij from the dynamics part.

@ocots @jbcaillau @joseph-gergaud this version will be the 1.0.9 release

[PierreMartinon]

@ocots @jbcaillau Manual sparsity patterns added, see #592

[jbcaillau]

@PierreMartinon @vmerc great job, thanks to both of you 👍🏽

A big thanks to @vmerc !

• stagewise control is now the default for GL2 and GL3 scheme, with the stepwise control still available through options gauss_legendre_2|3_constant_control
• [DONE] manual sparsity patterns are not available yet, but should not be too difficult to adapt from the stepwise control formulation
• solution currently outputs the 'averaged' control at each step, not the individual stage controls (which could be retrieved with a bit of work in the high-level generic getter). The averaged control uses the same weights b_j as the slopes in the discretized dynamics x_i+1 = x_i + h_i sum_j b_j f_ij. Recall that the f_ij are part of the NLP unknown, and subject to the equations f_ij = f(t_ij, x_ij, u_ij). The state values at the stages are treated as internal variables and computed via x_ij = x_i + sum_k a_jk f_ik. Alternately, it may be possible to use the x_ij as NLP unknowns and compute the f_ij instead with proper buffering. To be investigated.

you mean use x_ij instead of f_ij as NLP unknowns? looks a linear change of variables in view of x_ij = x_i + sum_k a_jk f_ik so probably not crucial

• another related question would be the handling of path constraints, which are currently evaluated at time steps using the first stage controls, ie g(t_i, x_i, u_i1). We could use the averaged control instead of the first one, or conversely evaluate the constraint at the stage times, reusing the x_ij from the dynamics part.

right:

• the current version is not consistant with stages and using the average control is indeed what should be done
• evaluating the constraint at stage times would look even more consistant but has some price (same dimension augmentation as with the stages); we need some math here to examine the gain in doing so in terms of error, convergence... check in Hager's book?

@ocots @jbcaillau @joseph-gergaud this version will be the 1.0.9 release

[ocots]

Is it interesting to reformulate the problem to have the zi as unknowns instead of ki?

Réf. Hairer et al