CC BY
22
УДК 519.246
Б.Г. Кухаренко Институт машиноведения им. А.А. Благонравова РАН Московский физико-технический институт (государственный университет)
Исследование динамики нелинейных систем на основе переходных
временных рядов
Рассматриваются переходные временные ряды для переменных нелинейных динамических систем. На примере численных решений модели Даффинга продемонстрированы особенности спектрального оценивания переходных временных рядов по методу Прони и на основе преобразования Гильберта. Показано, как определяется зависимость частоты от амплитуды численного решения, представляющего свободное нелинейное колебание. Определение этой зависимости позволяет идентифицировать нелинейную динамическую систему.
Ключевые слова: нелинейные динамические системы, переходные временные ряды, спектральный анализ, метод Прони, преобразование Гильберта, модель Даффин-га, скелетные кривые, идентификация динамических систем.
I. Спектральное оценивание временных рядов для переменных нелинейных
динамических систем
Х.1. Модель локально переходных временных рядов с колебательным изменением и
метод Прони
Нелинейная система описывается моделью с конечным числом степеней свободы, и в ее пространстве состояний имеется несколько положений равновесия. Если в этом пространстве состояний нелинейная система совершает хаотическое движение, то оно представляется последовательностями локальных орбит возле нескольких ее положений равновесия [1-3]. Причем каждое из локальных движений возле одного из неустойчивых положений равновесия завершается в результате перескока на орбиту возле другого положения равновесия. В результате изменение переменных нелинейной системы во времени имеет гладкий колебательный характер только ограниченное число периодов, после которых с вероятностью единица происходит скачек. То есть если и существует некоторое глобальное колебательное изменение во времени конкретной нелинейной системы, то, как правило, оно представляется только его частичными реализациями. Для нелинейных временных рядов, представляющих хаотическую динамику, естественной является модель локально-переходных временных рядов [4].
Стандартные процедуры локального (во времени) спектрального анализа временного ряда ж[1,Ло] (определенного с шагом дискретизации времени ДЬ) используют его сегментирование посредством сдвига временного окна фиксированной длины N ^ N0 (^Дt определяет временной масштаб для локального анализа этого нестационарного временного ряда) [5]. Для локально переходного временного ряда с колебательным изменением х[1,^] имеется временной масштаб NДЬ, такой, что в пределах последовательных сегментов
x[nJ,(nJ + N — 1)],nJ = N ■ 3,3 = [0 : топий^0^)]
огибающая амплитуды колебаний меняется монотонно. Такие сегменты (а значит, и сам ряд локально во времени) характеризуются переходным колебательным изменением.
Если необходимо определить зависимость частоты от амплитуды нелинейных колебаний, то требования к методу спектрального анализа нелинейных временных рядов (представляющих хаотическую динамику) основываются на используемой модели локально переходного колебательного изменения этих временных рядов. То есть сегменты ряда — это не квазистационарные составляющие нестационарного временного ряда (которые для наблюдаемого ряда могли бы оказаться
очень малыми по длине), а переходные, то есть с неизменным знаком фактора демпфирования (логарифмического декремента) главных спектральных составляющих сегмента. Использование спектрального анализа, соответствующего модели временного ряда, делает его результаты информативными (то есть имеющими физическую интерпретацию). Поэтому используется метод Прони (Prony R., 1796) [6]. Исторически — это первый метод спектрального анализа временных рядов, но он не имеет непрерывного аналога. Его широкое использование откладывалось из-за отсутствия устойчивого с вычислительной точки зрения алгоритма этого метода. Декомпозиция Прони последовательного сегмента x[1,N] временного ряда x[1,No] имеет вид
p
x[i] = ^ r [1](z[1])i—1 + n[i],i = 1,N,
1=1
где p — число определяемых полюсов сегмента временного ряда; z[l] = exp(At(—¿[l] + j2nf[l])),l = 1,p — полюса, которые определяются для сегмента временного ряда, где ¿[1] и f[l] — соответственно фактор демпфирования (логарифмический декремент) и частота; r[l] = a[l] ■ exp(j^>[l]),l = 1,p — вычеты в этих полюсах, где a[l] и ip[l] — соответственно амплитуда и фаза; n[l] — аддитивный шум. Длина N последовательного сегмента может быть любой, но лучшие результаты спектрального оценивания получаются, когда размер скользящего временного окна N At примерно равен двум средним периодам временного ряда x[1,N]. Более плотная выборка значений спектральных параметров временного ряда получается за счет использования дробных сдвигов временного окна.
Создание устойчивого с вычислительной точки зрения алгоритма метода Прони является результатом работ ряда авторов [7-12]. Матричные алгоритмы метода Прони основаны на специальной матричной форме — пучке матриц (matrix pencil) [13]. Для ганкелевой матрицы данных Y размерности (N — L) х (L + 1)
Y[1,(N - L),I] = x[I,(I + N — L)]T,I = 1,(L + 1),
где L — размерность вложения, определяются редуцированные матрицы
Y(0) = Y[1,(N — L) ,1L],Y(1) = Y[1,(N — L) ,2,(L + 1)]
и используется свойство понижения ранга пучка матриц
(Y(1) — zY(0))
при z = z[l],l = 1,p. Таким образом, полюса z = z[l],l = 1,p определяются как обобщенные собственные числа матриц Y(0) и Y(1). Критерий, позволяющий установить число p главных полюсов сегмента временного ряда, использует разложение сингулярных чисел (SVD):
Y(v) = U(v) ■ S(v) ■ (V(v))T,v = [0,1],
где
S(v) = diag(s(v)[1],s(v)[2],...),s(v)[1] > s(v)[2] > ... > 0
и
(U(V))T U(V) = I(n-L)x(N-L),(V(V))TV(V) = I(l+1)x(L+1),
где I(n-l)x(n-l) и I(l+1)x(l+1) — тождественные матрицы размерности (N — L) х (N — L) и (L + 1) х (L + 1) соответственно, а верхний индекс «T» обозначает транспонирование [14]. Определяются усеченные матрицы
S(v) = diag(s(v)[1],...s(v)[p],...,0),v = [0,1].
Использование матриц S(v),v = [0,1] позволяет определить матрицы Y(v),v = [0,1] с фиксированным рангом p. Определение полюсов zi,l = 1,p как обобщенных собственных чисел матриц
= [0,1] гарантирует вычислительную стабильность определения этих полюсов в пределе временного ряда без шума, поскольку матрицы Y(v),v = [0,1] имеют фиксированный ранг р. При минимальном числе полюсов р = 2 аппроксимация по методу Прони последовательных сегментов (гладкого) временного ряда х[1,Жо] позволяет оценить временную зависимость частоты этого временного ряда. Работы по совершенствованию алгоритмов на основе метода Прони для временных рядов, представляющих переменные динамических систем с целью их идентификации, начались в 1990-х и продолжаются до сих пор [15-18].
При определении зависимостей частоты локальных переходных колебаний от их амплитуды для координатно-временных зависимостей моделей с хаотической динамикой методу Прони нет альтернативы. В [4] полученные по методу Прони частичные зависимости частоты временных рядов от их амплитуды (с перекрытием диапазонов изменения амплитуды и частоты) объединяются в полную функциональную зависимость. Результат проверяется в основном (с точки зрения его асимптотических свойств) при малых значениях амплитуды колебаний (малые колебания описываются линеаризованными динамическими моделями, поэтому частоты малых колебаний определяются точно). Это составляет основу для метода выделения в пространстве состояний динамической модели как области переходных последовательностей неустойчивых локальных орбит порядка 1, так и областей устойчивых и неустойчивых глобальных орбит [4]. Однако для более простых динамических моделей переходное колебательное изменение характеризует временной ряд в целом, то есть колебание является непрерывным глобально во времени. В этом случае альтернативой методу Прони при оценке динамических характеристик систем на основе временных рядов является преобразование Гильберта [19]. Это позволяет сравнить результаты двух методов спектрального оценивания для переходных временных рядов с колебательным изменением.
Х.2. Спектральное оценивание переходных временных рядов с колебательным изменением на основе преобразования Гильберта
Преобразование Гильберта функции х = х(*) имеет вид
х(т)—■
СХЭ
1
ж(£) = — п
t — Т
— ОС
и представляет собой фильтр, который сдвигает фазы всех частотных составляющих х(*) на —п/2. Преобразование Гильберта позволяет ввести медленно меняющиеся функции, называемые огибающей а = а(*) и фазой ^ = ^(¿), по формуле
х(*) + ¿Ж(г) = а(*)ехр(,?Х*)).
Таким образом,
а(*) = л/(х(*))2 + (Ж(*))2
и
^(¿) = аг^(Ж(*)/х(*)).
Формула для мгновенной угловой частоты имеет вид
х(*)(—Ж(*)/^) — (^х(*)/^*)Ж(*)
^(¿) = ^(¿)/^* =
(а(*))2
Преобразование Гильберта координатно-временной зависимости х = х(*) динамической системы позволяет определить модальные параметры этой динамической системы. Например, в случае квазилинейной динамической системы
формулы для временных зависимостей ее модальных параметров имеют вид
5(^ а(*) + 2ш(*)
то есть зависят от производных огибающей и мгновенной частоты. В линейной системе (то есть -ш/-* = 0) фактор демпфирования (логарифмический декремент) 5(*) пропорционален логарифму огибающей а(*).
В случае вынужденных колебаний квазилинейной динамической системы формулы для модальных параметров более громоздкие, поскольку зависят от внешнего возбуждения у = у(*). Формулы для модальных параметров различных динамических систем с демпфированием и описание численной реализации метода спектрального оценивания временных рядов на основе преобразования Гильберта приведены в [19]. В качестве примера ниже приводится сравнение результатов спектрального оценивания на основе преобразования Гильберта (по методу [19]) и по методу Прони для (глобально переходных) численных решений для свободных и вынужденных колебаний в модели Даффинга с демпфированием [20].
II. Спектральное оценивание численных решений модели Даффинга
П.1. Спектральное оценивание численного решения модели Даффинга без
внешнего воздействия
Модель Даффинга имеет вид
где ш0 = 2п/0 — угловая частота свободных колебаний (при у(*) = 0) с малой амплитудой (при условии х ^ 1 колебания описываются линейным уравнением с 7 = 0), 5 > 0 — демпфирование и 7 > 0 — параметр нелинейности этой модели. В [19] в качестве примера рассматривается численное решение модели Даффинга при 5 = 2,5 (сильное демпфирование), собственной частоте /о = 15 и 7 = 2000 (сильная нелинейность), полученное при начальных условиях х(0) = 10 и (-/-*)х(0) = 0 (рис. 1). Такое решение не может быть оценено методами теории возмущений [20]. На рис. 2 показана оценка огибающей амплитуды колебания на рис. 1 посредством преобразования Гильберта. На рис. 3 показана оценка временной зависимости частоты колебания / = /(*) на рис. 1 посредством преобразования Гильберта. Для * £ [0,0,1] эти оценки (рис. 2 и 3) искажены, а для * £ [0,1,0,4] оценка огибающей (рис. 2) оказывается слегка заниженной. На рис. 4 сплошной линией показана оценка огибающей амплитуды а = а(£) колебания на рис. 1 по методу Прони. На рис. 5 показана оценка временной зависимости частоты / = /(*) колебания на рис. 1 по методу Прони. При убывании амплитуды колебания она асимптотически приближается к значению / = 15 (пунктирная линия на рис. 1).
Для оценок временных зависимостей рис. 4 и 5 по методу Прони отсутствует «краевой» эффект для * £ [0,0,1], характерный для оценок таких зависимостей на рис. 2 и 3. Для * £ [0,0,4] оценка огибающей амплитуды а = а(£) (рис. 4) оказывается слегка заниженной. Действительно, в основе метода Прони лежит линейная модель колебания х(£) = аexp(—5 ■ *) еов(2п/ ■ * + ^>), где а — амплитуда, ^ — фаза, 5 — фактор демпфирования (логарифмический декремент), / — частота. Для * £ [0,0,4] профиль нелинейного колебания более крутой, чем у синусоиды, поэтому при аппроксимации этого нелинейного колебания единственной косинусоидой (рис. 4) происходит заниженная оценка его амплитуды. При * ~ 0,4 оценка амплитуды колебаний по методу Прони становится практически точной. Интересно, что такой же эффект имеет место при оценке амплитуды свободных колебаний (рис. 1) на основе преобразования Гильберта (рис. 2). В остальном
— х о о . .
-2 +25^ + шох + 7х = уй,
оценки огибающей амплитуды а = а(і) и временной зависимости частоты / = /(і) колебания на рис. 1, полученные посредством преобразования Гильберта и по методу Прони, оказываются подобными. На рис. 6 представлена зависимость частоты от амплитуды / = /(а) колебания на рис. 1, соответствующая временным зависимостям огибающей амплитуды а = а(і) (рис. 4) и частоты / = /(і) (рис. 5), полученным по методу Прони. Обратная зависимость а = а(/) дает скелетную кривую для модели Даффинга с кубической нелинейностью упругой силы [20]. Скелетная кривая свободных колебаний является характерной для определенного типа нелинейности упругой силы, и определение этой зависимости позволяет идентифицировать нелинейную динамическую систему [19].
Рис. 1. Координатно-временная зависимость Рис 2 Оценка огибающей амплитуды к°лебания х = х(*) для свободного нелинейного колебания на рис. 1 посредством преобразования Гильберта
Рис. 3. Оценка временной зависимости частоты колебания на рис. 1 посредством преобразования Рис. 4. °ценка огибающей амплитуды колебания
Гильберта
на рис. 1 по методу Прони
Рис. 5. Оценка временной зависимости частоты колебания на рис. 1 по методу Прони
Рис. 6. Оценка зависимости частоты от амплитуды колебания на рис. 1 по методу Прони
II.2. Решения модели Даффинга с внешним воздействием
В [19] в качестве примера вынужденного колебания рассматривается численное решение модели Даффинга (рис. 7) при 5 = 2,5 (сильное демпфирование), собственной частоте /о = 30 и 7 = 2000 (сильная нелинейность) при начальных условиях х(0) = 0 и (-/-*)х(0) = 0 и синусоидальном возбуждении у = у(*) с линейной зависимостью частоты от времени (рис. 8). При локальном во времени синусоидальном возбуждении у = у(*) (рис. 8) вынужденное колебание (рис. 7) также имеет локально синусоидальный характер с переменной частотой. Тем не менее на рис. 9 оценка огибающей амплитуды а = а(*) вынужденного колебания (рис. 7) посредством преобразования Гильберта во временном интервале наступления «резонанса» * £ [1,1,1,3] оказывается заниженной. На рис. 10 показаны оценки временной зависимости частоты / = /(*)
(пунктирная линия) и фазы = ^(¿) (точечная линия) для колебания на рис. 7 и частоты возбуждения на рис. 8 (сплошная линия) посредством преобразования Гильберта, которые отражают «резонансное» поведение решения модели Даффинга при линейном возрастании частоты возбуждения. Однако оценки на рис. 11 и на рис. 12 огибающей амплитуды а = а(і) и частоты / = /(¿) вынужденного колебания (рис. 7) по методу Прони являются более точными, чем оценки на рис. 9 и 10 соответственно (отсутствуют эффекты сглаживания, характерные для оценок аналитических методов). Таким образом, и для глобально переходных нелинейных временных рядов спектральное оценивание по методу Прони оказывается наиболее точным.
Рис. 7. Координатно-временная зависимость Рис. 8. Временная зависимость у = у(*) для х = х(*) для вынужденного колебания возбуждения
Рис. 9. Оценка огибающей амплитуды колебания на рис. 7 посредством преобразования Гильберта
Рис. 10. Оценка временной зависимости частоты колебания на рис. 7 (пунктирная линия) и частоты возбуждения на рис. 8 (сплошная линия) посредством преобразования Гильберта
Рис. 11. Оценка огибающей амплитуды колебания на рис. 7 по методу Прони
Рис. 12. Оценка временной зависимости частоты колебания на рис. 7 (пунктирная линия) и частоты возбуждения на рис. 8 (сплошная линия) по методу Прони
Литература
1. Lichtenberg A.J., Liberman A.P. Regular and stochastic motion / Applied Mathematical
Science. — V. 38. — New York, Heidelberg, Berlin: Springer-Verlag, 1983.
2. Universality of Chaos: 2nd edition / Ed. by Cvitanovic P. — Copenhagen, Denmark: Niels Bohr Institute, 1989.
3. Cvitanovic P., Garpard P., Schreiber T. Investigating of the Lorenz gas in term of periodic orbits // Chaos: An Interdisciplinary Journal of Nonlinear Science. — 1992. — V. 2, I. 1. — P. 85-90.
4. Кухаренко Б.Г. Исследование по методу Прони динамики систем на основе временных
рядов. Труды МФТИ. — 2009. — Т. 1, № 2. — С. 176-192.
5. Kay S.M., Marple S.L. Spectrum analysis — A modern perspective // Proceedings of the IEEE. —
1981. — V. 69, N 11. — P. 1380-1419.
6. Prony R. Essai experimental et analytique // Journal de l’Ecole Polytechnique. — Paris, 1796. — V. 1, I. 2. — P. 24-76.
7. Pisarenko V.F. On the estimation of spectra by means of nonlinear function of covariance matrix // Geophysics Journal of Royal Astronomical Society. — 1970. — V. 28. — P. 511-531.
8. Pisarenko V.F. The retrieval of harmonics from a covariance function // Geophysics Journal of Royal Astronomical. Society. — 1973. — V. 33. — P. 347-366.
9. Kumaresan R., Tufts D.W. Estimating the parameters of exponentially damped sinusoids and pole-zero modeling in noise // IEEE Transactions on Acoustics, Speech, and Signal Processing. —
1982. — V. ASSP-30, N 4. — P. 833-840.
10. Tufts D.W., Kumaresan R. Singular value decomposition and improved frequency estimation using linear prediction // IEEE Transactions on Acoustics, Speech, and Signal Processing. — 1982. — V. ASSP-30, N 4. — P. 671-675.
11. Hua Y., Sarkar T.K. Matrix pencil method for estimating parameters of exponentially damped / undamped sinusoids in noise // IEEE Transactions on Acoustics, Speech, and Signal Processing. — 1990. — V. ASSP-38, N 5. — P. 814-824.
12. Hua Y., Sarkar T.K. On SVD for estimating generalized eigenvalues of singular matrix pencil in noise // IEEE Transactions on Acoustics, Speech, and Signal Processing. — 1991. — V. ASSP-39, N 4. — P. 892-900.
13. Гантмахер Ф.Р. Теория матриц. — М.: Наука, 1966.
14. Scharf L.L. The SVD and reduced rank signal processing // Signal processing. — 1991. —
V. 25, N 2. — P. 113-134.
15. Barone P., Massaro E., Polichetti A. The segmented Prony method for the analysis of nonstationary time-series // Astronomy and Astrophysics. — 1989. — V. 209, N 1-2. — P. 435-444.
16. Sarkar T.K., Pereira O. Using the matrix pencil method to estimate the parameters of a sum of complex exponentials // IEEE Antenna and Propagation Magazine. — 1995. — V. 37, N 1. — P. 48-55.
17. Kukharenko B.G. Use of the Prony method for modal identification of slow-evolutionary linear
structures // Journal of Structural Control. — 2000. — V. 7, N 2. — P. 203-218.
18. Кухаренко Б.Г. Технология спектрального анализа на основе быстрого преобразования
Прони // Информационные технологии. — 2008. — № 4. — С. 38-42.
19. Feldman M. Hilbert Transform Applications in Mechanical Vibration. — Wiley, 2011.
20. Rand R.H. Lecture Notes on Nonlinear Vibrations. — Cornell University, Ithaca, NY: Dept. Theoretical & Applied Mechanics, 2005.
Поступила в редакцию 27.06.2011.
