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

УДК 004.021, 681.3.069

ИДЕНТИФИКАЦИЯ ПАРАМЕТРИЧЕСКИ СВЯЗАННЫХ МОДЕЛЕЙ СЛОЖНЫХ СИСТЕМ С.В. Иванов

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

Введение

Развитие методов и технологий компьютерного моделирования стимулирует интерес исследователей к изучению так называемых сложных систем (complex systems). Под сложной системой подразумевается [1] система, которая (а) состоит из большого числа компонентов; (б) допускает «дальние» связи между компонентами; (в) обладает многомасштабной (в том числе пространственно-временной) изменчивостью. Перечисленные свойства порождают ряд специфических проблем компьютерной реализации моделей сложных систем. Так, большое число компонентов отражается на ресурсоем-кости вычислительных процедур и повышенных требованиях к объемам оперативной памяти для хранения структур данных. Многомасштабность сложных систем требует использования параметрически связанных моделей, описывающих иерархию смежных диапазонов изменчивости. Учет дальних связей, когда каждый компонент системы может взаимодействовать с каждым, приводит к необходимости предварительного решения проблемы связности при применении параллельных вычислений, что затрудняет или даже приводит к невозможности распараллеливания формальными методами. В силу условности выделения границ диапазонов изменчивости процедуры связывания или замыкания таких моделей традиционно используют эмпирические закономерности и коэффициенты, значения которых определяются посредством идентификации. Задача идентификации моделей заключается в определении их структуры и параметров по результатам наблюдений над входными и выходными переменными реальной системы и решается методами оптимизации [2].

Особенности идентификации моделей сложных систем

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

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

- неоднородность, которая возникает в силу различий в постановках экспериментов, методов и средств измерения.

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

В итоге данные измерений можно описать в виде набора векторов У = (У/ (у)) / Ё1к N, у ё1кМ, где Уу - разрозненные данные измерений по N источникам и М масштабам изменчивости, а V - пространственная и (или) временная координата. При этом особенности представления моделей связаны со способом представления их параметров. Так, вектор неизвестных (настраиваемых) параметров удобно представить как Е = (Е1,Е2,...,Е), где Е; - параметры /-ой модели, которые в общем

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

' УN = ^ (У, Е N X

yN-1 = FN-1(v, yN, E N-1)

(1)

. У1 = Fi(v, yN y2, Ei).

В итоге общая задача идентификации модели сложной системы состоит в минимизации невязки между выходами моделей разного уровня F1(E1),..., FN (E N) и разрозненными наблюдениями У за поведением реальной системы J(F1v..,FN,У) —min при заданном наборе ограничений G(Ei) еЦ,i е 1...N .

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

\

/ N

I Внешняя \ к

ч среда р?

E * е (E ^, E s )P

i \ i ? i SP

Рис. 1. Функциональная схема процедуры идентификации модели сложной системы

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

Идентификация параметров модели приводного ветра

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

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

к, ^)=

1

л р

дР дР ду' дх / параметр Кориолиса; р

(2)

где /к = 2О зт(ф) - параметр Кориолиса; р - плотность воздуха, О - угловая скорость вращения Земли; ф - широта места, и и vg - компоненты геострофического

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

Как правило, атмосферные течения не являются прямыми, а двигаются вдоль криволинейных траекторий. Это подразумевает дополнительные ускорения вдоль радиусов кривизны, т.е. дополнительную центростремительную силу, что порождает градиентный ветер. Модуль градиентного ветра можно вычислить по формуле

От =

л*

2

-1 +

1 +

40

л*

, где /к - параметр Кориолиса, О = ^и2§ + ^ - скорость гео-

строфического ветра, * - радиус кривизны изобары в интересующей точке.

Для перехода от скорости градиентного ветра От к скорости приводного ветра на

высоте 10 м Ол

о

10 необходимо знать коэффициенты в = —10

От

и а - угол отклонения на-

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

Идентификация параметров модели, по сути, разбивается на два этапа (рис. 2).

Рис. 2. Общая схема идентификации модели приводного ветра в синоптическом и климатическом диапазонах изменчивости

На первом этапе определяются параметры расчета прямоугольной сеточной области (5х, 5у). В качестве критерия качества идентификации предлагается использовать векторный коэффициент корреляции между геострофическим ветром и измерениями приводного ветра. Вторым этапом как в алгоритме моделирования, так и в процедуре идентификации является оценка коэффициентов перехода от градиентного ветра к приводному ветру, а также углов отклонения направления ветра от изобары в зависимости от стратификации атмосферы. При расчете этих коэффициентов необходимо, помимо значения модуля и направления смоделированного градиентного ветра, использовать значения температуры воды и воздуха в данной точке акватории (стратификация атмосферы). Коэффициенты перехода (а, в) определяются через средние многолетние значения отношения измеренного и смоделированного градиентного ветра в соответствующих диапазонах разностей температур воды и воздуха и силы градиентного ветра. Предложенный подход к идентификации параметров модели приводного ветра был использован для расчета ветра в акватории Атлантического океана. В качестве измерений были использованы данные с буя 44004, имеющего координаты 38,47° (с.ш.) и 70,56° (з.д.).

25

1 20 39 58 77 96 115 134 153 172 191 210 229 248 267 286 305 324 343 362 381 400 419 438 457 476 495

часы [2005.3.24-2005.3.25 ]

Модельный ветер .....................Измерения

Рис. 3. Фрагменты реализаций моделей ветра для акватории Атлантического океана

Уточнение параметров модели посредством идентификации позволило повысить качество результирующих данных (векторный коэффициент корреляции составляет 0,82 по отношению к измерениям; для данных реанализа КСБР/КСЛЯ эта величина составляет 0,74).

Идентификация параметров модели климатических спектров

Под климатическим спектром волн понимают характерный элемент ансамбля спектров, имеющий определенную вероятность и обусловливающий некоторые характерные волновые условия на заданной акватории. Функционально подобные классы спектров могут состоять из нескольких волновых систем и представлять собой хорошо известные аппроксимации, что позволяет представить спектр £ (ш, 9) в виде £ (ш, 9,2), где 2 - набор параметров. В данной работе рассматривается аппроксимация спектра вида JONSWAP [3], классическая запись которого в частотной области ш имеет вид

£ (ш) = а^рехр

ш

-1,25( ш)-

ш р

у5, 5 = ехр

(ш-ш Р )2 2а2шР

(3)

Здесь а - так называемый параметр Филипса (параметр масштаба спектра), у -параметр пиковатости, а - параметр формы, ш р - частота максимума спектра.

4

Даже при постоянном и устойчивом ветре волны распространяются не в одном фиксированном направлении, а в некотором секторе. Частотно-направленные спектры ветрового волнения и зыби можно представить в форме 5(ш, 0) = 5(0), где 5(ш) -

частотный спектр волнения, а Q (0) - функция углового (0) распределения энергии

_ |- /2 -,-1

* 1008(0) 5 00

Q(0) = С(*)[еоз(-—)]*, С(*) =

2

-п/2

(4)

Поскольку объектом анализа является спектральная плотность случайного процесса волнения (как квадратичная характеристика энергии волн), критерием качества для процедуры идентификации является сумма квадратов невязки модельного спектра

и измерений 5(ш, 0). Таким образом, расчетная процедура сводится к нелинейной оптимизации функционала

<-*-) 71 N

3(31 кЕN) ={{ 5(ш,0)-Ё5,(ш,0,Е,)

dшd0

(5)

0 0 1_ I=1

относительно набора параметров Е1 N (N - неизвестное число волновых систем).

В результате идентификации параметров модели спектра выполняется снижение мерности массива данных - сведение 5(ш, 0, х, у, ^) к набору параметров Е(х, у, ^) существенно меньшей мерности, чем исходные данные, что позволяет производить классификацию спектров относительно полученных параметров.

На рис. 4 приведена схема идентификации модели климатических спектров морского волнения. Ее отличием от аналогичной схемы для расчета полей приводного ветра (рис. 2) является применение процедуры идентификации к наблюдениям синоптического диапазона для получения характеристик в климатическом диапазоне, в то время как для модели расчета полей ветра наблюдения как раз климатического диапазона использовались для настройки модели в синоптическом диапазоне.

База данных спектров

Рис. 4. Общая схема идентификации модели климатических спектров

Модель (3) различает одно-, двух- и многопиковые (по переменным ш и 0) спектры. Выделены пять классов климатических спектров: ветровые волны; зыбь (к=2); ветровые волны и близкая по частоте «свежая» зыбь (к=3); ветровые волны и «свежая» зыбь, различающиеся и по частоте, и по направлению (=4); ветровые волны и зыбь, близкие (плохо различающиеся) по частоте и направлению (к=5). Каждый класс соответствует устойчивому состоянию к, следовательно, синоптическая изменчивость вол-

2

нения может быть представлена как марковская цепь к = к(^) с матрицей переходных вероятностей р(д+1) = р{к+1) = г | к) = 7} г,7 = 1, т , и вектором предельной вероятности п7 = р{к^) = 7} 7 = 1, т . Оценка коэффициентов в матрице переходных вероятностей и векторе предельной вероятности также является частью процедуры идентификации параметров модели.

В табл. 1 показана перемежаемость климатических спектров в течение года для центральной точки Тирренского моря.

Таблица 1. Вероятностные характеристики (%) перемежаемости климатических спектров морского волнения. ВЕСЬ ГОД

Класс Переходные вероятности из класса в класс (%) Повторяемость по классам (%)

I II III IV V

I 61 - 36 - 3 11

II - 66 30 3 1 3

III 8 2 75 1 14 46

IV 1 - 23 51 25 4

V 1 - 14 4 81 36

Идентификация параметров эпидемиологической модели ВИЧ

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

Эпидемиологическая динамика ВИЧ описывается моделью комплексной сети, в которой каждый узел представляет собой отдельного человека, находящегося в одном из трех основных состояний: восприимчивый к вирусу (£), зараженный (I), удаленный из сети в результате диагностики СПИД или смерти (Я). Сетевые связи (ребра графа) указывают на наличие контактов между людьми, которые могут приводить к новым случаям заражения. Для корректного описания динамики вируса внутри отдельного человека между состояниями зараженный и удаленный вводится ряд дополнительных промежуточных состояний, описываемых нестационарной марковской моделью специального вида [4]. Нестационарность модели определяется наличием факторов, оказывающих влияние на вероятность заражения, эффектом диагностики и лечения.

Комплексная сеть задается набором {О, З}, где О - ориентированный граф с множествами вершин V и ребер Е, а З - эволюционный оператор, определяющий изменение сети за интервал времени ^: {V, Е^ + =ЗV, Е^. Эволюционный оператор

может быть представлен как композиция нескольких операторов З = ® Зк, соответст-

к

вующих различным динамическим факторам (процессу распространения инфекции, динамике ребер и динамике узлов). В общем случае такие операторы являются некоммутативными.

Процедура моделирования эпидемиологической динамики ВИЧ задается следующей последовательностью шагов.

1. Сгенерировать сеть в соответствии с выбранным распределением степеней вершин и начальным числом инфицированных р0.

2. Заразить узлы, имеющие общие ребра с инфицированными узлами, с вероятностью X для каждого узла независимо.

3. Для каждого инфицированного узла применить правило развития СПИД через коэффициенты переходных вероятностей р j. Узлы, достигшие стадии СПИД, из сети

удаляются.

4. Применить демографическое правило (часть узлов удаляется из сети, а часть добавляется в соответствии с демографическим балансом).

5. Зафиксировать текущее состояние узлов и сгенерировать новую сеть.

6. Повторять (2)-(5) необходимое число раз.

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

коэффициент d. Это позволяет сформулировать задачу идентификации модели комплексной сети как задачу оптимизации, на каждом шаге которой выполняются шаги алгоритма (2)-(5),

2007

2 [r (t) - AIDS, ]2 ^ min. (6)

t=1981 (X, d,p0)

Здесь AIDSt - ежегодное число официально регистрируемых случаев СПИД в

конкретной стране; в силу специфики AIDSt как интегральной характеристики по

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

Начальные параметры

{p(0), d, XX

Идентификация

p(0), d, {X}, Pj

Рис. 6. Общая схема идентификации параметров модели эпидемиологической динамики ВИЧ

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

ны данные по случаям СПИД в США (доступны с 1981 г.). На рис. 7, а-б, приведены интегральные характеристики эпидемиологической ситуации - количество заболевших СПИД и количество ВИЧ-инфицированных.

а б

Рис. 7. Результат моделирования динамики ВИЧ в США: а - результат моделирования и данные измерений пол случаям СПИД; б - модельная кривая случаев ВИЧ (для гомосексуальной популяции)

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

Выводы

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

Литература

1. Boccara N. Modeling complex systems. - Springer, 2004. - 397 p.

2. Эйкхофф П. Основы идентификации систем управления. - М.: Мир, 1975.

3. Hasselmann K., et al. Measurements of wind-wave growth and swell decay during the Joint North Sea Wave Project (JONSWAP) // Dtsch. Hydrogr. Z. - 1973. - Р. 1-95.

4. Sloot P.M.A., Ivanov S.V., Boukhanovsky A.V., Van De Vijver D.A.M.C. and Boucher C.A.B. Stochastic simulation of HIV population dynamics through complex network modeling // International Journal of Computer Mathematics. - 2008. - Р. 1175-1187.

Иванов Сергей Владимирович - Санкт-Петербургский государственный универси-

тет информационных технологий, механики и оптики, м.н.с, svivanov@mail.ifmo.ru