авторефераты диссертаций БЕСПЛАТНАЯ БИБЛИОТЕКА РОССИИ

КОНФЕРЕНЦИИ, КНИГИ, ПОСОБИЯ, НАУЧНЫЕ ИЗДАНИЯ

<< ГЛАВНАЯ
АГРОИНЖЕНЕРИЯ
АСТРОНОМИЯ
БЕЗОПАСНОСТЬ
БИОЛОГИЯ
ЗЕМЛЯ
ИНФОРМАТИКА
ИСКУССТВОВЕДЕНИЕ
ИСТОРИЯ
КУЛЬТУРОЛОГИЯ
МАШИНОСТРОЕНИЕ
МЕДИЦИНА
МЕТАЛЛУРГИЯ
МЕХАНИКА
ПЕДАГОГИКА
ПОЛИТИКА
ПРИБОРОСТРОЕНИЕ
ПРОДОВОЛЬСТВИЕ
ПСИХОЛОГИЯ
РАДИОТЕХНИКА
СЕЛЬСКОЕ ХОЗЯЙСТВО
СОЦИОЛОГИЯ
СТРОИТЕЛЬСТВО
ТЕХНИЧЕСКИЕ НАУКИ
ТРАНСПОРТ
ФАРМАЦЕВТИКА
ФИЗИКА
ФИЗИОЛОГИЯ
ФИЛОЛОГИЯ
ФИЛОСОФИЯ
ХИМИЯ
ЭКОНОМИКА
ЭЛЕКТРОТЕХНИКА
ЭНЕРГЕТИКА
ЮРИСПРУДЕНЦИЯ
ЯЗЫКОЗНАНИЕ
РАЗНОЕ
КОНТАКТЫ


Pages:     | 1 |   ...   | 6 | 7 || 9 | 10 |

«Б.П. Безручко, Д.А. Смирнов МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ И ХАОТИЧЕСКИЕ ВРЕМЕННЫЕ РЯДЫ Издательство ГосУНЦ «Колледж» Саратов, 2005 УДК ...»

-- [ Страница 8 ] --

Это оценка смещенная: она занижена, т.к. знаменатель больше числителя. На рис.8.2 показана зависимость этого асимптотического значения от отношения шум/сигнал для нормально распределенного шума. Оно практически не отличается от истинного значения только при малых шумах (смещение меньше 1%, если уровень шума меньше 0.05). Смещение оценки тем больше, чем больше уровень шума. Рис.8.2. Оценки параметра c в Аналогичные свойства наблюдаются и при (8.1), полученные с помощью других распределениях шума [272, 50]. Но обычного МНК, при N =, c0 = 2 и дисперсии шума.

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

8.1.2.2. Как повысить точность оценок в случае значительного измерительного шума? Отчасти это удается при использовании так Часть II. Моделирование по временным рядам называемого метода полных наименьших квадратов [246], когда минимизируется сумма квадратов ортогональных расстояний от точек ( n, n+1 ) до графика f ( xn, c), рис.8.1,в. При этом учитывается то обстоятельство, что отклонение наблюдаемых точек с координатами ( n, n+1 ) от графика искомой функции f ( xn, c0 ) вызвано действием шума на обе координаты. Поэтому отклонение может происходить в любом направлении, а не только по вертикали. Использование именно ортогональных расстояний обосновано в [246] как приближенный вариант ММП.

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

( n+1 ) N 1 f ( n ) ( x1, c) min, S (c, x1 ) = (8.4) n = где f (n ) – n-я итерация отображения x n+1 = f ( x n, c), f ( 0) ( x, c) = x, и в число оцениваемых величин включено начальное состояние модели x1.

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

Поэтому дисперсия оценок (8.4) убывает в таком случае очень быстро с ростом N, иногда даже Рис.8.3. Целевые функции для квадратичного экспоненциально [287, отображения (8.1) при N = 20, c0 = 1.85, x1 = 0.3 : 239]. Это желательное слева – целевая функция для итераций в прямом свойство, но оно имеет времени (8.4), справа – целевая функция для место на практике только итераций в обратном времени (8.5) при условии, что всегда удается находить глобальный минимум (8.4). Даже при умеренном N график функции S в случае хаотической системы становится сильно «изрезанным» (рис.8.3,а), так что найти глобальный минимум численными методами [73] практически невозможно. Для этого потребовалось бы очень «удачное» задание стартовых догадок для c и x1. Об асимптотических свойствах оценок тоже говорить трудно, т.к. в пределе N целевая Глава 8. Модельные уравнения: оценка параметров функция становится негладкой. Поэтому разрабатываются модификации ММП в приложении к задаче оценки параметров модели по хаотическому ряду [287, 316].

Так, предлагалось делить исходный ряд на сегменты небольшой длины L, по каждому из которых можно найти минимум (8.4), и усреднять полученные оценки (кусочный метод). Это целесообразный подход, но итоговая оценка может оставаться асимптотически смещенной, а ее дисперсия убывает опять только как N 1. Некоторые варианты улучшения свойств оценок описаны ниже (п. 8.2).

Здесь мы отметим только специальный подход, предложенный в [158] для случая одномерных отображений. Он основан на использовании того свойства, что единственный ляпуновский показатель одномерного хаотического отображения в обратном времени становится отрицательным и траектория отображения не столь чувствительна к параметрам и «начальному» условию. Поэтому минимизируется величина ( N n ) N 1 f ( n ) ( x N, c) min, S (c, x1 ) = (8.5) n = где f ( n ) – n-ая итерация отображения x n+1 = f ( x n, c) в обратном времени.

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

При малых и умеренных уровнях шума ( x до 0.05-0.15) погрешности метода (8.5) оказываются меньше, чем для кусочного. Причем для малых шумов (8.5) дает асимптотически несмещенные оценки и дисперсия убывает в типичном случае как N 2. Эта высокая скорость обусловлена возвратами траектории в окрестность экстремума функции f [316].

8.2. Скрытые переменные Когда уровень измерительного шума значителен, переменную состояния x часто считают «скрытой», т.к. ее значения, строго говоря, не известны. «Еще более скрытыми» являются переменные, зашумленные значения которых нельзя ни непосредственно измерить, ни вычислить из рядов наблюдаемых величин, что часто имеет место на практике. В этом случае получение оценок параметров гораздо более проблематично, чем в задачах п. 8.1. Но если удается успешно провести эту процедуру, то, как ее побочный продукт, появляется дополнительная возможность – получить временные ряды скрытых переменных. Тогда процедура моделирования опять выступает как измерительный прибор, но уже и в отношении динамических переменных, а не только параметров.

Часть II. Моделирование по временным рядам 8.2.1. Методы Мы проиллюстрируем их на примере оценки параметров обыкновенных дифференциальных уравнений без динамического шума.

Объект – классическая хаотическая система (система Лоренца):

x1 = с1 ( x 2 x1 ), & x 2 = x 2 + x1 (с3 x3 ), (8.6) & x3 = с 2 x3 + x1 x 2, & с параметрами c1 = 10, c2 = 8 3, c3 = 46. Наблюдается только зашумленная реализация величины x1 : n = x1 (t n ) + n, переменные x2 и x3 – скрытые.

Модель имеет вид (8.6), все три параметра сi считаются неизвестными.

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

Непосредственное решение задачи типа (8.4) называют методом подгонки начального условия [238, 334]. Согласно сказанному в п.8.1.1.2, он неприменим для длинных хаотических рядов. Простое деление ряда на сегменты с последующим усреднением оценок дает низкую точность итоговой оценки, а итерации в обратном времени для многомерной диссипативной системы не пригодны.

8.2.1.2. Метод множественной стрельбы. Отчасти обойти трудности и использовать более длинные ряды позволяет алгоритм Бока [196, 187], который называют еще методом множественной стрельбы, т.к. при этом от решения задачи Коши для получения траектории модели на всем интервале наблюдения переходят к решению набора краевых задач. А именно, разбивают исходный ряд {1, 2,..., N } на L неперекрывающихся сегментов длины M и рассматривают начальные состояния модели x (i ) на каждом их них (т.е. в моменты времени t (i 1) M +1, i = 1,..., L ) как оцениваемые величины, но не свободные. Решается задача условной минимизации, которая в случае скалярной наблюдаемой = h(x) + (в нашем примере = x1 + ) записывается в виде [ )] ( L M 1 S (c, x (1), x ( 2),..., x ( L ) ) = (t (i 1) M + n ) h x(t (i 1) M + n, x (i ), c) min, (8.7) i =1 n = x(t M, x (i ), c) = x (i +1), i = 1,..., L 1.

Здесь выражение x(t, x (i ), c) задает траекторию модели (решение модельных уравнений), т.е. состояние модели x в момент времени t при Глава 8. Модельные уравнения: оценка параметров начальном условии x (i ) и параметрах c. Первая строка означает минимизацию отклонения реализации модели от ряда для всего интервала наблюдения, вторая – «сшивание» сегментов, чтобы в итоге получить непрерывную траекторию модели на всем интервале наблюдения.

«Сшивание» накладывает L 1 условие типа равенства на L искомых векторов x (i ), т.е. только один из них можно считать «свободным параметром задачи».

Как обычно, задача решается численными итерационными методами, начиная со стартовых догадок для искомых величин c, x (1), x ( 2),..., x ( L ).

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

Рис.8.4. Оценка параметров по хаотической реализации координаты x = x1 системы Лоренца (8.6), N = 100 точек, выборочный интервал 0.04, нормальный измерительный шум со стандартным отклонением 0.2 от стандартного отклонения сигнала [238]:

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

б) метод множественной стрельбы. Процесс подгонки сходится к глобальному минимуму, где траектория модели и оценки параметров близки к истинным Часть II. Моделирование по временным рядам На рис.8.4 представлен пример применения методов для оценки параметров системы (8.6). Метод множественной стрельбы «находит»

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

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

8.2.1.3. Модификация метода множественной стрельбы. Как показано в [41], снизить остроту проблемы можно, допустив разрывность итоговой траектории модели в некоторые моменты времени на интервале наблюдения, т.е. отказавшись от некоторых равенств во второй строке (8.7). При этом легче отыскать глобальный минимум целевой функции S в (8.7).

Такая модификация позволяет использовать сколь угодно длинные хаотические ряды, а расплата состоит в том, что иногда модель с неадекватной структурой может быть признана «хорошей» за счет способности воспроизвести короткий участок ряда. Поэтому требуется осторожный выбор числа и размера сегментов непрерывности траектории модели. Отметим также дополнительную трудность ситуации со скрытыми переменными: кроме удачных стартовых догадок для параметров c очень важно найти и удачные стартовые догадки для скрытых переменных – компонент x (i ) – в отличие от чрезмерно оптимистичных ранних утверждений [187]. Часто приходится действовать наугад, но полезную информацию можно получить из предварительного изучения свойств модельных уравнений при пробных значениях параметров [41].

8.2.1.4. Заключительные комментарии. Существуют и развиваются способы оценки параметров и скрытых переменных, пригодные в случае одновременного наличия динамического и измерительного шумов. Они основаны на байесовском подходе [274, 217, 198] и использовании модифицированного фильтра Калмана [307, 334]. На эту широкую область исследований мы обращаем внимание читателя, но рассмотреть не имеем возможности.

Верификация моделей в рассмотренных в данной главе задачах о прозрачном ящике проводится с помощью двух основных подходов. Во первых, это анализ остатков (невязок) модели, т.е. проверка их соответствия предполагаемым свойствам шумов (см. п. 7.3) [45]. Во вторых, это расчет динамических, геометрических и топологических Глава 8. Модельные уравнения: оценка параметров характеристик аттрактора модели и их сравнение с соответствующими свойствами объекта [233], см. п. 10.4.

8.2.2. Что дают успехи и какая польза от неудач моделирования?

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

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

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

8.2.2.1. Пример из клеточной биологии. Для многих приложений необходимо выяснить, изменение каких свойств клеток ведет к наиболее эффективному подавлению нежелательного процесса в них, и как целенаправленно воздействовать именно на эти свойства.2 Для ответа на Так, рост раковых клеток обусловлен тем, что они производят вещества «неадекватно» окружающей обстановке. Вариант борьбы с заболеванием, который пока только в предположениях, мог бы опираться на эмпирическое моделирование, подобное описанному в п. 8.2.2.1. Здесь свойства биохимического процесса в клетках типа BaF3 [321] и возможность управления ими выяснены с опорой на экспериментальный временной ряд и привлечение биохимических соображений для формирования структуры модели.

Часть II. Моделирование по временным рядам поставленные вопросы может быть достаточно получения адекватной математической модели. Рассмотрим пример успешного моделирования сложного биологического процесса по временным рядам.

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

В частности, они обеспечивают размножение, дифференциацию и выживание клеток. Здесь рассмотрен так называемый сигнальный путь JAK-STAT [321], который преобразует внешний химический сигнал в активацию транскрипции с соответствующего гена внутри ядра клетки, см. схему на рис.8.5. Одна из Рис.8.5. Схема биохимического самых простых математических сигнального процесса в клетке [321] моделей этого биохимического процесса может быть записана, исходя из закона действующих масс (обычный подход в химической кинетике), и имеет вид:

dx1 dt = k1 x1 E (t ), dx2 dt = k 2 x2 + k1 x1 E (t ), (8.8) dx3 dt = k3 x3 + k 2 x2 2, dx4 dt = k3 x3.

Здесь ki – параметры скоростей реакций, E (t ) (Epo на рис.8.5) – концентрация вещества эритропойетина в окружающей среде, в ответ на изменение которой активируются соответствующие рецепторы в мембране клетки. Рецепторы связываются с тирозинкиназами типа JAK-2, которые имеются в цитоплазме клетки. Последние вступают в реакцию с молекулами вещества STAT5, концентрация которых обозначена x1. В результате реакции последние фосфорилируются. Получаются мономерные фосфорилированные молекулы STAT5, концентрация которых обозначена x2. Эта реакция ведет к убыванию концентрации x (первое уравнение (8.8)) и увеличению концентрации x2 (положительное слагаемое во втором уравнении (8.8)). Мономерные молекулы, встречаясь друг с другом, димеризуются. Концентрация димеризованных молекул – x3. Эта реакция ведет к убыванию x2 (отрицательное слагаемое во втором Глава 8. Модельные уравнения: оценка параметров уравнении (8.8)) и увеличению x3 (положительное слагаемое в третьем уравнении (8.8)). Димеризованные молекулы проникают в ядро, их концентрация в ядре – x4. Этот процесс ведет к убыванию x (отрицательное слагаемое в третьем уравнении (8.8)) и увеличению x (четвертое уравнение (8.8)). В ядре димеризованные молекулы активируют транскрипцию гена, обозначаемого CIS (target gene на рис.8.5). В результате производится некоторый белок. После этого димеризованные молекулы распадаются на мономерные, которые согласно гипотезе, лежащей в основе модели (8.8), деградируют в ядре.

Но есть и другая версия: мономерные молекулы STAT возвращаются с некоторой задержкой из ядра в цитоплазму клетки. При таком предположении математическая модель несколько меняется и принимает вид:

dx1 (t ) dt = k1 x1 (t ) E (t ) + 2k 4 x3 (t ), dx 2 (t ) dt = k 2 x 2 (t ) + k1 x1 (t ) E (t ), (8.9) dx3 (t ) dt = k 3 x3 (t ) + k 2 x 2 (t ) 2, dx 4 (t ) dt = k 4 x3 (t ) + k 3 x3 (t ), где появляются слагаемые с запаздыванием в первом и четвертом уравнениях, – время запаздывания. Какая именно из двух моделей (т.е.

какая из двух гипотез) более адекватна – неизвестно. Возможности наблюдения крайне ограничены, техника весьма дорогостояща и процесс измерений очень сложен [321]. В эксперименте можно было измерять с точностью до постоянного множителя только суммарную массу STAT5 в цитоплазме 1 и суммарную массу фосфорилированного STAT5 в цитоплазме 2 :

1 = k 5 ( x2 + 2 x3 ), (8.10) 2 = k 6 ( x1 + x2 + 2 x3 ), где k 5, k 6 – неизвестные коэффициенты пропорциональности. Кроме этих двух наблюдаемых измерялась и величина концентрации E(t), и тоже с точностью до коэффициента пропорциональности (рис.8.6,а). Подчеркнем, что все переменные моделей – скрытые, есть лишь две наблюдаемых, которые связаны с четырьмя переменными известным образом (8.10) Часть II. Моделирование по временным рядам Рис.8.6. Результаты моделирования биохимического сигнального процесса в клетке [321]: а) временной ряд внешнего воздействия – изменения концентрации эритропойетина;

б) результаты моделирования с помощью (8.8) по двухканальному ряду;

в) результаты моделирования с помощью (8.9) В работе [321] по описанным экспериментальным данным оценивались параметры обеих моделей (8.8) и (8.9). Была показана неадекватность первой и хорошее соответствие с экспериментом – второй.

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

Модель (8.8) оказалась не способной воспроизвести сигнал (рис.8.6,б), а модель (8.9) вполне адекватна в этом отношении (рис.8.6,в), оценка времени запаздывания составила 6 минут, что согласуется по порядку величины с результатами других авторов, полученными другими методами для похожих объектов.

Таким образом, только с помощью моделирования по рядам авторы смогли сделать нетривиальный вывод о том, что возврат молекул STAT5 в цитоплазму играет существенную роль в изучаемом процессе. Выяснены и Глава 8. Модельные уравнения: оценка параметров некоторые детали процесса, которые нельзя наблюдать непосредственно, например, пребывание молекул STAT5 в ядре в среднем около 6 минут.

Исследуя модель (8.9), можно предсказать, что произойдет при изменении тех или иных параметров. Например, авторы исследовали, как изменится суммарное количество произведенного клеткой вещества (оно пропорционально общему количеству STAT5, принявшему участие в процессе) при изменении разных параметров модели. Оказалось, что оно очень слабо зависит от параметров k1,k 2 и сильно – от k 3, k 4,. Т.е.

процессы в ядре клетки играют первостепенную роль. Согласно модели (8.9) уменьшение k 4 до нуля (подавление экспорта из клетки) ведет к уменьшению произведенного вещества на 55 %. В эксперименте с помощью добавления лептомицина B снижали экспорт из клетки (аналог параметра k 4 ) только на 60 %. Согласно модели это должно привести к уменьшению количества активированного STAT5 на 40 %. В эксперименте наблюдалось уменьшение на 37 %. Таким образом, модельный прогноз отлично подтвержден, что еще более повышает доверие к модели и позволяет использовать ее для более детального изучения процесса и управления им в случае необходимости. В связи с этим достижением уже более предметно можно говорить о возможностях использования эмпирических моделей для медицинских целей, упомянутых на с. 239.

Однако следует заметить, что, несмотря на указанные успехи, задача моделирования часто может оказаться технически неразрешимой даже в рассмотренной относительно простой постановке «прозрачного ящика»

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

Глава 9. Модельные уравнения: восстановление нелинейных характеристик В случае только частичной осведомленности о структуре модельных уравнений x n +1 = f (x n, c) или dx dt = f (x, c) задача становится сложнее, чем просто оценка параметров, рассмотренная в предыдущей главе. На наш взгляд, это наиболее плодотворное и практически полезное направление исследований в области моделирования по временным рядам.

В ситуации, которую мы договорились в п. 5.2 называть «серым ящиком», неизвестны некоторые компоненты функции f.

Проиллюстрируем особенность постановки задачи на простом примере.

Пусть объект – нелинейный диссипативный осциллятор, который описывается уравнениями dx1 dt = x2, (9.1) dx2 dt = 0 x2 + F ( x1 ) + A0 cos( 0t + 0 ), где 0, A0, 0, 0 – параметры, F – нелинейная возвращающая сила, вид которой неизвестен. Пусть в качестве наблюдаемого ряда служит временная реализация переменной x1 : = x1. Отметим, что функция F – это только одна компонента, которая задает поле скоростей динамической системы (9.1) наряду с другими членами в правой части. Это характеристика объекта, которая имеет ясный физический смысл и может представлять сама по себе значительный интерес. По условиям эксперимента ее значения могут быть недоступны прямому измерению, т.е.

получить экспериментальные точки на плоскости ( x1, F ) невозможно. Но информация о функции F содержится в наблюдаемом ряде, т.к. F влияет на динамику объекта. «Выделить» значения F можно косвенным путем: через построение эмпирической модели объекта, в структуре которой заложена модельная функция, соответствующая F. А именно, модель следует строить в виде dx1 dt = x2, (9.2) dx2 dt = x2 + f ( x1, c) + Acos(t + ), где f ( x1, c) должна аппроксимировать F.1 Если удастся получить «хорошую» модель (9.2), то будет подтверждена справедливость идей, лежащих в основе этой модели, и одновременно восстановлена нелинейная Аппроксимация функции одного переменного много легче общей задачи аппроксимации функции многих переменных, возникающей обычно в случае «черного ящика» (глава 10).

Часть II. Моделирование по временным рядам характеристика F в виде f ( x1, c). Еще раз подчеркнем, что другого пути измерения характеристики F может не быть.

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

многие общие законы имеют такую форму.

Как и в главе 8, рассматриваемые модели задают зависимости будущего состояния объекта x n+1 от текущего x n или скорости изменения состояния dx dt от самого состояния x. Отличие состоит лишь в том, что перед процедурой оценки параметров надо задать функциональную форму искомых характеристик (например, разложение в каком-либо функциональном базисе, п. 7.2.4), см. п. 9.1. Из-за этой специфики здесь очень важны вопросы оптимизации структуры модели (п. 9.2) или ее специального подбора для избранных объектов (п. 9.3) или классов объектов (п. 9.4).

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

Модельные отображения. Пусть объект – xn+1 = F ( xn ). Вид одномерное отображение функции F неизвестен, наблюдаемая n = xn.

Считая размерность исходной системы известной, будем строить модель в виде одномерного отображения xn+1 = f ( xn, c). Из-за простоты примера в данном случае точки на плоскости ( n, n+1 ) фактически представляют график Рис.9.1. Построение одномерного функции F. Задача – подобрать вид функции модельного f ( x, c) и параметры c так, чтобы она отображения – поиск аппроксимировала эти точки наилучшим образом зависимости (рис.9.1). Подчеркнем, что не каждый отдельный следующего значения параметр, а только вся модельная функция f ( x, c) от предыдущего по экспериментальным теперь может иметь физический смысл.

( ) Мы получили задачу, почти такую же, как в главе 7 (см. рис.7.1,б и Глава 9. Модельные уравнения: восстановление нелинейных характеристик (7.24)). Отличие только в том, что по осям отложены не (t, ), а ( n, n+1 ).

Поэтому и методики используются те же, что в п. 7.2, с заменой пары величин (t, ) на ( n, n+1 ). Так, для аппроксимации функции одной переменной (рис.9.1) удобно использовать алгебраические многочлены невысокого порядка или кубические сплайны, см. п. 7.2.4. Параметры оценивают чаще всего с помощью обычного МНК (8.3).

При наличии динамического шума в объекте ( xn+1 = F ( xn ) + n ) в процедуре построения модели ничего не меняется. Измерительный шум ( n = xn + n ) существенно усложняет задачу оценки параметров (п. 8.1.2).

При низком его уровне еще пригоден обычный МНК (п. 8.1.2.1). При большом уровне желательно было бы использовать более сложные методики (пп. 8.1.2.2, 8.2.1), но в случае «серого ящика» модель, как правило, содержит довольно большое число оцениваемых параметров, что многократно затрудняет использование этих методик.

Модельные дифференциальные уравнения. Для описания сложных, в частности, хаотических, движений обычно используются нелинейные модельные ОДУ dx dt = f (x, c) с числом переменных не менее трех. В рассматриваемой постановке неизвестны некоторые компоненты поля скоростей f, как в примере (9.1) и (9.2). Эти нелинейные характеристики как бы «зашиты» внутри структуры модели и их можно найти, только если построить всю модель. Рассмотрим некоторые детали методики.

Первый случай – наблюдаются все динамические переменные xk :

k (ti ) = xk (ti ) + k (ti ), k = 1,..., D. Для построения модели нужно аппроксимировать зависимость производной xk (t i ) от состояния x(t i ) с & f k (x, c k ) для каждого k = 1,..., D. Значения помощью функции производных x k (t i ) обычно необходимо получить по наблюдаемым & данным k (t i ) с помощью численного дифференцирования (п.7.4.2).

Обозначим полученные оценки xk (t i ). «Сглаженные» значения самих & наблюдаемых, обычно получаемые как «побочный продукт» при дифференцировании, обозначим x k (t i ). По полученным рядам величин xk (ti ), xk (ti ) параметры модельных функций f k (x, c k ), которые включают & в себя и искомые нелинейные характеристики, рассчитываются с помощью обычного МНК:

S (с k ) = ( xk (ti ) f k (x(ti ), с k ) ) min, k = 1,..., D.

(9.3) & i Второй случай – наблюдается единственная переменная:

(ti ) = x1 (ti ) + (ti ). Шансы на успешное моделирование появляются, если уравнения исходного объекта с входящими в них нелинейными характеристиками заданы в стандартном виде (3.27) (п. 3.5.3), где Часть II. Моделирование по временным рядам компоненты вектора состояния x – последовательные производные переменной x1 (ti ). Модель тогда следует строить также в стандартном виде d D x(t ) dt D = f ( x(t ), dx(t ) dt, d D 1 x(t ) dt D 1, c), (9.4) Для этого по наблюдаемой (ti ) путем численного дифференцирования (п. 7.4.2) получают оценки x, x (1), x ( 2),..., x ( D1) и x ( D ) (верхний индекс означает порядок производной). Параметры функции f, которая включает в себя нелинейные характеристики, оценивают с помощью обычного МНК:

( ) S (с) = x ( D ) (ti ) f ( x(ti ), x (1) (ti ),..., x ( D1) (ti ), с) min. (9.5) i В обоих случаях методики эффективны при малом измерительном шуме. Желательно также, чтобы функции f k в (9.3) и f в (9.5) линейно зависели от параметров (псевдолинейные модели, см. п. 7.2.4). При большом шуме моделирование крайне затрудняется, т.к. численное дифференцирование еще усиливает его, особенно при расчете D производных в (9.5). Обычный МНК становится непригодным, а использование более сложных методик в случае многомерных моделей со значительным числом неизвестных параметров нереалистично (см. с. 242).

9.2. Оптимизация структуры модели К задаче о восстановлении характеристик в полной мере относится утверждение о важности выбора структуры модели (пп. 7.2.3, 7.2.4). Для выбора ее размера, например, числа слагаемых многочлена, пригодны критерии, описанные в п. 7.2.3. Но как выбрать функции-слагаемые из большого набора, чтобы обеспечить наилучшую модель заданного размера?

Рассмотрим эффективный подход к подбору функций-слагаемых [25] на примере реконструкции уравнений осциллятора Ван дер Поля – Тоды dx1 dt = x2, dx2 dt = (1 x12 ) x2 1 + e x, (9.6) = x1.

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

Рис.9.2. Фазовая траектория осциллятора Ван дер Поля – Тоды (9.6), включающая переходный процесс. Аттрактор – предельный цикл. Цифры – номера отсчетов, интервал Глава 9. Модельные уравнения: восстановление нелинейных характеристик dx1 dt = x2, (9.7) dx2 dt = f ( x1, x2, c), K ci, j x i y j, i + где f ( x, y, c) = j K. Многие слагаемые этого многочлена i, j = «лишние» – не имеют аналогов в уравнении (9.6): c0,0, c1,1 x1 x2, c0, 2 x2 и др.

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

Их желательно удалить из модельных уравнений.

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

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

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

Часть II. Моделирование по временным рядам Рис.9.3. Результаты реконструкции осциллятора Ван дер Поля – Тоды, начиная с полного многочлена порядка K = 7. а) зависимость коэффициентов при указанных сверху слагаемых многочлена в (9.7) от стартового момента сегмента реконструкции (в единицах интервала выборки);

б) погрешность минимальна при 20 исключенных слагаемых Будем строить модель (9.7) с многочленом высокого порядка K по W:

последовательным сегментам временного ряда длиной { ( k 1)W +1,..., ( k 1)W +W }, k = 1, 2,..., L. Получим набор оценок коэффициентов сi(,kj) (рис.9.3,а). Степень стабильности оценки каждого коэффициента сi, j определим как модуль отношения эмпирического среднего (1 L ) (ck, j )2.

L L сi, j = (1 L ) сi(,kj) к стандартному отклонению ck k =1 j = Исключим слагаемое, соответствующее наименее стабильному коэффициенту. Для упрощенной структуры модели повторим всю процедуру. Продолжая ее дальше, будем последовательно удалять «нестабильные слагаемые».

Исключение слагаемых прекращается, когда качество модели перестает улучшаться. На рис.9.3,б мы использовали как критерий качества погрешность описания объекта в широкой области V фазового [ ] 2 = {f ( x1, x2, c) (1 x12 ) x2 1 + e x } dx1dx2.

пространства После V исключения 20 слагаемых из модельного многочлена порядка K = погрешность уменьшается на порядок по сравнению со стартовой структурой модели [190].

9.3. Восстановление характеристик нелинейного элемента электрической цепи Рассмотрим показательный пример из области радиофизики, когда хаотическое поведение достаточно простой нелинейной системы успешно моделируется в постановке «серый ящик». Объект – RLC-контур с переключаемыми конденсаторами, находящийся под гармоническим внешним воздействием амплитуды U 0 и частоты 0. Его схема представлена на рис.9.4,а. K – электронный ключ: микросхема, содержащая десятки транзисторов и других пассивных элементов, которая питается от специального источника постоянного напряжения. При малых значениях напряжения U на емкости C1 происходят линейные колебания в контуре RLC1 (сопротивление ключа очень велико), когда напряжение U достигает порогового значения U пор, сопротивление ключа резко падает, он замыкает цепь и подключает емкость C 2. Обратное переключение Глава 9. Модельные уравнения: восстановление нелинейных характеристик происходит при понижении U относительно U пор (ключ обладает гистерезисом). Наличие нелинейности приводит к возможности хаотических колебаний в контуре.

При U пор 0 модель этой системы, полученная на основе законов Кирхгофа, примет вид уравнения неавтономного нелинейного осциллятора (9.1), где безразмерные переменные t = t / LC1, x = q / C 2 U пор (t – время, q – суммарный заряд на емкостях C1 и C2 ), а F – кусочно-линейная функция, связанная с емкостями, т.е. с вольт-фарадной характеристикой нелинейного элемента (переключающихся конденсаторов).

В эксперименте измерялся хаотический временной ряд силы тока через резистор R, это величина dx dt в уравнениях (9.1). С помощью численного интегрирования наблюдаемого ряда был получен временной ряд значений x, с помощью численного дифференцирования – оценки d 2 x dt 2. Будем считать, что информации о кусочно-линейном виде функции F нет, тем более что это лишь приближение, игнорирующее гистерезис при переключении ключа и т.п. Модели строились в виде (9.2) с многочленом f различных порядков K. На рис.9.4,в представлены результаты для наилучшей модели с K = 9, которая хорошо воспроизводит наблюдаемое хаотическое поведение (рис.9.4,б). Теоретическая кусочно линейная «возвращающая сила» и модельный многочлен f с хорошей точностью совпадают в области наблюдаемого движения, показанной штриховой линией на рис.9.4,г. Отметим, что без информации о структуре уравнений (9.2) получить адекватную эмпирическую модель с физическим смыслом невозможно [34].

Рис.9.4. Моделирование радиофизической системы: а) схема экспериментальной установки;

б) экспериментальная хаотическая траектория на плоскости заряд – сила тока, измерения I проводились с помощью 12-разрядного АЦП при t = 4 мкс, C1 = 0. мкФ, C2 = 4.4 мкФ, L = 0.02 Гн, R = 10 Ом, Uпор = -0.2 В, периоде воздействия T 84.02, U0 = 2.344 В (т.е. 0.02, 1.02, C1/(C1+C2) 1/45), единицы по осям безразмерные [34];

в) траектория восстановленной модели (9.2) с многочленом f 9-го порядка, г) график f с обратным знаком (жирная линия) и ожидаемая кусочно-линейная зависимость (тонкая) Часть II. Моделирование по временным рядам Рассмотренный простой пример иллюстрирует возможность определения характеристик элемента нелинейной системы с помощью процедуры моделирования по временному ряду, причем даже в режимах больших амплитуд и хаоса, когда эти характеристики могут быть недоступны измерениям с помощью обычных средств. Он успешно применялся для исследования динамических характеристик конденсатора с сегнетоэлектриком [236], полупроводниковых диодов [322], волоконно оптических систем [333]. Подробнее познакомиться с объектом и процедурой моделирования на примерах эталонных теоретических осцилляторов, контура с переключаемыми конденсаторами и контура с полупроводниковым диодом с p-n переходом можно в лабораторных работах [30].

9.4. Специфические подходы к выбору структуры модели Неопределенность в структуре модели может быть не столь мала, как в примерах, описанных выше. «Ящик» может быть «темно-серым», а не «светло-серым» (п. 5.2), что сильно усложняет моделирование. Но в ряде случаев даже небольшой объем априорной информации о свойствах системы или результаты предварительного анализа рядов наблюдаемых могут привести к успеху моделирования при их учете в структуре модели.

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

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

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

Если предположение о наличии воздействия сформировано, то в уравнения модели полезно включить функции, явно зависящие от времени, для описания этого воздействия. Например, для описания гармонического аддитивного воздействия целесообразна структура модели [34, 35, 194]:

d D x dt D = f ( x, dx dt,..., d D 1 x dt D 1, c) + a cos t + b sin t, (9.8) где x – наблюдаемая, f – алгебраический многочлен. В этом случае достаточно использовать меньшее число переменных D, чем было бы необходимо при использовании стандартной структуры (9.4), что и дает преимущества при использовании специальной структуры (9.8), см. также Глава 9. Модельные уравнения: восстановление нелинейных характеристик п. 10.2.2. Уравнение осциллятора (9.2) – частный случай (9.8) при D = 2 и неполном многочлене от двух переменных.

Отметим, что для успеха моделирования кроме выбора структуры уравнений необходимо еще преодолеть специфическую техническую проблему. Это проблема оценки частоты воздействия, которая входит как нелинейный параметр в уравнения модели (9.8). Как обычно, задают стартовую догадку и решают задачу минимизации целевой функции типа (9.5) итерационным методом. Но особенность здесь состоит в том, что правая часть модельных уравнений (9.8) очень чувствительна к значению при больших t, аналогично примеру (7.19) в п. 7.1.2.2. Поэтому и целевая функция S типа (9.5) чувствительна к при большой длине ряда N. Это ведет к тому, что в случае отыскания глобального минимума S дисперсия оценки быстро уменьшается с ростом длины ряда: как N 3, аналогично примеру (7.19). С одной стороны, это дает возможность очень точно определить. С другой стороны, возникают трудности отыскания глобального минимума, т.к. требуется очень удачная стартовая догадка для. Учитывая это, следует тщательно перебирать различные стартовые догадки.

Если же по условиям эксперимента известна заранее, но с погрешностью, то важно помнить, что в случае очень длинного ряда малая погрешность может привести к плохому описанию «истинного»

воздействия соответствующими слагаемыми в (9.8) из-за «набега разности фаз» между ними [34] и к бесполезности использования структуры (9.8).

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

В случае произвольного регулярного воздействия (сложного периодического и квазипериодического) более адекватна структура модели d D x dt D = f ( x, dx dt,..., d D 1 x dt D1, c) + g (t, c), (9.9) где функция g(t) описывает воздействие и может иметь вид тригонометрического многочлена [159]:

k Ki g (t ) = ci, j cos( j i t + i, j ). (9.10) i =1 j = Причем адекватные модели с тригонометрическими многочленами могут быть получены и при очень большом числе (сотни) используемых гармоник K i, тогда как при аппроксимации алгебраическими многочленами увеличение их порядка чревато неустойчивостью траекторий модели.

Часть II. Моделирование по временным рядам Чтобы учесть и мультипликативное, и более сложное воздействие, зависимость от времени можно вводить во все коэффициенты многочлена f в (9.8) [194]. Так, на рис.9.5 показан пример реконструкции неавтономного осциллятора Тоды при комбинированном гармоническом воздействии d 2 x dt 2 = 0.45 dx dt + (5 + 4 cos t )(e x 1) + 7 sin t. (9.11) Рис.9.5. Реконструкция уравнений нелинейного осциллятора Тоды (9.11) по ряду переменной x: а) аттрактор объекта, б) аттрактор неавтономной модели с зависимостью от времени во всех коэффициентах многочлена (D = 2, K = 9, см. текст), в) траектория стандартной модели (9.4) с D = 4, K = Модель строилась в виде (9.4), где воздействие вводилось во все коэффициенты многочлена f, т.е. вместо сk подставлялись выражения сk + ak cost + bk cost. Лучшая модель получена при D = 2, K = 9, ее фазовая траектория очень схожа с экспериментальной (рис.9.5,а,б).

Типичная неустойчивая траектория эмпирической стандартной модели (9.4), размерность которой должна быть не менее 3 для описания хаотического режима, показана на рис.9.5,в.

Эффективность всех упомянутых подходов была показана на численных примерах реконструкции уравнений по зашумленным хаотическим реализациям эталонных осцилляторов при различных видах воздействия: периодическом импульсном, периодическом с субгармониками, квазипериодическом [159, 193, 194].

9.4.2. Системы с запаздыванием Методы построения моделей систем с запаздыванием активно развиваются в последнее время [202, 203, 331, 332, 241, 191, 290, 139, 140].

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

Рассмотрим пример, когда моделирование системы с запаздыванием попадает под раздел «серый ящик» и может быть проведено с помощью подходов, близких к описанным выше. Речь пойдет о системах вида:

Глава 9. Модельные уравнения: восстановление нелинейных характеристик 0 x(t ) = x(t ) + F ( x(t )), (9.12) & где наблюдаемая = x, возможно с измерительным шумом.

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

0 x(t ) = x(t ) + µ 0 sin( x(t 0 ) x0 ). (9.13) & Рис.9.6. a) Временная реализация уравнения Икеды (9.13), 0 = 1.0, µ 0 = 20.0, 0 = 2.0.

б) Число пар экстремумов M(), нормированное на общее число экстремумов ряда.

Mmin() = M(2.0). в) Восстановленная нелинейная функция. Численные эксперименты с добавлением измерительного шума показывают возможность реконструкции при уровнях шума до 20% (по отношению среднеквадратичных амплитуд шума и сигнала) Модели строятся в виде:

x(t ) = x(t ) + f ( x(t ), c). (9.14) & Аналогично выше рассмотренным примерам, можно решать задачу ( x(t n ) + x(t n ) f ( x(t n ), c) ) min, рассматривая параметр & n инерционности и время задержки как дополнительные неизвестные параметры [203]. Но есть и специальный, более эффективный подход [191], который состоит в следующем.

Статистический анализ временных интервалов, разделяющих экстремумы во временных реализациях различных модельных и реальных систем с запаздыванием (9.12), позволяет установить, что зависимость числа M пар экстремумов временной реализации, удаленных друг от друга на время, от величины, имеет четкий минимум при, соответствующем времени запаздывания системы, рис.9.6,б. Это наблюдение дает способ оценки времени запаздывания и диагностический признак принадлежности объекта к системам этого класса. Имея оценку 0, можно найти оценку параметра инерционности и аппроксимирующую функцию f путем перебора пробных значений и выбора такого, при котором экспериментальные точки на плоскости ( x(t ), x(t ) + x(t ) ) наилучшим & образом ложатся на некоторую гладкую кривую. Последняя и представляет собой график искомой функции f. Рис.9.6,в иллюстрирует Часть II. Моделирование по временным рядам такое восстановление «истинной» функции F (x) для системы (9.13). По полученному графику можно найти и формулу для аппроксимирующей функции f в некотором функциональном базисе или в специальном виде.

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

Глава 10. Реконструкция уравнений: «черный ящик»

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

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

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

Его теоретическому обоснованию посвящены известные теоремы Такенса, которые упоминаются по необходимости и по традиции весьма часто (п. 10.1).

Не менее важен и труден в постановке «черного ящика» следующий этап – аппроксимация зависимости следующего значения вектора состояния от текущего (при построении модельных отображений x n +1 = f ( x n, c) ) или скорости изменения вектора состояния от самого вектора (при построении модельных дифференциальных уравнений dx dt = f ( x, c) ). На практике обычно удается добиться успеха, если оказывается достаточным использование не очень большой размерности модели, грубо говоря, не более 5-6. Для построения моделей большей размерности требуются огромные объемы экспериментальных данных, а аппроксимация функций многих переменных намного сложнее приближения одномерных зависимостей (пп. 7.2, 9.1, 9.3). Причем трудности быстро нарастают с ростом размерности модели – это так называемое «проклятие размерности» [254], основное препятствие при моделировании множества реальных процессов.


Говорят также о «реконструкции фазовой траектории» и «реконструкции векторов состояния».

Глава 10. Реконструкция уравнений: «черный ящик»

Чтобы развеять пессимистические нотки предыдущих абзацев, отметим, что и в такой тяжелой постановке, как «черный ящик», тоже удается добиться успеха при моделировании реальных сложных объектов.

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

10.1. Реконструкция фазовой траектории В соответствии со сказанным в п. 6.1.2, при построении эмпирических моделей по временному ряду { (t1 ), (t 2 ),..., (t N )} можно использовать в качестве недостающих переменных последовательные значения наблюдаемой (вектор состояния модели x(ti ) = [ (ti ), (ti + ),..., (ti + ( D 1) )], где – время задержки), ее последовательные производные (вектор состояния x(ti ) = [ (ti ), d (ti ) dt,..., d D 1 (ti ) dt D 1 ] ) и т.д. Эти подходы часто применялись с давних пор даже без специального обоснования. Так, первый из них используется с 1927 года [338] при построении широко известных авторегрессионных моделей (см. (4.12) в п. 4.4), где будущее значение наблюдаемой предсказывается по нескольким предыдущим. Это выглядит вполне разумно: если нет другой информации, кроме наблюдаемого ряда, то на что же еще опираться при прогнозе, как не на предыдущие значения наблюдаемой или какие-либо их комбинации?

В начале 1980-х годов была показана связь обоих этих эмпирических подходов с теорией динамических систем. Было доказано, что при реконструкции по скалярной временной реализации динамической системы2 и метод временных задержек, и метод последовательных производных гарантируют, что в новых переменных будет получено эквивалентное3 описание исходной динамики при достаточно большой размерности восстановленных векторов D. А именно, должно выполняться условие D 2d, где d – размерность множества M в фазовом пространстве исходной системы, на котором происходит моделируемое движение.4 Эти При некоторых условиях гладкости, см. ниже п. 10.1.1.2.

Точнее, диффеоморфное описание, п. 10.1.1.2.

Множество M – компактное гладкое многообразие, а величина d – его топологическая размерность, см. п. 10.1.1.2. Есть обобщения теоремы на случай негладких множеств M и фрактальной размерности d, которые мы здесь не обсуждаем [299]. Заметим, что множество M, о котором идет речь в теоремах, не обязательно соответствует движению на аттракторе. Например, пусть исследуемая система имеет аттрактор – предельный цикл C, который «намотан» на тор M. Если нас интересует описание только установившегося периодического движения на цикле C, то, согласно теореме Часть II. Моделирование по временным рядам утверждения составляют содержание знаменитых теорем Такенса [323], см. п. 10.1.1.

Подчеркнем, что теоремы относятся только к случаю, когда исходный объект – динамическая система (см. п. 2.2.1, с. 44). При моделировании реальных объектов использование упомянутых методов возможно и без упоминания о теоремах (как и делалось на протяжении не одного десятка лет), т.к. на практике нельзя проверить, выполняются ли условия теорем,5 а размерность d неизвестна. Но значение теоретических результатов Такенса все же достаточно велико. Во-первых, после них стало ясно, что для достаточно широкого класса систем упомянутые методы с гарантией дают недостающие переменные, пригодные для построения динамической модели. Так что теоремы «благословляют» практическое применение методов, особенно если есть какие-либо соображения в поддержку того, что условия теорем для данной ситуации выполняются. Во-вторых, с опорой на идеи теории динамических систем впоследствии были развиты новые полезные подходы к подбору параметров алгоритма реконструкции – времени задержки (например, первый минимум функции взаимной информации), размерности модели D (например, метод ложных ближайших соседей) и т.д., см. п. 10.1.2.

10.1.1. Теоремы Такенса Поясним сначала нестрого содержание теорем на простом примере (п. 10.1.1.1), а затем приведем их математическую формулировку и обсудим некоторые детали на более строгом языке (п. 10.1.1.2). В данной главе мы будем обозначать вектор состояния исходной системы y в отличие от восстановленных векторов – x. Обозначение d относится всегда к размерности множества M в фазовом пространстве исходной системы (это не обязательно размерность всего пространства, т.е. векторов y), а D – это размерность восстановленных векторов x и, следовательно, модели.

10.1.1.1. Пояснительный пример. Пусть объект представляет собой трехмерную динамическую систему с непрерывным временем. Вектор состояния y = ( y1, y 2, y3 ). Пусть движение происходит на предельном цикле (рис.10.1,а), т.е. на множестве M размерности d = 1.

Такенса, для реконструкции достаточно использовать D 2 модельных переменных.

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

Насколько вообще можно говорить о выполнении условий математической теоремы применительно к реальному объекту.

Глава 10. Реконструкция уравнений: «черный ящик»

Если бы наблюдались все три переменных y1, y 2, y3 (все компоненты вектора y), то можно было бы сразу приступать к аппроксимации зависимости следующего состояния y (t + ) от текущего y (t ). Эта зависимость однозначна, т.к. y – вектор состояния. Последнее означает, что каждый раз, когда наблюдается некоторое значение вектора y = y *, через фиксированный интервал времени каждый раз будет наблюдаться одно и то же событие (одно и то же будущее). Если же наблюдаются не все компоненты вектора состояния, то ситуация может быть сложнее.

Поставим вопрос: сколько и каких переменных достаточно использовать для эквивалентного описания исходной динамики, а какие варианты не подходят для этого?

Поскольку множество M, на котором происходит рассматриваемое движение, одномерно (d = 1), то должна существовать такая скалярная динамическая переменная, с помощью которой можно описать это движение. Например, замкнутую кривую M (рис.10.1,а) можно отобразить на окружность (рис.10.1,б). Важно, что при этом можно взаимно однозначно связать векторы y (t ) на цикле M с углом Рис.10.1. Одномерное представление динамики на предельном цикле (а) с поворота точки на окружности помощью отображения на окружность (б) (t ).6 Благодаря взаимной и проекции на координатную ось (в).

однозначности, переменная Размерность фазового пространства полностью задает состояние исходной системы равна 3, размерность многообразия M – замкнутой кривой d = 1, системы: значение фазы * всегда размерность «восстановленных» векторов соответствует одному и тому же в обоих случаях D = 1. Два разных значению вектора y * в тот же исходных состояния (черные кружки на рисунке а) соответствуют двум разным момент времени. Располагая точкам на окружности (рисунок б) и одной наблюдаемой, можно строить и той же точке на оси (рисунок в).

одномерную ( D = 1 ) динамическую Отображение цикла на окружность взаимно однозначно, а на отрезок – нет модель с x1 =.

Но далеко не каждая переменная подходит для одномерного динамического описания. Так, если наблюдается лишь одна из компонент вектора y, например y1, то замкнутая кривая отображается на отрезок (это Эта переменная – «свернутая» фаза колебаний, рассматривавшаяся в п. 6.4.3.

Часть II. Моделирование по временным рядам просто проекция на ось y1 ). Это отображение не взаимно однозначно:

* каждой точке отрезка с координатой y1 (t ) соответствуют, по меньшей мере, два вектора состояния y (t ), отличающиеся направлением дальнейшего движения – в сторону увеличения или уменьшения y1, т.е.

влево или вправо на оси y1 (t ), см. рис.10.1,а,в. Таким образом, задание y не определяет однозначно состояния системы. Если наблюдается * некоторое значение y1 = y1, то через фиксированный интервал времени может реализоваться один из двух (или более) вариантов будущего.

Поэтому одномерное описание наблюдаемого движения при использовании переменной x1 = y1 невозможно.

Итак, при использовании размерности модели D, равной размерности наблюдаемого движения d, построение динамической модели может оказаться как возможным (если «повезло»), так и невозможным. Оба варианта типичны.7 Это утверждение мы проиллюстрировали на простом примере, но оно справедливо и в общем случае.

Что меняется если использовать двухмерное описание? Здесь оказываются возможными и типичными те же две ситуации. Они проиллюстрированы на рис.10.2. Рассмотрим в качестве наблюдаемых две компоненты исходного вектора состояния – y1 и y 2, т.е. модельный вектор x = ( y1, y 2 ). Это соответствует проекции замкнутой кривой на плоскость ( y1, y 2 ). В проекции можно наблюдать замкнутую кривую без самопересечений (рис.10.2,а) или с самопересечениями (рис.10.2,б), в зависимости от Рис.10.2. Два случая при проецировании формы замкнутой кривой и ее одномерного многообразия из трехмерного ориентации в пространстве. В пространства в двухмерное: а) взаимно первом случае связь исходной однозначное отображение, б) не взаимно замкнутой кривой и ее однозначное (в проекции появилась точка проекции взаимно однозначна, самопересечения кривой) т.е. двухмерные векторы x полностью определяют состояние и могут использоваться для построения Типичность здесь и ниже понимается в смысле «грубости», «случая общего положения», т.е. ситуация не меняется при «малом шевелении» исходной системы и наблюдаемой. Более строгую формулировку см. в п. 10.1.1.2.

Глава 10. Реконструкция уравнений: «черный ящик»


динамической модели. Во втором случае это не так. Точка * * самопересечения y1, y 2 на плоскости ( y1, y 2 ), рис.10.2,б, соответствует двум различным состояниям исходной системы (т.е. связь x и y не взаимно * * однозначна). Поэтому, когда наблюдаются значения y1, y 2, то будущее по ним однозначно предсказать нельзя. Значит, с помощью такого вектора x не удастся осуществить динамическое описание процесса. Аналогичная ситуация имеет место, если в качестве двух переменных модели используются не y1 и y 2, а две переменных, полученных любым из упомянутых выше методов. Пусть наблюдаемая = h(y ), где h – произвольная гладкая (измерительная) функция, а компоненты вектора x получены методом временных задержек x(t ) = ( (t ), (t + )). При этом, в зависимости от h и, может иметь место ситуация, аналогичная рис.10.2,а, – нет самопересечений на плоскости ( x1, x2 ), но столь же типична и ситуация рис.10.2,б, когда динамическое моделирование невозможно.

Таким образом, даже число переменных большее, чем размерность исходного движения, D = 2 d = 1, не гарантирует возможности динамического описания.

Наконец, рассмотрим трехмерное представление. Например, когда для = h(y ) наблюдаемой векторы восстанавливаются в виде x(t ) = ( (t ), (t + ), (t + 2 )). Образ исходной замкнутой кривой в трехмерном пространстве ( x1, x2, x3 ) – также замкнутая кривая. И в типичном случае она тоже не имеет самопересечений. Это означает, что имеет место взаимно однозначное соответствие между векторами x и y.

Исходное движение на предельном цикле может быть эквивалентно описано с помощью векторов x. Вообще говоря, самопересечение кривой в пространстве x1, x2, x3 может иметь место, но это ситуация вырожденная, нетипичная. Интуитивно не трудно согласиться с тем, что самопересечение кривой в трехмерном пространстве, грубо говоря, крайне маловероятно.

Итак, в рассмотренном примере эквивалентное описание динамики достигается с гарантией только при восстановлении в пространстве размерности D 2d.9 Это и есть основное утверждение теорем Такенса.

Подчеркнем, что это – достаточное условие. В отдельных случаях может оказаться возможным эквивалентное описание даже при D = d (при удачной наблюдаемой). При реконструкции уравнений на практике теорема Такенса служит лишь для дополнительной психологической Это возможно только локально – в удалении от точки самопересечения.

Эквивалентное описание движения на предельном цикле гарантируется при D = 3, даже если этот цикл «живет» в бесконечномерном фазовом пространстве.

Часть II. Моделирование по временным рядам уверенности, ведь согласно ей существует конечная размерность модели D, при которой вероятен успех динамического описания! Реально же приходится пробовать различные значения D, начиная с малых, и стараться получить «хорошую» модель как можно меньшей размерности, чтобы избежать трудностей, связанных с упомянутым выше «проклятием размерности».

10.1.1.2. Математические детали. Сформулируем теоремы более строго. Для этого потребуется ввести ряд терминов и обозначений. Пусть объект представляет собой динамическую систему:

y (t 0 + t ) = t (y (t 0 )), (10.1) (t ) = h(y (t )), где y – вектор состояния, t – оператор эволюции, h – измерительная функция.10 Вектор наблюдаемых конечномерен R m, причем мы будем говорить только о наиболее распространенном (и наиболее сложном для моделирования) случае скалярного временного ряда (ti ), m = 1.

Многообразие. Пусть движение системы происходит на многообразии M конечной размерности d.11 Многообразие – это обобщение понятия гладкой поверхности в евклидовом пространстве (подробности см. в [58, 112, 299] и [115, с. 244-248]). Грубо, d-мерное многообразие M – это поверхность, которую локально в окрестности каждой точки можно параметризовать с помощью d евклидовых координат. Другими словами, любую точку p M и ее локальную окрестность U ( p ) можно взаимно однозначно и взаимно непрерывно отобразить на d-мерный фрагмент (например, шар) пространства R d. Полученный образ : U (U ) называют картой окрестности, а само непрерывное отображение – гомеоморфизмом. Примерами двухмерных многообразий в трехмерном Если объект – отображение y (t n +1 ) = F(y (t n )), то оператор эволюции t (y (t 0 )) есть просто функция F. Если объект – система ОДУ dy dt = F ( y (t )), то функция t (y (t 0 )) – результат интегрирования ОДУ на интервале шириной t. Если исходная система – уравнение в частных производных y t = F ( y, y r, 2 y r 2,...), где r – пространственная координата, то y есть вектор бесконечномерного пространства функций y(r), а t – оператор, действующий в этом пространстве.

Это может иметь место даже для бесконечномерной системы. Для описания многомерных систем вводится полезное понятие инерциального многообразия [115]. Не строго, это многообразие минимальной размерности d, каждая точка которого полностью определяет состояние исходной системы. Величины, определяющие положение точек на многообразии (d штук), называют параметрами порядка.

Глава 10. Реконструкция уравнений: «черный ящик»

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

(рис.10.3,а), но не (двойной) конус (рис.10.3,б).

Рис.10.3. Примеры множеств, которые являются (а) и не являются (б) многообразиями Если – n раз дифференцируемое отображение с таким же обратным, то говорят, что M принадлежит классу С n. В случае n отображение называют диффеоморфизмом. Если многообразие M отображается с помощью диффеоморфизма на многообразие S R D, D d, то M и S называют диффеоморфными друг другу. При этом говорят, что S – вложение многообразия M в евклидово пространство R D.

Ниже потребуется еще, чтобы многообразие M было ограничено и замкнуто. Первое условие означает, что его можно заключить в шар конечного радиуса, а второе – то, что все предельные точки многообразия M принадлежат ему же. Такое многообразие M в конечномерном пространстве называют компактным.

Постановка вопроса и обозначения. Каждой фазовой траектории системы (10.1) y (t ), 0 t, на многообразии M соответствует временная реализация наблюдаемой величины : (t ) = h(y (t )), 0 t.

Вектор y (t 0 ) однозначно определяет будущее поведение системы (10.1), в частности, всю реализацию (t ), t t 0. Можно ли по фрагменту временной реализации (t ) в окрестности момента t 0 однозначно определить положение системы на многообразии M в момент t 0 и, следовательно, всю будущую эволюцию? Другими словами, можно ли по значениям (t ) на ограниченном интервале времени «восстановить» состояние системы? Это ключевой вопрос, и теоремы Такенса дают положительный ответ на него при выполнении ряда условий.

Сначала введем необходимые обозначения для строгого рассмотрению метода задержек. Поставим в соответствие вектору y(t) D мерный вектор x(t ) = [ (t ), (t + ),..., (t + ( D 1) )]. Зависимость x от y (имеются в виду значения векторов в один и тот же момент времени t) описывается однозначным отображением : M R D, которое выражается через оператор эволюции t : M M и функцию h : M R следующим образом:

Часть II. Моделирование по временным рядам 1 (y (t )) h(y (t )) h( (y (t ))) (y (t )) =.

x(t ) = (y (t )) 2 (10.2)......

h ( D (y (t )) ( D 1) ( y (t ))) Гладкость (непрерывность, дифференцируемость, существование и непрерывность производных высокого порядка) определяется гладкостью и h. Образом многообразия M под действием является некоторое множество S R D.

Поставленный выше вопрос формулируется теперь так: является ли диффеоморфизмом? Если да, то S – это вложение M и любому вектору x из S соответствует только один вектор y из M.12 В таком случае x(t) может служить вектором состояния для описания динамики на M и уравнения (10.1) можно переписать в терминах новых переменных:

x(t 0 + t ) = t (x(t 0 )), (10.3) где новый оператор эволюции t (x) = ( t ( 1 (x))). Причем в силу диффеоморфизма локальные свойства динамики (устойчивость и типы неподвижных точек и т.д.) сохраняются. Любой фазовой траектории y(t) на M взаимно однозначно соответствует траектория x(t) на S. Если система (10.1) имеет аттрактор в M, то система (10.3) имеет аттрактор в S;

такие характеристики, как фрактальная размерность и ляпуновские показатели для этих аттракторов одинаковы. Другими словами, система (10.3) на многообразии M и система (10.1) на многообразии S – это два «почти одинаковых» представления одной и той же динамической системы.

Ясно, что отображение (10.2) является диффеоморфизмом не всегда. Так, на рис.10.2,б представлен пример, когда гладкое отображение (проецирование) одномерного многообразия M имеет неоднозначное обратное отображение (рис.10.2,б). Возможна и другая нежелательная ситуация, представленная на рис.10.4: Другими словами, если в разные моменты времени t встречается один и тот же фрагмент ряда [ (t ), (t + ),..., (t + (d 1) ) ], то за ним всегда следует одно и то же будущее. Это дало бы обоснование прогностическому методу аналогов, Рис.10.4. При применявшемуся еще Э. Лоренцем. Метод основан на поиске фрагментов временного проецировании ряда в прошлом, которые «похожи» на текущий (с точностью до конечной величины одномерного многообразия метрике), и использовании какой-либо комбинации следовавших за в некоторой M ними фрагментов место может иметь в качестве прогноза. В более современной трактовке – это локальные модели, см. п.10.2.1.4.

Глава 10. Реконструкция уравнений: «черный ящик»

однозначно, но не дифференцируемо. Последнее имеет место в точке возврата на множестве S (в проекции), в окрестности которой двухмерный вектор ( y1, y 2 ) не может использоваться для описания динамики с помощью ОДУ, т.к. точка возврата оказывается неподвижной точкой, и кривая S не может быть предельным циклом. Дифференциальные свойства M и S оказываются различными из-за недифференцируемости отображения 1.

Формулировка теоремы. Возвращаясь к системе (10.1) и отображению (10.2), можно сказать, что для разных, M, h, d, D и можно встретиться с любой из упомянутых ситуаций. Множества самопересечений и множества возврата на S = (M) могут быть весьма обширными, что крайне нежелательно. Но можно встретить и «хорошую» ситуацию вложения (рис.10.2,а). Приведенный ниже результат, впервые был строго получен голландским математиком Флорисом Такенсом [323], а затем обобщен в работе [299]. Он говорит о том, при каких условиях обеспечивается вложение исходного компактного d-мерного многообразия M в пространство R D с помощью отображения (10.2). Теорема 1. Пусть M – компактное d-мерное многообразие класса C 2.

Для почти любой 14 пары t, h (где t и h дважды непрерывно дифференцируемы на M) отображение : M R D, определяемое (10.2), есть диффеоморфизм15 для почти любого16 0 и D 2d.

Обсуждение. Итак, если взять достаточно большую размерность векторов временных задержек x (10.2), то получим (за исключением вырожденных случаев) вложение многообразия M и сможем использовать x в качестве векторов состояния динамической модели. Возможна следующая наглядная интерпретация того, что размерность D должна быть больше именно величины 2d [115]. Чтобы установить неоднозначность Теорема Такенса тесно связана с теоремой Уитни, которая касается произвольных отображений, из курсов дифференциальной геометрии. Первая отличается тем, что относится к специальному случаю отображений (10.2), определяющихся оператором эволюции динамической системы [299].

Термин «почти любая пара» понимается у Такенса в смысле типичности. Например, если для какого-либо t отображение (10.2) не дает вложения, то найдется такая сколь угодно малая вариация t + t, что вложение будет достигнуто. Более строго, типичные свойства выполняются на пересечении открытых и всюду плотных множеств. Метрическим аналогом типичности является превалентность (prevalence) [299].

Следовательно, отображение (10.2) дает вложение многообразия M. Пространство D R, в котором содержится образ S = (M) называют пространством вложения.

Например, в случае наличия предельного цикла внутри M величина не должна быть равна периоду колебаний на этом цикле. Подробнее см. [299].

Часть II. Моделирование по временным рядам отображения 1, надо на многообразии M найти такие два вектора y 1 и y 2, для которых ( y 1 ) = ( y 2 ). Последнее равенство представляет собой систему D уравнений с 2d переменными (по d компонент каждого вектора y 1 и y 2, задающих положение на M). Грубо, эта система не имеет решений в типичном случае, если число уравнений больше числа неизвестных, т.е.

D 2d. Это и есть утверждение теоремы Такенса.

Еще раз подчеркнем, что условие D 2d, наложенное на размерность вектора x, – достаточное, но не необходимое. Т.е. при соблюдении условия диффеоморфизм гарантирован с точностью до типичности, но если «повезет», то можно получить «хорошую реконструкцию» и при меньшей размерности, см. рис.10.1,а,б, где вложение одномерного многообразия достигнуто при D = 1 и является случаем общего положения.

Какими могут быть те нетипичные случаи, когда теоремы не справедливы? Упомянем два примера [115].

1) Измерительная функция – константа: h(y ) = a. Это гладкая функция, но она отображает всю динамику в одну точку. Малым шевелением измерительной функции это устраняется – добавкой к a «малой» функции от y (конечно, не просто константы).

2) Система, состоящая из двух подсистем с однонаправленной связью dy 1 dt = F ( y 1, y 2 ), dy 2 dt = G ( y 2 ), при наблюдении только за ведущей подсистемой – = h( y 2 ). В наблюдаемой нет информации об y1, поэтому вложение невозможно. Устраняется добавкой зависимости от y1 к.

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

(t ) d (t ) dt x(t ) =, (10.4)...

D 1 D d (t ) dt где D 2d. Она формулируется полностью аналогично теореме 1, но с более жесткими требованиями к гладкости функций t и h. А именно, требуется существование непрерывных производных D-го порядка каждой из них, чтобы существовали производные, используемые в (10.4).

Глава 10. Реконструкция уравнений: «черный ящик»

Если производные аппроксимировать конечными разностями, то (10.4) сводится к фильтрованному вложению [229].

На практике всегда есть шумы, а к такому случаю теоремы Такенса не относятся, хотя имеются некоторые их обобщения [319, 207]. Тем не менее, теоремы представляют ценность для решения практических задач моделирования, о чем уже говорилось в начале п. 10.1 (см. с. 256).

10.1.2. Практические алгоритмы реконструкции 10.1.2.1. Метод временных задержек. Это наиболее популярный из методов реконструкции. Из скалярного ряда { i }i =1 формируется ряд N {x i = (i,i+l,...,i+( D1)l )}iN=1( D1)l. = lt векторов Величина задержки теоретически может быть почти любой, но на практике часто нежелательны как слишком малые l, которые ведут к сильной корреляции компонент17 вектора состояния, так и слишком большие l, при которых может слишком усложниться структура восстановленного аттрактора.

Поэтому предлагалось выбирать величину, равной первому нулю автокорреляционной функции [229], первому минимуму функции взаимной информации [227], и пр. [269]. Используется неравномерное вложение, когда временные интервалы между компонентами вектора состояния не одинаковы, что уместно для систем с несколькими характерными временными масштабами [220, 249]. Для неоднородной динамики (например, с чередующимися интервалами почти периодического и очень сложного движений) предложено переменное вложение, при котором набор временных задержек и даже размерность вектора x зависят от его положения в пространстве состояний [249].

Каждый из подходов хорош для своего круга систем, нередко ведет к сложным вычислениям и не гарантирует успеха в общем случае [115].

Как выбрать размерность модели D, опираясь на анализ наблюдаемого ряда? Есть различные методы оценки: метод ложных соседей [258], метод главных компонент [199], метод Грассбергера – Прокаччиа [234], метод хорошо приспособленного базиса [108]. Два первых мы рассмотрим ниже (пп. 10.1.2.2, 10.1.2.3). Однако часто приходится просто перебирать разные пробные значения размерности D, начиная с малых и постепенно увеличивая, и строить модельные уравнения для каждого D, пока не будет получена «хорошая» модель. В таком случае подбор размерности и даже временных задержек может становиться частью единого процесса моделирования, а не отдельной замкнутой первой стадией.

При малом интервале выборки и l = 1 восстановленная фазовая траектория вытягивается вдоль главной диагонали x1 = x 2 =... = x D (практически совпадает с ней).

Часть II. Моделирование по временным рядам 10.1.2.2. Метод ложных соседей. Это метод целочисленной оценки размерности аттрактора сверху. Он основан на проверке того свойства, что фазовая траектория, восстановленная в пространстве достаточной размерности не должна иметь самопересечений. Проиллюстрируем метод на простом примере восстановления фазовой траектории по временной реализации синусоиды (t ) = sin t, рис.10.5,а.

При D = 1 (т.е. x(t ) = (t ) ) восстановленное множество лежит на отрезке прямой, рис.10.5,б. На нем точка с временным индексом k имеет соседей – точки l и s, соответствующие двум разным состояниям системы (разные по знаку производные (t ) ). В двухмерном пространстве ( x(t ) = [ (t ), (t + )] ) все точки «разойдутся», но точки k и l слабо, а k и s – сильно, рис.10.5,в. По этому признаку k и l называют «истинными», а k и s – «ложными» соседями.

Один из вариантов алгоритма таков. При пробной размерности D для каждого восстановленного вектора x k отыскивают одного (самого близкого) соседа;

увеличив D на 1, определяют, какие из соседей оказались ложными (сильно разошлись), а какие – истинными. Подсчитывают отношение числа ложных соседей к общему числу восстановленных векторов. Откладывают это число в зависимости от D, рис.10.5,г. Если при увеличении D это относительное число самопересечений уменьшается до нуля при некотором значении D *, то последнее и есть оценка размерности пространства, в котором достигается вложение траектории моделируемого движения. На практике, начиная с некоторой «правильной» размерности D *, число ложных соседей становится достаточно малым (но нестрого нулевым из-за шумов и т.д.) и далее не уменьшается. Такое значение D * принимают в качестве оценки размерности модели. В случае с синусоидой – это число 2, рис.10.5,г. Подробности см., например, в [115].

Глава 10. Реконструкция уравнений: «черный ящик»

Рис.10.5. Иллюстрация метода ложных соседей. а) временная реализация (t), значками помечены точки (k), (s), (l), соответствующие моментам времени k, s, l и близким значениям, а также точки, сдвинутые относительно них на = 3t. б) траектория, восстановленная в одномерном пространстве, в) траектория, восстановленная в двухмерном пространстве, г) график зависимости отношения числа ложных соседей к общему числу точек в восстановленной траектории от пробной размерности восстановленных векторов D 10.1.2.3. Метод главных компонент. Он может использоваться как для оценки размерности, так и для восстановления векторов состояния.

Метод используется в разных областях и имеет много названий, в приложении к реконструкции был предложен в [199]. Идея состоит в повороте осей координат в многомерном пространстве и выборе малого числа направлений, вдоль которых, в основном, развивается движение.



Pages:     | 1 |   ...   | 6 | 7 || 9 | 10 |
 





 
© 2013 www.libed.ru - «Бесплатная библиотека научно-практических конференций»

Материалы этого сайта размещены для ознакомления, все права принадлежат их авторам.
Если Вы не согласны с тем, что Ваш материал размещён на этом сайте, пожалуйста, напишите нам, мы в течении 1-2 рабочих дней удалим его.