Введение  Первичная разработка, модификация и прототипирование различных динамических устройств связаны с необходимостью проведения большого количества экспериментов. В ходе этих экспериментов характеристики устройств могут быть получены посредством изучения отклика на заранее заданные тестовые воздействия, применяемые к устройствам. Устройства, обладающие рядом заранее известных характеристик, объединяются в соответствующий класс. В настоящей работе рассматривается задача идентификации динамической системы, которой соответствует класс устройств, имеющих в своем составе упругие элементы с зоной нечувствительности и с сухим трением.Актуальность работы обусловлена необходимостью проведения экспериментальных исследований на этапе разработки, прототипирования и внедрения в производство различных устройств в машиностроении и на транспорте. В настоящее время в этих отраслях в нашей стране возрастает активность коллективов разработчиков во взаимодействии с производителями различного оборудования. Кроме того, в ряде случаев требуется доработка уже эксплуатирующегося оборудования. В частности, в хлопкопрядении к устройствам рассматриваемого класса относятся следящие приводы, которые управляют линейной плотностью выпускаемого полуфабриката (чесальной ленты) [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), получим частотную передаточную функцию вида:   Wjω=1k- sAω-a2ω2+vA(ω)j  .                         (10)   Параметры k, s, a2, v частотной передаточной функции (10), можно получить, решив соответствующую задачу идентификации [5-8]. Функцию (10) можно переписать в альтернативном виде:W(jω) = P(ω) + jQ(ω),где P(ω) и jQ(ω) - вещественная и мнимая части частотной передаточной функции соответственно, определяемые следующими выражениями:                                           Pjω= k- sA(ω) -a2ω2   (k- sAω -a2ω2)2+(vAω)2 ,                                                                                                                                                         (11)                                         Qjω=-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мjω=1km- smAω-e2ω2+vmA(ω)j  .или в видеWм(jω) = α/( γ + jδ),где                                    α = 1, γ = km -  sm /A(ω) - e2 ω2,   δ = vm /A(ω).                                                         Идентификацию рассматриваемой системы можно свести к решению системы линейных уравнений относительно неизвестных параметров km,,   sm,  e2,   vm    [9]         Φe=u  ,                              (13)vmrexp=-i=1rexpQiPi2+Qi2    .где         Φ=Φ00Φ01Φ02Φ10Φ11Φ12Φ20Φ21Φ22 ,                        e=kmsme2  ,       u=u0u1u2  . Здесь элементы матрицы   Φ  и элементы вектора u    имеют значения    Φ00=  i=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,,   sm,  e2   в виде формул Крамера    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 - km,  errb = b - bm,  erra2 = a2 - e2,  errg = g - gm, определения параметров k, b, a2 , g  и среднеквадратические отклонения sko для случайных величин errk, errb, erra2, errg:                                    sko_errk =0,0279;  sko_errb = 0,0234;  sko_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 определения параметра bFig. 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=kpar=bpar=gsko_errpar (par=k,b,g)0.02790.02340.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]. Предложенный в настоящей статье алгоритм является развитием классических методов анализа динамических систем и может применяться для идентификации реальных устройств в процессе их разработки, модификации и прототипирования, а также в качестве компонента той или иной гибридной модели, предназначенной, например, для решения социальных дилемм со множеством агентов. Предварительная идентификация параметров системы позволит сформировать обоснованные критерии исключения из массива накапливаемых данных заведомо ложных элементов, обусловленных различными сбоями и помехами при их измерении и записи.