Development of identification methods for fractional differential equations with riemann-liouville fractional derivative

Full Text

Введение. До настоящего времени актуальной остаётся задача разработки эффективных методов параметрической идентификации систем, описываемых при помощи дифференциальных уравнений c дробными производными. Несмотря на обилие литературы по дробному исчислению (см., например, библиографию в [1-11]) и множество приложений в задачах реологии, вязкоупругости, аномальной диффузии в пористых (фрактальных) структурах, теории автоматического управления, физической химии и биологии, исследования в области параметрической идентификации таких систем практически отсутствуют. Целью данной работы является разработка новых методов определения параметров задачи типа Коши для дифференциальных уравнений c дробными производными Римана-Лиувилля порядка α ∈ (1, 2). В работах [12, 13] рассмотрен метод параметрической идентификации дифференциальных ISSN: 2310-7081 (online), 1991-8615 (print); doi: http://dx.doi.org/10.14498/vsgtu1272 © 2014 Самарский государственный технический университет. Образец цитирования: А. С. О в с и е н к о, “Разработка методов идентификации параметров дифференциальных уравнений с дробной производной Римана-Лиувилля” // Вестн. Сам. гос. техн. ун-та. Сер. Физ.-мат. науки, 2014. № 1 (34). С. 134-144. 134 Разработка методов параметрической идентификации дифференциальных уравнений . . . уравнений дробного порядка α: α ∈ (0, 1) на основе построения аппроксимационного решения для исследуемого процесса. Однако в случае идентификации фрактальных осцилляторов порядка α ∈ (1, 2) оценки, вычисленные при помощи метода, описанного в [13], имеют достаточно высокую погрешность, что приводит к необходимости разработки новых процедур оценивания. Постановка задачи. Рассмотрим дробное осцилляционное уравнение следующего вида: α D0t y - µy = 0, (1) где α D0t y = t 1 d Γ(1 - α) dt 0 y(x)dx , (t - x)α t > 0, 1<α<2 - левосторонний дробный дифференциальный оператор Римана-Лиувилля порядка α [9]. Для уравнения (1) рассмотрим условия типа Коши: α-1 lim D0t y = a1 , t→0 α-2 lim D0t y = a2 . t→0 (2) Уравнение (1) является широко известным уравнением Баррета [14], которое неоднократно встречалось в литературе в связи с множеством приложений. Уравнением (1), в частности, может быть описано дробное обобщение закона Мальтуса, который хорошо известен в биологии и демографии [3, 10]. Это уравнение возникает также в задачах дробных реологических моделей [15, 16]. Решение задачи типа Коши (1), (2) можно представить в виде [14, 8] y (t) = a1 tα-1 Eα (µtα ; α) + a2 tα-2 Eα (µtα ; α - 1), ˜ где ∞ Eα (µtα ; ν) = n=0 (µtα )n , Γ(nα + ν) (3) µ, ν ∈ R - функция типа Миттаг-Леффлера. Согласно [17], функция типа Миттаг- Леффлера, входящая в выражение (3), может быть представлена в виде Eα (z; α) = при z > 0, | arg z| 1 1 1 (1-α) exp z α + H(z) zα α γ, πα/2 < γ < min{π, πα}, α < 2, |z| → ∞, или Eα (z; α) = H(z) при z < 0, γ (4) | arg z| (5) π, где M H(z) = - k=1 z -k + O(|z|-1-M ). Γ(α - kα) Представления (4) и (5) имеют место при α ∈ (0, 2); для α > 2 необходимо использовать другие соотношения [4], что значительно усложняет дальнейшую процедуру идентификации. В рамках данной работы ограничимся рассмотрением случая α ∈ (1, 2). 135 А. С. О в с и е н к о Задачу идентификации, решаемую в данной работе, можно сформулировать следующим образом: необходимо вычислить оценки параметров µ, a1 , a2 процесса (3) при известном порядке α ∈ (1, 2) дробного дифференциального оператора, входящего в уравнение (1). Идентификация параметров на основе линейно-параметрических дискретных моделей. Пусть µ > 0 и имеет место формула (4). Решение (3) задачи типа Коши (1), (2) после элементарных преобразований можно представить в виде 1 1 1 -1 µ α exp (µtα ) α + α m2 Γ(kα + 1) sin(πkα) -1-kα Γ(kα + 2) sin(π + πkα) -2-kα t + a2 t + k+1 πµ πµk+1 1 y (t) = a1 + a2 µ α -1 ˜ m1 + a1 k=1 k=1 + O(|µtα |-1-m ), (6) где m = min{m1 , m2 }. На основе представления (6) можно построить эффективную аппроксимацию процесса, описываемого при помощи (1): y (t) = c0 ept + c1 ts1 + d1 ts2 , ˆ (7) где 1 a1 a2 a1 + a2 µ α -1 1 -1 µ α , c1 = , d1 = , c0 = 2 α Γ(-α)µ Γ(-1 - α)µ2 1 p = µα , s1 = -1 - α, s2 = -2 - α. Абсолютная погрешность аппроксимации (7) составляет O(|µ-3 t-2α-1 |). Задача параметрической идентификации процесса, описываемого выражением (7), может быть решена при помощи линейно-параметрических дискретных моделей (ЛПДМ) [18]. Переобозначая в выражении (7) параметры s1 = s и s2 = s - 1, получим для функции (7) ЛПДМ следующего вида: y1 = λ6 , ˆ yk = λ1 yk-1 + λ2 (k + l)s + λ3 (k + l - 1)s + λ4 (k + l)s-1 + ˆ ˆ +λ5 (k + l - 1)s-1 , k = 2, 3, . . . , N, где λ1 = exp(pτ l), λ2 = c1 τ s , λ3 = -λ1 λ2 , λ4 = d1 τ s-1 , λ5 = -λ1 λ4 , λ6 = λ1+l c0 + λ2 (1 + l)s + λ4 (1 + l)s-1 . 1 При использовании среднеквадратичного критерия аппроксимации [18] функции (3) моделью (7) на конечном множестве точек tk = τ (k + l), k = = 1, 2, . . . , N , где N - объём выборки, τ - период дискретизации функции (7), l - целое число, минимизируется функционал e 2 = y-y ˆ 2 → min или в развёрнутой форме N N 2 (ek )2 → min, (yk - yk ) = ˆ k=1 136 k=1 Разработка методов параметрической идентификации дифференциальных уравнений . . . где yk - значения функции (3), вычисленные при помощи полученных оценок ˆ коэффициентов, yk - значения, полученные в результате эксперимента, e - вектор остатков. С учётом этого ЛПДМ в форме стохастических разностных уравнений можно представить в виде   y1 = λ 6 + e 1 ,  yk = λ1 yk-1 + λ2 (k + l)s + λ3 (k + l - 1)s + +λ4 (k + l)s-1 + λ5 (k + l - 1)s-1 + ηk ,   ηk = ek - λ1 ek-1 , k = 2, 3, . . . , N. В матричной форме ЛПДМ будет иметь вид b = F λ + η, η = Pλ e; где η = (η1 , η2 , . . . , ηN ) - N -мерный вектор эквивалентного возмущения в стохастическом разностном уравнении; Pλ - матрица эквивалентного возмущения; b = (y1 , y2 , . . . , yN ) - N -мерный вектор правой части; F - (N ×6)матрица регрессоров; ek = yk - yk - вектор остатков; λ = (λ1 , λ2 , . . . , λ6 ) - ˆ вектор коэффициентов ЛПДМ. Оценки коэффициентов ЛПДМ на начальном этапе реализации итерационной процедуры [19] вычисляются при помощи метода наименьших квадратов по формуле ˆ λ(0) = (F T F )-1 F b, ˆ где λ - вектор оценок коэффициентов стохастического разностного уравнения. В соответствии с алгоритмом данной процедуры формируется последовательность приближённых решений ˆ λ(s) = (F Ω-1 F )F Ω-1 b, ˆ (s-1) ˆ (s-1) λ λ s = 1, 2, . . . , 6, (8) ˆ где λ - вектор оценок коэффициентов стохастического разностного уравнения, вычисленный на s-том шаге, -1 Ω-1 ˆ ˆ ˆ (s-1) = (Pλ(s-1) Pλ(s-1) ) λ - квадратная симметричная матрица размера N ×N . Последовательность (8) сходится к λ - решению матричного уравнения F (Pλ Pλ )-1 F λ = F (Pλ Pλ )-1 b ˆ ˆ ˆ ˆ ˆ (см., например, [18]). По найденным оценкам λj , j = 1, 2, . . . , 6, могут быть вычислены параметры процесса (3). Для того чтобы привести в соответствие размерность вектора коэффициентов ЛПДМ и число подлежащих оцениванию параметров процесса, предложено использовать две дополнительные итерационные процедуры:  ˆ ˆ  y1 = λ1 λl c0 + λ2 (1 + l)s + λ4 (1 + l)s-1 + e1 ,  1  yk = λ1 yk-1 + λ2 (k + l)s + λ3 (k + l - 1)s + +λ4 (k + l)s-1 + λ5 (k + l - 1)s-1 + ηk ,    ηk = ek - λ1 ek-1 , k = 2, 3, . . . , N, 137 А. С. О в с и е н к о где λ1 = λ1 , λ2 = λ2 , λ3 = λ3 , λ4 = λ4 , λ5 = λ5 , и   y1 = λ4 + e1 ,   ˆ ˆ yk = λ1 (yk-1 - λ2 (k + l - 1)s - λ4 (k + l - 1)s-1 )+ s + λ (k + l)s-1 + η , +λ2 (k + l)  k 3   ηk = ek - λ1 ek-1 , k = 2, 3, . . . , N, где λ1 = λ1 , λ2 = λ2 , λ3 = λ4 , λ4 = λ6 . Выражения для вычисления оценок коэффициентов модели (7) с учётом дополнительных итерационных процедур будут иметь следующий вид: 1 ˆ ˆ ˆ ˆ ln λ1 , c1 = λ2 τ -s+1 , d1 = λ3 τ -s , ˆ τ ˆ ˆ ˆ λ - λ2 (1 + l)s - λ3 (1 + l)s-1 . c0 = 4 ˆ ˆ λ1 (1 + l) p= ˆ (9) Таким образом, оценки параметров исследуемого процесса (3) могут быть вычислены следующим образом: µ = pα , ˆ ˆ a1 = ˆ πˆ1 µ2 c ˆ , Γ(α + 1) sin(πα) a2 = ˆ c0 α ˆ µ ˆ 2 -1 α - a1 ˆ 1 . µα ˆ Применение метода минимизации невязки при определении параметров процесса. Однако применение описанного выше метода при решении задач идентификации дробно-дифференциальных операторов порядка α ∈ (1, 2) обеспечивает хорошие результаты при определении оценки параметра µ, но не может быть использовано при идентификации параметров a1 и a2 . В связи с этим необходимы другие подходы к решению поставленной задачи идентификации. Одним из таких подходов может быть описываемый далее метод определения параметров дробно-дифференциального уравнения на основе применения процедуры минимизации невязки. Сущность данного метода заключается в применении оптимизационных процедур не только для вычисления оценок коэффициентов ЛПДМ, но и в дальнейшем при определении аналитических соотношений между параметрами задачи (1), (2) и вычисленными оценками коэффициентов аппроксимационного решения, что позволяет рассматривать предлагаемую методику в целом в качестве нового подхода к идентификации дробных дифференциальных уравнений порядка α ∈ (1, 2). Оценки коэффициентов дробного осцилляционного уравнения (1) могут быть вычислены при помощи метода минимизации невязки [20]. Суть метода заключается в следующем: ставится задача нахождения параметров функции невязки α Ψ(µ, a1 , a2 ) = D0t y - µˆ ˆ y таким образом, чтобы невязка в некотором смысле имела минимальное значение; в качестве такого критерия выбран N α (D0t yk - µˆk )2 → min, ˆ y I(µ, a1 , a2 ) = k=1 138 (10) Разработка методов параметрической идентификации дифференциальных уравнений . . . при этом необходимо также выполнение условий типа Коши: α-1 lim D0t y = a1 , ˆ t→0 α-2 lim D0t y = a2 . ˆ t→0 (11) Очевидно, что для соблюдения критерия (10) необходимо выполнение условия (12): ∂I =2 ∂µ N α (D0t yk - µˆk )(-ˆk ) = 0, ˆ y y (12) k=1 откуда µ= N αˆ ˆ k=1 D0t yk · yk . N y 2 k=1 (ˆk ) Вычисляя производные дробного порядка, получаем α-2 ˆ D0t y (t) = 1 Γ(2 - α) α-1 ˆ D0t y (t) = t 0 y (x)dx ˆ = c0 t2-α E1 (pt; 3 - α)+ (t - x)α-1 Γ(s) Γ(s + 1) + d1 t1-α+s , + c1 t2-α+s Γ(3 - α + s) Γ(2 - α + s) d α-2 c0 Γ(s + 1) ˆ D0t y (t) = α-1 E1 (pt; 2 - α) + c1 t1-α+s , dt t Γ(2 - α + s) d α-1 α ˆ D0t y = D0t y = c0 t-α E1 (pt; 1 - α). ˆ dt Используя оценки коэффициентов аппроксимирующей функции, вычисленные по формулам (9), из соотношений (11) и (12) получим a1 = c1 Γ(α), ˆ ˆ µ= ˆ ˆ a2 = d1 Γ(α - 1), ˆ N s s-1 )(τ k)-α ˆ ˆ p c p ˆ k=1 c0 E1 (ˆτ k; 1 - α)(ˆ0 exp(ˆτ k) + c1 (τ k) + d1 (τ k) . N s s-1 )2 ˆ c p ˆ k=1 (ˆ0 exp(ˆτ k) + c1 (τ k) + d1 (τ k) Задача идентификации параметров (1), (2) в случае µ < 0 может быть решена аналогично описанной выше методике. С учётом (4) и (5) решение (3) может быть аппроксимировано при помощи функции следующего вида: y (t) = c1 ts1 + c2 ts2 + c3 ts3 , ˆ (13) где s1 = α - 1, s2 = α - 2. Методика построения ЛПДМ для аппроксимации (13) аналогична описанной выше. Производные дробного порядка вычисляются следующим образом: Γ(s1 + 1) + Γ(3 - α + s1 ) Γ(s2 + 1) Γ(s3 + 1) + c2 t2-α+s2 + c3 t2-α+s3 , Γ(3 - α + s2 ) Γ(3 - α + s3 ) α-2 2-α D0t y (t) = I0t y (t) = c1 t2-α+s1 ˆ ˆ 139 А. С. О в с и е н к о α-1 D0t y (t) = ˆ Γ(s3 + 1) d α-2 Γ(s1 + 1) D0t y (t) = c1 t1-α+s1 ˆ + c3 t1-α+s3 , dt Γ(2 - α + s1 ) Γ(2 - α + s3 ) d α-1 Γ(s3 + 1) α D0t y (t) = D0t y (t) = c3 ts3 -α ˆ ˆ . dt Γ(1 - α + s3 ) С помощью итерационной процедуры могут быть вычислены оценки коэффициентов λj , j = 1, 2, 3, которые известным образом связаны с параметрами модели (13): ˆ ˆ c1 = (λ2 - λ1 )τ -s1 , ˆ ˆ c2 = λ1 τ -s2 , ˆ ˆ c3 = λ2 τ -s3 . ˆ Функциональная зависимость между параметрами решения и параметрами аппроксимирующей модели может быть получена при помощи метода минимизации невязки [20]: a1 = c1 Γ(α), ˆ ˆ µ = c3 ˆ ˆ a2 = c2 Γ(α - 1), ˆ ˆ N s3 -α Γ(ˆ + 1)(ˆ (τ k)s1 + c (τ k)s2 + c (τ k)s3 ) ˆ ˆ ˆ ˆ s3 c1 ˆ2 ˆ3 k=1 (τ k) . N s1 ˆ s2 ˆ s3 2 ˆ s c ˆ ˆ k=1 Γ(ˆ3 + 1 - α)(ˆ1 (τ k) + c2 (τ k) + c3 (τ k) ) Результаты. Для реализации предложенных алгоритмов был разработан программный комплекс на языке Object Pascal, позволяющий на основе экспериментальных данных с учётом случайной аддитивной помехи в результатах наблюдений производить оценку параметров исследуемого процесса, описываемого задачей (1), (2). Проведены численные исследования эффективности предложенных методов, основанные на определении погрешности вычисления оценок параметров задачи (1), (2). Целью исследований являлось установление зависимости смещения полученных оценок параметров и среднеквадратического отклонения модели от периода дискретизации τ , порядка дробного дифференциального оператора α и величины случайной помехи ε в результатах измерений. Результаты численных исследований представлены на рис. 1-3. Цифры на графиках обозначают методы, для которых получены результаты: 1 - метод асимптотического разложения, 2 - метод идентификации на основе процедуры минимизации невязки. По результатам исследований можно сделать вывод о том, что применение метода минимизации невязки при решении задачи определения параметров дифференциальных уравнений с дробной производной позволяет существенно снизить погрешность определения оценок параметров a1 и a2 , являющихся начальными условиями задачи типа Коши. В то же время при реализации процедуры оценивания необходимо принимать во внимание и результаты, полученные с помощью других методов, что позволяет минимизировать величину среднеквадратического отклонения модели от входных данных. Между тем эффективное оценивание параметров исходного процесса при µ < 0 становится возможным лишь с использованием процедуры минимизации невязки. Таким образом, разработаны новые методы определения параметров дифференциального уравнения с дробными производными на основе разностных уравнений. Предложена аппроксимация решения дифференциального уравнения с дробными производными порядка α ∈ (1, 2), позволяющая свести 140 Разработка методов параметрической идентификации дифференциальных уравнений . . . a b Рис. 1. Зависимость среднеквадратического отклонения результатов, полученных по методу асимптотического разложения (кривая 1) и методу идентификации на основе процедуры минимизации невязки (кривые 2), от порядка дробного дифференциального оператора при µ > 0 (a) и µ < 0 (b) Figure 1. Standard deviation of results obtained from the asymptotic expansions method (the curve with the number 1), and the identification method based on residual minimization procedure (the curves with the number 2) in dependence on order of fractional derivative for µ > 0 (a), and µ < 0 (b) a b Рис. 2. Зависимость смещения оценок параметров µ (a) и a2 (b) от порядка дробного ˆ ˆ дифференциального оператора при µ > 0: 1 - метод асимптотического разложения, 2 - метод идентификации на основе процедуры минимизации невязки Figure 2. Bias of an estimator for parameter µ (a), and parameter a2 (b) in dependence on order ˆ ˆ of fractional derivative when µ > 0; the curve with the number 1 corresponds to the asymptotic expansions method; the curves with the number 2 correspond to the identification method based on residual minimization procedure 141 А. С. О в с и е н к о a b Рис. 3. Зависимость относительной дисперсии оценок параметров a1 (a) и µ (b) от величины случайной аддитивной помехи при µ < 0 для метода идентификации на основе процедуры минимизации невязки Figure 3. Variance of parameter a1 (a), and parameter µ (b) estimations in dependence ˆ ˆ on additional stochastic error when µ < 0 for the identification method based on residual minimization procedure исходную задачу идентификации к определению параметров аппроксимирующей функции. Построены ЛПДМ, связывающие в форме разностных уравнений результаты наблюдений мгновенных ординат процесса, описываемого дифференциальным уравнением с дробными производными. Предложенный подход позволяет вычислить оценки коэффициентов решения на основе среднеквадратического оценивания коэффициентов ЛПДМ. Получены соотношения между коэффициентами аппроксимирующей функции и параметрами аналитического решения дробного дифференциального уравнения на основе применения метода минимизации невязки. Численные исследования показывают, что описанные методы обладают достаточно высокой эффективностью и могут быть использованы при решении задач параметрической идентификации систем, содержащих дробные дифференциальные операторы.

About the authors

Anna S Ovsienko

Samara State Technical University

Email: sanabella@yandex.ru
Postgraduate Student, Dept. of Applied Mathematics & Computer Science 244, Molodogvardeyskaya st., Samara, 443100, Russian Federation

References

Supplementary files