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

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

УДК 519.246

В. Е. Зотеев

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

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

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

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

Рассмотрим построение линейно параметрической дискретной модели, связывающей в виде рекуррентной формулы несколько последовательных мгновенных значений кривой ползучести. В соответствии с разработанной теорией неупругого реологического деформирования материалов уравнение кривой ползучести y (t) при постоянном напряжении может быть представлено в виде суммы нескольких (обычно от двух до четырёх) экспоненциальных составляющих [4]:

При равномерной с периодом т дискретизации непрерывной функции (1) получаем её дискретный

p

{l)

i=1

аналог Уk = 52 ai [І — exp (—aiTk)] или

i=1

p

Уk

as — a

i=1

Yj ai^k’

{2)

где as = ai, k Є {0} U N, Hi = exp (—aiT).

i=1

Применяя к {2) Z-преобразование, имеем

Z{Vk} =

p ai

1 - г-1 4-1 г 1 - ^-Г

г=1

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

Ap (z ^ Z {yk} —

asAp z-1

sp

T^z

1

= Bp-l (z

1

{З)

p-l

где Ар (г -1) = П (1 — -1) = 1 — 52 Агг г и Вр_і (г *) = 52 Ьгг г — многочлены степени р и р — 1

І=1 І=1 І=0

относительно переменной г-1. Можно показать, что коэффициенты многочлена Ар (г-1) описываются формулами:

Л1 = 52 ^i> Л2 = — X] № і

i=1

Л3 — 52 HiHj Hk j

i, j, k i<j<k

Лp-l =(—1)p S ^іі №2 • • • №„-11 Лp =(—1)P+1 П №•

i, j i<j

І1 ,І2, ...,ip-1 І1 <І2 <...<І„-1

i=1

{4)

Из (3) с помощью обратного ^-преобразования можно получить

— А1 ук-1 — А2ук;-2 — • • • — АрЩ:-р — Ар+1 = + Ь1^й-1 + Ь2^й-2 + • • • + Ьр-А-(р-1), (5)

где йд —символ Кронекера (к € {0} и М), а коэффициенты Ар+1 определяются формулой

Лp+1 = as 1 1 — ^ Лі j=1

С учётом (2) введём обозначения:

{б)

Ap+2 = У1 = as — X] aiHi = 52 ai (1 — Hi^

i=1 i=1

pp 2

Лp+з = У2 = as — Ё aiH2 = 52 ai (1 — Hi),

i=1 i=1

pp j-1

ЛP+j = 1 -1 = as S ai Hi = S ai ( 1 Hi

i=1

i=1

j-1

Л2p =

pp Уp-l = as — 52 aiHp = £ ai (1 — Hp ) •

i=1 i=1 v 7

І=1

Используя (5) и (7), коэффициенты многочлена Вр_1 (г-1) можно представить в виде

Ь0 = 0, Ь1 = Ар+2) Ь2 = Ар+3 — А1 Ар+2) Ь3 = Ар+4 — А1 Ар+3 — А2 Ар+2) •

3 _1

= ар+з+1 — а1ар+з — А2Ар+з_ 1— • • •— Аз_ 1 Ар+2 = Ар+з+1 — 52 АгАр+з+1_^ • • •,

i=1

p-2

Ьp-1 = ^p — Л1Л2p—1 — Л2Л2p—2 — • • • — Лp-2Лp+2 = ^p — 52 ЛІЛ2p—i

i=1

{7)

1

s

При этом линейно параметрическая дискретная модель, связывающая в виде рекуррентной формулы р +1 мгновенное значение дискретной функции (2), принимает вид

У к =

(8)

Ар+1+к, если к = 1, 2, 3, . . . , р — 1,

А1 Ш:-1 + А2Щ:-2 + • • • + Арук;-р + Ар+Ъ если к = p, р + 1 . . . , N — 1,

где коэффициенты А? =1, 2, 3, .. •, 2р) определяются соотношениями (4), (6) и (7).

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

(г = 1, 2, • • •, р) алгебраического уравнения степени р:

Ир — А1^р-1 — А2^р-2 — ... — Ар-1^ — Ар = 0,

а затем воспользоваться формулами

(іі — — -

1п Ні

(9)

(10)

При параметрической идентификации кривой ползучести в процессе эксперимента формируется выборка результатов измерений уд = уд + ед, которые содержат случайную помеху ед, к = 0, 1, ..., N — 1, где N — объём выборки. В этом случае линейно параметрическая дискретная модель (8) принимает вид стохастического разностного уравнения:

єо,

Ар+1+к + єк,

если к = 0,

р — 1,

_ -•р+1+^ если к = 1, 2,

р р

£ А;Ш:-;' + Ар+1 — £ А?ек-.? + ек, если к = p, р +1, •••,N — 1 ;=1 ?=1

и может быть представлена в форме обобщённой регрессионной модели:

Ь = ^А + п,

П = Ра£.

(11)

(12)

Здесь А = (Аі, А2, ..., А2р)Т — вектор неизвестных коэффициентов линейно параметрической дискретной модели; є = (єо, єі, •••, єм-і)Т — N-мерный вектор случайной помехи в результатах наблюдений; п = (пі, П2, ..., Пм)Т — N-мерный вектор эквивалентного случайного возмущения в стохастическом разностном уравнении; Ь = (уо, Уі,

ум-і)Т — N-мерный вектор правой части; ^

= /ь/2: . :/р:/р+ь . :/2р

формулами:

— матрица регрессоров размера N х 2р, столбцы которой описываются

/і = (0, 0, . . . , 0, ур-і, Ур, . . ., Ум-2 )Т

/2 = ( 0, 0, . . . , 0, Уp-2, ур-і, . ■ ., Ум-3 )Т

/р = ( 0, 0, . . . , 0, Уо, У1, .. У 1 р і )Т

/р+і : = ( 0, 0, . . . , 0, 1, 1, . . . , 1 )Т

/р+2 : = ( 0, 1, . . . , 0, 0, 0, . . . , 0 )Т

2р = ( 0, 0, . . . , 1, 0, 0, . . . , 0 )Т

Матрица Р размера N х N в стохастическом разностном уравнении эквивалентного возмущения нижняя треугольная. Первые р строк матрицы описываются формулой

Ріі

если ] = і, если j = і.

і = 1, 2,

Р, І = 1, 2,

N.

Остальные строки матрицы Р имеют вид:

рр+і

рр+2

— ( Аp, Ар-Ъ Ар-

= (0, Рм = (0,

Ар

0,

р-2, р- і ,

..., —А2, -Аі , 1, 0, 0, ..., 0),

Ар-2, ..., -А2, — АЪ 1, 0, ..., 0),

0, -Ар, —Аp-1, — Аp-2, ... , -А2, — АЪ 1).

Помехоустойчивое определение параметров а и а кривой ползучести (1) сводится к среднеквадратичному оцениванию коэффициентов А? линейно параметрической дискретной модели (11), при котором могут быть использованы различные схемы вычислений. Наиболее простым является алгоритм, в основе которого лежит минимизации функционала ||п||2 = ||Ь — РА||2 —► тіп. В этом случае оценки коэффициентов стохастического разностного уравнения вычисляются по формуле

Л = (РТР) і РТЬ.

(13)

Однако такой подход мало эффективен из-за большого смещения оценок, обусловленного корреляцией между элементами матрицы Р и эквивалентным случайным возмущением Пк. Можно показать, что величина смещения в оценке вектора коэффициентов разностного уравнения приближённо равна — ^ — р) (РТР) А*, где А* = ( А1, А2, • • •, Ар, 0, ..., 0)т, —дисперсия случайного

возмущения в результатах наблюдений. Другой, более эффективный и помехоустойчивый, метод среднеквадратичного оценивания коэффициентов стохастического разностного уравнения использует минимизацию функционала ||е||2 = ||у — у||2 —► ш1п. В основе такого подхода лежит применение итерационной процедуры уточнения среднеквадратичных оценок коэффициентов А?, что позволяет практически устранить смещение в оценках и тем самым добиться высокой точности вычисления параметров исследуемой функции [5].

Из формулы (12) получаем Р-1Ь = Р;-1^1 А + е. Элементы р-1 обратной матрицы Р-1 зависят от коэффициентов А? и описываются формулами

р

-і \ 1, если і = і,

? | 0, если і = і,

при і = 1, 2, ..., р, І = 1, 2, ..., N

р

-і

і?

к=і

при і = 1, 2, ..., р, і = 1, 2,

N.

Е Акр,_к ?, если І ^

к=і

В основе итерационной процедуры лежит минимизация функционала

ІєІІ2^

Р-! Ь — Р-і Р А

л(к) л(к)

тіп,

где Рд(к) —обратная матрица, при вычислении элементов которой используются оценки А?к) коэффициентов разностного уравнения, найденные на к-ой итерации. Очевидно, что данный функционал представляет собой квадратичную форму относительно искомых коэффициентов А?. Следовательно, он достигает своего неотрицательного минимума. При этом нетрудно показать, что минимум данного функционала достигается в точке

Л(к+і) = (РТ^Л(к)Р)-і РТ^Л(к)Ь

(14)

где ^л(к) = (Р-1))Т (Р-1)).

На первом шаге итерационной процедуры среднеквадратичного оценивания коэффициентов разностного уравнения по формуле (13) находится начальное приближение:

д(0) = (>т Р) РТЬ.

Затем, полагая в формуле (14) к = 0 и Л(к) = Л(0), вычисляется следующее приближение:

А(1) = (РТОл(0) Р )-1 РТОл(0) Ь.

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

дМ - Л(к) ^ 0,01 Л(к)

и

2

Результаты численно-аналитических исследований показали не только хорошую сходимость итерационной процедуры (три-четыре итерации), но и её высокую эффективность по сравнению с обычным МНК-оцениванием. За счёт устранения смещения в оценках А^ погрешность вычисления характеристик кривой ползучести уменьшается в десятки и сотни раз в широком диапазоне изменения параметров численно-аналитического эксперимента.

Полученные среднеквадратические оценки коэффициентов линейно параметрической дискретной модели (11) используются при вычислении помехоустойчивых оценок параметров кривой ползучести а и аг (г = 1, 2, ..., р). Для этого сначала находятся корни Кг (г = 1, 2, ..., р) алгебраического уравнения (9):

КР + А1КР 1 + А2КР 2 + • • • + Ар—1К + Ар = 0.

Затем по формулам (10) вычисляются оценки параметров сц = т—11пКг (г = 1, 2, ..., р). На заключительном этапе по найденным Аг и Кг посредством решения системы линейных алгебраических уравнений

а1

+

А2

+ ... + Ар — Ар+1 р .

1— Л

г=1

+ ... + а - т-Н Аp+2,

+ ... + (1 - ар = Аp+3,

+ ... + 1 1 К рр 1 раА £ А

2, . . . , р) в функции (1).

(1 - ) А1 + (1 - ) а2 +

вычисляются оценки коэффициентов аг (г = 1,

Проведена апробация разработанного метода идентификации кривой ползучести на основе стохастических разностных уравнений в научно-техническом эксперименте. Исследовался образец по-лихлорвинилового пластиката длиной 1000 мм, площадью сечения 1,2 мм2 при температуре 20 °С и нагрузке Р = 5,6 Н. Результаты эксперимента по деформации ползучести уэксп (£) за 8 часов представлены в таблице, а также изображены на рисунке.

Штриховая линия на рисунке соответствует кривой ползучеЭкспериментальные значения (точки) и кривые ползучести, постро-

сти, полученной известным мето-енные методом выделения экспоненциальных составляющих (штри- ’ •’

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

(сплошная линия) составляющих [1], а сплошная —

методом, в основе которого лежит

2

Экспериментальные данные для деформации ползучести уэксп и результаты вычислений методом выделения экспоненциальных составляющих ух и на основе стохастических разностных уравнений у2

t, час 0,000 0,00417 0,01667 0,03333 0,08333 0,25 0,50 1,0

Уэксп 0 0,001 0,003 0,0045 0,007 0,0095 0,0125 0,0135

VI 0 0,001 0,00296 0,0045 0,00667 0,0095 0,01158 0,01337

У2 0 0,00065 0,00231 0,00399 0,00685 0,00985 0,01184 0,01400

1, час 2,0 3,0 4,0 5,0 6,0 7,0 8,0

Уэксп 0,016 0,017 0,018 0,0185 0,019 0,0195 0,02

VI 0,01517 0,01648 0,01751 0,01831 0,01894 0,01944 0,01982

У2 0,01595 0,01709 0,01793 0,01858 0,01908 0,01946 0,01976

среднеквадратичное оценивание коэффициентов стохастического разностного уравнения, описывающего результаты эксперимента [9]. Очевидно, что в среднеквадратичном смысле сплошная кривая более близка к экспериментальным точкам, чем штриховая в интервале от 1 до 6 часов.

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

У1 (*) = 8,47 ■ 10—3 -1 - е—°’24*) + 1,39 ■ 10—3 -1 - е—°’25*) + 6,54 ■ 10—3 <1 - е—3’92*) +

+ 4,34 ■ 10—3 -1 - е—37>96*< + 0,53 ■ 10—3 -1 - е—157>82*< ,

из которых первые две практически одинаковы, а доля последней составляет всего 3,6 % от общей мощности «сигнала», что статистически незначимо. Построенная на основе стохастических разностных уравнений модель кривой ползучести аппроксимирует экспериментальные данные с погрешностью 2,5 % и содержит только три экспоненциальные составляющие:

У2 (*) = 7,99 ■ 10—3 -1 - е—°’26*) + 6,06 ■ 10—3 -1 - е—2’33*) + 6,70 ■ 10—3 -1 - е—21’82*)

При этом величины а (г = 1, 2, 3) различаются примерно на порядок, что полностью соответствует теории построения уточнённых моделей неупругого реологического деформирования материалов и элементов конструкций.

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

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

1. Работное, Ю. Н. Ползучесть элементов конструкций [Текст] / Ю. Н. Работнов. —М.: Наука, 1966. —752 с.

2. Радченко, В. П. Экспериментальное исследование и анализ полей неупругих микро- и макро неоднородностей сплава АД-1 [Текст] / В. П. Радченко, С. А. Дудкин, М. И. Тимофеев // Вестн. Сам. гос. техн. ун-та. Сер.: Физ.-мат. науки. — 2002. — № 16. — С. 111-117. — ISBN 5-7964-0355-9

3. Радченко, В. П. Энергетический подход к прогнозированию ползучести и длительной прочности материалов в стохастической постановке [Текст] / В. П. Радченко // Проблемы прочности. — 1992. — № 2. — С. 34-40.

4. Самарин, Ю. П. Построение экспоненциальных аппроксимаций для кривых ползучести методом последовательного выделения экспоненциальных слагаемых [Текст] / Ю. П. Самарин // Проблемы прочности. — 1974. — № 9. — С. 2427.

5. Радченко, В. П. Определение динамических характеристик механической системы на основе стохастических разностных уравнений колебаний [Текст] / В. П. Радченко, В. Е. Зотеев // Известия вузов. Машиностроение. — 2007. — № 1. — С. 3-10.

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

[email protected]

V. E. Zoteev

PARAMETRICAL IDENTIFICATION OF CREEP’S CURVE ON THE BASIS OF STOCHASTIC DIFFERENCE EQUATIONS

The numerical method of parametrical identification of the creep’s curve is considered, allowing to increase accuracy of forecasting of processes not elastic deformations in tasks of an estimation of individual behavior of a concrete element of a structure. The method is based on linear parametrical discrete model, describing in the form of stochastic difference equations results of supervision of creep’s curve during experiment lays.

Samara State Technical University, Russia Received 11.10.2007

[email protected]