Задача о сумме подмножеств11.04.2024 20:15

Задача о сумме подмножеств в общей формулировке звучит так:

Существует множество S чисел. Вопрос состоит в том, будет ли сумма некоторого подмножества от S равна заданному числу Т.

Известно, что данная задача NP-полная.

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

Частным случаем задачи о сумме подмножеств является задача разбиения множества чисел:

Множество чисел S необходимо разбить на два подмножества S1 и S2, где сумма S1 равна сумме S2.

(От задачи о сумме подмножеств текущая отличается только тем, что T = Sum (S1) / 2 = Sum (S2) / 2)

Хочу предложить вам простой и элегантный способ относительно быстрого решения обеих задач методом целочисленного линейного программирования (ЦЛП). Мы получим не только точный ответ на вопрос, но и найдём искомое подмножество.

И так, у нас есть множество S натуральных чисел:

$S = \left\{ S_{i} \right\}; S_{i} \in \mathbb{Z};$

Нам необходимо ответить на вопрос: Возможно ли используя элементы данного множества получить сумму равную T?

Для нахождения ответа нужно решить линейное уравнение:

$\sum_{i=1}^{n}S_{i}*X_{i}=T; X_{i}\in\{0,1\};\ S_{i}\in \mathbb{Z} ;T\in \mathbb{Z};$

Если решение будет найдено, то ответ на вопрос положительный, если же решение невозможно, то отрицательный. Если решений несколько, нас устроит любое.

Используя Python и библиотеку PuLP, проиллюстрируем на примере.

import pulp as pl
src = [59, 63, 15, 43, 75, 72, 61, 48, 71, 55]
n = len(src)

# Искомое значение
val = sum(src) // 2
# Объявляем модель
model = pl.LpProblem(name="Subset sum problem", sense=pl.LpMinimize)       
# Подключаем решатель CBC
solver = pl.PULP_CBC_CMD(msg=True)
# Объявляем переменные
x = [pl.LpVariable(name=f'x{i:05}', cat='Binary') for i in range(n)]    
# Целевая функция
model += pl.lpSum([x[i] for i in range(n)])    
# Основное ограничение
model += pl.lpSum([x[i] * src[i] for i in range(n)]) == val

print(model)

В результате мы получим следующую модель решения линейного уравнения:

Subset sum problem:
MINIMIZE
1*x00000 + 1*x00001 + 1*x00002 + 1*x00003 + 1*x00004 + 1*x00005
 + 1*x00006 + 1*x00007 + 1*x00008 + 1*x00009 + 0
SUBJECT TO
_C1: 59 x00000 + 63 x00001 + 15 x00002 + 43 x00003 + 75 x00004 + 72 x00005
 + 61 x00006 + 48 x00007 + 71 x00008 + 55 x00009 = 281
VARIABLES
0 <= x00000 <= 1 Integer
0 <= x00001 <= 1 Integer
0 <= x00002 <= 1 Integer
0 <= x00003 <= 1 Integer
0 <= x00004 <= 1 Integer
0 <= x00005 <= 1 Integer
0 <= x00006 <= 1 Integer
0 <= x00007 <= 1 Integer
0 <= x00008 <= 1 Integer
0 <= x00009 <= 1 Integer

Запускаем решатель и выводим решение

# Находим решение
status = model.solve(solver)
res = []
if status == 1:    
    for j, v in enumerate(model.variables()):
        if v.value() > 0:
            res.append(src[j])
print('Результат:', sum(res), res)

Результат: 281 [63, 75, 72, 71]

Что и требовалось получить.

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

Например для поиска решения следующего списка из 24 элементов потребовалось порядка двадцати секунд.

source = [1485498, 1643002, 1391445, 1280110, 1145269, 1195187, 1908655, 1709512, 
          1006747, 1354760, 1527205, 1486247, 1941933, 1634072, 1084740, 1350249, 
          1581194, 1981871, 1646604, 1734106, 1042882, 1763564, 1397430, 1177643]

Возможно ли как-то ускорить процесс нахождение решения?

Именно ответ на этот вопрос и побудил меня написать данный пост. Так как задача ЦЛП является NP-сложной, то считается что алгоритм, который находит решение является экспоненциальным от числа множеств.

Поэтому предлагаю «есть слона по кусочкам» (разбить сложную задачу на более простые и менее объемные части).

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

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

Y = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24]

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

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

Y = [1, 24, 2, 23, 3, 22, 4, 21, 5, 20, 6, 19, 7, 18, 8, 17, 9, 16, 10, 15, 11, 14, 12, 13]

Читать нужно следующим образом: Сумма значений любых Y элементов множества должна быть равна искомому T, где Y — значение из списка выше.

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

Отличие от кода предыдущего решения будет только в том, что мы будем искать не одно решение, а несколько:

loop = list(range(1, n + 1))
loop = [loop.pop(0) if i % 2 == 0 else loop.pop(-1) for i in range(n)]
for val in loop:
    # Устанавливаем дополнительное ограничение шага цикла
    model.constraints['Y'] = pl.lpSum([x[j] for j in range(n)]) == val 
    print(model)
    
    # Находим решение
    status = model.solve(solver)  
    if status == 1:
        res = []
        for j, v in enumerate(model.variables()):
            if v.value() > 0:
                res.append(src[j])
            return res
return []

Таким образом мы последовательно будем строить ряд моделей и пытаться найти решение.

Subset sum problem:
MINIMIZE
1*x00000 + 1*x00001 + 1*x00002 + 1*x00003 + 1*x00004 + 1*x00005 + 1*x00006
 + 1*x00007 + 1*x00008 + 1*x00009 + 1*x00010 + 1*x00011 + 1*x00012 + 1*x00013
 + 1*x00014 + 1*x00015 + 1*x00016 + 1*x00017 + 1*x00018 + 1*x00019 + 1*x00020
 + 1*x00021 + 1*x00022 + 1*x00023 + 0
SUBJECT TO
_C1: 1485498 x00000 + 1643002 x00001 + 1391445 x00002 + 1280110 x00003
 + 1145269 x00004 + 1195187 x00005 + 1908655 x00006 + 1709512 x00007
 + 1006747 x00008 + 1354760 x00009 + 1527205 x00010 + 1486247 x00011
 + 1941933 x00012 + 1634072 x00013 + 1084740 x00014 + 1350249 x00015
 + 1581194 x00016 + 1981871 x00017 + 1646604 x00018 + 1734106 x00019
 + 1042882 x00020 + 1763564 x00021 + 1397430 x00022 + 1177643 x00023
 = 17734962

Y: x00000 + x00001 + x00002 + x00003 + x00004 + x00005 + x00006 + x00007
 + x00008 + x00009 + x00010 + x00011 + x00012 + x00013 + x00014 + x00015
 + x00016 + x00017 + x00018 + x00019 + x00020 + x00021 + x00022 + x00023 = 1

VARIABLES
0 <= x00000 <= 1 Integer
0 <= x00001 <= 1 Integer
0 <= x00002 <= 1 Integer
0 <= x00003 <= 1 Integer
0 <= x00004 <= 1 Integer
0 <= x00005 <= 1 Integer
0 <= x00006 <= 1 Integer
0 <= x00007 <= 1 Integer
0 <= x00008 <= 1 Integer
0 <= x00009 <= 1 Integer
0 <= x00010 <= 1 Integer
0 <= x00011 <= 1 Integer
0 <= x00012 <= 1 Integer
0 <= x00013 <= 1 Integer
0 <= x00014 <= 1 Integer
0 <= x00015 <= 1 Integer
0 <= x00016 <= 1 Integer
0 <= x00017 <= 1 Integer
0 <= x00018 <= 1 Integer
0 <= x00019 <= 1 Integer
0 <= x00020 <= 1 Integer
0 <= x00021 <= 1 Integer
0 <= x00022 <= 1 Integer
0 <= x00023 <= 1 Integer

Отличаться они будут только значением правой части равенства в 18 строке.

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

Примечание

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

В моих тестах только CBC от COIN-OR Foundation и OR-Tools от Google, справились с данной задачей верно.

Заключение:

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

Спасибо за внимание!

Полный код

import pulp as pl
import random as rnd
from datetime import datetime
# ==========================================================================
def ssp(src: list, val: int):
    """
    Задача о сумме подмножеств натуральных чисел (Subset sum problem - SSP)
    """
    n = len(src)
    # Объявляем модель
    model = pl.LpProblem(name="ssp", sense=pl.LpMinimize)       
    # Подключаем решатель CBC
    solver = pl.PULP_CBC_CMD(msg=True)
    # Объявляем переменные
    x = [pl.LpVariable(name=f'x{i:05}', cat='Binary') for i in range(n)]    
    # Целевая функция
    model += pl.lpSum([x[i] for i in range(n)])    
    # Основное ограничение
    model += pl.lpSum([x[i] * src[i] for i in range(n)]) == val   

    # Находим решение
    status = model.solve(solver)
    if status == 1:
        s = []
        for j, v in enumerate(model.variables()):
            if round(v.value()) > 0:
                s.append(src[j])
        return s
    return []
# ==========================================================================
def ssp_optimized(src: list, val: int):
    """
    Задача о сумме подмножеств натуральных чисел оптимизированный алгоритм
    (Subset sum problem optimized - SSP)
    """
    n = len(src)
    # Объявляем модель
    model = pl.LpProblem(name="ssp", sense=pl.LpMinimize)    
    # Подключаем решатель CBC
    solver = pl.PULP_CBC_CMD(msg=True)
    # Объявляем переменные
    x = [pl.LpVariable(name=f'x{i:05}', cat='Binary') for i in range(n)]
    # Целевая функция
    model += pl.lpSum([x[i] for i in range(n)])
    # Основное ограничение
    model += pl.lpSum([x[i] * src[i] for i in range(n)]) == val

    loop = list(range(1, n + 1))
    loop = [loop.pop(0) if i % 2 == 0 else loop.pop(-1) for i in range(n)]
    for val in loop:
        # Устанавливаем дополнительное ограничение шага цикла
        model.constraints['Y'] = pl.lpSum([x[j] for j in range(n)]) == val 

        # Находим решение
        status = model.solve(solver)  
        if status == 1:
            s = []
            for j, v in enumerate(model.variables()):
                if round(v.value()) > 0:
                    s.append(src[j])
            return s
    return []
# ==========================================================================
n = 23
rnd.seed(0)
# Генерируем случайный ряд чисел
src = [rnd.randint(1000000, 2000000) for i in range(n)]
t = sum(src) // 2
print(f'source = {src}')
print()

start_time = datetime.now()
res = ssp(src, t)
date_diff = (datetime.now() - start_time).total_seconds()
print('solution =', res)
print('time =', date_diff)
print()

start_time = datetime.now()
res = ssp_optimized(src, t)
date_diff = (datetime.now() - start_time).total_seconds()
print('solution optimized =', res)
print('time =', date_diff)