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

Abstract


The methods for parametric identification of fractional differential operators with $\alpha \on (1, 2)$ degree according to instantaneous values of experimental observations based on the Barrett differential equation example are suggested. The methods are based on construction of the linear parametrical discrete model for fractional differential equation. The coefficients of the model are associated with the required parameters of differentiation equation of fractional order. Different approaches to the determination of the relationships between the parameters of the differential equation and the discrete model coefficients are considered. Connection expressions for coefficients of linear parametrical discrete model and Cauchy type problem parameters to be identified are obtained. The algorithm of the method which let us reduce the problem to computation of mean-square estimates for coefficients of linear parametrical discrete model is described. Numerical investigations have been done; furthermore, their results let us conclude high efficiency of the methods.

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
244, Molodogvardeyskaya st., Samara, 443100, Russian Federation
Postgraduate Student, Dept. of Applied Mathematics & Computer Science

References

  1. K. B. Oldham, J. Spanier, The fractional calculus. Theory and applications of differentiation and integration to arbitrary order, Mathematics in Science and Engineering, vol. 111, New York, London, Academic Press, 1974, xiii+234 pp.
  2. K. S. Miller, B. Ross, An introduction to the fractional calculus and fractional differential equations, New York, Jon Wiley & Sons. Inc., 1993, xiii+366 pp.
  3. А. М. Нахушев, Уравнения математической биологии, М.: Высшая школа, 1995. 301 с.
  4. R. Gorenflo, F. Mainardi, “Fractional calculus, integral and differential equations of fractional order”, Fractals and Fractional Calculus in Continuum Mechanics, Courses and Lecture Notes, vol. 378, Wien, New York, Springer Verlag, 1997, pp. 223-276, arXiv: 0805.3823 [math-ph]. doi: 10.1007/978-3-7091-2664-6_5.
  5. I. Podlubny, Fractional differential equations. An introduction to fractional derivatives, fractional differential equations, to methods of their solution and some of their applications, Mathematics in Science and Engineering, vol. 198, San Diego, Academic Press, 1999, xxiv+340 pp.
  6. В. В. Учайкин, Метод дробных производных, Ульяновск: Артишок, 2008. 512 с.
  7. А. М. Нахушев, Элементы дробного исчисления и их применение, Нальчик: КБНЦ РАН, 2003. 299 с.
  8. А. М. Нахушев, Дробное исчисление и его применение, М.: Физматлит, 2003. 271 с.
  9. A. A. Kilbas, H. M. Srivastava, J. J. Trujillo, Theory and applications of fractional differential equations, North-Holland Mathematics Studies, Amsterdam, Elsevier, 2006, xv+523 pp.
  10. В. В. Васильев, Л. А. Симак, Дробное исчисление и аппроксимационные методы в моделировании динамических систем, Киев: НАН Украины, 2008. 256 с.
  11. S. Das, Functional Fractional Calculus, Berlin, Heidelberg, Springer-Verlag, 2011, ix+612 pp. doi: 10.1007/978-3-642-20545-3.
  12. А. С. Овсиенко, “Разработка метода идентификации параметров дробного дифференциального оператора” / Труды Восьмой Всероссийской научной конференции с международным участием. Часть 2 / Математическое моделирование и краевые задачи, Самара: СамГТУ, 2011. С. 195-201.
  13. А. С. Овсиенко, “Параметрическая идентификация задачи типа Коши для одного дробного дифференциального уравнения” // Вестн. Сам. гос. техн. ун-та. Сер. Физ.-мат. науки, 2012. No 1(26). С. 157-165. doi: 10.14498/vsgtu997.
  14. J. H. Barrett, “Differential equations of non-integer order”, Canad. J. Math., 1954, vol. 6, no. 4, pp. 529-541. doi: 10.4153/cjm-1954-058-2.
  15. Е. Н. Огородников, В. П. Радченко, Н. С. Яшагин, “Реологические модели вязкоупругого тела с памятью и дифференциальные уравнения дробных осцилляторов” // Вестн. Сам. гос. техн. ун-та. Сер. Физ.-мат. науки, 2011. No 1(22). С. 255-268. doi: 10.14498/vsgtu932.
  16. Е. Н. Огородников, “Некоторые аспекты теории начальных задач для дифференциальных уравнений с производными Римана-Лиувилля” // Вестн. Сам. гос. техн. ун-та. Сер. Физ.-мат. науки, 2010. No 5(21). С. 10-23. doi: 10.14498/vsgtu812.
  17. М. М. Джрбашян, Интегральные преобразования и представления функций в комплексной области, М.: Наука, 1966. 672 с.
  18. В. Е. Зотеев, Параметрическая идентификация диссипативных механических систем на основе разностных уравнений, ред. В. П. Радченко, М.: Машиностроение-1, 2009. 344 с.
  19. А. С. Овсиенко, “Идентификация параметров процесса аномальной диффузии на основе разностных уравнений” // Вычислительные технологии, 2013. Т. 18, No 1. С. 65-73.
  20. Н. Н. Калиткин, Численные методы, М.: Наука, 1978. 512 с.

Statistics

Views

Abstract - 14

PDF (Russian) - 1

Cited-By


PlumX

Dimensions

Refbacks

  • There are currently no refbacks.

Copyright (c) 2014 Samara State Technical University

Creative Commons License
This work is licensed under a Creative Commons Attribution 4.0 International License.

This website uses cookies

You consent to our cookies if you continue to use our website.

About Cookies