CC BY
77
ВЕСТНИК ТОМСКОГО ГОСУДАРСТВЕННОГО УНИВЕРСИТЕТА
2011 Управление, вычислительная техника и информатика № 3(16)
УДК 517.511
В.И. Смагин, С.В. Смагин
ФИЛЬТРАЦИЯ В ЛИНЕЙНЫХ ДИСКРЕТНЫХ НЕСТАЦИОНАРНЫХ СИСТЕМАХ С НЕИЗВЕСТНЫМИ ВОЗМУЩЕНИЯМИ
Рассматривается алгоритм синтеза оптимального фильтра, определяющего оценку вектора состояния дискретной линейной нестационарной динамической системы с аддитивными возмущениями, содержащими неизвестную постоянную составляющую. Приводятся результаты вычислительного эксперимента.
Ключевые слова: линейные дискретные нестационарные системы, фильтр Калмана, неизвестные возмущения.
В работах многих авторов большое внимание уделяется разработке алгоритмов калмановской фильтрации для класса систем с неизвестными аддитивными возмущениями и параметрами, которые могут использоваться в качестве моделей реальных физических систем, моделей объектов с неизвестными сбоями.
Известные методы вычисления оценок вектора состояния базируются на алгоритмах, использующих оценки неизвестного возмущения [1 - 11]. В работах [1, 2] рассматриваются алгоритмы расширения пространства состояний (к основной модели объекта добавляется модель ненаблюдаемого возмущения) и алгоритм двухэтапной фильтрации, уменьшающий вычислительные затраты за счет декомпозиции задачи. В работах [3 - 11] изучены алгоритмы рекуррентной оптимальной фильтрации, использующие оценки неизвестного возмущения, имеющие достаточно жесткие условия их разрешимости.
В настоящей работе для дискретного нестационарного объекта с неизвестной постоянной составляющей возмущений предлагается метод оптимальной фильтрации, не использующий оценки неизвестного возмущения. Метод базируется на преобразовании модели и сведении к задаче линейной калмановской фильтрации [12 - 14]. В настоящей статье обобщаются результаты [15] на случай решения задачи для нестационарного дискретного объекта.
1. Постановка задачи
Рассматривается дискретная система, которая описывается следующими разностными уравнениями:
x(k +1) = A(k)x(k) + f + q(k), x(0) = x0 , (1)
где x(k) e Rn - вектор состояния; A(k) - nxn-матрица; f - неизвестный постоянный вектор; q(k) - белая гауссовская случайная последовательность с характери-
стиками
M {q(k)} = 0 , M{q(k)qT ( j)} = Q(k)bk ] . (2)
Канал наблюдений имеет вид
y(k ) = S (k ) x(k ) + v(k ), (3)
y(k) e R1 - вектор измерений; S(k) - матрица размерности l x n ; v(k) - белая гаус-
совская случайная последовательность ошибок измерений, с характеристиками:
М{у(к)} = 0 , М{д(к)ут (])} = 0, М{у(к)ут (у)} = V(к)81} ; (4)
для матриц (£(к), А(к)) выполняются условия наблюдаемости. Вектор х0 является случайным и не зависит от от процессов д(к) и у(к), при этом
М{х(0)} = хо, М {(х(0) - Хо )(х(0) - Хо)т } = Ро .
Для системы (1) и канала наблюдений (3) требуется синтезировать фильтр, вычисляющий оценку вектора состояния, не использующий оценки неизвестной постоянной составляющей возмущений.
2. Синтез фильтра
Преобразуем дискретную систему (1). Исключаем постоянную составляющую возмущений / из описания объекта посредством вычитания из уравнения (1) такого же уравнения, но со сдвигом на один такт:
х(к) = А (к -1) х(к -1) + / + q(k -1). (5)
В результате получаем следующее уравнение:
х(к +1) = (А (к) + Еп) х(к) - А (к -1) х(к -1) + q(k) - q( к -1). (6)
Расширим пространство состояний системы путем добавления к уравнению (6) тождества х(к) = х(к). Обозначим
X (к) = ( хГ„ ) • «к) = ( q'k)- ■>). ()
Систему (1) представим в векторно-матричной форме
X(к +1) = А (к)X (к) + д (к), X (0) = Х0, (8)
где А (к) - 2п х 2п -матрица имеет следующую блочную структуру:
-к) = ( А<кЕ+ Еп -А<0 - 0 ^. (
Случайный вектор X0 = (х^ х-1 )т имеет следующие характеристики:
М{X(0)} = X0, М {(X0 -X0)^0 -X0)т} = Р0, (10)
где X0 = (х0т х- )т . Отметим, что здесь дополнительно вводится п-мерный вектор х-1, который является независимым от д(к) и у(к), а характеристики (10) могут быть получены по априорной информации об объекте (1).
Отметим, что в рассмотренной модели (8) процесс д (к) не является белой гауссовской последовательностью, процессы д (к) и д (к -1) будут коррелированны:
Q(k), если у = к,
Q(k -1), если ] = к -1, (11)
0, если 0 < ] < к -1,
М{д (к) д т (у)} =
где Q{k) = |'«к> + «к-1) 0), Q<k-1) =(-«0-■> ¡¡). (12)
Представим канал наблюдений для расширенной системы (8) в виде
у(к) = 5 (к) X (к) + v(k), (13)
где 5 (к) = (5(к) 0), v(k) - случайная последовательность ошибок измерений с характеристиками (4).
В качестве уравнения для вычисления оценки вектора состояния расширенной системы выберем уравнение, по своей структуре совпадающее с фильтром Кал-мана:
Х(к +1) = Л(к )Х(к) + К (к )(y(k +1) - 5 (к +1) Л (к )Х(к)), Х(0) = Х0. (14)
Учитывая (8) и (14), получим следующее уравнение для ошибки е(к) = Х(к) - X(к):
е(к + 1) = (Л(к) - К (к) 5 (к +1) Л(к ))е(к) + К (к )м(к + 1) + (К (к) 5 (к + 1) - Е2п )д (к). (15)
В силу (11) и (15), матрица Р (к) = М{е(к )ет (к)} определится из следующего разностного уравнения:
Р(к +1) = (Л(к) - К (к)5(к +1)Л (к))Р(к)(Л(к) - К (к)5(к +1)Л (к))т +
+(К (к) 5 (к +1) - Е2п Йк)(К (к )5 (к +1) - Е2п )т + К (к) V (к +1) Кт (к) +
+(Л(к) - К (к )5 (к +1) Л(к))(К (к -1)5 (к) - Е2п) х
х0(к -1)( К (к) 5 (к +1) - Е2 п )т + (К (к )5 (к +1) - Е2п) х
х0(к - 1)(К (к -1) 5 (к) - Е2п )т (Л(к) - К (к)5 (к + 1)Л(к))т , Р(0) = Р0. (16)
Оптимизируемый критерий зададим в виде
3 (к +1) = ЪР (к +1). (17)
Оптимальные коэффициенты передачи фильтра К(к) определяются из условия
3 (к +1)
dK (к)
= 0. (18)
Учитывая (17) и правую часть уравнения (16), применяя правила матричного дифференцирования следа от матрицы [16], получим из условия (18) уравнение для определения матрицы К(к):
- Л (к) Р (к) Л(к )т 5 (к + 1)т + К (к) 5 (к +1) Л (к) Р (к) Л(к )т 5 (к + 1)т +
+К (к) 5 (к + 1)^(к )5 (к )т - &(к) 5 (к + 1)т - К (к) 5 (к + 1)0(к -1) х х5 (к )т К (к - 1)т Л(к )т 5 (к + 1)т + К (к) 5 (к + 1)0(к -1) Л(к )т 5 (к + 1)т -
-К (к) 5 (к +1) Л(к) К (к -1) 5 (к )0(к -1)5 (к + 1)т +
+К (к )5 (к +1) Л(к )0(к -1) 5 (к + 1)т + 0(к -1)5 (к )т К (к - 1)т х хЛ(к)т 5(к + 1)т - 0(к -1)Л(к)т 5(к + 1)т -Л(к)0(к -1)5(к + 1)т +
+Л( к) К (к -1) 5 (к)0(к -1) 5 (к + 1)т + К (к V (к +1) = 0. (19)
Решение последнего уравнения относительно К(к) дает следующий результат:
К (к) = Р(к)5(к + 1)т (5(к +1)Р(к)5(к + 1)т + V(к +1))-1, (20)
где P(к) = Л (к)P(к)Л(к)т + Q(k - 1)(E2n - S(к)т K(к - 1)т )Л(к)т +
+A(k)(Eln - K(к -1)5(к))Q(k -1) + Q(k). (21)
Отметим, что для вычисления коэффициентов передачи (20), в силу (21), необходимо задать начальные значения коэффициентов K(-1).
Подставив в уравнение (16) выражение для оптимального коэффициента передачи (20), получим уравнение
P(к +1) = (E2n - K (к)S (к +1))P(к), P(0) = Р0. (22)
Основной результат сформулируем в виде теоремы, учитывая симметричность и блочное представление матриц P(к) и P(к):
P(к) = f p (к) (к) 1, P(k) = f p1(к) p2T (к) 1, (23)
IР 2 (к) p з(к)) У Р2(к) Рз(к))
блочные структуры матриц Л(к), Q(k), Q(k), S(к) и представление матрицы K(к) в виде
k (к >=( K%). <24)
Теорема. Пусть процесс с неизвестным постоянным возмущением определяется уравнениями (1) и канал наблюдений имеет вид (3). Тогда оптимальный алгоритм фильтрации определится следующими разностными уравнениями: x(k +1) = (A (k) + En) x(k) - Л (к -1) x(k -1) + K1 (к)(y(k +1) -
- S (к +1)[( Л (к) + En) x(k) - Л (к -1) x(k -1)] (25)
с начальными условиями
x(0) = x0, x(1) = M{x(1)} = x . (26)
Матрица Kx (k) в (25) определяется по формуле
K (к) = рх (к) S (к + 1)т (S (к +1) Р (к) S (к + 1)т + V (к +1))-1, (27)
где матрица р (к) в^гчисляется из системы уравнений
Р(к) = (Л (к) + En)р1(к)(Л (к) + En )т - Л(к -1) Р2(к)(Л (к) + En )т -
-(Л( к) + En) рТ (к) Л( к - 1)т + Л (к -1) р3(к) Л (к - 1)т + Q( к -1) S (к )т Kj(k - 1)т х х( Л (к) + En )т - Q(k -1) S (к )т K2 (к - 1)т Лт (к -1) +
+(Л( к) + En) Kj (k -1) S (k) Q( k -1) - Л( k -1) K2 (k -1) S (k) х xQ (k -1) - (Л(к) + En )Q(k -1) - Q(k -1)( Л(к) + En )T + Q(k) + Q(k -1),
Р2 (к) = #(k)(Л(к) + En )T - p2 (к)Л(к - 1)T +
+K1 (k - 1)S(k)Q(k -1) - Q(k -1), ^3 (k) = p1 (k),
Р1 (k +1) = (En - K (k)S(k +1))p (к), Р1 (0) = Р10,
Р2 (k +1) = - K 2 (k) S (k +1) P (k) + p 2 (к), Р2 (0) = Р20,
Р3 (k +1) = -K2 (k)S(k +1) p2 (к) + Р3 (к), Р3 (0) = p^ ,
K2 (k) = p 2 (k )S (k + 1)T (S (k +1) p (k) S (к + 1)T + V (к +1))-1. (28)
В (28) начальные условия р10, р2 0, р3 0, являются соответствующими блоками матрицы Р0. Отметим, что для выполнения расчетов в (28) необходимо задать начальные условия для КД-1) и К2(-1).
Замечание. Управляемый объект
х(к +1) = А(к)х(к) + В(к)и(к) + / + д(к), х(0) = х0, (29)
при исключении неизвестного постоянного возмущения / объекта, необходимо преобразовать к виду, который будет отличаться от (8) одним слагаемым:
X (к +1) = А (к) X (к) + В(к )(и(к) - и(к -1) + д (к), X (0) = Х0, (30)
где матрица А(к) приведена в формуле (9), д (к) имеет характеристики (11), (12). В (30) матрица В (к) имеет вид
В(к) =(В0к)). (31)
Тогда уравнения фильтра будут следующими:
Х(к +1) = (А (к) + Еп) Х(к) - А (к -1) Х(к -1) + В(к )(и(к) - и (к -1)) + К1 (к)(у(к +1) --Б(к +1)[(А(к) + Еп)Х(к) - А(к -1)Х(к -1) + В(к)(и(к) -и(к -1))], (32)
с начальными условиями (26), а матрица К1(к) определяется в соответствии с (27) и (28).
3. Результаты вычислительного эксперимента
Рассмотрим применение алгоритма фильтрации для модели второго порядка вида (1) и канала наблюдений (3) со следующими значениями параметров:
( 0 1 А _ (0,01 0 А ТЛ
() = ^0,05 0,925 + 0,Ыи(0,01к)) ’ ® = [ 0 0,02} ’ = , ’
х = (1 1); Х0 =(иР0 =(100 100} (
Вычисление оценок вектора Х(к) можно выполнить, используя двухэтапный алгоритм фильтрации [2]. Модель измерений в этом случае с учетом (1) представляется в виде
у(к +1) = 5Х(к +1) + у(к +1) = £А(к)Х(к) + Б/ + 5д(к) + у(к +1). (34)
Рекуррентные уравнения оценивания неизвестного вектора / имеют вид /(к +1) = /(к) + К/ (к)(у(к +1) - £А(к)Х(к) - Б/(к)), Д0) = /0,
Кг (к) = Рг (к) Бт (БРГ (к) Бт + ^т + V )-1,
Р/ (к +1) = (Е2 - К/ (к)Б)Р/ (к), Р/ (0) = Р/0, (35)
где М{/} = /0, М{(/ - /0)(/ - /0)т } = Р/0. (36)
Оценка вектора состояния для объекта с неизвестным постоянным входом задается уравнением:
х(к +1) = Л(к)х(к) + /(к) + Кх (к)(у(к +1) - БЛ(к)х(к) - Б/(к)) , (37)
где матрица Кх (к) определяет коэффициенты передачи фильтра Калмана. При моделировании используем
(01 Р =Г1,0 01 ,0 ] , Л { 0 1,0 ].
/ =
(38)
Применение расширенного фильтра Калмана [1] для данного примера (в этом случае уравнение (1) расширяется путем добавления уравнения /(к+1) = /(к)) приводит к необходимости построения фильтра Калмана для дискретной системы со следующими матрицами динамики, канала наблюдений и интенсивностей аддитивных возмущений:
( Л (к) К21 ( (б 01
0 Е2
0 0
(39)
Использование в данном примере методов, описанных в работах [4, 5, 8, 9], невозможно в силу невыполнения условий существования оптимальных оценок неизвестного входного вектора [4, 5, 8, 9]:
п > т и I > т . (40)
В [4, 5, 8, 9] неизвестное возмущение определяется в виде / = Ой, где й - неизвестный т-мерный вектор, О - п х т -известная матрица. В рассмотренном примере О = Е2, п = 2, т = 2, I = 1, а это означает, что условия (40) не выполняются.
Применение алгоритма фильтрации исследовалось также для неизвестного переменного возмущения с тремя возможными значениями компонент вектора /:
1, если 0 < к < 9,
/1(к) = /2(к) = < -1, если 9 < к < 25,
1, если 25 < к < 50.
На рис. 1 приведены реализации процессов и их оценок для трех сравниваемых фильтров. Отметим, что при реализации алгоритма фильтрации (25), начальные значения К1(-1) и К2(-1) задавались нулевые.
-10
0
10
20
30
40
-10
к 0
10
20
30
40
к
Рис. 1. Реализации процессов и оценок (1 - реализация х(к); 2 - оценка, построенная по алгоритму (25); 3 - оценка, построенная по двухэтапному алгоритму; 4 - оценка для расширенного фильтра Калмана)
На рис. 2 приведены ошибки оценивания компонент вектора состояния.
Рис. 2. Графики ошибок фильтрации (1 - ошибка для оценки, построенной по алгоритму (25); 2 - ошибка для оценки, построенной по двухэтапному алгоритму; 3 - ошибка для расширенного фильтра Калмана)
Как видно из рисунков для рассмотренного примера, качество оценок, полученных с помощью фильтра (25), лучше, чем для двухэтапного алгоритма фильтрации и расширенного фильтра Калмана, использующих оценки неизвестного возмущения. Отметим также, что для алгоритма фильтрации (25) не нужно задавать априорную информацию о характеристиках распределения начальных значений 7 и Р/0 .
Ниже, в таблице, приведены средние значения среднеквадратических ошибок оценивания для трех рассматриваемых методов, рассчитанных по 50 реализациям. Как видно из таблицы, предложенный метод фильтрации (25) обеспечивает среднюю ошибку в 3 - 4 раза меньшую, чем другие методы.
Средние значения среднеквадратических ошибок для компонент вектора состояния
Алгоритм (25) Двухэтапный алгоритм Расширенный фильтр Калмана
е1>Ср = 0,0912 е1,ср = 0,3128 Єі,ср = 0,4103
Є2,ср = 0,0945 е2,ср = 0,2917 е2,ср = 0,4296
Заключение
Разработан алгоритм синтеза дискретного оптимального нестационарного фильтра для объекта, возмущения которого содержат неизвестную постоянную составляющую. Алгоритм построен на основе расширения пространства состояния и исключения из модели неизвестной составляющей. В отличие от классического фильтра Калмана, предложенный фильтр использует рекуррентные оценки, построенные на двух предыдущих тактах. Как показали результаты вычислительного эксперимента, алгоритм может быть применен для кусочно-постоянной неизвестной аддитивной составляющей возмущений.
ЛИТЕРАТУРА
1. Astrom K., EykhoffP. System identification. A survey // Automatica. 1971. V. 7. P. 123-162.
2. FriedlandB. Treatment of bias in recursive filtering // IEEE Trans. on Automat. Contr. 1969. V. AC-14. P. 359-367.
3. Chen J., Patton R. J. Optimal filtering and robust fault diagnosis of stochastic systems with unknown disturbances // IEE Proc. Control Theory Appl. 1996. V. 143. P. 31-36.
4. Darouach M., Zasadzinski M. Unbiased minimum variance estimation for systems with unknown exogenous inputs // Automatica. 1997. V. 33. P. 717-719.
5. Darouach M., Zasadzinski M., Xu S. J. Full-order observers for linear systems with unknown inputs // IEEE Trans. on Automat. Contr. 1999. V. AC-39. P. 606.
6. Gillijns S., Moor B. Unbiased minimum-variance input and state estimation for linear discrete-time systems // Automatica. 2007. V. 43. P. 111-116.
7. Hou M., Patton R. Optimal filtering for systems with unknown inputs // IEEE Trans. on Automat. Contr. 1998. V. AC-43. P. 445-449.
8. Hsieh C.-S. A unified solution to unbiased minimum-variance estimation for systems with unknown inputs // Proc.17th World Congress The International Federation of Automatic Control. Seoul. Korea. July 6 - 11, 2008. P. 14502-14509.
9. Hsieh C.-S. Robust two-stage Kalman filters for systems with unknown inputs // IEEE Trans. on Automat. Contr. 2000. V. AC-45. P. 2374-2378.
10. Hsieh C.-S. Extension of the optimal unbiased minimum-variance filter for systems with unknown inputs // Proc. 15th IEEE International Workshop on Nonlinear Dynamics of Electronic Systems. Tokushima. Japan. 2007. P. 217-220.
11. Hsieh C.-S. Robust parameterized minimum variance filtering for uncertain systems with unknown inputs // Proc. American control conference. New York. 2007. P. 5118-5123.
12. Kalman R.E., Busy R. A new results in linear filtering and prediction theory // Trans. ASME J. Basic Engr. 1961. V. 83. P. 95-108.
13. Браммер К., ЗиффлингГ. Фильтр Калмана - Бьюси. М.: Наука, 1972. 200 с.
14. Пугачев В.С., Синицин И.Н. Стохастические дифференциальные уравнения М.: Наука, 1990. 630 с.
15. Смагин С.В. Фильтрация в линейных дискретных системах с неизвестными возмущениями // Автометрия. 2009. Т. 45. № 6. C. 29-37.
16. Амосов А.А., Колпаков В.В. Скалярно-матричное дифференцирование и его применение к конструктивным задачам теории связи // Проблемы передачи информации. 1972. № 1. С. 3-15.
Смагин Валерий Иванович
Смагин Сергей Валерьевич
Томский государственный университет
E-mail: vsm@mail.tsu.ru; ssv@fpmk.tsu.ru Поступила в редакцию 6 декабря 2010 г.
