CC BY
14
Прикладные задачи
^^^^^^^^^^»нелинейной теории колебаний и вслн
УДК 537.86
ВОССТАНОВЛЕНИЕ ПО ВРЕМЕННЫМ РЯДАМ АРХИТЕКТУРЫ СВЯЗЕЙ И ПАРАМЕТРОВ ЭЛЕМЕНТОВ В АНСАМБЛЯХ СВЯЗАННЫХ ОСЦИЛЛЯТОРОВ С ЗАДЕРЖКОЙ
И. В. Сысоев1'2, Д. Д. Кульминский2'1, В. И. Пономаренко2'1, М.Д. Прохоров2
1 Национальный исследовательский Саратовский государственный университет им. Н.Г. Чернышевского Россия, 410012 Саратов, Астраханская, 83 2Саратовский филиал Института радиотехники и электроники им. В. А. Котельникова РАН Россия, 410019 Саратов, Зеленая, 38 E-mail: ivssci@gmail.com; kulminskydd@gmail.com; ponomarenkovi@gmail.com; mdprokhorov@yandex.ru
Цель. Предложить новый подход к восстановлению архитектуры связей и параметров элементов в ансамблях связанных осцилляторов, описываемых дифференциальными уравнениями первого порядка с запаздыванием, по временным рядам их колебаний.
Метод. Метод основан на минимизации целевой функции, характеризующей расстояние между точками реконструируемой нелинейной функции данного элемента, и разделении восстановленных коэффициентов связи на значимые и незначимые. Минимизация целевой функции осуществляется методом наименьших квадратов. Время запаздывания определяется как соответствующее минимуму целевой функции по всем пробным временам запаздывания.
Результаты. Эффективность предложенного метода продемонстрирована в численном эксперименте на примере хаотических временных рядов ансамбля, состоящего из диффузионно связанных неидентичных уравнений Маккея-Гласса в присутствии шума, а также в натурном эксперименте на примере временных рядов резистивно связанных радиотехнических генераторов с запаздывающей обратной связью. Метод обеспечивает более высокую, чем ранее предложенные подходы, вычислительную эффективность за счёт использования неитерационных алгоритмов минимизации целевой функции и отбора значимых коэффициентов. При этом оценки коэффициентов связи и параметра инерционности являются несмещёнными.
Обсуждение. Метод может быть полезен для восстановления параметров элементов и архитектуры связей в системах различной природы: радиотехнических, биологических и иных, описываемых уравнениями первого порядка с запаздыванием.
Ключевые слова: Анализ временных рядов, реконструкция уравнений, ансамбли осцилляторов, системы с запаздыванием. DOI: 10.18500/0869-6632-2016-24-3-21-37
Ссылка на статью: Сысоев И.В., Кульминский Д.Д., Пономаренко В.И., Прохоров М.Д. Восстановление по временным рядам архитектуры связей и параметров элементов в ансамблях связанных осцилляторов с задержкой // Известия вузов. Прикладная нелинейная динамика. 2016. Т. 24, № 3. С. 21-37.
Введение
В последние годы большое внимание исследователей привлекает задача определения структуры и количественной оценки связей в ансамблях, состоящих из большого числа взаимодействующих между собой элементов, по временным рядам их колебаний. Важность этой задачи определяется тем, что архитектура и интенсивность связей во многом определяют особенности коллективной динамики элементов ансамбля и возможность их синхронизации [1-3]. Для реконструкции связей в многоэлементных системах используются методы, основанные на моделировании фазовой динамики [4-6], анализе причинности по Грейнджеру [7, 8], методы адаптивного управления [9-11] и другие методы [12-14]. Однако в большинстве описанных случаев элементы исследуемых ансамблей либо не имеют запаздывающей обратной связи, либо ее величина предполагается известной. Вместе с тем, системы с запаздыванием чрезвычайно широко распространены в природе и технике. Их повсеместность обусловлена такими фундаментальными свойствами, как конечная скорость распространения сигнала и наличие запаздывающей обратной связи, присущими многим физическим, химическим, климатическим и биологическим системам и процессам [15,16]. Ансамбли, состоящие из уравнений с запаздыванием, широко используются для моделирования и описания процессов в различных многоэлементных системах с задержкой [17-19].
Методы, позволяющие восстановить архитектуру и величину связей и одновременно оценить собственные параметры элементов в ансамблях систем с запаздыванием, были предложены в [20,21]. Наряду с рядом достоинств, эти методы имеют и недостатки. Например, для реализации метода [20] требуется обратимость функций, описывающих собственную динамику элементов ансамбля, отсутствие шума и задание стартовых догадок для времен запаздывания вблизи их истинных значений, что не всегда оказывается возможным. В методе, предложенном нами в [21], для восстановления времени запаздывания элементов используется отдельная процедура [22-24], а использование итерационного алгоритма для реконструкции архитектуры связей приводит к относительно большому времени работы, при этом полученный результат может зависеть от стартовых догадок архитектуры связей.
В настоящей работе предлагается новый метод восстановления ансамблей, состоящих из систем с запаздыванием, свободный от указанных недостатков. Метод основан на минимизации целевой функции для каждого элемента ансамбля, характеризующей расстояние между точками реконструируемой нелинейной функции, отсортированными по величине абсциссы, и использовании различных алгоритмов для разделения восстановленных коэффициентов связи на значимые и незначимые.
1. Метод
Рассмотрим ансамбль, состоящий из диффузионно связанных систем с запаздыванием, каждая из которых описывается уравнением следующего вида:
Б
егХг(г) = -хг(г) + f (хг(г - т^) + ^ к^ (хз (г) - х^г)), (1)
3 = 1,3 =
где г = 1... Б, Б - число элементов в ансамбле; ег - параметр, характеризующий инерционные свойства г-го элемента ансамбля; т% - время запаздывания; /\ - нелинейная функция; кг,з - коэффициенты связи, характеризующие воздействие ] -го элемента на г-й. В наиболее общем случае между любыми двумя элементами ансамбля существует взаимная связь.
Пусть у нас имеются временные ряды х* = {хг}1-!=1 длиною N точек всех Б осцилляторов, измеренные с шагом выборки АЬ. Пусть также все функции ¡г непрерывны, а длина временного ряда N достаточна для того, чтобы даже участкам самого быстрого на всём отрезке [шт(хг); шах(хг)} изменения ¡г(хг) во временном ряде {х}^ соответствовало несколько десятков точек.
Введем дискретное время запаздывания вг = тг/АЬ и перепишем уравнение (1) в виде
Б
¡г (хг(п)) = £гХг(н + в г) + хг(п + в г) — ^ кг,з (хз (п + вг) - Хг(п + вг)) , (2)
3=1,3=г
где п = 1.. .п — вг. Величины х*(п + вг) могут быть найдены численным дифференцированием исходного ряда. Для этого необходимо, чтобы шаг выборки АЬ был достаточно мал, чтобы позволить разрешить все существенные временные масштабы. Для каждого осциллятора отсортируем значения хг(п) по возрастанию, обозначив такую сортировку как преобразование Р, сопоставляющее точке с номером п в исходном ряде точку с номером Q(xг,n) в отсортированном ряде. Обратное преобразование, сопоставляющее точке с номером Q(xг, п) в отсортированном ряде точку с номером п в исходном, обозначим Q-1. Тогда п = Q-1(Q(n)). Для краткости обозначений здесь и далее зависимость Q и Q-1 от хг указывать не будем.
Пусть некоторая точка имеет номер п в исходном ряде и номер Q(n) в отсортированном ряде. Тогда её сосед справа в отсортированном ряде будет иметь номер Q(n) + 1, а в исходном ряде его номер будет рп = Q-1(Q(n) + 1), причем номера п и рп в общем случае не будут близки. Так как точки с номерами Q(n) и Q(n) + 1 являются соседними в отсортированном ряде, значения динамической переменной в этих точках будут близки. А значит, будут близки и значения функции ¡г от этих переменных поскольку все функции ¡г, как уже было сказано выше, изначально предполагаются непрерывными. Обозначим разность значений функции ¡г в этих точках как Ьг (п)
Ьг(п) = ¡г (хг(рп)) — ¡г (хг(п)). (3)
Используя уравнение (2), запишем уравнение (3) как
Б
Ьг(п) =(хг(рп + вг) — ^ кг ,з{х^(рп + вг) — хг(рп + вг)) + егхг(рп + в*))
з=1, з=г
Б
— (хг(п + вг) — ^ кг,з {х^(п + вг) — хг(п + вг)) + егхг(п + вг)). з=1,з=г
(4)
Введём новые обозначения и перепишем выражение (4) в следующем виде:
Б
Ъг(п) = АХг(п) - Е кг,] (АХ] (п) - АХг(п)) + £гАхг(п), (5)
АХг(п) = Хг(рп + бг) - Хг(п + бг), (6)
АХг(п) = Хг(рп + бг) - Хг(п + бг). (7)
Обозначим через Ь2 сумму Ъ2(п)
N-вi-l N-вi-l
<2,
— бi — L N —бi — L,
Ь = Е Ь2(п)= Е (АХг(п)-п=1 п=1 ^
Б \2 У] кг,з (АХ](п) - АХг(п)) - (-ег)АХг(п)) .
(8)
3=1,3=г
Величину Ь2 можно рассматривать как функцию от параметров бг, кг,з и ег, величина которых заранее не известна. При правильном выборе этих параметров Ь2 будет меньше, чем при ошибочном. Это объясняется тем, что при неправильном выборе вг, кг,з и ег расстояния (4) не будут малы даже для соседних точек в отсортированном ряде.
Заметим, что предложенная мера в некотором смысле характеризует длину восстановленной нелинейной функции, а именно, представляет собою сумму квадратов только вертикальных компонент расстояний между точками нелинейной функции, в то время как горизонтальные компоненты (Хг(рп) - Хг(п)) не поддаются оптимизации по параметрам бг, кг,з и ег и потому не включены в целевую функцию (8).
Поскольку (6) и (7) зависят только от вг, ег и кг,з при данном фиксированном г и не зависят от иных бт, ет или кт,з, для которых т = г, минимизацию целевой функции (8) возможно проводить отдельно для каждого осциллятора. При этом фактически решается задача реконструкции неавтономной системы по векторному ряду {Хг(п + вг), Хг(п), Хг(п + бг)^—6, одна компонента которого измеряется, другая восстанавливается методом задержек, а третья - методом дифференцирования. Временные ряды внешнего воздействия вычисляются явно по известным рядам переменных {Хт(п)}^^=в.+1. Таким образом, нет необходимости восстанавливать фазовое пространство большой размерности, пропорциональной числу осцилляторов в ансамбле.
Если зафиксировать вг, то задачу минимизации (8) можно рассматривать как линейную задачу о наименьших квадратах, где кг,з и (-ег) суть искомые коэффициенты (их всего В штук для каждого г-го элемента), АХг(п) - аппроксимируемые величины, а матрица значений базисных функций состоит из (Ах3-(п) - АХг(п)) и АХг(п). В такой постановке Ь2 представляет собою целевую функцию. Поиск её экстремума является стандартной задачей, которая может быть решена нерекурсивно.
Поскольку время запаздывания бг заранее не известно, минимизацию целевой функции (8) можно провести для различных пробных дискретных времен запазды-
вания вг, перебираемых из некоторого интервала. Минимум зависимости Ь2(т'г), где т'г = в'АЬ, будет наблюдаться при истинном времени запаздывания т*.
С увеличением числа N точек во временном ряде хг число членов суммы (8) будет расти пропорционально N. Вместе с тем, с ростом N будут уменьшаться расстояния между точками в отсортированном ряде и, как следствие, будут уменьшаться величины |6г(п)|. В среднем это уменьшение пропорционально 1/N. То есть каждый член Ь"1(п) суммы (8) будет убывать пропорционально 1/N с ростом N. Следовательно, Ь2 — 0 при N -ж. Значит, при N — ж предложенный метод является асимптотически точным, а полученные с его помощью оценки параметров являются асимптотически несмещенными.
Предложенный подход сформулирован для систем, описываемых уравнениями с запаздыванием, содержащими только непрерывные нелинейные функции. Однако если такие функции будут содержать на всём отрезке изменений аргумента небольшое число разрывов первого рода и сумма квадратов скачков функции будет много меньше общей суммы квадратов вертикальных расстояний (8), метод, по всей видимости, сохранит свою работоспособность при конечном N, несколько потеряв в точности. При этом свойства асимптотической несмещённости оценок будут утеряны.
Предложенный нами алгоритм описан для общего случая, при котором между любыми двумя элементами ансамбля существует двунаправленная связь, что редко бывает на практике. Если в действительности воздействие ]-го элемента на г-й отсутствует, то соответствующий коэффициент связи кг з в модельном уравнении (1) следует положить равным нулю. Однако в результате описанной выше процедуры реконструкции мы всегда получаем для каждого элемента ансамбля Б — 1 ненулевых коэффициентов связи кг 3; некоторые из них при отсутствии соответствующих связей являются лишними. Разделить восстановленные коэффициенты связи на значимые и незначимые и, таким образом, отбросить лишние связи можно с помощью метода К-средних [25].
Проведем для этого кластеризацию восстановленных коэффициентов к'г 3 в одномерном пространстве, поделив их на 2 кластера: значимых и незначимых коэффициентов. В качестве начальных положений центров кластеров зададим максимальное и минимальное из восстановленных значений к'г 3. Так как в общем случае значимые коэффициенты по абсолютной величине много больше незначимых, кластеризацию удобно проводить в логарифмическом масштабе. Определив незначимые коэффициенты связи, положим их равными нулю и повторно восстановим все кг з для повышения точности реконструкции.
Такой подход позволяет восстановить архитектуру связей в ансамбле. Следует отметить, что рассмотренный метод имеет существенно более высокое быстродействие по сравнению с методом реконструкции модельных уравнений элементов ансамбля и диагностики значимости связей с помощью последовательного пробного исключения или добавления коэффициентов связи в модель, предложенном нами в [21], поскольку в отличие от итерационного метода [21] реконструкция параметров проводится лишь дважды.
Подход к определению лишних связей, основанный на методе К-средних, хорошо работает в отсутствие шумов и соизмеримом количестве как значимых, так и незначимых коэффициентов связи. Однако при невыполнении этих условий границы
соседних кластеров оказываются близки друг к другу. В результате точность метода снижается, и он может находить ложные связи или пропускать часть имеющихся связей. В таких случаях для точной реконструкции архитектуры связей требуется использовать более двух кластеров в методе К-средних. К сожалению, при анализе экспериментальных данных, когда число связей между элементами ансамбля заранее не известно, дать точные рекомендации по выбору количества кластеров оказывается затруднительно.
Для разделения восстановленных коэффициенты связи кгз полной модели на значимые и незначимые можно использовать другой подход. Отсортируем все коэффициенты кг,з по абсолютной величине в порядке убывания. Затем рассчитаем величину
Б
Л = Е Ь2/(М - бг - 1) (9)
г=1
при введении в ансамбль, состоящий из В модельных уравнений (1), только одного коэффициента связи, самого большого по абсолютной величине. Нормировка компонентов, входящих в Л, на величину (К - бг - 1) необходима, поскольку иначе разные осцилляторы внесут разный вклад в меру (9) из-за отличающихся времен запаздывания бг и, как следствие, разного числа слагаемых в формуле для целевой функции (8). Далее будем добавлять в модель по одному коэффициенту связи в порядке уменьшения их абсолютного значения и снова рассчитывать Л. Наконец, построим зависимость Л от числа М коэффициентов связи в модели, где М = 1,...,В(В - 1).
По мере учета в модели все большего числа реально существующих связей она будет становиться все более точной. При этом величины (5), характеризующие расстояние между точками реконструируемых нелинейных функций элементов, а вместе с ними и величины (8) и (9) будут уменьшаться. Уменьшение величины Л на графике Л(М) практически остановится при М = Мгеа\, где Мгеа\ - число реальных связей. Последующее добавление в модель незначимых коэффициентов, соответствующих отсутствующим связям, почти не повлияет на величину Л. При этом, для М ^ Мгеа\ зависимость Л(М) будет оставаться почти постоянной. Графики, наглядно иллюстрирующие такой подход, приведены в следующем разделе. Отметим, что этот метод требует существенно больше вычислительных затрат, чем способ разделения восстановленных коэффициентов связи на значимые и незначимые на основе метода К-средних.
На практике для контроля точности восстановления априорно неизвестной архитектуры связей в ансамбле можно использовать оба рассмотренных подхода, которые в идеале должны дать одинаковые результаты.
2. Восстановление ансамбля, состоящего из связанных уравнений Маккея-Гласса
В качестве первого примера восстановим параметры элементов и архитектуру связей в ансамбле диффузионно связанных систем Маккея-Гласса [26], описываемых
Рис. 1. а - архитектура связей в ансамбле из 20 элементов вида (10); б - временной ряд переменной х1(Ь) ансамбля из уравнений Маккея-Гласса в присутствии нормального белого шума со среднеквадратичным отклонением ^ = 0.004
уравнением (1) с функцией
/г{Хг(г — Тг)) =
агхг(~^ Тг) Ъ (1+ Хг!°(* - Тг))
(10)
и £г = 1/Ъг. Уравнение Маккея-Гласса, описывающее процесс выработки организмом красных кровяных клеток, является эталонным уравнением с запаздыванием, широко используемым при численных исследованиях систем с задержкой. На рис. 1, а приведена архитектура случайно выбранных связей в ансамбле из 20 элементов. Из 380 возможных связей между элементами ансамбля имеется 60 связей, среди которых есть как однонаправленные, так и взаимные. Все элементы ансамбля являются неидентичными. Их параметры принимают случайные значения в следующих интервалах: т € [250; 400], £ € [7.5; 12.5], ц € [0.20; 0.25], кг^ € [0.02; 0.06]. При этом все элементы колеблются хаотически. Длина временных рядов N = 104 при шаге выборки ДЬ = 0.5. Рассмотрим случай отсутствия шума и случай, когда к временному ряду каждого элемента добавлен некоррелированный нормальный шум ^¿(¿) с нулевым средним и среднеквадратичным отклонением аг = 0.004. На рис. 1, б приведен зашумленный хаотический временной ряд колебаний первого элемента при Т! = 263, £! = 12.32, (! = 0.218, км = 0.0475, кц5 = 0.0294, кМ8 = 0.0292.
На рис. 2 приведены зависимости Ь2^) для всех 20 элементов ансамбля для случаев отсутствия и присутствия шумов. Глобальные минимумы всех Ь2^) наблюдаются в точности при истинных временах запаздывания. Величины Ь2 нормированы на величину N — 9г — 1) для удобства сравнения. Отметим, что для оценки временного ряда производных {хг(п)} временному ряду {хг(п)} использовалась аппроксимация со сглаживанием параболой.
При реконструкции модельного уравнения (1) для каждого элемента ансамбля получаются 19 ненулевых коэффициентов связи к'^, часть из которых являются лишними. Лишние коэффициенты по модулю гораздо на несколько порядков меньше действительных и поэтому их можно выявить с помощью кластеризации в логариф-
Рис. 2. Зависимости для каждого из 20 элементов ансамбля из уравнений Маккея-Гласса в
отсутствие шума (а) и в присутствии нормального шума с = 0.004 (б)
гтта|—I 111111П|—I 11 п|||||—I 11 п|||||—I I I п|||||—I 11 п|||||—I 11 п|||||
а,— 0
О
а,= 0.04
© ©
шщ|..................................................................
10"7 10"6 10"5 10"4 10"3 10-2 I и
р
Рис. 3. Распределение значений модулей оценок коэффициентов связи для всех элементов ансамбля из уравнений Маккея-Гласса в отсутствие шума (сверху) и в присутствии шума (снизу). Значимые коэффициенты показаны крестиками, расположенными внутри окружности, а незначимые коэффициенты показаны точками, расположенными внутри прямоугольника
мическом масштабе, как это показано на рис. 3. Видно, что модули всех значений к'г з хорошо делятся на 2 кластера, состоящие из действительных (справа) и лишних (слева) коэффициентов. Наличие шума приводит к сближению кластеров. Очевидно, что, начиная с некоторого критического уровня шума, метод начинает выдавать ошибки, пропуская часть имеющихся связей или находя ложные.
Члены суммы в (5), соответствующие лишним коэффициентам связи, были удалены из модели и реконструкция была проведена повторно уже без них.
В результате архитектура связей в ансамбле оказалась восстановленной правильно, в точном соответствии с рис. 1, а, как для случая отсутствия шумов, так и при их наличии. В отсутствие шума точность реконструкции коэффициентов модели составила 4 значащих цифры. Неидеальная точность в таком случае объясняется конечной длиной ряда и необходимостью численной оценки производной. При наличии шума значения коэффициентов
также были восстановлены с удовлетворительной точностью. В частности, для первого осциллятора они составили к[ 4 = 0.0467 при истинном значении к\,4 = 0.0475, к'1 15 = 0.0288 при истинном значении к1 ,15 = 0.0294, к'1 18 = 0.0287 при истинном значении к1}18 = 0.0292.
Чтобы охарактеризовать точность восстановления коэффициентов связи в среднем по ансамблю, была рассчитана средняя относительная погрешность коэффициентов по формуле (11)
к - - кг- ■ (11)
* = < ^\) - •
где усреднение проводилось только для действительных коэффициентов. Величина Ек составила 0.023, колеблясь от 0.002 до 0.058. Параметры инерционности были восстановлены с ещё большею точностью, так, для первого осциллятора е1 = 12.29 при истинном значении в1 = 12.32, аналогичным образом рассчитанная средняя относительная погрешность реконструкции составила 0.007 от абсолютной величины. Большая точность реконструкции параметров е% является, скорее всего, следствием того, что они вносят более уникальный вклад в целевую функцию (8), будучи до-множены на базисные функции АЖ%(п), в то время как параметры связи к%, - стоят перед однотипными базисными функциями вида (Аж- (п) — Аж%(п)), вклад которых может частично взаимно компенсироваться, особенно при частичной синхронизации в ансамбле.
На рис. 4, а приведена восстановленная нелинейная функция /1 первого элемента для случая отсутствия шума. Она очень точно совпадает с истинной функцией уравнения Маккея-Гласса (настолько, что на одном графике их невозмоно различить в использованном масштабе), поэтому приводить исходную функцию мы не стали. На рис. 4, б приведена та же нелинейная функция /1, построенная при указанных выше восстановленных значениях и при наличии шума. Точность восстановления параметров и нелинейных функций для остальных элементов ансамбля примерно такая же.
Очевидно, что лишние и действительные связи вносят различный по абсолютной величине вклад в целевую функцию (8). На этом основании нами был апробиро-
Рис. 4. Восстановленная функция /1(х1) в отсутствие шума (а) и в присутствии нормального шума со среднеквадратичным отклонением о = 0.004 (б)
ван альтернативным кластеризации методом К-средних алгоритм разделения действительных и лишних коэффициентов среди всех оценок к'г 3, основанный на построении зависимости (9) (рис. 5). При М = Мгеа! величина Л(М) практически достигает своего минимального значения и с увеличением М меняется очень слабо (менее чем на 1%), что видно на графике. В отсутствие шума Л(Мгеа}) очень близка к нулю и составляет 1.3 • 10-5, в то время как потеря даже одного коэффициента приводит к её резкому росту: Л(Мгеа[ — 1) = 55.7 • 10"5, то есть примерно в 43 раза. Таким образом, оба предложенных подхода дали в рассмотренном примере одинаковый результат и позволили точно реконструировать архитектуру связей.
Рис. 5. Зависимости Л(М) для ансамбля из уравнений Маккея-Гласса в отсутствие шума (серый цвет) и в присутствии нормального шума со среднеквадратичным отклонением = 0.004 (чёрный цвет)
3. Восстановление цепочки экспериментальных радиотехнических генераторов с запаздывающей обратной связью
Мы применили метод к экспериментальным временным рядам цепочки, состоящей из В = 10 однонаправленно связанных радиотехнических генераторов с запаздывающей обратной связью. Каждый генератор представляет собой кольцевую систему, состоящую из линии задержки, нелинейного элемента и низкочастотного КО-фильтра первого порядка (рис. 6, a). Нелинейные элементы и линии задержки генераторов выполнены на микроконтроллерах, а фильтры - на аналоговых элементах. Аналоговые и цифровые элементы схемы сопрягались с помощью аналого-цифровых и цифро-аналоговых преобразователей. Связь генераторов осуществлялась с помощью резисторов К&. Модельное уравнение, которым описывается 1-й элемент такой цепочки, имеет вид
Б
КгОгШ = —Уг(1) + ¡г{Уг(1 — Тг)) + ^ ^ (У — У^) , (12)
з=1, з=г
где Уг(Ь) и Уг(Ь — т*) суть напряжения на входе и выходе линии задержки; т* - время запаздывания; Кг и Ог - сопротивление и ёмкость элементов фильтра; ¡г - передаточная характеристика нелинейного элемента. Уравнение (12) сводится к виду (1), если положить ег = КгОг.
Цепочка состоит из неидентичных элементов, параметры которых принимают значения в следующих интервалах: тг € [2.50; 4.75] мс, ег € [0.203; 0.536] мс, кг,г+1 € [0.10; 0.23], кг,з = 0 = г + 1. Все нелинейные элементы имели квадратичную передаточную характеристику. Хаотические сигналы Уг(Ь) записывались с помощью 10-канального аналого-цифрового преобразователя с частотой выборки ¡8 = 100 кГц. На рис. 6, б приведен фрагмент временной реализации сигнала У1{Ь) в первом генераторе, имеющем параметры т1 = 2.5 мс, е1 = 0.203 мс, к1, 2 = К1/КС1 = 0.21.
На рис. 7 приведены зависимости Ь1(х'%) для всех 10 элементов цепочки. Глобальные минимумы девяти зависимостей наблюдаются при истинных временах запаздывания генераторов. Лишь для одного из генераторов минимум оказался смещён относительно истинного времени запаздывание на время, равное одному интервалу выборки.
а
АЦП Линия Нелинейный ЦАП
задержки —* элемент —>
К
С\
Рис. 6. a - блок-схема кольцевого генератора с запаздывающей обратной связью, где «АЦП» -аналого-цифровой преобразователь, «ЦАП» - цифро-аналоговый преобразователь. б - хаотическая временная реализация У1 (£) первого генератора
Рис. 7. Зависимости Ь2 (х^) для всех экспериментальных генераторов в цепочке
Для каждого генератора цепочки мы получаем 9 восстановленных коэффициентов связи к з. Из этих 9 коэффициентов для 9 элементов цепочки только один коэффициент является значимым, а для элемента на краю цепочки, на который не действуют другие генераторы, вообще не должно быть к\ ^. Поскольку число действительных коэффициентов на порядок меньше числа лишних коэффициентов, использование метода К-средних для разделения восстановленных коэффициентов связи на значимые и незначимые оказывается не столь эффективным, как в рассмотренном выше примере. Для успеха алгоритма потребовалось делить всё множество полученных ^ на четыре кластера, считая значимыми только к[ ^, принадлежащие верхнему из четырёх. При использовании двух кластеров в архитектуре связей оставались несколько лишних коэффициентов связи.
Метод, основанный на использовании меры (9), оказался более эффективным. На рис. 8, а построена зависимость Л(М). Видно, что при М = 9 зависимость Л(М) практически достигает минимума и с дальнейшим увеличением М почти не меняется. Оставив только 9 самых больших по модулю восстановленных коэффициентов связи и удалив все остальные незначимые, удалось получить оценки коэффициентов связи и параметров инерционности, близкие к номиналам. Например, для первого генератора восстановленные параметры инерционности и связи имеют следующие значения: е\ = 0.204 мс, к' 2 = 0.22. Восстановленная при этих значениях нелинейная функция приведена на рис. 8 б. Она достаточно хорошо совпадает с истинной передаточной характеристикой нелинейного элемента первого генератора. Аналогичные результаты были получены для остальных элементов.
Л 0.08 0.06 0.04 0.02
0
0 20.0 40.0 60.0 80.0 М 0 0.8 1.6 2.4 х.,В
а о
Рис. 8. Зависимость Л(М) для цепочки экспериментальных генераторов с запаздывающей обратной связью (а). Восстановленная нелинейная функция первого генератора /1X1) (б)
Заключение
Нами предложен новый эффективный метод, позволяющий определить значения параметров элементов и архитектуру связей в ансамблях связанных систем, описываемых дифференциальными уравнениями с запаздыванием, по временным рядам их колебаний. Метод основан на минимизации методом наименьших квадра-
тов целевой функции, характеризующей расстояние между точками реконструируемой нелинейной функции, отсортированными по величине абсциссы, для каждого элемента ансамбля. Предложенный подход позволяет с высокой точностью восстановить времена запаздывания, параметры инерционности, нелинейные функции и коэффициенты связи всех элементов ансамбля. Отказ от использованного нами ранее итерационного алгоритма для реконструкции архитектуры связей [21] позволяет на порядок увеличить скорость работы метода.
Для разделения восстановленных коэффициентов связи на значимые и незначимые предложено использовать два алгоритма. Один из них основан на кластеризации восстановленных коэффициентов связи с помощью метода K-средних, а другой - на построении зависимости суммы целевых функций всех элементов от количества коэффициентов связи в модельных уравнениях.
Предложенный подход можно применять к ансамблям, состоящим из неидентичных систем с запаздыванием с произвольным числом однонаправленных и взаимных связей между ними. Эффективность метода продемонстрирована на примере хаотических временных рядов ансамблей диффузионно связанных модельных систем Маккея-Гласса, в том числе при наличии шума, а также на примере хаотических экспериментальных временных рядов резистивно связанных радиотехнических генераторов с запаздывающей обратной связью. Предложенный метод может быть также применён для реконструкции ансамблей с другими типами связи элементов, например, для ансамблей, состоящих из систем с запаздыванием, связанных через производную или общим полем.
Работа выполнена при поддержке РФФИ, грант № 16-02-00091.
Библиографический список
1. Afraimovich VS., Nekorkin V.I., Osipov G.V., Shalfeev V.D. Stability, Structures, and Chaos in Nonlinear Synchronization Networks. Singapore: World Scientific, 1995.
2. Пиковский А., Розенблюм М., Куртс Ю. Синхронизация: Фундаментальное нелинейное явление. М: Техносфера, 2003. 496 c.
3. Boccaletti S., Latora V., Moreno Y., Chavez M., Hwang D.U. // Phys. Rep. 2006. Vol. 424. P. 175.
4. Безручко Б.П., Смирнов Д.А. Математическое моделирование и хаотические временные ряды. Саратов: ГосУНЦ «Колледж», 2005.
5. Timme M. Revealing network connectivity from response dynamics // Phys. Rev. Lett. 2007. Vol. 98. 224101.
6. Smirnov D.A., Bezruchko B.P. Detection of couplings in ensembles of stochastic oscillators // Phys. Rev. E. 2009. Vol. 79. 046204.
7. Kammski M., Ding M., Truccolo W.A., Bressler S.L. Evaluating causal relations in neural systems: Granger causality, directed transfer function and statistical assessment of significance // Biol. Cybern. 2001. Vol. 85. P. 145.
8. Sysoev I.V., Sysoeva M.V. Detecting changes in coupling with Granger causality method from time series with fast transient processes // Physica D. 2015. Vol. 309. P. 9.
9. Liu H., Lu J.-A., LuJ., Hill D.J. Structure identification of uncertain general complex dynamical networks with time delay // Automatica. 2009. Vol. 45. P. 1799.
10. Xu Y., Zhou W, Fang /.Topology identification of the modified complex dynamical network with non-delayed and delayed coupling // Nonlinear Dynamics. 2012. Vol. 68. P. 195.
11. Yang X.L., Wei T. Revealing network topology and dynamical parameters in delay-coupled complex network subjected to random noise // Nonlinear Dynamics. 2015. Vol. 82. P. 319
12. Chen J., Lu J., Zhou /.Topology identification of complex networks from noisy time series using ROC curve analysis // Nonlinear Dynamics. 2014. Vol. 75. P. 761.
13. Zhang Z., Zheng Z., Niu H., Mi Y., Wu S., Hu G. Solving the inverse problem of noise-driven dynamic networks // Phys. Rev. E. 2015. Vol. 91. 012814.
14. Wens V. Investigating complex networks with inverse models: Analytical aspects of spatial leakage and connectivity estimation // Phys. Rev. E. 2015. Vol. 91. 012823.
15. Hale J.K., Lunel S.M.V. Introduction to Functional Differential Equations. New York: Springer, 1993.
16. Kuang Y Delay Differential Equations with Applications in Population Dynamics. Boston: Academic Press, 1993.
17. Bocharov G.A., Rihan F.A. Numerical modelling in biosciences using delay differential equations // J. Comp. Appl. Math. 2000. Vol. 125. P. 183.
18. Mincheva M., Roussel M.R. Graph-theoretic methods for the analysis of chemical and biochemical networks. II. Oscillations in networks with delays // J. Math. Biol. 2007. Vol. 55. P. 87.
19. Heiligenthal S., Jüngling T., D'Huys O., Arroyo-Almanza D.A., Soriano M.C., Fischer I., Kanter I., Kinzel W. Strong and weak chaos in networks of semiconductor lasers with time-delayed couplings // Phys. Rev. E. 2013. Vol. 88. 012902.
20. Wu X., Sun Z., Liang F., Yu C. Online estimation of unknown delays and parameters in uncertain time delayed dynamical complex networks via adaptive observer // Nonlinear Dynamics. 2013. Vol. 73. P. 1753.
21. Сысоев И.В., Прохоров М.Д., Пономаренко В.И., Безручко Б.П. Определение параметров элементов и архитектуры связей в ансамблях связанных систем с запаздыванием по временным рядам // ЖТФ. 2014. Т. 84, вып. 10. С. 16.
22. Пономаренко В.И., Прохоров М.Д., Караваев А.С., Безручко Б.П. Определение параметров систем с запаздывающей обратной связью по хаотическим временным реализациям // ЖЭТФ. 2005. Т. 127, вып. 3. С. 515.
23. Prokhorov M.D., Ponomarenko V.I. Estimation of coupling between time-delay systems from time series // Phys. Rev. E. 2005. Vol. 72. 016210.
24. Prokhorov M.D., Ponomarenko V.I. Reconstruction of time-delay systems using small impulsive disturbances // Phys. Rev. E. 2009. Vol. 80. 066206.
25. Мандель И.Д.Кластерный анализ. М: Финансы и статистика, 1988. 176 с.
26. Mackey M.C., Glass L. Oscillations and chaos in physiological control systems // Science. 1977. Vol. 197 P. 287.
Поступила в редакцию 29.04.2016 После доработки 13.05.2016
RECONSTRUCTION OF COUPLING ARCHITECTURE AND PARAMETERS OF TIME-DELAYED OSCILLATORS IN ENSEMBLES FROM TIME SERIES
I. V Sysoev1'2, D.D. Kulminskiy2'1, V.I. Ponomarenko2'1, M.D. Prokhorov2
1 National Research Saratov State University Astrahanskaya, 83, 410012 Saratov, Russia 2Kotel'nikov Institute of Radio-Engineering and Electronics of RAS, Saratov Branch
Zelenaya, 38, 410019 Saratov, Russia E-mail: ivssci@gmail.com; kulminskydd@gmail.com; ponomarenkovi@gmail.com; mdprokhorov@yandex.ru
Purpose. To suggest a new approach to reconstruction of couping architecture and individual parameters of first-order time-delayed oscillators from experimental series of their oscillations.
Method. The method is based on minimization of target function, which characterizes a distance between points of nonlinear function of a current oscillator, which is to be reconstructed. Then estimated coupling coefficients are split into significant and insignificant. Minimization of target function is processed with least squares routine. Delay time is estimated as a trial delay corresponding to a minimum of target function over all trial delays.
Results. Efficiency of the proposed method was demonstrated in numerical experiment from time series of an ensemble of diffusively coupled nonidentical Mackey-Glass oscillators in presence of noise. Also a hardware experiment was considered in which resistively coupled generators with delay line were studied. The method demonstrated higher computational efficiency than previously suggested approaches due to use of not iterative algorithms for target function minimization and significant coefficient selection. Herewith estimates of coupling coefficients and inertance parameter are asymptotically unbiased.
Discussion. The proposed approach may be useful for reconstruction of parameters of elements and coupling architecture in systems of different nature: radioengineering, biological or others, which can be described using first-order time-delay equations.
Keywords: Time-series analyses, reconstruction of equations, ensembles of oscillators, time-delayed systems.
DOI: 10.18500/0869-6632-2016-24-3-21-37
Paper reference: Sysoev I.V., Kulminskiy D.D., Ponomarenko V.I., Prokhorov M.D. Reconstruction of coupling architecture and parameters of time-delayed oscillators in ensembles from time series // Izvestiya VUZ. Applied Nonlinear Dynamics. 2016. Vol. 24. Issue 3. P. 21-37.
References
1. Afraimovich VS., Nekorkin V.I., Osipov G.V., Shalfeev V.D. Stability, Structures, and Chaos in Nonlinear Synchronization Networks. Singapore: World Scientific, 1995.
2. Pikovsky A., Rosenblum M., Kurths /.Synchronization: A Universal Concept in Nonlinear Sciences. Cambridge University Press, United Kingdom, 2003.
3. Boccaletti S., Latora V., Moreno Y., Chavez M., Hwang D.U. // Phys. Rep. 2006. Vol. 424. P. 175.
4. Bezruchko B., Smirnov D. Extracting Knowledge From Time Series. Springer: Complexity, 2010.
5. Timme M. Revealing network connectivity from response dynamics // Phys. Rev. Lett. 2007. Vol. 98. 224101.
6. Smirnov D.A., Bezruchko B.P. Detection of couplings in ensembles of stochastic oscillators // Phys. Rev. E. 2009. Vol. 79. 046204.
7. Kammski M., Ding M., Truccolo W.A., Bressler S.L. Evaluating causal relations in neural systems: Granger causality, directed transfer function and statistical assessment of significance // Biol. Cybern. 2001. Vol. 85. P. 145.
8. Sysoev I.V., Sysoeva M.V. Detecting changes in coupling with Granger causality method from time series with fast transient processes // Physica D. 2015. Vol. 309. P. 9.
9. Liu H., Lu J.-A., LuJ., Hill D.J. Structure identification of uncertain general complex dynamical networks with time delay // Automatica. 2009. Vol. 45. P. 1799.
10. Xu Y., Zhou W., Fang J. Topology identification of the modified complex dynamical network with non-delayed and delayed coupling // Nonlinear Dynamics. 2012. Vol. 68. P. 195.
11. Yang X.L., Wei T. Revealing network topology and dynamical parameters in delay-coupled complex network subjected to random noise // Nonlinear Dynamics. 2015. Vol. 82. P. 319.
12. Chen J., Lu J., Zhou J.Topology identification of complex networks from noisy time series using ROC curve analysis // Nonlinear Dynamics. 2014. Vol. 75. P. 761.
13. Zhang Z., Zheng Z., Niu H., Mi Y., Wu S., Hu G. Solving the inverse problem of noise-driven dynamic networks // Phys. Rev. E. 2015. Vol. 91. 012814.
14. Wens V. Investigating complex networks with inverse models: Analytical aspects of spatial leakage and connectivity estimation // Phys. Rev. E. 2015. Vol. 91. 012823.
15. Hale J.K., Lunel S.M.V. Introduction to Functional Differential Equations. New York: Springer, 1993.
16. Kuang Y Delay Differential Equations with Applications in Population Dynamics. Boston: Academic Press, 1993.
17. Bocharov G.A., Rihan F.A. Numerical modelling in biosciences using delay differential equations // J. Comp. Appl. Math. 2000. Vol. 125. P. 183.
18. Mincheva M., Roussel M.R. Graph-theoretic methods for the analysis of chemical and biochemical networks. II. Oscillations in networks with delays // J. Math. Biol. 2007. Vol. 55. P. 87.
19. Heiligenthal S., Jungling T., D'Huys O., Arroyo-Almanza D.A., Soriano M.C., Fischer I., Kanter I., Kinzel W. Strong and weak chaos in networks of semiconductor lasers with time-delayed couplings // Phys. Rev. E. 2013. Vol. 88. 012902.
20. Wu X., Sun Z., Liang F., Yu C. Online estimation of unknown delays and parameters in uncertain time delayed dynamical complex networks via adaptive observer // Nonlinear Dynamics. 2013. Vol. 73. P. 1753.
21. Sysoev I.V., Prokhorov M.D., Ponomarenko V.I., Bezruchko B.P. // Tech. Phys. 2014. Vol. 59. P. 1434.
22. Ponomarenko V. I., Prokhorov M. D., Karavaev A. S., Bezruchko B. P. Recovery of parameters of delayed-feedback systems from chaotic time series // Journal of Experimental and Theoretical Physics. 2005. Vol. 100. Issue 3. P. 457.
23. Prokhorov M.D., Ponomarenko V.I. Estimation of coupling between time-delay systems from time series // Phys. Rev. E. 2005. Vol. 72. 016210.
24. Prokhorov M.D., Ponomarenko V.I. Reconstruction of time-delay systems using small impulsive disturbances // Phys. Rev. E. 2009. Vol. 80. 066206.
25. Mandel I.D. Cluster Analysis. Moscow: Finance and Statistics, 1988. 176 p. (in Russian).
26. Mackey M.C., Glass L. Oscillations and chaos in physiological control systems // Science. 1977. Vol. 197. P. 287.
Сысоев Илья Вячеславович - родился в Саратове (1983), окончил факультет нелинейных процессов СГУ (2004), защитил диссертацию на соискание учёной степени кандидата физико-математических наук (2007). Доцент базовой кафедры динамического моделирования и биомедицинской инженерии, ответственный секретарь редакционной коллегии журнала «Известия вузов. ПНД». Научные интересы - исследование сигналов биологической природы методами нелинейной динамики, исследование эффективности и модернизация подходов к анализу сигналов. Автор более 40 публикаций.
410012 Саратов, Астраханская, 83
Саратовский государственный университет имени Н.Г. Чернышевского 410019 Саратов, ул. Зеленая, 38
Саратовский филиал Института радиотехники и электроники РАН E-mail: ivssci@gmail.com
Кульминский Данил Дмитриевич -родился в 1991 году в Саратове, окончил Саратовский государственный университет в 2014 году. После окончания СГУ работает в СФ ИРЭ им. В.А. Котельникова РАН младшим научным сотрудником. Аспирант и инженер кафедры динамического моделирования и биомедицинской инженерии СГУ им. Н.Г. Чернышевского. Автор 10 научных статей в отечественных и зарубежных изданиях. Стипендиат фонда «Династия».
410019 Саратов, ул. Зеленая, 38
Саратовский филиал ИРЭ им. В.А. Котельникова РАН 410012 Саратов, ул. Астраханская, 83
Саратовский государственный университет им. Н.Г. Чернышевского E-mail: kulminskydd@gmail.com
Пономаренко Владимир Иванович - родился в Саратове (1960). Окончил Саратовский государственный университет (1982). Защитил диссертацию на соискание ученой степени кандидата физико-математических наук (1992) и доктора физико-математических наук (2008). Ведущий научный сотрудник Саратовского филиала Института радиотехники и электроники РАН. Профессор кафедры динамических систем СГУ. Область научных интересов - статистическая радиофизика, анализ временных рядов, нелинейная динамика и ее приложения. Автор более 200 научных публикаций.
410019 Саратов, ул. Зеленая, 38
Саратовский филиал ИРЭ им. В.А. Котельникова РАН 410012 Саратов, ул. Астраханская, 83
Саратовский государственный университет им. Н.Г. Чернышевского E-mail: ponomarenkovi@gmail.com
Прохоров Михаил Дмитриевич - родился в Саратове (1968). Окончил Саратовский государственный университет (1992). Защитил диссертацию на соискание ученой степени кандидата физико-математических наук (1997) и доктора физико-математических наук (2008). Заведующий лабораторией Саратовского филиала Института радиотехники и электроники им. В.А. Котельникова РАН. Область научных интересов - нелинейная динамика и ее приложения, математическое моделирование, анализ временных рядов. Имеет более 200 научных публикаций.
410019 Саратов, ул. Зеленая, 38
Саратовский филиал ИРЭ им. В.А. Котельникова РАН E-mail: mdprokhorov@yandex.ru
