CC BY
7
Math-Net.Ru
А. Н. Кувшинова, А. В. Цыганов, Ю. В. Цыганова, Математическое моделирование процесса параметрической идентификации моделей конвективно-диффузионного переноса с применением БУБ-фильтра Калмана, Вестн. Сам. гос. техн. ун-та. Сер. Физ.-мат. науки, 2021, том 25, номер 4, 716— 737
001: https://doi.org/10.14498/vsgtu1876
Использование Общероссийского математического портала Math-Net.Ru подразумевает, что вы прочитали и согласны с пользовательским соглашением
http://www.mathnet.ru/rus/agreement
Параметры загрузки: IP: 109.252.159.147 30 марта 2022 г., 09:16:38
ISSN: 2310-7081 (online), 1991-8615 (print) d https://doi.org/10.14498/vsgtu1876
Математическое моделирование, численные методы и комплексы программ
УДК 519.254, 004.94
Математическое моделирование процесса параметрической идентификации моделей конвективно-диффузионного переноса с применением 8УЮ-фильтра Калмана
А. Н. Кувшинова1, А. В. Цыганов1, Ю. В. Цыганова2
1 Ульяновский государственный педагогический университет им. И. Н. Ульянова, Россия, 432071, Ульяновск, пл. Ленина, 4/5.
2 Ульяновский государственный университет, Россия, 432017, Ульяновск, ул. Льва Толстого, 42.
Аннотация
Рассматривается задача математического моделирования процесса идентификации коэффициентов уравнения в частных производных в моделях конвективно-диффузионного переноса по результатам зашум-ленных измерений значений искомой функции с применением нового метода, относящегося к классу рекуррентных методов параметрической идентификации на основе алгоритмов оптимальной дискретной фильтрации калмановского типа. Рассматриваются одномерные модели с постоянными коэффициентами, граничными условиями первого рода или смешанными граничными условиями первого и третьего рода.
Предлагаемый метод решения задачи основан на переходе от исходной непрерывной модели с уравнением в частных производных к модели, описываемой линейной дискретной динамической системой в пространстве состояний, и применении к ней метода максимального правдоподобия с построением критерия идентификации (функции правдоподобия)
Научная статья
© Коллектив авторов, 2021 © СамГТУ, 2021 (составление, дизайн, макет) 3 ©® Контент публикуется на условиях лицензии Creative Commons Attribution 4.0 International (https://creativecommons.org/licenses/by/4.0/deed.ru)
Образец для цитирования
Кувшинова А. Н., Цыганов А. В., Цыганова Ю. В. Математическое моделирование процесса параметрической идентификации моделей конвективно-диффузионного переноса с применением SVD-фильтра Калмана // Вестн. Сам. гос. техн. ун-та. Сер. Физ.-мат. науки, 2021. Т. 25, №4. С. 716-737. https://doi.org/10.14498/vsgtu1876. Сведения об авторах
Анастасия Николаевна Кувшинова © https://orcid.org/0000-0002-3496-5981 аспирант; каф. высшей математики; e-mail:kuvanulspu@yandex.ru
Андрей Владимирович Цыганов © https://orcid.org/0000-0002-4173-5199 кандидат физико-математических наук, доцент; профессор; каф. высшей математики; e-mail: andrew.tsyganov@gmail.com
Юлия Владимировна Цыганова А https://orcid.org/0000-0001-8812-6035 доктор физико-математических наук, доцент; профессор; каф. информационных технологий; e-mail: tsyganovajv@gmail .com
на основе величин, вычисляемых БУБ-модификацией фильтра Калма-на. Данный фильтр основан на сингулярном разложении ковариационной матрицы ошибок оценивания вектора состояния и устойчиво работает даже в тех случаях, когда она близка к вырожденной. БУБ-фильтр хорошо зарекомендовал себя при решении различных задач дискретной фильтрации и параметрической идентификации и обладает целым рядом преимуществ по сравнению с традиционно используемым стандартным фильтром Калмана, главным из которых является устойчивость к ошибкам машинного округления.
Приводятся результаты компьютерного моделирования процессов параметрической идентификации в системе МЛТЬЛБ с использованием специализированного программного комплекса. Результаты численных экспериментов подтверждают работоспособность предложенного метода и его преимущества по сравнению с аналогичным методом на основе стандартного фильтра Калмана.
Ключевые слова: модель конвективно-диффузионного переноса, параметрическая идентификация, фильтр Калмана, БУБ-фильтр.
Получение: 3 августа 2021 г. / Исправление: 7 декабря 2021 г. / Принятие: 21 декабря 2021 г. / Публикация онлайн: 28 декабря 2021 г.
Введение и постановка задачи. Математические модели, описываемые уравнениями в частных производных, в частности, модели конвективно-диффузионного переноса, широко используются для описания природных и технических процессов [1,2]. Для данных моделей актуальными являются задачи идентификации их параметров по результатам измерений значений искомой функции в отдельных точках рассматриваемой области. Такие задачи относятся к обратными задачам математической физики и в общем случае являются некорректно поставленными [3].
В ряде работ (см., например, [4-9] для решения задач параметрической идентификации моделей, описываемых уравнениями в частных производных, предложено использовать рекуррентные методы, основанные на алгоритмах дискретной фильтрации калмановского типа. Такой подход обладает рядом преимуществ, например для систем, работающих в реальном времени, однако следует заметить, что эффективность реализации рекуррентных методов параметрической идентификации может существенно зависеть от выбора соответствующих алгоритмов оптимальной дискретной фильтрации, поскольку недостатки классического фильтра Калмана, в частности неустойчивость к ошибкам машинного округления, широко известны и описаны в литературе [10]. В связи с этим актуальными являются вопросы разработки алгоритмов параметрической идентификации на основе численно эффективных модификаций фильтра Калмана.
Пусть дана математическая модель конвективно-диффузионного переноса, описываемая следующими уравнениями:
дс дс д2с
(1)
с(х, 0) = <р(х), с(а,1) = / (I),
с(М) = д(1),
(2)
(3)
(4)
или
dc(b,t)
дх
= —\[c(b,t) -9(f)],
(5)
где х £ [a; Ь] — пространственная координата; t £ [0; Т] — время; с(х, t) — искомая функция, например концентрация некоторого вещества в потоке жидкости в точке с координатой х в момент времени t; v — скорость конвекции; а — коэффициент диффузии; (2) — начальное условие; (3), (4), (5) — граничные условия. Таким образом, рассматриваются модели либо с граничными условиями первого рода, либо со смешанными граничными условиями первого и третьего рода.
Рассмотрим задачу параметрической идентификации, состоящую в определении коэффициентов v и а в уравнении (1) по результатам зашумленных измерений значений функции с(х, t) в отдельных точках рассматриваемого отрезка в различные моменты времени (функции ф(х), f (t), g(t) и коэффициентов Л, входящих в начальное и граничные условия предполагаются известными).
Данная задача ранее исследовалась авторами в работах [11-13], однако в них использовались методы идентификации на основе стандартного фильтра Калмана. Основной целью данной статьи является разработка нового метода параметрической идентификации для моделей конвективно-диффузионного переноса на основе SVD-модификации фильтра Калмана и исследование его преимуществ по сравнению со стандартным фильтром с применением методов и средств компьютерного моделирования.
1. Дискретизация исходной модели. Перейдем от исходной непрерывной модели к дискретной модели, описываемой линейной динамической системой в пространстве состояний. Следуя [12,13], зададим в рассматриваемой пространственно-временной области конечно-разностную сетку {(xi,tk) | i = 0,1,...,N, к = 0, 1,...,К}, где
Xi = а + iAx, tfc = к At, Ах =
Ъ — а N - 1:
At
Т
К- 1
Обозначим = с(хг^к), с0 = с(хг, 0) = <р(хг), /к = /), дк = д(Ьк). Заменяя частные производные в уравнении (1) их конечно-разностными аппроксимациями, в случае граничных условий (3), (4) получаем дискретную линейную динамическую систему
„fc
/->fc
сп-2 ph
^п—1 pfc
Ск
Ü2 а1 0
аз (12 а1
0 0 0
0
аз (12
0 0 0
Fk-
0 0 0
(12 a1 0
аз Ü2 a1
0 0 0
0
аз «2J
pfc-1 "2 „к-1
fc-1 Ln-2 fc-1 Ln-1 „к-1
Cfc-l
+
~a1 0
0 0
0 0
0 0
0 0
0 аз.
~в! -l
fk-1 -,k-1
«fc-i
к = 1, 2,..., К, (6)
с
с
з
l
а в случае граничных условий (3), (5) —
.-¡к
сп-2 рк
^п—1 рк
Ск
0,2 аз 0. .0 0 0
а1 (12 аз . .0 0 0
0 а1 й2 . .0 0 0
0 0 0. . й2 аз 0
0 0 0. . а1 02 аз
0 0 0. . 0,40,1 й4й2 04йз
„к-Г -1
1
„к-1
0к-1 "п—2 пк-1 "п— 1 „к-1
Ск-1
+
'а1 0
0 0
0 0
0 0
0 0
0 Оъ.
в^ 1
к-1 к
9
Ик-1
к = 1, 2,..., К. (7)
Системы (6) и (7) являются дискретными линейными динамическими системами, в которых граничные условия входят в двумерный вектор входных воздействий и^. Коэффициенты в матрицах обеих систем имеют следующий вид:
Й1 = П + Г2, 0,2 = 1 - 2Г2, аз = Г2 - Г1,
1
а4
1 + ХАх'
а5
ХАх 1 + ХАх,
где
П =
=
2Аж' 2 Аж2 а сами матрицы являются постоянными (^ = Р £
Вк = В е
ъпх2\
4 = н с\ 4 + ГеП
рк ек
В системе (6) п = N — 1 — вектор состояния состоит из всех внутренних узлов пространственной сетки, а в системе (7) п = N — вектор состояния состоит из всех внутренних узлов пространственной сетки и правой границы.
2. Выбор структуры измерителя. Добавим к уравнениям (6) и (7) модель измерителя в следующем виде:
к = 1, 2,..., К,
где Н £ — матрица измерений, определяющая структуру измерителя;
т — количество измеряемых компонент вектора состояния ; Ск £ — погрешность измерителя. Очевидным выбором является Н = I, где I — единичная матрица. Однако это решение избыточное, поскольку единичная матрица означает наличие п сенсоров для сбора данных измерений.
Обозначим через в неизвестный (в общем случае векторный) параметр линейной динамической системы (6) или (7), подлежащий идентификации. Поставим задачу выбора структуры измерителя с минимальным количеством сенсоров, при которой обеспечивается существование и единственность стационарного оптимального дискретного фильтра Калмана при истинном значении параметра в [14]. Достаточным условием является полная наблюдаемость линейной динамической системы. Данное условие накладывает ограничения
к
с
с
с
з
1
п
на область определения D(ff) параметра в, а также на возможность реализации процесса параметрической идентификации.
Для проверки условия полной наблюдаемости необходимо вычислить ранг матрицы наблюдаемости [15]
^rni = [Нт (HF )т (HF 2)т ... (HFn-1)T]T,
который должен быть равен п (размеру вектора состояния ск).
Следовательно, для решения задачи выбора подходящей структуры измерителя необходимо найти такую матрицу с минимальным количеством сенсоров, для которой rank Жвт1 = п.
В [16] проведен анализ полной наблюдаемости модели (6) и решена задача выбора структуры измерителя для п = 5. Доказано, что при выборе матрицы Н в виде
10000 0 0 0 0 1
Н
модель (6), (8) является полностью наблюдаемой при условии, что коэффициенты а\ и аз в модели (6) не равны нулю. Можно показать, что и в общем случае при выборе матрицы Н в виде
Н
еТ
ет
Л
где г = ^ и строка е1 = (0,..., 0,1,0,..., 0) (единица на г-м месте), модель (6), (8) является полностью наблюдаемой при условии, что коэффициенты а\ и а3 в модели (6) не равны нулю.
Аналогичные результаты получаются для модели (7). Таким образом, для идентификации параметра в в моделях (6) и (7) достаточно всего лишь двух сенсоров для сбора данных измерений.
3. Описание процесса параметрической идентификации. Предположим, что характеристики шума в измерителе известны, а шаги пространственно-временной сетки Ах и заданы, тогда неизвестными параметрами дискретных моделей (6) и (7), подлежащими идентификации, являются скорость конвекции V и коэффициент диффузии а, от которых зависят коэффициенты а\, а,2, Оз, входящие в элементы матриц Р и В.
Дискретные модели (6), (7) с моделью измерителя (8) в общем виде можно записать как
{ ск = ^(в)ск-1 + В(9)ик-г, (9)
\ гк = Нск + &, к = 1,2,..., К, (9)
где д = (у,а)т € К2, ск € Мп — вектор состояния, ик € К2 — вектор входных воздействий, хк € Мт — вектор измерений. Предположим, что погрешность измерителя £к € Мт является гауссовым белым шумом с нулевым математическим ожиданием и ковариационной матрицей К > 0. Матрицы моделей (6) и (7) ^ (в) € Мпхп, В (в) € Мпх2, зависят от параметра в.
Рассмотрим задачу параметрической идентификации модели (9) по доступным данным с целью оценки неизвестного (векторного) параметра в. Задача параметрической идентификации заключается в нахождении неизвестного параметра в по известным входным сигналам ик_1 = [ио,и\,..., ик_1}
и выходным данным измерений Z^ = [z\,z2,..., zk} в соответствии с выбранным критерием качества идентификации fí(в; Zk, ^Q^-1). В этом случае задача оценки неизвестного параметра требует решения задачи нелинейного программирования
§min = argmin fí (в; ZK, UQK-l), (10)
eeD(e)
где D(d) С R2 — область определения параметра в.
Одним из известных и наиболее популярных методов решения задачи параметрической идентификации является метод максимального правдоподобия [17], который заключается в отыскании экстремума функции правдоподобия. Для решения данной задачи можно использовать различные численные методы — градиентный метод, метод Ньютона, метаэвристические методы оптимизации (например генетический алгоритм или метод имитации отжига) и др. [18]. Готовые программные реализации данных методов, как правило, требуют от пользователя задания начального значения параметра в, ограничений на переменные и описания целевой функции. Вопрос выбора конкретного численного метода оптимизации также может быть предметом отдельного исследования.
В качестве целевой функции для реализации процедуры параметрической идентификации выберем критерий идентификации (10) в виде отрицательной логарифмической функции правдоподобия [17]
/KF_LR(e; ZK,U0K-l) =
= ln(2tf) + 2 E{\n[det(K,k(в))] + (0>k(0)}, (11)
k=l
где вектор невязки измерений щ (в) и его ковариационную матрицу
(0)= Е {vk (в)иТ(в)}
при заданных значениях параметра в вычисляют по известным уравнениям фильтра Калмана [10].
Алгоритм 1 (Стандартный фильтр Калмана для модели (9)).
® Начальные данные.
Положить Cq|q = со и Р0|о = По, где со — начальное значение оценки вектора состояния, По > 0 — начальное значение ковариационной матрицы ошибки оценивания.
® Рекуррентно обновлять величины (к > 0):
I. Экстраполяция.
1.1. Вычислить матрицу ковариации ошибки оценки
Pk+i = FPkik Fт.
1.2. Найти оценку вектора состояния Ck+i = Fck\k + Buk.
II. ФИЛЬТРАЦИЯ.
11.1. Вычислить матрицу Калмана
Кк+1 = Рк+1НтЪ-1к+1, где Х^+1 = К + НРк+1Нт.
11.2. Найти оценку вектора состояния
= ск+1 + Kk+lЬ'k+l, где ик+1 = хк+1 — Нск+1.
11.3. Вычислить матрицу ковариации ошибки оценки
Рк+11к+1 = Рк+1 — Кк+1НРк+1.
® КОНЕЦ.
Отметим, что решение задачи параметрической идентификации параметра в = (у, а)1 для дискретной модели (6) методом максимального правдоподобия на основе стандартного фильтра Калмана получено в [13]. Решение задачи численной идентификации скорости конвекции в модели конвективно-диффузионного переноса с помощью метаэвристических алгоритмов рассмотрено в [11] также с применением стандартного алгоритма Калмана.
Как уже было отмечено ранее, основной целью данной работы является разработка нового метода параметрической идентификации модели конвективно-диффузионного переноса, обладающего улучшенными вычислительными свойствами. Хорошо известно, что стандартная реализация фильтра (алгоритм 1) является неустойчивой по отношению к ошибкам машинного округления (см., например, подробный анализ и обсуждение в [10]). Для повышения точности оценки вектора состояния вместо стандартного фильтра Калмана предпочтительнее использовать его численно устойчивые модификации [19].
4. Идентификация параметров дискретной модели с применением БУЮ-фильтра Калмана. В данной работе мы предлагаем новый подход к идентификации параметров дискретных моделей конвективно-диффузионного переноса на основе численно устойчивого БУБ-фильтра Калмана.
Рассмотрим БУБ-факторизацию [20, теорема 1.1.6]. Любую матрицу А € Стхп ранга г можно представить в виде
А = Ж ХУ *, Х =
5 0 00
€ Стхп, 5 = &\щ{оъ...,ог},
где 'Ж € Стхт, У € Спхп — унитарные матрицы, У * —означает сопряженную и транспонированную к У, 5 € Мгхг — вещественная неотрицательная диагональная матрица. Величины И1 ^ 02 ^ ... ^ ог > 0 являются сингулярными значениями матрицы А. Заметим, что если г = п и/или г = т, некоторые из нулевых подматриц в Х отсутствуют.
Метод сингулярного разложения (или БУБ-факторизация) известен как наиболее точный метод факторизации матриц, особенно для матриц, близких к вырожденным. Кроме того, сингулярное разложение существует для любой матрицы, чего нельзя сказать, например, о разложении Холецкого. Поэтому модификации фильтра Калмана на основе БУБ-факторизации обладают такой же улучшенной численной устойчивостью к ошибкам машинного округления, что и все известные квадратно-корневые модификации [10,19]. Кроме того, БУБ-фильтры имеют дополнительные преимущества, такие как:
(1) все собственные значения матриц ковариации ошибок автоматически вычисляются на каждом шаге работы алгоритма фильтрации и могут использоваться для автоматического анализа и/или сокращения исходной модели;
(2) информационные матрицы (обратные к ковариационным) легко вычисляются путем инверсии диагональных факторов в БУБ-разложении, что создает элегантный способ построения алгоритмов информационного типа и фильтров смешанного типа с автоматическим переключением с ковариационного режима фильтрации на информационный.
Насколько известно авторам данной статьи, впервые идею построения численно устойчивой модификации фильтра Калмана с применением сингулярного разложения предложили Ошман и Бар-Ицхак [21]. Авторы назвали свой вариант БУБ-фильтра как У-Л-фильтр. Для реализации алгоритма необходимо выполнить как БУБ-разложение, так и разложение Холецкого, а также минимум три операции матричного обращения. Затем Ошманом была предложена информационная форма V-Л-фильтра [22]. В [23] У-Л-фильтр был применен для решения задачи параметрической идентификации линейной дискретной стохастической системы.
Позднее другие авторы предложили свои варианты БУБ-фильтра, соответствующие стандартному фильтру Калмана [24] и расширенному фильтру Калмана [25]. Указанные модификации во многом схожи с У-Л-фильтром. Ограничением их применения является требование положительной определенности матриц ковариации шумов в объекте и измерителе на каждой итерации работы алгоритма, поскольку в нем необходимо применять разложение Холецкого для вычисления квадратного корня ковариационной матрицы. Также необходимо выполнить минимум три матричных обращения.
С целью устранения указанных недостатков предыдущих версий БУБ-фильтра в [26] авторами был предложен новый, улучшенный вариант БУБ-фильтра Калмана. Его отличие от других известных нам вариантов в том, что данная модификация фильтра Калмана свободна от выполнения условий Qk ^ 0 и К > 0, требуемых как в стандартном фильтре Калмана, так и во всех его квадратно-корневых модификациях [19]. Другим значимым преимуществом указанной модификации является наличие всего лишь одного обращения диагональной матрицы в уравнениях фильтра. Согласно [26], по результатам сравнительного анализа данный вариант БУБ-фильтра показал наилучшие результаты в плане численной устойчивости к ошибкам машинного округления на примере решения плохообусловленных задач.
Следует также отметить, что БУБ-фильтр подтвердил свою эффективность при решении задач параметрической идентификации [27], задачи оценки состояния и параметров полета летательного аппарата [25], задачи калма-новской фильтрации показаний инерциального измерительного модуля (1Ми) [28] и др.
Учитывая изложенное выше, в данной работе для построения процедуры параметрической идентификации мы выбрали улучшенный вариант БУБ-фильтра [26].
Рассмотрим представление ковариационных матриц ошибок оценивания в виде Рк = ОркОрк ©^, где ©рк — ортогональная матрица и Ирк — диагональная матрица, содержащая сингулярные значения матрицы Рк. Уравнения БУБ-фильтра позволяют рекуррентно обновлять БУБ-факторы {© Рк ,0Рк}
матрицы Рк с помощью сингулярного разложения (процедуры БУБ-факто-ризации). Далее запишем уравнения БУБ-фильтра с учетом модели (9).
Алгоритм 2 (ЯУБ-фильтр для модели (9)).
® Начальные данные.
Выполнить SVD-разложение матриц По = вп0^п0@п0, R = QrDrQr. Положить Со|о = со, ©р0|0 = впо, Dр/0
1/2
D
1/2 п0.
® Рекуррентно овновлять величины (к > 0): I. Экстраполяция.
1.1. Выполнить БУБ-разложение матрицы
D%2 вП Fk
pk|k Рк|к к
т
= Wtu
D
1/2
Рк+l
в
п
Рк+l ■
1.2. Получить БУБ-факторы |врй+1, Орк+1} матрицы Рк+1.
1.3. Найти оценку вектора состояния Ск+\ = Рсщк + Ви,к. II. Фильтрация.
11.1. Выполнить БУБ-разложение матрицы
dr/2 вп Dpk+ieik+i^ П
=
(1)
ми
D
1/2
!> ,к + 1
0
вп
v,k + 1
11.2. Найти
Кк+1 = Рк+\Н ТвЕ1/,к+1, Кк+\ =
11.3. Выполнить БУБ-разложение матрицы
-1
v,k+1
в
Т
',к + 1
D
Р/2+1 вПк+1 V -K^+iH)
Щ2в1ктк+1
т
=
(2)
ми
D
1/2
Рк + 1|к + 1 0
в Рк
к+1| к+1
II.4. Получить SVD-факторы {врк+1|к+1 ,Dpk+1|k+1} матрицы Рк+1|к+1.
11.5. Найти оценку вектора состояния Ск+1|к+1 = Ск+1 + Кк+\0-1 ^ ^к+1, гДе "к+1 = к+1 (Хк+1 - Яск+1). ® Конец.
Теперь для построения процедуры параметрической идентификации необходимо переписать выражение для вычисления логарифмической функции правдоподобия (11) в терминах БУБ-фильтра. Учитывая, что
^Е^к) = ) и "к К,к "к = "к,
запишем
_lr(0; zi , Uq 1) =
К
K^ ln(2vr) + 2E {ln[det(Dsl/,k)] + D-^ 4 }, (12)
к=1
К
где диагональная матрица к и вектор ^ доступны на каждом шаге работы алгоритма 2.
5. Численный анализ эффективности предложенного подхода.
Рассмотрим на численных примерах работоспособность и преимущества предложенного подхода. Пусть требуется идентифицировать параметры V и а следующей модели:
дс дс д2с дЪ дх д 2х'
с(х, 0) = 0, с(0,£) = 4*| 8ш
с(м) = г,
(13)
(14)
(15)
(16)
или
дс(1,г)
дх
= -[с(1,1) -*],
(17)
где с(х, ¿) — концентрация вещества в одномерном потоке; х € [0; 1]; £ € [0; 2]; начальное условие (14) соответствует концентрации вещества в начальный момент времени; граничное условие (15)—периодическому изменению концентрации вещества с возрастающей амплитудой на левом конце отрезка; граничное условие (16) —линейному возрастанию концентрации вещества на правом конце отрезка; граничное условие (17) —линейному возрастанию концентрации вещества в окружающей среде на правом конце отрезка.
Пусть в уравнении (13) V = 2, а = 1. Процесс идентификации будем моделировать в системе МЛТЬЛБ при помощи специализированного программного комплекса авторской разработки. Зададим на отрезке [0; 1] пространственную сетку с 6 узлами по оси Ох (Ах = 0.2). Шаг по оси ОЬ выберем
(Ах)2
Аt €
= 0.02, что
из условия устойчивости конечно-разностной схемы. 2
соответствует К = 101. Получим решения рассматриваемой задачи методом конечных разностей для различных комбинаций граничных условий. На рис. 1 приведены графики соответствующих решений.
Рассмотрим процесс идентификации коэффициентов уравнения (13) по результатам зашумленных измерений в узлах пространственной сетки, соответствующих первой и последней компонентам вектора состояния. В этом случае для задачи с граничными условиями (15), (16) измеритель будет иметь вид
+
к = 1, 2,
,К,
а для задачи с граничными условиями (15), (17):
1 0 0 0 0 00001
+
к = 1, 2,
е
г
с
2
с
с
с
с
г
2
О О
b
Рис. 1. Графики решения задачи с граничными условиями (15), (16) (a)
и (15), (17) (b)
[Figure 1. Plots of solutions of problem with boundary conditions (15), (16)
(a) and (15), (17) (b)]
a
b
Рис. 2. Графики измерений для задачи с граничными условиями (15), (16) (a)
и (15), (17) (b)
[Figure 2. Plots of measurements for problem with boundary conditions (15), (16) (a)
and (15), (17) (b)]
Таблица 1
Средние значения параметров и ошибки идентификации для R — Rq и граничных условий (15), (16) [Mean values and identification errors for R — RQ and boundary
conditions (15), (16)]
Algorithm Mean values RMSE MAPE
V a V a V a
CKF SVD 1.9951 1.9951 0.9962 0.9962 4.63 • 10-2 4.63 • 10-2 2.23 • 10-2 2.23 • 10-2 1.8242 1.8242 1.8325 1.8325
Таблица 2
Средние значения параметров и ошибки идентификации для R — Rq и граничных условий (15), (17) [Mean values and identification errors for R — Rq and boundary
conditions (15), (17)]
Algorithm Mean values RMSE MAPE
V a V a V a
CKF SVD 1.9996 1.9996 0.9984 0.9984 4.99 • 10-2 4.99 • 10-2 2.31 • 10-2 2.31 • 10-2 1.9539 1.9539 1.8660 1.8660
Таблица 3
Средние значения параметров и ошибки идентификации для R — Ri и граничных условий (15), (16) [Mean values and identification errors for R — R1 and boundary
conditions (15), (16)]
5 Algorithm Mean values RMSE MAPE
V a V a V a
10-1Q CKF SVD 2.0000 2.0000 1.0000 1.0000 4.55 • 10-6 4.55 • 10-6 2.34 • 10-6 2.34 • 10-6 1.80 • 10-4 1.80 • 10-4 1.80 • 10-4 1.85 • 10-4
10-11 CKF SVD 2.0000 2.0000 1.0000 1.0000 1.47 • 10-6 1.47 • 10-6 7.09 • 10-7 7.10 • 10-7 5.91 • 10-5 5.91 • 10-5 5.78 • 10-5 5.78 • 10-5
10-12 CKF SVD 2.0000 2.0000 1.0000 1.0000 5.09 • 10-7 5.09 • 10-7 1.94 • 10-7 1.94 • 10-7 1.99 • 10-5 1.99 • 10-5 1.56 • 10-5 1.55 • 10-5
10-13 CKF SVD 2.0000 2.0000 1.0000 1.0000 1.65 • 10-7 1.65 • 10-7 6.73 • 10-8 6.74 • 10-8 6.60 • 10-6 6.60 • 10-6 5.40 • 10-6 5.40 • 10-6
10-14 CKF SVD 2.0909 2.0000 1.1295 1.0000 3.72 • 10-1 4.95 • 10-8 2.45 • 10-1 2.44 • 10-8 13.3856 2.00 • 10-6 18.2959 1.97 • 10-6
10-15 CKF SVD 2.4525 2.0000 2.4007 1.0000 4.92 • 10-1 2.84 • 10-8 1.4422 1.38 • 10-8 24.2181 1.24 • 10-6 141.0602 1.21 • 10-6
10-16 CKF SVD 2.0000 1.0000 2.77 • 10-8 1.29 • 10-8 1.36 • 10-6 1.26 • 10-6
Таблица 4
Средние значения параметров и ошибки идентификации для R = R2 и граничных условий (15), (16) [Mean values and identification errors for R = R2 and boundary
conditions (15), (16)]
5 Algorithm Mean values RMSE MAPE
V a V a V a
10-10 CKF SVD 1.9998 1.9998 1.0002 1.0002 8.53 • 10-3 8.53 • 10-3 3.97 • 10-3 3.97 • 10-3 3.45 • 10-1 3.45 • 10-1 3.09 • 10-1 3.09 • 10-1
10-11 CKF SVD 1.9996 1.9996 1.0003 1.0003 8. 8. 77 to to i—1 i—1 о о со со 3.81 • 10-3 3.81 • 10-3 3.44 • 10-1 3.44 • 10-1 2.98 • 10-1 2.98 • 10-1
10-12 CKF SVD 1.9998 1.9999 1.0003 1.0003 9.17 • 10-3 9.17 • 10-3 3.58 • 10-3 3.58 • 10-3 3.70 • 10-1 3.70 • 10-1 2.91 • 10-1 2.91 • 10-1
10-13 CKF SVD 2.0008 2.0008 0.9999 0.9999 9.03 • 10-3 9.04 • 10-3 3.80 • 10-3 3.80 • 10-3 3.58 • 10-1 3.58 • 10-1 3.06 • 10-1 3.07 • 10-1
10-14 CKF SVD 2.0011 1.9998 1.1824 1.0000 3.80 • 10-1 8.23 • 10-3 3.53 • 10-1 3.70 • 10-3 8.5127 3.36 • 10-1 18.6786 2.96 • 10-1
10-15 CKF SVD 2.0556 1.9998 1.2045 1.0000 6.08 • 10-1 8.67 • 10-3 3.75 • 10-1 3.25 • 10-3 21.0730 3.30 • 10-1 23.0515 2.64 • 10-1
10-16 CKF SVD 2.1827 2.0002 1.2980 1.0001 7.79 • 10-1 9.15 • 10-3 6.69 • 10-1 3.57 • 10-3 29.2086 3.59 • 10-1 41.3944 2.79 • 10-1
Таблица 5
Средние значения параметров и ошибки идентификации для R = Дз и граничных условий (15), (16) [Mean values and identification errors for R = R3 and boundary
conditions (15), (16)]
5 Algorithm Mean values RMSE MAPE
V a V a V a
10-10 CKF SVD 2.0000 2.0000 1.0000 1.0000 5.88 • 10-6 4.79 • 10-6 2.97 • 10-6 2.97 • 10-6 2.40 • 10-4 2.40 • 10-4 2.40 • 10-4 2.40 • 10-4
10-11 CKF SVD 2.0000 2.0000 1.0000 1.0000 1.67 • 10-6 1.67 • 10-6 9.22 • 10-7 9.21 • 10-7 6.49 • 10-5 6.48 • 10-5 7.50 • 10-5 7.49 • 10-5
10-12 CKF SVD 2.0000 2.0000 1.0000 1.0000 5.63 • 10-7 5.63 • 10-7 2.94 • 10-7 2.94 • 10-7 2.25 • 10-5 2.25 • 10-5 2.41 • 10-5 2.42 • 10-5
10-13 CKF SVD 2.0000 2.0000 1.0000 1.0000 2.04 • 10-7 2.05 • 10-7 9.85 • 10-8 9.85 • 10-8 8.06 • 10-6 8.08 • 10-6 7.78 • 10-6 7.78 • 10-6
10-14 CKF SVD 2.4218 2.0000 2.3225 1.0000 5.04 • 10-1 6.59 • 10-8 1.3914 3.38 • 10-8 23.2788 2.60 • 10-6 132.5812 2.56 • 10-6
10-15 CKF SVD 2.0530 2.0000 1.2871 1.0000 6.02 • 10-1 4.32 • 10-8 3.99 • 10-1 2.13 • 10-8 24.6737 1.94 • 10-6 30.9475 1.93 • 10-6
10-16 CKF SVD 3.0785 2.0000 1.6264 1.0000 1.1877 3.54 • 10-8 0.6924 1.79 • 10-8 55.5319 1.74 • 10-6 63.5329 1.76 • 10-6
Таблица 6
Средние значения параметров и ошибки идентификации для R — Ri и граничных условий (15), (17) [Mean values and identification errors for R — Ri and boundary
conditions (15), (17)]
5 Algorithm Mean values RMSE MAPE
V a V a V a
10-10 CKF SVD 2.0000 2.0000 1.0000 1.0000 5.06 • 10-6 5.06 • 10-6 2.23 • 10-6 2.23 • 10-6 2.03 • 10-4 2.03 • 10-4 1.81 • 10-4 1.81 • 10-4
10-11 CKF SVD 2.0000 2.0000 1.0000 1.0000 1.47 • 10-6 1.47 • 10-6 7.24 • 10-7 7.25 • 10-7 5.80 • 10-5 5.80 • 10-5 5.89 • 10-5 5.90 • 10-5
10-12 CKF SVD 2.2732 2.0000 1.8860 1.0000 5.42 • 10-1 4.81 • 10-7 1.1367 2.51 • 10-7 19.9823 1.80 • 10-5 88.8022 1.97 • 10-5
10-13 CKF SVD 2.0295 2.0000 1.5567 1.0000 4.69 • 10-1 1.54 • 10-7 8.75 • 10-1 7.30 • 10-8 15.5471 6.09 • 10-6 56.4061 5.92 • 10-6
10-14 CKF SVD 1.7208 2.0000 1.8216 1.0000 9.00 • 10-1 5.17 • 10-8 1.06342 2.34 • 10-8 36.9710 2.06 • 10-6 85.6284 1.85 • 10-6
10-15 CKF SVD 2.0000 1.0000 1.86 • 10-8 8.75 • 10-9 7.55 • 10-7 7.05 • 10-7
10-16 CKF SVD 2.0000 1.0000 1.02 • 10-8 6.07 • 10-9 4.54 • 10-7 5.55 • 10-7
Таблица 7
Средние значения параметров и ошибки идентификации для R — R2 и граничных условий (15), (17) [Mean values and identification errors for R — R2 and boundary
conditions (15), (17)]
5 Algorithm Mean values RMSE MAPE
V a V a V a
10-10 CKF SVD 1.9998 1.9998 1.0000 1.0000 1.12 • 10-2 1.12 • 10-2 6.20 • 10-3 6.20 • 10-3 4.48 • 10-1 4.48 • 10-1 5.00 • 10-1 5.00 • 10-1
10-11 CKF SVD 2.0011 2.0011 1.0004 1.0004 1.20 • 10-2 1.20 • 10-2 6.30 • 10-3 6.30 • 10-3 4.82 • 10-1 4.82 • 10-1 5.02 • 10-1 5.02 • 10-1
10-12 CKF SVD 2.0867 2.0008 1.1835 1.0000 4.46 • 10-1 1.18 • 10-2 4.78 • 10-1 5.00 • 10-3 8.9741 4.63 • 10-1 19.0136 4.60 • 10-1
10-13 CKF SVD 1.7084 2.0000 1.4499 1.0000 8.48 • 10-1 1.14 • 10-2 6.67 • 10-1 6.37 • 10-3 28.9005 4.53 • 10-1 45.2698 5.07 • 10-1
10-14 CKF SVD 1.3955 1.9994 1.3397 1.0004 8.62 • 10-1 1.19 • 10-2 5.20 • 10-1 6.03 • 10-3 32.7939 4.54 • 10-1 34.9237 4.66 • 10-1
10-15 CKF SVD 1.5160 1.9988 1.2420 1.0006 8.89 • 10-1 1.12 • 10-2 4.28 • 10-1 6.11 • 10-3 37.9516 4.40 • 10-1 30.0707 4.85 • 10-1
10-16 CKF SVD 2.0004 1.0001 1.04 • 10-2 6.32 • 10-3 4.18 • 10-1 5.09 • 10-1
Таблица 8
Средние значения параметров и ошибки идентификации для Д = Д3 и граничных условий (15), (17) [Mean values and identification errors for R = Д3 and boundary
conditions (15), (17)]
5 Algorithm Mean values RMSE MAPE
V a V a V а
10-10 CKF SVD 2.0000 2.0000 1.0000 1.0000 5.90 • 10-6 5.90 • 10-6 2. 2. 88 to to 11 о о СП СП 2.38 • 10-4 2.38 • 10-4 2.24 • 10-4 2.24 • 10-4
10-11 CKF SVD 2.0000 2.0000 1.0000 1.0000 1.76 • 10-6 1.76 • 10-6 8.62 • 10-7 8.62 • 10-7 7.07 • 10-5 7.07 • 10-5 6.86 • 10-5 6.86 • 10-5
10-12 CKF SVD 2.2020 2.0000 1.7162 1.0000 4.04 • 10-1 6.42 • 10-7 1.0238 2.59 • 10-7 13.9583 2.65 • 10-5 71.8292 2.05 • 10-5
10-13 CKF SVD 2.4736 2.0000 2.4251 1.0000 4.86 • 10-1 1.64 • 10-7 1.4618 9.01 • 10-8 23.6781 6.53 • 10-6 142.5113 7.16 • 10-6
10-14 CKF SVD 1.9462 2.0000 2.2192 1.0000 7.26 • 10-1 5.46 • 10-8 1.3157 2.38 • 10-8 29.1861 2.21 • 10-6 121.9198 1.91 • 10-6
10-15 CKF SVD 2.3016 2.0000 2.2533 1.0000 6.05 • 10-1 1.88 • 10-8 1.3677 9.76 • 10-9 27.1425 7.56 • 10-7 128.7952 7.76 • 10-7
10-16 CKF SVD 2.0000 1.0000 1.12 • 10-8 5.98 • 10-9 4.86 • 10-7 5.43 • 10-7
Решение задачи получим при следующих вариантах матрицы ковариации шума Д в измерителе:
Rq =
10-2 0
0
10
-2
Ri =
S 0
0 s
R2 =
'S 0 0 10
-4
R3 =
10
4
0
0 s
где 5 = 10-10,10-11,..., 10-16. При уменьшении 5 погрешность измерений соответствующего сенсора уменьшается, что соответствует случаю «почти точных» измерений, которые могут приводить к вырожденности и плохой обусловленности ковариационной матрицы ошибок оценивания в стандартном фильтре Калмана.
На рис. 2 приведены примеры графиков смоделированных измерений с матрицей ковариации шума Д = До для рассматриваемой задачи с различными комбинациями граничных условий.
В качестве критериев идентификации возьмем логарифмические функции правдоподобия (11) и (12).
В табл. 1-8 приведены результаты численных экспериментов по идентификации коэффициентов V и а для различных значений матрицы Д и граничных условий. Минимизация критериев идентификации выполнялась в области = [в = (-и, а)т | V € [0; 5], а € [0; 5]} при помощи функции £ш1псоп системы ИАТЬАБ. В качестве начальной точки для функции £ш1псоп выбирался центр области Для каждого варианта матрицы Д и каждого значения 5 выполнялось усреднение найденных параметров и вычисление ошибок ЯМБЕ и МАРЕ по результатам 200 запусков.
Полученные данные численных экспериментов показывают, что сначала для обоих алгоритмов с увеличением точности измерений средние значения
идентифицированных параметров v и а стремятся к истинным (2 и 1 соответственно), а ошибки RMSE и MAPE уменьшаются. Однако начиная со значения 5 = 10-14 для граничных условий (15), (16) и 5 = 10-13 для граничных условий (15), (17) ошибки процедуры идентификации на основе стандартного фильтра Калмана начинают быстро нарастать вплоть до аварийного завершения работы функции минимизации fmincon (прочерки в таблицах). Данный факт обусловлен расходимостью стандартного фильтра Калмана (резким ухудшением обусловленности ковариационной матрицы при увеличении точности измерений) и, как следствие, некорректным вычислением значений критерия идентификации (11). В то же время процедура идентификации на основе SVD-фильтра выполняется корректно для всех значений
Таким образом, при использовании рекуррентных методов параметрической идентификации непосредственное влияние на качество идентификации могут оказывать вычислительные свойства алгоритма оптимальной дискретной фильтрации, на основе которого вычисляются значения критерия идентификации.
Заметим, что все вычисления в системе MATLAB по умолчанию проводятся с двойной точностью, однако на практике встречаются ситуации, когда вычисления необходимо проводить с одинарной точностью, например, во встраиваемых системах [28]. В этом случае использование SVD-фильтра, несмотря на больший объем вычислений, будет более предпочтительным.
Заключение. Рассмотрена задача математического моделирования процесса параметрической идентификации моделей конвективно-диффузионного переноса, описываемых уравнениями в частных производных с начальным и граничными условиями, с применением SVD-модификации фильтра Кал-мана. Рассматриваются одномерные модели с постоянными коэффициентами, граничными условиями первого рода или смешанными граничными условиями первого и третьего рода. Начальное и граничные условия предполагаются известными.
Метод решения задачи основан на переходе от исходной непрерывной модели с уравнением в частных производных к модели, описываемой линейной дискретной динамической системой в пространстве состояний, с известными входными воздействиями, двумя неизвестными параметрами, соответствующими скорости конвекции v и коэффициенту диффузии а, и с зашумленными измерениями, моделирующими экспериментальные данные.
В работе показано, что для выполнения условия полной наблюдаемости моделей достаточно, чтобы матрица измерений Н состояла всего из двух строк, соответствующих наличию в модели измерителя двух сенсоров.
Основным результатом работы является новый рекуррентный метод параметрической идентификации моделей конвективно-диффузионного переноса, основанный на использовании метода максимального правдоподобия с построением критерия идентификации (функции правдоподобия) на основе величин, вычисляемых SVD-модификацией фильтра Калмана. Оптимизация критерия идентификации выполнена с помощью функции fmincon системы MATLAB.
С помощью компьютерного моделирования в системе MATLAB проведен сравнительный анализ численной эффективности алгоритма параметрической идентификации с применением стандартного фильтра Калмана и SVD-фильтра. Результаты численных экспериментов подтверждают работоспособность предложенного метода и его надежность в плане численной устойчиво-
сти к ошибкам машинного округления при решении задачи параметрической идентификации.
Конкурирующие интересы. Конкурирующих интересов не имеем. Авторский вклад и ответственность. Все авторы принимали участие в разработке концепции статьи и в написании рукописи. Авторы несут полную ответственность за предоставление окончательной рукописи в печать. Окончательная версия рукописи одобрена всеми авторами.
Финансирование. Исследование выполнено при финансовой поддержке РФФИ и Ульяновской области в рамках научного проекта № 19-41-730009.
Библиографический список
1. Леонтьев А. И., Кожинов И. А., Исаев С. И. и др. Теория тепломассообмена. М.: МГТУ им. Н. Э. Баумана, 2018. 462 с.
2. Farlow S. J. Partial Differential Equations for Scientists and Engineers. New York: Dover Publ., 1982. ix+414 pp.
3. Денисов А. М. Введение в теорию обратных задач. М.: МГУ, 1994. 208 с.
4. Мацевитый Ю. М., Мултановский А. В. Идентификация параметров теплообмена методом оптимальной динамической фильтрации // Теплофизика высоких температур, 1979. Т. 17, №5. С. 1053-1060.
5. Карпов А. А., Тихонова Т. А. Восстановление нестационарных тепловых потоков по экспериментальным данным // Матем. модел, 2000. Т. 12, №5. С. 101-106.
6. Симбирский Г. Д., Лантрат В. К. Применение цифрового фильтра Калмана для параметрической идентификации высокотемпературного термопреобразователя // Автомобиль и электроника. Современные технологии, 2017. №11. С. 68-75.
7. Пилипенко Н. В. Применение фильтра Калмана в нестационарной теплометрии. СПб.: Унив. ИТМО, 2017. 36 с.
8. Матвеев М. Г., Копытин А. В., Сирота Е. А. Комбинированный метод идентификации параметров распределенной динамической модели / Информационные технологии и нанотехнологии (ИТНТ-2018): сборник трудов 4-й Международной конференции и молодежной школы. Самара, 2018. С. 1651-1657.
9. Пилипенко Н. В., Заричняк Ю. П., Иванов В. А., Халявин А. М. Параметрическая идентификация дифференциально-разностных моделей теплопереноса в одномерных телах на основе алгоритмов фильтра Калмана // Научно-технический вестник информационных технологий, механики и оптики, 2020. Т. 20, № 4. С. 584-588. https:// doi.org/10.17586/2226-1494-2020-20-4-584-588.
10. Grewal M. S., Andrews A. P. Kalman Filtering. Theory and Practice with MATLAB. Hoboken, NJ: John Wiley and Sons, 2015. xvii+617 pp. https://doi.org/10.1002/ 9781118984987.
11. Tsyganov A. V., Tsyganova Yu. V., Kuvshinova A. N., Tapia Garza H. R. Metaheuristic algorithms for identification of the convection velocity in the convection-diffusion transport model/ CEUR Workshop Proceedings. vol.2258, 2018. pp. 188-196. http://ceur-ws.org/ Vol-2258/paper24.pdf.
12. Кувшинова А. Н. Динамическая идентификация смешанных граничных условий в модели конвективно-диффузионного переноса в условиях зашумленных измерений // Журнал СВМО, 2019. Т. 21, №4. С. 469-479. https://doi.org/10.15507/2079-6900.21. 201904.469-479.
13. Kuvshinova A. N., Tsyganov A. V., Tsyganova Yu. V., Tapia Garza H. R. Parameter identification algorithm for convection-diffusion transport model// J. Phys.: Conf. Ser., 2021. vol. 1745, 012110. https://doi.org/10.1088/1742-6596/1745/1/012110.
14. Фомин В. Н. Рекуррентное оценивание и адаптивная фильтрация. М.: Наука, 1984. 288 с.
15. Maybeck P. S. Stochastic Models, Estimation, and Control. vol. 1 / Mathematics in Science and Engineering. vol. 141. New York, San Francisco, London: Academic Press, 1979. xix+423 pp.
16. Кувшинова А. Н. Анализ дискретной линейной стохастической модели конвективно-диффузионного переноса // Ученые записки УлГУ. Сер. Математика и информационные технологии, 2019. №1. С. 65-69.
17. Âstrôm K. J. Maximum likelihood and prediction error methods // Automatica, 1980. vol. 16, no. 5. pp. 551-574. https://doi.org/10.1016/0005-1098(80)90078-3.
18. Васильев В. П. Численные методы решения экстремальных задач. М.: Мир, 1982. 520 с.
19. Цыганова Ю. В., Куликова М. В. О современных ортогонализованных алгоритмах оптимальной дискретной фильтрации // Вестн. ЮУрГУ. Сер. Матем. моделирование и программирование, 2018. Т. 11, №4. С. 5-30. https://doi.org/10.14529/mmp180401.
20. Bjôrck Â. Numerical Methods in Matrix Computations / Texts in Applied Mathematics. vol.59. Cham: Springer, 2015. xvi+800 pp. https://doi.org/10.1007/ 978-3-319-05089-8.
21. Oshman Y., Bar-Itzhack I. Y. Square root filtering via covariance and information eigenfactors// Automatica, 1986. vol.22, no. 5. pp. 599-604. https://doi.org/10.1016/ 0005- 1098(86)90070- 1.
22. Oshman Y. Square root information filtering using the covariance spectral decomposition/ Proc. of the 27th Conf. on Decision and Control, 1988. pp. 382-387. https://doi.org/ 10.1109/CDC.1988.194335.
23. Oshman Y. Maximum likelihood state and parameter estimation via derivatives of the V-Lambda filter // J. Guid. Control Dyn., 1992. vol. 15, no. 3. pp. 717-726. https://doi.org/ 10.2514/3.20896.
24. Wang L., Libert G., Manneback P. Kalman filter algorithm based on singular value decomposition/ Proc. of the 31st Conf. on Decision and Control, 1992. pp. 1224-1229. https:// doi.org/10.1109/IECQN.1992.254406.
25. Zhang Y., Dai G., Zhang H., Li Q. A SVD-based extended Kalman filter and applications to aircraft flight state and parameter estimation/ Proc. of 1994 American Control Conf., 1994. pp. 1809-1813. https://doi.org/10.1109/ACC.1994.752384.
26. Kulikova M. V., Tsyganova J. V. Improved discrete-time Kalman filtering within singular value decomposition // IET Control Theory Appl., 2017. vol. 11, no. 15. pp. 2412-2418, arXiv: 1611.03686 [math.OC]. https://doi.org/10.1049/iet- cta.2016.1282.
27. Tsyganova J. V., Kulikova M. V. SVD-based Kalman filter derivative computation// IEEE Trans. Autom. Control, 2017. vol. 62, no. 9. pp. 4869-4875, arXiv: 1612.04777 [cs.SY]. https://doi.org/10.1109/TAC.2017.2694350.
28. Alessandrini M., Biagetti G., Crippa P., Falaschetti L., Manoni L., Turchetti C. Singular value decomposition in embedded systems based on ARM Cortex-M architecture // Electronics, 2021. vol.10, no. 1, 34. https://doi.org/10.3390/electronics10010034.
Vestn. Samar. Gos. Tekhn. Univ., Ser. Fiz.-Mat. Nauki
[J. Samara State Tech. Univ., Ser. Phys. Math. Sci.], 2021, vol. 25, no. 4, pp. 716-737 ISSN: 2310-7081 (online), 1991-8615 (print) d https://doi.org/10.14498/vsgtu1876
MSC: 93A30, 65C20
Mathematical modeling of parameter identification process of convection-diffusion transport models using the SVD-based Kalman filter
A. N. Kuvshinova1, A. V. Tsyganov1, Yu. V. Tsyganova2
1 Ilya Ulyanov State Pedagogical University,
4/5, Lenin Square, Ulyanovsk, 432071, Russian Federation.
2 Ulyanovsk State University,
42, L. Tolstoy st., Ulyanovsk, 432017, Russian Federation.
Abstract
The paper addresses a problem of mathematical modeling of the process of identifying the coefficients of a partial differential equation in convection-diffusion transport models based on the results of noisy measurements of the function values. Identification process is performed using a new method belonging to the class of recurrent parameter identification methods based on optimal discrete Kalman-type filtering algorithms. One-dimensional models with constant coefficients, boundary conditions of first kind, or mixed boundary conditions of first and third kind are considered.
The proposed method is based on the transition from the initial continuous model with a partial differential equation to the model described by the state-space linear discrete-time dynamic system and the application of the maximum likelihood method to it with construction of an identification criterion (likelihood function) based on the values calculated by the SVD algorithm of the Kalman filtering. This filter is based on the singular value decomposition of error covariance matrix and works stably even in cases when it is close to singular. The SVD filter has proven itself well in solving various problems of discrete filtering and parameter identification. It has several advantages over the traditionally used conventional Kalman filter. The main of which is robustness against machine roundoff errors.
Research Article
© Authors, 2021
© Samara State Technical University, 2021 (Compilation, Design, and Layout)
3 ©® The content is published under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/) Please cite this paper in press as:
Kuvshinova A. N., Tsyganov A. V., Tsyganova Yu. V. Mathematical modeling of parameter identification process of convection-diffusion transport models using the SVD-based Kalman filter, Vestn. Samar. Gos. Tekhn. Univ., Ser. Fiz.-Mat. Nauki [J. Samara State Tech. Univ., Ser. Phys. Math. Sci.], 2021, vol. 25, no. 4, pp. 716-737. https://doi.org/10.14498/vsgtu1876 (In Russian). Authors' Details:
Anastasia N. Kuvshinova © https://orcid.org/0000-0002-3496-5981 Postgraduate Student; Dept. of Higher Mathematics; e-mail:kuvanulspu@yandex.ru Andrey V. Tsyganov © https://orcid.org/0000-0002-4173-5199 Cand. Phys. & Math. Sci; Professor; Dept. of Higher Mathematics; e-mail: andrew.tsyganov@gmail.com
Yulia V. Tsyganova A https://orcid.org/0000-0001-8812-6035 Dr. Phys. & Math. Sci.; Professor; Dept. of Information Technology; e-mail: tsyganovajv@gmail .com
Computer modeling of parameter identification has been processed with the MATLAB system using a specialized software package. The results of numerical experiments confirm the efficiency of the proposed method and its advantages compared to the similar one based on the conventional Kalman filter.
Keywords: convection-diffusion transport model, parameter identification, Kalman filter, SVD filter.
Received: 3rd August, 2021 / Revised: 7th December, 2021 / Accepted: 21st December, 2021 / First online: 28th December, 2021
Competing interests. We have no competing interests.
Authors' contributions and responsibilities. Each author has participated in the article concept development and in the manuscript writing. The authors are absolutely responsible for submit the final manuscript to print. Each author has approved the final version of manuscript.
Funding. The reported study was funded by RFBR and Ulyanovsk region, project no. 19-41-730009.
References
1. Leont'ev A. I., Kozhinov I. A., Isaev S. I., et al. Teoriia teplomassoobmena [Theory of Heat and Mass Transfer]. Moscow, Bauman Moscow State Techn. Univ., 2018, 462 pp. (In Russian)
2. Farlow S. J. Partial Differential Equations for Scientists and Engineers. New York, Dover Publ., 1982, ix+414 pp.
3. Denisov A. M. Vvedenie v teoriiu obratnykh zadach [Introduction to the Theory of Inverse Problems]. Moscow, Moscow State Univ., 1994, 208 pp. (In Russian)
4. Matsevityi Yu. M., Multanovskii A. V. Identification of heat transfer parameters by the method of optimal dynamic filtering, High Temperature, 1979, vol. 17, no. 5, pp. 1053-1060 (In Russian).
5. Karpov A. A., Tikhonova T. A. Recovery of non-stationary heat flows from experimental data, Matem. Mod., 2000, vol. 12, no. 5, pp. 101-106 (In Russian).
6. Simbirskiy G. D., Lantrat V. K. Application of the Kalman digital filter for parametric identification high-temperature thermocouple, Autom. Electron. Modern Technology, 2017, no. 11, pp. 68-75 (In Russian).
7. Pilipenko N. V. Primenenie fil'tra Kalmana v nestatsionarnoi teplometrii [Applying the Kalman Filter in Non-Stationary Heat Metering]. St. Petersburg, ITMO Univ., 2017, 36 pp. (In Russian)
8. Matveev M. G., Kopytin A. V., Sirota E. A. Combined method for identifying the parameters of a distributed dynamic model, In: Proc. IV Int. Conf. (ITNT, 2018). Samara, 2018, pp. 1651-1657 (In Russian).
9. Pilipenko N. V., Zarichnyak Yu. P., Ivanov V. A., Khalyavin A. M. Parametric identification of differencial-difference models of heat transfer in one-dimensional bodies based on Kalman filter algorithms, Scientific and Technical Journal of Information Technologies, Mechanics and Optics, 2020, vol.20, no. 4, pp. 584-588 (In Russian). https://doi.org/10.17586/ 2226-1494-2020-20-4-584-588.
10. Grewal M. S., Andrews A. P. Kalman Filtering. Theory and Practice with MATLAB. Hoboken, NJ, John Wiley and Sons, 2015, xvii+617 pp. https://doi.org/10.1002/ 9781118984987.
11. Tsyganov A. V., Tsyganova Yu. V., Kuvshinova A. N., Tapia Garza H. R. Metaheuristic algorithms for identification of the convection velocity in the convection-diffusion transport
model, In: CEUR Workshop Proceedings, vol.2258, 2018, pp. 188-196. http://ceur-ws. org/Vol-2258/paper24.pdf.
12. Kuvshinova A. N. Dynamic identification of boundary conditions for convection-diffusion transport model in the case of noisy measurements, Zhurnal SVMO, 2019, vol.21, no. 4, pp. 469-479 (In Russian). https://doi.org/10.15507/2079-6900.21.201904.469-479.
13. Kuvshinova A. N., Tsyganov A. V., Tsyganova Yu. V., Tapia Garza H. R. Parameter identification algorithm for convection-diffusion transport model, J. Phys.: Conf. Ser., 2021, vol. 1745, 012110. https://doi.org/10.1088/1742-6596/1745/1/012110.
14. Fomin V. N. Rekurrentnoe otsenivanie i adaptivnaia fil'tratsiia [Recurrent Estimation and Adaptive Filtering]. Moscow, Nauka, 1984, 288 pp. (In Russian)
15. Maybeck P. S. Stochastic Models, Estimation, and Control, vol.1, Mathematics in Science and Engineering, vol. 141. New York, San Francisco, London, Academic Press, 1979, xix+423 pp.
16. Kuvshinova A. N. Analysis of discrete linear stochastic model of convection-diffusion transport, Uchenye Zapiski UlGU. Ser. Matem. Inform. Tekhn., 2019, no. 1, pp. 65-69 (In Russian).
17. Âstrôm K. J. Maximum likelihood and prediction error methods, Automatica, 1980, vol. 16, no. 5, pp. 551-574. https://doi.org/10.1016/0005-1098(80)90078-3.
18. Vasil'ev V. P. Chislennye metody resheniia ekstremal'nykh zadach [Numerical Methods for Solving Extreme Problems]. Moscow, Mir, 1982, 520 pp. (In Russian)
19. Tsyganova Yu. V., Kulikova M. V. On modern array algorithms for optimal discrete filtering, Vestn. YuUrGU. Ser. Mat. Model. Progr., 2018, vol.11, no. 4, pp. 5-30 (In Russian). https://doi.org/10.14529/mmp180401.
20. Bjôrck Â. Numerical Methods in Matrix Computations, Texts in Applied Mathematics, vol.59. Cham, Springer, 2015, xvi+800 pp. https://doi.org/10.1007/ 978-3-319-05089-8.
21. Oshman Y., Bar-Itzhack I. Y. Square root filtering via covariance and information eigenfactors, Automatica, 1986, vol.22, no. 5, pp. 599-604. https://doi.org/10.1016/ 0005- 1098(86)90070- 1.
22. Oshman Y. Square root information filtering using the covariance spectral decomposition, In: Proc. of the 27th Conf. on Decision and Control, 1988, pp. 382-387. https://doi.org/ 10.1109/CDC.1988.194335.
23. Oshman Y. Maximum likelihood state and parameter estimation via derivatives of the V-Lambda filter, J. Guid. Control Dyn., 1992, vol.15, no. 3, pp. 717-726. https://doi.org/ 10.2514/3.20896.
24. Wang L., Libert G., Manneback P. Kalman filter algorithm based on singular value decomposition, In: Proc. of the 31st Conf. on Decision and Control, 1992, pp. 1224-1229. https:// doi.org/10.1109/IEC0N.1992.254406.
25. Zhang Y., Dai G., Zhang H., Li Q. A SVD-based extended Kalman filter and applications to aircraft flight state and parameter estimation, In: Proc. of 1994 American Control Conf., 1994, pp. 1809-1813. https://doi.org/10.1109/ACC.1994.752384.
26. Kulikova M. V., Tsyganova J. V. Improved discrete-time Kalman filtering within singular value decomposition, IET Control Theory Appl., 2017, vol. 11, no. 15, pp. 2412-2418, arXiv: 1611.03686 [math.OC]. https://doi.org/10.1049/iet- cta.2016.1282.
27. Tsyganova J. V., Kulikova M. V. SVD-based Kalman filter derivative computation, IEEE Trans. Autom. Control, 2017, vol.62, no. 9, pp. 4869-4875, arXiv: 1612.04777 [cs.SY]. https://doi.org/10.1109/TAC.2017.2694350.
28. Alessandrini M., Biagetti G., Crippa P., Falaschetti L., Manoni L., Turchetti C. Singular value decomposition in embedded systems based on ARM Cortex-M architecture, Electronics, 2021, vol.10, no. 1, 34. https://doi.org/10.3390/electronics10010034.
