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

Управление ЛА по заданному лучу, это метод навигации, при котором управление движением осуществляется посредством РЛС или другого устройства способного построить направление на заданную точку. Такой метод навигации широко используется в ЛА разных видов (в том числе и БПЛА мультироторного и самолётного типов).

Основные принципы этого метода:

1. Формирование вектора управления:

Вектор формируется в заранее оговоренной системе координат. Обычно этот вектор задаётся в виде вектора скорости или направления, например, по азимуту и углу наклона. В данной статье предполагается, что вектор формируется с помощью РЛС установленной на земле.

2. Планирование управления:

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

Структурные схемы используемые в статье изначально строились с помощью Simulink также с помощью matlab были реализованы алгоритмы минимизации для определения коэффициентов передаточных функций. В рамках работы по импортозамещению, данная САУ ЛА была построена в среде Engee, были получены частотные и временные характеристики контуров управления и построена общая САУ, включающая в себя все исследуемые контура управления.

В данной статье используется осесимметричный ЛА, поэтому каналы продольного и бокового движения идентичны.

Передаточная функция угловой скорости тангажа представлена в уравнении:

Передаточная функция по перегрузке представлена в уравнении:

Где:

Константы, характеризующие динамику ЛА:

Уравнение изолированного движения по крену:

b1 = 50, b3 = 6

Математические модели исполнительных приводов:

Структурная схема астатической САУ нормальной перегрузкой представлена на рисунке 3, для боковой перегрузки схема имеет аналогичный вид.

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

В данном законе управления входит интеграл по ошибке управления (интеграл от разности фактического и заданного ускорения), это позволяет устранить статические ошибки при отработке команды и при действии возмущающего момента Mzb

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

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

Числитель:

Знаменатель:

Коэффициенты знаменателя передаточной функции:

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

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

1. Простота реализации (например, для дальнейшей доработки и модификации)

2. Меньшие требования к формулировке задачи (не обязательно задавать точную границу или сложные ограничения)

3. Быстродействие (в случаях, когда функция легко дифференцируема и не имеет сложных ограничений, градиентный спуск находит решение существенно быстрее встроенного метода)

Для того чтобы использовать методы оптимизации необходимо на базе полученных коэффициентов полиномов для числителя и знаменателя передаточной функции САУ записать целевую функцию, которая будет минимизироваться и функцию ограничений, внутри которых будет происходить поиск минимума. В частности, для контура управления нормальной перегрузкой будем минимизировать функцию трех неизвестных которыми будут коэффициенты усиления

Минимизация времени нарастания переходного процесса будет происходить при максимизации величины z = a0*(a1)^-1. Ограничение накладываемое на перерегулирование при переходном процессе

можно заменить неравенством:

а устойчивость САУ гарантируется выполнением достаточных условий:

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

min(-z) = min(-a0*(a1)^-1)

При ограничениях:

Engee скрипт, реализующий метод градиентного спуска представлен в листинге 1

Листинг 1 – скрипт Engee для оптимизации коэффициентов передаточной функции

# подключение библиотек import Pkg Pkg.add("ControlSystems") ## библитека для анализа САУ Pkg.add("Plots") ## библитека для построения графиков using ControlSystems using Plots using LinearAlgebra # библиотека для операций с векторами и матрицами (используется в методе градиентного спуска) # Ограничения верхнего уровня (для расчета оптимальных коэффициентов передаточной функции) lu = fill(300.0, 3) # Ограничения нижнего уровня (для расчета оптимальных коэффициентов передаточной функции) lb = fill(0.3, 3) # Начальная точка с которой начинаем итерации по методу градиентного спуска x0 = [0.5, 0.5, 0.5] # Функция целей function myfun(x::Vector{T}) where {T} return (150850000*x[1]) / (150850000*x[2] - 168000*x[1] + 67882500*x[3] + 49005000) end # Функция ограничения неравенства function mycon(x::Vector{T}) where {T} c1 = 2.15 - (150850000*x[1])^-1 (150850000x[2] - 168000*x[1] + 67882500*x[3] + 49005000) (73687500x[3] - 150000*x[1] - 168000*x[2] + 2009250) inv(588909 - 150000x[2]) c2 = 2.15 - (150850000*x[2] - 168000*x[1] + 67882500*x[3] + 49005000)^-1 (73687500x[3] - 150000*x[1] - 168000*x[2] + 2009250) (588909 - 150000x[2]) inv(9534) c3 = 2.15 - (73687500x[3] - 150000*x[1] - 168000*x[2] + 2009250)^-1 (588909 - 150000x[2]) 9534 inv(75) return [c1, c2, c3] end # Градиент функции цели (оцениваем численно) function grad_f(x::Vector{T}, eps=1e-6)::Vector{T} where {T} n = length(x) # расчет длины вектора х gradient = zeros(T, n) # вектор для хранения значений градиента for i in 1:n dx = copy(x); dx[i] += eps # перебор значений х df = myfun(dx) - myfun(x) # расчет разности значений функции в ближайщих точках х gradient[i] = df / eps # расчет значений градиента end return gradient end # Алгоритм градиентного спуска function minimize_gradient_descent(f, x0, max_iter=100, learning_rate=0.1, tol=1e-8) x = x0 # начальное значение х current_val = f(x0) # начальное значение функции prev_val = Inf # начальное значение критерия останова iter_count = 0 ## начальное значение счетчика while true ## бесконечный цикл current_val = f(x) # текущее значение функции if abs(current_val - prev_val) < tol || iter_count >= max_iter ## критерий остановки итераций break end g = grad_f(x) # расчет градиента функции в точке step = -learning_rate .* g # расчет длины шага # Проверяем границы (ограничения) proposed_x = clamp.(x .+ step, lb, lu) # проверка того попадает ли новая точка в интервал # между заданными ограничениями x .= proposed_x # обновление значения х prev_val = current_val # обновление значения функции iter_count += 1 # увеличение счетчика итераций end return x, current_val end # Запускаем процесс минимизации x_optimal, fval = minimize_gradient_descent(myfun, x0) println("Оптимизированный вектор х:", x_optimal) println("Значение целевой функции:", fval) Kwz = x_optimal[1] Kny = x_optimal[2] Knyi = x_optimal[3] K_iny=0,9683; K_ny=1,54; K_ωz=0.3. c1 = 1.12; c2 = 86; c3 = 131; c4 = 1; c5 = 0; c6 = 21.8; c9 = 0.12; g = 9.81; a0 = c2+c1*c4; a1 = (c1+c4+c5) b1 = 50 b3 = 6 K = b3/b1 T = 1/b1 K1 = 1 K2 = 1 a5 = 75 a4 = 9534 a3 = (588909 - 150000*Kny) a2 = (73687500*Kwz - 150000*Knyi - 168000*Kny + 2009250) a1 = (150850000*Kny - 168000*Knyi + 67882500*Kwz + 49005000) a0 = 150850000*Knyi W2 = tf([-150000*Knyi, -168000*Knyi, 150850000*Knyi],[a5, a4, a3, a2, a1, a0])

Результат расчета передаточной функции:

TransferFunction{Continuous, ControlSystemsBase.SisoRational{Float64}}

-145256.38691773862s^2 - 162687.15334786725s + 1.460795064436058e8

-----------------------------------------------------------------------------------------------------------------

75.0s^5 + 9534.0s^4 + 357783.7499868162s^3 + 2.3711383333067495e7s^2 + 3.016420226099106e8s + 1.460795064436058e8

y, t = step(W2,0:0.1:10) # расчет переходного процесса t = collect(0:0.1:10) # вектор времени для графика переходного процесса plot(t,y',title="Переходный процесс") ## построение графика переходного процесса

По переходному процессу САУ перегрузкой можно сделать вывод что перерегулирование при переходном процессе отсутствует, а время переходного процесса равно 6 с.

y2, t = impulse(W2,0:0.1:10) # расчет точек импульсной характеристики t = collect(0:0.1:10) plot(t,y2',title="Импульсная характеристика") # построение графика импульсной характеристики

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

bodeplot(W2) # построение графика частотной характеристики

nyquistplot(W2) # годограф Найквиста

На рисунке 6 видно как годограф Найквиста не охватывает точку -1 на действительной оси, из этого можно сделать вывод о том что система устойчива.

pzmap(W2) # диаграмма нулей и полюсов

По расположению полюсов видно, что положительные полюсы в системе отсутствуют (система устойчива) и колебательная составляющая полюсов сильно меньше апериодической

marginplot(W2) # диаграмма запасов устойчивости

На диаграмме запасов устойчивости показано значение запаса устойчивости по амлитуде (в абсолютных единицах и ниже пересчитано в дБ), запас устойчивости по фазе не ограничен.

Пересчитаю значение запаса устойчивости в логарифмический масштаб

L = 20*log10(48)

print("Запас устойчивости: ")

print(L)

print(" Дб")

Запас устойчивости:33.624824747511745 Дб

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

Также с САУ управления по лучу присутствует контур стабилизации угловой скорости вращения ЛА, в рамках данной статьи рассматривается один канал измерения угловой скорости (продольная ось) в рамках данной статьи контур управления угловой скоростью синтезироваться не будет, но будет проведен анализ его частотных и временных характеристик

Листинг 2– скрипт Engee для частотного и временного анализа контура управления угловой скоростью ЛА

using ControlSystems W1 = tf(7500,[1, 125, 7500]) W2 = tf(0.12,[0.02, 1]) W3 = 9.26*W1*W2 W4 = feedback(W3,1) W5 = tf(8.21,[1, 0])*W4 W6 = feedback(W5,1) y, t = step(W6,0:0.1:2) # расчет переходного процесса t = collect(0:0.1:2) plot(t,y',title="Переходный процесс") ## построение графика переходного процесса y2, t = impulse(W6,0:0.1:10) # расчет точек импульсной характеристики t = collect(0:0.1:10) plot(t,y2',title="Импульсная характеристика") # построение графика импульсной характеристики nyquistplot(W6) # годограф Найквиста pzmap(W6) # диаграмма нулей и полюсов marginplot(W6) # диаграмма запасов устойчивости

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

TransferFunction{Continuous, ControlSystemsBase.SisoRational{Float64}}

68422.14000000001

----------------------------------------------------------

0.02s^4 + 3.5s^3 + 275.0s^2 + 15834.0s + 68422.14000000001

Также стоит заметить что по умолчанию в Engee для интегрирования моделей используется метод Эйлера с постоянным шагом. Для некоторых моделей при данном методе интегрирования решение получается расходящимся (как например и для данной схемы). Для того чтобы добиться устойчивости при интегрировании структурной схемы шаг интегрирования выбрал переменный и метод интегрирования Feagin14 (метод Фиджина четырнадцатого порядка, вариант метода Рунге — Кутты)

Для построения САУ ЛА по лучу требуется в структурную схему рис. 19 встроить контура управления перегрузкой и угловым движением относительно продольной оси. Структурная схема, в которой присутствуют вышеуказанные контура показана на рисунке 20.

В составе денной схемы присутствуют контура управления угловым движением ЛА и управления нормальной перегрузкой. Переходный процесс для данной схемы занял 20 с, перерегулирование в данном переходном процессе составляет 40%. Данные значения связаны с высокой массой исследуемого ЛА и необходимостью в уточнении коэффициентов регулятора после объединения разных контуров управления в одну САУ.