|
4 | 4 |
|
5 | 5 | using OptimalControl
|
6 | 6 | using Plots
|
7 |
| -unicodeplots() |
| 7 | +#unicodeplots() |
| 8 | +ENV["GKSwstype"]="nul" # no plot display on stdout |
8 | 9 |
|
9 | 10 | # Parameters
|
10 | 11 | const Cd = 310
|
@@ -51,8 +52,10 @@ f(x, u) = F0(x) + u[1]*F1(x)
|
51 | 52 | constraint!(ocp, :dynamics, f)
|
52 | 53 |
|
53 | 54 | # Solve
|
54 |
| -sol = solve(ocp) |
| 55 | +N = 30 |
| 56 | +sol = solve(ocp, grid_size=N) |
55 | 57 | plot(sol)
|
| 58 | +savefig("sol-direct.pdf") |
56 | 59 |
|
57 | 60 | ## Indirect solve
|
58 | 61 |
|
@@ -101,7 +104,7 @@ function shoot!(s, p0, t1, t2, t3, tf) # B+ S C B0 structure
|
101 | 104 | end
|
102 | 105 |
|
103 | 106 | # Initialisation from direct solution
|
104 |
| -t = sol.T; t = (t[1:end-1] + t[2:end]) / 2 |
| 107 | +t = sol.T; N = length(t)-1; t = (t[1:end-1] + t[2:end]) / 2 |
105 | 108 | x = [ sol.X[i, 1:3] for i ∈ 1:N+1 ]; x = (x[1:end-1] + x[2:end]) / 2
|
106 | 109 | u = [ sol.U[i, 1 ] for i ∈ 1:N+1 ]; u = (u[1:end-1] + u[2:end]) / 2
|
107 | 110 | p = [ sol.P[i, 1:3] for i ∈ 1:N ]
|
@@ -143,4 +146,5 @@ m_plot = plot(sol, idxs=(0, 3), xlabel="t", label="m")
|
143 | 146 | pr_plot = plot(sol, idxs=(0, 4), xlabel="t", label="p_r")
|
144 | 147 | pv_plot = plot(sol, idxs=(0, 5), xlabel="t", label="p_v")
|
145 | 148 | pm_plot = plot(sol, idxs=(0, 6), xlabel="t", label="p_m")
|
146 |
| -plot(r_plot, pr_plot, v_plot, pv_plot, m_plot, pv_plot, layout=(3, 2)) |
| 149 | +plot(r_plot, pr_plot, v_plot, pv_plot, m_plot, pm_plot, layout=(3, 2)) |
| 150 | +savefig("sol-indirect.pdf") |
0 commit comments