Сравнительный анализ методов прогнозирования процессов потребления природного газа Текст научной статьи по специальности «Математика»

У cmammi розглядаються, аналiзуються i порiвнюються методи прогнозування взаемозв'язаних процеыв, як базуються на моделяхавторегресп - ттегрованого ковзно-го середнього з екзогенними змтними, метод групового врахування аргументiв, метод "Гусениця''-SSA на прикладi прогнозування процеыв споживання природного газу з ура-хуванням хронологiчних i метеорологiчних факторiв

Ключовi слова: процеси споживання, при-родний газ, АР1КСЕ, МГВА, "Гусениця''-SSA, порядок сезонностi, спектральнерозкладан-ня, вектор вкладень, траекторна матриця, полтом Колмогорова-Габора, зовтшне допо-внення

□-□

В статье рассматриваются, анализируются и сравниваются методы прогнозирования взаимосвязанных процессов, основанные на моделях авторегрессии - проинтегрированного скользящего среднего c экзогенными переменными, метод группового учета аргументов, метод "Гусеница''-SSA на примере прогнозирования процессов потребления природного газа с учетом хронологических и метеорологических факторов

Ключевые слова: процессы потребления, природный газ, АРПССЭ, МГУА, "Гусеница''-SSA, порядок сезонности, спектральное разложение, вектор вложения, траекторная матрица, полином Колмогорова-Габора,

внешнее дополнение

□-□

In article are considered, analyzed and compared interconnected processes forecasting methods, founded on the autoregressive -integrated moving average with exogenous variables model, group method of data handling, "Caterpillar"-SSA method on example of natural gas consumption processes forecasting with provision for both chronological and meteorological factors

The Keywords: consumption processes, natural gas, ARIMAX, GMDH, "Caterpillar"-SSA, seasonal order, spectral decomposition, vector of the embedding, trajectory matrix, polynomial of Kolmogorov-Gabor, external addition

УДК 681.513

сравнительным

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

А.Д. Тевяшев

Доктор технических наук, профессор, заведующий

кафедрой*

В.Н. Щелкалин

Аспирант*

*Кафедра прикладной математики Харьковский Национальный Университет радиоэлектроники пр. Ленина, 14, г. Харьков, Украина, 61200

1. Введение

Газотранспортная система (ГТС) Украины является пространственно распределенной системой, которая охватывает практически всю территорию Украины. Скорость движения газа в ГТС относительно неболь-

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

из необходимых условий повышения эффективности системы поддержки принятия решений по оперативно-диспетчерскому управлению режимами транспорта и распределения природного газа в ГТС Украины является создание и внедрение системы прогнозирования процессов потребления природного газа для всех категорий потребителей ГТС Украины с учётом влияния на них метеорологических и хронологических факторов.

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

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

2. рассматриваемые модели

а" =№ )■ )X...хащ-.. )=ГК. № )

- обобщенный оператор авторегрессии порядка п^

сП =с. №)■ с.2№ )X...хс-п. )=Пс. (qs' ) -

-с -с -с -с 4=1 -с

- обобщенный оператор скользящего среднего порядка

-с;

с0 - параметр, учитывающий детерминированный тренд;

ю|; = У У°2 -"У^гУ - стационарный случайный процесс, получаемый путем трансформации преобразованного исходного нестационарного процесса гу ; <О - порядок взятия разности; У. и qs' - упрощающие операторы такие, что

У^у =(1-qs' )■ гу = гу - г^.

Прогноз величины гу, вычисленный в момент 1 с упреждением 1 определяется по выражению

^У(1)=со гу+н -^ГсП

1=1

= |г 1 (1 -1),1 > 1,

2ы-,Л < 1,

(3)

0,1 > 1,

1 |гы_, - г 1+1-1-1 ()Д <1,

Далее будут приведены методы прогнозирования в обобщенном варианте, предполагая, что имеется прогнозируемый временной ряд yí и N экзогенных временных рядов х|., 1 = 1, N , однако для прогнозирования процессов газопотребления за N = 1 можно принять временной ряд изменения температуры окружающей среды. Также в роли экзогенного временного ряда при прогнозировании потребления природного газа может выступать временной ряд изменения скорости ветра.

2.1. модели авторегрессии - проинтегрированного скользящего среднего

В наиболее общем виде мультипликативная модель авторегрессии - проинтегрированного скользящего среднего, содержащая п. периодических компонент представляется в виде [1]:

а-,№)■ а^2)х...хУ22 ...У^гу =

= с1,^ с^2) х... х с-п.^4) ^

(1)

где аЦ1 (qЭl) = 1 - а^'' - a2q2■s' - ... - ап^"^1 - полином qs' степени па, определяющий авторегрессионную составляющую периодической компоненты с периодом Б'; с!'^5') = 1 - с^'' - c2q2■s' -... - c"'q"c - полином qs' степени пс, определяющий составляющую скользя-

щего среднего периодической компоненты с периодом

гу

Б'; et - остаточные ошибки модели, 1 - пронормиро-

временной

ванный от 0 до 1 по формуле гу = ———

ряд .

- У1

а„г (q)

где п;: = Х па ■, п; = Х пс ■.

где пу = пС + Х<О ■.'; с0 - общая константа, учи-

'=1 V • А V

тывающая детерминированный тренд; а' , ' = 1,па -коэффициенты обобщенного оператора авторегрессии, определяемые в два этапа согласно выражению аУУ Ы = а^С (q)У^1 У°2 . ..УО% , т.е вычисляются коэффициенты обобщенного оператора авторегрессии ""з п2 па

а П+,2 +-+К = (-1)п.-1 X X "Ха1-2 х ... х а-, (4)

£=0 '=0 1=0

причем

а0 = а2 = ... = а^ =-1; величина хоть и обозначается как период, но обычно считается равной 1, т.е. фактически имеем (п. -1) периодических компонент или считаем за период равный 1.

Затем вычисляются параметры а4У, аУ, ..., ауу оператора в результате выполнения п. итераций, осуществляемых по формуле п

аУ =

а = 011

^ = , па

(5)

Выражение (1) может быть записано в более компактной форме:

с0 + спС (q)

аП;(q)^ = с0 + ■ или =-^^— (2)

-а" ,1 = пх +1,пс + 8,,

' •> а 'а ''

причем вычисление по (5) для каждого повто-

1 П У

ряется на каждой итерации О раз с заменой а на а и пС на пС; + Э' . Для сокращения времени получения параметров аУ , аУ, ..., аУУ необходимо начинать вычисления по (5) с максимального значения периода .'_._

Неизвестные коэффициенты а', ) = 1,п. , ' = 1,па ; с', ) = 1,п. , ' = 1,п:с находятся при помощи алгоритма Марквардта для нелинейного метода наименьших квадратов.

Представим модель (2) = -

то же, гу =

с

с0+с„; (q)

■ или, что одно и

аУУ

в виде гу = а^= Ха' ■ е^, тогда дис-

персия прогнозов с упреждением 1 равна V (1) = о2 Ха2, где

коэффициенты а', ' = 0,1 -1 находятся по формуле

=1

=1

=1

а =

14=0,

УХ-а ,-Сп,i>о,

,=1

где а; = 0 , 1 > па ; с = 0 , 1 > пс ; ое - дисперсия

п!; сп= 0, 1 > пС ; о2 остаточных ошибок модели, определяемая в резуль тате решения задачи параметрической оптимизации модели.

2.2. метод авторегрессии - проинтегрированного скользящего среднего с экзогенными переменными

Математическая модель взаимосвязанных процессов, зависящих от метеорологических и хронологических факторов в операторной форме может быть представлена в виде [2]:

N ьПЬ1 (д) х

1м аП (д)- ^ Ьуу(д)■(6)

а1 Ь

N - количество экзогенных переменных, т1 - задержка 1-того экзогенного временного ряда по времени относительно прогнозируемого временного ряда т1; аП (д), ЬП (д) - полиномы д степеней па1 и пь1 соответственно; Ь^у (д), gnпс (д) - полиномы д степеней пу и п? соответственно; - остаточные ошибки модели;

- пронормированные от 0 до 1 экзогенные временные ряды х1, 1 = 1, N .

Выражение (6) в более компактной форме можно представить в виде:

N

а (д К=1 ь 1 (д)-2Х--т1 +с (дК,

1=1

где

а (д )=ьуу(д )-ГГаП1 (д)

Ь 1=1 ■

ь1 (д) = ЬПЬ1 (д)■ ЬУу(дМ1апа, (д)

С (д)=gПс(q)■П аП1(д), 1=

с упреждением :

N

пу+Еп,1

2У (1)= I ^<4 + ,=1

Ь0 ^х

п!+пь1 +1 п.к

1- I ь1-,=1

»с+х»,<

- I С,-

1=1

'1+1-1

2+,,, < 0,

< 0,

числяют в соответствии с уравнением (3). Значения е^,, , < 0 могут быть вычислены в три этапа:

N

1. Определяются значения 1У для 1 ^!(т1 + пь1) +1

N

г=1 ^=

м аПа1 (д)

N / ,

= 1(а1 ■ 1у-1+а2 ■ 1у-2+-+^ ■ ^)+

N / ■ ■ ■ \

+1(ь0<т1 - ь1 ■ 1Х-т1 -1 - - - ьПь1 <т1 )

Предшествующие значения 1У приравниваются к нулю.

2. Вычисляются значения

"У у ~У

11 = 1У -11 .

3. Непосредственно вычисление = еП осуществляется итерационно по формуле

пЬ ng

е 1 = е;-1 -1 ь^ц+1 ,

,=1

,=1

(7)

где е° = 1У-ц , 1 ^1(пЬ + d1 )) , ц - среднее значе-1=1

"У

ние временного ряда 11.

Предположим, что адекватная статистическая модель процесса 1Х может быть экономично представлена в классе общей мультипликативной модели типа (2), т.е.

1 сД(д)

Х1_ „х

еХ , 1 = 1,М

(10)

(8)

где еХ - процесс типа белого шума. Подставляя (10) в (6) и обозначая

1 ьПь1 (д) ■ сх1„пс (д) ■ дт1 ь1 (а) а (д) = —ь-1———СХ!——— = ^, (11)

Таким образом, получаем выражение для прогноза

в(д)=

аП,1 ЫЛч.Ы а1 (д)

ьУу(д),

(9)

представим модель (6) в виде

1у = 1ах1 (дК + Р(д)■ е4,

(12)

(13)

2У = •{ "" 1х=^ ' =Jet+l,j <0,

[1У ((),, > 0, +) [1х((),, > 0, ^ {0,, > 0,

где а,, ь , с, определяются из уравнений (8).

х1

Прогнозы 11 (1) получают либо непосредственно от метеостанций (прогнозы температуры воздуха в случае прогнозирования газопотребления), либо вы-

где еХ1 и статистически независимы 1 = 1,К. Тогда дисперсии прогнозов с упреждением 1 равны

V(l) = 11(а 1)2 ^о2х1 +о21в2 , где аХ1 аХ1 и в, - весо-

1=1 ,=0 е 2 ,=0 2

вые коэффициенты; о 1 , ое - дисперсии остаточных

ех

ошибок соответствующих моделей. Для вычисления весовых коэффициентов ах1 и в, сначала определяются параметров К , Ы , ..., К? +п 1 оператора

ь1 (д )=ьПь1 (д К^д )дт1 =

= Ыд^ - Ьдт1+1 - Ыдт1+2 - — - К?

Ь' Ы = ЬПЬ, (q=

= boqШ' - Ь^Ш' +1 - Ь12qm'+2 -... - Ь^ +п' q

х' Ь

Для этого используется формула

(14)

ЬИ

Ь0,1=о,

Ь0 ■ Сх' П-!Ь1 ■ с^М = ^

к=1

п'

(15)

ЬО ■ С'"-У Ьк ■ С' "к,1 = П' + 1,пх + п,',

0 х'1 ¿—I к х'Л-^ Ь' ' СХ' Ь"

где Сх'П = -1. ^ ^

Далее определяются параметры а1, а2, оператора

а' ^) = ^^Н'^, (q) =

п7 +п

= 1 - а^ - a2q2 - ... - а^. +п' q *х' "' Для этого используется формула

а = -У ак ■ ах'^к , к=0

гДе а0 =-1, ах'0=-1, а' =0, 1

V х -1

1 > п0 + п ' .

а<

>п

(16)

(17)

а'0=О,

х 1

Подставляя (14) и (16) в (11), имеем

1 - а^ - a2q -... - а^.+п^ х' а |( + + а2q + ...) =

= (^т' - Ь^1 - b2qm'+2 -... - Ы^Л^' +"Ь'|. (18) Если Ш' = 0 , то получаем

Ьо,1 = О, 1 „' „'

У ак а'-к - Ь',1 > 1,

(19)

где а'= 0, 1 >п°0' + па' ; Ь = 0, 1 >п^ + пЬ ' = 1,N. Если Ш' > 0 , то

0,1 < "', а' = - Ь0,1 = Ш',

У ак - Ь'-Ш',.1 >

(20)

оценивания оптимальных параметров модели - алгоритме Марквардта [2].

Автоматически оптимальные параметры структуры модели па' , пЬ', п^, п^ , Ш' можно определить методом случайного поиска в конечном множестве наиболее вероятных структур модели или методом целенаправленного перебора. Указанное множество формируется по значениям параметров п5, б1, б2 , ..., зп и ограничений на диапазон изменения параметров структур па' , пЬ', п^ , , Ш'. Параметры п^, п^ , б1, б2, ..., 8п могут также определяться путем анализа автокорреляционных и частных автокорреляционных функций, а параметры па' , пЬ', Ш' могут определяться путем анализа функции отклика на единичный импульс.

Таким образом, при помощи метода АРПСС можно вычислять прогнозы температуры воздуха, которые затем могут, в случае необходимости, использоваться для вычисления по моделям АРПССЭ прогнозов газопотребления.

2.3. Метод "Гусеница"-88А

Базовый метод 'Tусеница"-SSA (SSA) [4] состоит в преобразовании исходного ряда в многомерный, путем составления траекторной матрицы из так называемых векторов вложения, группировки членов сингулярного разложения получившейся траекторной матрицы и последующем её восстановлении. При этом часто оказывается возможным выделить отдельные аддитивные составляющие исходного ряда, такие как тренд, различные колебательные и периодические компоненты, а также шумовую компоненту. Также этим методом можно получать прогнозы последующих значений временных рядов по рекуррентному или векторному алгоритму.

Процедура вложения заключается в преобразовании исходного ^+1)-мерного временного ряда ^ - количество экзогенных временных рядов и один прогнозируемый) в последовательность L - мерных векторов ^ - ширина окна), число которых равно ( +1) К , К=п^+1, где п - длина временных рядов:

V (Уц,У-2,"-,У^2)Т , Х',1= (Х'_1,Х'_2,..,Х^_2)Т,

' = ^ 1 = 1Д

В результате формируем траекторную матрицу

х = [Х1:Х2:. :Хк:

= [Х1,1:Х1,2: . :Х1,К

Y] =

Х2,1 : Х2,2 :

: X,

•:ХК,1 :ХК,2 : . :Х№:^2: . : YK].(22)

где а' = 0, 1 > п! + па'; 1 0, 1 > < + пЬ'.

Этап сингулярного разложения и группировки заключается в формировании матриц разложения [5],

Веса оператора (12) находятся по формуле [1,1 = 0,

(21)

Z' =

и' ■ (и' )Т ■ Х1 :и' ■ (и' )Т ■ Х2:... :и' ■ (и')Т ■ Хк :и' ■ (и')Т ■ Y

,(23)

У - 8 м > o,

где Ь0 = 0, ' > п0; 8- = 0, 1 > п^ .

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

где и' - ьтый собственный вектор сингулярного разложения траекторной матрицы и сложении отобранных матриц разложения Z = Z' +... + Zj. Таким образом, для каждого из рядов получено разложение столбцов их траекторных матриц по общему базису

Этап диагонального усреднения переводит матрицы Z , ' = +1 состоящие из К столбцов от (' -1) К -ого до ' К - 1-ого матрицы Z в ряды Х0, Х1, ..., Хп-1,

к=1

1

к=1

к=1

подлежащие в дальнейшем прогнозированию, по формуле

1

~к X i =

—X z а-]+2,1 = 0, тт (L, К) - 2;

]=1 1

. „ и X Zj,i-j+2,1 = тт(L,K)- 1,тах(L,K)-1;(24) тт (L,K)

1

|х^,К)+1

-г X Zj,i-j+2,i = тах ),п -1,

j=k-K*+2

к = 1,К +1, причем ряды хо, Х1, ..., хп-1, 1 = 1,К соответствуют преобразованиям 1-ых временных рядов, соответствующих экзогенным факторам, при помощи сингулярного спектрального анали-

-N+1 -N+1 -N+1

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

-N+1 -N+1 -N+1

хо , Х1 , ..., Хп-1 будет соответствовать трендо-вой составляющей исходного временного ряда уо , у1, ..., уп-1. Если же базис (,...,№) выбирается из пар собственных векторов, соответствующих близким по величине собственным значениям тра-

-N+1 -N+1

екторной матрицы, то временной ряд хо , х1 , ...,

-N+1

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

Алгоритм рекуррентного прогноза [4] временного

ряда заключается в вычислении линейной ком, -N+1

бинации последних элементов временного ряда хо ,

- последний элемент вектора и1; г - количество элементов разложения; 1 - глубина прогноза.

Модификацией рекуррентного метода SSA-прогнозирования является векторное SSA-прогнозирование [4], которое в ряде случаев позволяет получать более точные прогнозы. Коротко векторный алгоритм SSA-прогнозирования временного ряда заключается в том, что, перед этапом диагонального усреднения, необходимо рекуррентно продолжить

-N+1

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

^N+14 Zl,i

»N+1 Z 2,1

»N+1 ZL,i

-N+1

х1 ,

-N+1

х п-1:

у 1=Х V

j=l

где а1 - элементы вектора т 1

(aL-1,aL-2,...,al) ="

, 1 = п,п +1 -1; х1 =

X UL

(25)

1 -х(

где

вектор, состоящий из пер-

(и1 и2 ... UL-:

вых (L-1) элементов 1-того собственного вектора и1 сингулярного разложения траекторной матрицы; ^

' и1

и?

+11 -Х(uL)? ■ ^^...аД

^N+1 Z 2,1-1

^N+1 Z 3,1-1

^N+1 ZL,i-1

(^-1,^-2,...,а1)

/«N+1 Л Z 2,1-1

^N+1 Z 3,1-1

^N+1 ZL,i-1

1=К- + 1, ...,К

N+1^ N+1

Zl,i Zl,i

Л N+1 -N+1

Z 2,1 = Z2,i

~ N+1 - N+1

1 ZL,i ^ ZL,i

1=1,...,К

(26)

После операции диагонального усреднения продол-

-N+1

женной матрицы Ъ получаем ряд 1 последних значений которого у ., 1 = п,п +1 -1 будут соответствовать искомым прогнозным значениям временного ряда yí.

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

Аналогично при помощи методов рекуррентного или векторного SSA-прогнозирования можно получать прогнозы экзогенных временных рядов температуры воздуха. Однако в этом случае траекторную матрицу следует формировать из векторов вложения только конкретного экзогенного временного ряда. 2.4. Метод группового учета аргументов В наиболее общем виде функцию, аппроксимирующую зависимость одного временного ряда от нескольких М других можно представить следующим образом:

у< = F (2,..,хМ).

В качестве такого аппроксиматора часто выступает полином Колмогорова-Габора [6]:

и

L-1

и

L-1

+

а

L-1

L-2

а

1

=1

М ММ МММ

Уt=аоа. ■х'+ХХVх' ■ х1+ХХХ а

¡=1 Л—1

=1 Л—1 к—1

¡Лк

х' ■ Х{ ■ хк + ■■■,(21)

так как с помощью такого полинома можно добиться достаточно точной аппроксимации любой дифференцируемой функции.

В методе группового учета аргументов эта сложная зависимость заменяется множеством простых функций:

У1« = f (,Х?); У?» = f (х^);.., у^1 = f (^М) ,(28)

где С = СМ, причем функция f всюду одинакова.

Очень часто в качестве функции f выбираются простые зависимости вида:

У1 (х1,х1 ) = а0 + а1 ■ х1 + а2 ■ х1 + а3 ■ х1 ■ х1

У (х1,х1 ) =

— а« + а ■ х^. + а о ■ х,. + а? ■ х,. ■ х^ + а,

(29)

(30)

+ а ■(х

связывающие только две переменные.

Модели (28) составляют первый ряд метода, из которых выбираются несколько Я наилучших по комбинированному критерию эффективности и внешнего дополнения, т.е. с одной стороны наилучшим образом аппроксимирующие обучающие данные, а с другой стороны показывающие наилучшие результаты на проверочной выборке - прогнозирующие наилучшим образом поведение функции на выборке, не охваченной экспериментом. Во втором ряду алгоритма полученные и отобранные на обучающей выборке значения у'/1 перенумеровываются по порядку и рассматриваются в качестве аргументов нового второго ряда:

у;<2) — f(у1<1>,у2<1>); у2<2> — f(«у*»);...;

у«« — f (.

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

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

Параметры частичных описаний (29) или (30) обычно вычисляются при помощи метода наименьших квадратов, однако ниже предложено использовать ро-

бастный метод оценивания для получения моделей адекватных процессу, который имеет следующий алгоритм [1]:

е1+1 = e¡+а ■ р1/2 ■ н-1 (е1 )ь (е1),

где 1 - номер итерации, е - вектор параметров, такой что е = (а0,а1,а2,а^ для частичных описаний вида (29) и е = (а0,а1,а2,а3,а4,а^ -для описаний вида (30),

П "1

ь (е)—|> ( (е))-7'-1, н (е)=£ ^Т^ (V' (

' ^—(^«—йх, I: «(х )=

1,х I С, х,|х| < С, —1, х < с;

где а - постоянная, обычно выбираемая равной 0,5; с - любое число между 1 и 2; п1 < п , так как объем выборок уменьшается за счет взятия задержек по времени. Начальными значениями оценок для этого алгоритма служат оценки е', получаемые по следующему алгоритму, принадлежащему семейству алгоритмов фильтрации Калмана [1]:

е 1—е 1+St ■ 7'—1 ■(у'-(е 1—1 )ТЛ-1);

S — S --

■ гт ■S

-1 п-1 >-У1 t —

1+ 2''-1 ■ 31-4 ■

, 1 — 1,п1. 30 — Е.

В этом алгоритме zí - вектор следующего вида

—(1,х'1,х'1,х1 ■ х{) для частичных описаний вида (29) / 2 2\Т

и —(1,х1,х'1,х'1 -х',(х') ,(х') - для частичных описаний типа (30).

Таким образом, метод группового учета аргументов позволяет определять прогнозы как отдельно экзогенных временных рядов (временных рядов температуры воздуха), так и прогнозы зависимого временного ряда (ряда газопотребления) совместно с экзогенными. При этом М>^ так как помимо экзогенных переменных можно использовать их задержки по времени (лаговые эффекты) на несколько единиц, задержки по времени самого прогнозируемого зависимого временного ряда, а также временные ряды, составленные из исходных зависимого и экзогенных с элементами, взятыми с задержками равными периодам , 1 — 1,п5 .

3. Сравнительный анализ результатов прогнозирования

Рассмотрим основные этапы построения АРПССЭ модели процесса суточного газопотребления с учетом изменения температуры воздуха.

Средняя ошибка такого оперативного прогноза не должна превышать 3%. На рис. 1 представлены временные ряды среднесуточных изменений температуры воздуха и потребления газа, взятые с 1 сентября по 31 мая.

1—1

Рис. 1. Графики изменений температуры воздуха и потребления газа

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

Рис. 2. Автокорреляционная и частная автокорреляционная функции изменения температуры воздуха

Рис. 3. позволяет говорить, что для выравнивания спектра, необходимо брать модель процесса изменения температуры воздуха вида АРПСС (3, 1, 2)

Рис. 3. Автокорреляционная и частная автокорреляционная функции изменения температуры воздуха после взятия первых разностей

Автокорреляционная и частная автокорреляционная функции "выбеленного" временного ряда изменения температуры говорят об адекватности модели выравнивания спектра [3].

Рис. 4. Выбеленный временной ряд изменения температуры воздуха и его автокорреляционная и частная автокорреляционная функции

По функции отклика на единичный импульс в данном случае сложно судить о структуре передаточной функции, поэтому перебором всех моделей кандидатов со сложностью не превышающей 2 и анализируя СКО и Ц-статистику остаточных ошибок, их корреляцию с экзогенными данными была получена модель АРПССЭ

^ -0,447-0,155q и А ооп ч

= —1 -0 425q— +(1 + 0,222^е4. Для прогнозирования температуры воздуха модель получилась вида: 1 + 0,^ - 0,635q2

Ух =-

-е, . Методом целена-

1 + 0,257q-0,377q2 + 0,207q3 4 правленного перебора вариантов итеративно осуществляется выбор оптимального объема обучающей выборки, который, для рассматриваемого примера, составил 101 значение.

Рис. 5. Функция отклика на единичный импульс

С середины марта, когда среднесуточная температура становится более 4°С, связь процесса потребления газа с изменением среднесуточной температуры воздуха теряется (рис. 6, после значения 93) и изменяется структура модели, что видно из графиков изменения параметров модели, поэтому прогнозировался участок времени до возникновения данной разладки.

Л ,Л

™1ЙЯ

Рис. 6. Графики изменения параметров модели

Средняя ошибка прогнозирования составила 2,018%.

Рис. 7. Процентные ошибки прогнозов

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

№

! йй

нФ

Рис. 8. АРПССЭ Прогнозы потребления природного газа с учетом изменения температуры воздуха и 95%-ые доверительные интервалы

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

Средняя ошибка прогнозирования составила 6.01%.

Рис. 9. Прогнозирование потребления природного газа с учетом изменения температуры воздуха методом "Гусеница"-SSA

При прогнозировании методом МГУА, использовался тот же объем выборок, только 30% выборки использовалось в качестве внешнего дополнения. В результате были получены следующие три слоя частичных описаний:

у!<3> —-0.012 +1.32 ■ у'(2) - 0.334 ■ у9(2) +

+5.46-(у1<2>)2 -9.945■ у'<2> ■ у9<2> + 4.484-(у9<2>) ;

у'(2> — 0.016 + 0.618 ■ У2^ + 0.312 ■ у4(1> -

-9.5-(у2«)2 +18.308■ у2«■ у4« -8.166(

у92> — 0.011 + 0.894 ■ у5(1> + 0.104 ■ у^ -

-3.35()2 + 6.493■ у?«■ у™ -3.143-(у?«);

у2(1 — -0.024 + 0.411 ■ х'-1 +1.212 ■ у1-2 --0.011 ■ х2-1 + 0.248 ■ х1-1 ■ у1-2 + 0.231 ■ у2-2; у4(1 — 0.011 - 0.105 ■ х1-2 + 0.819 ■ у1-2 --0.19 ■ х2-2 -1.412 ■ х'-2 ■ У'-2 - 0.112 ■ у2-2; у5(1> — 0.0030 - 0.238 ■ у1-1 +1.182 ■ у1-2 --0.194 ■ у2-1 + 0.118 ■ У'-1 ■ У'-2 - 0.042 ■ у2-2; у3(1> — 0.032 - 0.488 ■ х1-2 + 0.586 ■ у1-1 --0.135■ х2-2 - 0.411. х^ ■ У1-1 -0.341 ■ у2-±.

Рис. 10. Прогнозы и средние проценты ошибок прогнозирования по МГУА

Метод показал хорошие результаты, однако средний процент ошибки прогнозирования составил 2.6%.

4. Быводы

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

Метод "Гусеница"-SSA, несмотря на определенное количество публикаций, для прогнозирования потребления природного газа получил неудовлетворительные результаты. Основной проблемой, с которой можно столкнуться при прогнозировании газопотребления данным методом заключается в том, что метод разлагает временные ряды на различные компоненты адекватно не для всех временных рядов. Малые параметры ширины окна - 2 и порядка разложения - 1 вызвано тем, что модель определила псевдосезонную компоненту с периодом порядка 8, и, при больших значениях параметра разложения прогнозировала эту составляющую, однако анализ автокорреляционной и частной автокорреляционной функций показал отсутствие даже недельной сезонной компоненты. Данный метод может достаточно хорошо применяться для анализа временных рядов с ярко выраженными трендовой и сезонными составляющими.

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

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

Необходимо продолжить исследования применения АРПССЭ-моделей для прогнозирования процес-

сов потребления природного газа с учётом влияния на них метеорологических и хронологических факторов на значительно большем объёме экспериментальных данных (от трёх до пяти лет) для всех категорий потребителей ГТС Украины.

Литература

1. Евдокмимов А.Г., Тевяшев А.Д. Оперативное управле-

ние потокораспределением в инженерных сетях. - К.: 1979.

2. Бокс Дж., Дженкинс Г. Анализ временных рядов. Прогноз

и управление. Вып. 2. Пер. с англ. - М.: "Мир", 1974. -197 с.

3. Бэнн Д.В., Фармер Е.Д. Сравнительные модели прогно-

зирования электрической нагрузки: Пер. с англ. - М.: Энергоатомиздат, 1987. - 200 с.: ил.

4. Голяндина Н.Э. Метод «Гусеница»-SSA: прогноз времен-

ных рядов: Учеб. пособие. СПб: Изд-во СПбГУ, 2004. 52 с.

5. Н.Голяндина, В.Некруткин, Д.Степанов. Варианты метода

"Гусеница"-SSA для анализа многомерных временных рядов. Труды II Международной конференции "Идентификация систем и задачи управления" SICPRO'03. Москва, 2003, с. 2139-2168.

6. Ивахненко А.Г. Долгосрочное прогнозирование и управ-

ление сложными системами. - К.: "Техшка", 1975. - 311 с.

7. Кашьяп Р.Л., Рао А.Р. Построение динамических стохасти-

ческих моделей по экспериментальным данным. - М.: Наука., 1983. - 384 с.