CC BY
24
УДК 004.94 Дата подачи статьи: 26.07.21, после доработки: 13.09.21
DOI: 10.15827/0236-235X.136.639-648 2021. Т. 34. № 4. С. 639-648
Программный комплекс для компьютерного моделирования процессов параметрической идентификации математических моделей конвективно-диффузионного переноса
А.Н. Кувшинова 1, аспирант, kuvanulspu@yandex.ru
А.В. Цыганов 1, к.ф.-м.н., профессор, andrew.tsygan.ov@gm.ail.com
1 Ульяновский государственный педагогический университет им. И.Н. Ульянова, г. Ульяновск, 430071, Россия
В работе описан программный комплекс для компьютерного моделирования процессов параметрической идентификации одномерных математических моделей конвективно-диффузионного переноса с постоянными коэффициентами в условиях зашумленных измерений. Реализованные в программном комплексе методы параметрической идентификации основаны на переходе от математической модели, описываемой уравнениями в частных производных с начальным и граничными условиями, к модели, описываемой линейной дискретной стохастической системой в пространстве состояний с зашумленными измерениями, с последующим применением к данной системе методов оптимальной дискретной фильтрации и параметрической идентификации.
Программный комплекс написан на языке МАТЬАВ и представляет собой набор функций и скриптов, предназначенных для построения конечно-разностной сетки, дискретизации исходной задачи, построения решения исходной задачи методом конечных разностей, моделирования экспериментальных данных, численной идентификации граничных условий, идентификации скорости конвекции и коэффициента диффузии, а также выполнения различных вспомогательных операций, таких как визуализация результатов, проведение численных экспериментов и др.
Идентификация граничных условий выполняется при помощи алгоритма С. Гиллийнса и Б. Мура одновременного оценивания вектора состояния и неизвестного вектора входных воздействий линейной динамической системы. Идентификация скорости конвекции и коэффициента диффузии выполняется при помощи минимизации критерия идентификации, взятого в виде логарифмической функции правдоподобия, которая строится на основе величин, вычисляемых фильтром Калмана. Минимизация критерия идентификации выполняется с помощью функций системы МАТЬАВ для безградиентной и градиентной оптимизации.
Приводятся подробное описание программного комплекса и результаты компьютерного моделирования, подтверждающие работоспособность алгоритмов.
Ключевые слова: конвективно-диффузионный перенос, граничные задачи, дискретная линейная стохастическая модель, параметрическая идентификация, оптимальное оценивание, MATLAB.
Для описания различных природных и технических процессов широко используются математические модели тепломассопе-реноса [1, 2]. Для данных моделей актуальными являются задачи идентификации их параметров по результатам измерений значений искомой функции в отдельных точках рассматриваемой области. Такого рода задачи называются обратными и в общем случае относятся к некорректно поставленным задачам математической физики [3]. Для практических приложений традиционными методами решения этих задач являются экстремальные постановки с последующей параметрической идентификацией, однако возможности применения данных методов в системах, работающих в реальном времени, существенно
ограничены по причине большого объема вычислений.
В ряде работ отечественных и зарубежных авторов для решения обратных задач тепло-массопереноса предлагается использовать рекуррентные методы параметрической идентификации, основанные на алгоритмах оптимальной дискретной фильтрации (преимущественно калмановского типа).
Данные методы показали свою эффективность при решении задач температурной диагностики двигателей и паровых турбин [4-7], нестационарной теплометрии [8-11], конвективной диффузии [12]. В работах [13, 14] приводятся примеры применения различных модификаций фильтра Калмана и фильтра частиц к задачам теплопереноса. Применение
данных методов требует проведения исследований относительно устойчивости, сходимости и оценки погрешностей получаемых результатов. В связи с этим актуальной является задача разработки ПО для проведения соответствующих исследований. Следует также заметить, что численная эффективность реализации рекуррентных методов параметрической идентификации в определенных случаях может существенно зависеть от эффективности программной реализации соответствующих методов оптимальной дискретной фильтрации.
В данной работе описывается программный комплекс, предназначенный для моделирования процессов параметрической идентификации математических моделей конвективно-диффузионного переноса, в котором реализованы алгоритмы параметрической идентификации, предложенные в работах [15-18].
Постановка задачи
Рассмотрим одномерную модель конвективно-диффузионного переноса с постоянными коэффициентами, заданную следующими уравнениями:
де де
--+ v— = а
dt дх
д 2е д2 х'
a < х < b, 0 < t <
е(х, 0) = ф(х), a < х < b, е (a, t)= f (t), 0 < t <+x>, c (b, t) = g (t), 0 < t <+x
или
де_ дх
(b, t) = -X[c(b, t)-g(t)], 0 <t <+»,
(1)
(2)
(3)
(4)
(4')
где х - пространственная координата; t -время; а и Ь - границы рассматриваемого отрезка; с(х, () - искомая функция, например, концентрация загрязняющего вещества в потоке жидкости в точке с координатой х в момент времени V, V - скорость конвекции; а -коэффициент диффузии; (2) - начальное условие; (3), (4), (4') - граничные условия. Таким образом, рассматривается задача либо с граничными условиями первого рода, либо со смешанными граничными условиями первого и третьего рода.
Задача параметрической идентификации модели (1)-(4) ((1)-(4')) состоит в определении коэффициентов V и а в уравнении (1), а также функций Д/') и g(t) в граничных условиях по результатам зашумленных измерений искомой функции в отдельных точках рас-
сматриваемого отрезка в различные моменты времени.
Дискретизация модели
Для решения поставленной задачи перейдем от исходной непрерывной модели к дискретной, описываемой линейной динамической системой в пространстве состояний с за-шумленными измерениями, которая в общем случае имеет следующий вид:
|Ч = -А-1 + Бк-1ик-1 +
гк = Нкск +4, к =1, 2,..., где ск е К" - вектор состояния системы; ик е Кг - вектор входных воздействий; 2к е К™ -вектор измерений; Ек е К"х", Бк е К"Хг, Ок е К"х«, Нк е К™х" - матрицы; ^к е К« и £к е К™ - шумы, образующие независимые нормально распределенные последовательности с нулевым математическим ожиданием и ковариационными матрицами Q > 0 и Я > 0.
Дискретизируя исходную модель в пространственно-временной области на сетке X = а + ¿Ах, ¿к = кЬ1, ¿ = 0,1,...,И, к = 0,1, ...,К, в случае граничных условий (3) и (4) получаем дискретную линейную динамическую систему:
"1
е c2
„k-1
n-2 cn-1
" ck "
ck
<t
cl2
ck-1
<
a2 a3 0 0 0 0
a a2 аъ • 0 0 0
0 a a2 0 0 0
а, а.
f-1
0 0
0 0 Ц,- 1
0 a3
k = 1, 2,...,К,
(5)
а в случае граничных условий (3) и (4'):
" ck" a2 a 0 0 0 0 1
t c 2 a a аъ 0 0 0 k-1 c2
i c3 0 a a 0 0 0 k-1 c3
i cn - 2 0 0 0 a a 0 i-1 cn-2
t cn -1 0 0 0 aj a a k-1 c-
i . Cn . 0 0 0 a4a2 a4a,_ i-1 Cn
ck ck-i
к
с
k-1
k-1
+
в
a 0"
0 0
0 0 y-r
0 0
0 0 uk-1
0 a5
k = 1, 2,
(5')
Системы (5) и (5') являются детерминированными (шум в объекте отсутствует) дискретными линейными динамическими системами, в которых граничные условия входят в двумерный (г = 2) вектор входных воздействий ик. Коэффициенты в матрицах систем (5) и (5') имеют следующий вид: а\ = п + Г2, а2 = 1 - 2г2, аз = Г2 - Г1, 1 ХАх
.а. =-
где г1 = -
1 + ХАх aAt .
—-, Ах - пространственный
1 + ХДх vДt
2Дх' 2 Дх шаг сетки; At - временной шаг сетки, а сами матрицы являются постоянными (Ек = Е, Вк = В).
В системе (5) п = N — 1 вектор состояния ск состоит из всех внутренних узлов пространственной сетки, а в системе (5') п = N вектор состояния ск состоит из всех внутренних узлов пространственной сетки и правой границы.
Добавим к уравнениям (5) и (5') модель измерителя в виде уравнения зашумленных измерений:
Г k 1 Z1 Г kl C1 "Si "
А = H < + £
k k Cn _ M.
к = I, 2,...,К,
(6)
где Н е КтХп - матрица измерений (структура измерителя); т - количество измеряемых компонент вектора состояния.
Для решения задачи параметрической идентификации полученных дискретных моделей (5), (6) и (5'), (6) может быть применен богатый арсенал методов оптимальной дискретной фильтрации и параметрической идентификации.
Описание программного комплекса
Программный комплекс реализован на языке МЛТЬЛБ в виде набора скриптов и функций, позволяющих решать задачи
- дискретизации модели и моделирования экспериментальных данных;
- идентификации граничных условий;
- идентификации коэффициентов конвекции и диффузии.
Рассмотрим эти задачи подробнее.
Дискретизация модели и моделирование экспериментальных данных. Для представления исходной модели используются следующие структуры: pde - коэффициенты уравнения, domain - область изменения переменных, ic - начальное условие, bc - граничные условия, bct - типы граничных условий, а для ее дискретизации: grid - пространственно-временная сетка, Ids - линейная динамическая система.
В структуре grid хранится информация о количестве узлов, узлах и шагах сетки, а в структуре Ids - матрицы F и B.
Построение сетки выполняется с помощью функции grid = compute_grid(pde, domain, nX, nT), где nX, nT - число узлов по соответствующим переменным (в случае отсутствия аргумента nT число узлов по переменной t вычисляется автоматически из условия сходимости конечно-разностной схемы).
Дискретизация исходной модели выполняется при помощи функции Ids = = pde2lds(pde, bc, bct, grid).
Для моделирования процесса измерений используется структура sens, хранящая информацию о матрицах H и R. Измеритель задается с помощью функции sens = sensor(bct, nX, variance, problemtype), где variance - дисперсия шума в измерителе; problemtype -строковый параметр, определяющий вид матрицы H.
Моделирование процесса измерений осуществляется путем последовательного решения исходного уравнения конечно-разностным методом и зашумления значений решения в выбранных точках. Для этой цели используются следующие функции:
solution = solve_findiffpde, ic, bc, bct, grid),
Z = measure(solution, bct, sens), в результате чего получаются численное решение исходной задачи solution в виде матрицы размера KxN и матрица результатов измерений Z е RmxK, в которой по столбцам хранятся векторы измерений zk, k = 1, 2, ..., K.
Рассмотрим процесс дискретизации исходной модели и моделирования процесса измерений на примере следующей задачи:
dc „ dc d2c
— + 2— =-z-,
dt dx d x
0 < x < 1,0 < t < +TO,
(7)
+
в
k-1
c (x, 0) = 0,0 < x < 1, (8)
c(0, t) = |sin3rct| ,0 < t < +», (9)
^(1, t) = -[c(1, t)-1],0< t <+». (10)
В данном случае рассматривается задача со смешанными граничными условиями вида (1), (2), (3), (4'), где v =2, а =1, ф(х) = 0, ft) = \ sin 3%t\, g(t) = t, l = 1.
Приведем содержимое скрипта run_prob-lem_demo с описанием модели (7)—(10):
pde.v = 2; pde.alpha = 1;
domain.a = 0; domain.b = 1; domain.Tmin = 0; domain.Tmax = 1;
ic.f = @ic_f; bc.fl = @bc_fl; bc.cl = 1.0; bc.f2 = @bc_f2; bc.c2 = 1.0;
bct.tl = 1; bct.t2 = 3;
nX = б; I
variance = 0.02A2;I
convection speed diffusion coefficient
left x bound right x bound start time % stop time
% initial condition function % left boundary condition function % left boundary condition constant right boundary condition function right boundary condition constant
% left boundary condition type % right boundary condition type
% number of x-grid points noise variance
grid = compute_grid(pde, domain, nX); lds = pde2lds(pde, bc, bct, grid); sens = sensor(bct, nX, variance, 'bc'); solution = solve_findiff(pde, ic, bc, bct, grid); Z = measure(solution, bct, sens);
print_grid(grid);
print_lds(lds);
print_sensor(sens);
plot_conditions(ic, bc, bct, grid); plot_solution(solution, grid); plot_measurements(bct, grid, sens, Z);
В скрипте выполняется дискретизация задачи на сетке с 6 пространственными узлами (количество узлов по переменной t вычисляется автоматически и равно 101), а затем моделируется процесс зашумленных измерений решения с дисперсией шумов 0.022.
Для корректной работы скрипта пользователем должны быть в отдельных файлах предварительно реализованы функции ic_f, bcf1 и bcf2 с описанием начального и граничных условий. Например, функция, описывающая правое граничное условие, реализована в файле bc_f1.m следующим образом:
function c = bc_f1(t) c = abs(sin(3*pi*t)); end
В результате работы скрипта построены: - конечно-разностная сетка
nX = 6 nT = 101 xGrid = 0:0.2:1 tGrid = 0:0.01:1
- линеиная динамическая система
F (5x5) :
0.5000 0.3000 0 0 0
B(5x2):
0.3000 0 0 0 0
0.2000 0.5000 0.3000 0 0
0 0 0 0
0.1бб7
0
0.2000 0.5000 0.3000 0.2500
0 0
0.2000 0.5000 0 ,4167
0 0 0
0.2000 0.1667
измеритель
R(2x2):
l.0e-03 *
0 0.4000
Графики решения, полученного методом конечных разностей и смоделированных измерений, приведены на рисунках 1 и 2 соответственно.
Идентификация граничных условий. В программном комплексе реализован метод численной идентификации функцийД?) и g(t),
0
входящих в граничные условия, на основе алгоритма С. Гиллийнса и Б. Мура одновременного оптимального оценивания вектора состояния и неизвестного вектора входных воздействий [19]. Особенностью данного алгоритма является то, что ему не требуется априорная информация о виде входных воздействий. Подробно применение данного алгоритма к задаче идентификации граничных условий модели конвективно-диффузионного переноса описано в работах [16, 17].
Алгоритм идентификации граничных условий реализован в виде функции [estSol, estBCF2, PHistory, PUHistory] = bc_id_gm(ic, bct, grid, Ids, sens, Z).
Для работы алгоритма необходимы предварительная дискретизация задачи и измерение значений решения в узлах пространственной сетки, соответствующих первому и последнему элементам вектора состояния (в случае условий первого рода это будут второй и предпоследний узлы, а в случае условий первого и третьего рода - второй и последний узлы сетки).
Выходными аргументами функции являются: estSol - оценка решения в узлах пространственно-временной сетки; estBCF2 -оценка значений функции в граничном условии третьего рода; PHistory - ковариации ошибок оценивания вектора состояния; PUHistory - ковариации ошибок оценивания вектора входных воздействий.
estSol представляет собой матрицу того же размера, что и матрица численного решения solution, полученного по методу конечных разностей. В случае, когда решается задача идентификации с двумя граничными условиями первого рода, их оценка сохраняется в первом и последнем столбцах estSol (вектор estBCF2 при этом является пустым), а в остальных столбцах хранятся оценки вектора состояния. Если решается задача с условиями первого и третьего рода, их оценки сохраняются соответственно в первом столбце estSol и векторе estBCF2, а в остальных столбцах хранятся оценки вектора состояния. В матрицах PHistory и PUHistory по столбцам хранятся ковариации ошибок оценивания вектора состояния и измерений соответственно, характеризующие качество получаемых оценок.
Для запуска процесса идентификации граничных условий имеется скрипт run_bcjd_gm, содержимое и параметры которого для задачи (7)-(10) аналогичны приведенному ранее
скрипту, но с добавлением после моделирования измерений вызова функции bcjdgm и визуализацией получаемых результатов. На рисунках 3-6 приведены некоторые результаты его работы.
Идентификация коэффициентов конвекции и диффузии. В программном комплексе реализованы несколько методов идентификации коэффициентов уравнения (1), описанных в работах [15] и [18]. Данные методы основаны на поиске значений параметров v и а, минимизирующих критерий идентификации, который задается в виде логарифмической функции правдоподобия и строится на величинах, вычисляемых фильтром Калмана: векторе невязок измерений и на его ковариационной матрице (см., например, [20]).
Для вычисления значений критерия идентификации реализованы две функции: LR = kf_lr(theta, bc, bct, grid, sens, phi, f1, f2, Z), [LR, LRG] = kf_lr_grad(theta, bc, bct, grid, sens, phi, f1, f2, Z).
Первый аргумент обеих функций - вектор значений параметров 0 = [v, а]. Вторая функция отличается от первой дополнительным вычислением градиента функции правдоподобия.
На рисунке 7 для задачи (7)-( 10) приведен результат выполнения скрипта exp_plot_kf_lr, предназначенного для построения и визуализации критерия идентификации в заданной пользователем области изменения переменных.
Минимизация критерия идентификации может выполняться различными функциями MATLAB, предназначенными для минимизации функций нескольких переменных. В программном комплексе для этих целей используются функции simulannealbnd (Simulated Annealing - метод имитации отжига), ga (Genetic Algorithm - генетический алгоритм) и fmincon (условная минимизация). Первые две функции являются безградиентными ме-таэвристическими методами оптимизации, и в них в качестве целевой функции используется функция kflr. Функция fmincon может использоваться как с градиентной (kf_lr_ grad), так и с безградиентной (kf lr) версией функции вычисления значений критерия идентификации (во втором случае градиент целевой функции аппроксимируется самой функцией fmincon).
Для идентификации коэффициентов уравнения в программном комплексе реализованы функции coef_id_ga, coef_id_sa,
Рис. 3. Оценка решения (estSol) Fig. 3. The solution estimate (estSol)
Рис. 4. Функция f(t) и ее оценка (estSol(:,1)) Fig. 4. f(t) and its estimate (estSol(:,1))
1 -pun
4 ---PU22
1
i '
! \
1 1.
!-I-i4 1
\
10 20 30 40
50 60 k
70 BO 90 100
Рис. 5. Функция g(t) и ее оценка (estBCF2) Fig. 5. g(t) and its estimate (estBCF2)
Рис. 6. Ковариации ошибок оценивания
граничных значений (вектора uk) Fig. 6. Covariances of errors in estimating boundary values (vector uk)
coefidfmincon и coefidfmincon grad с унифицированным интерфейсом. Приведем пример реализации функции coef_id_sa:
function theta = coef_id_sa(ic, bc, bct, grid, sens, saOptions, LB, UB, Z, theta0)
vLB = LB(1); vUB = UB(1); alphaLB = LB(2); alphaUB = UB(2);
phi = ic.f(grid.xGrid); % phi(x) f1 = bc.f1(grid.tGrid); % f(t) f2 = bc.f2(grid.tGrid); % g(t)
if nargin < 10
v0 = vLB + 0.5*(vUB-vLB);
alpha0 = alphaLB + 0.5*(alphaUB-alphaLB); theta0 = [v0 alpha0];
end
objFnc = @(theta) kf_lr(theta, bc, bct, grid, sens, phi, f1, f2, Z);
[theta, fval, exitflag, output] = simulan-nealbnd(objFnc, theta0, [vLB alphaLB], [vUB alphaUB], saOptions);
Описание исходной задачи, настройка параметров и вызов нужного алгоритма выполняются в скриптах ru"_coef_id_ga, ru"_coefJd_sa, ru"_coefJd_fmmco" и ru"_coefJdJmmco"_grad. Представим фрагмент скрипта ru"_coef_id_sa и результат его выполнения (идентифицированные значения параметров V и а):
vLB = 0; vUB = 5; alphaLB = 0; bound
alphaUB = 5; bound
convection speed lower bound convection speed upper bound diffusion coefficient lower
% diffusion coefficient upper saOptions = saoptimset('simulannealbnd');
saOptions = saoptimset(saOptions, 'TimeLimit', Inf);
saOptions = saoptimset(saOptions, 'MaxIter', Inf);
saOptions = saoptimset(saOptions, 'Reanneal-Interval', 100);
saOptions = saoptimset(saOptions, 'StallIterLim-it', 100);
saOptions = saoptimset(saOptions, 'MaxFunEvals', Inf);
saOptions = saoptimset(saOptions, 'Display', 'off') ;
saOptions = saoptimset(saOptions, 'DisplayInter-val', 1);
saOptions = saoptimset(saOptions, 'OutputFcn', @sa_out);
grid = compute_grid(pde, domain, nX); solution = solve_findiff(pde, ic, bc, bct, grid);
sens = sensor(bct, nX, variance, 'coef'); Z = measure(solution, bct, sens);
theta = coef_id_sa(ic, bc, bct, grid, sens, saOptions, [vLB alphaLB], [vUB alphaUB], Z);
disp(['v = ', num2str(theta(1),15) ', num2str(theta(2),15)]);
alpha
>> run_coef_id_sa v = 1.94311777853335, alpha
0.989708385847537
Заключение
В статье описан программный комплекс, предназначенный для моделирования процессов параметрической идентификации одномерных моделей конвективно-диффузионного переноса с постоянными коэффициентами, в котором реализованы методы параметрической идентификации коэффициентов и граничных условий математической модели.
Помимо основных функций параметрической идентификации и различных вспомога-
тельных функций, в программном комплексе реализованы скрипты, упрощающие решение задач параметрической идентификации, а также используемые для проведения серий экспериментов по исследованию статистических свойств реализованных алгоритмов (например, точности и времени идентификации).
Данный программный комплекс может служить удобным средством для разработки и исследования новых алгоритмов параметрической идентификации и оптимального оценивания моделей конвективно-диффузионного переноса на основе алгоритмов оптимальной дискретной фильтрации. В настоящее время авторами ведется работа по созданию и реализации новых алгоритмов параметрической идентификации.
Исследование выполнено при финансовой поддержке РФФИ и Ульяновской области в рамках научного проекта № 19-41-730009.
Литература
1. Леонтьев А.И. Теория тепломассообмена. М.: Изд-во МГТУ, 1997. 683 с.
2. Фарлоу С. Уравнения с частными производными для научных работников и инженеров; [пер. с англ.]. М.: Мир, 1985. 384 с.
3. Денисов А.М. Введение в теорию обратных задач. М.: Изд-во МГУ, 1994. 208 с.
4. Симбирский Д.Ф. Температурная диагностика двигателей. Киев: Техника, 1976. 208 с.
5. Мацевитый Ю.М., Мултановский А.В. Идентификация параметров теплообмена методом оптимальной динамической фильтрации // ТВТ. 1979. Т. 17. № 5. С. 1053-1060.
6. Карпов А.А., Тихонова Т.А. Восстановление нестационарных тепловых потоков по экспериментальным данным // Матем. моделирование. 2000. Т. 12. № 5. С. 101-106.
7. Симбирский Г.Д., Лантрат В.К. Применение цифрового фильтра Калмана для параметрической идентификации высокотемпературного термопреобразователя // Автомобиль и электроника. Современные технологии. 2017. № 11. С. 68-75.
8. Пилипенко Н.В. Применение фильтра Калмана в нестационарной теплометрии. СПб: Изд-во Университет ИТМО, 2017. 36 с.
9. Пилипенко Н.В., Заричняк Ю.П., Иванов В.А., Халявин А.М. Параметрическая идентификация дифференциально-разностных моделей теплопереноса в одномерных телах на основе алгоритмов фильтра Калмана // Научно-технический вестник информационных технологий, механики и оптики. 2020. Т. 20. № 4. С. 584-588. DOI: 10.17586/2226-1494-2020-20-4-584-588.
10. Daouas N., Radhouani M.-S. A new approach of the Kalman filter using future temperature measurements for nonlinear inverse heat conduction problems. Numerical Heat Transfer Fundamentals, 2004, vol. 45, no. 6, pp. 565-585. DOI: 10.1080/10407790490430598.
11. Pacheco C.C., Orlande H.R.B., Colago M.J., Dulikravich G.S. Identification of a position and time dependent heat flux by using the Kalman filter and improved lumped analysis in heat conduction. Proc. ICCM2014, 2014, Cambridge, England, pp. 801-809.
12. Матвеев М.Г., Копытин А.В., Сирота Е.А. Комбинированный метод идентификации параметров распределенной динамической модели // Сб. тр. IV Междунар. конф. ИТНТ. 2018. С. 1651-1657.
13. Myers M.R., Jorge A.B., Mutton M.J., Walker D.G. A comparison of extended Kalman filter, ultrasound time-of-flight measurement models for heating source localization. Inverse Problems in Science and Engineering, 2012, vol. 20, no. 7, pp. 991-1016. DOI: 10.1080/17415977.2012.669272.
14. Orlande H.R.B., Colago M.J., Dulikravich G.S., Vianna F.L.V. et al. State estimation problems in heat transfer. Int. J. for Uncertainty Quantification, 2012, vol. 2, no. 3, pp. 239-258. DOI: 10.1615/Int.J. UncertaintyQuantification.2012003582.
15. 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. Proc. II Int. Sci. and Pract. Conf. FTI, 2018, pp. 188-196.
16. Цыганов А.В., Цыганова Ю.В., Кувшинова А.Н. Динамическая идентификация граничных условий в модели конвективно-диффузионного переноса в условиях зашумленных измерений // Сб. тр. V Междунар. конф. ИТНТ. 2019. Т. 3. С. 169-177.
17. Кувшинова А.Н. Динамическая идентификация смешанных граничных условий в модели конвективно-диффузионного переноса в условиях зашумленных измерений // Журнал Средневолж-ского математического общества. 2019. Т. 21. № 4. С. 469-479. DOI: 10.15507/2079-6900.21.201904. 469-479.
18. Кувшинова А.Н., Цыганов А.В., Цыганова Ю.В., Тапиа Гарса У.Р. Алгоритм численной идентификации параметров в модели конвективно-диффузионного переноса // Сб. тр. VI Междунар. конф. ИТНТ. 2020. С. 825-832.
19. Gillijns S., Moor B.D. Unbiased minimum-variance input and state estimation for linear discrete-time systems. Automatica, 2007, vol. 43, no. 1, pp. 111-116. DOI: 10.1016/j.automatica.2006.08.002.
20. Astrom K.J. Maximum likelihood and prediction error methods. Automatica, 1980, vol. 16, no. 5, pp. 551-574.
Software & Systems Received 26.07.21, Revised 13.09.21
DOI: 10.15827/0236-235X.136.639-648 2021, vol. 34, no. 4, pp. 639-648
A software package for computer modeling of parametric identification processes for mathematical models of convection-diffusion transfer
A.N. Kuvshinova 1, Postgraduate Student, kuvanulspu@yandex.ru
A.V. Tsyganov 1, Ph.D. (Physics and Mathematics), Professor, andrew.tsyganov@gmail.com
1 Ulyanovsk State Pedagogical University named after I.N. Ulyanov, Ulyanovsk, 430071, Russian Federation
Abstract. The paper describes a software package for computer modeling of parameter identification processes of one-dimensional mathematical models of convection-diffusion transfer with constant coefficients under noisy measurement conditions. The parameter identification methods implemented in the software package are based on the transition from a mathematical model described by partial differential equations with initial and boundary conditions to a model described by a linear discrete stochastic system in the state space with noisy measurements followed by the application of optimal discrete filtering and parameter identification methods to this system.
The software package is written in MATLAB and is a set of functions and scripts for constructing a finite-difference grid, discretizing the original problem, constructing a solution to the original problem by
the finite difference method, modeling experimental data, numerical identification of boundary conditions, identification of the convection velocity and diffusion coefficient, as well as performing various auxiliary operations, such as: visualization of results, conducting numerical experiments, etc.
The identification of boundary conditions is performed using the algorithm of S. Gillijns and B. Moore for simultaneous estimation of the state vector and the unknown control vector of a linear dynamic system. The identification of the convection velocity and the diffusion coefficient is performed by minimizing the identification criterion taken in the form of a logarithmic likelihood function based on the values calculated by the Kalman filter. The minimization of the identification criterion is performed using the functions of the MATLAB system for gradient-free and gradient optimization.
The paper gives a detailed description of the software package and the results of computer modeling confirming the operability of the algorithms.
Keywords: convection-diffusion transfer, boundary problems, discrete linear stochastic model, parameter identification, optimal estimation, MATLAB.
Acknowledgements. The study was carried out with the financial support of the RFBR and the Ulyanovsk region within the framework of the scientific project no. 19-41-730009.
References
1. Leontyev A.I. Theory of Heat and Mass Transfer. Moscow, 1997, 683 p. (in Russ.).
2. Farlow S.J. Partial Differential Equations for Scientists and Engineers. John Wiley and Sons Publ., 1982, 429 p. (Russ. ed.: Moscow, 1985, 384 p.).
3. Denisov A.M. Introduction to the Theory of Inverse Problems. Moscow, 1994, 208 p. (in Russ.).
4. Simbirskiy D.F. Temperature Diagnostics of Engines. Kiev: Tekhnika, 1976, 208 p. (in Russ.).
5. Matsevity Yu.M., Multanovsky A.V. Identification of heat transfer parameters using the optimal dynamic filtration method. High Temperature, 1979, vol. 17, no. 5, pp. 1053-1060 (in Russ.).
6. Karpov A.A., Tikhonova T.A. Recovery of non-stationary heat flows from experimental data. Ma-tem. Mod., 2000, vol. 12, no. 5, pp. 101-106 (in Russ.).
7. Simbirsky G.D., Lantrat V.K. Application of the Kalman digital filter for parametric identification high-temperature thermocouple. Avtomobil i Elektronika. Sovremennye Tekhnologii, 2017, no. 11, pp. 68-75 (in Russ.).
8. Pilipenko N.V. Applying the Kalman Filter in Non-Stationary Heat Metering. St. Petersburg, 2017, 36 p. (in Russ.).
9. Pilipenko N.V., Zarichnyak Yu.P., Ivanov V.A., Halyavin A.M. parametric identification of differencial-difference models of heat transfer in on e-dimensional bodies based on Kalman filter algorithms. Sci. Tech. J. Inf. Technol. Mech. Opt., 2020, vol. 20, no. 4, pp. 584-588 (in Russ.). DOI: 10.17586/2226-14942020-20-4-584-588.
10. Daouas N., Radhouani M.-S. A new approach of the Kalman filter using future temperature measurements for nonlinear inverse heat conduction problems. Numerical Heat Transfer Fundamentals, 2004, vol. 45, no. 6, pp. 565-585. DOI: 10.1080/10407790490430598.
11. Pacheco C.C., Orlande H.R.B., Colago M.J., Dulikravich G.S. Identification of a position and time dependent heat flux by using the Kalman filter and improved lumped analysis in heat conduction. Proc. ICCM2014, 2014, Cambridge, England, pp. 801-809.
12. Matveev M.G., Kopytin A.V., Sirota E.A. Combined method for identifying the parameters of a distributed dynamic model. Proc. IVInt. Conf. ITNT, 2018, pp. 1651-1657 (in Russ.).
13. Myers M.R., Jorge A.B., Mutton M.J., Walker D.G. A comparison of extended Kalman filter, ultrasound time-of-flight measurement models for heating source localization. Inverse Problems in Science and Engineering, 2012, vol. 20, no. 7, pp. 991-1016. DOI: 10.1080/17415977.2012.669272.
14. Orlande H.R.B., Colago M.J., Dulikravich G.S., Vianna F.L.V. et al. State estimation problems in heat transfer. Int. J. for Uncertainty Quantification, 2012, vol. 2, no. 3, pp. 239-258. DOI: 10.1615/Int.J. UncertaintyQuantification.2012003582.
15. 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. Proc. II Int. Sci. andPract. Conf. FTI, 2018, pp. 188-196.
16. Tsyganov A.V., Tsyganova Yu.V., Kuvshinova A.N. Dynamic identification of boundary conditions for convection-diffusion transport model in the case of noisy measurements. Proc. VInt. Conf. ITNT, 2019, vol. 3, pp. 169-177 (in Russ.).
17. Kuvshinova A.N. Dynamic identification of mixed boundary conditions for convection-diffusion transport model in the case of noisy measurements. Middle Volga Mathematical Society J., 2019, vol. 21, no. 4, pp. 469-479 (in Russ.). DOI: 10.15507/2079-6900.21.201904.469-479.
18. Kuvshinova A.N., Tsyganov A.V., Tsyganova Yu.V., Tapia Garza H.R. Algorithm for numerical identification of parameters for convective-diffusion transport model. Proc. VI Int. Conf. ITNT, 2020, pp. 825-832 (in Russ.).
19. Gillijns S., Moor B.D. Unbiased minimum-variance input and state estimation for linear discrete-time systems. Automatica, 2007, vol. 43, pp. 111-116. DOI: 10.1016/j.automatica.2006.08.002.
20. Astrom K.J. Maximum likelihood and prediction error methods. Automatica, 1980, vol. 16, no. 5, pp. 551-574.
Для цитирования
Кувшинова А.Н., Цыганов А.В. Программный комплекс для компьютерного моделирования процессов параметрической идентификации математических моделей конвективно-диффузионного переноса // Программные продукты и системы. 2021. Т. 34. № 4. С. 639-648. DOI: 10.15827/0236-235X. 136.639-648.
For citation
Kuvshinova A.N., Tsyganov A.V. A software package for computer modeling of parametric identification processes for mathematical models of convection-diffusion transfer. Software & Systems, 2021, vol. 34, no. 4, pp. 639-648 (in Russ.). DOI: 10.15827/0236-235X.136.639-648.
