CC BY
5Нелинейная динамика. 2013. Т. 9. № 4. С. 613-625. Полнотекстовая версия в свободном доступе http://nd.ics.org.ru
УДК: 537.86 М8С 2010: 37М10
Реконструкция систем с запаздыванием под внешним периодическим воздействием
М.В.Сысоева, В. И. Пономаренко, М.Д.Прохоров, И.В.Сысоев
Предложен метод восстановления систем первого порядка с запаздыванием, находящихся под внешним периодическим воздействием, по их временным рядам. Особенностью метода является учет структуры уравнения системы при построении регрессионной модели. Метод позволяет восстановить время запаздывания, параметр инерционности, нелинейную функцию системы, а также амплитуду и частоту внешнего периодического воздействия. Эффективность метода продемонстрирована на численных примерах при реконструкции различных неавтономных систем с задержкой.
Ключевые слова: реконструкция уравнений, системы с запаздыванием, анализ временных рядов
Получено 21 мая 2013 года
После доработки 9 октября 2013 года
Работа выполнена при поддержке РФФИ, гранты №13-02-00227 и №12-02-00377, а также гранта Президента РФ, МК-4435.2012.8.
Сысоева Марина Вячеславовна bobrichek@mail.ru Сысоев Илья Вячеславович ivssci@gmail.com
Саратовский государственный университет им. Н.Г.Чернышевского 410012, Россия, г. Саратов, ул. Астраханская, д. 83
Пономаренко Владимир Иванович ponomarenkovi@gmail.com
Саратовский государственный университет им. Н.Г.Чернышевского 410012, Россия, г. Саратов, ул. Астраханская, д. 83
Саратовский филиал Института радиотехники и электроники им. В. А. Котельникова РАН 410019, Россия, г. Саратов, ул. Зеленая, д. 38
Прохоров Михаил Дмитриевич mdprokhorov@yandex.^
Саратовский филиал Института радиотехники и электроники им. В. А. Котельникова РАН 410019, Россия, г. Саратов, ул. Зеленая, д. 38
1. Введение
Задача реконструкции динамических моделей (такой термин принят в нелинейной динамике [1]), или идентификации систем (в терминах математической статистики [2]), происходит из желания исследователей по доступным в ходе наблюдения или активного эксперимента данным подобрать оператор эволюции для изучаемого объекта. Чаще всего данные представляются в виде временных рядов — дискретных последовательностей значений одной (скалярный ряд) или нескольких (векторный ряд) наблюдаемых величин от времени. Если наблюдение представляет собой измерение, то обычно такие ряды снимаются с помощью аналого-цифровых преобразователей через равные промежутки времени (эквидистантный ряд); такие ряды получаются при работе некоторых радиофизических, медицинских и многих других измерительных и диагностических приборов.
Задача реконструкции, как и всякая обратная задача, в общем случае некорректна и может иметь множество различных решений, поскольку наблюдаемые ряды не несут полной информации о породившей их системе, а кроме того искажены помехами и шумами. Поэтому часто задачу реконструкции классифицируют по степени сложности, или по объему доступной об объекте моделирования дополнительной информации [3]. Можно условно выделить следующие варианты: 1) задача о «черном ящике», когда кроме временного ряда ничего более о системе не известно, 2) задача о «сером ящике», когда имеется дополнительная информация, как правило, помогающая определить принадлежность исследуемой системы к некоторому более или менее специальному классу, конкретизировать, сузить задачу, 3) задача о «белом ящике», когда оператор эволюции известен из первых принципов, так что остается только определить параметры и, возможно, измерить скрытые переменные. Следует отметить, что задача о черном ящике, как правило, не решаема, а задача о белом ящике встречается на практике достаточно редко (один из успешных примеров см. [4]), хотя оба типа задач и представляют большой практический интерес.
К настоящему времени накоплено достаточно большое число результатов по реконструкции в различных специальных случаях, в частности: для коротких сигналов [5], для защиты информации [6], для неавтономных систем, описываемых обыкновенными дифференциальными уравнениями, под периодическим [7] и квазипериодическим [8] воздействием, для стохастических систем на основе байесова подхода [9, 10], при наличии скрытых переменных [11, 12].
В данной работе рассматривается задача реконструкции уравнений состояния в случае, когда изучаемый объект может быть описан уравнениями с запаздывающей обратной связью. Системы, для которых уравнения с запаздыванием могут послужить подходящей моделью, достаточно распространены в природе [13-15]. Например, реальные системы, описываемые уравнениями с запаздывающим аргументом, встречаются в биофизике (уравнение Маккея-Гласса, описывающее регуляцию состава крови [16, 17], уравнение барорефлектор-ной регуляции артериального давления [18]), в оптике (уравнение Икеды, используемое для описания динамики пассивного оптического резонатора [19]), в радиофизике (радиотехнический генератор с запаздывающей обратной связью [20]).
Неавтономные модели могут включать в себя также некоторое внешнее воздействие, в качестве которого может выступать дыхание (в задаче описания регуляции артериального давления [18]), принимаемый сигнал (в радиофизических или оптических системах передачи информации [21]) или сигнал управляющей системы (в задачах исследования синхронизации внешним сигналом [22]).
Так как системы с запаздыванием обладают бесконечномерным фазовым пространством, это затрудняет решение задачи реконструкции по временному ряду с помощью уни-
версальных алгоритмов, обзор которых сделан в [1]. Действительно, при больших временах задержки даже скалярные дифференциальные уравнения с запаздыванием могут демонстрировать хаотические движения очень высокой размерности. Поэтому для реконструкции систем с запаздыванием разрабатываются специальные методы [23-36]. Однако задача восстановления модельных уравнений систем с запаздыванием по временным рядам в случае, когда эти системы находятся под внешним воздействием, до сих пор остается малоисследованной. Вообще говоря, такая ситуация типична для многих практически важных задач [37, 38]. Следует отметить, что специализированные подходы активно разрабатываются также для реконструкции систем, описываемых неавтономными обыкновенными дифференциальными уравнениями [8].
В данной работе мы предлагаем метод восстановления систем с запаздывающей обратной связью, находящихся под силовым внешним периодическим воздействием, по их временным рядам. Метод основан на использовании априорной информации о виде модельного уравнения системы, структура которого, в том числе наличие внешнего воздействия, учитывается при построении модели.
2. Описание метода 2.1. Идея метода
Рассмотрим уравнение системы с запаздыванием порядка Б в достаточно общем виде:
где То — время запаздывания, f — нелинейная функция. К уравнению вида (2.1) сводятся, например, уравнение Маккея - Гласса [16] и уравнение Икеды [19], ставшие эталонными при исследовании систем с запаздыванием.
Пусть на систему (2.1) действует внешний периодический сигнал х(1), а модельное уравнение при этом принимает вид
Далее будем полагать, что порядок системы с запаздыванием Б и способ введения внешнего воздействия известны априорно, а параметры е и т, нелинейную функцию f и форму внешнего воздействия г(Ь) — его амплитуду, частоту, спектральный состав — требуется найти. Тогда для восстановления по временному ряду системы вида (2.2) мы предлагаем следующий подход:
• вместо задачи о реконструкции дифференциальных уравнений будем решать задачу о реконструкции их разностной схемы (то есть фактически отображения последова-ния), построенной методом Эйлера;
• будем аппроксимировать нелинейную функцию f степенным полиномом;
• время запаздывания будем оценивать одним из методов, предложенных в работах [29, 33], прежде чем восстанавливать остальные параметры;
• внешнее воздействие будем аппроксимировать тригонометрическим полиномом по аналогии с [8].
(2.1)
(2.2)
Далее опишем алгоритм подробнее на примере системы первого порядка (Б — 1); для систем более высоких порядков процедура аналогична с поправкою на большую сложность соответствующих формул. Ограничимся случаем аддитивного воздействия, имея в виду, что в целом аналогично [7] можно обобщить на случаи произвольного воздействия метод, разработанный для аддитивного воздействия. Уравнение системы с запаздыванием первого порядка может быть получено из (2.2) и традиционно записывается следующим образом:
е^- = ¡{х,(1),х(1-то)) +г(1), (2.3)
где е — параметр, характеризующий инерционные свойства системы, выделенный из функции f.
Перейдем от дифференциального уравнения (2.3) к разностному уравнению, для систем первого порядка получится следующее выражение:
^д)"-^ = ¡(х(1),ха-то)) +г(1), (2.4)
где АЬ — малое время, равное интервалу выборки. Уравнение (2.4) можно переписать в виде
х{1 + А1) = + ^¡(х,{1),х {I - т0) ) + (2-5)
Запишем уравнение (2.5) в виде отображения, включив слагаемое и коэффициент
в состав нелинейной функции и функции внешнего воздействия:
хп+1 — f (уХп+ •> (2-6)
где п — Ь/АЬ — дискретное время, т — То/АЬ — дискретное время задержки, штрих у Хп+1 означает, что имеется в виду не наблюдаемое, а прогнозируемое по наблюдаемым значение (в идеале они, конечно, должны совпасть). Функция f' представляется в виде полинома степени N: если нет дополнительной информации о свойствах исходной системы — полинома общего вида со всеми перекрестными членами, если такая информация имеется — специализированного полинома или суммы полиномов более простого вида.
Периодическое внешнее воздействие представим в виде тригонометрического полинома:
м
х'п — ^ (вк сов(кшпАЬ) + dk вт(кшпАЬ)), (2.7)
к=1
где вк и dк — коэффициенты тригонометрического полинома степени М, аппроксимирующего внешнее периодическое воздействие.
2.2. Оценка коэффициентов
Для нахождения коэффициентов модели (2.6) используется критерий наименьших квадратов, минимизирующий средний квадрат одношаговой ошибки прогноза построенных моделей:
— ((жп+1 - х'п+1(т))2), (2.8)
где угловые скобки обозначают усреднение по времени. Такая мера показывает, насколько хорошо модель описывает наблюдаемую реализацию, если начальные условия для одноша-гового прогноза выбираются из исходного временного ряда.
Поскольку модельное уравнение (2.6) имеет нелинейную зависимость от частоты и внешнего воздействия, необходимо использовать нелинейный метод наименьших квадратов, то есть прибегать к итерационным алгоритмам. Нами использован метод Ньютона, заключающийся в следующем: задаются начальные догадки для искомых коэффициентов модели, далее система в окрестности них линеаризуется и с помощью линейного метода находятся поправки к начальным догадкам, эти поправки добавляются к значениям параметров, после чего процедура повторяется до тех пор, пока значения поправок на очередном шаге не станут пренебрежимо малы.
В нашем случае задача упрощается тем, что при использовании полиномиальной аппроксимации f только один коэффициент (частота и) входит нелинейно, поэтому метод был модифицирован таким образом, что начальное условие задавалось и поправки вычислялись только для и, остальные коэффициенты рассчитывались на каждом шаге линейным методом. В качестве начального приближения для и была использована оценка по спектру мощности колебаний системы (2.3). Поскольку в спектре мощности неавтономной системы с запаздыванием в общем случае сложно выделить пик, соответствующий периодическому внешнему воздействию, особенно если амплитуда воздействия мала, начальное приближение может быть недостаточно точным, что приведет итерационную процедуру к нахождению локального минимума функции (2.8) вместо глобального. Эту проблему можно частично решить, если перебирать начальные догадки из некоторого достаточно широкого интервала с заданным шагом и проводить процедуру реконструкции для каждой начальной догадки, после чего выбрать результат, соответствующий наименьшему значению функции (2.8).
Следует отметить, что в случае длинного временного ряда, содержащего большое число периодов внешнего воздействия, даже малая погрешность в определении и неизбежно ведет к «набегу фазы» и плохому описанию внешнего воздействия соответствующими слагаемыми в модели (2.6). При реконструкции это сказывается таким образом, что область сходимости глобального минимума функции (2.8) оказывается очень узка. Поэтому реконструкцию следует проводить для достаточно коротких по сравнению с характерным периодом воздействия временных рядов (например, 5-10 периодов). При этом, однако, теряется преимущество, которое дает более длинный ряд с точки зрения точности и достоверности оценок. Эти потери можно частично скомпенсировать, если разбить исходный длинный временной ряд на короткие участки, провести реконструкцию для каждого из них, а затем усреднить полученные значения коэффициентов. Аналогичный подход был использован в работе [39] при реконструкции неавтономных нелинейных осцилляторов без запаздывания по временным рядам. Недостатком такого подхода можно считать довольно посредственные статистические свойства оценок коэффициентов (по сравнению тем, как если бы они были получены по всему ряду целиком).
Альтернативный подход состоит в том, чтобы сначала грубо оценить частоту внешнего воздействия по коротким временным рядам длиною в 5-10 характерных периодов внешнего воздействия, перебирая начальные приближения для и с большйм шагом из широкого интервала. Далее реконструкция проводится уже для всего имеющегося временного ряда, причем полученное на предыдущем этапе значение и используется в качестве начального приближения.
С увеличением степени полиномов N и М средний квадрат одношаговой ошибки прогноза (2.8) моделей уменьшается. Мы ограничиваемся такими значениями N и М, при которых величина (2.8) достигает насыщения и при дальнейшем их увеличении меняется слабо.
2.3. Использование априорной информации
Для многих практически важных систем [16, 19] уравнения вида (2.1) и даже (2.3) являются слишком общими, так что на основе априорной информации можно конкретизировать вид функции /. Например, для большинства известных генераторов с запаздывающей обратной связью можно функцию / представить в виде суперпозиции члена, ответственного за мгновенную динамику (как правило, это линейное слагаемое —х(Ь)), и нелинейной функции ф(х (Ь — т)), отвечающей за задержанную динамику:
/(Ь) — — х(Ь)+ф(х (Ь — т)). (2.9)
Учет априорной информации о природе системы крайне желателен, поскольку, как неоднократно подчеркивалось [3], использование универсальных моделей (в том числе с полиномами общего вида в качестве аппроксимирующих функций) приводит к необходимости оценки очень большого числа коэффициентов, что сложно реализовать на практике.
При аппроксимации внешнего воздействия тригонометрическим полиномом (2.7), уравнение (2.6) для систем с нелинейной функцией вида (2.9) может быть представлено следующим образом:
N М
ХП+1 — ахп + У] ЪихП-т + ^ [ск со8(кииАЬ) + йк вт(кипАг)). (2.10)
к=0 к=1
Значение т, при котором наблюдается минимальное значение Б2, принимается в качестве оценки дискретного времени запаздывания, которое связано с истинным временем запаздывания соотношением то — тАЬ. Определив коэффициенты модели (2.10), можно восстановить параметр инерционности и нелинейную функцию системы с запаздыванием (2.9), а также амплитуду внешнего воздействия х(Ь). В качестве оценки параметра е будем использовать величину ' — Аг/(1 — а). Нелинейную функцию / можно восстановить по коэффициентам Ък аппроксимирующего ее полинома, причем коэффициенты Ък следует умножить на е/АЬ. Амплитуду А внешнего воздействия можно оценить по коэффициентам Сд. и (4- тригонометрического полинома: А = \ZSfcli (с| + Частота и внешнего воздействия находится с помощью метода Ньютона описанным выше способом.
3. Тестирование метода
3.1. Восстановление неавтономной системы с запаздыванием с квадратичной нелинейностью
Рассмотрим сначала применение метода к временным рядам системы с запаздыванием (2.3), функция / которой является квадратичной, а внешнее воздействие гармоническим:
ех(Ь) — —х(г) + Л — х2(г — т0) + А вт(иЬ), (3.1)
где Л — управляющий параметр. Параметры уравнения (3.1) выберем следующие: Л — 1.85, то — 100, е — 20, А — 0.7, и — 0.02п. При этих параметрах неавтономная система (3.1) демонстрирует хаотические колебания (рис. 1а). Для построения модели будем использовать временной ряд длиной Ь — 20 000 точек при интервале выборки АЬ — 1. При построении модельного уравнения (2.6) ограничимся степенным полиномом второй степени N — 2 и тригонометрическим полиномом первой степени М — 1.
(а) (Ь)
'п—Г
(с)
Рис. 1. (а) Временной ряд системы с запаздыванием (3.1) в режиме хаотических колебаний. (Ъ) Зависимость среднего квадрата одношаговой ошибки прогноза модели от пробного времени запаздывания при N — 2 и М — 1. (с) Восстановленная нелинейная функция.
На рисунке 1Ь приведена зависимость погрешности (2.8) модели (2.6) от пробного времени задержки т , перебираемого из интервала от 0 до 500 с шагом 1, при N — 2 и М — 1. График демонстрирует отчетливый минимум при т — То — 100.
С помощью метода Ньютона найдем частоту воздействия /о — и/2п. Восстановленная частота воздействия /о совпадает с истинной частотой: /о — /о — 0.01. Коэффициенты модели (2.6) принимают при этом следующие значения: а — 0.95, Ьо — 0.0925, Ь1 ~ 10-16, Ь2 — -0.05, в1 ~ 10-17, dl — 0.035. По ним легко получить оценку значений е, Л и А: Л — 20, Л — Ьо/(1 — а1) — 1.85, А — 0.7. Заметим, что все параметры системы (3.1) восстановлены точно. Восстановленная нелинейная функция приведена на рисунке 1с. Она практически совпадает с истинной квадратичной функцией системы (3.1).
Мы исследовали эффективность метода в присутствии динамического и измерительного шума. Было показано, что при добавлении в динамику системы нормального белого шума (с нулевым средним и среднеквадратичным отклонением 20 % от среднеквадратичного отклонения сигнала системы без шума) положение минимума Б2(т) позволяет точно восстановить То, а метод Ньютона — точно восстановить /о. Остальные параметры системы удается восстановить с очень хорошей точностью: Л — 20.01, Л — 1.851, А — 0.701. С увеличением уровня шума минимум Б2(т) становится все менее выраженным, а потом исчезает.
Аналогично, метод оказывается работоспособен при добавлении измерительного шума с уровнем порядка 5 %—20 %. Измерительный и динамический шумы одинаковой амплитуды вносят примерно одинаковые искажения в оценки коэффициентов, что обусловлено тем, что мы восстанавливали отображение последования (разностную схему), а не дифференциальное уравнение, и, таким образом, избежали необходимости численно оценивать производную наблюдаемого временного ряда (именно измерительный шум вносит наибольшие искажения в оценку производной).
Работоспособность метода исследована при различных частотах внешнего периодического воздействия. Установлено, что метод остается эффективным в широком диапазоне
значений и. Нижняя граница допустимой частоты воздействия определяется длиной ряда, который должен содержать несколько периодов внешнего воздействия. Верхняя граница определяется для рассмотренных систем частотой дискретизации временного ряда: необходимо иметь хотя бы несколько точек на период воздействия; при этом работоспособность метода не зависит от соотношения собственной частоты колебаний, определяемой временем запаздывания, и частоты воздействия.
Возьмем теперь в качестве внешнего воздействия периодический сигнала осциллятора ван дер Поля, уравнение которого имеет вид
¿>(t) - (§ - и2(t)) v(t) + П2и(t) = 0, (3.2)
где § — параметр нелинейности, а Q — частота линейных колебаний. Пусть § = 1, Q = 0.02п. При этих значениях параметров система (3.2) демонстрирует сильно несинусоидальные автоколебания с частотой fo ^ 0.0022 и амплитудой V ~ 2. Отметим, что разница между частотой автоколебаний fo и параметром Q/2n = 0.01, определяющим частоту линейных колебаний, обусловлена нелинейностью системы.
Уравнение неавтономной системы с запаздыванием принимает при этом вид
£x(t) = -x(t) + Л - x2(t - то) + Kv(t), (3.3)
где K — коэффициент связи, определяющий амплитуду внешнего воздействия A = KV. Пусть K = 0.1, при этом A = 0.2.
Параметры уравнения (3.3) оставим такими же, как в уравнении (3.1): Л = 1.85, то = = 100, £ = 20. Поскольку сигнал осциллятора ван дер Поля является более сложным, чем гармонический сигнал, для его аппроксимации требуется тригонометрический полином более высокого порядка. При построении модельного уравнения (2.6) мы использовали N = 2 и M = 9.
Зависимость S2(т), как и для случая гармонического воздействия, демонстрирует минимум при т = 100. Восстановленная частота воздействия fo = 0.002224. Остальные параметры системы удается восстановить с высокой точностью: е = 19.985, Л = 1.8499, A = = 0.230. Восстановленная нелинейная функция практически совпадает с функцией, приведенной на рисунке 1c.
Присутствие шума оказывает качественно такое же влияние на результаты реконструкции системы (3.3), как и на результаты реконструкции системы (3.1), находящейся под гармоническим внешним воздействием.
Большинство известных методов реконструкции систем с запаздыванием ориентировано на использование хаотических временных рядов. Предложенный в данной работе метод может быть применен для восстановления неавтономных систем с запаздыванием, совершающих как хаотические, так и периодические колебания. При применении метода к системам (3.1) и (3.3), находящимся в режиме периодических колебаний, нами были получены результаты, качественно подобные представленным выше для хаотических режимов.
3.2. Восстановление неавтономного уравнения Икеды
Далее метод применялся к хаотическим временным рядам уравнения Икеды
x(t) = -x(t) + ¡j, sin (x(t - то) - xo), (3.4)
находящегося под внешним периодическим воздействием в отсутствие и в присутствии динамического и измерительного шума. Уравнение (3.4) относится к системам с запаздыванием, описываемым уравнением (2.3) с е = 1. Как и в разделе 3.1, в качестве внешнего воздействия использовались гармонический сигнал и сигнал осциллятора ван дер Поля.
Поскольку число рассмотренных ситуаций достаточно велико, а полученные для них результаты качественно похожи, ограничимся иллюстрацией лишь одного случая. Применим метод к хаотическому временному ряду уравнения Икеды, находящегося под внешним гармоническим воздействием в присутствии динамического шума,
x(t) = —x(t) + j sin (x(t — т0) — x0) + A sin(ut) + £(t) (3.5)
при j = 20, то = 2, Хо = n/3, A = 1, и = 0.2п. На рисунке 2 a приведен фрагмент временного ряда для случая, когда £(t) — это 20 %-ный гауссовский белый шум с нулевым средним значением. Весь ряд состоял из 20 000 точек при интервале выборки At = 0.01.
При построении модельного уравнения (2.10) мы использовали полиномы различной степени. Ограничимся степенным полиномом с N = 15 и гармоническим полином первой степени M = 1, поскольку при дальнейшем увеличении N и M величина Б2(т) меняется слабо. Зависимость S2(t), построенная при т, изменяющемся от 0 до 400 с шагом 1, имеет минимум при т = m = 20 (рис. 2b), обеспечивая точное восстановление то = mAt = 2.
0.0012
S2(t)
0.0006
2000
~ с—' ^ V/ 1 -
1
100
200
г
300
400
(Ь)
fiVn-r)
Рис. 2. (а) Хаотический зашумленный временной ряд неавтономного уравнения Икеды (3.5). (Ъ) Зависимость одношаговой ошибки прогноза модели от пробного времени запаздывания при N =15 и М =1. (с) Восстановленная нелинейная функция.
Восстановленная с помощью метода Ньютона частота воздействия fo = 0.1003 близка к истинной частоте fo = и/2п = 0.1. Коэффициент a = 0.9898 модели (2.10) дает оценку параметра инерционности е = 0.98, а коэффициенты в\ = —0.0029 и di = 0.0094 дают оценку амплитуды воздействия A = 0.96.
Восстановленная нелинейная функция приведена на рисунке 2c. Ее амплитуда позволяет получить оценку J1 = 20.059. Параметр xo может быть найден по формуле Хо =
= — arcsin (ex(t) + x(t))^j, где ex(t) + x(t) равно значению восстановленной функции /
при x(t — т0) = 0. Получаем Х0 = 1.063 (п/3 w 1.047).
3.3. Восстановление неавтономного уравнения Маккея — Гласса
В качестве последнего примера рассмотрим реконструкцию уравнения Маккея - Гласса
0.2х(Ь — т0) ,
х = -0.1хН) +-^-—, 3.6
1 + - го)
которое делением на 0.1 может быть приведено к виду (2.3) с е =10. Как и в предыдущих разделах, были исследованы случаи гармонического внешнего воздействия и воздействия сигналом осциллятора ван дер Поля как в отсутствие, так и в присутствии шума.
Так как полученные результаты опять качественно похожи, ограничимся иллюстрацией лишь одного случая. Применим метод к хаотическому временному ряду уравнения Маккея-Гласса, находящегося под воздействием сигнала осциллятора ван дер Поля (3.2) в присутствии измерительного шума:
0.2х{1 - т0) 1 + ж10С£-т0)
х = —0.1х(1) + ; + Кту{1). (3.7)
Пусть параметры осциллятора ван дер Поля (3.2) определены как 5 = 1, О = 0.05^. При этих значениях параметров система (3.2) демонстрирует сильно несинусоидальные автоколебания с частотой /о ~ 0.0115 и амплитудой V ~ 2. Возьмем то = 300, К =1, при этом А = 2.
Добавим к временному ряду уравнения (3.7) 5 %-ный гауссовский белый шум с нулевым средним значением. На рисунке 3а приведен фрагмент полученного временного ряда. Весь ряд состоял из 20 000 точек при интервале выборки АЬ = 1.
Рис. 3. (а) Хаотический зашумленный временной ряд неавтономного уравнения Маккея - Гласса (3.7). (Ь) Зависимость одношаговой ошибки прогноза модели от пробного времени запаздывания при N = 5 и М = 9. (с) Восстановленная нелинейная функция.
При построении модельного уравнения (2.10) мы использовали N = 5 и М = 9. Зависимость Б2 (т), построенная при т, изменяющемся от 0 до 500 с шагом 1, имеет минимум при т = то = 300 (рис. 3Ь). Восстановленная частота воздействия /о = 0.0115. Остальные параметры системы удается восстановить с хорошей точностью: е~ = 10.641 (истинное е = 10),
A = 2.307. Восстановленная нелинейная функция (рис. 3c) хорошо совпадает с истинной нелинейной функцией уравнения (3.7).
4. Заключение
Предложен подход к реконструкции неавтономных систем с запаздыванием, описываемых дифференциальным уравнением с задержкой. Подход основывается на реконструкции по экспериментальным данным отображения последования, представляющего собой разностную схему для исходных уравнений; при этом априорные знания о структуре системы используются для определения ее порядка, способа введения внешнего воздействия и способа аппроксимации нелинейной функции. Оценка коэффициентов модели, в том числе частоты воздействия и, производится нелинейным методом Ньютона, для достижения глобального минимума используется перебор начальных догадок для и, а также варьирование длины используемого ряда.
Метод позволяет восстановить управляющие параметры и нелинейную функцию системы с запаздыванием, а также частоту и амплитуду внешнего периодического воздействия. Метод может быть применен для реконструкции неавтономных систем с запаздыванием, совершающих как хаотические, так и периодические колебания, и остается эффективным при высоких уровнях шума. Следует отметить, что возможность восстановления систем в периодических режимах является существенным плюсом метода по сравнению с методами, разработанными для обыкновенных дифференциальных уравнений [8], где для успеха реконструкции временной ряд обязательно должен быть хаотическим. Работоспособность метода показана на различных модельных системах с запаздыванием при различных видах периодического внешнего воздействия.
Несмотря на то, что в работе метод тестировался на системах первого порядка с аддитивным воздействием, он может быть использован и для систем с другим способом введения воздействия (например, аналогичные результаты для систем в виде обыкновенных дифференциальных уравнений были получены в [7]) и для систем более высокого порядка. В этих случаях также необходимо учитывать априорную информацию о структуре модельного уравнения системы, вводя соответствующие члены в восстанавливаемую модель и стараться максимально сократить сложность нелинейной функции. При реконструкции системы предложенный подход опирается на уже известные методы восстановления времени запаздывания [29, 33].
Список литературы
[1] Gouesbet G., Meunier-Guttin-Cluzel G., Menard O. Chaos and its reconstruction. New York: Nova Science, 2003. 320 pp.
[2] Eykhoff P. System identification: Parameter and state estimation. London-New York: Wiley-Interscience, 1974. 555 pp.
[3] Безручко Б. П., Смирнов Д. А. Математическое моделирование и хаотические временные ряды. Саратов: ГосУНЦ «Колледж», 2005. 320 с.
[4] Horbelt W., Timmer J., Bünner M. J., Meucci R., Ciofini M. Identifying physical properties of a CO2 laser by dynamical modeling of measured time series // Phys. Rev. E, 2001, vol.64, no. 1, 016222, 7 pp.
[5] Павлов А. Н., Янсон Н. Б., Капитаник Т., Анищенко В. С. Реконструкция динамических систем по сигналам малой длительности // Письма в ЖТФ, 1999, т. 25, № 11, с. 7-13.
[6] Анищенко В. С., Павлов А. Н., Янсон Н. Б. Реконструкция динамических систем в приложении к решению задачи защиты информации // ЖТФ, 1998, т. 68, № 12, с. 1-8.
[7] Besruchko B.P., Smirnov D.A. Constructing nonautonomous differential equations from experimental time series // Phys. Rev. E., 2000, vol. 63, no. 1, 016207, 7 pp.
[8] Безручко Б. П., Смирнов Д. А., Сысоев И. В., Селезнев Е. П. Реконструкция моделей неавтономных систем с дискретным спектром воздействия // Письма в ЖТФ, 2003, т. 29, №19 с. 69-76.
[9] Куликов М.Ю., Мухин Д. Н., Фейгин А. М. Байесова стратегия оценки точности характеристик, извлекаемых из экспериментальных данных с помощью базовых динамических моделей атмосферных фотохимических систем // Изв. вузов. Радиофизика, 2009, т. 52, №9, с. 690-699.
[10] Яхно Ю.В., Мольков Я. И., Мухин Д. Н., Лоскутов Е. М., Фейгин А. М. Реконструкция оператора эволюции как способ анализа электрической активности мозга при эпилепсии // Изв. вузов. ПНД, 2011, т. 19, №6, с. 156-172.
[11] Baake E., Baake M., Bock H. G., Briggs K.M. Fitting ordinary differential equations to chaotic data // Phys. Rev. A, 1992, vol.45, no. 8, pp. 5524-5529.
[12] Bezruchko B. P., Smirnov D. A., Sysoev I. V. Identification of chaotic systems with hidden variables (modified Bock's algorithm) // Chaos Solitons Fractals, 2006, vol. 29, no. 1, pp. 82-90.
[13] Kuang Y. Delay differential equations with applications in population dynamics. Boston: Academic Press, 1993. 398 pp.
[14] Bocharov G. A., Rihan F. A. Numerical modelling in biosciences using delay differential equations // J. Comput. Appl. Math., 2000, vol. 125, nos. 1-2, pp. 183-199.
[15] Mokhov 1.1., Smirnov D.A. El Niño - Southern Oscillation drives North Atlantic Oscillation as revealed with nonlinear techniques from climatic indices // Geophys. Res. Lett., 2006, vol. 33, no. 3, L03708, 4 pp.
[16] Mackey M.C., Glass L. Oscillations and chaos in physiological control systems // Science, 1977, vol. 197, pp. 287-289.
[17] Гласс Л., Мэки М. От часов к хаосу: Ритмы жизни. Москва: Мир, 1991. 248 с.
[18] Ringwood V., Malpas S.C. Slow oscillations in blood pressure via a nonlinear feedback model // Am. J. Physiol. Regul. Integr. Comp. Physiol., 2001, vol.280, R1105-R1115.
[19] Ikeda K. Multiple-valued stationary state and its instability of the transmitted light by a ring cavity system // Opt. Commun., 1979, vol.30, no. 2, pp. 257-261.
[20] Дмитриев А. С., Кислов В. Я. Стохастические колебания в радиофизике и электронике. Москва: Наука, 1989. 280 с.
[21] Prokhorov M. D., Ponomarenko V. I. Encryption and decryption of information in chaotic communication systems governed by delay-differential equations // Chaos Solitons Fractals, 2008, vol. 35, no. 5, pp. 871-877.
[22] Пиковский А. С., Розенблюм М.Г., Куртс Ю. Синхронизация: Фундаментальное нелинейное явление. Москва: Техносфера, 2003. 508 с.
[23] 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, no. 5, pp. 336-344.
[24] Tian Y.-C., Gao F. Extraction of delay information from chaotic time series based on information entropy // Phys. D, 1997, vol. 108, nos. 1-2, pp. 113-118.
[25] Ellner S.P., Kendall B.E., Wood S.N., McCauley E., Briggs C.J. Inferring mechanism from time-series data: Delay differential equations // Phys. D, 1997, vol. 110, nos. 3-4, pp. 182-194.
[26] Hegger R., Bünner M.J., Kantz H., Giaquinta A. Identifying and modeling delay feedback systems // Phys. Rev. Lett., 1998, vol.81, no.3, pp. 558-561.
[27] Buünner M. J., Ciofini M., Giaquinta A., Hegger R., Kantz H., Meucci R., Politi A. Reconstruction of systems with delayed feedback: 1. Theory // Eur. Phys. J. D, 2000, vol. 10, no. 2, pp. 165-176.
[28] Horbelt W., Timmer J., Voss H.U. Parameter estimation in nonlinear delayed feedback systems from noisy data // Phys. Lett. A, 2002, vol.299, nos. 5-6, pp. 513-521.
[29] Пономаренко В. И., Прохоров М.Д., Караваев А. С., Безручко Б. П. Определение параметров систем с запаздывающей обратной связью по хаотическим временным реализациям // ЖЭТФ, 2005, т. 127, №3, с. 515-527.
[30] Ortin S., Gutierrez J. M., Pesquera L., Vasquez H. Nonlinear dynamics extraction for time-delay systems using modular neural networks synchronization and prediction // Phys. A, 2005, vol.351, no. 1, pp. 133-141.
[31] Siefert M. Practical criterion for delay estimation using random perturbations // Phys. Rev. E, 2007, vol. 76, no. 2, 026215, 5 pp.
[32] Yu D., Frasca M., Liu F. Control-based method to identify underlying delays of a nonlinear dynamical system // Phys. Rev. E, 2008, vol.78, no. 4, 046209, 9 pp.
[33] Prokhorov M.D., Ponomarenko V.I. Reconstruction of time-delay systems using small impulsive disturbances // Phys. Rev. E, 2009, vol.80, no. 6, 066206, 6pp.
[34] 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, no.4, 046212, 9 pp.
[35] Ma H., Xu B., Lin W., Feng J. Adaptive identification of time delays in nonlinear dynamical models // Phys. Rev. E, 2010, vol.82, no. 6, 066210, 11pp.
[36] Dai Ch., Chen W., Li L., Zhu Y., Yang Y. Seeker optimization algorithm for parameter estimation of time-delay chaotic systems // Phys. Rev. E, 2011, vol.83, no.3, 036203, 11pp.
[37] Shigaev A. M., Dmitriev B. S., Zharkov Y. D., Ryskin N. M. Chaotic dynamics of delayed feedback klystron oscillator and its control by external signal // IEEE Trans. Electron Devices, 2005, vol. 52, no. 5, pp. 790-797.
[38] Marchewka C., Larsen P., Bhattacharjee S., Booske J., Sengele S., Ryskin N.M., Titov V.N. Generation of chaotic radiation in a driven traveling wave tube amplifier with time-delayed feedback // Phys. Plasmas, 2006, vol. 13, no. 1, 013104, 20 pp.
[39] Безручко Б. П., Селезнев Е. П., Смирнов Д. А. Реконструкция уравнений неавтономного нелинейного осциллятора по временному ряду: Модели, эксперимент // Изв. вузов. ПНД, 1999, т. 7, №1, с. 49-68.
Reconstruction of time-delay systems under external periodic driving
Marina V. Sysoeva1, Vladimir I. Ponomarenko2, Mikhail D. Prokhorov3, Ilya V. Sysoev4
1,2'4Saratov State University named after N. G. Chernyshevsky Astrahanskaya st. 83, Saratov, 410012, Russia
2'3Saratov Branch, Kotel'nikov Institute of Radio-Engineering and Electronics, Russian Academy of Sciences Zelenaya 38, Saratov, 410019, Russia
1bobrichek@mail.ru, 2ponomarenkovi@gmail.com, 3mdprokhorov@yandex.ru, 4ivssci@gmail.com
A method is proposed for the reconstruction of first-order time-delay systems under external periodic driving from their time series. The method takes into account the structure of the model equation of the system, while constructing the autoregressive model. The proposed method allows one to reconstruct the delay time, the parameter characterizing the system inertial properties, the nonlinear function, and the amplitude and frequency of the external periodic driving. The method efficiency is demonstrated in a numerical experiment by reconstructing a number of different nonautonomous time-delay systems.
MSC 2010: 37M10
Keywords: reconstruction of model equations, time-delay systems, time series analysis Received May 21, 2013, accepted October 9, 2013
Citation: Rus. J. Nonlin. Dyn., 2013, vol. 9, no. 4, pp. 613-625 (Russian)
