Расчет границы области асимптотической устойчивости динамической системы Текст научной статьи по специальности «Физика»

УЧЕНЫЕ ЗАПИСКИ ЦАГИ Т ом ХХ/ 19 90

мз

УДК 629.735.33.015.017.22/23

РАСЧЕТ ГРАНИЦЫ ОБЛАСТИ АСИМПТОТИЧЕСКОЙ УСТОЙЧИВОСТИ ДИНАМИЧЕСКОЙ СИСТЕМЫ

М. Г. Гоман, А. В. Храмцовский

Предложен алгоритм расчета двумерных сечений границы области асимпт0ти4еск0й уст0йчив0сти установившегоСя решения многомерной ди-наМИЧеСк0й системы. Рассматривается СЛучай, когда граница образована сепаратрисной поверхностью, проходящей через седловое стационарное решение. Алгоритм включает в себя решение краевой задачи посредством Численного интегрирования динамических уравнений и использования метода непрерывного продолжения. В качестве примера рассмотрена задача бокового движения самолета, которое описывается нелинейными уравнениями третьего порядка. Проведено качественное исследование особых решений и расчет сечений области асимптотической устойчивости при изменении параметра управления.

Можно выделить два основных подхода к задаче об определена областей асимптотической устойчивости нелинейной динамической системы. Широко используемым инструментом аналитического исследования устойчивости «в большом» являются функции Ляпунова и прямой метод Ляпунова, которым посвящена обширная литература [1, 2]-Имеется большое число частных способов построения функций Ляпунова, в том числе разработанные в последнее время машинные методы [3]. Получаемые с их помощью оценки области притяжения зависят от вида выбранной функции Ляпунова и не всегда достаточно хорошо совпадают с искомой областью.

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

Вторая группа методов основывается на качественных представлениях о структуре особых решений динамической системы и на численном восстановлении границы области асимптотической устойчивости [4, 5].

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

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

Х=Г(Х), 0)

где Х е Яп; F — вектор-функция, задающая отображение Яп-+Яп 11 удовлетворяющая известным достаточным условиям существования и единственности решения Х ^, Хо) при начальном условии Х (О, Х0) =Хо.

Пусть начало координат является изолированным и асимптотически устойчивым состоянием равновесия системы (1). При этом система уравнений (1) может иметь и другие особые точки (как устойчивые, так и неустойчивые), а также более сложные установившиеся решения.

Область асимптотической устойчивости, или область притяжения й точки Х= О определяется как множество всех точек Хо, таких, что: Нт Ху, Хо) =0. Построение области й является одной из центральных

проблем в теории устойчивости. Граница Ев области асимптотической устойчивости й при определенных условиях является гладкой и плотно заполненной целыми фазовыми траекториями, при этом ее структура может быть простой и доступной для приближенного вычисления [6].

У систем второго порядка (п=2) область асимптотической устойчивости ограничивается либо неустойчивой замкнутой траекторией (неустойчивым циклом), либо траекториями, входящими в седловые особые точки, либо замкнутой кривой критических точек. В системах большей размерности (п:;3) возможные ситуации становятся значительно разнообразнее и сложнее [6].

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

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

Аналогичным способом можно строить сепаратрисную поверхность и для систем более высокого порядка, в частности при п=3 и 4 [5]. Однако метод построения сепаратрисной поверхности путем интегрирования в обратном времени требует большого объема вычислений, для того чтобы достаточно плотно заполнить фазовыми траекториями поверхность, особенно при большом удалении от седловой особой точки.

Ограничимся случаем, когда граница области асимптотической устойчивости установившегося решения многомерной динамической системы полностью или частично образована сепаратрисной поверхностью 5;;-1, проходящей через седловую особую точку типа О”-1’1 (обозначения особых решений взяты из работы [6]).

Размерность сепаратрисной поверхности на единицу меньше,

чем размерность полного пространства. Сепаратрисная поверхность состоит из фазовых траекторий, которые при 1-+ 00 стремятся в седловую особую точку Оп-м, имеющую одно положительное собственное число. Из седловой особой точки выходят две полутраектории , одна из которых захватывается исследуемым устойчивым решением Оп0, а вторая — каким-либо другим устойчивым решением (рис. 1).

Восстановлению таких сепаратрисных поверхностей должен предшествовать достаточно полный качественный анализ всех особых решений динамической системы (1). Требуется, в частности, найти и провести локальное исследование всех седловых особых точек типа °п“1-1. Расчет траекторий выходящих из этих седловых точек, помогает

Рис. 1 Рис. 2

упорядочить устойчивые установившиеся решения, такие, как стационарные °п’ о или периодические Г”’1, области устойчивости которых разделяются сепаратрисными поверхностями типа 5 :-1..

2. Алгоритм расчета двумерных сечений сепаратрисных поверхностей. Наглядное представление о пространственной области устойчивости можно получить по набору двумерных сечений сепаратрисной поверхности 5;-1 различными плоскостями Рг- Линии Sf, ограничивающие сечение области устойчивости, являются пересечениями выбираемых двумерных плоскостей и сепаратрисной поверхности: 5? = Р2 п (рис. 2).

Фазовые траектории, проходящие через точки линии 5^, принадлежат и при ^ 00 стремятся к особой точке О”-1’1. Смещение фа-

зовой точки в плоскости Рг с линии 5^ приводит к тому, что выходящая из нее фазовая траектория уже не попадет в особую точку О”-1’1.

В плоскости Р2 введем координаты (|, '1']).

Далее, введем меру «промаха» для фазовой траектории, проходящей вблизи границы 5^, относительно особой точки О”-1’1.

Для этого рассмотрим гиперплоскость £”-ь проходящую через особую точку О”-1’1, и касательную к сепаратрисной поверхности Эта

гиперплоскость представляет собой линейное подпространство размерности (п—1) с базисом, образованным собственными векторами матрицы Якоби динамической системы в особой точке, которые соответствуют (п—1) собственным числам с отрицательными действительными частями. Далее, зададим цилиндрическую поверхность С я, являющуюся геометрическим местом точек, равноудаленных на расстояние Я от прямой, проходящей через особую точку О”-1’1 и ортогональной к гиперплоскости £”-1 (см. рис. 2). Направление оси цилиндра Ся совпадает, таким образом, с направлением вектора, ортогонального тем собственным векторам, на которые «натянута» гиперплоскость £л-ь В ряде случаев ось цилиндра может совпадать с направлением выхода из особой точки траекторий

В окрестности линии 5^ на плоскости Р2 существует некоторая область точек, такая, что начинающиеся из нее фазовые траектории пересекают поверхность цилиндра Ся- Размер этой области зависит от радиуса цилиндра Я. Мерой «промаха» для фазовой траектории может служить расстояние от фазовой точки до гиперплоскости £п-1 в момент пересечения фазовой траекторией цилиндра Ся. Для вычисления меры промаха будем использовать скалярную функцию Л:

Л (Х) = (Х -Хо)т -По,

где вектор Хо — задает координаты особой точки О”-м, а по — единичный вектор нормали к гиперплоскости £”-1.

6 -- Ученые записки N2 3

81

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

г(Х) = р(Х — Хо-й(Х). по) - Л, (2)

где р (г) — евклидова норма вектора г.

Отображение точки плоскости Рг на поверхность цилиндра Ся, задаваемое условием г(Х) =0, позволяет поставить в соответствие каждой точке из некоторой окрестности линии 5^ на плоскости Рг скалярную величину Ля- Это отображение и, соответственно, определение функции Ля (|, 11) выполняется посредством численного интегрирования системы уравнений (1), с одновременным вычислением и контролем знака функции (2).

Задача отыскания линии 5^ — сечения сепаратрисной поверхности плоскостью Рг — сводится к нахождению однопараметрического семейства решений скалярного уравнения, зависящего от двух переменных | и 11 — координат на плоскости Рг начальной точки отображения:

АЛ(6, Ю = 0. (3)

Для того чтобы алгоритм работал эффективно, важно правильно выбрать размер цилиндра Я. С одной стороны, этот размер должен быть достаточно мал, для того чтобы касательная гиперплоскость Ьп-1 хорошо приближала сепаратрисную поверхность внутри и на поверхности цилиндра; а с другой, — достаточно велик, чтобы размер окрестности линии 5^, из которой выходящие фазовые траектории пересекают цилиндр, был не слишком мал.

Условие (3) определяет приближенную границу области устойчивости, которая при Я;!- 0 стремится к сечению 5^ сепаратрисной поверхности 5+1-

Построение однопараметрического семейства решений уравнения

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

(4)

Требуемая точность решения обеспечивается соответствующим выбором коэффициента k и шага интегрирования М в (4). Величины и А’ — частные производные функции Ля по переменным | и 11 — определяются численно.

3. Качественное исследование нелинейной динамики бокового движения самолета. При изучении особенностей пространственного движения самолета часто приходится решать нелинейные динамические задачи. Большую помощь в этом оказывают качественные методы [8].

Уравнения пространственного движения самолета можно существенно упростить, если собственная частота в продольном движении

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

/) Юу ю2 +

/у шу = ^х “х Ю + Муа,

Р == (В соэ а + Юх 51па

(5)

где проекция угловой скорости на боковую ось определяется из кинематической схемы движения: мг=ЮхР (угол скольжения р считается малой величиной); /х, /у, т — главные моменты инерции и масса самоле-

та; V — скорость полета.

Аэродинамические моменты Мха,

Муа и сила могут быть пред-

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

®х = а1 Р + 02^ + аз Юу + Мх,

шу=(Ьі + *2 Шж)Р + Ь3Шу + Ь4Шд., Р — Фу ^ а0 + Юх sln а0 + СіР"

(6)

Параметры, входящие в уравнения (6), зависят от режима полета, аэродинамических и инерционных характеристик самолета [8, 9]. В настоящей работе при проведении расчетов эти параметры принимались равными:

01 = -14с-2, а2 = — 1,4 с-1, а3=О, &1 = -2с~2,

Ь2=0,6, Ь3 = — 0,2 с-1, Ь( = О, ао = 0,1, с, = 0.

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

[Мх] = с-2, [«х, Фу] = с-1.

Полагая правые части уравнений (6) равными нулю и решая полученную систему нелинейных уравнений, можно определить координаты особых точек системы уравнений (6). При всех значениях параметра управления Мх система (6) имеет три особые точки. На рис. 3 приведена зависимость одной из координат Ро всех трех особых точек от параметра управления. Характер устойчивости возмущенного движения и структура особых точек могут быть определены по собственным числам линеаризованной системы уравнений (6).

Две особые точки — О|л и О|л с ненулевыми значениями координат при Мх = О при всех значениях параметра управления неустойчивы. Эти особые точки имеют седло-фокусную структуру, что определяется расположением собственных чисел на комплексной плоскости — одно

положительное действительное число и пара комплексно-сопряженных с отрицательной действительной частью.

Ветвь решений, проходящая через начало координат, соответствует устойчивым особым точкам О^0 при Мж<3,55. При больших значениях

параметра управления А1ж>3,55 возникает неустойчивость колебательного характера О®'0 -+ 0\л + Г3’1— комплексно-сопряженная пара собственных чисел переходит в правую полуплоскость. Особая точка становится неустойчивым пространственным фокусом.

В момент изменения характера устойчивости особой точки при •/Иж~3,55 появляются устойчивые замкнутые траектории Г3’1 — предельные циклы, амплитуда которых быстро растет с увеличением параметра Мх. На рис. 3 изображены амплитуды возникающих колебаний по. углу скольжения. Эти колебания существуют в узком диапазоне значений параметра управления.

Форма цикла изменяется по мере роста его амплитуды. На рис. 4 приведены проекции замкнутых фазовых траекторий на плоскость фазовых координат Юх, юу при нескольких значениях параметра управления Мх. Наряду со значением самого параметра, указаны велечины периода колебаний Т и значения одного из мультипликаторов — собственных чисел якобиана точечного отображения в неподвижной точке. Расчет периодических решений, представленных на рис. 3, 4 сделан с помощью алгоритма поиска неподвижной точки точечного отображения, определенного на некоторой выбранной секущей плоскости [10]. С ростом амплитуды период колебаний увеличивается, особенно значительно — при приближении цикла к неустойчивой особой точке О^1. На рис. 4 видна деформация цикла и вытягивание его в направлении неустойчивой особой точки. Изменение положения особых точек при этом незначительно, в частности, на рис. 4 они изображены при различных значениях параметра управления: внутренняя особая точка О^о при Мж = 3,55 (в момент возникновения замкнутых траекторий), а седловая особая точка 0%1 соответствует моменту слияния с ней замкнутой траектории

при Мх",3,9.

При значениях параметра управления Мж>3,9 замкнутые фазовые траектории исчезают, и уравнения (6) уже не имеют никаких устойчивых установившихся решений, поскольку все три особые точки при этом неустойчивы.

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

(6), которые существуют при данном значении параметра управления. При Мх<3,55 установление параметров движения носит практически апериодический характер. При приближении к Мх-+■ 3,55 возникают плохо затухающие колебания, а при 3,55<Мх<3,9 переходный процесс приобретает установившийся колебательный характер. При превышении критической величины Мх;;:3,9 параметры движения после нескольких колебаний резко и необратимо возрастают [9].

Существование устойчивых решений — стационарных при Мх<3,55 и периодических при 3,55<Мх<3,9 — еще не гарантирует того, что они будут реализованы при управляемом движении. При каждом значении параметра управления имеется определенная область начальных значений фазовых переменных, из которой происходит выход на устойчивый установившийся режим движения. Определение области устойчивости в данной задаче связано с расчетом положения сепаратрисных поверхностей 5+.

В рассматриваемой задаче область притяжения устойчивой особой точки О®0, или устойчивой замкнутой траектории Г3’1 ограничивается двумя сепаратрисными поверхностями, проходящими через неустойчивые особые точки О^1 и О}1. Такая картина является наиболее характерным случаем и при рассмотрении других задач пространственного движения самолета большей размерности [8], в которых сепаратрис-ные поверхности проходят через особые точки типа 0я-1'1.

На рис. 5, а, б приведены различные сечения двух сепаратрисных поверхностей уравнений (6) при нулевом значении параметра управления Мх=О. Результаты получены с помощью описанного выше алгоритма. Плоскости сечений Р2 в этих случаях выбирались перпендикулярными оси Шх, в случае (а) плоскость сечения Р2 проходит через особую точку О®'0 в начале координат, в случае (б) —через седловую особую точку Ор. Области плоскости Р, из которых фазовые траектории захватываются устойчивой особой точкой О3-0, выделены штриховкой.

На рис. 5, б в проекции на выбранную плоскость сечения изображены характерные фазовые траектории, выходящие из особой точки 0«

Изменения сечений се-паратрисных поверхностей и соответствующая деформация областей захвата особой точки О3’0 при вариации значений параметра управления Мх приведены на рис. 6, а,

.6, в. Плоскость сечения при этом выбрана другим образом. Она фиксирована и определяется условием Шу = О.

Уже при Мх=1,5 (см. рис. 6, б) сечение области «захвата» становится двусвязным, что происходит из-за пространственной деформации и подгибания се-

Шх,=0; “Л2=-7,9/; ш,л- 1,91

Ллоскость сечения Р,

Плоскость сечения Р,

проходит через начало проходит через сеИлобт координат (0?'°ъРг) точку (0}’'е ^2)

о)

б)

Рис. 5

Мх=О

шу,—0,09е шу’г=0.189

Шу=-0,202

Мх=2,б

Ыу2=0,188

М,х=2,925

шу,=-';.155

Ыу. -0,187

сечение области захвата устойчивой особой точки о]'0

Рис. 6

паратрисной поверхности, проходящей через неустойчивую особую точку Оз'1 , к которой приближается устоичивая особая точка 0\ ,0-При дальнейшем увеличении параметра Мж=2,6 (см. рис. 6, в) появляется еще одна внутренняя область «неустойчивости» при |3<0, ю*<0, а при Мх = 2,925 (см. рис. 6, г) происходит слияние внутренних областей, «неустойчивости», в результате чего образуются узкие полоски «захвата».

Если рождение периодических решений — замкнутых фазовых траекторий у уравнений (6), семейство проекций которых приведено на рис. 4, связано с бифуркацией особой точки О3,о -+ О}’2 + Г3-1 при

Мх::3,55, то исчезновение замкнутой фазовой траектории связано с бифуркацией сепаратрисной поверхности 5^ и особой траектории 5^,, связанных с особой точкой 0|л. В момент исчезновения предельного цикла при Мх::3,9 особая траектория 5^ становится двояко асимтотической к особой точке О24 и совпадает с предельным циклом.

На рис. 7, а, б приведены сечения сепаратрисной поверхности плоскостью, проходящей через особую точку 0|’1 ,перпендикулярно оси Юу, На рис, 7, а сечение сепаратрисной поверхности при Мх = 3,875 охватывает небольшую область, из которой траектории захватываются предельным циклом Г3’1,

Превышение критического значения параметра М*::4,1 приводит к перестроению поверхности 5^, что видно по рассчитанному ее сечению, — она начинает спирально накручиваться на фазовую траекторию 5]+> асимптотически входящую в особую точку 0{’2.

Таким образом, поверхность S+(O<1’i) , являющаяся многообразием

траекторий, стремящихся при t-+ 00 в особую точку Оз'1, уже перестает быть границей области устойчивости, при этом, как уже отмечалось, исчезают устойчивые установившиеся решения.

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

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

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

ЛИТЕРАТУРА

1. Б а р б а ш и н Е. А. Функция Ляпунова. — М.: Наука, 1974.

2. К и р и н Н. Е., Н е л е п и н Р. А., Б а й д а е в В. Н. Построение области притяжения по методу Зубова. — Дифференциальные уравнения, 1981, т. 17, №8.

3. Метод функций Ляпунова в динамике нелинейных систем. /Под ред. В. М. Матросова, Р. И. Козлова. — Новосибирск.: Наука, Сиб. отд.,

1983.

4. G е n е s i о Р., Т а r t а g 1 i а М., V i с i n о А. Оп the estimation of asymptotic stability regions: State of the art and new proposals, IEEE Trans. Оп aut. contr., V. АС — 30, 1985, N 8.

5. Т о н д л А. Нелинейные колебания механических систем. — М.:

Мир, 1973.

6. Б у т е н и н Н. В., Н е й м а р к Ю. И., Фу ф а е в Н. А. Введение в теорию нелинейных колебаний. — М.: Наука, 1976.

7. Г о м а н М. Г. Дифференциальный метод продолжения решения , систем конечных уравнений, зависящих от параметров. — Ученые записки ЦАГИ, 1986, т. 17, .N2 5.

8. Б ю ш г е н с Г. С., С т у д н е в Р. В. Динамика самолета. Пр ост-ранственное движение. — М.: Машиностроение, 1983.

9. Г о м а н М. Г., С у х а н о в В. Л. Автоколебательные режимы движения самолета при вращении по крену. — Ученые записки ЦАГИ,

1985, т. 16, .N2 2.

10. Z а g а у n о v G. 1., G о та n М. G. Bifurcation analysis of cri-Нса1 aircraft flight regimes. — ICAS Proceedings, Toulouse, 1984.

Рукопись поступила 24/JV 1989 г.