CC BY
58
Изв. вузов «ПНД», т. 19, № 5, 2011
УДК 537.86
ВОССТАНОВЛЕНИЕ СИСТЕМ НЕЙТРАЛЬНОГО ТИПА С ЗАПАЗДЫВАНИЕМ
А.С. Караваев, В.И. Пономаренко, М.Д. Прохоров
Предложены методы реконструкции систем с задержкой, моделируемых дифференциальными уравнениями нейтрального типа с запаздыванием, по временным рядам. Эффективность методов продемонстрирована на численных примерах при восстановлении обобщенного уравнения Маккея-Гласса и модельных уравнений, описывающих качку корабля и колебания тела вертикально стоящего человека.
Ключевые слова: Реконструкция уравнений, системы с запаздыванием, анализ временных рядов.
Введение
Моделирование и исследование многомерных колебательных систем со сложной динамикой по зашумленным временным реализациям экспериментально наблюдаемых величин является актуальной междисциплинарной задачей, имеющей большое значение при изучении многих реальных объектов, описание которых «из первых принципов» затруднено или невозможно. Решение этой задачи важно не только для описания динамики системы и прогноза ее поведения во времени и при изменении параметров. Оно позволяет также оценить адекватность заложенных в модели представлений об объекте, осуществить классификацию систем и режимов их функционирования, восстановить скрытые переменные и определить значения параметров, недоступных непосредственному измерению в эксперименте.
Использование при реконструкции модельных уравнений по временным рядам универсальных методик, не учитывающих особенности объекта, как правило, не приводит к успеху. На хороший результат обычно можно рассчитывать лишь при использовании специальных технологий реконструкции для определенных классов объектов. Такой класс, в частности, составляют системы с запаздыванием. Обычно они моделируются дифференциальными уравнениями с запаздывающим аргументом. Такие модели успешно применяются во многих разделах физики, биологии и
3
химии [1-4]. В достаточно общем случае системы с запаздывающей обратной связью описываются уравнением следующего вида:
SnX^Kt) +en-ix(n-r) (t) +---+ eix(t) = F (x(t),x(t - xi),...,x(t - xk)), (1)
где x(t) - состояние системы в момент времени t; x(n)(t) - производная по времени порядка n; x1,...,Xk - времена запаздывания; ei,...,en - параметры, характеризующие инерционные свойства системы; F - некоторая функция.
Так как системы с запаздыванием обладают бесконечно большим числом степеней свободы, это вносит трудности в решение задачи реконструкции по временному ряду. Даже уравнения первого порядка с запаздыванием могут демонстрировать хаотические колебания, которым в фазовом пространстве соответствуют аттракторы очень высокой размерности. Поэтому для реконструкции систем вида (1) с задержкой разрабатывают специальные методы [5-15].
Помимо систем, моделируемых уравнением (1), динамика которых зависит от состояния динамической переменной в задержанные моменты времени, существует класс систем с запаздыванием, динамика которых зависит от скорости изменения динамической переменной в задержанные моменты времени. Такие системы достаточно распространены в физике, биологии, инженерных науках [4,16-20]. Они описываются дифференциальными уравнениями нейтрального типа с запаздыванием. Наиболее часто используемыми моделями нейтрального типа являются дифференциальные уравнения второго порядка с запаздыванием
e2x(t) + eix(t) = F (x(t),x(t - x{),x(t - T2)). (2)
К таким уравнениям относится, например, обобщенное уравнение Маккея-Гласса [8]
x(t) = -уx(t) - U)lx(t) + о}§/ (x(t - Ti)) + - f), (3)
где нелинейная функция f имеет вид
/(*(«-*■» = ^(7-!,)■ <4>
К уравнению вида (2) относится также уравнение, описывающее малые колебания вертикально стоящего человека
x(t) = -ycc(t) - w2x(t) - k1x(t - т1) - k2x(t - t1), (5)
где у, , k1, k2 - положительные параметры [21]. Время запаздывания т1 соответ-
ствует нейрофизиологической задержке, возникающей вследствие конечной скорости распространения и обработки нервных сигналов в организме человека. Принято считать, что в уравнении (5) члены без задержки описывают пассивные компоненты системы управления движением, а члены с запаздыванием описывают активные компоненты [21].
При малой амплитуде колебаний уравнением вида (2) описывается также качка корабля, оборудованного для предотвращения качки на волнах специальными резервуарами, частично заполненными водой [20]
x(t) = -ycc(t) - rn0x(t) - kx(t - т1), (6)
4
где x(t) - угол отклонения судна от вертикального положения; ю > 0 - частота собственных колебаний; у > 0 - коэффициент затухания; к > 0 - коэффициент трения. Заметим, что уравнение (6) не зависит от x(t — х\).
В отличие от уравнения (3), способного демонстрировать хаотическую динамику, уравнения (5) и (6) описывают затухающие колебания. Однако под действием внешней силы, в том числе под действием шума, системы (5) и (6) могут демонстрировать сложную динамику. Задача реконструкции времени запаздывания таких систем по их временным рядам до настоящего времени не исследовалась. В данной работе мы предлагаем методы восстановления систем, описываемых дифференциальными уравнениями нейтрального типа с запаздыванием, и апробируем их на различных модельных уравнениях.
1. Описание методов
1.1. Восстановление времени запаздывания с помощью статистического анализа экстремумов временного ряда и их производных. Ранее нами было установлено, что во временных реализациях систем с запаздыванием вида (1) с одним временем задержки практически отсутствуют экстремумы, удаленные друг от друга на время запаздывания [9]. Если система (1) совершает хаотические колебания, то экстремумы в ее временном ряде расположены нерегулярно и расстояние между ними принимает различные значения. На основе этого свойства нами был предложен метод определения %\, использующий статистический анализ временных интервалов между экстремумами хаотического временного ряда системы с запаздыванием. Определив для различных значений т число N ситуаций, при которых точки временного ряда, разделенные интервалом времени т, одновременно являются экстремальными, и построив зависимость N(т), легко найти время задержки Ti как значение, при котором наблюдается абсолютный минимум этой зависимости [9].
Покажем, что этот метод определения времени задержки можно развить на системы нейтрального типа с запаздыванием. Рассмотрим сначала систему (6). Продифференцировав ее по t и перенеся член со второй производной в левую часть уравнения, получим
x(t) + yX(t) = —w>^Xc(t) — kX(t — т1). (7)
Как видно из уравнения (7), если при X(t) = 0 одна из производных X(t) или x (t) отлична от нуля или обе эти производные отличны от нуля и имеют одинаковый знак, то должно выполняться условие X(t — Ti) = 0. Для различных значений т построим зависимость M (т), где M - число пар точек, для которых одновременно выполняются два условия: X(t) = 0 и X(t — т) = 0. Причем при построении зависимости M(т) будем использовать только такие экстремальные точки X(t) = 0, в которых вторая и третья производная имеют одинаковый знак. Зависимость M (т) будет иметь минимум при времени т, соответствующем времени запаздывания системы.
Рассмотрим теперь уравнение (5). Его производная по времени имеет вид
X (t) + yX(t) = — rn0X(t) — k1ic(t — т1) — k2 X(t — т1). (8)
5
Аналогичный вид имеет производная по времени уравнения (3). При x(t) = 0 и одинаковом знаке x(t) и (t) сумма второго и третьего слагаемых в правой части (8) отлична от нуля. А это возможно, когда x(t — т\) = 0 и/или x(t — ti) = 0. В этом случае график M (т) будет иметь выраженный минимум при истинном значении времени запаздывания.
Отметим, что для восстановления времени запаздывания Ti уравнений (3) и (5) можно использовать зависимость N(т), предложенную нами для систем с запаздыванием (1) [9]. Это объясняется тем, что в уравнениях (3) и (5) запаздывающая обратная связь такова, что динамическая переменная и ее производная задержаны на одно и то же время Ti.
Методы реконструкции времен запаздывания, основанные на построении и анализе зависимостей N(т) и M(т), предназначены для восстановления систем вида (1) и (2), находящихся в режиме хаотических колебаний. Однако, в отличие от уравнения (3), способного демонстрировать хаотическую динамику, уравнения (5) и (6) описывают затухающие колебания. Воздействие на системы (5) и (6) шумом способно индуцировать в них колебания, при этом системы могут демонстрировать сложное поведение [20]. В этом случае может быть использован предложенный метод восстановления времени запаздывания, основанный на статистическом анализе экстремумов временного ряда и их производных. Отметим, что включение шума в уравнения нарушает связи между производными переменных по времени. Кроме того, оценка самих производных (особенно высокого порядка) по временному ряду становится все менее точной с увеличением уровня шума. Эти обстоятельства накладывают ограничения на применимость метода, который перестает быть эффективным, начиная с некоторого уровня шума.
1.2. Восстановление времени запаздывания по отклику системы на периодическое внешнее воздействие прямоугольным импульсным сигналом. Для
систем с запаздыванием, описываемых уравнением вида (1), нами недавно был предложен метод восстановления времени задержки, основанный на возмущении системы периодическим внешним воздействием и анализе отклика [13]. Если воздействие y(t) имеет вид прямоугольных импульсов, то траектория x(t) системы (1) будет претерпевать возмущения через время т (i = 1,...к) после начала и после окончания каждого импульса. При этом на реализации x(t) появляются изломы, которые при малой амплитуде воздействия практически не заметны. Изменение динамики системы становится более заметным, если численно продифференцировать x(t) по времени. Зависимость x(t) демонстрирует скачок через время т после прохождения переднего и заднего фронтов прямоугольного импульса. Наиболее удобна для анализа отклика системы с запаздыванием (1) на внешнее воздействие вторая производная x(t). Ее временная реализация демонстрирует узкие пики или провалы через время т после начала и окончания импульса, хорошо заметные даже при малых амплитудах воздействия. Взаимная корреляционная функция
C (s)
(№)\|&(t + s)|) '(iSWI2) (l*(i)|2)
(9)
где угловые скобки обозначают усреднение по времени, имеет четко выраженные максимумы при s = т^, i = 1,...к [13].
6
Аналогичные рассуждения остаются справедливыми и для случая периодического воздействия прямоугольными импульсами y(t) на системы нейтрального типа с запаздыванием. Временные реализации X(t) таких систем, описываемых уравнениями (3), (5), (6), имеют хорошо заметные пики и провалы через время Ti после начала и окончания импульса. Следовательно, для систем нейтрального типа с запаздыванием функция (9) также будет демонстрировать максимумы при s = Ti и может быть использована для восстановления времени задержки.
Помимо пиков при s = Ti функция C(s) имеет другие пики, положение которых определяется отношениями T/т1 и Q/T , где T - период, а Q - длительность импульсного воздействия [13]. То есть при разных T пики C(s) наблюдаются при разных s. Неизменным остается лишь положение пика, соответствующего времени запаздывания. Поэтому для определения времени запаздывания следует подействовать на исследуемую систему сначала импульсным сигналом с T = Ti, а затем с T = T2, близким к Ti, и сравнить функции C(s) для первого и второго случаев. Обнаружив пик C(s), положение которого не меняется при изменении T, мы найдем Ti.
1.3. Восстановление линейных и нелинейных функций и параметров инерционности. Определив тем или иным способом время задержки Ti исследуемой системы нейтрального типа с запаздыванием, мы можем восстановить параметры инерционности ei и £2 и функции от задержанных переменных. Для этого будем использовать следующий подход. Рассмотрим сначала уравнение (6). Поделим его на и перепишем в виде
£2X(t) + £iX(t) + x(t) = f (X(t - Ti)) , (10)
где e2 = l/, ei = у/ш^, f - линейная функция. Из уравнения (10) следует, что, если построить на плоскости множество точек с координатами (X(t — Ti),£2X(t) + +eiX(t) +X(t)), то оно воспроизведет функцию f. Поскольку заранее величины ei и е2 не известны, будем строить зависимости e2X(t) + £iX(t) + X(t) от X(t — Ti) для различных значений e i и £2, добиваясь однозначности, которая возможна только при £ i = £Ъ £2 = £2.
В качестве количественного критерия однозначности при таком поиске ei и £2 будем использовать минимальную длину L (£i, £2) линии, последовательно соединяющей точки на перебираемых плоскостях вложения, упорядоченные по величине абсциссы. Минимум зависимости L (£i ,£2) будет наблюдаться при правильном выборе параметров, а построенное при этом значении множество точек на плоскости воспроизведет функцию f. При ошибочном выборе значений £i и £2 мы получаем на плоскости набор точек, не связанных между собой функционально. Чем менее точно определены параметры, тем более беспорядочно расположены точки, а соединяющая их ломаная линия имеет большую длину, чем в случае, когда множество точек ложится на одномерную кривую. Похожий подход использовался нами в [9] при реконструкции систем с запаздыванием вида (1). Отметим, что применение такого подхода возможно только в том случае, если нам априорно известен вид модельного уравнения, а неизвестными являются параметры.
Восстановив параметры £i и £2, можно найти шо и у: шо = л/1/£2 , у = £i/£2. Параметр к уравнения (6) можно оценить по углу наклона восстановленной линейной функции f.
7
Процедура реконструкции параметров уравнения (5) несколько сложнее. Перепишем уравнение (5) в виде
£2x(t) + eix(t) + x(t) = /1 (x(t — ti)) + /2 (x(t — ti)) , (11)
где ei и e2 имеют тот же вид, что и в уравнении (10), /1 и /2 - линейные функции. Уравнение (11) описывает двумерную поверхность, на которой располагаются точки, посещаемые системой (5) в трехмерном пространстве вложения (x(t — тi),x(t — t1),
e2x(t) + ei x(t) + x(t)).
Сечение этой поверхности плоскостью x(t — ti) = const позволяет с точностью до константы определить вид функции /1, так как точки сечения удовлетворяют уравнению e2x(t) + e1x(t) + x(t) = /1 (x(t — t1)) + c1, где c1 = /2 (x(t — t1)) при выбранном постоянном значении x(t — Т1). Аналогично можно восстановить с точностью до константы функцию /2, построив сечение x(t — Т1) = const, точки которого описываются уравнением e2x(t) + e1x(t) + x(t) = /2 (x(t — t1)) + c2, где c2 = /1 (x(t — t1)) при фиксированном x(t — t1).
Так как параметры e1 и e2 априорно не известны, будем перебирать их из некоторого интервала, находя в сечении x(t — t1) = const или x(t — t1) = const зависимость, наиболее близкую к однозначной. Как и в рассмотренном выше случае, в качестве количественного критерия однозначности удобно использовать минимальную длину L (£1,62) линии, соединяющей упорядоченные по величине абсциссы точки сечения, где e1 и ^2 - пробные значения параметров e1 и e2. Минимум L (£1, e 2) будет наблюдаться при e 1 = e1, £2 = e2, при этом в указанных сечениях получим функции /1 и /2. После этого описанным выше способом можно вычислить юо и у и оценить параметры к1 и k2.
Аналогичным образом можно восстановить параметры уравнения (3), которое тоже можно записать в виде (11). Отличие от уравнения (5) будет в том, что функции /1 и /2, входящие в (11), являются нелинейными.
2. Применение методов
2.1. Восстановление обобщенного уравнения Маккея-Гласса. Применим описанные в разделе 1 методы для восстановления по наблюдаемому временному ряду обобщенного уравнения Маккея-Гласса (3) в присутствии гауссова S-коррелированного шума Z(t). Уравнение исследуемой системы имеет вид
x(t) = -уx{t) - (£>lx(t) + СОо/ (x(t - Ti)) + - Ti) + (12)
где {Z(t)Z(tr)) = 2DS(t — t'), D - интенсивность шума. Зададим следующие значения параметров: у = 1.5, ю0 = 1.0, т1 = 100.0, a = 3.0, c = 10.0. При этих параметрах в отсутствие шума система демонстрирует движение на хаотическом аттракторе высокой размерности [8]. На рис. 1, а приведен фрагмент временного ряда колебаний системы (12) при D = 0.014. Весь ряд состоял из 106 точек при шаге интегрирования h = 0.01 и содержал около 6700 экстремумов.
8
Подсчитав число M одновременных обращений в нуль x(t) и x(t — т) для различных значений т, перебираемых с шагом 0.1, построим зависимость M (т), введя нормировку M на общее число экстремумов в ряде (рис. 1, б). Напомним, что для построения M (т) используются только такие экстремальные точки x(t) = 0, в которых вторая и третья производные имеют одинаковый знак. Для оценки производных по временному ряду мы использовали локальную аппроксимацию полиномами 4 порядка, подгоняемыми методом наименьших квадратов. Отметим, что для сглаживания временного ряда и уменьшения числа обусловленных шумом экстремумов при оценке производной в процедуре локальной аппроксимации использовалось 7 соседних точек. Несмотря на присутствие шума, зависимость M (т) демонстрирует четко выраженный минимум при т = 100.2, давая хорошую оценку времени запаздывания (см. рис. 1, б). Относительная погрешность А при восстановлении т составляет 0.2%.
Для иллюстрации метода восстановления времени запаздывания, основанного на анализе отклика системы на внешнее воздействие, подействуем на систему (12) внешним сигналом y(t), представляющим собой прямоугольные импульсы. На рис. 1, в построены взаимные корреляционные функции (9) для случая, когда возмущающие импульсные сигналы имеют амплитуду A = 5, периоды Тх = 250, Т2 = 300 и длительность Q = 0.1, а параметры системы (12) такие же, как в рассмотренном выше случае. При шаге изменения s, равном 0.01, C(s) имеет высокий пик при s = тх = 100.00, который не смещается при изменении периода воздействующих импульсов.
Восстановим теперь параметры у и w0 и нелинейную функцию уравнения (12) по ее временному ряду в отсутствие внешнего воздействия. Возьмем точки экстремального сечения X(t — тх) = 0 при восстановленном тх = 100 и будем проецировать их на плоскости (x(t — тх), x(t)/w2 + yx(t)/w2 + x(t)) при различных значениях у и w, перебираемых с шагом 0.01 из некоторого интервала, добиваясь зависимости, наиболее близкой к однозначной. В отсутствие шума (D = 0) длина L (у, w) линии, последовательно соединяющей точки на перебираемых плоскостях вложения, упорядоченных по величине абсциссы, имеет минимум при у = 1.50 и w = 1.00, то есть при истинных значениях у и w0 .На рис. 2, а приведена восстановленная нелинейная функция системы. При D = 1.4 • 10-4 минимум L (у, w) наблюдается при у = 2.33
Рис. 1. а - временной ряд обобщенного уравнения Маккея-Гласса (12) в режиме хаотических колебаний; б - зависимость M(т), нормированная на общее число экстремумов в ряде, Mmin (т) = = M(100.2); в - взаимные корреляционные функции (9) при Т1 = 250 (сплошная линия) и Т2 = 300 (пунктирная линия). Пики C(s) при s = 100.00 совпадают для Ti и Т2
9
Рис. 2. Восстановленная нелинейная функция обобщенного уравнения Маккея-Гласса (12): а - у = 1.50, со = 1.00, D = 0; 6 - у = 2.33, со = 1.31, D = 1.4 • 1СГ4
0 40 80 120 160 t ^ 0 10 20
Рис. 3. а - временной ряд обобщенного уравнения Маккея-Гласса (3) в режиме периодических колебаний; б - взаимные корреляционные функции (9) при Т = 20 (сплошная линия) и Т2 = 25 (пунктирная линия). При s = 10.0 пики C(s) для Т =20 и Т2 = 25 совпадают
(А = 55% ) и w = 1.31 (А = 31% ). При этих параметрах восстановленная нелинейная функция имеет вид, представленный на рис. 2, б. С увеличением уровня шума качество реконструкции параметров у и wo и нелинейной функции ухудшается.
Следует отметить, что предложенный нами метод реконструкции времени задержки, основанный на анализе отклика системы на внешнее воздействие, можно применять к системам нейтрального типа, находящимся как в режиме хаотических, так и в режиме периодических колебаний. Пусть параметры уравнения (3) у = 1.5, w0 = 1.0, Ti = 10.0, a = 0.9, c =10 таковы, что система демонстрирует периодические колебания (рис. 3, а). Подействуем на исследуемую систему периодическими импульсными сигналами с амплитудой A = 1, периодами Ti = 20, T2 = 25 и длительностью Q = 0.01. Из рис. 3, б видно, что C(s) имеет высокий пик при s = Ti = 10.0, положение которого не изменяется при изменении периода воздействующих импульсов, то есть время запаздывания системы восстанавливается точно.
2.2. Восстановление уравнения, моделирующего качку корабля. Как было отмечено выше, модельное уравнение (6), описывающее при малой амплитуде колебаний качку корабля на волнах, демонстрирует в отсутствие внешнего воздействия затухающие колебания. Однако воздействие на систему шумом способно индуцировать в ней сложную колебательную динамику [20]. Подействуем на систему (6) гауссовым S-коррелированным шумом Z(t) и запишем уравнение возмущаемой системы в виде
X(t) = —ycfc(t) — w0x(t) — kx(t — T1) + Z(t). (13)
10
На рис. 4, а приведен фрагмент временного ряда колебаний системы (13) при Y = 2, ю0 = 0.8 , к = 1, т1 = 50, D = 0.1, h = 0.01. Весь ряд состоял из 105
точек и содержал около 4100 экстремумов. Зависимость M(т), построенная при шаге изменения т, равном 0.1, приведена на рис. 4, б. Абсолютный минимум M (т) наблюдается при т = 50.1 (А = 0.2%).
Исследования показывают, что с увеличением значений параметров у и к абсолютный минимум M(т) становится более выраженным, и при у и к, больших 10, метод позволяет точно восстановить т1. Зависимость M(т), построенная при у = 20, к = 20 и тех же юо, т1 и D, что и рис. 4, б, демонстрирует хорошо заметный абсолютный минимум при т = 50.0 (рис. 4, в).
Рис. 4, г иллюстрирует восстановление времени запаздывания по отклику системы на периодическое внешнее воздействие в виде прямоугольного импульсного сигнала. Взаимные корреляционные функции (9) построены на рис. 4, г при шаге изменения s, равном 0.01, для случаев, когда возмущающие импульсные сигналы имеют амплитуду A = 1, периоды Ti = 200, T2 = 250 и длительность Q = 0.01. Из рис. 4, г видно, что C(s) имеет высокий максимум при s = т1 = 50.00, положение которого не изменяется при изменении периода воздействующих импульсов, то есть время запаздывания системы восстанавливается точно. Отметим, что метод определения т1 с помощью возбуждающих импульсов можно применять непосредственно к системе, моделируемой уравнением (6), при D = 0, поскольку внешний импульсный сигнал выводит систему из положения равновесия и вызывает в ней колебания.
Проиллюстрируем реконструкцию остальных параметров уравнения (13) по временному ряду, полученному при следующих значениях параметров: у = 2, юо = 0.8, к = 1, т1 = 50, D = 0.1, для случая, когда система возбуждается импульсами с амплитудой A =1, периодом T = 250 и длительностью Q = 0.01. 11
Мт;п(т) = М(50.1); в - зависимость М(т), нормированная на общее число экстремумов в ряде, при Y = 20, юо = 0.8, к = 20, тх = 50, D = 0.1, Мт;п(т) = M(50.0); г - взаимные корреляционные функции (9) при Ti = 200 (сплошная линия) и Т2 = 250 (пунктирная линия). При s = 50.00 пики C(s) для Т = 200 и Т2 = 250 совпадают
11
Для восстановления параметров у и юо системы мы проецировали ее траекторию в промежутках между внешними импульсами на плоскости (x(t — ti),X + yx + ю2х) для различных значений у и ю, перебираемых с шагом 0.01 из некоторого интервала, а значение времени запаздывания Ti считали уже восстановленным одним из описанных выше способов. При этом мы искали однозначную зависимость в перебираемых пространствах вложения. Длина L (у, ю) линии, последовательно соединяющей точки на перебираемых плоскостях вложения, упорядоченных по величине абсциссы, имеет минимум при Y = у = 2.00 и ю = ю0 = 0.80. При этих параметрах проекция траектории системы имеет вид, представленный на рис. 5. Определив угловой коэффициент прямой, аппроксимирующей набор точек на рисунке, получаем оценку к = 1.0, что в точности совпадает с истинным значением.
2.3. Восстановление модельного уравнения колебаний вертикально стоящего человека. Уравнение (5), описывающее малые колебания вертикально стоящего человека, так же как и уравнение (6), моделирующее колебания корабля при качке, демонстрирует затухающие в отсутствие внешнего воздействия колебания. Поэтому, как и в рассмотренном выше случае, подействуем на систему (5) гауссовым S-коррелированным шумом Z(t), индуцирующим в ней сложные колебания, и запишем уравнение возмущаемой системы в виде
x(t) = —yA(t) — ю^х(к) — k1x(t — т1) — k2x(t — т1) + Z(t). (14)
На рис. 6, а приведен временной ряд колебаний системы (14) при у = 10, ю0 = 0.9, к1 = 0.5, к2 = 10, т1 = 50.0, D = 1, h = 0.01. Весь ряд состоял из
105 точек и содержал около 3700 экстремумов. Зависимость M(т), построенная при шаге изменения т, равном 0.1, приведена на рис. 6, б. Абсолютный минимум M(т) наблюдается при т = 50.0, позволяя точно определить время задержки.
Подействуем на систему (14) внешним прямоугольным импульсным сигналом y(t) и построим функцию (9) (рис. 6, в). При A = 1, T1 = 200 и T2 = 250, Q = 0.01, D = 0.1 и шаге изменения s, равном 0.01, C(s) имеет высокий пик при s = T1 = = 50.00, который не смещается при изменении периода воздействующего импульсного сигнала. Таким образом, T1 восстанавливается точно.
Восстановление остальных параметров модельного уравнения (14) по временному ряду проиллюстрируем при их значениях, указанных выше, для случая, когда система возбуждается в отсутствие шума (D = 0) импульсами с A = 1, T = 250 и Q = 0.01. Процедура реконструкции состоит из нескольких этапов. Сначала построим экстремальное сечение x(t — т1) =0 реализации, сгенерированной исследуемой системой. При этом полагаем, что время запаздывания T1 = 50 мы уже восстановили описанным выше способом. В отсутствие внешнего воздействия точки этого 12
Рис. 5. Проекция траектории системы (13) на плоскость (X(t — ti),Х + уХ + ю2х) при у = 2.00, ю = 0.80
12
сечения удовлетворяют уравнению x(t) + yx(t) + u>2x(t) = —kix(t — xi). Но у, ю0 и ki нам заранее не известны, поэтому мы проецируем точки экстремального сечения на плоскости (x(t — xi),X(t) + yX(t) + w2x(t)) для различных значений у и ю, перебираемых с шагом 0.01 из некоторого интервала, и ищем такие у и ю, при которых множество точек на плоскости ложится на прямую линию. То есть в качестве оценок у и юо мы берем такие у и ю, при которых длина L (у, ю) линии, соединяющей упорядоченные по величине абсциссы точки проекции, оказывается минимальной. Полученные оценки у = 9.90 (А = 1%) и ю = 0.88 (А = 2%) близки к истинным у и юо. Отметим, что для уменьшения влияния на результаты реконструкции переходных процессов, обусловленных воздействием на систему импульсов, мы не используем при проецировании точки реализации, попадающие в интервал длительностью ti/5 после каждого импульса.
Рис. 6. а - временной ряд уравнения (14), моделирующего колебания вертикально стоящего человека; б - зависимость М(т), нормированная на общее число экстремумов в ряде, МШт(т) = = М(50.0); в - взаимные корреляционные функции (9) при Т1 = 200 (сплошная линия) и
Т2 = 250 (пунктирная линия). При s = 50.00 пики C(s) для Т1 = 200 и Т2 = 250 совпадают
Аппроксимируем множество точек проекции, построенной на рис. 7, а при у = 9.90 и ю = 0.88, прямой линией. Угловой коэффициент аппроксимирующей прямой у1 = 0.46 (А = 8% ) дает достаточно хорошую оценку параметра ki.
Для оценки параметра k2 спроецируем траекторию системы (14) на плоскость (x(t — т1), x(t)+yx(t)+ ^:>2x(t) + kix(t — т1^ при у = 9.90, ю = 0.88 и
Рис. 7. Реконструкция параметров к1 и к2 системы (14): а - проекция точек экстремального сечения x(t — т1) = 0 на плоскость (x(t — T1),x(t) + yx(t) + ю2х(Ь)) при у = 9.90, ю = 0.88; б - проекция
траектории системы (14) на плоскость I X(t — T1),X(t) + yx(t) + rn2x(t) + k1x(t — т1)) при у = 9.90,
ю = 0.88, к1 = 0.46
13
к1 = 0.46 (рис. 7, б). Затем аппроксимируем множество точек проекции прямой линией и найдем ее угловой коэффициент. Его значение к2 = 9.9 (А = 1% ) близко к истинному к.2 = 10.
Заключение
Нами предложены методы восстановления модельных дифференциальных уравнений с задержкой для систем нейтрального типа с запаздыванием по их временным рядам. Методы могут быть применены для реконструкции систем нейтрального типа, находящихся как в режиме хаотических, так и в режиме периодических колебаний. Для восстановления времени запаздывания использован метод, основанный на статистическом анализе экстремумов временного ряда и их производных, и метод, основанный на анализе отклика системы на периодическое внешнее воздействие в виде прямоугольных импульсов. Для восстановления линейных и нелинейных функций и параметров инерционности систем нейтрального типа с запаздыванием использован метод, основанный на проецировании траектории системы в специальным образом выбранные подпространства малой размерности.
Эффективность методов продемонстрирована на численных примерах при восстановлении обобщенного уравнения Маккея-Гласса в хаотическом и периодическом режимах и при реконструкции модельных уравнений, описывающих качку корабля и колебания тела вертикально стоящего человека.
Работа выполнена при поддержке РФФИ, грант № 10-02-00980, аналитической ведомственной целевой программы «Развитие научного потенциала высшей школы (2009-2011 годы)», проект № 2.1.1/1738, программы РАН и гранта Президента РФ для поддержки молодых кандидатов наук.
Библиографический список
1. Mackey M.C., Glass L. Oscillations and chaos in physiological control systems // Science. 1977. Vol. 197. P. 287.
2. Ikeda K. Multiple-valued stationary state and its instability of the transmitted light by a ring cavity system // Opt. Commun. 1979. Vol. 30. P. 257.
3. Epstein I.R. Delay effects and differential delay equations in chemical-kinetics // Int. Rev. in Phys. Chem. 1992. Vol. 11. P. 135.
4. Kuang Y. Delay Differential Equations with Applications in Population Dynamics. Boston: Academic Press, 1993.
5. Voss H., Kurths J. Reconstruction of non-linear time delay models from data by the use of optimal transformations // Phys. Lett. A. 1997. Vol. 234. P. 336.
6. Tian Y.-C., Gao F. Extraction of delay information from chaotic time series based on information entropy // Physica D. 1997. Vol. 108. P. 113.
7. Hegger R., BUnner M.J., Kantz H., Giaquinta A. Identifying and modeling delay feedback systems // Phys. Rev. Lett. 1998. Vol. 81. P. 558.
8. BUnner M.J., Ciofini M., Giaquinta A., Hegger R., Kantz H., Meucci R., Politi A. Reconstruction of systems with delayed feedback: (I) Theory // Eur. Phys. J. D. 2000. Vol. 10. P. 165. 14
14
9. Пономаренко В.И., Прохоров М.Д., Караваев А.С., Безручко Б.П. Определение параметров систем с запаздывающей обратной связью по хаотическим временным реализациям // ЖЭТФ. 2005. Т. 127. Вып. 3. С. 515.
10. Ortin S., Gutierrez J.M., Pesquera L., Vasquez H. Nonlinear dynamics extraction for time-delay systems using modular neural networks synchronization and prediction // Physica A. 2005. Vol. 351. P. 133.
11. Siefert M. Practical criterion for delay estimation using random perturbations // Phys. Rev. E. 2007. Vol. 76. 026215.
12. Yu D., Frasca M., Liu F. Control-based method to identify underlying delays of a nonlinear dynamical system // Phys. Rev. E. 2008. Vol. 78. 046209.
13. Prokhorov M.D., Ponomarenko V.I. Reconstruction of time-delay systems using small impulsive disturbances // Phys. Rev. E. 2009. Vol. 80. 066206.
14. Zunino L., Soriano M.C., Fischer I., Rosso O.A., Mirasso C.R. Permutation-information-theory approach to unveil delay dynamics from time-series analysis // Phys. Rev. E. 2010. Vol. 82. 046212.
15. Ma H., Xu B., Lin W., Feng J. Adaptive identification of time delays in nonlinear dynamical models // Phys. Rev. E. 2010. Vol. 82. 066210.
16. Gopalsamy K. Oscillations in neutral delay-differential equations // J. Math. Phys. Sci. 1987. Vol. 21. P. 23.
17. Gopalsamy K. Stability and Oscillations in Delay Differential Equations of Population Dynamics. Dordrecht: Kluwer, 1992.
18. Hale J.K., Lunel S.M.V Introduction to Functional Differential Equations. New York: Springer, 1993.
19. Bocharov G.A., Rihan F.A. Numerical modelling in biosciences using delay differential equations // J. Comp. Appl. Math. 2000. Vol. 125. P. 183.
20. Patanarapeelert K., Frank T.D., Friedrich R., Beek P.J., Tang I.M. A data analysis method for identifying deterministic components of stable and unstable time-delayed systems with colored noise // Phys. Lett. A. 2006. Vol. 360. P. 190.
21. Peterka R.J.Sensorimotor integration in human postural control // J. Neurophysiol. 2002. Vol. 88. P. 1097.
Саратовский госуниверситет Поступила в редакцию 14.06.2011
Саратовский филиал ИРЭ После доработки 19.10.2011
им. В.А. Котельникова РАН
RECONSTRUCTION OF NEUTRAL TIME-DELAY SYSTEMS
A.S. Karavaev, V.I. Ponomarenko, M.D. Prokhorov
The methods are proposed for the reconstruction of time-delay systems modeled by neutral delay-differential equations from their time series. The methods are successfully applied to the recovery of generalized Mackey-Glass equation and equations modeling ship rolling and human movement from simulated data.
Keywords: Reconstruction of equations, time-delay systems, time series analysis. 15
15
Караваев Анатолий Сергеевич - родился в 1981 году в Саратове. Окончил Саратовский государственный университет в 2004 году. Кандидат физикоматематических наук (2007). Доцент кафедры динамического моделирования и биомедицинской инженерии СГУ. Область научных интересов - моделирование по временным рядам, нелинейная динамика и ее приложения. Имеет более 50 научных публикаций.
410012 Саратов, ул. Астраханская, 83 Саратовский госуниверситет E-mail: karavaevas@gmail.com
Пономаренко Владимир Иванович - родился в 1960 году в Саратове. Окончил Саратовский государственный университет в 1982 году. Защитил диссертацию на соискание ученой степени кандидата физико-математических наук (1992) и доктора физико-математических наук (2008). Ведущий научный сотрудник Саратовского филиала Института радиотехники и электроники РАН. Область научных интересов - статистическая радиофизика, анализ временных рядов, нелинейная динамика и ее приложения. Автор более 130 научных публикаций.
410019 Саратов, ул. Зеленая, 38
Саратовский филиал Института радиотехники и электроники РАН E-mail: ponomarenkovi@gmail.com
Прохоров Михаил Дмитриевич - родился в 1968 году в Саратове. Окончил Саратовский государственный университет в 1992 году. Защитил диссертацию на соискание ученой степени кандидата физико-математических наук (1997) и доктора физико-математических наук (2008). Ведущий научный сотрудник Саратовского филиала Института радиотехники и электроники РАН. Область научных интересов - нелинейная динамика и ее приложения, математическое моделирование, анализ временных рядов. Имеет более 100 научных публикаций.
410019 Саратов, ул. Зеленая, 38
Саратовский филиал Института радиотехники и электроники РАН E-mail: mdprokhorov@yandex.ru 16
16
