Skip to content

Commit 94d2f33

Browse files
committed
foo
1 parent a4ad9d7 commit 94d2f33

8 files changed

+94
-63
lines changed

docs/make.jl

+2-2
Original file line numberDiff line numberDiff line change
@@ -36,8 +36,8 @@ makedocs(
3636
"application-batch.md",
3737
"application-orbit.md",
3838
"application-sail.md",
39-
"application-surface-revolution.md",
40-
"application-mri-saturation.md",
39+
"Catenoid solution" => "application-surface-revolution.md",
40+
"MRI: saturation problem" => "application-mri-saturation.md",
4141
],
4242
"FGS 2024" => "fgs-2024.md",
4343
"JuliaCon 2024"=> "juliacon2024.md",

docs/src/application-orbit.md

+26-25
Original file line numberDiff line numberDiff line change
@@ -171,7 +171,7 @@ p0 = p0 / norm(p0) # Normalization |p0|=1 for free final time
171171
jshoot(ξ) = ForwardDiff.jacobian(shoot, ξ)
172172
shoot!(s, ξ) = (s[:] = shoot(ξ); nothing)
173173
jshoot!(js, ξ) = (js[:] = jshoot(ξ); nothing)
174-
bvp_sol = fsolve(shoot!, jshoot!, ξ; show_trace=true); println(bvp_sol)
174+
#bvp_sol = fsolve(shoot!, jshoot!, ξ; show_trace=true); println(bvp_sol)
175175
nothing # hide
176176
```
177177

@@ -199,36 +199,37 @@ tf = 1.210e3; p0 =-[-2.215319700438820e+01, -4.347109477345140e+01, 9.6131888072
199199
p0 = p0 / norm(p0) # Normalization |p0|=1 for free final time
200200
ξ = [tf; p0]; # Initial guess
201201
202-
bvp_sol = fsolve(shoot!, jshoot!, ξ; show_trace=true); println(bvp_sol)
202+
#bvp_sol = fsolve(shoot!, jshoot!, ξ; show_trace=true); println(bvp_sol)
203203
nothing # hide
204204
```
205205

206206
## Plots
207207

208208
```@example main
209-
tf = bvp_sol.x[1]
210-
p0 = bvp_sol.x[2:end]
211-
ode_sol = fr((0, tf), x0, p0)
212-
t = ode_sol.t; N = size(t, 1)
213-
P = ode_sol[1, :]
214-
ex = ode_sol[2, :]
215-
ey = ode_sol[3, :]
216-
hx = ode_sol[4, :]
217-
hy = ode_sol[5, :]
218-
L = ode_sol[6, :]
219-
cL = cos.(L)
220-
sL = sin.(L)
221-
w = @. 1 + ex * cL + ey * sL
222-
Z = @. hx * sL - hy * cL
223-
C = @. 1 + hx^2 + hy^2
224-
q1 = @. P *((1 + hx^2 - hy^2) * cL + 2 * hx * hy * sL) / (C * w)
225-
q2 = @. P *((1 - hx^2 + hy^2) * sL + 2 * hx * hy * cL) / (C * w)
226-
q3 = @. 2 * P * Z / (C * w)
227-
228-
plt1 = plot3d(1; xlim = (-60, 60), ylim = (-60, 60), zlim = (-5, 5), title = "Orbit transfer", legend=false)
229-
@gif for i = 1:N
230-
push!(plt1, q1[i], q2[i], q3[i])
231-
end every N ÷ min(N, 100)
209+
# tf = bvp_sol.x[1]
210+
# p0 = bvp_sol.x[2:end]
211+
# ode_sol = fr((0, tf), x0, p0)
212+
# t = ode_sol.t; N = size(t, 1)
213+
# P = ode_sol[1, :]
214+
# ex = ode_sol[2, :]
215+
# ey = ode_sol[3, :]
216+
# hx = ode_sol[4, :]
217+
# hy = ode_sol[5, :]
218+
# L = ode_sol[6, :]
219+
# cL = cos.(L)
220+
# sL = sin.(L)
221+
# w = @. 1 + ex * cL + ey * sL
222+
# Z = @. hx * sL - hy * cL
223+
# C = @. 1 + hx^2 + hy^2
224+
# q1 = @. P *((1 + hx^2 - hy^2) * cL + 2 * hx * hy * sL) / (C * w)
225+
# q2 = @. P *((1 - hx^2 + hy^2) * sL + 2 * hx * hy * cL) / (C * w)
226+
# q3 = @. 2 * P * Z / (C * w)
227+
228+
# plt1 = plot3d(1; xlim = (-60, 60), ylim = (-60, 60), zlim = (-5, 5), title = "Orbit transfer", legend=false)
229+
# @gif for i = 1:N
230+
# push!(plt1, q1[i], q2[i], q3[i])
231+
# end every N ÷ min(N, 100)
232+
nothing
232233
```
233234

234235
## References

docs/src/application-surface-revolution.md

+4-3
Original file line numberDiff line numberDiff line change
@@ -1,8 +1,9 @@
11
# The surface of revolution of minimum area
22

3-
We consider the well-known surface of revolution of minimum area problem which dates
4-
back to Euler[^1] [^2]. This is a problem from
5-
[calculus of variations](https://en.wikipedia.org/wiki/Calculus_of_variation) but we consider
3+
We consider the well-known surface of revolution of minimum area problem which dates back to Euler[^1] [^2].
4+
We already know that the solution is a [catenoid](https://en.wikipedia.org/wiki/Catenoid) and we propose to retrieve
5+
this numerically and compute conjugate points, in relation with second order conditions of local optimality.
6+
This is a problem from [calculus of variations](https://en.wikipedia.org/wiki/Calculus_of_variation) but we consider
67
its optimal control version. We minimise the cost integral
78

89
```math

docs/src/tutorial-basic-example-f.md

+2-1
Original file line numberDiff line numberDiff line change
@@ -28,7 +28,8 @@ starting from the condition $x(0) = (-1, 0)$ and with the goal to reach the targ
2828
[Double integrator: energy minimisation](https://control-toolbox.org/docs/ctproblems/stable/problems/double_integrator_energy.html#DIE)
2929
for the analytical solution and details about this problem.
3030

31-
First, we need to import the `OptimalControl.jl` package to define and solve the optimal control problem. We also need to import the `Plots.jl` package to plot the solution.
31+
First, we need to import the `OptimalControl.jl` package to define the optimal control problem and `NLPModelsIpopt.jl` to solve it.
32+
We also need to import the `Plots.jl` package to plot the solution.
3233

3334
```@example main
3435
using OptimalControl

docs/src/tutorial-basic-example.md

+2-1
Original file line numberDiff line numberDiff line change
@@ -28,7 +28,8 @@ starting from the condition $x(0) = (-1, 0)$ and with the goal to reach the targ
2828
[Double integrator: energy minimisation](https://control-toolbox.org/docs/ctproblems/stable/problems/double_integrator_energy.html#DIE)
2929
for the analytical solution and details about this problem.
3030

31-
First, we need to import the `OptimalControl.jl` package to define and solve the optimal control problem. We also need to import the `Plots.jl` package to plot the solution.
31+
First, we need to import the `OptimalControl.jl` package to define the optimal control problem and `NLPModelsIpopt.jl` to solve it.
32+
We also need to import the `Plots.jl` package to plot the solution.
3233

3334
```@example main
3435
using OptimalControl

docs/src/tutorial-double-integrator.md

+2-1
Original file line numberDiff line numberDiff line change
@@ -19,7 +19,8 @@ a line without fricton.
1919
<img src="./assets/chariot.png" style="display: block; margin: 0 auto 20px auto;" width="300px">
2020
```
2121

22-
First, we need to import the `OptimalControl.jl` package to define and solve the optimal control problem. We also need to import the `Plots.jl` package to plot the solution.
22+
First, we need to import the `OptimalControl.jl` package to define the optimal control problem and `NLPModelsIpopt.jl` to solve it.
23+
We also need to import the `Plots.jl` package to plot the solution.
2324

2425
```@example main
2526
using OptimalControl

docs/src/tutorial-initial-guess.md

+27-15
Original file line numberDiff line numberDiff line change
@@ -4,14 +4,20 @@
44
CurrentModule = OptimalControl
55
```
66

7-
We present in this tutorial the different possibilities to provide an initial guess to solve an optimal control problem using the [`solve`](@ref) command. For the illustrations, we define the following optimal control problem.
7+
We present in this tutorial the different possibilities to provide an initial guess to solve an
8+
optimal control problem using the [`solve`](@ref) command.
9+
10+
First, we need to import the `OptimalControl.jl` package to define the optimal control problem and `NLPModelsIpopt.jl` to solve it.
11+
We also need to import the `Plots.jl` package to plot the solution.
812

913
```@example main
1014
using OptimalControl
1115
using NLPModelsIpopt
1216
using Plots
1317
```
1418

19+
For the illustrations, we define the following optimal control problem.
20+
1521
```@example main
1622
t0 = 0
1723
tf = 10
@@ -36,8 +42,6 @@ This will default to initialize all variables to 0.1.
3642
```@example main
3743
# solve the optimal control problem without initial guess
3844
sol = solve(ocp, display=false)
39-
40-
# print the number of iterations
4145
println("Number of iterations: ", sol.iterations)
4246
nothing # hide
4347
```
@@ -48,10 +52,12 @@ Let us plot the solution of the optimal control problem.
4852
plot(sol, size=(600, 450))
4953
```
5054

51-
Note that the following formulations are equivalent
55+
Note that the following formulations are equivalent to not giving an initial guess.
56+
5257
```@example main
5358
sol = solve(ocp, display=false, init=nothing)
5459
println("Number of iterations: ", sol.iterations)
60+
5561
sol = solve(ocp, display=false, init=())
5662
println("Number of iterations: ", sol.iterations)
5763
nothing # hide
@@ -65,20 +71,23 @@ Except when initializing from a solution, the arguments are to be passed as a na
6571
We first illustrate the constant initial guess, using vectors or scalars according to the dimension.
6672

6773
```@example main
68-
# solve the optimal control problem with initial guess
74+
# solve the optimal control problem with initial guess with constant values
6975
sol = solve(ocp, display=false, init=(state=[-0.2, 0.1], control=-0.2, variable=0.05))
70-
71-
# print the number of iterations
7276
println("Number of iterations: ", sol.iterations)
7377
nothing # hide
7478
```
7579

7680
Partial initializations are also valid, as shown below. Note the ending comma when a single argument is passed (tuple).
7781
```@example main
82+
# initialisation only on the state
7883
sol = solve(ocp, display=false, init=(state=[-0.2, 0.1],))
7984
println("Number of iterations: ", sol.iterations)
85+
86+
# initialisation only on the control
8087
sol = solve(ocp, display=false, init=(control=-0.2,))
8188
println("Number of iterations: ", sol.iterations)
89+
90+
# initialisation only on the state and the variable
8291
sol = solve(ocp, display=false, init=(state=[-0.2, 0.1], variable=0.05))
8392
println("Number of iterations: ", sol.iterations)
8493
nothing # hide
@@ -92,10 +101,7 @@ For the state and control, we can also provide functions of time as initial gues
92101
x(t) = [ -0.2t, 0.1t ]
93102
u(t) = -0.2t
94103
95-
# solve the optimal control problem with initial guess
96104
sol = solve(ocp, display=false, init=(state=x, control=u, variable=0.05))
97-
98-
# print the number of iterations
99105
println("Number of iterations: ", sol.iterations)
100106
nothing # hide
101107
```
@@ -107,49 +113,55 @@ For the values to be interpolated both matrices and vectors of vectors are allow
107113
Simple vectors are also allowed for variables of dimension 1.
108114

109115
```@example main
116+
# initial guess as vector of points
110117
time_vec = LinRange(t0,tf,4)
111118
x_vec = [[0, 0], [-0.1, 0.3], [-0.15,0.4], [-0.3, 0.5]]
112119
u_vec = [0, -0.8, -0.3, 0]
113120
114121
sol = solve(ocp, display=false, init=(time=time_vec, state=x_vec, control=u_vec, variable=0.05))
115-
116122
println("Number of iterations: ", sol.iterations)
117123
nothing # hide
118124
```
119125

120126
Note: in the free final time case, the given time grid should be consistent with the initial guess provided for the final time (in the optimization variables).
121127

122128
## Mixed initial guess
129+
123130
The constant, functional and vector initializations can be mixed, for instance as
131+
124132
```@example main
133+
# we can mix constant values with functions of time
125134
sol = solve(ocp, display=false, init=(state=[-0.2, 0.1], control=u, variable=0.05))
126135
println("Number of iterations: ", sol.iterations)
127-
nothing # hide
128136
137+
# wa can mix every possibility
129138
sol = solve(ocp, display=false, init=(time=time_vec, state=x_vec, control=u, variable=0.05))
130139
println("Number of iterations: ", sol.iterations)
131140
nothing # hide
132141
```
133142

134143
## Solution as initial guess (warm start)
144+
135145
Finally, we can use an existing solution to provide the initial guess.
136146
The dimensions of the state, control and optimization variable must coincide.
137147
This particular feature allows an easy implementation of discrete continuations.
148+
138149
```@example main
139150
# generate the initial solution
140151
sol_init = solve(ocp, display=false)
141152
142153
# solve the problem using solution as initial guess
143154
sol = solve(ocp, init=sol_init, display=false)
144-
145-
# print the number of iterations
146155
println("Number of iterations: ", sol.iterations)
147156
nothing # hide
148157
```
149158

150-
Note that you can also manually pick and choose which data to reuse from a solution, by recovering the functions ```sol.state```, ```sol.control``` and the values ```sol.variable```.
159+
Note that you can also manually pick and choose which data to reuse from a solution, by recovering the
160+
functions ```sol.state```, ```sol.control``` and the values ```sol.variable```.
151161
For instance the following formulation is equivalent to the ```init=sol``` one.
162+
152163
```@example main
164+
# use a previous solution to initialise picking information
153165
sol = solve(ocp, display=false, init=(state=sol.state, control=sol.control, variable=sol.variable))
154166
println("Number of iterations: ", sol.iterations)
155167
nothing # hide

docs/src/tutorial-nlp.md

+29-15
Original file line numberDiff line numberDiff line change
@@ -6,9 +6,10 @@ CurrentModule = OptimalControl
66

77
We describe here some more advanced operations related to the discretized optimal control problem.
88
When calling ```solve(ocp)``` three steps are performed internally:
9-
- first, the OCP is discretized into a DOCP (a nonlinear optimization problem) with *direct_transcription*
10-
- then, this DOCP is solved, also with the method *solve*
11-
- finally, a functional solution of the OCP is rebuilt from the solution of the discretized problem, with *build_solution*
9+
10+
- first, the OCP is discretized into a DOCP (a nonlinear optimization problem) with `direct_transcription`,
11+
- then, this DOCP is solved, also with the method `solve`,
12+
- finally, a functional solution of the OCP is rebuilt from the solution of the discretized problem, with `build_solution`.
1213

1314
These steps can also be done separately, for instance if you want to use your own NLP solver. Let us load the modules
1415

@@ -33,43 +34,56 @@ end
3334
nothing # hide
3435
```
3536

36-
First let us discretize the problem
37+
First let us discretize the problem.
38+
3739
```@example main
3840
docp = direct_transcription(ocp)
3941
nothing # hide
4042
```
41-
The DOCP contains a copy of the original OCP, and the resulting discretized problem, in our case an *ADNLPModel*.
42-
You can extract this raw NLP problem with the *get_nlp* function
43+
44+
The DOCP contains a copy of the original OCP, and the resulting discretized problem, in our case an `ADNLPModel`.
45+
You can extract this raw NLP problem with the `get_nlp` function.
46+
4347
```@example main
4448
nlp = get_nlp(docp)
4549
```
50+
4651
You could then use a custom solver that would return the solution for the NLP problem, such as
52+
4753
```
4854
nlp_sol = MySolver(get_nlp(docp))
4955
```
50-
For illustrative purpose we can mimick this by using the *solve* from *CTDirect* on the DOCP, and extract the NLP solution and multipliers
56+
57+
For an example we use the `ipopt` solver from `NLPModelsIpopt.jl` package to solve the NLP problem.
58+
5159
```@example main
52-
using CTDirect
53-
nlp_sol = CTDirect.solve(docp, display=false)
60+
using NLPModelsIpopt
61+
nlp_sol = ipopt(get_nlp(docp); print_level=4, mu_strategy="adaptive", tol=1e-8, sb="yes")
5462
nothing # hide
5563
```
56-
Then we can rebuild and plot an OCP solution (note that the multipliers are optional, but the OCP costate will not be retrieved if the multipliers are not provided)
64+
65+
Then we can rebuild and plot an OCP solution (note that the multipliers are optional, but the OCP costate will not be retrieved if the multipliers are not provided).
66+
5767
```@example main
5868
sol = build_solution(docp, primal=nlp_sol.solution, dual=nlp_sol.multipliers)
5969
plot(sol)
6070
```
6171

62-
An initial guess, including warm start, can be passed to *direct_transcription* the same way as for *solve*
72+
An initial guess, including warm start, can be passed to `direct_transcription` the same way as for `solve`.
73+
6374
```@example main
6475
docp = direct_transcription(ocp, init=sol)
6576
nothing # hide
6677
```
67-
and it can also be changed after the transcription is done, with *set_initial_guess*
78+
79+
It can also be changed after the transcription is done, with `set_initial_guess`.
80+
6881
```@example main
6982
set_initial_guess(docp, sol)
7083
nothing # hide
7184
```
7285

73-
Back to the function *solve*, passing an explicit initial guess to
74-
- *solve(ocp)* will transmit it to *direct_transcription*, therefore the resulting DOCP will have the initial guess embedded in it.
75-
- *solve(docp)* will override the initial guess in the DOCP and used the given one instead, without modifying the DOCP.
86+
Back to the function `solve`, passing an explicit initial guess to
87+
88+
- `solve(ocp)` will transmit it to `direct_transcription`, therefore the resulting DOCP will have the initial guess embedded in it.
89+
- `solve(docp)` will override the initial guess in the DOCP and used the given one instead, without modifying the DOCP.

0 commit comments

Comments
 (0)