Introdução a Constraint Programming usando Gecode

Agosto 19, 2012

Neste post vamos falar sobre constraint programming (CP) ou programação com restrições.

Constraint Programming é uma técnica em que se define um espaço de busca, um conjunto de restrições e deseja-se encontrar soluções viáveis, isto é, pontos que pertencem ao espaço de busca e que satisfazem as restrições.

Segundo a Wikipedia, surgiu a partir da programação por restrições lógicas (constraint logic programming). Essa variante de programação lógica foi criada por Jaffar e Lassez, que introduziram restrições ao Prolog II, sendo que Prolog III possui suporte nativo a programação por restrições lógicas.

É possível trabalhar com CP em linguagens imperativas através de frameworks. Algumas das bibliotecas mais famosas são Gecode (C++), IBM ILOG CP (C++ proprietário), Choco (Java), JaCoP (Java) e Google OR-Tools (C++, com bindings para Java, Python e C#). A biblioteca que utilizaremos para implementar as soluções para os exemplos é a Gecode (v 3.7.3).

Em geral os resolvedores de CP fazem uso de backtracking para encontrar essas soluções e portanto além da modelagem em si, outro fator que influencia o desempenho é a maneira como se faz as ramificações (branching) na árvore de busca.

Um breve tutorial de Gecode

Gecode é uma biblioteca open-source escrita em C++. É disponibilizada com licença MIT, é bastante eficiente, tem suporte a paralelismo e é bem documentada [1].

Modelos no Gecode devem definir um “espaço”. Dentro desse espaço ficam as variáveis, os propagadores (implementações de restrições) e branchers (implementações de branchings).

Um espaço é representado pela classe Space, que o modelo deve estender.

Vamos agora detalhar alguns conceitos utilizados pelo Gecode. Após isso, apresentaremos alguns exemplos clássicos que pode ser atacados com programação por restrição.

Variáveis

O Gecode define novos tipos para especificar o domínio das variáveis, como por exemplo IntVar, para variáveis inteiras e BoolVar para variáveis binárias. Também há o tipo representando um vetor desses tipos, IntVarArray e BoolVarArray.

Podemos construir um vetor de inteiros da seguinte maneira:

IntVarArray (Space &home, int n, int min, int max)

Sempre precisamos passar uma referência de um Space como parâmetro para informar em qual espaço criaremos a variável. O parâmetro n define o número de elementos no vetor, min é o limite inferior e max o liminte superior dos elementos nesse vetor.

Restrições

As restrições delimitam o domínio de busca de uma solução. No Gecode existem funções que criam propagadores implementando essas restrições. Algumas das principais são:

Restrições de relação:

void rel (Home home, IntVar x0, IntRelType r, IntVar x1);

Aqui, Home é equivalente a &Space, pois, novamente, precisamos especificar em que espaço a restrição será gerada. Essa restrição define uma relação entre a variável x0 e x1 usando uma relação r.

Um exemplo de uso é

rel(*this, x, IRT_NQ, 0)

Que modela a restrição x!= 0. Esse exemplo pode ser escrito de maneira mais legível com o uso de sobreposição de operadores:

rel(*this, x != 0)

Restrições de distinção:

Uma restrição bastante útil é a distinct, que exige que todos os valores de um vetor tenham valores diferentes:

void distinct (Home home, IntVarArray x0);

Restrições com expressões regulares:

É possível especificar que o vetor de solução satisfaça uma expressão regular através da restrição extensional.

extensional(Home home, const IntVarArgs &x, DFA d)

Essa restrição espera um vetor de inteiro e um DFA (Deterministic Finite Automata ou Autômato Determinístico Finito), que pode ser criado a partir de uma expressão regular. Essa função restringe os valores de x de forma que seja aceitável pelo DFA.

O tipo REG é uma maneira conveniente de se construir uma expressão regular a partir de expressões regulares menores.

Podemos criar a “base” de uma expressão regular com um único valor ou um conjunto possível de valores:

REG(int i);
REG (const IntArgs &x);

No exemplo abaixo, r corresponde a “1” e s corresponde a “(0|1)”.

REG r(1);

int v[] = {0, 1};
REG s(IntArray(vector<int>(v, v+2)));

Os operadores * e + foram convenientemente sobrescritos para modelar respectivamente “0 ou mais repetições” e “1 ou mais repetições”, por exemplo

REG r(1);
REG s = *r;
REG t = +r;

A variável s corresponde à expressão regular “1*” e t à “1+”.

O operador () também foi sobrescrito com duas versões: recebendo um inteiro n como parâmetro, que indica que a expressão da variável deve aparecer pelo menos n vezes; e a outra versão recebendo n e m, indicando que a expressão deve aparecer entre n e m (inclusive) vezes. Assim, no exemplo

REG r(1);
REG s = r(10);
REG t = r(3, 4);

A variável s corresponde à “1{10,}*” e t corresponde a “1{3,4}”.

Note que uma expressão regular seguida por “{n}” deve aparecer exatamente n vezes (e não pelo menos n vezes). Portanto, para modelar esse caso usando REG, devemos fazer

REG s = r(n, n);

Veremos a utilidade desse tipo de restrição no exemplo do Nonograma!

Cópia e clone

Como o algoritmo usa backtrack, ele precisa manter cópias para eventualmente voltar a uma solução anterior. Para isso precisamos implementar o método puramente virtual da classe Space, o copy(), retornando um ponteiro para a cópia criada.

virtual Space* copy(bool share) = 0;

A flag share está relacionada com o uso de múltiplas threads, mas não vamos nos preocupar com ela nesse post introdutório. Na cópia, o seguinte construtor da classe Space deve ser chamado:

Space::Space(boolean clone, Space &space);

Para fazer o clone de uma variável, podemos usar o método update(), presente na classe base dos tipos IntVar e IntVarArray:

template<class VarImp>
void Gecode::VarImpVar< VarImp >::update(Space &home, bool share, VarImpVar< VarImp > &y) 

Ele equivale a um clone. Existe um construtor de IntVar que recebe uma referência para outra IntVar, mas ele não cria um cópia e sim uma referência. Por exemplo,

IntVar x(y);

Aqui x e y referem-se à mesma variável.

Branching

Precisamos especificar o modo como será feito o backtracking através do método branch().

branch (Home home, const IntVarArgs &x, IntVarBranch vars, IntValBranch vals)

O primeiro parâmetro especifica o espaço onde a regra de branching será criada, x é um vetor de variáveis sobre o qual faremos branch, vars indica define como será feita a escolha das variáveis. O parâmetro vals indica como será feita a escolha dos valores.

Nos exemplos usaremos sempre a seguinte configuração:

branch(*this, sol, INT_VAR_SIZE_MIN, INT_VAL_MIN);

Aqui, sol representa a solução, INT_VAR_SIZE_MIN indica que o backtracking sempre fará branching sobre as variáveis com menor domínio primeiro e o parâmetro INT_VAL_MIN indica que os valores “chutados” serão feitos em ordem crescente.

Resolvedores

No Gecode há três tipos de motores de busca (search engines), um que usa busca em profundidade (chamado de DFS), outro que usa Branch-and-Bound (BAB) e um chamado Depth First Restart, que não sei exatamente o que é, mas talvez seja uma busca em profundidade iterativa (Iterative deepening depth search).

Os dois últimos motores são usados para quando se deseja otimizar uma função objetivo. Nos nosso exemplos estamos preocupados apenas em encontrar uma solução para o problema, então vamos usar apenas o DFS.

O resolvedores são um template para um tipo genérico e podem ser instanciados usando uma instância desse tipo. No nosso primeiro exemplo vamos implementar a classe SendMoreMoney e portanto podemos instanciar o motor de busca fazendo:

SendMoreMoney *m = new SendMoreMoney();
DFS<SendMoreMoney> e(m);

Os resolvedores definem o método next(), que retorna uma solução do problema (que é uma instância do modelo) e avança para a próxima. Se nosso modelo definir a função de imprimir uma solução, print(), podemos imprimir todas as soluções encontradas da seguinte maneira:

while(SendMoreMoney* s = e.next()){
    s->print();
    delete s;
}

A seguir, vamos comentar sobre 4 problemas clássicos e apresentar soluções usando CP. Os códigos-fontes estão disponíveis no Github e são adaptações dos códigos providos pela própria Gecode.

1. Puzzles Alpha-Numéricos

Um puzzle que pode ser resolvido com programação por restrições é aquele em que são dadas equações de palavras e cada letra da palavra deve ser substituída por um algarismo de modo que a equação seja satisfeita.

Letras iguais devem ser substituídas pelo mesmo algarismo, letras diferentes devem ter algarismos diferentes e a primeira letra de cada palavra não pode ser substituída por 0.

Um exemplo clássico é:

SEND + MORE = MONEY

Podemos resolver esse exemplo definindo 8 variáveis D, E, M, N, O, R, S, Y sujeitas às seguintes restrições:

1) As variáveis devem pertencer ao conjunto de inteiros [0-9].
2) M e S devem ser diferentes de 0.
3) As variáveis devem ter valores diferentes
4) A seguinte igualdade deve ser satisfeita

S*1000 + E*100 + N*10 + D +
M*1000 + O*100 + R*10 + E =
M*10000 + O*1000 + N*100 + E*10 + Y

Podemos implementar essas restrições da seguinte forma:

// Referencias para melhor legibilidade
IntVar 
s(sol[0]), e(sol[1]), n(sol[2]), d(sol[3]),
    m(sol[4]), o(sol[5]), r(sol[6]), y(sol[7]);

// Números nao podem comecar com 0
rel(*this, s != 0);
rel(*this, m != 0);

// Todas as variaveis devem ter valores distintos
distinct(*this, sol);

// Vetor de inteiros
IntArgs c(4+4+5);
// Vector de variaveis inteiras
IntVarArgs x(4+4+5);

// SEND
c[0]=1000; c[1]=100; c[2]=10; c[3]=1;
x[0]=s; x[1]=e; x[2]=n; x[3]=d;
// MORE
c[4]=1000; c[5]=100; c[6]=10; c[7]=1;
x[4]=m; x[5]=o; x[6]=r; x[7]=e;
// MONEY
c[8]=-10000; c[9]=-1000; c[10]=-100; c[11]=-10; c[12]=-1;
x[8]=m; x[9]=o; x[10]=n; x[11]=e; x[12]=y;

linear(*this, c, x, IRT_EQ, 0);

2. O problema das N-damas

O problema das N-damas (do inglês N-queens) é também bastante conhecido. Dado um tabuleiro de xadrez NxN, encontrar uma solução em que N damas sejam dispostas neste tabuleiro de forma que nenhuma esteja em posição de ataque com a outra. Segundo Hoffman, Loessi e Moore [2], sempre existe solução para N maior ou igual a 4.

Podemos definir a solução em outras palavras: queremos uma solução em que nenhum par de damas esteja na mesma coluna, na mesma linha ou na mesma diagonal. Claramente cada coluna deverá conter exatamente uma dama, de modo que podemos modelar a solução através de um vetor linhas de N posições, onde a linhas[i] representa a linha da rainha na i-ésima coluna.

Se colocarmos uma restrição de que os valores de linhas devam ser diferentes, garantimos que as damas estão em linhas diferentes. Entretanto, isso ainda não impede o posicionamento na mesma diagonal.

Se i e j estão na mesma diagonal principal, então dada a célula da dama i é dada por (linha[i], i) e a de j é (linha[i]+k, i+k). Se estão na diagonal invertida, então a posição de j é (linha[i]+k, i-k).

É fácil ver que i e j compartilham alguma diagonal se e somente se linha[i]-i = linha[j]-j ou se linha[i]+i = linha[j]+j.

Podemos implementar essas restrições da seguinte forma:

for (int i = 0; i<n; i++)
    for (int j = i+1; j<n; j++) {
        rel(*this, sol[i] != sol[j]);
        rel(*this, sol[i]+i != sol[j]+j);
        rel(*this, sol[i]-i != sol[j]-j);
    } 

3. Sudoku

No Sudoku clássico tem-se um tabuleiro 9×9, formado por 9 sub-tabuleiros 3×3. O objetivo é preencher o tabuleiro com algarismos de 1 a 9 de forma que nenhuma coluna, nenhuma linha e nenhum sub-tabuleiro contenha números repetidos. Além disso, na grande maioria das vezes os Sudokus já vêm com alguma de suas células preenchidas.

O Gecode oferece uma estrutura bastante simples de se trabalhar com matrizes, que é a classe Matrix. Temos os métodos row(i) e col(j) que retorna respectivamente a i-ésima linha e a j-ésima coluna.

Além disso, temos o método slice(i1, i2, j1, j2) que retorna a submatriz correspondendo a [i1, i2-1] e [j1, j2-1], o que é útil para modelar a restrição de sub-tabuleiros.

As restrições podem ser implementadas praticamente a partir da definição:

// Linhas e colunas devem ter valores diferentes
for (int i=0; i<N; i++) {
    distinct(*this, m.row(i));
    distinct(*this, m.col(i));
 }

// Restricoes dos sub-quadrados
for (int i=0; i<N; i+=n) {
    for (int j=0; j<N; j+=n) {
        distinct(*this, m.slice(i, i+n, j, j+n));
    }
 }

// Celulas ja preenchidas
for (int i=0; i<N; i++)
    for (int j=0; j<N; j++)
        if (mat[i][j] != '.')
            rel(*this, m(i,j), IRT_EQ, mat[i][j] - '0');

4. Nonograma

No nonograma é dado um tabuleiro com n linhas e m colunas. Para cada linha e coluna é dado um conjunto de números. Isso indica que naquela linha ou coluna deve blocos com exatamente aqueles comprimentos. Considere o exemplo abaixo:

Já tinha falado sobre esse puzzle anteriormente e à época não tinha conseguido pensar em uma solução para ele.

Vimos que o Gecode oferece uma maneira de especificar as restrições usando DFA’s. Se usarmos uma matriz de variáveis binárias como nosso vetor de soluções, então dada uma especificação de linha ou coluna dada por

3 4 5

Podemos criar a seguinte expressão regular:

0* 1{3} 0+ 1{4} 0+ 1{5} 0*

Basicamente, nossa linha ou coluna deve começar com 0 ou mais ‘0‘s, seguido por um bloco de exatamente 3 ‘1‘s. Entre dois blocos consecutivos de ‘1‘s deve haver pelo menos um ‘0‘.

Podemos transformar uma especificação em uma expressão regular da seguinte maneira:

// Expressao regular a partir de uma linha
DFA regexp(std::vector<int> &v) {
    REG r0(0), r1(1);
    REG border = *r0;
    REG between = +r0;
    REG r = border;
    for(size_t i = 0; i < v.size(); i++){
        if(i > 0) r += between;
        r += r1(v[i], v[i]);
    }
    return r + border;
}

Agora basta adicionar uma expressão regular dessas para cada coluna e linha. Nossa entrada é dada por um vetor de vetores com as especificações das linhas (_rows) e das colunas (_cols):

Matrix<BoolVarArray> m(sol, _cols.size(), _rows.size());
    
// Restricoes das colunas
for (size_t w=0; w < cols.size(); w++)
    extensional(*this, m.col(w), regexp(cols[w]));
// Restricao das linhas
for (size_t h=0; h < rows.size(); h++)
    extensional(*this, m.row(h), regexp(rows[h]));
// Branch   
for (int w = 0; w < cols.size(); w++)
    branch(*this, m.col(w), INT_VAR_NONE, INT_VAL_MIN);

Constraint Programming vs. Programação linear Inteira

Muitos dos problemas apresentados aqui também podem ser modelados como programação linear inteira e resolvidos através de um algoritmo de branch-and-bound.

Lembrando, no algoritmo de branch-and-bound resolve-se a relaxação linear do modelo em cada nó da árvore de busca. O valor dessa relaxação serve como limitante inferior (superior) para um problema de minimização (maximização). Com isso o algoritmo pode descartar soluções viáveis do problema sem analisá-las se o valor da relaxação em um nó for pior do que o valor da melhor solução encontrada até o momento.

Em CP não temos o valor desse limitante e portanto em um problema com função objetivo, em geral não se pode afirmar que uma solução é ótima a menos que se tenha analisado todas as soluções viáveis do problema.

Por outro lado, existem algumas modelagem que são complicadas de se fazer em programação linear inteira. A mais comum delas é a restrição de que duas variáveis devam ter valores diferentes.

Por exemplo no caso das N-damas, uma das restrições é de que cada elemento de linha tenha um valor diferente. Uma abordagem que pode ser usada para que variáveis x e y inteiras tenham valores diferentes é

(1) x >= y + 1

ou

(2) x <= y - 1

Lembrando que restrições disjuntas podem ser modeladas através da inclusão de uma variável binária, z que vale 0 se (1) for utilizada e 1 se (2) for utilizada. Dado uma constante M suficientemente grande temos:

(3) x \ge y + 1 - M \cdot z
(4) x \le y - 1 + M \cdot (1 - z)

Se z é 0, (3) vira (1) e (4) torna-se redundante. Se z é 1, (3) torna-se redundante e (4) vira (2).

Uma alternativa mais interessante é modelar os valores das variáveis inteiras como variáveis binárias. No caso do problema das N-damas, ao invés de trabalhar com linhas, poderíamos definir a variável c_{i,j} para cada célula do tabuleiro, sendo c_{i,j} = 1 se a dama está na posição (i,j) do tabuleiro e 0 caso contrário.

As restrições lineares a ser adicionadas são tais que a soma em cada linha e cada coluna seja exatamente 1 e que a soma em cada diagonal seja menor ou igual a 1.

Conclusão

Neste post fizemos uma pequena introdução a programação por restrições usando Geocode. Há muitos detalhes que faltam ser explorados nessa biblioteca como a customização de propagadores e branchers, bem como a comparação do desempenho de diferentes modelos.

Além disso, os exemplos empregados aqui são relativamente simples, com poucas restrições. A técnica de CP vem sendo aplicada a problemas do mundo real onde existem bastante detalhes que são modelados na forma de restrições. Pretendo estudar alguns papers e escrever sobre eles.

Ao decidir estudar CP, fiquei na dúvida sobre qual biblioteca utilizar. Eu pretendia usar a OR-Tools usando Python porque deve ser mais simples escrever os modelos. Entretanto, a documentação dessa bilioteca ainda está muito escassa e por isso acabei escolhendo a Gecode. De qualquer maneira Gecode também tem binding para uma linguagem simples de se escrever, que é o Ruby. Talvez seja uma boa oportunidade para aprender um pouco dessa linguagem.

Referências

[1] Gecode – GEneric COnstraint Development Environment
[2] Constructions for the Solution of the m Queens problem – E. Hoffmann , J. Loessi e R. Moore (Mathematics Magazine, Vol. 4, No. 2)

Anúncios

Mapas de Símbolos Fisicamente Realizáveis: a função Max-Min

Junho 5, 2012

Já escrevi sobre o tema de pesquisa do meu mestrado — otimização da visualização de mapas de símbolos proporcionais — em outros posts [1], [2] e [3].

Desta vez gostaria de falar sobre uma das variantes desse problema proposta por Cabello et al [4] onde o objetivo é basicamente dispor os símbolos no mapa de modo a maximizar a visibilidade do disco menos visível. Mais formalmente, dado um conjunto de discos S e uma disposição deles no plano, denotamos por b_i o comprimento da borda visível do disco i \in S.

O objetivo é encontrar uma disposição dos discos que maximize \min \{b_i\}. Não por acaso tal função é chamada de Max-Min.

Cabello et al. apresentaram um algoritmo guloso que maximiza essa função objetivo quando os discos só podem ser empilhados (não podem ser entrelaçados). A ideia é bem simples: dado um conjunto de n discos S, escolha o disco i^* que se posto por baixo de todos os outros teria o maior comprimento de borda visível (ou seja, b_i^* é máximo). Coloque i^* na base e determine a ordem dos outros discos recursivamente.

A versão “força-bruta” desse algoritmo é O(n^3), mas eles apresentam uma versão O(n^2 \log n) usando árvores de segmentos.

Se permitirmos que os discos possam ser entrelaçados, conforme a figura abaixo à esquerda, o problema se torna NP-difícil!

Disposição com entrelaçamento vs. disposição em pilha

Já havíamos atacado a versão que objetiva maximizar a soma das bordas visíveis considerando todos os discos (chamada de Max-Total) e permitindo entrelaçamento usando Programação Linear Inteira. Esse trabalho foi apresentado no Sibgrapi de 2011.

Após esse congresso, fomos convidados a estender nosso trabalho com pelo menos 30% de novo conteúdo e submeter para a publicação em um jornal científico, o The Visual Computer.

Decidimos então atacar a versão com função objetivo Max-Min usando PLI também. A formulação é bastante parecida com aquela que apresentamos no Sibgrapi, mas resolvedores de PLI têm bastante dificuldade em lidar com formulações do tipo max-min, onde se quer maximizar o menor valor de uma função.

Em nossos experimentos, instâncias para as quais o resolvedor achava a solução ótima da função Max-Total em segundos, não eram resolvidas mesmo rodando durante algumas horas para a versão Max-Min.

Uma possível causa para esse comportamento é que para formulações max-min temos muitas soluções com mesmo valor (platôs) e isso faz com que os otimizadores se percam.

Para um exemplo, considere as soluções na figura abaixo. Embora sejam bastante diferentes, elas possuem o mesmo valor da função Max-Min porque o disco destacado em vermelho é tão pequeno, que é ele que define o valor da função objetivo não importando a disposição entre o restante dos discos.

Soluções com mesmo valor da função Max-Min (definida pela borda do disco vermelho)

Como muitas das nossas instâncias podem ser decompostas em componentes menores, usamos a seguinte ideia para tentar melhorar o desempenho do nosso algoritmo. Considere as componentes decompostas c_1, c_2, \cdots, c_k. Seja s(c_i) o valor da solução ótima de c_i para Max-Min limitado a desenhos em pilha e p(c_i) permitindo também desenhos entrelaçados.

Suponha que as componentes estão numeradas de forma que s(c_1) \le s(c_2) \le \cdots \le s(c_k). Assim, a solução ótima da instância completa para desenhos em pilha é igual a s(c_1).

Temos também que s(c_i) \le p(c_i), pois estamos maximizando uma função e p considera tanto desenhos em pilha quanto entrelaçados.

Infelizmente não podemos afirmar que p(c_i) \le p(c_{i+1}). Porém, seja i^* = \mbox{argmin}_i\{p(c_i)\}. Se p(c_{i'}) \le s(c_{i'+1}) para algum i' \le k-1, então temos que p(c_{i'}) \le s(c_{j}) para i' < j e portanto p(c_{i'}) \le s(c_{j}) para i' < j, ou seja, i^* \le i'.

A ideia do algoritmo é resolver cada componente com o algoritmo polinomial de Cabello e ordená-las pelo valor da solução, obtendo s(c_1) \le s(c_2) \le \cdots \le s(c_k). Em seguida resolvemos as componentes nessa ordem usando a formulação do Max-Min permitindo desenhos entrelaçados, obtendo a cada passo p(c_i). Se em algum passo encontramos p(c_{i}) \le s(c_{i+1}), podemos parar pois já encontramos o valor ótimo.

Esse nosso trabalho foi aceito, tendo sido publicado online recentemente [5].

Referências

[1] Blog do Kunigami: Mapas de Símbolos Proporcionais
[2] Blog do Kunigami: Mapas de símbolos proporcionais e a meia integralidade da cobertura de vértices
[3]Blog do Kunigami: Visualização ótima de desenhos fisicamente realizáveis
[4] S. Cabello, H. Haverkort, M. van Kreveld, and B. Speckmann, “Algorithmic aspects of proportional symbol maps”, Algorithmica, 2010.
[5] G. Kunigami, P. de Rezende, C. de Souza and T. Yunes, “Generating Optimal Drawings of Physically Realizable Symbol Maps with Integer Programming”, 2012


Relaxação Lagrangeana – Prática

Março 11, 2012

Em um post anterior, expliquei brevemente o conceito de relaxações, focando em relaxação lagrangeana. Hoje vamos apresentar uma heurística que encontra multiplicadores lagrangeanos.

Aplicaremos o algoritmo no problema da localização de instalações.

Algoritmo

O algoritmo mais conhecido para encontrar bons multiplicadores lagrangeanos é dado por Beasley [1], e é conhecido por método do subgradiente.

Geometricamente, se tivéssemos dois multiplicadores lagrangeanos, estes representariam pontos em um plano e os valores das soluções ótimas do dual lagrangeano para um desses pontos (multiplicadores) formariam curvas de níveis. O que o algoritmo faz é procurar um ponto com o valor mais próximo do da solução ótima primal e para isso anda em uma direção (gradiente) em que há melhoria dessa função.

Se não houver melhoras por um longo tempo, pode ser que o “passo” que estamos dando a cada iteração esteja muito grande e pode ser que haja uma solução melhor do meio desse passo. Por isso, vamos diminuímos gradativamente o passo, até que eventualmente este fique menor que um tamanho mínimo.

Onde a função Atualiza_multiplicadores é dada por:

Neste algoritmo é usado o valor de uma solução primal, que é obtida a partir da relaxação através de uma heurística.

Aplicação ao problema da Localização de Instalações

Vamos usar a relaxação das desigualdades que forçam os clientes a serem cobertos por pelo menos uma instalação (como descrito em um post anterior).

Heurística lagrangeana. Como estamos trabalhando com a versão sem restrição de capacidade, dado um conjunto de fábricas abertas é trivial encontrar a solução de melhor custo, bastando ligar cada cliente à instalação aberta mais barata. É essa heurística que usaremos para obter uma solução primal a partir de uma solução relaxada.

Multiplicadores lagrangeanos. Note que o algoritmo requer um conjunto de multiplicadores lagrangeanos inicial. No nosso caso, fizemos u_i = \frac{1}{n}, onde n é o número de desigualdades relaxadas.

Para os outros parâmetros do algoritmo usamos: \pi = 2, T_\pi = 0.2 e N_{stuck} = 30.

Fiz uma implementação em C++, tentado separar o problema (classe FacilityLocation implementando a classe virtual LagrangeanModel) do algoritmo (classe GradientMethod). Como de costume, o código está no github.

Resultados computacionais

Para testar o algoritmo, usei as instâncias do problema Capacitated Facility Location da OR-Library [3]. Como estamos lidando com a versão não restrita (uncapacitated), ignorei os dados de capacidade.

Há também este arquivo com os valores ótimos das soluções de forma que podemos observar a qualidade das soluções obtidas via relaxação empiricamente. A tabela a seguir contém o nome da instância, o número de instalações, o número de clientes, a distância percentual entre o ótimo e o melhor valor primal obtido (Gap), a distância percentual entre os melhores valores valores primal e dual (Gap DP) e finalmente o tempo.

Observamos que a solução ótima foi encontrada para boa parte das instâncias. O pior resultado foi para a instância capb, devido à dificuldade do algoritmo em encontrar bons multiplicadores lagrangeanos (note a distância entre os valores primal e dual).

Note que mesmo para instâncias grandes, provavelmente inviáveis de serem resolvidas de maneira exata, são resolvidas em menos de 6 segundos (processador i7 2.2GHz, sem paralelismo).

Conclusão

Embora tenhamos que lidar com uma formulação do problema para decidir quais restrições relaxar e obter o problema relaxado, geralmente para resolvê-lo não precisamos trabalhar com a formulação explicitamente (em geral resolver PL’s, principalmente de formulações grandes, não é muito rápido). Assim, dependendo da complexidade do problema resultante, a relaxação lagrangeana pode ser bastante eficiente.

Vimos que na prática as soluções obtidas via heurística lagrangeana podem ser muito boas. O principal problema é que a escolha dos parâmetros podem influenciar bastante a qualidade dessas soluções.

Eu já havia feito um trabalho sobre relaxação lagrangeana durante uma disciplina na faculdade, sobre um problema chamado Axis Parallel Minimum Stabbing Number Spanning Tree [2].

Referências

[1] J. E. Beasley – Modern Heuristics for Combinatorial Optmization Problem, Chapter 6 (Lagrangean Relaxation). Blackwell Scientific Publications, 1993.
[2] Relatório Programação Linear Inteira (MO420) – Guilherme Kunigami
[3] OR-Library


Relaxação Lagrangeana – Teoria

Fevereiro 5, 2012

Neste post vamos comentar brevemente sobre relaxações de formulações lineares inteiras e focar especialmente em relaxações lagrangeanas. Depois, vamos demonstrar a aplicação dessa técnica no problema da localização de instalações (facility location). A menos que dito o contrário, estamos falando de problemas de minimização.

Relaxações

Em programação linear inteira, usamos o termo relaxação para indicar que algumas das restrições do problema original foram relaxadas. Provavelmente o exemplo mais conhecido de relaxação é o da relaxação linear. Nesse caso em particular, removemos as restrições de integralidade do modelo.

Mais formalmente, se tivermos uma formulação F representada por

z = \min \{c(x): x \in P\}

sendo P o conjunto de pontos que satisfazem as restrições dessa formulação, c(x) a função objetivo e z o valor da solução ótima. Dizemos que uma formulação F_R, dada por

z_R = \min \{f(x): x \in R\}

é uma relaxação se todo ponto de P está em R, ou seja P \subseteq R e se f(x) \le c(x) para x \in P.

Note que devido a essas propriedades, a solução ótima da relaxação é um limitante dual (ou seja, no nosso caso é um limite inferior e superior no de maximização) para a formulação original ou seja, z_R \le z. Limitantes duais podem ser utilizados, por exemplo, para melhorar o desempenho do algoritmo de Branch and Bound.

A vantagem de se usar relaxações é que em certos casos elas são mais fáceis de serem resolvidas. Por exemplo, as relaxações lineares podem ser resolvidas em tempo polinomial, enquanto a versão restrita a inteiros, em geral, é NP-difícil.

Podemos também relaxar uma formulação removendo algumas de suas restrições. Note, entretanto, que quanto mais relaxamos a formulação pior tende a ser a qualidade do limitante dual (considere o caso extremo onde removemos todas as desigualdades do modelo: não deixa de ser uma relaxação, mas o limitante dual obtido não serve para nada).

O ideal é que encontremos uma relaxação fácil de resolver e que dê limitantes duais de qualidade, ou seja, bem próximos do valor ótimo da formulação original.

Relaxação Lagrangeana

Existe uma técnica denominada relaxação lagrangeana [1], que consiste em remover algumas das restrições da formulação original, mas tenta embutir essas desigualdades na função objetivo. A ideia é penalizar a função objetivo quando as restrições removidas forem violadas. O “peso” dessas penalidades é controlado por coeficientes chamados multiplicadores lagrangeanos.

Em termos de matemáticos, suponha que tenhamos uma formulação com m restrições e n variáveis:

z^* = \min cx

Sujeito a

\begin{array}{llcl}    & A_1x & \le & b_1 \\   & A_2x & \le & b_2 \\   & x & \ge & 0 \\   & x & \in & Z^n  \end{array}

Onde A_1x \le b_1 representa as desigualdades que queremos relaxar e A_2x \le b_2 as desigualdades que vamos manter. A relaxação lagrangeana é então dada por

z(u) = \min cx - u(b_1 - A_1x)

Sujeito a

\begin{array}{llcll}   & A_2x & \le & b_2 & \\  & x & \ge & 0 & \\  & x & \in & Z^n & \\  & u & \ge & 0, & u \in R^m_1  \end{array}

Onde u é um vetor representando os multiplicadores lagrangeanos (um por cada desigualdade removida). Por exemplo, considere a desigualdade j que será removida:

\sum_{i = 1}^n a_{ij} x_i \le b_j

Ao remover (11), o seguinte fator é subtraído da função objetivo

(1) u_j (b_j - \sum_{i = 1}^n a_{ij} x_i)

Note que quanto mais “violada” a desigualdade estiver, mais negativo o termo (1) vai ficar. Como a função é de minimização e estamos subtraindo, intuitivamente a solução ótima para o problema tenderá a não violar as desigualdades, melhorando o resultado do limitante dual.

Escolhe-se remover as desigualdades de forma que o problema resultante seja mais fácil resolver para um valor fixo de multiplicadores lagrangeanos. Porém, agora temos que resolver um outro problema: encontrar um conjunto de multiplicadores lagrangeanos u^* que minimize z(u).

Aplicação ao problema da localização de instalações

O problema da localização de instalações é bastante estudado na literatura. Uma de minhas iniciações científicas consistia em estudar algoritmos aproximados para esse problema.

Basicamente temos um grafo bipartido onde uma partição corresponde ao conjunto de clientes C e a outra ao conjunto de fábricas/instalações F. Cada cliente i pode ser potencialmente atendido pela fábrica j, a custo um c_{i,j}. Além disso, paga-se um custo c_j por abrir a fábrica j, dado por y_j. O objetivo é atender todos os clientes com o menor custo.

Vamos apresentar uma formulação linear inteira para esse problema. Seja x_{i,j} uma variável binária que vale 1 sse o cliente i é atendido pela fábrica j e a variável y_j uma variável binária que vale 1 sse a fábrica j foi aberta. Temos a seguinte formulação:

z = \min \sum_{i \in C, j \in F} x_{ij} c_{ij} + \sum_{j \in F} y_j c_j

Sujeito a

\begin{array}{llcll}   (2) & \sum_{j \in F} x_{ij} & \ge & 1 & \forall \, i \in C \\  (3) & x_{ij} & \le & y_j & \forall i \in C, \forall j \in F \\  & x_{i,j}, y_{j} & \in & \{0, 1\} & \forall i \in C, \forall j \in F  \end{array}

A desigualdade (2) diz que cada cliente deve ser atendido por pelo menos uma fábrica e a desigualdade (3), que uma fábrica deve estar aberta para que possa atender a algum cliente.

Essa é uma formulação bem simples e temos basicamente dois grupos de desigualdades. Vamos tentar relaxar o grupo de desigualdades indicado por (2). Teremos:

z_1(u) = \min \sum_{i \in C, j \in F} x_{ij} c_{ij} +
\sum_{j \in F} y_j c_j - \sum_{i \in C} u_i (\sum_{j \in F} x_{ij} - 1)

Sujeito a

\begin{array}{llcll}   & x_{ij} & \le & y_j & \forall i \in C, \forall j \in F \\  & u_i & \ge & 0 & \forall i \in C\\   & x_{i,j}, y_{j} & \in & \{0, 1\} & \forall i \in C, \forall j \in F  \end{array}

Com um pouco de álgebra, é possível “fatorar” o somatório da função objetivo:

z_1(u) = \min \sum_{j \in F} (\sum_{i \in C} (x_{ij} (c_{ij} - u_{i})) + c_j y_j)) + \sum_{i \in C} u_i

Como o termo \sum_{i \in C} u_i é constante, podemos resolver o problema sem ele e depois somá-lo ao valor da solução obtida. Além disso, note que cada fábrica contribui de maneira independente à função objetivo e em cada restrição só aparecem variáveis associadas a uma dada fábrica.

Desta forma, podemos decompor essa formulação por cada fábrica j:

\min \sum_{i \in C} (x_{ij} (c_{ij} - u_{i})) + c_j y_j

Sujeito a:

\begin{array}{llcll}   & x_{ij} & \le & y_j & \forall i \in C\\  & x_{i,j}, y_{j} & \in & \{0, 1\} & \forall i \in C  \end{array}

Cada um desses problemas pode ser resolvido por inspeção. Seja c'_{ij} = c_{ij} - u_{i}. Vamos considerar duas possibilidades. Se fizermos y_j = 0, não atenderemos nenhum cliente, levando a uma solução de custo 0. Se y_j = 1, só compensa atender aos clientes com c'_{ij} < 0. Assim, o custo ótimo dessa formulação é:

\min\{0, c_j + \sum_{c'_{ij} < 0} c'_{ij} \}

A outra alternativa consiste em relaxar a desigualdade (3). Nesse caso, depois de alguma álgebra, teremos:

z_2(u) = \min \sum_{i \in C, j \in F} x_{ij} (c_{ij} - u_{ij}) + \sum_{j \in F} y_j (c_j + (\sum_{i \in C} u_{ij}))

Sujeito a:

\begin{array}{llcll}   & \sum_{j \in F} x_{ij} & \ge & 1 & \forall \, i \in C \\  & u_{ij} & \ge & 0 & \forall i \in C, j \in F\\  & x_{i,j}, y_{j} & \in & \{0, 1\} & \forall i \in C, \forall j \in F  \end{array}

Note que como não há restrições misturando x e y, podemos resolver o modelo para cada tipo de variável separadamente. No caso da variável y, como c_j \ge 0 e u \ge 0, não compensa abrir fábrica nenhuma. No caso de x, podemos resolver de maneira independente para cada x_i, por inspeção. Seja c'_{ij} = c_{ij} - u_{ij}. Se c'_{ij} > 0, podemos fazer x_{ij} = 1. Caso contrário x_{ij} = 0, exceto no caso em que c'_{ij} \ge 0 para todo j \in F. Nesse caso, devido à restrição (25), temos que fazer x_{ij^*} = 1, onde j^* = \mbox{argmin}_{j \in F} \{c'_{ij} \}.

Temos duas formulações que podem ser facilmente resolvidas. Entretanto, devemos levar em conta qual possui a formulação mais forte. Em outras palavras, dado z_1^* = \max \{ z_1(u) : u \ge 0 \} e z_2^* = \max \{ z_2(u) : u \ge 0 \}, queremos o limitante dual que mais se aproxima de z. Entretanto, não sei dizer qual das duas relaxações apresentadas é a mais forte.

Conclusão

A relaxação Lagrangeana é uma técnica poderosa para se obter limitantes duais de problemas combinatórios que podem ser modelados como programas lineares inteiros. Esses limitantes podem ser usados na tentativa de melhorar o desempenho de algoritmos branch-and-bound.

Além do mais, para certos problemas, existe a possibilidade de ajustar a solução obtida via relaxação de modo a obter uma solução viável para o problema original, esperando que a qualidade da solução não seja muito pior.

Esse post apresentou a teoria por trás da relaxação lagrangeana. Num próximo post vou apresentar uma heurística que visa encontrar bons multiplicadores lagrangeanos.

Referências

[1] The Lagrangian relaxation method for solving integer programming problems – M.L. Fisher Management Science, 1981.


Linearização de desigualdades não-lineares

Dezembro 18, 2011

Há algum tempo atrás, assisti uma apresentação em que o palestrante chegou a uma equação não linear e disse que era possível linearizá-la, entretanto sem dizer como, apenas passando uma referência.

Não consegui encontrar a referência na época, mas recentemente alguém postou no OR Exchange (sistema semelhante ao StackOverflow, mas voltado à Pesquisa Operacional) uma dúvida relacionada. Uma das respostas passou uma referência para um documento do AIMMS chamado Integer Programming Tricks [1].

Nesse documento há algumas técnicas de modelagem de programação linear inteira. Algumas delas eu já conhecia como por exemplo modelar restrições disjuntas (quando exatamente uma de duas restrições deve ser satisfeita).

A novidade é a técnica usada para transformar um produto de duas variáveis, por exemplo x_1 \cdot  x_2, em uma equação linear. Há três casos que devemos considerar: quando ambas variáveis são binárias; quando uma das variáveis é binária e a outra é contínua; e quando as duas variáveis são contínuas.

Caso 1: Ambas variáveis binárias

No caso em que ambas x_1 e x_2 são binárias, a modelagem é fácil:

(1) \, y \le x_1

(2) \, y \le x_2

(3) \, y \ge x_1 + x_2 - 1

(4) \, y \in \{0, 1\}

Fazendo as 4 possíveis combinações de valores, é fácil ver que y = x_1 \cdot x_2.

Caso 2: Uma das variáveis é contínua

A formulação nesse caso não é muito diferente. Vamos supor que é x_2 a variável contínua e é limitada por 0 \le x_2 \le u. Temos:

(5) \, y \le u \cdot x_1

(6) \, y \le x_2

(7) \, y \ge x_2 - u(1 - x_1)

(8) \, y \ge 0

Se x_1 = 0, temos que y = 0. Se x_1 = 1, então (7) vira y \ge x_2, que combinado com (6), implica y = x_2.

Caso 3: Ambas as variáveis contínuas

Esse caso é bem mais complicado e requer a introdução de um novo conceito: formulações lineares por partes (numa tradução livre de piecewise linear formulation).

Formulação linear por partes

Essa técnica consiste em aproximar uma dada função/restrição não-linear por uma função linear por partes, e é bastante comum em métodos numéricos.

Uma função é linear por partes (piecewise linear function) se pudermos subdividir o domínio em intervalos nos quais a função é linear.

x

Função não linear aproximada por uma função linear por partes

As funções que podem ser aproximadas por essa técnica devem ser separáveis. Uma função é separável se ela pode ser escrita como a soma de funções que só dependem de uma variável. Por exemplo, x_1^3 + x_2^2 é separável, mas x_1 \cdot x_2 não é.

Assim, dada uma função separável, considere uma função g do somatório. Suponha que aproximamos g usando uma função linear por partes f, com intervalos definidos pelos pontos x_1, x_2, \cdots, x_n, sendo que a função f entre x_i e x_{i+1} é linear.

Podemos modelar a função f(x) através de programação linear. Para isso, considere as variáveis contínuas \lambda_i, tal que 0 \le \lambda_i \le 1. Temos então

(9) \, \sum_{i=1}^n \lambda_i = 1

(10) \, f(x) = \sum_{i=1}^{n} f(x_i) \lambda_i

(11) \, x = \sum_{i=1}^{n} x_i \lambda_i

Com a restrição adicional de que no máximo dois \lambda‘s sejam não nulos e no caso de serem dois, esses \lambda‘s devem ser consecutivos. A ideia é que se só um \lambda_i for não nulo, ele representa x_i e f(x_i). No caso de \lambda_i e \lambda_{i+1} serem não nulos, eles representam x' = x_i \lambda_i + x_{i+1} \lambda_{i+1} e f(x') =  f(x_i) \lambda_i + f(x_{i+1}) \lambda_{i+1}.

Note que como \lambda_i + \lambda_{i+1} = 1, [x', f(x')] é uma combinação convexa de [x_i, f(x_i)] e [x_{i+1}, f(x_{i+1})] e portanto pertence à reta que une esses dois pontos.

Essa restrição pode ser embutida no algoritmo do simplex. Para isso, os resolvedores de programação linear inteira usam o SOS2 (Special Order Sets 2), no qual devemos especificar o conjunto de variáveis x_i e f(x_i) (os SOS2 não tem nada a ver com SOS1, sobre o qual falei num post anterior).

Entretanto, alguns resolvedores como o GLPK não possuem essa funcionalidade embutida. Em [2], eles apresentam uma formulação de PLI para modelar essa restrição.

São introduzidas novas variáveis binárias, z_i para i = 1, \cdots, n-1. Temos que z_i = 1 se e somente se, o valor de x está entre x_i e x_{i+1}. Também introduzimos uma variável contínua s_i para i = 1, \cdots, n-1, onde 0 \le s_i \le 1. Nesse caso s_i representa a combinação convexa entre x_{i} e x_{i+1}.

Definindo y_i = f(x_i) para i = 1, \cdots, n e y = f(x), temos o seguinte conjunto de restrições:

(12) \, \sum_{i = 1}^{n-1} z_i = 1

(13) \, s_i \le z_i

(14) \, x = \sum_{i = 1}^{n-1} x_i z_i + (x_{i+1} - x_{i}) s_i

(15) \, y = \sum_{i = 1}^{n-1} y_i z_i + (y_{i+1} - y_{i}) s_i

Como z_i‘s são binárias, exatamente um z_j será igual a 1 (e o resto zero). Nesse caso, (14) reduzirá para:

x = x_j + (x_{j + 1} - x_{j}) s_j,

que podemos reescrever como

x = x_j (1 - s_j) + x_{j+1} s_j

Para a equação (15) teremos, de maneira análoga:

y = y_j (1 - s_j) + y_{j+1} s_j

Ficando mais evidente ver que [x, y] é combinação convexa de [x_j, y_j] e [x_{j+1}, y_{j+1}].

Eliminação do produto de variáveis diferentes

Agora que vimos como fazer a formulação linear por partes, vamos ver como linearizar x_1 \cdot x_2 supondo que ambas são contínuas.

Para esse fim, definimos duas variáveis y_1 = \frac{1}{2} (x_1 + x_2) e y_2 = \frac{1}{2} (x_1 - x_2). Agora é fácil ver que y_1^2 - y_2^2 = x_1 \cdot x_2, e é uma função separável.

Agora podemos usar a formulação apresentada anteriormente da seguinte forma: aproximar y_1^2 usando k pontos (note o tradeoff entre precisão e tempo de execução), gerando pares [{y_1}_i, {y_1}_i^2] para i = 1, \cdots, k. Agora basta substituir ocorrência de y_1 pelo ‘x’ de (14) e y_1^2 pelo ‘y’ da equação (15). Repetimos o mesmo procedimento para a função -y_2^2.

Pronto! Agora todas as restrições do modelo são lineares!

Conclusão

Aprendi uma técnica muito útil para atacar alguns problemas com funções/restrições não-lineares. Ainda pretendo aprender programação não-linear, quer cursando uma disciplina, quer estudando um livro por conta.

Fiquei na dúvida sobre como aproximar uma função com domínio muito grande. Podem ser necessários muitos intervalos para uma boa aproximação e o número de variáveis cresce proporcionalmente.

Outra pergunta que não sei responder é o que dá para fazer quando ambas as variáveis são inteiras. É possível transformar uma variável inteira em binária, mas novamente, dependendo do tamanho do domínio das variáveis inteiras, a formulação pode ficar com um tamanho demasiadamente grande. Será que existe um jeito mais eficiente de linearizar x_1 \cdot x_2 nesse caso?

Só descobri esse material graças ao tópico no OR Exchange. Por sua vez, só fiquei sabendo do tópico graças ao twitter do OR Exchange. Aliás, tenho conseguido usar o twitter como ferramenta para obtenção de informações (úteis!). A maior parte das contas que passei a seguir ultimamente são agregadores de notícias/links relacionados à pesquisa operacional.

Referências

[1] AIMMS Modelling Guide – Integer Programming Tricks
[2] SOS2 constraints in GLPK


Visualização ótima de desenhos fisicamente realizáveis

Outubro 23, 2011

Hoje vou descrever o trabalho que apresentei no Sibgrapi 2011, entitulado “Determining an Optimal Visualization of Physically Realizable Symbol Maps“, escrito em conjunto com meus orientadores (professores Pedro e Cid) e o professor Tallys da Universidade de Miami. Há um

Trata-se de uma extensão do estudo que realizamos para o problema de mapas de símbolos proporcionais usando desenhos em pilha, sobre o qual comentei em um post anterior.

Desenho fisicamente realizável

Um desenho fisicamente realizável é uma generalização de um desenho em pilha, pois este permite ordens cíclicas, como no exemplo abaixo.


Desenho fisicamente realizável e desenho em pilha

Podemos pensar que desenhos fisicamente realizáveis são aqueles feitos a partir de discos de papel, desde que estes não sejam dobrados nem cortados.

O problema é encontrar um desenho fisicamente realizável que maximize a soma dos comprimentos das bordas visíveis. Esse problema é chamado Max-Total para desenhos fisicamente realizáveis.

Tínhamos desenvolvido um modelo de programação linear inteira para o problema Max-Total para desenhos em pilha e mostramos que um subconjunto das desigualdades desse modelo, resolvia o problema Max-Total para desenhos fisicamente realizáveis.

Discutimos aspectos teóricos sobre as desigualdades utilizadas no modelo, argumentando que elas tornam a formulação apertada. Na prática um modelo com formulação apertada tende a ser resolvido mais rapidamente pelo algoritmo de branch and bound.

Técnica de decomposição

Desenvolvemos uma técnica de decomposição nova, que só serve para esse tipo de desenho. A ideia básica das decomposições é quebrar o conjunto de discos de entrada em componentes menores e mostrar que podemos resolver a cada componente de maneira independente e então combiná-las em tempo polinomial para gerar a solução ótima da instância completa. Uma decomposição óbvia é considerar conjuntos de discos disjuntos no plano.

Para desenhos fisicamente realizáveis, podemos ignorar a interseção de discos em alguns casos. Isso por sua vez pode gerar novas componentes disjuntas. No exemplo abaixo, mostramos que é possível ignorar a interseção dos discos nas faces amarelas. Isso quebra a componente em duas que podem ser resolvidas isoladamente.

Exemplo de decomposição

Junto com as outras decomposições que já havíamos desenvolvido para desenhos em pilha, conseguimos diminuir o tamanho das instâncias de teste.

Resultados computacionais

Após resolvermos a formulação do problema para desenhos fisicamente realizáveis, constatamos que na grande maioria dos casos, o valor da solução era exatamente igual ao da solução para desenhos em pilha.

Mesmo para as instâncias onde houve diferença, esse valor foi muito pequeno e visualmente é muito difícil distinguir as soluções, conforme pode ser visto na figura a seguir.

Zoom de uma solução ótima usando desenho em pilha e desenho fisicamente realizável

Entretanto, também tivemos uma boa surpresa, que foi o tempo de execução desse novo modelo. Em média, ele foi aproximadamente 2 ordens de magnitude mais rápido do que modelo para desenhos fisicamente realizáveis.

Conclusão

Devido aos tempos obtidos com o novo modelo e o fato de as soluções obtidas serem muito próximas a desenhos em pilha, uma ideia é usarmos esse modelo para resolver o problema de desenhos em pilha, adicionando restrições para remover os eventuais ciclos que apareçam.

Podemos adicionar essas restrições de remoção de ciclo através de um algoritmo de planos de corte, de um modo similar ao usado para resolver o problema do caixeiro viajante.


Mapas de símbolos proporcionais e a meia integralidade da cobertura de vértices

Julho 31, 2011

O problema que ataquei no mestrado consistia em achar uma ordem de empilhamento entre símbolos em um mapa, de modo a maximizar uma dada função objetivo. Não se sabe a complexidade desse problema e por isso usei programação linear inteira.

Há uma variante desse problema, que consiste em maximizar uma outra função objetivo, para a qual existe um algoritmo polinomial.

Depois de obter as soluções para o meu problema, precisei fazer uma comparação entre as soluções desses dois problemas. Aconteceu que a diferença visual entre as duas soluções não é tão perceptível, sendo que apenas alguns discos estão em posição diferente. Para se ter uma idea, veja um dos casos


Jogo dos 7 erros ;) (clique na imagem para aumentar)

Decidi fazer um programa para destacar apenas os discos que mudaram de lugar. A minha primeira ideia foi marcar todos os discos que estavam em níveis diferentes entre as duas soluções.

Isso não funciona bem quando temos duas componentes disjuntas de discos, por exemplo na figura abaixo, com duas componentes disjuntas {A, B, C} e {D, E}. Note que as ordens {A > B > C > D > E} e {D > E > A > B > C} resultam visualmente na mesma figura, embora nenhum disco esteja exatamente na mesma posição na ordenação e portanto o programa vai marcar todos os discos.

Desenho que pode ser representado por mais de uma ordem dos discos

Uma segunda ideia foi marcar, considerando todos os pares de discos que se interceptam, apenas aqueles em que a ordem relativa do par foi invertida. Com essa estratégia, nenhum disco do exemplo acima seria destacado, o que está de acordo com nossas expectativas.

Porém, essa abordagem ainda apresenta um efeito indesejável. Se tivermos uma pilha de discos, onde todos se interceptam e um disco que está no topo em uma solução aparece na base na outra, sendo que todos os outros discos permaneceram intactos, como por exemplo {A > B > C > D} e {B > C > D > A}, ainda assim todos os discos serão destacados pois todos os pares formados com A tiveram sua ordem relativa invertida. Note que seria suficientemente informativo destacar apenas o disco A.

Assim, para todo par de discos que se intercepta, basta destacarmos um deles. Como então destacar o menor número possível para diminuir a poluição visual?

Se olharmos para o grafo de discos desse problema, existem algumas arestas que queremos destacar (cobrir), usando pelo menos um dos vértices mas tentando minimizar a quantidade desses vértices. Soa familiar?

Sim, reduzimos para o problema da cobertura de vértices (vertex cover).

Cobertura de vértices

O problema da cobertura de (arestas por) vértices é sabidamente NP-difícil. Assim, é pouco provável um algoritmo polinomial para resolvê-lo de maneira ótima. No meu trabalho, usei uma heurística simples: percorra os vértices começando por aqueles de maior grau. Se todas as arestas incidentes a esse vértice já foram cobertas por algum vértice da cobertura, ignore-o. Caso contrário adicione-o à cobertura.

Aplicando esse algoritmo nas figuras do começo do post, ficamos com o seguinte resultado:

Ficou mais fácil encontrar as diferenças entre as soluções (clique na imagem para aumentar)

Esse método ainda tem um problema que é quando um par de discos tem sua ordem relativa trocada, mas a parte que foi visualmente modificada está encoberta por um outro disco. Um exemplo disso pode ser visto no canto inferior direito da figura acima. Para consertar esse problema teríamos que fazer um pré-processamento e só considerar vértices cujos discos têm alguma parte do seu borda diferente entre uma solução e outra.

Modelo de programação linear

Como não poderia deixar de ser, apresento aqui o modelo de programação linear inteira da versão mais geral desse problema, na qual os vértices têm um custo associado. Dado um grafo G = (V, E) com uma função de custo c : V \rightarrow \mathbb{R}^+, definimos a variável binária x_v para cada v \in V, que vale 1 se o vértice v está na cobertura e 0 caso contrário. O modelo então fica:

Sujeito a

(1) Para todo (u, v) \in E,

(2) Para todo v \in V,

Meia integralidade

O modelo apresentado tem uma propriedade interessante chamada de meia-integralidade. Isso quer dizer que, dada a relaxação linear desse modelo, existe uma solução ótima onde os valores das variáveis x_v são 0, 1/2 ou 1. A prova é relativamente curta e está descrita na Seção 14.3 de [1].

Assim, dada uma solução fracionária do modelo, podemos usar a técnica de arredondamento, criando uma solução x'_v com x'_v = 1 sempre que x_v \ge 1/2. É possível mostrar que essa solução é uma cobertura de vértices válida e que seu custo não é maior do que 2 vezes o valor da solução ótima, resultando em um algoritmo 2-aproximado.

Referências

[1] Approximation Algorithms – Vijay V. Vazirani (2003)