F-16 Fighting Falcon — нелинейная продольная динамика¶

General Dynamics F-16 Fighting Falcon — американский многофункциональный лёгкий истребитель 4-го поколения. Данный модуль предоставляет нелинейную модель продольного канала полёта, реализованную на чистом Python/NumPy. Аэродинамические коэффициенты интерполируются кубическими сплайнами из таблиц аэродинамической трубы, обеспечивая высокоточную динамику в широком диапазоне углов атаки и отклонений управляющих поверхностей. В комплекте идёт среда Gymnasium для обучения агентов управления.

Быстрый старт

Запустите среду или модель в пару строк кода.

К примеру
API модели

Документация Python-класса нелинейной продольной динамики F-16.

К API
Среда Gymnasium

Готовая среда для обучения агентов управления.

К среде
Теория

Нелинейные уравнения состояния и структура аэродинамических таблиц.

К модели

ЛТХ (справочно)¶

Параметр	Значение
Модификация	F-16A Block 10
Размах крыла, м	9.45
Длина самолёта (со штангой ПВД), м	15.03
Высота самолёта, м	5.09
Площадь крыла, м²	27.87
Угол стреловидности, °	40.0
Нормальная взлётная масса, кг	11467

Как устроен объект управления¶

В отличие от линейной модели F-16, использующей постоянные матрицы коэффициентов \(A\) и \(B\), нелинейная модель вычисляет аэродинамические силы и моменты из табличных данных на каждом шаге. Уравнения движения интегрируются численно (методом Эйлера или RK4).

Вектор состояния и управляющее воздействие:

\[ x = \begin{bmatrix} \alpha & \omega_z & \delta_{\text{stab}} & \dot{\delta}_{\text{stab}} \end{bmatrix}^{\top}, \quad u = \delta_{\text{stab,act}} \]

ПеременныеПараметры по умолчанию

\(\alpha\): угол атаки, рад
\(\omega_z\): угловая скорость тангажа, рад/с
\(\delta_{\text{stab}}\): отклонение стабилизатора, рад
\(\dot{\delta}_{\text{stab}}\): скорость отклонения стабилизатора, рад/с
\(\delta_{\text{stab,act}}\): команда на стабилизатор (управляющее воздействие), рад

Параметр	Обозначение	Значение
Масса самолёта	\(m\)	9295.44 кг
Площадь крыла	\(S\)	27.87 м²
Средняя аэродинамическая хорда	\(b_A\)	3.45 м
Момент инерции тангажа	\(J_z\)	75673.6 кг м²
Постоянная времени стабилизатора	\(T_{\text{stab}}\)	0.03 с
Коэффициент демпфирования стабилизатора	\(\xi_{\text{stab}}\)	0.707
Высота полёта	\(H\)	3000 м
Скорость полёта	\(V\)	150 м/с
Ускорение свободного падения	\(g\)	9.80665 м/с²

О единицах измерения

Внутри модели все углы и угловые скорости задаются в радианах. Среда Gymnasium принимает действия в градусах (диапазон \(\pm 25°\)) для совместимости с существующими агентами и преобразует их в радианы внутренне.

Математическая модель¶

Модель решает следующую систему нелинейных ОДУ на каждом шаге:

\[ \dot{\alpha} = \omega_z - \frac{R_y - mg}{mV} \]

\[ \dot{\omega}_z = \frac{M_{Rz}}{J_z} \]

где аэродинамическая сила и момент:

\[ R_y = qS \cdot C_y(\alpha, \beta, \delta_{\text{stab}}, \delta_{\text{lef}}, \omega_z, V, b_A, \delta_{\text{sb}}) \]

\[ M_{Rz} = qS b_A \cdot m_z(\alpha, \beta, \delta_{\text{stab}}, \delta_{\text{lef}}, \omega_z, V, b_A, \delta_{\text{sb}}) + x_{\text{cg}} \cdot R_y \]

Скоростной напор \(q\) вычисляется по модели стандартной атмосферы (МСА):

\[ q = \frac{1}{2} \rho(H) V^2, \qquad \rho = \rho_0 \left(\frac{T_0 - LH}{T_0}\right)^{\frac{g}{LR} - 1} \]

Аэродинамические таблицы¶

Коэффициенты \(C_y\) и \(m_z\) не являются постоянными — это многомерные функции, интерполируемые из данных аэродинамической трубы, хранящихся в файлах .npz. Таблицы интерполируются кубическими сплайнами (csaps).

Модель привода¶

Стабилизатор моделируется как колебательное звено второго порядка с ограничениями по положению и скорости:

\[ \ddot{\delta}_{\text{stab}} = \frac{-2 T_{\text{stab}} \xi_{\text{stab}} \dot{\delta}_{\text{stab}} - \delta_{\text{stab}} + \delta_{\text{stab,act}}}{T_{\text{stab}}^2} \]

Ограничения привода

По умолчанию в модели применяются предельные значения управления:

Максимальное отклонение стабилизатора: \(\pm 25°\)
Максимальная скорость отклонения: \(\pm 60°/\text{s}\)

Методы интегрирования¶

Доступны два численных интегратора:

Euler (по умолчанию) — метод Эйлера первого порядка.
RK4 — метод Рунге-Кутты 4-го порядка, обеспечивает более высокую точность при том же шаге.

Источник данных¶

Stevens & Lewis, "Aircraft Control and Simulation".
Аэродинамические таблицы F-16A по данным аэродинамической трубы, сохранённые в формате NumPy .npz.

Быстрый старт¶

GymnasiumТолько модель

import gymnasium as gym
import numpy as np

from tensoraerospace.envs import NonlinearLongitudinalF16
from tensoraerospace.utils import generate_time_period
from tensoraerospace.signals.standard import sinusoid

dt = 0.01
tp = generate_time_period(tn=20, dt=dt)
number_time_steps = len(tp)
reference_signal = sinusoid(
    degree=3, tp=tp, frequency=0.1, output_rad=True
).reshape(1, -1)

env = gym.make(
    'NonlinearLongitudinalF16-v0',
    number_time_steps=number_time_steps,
    initial_state=np.array([0.0, 0.0]),
    reference_signal=reference_signal,
    dt=dt,
    integrator="euler",
)

state, info = env.reset()
for _ in range(number_time_steps - 1):
    action = np.array([0.0])  # градусы
    state, reward, terminated, truncated, info = env.step(action)
    if terminated or truncated:
        break

import numpy as np
from tensoraerospace.aerospacemodel.f16.nonlinear.longitudinal import LongitudinalF16

dt = 0.01
number_time_steps = 200

# Состояние: [alpha, wz, stab, dstab] (рад)
x0 = np.array([0.0, 0.0, 0.0, 0.0])

model = LongitudinalF16(
    x0=x0,
    selected_state_output=["alpha", "wz"],
    dt=dt,
    integrator="rk4",
)

for t in range(number_time_steps - 1):
    u = np.array([np.radians(-2.0)])  # команда на стабилизатор (рад)
    state_next = model.run_step(u)

Python API¶

МодельСреда GymnasiumПараметры

`LongitudinalF16(x0, selected_state_output=None, t0=0, dt=0.01, integrator='euler')` ¶

Bases: ModelBase

F-16 in isolated longitudinal channel.

Action: stab_act (elevator command, rad). State: alpha, wz, stab, dstab (rad / rad·s⁻¹).

`current_state` `property` ¶

Most recent state as a flat 1-D ndarray (alpha, wz, stab, dstab).

`NonlinearLongitudinalF16(initial_state, reference_signal, number_time_steps, tracking_states=None, state_space=None, control_space=None, output_space=None, reward_func=None, use_reward=True, dt=0.01, integrator='euler', control_bias=0.0, feedforward_fn=None, airspeed=200.0, render_mode=None, chart_states=('alpha', 'wz', 'stab'), trail_length=None, initial_pitch=0.0, damage_profile=None, damage_observable=False, damage_event_callback=None)` ¶

Bases: Env

Gymnasium env over the pure-numpy nonlinear F-16 longitudinal model.

Parameters:

Name	Type	Description	Default
`initial_state`	`ndarray`	Initial state vector. Either the full 4-element model state `[alpha, wz, stab, dstab]` or a shorter vector matching `state_space` (missing components are zero-filled). Values are interpreted in radians (matching the underlying model).	required
`reference_signal`	`ndarray`	Reference trajectory of shape `(1, T)`, radians.	required
`number_time_steps`	`int`	Episode length.	required
`tracking_states`	`list[str] \| None`	Names of tracked states. Defaults to `["alpha"]`.	`None`
`state_space`	`list[str] \| None`	Subset of model states to expose as the observation. Observations are returned in radians. Defaults to `["alpha", "wz"]`.	`None`
`control_space`	`list[str] \| None`	Names of control channels. Defaults to `["stab"]`.	`None`
`output_space`	`list[str] \| None`	Compatibility alias of `state_space`.	`None`
`reward_func`	`Callable[[ndarray, ndarray, int], ndarray \| float] \| None`	Custom reward callable `(state, ref_signal, ts) -> float`.	`None`
`use_reward`	`bool`	If False, reward is fixed at 1.0 each step.	`True`
`dt`	`float`	Discretisation step (s). Defaults to 0.01.	`0.01`
`integrator`	`Literal['euler', 'rk4']`	`"euler"` (default, matches matlab) or `"rk4"`.	`'euler'`
`control_bias`	`float`	Constant offset (degrees) added to every action before clipping and conversion to radians. Use this to operate the agent in "delta around trim" space when the trim control is non-zero (e.g., on the nonlinear F-16, trim elevator is ~-4.45°). Default 0.	`0.0`
`feedforward_fn`	`Callable[[int, ndarray], float] \| None`	Optional callable `(time_step, reference_signal) ->ff_action_deg`. Whenever set, the env adds the returned offset (degrees) to the agent's action before clipping. Use this to inject a precomputed feedforward map (e.g., trim elevator as a function of reference angle of attack) so a reactive agent like IHDP only has to learn the small disturbance correction instead of the full reference-tracking trajectory. `time_step` is the current step index; `reference_signal` is the same array given at construction time. The callable can return a scalar or an array of length `len(control_space)`.	`None`

Action units: by convention IHDP and most existing tensoraerospace agents were tuned with the linear F-16 env, whose elevator action is interpreted in degrees with magnitude limit 25 and rate limit 60. This env keeps the same convention so those agent settings transfer: action is in degrees, range [-25, 25], and is converted to radians internally before being handed to the underlying numpy model.

`default_reward(state, ref_signal, ts)` `staticmethod` ¶

Tracking-error reward with a small angular-rate penalty.

`F16LongParameters(m=9295.44, S=27.87, bA=3.45, Jz=75673.6, Tstab=0.03, Xistab=0.707, maxabsstab=(lambda: math.radians(25))(), maxabsdstab=(lambda: math.radians(60))(), lef=0.0, sb=0.0, g=9.80665, Oy=3000.0, V=150.0, damage_state=None, damage_geometry=None)` `dataclass` ¶