Автоматика и телемеханика, № 4, 2023

Стохастические системы

Автор посвящает эту статью

светлой памяти своего учителя и друга

Бориса Теодоровича Поляка

(Институт проблем управления им. В.А. Трапезникова РАН, Москва;

Национальный исследовательский университет

“Московский физико-технический институт”)

СРАВНЕНИЕ ГАРАНТИРУЮЩЕГО

И КАЛМАНОВСКОГО ФИЛЬТРОВ¹

Предлагается новый подход к задаче фильтрации при произвольных

ограниченных внешних возмущениях, основанный на ее сведении к за-

даче оптимизации. Подход обладает невысокой вычислительной сложно-

стью, предполагая на каждом итерационном шаге лишь решение уравне-

ний Ляпунова. При этом его преимуществами, существенными с инже-

нерно-практической точки зрения, являются возможность ограничения

величины матрицы фильтра, а также возможность строить оптимальные

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

отдельности.

Выписан градиентный метод для отыскания матрицы фильтра. Как

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

весьма эффективной и приводящей к вполне удовлетворительным резуль-

татам. Статья продолжает серию работ, посвященную синтезу обратной

связи в задачах управления с позиций оптимизации.

Ключевые слова: линейная система, внешние возмущения, фильтрация,

фильтр Калмана, наблюдатель Люенбергера, оптимизация, уравнение

Ляпунова, градиентный метод, метод Ньютона, сходимость.

DOI: 10.31857/S0005231023040050, EDN: QIBUNX

1. Введение

Классическая постановка задачи фильтрации (т.е. оценки состояния ди-

намической системы по измерениям) при случайных возмущениях допускает

практически исчерпывающее решение с помощью фильтра Калмана [1]; см.

подробнее монографию [2], а также [3]. Однако часто известно лишь, что все

возмущения являются ограниченными, а в остальном произвольными; в этом

¹ Исследование выполнено при частичной финансовой поддержке Российского научного

фонда (проект № 21-71-30005), https://rscf.ru/project/21-71-30005/

случае можно строить гарантирующие (а не вероятностные) оценки состоя-

ний. В начале 1970-х гг. в работах Швеппе [4], Куржанского [5] и Черно-

усько [6] была развита эллипсоидальная техника фильтрации. Позже в [7, 8]

рассматривалась проблема фильтрации с ограниченными неслучайными воз-

мущениями для стационарных задач: искалась оценка состояния такая, что

ее ошибка гарантированно заключена в единый эллипсоид для всех моментов

времени, т.е. оценка является равномерной, а сам фильтр искался в классе

линейных стационарных фильтров. В этом классе задач и оценок проблема

оказалась полностью разрешимой: удалось построить оптимальный фильтр

и оценку состояния. С технической точки зрения в [7, 8] был применен аппа-

рат линейных матричных неравенств (ЛМН) [9], а исходная задача сведена к

параметрической задаче полуопределенного программирования. Системати-

ческое изложение этой техники можно найти в монографии [10].

С другой стороны, в последнее время стал очень популярным подход к за-

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

обоснование подобных методов появилось лишь недавно, см. [11-15]. В [16] по-

добный подход был впервые применен к задачам с внешними возмущениями,

в [17] — к задаче синтеза обратной связи по выходу при помощи наблюдателя,

а в [18] — к задаче синтеза ПИД-регуляторов.

Настоящая статья, цель которой состоит в упрощении алгоритма вычис-

ления гарантирующего фильтра и его численном сравнении с фильтром Кал-

мана, продолжает обе эти линии исследований. В ней предлагается оптимиза-

ционный алгоритм решения задачи фильтрации при неслучайных ограничен-

ных внешних возмущениях. Он обладает невысокой вычислительной сложно-

стью, предполагая на каждом итерационном шаге лишь решение уравнений

Ляпунова. При этом существенным с инженерно-практической точки зрения

преимуществом предлагаемого подхода по сравнению с ЛМН-подходом явля-

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

ют примеры, предлагаемая рекуррентная процедура является эффективной

и приводящей к вполне удовлетворительным результатам.

При этом важно отметить, что в отличие от фильтра Калмана в рамках

предлагаемого подхода оказывается возможным строить оптимальные мат-

рицы фильтра по каждой из координат вектора состояния системы в отдель-

ности.

Структура статьи следующая: раздел 2 содержит постановку задачи, в

разделе 3 обсуждается предлагаемый подход к построению гарантирующего

фильтра, в разделе 4 приводится и обосновывается алгоритм вычисления

оптимальной матрицы фильтра, раздел 5 посвящен непрерывной постановке

задачи, раздел 6 содержит описание и обсуждение результатов вычислений

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

полученных результатов.

Всюду далее | · | — евклидова норма вектора, ∥ · ∥ — спектральная нор-

ма матрицы, ∥ · ∥_F — фробениусова норма матрицы,^T — символ транспо-

нирования, tr — след матрицы, 〈·, ·〉 — скалярное произведение Фробениуса

для матриц, I — единичная матрица соответствующей размерности, λ_i(A) —

собственные значения матрицы A, σ_i(A) — сингулярные числа матрицы A,

а σ(A) = - maxRe (λi(A)) — степень устойчивости гурвицевой матрицы A,

ρ(A) = max|λi(A)| — спектральный радиус шуровской матрицы A. Все мат-

ричные неравенства понимаются в смысле знакоопределенности матриц.

2. Постановка задачи

Рассмотрим систему в дискретном времени

xk+1 = Ax_k + B₁u_k + D₁w_k,

(1)

y_k = Cx_k + B₂u_k + D₂w_k,

z_k = C₁x_k,

где A ∈ R^n×n, B₁ ∈ R^n×p, B₂ ∈ R^ℓ×p, C ∈ R^ℓ×n, C₁ ∈ R^r×n, D₁ ∈ R^n×m,

D₂ ∈ R^ℓ×m, с состоянием x_k ∈ Rⁿ, начальным условием x₀, входом u_k ∈ R^p,

наблюдаемым выходом y_k ∈ R^ℓ, оцениваемым выходом z_k ∈ R^r и внешним

возмущением (шумом) w_k ∈ R^m, ограниченным в каждый момент времени:

|w_k| ≤ 1

для всех k = 0, 1, . . . ;

пара (A, D₁) управляема, пара (A, C) наблюдаема.

Пусть состояние x_k системы недоступно измерению и информация о систе-

ме предоставляется ее выходом y_k. Для оценивания выхода z_k будем использо-

вать фильтр, описываемый линейным разностным уравнением относительно

оценки состояния x_k:

(2)

xk+1 = Axk + B1uk + L(yk - Cxk - B2uk),

x₀

= 0,

где L ∈ R^n×ℓ. Подчеркнем, что структура фильтра задается заранее — он

является линейным стационарным, подлежит выбору лишь постоянная мат-

рица L. Эта структура фильтра такая же, как и в известном наблюдателе

Люенбергера [19, 20]. По сути, можно рассматривать этот фильтр как обоб-

щение наблюдателя Люенбергера на задачи с помехами.

Задачей является минимизация ошибки оценки

z_k - z_k = C₁(x_k - x_k) = C₁e_k,

где e_k = x_k - x_k — невязка, удовлетворяющая согласно (1), (2) разностному

уравнению

(3)

ek+1 = (A - LC)e_k + (D₁ - LD₂)w_k, e₀ = x₀.

Обратим внимание, что допустимая матрица фильтра L стабилизирует си-

стему (3), обращая матрицу A-LC в шуровскую. Ее существование вытекает

из свойства наблюдаемости исходной системы.

Важно отметить, что здесь рассматривается случай неслучайных ограни-

ченных помех. Для случайного, гауссовского шума вполне естественно приме-

нять калмановскую фильтрацию, однако одна из целей настоящей статьи —

привлечь внимание к фильтрации с ограниченными шумами. Также на де-

монстрационных примерах будет рассмотрено, как работает фильтр Калмана

с ограниченными помехами и наоборот — к чему приведет применение к слу-

чайным помехам рассматриваемой модели с ограниченным шумом.

Заметим также, что в настоящей статье рассматривается более общая по-

становка задачи, чем в [7]: в описании системы присутствует вход u_k.

Наконец, учет ограничения на внешнее возмущение вида

|w_k| ≤ γ для всех k = 0, 1, . . .

производится очевидным путем масштабирования матриц: D₁ := γD₁, D₂ :=

:= γD₂.

3. Гарантирующий фильтр

В настоящей статье предлагается гарантирующий подход к решению за-

дачи фильтрации при ограниченных шумах, для которого можно явным об-

разом выписать прямые формулы, основанные на градиентом спуске. Он ос-

новывается на концепции инвариантных эллипсоидов, см. подробнее [10, 21].

Один из основных результатов статьи представлен следующим утвержде-

нием.

Теорема 1. Пусть L^∗, P^∗ — решение оптимизационной задачи

(4)

min f(L, α), f(L, α) = tr C₁P C^T1 + ρ∥L∥2F ,

при ограничении

(5)

(A - LC)P (A - LC)^T - P +

(D₁ - LD₂)(D₁ - LD₂)^T

1-α

относительно матричных переменных P = P^T ∈ R^n×n, L ∈ R^n×ℓ и скаляр-

ного параметра 0 < α < 1.

Тогда ошибка оценки z_k - z_k выхода системы (1) с нулевым начальным

условием при помощи наблюдателя (2) с матрицей L^∗ заключена в мини-

мальный ограничивающий эллипсоид с матрицей

C₁P^∗C^T1.

Переходя к обоснованию теоремы, напомним следующий результат [22].

Лемма 1. Пусть матрица A шуровская, ρ = max|λi(A)| < 1, пара (A, D)

управляема, а матрица P(α) ≻ 0, ρ² < α < 1, удовлетворяет дискретному

уравнению Ляпунова

(6)

AP A^T - P +

DD^T

= 0.

1-α

Тогда:

1) задача об оптимальном ограничивающем эллипсоиде для выхода систе-

мы

xk+1 = Ax_k + Dw_k,

z_k = Cx_k,

с начальным условием x₀ и ограниченными внешними возмущениями |w_k| ≤

≤ 1 сводится к минимизации одномерной функции f(α) = trCP(α)C^T на

интервале ρ² < α < 1;

2) если α^∗ — точка минимума и x₀ удовлетворяет условию

x^T0P-1(α^∗)x₀ ≤ 1,

то гарантируется оценка

|z_k|² ≤ f(α^∗), k = 0, 1, . . .

Итак, рассмотрим величину C₁e_k в качестве линейного выхода системы (3).

Тогда если заключить невязку e_k в инвариантный эллипсоид

{

}

E =

e∈Rⁿ: e^TP-1e≤1

P ≻ 0,

то C₁e_k будет содержаться в ограничивающем эллипсоиде

{

}

(7)

E_z =

e_z ∈ R^r : eTz(C₁PC^T1)-1e_z ≤ 1

размер которого и будем минимизировать. Таким образом, оценивается

асимптотическая (а при малых уклонениях и равномерная по k) точность

фильтрации.

В соответствии с леммой 1 исходная задача свелась к матричной оптимиза-

ционной задаче (4)-(5). Обратим внимание, что в минимизируемую функцию

f (L, α) помимо компоненты, определяющей размер ограничивающего эллип-

соида (7) по критерию следа, введен штраф за величину матрицы фильтра

(при этом коэффициент ρ > 0 регулирует его важность). В то же время его

наличие гарантирует коэрцитивность минимизируемой функции по L. Запись

f (L, α) подчеркивает, что при заданных L и α матрица P находится из урав-

нения Ляпунова (5); тем самым независимыми переменными являются L и α.

Сделаем важное

Замечание 1. Рассматриваемый гарантирующий подход позволяет

строить оптимальные матрицы фильтра по каждой из координат вектора

состояния системы в отдельности (этой возможности нет в фильтре Калма-

на). Действительно, положив в задаче из теоремы 1 в качестве матрицы C₁

транспонированный i-й координатный вектор, приходим к матрице фильтра,

минимизирующей невязку x(i)k - x(i)k.

4. Вычисление оптимальной матрицы фильтра

Напомним, что в [7] был предложен вариант гарантирующего подхода к

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

технике линейных матричных неравенств и предполагает решение парамет-

рической задачи полуопределенного программирования. В рамках рассмат-

риваемой оптимизационной задачи (4)-(5) нет необходимости применять этот,

технически сравнительно сложный, аппарат (несмотря на то, что и мини-

мизирумая функция, и ограничение невыпуклы по совокупности перемен-

ных P , L и α). В этом разделе будет предложен регулярный итеративный

подход к ее решению, в основе которого лежит применение градиентного ме-

тода по переменной L и минимизации по α по методу Ньютона. Приведем

принципиальную схему алгоритма.

Алгоритм 1 для минимизации f(L,α):

1. Задаемся параметрами ε > 0, γ > 0, 0 < τ < 1 и начальным допустимым

(

)

приближением L₀. Вычисляем величину α₀ =

1 + ρ²(A - L₀C)

/2.

2. На j-й итерации, имея величины L_j и α_j, находим градиент H_j =

= ∇_Lf(L_j,α_j). Если ∥H_j∥ ≤ ε, то L_j принимаем за приближенное решение.

3. Делаем шаг градиентного метода

Lj+1 = L_j - γ_jH_j,

при этом длину шага γ_j > 0 подбираем дроблением γ до выполнения условий:

а. Lj+1 обращает матрицу (A - LC)/√αj в шуровскую;

б. f(Lj+1) ≤ f(L_j) - τγ_j ∥H_j ∥².

4. Для полученного Lj+1 решаем задачу минимизации f(Lj+1, α) по α по

методу Ньютона:²

f^′(α_j)

αj+1 = α_j -

f^′′(α_j)

и получаем αj+1. Переходим к п. 2.

В Алгоритме 1 величины ∇_Lf(L,α), f^′(α) и f^′′(α) определяются следую-

щим образом:

(

)

∇_Lf(L,α) = 2 ρL -

Y (A - LC)P C^T -

Y (D₁ - LD₂)D^T

1-α

(

)

f^′(α) = tr Y

(D₁ - LD₂)(D₁ - LD₂)^T -

(A - LC)P (A - LC)T ,

(1 - α)²

α²

(

f^′′(α) = 2tr Y

(D₁ - LD₂)(D₁ - LD₂)^T +

(1 - α)³

)

(A - LC)(P - X)(A - LC)^T

α³

² Реально требуется не более трех-четырех итераций для получения решения с

большой точностью, если начальная точка не слишком близка к границам интервала

(

)

ρ²(A - Lj+1C), 1

где матрицы P , Y и X являются решениями дискретных уравнений Ляпуно-

ва (5),

(8)

(A - LC)^TY (A - LC) - Y + C^T1C₁

(A - LC)X(A - LC)^T - X +

(D₁ - LD₂)(D₁ - LD₂)^T-

(1 - α)²

(A - LC)P (A - LC)^T = 0

α²

соответственно.

Обоснование Алгоритма 1 с выводом соответствующих формул помещено

в Приложение 1.

Важным моментом является выбор пробного шага градиентного метода.

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

дем для некоторого допустимого L решение P дискретного уравнения Ляпу-

нова

(A - LC)P (A - LC)^T - P = -I.

Рассмотрим приращение по L:

L → L - γH, H = ∇_Lf(L,α),

и найдем, для каких γ матрица P останется матрицей квадратичной функции

Ляпунова для A - (L - γH)C, т.е.

(A - (L - γH)C) P (A - (L - γH)C)^T - P ≺ 0

или по лемме Шура

(

)

A - (L - γH)C

≻ 0.

(A - (L - γH)C)^T

P-1

Полученное матричное неравенство, которому можно придать вид

(

)

(

)

A - LC

+γ

≻ 0,

(A - LC)^T P-1

(HC)^T

выполняется при

((

)

(

))

A - LC

γ<λ-1

^max (HC)^T

(A - LC)^T P-1

5. Непрерывный случай

Рассмотрим непрерывную систему

x = Ax + B₁u + D₁w, x(0) = x₀,

(9)

y=Cx+B₂u+D₂w,

z = C₁x,

где A ∈ R^n×n, B₁ ∈ R^n×p, B₂ ∈ R^ℓ×p, C ∈ R^ℓ×n, C₁ ∈ R^r×n, D₁ ∈ R^n×m,

D₂ ∈ R^ℓ×m, с состоянием x(t) ∈ Rⁿ, входом u(t) ∈ R^p, наблюдаемым выходом

y(t) ∈ R^ℓ, оцениваемым выходом z(t) ∈ R^r и внешним возмущением (шумом)

w(t) ∈ R^m, ограниченным в каждый момент времени:

|w(t)| ≤ 1 для всех t ≥ 0;

пара (A, D₁) управляема, пара (A, C) наблюдаема.

Для оценивания выхода z используется фильтр, описываемый линейным

диференциальным уравнением относительно оценки состояния x:

(10)

x = Ax + B₁u + L(y - Cx- B₂

u),

x(0) = 0,

где L ∈ R^n×ℓ.

Как и в дискретном случае, задачей является минимизация ошибки оценки

z - z = C₁(x - x) = C₁e,

где e(t) = x(t) - x(t) — невязка, удовлетворяющая согласно (9), (10) диффе-

ренциальному уравнению

(11)

ė = (A - LC)e + (D₁ - LD₂)w, e(0) = x₀.

При этом допустимая матрица фильтра L стабилизирует систему (11), обра-

щая матрицу A - LC в гурвицевую. Ее существование вытекает из свойства

наблюдаемости исходной системы.

Следующая лемма является непрерывным аналогом леммы 1.

Лемма 2 [9, 10]. Пусть матрица A гурвицева, σ = -maxRe (λi(A)) > 0,

пара (A, D) управляема, а матрица P (α) ≻ 0, 0 < α < 2σ, удовлетворяет

уравнению Ляпунова

(

)

(

)_T

A+

I P +P A+

DD^T = 0.

Тогда:

1) задача об оптимальном ограничивающем эллипсоиде для выхода систе-

мы

x = Ax + Dw, x(0) = x₀,

z=Cx

с ограниченными внешними возмущениями |w(t)| ≤ 1 сводится к минимиза-

ции одномерной функции f(α) = tr CP (α)C^T на интервале 0 < α < 2σ;

2) если α^∗ — точка минимума и x(0) удовлетворяет условию

x^T(0)P-1(α^∗)x(0) ≤ 1,

то гарантируется оценка

|z(t)|² ≤ f(α^∗),

0 ≤ t < ∞.

Рассуждая аналогично дискретному случаю и воспользовавшись лем-

мой 2, приходим к следующему результату.

Теорема 2. Пусть L^∗, P^∗ — решение оптимизационной задачи

min f(L, α), f(L, α) = tr C₁P C^T1 + ρ∥L∥2F ,

при ограничении

(

)

(

)_T

(12)

A-LC +

I P +P A-LC +

(D₁ - LD₂)(D₁ - LD₂)^T

относительно матричных переменных P = P^T ∈ R^n×n, L ∈ R^n×ℓ и скаляр-

ного параметра α > 0.

Тогда ошибка оценки z - z выхода системы (9) с нулевым начальным

условием при помощи наблюдателя (10) с матрицей L^∗ заключена в ми-

нимальный ограничивающий эллипсоид с матрицей

C₁P^∗C^T1.

Свойства минимизируемой функции и ее производных, установленные в

Приложении 2, позволяют построить метод минимизации и обосновать его

сходимость.

Алгоритм 2 для минимизации f(L,α):

1. Задаемся параметрами ε > 0, γ > 0, 0 < τ < 1 и начальным допустимым

приближением L₀. Вычисляем величину α₀ = σ(A - L₀C).

2. На j-й итерации, имея величины L_j и α_j, находим градиент H_j =

= ∇_Lf(L_j,α_j). Если ∥H_j∥ ≤ ε, то L_j принимаем за приближенное решение.

3. Делаем шаг градиентного метода

Lj+1 = L_j - γ_jH_j,

при этом длину шага γ_j > 0 подбираем дроблением γ до выполнения условий:

а. Lj+1 обращает матрицу A - LC +αj2Iвгурвицевую;

б. f(Lj+1) ≤ f(L_j) - τγ_j ∥H_j ∥².

4. Для полученного Lj+1 решаем задачу минимизации f(Lj+1, α) по α и по-

лучаем αj+1. Переходим к п. 2.

В Алгоритме 2 величины ∇_Lf(L,α), f^′(α) и f^′′(α) определяются следую-

щим образом:

(

)

∇_Lf(L,α) = 2 ρL - Y PC^T -

Y (D₁ - LD₂)D^T

(

)

f^′(α) = tr Y P -

(D₁ - LD₂)(D₁ - LD₂)T

α²

(

)

f^′′(α) = 2tr Y X +

(D₁ - LD₂)(D₁ - LD₂)T

α³

где матрицы P , Y и X являются решениями дискретных уравнений Ляпуно-

ва (12),

(

)_T

(

)

(13)

A - LC +

Y + Y A - LC +

+C^T1C₁

(

)

(

)_T

(14)

A - LC +

I X + X A - LC +

+P -

(D₁ - LD₂)(D₁ - LD₂)^T = 0.

α²

Обоснование Алгоритма 2 с выводом соответствующих формул помещено

в Приложение 2.

Выбор пробного шага градиентного метода предлагается проводить из сле-

дующих соображений. Найдем для некоторого допустимого L решение P

уравнения Ляпунова

(A - LC)P + P (A - LC)^T = -I.

Рассмотрим приращение по L:

L → L - γH, H = ∇_Lf(L,α),

и найдем, для каких γ матрица P останется матрицей квадратичной функции

Ляпунова для A - (L - γH)C, т.е.

(A - (L - γH)C) P + P (A - (L - γH)C)^T ≺ 0.

С учетом исходного уравнения имеем

(

)

HCP + P(HC)^T

≺I,

откуда

(

)

γ<λ-1max

HCP + P(HC)^T

6. Примеры и обсуждение

Рассмотрим дискретную модель объекта в пространстве состояний:

xk+1 = Ax_k + B₁u_k + Gw_k,

(15)

y_k = Cx_k + B₂u_k + v_k,

с состоянием x_k ∈ Rⁿ, начальным условием x₀, входом u_k ∈ R^p, наблюдаемым

выходом y_k ∈ R^ℓ, шумом w_k ∈ R^m и ошибкой измерений v_k ∈ R^ℓ; здесь A, B₁,

B₂, C и G — известные матрицы соответствующих размерностей; величи-

ны w_k и v_k предполагаются независимыми.

Применительно к системе (15) фильтр Калмана имеет следующий вид

(в предположении, что величины w_k и v_k распределены нормально с нулевым

математическим ожиданием и ковариационными матрицами Q и R соответ-

ственно).

Этап экстраполяции:

xk+1|k = Axk|k + B1uk,

Pk+1|k = AP_k|kA^T + GQG^T.

Этап коррекции:

K_k = Pk|k-1C^T(CPk|k-1C^T + R)-1,

xk|k = xk|k-1 + Kk(yk - Cxk|k-1 - B2uk),

P_k|k = (I - K_kC)Pk|k-1.

Будем рассматривать три следующие постановки задачи.

1. Модель M₁ со случайными помехами: шум w_k распределен по нормаль-

ному закону с нулевым математическим ожиданием и ковариационной мат-

рицей Q, а погрешность измерений v_k имеет нормальное распределение с

нулевым математическим ожиданием и ковариационной матрицей R:

w_k ∼ N(0,Q), v_k ∼ N(0,R).

2. Модель M₂ со случайными ограниченными помехами: величина w_k рав-

номерно распределена на кубе [-w, w]^m, а величина v_k равномерно распреде-

лена на кубе [-v, v]^ℓ:

w_k ∼ U([-w,w]^m), v_k ∼ U([-v,v]^ℓ).

3. Модель M₃ с неслучайными ограниченными помехами: величины w_k

и v_k могут принимать произвольные значения на кубах [-w,w]^m и [-v,v]^ℓ

соответственно:

|w_k|_∞ ≤ w,

|v_k|_∞ ≤ v.

В рамках этих моделей сравним работоспособность фильтра Калмана и

гарантирующего подхода на следующих примерах (несмотря на то, что га-

рантирующие оценки не обоснованы для модели M₁, а фильтр Калмана —

для модели M₃).

На приводимых ниже рисунках показаны истинная траектория системы,

ее наблюдение (если оно имеется) и оценки, предоставляемые калмановским

(ФК) и гарантирующим (ГФ) фильтром; в последнем случае также показана

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

Пример 1. Рассмотрим вагонетку, которая может двигаться без трения

по бесконечным рельсам [23]. В начальный момент вагонетка находится в

нулевом положении и на нее воздействуют внешние возмущения. Положе-

ние вагонетки измеряется каждые Δt секунд, при этом измерения неточны.

Задача состоит в отслеживании положения s и скорости s = v вагонетки.

Соответствующая система представима в формате (15) при

)

⁽s

x_k =

x₀ = 0,

)

⁽1 Δt

(

)

⁽(Δt)²/2⁾

B₁ = B₂ = 0, C =

Δt

1. В рамках модели M₁ предполагается, что на k-м такте вагонетка дви-

жется с постоянным ускорением, распределенным по нормальному закону

с нулевым математическим ожиданием и среднеквадратическим отклонени-

ем σ_x, а погрешность измерений имеет нормальное распределение с нулевым

математическим ожиданием и среднеквадратическим отклонением σ_y:

w_k ∼ N(0,σ_x), v_k ∼ N(0,σ_y).

Построим фильтр Калмана, а также воспользуемся гарантирующим под-

ходом, объединив помехи w_k и v_k в общий вектор возмущений

)

⁽w_k

(16)

w_k =

v_k

При этом матрицы D₁ и D₂ в системе (1) примут вид

√

(

)

√

(

)

D₁ = 3σ_x

G 0

D₂ = 3σ_y

что позволит варьировать величины w_k и v_k независимо друг от друга на

интервалах

|w_k| ≤ 3σ_x,

|v_k| ≤ 3σ_y.

Гарантирующий подход привел к оптимальным матрицам фильтра

)

⁽0,2359

L∗1 =

0,1412

2,0

Траектория

Наблюдение

ГФ: оценка

1,5

ГФ: оценка

ФК: оценка

ГФ: границы

1,0

0,5

0.5

1,0

1,5

2,0

t, сек

Рис. 1. Динамика координат и их оценок в примере 1 (модель M₁).

для координаты x⁽¹⁾ = s и

)

⁽0,1122

L∗2 =

0,0386

для координаты x⁽²⁾ = s.

Результаты сравнения с фильтром Калмана при

Δt = 0,1, σ_x = 0,1, σ_y = 0,5

показаны на рис. 1,а и 1,б .

2. В рамках модели M₂ будем предполагать, что на k-м такте вагонетка

движется с постоянным ускорением, равномерно распределенным на отрезке

[-a, a], а погрешность измерений равномерно распределена на отрезке [-v, v]:

w_k ∼ U(-a,a), v_k ∼ U(-v,v).

Построим фильтр Калмана (при σ_x = a/3, σ_y = v/3), а также воспользу-

емся гарантирующим подходом — для возмущения (16) и матриц

√

(

)

√

(

)

D₁ = a

G 0

D₂ = v

Гарантирующий подход привел к оптимальным матрицам фильтра

)

⁽0,1393

L∗1 =

0,0489

для координаты x⁽¹⁾ = s и

)

⁽0,0574

L∗2 =

0,0101

для координаты x⁽²⁾ = s.

2,0

Траектория

Наблюдение

ГФ: оценка

1,5

ФК: оценка

ГФ: границы

1,0

0,5

0.5

1,0

1,5

t, сек

Рис. 2. Динамика координат и их оценок в примере 1 (модель M₂).

Результаты сравнения с фильтром Калмана при

Δt = 0,1, a = 0,1, v = 2

показаны на рис. 2,а и 2,б .

3. В рамках модели M₃ будем предполагать, что w_k и v_k принимают про-

извольные значения на отрезках [-a, a] и [-v, v] соответственно:

|w_k| ≤ a,

|v_k| ≤ v.

Применим фильтр Калмана (для σ_x = a, σ_y = v), а также воспользуемся

гарантирующим подходом — для возмущения (16) и матриц

√

(

)

√

(

)

D₁ = a

G 0

D₂ = v

Гарантирующий подход привел к оптимальным матрицам фильтра

)

⁽0,1397

L∗1 =

0,0492

для координаты x⁽¹⁾ = s и

)

⁽0,0574

L∗2 =

0,0101

для координаты x⁽²⁾ = s.

Результаты сравнения с фильтром Калмана при

Δt = 0,1, a = 0,1, v = 2

показаны на рис. 3,а и 3,б .

Траектория

Наблюдение

ГФ: оценка

ФК: оценка

ГФ: границы

t, сек

Рис. 3. Динамика координат и их оценок в примере 1 (модель M₃).

Пример 2. Следующий пример [24] связан с оцениванием движения сна-

ряда по баллистической траектории при наличии внешних возмущений и до-

ступных наблюдению (зашумленных) координатах. Соответствующая систе-

ма имеет вид (15), где

⎞

^⎛s_x

⎜sy⎟

⎜

⎟

⎝vx⎠

v_y

есть вектор состояния системы, состоящий из проекций координаты и скоро-

сти снаряда на горизонтальную и вертикальную оси, а

⎞

⎛

⎞

^⎛1 0 Δt

⎜₀₁

Δt

⎟

⎜

⎟

⎜

⎟

⎜

⎟

⎝001-b

⎠^,

⎝

⎠^,

1-b

-gΔt

(

)

B₁ = G = I, B₂ = 0, C =

Здесь Δt — интервал между измерениями, 0 < b ≪ 1 — коэффициент со-

противления воздуха, g — гравитационная постоянная, шум w_k распределен

нормально с нулевым математическим ожиданием и ковариационной мат-

рицей Q_k ≻ 0, а помехи измерения v_k распределены нормально с нулевым

математическим ожиданием и ковариационной матрицей R_k ≻ 0. При этом

Δt = 0,1, b = 10-4, g = 9,8, Q = 0,1I, R = 500I.

1. В рамках “стандартной” модели M₁ предполагается, что

w_k ∼ N(0,σ_xI), σ2x = 0,1,

v_k ∼ N(0,σ_yI), σ2y = 500.

12 000

700

Траектория

Наблюдение

ГФ: оценка

10 000

ГФ: оценка

600

ФК: оценка

ГФ: границы

500

8000

400

6000

300

4000

200

2000

100

2000

100

t, сек

Рис. 4. Динамика координат и их оценок в примере 2 (модель M₁).

Как и в предыдущем примере, для того чтобы воспользоваться гарантирую-

щим подходом, объединим помехи w_k и v_k в общий вектор возмущений, а

матрицы D₁ и D₂ в системе (1) примут вид

√

(

)

√

(

)

D₁ = 3σ_x

G 0

D₂ = 3σ_y

что позволит варьировать величины w_k и v_k независимо друг от друга на

интервалах

|w_k| ≤ 3σ_x,

|v_k| ≤ 3σ_y.

Гарантирующий подход доставил следующие оптимальные матрицы филь-

тра для каждой из четырех координат (здесь и далее в матрицах фильтра

обнулены элементы, по абсолютной величине не превосходящие 10-6):

⎞

⎛

⎞

^⎛0,5946

0,7277

⎜

0,6822⎟

⎜

0,6340⎟

⎜

⎟

⎜

⎟

L∗1 =

L∗2 =

⎝0,8467

⎠^,

⎝1,2376

⎠^,

1,0590

0,9879

⎞

⎛

⎞

^⎛0,0971

0,1393

⎜

⎟

⎜

⎟

0,1389

0,0975

⎜

⎟

⎜

⎟

L∗3 =

L∗4 =

⎝0,0284

⎠^,

⎝0,0459

⎠^.

0,0456

0,0285

Результаты сравнения для координат s_y и v_y показаны на рис. 4,а и 4,б.

Обратим внимание, что поведение траекторий оценок на начальном участ-

ке обусловлено тем, что фактическая начальная точка

⎛

⎞

⎜₀

⎟

⎜

⎟

x₀ =

⎝₃₀₀⎠

600

8000

800

Траектория

7000

Наблюдение

700

ГФ: оценка

ФК: оценка

ГФ: границы

6000

ГФ: границы

600

5000

500

4000

400

3000

300

2000

200

1000

100

1000

100

t, сек

Рис. 5. Динамика координат и их оценок в примере 2 (модель M₂).

14 000

700

Траектория

Наблюдение

ГФ: оценка

12 000

600

ГФ: оценка

ФК: оценка

ГФ: границы

10 000

ГФ: границы

500

8000

400

6000

300

4000

200

2000

100

2000

100

t, сек

Рис. 6. Динамика координат и их оценок в примере 2 (модель M₃).

оказалась вне минимального инвариантного эллипсоида для невязки. Однако

в силу того, что инвариантный эллипсоид обладает свойством притягиваемо-

сти, после нескольких шагов гарантирующие оценки начинают охватывать

истинную траекторию.

2. В рамках модели M₂ будем предполагать, что

w_k ∼ U(-w,w), v_k ∼ U(-v,v),

где w = 3σ_x, v = 3σ_y.

Поскольку в этом случае

√

(

)

√

(

)

D₁ = w

G 0

D₂ = v

то гарантирующий подход приводит к тем же матрицам фильтра, что и в

модели M₂. Результаты его сравнения с фильтром Калмана для координат s_y

и v_y показаны на рис. 5,а и 5,б.

3. В рамках модели M₃ будем предполагать, что

|w_k| ≤ σ_x,

|v_k| ≤ σ_y.

Гарантирующий подход при

√

(

)

√

(

)

D₁ = σ_x

G 0

D₂ = σ_y

привел к оптимальным матрицам фильтра

⎞

⎛

⎞

^⎛0,6362

0,6697

⎜

⎟

⎜

⎟

0,7256

0,5747

L∗1 =

⎜

⎟

L∗2 =

⎜

⎟

⎝0,9779

⎠^,

⎝1,0517

⎠^,

1,2126

0,8206

⎞

⎛

⎞

^⎛0,0971

0,1393

⎜

0,1389⎟

⎜

0,0975⎟

⎜

⎟

⎜

⎟

L∗3 =

L∗4 =

⎝0,0284

⎠^,

⎝0,0459

⎠^.

0,0456

0,0285

Результаты его сравнения с фильтром Калмана для координат s_y и v_y пока-

заны на рис. 6,а и 6,б .

Как видно из примеров, и фильтр Калмана, и гарантирующий фильтр,

с одной стороны, не слишком сильно отличаются по своим результатам,

а с другой — во всех трех моделях они вполне работоспособны. При этом, как

ожидалось, для гауссовских помех фильтр Калмана дает немного (но не кар-

динально) лучшее поведение оценки по сравнению с гарантирующим филь-

тром, тогда как для неслучайных ограниченных помех преимущество оста-

ется за гарантирующим фильтром.

Отметим также, что ширина трубки, в которую гарантированно заключе-

на соответствующая оценка, довольно сильно завышена, что хорошо видно

на рисунках; это — характерное поведение методов гарантированного оцени-

вания, направленных на противодействие “наихудшей” из возможных реали-

заций неопределенности.

Наконец, гарантирующий фильтр позволяет получить равномерную оцен-

ку точности фильтрации; здесь уместно обратить внимание на поведение га-

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

этом последняя также имеет явно выраженный всплеск).

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

Предложен новый подход к задаче синтеза гарантирующего фильтра, ос-

нованный на сведении проблемы к задаче матричной оптимизации, где пе-

ременной является матрица фильтра. Далее эта задача решается градиент-

ным методом; его сходимость теоретически обосновывается для ряда важных

частных случаев.

Рассмотренные примеры демонстрируют работоспособность и эффектив-

ность предлагаемого алгоритма; проведено его сравнение с фильтром Кал-

мана на трех различных постановках задач.

В статье рассматривалась задача построения стационарного фильтра;

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

задачу фильтрации в духе работ [25, 26] с помощью инструментария [27, 28].

ПРИЛОЖЕНИЕ 1

Лемма П.1.1. Пусть X и Y — решения двойственных дискретных урав-

нений Ляпунова с шуровской матрицей A:

A^TXA - X + W = 0 и AY A^T - Y + V = 0.

Тогда tr (XV ) = tr (Y W ).

Доказательство леммы П.1.1. В самом деле, прямым вычислением

имеем

(

)

tr (XV ) = tr

X(Y - AY A^T)

= tr (XY ) - tr (XAY A^T) =

(

)

= tr (Y X) - tr (Y A^TXA) = tr

Y (X - A^TXA)

= tr (Y W ).

Лемма П.1.1 доказана.

Лемма П.1.2. Для решения P дискретного уравнения Ляпунова

AP A^T - P + Q = 0

с шуровской матрицей A и Q ≻ 0 справедливы оценки:

λ_min(Q)

(П.1.1)

λ_max(P) ≥

λ_min(P) ≥

1-ρ²

1 - σ^2min(A)

где ρ = max|λi(A)|, а σmin(A) — наименьшее сингулярное число матрицы A.

Если же Q = DD^T и пара (A, D) управляема, то

∥u^∗D∥

(П.1.2)

λ_max(P) ≥

> 0,

1-ρ²

где

u^∗A = λu^∗,

|λ| = ρ,

∥u∥ = 1,

т.е. u — левый собственный вектор матрицы A, отвечающий собствен-

ному значению λ матрицы A с наибольшим модулем. Вектор u и число λ

могут быть комплексными; здесь u^∗ означает комплексное сопряжение и

транспонирование.

Доказательство леммы П.1.2. Оценки (П.1.1) хорошо известны, см.,

например, [29]. Докажем справедливость оценки (П.1.2). Явное решение дис-

кретного уравнения Ляпунова для шуровской матрицы A имеет вид

∑

P = A^kDD^T(A^T)^k.

k=0

Умножая это равенство справа на u и слева на u^∗ и учитывая, что u^∗A^k =

= λ^ku^∗, (A^T)^ku = (λ^∗)^ku, получаем, что

∑

u^∗D∥²

λ_max(P) ≥ u^∗Pu = u^∗A^kDD^T(A^T)^ku = (λλ^∗)^ku^∗DD^Tu =^∥

1-ρ²

k=0

причем ∥u^∗D∥ > 0 в силу управляемости пары (A, D), см., например, [10, тео-

рема Д.1.5]. Лемма П.1.2 доказана.

Итак, начнем с оптимизации функции f(α) и рассмотрим оптимизацион-

ную задачу

min f(α), f(α) = tr CP C^T

при ограничении

AP A^T - P +

DD^T = 0

1-α

относительно матричных переменных P = P^T ∈ R^n×n и скалярного парамет-

ра 0 < α < 1.

Здесь наложены более жесткие требования к постановке задачи, чем обыч-

но: будем предполагать, что матрица C выхода системы — квадратная невы-

рожденная. Это предположение можно было бы ослабить, но цель сейчас —

получить наиболее простые и наглядные результаты.

Лемма П.1.3. Пусть матрица A шуровская, ρ — спектральный радиус

матрицы A, ρ² < α < 1, пара (A, D) управляема, а матрица C такова, что

C^TC ≻ 0. Тогда функция f(α) = tr CP(α)C^T обладает следующими свой-

ствами:

а) функция f(α) определена, положительна и сильно выпукла на интер-

вале ρ² < α < 1, а ее значения стремятся к бесконечности на концах ин-

тервала, причем существует c > 0 такое, что

(П.1.3)

f (α) ≥

c, ρ²

< α < 1;

(1 - α)(α - ρ²)

б) производная функции f(α) имеет вид

(

)

f^′(α) = tr Y

DD^T -

AP A^T

(1 - α)²

α²

где P и Y — решения дискретных уравнений Ляпунова

(П.1.4)

AP A^T - P +

DD^T

1-α

(П.1.5)

A^TY A - Y + C^T

C = 0.

в) вторая производная функции f(α) определяется формулой

(

)

f^′′(α) = 2tr Y

DD^T +

A(P - X)A^T

(1 - α)³

α³

где P , Y и X — решения дискретных уравнений Ляпунова (П.1.4) и (П.1.5)

AXA^T - X +

DD^T -

AP A^T = 0,

(1 - α)²

α²

причем f^′′(α^∗) > 0 и f^′′(α) монотонно возрастает слева и справа от α^∗.

Доказательство леммы П.1.3.

а. Уравнение (6) представимо в виде

(

)

(

)_T

-P =-

DD^T

√αAP

√αA

1-α

и согласно [10, лемма 1.2.6] имеет единственное решение тогда и только тогда,

когда матрица

√

A шуровская: |λ_i

√αA)|<1,т.е.приρ2 <α<1.

Оценим величину f(α) = tr CP (α)C^T, используя лемму П.1.2 с очевидны-

ми заменами:

f (α) = tr CP (α)C^T ≥ λ_min(C^TC)λ_max (P (α)) ≥

∥u^∗D∥²λ_min(C^TC)

≥

∥u^∗D∥²λ_min(C^TC),

(1 - α) (1 - ρ²(A/√α))=

(1 - α)(α - ρ²)

где u имеет тот же смысл, что и в лемме П.1.2, а величина ∥u^∗D∥² положи-

тельна в силу предположения об управляемости пары (A, D), а тем самым и

пары (A/√α,D).

Покажем теперь, что функция f(α) = tr CP (α)C^T строго выпукла на ин-

тервале (ρ², 1). В соответствии с [10, лемма 1.2.6] решение уравнения (П.1.4)

представимо в явном виде как

(

)_k

(

)_k

∑

P (α) =

DD^T

A^T

A^kDD^T(A^T)^k

√αA

1-α

√α

(1 - α)α^k

k=0

H_k

g(α,k)

Но H_k ≻ 0 и g(α, k) > 0 при 0 < α < 1, поэтому на интервале (ρ², 1) имеем

∑

P (α) = g(α, k)H_k ≻ 0

k=0

f (α) = tr P (α)C^TC > 0.

Прямым вычислением получаем

(

)

g^′(α,k) =

g(α, k),

1-α

((

)₂

)

g^′′(α,k) =

g(α, k) ≥

g(α, k)

1-α

(1 - α)²

α²

(1 - α)²

(здесь дифференцирование производится по α), так что

∑

f^′′(α) =

g^′′(α,k)tr CH_kC^T ≥

f (α) ≥

f (α^∗) > 0.

(1 - α)²

(1 - ρ²)²

k=0

Таким образом, вторая производная функции f(α) положительна и стре-

мится к бесконечности на концах интервала (ρ², 1).

Далее, прямым вычислением четвертой производной получаем

∑

k(k + 1)(k + 2)(k + 3)

g(IV)(α,k) =

≥

αk+4

(1 - α)⁴

k=0

так что

∑

f(IV)(α) =

g(IV)(α,k)tr CH_kC^T ≥

k=0

∑

≥

tr CH_kC^T >

tr CH_kC^T > 0,

(1 - α)⁴

(1 - ρ²)⁴

k=0

т.е. вторая производная f^′′(α) сама является выпуклой и растет на границах

интервала.

б. Выведем теперь формулу для производной функции f(α). В уравне-

нии (П.1.4) решение P является функцией от α. Продифференцируем это

уравнение; под P^′ будем понимать производную по α:

(П.1.6)

AP^′A^T - P^′ +

DD^T -

AP A^T

= 0.

(1 - α)²

α²

Применяя лемму П.1.1 к двойственным уравнениям (П.1.6) и (П.1.5), полу-

чаем желаемую формулу

(

)

f^′(α) = tr CP^′C^T = tr P^′C^TC = tr Y

DD^T -

AP A^T

(1 - α)²

α²

в. Аналогично получим выражение для второй производной f(α). Диффе-

ренцируя уравнение (П.1.6) по α, получаем

AP^′′A^T - P^′′ +

DD^T +

AP A^T -

AP^′A^T = 0.

(1 - α)³

α³

Вновь применяя лемму П.1.1 к этому уравнению и уравнению (П.1.5) (и имея

в виду, что X = P^′), получаем

(

)

f^′′(α) = tr CP^′′C^T = tr P^′′C^TC = 2tr Y

DD^T +

A(P - X)A^T

(1 - α)³

α³

Лемма П.1.3 доказана.

Заметим, что для вычисления функции f(α) и двух ее производных до-

статочно решить три дискретных уравнения Ляпунова.

Указанные свойства функции позволяют применить метод Ньютона для ее

минимизации. Задаемся начальным приближением ρ²(A) < α₀ < 1, например

(

)

α₀ =

1 + ρ²(A)

/2, и применяем итерационный процесс

f^′(α_j)

(П.1.7)

αj+1 = α_j -

f^′′(α_j)

Следующая теорема гарантирует глобальную сходимость алгоритма; ее

справедливость устанавливается аналогично сходному результату в [16].

Теорема П.1.1

[16]. В методе (П.1.7) справедливы оценки

f^′′(α₀)

|α_j - α^∗| ≤

|α₀ - α^∗|,

|αj+1 - α^∗| ≤ c|α_j - α^∗|²,

2j f^′′(α^∗)

где c > 0 — некоторая константа (она может быть выписана явно).

Первая оценка гарантирует глобальную сходимость метода (быстрее, чем

геометрическая прогрессия с коэффициентом 1/2), а вторая — квадратич-

ную сходимость в окрестности решения. Реально требуется не более трех -

четырех итераций для получения решения с большой точностью (если только

начальная точка не слишком близка к границам интервала).

Вернемся к оптимизационной задаче (4)-(5) и займемся минимизацией

функции

f (L) = min f(L, α),

предварительно исследовав ее свойства.

Лемма П.1.4. Функция f(L) определена и положительна на множе-

стве S допустимых матриц фильтра.

Действительно, если матрица A - LC шуровская, то ρ(A - LC) < 1 и для

ρ²(A - LC) < α < 1 существует решение P ≽ 0 дискретного уравнения Ля-

пунова (5). Тем самым определена строго положительная функция f(L, α),

при этом f(L) > 0 в силу (П.1.3). Как и в непрерывном случае, множество

ее определения S может быть невыпуклым и несвязным, а его границы —

негладкими.

Лемма П.1.5. На множестве S функция f(L) коэрцитивна (т.е. стре-

мится к бесконечности на границе области), причем справедливы следующие

оценки:

λ_min(C₁C^T1)

(П.1.8)

f (L) ≥

∥D₁ - LD₂∥2F ,

1 - ρ²(A - LC) 1 - σ^2min(A - LC)

f (L) ≥ ρ∥L∥².

Доказательство леммы П.1.5. Рассмотрим последовательность до-

пустимых матриц {L_j } ⊆ S такую, что L_j → L ∈ ∂S, т.е. ρ(A - LC) = 1. Это

означает, что для любого ε > 0 найдется число N = N(ε) такое, что неравен-

ство

|ρ(A - L_j C) - ρ(A - LC)| = 1 - ρ(A - L_jC) < ε

справедливо для всех j ≥ N(ε).

Пусть P_j — решение уравнения (5), ассоциированного с матрицей филь-

тра L_j :

(A - L_jC)P_j (A - L_jC)^T - P_j +

(D₁ - L_j D₂)(D₁ - L_jD₂)^T = 0,

α_j

1-α_j

а Y_j — решение двойственного к нему дискретного уравнения Ляпунова

(A - L_jC)^TY_j(A - L_jC) - Y_j + C₁C^T1 = 0.

α_j

Тогда с учетом леммы П.1.2 имеем:

(

)

f (L_j) = tr (C₁P_j C^T1) + ρ∥L_j ∥2F ≥ tr

P_jC₁C^T1

(

)

= tr Y_j

(D₁ - L_jD₂)(D₁ - L_j D₂)T

≥

1-α_j

λ_min(C₁C^T1)

≥

λ_min(Y_j)∥D₁ - L_jD₂∥2F ≥

∥D₁ - L_jD₂∥2F ≥

1-α_j

1 - α_j 1 - σ^2min(A - L_jC)

λ_min(C₁C^T1)

≥

∥D₁ - L_jD₂∥2F ≥

1 - ρ²(A - L_jC) 1 - σ^2min(A - L_jC)

λ_min(C₁C^T1)

≥

∥D₁ - L_jD₂∥2F

-−-→ +∞,

ε 1 - σ^2min(A - L_jC)

ε→0

поскольку ρ²(A - L_jC) < α_j < 1.

C другой стороны,

f (L_j) = tr (C₁P_j C^T1) + ρ∥L_j ∥2F ≥ ρ∥L_j ∥2F ≥ ρ∥L_j ∥²

-−-----→ +∞.

∥L_j ∥→+∞

Лемма П.1.5 доказана.

Введем в рассмотрение множество уровня

S₀ = {L ∈ S : f(L) ≤ f(L₀)}.

Из леммы П.1.5 вытекает очевидное

Следствие П.1.1. Для любых L₀ ∈ S множество S₀ ограничено.

C другой стороны, у функции f(L) на множестве S₀ существует точка

минимума (как у непрерывной — в силу свойств решения дискретного урав-

нения Ляпунова — функции на компактном множестве), но множество S₀ не

имеет общих точек с границей S в силу (П.1.8). Далее будет показано, что

f (L) дифференцируема на S₀. Следовательно, справедливо

Следствие П.1.2. Существует точка минимума L_∗ на множестве S,

и в ней градиент функции f(L) обращается в нуль.

Перейдем к свойствам градиента функции f(L, α).

Лемма П.1.6. Функция f(L,α) определена на множестве стабилизи-

рующих L и для ρ²(A - LC) < α < 1. На этом допустимом множестве она

дифференцируема, причем градиент дается выражениями

(

(П.1.9)

∇_αf(L,α) = tr Y

(D₁ - LD₂)(D₁ - LD₂)^T -

(1 - α)²

)

(A - LC)P (A - LC)T ,

α²

(

)

(П.1.10)

∇_Lf(L,α) = 2 ρL -

Y (A - LC)P C^T -

Y (D₁ - LD₂)D^T

1-α

где матрицы P и Y являются решениями дискретных уравнений Ляпуно-

ва (5) и (8).

Минимум f(L,α) достигается во внутренней точке допустимого мно-

жества и определяется условиями

∇_Lf(L,α) = 0,

∇_αf(L,α) = 0.

При этом f(L, α) как функция от α строго выпукла на ρ²(A - LC) < α < 1

и достигает минимума во внутренней точке этого интервала.

Доказательство леммы П.1.6. Имеем задачу:

min f(L, α), f(L, α) = tr C₁P C^T1 + ρ∥L∥2F

при ограничении в виде дискретного уравнения Ляпунова (5) относительно

матрицы P инвариантного эллипсоида.

Следуя лемме П.1.3, дифференцирование по α производится в соответ-

ствии с соотношениями (П.1.9), (5) и (8). Для дифференцирования по L да-

дим ему приращение ΔL и обозначим соответствующее приращение P

че-

рез ΔP ; в результате (5) примет вид

(A - (L + ΔL)C) (P + ΔP ) (A - (L + ΔL)C)^T - (P + ΔP )+

(D₁ - (L + ΔL)D₂) (D₁ - (L + ΔL)D₂)^T = 0.

1-α

Оставляя обозначение ΔP для главной части приращения, получаем

(

(A - LC)P (A - LC)^T - ΔLCP (A - LC)^T - (A - LC)P (ΔLC)^T+

)

+(A - LC)ΔP (A - LC)^T

- (P + ΔP )+

(

(D₁ - LD₂)(D₁ - LD₂)^T - ΔLD₂(D₁ - LD₂)^T -

1-α

)

- (D₁ - LD₂)(ΔLD₂)^T

= 0.

После вычитания уравнения (12) из этого уравнения, имеем:

(П.1.11)

(A - LC)ΔP (A - LC)^T -

(

)

- ΔP -

ΔLCP(A - LC)^T + (A - LC)P(ΔLC)^T

(

)

−

ΔLD₂(D₁ - LD₂)^T + (D₁ - LD₂)(ΔLD₂)^T

= 0.

1-α

Вычислим приращение функционала f(L), линеаризуя соответствующие

величины:

Δf(L) = tr C₁ΔPC^T1 + ρtr L^TΔL + ρtr (ΔL)^TL = tr ΔPC^T1C₁ + 2ρtr L^TΔL.

По лемме П.2.1 из двойственных уравнений (П.1.11) и (8) имеем

⁽1

Δf(L) = -2tr Y

ΔLCP(A - LC)^T +

)

ΔLD₂(D₁ - LD₂)^T

+ 2ρ tr L^TΔL =

1-α

(

)

= 2tr ρL^TΔL -

CP(A - LC)^TY -

D₂(D₁ - LD₂)^TY ΔL =

1-α

9 (

)

ρL -

Y (A - LC)P C^T -

Y (D₁ - LD₂)D^T

,ΔL

1-α

Таким образом, приходим к соотношению (П.1.10). Лемма П.1.6 доказана.

Градиент функции f(L) не является липшицевым на множестве S, однако

можно показать, что он обладает этим свойством на его подмножестве S₀

(аналогично тому, как это было сделано в [16]).

Установленные свойства минимизируемой функции и ее производных обос-

новывают метод минимизации, реализованный в виде Алгоритма 1.

ПРИЛОЖЕНИЕ 2

Лемма П.2.1

[16]. Пусть X и Y — решения двойственных уравнений

Ляпунова с гурвицевой матрицей A:

A^TX + XA + W = 0 и AY + Y A^T + V = 0.

Тогда tr (XV ) = tr (Y W ).

Свойства функции f(α), установленные в [16], полностью переносятся на

рассматриваемый случай. В частности, функция f(α) определена, положи-

тельна и сильно выпукла на интервале 0 < α < 2σ(A - LC), а ее значения

стремятся к бесконечности на концах интервала, причем существует c > 0

такое, что

(П.2.1)

f (α) ≥

0 < α < 2σ(A - LC).

α(2σ - α)

Минимизацию функции f(α) можно эффективно осуществлять при по-

мощи метода Ньютона. Зададимся начальным приближением

0<α₀ <

< 2σ(A - LC), например α₀ = σ(A - LC), и применим итерационный процесс

f^′(α_j)

αj+1 = α_j -

f^′′(α_j)

где согласно [16]

(

)

(П.2.2)

f^′(α) = tr Y P -

(D₁ - LD₂)(D₁ - LD₂)T

α²

(

)

f^′′(α) = 2tr Y X +

(D₁ - LD₂)(D₁ - LD₂)T

α³

а P, Y и X — решения уравнений Ляпунова (12), (13) и (14). При этом тео-

рема П.1.1 сохраняет свою силу.

Следующая лемма является непрерывным аналогом леммы П.1.4.

Лемма П.2.2. Функция f(L) определена и положительна на множе-

стве S допустимых матриц фильтра.

Действительно, если матрица A - LC гурвицева, то σ(A - LC) > 0 и для

0 < α < 2σ(A - LC) существует решение P ≽ 0 уравнения Ляпунова (12).

Тем самым определена строго положительная функция f(L, α), при этом

f (L) > 0 в силу (П.2.1). Множество ее определения S может быть невыпук-

лым и несвязным, причем его границы могут быть негладкими, см. [16].

Лемма П.2.3. На множестве S допустимых матриц функция f(L) ко-

эрцитивна (т.е. стремится к бесконечности на границе области), причем

справедливы следующие оценки:

λ_min(C₁C^T1)∥D₁ - LD₂∥2F

(П.2.3)

f (L) ≥

4σ(A - LC) (∥A - LC∥ + σ(A - LC))

f (L) ≥ ρ∥L∥².

Доказательство леммы П.2.3. Рассмотрим последовательность до-

пустимых матриц {L_j } ⊆ S такую, что L_j → L ∈ ∂S, т.е. σ(A - LC) = 0. Это

означает, что для любого ε > 0 найдется число N = N(ε) такое, что неравен-

ство

|σ(A - L_jC) - σ(A - LC)| = σ(A - L_jC) < ε

справедливо для всех j ≥ N(ε).

Пусть P_j — решение уравнения (12), ассоциированного с матрицей филь-

тра L_j :

(

)

(

)_T

α_j

A-L_jC+

I P_j +P_j A-L_jC+

(D₁ - L_jD₂)(D₁ - L_jD₂)^T = 0,

α_j

а Y_j — решение двойственного к нему уравнения Ляпунова

(

)_T

(

)

α_j

A-L_jC+

Y_j + Y_j A - L_jC +

+ C₁C^T1 = 0.

Тогда с учетом [16, лемма П.3] имеем:

f (L_j) = tr (C₁P_j C^T1) + ρ∥L_j ∥2F ≥

(

)

≥ tr (P_j C₁C^T) = tr Yj (D1 - LjD2)(D1 - LjD2)T

≥

α_j

λ_min(C₁C^T1)

≥

λ_min(Y_j)∥D₁ - L_jD₂∥2F ≥

α_j

α_j 2∥A - L_jC +αj2I∥∥D1 -LjD²∥^F ≥

λ_min(C₁C^T1)

≥

4σ(A - L_jC) ∥A - L_jC +αj2I∥∥D1 -Lj D²∥^F ≥

λ_min(C₁C^T1)

≥

∥D₁ - L_jD₂∥2F

-−-→ +∞,

4ε(∥A - L_j C∥ + ε)

ε→0

поскольку 0 < α_j < 2σ(A - L_j C) и

α_j

∥A - L_j C +

I∥ ≤ ∥A - L_jC∥ +

C другой стороны,

f (L_j) = tr (C₁P_j C^T1) + ρ∥L_j ∥2F ≥ ρ∥L_j ∥2F ≥ ρ∥L_j ∥²

-−-----→ +∞.

∥L_j ∥→+∞

Лемма П.2.3 доказана.

Введем в рассмотрение множество уровня

S₀ = {L ∈ S : f(L) ≤ f(L₀)}.

Из леммы П.2.3 вытекает очевидное

Следствие П.2.3. Для любых L₀ ∈ S множество S₀ ограничено.

C другой стороны, у функции f(L) на множестве S₀ существует точка ми-

нимума (как у непрерывной — в силу свойств решения уравнения Ляпунова —

функции на компактном множестве), но множество S₀ не имеет общих точек с

границей S в силу (П.2.3). Далее будет показано, что f(L) дифференцируема

на S₀. Следовательно, справедливо

Следствие П.2.4. Существует точка минимума L_∗ на множестве S,

и в ней градиент функции f(L) обращается в нуль.

Перейдем к свойствам градиента функции f(L, α).

Лемма П.2.4. Функция f(L,α) определена на множестве допустимых L

и для 0 < α < 2σ(A - LC). На этом допустимом множестве она дифферен-

цируема, причем градиент дается выражениями

(

)

∇_αf(L,α) = tr Y P -

(D₁ - LD₂)(D₁ - LD₂)T

α²

(

)

(П.2.4)

∇_Lf(L,α) = 2 ρL - Y PC^T -

Y (D₁ - LD₂)D^T

где матрицы P и Y являются решениями уравнений Ляпунова (12) и (13).

Минимум f(L,α) достигается во внутренней точке допустимого мно-

жества и определяется условиями

∇_Lf(L,α) = 0,

∇_αf(L,α) = 0.

При этом f(L, α) как функция от α строго выпукла на 0 < α < 2σ(A - LC)

и достигает минимума во внутренней точке этого интервала.

Доказательство леммы П.2.4. Имеем задачу:

min f(L, α), f(L, α) = tr C₁P C^T1 + ρ∥L∥2F

при ограничении в виде уравнения Ляпунова (12) относительно матрицы P

инвариантного эллипсоида.

Дифференцирование по α производится в соответствии с соотношения-

ми (П.2.2), (12) и (13). Для дифференцирования по L дадим ему приращение

ΔL и обозначим соответствующее приращение P через ΔP; в результате (12)

примет вид

(

)

(

)_T

A - (L + ΔL)C +

I (P + ΔP ) + (P + ΔP ) A - (L + ΔL)C +

(D₁ - (L + ΔL)D₂) (D₁ - (L + ΔL)D₂)^T = 0.

Оставляя обозначение ΔP для главной части приращения, получаем

(

)

(

)_T

A - (L + ΔL)C +

I P + P A - (L + ΔL)C +

(

)

(

)_T

+ A - LC +

I ΔP + ΔP A - LC +

(

(D₁ - LD₂)(D₁ - LD₂)^T - ΔLD₂(D₁ - LD₂)^T -

)

- (D₁ - LD₂)(ΔLD₂)^T

= 0.

После вычитания уравнения (12) из этого уравнения, имеем:

(

)

(

)_T

(П.2.5)

A - LC +

I ΔP + ΔP A - LC +

(

)

−ΔLCP - P(ΔLC)^T -

ΔLD₂(D₁ - LD₂)^T + (D₁ - LD₂)(ΔLD₂)^T

= 0.

Вычислим приращение функционала f(L), линеаризуя соответствующие

величины:

Δf(L) = tr C₁ΔPC^T1 + ρtr L^TΔL + ρtr (ΔL)^TL = tr ΔPC^T1C₁ + 2ρtr L^TΔL.

По лемме П.2.1 из двойственных уравнений (П.2.5) и (13) имеем

(

)

Δf(L) = - tr 2Y ΔLCP +

ΔLD₂(D₁ - LD₂)T

+ 2ρ tr L^TΔL =

(

)

= 2tr ρL^TΔL - CPY ΔL -

D₂(D₁ - LD₂)^TY ΔL

9 (

)

ρL - Y P C^T -

Y (D₁ - LD₂)D^T

,ΔL

Таким образом, приходим к соотношению (П.2.4). Лемма П.2.4 доказана.

СПИСОК ЛИТЕРАТУРЫ

1. Kalman R.E. A New Approach to Linear Filtering and Prediction Problems // J.

Basic Engineer. 1960. V. 82. I. 1. P. 35-45.

2. Kailath T., Sayed A.H., Hassibi B. Linear Estimation. N.J.: Prentice Hall, 2000.

Матасов А.И. Основы теории фильтра Калмана. М.: Изд-во МГУ, 2021.

Schweppe F.C. Uncertain Dynamic Systems. N.J.: Prentice Hall, 1973.

Куржанский А.Б. Управление и наблюдение в условиях неопределенности. М.:

Наука, 1977.

Черноусько Ф.Л. Оценивание фазового состояния динамических систем. М.: На-

ука, 1988.

Поляк Б.Т., Топунов М.В. Фильтрация при неслучайных возмущениях: метод

инвариантных эллипсоидов // Докл. РАН. 2008. Т. 418. № 6. С. 749-753.

Polyak B.T., Topunov M.V. Filtering under Nonrandom Disturbances: The Method

of Invariant Ellipsoids // Dokl. Math. 2008. V. 77. No. 1. P. 158-162.

Хлебников М.В., Поляк Б.Т. Фильтрация при произвольных ограниченных

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

13-я

Мультиконференция по проблемам управления (МКПУ-2020). Матер. XXXII

Конференции памяти выдающегося конструктора гироскопических приборов

Н.Н. Острякова. Санкт-Петербург, 6-8 октября 2020 г. СПб.: Концерн “ЦНИИ

“Электроприбор”, 2020. С. 291-294.

Boyd S., El Ghaoui L., Feron E., Balakrishnan V. Linear Matrix Inequalities in

System and Control Theory. Philadelphia: SIAM, 1994.

10.

Поляк Б.Т., Хлебников М.В., Щербаков П.С. Управление линейными система-

ми при внешних возмущениях: Техника линейных матричных неравенств. М.:

ЛЕНАНД, 2014.

11.

Fazel M., Ge R., Kakade S., Mesbahi M. Global Convergence of Policy Gradient

Methods for the Linear Quadratic Regulator // Proc. 35th Int. Conf. Machine

Learning. Stockholm, Sweden, July 10-15, 2018. V. 80. P. 1467-1476.

12.

Mohammadi H., Zare A., Soltanolkotabi M., Jovanović M.R. Global Exponential

Convergence of Gradient Methods Over the Nonconvex Landscape of the Linear

Quadratic Regulator // Proc. 2019 IEEE 58th Conf. Decision Control. Nice, France,

December 11-13, 2019. P. 7474-7479.

13.

Zhang K., Hu B., Basar T. Policy Optimization for H₂ Linear Control with H_∞

Robustness Guarantee: Implicit Regularization and Global Convergence // Proc. 2nd

Conference on Learning for Dynamics and Control (2nd L4DC). Zürich, Switzerland,

June 11-12, 2020. P. 179-190.

14.

Bu J., Mesbahi A., Fazel M., Mesbahi M. LQR through the Lens of First Order

Methods: Discrete-Time Case // arXiv:1907.08921, 2019.

15.

Fatkhullin I., Polyak B. Optimizing Static Linear Feedback: Gradient Method //

SIAM J. Control Optim. 2021. V. 59. No. 5. P. 3887-3911.

16.

Поляк Б.Т., Хлебников М.В. Синтез статического регулятора для подавления

внешних возмущений как задача оптимизации // АиТ. 2021. № 9. С. 86-115.

Polyak B.T., Khlebnikov M.V. Static Controller Synthesis for Peak-to-Peak Gain

Minimization as an Optimization Problem // Autom. Remote Control. 2021. V. 82.

No. 9. P. 1530-1553.

17.

Поляк Б.Т., Хлебников М.В. Синтез обратной связи по выходу при помощи на-

блюдателя как задача оптимизации // АиТ. 2022. № 3. С. 7-32.

Polyak B.T., Khlebnikov M.V. Observer-Aided Output Feedback Synthesis as an

Optimization Problem // Autom. Remote Control. 2022. V. 83. No. 3. P. 303-324.

18.

Поляк Б.Т., Хлебников М.В. Новые критерии настройки ПИД-регуляторов //

АиТ. 2022. № 11. С. 62-82.

Polyak B.T., Khlebnikov M.V. New Criteria for Tuning PID Controllers // Autom.

Remote Control. 2022. V. 83. No. 11. P. 1724-1741.

19.

Luenberger D.G. Observing the State of a Linear System // IEEE Transactions on

Military Electronics. 1964. V. 8. P. 74-80.

20.

Luenberger D.G. An Introduction to Observers // IEEE Trans. Autom. Control.

1971. V. 35. P. 596-602.

21.

Поляк Б.Т., Хлебников М.В., Щербаков П.С. Линейные матричные неравенства

в системах управления с неопределенностью // АиТ. 2021. № 1. С. 3-54.

Polyak B.T., Khlebnikov M.V., Shcherbakov P.S. Linear Matrix Inequalities in

Control Systems with Uncertainty // Autom. Remote Control. 2021. V. 82. No. 1.

P. 1-40.

22.

Назин С.А., Поляк Б.Т., Топунов М.В. Подавление ограниченных внешних воз-

мущений с помощью метода инвариантных эллипсоидов // АиТ. 2007. № 3.

С. 106-125.

Nazin S.A., Polyak B.T., Topunov M.V. Rejection of Bounded Exogenous Distur-

bances by the Method of Invariant Ellipsoids // Autom. Remote Control.

2007.

V. 68. No. 3. P. 467-486.

23.

URL en.wikipedia.org/wiki/Kalman_filter

24.

Humpherys J., Redd P., West J. A Fresh Look at the Kalman Filter // SIAM Rev.

2012. V. 54. Iss. 4. P. 801-823.

25.

Tang W., Zhang Q., Wang Z., Shen Y. Ellipsoid Bundle and its Application to Set-

Membership Estimation // IFAC-PapersOnLine. 2020. V. 53. I. 2. P. 13688-13693.

26.

Tang W., Zhang Q., Wang Z., Shen Y. Set-Membership Filtering with Incomplete

Observations // Inform. Sci. 2020. V. 517. P. 37-51.

27.

Polyak B.T., Nazin S.A., Durieu C., Walter E. Ellipsoidal Parameter or State Esti-

mation under Model Uncertainty // Automatica. 2004. V. 40. I. 7. P. 1171-1179.

28.

Durieu C., Walter E., Polyak B. Multi-Input Multi-Output Ellipsoidal State Bound-

ing // J. Optim. Theory Appl. 2001. V. 111. No. 2. P. 273-303.

29.

Kwon W.H., Moon Y.S., Ahn S.C. Bounds in Algebraic Riccati and Lyapunov Equa-

tions: A Survey and Some New Results // Int. J. Control. 1996. V. 64. P. 377-389.

Статья представлена к публикации членом редколлегии Е.Я. Рубиновичем.

Поступила в редакцию 23.09.2022

После доработки 21.12.2022

Принята к публикации 29.12.2022