I am coding the resolution of my optimal control problem on Julia.
u and v are my two controls and x=[s,i,m]
I implemented this problem on Julia, and want to solve it with IPOPT. I approached this integral by using the Trapezoidal rule, like so : @expression(sys, J1, 0.5 * Δt * sum((A * x[t] + B * u[t] + C * v[t] + A * i[t+1] + B * u[t+1] + C * v[t+1]) for t in 1:N)) . N being the size of my discretisation.
I defined my constraints and dynamics, then used the Cranck-Nickolson scheme.
My problem is that the solution behaves weirdly; for N \in [200, 7000] my objective is ascendant, i obtain : J1=26 for N=200 and J1=300 for N=7000. Starting from around N=7200, my objective is close to 13, by N=10000, J1=11.
It seems really weird to me to have this change of behaviour, with u that had a singular arc for N \in [200, 7000] and that is now (N=10000) bang bang.
I suspect i have a problem in my code that gives me these weirdlooking results but i dont know where.
I tried using the Runge-kutta method by replacing u and v by their optimal values (found by the solver ) in my dynamics and trying to solve the system to see if s, i and m given by Runge-Kutta coincide with the values initially given by the solver, but i can say that my results do not show that x_solver = x_runge-kutta.
One of the possibilities would be that my approximation of the integral is incorrect,but i'm not sure.
I thank you all for your help !