CC BY
23
структуры и моделирование 2017. №2 1(41). С. 55-74
УДК 519.16
О ДОВЕРИТЕЛЬНЫХ ИНТЕРВАЛАХ ДЛЯ ЧИСЛА ЛОКАЛЬНЫХ ОПТИМУМОВ
А.В. Еремеев1
доцент, д.ф.-м.н., e-mail: eremeev@onm.oscsbras.ru К.Р. Ривс2
профессор, e-mail: c.reeves@coventry.ac.uk
1 Институт математики им. С.Л.Соболева Сибирского отделения РАН
2Университет Ковентри
Аннотация. Число локальных оптимумов является важным показателем сложности задач оптимизации для алгоритмов локального поиска. В настоящей работе обсуждаются некоторые методы нахождения доверительных интервалов для этого параметра в задачах, где мощность пространства поиска не позволяет выполнить полный перебор решений. Представлены результаты вычислительных экспериментов, полученные с использованием предложенных методов для модели NK ландшафтов С. Кауффмана, задачи отыскания двоичной последовательности с минимальной автокорреляцией, задачи о вершинном покрытии и задачи распределения бункерных ёмкостей в производственных линиях.
Ключевые слова: локальный поиск, задача комбинаторной оптимизации, перепись Шнабеля, консервативный доверительный интервал.
Введение
Значительное число эвристических методов решения задач комбинаторной оптимизации основано на поиске улучшающих решений в окрестности некоторой текущей точки. При этом могут быть использованы различные стратегии выбора новой текущей точки из окрестности, в которую будет сделан очередной переход. На практике такие методы быстро приходят к локальному оптимуму, после чего требуется либо начинать поиск с новой точки (см., например, [20,30,34]), либо прибегать к некоторому правилу постепенного выхода из окрестности, как это делается в алгоритмах поиска с запретами [14] и поиска с переменными окрестностями [16]. Как показывает практика, для успешного применения таких методов большое значение имеет удачный выбор системы окрестностей. Сложность индивидуальной задачи с заданной системой окрестностей характеризуется такими параметрами, как трудоёмкость выбора улучшающего решения в окрестности, среднее и максимальное число шагов до локального оптимума, гарантированная относительная погрешность локальных оптимумов и общее их число. Исследование генетических алгоритмов также
приводит к необходимости анализа структуры множества локальных оптимумов [26].
В настоящей работе предлагаются и обосновываются статистические методы построения доверительных интервалов для числа локальных оптимумов на основе результатов многократного выполнения локального поиска. Методы предназначены для исследования индивидуальных задач комбинаторной оптимизации, в которых большая размерность пространства решений не позволяет осуществить полный перебор всех его элементов за приемлемое время. Будем далее обозначать число локальных оптимумов задачи через V.
Из литературы известно несколько подходов к оцениванию величины V на основе результатов многократного выполнения локального поиска из случайно выбранных начальных точек. Один из подходов основывается на упрощающем предположении о том, что все локальные оптимумы имеют равные вероятности попадания в них (предположение изотропности) [27]. Другие подходы состоят в подборе параметров некоторого распределения вероятностей (экспоненциального, гамма, логнормального и т.д.) [12,28] для описания вероятности попадания в различные локальные оптимумы или использовании таких непараметрических методов, как бутстреп-метод или метод складного ножа [28,29]. Некоторые из этих подходов основываются на методах, разработанных в статистике и экологии для оценивания численности популяций (см., например, [24,33]).
Настоящая работа представляет собой расширение краткой статьи [11] и содержит более детальный теоретический анализ и более широкое экспериментальное исследование. В п. 1 рассматриваются методы построения консервативных доверительных интервалов, основанные непосредственно на функции распределения наблюдаемой целочисленной случайной величины. Здесь же предложены методы их вычисления, основанные на первом повторении локального оптимума и на числе различных локальных оптимумов, полученных за время испытаний. В п. 2 приводятся описание тестовых примеров и результаты вычислительных экспериментов, проведённых с предложенными методами построения доверительных интервалов и с непараметрическим бутстреп-методом из [29]. Описание последнего содержится в приложении.
1. Методы построения доверительных интервалов
Рассмотрим сначала достаточно общую ситуацию, когда при любом значении представляющего интерес неизвестного параметра V функция распределения наблюдаемой целочисленной случайной величины Т известна. Пусть ¥(Ь; V) = Р{Т ^ ¿} — функция распределения величины Т при указанном значении V. Наблюдаемое значение величины Т, полученное в фактически проведённом эксперименте, будем обозначать через т, т € N.
Если упрощённо считать, что V может принимать любое вещественное значение и ¥(Ь; V) непрерывно зависит от V, то при заданных вероятностях аьа2, таких, что а + а2 < 1, доверительный интервал ^\(т),^(т)] уровня
(1 — а1 — а2) ■ 100% может быть построен решением уравнений
1 — Е(т; VI) = «1, Е(т; = «2 (1)
относительно неизвестных v1 и ^ (см., например, [4,22]).
Далее интервал v1(т) ^ V ^ ^(т) будем называть консервативным доверительным интервалом с доверительным уровнем (1 — а) ■ 100%, если с вероятностью не менее 1 — а выполняются неравенства v1(т) ^ V ^ v2(r). В случае ^(т) = консервативный доверительный интервал будем называть консервативной нижней границей.
Поскольку в поставленной задаче параметр V принимает только натуральные значения, необходимо соответственно адаптировать (1).
Утверждение 1. Пусть
V!(т) = € N : 1 — Е(т — 1; V7) ^ а1>,
v2(т) = шax{v7 € N : Е(т; V7) ^ а2},
тогда [v1(т),v2(т)] является консервативным доверительным интервалом целочисленного параметра V с уровнем значимости (1 — а1 — а2) ■ 100%, причём P{v < v1 (Т)} ^ аь Р{^(Т) < V} ^ а2.
Доказательство. Пусть для некоторого значения V утверждение неверно. Тогда имеет место по крайней мере одно из двух неравенств: а1 < Р^ < v1(т)}, а2 < Р^(т) < V}. В первом случае
а1 < Р^ < шш( V7 : 1 — Е(Т — 1; V7) ^ а1)} ^ Р{1 — Е(Т — 1; V) < а1} =
= Р{Т — 1 ^ ш1п(0 : 1 — Е(0; V) < а1)} = Р{Т ^ а1)},
где 0т1п(V, а1) = шт{07 : 1 — Е(07 — 1; V) < а1}. Однако
Р{Т ^ 0тт(V, а1)} = 1 — Е(0min(v, а1) — 1; V) < а1 (2)
по определению 0min, что приводит к противоречию. Во втором случае
а2 < P{шax{v7 : Е(Т; V7) ^ а2} < V} ^
^ Р{Е(Т; V) < а2} = Р{Т < ^п^)},
где $тт(V, а2) = шт{0 : Е(0; V) ^ а2}. Однако Р{Т < б^т^, а2)} < а2. □
Модифицируем предыдущее утверждение с целью нахождения консервативной нижней границы заданного доверительного уровня а в случае, когда распределение вероятностей случайной величины Т зависит не только от V, но и от дополнительного вектора параметров у, принадлежащего некоторому множеству У(V) С . Обозначим в таком случае функцию распределения для Т через Е(£; V, у).
Утверждение 2. Пусть
VI (т )=шт^ V е N : 1- т£ Е(т - 1; V', у) ^ а>, \ уел*') /'
тогда P{v < VI(Т)} ^ а.
Доказательство. Достаточно в доказательстве предыдущего утверждения заменить определение 0тщ^, а) на
©тт^,а) = тт|©' € N : 1 - у т£ ^(©' - 1; V,у) < а| ,
и соответствующим образом модифицировать неравенство (2): Р{Т ^ 0min(v,а)} = 1 - Е(0т1п(^а) - 1; V, у) ^ ^ 1 - т£{F(0min(v,а) - 1; V, у) : у €?(V)} < а.
□
Если при любом Ь функция F(t; V) (функция Fy(Ь, V) = т£уеуF(т --1; V', у)) монотонна по V, то доверительный интервал из утверждения 1 (утверждения 2) может быть вычислен с помощью метода дихотомии. Ниже будет показано, что для рассматриваемых в настоящей работе оценок величины V монотонность по V имеет место.
1.1. Ожидание первого повторения
Пусть х^ — вероятность попадания локального поиска в локальный оптимум г, г = 1,..., V, и пусть Т — номер испытания локального поиска, при котором впервые получен ранее встречавшийся локальный оптимум. Функцию распределения случайной величины Т, соответствующую вектору р = (р1,...,р^), будем обозначать через Fp(t; V). В частном случае изотропного распределения, когда вероятности попадания во все локальные оптимумы одинаковы, введём обозначение р = (1^,..., 1^).
Следующее утверждение показывает, что изотропное распределение является в некотором смысле экстремальным.
Утверждение 3. Для любых Ь, V € N Ь ^ 2, функция Fp(t; V) достигает своего единственного минимума при р = р.
Доказательство. Пусть 5 — множество всех последовательностей длины Ь из элементов множества {1,...^}, не содержащих повторений, и пусть а(г) — г-й элемент последовательности а. Тогда
г
Fp(t; V) = 1 - Р{Т > Ь} = 1 - ^ Д р^).
аеяi=1
Fp(t; v) — непрерывная функция на компактном множестве
л = j(pi,P2,...,Pv) ^ 0 : ^pi = 1
когда t фиксировано. Таким образом, существует распределение вероятностей попадания в локальные оптимумы, при котором достигается минимум функции Fp(t; v) на Л. Покажем, что никакое распределение p, кроме p, таким свойством не обладает.
Предположим противное: пусть минимум достигается на векторе p, где pi = pj для некоторых i и j. Не теряя общности, предположим, что i =1, j = 2, p2 = Yp1 и y < 1. Рассмотрим распределение вероятностей
, , pi + p2 , ,
pi = p2 = -2—, P3 = p3, ... ,Pv = pv
и соответствующую ему функцию распределения Fp(t; v). Достаточно показать, что Fp/ (t; v) < Fp(t; v), т.к. это приведёт к противоречию с оптимальностью р. Обозначим через S(j) подмножество S, состоящее из тех последовательностей, в которые входит элемент j, j = 1,2. Тогда Fp(t; v) = = 1 - а(р) - b(p) - c(p) - d(p), где
t t
a(p) = Y1 П^ф b(p) = Y1 E[pff(i)'
CT€S(1)nS(2) i=1 ct€S(1)\S(2) i=1
tt c(p) = Y1 E[pff(i)' d(p) = Y1 IW
ct€S(2)\S(1) i=1 ct€S\(S(1)US(2)) i=1
Аналогично Fp/(t; v) = 1 - a(p') - b(p') - c(p') - d(p') и d(p') = d(p). Заметим, что если элемент 2 во всех последовательностях из S(2)\S(1) заменить на 1, то будет получено множество последовательностей S(1)\S(2). Таким образом, с одной стороны,
tt b(p)+c(p) = Y^ np^(i)+y ( ^ np^(i)
CT€S(1)\S(2) i=1 \o-eS(1)\S(2) i=1
= (1+y)p1 • Y1 П p^(i).
ct€S(1)\S(2) i:o-(i)=1
С другой стороны, p1 = p2 = (p1 + p2)/2, поэтому
b(p7) + C(p;) = (p1 + p2) П p^(i)
ct€S(1)\S(2) i:a(i) = 1
и b(p) + c(p) = b(p') + c(p'). Далее
t
Kp) = YP2 • Y1 П P^
CTGS(1)nS(2) i:<r(i)=1,2
2 г
^ П Р*)-
Однако (1 + 7)2 > 47, поэтому а(р') > а(р) и Ер/(¿; V) < Ер(£; V). □
Как следствие из доказанного утверждения вытекает свойство монотонности функции Ер(£; V).
Следствие 1. Ер(£; V) является строго убывающей функцией от V.
Доказательство. Достаточно сравнить оптимальное по утверждению 3 значение Ер(£; V) с величиной Ер(£; V), где р =(1/^ - 1), 1/(v - 1),..., 1/(v - 1), 0). □
Ввиду монотонности Ер(£; V) по V, в изотропном случае доверительный интервал для V может быть вычислен методом дихотомии с использованием утверждения 1, при этом Е(£; V) = 1 — Консервативная нижняя гра-
ница v1 с заданным доверительным уровнем 100 ■ (1 — а)%, справедливая при любом распределении р, вытекает из утверждения 2. Далее будем обозначать эту границу через 1 — а).
Если задача оптимизации имеет чрезвычайно большое число локальных оптимумов, а число испытаний локального поиска г оказалось недостаточным, чтобы встретить хотя бы один локальный оптимум более одного раза, то консервативная нижняя граница доверительного уровня 100(1 — а)% для V может быть вычислена в предположении о том, что следующее (гипотетическое) испытание приведёт к первому повторению. Оценим скорость роста оценки Ь(г + 1,1 — а) как функции от г. Заметим, что согласно утверждению 2, Ь(г + 1,1 — а) есть минимальное значение при котором г!/^г ^ а. Однако
(г)^ = (^ — г)!^ = 6ХР{^— ^
г — Н ( г(г + 1)1 Г г2 1
= --^ — 2^1 ,
поэтому Ь(г + 1,1 — а) ^ 0.5г2/ 1п(а-1), и, к примеру, при а = 0.05 имеем 0.5г2/ 1п(а-1) ~ 0.167г2. При больших значениях г эта простая оценка хорошо аппроксимирует величину Ь(г + 1,1 — а).
1.2. Метод переписи Шнабеля
Один из методов, используемых в экологии для оценки численности популяций, носит название переписи Шнабеля [32,33]. Из популяции с неизменным составом производится отлов заранее заданного числа особей, каждая из которых имеет константную вероятность отлова. При отлове особь помечается (если она не была помечена ранее) и возвращается обратно в популяцию. Статистические оценки численности популяции вычисляются на основе данных о числе помеченных особей среди отловленных.
а(р')
Р1(1 + 7)
2
Для оценки числа локальных оптимумов метод переписи Шнабеля адаптируется следующим образом. Выполняются г испытаний локального поиска из начальных точек, выбранных случайным образом с одинаковым распределением (число г выбирается априори). Пусть К — случайная величина, равная числу различных локальных оптимумов, полученных за г испытаний. Обозначим через Ер (к; V, г) = Р{К ^ к} функцию распределения величины К в предположении о том, что р = ... /ри) — вектор вероятностей попадания в локальные оптимумы.
Для вычисления доверительных интервалов на параметр V известны приближенные методы, в которых используется аппроксимация распределения величины К нормальным распределением (см., например, [27,33]), а также точные методы [13,23], основанные на непосредственном вычислении функции распределения, подобно выражениям из утверждения 1.
Как было замечено в [7], в изотропном случае распределение вероятностей случайной величины К имеет вид:
Р{К = к} = -- ^ — к + 1) *(г-к) = V 5(^ , (3)
V7- (V — к)! V7
где 5(г, к) — число Стирлинга второго рода, равное количеству способов разбить г-элементное множество на к непустых подмножеств без учета порядка подмножеств, 5(г, к) = £ £*=о(—1)^ (*)37.
Для вычисления оценки максимального правдоподобия г>(г, к) может быть рассмотрена величина
г, к) = 1п(Р{К = к}) = 1п V! — 1п^ — к)! + 1п 5(г, к) — г 1п V.
В [8] Дж. Даррохом предложено приближённо (с погрешностью не более единицы) вычислять г>(г, к), как корень уравнения г, к) = — 1,г, к), т.е.
(V — 1)г = vr-1(v — к), (4)
т.к. 1пV! — 1п^ — 1)! = 1пV. Корень уравнения (4) может быть вычислен с помощью метода Ньютона или методами оптимизации унимодальной функции на прямой. Заметим, что при к = г уравнение г, к) = 0 не имеет конечного
решения, и оценка £(г, к) лишена смысла.
Найденные в [8] асимптотические оценки дисперсии позволяют построить консервативный доверительный интервал с центром в точке £(г, к), используя упрощающее предположение о нормальности распределения данной оценки. Как показано в [13,23], более точно вычислить консервативный доверительный интервал удаётся при использовании нормальной аппроксимации распределения К. Однако оба подхода могут давать существенные ошибки в ситуациях, когда предположение нормальности оказывается неуместным. Как видно из рис. 1, доверительный интервал из [8] даёт корректные оценки при достаточно больших значениях г. При малых значениях г возникают погрешности из-за нарушения гипотезы нормальности.
12000 10000 8000 6000 4000 2000
/ _ / V / / /
/ / / / к/г = 0.95
/- / / / / / / / /
\ / к/г = 0.8 х ----
у У У
____________.!______________________" kr = 0.7
500
1000
1500
Г
2000
Рис. 1. Доверительные интервалы для V, построенные симметрично к оценке т)(т,к) по методу Дж. Дарроха [8] и по точным формулам из утверждения 1. Пунктирные линии показывают 95% доверительные интервалы, полученные в предположении нормальности К. Сплошные линии — оценки максимального правдоподобия с точными 95% доверительными интервалами.
Рассмотрим следствие основного результата А.М. Зубкова и Н.Н. Попова из [2], на основе которого будет получен аналог утверждения 3 для переписи Шнабеля. Пусть имеется r частиц, размещаемых случайным образом, независимо, с одинаковым распределением по ячейкам 1, 2,...,v. Вероятность попадания в j-ю ячейку равна pj, j = 1,...,v, Yl'j=1 Pj = 1. Как показано в [2], в случае изотропного распределения вероятностей (p1,... ,pv) = p число занятых ячеек ^(r, p) «стохастически больше» числа занятых ячеек ^(r, p') при любом другом распределении вероятностей p', т.е. P{^(r,p) ^ k} ^ P{^(r,p') ^ k} для всех r G N, k = 1,..., v.
Обозначая через Fp(k; r, v) функцию распределения случайной величины K при вероятностях попадания в локальные оптимумы p1,... ,pv, p = (p1,... ,pv), указанный выше результат из [2] записывается как
Утверждение 4. При любых фиксированных v G N, r G N, k = 0,..., r, функция Fp(k; r, v) достигает минимума в точке p = p.
Утверждение 4 позволяет вычислять консервативную нижнюю границу L'(r, k, p) для v с заданным доверительным уровнем p ввиду утверждения 2. Корректность такой нижней границы при условии справедливости утверждения 4 была отмечена нами в [11]. Аналогичная нижняя граница предложена Ч. Мао в [21]. Заметим, что Fp(r; v, r) = 1 при любом v, поэтому при k = r нижняя оценка L'(r, k, p) не определена. Подобно доказательству следствия 1, из утверждения 4 получаем
Следствие 2. Fp(k; r, v) является невозрастающей функцией от v при любых r G N, k = 0, . . . , r.
Последнее следствие было независимо доказано в [11,23] без использования результата [2]. Это следствие даёт основание применению метода дихотомии для построения нижней границы ¿'(г, к, р) в соответствии с утверждением 2, а в изотропном случае - консервативного доверительного интервала в соответствии с утверждением 1.
2. Экспериментальное исследование 2.1. Используемые задачи
Изложенные выше методы оценивания числа локальных оптимумов были опробованы на ряде комбинаторных задач оптимизации. Среди них задачи с целевыми функциями, построенными случайным образом по модели ЛК-ландшафтов С. Кауффмана [19], задачи отыскания двоичной последовательности с минимальной автокорреляцией [15], задачи о вершинном покрытии и задачи распределения бункерных ёмкостей в производственных линиях [9]. Перечисленные задачи формулируются следующим образом.
Модель ЛК-ландшафтов предложена в [19] для генерации индивидуальных задач безусловной оптимизации с настраиваемой «изрезанностью ландшафта» целевой функции ^ : {0,1}п ^ при размерности пространства решений п = N. Значение этой функции в точке x е {0,1}п определяется как среднее арифметическое значений вспомогательных функций ^¿(ж^ж^, ...,ж^), г = 1,...,п, каждая из которых зависит от соответствующей компоненты ж^ а также от других К компонент ж^ , ...,ж^. Наборы номеров {г1, ...,гк} на этапе генерации исходных данных выбираются случайным образом среди К-элементных подмножеств множества {1,..., г — 1, г + 1,..., п} с равномерным распределением. Значения функции г = 1,...,п, для каждого из 2К+1 возможных значений её аргумента назначаются датчиком случайных чисел, выбираемых из интервала [0,1). Эти значения также входят в исходные данные задачи.
Задача отыскания двоичной последовательности с минимальной автокорреляцией состоит в отыскании вектора у е {—1, +1}п, минимизирующего
Его— 1 2
с2, где с, — величина автокорреляции для последовательности у при смещении ], которая в данном случае имеет вид с, = ^¿+.7. Ин-
дивидуальная задача отыскания двоичной последовательности с минимальной автокорреляцией полностью описывается числом п.
Невзвешенная задача о вершинном покрытии (ЗВП). Данная задача является ЫР-трудной в сильном смысле (см., например, [1]) и имеет следующую постановку. Пусть дан граф С = (V, Е) с множеством вершин V и множеством рёбер Е. Подмножество С с V называется вершинным покрытием в С, если каждое ребро инцидентно хотя бы одной вершине из С. Требуется найти вершинное покрытие С* минимальной мощности.
Пусть в графе С = (V,Е) множество вершин имеет вид V = ,...,гп}. Для представления ЗВП в форме задачи безусловной оптимизации на множестве {0,1}п предположим, что элемент ж^-, ] = 1,...,п, равен 1 тогда и только тогда, когда вершина г входит в искомое множество. Обозначим через Смножество вершин, соответствующее вектору x. Пусть целевая функция /= £™=1 ж^ в случае, если С— вершинное покрытие; иначе вводится штраф за нарушение ограничений задачи и значение целевой функции полагается равным /= £™=1 ж^ + и^) ■ IV|, где и^) — число рёбер в Е, которым не инцидентна ни одна вершина из С
Для экспериментов с ЗВП далее используются тестовые примеры для задачи о наибольшей клике из коллекции ШМАСБ [18], преобразованные в индивидуальные задачи о вершинном покрытии с помощью известного преобразования (см., например, [1]).
Задача распределения бункерных ёмкостей в производственной линии.
При управлении производственными линиями, в которых детали перемещаются от одной единицы оборудования (ЕО) к другой с помощью некоторого транспортного механизма, возникает следующая задача оптимизации объёмов бункеров. Предполагается, что вследствие отказов оборудования в процессе работы линии возникают остановки ЕО, случайные по моменту возникновения и длительности. Последствия отказов распространяются на смежные операции из-за невозможности передать деталь на следующую операцию или отсутствия деталей на входе ЕО. Наличие бункеров (ёмкостей) для складирования деталей между ЕО позволяет снизить влияние отказов на соседние операции и повысить среднюю производительность линии. Однако установка бункеров связана с дополнительными капитальными затратами и увеличением числа складируемых деталей. Задача состоит в выборе объёмов бункеров с учётом средней производительности линии V(Ь) в стационарном режиме, капитальных затрат на установку бункеров /(Ь) и средней стоимости хранения деталей ^(Ь) в стационарном режиме.
Любая ЕО может находиться в состояниях «работа», «отказ» (идёт восстановление ЕО после её поломки), «блокировка» (невозможно передать обработанную деталь на следующую операцию) и «простой» (отсутствие деталей на входе ЕО). Пусть п — количество бункеров, т — число ЕО в линии. Через Z+ обозначим множество неотрицательных целых чисел. Введём вектор переменных Ь = ,...,Л,га) е где ^ - объем г-го бункера, г = 1,...,п. Структура линии представляется ориентированным последовательно-параллельным графом С, у которого вершины 60,..., 6га+1 отвечают бункерам, а дуги а1,..., ат — единицам оборудования. Каждая дуга направлена от вершины входного бункера соответствующей ЕО к вершине его выходного бункера. Возможны кратные дуги, т.е. параллельно работающие ЕО с общими бункерами. Единственная вершина 60, не имеющая входящих дуг, соответствует входному бункеру линии, а единственная вершина 6га+1 без выходящих дуг - выходному. Каждой вершине приписано максимальное число вмещаемых деталей в бункер ], причём Л,0 = Л,га+1 = то. Предполагается, что на входе линии всегда имеется до-
статочное число деталей и обработанные детали всегда могут быть помещены в выходной бункер, поскольку он имеет неограниченный объём.
Каждая дуга а графа С характеризуется тройкой параметров (Т°,Тв,и) ^ Находясь в состоянии «работа», ЕО с номером г имеет постоянное время обработки детали и, г = 1,...,т. Отказы и восстановления различных ЕО происходят независимо, и время наработки на отказ, также как и время восстановления ЕО, имеет экспоненциальное распределение. Для каждой ЕО ТО — среднее время наработки на отказ, а ТВ — среднее время восстановления. В состояниях «блокировка» и «простой» ЕО не отказывает.
Обозначим среднее число деталей в стационарном режиме через qj•(Ь). Тогда критерий максимизации имеет вид
п— 1
/(V(Ь),В(ЬШЬ)) = Тат£(V(Ь)) - 3(В(Ь)) - Тат ^ с,qj(Ь),
j=1
где Тат — период амортизации линии; 5(V(Ь)) — доход от продажи готовой продукции за единицу времени; с, — затраты на хранение одной детали в буфере ] в единицу времени; 3(В(Ь)) — капитальные затраты на установку бункеров. В качестве метода приближенного вычисления средней производительности линии V(Ь) и среднего числа деталей qj(Ь) в каждом бункере ] используется метод из [6], основанный на замене участков линии эквивалентными ЕО с применением марковских моделей для пар последовательных и параллельных ЕО.
2.2. Результаты экспериментов
В экспериментах с описанными выше задачами используется система окрестностей радиуса 1, порождённая метрикой Хэмминга. Для задач, построенных по модели ЫК-ландшафтов и задач отыскания двоичной последовательности с минимальной автокорреляцией и задач распределения ёмкостей бункерных устройств выполнялось по 1000 независимых испытаний алгоритма локального поиска. Для большинства задач вершинного покрытия этот параметр был увеличен до 10000.
В некоторых задача малой размерности удавалось осуществить полный перебор элементов пространства решений — в таких случаях построенные доверительные интервалы сравниваются с истинным значением V. В случаях, когда полный перебор оказывался чрезмерно трудоёмким, доверительные интервалы сравнивались с бутстреп-оценкой из [10,29], показавшей хорошие результаты в экспериментах [10, 17]. Описание бутстреп-метода из [29] приводится в приложении. При расчётах бутстреп-методом число независимых реализаций В всюду полагалось равным 1000.
Модель ЫК-ландшафтов. Эксперименты по оценке числа локальных оп-тимумов в задачах, построенных по модели ЫК-ландшафта, проводились при
различной «изрезанности ландшафта», определяемой величиной К в диапазоне от 2 до 12, в то время как размерность п =15 оставалась неизменной. Рис. 2 показывает истинное число локальных оптимумов и двухсторонние 95%-доверительные интервалы, вычисленные методом переписи Шнабе-ля в предположении изотропности р = р с использованием утверждения 1. Доверительные интервалы вокруг значений максимума правдоподобия и(т, к) рассчитаны в предположении изотропности. Рисунок также содержит нижние границы Ь(Ь, 0.95) и Ь'(Ь, 0.95) доверительного уровня 95%.
Рис. 2. Доверительные интервалы для числа локальных оптимумов в ^К-ландшафтах
График значений V на рис. 2 существенно отклоняется от оценки максимального правдоподобия по методу Шнабеля, что показывает неадекватность предположения об изотропности применительно к задачам данного класса. Нижние границы, основанные на ожидании первого повторения, оказываются даже ниже числа г обнаруженных оптимумов за к испытаний локального поиска.
Задачи отыскания двоичной последовательности с минимальной автокорреляцией. Рис. 3, подобно рис. 2, представляет результаты экспериментов с задачами отыскания двоичной последовательности с минимальной автокорреляцией. Длина последовательности п варьируется от 7 до 29. Для задач с п > 24 полный перебор пространства решений оказался чрезмерно трудоёмким, поэтому точное значение величины V для задач большой размерности не изображено. Как видно из рисунка, бутстреп-метод даёт оценку наиболее близкую к истинному значению V.
Из приведённого графика видно, что при малой размерности задач (п < 12) оценка максимального правдоподобия основанная на методе переписи Шна-беля, даёт достаточно близкие к истинным значения. Далее с увеличением раз-
Рис. 3. Оценки числа локальных оптимумов в задачах отыскания двоичной последовательности с минимальной автокорреляцией
мерности задач предположение изотропности, по-видимому, становится неадекватным. Оценка Ь(Ь, 0.95), в отличие от Ь'(т,к, 0.95), ведёт себя очень нестабильно.
Задачи о вершинном покрытии из коллекции DIMACS. Результаты экспериментов по оцениванию числа локальных оптимумов для задач вершинного покрытия представлены в табл. 1. Здесь столбец «95% интервал» содержит 95% доверительные интервалы, построенные с использованием утверждения 1 в предположении изотропности.
Как выяснилось, задачи большой размерности в коллекции DIMACS характеризуются большими значениями величины v, в связи с этим для получения содержательных результатов оценивания числа локальных оптимумов рассматривались только графы с числом вершин не более 200. При r = 1000 на всех тестовых примерах, кроме brock200-2 и brock200-4 не обнаруживалось ни одного повторения ранее найденных локальных оптимумов, а значит оценка, основанная на ожидании первого повторения, была не определена. Для других задач число испытаний r было увеличено до 10 000.
Задачи gn200p0.9-44, mann a27, san200-0.9-1 и san200-0.9-3 не представлены в табл. 1 в связи с тем, что и при r = 10000 вновь было получено k = r. В этих случаях нижняя граница может быть построена по формуле L(r + 1,0.95) « 0.167r2 = 1.67 ■ 108 (см. п. 1.1). Как видно из табл. 1, нижняя
Таблица 1. Оценки числа локальных оптимумов в ЗВП
Задача |V | r k L(t, 0.95) L'(r,k, 0.95) бутстреп 95% интервал
john8-4-4 70 104 4509 572 5260 7145 5248-5398
brock200-2 200 103 972 346 13263 33761 12214-26212
mann a9 45 104 9406 7824 7.6-104 1.5105 74808-87445
keller4 171 104 9805 1.4105 2.3 -105 4.9-105 2.2-105-2.9-105
brock200-4 200 103 999 8.1105 1.1105 8.3-105 9.0-104-2.0-107
john16-2-4 120 104 9939 1.6105 6.6-105 1.6106 6.4105-1.1106
san200-0.72 200 104 9945 1.2106 7.2-105 1.8106 7.0-105-1.2-106
san200-0.92 200 104 9995 4.1106 4.8-106 1.2107 4.3106-3.1 107
c125.9 125 104 9999 1.5107 1.1107 108 9.0106-1.1 109
gn200p0.9-55 200 104 9999 1.1 106 1.1107 108 9.0106-1.1 109
оценка, основанная на времени первого повторения, и в данной серии даёт неустойчивые результаты, в ряде случаев оказываясь ближе к оценке бутстре-па чем нижняя оценка, построенная по методу переписи Шнабеля (brock200-4, san200-0.72, c125.9).
Полученные данные показывают, что число локальных оптимумов в ЗВП близкой размерности может существенно различаться и зависит от структуры графа. Доверительный интервал, построенный по методу переписи Шнабе-ля в предположении изотропности, во многих случаях не содержит значения бутстреп-оценки. Таким образом, предположение изотропности в задачах о вершинном покрытии не оправдывается.
Задача распределения ёмкостей бункеров. Результаты оценивания числа локальных оптимумов в задачах распределения ёмкостей бункеров, созданных на основе реальных данных технологических линий автозавода Renault [6], показаны на рис. 4. Число промежуточных бункеров n и мощность множества допустимых решений в этих задачах приведена в табл. 2. Среди тестовых примеров 4 линии имеют последовательную структуру, а 2 линии AS7, AS8 имеют последовательно-параллельную структуру из 10 ЕО. Полный перебор пространства решений во всех задачах кроме AS3 оказался чрезмерно трудоёмким из-за значительного времени вычисления целевой функции. В связи с этим точное число локальных оптимумов для этих задач не вычислялось.
Как видно из рисунка, на этих задачах нижняя граница L'(r,k, 0.95), соответствующая равновероятному попаданию в локальные оптимумы, оказывается близка к бутстреп-оценке. Можно предположить, что структура локальных оптимумов задачи распределения ёмкостей бункерных устройств близка к изотропной. Причиной тому может служить симметричность «ландшафта» целевой функции в этих задачах, отмеченная в работе [5].
ав.З авЛ ав.б аБ.2 ав.4
«—к -•--|_'(г,к,0.95) -•--Щ,0.95) -■- бутстреп
Рис. 4. Оценки числа локальных оптимумов в задачах распределения ёмкостей бункеров
Таблица 2. Тестовые задачи распределения ёмкостей бункеров
п допустимых решений
AS1 4 4 х 105
AS2 9 9 х 1011
AS4 3 4 х 103
AS6 13 2 х 1022
AS7 7 9 х 107
AS8 7 5 х 108
3. Заключение
Проведённые вычислительные эксперименты на известных тестовых примерах средней размерности показывают адекватность предложенных нижних оценок и бутстреп-метода для приближенного оценивания числа локальных оптимумов. В то же время построенные двусторонние доверительные интервалы зачастую не покрывают истинное значение оцениваемого параметра, что говорит о неприменимости к рассматриваемым задачам предположения о равновероятном попадании в локальные оптимумы алгоритма локального поиска.
Для оценивания числа локальных оптимумов снизу в большинстве случаев предпочтительнее использовать нижнюю границу, основанную на методе переписи Шнабеля, т.к. нижняя граница, основанная на времени первого повторения локальных оптимумов, даёт нестабильные результаты. Тем не менее последняя оценка может оказаться единственным приемлемым методом, когда локальные оптимумы повторяются чрезвычайно редко и априори сложно выбрать число испытаний, при котором будет получено хотя бы одно повторение ранее полученного локального оптимума.
Дальнейшие исследования могут быть направлены на использование предложенных методов в правилах остановки для эвристик локального поиска с
многократным перезапуском. В частности, необходимое условие для остановки алгоритма может формулироваться как условие превышения числа найденных локальных оптимумов над нижней границей заданного доверительного уровня. Также предложенные методы могут использоваться для классификации тестовых примеров по их сложности для алгоритмов, основанных на локальном поиске. Методы оценивания числа локальных оптимумов могут также быть использованы в биоинформатике (см., например, [31]).
Благодарности
Работа выполнена при поддержке Программы фундаментальных научных исследований государственных академий наук на 2013-2020 годы, п. 1.5.1.5. «Исследование и решение задач комбинаторной оптимизации с использованием целочисленного программирования». Авторы благодарны С.А. Клокову за помощь в доказательстве утверждения 3.
Приложение
Настоящее приложение содержит описание бутстреп-метода для оценки числа локальных оптимумов, предложенного в [10,29]. Бутстреп-метод предложен Б. Эфроном для решения целого ряда задач статистического анализа [3], и в частности, для непараметрической оценки смещения. Пусть Xi,...,Xn -независимые одинаково распределённые случайные величины с неизвестным распределением F на некотором пространстве X. Предположим, что на основе наблюдаемых значений
Xi = x i, ..., Xr = xn
может быть вычислена оценка 9(x1,... ,xn) некоторого интересующего нас функционального параметра распределения 9(F) е R (например, математического ожидания, дисперсии или |supp F|, как в случае оценивания числа локальных оптимумов или особей в популяции). Интерес представляет распространённый случай, когда оценка имеет вид
9(xi ,...,xn) = 9(F),
где F — эмпирическое распределение, приписывающее массу n каждой из точек имеющейся выборки x1,...,xn. Требуется оценить величину смещения оценки 9 = 9(F):
bias = EF[9(F) - 9(F)], (5)
где Ef обозначает математическое ожидание в предположении, что X1,... ,Xn — независимые одинаково распределённые случайные величины с распределением F.
В бутстреп-методе для оценки смещения bias при фиксированном F рассматривается случайная выборка объема n:
X * = г* X * = г*
XI x1, . . . , Xn xn,
где X*,... ,X* — независимые одинаково распределённые случайные величины с вероятностным распределением F. Выборочное распределение интересующей нас случайной величины 9(F) — 9(F) аппроксимируется бутстреп-распределением случайной величины 9(X*,..., X*) — 9(F). Для получения бутстреп-оценки смещения BIAS* в формуле (5) производится замена неизвестного распределения случайной величины 9(F) — 9(F) имеющимся бутстреп-распределением:
BIAS* = E*[(9(Xi,...,X*) — 9(F)],
где E* обозначает математическое ожидание в предположении, что X*,...,X* — независимые одинаково распределённые случайные величины с вероятностным распределением F, зафиксированным на его наблюдаемом значении.
Величину BIAS* в некоторых случаях удаётся вычислить точно или приближённо теоретическими методами (см., например, гл. 1 в [3]), однако, наиболее часто для этого используется аппроксимация бутстреп-распределения методом Монте-Карло, где генерируется достаточно большое число B независимых реализаций x 1b,..., x*nb, b = 1,..., B, случайной выборки X**,..., X* из распределения F. В таком случае бутстреп-оценка смещения аппроксимируется величиной
1 B
BIAS* = -^Y, 9(x*b,..., x*nb) — 9(F). b= 1
Алгоритм 1. Бутстреп-оценка смещения
1. Построить эмпирическое распределение F, приписывающее массу 1 каждой из точек имеющейся выборки x1, . . . , xn.
2. Для b = 1,..., B независимо выполнять:
2.1. Построить бутстреп-выборку из F:
X * = x* X * = x*
XI xib, . . . , Xn xnb.
2.2. Вычислить 9* = 9(x1b,... ,xnb).
3. Вычислить бутстреп-оценку смещения
1 B
BIAS* = B E
b=1
4. Конец работы.
С использованием бутстреп-оценки смещения, найденной алгоритмом 1, может быть получена скорректированная бутстреп-оценка:
9* = 9 — BIAS*. (6)
Применительно к задаче оценки числа локальных оптимумов описанный алгоритм позволяет оценивать смещение статистики V(г, к), вычисленной по
методу максимального правдоподобия. В данном случае в = v, X = {1,...,v} - множество всех локальных оптимумов, распределение F описывает случайную величину, равную номеру локального оптимума, получаемого алгоритмом локального спуска из случайно выбранной начальной точки. На шаге 1 алгоритма 1 по n = r испытаниям локального спуска строится эмпирическое распределение F = (vl/r,...,vk/r), где Vj — число испытаний, при которых был получен локальный оптимум с номером j, j = 1,...,k.
На шаге 2.2. полагаем V* = v(r, k*), где k* — число различных значений в бутстреп-выборке (j,.-,j*) на итерации b. На шаге 3 = вb = 1,...,B и оценка смещения вычисляется по формуле
в
BIAS* = - v (r, k))/B.
b=l
При нарушении условия изотропности локальные оптимумы, куда вероятность попадания локального спуска мала, имеют мало шансов быть представленными в эмпирическом распределении F. В таком случае F оказывается ближе к равномерному распределению по сравнению с истинным распределением вероятностей F. В результате скорректированная бутстреп-оценка (6) отклоняется в сторону оценки максимального правдоподобия û(r,k). Приведённые в [10] результаты эксперимента подтверждают это.
Литература
1. Гэри М., Джонсон Д. Вычислительные машины и труднорешаемые задачи. М. : Мир, 1982.
2. Зубков А.М., Попов Н.Н. Отношение частичного порядка, порождённое распределениями числа занятых ячеек // Матем. заметки. 1982. T. 32, № 1. С. 97-102.
3. Эфрон Б., Тибширани Г.Г.Р. Нетрадиционные методы многомерного статистического анализа. М. : Финансы и статистика, 1988.
4. Севастьянов Б.А. Курс теории вероятностей и математической статистики. М. : Наука, 1982.
5. Сигаев В.С. О свойствах локальных оптимумов задачи размещения буферных накопителей // Вестник Омского университета. 2006. № 4. С. 1-3.
6. Ancelin B., Semery A. Calcul de la productivité d'une ligne integrée de fabrication // RAIRO Automatiq., Productiq. Informatiq. Industrielle. 1987. V. 21. P. 209-238.
7. Craig C.C. Use of marked specimens in estimating populations // Biometrika. 1953. V. 40, N. 1-2. P. 170-176.
8. Darroch J.N. The multiple-recapture census. I: estimation of a closed population // Biometrika. 1958. V. 45, N. 3-4. P. 343-359.
9. Dolgui A., Eremeev A., Kolokolov A., Sigaev V. A genetic algorithm for the allocation of buffer storage capacities in a production line with unreliable machines // Journal of Mathematical Modelling and Algorithms. 2002. V. 1, N. 2. P. 89-104.
10. Eremeev A.V., Reeves C.R. Non-parametric estimation of properties of combinatorial landscapes // Applications of Evolutionary Computing: Proceedings of EvoWorkshops
2002 / Cagnoni S., Gottlieb J., Hart E., Middendorf M., Raidl G. Lecture Notes in Computer Science, Vol. 2279. Berlin : Springer-Verlag, 2002. P. 31-40.
11. Eremeev A.V., Reeves C.R. On confidence intervals for the number of local optima // Applications of Evolutionary Computing, EvoWorkshop 2003: EvoBIO, EvoCOP, EvolASP, EvoMUSART, EvoROB, and EvoSTIM / Raidl G.R., Meyer J.-A., Middendorf M., Cagnoni S., Cardalda J.J.R., Corne D., Gottlieb J., Guillot A., Hart E., Johnson C.G., Marchiori E. Lecture Notes in Computer Science, Vol. 2611. Heidelberg : Springer, 2003. P. 224-235.
12. Garnier J., Kallel L. How to detect all maxima of a function? // Proceedings of the Second EVONET Summer School on Theoretical Aspects of Evolutionary Computing (Anvers, 1999). Berlin : Springer, 2001. P. 343-370.
13. Garthwait P.H., Buckland S.T. Analysis of multiple recapture census by computing conditional probabilities // Biometrics. 1990. V. 46, N. 1. P. 231-238.
14. Glover F. and Laguna M. Tabu Search. Norwell, MA: Kluwer, 1997.
15. Golay M.J.E. Series for low-autocorrelation binary sequences. IEEE Trans. Inform. Theory. 1977. V. 23. P. 43-51.
16. Hansen P., Mladenovic N. Variable neighborhood search: Principles and applica- tions // European Journal of Operational Research. 2001. V. 130. P. 449-467.
17. Hernando L., Mendiburu A., Lozano J.A. An evaluation of methods for estimating the number of local optima in combinatorial optimization problems // Evol. Comput. 2013. V. 21, N. 4. P. 625-658.
18. Johnson D.S., Trick M.A. Introduction to the second DIMACS challenge: cliques, coloring, and satisfiability // DIMACS Series in Discrete Math. and Theoretical Comput. Sci. / Johnson D., Trick. M. Providence, RI : American Math. Soc., 1996. Vol. 26. P. 1-10.
19. Kauffman S.A. Adaptation on rugged fitness landscapes // Lectures in the Sciences of Complexity, Vol. I of SFI Studies. Addison-Wesley, 1989. P. 619-712.
20. Kochetov Yu., Mikhailova A., Plyasunov A. A Genetic Local Search Algorithm for the Graph Partitioning Problem with Cardinality Constraints // Preprints of the 13th IFAC Symposium on Information Control Problems in Manufacturing (INCOM09), Moscow, Russia, June 3-5, 2009. P. 1991-1996.
21. Mao C.X. On the nonidentifiability of population sizes // Biometrics. 2008. V. 64. P. 977-981.
22. Mood A.M., Graybill F.A., Boes D.C. Introduction to the theory of statistics. Third Ed. New York : McGraw-Hill, 1973.
23. Pickands J., Raghavachari M. Exact and asymptotic inference for the size of a population // Biometrika. 1987. V. 74, N. 2. P. 355-363.
24. Pledger S. The performance of mixture models in heterogeneous closed population capture-recapture // Biometrics. 2005. V. 61, N. 3. P. 868-876.
25. Reeves C. Genetic algorithm for the operations researcher // INFORMS Journal on Computing. 1997. V. 9, N. 3. P. 231-250.
26. Reeves C.R. The "Crossover Landscape"and the Hamming Landscape for Binary Search Spaces // Foundations of Genetic Algorithms 7 / De Jong K., Poli R., Rowe J. San Francisco, CA: Morgan Kaufmann. 2003. P. 81-98.
27. Reeves C.R. Estimating the Number of Optima in a Landscape, Part I: Statistical Principles. Coventry University Technical Report SOR#01-03. 2001.
28. Reeves C.R. Estimating the Number of Optima in a Landscape, Part II: Experimental Investigations. Coventry University Technical Report SOR#01-04. 2001.
29. Reeves C.R., Eremeev A.V. Statistical analysis of local search landscapes // J. Oper. Res. Soc. 2004. V. 55, N. 7. P. 687-693.
30. Resende M.G.C. and Ribeiro C.C. Greedy randomized adaptive search procedures // Handbook of metaheuristics / Glover F. and Kochenberger G. N.Y.: Kluwer Academic Publishers, 2003. P. 219-249.
31. Sahoo S., Albrecht A.A. Approximating the set of local minima in partial RNA folding landscapes // Bioinformatics. 2012. V. 28, Issue 4. P. 523-530.
32. Schnabel Z.E. The estimation of the total fish population of a lake // American Math. Monthly. 1938. V. 45. P. 348-352.
33. Seber G.A.F. The estimation of animal abundance. London : Charles Griffin, 1973.
34. Yagiura M., Kishida M., Ibaraki. T. A 3-flip neighborhood local search for the set covering problem // Eur. J. Oper. Res. 2006. V. 172. P. 472-499.
ON CONFIDENCE INTERVALS FOR THE NUMBER OF LOCAL OPTIMA
A.V. Eremeev1
Associate Professor, Dr. Sci. (Phys.-Math.), e-mail: eremeev@ofim.oscsbras.ru
C.R. Reeves2 Professor Emeritus, e-mail: c.reeves@coventry.ac.uk
1Omsk Branch of Sobolev Institute of Mathematics SB RAS 2Coventry University
Abstract. The number of local optima is an important indicator of optimization problem difficulty for local search algorithms. Here we will discuss some methods of finding the confidence intervals for this parameter in problems where the large cardinality of the search space does not allow exhaustive investigation of solutions. Computational results are reported that were obtained by using these methods for N K landscapes model of S.A. Kauffman, for the low autocorrelation binary sequence, for buffer allocation problems in production line, and vertex cover problems.
Keywords: local search, combinatorial optimization problem, Schnabel census, conservative confidence interval.
Дата поступления в редакцию: 02.02.2017
