Skip to content
Snippets Groups Projects
Commit c98ffa6b authored by Olivier Cots's avatar Olivier Cots
Browse files

rm tp

parent 34ae4bad
Branches
No related tags found
No related merge requests found
Expand all Hide whitespace changes
Inline Side-by-side
tp/1d-regulator-etu.ipynb deleted 100644 → 0
+ 0
− 595
View file @ 34ae4bad
%% Cell type:markdown id: tags:
![n7](http://cots.perso.enseeiht.fr/figures/inp-enseeiht.png)
# A one dimensional regulator problem
We consider the following optimal control problem:
$$
\left\{
\begin{array}{l}
\min \displaystyle J(x, u) := \displaystyle \int_0^{t_f} x^2(t) \, \mathrm{d}t \\[1.0em]
\dot{x}(t) = f(x(t), u(t)) := \displaystyle u(t), \quad |u(t)| \le 1, \quad t \in [0, t_f] \text{ a.e.}, \\[1.0em]
x(0) = 1, \quad x(t_f) = 1/2.
\end{array}
\right.
$$
To this optimal control problem is associated the stationnary optimization problem
$$
\min_{(x, u)} \{~ x^2 ~ | ~ (x, u) \in \mathrm{R} \times [-1, 1],~ f(x,u) = u = 0\}.
$$
The static solution is thus $(x^*, u^*) = (0, 0)$. This solution may be seen as the static pair $(x, u)$ which minimizes the cost $J(x, u)$ under
the constraint $u \in [-1, 1]$.
It is well known that this problem is what we call a *turnpike* optimal control problem.
Hence, if the final time $t_f$ is long enough the solution is of the following form:
starting from $x(0)=1$, reach as fast as possible the static solution, stay at the static solution as long as possible before reaching
the target $x(t_f)=1/2$. In this case, the optimal control would be
$$
u(t) = \left\{
\begin{array}{lll}
-1 & \text{if} & t \in [0, t_1], \\[0.5em]
\phantom{-}0 & \text{if} & t \in (t_1, t_2], \\[0.5em]
+1 & \text{if} & t \in (t_2, t_f],
\end{array}
\right.
$$
with $0 < t_1 < t_2 < t_f$. We say that the control is *Bang-Singular-Bang*. A Bang arc corresponds to $u \in \{-1, 1\}$ while a singular control corresponds to $u \in (-1, 1)$. Since the optimal control law is discontinuous, then to solve this optimal control problem by indirect shooting methods and find the *switching times* $t_1$ and $t_2$, we need to give the structure to the algorithm. Hence, assuming we do not know the optimal structure, we need an alternative. We thus consider a direct method which consists roughly speaking in a full discretization of the state and control variables on a given time steps grid.
<div class="alert alert-warning">
**Main goals**
1. Retrieve the optimal structure by a direct method.
2. Find the switching times $t_1$ and $t_2$ by indirect shooting.
</div>
%% Cell type:markdown id: tags:
## Direct method
%% Cell type:code id: tags:
``` julia
# For direct methods
using JuMP, Ipopt
# To plot solutions
using Plots
```
%% Cell type:code id: tags:
``` julia
# See the following links to get examples:
# https://ct.gitlabpages.inria.fr/gallery/goddard-j/goddard.html
# https://ct.gitlabpages.inria.fr/gallery/turnpike/2d/turnpike.html
# Create JuMP model, using Ipopt as the solver
sys = Model(optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 5))
set_optimizer_attribute(sys,"tol",1e-8)
set_optimizer_attribute(sys,"constr_viol_tol",1e-6)
set_optimizer_attribute(sys,"max_iter",1000)
#######
####### DEBUT A MODIFIER
#######
####### Les ... sont a remplacer !
#######
####### Attention, on ecrit le probleme sous la forme d'un probleme de Mayer
#######
# Parameters
t0 = 0. # initial time
tf = 0. # final time
c0 = 0. # Initial cost
x0 = 0. # Initial position
xf = 0. # Final position
N = 50 # Grid size
Δt = (tf-t0)/N # Time step
@variables(sys, begin
... # Cost
... # Position
... # Control
end)
# Objective
@objective(sys, Min, ...)
# Boundary constraints
@constraints(sys, begin
con_c0, ... # constraint on inital cost
con_x0, ... # constraint on inital x
con_xf, ... # constraint on final x
end)
# Dynamics with Crank-Nicolson scheme
@NLconstraints(sys, begin
con_dc[j=1:N], ... # differential constraint on the dynamics of the cost
con_dx[j=1:N], ... # differential constraint on the dynamics of the state x
end)
#######
####### FIN A MODIFIER
#######
```
%% Cell type:code id: tags:
``` julia
# Solve for the control and state
println("Solving...")
status = optimize!(sys)
println()
# Display results
if termination_status(sys) == MOI.OPTIMAL
println(" Solution is optimal")
elseif termination_status(sys) == MOI.LOCALLY_SOLVED
println(" (Local) solution found")
elseif termination_status(sys) == MOI.TIME_LIMIT && has_values(sys)
println(" Solution is suboptimal due to a time limit, but a primal solution is available")
else
error(" The model was not solved correctly.")
end
println(" objective value = ", objective_value(sys))
println()
# Retrieves values (including duals)
c = value.(c)[:]
x = value.(x)[:]
u = value.(u)[:]
t = (1:N+1) * value.(Δt)
pc0 = dual(con_c0)
px0 = dual(con_x0)
pxf = dual(con_xf)
if(pc0*dual(con_dc[1])<0); pc0 = -pc0; end
if(px0*dual(con_dx[1])<0); px0 = -px0; end
if(pxf*dual(con_dx[N])<0); pxf = -pxf; end
if (pc0 > 0) # Sign convention according to Pontryagin Maximum Principle
sign = -1.0
else
sign = 1.0
end
pc = [ dual(con_dc[i]) for i in 1:N ]
px = [ dual(con_dx[i]) for i in 1:N ]
pc = sign * [pc0; pc[1:N]]
px = sign * [px0; (px[1:N-1]+px[2:N])/2; pxf]; # We add the multiplier from the limit conditions
```
%% Output
Solving...
******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
Ipopt is released as open source code under the Eclipse Public License (EPL).
For more information visit https://github.com/coin-or/Ipopt
******************************************************************************
This is Ipopt version 3.13.4, running with linear solver mumps.
NOTE: Other linear solvers might be more efficient (see Ipopt documentation).
Number of nonzeros in equality constraint Jacobian...: 403
Number of nonzeros in inequality constraint Jacobian.: 0
Number of nonzeros in Lagrangian Hessian.............: 100
Total number of variables............................: 153
variables with only lower bounds: 0
variables with lower and upper bounds: 51
variables with only upper bounds: 0
Total number of equality constraints.................: 103
Total number of inequality constraints...............: 0
inequality constraints with only lower bounds: 0
inequality constraints with lower and upper bounds: 0
inequality constraints with only upper bounds: 0
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
0 0.0000000e+00 1.00e+00 0.00e+00 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0
1 0.0000000e+00 3.97e-02 1.23e-01 -1.7 1.00e+00 - 7.83e-01 1.00e+00h 1
2 1.0516638e+00 7.88e-05 1.93e-02 -1.7 1.05e+00 - 9.51e-01 1.00e+00h 1
3 5.9684864e-01 2.59e-03 7.32e-03 -2.5 9.01e-01 - 7.46e-01 1.00e+00f 1
4 4.3476735e-01 1.70e-03 4.68e-03 -3.8 4.27e-01 - 7.74e-01 1.00e+00h 1
5 3.9437257e-01 6.32e-04 5.67e-04 -3.8 3.27e-01 - 9.50e-01 1.00e+00h 1
6 3.7872229e-01 2.26e-04 6.96e-04 -5.7 1.62e-01 - 8.28e-01 1.00e+00h 1
7 3.7635925e-01 4.48e-05 5.14e-04 -5.7 3.42e-01 - 6.87e-01 1.00e+00h 1
8 3.7566913e-01 8.63e-06 2.06e-04 -5.7 3.96e-01 - 7.23e-01 1.00e+00h 1
9 3.7550672e-01 2.24e-06 5.23e-05 -5.7 3.66e-01 - 8.08e-01 1.00e+00h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
10 3.7547509e-01 4.68e-07 1.11e-16 -5.7 2.70e-01 - 1.00e+00 1.00e+00h 1
11 3.7540196e-01 2.52e-07 3.34e-06 -8.6 1.90e-01 - 8.57e-01 1.00e+00h 1
12 3.7540050e-01 7.12e-08 1.28e-06 -8.6 2.52e-01 - 8.33e-01 1.00e+00h 1
13 3.7540036e-01 2.43e-08 1.37e-07 -8.6 2.05e-01 - 9.60e-01 1.00e+00h 1
14 3.7540039e-01 5.99e-09 1.11e-16 -8.6 8.66e-02 - 1.00e+00 1.00e+00h 1
15 3.7540034e-01 9.81e-10 1.11e-16 -9.0 3.11e-02 - 1.00e+00 1.00e+00h 1
Number of Iterations....: 15
(scaled) (unscaled)
Objective...............: 3.7540033633411968e-01 3.7540033633411968e-01
Dual infeasibility......: 1.1102230246251565e-16 1.1102230246251565e-16
Constraint violation....: 9.8117847180390072e-10 9.8117847180390072e-10
Complementarity.........: 2.6817939297375691e-09 2.6817939297375691e-09
Overall NLP error.......: 2.6817939297375691e-09 2.6817939297375691e-09
Number of objective function evaluations = 16
Number of objective gradient evaluations = 16
Number of equality constraint evaluations = 16
Number of inequality constraint evaluations = 0
Number of equality constraint Jacobian evaluations = 16
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations = 15
Total CPU secs in IPOPT (w/o function evaluations) = 2.366
Total CPU secs in NLP function evaluations = 2.017
EXIT: Optimal Solution Found.
(Local) solution found
objective value = 0.3754003363341197
%% Cell type:code id: tags:
``` julia
#c_plot = plot(t, c, xlabel = "t", ylabel = "c", legend = false, fmt = :png)
x_plot = plot(t, x, xlabel = "t", ylabel = "x", legend = false, fmt = :png)
u_plot = plot(t, u, xlabel = "t", ylabel = "u", legend = false, fmt = :png, size=(800,400), linetype=:steppre)
#pc_plot = plot(t, pc, xlabel = "t", ylabel = "pc", legend = false, fmt = :png)
px_plot = plot(t, px, xlabel = "t", ylabel = "px", legend = false, fmt = :png)
#display(plot(c_plot, x_plot, pc_plot, px_plot, layout = (2,2), size=(800,400)))
display(plot(x_plot, px_plot, layout = (1,2), size=(800,400)))
display(u_plot)
```
%% Output
%% Cell type:markdown id: tags:
## Indirect method
%% Cell type:code id: tags:
``` julia
# For indirect methods
using DifferentialEquations, NLsolve, ForwardDiff
```
%% Cell type:code id: tags:
``` julia
# Function to get the flow of a Hamiltonian system
function Flow(hv)
function rhs!(dz, z, dummy, t)
n = size(z, 1)÷2
dz[:] = hv(z[1:n], z[n+1:2*n])
end
function f(tspan, x0, p0; abstol=1e-12, reltol=1e-12, saveat=0.01)
z0 = [ x0 ; p0 ]
ode = ODEProblem(rhs!, z0, tspan)
sol = DifferentialEquations.solve(ode, Tsit5(), abstol=abstol, reltol=reltol, saveat=saveat)
return sol
end
function f(t0, x0, p0, tf; abstol=1e-12, reltol=1e-12, saveat=[])
sol = f((t0, tf), x0, p0, abstol=abstol, reltol=reltol, saveat=saveat)
n = size(x0, 1)
return sol[1:n, end], sol[n+1:2*n, end]
end
return f
end;
```
%% Cell type:markdown id: tags:
The pseudo-Hamiltonian is
$$
H(x, p, u) = pu - x^2
$$
and the pseudo-Hamiltonian system is
$$
\vec{H}(x, p, u) = (\partial_p H(x, p, u), -\partial_x H(x, p, u)) = (u, 2x).
$$
The maximizing control is given by
$$
u(t)~ \left\{~
\begin{array}{lll}
= -1 & \text{if} & p(t) < 0, \\[0.5em]
\in [-1, 1] & \text{if} & p(t) = 0, \\[0.5em]
= +1 & \text{if} & p(t) > 0,
\end{array}
\right.
$$
A *singular arc* is a part of an extremal such that $p(t) = 0$. Here, we have a singular arc for $t \in [t_1, t_2]$.
On the interval $[t_1, t_2]$, we have $p(t)=0$. Hence, along this arc $\dot{p}(t)=2x(t)=0$. Derivating once again,
we have $\dot{x}(t) = u(t) = 0$. The singular control is thus the null control.
%% Cell type:code id: tags:
``` julia
# Pseudo-Hamiltonian system
function hv(x, p, u)
r = zeros(eltype(x), 2)
r[1] = u
r[2] = 2.0*x[1]
return r
end
hv_min(x, p) = hv(x, p, -1.0)
hv_max(x, p) = hv(x, p, 1.0)
hv_sin(x, p) = hv(x, p, 0.0)
fmin = Flow(hv_min)
fmax = Flow(hv_max)
fsin = Flow(hv_sin);
```
%% Cell type:markdown id: tags:
We have determined that the optimal control follows the strategy:
$$
u(t) = \left\{
\begin{array}{lll}
-1 & \text{if} & t \in [0, t_1], \\[0.5em]
\phantom{-}0 & \text{if} & t \in (t_1, t_2], \\[0.5em]
+1 & \text{if} & t \in (t_2, t_f],
\end{array}
\right.
$$
with $0 < t_1 < t_2 < t_f=2$.
We define the following shooting function:
$$
S(p_0, t_1, t_2) :=
\begin{pmatrix}
x(t_1, t_0, x_0, p_0, u_-) \\
p(t_1, t_0, x_0, p_0, u_-) \\
x(t_f, t_2, x_2, p_2, u_+) - 1/2
\end{pmatrix},
$$
where $z_2 := (x_2, p_2) = z(t_2, t_1, x_1, p_1, u_0)$, $z_1 := (x_1, p_1) = z(t_1, t_0, x_0, p_0, u_-)$ and where z(t, s, a, b, u) is the solution at time $t$ of the Hamiltonian system associated to the control u starting at time $s$ at the initial condition $z(s) = (a,b)$.
We have introduced the notation $u_-$ for $u\equiv -1$, $u_0$ for $u\equiv 0$ and $u_+$ for $u\equiv +1$.
%% Cell type:code id: tags:
``` julia
# See the following link to get more details about this problem:
# https://ct.gitlabpages.inria.fr/gallery/shooting_tutorials/multiple_shooting_homotopy.html
#######
####### DEBUT A MODIFIER
#######
####### Les ... sont a remplacer !
#######
# Shooting function
x0 = [x0]
function shoot(p0, t1, t2)
# integration
x1, p1 = ...
x2, p2 = ...
x3, p3 = ...
# conditions
s = zeros(eltype(p0), 3)
#
s[1] = ...
s[2] = ...
s[3] = ...
return s
end;
#######
####### FIN A MODIFIER
#######
#######
```
%% Cell type:code id: tags:
``` julia
# Find the initial guess
η = 1e-3
tsin = t[ abs.(px) .≤ η ]
p0 = px[1]
t1 = min(tsin...)
t2 = max(tsin...)
y = [ p0 ; t1 ; t2]
println("Initial guess:\n", y)
```
%% Output
Initial guess:
[-1.0000353948050356, 1.04, 1.56]
%% Cell type:code id: tags:
``` julia
# Jacobian of the shooting function
foo(y) = shoot(y[1], y[2], y[3])
jfoo(y) = ForwardDiff.jacobian(foo, y)
# Solve shoot(p0, t1, t2) = 0.
nl_sol = nlsolve(foo, jfoo, y; xtol=1e-8, method=:trust_region, show_trace=true);
# Retrieves solution
if converged(nl_sol)
p0 = nl_sol.zero[1]
t1 = nl_sol.zero[2]
t2 = nl_sol.zero[3];
println("\nFinal solution:\n", nl_sol.zero)
else
error("Not converged")
end
```
%% Output
Iter f(x) inf-norm Step 2-norm
------ -------------- --------------
0 1.000000e-01 NaN
1 1.600000e-03 7.212800e-02
2 4.440892e-16 1.600000e-03
Final solution:
[-1.0000000000000013, 1.0000000000000004, 1.5000000000000002]
%% Cell type:code id: tags:
``` julia
# Plots
ode_sol = fmin((t0, t1), x0, p0)
tt0 = ode_sol.t
xx0 = [ ode_sol[1, j] for j in 1:size(tt0, 1) ]
pp0 = [ ode_sol[2, j] for j in 1:size(tt0, 1) ]
uu0 = -ones(size(tt0, 1))
ode_sol = fsin((t1, t2), xx0[end], pp0[end])
tt1 = ode_sol.t
xx1 = [ ode_sol[1, j] for j in 1:size(tt1, 1) ]
pp1 = [ ode_sol[2, j] for j in 1:size(tt1, 1) ]
uu1 = zeros(size(tt1, 1))
ode_sol = fmax((t2, tf), xx1[end], pp1[end])
tt2 = ode_sol.t
xx2 = [ ode_sol[1, j] for j in 1:size(tt2, 1) ]
pp2 = [ ode_sol[2, j] for j in 1:size(tt2, 1) ]
uu2 = ones(size(tt2, 1))
t_shoot = [ tt0 ; tt1 ; tt2 ]
x_shoot = [ xx0 ; xx1 ; xx2 ]
p_shoot = [ pp0 ; pp1 ; pp2 ]
u_shoot = [ uu0 ; uu1 ; uu2 ]
x_plot = plot(t_shoot, x_shoot, xlabel = "t", ylabel = "x", legend = false, fmt = :png)
u_plot = plot(t_shoot, u_shoot, xlabel = "t", ylabel = "u", legend = false, fmt = :png, size=(800,400), linetype=:steppre)
px_plot = plot(t_shoot, p_shoot, xlabel = "t", ylabel = "px", legend = false, fmt = :png)
display(plot(x_plot, px_plot, layout = (1,2), size=(800,400)))
display(u_plot)
```
tp/animation-etu.ipynb deleted 100644 → 0
+ 0
− 638
View file @ 34ae4bad
This diff is collapsed.
tp/ode-etu.ipynb deleted 100644 → 0
+ 0
− 8940
View file @ 34ae4bad
This diff is collapsed.
tp/sandbox.ipynb deleted 100644 → 0
+ 0
− 96
View file @ 34ae4bad
%% Cell type:markdown id: tags:
![n7](http://cots.perso.enseeiht.fr/figures/inp-enseeiht.png)
# An empty notebook to play with
%% Cell type:markdown id: tags:
## Direct method
%% Cell type:code id: tags:
``` julia
# For direct methods
using JuMP, Ipopt
# To plot solutions
using Plots
```
%% Cell type:markdown id: tags:
## Indirect method
%% Cell type:code id: tags:
``` julia
# For indirect methods
using DifferentialEquations, NLsolve, ForwardDiff
```
%% Cell type:code id: tags:
``` julia
# Function to get the flow of a Hamiltonian system
function Flow(hv)
function rhs!(dz, z, dummy, t)
n = size(z, 1)÷2
dz[:] = hv(z[1:n], z[n+1:2*n])
end
function f(tspan, x0, p0; abstol=1e-12, reltol=1e-12, saveat=0.01)
z0 = [ x0 ; p0 ]
ode = ODEProblem(rhs!, z0, tspan)
sol = DifferentialEquations.solve(ode, Tsit5(), abstol=abstol, reltol=reltol, saveat=saveat)
return sol
end
function f(t0, x0, p0, tf; abstol=1e-12, reltol=1e-12, saveat=[])
sol = f((t0, tf), x0, p0, abstol=abstol, reltol=reltol, saveat=saveat)
n = size(x0, 1)
return sol[1:n, end], sol[n+1:2*n, end]
end
return f
end;
```
tp/simple-shooting-etu.ipynb deleted 100644 → 0
+ 0
− 176
View file @ 34ae4bad
%% Cell type:markdown id: tags:
![tata](http://cots.perso.enseeiht.fr/figures/inp-enseeiht.png)
# Tir simple indirect : une introduction
Le but est de résoudre par du tir simple indirect, deux problèmes simples de contrôle optimal dont le contrôle maximisant est lisse.
Le premier problème possède des conditions aux limites simples tandis que le second aura des conditions de transversalités.
**Packages.**
%% Cell type:code id: tags:
``` julia
using DifferentialEquations, Plots, NLsolve, ForwardDiff
```
%% Cell type:markdown id: tags:
## Exercice 1
On considère le problème de contrôle optimal suivant :
$$
\left\{
\begin{array}{l}
\displaystyle \min\, \frac{1}{2} \int_0^{1} u^2(t) \, \mathrm{d}t \\[1.0em]
\dot{x}(t) = -x(t) + \alpha x^2(t) + u(t), \quad u(t) \in \mathrm{R}, \quad t \in [0, 1] \text{ p.p.}, \\[1.0em]
x(0) = -1, \quad x(1) = 0.
\end{array}
\right.
$$
Pour $\alpha \in \{0, 1.5\}$, faire :
1. Afficher le flot du système hamiltonien associé et afficher 5 fronts aux temps $0$, $0.25$, $0.5$, $0.75$ et $1$. On pourra utiliser le TP précédent.
2. Coder la fonction de tir
$$
S(p_0) = \pi(z(t_f, x_0, p_0)) - x_f,
$$
où $t_f = 1$, $x_0=-1$, $x_f = 0$. On pourra se référer au polycopié de cours, cf. chapitre 5, pour les notations et la définition de la fonction de tir.
3. Afficher la fonction de tir pour $p_0 \in [0, 1]$.
4. Résoudre l'équation $S(p_0) = 0$ et tracer l'extrémale dans le plan de phase, cf. question 1.
Remarques.
* Voici quelques appels de fonctions pouvant aider au calcul du flot d'un système hamiltonien~:
```julia
# Hamiltonian system
function hv(x, p)
r = [0; 0]
return r
end
# Makes flow from a Hamiltonian system
f = Flow(hv)
```
* Voici quelques appels de fonctions pouvant aider à la résolution des équations de tir~:
```julia
# jacobienne de la fonction de tir
jshoot(p0) = ForwardDiff.jacobian(shoot, p0)
# résolution de shoot(p0) = 0. L'itéré initial est appelé p0_guess
sol = nlsolve(shoot, jshoot, p0_guess; xtol=1e-8, method=:trust_region, show_trace=true);
```
%% Cell type:code id: tags:
``` julia
function Flow(hv)
function rhs!(dz, z, dummy, t)
n = size(z, 1)÷2
dz[:] = hv(z[1:n], z[n+1:2*n])
end
function f(tspan, x0, p0; abstol=1e-12, reltol=1e-12, saveat=[])
z0 = [ x0 ; p0 ]
ode = ODEProblem(rhs!, z0, tspan)
sol = solve(ode, Tsit5(), abstol=abstol, reltol=reltol, saveat=saveat)
return sol
end
function f(t0, x0, p0, tf; abstol=1e-12, reltol=1e-12, saveat=[])
sol = f((t0, tf), x0, p0, abstol=abstol, reltol=reltol, saveat=saveat)
n = size(x0, 1)
return sol[1:n, end], sol[n+1:2*n, end]
end
return f
end
```
%% Cell type:code id: tags:
``` julia
```
%% Cell type:markdown id: tags:
## Exercice 2
On considère le problème de contrôle optimal suivant :
$$
\left\{
\begin{array}{l}
\displaystyle \min\, \frac{1}{2} \int_0^{t_f} ||u(t)||^2 \, \mathrm{d}t \\[1.0em]
\dot{x}(t) = (0, \alpha) + u(t), \quad u(t) \in \mathrm{R}^2, \quad t \in [0, t_f] \text{ p.p.}, \\[1.0em]
x(0) = a, \quad x(t_f) \in b + \mathrm{R} v,
\end{array}
\right.
$$
avec $a = (0,0)$, $b = (1,0)$, $v = (0, 1)$ et $t_f = 1$.
Pour différente valeur de $\alpha$, faire :
1. Résoudre le problème de contrôle optimal par du tir simple.
2. Tracer la trajectoire solution dans le plan $(x_1, x_2)$.
3. Tracer pour la solution, attaché au point $x(t_f)$, les vecteurs $\dot{x}(t_f)$ et $p(t_f)$. Commentaire ?
%% Cell type:code id: tags:
``` julia
```
tp/transfert-orbital-temps-min-etu.ipynb deleted 100644 → 0
+ 0
− 447
View file @ 34ae4bad
%% Cell type:markdown id: tags:
![n7](http://cots.perso.enseeiht.fr/figures/inp-enseeiht.png)
# Transfert orbital en temps minimal
Ce sujet (ou TP-projet) est à rendre (voir la date sur moodle et les modalités du rendu) et sera évalué pour faire partie de la note finale de la matière Contrôle Optimal.
----
On considère le problème du transfert d'un satellite d'une orbite initiale à l'orbite géostationnaire à temps minimal. Ce problème s'écrit comme un problème de contrôle optimal sous la forme
$$
\left\lbrace
\begin{array}{l}
\min J(x, u, t_f) = t_f \\[1.0em]
\ \ \dot{x}_{1}(t) = ~ x_{3}(t) \\
\ \ \dot{x}_{2}(t) = ~ x_{4}(t) \\
\ \ \dot{x}_{3}(t) = -\dfrac{\mu\, x_{1}(t)}{r^{3}(t)} + u_{1}(t) \\
\ \ \dot{x}_{4}(t) = -\dfrac{\mu\, x_{2}(t)}{r^{3}(t)}+u_{2}(t), ~~ ||u(t)|| \leq \gamma_\mathrm{max}, ~~ t \in [0,t_f] ~~ \text{p. p}, ~~ u(t) = (u_1(t),u_2(t)), \\[1.0em]
\ \ x_{1}(0)=x_{0,1}, ~~ x_{2}(0)=x_{0,2}, ~~ x_{3}(0)=x_{0,3}, ~~ x_{4}(0)=x_{0,4}, \\
\ \ r^2(t_f) = r_{f}^2, ~~ x_{3}(t_f)=-\sqrt{\dfrac{\mu}{r_{f}^{3}}}x_{2}(t_f), ~~ x_{4} (t_f)= \sqrt{\dfrac{\mu}{r_{f}^{3}}}x_{1}(t_f), \\
\end{array}
\right.
$$
avec $r(t)=\sqrt{x_{1}^{2}(t)+x_{2}^{2}(t)}$.
Les unités choisies sont le kilomètre pour les distances et l'heure pour les temps. On donne les paramètres suivants :
$$
\mu=5.1658620912 \times 10^{12} \ \mathrm{km}^{3}.\mathrm{h}^{-2}, \quad r_{f} = 42165 \ \mathrm{km}.
$$
Le paramètre $\gamma_\mathrm{max}$ dépend de la poussée maximale $F_\mathrm{max}$ suivant la relation :
$$
\gamma_\mathrm{max} = \frac{F_\mathrm{max}\times 3600^2}{m}
$$
où m est la masse du satellite qu'on fixe à $m=2000\ \mathrm{kg}$.
%% Cell type:markdown id: tags:
## Résolution via du tir simple indirect
%% Cell type:code id: tags:
``` julia
# Packages utiles
using DifferentialEquations, NLsolve, ForwardDiff, Plots, LinearAlgebra
```
%% Cell type:code id: tags:
``` julia
# Les constantes du pb
x0 = [-42272.67, 0, 0, -5796.72] # état initial
μ = 5.1658620912*1e12
rf = 42165.0 ;
F_max = 100.0
γ_max = F_max*3600.0^2/(2000.0*10^3)
t0 = 0.0
rf3 = rf^3 ;
α = sqrt(μ/rf3);
```
%% Cell type:code id: tags:
``` julia
# Function to get the flow of a Hamiltonian system
function Flow(hv)
function rhs!(dz, z, dummy, t)
n = size(z, 1)÷2
dz[:] = hv(z[1:n], z[n+1:2*n])
end
function f(tspan, x0, p0; abstol=1e-12, reltol=1e-12, saveat=0.1)
z0 = [ x0 ; p0 ]
ode = ODEProblem(rhs!, z0, tspan)
sol = DifferentialEquations.solve(ode, Tsit5(), abstol=abstol, reltol=reltol, saveat=saveat)
return sol
end
function f(t0, x0, p0, t; abstol=1e-12, reltol=1e-12, saveat=[])
sol = f((t0, t), x0, p0; abstol=abstol, reltol=reltol, saveat=saveat)
n = size(x0, 1)
return sol[1:n, end], sol[n+1:2*n, end]
end
return f
end;
```
%% Cell type:markdown id: tags:
<div class="alert alert-info">
**_Question 1:_**
* Donner le pseudo-hamiltonien associé au problème de contrôle optimal.
* Calculer le contrôle maximisant dans le cas normal ($p^0=-1$). On supposera que $(p_3, p_4)\neq (0,0)$).
* Donner le pseudo système hamiltonien
$$
\vec{H}(x, p, u) = \left(\frac{\partial H}{\partial p}(x, p, u),
-\frac{\partial H}{\partial x}(x, p, u) \right).
$$
</div>
%% Cell type:markdown id: tags:
**Réponse** (double cliquer pour modifier le texte)
%% Cell type:markdown id: tags:
<div class="alert alert-info">
**_Question 2:_** Compléter le code suivant pour calculer la trajectoire optimale et l'afficher.
</div>
%% Cell type:code id: tags:
``` julia
#####
##### A COMPLETER
# Contrôle maximisant
function control(p)
u = zeros(eltype(p),2)
u[1] = ...
u[2] = ...
return u
end;
# Hamiltonien maximisé
function hfun(x, p)
h = ...
return h
end
# Système hamiltonien
function hv(x, p)
n = size(x, 1)
hv = zeros(eltype(x), 2*n)
...
return hv
end
f = Flow(hv);
```
%% Cell type:markdown id: tags:
On note
$
\alpha := \sqrt{\frac{\mu}{r_f^3}}.
$
La condition terminale peut se mettre sous la forme $c(x(t_f)) = 0$, avec $c \colon \mathbb{R}^4 \to \mathbb{R}^3$
donné par
$$
c(x) := (r^2(x) - r_f^2, ~~ x_{3} + \alpha\, x_{2}, ~~ x_{4} - \alpha\, x_{1}).
$$
Le temps final étant libre, on a la condition au temps final
$$
H(x(t_f), p(t_f), u(t_f)) = -p^0 = 1. \quad \text{(cas normal)}
$$
De plus, la condition de transversalité
$$
p(t_f) = c'(x(t_f))^T \lambda, ~~ \lambda \in \mathbb{R}^3,
$$
conduit à
$$
\Phi(x(t_f), p(t_f)) := x_2(t_f) \Big( p_1(t_f) + \alpha\, p_4(t_f) \Big) - x_1(t_f) \Big( p_2(t_f) - \alpha\, p_3(t_f) \Big) = 0.
$$
En considérant les conditions aux limites, la condition finale sur le pseudo-hamiltonien et la condition de transversalité, la fonction de tir simple est donnée par
\begin{equation*}
\begin{array}{rlll}
S \colon & \mathbb{R}^5 & \longrightarrow & \mathbb{R}^5 \\
& (p_0, t_f) & \longmapsto &
S(p_0, t_f) := \begin{pmatrix}
c(x(t_f, x_0, p_0)) \\[0.5em]
\Phi(z(t_f, x_0, p_0)) \\[0.5em]
H(z(t_f, x_0, p_0), u(z(t_f, x_0, p_0))) - 1
\end{pmatrix}
\end{array}
\end{equation*}
où $z(t_f, x_0, p_0)$ est la solution au temps de $t_f$ du pseudo système hamiltonien bouclé par
le contrôle maximisant, partant au temps $t_0$ du point $(x_0, p_0)$. On rappelle que l'on note
$z=(x, p)$.
%% Cell type:code id: tags:
``` julia
#####
##### A COMPLETER
# Fonction de tir
function shoot(p0, tf)
s = zeros(eltype(p0), 5)
...
return s
end;
```
%% Cell type:code id: tags:
``` julia
# Itéré initial
y_guess = [1.0323e-4, 4.915e-5, 3.568e-4, -1.554e-4, 13.4] # pour F_max = 100N
# Jacobienne de la fonction de tir
foo(y) = shoot(y[1:4], y[5])
jfoo(y) = ForwardDiff.jacobian(foo, y)
# Résolution de shoot(p0, tf) = 0
nl_sol = nlsolve(foo, jfoo, y_guess; xtol=1e-8, method=:trust_region, show_trace=true);
# On récupère la solution si convergence
if converged(nl_sol)
p0_sol_100 = nl_sol.zero[1:4]
tf_sol_100 = nl_sol.zero[5]
println("\nFinal solution:\n", nl_sol.zero)
else
error("Not converged")
end
```
%% Cell type:code id: tags:
``` julia
# Fonction d'affichage d'une solution
function plot_solution(p0, tf)
# On trace l'orbite de départ et d'arrivée
gr(dpi=300, size=(500,400), thickness_scaling=1)
r0 = norm(x0[1:2])
v0 = norm(x0[3:4])
a = 1.0/(2.0/r0-v0*v0/μ)
t1 = r0*v0*v0/μ - 1.0;
t2 = (x0[1:2]'*x0[3:4])/sqrt(a*μ);
e_ellipse = norm([t1 t2])
p_orb = a*(1-e_ellipse^2);
n_theta = 101
Theta = range(0, stop=2*pi, length=n_theta)
X1_orb_init = zeros(n_theta)
X2_orb_init = zeros(n_theta)
X1_orb_arr = zeros(n_theta)
X2_orb_arr = zeros(n_theta)
# Orbite initiale
for i in 1:n_theta
theta = Theta[i]
r_orb = p_orb/(1+e_ellipse*cos(theta));
X1_orb_init[i] = r_orb*cos(theta);
X2_orb_init[i] = r_orb*sin(theta);
end
# Orbite d'arrivée
for i in 1:n_theta
theta = Theta[i]
X1_orb_arr[i] = rf*cos(theta) ;
X2_orb_arr[i] = rf*sin(theta);
end;
# Calcul de la trajectoire
ode_sol = f((t0, tf), x0, p0)
t = ode_sol.t
x1 = [ode_sol[1, j] for j in 1:size(t, 1) ]
x2 = [ode_sol[2, j] for j in 1:size(t, 1) ]
u = zeros(2, length(t))
for j in 1:size(t, 1)
u[:,j] = control(ode_sol[5:8, j])
end
px = plot(x1, x2, color = :red, legend = false,
border=:none, size = (800,400), aspect_ratio=:equal)
plot!(px, X1_orb_init, X2_orb_init, color = :green)
plot!(px, X1_orb_arr, X2_orb_arr, color = :blue)
plot!(px, [x0[1]], [x0[2]], seriestype=:scatter, color = :red)
plot!(px, [0.0], [0.0], color = :blue, seriestype=:scatter, markersize = 20)
pu1 = plot(t, u[1,:], color = :red, xlabel = "t", ylabel = "u1", legend = false)
pu2 = plot(t, u[2,:], color = :red, xlabel = "t", ylabel = "u2", legend = false)
display(plot(pu1, pu2, layout = (1,2), size = (800,400)))
display(px)
end;
```
%% Cell type:code id: tags:
``` julia
# On affiche la solution pour Fmax = 100
plot_solution(p0_sol_100, tf_sol_100)
```
%% Cell type:markdown id: tags:
# Etude du temps de transfert en fonction de la poussée maximale
<div class="alert alert-info">
**_Question 3:_**
* A l'aide de ce que vous avez fait précédemment, déterminer $t_f$ pour
$$
F_\mathrm{max} \in \{100, 90, 80, 70, 60, 50, 40, 30, 20 \}.
$$
* Tracer ensuite $t_f$ en fonction de $F_\mathrm{max}$ et commenter le résultat.
</div>
%% Cell type:code id: tags:
``` julia
# Les différentes valeurs de poussées
F_max_span = range(100, stop=20, length=11);
γ_max_span = [F_max_span[j]*3600^2/(2000*10^3) for j in 1:size(F_max_span,1)];
```
%% Cell type:code id: tags:
``` julia
#####
##### A COMPLETER
# Solution calculée précédemment
y_guess = [p0_sol_100; [tf_sol_100]]
tf_sols = zeros(length(γ_max_span))
p0_sols = zeros(4, length(γ_max_span))
...
```
%% Cell type:code id: tags:
``` julia
# Affichage de tf en fonction de la poussée maximale
plot(F_max_span, tf_sols, aspect_ratio=:equal, legend=false, xlabel="Fmax", ylabel="tf")
```
%% Cell type:code id: tags:
``` julia
# Plots sol pour F_max = 20N
plot_solution(p0_sols[:, end], tf_sols[end])
```
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment