Sistemas Lineares Aproximados Derivados de Problemas de Fluxo Multiproduto em Métodos de Pontos Interiores

RESUMO

Uma das abordagens utilizadas para resolver o sistema linear que surge a cada iteração nos métodos de pontos interiores primal-dual é reduzi-lo a um sistema linear equivalente simétrico definido positivo, conhecido como sistema de equações normais, e aplicar a fatoração de Cholesky na matriz do sistema. A grande desvantangem desta abordagem é o preenchimento gerado durante a fatoração, o que pode tornar seu uso inviável, por limitação de tempo e memória. Com o intuito de contornar o problema de preenchimento gerado na fatoração de Cholesky, neste trabalho, estamos propondo uma abordagem que resolve de forma direta sistemas lineares aproximados do sistema de equações normais derivados de problemas de fluxo multiproduto e que exerce um certo controle sobre o preenchimento.

Palavras-chave:
Método de pontos interiores primal-dual; Fatoração de Cholesky; Sistema de Equações Normais

ABSTRACT

One of the approaches used to solve the linear system that arises at each iteration in the primal-dual interior points methods is reduce it to an equivalent symmetrical positive defined linear system, known as normal equation system and apply the factorization in the system matrix. The major disavantage of this approach is the filling generated during factorization, which may turn their use impractical for limited time and memory. In order, to deal with the fill-in problem in the Cholesky factorization, in this paper, we proposed an approach that solve directly approximate normal equations linear systems derived from the multiproduct commodity problems and apply some control over the filling.

Keywords:
Primal-dual; interior point methods; Cholesky factorization; Normal equations systems

1 INTRODUÇÃO

Um Problema de Programação Linear (PPL) consiste em encontrar uma solução que satisfaz equações e/ou inequações lineares e que simultaneamente maximize ou minimize uma função linear denominada função objetivo.

Em sua forma padrão este problema é dado por

(1)

em que A é uma matriz de ordem m × n, x ∈ ℝⁿ , b ∈ ℝ^m , c ∈ ℝⁿ e _{^{cT x}} é a função objetivo do problema que deve ser minimizada.

A todo PPL na forma padrão, pode-se associar um problema dual que possui a seguinte forma

(2)

em que z ∈ ℝⁿ e _{^{bT y}} é a função objetivo do problema dual que deve ser maximizada.

Por estarem associados desta forma, denominamos o par de PPLs (1) e (2) de par primal-dual, sendo (1) o problema primal e (2) o problema dual.

Uma solução factível para o par primal-dual é um vetor (x, y, z) que satisfaz as retrições Ax =b, x ≥ 0 _{^{AT y}} + z = c e z ≥ 0. Uma solução ótima (x∗, y∗, z∗) para o par primal-dual é uma solução factível com o menor e maior valor possível para as funções objetivos dos problemas primal e dual, respectivamente.

O Teorema 1, enunciado a seguir, apresenta as condições de Karush-Kuhn-Tucker (KKT) para problemas de progração linear, que oferece uma condição necessária e suficiente para que um vetor (x, y, z) seja uma solução ótima para o par primal-dual.

Theorem 1. (Condicões de otimalidade) ¹⁹. (x, y, z) são solucões ótimas para o par primaldual se, e somente se, as seguintes condicões são satisfeitas

(3)

(4)

(5)

(6)

1.1 Método de Pontos Interiores Primal-Dual

Os métodos de pontos interiores (MPI) para problemas de programação linear, como a maioria dos algoritmos iterativos em otimização, possuem dois pontos-chave: um procedimento quedetermina uma direção do passo ou direção de busca em cada iteração e uma medida conveniente deste passo.

Nos métodos de pontos interiores do tipo primal-dual a direção de busca é calculada por meiode um sistema linear obtido pela aplicação de variações do método de Newton ¹⁷ às três igualdades (3), (4) e (5) ¹⁹.

Logo, a direção de busca (dx, dy, dz) em cada iteração pode ser calculada por meio do sistema linear

(7)

em que X e Z são matrizes diagonais, cujos elementos da diagonal são dados ordenadamente pelas entradas dos vetores x e z, repsectivamente, os resíduos são da forma _^rp = b − Ax,_^rd = c − _{^{At y}} −z e _^rc = −Xze + σ μe, com μ = _{^{xt z/n}} (gap de dualidade) e σ é um parâmetro de centralização entre 0 e 1.

O tamanho do passo é dado de forma que as variáveis x e z sejam estritamente positivas.

Em aplicações reais o sistema linear (7) quase sempre possui dimensões elevadas e um alto grau de esparsidade e sua resolução é a etapa do método que requer mais esforço computacional.

Por meio de operações algébricas, o sistema pode ser reduzido a um sistema linear equivalente cuja matriz é simétrica definida positiva dado por

(8)

em que D = XZ ⁻¹ e

(9)

Esse sistema linear é conhecido como sistema de equações normais. Uma das abordagens utilizadas para resolvê-lo é a aplicação da fatoração de Cholesky na matriz dos coeficientes, já que ela é simétrica definida positiva.

A fatoração de Cholesky de uma matriz simétrica definida positiva A de ordem m consiste em calcular uma matriz L triangular inferior, denominada fator de Cholesky, tal que A = _^LLT . Os elementos da diagonal de L são denominados pivôs.

A solução de um sistema linear com a matriz A, pode então, ser calculada por meio da resolução de dois sistemas triangulares com as matrizes L e _^LT . Esta abordagem é uma forma bastante estável para calcular soluções de sistemas lineares com a matriz dos coeficientes simétrica definida positiva, por isso, ela é utilizada nos MPI primal-dual sempre que a decomposição não é extremamente cara.

A grande desvantagem da fatoração de Choleky é o preenchimento gerado durante a decomposição, que consiste no surgimento de elementos não nulos no fator de Cholesky nas posições em que na matriz original os elementos são nulos. A Observação 1 mostra como o preenchimento ocorre.

Observação 1.Considere a matriz simétrica definida positiva de ordem 5×5, com o padrão de esparsidade dado à esquerda, em que × são elementos não nulos. Apenas a parte triangular inferior da matriz está representada. A parte superior pode ser obtida pela transposição da inferior.

O cálculo da segunda coluna de L gera um preenchimento na posição (4,2), como ilustrado na matriz à direta, indicado com o sinal + nesta posição. Uma vez calculada a segunda coluna, o cálculo da terceira coluna produzirá um elemento não-nulo na posição (4,3), também indicado pelo símbolo +.

O preenchimento gerado na fatoração de Cholesky leva à necessidade de mais espaço de armazenamento e maior tempo computacional para a resolução dos sistemas lineares. Logo, seu uso pode se tornar inviável por limitação de tempo e memória, fazendo com que a utilização de métodos iterativos se torne mais apropriada. Neste contexto, a escolha mais natural de um método para resolver o sistema de equações normais é o método dos gradientes conjugados precondicionado (GCP).

Em geral, a matriz do sistema de equações normais sofre mudanças significativas entre as iterações dos MPI primal-dual e se torna extremamente mal condicionada nas iterações finais, por isso, junto ao método GCP nos MPI é mais adequado a utilização de precondicionadores híbridos compostos por mais de um tipo de precondicionador, em que cada um é mais apropriado para uma determinada etapa da resolução do problema ^9),(18.

Um precondicionador que se mostra eficiente para as iterações iniciais dos MPI é construído pela fatoração controlada de Cholesky (FCC)⁶, que é um tipo de fatoração incompleta de Cholesky.

Neste tipo de fatoração, ao fatorar uma matriz simétrica definida positiva, é obtida uma matriz triangular inferior é mais esparsa que o fator de Cholesky. Na seção 2, a fatoração controlada de Cholesky será descrita mais detatalhadamente., denominada fator incompleto de Cholesky, cujos elementos são calculados como na fatoração de Cholesky, porém devido a algum tipo de controle de preenchimento,

A consequência de resolver os sistemas lineares que surgem dos métodos de pontos interiores primal-dual de forma iterativa, é que a direção de busca em cada iteração é calculada aproximadamente. A direção calculada satisfará um sistema da forma

(10)

em que e , são componentes de um erro residual.

Se o sistema de equações normais é resolvido por um método iterativo então a direção dy calculada satisfaz

em que h é como em (9) e (10) será dado por: é um erro residual. Desta forma, o erro residual de

(11)

(12)

(13)

Os métodos de pontos interiores que resolvem o sistema de Newton de forma aproximada são conhecidos como métodos de pontos interiores inexatos. Em geral, nestes métodos, uma direção de busca calculada é aceitável se o erro residual satisfaz a condição:

em que _^yk é uma sequência uniformemente menor que um.

As análises dos métodos de pontos interiores inexatos mostram que se o erro residual é controlado, ainda é possível calcular uma boa direção de busca em cada iteração ^1)-(3),(15.

2 FATORAÇÃO CONTROLADA DE CHOLESKY

Em uma fatoração incompleta de Cholesky de uma matriz simétrica definida postiva A deordem m, uma vez calculada uma coluna de R, em que R é uma matriz denominada matriz resto., tal que A = como na fatoração de Cholesky, o controle sobre o preenchimento é feito descartando por algum critério um certo número de elementos não-nulos da coluna. Após este descarte, passa-se para o cálculo da próxima coluna, realizando o mesmo procedimento. Desta forma, é obtida uma matriz triangular inferior −

Em geral, as fatorações incompletas de Cholesky são utilizadas como precondicionadores para os métodos de gradientes conjugados precondicionados. Em ⁸ é mostrado que o número de iterações do GCP está quase que diretamente relacionado com a norma da matriz R. Isto motivou Campos ⁶ a construir a FCC baseada na minimização da norma de Frobenius da matriz E = L, em que L é o fator de Cholesky da matriz A. De fato, reescrevendo R da forma -

pode-se observar que, ||R|| → 0 quando ||E|| → 0.

O problema de minimização da norma de Frobenius da matriz E pode ser escrito como

(14)

com _^cj = ∑_{^mi =1}|_{^{li j}} − , i,j=1,...,m são os elementos de L e |2, em que li j e , respectivamente. Os elementos cj, por sua vez, podem ser reescritos como soma de dois somatórios

(15)

em que _^nj é o número de elementos não nulos da coluna j de A e η é um parâmetro que determina a quantidade de preenchimento adicional permitida em relação a _^nj , na coluna j de .

A primeira parcela de (15) contém todos os elementos não nulos da coluna j de L, cujo correspondente em _{^{li j}} , uma redução da norma de Frobenius de E pode ser obtida por meio da seguinte heurística: e a segunda contém apenas os elementos restantes do fator completo ≈ é nulo. Se consideramos que

• Aumentar o valor de η permitindo mais preenchimento. Isto faz com que _^cj , j=1, ..., m e consequentemente ||E||_F decresça, pois o primeiro somatório em (15) conterá mais elementos e o segundo menos elementos.

• Para um η fixo, escolher os _^ηj + η maiores elementos da coluna j de L, para j=1, ..., m, em valor absoluto. Assim, para uma determinada quantidade de armazenamento, um fator incompleto ótimo é obtido, já que os maiores elementos estarão no primeiro somatório e os menores no segundo.

A FCC combina os dois itens acima. Uma vez calculada a coluna j de , ela propõe utilizar um parâmetro η tomando o número máximo de elementos não-nulos desta coluna como sendo kj = η + ηj. Os kj elementos maiores em valor absoluto calculados, serão mantidos na coluna j de e os demais descartados, isto é, serão nulos.

Esta heurística permite prever o armazenamento necessário na fatoração e exercer um controle de preenchimento de acordo com a memória disponível.

Como para um dado problema, a ordem m e o número de elementos não nulos de A são fixos, o número de elementos não nulos de η, que por sua vez, pode assumir qualquer valor de -m a m. Se η =-m, a matriz η =0, ela requer o mesmo espaço de armazenamento da matriz A e se η =m, a matriz será diagonal, se depende apenas do parâmetro é o fator de Cholesky completo.

2.1 Tratamento de pivôs inadequados na fatoração incompleta de Cholesky

A existência de uma fatoração incompleta foi estabelecida apenas para certas classes de matrizes, como a classe das M-matrizes e a classe das H-matrizes, na qual as matrizes diagonalmente dominantes estão incluídas ⁴. Para matrizes simétricas definidas mais gerais, fatorações incompletas de Cholesky podem falhar devido a ocorrência de pivôs nulos ou negativos calculados durante a fatoração. Infelizmente, muitas aplicações importantes lidam com matrizes que não são M-matrizes ou H-matrizes.

Pivôs inadequados, são reconhecidos como um problema desde o início do precondicionamento por fatorações incompletas, e inúmeras estratégias foram propostas para remediá-lo.

Uma boa estratégia, foi sugerida por Manteuffel ¹⁴. Nela, se uma fatoração incompleta de Cholesky de A falhar devido a pivôs inadequados, é realizado um incremento global na diagonal de A, antes de se tentar uma nova fatoração. Assim, a fatoração incompleta é aplicada na matriz A + α diag(A), em que α > 0 e diag(A) denota a parte diagonal de A. Se a fatoração incompleta falhar novamente, é realizado um incremento em α e a fatoração é reiniciada. O processo é repetido até que uma fatoração incompleta seja calculada com sucesso. =

Esta estratégia é baseada no fato de que existe um valor α* para o qual a fatoração incompleta de α* como sendo o menor α para o qual A + α diag(A) é diagonalmente dominante. Contudo, bons resultados são geralmente obtidos para valores de α bem menores que α* e maiores que o menor α para o qual = existe. Pode-se, por exemplo, tomar admite uma fatoração incompleta. A desvantagem desta abordagem é a dificuldade em se determinar um valor de αque não seja muito maior que o necessário para evitar novos reinícios.

O código da fatoração controlada de Cholesky de Campos utiliza a estratégia de Manteuffel. Bocanegra, em seu precondicionador híbrido ⁵, utilizou este código, com α determinado pela fórmula

(16)

tomando λ = 5 · 10⁻⁴ e i =1,2,...,15, em que i é o número de reinícios da fatoração.

3 NOVA ABORDAGEM NA SOLUÇÃO DIRETA DO SISTEMA DE EQUAÇÕES NORMAIS

Se por um lado a fatoração de Cholesky do sistema de equações normais oriundos dos MPI primal-dual gera muitos preenchimentos, por outro, a resolução deste sistema pelo método dos gradientes conjugados, exige a utilização de precondicionadores sofisticados. Neste trabalho, com o intuito de contornar o problema de preenchimento gerado na fatoração de Cholesky e com a motivação de que uma direção aproximada com boas propriedades pode ser calculada no método de pontos interiores, é proposta uma abordagem que de certa forma combina características dos métodos acima. A abordagem resolve de forma direta sistemas lineares aproximados do sistema de equações normais, mas exercendo um certo controle sobre o preenchimento. Isto é feito substituindo a fatoração de Cholesky da matriz do sistema de equações normais, pela fatoração controlada de Cholesky.

Desta forma, ao invés de resolver em cada iteração o sistema linear _^LLTdy = h, onde L é o fator de Cholesky da matriz _^ADAT , é resolvido ainda diretamente, um sistema linear da forma

(17)

em que _^ADAT , calculado por meio da fatoração controlada de Cholesky. Assim, a cada iteração do método é encontrada uma solução satisfaz o sistema que é uma aproximação da solução original, isto é, é o fator incompleto da matriz

(18)

em que .

Se (11) satisfaça a condição de convergência, o método proposto realiza uma única vez substituição regressiva e progressiva para resolver os sistemas lineares com as matrizes triangulares é controlado, o método proposto é de certa forma um método de pontos interiores inexato, no sentido de que uma direção aproximada é calculada. Entretanto, difere dos métodos inexatos usuais que resolvem o sistema de equações normais. Os métodos inexatos usuais utilizam um método iterativo para resolver o sistema de equações normais “exato”, já a abordagem proposta resolve de forma direta um sistema aproximado do sistema de equações normais. Assim, enquanto a cada iteração dos métodos inexatos usuais, muitas iterações do método iterativo são feitas, até que o resíduo .

O sucesso de um MPI inexato depende de tornar o resíduo que surge do lado direito da equação de Newton menor do que uma tolerância aceitável para atingir a convergência.

Observamos que, tomando a norma dois de , obtém-se

Então, se ||L|| F, ||||F→0. Logo, na abordagem proposta é desejável minimizar a norma de Frobenius de E = L - ||F e ||||2→0, quando ||L - ,o que justifica a escolha da FCC para a construção do fator incompleto.||2 forem limitadas, temos que ||

Além disso, o parâmetro de preenchimento η da FCC pode ser controlado convenientemente, de forma a iniciar com um valor que permita na primeira iteração, uma matriz o mais esparsa possível. À medida que o MPI avança nas iterações, o valor de η aumenta progressivamente, permitindo um aumento gradual da densidade de ao longo das iterações, até atingir a fatoração de Cholesky completa. Ao mesmo tempo espera-se que a convergência do método seja alcançada. Desta forma, a resolução do sistema se torna mais rápida nas iterações iniciais, já que haverá uma economia de tempo tanto no cálculo dos elementos do fator incompleto, quanto na resolução dos dois sistemas triangulares.

4 EXPERIMENTOS COMPUTACIONAIS

A abordagem proposta foi incorporada ao código PCx ⁷ que é uma implementação do método primal-dual preditor-corretor. Todas as rotinas deste código são implementadas em Linguagem C, exceto a parte responsável pela solução dos sistemas lineares que é feita aplicando a fatoração de Cholesky no sistema de equações normais, utilizando o código da fatoração de Cholesky esparsa projetada por Ng e Peyton ¹⁶, em FORTRAN77. Na implementação do código com a nova abordagem, que será denominada código PCx modificado (PCx_mod ), foi admitida a existência de duas fases. Na primeira fase são calculados fatores de Cholesky incompletos, utilizando um código da FCC implementado em FORTRAN77 que foi fornecido pelo seu autor ⁶, mas com algumas adaptações. Na segunda fase, são calculadas fatorações de Cholesky completas utilizando o código de Ng e Pyton. Isto porque, o código da FCC não é eficiente para o cálculo do fator completo, pois necessita de operações adicionais e não explora o fato de que a estrutura esparsa da fatoração final é conhecida.

A mudança para a segunda fase é feita quando a quantidade de elementos não nulos do fator incompleto calculado é 95% da quantidade de elementos não nulos do fator de Cholesky completo³ ou quando não está sendo obtido progresso com a fatoração controlada de Cholesky. O progresso de uma iteração para outra é medido através do quociente

são os gaps de dualidade da iteração atual e anterior, respectivamente. Na implementação foi considerado que se ρ ≥ 0,99, nenhum progresso está sendo obtido, assim a mudança para a segunda fase é feita na iteração imediatamente seguinte.

Ainda na primeira fase, uma vez determinado um valor de η inicial, a partir da segunda iteração, seu valor é atualizado ou não, em função de ρ, da seguinte forma:

em a e b são constantes, com a < b.

Desta forma, se um bom progresso é obtido de uma iteração para outra, isto é, se há uma redução significativa do gap de dualidade, o valor do parâmetro η não muda de uma iteração para outra. E se pouco progresso é obtido na redução do gap de dualidade, é realizado um incremento maior em η.

Na seção 2.1, foi mencionado o problema do surgimento de pivôs inadequados nas fatorações incompletas de Cholesky. No método proposto, para contornar este problema, foi utilizada a mesma técnica que o código PCx utiliza quando pivôs inadequados surgem na fatoração de Cholesky. Um pivô ^-8 é substituído por 10¹²⁸ e a fatoração simplesmente continua a partir daí, evitando reiniciar a fatoração. Isso faz com que os elementos da coluna i de ¹⁰. sejam essencialmente nulos, replicando técnicas utilizadas para fatoração de matrizes semidefinidas positivas menor que 10

Os testes apresentados nesta seção foram realizados com os problemas de fluxo multiproduto do tipo PDS da biblioteca Kennington ¹¹. A Tabela 1 exibe em suas colunas, nesta ordem, o nome dos problemas, o número de linhas, o número de colunas, a quantidade de elementos não nulos da matriz _^AAT (|_^AAT| ) e a densidade do fator de Cholesky desta última matriz.

Tabela 1
Problemas testes.

Problema	Linhas	Colunas	|AAT |	Densidade Chol
pds-10	16558	49932	77866	0,014
pds-20	33874	108175	170973	0,130
pds-30	49944	158489	251580	0,013
pds-40	64265	214385	342483	0,013
pds-50	83060	275814	432371	0,013
pds-60	96503	332862	524563	0,012
pds-70	114944	390005	607375	0,013
pds-80	126120	430800	677160	0,012
pds-90	142823	475448	741153	0,011
pds-100	156243	514577	788421	0,010

Nos testes foram utilizados os valores 0, para o valor inicial de η, a = 10 e b = 25, sendo estes os melhores valores obtidos para os problemas PDS, dentre alguns valores testados.

A Tabela 2, exibe uma comparação, na segunda e na terceira coluna, entre o número de iterações dos códigos PCx e PCx_mod , respectivamente. Na quarta coluna temos o número de iterações da primeira fase do código, ou seja, o número de iterações em que se utilizou um fator incompleto para calcular uma direção de busca. Nas duas últimas colunas da tabela, pode-se comparar o tempo de processamento, em segundos, para a resolução dos problemas pelos dois códigos. A última linha da tabela exibe a soma total de cada coluna.

Tabela 2
Comparação do número de iterações e do tempo de processamento.

Problema	Iteração			Tempo
	PCx	PCxmod	FCC	PCx	PCxmod
pds-10	40	55	20	30,61	17,79
pds-20	52	66	29	177,01	140,09
pds-30	62	66	26	563,14	338,26
pds-40	69	67	25	1826,69	978,52
pds-50	68	72	30	2760,82	1700,20
pds-60	69	77	34	4495,51	2883,73
pds-70	72	82	39	7425,77	4562,90
pds-80	74	75	33	10631,58	6083,00
pds-90	75	83	45	13348,96	6845,25
pds-100	80	84	11	15724,24	14222,93
Total	661	727	292	56984,33	37772,67

Com exceção do problema pds-40, houve um aumento no número de iterações dos demaisproblemas, somando ao todo uma aumento de 66 iterações, como pode ser visto na Tabela 2. Entretanto, em contrapartida, houve uma redução significativa no tempo de processamento de todos os problemas, totalizando uma redução de 19211,66 segundos.

A FCC foi utilizada em 40% do total de iterações e a redução total do tempo de processamento foi de 33,7%. A grande redução obtida no tempo se deveu à resolução dos sistemas lineares mais esparsos nas iterações iniciais.

A Tabela 3 permite analisar, para cada problema, o quanto o fator incompleto de Cholesky calculado na primeira iteração, é menos denso que o fator de Cholesky. Isto pode ser observado, comparando o número de elementos não nulos do fator de Cholesky L e do fator incompleto da primeira iteração, denotado por _^AAT . Os valores do quociente, exibidos na segunda e terceira coluna da tabela, respectivamente. Note que, neste caso, como o η inicial é igual a zero, o número de elementos não nulos da matriz é igual ao número de elementos não nulos da matriz

exibidos na quarta coluna, nos fornecem uma melhor visão da esparsidade de L, pois, se κ é próximo de zero, então L e quando κ é próximo de um, é muito mais esparso que é quase o fator de Cholesky. em relação a

Tabela 3
Comparação entre Cholesky completa e o incompleta.

Problema	|L|	|0 |	κ
pds-10	1687660	79866	0,047
pds-20	7089645	170973	0,024
pds-30	14643812	251580	0,017
pds-40	28195225	342483	0,012
pds-50	41767191	432371	0,010
pds-60	58118583	5224563	0,009
pds-70	78724946	607375	0,007
pds-80	94270275	677160	0,008
pds-90	110971961	741153	0,007
pds-100	120240013	788421	0,006

Pode-se observar na Tabela 3 que os primeiros fatores incompletos das matrizes dos sistemas de equações normais dos problemas do tipo PDS são extremamente esparsos em comparação ao fator de Cholesky completo, o que justifica a redução no tempo de processamento total apesar do aumento do número de iterações.

4.1 Erro residual

Na prática, os métodos de pontos interiores inexatos usuais, que resolvem o sistema de equações normais, utilizam um método iterativo para resolver este sistema. O critério de parada do método iterativo é baseado no controle do erro residual, assim, as iterações são realizadas até que o erro residual atinja uma condição pré-estabelecida para que a convergência ocorra.

Na abordagem proposta, em uma iteração de métodos de pontos interiores, não é possível a priori, controlar o erro residual, podendo ocorrer de ele ser maior do que a condição pré-estabelecida para a convergência. Entretanto, dado que erro residual no lado direito do sistema de equações normais satisfaz a condição

caso seja necessário, pode-se aumentar a densidade da matriz L, até que ||₂ satisfaça a condição de convergência. Na pior das hipóteses, pode ser tomado L. de forma que ela se aproxime mais da matriz = ||

Assim, calculada uma direção aproximada em uma iteração de método de pontos interiores, é verificado se ela satisfaz a condição de convergência. Se isto não ocorre, descarta-se a direção calculada, calcula-se uma nova matriz permitindo mais preenchimentos e uma nova direção é calculada utilizando esta matriz.

Em ¹³, Korzak propõe um método de pontos interiores primal-dual infactível e inexato baseado no método proposto por Kojima, Mizuno e Yoshise em ¹². No método de Korzak a convergência é garantida se as componentes do erro residual no sistema (10) satisfazem as condições:

(19)

(20)

(21)

Com τ₁ ∈ (0, 1>, τ₂ ∈ (0, 1> e τ₃ ∈ (0, 1>. Na prática, são utilizados τ₁ = 0,95, τ₂ = 0,95 e τ₃ = 0,049.

Baseado no trabalho de Korzak, foram realizados alguns testes computacionais para verificar o comportamento do erro residual na abordagem proposta. Neste caso, as condições (20) e (21) são atendidas, já que Tabela 4 exibe estes valores para as quinze primeiras iterações de cinco dos problemas testes. = 0. Assim, como , é verificado apenas o valor da norma no processamento de alguns problemas, nas iterações em que a fatoração incompleta de Cholesky foi realizada. A = 0 e

Tabela 4
Valor da norma nas quinze primeiras iterações.

Iteração	pds-10	pds-30	pds-50	pds-70	pds-90
1	4∙10-2	4∙10-2	3∙10-2	3∙10-2	2∙10-2
2	1∙10-4	6∙10-3	1∙10-3	4∙10-3	5∙10-3
3	3∙10-4	3∙10-4	5∙10-4	1∙10-4	2∙10-4
4	2∙10-5	6∙10-5	2∙10-5	3∙10-5	2∙10-5
5	3∙10-4	1∙10-4	1∙10-4	1∙10-5	2∙10-5
6	3∙10-4	2∙10-4	1∙10-5	1∙10-5	5∙10-5
7	9∙10-4	3∙10-4	1∙10-4	2∙10-5	1∙10-4
8	2∙10-4	4∙10-4	2∙10-4	1∙10-4	1∙10-4
9	1∙10-4	4∙10-4	5∙10-4	1∙10-3	5∙10-4
10	1∙10-4	4∙10-4	1∙10-3	1∙10-4	6∙10-4
11	1∙10-4	1∙10-4	9∙10-4	7∙10-4	7∙10-4
12	1∙10-3	1∙10-4	1∙10-3	5∙10-4	6∙10-4
13	-	1∙10-4	1∙10-3	5∙10-4	4∙10-3
14	-	1∙10-4	1∙10-4	9∙10-4	1∙10-3
15	-	1∙10-3	8∙10-4	1∙10-3	1∙10-3

Na Tabela 4 pode ser verificado que o valor de (19) é atendida nestas iterações. é menor do que 0,05 em todas as iterações. Assim a condição de convergência

A desvantagem da abordagem proposta é que uma direção aproximada calculada pode não satisfazer a condição de convergência, entretanto nos testes computacionais verificamos que nas quinze primeiras iterações dos problemas testados, nenhuma direção inadequada foi calculada. Mas caso isso ocorra, a direção inadequada é descartada e aumentando a densidade da matriz calculada pela fatoração incompleta, é possível calcular uma nova direção mais adequada.

5 CONCLUSÕES

Neste trabalho foi apresentado uma nova abordagem para a solução do sistema de equações normais oriundos dos problemas de fluxo multiproduto em métodos de pontos interiores.

O método proposto é de certa forma um método de pontos interiores inexato, no sentido de que uma direção aproximada é calculada em cada iteração. Mas diferente dos métodos inexatos usuais que resolvem o sistema de equações normais “exato” por meio de um método iterativo, a abordagem proposta resolve de forma direta um sistema aproximado do sistema de equações normais exercendo um controle sobre o preenchimento.

Os testes computacionais, mostraram que a abordagem proposta é bastante eficiente para a resolução dos problemas de fluxo multiproduto do tipo PDS, sendo competitivo com o código PCx.

A desvantagem do método é que não é possível prever a priori, se uma direção calculada satisfaz as condições de convergência, entretanto esta direção pode ser descartada e uma uma nova direção, mais adequada, pode ser eventualmente calculada.

Como trabalho futuro iremos extender este método para outros tipos de problemas de programação linear.

6 AGRADECIMENTOS

Os autores agradecem à FAPESP e ao CNPq pelo suporte financeiro deste trabalho.

REFERÊNCIAS

• ¹
  G. Al-Jeiroudi & J. Gondzio. About the “Convergence Analysis of the Inexact Infeasible Interior-Point Method for Linear Optimization”, Journal of Optimization Theory and Applications, 141 (2009), 231-247.
• ²
  V. Baryamureeba & T. Steihaug. “On the convergence of an inexact primal-dual interior point method for linear programming”, in Lecture Notes in Computer Science, Springer, Berlin/Heidelberg (2006).
• ³
  S. Bellavia. “An inexact interior point method”, Journal of Optimization Theory and Applications,96 (1988), 199-121.
• ⁴
  M. Benzi. “Preconditioning techniques for large linear systems: A survey”. Journal of Computational Physics, 182(2) (2002), 418-477.
• ⁵
  S. Bocanegra, F.F. Campos & A.R.L. Oliveira. “Using a hybrid preconditioner for solving large-scale linear systems arising from interior point methods”, Computational Optimization and Applications, 36 (2007), 149-164. Special issue on “Linear Algebra Issues arising in Interior Point Methods”.
• ⁶
  F.F. Campos. “Analysis of Conjugate Gradients - type methods for solving linear equations”, PhD thesis, Oxford University Computing Laboratory, Oxford (1995).
• ⁷
  J. Czyzyk, S. Mehrotra, S.J. Wright & M. Wagner. “PCx: An interior-point code for linear programming”, Optimization Methods and Software, 11(1) (1999), 397-430.
• ⁸
  I.S. Duff & G.A. Meurant. “The effect of ordering on preconditioned conjugate gradients”, BIT, 29 (1989), 635-657.
• ⁹
  C.T.L.S. Ghidini, A.R.L. Oliveira & D.C. Sorensen. “Computing a hybrid preconditioner approach to solve the linear systems arising from interior point methods for linear programming using the gradient conjugate method”, Annals of Management Science, (2013).
• ¹⁰
  N. Highan. “Acuraccy and Stability of Numerical Algorithms”, SIAM Publications, Philadelphia, (1996).
• ¹¹
  J.L. Kennington, S. Niemi & S.J. Wichmann. “An empirical evaluation of the korbx algorithms for military airlift applications”, Operations Research, 28, (1990), 240-248.
• ¹²
  M. Kojima, S. Mizuno & A. Yoshise. “A primal-dual interior-point method for linear programming. Progress in mathematical programming”, interior-point and related methods, Spring-Verlag, New York, 61 (1989), 29-47.
• ¹³
  J. Korzak. Convergence analysis of inexact infeasible-interior-point algorithms for solving linear programming problems”, Siam J. Optim, 11 (2000), 133-148.
• ¹⁴
  T.A. Manteuffel. “An incomplete factorization technique for positive definite linear systems”, Math. Comp., 34 (1980), 473-497.
• ¹⁵
  R.D.S. Monteiro & J.W. O’Neal. “Convergence analysis of long-step primal-dual infeasible interior point LP algorithm based on iterative linear solvers”, Georgia Institute of Technology (2003).
• ¹⁶
  E. Ng & B.P. Peyton. “Block sparse Cholesky algorithms on advanced uniprocessors computers”, SIAM Journal on Scientific Stat. Computing, 14 (1993), 1034-1056.
• ¹⁷
  R.J. Vanderbei. “Linear Programming - Foundations and Extensions”, Kluwer Academics Publishers, Boston, USA (1996).
• ¹⁸
  M.I. Velazco, A.R.L. Oliveira & F.F. Campos. “Heuristics for implementation of a hybrid preconditioner for interior-point methods”, Pesquisa Operacional, 34 (2011), 2553-2561.
• ¹⁹
  S.J. Wright. “Primal-Dual Interior-Point Methods”, SIAM Publications, SIAM, Philadelphia, PA, USA, (1996).

• 3
  O número de elementos não nulos do fator de Cholesky completo é calculado pelo código de Ng e Peyton na fase incial, antes de iniciar as iterações.

Datas de Publicação

• Publicação nesta coleção
  Jan-Apr 2017

Histórico

• Recebido
  20 Nov 2016
• Aceito
  18 Mar 2017