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

УДК 519.248:681.5.001.3

ОБОБЩЕННЫЕ ВЕРОЯТНОСТНЫЕ КРИТЕРИИ В ЗАДАЧЕ НЕЛИНЕЙНОЙ ПАРАМЕТРИЧЕСКОЙ ИДЕНТИФИКАЦИИ

C.B. Соколов, П.А. Кучеренко

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

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

ВВЕДЕНИЕ

Существующие методы решения задачи стохастической параметрической идентификации требуют, как правило, для своей удовлетворительной реализации принятия ряда крайне упрощающих ограничений — линейности модели измерителя относительно параметров, необходимости нормального вида распределения (с нулевыми математическими ожиданиями и конечными дисперсиями) аддитивных помех наблюдаемых сигналов и др. В большинстве реальных ситуаций это оказывает значительное отрицательное влияние на качественные характеристики процедуры идентификации и, как следствие, снижает потенциально возможную точность получаемых оценок параметров [1—8]. Анализ литературы по данной тематике показал, что работы в этой области в основном посвящены вопросам, связанным с идентификацией параметров модели наблюдаемого стохастического объекта, в то время как вопросы, касающиеся определения неизвестных параметров самого наблюдателя (измерителя параметров состояния объекта), практически не освещены. Вместе с тем, задачи подобного рода очень часто возникают в самых различных областях телекоммуникации и связи, радионавигации, метрологии и пр. К числу наиболее распространенных из них можно отнести задачи определения неизвестных характеристик тракта передачи (времени распространения сиг-

налов, параметров среды передачи и др.), определения неизвестных или «шумящих» коэффициентов усиления (или других характеристик) приемника и пр.

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

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

Вопросы идентификации по сложным статистическим критериям уже рассматривались в работах различных авторов [7, 8], но принятые в них критерии, основанные на минимизации средне-квадратической ошибки оценивания и определении апостериорного математического ожидания оцениваемых параметров (в том числе, с использованием статистической и дисперсионной линеаризации) или максимизации апостериорной плотности вероятности (АПВ) параметров наблюдаемого объекта, как это будет видно из последующих

24

CONTROL SCIENCES № 4 • 2008

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

1. ПОСТАНОВКА ЗАДАЧИ НЕЛИНЕЙНОЙ ПАРАМЕТРИЧЕСКОЙ ИДЕНТИФИКАЦИИ

Пусть Л-мерный вектор состояния нелинейного дискретного объекта хк задан в самом общем случае разностным уравнением

хк = ^ - 1, п), (1)

где — известная нелинейная вектор-функция с дифференцируемыми компонентами, для которых существуют обратные функции; хк _ 1 — Л-мер-ный вектор переменных состояния на (к — 1)-м шаге времени; п — Л-мерный вектор шума с известной Л-мерной невырожденной функцией плотности распределения вероятности #(п).

Вектор наблюдения размерности М для к-го момента времени zk описывается следующим уравнением (в общем случае также нелинейным):

% = Х(С^ % ^ (2)

где с(...) — известная нелинейная вектор-функция наблюдения с дифференцируемыми компонентами, для которых существуют обратные функции; Ск — вектор (или матрица) параметров наблюдателя соответствующей размерности, в общем случае нестационарный; w — М-мерный вектор шума с известной М-мерной невырожденной функцией плотности распределения вероятности

Для сокращения дальнейшей записи набор векторов сигналов наблюдения zг, i = 1, ..., к, обозна-к

чим через z1.

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

вероятности переменных состояния р(хк| zk; Ск) в общем случае нелинейна:

J = 1 Ф(р(хК; Ск))^к = Ж(Ск),

(3)

где Ф — известная нелинейная аналитическая функция, X — заданная область пространства состояний; р(хк| zk; Ск) — АПВ вектора переменных

состояния объекта, W — известная нелинейная функция.

Различные вариации вида критериальной функции Ф позволяют охватить достаточно широкий класс критериев оптимальности [9]:

— для критерия минимума отклонения АПВ от

заданной функции 5: Ф(р) = (р — 5) , Ф(р) = —р1п(5/р) (критерий Кульбака) и др.;

— для критерия максимума (минимума) вероятности существования параметров состояния в некоторой области пространства состояний: Ф(р) = ±р;

— для критерия максимума информации о параметрах состояния: Ф(р) = —рГГ!(кри-

^ дхк' ^дхк'

терий Фишера), Ф(р) = —р ^р (критерий Шеннона) и др.

Очевидно, что задача в данной постановке сводится к отысканию текущей апостериорной плотности р(хк| z1; Ск) и последующему определению текущего значения искомого параметра Ск, удовлетворяющего выбранному критерию оптимальности J.

2. СИНТЕЗ АЛГОРИТМА НЕЛИНЕЙНОЙ ПАРАМЕТРИЧЕСКОЙ ИДЕНТИФИКАЦИИ

Известно [10], что многомерная апостериорная плотность вероятности вектора состояния х для

к-го момента времени р(хк| zk; Ск) определяется выражением

+ да + да

Р(хк14; Ск) = Щ^) 1... 1 Р(хк _ 11 zk ^1;

—да —да

Ск _ 1МХК _ 1)йХк _ ЛМ,

+ да +да / + да + да

МС) = 1... 1 [ 1... 1 р(Хк _ 11 zí-1;

—да —да —да —да

Ск _ 1^Р(хк|хк _ 1)йХк _ lP(zk|xk) I ^ (4)

где р(хк - 11 zk 1; Ск _ 1) — определенная на (к — 1) -м шаге апостериорная плотность вероятности вектора состояния объекта; р(хк|хк _ 1) — определяемая на текущем шаге алгоритма Л-мерная условная плотность вероятности вектора хк; p(zk|xk) — определяемая на текущем шаге алгоритма функция правдоподобия для многомерного наблюдения.

Многомерная условная плотность р(хк|хк _ 1) может быть получена из исходного уравнения объекта (1) при известном виде плотности распределения вероятности значений шума #(п) (в предпо-

ложении их взаимной статистической независимости) [11]:

Р(хк|хк - 1) = ^У^ Хк - 1))|YУ|,

дУ( 1 ) ( Хк. Хк - 1 ) д/( 1 ) ( Хк. Хк - 1 )

где

где Yy =

дх,

к( 1)

дх

к (Л)

дУ( л-) (

Хк> Хк - 1)

дУ ( л) (

Хк> Хк - 1)

дх

к( 1)

дх

к (Л)

У(хк, Хк - 1) =

У( 1)(Хк, Хк - 1)

7( Л)(

Хк> Хк - Ь

где У^ — якобиан преобразования от вектора переменных п к вектору переменных хк, зависимость которых в общем случае нелинейная и определяется выражением (1); У(г)(хк, хк - 1), i = 1, ..., N — полученные в результате обратного преобразования соответствующих компонент однозначно определенные функции; х^, i = 1, ..., N — компоненты вектора состояния объекта хк.

Аналогичным образом из уравнения (2) можно определить входящую в выражение (4) функцию правдоподобия для многомерного наблюдения:

Р(2к|хк) = (У(2к> ск Хк^1,

дУ (1 ) ( 2 к. с к. Х к) дУ (1 ) ( 2к. с к. Хк)

Y =

V

дг,

к (1)

дг

•к(М)

дУ(М)(2к> ск, Хк) дУ(М)(2к> ск, Хк)

дг,

к (1)

ск> Хк) =

дг

•к(М)

У(1)(2к> ск> Хк) У(М)(2к> ск> Хк)

где Yv — якобиан преобразования от вектора переменных w к вектору 2к; у(/)(гк, ск, хк), у = 1, ..., М — полученные в результате обратного преобразования соответствующих компонент х(-.) однозначно определенные функции; г^, у = 1, ..., М — компоненты вектора наблюдения 2к.

Так как АПВ р(хк _ 11 2к -1; с к - 1) в правой части

равенства (4) — известная функция (определенная на предыдущем шаге), рекуррентный алгоритм определения АПВ переменной состояния для к-го момента времени при наличии дискретных отсчетов сигнала наблюдения 21 принимает вид:

р(хк| 21; ск) =

_ Л(Хк, ск)

h* (ск)

(5)

л(Хк ск) = |... |р(Хк _ К-1; ск _ 1) х

-да -да

X 4(7^ Хк _ ^ЭД^к - 1* Ck, х*))^

Ч - 1>>\ Ау1 к -

+ да + да

^(ск) = { ... { л(х^ ск)^Хк.

-да -да

Тогда, с учетом выражения (5), обобщенный вероятностный критерий (3) можно окончательно представить следующим образом:

/ = | Ф

Л ( Хк, с к)

^(ск)

^к.

(6)

Идентификация параметра ск предполагает минимизацию (максимизацию) критерия (6). Для этой цели возможно применение известных и широко применяемых методов оптимизации функций многих переменных, выбор которых определяется особенностями исследуемого объекта и его наблюдателя [12—14]. Другими словами, в конечном счете, выбор зависит от конкретного вида критериального выражения — аналитических свойств получаемой целевой функции (выпуклости, непрерывности, повторной дифференцируемости и др.), а также необходимого для реализации самого метода оптимизации уровня информации об этой функции. Так, в случае вида критериального выражения, близкого к квадратичному, целесообразно применять методы, использующие свойства сопряженных направлений (метод Давидона— Флетчера—Пауэлла, обеспечивающий в подобном случае относительно быструю сходимость за конечное число шагов, и т. п.), тогда как для случаев недифференцируемых функций неизвестной формы — методы прямой оптимизации (метод Нелде-ра—Мида и его модификации, метод Хука—Джив-са и др.).

Эффективность предложенного подхода проиллюстрируем на следующем примере, имеющем самостоятельное теоретическое и практическое значение.

3. РЕШЕНИЕ ЗАДАЧИ НЕЛИНЕЙНОЙ ПАРАМЕТРИЧЕСКОЙ ИДЕНТИФИКАЦИИ НА ОСНОВЕ КРИТЕРИЯ МИНИМУМА ВЕРОЯТНОСТИ ОШИБКИ ОЦЕНИВАНИЯ

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

+да +да

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

Итак, в качестве возможного варианта критерия J примем далее критерий минимума АПВ текущей ошибки оценивания стк переменных состояния объекта на выбранном интервале ее предельно допустимого изменения — от а • до а „ . т. е.

^л.1п тах

тт J = тт 1 ... 1 р(стк| zí )<1а1

к

где стк = х

к

хк — вектор текущей ошибки оце-

нивания, хк — вектор текущих оценок переменных состояния объекта, р(стк^к) — многомерная

условная апостериорная плотность вероятности ошибки оценивания.

Учитывая линейную зависимость значений текущей ошибки ак и переменной состояния хк:

а к = хк — хк, выразим АПВ ошибки оценивания р(стк| z1) через определенную ранее АПВ переменной состояния p(xk| zk; Ск):

р(ак1 zí) = РК + хк | z1; Ск).

В этом случае минимизацию критерия идентификации можно представить следующим образом:

тт J = min 1 ... 1 р(стк| zí )<1ак

тс[п 1 ... | p(стk + хк|zk; С^ак

(7)

Производя соответствующую замену переменных в полученном ранее выражении для АПВ параметров состояния (5) и обозначив критериальное выражение через 0(Ск), поиск минимума критерия (7) можно представить в следующем компактном виде:

тт J = min 0(Ск),

(8)

где

ашах ашах

о(Ск) = 1 ... 1 р(^к + хк^к; Ск)^ак :

ашт ашт

шах шах

1 ... 1

ашт ашт

л ( ^к + хк. С к) А*(Ск)

Лст,,

Как было отмечено, идентификация текущего значения параметра, удовлетворяющего критерию (8), предполагает минимизацию функции ^(Ск) на основе известных методов оптимизации функций многих переменных, конкретное применение которых рассмотрим далее в процессе численного моделирования конкретной процедуры идентификации.

4. ПРИМЕР

Проиллюстрируем эффективность предложенного метода идентификации следующим примером.

Пусть стохастический дискретный объект задан нелинейным разностным уравнением

хк = /(хк _ 1) + П = хк - 1 + АХ2 - 1 + П,

Х1 = 0,5, (9)

где п — белый гауссовский шум с нулевым средним и дисперсией Бп = 0,025, а = —0,4 — коэффициент уравнения объекта.

Переменные состояния заданного объекта наблюдаются измерителем, описываемым уравнением

2

^ = х(с, Хк) + я = ЬсХк + я, (10)

где Ь = 0,35 — коэффициент уравнения объекта, с — неизвестный искомый параметр наблюдателя (для рассматриваемого далее модельного примера выберем исходное значение с = 1 для всех к), я — белый гауссовский шум с нулевым средним и дисперсией = 0,15.

Для последующего определения оценок воспользуемся достаточно эффективным в смысле вычислений рекуррентным нелинейным фильтром Калмана (предполагающим линеаризацию нелинейных функций /(хк _ 1) и х(с, хк) в окрестности оценки). Линеаризованные уравнения объекта (9) и измерителя (10) принимают следующий вид:

хк = - 1 + вк + п Х1 = 0,5 ^ = с(ЕкХк + рк) + ^ где Ак = 1 - 0,8 хк -1, Вк = -0,4( хк -1 )2, Ек = 0,7 хк -1, Рк = -0,35( хк -1) — коэффициенты, полученные в результате линеаризации функций /(х_1) и %(с, хк) в окрестностях оценок переменной состояния объекта для (к — 1)-го шага хк -1.

Текущие значения оценки переменной состояния объекта определялись посредством субоптимального фильтра Калмана [10], имеющего в рассматриваемом случае вид

хк = и(с) = Ак хк -1 + Вк + + Кк(% — с(Ек(Акхк -1 + Вк) + Рк)^ х* = 0,7,

К = 4, Я, =

1

2 ч -1

як - 1 + Бп

+ с-) , я1 = 0,5.

а

а

шах

шах

а

а

а

а

шах

шах

а

а

а

а

шах

шах

а

а

к

к

Соответственно, функция О(с) для к-го шага алгоритма примет вид

сттах

п(с) = | р(СТк + Хк |гк1; с)^ =

атт атах / +да

= / I ^ / ^(Хк - 11 г1-1) ^ Х

а . ^ h (с) -да 2П

(стк + и( с) -/(Хк - 1))2

х е

2Д„

аХ

к - 1

х е

( гк - х ( с, а к + - с ) ) ) 2 2 Д..

А--'

(11)

П 1, щ)

5 / /

4 / / \ \

3 / / 1 \

\ 1 \ \ \

/2 У'

Рис. 2. Сечение вдоль оси стк функции У(с, стк) при с = 1 и гауссовская плотность распределения с математическим ожиданием 0,29 и дисперсией 0,064 (к = 150)

к -1

Г(с) = | | р(Хк - 1 г1-1)

( а к + Д( с ) -/( _ , ) ) 2

2 D „

ял72П

^к - 1 х

(гк-х(с, ак + и(с)))

2Д„

аХк.

Априорную плотность вероятности для первой итерации алгоритма выберем нормальной с дисперсией = 0,45 и математическим ожиданием 0,7. Здесь интересно отметить, что отклонения среднего значения априорной плотности от начального значения переменной состояния не оказывают в дальнейшем существенного влияния на качество процедуры идентификации (алгоритм идентификации к ним устойчив).

Рис. 1. График АПВ текущей ошибки оценивания для к-го шага алгоритма (к = 150)

На рис. 1 представлен полученный в результате моделирования график входящей в формулу (11) АПВ текущей ошибки оценивания для к-го шага алгоритма (к = 150), которая, являясь функцией текущей ошибки стк, зависит в то же время от значений искомого параметра с: (р(стк + Х* | гк) =

= р(^к + и(с)|гк) = Г(с, Стк)).

Интересно, что характер получаемых в результате моделирования АПВ (как для переменной состояния нелинейного объекта, так и для ошибки оценивания <гк) существенно негауссовский (ввиду нелинейности исходных уравнений объекта и измерителя). Отмеченное обстоятельство иллюстрируется приведенными на рис. 2 графиками, соответствующими сечению вдоль оси стк функции Г(с, стк) при с = 1 (непрерывная линия) и гауссов-ской плотности распределения с математическим ожиданием 0,29 и дисперсией 0,064 (штриховая линия).

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

Интегралы в выражении (11) здесь и в дальнейшем определялись численно с помощью квадратурных формул с шагом А = 0,025. Бесконечные пределы интегрирования по переменной состояния Х были заменены на конечные, удовлетворяющие точностным требованиям к алгоритму оценки

(Хшт = -0,5; Хтах = 1,5).

Эффективность работы алгоритма иллюстрируется приведенным на рис. 3 графиком зависимости критериального выражения О(с) от искомого параметра.

Очевидно, что минимум критериального выражения (и, следовательно, сечение построенной

X

+ да + да

ли —со

х

е

х

е

Рис. 3. Зависимость функции W(c) от искомого параметра (к = 150)

функции Г(с, ак) вдоль оси ак с наименьшей площадью) располагается около истинного значения искомого параметра с = 1.

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

Для минимизации выбранного критерия на очередном шаге алгоритма — т. е. для однозначного определения численного значения текущей оценки искомого параметра — целесообразно, задав некоторый интервал возможных значений параметра с (в примере был выбран интервал 0 < с < 6), воспользоваться одним из методов прямой минимизации. В данном случае применялся модифицированный симплексный метод прямой минимизации Нелдера—Мида — метод Бокса, обладающий достаточной вычислительной эффективностью и удобной программной реализацией [14].

Результаты компьютерного моделирования процедуры нелинейной параметрической идентификации показали, что при выборе числа дискретных значений сигнала наблюдения к 1 280 отклонение оценки параметра с наблюдателя от его истинного значения не превышает 17 % от его значения.

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

ЗАКЛЮЧЕНИЕ

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

— более высокий по сравнению с существующими методами уровень потенциальной точности

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

— инвариантность к виду плотности распределения вероятности шума как объекта, так и наблюдателя;

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

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

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

ЛИТЕРАТУРА

1. Гроп Д. Методы идентификации систем. — М.: Мир, 1979. — 302 с.

2. Теория моделей в процессах управления / Б.Н. Петров, Г.М. Уланов, И.И. Гольденблат, С.В. Ульянов. — М.: Наука, 1978. — 223 с.

3. Справочник по теории автоматического управления / Под ред. А.А. Красовского. — М.: Наука, 1987. — 712 с.

4. Штейнберг Ш.Е. Идентификация в системах управления. — М.: Энергоатомиздат, 1987. — 80 с.

5. Эйкхофф П. Основы идентификации систем управления. — М.: Мир, 1975. — 683 с.

6. Ли Р. Оптимальные оценки, определение характеристик и управление. — М.: Наука, 1966. —176 с.

7. Пугачев В.С., Казаков И.Е., Евланов Л.Г. Основы статистической теории автоматических систем. — М.: Машиностроение, 1974. — 400 с.

8. Пащенко Ф.Ф. Введение в состоятельные методы моделирования систем. Идентификация нелинейных систем. — М.: Финансы и статистика, 2007. — 288 с.

9. Хуторцев В.В., Соколов С.В, Шевчук П.С. Современные принципы управления и фильтрации в стохастических системах. — М.: Радио и связь, 2001. — 807 с.

10. Тихонов В.И., Харисов В.Н. Статистический анализ и синтез радиотехнических устройств и систем. — М.: Радио и связь, 1991. — 608 с.

11. Левин Б.Р. Теоретические основы статистической радиотехники. — М.: Радио и связь, 1989. — 455 с.

12. Измаилов А.Ф., Солодов М.В. Численные методы оптимизации. — М.: Физматлит, 2003. — 304 с.

13. Васильев Ф.П. Методы оптимизации. — М.: Факториал, 2002. — 824 с.

14. Банди Б. Методы оптимизации: Вводный курс. — М.: Радио и связь, 1988. — 128 с.

e-mail: [email protected]

Статья представлена к публикации членом редколлегии Ф.Ф. Пащенко. □