Reconstruction of the spatiotemporal reaction rate function in high-dimensional problems

Full Text

ВВЕДЕНИЕ

Обратные задачи для уравнений параболического типа, связанные с восстановлением неизвестных коэффициентов и функций по ограниченным данным наблюдений, являются классическими. Такие задачи возникают при моделировании процессов тепломассопереноса [1; 19; 20; 22], фильтрации, переноса примесей, а также в экологических приложениях [9; 13; 26], где необходимо определить функцию скорости реакции по данным о концентрации вещества.

Вопросы определения младшего коэффициента в параболических уравнениях исследовались в литературе в различных постановках. Существенная часть работ посвящена случаям, в которых неизвестный коэффициент зависит только от пространственных переменных и восстанавливается по условиям финального переопределения. Для таких задач получены результаты существования и единственности решений, а также установлены условия корректности, в ряде случаев доказательства опираются на принцип максимума и условия знакоопределенности искомого коэффициента [3; 4; 11; 12; 17; 21]. Наряду с этим рассматривались постановки, в которых коэффициент представляется в виде суммы известных функций с неизвестными параметрами, а дополнительная информация задается в виде значений решения в фиксированные моменты времени или в отдельных точках области [3; 14].

Интерес представляют задачи с локальными и интегральными условиями переопределения. В частности, исследовались одномерные постановки с данными типа Коши на части границы, а также задачи определения пространственно-зависящего коэффициента по интегральным наблюдениям [10; 16; 24]. Часть работ относится к случаям, когда младший коэффициент зависит только от времени. Для таких задач разработаны как теоретические подходы, так и численные методы восстановления, в том числе при одновременном определении нескольких неизвестных коэффициентов и источников [2; 5; 15; 18; 23; 25].

Задача идентификации коэффициента поглощения и скорости потребления метана по точечным данным уже была рассмотрена нами ранее [6; 8]. Однако рассмотренные ранее постановки не охватывают случай пространственно-временной зависимости искомого реакционного члена в задачах высокой размерности. При переходе к математическим моделям высокой размерности возрастает как вычислительная сложность решения прямой задачи, так и сложность идентификации неизвестной функции по ограниченным данным наблюдений. Данный подход представляет интерес для более полного описания сложных процессов. Задача восстановления пространственно-временной функции скорости реакции в четырехмерной области, в которой искомый реакционный член зависит от пространственно-временных переменных, является актуальной.

В работе [7] были получены результаты о корректности обратной задачи определения коэффициента поглощения по точечным данным в параболическом уравнении, доказаны существование и единственность решения в пространствах Соболева, получены априорные оценки. В настоящей работе этот подход развивается на случай восстановления пространственно-временной функции скорости реакции в конечномерном пространстве и построения численного алгоритма.

Рассматривается обратная задача восстановления пространственно-временной функции скорости реакции g = g(x, y, z, t), входящей в параболическое уравнение как младший коэффициент. Пусть G ⊂ R³ – ограниченная пространственная область с границей Г,

$Q=G\times \left(\mathrm{0,}T\right), S=Г\times \left(\mathrm{0,}T\right).$

Искомая функция u = u(x, y, z, t) удовлетворяет уравнению

$\begin{array}{c}{u}_t-Lu+g\left(x,y,z,t\right)u=f\left(x,y,z,t\right), \left(x,y,z,t\right)\in Q,\end{array}$ (1)

где оператор L имеет вид

$Lu=\sum _{i,j=1}^3{a}_{ij}\left(x,y,z,t\right)u_{x_i{x}_j}-\sum _{i=1}^3{a}_i\left(x,y,z,t\right)u_{x_i}-a_0\left(x,y,z,t\right)u.$

Здесь g(x, y, z, t) – неизвестная пространственно-временная функция скорости реакции, подлежащая восстановлению.

Уравнение (1) дополняется начально-краевыми условиями

$\begin{array}{c}{\left.u\right|}_{t=0}=u_0\left(x,y,z\right), {\left. Bu\right|}_S= g_0\left(x,y,z,t\right),\end{array}$ (2)

где

$Bu=\sum _{i=1}^3{\gamma }_i\left(x,y,z,t\right)\frac{\partial u}{\partial x_i}+\sigma \left(x,y,z,t\right)u.$

Условия переопределения имеют вид

$\begin{array}{c}u\left(x_i, y_i,z_i,t_i\right)={\psi }_i, i=\mathrm{1,2,}\dots ,r,\end{array}$ (3)

где (x_i, y_i, x_i, t_i) ∈ Ǭ.

Целью настоящей работы является построение и исследование численного алгоритма восстановления пространственно-временной функции скорости реакции g(x, y, z, t) по точечным данным. Неизвестная функция аппроксимируется конечномерной линейной комбинацией базисных функций, что позволяет свести исходную обратную задачу (1)–(3) к задаче определения конечного набора коэффициентов.

РЕЗУЛЬТАТЫ И ОБСУЖДЕНИЕ

Численный алгоритм решения прямой и обратной задачи

Рассматривается область

$\text{Ω}={\left(\mathrm{0,1}\right)}^3, Q=\text{Ω}\times \left(\mathrm{0,}T\right], T=1.$

В вычислительном эксперименте используется частный случай уравнения (1)

$u_t-\sum _{i=1}^3\frac{\partial }{\partial x_i}\left(k\left(x,y,z\right)\frac{\partial u}{\partial x_i}\right)+g\left(x,y,z,t\right)u=s\left(x,y,z,t\right), \left(x,y,z,t\right)\in Q.$

При этом k(x, y, z) = 1+x – коэффициент диффузии, в рассматриваемой модели он зависит только от координаты x, оставаясь постоянным по переменным y и z. На левой грани x = 0 задается условие Дирихле u(0, y, z, t) = 1, на правой грани x=1 задается краевое условие $k\left(x,y,z\right)\frac{\partial u}{\partial x}=1$

Дополнительная информация задается точечными измерениями

$d_i=u\left(x_i, y_i,z_i,t_i\right), i=\mathrm{1,2,}\dots ,M,$

где M = 4000. Точечные измерения определяются как значения решения прямой задачи. Искомая функция скорости реакции восстанавливается в виде конечной суммы радиально-базисных функций

$g_N\left(x,y,z,t\right)=\sum _{j=1}^N{c}_j{Ф}_j\left(x,y,z,t\right),$

где c_j ∈ R – неизвестные коэффициенты, а Ф_j – базисные функции. Используются радиально-базисные функции вида

${Ф}_j\left(X\right)=\max \left(\mathrm{0,} R^2-X-{{\text{Θ}}_j}_2^2\right), X=\left(x,y,z,t\right)\in {\left[\mathrm{0,1}\right]}^4,$

где Θ_j – центр базисной функции, R > 0 – радиус. Эмпирическим путем был получен базис с пространственными центрами ${\left\{\frac{1}{3},\frac{2}{3}\right\}}^3$, временными центрами {0,25, 0,5, 0,75}.

В рассматриваемом эксперименте используется N = 24 базисных функций с радиусом R = 0,6.

Целью эксперимента является построение и исследование алгоритма восстановления коэффициентов c_j и всей пространственно-временной функции скорости реакции по точечным данным.

Прямая задача дискретизировалась методом конечных элементов на регулярной сетке (6+1)×(6+1)×(6×1). Приближенное решение представляется в виде

$u_h\left(x,y,z,t_n\right)=\sum _{l=1}^{N_h}{U}_l\left(t\right){\phi }_l\left(x,y,z\right),$

где {φ_l} – базисные функции, а U(t) = (U₁(t), U₂(t), …, U_Nh(t))^T – вектор неизвестных коэффициентов.

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

$M\frac{dU}{dt}+\left(K+G\left(t,c\right)\right)U\left(t\right)=F\left(t\right),$

где M – матрица масс c элементами

$M_{il}=\left({\phi }_l, {\phi }_i\right)={\int }_{\text{Ω}}{\phi }_l{\phi }_{i}dx, i,l=\mathrm{1,2,}\dots ,N_h,$

K – матрица диффузии с элементами

$K_{il}=\sum _{m=1}^3{\int }_{\text{Ω}}k\left(x,y,z\right)\frac{\partial {\phi }_l}{\partial x_m}\frac{\partial {\phi }_i}{\partial x_m}dx, i,l=\mathrm{1,2,}\dots ,N_h,$

G(t,c) – матрица реакционного члена с элементами

$G_{il}\left(t,c\right)={\int }_{\text{Ω}}g\left(x,y,z,t\right){\phi }_l{\phi }_{i}dx, i,l=\mathrm{1,2,}\dots ,N_h,$

а F(t) = (F₁(t), F₂(t), …, F_Nh(t))^T – вектор правой части с компонентами

$F_i\left(t\right)={\int }_{\text{Ω}}s\left(x,y,z,t\right){\phi }_{i}dx, i=\mathrm{1,2,}\dots ,N_h.$

Для аппроксимации по времени использовалась неявная схема Эйлера с шагом ∆t = 0,05, что соответствует N_t = 20 временным слоям на интервале 0 ≤ t ≤ 1. Обозначая через Uⁿ = (U₁ⁿ, U₂ⁿ, ..., Uⁿ_Nh) ≈ U(t_n) вектор коэффициентов на временном слое t_n, получаем полностью дискретную систему

$M\frac{U^n-U^{n-1}}{∆t}+\left(K+G_n\left(c\right)\right)U^n=F^n, n=\mathrm{1,2,}\dots ,N_t.$

В качестве g использовалась гладкая функция

$g_{true}\left(x,y,z,t\right)=2\sin \left(\pi x\right)\sin \left(\pi y\right)\sin \left(\pi z\right)\left(1+\mathrm{0,5}\sin \left(\pi t\right)\right)+\mathrm{0,5,}$

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

Вводится априорный вектор коэффициентов c ̂. Далее в качестве искомой величины рассматривается отклонение неизвестного вектора c от априорного приближения c ̂

$\delta =c-\widehat{c}.$

При каждом значении δ строится очередное приближение коэффициента реакции

$g_{\delta }\left(x,y,z,t\right)=\sum _{j=1}^N\left({\widehat{c}}_j+{\delta }_j\right){\text{Φ}}_j\left(x,y,z,t\right),$

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

${\widehat{d}}_i\left(\delta \right)=u\left(t_i, x_i, y_i, z_i;g_{\delta }\right), i=\mathrm{1,2,}\dots ,M.$

Тем самым задача восстановления коэффициента сводится к конечномерной задаче оптимизации по вектору δ ∈ R_N.

В качестве целевого функционала используется взвешенный функционал

$J\left(\delta \right)=\sum _{i=1}^M[\sqrt{w_i^{val}}\frac{b_i}{s_i}{\left({\widehat{d}}_i\left(\delta \right)-d_i)\right]}^2+\lambda {\delta }_2^2.$

Здесь $w_i^{val}={\left(\frac{\left|d_i\right|}{ma{x}_k\left|d_k\right|}\right)}^2$ – вес, отражающий относительную значимость измерения, $b_i={\left(\sum _{j\in L\left(t_i\right)}{w}_j^{val}\right)}^{-1/2}$ – множитель балансировки по временным слоям, $s_i=\text{max}\left(\left|d_i\right|{\mathrm{,10}}^{-6}ma{x}_k\left|d_k\right|\right)$ – нормирующий коэффициент, предотвращающий неустойчивость при малых значениях данных, L(t_i) – множество измерений на том же временном слое, что и точка t_i, а λ = 0,05 – параметр регуляризации.

Минимизация функционала J(δ) выполняется методом L-BFGS-B (Limited-memory Broyden – Fletcher – Goldfarb – Shanno with Bounds), относящимся к классу квазиньютоновских методов оптимизации. Выбор данного метода обусловлен его эффективностью при минимизации гладких функционалов большой размерности, а также возможностью учитывать априорные ограничения на параметры. Метод использует градиентную информацию и строит приближение матрицы вторых производных, что обеспечивает более высокую сходимость по сравнению с классическими градиентными методами.

После завершения оптимизации восстанавливаются коэффициенты

$c^{*}=\widehat{c}+{\delta }^{*}$

и искомая функция

$g_{rec}\left(x,y,z,t\right)=\sum _{j=1}^N{c}_j^{*}{Ф}_j\left(x,y,z,t\right).$

Численные эксперименты

Численный эксперимент был проведен на компьютере со следующими характеристиками: процессор AMD Ryzen 5 8600G, оперативная память 32 ГБ, операционная система Windows 10.

Вычислительный эксперимент проводился на регулярной сетке (6+1)×(6+1)×(6+1) при использовании M = 4000 искусственных измерений, распределенных по временным слоям и внутренним слоям по глубине.

Для иллюстрации используемой пространственно-временной радиально-базисной параметризации на рисунке 1(a) приведена геометрия расположения базисных функций на фиксированном пространственном срезе. На рисунке 1(б) показан пример одной базисной функции.

Рисунок 1. Радиально-базисная параметризация: a) пространственная геометрия радиально-базисной параметризации (срез z = 0,5); б) пример базисной функции

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

В результате проведенного численного эксперимента было получено восстановление пространственно-временной функции скорости реакции g(x, y, z, t) и концентрации на координатной сетке.

В таблице 1 приведены показатели ошибки восстановления концентрации в точках измерений на каждом временном слое. Минимальные ошибки восстановления концентрации достигаются на промежуточных временных слоях, тогда как к финальному моменту времени ошибка возрастает. Это можно объяснить накоплением погрешности при пошаговом продвижении по времени, а также уменьшением чувствительности концентрации на координатной сетке к локальным изменениям искомой функции скорости реакции на поздних стадиях процесса. Для количественной оценки точности использовалась абсолютная взвешенная ошибка E_L2w и регуляризированная взвешенная относительная ошибка rMAPE. Обе метрики применялись для сопоставления истинной концентрации c и восстановленной концентрации c ̂. Абсолютная взвешенная ошибка вычислялась следующим образом

$E_{L_{2}w}\left(c,\widehat{c}\right)={\left({\sum }_{i=1}^M{w}_i{\left(c_i-{\widehat{c}}_i\right)}^2\right)}^{\frac{1}{2}},$

здесь w_i ≥ 0 – вес соответствующей точки сравнения. Данная величина характеризует абсолютное отклонение восстановленного решения от эталонного в среднеквадратичном смысле с учетом возможного перевзвешивания отдельных наблюдений. Регуляризированная взвешенная средняя относительная ошибка rMAPE вычислялась следующим образом:

$rMAP{E}_w\left(c, \widehat{c}\right)=100\frac{{\displaystyle {\sum }_{i=1}^M{w}_i\frac{\left|c_i-{\widehat{c}}_i\right|}{\left|c_i\right|+\alpha s_c}}}{{\displaystyle {\sum }_{i=1}^M{w}_i}}, \alpha =\mathrm{0,05}$

где величина s_c = median{|c_i| : |c_i| > ε} – медиана ненулевых по модулю значений эталонной концентрации. Введение члена αs_c добавляет регуляризацию, предотвращая неустойчивое поведение относительной ошибки в тех областях, где c_i близко к нулю.

Таблица 1. Ошибки восстановления концентрации в точках измерений

Суммарная ошибка по концентрации в точках измерений составляет 2,6845 %. Двумерное отображение концентрации u(x, y, z, t) показано на рисунке 2.

Рисунок 2. Двумерное отображение концентрации

Представлено отображение концентрации из четырехмерного пространства переменных в двумерное пространство признаков. Это изображение следует интерпретировать как низкоразмерную проекцию полного 4D-распределения для анализа структуры данных, близости состояний и взаимного расположения пространственно-временных характеристик. Понижение размерности выполнялось с помощью метода UMAP (Uniform Manifold Approximation and Projection). Данный метод предназначен для построения низкоразмерного представления данных с сохранением их локальной геометрической структуры. В пространстве меньшей размерности ищется такое расположение точек, которое в наибольшей степени сохраняет структуру исходных. Следует отметить, что наименование осей не имеет непосредственного физического смысла и не соответствует напрямую пространственным координатам или времени. Они представляют собой координаты низкоразмерного вложения, возникающего в результате работы алгоритма UMAP. Интерпретация изображения основывается не на абсолютных значениях этих координат, а на взаимном расположении точек, форме и степени совпадения структур на изображениях для исходного и восстановленного распределений.

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

Относительная ошибка восстановления функции скорости реакции g(x, y, z, t) составила rMAPE=4,4 %. Двумерное отображение функции скорости реакции g(x, y, z, t) показано на рисунке 3.

Рисунок 3. Двумерное отображение функции скорости реакции

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

ЗАКЛЮЧЕНИЕ И ВЫВОДЫ

В настоящей работе рассмотрена обратная задача восстановления пространственно-временной функции скорости реакции, входящей в параболическое уравнение в качестве младшего коэффициента. Для численного решения задачи искомая функция аппроксимировалась конечномерной линейной комбинацией радиально-базисных функций в четырехмерном пространстве переменных (x, y, z, t), что позволило свести исходную обратную задачу к задаче определения конечного набора коэффициентов.

Построен численный алгоритм, включающий решение прямой задачи методом конечных элементов, параметризацию искомой функции с помощью радиально-базисных функций. Проведена минимизация методом L-BFGS-B. Проведенный эксперимент показал, что предложенная схема позволяет восстановить как значение концентрации в точках измерений, так и пространственно-временную функцию скорости реакции. Суммарная относительная ошибка восстановления концентрации составила 2,6845 %, а относительная ошибка восстановления функции g(x, y, z, t) – 4,4 %.

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

About the authors

Julia A. Tukmacheva

Author for correspondence.
Email: y_tukmacheva@ugrasu.ru

Postgraduate student, Engineering School of Digital Technologies

References

Supplementary files

Supplementary Files

Action

1. JATS XML

Download

2. Figure 1. Radial-basis parameterization: a) spatial geometry of radial-basis parameterization (cut z = 0.5); b) example of basis function

Download (304KB)

Indexing metadata

3. Figure 2. Two-dimensional display of concentration

Download (138KB)

Indexing metadata

4. Figure 3. Two-dimensional display of the reaction rate function

Download (133KB)

Indexing metadata