Вычислительные технологии
Том 18, № 4, 2013
Глобальная оптимизация на основе гибридного метода усреднения координат и метода роя частиц
В. Д. Кошур
Институт космических и информационных технологий, Сибирский федеральный университет, Красноярск, Россия e-mail: [email protected]
Представлен гибридный алгоритм глобальной оптимизации, в котором взаимно дополняющими компонентами являются метод усреднения координат и метод роя частиц. Для ряда тестовых функций показаны его преимущества и особенности работы.
Ключевые слова: глобальная оптимизация, метод усреднения координат, метод роя частиц.
Введение
Поиск решений нелинейных оптимизационных задач и особенно задач глобальной оптимизации (ГО) является одной из широко востребованных проблем вычислительной математики. В прикладных задачах целевая функция (ЦФ), как правило, имеет большое число переменных, не задана в аналитической форме и вычисляется как некоторая интегральная характеристика сложного динамического процесса (ЦФ-СДП), например, функция аэродинамического сопротивления проектируемого летательного аппарата.
Разработка эффективных методов, в определённой степени адаптивных к характеру изменяемости ЦФ, особенно актуальна в связи с развитием вычислительной техники и возможности использования параллельных вычислительных систем. Общие базовые положения по данной тематике приведены в работах [1-6], где рассмотрены оригинальные подходы и обзоры различных численных методов и их модификаций для решения задач оптимизации и ГО, а также в [7, 10], в которых представлены методы на основе неравномерных покрытий, реализованные как параллельные вычислительные алгоритмы ГО. Кроме того, необходимо отметить относительно новое направление исследований и построения алгоритмов ГО на основе методов интервального анализа, подробно изложенное в монографии [11].
В настоящей работе основное внимание уделено ряду методов "нулевого порядка", согласно которым путём вычислительного алгоритма определяются значения ЦФ только в пробных точках, с ориентацией в дальнейшем на решение задач с ЦФ-СДП при минимальных требованиях к непрерывности и ограниченности ЦФ. Для приближённых оценок изменяемости ЦФ используются значения максимумов величин, полученных как отношения разности значений ЦФ к расстоянию для всех пар пробных точек (нижняя оценка константы Липшица).
Наиболее востребованными в практических приложениях, по мнению автора, являются следующие детерминированные, стохастические и эвристические методы:
— методы на основе различных вариантов генетического алгоритма (ГА), эволюционных вычислений и их модификаций [12-18];
— методы на основе роя частиц — Particle Swarm Optimization (PSO) с введением адаптационных модификаций [19-22];
— методы случайного поиска и моделирования отжига [23-26];
— метод усреднения координат [5] и методы на основе инверсных регрессий [27-31];
— методы, основанные на оценках константы Липшица, неравномерных покрытий [6-10] и диагональные методы [32-36].
Перечисленные методы имеют как сильные, так и слабые стороны. В настоящей работе предлагается гибридный метод ГО на основе сочетания метода усреднения координат и PSO. Это позволяет в методе усреднения координат перейти от случайного выбора пробных точек к использованию текущих координат роя частиц, коллективное движение которых происходит адаптивно, подстраиваясь под характер изменения ЦФ, а при движении роя частиц учитывать их смещения в направлении найденного усреднённого центра на основе метода усреднения координат. Таким образом, осуществляется взаимодействие двух методов. Дополнительным приёмом, ускоряющим процесс сходимости гибридного алгоритма, является включение в процесс вычислений нескольких шагов процедуры Хука — Дживса, уточняющих текущие координаты лучшей и/или худшей частицы в рое.
Следует отметить, что предлагаемый гибридный метод, как и одна из его составляющих, PSO, не дает 100%-й гарантии нахождения глобального минимума ЦФ. Его результаты носят вероятностный характер. Данный метод пока не претендует на "промышленное" использование. Для этого необходимо проведение большого объёма вычислений, сравнительного анализа, статистической обработки результатов с использованием стандартизованных каталогов тестовых функций. Ниже описываются один из вариантов реализации гибридного метода и результаты проведённых вычислительных экспериментов для трёх специально выбранных тестовых функций. Эти численные результаты носят главным образом иллюстративный характер и позволяют в определённой степени оценить работоспособность метода.
1. Гибридный метод
Рассматривается ограниченная непрерывная функция f (х): П ^ К, где х = (х\,х2,... ,хп) Е П С Кп. Множество П является областью допустимых значений переменных и в простейшем случае представляет п-мерный параллелепипед с заданными
сторонами
I = 1, 2,... , п. Необходимо найти приближённое значение глобального минимума f * и хотя бы одну точку х*, в которой это значение достигается с заданной допустимой погрешностью £f для значений целевой функции:
^п = ттf (х), х Е П, f* - eí < ^ < f*, f* = f (х*). (1)
Вычислительная процедура нахождения приближённого значения координат точки х* в методе усреднения координат строится на основе итерационного процесса, который в непрерывной форме имеет вид [5]
xf+1] = J xtpjfe](x)dx, i = 1, 2,..., n, (2)
пи
pjk](x) = P?](/(x)) // Pjk](/(x))dx, (3)
пи
здесь k — номер шага вычислительного процесса, — область усреднения координат на шаге k. При этом вводится последовательность непрерывных функций Pj(y), s = 1, 2, 3,..., таких, что Vy G R значения Pj(y) > 0, и для последовательности вида Ps(y)/Ps(z) выполняется условие монотонного неограниченного возрастания при увеличении параметра селективности s и любых фиксированных значениях y, z с условием y < z. Примерами функций Pj(y), в частности, являются функции exp (—sy), s-y, y-j, а также класс функций вида (1 — yr )J при y G [0,1], r = 1, 2, 3,... Для приведённых ниже примеров численной минимизации использована функция (1 — y2)j.
По мере роста s крутизна ядер pjk](x) увеличивается, что в свою очередь приводит к увеличению веса координат, которые соответствуют лучшим значениям ЦФ, и в предельном случае последовательность усреднённых координат сходится к глобальному минимуму (соответствующая теорема сходимости доказана в работе [5]).
Для численной реализации метода усреднения координат одним из эффективных способов повышения точности вычисления интегралов (2) является последовательное увеличение количества пробных точек x(j)[fc], j = 1, 2,..., M[k], на k-м шаге итерационного процесса, т. е. M[k] > M[k-1]. Чтобы случайно не исключить точку глобального минимума, область усреднения Q[fc] в данном случае может рассматриваться как адаптивно изменяющейся [5], так и неизменной.
В предлагаемой модификации алгоритма итерационного вычисления усреднённых координат введено адаптивное смещение пробных точек, реализуемое как движение роя частиц в методе PSO с модификацией FDR (fitness-distance ratio based PSO) [22]. При этом вводится дополнительное смещение частиц роя к найденному центру усреднённых координат, что вносит в алгоритм PSO новый фактор обмена информацией между частицами и дополнительную стабилизацию процесса коллективного поиска роем частиц глобального минимума ЦФ.
При вычислении интегралов в формулах (2), (3) используется суммирование значений подынтегральных выражений на множестве пробных точек с учётом объёмов подобластей дискретизации области интегрирования Л В предлагаемом гибридном алгоритме координаты пробных точек отождествляются с координатами роя частиц, которые изменяются в процессе коллективного поиска глобального минимума и в начале вычислительного процесса задаются в узлах некоторой расчётной сетки или генерируются датчиком случайных чисел с равномерным распределением на интервалах i = 1, 2,... , n. В дискретной форме соотношения (2), (3) принимают
вид
M и
f+1] = £ xPW (xj)[fc]) Vj)[fc], i = 1, 2,..., n, j=i
pJfcl(gM(xWM))
PJJ x
m и
E Pjfc](g[fci(x(j)[fc]))V(j)[k] j=1
где V(j)[k] соответствует n-мерному объёму при разбиении области на подобласти, связанные с семейством точек интегрирования x(j)[fc], j = 1, 2,... , M[k]. Здесь g[k](x) — вспомогательные функции, масштабирующие значения целевой функции f (x(j)[fc]) к диапазону [0,1], которые определяются как
f (x(i)M 1 _ f [k]
gW(x«j,l4) = f f, (6)
fmax f min
где
fJSL = mm (f (x(1)[k]),f (x(2)[k]),... ,f (x(Mk)[k])), fmax = max (f (x(1)[k]), f (x(2)[k]),... , f (x(MW)[k])).
Предполагая, что для пробных точек — текущих координат роя частиц — в формулах (4), (5) при реализации приближённого интегрирования можно выбрать подобласти так, чтобы величины их объёмов V(j)[k], j = 1, 2,... , M[k], были приблизительно равными, расчётные формулы упрощаются к виду
M И
.,Jk+1| = £ jpW^W ), i =1,2,...,n, (7)
' (ДН) = (8)
р Vх ) м [к]
Е рР^КхСом))
3=1
Следует отметить, что для вычислительного алгоритма усреднения координат [5] конкретный вид подобластей дискретизации области интегрирования является не существенным, поэтому использование более простых при численной реализации соотношений (7), (8) вместо (4), (5) вполне оправданно.
Адаптивное изменение координат роя частиц при переходе на (к + 1)-й шаг проводится по схеме РБО [19-22], и дополнительно вводится составляющая движения к усреднённому центру координат х[к] для каждой ]-й частицы роя в следующей форме:
х(3)[к+1] = х(3)[к] + аУ(3)И + Б3)[к] + и[0, во] 0 (х[к] - х3)[к]). (9)
Здесь У(3)[к] — инерционная составляющая движения ]-й частицы, Б(3)[к] — вектор адаптивного смещения ^'-й частицы, определяемый тремя слагаемыми случайного смещения этой частицы [22]:
вшм = ар] + ар] + ар], (10)
где
3] = и[0,в1] 0 (х!з)[к]- х3)[к]),
аР] = и[0,в2] 0 (х^ - х(з)[к]), а(з)[к] = и[о, вз] 0 (х(^))[к] - х(з)[к]). (11)
В этих формулах использованы обозначения: х^7^ — лучшие координаты ]-й частицы
за к итераций, определяемые по значению ЦФ — когнитивная составляющая сме-
Л [к]
щения частицы); х^ — координаты лучшей частицы в рое с минимальным значением ЦФ за к итераций (в2?)[к] — социальная составляющая смещения частицы); х(9(:?'))[к] — координаты частицы с номером ), в направлении которой скорость убывания целевой функции наибольшая (в3?)[к] — составляющая изменяемости ЦФ по локальной оценке константы Липшица); и[0,в] — вектор с компонентами равномерно распределённых случайных чисел в интервале [0,в]; символом ® обозначено покомпонентное умножение векторов; коэффициенты а, вт, т = 0,1, 2, 3, являются настраиваемыми параметрами гибридного вычислительного алгоритма.
Таким образом, соотношения (4)—(11) после конкретизации вида ядер Р8[к](у) с возрастающим параметром селективности в и задания конкретных значений коэффициентов а, вт, т = 0,1, 2, 3, полностью определяют гибридный вычислительный алгоритм ГО на основе методов усреднения координат и роя частиц. Выбор коэффициентов численного алгоритма ГО может быть осуществлен путём мета-оптимизации [37] и выходит за рамки статьи.
2. Результаты вычислительных экспериментов
Для наглядного представления работы предлагаемого гибридного алгоритма ниже в графической форме приведены результаты вычислительного процесса поиска глобального минимума для трёх тестовых функций.
На рис. 1 показаны результаты процесса минимизации функции двух переменных вида
Рис. 1. Результаты процесса минимизации целевой функции (12) гибридным алгоритмом
/ (ж,у) = (ж - у)2 С08(ж)е0^ -У=]+2 (12)
в прямоугольной области [-12, 8] х [-15, 5] после 100 итераций (здесь точки и пунктирные линии — соответственно промежуточные положения и траектории 25 частиц в рое). Стартовые координаты частиц задавались в узлах регулярной сетки 5 х 5 расчётной области, значения коэффициентов а = 0, в0 = 0.1, в1 = 0, в2 = 0, в3 = 0.2; при этом достигнутое значение параметра селективности в = 212. Сплошной ломаной линией с маркерами показаны положения точек на поверхности функции для последовательно уточняемых значений усреднённых координат х[к]. Проведено 100 запусков программы, приемлемая точность (порядка 10-2) по значениям координат лучшей частицы достигалась за 10-15 итераций, в дальнейшем происходила последовательная концентрация частиц в окрестности глобального минимума. Следует отметить, что метод является статистическим и для получения "надежного результата по вероятности" необходимо проведение многократных запусков программы с измененными значениями случайных векторов и[0, в].
На рис. 2 представлены карта изолиний целевой функции (12) и результаты поиска глобального минимума при аддитивной случайной помехе с равномерным распределением, составляющей 20 % амплитуды изменения значений ЦФ (здесь пунктирные линии — траектории частиц, сплошная линия — изменение усреднённых координат). В расчёте использовано 36 частиц со стартовыми значениями в узлах регулярной сетки 6 х 6 расчётной области.
Рис. 2. Карта изолиний целевой функции (12) и результаты поиска глобального минимума при аддитивной случайной помехе
На рис. 3 приведён вариант вычислительного процесса минимизации классической тестовой функции Розенброка
/(х, у) = 100(у - X2)2 + (1 - х)2. (13)
Для данного расчёта было использовано 36 частиц в узлах сетки 6 х 6, параметры гибридного вычислительного алгоритма не изменялись. Следует отметить, что движение роя частиц в этом случае напоминает физический процесс постепенного скатывания капель воды в пологом желобе.
Переход к большим размерностям варьируемых переменных приводит к необходимости экспоненциального роста количества пробных точек или частиц в рое. Возможность проведения параллельных вычислений указывает на целесообразность использования в алгоритме нескольких семейств роев частиц. На рис. 4 показаны результаты итерационного процесса вычисления глобального минимума после 10 итераций для размерности п = 50 негладкой четырехэкстремальной функции вида
/(х) = —5 ехр (-351) - 10ехр(-2^2) - 7ехр(-2.5£з) - Зехр(-£4), (14)
где
п п п п
£ = ^ |х + 1|0'6 , S2 = ^2 |х»| , £3 = ^2 |Хг - 1|0'8 , £4 = ^2 |Хг - 2|°'9 . ¿=1 ¿=1 ¿=1 ¿=1
При этом были использованы 50 роев по 100 частиц в каждом, усреднение координат проводилось по всем частицам, значения координат лучшей частицы в каждом рое
Рис. 3. Процесс минимизации функции Розенброка
уточнялись процедурой Хука — Дживса с начальным шагом Л0 = 1.1 и числом внутренних итераций, равным пяти, параметры гибридного алгоритма приняты следующими: а = 0.4, во = 1.0, в1 = 0, в2 = 1.0, вз = 2.0.
Для данного варианта расчёта на рис. 5 представлены распределение значений целевой функции по частицам одного роя и значения координат лучшей частицы (с номером 37), которая достигла глобального минимума в точке х = 0.
Рис. 4. Минимизация функции (14) при п = 50: 1 — изменение максимальных значений ЦФ в списке частиц, 2 — изменение значений ЦФ в улучшаемых усреднённых координатах, 3 — изменение минимальных значений ЦФ в списке частиц
-2-
-4
Ш ж
в;;® йМв
J_1_
0 10 20 30 40 50 60 70 80 90 100
а
б
х 10"3
1 1 ОЙГ. 1 ни 1 1 1 1 Т а о- <?■ т # л г!г „.-i.fr Т П*
1 1 1 1 « 1 1 1 • | ш 1 1 1 1 А Г 1 I 4 4- J 1 1 1
0 10 20 30 40 50
Рис. 5. Распределение значений целевой функции по частицам одного роя (а) и значения 50 координат лучшей частицы, которая достигла глобального минимума функции (14) (б)
Рис. 6. Движение 36 частиц одного роя и изменение усреднённых координат в вычислительном процессе минимизации функции (14) для двух переменных
Проведённые вычислительные эксперименты по минимизации функции (14) при п = 100 показали, что если общее количество частиц не увеличивать, то гибридный алгоритм приводит к одному из локальных минимумов, чаще всего попадая в точку, где все координаты равны 2. Это объясняется тем, что данный локальный минимум имеет наиболее широкую область притяжения (последнее слагаемое в формуле (14)), в силу чего с большей вероятностью хотя бы одна из частиц попадает в указанную область и оказывает максимальное влияние на дальнейшее поведение роя частиц.
Для наглядности вид поверхности ЦФ (14) и вычислительный процесс минимизации для двух переменных показан на рис. 6. В этом вычислительном эксперименте было использовано 36 частиц одного роя с начальным положением частиц в узлах регулярной сетки 6 х 6.
Следует отметить, что использование нескольких семейств роев частиц позволяет в ряде случаев находить одновременно как глобальный минимум, так и локальные минимумы целевой функции, что может представлять интерес при решении прикладных задач.
Заключение
Представлены гибридный метод глобальной оптимизации на основе сочетания методов роя частиц и усреднения координат и его модификации с использованием нескольких роев частиц и включением процедуры Хука — Дживса. Для целевых функций с большим количеством переменных и высоким значением константы Липшица, как правило,
требуется эффективный выбор параметров алгоритма [37]. Предполагается дальнейшая разработка гибридного метода с нейросетевой подстройкой изменения параметров алгоритма, т. е. адаптивное управление этими параметрами, в частности, когда каждый из роев частиц осуществляет поиск самостоятельно и в то же время периодически обменивается информацией об успешном поиске с другими роями частиц по аналогии с компьютерной технологией мультиагентных систем [38]. Это направление исследований может стать определённой альтернативой экспоненциальному увеличению количества пробных точек при решении оптимизационных задач большой размерности.
Список литературы
[1] РАстригин Л.А. Статистические методы поиска. М.: Наука, 1968.
[2] Сухарев А.Г. Глобальный экстремум и методы его отыскания. М.: Изд. МГУ, 1981.
[3] Жиглявский А.А., ЖилинскАс А.Г. Методы поиска глобального экстремума. М.: Наука, 1991.
[4] Horst R., Pardalos P.M., et al. Handbook of Global Optimization. Dordrecht: Kluwer, 1995.
[5] Рубан А.И. Глобальная оптимизация методом усреднения координат. Красноярск: ИПЦ КГТУ, 2004. 302 с.
[6] Евтушенко Ю.Г. Методы решения экстремальных задач и их применение в системах оптимизации. М.: Наука, 1982. 432 с.
[7] Евтушенко Ю.Г., Малкова В.У., Станевичюс А.А. Параллельный поиск глобального экстремума функций многих переменных // Журн. вычисл. математики и матем. физики. 2009. Т. 49, № 2. С. 255-269.
[8] Евтушенко Ю.Г., Жадан В.Г. Об одном подходе реализации численных методов нелинейного программирования // Изв. АН СССР. Техническая кибернетика. 1983. № 1. С. 47-59.
[9] Евтушенко Ю.Г., Посыпкин М.А. Применение метода неравномерных покрытий для глобальной оптимизации частично целочисленных нелинейных задач // Журн. вычисл. математики и матем. физики. 2011. Т. 51, № 8. С. 1376-1389.
[10] Евтушенко Ю.Г. Численный метод поиска глобального экстремума функций (перебор на неравномерной сетке) // Там же. 1971. Т. 11, № 6. С. 1390-1403.
[11] Хансен Э., Уолстер Дж.У. Глобальная оптимизация с помощью методов интервального анализа. М.; Ижевск: НИЦ "Регулярная и хаотическая динамика", Институт компьютерных исследований, 2012. 516 с.
[12] Holland J.H. Adaptation in Natural and Artificial System. Ann Arbor: Univ. of Michigan Press, 1975.
[13] Goldberg D.E. Genetic Algorithms in Search, Optimization, and Machine Learning. Addison-Wesley, 1989.
[14] Курейчик В.М. Генетические алгоритмы. Состояние. Проблемы. Перспективы // Изв. РАН. Теория и системы управления. 1999. № 1. С. 144-160.
[15] Michalewicz Z. Genetic Algorithm + Data Structures = Evolution Programs. 3rd Edn. New York: Springer-Verlag, 1996.
[16] Stom R., Price K. Differential evolution — a simple and efficient heuristic for global optimization over continuous spaces // J. of Global Optimizat. 1997. Vol. 11. P. 341-359.
[17] Кошур В.Д., Ильин В.А. Нейронная сеть Хопфилда как кроссовер генетического алгоритма //V Всероссийская научно-техн. конф. "Нейроинформатика-2003". Сб. науч. тр. Ч. 1. М.: МИФИ, 2003. С. 92-100.
[18] Рутковскля Д., Пилинский М., Рутковский Л. Нейронные сети, генетические алгоритмы и нечеткие системы. М.: Горячая линия — Телеком, 2004.
[19] Kennedy J., Eberhard R. Particle swarm optimization // Proc. of IEEE Intern. Conf. on Neural Networks. 1995. Vol.4. P. 1942-1948.
[20] Mendes R., Kennedy J., Neves J. The fully informed particle swarm: Simpler, maybe better // IEEE Trans. on Evolutionary Comput. 2004. Vol. 8. P. 204-210.
[21] Kennedy J., Mendes R. Neighborhood topologies in fully informed and best-of-neigh-borhood particle swarms // Systems, Man, and Cybernetics. 2006. Vol. 36. P. 515-519.
[22] Карпенко А.П., Селиверстов Е.Ю. Обзор методов роя частиц (particle swarm optimization) // Электронное научно-техн. издание "Наука и образование". 2009. № 3. http://technomag.edu.ru/doc/116072.html.
[23] Aarts E.N.L., Laarhoven P.J.M. Simulated Annealing: Theory and Applications. London: Kluwer, 1987.
[24] Ali M., Storey C. Aspiration based simulated annealing algorithm // J. of Global Optimizat. 1996. No. 11. P. 181-191.
[25] Hamma B., Viitanen S., Torn A. Parallel continuous simulated annealing for global optimization // Optimization Methods and Software. 2000. Vol. 13, No. 2. P. 93-116.
[26] Орлянскля И.В. Современные подходы к построению методов глобальной оптимизации // Электронный журнал "Исследовано в России". 2002. С. 2097-2108. http://zhurnal.ape.relarn/ru/articles/2002/189.pdf.
[27] Кошур В.Д. Адаптивный алгоритм глобальной оптимизации на основе взвешенного усреднения координат и нечётко-нейронных сетей // Нейроинформатика. Электронный рецензируемый журнал. 2006. Т. 1, № 2. С. 106-124. http://www.niisi.ru/iont/ni/Journal/N2/Koshur.pdf.
[28] Кошур В.Д., Пушкарёв К.В. Глобальная оптимизация на основе инверсных соотношений и обобщённо регрессионных нейронных сетей // X Всероссийская научно-техн. конф. "Нейроинформатика-2008". Сб. науч. тр. Ч. 2. М.: МИФИ, 2008. С. 182-192.
[29] Koshur V., Kuzmin D., Legalov A., Pushkaryov K. Solution of large-scale problems of global optimization on the basis of parallel algorithms and cluster implementation of computing processes // Parallel Comput. Technolog. 10th Intern. Conf., PaCT 2009. Springer, 2009. P. 121-125.
[30] Кошур В.Д., Пушкарёв К.В. Дуальные обобщённо-регрессионные нейронные сети для решения задач глобальной оптимизации // XII Всероссийская научно-техн. конф. "Нейроинформатика-2010". Сб. науч. тр. Ч. 2. М.: НИЯУ МИФИ, 2010. С. 219-227.
[31] Кошур В.Д., Пушкарёв К.В. Глобальная оптимизация на основе нейросетевой аппроксимации инверсных зависимостей // XIII Всероссийская научно-техн. конф. "Нейроинформатика-2011". Сб. науч. тр. Ч. 1. М.: НИЯУ МИФИ, 2011. С. 89-98.
[32] Стронгин Р.Г. Поиск глобального минимума. М.: Знание, 1990.
[33] Нефёдов В.Н. Некоторые вопросы решения липшицевых задач глобальной оптимизации с использованием метода ветвей и границ // Журн. вычисл. математики и матем. физики. 1992. Т. 32, № 4. С. 512-529.
[34] Hansen P., Jaumard B. Lipschitz optimization // Handbook of Global Optimization / Eds. R. Horst, P.M. Pardalos. 1995. Vol. 1. P. 407-493.
[35] Sergeyev Y.D. On convergence of "Divide the Best" global optimization algorithms // Optimization. 1998. Vol. 44, No. 3. P. 303-325.
[36] Сергеев Я.Д., Квасов Д.Е. Диагональные методы глобальной оптимизации. М.: ФИЗ-МАТЛИТ, 2008.
[37] Pedersen M.E.H. Tuning & Simplifying Heuristical Optimization: Thesis for the degree of Doctor of Philosophy. Univ. of Southampton, School of Eng. Sci., Comput. Eng. and Design Group, 2010. 204 p.
[38] Shoham Y., Leyton-Brown K. Multiagent Systems. Algorithmic, Game-Theoretic, and Logical Foundations. Cambridge Univ. Press, 2009.
Поступила в 'редакцию 29 ноября 2011 г., с доработки — 10 июня 2013 г.