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

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

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


Pages:     | 1 || 3 | 4 |

«Министерство образования и науки РФ Научно-исследовательский Иркутский государственный технический университет ...»

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

Yk 1/ k – оценка вектора Y на момент k-1, вычисленная с учетом k имеющихся наблюдений.

Соотношения (2.33) носят название дискретного фильтра Калмана.

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

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

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

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

При построении фильтра учитываются особенности рассматриваемой динамической системы. Могут использоваться различные классы математических моделей: полиномиальные, экспоненциальные и др. Для рассматриваемого случая групповых эталонов времени и частоты наиболее естественным способом является применение уравнений авторегрессии – проинтегрированного скользящего среднего (АРПСС), так как этот класс моделей предназначен для обработки временных рядов [18]. Модели АРПСС относятся к динамическим стохастическим моделям [30].

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

Стандартные преобразования (2.19) применимы к системам с одним входом и одним выходом. Обобщение стандартного преобразования для случая многомерных систем выполнено D. Graup [19]. При этом в уравнениях (2.18) скаляры yk и zk становятся векторными величинами Yk, Zk, векторы g, h* блочными матрицами и H*, а матрица * – блочной матрицей ** * 1 * 0, ** (2.34) 0 0 * 0 0 n где каждое * имеет вид * из уравнения (2.18). Уравнения (2.18) i примут в таком случае вид:

Yk **Yk 1 Ak (2.35) Z k H *Yk Элементы Hij блочной матрицы H* отражают связи между различными входами и выходами системы. На главной диагонали этой матрицы расположены элементы вектора H H, H ** T, H n, 1 где каждое Hi соответствует вектору H* в уравнениях (2.18).

В групповых эталонах, как и в иных системах с независимыми друг от друга элементами, взаимодействие между различными входами и выходами происходит только на уровне оператора H*. По существу, речь идёт об одномерных системах, описываемых уравнениями динамики (2.18) и объединёнными в многомерную измерительную систему оператором H*. Из-за независимости шумов различных стандартов частоты в эталонных комплексах эталонов времени и частоты (такая независимость обеспечивается техническими мерами – экранированием, развязкой цепей питания и др.) и вообще, при независимых элементах любых других многомерных систем, матрица тоже диагональная: элементы матрицы, расположенные на главной диагонали являются составляющими расширенного вектора g* g g, g,, g, *T 1 2 n где каждое gi определёно формулой (2.19).

Оператор H* представлен блочной матрицей 0 0 0 0 0 0 0 E E E1 H 1, 0 0 E1 E1 Первая строка матрицы H* представлена нулевыми элементами, поскольку отражает “фиктивное” измерение y1-y1. Матрица Q также блочная Q1 0 E Q Q 1, Qn 0 где Qi определены соотношениями Q a I (I – единичная матрица), поскольку предполагается, что шумы не коррелированны.

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

Субоптимальные фильтры, применяемые в системах автоматического управления, основываются на методах понижения размерности вектора состояния и методах расщепления (декомпозиции) фильтра [20,73].

При понижении размерности вектора состояния системы получение оценок разбивается на две подсистемы. Первая подсистема – динамическая (дополнительный канал наблюдения), описывается конечно-разностным уравнением t 1 f X t D f Z t, (2.36) где f и D f - матрицы размером N M N M и N M M, N – размерность вектора состояния Y в фильтре Калмана.

M – размерность вектора измерений Z.

Вторая подсистема завершает построение оценки вектора состояния.

Выходной сигнал представляется в виде [20] Y t X t T Z t (2.37) Метод декомпозиции вектора состояния основывается на представлении его в виде суммы k Y k F j yi, (2.38) j где y j D j Y i jk, F j и D j - матрицы, удовлетворяющие условию k F D I, (2.39) j j j где I – единичная матрица.

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

Используем метод понижения размерности для получения оценки состояния опорного элемента Y1. Вектор состояния системы найдётся из выражения (2.37). Однако, в отличие от рассмотренного в [20] метода, подвектор Y1 найдём в соответствии с принципом декомпозиции, используя соотношение (2.38). При этом уравнение для оценки подвектора yi примет вид y *j t y j t k j Gi Z HY, (2.40) где y j t - прогноз подвектора y j, kj – коэффициент усиления невязки, Gj – матрица, определяющая, как вектор Z используется для коррекции оценки yj.

Структура такого субоптимального фильтра представлена на рис. 2.3.

G1 K1 F T * … * Z Gj Kj Fj P … B Gk Kk Fk Eз H Рис. 2.3 – Структурная схема субоптимального фильтра Операторы,,,, F, K соответствуют рассмотренным выше матрицам. E3 – оператор, выделяющий подвектор состояния опорного элемента из вектора Y, B – оператор задержки на один такт.

Декомпозицию вектора состояния рекомендуется производить таким образом, чтобы подвекторы не были связаны между собой ни через уравнения объекта, ни через помехи [20]. В соответствии с этой рекомендацией выполним декомпозицию так, чтобы подвекторы были равными yi. При одинаковой размерности всех подвекторов матрицы Dj имеют вид D j 0 0 I 0.

p p Подматрица I размерности расположена в j-м столбце блочной матрицы Dj (p – размерность подвекторов). Матрицы Gj при таком разбиении будут равны матрицам Dj. Для выполнения требования восстановимости подвектора состояния Y1 в рассматриваемом фильтре условие (2.39) должно быть заменено условием (2.41) n D j I 0 0 n k.

F (2.41) j j Из этого условия следует, что матрицы Fj должны быть заменены операторами, осуществляющими перенос содержимого j-го блока в первый с инвертированием знака и умножением на скалярную величину i, такую n что I.

j j Коэффициент усиления субоптимального фильтра, построенного по принципу декомпозиции, определяется соотношениями [20] n K t F j K j t G j, (2.42) j K j t Pj t H j H j Pj t H T R j t, (2.43) j Pj t j Pj t 1 Tj j Q j T, (2.44) j где Pj t - экстраполированная ковариационная матрица подвектора Yj, Pj t - ковариационная матрица, вычисленная с помощью измерения, выполненного в момент t.

Матрицы Pj t и Pj t связаны соотношением [9] Pj t I K t H Pj t (2.45) Матрицы Qj, Hj, Rj, j в (2.43), (2.44) находятся из соответствующих матриц полного фильтра [20] Q j D j Q j QT ;

H j G j H D ;

R j G j R RT ;

j D j D j j j j D Псевдообратная матрица находится с помощью скелетного j разложения матрицы D [10] D j B C ;

D C T CC T BB T B T 1 (2.46) j Нетрудно убедиться, что D D T (матрица B – единичная, размерности j j p p, матрица C равна матрице D).

Если элементы системы и шумы измерений не зависят друг от друга (в групповых эталонах это справедливо), то матрицы Pj, Qj, Rj, j являются j-ми диагональными блоками соответствующих матриц полного фильтра. Матрица Hj имеет размерность p p. Её вид определяется матрицей H из выражений (2.35) и вектором измерений в стандартном каноническом преобразовании [20] 0 0 0 Hj 0 Подвектор Y1*, являющийся оценкой состояния опорного элемента, будет определяться выражениями n Y1* j Y j Z j (2.47) j n j 2,, p ;

1 I j (2.48) j Fj K j j Y j Y1 G1 H Y ;

(2.49) Y1 E3 Y ;

(2.50) (2.51) Z j Gj Z Y B Y* (2.52) Выходной сигнал фильтра Y * формируется в соответствии с формулой (2.37) (метод понижения размерности). Выбрав матрицы P и T так, что T I I I, T получим T Y * Y1* Y1* Z 2 Y1* Z 3 Y1* Z n (2.53) Выражения (2.47), (2.53) аналогичны соответствующим выражениям для оценок относительных отклонений частот генераторов, полученных на основе применения моделей АРСС (авторегрессии – скользящего среднего) [15].

Следовательно, алгоритмы обработки данных, опирающиеся на использование моделей АРСС, при определённых условиях являются алгоритмами субоптимальной фильтрации. Такими условиями являются, очевидно, равенства весов gi в алгоритме фильтрации соответствующим элементам матричных коэффициентов j в (2.48)-(2.50). Коэффициенты усиления kj определены выражениями (2.42)-(2.44). При принятых выше представлениях матриц, входящих в эти выражения, получаем j P11 0 P11j r11j j P21 j 0 P11 r11j, Kj j Ppj 0 P11j r11j где P j и r11 - левые верхние элементы ковариационных матриц Pj и R j.

j В соответствии с условием восстановимости вектора состояния, матрица i равна матрице Kj, умноженной на скалярную величину j (операция перемещения блока Kj размерности p p в первую клетку блочной матрицы размерности p p не меняет численных значений элементов блока). Таким j образом, элемент 11, соответствующий весовому коэффициенту g i равен Pj j j 11 j j. (2.54) P r n Выполняя условие нормировки 1 получим j j 1 P11 r j j j. (2.55) n 1 P j j r 11 j В результате получаем P11 P11j r j j j j 2,, n (2.56) 1 P11j r11j n 1 j 11 1 11. (2.57) j Так как элементы P j экстраполированной ковариационной матрицы Pj являются дисперсиями прогноза состояния j-го элемента системы, то при отсутствии шумов измерений ( Rj 0 ввиду их малости) весовые j коэффициенты 11, найденные из условия субоптимальной фильтрации, совпадают с весами gi, используемыми в алгоритмах оценивания состояния, опирающихся на динамические стохастические модели (АРСС) [10,15,21].

Следует отметить, что на связь фильтров Калмана с моделями временных рядов указывал академик В.С. Пугачёв [18].

Таким образом, использование фильтра Калмана является достаточно формальным с математической точки зрения способом построения процедуры динамической обработки измерительной информации в эталонах времени и частоты, основанной на использовании прогнозов относительных отклонений частот, вычисленных посредством моделей АРСС. Показано, что предлагаемая процедура оценивания является процедурой субоптимальной Калмановской фильтрации.

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

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

Модели АРПСС относятся к динамическим стохастическим моделям и могут описывать стационарные и стационарно-разностные процессы. В случае стационарных процессов предполагается, что процесс остается в равновесии относительно постоянного среднего уровня. Нестационарные процессы (и их модели) не имеют естественного среднего значения. Для описания стационарных процессов используются модели авторегрессии – скользящего среднего (АРСС), нестационарные же процессы описываются моделями авторегрессии — проинтегрированного скользящего среднего (АРПСС) [30].

Существует большое количество различных расширений моделей АРПСС [82].

Временные ряды, в которых последовательные значения существенно зависимы, можно рассматривать как генерируемые последовательностью независимых импульсов at. Эти импульсы – реализации случайных некоррелированных величин с фиксированным распределением, которое обычно предполагается нормальным с нулевым средним значением и дисперсией a2 (“белый шум”). В модели авторегрессии текущее значение процесса выражается как конечная линейная совокупность предыдущих значений процесса и импульса at [14].

Обозначим значения процесса в равноотстоящие моменты времени t, t 1, t 2... как yt, yt 1, yt 2.... Обозначим как ~t, ~t 1, ~t 2,... отклонения от yy y среднего значения, так что ~t yt. Тогда процесс авторегрессии (АР) y порядка p будет описываться выражением.

~ ~ ~... ~ a (2.58) yt 1 y t 1 2 yt 2 p yt p t Если ~t линейно зависит от конечного числа q предыдущих a, – такой y процесс называется процессом скользящего среднего (СС) порядка q :

~ a a a... a (2.59) yt t 1 t 1 2 t 2 q t q Для достижения большей гибкости в подгонке моделей к наблюдаемым временным рядам во многих случаяхв используют комбинированную модель авторегрессии – скользящего среднего (АРСС):

~ ~... ~ a a... a (2.60) yt 1 y t 1 p yt p t 1 t 1 q t q При порядке разности d=0 модель описывает стационарный процесс.

Модели с порядком разности d, отличным от 0 (модели авторегрессии – проинтегрированного скользящего среднего - АРПСС) применяются для описания нестационарных процессов, содержащих детерминированные составляющие. Мнения специалистов по поводу использования разностных авторегрессионных моделей расходятся [23]. Это связано с рядом сложностей в отношении интерпретации получаемых результатов при порядке разности d0.

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

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

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

Классический подход, предлагаемый Боксом и Дженкинсом [14], опирается на изучение автокорелляционной и частной автокорелляционной функций исходного временного ряда, в результате которого определяется характер и порядок процесса, а, следовательно, количество включаемых в модель членов АР и СС. После этого тем или иным способом осуществляется подгонка параметров модели.

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

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

n N I z j, i y j,1 (t ) y j, i 1 (t ) (2.61) j 1 i Данный функционал зависит от величины ошибки прогноза моделей АРПСС, а, следовательно, и от качества моделей.

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

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

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

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

Из практики применения моделей АРПСС известно, что полезные значения количества членов АР и СС (p и q), при которых достигается адекватное описание наблюдаемых временных рядов лежат в интервале [0;

3].

В частности, для водородных стандартов частоты, эксплуатирующихся в ГСВЧ, успешно применялись модели АРПСС, порядки которых (p, d, q) равнялись (2, 1, 2) [15]. Вместе с тем, имеются предположения о том, что для более точного моделирования процессов, возникающих в выборке при выполнении наблюдений на малых интервалах может оказаться оправданным применение моделей более высоких порядков. Так, в [48] для моделирования каждого из водородных генераторов применялась линейная комбинация из 4 марковских процессов (моделей авторегрессии порядка 1) и сделано предположение о возможности замены их одной моделью авторегрессии порядка 4.

Полагая порядки АР и СС минимальными для всех элементов группового эталона (т.е. p=0, q=0), можно последовательно повышать их, решая при этом оптимизационную задачу. Очевидно, что, как и в регрессионном анализе, увеличение p и q не должно вести к росту функционала I, где обобщённый вектор параметров моделей [21]. Однако, как нетрудно заметить, pmax qmax 3, даже для при ограничении максимального числа членов группового эталона, включающего всего два хранителя, оптимизационную задачу необходимо решать 256 раз. При увеличении числа мер в групповом эталоне трудоёмкость подобного алгоритма стремительно возрастает. Поэтому, несмотря на теоретическую значимость, практическая ценность такого алгоритма невелика, и имеет смысл исследовать только упрощённый алгоритм с одновременным наращиванием порядков для всех моделей. Множество возможных моделей в таком случае позволяет перебрать все практически целесообразные варианты, вычислить значения функционала (проведя параметрическую оптимизацию), а затем ранжировать эти модели по какому либо признаку, и определить “наилучшую” (в заданном смысле). Если порядок разности в модели d принимается равным нулю (отсутствует тренд), потребуется перебрать 24 модели (модель с нулевым количеством членов рассматривать нецелесообразно), что хотя и требует достаточно продолжительных вычислений, но может быть выполнено с привлечением разумных вычислительных ресурсов.

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

Как упоминалось выше, исходные ряды могут наилучшим образом описываться моделями АРПСС различной структуры (т.е. процессы, порождаемые каждым из хранителей, могут иметь различные типы и структуры), следовательно, для наиболее точного результата следовало бы перебирать модели для каждого ряда, что приводит к чрезмерному росту вычислительной ёмкости задачи. При выборе таким способом модели в случае разностных рядов, полученных из временных рядов, описываемых моделями различной структуры, поведение функционала I может иметь более сложный характер. Таким образом, функционал I не может быть непосредственно использован для выбора структуры модели.

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

I I2, где: (2.62) N (n k ) N (n ( p q)) - число степеней свободы, N – число исходных рядов, n – число членов в каждом ряду, p – число членов АР, включенных в модель, q – число членов СС, включенных в модель.

Выражение I 2 представляет собой функционал I, разделенный на количество степеней свободы, которое учитывает число членов ряда и оцениваемых параметров. Таким образом, знаменатель имеет тенденцию к уменьшению по мере роста сложности модели, а числитель представляет собой тот же функционал I. По существу, вводится своеобразный “штраф” за каждый дополнительный член АР или СС, включенный в модель.

Результаты моделирования, приведённые в соответствующем разделе, показали, что хотя поведение критерия существенным образом отличается от поведения функционала I, в результате чего в числе “наилучших” оказываются модели с менее сложной структурой, однако представляющиеся наиболее эффективными модели в их число всё равно не попали. Такой результат оказался обусловлен недостаточно большим “штрафом” за увеличение сложности модели: сумма ( p q ) нарастает недостаточно динамично и мало влияет на знаменатель при большом количестве наблюдений n. В качестве выхода из ситуации было предложено увеличение такого “штрафа”, путём введения дополнительного коэффициента. Домножением суммы количества членов авторегрессии и скользящего среднего на некоторый коэффициент k было получено выражение:

I, (2.63) I N (n k ( p q)) где k1 – добавочный коэффициент.

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

Таким образом, данный критерий заведомо не является универсальным.

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

При этом высокая трудоёмкость вычислений не является препятствием.

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

Следует отметить, что узким местом описанного способа структурной идентификации является процедура “подгонки” модели. Будучи реализованной методом градиентного спуска с дроблением шага, она требует вычисления частных производных функционала I по каждому из коэффициентов, что, при численном дифференцировании, в свою очередь, требует пересчёта функционала для каждого приращения параметра. Таким образом, на каждом шаге спуска для модели структуры (p,d,q) требуется вычислить функционал I некоторое количество раз K, при этом K может быть определено, как:

K=1+N·(p+q)+1= 2+N·(p+q), (2.64) где N – число исходных рядов, p – число членов авторегрессии, включенных в модель, q – число членов скользящего среднего, включенных в модель.

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

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

В качестве альтернативного подхода к решению задачи построения моделей АРПСС при неизвестной структуре может быть предложено использование принципа “дополнительной суммы квадратов”. Вопрос о том, стоит ли включать в регрессионную модель те или иные члены может решаться анализом дополнительной части суммы квадратов, порождённой регрессией, которая связана с включением в модель рассматриваемых членов. Средний квадрат, который получается из этой дополнительной суммы может быть далее сопоставлен с оценкой s2 параметра 2, чтобы выяснить, является ли значимым различие между ними. Если средний квадрат значимо превышает оценку s2, то такие члены следует включить в модель. В противном случае их следует рассматривать как излишние, и не включать в модель [25].

Применительно к рассматриваемой задаче, использование принципа “дополнительной суммы квадратов” позволяет, с одной стороны, отказаться от жёсткого ограничения на порядок модели, а с другой – избежать перебора всего множества возможных моделей с последующей их подгонкой и вычислением функционала I: наращивание порядка модели прекратится после признания дополнительной суммы квадратов, внесённой добавочными членами вектора параметров, незначимой. Проблемой при использовании такого метода, однако, может явиться то, что из-за действия случайных факторов предыдущая модель может оказаться значимо хуже последующих, которые будут отсечены ввиду незначительности дополнительной суммы квадратов предыдущего шага.

Ещё один альтернативный подход заключается в построении предварительных оценок относительных отклонений частоты элементов группового эталона с последующим их уточнением на основе соотношения (2.20). При этом для получения начальных оценок находится решение уравнений (2.18) с помощью псевдообратной матрицы A. Известно, что такое решение доставляет минимум сумме квадратов отклонения оценок составляющих вектора состояния Y (t ) при минимальной его норме.

Псевдообратная матрица A может быть найдена с помощью скелетного разложения матрицы [10] b1r b c11 c12 c1n b b2 r, A B C 21 (2.65) c r 1 c r 2 c 2 n bmr bm где ранги сомножителей B и C равны рангу произведения, rB rC r.

Тогда 1 A C B C C C B B B. (2.66) В качестве столбцов матрицы B берутся первые r столбцов матрицы A ( r k 1 ).

1 1 B 1 0 Элементы матрицы C найдутся из решения матричного уравнения (2.65).

Например, при размерности матрицы A, равной (2x3) получаем матрицы 1 1 -1 1 0 1 1 C 0 1 1 и A 3 2 1.

B (2.67) 1 0 1 С учётом того что z1 y1 y2 (2.68) z 2 y1 y имеем y1 y2 y1 y Y A Z 2 y1 y2 y1 y3.

(2.69) 3 y1 y2 2 y1 y Таким образом, оценка состояния опорного элемента находится как:

2 y1 y 2 y, (2.70) y или, с учётом (2.68) y1 zi.

3 i Для произвольной размерности матриц выражение имеет вид:

1 n zi (2.71) y N i Оценки состояния других элементов могут быть найдены либо непосредственно из соотношения (2.68), либо с помощью простого пересчёта результатов измерений для данного момента времени t.

yi(t) y1(t) zi(t) (2.72) Полученные результаты легко обобщаются для произвольной размерности матрицы A.

Временные ряды, построенные по формулам (2.71), (2.72) можно использовать для идентификации структуры моделей АРСС, т.е. для определения параметров p и q. Для этого требуется построить для каждого полученного таким образом временного ряда “грубых” оценок автокорреляционные (АКФ) и частные автокорреляционные (ЧАКФ) функции, а затем использовать методику Бокса-Дженкинса, основанную на сопоставлении выборочных АКФ и ЧАКФ их теоретическим значениям при конкретных p и q. При этом необходимо помнить, что исследуемые временные ряды (ряды оценок yi(t) ) содержат погрешности, которые неизбежно повлекут за собой погрешности при оценивании АКФ и ЧАКФ.

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

N N 1 y1 i.

y i y1 y1 (2.73) N N Если рассматриваемые линейные процессы некоррелированы, АКФ суммы процессов равна сумме автокорреляционных функций [26]. Первый член в выражении (2.73) определяет смещение оценок y1(t), что не вносит изменений в АКФ временного ряда. Второй член с ростом N при независимых временных рядах стремится к нулю и, следовательно, также не вносит существенной погрешности в определение порядков p и q. Однако, при малых N, а в реальности приходится иметь дело именно с такой ситуацией, влияние этой составляющей погрешности может существенно исказить оценку порядков p и q. Впрочем, сама методика Бокса-Дженкинса, основанная на понятии “близости” теоретических и выборочных автокорреляционных и частных автокорреляционных функций, не даёт точного ответа на вопрос о структуре модели, а предполагает лишь нахождение приближённых значений порядков АР и СС. Поэтому зачастую оценивают порядок p и q приближённо, исследуя несколько альтернативных структур моделей.

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

1. По имеющимся результатам взаимных измерений элементов группового эталона находятся “грубые” (предварительные) оценки относительных отклонений их частот от приписанных значений по формулам (2.71), (2.72) для каждого момента t t 1,n.

2. По найденным временным рядам строятся их динамические модели (модели АРСС). Структурная идентификация выполняется с использованием в качестве “исходных” рядов “грубых” оценок. Также на этом этапе могут быть с минимальными трудозатратами получены начальные приближения для параметров моделей АРСС.

3. С использованием полученных структур моделей (и, опционально, полученных начальных приближений оценок параметров моделей) путём минимизации функционала (2.61) строятся уточнённые оценки на основе соотношения (2.20).

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

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

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

В общем случае временные ряды могут быть представлены в виде [27] Y t A f тр t Б t В t t, (2.74) 1, если факторы типа c участвуют в формировании Y t где (c) 0 в противном случае c=A, Б или В;

f тр t - функция тренда (долговременная составляющая временного ряда – полиномиальный тренд);

t - динамическая составляющая (динамический тренд);

t - сезонная составляющая;

t - случайная составляющая.

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

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

Таким образом, учёт детерминированных функций f тр t в построенных моделях принципиально может быть выполнен двумя способами:

1. Непосредственной оценкой параметров линейных функций f i t b0 b1i t i (или более сложных) по результатам измерений и последующим исключением влияния этих функций из результатов измерений zi t ;

2. Использованием разностных моделей y t y t y t-1, включающих в правой части постоянный член 0 [22,27] – в этом случае речь идёт о моделях авторегрессии-проинтегрированного скользящего среднего (АРПСС).

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

При таком подходе на первом этапе построения моделей временных рядов требуется оценить параметры линейных функций. Затем необходимо устранить влияние этих функций из рядов измерений zi t и таким образом свести задачу построения моделей временных рядов к описанной ранее процедуре.

Поскольку при формировании “автономной шкалы времени”, т.е. шкалы времени, построенной только по результатам внутренних сличений отдельного эталона, в обязательном порядке производится начальная “привязка” шкалы эталона к шкале эталонов более высокого уровня (например, шкалы вторичных эталонов “привязываются” к шкале Государственного эталона), то можно по результатам внешних сличений установить начальные отклонения частоты, характеризующиеся членами b0i и внести соответствующие поправки в ряды измерений zi t. Иными словами, при оценивании линейных трендов можно полагать все коэффициенты b0i равными нулю. Строя затем модели парной регрессии для рядов zi t, где i 1,2,, N, получим систему линейных уравнений.

b b1i t d i t, i 1,2,,N (2.75) где di – угол наклона линейной функции для ряда zi.

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

1. Имеется априорная информация о значении (полученная по b результатам “внешних сличений” эталона – “внешняя привязка” группового эталона). В этом случае из системы уравнений (2.75) немедленно получаем оценки b i, т.е. находим решение поставленной задачи.

2. Априорная информация о коэффициентах b i отсутствует, но для какого либо из рядов zi подтверждается гипотеза H0: b1i 0. В этом случае с большой долей вероятности можно считать, что b11 и b1i равны нулю, т.к. в противном случае придётся признать, что углы наклона линейных трендов частоты у опорного и i-го элементов группового эталона равны, а это крайне маловероятно (следует напомнить, что конструктивное исполнение и условия эксплуатации сводят вероятность такой ситуации практически к нулю).

В случае неприятия гипотезы H0 для всех рядов измерений, можно предположить что имеет место тренд частоты у опорного элемента. Тогда следует пересчитать ряды измерений на “новый опорный элемент”, что выполняется простыми линейными операциями.

i 1,2,, N z zi 1 z1 (2.76) i и вновь вернуться к процедуре проверки гипотез, описанной выше.

3. Процедура проверки гипотез не дала положительных результатов (обнаружены линейные тренды во всех временных рядах). Единственный возможный подход в этом случае – использовать метод наименьших квадратов для решения системы уравнений (2.75) с матрицей наблюдений H. Это выполняется посредством минимизации функционала вида n N I Т z i,t b11t b01 b1i t b0i (2.77) t 1 i Как уже упоминалось выше, оценка b1 может быть найдена как среднее значение результатов di.

N 1 d (2.78) b1 i n i В таком случае найденный вектор имеет минимальную норму [10].

N b Погрешность решения равна нулю только при выполнении условия 0, i i что было проверено при моделировании описанных процедур.

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

1. Выявляются детерминированные (линейные) тренды рядов относительных отклонений частоты водородных генераторов.

1.1. Если имеются результаты внешних сличений, то оценки параметров линейных функций ( bi0, bi1 ) находятся из анализа этих данных.

1.2. Если такой информации нет, проводится регрессионный анализ результатов взаимных измерений, т.е. строятся уравнения парной y1 (t ) yi (t ).

регрессии для всех рядов Проверяется гипотеза о равенстве нулю коэффициентов bi1. Если гипотеза не отвергается хотя бы для одного из рядов с индексом j, считается, что b1 0, b1 0. В противном случае считается, что линейная регрессия j отсутствует у элемента группового эталона, для которого модуль коэффициента b1j минимален.

1.3. Из недоопределённой системы линейных уравнений b11 bi1 b11i i 1,, N 1, (2.79) с учётом полученной на этапе 1.2 оценки b1 0, находятся оценки j линейных трендов частоты для всех элементов эталона.

2. Из рядов взаимных измерений устраняется влияние трендов.

3. По полученным скорректированным рядам измерений находятся МНК оценки значений относительных отклонений частоты опорного элемента для всех t, как это было описано ранее. Эти оценки совпадают со средним для каждого момента t. Оценки относительных значением z (t ) отклонений частоты для всех остальных элементов группового эталона (водородных стандартов) находятся по формуле (2.21).

4. Найденные ряды МНК-оценок относительных отклонений частот генераторов используются для построения моделей АРСС по методике Бокса-Дженкинса. Построенные модели используются в качестве первого приближения в процедуре минимизации функционала (2.22).

5. Полученные ряды оценок корректируются с учётом трендов, удалённых на этапе 2.

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

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

2.6. Выводы Задача оценивания состояния в системах с неполной матрицей наблюдения может быть решена с привлечением в качестве дополнительного канала наблюдений прогнозов, построенных на основе моделей процессов. В качестве таких моделей целесообразно применить стохастические динамические модели – модели АРСС. Подобный алгоритм фильтрации можно рассматривать как субоптимальный фильтр Калмана, построенный методами понижения порядка фильтра и декомпозиции вектора состояния системы.

Модели АРСС достаточно распространены и методология их применения в общем случае хорошо разработана. Однако применительно к системам с неполной матрицей наблюдения в целом и рассматриваемой задаче в частности, возникает целый ряд сложностей, не позволяющих реализовать классическую методику Бокса-Дженкинса. Это связано с отсутствием исходных рядов наблюдений для построения АКФ и ЧАКФ. Предложен оригинальный подход к решению этой проблемы, основанный на использовании для предварительной идентификации моделей АРСС рядов “грубых” МНК-оценок, полученных по алгоритму среднего арифметического.

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

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

3. Программный комплекс оценивания состояния эталонов времени и частоты по результатам взаимных измерений 3.1. Архитектура приложения и применяемые технологии В предыдущей главе была предложена методика оценивания относительных отклонений частоты в групповых эталонах времени и частоты, основанная на использовании в процедуре оценивания прогнозирующих моделей – моделей авторегрессии – проинтегрированного скользящего среднего (АРПСС). В процессе разработки данной методики, возникла необходимость её тестирования и испытания - как в режиме моделирования, так и при работе с реальными данными, полученными от эталона. Ввиду весьма существенной трудоёмкости производимых вычислений (поскольку речь идёт о численной оптимизации функций многих переменных), реализация такого алгоритма возможна только с использованием ЭВМ. Возникла необходимость разработки программного комплекса оценивания состояния эталонов времени и частоты.

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

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

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

На практике, среди основных способов обеспечения кросс платформенности можно выделить такие, как:

применение интерпретируемых и байт-код ориентированных языков программирования (Perl, Python, Java и др.) применение традиционных языков программирования в связке с обеспечивающими абстракцию от системных средств кросс платформенными библиотеками и фреймворками (например, Qt, GTK) Второй путь весьма популярен в коммерческих проектах и проектах, ориентированных на конечного потребителя. Он обеспечивает высокую производительность (в частности, за счёт отсутствия накладных расходов на интерпретацию) и отсутствие потребности в установке дополнительных программных средств (интерпретаторов, виртуальных машин) пользователем.

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


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

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

Одним из наиболее популярных языков такого типа является язык Perl.

Несмотря на изначальную ориентированность на обработку текстовых данных, Perl имеет очень развитую функциональность и широкую сферу применения, в том числе, и для решения научных задач [28]. Интерпретатор Perl портирован на большинство распространённых платформ, включая, конечно, ОС семейства Windows и *nix-системы, что позволяет выполнять написанные на нём скрипты (программы) под управлением любой из них без каких-либо дополнительных усилий. Разработка и отладка приложений на Perl занимает сравнительно мало времени, поскольку не требует перекомпиляции кода во время экспериментов или исправления ошибок, а сам язык обладает очень высокой смысловой ёмкостью (конструкции Perl выражают на единицу программного кода гораздо больше выполняемых ЭВМ действий, чем традиционные языки программирования, такие как Pascal, C++) и богатой встроенной функциональностью (в т.ч. развитой поддержкой динамической типизации и структур данных). При всём этом, Perl обладает очень высокой производительностью, в ряде случаев не слишком сильно уступая даже компилируемым языкам, таким как C++. Использование интерпретатора Perl бесплатно (он представляет собой СПО), ограничения в распоряжении написанными с его использованием программами отсутствуют. В силу развитых средств взаимодействия с системным окружением, программы на Perl хорошо интегрируются в существующую программную инфраструктуру.

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

Исходя из требований минимизации трудозатрат на разработку, было принято решение реализовать текстовый (“консольный”) интерфейс программы. Разработка графического интерфейса хотя и возможна, но не оправдана, т.к. для исследовательской работы он не требуется, а при внедрении программного комплекса в ИТ-инфраструктуру того или иного эталона могут возникнуть новые требования, предусмотреть которые заранее не представляется возможным. Возникающую же в период проведения исследований задачу визуализации результатов экспериментов (построение графиков, автокорреляционных и частных автокорреляционных функций и т.п.), а также ряд типовых статистических задач (проверка выборок на нормальность и др.) было решено возложить на внешнее средство. Поскольку в процессе проводимых исследований уже использовался пакет StatSoft Statistica версий 6.0 и 8.0, было решено выполнить интеграцию с ним.

Данный пакет имеет поддержку технологии OLE, позволяющий организовать взаимодействие с экземпляром процесса извне, в частности, из другой программы. В свою очередь, модуль Win32::OLE позволяет Perl приложениям в среде Windows получать доступ к функциональности программ, поддерживающих такой интерфейс. Возникающее при этом ограничение, связанное с платформой проведения исследований, не является проблемой, поскольку StatSoft Statistica и так является Windows-приложением, необходимым для проведения исследований.

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

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

1. Оценивание детерминированных трендов и их удаление из рядов измерений.

2. Структурная идентификация моделей авторегрессии – скользящего среднего (АРСС).

3. Идентификация параметров моделей и оценивание состояния.

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

Оценка относительного отклонения частоты опорного генератора на каждом такте работы алгоритма находится как n n yоп g i zi yi( 1 ) g i yоп yi yi( 1 ) (3.1) i 1 i где yi (1) – прогноз, вычисленный на предыдущем такте, gi 1 - вес i-го хранителя.

N Оценки для остальных генераторов вычисляются как yi yоп zi (3.2) Идентификация параметров моделей осуществляется путём оптимизации функционала вида n N I z t,i yt,1 yt,i 1.

(3.3) t 1 i Прогнозы yt,i для каждого i-го генератора вычисляются как:

yt 1 yt 1 2 yt 2... t yt p 1 at 1... q at q, (3.4) где – коэффициенты авторегрессии, p – количество включенных членов авторегрессии, – коэффициенты скользящего среднего, q – количество включенных членов скользящего среднего.

Ошибка прогноза предыдущего шага расценивается, как at-1 на следующем.

Фактически, при этом производится минимизация ошибки прогноза.

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

Оптимизация функционала осуществляется посредством метода градиентного спуска с дроблением шага [29]. Последовательно производится численное дифференцирование вперёд методом двух точек (как наиболее быстрым – высокая точность вычисления производных для данной задачи не требуется) по каждому коэффициенту моделей АРСС (с шагом x =0.0001), для чего вычисляется исходное значение функционала (3.3) и значение функционала после приращения по той переменной, по которой в данный момент ведётся дифференцирование, и находится их разность:

I I x x I x. (3.5) x x Данная процедура является достаточно ресурсоёмкой, поскольку ресурсоёмкой задачей является вычисление функционала (3.3) вообще:

требуется пересчёт оценок для каждого ряда i на все моменты времени t.

В результате, находим градиент функционала I:

I I I I I 1, p, 1, q.

,,,, (3.6), q p 1 Параллельно с этим, вычисляется евклидова норма градиента:

2 p q I I.

I (3.7) k 1 k l 1 l Имея градиент и норму градиента, выполняется шаг градиентного спуска:

I i I i1 I i k, (3.8) I i где 1 - шаг градиента (задан равным 0.025) Проверяется, выполняются ли условия:

I i I i1 0.1 k I i, (3.9) I i I i1 0, (3.10) k 1 109, (3.11) Условия (3.10) и (3.11) связаны, преимущественно, с точностью представления чисел в ЭВМ и исключают выполнение не имеющих смысла итераций.

Если выполняются все три условия, то происходит дробление шага k k 1. (3.12) Коэффициент дробления шага принят равным 0.5, т.е. шаг делится пополам.

После дробления шага вновь выполняется вычисление значения функционала (3.3), но уже с новым значением k и вновь проверяется соответствие условиям (3.9)-(3.11).

Если значение является приемлемым (условия перестают выполняться), то итерация завершается. Проверяется, является ли снижение функционала I I i достаточно существенным, в данном случае условие не выполняется, i если модуль разности меньше 0.005.

I i1 I i 0.005. (3.13) Если условие выполняется, то считается, что заданная точность достигнута и процесс оптимизации завершается. Также процесс завершается, если количество итераций превышает 80 – считается, что если условие не выполнено за такое количество итераций, то структура модели выбрана неудачно (не обеспечивает сходимости алгоритма) и дальнейшие затраты вычислительных ресурсов на неё не оправданы.

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

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

Кроме идентификации параметров моделей, имеет место задача идентификации структуры моделей. Ранее были описаны два предлагаемых способа решения этой задачи – способ, основанный на переборе наиболее вероятных структур моделей и последующем выборе оптимальной, и способ, базирующийся на использовании рядов предварительных оценок. В программном комплексе реализованы обе данные методики. Реализация первой из них, основанной на последовательном наращивании порядков модели, строится на том, что для моделей с количеством членов обоих типов от 0 до или 4 выполняется оптимизация функционала I, после чего вычисляются значения критерия вида:

I I3, (3.14) N (n k ( p q)) где N – число исходных рядов, n – число членов в каждом ряду, k=3 – добавочный коэффициент.


Результаты помещаются в массив, который сортируется по значению выражения (3.14), выбирается наименьший результат и используется для дальнейшей работы.

Альтернативный вариант основывается на том, что вычисляются как среднее арифметическое измерений предварительные оценки состояния вида:

1 N zi.

(3.15) y N i Оценки остальных генераторов вычисляются как (3.2).

После этого строятся ЧАКФ и АКФ полученных рядов, которые исследуются с целью определения структуры моделей по методике Бокса Дженкинса. С целью экономии времени на разработку, было решено использовать для визуализации АКФ и ЧАКФ средства StatSoft Statistica, как это описывалось ранее. Таким образом, при работе в таком варианте, комплекс запрашивает пользователя о необходимости экспорта в Statistica. Если пользователь отказывается, то ряды экспортируются в текстовые файлы (для визуализации иными средствами). В случае согласия пользователя, результаты экспортируются в Statistica, где по ним строятся ЧАКФ и АКФ. Пользователь анализирует полученные графики функций, идентифицирует порядки АР и СС и указывает их в программе оценивания состояния. Далее определённая таким образом структура используется для построения моделей и получения уточнённых оценок состояния системы.

Другой важной задачей является оценка и удаление линейных трендов.

Реализацией методики оценивания предусматривается оценивание и исключение из рядов взаимных измерений линейных трендов по заданным априорным сведениям о них (полученным в процессе “привязки” эталона постоянным составляющим и оцененным, например, через ряды внешних сличений эталона коэффициентам угла наклона), в том числе с использованием процедуры “дооценивания” тренда, весьма похожей по реализации на процедуру минимизации функционала I. Отличие состоит в виде используемой при этом модели – вместо моделей АРСС, используется линейная модель вида yiтр (t ) bi1t bi0. Оптимизируется методом градиентного спуска с дроблением шага функционал вида.

n N I Т z i,t b11t b01 b1i t b0i (3.16) t 1 i При этом шаг градиентного спуска выбирается значительно меньшим, чем в случае функционала (3.3), порядка 0.0005. Также меньшим выбирается и шаг дифференцирования x =0.00001. Это объясняется тем, что при больших t приращение коэффициента угла наклона даёт очень большой рост значения функции, таким образом, для обеспечения сходимости алгоритма требуется давать очень малое приращение. После определения параметров модели тренда, тренды удаляются из рядов измерений ziсл zi y1тр yiтр.

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

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

Случайная составляющая ряда генерируется моделями вида ytсл 1 yt 1... pd yt pd 1at 1... q at q at. (3.18) Сложность модели задаётся для каждого генератора индивидуально путём изменения параметров pi и qi. Генерация рядов выполняется в порядке возрастания номера члена, на каждом шаге t генерируется нормально at распределённая случайная величина с заданными характеристиками (нулевым математическим ожиданием и заданным с.к.о.), которая запоминается для следующих шагов (в зависимости от выбранного количества членов скользящего среднего).

Детерминированный тренд ряда моделируется при помощи линейных моделей ytтр b1t b 0. (3.19) Далее тренд накладывается на сгенерированную ранее стохастическую составляющую:

ytген ytтр ytсл. (3.20) Измерительная система, которая в принятой модели не вносит собственных шумов, моделируется в каждый момент времени t вычитанием соответствующих членов рядов ziген y1ген yiген. (3.21) В результате, на выходе блока генерации данных получаются ряды “измерений”, в известной мере соответствующие по своим характеристикам рядам, поступающим от измерительной системы эталонов времени и частоты.

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

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

Следует отметить, что размерность обрабатываемых данных (количество временных рядов и объём выборки) явных ограничений не имеет и, таким образом, ограничена только вычислительными возможностями системы и особенностями текущей версии интерпретатора Perl.

Результаты работы программы также выводятся в файлы сходного формата, кроме того, часть сервисной информации выводится в поток стандартного вывода (в консоль ОС). Для оперативной обработки и отображения результатов используется интеграция с пакетом Statistica компании StatSoft посредством технологии OLE. Естественно, что такая возможность доступна только при наличии указанной программы на ПК пользователя.

3.3. Результаты В результате реализации предлагаемой методики оценивания состояния был разработан программный комплекс, основными функциями которого является:

1. Генерация (синтез) данных по заданной модели АРСС и характеристикам случайной составляющей ряда.

2. Ввод реальных или ранее сгенерированных данных (рядов) из текстового файла.

3. Выявление и исключение линейных трендов из исходных данных.

4. Построение оценок среднего арифметического.

5. Идентификация порядка моделей АРСС:

способом перебора моделей, через предварительные оценки.

6. Построение оценок с использованием моделей АРСС:

подгонка моделей методом градиентного спуска:

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

оценивание состояния объекта.

7. Вывод промежуточных и окончательных результатов в текстовые файлы и на экран (в терминал).

8. Экспорт промежуточных и окончательных результатов в ППП StatSoft Statistica.

Приложение, реализованное с использованием языка программирования Perl версии 5.8 имеет простой текстовый (“консольный”) интерфейс и работает в диалоговом режиме.

Упрощённая схема функционирования программного комплекса приведена на рис. 3.1.

Ввод N-1 разностных рядов из файла Генерация N рядов по моделям АРПСС Получение N-1 разностных рядов Предварительная обработка данных Оценка и удаление трендов Построение рядов оценок по среднему Построение и вывод АКФ, ЧАКФ Пользователь Задание структуры модели Построение моделей рядов (Градиентный метод с дроблением шага) Оценивание состояния системы Постобработка и вывод Рис. 3.1 – Упрощённая схема функционирования программного комплекса Приложение реализовано с использованием структурного подхода к программированию и состоит из ряда модулей:

model_arima.pl – функциональность, связанная с построением моделей АРСС и оцениванием с их использованием;

model_average.pl – функциональность, связанная с построением оценок методом среднего арифметического;

model_generate.pl – генерация синтетических рядов;

model_io.pl – ввод и вывод данных;

model_main.pl – главный модуль;

model_single_arima.pl – функциональность, связанная с подгонкой моделей АРСС к одиночным рядам;

model_trends.pl – функциональность, связанная с исключением и наложением детерминированных трендов;

model_utl.pl – вспомогательные функции;

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

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

Особенно интересными являются результаты, связанные с выбором начальных приближений для процедуры оптимизации функционала I. Кроме наиболее простого варианта с выбором нулевых начальных приближений, была рассмотрена альтернативная вычислительная схема, основанная на использовании рядов предварительных оценок, выполненных по методу среднего арифметического. Было предложено изначально подгонять (раздельно) модели АРСС к рядам “грубых” оценок, получая, таким образом, начальные приближения для процедуры оптимизации, а уже после этого уточнять в процессе оптимизации функционала I. Вопреки имевшемуся мнению об избыточности подобных вычислений, такой подход оказался весьма эффективным. При одинаковых исходных данных (5 рядов реальных данных длиной 90 значений) проводились вычисления тем и другим способом.

Результаты в обоих случаях оказались полностью идентичны, однако суммарное время вычислений в случае использования начальных приближений, полученных через “грубые” оценки, оказалось существенно ниже (1,44 секунды против 4,08 секунды, при этом время, затраченное собственно на процедуру оценивания снизилось ориентировочно в 10 раз). Следует отметить, что внимательное изучение того, как использовалось время, показало, что в случае использования предварительных оценок процедура может быть ещё более оптимизирована, если свести к минимуму копирование данных в оперативной памяти.

Вывод приложения-профайлера приведён на рис. 3.2-3.3. Как можно model_params_estimation увидеть, функция оценивания состояния с вложенными в неё функциями отнимает львиную долю процессорного времени в случае использования нулевых начальных приближений. Из вызываемых ею функций наиболее значимы по затратам функция прогнозирования gen_next_observation и функция построения оценки частоты на момент времени t freq_estimation_on_t. Это является логичным, хотя возможна дополнительная оптимизация, например, путём использования inline-функций.

В случае с использованием предварительной подгонки моделей для получения начальных приближений картина отличается. Доля затрат на model_params_estimation функцию резко снижается, как снижается и абсолютное время её выполнения. Зато появляются затраты на вызовы функции индивидуальной подгонки моделей АРСС model_params_estimation_sing. При этом непропорционально увеличивается время выполнения функции copy_series_struct, копирующей структуру, содержащую временной ряд. Это связано с неоптимальной опытной реализацией model_params_estimation_sing и может быть устранено в целях оптимизации исключением потребности в полном копировании содержимого структуры временного ряда. Суммарное же время выполнения оказалось ниже более чем в два раза.

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

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

Разработанный программный комплекс для оценивания состояния эталонов времени и частоты по результатам взаимных измерений внедрён в деятельность ВСФ ВНИИФТРИ. Пройдена процедура регистрации созданного программного средства, в результате чего получено Свидетельство о государственной регистрации программы для ЭВМ №2012617062 (см.

Приложение 2).

Вместе с тем, имеется целый ряд направлений для доработки и дальнейшего развития данного продукта. В частности, исследование реальных данных, поступивших от эталона ГЭТ1-5 ВСФ ВНИИФТРИ, показало, что в ряде случаев ряды могут иметь не только линейный тренд, но и квадратичный.

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

4. Моделирование процессов оценивания состояния эталона 4.1. Структура специализированной системы моделирования Сложность проведения исследований в таких системах, как эталоны времени и частоты состоит в том, что возможность прямого контроля правильности результатов оценивания отсутствует. Действительно, непосредственное измерение оцениваемых величин невозможно: сущность эталона такова, что он хранит физическую величину с точностью, превосходящей точность других приборов. Для вторичных эталонов возможно использование результатов внешних сличений с Государственным эталоном [4], однако они сами по себе также содержат существенные погрешности и шумы, неизбежно возникающие в системе сличения.

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

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

Сформулируем требования к модели.

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

Требуемая математическая модель должна позволять получить N рядов истинных значений величин, а также на их основе – N-1 разностных рядов – результатов “взаимных измерений”. После обработки данных оценки могут быть сравнены с исходными рядами и таким образом определена точность алгоритма оценивания. Кроме того, должна обеспечиваться возможность наложения линейного тренда с заданным постоянным уровнем и углом наклона на полученные исходные ряды, с целью моделирования оценивания при наличии в данных линейных трендов. Вполне логичным является предположить, что модель эталона времени и частоты может быть декомпозирована на N моделей генераторов и модель измерительной системы.

При выборе класса моделей, описывающих водородные стандарты, следует исходить из следующих соображений. Во-первых, имевшие место ранее исследования [15,21], говорят о том, что реальные водородные стандарты частоты удовлетворительно описываются моделями авторегрессии – проинтегрированного скользящего среднего. Во-вторых, используемый математический аппарат опирается на применение моделей АРСС. Таким образом, в идеальном случае и при отсутствии детерминированного тренда наиболее точная оценка будет получена в случае, если при оценивании будет использоваться точно такая же модель генератора, как и при генерации синтетических данных. Исходя из всего сказанного, логично использовать для отладки методики оценивания состояния, основанной на авторегрессионных моделях, в качестве моделей генераторов модели АРСС вида:

ytсл 1 yt 1... p d yt p d 1at 1... q at q at (4.1) Сложность модели задаётся для каждого генератора индивидуально путём изменения параметров pi и qi. Генерация рядов выполняется в порядке возрастания номера члена, на каждом шаге t генерируется нормально распределённая случайная величина с заданными характеристиками at (нулевым математическим ожиданием и заданным с.к.о.), которая запоминается для следующих шагов (в зависимости от выбранного количества членов скользящего среднего).

Для моделирования детерминированного тренда можно использовать модели АРПСС, в которых порядок разности больше нуля. Однако наиболее просто детерминированный тренд ряда может быть смоделирован при помощи линейных моделей:

ytтр b1t b 0 (4.2) Член b 0 в таком случае описывает постоянную составляющую ряда (её можно было бы включить и в модели АРСС, однако это менее удобно с точки зрения интерпретации процесса).

Наложение тренда в таком случае производится простым суммированием детерминированной и стохастической составляющей:

ytген ytтр ytсл (4.3) Измерительная система, которая в принятой модели не вносит собственных шумов, моделируется в каждый момент времени t простым вычитанием:

ziген y1 yiген ген (4.4) Сгенерированные моделью измерения поступают на вход алгоритма оценивания так, как будто это – реальные данные. Алгоритм выполняет все этапы оценивания, включенные в него, и осуществляет построение оценок состояния элементов модели эталона времени и частоты yi для каждого момента времени t. На выходе алгоритма, таким образом, имеется N рядов оценок.

Имея N рядов исходных наблюдений yiген t и столько же рядов оценок yi t, подаём эти ряды на элемент сравнения. Сравнение может выполняться различными способами – как путём непосредственного сравнения значений рядов между собой (в т.ч. с вычислением рядов остатков и их последующим анализом), так и путём вычисления и последующего сравнения тех или иных статистик по рядам.

Структурная схема системы моделирования приведена на рис. 4.1.

Генерация стохастических составляющих y1 t y2 t y N 1 t y N t сл сл сл сл … Модель Наложение линейных трендов y1 t y N t y 2 t y N 1 t … Вычисление разностных рядов z1 t z 2 t z N 1 t … Алгоритм оценивания:



Pages:     | 1 || 3 | 4 |
 





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

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