Ю. П. Драница Судоводительский факультет мгту, кафедра высшей математики и программного обеспечения ЭВМ

Вид материала

Документы

Содержание 2.2. Основные проблемы моделирования 2.3. Современные методы оценки тренда 3. Предлагаемый метод решения задачи М., Недра

Подобный материал:

• Н. И. Лобачевского Факультет Вычислительной математики и кибернетики Кафедра Математического, 169.45kb.
• Н. И. Лобачевского Факультет Вычислительной математики и кибернетики Кафедра Математического, 172.6kb.
• Н. И. Лобачевского Факультет Вычислительной математики и кибернетики Кафедра Математического, 123.69kb.
• Н. И. Лобачевского Факультет Вычислительной математики и кибернетики Кафедра Математического, 132.68kb.
• Н. Э. Баумана Факультет "Инженерный бизнес и менеджмент" Кафедра "Менеджмент", 786.11kb.
• И. И. Мечникова Институт математики, экономики и механики Кафедра математического обеспечения, 900.66kb.
• По ЭВМ перечень примерных контрольных вопросов и заданий для текущей работы, 40.23kb.
• Кафедра системного анализа, 34.57kb.
• Программа дисциплины «компьютерная подготовка», 147.49kb.
• Правительство Российской Федерации Национальный исследовательский университет Высшая, 320.55kb.

Об одном методе моделирования нестационарных динамических систем и процессов


Ю.П. Драница

Судоводительский факультет МГТУ, кафедра высшей математики
и программного обеспечения ЭВМ


Аннотация. Разработан корректный метод оценки нестационарности многомерных динамических процессов и ее учета при построении линейных математических моделей. Теория метода основана на аппроксимации процесса системой линейных разностных уравнений в прямом и обратном времени, при этом не требуется знания ковариационной матрицы. Метод является обобщением теории Дж. Бурга, разработанной для одномерного и стационарного случая. Выполнена оптимизация алгоритма за счет использования свойств матриц Теплица и рекурсии Левинсона.


Abstract. The correct method of an estimation of nonstationarity of the multidimensional dynamic processes and its taking into account when constructing linear mathematical models has been worked out. The theory of the method is based on approximation of the process by the linear differential equations system in direct and inverse time, and the knowledge of a covariance matrix is not required. The method is the generalization of the J. Burg theory, worked out for the one-dimensional and stationary case. The optimization of the algorithm has been executed by using properties of Teplitz matrixes and Levinson recursion.


1. Введение

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

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

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


2. Современное состояние проблемы

2.1. Исходные положения

Пусть наблюдения за некоторой динамической системой или процессом проводятся в дискретные моменты времени t_i, i = 1,2,...,N, где N - общее число наблюдений. Измерения за состоянием системы характеризуется вектором S(t_i) c размерностью m. Следовательно, совокупность наблюдений образует матрицу S размерностью m´N. Отметим, что вместо независимой переменной типа времени могут использоваться другие упорядоченные параметры, например, глубины, высоты и т.п., однако без ограничения общности будем считать, что независимой переменной является время.

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


S(t_i) = U(t_i) + Е(t_i), (1)


где U(t_i) - вектор истинных значений параметров состояния системы в момент времени t_i; Е(t_i) - суммарный вектор ошибок измерений в тот же момент времени; размерность вектора Е(t_i) равна m.

Выражение (1) представляет собой линейную модель процесса наблюдения и указывает на независимость сигнальной части измерений U(t_i) и помехи Е(t_i). Сигнальная часть измерений, в зависимости от физической постановки задачи, может считаться либо неизвестным детерминированным вектором, либо набором случайных величин. Помехи, с одной стороны, представляют ошибки измерений, а с другой - отражают неполноту и ограниченность модели данных (1). В общем случае под помехами будем понимать любую информацию, которая не является полезной для отыскания вектора истинных значений параметров состояния системы в рамках принятой модели. Иногда помехи называют "шумом" наблюдений. Неполнота и ограниченность модели (1) может, в частности, заключаться в том, что реальные данные S(t_i) содержат также систематические компоненты, или тренды параметров. В этом случае более адекватной моделью наблюдений будет следующая


S(t_i) = A(t_i) + V(t_i) + Е(t_i), (2)


где A(t_i), V(t_i) - соответственно вектора трендовой и "сигнальной" составляющих данных размерности m; E(t_i) - по-прежнему вектор шума.

Относительно суммарного вектора ошибок Е(t_i) принято считать, что он образует квазислучайную (коррелированную) совокупность данных. Трендовая составляющая отражает временную нестационарность процесса измерений и может иметь достаточно сложную природу. Например, она может состоять из систематического ухода показаний измерительных приборов, связанных с техническими особенностями этих устройств (приборные ошибки), и смещений, которые отражают реальный тренд наблюдаемого объекта. В этом случае вместо (2) должна быть использована более сложная модель, учитывающая данную реальность. Будем, однако, считать, что все приборные ошибки не имеют тренда и сосредоточены в "шумовой" компоненте сигнала. Это не умаляет общности рассуждений, так как при необходимости всегда может быть применена более адекватная модель наблюдений.


2.2. Основные проблемы моделирования

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

Допустим, что по наблюдениям типа (1) или (2) требуется построить математическую модель динамической системы. Наиболее простой путь решения этой задачи заключается в использовании линейного метода. Математической основой такого подхода являются хорошо развитая теория линейных стационарных операторов и преобразование Фурье, которые позволяют ввести понятия передаточной функции, спектра, эргодичности и пр. Известно, что исчерпывающей информацией о линейной динамической системе является матрица ковариаций вектора состояний, позволяющая синтезировать ее математическую модель, которая ясно отражает такие ее физические свойства, как распределение энергии по частотам, устойчивость, наличие периодичностей, характер их затухания или возбуждения и другие параметры. Следовательно, основной задачей при линейном моделировании является корректная оценка ковариационной матрицы вектора состояний, в том числе его трендовой, "сигнальной" и "шумовой" компонент.

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

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

Рассмотрим, к чему может приводить неучет тренда.

Эта задача тесно связана с проблемой влияния скрытых систематических ошибок на параметры идентифицируемой модели. Пусть полная модель некоторой совокупности данных описывается, например, выражением (2). Будем считать, что часть этой модели нам неизвестна и представляет собой скрытую систематическую ошибку. Это означает, что вместо соотношения (2) мы используем неполную модель (1), где U(t_k) = A(t_k) + V(t_k). В работе (Макузе, 1989) показано, что скрытые ошибки приводят к систематическому смещению оценок параметров модели относительно истинных значений. Практические следствия этого могут быть весьма существенными, способными создать совершенно неадекватные представления о свойствах и параметрах изучаемого явления. Таким образом, задача объективного выявления нестационарности и ее учета является актуальной как с точки зрения разработки моделей с переменными коэффициентами, так и с постоянными для выработки критерия адекватности модели натурным данным.

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

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


2.3. Современные методы оценки тренда

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

Основная возникающая при этом проблема заключается в том, что для того, чтобы сделать эти оценки с достаточной для практических нужд точностью, требуются большие объемы информации. Так, например, в (Котюк и др., 1967) показано, что в случае разложения функции автокорреляции (ФАК) по системе затухающих косинусоид ошибки ее оценивания представляются суммой ошибок по каждой гармонике. При этом надежное оценивание ординат ФАК возможно лишь для тех составляющих процесса, период которых в 5-10 раз меньше продолжительности наблюдений. Кроме того, ошибка оценивания возрастает с увеличением временного запаздывания между функциями, следовательно, достоверные оценки ФАК могут быть получены только при небольших сдвигах. Это означает, что по ограниченным по длине данным нельзя с достаточной для практических целей точностью оценить не только тренд процесса, но и его наиболее низкочастотные компоненты. В то же время синтез модели только по достоверным оценкам ковариационной функции неизбежно приводит к ее недостаточному разрешению в области низких частот, что снижает ценность получаемых решений. Однако построение линейной динамической модели возможно и без знания ковариационной функции. Рассмотрим эту задачу более подробно.


3. Предлагаемый метод решения задачи

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


X(t_k) = S_i A_bi(t_k) X(t_k-i) + f(t_k), i = 1,2,…,l , (3)

1. 
   
2. где X(t_k) - m-мерный вектор состояния системы (фазовые переменные); A_bi(t_к) - заданные матрицы размерности m´m (переходные матрицы); f(t_k) - неизвестный m-мерный возмущающий вектор.
   

В общем случае фазовый вектор X определяется как множество таких переменных, которые вместе с известными входными сигналами необходимы и достаточны для нахождения отклика системы. Матрицы A являются матрицами частных производных вида ¶X(t_k)/¶X(t_k-i) i = 1,2, ... ,l, которые, возможно, зависят от времени. Возмущающий вектор f(t_к) представляет, с одной стороны, воздействие на систему внешней среды, а с другой - ошибки, связанные с неадекватностью модели натуре.

Уравнение (3) называется дискретной моделью состояния системы. Это линейное разностное уравнение порядка m´l. С другой стороны, эта модель является уравнением многомерной линейной авторегрессии, или одношаговым линейным предсказателем. Это уравнение позволяет предсказать значение фазового вектора X(t_k) в момент времени t_k по его значениям в прошлые моменты времени t_k-1, t_k-2, ... , t_k-l. Число участвующих в прогнозе прошлых значений фазового вектора называется лагом модели. Модель (3) представляет так называемый физически реализуемый линейный многомерный цифровой фильтр. Его мерность определяется числом компонент фазового вектора.

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


X(t_k) = S_i A_pi(t_k) X(t_k+i) + f(t_k), i = 1,2,…,l. (4)

3. 
   

В частном случае, когда матрицы A_bi, A_pi не меняются во времени, модели (3-4) описывают стационарную во времени систему. При одномерном фазовом векторе (m = 1), модели (3-4) упрощаются. Векторы X(t_k) превращаются в скалярные величины X(t_k), вектора f(t_k) - в случайные функции времени f(t_k), а совокупность матриц A - в вектор строки. Модель (3) часто называют перспективным предсказанием, а (4) - ретроспективным.

Задача идентификации моделей (3-4) заключается в определении матриц A_i(t_k) и лага l по априорным или экспериментальным данным. При известной матрице ковариаций вектора X, решение этой задачи хорошо известно (Клаербоут, 1981). Если фазовый вектор непосредственно измеряем, то его матрица автоковариаций может быть оценена и по экспериментальным данным. Но, как показано выше, существенным ограничением такого подхода является то, что достигнутая при этом точность является достаточно высокой только при небольших лагах модели. Однако практически лаг модели заранее, как правило, неизвестен, и должен быть оценен по экспериментальным данным. Все эти проблемы сильно ограничивают возможности определения параметров моделей (3-4) по методу непосредственной эмпирической оценки ковариационной матрицы.

Для одномерного случая и в стационарном приближении синтез моделей (3-4) возможен и при неизвестной автоковариационной матрице (Сильвия, Робинсон, 1983; Клаербоут, 1981). Основная идея этого метода была разработана Дж. Бургом для расчета спектров высокого разрешения и заключалась в симметричном использовании как перспективных, так и ретроспективных предсказаний. Важным достоинством этого подхода, по сравнению с классическими схемами, является возможность построения моделей с лагами, превышающими область достоверных оценок ФАК, рассчитанных по общепринятым методам. Последнее достигается благодаря статистически более достоверной экстраполяции параметров модели. Преимущества этого метода особенно полезны при работе с короткими выборками, на которых общепринятые статистические методы часто оказываются неприемлемыми геофизически. Эта ситуация хорошо была изучена Дж. Бургом, в результате чего он использовал некоторые специфические свойства матриц Теплица для разработки своей процедуры максимально-энтропийного оценивания спектра и параметров одномерной линейной и стационарной модели. Вычислительная эффективность метода основана на рекурсии Левинсона, которая позволяет строить очень экономичные алгоритмы. Отметим, что рекурсия Левинсона для временной области имеет такое же важное значение, как и быстрое преобразование Фурье для частотной. Полученная в ходе вычислений информация может быть использована затем для более корректной эмпирической оценки корреляционной матрицы фазового вектора.

Распространим идею Бурга на синтез моделей (3)-(4) по экспериментальным данным на многомерный случай. При этом предполагается, что ковариационная функция неизвестна, а моделируемый процесс нестационарен. Синтез модели будем осуществлять методом последовательных итераций, с увеличением на каждом шаге лага модели на единицу. Рассмотрим перспективное предсказание (экстраполяцию вперед) и ретроспекцию (экстраполяцию назад). Эти экстраполяции будем помечать символом "€". Согласно рекурсивному методу Теплица, сначала будем находить уравнения авторегрессии (3)-(4) единичного лага, затем, используя найденные соотношения, строим модель с задержкой на 2 единицы, и т.д. пока не будет достигнута желаемая точность аппроксимации, или получаемые приближения уже не будут дальше улучшаться. Достигнутый при этом шаг итерации примем за эффективный лаг модели. Для реализации этой идеи на каждом шаге итерации будем вычислять достигнутую погрешность прогноза. В заключение отметим также, что рекурсия Теплица является классической рекурсией, встречающейся в теории полиномов, ортогональных на единичном круге. Итак, на i-ом итерационном шаге имеем


X€_bi(t_k) = X_pi(t_k-i+1) C_bi + B_i(t_k),

X€_pi(t_k-i+1) = X_bi(t_k) C_pi + P_i(t_k-i+1). (5)


Здесь C_bi, C_pi - матрицы частных производных размерности m´m (переходные матрицы); значки "_b", "_p" обозначают принадлежность матрицы соответственно для перспективного и ретроспективного прогноза; второй нижний индекс у этих матриц означает номер итерации; B_i(t_k), P_i(t_k) – вектора квазислучайных невязок с матожиданием Е(B_i) = E(P_i) = 0 и размера m; k = i,(i+1),…,N; i = 2,3,…,l - номер итерации; 0 - нулевой вектор размера m.

Из формулы (5) видно, что число отсчетов между прогнозируемой величиной и предиктором равно |i-1|, поэтому начальная итерация имеет номер i_н = 2. Определим начальные условия (при i = 1) следующим образом: C_b1 = C_p1 = E, где E – единичная матрица размера m´m. Выражение (5) является одношаговым предсказателем, использующим для прогноза текущего значения вектора одно его прошлое (будущее) значение, сдвинутое на (i-1) временных тактов. Представим соотношения (5) в матричном виде


X€_bi = X_piC_bi + B_i,

X€_pi = X_biC_pi + P_i, (6)


где X_bi, X_pi, B_i, P_i – матрицы размера (N-i)´m.

Проблема состоит в определении оптимальных (по некоторому критерию) значений матриц C_bi, C_pi. При (N-i) < m проблему (6) можно рассматривать как типичную задачу метода наименьших квадратов (МНК) с критерием оптимальности в виде суммы квадратов невязок.


J1 = min{(X_bi - X_pi C_bi)¢ (X_bi - X_pi C_bi)},

J2 = min{(X_pi - X_bi C_pi)¢ (X_pi - X_bi C_pi)}, (7)


где символом " ¢ " обозначена операция транспонирования.

Нам необходимо найти такие оценки C_bi, C_pi параметров модели (6), которые удовлетворяли бы этому условию. Для раскрытия экстремума функций J1(C_bi), J2(C_pi) воспользуемся правилами векторного дифференцирования (Рао, 1968). Пусть Z, A, B – матрицы или вектора. Тогда имеют место следующие соотношения:


¶(B¢Z)/ ¶Z = ¶(Z¢B)/ ¶Z = B, ¶(Z¢Z)/ ¶Z = 2Z,

¶(B¢AZ)/ ¶Z = A¢B, ¶(Z¢AZ)/ ¶Z = 2AZ. (8)


Используя эти соотношения, находим необходимое условие экстремума


¶(J1)/ ¶Z = - X¢_biX_pi + X¢_piX_piC_bi = 0,

¶(J2)/ ¶Z = - X¢_piX_bi + X¢_biX_biC_pi = 0,


откуда имеем

X¢_piX_piC_bi = X¢_biX_pi,

X¢_biX_biC_pi = X¢_piX_bi. (9)


Условием минимума квадратичных форм J1, J2 означает, что матрицы


¶²(J1)/¶Z² = X¢_piX_pi, ¶²(J2)/ ¶Z² = X¢_biX_bi


должны быть положительно определены.

Выражение (9) представляет собой матричную запись так называемой нормальной системы уравнений. Если матрицы X_bi, X_pi неособенные, то системы (9) имеют единственные решения


C_bi = (X¢_piX_pi)^-1X¢_biX_pi,

C_pi = (X¢_biX_bi)^-1X¢_piX_bi, (10)


где символом "^-1" обозначена операция вычисления обратной матрицы.

Множители правой части соотношения (10) представляют матрицы дисперсий и взаимных ковариаций между будущими и прошлыми значениями вектора состояний. Как следует из формулы (5), дисперсии, вычисляемые по предыдущему выражению, рассчитываются на разных временных базах (для перспективного прогноза используется база t₁ £ t_k £ t_N-i+1, а для ретроспективного - t_i+1 £ t_i £ t_N). Временная несовмещенность этих баз зависит от номера итерации и в целом составляет 2i отсчетов данных. Взаимные ковариации, вычисляемые для различных временных задержек, также пропорциональны номеру текущей итерации. Следовательно, при увеличении лага модели уменьшается число используемых в вычислениях данных.

Упростим выражение (10). Для этого введем следующие обозначения, смысл которых легко установить при анализе предыдущих формул:


X¢_biX_pi = R_bpi,

X¢_piX¢_bi = R_pbi,

X¢_piX_pi = D_pi,

X¢_biX_bi = D_bi . (11)


Тогда соотношения (10), с учетом обозначений (11), могут быть переписаны в более компактном виде


C_bi = (D_pi)^-1 R_bpi,

C_pi = (D_bi)^-1 R_pbi. (12)


Отметим, что s-ая строка матриц C_bi, C_pi является строкой весовых коэффициентов для прогноза s-ой составляющей текущего фазового вектора по информации, задержанной на (i-1) временных тактов.

Известно, что решение (12) может иметь вычислительные трудности при обращении матриц D, если они вырождены или плохо обусловлены. Это является достаточно серьезным ограничением, которое может разрушить сходимость итерационного процесса. В работе (Химмерблау, 1975) изложен способ, позволяющий улучшить обусловленность матрицы D, в котором предлагается обращать матрицы вида (D + lE), где E – единичная матрица, l - параметр регуляризации. Известна также модификация приведенного выше метода, в котором обращается следующая матрица (D + lD_d), где D_d – диагональная матрица, составленная из соответствующих элементов матрицы D. Величина параметра l находится либо путем оптимизации его значения, либо задается априорно. В последнем случае параметр регуляризации l обычно принимается равным 0.01¸0.05, т.е. порядка нескольких процентов. Другим способом решения этой проблемы является сингулярный анализ (Троян, Соколов, 1989).

Имея оценки переходных матриц (12), можно выполнить прогноз фазового вектора по следующим формулам:


B_ki = P_(k-i+1)(i-1)C_bi,

P_(k-i+1)i = B_k(i-1)C_pi, (13)


где первый индекс у векторов B, P означает дискретное время, а второй – номер итерации. В результате получим ошибки предсказания и ретроспекции, которые следуют из системы уравнений (5):

B_ki = B_k(i-1) - P_{(k-i+1) (i-1)} C_bi,

P_(k-i+1)i = P_(k-i+1)(i-1) - B_k(i-1) C_pi (14)


при начальных условиях B_k1 = X(t_k), P_k1 = X(t_k), k = 1,2,...,N.

Поясним физическую сущность введенных выше соотношений. С учетом формулы (11) и невязок (14), матрицы D_bi, D_pi интерпретируются как ковариационные матрицы ошибок предсказания и ретроспекции (с размерностями m´m). Вектора B_i и P_i интерпретируются, соответственно, как ошибки предсказания и ретроспекции на i-ом шаге после учета линейных связей предыдущих (i-1) итераций. Поэтому матрицы R_bpi и R_pbi представляют взаимные функции частных автоковариаций составляющих фазового вектора (ВФЧАК) предсказания и ретроспекции на (i-1)-ом сдвиге. И в этом случае, как видно из формул (12), матрицы C_bi, C_pi можно трактовать как нормированные значения ВФЧАК. Мы используем понятие частных автоковариаций, чтобы подчеркнуть тот факт, что на i-ой временной задержке учитываются только те линейные связи, которые еще остались после устранения линейности для более ранних запаздываний.

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

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

Для стационарного процесса теоретические дисперсии являются постоянными, а между ВФЧАК существуют простые зависимости (Клаербоут, 1981). В целом все эти условия можно описать следующими соотношениями:

D_b = D_p, (15¢)

R_bp = R'_pb, (15¢¢)

C_b = C ¢_p. (15¢¢¢)


Рассмотрим одномерный вариант метода. В этом случае матрицы из выражений (15¢-15¢¢¢) вырождаются в числа, при этом ВФЧАК трансформируется в ФАК. Известно, что ФАК в одномерном случае является четной функцией, и, следовательно, для любого процесса (с трендом или без него) выполняется равенство R_b = R_p. При практической реализации алгоритмов в классическом варианте выполнение соотношения (15¢) достигается за счет усреднения дисперсий ошибок предсказания и ретроспекции. Анализ формулы (12) показывает, что при выполнении равенств (15¢-15¢¢), автоматически выполняется и соотношение (15¢¢¢). Многомерный стационарный случай не имеет в этом отношении принципиального отличия от одномерного, следовательно, и для него выполняются условия (15¢-15¢¢¢).

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

Этот вывод еще в большей степени применим к многомерному случаю, как к объекту более сложной природы. Оценка и учет этого расхождения, по нашему мнению, позволяют смоделировать параметры возможного тренда. Практическая реализация этой идеи заключается в раздельной оценке матриц D_b и D_p по не полностью совмещенным временным интервалам. Степень разноса данных зависит как от лага модели l, так и от длины базы исходных данных N.

Итак, после второй итерации получаем матрицы C_b2, C_p2 и вектора остатков B₂ и P₂, а также оценки их дисперсий D_b2 и D_p2 и ковариаций R_bp2, R_pb2. При необходимости выполняем третью итерацию. Для этого рассчитаем по формуле (12) матрицы C_b3, C_p3. При этом вместо векторов X_b2, X_p2 используем вектора остатков B_k2 и P_k2 соответственно, полученные по формулам (14). После применения выражения (13), а затем (14) к этим векторам вновь получаем остатки B_k3 и P_k3, но уже после третьей итерации и т.д. На l-ом шаге получаем матрицы С_bl и С_pl, вектора невязок B_kl, P_kl, оценки их дисперсий и ковариаций.

Совокупность матриц C_bi, C_pi является достаточной для построения оценок фазового вектора, однако практическим неудобством при этом является необходимость использовать не исходные данные, а их разности. Для явного определения параметров моделей (3)-(4) требуется вычислить матрицы A_bi, A_pi по ранее полученным матрицам C_bl, C_pl, т.е. нужно перейти от частных моделей к обобщенным. Для этого воспользуемся рекуррентной формулой Теплица для одномерного случая (Никитин, 1979), которая позволяет по коэффициентам модели (l-1) порядка получить коэффициенты модели l-ого порядка, если известны матрицы C_bl, C_pl. Обозначим матрицы моделей (l-1) порядка следующим образом: A_bi(l-1), А_pi(l-1), i = 1,2,…,(l-1). В этом случае матрицы моделей (3-4) l-ого порядка, по аналогии с одномерным случаем, могут быть вычислены по следующим формулам:


A_bl = C_bl  - C_blA_bi(l-1) i = 1,2,…,(l-1),

A_pl = C_pl  - C_plA_pi(l-1) (членов нет, если l = 1). (16)


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

Рассмотрим этот вопрос более подробно.

Будем анализировать эволюцию информации по мере увеличения лага модели. Если принятая к рассмотрению модель адекватна анализируемому процессу, то при увеличении ее порядка l элементы матриц C_bl, C_pl (а также A_bl, A_pl) должны приближаться к нулю, а след матриц D_bl, D_pl – стремиться к постоянным величинам. Большие остаточные значения составляющих матриц С_bl, C_pl свидетельствуют или о малости принятого лага модели, или о нелинейности связей между исходными данными, или о неполном информационном базисе. Последовательно увеличивая лаг модели, можно уменьшить эту неопределенность. Матрицы D_bl, D_pl характеризуют немоделируемую часть исходных данных, другими словами, отображают шум модели. Стабилизация оценок этих матриц свидетельствует о том, что достигнута предельная точность выбранной модели, и дальнейшее увеличение ее лага уже не может привести к существенному улучшению оценок, а иногда и ухудшает их. Таким образом, анализ поведения матриц в зависимости от лага модели позволяет определить оптимальное значение порядка разностных уравнений (3-4).

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

Аналогичным образом можно проинтерпретировать стремление к нулю элементов матриц C, A по мере увеличения лага модели. Известно, что белый шум имеет ковариационную зависимость типа дельта-функции, другими словами, два его различных временных отсчета статистически не связаны, и, следовательно, он не может нести какой-либо прогностической информации, что и отражается в поведении матриц C, A. Рассмотрим эту проблему с несколько другой стороны. Будем считать, что выражения (3-4) представляют модели с постоянными коэффициентами, т.е. предположим, что матрицы A_bl, A_pl не меняются во времени. Это предположение соответствует стационарному приближению, для которого будет выполняться равенство A_b = A¢_pl. Таким образом, для классического метода наблюдается равноправность, с точки зрения качества прогноза, прямого и обратного времен. Для рассматриваемой методики, как видно из приведенных выше замечаний, эта равноправность должна нарушаться даже для модели с постоянными коэффициентами. Это нарушение будет проявляться в виде некоторого рассогласования аналогичных матриц, построенных для прямого и обратного времен. Следовательно, будут несколько отличаться и соответствующие прогностические модели и качество их работы. Степень рассогласования аналогичных матриц, или качества прогноза, может являться характеристикой нестационарности процесса. Но особенно хорошим индикатором тренда являются нормированные ФАК, имеющие ясный физический смысл и наиболее простые и наглядные для анализа. Поэтому рассмотрим задачу расчета корреляционных матриц по полученной выше информации.

Покажем далее, что знание матриц A_bl, A_pl позволяет оценить функции автокорреляции фазового вектора. Для одномерного случая Дж. Бургом показано, что эти параметры связаны следующим матричным соотношением (Никитин, 1979):

R(0) R(1) … R(l)		1		Dl
R(1) R(0) … R(l-1)		A1l	=	0
… … … …		…		…
R(l) R(l-1) … R(0)		All		0	(17)