Skip to content

Commit 419639c

Browse files
ocotsocots
authored
Merge pull request #104 from control-toolbox/develop
Develop
2 parents fb2ed6b + a25fc1d commit 419639c

Some content is hidden

Large Commits have some content hidden by default. Use the searchbox below for content that may be hidden.

57 files changed

+1669
-2852
lines changed

‎examples/basic_f.ipynb

+2-226
Large diffs are not rendered by default.

‎examples/callbacks-print.jl

+19-15
Original file line numberDiff line numberDiff line change
@@ -4,18 +4,22 @@ using Printf
44
using LinearAlgebra
55

66
# description of the optimal control problem
7-
t0 = 0.0 # t0 is fixed
8-
tf = 1.0 # tf is fixed
9-
x0 = [-1.0; 0.0] # the initial condition is fixed
10-
xf = [0.0; 0.0] # the target is fixed
11-
A = [0.0 1.0
12-
0.0 0.0]
13-
B = [0.0; 1.0]
14-
f(x, u) = A * x + B * u[1]; # dynamics
15-
L(x, u) = 0.5 * u[1]^2 # integrand of the Lagrange cost
16-
17-
# problem definition
18-
ocp = OptimalControlProblem(L, f, t0, x0, tf, xf, 2, 1)
7+
t0 = 0
8+
tf = 1
9+
x0 = [ -1, 0 ]
10+
xf = [ 0, 0 ]
11+
ocp = Model()
12+
state!(ocp, 2) # dimension of the state
13+
control!(ocp, 1) # dimension of the control
14+
time!(ocp, [t0, tf])
15+
constraint!(ocp, :initial, x0)
16+
constraint!(ocp, :final, xf)
17+
A = [ 0 1
18+
0 0 ]
19+
B = [ 0
20+
1 ]
21+
constraint!(ocp, :dynamics, (x, u) -> A*x + B*u[1])
22+
objective!(ocp, :lagrange, (x, u) -> 0.5u[1]^2)
1923

2024
# replace default callback
2125
function myprint(i, sᵢ, dᵢ, Uᵢ, gᵢ, fᵢ)
@@ -27,15 +31,15 @@ function myprint(i, sᵢ, dᵢ, Uᵢ, gᵢ, fᵢ)
2731
@printf("%16.8e", norm(Uᵢ) > 1e-14 ? norm(sᵢ * dᵢ) / norm(Uᵢ) : norm(sᵢ * dᵢ)) # Stagnation
2832
end
2933
cbs = (PrintCallback(myprint),)
30-
sol = solve(ocp, callbacks=cbs);
34+
sol = solve(ocp, :direct, :shooting, callbacks=cbs);
3135

3236
# add text to default print callback
3337
function myprint2(i, sᵢ, dᵢ, Uᵢ, gᵢ, fᵢ)
3438
symbols = ("--", "\\", "|", "/")
3539
print(" ", symbols[mod(i, 4)+1])
3640
end
3741
cbs = (PrintCallback(myprint2, priority=0),)
38-
sol = solve(ocp, callbacks=cbs);
42+
sol = solve(ocp, :direct, :shooting, callbacks=cbs);
3943

4044
# add text to default print callback, saving old value of f
4145
global old_f = 0.0
@@ -48,4 +52,4 @@ function myprint3(i, sᵢ, dᵢ, xᵢ, gᵢ, fᵢ)
4852
old_f = fᵢ
4953
end
5054
cbs = (PrintCallback(myprint3, priority=0),)
51-
sol = solve(ocp, callbacks=cbs);
55+
sol = solve(ocp, :direct, :shooting, callbacks=cbs);

‎examples/callbacks-stop.jl

+19-15
Original file line numberDiff line numberDiff line change
@@ -4,18 +4,22 @@ using Printf
44
using LinearAlgebra
55

66
# description of the optimal control problem
7-
t0 = 0.0 # t0 is fixed
8-
tf = 1.0 # tf is fixed
9-
x0 = [-1.0; 0.0] # the initial condition is fixed
10-
xf = [0.0; 0.0] # the target is fixed
11-
A = [0.0 1.0
12-
0.0 0.0]
13-
B = [0.0; 1.0]
14-
f(x, u) = A * x + B * u[1]; # dynamics
15-
L(x, u) = 0.5 * u[1]^2 # integrand of the Lagrange cost
16-
17-
# problem definition
18-
ocp = OptimalControlProblem(L, f, t0, x0, tf, xf, 2, 1)
7+
t0 = 0
8+
tf = 1
9+
x0 = [ -1, 0 ]
10+
xf = [ 0, 0 ]
11+
ocp = Model()
12+
state!(ocp, 2) # dimension of the state
13+
control!(ocp, 1) # dimension of the control
14+
time!(ocp, [t0, tf])
15+
constraint!(ocp, :initial, x0)
16+
constraint!(ocp, :final, xf)
17+
A = [ 0 1
18+
0 0 ]
19+
B = [ 0
20+
1 ]
21+
constraint!(ocp, :dynamics, (x, u) -> A*x + B*u[1])
22+
objective!(ocp, :lagrange, (x, u) -> 0.5u[1]^2)
1923

2024
# replace default stop callbacks
2125
function mystop(i, sᵢ, dᵢ, xᵢ, gᵢ, fᵢ, ng₀, oTol, aTol, sTol, iterations)
@@ -29,11 +33,11 @@ function mystop(i, sᵢ, dᵢ, xᵢ, gᵢ, fᵢ, ng₀, oTol, aTol, sTol, iterat
2933
return stop, stopping, message, success
3034
end
3135
cbs = (StopCallback(mystop),)
32-
sol = solve(ocp, callbacks=cbs);
36+
sol = solve(ocp, :direct, :shooting, callbacks=cbs);
3337

3438
# add stop callback to existing ones
3539
cbs = (StopCallback(mystop, priority=0),)
36-
sol = solve(ocp, callbacks=cbs);
40+
sol = solve(ocp, :direct, :shooting, callbacks=cbs);
3741

3842
# print (not replacing) and stop (replacing) callbacks
3943
global old_f = 0.0
@@ -46,4 +50,4 @@ function myprint(i, sᵢ, dᵢ, xᵢ, gᵢ, fᵢ)
4650
old_f = fᵢ
4751
end
4852
cbs = (PrintCallback(myprint, priority=0), StopCallback(mystop, priority=1))
49-
sol = solve(ocp, callbacks=cbs);
53+
sol = solve(ocp, :direct, :shooting, callbacks=cbs);

‎examples/direct-shooting.ipynb

+16-15
Original file line numberDiff line numberDiff line change
@@ -36,21 +36,22 @@
3636
],
3737
"source": [
3838
"# description of the optimal control problem\n",
39-
"t0 = 0.0 # t0 is fixed\n",
40-
"tf = 1.0 # tf is fixed\n",
41-
"x0 = [-1.0; 0.0] # the initial condition is fixed\n",
42-
"xf = [0.0; 0.0] # the target is fixed\n",
43-
"A = [0.0 1.0\n",
44-
" 0.0 0.0]\n",
45-
"B = [0.0; 1.0]\n",
46-
"f(x, u) = A * x + B * u[1]; # dynamics\n",
47-
"L(x, u) = 0.5 * u[1]^2 # integrand of the Lagrange cost\n",
48-
"\n",
49-
"# problem definition\n",
50-
"ocp = OptimalControlProblem(L, f, t0, x0, tf, xf, 2, 1, :autonomous)\n",
51-
"\n",
52-
"# display the formulation of the problem\n",
53-
"display(ocp)\n",
39+
"t0 = 0\n",
40+
"tf = 1\n",
41+
"x0 = [ -1, 0 ]\n",
42+
"xf = [ 0, 0 ]\n",
43+
"ocp = Model()\n",
44+
"state!(ocp, 2) # dimension of the state\n",
45+
"control!(ocp, 1) # dimension of the control\n",
46+
"time!(ocp, [t0, tf])\n",
47+
"constraint!(ocp, :initial, x0)\n",
48+
"constraint!(ocp, :final, xf)\n",
49+
"A = [ 0 1\n",
50+
" 0 0 ]\n",
51+
"B = [ 0\n",
52+
" 1 ]\n",
53+
"constraint!(ocp, :dynamics, (x, u) -> A*x + B*u[1])\n",
54+
"objective!(ocp, :lagrange, (x, u) -> 0.5u[1]^2)\n",
5455
"\n",
5556
"# solution of the optimal control problem\n",
5657
"N = 1001\n",

‎examples/direct-shooting.jl

+16-14
Original file line numberDiff line numberDiff line change
@@ -1,21 +1,23 @@
11
using OptimalControl
22

3-
# description of the optimal control problem
4-
t0 = 0.0 # t0 is fixed
5-
tf = 1.0 # tf is fixed
6-
x0 = [-1.0; 0.0] # the initial condition is fixed
7-
xf = [0.0; 0.0] # the target is fixed
8-
A = [0.0 1.0
9-
0.0 0.0]
10-
B = [0.0; 1.0]
11-
f(x, u) = A * x + B * u[1]; # dynamics
12-
L(x, u) = 0.5 * u[1]^2 # integrand of the Lagrange cost
3+
t0 = 0
4+
tf = 1
5+
x0 = [ -1, 0 ]
6+
xf = [ 0, 0 ]
137

14-
# problem definition
15-
ocp = OptimalControlProblem(L, f, t0, x0, tf, xf, 2, 1, :autonomous)
8+
ocp = Model()
169

17-
# display the formulation of the problem
18-
display(ocp)
10+
state!(ocp, 2) # dimension of the state
11+
control!(ocp, 1) # dimension of the control
12+
time!(ocp, [t0, tf])
13+
constraint!(ocp, :initial, x0)
14+
constraint!(ocp, :final, xf)
15+
A = [ 0 1
16+
0 0 ]
17+
B = [ 0
18+
1 ]
19+
constraint!(ocp, :dynamics, (x, u) -> A*x + B*u[1])
20+
objective!(ocp, :lagrange, (x, u) -> 0.5u[1]^2)
1921

2022
# initial iterate
2123
u_sol(t) = 6.0-12.0*t # solution

‎examples/goddard.ipynb

+20-111
Original file line numberDiff line numberDiff line change
@@ -30,20 +30,26 @@
3030
"m0 = 1\n",
3131
"mf = 0.6\n",
3232
"\n",
33-
"ocp = Model()\n",
34-
"\n",
35-
"@variable ocp tf\n",
36-
"@time ocp t ∈ [ t0, tf ]\n",
37-
"@state ocp x[3]\n",
38-
"@control ocp u\n",
39-
"\n",
40-
"@constraint ocp x(t0) == [ r0, v0, m0 ]\n",
41-
"@constraint ocp 0 ≤ u(t) ≤ 1\n",
42-
"@constraint ocp r0 ≤ x(t)[1] :state_con1\n",
43-
"@constraint ocp 0 ≤ x(t)[2] ≤ vmax :state_con2\n",
44-
"@constraint ocp mf ≤ x(t)[3] ≤ m0 :state_con3\n",
33+
"ocp = @def begin\n",
34+
"\n",
35+
" variable: tf\n",
36+
" time: t ∈ [ t0, tf ]\n",
37+
" state: x[3]\n",
38+
" control: u\n",
39+
" \n",
40+
" r = x₁\n",
41+
" v = x₂\n",
42+
" m = x₃\n",
43+
" \n",
44+
" x(t0) == [ r0, v0, m0 ]\n",
45+
" 0 ≤ u(t) ≤ 1\n",
46+
" 0 ≤ r(t) => eq1\n",
47+
" 0 ≤ v(t) ≤ vmax => eq2\n",
48+
" mf ≤ m(t) ≤ m0 => eq3\n",
4549
" \n",
46-
"@objective ocp max x(tf)[1]\n",
50+
" r(tf) -> max\n",
51+
" \n",
52+
"end\n",
4753
"\n",
4854
"function F0(x)\n",
4955
" r, v, m = x\n",
@@ -60,108 +66,11 @@
6066
"\n",
6167
"f(x, u) = F0(x) + u*F1(x)\n",
6268
"\n",
63-
"@constraint ocp ẋ(t) == f(x(t), u(t))\n",
69+
"@def! ocp ẋ(t) == f(x(t), u(t))\n",
6470
"\n",
6571
"sol = solve(ocp)\n",
6672
"plot(sol)"
6773
]
68-
},
69-
{
70-
"cell_type": "markdown",
71-
"metadata": {
72-
"tags": []
73-
},
74-
"source": [
75-
"## Indirect solve"
76-
]
77-
},
78-
{
79-
"cell_type": "code",
80-
"execution_count": null,
81-
"metadata": {
82-
"tags": []
83-
},
84-
"outputs": [],
85-
"source": [
86-
"using MINPACK\n",
87-
"\n",
88-
"# Bang controls\n",
89-
"u0(x, p) = 0.\n",
90-
"u1(x, p) = 1.\n",
91-
"\n",
92-
"# Computation of singular control of order 1\n",
93-
"H0(x, p) = p' * F0(x)\n",
94-
"H1(x, p) = p' * F1(x)\n",
95-
"H01 = Poisson(H0, H1)\n",
96-
"H001 = Poisson(H0, H01)\n",
97-
"H101 = Poisson(H1, H01)\n",
98-
"us(x, p) = -H001(x, p) / H101(x, p)\n",
99-
"\n",
100-
"# Computation of boundary control\n",
101-
"@remove_constraint ocp :state_con1\n",
102-
"@remove_constraint ocp :state_con3\n",
103-
"@constraint ocp x(tf)[3] == mf :final_con # changed to equality constraint for shooting\n",
104-
"g(x) = constraint(ocp, :state_con2, :upper)(x, 0) # g(x, u) ≥ 0 (cf. nonnegative multiplier)\n",
105-
"ub(x, _) = -Ad(F0, g)(x) / Ad(F1, g)(x)\n",
106-
"μb(x, p) = H01(x, p) / Ad(F1, g)(x)\n",
107-
"\n",
108-
"# Flows\n",
109-
"f0 = Flow(ocp, u0)\n",
110-
"f1 = Flow(ocp, u1)\n",
111-
"fs = Flow(ocp, us)\n",
112-
"fb = Flow(ocp, ub, (x, _) -> g(x), μb)\n",
113-
"\n",
114-
"# Shooting function\n",
115-
"function shoot!(s, p0, t1, t2, t3, tf) # B+ S C B0 structure\n",
116-
"\n",
117-
" x1, p1 = f1(t0, x0, p0, t1)\n",
118-
" x2, p2 = fs(t1, x1, p1, t2)\n",
119-
" x3, p3 = fb(t2, x2, p2, t3)\n",
120-
" xf, pf = f0(t3, x3, p3, tf)\n",
121-
" s[1] = constraint(ocp, :final_con)(t0, x0, tf, xf)\n",
122-
" s[2:3] = pf[1:2] - [ 1, 0 ]\n",
123-
" s[4] = H1(x1, p1)\n",
124-
" s[5] = H01(x1, p1)\n",
125-
" s[6] = g(x2)\n",
126-
" s[7] = H0(xf, pf) # free tf\n",
127-
"\n",
128-
"end\n",
129-
"\n",
130-
"# Initialisation from direct solution\n",
131-
"t = sol.time\n",
132-
"x = sol.state\n",
133-
"p = sol.adjoint\n",
134-
"H1_plot = plot(t, H1.(x, p), xlabel = \"t\", ylabel = \"H1\", legend = false)\n",
135-
"g_plot = plot(t, g.(x), xlabel = \"t\", ylabel = \"g\", legend = false)\n",
136-
"display(plot(u_plot, H1_plot, g_plot, layout = (3,1)))\n",
137-
"η = 1e-3\n",
138-
"t13 = t[ abs.(H1.(x, p)) .≤ η ]\n",
139-
"t23 = t[ 0 .≤ g.(x) .≤ η ]\n",
140-
"p0 = p[1]\n",
141-
"t1 = min(t13...)\n",
142-
"t2 = min(t23...)\n",
143-
"t3 = max(t23...)\n",
144-
"tf = t[end]\n",
145-
"ξ = [ p0 ; t1 ; t2 ; t3 ; tf ]\n",
146-
"\n",
147-
"println(\"Initial guess:\\n\", ξ)\n",
148-
"\n",
149-
"# Solve\n",
150-
"nle = (s, ξ) -> shoot!(s, ξ[1:3], ξ[4], ξ[5], ξ[6], ξ[7])\n",
151-
"sol = fsolve(nle, ξ, show_trace=true)\n",
152-
"println(sol)\n",
153-
"\n",
154-
"# Plots\n",
155-
"p0 = sol.x[1:3]\n",
156-
"t1 = sol.x[4]\n",
157-
"t2 = sol.x[5]\n",
158-
"t3 = sol.x[6]\n",
159-
"tf = sol.x[7]\n",
160-
"\n",
161-
"f = f1 * (t1, fs) * (t2, fb) * (t3, f0) # concatenation of the four Hamiltonian flows\n",
162-
"sol = f((t0, tf), x0, p0)\n",
163-
"plot(sol)"
164-
]
16574
}
16675
],
16776
"metadata": {

0 commit comments

Comments
 (0)