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

Математическое моделирование

УДК 519.246

ПАРАМЕТРИЧЕСКАЯ ИДЕНТИФИКАЦИЯ ДИФФЕРЕНЦИАЛЬНЫХ ОПЕРАТОРОВ ДЛЯ СИСТЕМ С ТУРБУЛЕНТНЫМ ТРЕНИЕМ НА ОСНОВЕ РАЗНОСТНЫХ УРАВНЕНИЙ

В. Е. Зотеев, М. А. Заусаева, A.A. Егорова

Самарский государственный технический университет,

443100, Самара, ул. Молодогвардейская, 244.

E-mail: zoteev-ve@mail. ru

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

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

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

Владимир Евгеньевич Зотеев (д.ф.-м.н., доцент), доцент, каф. прикладной математики и информатики. Мария Анатольевна Заусаева, аспирант, каф. прикладной математики и информатики. Александра Арсеновна Егорова, аспирант, каф. прикладной математики и информатики.

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

Динамические процессы в системах с турбулентным (гидродинамическим трением) описываются нелинейным дифференциальным уравнением

где у(Ь)—реакция системы на некоторое входное возмущение /(¿) (например, типовое тестовое воздействие); т и с — масса и коэффициент жесткости системы; Ъу'{Щу1 {Щ —диссипативная сила (сила трения), обуславливающая рассеяние энергии колебаний, причём Ь2 <С тс [1].

Рассмотрим задачу параметрической идентификации нелинейного дифференциального оператора Ь по результатам измерений реакции системы с турбулентным трением на импульсное (ударное) воздействие, которое может быть описано специальной ¿-функцией:

Известно, что в этом случае весовая функция может быть найдена из решения однородного нелинейного дифференциального уравнения

Ьу(Ь) = ту” (*) + + су (г) = /(*),

при £ = О, при £ ф О,

ту" (Ь) + Ъу'Щу'Щ + су (г) = О с начальными условиями у(0) = 0, у'{0) = 1.

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

-1 _|_ СИ ,

37ТГПл/ш

где а и ф — произвольные постоянные.

С учётом начальных условий произвольные постоянные находятся из системы уравнений

: асовф = 0;

= —ахГсТтътф — а2 „со$ф = 1.

у ' Т бтгтпу'т Т

Отсюда получаем, что а = \/т/с и ф = —тг/2 + 2тгк, к Є 2.

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

, . (1/^) , ч

т = і + даг)(8Ш1^ <2)

где 5 = 86/(3л/тс), ш = \fcjm, Т = 27г/о; — декремент, частота и период колебаний соответственно. Параметры дифференциального оператора в относительных к массе системы единицах могут быть найдены по следующим формулам:

с 2 Ь 3 .

— = иг, — = -оиз.

т т 8

При измерении мгновенных значений Ук, к = 0,1,... , импульсной характеристики (2) первый отсчёт уо в выборке результатов наблюдений может не совпадать по времени с моментом і = 0, соответствующим приложению импульсного воздействия на систему. Обозначим момент времени измерения первого отсчета уо относительно момента приложения импульсного воздействия через ¿о- Тогда, переходя к новой системе временных координат: і = = + ¿0) где момент времени = 0 соответствует началу формирования

выборки результатов измерений, получаем следующее уравнение, описывающее импульсную характеристику системы с турбулентным трением:

ао

у (¡l) = і + (/¡„теЮ5И'+м'

где 5о = 5/{1 + (¿/Т)іо)—декремент колебаний, соответствующий моменту времени начала формирования выборки результатов измерений;

(І/о;) 2-7Г

ао = —

l + (5/T)to 2iru) + u)25to

и фо = coto —7г/2 — начальные амплитуда и фаза колебаний, соответствующие времени ¿0-

Таким образом, в системе временных координат і Є [0; оо], начало которой і = 0 соответствует времени первого отсчета в выборке результатов наблюдений, уравнение импульсной характеристики системы с турбулентным трением имеет вид

у® = -і, (з)

1 "Т" 2тг

где

П~ X 2уг , , ,.ч

ш = л/с,/тп, 5 = —, 6 = ,—, а0 = ----------■———, фо = иЛо - тг/2. 4

Зл/тс Зл/тс 2пш +

Очевидно, что параметры ао и в уравнении (3) функционально связаны между собой:

2тг

а° 2тгш + со5 (фо + 7г/2)

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

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

Ук =-----о!чт, с°8 (штк + ,Фо) + £к, к = 0,1,2,... , N — 1, (5)

1 + к

где т— период дискретизации функциональной зависимости (3); N — объём выборки результатов наблюдений; —случайная помеха в результатах наблюдений. Оценки динамических характеристик со и 5 находятся из условия минимизации функционала

N-1 N-1

■] 5) = ^ (Ук ~ Ук)2 ш1п’

к=0 к=О

где Ук — предсказанные по построенной модели значения отклика системы.

С этой целью целесообразно использовать численные методы, в основе которых лежат разностные уравнения [1]. Применяя алгоритмы построения линейно-параметрических дискретных моделей, описанные в [1], для временной последовательности результатов наблюдений (5) можно получить следующую систему стохастических разностных уравнений:

Уо — А 2 + £о]

< Уі = Аз + Єї]

Ук + Ук-2 = АоУй-і - Лі [кук - Ао (к - 1) ук-\ + (к -2) ук-2} + Щ]

^Лк = [1 + (к — 2) Лі] Єк~2 — Ао [1 + (к — 1) Аі] Єк-і + [1 + кХ\] Єк,

к = 2,3,...,М-1, (6)

коэффициенты которых связаны с параметрами функциональной зависимости (3) соотношениями: Ао = 2coswr, Ai = аоС026т/(2n), Х2 = аосовфо, A3 = = a0( 1 + Ai)-1 cos (wt + фо).

Задача вычисления среднеквадратичных оценок коэффициентов Ао, Ai, А2, A3 разностных уравнений (6) решается в формате прикладного линейного регрессионного анализа на основе обобщённой регрессионной модели, которая в матричной форме имеет вид

b = F А + г]; Г/ = РХ£.

(7)

Элементы векторов и матриц, входящих в систему (7), описываются в соответствии с уравнениями (6) [1]. Применение итерационных процедур среднеквадратичного оценивания коэффициентов А^ позволяет обеспечить минимизацию функционала .1(00, 5) и за счёт этого практически устранить смещение оценок А.,- по сравнению с классическим методом наименьших квадратов [1].

Оценки параметров импульсной характеристики (5) и коэффициентов дифференциального оператора в уравнении (1) могут быть вычислены по следующим формулам:

1 Ао ' г,

со = — arccos —, ао = 2 д

т 2 1

А2 + А2 + AfA2 + 2А1А2 — А0А2А3 — А0А1А2А3 4^2 ’

S =

2ттХ\ аош2т’

ао sin фо =

А0А2 — 2А3 — 2А1А3

4-Ag

■Фо =

arctg f 0,0 д” J , при Х2 ^ 0;

arctg f а° j +7г, при Х2 < 0, aosin^o ^ 0;

arctg (а° j — 7г, при Х2 < 0, ао sin^o < 0,

фо = Фо + 27Г • int

1 — ао со

OoLOÖ

ФІ

2тг

с = со2т, Ъ = ^-бсот, 8

где int [х\ — округление числа х до ближайшего целого. В соответствии с (4) время запаздывания можно оценить по формуле ¿о = {фо + 7г/2) /со.

Существенным недостатком модели (6) является то, что она не учитывает существующей функциональной зависимости между параметрами ао и фо в решении (5). Это находит отражение в использовании двух независимых коэффициентов модели А2 и A3. С учётом начальных условий эти коэффициенты вполне определяются параметрами со и 5 весовой функции (2):

Х2 = ао cos фо =

2тг sin coto СО (2-7Г + COÖto)

т — Aiio

сот

sinwto,

Лз = ТТТ~ COS (WT + ^ = , ,Т (л У\° N Sin (Í0 + Т)] •

I + Al LOT (1 + Лі)

С учётом этого недостатка при известной величине to построена модель, описывающая в форме стохастических разностных уравнений результаты наблюдений импульсной характеристики системы с турбулентным трением:

Уо-^Т) sinw«í0 = —¿fcf sin (ш^Но) Al + є0;

Уі - Ш sin [й(І) (r + ÍQ)] =- У1 + Wh sin [^(i) (r + ÍQ)] Лі + (! + Лі)

Ук + Ук-2 = ХоУк-1 - Аі [(к - 2) ук_2 - Ао (к- 1) ук-1 + кук\ + щ] т)к = [1 + (к — 2) Ai] £fc_2 — Ао [1 + (А: — 1) Аі] єк-і + [1 + к\\] єк,

к = 2,3, ...,N-1. (8)

Алгоритм вычисления среднеквадратичных оценок коэффициентов модели (8) включает дополнительную итерационную процедуру уточнения величины и)(г\ г = 0,1,... , первоначальная оценка которой находится на основе модели (6). Для обеспечения сходимости этой итерационной процедуры реализуется алгоритм вычислений

Ао} = (1 + с) Aq 1} - с(р(Хо 1}), і = 1,2,

где <р(\о 1 ) — оператор, непосредственно описывающий процедуру уточнения Ао на основе численных методов [1] среднеквадратичного оценивания коэффициентов модели (8):

В общем случае параметр с выбирается из следующих условий: с = — 1, если -1 < <р' (Ао) < 1; ^(Л2о)_1 < с < 0, если у/ (Ао) < -1, и < с < i+^(a0)>

если (р1 (Ао) > 1. На практике параметр с подбирается рекуррентно:

Cj = Cj~i/2, j = 1,2,...; со = -1,

до тех пор, пока итерационная процедура уточнения величины u)W не будет сходиться. Оценки параметров импульсной характеристики (2), а также коэффициентов дифференциального оператора, описывающего динамику систем с турбулентным трением, можно найти по следующим формулам:

1 Ао i 27Г А1 Л «2 i

lü = -arccos—, д =--------------------------------------------, с = шт., b = -oujm.

т 2 Cj(j — Aiio) 8

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

можно, если ввести в описание модели ещё один неизвестный коэффициент Аг, зависящий от величины ¿о- Это позволит не только найти параметры импульсной характеристики с учётом их зависимости от начальных условий, но также оценить время запаздывания момента начала формирования выборки результатов эксперимента. В этом случае система стохастических разностных уравнений принимает вид:

Алгоритм среднеквадратичного оценивания коэффициентов А^ модели (9) базируется на методах прикладного регрессионного анализа в формате обобщённой регрессионной модели (7), включает три итерационные процедуры и аналогичен алгоритму вычислений, описанному для модели (8). Оценки декремента колебаний 5 и времени запаздывания ¿о находятся из решения системы нелинейных уравнений

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

Компьютерное моделирование проводилось при следующих значениях коэффициентов дифференциального оператора: т = 1, Ь = 0,5, с = 10, что соответствует параметрам весовой функции (2): ш = 3,162 и 5 = 0,422. Формировалась выборка мгновенных значений ук, к = 0,1,2,... , Л/- — 1, импульсной характеристики с шагом дискретизации т = 0,2 и объёмом N = 21. К смоделированным в соответствии с формулой (5) дискретным значениям у к добавлялась аддитивная случайная помеха єк, величина которой

изменялась от 0 до 10%. На основе каждой из описанных выше моделей в форме стохастических разностных уравнений (6), (8) и (9) вычислялись среднеквадратичные оценки параметров динамического процесса со и 5 и дифференциального оператора — Ь и с. В качестве оценки погрешности результатов

'уо — А 2 + єо',

{ йіп ІР(І) (т + іо)] = “ г/і + йіп [и(і) (т + ¿о)] Аі + (1 + Аі) єі;

Ук + У к-2 = АоУй-і - Аі [(к - 2) ук_ 2 - А0 (к - 1) ук_г + кук\ + г/к] іЛк = [1 + (к — 2) Аі] єк-2 — Ао [1 + (к — 1) Аі] єк-\ + [1 + кХ\] єк,

к = 2,3,...,М-1, (9)

где коэффициенты модели описываются такими формулами:

со5т 27г віп 27гАі віп

Ао = 2 сов и;т, Аі

Ао 2 сов СОТ, Аі о і /г) і г# \ 2Х+

2тг + шдіо оо (2-7Г + содіо) согдіо

где со = т 1 агссов Ао/2.

^N-1 2

Н=° * • юо%

вычислении использовался второй центральный момент относительно точного значения параметра:

2

М

(ъ - ъ)2 = б [¿] + (м [6] - ь) .

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

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

Точки 1 соответствуют результатам вычислений с использованием разностных уравнений (6), а точки 2 и 3 — линейно-параметрическим дискретным моделям (8) и (9), учитывающим функциональную связь между параметрами динамического процесса, обусловленную начальными условиями в дифференциальном уравнении. Видно, что применение численных методов среднеквадратичного оценивания коэффициентов разностных уравнений, учитывающих в своей структуре начальные условия для дифференциального уравнения (1), позволяют существенно, в ряде случаев почти в два раза, повысить точность вычисления параметров дифференциального оператора для систем с турбулентным трением.

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

А6, АЬ, % Ас, %

Рис. 1. Зависимости погрешности вычисле- Рис. 2. Зависимости погрешности вычисления декремента колебаний <5 и параметра Ь ния параметра с дифференциального опера,-дифференциального оператора от величины тора от величины случайной помехи в ре-случайной помехи в результатах наблюдений зультатах наблюдений

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

Работа выполнена в рамках Аналитической ведомственной целевой программой «Развитие научного потенциала высшей школы» (проект № РНП.2.1.1/745) и при поддержке РФФИ (проект № 10-01-00644-а).

БИБЛИОГРАФИЧЕСКИЙ СПИСОК

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

2. Зотеев В. Е. Параметрическая идентификация линейной динамической системы на основе стохастических разностных уравнений // Матем. моделирование, 2008. — Т. 20, №9. - С. 120-128.

3. Зотеев В. Е. Итерационный метод среднеквадратичного оценивания коэффициентов стохастического разностного уравнения колебаний систем с турбулентным трением // Вестн. Сам. гос. техн. ун-та. Сер. Физ.-мат. науки, 2005. — №38. — С. 100-109.

Поступила в редакцию 01/11/2010; в окончательном варианте — 03/Х/2010.

MSC: 65C20, 65P40, 34C15, 37M05

PARAMETRIC IDENTIFICATION OF THE DIFFERENTIAL OPERATORS FOR THE SYSTEMS WITH TURBULENT FRICTION ON THE BASE OF FINITE-DIFFERENCE EQUATIONS

V.E. Zoteev, M. A. Zausaeva, A. A. Egorova

Samara State Technical University,

244, Molodogvardeyskaya st., Samara, 443100, Russia.

E-mail: zoteev-ve®mail.ru

Numerical technique for determination of the differential operator paramaters for the systems with turbulent friction is considered. Mean-square estimation of coefficients of stochastic difference equations describing observation results of the impulse characteristic of the system form, the basis of this method. The results of numerical analytical treatments of recommended, symbolic models and, algorithms for calculation on their basis are presented. These results confirm high efficiency of the developed, method for the parametric identification on the base of finite-difference equations.

Key words: systems with, turbulent friction, parametric identification, finite-difference equations, mean-square estimation.

Original article submitted 01/11/2010; revision submitted 03/X/2010.

Vladimir E. Zoteev (Dr. Sci. (Phys. & Math.)), Associate Professor, Dept, of Applied Mathematics & Computer Science. Maria A. Zausaeva, Postgraduate Student, Dept, of Applied Mathematics & Computer Science. Alexandra A. Egorova, Postgraduate Student, Dept, of Applied Mathematics & Computer Science.