Some aspects of the implementation of the PRPHMM 1.0 software package for refining the parameters of hereditary mathematical models of radon transfer in a storage chamber

Full Text

Введение

В последние годы возрос интерес к исследованию так называемых дифференциальных уравнений дробного порядка [1, 2], в которых неизвестная функция содержится под знаком производной дробного порядка. Это обусловлено как развитием самой теории дробного интегрирования и дифференцирования, так и приложениями таких конструкций в различных областях науки [3, 4].

Рассматривается дробное уравнение вида:

${\partial }_{0t}^{\alpha \mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}}x\mathrm{<mo\; stretchy="false">(</mo>}\phi \mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}F\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>,}t\mathrm{<mo\; stretchy="false">)</mo>,}x\mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo><mo>=</mo>}{x}_0\mathrm{,}$ (1)

где, дробная производная типа Герасимова-Капуто переменного 0 < α(t) < 1 порядка имеет вид:

${\partial }_{0t}^{\alpha \mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}}x\mathrm{<mo\; stretchy="false">(</mo>}\phi \mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}\frac{1}{\Gamma \mathrm{<mo\; stretchy="false">(</mo>1}-\alpha \mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo>}}{\displaystyle {\int }_0^t}\frac{dx\mathrm{<mo\; stretchy="false">(</mo>}\phi \mathrm{<mo\; stretchy="false">)</mo>}}{dt}\frac{1}{{\mathrm{<mo\; stretchy="false">(</mo>}t-\phi \mathrm{<mo\; stretchy="false">)</mo>}}^{\alpha \mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}}}d\phi \mathrm{,}$ (2)

где, $\Gamma \mathrm{<mo\; stretchy="false">(</mo>}\cdot \mathrm{<mo\; stretchy="false">)</mo>}$ $–$ известная гамма-функция Эйлера.

Уравнение (1) лежит в основе некоторых разработанных авторами математических моделей динамических процессов: солнечной активности [5], динамики смертности от COVID-19 [6] а также моделировании объемной активности радона в накопительной камере [7]. Это достигается за счёт введения эффекта памяти [8], с помощью (2), в модель того или иного процесса, где порядок производной будет связан с интенсивностью процесса. В тоже время, выбирая в модели (1) вид функции F(A(t), t) описывающие различные механизмы исследуемого динамического процесса, можно моделировать различные динамические режимы.

Модельные уравнения на основе (1) решаются с помощью численных методов. В работах [9, 10] подробно исследуется математический аппарат модельных уравнений типа (1).

Отметим, что в упомянутых исследованиях параметры моделей подбирались вручную по максимуму коэффициента детерминации и коэффициента корреляции Пирсона с экспериментальными данными объемной активности радона (ОАР), учитывая имеющееся понимание природы моделируемого процесса. Такой подход является трудоемким, что приводит к идеям применения способов автоматизации подбора оптимальных параметров за счет решения соответсвующей обратной задачи.

В настоящей статье мы будем рассматривать случай, когда α(t) в уравнении (1) является константой (эредитарная α-модель). В статье авторов [11] были приведены результаты применения программного комплекса PRPHMM 1.0 для восставновления значений проницаемости среды и коэффициента воздухообмена на пункте радонового мониторинга Карымшина. В настоящей статье дается описание программного комплекса PRPHMM 1.0, приводится листинг модулей и пример применения.

Эредитарная α-модель ОАР

Для решения проблемы поиска оптимальных параметров можно решать обратную задачу $—$ распространенный тип задач во многих научных областях, где необходимо определить значения параметров модели на основе наблюдаемых данных [12], но невозможно провести прямые измерения этих параметров. Необходимость такого подхода часто возникает при работе с геологическими данными, в геофизике и сейсмологии [13] и др.

Далее, для иллюстрации кода и того что он реализует, рассмотрим конкретный пример а именно эредитарную α-модель ОАР в накопительной камере. Подробнее о процессе преноса радона можно узнать из работ [14, 15], а о α-модели [7, 11].

Анализ данных ОАР и сопутствующих параметров получаемых в ходе непрерывного мониторинга, является одним из методов поиска предвестников землетрясений. Это связано с тем, что на ОАР влияют изменения напряженно-деформированного состояния среды через которую подпочвенный газ выходит на поверхность [16], поэтому радон считается известным и хорошо себя зарекомендовавшим индикатором процессов протекающий в такой среде. Мониторинг ²²²Rn как метод поиска предвестников сейсмических событий, за последние годы прекрасно себя показал, особенно как краткосрочный предвестник (до 15 суток) [17].

Эредитарная $\alpha$ -модель ОАР представляется следующей задачей Коши:

${\partial }_{\mathrm{0,}t}^{\alpha }A\mathrm{<mo\; stretchy="false">(</mo>}\phi \mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}-{\lambda }_{0}A\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}+{\lambda }_0\mathrm{,}A\mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo><mo>=</mo>}{A}_0\mathrm{,}$ (3)

где, $A\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}$ $–$ функция решения, зависимость ОАР от времени в камере; ${\partial }_{\mathrm{0,}t}^{\alpha }A\mathrm{<mo\; stretchy="false">(</mo>}\phi \mathrm{<mo\; stretchy="false">)</mo>}$ $–$ дробная производная Герасимова-Капуто [18, 19] постоянного порядка 0 < α < 1, частный случай (2); α $–$ интенсивность переноса радона, порядок дробной производной; λ₀ $–$ коэффициент воздухообмена (КВО) в камере.

Задача (3) решается численно в равномерной сеточной области:

$\begin{array}{c}h\mathrm{<mo>=</mo>}T\mathrm{<mo>/</mo>}N\mathrm{,}\widehat{\Omega }\mathrm{<mo>=</mo>}\left\{\mathrm{<mo\; stretchy="false">(</mo>}{t}_i\mathrm{<mo>=</mo>}ih\mathrm{<mo\; stretchy="false">)</mo><mo>:</mo>0}\le i\mathrm{<mo><</mo>}N\right\}\mathrm{,}\widehat{\mathbb{A}}\in \widehat{\Omega }\mathrm{,}\\ A\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{A}_i\mathrm{,}\mathrm{0<mo><</mo>}{A}_i\mathrm{<mo><</mo>1.}\end{array}$ (4)

а для решения (3) на сетке (4) использоваться безусловно устойчивая численная схема IFDS апробированная в ряде тестовых и прикладных задач [5, 10].

Методика решения обратной задачи для $\alpha$ -эредитарной модели

Пусть $A_i\in \widehat{\mathbb{A}}$ $–$ функция некоего известного класса сеточных функций, но её решение зависит от набором параметров $\overleftarrow{X}\mathrm{<mo>=</mo>}\left[X_0\mathrm{,...,}{X}_{K-1}\right]$, где $K\mathrm{<mo>=</mo>2}$ $–$ число восстанавливаемых параметров, а X₀ = α, X₁ = λ₀. Пусть значения дискретной функции решения $A_i\in \widehat{\mathbb{A}}$ неизвестны, но известна дополнительная информация (экспериментальные данные RVA) $A_i\mathrm{<mo>=</mo>}{\theta }_i\mathrm{<mo>=</mo>}\overleftarrow{\theta }$ о решении разностной прямой задачи Коши для эредитарной $\alpha$ -модели RVA.

Тогда разностная обратная задача для (3) на сетке (4) $–$ это восстановление значений $\overleftarrow{X}\mathrm{<mo>=</mo>}\left[X_0\mathrm{,}{X}_1\right]$ по известным A_i = θ_i экспериментальным данным:

$\begin{array}{c}{A}_i\mathrm{<mo>=</mo>1}-\frac{\widehat{{\partial }_{\mathrm{0,}ih}^{X_0}}{A}_i}{X_1}\mathrm{,}{A}_i\mathrm{<mo>=</mo>}{\theta }_i\mathrm{,}1\le i\mathrm{<mo><</mo>}N\mathrm{,}\\ \widehat{{\partial }_{\mathrm{0,}ih}^{X_0}}{A}_i\mathrm{<mo>=</mo>}\frac{h^{-X_0}}{\Gamma \mathrm{<mo\; stretchy="false">(</mo>2}-X_0\mathrm{<mo\; stretchy="false">)</mo>}}{\displaystyle \sum _{j\mathrm{<mo>=</mo>0}}^{i-1}}\left({\mathrm{<mo\; stretchy="false">(</mo>}j+\mathrm{1<mo\; stretchy="false">)</mo>}}^{1-X_0}-j^{1-X_0}\right)\left(A_{i-j}-A_{i-j-1}\right)\mathrm{.}\end{array}$ (5)

Для решения (5) обратимся к теории безусловной оптимизации [20]. Для этого необходимо минимизировать функционал невязки:

$\begin{array}{c}\overleftarrow{\eta }\mathrm{<mo>=</mo>}\overleftarrow{\theta }-\omega \mathrm{<mo\; stretchy="false">(</mo>}\overleftarrow{X}\mathrm{<mo\; stretchy="false">)</mo>,}\min \left(\Psi \left(\overleftarrow{X}\right)\right)\mathrm{<mo>=</mo>}\frac{1}{2}{\displaystyle \sum _{i\mathrm{<mo>=</mo>0}}^{N-1}}{\eta }_i^2\mathrm{<mo>=</mo>}\frac{1}{2}{\displaystyle \sum _{i\mathrm{<mo>=</mo>0}}^{N-1}}{\left({\theta }_i-{\omega }_i\right)}^2\mathrm{,}\end{array}$ (6)

где, $\overleftarrow{\eta }$ $–$ вектор невязки размерности N > K, а вектор $\omega \mathrm{<mo\; stretchy="false">(</mo>}\overleftarrow{X}\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}\left[{\omega }_0\mathrm{,...,}{\omega }_N\right]$ $–$ вектор модельных данных, т. е. решение прямой задачи относительно некоторого приближения $\overleftarrow{X}$, получаемого в ходе решения обратной задачи.

Разностная обратная задача решается методом безусловной оптимизации ньютоновского типа [21], а именно итерационным методом Левенберга-Марквардта [22, 23], представимого в виде:

$\Delta X\mathrm{<mo>=</mo>}\left(-H^{-1}\right)\times \left(J^T\times \overleftarrow{\eta }\right)\mathrm{,}H\mathrm{<mo>=</mo>}{J}^T\times J+\gamma E\mathrm{,}$ (7)

где ΔX $–$ оптимальное приращение $\overleftarrow{X}$ для следующей итерации; E $–$ единичная матрица размерности K×K; $J\mathrm{<mo>=</mo>}J\left(\overleftarrow{X}\right)$ $–$ матрица Якоби размерности N×K с элементами вычисляемыми во формуле: $J_{i\mathrm{,}k}\mathrm{<mo>=</mo>}\frac{\partial {\eta }_i}{\partial X_k}\mathrm{,}i\mathrm{<mo>=</mo>0..}N-\mathrm{1,}k\mathrm{<mo>=</mo>0..}K-1$; производная $\frac{\partial {\eta }_i}{\partial X_k}$ аппроксимируется разностным оператором $J_{i\mathrm{,}k}\mathrm{<mo>=</mo>}\frac{{\eta }_i^{\delta }-{\eta }_i}{\delta X_k}$, где δX $–$ заданное малое приращение $\overleftarrow{X}$; γ $–$ параметр регуляризации метода. Если γ∈ℝ>0, а также матрица Гессе H положительно определена, то тогда ΔX является направлением спуска для оптимального шага метода; Стартовое значение: ${\gamma }^{\mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo>}}\mathrm{<mo>=</mo>}v\cdot {{\displaystyle \max }}_i\left(diag\left(J{\left(X^{\mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo>}}\right)}^T\times J\left(X^{\mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo>}}\right)\right)\right)$, где v $–$ заданная стартовая константа.

Решение обратной задачи (5) методом Левенберга-Марквардта (7), далее (IP-LB), сводится к тому, чтобы в ходе цикла, начиная с заданных постоянных X⁽⁰⁾, δX, v а также c $–$ константы для пересчёта γ, многократно вычисляя решение прямой задачи при приближениях $\overleftarrow{X}$, получаемых в ходе решения обратной задачи, вычислить оптимальные значения $\overleftarrow{X}$. Подробнее с алгоритмом для реализующий оптимизацию вектора $\overleftarrow{X}$ можно ознакомиться в работе [24].

Критерием того, что метод сходится к оптимальному решению является ε ≤ Σ, где Σ $–$ заданная точность решения IP-LB, $\epsilon \mathrm{<mo>=</mo>}\frac{1}{N}{\displaystyle {\sum }_{i\mathrm{<mo>=</mo>0}}^{N-1}}{\left[{\eta }_i^{\Delta }\right]}^2$ $–$ среднеквадратичная ошибка между экспериментальными и модельными данными RVA.

Критерий того, что метод попал в "ловушку" локального минимума $–$ это отсутствие существенного изменения значений $\overleftarrow{\Delta X}$ в ходе итераций. Иначе говоря $\Delta \left(\overleftarrow{\Delta X}\right)\to 0$, где $\Delta \left(\overleftarrow{\Delta X}\right)$ $–$ оценка скорости «приращения приращений» определяемая:

$\Delta \left(\overleftarrow{\Delta X}\right)\mathrm{<mo>=</mo>}\underset{n\to \infty }{{\displaystyle \lim }}\left(\frac{1}{K}{\displaystyle \sum _{k\mathrm{<mo>=</mo>0}}^{K-1}}\left|\Delta X_k^{\mathrm{<mo\; stretchy="false">(</mo>}n\mathrm{<mo\; stretchy="false">)</mo>}}-\Delta X_k^{\mathrm{<mo\; stretchy="false">(</mo>}n-\mathrm{1<mo\; stretchy="false">)</mo>}}\right|\right)\to 0.$

Программная реализация алгоритма решения обратной задачи

Далее в (лист. 1$–$25) представлен код реализующий описанный итерационный метод Левенберга-Марквардта, для решения обратной задачи в рамках программного комплекса PRPHMM 1.0 (рис. 1) на языке MATLAB [25]. Решение основного матричного уравнения (7) представлено в коде (лист. 14$–$15).

 

Рис. 1. Файловая структура программного комплекса PRPHMM 1.0, в частности ветвь содержащая подпрограммы для решения обратной задачи

[Figure 1. File structure of the PRPHMM 1.0 software package, in particular the branch containing subroutines for solving the inverse problem ]

 

Листинг 1. Calc_Coeff_and_Add_param_CUSTOM.m

1 % пере-вычисление значений коэф. на основе описанных пользов. формул

2 function [ model ] = ...

3 Calc_Coeff_and_Add_param_CUSTOM( model, IM )

4 if( isfield( model( IM ), ’derivative’ ) == 1 )

5 if( isfield( model( IM ).derivative,’const_of_process’) == 0 )

6 disp( [’ параметр [model(’ num2str(IM) ...

7 ’).derivative.const_of_process] НЕ определён ’ ...

8 ’ -->> поэтому, будет создано = 1 ’] );

9 model( IM ).derivative.const_of_process = 1;

10 end

11 end

12 if( isfield( model( IM ), ’parameter’ ) == 1 )

13 names_fields = fieldnames( model( IM ).parameter );

14 num_fields_coeff = length( names_fields );

15 for sub_ind = 1 : num_fields_coeff

16 name_curr_field = char(names_fields(sub_ind));

17 formula = ’ ’;

18 recive_formula = sprintf(...

19 ’formula = model(%d).parameter.%s.define;’,...

20 IM, string(name_curr_field) );

21 eval( recive_formula );

22 model(IM).parameter.(name_curr_field).define_start_type=class(formula);

23 replace_define_field = sprintf(...

24 ’model(%d).parameter.%s.define = ’’%s’’;’,...

25 IM, string(name_curr_field), string(formula) );

26 eval( replace_define_field );

27 t = model( IM ).t;

28 h = model( IM ).h;

29 N = model( IM ).N;

30 T = model( IM ).T;

31 ui = model( IM ).ui;

32 if (isa(formula, ’char’) == 1)

33 array_values = eval(formula);

34 else

35 array_values(1: N + 1) = formula;

36 end

37 eval( sprintf( ...

38 ’model(%d).parameter.%s.values = array_values;’,...

39 IM, string( name_curr_field )) );

40 end

41 end

42 end

 

Листинг 2. set__parameter_equation.m

1 % пере-расчет параметров модели для решения обратной задачи

2 function [ model_UPDATE ] = set__parameter_equation( model, IM, X_n,X_n_m1)

3 model_UPDATE = model;

4 model_UPDATE( IM ).parameter.alpha.define = X_n( 1 );

5 model_UPDATE( IM ).lambda = X_n( 2 );

6 [ model_UPDATE ] = Calc_Coeff_and_Add_param_CUSTOM( model_UPDATE, IM );

7 disp( [ ’| X_0 = ’ num2str( X_n_m1(1) ) ’ -> | ’ num2str(X_n(1))]);

8 disp( [ ’| X_1 = ’ num2str( X_n_m1(2) ) ’ -> | ’ num2str(X_n(2))]);

9 end

 

Листинг 3. Vector__bias_value.m

1 % вычисление невязки--значения

2 function [ bias ] = Vector__bias_value( X, divisor )

3 bias = 0.0;

4 X_size = size(X,1);

5 for i = 1 : X_size

6 bias = bias + X( i )^2;

7 end

8 bias = bias / divisor;

9 end

 

Листинг 4. get__appoximation_X_n_delta.m

1 % присвоение заданого малого приращения

2 function [ X_n_delta ] = get__appoximation_X_n_delta( X_n, delta_X )

3 X_n_delta = X_n + delta_X;

4 end

 

Листинг 5. get__appoximation_X_n__Delta.m

1 % присвоение ОПТИМАЛЬНОГО (вычисленного) приращения

2 function [ X_n__Delta ] = get__appoximation_X_n__Delta( X_n, Delta__X )

3 X_n__Delta = X_n + Delta__X;

4 end

 

Листинг 6. get__delta_X.m

1 % СТАРОТОВОЕ - ПО умолчанию, приращение для параметров функции

2 function [ delta_X ] = get__delta_X( model, IM )

3 delta_X( 1 ) = model( IM ).type_task.par_rest.X_0_inc_start;

4 delta_X( 2 ) = model( IM ).type_task.par_rest.X_1_inc_start;

5 end

 

Листинг 7. get__Delta__X.m

1 % получение ОПТИМАЛЬНОГО (вычисленного) приращения

2 function [ Delta__X ] = get__Delta__X( result_M_eq )

3 n_d = size( result_M_eq, 1 );

4 for i = 1 : n_d

5 Delta__X( i ) = result_M_eq( (n_d + 1) - i );

6 end

7 end

 

Листинг 8. Get_Experemental_Data.m

1 % получение вектора эксперементальных данных

2 function [expr_data] = Get_Experemental_Data( model, IM, data_set )

3 ind_ED = model( IM ).type_task.INDEX_exper_data;

4 expr_data = data_set( ind_ED ).DATA_.VALUE;

5 end

 

Листинг 9. get__X_0.m

1 % стартовая итерация, для 1-го решения прямой задачи

2 function [ X_0 ] = get__X_0( model, IM )

3 X_0( 1 ) = model( IM ).type_task.par_rest.X_0_start;

4 X_0( 2 ) = model( IM ).type_task.par_rest.X_1_start;

5 end

 

Листинг 10. get__X_default.m

1 % получить и запомнить значения X_0 (ПО умолчанию)

2 function [ X_default ] = get__X_default( model, IM )

3 X_default( 1 ) = model( IM ).parameter.alpha.values( 1 );

4 X_default( 2 ) = model( IM ).lambda;

5 end

 

Листинг 11. Matrix_Jacobi.m

1 % вычисление матрицы Якоби

2 function [ J ] = Matrix_Jacobi( X_n_p1, X_n, delta_X, type_finite_diff )

3 size__X_n_p1 = size( X_n_p1 , 1 );

4 size__X_n = size( X_n , 1 );

5 size__delta_X = size( delta_X, 2 );

6 if( size(X_n_p1,1) ~= size(X_n,1) )

7 locat = sprintf( ’Error ! dimensions wrong, %d != %d’, ...

8 size__X_n_p1, size__X_n );

9 error( char(locat) );

10 end

11 if( type_finite_diff == "upward f-d" )

12 for i = 1 : size__X_n_p1

13 for j = 1 : size__delta_X

14 J(i,j) = ( X_n_p1(i) - X_n(i) ) / delta_X(j);

15 end

16 end

17 end

18 end

 

Листинг 12. Matrix_Jacobi_custom.m

1 % вычисление матрицы Якоби

2 function [J_out] = Matrix_Jacobi_custom( J, X, delta_X,type_finite_diff,...

3 type_calc_Jacobi )

4 J_out = J;

5 X_size = size(X,2);

6 if( type_finite_diff == "upward f-d" )

7 if( type_calc_Jacobi == "classicaly" )

8 for i = 1 : X_size

9 for j = 1 : X_size

10 J_out(i,j) = ( ( X(i) + delta_X(j) ) - X(i) ) / delta_X(j);

11 end

12 end

13 else

14 error( [’НЕИЗВЕСТНЫЙ : [type_calc_Jacobi]’ ] );

15 end

16 else

17 error( [’НЕИЗВЕСТНЫЙ : [type_finite_diff]’ ] );

18 end

19 end

 

Листинг 13. calc__J_i_j.m

1 % вычисление основной матрицы Якоби

2 function [ J_i_j ] = calc__J_i_j( eta_n_delta, eta_n, delta_X, n )

3 locat_nd = sprintf( ’eta^( X^(%d) + delta X )’, n );

4 locat = sprintf( ’eta^( X^(%d) )’, n );

5 [J_i_j] = Matrix_Jacobi( eta_n_delta, eta_n, delta_X, ’upward f-d’ );

6 display_control_sum( locat_nd, eta_n_delta );

7 display_control_sum( locat, eta_n );

8 display_control_sum( ’J_i_j ’, J_i_j );

9 end

 

Листинг 14. calc__Gesse_Matrix.m

1 % вычисоение матрицы Гессе - ключевой для метода реш. обратной задачи

2 function [ Gesse_Matrix ] = calc__Gesse_Matrix( J_i_j, gamma_reg )

3 size_J_i_j = size(J_i_j,2);

4 [ gamma_E ] = Unitary_Matrix( size_J_i_j, size_J_i_j );

5 gamma_PART = gamma_E * gamma_reg;

6 Gesse_Matrix = (J_i_j’ * J_i_j) + gamma_PART;

7 end

 

Листинг 15. calc__Matrix_equation.m

1 % вычисление решения основного матричного уравнения

2 function [ result_M_eq ] = calc__Matrix_equation(n, J_i_j,gamma_reg, Theta)

3 [ Gesse_Matrix ] = calc__Gesse_Matrix( J_i_j, gamma_reg );

4 Gesse_Matrix_Inverse = inv( Gesse_Matrix );

5 J_i_j_T_mul_Theta = - 1 * J_i_j’ * Theta;

6 result_M_eq = Gesse_Matrix_Inverse * J_i_j_T_mul_Theta;

7 end

 

Листинг 16. Solve_Inverse_Problem.m

1 % решение обратной задачи

2 function [ model_OUT ] = Solve_Inverse_Problem( model, IM, expr_data )

3 model_OUT = model;

4 tic;

5 switch( model( IM ).type_task.method )

6 case( "Levenberg-Marquardt" )

7 [ model_OUT ] = IP_Levenberg_Marquardt( model, IM, expr_data );

8 otherwise

9 error( [’НЕИЗВЕСТНЫЙ МЕТОД : [model(’ num2str( IM ) ...

10 ’).only_IP.method] -- решения Обратной задачи’ ] );

11 end

12 model_OUT( IM ).time_calc = toc;

13 model_OUT( IM ).RES_.var.x_raw = model_OUT( IM ).RES_.var.x_raw’;

14 model_OUT(IM).RES_.diff.x_dot_raw = model_OUT(IM).RES_.diff.x_dot_raw’;

15 model_OUT( IM ).RES_.var.x = model_OUT( IM ).RES_.var.x’;

16 model_OUT( IM ).RES_.diff.x_dot = model_OUT( IM ).RES_.diff.x_dot’;

17 end

 

Листинг 17. RP_LM_seq___step_1.m

1 function [ gamma_reg ] = RP_LM_seq___step_1( k, n_d, n,count,v,X_0,delta_X)

2 display__step( 1, n, count );

3 J_X_0 = zeros( k, n_d );

4 [ J_X_0 ]=Matrix_Jacobi_custom(J_X_0,X_0,delta_X,"upward f-d","classicaly");

5 J_X_0_T = J_X_0’;

6 B_0 = J_X_0_T * J_X_0;

7 md_B_0( 1 : size(B_0,2) ) = 0;

8 for i = 1 : size(B_0,2)

9 md_B_0( i ) = B_0( i,i );

10 end

11 max_md_B_0 = max( md_B_0 );

12 gamma_reg = v * max_md_B_0;

13 end

 

Листинг 18. RP_LM_seq___step_2.m

1 function [ s_0, eta_0 ] = RP_LM_seq___step_2( model, IM, n, count, X_0, ...

2 X_default, expr_data, s_0_in)

3 model_UPDATE = model;

4 display__step( 2, n, count );

5 [ model_UPDATE ] = ...

6 set__parameter_equation( model, IM, X_0, X_default );

7 [ model_UPDATE ] = Solve_Direct_Problem( model_UPDATE, IM );

8 local_res = model_UPDATE( IM ).RES_.var.x;

9 [ eta_0 ] = expr_data - local_res;

10 display_control_sum( ’eta_0’, eta_0 );

11 s_0_previous = s_0_in;

12 [ s_0 ] = Vector__bias_value( eta_0, 2 );

13 display__s_bias( n, s_0, s_0_previous, 0 );

14 end

 

Листинг 19. RP_LM_seq___step_3.m

1 function [ J_i_j, eta_n_delta ] = RP_LM_seq___step_3( model, IM,n,count,...

2 X_n, delta_X, expr_data, eta_n )

3 model_UPDATE = model;

4 display__step( 3, n, count );

5 [ X_n_delta ] = get__appoximation_X_n_delta( X_n, delta_X );

6 [ model_UPDATE ] = set__parameter_equation( model, IM, X_n_delta, X_n);

7 [ model_UPDATE ] = Solve_Direct_Problem( model_UPDATE, IM );

8 local_res = model_UPDATE( IM ).RES_.var.x;

9 [ eta_n_delta ] = expr_data - local_res;

10 [ J_i_j ] = calc__J_i_j( eta_n_delta, eta_n, delta_X, n );

11 end

 

Листинг 20. RP_LM_seq___step_4.m

1 function [ s_1, X_n__Delta, model_UPDATE ] = ...

2 RP_LM_seq___step_4( model, IM, n, count, J_i_j, ...

3 gamma_reg, Theta, X_n, expr_data, s_1_in )

4 model_UPDATE = model;

5 display__step( 4, n, count );

6 [ result_M_eq ] = calc__Matrix_equation( n, J_i_j, gamma_reg, Theta );

7 [ Delta__X ] = get__Delta__X( result_M_eq );

8 [ X_n__Delta ] = get__appoximation_X_n__Delta( X_n, Delta__X );

9 [ model_UPDATE ] = set__parameter_equation( model, IM, X_n__Delta,X_n);

10 [ model_UPDATE ] = Solve_Direct_Problem( model_UPDATE, IM );

11 local_res = model_UPDATE( IM ).RES_.var.x;

12 [ eta_n__Delta ] = expr_data - local_res;

13 display_control_sum( ’eta_n__Delta ’, eta_n__Delta );

14 s_1_previous = s_1_in;

15 [ s_1 ] = Vector__bias_value( eta_n__Delta, 2 );

16 display__s_bias( n, s_1, s_1_previous, 1 );

17 end

 

Листинг 21. RP_LM_seq___step_5__exit.m

1 function [ is_exit, X_n__SUMM, count_summ ] = ...

2 RP_LM_seq___step_5__exit( model, IM, n, count, ...

3 n_d, model_OUT, expr_data, ...

4 X_n__SUMM, X_n,count_summ, ...

5 count_summ_max )

6 display__step( 5, n, count );

7 is_exit = 0;

8 SIGMA = model( IM ).type_task.par_.SIGMA;

9 local_res = model_OUT( IM ).RES_.var.x;

10 if( try__exit__by_MSE( expr_data, SIGMA, n, local_res ) == 1 )

11 is_exit = 1;

12 return;

13 end

14 if( count_summ < count_summ_max )

15 X_n__SUMM = X_n__SUMM + X_n;

16 count_summ = count_summ + 1;

17 end

18 if( count_summ == count_summ_max )

19 if( try__freeze_in_local_optimum( X_n__SUMM, X_n, SIGMA, n,...

20 count_summ_max ) == 1 )

21 is_exit = 1;

22 return;

23 end

24 X_n__SUMM( 1 : n_d ) = 0;

25 count_summ = 0;

26 end

27 end

 

Листинг 22. RP_LM_seq___step_6__continue.m

1 function [ n, X_n, s_0, s_1, X_n__Delta, gamma_reg, model_OUT ] =...

2 RP_LM_seq___step_6__continue( model, IM, c, n, count, ...

3 s_0, s_1, gamma_reg, ...

4 model_OUT_last, expr_data, ...

5 X_n__Delta, delta_X )

6 model_OUT = model_OUT_last;

7 previous_gamma = -1;

8 display__step( 6, n, count );

9 if( s_0 > s_1 )

10 display__exit( n, ’s_0’, ’>’, ’s_1’, s_0, s_1,...

11 ’ n++, Пересчёт шагов 3 4 ’ );

12 s_0 = s_1;

13 X_n = X_n__Delta;

14 previous_gamma = gamma_reg;

15 gamma_reg = gamma_reg / c;

16 n = n + 1;

17 display__gamma_ext( n, ’gamma_reg’, ’/’, ’gamma_reg’, ’c’,...

18 gamma_reg, previous_gamma, c );

19 if( n >= 1 )

20 local_res = model_OUT_last( IM ).RES_.var.x;

21 eta_n = expr_data - local_res;

22 end

23 [ J_i_j, eta_n_delta ] = ...

24 RP_LM_seq___step_3( model, IM, n, count, ...

25 X_n, delta_X, expr_data, eta_n );

26 [ s_1, X_n__Delta, model_OUT ] = ...

27 RP_LM_seq___step_4( model, IM, n, count, J_i_j, ...

28 gamma_reg,eta_n_delta, X_n, ...

29 expr_data, s_1 );

30 else

31 previous_gamma = gamma_reg;

32 gamma_reg = gamma_reg * c;

33 display__gamma_ext( n, ’gamma_reg’, ’*’, ’gamma_reg’, ’c’,...

34 gamma_reg, previous_gamma, c );

35 display__exit( n, ’s_0’, ’<=’, ’s_1’, s_0, s_1,...

36 ’ Пересчёт 4 ’ );

37 [ s_1, X_n__Delta, model_OUT ] = ...

38 RP_LM_seq___step_4( model, IM, n, count, J_i_j, ...

39 gamma_reg,eta_n_delta, ...

40 X_n, expr_data, s_1 );

41 end

42 end

 

Листинг 23. IP_Levenberg_Marquardt.m

1 function [ model_OUT ] = IP_Levenberg_Marquardt( model, IM, expr_data )

2 model_OUT = model;

3 N = model( IM ).N;

4 n = 0;

5 count = 0;

6 COUNT__TIME_OUT_RP = model( IM ).type_task.par_.TIME_OUT;

7 k = model( IM ).type_task.par_rest.k;

8 n_d = k;

9 c = model( IM ).type_task.par_.c;

10 v = model( IM ).type_task.par_.v;

11 gamma_reg = 0;

12 s_0 = 0.0; s_1 = 0.0;

13 s_0_previous = 0.0; s_1_previous = 0.0;

14 % ======= обявление матриц и векторов необходимых для алгоритма =======

15 X_default = get__X_default( model, IM );

16 X_0( 1 : k ) = 0;

17 X_n( 1 : k ) = 0;

18 delta_X( 1 : n_d ) = 0;

19 J_i_j = zeros( N, n_d );

20 X_n__SUMM( 1 : n_d ) = 0;

21 % =================== СТАРТ решения обратной задачи ===================

22 [ X_0 ] = get__X_0( model, IM );

23 [ delta_X ] = get__delta_X( model, IM );

24 [ gamma_reg ] = RP_LM_seq___step_1( k, n_d, n, count, v, X_0, delta_X);

25 disp( [’gamma_reg = ’ num2str(gamma_reg)]);

26 [ s_0, eta_0 ] = RP_LM_seq___step_2( model, IM, n, count, X_0, ...

27 X_default, expr_data, s_0 );

28 if( n == 0 )

29 X_n = X_0;

30 eta_n = eta_0;

31 end

32 [ J_i_j, eta_n_delta ] = RP_LM_seq___step_3( model, IM, n,count, X_n,...

33 delta_X, expr_data, eta_n );

34 [ s_1, X_n__Delta, model_OUT ] =

35 RP_LM_seq___step_4( model, IM, n, count, J_i_j, gamma_reg, ...

36 eta_n_delta, X_n,expr_data, s_1 );

37 % ========================= ГЛАВНЫЙ ЦИКЛ ==============================

38 count_summ = 0;

39 count_summ_max = 5;

40 while( count < COUNT__TIME_OUT_RP )

41 [ is_exit, X_n__SUMM, count_summ ] = ...

42 RP_LM_seq___step_5__exit( model, IM, n, count, n_d, model_OUT,...

43 expr_data, X_n__SUMM, X_n, count_summ, count_summ_max );

44 if( is_exit == 1 )

45 break;

46 end

47 [ n, X_n, s_0, s_1, X_n__Delta, gamma_reg, model_OUT ] = ...

48 RP_LM_seq___step_6__continue( model, IM, c, n, count, s_0, s_1, ...

49 gamma_reg, model_OUT, expr_data, X_n__Delta, delta_X );

50 count = count + 1;

51 end

52 end

 

Листинг 24. try__exit__by_MSE.m

1 % проверка : выход по достижению заданого среднего квадратического ?

2 function [bool] = try__exit__by_MSE( expr_data,SIGMA, n,result__X_n__Delta)

3 bool = 0;

4 MSE = - 1.0;

5 M = size( expr_data, 1 );

6 eta_mse = expr_data - result__X_n__Delta;

7 MSE = Vector__bias_value( eta_mse, M );

8 if( MSE <= SIGMA )

9 bool = 1;

10 display__exit( n, ’MSE’, ’<=’, ’SIGMA’, MSE, SIGMA, ’ EXIT ’ );

11 else

12 display__exit( n,’MSE’,’>’,’SIGMA’,MSE,SIGMA,’ ВЫЧИСЛЯЕМ ДАЛЬШЕ ’);

13 end

14 end

 

Листинг 25. try__freeze_in_local_optimum.m

1 % проверка : застревания решения возле некоторого состояния

2 function [ bool ] = try__freeze_in_local_optimum( X_n__SUMM, X_n__last, ...

3 SIGMA,n, count_summ_max )

4 bool = 0;

5 X_n_size = size(X_n__SUMM,2); SIGMA_pow_2 = SIGMA^2;

6 X_n__SUMM_mean( 1 : X_n_size ) = 0;

7 X_n__SUMM_mean = X_n__SUMM * ( 1 / count_summ_max );

8 X_n__SUMM_mean_bias( 1 : X_n_size ) = 0;

9 X_n__SUMM_mean_bias = X_n__SUMM_mean - X_n__last;

10 disp( [ ’ try__freeze_in_local_optimum : X_n__SUMM_mean_bias = ’...

11 num2str(X_n__SUMM_mean_bias) ] );

12 X_n_MSE = Vector__bias_value( X_n__SUMM_mean_bias, X_n_size );

13 disp( [ ’ try__freeze_in_local_optimum : MSE ’ ...

14 ’( X_n__SUMM_mean, X_n__last ) = ’ num2str(X_n_MSE) ] );

15 if( X_n_MSE <= SIGMA_pow_2 )

16 bool = 1;

17 display__exit( n, ’MSE( X_n )’, ’<=’, ’SIGMA ^2’, ...

18 X_n_MSE, SIGMA_pow_2, ’ ВЫХОД ( Застыл !!! ) ’);

19 else

20 display__exit( n, ’MSE( X_n )’, ’>’, ...

21 ’SIGMA ^2 (проверка застревания в локальном оптимуме)’, ...

22 X_n_MSE, SIGMA_pow_2, ’ ПЕРЕПРОВЕРКА в следующем ОКНЕ ’ );

23 end

24 end

 

Тестовый пример

Листинг 26. Файл управляющих параметров – mod_1

1 RESET;

2 COLOR;

3 DEF_COUNTERS;

4 TVHN;

5 % -----------------< Определение ... МОДЕЛЕЙ >-----------------------

6 model = [];

7 % -------------------------------------------------------------------------

8 % classic KRMR-2_(24_aug-27_aug)_(h_1_6)

9 IND = IND + 1;

10 model(IND).name = "Классическая модель RVA";

11 model(IND).N = 10000;

12 model(IND).h = 1/6;

13 model(IND).start_point = 0.9;

14 model(IND).A_max = 1;

15 model(IND).lambda = 0.06;

16 model(IND).parameter.a.define = 0;

17 model(IND).parameter.b.define = - model(IND).lambda;

18 model(IND).parameter.c.define = model(IND).A_max * model(IND).lambda;

19 model(IND).derivative.type = ’VOGC’;

20 model(IND).parameter.alpha.define = 0.9999;

21 % -- подстройка парам. [model().] под данные

22 model(IND).impact_data.INDEX_data_set = 1;

23 model(IND).impact_data.re_define_N = 1;

24 model(IND).impact_data.re_define_start_point = 1;

25 % -- численный метод решения

26 model(IND).num_method.result_last_cutt = 1;

27 model(IND).num_method.type = ’IFDS (MNM)’;

28 model(IND).num_method.IFDS_accuracy = 1e-3;

29 model(IND).num_method.IFDS_start_iter = ’EFDS u(N)’;

30 % -- тип решения Direct problem

31 model(IND).type_task.type = ’DIRECT problem’;

32 % -- отрисовка

33 model(IND).display.is_plot_proc = 1;

34 model(IND).display.parameter.style = ’--’;

35 model(IND).display.parameter.color = COLOR_.Blue;

36 model(IND).display.result.style = ’-’;

37 model(IND).display.result.color = COLOR_.Blue;

38 % -------------------------------------------------------------------------

39 % KRMR-2_(24_aug-27_aug)_(h_1_6)

40 IND = IND + 1;

41 model(IND).name = "Эредитарная \alpha-модель RVA";

42 model(IND).N = 10000;

43 model(IND).h = 1/6;

44 model(IND).start_point = 0.9;

45 model(IND).A_max = 1;

46 model(IND).lambda = 0.055;

47 model(IND).parameter.a.define = 0;

48 model(IND).parameter.b.define = - model(IND).lambda;

49 model(IND).parameter.c.define = model(IND).A_max * model(IND).lambda;

50 model(IND).derivative.type = ’VOGC’;

51 model(IND).parameter.alpha.define = 0.818;

52 % -- подстройка парам. [model().] под данные

53 model(IND).impact_data.INDEX_data_set = 1;

54 model(IND).impact_data.re_define_N = 1;

55 model(IND).impact_data.re_define_start_point = 1;

56 % -- численный метод решения

57 model(IND).num_method.result_last_cutt = 1;

58 model(IND).num_method.type = ’IFDS (MNM)’;

59 model(IND).num_method.IFDS_accuracy = 1e-3;

60 model(IND).num_method.IFDS_start_iter = ’EFDS u(N)’;

61 % -- тип решения Direct problem

62 model(IND).type_task.type = ’DIRECT problem’;

63 % -- отрисовка

64 model(IND).display.is_plot_proc = 1;

65 model(IND).display.parameter.style = ’--’;

66 model(IND).display.parameter.color = COLOR_.Bright_green;

67 model(IND).display.result.style = ’-’;

68 model(IND).display.result.color = COLOR_.Bright_green;

69 % -------------------------------------------------------------------------

70 % INVERSE KRMR-2_(24_aug-27_aug)_(h_1_6)

71 IND = IND + 1;

72 model(IND).name = "Обратная задача";

73 model(IND).N = 10000;

74 model(IND).h = 1/6;

75 model(IND).start_point = 0.9;

76 model(IND).A_max = 1;

77 model(IND).lambda = 0.06;

78 model(IND).parameter.a.define = ’0 * ui’;

79 model(IND).parameter.b.define = ’- model(IM).lambda * ui’;

80 model(IND).parameter.c.define = ’model(IM).A_max * model(IM).lambda * ui’;

81 model(IND).derivative.type = ’VOGC’;

82 model(IND).parameter.alpha.define = ’0.9999 * ui’;

83 % -- подстройка парам. [model().] под данные

84 model(IND).impact_data.INDEX_data_set = 1;

85 model(IND).impact_data.re_define_N = 1;

86 model(IND).impact_data.re_define_start_point = 1;

87 % -- численный метод решения

88 model(IND).num_method.result_last_cutt = 1;

89 model(IND).num_method.type = ’IFDS (MNM)’;

90 model(IND).num_method.IFDS_accuracy = 1e-3;

91 model(IND).num_method.IFDS_start_iter = ’EFDS u(N)’;

92 % -- тип решения Inverse problem

93 model(IND).type_task.type = ’INVERSE problem’;

94 model(IND).type_task.method = ’Levenberg-Marquardt’;

95 model(IND).type_task.INDEX_exper_data = 1;

96 model(IND).type_task.par_.TIME_OUT = 50;

97 model(IND).type_task.par_.v = 100;

98 model(IND).type_task.par_.c = 2;

99 model(IND).type_task.par_.SIGMA = 7e-3;

100 model(IND).type_task.par_rest.k = 2;

101 model(IND).type_task.par_rest.X_0_start = 0.05;

102 model(IND).type_task.par_rest.X_1_start = 0.0025;

103 model(IND).type_task.par_rest.X_0_inc_start = 0.01;

104 model(IND).type_task.par_rest.X_1_inc_start = 0.0005;

105 % -- отрисовка

106 model(IND).display.is_plot_proc = 1;

107 model(IND).display.parameter.style = ’--’;

108 model(IND).display.parameter.color = COLOR_.Red;

109 model(IND).display.result.style = ’-’;

110 model(IND).display.result.color = COLOR_.Red;

111 % --------------< Определение ... НАБОРОВ ДАННЫХ >--------------------

112 data_set = [];

113 IND = 0;

114 % -------------------------------------------------------------------------

115 % KRMR-2_(24_aug-27_aug)_(h_1_6)

116 IND = IND + 1;

117 data_set(IND).file.name = ’KRMR-2.xls’;

118 data_set(IND).file.name_add = ...

119 "Эксперементальные данные RVA (KRMR-2: 24 aug-27 aug) (h 1/6)";

120 data_set(IND).file.name_column_Data = "middle (Radon1 (Bq.m3))";

121 data_set(IND).file.index_column_date = [1];

122 % -- выборка по календарю -

123 data_set(IND).date_bound_op.index_type_datatime = 4;

124 data_set(IND).date_bound_op.filler_Second = 0;

125 data_set(IND).date_bound.start.Year = 2020;

126 data_set(IND).date_bound.start.Month = 08;

127 data_set(IND).date_bound.start.Day = 24;

128 data_set(IND).date_bound.start.Hour = 13;

129 data_set(IND).date_bound.start.Minute = 40;

130 data_set(IND).date_bound.end.Year = 2020;

131 data_set(IND).date_bound.end.Month = 08;

132 data_set(IND).date_bound.end.Day = 27;

133 data_set(IND).date_bound.end.Hour = 03;

134 data_set(IND).date_bound.end.Minute = 30;

135 % -- обработка

136 data_set(IND).processing.normalize_on_max = true;

137 % -- отрисовка

138 data_set(IND).display.data.style = ’-’;

139 data_set(IND).display.data.color = COLOR_.Black;

140 % ---------------< ЗАПУСК ПРОГРАММНОГО КОМПЛЕКСА >--------------------

141 main_PRPHMM_1_0;

 

На рис. 2 приведены графики ОАР, полученные по результам модельных расчетов (классическая модель и α-модель) в программном комплексе PRPHMM 1.0 и результатам наблюдений с пункта Камчатского филиала федерального исследовательского центра «Единая геофизическая служба РАН» Карымшина (KRMR-2) 24-27 августа 2024 г. (предоставлены к.ф.-м.н., Макаровым Е.О.). Мы видим, что решение обратной задачи позволяет давать в автоматическом режиме приемлемые оценки порядка дробной производной α и коэффициента воздухообмена λ₀ в модельном уравнении. Коэффициент воздухообмена достаточно тяжело измерить с помощью геофизической аппаратуры. В тоже время зная его можно расчитать плотность потока радона в накопительной камере, что является важным при исследований аномалий в ОАР, в том числе предшествующих сейсмическим событиям.

 

Рис. 2. Результат работы программного комплекса PRPHMM 1.0. Визуализация экспериментальных и модельных данных [11].

[Figure 2 The result of the PRPHMM 1.0 software package. Visualisation of experimental and model data [11].]

 

Заключение

Программный комплекс PRPHMM 1.0 позволяет проводить автоматизацию расчетов для уточнения значений параметров эредитарных математических моделей переноса радона в накопительной камере. В статье с помощью программного комплекса PRPHMM 1.0 и экспериментальных данных ОАР на пункте Карымшина для эредитарной α-модели переноса радона в накопительной камере были уточнены параметры порядка дробной производной α (проницаемость среды) и коэффициента воздухообмена λ₀. Расчеты показали, что значения оценимаемых параметров является адекватными. Программный комплекс PRPHMM 1.0 успешно прошел тестирование и в настоящее время ведется работа по восставновлению значний параметров функций α(t) для эредитарной α(t)-модели.

Благодарности. Авторы благодарят д.ф.-м.н. Паровика Р.И. за ценные замечания при обсуждении результатов, полученных в статье.

About the authors

Dmitrii A. Tverdyi

Institute of Cosmophysical Research and Radio Wave Propagation FEB RAS

Author for correspondence.
Email: tverdyi@ikir.ru
ORCID iD: 0000-0001-6983-5258

PhD (Phys. & Math.), Researcher, Electromagnetic Radiation Laboratory

Russian Federation, 684034, Paratunka village, Mirnaya str., 7

Evgeny O. Makarov

Federal Research Center «Unified Geophysical Service of the Russian Academy of Sciences»

Email: tverdyi@ikir.ru
ORCID iD: 0000-0002-0462-3657

PhD (Phys. & Math.), Senior Researcher, Acoustic and Radon Monitoring Laboratory, Kamchatka Branch

Russian Federation, 683023, Petropavlovsk-Kamchatsky, Piipa Boulevard st., 9

References

Supplementary files