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)


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)


Mapas de Símbolos Proporcionais

Abril 17, 2011

Sexta-feira passada apresentei uma palestra na série de seminários do LOCo. Falei sobre o problema que estou atacando no mestrado.

Um mapa de símbolos proporcionais emprega símbolos para exibir eventos associados a alguma localidade e intensidade. O símbolo é posicionado no local de ocorrência do evento e a sua área fica sendo proporcional à intensidade desse evento. Focamos em mapas que usam círculos opacos como símbolos. A figura abaixo é um exemplo representando as populações das maiores cidades do Brasil.

Note que há sobreposição entre os discos. Dependendo da ordem em que os empilhamos, pode ser que haja mais ou menos porções dos discos visíveis. Há casos ruins em que a borda do disco fica totalmente encoberta, como na figura abaixo.


Não podemos inferir o centro nem o raio do disco com bordo escondido.

Para tentar fazer com que o desenho tenha bastante borda visível, usamos uma métrica que consiste em maximizar o comprimento visível total dos discos. Com isso em mente, é possível definir um modelo de programação linear inteira, com um modelo que atribui cada disco a um nível.

Além do modelo básico, desenvolvemos algumas desigualdades adicionais para tornar o modelo mais forte, nos baseando em propriedades geométricas, que impedem certas configurações de acontecerem.

Também desenvolvemos duas técnicas de decomposição de instâncias. Um jeito trivial de decompor instâncias é considerar componentes de discos conexas de maneira independente.

Nossa primeira técnica de decomposição vem da observação que um disco que está contido no outro sempre pode ser desenhado por cima em uma solução ótima. Dessa forma, em uma instância como a da figura abaixo, a componente {a, b} pode ser desenhada de maneira independente de {c, d, e, f} e na hora de construir a solução ótima, basta desenhar a solução obtida para {a, b} em cima da solução {c, d, e, f}.

Componente {a,b} está contida em {f} e pode ser resolvida separadamente.

Podemos definir um grafo Gs sobre os discos de S, com conjunto de vértices correspondendo aos discos e há uma aresta (i, j) em Gs se os discos correspondentes se interceptam. Falei sobre esse tipo de grafo em um post anterior.

A outra técnica consiste em remover um ponto de articulação de Gs e replicá-lo nas componentes conexas resultantes, como na figura abaixo.


O disco f representa um ponto de articulação em Gs. As figuras (ii), (iii) e (iv) são as componentes resultantes de sua remoção, com f replicado.

Mostramos que é possível resolver cada componente com o ponto de articulação replicado de maneira independente. Depois, basta juntar os discos mantendo a ordem relativa de cada sub-solução. Isso pode ser feito através de uma ordenação topológica.

Implementamos todas essas ideias e reportamos os resultados. As instâncias testadas foram as mesmas de um trabalho anterior no qual nos baseamos. As técnicas de decomposição foram efetivas na redução do tamanho das instâncias e o resolvedor XPRESS conseguiu resolver nosso modelo para quase todas as instâncias. Porém, ficaram algumas componentes sem serem resolvidas, o que nos motiva a procurar novos modelos e novas desigualdades.

Esse nosso trabalho foi aceito na Computational Geometry and Applications 2011.


Formulações do TSP

Março 4, 2011

O problema do caixeiro viajante (do inglês Traveling Salesman Problem, abreviado como TSP) é o problema mais conhecido e estudado em otimização combinatória. Creio que alguns dos motivos para isso seja sua definição simples, sua dificuldade e sua aplicabilidade prática.

A definição é de fato simples. Imagine que você tem n cidades e entre cada par de cidades existe uma rodovia com uma dada distância. Você deseja partir de uma das cidades, percorrer as outras exatamente uma vez e então retornar para a cidade de origem de forma a minimizar a distância percorrida. Em termos de teoria de grafos, dado um grafo completo de n vértices e custo nas arestas, encontrar um ciclo com exatamente n vértices de menor custo.

O TSP é conhecidamente NP-difícil. Porém, entre os problemas NP-difíceis, há ainda uma diferença na dificuldade entre eles. Algoritmos de aproximação podem dar uma noção dessa diferença. Alguns problemas admitem algoritmos com fatores de aproximação arbitrariamente próximos de 1, enquanto outros não admitem sequer fator de aproximação constante.

O problema da mochila por exemplo, é um dos problemas NP-difíceis mais “tratáveis”. Ele pertence à classe dos PTAS (Polynomial Time Approximation Scheme), o que significa que ele pode ser aproximado por um fator tão próximo de 1 quanto se queira, embora o tempo de execução aumente em função dessa proximidade. Por outro lado, a versão geral do TSP não pode ser aproximada por um fator constante a menos que P = NP. Por isso, o algoritmo do TSP é um dos problemas mais difíceis da computação, embora afirma-se que abelhas saibam resolvê-lo.

Formulações

Há pelo menos três formulações de programação linear inteira para esse problema. Uma delas foi publicada em 1960, por Miller, Tucker e Zemlin e que denominaremos MTZ, é baseada em fluxos em redes. Considere um grafo direcionado completo D(V, A) de n nós e custo nas arestas. Seja uma variável binária que vale 1 se a aresta (i,j) está na solução, e 0 caso contrário. Denotando o custo de uma aresta (i,j) por , a função objetivo é simplesmente

Como queremos encontrar um ciclo que passe por todos os vértices, cada vértice deve ter exatamente uma aresta entrando e outra saindo, o que é gatantido com as seguintes restrições

(1) Para todo

(2) Para todo

Somente com essas restrições, as soluções serão coberturas de vértices por ciclos, ou seja, teremos um ou mais ciclos de forma que cada vértice estará contido em exatamente um ciclo. A ideia dos autores de [1] é introduzir variáveis representando potenciais de um vértice v. Sempre que uma aresta (i,j) é usada, o potencial do vértice j tem que ser maior do que i. Isso criará um impedimento na criação de ciclos, da mesma maneira do modelo para o caminho mais longo, conforme esse post anterior.

Porém, isso também impedirá a formação do ciclo que estamos procurando. O truque é não usar essa restrição sobre um dos vértices, por exemplo o vértice 1. Dessa forma as desigualdades impedirão quaisquer soluções que contenham subciclos que não incluam o vértice 1. As únicas soluções que satisfarão isso, além de (1) e (2), são os ciclos de tamanho n. A desigualdade fica,

(3) Para todo

Onde M é uma constante suficientemente grande. Se , a desigualdade fica redundante. Se , forçamos .

A segunda que vamos apresentar é devido a Dantzig, Fulkerson e Johnson, ou DFJ. O modelo possui as variáveis e as restrições (1) e (2). A diferença fica por conta da eliminação de subciclos.

(4) Para todo

É possível mostrar que se, para um subconjunto não-próprio S de V, existem |S| ou mais arestas, então existe pelo menos um ciclo. O problema é que o número de desigualdades (4) é exponencial e para isso um algoritmo de planos de corte, como o Branch and Cut, deve ser usado. Uma maneira de reduzir o número dessas desigualdades é observar que se uma solução satisfaz (1) e (2) e existe um subciclo de tamanho maior do que |S|/2, então existe um subciclo de tamanho menor do que |S|/2, de forma que podemos considerar as desigualdades (4) apenas para

Uma outra maneira de se eliminar subciclos é através das desigualdade chamadas cut set constraints, dadas por

(5) Para todo

A ideia aqui é que dada uma cobertura por ciclos contendo subciclos, existe uma maneira de separar o conjunto de vértices em 2, de forma que cada subciclo fique exatamente de um dos lados, ou seja, nenhum aresta “cruza” essa partição. Note que o número de desigualdades (5) também é exponencial.

Conclusão

Embora o modelo MTZ tenha sido desenvolvido depois do DFJ, e tenha um número polinomial de desigualdades, na prática o DFJ é bem mais eficiente. A explicação pode ser dada pela combinatória poliédrica, já que é possível mostrar que o modelo DFJ está contido no MTZ. Em um post futuro, pretendo apresentar um esboço dessa prova.

Uma questão que surgiu ao estudar as formulações é a relação entre os modelos com a desigualdade (4) e (5). Serão eles equivalentes? Um está contido no outro? Essa é outra pergunta que pretendo pesquisar.

Referências

[1] Miller, Tucker and Zemlin — Integer Programming Formulation of Traveling Salesman Problems
[2] Dantzig, Fulkerson and Johnson — Solution of a Large-Scale Traveling-Salesman Problem


Lifted Cover Inequalities

Fevereiro 4, 2011

Nesse post falarei sobre as desigualdades de cobertura erguidas (em uma horrorosa tradução livre de lifted cover inequalities). A minha ideia inicial era comentar sobre cortes automáticos que resolvedores como o XPRESS inserem no programa para apertar o modelo. Há dois tipos de desigualdades usadas como cortes: as de cobertura erguidas e aquelas resultantes do procedimento de Chvátal-Gomory.

Como esse segundo método faz uso direto do tableau do simplex, antes eu gostaria de escrever um post revisando o primal simplex e o dual simplex. Por isso, ater-me-ei às coberturas erguidas.

Essas desigualdades surgiram a partir do problema da mochila binária. Vamos relembrar o problema: dado um conjunto de n itens — cada qual com um peso e valor — e uma mochila de capacidade b, queremos escolher um subconjunto de itens cuja soma dos pesos respeite à capacidade da mochila e cuja soma dos valores seja máxima.

Vamos supor que os valores dos itens sejam dados por e os pesos por (para facilitar, suponha que os itens estão ordenados em ordem não crescente de peso, ou seja ). Para cada item i, definimos a variável binária que vale 1 se vamos colocar o item na mochila e 0 caso contrário. O modelo de programação linear inteira é dado por:

Sujeito a:

(1)

(2)

Definições

Definimos qualquer subconjunto de itens que cabem na mochila como conjunto independente. Assim, se S é um conjunto independente, temos que

Se, por outro lado, o subconjunto não cabe na mochila, denominamo-no uma cobertura.

Adicionalmente, uma cobertura é dita minimal se a remoção de qualquer item do subconjunto torna-o independente. Dada uma cobertura C, uma desigualdade de cobertura é definida como:

Ela é claramente válida pois por definição incluir todos os itens de uma cobertura violará a capacidade da mochila.

Cobertura Estendida

Podemos estender essa desigualdade da seguinte forma: Seja j o item mais pesado dessa cobertura. Se colocarmos um item mais pesado que j, ainda assim não podemos ter mais do que |C|-1 itens. Aplicando esse argumento várias vezes, podemos inserir todos os itens j’ < j (que são mais pesados que j por causa da ordenação inicial) na cobertura, para deixar a desigualdade mais apertada. Definimos então a cobertura estendida E(C) como sendo a cobertura C mais os itens mais pesados que o item j. A desigualdade abaixo é então válida:

Lifting das Desigualdades de Cobertura

Uma outra maneira de melhorar uma desigualdade de cobertura é através de lifting, que consiste em aumentar os coeficientes das variáveis do lado esquerdo de uma desigualdade na forma (<=), desde que ela continue viável. Obviamente só conseguimos fazer lifting em desigualdades que não representam facetas, uma vez que a desigualdade que sofreu lifting domina a original.

As desigualdades de cobertura estendida é uma forma de lifting, pois estamos aumentando de 0 para 1 o coeficiente dos itens mais pesados que j. Porém, usando um método diferente (porém mais caro), podemos conseguir aumentar os coeficientes para valores maiores do que 1. Vamos considerar um exemplo de mochila, extraído de [1]:

Uma desigualdade de cobertura válida é

(3)

A cobertura estendida consiste em incluir os itens e , os itens mais pesados que , resultado na desigualdade de cobertura estendida: .

Ao invés de aumentar o valor de para 1 em (3), vamos usar uma variável , levando a

(4)

A pergunta é: qual o maior valor que pode assumir, desde que (4) se mantenha válida? Para isso, temos que encontrar o valor de considerando todas as possibilidades e escolher o menor. Primeiro, se então é ilimitado, o que não ajuda muito. Vamos considerar então que . Com isso, podemos reduzir o problema:

Agora queremos encontrar o maior valor de em (4) que respeite a desigualdade. Para isso, devemos considerar o pior caso, ou seja, em que a soma dos outros termos é máxima,

(5)

A soma máxima pode ser determinada resolvendo-se outro problema da mochila:

Sujeito a

A solução ótima desse problema é 1, o que significa que por (5), . Temos então a desigualdade

(6)

Repetindo esse procedimento para fazer um lifting em (6) aumentando o coeficiente de para , concluiremos que , levando a desigualdade

(7)

Que domina inclusive a desigualdade de cobertura estendida. Porém, isso tem seu preço: para estender uma desigualdade de cobertura só precisamos inserir os itens mais pesados, enquanto nesse último procedimento tivemos que resolver dois novos problemas da mochila. Ou seja, para resolver um problema NP-difícil precisamos resolver sub-problemas que também são NP-difíceis? Veremos…

Programação Dinâmica

Como se sabe, o problema da mochila pode ser resolvido através de um algoritmo clássico de programação dinâmica que possui complexidade pseudo-polinomial de O(nb). Basicamente define-se uma submochila m(i, j), que quer representa o maior valor que se consegue com uma mochila usando os i primeiros itens e com capacidade j. A recorrência é

Porém, é possível indexar o valor da mochila ao invés da capacidade. Nesse sentido definimos m'(i, v) como o menor peso que conseguimos usando os i primeiros itens e com valor v. A recorrência é semelhante:

Com a restrição de que os valores dos itens sejam inteiros. Para encontrar maior valor da mochila, basta percorrer a matriz m’ e determinar o maior v, tal que m'(n, v) <= b.

A complexidade desse segundo algoritmo é O(max_v n), onde max_v é o maior valor que a mochila pode ter (ou seja, um limite superior para o custo da função objetivo). No nosso caso, nossa função objetivo corresponde a um pedaço do lado esquerdo da desigualdade de cobertura. Como tal desigualdade é válida, sabemos que o lado esquerdo nunca ultrapassará |C|-1, logo max_v . Concluímos que os subproblemas da mochila podem ser resolvidos em .

Desigualdades de Cobertura para outros problemas

As desigualdades levantadas para coberturas podem ser usadas para outros problemas. Se temos um PLI na forma para i=1,…,m onde representa a i-ésima linha de A, então podemos particionar P em vários problemas da mochila na forma e encontrar desigualdades fortes para esses subproblemas.

Como temos que

A esperança é que as desigualdades fortes para os subproblemas (em especial as desigualdades de cobertura) sejam úteis para P na prática. Mais detalhes sobre essa técnica, ver em [1].

Bibliografia

[1] Crowder, H. and Johnson, E.L. and Padberg, M. – Solving large-scale zero-one linear programming problems