Skip to content

Commit e37a793

Browse files
committed
free initial time ok
1 parent b686a1f commit e37a793

File tree

1 file changed

+117
-132
lines changed

1 file changed

+117
-132
lines changed

docs/src/tutorial-free-times-initial.md

Lines changed: 117 additions & 132 deletions
Original file line numberDiff line numberDiff line change
@@ -4,228 +4,213 @@ Draft = false
44

55
# [Optimal control problem with free initial time](@id tutorial-free-times-initial)
66

7-
In this tutorial, we explore an optimal control problem with free initial time `t0`.
7+
In this tutorial, we explore a minimum-time optimal control problem where the **initial time $t_0$ is free** (but negative). We aim to determine the latest possible starting time for the system to still reach the goal at a fixed final time $t_f = 0$.
88

9-
## Definition of the problem
9+
## Problem definition
1010

11-
We consider the double integrator in minimum time. We start from the initial condition $x(t_0) = [0, 0]$ and we aim to reach at the final time $t_f = 0$ the target $x(t_f) = [1, 0]$. The final time is fixed but the initial time is free. The objective is thus to maximise `t0` which is negative.
11+
We consider the classic **double integrator** system. The goal is to move from rest at position $0$ to position $1$ (also at rest), using bounded control:
12+
13+
- Initial condition: $x(t_0) = [0, 0]$,
14+
- Final condition: $x(t_f) = [1, 0]$ with $t_f = 0$ fixed,
15+
- Control bounds: $u(t) \in [-1, 1]$,
16+
- The objective is to **maximize $t_0$**, i.e., minimize $-t_0$, which corresponds to a **minimum-time formulation** with fixed arrival time and free departure time.
1217

1318
```@example initial_time
1419
using OptimalControl
1520
16-
tf = 0 # final time
21+
tf = 0 # Final time
22+
x0 = [0, 0] # Initial state
23+
xf = [1, 0] # Final state
1724
1825
@def ocp begin
1926
20-
t0 ∈ R, variable # the initial time is free
27+
t0 ∈ R, variable # t₀ is free and will be optimized
2128
t ∈ [t0, tf], time
2229
x ∈ R², state
2330
u ∈ R, control
2431
25-
-1 ≤ u(t) ≤ 1
32+
-1 ≤ u(t) ≤ 1 # Control bounds
2633
27-
x(t0) == [0, 0]
28-
x(tf) == [1, 0]
34+
x(t0) == x0 # Initial constraint
35+
x(tf) == xf # Final constraint
2936
30-
0.05 ≤ -t0 ≤ Inf # t0 is negative
37+
0.05 ≤ -t0 ≤ Inf # Ensure t₀ is negative and bounded
3138
32-
ẋ(t) == [x₂(t), u(t)]
39+
ẋ(t) == [x₂(t), u(t)] # Dynamics
3340
34-
-t0 → min # or t0 → max
41+
-t0 → min # Equivalent to maximizing t₀
3542
3643
end
3744
nothing # hide
3845
```
3946

40-
# Direct resolution
47+
## Direct solution
4148

42-
Let us solve the problem with a direct method.
49+
We solve the problem using a **direct transcription method**.
4350

4451
```@example initial_time
4552
using NLPModelsIpopt
46-
sol = solve(ocp)
53+
direct_sol = solve(ocp)
4754
nothing # hide
4855
```
4956

50-
The initial time obtained with the direct method is:
57+
The optimal initial time obtained is:
5158

5259
```@example initial_time
53-
t0 = variable(sol)
60+
t0 = variable(direct_sol)
5461
```
5562

56-
We can plot the solution.
63+
We can now plot the solution:
5764

5865
```@example initial_time
5966
using Plots
60-
plot(sol; label="direct", size=(700, 700))
61-
```
62-
63-
## Mathematical verification
64-
65-
Let us compute the solution analytically, thanks to the Pontryagin's Maximum Principle (PMP). First, we define the pseudo-Hamiltonian
66-
67-
```math
68-
H(x, p, u) = p_1 x_2 + p_2 u.
67+
plt = plot(direct_sol; label="direct", size=(700, 700))
6968
```
7069

71-
If $(x(\cdot), u(\cdot), t_0)$ is a solution, then, there exists $(p(\cdot), p^0) \ne 0$, $p^0 \le 0$, satisfying the conditions given by the (PMP). First, since the initial time is free, we have $H(x(t_0), p(t_0), u(t_0)) = -p^0$. Then, the costate $p(\cdot)$ satisfies the adjoint equations
72-
73-
```math
74-
\begin{aligned}
75-
\dot{p}_1(t) &= 0, \\
76-
\dot{p}_2(t) &= -p_1(t),
77-
\end{aligned}
78-
```
70+
## Pontryagin's Maximum Principle
7971

80-
which gives $p_1(t) = \alpha$, $p_2(t) = -\alpha t + \beta$. The maximisation condition from the PMP implies that $u(t) = 1$ if $p_2(t)>0$ and $u(t)=-1$, if $p_2(t)<0$. We switch at time $t_s$ from one control to the other if $p_2(t_s) = 0$, that is when $\beta = \alpha\, t_s$. Note that $p_2$ cannot vanish on an interval of non-empty interior since otherwise, we would have $\alpha = \beta = 0$ and so $H(x(t_0), p(t_0), u(t_0)) = -p^0 = 0$ but $(p(\cdot), p^0) \ne 0$. Hence, there is at most one switching time $t_s$. We can demonstrate that in the setting above, there is exactly one switching time and we switch from $u=-1$ to $u=+1$. Thus, the optimal control is of the form $u(t) = 1$ for $t \in [t_0, t_s)$ and $u(t) = -1$ for $t \in [t_s, 0]$.
72+
The optimal control appears to be **bang-bang**, switching once from $u = 1$ to $u = -1$. This motivates us to apply the Pontryagin's Maximum Principle (PMP) analytically to verify and interpret the solution.
8173

82-
Let us find now $\alpha$, $\beta$, $t_s$ and $t_0$. We can integrate the system: On the interval $[t_0, t_s]$ we have $\dot{x}_2(t) = u(t) = 1$, hence, $x_2(t) = t - t_0$. Since $\dot{x}_1(t) = x_2(t)$, we get
74+
We define the **pseudo-Hamiltonian**:
8375

8476
```math
85-
x_1(t) = \frac{(t - t_0)^2}{2}.
86-
```
87-
88-
At switching time $t_s$ whe have:
89-
90-
```math
91-
x_2(t_s) = t_s - t_0, \quad x_1(t_s) = \frac{(t_s - t_0)^2}{2},
92-
```
93-
94-
that are the initial conditions of the integrations on the interval $[t_s, t_f]$:
95-
96-
```math
97-
\dot{x}_2(t) = u(t) = -1 \quad \Rightarrow \quad x_2(t) = x_2(t_s) - (t - t_s) = 2\, t_s - t_0 - t
98-
```
99-
100-
and
101-
102-
```math
103-
\dot{x}_1(t) = x_2(t) \quad \Rightarrow \quad x_1(t) = x_1(t_s) + \int_{t_s}^t x_2(\sigma)\, \mathrm{d}\sigma =
104-
\frac{(t_s - t_0)^2}{2} + (2\, t_s - t_0) (t - t_s) - \frac{t^2}{2} + \frac{t_s^2}{2}.
77+
H(x, p, u) = p_1 x_2 + p_2 u.
10578
```
10679

107-
The final condition on $x_2$ gives ($t_f = 0$):
80+
Let $(x(\cdot), u(\cdot), t_0)$ be the solution. By PMP, there exists $(p(\cdot), p^0) \neq 0$ with $p^0 \leq 0$ and $p(\cdot)$ absolutely continuous, such that $H(x(t_0), p(t_0), u(t_0)) = -p^0$. The adjoint equations are:
10881

10982
```math
110-
x_2(0) = 0 \quad \Rightarrow \quad t_0 = 2\, t_s.
83+
\dot{p}_1(t) = 0, \quad \dot{p}_2(t) = -p_1(t),
11184
```
11285

113-
The final condition of $x_1$ gives:
86+
so $p_1(t) = \alpha$, $p_2(t) = -\alpha (t-t_0) + \beta$. The maximization condition gives:
11487

11588
```math
116-
x_1(0) = t_s^2 = 1 \quad \Rightarrow \quad t_s = -1.
117-
```
118-
119-
We deduce:
120-
89+
\left\{
90+
\begin{array}{lll}
91+
u(t) = 1 & \text{ if } & p_2(t) > 0, \\[0.5em]
92+
u(t) = -1 & \text{ if } & p_2(t) < 0, \\[0.5em]
93+
u(t) \in [-1, 1] & \text{ if } & p_2(t) = 0.
94+
\end{array}
95+
\right.
96+
```
97+
98+
We assume one switch at $t_1$ with $p_2(t_1) = 0 \Leftrightarrow \beta = \alpha (t_1-t_0)$.
99+
100+
!!! note
101+
If $p_2(t) = 0$ on a nontrivial interval, then $p \equiv 0$ and $p^0 = 0$, violating PMP.
102+
103+
We integrate the system in two phases:
104+
- On $[t_0, t_1]$ with $u = 1$:
105+
$x_2(t) = t - t_0$ and $x_1(t) = \frac{1}{2}(t - t_0)^2$.
106+
- At switch $t = t_1$:
107+
```math
108+
x_2(t_1) = t_1 - t_0, \quad x_1(t_1) = \frac{(t_1 - t_0)^2}{2}
109+
```
110+
- On $[t_1, 0]$ with $u = -1$:
111+
$x_2(t) = x_2(t_1) - (t - t_1) = 2 t_1 - t_0 - t$.
112+
Integrating gives:
113+
```math
114+
x_1(t) = x_1(t_1) + (2t_1 - t_0)(t - t_1) - \frac{t^2}{2} + \frac{t_1^2}{2}
115+
```
116+
117+
Final constraints:
121118
```math
122-
t_0 = 2\, t_s = -2.
119+
x_2(0) = 0 \Rightarrow t_0 = 2 t_1, \quad x_1(0) = 1 \Rightarrow t_1 = -1
123120
```
124121

125-
To get $\alpha$ and $\beta$ we use the initial condition on the pseudo-Hamiltonian and the swtiching condition $p_2(t_s) = 0$. The switching condition gives $\beta = \alpha\, t_s = -\alpha$. The initial condition on the pseudo-Hamiltonian gives:
126-
127-
```math
128-
H(x(t_0), p(t_0), u(t_0)) = p_1(t_0) x_2(t_0) + p_2(t_0) u(t_0) = -\beta = -p^0 = 1.
129-
```
122+
Therefore: $t_0 = -2$.
130123

131-
The dual variable $p^0 \ne 0$ since otherwise $(p(\cdot), p^0) = 0$ which is not. Then, $p^0 < 0$ and we can normalise it to $p^0 = -1$. At the end, we have:
124+
To find $\alpha$ and $\beta$, use:
125+
- Switching: $\beta = \alpha (t_1-t_0) = \alpha$,
126+
- PMP: $H(x(t_0), p(t_0), u(t_0)) = -p^0 = \beta = 1$.
132127

128+
So the full solution is:
133129
```math
134-
\alpha = 1, \quad \beta = -1, \quad t_s = -1, \quad t_0 = -2.
130+
p(t_0) = (1, 1), \quad t_1 = -1, \quad t_0 = -2.
135131
```
136132

137-
## Indirect method
133+
## Indirect method
138134

139-
Let us compute numerically the solution of the PMP using the solution from the direct method. We define the pseudo-Hamiltonian and the flows associated to the two controls.
135+
We now solve the PMP system **numerically** using the direct solution as an initial guess.
140136

141137
```@example initial_time
142-
# pseudo-Hamiltonian
143138
H(x, p, u) = p[1]*x[2] + p[2]*u
144139
145-
# Define flows for u = +1 and u = -1
146-
const u_pos = 1
147-
const u_neg = -1
140+
# Define flows for constant controls
141+
u_pos = 1
142+
u_neg = -1
148143
149144
using OrdinaryDiffEq
150145
f_pos = Flow(ocp, (x, p, t0) -> u_pos)
151146
f_neg = Flow(ocp, (x, p, t0) -> u_neg)
147+
nothing # hide
152148
```
153149

154-
Shooting function adapted for free initial time
155-
```@example initial_time
150+
The **shooting function** encodes:
151+
- continuity of state and costate,
152+
- satisfaction of boundary conditions,
153+
- switching condition ($p_2(t_1) = 0$),
154+
- pseudo-Hamiltonian condition at $t_0$.
156155

156+
```@example initial_time
157157
function shoot!(s, p0, t1, t0)
158-
# Initial conditions
159-
x0 = [0, 0] # fixed initial state
160-
tf = 0 # fixed final time
161-
162-
# The flows (time runs from t0 to tf = 0)
163-
x1, p1 = f_pos(t0, x0, p0, t1) # from t0 to t1
164-
xf, pf = f_neg(t1, x1, p1, tf) # from t1 to tf = 0
165-
166-
# Final control
167-
uf = -1
168-
169-
# Target final state
170-
xf_target = [1, 0]
171-
172-
# Shooting conditions
173-
s[1:2] = xf .- xf_target # reach the target at tf = 0
174-
s[3] = H(xf, pf, uf) - 1 # final condition on pseudo-Hamiltonian
175-
s[4] = H1(x1, p1) # switching condition at t1
158+
x1, p1 = f_pos(t0, x0, p0, t1)
159+
x2, p2 = f_neg(t1, x1, p1, tf)
160+
161+
s[1:2] = x2 - xf # Final state match
162+
s[3] = H(x0, p0, u_pos) - 1 # Hamiltonian normalization
163+
s[4] = p1[2] # Switching condition
176164
end
165+
nothing # hide
177166
```
178167

179-
Get initial guess from direct solution
168+
Get an initial guess from the direct solution:
169+
180170
```@example initial_time
181-
t = time_grid(sol)
182-
x = state(sol)
183-
p = costate(sol)
184-
φ(t) = H1(x(t), p(t)) # switching function
171+
t = time_grid(direct_sol)
172+
x = state(direct_sol)
173+
p = costate(direct_sol)
174+
φ(t) = p(t)[2] # Switching function
185175
186-
p0 = p(t[1]) # initial costate (at t0)
187-
t1 = t[argmin(abs.(φ.(t)))] # switching time
188-
t0 = t[1] # initial time (negative value)
176+
t0 = variable(direct_sol)
177+
p0 = p(t0) # Initial costate
178+
t1 = t[argmin(abs.(φ.(t)))] # Approximate switching time
189179
190-
println("p0 = ", p0)
191-
println("t1 = ", t1)
192-
println("t0 = ", t0)
180+
println("Initial guess:\n\np0 = ", p0, "\nt1 = ", t1, "\nt0 = ", t0)
193181
194-
#Test shooting function at initial guess
195182
s = zeros(4)
196183
shoot!(s, p0, t1, t0)
197184
using LinearAlgebra: norm
198-
println("\nNorm of the shooting function: ‖s‖ = ", norm(s), "\n")
185+
println("\n‖s‖ = ", norm(s))
186+
```
187+
188+
Solve the shooting equations:
199189

200-
# Solve the nonlinear system
190+
```@example initial_time
201191
using NonlinearSolve
202192
203-
# Aggregated function
204193
nle! = (s, ξ, λ) -> shoot!(s, ξ[1:2], ξ[3], ξ[4])
205-
206-
# Initial guess: [p0[1], p0[2], t1, t0]
207194
ξ_guess = [p0..., t1, t0]
208195
prob = NonlinearProblem(nle!, ξ_guess)
209196
210-
# Solve the problem
211-
indirect_sol = solve(prob; show_trace=Val(true), abstol=1e-8, reltol=1e-8)
197+
nle_sol = solve(prob; show_trace=Val(true), abstol=1e-8, reltol=1e-8)
198+
212199
# Extract solution
213-
p0_opt = indirect_sol.u[1:2]
214-
t1_opt = indirect_sol.u[3]
215-
t0_opt = indirect_sol.u[4]
216-
```
200+
p0 = nle_sol.u[1:2]
201+
t1 = nle_sol.u[3]
202+
t0 = nle_sol.u[4]
217203
218-
We can now plot and compare with the direct method.# Compute and plot the optimal trajectory
204+
println("\nIndirect solution:\n\np0 = ", p0, "\nt1 = ", t1, "\nt0 = ", t0)
219205
220-
```@example initial_time
221-
x0 = [0, 0]
222-
tf = 0
206+
shoot!(s, p0, t1, t0)
207+
println("\nFinal shooting residual: ‖s‖ = ", norm(s))
208+
```
223209

224-
# Concatenate flows: u=+1 from t0 to t1, then u=-1 from t1 to tf
225-
f = f_pos * (t1_opt, f_neg)
226-
flow_sol = f((t0_opt, tf), x0, p0_opt; saveat=range(t0_opt, tf, 100))
210+
Compare with the direct solution:
227211

228-
# Plot comparison
229-
plt = plot(sol; label="direct", size=(800, 800))
230-
plot!(plt, flow_sol; label="indirect", color=:red)
231-
```
212+
```@example initial_time
213+
f = f_pos * (t1, f_neg) # Concatenate bang-bang flow
214+
indirect_sol = f((t0, tf), x0, p0; saveat=range(t0, tf, 200))
215+
plot!(plt, indirect_sol; label="indirect", color=:red, linestyle=:dash)
216+
```

0 commit comments

Comments
 (0)