Filtration and diagnostics in discrete stochastic systems with jump parameters and multiplicative perturbations.pdf Задачи построения оценок и управлений для систем со случайными скачкообразными параметрами актуальны для различных реальных объектов. В качестве примера таких объектов можно рассмотреть, например, энергетические системы [1, 2], системы управления летательными аппаратами [3], системы связи [4, 5], задачи обнаружения неисправностей [6-8] и экономические системы [9, 10]. В работах [10-16] задачи фильтрации и управления рассматриваются для систем со случайными скачкообразными параметрами и аддитивными мультипликативными возмущениями. В настоящей работе рассмотрена задача фильтрации вектора состояния и идентификации состояния марковской цепи, которая входит в описание линейной стохастической системы с мультипликативными возмущениями. Решение получено с использованием принципа разделения, фильтрации Калмана и оценок неизвестного входного сигнала [17-22]. Предлагается выбрать матрицу коэффициентов передачи фильтра на основе минимизации суммы квадратичных форм ошибок оценивания. Приведен численный пример решения задачи фильтрации и идентификации для линейной стохастической системы с мультипликативными возмущениями и скачкообразными параметрами. 1. Постановка задачи Рассмотрим следующую линейную дискретную стохастическую систему с мультипликативными возмущениями и скачкообразным параметром: m x(k + 1) = (k)x(k) + Bi(k)u(k) + £ AS1 (k)x(k) 0s,Y(k)(k) + qy(k)(kX x(0) = ^ (1) s=1 где x(k) e Rn - вектор состояние системы, x0 - случайный вектор с известным математическим ожиданием и ковариацией x = М^}, N0. =М{(x0 - ^Хx0 - x0)T /у = у,.}, i = 1,r; Ay(kBy^), As,y(k) - известимте матрицы; u(k) e Rm - известный вектор; у = y(k) - марковская цепь с r состояниями (у1, y2,...,yr); qY(k)(k), 0sy(k) - гауссовские случайные векторы с характеристиками М{ qy№)(k)} = 0, М{0у№)(k)} = 0, М{qY(k)(k)qyT(k)(j) / уф = 1 (k),k < £ < j} = Ql(k)5^., M{0y(k)(k)0^)(j) /1 © = 1 (k),k < £ < j} = Hy№^ (M{} - математическое ожидание, T - символ матричного транспонирования, 5kj- - символ Кроне-кера). Вероятность состояний скачкообразного процесса pj(k) = P{y(k) = j}, j = 1,r, удовлетворяет уравнению Pj(k +1) = gp(k)р^, Pj(0) = Pj, j = й (2) i=i Здесь pt j - вероятность перехода из состояния , в состояние j за один шаг, pj - начальная вероятность j-го состояния. Вектор наблюдения: y(k) = Sy№) x(k) + vy(k )(k), (3) где vy(k) (k) - гауссовская случайная последовательность, не зависящая от qy(k) (k), с характеристиками: М{vy(k)(k)} = 0, М{vy(k)(k)vyT(k)(j)/у(S) = Y(k),k < S < j} = Vj(k)8y . Предполагается, что пара матриц Ay , Sy детектируема. По наблюдениям (3) требуется найти оценку параметра скачкообразного параметра y(k) (проблема диагностики) и оценку вектора состояния X(k), которую определим из условия минимума следующего критерия: J(Q,T,,) = М jg eT(k)He(k)/у(0) = у, j, (4) где H > 0 - весовая матрица, e(k) = x(k) - X(k), k e [Q, T]. 2. Синтез оптимального фильтра Нетрудно убедиться, что при ошибочном определении значения скачкообразного параметра y(k) модель (1) можно представить в эквивалентном виде, но с дополнительным вектором входа. Действительно, когда система (1) находится в j-м состоянии (у = уj) и это состояние ошибочно идентифицировано как ,-е (j Ф,), уравнение (1) представляется как модель с дополнительным входом: m x(k + 1) = A,x(k)+Bu(k)+f (k)+gAs ix(k) 0si (k) + qt (k), x(0) = xq, (5) s=1 где вектор дополнительного входа определяется по формуле mm f (k) = (Aj - A,)x(k) + (Bj - Bt )u(k) + g As jx(k)0s,j (k) - g Asix(k)0„ (k) + q] (k) - qt (k). (6) s=1 s=1 Здесь и далее введены обозначения для матриц AY(k), A5,y(k), BY(k^ Sy(k^ QY(k^ SY(k), VY(k) при Y(k) = Y, : A, A ., B , S., О , S,, V (, = й). Заметим также, так как j-е состояние нам неизвестно, то дополнительный вход f (k) также будет неизвестным. Поэтому для построения оценки вектора состояния воспользуемся рекуррентным оценивателем, совпадающим по структуре с фильтром Калмана и вычисляющим оценку состояния x(k) при наличии в модели объекта неизвестного входа [17-22]: x(k +1) = A,x(k) + f (k) + K(k)(y(k +1) -St(A,x(k) + f(k)), x(0) = xq, (6) где (k) (, = 1, r ) - матрица коэффициентов передачи фильтра. При этом мы воспользуемся принципом разделения, который означает, что сначала мы построим оценки вектора x(k) в предположении, что вектор f (k) известен, затем строятся оценки вектора f (k) в предположении, что оценка вектора состояния x(k) известна. Матрицы (k) определяются на основе минимизации критерия (4). Запишем критерий (4) в виде суммы: J (0,T, О = 2 tr Nt (k) H, (8) k=0 где tr - след квадратичной матрицы, матрица N, (k) = M{e(k )e(k )Т / у = yi} определяется из уравнения N(k +1) = (A - k(k)SA)(]>,jNj (k))(a - K(k)SA)T + к(k)viki (k)T + j=1 m r m /(\\\\ +(E - K (k )Si )[£ A, (X puNj. (k)) AT + 2 AAk) xT (k) A\\ + Qt ](E - Kt (k )S, )Т, N (0) = N0, (9) s=1 j=1 s=1 где Е - единичная матрица соответствующей размерности. Введем функцию Ляпунова следующего вида: W (k, Nt (k)) = tr Nt (k) H + tr i (t)L (t), (10) t=k где (t) = (E - K(k)S,)EtE - K(k)S,)T + Kt(t)VtKi(t)T + T,(t) (T,, (t) - некоторые положительно опре s.i / , s,i s=1 j=1 s=1 деленные матрицы, Ei =2Asi(2PijNj(k))АТ + 2Asixc(k)xT (k)ATsi + Q,). Здесь предполагается, что в (10) матрицы Li (t) положительно определены и удовлетворяют уравнению L (t) = (a - K (t)SA )Т (2 р,,jLj (t+1))(A - K (t)SA)+н + j=1 +(e - K (k)s, )t (2 aS,, (2 Pi,jLj (t+1))aS,,. )(e - K (k)s, ), L (T) = lT , (11) s=1 j=1 где LT - некоторая положительно определенная матрица. Суммируя по k = t,T -1 конечные разности функции W(k,N,. (k)) и учитывая формулу (11), получим T-1 T-1 2 a W (k, Nt (k)) =2 [W (k +1, Nt (k +1)) - W (k, Nt (k))] = k=t k=t = 2 tr[ Nt (k + 1)L, (k +1) - N, (k )L (k) - (k) L (k)]. k=t С другой стороны, это выражение может быть представлено следующим образом: T-1 2 AW (k, Nt (k)) = W (t +1, Nt (t +1)) - W (t, Nt (t)) +... + W (T, Nt (T)) k=t (12) (13) - W(T -1, N (T -1)) = tr N (T)Lt - tr N (t)L (t) - tr2 ^ (k)L (k). k=t Добавим к формуле (8) разность правых частей (12) и (13). В результате критерий (8) примет вид: j(0, t, i) = 2 tr ni (k)H - 2 tr n (k)l (k) + 2 tr [(at - K (k)s,a, )(2 p , N (k))(A, - K, (k)S,A, ) k=0 k=1 k=0 j=1 + (E - K (k)S,, )Hi (E - K (k)S,, )t,. + K (k)VK (k)T ]Li (k +1). (14) Применяя правила матричного дифференцирования функции tr [23], вычислим производные J (0,T, 0 {222tr N (k) н - 2tr N (k) L (k) + ж, (k) ж, (k) tf +2tr[(A - K{k)SA)(ZPijNj (k))(a, - к(k)SA,)T + (E - K(k)S,)S, (E - к(k)S,)T + k =0 j=1 + +к, (k)vlKl (k)T (k +1)} = X 2[-L, (k+1)Д. (X p,N (k)) aJsJ + L (k + 1)K (k)SAx r *=0 ^ (15) x(XPuN-(k))aJsJ -L(k+1)stEST + L(k+1)K(k)stEST + L(k+1)k,(k)v]. j=1 Приравнивая эту производную к нулю, получаем формулу для определения матриц K. (k) (t = 1, r): K(k) = (A,(]Tp,j.Nj.(k))AjSj + E.Sj)[S . A(Zp-Nj(k))AjSj + St E,Sj + Vt]-1. (16) j=1 j=1 Учитывая (12) и (13), можно показать, что конечная разность функции Ляпунова (10) определяется по формуле AW (k, Nt(k)) = W (k +1, Nt (k +1)) - W (k, Nt (k)) = = trN(k +1)H -trN(k)H - tr[(E -K(k)S,)E ,(E -K,(k)S,)T + K,(k)VtKt(k)T+^,(t)]L(k). ^ В силу того, что решения матричных разностных уравнений (9), (11) положительно определены, функция Ляпунова (10) будет положительной. Устойчивость по Ляпунову динамики фильтра гарантируется тем, что матрицу (t) всегда можно задать так, чтобы выражение (17) стало отрицательным. В стационарном случае матрицы переноса Kt будут постоянными и определяются из следующих матричных алгебраических уравнений: N = (A - KS A )(Z PuN)(a, - KSA )T + (18) -=1 m +(E - KSt )[Z As,t(Z P-Nj) AT +Z А^, , XXj Ajt + a ](E - KS, )T + KyK m K , = (At (Z Pt,jNj )Aj Sj + [Z As, t(Z Pt,jNj) AJ, + X A. ,xxT Aj, + Q, ]Sj) x j=1 s=1 j=1 s=1 r m r m х^,Аг(ZPtJNj)ATSj + S, [X4,t(XJ)Asj, +ZAs,XXj4* + Qt]Sj + Vt]-1. j=1 s=1 j=1 s=1 s,t V^r t,j j s,t ^^ s, t s,t s=1 j=1 s=1 (19) Отметим, что если существуют положительно определенные решения N, (t = 1, r) матричного m r m уравнения (18), то из условия (E - KS,)[£ АЛХPuN)Aj +Z 4,^Aj + Q,](E - KS,)j + KyK- > 0 s=1 j=1 s=1 следует справедливость теоремы 1.6 [24], а это означает устойчивость стационарного фильтра: X(k +1) = AX(k) + /Д) + K, (y(k +1) - S, (4X(k) + f,(k)), X(0) = X, (20) В качестве алгоритма оценки неизвестного входа мы будем использовать МНК-оценки; в этом случае оценка может быть построена на основе минимизации дополнительного критерия [17] в предположении, что значение скачкообразного параметра известно (здесь предполагается, что у = у,): k+11 2 2 JКLMS> (/ (k)) = X {||y(t) - St (4 X(t -1) + Btu(t -1) + / (t -1))||W + || / (t -1) ||W }, (21) где W, W - положительно определенные весовые матрицы. Минимизируя (21), мы получаем оценки неизвестного входа: /(k) = (S,T WSt + W)1 Sj W[y(k +1) - S, (At X(k) + Bt u(k))]. (22) 3. Оценка скачкообразного параметра Из (6) видно, что дополнительный вход в модели (5) будет отсутствовать, если состояние скачкообразного параметра идентифицировано правильно, т.е. j = t. Поэтому предлагается алгоритм оценивания параметра у(k) строить на основе определения номера t по оценке (22), для которой евклидова норма минимальна. Отметим, что для повышения точности определения значения оценки скачкообразного параметра y(k) предлагается использовать предварительное экспоненциальное сглаживание значения нормы Значение параметра определяется по следующему алгоритму: 1) вычисляется норма оценки неизвестного вектора y(i,k) = || f (k) || для всех у = у. (i = 1,r); 2) выполняется экспоненциальное сглаживание для значений y(i, k) у^ (i,k +1) = a^(i,k) + (1 -a)y^ (i, k), i = 1, r, (23) где a - коэффициент сглаживания; 3) для каждого момента времени k определяется значение i, для которого сглаженное значение нормы у (i, k) минимально (i = argmin{y (i, k)}). В качестве оценки y(k) выбирается уt, у которого значение индекса i определено в п. 3 алгоритма. f (k ) 4. Результаты моделирования Рассмотрим задачу моделирования алгоритма фильтрации и диагностики для дискретной системы с мультипликативными возмущениями и тремя состояниями цепи y(k) (у1 = 1, у2 = 2, у3 = 3) со следующей матрицей вероятностей перехода: (0,4 0,2 0,4 1 [ Ри ] = 0,3 0,5 0,2 v0,3 0,3 0,4, Моделирование выполнено на интервале времени k е [0, 400] для следующих исходных данных: ( 0,85 0,1 Л v-0,05 0,94, ( 0,6 0,051 ( 0,7 0,121 (1 01 v0 1, A = A = A3 = B1 = B2 = B3 = v-0,02 0,45, v-0,04 0,62, (0,31 v 0,3 (0,005 0 1 (0,001 0 1 0 0,001 (0,001 0 1 0 0, 0005 u(k) = Q1 = Q2 = Q3 = v 0 0,005, s = s2 = s3 = (1 1), v = v2 = v3 = 0,01. (0,1 0 1 W = 1, W2 = 0 0,1, Оценки фильтрации рассчитываются по формулам (18)-(20), (22). Параметр y(k) оценивался с помощью алгоритма, описанного в разд. 3. Результаты моделирования приведены на рис. 1-4, которые иллюстрируют качество оценивания и идентификации. Рис. 2. Сглаженные величины норм y(i,k) (i = 1,3 ) Fig. 2. Smoothed values of norms y(i,k) (i = 1,3 ) Рис. 1. Значения параметра y(k) и оценки y(k) Fig. 1. Parameter value y(k) and estimate y(k) На рис. 1 видно, что ошибки идентификации возникают в момент изменения значения параметра перехода. На рис. 2 показано, что минимальному значению нормы у(/, k) соответствует номер i, который принимается за оценку скачкообразного параметра y(k). Рис. 3. Графики компоненты вектора состояния x1(k) и ее оценки x1(k) Fig. 3. Graphs of the component of state vector x1(k) and its estimate x1 (k) Рис. 4. Графики компоненты вектора состояния x2(k) и ее оценки x2(k) Fig. 4. Graphs of the component of state vector x2(k) and its estimate x2 (k) С помощью метода статистического моделирования (для 100 реализаций) получили, что процент ошибочных оценок значений параметров y(k) составляет 5,75%. Заключение Получено решение задачи синтеза алгоритма фильтрации вектора состояния и идентификации значения скачкообразного параметра, входящего в описание линейной дискретной системы с аддитивными и мультипликативными возмущениями. Задача решена посредством интерпретации модели в виде модели системы с неизвестным входным вектором, который появляется при ошибке идентификации скачкообразного параметра. Значение оценки скачкообразного параметра предлагается определять из условия минимума сглаженной нормы оценки вектора неизвестного входа.

Скачать электронную версию публикации