@@ -7,9 +7,8 @@ const t0 = 0
7
7
const tf = 1
8
8
const x0 = - 1
9
9
const xf = 0
10
- const α = 1.5
10
+ const α = 1.5
11
11
ocp = @def begin
12
-
13
12
t ∈ [t0, tf], time
14
13
x ∈ R, state
15
14
u ∈ R, control
@@ -19,8 +18,7 @@ ocp = @def begin
19
18
20
19
ẋ (t) == - x (t) + α * x (t)^ 2 + u (t)
21
20
22
- ∫ ( 0.5 u (t)^ 2 ) → min
23
-
21
+ ∫ (0.5 u (t)^ 2 ) → min
24
22
end
25
23
nothing # hide
26
24
@@ -29,7 +27,7 @@ u(x, p) = p
29
27
30
28
π ((x, p)) = x
31
29
32
- S (p0) = π ( φ (t0, x0, p0, tf) ) - xf # shooting function
30
+ S (p0) = π (φ (t0, x0, p0, tf)) - xf # shooting function
33
31
34
32
ξ = [0.1 ] # initial guess
35
33
@@ -40,19 +38,20 @@ function fsolve(f, j, x; kwargs...)
40
38
catch e
41
39
println (" Erreur using MINPACK" )
42
40
println (e)
43
- println (" hybrj not supported. Replaced by hybrd even if it is not visible on the doc." )
41
+ println (
42
+ " hybrj not supported. Replaced by hybrd even if it is not visible on the doc."
43
+ )
44
44
MINPACK. fsolve (f, x; kwargs... )
45
45
end
46
46
end
47
47
48
48
using DifferentiationInterface
49
- import ForwardDiff
49
+ using ForwardDiff : ForwardDiff
50
50
backend = AutoForwardDiff ()
51
51
52
- nle! = ( s, ξ) -> s[1 ] = S (ξ[1 ]) # auxiliary function
52
+ nle! = (s, ξ) -> s[1 ] = S (ξ[1 ]) # auxiliary function
53
53
jnle! = (js, ξ) -> jacobian! (nle!, similar (ξ), js, backend, ξ) #
54
54
55
-
56
55
indirect_sol = fsolve (nle!, jnle!, ξ; show_trace= true ) # resolution of S(p0) = 0
57
56
p0_sol = indirect_sol. x[1 ] # costate solution
58
57
println (" \n costate: p0 = " , p0_sol)
@@ -61,71 +60,85 @@ println("shoot: |S(p0)| = ", abs(S(p0_sol)), "\n")
61
60
sol = φ ((t0, tf), x0, p0_sol)
62
61
plot (sol)
63
62
64
- using Plots. PlotMeasures
65
- function Plots. plot (S:: Function , p0:: Float64 ; Np0= 20 , kwargs... )
66
-
63
+ using Plots. PlotMeasures
64
+ function Plots. plot (S:: Function , p0:: Float64 ; Np0= 20 , kwargs... )
65
+
67
66
# times for wavefronts
68
- times = range (t0, tf, length= 3 )
67
+ times = range (t0, tf; length= 3 )
69
68
70
69
# times for trajectories
71
- tspan = range (t0, tf, length= 100 )
70
+ tspan = range (t0, tf; length= 100 )
72
71
73
72
# interval of initial covector
74
- p0_min = - 0.5
75
- p0_max = 2
73
+ p0_min = - 0.5
74
+ p0_max = 2
76
75
77
76
# covector solution
78
- p0_sol = p0
79
-
77
+ p0_sol = p0
78
+
80
79
# plot of the flow in phase space
81
- plt_flow = plot ()
82
- p0s = range (p0_min, p0_max, length= Np0)
83
- for i ∈ eachindex (p0s)
80
+ plt_flow = plot ()
81
+ p0s = range (p0_min, p0_max; length= Np0)
82
+ for i in eachindex (p0s)
84
83
sol = φ ((t0, tf), x0, p0s[i])
85
84
x = state (sol).(tspan)
86
85
p = costate (sol).(tspan)
87
- label = i== 1 ? " extremals" : false
88
- plot! (plt_flow, x, p, color= :blue , label= label)
89
- end
90
-
86
+ label = i== 1 ? " extremals" : false
87
+ plot! (plt_flow, x, p; color= :blue , label= label)
88
+ end
89
+
91
90
# plot of wavefronts in phase space
92
- p0s = range (p0_min, p0_max, length= 200 )
93
- xs = zeros (length (p0s), length (times))
94
- ps = zeros (length (p0s), length (times))
95
- for i ∈ eachindex (p0s)
96
- sol = φ ((t0, tf), x0, p0s[i], saveat= times)
97
- xs[i, :] .= state (sol).(times)
98
- ps[i, :] .= costate (sol).(times)
99
- end
100
- for j ∈ eachindex (times)
101
- label = j== 1 ? " flow at times" : false
102
- plot! (plt_flow, xs[:, j], ps[:, j], color= :green , linewidth= 2 , label= label)
103
- end
104
-
91
+ p0s = range (p0_min, p0_max; length= 200 )
92
+ xs = zeros (length (p0s), length (times))
93
+ ps = zeros (length (p0s), length (times))
94
+ for i in eachindex (p0s)
95
+ sol = φ ((t0, tf), x0, p0s[i]; saveat= times)
96
+ xs[i, :] .= state (sol).(times)
97
+ ps[i, :] .= costate (sol).(times)
98
+ end
99
+ for j in eachindex (times)
100
+ label = j== 1 ? " flow at times" : false
101
+ plot! (plt_flow, xs[:, j], ps[:, j]; color= :green , linewidth= 2 , label= label)
102
+ end
103
+
105
104
#
106
- plot! (plt_flow, xlims= (- 1.1 , 1 ), ylims= (p0_min, p0_max))
107
- plot! (plt_flow, [0 , 0 ], [p0_min, p0_max], color= :black , xlabel= " x" , ylabel= " p" , label= " x=xf" )
108
-
105
+ plot! (plt_flow; xlims= (- 1.1 , 1 ), ylims= (p0_min, p0_max))
106
+ plot! (
107
+ plt_flow,
108
+ [0 , 0 ],
109
+ [p0_min, p0_max];
110
+ color= :black ,
111
+ xlabel= " x" ,
112
+ ylabel= " p" ,
113
+ label= " x=xf" ,
114
+ )
115
+
109
116
# solution
110
117
sol = φ ((t0, tf), x0, p0_sol)
111
118
x = state (sol).(tspan)
112
119
p = costate (sol).(tspan)
113
- plot! (plt_flow, x, p, color= :red , linewidth= 2 , label= " extremal solution" )
114
- plot! (plt_flow, [x[end ]], [p[end ]], seriestype= :scatter , color= :green , label= false )
115
-
120
+ plot! (plt_flow, x, p; color= :red , linewidth= 2 , label= " extremal solution" )
121
+ plot! (plt_flow, [x[end ]], [p[end ]]; seriestype= :scatter , color= :green , label= false )
122
+
116
123
# plot of the shooting function
117
- p0s = range (p0_min, p0_max, length= 200 )
118
- plt_shoot = plot (xlims= (p0_min, p0_max), ylims= (- 2 , 4 ), xlabel= " p₀" , ylabel= " y" )
119
- plot! (plt_shoot, p0s, S, linewidth= 2 , label= " S(p₀)" , color= :green )
120
- plot! (plt_shoot, [p0_min, p0_max], [0 , 0 ], color= :black , label= " y=0" )
121
- plot! (plt_shoot, [p0_sol, p0_sol], [- 2 , 0 ], color= :black , label= " p₀ solution" , linestyle= :dash )
122
- plot! (plt_shoot, [p0_sol], [0 ], seriestype= :scatter , color= :green , label= false )
123
-
124
+ p0s = range (p0_min, p0_max; length= 200 )
125
+ plt_shoot = plot (; xlims= (p0_min, p0_max), ylims= (- 2 , 4 ), xlabel= " p₀" , ylabel= " y" )
126
+ plot! (plt_shoot, p0s, S; linewidth= 2 , label= " S(p₀)" , color= :green )
127
+ plot! (plt_shoot, [p0_min, p0_max], [0 , 0 ]; color= :black , label= " y=0" )
128
+ plot! (
129
+ plt_shoot,
130
+ [p0_sol, p0_sol],
131
+ [- 2 , 0 ];
132
+ color= :black ,
133
+ label= " p₀ solution" ,
134
+ linestyle= :dash ,
135
+ )
136
+ plot! (plt_shoot, [p0_sol], [0 ]; seriestype= :scatter , color= :green , label= false )
137
+
124
138
# final plot
125
- plot (plt_flow, plt_shoot; layout= (1 ,2 ), leftmargin= 15 px, bottommargin= 15 px, kwargs... )
126
-
139
+ plot (plt_flow, plt_shoot; layout= (1 , 2 ), leftmargin= 15 px, bottommargin= 15 px, kwargs... )
127
140
end
128
141
129
142
p0_sol
130
143
131
- plot (S, p0_sol; size= (800 , 450 ))
144
+ plot (S, p0_sol; size= (800 , 450 ))
0 commit comments