Разработка математической модели системы управления углом тангажа летательного аппарата

Система управления летательным аппаратом (ЛА) включает в себя несколько контуров управления для управления скоростью, высотой и навигацией самолёта.

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

В данном практическом руководстве:

• Представлен процесс разработки математической модели продольного движения самолёта[1] с помощью блока Jython;

• Приведён процесс создания модели управления углом тангажа ЛА с использованием библиотеки блоков Автоматика;

• Представлена Python программа для автоматизированной настройки коэффициентов ПИД-регулятора системы управления углом тангажа.

1.  Математическая модель продольного движения самолёта

В данном практическом руководстве будет выполнено моделирование продольного движения гражданского самолёта, схематично представленного на рисунке ниже[1].

Система координат гражданского самолёта

На рисунке выше приняты следующие обозначения.

Обозначения	Описание
Чёрные оси координат	Оси координат: \(\mathbf{x, z}\) - оси связанной СК; \(\mathbf{x', z'}\) - оси земной СК.
Красные векторы сил	Основные силы, действующие на ЛА при продольном движении: Thrust - сила тяги, Lift - подъёмная сила, Drag - сила сопротивления, Weight - вес ЛА.

В модели гражданского самолёта используются следующие допущения:

1. ЛА находится в устойчивом состоянии на постоянной высоте и скорости;

2. Изменение угла тангажа не приводит к изменению скорости самолета.

Принятые допущения позволяют сформировать модель продольного движения ЛА в виде системы уравнений (1).

(1)\[\begin{split}\left\{ \begin{aligned} \dot{\alpha} &= \mathrm{\mu} \cdot \mathrm{\Omega} \cdot \Bigg[ -\left(\mathrm{C}_{\mathrm{L}} + \mathrm{C}_{\mathrm{D}}\right) \cdot \mathrm{\alpha} + \frac{q}{\left(\mathrm{\mu} - \mathrm{C}_{\mathrm{L}}\right)} - \left(\mathrm{C}_{\mathrm{W}} \cdot \sin\mathrm{\gamma}\right) \cdot \mathrm{\theta} + \mathrm{C}_{\mathrm{L}} \Bigg]; \\ \dot{q} &= \frac{\mathrm{\mu} \cdot \mathrm{\Omega}}{2 \cdot i_{yy}} \Bigg[ \left(\mathrm{C}_{\mathrm{M}} - \mathrm{\eta} \cdot \left(\mathrm{C}_{\mathrm{L}} + \mathrm{C}_{\mathrm{D}}\right)\right) \cdot \alpha + \left(\mathrm{C}_{\mathrm{M}} + \mathrm{\sigma} \cdot \mathrm{C}_{\mathrm{M}} \cdot \left(1 - \mathrm{\mu} \cdot \mathrm{C}_{\mathrm{L}}\right)\right) \cdot q \\ &\quad + \left(\mathrm{\eta} \cdot \mathrm{C}_{\mathrm{W}} \cdot \sin\mathrm{\gamma}\right) \cdot \delta \Bigg]; \\ \dot{\theta} &= \Omega \cdot q. \end{aligned} \right.\end{split}\]

В уравнении (1) используются следующие множители:

(2)\[\begin{split}\begin{aligned} \mu &= \frac{\rho \cdot S \cdot \dot{c} }{4 \cdot m}; \\ \Omega &=\frac{2 \cdot U}{\dot{c}}; \\ \sigma &=\frac{1}{1+\mu \cdot C_{L}}; \\ \eta &=\mu \cdot \sigma \cdot C_{M}. \end{aligned}\end{split}\]

В уравнениях (1) и (2) приняты следующие обозначения:

Обозначение	Описание
\(\alpha\)	Угол атаки, приведён на рисунке выше
\(\dot{\alpha}\)	Производная первого порядка угла атаки
\(\theta\)	Угол тангажа (см. рисунок выше)
\(\dot{\theta}\)	Производная первого порядка угла тангажа
\(q\)	Угловая скорость тангажа
\(\dot{q}\)	Производная первого порядка угловой скорости тангажа
\(\delta\)	Угол руля высоты (см. рисунок выше)
\(\rho\)	Плотность воздуха
\(S\)	Площадь крыла
\(\dot{c}\)	Средняя длина хорды
\(m\)	Масса ЛА
\(U\)	Скорость установившегося полета
\(C_{T}\)	Коэффициент тяги
\(C_{D}\)	Коэффициент сопротивления
\(C_{L}\)	Коэффициент подъёмной силы
\(C_{W}\)	Коэффициент веса
\(C_{M}\)	Коэффициент момента тангажа
\(\gamma\)	Угол направления полета (см. рисунок выше)
\(i_{yy}\)	Нормализованный момент инерции

В случае использования коэффициентов гражданского самолёта уравнение (1) принимает следующий вид (3).

(3)\[\begin{split}\left\{ \begin{aligned} \dot{\alpha} &= -0.313 \cdot \alpha + 56.7 \cdot q + 0.232 \cdot \delta; \\ \dot{q} &= -0.0139 \cdot \alpha -0.426 \cdot q + 0.0203 \cdot \delta; \\ \dot{\theta} &= 56.7 \cdot q. \end{aligned} \right.\end{split}\]

Для последующей реализации модели продольного движения ЛА в блоке Jython потребуется матричное представление уравнения (3).

(4)\[\begin{split}\left[\begin{array}{c}\dot{\alpha} \\ \dot{q} \\ \dot{\theta} \end{array}\right] = \left[\begin{array}{ccc} -0.313 & 56.7 & 0\\ -0.0139 & -0.426 & 0 \\ 0 & 56.7 & 0 \end{array}\right] \cdot \left[\begin{array}{c} \alpha \\ q \\ \theta \end{array}\right]+\left[\begin{array}{c} 0.232 \\ 0.0203 \\ 0 \end{array}\right] \cdot \delta \end{split}\]

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

(5)\[\begin{split}y=\left[\begin{array}{ccc}0& 0& 1\end{array}\right] \cdot \left[\begin{array}{c}\alpha \\ q\\ \theta \end{array}\right] \end{split}\]

Для решения задачи по определению угла тангажа \(\theta\) выполняется интегрирование \(\dot{\theta}\) производной угла тангажа численным методом Эйлера в соответствии с уравнением (6).

(6)\[\begin{split}\left[\begin{array}{c}\alpha _{t+\mathrm{\Delta }t}\\ q_{t+\mathrm{\Delta }t}\\ \theta _{t+\mathrm{\Delta }t}\end{array}\right]=\left[\begin{array}{c}\alpha _{t}\\ q_{t}\\ \theta _{t}\end{array}\right]+\mathrm{\Delta }t\cdot \left[\begin{array}{c}\dot{\alpha }\\ \dot{q}\\ \dot{\theta }\end{array}\right]\end{split}\]

В уравнении (6) приняты следующие обозначения:

Обозначение	Описание
\(\Delta t\)	Шаг интегрирования. Используемое значение шага для настройки проекта приведено в таблице ниже.
\(\alpha_{t}\)	Значение угла атаки на предыдущем шаге
\(q_{t}\)	Значение угловой скорости тангажа на предыдущем шаге
\(\theta_{t}\)	Значение угла тангажа на предыдущем шаге
\(\alpha_{t + \Delta t}\)	Значение угла атаки на текущем шаге
\(q_{t + \Delta t}\)	Значение угловой скорости тангажа на текущем шаге
\(\theta_{t + \Delta t}\)	Значение угла тангажа на текущем шаге

2.  Обзор демонстрационной модели управления углом тангажа самолёта в приложении REPEAT VISION

В наборе демонстрационных проектов приложения REPEAT VISION доступен пример модели управления углом тангажа самолёта. Далее выполним анализ данного примера, который нам поможет описать все необходимые шаги по реализации модели управления углом тангажа ЛА на платформе REPEAT VISION.

В приложении REPEAT VISION добавьте себе проект [DEMO] Система управления углом тангажа самолёта из папки Демо проекты и откройте его для последующего выполнения обзора.

Демонстрационный проект системы управления углом тангажа самолёта

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

Настройка	Значение
Синхронизация с реальным временем	Выключена
Время моделирования	40 \([с]\)
Шаг интегрирования	1 \([мс]\)
Шаг дискретизации	100 \([мс]\)

2.1.  Реализация модели продольного движения ЛА

Для расчёта угла тангажа \(\theta\) в модели продольного движения ЛА используется два блока:

Блок	№ блока в проекте	Назначение для модели
Константа	Блок №1	Необходим для назначения \(\theta\), уставки угла тангажа (см. схему ЛА), величиной \(0.2 [рад]\)
Jython	Блок №8	Необходим для непосредственной реализации уравнения (3) с использованием матричного вида (4)

Модель продольного движения ЛА реализована на языке программирования Python с применением блока Jython.

Для обзора модели продольного движения ЛА требуется в свойствах блока №8 Jython нажать на кнопку Редактировать код блока.

В блоке Jython доступны три вкладки для реализации кода Python программы.

Инициализация

Код данного этапа выполняется 1 раз до начала моделирования с целью подготовки к процессу расчёта. На данном этапе необходимо выполнить:

1. Импорт необходимых библиотек;

2. Создать объекты программы, у которых требуется сохранять состояние между выполнением шагов интегрирования в процессе расчёта;

3. Создать необходимые функции и классы для последующего вызова на этапе моделирования, в Коде расчета.

Код расчета

Код данного этапа выполняется на каждом шаге интегрирования от начала расчёта до его остановки. На данном этапе необходимо:

1. Вызывать функции и методы классов Python программы, отвечающих за моделирование;

2. Использовать значения in входных портов блока Jython и возвращать результаты на out выходные порты блока Jython;

Объекты, созданные в рамках кода расчёта не сохраняют своего состояния после завершения шага интегрирования.

Если вашим объектам необходимо сохранять состояние между шагами интегрирования, то их необходимо создавать и настраивать на этапе инициализации, а затем использовать в коде расчёта.

Завершение

Код данного этапа выполняется 1 раз при остановке расчёта.

При завершении расчёта может выполняться обработка данных объектов Python программы, у которых требовалось сохранять состояние между выполнением шагов интегрирования в процессе расчёта.

Далее в текстовом редакторе кода блока Jython требуется перейти на вкладку Инициализация для анализа кода подготовки и настройки модели продольного движения ЛА.

В рамках инициализации создаётся матрица и вектор коэффициентов, а также реализуется функция для решения системы уравнений (4) с целью дальнейшего получения вектора угловых значений атаки и тангажа при продольном движении ЛА на следующем этапе выполнения моделирования.

Инициализация

from org.apache.commons.math3.linear import Array2DRowRealMatrix, ArrayRealVector


COEFF_MATRIX = Array2DRowRealMatrix([
    [-0.313, 56.7, 0.],
    [-0.0139, -0.426, 0.],
    [0., 56.7, 0.]
])
COEFF_VECTOR = ArrayRealVector([0.232, 0.0203, 0.]) 
STATE_VECTOR = ArrayRealVector([0., 0., 0.])
INTEGRATION_STEP = 1e-3

def calculate_angles(input_value):
    """
    Функция выполняет расчёт вектора угловых значений атаки и тангажа
    при продольном движении ЛА.

    Входным аргументом input_value функции является 
    вещественное значение угла руля высоты.
    """

    # В левой части выполняется 
    # умножение матрицы коэффициентов на вектор состояния
    left_part = COEFF_MATRIX.operate(STATE_VECTOR)

    # В правой части выполняется 
    # умножение вектора коэффициентов на значение угла руля высоты 
    right_part = COEFF_VECTOR.mapMultiply(input_value)
    
    # В результате выполняется сложение левой и правой частей
    angles = left_part.add(right_part)

    return angles

• Строка 1. Импорт необходимых структур данных из модуля линейной алгебры библиотеки Apache Commons Math. Матрица Array2DRowRealMatrix и вектор ArrayRealVector вещественных значений будут использоваться далее для расчёта системы уравнений (3), представленных в матричном виде (4);

• Строки 4, 8. Создание объекта матрицы коэффициентов уравнений (4) продольного движения гражданского самолёта. Матрица далее будет использована в функции calculate_angles;

• Строка 9. Создание вектора коэффициентов правой части уравнений (4);

• Строка 10. Создание вектора состояний модели с начальными значениями угла атаки \(\alpha\), угловой скорости тангажа \(q\) и угла тангажа \(\theta\). Далее на этапе моделирования вектор STATE_VECTOR будет обновлять состояние при завершении каждого шага интегрирования;

• Строка 11. Создание значения шага интегрирования (в секундах), ранее указанного при настройке проекта;

• Строки 13, 33. Создание функции calculate_angles с целью решения системы уравнений (4) для расчёта вектора угловых значений атаки и тангажа при продольном движении ЛА. Аргументом функции input_value является вещественное значение \(\delta\) угла руля высоты (см. рисунок выше).

Далее в текстовом редакторе кода блока Jython требуется перейти на вкладку Код расчета для анализа этапа программы, которая выполняется на этапе моделирования.

В рамках этапа моделирования рассчитывается вектор угловых значений атаки и тангажа на каждом шаге интегрирования с начального времени моделирования и до остановки расчёта. Данный вектор состояния обновляется с применением явного метода Эйлера.

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

Код расчета

# Выполняется расчёт вектора угловых значений атаки и тангажа
# при продольном движении ЛА
angles = calculate_angles(in_1)

# Обновляется вектор состояния с учётом шага интегрирования
angles_step = angles.mapMultiply(INTEGRATION_STEP)
STATE_VECTOR = STATE_VECTOR.add(angles_step)

# На выходной порт передаётся значение угла тангажа 
# по индексу №2 вектора состояния
out_1 = STATE_VECTOR.getEntry(2)

• Строка 3. На каждом шаге интегрирования, равном \(1 [мс]\), выполняется расчёт вектора угловых значений атаки и тангажа для актуального входного значения \(\delta\) угла руля высоты, полученного из объекта in_1 входного порта;

• Строки 6, 7. Выполняется обновление вектора состояния STATE_VECTOR с применением явного метода Эйлера.

• Строка 11. В результате выполнения Кода расчета на каждом шаге интегрирования создаётся объект out_1 выходного порта блока Jython, в который выводится актуальное значение угла тангажа \(\theta\).

2.2.  Реализация модели системы управления углом тангажа ЛА

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

Блок	№ блока в проекте	Назначение для модели
Сравнивающее устройство	Блок №2	Выполняет вычитание из уставки угла тангажа расчётное значение \(\theta\) угла тангажа модели, реализованной в блоке Jython. Разность передаётся в ПИД-регулятор системы управления через разветвитель №3.
Разветвитель	Блок №3	Передаёт сигнал от сравнивающего устройства в передаточные звенья №4, 5, 6
Усилитель	Блок №4	Необходим для регулирования пропорционального коэффициента до достижения приемлемой скорости нарастания сигнала
Интегрирующее звено	Блок №5	Предназначен для усиления интегрального коэффициента для снижения установившейся ошибки
Дифференцирующее звено	Блок №6	Необходим для настройки дифференциального коэффициента для уменьшения критических осцилляций в системе
Сумматор	Блок №7	Суммирует результаты передаточных звеньев в рамках реализации ПИД-регулятора

В блоке №6 дифференцирующего звена параметру Таu присваивается малое значение времени, равное \(1\mathrm{e}{-5}\).

Проект содержит набор глобальных переменных в соответствии с таблицей ниже.

Параметр	Имя глобальной переменной	Применяется к блоку
Коэффициент пропорционального звена	PROPORTIONAL_COEFF	блок №4
Коэффициент интегрирующего звена	INTEGRAL_COEFF	блок №5
Коэффициент дифференцирующего звена	DIFFERENTIAL_COEFF	блок №6

Подсказка

Глобальные переменные позволяют удобно управлять настройками проекта в процессе адаптации и отладки модели.

Также глобальные переменные необходимы для взаимодействия внешних программ пользователя с моделью, реализованной на платформе REPEAT VISION. Подробнее о взаимодействии с платформой REPEAT VISION вы можете ознакомиться в разделе уроков.

Результаты первого моделирования системы управления углом тангажа ЛА

Для получения первых результатов моделирования необходимо присвоить значения глобальным переменным ПИД-регулятора в соответствии с таблицей ниже.

Параметр	Имя глобальной переменной	Величина
Коэффициент пропорционального звена	PROPORTIONAL_COEFF	1.0
Коэффициент интегрирующего звена	INTEGRAL_COEFF	0.0
Коэффициент дифференцирующего звена	DIFFERENTIAL_COEFF	0.0

Далее выполните запуск проекта на расчёт и создайте новый график на панели Визуализация.

На график добавьте параметр порта out_1 Jython блока №8 с целью анализа значений \(\theta\) угла тангажа.

(png, hires.png)

Зависимость изменения угла тангажа от времени.

На графике выше приняты следующие обозначения.

Линии	Описание
Синяя линия	Зависимость изменения угла тангажа \(\theta [рад]\) от времени \([c]\)

Для анализа результатов расчётов и адаптации модели системы управления углом тангажа ЛА используются следующие критерии качества отклика:

№	Критерий	Максимальное допустимое значение
1	Перерегулирование \(\sigma\)	0.15
2	Время регулирования \(t_{рег}\) с коридором точности в 0.05	\(10 [c]\)
3	Относительная ошибка стационарного значения угла тангажа \(\varepsilon_{S}\)	0.05

Критерии качества отклика определяются по формулам (7).

(7)\[\begin{split}\begin{aligned} \varepsilon_{S} &= \frac{\theta_{S} - \theta_{T}}{\theta_{T}} \\ \sigma &= \frac{\theta_{max} - \theta_{S}}{\theta_{S}} \end{aligned}\end{split}\]

В уравнении (7) приняты следующие обозначения:

Обозначение	Описание
\(\varepsilon_{S}\)	Относительная ошибка стационарного значения угла тангажа
\(\theta_{S}\)	Расчётное стационарное значение угла тангажа \(\theta\)
\(\theta_{T}\)	Уставка угла тангажа \(\theta = 0.2 [рад]\)
\(\sigma\)	Перерегулирование
\(\theta_{max}\)	Максимальное расчётное значение угла тангажа \(\theta\)

Критерий времени регулирования \(t_{рег}\) с коридором точности в 0.05 изображён на графике ниже.

(png, hires.png)

Критерий времени регулирования с коридором точности.

Далее необходимо провести анализ результатов расчёта модели системы управления углом тангажа ЛА, выполненного с исходными данными, приведённых в таблице выше.

Применяя формулу (7), получаются следующие значения критериев: \(\sigma = 0\), \(t_{рег} = 22.2 [с]\), \(\varepsilon_{S} = -0.012\).

Критерии качества не выполнены. Модель требует настройки коэффициентов системы управления.

Совет

Эмпирический метод настройки коэффициентов модели системы управления может занимать длительное рабочее время.

С целью снижения трудозатрат настройки коэффициентов системы управления и для повышения точности по адаптации модели необходимо воспользоваться оптимизатором.

Описание автоматизированного процесса адаптации модели системы управления углом тангажа ЛА через решение оптимизационной задачи приведено в следующем разделе.

2.3.  Адаптация модели системы управления углом тангажа ЛА

Для решения оптимизационной задачи определения оптимальных коэффициентов ПИД-регулятора системы управления углом тангажа ЛА в проекте используются следующие блоки:

Блок	№ блока в проекте	Назначение для модели
Выходной порт	блок №10	Необходим для передачи в Python программу результатов расчёта \(\theta\) угла тангажа ЛА

В блоке №10 назначено пользовательское название pitch для удобства получения результатов расчёта.

Так же для решения оптимизационной задачи в проекте используются глобальные переменные, подготовленные ранее.

Глобальные переменные используются оптимизатором как параметры целевой функции для решения оптимизационной задачи.

См. также

Подробную информацию о настройке проекта платформы REPEAT VISION для взаимодействия с внешними программами можно получить в разделе уроков.

Постановка оптимизационной задачи

Задача оптимизации Минимизировать целевую функцию вида (8)

(8)\[F(K_p, K_i, K_d) = w_{\sigma} \cdot \sigma^2 + w_S \cdot \varepsilon_S^2 + w_{t_{рег}} \cdot t_{рег}^2,\]

В функции (8) применяются следующие обозначения \( w_{\sigma}, w_S, w_{t_{рег}} \) — весовые коэффициенты для перерегулирования, стационарной ошибки и времени регулирования, позволяющие по-разному учитывать эти компоненты ошибок.

Минимизация целевой функции означает уменьшение ошибок и достижение параметров требуемой работы системы управления.

Дополнительно, после проведения оптимизации проверяются требования критериев качества модели:

\[\sigma \leq 0.15, \quad \varepsilon_S \leq 0.05, \quad t_{рег} \leq 0.4. \]

При решении оптимизационной задачи используются следующие ограничения на переменные:

(9)\[K_p, K_i, K_d \in [0, 2],\]

В итоге имеем задачу оптимизации:

(10)\[\begin{split}\left\{ \begin{aligned} & \min_{K_p, K_i, K_d} F(K_p, K_i, K_d), \\ & \sigma \leq 0.15, \\ & \varepsilon_S \leq 0.05, \\ & t_{рег} \leq 0.4, \\ & K_p, K_i, K_d \in [0, 2]. \end{aligned} \right.\end{split}\]

Реализация оптимизационной задачи на платформе REPEAT VISION

Сначала откройте ноутбук и в новой тетрадке создайте скрипт для подключения к модели через платформу REPEAT VISION.

Далее в первую ячейку добавьте код загрузки Python пакета repeat, необходимого для взаимодействия с платформой REPEAT VISION и для решения оптимизационной задачи.

import piplite
await piplite.install("repeat")

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

import numpy as np
import pandas as pd
import scipy.optimize as so
import matplotlib.pyplot as plt

В следующей ячейке выполните импорт Python пакета repeat.

import repeat as rp

Далее требуется настроить подключение к платформе REPEAT VISION и выбрать модель для решения оптимизационной задачи.

user = rp.User(
  "enter a unique API key from REPEAT VISION user cabinet instead of this comment"
  )
app = rp.Application(user)
# enter a project ID from REPEAT VISION application instead of 12345
project = 12345
model = app.get_exploration_model(project)

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

print(model, model.existing_variables, sep = "\n")

# ExplorationModel(project = 12345, 
# name = '[DEMO] Система управления углом тангажа самолёта', 
# time_interval = None)
# DIFFERENTIAL_COEFF    0.0
# INTEGRAL_COEFF        0.0
# PROPORTIONAL_COEFF    1.0
# dtype: float64

В следующей ячейке выполним подготовку критериев качества отклика в соответствии с формулами (7) и в соответствии с таблицей, содержащей допустимые значения критериев. Критерии качества реализуются в виде методов класса СontrolСriterion.

При инициализации экземпляра класса требуется передать в метод __init__() временной ряд результата расчёта модели, предварительно полученной от платформы REPEAT VISION.

class СontrolСriterion:
    """
    Класс предназначен для определения критериев качества отклика
    
    Parameters
    ----------
    results : pandas.Series
      Временной ряд результата расчёта модели платформы REPEAT VISION
    """
    def __init__(self, results):
        self.results = results
        self.steady_state = results.iloc[-1]
        self.maximum = results.max()
        
    def get_stabilization_time(self, relative_width):
        """
        Метод рассчитывает время регулирования
  
        Parameters
        ----------
        relative_width : float
          Значение доли (от 0 до 0.5) установившейся величины результата расчёта.
  
        Returns
        -------
        stabilization_time : float
          Время регулирования.
        """
        up_bound = 1 + relative_width
        low_bound = 1 - relative_width
        values_up = self.results[self.results > self.steady_state * up_bound]
        values_low = self.results[self.results < self.steady_state * low_bound]
        stabilization_time = np.nanmax([
            values_up.index.max(), 
            values_low.index.max()
        ])
        return stabilization_time

    def get_overshoot_err(self):
        """
        Метод рассчитывает значение перерегулирования
  
        Returns
        -------
        relative_error : float
          Величина перерегулирования.
        """
        relative_error = (
            self.maximum - self.steady_state
        ) / self.steady_state
        
        return relative_error

    def get_steady_state_err(self, target):
        """
        Метод рассчитывает относительную ошибку стационарного значения
  
        Parameters
        ----------
        target : float
          Уставка значения.
  
        Returns
        -------
        relative_error : float
          Относительная ошибка стационарного значения.
        """
        relative_error = (
            self.steady_state - target
            ) / target
        
        return relative_error

• Строки 15, 37. Реализация метода расчёта времени регулирования. Аргументу relative_width передаётся значение 0.05 коридора точности в соответствии с таблицей критериев качества;

• Строки 39, 52. Реализация метода расчёта перерегулирования в соответствии с формулой (7);

• Строки 54, 72. Реализация метода расчёта относительной ошибки стационарного значения в соответствии с формулой (7). Аргументу target передаётся значение уставки угла тангажа в соответствии с настройками блока №1 Константа.

Методы расчёта критериев качества класса СontrolСriterion далее будут использоваться при реализации целевой функции objective_kernel далее.

def objective_kernel(angles, target, relative_width):
    """
    Целевая функция суммы квадратов нормализованных на единицу
    критериев качества системы управления с учётом весов критериев.

    Parameters
    ----------
    angles : pandas.Series
      Временной ряд результатов расчёта угла [рад] тангажа ЛА.
    target : float
      Величина уставки угла [рад] тангажа ЛА.
    relative_width : float
      Значение доли (от 0 до 0.5) установившейся величины результата расчёта.
  
    Returns
    -------
    y : float
      Значение целевой функции.
    """
    angles.index = angles.index / 1e3 # [мс] в [с]
    сrit = СontrolСriterion(angles)
    overshoot_error = сrit.get_overshoot_err()
    steady_state_error = сrit.get_steady_state_err(target)
    stabilization_time = сrit.get_stabilization_time(relative_width)
    stabilization_time_relative = stabilization_time / angles.index[-1]
    weight_overshoot, weight_steady_state, weight_time = 5, 1, 10
    y = (
        weight_overshoot * overshoot_error**2
        + weight_steady_state * steady_state_error**2
        + weight_time * stabilization_time_relative**2
    )
    return y

• Строка 20. Выполняется перевод величины времени из миллисекунд в секунды, поскольку API платформы REPEAT VISION формирует временные ряды с размерностью времени \([мс]\);

• Строки 22, 24. Вычисление критериев качества регулятора в соответствии с формулами (7);

• Строка 25. Нормализация критерия времени на единицу для использования в целевой функции (10);

• Строка 26. Создание постоянных значений весов критериев для успешного решения задачи оптимизации. Веса критериев назначаются инженером эмпирически;

• Строки 27, 31. Вычисление значения целевой функции в соответствии с уравнением (10).

В следующей ячейке необходимо подготовить вспомогательную функцию run() для взаимодействия с платформой REPEAT VISION для передачи variables значений глобальных переменных модели и для последующего получения результатов расчёта временного ряда значений угла тангажа.

def run(model, variables):
    """
    Функция подготовки временного ряда результатов одного расчёта
     
    Parameters
    ----------
    model : repeat.ExplorationModel
      Экземпляр класса ExplorationModel для взаимодействия с платформой REPEAT VISION.
    variables : pandas.Series
      Набор глобальных переменных модели.
  
    Returns
    -------
    angles : pandas.Series
      Временной ряд результатов расчёта угла [рад] тангажа ЛА с 0 до 40 [с].
  
    """
    with model as md:
        md.run(variables)
        angles = md.get_results("pitch",  
                                rp.TimeInterval(start = 0, end = 40.0e3))
        
    return angles

Функция run() возвращает временной ряд значений угла тангажа из блока №10 Выходной порт. Пример результатов расчёта приведён на графике результатов первого расчёта.

В следующей ячейке выполняется подготовка вспомогательной функции для обеспечения взаимодействия с интерфейсом пакета scipy для дальнейшего процесса решения оптимизационной задачи.

def objective(pid_terms, model, target, relative_width):
    """
    Вспомогательная внешняя функция для целевой функции objective_kernel и
    и для функции run. Функция objective обеспечивает взаимодействие
    непосредственно с оптимизационными алгоритмами пакета scipy

    Parameters
    ----------
    pid_terms : numpy.array
      Вектор трёх значений для глобальных переменных модели REPEAT VISION.
    model : repeat.ExplorationModel
      Экземпляр класса ExplorationModel для взаимодействия с платформой REPEAT VISION.
    target : float
      Величина уставки угла [рад] тангажа ЛА.
    relative_width : float
      Значение доли (от 0 до 0.5) установившейся величины результата расчёта.
  
    Returns
    -------
    y : float
      Значение целевой функции.
    """
    variables = pd.Series(dict(
        DIFFERENTIAL_COEFF = pid_terms[0], 
        INTEGRAL_COEFF = pid_terms[1],
        PROPORTIONAL_COEFF = pid_terms[2],
    ))
    angles = run(model, variables)
    y = objective_kernel(angles, target, relative_width)
    return y

• Строки 23, 27. Создаётся объект variables pandas.Series, необходимый для передачи глобальных переменных в платформу REPEAT VISION.

В следующей ячейке настраивается взаимодействие с оптимизатором differential_evolution из библиотеки scipy.optimize и выполняется процесс решения оптимизационной задачи.

solution = so.differential_evolution(
    objective,
    x0 = [0, 0, 1],
    args = (model, 0.2, 0.05),
    popsize = 8,
    maxiter = 4,
    bounds = so.Bounds(
      [0, 0, 0,], 
      [2, 2, 2,]
                  ),
    seed = 5,
)

# Время расчёта 0:00:02.170500 ...

• Строка 2. Первым аргументом оптимизатору передаётся целевая функция objective, которую он будет минимизировать;

• Строка 3. Следующему именованному аргументу x0 передаётся список начальных значений глобальных переменных в строгом порядке, использованном в функции objective: DIFFERENTIAL_COEFF, INTEGRAL_COEFF, PROPORTIONAL_COEFF;

• Строка 4. Аргументу args передаётся набор постоянных параметров и объект model в порядке, назначенном в функции objective: model, target, relative_width;

• Строки 7, 10. Аргументу bounds назначаются ограничения на область вещественных значений поиска оптимальных значений переменных DIFFERENTIAL_COEFF, INTEGRAL_COEFF, PROPORTIONAL_COEFF в соответствии с функциями ограничения целевой функцией (9);

• Строка 11. Аргументу seed назначается величина для обеспечения воспроизводимости результата во время отладки решения задачи оптимизации.

Начальные значения глобальных переменных, переданные аргументу x0, формируются инженером эмпирически. Начальное приближение снижает время поиска оптимального решения.

Ограничения на область поиска оптимальных значений, переданные аргументу bounds, позволяют сосредоточить алгоритм на указанной области пространства параметров. Выбор области поиска оптимальных значений должен опираться на физические ограничения модели и технические требования к модели.

2.4.  Анализ результатов адаптированной модели системы управления

Для выполнения анализа результатов необходимо в следующей ячейке запустить модель системы управления с оптимальными значениями коэффициентов ПИД-регулятора для получения актуальных результатов изменения угла тангажа ЛА.

variables = pd.Series(dict(
        DIFFERENTIAL_COEFF = solution.x[0], 
        INTEGRAL_COEFF = solution.x[1],  
        PROPORTIONAL_COEFF = solution.x[2]  
))
results = run(model, variables)
results.index = results.index / 1e3
axis = result.plot()
axis.set_xlabel(r"Time, $[s]$")
axis.set_ylabel(r"Pitch angle, $\theta, [rad]$")
axis.grid(True)

• Строки 1, 5. Создаётся объект variables с оптимальными значениями глобальных переменных модели, передаваемых из вектора x объекта solution, который был создан в предыдущей ячейке по результату решения оптимизационной задачи;

• Строка 7. Миллисекунды переводятся в секунды для удобства анализа временного ряда угла тангажа;

В результате работы ячейки выше выполняется построение графика,приведённого ниже, для анализа актуальных результатов.

(png, hires.png)

Изменение угла тангажа от времени в адаптированной модели.

Линии	Описание
Синяя линия	Результаты первого расчёта зависимости угла тангажа от времени
Зелёная линия	Изменение угла тангажа от времени в адаптированной модели

Результаты анализа критериев качества системы управления представлены в таблице ниже.

№	Критерий	Результат оптимизации	Максимальное допустимое значение
1	Перерегулирование	0.1293	0.15
2	Время регулирования с коридором точности в 0.05	4.80	\(10 [c]\)
3	Относительная ошибка стационарного значения угла тангажа	0.0034	0.05

3.  Выводы

В рамках данного практического руководства вам удалось освоить:

• Процесс разработки математической модели продольного движения ЛА с помощью блока Jython;

• Процесс создания модели управления углом тангажа ЛА с использованием библиотеки блоков Автоматика;

• Процесс разработки Python программы для автоматизированной настройки коэффициентов ПИД-регулятора системы управления углом тангажа.

4.  Список литературы

---

[1] (1,2)

Zakarya Mohammed Motea, Herman Wahid, Junaid Zahid, Swan Htet Lwin and Abbas Musa Hassan, A Comparative Analysis of Intelligent and PID Controllers for an Aircraft Pitch Control System.