Построение кусочно-линейной функции Ляпунова для динамических систем второго порядка Текст научной статьи по специальности «Математика»

Математика к Математическое

моделирование

Сетевое научное издание УДК 517.938

Ссылка на статью:

// Математика и математическое моделирование. МГТУ им. Н.Э. Баумана. Электрон. журн. 2017. №3. С. 32-43.

Б01: 10.24108/шаШш.0317.0000073

Представлена в редакцию: 13.06.2017 © МГТУ им. Н.Э. Баумана

Построение кусочно-линейной функции Ляпунова для динамических систем второго порядка

Касаткина Т. С.1'*

kasatkina_t_s@mail.ru 1МГТУ им. Н.Э. Баумана, Москва, Россия

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

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

Введение

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

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

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

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

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

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

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

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

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

1. Постановка задачи

Рассмотрим автономную систему дифференциальных уравнений

х = / (х), (1)

где х Е Ега — вектор состояния; /: Б ^ Ега, Б С Ега, — непрерывно-дифференцируемое отображение. Считаем, что точка х = 0 — положение равновесия рассматриваемой системы, т.е. /(0) = 0.

Функция V: D ^ R называется положительно определенной, если V(0) = 0 и V(x) > 0 при x = 0. Положительно определенную функцию V называют функцией Ляпунова, если V(x) < 0, x G D, где V — производная функции V в силу системы (1).

Прямой метод Ляпунова исследования положения равновесия на устойчивость основан на следующей теореме [9].

Теорема 1. Пусть x = 0 — положение равновесия системы (1), область D содержит точку x = 0, V — непрерывно дифференцируемая положительно определенная в D функция. Тогда:

- если V(x) < 0 в D \ {0}, то положение равновесия x = 0 устойчиво;

- если V(x) < 0 в D \ {0}, то положение равновесия x = 0 асимптотически устойчиво;

-если V(x) > Ci||x|| и V(x) < —C2||x|| в D, где ||x|| — какая-либо норма в Rn, а Сь

C2, а — положительные постоянные, то положение равновесия x = 0 экспоненциально устойчиво.

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

Величину

D+V(xo) = lim V(x0 +hf(xo) h^+o h

называют верхней производной Дини функции V в точке x0 по направлению векторного поля f, соответствующего системе (1). Оказывается, что утверждение теоремы 1 остается в силе, если в формулировке производную в силу системы заменить верхней производной Дини [4, 8].

Мы ограничимся рассмотрением случая положения равновесия двумерной системы, которое является экспоненциально устойчивым. В этом случае функция Ляпунова V(x) должна удовлетворять условиям [6]

V(x) > Ci11x112, D+V(x) < -C2||x||2. (2)

Функцию Ляпунова будем строить на основе некоторой триангуляции, покрывающей некоторую область X, близкую к исходной области D: достаточно выбрать некоторую многоугольную окрестность точки x = 0, целиком расположенную в D. Тогда X представляют собой объединение конечного числа треугольников, которые либо не пересекаются, либо стыкуются только по границе. Определив значения функции V в вершинах треугольников, внутри треугольников задаем значения этой функции путем линейной интерполяции.

2. Алгоритм построения функции Ляпунова

Любой треугольник Т определяется тремя своими вершинами х0, X, х2 и является выпуклой оболочкой множества {х0, х, х2}. В дальнейшем будем это отражать обозначением Т = со{х0, х^ х2}. Если в вершинах х0, х, х2 треугольника Т заданы значения функции ^0, Уь У2, то линейная функция V(х), определяемая этими значениями, может быть записана в виде

V (х) = шх + а,

где градиент ш = (ш, ш2) и смещение а этой функции могут быть найдены из соотношений

Ш1 ^ х1 х0 2 х0 1 -1

Ш2 = х11 2 1 1 Vl

а ! х11 2 1 1)

(3)

V

(X, X), ^ = 0, 1, 2, — координатное представление точек х0, х1, х2.

Здесь х^-

В процессе построения будем использовать евклидову норму ||х||2. Контролировать условия экспоненциальной устойчивости для треугольника Т можно по значениям функции и ее производной в вершинах треугольника. Действительно, если

V > С1 ||х ||2, ] = 0, 1, 2, то

V(х) = V(«0^0 + «1^1 + «2X2) = «0^ + «1^1 + a2V2 >

> 11Х012 + а11|Ж11|2 + «2^X2!^ > СЦаХ0 + «1X1 + «2X2Ц2 = СЩх^.

Согласно, [6], можно подобрать постоянные Е1 и Е2 так, что при выполнении неравенств

ш/(х0) < ||х01|2, ш/(х1) + Ех(|тх| + |Ш2|) < -Ш/(х2)+ ^(Ы + Ы) < -будет выполняться неравенство

ш/(х) < — |х|2, х € Т.

| X112, 2

(4)

В качестве постоянных Е1 и Е2 можно выбрать величины [6]

Е1 = ||х1 — х0 || 2 ^21 х1 — х0^2 + ||х1 — х0 || 2 )В, Е = 11X 2 — х0||2( ||х1 — х0 || 2 + 2|х1 — х0|^ В,

где

В > шах вир

Ш,г,з хет

д2/

дх^дхс

(5)

(6)

Пусть имеется триангуляция Т = {Т1, Т2, ..., Т,} заданной области, обладающая тем свойством, что точка х = 0 является вершиной в каждом треугольнике Т„, в который эта

точка входит, причем будем считать, что в перечне вершин такого треугольника точка х0 имеет индекс 0 (рис. 1). Вершины треугольника Ту обозначим , х^д, . Выберем значения У,,0, Кд, У/,2 кусочно-линейной функции V(х) в вершинах треугольников Т, г = 1, 2, ..., V, так, чтобы выполнялись следующие условия:

У,0 > ||Х,0112, У,1 > ||Х,1 |12, Уг,2 > ||Х,21Ь; (7)

К,1| < С,1, К,2| < С,2, (8)

^¿,1 /1 (х,о) + ^¿,1 /1 (х,о) < -||х,0 112, (9)

^¿,1 /1 (х,1)+ ^¿,1 /1 (Жг,1)+ E¿д (Сд + С^ ) < -||Хд Ц2, (10)

^¿,1 /1 (х,2)+ ^¿,2/1 (Хд )+ E¿,2 (С*д + С^ ) < 11 Х,2 112. (11)

Здесь w¿,1, ^¿,2 — координаты вектора ^, вычисленные для треугольника T¿ согласно (3) (они являются линейными функциями неизвестных У,7); E¿,1, E¿,2 — коэффициенты, вычисленные для треугольника T¿ согласно формулам (5); Сд, С,2 — коэффициенты, подлежащие определению.

У;

/и дзА --- А1Х х6,1

х3,0 хгК х5,1 /

\ ^5,0 А ^ч \ X

\ Аз,г/

^^ \х5,2/

Рис. 1. Триангуляция окрестности положения равновесия

Неравенства (7)—(11) обеспечивают выполнение неравенств (4), что в свою очередь гарантирует, что кусочно-линейная функция, заданная значениями У,7, будет удовлетворять соотношениям (3), т.е. будет функцией Ляпунова.

Неравенства (7)—(11) являются линейными относительно неизвестных У,7, Сд, С,2, г = 1, 2, ..., V, 3 = 0, 1, 2. При этом надо иметь в виду, что эти переменные повторяются, поскольку некоторые треугольники имеют одинаковые вершины. Исключив повторения

путем переименования уникальных вершин х^-, получим последовательности хк, Ук, к = 1, 2, ..., К. Добавив целевую функцию, например

Н(У11 У,..., Ук, Со,1, Со,2, С\,\, С\,2,..., Ои,1, а,2) = VI + У> + ... + ^,

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

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

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

Для построения стартовой триангуляции 70 каждый квадрат [т, т +1] х [п, п +1] разбиваем на два треугольники диагональю согласно следующему правилу. Если т > 0, п > 0 или т < 0, п < 0, то разбиение осуществляется диагональю из левого нижнего угла в правый верхний, в остальных случаях разбиение осуществляется другой диагональю, идущей из левого верхнего угла в правый нижний. Таким образом, триангуляцию составляют треугольники [гт,п, ¿ш+1,и, ¿Ш+1,и+1} и {^т,п, ^т,п+1, ¿ш+1,и+1}, если т > 0, п > 0 или

т < ° п < ° и треугольники {zm,n, 2ш+1,и, 2-ш,и+1} и {^ш+1,и+1, ¿ш+1,и, гт,и+1}, если

т > 0, п < 0 или т < 0, п > 0 (рис. 2).

4 -

з -2 -1 -х™ о--1 --2 --з -

-4 □_I_I_I_I_I_I_I_!_

-4-3-2-1 0 1 2 3 4

Х1

Рис. 2. Триангуляция 70

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

Рис. 3. Триангуляция T1

треугольника. Каждый треугольник с вершиной ^0 0 разбиваем на два, соединяя вершину £0 0 с серединой противоположной стороны (рис. 3).

Таким образом, каждый треугольник {ж0, ж1, ж2} триангуляции 70 в случае ж0 = 0 порождает два треугольника

3 3

3 xi + x2

4Х0' 4Хь 4 2

а в случае x0 = 0 — четыре треугольника

3 3 x0 + xi 3 x0 + x2

Xo

44

42

3 3 x2 + xi 3 x0 + x2

X2

44

2

42

33

4xo, 4x2, 4

3 xi + x2

2

3 3 x0 + xi 3 (xi + x2) 4xi ' 4 2 ' 4 2

3 x0 + xi 3 x0 + x2 3 xi + x2

4 2 ' 4 2 ' 4 2

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

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

Измельчение триангуляций можно реализовывать разными способами. Однако такое измельчение само по себе не гарантирует, что на каком-то шаге очередная задача линейного программирования будет иметь решение. Согласно [6] предложенный алгоритм построения последовательности триангуляций позволяет на каком-то шаге получить совместную задачу линейного программирования и тем самым получить функцию Ляпунова.

Система второго порядка

3. Пример

X 1 = — (Ьх2 + а)х 1, ,2

х 2 = 7X1 — Х2

имеет нулевое положение равновесия. Якобиан системы в нуле имеет вид

3(0, 0) =

—а 0 0

1

При а > 0 собственные значения этой матрицы отрицательны. Поэтому нулевое положение экспоненциально устойчиво [9]. Выберем значения параметров а = 2, Ь = 0.02, 7 = 0.01 и рассмотрим задачу построения кусочно-линейной функции Ляпунова для рассматриваемой системы.

Минимально значение В, вычисленное согласно (6), составляет 27 = 0.02. При В = 0.03 рассматриваемый алгоритм, описанный в разд. 2, приводит к решению задачи на множестве [—3; 3]2 за одну итерацию. На рис. 2 показана триангуляция Т0. На рис. 4 представлен график кусочно-линейной функции Ляпунова, полученный при В = 0.03.

Рис. 4. График функции Ляпунова при В = 0.03

При В = 0.32 построение решения происходит за две итерации. На рис. 3 показана триангуляция Т1. На рис. 5 представлен график соответствующей кусочно-линейной функции Ляпунова.

Дальнейшее увеличение параметра В приводит к увеличению количества итераций и росту размерности задачи линейного программирования. Так, при В = 1 для построения

Рис. 5. График функции Ляпунова при B = 0.32

функции Ляпунова двух итераций уже недостаточно. На третьей итерации матрица системы ограничений в задаче линейного программирования имеет размеры 23665 х 7505.

Заключение

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

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

Работа выполнена при поддержке РФФИ (гранты 17-07-00653, 16-07-00927).

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

1. Зубов В.И. Методы A.M. Ляпунова и их применение. Л.: Изд-во Ленингр. ун-та, 1957. 241 с.

2. Зубов В.И. Устойчивость движения (Методы Ляпунова и их применение): учеб. пособие.

2-е изд. М.: Высшая школа, 1984. 232 с.

3. Giesl P. Construction of global Lyapunov functions using radial basis functions. B.; N.Y.: Springer, 2007. 166 p. (Lecture Notes in Mathematics; vol. 1904). DOI: 10.1007/ 978-3-540-69909-5

4. Giesl P., Hafstein S. Existence of piecewise affine Lyapunov functions in two dimensions // J. of Mathematical Analysis and Applications. 2010. Vol.371, iss. 1. Pp. 233-248. DOI: 10.1016/j.jmaa.2010.05.009

5. Крищенко А.П. Исследование асимптотической устойчивости положений равновесия методом локализации инвариантных компактов // Автоматика и телемеханика. 2017. №6. С. 36-56.

6. Giesl P., Hafstein S. Construction of Lyapunov functions for nonlinear planar systems by linear programming // J. of Mathematical Analysis and Applications. 2012. Vol.338, iss. 1. Pp. 463-479. DOI: 10.1016/j.jmaa.2011.10.047

7. Marinosson S.F. Lyapunov function construction for ordinary differential equations with linear programming // Dynamical Systems. 2002. Vol.17, iss. 2. Pp. 137-150. DOI: 10.1080/0268111011011847

8. Marinosson S.F Stability analysis of nonlinear systems with linear programming: A Lyapunov functions based approach: PhD thesis. Duisburg, 2002. 103 p.

9. Халил Х.К. Нелинейные системы: пер. с англ. М.-Ижевск: Регулярная и хаотическая динамика, Институт компьютерных исследований, 2009. 832 с. [Khalil H. Nonlinear Systems. 3rd ed. Upper Saddle River: Prentice Hall, 2002.].

Mathematics i Mathematical Modelling

Electronic journal of the Bauman MSTU

Mathematics and Mathematical Modelling of the Bauman MSTU, 2017, no. 3, pp. 32-43.

DOI: 10.24108/mathm.0317.0000073

Received: 13.06.2017

© Bauman Moscow State Technical University

The Piecewise-Linear Lyapunov Function Construction for Dynamical Systems of the Second Order

Kasatkina T. S.1*

* kasatkina_t_s@mail.ru 1 Bauman Moscow State Technical University, Russia

Keywords: Lyapunov function, stability of the equilibrium point, dynamical systems, linear programming problem

The paper considers a stability problem of the equilibrium points of dynamical systems of the second order and describes an investigative technique based on the Lyapunov function construction. This method is useful for exploring the equilibrium points, which are exponentially stable.

One of the methods to analyse the equilibrium point stability of the systems of ordinary differential equations is search for the Lyapunov function. It is mandatory for the found function to meet the specific conditions. So far, there is no universal technique to construct the Lyapunov function. There are developed approaches, which allow constructing the Lyapunov functions in different forms. One of well-studied technique to construct the Lyapunov function is based on the Zubov's equation solution [1,2]. Its disadvantage is that it is necessary to solve a partial differential equation. A number of articles offer piecewise specifying of Lyapunov function at different meshes of the state space containing equilibrium point. Publications present studies of the Lyapunov function construction for a diversity of techniques to specify the triangle meshes.

The paper describes a construction method of the piecewise-linear Lyapunov function on the triangle mesh. Presents an algorithm of triangulation of the state space within which there is the equilibrium point. The suggested method is based on the solution of the linear programming problem. The variables of this problem are the values of the Lyapunov function in the nodes of the mesh and the additional constants, which ensure that function satisfies specific conditions. This method is iterative. The initial triangulation is built. If there is no solution on the current iteration, another triangulation is built, and finding a new problem solution with the new amount of the variables takes place.

Implementation of this method for different values of the parameter is analyzed for a dynamical system of the second order. It is shown that the number of required iterations depends on the parameter value. A lot of iterations involve essential calculus problems. The algorithm efficiency

considerably varies with respect to the eigenvalues of the system Jacobian matrix in the equilibrium point. Development of technique to construct the Lyapunov function, which is useful for wider class of dynamical systems, is planned.

References

1. Zubov V.I. Metody A.M. Liapunova i ikhprimenenie [Liapunov's methods and its applications]. Leningrad: Leningrad State Univ. Publ., 1957. 241 p. (in Russian).

2. Zubov V.I. Ustojchivost' dvizheniia (Metody Liapunova i ikh primenenie) [The stability of motion (Lyapunov's methods and their application)]: a textbook. 2nd ed. Moscow: Vysshaia shkola Publ., 1984. 232 p. (in Russian).

3. Giesl P. Construction of global Lyapunov functions using radial basis functions. B.; N.Y.: Springer, 2007. 166 p. (Lecture Notes in Mathematics; vol. 1904). DOI: 10.1007/ 978-3-540-69909-5

4. Giesl P., Hafstein S. Existence of piecewise affine Lyapunov functions in two dimensions. J. of Mathematical Analysis and Applications, 2010, vol.371, iss. 1, pp. 233-248. DOI: 10.1016/j.jmaa.2010.05.009

5. Krishchenko A.P. Investigation of asymptotic stability of equilibria by localization of the invariant compact sets. Automation and Remote Control, 2017, vol. 78, iss. 6, pp. 989-1005. DOI: 10.1134/S0005117917060030

6. Giesl P., Hafstein S. Construction of Lyapunov functions for nonlinear planar systems by linear programming. J. of Mathematical Analysis and Applications, 2012, vol.388, iss. 1, pp. 463-479. DOI: 10.1016/j.jmaa.2011.10.047

7. Marinosson S. Lyapunov function construction for ordinary differential equations with linear programming. Dynamical Systems, 2002, vol.17, iss. 2, pp. 137-150. DOI: 10.1080/ 0268111011011847

8. Marinosson S. Stability analysis ofnonlinear systems with linear programming: A Lyapunov functions based approach. PhD thesis. Duisburg, 2002. 103 p.

9. Khalil H. Nonlinear Systems. 3rd ed. Upper Saddle River: Prentice Hall, 2002. 750 p. (Russ. ed.: Khalil H. Nelinejnye sistemy. M.-Izhevsk: Reguljarnaja i haoticheskaja dinamika, Institut komp'juternyh issledovanij, 2009. 832 p.).