Numerical method of estimation of parameters of the nonlinear differential operator of the second order

Full Text

Введение. Одной из основных проблем, возникающих при исследовании объектов, явлений или систем различной физической природы, является нелинейность математической модели. В этих случаях достоверная оценка параметров модели по результатам эксперимента или промышленных испытаний представляет собой сложную вычислительную задачу. Обычно эта проблема решается на основе линеаризации модели тем или иным способом, однако возникающая при этом систематическая погрешность в результатах вычислений параметров может оказаться недопустимо большой. В ряде случаев, например, когда характеристики нелинейности исследуемой системы являются важнейшим диагностическим признаком ее технического состояния и линеаризация принципиально недопустима, приходится применять известные методы нелинейного оценивания, которые не всегда дают положительный результат. В данной работе рассматривается нелинейный дифференциальный оператор второго порядка L ≡ my (t) + br y(t), y (t) + cy(t), описывающий широкий класс систем с одной степенью свободы со слабой нелинейностью общего вида, где m, b и c - коэффициенты пропорциональности; функция r y(t), y (t) описывает нелинейность общего вида. В частности, такой оператор широко используется при описании затухающих свободных колебаний различного рода механических систем с диссипативными силами, пропорциональными n-ной степени скорости движения (внутреннее трение в материале; конструкционное трение в опорах, шарнирах, сочленениях; силы сопротивления жидкой и газообразной среды; силы, возникающие с нагружением поглотителей энергии и т.п.): L ≡ my (t) + by (t)|y (t)|n-1 + cy(t) = 0, (1) где m - масса системы; cy(t) - линейная восстанавливающая сила; слагаемое by (t)|y (t)|n-1 описывает нелинейную диссипативную силу, пропорциональную n-ной степени скорости движения [1-4]; c и b - коэффициенты пропорциональности, n - параметр, характеризующий нелинейность диссипативной силы. Одной из важнейших задач при построении математической модели технической системы с использованием этого нелинейного дифференциального оператора является задача достоверной оценки его параметров по результатам наблюдений реакции системы на некоторое типовое тестовое воздействие [5-7]. Широко используемый в практике эксперимента цифровой спектральный анализ и высокоэффективные методы вибродиагностики, а также методы корреляционного анализа [8-13] по своей природе не позволяют обеспечить требуемую точность и достоверность оценок характеристик нелинейности диссипативной силы, которые являются важнейшим диагностическим признаком технического состояния механической системы. В то же время известные и традиционно применяемые методы определения характеристик рассеяния энергии колебаний различных механических конструкций и демпфирующих свойств материалов уже не вписываются в формат современных информационных технологий [2, 4, 5]. В данной работе проблема параметрической идентификации нелинейного дифференциального оператора второ557 З о т е е в В. Е., С т у к а л о в а Е. Д., Б а ш к и н о в а Е. В. го порядка решается на основе среднеквадратичного оценивания параметров импульсной характеристики системы по результатам наблюдений, полученных в ходе научно-технического эксперимента. Построение импульсной переходной характеристики методом Ван-дер-Поля. В [14] отмечается, что импульсная переходная функция (характеристика) g(t) системы, описываемой нелинейным дифференциальным оператором, может быть получена из решения однородного дифференциального уравнения my (t) + br y(t), y (t) + cy(t) = 0 с начальными условиями вида y(0) = 0, y (0) = 1/m. Рассмотрим построение приближенного решения этого нелинейного дифференциального уравнения методом Ван-дер-Поля с учетом малости параметра b [15]. В соответствии с этим методом решение ищется в виде y(t) = a(t) cos φ(t), φ(t) = ωt + β(t), (2) где a(t) и β(t) - некоторые функции, мало изменяющиеся за период колебаний T = 2π/ω. При этом предполагается выполнение равенства [15] y (t) = -ωa(t) sin φ(t), (3) то есть, чтобы y(t) и y (t) имели тот же вид, что и при отсутствии в системе слабой нелинейности общего вида br y(t), y (t) , когда b = 0. Дифференцируя (2) и сравнивая результат с (3), получаем равенство a (t) cos φ(t) - a(t)β (t) sin φ(t) = 0. (4) Вычисляя вторую производную y (t) из (3), подставляя полученный результат и функцию (2) в исходное дифференциальное уравнение, с учетом малости величин a (t), β (t) и b получаем два соотношения: mω 2 a(t) cos φ(t) = c cos φ(t), a (t) sin φ(t) + a(t)β (t) cos φ(t) = b r a(t) cos φ(t), -ωa(t) sin φ(t) . mω Из первого равенства следует, что ω = c/m, а второе вместе с (4) образует систему уравнений, из решения которой выражаются функции a (t) и a(t)β (t):    a (t) = b sin φ · r a cos φ, -ωa sin φ , mω (5) b   a(t)β (t) = cos φ · r a cos φ, -ωa sin φ . mω В правых частях (5) и далее для краткости обозначено a = a(t), ω = ω(t). С учетом слабой нелинейности (малости параметра b) функции a(t) и β(t) и их производные за время периода колебаний T практически не изменяются. Усредняя (интегрируя по переменной t) на промежутке времени [0; T ] левую и правую части равенств (5), с учетом dt = dφ/ω и соответствующей 558 Численный метод оценки параметров нелинейного дифференциального оператора. . . замены промежутка интегрирования в правой части с [0; T ] на [0; 2π], получаем систему уравнений, из решения которой можно найти функции a(t) и β(t), входящие в уравнение (2), описывающее импульсную характеристику:  2π b   r(a cos φ, -ωa sin φ) sin φ dφ,  a (t) = 2πmω 0 (6) 2π b    a(t)β (t) = r(a cos φ, -ωa sin φ) cos φ dφ. 2πmω 0 В рассматриваемом нами частном случае при нелинейной диссипативной силе, пропорциональной n-ной степени скорости движения by (t)|y (y)|n-1 , имеем: r y(t), y (t) = -aω sin φ| - aω sin φ|n-1 = -an ω n sin φ| sin φ|n-1 . При этом интегралы в (6) имеют вид 2π 2π | sin φ|n+1 dφ r(a cos φ, -ωa sin φ) sin φ dφ = -an ω n 0 0 и 2π 2π r(a cos φ, -ωa sin φ) cos φ dφ = -an ω n 0 sin φ| sin φ|n-1 cos φ dφ. 0 Очевидно, что в последнем выражении интегральная функция нечетная и имеет период равный 2π. Следовательно, 2π sin φ| sin φ|n-1 cos φ dφ = 0. 0 Обозначая 2π | sin φ|n+1 dφ, Jn = 0 с учетом последнего замечания систему (6) запишем в виде  bω n-1 n  a (t) = - a (t)Jn ; 2πm  a(t)β (t) = 0. (7) Обозначая a(0) = a0 и β(0) = β0 , из первого и второго уравнений системы (7) соответственно получаем a0 a(t) = , β(t) = β0 . bω n-1 n-1 n-1 1 + (n - 1) a Jn t 2πm 0 Отсюда с учетом (2) общее решение однородного нелинейного дифференциального уравнения (1) можно записать в виде a0 cos(ωt + β0 ) y(t) = n-1 bω n-1 n-1 1 + (n - 1) a Jn t 2πm 0 . (8) 559 З о т е е в В. Е., С т у к а л о в а Е. Д., Б а ш к и н о в а Е. В. Отметим, что полученное методом Ван-дер-Поля уравнение (8) свободных колебаний нелинейной диссипативной системы с силой трения, пропорциональной n-ной степени скорости движения, полностью совпадает с результатом, полученным на основе метода энергетического баланса в работах [1-4, 7]. С учетом начальных условий y(0) = 0, y (0) = 1/m из (8) получаем частное решение, описывающее импульсную характеристику системы с нелинейным дифференциальным оператором второго порядка (1): sin ωt g(t) = ωm n-1 b 1 + (n - 1) Jn t 2πmn . (9) Построенная математическая модель (9) импульсной характеристики системы с нелинейным дифференциальным оператором (1) содержит четыре параметра: m, n, b и ω = c/m, которые известным образом связаны с параметрами дифференциального оператора, что позволяет свести задачу параметрической идентификации дифференциального оператора к оценке параметров импульсной характеристики. Величину Jn при различных значениях n можно найти из таблиц, приведенных в работах [1-4, 7]. В частности, в табл. 1 представлены значения Jn для некоторых n. Значения n 0.0 Jn 4.00 n 1.6 Jn 2.83 интеграла Jn 0.2 0.4 3.77 3.58 1.8 2.0 2.75 2.67 Таблица 1 [The value of the integral Jn ] 0.6 0.8 1.0 1.2 1.4 3.42 3.27 3.14 3.03 2.92 2.2 2.4 2.6 2.8 3.0 2.60 2.53 2.47 2.41 2.37 При формировании выборки результатов измерений yk свободных колебаний нелинейной диссипативной системы в дискретные моменты времени, на основе которой вычисляются оценки параметров, момент первого измерения t0 , как правило, не совпадает с моментом t = 0 приложения импульса силы к системе. Поэтому преобразуем построенную модель (9) для системы координат, начало которой совпадает с моментом времени t0 первого отсчета в выборке результатов наблюдений: y(t) = g(t + t0 ). Заменяя в (9) переменную t на t + t0 , после простых алгебраических преобразований получаем математическую модель в форме уравнения свободных колебаний нелинейной диссипативной системы с силой сопротивления, пропорциональной n-ной степени скорости движения: y(t) = a0 cos(ωt + ψ0 ) n-1 где ω = δ0 1 + (n - 1) t T (10) c/m - частота колебаний; T = 2π/ω - период колебаний; δ0 = δ(a0 ) = 560 , bω n n-1 a Jn c 0 (11) Численный метод оценки параметров нелинейного дифференциального оператора. . . - логарифмический декремент колебаний; 1 a0 = ωm n-1 b Jn t0 1 + (n - 1) 2πmn , ψ0 = ωt0 - π 2 (12) - амплитуда и фаза колебаний, соответствующие моменту времени t0 первого отсчета в выборке результатов наблюдений относительно момента приложения импульса силы (t0 - время запаздывания). Из соотношений (11) и (12) можно получить формулы, которые позволяют по оценкам параметров n, δ0 , ω и a0 вычислить параметры m, b и c для нелинейного дифференциального оператора (1): m= n-1 2π - (n - 1)δ0 ωt0 , ω n-1 2πan-1 0 b= δ0 m , n-1 n-2 Jn a0 ω c = mω 2 . (13) Отметим, что формулы (13) предполагают известной величину t0 - время запаздывания момента первого отсчета y0 в выборке результатов наблюдений относительно момента времени приложения импульса силы. При неизвестной величине t0 ее можно найти из формулы (12) только при условии t0 ∈ [0; T ). Таким образом, задача параметрической идентификации нелинейного дифференциального оператора, описывающего диссипативную механическую систему с силами трения, пропорциональными n-ной степени скорости движения, может быть решена на основе оценки параметров свободных затухающих колебаний системы. Анализ известных методов оценки параметров нелинейной диссипативной системы. Известные методы оценки динамических характеристик нелинейной диссипативной системы по результатам наблюдений ее свободных колебаний описаны в [3-5, 7, 10, 12, 16]. К основным недостаткам большинства из описанных методов можно отнести, во-первых, то, что практически все применяемые при экспериментальных исследованиях методы предполагают линеаризацию математической модели, описывающей в той или иной форме динамический процесс в диссипативной механической системе. Это является одним из источников систематической погрешности в оценках динамических характеристик и для систем с сильной нелинейностью вообще недопустимо. Во-вторых, большинство из известных методов за время одной реализации затухающих колебаний не предусматривает получения избыточного числа измерений, используемых для вычисления одной оценки динамических характеристик. Тем самым исключается возможность обеспечения помехозащищенности этих оценок за счет применения статистических методов обработки экспериментальных данных. В-третьих, алгоритмы вычислений предполагают достаточно громоздкие промежуточные графические построения или функциональные дополнительные преобразования колебаний с помощью различных технических устройств [4, 10, 12]. Во всех случаях это вносит дополнительные погрешности в результаты вычислений оценок динамических характеристик и не способствует эффективному применению вычислительной техники при экспериментальных исследованиях [7]. Устранить эти недостатки можно только на основе применения статистических методов 561 З о т е е в В. Е., С т у к а л о в а Е. Д., Б а ш к и н о в а Е. В. обработки результатов эксперимента на базе компьютеризации экспериментальных исследований. Наиболее эффективно эта задача может быть решена на основе вычисления среднеквадратичных оценок параметров модели (10), минимизирующих величину отклонения этой модели от данных эксперимента в евклидовом нормированном пространстве: N -1 (yk - yˆk )2 → min, y - yˆ = k=0 где y - N -мерный вектор результатов эксперимента, yˆ - N -мерный вектор результатов расчета по модели (10). В такой постановке задачу параметрической идентификации дифференциального оператора можно интерпретировать как задачу нелинейного оценивания (нелинейной регрессии). Однако результаты научных исследований показали, что известные методы нелинейного оценивания: метод Ньютона-Гаусса, метод Хартли, метод Левенберга- Маркуардта, а также такие общие методы минимизации функций, как градиентный и сопряженных градиентов [17-20], не позволяют находить устойчивые оценки параметров модели (10), обеспечивающие минимум квадрата нормы ее отклонения от данных эксперимента. Основными причинами этого являются проблема выбора вектора начальных оценок параметров и требование удовлетворения на каждой итерации уточнения параметров ограничениям, накладываемым видом функциональной зависимости (10). В работе [7] эта проблема решалась на основе аппроксимации модели (10) функциональной зависимостью вида yˆ(t) = a0 cos(ωt + ψ0 ) δ0 n δ0 1+ t+ 1- T 2 T 2 . t2 На основе этой зависимости были построены разностные уравнения, описывающие результаты эксперимента, и получена линейная регрессионная модель, коэффициенты которой известным образом связаны с параметрами свободных колебаний диссипативной системы [7, 22]. Такой подход позволил обеспечить высокую адекватность аппроксимативной модели результатам наблюдений, и, как следствие, достоверные оценки ее параметров. Однако аппроксимация модели (10) неизбежно приводит к методической систематической погрешности в результатах расчета параметров режима свободных колебаний диссипативной системы, а значит, и в оценках параметров нелинейного дифференциального оператора (1). На основе численно-аналитического моделирования проведены исследования погрешностей, обусловленных аппроксимацией модели (10). Численный эксперимент проводился следующим образом. При времени наблюдения tH = 10 и периоде дискретизации τ = 0.2 по формуле a0 ak = , δ0 n-1 1 + (n - 1) τ k T где a0 = 1, T = 1, k = 0, 1, 2, . . . , N - 1, формировались выборки результатов расчета значений огибающей амплитуд колебаний (объем выборок составил 562 Численный метод оценки параметров нелинейного дифференциального оператора. . . N = 51). Были установлены следующие значения уровней затухания колебаний δa = a(tH )/a0 за время tH = 10: δa = 0.1; 0.2; 0.25; 0.33; 0.5; 0.8. Исследования проводились для 17-ти значений параметра нелинейности n в диапазоне от 0 до 2.5. Величина декремента колебаний в уравнении огибающей амплитуд колебаний рассчитывалась по формуле δ0 = 1 - an-1 (1 - an-1 )(n - 1)tH таким образом, чтобы обеспечить при заданном значении n требуемую величину уровня затухания δa. Для оценки погрешности вычисления параметров n и δ0 использовалась аппроксимация уравнения огибающей амплитуд колебаний в виде a0 , a0 = 1, T = 1. (14) δ0 n δ0 2 2 1+ t+ 1- t T 2 T Оценки параметров модели (14) находились из условия минимизации остаточной суммы квадратов a ˆ(t) = N -1 (ak - a ˆk )2 → min, a-a ˆ = k=0 где N = 51, a ˆk = a ˆ(tk ), k = 0, 1, . . . , N -1. Величина ∆a = a - a ˆ / a , характеризующая среднеквадратическое отклонение модели огибающей амплитуд колебаний от ее аппроксимации (14), в диапазоне изменения n ∈ [0.0; 2.5] и δa ∈ [0.1; 0.8] не превышала 0.05. Результаты исследований в виде кривых зависимостей относительных погрешностей ∆n = |(n - n ˆ )/n| и ∆δ0 = |(δ0 - δˆ0 )/δ0 | представлены на рис. 1. Очевидно, что только в относительно узком диапазоне изменения параметра n от 1.3 до 2.1 погрешности ∆n и ∆δ0 могут считаться приемлемыми. Причем с уменьшением уровня δa 0.33 погрешности резко возрастают. Расчет величины отношения числа результатов вычисления параметров n и δ0 , относительная погрешность которых не превышает 0.05, к общему числу 150 результатов вычислений (при n ∈ [0.1; 2.5] и δa ∈ [0.1; 0.8]), сделанный на основе таблицы данных исследования, показал следующие результаты. Для параметра n эта величина составляет всего 0.3, а для параметра δ0 - 0.76. Таким образом, можно сделать вывод о том, что область применения численного метода оценки параметров режима свободных колебаний механической системы, предложенного в [7], существенно ограничена величиной и степенью нелинейности (параметры b и n) сил трения, действующих в диссипативной системе. Численный метод оценки параметров на основе разностных уравнений. В данной работе предлагается новый численный метод оценки параметров математической модели (10) по результатам эксперимента. Этот метод состоит из двух этапов. На первом этапе модель (10) аппроксимируется функциональной зависимостью вида yˆ(t) = c0 cos(ωt + ψ0 ) , 1 + c1 t + c2 t2 563 З о т е е в В. Е., С т у к а л о в а Е. Д., Б а ш к и н о в а Е. В. Рис. 1. Зависимости относительной погрешности ∆n и ∆δ0 от параметра нелинейности n дифференциального оператора при различных уровнях затухания колебаний δa за промежуток времени наблюдения tH [Figure 1. Dependences of a relative error ∆n and ∆δ0 from the parameter of nonlinearity of n the differential operator at various levels of attenuation of fluctuations δa for a period of observation of tH ] дискретные значения которой в точках tk = τ k, k = 0, 1, 2, . . . , N - 1, описываются формулой c0 cos(ωτ k + ψ0 ) . (15) 1 + c1 τ k + c2 τ 2 k 2 Здесь же по найденным среднеквадратичным оценкам cˆ0 , cˆ1 и cˆ2 параметров этой модели вычисляются значения огибающей амплитуд колебаний: yˆk = ak = cˆ0 , 1 + cˆ1 τ k + cˆ2 τ 2 k 2 k = 0, 1, 2, . . . , N - 1. (16) На втором этапе с учетом сформированной выборки результатов вычислений ak по формуле (16) на основе минимизации функционала N -1 a-a ˆ 2 (ak - a ˆk )2 → min, = k=0 где a0 , δ0 ω ˆ 1 + (n - 1) τk 2π находятся среднеквадратичные оценки параметров n и δ0 . Для этого при формировании линейной регрессионной модели, коэффициенты которой известным образом связаны с параметрами нелинейной зависимости n и δ0 , используется дифференцирование зависимостей, описывающих огибающую амплитуд колебаний a(t) в модели (10) и ее аппроксимацию a ˆ(t). a ˆk = n-1 564 Численный метод оценки параметров нелинейного дифференциального оператора. . . На первом этапе численного метода вначале на основе дискретной модели (15) строится система разностных уравнений, описывающих результаты эксперимента yk , k = 0, 1, 2, . . . , N - 1, где N - объем выборки результатов наблюдений с периодом дискретизации τ . Из формулы (15) можно получить аналогичные соотношения для yˆk-1 и yˆk-2 , и на их основе с помощью простых алгебраических преобразований построить систему разностных уравнений вида  yˆ0 = λ4 ;    yˆ1 = λ5 ;  yˆk-2 + yˆk = λ1 yˆk-1 + λ2 [(k - 2)ˆ yk-2 - λ1 (k - 1)ˆ yk-1 + k yˆk ]+  2 2 2  +λ3 [(k - 2) yˆk-2 - λ1 (k - 1) yˆk-1 + k yˆk ],   k = 2, 3, . . . , N - 1, коэффициенты которой связаны с параметрами модели (15) соотношениями λ2 = -c1 τ, λ3 = -c2 τ 2 , 1 λ5 = (λ4 cos ωτ - c0 sin ψ0 sin ωτ ). 1 - λ2 - λ3 λ1 = 2 cos ωτ, λ4 = c0 cos ψ0 , С учетом случайного разброса εk в результатах эксперимента yk относительно модели yˆk : yk = yˆk + εk построенную систему уравнений можно привести к виду  y0 = λ4 + ε0 ;    y1 = λ5 + ε0 ;    (i)   yk-2 + yk = λ1 yk-1 + λ2 [(k - 2)yk-2 - λ1 (k - 1)yk-1 + kyk ]+    (i) +λ3 [(k - 2)2 yk-2 - λ1 (k - 1)2 yk-1 + k 2 yk ] + ηk+1 ,  k = 2, 3, . . . , N - 1;   2 ]ε  η = [1 - λ (k - 2) - λ (k - 2) -  2 3 k+1 k-2   (i)   -λ1 [1 - λ2 (k - 1) - λ3 (k - 1)2 ]εk-1 + [1 - λ2 k - λ3 k 2 ]εk ,   k = 2, 3, . . . , N - 1, (17) (i) ˆ при i = 0, 1, 2, . . . является некоторой известной оценкой где величина λ 1 коэффициента λ1 . При такой интерпретации на основе системы разностных уравнений (17) можно построить линейную обобщенную регрессионную модель вида b = Fλˆ (i) λ + η; η = Pλˆ (i) ε, (18) элементы которой описываются следующими соотношениями: λ = (λ1 , λ2 , λ3 , λ4 , λ5 ) - вектор неизвестных коэффициентов; b = (y0 , y1 , y0 + y2 , y1 + y3 , . . . , yN -3 + yN -1 ) , - вектор свободных членов; ε = (ε0 , ε1 , ε2 , . . . , εN -1 ) - вектор случайной помехи в результатах эксперимента; 565 З о т е е в В. Е., С т у к а л о в а Е. Д., Б а ш к и н о в а Е. В.  Fλˆ(i) 0 0 y1 y2 y3 .. .        =    yN -2   0 0 (i) -λ1 y1 + 2y2 (i) y1 - λ1 2y2 + 3y3 (i) 2y2 - λ1 3y3 + 4y4 .. . 0 0 (i) -λ1 y1 + 4y2 (i) y1 - λ1 4y2 + 9y3 (i) 4y2 - λ1 9y3 + 16y4 .. . (N - 3)yN -3 - (i) -λ1 (N - 2)yN -2 + +(N - 1)yN -1 (N - 3)2 yN -3 - (i) -λ1 (N - 2)2 yN -2 + +(N - 1)2 yN -1 1 0 0 0 0 .. .  0 1  0  0  0  ..  .  0 0   - матрица регрессоров; η = (η1 , η2 , η3 , . . . , ηk , . . . , ηN ) - вектор невязки;  Pλˆ(i)      =     1 0 1 0 0 .. . 0 0 0 0 1 0 0 (i) (i) (i) (i) (i) -λ1 (1 - λ2 - λ3 ) 1 - λ2 - λ3 0 (i) (i) (i) (i) (i) (i) (i) 1 - λ2 - λ3 -λ1 (1 - 2λ2 - 4λ3 ) 1 - 3λ2 - 9λ3 (i) (i) (i) (i) (i) 0 1 - 2λ2 - 4λ3 -λ1 (1 - 3λ2 - 9λ3 ) .. .. .. . . . 0 0 0  ... 0 ... 0   ... 0   ... 0  (19)  ... 0   ..  ... . ... (i) (i) 1 - (N - 1)λ2 - (N - 1)2 λ3 - матрица линейного преобразования вектора случайной помехи в результатах эксперимента. Применяя алгоритм численного метода нелинейного оценивания на основе разностных уравнений [7, 21-23], с учетом второго уравнения в системе (18) преобразуем первое уравнение к виду P ˆ-1 b = P ˆ-1 F ˆ (i) λ + ε. (i) (i) λ λ λ (20) (i) Поскольку величина λˆ1 в левой части и в матрице уравнения регрессии (19) известна, это уравнение линейно относительно вектора коэффициˆ этого вектора, миниентов λ. Следовательно, среднеквадратичную оценку λ мизирующего сумму квадратов величин разброса данных эксперимента относительно модели (15): N -1 ε2k → min, k=0 можно найти из решения нормальной системы линейных алгебраических уравнений -1 ˆ Fλˆ (i) Ω-1 ˆ (i) λ = Fλ ˆ (i) Ω ˆ (i) b, ˆ (i) Fλ λ 566 λ Численный метод оценки параметров нелинейного дифференциального оператора. . . где Ωλ(i) = Pλ(i) Pλ(i) . Отсюда получаем ˆ i+1 = F (i) Ω-1 F ˆ (i) λ ˆ ˆ (i) λ λ λ -1 Fλˆ (i) Ω-1 ˆ (i) b, λ i = 0, 1, 2, . . . (21) Формула (20) описывает итерационный процесс уточнения оценки вектора коэффициентов обобщенной регрессионной модели (18). В качестве начального приближения оценки коэффициента λ1 можно выбрать либо величину ˆ (0) = λ 1 N -1 k=2 (yk-2 + yk )yk-1 , N -1 2 k=2 yk-1 которая минимизирует сумму квадратов отклонений N -1 [yk - c0 cos(ωτ k + ψ0 )]2 → min, k=0 ˆ полученный из уравнения регрессии либо первый элемент в оценке вектора λ, вида yk-1 + yk = λ1 yk-1 + λ2 [(k - 2)yk-2 + kyk ] + λ3 (k - 2)2 yk-2 + k 2 yk + + λ4 (k - 1)yk-1 + λ5 (k - 1)2 yk-1 + ηk , k = 2, . . . , N - 1, -1 2 при минимизации остаточной суммы квадратов N k=2 ηk → min . Вычисления по формуле (21) следует продолжать до выполнения условия ˆ (i+1) - λ ˆ (i) = max λ ˆ i+1 - λ ˆ i < θ, λ j j (22) j=1,3 где θ - произвольно малая положительная величина. Полученные при выполнении условия (22) оценки коэффициентов обобщенной регрессионной моˆ j , j = 1, 2, 3, 4, 5, минимизирующие дели (18) следует принять за оценки λ величину среднеквадратического отклонения модели (15) от данных эксперимента. ˆ j с учеПо найденным среднеквадратическим оценкам коэффициентов λ том соотношений (16) оценки параметров модели (15) можно найти по формулам ω ˆ= cˆ0 = ˆ1 ˆ2 ˆ3 1 λ λ λ arccos , cˆ1 = - , cˆ2 = - 2 , τ 2 τ τ ˆ ˆ ˆ ˆ ˆ λ1 λ4 - 2λ5 (1 - λ2 - λ3 ) , u= ˆ2 4-λ 1 ˆ1λ ˆ 4 - 2λ ˆ 5 (1 - λ ˆ2 - λ ˆ 3 )]2 [λ ˆ2, +λ 4 ˆ2 4-λ 1  u ˆ 4 > 0,  arctg , при λ   ˆ4 λ   u ˆ 4 < 0, + sign(u) · π, при λ ψˆ0 = arctg ˆ  λ4    π  ˆ 4 = 0. sign(u) · , при λ 2 ˆ 2 + u2 = λ 4 567 З о т е е в В. Е., С т у к а л о в а Е. Д., Б а ш к и н о в а Е. В. Завершается первый этап вычислением значений аппроксимации огибающей амплитуд колебаний в (15) по формуле (16): ak = cˆ0 . 1 + cˆ1 τ k + cˆ2 τ 2 k 2 Заметим, что по результатам вычислений на первом этапе можно найти приближенную оценку логарифмического декремента колебаний δ0 в (10), соответствующего начальной амплитуде колебаний a0 : 4π 2 λ2 + 2πλ2 arccos 2π(c1 ω 2 + 2πc2 ) δ0 = 2 = ω + 2πc1 ω + 4π 2 c2 4π 2 λ λ1 2 λ1 λ1 - arccos2 2 + 2πλ2 arccos 2 2 . На втором этапе предлагаемого численного метода находятся оценки параметров нелинейности системы n и δ0 . В основе вычисления оценок лежит минимизация функционала N -1 a-a ˆ 2 (ak - a ˆk )2 → min, = k=0 где элементы вектора a ˆ находятся по формуле a0 a ˆk = n-1 δ0 ω ˆ 1 + (n - 1) τk 2π . (23) Очевидно, что задача среднеквадратического оценивания параметров модели (23) по результатам расчета по формуле (16) относится к достаточно сложной задаче нелинейной регрессии [17-20]. В данной работе предлагается перейти к линейной регрессионной модели, коэффициенты которой известным образом связаны с параметрами n и δ0 . С этой целью предлагается вначале провести дифференцирование непрерывных в области своего определения функций: a(t) = cˆ0 , 1 + cˆ1 t + cˆ2 t2 a0 a ˆ(t) = n-1 δ0 ω ˆ 1 + (n - 1) t 2π . Находя производные, получаем cˆ0 (ˆ c1 + 2ˆ c2 t) a (t) = - , (1 + cˆ1 t + cˆ2 t2 )2 a ˆ (t) = - Очевидно, что имеет место равенство a(t) = a ˆ(t) + ε(t), 568 δ0 ω ˆ 2π δ0 ω ˆ 1 + (n - 1) t 2π a ˆ(t). Численный метод оценки параметров нелинейного дифференциального оператора. . . где непрерывная функция ε(t) описывает ошибку аппроксимации зависимости a ˆ(t) построенной моделью a(t). Дифференцируя это равенство, получаем a (t) = a ˆ (t) + ε (t). Теперь, переходя в область дискретного изменения аргумента t: tk = τ k, k = 0, 1, 2, . . . , N -1, с учетом полученных выше соотношений имеем δ0 ω ˆ 2π ak = - (ak - εk ) + εk , δ0 ω ˆ τk 2π где значения ak , k = 0, 1, 2, . . . , N -1, вычисляются по формуле (16), а зачения ak по формуле 1 + (n - 1) ak = cˆ0 (ˆ c1 + 2ˆ c2 τ k) , (1 + cˆ1 τ k + cˆ2 τ 2 k 2 )2 k = 0, 1, 2, . . . , N - 1. Отсюда после простых алгебраических преобразований, применяя формулы численного дифференцирования второго порядка εk = εk+1 - εk-1 , 2τ k = 1, 2, . . . , N - 2, и εN -3 - 4εN -2 + 3εN -1 , 2τ получаем систему уравнений вида    a0 = λ3 + ε0 ;   1 + kλ1 1 + kλ1    ak = -ak kλ1 - ak λ2 - 2τ εk-1 + λ2 εk + 2τ εk+1 ,   k = 1, 2, . . . , N - 2; 1 + (N - 1)λ1    aN -1 = -aN -1 (N - 1)λ1 - aN -1 λ2 + εN -3 -   2τ   3 2   [1 + (N - 1)λ1 ] εN -1 , - [1 + (N - 1)λ1 ] εN -2 + λ2 + τ 2τ εN -1 = (24) коэффициенты которой связаны с параметрами зависимости a(t) формулами λ1 = (n - 1) δ0 ω ˆ τ, 2π λ2 = δ0 ω ˆ , 2π λ3 = a ˆ0 . (25) Соответствующая системе уравнений (24) обобщенная регрессионная модель принимает вид b = F λ + η, (26) η = Pλˆ (i) ε, элементы которой описываются следующими соотношениями: λ = (λ1 , λ2 , λ3 ) - вектор неизвестных коэффициентов; ε = (ε0 , ε1 , ε2 , . . . , εN -1 ) - вектор отклонений модели a ˆ(t) от результатов расчета a(t); η = (η1 , η2 , . . . , ηN ) - вектор невязки; b = (a0 , a1 , a2 , . . . , aN -1 ) (27) 569 З о т е е в В. Е., С т у к а л о в а Е. Д., Б а ш к и н о в а Е. В. - вектор свободных членов;   1 0  0 0  ..  . 0 0 0  -a1 -a1  -2a2 -a2  F = -3a -a3  3  . .. ..  . -(N - 1)aN -1 -aN -1 (28) - матрица регрессоров; 1  1 + λ(i) 1  -  2τ    0   Pλ =   0   ..   .  0   0 0  0 (i) 1 + λ1 - 2τ (i) λ2 - (i) 2λ1 1+ 2τ 0 0 0 (i) (i) - λ2 (i) - 0 .. . 0 0 0 1 + 3λ1 2τ .. . 0 0 0 ... ... ... ... .. . 0 0 0 0 .. . ... 1 + (N - 3)λ1 2τ ... λ2 1 + 2λ1 2τ (i) λ2 .. . 0 0 0  0 0 0 0 .. . (i) 0 (i) (i) (i) ... - 2 1 + (N - 2)λ1 τ 1 + (N - 2)λ1 2τ (i) 3 1 + (N - 2)λ1 (i) λ2 + 2τ           (29)         - матрица линейного преобразования вектора ε отклонения момента от результатов расчета ak . Применяя алгоритм численного метода нелинейного оценивания на основе разностных уравнений [7, 21-23] с целью минимизации суммы квадратов отклонения модели (23) от результатов расчета по формуле (16): N -1 ε2k → min, k=0 можно получить рекуррентную формулу ˆ i+1 = [F (i) Ω-1 λ F (i) ]-1 Fλ(i) Ω-1 b, λ λ(i) λ λ(i) 570 i = 0, 1, 2, . . . , (30) Численный метод оценки параметров нелинейного дифференциального оператора. . . где матрица Ωλ(i) = Pλ(i) Pλ(i) находится с учетом формулы (29), а вектор b и матрица F определяются формулами (27) и (28) соответственно. Формула (30) описывает итерационный процесс уточнения оценки вектора коэффициентов регрессионной модели (26). Начальное приближение вектора оценок ˆ (0) можно найти из условия минимизации функционала λ η из соотношения 2 = b - Fλ 2 → min ˆ (0) = (F F )-1 F b. λ Вычисления по формуле (30) следует продолжать до выполнения условия (22): ˆ i+1 - λ ˆ i < θ, ˆ (i+1) - λ ˆ (i) = max λ λ j j j=1,3 где θ - заданная малая величина. По найденным среднеквадратическим оценкам коэффициентов с учетом соотношений (25) определяются оценки параметров модели (23): 2π ˆ δˆ0 = λ2 , ω ˆ n ˆ= ˆ1 λ + 1, ˆ2τ λ ˆ3. a ˆ0 = λ Второй этап численного метода завершается вычислением оценок параметров нелинейного дифференциального оператора (1) по формулам m= n ˆ -1 2π - (ˆ n - 1)δˆ0 ω ˆ t0 , n ˆ -1 n ˆ -1 ˆ 2πˆ a0 ω b= δˆ0 m ˆ , n ˆ -1 n ˆ ˆ -2 Jn a ˆ0 ω c = mˆ ˆ ω2. (31) Результаты численно-аналитических исследований. Для проверки эффективности предлагаемого численного метода проведены численно-аналитические исследования, а также выполнен сравнительный анализ точности результатов вычислений, полученных этим методом и численным методом, описанным в [7, 21, 22]. Имитационное компьютерное моделирование и реализация алгоритмов этих методов проводились при следующих параметрах нелинейного дифференциального оператора (1): m = 0.4, b = 0.2, n = 3 и c = 12. Такие значения параметров соответствуют достаточно высокой степени нелинейности и величине диссипативной силы, обуславливающей затухание свободных колебаний в системе. В соответствии с формулой (10) с периодом дискретизации τ = 0.05 формировались выборки результатов наблюдений yˆk , объем которых был равен N = 40. К результатам моделирования yˆk добавлялась случайная некоррелированная помеха εk , имеющая нормальный закон распределения с математическим ожиданием M [εk ] = 0 и дисперсией D[εk ] = σε2 : εk ∈ N (0, σε2 ), которая генерировалась таким образом, чтобы ее величина удовлетворяла равенству γ= ε = yˆ N -1 2 k=0 sk , N -1 2 ˆk k=0 y 571 З о т е е в В. Е., С т у к а л о в а Е. Д., Б а ш к и н о в а Е. В. где значения γ выбирались из диапазона от 0 до 0.05 равномерно с шагом 0.005. По результатам суммирования yk = yˆk + εk , k = 0, 1, 2, . . . , N - 1, в соответствии с алгоритмами двух различных методов вычислялись оценки m, ˆ ˆb, n ˆ и cˆ искомых параметров. При одном и том же значении мощности помехи γi формирование выборки результатов наблюдений yk повторялось 100 раз. Затем для каждого значения γi из диапазона от 0 до 0.05 вычислялись среднеквадратичные оценки погрешности. В качестве оценки погрешности был выбран второй центральный момент оценки параметра дифференциального оператора pˆ относительно его точного значения p: M [(p - pˆ)2 ] = D[ˆ p] + (M [ˆ p] - p)2 . Такая оценка погрешности наиболее информативна, так как она интегрировано включает в себя дисперсию оценки параметра D[ˆ p] = M [(ˆ p - M [ˆ p])2 ] и квадрат смещения математического ожидания этой оценки M [ˆ p] относительно его точного значения. При приближенном вычислении математического ожидания M [ˆ p] использовалось среднее арифметическое результатов расчета параметра в 100 реализациях при одном и том же значении γi мощности случайной помехи: 100 1 pˆj . M [ˆ p] ≈ 100 j=1 Результаты исследований в форме зависимостей относительной суммы квадратов ошибки 100 1 ∆p = (pj - pˆj )2 , j=1 10|p| где p любой из параметров m, c, b, n, от величины случайной помехи в данных эксперимента, представлены на рис. 2. По результатам исследований можно сделать следующие выводы. Погрешности оценок параметров следует анализировать дифференцированно в зависимости от этапа численного метода. На первом этапе вычисляются параметры, связанные с частотой колебаний ω, а также с характеристикой затухания этих колебаний δ0 , но без учета вида модели by (t)|y (t)|n-1 нелинейной зависимости диссипативной силы от скорости движения. Из кривых, представленных на рис. 2 (сверху), следует, что суммы квадратов ошибок ∆m и ∆c, во-первых, существенно зависят от величины случайной помехи в результатах наблюдений, и, во-вторых, с увеличением этой величины кривые 1 и 2 практически совпадают. Это объясняется тем, что алгоритмы вычислений оценок на первом этапе предлагаемого метода и известным методом [7, 21, 22] похожи и дают результаты с близкими дисперсиями. При этом проведенные исследования показали, что доминирующей составляющей в суммарной среднеквадратичной ошибке M [(p - pˆ)2 ], p ∈ {m, c}, является дисперсия D[ˆ p], а квадрат смещения (M [ˆ p]-p)2 математического ожидания оценки параметра относительно точного значения практически не зависит от ε и не превышает величины 0.26 для известного метода и 0.07 для нового численного метода от общей суммы квадратов ошибок (при γ = 0.05). 572 Численный метод оценки параметров нелинейного дифференциального оператора. . . Рис. 2. Сравнительный анализ относительных погрешностей вычисления параметров нелинейного дифференциального оператора (1 - результаты вычисления известным методом [7, 21, 22], 2 - результаты предлагаемым численным методом) [Figure 2. The comparative analysis of relative errors of calculation of parameters of the nonlinear differential operator (1-the results of calculation by the known method [7, 21, 22], 2-the results by an estimated numerical method)] Очевидно, что по графикам, представленным на рис. 2 (сверху), величину смещения оценок параметров m и c в относительных единицах можно оценить по точкам, соответствующим ε = 0, в которых D[ˆ p] = 0, p ∈ {m, c}. В целом, как видно из рис. 2 (сверху), погрешности вычисления параметров m и c невелики и имеют порядок случайной помехи в результатах наблюдений. Иной характер имеют зависимости суммарной среднеквадратичной погрешности ∆b и ∆n оценок параметров. В известном методе эти оценки вычисляются непосредственно по формулам, построенным при разработке алгоритма метода [7, 21, 22]. При этом аппроксимация нелинейной зависимости огибающей амплитуд колебаний в (10) функцией вида (14) обуславливает заметное смещение оценок параметров b и n, что приводит к недопустимо большим погрешностям в результатах вычислений (см. рис. 1). В отличие от оценок параметров m и c в суммарной среднеквадратичной погрешности параметров b и n: M [(p - pˆ)2 ], p ∈ {b, n}, доминирующим является квадрат смещения (M [ˆ p] - p)2 , а не дисперсия D[ˆ p], которая составляет величину, не 573 З о т е е в В. Е., С т у к а л о в а Е. Д., Б а ш к и н о в а Е. В. превышающую 0.02 от величины общей суммы квадратов ошибок. В алгоритме предлагаемого численного метода оценки параметров b и n находятся на втором этапе из условия минимизации величины среднеквадратического отклонения аппроксимации (16), построенной на предыдущем этапе, от дискретной модели (23). В новом предложенном численном методе, как показали результаты исследований, смещение M [ˆ p] - p мало, соизмеримо со среднеквадратическим отклонением σ = D[ˆ p], и составляет величину от 0.17 до 0.24 от смещения соответствующих оценок, полученных известным методом [7, 21, 22]. На рис. 2 (снизу) представлены результаты вычислений суммарной среднеквадратичной погрешности оценок параметров b и n: кривые 1 соответствуют результатам вычислений известным методом [7, 21, 22]; кривые 2 - результатам вычислений новым численным методом. Видно, что в относительных единицах погрешность оценки параметра b уменьшилась в 5 раз, а параметра n - в 3 раза. Таким образом, по результатам проведенных исследований можно сделать общий вывод о том, что применение разработанного численного метода существенно (в несколько раз) повышает точность вычисления оценок параметров нелинейного дифференциального оператора за счет устранения смещения в результатах вычислений. Результаты апробации численного метода при обработке данных эксперимента. Разработанный численный метод оценки параметров нелинейного дифференциального оператора (1) был апробирован при обработке результатов эксперимента, полученных на основе компьютерного моделирования. Значения параметров нелинейного дифференциального оператора (1) при имитационном моделировании были выбраны следующим образом: m = 0.4; b = 0.2; n = 3.0 и c = 12.0. На основе этих значений были рассчитаны параметры математической модели (10), которая описывает свободные колебания диссипативной системы: ω = c/m, затем при выбранном времени запаздывания t0 = 0.1 по формулам (12) вычислялись амплитуда a0 и фаза колебаний ψ, соответствующие моменту времени t0 , и в соответствии с (11) рассчитывалась величина декремента колебаний δ0 . Результаты этих расчетов представлены в последних столбцах второй строки табл. 2. Таким образом, имитационная модель колебаний диссипативной системы имела вид y˜(t) = 0.411 cos(5.477t - 1.023) √ . 1 + 1.889t (32) По формуле (32) при tk = τ k, k = 0, 1, 2, . . . , N - 1, с периодом дисТаблица 2 Результаты вычислений оценок параметров [The results of calculations of the parameter estimates] Parameter Differential operator Impulse characteristic m b n c δ0 ω n a0 ψ0 Exact value 0.4 0.2 3.0 12.0 1.09 5.48 3.0 0.41 -1.02 Parameter estimation 0.43 0.299 2.64 12.85 1.12 5.45 2.64 0.42 -0.99 Relative error 0.084 0.493 0.12 0.071 0.125 0.006 0.12 0.031 0.03 574 Численный метод оценки параметров нелинейного дифференциального оператора. . . кретизации τ = 0.05 были вычислены значений дискретной функции y˜k , tk = τ k, k = 0, 1, 2, . . . , 39. К этим значениям была добавлена случайная помеха (некоррелированные случайные величины εk ∈ N (0, σε2 ) с нулевым математическим ожиданием и одинаковыми дисперсиями, распределенные по нормальному закону), величина которой была сгенерирована таким образом, чтобы отклонение по евклидовой норме в относительных единицах ∆ε = y - y˜ y˜ составило величину ∆ε = 0.2. Полученные таким образом величины yk = y˜k + εk , k = 0, 1, 2, . . . , N - 1, (33) использовались далее в качестве смоделированной выборки результатов эксперимента. Затем на основе описанного выше алгоритма численного метода были найдены оценки параметров математической модели (10), которые представлены в последних столбцах третьей строки табл. 2. В результате была построена математическая модель колебаний диссипативной системы вида yˆ(t) = 0.424 cos(5.445t - 0.992) √ . 1.641 1 + 1.742t (34) На рис. 3 изображены результаты апробации метода на основе компьютерного моделирования: маркеры 1 соответствуют смоделированным результатам эксперимента yk со случайной помехой ∆ε = 0.2 (см. формулу (33)); Рис. 3. 1 - результаты модельного эксперимента yk (см. формулу (33)); 2 - динамический процесс y˜(t), смоделированный (точный) по формуле (32); 3 - динамический процесс yˆ(t), восстановленный по построенной модели (34) [Figure 3. 1-the results of the model experiment yk (see Eq. (33)); 2-the dynamic process y˜(t) simulated (exact) on Eq. (32); 2-the dynamic process yˆ(t) restored on the constructed model (see Eq. (34))] 575 З о т е е в В. Е., С т у к а л о в а Е. Д., Б а ш к и н о в а Е. В. кривая 2, построенная по формуле (32), описывает смоделированный «точный» динамический процесс y˜(t); кривая 3 описывает динамический процесс yˆ(t), соответствующий построенной модели (34). Полученные данные показывают, что даже при достаточно высоком уровне случайной помехи ∆ε = 0.2 восстановленная по результатам эксперимента кривая 3, описывающая динамический процесс в системе, практически совпадает с «точной» кривой 2, смоделированной при заданных значениях параметров нелинейного дифференциального оператора (1). Количественно это расхождение между построенной моделью (34) и имитационной «точной» моделью (32) можно оценить, используя евклидову норму в векторном пространстве: y˜ - yˆ = 0.051. yˆ Очевидно, что по сравнению с величиной ∆ε = 0.2 расхождение между данными эксперимента и результатами расчета по модели (34) незначительно. Найденные по формулам (31) оценки параметров нелинейного дифференциального оператора (1) приведены в первых столбцах третьей строки табл. 2. В последней строке табл. 2 представлены относительные погрешности результатов вычислений параметров нелинейного дифференциального оператора и параметров построенной модели (34) динамического процесса в системе, описываемой с помощью этого дифференциального оператора. Заключение. По результатам проведенных исследований численного метода оценки параметров нелинейного дифференциального оператора можно сделать следующие выводы. Разработан новый численный метод оценки параметров нелинейного дифференциального оператора с диссипативной силой, пропорциональной n-ной степени скорости движения. В основе метода лежит среднеквадратичное оценивание коэффициентов обобщенной регрессионной модели, построенной с учетом разностных уравнений, описывающих результаты измерений импульсной характеристики системы. Реализованная в методе двухэтапная процедура дифференцированного оценивания параметров динамического процесса позволяет обеспечить высокую адекватность построенной модели данным эксперимента за счет использования уравнения скорости затухания колебаний при разработке алгоритма вычислений на втором этапе метода. Применение разработанного численного метода позволяет существенно (в несколько раз) повысить точность оценок параметров нелинейного дифференциального оператора по сравнению с известными методами за счет устранения смещения в оценках, обусловленного применением аппроксимации при моделировании огибающей амплитуд колебаний. Конкурирующие интересы. Мы не имеем конкурирующих интересов. Авторский вклад и ответственность. Все авторы принимали участие в разработке концепции статьи и в написании рукописи. Авторы несут полную ответственность за предоставление окончательной рукописи в печать. Окончательная версия рукописи была одобрена всеми авторами.

About the authors

Vladimir E Zoteev

Samara State Technical University

Email: zoteev-ve@mail.ru
Dr. Tech. Sci.; Professor; Dept. of Applied Mathematics & Computer Science 244, Molodogvardeyskaya st., Samara, 443100, Russian Federation

Ekaterina D Stukalova

Samara State Technical University

Email: stukalova.ed@samgtu.ru
Graduate Student; Dept. of Applied Mathematics & Computer Science 244, Molodogvardeyskaya st., Samara, 443100, Russian Federation

Elena V Baskinova

Samara State Technical University

Email: bashkinova-ev@yandex.ru
Cand. Phys. & Math. Sci., Associate Professor; Associate Professor; Dept. of Applied Mathematics & Computer Science 244, Molodogvardeyskaya st., Samara, 443100, Russian Federation

References

Supplementary files