0/1 Knapsack Problem

O problema da mochila é um velho conhecido dos que se preocupam com questões de optimização. Na sua formulação binária consiste em encontrar a melhor selecção de entre um conjunto de objectos, cada um com um peso e um valor, que uma mochila de capacidade limitada consegue suportar, procurando maximizar o valor total dos objectos seleccionados. Para os que têm uma inclinação matemática a definição rigorosa é a seguinte: $$ maximizar \quad\quad z = \sum_{j=1}^n v_j x_j$$ $$sujeito \quad a \quad\quad \sum_{j=1}^n p_j x_j \le c$$ $com \quad\quad x_j = 0 \quad ou \quad x_j = 1$ (vale 1 quando o item de ordem $j$ vai para a mochila). $v_j$ é o valor, $p_j$ é o peso, e $c$ a capacidade.

A dificuldade deste problema aumenta com o aumento da dimensão do conjunto dos objectos, como é natural, mas também é função da relação entre os pesos, os valores e a capacidade.

Vamos procurar resolver esta questão recorrendo a um dos métodos estocásticos de estado único dados nas aulas. Vamos ter que resolver várias questões, nomeadamente, a representação do problema, o modo como calculamos a qualidade das soluções candidatas, o algoritmo a usar e a forma como determinamos soluções alternativas candidatas a melhor solução.

Dada a natureza do problema uma representação binária é adequada: numa cadeia binária um 1 numa posição $j$ significa que o item $j$ está na mochila. Quanto à qualidade o mais simples é fazer a qualidade igual à soma dos valores de todos os objectos colocados na mochila. Vamos escolher como algoritmo o trepa-colinas básico. Já o candidato alternativo obtém-se trocando um bit uma posição escolhida aleatoriamente na cadeia candidata corrente. A qualidade de um indivíduo será dada simplesmente pelo valor total dos objectos seleccionados como decorre da definição do problema. Mas temos um ponto prévio a resolver: a escolha dos dados para testar a nossa solução. Os pesos e os valores podem ser não correlacionados, fracamente relacionados ou fortemente relacionados. A capacidade pode ser escolhida em função do peso total dos objectos, por exemplo, considerando metade desse valor. Podemos implementar as três situações de modo muito simples:

def generate_uncor(size_items,max_value):
    weights = [random.uniform(1,max_value) for i in range(size_items)]
    values = [random.uniform(1,max_value) for i in range(size_items)]
    capacity = int(0.5 * sum(weights))
    return weights, values, capacity

def generate_weak_cor(size_items,max_value, amplitude):
    """No negative weight allowed!"""
    weights = [random.uniform(1,max_value) for i in range(size_items)]
    values = []
    for i in range(size_items):
        value = weights[i] + random.uniform(-amplitude,amplitude)
        while value <= 0:
            value = weights[i] + random.uniform(-amplitude,amplitude)
        values.append(value)
    capacity = int(0.5 * sum(weights))
    return weights, values, capacity

def generate_strong_cor(size_items,max_value,amplitude):
    weights = [random.uniform(1,max_value) for i in range(size_items)]
    values = [weights[i] + amplitude for i in range(size_items)]
    capacity = int(0.5 * sum(weights))
    return weights, values, capacity

Como se pode ver, a geração dos dados depende de dois parâmetros: o valor máximo para os pesos, e o que designamos por amplitude (que determina o afastamento entre os pesos e os valores). Na listagem abaixo estes parâmetros estão definidos. Pode sempre testar com outros valores e observar a existência ou não de diferenças na qualidade dos resultados.

. Posto isto apresentamos a solução simples para o 0/1 KP.

import random
import matplotlib.pyplot as plt

# Data Sets
def generate_uncor(size_items,max_value):
    weights = [random.uniform(1,max_value) for i in range(size_items)]
    values = [random.uniform(1,max_value) for i in range(size_items)]
    capacity = int(0.5 * sum(weights))
    return weights, values, capacity

def generate_weak_cor(size_items,max_value, amplitude):
    """No negative weight allowed!"""
    weights = [random.uniform(1,max_value) for i in range(size_items)]
    values = []
    for i in range(size_items):
        value = weights[i] + random.uniform(-amplitude,amplitude)
        while value <= 0:
            value = weights[i] + random.uniform(-amplitude,amplitude)
        values.append(value)
    capacity = int(0.5 * sum(weights))
    return weights, values, capacity

def generate_strong_cor(size_items,max_value,amplitude):
    weights = [random.uniform(1,max_value) for i in range(size_items)]
    values = [weights[i] + amplitude for i in range(size_items)]
    capacity = int(0.5 * sum(weights))
    return weights, values, capacity


#  Hill Climbing  
#  Binary representation
#  Random neighbor

def kp_hc(max_iter,problem_size, weights,values,capacity):
    candidate = random_indiv(problem_size)
    cost_candi = fitness(candidate,weights,values,capacity)
    for i in range(max_iter):
        next_neighbor = random_neighbor(candidate)
        cost_next_neighbor = fitness(next_neighbor,weights,values,capacity)
        if cost_next_neighbor >= cost_candi: 
            candidate = next_neighbor
            cost_candi = cost_next_neighbor        
    return candidate,cost_candi/sum(values)

# ------------------------------------------------------------------------------
def random_indiv(size):
    return [random.randint(0,1) for i in range(size)]

 
def random_neighbor(individual):
    new_individual = individual[:]
    pos = random.randint(0,len(individual) - 1)
    gene = individual[pos]
    new_gene = (gene + 1) % 2
    new_individual[pos] = new_gene
    return new_individual


# Fitness for a binary representation. Unfeasible individuals have fitness 0

def fitness(individual, weights,values, capacity):
    items = [ i for i in range(len(individual)) if individual[i] == 1]
    total_weight = sum([ weights[i] for i in items])
    if total_weight > capacity:
        return 0
    else:
        total_value = sum([ values[i] for i in items])
        return total_value
 
    
if __name__ == '__main__':
    size_1 = 100
    size_2 = 250
    size_3 = 500
    max_value = 10
    amplitude = 5
    weights_nc, values_nc, capacity_nc = generate_uncor(size_1, max_value)
    weights_wc, values_wc, capacity_wc = generate_weak_cor(size_2, max_value,amplitude)
    weights_cor, values_cor, capacity_cor = generate_strong_cor(size_3, max_value, amplitude)
    
 
    print(kp_hc(1000,size_1,weights_nc,values_nc,capacity_nc))
    print(kp_hc(1000,size_2,weights_wc,values_wc,capacity_wc))
    print(kp_hc(1000,size_3,weights_cor,values_cor,capacity_cor))

Podemos explorar variantes? Claro. Por exemplo, visualizar os resultados.

def kp_hc_visual(max_iter,problem_size, weights,values,capacity):
    candidate = random_indiv(problem_size)
    cost_candi = fitness(candidate,weights,values,capacity)
    quality = [cost_candi]
    for i in range(max_iter):
        next_neighbor = random_neighbor(candidate)
        cost_next_neighbor = fitness(next_neighbor,weights,values,capacity)
        if cost_next_neighbor >= cost_candi: 
            candidate = next_neighbor
            cost_candi = cost_next_neighbor
        quality.append(cost_candi)
    display_kp(quality)
    return candidate,cost_candi/sum(values)

def display_kp(data):
    """Plot the data"""
    x = list(range(len(data)))
    plt.title('0/1 Knapsack')
    plt.xlabel('Generations')
    plt.ylabel('Value')
    plt.grid(True)
    plt.plot(x,data, 'r')
    plt.show()      

Podemos também explorar outras funções de qualidade:

def fitness_2(individual, weights,values, capacity):
    items = [ i for i in range(len(individual)) if individual[i] == 1]
    total_weight = sum([ weights[i] for i in items])
    total_value = sum([ values[i] for i in items])
    if total_weight > capacity:
        return total_value - (total_weight - capacity)
    else:
        return total_value  
    
def fitness_3(individual, weights,values, capacity):
    items = [ i for i in range(len(individual)) if individual[i] == 1]
    total_weight = sum([ weights[i] for i in items])
    total_value = sum([ values[i] for i in items])
    if total_weight == 0:
        return 0
    return total_value/sum(values)

A primeira alternativa, no caso de haver excesso de peso, penaliza em função da amplitude desse excesso. A segunda alternativa, usa como medida de qualidade o valor dos objectos que conseguimos colocar no saco em função do valor total dos objectos. É claro que agora precisamos de alterar ligeiramente o programa principal para que a função de qualidade passe a ser um parâmetro.

def kp_hc_visual_fit(max_iter,problem_size, weights,values,capacity,fitness):
    candidate = random_indiv(problem_size)
    cost_candi = fitness(candidate,weights,values,capacity)
    quality = [cost_candi]
    for i in range(max_iter):
        next_neighbor = random_neighbor(candidate)
        cost_next_neighbor = fitness(next_neighbor,weights,values,capacity)
        if cost_next_neighbor >= cost_candi: 
            candidate = next_neighbor
            cost_candi = cost_next_neighbor
        quality.append(cost_candi)
    display_kp(quality)
    return candidate,cost_candi/sum(values)

Agora só falta explorar o que foi apresentado. Pode também recorrer a outro algoritmo que não o trepa colinas.

Números de João Brandão

Existem problemas de enunciado simples mas que se revelam bem complexos quando se trata de os resolver. Um exemplo disso é o problema conhecido pelo nome de Os números de João Brandão. Reza assim. Consideremos o conjunto dos números inteiros entre 1 e um dado n. Formar o seu sub-conjunto de maior cardinalidade em que nenhum dos elementos é igual à média da soma de outros dois membros presentes no sub-conjunto. Por exemplo, para n=5 uma solução possível é {1,2,4,5}. É claro quanto maior for o valor de n mais complexa é a questão.

Qual será a melhor forma de o resolver? Um método simples é o de enumerar todos os sub-conjuntos de inteiros que é possível construir e seleccionar os que não violam a restrição indicada. Depois é fácil: basta escolher os de maior cardinalidade. Mas sabemos quantos sub-conjuntos podem existir para um dado n, para determinar se este método de força bruta método é exequível? Se pensarmos um pouco chegamos à conclusão que esse número é igual à soma de todas as combinações de n tomadas i a i, com i a variar entre 1 e n. Feitas as contas esse valor é igual a 2 levantado a n. Um número muito grande para valores de n razoáveis...

Um alternativa a esta questão da dimensão é usar um método heurístico, por exemplo, um dos métodos estocásticos de estado único dados nas aulas. Não temos a certeza de encontrar uma solução óptima, mas pelo menos teremos um valor ... antes de morrermos. Dito isto, uma primeira questão que se coloca é a de como representar o problema e de como definir a qualidade de uma solução. A primeira resolve-se facilmente: um vector binário, em que um 1 na posição i significa que o número i faz parte da solução. Quanto à qualidade temos duas coisas a considerar: (1) não existirem números que violem a restrição enunciada e, (2) ter o tamanho maior possível. Na solução apresentada iremos considerar estes dois aspectos. Vejamos como.

Dada uma solução na forma de um vector binário a qualidade será tanto melhor quanto maior for o número de números presentes, e isso é dado pela quantidade de uns no vector. Por outro lado, temos que calcular a quantidade de violações da restrição. Daí o código:

def evaluate(indiv,comp):
    alfa = 1
    beta = 1.5
    return alfa * sum(indiv) - beta * viola(fenotipo(indiv),comp)

Esta solução faz intervir dois parâmetros, alfa e beta, que pesam do ponto de vista relativo os dois aspectos (tamanho e violações). Notar que quanto maior for o valor achado, melhor!

Para calcular o número de violações vamos usar o conjunto dos inteiros correspondente à cadeia binária (função fenótipo):

def fenotipo(indiv):
    fen = [i+1 for i in range(len(indiv)) if indiv[i] == 1]
    return fen

O algoritmo que calcula o número de violações parte de uma ideia simples: um número está a mais se os inteiros à mesma distância dele também estiverem no conjunto. Por exemplo, se eu tiver o número 4, então a presença simultânea do par (3,4), (2,6) ou (1,7) provoca uma violação cada. Usando esta ideia chegamos ao código:

def viola(indiv,comp):
    # count violations of constraints
    v = 0
    for elem in indiv:
     limite = min(elem-1,comp - elem)
     vi = 0
     for j in range(1,limite+1):
      if ((elem - j) in indiv) and ((elem+j) in indiv):
       vi += 1
     v += vi
    return v

Resolvido esta primeira parte da busca de uma solução, a seguir escolhemos um algoritmo, por exemplo, o Trepa Colinas em que o vizinho que consideramos para comparar com a solução candidata corrente é o melhor vizinho deste último:

def jb_hc(problem_size, max_iter,fitness):
    candidate = random_indiv(problem_size,problem_size)
    cost_candi = fitness(candidate)
    for i in range(max_iter):
        next_neighbor = best_neighbor(candidate,fitness)
        cost_next_neighbor = fitness(next_neighbor,problem_size)
        if cost_next_neighbor >= cost_candi: 
            candidate = next_neighbor
            cost_candi = cost_next_neighbor        
    return candidate

Notar como é feita a comparação, pois, como dissemos, quanto maior o valor dado pela função que calcula o mérito melhor. E como calculamos o melhor vizinho? em problemas discretos em que o número de vizinhos é finito e pequeno a resposta é fácil:

 def best_neighbor(individual, fitness):
    best = individual[:]
    best[0] = (best[0] + 1) % 2
    for pos in range(1,len(individual)):
        new_individual = individual[:]
        new_individual[pos]= (individual[pos] + 1) % 2
        if fitness(new_individual) > fitness(best):
            best = new_individual
    return best

Na solução apresentada vamos fabricando de modo ordenado os vizinhos de uma solução candidata, considerando que estes são aqueles que diferem em apenas uma posição dessa solução. à medida que geramos uma nova solução vamos actualizando, se for caso disso, o melhor vizinho. Juntando tudo, aqui fica o código completo. Agora é só experimentar para diferentes valores de n. Uma nota: quando n é igual a 100, se conseguir uma solução de tamanho 27 ... diga-me!

import random

# Basic Hill Climbing

def jb_hc(problem_size, max_iter,fitness):
    candidate = random_indiv(problem_size,problem_size)
    cost_candi = fitness(candidate)
    for i in range(max_iter):
        next_neighbor = best_neighbor(candidate,fitness)
        cost_next_neighbor = fitness(next_neighbor,problem_size)
        if cost_next_neighbor >= cost_candi: 
            candidate = next_neighbor
            cost_candi = cost_next_neighbor        
    return candidate
 
    
# Random Individual
def random_indiv(size):
    return [random.randint(0,1) for i in range(size)]

# Best neighbor
def best_neighbor(individual, fitness):
    best = individual[:]
    best[0] = (best[0] + 1) % 2
    for pos in range(1,len(individual)):
        new_individual = individual[:]
        new_individual[pos]= (individual[pos] + 1) % 2
        if fitness(new_individual) > fitness(best):
            best = new_individual
    return best


# Fitness for JB
def evaluate(indiv,comp):
    alfa = 1
    beta = 1.5
    return alfa * sum(indiv) - beta * viola(fenotipo(indiv),comp)

def viola(indiv,comp):
    # count violations of constraints
    v = 0
    for elem in indiv:
     limite = min(elem-1,comp - elem)
     vi = 0
     for j in range(1,limite+1):
      if ((elem - j) in indiv) and ((elem+j) in indiv):
       vi += 1
     v += vi
    return v
    
# Auxiliar

def fenotipo(indiv):
    fen = [i+1 for i in range(len(indiv)) if indiv[i] == 1]
    return fen

    
if __name__ == '__main__':
 # For test purposes: beware of the time it may takes...
    res = fenotipo(jb_hc(100,10000,evaluate))
    quali = viola(res, len(res))
    print(quali,'\n', len(res), '\n',res)
 

É tempo de retomar o caminho

Este blog tem estado um pouco adormecido. Vamos retomar a escrita para acompanhar a edição de 2013-2014 da cadeira de Computação Evolucionária. Como no passado aqui irão aparecer ideias, acrescentos e soluções que o tempo for justificando.