@@ -90,14 +90,14 @@ nothing # hide
9090Let us solve it with a direct method (we set the number of time steps to 200):
9191
9292``` @example main
93- sol = solve(ocp; grid_size=200)
93+ direct_sol = solve(ocp; grid_size=200)
9494nothing # hide
9595```
9696
9797and plot the solution:
9898
9999``` @example main
100- plt = plot(sol ; label="Direct", size=(800, 600))
100+ plt = plot(direct_sol ; label="Direct", size=(800, 600))
101101```
102102
103103!!! note "Nota bene"
@@ -161,15 +161,24 @@ We are now ready to solve the shooting equations:
161161# in-place shooting function
162162nle!(s, ξ, λ) = shoot!(s, ξ[1:2], ξ[3], ξ[4])
163163
164- # initial guess: costate and final time
165- ξ_guess = [0.1, 0.1, 0.5, 1]
164+ # initial guess for the Newton solver from direct method
165+ t = time_grid(direct_sol) # the time grid as a vector
166+ p = costate(direct_sol) # the costate as a function of time
167+ p0_guess = p(t0) # initial costate
168+ tf_guess = variable(direct_sol) # final time
169+
170+ # find switching time t1 where p2(t) changes sign
171+ p2_values = [p(ti)[2] for ti in t]
172+ t1_guess = t[findfirst(i -> i > 1 && p2_values[i] * p2_values[i-1] < 0, 2:length(p2_values))]
173+
174+ ξ_guess = [p0_guess[1], p0_guess[2], t1_guess, tf_guess]
166175
167176# NLE problem
168177prob = NonlinearProblem(nle!, ξ_guess)
169178
170179# resolution of the shooting equations
171- sol = solve(prob; show_trace=Val(true))
172- p0, t1, tf = sol .u[1:2], sol .u[3], sol .u[4]
180+ shoot_sol = solve(prob; show_trace=Val(true))
181+ p0, t1, tf = shoot_sol .u[1:2], shoot_sol .u[3], shoot_sol .u[4]
173182
174183# print the solution
175184println("\np0 = ", p0, "\nt1 = ", t1, "\ntf = ", tf)
@@ -182,10 +191,10 @@ Finally, we reconstruct and plot the solution obtained by the indirect method:
182191φ = f_max * (t1, f_min)
183192
184193# compute the solution: state, costate, control...
185- flow_sol = φ((t0, tf), x0, p0; saveat=range(t0, tf, 200))
194+ indirect_sol = φ((t0, tf), x0, p0; saveat=range(t0, tf, 200))
186195
187196# plot the solution on the previous plot
188- plot!(plt, flow_sol ; label="Indirect", color=2, linestyle=:dash)
197+ plot!(plt, indirect_sol ; label="Indirect", color=2, linestyle=:dash)
189198```
190199
191200!!! note
0 commit comments