Экспериментируем с кодом модели самолета из книги [1]. В книге рассмотрен отклик трехстепенной модели самолета на импульс руля высоты. Здесь мы заменим код из книги на встроенную функцию Matlab, которая называется ode45. Построим тот же график зависимости угла атаки от времени. Также, я добавил возможность задавать возмущение балансированного полета по точкам.
Отклик трехстепенной модели самолета: набор программ
Рассмотрим по порядку, как работает тот код, что приведен в книге [1]. В начале замечу, что для работы модели нужны те программы, которые рассматривались ранее в статье о балансировке. Теперь, рассмотрим новые скрипты и функции, добавленные для моделирования процесса во времени. Итак есть следующие программы:
nlsim — главный скрипт, моделирующий отклик трехстепенной модели самолета. Выполнение начинается с его запуска. Этот скрипт задает начальные векторы состояния и управления по ранее рассчитанному сбалансированному состоянию самолета. Далее, на каждом шаге по времени мы вычисляем новый вектор состояния при помощи функции RK4. В начале каждого шага по времени, текущие переменные состояния x сохраняются в выходной вектор Y. Вектор Y мы в конце работы скрипта выводим на график. Импульс руля высоты, тоже задаваемый скриптом nlsim, представляет собой установку определенных значений во вторую ячейку вектора управления u(2). Это делается на нужных шагах по времени (между 1-й и 2-й секундами) при помощи набора условных операторов if-else. Установку времени моделирования (50 секунд) содержит переменная runtime.
RK4 — это функция. Функции отличаются от скрипта тем, что имеют входные аргументы и имя функции указанное в файле должно совпадать с именем файла без расширения. Данная функция производит вычисление одного шага метода Рунге-Кутты [2, стр.226]. Каждый вызов RK4 четыре раза определяет производные состояния, это как раз делается путем вызова модели transp() функцией feval(). И, в завершении, определяется значение вектора состояния xnew для следующего временного шага
Эксперимент
Эксперимент заключается в том, что я поменяю функцию RK4 на матлабовскую ode45. Это реализация метода Дормана-Принса, являющегося разновидностью метода Рунге-Кутты. Для ode45, Matlab имеет систему настроек, которые в простенькой реализации RK4 конечно отсутствуют. Некоторые из этих настроек введены в новый вариант скриптов, моделирующих отклик трехстепенной модели самолета. Новые скрипты имеют суффикс «ode45».
nlsimODE45 — функция этого файла та же, что у nlsim. Итерации по всем временным шагам здесь выполняются в функции ode45. Поэтому цикла по временным шагам нет. Однако, нам нужно устанавливать управление. И начальное и возмущающее. Кроме того, ode45 в качестве параметра odefun принимает функцию со входным аргументом x и выходным dx (производных). Размерность должна совпадать.
transpOde45 — целевая функция для применения в качестве odefun в ode45. Работа с вектором управления инкапсулирована полностью в эту функцию. Импульс руля высоты задается здесь. Это сделано из-за требований ode45 к аргументу odefun. Применить в ode45 напрямую transp(time, x, u) не получится. Слишком много аргументов.
Кроме того, я выделил возмущения вводимые в вектор управления в текстовые файлы. ControlIncrementsElev01 — для возмущения в установке руля, ControlIncrementsThtl01 — для возмущения по тяге. Нужное возмущение выбирается при чтении из файла в ctrlData (см. nlsimODE45 ). Теперь можно задавать по точкам величину отклонения от сбалансированного состояния и рассчитывать отклик. При вычислениях, когда программа ode45 вызывает transpOde45 происходит вызов getU_3Dof. Эта функция интерполирует временные шаги, попавшие между точками.
Работа
В архиве к этой статье приложена вся рабочая папка. Установки в скриптах при простом запуске команды nlsimODE45, обеспечат моделирование отклика на знакопеременное отклонение на 2 градуса руля высоты при балансированном полете со скоростью 250 футов в секунду. Далее, рассмотрим, как менять модель и условия задачи.
Характеристики модели менять нужно в файле transp, рассмотренном в предыдущей статье. Подробнее рассмотрим, как поменять скорость полета или тип возмущения.
Сперва, надо trimTable скорректировать, установив нужную скорость (пусть 250 футов в секунду):
x(1)=250; % TAS (Vt), ft/s.
Далее, запускаем этот скрипт и получаем тягу, отклонение руля высоты и угол атаки в градусах для состояния балансировки.
trimTable TAS(ft/s) = 1.844959e-01, H(ft) = 0 Thrl(0...1) Elev(deg) AoA(deg) 0.184496 -9.218456 9.277302
TAS(ft/s) = 1.844959e-01, H(ft) = 0 Thrl(0…1) Elev(deg) AoA(deg) 0.184496 -9.218456 9.277302[/code]
Мы решили задачу предыдущей статьи о балансировке для заданной скорости. Полученный угол атаки входит в вектор состояния, а тяга и отклонение руля входят в вектор управления. Именно к этим значениям управляющих величин (тяга и/или угол отклонения руля) программа nlsimODE45 будет прибавлять инкременты для введения возмущения в моделируемый процесс. Полученное состояние балансировки мы прописываем в скрипт nlsimODE45. Точнее, учитываем его в глобальных переменных x0 и u_trim. Также прописываем в матлабовскую «читалку» dlmread имя файла с возмущениями по управлению.
ctrlData = dlmread('ControlIncrementsElev01.txt',' ',1,0);
% Вектор состояния балансировки 'x0' и управление балансировки 'u_trim'.
global x0;
x0=[250 0.16192 0.16192 0 0]; % Стр.195 [TAS AoA Pitch Pitch_Rate Alt]
global u_trim;
xcg=0.25; land=0;
u_trim=[0.1845 -9.2184 xcg land]; % [thtl elev xcg land]
Обратите внимание, 0.16192 — это угол атаки 9.277302, только не в градусах, а в радианах. Он равен тангажу.
Далее — выполняем вызов скрипта для расчета отклика на возмущения относительно состояния балансировки, заданные в указанном файле:
nlsimODE45 506 successful steps 6 failed attempts 3073 function evaluations
Полученный график зависимости угла атаки от времени хорошо совпадает с приведенным в книге[1,стр.196]. На графике угла атаки видно короткопериодическую и длиннопериодическую составляющие продольного движения самолета после отклонения руля.
Напоминаю, что теперь, при помощи текстовых файлов, в которых задается возмущение вектора управления, можно задавать в том числе возмущения протяженные по времени на несколько временных шагов. Также можно использовать не линейную, а например кубическую интерполяцию между точками управляющего вектора. Это настраивается в функции getU_3Dof(time). Точнее, в вызове u_out = interp1(t,u’,time,’linear’, ‘extrap’);
Скачать
Литература
[1] — B.Stevens, F.Lewis, E.Johnson — Aircraft Control and Simulation (3ed).
[2] — В.М. Вержбицкий — Численные методы. Математический анализ и обыкновенные дифференциальные уравнения. М. «Высшая школа» 2001 г.