Реконструкция модельных уравнений сетей осцилляторов с запаздыванием в динамике узлов и связях между ними: Обзор Текст научной статьи по специальности «Математика»

Обзоры актуальных

проблем нелинейной динамики

УДК 530.182 https://doi.org/10.18500/0869-6632-2019-27-4-13-51

Реконструкция модельных уравнений сетей осцилляторов с запаздыванием в динамике узлов и связях между ними:

Обзор

И. В. Сысоев1'2, В. И. Пономаренко1'2, М.Д. Прохоров2

Саратовский национальный исследовательский государственный университет имени Н.Г. Чернышевского

Россия, 410012 Саратов, Астраханская, 83 2Саратовский филиал Института радиотехники и электроники им. В.А. Котельникова РАН Россия, 410019 Саратов, Зеленая, 38 E-mail: ivssci@gmail.com, ponomarenkovi@gmail.com, mdprokhorov@yandex.ru Автор для переписки Илья Вячеславович Сысоев, ivssci@gmail.com Поступила в редакцию 22.07.2019

Цель данного обзора - показать современный уровень исследований в области реконструкции по имеющимся временным рядам моделей сетей, в которых отдельные узлы описываются уравнениями с запаздыванием, либо запаздывание присутствует в функциях связи.

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

Представлены результаты реконструкции по временным рядам модельных осцилляторов с различными нелинейными функциями, различными функциями связи, при вариации числа узлов в сетях в широких пределах: от трёх до нескольких десятков. Также приведены результаты реконструкции моделей по временным рядам различных радиофизических экспериментов.

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

Ключевые слова: уравнения с запаздыванием, сети связанных осцилляторов, реконструкция по временным рядам, коллективная динамика.

Образец цитирования: Сысоев И.В., Пономаренко В.И., Прохоров М.Д. Реконструкция модельных уравнений сетей осцилляторов с запаздыванием в динамике узлов и связях между ними: Обзор//Известия вузов. ПНД. 2019. T. 27, № 4. С. 13-51. https://doi.org/10.18500/0869-6632-2019-27-4-13-51

Финансовая поддержка. Работа выполнена при поддержке Российского фонда фундаментальных исследований, грант № 19-02-00071 (реконструкция модельных уравнений сетей) и в рамках государственного задания (постановка радиофизического эксперимента).

https://doi.org/10.18500/0869-6632-2019-27-4-13-51

Reconstruction of model equations of networks of oscillators with delay in node dynamics and couplings between them:

Review

I. V. Sysoev1,2, V. I. Ponomarenko1,2, M. D. Prokhorov2

1 Saratov State University 83, Astrakhanskaya, Saratov 410012, Russia 2 Kotelnikov Institute of Radioengineering and Electronics of RAS, Saratov Branch

38, Zelenaya, Saratov 410019, Russia E-mail: ivssci@gmail.com, ponomarenkovi@gmail.com, mdprokhorov@yandex.ru Correspondence should be addressed to Ilya V. Sysoev, ivssci@gmail.com Received 22.07.2019

The aim of this review is to show the modern level of research in the area of reconstruction of network models from measured time series, for which individual nodes are described by time-delayed equations or there is a delay in coupling.

Methods are described for reconstruction of coupling coefficients and functions, nonlinear functions of individual nodes and for detection of superfluous couplings. The techniques for delay time detection are considered separately due to their choice is crucial for success of entire reconstruction procedure.

There presented the results of reconstruction from times series of model oscillators with different nonlinear functions, coupling fucntions, with number of nodes in a netwoks ranging widely (from 3 to tens of nodes). In addition, the results of reconstruction of models from different radiophysical experiments are presented.

The advantages and shortcomings of proposed approaches are discussed in comparison with other known from literature methods of coupling estimation. The effects of time series length, the amount of a priori information, measurement noise, calculation errors on method efficiency are considered.

Key words: time-delayed equations, networks of coupled oscillators, reconstruction from time series, collective dynamics.

Reference: Sysoev I.V., Ponomarenko V.I., Prokhorov M.D. Reconstruction of model equations of networks of oscillators with delay in node dynamics and couplings between them: Review. Izvestiya VUZ. Applied Nonlinear Dynamics, 2019, vol. 27, no. 4, pp. 13-51. https://doi.org/10.18500/0869-6632-2019-27-4-13-51

Acknowledgements. This work was supported by Russian Foundation for Basic Research, grant No 19-02-00071 (reconstruction of model equations of networks) and within the framework of the state task (radiophysical experiment).

Введение

Ансамбли связанных дифференциальных уравнений с запаздыванием широко используются для моделирования и описания процессов в сетях различных физических [1-3], химических [4] и биологических [5-7] колебательных систем. Элементы этих сетей, будучи описаны достаточно простыми уравнениями, способны и сами по себе демонстрировать огромное богатство режимов [8,9], в том числе большое число сосуществующих аттракторов [10], что даёт надежду описывать с их использованием объекты, демонстрирующие самые различные типы колебаний. Введение связей дополнительно обогащает систему. В то же время, объём данных, которые нужно иметь об объектах, описываемых уравнениями с запаздыванием, для оценки их параметров или восстановления уравнений по измеренным временным рядам, не сильно отличается от того, который необходим для реконструкции моделей в виде обыкновенных дифференциальных уравнений.

Процесс реконструкции ансамблей систем с запаздыванием может быть условно разделён на три задачи, которые могут решаться как по отдельности, так и одновременно:

• реконструкция времени или времён запаздывания;

• реконструкция уравнений для отдельных элементов сети;

• реконструкция архитектуры связей.

Все три эти задачи нетривиальны ввиду того, что даже просто одиночная система с запаздыванием может демонстрировать бесконечномерную хаотическую динамику, и прямая реконструкция таких систем с использованием метода восстановления вектора состояния методом задержек [11] часто терпит неудачу. Многие из более продвинутых методов основаны на проекции бесконечномерного фазового пространства систем с запаздыванием на низкоразмерные подпространства [12-16]. Различные критерии используются для оценки качества реконструкции: минимальная ошибка прогноза построенной модели [12-14], минимальное значение информационной энтропии [15], или различные меры сложности прогнозируемых временных рядов [16]. В литературе также предложены некоторые другие методы и подходы к оценке параметров систем с задержкой по временным рядам: регрессионный анализ [17,18], статистический анализ временных интервалов между экстремумами временного ряда [19], анализ ближайших соседей [20], информационно-теоретические подходы [21,22], метод множественной стрельбы [23], алгоритм оптимизации искателя [24] и адаптивная синхронизация [25,26]. Отдельная группа методов восстановления систем с запаздыванием основана на анализе реакции объекта на внешние возмущения [27-30]. Однако большинство этих подходов могут быть использованы только для восстановления уравнений автономных систем с запаздыванием.

Проблема реконструкции становится ещё сложнее при наличии взаимодействия между отдельными системами с запаздыванием и требует разработки альтернативных методов [31]. Архитектура и сила взаимодействий между элементами сети определяют возможность их синхронного поведения [32,33]. Подход к одновременной оценке связанности в сети и параметров отдельных узлов, включая время запаздывания для систем с запаздыванием, был не так давно предложен в работе [34]. Однако предложенное решение имеет значительное число ограничений: возможность искажения оператора эволюции или уже измеренных рядов шумом игнорировалась, нелинейные функции элементов сети предполагались обратимыми, а стартовые догадки для времён запаздывания выбирались очень близкими к истинным значениям.

Ещё один подход был недавно предложен в работах [35-38], в основном, для систем без запаздывания, но может быть использован и для реконструкции осцилляторов с запаздыванием, если сама величина запаздывания известна. Предполагается, что все переменные всех осцилляторов измеряются, структура уравнений определена, но нелинейные функции неизвестны. Дополнительно предполагается, что все эти функции могут быть представлены в виде небольшого числа широко известных элементарных функций (см.,например, [39]).

Далее нами описывается подход, представляющий собою семейство близких по сути методов реконструкции ансамблей осцилляторов с запаздыванием в собственной динамике отдельных узлов сети и с запаздыванием в связях между ними. Узлы описываются дифференциальными уравнениями первого порядка. Основная идея предлагаемого подхода состоит в том, чтобы в качестве целевой функции использовать длину описания одной из нелинейных функций, входящих в уравнения для динамики отдельного узла. Таким образом удаётся существенно увеличить универсальность подхода и его грубость, поскольку угадывать аппроксимацию для соответствующей функции не нужно, а число оцениваемых коэффициентов также сокращается по сравнению с методами типа [39]. Представленные методы дают возможность восстановить по временным рядам нелинейные функции всех подсистем и все связи (и отличить реально присутствующие связи от ложных). Работоспособность методов демонстрируется на численных примерах и в радиофизическом эксперименте. Методы реконструкции уравнений с запаздыванием были предложены в работах [40,41], методы реконструкции уравнений с запаздыванием в связях - в работах [42-45]. Для обобщённых осцилляторов ван дер Поля без запаздывания аналогичный подход предложен в [46].

1. Уравнения ансамбля систем первого порядка с запаздыванием

Рассматриваются сети колебательных элементов, каждый из которых описывается уравнением следующего вида:

Б

£гХг(г) = -Хг (¿) + /¿(Х*^ - Т*)) + ^ к^ (ж^ (¿) - Х*(£)) , (1)

3 = 1,3 =

где г = 1,... , В - номера узлов в сети; параметры ег характеризуют инерционные свойства элементов; т* - время запаздывания; /г - произвольная гладкая нелинейная функция; кг,3 - коэффициенты связи, характеризующие силу воздействия ] -го элемента сети на г-й.

Далее предполагаем, что в нашем распоряжении имеются временные ряды {Хг(1п)}1^=1 всех В осцилляторов ансамбля, имеющие длину N, измеренные эквидистантно с шагом выборки Д£, и будем для простоты обозначать Хг(п) = Хг(£п). Также введём дискретное время запаздывания 9г = |_тг/Д£_|. Временные ряды производных {ХСг(Ьп)}1^=1 и, при необходимости, вторых производных {Хг(^п)}^=1 рассчитаем численно с использованием фильтра Савитцки-Голэя для сглаживания [47].

Следует отметить, что время запаздывания является параметром, точность восстановления которого имеет наибольшее влияние на качество реконструкции системы. Даже небольшая погрешность его определения приводит, как правило, к неправильному восстановлению архитектуры связей в сети и большим ошибкам оценки остальных параметров. Поэтому далее предлагается подход к восстановлению параметров элементов и архитектуры связей в ансамбле систем с запаздыванием (1), который состоит из двух этапов. Сначала восстанавливаем по временным рядам время запаздывания т* каждого элемента, а затем, зная т*, восстанавливаем остальные параметры (е* и кг,з) и нелинейные функции /¿. Это позволяет существенно упростить решение задачи и за счёт этого добиться более высокой точности оценки параметров.

1.1. Реконструкция времени запаздывания элементов. Ранее в [19] было установлено, что во временных реализациях изолированных (к*,3 = 0) систем с запаздыванием вида ( ) практически отсутствуют экстремумы, удаленные друг от друга на время запаздывания. Если такая система совершает хаотические колебания, то экстремумы в её временном ряде расположены нерегулярно и расстояние между ними принимает различные значения. На основе этого свойства был предложен метод определения времени задержки т*, использующий статистический анализ временных интервалов между экстремумами хаотического временного ряда системы с запаздыванием. Определив для различных значений т число V* ситуаций, при которых точки временного ряда, разделенные интервалом времени т, одновременно являются экстремальными, и построив зависимость Vi(тi), можно найти время задержки т* как значение, при котором наблюдается абсолютный минимум этой зависимости [19].

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

егхг(1) = -Хг(1) + /Х'и!^ Х*(* - т*) + Е кг,з (Х,- (*) - Хг(1)). (2)

* * 3 = 1,3=г

При наличии инерционности (еi > 0), что соответствует реальным ситуациям, экстремумы во временной реализации Жг(£) близки к квадратичным, а следовательно, в экстремальных точках х= 0, Xi(t) = 0. Тогда, если при х¿(¿) = 0 в типичном случае х^) = 0, то при еi = 0 должно выполняться условие

^"тТ))) XX i(t - *) + Е ^ XX ^ (t) = 0. (3)

i i 3 = 1,3 =

Для выполнения условия (3) необходимо, чтобы либо х^ — Тi) = 0, либо выполнялось условие

Б

Е ^ хХ 3 (t) = 0. (4)

3 = 1,3 =

Условие ( ) никогда не выполняется в случае отсутствия связи (кг^ = 0) и в случае сильной связи, обеспечивающей синхронизацию элементов, так как при этом х3 (^ = ха хi (t) = 0 при выводе условия (3). Следовательно, в этих пограничных случаях первое слагаемое в (3) отлично от нуля, а значит, производные хи х^ — тi) одновременно в нуль не обращаются, то есть на удалении т от квадратичного экстремума во временном ряде хне должно быть другого экстремума. В промежуточных ситуациях слабой и умеренной связи существует вероятность обнаружить пару экстремумов на удалении т друг от друга. Однако, как показывают численные эксперименты, в общем случае эта вероятность меньше, чем вероятность встретить пару экстремумов на удалении т = т^ В результате число разделенных интервалом времени т экстремумов будет меньше числа экстремумов, разделенных другими значениями времени т, и график у^т) будет иметь минимум при т = т . Таким образом, можно ожидать, что качественные особенности зависимости уДт) сохраняются для систем (1) в широком диапазоне значений коэффициентов связи. Отметим, что такой метод определения времени запаздывания обладает высоким быстродействием, поскольку использует только операции сравнения и сложения, не требуя вычисления каких-либо мер сложности движения или ошибки аппроксимации данных.

Однако, если система с запаздыванием совершает периодические колебания, такой подход оказывается неработоспособным, так как экстремумы во временном ряде расположены регулярным образом. Для изолированных систем с запаздыванием, находящихся в режиме периодических автоколебаний, не так давно был предложен метод восстановления времени задержки, основанный на возмущении динамики системы внешним воздействием и анализе отклика [30]. Если на переменную х^) изолированной системы с запаздыванием подействовать внешним сигналом у^), представляющим собою прямоугольные импульсы, и построить взаимную корреляционную функцию

с = (ЩЩх^ + $)рг i '

то С^в) будет иметь четко выраженный максимум при в = т^ Недостатком такого подхода по сравнению с рассмотренным выше методом является необходимость активного воздействия на систему. С другой стороны, метод позволяет использовать импульсы малой амплитуды, что позволяет свести воздействие на систему к минимуму.

Исследуем возможность применения такого метода для определения времен задержки в ансамбле связанных систем с запаздыванием. Рассмотрим такой способ возбуждения элементов ансамбля внешними сигналами yiпри котором их модельные уравнения имеют следующий вид:

Б

еiXi(t) = — х^) + ^ (х^ — т) + у^ — *))+ £ к^ (х^(t) — х^)), (6)

3 = 1,3 =

где сигнал у* (¿) представляет собою прямоугольные импульсы с амплитудою А*, периодом Т* и длительностью Тг'. Для восстановления времени задержки т* только одного г-го элемента ансамбля достаточно подействовать внешним сигналом только на этот г-й элемент. Как уже было отмечено выше, наличие взаимодействия между системами с запаздыванием приводит к возмущению фазовых траекторий. Эти возмущения снижают чувствительность взаимной корреляционной функции (5) как характеристики определения времени задержки. В результате для восстановления т* в общем случае необходимо увеличить амплитуду внешнего воздействия А* по сравнению со случаем несвязанных систем с запаздыванием. Метод можно использовать при любых значениях коэффициентов связи кг,3. Кроме того, этот метод может быть применен к системам (1), совершающим как периодические, так и хаотические колебания. Ещё одним достоинством метода является то, что он остается эффективным при высоких уровнях шума, в несколько раз превышающих допустимый уровень шума для метода определения т*, основанного на статистическом анализе экстремумов временного ряда.

2. Реконструкция параметров инерционности, нелинейных функций и архитектуры связей

Определив т*, можно восстановить параметр е*, нелинейную функцию /г и коэффициенты связи кг,3 г-й системы с запаздыванием ( ), имея временные ряды колебаний всех элементов ансамбля. Для этого предлагается следующий подход. Уравнение (1) перепишем в виде

Б

/г(Хг(£ - Т*)) = е*Хг(¿) + Х*(£) - ^ кг,3 (Х3(¿) - Х*(£)) . (7)

3=1,3=г

Если построить зависимость части уравнения (7) от х*(£ - т*), то она воспроизведет функцию /¿. Поскольку заранее величины е* и кг,3 не известны, будем искать их, минимизируя некоторую величину характеризующую расстояние между точками на плоскости (хг, /¿(х)), отсортированными по возрастанию х*. Пусть рп - это номер в исходном ряде момента измерения, соответствующего в отсортированном по возрастанию х* ряде значению, непосредственно следующему за Хп. Тогда можно записать в виде

N — 0^—1

¿г(е*, кг,3) = Е ((Хг(Рп) - Х*(п))2 + (/(Хг(рп)) - /(х*(п)))2) =

п=1 N — 1

Е ((Хг(Рп - бг) - Хг(п - б*))2 +

N—е4—1

// / ч / ч\2 . /

-ч-^--/ -ч // .„ чХг

п=1 N — 1

Хг(/7п

п=0;+1

+ (ег(Хг(Рп) - Х*(п)) + Х*(рп) - Х*(п) - (8)

Б

кг,3 Е (Х3(Рп) - Хг(рп) - Х3(п) + Х*(п)))^ . 3=1,3=г

При ошибочном выборе значений е* и кг3 точки на плоскости (х*, /¿) не ложатся на одномерную кривую, образуя «облако», величина £г(ег,кг)3) будет многократно больше, чем при истинных е* и кг>3-.

Для ег и кг,3 задаются стартовые догадки, которые затем уточняются симплекс-методом [48], минимизируя величину ( ), минимум которой обозначим . При В ^ 4 и отсутствии шума все параметры восстанавливаются с высокою точностью. Однако уже при В > 4 типичной является ситуация, когда метод не позволяет выявить отсутствующие связи между элементами

ансамбля. Такие связи диагностируются как слабые из-за наличия опосредованных связей через другие элементы. Избавиться от незначимых связей позволяет метод последовательного пробного исключения коэффициентов fc¿j из модели ( ). Для этого выдвигается гипотеза отсутствия связи j ^ i между двумя элементами, соответствующее слагаемое с коэффициентом связи исключается из модели, и остальные параметры модели восстанавливаются, как и ранее, по минимуму целевой функции L( ). Затем для зафиксированного i процедура повторяется, исключая другой kjj, и так далее для всех j = i. При этом на каждом шаге предполагается отсутствие воздействия на i-й элемент лишь со стороны одного из остальных j-х элементов ансамбля. Наконец, определяется, при исключении какого из fc¿j получается minj L¿,d-i, и оценивается статистическая значимость величины Л = L¿;d_i/L¿;d из следующих соображений.

При больших N разности, входящие в (8) в квадрате, распределены по закону, близкому к нормальному, причём N/2 из них можно считать независимыми, поскольку они не имеют общих координат. Кроме того, LÍ;d зависит от D параметров модели ( ), что уменьшает общее число независимых величин в (8) до N/2 — D. Тогда, учитывая, что сумма квадратов K независимых нормально распределенных величин распределена по закону х2 с K степенями свободы [49], получим, что распределение величин L, полученных при разных значениях параметров и/или наличии шума, подчиняется закону х2 с N/2 — D степенями свободы, а величин L¿,d_i - закону X2 с N/2 — D + 1 степенями свободы.

Известно, что величина X , являющаяся отношением двух независимых случайных величин, распределённых по закону х2, имеет распределение Фишера-Снедекора с функцией распределения

(X) = 2, (9)

где Bd - неполная бета-функция, а d = vX/(vX + w) [50]. Следовательно, величина Л имеет функцию распределения ( ) с параметрами v = N/2 — D + 1 и w = N/2 — D. Обозначим через Л1_р такое значение Л, при котором (L, 1 — p) = 1 — p, где p - уровень статистической значимости. Тогда, если Л > Л1_р, то на уровне значимости p можно сделать вывод о наличии связи между элементами, а значит, все = 0. В противном случае делается вывод об отсутствии связи j ^ i между соответствующими элементами и далее проверяется значимость остальных связей, последовательно исключая из оставшихся связей i-го элемента по одной. Процедура повторяется, пока все связи не окажутся значимыми. Такой подход позволяет восстановить архитектуру связей, параметры всех элементов и их нелинейные функции.

Если известно, что число связей между элементами ансамбля мало, то для восстановления архитектуры и величины связей предпочтительнее использовать метод последовательного пробного добавления коэффициентов fc¿j в модель ( ). Сначала найдем минимум L¿ i функции ( ), предположив, что в уравнении ( ) отсутствуют все (связей нет). Затем будем вводить в ( ) по одному fcij, находя минимум Li 2 функции ( ). Перебрав все j = i, найдём Lmin;i)2 = minj(Lг,2). Если Л > Л1-р, где Л = L¿д/Lmin,i)2, а построена при v = N/2 — 1, w = N/2 — 2, то введённая связь отлична от нуля на уровне значимости не ниже p. Процедура повторяется, пока очередная добавленная в модель связь не окажется незначимою.

3. Апробация метода в численном эксперименте

3.1. Восстановление уравнений цепочки однонаправленно связанных уравнений Икеды. В качестве первого примера восстановим параметры элементов и архитектуру связей в цепочке однонаправленно связанных элементов (рис. 1, a), описываемых уравнением Икеды [51]

¿j(í) = —x¿(í) + Цг Sin (x¿(í — Ti) — Жо,г) + fc^-l(x¿_i — £¿(í)). (10)

Число элементов в цепочке D = 10, и использовано кольцевое граничное условие xi = xd+i. Уравнение Икеды описывает сдвиг x фазы электрического поля в нелинейной поглощающей среде кольцевого резонатора. Параметр ^ характеризует интенсивность лазерного излучения, подаваемого на вход резонатора, т — время запаздывания, x0;i — постоянный фазовый сдвиг. Уравнение ( ) является частным случаем уравнения ( ) с = 1. Пусть замкнутая в кольцо цепочка состоит из неидентичных элементов, параметры которых принимают случайные значения в следующих интервалах: т € [2, 5], ^ € [15, 25], x0)i € [0, 2п], ki,i-1 € [0.1, 0.5], и на каждый элемент действует независимый нормальный белый шум ^¿(i) с нулевым средним и дисперсией о2 = 0.01. При этом все элементы колеблются хаотически. Проиллюстрируем результаты применения предложенной методики восстановления параметров на примере одного из элементов цепочки с параметрами т7 = 2.15, ц7 = 21.67, x0,7 = 3.88, fc7,6 = 0.284. На рис. , b приведён фрагмент временного ряда колебаний в этом элементе цепочки. Подсчитав число V7 одновременных обращений в нуль x7(t) и x7(t — Т7) для различных значений Т7, перебираемых с шагом, равным выборке интегрирования At = 0.01, построим зависимость V7(t), введя нормировку V7 на общее число экстремумов в ряде (рис. 1, с). Для оценки производной по временному ряду была использована локальная аппроксимация квадратичным полиномом. Абсолютный минимум V7(t) наблюдается при истинном времени запаздывания т = Т7 = 2.15. Для построения графика V7(t) был использован временной ряд длиною в N = 40000 отсчётов, содержащий около 1600 экстремумов.

На рис. 1, d приведена нелинейная функция /7, реконструированная с помощью метода последовательного пробного добавления коэффициентов связи в модель при p = 0.05. Функция /7

5 0 -5 -10

a

b

"Г

10

г

30

12 3 t

f 20

0

-20

d

-10 0 10 X,

12 3 4 5 6789 i

Рис. 1. a — замкнутая в кольцо цепочка из 10 однонаправленно связанных систем с запаздыванием; b — хаотический временной ряд переменной x7 (t) цепочки (10) в присутствии нормального шума с дисперсией а2 = 0.01; c — зависимость v7(t), v7 min = v7(2.15); d — функция f7, восстановленная на плоскости (x7, f7(x)); e — диаграмма результата восстановления архитектуры связей: чёрным цветом показаны правильно восстановленные связи, белым - правильно диагностированные отсутствующие связи

Fig. 1. a - ten unidirectionally coupled in ring time-delayed systems; b - chaotic time series of variable x7(t) from ensemble (1 ) in presence of normal noise with varience а2 = 0.01; c - dependence v7(t), v7min = v7(2.15); d - function f7 reconstructed on the plane (x7,f7(x)); e - diagram of of reconstructed coupling architecture, where correctly revealed couplings are shown in black and couplings rightly detected as absent are shown in white

0

e

c

построена при восстановленных значениях параметров т7 = 2.15, к7в = 0.276. Она достаточно хорошо совпадает с истинной функцией уравнения Икеды. Аппроксимация восстановленной нелинейной функции гармоническою функцией позволила получить следующие оценки параметров: ц'7 = 21.80 и 7 = 3.97, которые находятся как амплитуда и начальная фаза гармонической функции, соответственно. Аналогичным образом проводилось восстановление параметров и коэффициентов связи остальных элементов цепочки.

Результат реконструкции архитектуры связей в цепочке, полученный с помощью метода последовательного пробного добавления коэффициентов связи в модель, приведен на рис. 1, e. Клетка с координатами (г, ]) показывает влияние ] ^ г, кроме клеток на диагонали, не несущих никакой информации. На уровне значимости р = 0.05 найдены все десять реально имеющихся связей, ложных связей не обнаружено. Отметим, что для построения рис. , d,e использованы временные ряды длиною в 10000 отсчётов, то есть существенно короче, чем для поиска величины запаздывания.

Далее параметры системы (10) задавались таким образом, чтобы все элементы цепочки совершали периодические колебания в отсутствии связи, при этом связи приводили к слабому возмущению этого режима. Для этого параметры элементов случайным образом генерировались из равномерного распределения на следующих интервалах: т € [0.5,1.2], ц € [4.5, 5.5], ж0,г € [0, 2п], к ,г_1 € [0.1, 0.5]. Кроме того, на каждый элемент системы действовал независимый нормальный шум с нулевым средним и дисперсией о2 = 0.001. На рис. , a приведён фрагмент временного ряда колебаний в седьмом элементе цепочки, параметры которого: т7 = 0.94, ц7 = 5.17, ж0,7 = 3.32, к7 в = 0.284. Поскольку метод оценки времени запаздывания, основанный на статистическом анализе экстремумов временного ряда, неприменим к периодическим временным рядам, для определения времени запаздывания Т7 подействуем на переменную Ж7(£) слабым внешним сигналом ут(£), представляющим собой прямоугольные импульсы, таким образом, чтобы динамика возбуждаемого элемента описывалась уравнением (6). На рис. 2, Ь построена взаимная корреляционная функция (5) для случая, когда возмущающий импульсный сигнал имеет амплитуду А 7 = 0.1, период Т7 = 3 и длительность Т7 = Т7/2. При шаге изменения в, равном 0.01, С7(в) имеет высокий максимум при в = 0.94, то есть время запаздывания восстанавливается точно.

С помощью метода последовательного пробного добавления коэффициентов связи в модель на уровне значимости р = 0.05 были получены следующие оценки параметров седьмого

х7 2 -

1 -

-1

-2 -

а

t

C7 0.03

0.02

0.01

0.00

Ъ

0.0

0.4

0.8

1.2

Рис. 2. a — временной ряд переменной x7(t) цепочки ( ) в периодическом режиме в присутствии нормального динамического шума с а2 = 0.001; b — взаимная корреляционная функция C7(s), C7,max(s) = C7(0.94)

Fig. 2. a — time series of variable x7(t) of the chain (1С) in periodic regime in presence of normal white noise with а2 = 0.001; b — cross-correlation function C7(s), C7,max(s) = CV(0.94)

0

5

элемента: ц7 = 5.21, 7 = 3.70, k7 6 = 0.303. Для остальных элементов погрешности восстановления параметров имели те же значения по порядку величины. Результат реконструкции архитектуры связей в цепочке совпадает с результатом, полученным в случае хаотической динамики (см. рис. 1, e): все связи были найдены правильно.

3.2. Восстановление ансамбля произвольно связанных уравнений Икеды. Рассмотрим ситуацию, когда элементы ансамбля, описываемые уравнениями Икеды вида (10), связаны не в кольцо, а достаточно произвольным образом так, что они описываются уравнениями (11).

D

¿i(í) = -Xi(t) + sin - Ti) - x0 ,i) + E ki (xj(t) - Xi(t)). (11)

j=1 , j=i

На рис. 3, a приведена архитектура случайно выбранных связей в ансамбле из 10 элементов. Из 90 возможных связей между элементами ансамбля реально присутствуют 40, причём некоторые элементы связаны только однонаправленно, а другие - взаимно. Параметры неидентичных элементов заданы таким же образом, как в первом рассмотренном примере для цепочки (10).

Рис. 3. Архитектура связей в ансамбле из 10 произвольно связанных элементов, описываемых уравнениями (11) (a); зависимость v7(t), v7,min(T) = v7(2.15) (b); функция f7, восстановленная на плоскости (x7,f7) (c); диаграмма результата реконструкции архитектуры связей (d)

Fig. 3. Coupling architecture in the ensemble of 10 arbitrary coupled elements described with equations (11) (a); dependence v7(t), v7,min(t) = v7(2.15) (b); function f7 (c); diagram of coupling reconstruction results (d), where colors follow the pattern from Fig. 1

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

Как и для ранее рассмотренных примеров, далее проиллюстрированы результаты восстановления седьмого элемента с параметрами т7 = 2.15, ц7 = 21.67, х0,7 = 3.88, к7;4 = 0.445, к7,б = 0.172, ^7,9 = 0.311, ^7,10 = 0.435, ^ = 0, при э = 1, 2, 3, 5, 7, 8. На рис. 3, Ь приведена зависимость У7(т) при шаге изменения т, равном 0.01. Минимум У7(т) наблюдается при истинном времени запаздывания т = т7 = 2.15.

На рис. 3, c приведена нелинейная функция /7, реконструированная с помощью метода последовательного пробного добавления коэффициентов связи в модель на уровне значимости р = 0.05. Функция /7 построена при восстановленных значениях параметров т7 = 2.15, А7 4 = 0.517, А7,6 = 0.188, А7,9 = 0.355, А7,10 = 0.490, ^ = 0, для всех э = 1,2,3, 5, 7,8. Из-за погрешности определения коэффициентов связи и более высокого уровня шума качество восстановления /7 несколько хуже, чем на рис. 1, d. Аппроксимация реконструированной функции /7 гармоническою функцией позволяет оценить и параметры ^ в исходных уравнениях (11), в частности, для рассматриваемого элемента р/7 = 22.00.

Параметры связи и индивидуальных узлов были аналогично оценены для всех 10 элементов ансамбля. Результат реконструкции архитектуры связей приведён на рис. 3, d. На уровне значимости р = 0.05 все 40 связей определяются правильно при использовании как метода добавления связей, так и метода последовательного пробного исключения коэффициентов связи из модели.

Далее рассмотрен случай, когда на элементы ансамбля (11) с теми же значениями параметров и той же конфигурацией и величиной связей, что на рис. 3, a, действует сильный независимый нормальный шум с нулевым средним и дисперсией о2 = 0.36. При таком большом шуме метод восстановления времени запаздывания, основанный на статистическом анализе экстремумов временного ряда, оказывается неэффективным и приходится прибегать к активному эксперименту. Для определения времени запаздывания г-го элемента подействуем на переменную Жг(£) слабым прямоугольным импульсным сигналом Уг(^) таким образом, что динамика возбуждаемого элемента описывается уравнением (6). На рис. 4, a построена взаимная корреляционная

Рис. 4. Случай высокого уровня шума: a — взаимная корреляционная функция C7(s), C7,max (s) = C7(2.15); b - диаграмма результата восстановления архитектуры связей: чёрным цветом показаны правильно восстановленные связи, белым — правильно диагностированные отсутствующие связи, серым цветом — пропущенные связи

Fig. 4. Reconstruction from series significantly contaminated with noise: a - cross-correlation function C7(s), C7,max(s) = = C7(2.15); b - diagram of coupling reconstruction: rightly reconstructed couplings are colored in black, missed couplings are colored in gray, absent couplings are colored in white

Рис. 5. Случай высокого уровня шума: a — взаимная корреляционная функция C7(s), C7,max(s) = C7(2.15); b — диаграмма результата восстановления архитектуры связей: белым цветом показаны правильно восстановленные связи, крестиками — правильно детектированные отсутствующие

Fig. 5. Reconstruction from series significantly contaminated with noise: a - cross-correlation function C7(s), C7,max(s) = = C7(2.15); b - diagram of coupling reconstruction: rightly reconstructed couplings are colored in white, correctly determined absent couplings are marked with crosses

функция ( ) для случая, когда возмущающий импульсный сигнал с амплитудой А 7 = 0.15, периодом Т7 = 5 и длительностью Т7 = Т7/2 действует на переменную Ж7(£). При шаге изменения в, равном 0.01, С7(в) имеет выраженный максимум при в = Т7 = 2.15. Отметим, что этот результат не уникален и, несмотря на высокий уровень шума, время запаздывания удается точно восстановить для всех элементов.

Метод добавления связей, так же как и метод исключения связей, даёт на уровне значимости р = 0.05 следующую оценку параметров седьмого элемента: ^7 4 = 0.396, ^7 б = 0.142, ^7,9 = 0.262, ^7 10 = 0.392, ^7 ^ = 0 при ] = 1, 2, 3, 5, 7, 8, ц'7 = 21.52, ж0 , 7 = 3.97. Истинные значения параметров такие же, как приведённые выше для рис. 3. Результат реконструкции архитектуры связей в ансамбле приведён на рис. , Ь. На уровне значимости р = 0.05 из 40 существующих связей удалось найти 35. Таким образом, из-за высокого уровня шума пропущенными оказались 5 связей. Отметим, что, увеличивая р, можно уменьшить количество пропущенных связей, однако при этом повышается вероятность обнаружения ложных.

Все ранее рассмотренные варианты характеризовались относительно редкими связями в ансамбле: 10 или 40 из 90. Поэтому далее рассмотрен случай, когда из 90 возможных связей между элементами ансамбля (11) имеется 80 связей. Параметры элементов заданы случайным образом в тех же интервалах, что и в рассмотренном выше примере с 40 связями. На каждый элемент действует независимый нормальный шум с нулевым средним и дисперсией о2 = 0.01. Положение абсолютного минимума зависимостей ^ (т) позволяет точно восстановить истинное время запаздывания для каждого элемента ансамбля. Для примера на рис. 5, a приведена зависимость У7(т), построенная при шаге изменения т, равном 0.01. Она имеет минимум при т = Т7 = 2.40. Результат реконструкции архитектуры связей, полученный с помощью метода последовательного пробного исключения коэффициентов связи из модели, приведен на рис. 5, Ь. Для удобства восприятия цвета были изменены по сравнению с другими диаграммами (иначе почти всё поле было бы закрашено чёрным). На уровне значимости р = 0.05 были правильно обнаружены все 80 связей, ложно положительных результатов метод не дал.

3.3. Восстановление ансамбля связанных уравнений Маккея-Гласса. Чтобы охарактеризовать общность полученных результатов, были рассмотрены ансамбли систем с запаздыва-

нием с другою функцией нелинейности - связанные системы Маккея-Гласса [52], для которых / представляется в виде ( ), а параметр инерционности = 1/6г-

fi(x) =

bi(1+ )

(12)

Уравнение Маккея-Гласса, описывающее процесс выработки организмом красных кровяных клеток, является эталонным уравнением с запаздыванием, широко используемым при численных исследованиях систем с задержкою.

Будем считать, что ансамбль состоит из неидентичных элементов, параметры которых случайно выбираются из следующих равномерных распределений: т € [300, 400], £ € [8,12], a € [0.2, 0.25], kij € [0.01, 0.05], и на каждый элемент дополнительно действует независимый нормальный шум ^¿(i) с нулевым средним и дисперсией о2 = 10-4. При этом все элементы при отсутствии шума и связи колеблются хаотически.

На рис. 6, a приведён фрагмент временного ряда колебаний седьмого элемента при т7=335, £7 = 10.2, Ä7,1 = 0.011, k7,2 = 0.046, Й7,3 = 0.043, ^7,4 = 0.016, k7j = 0 для всех 5 < j < 10. На рис. 6, b приведена зависимость V7(t), построенная по временному ряду переменной X7(t) длиною 40 000 отсчётов, содержащему около 2600 экстремумов. При шаге изменения т, равном 1, минимум V7(т) наблюдается при истинном времени запаздывания т = Т7 = 335.

На рис. 6, c серым цветом показана функция /7, полученная при реконструкции уравнения (1) в предположении, что все элементы ансамбля не связаны между собою. Эта функция построена при восстановленных значениях параметров т7 = 335, = 8.4. Чёрным цветом на

a

1000

2000

3000

4000

0.10

0.05

0.00 b

О 200 400 X

f7 1.2

0.0 0.6 1.2

c

J 9 8 7 6 5 4 3 2 1

d

1 2 3 4 5 678 9 i

Рис. 6. a — временной ряд переменной x7(t) ансамбля ( ) с функцией ( ) при добавлении в динамику нормального шума с a2 = 10-4; b — зависимость v7(t), v7>min(x) = v7(335); c — функция /7, восстановленная на плоскости (x7, /7) в предположении отсутствия связей (серый цвет) и при учёте связей (чёрный цвет); d — диаграмма результата восстановления архитектуры связей: чёрным цветом показаны правильно восстановленные связи, белым — правильно диагностированные отсутствующие связи, точками — ложные связи

Fig. 6. a — time series of variable x7(t) of the ensemble described by equations ( ) with nonlinear functions of type (12) and with white noise with varience a2 = 10-4 introduced into dynamics; b — dependence v7(t), v7,min(x) = v7(335); c — function /7, reconstructed on the plane (x7,/7) assuming no couplings (gray) and after reconstruction of couplings (black); d — diagram of coupling reconstruction: rightly reconstructed couplings are colored in black, superfluous couplings are shown by dotted bars, absent couplings are colored in white

t

7

рис. 6, c показана функция /7, восстановленная с помощью метода последовательного пробного добавления коэффициентов связи в модель с контролем значимости на уровне р = 0.05. Эта функция построена при восстановленных значениях параметров т^ = 335, е7 = 10.0, Л7,1 = 0.012, ^7,2 = 0.047, ¿4,3 = 0.044, ¿7,4 = 0.017, ¿7. = 0 для всех 5 < ; < 10. Учёт архитектуры связей существенно улучшает качество восстановления нелинейной функции и точность оценки параметров модели. Погрешности их восстановления вызваны преимущественно присутствием шума. Как и в остальных рассмотренных примерах, для восстановления параметров использованы временные ряды длиною 10 000 отсчётов.

Аналогичным образом проводится восстановление параметров и коэффициентов связи остальных элементов ансамбля. Результат реконструкции архитектуры связей в ансамбле, полученный с помощью метода добавления связей, приведён на рис. 6, d. На уровне значимости р = 0.05 мы правильно нашли все 40 имеющихся связей. Дополнительно метод показал наличие ещё двух связей, которые являются ложными. Однако доля этих ложных связей в общем количестве связей не превышает р. Как и ранее для сязанных систем Икеды, уменьшая р, можно избавиться от ложных связей, однако при этом некоторые связи окажутся пропущенными. При восстановлении архитектуры связей по тем же временным рядам с помощью метода исключения связей было получено большее количество ложных связей при том же р.

4. Восстановление уравнений ансамбля связанных экспериментальных радиотехнических генераторов с запаздывающей обратной связью

После апробирования на различных численных примерах предложенный метод был применён к экспериментальным временным рядам трёх связанных неидентичных радиотехнических генераторов с запаздывающею обратною связью. Блок-схема экспериментальной установки приведена на рис. 7, a. Она состоит из трёх связанных между собой кольцевых генераторов, каждый из которых включает линию задержки, нелинейный элемент и низкочастотный ЯС-фильтр первого порядка. Нелинейные элементы и линии задержки были выполнены на микроконтроллерах, а фильтры - на аналоговых элементах. Аналоговые и цифровые элементы схемы сопрягались с помощью аналого-цифровых и цифро-аналоговых преобразователей. Связь генераторов осуществлялась с помощью резисторов Яс. Модельное уравнение для г-го элемента ансамбля имеет вид (13), где Рг(£) и V» (£ — т») - напряжения соответственно на входе и выходе линии задержки г-го элемента, т» — время запаздывания, Я» и С» — сопротивление и ёмкость элементов фильтра, / — передаточная характеристика нелинейного элемента.

Б

ЯгСг%) = + / (^ — т))+ £ (V.(*) — ВД) . (13)

•7 = 1,.?=»

Уравнение ( ) может быть сведено к виду ( ) с использованием замен е» = Я»С», ж» = V».

Все три нелинейных элемента имели квадратичную передаточную характеристику /». Хаотические сигналы записывались с помощью трёхканального аналого-цифрового преобразователя с частотой выборки / = 10 кГц. На рис. 7, Ь приведён фрагмент временной реализации сигнала У1(£) в первом генераторе, имеющем параметры т1 = 13.6 мс, е1 = 2.88 мс, ¿1,2 = Я1/Яс = 0.1, ¿1,3 = 0.

При шаге изменения т, равном времени выборки точек 0.1 мс, абсолютный минимум зависимости У1(т) наблюдается при т = 13.6 мс (рис. 7, в). На рис. , d приведена функция /1, восстановленная по экспериментальным временным рядам с помощью метода добавления коэффициентов связи в модель при р = 0.05. Эта функция построена при восстановленных значениях

АГ)П DL-1 ND-1 ПАП

R

Vi(t)

ADC П1 -? ND-2 DAC

R,

v2(t)

Rc

Co

V3(t)

ADC DL-3 ND-3 DAC

V, V 4 -3 -

i H

0

N\

a

I

b

0 50 100 150 t, ms

0.03 -0.02 -0.01 -0.00

fi, V 4 -

2 -

c

0 10 20 x, ms 0 1 2 3 4 x1, V 1

d e

2 3 i

0

Рис. 7. a — блок-схема экспериментальной установки, где DL-1, DL-2 и DL-3 - линии задержки генераторов, ND-1, ND-2 и ND-3 — их нелинейные элементы, ADC — аналого-цифровые преобразователи, DAC — цифро-аналоговые преобразователи; b — хаотическая временная реализация V1(t) первого генератора; c — зависимость vi(x), v1>min(x) = = v1 (13.6 ms), d — восстановленная функция f1; e — диаграмма результата восстановления архитектуры связей

Fig. 7. a - block-scheme of the experimental device, where DL-1, DL-2 and DL-3 are delay lines of generators, ND-1, ND-2 and ND-3 are their nonlinear elements, ADC are analog-digital converters, DAC are digital-analog converters; b - chaotic time series V1(t) of the first generator; c - dependence v1(x), v1>min(x) = v1 (13.6 ms); d — reconstructed nonlinear function f1; e — diagram of results of coupling reconstruction

параметров т1 = 13.6 мс, е1 = 2.74 мс и 2 = 0.098. Незначимый член связи с коэффициентом з не был включён в модель. Восстановленная функция достаточно хорошо совпадает с истинною передаточною характеристикой /1 нелинейного элемента первого генератора. Метод исключения связей даёт для этой экспериментальной системы точно такие же результаты. Характеристики двух остальных элементов были восстановлены аналогичным образом. Результат реконструкции архитектуры связей в ансамбле приведен на рис. 7, e. Видно, что на уровне значимости р = 0.05 все связи найдены верно.

5. Модифицированный подход к реконструкции

Предложенный подход имеет два существенных недостатка: во-первых, он опирается на отдельную процедуру оценки времени запаздывания (либо метод, основанный на статистике экстремумов [19], либо метод, основанный на расчёте автокорреляционной функции [30]), во-вторых, он использует симплекс-метод для минимизации целевой функции (8) и итеративный подход для определения значимых связей. Оба эти недостатка существенно снижают простоту применения и многократно увеличивают вычислительное время. При этом для достаточно

больших ансамблей от 10 элементов и более процедуры добавления и исключения коэффициентов становятся недостаточно эффективными из-за частичной синхронизации между узлами, а симплекс-метод часто ведёт к локальным минимумам целевой функции. Поэтому в данном разделе предлагается модификация исходного подхода, которая в значительной степени лишена приведённых недостатков и может быть эффективно использована для ансамблей большего размера. Основная идея заключается в том, чтобы переформулировать целевую функцию так, что задача сводится к линейному методу наименьших квадратов. Это многократно снижает вычислительную сложность алгоритма и увеличивает надёжность результатов. В результате, время запаздывания может быть найдено перебором, так как оно фактически может быть выбрано только дискретно и должно быть кратно шагу выборки ДЬ.

5.1. Вывод целевой функции. Будем по-прежнему рассматривать систему (1) и сохраним все предположения и обозначения, дополнительно обозначив значения всех переменных ансамбля в момент времени Ь как вектор Хг(£); тогда, следуя ранее введённым упрощениям, значения всех переменных в момент времени Ьп можно обозначить как х^(п). Используя дискретное время запаздывания 9^, перепишем уравнение (1) в следующем виде:

э

/ (Жг(п)) = од(п + бг)+ жДп + — ^ ^,3 (X(п + 0г) - «¿(п + )) , (14)

где п = 1,..., N — 0г. Для каждого осциллятора отсортируем значения «¿(п) по возрастанию, обозначив такую сортировку как преобразование ставящее в соответствие точке с номером п в исходном ряде точку с номером ^г(хг,п) в отсортированном ряде. Обратное преобразование, ставящее в соответствие точке с номером ^(хг,п) в отсортированном ряде точку с номером п в исходном, обозначим Ф-1. Тогда п = ^-1(^г(п)). Для краткости обозначений здесь и далее зависимость ф и ф-1 от х^ указывать не будем, ограничившись индексом г.

Пусть некоторая точка имеет номер п в исходном ряде и номер ф (п) в отсортированном ряде. Тогда её сосед справа в отсортированном ряде будет иметь номер ^¿(п) + 1, а в исходном ряде его номер будет рп = ^¡^(^¿(п) + 1), причём номера п и рп в общем случае не будут близки. Так как точки с номерами ^¿(п) и ^¿(п) + 1 являются соседними в отсортированном ряде, значения динамической переменной в этих точках будут близки. А значит, будут близки и значения функции / от этих переменных, поскольку все функции /¿, как уже было сказано выше, изначально предполагаются непрерывными. Обозначим разность значений функции / в этих точках как бДп)

б;(п) = / («¿(Рп)) — / (жДп)). (15)

Запишем уравнение (15), используя уравнение (14)

э

6Дп) = ( «¿(Рп + 9») — Е (жз(Рп + 9^) — «¿(Рп + 01^ + ¿(Рп + 9»)

3=1,3=»

— | «¿(п + 01) — Е ^,3 (ж3 (п + 01) — «¿(п + 01^ + ¿(п + 0г) | .

3=1,3=«

(16)

Введём новые обозначения и перепишем выражение (16) в следующем виде:

э

61 (п) = Джг(n) — Е ^,3 (Д«3 (п) — ДЖг(n)) + бД«¿(п), (17)

3=1,3=1

Джг(n) = «¿(Рп + 01) — «¿(п + 01), (18)

Д« ¿(п) = « ¿(Рп + 01) — ¿¿(п + 01). (19)

Обозначим через сумму 6? (п)

N-0г-1 N-0г-1 , Э ч 2

¿2 = Е б2(п)= Е ^«¿(п) — Е ^¿,3(^3(п) — Д«;(п)) — (—¿(п)) . (20)

п=1 п=1 3=1,3 = ®

Величину можно рассматривать как функцию от параметров 02, ^¿,3 и б, значения которых заранее не известны. При правильном выборе этих параметров будет меньше, чем при ошибочном. Это объясняется тем, что при неправильном выборе 01, к,3- и б расстояния (17) не будут малы даже для соседних точек в отсортированном ряде.

Заметим, что в отличие от ранее использованной целевой функции (8) предложенная мера в некотором смысле характеризует длину восстановленной нелинейной функции, а именно, представляет собою сумму квадратов только вертикальных компонент расстояний между точками нелинейной функции, в то время как горизонтальные компоненты ^(рп) — «¿(п)) не поддаются оптимизации по параметрам 01, к,3 и б и потому не включены в целевую функцию (20).

Поскольку ( ) и (19) зависят только от 01, б и ^¿,3 при данном фиксированном г и не зависят от иных 0т, бт или кт,3, для которых т = г, минимизацию целевой функции ( 0) возможно проводить отдельно для каждого осциллятора. При этом фактически решается задача реконструкции неавтономной системы по векторному ряду {«¿(п + 01),«¿(п), А¿(п + 0г)}^=-|ег, одна компонента которого измеряется, другая восстанавливается методом задержек, а третья -методом дифференцирования. Временные ряды членов связи вычисляются явно по известным рядам переменных {«т,(п)}^=0.+1. Таким образом, нет необходимости восстанавливать фазовое пространство большой размерности, пропорциональной числу осцилляторов в ансамбле.

Если зафиксировать 01, то задачу минимизации (20) можно рассматривать как линейную задачу о наименьших квадратах, где к,3 и (—б) суть искомые коэффициенты (их всего В штук для каждого г-го элемента), Джг(n) - аппроксимируемые величины, а матрица значений базисных функций состоит из (Д«3 (п) — Д«!^)) и Д«¿(п). В такой постановке представляет собою целевую функцию. Поиск её экстремума является стандартною задачей, которая может быть решена нерекурсивно.

5.2. Аналитические свойства оценок. С увеличением числа N точек во временном ряде Хг число членов суммы (20) будет расти пропорционально N. Вместе с тем, с ростом N будут уменьшаться расстояния между точками в отсортированном ряде и, как следствие, уменьшаться величины ^¿(п)|. В среднем это уменьшение пропорционально . То есть каждый член 62(п) суммы (20) будет убывать пропорционально 1 с ростом N. Следовательно, ^ 0 при N ^ то. Значит, при N ^ то предложенный метод является асимптотически точным, а полученные с его помощью оценки параметров являются асимптотически несмещёнными.

Предложенный подход сформулирован для систем, описываемых уравнениями с запаздыванием, содержащими только непрерывные нелинейные функции. Однако, если такие функции будут содержать на всём отрезке изменений аргумента небольшое число разрывов первого рода и сумма квадратов скачков функции будет много меньше общей суммы квадратов вертикальных

расстояний (20), метод, по всей видимости, сохранит свою работоспособность при конечном N, несколько потеряв в точности. При этом свойства асимптотической несмещённости оценок будут утеряны.

Поскольку время запаздывания 9» заранее не известно, минимизацию целевой функции (20) можно провести для различных пробных дискретных времен запаздывания 9», перебираемых из некоторого интервала. Минимум зависимости Ь2(т»), где т» = 9»Д£, будет наблюдаться при истинном времени запаздывания т».

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

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

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

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

Для разделения восстановленных коэффициентов связи полной модели на значимые и незначимые можно использовать другой подход. Отсортируем все коэффициенты по абсолютной величине в порядке убывания. Затем рассчитаем величину

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

Б

(21)

»=1

величину (Ж — — 1) необходима, поскольку иначе разные осцилляторы внесут разный вклад в меру (21) из-за отличающихся времен запаздывания 9^ и, как следствие, разного числа слагаемых в формуле для целевой функции (20). Далее будем добавлять в модель по одному коэффициенту связи в порядке уменьшения их абсолютного значения и снова рассчитывать Л. Наконец, построим зависимость Л от числа М коэффициентов связи в модели, где М = 1,..., — 1).

По мере учета в модели всё большего числа реально существующих связей она будет становиться все более точною. При этом величины (17), характеризующие расстояние между точками реконструируемых нелинейных функций элементов, а вместе с ними и величины (20) и (21) будут уменьшаться. Уменьшение величины Л на графике Л(М) практически остановится при М = Мгеа1, где Мгеа1 - число реальных связей. Последующее добавление в модель незначимых коэффициентов, соответствующих отсутствующим связям, почти не повлияет на величину Л. При этом, для М ^ Мгеа зависимость Л(М) будет оставаться почти постоянной. Графики, наглядно иллюстрирующие такой подход, приведены далее в разделах, демонстрирующих применение метода в численном и радиофизическом эксперименте. Отметим, что этот метод требует существенно больше вычислительных затрат, чем способ разделения восстановленных коэффициентов связи на значимые и незначимые на основе метода К-средних, поскольку также является итерационным.

На практике для контроля точности восстановления априорно неизвестной архитектуры связей в ансамбле можно использовать оба рассмотренных подхода, которые в идеале должны дать одинаковые результаты.

6. Тестирование модифицированного подхода в численном эксперименте

6.1. Восстановление ансамбля, состоящего из связанных уравнений Маккея-Гласса.

В качестве первого примера восстановим параметры элементов и архитектуру связей в ансамбле диффузионно связанных систем Маккея-Гласса [52], описываемых уравнением (1) с функцией (12), который мы уже рассматривали ранее в разделе 3.3. Продемонстрируем преимущества модифицированного подхода перед оригинальным, увеличив количество элементов в ансамбле с 10 до 20. Архитектура случайно выбранных связей в ансамбле из 20 элементов приведена на рис. 8, а. Из 380 возможных связей между элементами ансамбля реально имеется 60, среди которых есть как однонаправленные, так и взаимные. Все элементы ансамбля являлись неиден-

Рис. 8. a - архитектура связей в ансамбле из 20 осцилляторов с запаздыванием ( ) и нелинейною функцией ( ); b - временной ряд переменной xi(i) ансамбля из уравнений Маккея-Гласса при добавлении к измеренному ряду всех элементов нормального белого шума со среднеквадратичным отклонением Oi = 0.004

Fig. 8. a - coupling architecture in the ensemble of 20 time-delayed oscillators ( ) and nonlinear function ( ); b - time series of the variable x1(t) of the ensemble in presence of dynamical noise with Oi = 0.004

L, 10-2

10"6

10-8

200

a

300

400

10-1 ]

10-2

10-3

10-4

10-5 b

Рис. 9. Зависимости L2 (xi) для каждого из 20 элементов ансамбля из уравнений Маккея-Гласса в отсутствие шума (a) и в присутствии нормального шума с Oi = 0.004 (b)

Fig. 9. Dependences L2 (xi) for all 20 elements of the ensemble of Mackey-Glass oscillators in absence of dynamical noise (a) and in presence of the normally distributed white noise with Oi = 0.004 (b)

10-4

тичными. Их параметры принимали случайные значения из равномерного распределения в следующих интервалах: т € [250; 400], б € [7.5; 12.5], а € [0.20; 0.25], € [0.02; 0.06]. При этом все элементы колебались хаотически. Длина временных рядов составляла N = 104 значений при шаге выборки ДЬ = 0.5. Рассмотрим случай отсутствия шума и случай, когда к временному ряду каждого элемента добавлен некоррелированный нормальный шум ^ (Ь) с нулевым средним и среднеквадратичным отклонением о = 0.004. На рис. , Ь приведён зашумленный хаотический временной ряд колебаний первого элемента.

На рис. 9 приведены зависимости целевой функции от пробного времени запаздывания ¿2К) для всех 20 элементов ансамбля в случаях отсутствия и присутствия шумов. Глобальные минимумы всех кривых L2(тг) наблюдаются в точности при истинных временах запаздывания. Величины нормированы на величину ^ — 02 — 1) для удобства сравнения. Отметим, что для оценки временного ряда производных {«г (п)} по временному ряду {«¿(п)} использовалась аппроксимация со сглаживанием параболою.

При реконструкции модельного уравне-_ ния (1) для каждого элемента ансамбля по-

гптт1.........i.........i.........i.........i.........i.........i 4

лучаются 19 ненулевых коэффициентов связи £¿3, часть которых являются лишними. Лишние коэффициенты по модулю на несколько порядков меньше действительных и поэтому их можно выявить с помощью кластеризации в логарифмическом масштабе, как это показано на рис. 10. Видно, что модули всех значений ^ 3 хорошо делятся на два кластера, состоящих из действительных (справа) и лишних (слева) коэффициентов. Наличие шума приводит к сближению кластеров. Очевидно, что, начиная с некоторого критического уровня шума, метод начинает выдавать ошибки, пропуская часть имеющихся связей или находя ложные.

Члены суммы в (17), соответствующие лишним коэффициентам связи, были удалены из модели и реконструкция была проведена

стП018 = 0.004

.= 0

10"'

10

10

10"4

10

10"2 \к j

Рис. 10. Распределение значений модулей оценок коэффициентов связи для всех элементов ансамбля из уравнений Маккея-Гласса в отсутствие шума (снизу) и в присутствии шума (сверху). Значимые коэффициенты показаны черными крестиками, а незначимые - серыми точками

Fig. 10. Distribution of absolute values of estimates of coupling coefficients for all elements of ensemble of Mackey-Glass equations in ebsence of noise (lower) and in its presence (upper). Significant coefficients are plotted with black crosses, and insignificant - with gray circles

повторно уже без них. В результате архитектура связей в ансамбле оказалась восстановлена правильно, в точном соответствии с рис. 8, а, как для случая отсутствия шумов, так и при их наличии. В отсутствие шума точность реконструкции коэффициентов модели составила 4 значащих цифры. Неидеальная точность в таком случае объясняется конечною длиною ряда и необходимостью численной оценки производной. При наличии шума значения коэффициентов также были восстановлены с удовлетворительною точностью. В частности, для первого осциллятора они составили к[ 4 = 0.0467 при истинном значении к1?4 = 0.0475, к'15 = 0.0288 при истинном значении к1,15 = 0.0294, к[ 18 = 0.0287 при истинном значении к1?18 = 0.0292.

Чтобы охарактеризовать точность восстановления коэффициентов связи в среднем по ансамблю, была рассчитана средняя относительная погрешность коэффициентов

Ek =

k

hj

(22)

где усреднение проводилось только для действительных коэффициентов. Величина Е& составила 0.023, колеблясь от 0.002 до 0.058. Параметры инерционности были восстановлены с ещё большею точностью, так, для первого осциллятора = 12.29 при истинном значении е1 = 12.32, аналогичным с Е& образом рассчитанная средняя относительная погрешность реконструкции Ее = 0.007. Большая точность реконструкции параметров ег является, скорее всего, следствием того, что они вносят более уникальный вклад в целевую функцию (20), будучи домножены на базисные функции Ах г (п), в то время как параметры связи стоят перед однотипными базисными функциями вида (Аж^ (п) — Ахг (п)), вклад которых может частично взаимно компенсироваться, особенно при частичной синхронизации в ансамбле.

На рис. 11, а приведена восстановленная нелинейная функция /1 первого элемента для случая отсутствия шума. Она очень точно совпадает с истинной функцией уравнения Маккея-Гласса - настолько, что на одном графике их невозможно различить в использованном масштабе; поэтому для их сопоставления мы вынуждены были отобразить восстановленную функцию более жирным чёрным, а исходную, отрисованную серым, достроить за пределы, в которых она была реконструирована. На рис. 11, Ь приведена чёрным та же нелинейная функция /1, построенная при указанных выше восстановленных значениях и при наличии шума. Видно, что в этом случае она уже заметно отличается от истинной. Точность восстановления параметров и нелинейных функций для остальных элементов ансамбля примерно такая же, как и для первого элемента ансамбля.

Рис. 11. Восстановленная (чёрная линия) и истинная (серая тонкая линия) функция fi (xi) в отсутствие шума (а) и в присутствии нормального шума со среднеквадратичным отклонением Oi = 0.004 (b)

Fig. 11. Nonlinear function fi (xi) - reconstructed (black thick points) and original (gray thin line); a - without noise, and b in presence of normal white noise with oi = 0.004

Рис. 12. Зависимости Л(М) для ансамбля из уравнений Маккея-Гласса в отсутствие шума (серый цвет) и в присутствии нормального шума со среднеквадратичным отклонением о = 0.004 (чёрный цвет)

Fig. 12. Dependence Л(М) for an ensemble of Mackey-Glass equations in absence of noise (gray) and in presence of normal white noise with oi = 0.004 (black)

Очевидно, что лишние и действительные связи вносят различный по абсолютной величине вклад в целевую функцию (20). На этом основании нами был апробирован альтернативный кластеризации методом K-средних алгоритм разделения лишних и необходимых коэффициентов среди всех оценок kj j, основанный на построении зависимости (21) (рис. 12). При M = Mreai величина Л(М) практически достигает своего минимального значения и с увеличением M меняется очень слабо (мене чем на 1%), что видно на графике. В отсутствие шума Л(Мгеа}) очень близка к нулю и составляет 1.3 ■ 10-5, в то время как потеря даже одного коэффициента приводит к её резкому росту до Л(Мгеа[ — 1) = 55.7 ■ 10-5, то есть примерно в 43 раза. Таким образом, оба предложенных подхода дали в рассмотренном примере одинаковый результат и позволили точно реконструировать архитектуру связей.

В реальности значения времён запаздывания т неизвестны, но их можно восстановить предложенным методом. Для этого проще всего использовать простой перебор дискретных значений 9j (точнее, чем с шагом выборки, их восстановить всё равно невозможно). Так как реконструкция для каждого элемента ансамбля проводится отдельно, эта задача может быть решена за относительно приемлемое время. Для примеров, приведённых на рис. 9, в диапазоне пробных времён запаздывания 400 ^ 0j ^ 1000 (601 вариант), что соответствует 200 ^ xj ^ 500, реконструкция занимала не более часа машинного времени на компьютере с процессором Intel Core i5 (для ускорения процесса применялась параллелизация вычислений).

Если внимательно рассмотреть результаты, приведённые на рис. 9, можно видеть, что для всех двадцати элементов ансамбля на кривой имеются очень чётко выраженные минимумы как в случае «чистых» рядов (a), так и при наличии измерительного шума (b), хотя, конечно, во втором случае они более размыты. При этом детальное сопоставление показывает, что во всех случаях минимум соответствует истинному времени запаздывания, то есть ошибка во времени реконструкции в рассмотренном случае меньше шага выборки.

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

6.2. Восстановление ансамбля, состоящего из связанных уравнений Икеды. Чтобы сопоставить результаты двух подходов и дополнительно протестировать общность нового метода, рассмотрим снова ансамбль, описываемый уравнениями (11), где мы предполагаем, что

10"6

1Q-5

10"4

10-3

10-2

10-1 \К\

15 10 5

a

b

0

50

100

150

200

Рис. 13. Результаты применения двух алгоритмов разделения коэффициентов на действительные и ложные для ансамбля из уравнений Икеды с запаздыванием в присутствии измерительного нормального шума со среднеквадратичным отклонением Oi = 0.003: a - распределение значений модулей оценок коэффициентов связи для всех элементов ансамбля (значимые коэффициенты показаны черными крестиками, расположенными внутри окружности, а незначимые - серыми точками, расположенными внутри прямоугольника); b - зависимость Л(М)

Fig. 13. Results of application of two techniques for coupling coefficient separation for an ensemble of Ikeda time-delayed oscillators in presence of measurement noise with Oi = 0.003, where subplot (a) show distribution of absolute values of coupling estimates for all elements of the ensemble (significant coefficients are shown with black crosses inside the ellipse, insignificant ones are shown with gray circles inside the rectangle), and subplot (b) shows the dependence Л(М)

Ег = 1. Пусть ансамбль состоит из 16 осцилляторов с 50 связями между ними (примерно по 3 связи на каждый) из 240 возможных, причём связи заданы случайно. Все элементы ансамбля неидентичны, их параметры задаются случайно из равномерных распределений на следующих интервалах: тг € [2; 5], цг € [15; 25], х0,г € [0; 2п] и € [0.1; 0.5]. Все осцилляторы без связи демонстрировали хаотическую динамику, использовались временные ряды длиною в 10000 отсчётов, как и для других ранее рассмотренных примеров. Интервал выборки составлял, как и ранее для систем Икеды, АЬ = 0.01. Ко временным рядам добавлялся измерительный гауссов независимый шум ^¿(Ь) с нулевым средним и стандартным отклонением аг = 0.003, что составляет 0.1-0.2% от стандартного отклонения сигнала (отношение сигнал/шум порядка 54-60 дБ). Такой уровень шума следует считать, конечно, невысоким.

Поскольку значения Ег были в данном численном эксперименте фиксированы равными единице вследствие свойств системы (11), коэффициенты связи восстанавливались очень точно с ошибками менее 10-4 от их истинного значения. Более интересно рассмотреть результаты разделения связей на действительные и ложные связи, приведённые на рис. 13.

Из рис. 13, а видно, что ложные (серые кружки) и истинные (чёрные крестики) связи легко разделить по абсолютной величине и они прекрасно кластеризуются. Рис. 13, Ь показывает, что на зависимости Л(М) имеется чёткий излом, соответствующий значению М = 50, то есть истинному количеству связей в ансамбле. Таким образом, оба алгоритма позволили верно восстановить архитектуру связей.

Результаты оценки времени запаздывания для всех 16 восстановленных систем Икеды приведены на рис. 14. Как и для ансамбля из 20 систем Маккея-Гласса, рассмотренных в предыдущем подразделе, все времена запаздывания восстанавливаются точно, по крайней мере, с точностью до интервала выборки.

Рис. 14. Зависимости целевой функции от пробного времени запаздывания Lf(xi) для ансамбля из уравнений вида (11) в присутствии нормального шума со среднеквадратичным отклонением Oi = 0.003

Fig. 14. Dependence of target function on trial delay time L2(xi) for an ensemble of (11) oscillators in presence of mesuarement white gaussian noise with Oi = 0.003.

7. Восстановление цепочки экспериментальных радиотехнических генераторов

с запаздывающей обратной связью

Модифицированный метод был применён к экспериментальным временным рядам цепочки, состоящей из В = 10 однонаправленно связанных радиотехнических генераторов с запаздывающей обратной связью. По сравнению с оригинальным подходом число элементов в ансамбле было увеличено с трёх до 10. Каждый генератор представляет собою кольцевую систему, состоящую из линии задержки, нелинейного элемента и низкочастотного ЯС-фильтра первого порядка (см. рис. 7, a). Модельное уравнение, которым описывается г-й элемент такой цепочки, имеет вид

D

RiCiVj(t) = -Vj(t) + /j(V(t - Tj)) + E ki j (Vj(t) - v-(t)).

j=1, j=i

(23)

Уравнение ( ) сводится к виду ( ), если положить ^ = Я С^.

Цепочка состоит из неидентичных элементов, параметры которых принимают значения в следующих интервалах: т € [2.50; 4.75] мс, ^ € [0.203; 0.536] мс, А , ¿+1 € [0.10; 0.23], А = 0 Ч? = г + 1- Все нелинейные элементы имели квадратичную передаточную характеристику. Хаотические сигналы записывались с помощью 10-канального аналого-цифрового преобразователя с частотой выборки Д = 100 кГц.

На рис. 15, a приведены зависимости для всех 10 элементов цепочки. Глобальные

минимумы девяти зависимостей наблюдаются при истинных временах запаздывания ге-

нераторов. Лишь для одного из генераторов минимум оказался смещён относительно истинного времени запаздывания на время, равное одному интервалу выборки.

Для каждого генератора цепочки мы получаем 9 восстановленных коэффициентов связи ^. Из этих 9 коэффициентов для 9 элементов цепочки только один коэффициент является значимым, а для элемента на краю цепочки, на который не действуют другие генераторы, вообще не должно быть А ^. Поскольку число действительных коэффициентов на порядок меньше числа лишних коэффициентов, использование метода К-средних для разделения восстановленных коэффициентов связи на значимые и незначимые оказывается не столь эффективным, как в рассмотренном выше примере. Для успеха алгоритма потребовалось делить всё множество полученных А ^ на 4 кластера, считая значимыми только А ^, принадлежащие верхнему из четырёх. При использовании двух кластеров в архитектуре связей оставались несколько лишних.

L 100

10-1 10-2 ] 10-3 10-4

a

t,, ms

Рис. 15. Зависимости для всех экспериментальных генераторов в цепочке (a) и результаты реконструкции

нелинейной функции первого осциллятора (b)

Fig. 15. Dependence L2 (t^) for all experimental generators in the chain (a) and results of reconstruction of nonlinear function for the first generator (b)

Метод, основанный на использовании меры (21), оказался более эффективным. На рис. 16 построена зависимость Л(М). Видно, что при М = 9 зависимость Л(М) практически достигает минимума и с дальнейшим увеличением М почти не меняется. Оставив только 9 самых больших по модулю восстановленных коэффициентов связи и удалив все остальные незначимые, удалось получить оценки коэффициентов связи и параметров инерционности, близкие к номиналам. Например, для первого генератора восстановленные параметры инерционности и связи имеют следующие значения: е' = 0.204 мс, к' 2 = 0.22. Восстановленная при этих значениях нелинейная функция / приведена на рис. 15, Ь. Она достаточно хорошо совпадает с истинной передаточной характеристикой нелинейного элемента первого генератора. Аналогичные результаты были получены для остальных элементов.

Л 0.10 0.08 0.06 0.04 0.02 0.00

0 20 40 60 80 M

Рис. 16. Зависимость Л(М) для цепочки экспериментальных генераторов с запаздывающей обратной связью

Fig. 16. Dependence Л(М) for the chain of experimental generators with time delay

8. Модель динамики частоты бёрстинга нейронов

Первая основная рассматриваемая модель была предложена в [54] и принадлежит к классу моделей эволюции основной частоты бёрстинга сети нейронов. Каждый узел описывается одним уравнением для трансмембранного потенциала Xi, находящимся под нелинейным внешним воздействием Н(х.) со стороны других нейронов:

Б

Xi = -Xi + к , . Н(х.), (24)

.7 = 1 , 3 =

I = 1, 2,...,В, Л,(х) = Ш(дх),

где д параметр масштабирования; к . - элементы матрицы связанности К, характеризующей эффективность синаптической передачи, и распределённые нормально с нулевым средним и среднеквадратичным отклонением где Л есть параметр, обобщённо характеризующий

силу связей. В работе [54] было аналитически показано, что ансамбль (24) демонстрирует хаотическое поведение в термодинамическом пределе В ^ то при (Лд) > 2 для произвольной типичной матрицы К.

В работе [43] предложен подход к реконструкции ансамбля (24) по временным рядам. Модель (24) может быть слегка модифицирована и обобщена, если рассмотреть возможность наличия запаздывания в связи:

х^г) = -х^)+ £ к.Н(х.(г - Т.)), (25)

3=1, 3=

где т,1. - времена запаздывания, свои для каждого члена связи.

Введение задержек в связях увеличивает сложность модели, что даёт возможность получить хаотическое поведение для меньшего числа узлов.

8.1. Реконструкция модели с запаздывающими связями. Реконструкция модифицированной модели (25) может быть осуществлена примерно тем же способом, и модели без запаздывания, если времена запаздывания известны. Обозначим целочисленные времена запаздывания в шагах выборки как /ДЬ. Тогда уравнение ( ) может быть переписано как

D

Xi(n) = -Xi(n) + ki,jh(xj(n - Qi,j)). (26)

j=1,j=i

На следующем шаге для данного г отсортируем все измеренные значения |жг(п)}^=1, например, в порядке возрастания, но можно и в порядке убывания, это не имеет принципиального значения. Введём отображение Qi аналогично тому, как мы делали это для ансамблей систем с запаздыванием, ставящее в соответствие номеру п значения х^п) в исходном ряде его номер Qi(n) в отсортированном. Также рассмотрим обратное отображение Q-1, ставящее в соответствие номеру Qi (п) в отсортированном ряде номер п в исходном. Таким образом, Q-1(Qi(n)) = п.

Рассмотрим два ближайших значения в отсортированном ряде, имеющие в исходном номера п и рп = Q-1(Qi(п) — 1), пусть в отсортированном ряде х^рп) стоит непосредственно перед х^п). Тогда обозначим разницу между ними как 8^п)

6Дп) = х^п) — Xi(ri,j,n), (27)

i,j,n

= Q-1(Qi(n - 9i,j) - 1).

Также рассмотрим одновременные значения в рядах всех прочих узлов {xj}: Xj (п) и Xj(г^-;П). Можно переписать (27) с использованием (26), выразив 8^п) следующим образом:

D

ôi(n) = Е ki,j&hi,j(n) - ÀXi(n), (28)

j=1,j=i

ÀXi(n) = Xi(n) - Xi(riyjyn), (29)

Àhi,j(n) = h(Xi(n - 9i,j)) - h(Xi(n,j,n)). (30)

Поскольку амплитуда колебаний конечна, значения Xi(n) и Xi(pn), будучи соседними в отсортированном ряде, очень близки друг к другу, если N достаточно велико. Чтобы охарактеризовать эту близость в среднем, рассмотрим сумму квадратов ôi(n):

N N / D \ 2

Si2(ki)= £ ô2(n)= Е Е ki,jÀhi,j(n) - ÀXi(n) . (31)

n=9i+1 n=9i+1 \j=1,j=i /

Далее рассмотрим правую часть уравнения (31). Вследствие хаотического и несинхронного поведения различных узлов близость Xi(n) и Xi(pn) не ведёт к близости одновременных значений Xj (n) и Xj (pn), как и производных Xi (n) и Xi(pn). Поэтому Àhitj (n) и ÀXi(n) в общем случае не являются малыми. Следовательно, единственный путь добиться малости (31) состоит в том, чтобы найти правильные значения ki,j так, чтобы слагаемые ki,jÀhi,j (n) и ÀXi(n) взаимно компенсировались.

Если функции связи h известны, формула (31) может рассматриваться как постановка задачи на наименьшие квадраты, где S2 есть целевая функция, которую следует минимизировать. Здесь ÀXi(n) и Àhi,j (n) могут быть получены по временным рядам, а в результате минимизации S2 может быть получена оценка матрицы связей K.

Рассмотрим, как ведёт себя функция (31) в пределе N ^ то для того, чтобы охарактеризовать асимпотические свойства оценок. В таком случае при условии конечности амплитуды колебаний все 8,(п) спадают в среднем по закону ~ N-1, так что квадраты 82(п) спадают как ~ N-2. Поэтому сумма (3]) должна спадать как ~ N х N-2 = N-1. Это означает, что £2(к) стремится к нулю в пределе N ^ то. Таким образом, оценки ^^ коэффициентов связи к^, полученные представленным методом, асимптотически не смещены.

8.2. Реконструкция времён запаздывания. Чтобы реконструировать времена запаздывания, если они неизвестны, можно провести сканирование всех возможных комбинаций 9^ для каждого г отдельно, как это предлагалось для систем с запаздыванием в [41]. Можно провести реконструкцию для всех пробных запаздываний в заданном диапазоне, и настоящее время запаздывания должно соответствовать минимуму ( ). Поскольку величины могут принимать только неотрицательные целые значения, такой подход вполне возможен, но потребует колоссальных вычислительных мощностей, так как сканировать придётся в (В — 1)-мерном пространстве (для каждого реконструируемого осциллятора реконструкция проводится отдельно, поэтому удастся избежать работы в В(В — 1)-мерном пространстве).

Вместо этого мы предлагаем итерационный алгоритм, похожий на метод градиентного спуска. Сначала для фиксированного г задаются начальные пробные значения 9^. Затем проводится реконструкция для этого вектора, полученное значение целевой функции назовём £20. Далее для одной из связей увеличим дискретное время задержки на 1 и найдем минимум Si,+j функции (31), после чего вернем прежнее значение 9^. Повторим эту процедуру по очереди для каждой связи. Далее найдем минимумы Si-j функции (31) при поочередном уменьшении каждого 9^ на 1. Из (В — 1) минимумов Si,+j и (В — 1) минимумов Si-j найдем наименьший Si)min. Если Si)min < Si)o, то значения , соответствующие Si)min, принимаются за новые стартовые догадки и процедура повторяется, пока не достигнет минимума.

8.3. Численное тестирование алгоритма. Чтобы реконструировать описанные выше сети осцилляторов с запаздывающими связями, необходимо знать все времена запаздывания т^-или их дискретные аналоги 9^. Однако предполагать, что в реальном физическом или биологическом эксперименте они будут доступны или известны из первых принципов, - слишком оптимистично. Поэтому в разделе 8.2 был предложен итеративный подход к реконструкции времён запаздывания, тестируемый ниже. Чтобы изучить сходимость алгоритма, пробная реконструкция проводилась при различных стартовых догадках для времён запаздывания. В рассматриваемом случае при В ^ 16 узлов в сети (24) и от 240 связей и более рассмотреть все возможные комбинации стартовых догадок сложно даже несмотря на то, что их значения дискретны с шагом, равным интервалу выборки. Поэтому мы ограничились весьма реалистичным сценарием, когда стартовые догадки 9о для запаздывания во всех связях принимаются равными. Чтобы охарактеризовать качество реконструкции, для каждого рассмотренного набора 6о рассчитывались две меры:

1. количество уаь3 (и процент уге1) корректно восстановленных времён запаздывания;

2. максимальная из всех ошибка реконструкции времени запаздывания тах(Д9егг). Обе эти меры отложены на рис. 17.

Исходные значения времён запаздывания т^ были сгенерированы из равномерного распределения на отрезке [2.5; 3.5], то есть дискретные времена запаздывания лежали на отрезке 9^ € [50; 70]. Из рис. 17 можно видеть, что существует большая область значений 9о вокруг центра рассматриваемого диапазона (90 = 60), стартуя из которой, значения запаздываний восстанавливаются очень точно: 2 ошибочных значения с отклонением от истинных на ±1Д£ для В = 16 из 240 коэффициентов связи, одна ошибка, равная 2Д£ для d = 24 из 554 времён запаздывания, и 6 ошибок, не превышающих по модулю 6 для В = 32 из 992 связей. Если, однако,

Рис. 17. Максимальная среди всех абсолютная ошибка при оценке времён запаздывания в связях шах(А6еГг) (я); процент уГе1 корректно определённых времён запаздывания (Ь)

Fig. 17. The overall maximal absolute value of the error of coupling time delay estimation max(A6err) (a); and the per cent vrei of precisely detected time delays (b)

реконструкция начинается с очень больших или очень малых стартовых догадок, она проваливается при той же длине временного ряда. Для сети из D = 32 узлов максимальная ошибка max(A9err) соответствовала всегда одной и той же единственной связи для всех стартовых догадок в широком диапазоне от 9о = 32 до 6о = 69, в то время как остальные запаздывания восстанавливались одинаково хорошо для всего диапазона. Эта связь соответствует относительно малому коэффициенту связи по модулю в 2.5 раза меньше среднего значения.

Чтобы убедиться в том, что полученные результаты не являются частным случаем для удачно выбранного нами ансамбля, времена запаздывания были восстановлены для 20 различных ансамблей при D = 16, для которых режим был хаотическим. При этом мы ограничились случаем 9о = 60 для всех ансамблей и всех осцилляторов. В результате реконструкции из 240 времён запаздывания для первого ансамбля верно оказались восстановлены 237, ещё в трёх были допущены ошибки A9err = т'/At = 1 для осциллятора № 9, A9err = 3 для осциллятора № 6 и A9err = 4 для осциллятора № 10.

Общие результаты по 20 ансамблям при длине ряда в 212 точек приведены на рис. 18, a в виде гистограммы. Всего из 5000 времён запаздывания неверно были оценены 55 (1.1%), в том числе с ошибкою в ±1 время дискретизации - 35 (0.7%). Для трёх ансамблей все времена запаздывания были восстановлены точно.

0.004

0.003

0.002

0.001

0.000

_L

_L

_L

_L

_L

-15 -10 -5 0 5 10 15 At/At a b

Рис. 18. a - гистограмма плотности распределения ошибок оценок времён запаздывания (столбец, соответствующий точной оценке при равной нулю ошибке, вырезан, поскольку его высота слишком велика для наглядного представления других столбцов); b - усреднённая по всем 16 осцилляторам и по всем 20 ансамблям нормированная зависимость значения целевой функции (31) и его стандартного отклонения от длины временного ряда N

Fig. 18. a - histogram of probability dencity of errors in delay time reconstruction. The column corresponding to exact estimate (no error) is cut due to its height is too large to be fitted to the plot. b - the averaged over all 20 ensembles and 16 nodes in every ensemble normalized dependence of target function (31) value on time series length with its standard deviation

Результаты реконструкции задержек в связях сильно зависят от динамического режима. Для рис. 17 были выбраны только хаотические режимы с достаточно большим А1 > 0.12 старшим ляпуновским показателем. Для меньших А1 число ошибок может быть существенно больше, например, для ансамбля из В = 24 осцилляторов с А1 ~ 0.04 все коэффициенты связи могут быть корректно восстановлены с соответствующими им временами запаздывания только для 6 осцилляторов из 24, а число запаздываний, восстановленных с ошибками, превышает 60% даже для стартовых догадок 90 = 60.

8.4. Асимптотическая точность метода. При описании метода было показано, что он должен давать асимптотически точные оценки при N ^ то. Чтобы подтвердить это в численном эксперименте, следует рассмотреть все существующие источники ошибок. Их четыре:

1. ошибки, связанные с недостаточно точным представлением экспериментальных данных

(квантование, шумы измерения и проч.);

2. ошибки, возникающие при численном расчёте производной;

3. ошибки, обусловленные конечной точностью расчётов;

4. ошибки, вызванные неточностью метода при конечном N.

Чтобы показать, что метод действительно является асимптотически точным, необходимо понять, как ведут себя все типы ошибок при N ^ то. В натурном эксперименте при наличии шумов измерений основной вклад дают ошибки первых двух типов, их вклад в целевую функцию (31) пропорционален N. В численном моделировании ошибкою первого типа можно пренебречь, если числа записаны с максимальною (родною) точностью, так как она заведомо не превышает погрешностей, возникающих при расчётах (тип 3).

Ошибкою второго типа можно пренебречь при стремлении шага выборки к нулю ^ 0. Однако на практике даже в численном эксперименте это невозможно, поэтому при исследовании асимптотических свойств метода мы дополнительно записали временные ряды производных для всех элементов ансамбля и использовали их вместо рядов производных, оцененных численно. Влияние, которое численная оценка производной оказывает на точность восстановления коэффициентов, очень невелико, но становится заметным при N > 213.

Ошибки, обусловленные конечной точностью расчётов, можно оценить теоретически, но это довольно сложно, так как используемые алгоритмы, реализованные в популярных библиотеках вроде Ьараек или питру имеют множество оптимизаций по точности. В любом случае, из-за матричных операций они имеют порядок не менее N3/2. Оценить эти ошибки можно, если использовать различную точность расчётов. К сожалению, в настоящее время все основные библиотеки методов ограничиваются одинарною (32 бита) и двойною (64 бита) точностями при операциях с плавающей запятою. Численный эксперимент показал, что при одинарной точности даже при N = 210 ошибка точности расчётов вносит заметный вклад. В то же время при использовании двойной точности надёжно выявить её наличие при длинах ряда до N = 215 включительно не удалось.

Чтобы показать асимптотическую сходимость метода, была проведена реконструкция 20 ансамблей из 16 осцилляторов по хаотическим временным рядам различной длины: от N = 210 (данная длина как минимальная была обозначена N0) до N = 215. Поскольку из-за различий во временах запаздывания и коэффициентах связи итоговые значения целевой функции £2г (г соответствует номеру осциллятора, I - номеру ансамбля) сильно варьировали, они были отнорми-рованы на значение, соответствующее минимальной рассмотренной длине ряда £?г(N0), а затем усреднены по всем осцилляторам во всех ансамблях (см. 18, б).

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

что при выбранном способе нормировки для N = N0 разброс будет нулевым по определению. В то же время, падение, которое должно быть линейным по нашей гипотезе, не является таковым и замедляется при больших N, что может быть обусловлено конечной точностью расчётов.

9. Подход для случая произвольных сигмоидальных функций связи

9.1. Модернизация метода на случай неизвестных функций связи. Отметим два важных недостатка описанного выше подхода, препятствующих его практическому применению. Во-первых, уравнения для отдельных узлов либо слишком просты (24), слишком специфичны и должны быть известны априорно. Во-вторых, функции связи предполагаются линейными, как в рассмотренной в главе о системах с запаздыванием модели (1) или известными заранее, как ранее в данной главе. Далее рассматривается задача о реконструкции ансамбля вида (32), являющегося обобщением (24) и (25) и содержащего произвольные непрерывные нелинейные собственные функции /¿(ж) для каждого узла и сигмоидные функции связи с параметрами дг,у (масштабный коэффициент) и <С,г,з (сдвиг), входящие в уравнения ( ) нелинейно.

Б

Хг = /¿(жг) + ^ кг,у Л (д^- (ж у (£ - т^) - ^)) • (32)

г=1,г=3

Обозначим максимальное из 6г,3- для данного г как вг и выразим нелинейную функцию /¿(ж) из (32) Уп = 6г + 1, • • •, N:

Б

/г(жг) = Хг - Е кгу Ш (дг>у(ж3(п - вг,з) - ^^ • (33)

г=1,г=3

Поскольку функция /г по определению непрерывна, близким значениям Хг с номерами п и рп в отсортированном ряде соответствуют близкие же значения функции /¿(ж(п)) и /г(ж(рп)) соответственно. Теперь разницу между ними обозначим как 8г(п):

Б

8г(п) = /г(ж(п)) - /г(ж(рп)) = Ажг(п) - ^ АНг,у (П - бгу), (34)

3 = 1,3=1

АХ ¿(п) = ж г (п) - ж г(рп),

АНг,з (п) = кг,з th (дг,у ж у (п - ) - ^) - кг,у th (дг,у ж у (рп - ) - 1г,у) • (35)

Аналогично рассмотренным ранее уравнениям, в частности, (28), величина 8г(п) в уравнении (34) будет обращаться в 0 в пределе N ^ то только если будут верно выбраны все параметры дг,3, и кгз в уравнении (35). Поэтому по-прежнему сумму 8г(п)2 можно рассматривать как целевую функцию £2, являющуюся функцией всех дг з, и кг,у:

Б

£(кг,1, • • • ,кг,Б,дг,1, • • • ,дг,Б,, • • •,Х,г,п) = Е 8г(п) (36)

3=1,3=г

Минимизация (36) даёт возможность вычислить значения дг,у, и кг,у, после чего функции /г определяются таблично из формулы (33). Минимизацию необходимо проводить одним из существующих специализированных подходов, например, методом Левенберга-Марквардта [55, 56] или методом отображения доверительных областей [57], в численном эксперименте были успешно опробованы оба подхода.

9.2. Тестирование модернизированного метода. В численном эксперименте рассматривались ансамбли с кубическою функцией нелинейности вида

/г(ж) = УгЖ + а^ж3. (37)

Временной ряд первого осциллятора приведён на рис. 19, a. Параметры ^г , уг, $г и аг получались как случайные числа из равномерного распределения на отрезках [-1; 1], [0.5; 1.5], [0.5; 1.5] и [0; 1], соответственно, значения А , 3 - из нормального распределения с нулевым средним и среднеквадратичным отклонением Уравнения решались методом Рунге-Кутты 4-го порядка с адаптивным шагом и интервалом выборки Д£ = 0.05. Результаты реконструкции нелинейной функции первого осциллятора и одной из его функций связи приведены на рис. 19, Ь и рис. 19, c, соответственно.

Рис. 19. a - временной ряд колебаний первого осциллятора ансамбля ( ) из 8 узлов; b - восстановленная (серый цвет) и исходная (чёрный) нелинейная функция этого осциллятора; c - восстановленная (серый) и исходная (чёрный) функция связи, соответствующая воздействию второго осциллятора на первый; d - схема архитектуры связей: чёрными квадратами показаны реально существующие связи, а белыми — отсутствующие

Fig. 19. a - time series from oscillator No. 1 from the ensemble ( ) with 8 nodes; b - original (black) and reconstructed (thick gray) nonlinear function of this oscillator; c - original (black) and reconstructed (gray) coupling function for coupling from oscillator No. 2 to oscillator No. 1; d - coupling architecture, where black squares show rightly detected couplings and white squares show couplings correctly considered to be absent

В реальных сетях как естественного, так и искусственного происхождения далеко не все узлы связаны между собою. В главе 1 были предложены алгоритмы выявления и удаления незначимых (лишних) связей, основанные на кластеризации коэффициентов ансамбля методом К-средних. Эта процедура оказалась хорошо работоспособна и для предлагаемого обобщённого алгоритма. На рис. 19, d приведена архитектура связей в ансамбле: все связи восстановились верно. Всего было рассмотрено 14 ансамблей, из них в 7 случаях архитектура была восстановлена без ошибок, в 4 случаях имелась 1 лишняя связь, в остальных 3 случаях имелись как лишние, так и пропущенные связи.

Поскольку в реальности значения времён запаздывания неизвестны или могут быть определены лишь примерно, предлагаемую методику можно модифицировать следующим образом. Вместо истинных вг , 3 используются некоторые догадки вг ^ и вся описанная процедура проводится для них при фиксированном г. Поскольку вг , 3 суть по определению целые числа, реконструкция проводится дополнительно 2(В — 1) раз для наборов значений вг ^ таких, что в каждом наборе только одно вг ^ отличается от исходных догадок на 1 в плюс или минус (исходный набор в таком случае можно считать «центральным»). Затем выбирается такой набор

. -.я

вг Л , для которого целевая функция (36) минимальна. Если это исходный набор, работа

прекращается и полученные значения объявляются оценками времён запаздывания. Если же минимум соответствует какому-либо смещённому набору \ вг Л , то теперь он выбирается в

качестве центрального (догадки) и процедура повторяется уже вокруг него. Такой подход можно рассматривать как адаптацию метода градиентного спуска к дискретному случаю.

Предложенный алгоритм оказывается достаточно эффективен, если исходные догадки лежат вблизи истинных значений, например, отличаются от них на 5% по абсолютной величине. При больших ошибках в исходных догадках для вг 3 часть времён запаздывания может определяться с погрешностями, поскольку алгоритм сходится к локальному, а не глобальному минимуму функции (36). Гораздо более скромная сходимость модифицированного подхода по сравнению с исходным вызвана, по-видимому, использованием нелинейного МНК. Исправить такие ошибки можно, сканируя область догадок с некоторым шагом, либо разрешив смещение по догадкам сразу на 2 шага по вг 3. Оба подхода существенно увеличивают вычислительные затраты, но могут быть эффективно распараллелены.

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

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

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

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

В предложенных подходах важную роль играют методы реконструкции времени запаздывания. Фактически, их можно подразделить на три класса. Первый - это методы, перенесённые из подходов, ориентированных на реконструкцию уединённых осцилляторов. Такие подходы могут хорошо работать, если связь слабая и статистические свойства сигнала не затрагиваются радикально. Для систем, где запаздывание входит в функцию связи, они неприменимы. Второй класс - это подходы, основанные на внешнем активном воздействии. Такие методы, как правило, очень эффективны, но далеко не всегда применимы на практике. Оба этих класса подразумевают, что время (или времена) запаздывания определяются до того, как реконструируются все остальные параметры и функции. Третий класс - это методы, основанные на множественной реконструкции с пробными временами запаздывания. В этом случае критерием правильности служит всё та же целевая функция, что используется для оценки остальных параметров. Однако оценка времени запаздывания всё равно производится отдельно (только уже после реконструкции всех параметров, а не до неё), поскольку оно может быть детектировано только с точностью до шага выборки, то есть оно фактически дискретно. Для одного-двух времён запаздывания задача решается перебором. Для большего числа приходится использовать методы, основанные на пошаговой минимизации с дискретными шагами, что порождает проблему локальных экстремумов. Хотя третий класс методов не самый эффективный, дискретность времени запаздывания позволяет надеяться, что по мере совершенствования компьютеров возможности его применения будут всё более расширяться.

Для некоторых изложенных подходов нами показана асимптотическая точность получаемых догадок с увеличением длины ряда. С одной стороны, такое свойство является важным теоретическим обоснованием для их применения, а также страховкой на случай увеличения числа осцилляторов в ансамбле (можно надеяться на успех, просто увеличив время наблюдения). С другой стороны, зависимость результатов от ошибок счёта при больших размерах выборки и ансамбля и недоступность длинных стационарных рядов для большого числа объектов моделирования (в первую очередь, для биологических систем) указывают на то, что асимптотическая точность может быть избыточным требованием, не всегда существенным при практическом применении.

В ряде случаев нами применялся нелинейный метод наименьших квадратов. Хотя многие исследователи стараются избежать его, по возможности, либо максимально сократить число нелинейно входящих параметров [58], полученный нами опыт свидетельствует, что в рассмотренном классе задач можно рассчитывать на успех, даже имея значительное (до нескольких десятков) число нелинейно входящих параметров. По всей видимости, такая ситуация обусловлена типом рассмотренных функций - сигмоид, являющихся ограниченными на числовой оси, что означает малую чувствительность к ошибкам и большие области сходимости в глобальный и существующие локальные экстремумы. Однако чёткого ответа на вопрос, сколько и каких нелинейно входящих параметров может быть в задаче о реконструкции системы по временным рядам, чтобы рассчитывать на успех, у нас нет. Полагаем, что данный вопрос заслуживает отдельного рассмотрения в будущем.

Библиографический список

1. 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, no. 1. P. 012902.

2. BuriC N., Vasovic N. Global stability of synchronization between delay-differential systems with generalized diffusive coupling // Chaos, Solitons & Fractals. 2007. Vol. 31, no. 2. P. 336-342.

3. Bindu M. Krishna and Manu P. John and Nandakumaran V.M. Multi-user bidirectional communication using isochronal synchronisation of array of chaotic directly modulated semiconductor lasers // Physics Letters A. 2010. Vol. 374, no. 17. P. 1835-1842.

4. Mincheva M., Roussel M.R. Graph-theoretic methods for the analysis of chemical and biochemical networks. II. Oscillations in networks with delays // Journal of Mathematical Biology. 2007. Vol. 55, no. 1. P. 87-104.

5. Kuang Y. Delay Differential Equations with Applications in Population Dynamics. Boston: Academic Press, 1993. 398 p.

6. Bocharov G.A., Fathalla A.R. Numerical modelling in biosciences using delay differential equations // Journal of Computational and Applied Mathematics. 2000. Vol. 125, no. 1. P. 183-199.

7. Orosz G., Moehlis J., Murray R.M.Controlling biological networks by time-delayed signals // Philosophical Transactions of the Royal Society of London A: Mathematical, Physical and Engineering Sciences. 2010. Vol. 368, no. 1911. P. 439-454.

8. Glyzin S.D., Marushkina E.A. Complicated Dynamic Regimes in a Neural Network of Three Oscillators with a Delayed Broadcast Connection // Automatic Control and Computer Sciences. 2018. Vol. 52, no. 7. P. 885-893.

9. Glyzin S.D., Kolesov A.Yu., Rozov N.Kh. Self-Sustained Relaxation Oscillations in Time-Delay Neural Systems // Journal of Physics: Conference Series. 2016. Vol. 727, no. 1. P. 012004.

10. Glyzin S., Goryunov V., Kolesov A. Spatially inhomogeneous modes of logistic differential equation with delay and small diffusion in a flat area // Lobachevskii Journal of Mathematics. 2017. Vol. 38, no. 5. P. 898-905.

11. Packard N., Crutchfield J., Farmer J., Shaw R. Geometry from a Time Series // Phys. Rev. Lett. 1980. Vol. 45. P. 712-716.

12. Fowler A.C., Kember G. Delay recognition in chaotic time series // Physics Letters A. 1993. Vol. 175, no. 6. P. 402-408.

13. Hegger R., BUnner M.J., Kantz H., Giaquinta A. Identifying and Modeling Delay Feedback Systems // Phys. Rev. Lett. 1998. Vol. 81, no. 3. P. 558-561.

14. Bunner M.J., Ciofini M, Giaquinta A., Hegger R., Kantz H., Meucci R., Politi A. Reconstruction of systems with delayed feedback: I. Theory // The European Physical Journal D. 2000. Vol. 10, no. 2. P. 165-176.

15. Yu-Chu Tian, Furong Gao. Extraction of delay information from chaotic time series based on information entropy // Physica D: Nonlinear Phenomena. 1997. Vol. 108, no. 1. P. 113-118.

16. Bunner M.J., Meyer Th., Kittel A., Parisi J. Recovery of the time-evolution equation of time-delay systems from time series // Phys. Rev. E. 1997. Vol. 56, no. 5. P. 5083-5089.

17. Voss H., Kurths J. Reconstruction of non-linear time delay models from data by the use of optimal transformations // Physics Letters A. 1997. Vol. 234, no. 5. P. 336-344.

18. Ellner S.P., Kendall B.E., Wood S.N., McCauley E, Briggs Ch.J. Inferring mechanism from time-series data: Delay-differential equations // Physica D: Nonlinear Phenomena. 1997. Vol. 110, no. 3. P. 182-194.

19. Prokhorov M.D., Ponomarenko V.I., Karavaev A.S., Bezruchko B.P. Reconstruction of time-delayed feedback systems from time series // Physica D: Nonlinear Phenomena. 2005. Vol. 203, no. 3. P. 209-223.

20. Prokhorov M.D., Ponomarenko V.I., Khorev V.S. Recovery of delay time from time series based on the nearest neighbor method // Physics Letters A. 2013. Vol. 377, no. 43. P. 3106-3111.

21. Udaltsov VS., Larger L., Goedgebuer J.P., Locquet A., Citrin D.S. Time delay identification in chaotic cryptosystems ruled by delay-differential equations // J. Opt. Technol. 2005, no. 5. P. 373-377.

22. Zunino L., Soriano M.C., Fischer I., Rosso O.A., Mirasso C.R. Permutation-information-theory approach to unveil delay dynamics from time-series analysis // Phys. Rev. E. 2010. Vol. 82, no. 4. P. 046212.

23. Horbelt W., Timmer J., Voss H.U. Parameter estimation in nonlinear delayed feedback systems from noisy data // Physics Letters A. 2002. Vol. 299, no. 5. P. 513-521.

24. Dai Chaohua, Chen Weirong, Li Lixiang, Zhu Yunfang, Yang Yixian. Seeker optimization algorithm for parameter estimation of time-delay chaotic systems // Phys. Rev. E. 2011. Vol. 83, no. 3. P. 036203.

25. Sorrentino F. Identification of delays and discontinuity points of unknown systems by using synchronization of chaos // Phys. Rev. E. 2010. Vol. 81, no. 6. P. 066218.

26. Ma Huanfei, Xu Bing, Lin Wei, Feng Jianfeng. Adaptive identification of time delays in nonlinear dynamical models // Phys. Rev. E. 2010. Vol. 82, no. 6. P. 066210.

27. Siefert M. Practical criterion for delay estimation using random perturbations // Phys. Rev. E. 2007. Vol. 76, no. 2. P. 026215.

28. Yu Dongchuan, Frasca Mattia, Liu Fang. Control-based method to identify underlying delays of a nonlinear dynamical system // Phys. Rev. E. 2008. Vol. 78, no. 4. P. 046209.

29. Ponomarenko V.I., Prokhorov M.D. Recovery of systems with a linear filter and nonlinear delay feedback in periodic regimes // Phys. Rev. E. 2008. Vol. 78, no. 6. P. 066207.

30. Prokhorov M.D., Ponomarenko V.I. Reconstruction of time-delay systems using small impulsive disturbances // Phys. Rev. E. 2009. Vol. 80, no. 6. P. 066206.

31. Prokhorov M.D., Ponomarenko V.I. Estimation of coupling between time-delay systems from time series // Phys. Rev. E. 2005. Vol. 72. P. 016210.

32. Afraimovich V.S., Nekorkin V.I., Osipov G.V., Shalfeev V.D. Stability, Structures and Chaos in Nonlinear Synchronization Networks. WORLD SCIENTIFIC, 1995.

33. Пиковский А., Розенблюм М., Куртс Ю. Синхронизация. Фундаментальное нелинейное явление. M.: Техносфера. НИЦ «Регулярная и хаотическая динамика», 2003. 496 p.

34. Xiaoming Wu, Zhiyong Sun, Feng Liang, Changbin Yu. Online estimation of unknown delays and parameters in uncertain time delayed dynamical complex networks via adaptive observer // Nonlinear Dynamics. 2013. Vol. 73, no. 3. P. 1753-1768.

35. Wang W.X., Yang R., Lai Y.C., Kovanis V., Grebogi C. Predicting catastrophes in nonlinear dynamical systems by compressive sensing // Phys. Rev. Lett. 2011. Vol. 106. P. 154101.

36. Han X., Shen Z, Wang W.-X., Di Z. Robust Reconstruction of Complex Networks from Sparse Data // Phys. Rev. Lett. 2015. Vol. 114. P. 28701.

37. Brunton S.L., Proctor J.L., KutzJ.N. Discovering governing equations from data by sparse identification of nonlinear dynamical systems // Proc. Natl. Acad. Sci. U.S.A. 2016. Vol. 113. P. 3932-7.

38. Mangan N.M., Brunton S.L., Proctor J.L., Kutz J.N.Inferring biological networks by sparse identification of nonlinear dynamics// IEEE Trans. Mol. Biol. Multi-Scale Commun. 2016. Vol. 2. P. 52-63.

39. Jose Casadiego, Mor Nitzan, Sarah Hallerberg, Marc Timme. Model-free inference of direct network interactions from nonlinear collective dynamics // Nature Communications. 2017. Vol. 8. P. 2192.

40. Sysoev I.V., Ponomarenko V.I., Prokhorov M.D., Bezruchko B.P. Reconstruction of ensembles of coupled time-delay system from time series // Phys. Rev. E. 2014. Vol. 89. P. 062911.

41. Sysoev I.V., Ponomarenko V.I., Kulminsky D.D., Prokhorov M.D. Recovery of couplings and parameters of elements in networks of time-delay systems from time series//Phys. Rev. E. 2016. Vol. 94. P. 052207.

42. Сысоев И.В., Пономаренко В.И. Реконструкция матрицы связей ансамбля идентичных нейропо-добных осцилляторов с запаздыванием в связи // Нелинейная динамика. 2016. Vol. 12, no. 4. P. 567-576.

43. Sysoev I.V., Ponomarenko V.I., Pikovsky A. Reconstruction of coupling architecture of neural field networks from vector time series // Commun. Nonlinear Sci. Numer. Simulat. 2018. Vol. 57. P. 342-351.

44. Сысоев И.В., Пономаренко В.И., Прохоров М.Д. Реконструкция ансамблей осцилляторов с нелинейными запаздывающими связями // Письма в ЖТФ. 2018. Vol. 44, no. 22. P. 57--64.

45. Sysoev I.V., Ponomarenko V.I., Prokhorov M.D. Reconstruction of ensembles of nonlinear neurooscilla-tors with sigmoid coupling function // Nonlinear Dynamics. 2019. Vol. 95, no. 3. P. 2103-2116.

46. Sysoev I.V. Reconstruction of ensembles of generalized Van der Pol oscillators from vector time series // Physica D. 2018. Vol. 384-385, no. 1. P. 1-11.

47. Savitzky A., Golay M.J.E. Smoothing and Differentiation of Data by Simplified Least Squares Procedures // Analytical Chemistry. 1964. Vol. 38, no. 8. P. 1627-1639.

48. Nelder J.A., Mead R. A simplex for function minimization // Computer Journal. 1965. Vol. 7. P. 308-313.

49. Kendall M, Stuart A. The Advanced Theory of Statistics. New York: MacMillan, 1979.

50. Johnson N.L., Kotz S., Balakrishnan N. Continuous Univariate Distributions. New York: Wiley, 1995. Vol. 2. 752 p.

51. Ikeda K. Multiple-valued stationary state and its instability of the transmitted light by a ring cavity system // Optics Communications. 1979. Vol. 30, no. 2. P. 257-261.

52. Mackey M.C., Glass L. Oscillation and chaos in physiological control systems // Science. 1977. Vol. 197, no. 4300. P. 287-289.

53. Мандель И.Д. Кластерный анализ. М.: Финансы и статистика, 1988. 176 p.

54. Sompolinsky H., Crisanti A., Sommers H.E. Chaos in random neural networks // Phys. Rev. Lett. 1988. Vol. 61, no. 3. P. 259-262.

55. Levenberg K. A method for the solution of certain non-linear problems in least squares // Quarterly of Applied Mathematics. 1944. Vol. 2. P. 164-168.

56. Marquardt D. An algorithm for least-squares estimation of nonlinear parameters // SIAM Journal on Applied Mathematics. 1963. Vol. 11, no. 2. P. 431-441.

57. Coleman T.F., Li Y. An interior trust region approach for nonlinear minimization subject to bounds // SIAM J. Opt. 1996. Vol. 6. P. 418-445.

58. Bezruchko B.P., Smirnov D.A. Constructing nonautonomous differential equations from experimental time series // Phys. Rev. E. Vol. 63, no. 1. P. 016207.

References

1. Heiligenthal S., Jiingling 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, no. 1, p. 012902.

2. Buric N., Vasovic N. Global stability of synchronization between delay-differential systems with generalized diffusive coupling. Chaos, Solitons & Fractals, 2007, vol. 31, no. 2, pp. 336-342.

3. Bindu M. Krishna and Manu P. John and Nandakumaran V.M. Multi-user bidirectional communication using isochronal synchronisation of array of chaotic directly modulated semiconductor lasers. Physics Letters A, 2010, vol. 374, no. 17, pp. 1835-1842.

4. Mincheva M., Roussel M.R. Graph-theoretic methods for the analysis of chemical and biochemical networks. II. Oscillations in networks with delays. Journal of Mathematical Biology, 2007, vol. 55, no. 1, pp. 87-104.

5. Kuang Y. Delay Differential Equations with Applications in Population Dynamics. Boston: Academic Press, 1993. 398 p.

6. Bocharov G.A., Fathalla A.R. Numerical modelling in biosciences using delay differential equations. Journal of Computational and Applied Mathematics, 2000, vol. 125, no. 1, pp. 183-199.

7. Orosz G., Moehlis J., Murray R.M. Controlling biological networks by time-delayed signals. Philosophical Transactions of the Royal Society of London A: Mathematical, Physical and Engineering Sciences, 2010, vol. 368, no. 1911, pp. 439-454.

8. Glyzin S.D., Marushkina E.A. Complicated Dynamic Regimes in a Neural Network of Three Oscillators with a Delayed Broadcast Connection. Automatic Control and Computer Sciences, 2018, vol. 52, no. 7, pp. 885-893.

9. Glyzin S.D., Kolesov A.Yu., Rozov N.Kh. Self-Sustained Relaxation Oscillations in Time-Delay Neural Systems. Journal of Physics: Conference Series, 2016, vol. 727, no. 1, p. 012004.

10. Glyzin S., Goryunov V., Kolesov A. Spatially inhomogeneous modes of logistic differential equation with delay and small diffusion in a flat area. Lobachevskii Journal of Mathematics, 2017, vol. 38, no. 5, pp. 898-905.

11. Packard N., Crutchfield J., Farmer J., Shaw R. Geometry from a Time Series. Phys. Rev. Lett., 1980, vol. 45, pp. 712-716.

12. Fowler A.C., Kember G. Delay recognition in chaotic time series. Physics Letters A, 1993, vol. 175, no. 6, pp. 402-408.

13. Hegger R., BUnner M.J., Kantz H., Giaquinta A. Identifying and Modeling Delay Feedback Systems. Phys. Rev. Lett., 1998, vol. 81, no. 3, pp. 558-561.

14. BUnner M.J., Ciofini M., Giaquinta A., Hegger R., Kantz H., Meucci R., Politi A. Reconstruction of systems with delayed feedback: I. Theory. The European Physical Journal D, 2000, vol. 10, no. 2, pp. 165-176.

15. Yu-Chu Tian, Furong Gao. Extraction of delay information from chaotic time series based on information entropy. Physica D: Nonlinear Phenomena, 1997, vol. 108, no. 1, pp. 113-118.

16. Bünner M.J., Meyer Th., Kittel A., Parisi J. Recovery of the time-evolution equation of time-delay systems from time series. Phys. Rev. E, 1997, vol. 56, no. 5, pp. 5083-5089.

17. Voss H., Kurths J. Reconstruction of non-linear time delay models from data by the use of optimal transformations. Physics Letters A, 1997, vol. 234, no. 5, pp. 336-344ro

18. Ellner S.P., Kendall B.E., Wood S.N., McCauley E., Briggs Ch.J. Inferring mechanism from time-series data: Delay-differential equations. Physica D: Nonlinear Phenomena, 1997, vol. 110, no. 3, pp. 182-194.

19. Prokhorov M.D., Ponomarenko V.I., Karavaev A.S., Bezruchko B.P. Reconstruction of time-delayed feedback systems from time series. Physica D: Nonlinear Phenomena, 2005, vol. 203, no. 3, pp. 209-223.

20. Prokhorov M.D., Ponomarenko V.I., Khorev V.S. Recovery of delay time from time series based on the nearest neighbor method. Physics Letters A, 2013, vol. 377, no. 43, pp. 3106-3111.

21. Udaltsov V.S., Larger L., Goedgebuer J.P., Locquet A., Citrin D.S. Time delay identification in chaotic cryptosystems ruled by delay-differential equations. J. Opt. Technol., 2005, no. 5, pp. 373-377.

22. Zunino L., Soriano M.C., Fischer I., Rosso O.A., Mirasso C.R. Permutation-information-theory approach to unveil delay dynamics from time-series analysis. Phys. Rev. E, 2010, vol. 82, no. 4, p. 046212.

23. Horbelt W., Timmer J., Voss H.U. Parameter estimation in nonlinear delayed feedback systems from noisy data. Physics Letters A, 2002, vol. 299, no. 5, pp. 513-521.

24. Dai Chaohua, Chen Weirong, Li Lixiang, Zhu Yunfang, Yang Yixian. Seeker optimization algorithm for parameter estimation of time-delay chaotic systems. Phys. Rev. E, 2011, vol. 83, no. 3, p. 036203.

25. Sorrentino F. Identification of delays and discontinuity points of unknown systems by using synchronization of chaos. Phys. Rev. E, 2010, vol. 81, no. 6, p. 066218.

26. Ma Huanfei, Xu Bing, Lin Wei, Feng Jianfeng. Adaptive identification of time delays in nonlinear dynamical models. Phys. Rev. E, 2010, vol. 82, no. 6, p. 066210.

27. Siefert M. Practical criterion for delay estimation using random perturbations. Phys. Rev. E, 2007, vol. 76, no. 2, p. 026215.

28. Yu Dongchuan, Frasca Mattia, Liu Fang. Control-based method to identify underlying delays of a nonlinear dynamical system. Phys. Rev. E, 2008, vol. 78, no. 4, p. 046209.

29. Ponomarenko V.I., Prokhorov M.D. Recovery of systems with a linear filter and nonlinear delay feedback in periodic regimes. Phys. Rev. E, 2008, vol. 78, no. 6, p. 066207.

30. Prokhorov M.D., Ponomarenko V.I. Reconstruction of time-delay systems using small impulsive disturbances. Phys. Rev. E, 2009, vol. 80, no. 6, p. 066206.

31. Prokhorov M.D., Ponomarenko V.I. Estimation of coupling between time-delay systems from time series. Phys. Rev. E, 2005, vol. 72, p. 016210.

32. Afraimovich V.S., Nekorkin V.I., Osipov G.V., Shalfeev V.D. Stability, Structures and Chaos in Nonlinear Synchronization Networks. WORLD SCIENTIFIC, 1995.

33. Pikovsky A., Rosenblum M. and Kurths J. Synchronization: A Universal Concept in Nonlinear Sciences. London: Cambridge University Press, 2001.

34. Xiaoming Wu, Zhiyong Sun, Feng Liang, Changbin Yu. Online estimation of unknown delays and parameters in uncertain time delayed dynamical complex networks via adaptive observer. Nonlinear Dynamics, 2013, vol. 73, no. 3, pp. 1753-1768.

35. Wang W.X., Yang R., Lai Y.C., Kovanis V., Grebogi C. Predicting catastrophes in nonlinear dynamical systems by compressive sensing. Phys. Rev. Lett., 2011, vol. 106, p. 154101.

36. Han X., Shen Z., Wang W.-X., Di Z. Robust Reconstruction of Complex Networks from Sparse Data. Phys. Rev. Lett., 2015, vol. 114, p. 28701.

37. Brunton S.L., Proctor J.L., Kutz J.N. Discovering governing equations from data by sparse identification of nonlinear dynamical systems. Proc. Natl. Acad. Sci. U.S.A., 2016, vol. 113, p. 3932.

38. Mangan N.M., Brunton S.L., Proctor J.L., Kutz J.N. Inferring biological networks by sparse identification of nonlinear dynamics. IEEE Trans. Mol. Biol. Multi-Scale Commun., 2016, vol. 2, pp. 52-63.

39. Jose Casadiego, Mor Nitzan, Sarah Hallerberg, Marc Timme. Model-free inference of direct network interactions from nonlinear collective dynamics. Nature Communications, 2017, vol. 8, p. 2192.

40. Sysoev I.V., Ponomarenko V.I., Prokhorov M.D., Bezruchko B.P. Reconstruction of ensembles of coupled time-delay system from time series. Phys. Rev. E, 2014, vol. 89, p. 062911.

41. Sysoev I.V., Ponomarenko V.I., Kulminsky D.D., Prokhorov M.D. Recovery of couplings and parameters of elements in networks of time-delay systems from time series. Phys. Rev. E, 2016, vol. 94, p. 052207.

42. Sysoev Ilya V., Ponomarenko Vladimir I. Reconstruction of the coupling matrix in the ensemble of identical neuron-like oscillators with time delay in coupling. Russian Journal of Nonlinear Dynamics, 2016, vol. 12, no. 4, pp. 567-576.

43. Sysoev I.V., Ponomarenko V.I., Pikovsky A. Reconstruction of coupling architecture of neural field networks from vector time series. Commun. Nonlinear Sci. Numer. Simulat., 2018, vol. 57, pp. 342-351.

44. Sysoev I.V., Ponomarenko V.I., Prokhorov M.D. Reconstruction of Ensembles of Oscillators with Nonlinear Time-Delay Feedbacks. Technical Physics Letters, 2018, vol. 44, no. 11, pp. 1024-1027.

45. Sysoev I.V., Ponomarenko V.I., Prokhorov M.D. Reconstruction of ensembles of nonlinear neurooscilla-tors with sigmoid coupling function. Nonlinear Dynamics, 2019, vol. 95, no. 3, pp. 2103-2116.

46. Sysoev I.V. Reconstruction of ensembles of generalized Van der Pol oscillators from vector time series. Physica D, 2018, vol. 384-385, no. 1, pp. 1-11.

47. Savitzky A., Golay M.J.E. Smoothing and differentiation of data by simplified least squares procedures. Analytical Chemistry, 1964, vol. 38, no. 8, pp. 1627-1639.

48. Nelder J.A., Mead R. A simplex for function minimization. Computer Journal, 1965, vol. 7, pp. 308-313.

49. Kendall M., Stuart A. The Advanced Theory of Statistics. New York: MacMillan, 1979.

50. Johnson N.L., Kotz S., Balakrishnan N. Continuous Univariate Distributions. New York: Wiley, 1995, vol. 2 , 752 p.

51. Ikeda K. Multiple-valued stationary state and its instability of the transmitted light by a ring cavity system. Optics Communications, 1979, vol. 30, no. 2, pp. 257-261.

52. Mackey M.C., Glass L. Oscillation and chaos in physiological control systems. Science, 1977, vol. 197, no. 4300, pp. 287-289.

53. Mandel I.D. Cluster Analysis. Moscow: Finansy i Statistika, 1988, 176 p. (in Russian).

54. Sompolinsky H., Crisanti A., Sommers H.E. Chaos in random neural networks. Phys. Rev. Lett., 1988, vol. 61, no. 3, pp. 259-262.

55. Levenberg K. A method for the solution of certain non-linear problems in least squares. Quarterly of Applied Mathematics, 1944, vol. 2, pp. 164-168.

56. Marquardt D. An algorithm for least-squares estimation of nonlinear parameters. SIAM Journal on Applied Mathematics, 1963, vol. 11, no. 2, pp. 431-441.

57. Coleman T.F., Li Y. An interior trust region approach for nonlinear minimization subject to bounds. SIAM J. Opt, 1996, vol. 6, pp. 418-445.

58. Bezruchko B.P., Smirnov D.A. Constructing nonautonomous differential equations from experimental time series. Phys. Rev. E, vol. 63, no. 1, p. 016207.

Сысоев Илья Вячеславович - родился в Саратове (1983), окончил факультет нелинейных процессов СГУ (2004), защитил диссертацию на соискание учёной степени кандидата физико-математических наук (2007). Доцент базовой кафедры динамического моделирования и биомедицинской инженерии, ответственный секретарь редакционной коллегии журнала «Известия вузов. ПНД». Научные интересы - исследование сигналов биологической природы методами нелинейной динамики, исследование эффективности и модернизация подходов к анализу сигналов. Автор более 40 публикаций.

Россия, 410012 Саратов, Астраханская, 83

Саратовский государственный университет имени Н.Г. Чернышевского Россия, 410019 Саратов, ул. Зеленая, 38

Саратовский филиал Института радиотехники и электроники РАН E-mail: ivssci@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