Open araujolma opened 5 years ago
I have solved partially this issue. I ended up not implementing the _psiy inspection as originally planned, but opted for a more manual "override" in which the code of each problem can be altered (by adding a few lines to the init section) to take advantage of the changes to LMPBVP.
I do not consider this issue to be 100% solved, though. There are two more "sub-issues" that I can think of right now:
The particular solutions are avoided when the initial value of the state for the arc is prescribed. However, for some arcs, most of the values for the states at the end are prescribed. I believe that, theoretically, the particular solution could, in these cases, be propagated backwards in time in order to take advantage of that. Doing this would require a HUGE alteration to the structure of LMPBVP, so this should only be tried if we are REALLY grinding for performance.
I am pretty sure that the rest can be changed to take advantage of low values of P_psi, for example. Looking at the coefficients obtained during rest, it seems that these could also not be calculated. Maybe doing so will destabilise Rest; it must be tried to know for sure.
As we know, the Linear Multi-Point Boundary Value Problem is solved by the method of particular solutions exposed in Miele and Wang (2003). It works by solving many (2ns + p + 1) particular solutions for the same ODE, imposing unit variations in each state (for each arc), then for each costate (for each arc), then each parameter, then no variation (homogeneous solution). After solving all these ODEs, it then solves a linear system of algebraic equations to get the coefficients of the solutions, in order to compose the final solution as a linear (actually, affine) combination of particular solutions.
After solving the LMPBVP for probCart by hand, I noticed again that it makes no sense to propagate an particular solution 'j' that will most definitely not be used, only to then solve a linear system (that is larger that it needs to be) and then obtain 'k_j = 0', as expected. For example: if the boundary condition of the first state in a given arc is to impose a given value for the state, say, 'x=1', then the grad corrections for this state (A) and arc will be necessarily zero, because the boundary conditions are already satisfied. An entire 2ns system of ODEs can be skipped! This kind of boundary condition is actually very, very common.
Theoretically, it should be easy to inspect the psiy matrix and identify the particular solutions which won't be used. Coding this is not easy, but I believe the computational cost reduction would be great! This is a potentially game changing enhancement for the code.