Решение систем нелинейных уравнений методами интервального распространения ограничений Текст научной статьи по специальности «Математика»

Вычислительные технологии

Том 7, № 2, 2002

РЕШЕНИЕ СИСТЕМ НЕЛИНЕЙНЫХ УРАВНЕНИЙ МЕТОДАМИ ИНТЕРВАЛЬНОГО РАСПРОСТРАНЕНИЯ ОГРАНИЧЕНИЙ

М. Лоенко

Новосибирский филиал Российского научно-исследовательского института искусственного интеллекта, Россия e-mail: ml@iis.nsk.su

New heuristic algorithm for obtaining 3B-compatibility, i. e. the algorithm for solution correction (SC-algorithm) is presented. Standard 3B-algorithm does not complete its work in some problems in the reasonable time because of the peculiarities of the 2B-filtration function and fixed order of variable choice. In the presented algorithm 2B-filtration is stopped at exceeding of some limiting time and the order of variable choice is defined dynamically.

Введение

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

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

Известны несколько определений локальной совместности NCSP, в частности достаточно слабая 2В-совместность, более сильные 3B-, 4В-совместности и т.д. Алгоритмы достижения таких совместностей описаны в [3]. Чаще всего используется 2В-алгоритм, так как

© М. Лоенко, 2002.

для большинства задач он завершает свою работу достаточно быстро. Несмотря на то, что с помощью 3В-алгоритма можно получить более точную оценку решения, он используется реже, поскольку требует гораздо больше времени.

Целью данной работы является построение алгоритма быстрого достижения 3В-сов-местности, основанного на использовании эвристик.

Статья имеет следующую структуру. В разд. 1 даются понятие МСБР, определения 2В- и ЗВ-совместности для МСБР над непрерывными областями, здесь также продемонстрированы примеры алгоритмов, с помощью которых достигаются 2В- и ЗВ-совместности. В разд. 2 приводится разработанный автором БС-алгоритм, являющийся альтернативным методом достижения 3В-совместности, приведены результаты тестов и сравнений.

1. Основы теории распространения ограничений

1.1. Понятие N08?

Определение 1. Численной задачей удовлетворения ограничений ЫСЯР будем называть набор М = (Х,О, С), где X — множество переменных {х\, ...,хп}; О = О\ х ... х Оп (Ог — множество значений переменной хг); С — множество ограничений {С\,..., Ст}, С^ — отношение (уравнение, неравенство, таблица), которое связывает некоторые переменные х 1,Хп^ .

Определение 2. Решением ЫСЯР будем называть все векторы А = {а\, ...,ап} Е О, для которых выполняются все ограничения С).

Определение 3. ЫСЯР М = (Х,О,С) называется глобально совместной тогда и только тогда, когда любой вектор а Е О принадлежит решению задачи М.

Задача нахождения множества решений МСБР или, что то же самое, построения эквивалентной глобально совместной МСБР является МР-трудной, поэтому существует множество понятий более слабой локальной совместности — совместность по дугам, путям для дискретных задач [4], совместность по границам для интервальных задач [3]. Для задач достижения этих локальных совместностей, т. е. нахождения эквивалентных локально совместных МСБР, известны алгоритмы, имеющие полиномиальную сложность.

При работе с непрерывными переменными очень удобно представление множества значений переменных с помощью интервалов (оно может быть также использовано и для целочисленных переменных). Итак, пусть ЕР — множество чисел с плавающей точкой (чисел, представимых в компьютере), расширенное двумя элементами {-то, Тогда

любое подмножество множества вещественных чисел может быть представлено минимальным интервалом I = [а,Ь], где а Е ЕР, Ь Е ЕР, таким, что X С I. Методы решения численных задач, которые используют технику распространения ограничений при интервальном представлении множества значений переменных, называются методами интервального распространения ограничений.

Для задач с непрерывными переменными введены свои виды локальной совместности, которые являются аналогами совместностей по дугам в дискретном случае.

1.2. 2В-совместность

Приведем определение 2В-совместности [3]. Обозначим через а+ минимальное ЕР-число, строго большее а, и через а- — максимальное ЕР-число, строго меньшее а.

Определение 4. Пусть M = (X,D,C) есть некоторая NCSP. M называется 2B-совместной тогда и только тогда, когда для любого ограничения Cj,...,xjk) £ C и любой переменной xji £ {xj1 ,...,xjk}. Если [a,b] — интервал, соответствующий области значений Dji, то существуют значения a¡ £ D¡ (l = jl,..., ji-1, jí+l, ■■■,jk) и значение aji £ [a,a+) такие, что Cj(aj1 ,...,ajk) удовлетворяется, и существуют значения b¡ £ D¡ (l = ji,jí-l, jí+i,jk) и значения bji £ (b-,b] такие, что Cj(bj1 ,...,bjk) удовлетворяется.

Рассмотрим алгоритм достижения 2В-совместности. Поскольку результатом работы алгоритма является сужение начальных областей, алгоритм называют также 2В-фильт-рацией и обозначают через Ф2В (M), где M есть NCSP.

Сначала все ограничения разбиваются на простейшие. Формально это описано в [5], здесь же мы поясним на примере, как это делается. Рассмотрим неравенство x2-siny+xt — 4 > 0. Введением дополнительных переменных оно разбивается на несколько простейших вида

f (x)=l y, x +l y =l z, (1)

где =l£ {=, >}; +l £ {+, —, *,/}; f £ {sin, cos, tg, arcsin, arccos, arctg, exp, log, pow, sqrt}.

Получаем

Ql = x2, q2 = siny, q3 = xt, qA = Qi — q2, + Q3 > 4.

Заметим, что при таком разбиении в каждом ограничении каждая переменная очевидным образом выражается через остальные. Таким образом, мы получили новую NCSP, решив которую мы очевидным образом сможем получить решение начальной NCSP. Начиная с этого момента будем считать что все ограничения в NCSP имеют вид (1).

Ядром алгоритма 2В-фильтрации является набор фильтрующих функций. Аргументами фильтрующей функции являются некоторое ограничение и некоторая переменная, которая входит в это ограничение. Эта функция сужает интервал — область возможных значений выбранной переменной в зависимости от областей значений всех переменных, входящих в ограничение. Рассмотрим ее работу на примере x = y — z.

Function filter(< op, x, y, z >)

begin

if (op = '—'') then

ret := 0;

res- := min (Dy) - max(Dz);

res+ := max (Dy) - min(Dz);

res := [res-, res+];

Dnew • = Dx ^ res ;

if (Dnew = 0) then

print(INCONSISTENSY);

EXIT

end if

if ((Dx\Dnew)= 0) then

ret : = ret U {x};

Dx • = Dnew j

end if

res- := min (Dx) + min(Dz);

res+ := max (Dx) + max(Dz);

res := [res-, res+];

Dnew * Dy П VGS ,

if (Dnew = 0) then

print(INCONSISTENSY); EXIT end if

if ((Dy\Dnew) = 0) then

ret : = ret U {y};

Dy ■ = Dnew >

end if

res- : = min (Dy) - max(Dx); res+ := max (Dy) - min(Dx);

res := [res-, res+];

Dnew • = Dz П res j

if (Dnew = 0) then print(INCONSISTENSY); EXIT end if

if ((Dz\Dnew) = 0) then

ret := ret U {z};

Dz • = Dnew >

end if return ret; end if

end

Функция filter(< op,x,y,z >) возвращает набор изменившихся переменных. Примеры обработки остальных операций и элементарных функций приведены в [5].

Алгоритм 2В-фильтрации работает следующим образом. Сначала мы ставим в очередь Queue все отношения. Один шаг алгоритма состоит в том, что сначала выталкивается из очереди Queue некое отношение c, для него вызывается функция filter. В случае, если функция возвратит непустое множество изменившихся переменных, в очередь Queue добавляются все отношения, в которых участвует хотя бы одна из изменившихся переменных, кроме тех, которые уже стоят в очереди. Алгоритм заканчивается, если очередь Queue пуста:

Function 2B-Filter(X, D, C) begin

Queue := 0 for each c £ C

Queue. put (c); end for

while Queue = 0 do c := Queue .get (); Changed := filter(c); for each x £ Changed for each c £ C

if ((x £ c ) AND NOT(c £ Queue)) Queue ,put(c );

end if end for end for end while end

Более подробное описание 2В-алгоритма и доказательство его сходимости приведены в [3]. Здесь стоит отметить, что при вычислении функции сужения на некоторых задачах наблюдается эффект медленной сходимости [6]. Чтобы гарантировать, что процесс завершится за реальное время, мы только при существенном уточнении переменной будем ставить в очередь все связанные с ней функции фильтрации, т. е. нам потребуется введение некоторого абсолютного и (или) относительного ограничения для того, чтобы игнорировать малые уточнения переменных. Правда, в этом случае мы можем не получить в итоге 2В-совместную систему, поскольку требования из определения 4 будут несколько ослаблены и полученная оценка решения будет несколько грубее.

Таким образом, работа описанного алгоритма зависит от некоторого параметра скрупулезности. Изменяя его, мы можем уменьшить время работы алгоритма в ущерб качеству получаемой оценки решения.

1.3. 3В-совместность

Определение 5. Пусть M = (X,D, C) есть некоторая NCSP и x — переменная из X с областью возможных значений [a,b]. Пусть также:

M1 — NCSP получена из M сужением Dx в D на D\ = [а,а+]; M2 — NCSP получена из M сужением Dx в D на DX = [Ь-,Ь]; Dx SB-совместна, если и только если 1) Фхв(M1) = 0, 2) Фхв(M2) = 0.

NCSP (X,D,C) SB-совместна тогда и только тогда, когда все ее области из D SB-совместны.

По аналогии с 2В-фильтрацией существует стандартный алгоритм достижения 3В-совместности, который называется ЗВ-фильтрацией и обозначается через Фзв(M), где M есть NCSP. Рассмотрим этот алгоритм.

Ядром алгоритма являются функции-корректоры переменных, которые применяются к области возможных значений некоторых переменных и сужают ее до тех пор, пока она не станет ЗВ-совместной. Эти функции вызываются отдельно для левых границ, отдельно для правых границ и выглядят следующим образом: Function 3B-Left-Bound(x) begin

stack := 0; stack .push(D) ; while stack = 0 do D := stack.pop(); D" := Фхв(-X, D', C); if D" = 0 then if width(DX') > 0 then

(LD, RD) := bisect (D ', x) stack .push(RD); stack .push(LD) ;

else

Dnew = [ min(DX), max(Dx)]; changed := (Dx\Dnew) = 0; Dx := Dnew; return changed; end if end if end while

print(INCONSISTENCY); EXIT end

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

Рассмотрим теперь, как выглядит ЗВ-фильтрующая функция. Эта функция запускает функции 3B-Left-Bound() и 3B-Right-Bound() для всех переменных задачи до тех пор, пока не обнаружится несовместность или эти функции не станут выдавать < false >. Function 3B-Filter() begin do

changed := < false >; for each x £ X

changed := changed OR 3B-Left-Bound(x); changed := changed OR 3B-Right-Bound(x); end for while (changed); end

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

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

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

2. Алгоритм коррекции решения

Излагаемый в данной главе алгоритм будем называть алгоритмом коррекции решения или БС-алгоритмом, а построенную на его основе функцию фильтрации — БС-фильтрацией.

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

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

Для того чтобы сформулировать условия прерывания корректоров, нам понадобится определение текущей расчетной скорости корректора.

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

Теперь мы сможем сформулировать условия прерывания работы корректоров:

1) если время с момента его вызова превышает некоторое значение ¿сац;

2) если время работы вызванной им функции Ф2в () превышает некоторое значение ;

3) если текущая расчетная скорость корректора станет меньше некоторого значения

^СОГГ.

Значения ¿сац, и зСОГГ зависят от размера решаемой задачи и вычисляются по некоторым эмпирическим формулам.

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

2.1. Настроечные параметры алгоритма

Перечислим параметры, которые влияют на работу БС-алгоритма. Выбор эвристик зависит от конкретной реализации алгоритма и приоритетов пользователя. Изменяя параметры, можно влиять на скорость решения, качество оценки и даже на количество используемой памяти. Для получения приведенных в следующем разделе результатов автором использовались следующие эвристики (согласно определению 1 т — количество отношений в системе):

1) формула вычисления ¿сац (¿сац = 300т);

2) формула вычисления ¿2в (¿2в = 80т);

3) формула вычисления штрафа.

После возврата из корректора на него накладывается штраф, равный 10, если корректор не изменил значение корректируемой переменной; равный 1, если он изменил значение менее чем на 5 %, и равный 0 — в других случаях.

Эти формулы были получены в результате тестирования алгоритма на большом количестве задач разного типа и являются, на взгляд автора, наиболее универсальными.

2.2. Результаты экспериментов

В таблице приведено время работы 3В- и БС-алгоритмов в секундах на пяти примерах из Приложения А. Эти результаты получены на РеП;тт-166.

Пример 3B SC

1 46 3

2 1198 36

3 28 <1

4 1228 <1

5 25 20

Как видно из примеров, скорость работы SC-алгоритма существенно выше скорости 3В-алгоритма. Это происходит потому, что в 3В-алгоритме очень многое зависит от порядка выбора переменных. В SC-алгоритме эта зависимость сведена к минимуму.

Заключение

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

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

Список литературы

[1] BENHAMou F. Interval constraint logic programming // Constraint Programming: Basic and Trends / A. Podelski (Ed.). LNCS. 1995. Vol. 910. P. 1-21.

[2] Older W., Vellino A. Constraint arithmetic on real intervals // Constraint Logic Programming: Selected Res. / F. Benhamou and A. Colmerauer (Eds). MIT Press, 1993. P. 175-195.

[3] Lhomme O. Consistency techniques for numeric csp's // Proc. 13th IJCAI / R. Bajcsy (Ed.). IEEE Computer Society Press, 1993. P. 232-238.

[4] Tsang E. Foundations of Constraint Satisfaction. N. Y.: Acad. Press, 1993.

[5] Kashevarova T., Leshchenko A., Petunin D., Semenov A. Combining various techniques with the algorithm of subdefinite calculations // Proc. 3rd Intern. Conf. of PACT'97. London, 1997. P. 287-306.

[6] Gotlieb A., Lhomme O., Rueher M., Taillibert P. Boosting the interval narrowing algorithm // Proc. JICSLP'96. MIT Press, 1996. P. 378-392.

[7] Boizumault P., JuSSIEN N. Dynamic Backtracking with Constraint Propagation Application to Static and Numeric CSPs.

[8] Prosser P. MAC-CBJ: maintaining arc consistency with conflict-directed backjumping.

Приложение A

Пример 1. [7]

x1 + x2 + x3 + x4 = 1, x1 + x2 — x3 + x4 = 3,

xi H- x2 H- x2 H- x4 — 4, xi + x2 + x3 + x2 = 2xi + 3, x1 = [—10,10],

x2 = [0, 10], x3 = [—10,10],

x4 = [—10, 10].

Пример 2. [7]

—x3x1ox11 — x5x10x11 — x7x10x11 + x4x12 + x6x12 + x8x12 = 0.4077,

x2x4x9 + x2x6x9 + x2x8x9 + x1x10 = 1.9115,

x3x9 + x5x9 + x7x9 = 1.9791,

3x2x4 + 2x2x6 + x2 x8 = 4.0616,

3x1x4 + 2x1x6 + x1 x8 = 1.7172,

3x3 + 2x5 + x7 = 3.9701,

x^ -И x2+1 — 1, i odd,

x1 > 0,

x10 > 0,

x11 > 0.

Пример 3.

xy + t — 2z = 4, x sin(x) + y cos(t) = 0, x — y + cos(z)2 = sin(t)2, x H yz = 2t.

Пример 4.

2

3X — 2y =77,

3X/2 — 2У2/2 = 7

Пример 5. ^royden banded [8]) xi := [—1,1], i = 1...20,

x1(2 + 5x1) + 1 — x2(1 + x2) = 0,

x2(2 + 5x2) + 1 — x1(1 + x1) — x3(1 + x3) = 0,

x3(2 + 5x2) + 1 — x1(1 + x1) — x2(1 + x2) — x4(1 + x4) = 0,

x4(2 + 5x4) + 1 — x1(1 + x1) — x2 (1 + x2) — x3(1 + x3) — x5 (1 + x5) = 0,

x5 (2 + 5x5) + 1 — x1(1 + x1) — x2 (1 + x2) — x3(1 + x3) — x4 (1 + x4) — x6 (1 + xa) = 0,

x6(2 + 5x2) + 1 — x1(1 + x1) — x2 (1 + x2) — x3(1 + x3) — x4 (1 + x4) — x5 (1 + x5) — —x7(1 + x7) = 0,

х7(2 + 5х2) + 1 — Х2(1 + Х2) — Хз(1 + Хз) — Х4(1 + Х4) — Хб(1 + Х5) — Хб(1 + Хб) —

— Х8 (1 + Хз) = 0,

Хз(2 + 5х2) + 1 — Хз(1 + Хз) — Х4 (1 + Х4) — Х5 (1 + Х5) — Хб(1 + Хб) — Х7(1 + Х7) —

— Хд (1 + Хд) = 0,

Хд (2 + 5х2) + 1 — Х4 (1 + Х4) — Х5 (1 + Х5) — Хб(1 + Хб) — Х7(1 + Х7) — Хз(1 + Хз) —

—Х1о(1 + Х10 Х1о(2 + 5х1о) + 1

— ХИ(1 + Х11

хп(2 + 5х?1) + 1 —Х1о(1 + Х10 Х12(2 + 5х12) + 1

—Х1з(1 + Х13

Х1з(2 + 5х1з) + 1

— Х14(1 + Х14

Х14(2 + 5х14) + 1

— Х15(1 + Х15

— Х1д(1 + Х19

= 0,

— Х5 (1 + Х5) — Хб(1 + Хб) — Х7(1 + Х7) — Хз (1 + Хз) — Хд (1 + Хд ) —

= 0,

— Хб(1 + Хб) — Х7 (1 + Х7) — Хз(1 + Хз) — Хд (1 + Хд) —

— Х12 (1 + Х12) = 0,

— Х7(1 + Х7) — Хз(1 + Хз) — Хд(1 + Хд) — Хю(1 + Х10) — Хц (1 + Хц) —

= 0,

— Хз(1

= 0,

— Хд(1

= 0,

+ Хз) — Хд (1 + Хд) — ХШ(1 + Х10) — Хц(1 + Хц) — Х12(1 + Х12) — + Хд) — Хю(1 + Х10) — Хц(1 + Хц) — Х12 (1 + Х12)х1з(1 + Х1з) —

Х15(2 + 5хУ + 1 — Хю( + Х10) — ХИ( + Х11) — Х12 ( + Х12) — Х1з( + Х1з) — Х14 (1 + Х14)

—Х1б(1 + Х1б) = 0,

Х1б(2 + 5х1б) + 1 — Х11 ( + Х11) — Х12 ( + Х12) — Х1з ( + Х1з) — Х14 ( + Х14) — Х15 (1 + Х15)

— Х17(1 + Х17) = 0,

Х17(2 + 5х17) + 1 — Х12 ( + Х12) Х1з( + Х1з) — Х14(1 + Х14) Х15 ( + Х15) — Х1б(1 + Х1б)

—Х1з(1 + Х1з) = 0,

Х1з(2 + 5х1з) + 1 — Х1з( + Х1з) — Х14 ( + Х14) — Х15 ( + Х15) Х1б( + Х1б) — Х17(1 + Х17)

— Х1д(1 + Х1д) = 0,

Х1д(2 + 5х1д) + 1 — Х14 ( + Х14) Х15 ( + Х15) — Х1б ( + Х1б) — Х17 ( + Х17) — Х1з(1 + Х1з)

— Х20(1 + Х20) = 0,

Х20(2 + 5х20) + 1 — Х15 ( + Х15) Х1б( + Х1б) — Х17 ( + Х17) Х1з( + Х1з) —

-

-

-

-

-

0.

Поступила в редакцию 17 сентября 2001 г.