Código
using OptimalControl
using NLPModelsIpopt
using DataFrames6 Problema Alterno
6.1 Método Directo
Se formula un problema alterno en el que la función objetivo conciste en maximizar la masa final del cohete al alcanzar una altitud final fija. Este enfoque permite evaluar la eficiencia en el consumo de combustible. Para garantizar la coherencia en el análisis, se emplean los mismos datos y restricciones del problema original, preservando la dinámica del sistema y las condiciones operativas.
\[ \left\{ \begin{split} \text{max } &m(t_f)\\ \dot{h} &= v\\ \dot{v} &= \frac{T-D}{m}-\frac{1}{h^2}\\ \dot{m} &= -\frac{T}{c}\\ T &\in \left[0,T_{max}\right]\\ h(0) &= 1, v(0) = 0, m(0) = 1\\ h(t_f) &= 1.01257\\ t_f &\text{ libre} \end{split} \right. \]
Para obtener los resultados numéricos de este problema, resolvemos el problema mediante el método directo, implementado en el lenguaje de programación Julia.
Definimos los datos del problema
Código
t0_m = 0 # Tiempo Inicial
h0_m = 1 # Altitud Inicial
v0_m = 0 # Velocidad Inicial
m0_m = 1 # Masa Inicial
h_target = 1.01283 # Altitud Final Objetivo
mf_min = 0.2 # Masa Final Mínima Autorizada
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 motorCódigo
@def ocp_mass begin
tf_m ∈ R, variable
t ∈ [ t0_m, tf_m ], time
x = (h, v, m) ∈ R³, state
u ∈ R, control
x(t0_m) == [ h0_m, v0_m, m0_m ]
h(tf_m) == h_target, (1)
m(tf_m) ≥ mf_min
0 ≤ u(t) ≤ 1
h(t) ≥ h0_m
ẋ(t) == F0(x(t)) + u(t) * F1(x(t))
m(tf_m) → max
end;
F0(x) = begin
r, v, m = x
D = Cd * v^2 * exp(-β*(r - 1)) # Drag force
return [ v, -D/m - 1/r^2, 0 ]
end
F1(x) = begin
r, v, m = x
return [ 0, Tmax/m, -b*Tmax ]
endCódigo
direct_sol_mass = solve(ocp_mass; 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.: 1
Number of nonzeros in Lagrangian Hessian.............: 1111
Total number of variables............................: 405
variables with only lower bounds: 101
variables with lower and upper bounds: 101
variables with only upper bounds: 0
Total number of equality constraints.................: 304
Total number of inequality constraints...............: 1
inequality constraints with only lower bounds: 1
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.0000000e-01 9.00e-01 1.00e+00 0.0 0.00e+00 - 0.00e+00 0.00e+00 0
1 -5.2232164e-01 1.04e-02 5.64e+01 -0.5 9.00e-01 - 7.98e-01 9.90e-01f 1
2 -5.5285947e-01 5.46e-03 7.09e+02 -0.6 3.28e-01 2.0 9.98e-01 9.90e-01h 1
3 -6.4577962e-01 2.33e-03 6.46e+02 -1.1 2.23e-01 1.5 1.00e+00 1.00e+00h 1
4 -6.0124963e-01 1.62e-03 1.33e+07 -0.5 1.31e+00 1.0 1.00e+00 5.75e-01f 1
5 -5.4253619e-01 8.12e-04 6.06e+06 -0.6 5.78e-01 - 7.01e-01 4.81e-01h 1
6 -5.2059521e-01 5.98e-04 2.71e+06 -1.3 1.98e-01 - 5.63e-01 2.50e-01f 3
7 -4.9742555e-01 4.68e-04 7.92e+06 -2.0 1.17e-01 0.6 6.55e-01 2.50e-01h 3
8 -4.4368540e-01 1.85e-04 2.09e+01 -2.0 5.37e-02 1.0 1.00e+00 1.00e+00h 1
9 -4.9559548e-01 2.91e-04 3.46e+00 -2.3 1.90e-01 - 1.00e+00 1.00e+00f 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
10 -5.0093351e-01 2.44e-06 2.49e-02 -4.0 5.34e-03 0.5 1.00e+00 1.00e+00h 1
11 -5.0940839e-01 1.32e-05 1.79e-02 -5.5 1.63e-02 0.0 1.00e+00 1.00e+00f 1
12 -5.2708055e-01 5.86e-05 4.89e+01 -6.7 4.69e-02 -0.4 1.00e+00 8.30e-01h 1
13 -5.4994473e-01 1.12e-04 6.34e+00 -6.2 1.34e-01 -0.9 1.00e+00 4.89e-01h 1
14 -5.6065923e-01 1.10e-04 3.24e+01 -6.5 4.39e-01 -1.4 1.00e+00 8.76e-02h 1
15 -5.6809749e-01 1.04e-04 2.54e+01 -6.4 1.94e+00 -1.9 1.00e+00 2.15e-02f 1
16 -5.7634469e-01 8.77e-05 3.47e+01 -6.8 2.25e-01 -1.4 1.00e+00 2.36e-01h 1
17 -5.8121261e-01 8.55e-05 3.61e+01 -7.0 6.93e-01 -1.9 1.00e+00 1.03e-01h 1
18 -5.8419394e-01 8.01e-05 2.70e+01 -6.7 2.36e+00 -2.4 1.00e+00 3.62e-02h 1
19 -5.8792501e-01 7.02e-05 2.89e+01 -7.0 1.44e-01 -2.0 1.00e+00 2.01e-01h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
20 -5.8983597e-01 6.28e-05 2.71e+01 -7.1 4.17e-01 -2.4 1.00e+00 1.05e-01h 1
21 -5.9070972e-01 6.00e-05 2.63e+01 -7.1 1.21e+00 -2.9 1.00e+00 4.51e-02h 1
22 -5.9127966e-01 5.82e-05 1.08e+00 -6.5 2.56e+00 -3.4 1.00e+00 2.99e-02h 1
23 -5.9161565e-01 5.71e-05 2.40e+02 -5.6 8.45e+00 - 1.00e+00 1.94e-02h 1
24 -6.0422424e-01 1.91e-03 1.65e+01 -5.6 7.61e-01 - 1.00e+00 9.16e-01h 1
25 -6.0177706e-01 8.27e-04 5.46e+00 -5.6 6.45e-01 - 9.29e-01 5.67e-01h 1
26 -5.9948065e-01 6.25e-05 1.74e-02 -5.1 7.21e-01 - 1.00e+00 1.00e+00h 1
27 -5.9983571e-01 1.38e-05 3.94e-04 -5.6 2.93e-01 - 1.00e+00 1.00e+00h 1
28 -5.9999480e-01 5.28e-06 1.69e-04 -6.5 1.45e-01 - 1.00e+00 1.00e+00h 1
29 -6.0001951e-01 1.32e-06 4.56e-05 -7.6 1.01e-01 - 1.00e+00 1.00e+00h 1
Number of Iterations....: 29
(scaled) (unscaled)
Objective...............: -6.0001951350455163e-01 -6.0001951350455163e-01
Dual infeasibility......: 4.5634917827996978e-05 4.5634917827996978e-05
Constraint violation....: 5.0046210509080424e-09 1.3167889997001136e-06
Variable bound violation: 0.0000000000000000e+00 0.0000000000000000e+00
Complementarity.........: 2.6395316602666028e-07 2.6395316602666028e-07
Overall NLP error.......: 2.6395316602666028e-07 4.5634917827996978e-05
Number of objective function evaluations = 36
Number of objective gradient evaluations = 30
Number of equality constraint evaluations = 36
Number of inequality constraint evaluations = 36
Number of equality constraint Jacobian evaluations = 30
Number of inequality constraint Jacobian evaluations = 30
Number of Lagrangian Hessian evaluations = 29
Total seconds in IPOPT = 3.774
EXIT: Optimal Solution Found.OptimalControlSolution6.2 Análisis de Sensibilidad
Al igual que con el Capítulo 5, se realiza un análisis de sensibilidad del Capítulo 6 variando el parámetro de la alura final tomando los datos de la salida del Sección 5.2.3 del problema original. Este análisis se lleva a cabo 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 masa final alcanzada.
Código
using OptimalControl
using NLPModelsIpopt
using DataFrames
t0_ms = 0 # Tiempo Inicial
h0_ms = 1 # Altitud Inicial
v0_ms = 0 # Velocidad Inicial
m0_ms = 1 # Masa Inicial
h_targets = [1.038372, 1.026699, 1.018521, 1.012829] # Altitud Final Objetivo
mf_min = 0.2 # Masa Final Mínima Autorizada
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 motor
Op_1_ms = DataFrame(rand(100, 5), [:t, :x_h, :x_v, :x_m, :u_t])
Op_2_ms = DataFrame(rand(100, 5), [:t, :x_h, :x_v, :x_m, :u_t])
Op_3_ms = DataFrame(rand(100, 5), [:t, :x_h, :x_v, :x_m, :u_t])
Op_4_ms = DataFrame(rand(100, 5), [:t, :x_h, :x_v, :x_m, :u_t])
#Op_5_ms = DataFrame(rand(100, 5), [:t, :x_h, :x_v, :x_m, :u_t])
List_dataf_ms = [Op_1_ms, Op_2_ms, Op_3_ms, Op_4_ms]
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]
@def ocp_ms begin # definition of the optimal control problem
tf_ms ∈ R, variable
t ∈ [ t0_ms, tf_ms ], time
x = (h, v, m) ∈ R³, state
u ∈ R, control
x(t0_ms) == [ h0_ms, v0_ms, m0_ms ]
h(tf_ms) == h_targets[i], (1)
m(tf_ms) ≥ mf_min
0 ≤ u(t) ≤ 1
h(t) ≥ h0_ms
ẋ(t) == F0(x(t)) + u(t) * F1(x(t))
m(tf_ms) → max
end;
sol_ms = solve(ocp_ms; grid_size=100)
t_ms = sol_ms.time_grid
x_ms = sol_ms.state
u_ms = sol_ms.control
for j in 1:100
List_dataf_ms[i].t[j] = t_ms[j]
List_dataf_ms[i].x_h[j] = x_ms(t_ms[j])[1]
List_dataf_ms[i].x_v[j] = x_ms(t_ms[j])[2]
List_dataf_ms[i].x_m[j] = x_ms(t_ms[j])[3]
List_dataf_ms[i].u_t[j] = u_ms(t_ms[j])
end
end