CC BY
97
Вестн. Сам. гос. техн. ун-та. Сер. Физ.-мат. науки. 2017. Т. 21, № 3. С. 556-580 ISSN: 2310-7081 (online), 1991-8615 (print) d http://doi.org/10.14498/vsgtu1560
УДК 517.962.24+519.246
Численный метод оценки параметров нелинейного дифференциального оператора второго порядка
В. Е. Зотеев, Е. Д. Стукалова, Е. В. Башкинова
Самарский государственный технический университет,
Россия, 443100, Самара, ул. Молодогвардейская, 244.
Аннотация
Проблема нелинейного оценивания параметров систем различной физической природы, описываемых нелинейным дифференциальным оператором, является важнейшей проблемой математического моделирования. В статье рассматривается новый численный метод оценки параметров нелинейного дифференциального оператора второго порядка с дис-сипативной силой, пропорциональной n-ной степени скорости движения. В основе численного метода лежит среднеквадратичное оценивание коэффициентов обобщенной регрессионной модели, построенной с учетом разностных уравнений, описывающих результаты измерений импульсной характеристики системы. Реализованная в методе двухэтапная процедура дифференцированного оценивания параметров динамического процесса позволяет обеспечить высокую адекватность построенной модели данным эксперимента. Применение разработанного численного метода позволяет существенно (в несколько раз) повысить точность оценок параметров нелинейного дифференциального оператора по сравнению с известными методами за счет устранения смещения в оценках, обусловленного использованием аппроксимации при моделировании огибающей амплитуд колебаний.
Ключевые слова: нелинейный дифференциальный оператор, диссипа-тивная сила, разностные уравнения, обобщенная регрессионная модель, нелинейная регрессия, среднеквадратическое оценивание.
Получение: 2 августа 2017 г. / Исправление: 16 сентября 2017 г. / Принятие: 18 сентября 2017 г. / Публикация онлайн: 13 ноября 2017 г.
Статья
3 ©® Контент публикуется на условиях лицензии Creative Commons Attribution 4.0 International (https://creativecommons.org/licenses/by/4.0/deed.ru) Образец для цитирования
Зотеев В. Е., Стукалова Е. Д., Башкинова Е. В. Численный метод оценки параметров нелинейного дифференциального оператора второго порядка // Вестн. Сам. гос. техн. ун-та. Сер. Физ.-мат. науки, 2017. Т. 21, № 3. С. 556-580. doi: 10.14498/vsgtu1560. Сведения об авторах
Владимир Евгеньевич Зотеев Ä http://orcid.org/0000-0001-7114-4894 доктор технических наук, доцент; профессор; каф. прикладной математики и информатики; e-mail: [email protected]
Екатерина Дмитриевна Стукалова. http: //orcid. org/0000- 0002- 6064- 7725 магистрант; каф. прикладной математики и информатики; e-mail: [email protected]
Елена. Викторовна. Башкинова. http://orcid.org/0000-0001-6887-6162 кандидат физико-математических наук; доцент; каф. прикладной математики и информатики; e-mail: bashkinova- [email protected]
Введение. Одной из основных проблем, возникающих при исследовании объектов, явлений или систем различной физической природы, является нелинейность математической модели. В этих случаях достоверная оценка параметров модели по результатам эксперимента или промышленных испытаний представляет собой сложную вычислительную задачу. Обычно эта проблема решается на основе линеаризации модели тем или иным способом, однако возникающая при этом систематическая погрешность в результатах вычислений параметров может оказаться недопустимо большой. В ряде случаев, например, когда характеристики нелинейности исследуемой системы являются важнейшим диагностическим признаком ее технического состояния и линеаризация принципиально недопустима, приходится применять известные методы нелинейного оценивания, которые не всегда дают положительный результат.
В данной работе рассматривается нелинейный дифференциальный оператор второго порядка
Ь = шу"(Ь) + Ьг^у^)) + еу(Ь),
описывающий широкий класс систем с одной степенью свободы со слабой нелинейностью общего вида, где т, Ь и с — коэффициенты пропорциональности; функция г(у(£),у'(¿)) описывает нелинейность общего вида.
В частности, такой оператор широко используется при описании затухающих свободных колебаний различного рода механических систем с дисси-пативными силами, пропорциональными п-ной степени скорости движения (внутреннее трение в материале; конструкционное трение в опорах, шарнирах, сочленениях; силы сопротивления жидкой и газообразной среды; силы, возникающие с нагружением поглотителей энергии и т.п.):
Ь = ту" (*) + Ьу'(;)|у'(;)Г1 + су(1) = 0, (1)
где т — масса системы; су(Ь) —линейная восстанавливающая сила; слагаемое Ьу'(^)|у/(^)|га-1 описывает нелинейную диссипативную силу, пропорциональную п-ной степени скорости движения [1-4]; с и Ь — коэффициенты пропорциональности, п — параметр, характеризующий нелинейность диссипативной силы.
Одной из важнейших задач при построении математической модели технической системы с использованием этого нелинейного дифференциального оператора является задача достоверной оценки его параметров по результатам наблюдений реакции системы на некоторое типовое тестовое воздействие [5-7]. Широко используемый в практике эксперимента цифровой спектральный анализ и высокоэффективные методы вибродиагностики, а также методы корреляционного анализа [8-13] по своей природе не позволяют обеспечить требуемую точность и достоверность оценок характеристик нелинейности диссипативной силы, которые являются важнейшим диагностическим признаком технического состояния механической системы. В то же время известные и традиционно применяемые методы определения характеристик рассеяния энергии колебаний различных механических конструкций и демпфирующих свойств материалов уже не вписываются в формат современных информационных технологий [2, 4, 5]. В данной работе проблема параметрической идентификации нелинейного дифференциального оператора второ-
го порядка решается на основе среднеквадратичного оценивания параметров импульсной характеристики системы по результатам наблюдений, полученных в ходе научно-технического эксперимента.
Построение импульсной переходной характеристики методом
Ван-дер-Поля. В [14] отмечается, что импульсная переходная функция (характеристика) g(t) системы, описываемой нелинейным дифференциальным оператором, может быть получена из решения однородного дифференциального уравнения
my''(í) + bv(y(í),y'(í)) + cy(í) = 0
с начальными условиями вида y(0) = 0, y'(0) = 1/m.
Рассмотрим построение приближенного решения этого нелинейного дифференциального уравнения методом Ван-дер-Поля с учетом малости параметра b [15]. В соответствии с этим методом решение ищется в виде
y(í) = a(í) cos 0(í), 0(í) = wí + e(t), (2)
где a(í) и в(í) — некоторые функции, мало изменяющиеся за период колебаний T = 2n/w. При этом предполагается выполнение равенства [15]
y'(í) = —wa(í) sin 0(í), (3)
то есть, чтобы y(í) и y'(í) имели тот же вид, что и при отсутствии в системе слабой нелинейности общего вида bv(y(í), y'(í)), когда b = 0. Дифференцируя (2) и сравнивая результат с (3), получаем равенство
a'(í) cos 0(í) — а(£)в'(í) sin 0(í) = 0. (4)
Вычисляя вторую производную y''(í) из (3), подставляя полученный результат и функцию (2) в исходное дифференциальное уравнение, с учетом малости величин a'(t), в'(í) и b получаем два соотношения:
mw2a(í)cos 0(í) = c cos 0(í),
a'(í) sin 0(í) + a(í)e'(í)cos 0(í) = ——v(a(í)cos 0(í), —wa(í) sin 0(í)).
Из первого равенства следует, что w = д/c/m, а второе вместе с (4) образует систему уравнений, из решения которой выражаются функции a'(í) и a(í)e'(í):
a'(í) =-sin ф ■ via cos ф, — wa sin ф),
mw b (5)
a(í)e'(í) =-cos ф ■ via cos ф, — wa sin ф).
mw
В правых частях (5) и далее для краткости обозначено a = a(í), w = w(í). С учетом слабой нелинейности (малости параметра b) функции a(í) и в(í) и их производные за время периода колебаний T практически не изменяются. Усредняя (интегрируя по переменной í) на промежутке времени [0; T] левую и правую части равенств (5), с учетом dí = ^ф/w и соответствующей
замены промежутка интегрирования в правой части с [0; T] на [0; 2п], получаем систему уравнений, из решения которой можно найти функции a(t) и fi(t), входящие в уравнение (2), описывающее импульсную характеристику:
b ¡'2п
a'(t) =- r(acosф, —wasinф) sinфйф,
() 2птш J о ( ф, ф) ф ф (6)
b г2п Vй/
a(t)P'(t) = -- r(a cos ф, —ша sin ф)cos фйф.
2птш J о
В рассматриваемом нами частном случае при нелинейной диссипативной силе, пропорциональной n-ной степени скорости движения by'(t)jy'(y)|n-1, имеем:
r{y(t),y'(t)) = —aw sin ф| — aw sin(
n-1 = -anwn sin ф| sin ф|n-1.
При этом интегралы в (6) имеют вид
г2ж г2ж
/ r(a cos ф, — wa sin ф) sin ф^ф = —anwn j sin ф|п+1 0ф J о J о
f2n г2ж
/ r(a cos ф, —wa sin ф) cos ф^ф = —anwn sin ф| sin ф|п-1 cos оо
Очевидно, что в последнем выражении интегральная функция нечетная и имеет период равный 2п. Следовательно,
г- 2п
/ sin ф| sin ф|П-1 cos ф^ф = 0.
о
Обозначая
Jn = / | sin ф|п+1 #, о
с учетом последнего замечания систему (6) запишем в виде
bwn-1
a'(t) =— an(i)Jn; (7)
a(t)e'(t) = 0.
Обозначая a(0) = a0 и в(0) = во, из первого и второго уравнений системы (7) соответственно получаем
a(t) =-. ^ , в (t) = во.
bwn-1
i + (n — i) -^r^an-1Xnt
Отсюда с учетом (2) общее решение однородного нелинейного дифференциального уравнения (1) можно записать в виде
а0 cos(wt + во)
y(t) = , ......... /л ■ (8
1 + (" — 1) ^ an-Jnt
и
n-1
Отметим, что полученное методом Ван-дер-Поля уравнение (8) свободных колебаний нелинейной диссипативной системы с силой трения, пропорциональной п-ной степени скорости движения, полностью совпадает с результатом, полученным на основе метода энергетического баланса в работах [1-4,7].
С учетом начальных условий у(0) = 0, у'(0) = 1/т из (8) получаем частное решение, описывающее импульсную характеристику системы с нелинейным дифференциальным оператором второго порядка (1):
g(t) =
sin wt
wm
1 + (n - 1)
b
(9)
2nmn
. Jnt
Построенная математическая модель (9) импульсной характеристики системы с нелинейным дифференциальным оператором (1) содержит четыре параметра: m, n, b и w = л/c/m, которые известным образом связаны с параметрами дифференциального оператора, что позволяет свести задачу параметрической идентификации дифференциального оператора к оценке параметров импульсной характеристики. Величину Jn при различных значениях n можно найти из таблиц, приведенных в работах [1-4,7]. В частности, в табл. 1 представлены значения Jn для некоторых n.
Таблица 1
Значения интеграла Jn [The value of the integral Jn]
n 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4
Jn 4.00 3.77 3.58 3.42 3.27 3.14 3.03 2.92
n 1.6 1.8 2.0 2.2 2.4 2.6 2.8 3.0
Jn 2.83 2.75 2.67 2.60 2.53 2.47 2.41 2.37
При формировании выборки результатов измерений у^ свободных колебаний нелинейной диссипативной системы в дискретные моменты времени, на основе которой вычисляются оценки параметров, момент первого измерения ¿о, как правило, не совпадает с моментом £ = 0 приложения импульса силы к системе. Поэтому преобразуем построенную модель (9) для системы координат, начало которой совпадает с моментом времени ¿о первого отсчета в выборке результатов наблюдений:
у(£) = £(* + ¿о).
Заменяя в (9) переменную £ на £ + ¿о, после простых алгебраических преобразований получаем математическую модель в форме уравнения свободных колебаний нелинейной диссипативной системы с силой сопротивления, пропорциональной п-ной степени скорости движения:
y(t) =
ао cos(wt + 0о)
1 + (n - 1) 11
(10)
где ш = л/с/т — частота колебаний; Т = 2п/ш — период колебаний;
Ьш™
¿о = ¿(ао) =-«П-1 Л
n-1
n-1
— логарифмический декремент колебаний;
1 п
ао =---, ^о = ш£о - - (12)
- 1 / Ь 2
шш п д/1 + (п - 1)---Л¿о
V 2пш-
— амплитуда и фаза колебаний, соответствующие моменту времени ¿о первого отсчета в выборке результатов наблюдений относительно момента приложения импульса силы (¿о — время запаздывания).
Из соотношений (11) и (12) можно получить формулы, которые позволяют по оценкам параметров п, ¿о, ш и ао вычислить параметры ш, Ь и с для нелинейного дифференциального оператора (1):
ш = -Л2п - (п - 1)ДоШ*°, Ь = *°ш , с = шш2. (13) У 2па- ш--1 а- ш™-27„ У 7
Отметим, что формулы (13) предполагают известной величину ¿о — время запаздывания момента первого отсчета уо в выборке результатов наблюдений относительно момента времени приложения импульса силы. При неизвестной величине ¿о ее можно найти из формулы (12) только при условии ¿о € [0; Т).
Таким образом, задача параметрической идентификации нелинейного дифференциального оператора, описывающего диссипативную механическую систему с силами трения, пропорциональными п-ной степени скорости движения, может быть решена на основе оценки параметров свободных затухающих колебаний системы.
Анализ известных методов оценки параметров нелинейной дис-сипативной системы. Известные методы оценки динамических характеристик нелинейной диссипативной системы по результатам наблюдений ее свободных колебаний описаны в [3-5,7,10,12,16]. К основным недостаткам большинства из описанных методов можно отнести, во-первых, то, что практически все применяемые при экспериментальных исследованиях методы предполагают линеаризацию математической модели, описывающей в той или иной форме динамический процесс в диссипативной механической системе. Это является одним из источников систематической погрешности в оценках динамических характеристик и для систем с сильной нелинейностью вообще недопустимо. Во-вторых, большинство из известных методов за время одной реализации затухающих колебаний не предусматривает получения избыточного числа измерений, используемых для вычисления одной оценки динамических характеристик. Тем самым исключается возможность обеспечения помехозащищенности этих оценок за счет применения статистических методов обработки экспериментальных данных. В-третьих, алгоритмы вычислений предполагают достаточно громоздкие промежуточные графические построения или функциональные дополнительные преобразования колебаний с помощью различных технических устройств [4,10,12]. Во всех случаях это вносит дополнительные погрешности в результаты вычислений оценок динамических характеристик и не способствует эффективному применению вычислительной техники при экспериментальных исследованиях [7]. Устранить эти недостатки можно только на основе применения статистических методов
обработки результатов эксперимента на базе компьютеризации экспериментальных исследований.
Наиболее эффективно эта задача может быть решена на основе вычисления среднеквадратичных оценок параметров модели (10), минимизирующих величину отклонения этой модели от данных эксперимента в евклидовом нормированном пространстве:
где у — Ж-мерный вектор результатов эксперимента, у — Ж-мерный вектор результатов расчета по модели (10). В такой постановке задачу параметрической идентификации дифференциального оператора можно интерпретировать как задачу нелинейного оценивания (нелинейной регрессии). Однако результаты научных исследований показали, что известные методы нелинейного оценивания: метод Ньютона—Гаусса, метод Хартли, метод Левенберга— Маркуардта, а также такие общие методы минимизации функций, как градиентный и сопряженных градиентов [17-20], не позволяют находить устойчивые оценки параметров модели (10), обеспечивающие минимум квадрата нормы ее отклонения от данных эксперимента. Основными причинами этого являются проблема выбора вектора начальных оценок параметров и требование удовлетворения на каждой итерации уточнения параметров ограничениям, накладываемым видом функциональной зависимости (10).
В работе [7] эта проблема решалась на основе аппроксимации модели (10) функциональной зависимостью вида
На основе этой зависимости были построены разностные уравнения, описывающие результаты эксперимента, и получена линейная регрессионная модель, коэффициенты которой известным образом связаны с параметрами свободных колебаний диссипативной системы [7,22]. Такой подход позволил обеспечить высокую адекватность аппроксимативной модели результатам наблюдений, и, как следствие, достоверные оценки ее параметров.
Однако аппроксимация модели (10) неизбежно приводит к методической систематической погрешности в результатах расчета параметров режима свободных колебаний диссипативной системы, а значит, и в оценках параметров нелинейного дифференциального оператора (1). На основе численно-аналитического моделирования проведены исследования погрешностей, обусловленных аппроксимацией модели (10). Численный эксперимент проводился следующим образом. При времени наблюдения ¿я = 10 и периоде дискретизации т = 0.2 по формуле
где ао = 1, Т = 1, к = 0,1, 2,..., N — 1, формировались выборки результатов расчета значений огибающей амплитуд колебаний (объем выборок составил
N-1
1|у — у11 = (Ук — ук)2 ^ т^
ао сов^ + ^о)
N = 51). Были установлены следующие значения уровней затухания колебаний ¿а = a(tH)/а0 за время tH = 10: ¿а = 0.1; 0.2; 0.25; 0.33; 0.5; 0.8.
Исследования проводились для 17-ти значений параметра нелинейности n в диапазоне от 0 до 2.5. Величина декремента колебаний в уравнении огибающей амплитуд колебаний рассчитывалась по формуле
¿с =
1 - а"-1
(1 - ап-1)(п - 1)tn
таким образом, чтобы обеспечить при заданном значении п требуемую величину уровня затухания ¿а.
Для оценки погрешности вычисления параметров n и ¿о использовалась аппроксимация уравнения огибающей амплитуд колебаний в виде
a(t) =---а0-—-г-, ао = 1, T = 1. (14)
i+|t+(i - n)(f) t2
Оценки параметров модели (14) находились из условия минимизации остаточной суммы квадратов
N -1
11а — а|| = ^ (ак — ak)2 ^ min,
k=0
где N = 51, ak = a(tk), k = 0,1,..., N — 1. Величина Да = ||а — а||/||а||, характеризующая среднеквадратическое отклонение модели огибающей амплитуд колебаний от ее аппроксимации (14), в диапазоне изменения n £ [0.0; 2.5] и ¿а £ [0.1; 0.8] не превышала 0.05.
Результаты исследований в виде кривых зависимостей относительных погрешностей Дп = |(n — П)/п| и Д¿0 = ^¿0 — ¿(^/¿о| представлены на рис. 1. Очевидно, что только в относительно узком диапазоне изменения параметра n от 1.3 до 2.1 погрешности Дп и Д¿о могут считаться приемлемыми. Причем с уменьшением уровня ¿а ^ 0.33 погрешности резко возрастают.
Расчет величины отношения числа результатов вычисления параметров n и ¿о, относительная погрешность которых не превышает 0.05, к общему числу 150 результатов вычислений (при п £ [0.1; 2.5] и ¿а £ [0.1; 0.8]), сделанный на основе таблицы данных исследования, показал следующие результаты. Для параметра п эта величина составляет всего 0.3, а для параметра ¿о — 0.76. Таким образом, можно сделать вывод о том, что область применения численного метода оценки параметров режима свободных колебаний механической системы, предложенного в [7], существенно ограничена величиной и степенью нелинейности (параметры b и п) сил трения, действующих в диссипативной системе.
Численный метод оценки параметров на основе разностных уравнений. В данной работе предлагается новый численный метод оценки параметров математической модели (10) по результатам эксперимента. Этот метод состоит из двух этапов. На первом этапе модель (10) аппроксимируется функциональной зависимостью вида
со cos(wt + ^о)
1 + Cit + C2t2
п п
—t— 6a = 0.1 -i- 6a = 0.2 —6a = 0.25 6a = 0.33 —•- 6a = 0.5 —•— 6a = 0.8
Рис. 1. Зависимости относительной погрешности An и Д50 от параметра нелинейности n дифференциального оператора при различных уровнях затухания колебаний 5а за промежуток времени наблюдения tH [Figure 1. Dependences of a relative error An and A50 from the parameter of nonlinearity of n the differential operator at various levels of attenuation of fluctuations 5а for a period of
observation of tH ]
дискретные значения которой в точках tk = Tk, k = 0,1, 2,..., N — 1, описываются формулой
с = Со cos(^rk + )
" = 1+ CiTk + С2Т2k2 . (ibj
Здесь же по найденным среднеквадратичным оценкам "о, " и ¿2 параметров этой модели вычисляются значения огибающей амплитуд колебаний:
"k =1+ CiZ- C2T 2 k2 • k = 0 1' 2---"N — 1- (16)
На втором этапе с учетом сформированной выборки результатов вычислений "k по формуле (16) на основе минимизации функционала
N -1
||" — "||2 = ^ ("k — "k)2 ^ min, k=0
где
"о
"k = -
1 + (n - 1) ^тк 2п
находятся среднеквадратичные оценки параметров п и ¿о. Для этого при формировании линейной регрессионной модели, коэффициенты которой известным образом связаны с параметрами нелинейной зависимости п и ¿о, используется дифференцирование зависимостей, описывающих огибающую амплитуд колебаний а^) в модели (10) и ее аппроксимацию а^).
п-1
На первом этапе численного метода вначале на основе дискретной модели (15) строится система разностных уравнений, описывающих результаты эксперимента у^, к = 0,1, 2,..., N — 1, где N — объем выборки результатов наблюдений с периодом дискретизации т. Из формулы (15) можно получить аналогичные соотношения для и ук_2, и на их основе с помощью простых алгебраических преобразований построить систему разностных уравнений вида
уо = А4;
у1 = А5;
< У/с-2 + Ш = А1Ук-1 + А2[(к — 2)Ук-2 — А1(к — 1)^-1 + кШ ] + +Аз[(к — 2)2ук-2 — А1(к — 1)2ук-1 + к2ук ],
к = 2,3,..., N — 1,
коэффициенты которой связаны с параметрами модели (15) соотношениями
Ai = 2 cos wt, А2 = —cir, A3 = —c2r2,
А4 = со cos ^о, А5 = --г1-— (А4 cos wt — со sin ^o sin wt).
1 — A2 — A3
С учетом случайного разброса в результатах эксперимента относительно модели yk:
Ук = yAk + efc
построенную систему уравнений можно привести к виду
Уо = А4 + ео; yi = А5 + ео;
Ук-2 + Ук = А1Ук-1 + А2[(к — 2)yfc-2 — А^к — 1)yfc_i + kyfc ] +
+Аз[(к — 2)2yfc-2 — А^к — 1)2 ук-i + k2yfc ] + Пк+i, (17)
k = 2,3,...,N — 1; 1 ;
Пк+i = [1 — А2(к — 2) — Аз(к — 2)2]efc-2 —
—А? [1 — А2(к — 1) — Аз (к — 1)2]ек-1 + [1 — А2 к — Азк2]ек,
к = 2,3,..., N — 1,
а (г)
где величина а( при i = 0,1, 2,... является некоторой известной оценкой коэффициента А(.
При такой интерпретации на основе системы разностных уравнений (17) можно построить линейную обобщенную регрессионную модель вида
b = А + П; (18)
элементы которой описываются следующими соотношениями: А = (А(, А2, А3, А4, А5)т — вектор неизвестных коэффициентов; b = (Уо, yi, Уо + y2, yi + Уз, ..., yw-з + yw-i)T, — вектор свободных членов; е = (ео, е1, е2,..., ew-()Т — вектор случайной помехи в результатах эксперимента;
=
0 0
у1 у2 уз
0 0
—а14)у1 + 2у2
у 1 — А1 )2у2 + 3уз у1 — А1) 4у2 + 9уз
0 0
—А(/)у1 + 4у2 1( )
2у2 — А1) 3уз + 4у4 4у2 — АС) 9уз + 16у4
( )
ум-2 (Ж — 3)ум-з —
—а1)(Ж — 2)ум-2+ + — 1)ум-1
— матрица регрессоров; П = (П1,П2,Пз,...,Пк, ...,nN )Т
вектор невязки;
(Ж — 3)2ум-з— 0 0 —а1)(Ж — 2)2ум-2 + + (Ж — 1)2ум-1
Р,
А«
0 1
—а1)(1 — А« — А«) 1 — А« — А« 0
0 0
1 — А« — А« —А1)(1 — 2А20 — 4А30)
0 0
0 к(<)
«
1 — 2а24) — 4А:
( )
1 — 3А2') — 9А3
—а1)(1 — 3А2° — 9А30)
1 — (Ж — 1)А2г ) — (Ж — 1)2А3
(19)
— матрица линейного преобразования вектора случайной помехи в результатах эксперимента.
Применяя алгоритм численного метода нелинейного оценивания на основе разностных уравнений [7,21-23], с учетом второго уравнения в системе (18) преобразуем первое уравнение к виду
Р-1Ь = Р-1 **«)А +
(О
(20)
Поскольку величина в левой части и в матрице уравнения регрессии (19) известна, это уравнение линейно относительно вектора коэффициентов А. Следовательно, среднеквадратичную оценку А этого вектора, минимизирующего сумму квадратов величин разброса данных эксперимента относительно модели (15):
^ 1
У^ е2к ^ тт,
к=о
можно найти из решения нормальной системы линейных алгебраических уравнений
^л« А = р-« ^ь,
0
0
0
0
где QA(i) = P^) PX(i). Отсюда получаем
= Ко QÄ(1) FÄ(i) ]"1fJ) ^(1) b i = 0,1, 2,... (21)
Формула (20) описывает итерационный процесс уточнения оценки вектора коэффициентов обобщенной регрессионной модели (18). В качестве начального приближения оценки коэффициента Ai можно выбрать либо величину
i (о) = Efc1 (yfc-2 + yfc )yfc-i
Al = TN-1 y2 ,
которая минимизирует сумму квадратов отклонений N -1
У^ [yk — c0 cos(wrk + ^0)]2 ^ min, fc=0
либо первый элемент в оценке вектора А, полученный из уравнения регрессии вида
yfc-1 + yfc = A1yfc_1 + A2 [(k — 2)yfc_2 + kyfc] + A3 [(k — 2)2yfc-2 + k2yfc] +
+ A4(k — 1)yfc-1 + A5(k — 1)2 yfc-1 + , k = 2,...,N — 1,
N 1 2
при минимизации остаточной суммы квадратов fc=2 П2 ^ min.
Вычисления по формуле (21) следует продолжать до выполнения условия
||A(i+1)-A(i)||
= maxiAj+1 — AjI < 0, (22)
и и ^=1,3 j ji
где 0 — произвольно малая положительная величина. Полученные при выполнении условия (22) оценки коэффициентов обобщенной регрессионной модели (18) следует принять за оценки Aj, j = 1, 2, 3, 4, 5, минимизирующие величину среднеквадратического отклонения модели (15) от данных эксперимента.
По найденным среднеквадратическим оценкам коэффициентов Aj с учетом соотношений (16) оценки параметров модели (15) можно найти по формулам
1 A1 „ A2 „ A3
ca = - arccos —, a =--, c2 = —„,
т 2 т т2
A1A4 — 2 A5(1 — A2 — Аз) u =-
У4—12
„ /i2 . 2 [ A1A4 — 2 A5(1 — A2 — A3)]2 . A2
Ао = V A4 + u2 = W-43^-+
U А
arctg —, при a4 > 0, A4
u
tA0 = arctg A—+ sign(u) ■ п, при A4 < 0,
A4
sign(u) —, при A4 = 0.
Завершается первый этап вычислением значений аппроксимации огибающей амплитуд колебаний в (15) по формуле (16):
со
ан =
1 + C1 tk + c2t 2k2'
Заметим, что по результатам вычислений на первом этапе можно найти приближенную оценку логарифмического декремента колебаний ¿о в (10), соответствующего начальной амплитуде колебаний ао:
г 2п(С1^2 + 2ПС2) 4п2д2 +2пЛ2 агссо8 у
¿0 = ""
w2 + 2nciw + 4п2С2 , 2Х . 0 , Ai 2 Ai'
4п2 A2 + 2nA2 arceos —— arceos2 —
На втором этапе предлагаемого численного метода находятся оценки параметров нелинейности системы n и ¿о. В основе вычисления оценок лежит минимизация функционала
N -1
22
|а — а||2 = ^ (ak — ak)2 ^ min,
k=0
где элементы вектора а находятся по формуле
Сн =-■ а° , . (23)
'1 + (п - 1) ¿Птк
Очевидно, что задача среднеквадратического оценивания параметров модели (23) по результатам расчета по формуле (16) относится к достаточно сложной задаче нелинейной регрессии [17-20]. В данной работе предлагается перейти к линейной регрессионной модели, коэффициенты которой известным образом связаны с параметрами п и ¿о. С этой целью предлагается вначале провести дифференцирование непрерывных в области своего определения функций:
/ ,ч со _ ао
а(£) = -----—2 , а(^) =
1 + с1 £ + о,*2' п-^^ ,/осс '
V1 + (п -^
Находя производные, получаем
¿ой
а'Ю = _ Со(С1 +2С2^) 2, а'(<) =--а(<).
() (1 + + с2^2)2' () 1 + (п -1) ¿оа
2п
Очевидно, что имеет место равенство
а(£) = а(£) + е(£),
где непрерывная функция е(£) описывает ошибку аппроксимации зависимости а(£) построенной моделью а(£). Дифференцируя это равенство, получаем а'(£) = а'(¿) + е'(£). Теперь, переходя в область дискретного изменения аргумента ¿: = тк, к = 0,1, 2,..., N— 1, с учетом полученных выше соотношений имеем
ай =--2п А Л—(а& — ) + 4,
1 + (п — 1) ^ тк
где значения а^, к = 0,1, 2,..., N—1, вычисляются по формуле (16), а зачения а'д, по формуле
^ = (1 У* 2)2 , к = 0,1, 2,..., N — 1.
й (1 + С1тк + с2т 2к2)2
Отсюда после простых алгебраических преобразований, применяя формулы численного дифференцирования второго порядка
4 = е*+ — ей-1, к = 1, 2,..., N — 2,
Й 2т
и
/ _ ем-з — 4ем-2 + Зем-1
ем-1= 2Т ,
получаем систему уравнений вида
( а0 = Аз + ео;
а/ а/ к Л аЛ 1 + кД1 е + Л е +1 + Ы1 е ай = — ай к Л1 — ай Л 2--^-ей-1 + Л2ей +--2Т-
к = 1, 2,..., N — 2; п' п' 1) Л а Л е
1 = — — 1) Л1 — ам-1Л2 +----ем-з —
2т
(24)
2
— [1 + ^ — 1)Л1] ем-2 + Л 2 + — [1 + ^ — 1)Л1]
т
З
2т
коэффициенты которой связаны с параметрами зависимости а(£) формулами Л1 = (п — 1) М т, Л 2 = ^, Лз = ао. (25)
Соответствующая системе уравнений (24) обобщенная регрессионная модель принимает вид
'Ь = РЛ + П, (26)
элементы которой описываются следующими соотношениями: Л = (Л1, Л 2, Л3)т — вектор неизвестных коэффициентов;
е = (е0,е1 , е2,...,ем-1)т — вектор отклонений модели а(£) от результатов расчета а(£); п = (Пъ П2,..., )Т — вектор невязки;
Ь = (ао,а/1,а/2,...,а'м-1)т (27)
— вектор свободных членов;
0 0 1
- a1 -a1 0
-2a2 -a2 0
F = 3a3 -a3 0
-(N - 1)aN-1 -aN-1 0
матрица регрессоров; 1
1 + A(1i) 2т А2 1 + A« 2т 0
0 1 + 2A1i) 2т A2 1 + 2A 2т
P\ = 0 0 1 + 3A1i) 2т A2
(i) 1
0 0 0
0 0 0
0 0 0
0 0 0
0 0 0 0
(i)
1 + (N - 3)А1'
2т
A(i) А2
1 + (N - 2)А
1 + (N - 2)А1
(i)
A2i) +
2т
1 + (N - 2)A
(28)
2т
(29)
— матрица линейного преобразования вектора е отклонения момента от результатов расчета ak.
Применяя алгоритм численного метода нелинейного оценивания на основе разностных уравнений [7,21-23] с целью минимизации суммы квадратов отклонения модели (23) от результатов расчета по формуле (16):
N-1
У^ е| ^ min, k=0
можно получить рекуррентную формулу
Ai+1 = [F^ П-1 FX(i)Q-1)b i = 0 1, 2,...
(30)
0
0
0
0
3
где матрица QA(i) = PA'(i) PA(i) находится с учетом формулы (29), а вектор b и матрица F определяются формулами (27) и (28) соответственно. Формула (30) описывает итерационный процесс уточнения оценки вектора коэффициентов регрессионной модели (26). Начальное приближение вектора оценок Л(о) можно найти из условия минимизации функционала
||П||2 = ||b - F Л||2 ^ min
из соотношения
Л(0) = (FTF)-1FT b. Вычисления по формуле (30) следует продолжать до выполнения условия (22):
Ц^1) - ЛМ|| = max|V+1 - aj| < 0,
и и j ji
где 0 — заданная малая величина.
По найденным среднеквадратическим оценкам коэффициентов с учетом соотношений (25) определяются оценки параметров модели (23):
Л 2п л Л1 Л Л
öo = — Л2, n = ---+1, ао = Лз.
W Л2Т
Второй этап численного метода завершается вычислением оценок параметров нелинейного дифференциального оператора (1) по формулам
n-i/2п - (П - 1)Ло<Лto , Л от /оп т = \ -Цг—:—--, b = ^^-, c = mw2. (31)
V 2nan-1wn -1 an -1wn-2J«
Результаты численно-аналитических исследований. Для проверки эффективности предлагаемого численного метода проведены численно-аналитические исследования, а также выполнен сравнительный анализ точности результатов вычислений, полученных этим методом и численным методом, описанным в [7,21,22]. Имитационное компьютерное моделирование и реализация алгоритмов этих методов проводились при следующих параметрах нелинейного дифференциального оператора (1): т = 0.4, Ь = 0.2, п = 3 и с = 12. Такие значения параметров соответствуют достаточно высокой степени нелинейности и величине диссипативной силы, обуславливающей затухание свободных колебаний в системе. В соответствии с формулой (10) с периодом дискретизации т = 0.05 формировались выборки результатов наблюдений уд, объем которых был равен N = 40. К результатам моделирования уд добавлялась случайная некоррелированная помеха ед, имеющая нормальный закон распределения с математическим ожиданием М[ед] =0 и дисперсией ] = о2 ед € N(0, о^), которая генерировалась таким образом, чтобы ее величина удовлетворяла равенству
Y
\
-1 s2
-1 a2 ;
где значения 7 выбирались из диапазона от 0 до 0.05 равномерно с шагом 0.005. По результатам суммирования = с + ек, к = 0,1, 2,..., N — 1, в соответствии с алгоритмами двух различных методов вычислялись оценки т, Ь, п и с искомых параметров. При одном и том же значении мощности помехи 71 формирование выборки результатов наблюдений повторялось 100 раз. Затем для каждого значения ^г из диапазона от 0 до 0.05 вычислялись среднеквадратичные оценки погрешности.
В качестве оценки погрешности был выбран второй центральный момент оценки параметра дифференциального оператора р относительно его точного значения р:
М [(р — р)2] = £[р] + (М [р] — р)2.
Такая оценка погрешности наиболее информативна, так как она интегрировано включает в себя дисперсию оценки параметра
£[р] = М [(р — М [р])2]
и квадрат смещения математического ожидания этой оценки М[р] относительно его точного значения. При приближенном вычислении математического ожидания М[р] использовалось среднее арифметическое результатов расчета параметра в 100 реализациях при одном и том же значении ^г мощности случайной помехи:
1 100 М [р]^ щ £ р.
.7 = 1
Результаты исследований в форме зависимостей относительной суммы квадратов ошибки
1 /^-^100 Ар = 7=1 (р7 — р )2,
где р любой из параметров т, с, Ь, п, от величины случайной помехи в данных эксперимента, представлены на рис. 2.
По результатам исследований можно сделать следующие выводы. Погрешности оценок параметров следует анализировать дифференцированно в зависимости от этапа численного метода. На первом этапе вычисляются параметры, связанные с частотой колебаний ш, а также с характеристикой затухания этих колебаний ¿0, но без учета вида модели Ьу7(^)|у'(^)|п 1 нелинейной зависимости диссипативной силы от скорости движения. Из кривых, представленных на рис. 2 (сверху), следует, что суммы квадратов ошибок Ат и Ас, во-первых, существенно зависят от величины случайной помехи в результатах наблюдений, и, во-вторых, с увеличением этой величины кривые 1 и 2 практически совпадают. Это объясняется тем, что алгоритмы вычислений оценок на первом этапе предлагаемого метода и известным методом [7,21,22] похожи и дают результаты с близкими дисперсиями. При этом проведенные исследования показали, что доминирующей составляющей в суммарной среднеквадратичной ошибке М[(р — р)2], р е {т, с}, является дисперсия ^[р], а квадрат смещения (М[р] —р)2 математического ожидания оценки параметра относительно точного значения практически не зависит от е и не превышает величины 0.26 для известного метода и 0.07 для нового численного метода от общей суммы квадратов ошибок (при 7 = 0.05).
0.08
0.06
% 0.04
0.02
1
/
2
0.01 0.02 0.03 0.04 0.05
7
0.6 0.5 0.4 < 0.3 0.2 0.1
0.25
0.20
0.15
е <1
0.10
0.05
___i——Т-
0.01
0.02 0.03
7
0.04
0.05
0.01
0.02 0.03
7
0.04
0.05
Рис. 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 невелики и имеют порядок случайной помехи в результатах наблюдений.
Иной характер имеют зависимости суммарной среднеквадратичной погрешности Ab и An оценок параметров. В известном методе эти оценки вычисляются непосредственно по формулам, построенным при разработке алгоритма метода [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], которая составляет величину, не
превышающую 0.02 от величины общей суммы квадратов ошибок. В алгоритме предлагаемого численного метода оценки параметров b и n находятся на втором этапе из условия минимизации величины среднеквадратическо-го отклонения аппроксимации (16), построенной на предыдущем этапе, от дискретной модели (23). В новом предложенном численном методе, как показали результаты исследований, смещение M[p] — p мало, соизмеримо со среднеквадратическим отклонением a = \JD[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, затем при выбранном времени запаздывания to = 0.1 по формулам (12) вычислялись амплитуда ао и фаза колебаний соответствующие моменту времени to, и в соответствии с (11) рассчитывалась величина декремента колебаний ¿о. Результаты этих расчетов представлены в последних столбцах второй строки табл. 2.
Таким образом, имитационная модель колебаний диссипативной системы имела вид
. , 0.411 cos(5.477t — 1.023) . .
y(t) =-Л --. (32)
yv 7 V1 + 1.889t
По формуле (32) при tk = тк, k = 0,1,2,...,N — 1, с периодом дис-
Таблица 2
Результаты вычислений оценок параметров [The results of calculations of the parameter estimates]
Parameter Differential operator Impulse characteristic
m b n c do ш n ao ^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
кретизации т = 0.05 были вычислены значений дискретной функции уь, tk = тк, к = 0,1, 2,..., 39. К этим значениям была добавлена случайная помеха (некоррелированные случайные величины еь € N(0, а^) с нулевым математическим ожиданием и одинаковыми дисперсиями, распределенные по нормальному закону), величина которой была сгенерирована таким образом, чтобы отклонение по евклидовой норме в относительных единицах
Де = IIу - У11
составило величину Де = 0.2. Полученные таким образом величины
Ук = Ук + ек, к = 0,1, 2,...,Ж - 1, (33)
использовались далее в качестве смоделированной выборки результатов эксперимента. Затем на основе описанного выше алгоритма численного метода были найдены оценки параметров математической модели (10), которые представлены в последних столбцах третьей строки табл. 2. В результате была построена математическая модель колебаний диссипативной системы вида
, 0.424cos(5.445t - 0.992) .
уМ = -с°4у/1 + 1.742? ) ■ (34)
На рис. 3 изображены результаты апробации метода на основе компьютерного моделирования: маркеры 1 соответствуют смоделированным результатам эксперимента уь со случайной помехой Де = 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))]
кривая 2, построенная по формуле (32), описывает смоделированный «точный» динамический процесс у(£); кривая 3 описывает динамический процесс у(¿), соответствующий построенной модели (34).
Полученные данные показывают, что даже при достаточно высоком уровне случайной помехи Де = 0.2 восстановленная по результатам эксперимента кривая 3, описывающая динамический процесс в системе, практически совпадает с «точной» кривой 2, смоделированной при заданных значениях параметров нелинейного дифференциального оператора (1).
Количественно это расхождение между построенной моделью (34) и имитационной «точной» моделью (32) можно оценить, используя евклидову норму в векторном пространстве:
Очевидно, что по сравнению с величиной Де = 0.2 расхождение между данными эксперимента и результатами расчета по модели (34) незначительно.
Найденные по формулам (31) оценки параметров нелинейного дифференциального оператора (1) приведены в первых столбцах третьей строки табл. 2. В последней строке табл. 2 представлены относительные погрешности результатов вычислений параметров нелинейного дифференциального оператора и параметров построенной модели (34) динамического процесса в системе, описываемой с помощью этого дифференциального оператора.
Заключение. По результатам проведенных исследований численного метода оценки параметров нелинейного дифференциального оператора можно сделать следующие выводы.
Разработан новый численный метод оценки параметров нелинейного дифференциального оператора с диссипативной силой, пропорциональной п-ной степени скорости движения. В основе метода лежит среднеквадратичное оценивание коэффициентов обобщенной регрессионной модели, построенной с учетом разностных уравнений, описывающих результаты измерений импульсной характеристики системы.
Реализованная в методе двухэтапная процедура дифференцированного оценивания параметров динамического процесса позволяет обеспечить высокую адекватность построенной модели данным эксперимента за счет использования уравнения скорости затухания колебаний при разработке алгоритма вычислений на втором этапе метода.
Применение разработанного численного метода позволяет существенно (в несколько раз) повысить точность оценок параметров нелинейного дифференциального оператора по сравнению с известными методами за счет устранения смещения в оценках, обусловленного применением аппроксимации при моделировании огибающей амплитуд колебаний.
Конкурирующие интересы. Мы не имеем конкурирующих интересов.
Авторский вклад и ответственность. Все авторы принимали участие в разработке концепции статьи и в написании рукописи. Авторы несут полную ответственность за предоставление окончательной рукописи в печать. Окончательная версия рукописи была одобрена всеми авторами.
Финансирование. Работа выполнена при поддержке РФФИ (проект № 16—01— 00249_а).
0.051.
НУ I
Библиографический список
1. Божко А. Е., Голуб Н. М. Динамико-энергетические связи колебательных систем. Киев: Наук. думка, 1980. 188 с.
2. Пановко А. Г. Основы прикладной теории колебаний и удара. Л.: Машиностроение, 1976. 320 с.
3. Пановко Я. Г. Внутреннее трение при колебаниях упругих систем. М.: Физматгиз, 1960. 194 с.
4. Писаренко Г. С., Яковлев А. П., Матвеев В. В. Вибропоглощающие свойства конструкционных материалов: Справочник. Киев: Наук. думка, 1971. 376 с.
5. Марков С. И., Минаев В. М., Артамонов Б. И. Идентификация колебательных систем автоматического регулирования/ Библиотека по автоматике. Выпуск 543. Л.: Энергия, 1975. 96 с.
6. Штейнберг Ш. Е. Идентификация в системах управления / Библиотека по автоматике. Выпуск 668. М.: Энергоатомиздат, 1987. 80 с.
7. Зотеев В. Е. Параметрическая идентификация диссипативных механических систем на основе 'разностных уравнений. М.: Машиностроение, 2009. 344 с.
8. Bendat J., Piersol A. Engineering Application of Correlation and Spectral Analysis. New York: Wiley-Interscience, 1980. xiv+302 pp.
9. Явленский К. Н., Явленский А. К. Вибродиагностика и прогнозирование качества механических систем. Л.: Машиностроение, 1983. 239 с.
10. Вибрации в технике: Справочник. В 6 т. М.: Машиностроение; Т. 1, 1978. 352 с.; Т. 2, 1979. 351 с.; Т. 5, 1981. 496 с.
11. Jenkins G. M., Watts D. G. Spectral analysis and its applications / Holden-Day series in time series analysis. San Francisco: Holden-Day, 1968. xviii+525 pp.
12. Добрынин С. А., Фельдман М. С., Фирсов Г. И. Методы автоматизированного исследования вибраций машин: Справочник. М.: Машиностроение, 1987. 224 с.
13. Marple S. Lawrence, Jr. Digital Spectral Analysis: With Applications / Prentice-Hall Series in Signal Processing. New York: Prentice-Hall, 1987. xx+492 pp.
14. Деч Г. Руководство к практическому применению преобразования Лапласа и z-преоб-разования. М.: Наука, 1971. 288 с.
15. Nayfeh A. H. Introduction to perturbation techniques. New York: John Wiley & Sons, 1993. xiv+519 pp.
16. Писаренко Г. С., Матвеев В. А., Яковлев А. П. Методы определения характеристик колебаний упругих систем. Киев: Наук. думка, 1976. 88 с.
17. Draper N. R., Smith H. Applied Regression Analysis / Wiley Series in Probability and Statistics. New York: John Wiley & Sons, 1998. xix+716. doi: 10.1002/9781118625590
18. Демиденко Е. З. Линейная и нелинейная регрессии. М.: Финансы и статистика, 1981. 302 с.
19. Marquardt D. W. An Algorithm for Least-Squares Estimation of Nonlinear Parameters // J. Soc. Indust. Appl. Math., 1963. vol.11, no. 2. pp. 431-441. doi: 10.1137/0111030.
20. Hartley H. O., Booker A. Nonlinear Least Squares Estimation// Ann. Math. Statist, 1965. vol.36, no. 2. pp. 638-650. doi: 10.1214/aoms/1177700171.
21. Зотеев В. Е. Построение разностных уравнений для повышения точности параметрической идентификации колебательных систем со слабой нелинейностью общего вида// Вестн. Сам. гос. техн. ун-та. Сер. Физ.-мат. науки, 2000. №9. С. 169-173. doi:10.14498/vsgtu44.
22. Зотеев В. Е. Разработка и исследование линейных дискретных моделей колебаний дис-сипативных систем // Вестн. Сам. гос. техн. ун-та. Сер. Физ.-мат. науки, 1999. №7. С. 170-177. doi: 10.14498/vsgtu222.
23. Зотеев В. Е. Параметрическая идентификация линейной динамической системы на основе стохастических разностных уравнений // Матем. моделирование, 2008. Т. 20, № 9. С. 120-128.
Vestn. Samar. Gos. Tekhn. Univ., Ser. Fiz.-Mat. Nauki
[J. Samara State Tech. Univ., Ser. Phys. Math. Sci.], 2017, vol. 21, no. 3, pp. 556-580
d http://doi.org/10.14498/vsgtu1560
ISSN: 2310-7081 (online), 1991-8615 (print)
MSC: 65C20, 65P40
Numerical method of estimation of parameters
of the nonlinear differential operator of the second order
V. E. Zoteev, E. D. Stukalova, E. V. Bashkinova
Samara State Technical University,
244, Molodogvardeyskaya st., Samara, 443100, Russian Federation.
Abstract
The main problem of mathematical simulation is the problem of nonlinear estimation of parameters of the different physical systems. The article contains new numerical method of parameters estimation of the nonlinear differential operator of the second order with the dissipative force, proportional to n-motion speed level assessment. Mean square estimation of coefficients of the generalized regression model constructed taking into account the difference equations describing results of measurements of a pulse response of system is the cornerstone of the numerical method. Two landmark procedure of differentiated estimation of parameters of dynamic process realized in a method allow to provide high adequacy of the constructed model to data of an experiment. Application of the developed numerical method allows to increase significantly (several times) the accuracy of estimates of parameters of the nonlinear differential operator in comparison with the known methods due to elimination of the offset in estimates caused by use of approximation in case of simulation of an envelope of vibration amplitudes.
Keywords: nonlinear differential operator, dissipative force proportional, difference equations, generalized regression model, nonlinear regression, mean square estimation.
Received: 2nd August, 2017 / Revised: 16th September, 2017 / Accepted: 18th September, 2017 / First online: 13th November, 2017
Article
3 ©® The content is published under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/) Please cite this article in press as:
Zoteev V. E., Stukalova E. D., Bashkinova E. V. Numerical method of estimation of parameters of the nonlinear differential operator of the second order, Vestn. Samar. Gos. Tekhn. Univ., Ser. Fiz.-Mat. Nauki [J. Samara State Tech. Univ., Ser. Phys. Math. Sci.], 2017, vol. 21, no. 3, pp. 556-580. doi: 10.14498/vsgtu1560 (In Russian). Authors' Details:
Vladimir E. Zoteev http://orcid.org/0000-0001-7114-4894
Dr. Tech. Sci.; Professor; Dept. of Applied Mathematics & Computer Science;
e-mail: [email protected]
Ekaterina D. Stukalova http://orcid.org/0000-0002-6064-7725 Graduate Student; Dept. of Applied Mathematics & Computer Science; e-mail: [email protected]
Elena V. Baskinova http://orcid.org/0000-0001-6887-6162
Cand. Phys. & Math. Sci., Associate Professor; Associate Professor; Dept. of Applied Mathematics & Computer Science; e-mail:[email protected]
Competing interests. We have no competing interests.
Authors' contributions and responsibilities. Each author has participated in the article concept development and in the manuscript writing. The authors are absolutely responsible for submitting the final manuscript in print. Each author has approved the final version of manuscript.
Funding. This work was supported by the Russian Foundation for Basic Research (project no. 16-01-00249_a).
References
1. Bozhko A. E., Golub N. M. Dinamiko-energeticheskie sviazi kolebatel'nykh sistem [Dynamic-energy relations of oscillatory systems]. Kiev, Nauk. dumka, 1980, 188 pp. (In Russian)
2. Panovko A. G. Osnovy prikladnoi teorii kolebanii i udara [Fundamentals of applied theory of vibrations and shock]. Leningrad, Mashinostroenie, 1976, 320 pp. (In Russian)
3. Panovko Ya. G. Vnutrennee trenie pri kolebaniiakh uprugikh sistem [Internal Friction in Vibrations of Elastic Systems]. Moscow, Fizmatgiz, 1960, 194 pp. (In Russian)
4. Pisarenko G. S., Yakovlev A. P., Matveev V. V. Vibropogloshchaiushchie svoistva konstrukt-sionnykh materialov: Spravochnik [Vibration-Absorbing Properties of Structural Materials: A Handbook]. Kiev, Nauk. dumka, 1971, 376 pp. (In Russian)
5. Markov S. I., Minaev V. M., Artamonov B. I. Identifikatsiia kolebatel'nykh sistem av-tomaticheskogo regulirovaniia [Identification of oscillatory automatic control systems]. Leningrad, Energiia, 1975, 96 pp. (In Russian)
6. Shteinberg Sh. E. Identifikatsiia v sistemakh upravleniia [Identification in control systems]. Moscow, Energoatomizdat, 1987, 80 pp. (In Russian)
7. Zoteev V. E. Parametricheskaia identifikatsiia dissipativnykh mekhanicheskikh sistem na osnove raznostnykh uravnenii [Parametric identification of dissipative mechanical systems based on difference equations]. Moscow, Mashinostroenie, 2009, 344 pp. (In Russian)
8. Bendat J., Piersol A. Engineering Application of Correlation and Spectral Analysis. New York, Wiley-Interscience, 1980, xiv+302 pp.
9. Yavlenskii K. N., Yavlenskii A. K. Vibrodiagnostika i prognozirovanie kachestva mekhanich-eskikh sistem [Vibration diagnostics and prediction of quality mechanical systems]. Leningrad, Mashinostroenie, 1983, 239 pp. (In Russian)
10. Vibratsii v tekhnike [Vibrations in Engineering], Handbook in 6 Vols. Moscow, Mashinostroenie (In Russian); vol. 1, 1978, 352 pp.; vol. 2, 1979, 351 pp.; vol. 5, 1981, 496 pp.
11. Jenkins G. M., Watts D. G. Spectral analysis and its applications, Holden-Day series in time series analysis. San Francisco, Holden-Day, 1968, xviii+525 pp.
12. Dobrynin S. A., Fel'dman M. S., Firsov G. I. Metody avtomatizirovannogo issledovaniia vibratsii mashin: Spravochnik [Methods of Automated Vibration Research of Machines: A Handbook]. Moscow, Mashinostroenie, 1987, 224 pp. (In Russian)
13. Marple S. Lawrence, Jr. Digital Spectral Analysis: With Applications, Prentice-Hall Series in Signal Processing. New York, Prentice-Hall, 1987, xx+492 pp.
14. Dech G. Rukovodstvo k prakticheskomu primeneniiu preobrazovaniia Laplasa i z-preobrazovaniia [Guide to the practical application of the Laplace transform and Z-transform]. Moscow, Nauka, 1971, 288 pp. (In Russian)
15. Nayfeh A. H. Introduction to perturbation techniques. New York, John Wiley & Sons, 1993, xiv+519 pp.
16. Pisarenko G. S., Matveev V. A., Yakovlev A. P. Metody opredeleniia kharakteristik kolebanii uprugikh sistem [Methods of determining the characteristics of vibration damping in elastic systems]. Kiev, Nauk. dumka, 1976, 88 pp. (In Russian)
17. Draper N. R., Smith H. Applied Regression Analysis, Wiley Series in Probability and Statistics. New York, John Wiley & Sons, 1998, xix+716. doi: 10.1002/9781118625590
18. Demidenko E. Z. Lineinaia i nelineinaia regressii [Linear and Nonlinear Regression]. Moscow, Finansy i statistika, 1981, 302 pp. (In Russian)
19. Marquardt D. W. An Algorithm for Least-Squares Estimation of Nonlinear Parameters, J. Soc. Indust. Appl. Math., 1963, vol.11, no. 2, pp. 431-441. doi: 10.1137/0111030.
20. Hartley H. O., Booker A. Nonlinear Least Squares Estimation, Ann. Math. Statist, 1965, vol.36, no. 2, pp. 638-650. doi: 10.1214/aoms/1177700171.
21. Zoteev V. E. Construction of difference equations for increasing the accuracy of parametric identification of vibrational systems with a weak general nonlinearity, Vestn. Samar. Gos. Tekhn. Univ., Ser. Fiz.-Mat. Nauki [J. Samara State Tech. Univ., Ser. Phys. Math. Sci.], 2000, no. 9, pp. 169-173 (In Russian). doi: 10.14498/vsgtu44.
22. Zoteev V. E. Development and research in digital linear models of dissipative systems wavering, Vestn. Samar. Gos. Tekhn. Univ., Ser. Fiz.-Mat. Nauki [J. Samara State Tech. Univ., Ser. Phys. Math. Sci.], 1999, no. 7, pp. 170-177 (In Russian). doi: 10.14498/vsgtu222.
23. Zoteev V. E. Parametrical identification of linear dynamical system on the basis of stochastic difference equations, Matem. Mod., 2008, vol.20, no. 9, pp. 120-128 (In Russian).
