Código
using OptimalControl
using NLPModelsIpopt
using DataFrames5 Problema Original
5.1 Formulación
El problema se enfoca en maximizar la altitud final de un cohete en ascenso vertical, considerando los efectos del arrastre atmosférico y el inverso cuadrado del campo gravitacional, tomando el empuje (\(T\)) como el único control. Los parámetros empleados corresponden a los presentados por Seywald y Cliff (1993), mientras que los datos relativos al vehículo y su aerodinámica provienen de Zlatskii y Kiforenko (1975). Dichos datos describen las características del misil Soviet SA-2 surface-to-air missile, NATO.
A continuación, se presenta la formulación matemática del problema (vista como un problema de Cálculo Variacional), basada en los parámetros y datos previamente descritos.
\[ \text{max } h(t_f) \tag{5.1}\]
sujeto a \[ \begin{split} \dot{h} &= v\\ \dot{v} &= \frac{T-D}{m}-\frac{1}{h^2}\\ \dot{m} &= -\frac{T}{c} \end{split} \tag{5.2}\]
el control (Thrust) \[ T\in \left[0,T_{max}\right] \tag{5.3}\]
y las condiciones de frontera \[ \begin{split} r(0) &= 1\\ v(0) &= 0 \\ m(0) &= 1 \\ m(t_f) &= m_f \end{split} \tag{5.4}\]
5.2 Resultados Numéricos
La existencia de la solución es clara debido a las condiciones del problema y lo discutido en Capítulo 4. Resolver el problema de Goddard implica determinar la trayectoria óptima \(\textbf{x}^*(t)\) y el control \(u^*(t)\) que maximicen la altitud final \(h(t_f)\), cumpliendo con las ecuaciones de dinámica, las restricciones de frontera, y las limitaciones en el control. Esto se puede lograr mediante dos enfoques principales: el método directo y el método indirecto. En este trabajo, implementamos ambos enfoques utilizando el lenguaje de programación Julia y la paquetería OptimalControl.jl (s. f.), presentando las soluciones y resultados obtenidos.
5.2.1 Método Directo
El método directo consiste en discretizar las ecuaciones de estado y de control, convirtiendo el problema de un problema de Calculo Variacional en un problema de programación no lineal (PNL). Se utiliza IPOPT para resolver el NLP, con MUMPS como el solver lineal para hacer más eficiente el manejo de grandes sistemas dispersos.
Las paqueterías que se usarán son las siguientes
Definimos los datos del problema (las restricciones, condiciones de frontera y la dinámica)
Código
t0 = 0 # Tiempo Inicial
h0 = 1 # Altitud Inicial
v0 = 0 # Velocidad Inicial
m0 = 1 # Masa Inicial
vmax = 0.1 # Velocidad Máxima Autorizada
mf = 0.6 # Masa Final Objetivo
const Cd = 310 # Coeficiente de Arrastre
const Tmax = 3.5 # Empuje Máximo
const β = 500 # Tasa de caída de Densidad
const b = 2 # Constante positiva de acuerdo al motorDefinimos el problema de control optimo
Código
@def ocp begin
tf ∈ R, variable
t ∈ [ t0, tf ], time
x = (h, v, m) ∈ R³, state
u ∈ R, control
x(t0) == [ h0, v0, m0 ]
m(tf) == mf, (1)
0 ≤ u(t) ≤ 1
h(t) ≥ h0
0 ≤ v(t) ≤ vmax
ẋ(t) == F0(x(t)) + u(t) * F1(x(t))
h(tf) → max
end;
F0(x) = begin
h, v, m = x
D = Cd * v^2 * exp(-β*(h - 1)) # Drag
return [ v, -D/m - 1/h^2, 0 ]
end
F1(x) = begin
h, v, m = x
return [ 0, Tmax/m, -b*Tmax ]
endSe reuelve utilizando IPOPT (Métodos de Punto Interior)
Código
direct_sol = solve(ocp; grid_size=100)This is Ipopt version 3.14.17, running with linear solver MUMPS 5.7.3.
Number of nonzeros in equality constraint Jacobian...: 2104
Number of nonzeros in inequality constraint Jacobian.: 0
Number of nonzeros in Lagrangian Hessian.............: 1111
Total number of variables............................: 405
variables with only lower bounds: 101
variables with lower and upper bounds: 202
variables with only upper bounds: 0
Total number of equality constraints.................: 304
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 -1.0100000e+00 9.00e-01 2.00e+00 0.0 0.00e+00 - 0.00e+00 0.00e+00 0
1 -1.0090670e+00 8.99e-01 6.67e+01 1.3 1.67e+02 - 3.64e-03 5.93e-04f 1
2 -1.0000907e+00 8.74e-01 1.83e+02 1.0 6.64e+00 - 3.63e-02 2.83e-02h 1
3 -1.0023670e+00 8.37e-01 1.34e+04 1.0 6.91e+00 - 2.11e-01 4.19e-02f 1
4 -1.0025025e+00 7.70e-01 9.45e+03 1.5 4.04e+00 - 1.00e+00 8.09e-02f 1
5 -1.0033626e+00 7.16e-01 1.48e+05 2.3 3.56e+00 - 3.73e-01 6.94e-02f 1
6 -1.0142503e+00 9.62e-03 3.99e+04 2.3 7.16e-01 - 4.49e-01 9.90e-01h 1
7 -1.0101264e+00 4.21e-03 4.24e+05 1.8 5.32e-01 - 5.03e-01 9.90e-01h 1
8 -1.0068427e+00 3.20e-04 2.87e+06 0.9 2.44e-01 - 6.73e-01 9.91e-01h 1
9 -1.0067336e+00 2.64e-06 2.30e+07 0.1 7.39e-02 - 7.07e-01 1.00e+00f 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
10 -1.0067340e+00 1.13e-10 6.50e+05 -5.0 2.90e-04 - 9.89e-01 1.00e+00h 1
11 -1.0067350e+00 2.81e-10 7.20e+03 -7.0 4.26e-04 - 9.89e-01 1.00e+00h 1
12 -1.0078967e+00 9.07e-04 5.95e+03 -3.0 7.61e-01 - 6.55e-01 7.21e-01f 1
13 -1.0081866e+00 3.67e-06 9.38e+03 -9.0 1.31e-02 - 8.60e-01 1.00e+00h 1
14 -1.0091814e+00 1.57e-04 1.79e+02 -4.6 2.18e-01 - 1.00e+00 9.28e-01h 1
15 -1.0105115e+00 2.55e-04 7.17e+02 -4.4 3.09e-01 - 1.00e+00 5.75e-01h 1
16 -1.0114149e+00 2.34e-05 7.98e-04 -5.1 4.29e-02 - 1.00e+00 1.00e+00h 1
17 -1.0122891e+00 7.90e-05 1.24e+02 -5.8 1.25e-01 - 9.98e-01 7.92e-01h 1
18 -1.0125091e+00 2.83e-05 4.85e-04 -6.1 1.06e-01 - 1.00e+00 1.00e+00h 1
19 -1.0125585e+00 1.21e-05 8.67e-05 -6.8 1.01e-01 - 1.00e+00 1.00e+00h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
20 -1.0125675e+00 4.40e-06 3.01e-05 -7.3 9.80e-02 - 1.00e+00 1.00e+00h 1
21 -1.0125707e+00 1.75e-06 3.70e-02 -8.0 1.14e-01 - 1.00e+00 9.91e-01h 1
22 -1.0125714e+00 9.50e-07 4.98e-03 -8.6 1.37e-01 - 1.00e+00 9.94e-01h 1
Number of Iterations....: 22
(scaled) (unscaled)
Objective...............: -1.0125714428815755e+00 -1.0125714428815755e+00
Dual infeasibility......: 4.9826923912073525e-03 4.9826923912073525e-03
Constraint violation....: 3.5941849951814930e-07 9.4995894299454420e-07
Variable bound violation: 9.9621924876114321e-09 9.9621924876114321e-09
Complementarity.........: 8.1796266157923928e-09 8.1796266157923928e-09
Overall NLP error.......: 5.1484604087044977e-07 4.9826923912073525e-03
Number of objective function evaluations = 23
Number of objective gradient evaluations = 23
Number of equality constraint evaluations = 23
Number of inequality constraint evaluations = 0
Number of equality constraint Jacobian evaluations = 23
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations = 22
Total seconds in IPOPT = 3.820
EXIT: Optimal Solution Found.OptimalControlSolution5.2.2 Método Indirecto
El método indirecto se basa en la teoría del mínimo de Pontryagin y requiere resolver un sistema de ecuaciones diferenciales en términos de las variables de estado, de coestado y las condiciones de transversalidad. Para verificar que las condiciones de frontera y las restricciones del problema se cumplen, se define una función de disparo. Se resuelve un sistema no lineal para encontrar los valores óptimos de los parámetros usando métodos numéricos. Finalmente, se construye la trayectoria óptima mediante la concatenación de flujos de los diferentes regímenes de control.
Código
t = time_grid(direct_sol)
x = state(direct_sol)
u = control(direct_sol)
p = costate(direct_sol)
H1 = Lift(F1) # H1(x, p) = p' * F1(x)
φ(t) = H1(x(t), p(t)) # Función "switching"
g(x) = vmax - x[2] # Restricción del estado v ≤ vmaxCódigo
# Controles
u0 = 0 # Control Apagado
u1 = 1 # Control "Bang"
H0 = Lift(F0) # H0(x, p) = p' * F0(x)
H01 = @Lie { H0, H1 }
H001 = @Lie { H0, H01 }
H101 = @Lie { H1, H01 }
us(x, p) = -H001(x, p) / H101(x, p) # Control Singular
ub(x) = -(F0⋅g)(x) / (F1⋅g)(x) # Frontera del Control
μ(x, p) = H01(x, p) / (F1⋅g)(x) # Multiplicador asociado a
# la restricción de estado g
# Flujo
f0 = Flow(ocp, (x, p, tf) -> u0)
f1 = Flow(ocp, (x, p, tf) -> u1)
fs = Flow(ocp, (x, p, tf) -> us(x, p))
fb = Flow(ocp, (x, p, tf) -> ub(x), (x, u, tf) -> g(x), (x, p, tf) -> μ(x, p))
x0 = [ r0, v0, m0 ] # Estado Inicial
#Función de Disparo
function shoot!(s, p0, t1, t2, t3, tf)
x1, p1 = f1(t0, x0, p0, t1)
x2, p2 = fs(t1, x1, p1, t2)
x3, p3 = fb(t2, x2, p2, t3)
xf, pf = f0(t3, x3, p3, tf)
s[1] = constraint(ocp, :eq1)(x0, xf, tf) - mf # Restricción de Masa Final
s[2:3] = pf[1:2] - [ 1, 0 ] # Condiciones de Transversalidad
s[4] = H1(x1, p1) # H1 = H01 = 0
s[5] = H01(x1, p1) # Entrada del Arco Singular
s[6] = g(x2) # g = 0 al entrar en el Arco Limite
s[7] = H0(xf, pf) # Ya que el Tiempo Final es libre
end
# Suposición Inicial
η = 1e-3
t13 = t[ abs.(φ.(t)) .≤ η ]
t23 = t[ 0 .≤ (g ∘ x).(t) .≤ η ]
p0 = p(t0)
t1 = min(t13...)
t2 = min(t23...)
t3 = max(t23...)
tf = t[end]
println("p0 = ", p0)
println("t1 = ", t1)
println("t2 = ", t2)
println("t3 = ", t3)
println("tf = ", tf)
# Norma de la función de disparo en la solución
using LinearAlgebra: norm
s = similar(p0, 7)
shoot!(s, p0, t1, t2, t3, tf)
println("\nNorma de la función de disparo: ‖s‖ = ", norm(s), "\n")
ξ = [ p0 ; t1 ; t2 ; t3 ; tf ] # Suposición Inicial
# Función auxiliar con entradas agregadas
nle! = (s, ξ, λ) -> shoot!(s, ξ[1:3], ξ[4], ξ[5], ξ[6], ξ[7])
# Problema NLE con Suposición Inicial
prob = NonlinearProblem(nle!, ξ)
# Resolución de S(ξ) = 0
indirect_sol = solve(prob; abstol=1e-8, reltol=1e-8, show_trace=Val(true))
# Recuperamos la solución de coestado junto con los tiempos
p0 = indirect_sol.u[1:3]
t1 = indirect_sol.u[4]
t2 = indirect_sol.u[5]
t3 = indirect_sol.u[6]
tf = indirect_sol.u[7]
println("")
println("p0 = ", p0)
println("t1 = ", t1)
println("t2 = ", t2)
println("t3 = ", t3)
println("tf = ", tf)
# Norma de la función de disparo en la solución
s = similar(p0, 7)
shoot!(s, p0, t1, t2, t3, tf)
println("\nNorma de la función de disparo: ‖s‖ = ", norm(s), "\n")
f = f1 * (t1, fs) * (t2, fb) * (t3, f0) # concatenación de los flujos
flow_sol = f((t0, tf), x0, p0) # Calcular la solución:
# estado, coestado, control...A continuación se presenta una gráfica comparando el método directo (Sección 5.2.1) y el método indirecto (Sección 5.2.2)
5.2.3 Análisis de Sensibilidad
Se realiza un análisis de sensibilidad del problema variando el parámetro de la masa final en un intervalo \([0.45,0.65]\) tomando valores con distancia de \(0.05\). Este estudio permite determinar el límite máximo de altitud en función de diferente valores de masa final, evaluando como afectan las restricciones dinámicas y el consumo de combustible. A través de simulaciones númericas (haciendo uso del Sección 5.2.1), se puede examinar la relación entre la eficiencia del control y la altitud final alcanzada.
Código
using OptimalControl
using NLPModelsIpopt
using DataFramesCódigo
t0 = 0 # initial time
h0 = 1 # initial altitude
v0 = 0 # initial speed
m0 = 1 # initial mass
mass = [0.45, 0.5, 0.55, 0.6, 0.65] # final mass to target
Op_1 = DataFrame(rand(100, 5), [:t, :x_h, :x_v, :x_m, :u_t])
Op_2 = DataFrame(rand(100, 5), [:t, :x_h, :x_v, :x_m, :u_t])
Op_3 = DataFrame(rand(100, 5), [:t, :x_h, :x_v, :x_m, :u_t])
Op_4 = DataFrame(rand(100, 5), [:t, :x_h, :x_v, :x_m, :u_t])
Op_5 = DataFrame(rand(100, 5), [:t, :x_h, :x_v, :x_m, :u_t])
List_dataf = [Op_1, Op_2, Op_3, Op_4, Op_5]
const Cd = 310 # Drag Coeficient
const Tmax = 3.5 # Maximal Thrust
const β = 500 # density decay rate
const b = 2 # Positive real number depending on the engine
F0(x) = begin
h, v, m = x
D = Cd * v^2 * exp(-β*(h - 1)) # Drag force
return [ v, -D/m - 1/h^2, 0 ]
end
F1(x) = begin
h, v, m = x
return [ 0, Tmax/m, -b*Tmax ]
endCódigo
for i in [1,2,3,4,5]
@def ocp_as begin # definition of the optimal control problem
tf ∈ R, variable
t ∈ [ t0, tf ], time
x = (h, v, m) ∈ R³, state
u ∈ R, control
x(t0) == [ h0, v0, m0 ]
m(tf) == mass[i], (1)
0 ≤ u(t) ≤ 1
h(t) ≥ h0
ẋ(t) == F0(x(t)) + u(t) * F1(x(t))
h(tf) → max
end;
direct_sol_as = solve(ocp_as; grid_size=100)
t_as = direct_sol_as.time_grid
x_as = direct_sol_as.state
u_as = direct_sol_as.control
for j in 1:100
List_dataf[i].t[j] = t_as[j]
List_dataf[i].x_h[j] = x_as(t_as[j])[1]
List_dataf[i].x_v[j] = x_as(t_as[j])[2]
List_dataf[i].x_m[j] = x_as(t_as[j])[3]
List_dataf[i].u_t[j] = u_as(t_as[j])
end
end5.2.4 Interpretación de los resultados
La solución óptima obtenida revela una estructura de control caracterísitica en tres fases claramente diferenciadas. En la gráfica del empuje (control \(u(t)\)), observamos que el cohete tiene un inicio con empuje máximo durante el intervalo \([t_0,t_1]\), lo que corresponde a la fase de ascenso inicial donde se busca vencer la atracción gravitacional.
Posteriormente, en el siguiente intervalo \([t_1,t_2]\), el control muestra un comportamiento singular, donde \(u(t)\) toma valores intermedios entre 0 y 1. Esta fase corresponde a un arco singular, donde la función de conmutación \(\phi(t)\) se anula en un intervalo no trivial. Durante esta fase, el cohete mantiene un equilibrio entre el consumo de combustible y la ganancia de altitud.
En el intervalo \([t_2,t_3]\), observamos un arco límite donde la restricción de velocidad \(v(t)\leq v_{\max}\) se activa. El control \(u(t)\) ajusta automáticamente el empuje para mantener la velocidad en su límite superior sin sobrepasarla, evidenciado por el comportamiento plano de la curva de velocidad en este intervalo.
Finalmente, en \([t_3,t_f]\), el control se apaga por completo permitiendo que el cohete continúe su ascenso por inercia hasta alcanzar la altitud máxima, momento en el cual la velocidad se anula simultáneamente con el agotamiento del combustible hasta alcanzar la masa final \(m_f=0.6\).
La solución obtenida se caracteriza por la secuencia de arcos: Bang-Singular-Límite-Bang, denotada como (1-S-B-0) (Seywald y Cliff (1993)).
La comparación que vemos entre el método directo e indirecto (Figura 5.1) muestra una concordancia notable, validando ambos enfoques. Las pequeñas discrepancias observables en las transiciones entre fases se deben a la discretización empleada en el método directo (grid_size=100), mientras que el método indirecto captura estas transiciones con mayor precisión al resolver el sistema continuo de ecuaciones diferenciales.