Факторизованная форма рекуррентного алгоритма текущего регрессионного анализа Текст научной статьи по специальности «Математика»

УДК 658.012

ФАКТОРИЗОВАННАЯ ФОРМА РЕКУРРЕНТНОГО АЛГОРИТМА ТЕКУЩЕГО РЕГРЕССИОННОГО АНАЛИЗА

ЯКУНИН А.В.

T

Предлагается численно устойчивая UDU -факторизованная схема рекуррентного алгоритма ТРА, основанного на МНК, для задачи параметрической идентификации линейной регрессионной модели.

1. Введение

Пусть идентифицируемый объект представлен линейно-параметризованной регрессионной моделью

Уn = cTxn n , (1)

где Xn eRN — вектор входных сигналов; yn eR1 — выходной сигнал; c eR N — вектор оцениваемых

параметров; £n eR1 — помеха измерения выходного сигнала; n = 0,1,2, ... — дискретное время. Если параметры объекта медленно изменяются во времени, то для отслеживания их динамики применяются алгоритмы текущего регрессионного анализа (ТРА) [1 - 4], основанные на учёте в рамках метода наименьших квадратов (МНК) лишь ограниченного числа L последних наблюдений.

Характеру процесса оценивания вполне соответствует предложенная в [4] рекуррентная схема алгоритма ТРА. Однако ей присуща потенциальная численная неустойчивость, обусловленная, во-первых, выделением в отдельный этап процедуры исключения устаревшей строки из матрицы наблюдений [5], во-вторых — наличием в вычислениях плохо обусловленных промежуточных, связанных с ковариационной, матриц [5, 6]. Кроме того, в приведенной в [4] формуле рекуррентного пересчёта оценки вектора параметров отброшены слагаемые, непосредственно учитывающие устаревшее наблюдение. Это ухудшает динамические свойства получаемой субоптимальной оценки.

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

на базе UDU т -факторизации ковариационной матрицы [6], обеспечивающей лучшую обусловленность промежуточных матриц.

2. Рекуррентный алгоритм ТРА на базе схемы замещения

Если в процессе идентификации на каждом n -м шаге используется лишь некоторое фиксированное

число L > N последних наблюдений, то текущая МНК-оценка вектора параметров имеет вид [2, 4]

Cn -Pn|L Xn|L Yn|L , (2)

где Xn|L = (Xn-L+i,Xn-L+2,---,Xn)T - матрица наблюдений входных сигналов размерности L х N ; -1

—ковариационная матрица раз-

мерн°сти N X N; Yn|L =(yn-L+1 Xn-L+2 ’• • •’ Уп)Т eR-L — вектор измерений выходного сигнала.

Пусть при поступлении нового наблюдения из рассмотрения одновременно исключается устаревшее. Для получения рекуррентных формул перехода от предыдущих оценок вектора параметров c n_1 и ковариационной матрицы Pn_1|L к текущим cn

и Pn|L применяются изложенные в [5] идеи процесса удаления и добавления строки в задаче МНК-оценивания.

Используя блочное представление матриц и вводя мнимую единицу i, i2 = -1, текущие оценки можно представить в виде

р.1т —

p n|L

X

T_!|LXn -1|L + (ixn-L^n):

(ixn-L M^n)1

3-1

(3)

Cn = Pn|L [xT_1L Mxn-L ^n J X

*( YHLMy»-lM'„ )T . (4)

Исходя из формулы обращения матричной суммы [7], после несложных, но громоздких преобразований соотношение (3) принимает рекуррентную форму:

Pn|L = Pn-l|L + (V8)Pn-1|L + а) xn-L xT-L _

-Pxn-L xT -Pxn xT-L -11 -y)xn xT

Pn-l|L

(5)

где

a_ xn Pn-l|L xn ; P_ xn-L Pn-l|L xn ; (6)

Y = xT-LPn-1|Lxn-L ; 8=1 + a-y-ay+p2 (7)

— скалярные коэффициенты, причём a > 0 и у > 0. Следует отметить, что наличие в (5) отрицательных членов, соответствующих устаревшему наблюдению, компенсируется, в отличие от [4], присутствием положительных слагаемых, отвечающих новому наблюдению.

С учётом полученных формул (5)-(7) выражение (4) также преобразуется к рекуррентному виду:

62

РИ, 2000, № 4

cn _ cn-l (V8) Sn_l|l |(^ + a)(yn-L cn-l x x xn-0 - P (yn - cT-l xn )] + (VS) Sn|L [ І1 “ y) x

:(yn -cT-1 xn) + p|Уп-L

- c

n -1 An-L

(8)

где

Sn-L|L = Pn-l|L Xn-L ; Sn|L = Pn-l|L xn (9)

— векторные коэффициенты усиления размерности N х 1.

3. Факторизованная схема рекуррентного оценивания

Однако нет необходимости находить указанную промежуточную факторизацию в явном виде. Более эффективна построенная по аналогии с [6] на основе (5)-(9) факторизованная рекуррентная схема ТРА в целом. Она имеет вид (n -l) - шаговой рекурсии:

fn-L _Un_i|l xn-L ; gn-L _Dn_i|l fn-L ; (11) fn =Un-l|L Xn ; gn = Dn-l|Lfn ; (12)

a1 _l + fnl gnl ; Pl _ fn-L)l gnl ; (13)

y= fn-L)lg(n -l)i; у 1 = 1 -у; (14)

Представив матрицу ковариаций Pn|L в факторизованной форме

Pn|L _U n|L D n|L U n|L j (10)

можно построить рекуррентную процедуру пересчёта матриц Un|L , Dn|L и коэффициентов усилеНия Sn-L|L , Sn|L •

Переписав соотношение (5) в виде

Un|L Dn|LUT|L = Un-1L {D„-1L + (V^[(l + «)* x(Dn-l|L U T_1L Xn-L XDn-1L U T_1L x"-L ) “

-p(Dn-l|LUn_,|L Xn-L)(Dn-lLUT_lLXn)' "

Dn-l|LUn_,|L X»)(Dn-1LUT_,LX»-L^

5l-al-y ; fn-L)l g(n-L)l; (15)

rnl = gnl ; D(n|L)ll = D(n-l|L)ll/8l ; (16)

для k = 2,3, ...,n последовательно вычисляются:

fnk gnk ; P fn —L)k gnk ; (17)

Y = fn-L)k g(n-L)k ; ^ = (ak-l/8k-l)fn-L)k ; (18) p = _(P k-1/8 k-l)fn _L)k; p = -(P k-l/8 k-Ofnk; (19)

®=-(y k-ll8k-l)fnk ; ak =ak_l +a ; (20)

Pk =Pk-l + P ; уk =yk-l-y ; (21)

8 k =8 k-1 _a k-1 a+ 2P k-1 P+Y k-1 У ; (22)

D(nf)kk = (8k-1 /8k)D(n-l|L)kk ; (23)

(1 -Д( D„ -l|L U n-lL X» )

fn-L)k g(n~L)k ; rnk - gnk ; (24)

для j = 1,2,k -1 последовательно вычисляются:

x(Dn-fL U T_1l xn ) j UML U(n|L)jk = U(n -l|L)jk fn-L)j +4 rnj +

и обозначив через Bn L матрицу в фигурных + Ёfn_L)j +® rnj ; (25)

скобках, а также положив f n n-L = Un-l|L Xn-L , fn-L)j ■_fn_L)j + U(n-l^jk g(n-L)k ; (26)

gn-L = D n-l|L fn -L , fn = U T_1L x n , gn = D n-l|L fn , rnj : = rnj + U(n-l|b)jkgnk . (27)

можно получить Bn|L = Dn-1|L + (V8)[(l + a) n gn-L gn-L По окончании рекурсии полагается a=an ; P = Pn ; у=уn ; 8=8n ; (28)

- Pgn-L gT - Pgn gn-L -(l“Г)gn g^ Un|LDn|LU„T|L = Un-lLBn|LUnT_lL

Если произвести факторизацию Bn|L = U^l D^l U^, то из единственности разложения следует

Un|L _ Un -l|L Un|L ; Dn|L _ Dn|L •

Sn-L|L rn -L ; Sn|L rn •

Далее по формуле (8) находится оценка cn .

(29)

В соотношениях (11)-(29) fn-L , gn-L , fn , gn , rn-L , rn — n -мерные рабочие векторы; a, p, у, X , ц,p,ю ,an,pn,уn,5n - скалярные рабочие величи-

ны.

РИ, 2000, № 4

63

4. Рекуррентное оценивание невязки

О точности получаемой оценки можно судить по величине квадрата длины невязки 2

Є

n

Yn|L Xn|L cn

В соответствии с рекуррен-

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

^n-l|L Xn-l|L cn (yn-L cn xn-l) +

+ (yn - cT xn )

с учётом соотношений (6), (7), (9) и

(Yn-l|L _ Xn-l|L cn-l) Xn-l|L = Yn-l|L _ Xn_^L

є n =

-T '’n-l|L

X± - Xn-l|L

-l

T

XT-lL Yn-lL

X

n -l|L

= Y

n -l|L

E - X

n-l|L

Xn-l|L Xn-l|L

-l

X

n -l|L

Xn-l|L -0 ; Pn-l|L [XT_lLXn-l|LJPn-l|L -Pn-l|L ,

где e -единичная матрица размерности L x L . Окончательная расчётная формула имеет вид

Щ _sn-l (VS) Iі + a)(yn-L cn-l xn-L^

_2P(yn-L _ cn-l xn-Lj(yn _ cn-l xn j _ ^ _ y) x

x (yn -cT-lxn) . (30)

Здесь є n _l — предыдущая оценка квадрата длины невязки.

5. Заключение

Разработанный алгоритм (8), (10)-(30) позволяет, аналогично [4], не только находить текущие оценки вектора параметров и ковариационной матрицы, но и рекуррентно пересчитывать длину невязки, дающей возможность судить о качестве процесса идентификации. Для повышения быстродействия можно, аналогично [4], вместо формулы (8) использовать соотношение

cn = cn-l +[(l-y)/s] Sn|L (yn - cT-l xn) , (31)

где непосредственно не учитывается отбрасываемое наблюдение. При этом справедлива следующая

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

Доказательство аналогично [4].

Рекуррентные процедуры оценивания чувствительны к выбору начальных оценок вектора параметров и ковариационной матрицы. Поэтому рекомендуется для вычисления хороших начальных

приближений вектора cn и матриц-факторов U n|L ,

D n і l использовать численно устойчивые процедуры, основанные на ортогональных преобразованиях с регуляризацией [8].

Предложенный алгоритм теоретически предполагает полноту ранга оценки ковариационной матрицы, что обеспечивается условием L > N . В дальнейшем предусмотрено модернизировать данный подход для случая L < N .

Литература: 1. Перельман И.И. Оперативная идентификация объектов управления. М.: Энергоиздат, 1982. 272 с. 2. Руденко О.Г., Штефан А., Хюбенталь Ф. Рекуррентный алгоритм МНК со скользящим окном при коррелированных помехах // Радиоэлектроника и информатика. 1998. №1. С. 79 - 81. 3. Косаревская Г.Н. Рекуррентный метод псевдообращения матриц и его применение в задаче идентификации // Вестник ЛГУ. Сер. 1. Математика. Механика. Астрономия. 1985. №22. С. 100 - 102. 4. Тимофеев В. А. Об одной модификации алгоритма текущего регрессионного анализа // Радиоэлектроника и информатика. 1999. №3. С. 30 - 32. 5. Лоусон Ч, Хенсон Р. Численное решение задач метода наименьших квадратов. М.: Наука, 1986. 232 с. 6. Bierman G.J. Factorization methods for discrete sequential estimation. N. Y.: Academic Press, 1977. 241 p. 7. Воеводин В.В., Кузнецов Ю.А. Матрицы и вычисления. М.: Наука, 1984. 320 с. 8. Регуляризованные UTDU -факторизации симметричных матриц и их применение / В.И. Мелешко, В.М. Задачин, Т.В. Ткаченко, И.Ф. Шматько. К.: Инт кибернетики, 1986. 28 с.

Поступила в редколлегию 27.06.2000

Рецензент: д-р техн. наук, проф. Филипенко И.Г.

Якунин Анатолий Викторович, канд. техн. наук, доцент кафедры высшей математики Харьковской государственной академии городского хозяйства. Научные интересы: численный анализ, распознавание образов, нестационарная гидродинамика. Увлечения: музыцирование, шахматы. Адрес: Украина, 61002, Харьков, ул. Революции, 12, тел. 45-90-30.

64

РИ, 2000, № 4