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

лектуальных систем. СПб.: Питер, 2001. 384 с.

6. Копайгородский А.Н. Методы, модели и программные средства построения информационной инфраструктуры исследований в энергетике: авторефер. дис. ... канд. техн. наук. Иркутск: ИСЭМ СО РАН, 2008. 25 с.

7. Рубашкин В.Ш. Онтологическая семантика: Знания. Онтологии. Онтологически ориентированные методы информационного анализа текстов. М.: ФИЗМАТЛИТ, 2012. 346 с.

8. Копайгородский А.Н., Массель Л.В. Разработка и интеграция основных компонентов информационной инфраструктуры научных исследований // Вестник ИрГТУ. 2006. № 2 (26). С. 20-24.

9. Воропай Н.И., Массель Л.В. ИТ-инфраструктура системных исследований в энергетике и предоставление ИТ-

услуг // Известия АН - Энергетика. 2006. № 3. С. 86-93.

10. Массель Л.В., Копайгородский А.Н. Технологии и система хранения данных и знаний для исследований в энергетике: мат-лы Всерос. конф. «Современные информационные технологии для научных исследований». Магадан: СВНЦ ДВО РАН, 2008. С. 64-66.

11. Копайгородский А.Н. Интеграция данных в исследованиях энергетики на основе онтологий: труды XVII Байкальской Всерос. конф. Иркутск: ИСЭМ СО РАН, 2012. Т. III. С. 62-68.

12. Копайгородский А.Н. Проектирование и реализация системы графического моделирования: труды XV Байкальской Всерос. конф. Иркутск: ИСЭМ СО РАН, 2010. Ч. 3. С. 22-28.

УДК 621.311.1

СМЕШАННЫЕ АВТОРЕГРЕССИОННЫЕ МОДЕЛИ И ПРОГНОЗИРОВАНИЕ ПРОЦЕССА ВЫРАБОТКИ ПАРА

© В.Г. Хапусов1, А.В. Баев2

Иркутский государственный технический университет, 664074, Россия, г. Иркутск, ул. Лермонтова, 83.

Рассматривается применение известной методики Бокса Д.Ж. и Дженкинса Г. для идентификации процесса выработки пара, который носит нестандартный характер. Поэтому к нему была подобрана модель авторегрессии проинтегрированного скользящего среднего (АРПСС) порядка (1 1 1). Построенные для различных характеристик процесса модели имеют не только самостоятельное значение, но и могут быть использованы для краткосрочного прогноза, а следовательно, и для оперативного управления производством. Ил. 2. Табл. 2. Библиогр. 1 назв.

Ключевые слова: идентификации; оценивание; диагностическая проверка; прогнозирование.

MIXED AUTOREGRESSIVE MODELS AND STEAM PRODUCTION FORECAST V.G. Khapusov, A.V. Baev

Irkutsk State Technical University, 83 Lermontov St., Irkutsk, 664074, Russia.

The article deals with the application of a well-known Box-Jenkins Forecasting technique to identify the non-standard process of steam generation. For this reason an autoregressive integrated moving average model (ARIMA) of (1 1 1) order has been chosen. Models built for various characteristics of the process are self-sufficient, but also can be used for short-term forecasting, and consequently, for operating production management. 2 figures. 2 tables. 1 source.

Key words: identification; evaluation; diagnostic test; forecasting.

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

ся к точности поддержания температуры и давления пара.

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

1Хапусов Владимир Георгиевич, доктор технических наук, профессор кафедры автоматизации производственных процессов, тел.: 9148883081, e-mail: hapusov@yandex.ru

Khapusov Vladimir, Doctor of technical sciences, Professor of the Department of Automation of Production Processes, tel.: 9148883081, e-mail: hapusov@yandex.ru

2Баев Анатолий Васильевич, кандидат технических наук, зав. кафедрой автоматизации производственных процессов, тел.: (3952) 405243, e-mail: baev@istu.edu

Baev Anatoliy, Candidate of technical sciences, Head of the Department of Automation of Production Processes, tel.: (3952) 405243, e-mail: baev@istu.edu

тистические, а не детерминированные явления.

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

Zt = (pi Zt-1 + (р2Z t-2 +... + PpZ t-p + at-61 at_x

-6

6 q a

2 at-2 ".-6q at-q'

или (р (В) 11 = 0(В) а, (2)

где а - случайный процесс; В - оператор сдвига назад, т.е. В2( = Ъи; рьр2,..., рр - параметры авторегрессии; 01 ,02,-, -параметры авторегрессии.

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

ф (В) 21 = 0 (В) а{, (3)

где ф (В) = Vе р (В); V - оператор разности; V = -Ъи; р(В) - оператор авторегрессии. Предполагается, что этот оператор стационарен, то есть корни уравнения р(В) = 0 лежат вне единичного круга; ф(В) -обобщенный оператор авторегрессии. Это нестационарный оператор, у которого Ь корней уравнения ф (В) = 0 равны единице.

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

Ро=(Рю (о >••••(()), несезонных параметров скользящего среднего 0 = (010,02О,....,0?о), начальной оценке общего постоянного члена 0О , начальной оценке дисперсии белого шума с2.

В случае нестационарности временного ряда вычисляются разности до тех пор, пока временной ряд не окажется стационарным, и отсюда получают порядок операций взятия несезонных разностей й > 0 .

Для каждого набора значений временной ряд

преобразуется следующим образом:

[(г, + т )Я Л* 0,

1п ( г + т) Л = ° где Л и т - параметры преобразования. Если Л Ф 0, параметр т выбран так, что г + т положительно для всех 1 Если Л=1, т приравнивается нулю, так что zt остаются неизменными. Затем находится разностный ряд, такой, что

о^^1^, 1=1, 2,....., п,

(1) где Vщ=щ—щ_1, п = N - С - число значений ря-

да о .

Далее вычисляются следующие величины. Среднее значение и дисперсия

1 _п_

1 2 2 = , зю= с0, где зо — дисперсия ряда о

п t=1

Автоковариационная функция

п—к

С,

k

J n-k

- ЕЦ-®)(®t+k

n t=l

где к = 0,1......К , К - максимальная задержка.

Автокорреляционная функция

rk = — ■

(4)

Частная автокорреляционная функция Г, l = l

Pu =<

l-l

rl -

Ер--Л-j , l = 2,3.

■ , L,

(5)

l-l

1 — ^рг—1

з =1

где р = РI—1,1 —ри(1—и—1 , \ = 1 , 2...... I-1, Ь

< К - максимальная задержка.

Частная автокорреляционная функция может быть оценена путем вычисления оценки наименьших квадратов последнего из параметров авторегрессии ри в

последовательности процессов авторегрессии АР(1), I

= 1, 2,......, подгоняемых к данным.

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

c

c

0

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

параметров /,р,6 и (Г2а в модели.

р( В —л)=б( В) а,

где - преобразованный временной

ряд; /- среднее значение ряда С0{,

Р = (Рю,--;Рро) - начальные оценки несезонных параметров авторегрессии; 6 = (б10 ,62о ,...,6?о)

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

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

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

соответствующим оценкам наименьших квадратов, согласно формуле

^ (k)=с^ (k)/(0),

(6)

где

n—k

Caa ( k ) = 1 X( at — a )( at+k — a )

n

t=1

a

=Xa, k = 0,1,...,K,

n

t=1

наконец, статистика % вычисляется как

к

%2 = пУ (к) и сравнивается с % - распреде-

к=1

лением с у = К — ё — р — д степенями свободы. Если полученная модель удовлетворительна, то Хр < X2 (расчетная меньше табличной).

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

После получения оценок наименьших квадратов

параметров (60, р1,ра,....,рр, 6,6,.,6) общей

модели вычисляются:

[ Z t+| ] =[ Z t (l)] Z t+|.| ] + (p2 [ Z t+l-2 ]

+... + (p[ Z t+i-p ] + [at +|] - 91

[at +1-1] - 92[at +1-2] - ... - 9q[at +4

(7)

2 \ (I) - прогноз с минимальной среднеквадратичной ошибкой для упреждения I.

Верхний и нижний вероятностные пределы

2 м. (±) = 2 t (I) ± иЛ^(Г), (8)

где и =1,96 или 2,58, в зависимости от того, лежит ли будущее между этими пределами с вероятностью в интервале 0,95 или 0,99 соответственно. Функция дисперсии

V(/) = ^а 1 у1 ща ., а г\ 1 7 = 0

где

щ j — веса, равны

1, j = 0

j

i щj

i =1

— 9 ,J)l

j"

В качестве объекта исследования был выбран котельный агрегат БКЗ-420-140-6, оборудованный четырьмя пылеприготовительными установками. БКЗ-420-140-6 однобарабанный, вертикально-водотрубный с естественной циркуляцией, предназначен для сжигания азейских бурых углей. Данные, собранные в течение 2-х часов наблюдений за нормальным ходом процесса выработки пара с интервалом отсчета 10 с, были подвергнуты статистическому анализу.

В период пассивного эксперимента контролировались технологические параметры, приведенные в табл. 1.

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

„ _ ,0,936 ,0,557

= (+0,017) ' 1 — (+0,039) ' а—1 + а ■

Учитывая, что двойное значение стандартной ошибки (0,017) для коэффициента регрессии не превышает значение самого коэффициента (0,935), мож-

^ ,0,936 Л

но сказать, что коэффициент рр = (+0017) значим.

Аналогичный вывод можно сделать и для

= (0,557 ) 61 = (+0,039)'

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

Таблица 1

Технологические параметры процесса выработки пара_

Параметр Условное обозначение Единицы измерения

Давление перегретого пара Рпер. п. МПа

Температура дымовых газов (в поворотной камере) Тд. г.п.к ° С

Температура перегретого пара Тпер. п. ° С

Расход перегретого пара Спар м3/ч

Уровень воды в барабане ^в в.бар Мм

Содержание О2 в дымовых газах О2 %

Разряжения в топке Р г разр. Кгс/м2

Расход топлива ( обороты двигателя ПСУ) ОПСУ Обор/мин

Расход воздуха (ток двигателя ВДН) Ток ДВ А

Расход дымового газа (ток двигателя Д) Ток ДС А

Рис. 1. Временной ряд расхода перегретого пара вПара

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

Наконец, статистика % сравнивается с х2 - распределением с у = К- й- р- ц степенями свободы. Если полученная модель удовлетворительна, то

Х2р ^ X2 (расчетная меньше табличной). В табл. 2

приведено табличное значение - распределения и

2

статистика г = 7,09.

Таблица 2

Модели временных рядов_

Параметр Математическая запись модели Модель Число степеней свободы X p

С пара „ .0,936 , „ ,0,557 . ^ Уг( - (+0,017 >' У2г-1 -(+0,039 >' аг-1 + аг АРПСС (1 1 1) 12 7,09 Х2=21,0

О2 справа „ Л,618 , „ /-0,787 , „ ,0,063 . „ У2г = (+0,095 ) • У2г-1 + (+0,112 ) •У2г-2 + (+0,042 )• У2г-3 - 1,63^ -0,773 (+0,088 ) • аг-1 (+0,094 ) • аг-2 + аг АРПСС (3 1 2 ) 9 11,31 Х2=16,9

О2 слева „ ,0.991 . = ( + а г '+ 0.005) г -1 г АР (1 0 0 ) 14 6,64 Х2=23,7

Рразр. „ ,0,780 . „ ,-0,420 . „ ,-0,026 . „ У2г = (+0,183 ) • Угг-1 + (+0,093 ) • Угг-2 + (+0,050 ) • Угг-3 + АРПСС 8 5,58

справа ^ ,-0,194 , „ ,1,129 , ,-0,580, ^ +(+0,046)' У2г-4 -(+0,185^аг-1 -(+0,129 )'аг-2 + аг (4 1 2) Х2=15,5

Рразр. слева У7 ,0,717 . ^ У2г = +0,026 " аг-1 + аг ПСС (0 1 1) 13 15,00 Х2=22,4

Тд.г. в п.к. справа „ ,0,973 . „ ,1,227 , ,-0,476 . У2г =(+0,010 ) • У2г-1 - (+0,037 ) • аг-1 - (+0,036 )' аг-2 + аг АРПСС (1 1 2) 11 15,54 Х2=19,7

Тд.г. в п.к. слева „ ,0,957 . „ ,1,063 . ,-0,424 . ^ = (+0,012 ) • У2г-1 (+0,040 ) • аг-1 (+0,039 ) • аг-2 + аг АРПСС (1 1 2) 11 13,29 Х2=19,7

„ .-0,705 . „ .0,409 ^ ,0,954 , „ У2г = (±0,087) • У2г-1 + (±0,069) • Уг(-2 + (±0,077) • У2г-3 +

Ток ДС-А 0,130 -0,616 0,432 + (+0,042)' У2г-4 - (+0,081 ),аг-1 - (+0,080),аг-2 - АРПСС (4 1 3) 7 12,5 Х2=14,1

0,814 - (+0,088 ),аг-3 + аг

Ток ДС-Б „ ,0,758 , „ ^ ,0,031 , „ /0,157 , „ У2Г = (+0,045 ) • У2Г-1 + (+0,047 > • -2 + (+0,038 ) ' У2г-3 - /0,909 , ^ - (+0,027 )аг -1 + аг АРПСС (3 1 1) 10 15,78 Х2=18,3

Ток ДВ-А „ ,0,579 , ,-0,068 . У2г = +0,038 ) • аг-1 - (+0,038 ) • аг-2 + аг АРПСС (0 1 2) 12 20,95 Х2=21

„ ,0,098 . „ ,-0,302 . „ ,0,671 . „ = (+0,108 ) • У2Г-1 + (+0,064 ) • Угг-2 + (+0,083 ) • УгГ-3 +

Ток ДВ-Б ,0,164 , „ ,0,10^ , ,-0,479 , +(+0,043 ) • У2Г-4 - (+0,106 ) • аГ-1 - (+0,075 ) ' аг-2 - 0,717 - (+0,102 )' аг-3 + аг АРПСС (4 1 3) 7 13,57 Х2=14,1

ОПСУ-А „ -0,533 .0,206 . _ ,0,91^ = (+0,058) • У2г-1 +(+0,058) • Уг/-2 - (+0,042)' аг-1 + аг АРПСС (2 1 1) 11 14,57 Х2=19,7

ОПСУ-Б „ ,0,846 . „ ,0,658 . ^ =(+0,043 ) • -1 - (+0,059 ) • аг-1 + аг АРПСС (1 1 1) 12 19,83 Х2=21

ОПСУ-В „ ,0,549 . „ ,0,179 . „ ,0,371 . У2г = (+0,114) • Уг1-1 + (+0,057) • Уг1 -2 - (+0,114)' аг-1 + аг АРПСС (2 1 1) 11 17,96 Х2=19,7

ОПСУ-Г 2 г =(:±с9с)03) •-1 + аг АР (1 0 0) 14 4,80 Х2=23,7

^в. в бар. п / 0,792 \ /-0,312ч , / 0,317 \ ^/-05814 , АРПСС (3 1 1) 11 13,06 Х2=19,7

Р в бараб. п /0,94^^ , / 1,066 \ . 1 = \+0,017/ ^ + ч+0,042] ^ + /-0,287ч + ч+0,039/ ^ а'-2 + аt АРПСС (1 1 2) 11 16,46 Х2=19,7

Р пер. п. У21 = 1+0,103] У21-1 + (+0,098] У21-2 + , / 1,45^ , /-0,586ч , АРПСС (2 1 2) 10 13,27 I2 = 18,3

Т пер. п. = 1+0,177/ + (+0,087] У21-2 + ( 0,036 у7г +(-0,099)^г + (+0,099] 1-3 + (+0,071] У21-4 АРПСС 6 11,34

,/-0,11^ /-0,738ч + 1+0,088/ ^ а'-1 + 1+0,099/ ^ а'-2 -0,082 -0,654 + 1+0,109] •а'-з + 1+0,086/^1-4 + а^ (4 1 4) I2 = 12,6

Forecasts; MDdaisflJO.I} Seasonal lag 12 Input: ВПАРА : Df-lJ Start of origin: 2 End □[ origin 700

1

1 1 II

г 11

J А

1

т

500 525 550 575 GOD C25 ОБО 075 7DD 725 750 - Observed---Fwecest ...... ± ЭО.ЦООО»

CaseNo. Forecasts; Model:(1.0,1) Seasonal lag: 12 (Данные) Input: GПАРА : D(-1) Start of origin: 2 End of origin: 700

Forecast Lower 90.0000% Upper 90.0000% Std.Err.

701 -0,677284 -2,39421 1.039644 1.042431

702 -0.6338641 -2,46982 1.202093 1 114699

703 -0,593227 -2,52743 1,340978 1,174350

704 -0,555196 -0,519603 -0.486291 -0,455115 -0.425938 -0,398632 -0,373076 -2,57153 1.461135 1.224213

705 -2,60521 1,506006 1,266275

706 -2,63074 1.658159 1,302000

707 -2,64981 1,739578 1,332506

708 -2,66371 1.811837 1,358663

709 -2,67347 1,876208 1,381166

710 -2,67989 1.933739 1.400580

Рис. 2. График прогноза преобразованного временного ряда расхода перегретого пара впара и таблица

прогнозируемых значений

Табличное значение - распределения для 12 уровней свободы при уровне значимости 0,05 равно 21,0, следовательно, модель адекватна.

Аналогично были построены модели других временных рядов, которые приведены в табл. 2.

Во многих случаях необходимо знать прогноз на несколько шагов вперед. Используя формулу (7) для вычисления прогнозируемого значения, в табл. 3 показаны результаты расчета прогнозов для преобразо-

ванного ряда Спара на момент 1 = 700 с упреждением I = 1, 2, 3, 4, 5, 6, 7, 8, 9.

Дисперсия прогноза по расходу перегретого пара о-2 = 1,086.

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

Статья поступила 11.11.2014 г.

Библиографический список

1. Бокс Д.Ж., Дженкинс Г. Анализ временных рядов, прогноз и управление. М.: Мир, 1974. 604 с.