Идентификация параметров модели конвекции–диффузии–реакции и неизвестных граничных условий при наличии случайных помех в измерениях

Полный текст

Введение и постановка задачи

Задачи параметрической идентификации относятся к классу обратных задач математического моделирования и заключаются в определении неизвестных параметров математической модели объекта, принадлежащей выбранному классу моделей, по известным входным сигналам и выходным данным измерений. Одними из основных подходов к их решению являются метод подпространств (subspace identification method) и метод минимума ошибки предсказания (minimum prediction-error (MPE) method). Фундамент данных подходов был заложен еще в середине прошлого века в работах [1, 2]. Базовые концепции второго подхода определены в работах Л. Льюнга [3, 4], в рамках которого конкретные методы идентификации получаются как частные случаи в зависимости от выбора функции потерь, критерия качества, структуры модели и численного метода минимизации. В нашей стране большой вклад в развитие теории идентификации систем внес выдающийся советский ученый Я. З. Цыпкин [5].

В классе дискретных линейных стохастических систем с известными входными сигналами методы идентификации неизвестного модельного параметра, от которого зависят системные матрицы, хорошо развиты. Классическим подходом к решению задачи параметрической идентификации модельного параметра $\theta$ одновременно с оцениванием вектора состояния дискретной линейной стохастической системы по известным входным и выходным данным является применение методов адаптивной фильтрации [6]. При этом критерий идентификации определяет разницу между выходами реальной физической системы и дискретной линейной стохастической модели с присоединенным фильтром Калмана, параметры которого необходимо настроить таким образом, чтобы численное значение критерия идентификации стало оптимальным [3]. В настоящее время существуют различные подходы к построению алгоритмов параметрической идентификации, требующие построения адаптивных фильтров. В общий класс методов минимума ошибки предсказания входят хорошо известные метод максимального правдоподобия [7] и метод наименьших квадратов [8], а также методы минимаксного оценивания [9]. Отдельно отметим метод вспомогательного функционала качества [10, 11].

Сходимости оценок неизвестных параметров как для линейных, так и для нелинейных стохастических систем посвящены, в частности, работы [12–14]. Общие условия сходимости оценок неизвестных параметров при решении задач параметрической идентификации линейных стохастических систем рассмотрены в работе Л. Льюнга [12]. При решении задач параметрической идентификации предположения о выполнении условий сходимости, среди которых полная управляемость и наблюдаемость дискретной линейной стохастической системы, а также непрерывная дифференцируемость по модельному параметру элементов системных матриц и компактность множества допустимых значений модельного параметра, обеспечивают применимость численных методов для минимизации критерия параметрической идентификации.

Отметим, что для квадратических критериев качества, являющихся выпуклыми и монотонно возрастающими функциями по искомому параметру, существует глобальный минимум. Следовательно, отыскание точки минимума представляется возможным с помощью подходящих численных методов [15].

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

Для систем с дискретным временем самые ранние подходы основаны на включении неизвестного вектора входных сигналов в вектор состояния системы. Предполагалось, что модель динамики вектора входных сигналов известна. В этом случае для решения задачи применялся расширенный фильтр Калмана.

Большой вклад в развитие теории дискретной фильтрации стохастических систем с неизвестными входными сигналами внесли С. Гиллейнс и Б. Де-Мор. В [16] они расширили результаты, полученные в [17, 18], и предложили рекуррентный алгоритм одновременного оценивания вектора состояния системы и вектора неизвестных входных сигналов. При этом полученные оценки имеют минимальную дисперсию ошибки. Кроме того, они доказали оптимальность вычисляемых оценок. Перечисленные результаты относятся к решению задач оценивания вектора состояния и неизвестных входных воздействий. Во всех постановках задач предполагалось, что системные матрицы, определяющие математическую модель объекта и измерителя, точно известны.

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

Практическое применение полученного решения поставленной задачи показано на примере идентификации параметров математических моделей конвекции–диффузии–реакции, которые относятся к моделям тепломассопереноса и широко применяются при исследовании природных и техногенных процессов [19]. Для данного класса моделей актуальной является задача идентификации как параметров самой модели, так и входящих в нее граничных условий по результатам измерений значений искомой функции в отдельных точках рассматриваемой области. Усложняет задачу наличие случайных помех в измерениях.

Применение рекуррентных методов калмановского типа для решения задачи параметрической идентификации моделей, описываемых уравнениями в частных производных, рассмотрено в работах [20–23]. Рекуррентные методы решения задач параметрической идентификации моделей конвекции–диффузии предложены в работах [24–27].

Развивая и дополняя полученные ранее результаты [24–29], данная работа ставит целью построение комбинированного двухэтапного метода идентификации моделей тепломассопереноса, основанного на совместном применении методов минимизации критерия идентификации безградиентного типа и рекуррентного оценивания неизвестных входных сигналов.

Рассмотрим математическую модель конвекции–диффузии–реакции, заданную уравнением
\[ \begin{equation}
\dfrac{\partial c(x,t)}{\partial t}+v\dfrac{\partial c(x,t)}{\partial x} = \alpha\dfrac{\partial^{2} c(x,t)}{\partial x^{2}}-\beta c(x,t)
\end{equation} \tag{1} \]
с начальным условием
\[ \begin{equation}
c(x,0) = \varphi(x)
\end{equation} \tag{2} \]
и граничными условиями третьего рода
\[ \begin{equation}
\dfrac{\partial c(a,t)}{\partial x} = \lambda [c(a,t)-f(t)],\quad
\dfrac{\partial c(b,t)}{\partial x} = -\lambda [c(b,t)-g(t)]
\end{equation} \tag{3} \]
либо смешанными граничными условиями первого и третьего рода
\[ \begin{equation}
c(a,t) = f(t),\quad
\dfrac{\partial c(b,t)}{\partial x} = -\lambda [c(b,t)-g(t)]
\end{equation} \tag{4} \]
или
\[ \begin{equation}
\dfrac{\partial c(a,t)}{\partial x} = \lambda [c(a,t)-f(t)],\quad
c(b,t) = g(t),
\end{equation} \tag{5} \]
где $x\in[a;b]$ — пространственная координата; $t\in[0;T]$ — время; $c(x,t)$ — искомая функция; $v$ — скорость конвекции; $\alpha$ — коэффициент диффузии; $\beta$ — коэффициент реакции; $\varphi (x)$, $f(t)$, $g(t)$ — заданные функции.

В уравнении (1) $c(x,t)$ может иметь смысл температуры, концентрации вещества и т.д. Поставим задачу идентификации параметров данной модели, состоящую в определении коэффициентов $v$, $\alpha$, $\beta$ и $\lambda$ при наличии случайных помех в измерениях значений функции $c(x,t)$ в отдельных точках рассматриваемого отрезка в различные моменты времени. При этом функции $f(t)$ и $g(t)$, входящие в граничные условия (3), (4) или (5), считаются неизвестными и также подлежат идентификации. Функция $\varphi(x)$ в начальном условии (2) предполагается известной.

В первом разделе предложено новое решение задачи параметрической идентификации моделей тепломассопереноса, которое заключается в переходе к дискретной линейной стохастической модели в пространстве состояний и разработке комбинированного двухэтапного метода идентификации, основанного на последовательном применении методов минимизации критерия идентификации безградиентного типа и рекуррентного оценивания неизвестных входных сигналов. Во втором разделе построены новые дискретные линейные стохастические модели конвекции–диффузии–реакции для трех разных типов граничных условий. В третьем разделе изложен этап I — идентификация параметров дискретной стохастической модели. Четвертый раздел содержит описание этапа II — рекуррентное оценивание неизвестных граничных условий модели. В пятом разделе представлены результаты компьютерного моделирования на языке MATLAB, а также их обсуждение. Заключение, в котором перечислены полученные результаты, завершает статью.

1. Комбинированный метод идентификации неизвестных параметров и неизвестных граничных условий моделей тепломассопереноса

Математическая модель конвекции–диффузии–реакции, заданная уравнениями (1)–(5), имеет параметрическую неопределенность двух принципиально различных типов:

1. неизвестный векторный параметр $\theta=(v,\alpha,\beta,\lambda)^{\top}\in\mathbb{R}^4$, элементы которого являются числовыми константами;
2. неизвестные функции $f(t)$ и $g(t)$, входящие в граничные условия (3)–(5).

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

Рис. 1. Схема процесса идентификации [Figure 1. Identification process diagram]

Общая схема процесса идентификации состоит из одного подготовительного и двух основных этапов. Априорные данные о модели включают в себя непрерывную параметризованную модель $M(\theta)$, заданную уравнениями в частных производных, измерительную информацию $Z$, которая представляет собой измеренные датчиками отдельные значения целевой функции $c(x,t)$ в заданных точках пространства, причем измерения являются неточными ввиду наличия случайных помех.

На подготовительном этапе с помощью метода конечных разностей осуществляется переход от исходной непрерывной модели $M(\theta)$ к дискретной модели $M_D(\theta)$, представленной разностными уравнениями в пространстве состояний.

На первом этапе при сборе измерительной информации в качестве граничных условий используется известный тестовый сигнал. Затем по доступным данным измерений $Z$ и тестовым сигналам $U$ проводится идентификация векторного параметра $\theta$ с помощью численной минимизации критерия качества идентификации $J$. По завершении первого этапа найденная оценка $\hat\theta^\star$ используется для настройки модели $M_D(\theta)$.

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

Таким образом, по завершении процесса идентификации математической модели мы получаем:

1. оценку параметра $\theta$ в соответствии с выбранным критерием идентификации $J$;
2. численные значения функций $f(t)$ и $g(t)$, входящих в граничные условия;
3. оценки значений целевой функции $c(x,t)$ в каждой точке пространственно-временной сетки в рассматриваемой области.

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

2. Подготовительный этап: переход от исходной модели к дискретной линейной динамической модели в пространстве состояний

Перейдем от исходной непрерывной модели к дискретной модели, представленной линейной динамической системой в пространстве состояний. Следуя [25, 26], зададим в рассматриваемой пространственно-временной области конечно-разностную сетку $\{(x_i,t_k) \mid i=0,1,\ldots,N, \; k=0,1,\ldots,K\}$, где
\[ \begin{equation}
x_i=a+i\Delta x,\quad t_k=k\Delta t, \quad \Delta x = ({b-a})/{N}, \quad \Delta t = {T}/{K}.
\end{equation} \tag{6} \]
Введем обозначения: $c_i^k=c(x_i,t_k)$, $\varphi_i=\varphi(x_i)$, $f^k=f(t_k)$, $g^k=g(t_k)$. Заменим частные производные в уравнении (1) их конечно-разностными аппроксимациями, используя двухслойную по времени и симметричную по пространственной переменной схему. Предполагая, что соответствующие условия устойчивости для разностной схемы выполняются [30], получим следующую систему уравнений:
\[ \begin{multline}
\dfrac{c^{k}_i-c^{k-1}_i}{\Delta t}+v\dfrac{c^{k-1}_{i+1}-c^{k-1}_{i-1}}{2\Delta x} = {}
\\
{} = \alpha\dfrac{c^{k-1}_{i+1}-2c^{k-1}_i+c^{k-1}_{i-1}}{{\Delta x}^2}-\beta c^{k-1}_i,\quad
k=1,2,\ldots,K,
\end{multline} \tag{7} \]
где $i=0,1,\ldots,N$ для граничных условий (3); $i=1,2,\ldots,N$ для условий (4); $i =0,1,\ldots,N-1$ для (5).

Обозначим $r_1 = {\alpha\Delta t}({{\Delta x}^2}$, $r_2 = \beta\Delta t$, $r_3 = {v\Delta t}/({2\Delta x})$, $r_4 = {2\Delta t}/{{\Delta x}^2}$ и выразим из (7) $c^k_i$:
\[ \begin{equation*}
c^k_i = (r_1+r_3)c^{k-1}_{i-1} + (1-r_2-r_4)c^{k-1}_{i} + (r_1-r_3)c^{k-1}_{i+1} = a_1c^{k-1}_{i-1} + a_2c^{k-1}_{i} + a_3c^{k-1}_{i+1},
\end{equation*} \]
где $a_1 = r_1+r_3$, $a_2 = 1-r_2-r_4$, $a_3 = r_1-r_3$.

С учетом (6) для начального условия (2) получаем
\[ \begin{equation*}
c_i^0=\varphi_i,\quad i=0,1,\ldots,N.
\end{equation*} \]

Дискретизируя граничные условия (3), получаем
\[ \begin{equation*}
\frac{c^k_1-c^k_0}{\Delta x} = \lambda[c^k_0-f^k],\quad
\frac{c^k_N-c^k_{N-1}}{\Delta x} = -\lambda[c^k_N-g^k], \quad k=1,\ldots,K.
\end{equation*} \]
Обозначим $a_4 = {1}/({1+\lambda\Delta x})$, $a_5 = {\lambda\Delta x}/({1+\lambda\Delta x})$. Тогда из первого граничного условия получаем
\[ \begin{equation*}
c^k_0 = a_4c^k_1+a_5f^k = a_4a_1c^{k-1}_0+a_4a_2c^{k-1}_1+a_4a_3c^{k-1}_2+a_5f^k, \quad k=1,\ldots,K,
\end{equation*} \]
а из второго —
\[ \begin{equation*}
c^k_N = a_4a_1c^{k-1}_{N-2}+a_4a_2c^{k-1}_{N-1}+a_4a_3c^{k-1}_N+a_5g^k, \quad k=1,\ldots,K.
\end{equation*} \]

Для случая смешанных граничных условий (4) и (5) находим соответственно выражения
\[ \begin{equation*}
c^k_0 = f^k,\quad
\frac{c^k_N-c^k_{N-1}}{\Delta x} = -\lambda[c^k_N-g^k], \quad k=1,\ldots,K,
\end{equation*} \]
и
\[ \begin{equation*}
\frac{c^k_1-c^k_0}{\Delta x} = \lambda[c^k_0-f^k],\quad
c^k_N = g^k, \quad k=1,\ldots,K,
\end{equation*} \]
откуда окончательно получаем
\[ \begin{equation*}
c^k_0 = f^k,\quad
c^k_N = a_4a_1c^{k-1}_{N-2}+a_4a_2c^{k-1}_{N-1}+a_4a_3c^{k-1}_N+a_5g^k, \quad k=1,\ldots,K,
\end{equation*} \]
и
\[ \begin{equation*}
c^k_0 = a_4a_1c^{k-1}_0+a_4a_2c^{k-1}_1+a_4a_3c^{k-1}_2+a_5f^k,\quad
c^k_N = g^k, \quad k=1,\ldots,K.
\end{equation*} \]

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

1) дискретная линейная динамическая система для модели (1), (2), (3):
\[ \begin{multline}\label{eq-lds-33}
\begin{bmatrix}c^k_0\\c^k_1\\c^k_2\\\vdots\\c^k_{N-2}\\c^k_{N-1}\\c^k_{N}\end{bmatrix}
=
\begin{bmatrix}
a_4a_1 \!\! &\!\! a_4a_2 \!\! &\!\! a_4a_3 \!\!&\ldots&0&0&0\\
a_1&a_2&a_3&\ldots&0&0&0\\
0&a_1&a_2&\ldots&0&0&0\\
\vdots&\vdots&\vdots&\ddots&\vdots&\vdots&\vdots\\
0&0&0&\ldots&a_2&a_3&0\\
0&0&0&\ldots&a_1&a_2&a_3\\
0&0&0&\ldots& \!\!a_4a_1 \!\!& \!\!a_4 a_2 \!\!& \!\!a_4a_3
\end{bmatrix}
\hspace{-1mm}
\begin{bmatrix}c^{k-1}_0\\c^{k-1}_1\\c^{k-1}_2\\\vdots\\c^{k-1}_{N-2}\\c^{k-1}_{N-1}\\c^{k-1}_{N}\end{bmatrix}
+
\begin{bmatrix}a_5&0\\0&0\\0&0\\\vdots&\vdots\\0&0\\0&0\\0&a_5\end{bmatrix}
\hspace{-0.5mm}
\begin{bmatrix}f^{k}\\g^{k}\end{bmatrix},\\
k=1,2,\ldots,K;
\end{multline} \tag{8} \]

2) дискретная линейная динамическая система для модели (1), (2), (4):
\[ \begin{multline}
\begin{bmatrix}c^k_1\\c^k_2\\c^k_3\\\vdots\\c^k_{N-2}\\c^k_{N-1}\\c^k_{N}\end{bmatrix}
=
\begin{bmatrix}
a_2&a_1&0&\ldots&0&0&0\\
a_1&a_2&a_3&\ldots&0&0&0\\
0&a_1&a_2&\ldots&0&0&0\\
\vdots&\vdots&\vdots&\ddots&\vdots&\vdots&\vdots\\
0&0&0&\ldots&a_2&a_3&0\\
0&0&0&\ldots&a_1&a_2&a_3\\
0&0&0&\ldots&a_4a_1&a_4a_2&a_4a_3
\end{bmatrix}
\begin{bmatrix}c^{k-1}_1\\c^{k-1}_2\\c^{k-1}_3\\\vdots\\c^{k-1}_{N-2}\\c^{k-1}_{N-1}\\c^{k-1}_{N}\end{bmatrix}
+
\begin{bmatrix}a_1&0\\0&0\\0&0\\\vdots&\vdots\\0&0\\0&0\\0&a_5\end{bmatrix}
\begin{bmatrix}f^{k-1}\\g^{k}\end{bmatrix},\\
k=1,2,\ldots,K;
\end{multline} \tag{9} \]

3) дискретная линейная динамическая система для модели (1), (2), (5):
\[ \begin{multline}
\begin{bmatrix}c^k_0\\c^k_1\\c^k_2\\\vdots\\c^k_{N-3}\\c^k_{N-2}\\c^k_{N-1}\end{bmatrix}
=
\begin{bmatrix}
a_4a_1&a_4a_2&a_4a_3&\ldots&0&0&0\\
a_1&a_2&a_3&\ldots&0&0&0\\
0&a_1&a_2&\ldots&0&0&0\\
\vdots&\vdots&\vdots&\ddots&\vdots&\vdots&\vdots\\
0&0&0&\ldots&a_2&a_3&0\\
0&0&0&\ldots&a_1&a_2&a_3\\
0&0&0&\ldots&0&a_1&a_2
\end{bmatrix}
\begin{bmatrix}c^{k-1}_0\\c^{k-1}_1\\c^{k-1}_2\\\vdots\\c^{k-1}_{N-3}\\c^{k-1}_{N-2}\\c^{k-1}_{N-1}\end{bmatrix}
+
\begin{bmatrix}a_5&0\\0&0\\0&0\\\vdots&\vdots\\0&0\\0&0\\0&a_3\end{bmatrix}
\begin{bmatrix}f^{k}\\g^{k-1}\end{bmatrix},\\
k=1,2,\ldots,K.
\end{multline} \tag{10} \]

Системы (8), (9) и (10) являются дискретными линейными динамическими системами с постоянными матрицами, в которых значения функций $f(t)$ и $g(t)$ входят в двумерный вектор входных сигналов (воздействий).

Модель измерителя зададим в виде
\[ \begin{equation}
z_k =H c_k + \xi_k, \quad k=1,2,\ldots,K,
\end{equation} \tag{11} \]
где $H\in\mathbb{R}^{m\times n}$ — матрица измерений, определяющая структуру измерителя; $n$ — количество компонент вектора состояния $c_k$, а $m$ — количество его измеряемых компонент; $\xi_k\in\mathbb{R}^m$ — погрешность измерителя.

3. Этап I: Идентификация неизвестного вектора параметров $\boldsymbol \theta$ дискретной линейной стохастической модели по известным входным и выходным данным

Обозначим через $\theta$ неизвестный (в общем случае векторный) параметр линейной динамической системы (8), (9) или (10), который необходимо идентифицировать по доступным измерениям (11).

Дискретные модели конвекции–диффузии–реакции (8), (9) или (10) с моделью измерителя (11) можно представить в общем виде:
\[ \begin{equation}
\left\{
\begin{array}{l}
c_k = F(\theta) c_{k-1} + B(\theta)u_{k-1},\\
z_k =H c_k + \xi_k,\quad k = 1,2,\ldots,K,
\end{array}
\right.
\end{equation} \tag{12} \]
где $\theta=(v,\alpha,\beta,\lambda)^{\top}\in\mathbb{R}^{4}$ — неизвестный векторный параметр, $c_k\in\mathbb{R}^n$ — вектор состояния, $u_k\in\mathbb{R}^r$ — вектор входных сигналов (воздействий), $z_k\in\mathbb{R}^m$ — вектор измерений. Предположим, что шаги пространственно-временной сетки $\Delta x$ и $\Delta t$ заданы, в уравнении измерений аддитивный шум $\xi_k\in\mathbb{R}^m$ является последовательностью нормально распределенных независимых случайных векторов с нулевым математическим ожиданием и известной положительно определенной ковариационной матрицей $R\in\mathbb{R}^{m \times m}$. Матрицы $F(\theta)\in\mathbb{R}^{n \times n}$, $B(\theta)\in\mathbb{R}^{n \times r}$, входящие в уравнение состояния модели, зависят от параметра $\theta$.

Рассмотрим задачу параметрической идентификации модели (12) по доступным измерениям $z_k$ с целью оценки неизвестного (векторного) параметра $\theta$.

На этапе I предполагаем, что функции $f(t)$ и $g(t)$, входящие в граничные условия (3), (4) или (5), известны. Последнее означает, что в соответствующей дискретной модели (12) вектор $u_k$ представляет собой известный тестовый сигнал. Требуется вычислить оценку $\hat\theta^*$ неизвестного параметра $\theta$ по известным входным сигналам $U^{K-1}_0=\{u_0,u_1,\dots,u_{K-1}\}$ и выходным данным измерений $Z^{K}_1=\{z_1,\dots,z_{K}\}$ в соответствии с выбранным критерием качества идентификации $\mathcal{J}(\theta;Z^{K}_1,U^{K-1}_0)$.

Указанная задача параметрической идентификации сводится к решению задачи нелинейного программирования
\[ \begin{equation}
\hat{\theta}_{\min}=\operatorname*{argmin}\limits_{\theta\in D(\theta)}\mathcal{J}(\theta;Z^{K}_1,U^{K-1}_0),
\end{equation} \tag{13} \]
где $D(\theta)\subseteq\mathbb{R}^{4}$ — область определения параметра $\theta$. Таким образом, оценку $\hat\theta^*$ определим как $\hat\theta^*=\hat{\theta}_{\min}$.

Область $D(\theta)$ допустимых значений параметра $\theta$ определяется из условий сходимости его оценок при решении задачи (13). Для дискретных линейных стохастических систем вида (12) общие условия сходимости известны и определены в [12, Леммы 2.1 и 3.1]. В данной работе предполагаем, что область $D(\theta)$ является компактом в $\mathbb{R}^4$ и удовлетворяет указанным условиям сходимости для рассматриваемого класса моделей: $\forall \theta\in D(\theta)$ дискретная модель (12) является полностью управляемой и наблюдаемой, а элементы матриц $F(\theta)$ и $B(\theta)$ — непрерывно дифференцируемые по $\theta$ функции.
Например, свойство полной наблюдаемости обеспечивается выполнением следующего условия [31]:
\[ \begin{equation}
\operatorname{rank} \mathcal{M}_{DTI}(\theta)=n,
\end{equation} \tag{14} \]
где
\[ \begin{equation}
\mathcal{M}_{DTI}(\theta)=\left[
\begin{array}{ccccc}
H^{\top} & (HF(\theta))^{\top} & (HF^2(\theta))^{\top} & \cdots & (HF^{n-1}(\theta))^{\top}
\end{array}
\right]^{\top}
\end{equation} \tag{15} \]
есть матрица наблюдаемости модели (12). Условие (14) накладывает ограничения на выбор структуры измерителя, определяемой матрицей $H$.

Задача выбора оптимальной структуры измерителя с минимальным количеством сенсоров может быть решена путем анализа свойства полной наблюдаемости линейной динамической системы с учетом (14). В частности, в [27] показано, что для идентификации параметров моделей конвективно-диффузионного переноса достаточно всего лишь двух сенсоров. Из (15) следует, что очевидным является выбор $H$ в форме единичной матрицы, что означает наличие $n$ сенсоров для сбора данных измерений.

Для нахождения решения задачи (13) можно использовать известные численные методы оптимизации: градиентный метод, метод Ньютона, метаэвристические методы (например, генетический алгоритм или метод имитации отжига) и др. [15]. Готовые программные реализации данных методов, как правило, требуют от пользователя задания начального значения параметра $\theta$, ограничений на переменные и описания целевой функции. Вопрос выбора конкретного численного метода оптимизации зависит от специфики решаемой задачи.

В качестве целевой функции для решения задачи численной минимизации (13) выберем отрицательную логарифмическую функцию правдоподобия [7]
\[ \begin{equation}
\mathcal{J}_{LR}(\theta;Z^{K}_1,U^{K-1}_0)=\frac{Km}{2} \ln(2\pi)+\frac{1}{2}\sum_{k=1}^{K}
\bigl\{\ln|\Sigma_{\nu,k}(\theta)|+\nu^{\top}_k(\theta)\Sigma^{-1}_{\nu,k}(\theta) \nu_k(\theta)
\bigr\},
\end{equation} \tag{16} \]
где вектор невязки измерений $\nu_k(\theta)=z_{k}-H\hat c_{k}(\theta)$ и его ковариационную матрицу $\Sigma_{\nu,k}(\theta)=E \bigl\{\nu_k(\theta) \nu^{\top}_k(\theta)\bigr\}$ при заданных значениях параметра $\theta$ вычисляют по известным уравнениям фильтра Калмана [32].

Следует отметить, что ранее решение задачи параметрической идентификации моделей конвективно-диффузионного переноса с известными граничными условиями с помощью метода максимального правдоподобия на основе стандартного фильтра Калмана получено в [26]. Решение задачи параметрической идентификации скорости конвекции в модели конвективно-диффузионного переноса с помощью метаэвристических алгоритмов рассмотрено в [24] также с применением стандартного алгоритма Калмана. В [27] предложен новый подход к идентификации параметров дискретных моделей конвективно-диффузионного переноса на основе численно устойчивого SVD-фильтра Калмана. Предложенные ранее решения могут быть использованы для реализации первого этапа комбинированного метода идентификации моделей конвекции–диффузии–реакции.

4. Этап II: Численная идентификация неизвестных граничных условий

На втором этапе оценку параметра дискретной модели (12) считаем известной, т. е. ${\theta=\hat\theta^*}$, а вектор $u_k$, соответствующий граничным условиям (3),
(4) или (5), неизвестным.

Рассмотрим следующую модель дискретной стохастической системы
\[ \begin{equation}
\left\{
\begin{array}{l}
c_k = F(\hat\theta^*) c_{k-1} + B(\hat\theta^*)u_{k-1},\\
z_k =H c_k + \xi_k,\quad k = 1,2,\ldots,K,
\end{array}
\right.
\end{equation} \tag{17} \]
где $\hat\theta^*$ — оценка параметра модели, вычисленная на этапе I; $c_k\in\mathbb{R}^n$ — вектор состояния; $u_k\in\mathbb{R}^r$ — вектор неизвестных входных воздействий; $z_k\in\mathbb{R}^m$ — вектор доступных измерений. Предположим, что погрешность измерителя ${\xi_k\in\mathbb{R}^m}$, как и на первом этапе, имеет характеристики $\xi_k\sim \mathcal{N}(0,R)$.

Рассмотрим задачу идентификации вектора $u_k$ модели (17) по доступным измерениям $Z^K_1$ с одновременным оцениванием вектора состояния $c_k$. Поскольку мы предполагаем, что априорная информация об эволюции входного сигнала отсутствует, для решения задачи идентификации будем использовать подход, разработанный С. Гиллейнсом и Б. Де-Мором [16].

Алгоритм Гиллейнса – Де-Мора состоит из трех последовательных шагов, повторяемых в цикле:

1. обновление оценки вектора состояния по времени;
2. оценка вектора неизвестного входного воздействия;
3. обновление оценки вектора состояния по текущим измерениям.

В работе [16] рассмотрены два варианта алгоритма, в которых шаги 1 и 2 совпадают, а шаги 3 различаются. В первом варианте на шаге 2 получается MVU-оценка (MVU — minimum-variance unbiased) вектора $\hat u_{k-1}$, а на шаге 3 — несмещенная оценка вектора состояния $\hat c_k$. Во втором варианте алгоритма на шаге 3 за счет более сложных вычислений получается MVU-оценка вектора состояния. Алгоритмы требуют выполнения условия
\[ \begin{equation}
\mathop{\mathrm{rank}} HB(\hat\theta^*) = \mathop{\mathrm{rank}} B(\hat\theta^*) = r,
\end{equation} \tag{18} \]
где $r$ — размер вектора $u_k$.

Условие (18) является достаточным для существования несмещенной оценки вектора состояния в алгоритмах одновременного оценивания вектора состояния и неизвестных входных сигналов [16, 18].

Отметим, что ранее решение задачи идентификации неизвестных граничных условий в модели конвективно-диффузионного переноса на основе алгоритмов Гиллейнса – Де-Мора рассматривалось в [25, 29]. В [28] получена квадратно-корневая модификация алгоритма Гиллейнса – Де-Мора. Предложенные ранее решения могут быть использованы для реализации второго этапа комбинированного метода идентификации моделей конвекции–диффузии–реакции.

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

Рассмотрим идентификацию параметров модели вида (1), (2), (4):
\[ \begin{equation*}
\dfrac{\partial c(x,t)}{\partial t}+v\dfrac{\partial c(x,t)}{\partial x} = \alpha\dfrac{\partial^{2} c(x,t)}{\partial x^{2}}-\beta c(x,t), \quad x\in[0{;}1],\quad t\in[0{;}1];
\end{equation*} \]
\[ \begin{equation*}
c(x,0) = \varphi(x);
\end{equation*} \]
\[ \begin{equation*}
c(0,t) = f(t), \quad \dfrac{\partial c(1,t)}{\partial x} = -\lambda[c(1,t)-g(t)],
\end{equation*} \]
где $v=3$, $\alpha=0.8$, $\beta=4$, $\lambda=0.5$, $\varphi(x)\equiv 0$.

Процесс идентификации будем моделировать в системе MATLAB. Зададим в рассматриваемой области пространственно-временную сетку c 11 узлами по оси $Ox$ и 201 узлом по оси $Ot$ ($\Delta x = 0.1$, $\Delta t = 0.005$) и получим решение прямой задачи методом конечных разностей. Смоделируем зашумленные измерения в узлах пространственной сетки $x=0.1, 0.2, \ldots, 1$, соответствующих всем компонентам вектора состояния $c^k_1$, $c^k_2$, $\ldots$, $c^k_{10}$. Матрица измерений имеет вид $H=I_{10}$.

Пусть функции, входящие в граничные условия, известны: $f(t)=g(t)\equiv5$. На рис. 2 приведены графики решения прямой задачи и смоделированных измерений с матрицей ковариации шума $R=0.05^2I_{10}$.

Рис. 2. Графики решения задачи (a) и зашумленных измерений (b)
[Figure 2. Plots of the solution (a) and noisy measurements (b)]

Минимизация критерия идентификации (16) выполнялась при помощи функции fmincon. Проведенные эксперименты показали, что точность идентификации параметров зависит как от вида дискретизации пространственно-временной области и уровня шума, так и от выбора граничных условий, начального решения, настроек функции минимизации и других факторов. В таблице приведены результаты идентификации параметров для различных значений уровня шума $R$ по результатам 100 экспериментов ($\mu$ — среднее значение, $\sigma^2$ — дисперсия, RMSE — среднеквадратичная ошибка, MAPE — средняя абсолютная процентная ошибка). Поиск значений каждого параметра осуществлялся на отрезке $[0.01; 5]$. В качестве начального решения выбиралась середина отрезка. Основные параметры функции fmincon: SpecifyObjectiveGradient = false, Algorithm = sqp, MaxFunctionEvaluations = 500.

Результаты идентификации [Identification results]

$R$	Parameter	$\mu$	$\sigma^2$	RMSE	MAPE
$0.05^{2}I_{10}$	$v$	2.990654	0.018655	0.136219	3.527277
	$\alpha$	0.794360	0.000096	0.011265	1.131970
	$\beta$	2.628127	5.509623	2.708609	56.280028
	$\lambda$	0.500398	0.000097	0.009822	1.589091
$0.01^{2}I_{10}$	$v$	2.995260	0.001055	0.032662	0.863970
	$\alpha$	0.798047	0.000032	0.005990	0.570405
	$\beta$	3.516434	2.039151	1.500865	28.702054
	$\lambda$	0.500095	0.000004	0.001965	0.318911
$0.005^{2}I_{10}$	$v$	2.998785	0.000287	0.016914	0.445395
	$\alpha$	0.799324	0.000010	0.003165	0.313788
	$\beta$	3.834971	0.611877	0.795608	15.794188
	$\lambda$	0.500047	0.000001	0.000982	0.159357

Definition: $R$ is the noise covariance matrix; $\mu $ is the mean of data set; $\sigma$ is the standard deviation; RMSE is the root mean square error; MAPE is the mean absolute percentage error

Результаты экспериментов показывают, что с уменьшением уровня шума в измерителе точность идентификации параметров возрастает. При этом параметр $\beta$ идентифицируется хуже всех, что может быть объяснено малым вкладом слагаемого $r_2$ в коэффициент $a_2$ дискретной модели.

Используем идентифицированные значения параметров $v$, $\alpha$, $\beta$, $\lambda$ (в таблице выделены жирным) для оценивания значений неизвестных функций $f(t)$ и $g(t)$, входящих в граничные условия.

Пусть $f(t)=4\left|{3t-\lfloor3t+0.5\rfloor}\right|$ (треугольная волна с периодом 1/3 и амплитудой 2), $g(t)=t$.

Зададим в рассматриваемой области пространственно-временную сетку c 9 узлами по оси $Ox$ и 201 узлом по оси $Ot$ ($\Delta x = 0.125$, $\Delta t = 0.005$) и получим решение прямой задачи методом конечных разностей. Смоделируем зашумленные измерения в узлах пространственной сетки $x=0.125$ и ${x=1}$, соответствующих компонентам вектора состояния $c^k_1$ и $c^k_{8}$. В этом случае матрица измерений имеет вид
\[ \begin{equation*}
H=\begin{bmatrix}
1&0&0&0&0&0&0&0\\
0&0&0&0&0&0&0&1\\
\end{bmatrix}.
\end{equation*} \]

На рис. 3 приведены графики решения прямой задачи и смоделированных измерений с матрицей ковариации шума $R=0.01^2I_{2}$.

Рис. 3. Графики решения задачи (a) и зашумленных измерений (b)
[Figure 3. Plots of the solution (a) and noisy measurements (b)]

Рис. 4. Оценки $f(t)$ (a), $g(t)$ (b) и решения задачи (c)
[Figure 4. Estimates of $f(t)$ (a), $g(t)$ (b) and solution (c)]

Для одновременного оценивания значений решения, а также функций $f(t)$ и $g(t)$ используем алгоритм Гиллейнса – Де-Мора, результаты работы которого представлены на рис. 4. Из приведенных графиков видно, что значения функции $f(t)$ оцениваются точнее, чем значения функции $g(t)$ ($\mathrm{RMSE}_f =0.034908$, $\mathrm{RMSE}_g=0.175363$).

Для сравнения заменим в рассматриваемой модели левое граничное условие на условие третьего рода
\[ \begin{equation*}
\frac{\partial c(0,t)}{\partial x} = \lambda[c(0,t)-f(t)]
\end{equation*} \]
и снова воспользуемся алгоритмом Гиллейнса – Де-Мора. Соответствующие графики приведены на рис. 5, 6. В данном случае ошибки оценивания значений $f(t)$ и $g(t)$ имеют одинаковый порядок ($\mathrm{RMSE}_f=0.184948$, $\mathrm{RMSE}_g =0.175259$).

Рис. 5. Графики решения задачи (a) и зашумленных измерений (b)
[Figure 5. Plots of the solution (a) and noisy measurements (b)]

Рис. 6. Оценки $f(t)$ (a), $g(t)$ (b) и решения задачи (c)
[Figure 6. Estimates of $f(t)$ (a), $g(t)$ (b) and solution (c)]

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

В работе предложен новый комбинированный двухэтапный метод параметрической идентификации моделей тепломассопереноса, на примере идентификации математических моделей конвекции–диффузии–реакции с неизвестными граничными условиями при наличии случайных помех в измерениях значений искомой функции $c(x,t)$. Решение заключается в переходе от исходной модели, описываемой уравнениями в частных производных, к дискретной линейной стохастической модели в пространстве состояний с неизвестными входными сигналами и разработке комбинированного двухэтапного метода идентификации, основанного на последовательном применении метода минимизации критерия идентификации безградиентного типа и рекуррентного метода одновременного оценивания вектора состояния и неизвестных входных сигналов.

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

Основными результатами работы являются:

1. новые дискретные линейные стохастические модели конвекции–диффузии–реакции для разных типов граничных условий;
2. общая схема процесса параметрической идентификации, включающая поэтапную идентификацию неизвестных параметров математической модели и идентификацию неизвестных граничных условий.

С целью проверки работоспособности предложенного подхода построена компьютерная модель конвекции–диффузии–реакции и выполнена реализация всех алгоритмов на языке MATLAB. Идентификация параметров проводилась безградиентным численным методом минимизации критерия идентификации, а идентификация граничных условий — рекуррентным алгоритмом Гиллейнса – Де-Мора. Проведена серия вычислительных экспериментов, результаты которых подтверждают работоспособность предложенного решения.

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

Конкурирующие интересы. Конкурирующих интересов не имеем.
Авторский вклад и ответственность. Ю.В. Цыганова — общая концепция статьи; схема комбинированного двухэтапного метода идентификации; работа с черновиком и переработанным вариантом рукописи. А.В. Цыганов — программная реализация алгоритмов; подготовка первичного варианта рукописи. А.Н. Кувшинова — вычислительные эксперименты; подготовка первичного варианта рукописи. Д.В. Галушкина — дискретные модели конвекции–реакции–диффузии; подготовка первичного варианта рукописи. Авторы несут полную ответственность за предоставление окончательной рукописи в печать. Окончательная версия рукописи была одобрена всеми авторами.
Финансирование. Исследование выполнено за счет гранта Российского научного фонда № 23–21–00361, https://rscf.ru/project/23-21-00361/.

Об авторах

Юлия Владимировна Цыганова

Ульяновский государственный педагогический университет им. И. Н. Ульянова

Автор, ответственный за переписку.
Email: tsyganovajv@gmail.com
ORCID iD: 0000-0001-8812-6035
SPIN-код: 8259-4594
Scopus Author ID: 6507218923
ResearcherId: F-7169-2013
http://www.mathnet.ru/person69680

доктор физико-математических наук, доцент; профессор; каф. информационных технологий

Россия, 432071, Ульяновск, пл. Ленина, 4/5

Андрей Владимирович Цыганов

Ульяновский государственный университет

Email: andrew.tsyganov@gmail.com
ORCID iD: 0000-0002-4173-5199
SPIN-код: 2729-7659
Scopus Author ID: 35186570100
ResearcherId: C-2331-2014
http://www.mathnet.ru/person178940

кандидат физико-математических наук, доцент; профессор; каф. высшей математики

Россия, 432017, Ульяновск, ул. Льва Толстого, 42

Анастасия Николаевна Кувшинова

Ульяновский государственный университет

Email: kuvanulspu@yandex.ru
ORCID iD: 0000-0002-3496-5981
SPIN-код: 2849-0643
Scopus Author ID: 57204965949
http://www.mathnet.ru/person141068

кандидат физико-математических наук; доцент; каф. высшей математики

Россия, 432017, Ульяновск, ул. Льва Толстого, 42

Дарья Валерьевна Галушкина

Ульяновский государственный педагогический университет им. И. Н. Ульянова

Email: dgalushkina73@gmail.com
ORCID iD: 0000-0003-4041-0533
https://www.mathnet.ru/person178370

аспирант; каф. информационных технологий

Россия, 432071, Ульяновск, пл. Ленина, 4/5

Список литературы

Дополнительные файлы