Идентификация параметров процесса аномальной диффузии на основе разностных уравнений Текст научной статьи по специальности «Математика»

Вычислительные технологии

Том 18, № 1, 2013

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

А. С. ОвсиЕнко Самарский государственный технический университет, Россия e-mail: [email protected]

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

Ключевые слова: аномальная диффузия, параметрическая идентификация, линейно-параметрическая дискретная модель, разностное уравнение.

Введение

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

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

При определении параметров линейной динамической системы достаточно широко используются методы идентификации на основе результатов измерений мгновенных значений ук реакции системы на некоторое типовое тестовое воздействие, например, путём применения моделей авторегрессии или скользящего среднего [3]. Обычно эти модели не учитывают случайный характер возмущений ек в результатах наблюдений ук = Ук + ек и построенные на их основе численные алгоритмы не используют статистические методы обработки результатов эксперимента и обладают низкой помехозащищённостью. Устранить указанный недостаток позволяет применение линейно-параметрических дискретных моделей (ЛПДМ) [4, 5]. Особенностью данных моделей

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

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

1. Постановка задачи

Математическая модель процесса аномальной диффузии может быть выражена при помощи дробно-дифференциального уравнения с дробными производными по пространству и времени:

д7с(х,Ь) 1дас(х,Ь)

—^-= «—--, (1)

дП дха ' v у

здесь Ь > 0, 0 <7 < 1, 1 < а < 2, с(х,Ь) — объёмная концентрация диффундирующих

частиц, d(x) — коэффициент диффузии. Дробная производная произвольного порядка

в имеет вид

X

те Л(х) = ^в! (х) =_1_^ Г /(е) ^

1 )(х) dx[в] ^ 1} (х) г ([в] +1 - в) dxв+1 ] (х- е)в-[в],

0

где [в] — целая часть числа в. Начальные условия для уравнения (1) будут следующими:

с(х, 0) = со(х). (2)

Краевые условия по пространственной координате за пределами исследуемого интервала принимаются равными нулю.

Численное решение уравнения аномальной диффузии было получено [6, 7] при помощи формулы Грюнвальда — Летникова. Явная схема численного решения дробно-дифференциального уравнения аномальной диффузии имеет вид

1 к+1 л ¿+1 1 ^к—1+1_ ^

-- У^ 91,1°к + = уа За>пСг—п+1,

ла

1=0 п=0

где коэффициенты д7>1 и да,п определяются как

Г(в - в)

дв,

Г(-в)Г(з + 1)-

Порядок аппроксимации точного решения составляет О (к + т). Схема (3) условно т7 dmax 7

устойчива при ——— < —. Шаблон численной схемы представлен на рис. 1. Схема (3)

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

Рис. 1. Шаблон численной схемы (3)

2. Численный метод идентификации

Введём обозначения: N — объём выборки по времени, N — размерность пространственной выборки. При помощи элементарных преобразований выражение (3) можно представить в виде:

т ^ л 'у А 7 ^ 1 ^ 1

с?+1 = С+1 + ^ + АТ- (-а) Е М^-з + ^ + 7 £ М/с?-1-3. (5)

где

3=0 3=0

1 — V 2 — V 3 — V I + 1 — V М = — ■ — ■ — ■ ...■ -++-Л-, " = 7,а, I.» > 0.

В случае I = 0 или к = 0 соответствующие суммы в (5) принимают значения, равные нулю. Выражение (5) представляет собой рекуррентную зависимость значения концентрации в точке с(%Н, (к + 1)т) от её значений во всех "предыдущих" точках (см. рис. 1).

Рекуррентные соотношения (5) позволяют свести поставленную задачу идентификации процесса аномальной диффузии к построению линейно-параметрической дискретной модели и последующему оцениванию её коэффициентов, которые линейным образом связаны с подлежащими определению параметрами задачи типа задачи Ко-ши (1), (2). Данный алгоритм среднеквадратичного оценивания позволяет учитывать в результатах измерений и стохастическую составляющую, что обеспечивает помехозащищённость метода.

На основе разностной схемы (5) может быть построена линейно-параметрическая дискретная модель, которая в матричной форме имеет вид

{:=2,+* (6)

где Ь = (с0,с1,..,с^ж-1, с0,... с^_1,... с^-1)Т — ^Ж^-мерный вектор правой части; ^ — матрица регрессоров (ЖЖх) х (3Жж); е — ^Ж^-мерный вектор случайной аддитивной помехи в результатах наблюдений; Р — матрица размера ЖгЖх х ЖЖх в стохастическом уравнении эквивалентного возмущения; п = (П1, П2,... П^^)Т — Ж^х-мерный вектор эквивалентного возмущения в стохастическом разностном уравнении. Вектор коэффициентов ЛПДМ (6) Л размера 3Жх выражается через параметры решения дробно-диф-

ференциального уравнения (1) следующим образом:

\ \ ат7 ^i-1 -10 лт Ai = 7, Ai+i =--——, г = 1

Т 7 di_ 1

A1+Nx+i = , г = 1, 2,..,Nx - 1

Ai+2Nx+i = co(ih), г = 0,1, 2,..,Nx - 1. (7) Матрица регрессоров F имеет вид

Ук+М = c0-i, г = 1, 2,..,Nx,

1 — Y

f2Nx+i,1 = C0-i + c1-l, г = 1, 2,

j-4

f(j-i)Nx+i,i = ^ MMj-3-kck-i, г = 1, 2,..,Nx, j = 4, 5,..,Nt, fc=0

?- 2

/ü-1)Nx+1,2 = c0 , j = 2 3, ^ "2

i-3

/(j-i)Nx+i,i+i = ^ MM0_2-fcck-2, г = 3, 4,.., Nx, j = 2, 3,.., Nt, fc=0

/(j-1)Nx+2,3 = ^-c0 2 + Ci 2, j = 2, 3, -Л

/ü-1)Nx+i,1+Nx+i = ^ 2> г = 1, j = 2, ^»Л

/i,2Nx+i = 1, г = 1,2,..,Nx, j = 2,3,..,Nt. (8)

Элементы матрицы размера NtNx эквивалентного возмущения P в стохастическом

уравнении (6) выражаются в виде 1 г = j;

-Awx+fc+i, г = (m + 1)Nx + k, j = mNx + k + 1,

k = 1, 2,..,Nx - 1, m = 0,1,..,Nt - 2;

— A1 — Ak+1, г = (m + 1)Nx + k, j = mNx + k, k = 1, 2,.., Nx, m = 0,1,..,Nt - 2; PiJ = ^ -AiMsY-1, г = (m + 1 + s)Nx + k, j = mNx + k, (9)

k = 1, 2,..,Nx, m = 0,1,..,Nt - 2 - s, s = 1, 2,..,Nt - 2; -Afc+i+sMk-1, г = (m + 1)Nx + k + s, j = mNx + s,

k = 1, 2,.., Nx - 1, m = 0,1,..,Nt - 2, s = 1, 2,..,Nx - k; 0, иначе.

При использовании среднеквадратичного критерия [3] на конечном множестве точек {xi, tj}, г = 0,1,.., Nx - 1, j = 0,1,.., Nt - 1, где NtNx — объём выборки, минимизируется функционал

IH|2 = l|c - с||2 ^ min, (10)

или в развёрнутой форме

Nx-1 Nt-1 Nx-1 Nt-1

J] J] (cj - С)2 = ^ (£j+(i-1)Nt)2 ^ min, i=0 j=0 i=0 j=0

где е = cj — cj — вектор остатков, cj — дискретные значения функции концентрации, полученные экспериментально, cj — искомые оценки значений функции концентрации. Оценки коэффициентов ЛПДМ находятся из условия минимизации среднеквадратичного отклонения (10) построенной модели пространственно-временной функциональной зависимости cj от экспериментальных данных cj. С этой целью первое уравнение в (6) преобразуется к виду

Р-1Ъ = P-1FA + е . (11)

Таким образом, задача сводится к минимизации функционала

—> min . (12)

I -1|2

е

л 2

Р-1Ъ — P-1FA

Л Л

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

Оценки коэффициентов линейно-параметрической дискретной модели на начальном этапе реализации итерационной процедуры [3] вычисляются при помощи метода наименьших квадратов по формуле

Л(0) = (рт Р )-1РТ Ь, (13)

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

А(в) = (РтП-1-!)р)РтП-ЦЬ, в = 1, 2,.., ЭМ». (14)

Здесь Л — вектор оценок коэффициентов стохастического разностного уравнения, вычисленный на в-м шаге, П—1_1) = (Р-(в_1))-1 — квадратная симметричная матрица размера х ^¿Жк, сходящаяся к решению Л матричного уравнения Рт(Р-РТ)-1 РЛ =

(Р-РТ)-1Ь. По найденным оценкам Л^, ] = 1, 2,.., 3Жх, могут быть вычислены параметры функции (3).

Однако в выражении (2) содержатся параметры 7 и а, которые считаются неизвестными, в связи с чем описанный алгоритм идентификации требует существенного преобразования. Автором разработана итерационная процедура, позволяющая реализовать предложенный метод идентификации следующим образом: на первом шаге значения параметров 7 и а в формуле (2) принимаются равными некоторым начальным приближениям 7(0) и а(0), на следующем итерационном шаге значения параметров пере-считываются через оценки коэффициентов ЛПДМ, вычисленные на каждой итерации по формулам (13), (14), и с учётом (2) определяются следующим образом:

а« = - , т(г) = Л« (15)

После нескольких итераций оценки параметров 7 и а принимаются равными значениям 7(г) и а(г), вычисленным на последней итерации, а оценки остальных параметров

процесса (1) могут быть найдены по формулам

Л2+Мх+г

Лг

т 7

1 — -

0,1,..,Мх - 2,

На Л

1 + Мх

т «7

Ло(гк) — Л1+2Мх+г, г — 0,1,..,^ - 1.

'16)

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

3. Результаты численного моделирования

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

Д7, %

30 25 20 15 10 5 0

■ г

\

?

I

ч

Да, % 20 15 10 5 0

б

■ 1

» « «

£

» к

0 2 4 6 8 №

0 2 4 6 8 №

8,%-20 15

Ю

•Л:

п___

0 2

6 8 №

а

в

Рис. 2. Зависимость смещения оценки параметров 7 (а), а (б) и среднеквадратичного отклонения в (в) от числа итераций N

На рис. 2, а, б показаны результаты исследований сходимости итерационной процедуры для различных значений начальных приближений 7(0) и а(0) в зависимости от числа итераций. Исследования проводились на тестовом сигнале при 7 = 0.3, а = 1.3, с0 = 0.6, й = 0.1, г = 0,..,Ях, N = 100, N = 3, к = 0.1, т = 0.01, е = 0. Данные усреднялись по пяти реализациям. Тремя линиями на рисунках показаны результаты, полученные при различных значениях начальных приближений 7(0) и а(0). Видно, что

8,% 1.2 1.0 0.8 0.6 0.4 0.2 0

Д7,%

0.6 0.5 0.4 0.3 0.2 0.1 0

1.5 1.0 0.5

♦

4 ♦ ► +

4 Л ♦

< ♦— >*

0 *

0

4 <

♦ < ♦

4 ♦ 4>-

0 1 2 3 4 £, %

3 4 е,%

0 1 2 3 4 е, % Да, % 0.25 0.20 0.15 0.10 0.05 0

♦ ♦

♦

♦ < ♦

а1,%

0.9

0.6

0.3

0 1 2 3 4 е,%

д

0 *

0

4>

4 ♦ ♦ ♦ ►

4 ♦

4 ♦ И <►- ►

3 4 е, %

Рис. 3. Зависимости среднеквадратичного отклонения в (а), смещения оценок параметров 7 (б), а (в) и дисперсии оценок параметров 7 (г), а (д) от величины случайной аддитивной помехи в результатах наблюдений е

а

в

г

с ростом количества итераций величины смещения в оценках параметров Л и а, а также среднеквадратичного отклонения модели в (рис. 2, в) от входной выборки данных стремятся к нулю, и при любых начальных приближениях 7(0) Е (0,1) и а(0) Е (1, 2) итерационная процедура сходится к точному решению уже при числе итераций N — 6.

На рис. 3 представлены графики зависимостей среднеквадратичного отклонения в построенной на основе результатов идентификации модели от точного решения, а также смещения А и дисперсии а2 оценок параметров дробных дифференциальных операторов, входящих в уравнение (1). Исследования проводились при начальных условиях, описанных выше. Результаты усреднялись по 100 реализациям. В качестве начальных приближений были выбраны значения а(0) — 1.5 и 7(0) — 0.5 — середины соответствующих интервалов. Число итераций принималось равным N — 6. Случайная аддитивная помеха е в результатах наблюдений считалась удовлетворяющей нормальному распределению. Из приведенных графиков видно, что при величине случайной аддитивной помехи е — 5 % величина среднеквадратичного отклонения модели от экспериментальных данных и дисперсии порядков дробных дифференциальных операторов не превышают 2%, а смещения полученных оценок находятся в пределах 1 %. Аналогичная картина наблюдается и в случае оценок коэффициентов диффузии, где дисперсия не превышает 5 %, а смещение составляет менее 2 %. Полученные данные свидетельствуют о том, что погрешность оценивания предложенным методом не превышает величины случайной аддитивной помехи в результатах наблюдений, что позволяет говорить о высокой эффективности описанной процедуры.

Были проведены также численные исследования смещения и дисперсии оценок в зависимости от периодов дискретизации т и к, показавшие, что представленный метод имеет высокую скорость сходимости на интервалах т Е (0.1, 0.25) и к Е (0.09, 0.17). Применение метода в указанных интервалах периодов дискретизации позволяет минимизировать значения погрешности оценивания.

Выводы

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

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

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

[1] ЭйкхОФФ П., Ванечек А., Савараги Е. и др. Современные методы идентификации систем / Под ред. П. Эйкхоффа. М.: Мир, 1983. 400 с.

[2] ОвсиЕнко А.С. Метод параметрической идентификации процесса аномальной диффузии на границе застойной зоны и области радиального стока скважины // Тез. докл. V науч-но-практ. конф. по математическому моделированию и компьютерным технологиям в процессах разработки месторождений. М.: Нефтяное хозяйство, 2012. C. 37.

[3] Зотеев В.Е. Параметрическая идентификация диссипативных механических систем на основе разностных уравнений. М.: Машиностроение, 2009. 344 с.

[4] Зотеев В.Е. Разработка и исследование линейных дискретных моделей колебаний диссипативных систем // Вестник СамГТУ. Физико-математические науки. 1999. Вып. 7. С. 170-177.

[5] Радченко В.П., Зотеев В.Е. Определение динамических характеристик механической системы на основе стохастических разностных уравнений колебаний // Изв. вузов. Машиностроение. 2007. № 1. С. 3-10.

[6] Петухов А.А., РЕвизников Д.Л., Сластушенский Ю.В. Описание аномальной диффузии с использованием дробно-дифференциальных уравнений и метода дискретных элементов// Материалы XVII Междунар. конф. по вычислительной механике и современным прикладным программным системам. М.: МАИ-ПРИНТ, 2011. С. 601-602.

[7] Петухов А.А., РЕвизников Д.Л. Алгоритмы численного решения дробно-дифференциальных уравнений // Вестник МАИ. 2009. Т. 16, № 6. С. 228-234.

Поступила в 'редакцию 30 июля 2012 г., с доработки — 26 ноября 2012 г.