CC BY
152
УДК 681.51 ББК 32.811.1
ОБНАРУЖЕНИЕ СТРУКТУРНО-ПАРАМЕТРИЧЕСКИХ ИЗМЕНЕНИЙ В СТОХАСТИЧЕСКИХ СИСТЕМАХ В РЕАЛЬНОМ МАСШТАБЕ ВРЕМЕНИ АЛГОРИТМАМИ НЕПРЕРЫВНЫХ ДРОБЕЙ И СТРУКТУРНОГО АНАЛИЗА
1 2 Карташов В. Я. , Новосельцева М. А.
(ГОУВПО «Кемеровский государственный
университет», Кемерово)
В статье разработан алгоритм обнаружения структурнопараметрических изменений в стохастических процессах и системах в режиме реального времени. Алгоритм основан на совместном использовании структурного анализа и теории непрерывных дробей.
Ключевые слова: структурная функция, непрерывная дробь, идентифицирующая матрица.
1. Введение
Обработка и анализ временных рядов представляет собой интенсивно развивающееся направление в математических, технических, экономических, политических, социальных и других науках. Традиционным при анализе временных рядов является предположение о том, что статистические свойства наблюдаемого ряда или свойства порождающей его системы
1 Владимир Яковлевич Карташов, заведующий кафедрой автоматизации исследований и технической кибернетики, доктор технических наук, профессор (kartash@kemsu.ru).
2 Марина Александровна Новосельцева, доцент кафедры автоматизации исследований и технической кибернетики, кандидат технических наук, доцент (aanov@pochta.ru).
сохраняют определенное постоянство во времени. Однако многие практические задачи текущего контроля производства, экономики, медицинской и технической диагностики, гидроакустики, геофизики и др. связаны с оценкой риска наступления той или другой непредвиденной ситуации, вызванной некоторым устойчивым изменением свойств наблюдаемого временного ряда, происходящим в неизвестный момент времени.
Интерес к проблематике данных задач стал возрастать с середины 60-х годов, что вызывалось потребностями практических приложений и появлением ЭВМ, давшей возможность автоматизировать данный процесс. При этом основные усилия исследователей направлялись на то, чтобы разработать методы, использующие как можно меньше априорной информации о наблюдениях и свойствах временного ряда. Популярность данной проблематики несомненна.
2. Постановка задачи
Определим круг практических задач и их постановки, в которых возникает потребность обнаруживать изменение свойств случайных процессов [7, 10, 12, 26].
Первый класс рассматривает ретроспективные (апостериорные) задачи. Пусть представлена реализация случайного процесса. Корректность методов анализа основных свойств случайных процессов, а также интерпретация результатов анализа в значительной степени влияют на правильность получения модели, оценку ее параметров и прогнозирование будущих значений временного ряда. Поэтому статистическая обработка реализации случайного процесса должна быть основана на предположении, что его свойства в процессе сбора данных не изменялись. Предварительным этапом любой статистической обработки должен быть этап проверки подобной неизменчиво-сти. Если характеристики не изменялись, то следует заниматься статистической обработкой в зависимости от целей исследователя. В противном случае возникает задача обнаружения моментов изменения статистических характеристик и разбиения
исходной выборки на несколько интервалов стационарности и однородности. Такая задача решалась авторами в [19, 21, 27].
Второй класс непосредственно связан с ситуацией обработки измерений случайного процесса в реальном масштабе времени. Допустим, что в некоторый заранее неизвестный момент времени происходит изменение какой-либо статистической характеристики процесса. Как обнаружить произошедшее изменение скорейшим образом после того, как оно возникло? В данном случае запаздывание в обнаружении момента изменения свойств временных рядов может привести к наступлению рисковой ситуации в системе [17, 28]. Эта задача часто формулируется как задача обнаружения момента разладки процесса [9, 12,
26, 29]. Именно ее решению будет посвящено дальнейшее изложение.
3. Обзор методов обнаружения изменения свойств стохастических процессов и систем
Существующие методы решения данной задачи позволяют достаточно эффективно улавливать момент изменения характеристик случайного процесса. Однако развитая теория имеет своим основанием сложный математический аппарат, не самый простой с точки зрения практического применения, поскольку его можно использовать, имея дополнительную априорную информацию о данных измерений [4, 5, 7, 9, 22]. В качестве априорных знаний или предположений при обнаружении изменений свойств стохастических сигналов и систем обычно используют: знание распределения момента появления изменений [33]; знание законов или классов распределений данных измерений в нормальном состоянии и после появления изменений [22, 29, 33]; знание вида модели процесса до и после изменений с указанием меняющихся параметров процесса [29]; задание величины измененного параметра либо пороговых значений процесса [7, 22, 29]; стационарность данных измерений до и после разладки [4, 29] и т.д.
В связи с отсутствием какой-либо априорной информации при обнаружении момента изменения свойств временных рядов
желательно применение непараметрических критериев, не использующих априорную информацию о распределениях. Кроме того, ограничения области применимости критериев одним типом разлаженного состояния могут привести к запаздыванию в принятии решения по управлению стохастическими объектами и системами [18].
Алгоритм кумулятивных сумм (АКС) [26, 29] является наиболее популярным, теоретически обоснованным и исследованным методом обнаружения изменений в последовательности наблюдений. Данный метод обобщен на различные типы данных измерений и имеет несколько модификаций: для независимых отсчетов, для зависимых отсчетов, для линейных объектов, для процессов авторегресии проинтегрированного скользящего среднего и т.д. Однако методика АКС имеет ряд недостатков, заключающихся в следующем.
1. Использование разного рода априорной информации, часто недоступной на практике для малоизученных явлений и при работе с новыми уникальными объектами.
2. Выбор шага дискретизации случайного процесса не ясен, подчеркивается его важность, но рекомендаций и правил его выбора не дается.
3. Использование малоэффективных корреляционных методик визуального анализа стационарности методом Бокса-Дженкинса, недостатки которых подробно описаны в [21, 27].
4. Идентификация временных рядов на основе перебора пробных моделей методикой Бокса-Дженкинса. Подробный анализ этой методики приведен в [16, 21, 27]. Авторы сами отмечают [26], что на качество работы алгоритма влияет «...неточность восстановления модели авторегрессии», заключающаяся в неверном задании порядка или определении параметров модели.
5. Метод основан на вычислении логарифма правдоподобия авторегрессионной модели, на корреляционной теории случайных процессов, которые применимы лишь в стационарных условиях.
Необходимость наличия большого количества априорной информации приводит к ограничению использования АКС.
Кроме того, в некоторых условиях нет возможности проводить длительную настройку алгоритма по какому-то конкретному параметру даже с наименьшей информацией о нем и с уточняющими процедурами. Следовательно, актуальной становится задача применения таких методов обработки данных измерений, которые не включали бы в себя вышеописанные недостатки. Представляется, что в таких условиях наиболее эффективными окажутся непараметрические методы, не использующие значительное большинство из перечисленных выше априорных сведений.
Одним из таких критериев является инверсионный алгоритм разделения данных на интервалы стационарности [16, 24, 25]. Алгоритм является одним из наиболее эффективных непараметрических критериев обнаружения изменений характеристик случайных процессов в условиях априорной неопределенности. В [16, 27] данный алгоритм был применен для ретроспективного анализа стационарности случайных процессов при идентификации стохастических объектов. Метод может использоваться в реальном масштабе времени и позволяет обнаружить изменения среднего, дисперсии, а также появление нестацио-нарности. Недостатком метода, как отмечает сам автор [24, 25], является то, что данный метод не позволяет выявлять момент смены корреляционной функции процесса, т.е. появления структурной неоднородности случайного процесса [16, 27].
4. Использование структурного анализа и теории непрерывных дробей для обнаружения структурно-параметрических изменений стохастических сигналов и систем
Для разработки критерия обнаружения изменений структурно-параметрических свойств процессов и стохастических систем воспользуемся структурным анализом - разделом теории случайных процессов, связанным с изучением последних на основе структурных функций [16, 30]. Структурный анализ случайных процессов в большинстве случаев приводит к более устойчивым характеристикам по сравнению с корреляционны-
ми. Эффективность структурного анализа заключается в том, что параметры структурной функции обладают свойствами инвариантности относительно некоторых форм нестационарно-сти, проявляющихся, например, при смещенности по математическому ожиданию, а также в случае квазистационарного характера случайного процесса.
Структурная функция [16, 30], предложенная А.Н. Колмогоровым, определяется как математическое ожидание квадрата разности сечений случайного процесса, соответствующим значениям аргумента t и t + т:
(1) Сх ^, t + т) = М{ х^) - х^+т) }2,
где х(0 - некоторый случайный процесс. Очевидно, что функция всегда неотрицательна, четна и удовлетворяет условию Сх(^ t) = 0.
Рассмотрим частный случай, когда х{() представляет собой стационарный случайный процесс с нулевым средним. Тогда структурная функция будет иметь вид
Сх ^, t + т) = М{ х(^) - х^ + т) }2 =
(2) = М (x2(t)) - 2М (х^) х^ + т)) +
+М (х+ т)) = 2Яхх (0) - 2#хх (т) = Сх (т).
Как видно из (2), структурная функция стационарного случайного процесса не зависит от текущего момента времени. На основании (2) можно утверждать, что структурная функция стационарного случайного процесса с течением времени стремится к установившемуся значению:
(3) Сх (т) ——® 2 Яхх (0).
Практическое построение структурной функции более надежно по сравнению с корреляционной, поскольку на нее не влияют ошибки определения среднего значения процесса х{() [30]. Оно осуществляется по формуле
1 N -к
(4) Сх (к) = —— X (х(0 - х(/ + к))2 ,
N - к =
где N - число измерений процесса х(^.
Проведенные ранее исследования [16, 19, 21, 27] позволили использовать структурную функцию для анализа стационарности случайного процесса, т.е. для решения ретроспективной задачи. Поскольку данная методика ляжет в основу алгоритма обработки реализации случайного процесса в реальном масштабе времени, опишем ее более подробно.
Для анализа стационарности использовались следующие закономерности [16, 19, 27]:
- структурная функция стационарного процесса ограничена и с течением времени выходит на установившееся значение;
- структурная функция нестационарного процесса с течением времени неограниченно возрастает.
С целью обеспечения эффективной алгоритмической и программной реализации метода анализа стационарности, а также исключения визуального анализа кривых структурной функции был разработан алгоритм получения модели структурной функции сигнала на основании ее значений [16, 19, 21, 27]. Кроме того, свойство стационарности случайного процесса отождествлялось со свойством устойчивости полученной модели структурной функции.
Алгоритм основан на использовании нетрадиционного математического аппарата - теории непрерывных дробей [15]. В частности наибольшее развитие в этой теории получил класс правильных С-дробей следующего вида:
(5)
Ь • —п-
и0 5 1
- Ь0 +'
1 +---2----
1 +
1+...
где аь ..., ап, Ь0 - постоянные величины.
Наиболее приемлемым и простым способом перехода к непрерывной дроби является модифицированный алгоритм
В. Висковатова [15, 13].
При нахождении структурной функции используется виртуальный подход [16], заключающийся в том, что процесс определения структурной функции Сх(к) можно представить как реакцию динамического виртуального объекта на единичное
входное воздействие. В этом случае Сх(к) для фиксированных к порождает передаточную функцию виртуального объекта, знаменатель которой представляет собой его характеристическое уравнение. Исследуя это уравнение на устойчивость любым известным методом, делается вывод о стационарности или нестационарности случайного процесса.
Критерий проверки стационарности случайного процесса х(?) с помощью структурной функции включает в себя следующие этапы.
1. На основании значений случайного процесса х(?) вычисляются значения структурной функции по формуле (4).
2. На основании значений структурной функции процесса определяется идентифицирующая матрица:
(6)
(-1)строка 1 1 1 1 1 1
0 строка Сх(0) С(А) Сх(2А?) Сх(иА?)
1 строка а:(0) а^Л?) а^А?) а^иА?)
2 строка а2(0) а2(Л?) а2(2А?) а2(иА?)
т строка ат(0) ат(А?) ат(2А?) а ( > )
в которой (-1)-строка содержит значения единичной функции 1(0, а (0)-строка - значения структурной функции входного сигнала Сх(кА?) в моменты времени {иЛ?^ ; А? - шаг дискретизации, а элементы ат(иА?) последовательно определяются с помощью соотношения:
(7) ат (пЛ) =
ат _ 2 ((и + 1)Л?) ат _1((и + 1)Л?)
ат_2 (0) ат_1(0)
где а_1(и) = 1(иЛ?); а0(и) = Сх(иЛ?); т = 1, 2, 3 ,..., и = 0, 1, 2, ... Вычисление элементов ат (иЛ?) продолжается до появления нулевой строки.
3. Элементы первого столбца матрицы (6) порождают частные числители правильной С-дроби [13, 15], что позволяет получить модель структурной функции сигнала в форме дискретной передаточной функции (ДПФ) виртуального объекта:
(8) Осх (2) =
С. (1) 2-1
1 +
а1 (0)2
1 + а2(0) 2_1
1+
а 2(0) 2 -1
Если в некоторой /-й строке (/ = 0, 1, 2, ...) матрицы (6) конечное число к/ первых элементов равны нулю, а последующие элементы отличны от нуля, то необходимо осуществить сдвиг влево на к/ элементов до появления в нулевом столбце ненулевого элемента и далее продолжить определение других элементов матрицы (6) по правилу (7). Для /-й строки при восстановлении модели структурной функции элемент а/ умножается на г_к‘ .
4. Производится проверка устойчивости полученного виртуального объекта по математической модели ДПФ (8) с помощью любого из известных критериев устойчивости [31]. В случае устойчивости объекта следует утверждать, что данный сигнал стационарен. В противном случае сигнал является нестационарным.
Пример 1. Пусть х(0 - стационарный случайный процесс с нулевым средним и корреляционной функцией До;(0 = е-0’3 Структурная функция процесса имеет вид Сх(т) = 2 - 2е~°’3т. На основании значений структурной функции, взятых с шагом дискретизации А( = 1, построим идентифицирующую матрицу (6):
1
1
1 1 1 1 1 1
0,51836 0,90238 1,18686 1,39761 1,55374 1,66940
-0,74082 -1,28963 -1,69620 -1,99739 -2,22052
0 0 0 0
С учетом того, что в нулевой строке был осуществлен сдвиг влево на один элемент, получаем модель структурной функции в форме ДПФ виртуального объекта
^ , ч 0,51836г-1
^С (г) — т.
С^ 1 - 0,74082г-
ДПФ объекта имеет один полюс гп = 0,74082, который лежит внутри единичной окружности с центром в начале координат плоскости г. Следовательно, объект устойчив, а это означает, что случайный процесс х(0 является стационарным.
Чтобы убедиться в том, что структурная функция может служить индикатором структурной неоднородности, рассмотрим пример.
Пример 2. Имеется стохастический процесс у (0, взятый из [3]:
(9) у'(0 —10 + х^) + х(Г -1), t — 1,...,50,
где х(0 - белый шум с нулевым средним и единичной дисперсией.
В момент времени t = 51 происходит изменение параметров модели случайного процесса, и она принимает вид [3]:
(10) у2(0 —10 + х(0-х^-1), t — 51,...,100 .
Графики процессов (9) и (10) представлены на рис. 1. Нетрудно показать, что математические ожидания и дисперсии процессов у'(0 и у2(0 одинаковы:
М (у1^)) — М (у2^)) —10,
Б( у^)) — Б( у2^)) — 2, поэтому применение инверсионного критерия [16, 24, 25] не является эффективным.
14
12
10
8
6
4
2
0
у№
11 21 31 41 51 61 71 81 91
Рис. 1. Реализации случайных процессов (9) при t = 1, ..., 50 и
(10) при t = 51, ..., 100
Корреляционные функции процессов (9) и (10) соответственно равны
(11) куЛ%) =
% = 0, % = 1, % > 2;
(12) Яу 2(%) =
2, т — 0,
-1, т — 1,
0, т > 2.
Тогда структурные функции процессов (9) (рис. 2а) и (10) (рис. 2б) соответственно равны:
0, т — 0,
2, т — 1,
4, т > 2;
0, т — 0,
6, т — 1,
4, т > 2.
(13) Су1 (%) =
(14) С у 2(%) =
Рис. 2. а) Структурная функция (13), б) структурная функция (14)
Определим структурную функцию стационарного случайного процесса на всем интервале і = 1, ..., 100 (рис. 3.) по формуле (4). Это доказывает, что объединение двух стационарных реализаций с разными структурными функциями приводит к
стационарной реализации с отличной от исходных структурной функцией. Аналогичное явление будет постепенно происходить при добавлении значений случайного процесса в режиме реального времени.
Рис. 3. Структурная функция объединенного процесса г = 1, ..., 100
Методы экспериментального исследования характеристик объектов, явлений, процессов, как правило, предусматривают проведение определенного числа опытов с фиксацией результатов [11]. Лишь после окончания последнего опыта, используя тот или иной алгоритм, находится нужная характеристика. Это относится и к структурной функции, определяемой на практике по формуле (4). В системах реального времени с использованием мониторинга наблюдений временных рядов необходимо отслеживать изменения их свойств на основе итерационных характеристик. Если какая-либо характеристика изменяется, то итерационные алгоритмы «следят» за этими изменениями [11, 23]. Так, например, итерационная формула для вычисления начального 5-го момента [11] после Л-го измерения имеет вид:
(15) тЛ (х) = тЛ_ (х) + -1 (хЛ - тЛ_ (х)).
Получим рекуррентную формулу для вычисления структурной функции случайного процесса. Подставляя в данное
равенство выражение, задающее структурную функцию (1), получим:
(16) Сых (к) = Сых -1(к) + -1— ((( - х(N - к))2 - Сых -1(к)],
N - к
где С^!-1(к) - оценка структурной функции, полученная на основе N - 1) измерения, к = 0, 1, 2, ...
Для проверки полученного равенства используем основную формулу (4) практического определения структурной функции. Раскроем скобки в полученном выражении (16):
1 N-1-к
СN (к) = ^-ГТ £ () - х(г' + к) )2 +
N - к -1 г=1
1 Г 1 N-1-к 1
—-[(х(^) - х(М - к))2-—— £ (х(/) - х(2 + к))2 ] =
N — к N — к — 1 2 =1
1 N-1-к 1
= ^ТГ £ Х) - х(2 + к) )2 +—— (х( N) - х( N - к))2 =
N — к 2=1 N — к
1 N - к
= —— £ Х) - х(2 + к))2 .
N - к 2=1
Таким образом, итерационная формула структурной функции тождественна вычислению структурной функции по N измерениям.
Использование итерационной формулы расчета структурной функции дает возможность получать сведения об изменении характеристик процесса в реальном масштабе времени. Все преимущества структурного анализа продолжают иметь место, кроме того, структурная функция может служить индикатором изменения как математического ожидания, дисперсии, так и корреляционной функции случайного процесса, а также наличия в нем нестационарной компоненты.
Разовьем идею ретроспективного анализа стационарности случайного процесса [16, 19, 21, 27] на случай его применения в реальном масштабе времени.
Пусть имеется реализация случайного процесса в режиме нормального функционирования системы. Шаг дискретизации А( выбран исходя из условия структурно-параметрической идентифицируемости [16, 27]. Необходима предварительная
проверка данной реализации на стационарность и однородность. Кроме того, следует получить идентифицирующую матрицу и модель виртуального объекта для дальнейшего контроля и обнаружения изменений свойств данных измерений. Далее следует осуществлять контроль времени начала возможных изменений внутренних свойств виртуального объекта на основе значений выходной структурной функции. Для этого осуществляется перерасчет структурной функции случайного процесса по итерационной формуле (16). Поскольку наличие нулевой строки в идентифицирующей матрице связано с определением порядка функции и показывает механизм аппроксимации непрерывными дробями ДПФ объекта [13], то различные изменения, которые могут появляться в нулевой строке, будут свидетельствовать об изменениях, происходящих в объекте. Теория непрерывных дробей относится к разделу алгоритмической математики, поэтому теоретически доказать этот факт не представляется возможным. В [13] данное утверждение доказано для апериодического звена 1-го порядка. Для объектов 2-го порядка справедливость утверждения доказана на ЭВМ. Из-за чрезвычайной громоздкости доказательства в общем случае данный факт можно показать лишь на тестовых примерах.
Таким образом, для решения задачи обнаружения структурно-параметрических изменений в реальном масштабе времени будем использовать итерационный алгоритм, в котором при добавлении каждого 2-го измерения процесса производится перерасчет выходной структурной функции, т.е. уточняется ее значение. Если свойства объекта, формирующего на выходе структурную функцию, изменяются, то будет меняться и его модель, следовательно, нулевая строка идентифицирующей матрицы перестает быть таковой [13].
При решении практических задач важным моментом является формирование критериев определения нулевых строк матрицы (6), т.е. определение модели виртуального объекта. На основе модельных исследований были установлены качественные признаки выявления нулевых строк. Таковыми являются [14]:
- близость к нулю элементов строки матрицы;
- изменение элементов строки по величине и знаку;
- резкое увеличение по модулю элементов следующей строки.
Возможна также дополнительная проверка гипотезы о равенстве нулю среднего значения элементов нулевой строки идентифицирующей матрицы с помощью любого непараметрического критерия [2, 6].
Если при добавлении новых данных измерений в идентифицирующей матрице нулевая строка сохраняет свое местоположение, полученная модель виртуального объекта остается неизменной, следовательно, с некоторой заданной вероятностью можно утверждать, что структурно-параметрические изменения отсутствуют. Если же при добавлении новых измерений нулевая строка в матрице отсутствует либо некоторые (или все) ее элементы резко возрастают (по модулю), полученная ранее модель не является адекватной виртуальному объекту. Таким образом, с некоторой заданной вероятностью рд можно утверждать, что в объекте произошли параметрические и (или) структурные изменения, связанные с изменениями характеристик временного ряда.
5. Алгоритм обнаружения
структурно-параметрических изменений стохастических процессов и систем
Проведенные исследования позволили сформулировать алгоритм обнаружения изменения характеристик случайного процесса х(0 с помощью структурной функции, который включает в себя следующие этапы.
Подготовительный этап
А) Имеется N - 1) измерение случайного процесса х(/А0 в режиме нормального функционирования системы, где А2 - шаг дискретизации, 2 = 0, ..., N - 1. На данном этапе также необходимо осуществить проверку измерений х(/АО на стационарность с помощью указанных выше критериев.
Б) На основе измерений случайного процесса х(/А) осуществляется вычисление значений структурной функции для количества отчетов N - 1) по формуле (4).
В) На основании значений структурной функции сигнала рассчитывается идентифицирующая матрица (6). Вычисления продолжаются до появления нулевой строки в матрице (обозначим ее за т-ю). Затем определяется модель виртуального объекта вида (8).
Г) Осуществляется проверка значений нулевой строки на независимость с заданным уровнем значимости а. Если гипотеза о независимости не отвергается, то строится а-процентный доверительный интервал (е(к) - а; е(к) + а) для значений нулевой строки. Если гипотеза отвергается, то данная строка не является нулевой, и получить модель виртуального объекта не представляется возможным. В этом случае вопрос определения модели виртуального объекта и структурной функции на его выходе остается открытым.
Этап контроля возможных изменений случайного процесса
A) При добавлении Л-го наблюдения оценка структурной функции пересчитывается по рекуррентной формуле (16).
Б) Осуществляется перерасчет идентифицирующей матрицы (6).
B) Производится оценка покрытия элементов т-й строки идентифицирующей матрицы а-процентным интервалом. Если хотя бы одно из значений вышло за границы интервала, с вероятностью рд = 1 - а можно утверждать, что модель виртуального объекта изменилась. Последнее означает смену какой-либо характеристики случайного процесса (математического ожидания, дисперсии, корреляционной функции, появления нестацио-нарности). В случае покрытия всех элементов нулевой строки а-процентным интервалом можно утверждать, что характеристики случайного процесса не изменялись.
6. Примеры использования алгоритма для обнаружения структурно-параметрических изменений стохастических процессов и систем
Пример 3. Имеется 250 измерений стационарного случайного процесса, заданного конечно-разностным уравнением:
(17) х(к) = 0,9х(к _ 1) + п(к),
где п(к) - белый шум с нулевым математическим ожиданием и дисперсией, равной 0,25. Получим значения структурной функции случайного процесса по формуле (4). На основе этих значений рассчитаем идентифицирующую матрицу (6):
№ столбца 1 2 3 4 5 6
(-1)строка 1 1 1 1 1 1
0 строка 0,2821 0,5387 0,7588 0,9771 1,1426 1,2870
1 строка -0,9099 -1,6899 -2,4639 -3,0507 -3,5626 -4,1004
2 строка 0,0527 -0,0181 0,1112 0,1354 0,0561 -0,1124
№ столбца 7 8 9 10 11 12
(-1)строка 1 1 1 1 1 1
0 строка 1,4387 1,6200 1,7822 1,4387 1,6200 1,7822
1 строка -4,7432 -5,3179 -5,8010 -6,086 -6,3733
2 строка -0,1013 -0,0575 0,1124 0,0816
Вторая строка матрицы является нулевой. С учетом того, что в нулевой строке был осуществлен сдвиг влево на один элемент, получаем модель структурной функции в форме ДПФ виртуального объекта
^ (_ ч = 0,2821г-1
^С (г) — 1 .
С^ 1 - 0,9099г-
Произведем проверку гипотезы о независимости значений e(k), стоящих в нулевой строке и представляющих собой погрешности вычислений. Для этого используем критерий серий [2, 16]. Составим вариационный ряд ~(к) из е(к), к = 1, ..., 10, и вычислим выборочное значение медианы по формуле:
med = 2 (~(к/ 2) + е(к/2 +1)) =0,0544.
Сравнивая каждое e(k) с медианой, ставим знак «+», если e(k) > med и знак «-», если e(k) < med . Образуется последовательность вида: — + +-------+ +. Серией называется последова-
тельность однотипных наблюдений, после которой следует наблюдения противоположного типа или же вообще нет никаких наблюдений. В вышеприведенном примере число серий S = 4.
Принятие гипотезы Н0 о независимости значений производится, если верно условие 5(102;1 - )< 5 5Х1^; ), где
критические значения определяются нормальным распределе-78
нием S и приводятся в [2], а - уровень значимости. Для нашего случая при уровне значимости а = 0,05 условие верно: 2 < 4 < 9. Последнее свидетельствует об отсутствии корреляции остатков в нулевой строке.
Зададим пятипроцентный интервал (в(к) - а; е(к) + а), в который должны попадать значения, стоящие в нулевой строке, в случае отсутствия структурно-параметрических изменений (таблица 1). Если при добавлении нового измерения хотя бы одно из значений, стоящих во 2-й строке выйдет за границы этого интервала, то с вероятностью рд = 0,95 можно утверждать, что модель объекта изменилась.
Таблица 1. Интервал нулевой строки к примеру 3
№ столбца 1 2 3 4 5
Верхняя граница 0,1027 0,0319 0,1612 0,1854 0,1061
Нижняя граница 0,0027 -0,0681 0,0612 0,0854 0,0061
№ столбца 6 7 8 9 10
Верхняя граница -0,0624 -0,0513 -0,0075 0,1624 0,1316
Нижняя граница -0,1624 -0,1513 -0,1075 0,0624 0,0316
Добавим 251-е значение случайного процесса х(к) без изменения его характеристик. Пересчитав структурную функцию по формуле (16), построим идентифицирующую матрицу:
№ столбца 1 2 3 4 5 6
(-1)строка 1 1 1 1 1 1
0 строка 0,2810 0,5371 0,7561 0,9732 1,1380 1,2856
1 строка -0,9116 -1,6911 -2,4639 -3,0504 -3,5754 -4,1428
2 строка 0,0566 -0,0117 0,1178 0,1284 0,0311 -0,0826
№ столбца 7 8 9 10 11 12
(-1)строка 1 1 1 1 1 1
0 строка 1,4449 1,6193 1,7800 1,9105 1,9914 2,0713
1 строка -4,7636 -5,3357 -5,8000 -6,0879 -6,3725
2 строка -0,0893 -0,0265 0,1219 0,0977
Значения, стоящие во 2-й строке, попадают в интервал (таблица 1), следовательно, нулевая строка сохраняется, и характеристики процесса не изменились.
На 252-м шаге изменяется корреляционная функция случайного процесса. Модель процесса принимает вид:
(18) х(к) = -0,9х(к -1) + п(к).
Графики процессов (17) и (18) приведены на рис. 4. Пересчитаем структурную функцию по рекуррентной формуле (16). Матрица на основе полученных значений примет вид:_______________
№ столбца 1 2 3 4 5 6
(-1)строка 1 1 1 1 1 1
0 строка 0,2889 0,5490 0,7675 0,9823 1,1464 1,3046
1 строка -0,9006 -1,6570 2,4006 -2,9689 -3,5162 -4,0987
2 строка 0,0608 -0,0084 0,1042 0,0647 -0,0346 -0,1002
№ столбца 7 8 9 10 11 12
(-1)строка 1 1 1 1 1 1
0 строка 1,4728 1,6414 1,8013 1,9190 1,9942 2,0774
1 строка -4,6823 -5,3179 -5,2360 -5,9039 -6,1918
2 строка -0,1013 -0,1313 -0,0298 0,0880
Значения в 5, 8, 9 столбцах 2-й строки выходят за границы интервала. Следовательно, изменение модели обнаружено без запаздывания. Добавим дополнительно 253-е наблюдение и посмотрим, какие изменения произойдут в матрице. Она принимает вид:
№ столбца 1 2 3 4 5 6
(-1)строка 1 1 1 1 1 1
0 строка 0,3024 0,5435 0,7659 0,9805 1,1408 1,2863
1 строка -0,7974 -1,5329 -2,2456 -2,7729 -3,2541 -3,8297
2 строка -0,1249 -0,2794 -0,2347 -0,3080 -0,5486 -0,7766
№ столбца 7 8 9 10 11 12
(-1)строка 1 1 1 1 1 1
0 строка 1,4604 1,6542 1,8079 1,9303 1,9899 2,0600
1 строка -4,4706 -4,9789 -5,3836 -5,5807 -5,8124
2 строка -0,7732 -0,7724 -0,6148 -0,7084
Вторая строка матрицы претерпевает значительные изменения и перестает быть нулевой. Ни одно из ее значений не попадает в интервал. Таким образом, изменение корреляционной функции обнаружено без запаздывания.
Пример 4. Имеется 50 значений временного ряда х(к), математическое ожидание которого равно нулю. В 51-й момент времени происходит изменение среднего значения ряда, которое становится равным 2 (рис. 5).
Рис. 5. Значения временного ряда х(к) к примеру 4
На основе имеющихся 50 измерений рассчитаем структурную функцию (4) и идентифицирующую матрицу (6):
1 2 3 4 5 6 7
1 1 1 1 1 1 1
0,5389 0,5114 0,4726 0,5891 0,5458 0,5619 0,5973
0,0511 0,1230 -0,0931 -0,0127 -0,0426 -0,1084
Первая строка нулевая, ее значения представляют собой независимые ошибки вычислений. Пятипроцентный интервал представлен в таблице 2.
Таблица 2. Интервал нулевой строки к примеру 4
№ столбца 1 2 2 4 5 6
Верхняя граница 0,1011 0,1730 -0,0431 0,0373 0,0074 -0,0584
Нижняя граница 0,0011 0,0730 -0,1431 -0,0627 -0,0926 -0,1584
После добавления 51 измерения структурная функция пересчитывается по формуле (16) и заполняется матрица:
1 2 3 4 5 6 7
1 1 1 1 1 1 1
0,6352 0,5168 0,5062 0,6432 0,5875 0,5813 0,6560
0,1863 0,2030 -0,0126 0,0752 0,0849 -0,0327
Элементы 1-5 столбцов выходят за границы интервала. Изменение характеристики случайного процесса обнаружено без запаздывания.
Для сравнения работоспособности критериев данная последовательность была также обработана с помощью инверсионного алгоритма [16, 24, 25], который на 56-м измерении обнаружил возрастающий тренд (запаздывание 6 значений).
Пример 5. 350 измерений случайного процесса х(к)
(19) х(к) = 0,5х(к -1) + 2п(к -1)
производились в нормальном режиме, где п(к) - белый шум с нулевым математическим ожиданием и дисперсией, равной
0,09. На основании значений структурной функции была рассчитана идентифицирующая матрица:
№ столбца 1 2 3 4
(-1)строка 1 1 1 1
0 строка 0,5076 0,7593 0,8840 0,9731
1 строка -0,4960 -0,7415 -0,9171 -0,9313
2 строка 0,0011 -0,1074 0,0395 0,0012
№ столбца 5 6 7 8
(-1)строка 1 1 1 1
0 строка 0,9803 0,9935 0,9546 0,9891
1 строка -0,9574 -0,8807 -0,9486
2 строка 0,1820 -0,0317
Вторая строка матрицы нулевая, ее элементы представляют собой независимые погрешности вычислений. Пятипроцентный интервал для нулевой строки представлен в таблице 3.
Таблица 3. Интервал нулевой строки к примеру 5
№ столбца 1 2 3 4 5 6
Верхняя граница 0,0511 -0,0574 0,0895 0,0512 0,2320 0,0183
Нижняя граница -0,0489 -0,1574 -0,0105 -0,0488 0,1320 -0,0817
На 351-ом измерении появляется нестационарность в форме линейного тренда и значения стохастического процесса становятся равными:
(20) х(к) = 0,5х(к -1) + 2п(к -1) + 0,04(к - 290).
Графики процессов (20) и (19) представлены на рис. 6. Матрица после 351-го наблюдения имеет вид:
№ столбца 1 2 3 4
(-1)строка 1 1 1 1
0 строка 0,5063 0,7573 0,8826 0,9711
1 строка -0,4957 -0,7432 -0,9180 -0,9309
2 строка -0,0035 -0,1088 0,0401 0,0001
№ столбца 5 6 7 8
(-1)строка 1 1 1 1
0 строка 0,9777 0,9909 0,9566 0,9909
1 строка -0,9572 -0,8893 -0,9570
2 строка 0,1632 -0,0413
к примеру 5
Все значения нулевой строки попали в интервал. Пересчитав структурную функцию и идентифицирующую матрицу с поступлением 352-го наблюдения, получим:___________________
№ столбца 1 2 3 4
(-1) строка 1 1 1 1
0 строка 0,5099 0,7624 0,8869 0,9795
1 строка -0,4951 -0,7394 -0,9210 -0,9317
2 строка 0,0016 -0,1208 0,0391 0,0246
№ столбца 5 6 7 8
(-1) строка 1 1 1 1
0 строка 0,9850 0,9913 0,9616 1,0078
1 строка -0,9442 -0,8858 -0,9764
2 строка 0,1550 -0,0864
Значение, стоящее в 6-ом столбце матрицы, не попало в интервал. Добавим 353-е наблюдение:
№ столбца 1 2 3 4
(-1) строка 1 1 1 1
0 строка 0,5096 0,7713 0,8984 0,9904
1 строка -0,5133 -0,7633 -0,9435 -0,9657
2 строка 0,0265 -0,0745 0,0623 0,2264
№ столбца 5 6 7 8
(-1) строка 1 1 1 1
0 строка 1,0018 1,0068 0,9672 1,0200
1 строка -0,9755 -0,8979 -1,0014
2 строка -0,0527 -0,2060
Значения в 4, 5, 6 столбцах не попали в интервал. Изменение модели обнаружено после 352-го наблюдения (запаздывание -1 наблюдение).
Пример 6. Имеется 250 измерений стационарного случайного процесса, имеющего модель
(21) х(к) = 0,85х(к -1) + п(к -1) + 0,4п(к),
где п(к) - белый шум с нулевым математическим ожиданием и дисперсией, равной 0,04. На основе имеющихся наблюдений была получена идентифицирующая матрица для структурной функции:
1 2 3 4 5 6
1 1 1 1 1 1
0,0512 0,1287 0,1943 0,2514 0,2957 0,3308
-1,5163 -2,7969 -3,9141 -4,7804 -5,4665
0,6717 1,2156 1,7614 2,1753
0,0349 -0,0409 -0,0857
Третья строка матрицы нулевая. Пятипроцентный интервал для нулевой строки представлен в таблице 4.
Таблица 4. Интервал нулевой строки к примеру 6
№ 1 2 3
Верхняя граница 0,0849 0,0091 -0,0357
Нижняя граница -0,0151 -0,0909 -0,1357
На 251-м наблюдении происходит параметрическое изменение модели случайного процесса. Модель процесса принимает вид:
(22) х(к) = 0,95х(к -1) + 3п(к -1) + 0,8п(к).
График процессов (21) и (22) представлен на рис. 7.
х(к)
изменение модели процесса
Рис. 7. График процессов (21) и (22) к примеру 6
После пересчета структурной функции на 251 шаге идентифицирующая матрица имеет вид:
1 2 3 4 5 6
1 1 1 1 1 1
0,0510 0,1283 0,1935 0,2504 0,2945 0,3300
-1,5167 -2,7978 -3,9133 -4,7794 -5,4745
0,6720 1,2176 1,7621 2,1621
0,0328 -0,0419 -0,0777
Все значения 3-ей строки попадают в интервал. Пересчитаем матрицу после добавления 252 наблюдения:
1 2 3 4 5 6
1 1 1 1 1 1
0,0508 0,1278 0,1930 0,2496 0,2934 0,3286
-1,5145 -2,7973 -3,9130 -4,7750 -5,4679
0,6675 1,2136 1,7602 2,1645
0,0287 -0,0534 -0,0901
Нулевая строка по-прежнему третья. Добавим следующее 253 наблюдение:
1 2 3 4 5 6
1 1 1 1 1 1
0,0511 0,1275 0,1925 0,2487 0,2923 0,3275
-1,4948 -2,7670 -3,8673 -4,7202 -5,4087
0,6437 1,1798 1,7095 2,1018
0,0183 -0,0685 -0,0107
Все значения нулевой строки попали в интервал. Пересчитав структурную функцию и идентифицирующую матрицу на 254 шаге, получим:
1 2 3 4 5 6
1 1 1 1 1 1
0,0509 0,1276 0,1940 0,2494 0,2930 0,3274
-1,5067 -2,8103 -3,8979 -4,7548 -5,4296
0,6417 1,2235 1,7423 2,1514
-0,0416 -0,1284 -0,1972
Все значения 3-ей строки выходят за границы интервала. Третья строка перестает быть нулевой. Параметрическое изменение модели процесса зафиксировано на 254-м наблюдении. Запаздывание составило 3 измерения.
7. Выводы
Предложенный алгоритм обнаружения структурнопараметрических изменений в моделях динамических объектов, процессов и систем, основанный на использовании модели структурной функции случайного процесса в форме непрерывной дроби, позволяет улавливать изменения либо в момент их появления, либо с небольшим запаздыванием (1-3 значения), что подтверждено рядом модельных примеров. Метод является непараметрическим, не требует знания априорной информации и ввиду его формализованности может применяться для решения различных практических задач контроля над функционированием объектов и систем в режиме реального времени. В отличие от других существующих методов, предложенный метод не требует настройки на конкретный тип изменений, стационарности данных измерений, знания законов распределений данных. Кроме того, на сложность и структуру метода не оказывает влияния неизвестность момента (или его распределения) появления изменения. Метод позволяет эффективно отслеживать изменения всех характеристик случайного процесса: среднего, дисперсии, корреляционной функции и появления нестационар-ности. Развитая теория структурно-параметрической идентифи-
кации [16, 19-21, 27] на подготовительном этапе позволяет осуществлять эффективную проверку данных на стационарность и однородность, получать модели объектов и систем без их перебора. Помимо этого при проектировании системы дискретного контроля за функционированием непрерывного динамического объекта необходимо, чтобы ее неотъемлемой функцией была возможность изменения шага дискретизации. Принцип вариации шага дискретизации [16, 19-21, 27] является достаточным условием восстановления свойств и характеристик непрерывного объекта и процесса. Неизменность непрерывной модели объекта в случае выполнении принципа вариации шага дискретизации при обработке временных рядов может служить критерием достоверных оценок их свойств.
Литература
1. АНДЕРСОН Т. Статистический анализ временных рядов. -М.: Мир, 1976.
2. БЕНДАТ ДЖ., ПИРСОЛ А. Прикладной анализ случайных процессов. - М.: Мир, 1989.
3. БОКС ДЖ., ДЖЕНКИНС Г. Анализ временных рядов. Прогноз и управление. - М.: Мир, 1974. - Вып. 1., 1974. - Вып. 2.
4. БОРОДКИН ЛИ., МОТТЛЬ В.В. Алгоритм обнаружения моментов изменения параметров уравнения случайного процесса // Автоматика и телемеханика. - 1976. - №6. - С. 23-31.
5. БРАВЕРМАН Э.М., МУЧНИК И.Б. Структурные методы обработки эмпирических данных. - М.: Наука, 1983.
6. БРИЛЛИНДЖЕР Д.Р. Временные ряды: обработка данных и теория. - М.: Мир, 1980.
7. БРОДСКИЙ Б.Е., ДАРХОВСКИЙ Б.С. Непараметрический метод обнаружения моментов переключения двух случайных последовательностей // Автоматика и телемеханика. -1989. - №10. - С. 66-74.
8. ВАЛЬД А. Последовательный анализ. - М.: Физматгиз, 1960.
9. ВОРОБЕИЧИКОВ С.Э., КОНЕВ В.В. Об обнаружении разладок в динамических системах // Автоматика и телемеханика. - 1990.
- №3. - С. 56-58.
10. ДАРХОВСКИЙ Б.С. Непараметрический метод для апостериорного обнаружения момента «разладки» последовательности независимых случайных величин // Теория вероятностей и ее применения. - 1976. - Том 21. - С. 180-184.
11. ДУДНИКОВ Е.Г., БАЛАКИРЕВ В.С., КРИВСУНОВ ВС.
Построение математических моделей химико-
технологических объектов. - Ленинград: Химия, 1970.
12. ЖИГЛЯВСКИЙ А.А., КРАСНОВСКИЙ А.Е. Обнаружение разладки случайных процессов в задачах радиотехники. -Ленинград: Изд-во Ленинградского университета, 1988.
13. КАРТАШОВ В.Я. Анализ и исследование аппроксимационных свойств непрерывных дробей при решении задачи структурнопараметрической идентификации динамических объектов // Препринт №22. - Барнаул: Изд-во Алтайского госуниверситета, 1996.
14. КАРТАШОВА Л.В., КАРТАШОВ В.Я. Построение причинно-следственных моделей социально-экономических процессов: монография. - Томск: Издательство Томского государственного педагогического университета, 2008.
15. КАРТАШОВ В.Я. Непрерывные дроби (определения и свойства). - Кемерово: Изд-во Кемеровского госуниверси-тета, 1999.
16. КАРТАШОВ В.Я., НОВОСЕЛЬЦЕВА М.А. Идентификация стохастических объектов: учебное пособие. - Томск: Изд-во Томского государственного педагогического университета, 2008.
17. КАРТАШОВ В.Я., НОВОСЕЛЬЦЕВА М.А. Обнаружение изменений внутренних свойств динамических объектов // Материалы XII Международной конференции по мягким вычислениям и измерениям, Санкт-Петербург, 2009. -
С. 118-121.
18. КАРТАШОВ В.Я., НОВОСЕЛЬЦЕВА М.А. Определение риска по структурным особенностям временных рядов // Материалы X Международной конференции по мягким вычислениям и измерениям, Санкт-Петербург, 2007. - С. 73-76.
19. КАРТАШОВ В.Я., НОВОСЕЛЬЦЕВА М.А. Способ идентификации линейного объекта // Патент РФ №2189622. - 2002. -
Бюл. №26.
20. КАРТАШОВ В.Я., НОВОСЕЛЬЦЕВА М.А. Способ идентификации линейного объекта // Патент РФ №2146063. - 2000. -Бюл. №6.
21. КАРТАШОВ В.Я., НОВОСЕЛЬЦЕВА М. А. Структурнопараметрическая идентификация стохастических объектов с использованием непрерывных дробей // Управление большими системами. - 2008. - Вып. 21. - С. 27-48.
22. КОНЕВ В.В. Последовательные оценки параметров стохастических динамических систем. - Томск: изд-во Томского ун-та, 1985.
23. МАКС Ж. Метода и техника обработки сигналов при физических измерениях. - М.: Мир, 1983. - Том 1.; Том 2.
24. МикроЭВМ в информационно-измерительных системах. - М.: Машиностроение, 1987.
25. МИРОНОВ И.И., ОСИПОВ С.Н. Многоконтурные системы обработки информации и активного управления. - М.: Энерго-атомиздат, 1997.
26. НИКИФОРОВ И.В. Последовательное обнаружение изменения свойств временных рядов. - М.: Наука, 1983.
27. НОВОСЕЛЬЦЕВА М.А. Идентификация моделей совместных случайных процессов для систем контроля горной техники: дис. канд. тех. наук. - Кемерово, 2001.
28. НОВОСЕЛЬЦЕВА М.А. Оценка риска в стохастических системах // Материалы VII Всероссийской научнопрактической конференции А8’2009 «Системы автоматизации в образовании, науке и производстве». - Новокузнецк: изд-во СибГИУ, 2009. - С.441-446.
29. Обнаружение изменения свойств сигналов и динамических систем / Под ред. М. Бассвиль, А. Банвениста. - М.: Мир, 1989.
30. РОМАНЕНКО А.Ф. СЕРГЕЕВ Г.А. Вопросы прикладного анализа случайных процессов. - М.: Советское радио, 1968.
31. СОЛОДОВНИКОВ ВВ., ПЛОТНИКОВ В.Н., ЯКОВ-
ЛЕВ А.В. Основы теории и элементы систем автоматического регулирования: учебное пособие для вузов. - М.: Машиностроение, 1985.
32. ШИРЯЕВ А.Н. О марковских достаточных статистиках в неаддитивных байесовских задачах последовательного анализа // Теория вероятностей и ее применения. - 1964. -№IX, 4. - С. 670-686.
33. ШИРЯЕВ А.Н. Статистический последовательный анализ.
- М.: Наука, 1976.
34. PAGE E.S. Continuous Inspecting Schemes // Biometrica. -1954. - Vol. 41. - P. 100-114.
35. WILLSKY A.S. A Survey of Design Methods for Failure Detection in Dynamic Systems // Automatica. - 1976. - Vol. 12. -P. 601-611.
DETECTION OF STRUCTURALLY-PARAMETRICAL CHANGES IN STOCHASTIC SYSTEMS IN REAL TIME BY ALGORITHMS OF CONTINUED FRACTIONS AND STRUCTURE ANALYSIS
Vladimir Kartashov, head of the Department of Automation of Researches and Technical Cybernetics of the Kemerovo State University, Doctor of Engineering, professor (kartash@kemsu.ru). Marina Novoseltseva, associate professor of Department of Automation of Researches and Technical Cybernetics of the Kemerovo State University, Candidate of Engineering, (aanov@pochta.ru).
Abstract: An algorithm is developed for real-time detection of
structural and parametric changes in stochastic processes and systems. The algorithm is based on the combined use of structural analysis and theory of continued fractions.
Keywords: structural function, the continued fraction, identifying matrix.
Статья представлена к публикации членом редакционной коллегии Я. И. Квинто
