«РЕКОНСТРУКЦИЯ УРАВНЕНИЙ КОЛЕБАТЕЛЬНЫХ СИСТЕМ ПРИ НАЛИЧИИ СКРЫТЫХ ПЕРЕМЕННЫХ И ВНЕШНИХ ВОЗДЕЙСТВИЙ ...»

Саратовский государственный университет

им. Н.Г. Чернышевского

Саратовское отделение Института радиотехники и электроники

Российской Академии Наук

На правах рукописи

УДК: 530.182, 53.083

СЫСОЕВ Илья Вячеславович

РЕКОНСТРУКЦИЯ УРАВНЕНИЙ КОЛЕБАТЕЛЬНЫХ

СИСТЕМ ПРИ НАЛИЧИИ СКРЫТЫХ ПЕРЕМЕННЫХ И

ВНЕШНИХ ВОЗДЕЙСТВИЙ

01.04.03 Радиофизика Диссертация на соискание ученой степени кандидата физико-математических наук

Научный руководитель:

доктор физико-математических наук, профессор Безручко Б.П.

Саратов Оглавление Введение 1 Реконструкция при наличии скрытых переменных 1.1 Введение..................................... 1.2 Методы оценки параметров при наличии скрытых переменных...... 1.2.1 Метод начального условия....................... 1.2.2 Метод множественной стрельбы.................... 1.3 Сравнительный анализ в численном эксперименте............. 1.3.1 Методика сравнения........................... 1.3.2 Оценка параметров системы Лоренца................. 1.3.3 Оценка влияния измерительного шума на результат реконструкции 1.3.4 Оптимальный выбор параметров алгоритма: количества фрагментов и их длины.............................. 1.3.5 Преимущества модифицированного метода............. 1.3.6 Выбор оптимального количества разрывов.............. 1.3.7 Универсальность оптимального выбора числа подсегментов L при различном числе сегментов непрерывности............ 1.3.8 Реконструкция системы Рёсслера................... 1.3.9 Подбор стартовых догадок для скрытых переменных........ 1.4 Выводы...................................... 2 Реконструкция систем под регулярным воздействием 2.1 Введение..................................... 2.2 Методика реконструкции............................ 2.3 Численный эксперимент............................ 2.3.1 Реконструкция при гладком периодическом воздействии...... 2.3.2 Влияние шума на результат реконструкции............. 2.3.3 Реконструкция при треугольном периодическом воздействии... 2.3.4 Реконструкция при воздействии c субгармониками......... 2.3.5 Реконструкция при квазипериодическом воздействии........ 2.4 Выводы...................................... 3 Восстановление внешнего воздействия методами работы со скрытыми переменными 3.1 Введение..................................... 3.2 Методика реконструкции............................ 3.3 Численные примеры реконструкции...................... 3.3.1 Автономный режим периодический, воздействие хаотическое... 3.3.2 Реконструкция по зашумлённым данным............... 3.3.3 Автономный режим хаотический, воздействие хаотическое.... 3.3.4 Автономный режим хаотический, воздействие шумом....... 3.3.5 Ситуация большого числа скрытых переменных.......... 3.3.6 Реконструкция уравнений генератора с 1,5 степенями свободы.. 3.4 Возможные приложения метода........................ 3.4.1 Скрытая передача и кодирование информации........... 3.5 Выводы...................................... 4 Приложение методов реконструкции 4.1 Введение..................................... 4.2 Способ определения характеристик нелинейных устройств........ 4.2.1 Выбор объекта и постановка эксперимента.............. 4.2.2 Выбор эквивалентного представления................. 4.2.3 Реконструкция в режиме малых сигналов.............. 4.2.4 Реконструкция в режиме больших периодических сигналов. Зависимость ёмкости от частоты воздействия............... 4.2.5 Реконструкция в режиме больших периодических сигналов. Зависимость ёмкости от амплитуды воздействия............. 4.2.6 Использование метода множественной стрельбы при реконструкции характеристик диода........................ 4.2.7 Учёт сопротивления базы........................ 4.4.2 Реконструкция модели нефрона по модельным и экспериментальным реализациям............................ Введение Развитие концепции динамического хаоса продемонстрировало возможность описания сложных движений с помощью простых нелинейных моделей и усилило интерес к методам моделирования по дискретным последовательностям экспериментальных данных (временным рядам).

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

Исторически в основе реконструкции1 лежит задача аппроксимации точек на плоскости (x, y) функцией y = f (x) [2, 3, 4, 5]. Но в настоящее время речь идёт об описании сложных, часто хаотических процессов, поэтому появляется необходимость построения по экспериментальным данным разностных уравнений (дискретных отображений отображений последования):

Реконструкция термин, используемый в нелинейной динамике. В математической статистике принят другой термин идентификация систем [1].

Рис. 1: два подхода к обработке экспериментальных временных рядов непосредственная обработка и обработка на основе построения модели.

обыкновенных дифференциальных уравнений (ОДУ):

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

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

• прогноз дальнейшего поведения с целью предсказания последствий или контроля и управления объектом;

• определение наличия и направленности (или даже структуры связей) между двумя системами или между подсистемами одной системы;

• классификация систем [10].

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

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

• скрытая (конфиденциальная) передача, кодирование и декодирование информации [12, 13];

Восстановленные по временным рядам модели используются в настоящее время в радиофизике [14], лазерной физике [15, 16], биофизике и биологии (в частности, при изучении структуры и механизмов функционирования клетки [17]), метеорологии [18, 19], сейсмографии [20], экономике [21, 22], медицине и физиологии [23, 24, 25], астрофизике [26], и т.д.

Развитие алгоритмов реконструкции долгое время шло по пути создания универсальных методик, рассчитанных на широкий класс систем обладающих высокою степенью универсальности. К таковым относится, например, так называемый стандартный подход [27], когда по скалярной наблюдаемой восстанавливаются уравнения в виде:

где yi, i = 1,..., D динамические переменные, D размерность модели, f (y1, y2,..., yD ) единственная неизвестная функция, представляемая обычно в виде полинома. Претензии таких моделей на общность теоретически обоснованы, поскольку в виде (3) может быть представлена почти любая система ОДУ (подробнее см. [28]) и любая непрерывная функция f в (3) может быть сколь угодно точно равномерно приближена алгебраическим полиномом (теорема Вейерштрасса). К тому же они требуют минимума дополнительной информации об объекте. Однако при их практическом применении возникают многочисленные и трудноразрешимые проблемы, главная причина которых проклятие размерности. Так принято называть проблемы, связанные повышением размерности модели: быстрый рост количества коэффициентов, недостаточная длина экспериментального ряда для восстановления пространства состояний, а при использовании последовательного дифференцирования, как в 3, также увеличение влияния шумов. Большое число коэффициентов разложения в ряд по выбранному базису приводит к росту вычислительных затрат, а главное существенно снижает область сходимости используемых итерационных алгоритмов, так как подавляющее большинство их являются лишними. Существующие способы устранения лишних коэффициентов достаточно громоздки, так как связаны с многократным перебором, и всё равно позволяют решить эту задачу лишь частично. Важно также, что коэффициенты 3 не имеют понятного физического смысла и их сложно интерпретировать.

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

Одною из распространённый ситуаций, требующих разработки специализированных алгоритмов необходимо, является случай так называемых скрытых переменных. Скрытыми называются такие (D l) из общего числа D переменных модели (2), которые не могут быть измерены в принципе или полученные временные ряды этих переменных непригодны для моделирования, например, из-за высокой зашумлённости. Основная причина несовершенство аппаратуры, а также риск разрушить объект (это особенно актуально для биологических систем). Такая ситуация часто встречается на практике, если модель записывается из первых принципов и стоит задача проверки её адекватности или оценки неизвестных параметров и нелинейных функций.

Основной принцип методов работы со скрытыми переменными был выдвинут достаточно давно и заключается в том, что начальные условия для скрытых переменных включаются в число неизвестных параметров3. Для этих начальных условий задаются стартовые догадки yl+1 (t1 ),..., yD (t1 )4, задаются стартовые догадки для параметров cs,...cs. Затем формулируется критерий, на основе которого стремятся достигнуть максимально возможАльтернативный подход, использующий явление хаотической синхронизации см. в [29].

Индексом s сверху будем далее обозначать стартовые догадки для соответствующих величин, индексом 0 начальные условия.

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

Применение описанного подхода напрямую5 для длинных хаотических рядов малоэффективно, поскольку высокая чувствительность траектории к начальным условиям приводит к тому, что задачу глобальной оптимизации сложно решить из-за большого количества локальных минимумов вблизи глобального (см. [30]). Поэтому в [31] предложена модификация, состоящая в том, что для каждой переменной задаётся не одно, а сразу несколько L начальных условий yk (t1 ), yk (tn+1),..., yk t(L1)n+1 на разных участках временного ряда. Такой модифицированный метод называется методом множественной стрельбы по аналогии с одноимённым методом решения краевых задач для ОДУ. В случае, если дополнительно на траекторию накладывается условие непрерывности (4) траектории ( сшивания фрагментов), алгоритм часто называют методом Бока.

где n число точек в одном фрагменте (таким образом N = Ln), y tjn+1, y0 tj(n1)+1 значение вектора состояния модели, полученного интегрированием уравнений с начальными условиями y0 tj(n1)+1.

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

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

Такой подход часто называют методом начального условия.

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

Другой практически важный класс объектов, для которых в диссертации разрабатываются специализированные методики, неавтономные системы. Необходимость этого обусловлена тем, что, хотя неавтономную динамическую систему можно свести к виду (3), это ведёт росту размерности и сложности аппроксимирующей функции f (y1, y2,..., yD ). Поэтому неавтономность следует учитывать в структуре модели. Информация об этом может быть известна априорно, либо такое предположение можно сделать, анализируя сами временные ряды, например, по спектру, если в нём присутствуют один или несколько резких пиков.

Существуют несколько подходов к реконструкции неавтономных систем. Один из них предполагает возможность измерения рядов той же системы в автономном режиме, на основе которых восстанавливаются параметры модели, а затем уже по невтономным реализациям реконструируется само воздействие [32, 33, 34]. Преимущество такого подхода состоят в том, что на воздействие не накладываются никакие ограничения: оно может быть как регулярным, так и хаотическим или шумовым. Но доступность рядов автономной системы существенное и трудновыполнимое на практике требование.

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

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

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

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

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

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

Для достижение поставленной цели были решены следующие задачи:

• проведена оценка работоспособности современных методов работы со скрытыми переменными, основанных на алгоритме множественной стрельбы, для чего требуется ввести объективные количественные критерии;

• усовершенствован метод Бока на случай длинных рядов, а также разработаны рекомендации на оптимальную длину используемого ряда, его способ деления и вид стартовых догадок для скрытых переменных;

• осуществлено расширение известных алгоритмов реконструкции неавтономных систем на случай произвольного регулярного (сложного периодического и квазипериодического), хаотического и шумового воздействия;

• разработан способ определения нелинейных характеристик радиотехнических устройств в реальном режиме эксплуатации;

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

Объекты исследования. Представленные в работе методы, как заимствованные, так и оригинальные тестируются на эталонных динамических системах, таких как неавтономный осциллятор Тоды [37, 38], система Лоренца [39], система Рёсслера [40], а также ставших уже классическими моделях радиотехнических генераторов хаоса с 1,5 и 2,5 степенями свободы [41, 42], рекомендуемых авторами для использования в системах связи на хаотической несущей. Метод реконструкции нелинейных характеристик радиофизических устройств тестируется на примере полупроводникового диода с p-n переходом типа КД202Р, от которого измеряются ряды тока через него и напряжения на нём при различных видах воздействия. В биологических приложениях используются модель нефрона [43] и различные модели голосовых связок человека [44, 45, 46].

На защиту выносятся следующие положения и результаты:

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

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

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

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

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

Научная новизна работы состоит в следующем:

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

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

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

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

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

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

Теоретическая и практическая значимость результатов.

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

• Исследование эффективности и пределов применимости методов множественной стрельбы с помощью впервые введённых количественных критериев подтверждает ранее высказанное мнение об их высокой работоспособности и предлагает общие рекомендации по выбору некоторых параметров моделирования: длины используемого временного ряда, количества его сегментов, способу подбора стартовых догадок для скрытых переменных. Также показано, что мало популярный модифицированный метод множественной стрельбы, допускающий разрывы траектории модели, более эффективен, чем оригинальный алгоритм Бока, причём его преимущества проявляются для более длинных временных рядов. Это позволит в дальнейшем при моделировании сократить и упростить этап подбора параметров модели.

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

Построение системы скрытой передачи или кодирования информации.

• На основе метода реконструкции неавтономных систем разработан способ измерения нелинейных характеристик устройств, позволяющий измерять величины, недоступные прямому измерению, в произвольном режиме эксплуатации исследуемого устройства. Способ опробован на примере реконструкции вольтамперных и вольтфарадных характеристик полупроводникового диода с p-n переходом. Данный способ позволяет проводить разбраковку нелинейных устройств на основе произвольно выбранного критерия, что может быть использовано, в частности, при построении систем телекоммуникаций, где важно с высокой точностью соблюсти идентичность компонентов системы приёмник–передатчик.

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

Апробация результатов. Основные результаты диссертации были доложены на следующих конференциях:

• International Symposium Topical problems of nonlinear wave physics (Nizhny Novgorod, 2003).

• The second international conference on circuits and systems for communication (Moscow, 2004).

• 6th International School on Chaotic Oscillations and Pattern Formation (Saratov, 2001).

• VII Международная школа-конференция Хаотические автоколебания и образование структур (Саратов, 2004).

• XII и XIII Всероссийские школы-конференции Нелинейные волны 2004 и Нелинейные волны 2006 (Нижний Новгород).

• Международная научно-техническая конференция Радиотехника и связь (Саратов, 2005).

• VI и VII Всероссийские научные конференции Нелинейные колебания механических систем (Нижний Новгород, 2002, 2005).

• Межвузовская конференция Современные проблемы электроники и радиофизики СВЧ (Саратов, 2001).

• I Конференция молодых учёных Наноэлектроника, нанофотоника и нелинейная физика (Саратов, 2006).

• научные школы-конференции Нелинейные дни в Саратове для молодых (Саратов, 2000-2006).

Также результаты неоднократно обсуждались на научных семинарах:

• кафедры динамического моделирования и биомедицинской инженерии ФНиБМТ СГУ, • кафедры электроники, колебаний и волн ФНП СГУ (объединённые семинары с участием сотрудников кафедры нелинейной физики ФНП СГУ и сотрудников отделения физики нелинейных систем НИИ ЕН Исследования были поддержаны Российским фондом фундаментальных исследований (гранты №02-02-17578, №05-02-16305), Фондом некоммерческих программ Династия и Американским фондом гражданских исследований и развития для государств бывшего Советского Союза (CRDF REC-006).

По теме диссертации опубликованы 22 работы: 4 статьи в реферируемых журналах, входящих в перечень рекомендованных ВАК, 2 в тематических сборниках статей, 16 статей и тезисов в сборниках трудов конференций.

Структура и объём диссертации Диссертация состоит из введения, четырёх глав, заключения и списка литературы. В ней содержится 101 страница текста, 51 рисунок, библиография из 110 наименований. Общий объём диссертации 150 страниц.

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

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

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

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

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

Результаты диссертационной работы и выводы обобщаются и обсуждаются в заключении.

Глава Реконструкция при наличии скрытых переменных 1.1 Введение Задача построения математической модели, описывающей сложную динамику исследуемого объекта, по временным реализациям экспериментально наблюдаемых величин (t), может иметь различные постановки.

Мы рассматриваем такую, в которой структура модельных уравнений полностью известна из первых принципов или других соображений:

где y модельная переменная (D-мерный вектор состояния), а c Pмерный вектор параметров, который необходимо найти. Задача сводится к оценке неизвестных параметров c1,..., cP по наблюдаемому временному ряду. Такая постановка встречается, например, в аэродинамике (оценка коэффициентов аэродинамических сил при летном эксперименте с самолетом [47]), лазерной физике (оценка скоростей перехода между уровнями рабочего вещества газового лазера [15]), радиотехнике (моделирование электрических цепей с сегнетоэлектрическими или полупроводниковыми нелинейностями [11, 48]), клеточной биологии (оценка параметров модели нейрона [49], модели сигнального пути клеток определенного вида [17]).

Стандартная процедура реконструкции модели (1.1) предполагает знание векторного временного ряда значений y, то есть требуется иметь ряды всех D компонент, формируемых по данным эксперимента из ряда наблюдаемой {i}. Например, одна из переменных может просто совпадать с наблюдаемой, т.е. измеряется непосредственно, а другие получаются из неё дифференцированием, интегрированием, использованием временной задержки, весовым суммированием и другими способами. Однако нередко из-за невозможности измерить нужную величину вследствие зашумлённости или иных причин удаётся сформировать ряды лишь для l < D компонент y. (D l) переменных, ряды которых не могут быть получены из данных наблюдения, называют скрытыми. Наличие скрытых переменных сильно усложняет задачу, требуя использования специальных подходов.

В этой ситуации основная тяжесть задачи построения модели переносится на этап оценки параметров, причем в число искомых величин кроме параметров вводят начальные условия. Обычно для оценки параметров используется метод наименьших квадратов1, который формулируется в данной ситуации следующим образом. Нужно подобрать такие начальные условия модели y0 и параметры c, чтобы обеспечить минимальное отклонение решения модельных уравнений y от данных экспериментального наl блюдения y. Сравнение проводится только по l нескрытым компонентам y. То есть требуется минимизировать функцию:

где yi y0, c (1.1), y аналогичные значения, полученные из ряда наблюдаемой. Миi нимизацию (1.2) проводят с помощью численных итерационных методов, отправляясь от некоторых стартовых догадок для искомых величин y0, c.

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

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

чительном N весьма изрезан и имеет множество локальных минимумов, а область сходимости в глобальный минимум, соответствующий наиболее точным оценкам, очень узка, то есть требуются исключительно удачные стартовые догадки, чтобы туда попасть. Обойти трудности обещает изящный технический трюк алгоритм множественной стрельбы [31, 30]. Но, как показывает опыт, он требует дополнительных усилий для его успешной реализации. Эта проблема уже привлекала к себе внимание [15, 50], но систематический её анализ не проводился.

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

1.2 Методы оценки параметров при наличии скрытых 1.2.1 Метод начального условия Так называют метод оценки параметров, состоящий в минимизации непосредственно функции (1.2), где N длина всего наблюдаемого временного ряда [51]. На практике чем N больше, тем более достоверны получаемые оценки. Во-первых, без достаточного количества точек при наличии шумов достоверность полученных оценок оказывается низкою, во-вторых, важно, чтобы используемый ряд отражал все временные масштабы объекта. Но метод начального условия, как правило, неприменим для достаточно больших N, особенно в случае хаоса, поскольку из-за экспоненциальной чувствительности траектории модели к начальным условиям y0 область сходимости в глобальный минимум целевой функции (1.2) столь мала, что в неё почти невероятно попасть2. Недостатки этого подхода в сравнении с другими (см. п. 1.2.2) наглядно показаны в литературе, в частности, в [30] и [48], поэтому подробное рассмотрение его нецелесообразно.

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

Алгоритм Бока Увеличить допустимые значения N и погрешностей в стартовых догадках позволяет модификация метода начального условия3. Она направлена на то, чтобы хотя бы на промежуточных этапах минимизации целевой функции снизить чувствительность траектории модели к начальным условиям y0. Это достигается путем разбиения исходного ряда на L более коротких сегментов (n длина сегмента, N = Ln) и использования начальных условий на этих сегментах y1, yn+1,..., y(L1)n+1 в качестве дополнительных искомых величин аргументов целевой функции S:

где нижний индекс y(i1)n+1 соответствует первой точке i-того фрагмента. Чтобы избежать большого числа независимых неизвестных, увеличиЕсли временной ряд периодический, проблема разбегания траекторий, хотя и не по экспоненциальному закону, тоже может иметь место, например, в ситуации, когда какой-либо из параметров связан с периодом: тогда небольшая ошибка в его задании на больших временах даст существенное расхождение.

Иной подход к решению этой проблемы предлагают методы, основанные на хаотической синхронизации. Их принцип состоит в том, что наблюдаемая синхронизует модель, причём одновременно решаются эволюционные уравнения для параметров, которые таким образом подгоняются по экспериментальному временному ряду. Подробнее см. [29, 52, 53, 54].

вающего дисперсию оценок, вводится условие непрерывности траектории модели сшивания фрагментов:

Минимизация (1.3) при условии (1.4) это задача условной минимизации.

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

Но в процессе работы итерационного метода сегменты всё лучше удовлетворяют условию (1.4), так что в итоге траектория модели оказывается непрерывной. При L = 1, n = N алгоритм Бока превращается в метод начального условия.

Хотя в [30] утверждается, что для эффективности алгоритма Бока не требуется стартовых догадок, близких к истинным значениям, практика показывает, что это далеко не всегда так (см. ниже, пп. 1.3.5,1.3.6). Алгоритм только до некоторой степени расширяет возможности оценки параметров по хаотическому ряду по сравнению с наивным методом начального условия. Причина этого состоит в том, что итоговое условие (1.4) часто оказывается для хаотических систем слишком жёстким и при достаточно большой длине N удается попасть только в локальные минимумы целевой функции (1.3).

Кусочные методы Чтобы точнее оценить параметры проще всего разбить ряд на короткие сегменты (L штук) и проводить реконструкцию по каждому из сегментов отдельно, например, с помощью того же алгоритма Бока. Итоговая оценка получается как среднее арифметическое c = ci. Такой подход называL ют кусочным [50]; при использовании на отдельных сегментах алгоритма Бока будем говорить о кусочном методе множественной стрельбы.

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

Модифицированный метод Опираясь на статистические соображения [50], можно утверждать, что использование сразу всего объема данных позволяет получать оценки с меньшим смещением, чем при кусочных подходах. Поэтому мы предлагаем обратить внимание на модификацию алгоритма Бока (далее модифицированный метод), которая фактически уже использовалась в работе [15] для нехаотических рядов и упоминалась в [50]. Она состоит в отказе от непрерывности траектории модели в некоторые моменты времени внутри интервала наблюдения, т.е. условие (1.4) не требуется для ( 1) моментов времени. Это означает, что начальные условия модели в эти моментов, включая начальный, становятся независимыми искомыми величинами. Мы выберем их распределёнными равномерно внутри интервала наблюдения.

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

число сегментов с независимыми начальными условиями и L число подсегментов, на которые делится каждый сегмент для реализации алгоритма Бока, где N = Ln. При = 1 модифицированный метод сводится к алгоритму Бока. Этот подход, на наш взгляд незаслуженно, не получил до сих пор широкого распространения. Может быть, это связано с тем, что нарушение гладкости траектории на первый взгляд не сулит хорошей модели и настораживает.

1.3 Сравнительный анализ в численном эксперименте 1.3.1 Методика сравнения Чтобы охарактеризовать и сопоставить эффективность алгоритмов были введены специальные количественные меры. Во-первых, введены нормированные на истинные значения c0 параметры:

Это сделано потому, что отклонение на одну и ту же абсолютную величину для разных параметров имеет существенно разные последствия для режима системы в зависимости от абсолютной величины этого параметра. Точка b = 0 соответствует глобальному минимуму, если пренебречь шумами, включая погрешности вычислений.

Во-вторых, было использовано наглядное графическое представление:

на рис. 1.1 в сечении пространства параметров плоскостью bi2, bi3 в шкале серого отмечены результаты выбора соответствующих стартовых догадок.

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

Области хороших догадок в типичном случае сильно изрезаны, к тому же само по себе такое представление не даёт количественной меры успешности метода. В качестве таковой меры µ можно использовать относительную площадь светлой области внутри круга с центром в точке b = 0 и радиусом r. Фактически, µ можно интерпретировать как вероятность попадания в глобальный минимум, если ошибка в стартовых догадках не больше r.

Однако величина µ не всегда удобна в применении, поскольку сходимость метода характеризуется не одним числом, а функцией µ(r). Более Рис. 1.1: плоскость нормированных стартовых догадок для параметров системы Лоренца (сечение пространства параметров плоскостью b1 = 0 ). Белым обозначены точки, откуда достигается глобальный минимум, оттенками серого локальные минимумы (интенсивность зависит от их удалённости от глобального), чёрным алгоритм минимизации расходится: (a), (b) алгоритм Бока при L = 30, n = 35 (рис. (b) представляет собою фрагмент из центральной части рисунка (a)), (c) модифицированный метод с простой интегральный критерий эффективности получается, если рассмотреть обратную зависимость r(µ). Терпимость метода к ошибкам в стартовых догадках усреднённо характеризуется величиной радиуса rµ, при которой вероятность попадания в глобальный минимум не меньше заданного µ (см. рис. 1.1b,c). Для определенности далее будет в основном использоватся радиус 100%-ной сходимости4.

r 1.3.2 Оценка параметров системы Лоренца В качестве первого объекта для исследования эффективности методов (работоспособности при длинных рядах и терпимости к стартовым догадкам) была выбрана система Лоренца с параметрами c0 = 10, c0 = 8/3, c0 = 465 по решению, полученному при начальном условии y1 = 7.60, y2 = 12.37, y3 = 38.66, произвольно выбранном на хаотическом аттракторе (старший ляпуновский показатель 1 = 1.23 [30]). Уравнения интегрировались методом Рунге-Кутты 4-го порядка с шагом 0.001, интервал выборки составлял 0.002 (около 280 точек на характерном периоде T ). В качестве ряда скалярной наблюдаемой использовалась последовательность значений переменной y1. Для большей реалистичности к переменной аддитивно добавлялся нормальный шум: = y1 +. Две остальные переменные y2 и y3 считались скрытыми, их временные ряды неизвестными.

Угадать истинные (идеальные, наилучшие) стартовые догадки маловероятно, поэтому, следуя работам [30, 48], при оценке параметров в качестве стартовых догадок для всех переменных модели y1, yn+1,..., y(L1)n+1,..., y(1)Ln+(L1)n+1 использовались значения наблюдаемой в соответствующие моменты времени. Такой подход неплохо зарекомендовал себя, хотя Вообще говоря, зависимость r(µ) может быть неоднозначна, но при использовании r100 эта проблема снимается.

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

и неидеален6. Для минимизации целевой функции использовался обобщенный метод Гаусса-Ньютона [30].

На рис. 1.1a,b в разных масштабах представлены результаты исследования сходимости в глобальный минимум для алгоритма Бока (п. 1.3.1) при различных стартовых догадках для параметров. Картины получены в случае длины ряда, при которой алгоритм Бока даёт максимальный радиус стопроцентной сходимости. На рис. 1.1 приведены сечения пространства параметров плоскостью b1 = 0, поскольку в рассматриваемом случае ошибки в параметре b1 наименее критичны. Как видно из рис. 1.1a,b, область сходимости алгоритма Бока в случае незашумленного ряда достаточно широка диаметр окружности r100на плоскости нормированных параметров больше 0.7, т.е. допускаются ошибки в стартовых догадках больше 70%. Существуют и очень удалённые от истинных значений области, откуда также достигается сходимость в глобальный минимум (рис. 1.1а). Однако модифицированный метод ещё менее требователен к стартовым догадкам для параметров: сравните рис. 1.1b и рис. 1.1c, где значения r100, r90, r80 существенно больше, а соответствующая r100 белая область шире.

1.3.3 Оценка влияния измерительного шума на результат реконструкции На рис. 1.2 представлена зависимость меры µ(r) при разных уровнях измерительного шума, добавлявшегося к ряду аддитивно и имевшего нормальное распределение. Из рисунков видно, что даже весьма существенный по абсолютной величине шум до 20% от амплитуды сигнала тем не менее не приводит к качественному изменению зависимости µ(r). При некоторых значениях шума: noise/signal = 0.1, где noise и signal среднеквадратичные отклонения шума и сигнала, величина µ(r) может быть даже больше, чем в случае noise/signal = 0 для отдельных реализаций Вариант использования идеальных стартовых догадок для скрытых переменных рассмотрен отдельно далее на частном примере, показывающем, что принципиально результаты ничем не отличаются от представленных далее.

шума, как в приведённом на рис. 1.2 примере7.

Рис. 1.3.3, где построены графики r100 (N ) радиуса 100%-ной сходимости в зависимости от длины ряда, показывает что низкая чувствительность алгоритма Бока к шумам является его общим свойством, поскольку выполняется для большого количества разных длин ряда. Графики для случая зашумлённых рядов качественно ничем не отличается от графика, соответствующего реконструкции без шума. Кроме того, рис. 1.3.3, как и рис. 1.2b, показывает, что модифицированный метод имеет столь же хорошие свойства устойчивости к шумам, как и оригинальный, несмотря на более низкие требования к непрерывности траектории (на рисунках показаны результаты применения модифицированного метода при = 2 и = 4).

1.3.4 Оптимальный выбор параметров алгоритма: количества Общим правилом для большинства методов обработки экспериментальных данных, в том числе и методов реконструкции моделей по временным рядам, является стремление исследователя задействовать максимально большой объём имеющейся информации, то есть использовать максимально большой ряд. Однако, применяя метода Бока, мы сталкиваемся с необходимостью ограничить длину используемого ряда, чтобы увеличить вероятность попадания в глобальный минимум целевой функции (1.3).

На самом деле необходимо выбрать наилучшим образом не только длину ряда N, но и количество фрагментов L, на которые он делится. На и длину фрагмента. Чтобы охарактеризовать зависимость r100(L, n) была проведена большая серия экспериментов с различными L и n, в каждом эксперименте строилось пространство стартовых догадок, как на рис. 1.1, и подсчитывался радиус 100%-ной сходимости. Результаты представлены на Хотя радиус сходимости может немного увеличиваться с ростом шума, точность полученных оценок параметров всегда снижается, т.е. значение параметров, соответствующее глобальному минимуму удаляется от истинных значений c0.

Рис. 1.2: Зависимость меры µ от радиуса r при различных уровнях измерительного шума, (a) реконструкция методом Бока с параметрами L = 35, n = 30, (b) реконструкция при разных значениях числа сегментов непрерывности: = 1 метод Бока и=2 модифицированный метод.

Рис. 1.3: зависимость радиуса стопроцентной сходимости r100 от общей длины ряда N при разных уровнях шума. Реконструкция при разных значениях числа сегментов непрерывности: = 1 метод Бока и = 4 модифицированный метод.

рис. 1.4, где L и n отложены по осям, в величина r100 изображена оттенками серого: чем светлее тон, тем больше r100. Для удобства дополнительно проведены линии уровня r100 = 0.2, r100 = 0.4, r100 = 0.6, r100 = 0.8, r100 = 1.0.

Рис. 1.4 показывает, что искомая зависимость r100(L, n) устроена достаточно сложно. При этом, однако, легко выделить одно или несколько значений (L, n), при которых достигается наибольшее значение r100. При удалении от этих значений величина r100 плавно снижается, оставаясь достаточно большою в значительном интервале общих длин ряда N. Общность этих выводов для достаточно широкого класса режимов подтверждается сравнением графиков 1.4a и 1.4b, соответствующих различным значениям параметров.

1.3.5 Преимущества модифицированного метода Увеличить длину ряда без ущерба для статистических свойств получаемых оценок параметров c можно, как уже упоминалось в 1.2.2, применяя модифицированный метод множественной стрельбы. Чтобы показать, что при этом радиус 100%-ной сходимости не только не уменьшается, но даже в среднем немного увеличивается по сравнению с оригинальным методом бока, были построены аналогичные пункту 1.3.4 зависимости r100(L, n) при тех же значениях параметров по тем же рядам для случая = 2, т.е. с допуском одного разрыва в середине ряда (см. рис. 1.5).

О преимуществах модифицированного метода свидетельствует более светлый тон рис. 1.5 по сравнению с рис. 1.5. Также максимально достижимое значение r100 больше для модифицированного метода: не только область с r100 1 больше, но и внутри неё выделяется область с r100 1.

Кроме того, рассматривая белые гиперболы (линии N = const) на рис. 3, можно сделать вывод, что преимущества модифицированного метода сильнее проявляются для более длинных рядов. Поскольку для реконструкции использовался модифицированный метод множественной стрельбы с = 2, для приведённых зависимостей N = 2Ln.

Рис. 1.4: зависимость радиуса стопроцентной сходимости r100 от количества фрагментов L, на которые делится ряд, и от длины фрагмента n. Величина r100 показана оттенками серого, чем светлее тон, тем больше r100. Проведены линии уровня r100 = 0.2, r100 = 0.4, r100 = 0.6, r100 = 0.8, r100 = 1.0. Дополнительно на рисунок наложены белые гиперболы линии, вдоль которых общая длина ряда N постоянна, и надписана величина N для каждой из них. Рис. (a) соответствует параметрам c0 = 10, c0 = 8/3, c0 = 28, рис. (b) параметрами c0 = 10, c0 = 8/3, c0 = 46.

Рис. 1.5: зависимость радиуса стопроцентной сходимости r100 от количества фрагментов L, на которые делится ряд, и от длины фрагмента n. Величина r100 показана оттенками серого, чем светлее тон, тем больше r100. Проведены линии уровня r100 = 0.2, r100 = 0.4, r100 = 0.6, r100 = 0.8, r100 = 1.0. Дополнительно на рисунок наложены белые гиперболы линии, вдоль которых общая длина ряда N постоянна, и надписана величина N для каждой из них. Поскольку для реконструкции использовался модифицированный метод множественной стрельбы с = 2, для приведённых зависимостей N = 2Ln. Рис.

(a) соответствует параметрам c0 = 10, c0 = 8/3, c0 = 28, рис. (b) параметрами c0 = 10, c0 = 8/3, c0 = 46.

1.3.6 Выбор оптимального количества разрывов Из сопоставления рис. 1.4 и рис. 1.5 видно, что модифицированный метод при = 2 не всегда даёт результаты лучшие, чем оригинальный алгоритм Бока: для рядов коротких рядов (N < 1400 для рис. (a) и N < для рис. (b)) наблюдается обратная картина, а для очень длинных рядов (N > 4000 для рис. (a) и N > 1500 для рис. (b)) методы дают почти одинаково плохие результаты. Последний результат можно объяснить таким образом, что допуск одного разрыва для таких длинных рядов недостаточен. Поэтому были построены интегральные кривые r100 (N ) при различных : от 1 (метод Бока) до 4. При этом в каждом случае первоначально строилась такая же зависимость r100 (L, n), как на рис. 1.4 и рис. 1.5, а затем для каждой величины N (т.е. вдоль белой гиперболы) брались L и n, соответствующие наибольшему при этом N радиусу r100. Результаты представлены на рис. 1.6.

Холмоподобная форма зависимостей r100 (N ) определяется двумя факторами: при очень малых N не хватает точек для подгонки параметров, а при больших сказывается экспоненциальная чувствительность к начальным условиям за время порядка = 1/1 малые возмущения достигают макромасштабов. Кривые, соответствующие большим значениям, достигают бльших значений r100 и сдвинуты вправо, в сторону более длинных рядов. При этом шире становятся и интервалы по оси N с ненулевыми значениями r100, то есть с увеличением числа разрывов диапазон длин ряда, в котором можно получить хорошую модель, возрастает.

Чтобы показать общность этих выводов и типичность кривых на стартовых догадок для скрытых переменных, т.е. начальные условия y1, yn+1,..., y(L1)n+1 на первом же шаге алгоритма точно такие, какие соответствуют глобальному минимуму (1.3)8 (см. рис. 1.7). Заметим, что Тем не менее, эти величины подгоняются прежним способом и, поскольку стартовые догадки для параметров cs как правило выбираются далёкими от их истинных значений cs, на промежуточных этапах подгонки начальные условия будут отличаться от идеальных.

кривые на рис. 1.7 существенно выше, чем на рис. 1.6, так при = для реальных стартовых догадок максимум кривой r100 (N ) составляет 1.36, для идеальных 10. Также сами кривые являются существенно более плавными и максимум на них выражен более явно. Всё это является следствием более удачных стартовых догадок.

Несмотря на то, что рис. 1.7 и рис. 1.6 хорошо демонстрируют влияние длины ряда на оптимальный выбор, на основе их непосредственного анализа сложно сформулировать какие-либо практические рекомендации по выбору количества разрывов. Поскольку экспоненциальное разбегание хаотических траекторий с близкими начальными условиями является главною причиной неэффективности метода Бока, представляется разумным связать оптимальное значение с ляпуновским временем, где = 1/, старший ляпуновский показатель. Для этого графики рис. 1.7 и рис. 1. были перестроены таким образом, что теперь по горизонтальной оси отложена не общая длина ряда N, а длина сегмента непрерывности Ln = N/ (см. рис. 1.8). Для удобства введена вторая ось абсцисс сверху, по которой время отложено не в отсчётах, а в величинах ляпуновского времени.

Из рис. 1.8 видно, что для случаев = 1 4 оптимальные длины Ln соответствуют 1–3 ляпуновским временам, немного смещаясь в сторону меньших значений с ростом. Такой эффект можно объяснить со следующих позиций. Шансы на успех реконструкции зависят от допустимой длины сегмента Ln, на которой малые возмущения начальных условий ещё не слишком возрастают, а также от числа P +D свободных параметров, которые нужно подогнать по ряду. Увеличивая число разрывов, мы ослабляем чувствительность к начальным условиям, но одновременно увеличиваем количество подгоняемых величин. В приведённом примере при = 1 их всего 6, при = 2 9, а при = 4 уже 15. Поэтому с ростом общего количества подгоняемых величин результат реконструкции более критичен к ошибкам в начальных условиях и интервал значений Ln с ненулевым r сужается.

Рис. 1.6: зависимость радиуса стопроцентной сходимости r100 от общей длины используемого ряда N при различном числе разрывов ( 1).

Рис. 1.7: зависимость радиуса стопроцентной сходимости r100 от общей длины используемого ряда N при различном числе разрывов ( 1). Использованы идеальные стартовые догадки для скрытых переменных.

1.3.7 Универсальность оптимального выбора числа подсегментов L при различном числе сегментов непрерывности Решение на основе оценки ляпуновского времени вопроса об оптимальном выборе числа сегментов непрерывности для модифицированного метода множественной стрельбы не снимает вопрос об оптимальном выборе количества подсегментов L, на которые этот сегмент делится и для которых условие (1.4) выполняется. При построении графиков на рис. 1.6, 1. и 1.8 оптимальное L выбралось перебором, что весьма трудоёмко. В то же время, рис. 1.4 и рис. 1.5 показывают, что при одном и том же N успех моделирования решающим образом зависит от количества подсегментов L:

при изменении L всего в два раза величина r100 может измениться в 5 и более раз.

Составить для величины L какую-либо оценку из общих соображений о природе сигнала, подобно оценке для, затруднительно. Однако некоторые общие свойства алгоритма множественной стрельбы, относящиеся к величине L, можно выявить, если построить зависимость r100 (L) при различных и при оптимальном выборе длины сегмента непрерывности Ln (см. рис. 1.3.7).

Рис. 1.3.7 показывает, что оптимальное значение L для любого одно и то же (в данном случае равно 20). Аналогичные результаты были получены при других значениях параметров: c0 = 10, c0 = 8/3, c0 = 28, а также при использовании идеальных стартовых догадок для скрытых переменных. Это свойство можно использовать, например, в следующей ситуации.

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

Рис. 1.8: зависимости радиуса стопроцентной сходимости r100 от длины сегмента непрерывности Ln = N/ при различном числе разрывов ( 1): (a) при использовании реальных стартовых догадок для скрытых переменных, (b) при использовании идеальных стартовых догадок.

Рис. 1.9: зависимости радиуса стопроцентной сходимости r100 от длины подсегмента L при различном числе разрывов ( 1) и при оптимальном выборе длины сегмента непрерывности Ln.

1.3.8 Реконструкция системы Рёсслера Чтобы показать общность полученных в пп. 1.3.2–1.3.7 на примере системы Лоренца результатов аналогичные исследования были проведены для системы Рёсслера при значениях параметров c1 = 0.2, c2 = 0.15, c3 = 10, что соответствует хаотическому режиму, и начальных условиях y1 = 0.21, y2 = 6.5, y3 = 0.22.

Основной средний период колебаний9 составляет примерно 6.0, старший ляпуновский показатель = 0.1. Система 1.7 решалась численно методом Рунге-Кутты 4-го порядка с шагом 0.002 и интервалом выборки 0.01. В качестве единственной наблюдаемой использовалась переменная y1, остальные две переменные считались скрытыми.

Данная система была выбрана в качестве тестового примера, поскольку, являясь как и система Лоренца 1.6 эталонным объектом нелинейной динамики, она имеет принципиально иную форму аттрактора. В системе Лоренца колебания происходят вокруг двух симметрично расположенных в плоскости (y1, y2) неустойчивых точек [39]. Огрубляя динамику, можно сказать, что колебания происходят вблизи плоскости, перпендикулярной (y1, y2) и пересекающей её примерно по диагонали y1 = y2. Таким образом, одновременные значения y1, y2 весьма близки друг к другу. В то же время система Рёсслера имеет спиральный аттрактор с одним центром вращения, причём основная динамика происходит в плоскости (y1, y2), а переменные y1 и y2 сдвинуты друг относительно друга на примерно четверть характерного периода.

Проще всего восстановить систему 1.7 используя идеальные стартовые догадки для скрытых переменных. Рис. 1.10 аналогичен рис. 1.1 и показыХотя режим, соответствующий выбранным значениям параметров, хаотический, вследствие наличия явного пика в спектре [40] и возможности легко ввести фазу колебаний как угол поворота на плоскости y1, y2 основной период колебаний можно выделить.

вает область сходимости в глобальный минимум целевой функции (1.3) в пространстве нормированных стартовых догадок (1.5) для параметров (в сечении плоскостью b3 = 0) при использовании оригинального метода Бока (a) и модифицированного метода с допуском одного разрыва (b).

Видно, что при использовании модифицированного метода удаётся добиться существенного роста области сходимости (светлая область на рис. 1.10b существенно больше, чем на рис. 1.10a).

Чтобы более общим образом охарактеризовать влияние количества разрывов на сходимость метода, на рис. 1.11 построены зависимости радиуса стопроцентной сходимости от длины ряда r100(N ) при различных (от 1 до 3) и оптимальном подборе L. Сравнение рис. 1.11 и рис. 1.7 показывает, что качественно результаты для систем Рёсслера и Лоренца одинаковы и показывают преимущества использования модифицированного метода. Так, максимальный r100 составляет 0.9 при = 1 (метод Бока), 0.95 при = 2 и уже 1.05 при = 3. При этом преимущества большего числа допустимых разрывов проявляются для длинных рядов, как это было и в примерах из п. 1.3.6.

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

1.3.9 Подбор стартовых догадок для скрытых переменных При моделировании системы Лоренца в пп. 1.3.2–1.3.7 мы пользовались подходом к формированию стартовых догадок, предложенным в [30].

Согласно ему стартовые для скрытых переменных выбираются в точности равными значениям из ряда наблюдаемой: y2,j = j, y3,j = j 10, где Двойной индекс внизу обозначает первое число номер переменной, второе число номер момента времени.

Рис. 1.10: плоскость нормированных стартовых догадок для параметров системы Рёсслера (сечение пространства параметров плоскостью b3 = 0). Белым обозначены точки, откуда достигается глобальный минимум, оттенками серого локальные минимумы (интенсивность зависит от их удалённости от глобального), чёрным алгоритм минимизации расходится: (a), алгоритм Бока при L = 25, n = 30, (b) модифицированный метод с L = 25, n = 30, = 2. Использованы идеальные стартовые догадки для скрытых переменных.

Рис. 1.11: зависимость радиуса стопроцентной сходимости r100 от общей длины используемого ряда N при различном числе разрывов ( 1) и оптимальном подборе L реконструкция системы Рёсслера. Использованы идеальные стартовые догадки для скрытых переменных.

j = 1, (n + 1),... ((L 1)n + 1) начальные моменты каждого подсегмента. Как уже упоминалось в п. 1.3.8, для системы Лоренца переменные y1 и y2 очень близки по абсолютным величинам и практически синфазны.

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

Однако попытка переноса такого подхода на систему Рёсслера ведёт к полному провалу: радиус сходимости получается нулевым (см. рис. 1.12a).

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

Выход из такой ситуации заключается в том, чтобы при составлении стартовых догадок максимально учесть специфику модели. Для системы Рёсслера это достигается следующим образом. Переменные y1 и y2 в ней похожи, но сдвинуты примерно на четверть характерного периода колебаний. Следовательно, в качестве стартовых догадок для y2 нужно взять ряд, но не одновременный, а сдвинутый на примерно T /4 вперёд (в нашем случае T 6.0). Величину характерного периода можно определить по пику в спектре колебаний или введя фазу как угол поворота радиус вектора на плоскости (y1, y2) и подсчитав среднее время, за которое она нарастает на 2. В то же время для переменной y3 в качестве стартовых догадок следует взять тождественный ноль, поскольку основное движение на аттракторе происходит вблизи плоскости (y1, y2), при эпизодических и кратковременных выбросах по переменной y3. Хотя в некоторых редких точках погрешность в выборе y3,j будет велика, в большинстве точек и, следовательно, в среднем она будет достаточно мала.

Результат использования составленных таким образом стартовых догадок представлен на рис. 1.12c. Видно, что появилась довольно существенная область, где достигается сходимость в глобальный минимум целевой функции 1.3, хотя это область существенно меньше, чем при использовании идеальных догадок (см. рис 1.12b).

Теперь с использованием подобранных таким образом стартовых догадок легко повторить результаты п. 1.3.8. На рис. 1.13 показаны в сравнении результаты при использовании метода Бока при L = 25, n = 30 и модифицированного подхода при L = 25, n = 30, = 2. Модифицированный метод, как и ранее, даёт заметно лучшие результаты (светлая область на рис. 1.13b существенно шире, чем на рис. 1.13a).

На рис. 1.14(a) подтверждаются выводы, сделанные по рис. 1.6 и рис. 1.7 для системы Лоренца как при идеальных, так и при реальных стартовых догадках для скрытых переменных, а также по рис. 1.11 для системы Рёсслера при идеальных стартовых догадках. А именно, модифицированный метод даёт существенно больший максимальный радиус 100%-ной сходимости: для = 2 r100 = 0.097 против r100 = 0.068 для = 1.

Также преимущества модифицированного метода становятся существенными при больших длинах ряда: для = 1 максимуму r100 соответствует Рис. 1.14(b) соответствует рис. 1.8. Видно, что оптимальная длина сегмента непрерывности Ln примерно равна, как и для системы Лоренца, или несколько меньше его, что может быть обусловлено индивидуальными характеристиками аттрактора. При этом происходит небольшое смещение максимума зависимости r100(Ln) в сторону меньших Ln с ростом, что, как уже объяснялось в п. 1.3.6, обусловлено отрицательным влиянием увеличения количества независимых искомых параметров: их 6 (три собственно параметра c1, c2 и c3 и три начальных условия y1,1, y2,1 и y3,1) для = 1 и уже 9 (появляются дополнительно начальные условия y1,(Ln+1), y2,(Ln+1) и y3,(Ln+1)) для = 2.

1.4 Выводы С помощью введённого количественного критерия проведено сравнение работоспособности различных реализаций метода множественной Рис. 1.12: плоскость нормированных стартовых догадок для параметров системы Рёсслера (сечение пространства параметров плоскостью b3 = 0). Белым обозначены точки, откуда достигается глобальный минимум, оттенками серого локальные минимумы (интенсивность зависит от их удалённости от глобального), чёрным алгоритм минимизации расходится: (a), алгоритм Бока при L = 25, n = 30, (b) модифицированный метод с L = 25, n = 30, = 2.

Рис. 1.13: плоскость нормированных стартовых догадок для параметров системы Рёсслера (сечение пространства параметров плоскостью b3 = 0). Белым обозначены точки, откуда достигается глобальный минимум, оттенками серого локальные минимумы (интенсивность зависит от их удалённости от глобального), чёрным алгоритм минимизации расходится: (a), алгоритм Бока при L = 25, n = 30, (b) модифицированный метод с L = 25, n = 30, = 2. Использованы реалистичные специально составленные стартовые догадки для скрытых переменных.

Рис. 1.14: зависимость радиуса стопроцентной сходимости r100 от общей длины используемого ряда N (a) о от длины сегмента непрерывности Ln = N/ (b) при различном числе разрывов ( 1) и оптимальном подборе L реконструкция системы Рёсслера.

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

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

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

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

Модифицированный метод имеет ряд преимуществ по сравнению с алгоритмом Бока, поскольку он свободен от ограничений, связанных с высокой чувствительностью к начальным условиям: накапливающиеся вдоль фазовой траектории невязки отбрасываются в разрешенных этим подходом разрывах. Во-первых, он позволяет использовать ряды большой длины, что даёт возможность уточнить оценки параметров за счёт введения в рассмотрение дополнительного объёма данных, причём такие оценки статистически предпочтительнее, чем получаемые с помощью кусочного метода множественной стрельбы. Во-вторых, при использовании рядов фиксированной длины он в большинстве случаев предъявляет меньшие требования к стартовым догадкам для искомых параметров, а иногда позволяет получить достаточно точные оценки, когда с помощью алгоритма Бока это невозможно сделать ни при каких стартовых догадках, даже равных истинным значениям.

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

Показано, что влияние шума (вплоть до 0.2 от уровня сигнала) на процесс моделирования несущественно для обеих методик. Следует сказать, что точность оценки параметров при условии попадания в глобальный минимум целевой функции и для выбранной длины ряда наибольшая для алгоритма Бока, предъявляющего самые жёсткие требования при подгонке, меньше для модифицированного метода, и ещё меньше для кусочного метода множественной стрельбы. Если же учесть, что попасть в глобальный минимум гораздо проще при использовании модифицированного метода, то он в итоге имеет практические преимущества даже в смысле точности оценок. Но следует учитывать и то, что, так как модифицированный метод предъявляет меньшие требования к адекватности используемой модели, существует опасность при слишком малом размере сегментов успешно подогнать под наблюдаемый ряд чужую модель.

Основные результаты первой главы изложены в работах [84, 85, 86, 87, 88, 89] Глава Реконструкция систем под регулярным воздействием 2.1 Введение Системы под внешним воздействием широко распространённый в природе класс. Большинство изучаемых нами объектов каким-либо образом связаны с окружающею средою: одни либо могут принимать сигналы, подвергаясь её воздействию, либо сами влияют на неё, либо эта связь двунаправленная. Такая связь не всегда напрямую отражается в структуре модельных уравнений. Например, все представители весьма популярного в нелинейной динамике класса автогенераторных систем являются по форме автономными, хотя и получают энергию извне (через отрицательное трение ). Однако и неавтономные по форме модели также весьма распространены (см., например, в приложении к радиофизике [11, 57, 58]). В этой связи важным оказывается вопрос о выборе структуры модели при реконструкции по временным рядам, если есть сведения о неавтономности моделируемого объекта: предпочесть ли структуру, в явном виде содержащую внешнее воздействие, или автономную структуру.

Универсальные подходы, ориентируемые на возможно более широкий класс объектов, получили широкое распространение. В случае, когда измерению доступна лишь одна наблюдаемая величина, такой стандартной структурой обычно является (см. [27]) система дифференциальных уравнений (2.1):

где в качестве x1 берется сама скалярная наблюдаемая, а функция f представляется в виде степенного полинома порядка K:

Претензии таких моделей на общность теоретически обоснованы, поскольку в виде (2.1) может быть представлена почти любая система ОДУ (подробнее см. [28]) и любая непрерывная функция f в (2.1) может быть сколь угодно точно равномерно приближена алгебраическим полиномом (2.2) (теорема Вейерштрасса). Однако случаи успешного применения стандартного подхода на практике единичны, особенно если велики D и K, т.

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

Причина неудачи стандартного подхода применительно к неавтономным системам в значительной степени состоит в том, что если в нашем распоряжении имеется единственный скалярный временной ряд, необходимо восстанавливать пространство состояний каким-либо известным способом1 за изменение структуры приходится платить значительным ростом размерности: согласно теореме Таккенса [63] размерность до 2 раз.

Работоспособность стандартного подхода можно повысить за счет частичного отказа от универсальности и разработки методик (технологий), ориентированных на сравнительно узкий класс объектов. В данной работе такая методика предлагается для систем, находящихся под регулярным (имеющим дискретный спектр) внешним воздействием. Это могут быть произвольные по форме периодические или квазипериодические изменения параметров или внешней силы. Предпосылками для моделирования Можно использовать последовательное дифференцирование [59, 60], интегрирование [61] или метод временных задержек [10, 62].

объекта в виде системы ОДУ с регулярным внешним воздействием могут быть физические соображения, например, наличие дискретных пиков в спектре мощности наблюдаемого ряда2, или априорная информация. Методика опирается на введение явной зависимости от времени в структуру модельных уравнений, аналогично тому, как это делается для уменьшения D в [11, 36] при гармоническом воздействии, а также на использование, наряду со степенными полиномами (2.2), тригонометрических рядов.

2.2 Методика реконструкции Пусть имеется временной ряд наблюдаемой величины {i}, спектр мощности которой имеет дискретные пики на частотах 1, 2,... m и их комбинационные составляющие. Ограничимся для начала при моделировании предположением о силовом характере воздействия на объект и стандартною структурой (2.1), которую модифицируем путем замены:

где функция g(t) представляет внешнее воздействие. Если она известна, либо известна её форма и нужно определить только некоторые параметры, её можно непосредственно ввести в (2.3). Иначе её можно представить в виде суммы тригонометрических полиномов:

где Tj = 2/j. Для периодического случая m = 1 и формула (2.4) преобразуется к более простому виду:

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

А квазипериодическом m равно числу несоизмеримых частот.

Примем за переменную x1 саму скалярную наблюдаемую, а временdxD (tj ) ные ряды величин x2 (ti ),..., xD (ti ), получим путем последовательdt ного дифференцирования ряда (tj ) после фильтрации шума. Для нахождения коэффициентов полинома f и функции g используем метод наименьших квадратов (МНК), минимизируя средний квадрат ошибки аппроксимации3 :

Использование структуры (2.1, 2.3, 2.5) вносит специфику в процедуру расчета коэффициентов модели по сравнению со стандартным подходом в рассматриваемом случае параметры Tj входят в выражение для воздействия g нелинейно. Поэтому приходится искать минимум погрешности аппроксимации 2 как функции многих переменных с помощью итерационных методов. Мы использовали метод Левенберга–Марквардта, но при этом значения линейно входящих параметров на каждом шаге рассчитывали линейным МНК. Согласно этому методу выбираются начальные приближения для всех неизвестных коэффициентов; функции f и g линеаризуются по этим коэффициентам в окрестности этих приближений; решением задачи на линейный МНК находятся поправки к сделанным приближениям. Алгоритм повторяется, пока не будет достигнут минимум (2.6), в общем случае локальный. Чтобы найти глобальный минимум, необходимо перебирать различные начальные приближения. В качестве начального приближения значений Tj удобно использовать экспериментальную оценку этих величин по спектру.

В общем случае для оценки параметров по данным используют метод максимального правдоподобия, частным случаем которого является МНК. В работе [55] показано, что в общем случае МНКоценки могут быть сильно смещены, но при рассматриваемых нами малых уровнях данным эффектом можно пренебречь. Поэтому далее используется только МНК, при котором сильно упрощаются расчеты.

2.3 Численный эксперимент Для апробации метода в качестве тестового объекта был выбран нелинейный неавтономный осциллятор Тоды:

где r = 0.5 параметр диссипации, а (t) внешнее воздействие. Уравнения (2.7) решались численно методом Рунге-Кутты 4-ого порядка с шагом t = 0.01 и таким же интервалом выборки. При выбранном значении управляющего параметра r и рассмотренных ниже типах воздействия в системе реализовывался хаотический режим. В каждой ситуации восстанавливались модели как в стандартной форме (2.1, 2.2), так и в модифицированном виде (2.1, 2.3, 2.5).

Для изучения работоспособности изложенного алгоритма реконструкция проводилась при четырёх различных формах регулярного воздействия (рис. 2.1). При этом воздействие последовательно усложнялось: сначала было рассмотрено периодическое гладкое, затем периодическое треугольное, более сложное для разложения в ряд, затем при наличии субгармоник и, наконец, квазипериодическое.

2.3.1 Реконструкция при гладком периодическом воздействии Первое воздействие представляло собою гауссовы импульсы переменной полярности, описываемые формулою (2.8) где A амплитуда воздействия, T его период, эффективная ширина импульса, а вспомогательная функция (t) соотношением (2.9):

где t0 начальный сдвиг по времени (начальная фаза воздействия). Такое воздействие было выбрано, поскольку в отличие от гармонического, рассмотренного в [36], оно имеет богатый набор кратных гармоник основной Рис. 2.1: рассмотренные виды внешнего воздействия: периодические (a) гауссовы импульсы, (b) треугольные импульсы; (c) сигнал с субгармониками (получен решением системы Рёсслера (1.7) при значениях параметров c1 = 0.1, c2 = 0.2, c3 = 5.7, что соответствует режиму периода 4); (d) квазигармонический двухчастотный сигнал.

частоты, но, в то же время, также является гладким. Величины A, T и варьировались в различных экспериментах в широких пределах: A [10, 25], T [2, 5], [0.15, 0.2].

Величины D, K и k подбирались оптимальным образом. Критериями качества моделей выступали погрешности аппроксимации, дальность прогноза pred и сходство фазовых портретов модели и объекта. Дальность прогноза pred рассчитывалась как временной интервал, на котором ошибка прогноза не превышала 5% от стандартного отклонения; этот интервал сравнивался с характерным временным масштабом воздействия его периодом T (или одним из периодов в квазипериодическом случае) или предельно возможной дальностью прогноза limit, определяемой величиною старшего ляпуновского показателя. Для анализа влияния шума к исходному временному ряду добавлялся нормальный шум; тогда для расчета производных использовался m-точечный сглаживающий полином (фильтр Савицки–Голэя [56]) второго порядка, величина m подбиралась.

При моделировании весь имеющийся в распоряжении временной ряд разбивался на две неравные части: по одной из них строилась модель, этот ряд был назван тренировочным, его длина составляла Ntrn, другая часть называлась тестовым рядом с длиною Ntest, на нем проверялась адекватность полученной модели, подсчитывались ошибка аппроксимации 4 и дальность прогноза pred. Заметим, что хорошее соответствие модели объекту на тренировочном ряде вовсе не означает, что эта модель хорошая, поскольку она может описывать не процесс как таковой, а данную конкретную реализацию, с присущими только её особенностями и шумами (подробнее см. [36, 59]).

Стандартный алгоритм (2.1, 2.2) не дает удовлетворительных результатов даже в отсутствие шумов. Наилучшая из полученных моделей (при использовании D = 3 и K = 2) демонстрирует хаотический режим в той же области, что и объект, но фазовый портрет её не похож на портрет объекта (ср. рис. 2.2a и рис. 2.2b), ошибка аппроксимации при этом очень При вычислении тестовой ошибки аппроксимации использовался ряд такой же длины, как и тренировочный.

велика = 86.9%, а дальность прогноза очень низкая pred = 0.078T. При этом предельно возможная дальность прогноза с точностью 5%, определяемая старшим ляпуновским показателем составляет для рассматриваемой системы в данном режиме limit 60T, т.е. в 800 раз больше. Следует отметить, что данная модель получилась при минимально возможных значениях размерности модели D = 3 и степени аппроксимирующего полинома K = 2: взять меньшую размерность нельзя, поскольку исходный сигнал хаотический и, следовательно, не может описываться двумерной моделью, а при K = 1 получилась бы система линейных уравнений, которая не имеет автоколебательных решений. В то же время согласно теореме Такенса объекта, гарантирует точное и однозначное воспроизведение моделью реализации объекта5. В нашем случае D0 = 3, следовательно Dmax = 7, но попытки построить модели с размерностью D = 4 и выше неизбежно заканчивались неудачею: такие модели не только демонстрировали огромную ошибку аппроксимации > 60% по тренировочному ряду и очень низкое время предсказания pred < 0.03T, но их поведение даже качественно не согласовывалось с экспериментально наблюдаемым: в качестве решения получался предельный цикл, устойчивое состояние равновесия (как правило с переходным процессом, который и описывал тренировочный ряд) или вовсе расходилось на бесконечность. К аналогичным последствиям приводило увеличение K > 2. Причина этого кроется как в росте влияния шумов при формировании рядов переменных, что препятствует увеличению размерности, так и в быстром росте числа лишних коэффициентов с увеличением степени полинома, так как эти коэффициенты учитывают не динамику объекта в данном режиме вообще, а особенности конкретной реализации.

В то же время модернизированная структура (2.1, 2.3, 2.5) при D = и K = 8 и достаточном количестве учтённых гармоник k (уже начиная с k = 8) обеспечивает хорошее качество реконструкции (см. рис. 2.3). КачеУсловие Dmax = 2D0 + 1 является достаточным, но не необходимым, на практике часто более адекватными оказываются модели относительно невысокой размерности ственно фазовые портреты объекта и модели идентичны, ошибка аппроксимации составляет менее 0.1% и время предсказания pred = 19.24T (при k = 15). При наличии измерительных шумов амплитуды до 1% результаты сохраняются (ср. рис. 2.3b и 2.3c). Успех модифицированного метода обусловлен как уменьшением количества численных производных (с трёх и более для стандартной модели до двух), так и более удачным выбором аппроксимирующей функции (2.3) по сравнению с (2.2), что привело к снижению количества лишних коэффициентов.

Существенным фактором успешности модифицированного метода является правильный выбор количества членов тригонометрического полинома. На рис. 2.4 представлены графики зависимостей ошибки аппроксимации и времени предсказания pred от степени тригонометрического полинома k. Из графика видно, что существует оптимальное с точки зрения pred значение k. При меньших k количество членов ряда Фурье ещё недостаточно для описания формы воздействия (2.8), а при больших k новые гармоники оказываются уже лишними: точности данных не хватает, чтобы оценить их коэффициенты и они отражают только индивидуальные характеристики тренировочного ряда. Ступенчатый характер зависимостей (k) и pred (k) определяется видом исходного внешнего воздействия, которое не содержит чётных гармоник, следовательно все они являются лишними коэффициентами и их учёт несущественно улучшает модель.

2.3.2 Влияние шума на результат реконструкции Цель данного раздела показать влияние измерительного шума на результат реконструкции, что уже было начато в п. 2.3.1. Для этого мы добавляли к временным рядам наблюдаемой нормальный шум с нулевым средним и различным среднеквадратичным отклонением n. Модель восстанавливалась только в модифицированном виде, поскольку, как показано в п. 2.3.1, стандартная структура (2.1, 2.2) была неэффективна даже при нулевом шуме. Поскольку в наибольшей степени влияние шума сказывается на качестве модели вследствие необходимости получать все переменРис. 2.2: Фазовые портреты объекта (a) и модели, построенной с использованием стандартной структуры (2.1, 2.2.) с D = 3 и K = 2.

Рис. 2.3: Фазовые портреты объекта (a) и модели с использованием модифицированной структуры (2.1, 2.3, 2.5) k = 8 без добавленного шума (b) и с 1%-ным шумом (с) при k = 15.

ные кроме первой численным дифференцированием ряда наблюдаемой, мы использовали при дифференцировании сглаживающий полином второго порядка (фильтр Савицки–Голэя [56]). При этом количество точек, по которым проводилось усреднение, подбиралось для каждого уровня шума отдельно.

Чтобы охарактеризовать величину подаваемого шума мы ввели велиn чину = lg, где n и s среднеквадратичные отклонения шума и сигнала соответственно. Зависимости pred () и limit() приведены на рис. 2.5.

Из приведённых графиков видно, что, хотя с ростом шума дальность прогноза снижается (при = 4 pred = 7.1, а при = 2 уже pred = 1.5), предельно возможное время прогноза также снижается (limit = 18.4 и limit = 4.8) и относительное значение pred /limit остаётся примерно таким же: pred /limit = 0.38 при = 4 и pred /limit = 0.32 при = 2. Таким образом, можно говорить о том, что при разумных уровнях шума до 1% по амплитуде (отношению среднеквадратичных отклонений) модифицированный метод вполне работоспособен. Главное разумно выбрать количество точек при использовании сглаживающего полинома: малое усреднение приведёт к росту влияния шумов, а слишком сильное исказит сигнал.

2.3.3 Реконструкция при треугольном периодическом воздействии Второе рассмотренное нами воздействие аналогично первому, но использованы не гауссовы импульсы, а треугольные (см. рис. 2.1b). Такая угловатая функция существенно сложнее для аппроксимации тригонометрическим полиномом, поскольку приходится учитывать существенно большее, чем для (2.9) количество членов. Цель данного численного эксперимента подтвердить общность полученных в п. 2.3.1 выводов и выяснить, насколько предлагаемый модифицированный метод устойчив к росту степени тригонометрического полинома. В представленном примере амплиРис. 2.4: зависимости (в %) и pred (в периодах внешнего воздействия T ) от числа гармоник k в аппроксимации внешней силы, при воздействии гауссовыми импульсами с T = 3.6, A = 20, t0 = 0, = 0. Рис. 2.5: зависимости дальности прогноза pred построенной с использованием модифицированного алгоритма модели и предельно возможной дальности прогноза limit (в периодах внешнего воздействия T ) от уровня аддитивного шума = lg, где n и среднеквадратичные отклонения шума и сигнала соответственно.

туда воздействия A = 20 (это амплитуда униполярного импульса, соответственно общая амплитуда в 2 раза больше), период T = 3.6, общая ширина импульсов = 0.7.

Несмотря на то, что при использовании стандартного подхода как размерность, так степень полинома варьировались при построении модели в широких пределах: перебирались всевозможные D от 3 до 6 и K от 2 до 8, построить модель, которая демонстрировала бы хаотические поведение, не удалось. Все построенные модели имели в качестве аттрактора либо предельный цикл, либо устойчивую точку, либо вовсе не имели аттрактора, а тренировочный временной ряд описывался за счёт переходного процесса.

Модифицированный метод дал удовлетворительный результат: при оптимальных значениях D = 2, K = 8 и k = 40 ошибка аппроксимации составила = 0.2%, а дальность прогноза pred = 13.7T при максимально возможном limit 70T. Несколько худшие, чем в предыдущем примере результаты объясняются сложностью аппроксимации воздействия тригонометрическим полиномом.

На рис. 2.6 приведены аналогичные рис. 2.4 зависимости ошибки аппроксимации и дальности прогноза pred от количества использованных членов разложения k. В целом, она имеет аналогичный вид, но насыщение происходит при существенно больших значениях k: при k = 35 для величины pred против k = 15 в предыдущем случае и при k = 85 для (хотя (k = 35) и (k = 85) отличаются всего в 1.3 раза) при том же k = 15 для гауссовых импульсов. Разный порядок насыщения для pred и объясняется, по-видимому, тем, что первая величина вычисляется по существенно более длинному тренировочному ряду, т.е. аттрактор охватывается более полно и, как следствие, коэффициенты при высоких гармониках, вычисленные не достаточно достоверно из-за точности представления данных, приводят к большим ошибкам в более удалённых от тренировочного ряда областях аттрактора.

Модифицированный подход показал большую устойчивость к росту количества коэффициентов тригонометрического полинома. Как показано на рис. 2.6, при k = 100 полученная модель почти ничем не хуже модели, соответствующей оптимальному выбору k = 35. Мы использовали и более длинные полиномы до 500 гармоник (т.е. 1000 коэффициентов). При этом нет никаких проблем с неустойчивостью полученного решения, свойственных моделям в стандартной форме при высоких степенях полинома (2.2).

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

2.3.4 Реконструкция при воздействии c субгармониками В представленных в пп. 2.3.1 и 2.3.3 примерах спектр внешнего воздействия содержал только гармоники основной частоты, что обусловило простоту аппроксимации функции (t) тригонометрическим полиномом (2.5).

Но на практике может встретиться такое воздействие, которое имеет не только высшие гармоники, но и субгармоники. Мы сконструировали такой пример, в котором в качестве сигнала воздействия выступала переменная y1 системы Рёсслера (1.7) при значениях параметров c1 = 0.15, c2 = 0.2 и c3 = 5.7, что соответствует режиму периода 4 (см. рис. 2.1c).

Стандартный подход, как и в случае воздействия треугольными импульсами (см. п. 2.3.3) не позволяет построить удовлетворительную модель:

ни одна из полученных моделей не демонстрировала хаотическое поведение. Модифицированный подход позволяет добиться хороших результатов, но для этом в качестве основного периода в (2.5) получается период, соответствующий наименьшей из субгармоник. Величина T1 в (2.5) подгоняется по ряду, как и другие коэффициенты, но при этом она входит в аппроксимирующую функцию (2.3) нелинейно, следовательно, если неудачно задать стартовую догадку для T1, можно попасть в локальный минимум целевой функции (2.6). Именно это и происходит, если в качестве стартовой доs гадки T1 использовать значение частоты, соответствующее пику в спектре мощности наблюдаемого сигнала. Следовательно, если есть основания полагать, что сигнал воздействия может содержать субгармоники, необходимо пробовать различные стартовые догадки для периода воздействия, в первую очередь значения T1, кратные значению периода, соответствующего основному пику в спектре мощности наблюдаемой.