IDENTIFICATION OF DYNAMIC SYSTEM WITH DEAD ZONE OF THE ELASTIC ELEMENT AND DRY FRICTION ALONG CURVES
Abstract and keywords
Abstract (English):
The study objective is to increase the efficiency of identifying parameters of dynamic systems that have an elastic element with a dead zone and dry friction in their composition according to experimental data. The problem to which the paper is devoted is to develop an identification algorithm that allows defining the parameters of the frequency response function of a nonlinear dynamic system for all experimentally obtained points of the system frequency curve. The research methods are based on the theory of dynamical systems, in particular on harmonic linearization. The novelty of the work is in applying modules of modifying deviations of the curve points of the system model from the corresponding experimentally obtained system curve points as a measure of proximity for the system curves and the model of the sum of squares. The study result is an algorithm for identifying a nonlinear system, which is reduced to solving a system of equations linear to unknown parameters of the transfer function of the system model. Computing experiment showed that the errors in defining the values of the system parameters with a dead zone of an elastic element and a non-linearity of dry friction type are approximately equal to the measurement errors of the experimental readings of this system curve. The main conclusion of the study is that the developed algorithm should be used to solve practical problems of identifying nonlinear systems, which include an elastic element with a dead zone and dry friction.

Keywords:
dynamic system, dead zone of elastic element, dry friction, identification, frequency hodograph.
Text
Publication text (PDF): Read Download

Введение

 

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

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

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

 

 

a2xt+gsignx(t)+Fx(t)=qsinωt                              (1)

 

В этом уравнении F(x(t)) – зависимость силы, действующей на пружину от координаты x(t). Эта зависимость содержит зону нечувствительности и приведена на рис 1.

 

 

Рис. 1. Зависимость силы, действующей на пружину от координаты x(t)

Fig. 1. The dependence of the force acting on the spring on the x(t) coordinate

 

 

На рис. 1: k – коэффициент жесткости пружины, b – параметр, определяющий ширину зоны нечувствительности. Зависимость силы трения  от координаты   xt , определяющая сухое трение, приведена на рис. 2

 

 

Рис. 2. Зависимость силы трения   от координаты    xt ,

определяющая сухое трение

Fig. 2. The dependence of the friction force

on the coordinate x ̇(t) determining dry friction

 

На рис. 2: g - сила сухого трения.

 

Материалы, модели, эксперименты и методы.

 

Исходными данными для нахождения значений параметров   k, b, a2, g уравнения (1) служат экспериментально полученные амплитудно-частотная и фазо-частотная характеристики системы, представленные в виде годографа. Зависимости силы, действующей на пружину с зоной нечувствительности от координаты x(t) и зависимость силы трения g от координаты               xt   , определяющая сухое трение, можно заменить на зависимости, определяемые коэффициентами гармонической линеаризации [3,4]. При такой замене вынужденные колебания системы будут определяться уравнением

 

                             a2x(t)+a1(ω)x(t)+a0(ω)x(t)=qsinωt  ,                                         (2)

 

где, в соответствии с [2]

 

  a0ω=k-2kπ arcsinbAω   +bAω1-b2A(ω)2  ,                 (3)

 

a1(ω) = 4g / (πA(ω)ω).                                              (4)

 

 

Здесь A(ω) - амплитуда синусоидальной составляющей колебаний системы на частоте ω. Рассмотрим колебания системы, которые удовлетворяют условию 0<bAω<0.5  . С учетом допущения, что

 

                                  arcsinbAωbAω  ,    bAω1-b2A(ω)2   bAω  ,

 

 

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

 

a0(ω) = k - 4kb / (πA(ω))       (5)

с учетом обозначения         

s = 4kb / π                              (6)

вместо (5) можем записать

a0(ω) = k - s / A(ω).               (7)

С учетом обозначения

 

v = 4g / π                                (8)

вместо (4) можем записать

a1(ω) = 4v / (A(ω)ω).             (9)

 

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

 

 

 

W=1k- sAω-a2ω2+vA(ω)j  .                         (10)

 

 

Параметры k, s, a2, v частотной передаточной функции (10), можно получить, решив соответствующую задачу идентификации [5-8].

Функцию (10) можно переписать в альтернативном виде:

W(jω) = P(ω) + jQ(ω),

где P(ω) и jQ(ω) - вещественная и мнимая части частотной передаточной функции соответственно, определяемые следующими выражениями:

 

                                          P= k- sA(ω) -a2ω2   (k- sAω -a2ω2)2+(vAω)2 ,                       

                                                                                                                                 (11)

                                         Q=-vA(ω)   (k- sAω -a2ω2)2+(vAω)2 .                                                        

 

Амплитуду вынужденных колебаний динамической системы для частоты ω можно получить из соотношения A2(ω) = P2(ω) + Q2(ω) или с учетом (11) 

 

Aω=-R±R2-4TU2S .

 

Здесь R =  (k - a2 ω2)2,   T = -2(k - a2 ω2)s,    U = s2+ v2 – 1.   

Частотную передаточную функцию W(jω) можно представить в виде соответствующего годографа на комплексной плоскости.

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

 

 

P1 =P(ω1),..., Prexp =P(ωrexp),   Q1 = Q(ω1),..., Qrexp = Q(ωrexp)                    (12)

 

 

С учетом того, что вид передаточной функции известен заранее, соответствующий годограф будет принадлежать классу годографов, задаваемых моделью Wм(jω) с учетом (10)

 

Wм=1km- smAω-e2ω2+vmA(ω)j  .

или в виде

Wм(jω) = α/( γ + jδ),

где

 

                                   α = 1, γ = km -  sm /A(ω) - e2 ω2,   δ = vm /A(ω).                                                       

 

Идентификацию рассматриваемой системы можно свести к решению системы линейных уравнений относительно неизвестных параметров km,,   sme2,   vm    [9]

         Φe=u  ,                              (13)

vmrexp=-i=1rexpQiPi2+Qi2    .

где

 

 

       Φ=Φ00Φ01Φ02Φ10Φ11Φ12Φ20Φ21Φ22 ,                        e=kmsme2  ,       u=u0u1u2  .

 

Здесь элементы матрицы   Φ  и элементы вектора u    имеют значения

 

   Φ00i=1rexpPi2+Qi2  ,     Φ01  =-i=1rexpPi2+Qi2  ,     Φ02=  -i=1rexpPi2+Qi2ωi2 ,

 Φ10  = i=1rexpPi2+Qi2  ,     Φ11  = -rexp  ,      Φ12  =-i=1rexpPi2+Qi2ωi2 ,

   Φ20 = i=1rexpPi2+Qi2ωi2 , Φ21=-i=1rexpPi2+Qi2 ωi2  ,  Φ22 =-i=1rexpPi2++Qi2ωi4 ,

u0=i=1rexpPi  ,       u1=i=1rexpPiPi2+Qi2  ,           u2=i=1rexpPiωi2  . 

Решая систему уравнений (13), получим выражения для параметров  km,,   sme2   в виде формул Крамера

  

km, = H0 / H,   sm, = H1 / H,   e2, = H0 / H, где

 

            Η=Φ00Φ01Φ02Φ10Φ11Φ12Φ20Φ21Φ22

 

Η0=u0Φ01Φ02u1Φ11Φ12u2Φ21Φ22   ,  Η1=Φ00u0Φ02Φ10u1Φ12Φ20u2Φ22   ,  Η2=Φ00Φ01u0Φ10Φ11u1Φ20Φ21u2   ,

 

а для параметра vm

                                          vm=-1rexpi=1rexpQiPi2+Qi2  

С учетом обозначений (6) и (8) получим

                                bm, = π sm / (4km),              gm, = π vm / 4.

 

Вычислительный эксперимент по оценке ошибок в значениях параметров частотной передаточной функции системы, полученных с применением разработанного алгоритма, проведен для значений  k = 1,0;    b = 0,1; e2 = 1,0; g = 0,5.   Ошибки измерения отсчетов Pi, Qi  (12) экспериментально полученного годографа были выбраны случайно в соответствии с равномерным законом плотности распределения вероятностей в интервале [-0,1; 0,1]. Число вычислительных экспериментов rexp = 20. Проводилось rser=200 серий таких экспериментов в частотном диапазоне от 0.2 до 2.0 Гц. В сериях экспериментов вычислялись шибки errk = k - kmerrb = b - bmerra2 = a2 - e2errg = g - gm, определения параметров k, b, a2 , g  и среднеквадратические отклонения sko для случайных величин errk, errb, erra2, errg:

                                   sko_errk =0,0279;  sko_errb = 0,0234sko_errg =,.0150.

Также построены гистограммы для этих величин.  На рис. 3-5,   приведены примеры гистограмм для случайных величин errk, errb, errg.

 

14

 

 

Рис. 3. Гистограмма ошибки errk  определения параметра k

Fig. 3. Histogram of the error errk for determining the parameter k

 

Рис. 4. Гистограмма ошибки errb определения параметра b

Fig. 4. Histogram of the error errb of determining the parameter b

 

 

Рис. 5. Гистограмма ошибки errg определения параметра g

Fig. 5. Histogram of the error errg for determining the parameter g

Результаты

 

Проведенный вычислительный эксперимент показал, что ошибки определения значении параметров k, b, g     имеют близкий к нормальному закону распределения. На рис. 3-5 приведены соответствующие этим параметрам гистограммы ошибок. Для нормального закона распределения примерно 95% результатов определения значения параметра находится в доверительном интервале ± 2sko и примерно 99% результатов находится в доверительном интервале ±3sko [10]. Оценки погрешностей определения значений параметров k, b, g в вычислительном эксперименте в виде доверительных интервалов приведены в таблице.

 

 

Таблица

Доверительные интервалы результатов определения значений параметров

Table  

Confidence intervals of the results of determining the values of parameters

 

par=k

par=b

par=g

sko_errpar (par=k,b,g)

0.0279

0.0234

0.0150

± 2sko_errpar (par=k,b,g)

±0.0558

±0.0468

±0.0300

± 3sko_errpar (par=k,b,g)

±0.0737

±0.0702

±0.0450

 

 

Таким образом показано, что ошибки определения значения параметров k = 1,0; b = 0,1;  g = 0,5 в вычислительном эксперименте приближенно соответствуют интервалу ошибки измерений отсчетов годографа [-0,1; 0,1].

 

 

Обсуждение/Заключение

 

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

Несмотря на все более широкое внедрение в практику методов искусственного интеллекта, классические методы и алгоритмы анализа динамических систем остаются актуальными и развиваются, многие специалисты продолжают работы в этом направлении [11-13]. Наиболее перспективным по мнению авторов подходом при этом является гибридизация, которая с одной стороны предусматривает использование современных вычислительных возможностей по обработке «больших» экспериментальных данных с помощью нейросетевых [11,12], нечётких [13] и других моделей, а с другой предусматривает основу в виде классических методов и алгоритмов анализа таких систем [14]. Кроме того, существуют гибридные модели на основе теории игр, рассматривающие поведение агентов как динамическую систему в контексте решения социальных дилемм [15]. Предложенный в настоящей статье алгоритм является развитием классических методов анализа динамических систем и может применяться для идентификации реальных устройств в процессе их разработки, модификации и прототипирования, а также в качестве компонента той или иной гибридной модели, предназначенной, например, для решения социальных дилемм со множеством агентов. Предварительная идентификация параметров системы позволит сформировать обоснованные критерии исключения из массива накапливаемых данных заведомо ложных элементов, обусловленных различными сбоями и помехами при их измерении и записи.

References

1. Sevostyanov PA. A computer model of changes in the characteristics of a fibrous material in the production process. Textile Industry Technology. Proceedings of Higher Educational Institutions. 2016;1:361.

2. Meshcheryakov VN, Tolcheev VM. Mathematical model of the electromechanical system of the electric drive of the cold rolling mill. Electrotekhnika. 2015; 2(2):12-21.

3. Kulyabov DS, Korolkova AV, Veliyeva TR. Application of the harmonic linearization method to the study of the self-oscillatory mode of systems with control. Discrete and Continuous Models and Applied Computational Science. 2017;25(3).

4. Popov EP, Paltov IP. Approximate methods for the study of nonlinear automatic systems. Moscow: GIFML; 1960.

5. Leonov GA. On the method of harmonic linearization. Automation and Remote Control. 2009;5:65-75.

6. Spiridonov EG, Leshchenko MN. Development of algorithms for estimating parameters of dynamic systems of transport infrastructure facilities. Transport: Science, Education, Production (TRANSPORT-2021). Proceedings of the International Scientific and Practical Conference; Rostov n/A: 2021. p. 211-214.

7. Pavlov YuN, Nedashkovsky VM, Tikhomirova EA, Shavyrin IB. Method of harmonic linearization in the problem of identifying nonlinear dynamic systems. Science and Education. 2014;4: 382-394. doi: 10.7463/0414.0704613.

8. Shapkarin AV, Prosandeev AV, Kullo IG. Analysis of nonlinear automatic control systems by harmonic balance method in MATLAB. Caspian Journal: Management and High Technologies. 2013;1:077-085.

9. Boikov IV, Krivulin NP. Methods of identifying dynamic systems. Software systems: theory and applications. 2014;5 (23):79-96.

10. Wentzel ES. Probability theory: textbook for universities. 6th ed. Moscow: Visshaya Shkola; 1999.

11. Tiumentsev Y, Egorchev M. Neural network modeling and identification of dynamical systems. Academic Press; 2019.

12. Quaranta G, Lacarbonara W, Masri SF. A review on computational intelligence for identification of nonlinear dynamical systems. Nonlinear Dynamics. 2020;99(2):1709-1761.

13. Škrjanc I, Blažič S, Agamennoni O. Identification of dynamic systems with a robust interval fuzzy model. Automatica. 2005;41(2):327-332.

14. Semenov AD. A combined neural network system for regulating the linear density of the tape. Textile Industry Technology. Proceedings of Higher Educational Institutions. 2021;2:109-112.

15. Akiyama E, Kaneko K. Dynamical systems game theory II: A new approach to the problem of the social dilemma. Physica D: Nonlinear Phenomena. 2002;167(1-2):36-71.

Login or Create
* Forgot password?