Научная статья на тему 'ФОРМИРОВАНИЕ АЛГОРИТМОВ СИСТЕМЫ АВТОМАТИЧЕСКОГО УПРАВЛЕНИЯ ПРЕОБРАЗУЕМОГО БЕСПИЛОТНОГО ЛЕТАТЕЛЬНОГО АППАРАТА'

ФОРМИРОВАНИЕ АЛГОРИТМОВ СИСТЕМЫ АВТОМАТИЧЕСКОГО УПРАВЛЕНИЯ ПРЕОБРАЗУЕМОГО БЕСПИЛОТНОГО ЛЕТАТЕЛЬНОГО АППАРАТА Текст научной статьи по специальности «Механика и машиностроение»

CC BY
69
20
i Надоели баннеры? Вы всегда можете отключить рекламу.
Журнал
Труды МАИ
ВАК
Ключевые слова
СИСТЕМА АВТОМАТИЧЕСКОГО УПРАВЛЕНИЯ / ПИД-РЕГУЛЯТОР / ЛИНЕЙНО-КВАДРАТИЧНЫЙ РЕГУЛЯТОР / НЕЛИНЕЙНАЯ ЛОГИКА / ОБРАТНАЯ ЗАДАЧА ДИНАМИКИ

Аннотация научной статьи по механике и машиностроению, автор научной работы — Аполлонов Дмитрий Вадимович, Бибикова Кристина Игоревна, Шибаев Владимир Михайлович, Ефимова Ирина Евгеньевна

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

i Надоели баннеры? Вы всегда можете отключить рекламу.

Похожие темы научных работ по механике и машиностроению , автор научной работы — Аполлонов Дмитрий Вадимович, Бибикова Кристина Игоревна, Шибаев Владимир Михайлович, Ефимова Ирина Евгеньевна

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.
i Надоели баннеры? Вы всегда можете отключить рекламу.

CREATION OF ALGORITHMS FOR THE AUTOMATIC CONTROL SYSTEM OF THE CONVERTIBLE UNMANNED AERIAL VEHICLE

Convertible unmanned aerial vehicle called a tiltrotor is considered. Tiltrotor can flight airplane mode and helicopter mode by changing the propulsion thrust vector by tilting the axes of rotation of the engines. The tiltrotor has the ability to make a vertical take-off and land, it can also perform horizontal flight at high-speed cruise. It is difficult to develop tiltrotor control system due to the variable characteristics of the UAV during flight. This paper describes analysis of methods for selecting controller parameters of control system regulators. When analyzing fault tolerance ensuring controllability imposes requirements on the choice of control laws in the angular motion control loop and trajectory motion control loop. In addition to the widely used PD- and PID-controllers and the linear-quadratic controller, algorithms based on the concept of the invers dynamic problem and fuzzy logic are considered. Results of research on the creation of algorithms for the automatic control system of a convertible unmanned aerial vehicle are presented. All the described methods were analyzed to ensure the robustness of the automatic control system in order to ensure the necessary quality of control in cases where the control object differs from the calculated one. In case of failure of individual elements of the system the characteristics of control object change during operation, this problem is also considered in this paper. A comparative assessment of various methods for implementing control loops for UAV trajectory and angular motion using proportional-integro-differential, linear-quadratic, nonlinear controllers and the method of the inverse dynamic problem is carried out to identify the advantages and disadvantages of all these methods.

Текст научной работы на тему «ФОРМИРОВАНИЕ АЛГОРИТМОВ СИСТЕМЫ АВТОМАТИЧЕСКОГО УПРАВЛЕНИЯ ПРЕОБРАЗУЕМОГО БЕСПИЛОТНОГО ЛЕТАТЕЛЬНОГО АППАРАТА»

Труды МАИ. 2022. № 122 Trudy MAI, 2022, no. 122

Научная статья УДК 519.711.2

DOI: 10.34759/Ы-2022-122-23

ФОРМИРОВАНИЕ АЛГОРИТМОВ СИСТЕМЫ АВТОМАТИЧЕСКОГО УПРАВЛЕНИЯ ПРЕОБРАЗУЕМОГО БЕСПИЛОТНОГО ЛЕТАТЕЛЬНОГО АППАРАТА

Дмитрий Вадимович Аполлонов1, Кристина Игоревна Бибикова2®, Владимир Михайлович Шибаев3, Ирина Евгеньевна Ефимова4

1,2,3,4Центральный аэрогидродинамический институт имени Н.Е. Жуковского, ЦАГИ,

Жуковский, Московская область, Россия

2ccfstd@tsagi.ru:

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

Для цитирования. Аполлонов Д.В., Бибикова К.И., Шибаев В.М., Ефимова И.Е. Формирование алгоритмов системы автоматического управления преобразуемого беспилотного летательного аппарата // Труды МАИ. 2022. № 122. DOI: 10.34759/trd-2022-122-23

CREATION OF ALGORITHMS FOR THE AUTOMATIC CONTROL SYSTEM OF THE CONVERTIBLE UNMANNED AERIAL VEHICLE

Dmitry V. Apollonov1, Kristina I. Bibikova20, Vladimir M. Shibaev3, Irina E. Efimova4

U'3'4Central Aerohydrodynamic Institute named after N.E. Zhukovsky, TsAGI, Zhukovsky,

Moscow Region, Russia 2ccfstd@tsagi.com0

Abstract. Convertible unmanned aerial vehicle called a tiltrotor is considered. Tiltrotor can flight airplane mode and helicopter mode by changing the propulsion thrust vector by tilting the axes of rotation of the engines. The tiltrotor has the ability to make a vertical take-off and land, it can also perform horizontal flight at high-speed cruise. It is difficult to develop tiltrotor control system due to the variable characteristics of the UAV during flight. This paper describes analysis of methods for selecting controller parameters of control system regulators. When analyzing fault tolerance ensuring controllability imposes requirements on the choice of control laws in the angular motion control loop and trajectory motion control loop. In addition to the widely used PD- and PID-controllers and the linear-quadratic controller, algorithms based on the concept of the invers dynamic

problem and fuzzy logic are considered. Results of research on the creation of algorithms for the automatic control system of a convertible unmanned aerial vehicle are presented. All the described methods were analyzed to ensure the robustness of the automatic control system in order to ensure the necessary quality of control in cases where the control object differs from the calculated one. In case of failure of individual elements of the system the characteristics of control object change during operation, this problem is also considered in this paper. A comparative assessment of various methods for implementing control loops for UAV trajectory and angular motion using proportional-integro-differential, linear-quadratic, nonlinear controllers and the method of the inverse dynamic problem is carried out to identify the advantages and disadvantages of all these methods. Keywords: automatic control system, PID-controller, linear-quadratic controller, fuzzy logic, the inverse problem of dynamics

For citation: Apollonov D.V., Bibikova K.I., Shibaev V.M., Efimova I.E. Creation of algorithms for the automatic control system of the convertible unmanned aerial vehicle. Trudy MAI, 2022, no. 122. DOI: 10.34759/trd-2022-122-23

1. Введение

Беспилотные летательные аппараты все шире используются для решения

задач различного характера. В зависимости от предполагаемых сценариев

применения разрабатываются проекты беспилотных летательных аппаратов (БЛА)

различных типов, схем и аэродинамических компоновок. В частности, БЛА

преобразуемого типа предлагается использовать в тех случаях, когда БЛА традиционных схем не могут быть применены, например, если требуется обеспечить относительно высокие скорость и дальность полета при отсутствии взлетно-посадочных полос [12,22].

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

Решение задачи формирования сил и моментов, потребных для управления перемещением и ориентацией выбранного БЛА в зависимости от режима полета и конфигурации БЛА была решена в статье «Вопросы выбора архитектуры системы автоматического управления преобразуемого беспилотного летательного аппарата -«конвертоплана» [3].

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

2. Модель динамики полета БЛА

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

БЛА. При этом оси вращения двигателей направлены вверх перпендикулярно относительно продольной оси летательного аппарата.

Движение БЛА в выбранном режиме можно описать следующими уравнениями движения, аналогичными уравнениям движения квадрокоптера [18]:

X = (sin -ф sin у — COS -ф sin тд cos у) —

т.

те

Y = — q + (cos д cosy)—

т

тЕ

Z = (cos -ф sin у + sin -ф sin д cos у) —

т.

IZZ — ¡YY , jTP _ , Мх Ых=----+ — + —

lxx 1ХХ 1ХХ

¡XX — ¡ZZ , My ¡TjVy

0)у =----+ ---муушу

'YY lYY

. _ iYY — ¡XX jTP mZ

ú)z =----ШхШу--—<tixLL+ -—

v lzz lzz lzz

Подъемная сила TE = Y,f.=1 Ti создаётся за счёт работы винтомоторной группы, ее значение равно сумме тяг всех четырех винтокольцевых движителей вентиляторного типа на концах консолей крыльев (рис. 1)

Поперечное управление (по углу крена у) производится аналогично БЛА квадрокоптерной схемы путём формирования управляющей составляющей момента Мх относительно продольной оси за счет изменения соотношения суммарной тяги двигателей, установленных на левых (Т1 и Т4) и правых (Т2 и Т3) консолях крыльев[6,14,15,16,20]. Аналогично, для продольного управления (Mz) - по углу тангажа 0 - путём изменения соотношения суммарной тяги двигателей, установленных перед центром масс БЛА (Т1 и Т2) и позади него (Т3 и Т4).

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

(управляющая составляющая Му) в «вертолетном» режиме БЛА рассматриваемой

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

Рисунок 1 - Схема управления БЛА при выполнении "вертолётного" режима

полета

Линеаризуя приведенные выше уравнения, можно привести их к следующему

виду:

' X = -дд •• те

ч = --д

т

2 = ?д

м.

0)х =

X

'XX

м.

0)у = у

у

УУ

М,

0)2 =

Исходя из этих уравнений, можно сформировать законы управления суммарной тягой и моментами, потребными для стабилизации углового положения

и высоты в рассматриваемом режиме [7,21,23].

3. Алгоритмы системы управления

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

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

Рисунок 2 - Структурная схема САУ преобразуемого БЛА при выполнении

"вертолётного" режима

В контуре управления угловым движением в зависимости от заданных и измеренных значений углов и угловых скоростей формируются потребные значения моментов крена, тангажа и рыскания. Полученные значения моментов и суммарной тяга двигателей передаются в блок формирования потребных отклонений органов управления, которые используются в качестве входных управляющих воздействий для блока расчета динамики полета БЛА. Сигналы текущих значений параметров линейного и углового движения БЛА, полученные с помощью системы бортовых измерений, поступают в контуры управления траекторным и угловым движением БЛА [17].

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

стандартных ПД- и ПИД-регуляторов будут также рассмотрены линейно-квадратичный регулятор (ЛКР-регулятор), регулятор на базе правил нечеткой логики (НЛ-регулятор) и регулятор, сформированный на основе концепции обратных задач динамики (ОЗД-регулятор).

3.1. Контур стабилизации высоты

Формирование контура управления высотой выполним на основании линеаризованного уравнения вертикального движения БЛА: Исходное уравнение:

TE

Y = —g + (cos д cos y) —

m

заменим линеаризованным уравнением, считая отклонения углов крена и тангажа небольшими:

•• Те Y = — — д т

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

ТЕ = ^тек —"зад) + *Л + ф™ —"зад)^

9

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

3.2. Контур управления перемещением в горизонтальной плоскости

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

X = (sin Ф sin Y — cos ф sin д cos y) —

те

Z = (cos Ф sin Y + sin ф sin д cos y) —

Для определения законов управления параметрами линеаризуем данные уравнения движения, тогда при ф = 0 имеем:

1 = —дд

Z = Yd

В этом случае в качестве законов управления координатами X и Z может быть использован ПИД-регулятор, закон управления которого имеет следующий вид:

^зад (^зад ^тек)Ьр УхдЬй + ^ J С^зад

Ниже на примере канала управления перемещением вдоль оси ОХ будут рассмотрены другие варианты контроллеров.

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

Рисунок 3 - Структурная схема контура управления координатой X

На вход контуров поступают заданные значения координат X и 7, которые сравниваются с измеренными значениями в текущий момент времени. На выходе

контуров будут формироваться заданные углы тангажа и крена.

3.3. Пересчёт углов тангажа и крена при произвольной ориентации БЛА

В случае, если продольная ось ОХ связанной с БЛА системы координат не совпадает с осью ОХё земной системы координат, заданные значения углов тангажа дзад и крена узад пересчитываются с учетом угла рыскания

дзаА = дсоБУ + уБтУ,

узад = усоБ^ — дБтУ,

где д, у - углы тангажа и крена соответственно, рассчитанные в предположении, что угол рыскания Ф = 0.

Рисунок 4 - Влияние угла рыскания на углы крена и тангажа

3.4. Контуры управления углами тангажа и крена

Сформированное в контуре управления координатой X значение потребного угла тангажа поступает на вход внутреннего контура - контура стабилизации угла тангажа. Уравнения углового движения относительно поперечной оси 07:

. мг

* г г

д =

При этом на выходе контура будет формироваться соответствующее заданное значение вращающего момента Мгвокруг оси ОЪ связанной системы координат, которое впоследствии будет подано на вход блока динамики БЛА:

Мгзад = Ьр(дТек - дзад) + ЫгКй,

Где Кр, Кй - коэффициенты регулирования угла тангажа.

Сформированное в контуре управления координатой Ъ значение потребного угла крена поступает на вход внутреннего контура - контура стабилизации угла крена. Уравнения углового движения относительно поперечной оси ОХ:

Мх

и>х = ух х

У = ^х

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

Мхзад = Кр(уТек - Узад) + ™хКф

Кр, Кй - коэффициенты обратных связей.

3.5. Контур управления углом рыскания

Линеаризованные уравнения углового движения относительно оси 0У имеют

вид:

Му

1УУ

Закон управления углом рыскания формируется с помощью ПД -регулятора. На вход контура управления углом рыскания ¥ поступает заданное значение угла рыскания Фзад. На выходе регулятора формируется при этом потребный момент рыскания:

МУзал = Кр(утек — Узад) + ШуКа,

Где Кр, Кй - коэффициенты регулирования угла рыскания.

4. Выбор параметров ПИД-регулятора

Рассмотрим более подробно вопрос выбора закона управления и параметров контроллеров, использованных при формировании законов управления САУ БЛА.

Рисунок 5 - Блок-схема ПИД-регулятора

Среди вышеперечисленных законов управления были рассмотрены два случая: случай использования пропорциональной и дифференциальной обратных связей (ПД-регулятор) и случай использования пропорциональной, интегральной и дифференциальной обратных связей (ПИД-регулятор). Выбор значений коэффициентов обратных связей проводился с помощью эталонной модели передаточной функции звена 2-го порядка вида

у

' ' 'ГГЧГ

га"

тек „2

р2 + 2$шр+ш2 зад

характеристики переходного процесса, которого необходимо воспроизвести.

1. ПИД-регулятор

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

2

у —_у

Лтек (р2 + 2%шр + ш2)Лзад

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

рКр + К, Лтек р3 + Кар2 + крГ + к*зад

или, после представления знаменателя в виде произведения двух полиномов:

К-п

к^р + Ъ

V =_21_У

Лтек (р2 + 2&р + ы2)(р + ш±)Лзад

Сравнивая выражения в знаменателях

р3 + Кйр2 + Крр + К1

и

(р2 + 2^шр + ш2)(р + ш-1) = р3 + р2(ш1 + + р(ш2 + 2+ ш2ш1, можем рассчитать коэффициенты ПИД-регулятора:

Кй = +

Кр = (ш2 +

= ш2ш1

Значения частоты ш и коэффициента демпфирования % выбираем равными

значениям этих параметров для эталонного звена второго порядка, значение частоты

среза ш1 апериодического выбираем из условия ш1 = (2 .„5) • ш.

16

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

1 Кр

тгр + 1 ] к

2. ПД-регулятор

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

Х - кр Х

тек п I ы ^ I ы зад

р2 + Кар + Кр д

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

1

Кр -Т2

X, м

о.а

0.6

0.4

0.2

■РЮ ■ рр

10

Рисунок 6 - Переходной процесс для случаев использования ПИД- и ПД-

регуляторов

Из данных примеров следует, что можно представить динамическую модель таким образом, что она будет иметь характеристики аналогичные характеристикам звена 2-го порядка. При этом динамика процесса с использованием ПИД -регулятора и ПД-регулятора будет отслеживаться практически идентично. Однако при каком-либо внешнем возмущении, генерирующем статическую ошибку, интегральное звено ПИД-регулятора будет эту ошибку компенсировать [10].

В контуре управления перемещением вдоль оси координат ОХ, который был настроен с использованием ПИД-регулятора, был проведён анализ устойчивости с помощью логарифмической амплитудно-фазовой частотной характеристики. (ЛАФЧХ) (рис. 7)

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

4-135 -Я

а.

и

п ся

-270

10"2 10"1 10° 101 102 103

Частота, рад/с

Рисунок 7 - ЛАФЧХ при настройке контура с помощью ПИД -регулятора

Запас устойчивости по амплитуде: 15,3 дБ. Запас устойчивости по фазе: 58 град.

5. Линейно-квадратичный регулятор

Линейно-квадратичный регулятор (ЛКР) обеспечивает оптимальное

управление замкнутой системой при минимизации квадратичного функционала

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

19

качества [2]. В отличие от ПИД-регулятора, при формировании сигналов обратных связей квадратичный регулятор не вычисляет производные, а использует переменные состояния системы, отслеживаемые наблюдателем (получаемые с помощью системы бортовых измерений).

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

х = Ах + Ви,

где и -вектор управления; х - вектор состояния.

Критерий оптимальности имеет вид:

^Г. от

(xтQx + итНи)&, где

о

Слагаемое fОT(xтQx)dt является интегральной квадратической ошибкой и характеризует качество управления, здесь Q - неотрицательно-определенная весовая матрица;

г СО у

^ (и1Ки)& — это слагаемое, характеризующее взвешенную «энергию» управления, которое ограничивает управление. Требуемое ограничение на управление может быть обеспечено соответствующим выбором матричной функции Я, где Я - положительная весовая матрица.

Оптимальное управление при этом имеет вид:

и = —Я-1ВтКх ,

где матрица коэффициентов обратной связи К определяется из уравнения Рикатти:

К = -КА - АТК + КВИ-1ВТК - £

Решение уравнений Рикатти и расчет матрицы коэффициентов К линейно-квадратичного оптимального регулятора можно выполнить с помощью библиотечной функции МЛТЬЛВ:

к = ^г(А,в,д,К)

Рассмотрим подробнее настройку системы автоматического управления с помощью ЛКР на примере контура управления координатой X.

Линеаризованное уравнение движения квадрокоптера вдоль оси ОХ имеет

вид:

Х = -дд

или

X = V У = -дд

Запишем эту систему равнений движения в матричном виде:

х = Ах + Ви,

где

А =

0 1 00

, В =

0

а

Выберем весовые матрицы по правилу Брайсона [1], тогда получим:

Q =

-тг 0

X.

м

0 -4т

V,

м

, я =

1

Где в знаменателях диагональных элементов указаны характерные масштабы

диапазона изменения соответствующей переменной.

Ъ сек

Рисунок 8 - Настройки ЛКР при различных параметрах весовых матриц

1. Базовый случай, выбраны параметры Хм = 0,1 м, Ум = 0,1 м/с, Эм = 0,01 рад. Весовые матрицы, оцененные по правилу Брайсона, имеют вид:

Q =

100 0 0 100

, Я = |10000|.

2. Случай ограничения перемещения по оси ОХ, Хм = 0,01 м

Q =

10000 0 0 100

Я = |10000|.

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

3. Случай, когда ограничена скорость перемещения Ум = 0,01 м/с.

Q =

100 0 0 10000

Я = |10000|.

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

4. Случай ограничения управления тангажа, Эм = 0,001 рад. Весовые

матрицы будут равны Я = 110000001, Q =

100 0 0 100

. Переходной процесс в этом

случае будет выходить на заданное значение без перерегулирования с допустимым временем срабатывания.

Используя линейно-квадратичный регулятор, в контуре управления

координатой X был проведён анализ устойчивости с помощью ЛАФЧХ. (рис. 9)

50

Рисунок 9 - ЛАФЧХ при настройке контура с помощью ЛКР

6. Формирование алгоритма управления на основе концепции обратных задач

динамики

Задача формирования закона управления динамической системой может рассматриваться как обратная задача динамики (ОЗД) [8,13,24]: назначается «траектория» движения фазовой координаты х*(0, рассчитывается «сила», необходимая для того, чтобы система двигалась по заданной траектории, строится управление и*, обеспечивающее создание необходимой управляющей силы.

Так же, как и рассмотренные выше подходы к настройке САУ, формирование структуры и выбор параметров закона управления рассмотрим на примере контура управления координатой х.

Считаем, что процесс движения вдоль оси ОХ описывается уравнением второго порядка вида

х = Р(х,х,д)

где Р(х,х,д) - функция, точно нам не известная, описывающая зависимость силы, действующей на летательный аппарат вдоль оси Х, от его положения и скорости, а также от угла тангажа д, являющегося в данном случае управляющим параметром, определяющим проекцию суммарной силы тяги винтов на горизонтальную ось.

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

_ 1

% = р2 + ^р + 0О %зад

или дифференциальным уравнением

X = Р'(Х,Х) = Ро(Хзад -х)- ргх

Параметрами, определяющими силу, вызывающую движение эталонной модели, являются фактическое положение системы - координата х, желаемое положение хзад и горизонтальная скорость ¥х = х.

Задача формирования алгоритма управления - обеспечить значения ускорения (или силы) такое же, как у эталонной модели, в данном случае - звена второго порядка с заданными параметрами. Для этого формируется управляющее отклонение угла тангажа

АадЮ = к[р'(х,х)-х]

Или, после интегрирования обеих частей,

$зад(£) = к I Р'(х,х)& — к(х — х1=0)

Преобразуя с учётом уравнения движения эталонной модели получим:

= кРо I (*зад — х)<И — крг(х — хг=о) — к(х — х=), 0

Где к - коэффициент усиления,

Ро,^! - коэффициенты, определяющие характер протекания назначаемой траектории движения:

С у С у

тх - постоянная времени эталонной модели (звена второго порядка), { - коэффициент демпфирования эталонной модели (звена второго порядка). Главное требование, на основе которого должен выполняться расчет коэффициента усиления к, заключается в том, чтобы фактическая траектория

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

x(t) = x*(t).

Поскольку величина к определяет динамику контура ускорения, то условие возможно большего приближения процессов x(t) ^ x*(t) будет реализовано в том случае, когда быстродействие контура ускорения существенно выше быстродействия внешнего контура. Для оценки рекомендуемого значения величины к рассмотрим линеаризованное уравнение динамики движения управляемого объекта в виде

х = F0^ Ах + Ах + F^ • (-Ад), F-$ > 0, для рассматриваемого случая: = д.

Тогда остальные уравнения алгоритма управления представим в виде:

А-0зад(Х) = k[F'(Ax,Ax) - Ах], F'(Ax,Ax) = -р0 • Ах - ргАх.

При этом динамика контура ускорения может быть оценена по уравнению свободного движения

Аах + (к • F-$ — F1) • Аах = 0, Откуда получаем оценку для постоянной времени контура ускорения

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

Тп = -.

а KF-e- Fi

1

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

?а= N >1.

а

следовательно, коэффициент усиления может быть оценен, как

к =

Из условия устойчивости замкнутого контура управления линеаризованной динамической системой с выбранным коэффициентом усиления к (по критерию Гурвица) получаем

1 . гт

Для случая перемещения на режиме, близком к висению, считая Р0 =^=0,

получим N > V2. Это минимальное значение, обеспечивающее устойчивость линеаризованной системы. В общем случае увеличение значения N и, следовательно, коэффициента усиления к обеспечивает устойчивость замкнутого контура системы управления при условии, что допустимый диапазон изменения управляющего сигнала (в данном случае - тангажа) позволяет реализовать рассчитанное управляющее воздействие.

Используя вышеприведённые уравнения, создана математическая модель канала управления координатой X.

На основании вышеизложенного метода настройки САУ, в контуре управления координатой X был проведён анализ устойчивости с помощью ЛАФЧХ (рис. 10).

Рисунок 10 - ЛАФЧХ при настройке контура с помощью ОЗД

7. Нечёткий регулятор

В качестве ещё одного формирования системы управления летательного аппарата рассмотрим использование алгоритм, построенный на основе принципов нечёткой логики (НЛ) [4,5,9,11,19,25]. Будем использовать нечёткий регулятор, работающий по принципу алгоритма Мамдани.

Рисунок 11 - Функциональная схема системы автоматического управления на

базе нечёткой логики

Алгоритм Мамдани математически описывается следующим образом:

1. Нечеткость (процесс фаззификации): нахождение степени истинности для условий:

V1 Ы!), V2 (и1), V1 (и2), V2 (и2),

где ^.1(и1), ц.2(и1) - функции принадлежности для переменной щ, ^(и!), ц.2(и2) - функции принадлежности для переменной и2.

2. Нечеткий вывод: находятся степени истинности для условий каждого из правил:

А= ^(и1) - ^(и*2) В= ц2(и2) - ц.2(и2)

где через "-" обозначена операция логического минимума, затем находятся усеченные функции принадлежности для выходной переменной ис:

где ^(и), ^¿(и) - функции принадлежности для переменной ис.

3. Композиция (аккумуляция): объединение найденных усеченных функций, в результате получаем итоговое нечеткое множество для переменной выхода с результирующей функцией принадлежности:

дс(и) = ^1(и) V д2(и),

где через "V" обозначена функция логического максимума.

4. Приведение к четкости (дефаззификация): нахождение четкого значения выходной переменной и2 центроидным методом.

Рассмотрим настройку нечёткого регулятора в контуре управления координатой х:

1. Указываются лингвистические термы, в данном случае входные переменные Ах и АУХ и выходная переменная: тангаж.

д1(и) = А - ¡и1(ис ); .у.2(и) = В - р2(ис )

Рисунок 12 - Структурная схема нечеткого регулятора, показывающая связь

входа и выхода

2. Определяются диапазоны значений переменных. Диапазон выставляется таким образом, чтобы он охватывал все возможные значения для переменной.

Рисунок 13 - Функции принадлежности входа Ах для базы правил нечеткого

регулятора

3. Определяются термы множества для лингвистических переменных. Для Ах используем следующие три термы: Минус = {отрицательная ошибка}, Ровно = {нулевая ошибка}, Плюс = {положительная ошибка}. Термов может быть больше: чем больше лингвистических переменных, тем точнее настройка.

4. Выбираются функции принадлежности для термов. Функции принадлежности (ФП) задаем в виде кусочно-линейных функций. Для 2 термов данные функции задаем в трапециевидной форме, а для центральной в треугольной. Можно использовать и другие фигуры (7, Б-образные функции и т.д.), но треугольники и трапеции использовать удобней. ФП располагаются так, чтобы они перекрывали друг друга. Чем больше перекрытие, тем выше степень принадлежности. Мы можем настраивать ширину для понятия «ровно», которая влияет на заброс.

5. Создается база правил нечеткого вывода:

ПРАВИЛО 1: ЕСЛИ «Отрицательная ошибка» И «Отрицательная ошибка», ТО «Пикирование».

ПРАВИЛО 2: ЕСЛИ «Отрицательная ошибка» И «Нулевая ошибка», ТО «Пикирование»

ПРАВИЛО 3: ЕСЛИ «Отрицательная ошибка» И «Положительная ошибка», ТО «Ровно»

И т.д.

Рисунок 14 - Редактор базы правил

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

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

1. Линеаризовать систему и анализировать устойчивость с помощью критериев Найквиста или Михайлова.

2. Упростить систему, например, для обратной связи по скорости ¥х.

Для скорости рассмотреть устойчивость в линейном приближении (вблизи

нулевого значения), а для оставшейся части НЛ-регулятора (управление по отклонению координаты х) рассмотреть критерий абсолютной устойчивости (Попова) для определения влияния нелинейности.

3. Построить и проанализировать фазовые траектории по критерию устойчивости Ляпунова.

Способ 1: Для анализа устойчивости линеаризуем систему, рассчитав коэффициенты обратной связи по координате К^ и скорости К^ вблизи нулевых значений этих переменных по формулам:

3* д0

К = __

Рх 2* Ах0 3 * д0

Ках = о „м

х 2 * V,

0

где Дх0, УХо, д0 - полуширина основания функций принадлежности треугольной формы в описания правил принятия решения НЛ-регулятора.

Как и в предыдущих примерах, для определения запаса по фазе и амплитуде воспользуемся логарифмической амплитудно-фазовой частотной характеристикой (графиком Боде) для разомкнутой системы).

50

Рисунок 15 - ЛАФЧХ при настройке контура с помощью НЛ

Результаты анализа ЛАФЧХ показывает, что система управления с использованием предлагаемого алгоритма устойчива.

Способ 2:

Упростим систему (канал управления координатой х), линеаризовав обратную связь по скорости Vx. В нечетком регуляторе остается только управление отклонениям Лх. Определим влияние нелинейности нечеткого контроллера с помощью критерия абсолютной устойчивости (критерия Попова) (рис. 16):

Рисунок 16. Упрощенная структура канала управления координатой х

Определим влияние нелинейности нечеткого контроллера с помощью критерия абсолютной устойчивости (критерия Попова).

д.

зад

К

Рисунок 16 - Управляющая поверхность системы нечеткого логического вывода для формирования заданного значения тангажа по отклонениям координаты

и скорости от заданных значений

1. Строим функцию принадлежности и составляем уравнение прямой,

описывающее предельный темп роста нелинейной зависимости. В нашем случае это

уравнение прямой: $ = к * Ах. Определяем коэффициент наклона прямой к.

37

2. Определяем точку, через которую должна пройти прямая с заданным 1

наклоном - -, не пересекающая статическую характеристику системы управления. к

3. Как видно из рис. 17, статическая характеристика полностью расположена правее прямой, проходящей через точку (в рассмотренно примере -точку -4.167; /0), с некоторым угловым коэффициентом, следовательно, согласно критерию Попова данная система является абсолютно устойчивой.

Рисунок 17 - Критерий абсолютной устойчивости

Способ 3: В зависимости от начальных условий Vx и X был построен фазовый портрет системы и определена устойчивость по методу Ляпунова

Х=Ух

дзаА = ЫЬ(Х,Ух) тд =

1

т2

1-г

Ух

——Г

й

Г Й

Я

"1-г и в ±р ФР 1\ —1—к П—т~ -1—\-

X

—г~т

——г

1—г

——К

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

~гг

Рисунок 18 - Фазовый портрет системы с нечетким регулятором

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

Другими словами, система устойчива после начальных отклонений, следовательно, она устойчива в целом.

8. Результаты моделирования с использованием предлагаемых алгоритмов

На основе описанных алгоритмов настройки параметров САУ было проведено моделирование контура управления координатой X.

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

Рисунок 19 - Настройка контура управления координатой X с использованием ЛКР, нечёткого регулятора и метода обратной задачи динамики

В случае наличия постоянного внешнего возмущения система, настроенная с

помощью обратной задачи динамики, удовлетворяет требованиям к

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

Рисунок 20 - Настройка контура управления координатой X с использованием различных регуляторов с учётом влияния постоянного внешнего возмущения

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

Рисунок 21 - Настройка контура управления координатой X с использованием различных регуляторов в случае неточного вычисления эффективностей

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

О 5 10 сек 15

Рисунок 22 - Настройка контура управления координатой X с использованием ПИД-регулятора и метода обратной задачи динамики

X, м

1 ,

:

1'

Хгай эталон

ОЗД-рвгулятор Р! Р-регулятор

15

Ъ сек

Рисунок 23 - Настройка контура управления координатой X с использованием ПИД-регулятора и метода обратной задачи динамики при неточных значениях эффективностей и при влиянии постоянного внешнего возмущения

9. Выводы

Основываясь на результатах вышеописанных исследований, можно сделать вывод о том, что каждый из рассмотренных методов настройки параметров САУ имеет различные преимущества и недостатки, которые можно свести в Таблицу 1.

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

Таблица 1 - Преимущества и недостатки алгоритмов настройки параметров

САУ

Алгорит м настройки параметров САУ Просто та использования Малая чувствительность по отношению к внешним возмущениям Возможность настройки: Запасы устойчивости

В частотной области Во временной области о фазе По амплитуде

ПИД-регулятор + + + 8° 15, 3 дБ

ЛКР + 1,7° 24, 8 дБ

НЛ + 6,1° 19, 9 дБ

ОЗД + + + 1,8° 23 дБ

Список источников

1. А. Брайсон, Хо Ю-Ши. Прикладная теория оптимального управления. - М.: Мир, 1972. - 544 с.

2. Алхаддад Мухаммад. Моделирование и управление ориентацией квадрокоптера с использованием линейного квадратического регулятора // Актуальные проблемы авиации и космонавтики. 2016. Т.1. № 12. С. 883-886.

3. Аполлонов Д.В., Бибикова К.И., Гаврилова А.В., Шибаев В.М. Вопросы выбора архитектуры системы автоматического управления преобразуемого беспилотного летательного аппарата - «конвертоплана» // Труды МАИ. 2021. №121. URL: http://trudymai.ru/published.php?ID=162672. DOI: 10.34759/trd-2021-121-25

4. Верещиков Д.В., Волошин В.А., Ивашков С.С., Васильев Д.В. Применение нечеткой логики для создания имитационной модели управляющих действий летчика // Труды МАИ. 2018. № 99. URL: http://trudymai.ru/published.php?ID=91926

5. Виноградов С.С. Синтез нечёткого навигационного регулятора для малоразмерного вертолёта «Раптор» // Труды МАИ. 2012. № 73. URL: http://trudymai.ru/published.php?ID=48562

6. Дмитриев В.И., Звонарев В.В., Лисицын Ю.Е. Методика обоснования рациональных способов управления беспилотным летательным аппаратом // Труды МАИ. 2020. № 112. URL: http://trudymai.ru/published.php?ID=116566. DOI: 10.34759/trd-2020-112-16

7. Кожевников В.А. Автоматическая стабилизация вертолетов. - М.: Машиностроение, 1977. - 152 с.

8. Крутько П.Д. Обратные задачи динамики в теории автоматического управления. - М.: Машиностроение, 2004. - 576 с.

9. Пегат А. Нечеткое моделирование и управление. - М.: Лаборатория знаний, 2020. - 801 с.

10. Попов Е.П. Теория линейных систем автоматического регулирования и управления. - М.: Наука. Главная редакция физ.-мат. литературы, 1989. - 304 с.

11. Сторожев С.А., Хижняков Ю.Н. Новый метод адаптации регулятора состояний с применением нечеткой логики // Труды МАИ. 2021. № 118. URL: http://trudymai.ru/published.php?ID=158255. DOI: 10.34759/trd-2021-118-16

12. Тищенко М.Н., Артамонов Б.Л. Проблемы повышения крейсерской скорости полета вертолета и пути их решения // Труды МАИ. 2012. № 55. URL: http ://trudymai.ru/published.php?ID=30012

13. David J. Murray-Smith. A Review Of Inverse Simulation Methods And Their application // International Journal of Modelling and Simulation, 2014, no. 34(3), 120-125. DOI:10.2316/Journal.205.2014.3.205-5906

14. Gareth D. Padfield. Helicopter Flight Dynamics Including a Treatment of Tiltrotor Aircraft, Wiley, 2018, 856 p.

15. Jingxian Liao. Mathematical modelling and model predictive controller design of a quad tiltrotor UAV, 2020. D0I:10.1177/0954406220971330

16. Ke Lu. Flight Dynamics Modeling and Dynamic Stability Analysis of Tilt-Rotor Aircraft // International Journal of Aerospace Engineering, 2019, no. 2, pp. 1-15. DOI: 10.1155/2019/5737212

17. Kristi Marie Kleinhesselink. Stability and Control Modeling of Tiltrotor Aircraft, 2007.

18. Luis Rodolfo Garcia Carrillo et al. Quad Rotorcraft Control. Vision-Based Hovering and Navigation, 2013.

19. Saied M. BFA fuzzy logic based control allocation for fault-tolerant control of multirotor UAVs // Aeronautical Journal -New Series, 2019, no. 123 (1267), pp. 13561373. D0I:10.1017/aer.2019.58

20. Mark E. Dreier. Introduction to Helicopter and Tiltrotor Simulation, American Institute of Aeronautics and Astronautics, 2007, 637 p.

21. Masayuki Sato. Flight Controller Design and Demonstration of Quad-Tilt-Wing Unmanned Aerial Vehicle // Journal of Guidance, Control, and Dynamics, 2014, no. 38(6), pp. 1-12. DOI: 10.2514/1.G000263

22. Olena Matiychyk. Integration Of Tiltrotor Aircraft Into Modern Air Transport Systems // Conference "Science - Future of Lithuania. Transport Engineering and Management", 2017. URL: http://jmk.transportas.vgtu.lt/index.php/tran2017/tran2017/paper/view/106

23. Takahito Mikami. Design of Flight Control System for Quad Tilt-Wing UAV, 2015. D0I:10.1109/ICUAS.2015.7152364

24. Thomas Lombaerts. Nonlinear Dynamic Inversion based Attitude Control for a Hovering Quad Tiltrotor eVTOL Vehicle // AIAA Scitech 2019 Forum, 2019. D0I:10.2514/6.2019-0134

25. Xue Ping Zhu., Yong Hua Fan. Design of Tiltrotor Flight Control System Using Fuzzy Sliding Mode Control, 2010. D01:10.1109/ICMTMA.2010.211

References

1. A. Braison, Kho Yu-Shi. Prikladnaya teoriya optimal'nogo upravleniya (Applied optimal control), Moscow, Mir, 1972, 544 p.

2. Alkhaddad Mukhammad. Aktual'nye problemy aviatsii i kosmonavtiki, 2016, vol. 1, no. 12, pp. 883-886.

3. Apollonov D.V., Bibikova K.I., Gavrilova A.V., Shibaev V.M. Trudy MAI, 2021, no. 121. URL: http://trudymai.ru/eng/published.php?ID=162672. DOI: 10.34759/trd-2021-121-25

4. Vereshchikov D.V., Voloshin V.A., Ivashkov S.S., Vasil'ev D.V. Trudy MAI, 2018, no. 99. URL: http ://trudymai.ru/eng/published.php?ID=91926

5. Vinogradov S.S. Trudy MAI, 2012, no. 73. URL: http://trudymai.ru/eng/published.php?ID=48562

6. Dmitriev V.I., Zvonarev V.V., Lisitsyn Yu.E. Trudy MAI, 2020, no. 112. URL: http://trudymai.ru/eng/published.php?ID=116566. DOI: 10.34759/trd-2020-112-16

7. Kozhevnikov V.A. Avtomaticheskaya stabilizatsiya vertoletov (Automatic stabilization of helicopters), Moscow, Mashinostroenie, 1977, 152 p.

8. Krut'ko P.D. Obratnye zadachi dinamiki v teorii avtomaticheskogo upravleniya (Inverse problems of dynamics in the theory of automatic control), Moscow, Mashinostroenie, 2004, 576 p.

9. Pegat A. Nechetkoe modelirovanie i upravlenie (Fuzzy modeling and control), Moscow, Laboratoriya znanii, 2020, 801 p.

10. Popov E.P. Teoriya lineinykh sistem avtomaticheskogo regulirovaniya i upravleniya (Theory of linear automatic control and control systems), Moscow, Nauka. Glavnaya redaktsiya fiz.-mat. literatury, 1989, 304 p.

11. Storozhev S.A., Khizhnyakov Yu.N. Trudy MAI, 2021, no. 118. URL: http://trudymai.ru/eng/published.php?ID=158255. DOI: 10.34759/trd-2021-118-16

12. Tishchenko M.N., Artamonov B.L. Trudy MAI, 2012, no. 55. URL: http ://trudymai.ru/eng/published.php?ID=3 0012

13. David J. Murray-Smith. A Review Of Inverse Simulation Methods And Their application, International Journal of Modelling and Simulation, 2014, no. 34(3), 120-125. DQI:10.2316/Journal.205.2014.3.205-5906

14. Gareth D. Padfield. Helicopter Flight Dynamics Including a Treatment of Tiltrotor Aircraft, Wiley, 2018, 856 p.

15. Jingxian Liao. Mathematical modelling and model predictive controller design of a quadtiltrotor UAV, 2020. DOI:10.1177/0954406220971330

16. Ke Lu. Flight Dynamics Modeling and Dynamic Stability Analysis of Tilt-Rotor Aircraft, International Journal of Aerospace Engineering, 2019, no. 2, pp. 1-15. DOI: 10.1155/2019/5737212

17. Kristi Marie Kleinhesselink. Stability and Control Modeling of Tiltrotor Aircraft, 2007.

18. Luis Rodolfo Garcia Carrillo et al. Quad Rotorcraft Control. Vision-Based Hovering and Navigation, 2013.

19. Saied M. BFA fuzzy logic based control allocation for fault-tolerant control of multirotor UAVs, Aeronautical Journal -New Series, 2019, no. 123 (1267), pp. 13561373. D0I:10.1017/aer.2019.58

20. Mark E. Dreier. Introduction to Helicopter and Tiltrotor Simulation, American Institute of Aeronautics and Astronautics, 2007, 637 p.

21. Masayuki Sato. Flight Controller Design and Demonstration of Quad-Tilt-Wing Unmanned Aerial Vehicle, Journal of Guidance, Control, and Dynamics, 2014, no. 38(6), pp. 1-12. DOI: 10.2514/1.G000263

22. Olena Matiychyk. Integration Of Tiltrotor Aircraft Into Modern Air Transport Systems, Conference "Science - Future of Lithuania. Transport Engineering and Management", 2017. URL: http://imk.transportas.vgtu.lt/index.php/tran2017/tran2017/paper/view/106

23. Takahito Mikami. Design of Flight Control System for Quad Tilt-Wing UAV, 2015. D0I:10.1109/ICUAS.2015.7152364

24. Thomas Lombaerts. Nonlinear Dynamic Inversion based Attitude Control for a Hovering Quad Tiltrotor eVTOL Vehicle, AIAA Scitech 2019 Forum, 2019. D0I:10.2514/6.2019-0134

25. Xue Ping Zhu., Yong Hua Fan. Design of Tiltrotor Flight Control System Using Fuzzy Sliding Mode Control, 2010. DOI:10.1109/ICMTMA.2010.211

Статья поступила в редакцию 22.11.2021; одобрена после рецензирования 15.12.2021; принята к публикации 21.02.2022.

The article was submitted on 22.11.2021; approved after reviewing on 15.12.2021; accepted for publication on 21.02.2022.

i Надоели баннеры? Вы всегда можете отключить рекламу.