В.В. Александров, С.С. Лемак, Н.А. Парусников - Лекции по механике управляемых систем, страница 13

Число измерений на каждом участке кроме первогообозначим через m (поскольку ΔT = const ). Тогда первый участок[T1 , T2 ] содержит m + 1 измерений. Имеет место соотношениеm + 1 + (n − 2)m = N + 184илиm(n − 1) = N.Введем безразмерный интервал между измерениями:ΔTΔt,(= m).Δτ =ΔTΔtНа каждом участке [Ti−1 , Ti ] введем безразмерное время τ :t − Ti−1, τ ∈ [0, 1].τ=ΔTОсновной элемент сглаживания при помощи кубических сплайнов— определение базовых функций. Вводятся четыре базовые функции:13ϕ1 (τ ) = (1 − τ )3 ,ϕ2 (τ ) = 1 − (2 − τ )τ 2 ,4431 32ϕ3 (τ ) = 1 − (1 + τ )(1 − τ ) ,ϕ4 (τ ) = τ .44Линейная независимость функций ϕ1 (τ ), ϕ2 (τ ), ϕ3 (τ ), ϕ4 (τ ) следует,например, из однозначности определения коэффициентов c1 , c2 , c3 , c4в соотношенииc1 ϕ1 (τ ) + c2 ϕ2 (τ ) + c3 ϕ3 (τ ) + c4 ϕ4 (τ ) =a1 + a2 τ + a3 τ 2 + a4 τ 3 .Оценка g̃ функции g ищется в виде!ci ϕ1 (τ ) + ci+1 ϕ2 (τ ) + ci+2 ϕ3 (τ ) + ci+3 ϕ4 (τ ), t ∈ [Ti , Ti+1 ],g̃ =0, t ∈/ [Ti , Ti+1 ].Непосредственной проверкой легко убедиться в том, что g̃(t)непрерывна вместе со своей первой и второй производной в узловыхточках и, стало быть, всюду.Для определения коэффициентовc1 , c2 , .

. . , cn , cn+1 , cn+2 ,имеем следующую переопределенную систему, имеющую так называемую ленточную структуру.Интервал [T1 , T2 ]:z0 = c1 ϕ1 (0) + c2 ϕ2 (0) + c3 ϕ3 (0) + c4 ϕ4 (0),z1 = c1 ϕ1 (Δτ ) + c2 ϕ2 (Δτ ) + c3 ϕ3 (Δτ ) + c4 ϕ4 (Δτ ),z2 = c1 ϕ1 (2Δτ ) + c2 ϕ2 (2Δτ ) + c3 ϕ3 (2Δτ ) + c4 ϕ4 (2Δτ ),...zm = c1 ϕ1 (1) + c2 ϕ2 (1) + c3 ϕ3 (1) + c4 ϕ4 (1).85Интервал [T2 , T3 ]:zm+1 = c2 ϕ1 (Δτ ) + c3 ϕ2 (Δτ ) + c3 ϕ3 (Δτ ) + c5 ϕ4 (Δτ ),...z2m = c2 ϕ1 (1) + c3 ϕ2 (1) + c4 ϕ3 (1) + c5 ϕ4 (1)....Интервал [Tn−1 , Tn ]:...zN = cn−1 ϕ1 (1) + cn ϕ2 (1) + cn+1 ϕ3 (1) + cn+2 ϕ4 (1).Решение этой системы обычно находят по методу наименьшихквадратов.Заметим, что алгоритм сглаживания может успешно применяться,если некоторые из измерений или даже целые блоки таких измеренийотсутствуют.86Лекция 12Критерий ортогональности и критерийусловного среднегоДалее решается задача оценивания, которая служит основой многочисленных модификаций метода наименьших квадратов.

Из нее, вчастности, следует и предыдущий результат.Рассматриваются два случайных вектора x(n × 1) и z(m × 1) сизвестными характеристиками:μz = M [z],μx = M [x],◦◦Pxx = M [xx ],◦◦Pzz = M [z z ],◦◦Pxz = M [xz ].Указанные характеристики будем называть априорной информацией.Ставится задача определения (построения оценки) скалярной величины α = c x, где c(n × 1) — известный вектор.В рамках априорной информации в качестве оценки α̃ величины αестественно принять ее математическое ожидание, т.е. положить α̃ =c μx .

Ошибка оценки Δα = α − α̃. Мерой ошибки оценки служит◦◦◦дисперсия D[Δα] = M [(c x)2 ] = M [c xx c] = c Pxx c.Пусть теперь измеряется вектор z и результат измерения известен.Поскольку известна мера линейной связи двух векторов x и z (задана матрица Pxz ), то возникает задача получения оценки величины α сучетом измерения z.Оценку будем искать как линейную, несмещенную и оптимальнуюпри условии, что критерием оптимальности служит минимум дисперсии ошибки оценки Δα.1. Линейность. Оценка α̃ ищется в видеα̃ = γ 1 μx + γ 2 μz + γ 3 z,где векторы γ 1 ,γ 2 и γ 3 подлежат определению.2. Несмещенность.

Условие несмещенности означает, что в среднем оценивание приводит к точному результату, т.е. M [Δα] = 0. ИмеемM [Δα] = M [c x−(γ 1 μx +γ 2 μz +γ 3 z)] = (c−γ 1 ) μx −(γ 2 +γ 3 ) μz ,отсюда γ 1 = c, γ 2 = −γ 3 . Переобозначим γ 3 = κ, тогдаα̃ = c μx + κ (z − μz ),87◦◦Δα = c x − κ z.Вектор κ подлежит дальнейшему определению.3. Оптимальность.D[Δα] = M [(Δα)2 ] = c Pxx c + κ Pzz κ − 2c Pxz κ.Условие оптимальности имеет вид∂D[Δα]= 0,∂κоткуда следует−1κ = c Pxz Pzz= c K,(12.1)(12.2)−1Pxz Pzz.где K =Если в качестве оценки x̃ вектора x выбратьx̃ = μx + K(z − μz ),(12.3)то при всяком c оценка α̃ = c x̃ оптимальна.Мерой ошибки оценки Δx = x − x̃ вектора x служит ковариация◦◦ ◦◦PΔx = M [ΔxΔx ] = M [(x − K z)(x − K z) ].

После выкладок сучетом (12.2) получаем−1Pzx .(12.4)PΔx = Pxx − Pxz PzzМерой ошибки оценки величины α служит дисперсияD[Δα] = c PΔx c.Величины x̃ и P̃xx можно интерпретировать как апостериорные характеристики (условные, уточненные в результате учета измерения z).Результаты, сведенные в формулы (12.2), (12.3), позволяют критерийоптимальности (12.1) сформулировать в виде критерия ортогональности: оценка x̃ оптимальна, если ошибка оценки ортогональна (не коррелирована) измерениям:◦M [Δxz ] = 0 .(12.5)◦◦В самом деле, из (12.5) следует (12.1) с учетом того, что Δx = x − K z.Как уже говорилось, критерий (12.5) служит отражением простогоэвристического рассуждения: оптимальная оценка должна вобрать всебя всю информацию об оцениваемой величине, содержащуюся в измерениях, и, следовательно, ошибка оценки не должна быть связана сизмерениями.Ввиду важности полученного результата дадим ему несколькоиную интерпретацию, напрямую используя нормальность закона распределения случайных векторов x и z.

Оценка x̃ будет определятьсякак условное математическое ожидание.88Пусть случайный вектор xu=zраспределен по нормальному закону с характеристиками ◦ ◦Pxx PxzPu = M uu =),, (Pxz = PzxPzx Pzzμxμu = M [u] =.μz(12.6)Функция плотности вероятности f (u) = f (x, z) вектора u записывается в следующем виде:!"11 −1exp − (u − μu ) Pu (u − μu ) . (12.7)f (u) =m+n12(2π) 2 |Pu | 2Пусть вектор z доступен измерению. Требуется уточнить априорную характеристику M [x] случайного вектора x и определить соответствующую ковариацию.Для решения этой задачи естественно воспользоваться формулойусловной вероятности:f (x, z),f (x|z) =f (z)где функция распределения f (x, z) определяется соотношением(12.7), а функция f (z) имеет вид!"11 −1(z−μf (z) =exp−)P(z−μ).zzm1zz2(2π) 2 |Pzz | 2Обозначим через μx|z и Px|z условное математическое ожидание иусловную ковариационную матрицу.

Тогда!"11 −1− (x − μx|z ) Px|z (x − μx|z ) ,f (x|z) =n1 exp2(2π) 2 |Px|z | 2 Px|z = M x − μx|z x − μx|z z .Далее нам понадобится соотношение−1Pxx PxzNxx−1Pu ==N =Pzx PzzNzxгде−1−1Pzx,Nxx = Pxx − Pxz PzzNxzNzz,−1−1Nzz = Pzz − Pzx PxxPxz,89−1−1−1Nxz = − Pxx − Pxz PzzPzxPxz Pzz.Справедливость последних соотношений непосредственно следует изравенства Pu N = E.Применяя формулу Байеса и приравнивая коэффициенты при одинаковых степенях x, получаем уже известный результат−1μx|z = μx + Pxz Pzz(z − μz ),−1Pzx .Px|z = Pxx − Pxz Pzz(12.8)Рассмотрим частный случай, когда вектор z линейно зависит от вектора x, т.е.z = Hx + r,H — известная матрица, r — случайный вектор с нулевым средним,трактуемый как погрешность измерения.Будем считать заданными априорные характеристикиμx = M [x],◦◦Pxx = M [xx ],R = M [rr ],и, кроме того, предположим, что вектор r не коррелирован с оцениваемым вектором x:◦M [xr ] = 0.Целесообразно ввести обозначения, в которых априорность фиксируется верхним индексом −, апостериорность — верхним индексом +:x̃− = μx , Px− = Pxx .После очевидных выкладок имеемPz− = Pzz = HPx− H + R,Pxz = Px− H ,μz = Hμx = H x̃− .Взамен соотношений (12.8) получимx̃+ = x̃− + K(z − H x̃− ),K = Px− H (HPx− H + R)−1 ,Px+90=Px−−KHPx− .(12.9)Лекция 13Дискретный фильтр КалманаДальнейшее изложение посвящено динамическим задачам оценивания, когда поведение во времени оцениваемой величины x описывается линейными стохастическими (непрерывными или дискретными) уравнениями.

При этом будет использоваться только одна формакритерия оптимальности — критерий ортогональности. Обсуждениеуказанного круга задач начнем с дискретного случая.Термин фильтрация означает, что алгоритм может быть интерпретирован как алгоритм выделения полезной информации на фонепомех.1. Алгоритмы дискретного фильтра КалманаСначала опишем дискретный вариант задачи, используя при этомрезультаты лекций 10 и 12.Поведение динамической системы подчиняется уравнениюxj+1 = Φj xj + qj ,где xj — вектор состояния в момент tj ; qj — дискретный белый шумс известной интенсивностью — ковариационной матрицей Qj ; Φj —известная переходная матрица.Заданы x̃0 = M [x0 ], P0 = M [(x0 − x̃0 )(x0 − x̃0 ) ].

В моментывремени j = 0, 1, 2, . . . поступает информация z0 , z1 , . . . :zj = Hj xj + rj .Погрешность информации rj — некоррелированный во времени вектор с нулевым средним и известной интенсивностью Rj : M [rj ri ] =Rj δji . Кроме того, полагаем M [xi rj ] = 0.Требуется в каждый момент j получить оценку x̃j , удовлетворяющую условию линейности, несмещенности, оптимальности. Введемобозначения:x̃−j — оценка вектора x в момент времени j, использующая измерения z0 , .