Управление ракетой

Это руководство было изначально предоставлено Иэном Даннингом (Iain Dunning).

Цель данного руководства — продемонстрировать настройку и решение нелинейной задачи оптимизации.

В качестве примера используется задача оптимального управления нелинейной ракетой.

Для моделирования и решения задач оптимального управления также можно использовать расширение JuMP InfiniteOpt.jl.

В этом руководстве используются следующие пакеты:

using JuMP
import Ipopt
import Plots

Обзор

Наша цель — максимизировать конечную высоту полета вертикально запускаемой ракеты.

Мы можем контролировать тягу ракеты и должны учитывать ее массу, расход топлива, силу тяжести и аэродинамическое сопротивление.

Рассмотрим базовое описание модели (полное описание, включая параметры ракеты, см. в документе COPS3).

В нашей модели три переменные состояния:

скорость: ;
высота: ;
масса ракеты и оставшегося топлива: ;

и одна управляющая переменная:

тяга: .

Движение ракеты описывается тремя уравнениями:

Скорость подъема:
Ускорение:
Скорость уменьшения массы:

где сопротивление является функцией высоты и скорости, сила тяжести является функцией высоты, а является константой.

Эти силы определяются следующим образом:

and

Мы используем дискретную модель времени с фиксированным числом временных шагов .

Таким образом, наша цель — максимизировать .

Данные

Все параметры этой модели нормализованы до безразмерных и взяты из COPS3.

h_0 = 1                      # Начальная высота
v_0 = 0                      # Начальная скорость
m_0 = 1.0                    # Начальная масса
m_T = 0.6                    # Конечная масса
g_0 = 1                      # Сила тяжести на поверхности
h_c = 500                    # Используется для сопротивления
c = 0.5 * sqrt(g_0 * h_0)    # Тяга к массе топлива
D_c = 0.5 * 620 * m_0 / g_0  # Масштабирование сопротивления
u_t_max = 3.5 * g_0 * m_0    # Максимальная тяга
T_max = 0.2                  # Количество секунд
T = 1_000                    # Количество временных шагов
Δt = 0.2 / T;                # Время дискретного шага

Формулировка на языке JuMP

Сначала создадим модель и выберем оптимизатор. Так как это нелинейная программа, необходимо использовать нелинейный решатель, например Ipopt. Линейный решатель, например HiGHS, использовать нельзя.

model = Model(Ipopt.Optimizer)
set_silent(model)

Далее создадим переменные состояния и управляющую переменную, каждая из которых индексируется по t. В нелинейных программах желательно всегда предоставлять начальные решения для каждой переменной.

@variable(model, x_v[1:T] >= 0, start = v_0)           # Скорость
@variable(model, x_h[1:T] >= 0, start = h_0)           # Высота
@variable(model, x_m[1:T] >= m_T, start = m_0)         # Масса
@variable(model, 0 <= u_t[1:T] <= u_t_max, start = 0); # Тяга

Мы реализуем граничные условия, фиксировав значения переменных.

fix(x_v[1], v_0; force = true)
fix(x_h[1], h_0; force = true)
fix(x_m[1], m_0; force = true)
fix(u_t[T], 0.0; force = true)

Целью является максимизация высоты в конце полета.

@objective(model, Max, x_h[T])

x_h[1000]

Силы определяются как функции:

D(x_h, x_v) = D_c * x_v^2 * exp(-h_c * (x_h - h_0) / h_0)
g(x_h) = g_0 * (h_0 / x_h)^2

g (generic function with 1 method)

Уравнения движения реализованы как ограничения.

ddt(x::Vector, t::Int) = (x[t] - x[t-1]) / Δt
@constraint(model, [t in 2:T], ddt(x_h, t) == x_v[t-1])
@constraint(
    model,
    [t in 2:T],
    ddt(x_v, t) == (u_t[t-1] - D(x_h[t-1], x_v[t-1])) / x_m[t-1] - g(x_h[t-1]),
)
@constraint(model, [t in 2:T], ddt(x_m, t) == -u_t[t-1] / c);

Теперь оптимизируем модель и проверим, найдено ли решение:

optimize!(model)
@assert is_solved_and_feasible(model)
solution_summary(model)

* Solver : Ipopt

* Status
  Result count       : 1
  Termination status : LOCALLY_SOLVED
  Message from the solver:
  "Solve_Succeeded"

* Candidate solution (result #1)
  Primal status      : FEASIBLE_POINT
  Dual status        : FEASIBLE_POINT
  Objective value    : 1.01278e+00
  Dual objective value : 4.66547e+00

* Work counters
  Solve time (sec)   : 8.58810e-01
  Barrier iterations : 24

Наконец, построим график решения:

function plot_trajectory(y; kwargs...)
    return Plots.plot(
        (1:T) * Δt,
        value.(y);
        xlabel = "Time (s)",
        legend = false,
        kwargs...,
    )
end

Plots.plot(
    plot_trajectory(x_h; ylabel = "Altitude"),
    plot_trajectory(x_m; ylabel = "Mass"),
    plot_trajectory(x_v; ylabel = "Velocity"),
    plot_trajectory(u_t; ylabel = "Thrust");
    layout = (2, 2),
)

Дальнейшие действия

Поэкспериментируйте с различными значениями констант. Как меняется решение? В частности, что происходит при изменении T_max?
К членам в правой части уравнений движения применяется интегрирование по формуле прямоугольников. Измените уравнения так, чтобы вместо этого использовалось правило трапеций. (Например, x_v[t-1] станет 0.5 * (x_v[t-1] + x_v[t]).) Есть ли разница?