markdown:
# «Б.П. Безручко, Д.А. Смирнов МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ И ХАОТИЧЕСКИЕ ВРЕМЕННЫЕ РЯДЫ Издательство ГосУНЦ «Колледж» Саратов, 2005 УДК ...»
Это оценка смещенная: она занижена, т.к. знаменатель больше числителя. На рис.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.
Траектория хаотической системы очень чувствительна к начальным условиям и параметрам.
Так, предлагалось делить исходный ряд на сегменты небольшой длины 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]. Идея состоит в повороте осей координат в многомерном пространстве и выборе малого числа направлений, вдоль которых, в основном, развивается движение.

Content-Type: text/html; charset=WINDOWS-1251


markdown:
Từ những hình ảnh dễ thương của chú chó Su, đến những phản ứng của những con khỉ khi nhìn thấy mình trong gương hay chú ngựa thổi nến... là những clip gây sốt trên mạng trong những ngày vừa qua.
*Hình ảnh của người đàn ông và chú chó nhỏ trong clip Puppyhood* *Clip Puppyhood - Nguồn: YouTube*
Mặc dù là một quảng cáo thương mại của hãng
** Purina** nhưng clip Puppyhood đã làm tan chảy nhiều người khi nhìn thấy những hình ảnh dễ thương của một người đàn ông đã chăm sóc và nuôi nấng một chú chó con khi mang chú về nhà.
Clip đã làm nhiều người xúc động khi thấy tình bạn đẹp đẽ giữa người đàn ông với chú cho nhỏ khi cả hai tìm thấy niềm vui khi "sống" cùng nhau.
Clip này đã thu hút đến hơn 2 triệu lượt người xem chỉ trong vòng 4 ngày kể từ khi đăng tải vào ngày 29-5.
*Clip phản ứng hài hước của những chú khỉ - Nguồn: YouTube*
Một clip hài hước khác về những phản ứng khác nhau của những con vật khi nhìn thấy chính mình trong gương. Người ta thử đặt một chiếc gương lớn trong rừng và theo dõi phản ứng của những con vật như thế nào khi chứng kiến chính mình trong gương.
Có những chú khỉ thì sợ hãi tấn công chính mình, có những con vật thì khép nép theo dõi, có con vật thì thử chạm vào mặt mình...
Clip này hiện nay đã có hơn 4,5 triệu lượt người xem khi chỉ trong vòng 2 ngày đăng tải trên YouTube.
*Clip chú ngựa thổi nến - Nguồn: YouTube*
Trong khi đó mọi người như "té ngửa" trong clip chú ngựa bất ngờ thổi nến. Chú ngựa Southside đã thổi nến và nhe răng cười khi mọi người cầm bánh sinh nhật trước mặt Southside và hát bài chúc mừng sinh nhật.
N.XUÂN
Theo Tuổi Trẻ

WARC-Type: response

Content-Type: application/http; msgtype=response

WARC-Identified-Payload-Type: text/html

Connection: close

Pragma: no-cache

P3P: CP="NOI ADM DEV PSAi COM NAV OUR OTRo STP IND DEM"

Cache-Control: no-store, no-cache, must-revalidate, post-check=0, pre-check=0

Set-Cookie: f43ef07555d0dd7ad70bd7202342ca87=dqb2vgp0co09ig4ecutc5b02i7; path=/; domain=www.marketingchienluoc.com; HttpOnly

Vary: Accept-Encoding,User-Agent

Content-Type: text/html; charset=utf-8


markdown:
**Đa dạng quà tặng**
Càng gần đến ngày tôn vinh phụ nữ, trên nhiều tuyến phố chuyên về hoa tươi như Kim Mã, Nguyễn Khuyến, Trần Xuân Soạn (Hà Nội) đã có dáng dấp của không khí chuẩn bị.
Theo tiết lộ của một số cửa hàng, giá hoa năm nay sẽ tăng nhẹ so với năm ngoái. Lý do là thời tiết năm nay diễn biến thất thường khiến nguồn cung có phần hạn chế. Hiện giá hoa trung bình khoảng 150.000-200.000 đồng/bó.
![ ](http://xmedia.nguoiduatin.vn/thumb_x500x/2013/10/10/image-thumb1381378957.jpg)
Không chỉ hoa, các mặt hàng dành cho phái nữ như nước hoa, mỹ phẩm, đồ trang sức, quần áo, giầy dép, đồ dùng nhà bếp cũng được bán khá nhiều, với các mức giá vô dùng hấp dẫn. Khác với hoa tươi, những mặt hàng này mang tính thiết thực hơn, được nhiều nam giới và các chị em phụ nữ quan tâm.
Chị Trần Thị Thu Phương - chủ cửa hàng
[quà lưu niệm](/timkiem-qu%C3%A0%20l%C6%B0u%20ni%E1%BB%87m.html) trên đường Tôn Thất Tùng cho biết ngay từ đầu tháng, cửa hàng đã tăng 50% lượng hàng dự trữ.
Theo chị Phương, năm nay có rất nhiều sản phẩm để cho khách hàng lựa chọn như vật dụng trang trí, mỹ phẩm, hộp đựng trang sức, kẹp tóc, gấu bông, túi xách. Hiện sức mua bắt đầu tăng khoảng 20% so với tuần trước, cao điểm mua sắm sẽ rơi vào những ngày cận kề.
Tuy nhiên, do ảnh hưởng của kinh tế nên nhiều cửa hàng dự báo sức mua dịp lễ 20/10 năm nay không tăng đột biến, chủ yếu nhóm quà tặng có giá trị trung bình, từ 100.000-300.000 đồng.
Bên cạnh đó, du lịch cũng là sự lựa chọn của nhiều khách hàng. Năm nay, 20/10 rơi vào cuối tuần, nên kỳ nghỉ lễ ngắn ngày được du khách lựa chọn nhiều nhất, trong đó có các tour trong nước, đặc biệt là tuyến miền Trung.
Ngoài ra, các spa, trung tâm thẩm mỹ còn tung ra nhiều gói thư giãn, phiếu quà tặng để khách mua làm quà cho chị em nhân dịp này.
**Tưng bừng khuyến mãi**
Đến hẹn lại lên, 20/10 cũng là dịp để các doanh nghiệp đưa ra nhiều hình thức khuyến mãi thu hút khách hàng. Tại các siêu thị, quà tặng 20/10 chủ yếu là các hộp quà gói sẵn các sản phẩm dầu gội, sữa tắm, mỹ phẩm có giá từ 55.000-165.000 đồng/hộp. Chuỗi hệ thống siêu thị Co.opmart tổ chức chương trình Thân thương dáng Việt với nhiều chương trình khuyến mãi lớn giảm giá đồ dùng và may mặc lên tới 50%.
Các cửa hàng quần áo, trung tâm thẩm mỹ, yoga cũng có nhiểu hình thức tặng quà hấp dẫn.
Bên cạnh đó, các công ty du lịch cũng giảm giá nhiều tour nội địa, nước ngoài để thu hút khách mua tour đi chơi dịp lễ.
Cùng với tour trong nước, những chuyến du lịch ngắm Thu nước ngoài cũng được các hãng lữ hành giới thiệu nhân dịp khuyến mãi mừng 20/10.
Mặc dù vậy, một số hãng lữ hành cho biết, 20/10 năm nay rơi vào ngày cuối tuần nhưng đến thời điểm này, khách mua tour khởi hành dịp này hiện chưa nhiều.
Theo
* Thông tấn xã Việt Nam*
|
|
Mời bạn đọc cùng làm báo với Nguoiduatin.vn, viết tin bài, bày tỏ quan điểm về mọi mặt của cuộc sống gửi về hộp thư

WARC-Truncated: length

WARC-Identified-Payload-Type: application/xhtml+xml

Vary: Accept-Encoding

Content-Type: text/html


markdown:
Thêm bài hát này vào danh sách Playlist
Lời đăng bởi: siiiidoooo
Canve chiều buồn, một trời hoang xơ xác,
Canve chiều buồn, nịnh giận bóng nhân gian,
Canve chiều buồn, Một Mình Con Chúa trời,
gánh lấy tội đời, để giải thoái con người.
ĐK : Canve xưa ấy, một mình GIÊSU
hy sinh gánh lấy, tội đời ô nhơ, tình yêu thương đó, bao la vô bờ,
đã chết vì yêu , trần thế hững hờ.
con xin hối lỗi, tội đời bao năm ,
xin thương tha thứ, từng ngày đi hoang,
giờ đây xám hối, ăn năn 1 lòng,
nguyện Chúa rủ thương, đừng chấp tội con.
2. Bao nhiêu dập vùi, hình hài ôi đẫm máu,
hy sinh chuộc tội, chịu ngàn nỗi thương đau,
yêu thương loài người, để rồi thân rã rời,
ai hay tình Ngài! đã chịu chết đớn đau. :( ..
Con xin hối lỗi, tội đời bao năm ,
xin thương tha thứ, từng ngày đi hoang,
giờ đây xám hối, ăn năn 1 lòng, nguyện Chúa rủ thương,
đừng chấp tội con.
3. Thân con mỏi mòn, chẳng còn đâu hơi sức,
cho con nghị lực, để chịu những thương đau,
xin tha tội đời ,vi ngày bao lỗi lầm ,
xin vâng lời Ngài, trọn lời hứa hôm nao.
4. Cô đơn tủi sầu, một mình Con gánh vác,
tâm can nặng buồn, Vì đời đã vong an,
lê thân miệt mài, mở đường ơn cứu độ,
cho nhân loại được, hưởng nguồn sống vô bờ. ..ĐK( 2 lần) \" God bless to you\"
Gia Ân
V.A
KaTa
Cà Dễ Dãi
Hồ Quang Hiếu
Hilary Duff
xem thêm
rút gọn
Top 100 bài hát
Top 100 Nhạc Hàn Hay Nhất
Top 100 Nhạc Country Hay Nhất
Top 100 Nhạc Hoa Hay Nhất
Top 10 Blues/Jazz Hay Nhất
Top 100 Nhạc Nhật Hay Nhất
Đơn vị chủ quản: Công ty Cổ phần NCT
Địa chỉ: Tòa nhà HAGL Safomec, 7/1 Thành Thái, P14, Q10, TP.HCM
Người chịu trách nhiệm nội dung: Ông Nhan Thế Luân - Email: support@nct.vn - Tel: (08) 3868 7979
Giấy phép MXH số 499/GP-BTTTT do Bộ Thông Tin và Truyền thông cấp ngày 28/09/2015
## Bình luận

Server: nginx+

P3P: CP='NOI ADM DEV PSAi COM NAV OUR OTRo STP IND DEM'

Vary: Accept-Encoding, User-Agent

Content-Type: text/html; charset=UTF-8


markdown:
Dựa trên những chia sẻ của HopperH, một công cụ đăng tải hình ảnh lên Instagram theo lịch trình sẵn có, nhiều người sẽ biết cách để trang Instagram của mình trở nên thu hút hơn.
Công cụ lập kế hoạch đăng ảnh lên Instagram HopperH đã đưa ra một bộ dữ liệu phân tích đặc biệt về Instagram. Trong báo cáo chi tiết của HopperH, công cụ này đã đưa ra nhiều gợi ý về thời gian đăng Instagram tốt nhất, bộ lọc, caption hay nhiều thông tin khác tốt nhất để bạn có thể tạo được một bài đăng gây ấn tượng nhất với mọi người.
Theo đó, bộ lọc được ưa chuộng và nhận được nhiều like, lượt thích nhất là Gingham. Tiếp sau đó là Claredon.
Ngoài ra, các bài đăng từ người dùng nữ cũng nhận được nhiều lượt tương tác và tham gia của mọi người hơn nam giới. Trung bình các bài đăng của phái nữ nhận được 578 lượt like so với 117 like của phái nam.
Thời gian đăng Instagram tốt nhất là vào buổi tối, từ 6h tới 7h30 (áp dụng trong tuần). Nếu bạn đăng ảnh vào cuối tuần, hãy lựa chọn khung giờ tốt nhất vào 11h sáng.
Các lưu ý quan trọng khác như các bức ảnh ghi lại trải nghiệm cá nhân như du lịch, hành trình giảm cân sẽ nhận được nhiều quan tâm hơn.
Cách lựa chọn hashtag thông minh cũng là cách giúp các bức ảnh của bạn thêm nổi bật và phổ biến hơn. Hashtag phổ biến nhất hiện nay là #love với hơn 1 tỷ bài đăng sử dụng hashtag này. Đứng vị trí cuối là #me với khoảng 312,8 triệu bài đăng.
Một bài đăng Instagram luôn đi kèm một dòng caption (tiêu đề). Các phân tích đã chỉ ra, độ dài caption quyết định khá nhiều tới việc người dùng có bấm like hay không. Nếu họ phải nhấn vào chữ More (Thêm) để xem caption, họ sẽ ít có xu hướng nhấn Like cho bài đăng của bạn. Ngoài ra, một caption có chứa một emoji nào đó cũng rất thú vị và thu hút hơn.

Cache-Control: private

Set-Cookie: websitecloseoropen=True; expires=Mon, 29-May-2017 21:23:00 GMT; path=/


markdown:
Nhằm đáp ứng nhu cầu khách hàng ! ONB.VN triển khai dịch vụ quảng cáo Google với giá tốt nhất cho khách hàng ! ONB.VN luôn cam kết các dịch vụ chất lượng, đúng giá thị trường ! Hiện tại nhiều nơi với lời chào quảng cáo google giá rẻ không giới hạn số lượng click, hiện thị hoàn toàn là lừa gạt, không có thật quý khách nên cẩn thận lưu ý. Vì google adword chỉ tính tiền trên lượt click và ngay cả Google Adword trang chính không hề cung cấp 1 gói nào là không giới hạn lượt click cả.
Một ví dụ điển hình 1 số khách hàng làm web tại onb ( trước khi ONB.VN làm dịch vụ quảng cáo google ) thì chi phí họ đóng tầm 800k - 1 triệu 1 tháng cho gói không giới hạn ! mà quảng cáo so với trước khi quảng cáo, số lượt người có tầm 10 người nhấp 1 ngày ! Mà ngay cả thông tin thống kê thì bên quảng cáo cũng giấu ko giao cho khách hàng biết ! nên khách hàng càng đóng phí cao thì tiền càng chảy vào túi họ nhiều hơn.
Khi hợp tác với ONB.VN, quý khác sẽ tư vấn rõ hơn về dịch vụ, từ khoá đó mức độ cạnh tranh như thế nào, dự trù 1 click mất bao nhiệu, tổng chi phí 1 tháng tầm bao nhiêu ...... và ONB.VN cung cấp thống kê mỗi ngày cho khách hàng biết 1 ngày có bao nhiêu lượt nhấp và chi phí là bao nhiêu !
"Đây chính là ưu điểm của ONB.VN ! Cam kết làm đúng sự thật không lừa gạt !"
Quý khách hàng lưu ý:
Nếu quý khách muốn quảng cáo với chi phí lớn hơn, vui lòng liên hệ với chúng tôi qua email trancongnghiep@yahoo.com để được phục vụ tốt nhất.
Chú thích:

Set-Cookie: websitecloseoropen=True; expires=Mon, 29-May-2017 21:27:38 GMT; path=/


markdown:
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
*Hãy quyết định ngay và liên hệ với chúng tôi để có báo giá tốt nhất:* **CÔNG TY CP TM DV AN SƯƠNG**
Website :
[www.otoansuong.vn](http://www.otoansuong.vn/)
Gmail : anhototai@gmail.com
**Hotline : 0933 770 688 - 0908.541.699** *" Rất hân hạnh được phục vụ quý khách "*
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
*Hãy quyết định ngay và liên hệ với chúng tôi để có báo giá tốt nhất:* **CÔNG TY CP TM DV AN SƯƠNG**
Website :
[www.otoansuong.vn](http://www.otoansuong.vn/)
Gmail : anhototai@gmail.com
**Hotline : 0933 770 688 - 0908.541.699** *" Rất hân hạnh được phục vụ quý khách "*

Accept-Ranges: bytes

Cache-Control: private, no-cache, no-store, proxy-revalidate, no-transform

Set-Cookie: currency=USD; expires=Wed, 21-Jun-2017 21:14:15 GMT; path=/; domain=www.otoansuong.vn


markdown:
LONDON, November 21, 2012 /PRNewswire/ --
A special purpose vehicle has been formed by Mauricio Trevino Zambrano to bid for the entire share capital of Nordic Petroleum ASA.
After 2 years of protracted negotiations between Mercom Oil Sands and Nordic Petroleum ASA with disagreements and un-workable contracts Mercom Oil Sands PLC seems to have reached a crossroads.
NAC has been formed to launch a bid for the entire share capital of Nordic Petroleum and is actively seeking willing shareholders to join in creating an exciting opportunity with a team that can raise capital, create value and effectively manage the assets.
Mr. Zambrano was instrumental in bringing the single largest contribution to the IPO of Mercom Oil Sands PLC and has a huge desire to see the transaction reap rewards and be successful.
Should an offer be successful NAC intends to offer the CHARD asset to Mercom but no discussion at board level or agreements to this effect are in place at this time.
An offer has been submitted to the board of NOP for good measure but NAC now seeks to deal only with NOP shareholders directly.
SOURCE Nordic Acquisition Co Ltd

Server: Apache


markdown:
[0](/mua-hang/gio-hang)
# LiPo 1S 400mAh Novus Sea Cobra
Sản phẩm này được bán tại
[www.amazon.com](http://www.amazon.com/LiPo-400mAh-Novus-Sea-Cobra/dp/B0032U2RHE)
Thương hiệu:
[Great Planes](/thuong-hieu/great-planes_th6807)
Ngay sau khi quý khách đặt hàng và thanh toán, ShipTo.vn sẽ mua sản phẩm từ Mỹ và vận chuyển về Việt Nam
Sản phẩm tạm thời ngưng bán trên Amazon. ShipTo.vn sẽ theo dõi và tìm mua từ nhà cung cấp khác
Chọn số lượng:
[Xem chi tiết sản phẩm tại www.amazon.com](http://www.amazon.com/LiPo-400mAh-Novus-Sea-Cobra/dp/B0032U2RHE)
-
Giá giảm 10-50%
so với hàng xách tay
-
100% sản phẩm thật
được mua trực tiếp từ Mỹ
-
Hoàn tiền 100%
nếu hàng bị thất lạc hư hỏng
-
Miễn phí dịch vụ
nếu giao hàng trễ
#### Mua hàng theo yêu cầu
-
1.Đặt mua
Đặt mua những sản phẩm được ShipTo.vn tuyển chọn, hoặc gửi link của sản phẩm muốn mua để nhận báo giá.
-
2.Thanh toán
Thanh toán cho ShipTo.vn (tiền mặt, chuyển khoản, hoặc thẻ tín dụng).
-
3.Vận chuyển
ShipTo.vn mua hàng từ nước ngoài, chuyển về Việt Nam, và lo tất cả thủ tục hải quan.
-
4. Giao hàng
ShipTo.vn giao hàng tận nơi miễn phí trong vòng 10 - 20 ngày làm việc kể từ ngày thanh toán.
#### Vận chuyển theo yêu cầu
-
1.Đăng ký địa chỉ Mỹ
Đăng ký tài khoản ShipTo.vn có ngay địa chỉ riêng tại Mỹ.
-
2.Tự mua hàng
Tự mua tại bất kỳ website nào và chuyển hàng tới địa chỉ Mỹ của bạn.
-
3.Vận chuyển
ShipTo.vn nhận hàng, kiểm hàng và chuyển quốc tế.
-
4. Giao hàng
ShipTo.vn thông quan, giao hàng tận tay và thu tiền dịch vụ.
ĐĂNG KÝ EMAIL ĐỂ BẮT ĐẦU

Cache-Control: private, no-cache="Set-Cookie"

Set-Cookie: website_id=1; expires=Tue, 22-May-2018 17:00:00 GMT; path=/


markdown:
[Lâm Đồng: Cung đường đẹp hiếm có mê hoặc giới trẻ khắp Việt Nam](http://www.tintaynguyen.com/lam-dong-cung-duong-dep-hiem-co-hoac-gioi-tre-khap-viet-nam/214326/)
Tà Năng là một xã thuộc huyện Đức Trọng, tỉnh Lâm Đồng. Nơi này luôn là cái tên được các phượt thủ khao khát chinh phục, bởi đây được mệnh danh là...

Server: cloudflare-nginx

X-Content-Type-Options: nosniff

X-Frame-Options: SAMEORIGIN

Set-Cookie: __cfduid=df15264614e0e56e03c081e9425e85d101495489253; expires=Tue, 22-May-18 21:40:53 GMT; path=/; domain=.tintaynguyen.com; HttpOnly

X-XSS-Protection: 1; mode=block

Link: <https://goo.gl/dYUjAe>; rel=shortlink


markdown:
|
|
##
|
|![]()
|
|![]()
|
|
![]()
![]()
|
|
|
|
|
|
|
||
|
|
|![]()
|
|
## Trò chơi Xếp hình kiểu 375 - Mã số trò chơi: 13313
## Game thuộc chủ đề Phá khối hình
|
|
## Hướng dẫn chơi game
|
|![]()
|Tham gia những trò chơi vui nhộn với kho game lớn nhất hiện nay với trên 20.520 trò chơi hấp dẫn dành cho máy tính và điện thoại có hỗ trợ công nghệ Flash, chơi game miễn phí và cực nhanh trên TroChoiVui.com. Những ngày lễ hội đặc biệt trong năm, chúng tôi đều cập nhật nhiều trò chơi hấp dẫn cho lễ hội Valentine, Halloween, Noel,.. Tham gia với các nhân vật hoạt hình được nhiều người yêu thích như Doraemon, Mario, Sonic, Pikachu, Dora, Ben10, Scooby Doo, Simpson, Pou,.. Mỗi ngày, chúng tôi luôn cập nhật số lượng lớn các trò chơi đa dạng ở nhiều thể loại như game đánh nhau, đấu võ, dàn trận, kinh doanh, đua xe, bạn gái, văn phòng, trẻ em, tô màu, phiêu lưu, bắn súng, nấu ăn, thời trang, xếp hình, xếp gạch, trang điểm, xếp bi, cá lớn nuốt cá bé, bắn xe tăng, bắn máy bay, bắn súng, bắn cung, bắn ná, đấu kiếm, ném đồ vật, rắn săn mồi, mặt trăng ăn vàng, hứng đồ vật, thiết kế người mẫu, sơn móng tay chân, câu cá, trượt tuyết, lướt ván, khúc côn cầu, bóng đá, bóng chuyền, bóng rổ, bóng bàn, cầu lông, tennis, cờ vua, cờ tướng, domino,.. Chắc chắn bạn sẽ có những giây phút thú vị khi chơi game tại trochoivui.com hàng ngày.|
|
Kết nối với chúng tôi qua |
Bản quyền thuộc về
Giấy phép ICP số 181/GP-CBC. Cơ quan chủ quản: Công ty TNHH Không Gian Việt - Website:


markdown:
|
|
##
|
|![]()
|
|![]()
|
|
![]()
![]()
|
|
|
|
|
|
|
||
|
|
|![]()
|
|
## Trò chơi Nhanh tay ghép khối - Mã số trò chơi: 8642
## Game thuộc chủ đề Khả năng sắp xếp
|
|
## Hướng dẫn chơi game
|
|![]()
|Tham gia những trò chơi vui nhộn với kho game lớn nhất hiện nay với trên 20.520 trò chơi hấp dẫn dành cho máy tính và điện thoại có hỗ trợ công nghệ Flash, chơi game miễn phí và cực nhanh trên TroChoiVui.com. Những ngày lễ hội đặc biệt trong năm, chúng tôi đều cập nhật nhiều trò chơi hấp dẫn cho lễ hội Valentine, Halloween, Noel,.. Tham gia với các nhân vật hoạt hình được nhiều người yêu thích như Doraemon, Mario, Sonic, Pikachu, Dora, Ben10, Scooby Doo, Simpson, Pou,.. Mỗi ngày, chúng tôi luôn cập nhật số lượng lớn các trò chơi đa dạng ở nhiều thể loại như game đánh nhau, đấu võ, dàn trận, kinh doanh, đua xe, bạn gái, văn phòng, trẻ em, tô màu, phiêu lưu, bắn súng, nấu ăn, thời trang, xếp hình, xếp gạch, trang điểm, xếp bi, cá lớn nuốt cá bé, bắn xe tăng, bắn máy bay, bắn súng, bắn cung, bắn ná, đấu kiếm, ném đồ vật, rắn săn mồi, mặt trăng ăn vàng, hứng đồ vật, thiết kế người mẫu, sơn móng tay chân, câu cá, trượt tuyết, lướt ván, khúc côn cầu, bóng đá, bóng chuyền, bóng rổ, bóng bàn, cầu lông, tennis, cờ vua, cờ tướng, domino,.. Chắc chắn bạn sẽ có những giây phút thú vị khi chơi game tại trochoivui.com hàng ngày.|
|
Kết nối với chúng tôi qua |
Bản quyền thuộc về
Giấy phép ICP số 181/GP-CBC. Cơ quan chủ quản: Công ty TNHH Không Gian Việt - Website:


markdown:
# Tôi phải làm gì đây khi em gái chồng yêu mình?
**( Theo NLĐO) - Tôi như rơi vào hố sâu khi em gái của chồng mình thổ lộ : "Trái tim em đóng băng, không thể run lên với bất cứ người đàn ông nào? Nó chỉ có thể bồi hồi trong lồng ngực và neo lại ở mỗi một người, là ai không lẽ chị không hề cảm nhận được sao?"…**
Tôi và Hoa là bạn học với nhau từ thời cấp hai. Hoa là một người con gái đẹp nhưng rất lạnh lùng, kiêu kỳ với tất cả mọi người. Nhưng với riêng tôi thì lại hoàn toàn khác.
Vì rất thân với nhau, mỗi lần càu nhàu hay giận ai, Hoa lấy xe chạy vèo đến gặp tôi, xả hết bực dọc trong lòng cho tôi nghe. Bù lại, những lúc rảnh, tôi cũng thường sang nhà Hoa để "tám" đến hết những chuyện trên trời dưới đất mới thôi. Điều lạ là tôi chẳng thấy chàng trai nào bên cạnh Hoa. Có lẽ là do Hoa rất lạnh lùng, thờ ơ với con trai.
Khác hẳn với em gái của mình, anh trai Hoa lại rất dịu dàng và thân thiện. Ngay lần đầu tiên gặp lúc anh vừa mới đi du học nước ngoài trở về, tôi đã nhận ra mình thích anh. Rồi tôi và anh gặp nhau thường xuyên hơn.
Biết chuyện, Hoa tìm cách tránh mặt tôi. Ở Hoa có một nỗi buồn gì đó hoặc một tâm sự gì đó được giấu kỹ trong lòng không cho ai biết. Một hôm Hoa hẹn gặp riêng tôi, nghiêm nghị nói: "Bồ hãy suy nghĩ cho kỹ trước khi đến với anh mình nhé". Giọng điệu và ánh mắt của Hoa khiến tôi có cảm giác hình như Hoa không muốn tôi làm chị dâu.
Nhưng rồi tôi và anh Hoa cũng đến với nhau. Sau đám cưới, tôi về sống chung một nhà với Hoa. Những ngày đầu, Hoa vẫn cố tình né tránh tôi nhưng dần dà sau đó, Hoa vui vẻ hơn, bớt u buồn, thậm chí hay tìm cách gần gũi với tôi hơn. Nhưng Hoa vẫn không lấy chồng.
Từ lúc về ở chung tôi chưa bao giờ thấy Hoa dẫn người đàn ông nào về nhà, dù chỉ là bạn bình thường. Trái lại, sau ngày làm việc, Hoa dành hết thời gian cho gia đình, chính xác là cho tôi. Thường khi đi làm về, Hoa ôm chầm lấy tôi rồi hôn một cái chẳng ngại ngùng gì. Có khi còn bình thản ngồi vào lòng tôi ngay tại phòng khách xem ti vi. Nói thật, thân và sống với nhau bao nhiêu năm rồi nhưng tôi vẫn có phần dị ứng với điều này. Tôi có nhắc nhở về sự đụng chạm thì Hoa lại cười hề hề: "Tại em quý chị nên mới thế!". Có lần khác vào phòng tôi hứng chuyện, Hoa còn vô tình đòi ngủ chung với chị dâu, báo hại anh trai phải đi tìm chỗ khác mà ngủ…
Sau bao nhiêu chuyện như thế, mọi người nghiêm khắc góp ý, Hoa có giữ chừng mực. Nhưng chỉ được một thời gian ngắn, đâu lại vào đấy. Bi kịch thật sự chỉ bắt đầu từ một đêm đông rét mướt, chồng tôi đi công tác, còn Hoa trở về nhà trong trạng thái say xỉn và vô tình đột nhập vào phòng tôi. Bị mất giấc ngủ, bực mình tôi hét lên : "Sao suốt ngày em chỉ biết lo chơi bời thế, không biết rằng ba mẹ đang sốt ruột về chuyện chồng con của em sao?". Hoa lè nhè đáp : "Em có muốn như vậy đâu chứ. Nhưng tại sao trái tim em đóng băng, không thể run lên với bất cứ người đàn ông nào? Nó chỉ có thể bồi hồi trong lồng ngực và neo lại ở mỗi một người, là ai không lẽ chị không hề cảm nhận được sao? Em biết là mình đã sai khi cứ mải miết chạy theo ảo ảnh mà mãi mãi không thể nào thuộc về mình…".
Tôi như đang rơi xuống một cái hố sâu khi em chồng nói yêu mình. Tôi cảm nhận rõ ràng hơn bao giờ hết sự chơi vơi, chênh vênh về sự luyến ái của Hoa và hiểu ra vì sao bao năm nay Hoa vẫn không chịu lấy chồng.
Tôi trở nên hoảng sợ. Trong suy nghĩ của tôi, chúng tôi mãi mãi vẫn là chị dâu em chồng, là hai đường thẳng song song, chỉ ở cạnh nhau, nương tựa vào nhau nhưng không thể là giao điểm được.
Tôi phải làm gì đây? Có nên tiếp tục sống chung trong một mái nhà?…
Bạn thích truyện này?

X-Cache: MISS