УДК 681.5:681.3:519.6
ТОЧНОСТЬ АЛГОРИТМОВ ПАРАМЕТРИЧЕСКОЙ ИДЕНТИФИКАЦИИ
Ю.В. Журавлев
МГТУ им. Н.Э. Баумана, Москва, Российская Федерация e-mail: [email protected]
Развита методика оценки погрешности класса алгоритмов параметрической идентификации, допускающих структурную декомпозицию на блок-формирователь матричной системы и блок-решатель этой системы. Исследованы теоретико-вероятностные характеристики блока-формирователя в алгоритме производящих функций при экспериментальной идентификации коэффициентов дифференциального уравнения моментов колебательной модели продольного короткопериодического движения летательного аппарата. Принята гипотеза о наиболее неблагоприятном распределении знаков ошибок регистрации экспериментальных данных. Аналитически вычислены два первых вероятностных момента матричных ошибок, которые предопределяют наследственную погрешность блока-решателя. Рекомендована регуляризация процесса идентификации параллельным включением блока фильтрации. Обсуждена отделимость этапа фильтрации экспериментальных измерений входных-выходных процессов от этапа собственно идентификации линейной стационарной динамической системы.
Ключевые слова: идентификация, функции Эрмита, число обусловленности, аналитическая оценка погрешности, фильтрация измерений, летательный аппарат.
ACCURACY OF ALGORITHMS OF PARAMETRIC IDENTIFICATION
Yu.V. Zhuravlev
Bauman Moscow State Technical University, Moscow, Russian Federation e-mail: [email protected]
A technique is developed for error estimation of the class of algorithms for parametric identification, allowing the structure decomposition to a block-shaper of the matrix system and a block-solver of this system. Theoretical-probabilistic characteristics of the block-shaper in the algorithm ofgenerating functions are investigated during the experimental identification of coefficients of the differential equation of moments of the oscillatory model of the longitudinal short-period movement of the aircraft. A hypothesis of the most unfavorable distribution of signs of errors of the experimental data registration is accepted. Two first probabilistic moments of matrix errors, which predetermine the inherited error of the block-solver, are calculated analytically. The regularization of the identification process is recommended by switching of the filtering unit in parallel. The feasibility of separation of the stage of filtering of experimental measurements of input-output processes from the stage of the identification itself of linear stationary dynamical system is discussed.
Keywords: identification, Hermite functions, condition number, analytical estimation of error, measurement filtering, aircraft.
Качество синтеза систем управления зависит от уровня информационной достоверности моделей, что чрезвычайно актуализирует
проблему идентификации объектов управления. Разработка вычислительных алгоритмов идентификации по необходимости сопровождается исследованиями их точностных характеристик [1, 2]. При этом могут привлекаться такие подходы, как аналитический, алгоритмический, имитационное моделирование [3]. Аналитическим путем изучаются априорные свойства решений. Алгоритмический и имитационный подходы опираются на статистический анализ.
В настоящей работе развита методология получения гарантирующих оценок погрешностей для обширного класса алгоритмов параметрической идентификации, структурно декомпозируемых в двухуровневую систему: первый уровень — блок формирования матричной системы (формирователь), второй уровень — блок решения матричной системы (решатель). С помощью аналитического аппарата нормированных пространств [4] оценивается норма погрешности решателя, обусловленной присутствием вариаций выходных результатов формирователя. Точностные характеристики формирователя изучены на примере алгоритма идентификации коэффициентов дифференциального уравнения моментов в модели продольного короткопериодиче-ского движения летательного аппарата [5]. В основе алгоритма лежит интегральный метод производящих функций Эрмита [6-9]. Родственные проблемы возникают и изучаются во всех без исключения задачах прикладной математики, использующих аппарат линейной алгебры, в частности в задачах распознавания образов [10]. Содержательное изложение тематически связанных вопросов можно найти в работе [11]. Как и большинство обратных задач, задача идентификации относится к классу некорректных [12]. Специфика таких задач вызывает необходимость совершенствования методологических приемов исследования традиционных и вновь разрабатываемых алгоритмов [13, 14].
Постановка задачи. Допустим, что структура идентифицируемой математической модели функционально линейна по параметрам и на этой основе некоторый алгоритм первого уровня, в дальнейшем именуемый формирователем, формирует матричную систему линейных алгебраических уравнений относительно вектора параметров X,
АХ = в, А е Япхп, в, X е Япх\
а на втором уровне другой блок алгоритма, в дальнейшем именуемый решателем, отыскивает X.
Пусть формирователь передает решателю точные матрицы А = = А0 е Дпхп и В = В0 е Дпх1, и пусть решатель при условии невырожденности матрицы А0 идеально точно находит Х0 = А-1В0 — истинное значение вектора идентифицируемых параметров. Но когда от формирователя поступают возмущенные матрицы А = А0 + ДА,
B = Во + ДВ, то идеальный решатель по уравнению (Ao + ДА) х х Х = В0 + ДВ выдаст возмущенное решение X = Х0 + ДХ. Здесь ДА £ Rnxn и ДВ, ДХ £ Rn — абсолютные погрешности. Абсолютная погрешность идентификации ДХ в определенном смысле оказывается наследственной погрешностью решателя. Источниками нетривиальности ДА, ДВ могут являются флуктуационные составляющие в записях экспериментальных данных, а также неточности аппроксимации математических операторов в алгоритме формирователя, в том числе вычислительные погрешности округления.
Требуется исследовать точностные характеристики алгоритма производящих функций [6] применительно к параметрической идентификации математической модели короткопериодического продольного движения летательного аппарата [5].
Оценивание наследственной погрешности решателя. Вычитанием двух матричных уравнений (А0 + ДА)(Х0 + ДХ) = В0 + ДВ и А0Х0 = В0 с невырожденной матрицей А0 получаем матричное уравнение относительно ДХ:
ДХ = А-1 • (ДВ - ДА • Х0 - ДА • ДХ). (1)
Опираясь на условие согласованности матричной и векторной норм
1|А • Х||<||А|Н|ХII, (2)
применим к (1) неравенство треугольника и получим
IIДХ|| < ||А0-1|| • (||ДВ|| + ||ДА|| • ||Х0|| + IIДА|| • ||ДХ||), откуда
ЦДХII • (1 - IIА-1 II • ЦДАЦ) < ЦА-1 II • (ЦДВЦ + ЦДАЦ • ЦХ0Ц). (3) Введем в рассмотрение вариации норм
5Х = ЦДХЦ/ЦХ0II, 5В = IIДВII/ЦВ0II, 5А = IIДАЦ/IIА0Ц,
а также обозначим v — число обусловленности, d — дискриминант корректности, K — коэффициент усиления ошибок:
v = IIА0II IIа-1 II, d = 1 - v¿A, K = v/d.
Необходима положительность дискриминанта корректности
d = 1 - v¿A > 0, (4)
чтобы из (3) получить точную верхнюю оценку абсолютной погрешности решения
IIДХII < K(ЦХ0^А+ IIА0II-1 • ЦДВЦ),
а также точную верхнюю оценку вариации решения
5Х < K[¿А + ЦДВЦ^ (IIA0II • ЦХ0Ц)-1].
По условию (2) имеем ||A0|| • ||Xo|| > ||B0||, откуда
(||Ao|| • ||Xo||)-1 <||Bo||-1
и потому
||AB|| • (||Ao|| • ||Xo||)-1 < ||AB|| • ||Bo||-1 = SB.
Окончательно запишем точную оценку вариации погрешности алгоритма в виде усиленной суммы матричных вариаций, внесенных формирователем:
SX < K(SB + SA). (5)
В частности, если SA = 0, то K = v, и тогда
SX < vSB. (6)
Видно, что 1 % погрешности в формируемых матрицах может перейти в K % погрешности результата, естественно K > 1. Полученная аналитическая оценка справедлива лишь при соблюдении условия d > 0, из которого следует, что при v = 1 и вариации SX = 10 % величина SA должна быть не более 9%, но при v = 100 и вариации SX = 10 % величина SA должна быть не более 0,09%. Этот пример указывает на чрезвычайную роль числа обусловленности матрицы.
Относительная методическая погрешность решателя является усиленной суммой относительных погрешностей исходных матриц. Поскольку главная часть нормы вектора определяется доминирующими компонентами вектора, то оценка SX будет объективно близкой по отношению к доминирующим компонентам. Например, если Xo = (1, 50, 20)т и абсолютная погрешность AX = (0,2, 5, 2)т, то в сферической норме ||Xo|| ~ 54, ||AX|| ~ 5,4, SX « 0,1, при этом относительные погрешности Ax1/x1 ~ 0,2, Ax2/x2 ~ 0,1, Ax3/x3 « 0,1, что отчетливо демонстрирует явление неустойчивости оценки погрешности меньшего параметра x1 = 1 на фоне оценок больших параметров. Проблема идентификации малых параметров, вообще говоря, относится к классу некорректных задач [13]. Принципиально необходимо проводить исследование чувствительности оценок малых параметров к влиянию разброса оценок доминирующих параметров.
В изложенном способе оценивания нормы погрешности решателя по умолчанию принимается гипотеза о максимально неблагоприятном распределении знаков всевозможных ошибок, так что подобные оценки являются гарантирующими, в случае супремальности — минимаксными, что актуально в стохастических условиях с высоким уровнем неопределенности.
Замечание. При любой матричной норме ||A|| > max |Aj|, где {Л} —
1<г<п
спектр собственных чисел матрицы A, тогда IIA-1 II > 1/ min |Л«|. Отсю-
1<г<п
да v > 1. Норма матрицы ||A||, подчиненная векторной сферической норме
= \ х2 и одновременно согласованная с ней, выражается через спек-
тральный радиус р (АтА) = ^тах матрицы АтА, а именно ||Л|| = ^р (АтА),
Мтт
где ^тах, Мтт — максимальное и минимальное (оба неотрицательные) собственные числа матрицы АтА.
С ростом V — числа обусловленности матрицы А, точностные характеристики решателя ухудшаются, что выдвигает постановку задачи о выборе активного управления объектом идентификации, доставляющем минимум V.
Точность формирователя. Исследование точности формирователя системы АХ = В рассматривается на примере интегрального метода производящих функций [6-8] применительно к задаче идентификации трех постоянных параметров , ж2 , ж3 уравнения моментов в модели короткопериодического движения летательного аппарата [6]:
(7)
$ + ж2 $ + ж3 а = и; ж4 (а — $) + ж5 а = и.
Экспериментальными данными служат записи процессов угловой скорости тангажа $(£), угла атаки а(£) и отклонения руля высоты и(£) на отрезке времени [0,Т]. В качестве производящих функций для поставленной задачи используются функции Эрмита ^ (г) (7 = 0,1, 2, 3) (рисунок) следующего вида:
Go = e-r2/2 ,Gi = -re-r2/2, G2 = (r2 - 1)e-r2/2, G3 = -(r3
6r)e
,-r2/2.
(8)
Одно из свойств функций Эрмита состоит в том, что dGj (r)/dr = = Gj+i(r) [9].
Функция Эрмита порядка j приближенно финитная: имеется "носитель" в виде отрезка [—rm, rm], вне которого она монотонна и lim Gj (r) = 0. При
этом функции Эрмита порядков, низших по сравнению с j, имеют аналогичные носители, но вложенные в [—rm, rm]. В частности, G3(rm) ~ 10-6 при rm = 5,6.
Согласуем T, rm, положив r = m(t — -T/2), m = 2rm/T.
Домножим все члены дифференциального уравнения моментов из модели (7) на функцию Эрмита Gj (r(t)) порядка j и почленно проинтегрируем по t Графики функций Эрмита на отрезке [0,Т]. Проблемный интеграл
х
T
Gj (r(t)) i(t)dt c недоступной функцией - после интегрирования
по частям сводится к интегралу с доступной функцией $:
т т
[ в3 (г(*)) = в,- (г(*)) $(*) т - т [ ву+1 (г(*))
Поскольку в у (±гт) « 0, постольку в у (г(£)) $(£)
T
m} ^ XJ > HWVAWJlDrVJ' Uj
T T
f Gj (r(t)) K?(t)dt w-m / Gj+1 (r(t)) i(t)dt
0, так что
Вместо дифференциального уравнения моментов получим интегральную модель в форме алгебраического уравнения
—т , ву+1^ Х1 + $, в^ Х2 + (а, в у) Хз = (м, в у),
T
где (Р,5) = у (г(^)) ^
о
Последовательно повторив данную процедуру с ву (г) (3 = 0,1, 2), в итоге получим матричное уравнение АХ = В, где
A=
/ -m(i,Go) (a,Go>\
-m/tf /г? ,GA (a,G1>
V -m(-(г?(a,G2> у/
(м,во) \ / Х1
в= | (м,в1) I , X = I Х2 I .
(м, в2^ \ Хз
Помехи измерительного тракта моделируются аддитивным наложением случайных процессов:
$(*) = $(*)+ &(*), «(¿) = «(¿)+ &(*), «(*) = М(*)+ &(*).
Здесь (£), г = 1, 2, 3, — независимые между собой, а также от измеряемых процессов центрированные помехи типа гауссовых белых шумов с известными первыми двумя моментами
0 пРи г = 3;
M [^i(t)] = 0, M (т)] =
S(t - т) пРи г = j.
Тогда формирователь выдаст матричные возмущения
/ —т (6,^) (6, С> \
ДЛ= -т <^1,С2> (6,^ (6, С >
V —т <6,Сз> <?1,С2> (^2,^2>) (9)
( (6,Со> " ДВ = I (6,^ V (?з,С2>
Вычисляя математические ожидания для ДА, ДВ, учтем
М [(&,С,- >] = (М &],С,- > = 0,
и получим статистическую несмещенность матричных возмущений (9):
М[ДЛ] = 0, М[ДВ] = 0.
Вычислим теперь дисперсию матричных возмущений (9), она будет отражать степень разброса ошибок:
Б [(&, С,->] = М [((&, С> — М [&, С3>])2] = М [((&, С>)2] = = М | / &(*) • С (г(*)) ^ I • I / &(*) • С (г(*)) ^
= M
т / т
[ G3 (r(t)) • ( / Ш • &(t) • Gj (r(t)) dr I dt
0
ч0
Т / Т
= 1 С3 (г(*)) • И М [&(*) • &(*)] • (г(*)) ¿Т I А =
00
т / т \ т
= I Сз (г(*)) I / ^ — Т)Сз (г(*)) ¿т I Л =Б& / (Сз (г(£)))2 Л. 0 0 0
т
Для вычисления интеграла У (Сз (г(£)))2 заменим £ на г =
= m(t — T/2) и найдем
т
0
(Gj (r(t)))2 dt = - [ (Gj (r))2 dr «- [ (Gj (r))2 dr.
Подставив вместо Сз (г) их формулы (8), получим, что интегрирование по частям будет сведено к вычислению интеграла Пуассона
m
— Г
m
p
/ e-x2dx =л/П; x
J —ж J
[ x2n+1e" "x2 dx =0.
J
2n „—x
1 • 3 • 5 ••• (2n - 1)
^/п, n = 1, 2, 3,...;
Приведем получившиеся числовые значения
\2
(Gj (r)) dx = Ajv^n,
1 л 3 . 87 где Ао = 1, А = -, А = -, Аз = —.
2 4 8
Результаты представлены матрицами дисперсий элементов ДА, ДБ
/
Яда =
/п m
2
^ m2
D5l--
2
3m2
87m2
D,
5i
D
52
D?1 • - D,4 • -
D
5i
8
52
D51 • - D5, • -
52
\
1 2
3
4 /
(10)
Ядв =
' п m
( D53 \
D?3 • 1
2
_ 3
\ ^ • -4)
По найденным матрицам дисперсий (10) можно прогнозировать априорные вероятностные границы разброса ошибок идентификации, привлекая формулы (5) и (6).
Обсуждение результатов. На примере матричного алгоритма идентификации динамических параметров модели (7) продемонстрирована конструктивность методологии получения априорных гарантирующих оценок погрешности решения в матричных алгоритмах, представляемых двухуровневой структурой формирователь-решатель. Можно сделать обобщающий вывод о стратегии исследований:
а) указать способы минимизации вариаций £А, £Б и, следовательно, минимизировать элементы дисперсионных матриц Ддл, ^дв из (10);
б) указать способы минимизации числа обусловленности матрицы А.
Заключение. О минимизации вариаций дА, дБ. В задачах идентификации приходится обращаться к первичной обработке экперимен-
тальных данных, в частности фильтровать и сглаживать зашумленные выходные процессы датчиков первичной информации. Очищение сигналов входа-выхода идентифицируемого линейного объекта от наложенных высокочастотных флуктуационных помех (с частотами выше w, с-1) возможно посредством пропускания этих сигналов через идентичные фильтры с передаточной функцией ——- (т ~ w-1). Наря-
т S + 1
ду с блоками "формирователь" и "решатель" в структуру матричного алгоритма вводится блок "регуляризатор" с параметром т. Корректность подобной регуляризации при идентификации асимптотически устойчивых объектов подтверждена математическим анализом и вычислительным экспериментом [14] и трактуется как разделение этапов фильтрации и идентификации.
О минимизации числа обусловленности матрицы A. Малость числа v = ||Ao|| • ЦА-1!! на специальных режимах движения идентифицируемого объекта — это задача планирования эксперимента. Данному вопросу должно быть уделено особое внимание. Проводившееся имитационное моделирование показало, что в зависимости от числа обусловленности, а оно на различных участках выборки исходных данных может широко изменяться от 10 до 107 и более, наследственная погрешность оценивания коэффициентов модели (7) адекватно пропорциональна числу обусловленности. При скачкообразном входном сигнале число обусловленности будет наименьшим на участке, соответствующем полному переходному процессу по координатам ? и а. Эти наблюдения подтверждают необходимость планирования эксперимента с учетом критерия максимума математической обусловленности проблемы обработки экспериментальных данных.
О единстве стратегических задач. Единство их выражается в форме неравенства (4) для знака дискриминанта корректности.
Дополнительные соображения. Численное воплощение конкретного алгоритма должно сопровождаться анализом распространения вычислительных погрешностей. Следует также считаться с присутствием малых параметров математической модели на фоне доминирующих, например демпфирующего коэффициента в модели корот-копериодического движения маневренного летательного аппарата [5], устойчивость оценки которого низкая, и не только по причине использования off-line метода идентификации, но и on-line метода, такого как адаптивный, описанный в работе [13].
ЛИТЕРАТУРА
1. Эйкхофф П. Основы идентификации систем управления. М.: Мир, 1975. 193 с.
2. Льюнг Л.Идентификация систем. Теория для пользователя / пер. с англ. М.: Наука, 1991. 432 с.
3. Иванов В.В. Методы вычислений на ЭВМ. Киев: Наукова думка, 1986. 584 с.
4. Форсайт Дж., Молер К. Численное решение систем линейных алгебраических уравнений. М.: Мир, 1969. 168 с.
5. Лебедев А.А., Чернобровкин Л.С. Динамика полета. М.: Машиностроение, 1978. 615 с.
6. Журавлев Ю.В. Идентификация линейных стационарных систем методом производящих функций // В кн.: Исследование, проектирование и расчет технических средств кибернетики, радио-систем и установок летательных аппаратов. под ред. В.И. Козлова и А.С. Протопопова. Тр. МАИ. Вып. 348. М.: Изд-во МАИ, 1976. С. 9-16.
7. Loeb J.M. and Cahen G.M. Extraction, a partir des enregistrements de mesures, des parameters dynamiques d'un systeme //Automatisme. No. 12, Dec. 1963. P. 479-486.
8. Tokaya K. The use of Hermite functions for systems identification // IEEE Trans., 1968, AC-13, August, No. 4.
9. Никифоров А.Ф., Уваров В.Б. Специальные функции математической физики. М.: Наука, 1984. 344 c.
10. Журавлев Ю.И., Гуревич И.Б. Распознавание образов и распознавание изображений // Распознавание, классификация, прогноз. 1989. Т. 2. № 5. С. 5-73.
11. Галанин М.П., Савенков Е.Б. Методы численного анализа математических моделей. М.: Изд-во МГТУ им. Н.Э. Баумана, 2010. 591 c.
12. Тихонов А.Н., Арсенин В.Я. Методы решения некорректных задач. М.: Наука, 1986. 288 c.
13. Журавлев Ю.В. Синтез адаптивного следящего идентификатора прямым методом Ляпунова // Вестник МГТУ им. Н.Э. Баумана. Сер. Естественные науки. 2010. № 2. С. 3-15.
14. Журавлев Ю.В. О разделимости задач фильтрации и идентификации линейных стационарных динамических систем // В кн.: Вопросы кибернетики и радиотехники / под ред. В.И. Козлова. Труды МАИ. Вып. 426. М.: Изд-во МАИ, 1977. С. 24-28.
REFERENCES
[1] Eykhoff P. System identification. London, Wiley, 1974. 253 p. (Russ. ed.: Eykkhoff P. Osnovy identifikatsii sistem upravleniya. Moscow, Mir Publ., 1975. 193 p.).
[2] Ljung L. System identification: Theory for the user. New York, Prentice Hall, 1987. 520 p. (Russ. ed.: L'yung L. Identifikatsiya sistem. Teoriya dlya pol'zovatelya. Moscow, Nauka Publ., 1991. 432 p.).
[3] Ivanov V.V. Metody vychisleniy na EVM [Methods of computer calculations]. Kiev, Naukova Dumka Publ., 1986. 584 p.
[4] Forsythe G.E., Moler C.B. Compuer solution of linear algebraic systems. Englewood Cliffs, Prentice Hall, 1967. 167 p. (Russ. ed.: Forsayt Dzh., Moler K. Chislennoe reshenie sistem lineynykh algebraicheskikh uravneniy. Moscow, Mir Publ., 1969. 168 p.).
[5] Lebedev A.A., Chernobrovkin L.S. Dinamika poleta [Flight dynamics]. Moscow, Mashinostroenie Publ., 1978. 615 p.
[6] Zhuravlev Yu.V. Identification of linear time-invariant systems by generating functions. In: Kozlov V.I., Protopopov A.S. Issledovanie, proektirovanie i raschet tekhnicheskikh sredstv kibernetiki, radio-sistem i ustanovok letatel'nykh apparatov [The research, design, and calculation of cybernetics hardware, radio and aircraft systems]. Tr. MAI [Trans. Moscow Aviat. Inst.], vol. 348. Moscow, MAI Publ., 1976. pp. 9-16 (in Russ.).
[7] Loeb J.M., Cahen G.M. Extraction, a partir des enregistrements de mesures, des parameters dynamiques d'un systeme. Automatisme, 1963, no. 12, pp. 479-486.
[8] Tokaya K. The use of Hermite functions for systems identification. IEEE Trans., 1968, vol. AC-13, no. 4.
[9] Nikiforov A.F., Uvarov V.B. Spetsial'nye funktsii matematicheskoy fiziki [Special functions of mathematical physics]. Moscow, Nauka Publ., 1984. 344 p.
[10] Zhuravlev Yu.I., Gurevich I.B. Artificial perception and image recognition. Raspoznavanie, Klassifikatsiya, Prognoz [Recognit. Classif. Predict.], 1989, vol. 2, no. 5, pp. 5-73 (in Russ.).
[11] Galanin M.P., Savenkov E.B. Metody chislennogo analiza matematicheskikh modeley [Methods of numerical analysis of mathematical models]. Moscow, MGTU im. N.E. Baumana Publ., 2010. 591 p.
[12] Tikhonov A.N., Arsenin V.Ya. Metody resheniya nekorrektnykh zadach [Methods for solving ill-posed problems]. Moscow, Nauka Publ., 1986. 288 p.
[13] Zhuravlev Yu.V. Synthesis of adaptive tracking identifier by the direct Lyapunov method. Vestn. Mosk. Gos. Tekh. Univ. im. N.E. Baumana, Estestv. Nauki [Herald of the Bauman Moscow State Tech. Univ., Nat. Sci.], 2010, no. 2, pp. 3-15 (in Russ.).
[14] Zhuravlev Yu.V. Separability of filtration and identification problems for linear timeinvariant dynamic systems. In.: Kozlov V.I. Voprosy kibernetiki i radiotekhniki [Problems of cybernetics and radio engineering]. Tr. MAI [Trans. Moscow Aviat. Inst.], vol. 426. Moscow, MAI Publ., 1977. pp. 24-28 (in Russ.).
Статья поступила в редакцию 3.04.2013
Юрий Васильевич Журавлев — старший преподаватель кафедры "Вычислительная математика и математическая физика "МГТУ им. Н.Э. Баумана. Автор более 20 научных работ в области моделирования, идентификации и управления динамическими системами.
МГТУ им. Н.Э. Баумана, Российская Федерация, 105005, Москва, 2-я Бауманская ул., д. 5.
Yu.V. Zhuravlev — senior teacher of "Computational Mathematics and Mathematical Physics" department of the Bauman Moscow State Technical University. Author of more than 20 publications in the field of simulation, identification and control of dynamical systems.
Bauman Moscow State Technical University, Vtoraya Baumanskaya ul., 5, Moscow, 105005 Russian Federation.