156
AVALIAÇÃO DE ALGORITMOS HEURÍSTICOS DE OTIMIZAÇÃO EM PROBLEMAS DE ESTIMAÇÃO DE PARÂMETROS Marcio Schwaab DISSERTAÇÃO SUBMETIDA AO CORPO DOCENTE DA COORDENAÇÃO DOS PROGRAMAS DE PÓS-GRADUAÇÃO DE ENGENHARIA DA UNIVERSIDADE FEDERAL DO RIO DE JANEIRO COMO PARTE DOS REQUISITOS NECESSÁRIOS PARA A OBTENÇÃO DO GRAU DE MESTRE EM CIÊNCIAS EM ENGENHARIA QUÍMICA. Aprovada por: ________________________________________________ Prof. José Carlos Costa da Silva Pinto, D.Sc. ________________________________________________ Prof. José Luiz Fontes Monteiro, D.Sc. ________________________________________________ Prof. Evaristo Chalbaud Biscaia Jr., D.Sc. ________________________________________________ Prof. Marcelo Castier, Ph.D. RIO DE JANEIRO, RJ - BRASIL SETEMBRO DE 2005

Marcio Schwa Ab

Embed Size (px)

Citation preview

Page 1: Marcio Schwa Ab

AVALIAÇÃO DE ALGORITMOS HEURÍSTICOS DE OTIMIZAÇÃO EM

PROBLEMAS DE ESTIMAÇÃO DE PARÂMETROS

Marcio Schwaab

DISSERTAÇÃO SUBMETIDA AO CORPO DOCENTE DA COORDENAÇÃO DOS

PROGRAMAS DE PÓS-GRADUAÇÃO DE ENGENHARIA DA UNIVERSIDADE

FEDERAL DO RIO DE JANEIRO COMO PARTE DOS REQUISITOS

NECESSÁRIOS PARA A OBTENÇÃO DO GRAU DE MESTRE EM CIÊNCIAS EM

ENGENHARIA QUÍMICA.

Aprovada por:

________________________________________________ Prof. José Carlos Costa da Silva Pinto, D.Sc.

________________________________________________ Prof. José Luiz Fontes Monteiro, D.Sc.

________________________________________________ Prof. Evaristo Chalbaud Biscaia Jr., D.Sc.

________________________________________________ Prof. Marcelo Castier, Ph.D.

RIO DE JANEIRO, RJ - BRASIL

SETEMBRO DE 2005

Page 2: Marcio Schwa Ab

ii

SCHWAAB, MARCIO

Avaliação de Algoritmos Heurísticos de

Otimização em Problemas de Estimação de

Parâmetros [Rio de Janeiro] 2005.

VII, 149 p. 29,7 cm (COPPE/UFRJ, M.Sc.,

Engenharia Química, 2005).

Dissertação - Universidade Federal do Rio

de Janeiro, COPPE.

1. Estimação de Parâmetros.

2. Algoritmos Heurísticos de Otimização.

3. Enxame de Partículas.

4. Região de Confiança dos Parâmetros.

I. COPPE/UFRJ II. Título ( série )

Page 3: Marcio Schwa Ab

iii

A meus pais, Inácio e Clarinda, a meus

irmãos, Edison e Fernando, e à Elisa pela

amizade e apoio. À Deus por sempre

colocar as pessoas certas na hora certa em

minha vida. Muito que tenho hoje devo a

estas pessoas.

Page 4: Marcio Schwa Ab

iv

AGRADECIMENTOS

Aos meus orientadores José Carlos Pinto e José Luiz Monteiro, pela amizade,

pelos ensinamentos e pela oportunidade concedida.

À Geraldo Crossetti, à Luciana Silva e à Marcos Lopes Dias pelos valiosos

dados experimentais da reação de polimerização.

À Ariane, à Neuman Resende e à Vera Salim pelos valiosos dados experimentais

da reação de produção de gás de síntese.

A todos os colegas, funcionários e professores do PEQ/COPPE/UFRJ, em

particular ao pessoal do LMSCP e do NUCAT.

Aos professores José Eduardo Olivo e Nádia Regina Camargo Fernandes

Machado da Universidade Estadual de Maringá, pelo apoio e incentivo para a minha

vinda ao Rio de Janeiro.

Aos meus velhos amigos da época do Handebol, em particular ao grande José

Carlos Mendes, mais conhecido como Spock, pela grande oportunidade que ele me

concedeu.

À Neuza Patrício de Santanna e à João Carlos de Azevedo pelas inúmeras cópias

de artigos na biblioteca e pelo bom humor durante o atendimento.

Aos meus pais, Inácio e Clarinda, que muito me ajudaram e se esforçaram para

que a realização da minha graduação e a minha vinda ao Rio de Janeiro fossem

possíveis.

Ao meu irmão Édison, pela amizade e pela grande ajuda que ele deu a nossos

pais na mudança para Santa Catarina.

Ao meu irmão Fernando, pela amizade e pelas Ervas de Tereré.

À Elisa por seu carinho e amizade.

A todos do futebol do PEQ, sem eles nada disto seria possível.

Ao grande Estado do Paraná e ao grande Atlético Paranaense, o melhor time de

futebol do mundo.

Ao CNPq pelo apoio financeiro.

A todos que contribuíram de alguma forma para a realização deste trabalho.

A todos, o meu muito obrigado.

Page 5: Marcio Schwa Ab

v

Resumo da Dissertação apresentada à COPPE/UFRJ como parte dos requisitos

necessários para a obtenção do grau de Mestre em Ciências (M.Sc.)

AVALIAÇÃO DE ALGORITMOS HEURÍSTICOS DE OTIMIZAÇÃO EM

PROBLEMAS DE ESTIMAÇÃO DE PARÂMETROS

Marcio Schwaab

Setembro/2005

Orientadores: José Carlos Costa da Silva Pinto

José Luiz Fontes Monteiro.

Programa: Engenharia Química

Este trabalho avalia a utilização de algoritmos heurísticos de otimização na

estimação de parâmetros de modelo não lineares. Os métodos heurísticos usados foram

o Monte Carlo, o Algoritmo Genético, o Recozimento Simulado e o Enxame de

Partículas. Estes métodos apresentam diversas vantagens sobre os métodos tradicionais

de estimação, como a realização de uma busca global em uma região de busca pré-

definida, sem a necessidade de uma estimativa inicial dos parâmetros e sem a

diferenciação da função objetivo ou do modelo. Assim, estes métodos são mais

robustos, quando comparados aos métodos tradicionais, já que os métodos heurísticos

utilizam somente os valores da função objetivo e não são influenciados pela presença de

parâmetros não-significativos. O método do Enxame de Partículas, em particular,

apresentou resultados superiores aos demais algoritmos heurísticos utilizados neste

trabalho, tanto na minimização de funções-teste, como na estimação de parâmetros de

problemas reais. Também foi possível mostrar que o grande número de avaliações da

função objetivo realizadas pelos métodos heurísticos durante a minimização pode ser

utilizado na construção das regiões de confiança dos parâmetros, sem a necessidade de

aproximações, garantindo uma analise estatística rigorosa dos resultados obtidos.

Page 6: Marcio Schwa Ab

vi

Abstract of Dissertation presented to COPPE/UFRJ as a partial fulfillment of the

requirements for the degree of Master of Science (M.Sc.)

EVALUATION OF HEURISTIC OPTIMIZATION ALGORITHMS

IN PARAMETER ESTIMATION PROBLEMS

Marcio Schwaab

September/2005

Advisors: José Carlos Costa da Silva Pinto

José Luiz Fontes Monteiro.

Department: Chemical Engineering

This work evaluates the use of heuristic optimization algorithms in nonlinear

parameter estimation problems. The heuristic methods used here were the Monte Carlo,

the Genetic Algorithm, the Simulated Annealing and the Particle Swarm algorithms.

This methods show several advantages over traditional parameter estimation methods,

as the realization of a global minimization in a predefined search region, without the

need for initialization parameter guesses and without differentiation of the objective

function or the model. For these reasons, these methods are more robust, when

compared to traditional methods, since the heuristic methods make use only of the

objective function values and are not influenced by presence of non significant

parameters. The Particle Swarm method, in particular, showed superior results over

others heuristic optimization methods used in this work, both for the minimization of

test functions and for estimation of kinetic parameters in kinetic models. It is also

shown that the large number of objective function evaluations performed by heuristic

methods during minimization can be used for construction of confidence regions of

parameter estimates, without need of approximations, ensuring a more rigorous

statistical analysis of the final obtained results.

Page 7: Marcio Schwa Ab

vii

ÍNDICE GERAL

1. INTRODUÇÃO .................................................................................................................1

2. REVISÃO BIBLIOGRÁFICA.............................................................................................5

2.1. A Função Objetivo.......................................................................................6

2.2. A Minimização...........................................................................................11

2.3. A Interpretação Estatística dos Resultados ..............................................27

3. OS ALGORITMOS HEURÍSTICOS .................................................................................35

3.1. Monte Carlo ...............................................................................................39

3.2. Algoritmo Genético....................................................................................49

3.3. Recozimento Simulado ..............................................................................62

3.4. Enxame de Partículas................................................................................71

3.5. Conclusões .................................................................................................85

4. PROBLEMAS DE ESTIMAÇÃO DE PARÂMETROS..........................................................88

4.1. Determinação da Região de Confiança em Modelos Não Lineares........88

4.1.1. Modelo de Michaelis-Menten .......................................................89

4.1.2. Modelo Cinético de 1ª Ordem .......................................................94

4.1.3. Modelo de Exponencial Dupla .....................................................99

4.2. Exemplo: Cinética da Reação de Polimerização de Etileno ..................103

4.3. Exemplo: Cinética da Reação de Produção de Gás de Síntese .............117

5. CONCLUSÕES E SUGESTÕES ......................................................................................132

REFERÊNCIAS BIBLIOGRÁFICAS...................................................................................136

APÊNDICE ......................................................................................................................143

Page 8: Marcio Schwa Ab

1

CAPÍTULO 1: INTRODUÇÃO

A aplicação de modelos matemáticos para a representação e interpretação de

processos químicos, tanto em escala laboratorial quanto industrial, é uma atividade

muito útil e bem difundida no campo da engenharia química. Os modelos matemáticos

são compostos por equações algébricas e/ou diferenciais que relacionam as variáveis do

processo, de modo que seja possível obterem-se previsões baseadas na resposta do

modelo, como, por exemplo, prever a conversão de uma reação a uma dada temperatura

e composição na alimentação do reator. Estas previsões fornecidas pelo modelo

permitem a simulação, o projeto e a otimização do processo que o modelo representa.

Entretanto, para que as previsões realizadas a partir de um modelo possam ser utilizadas

com segurança, é necessário que estas previsões representem o processo real com

precisão.

Assim, durante o desenvolvimento e aperfeiçoamento de um determinado

modelo, a etapa de estimação de parâmetros assume uma importância fundamental, pois

é através desta etapa que o modelo é comparado a dados reais (ou experimentais),

permitindo a avaliação da adequação do modelo. Os parâmetros são variáveis que

surgem durante o desenvolvimento dos modelos e que não podem ser medidas

experimentalmente (por exemplo, o fator pré-exponencial e a energia de ativação de

uma reação na equação de Arrhenius) e cujos valores devem ser estimados de tal forma

que o modelo represente os dados experimentais com a maior precisão possível.

O procedimento de estimação de parâmetros consiste na minimização de uma

função objetivo através do ajuste dos parâmetros do modelo, sendo a função objetivo

uma medida da distância entre os dados experimentais e as predições do modelo. Uma

função objetivo muito utilizada é a função de mínimos quadrados, que mede a soma dos

desvios quadráticos do modelo em relação ao experimento. Porém, quando a diferença

na grandeza das variáveis medidas é significante, ou mesmo quando os erros associados

às medições experimentais são diferentes, a utilização da função de mínimos quadrados

se torna inadequada. Neste caso, a utilização da função de máxima verossimilhança, que

maximiza a probabilidade de serem observados os valores experimentais efetivamente

obtidos, deve ser preferida. Essa função leva em consideração os erros associados a

cada medida e, conseqüentemente, a diferença relativa da ordem de grandeza das

Page 9: Marcio Schwa Ab

2

variáveis. Além disso, a utilização da função de máxima verossimilhança permite uma

análise estatística rigorosa dos resultados da estimação.

Uma vez definida a função objetivo, o passo seguinte é a sua minimização

através do ajuste dos parâmetros para, assim, obter-se um modelo com a melhor

predição possível. Os métodos de minimização tradicionalmente utilizados na estimação

de parâmetros são métodos determinísticos, nos quais, a partir de uma estimativa inicial

dos parâmetros, busca-se o mínimo da função objetivo. Dentre os métodos

determinísticos, destacam-se os métodos de Newton, os quais fazem uso do Gradiente e

da matriz Hessiana da função objetivo. As principais qualidades destes métodos são a

rápida convergência e a alta precisão dos valores obtidos.

Entretanto, devido à freqüente não linearidade encontrada nos modelos

matemáticos utilizados na engenharia química, a minimização da função objetivo pode

apresentar uma série de dificuldades, dentre as quais destacam-se, entre outros, a

presença de mínimos locais, a alta correlação entre os parâmetros, o desconhecimento

de uma boa estimativa inicial dos parâmetros e a presença de parâmetros não

significativos.

Assim, a utilização dos métodos tradicionais fica muito prejudicada, pois a sua

eficiência depende principalmente de uma boa estimativa inicial. Por outro lado, estes

métodos não são capazes de realizar uma busca global e para que o ótimo global seja

encontrado é necessária, mais uma vez, uma boa estimativa inicial dos parâmetros.

Para superar estas dificuldades, os métodos heurísticos de otimização surgem

como uma alternativa eficiente e robusta. Estes métodos são caracterizados por um

caráter aleatório na busca do ótimo, um grande número de avaliações da função objetivo

e por realizarem uma busca global em toda a região de interesse. Além disso, estes

métodos não dependem de uma boa estimativa inicial dos parâmetros e não utilizam

derivadas (Gradiente e matriz Hessiana) durante a otimização, o que garante a robustez

do método, mesmo em problemas com parâmetros não significativos.

Dentre os métodos heurísticos, o método de Monte Carlo, o Algoritmo Genético,

o Recozimento Simulado e o Enxame de Partículas serão avaliados neste trabalho. O

método de Monte Carlo apresenta um alto caráter aleatório, sendo que a busca consiste

em extensivas tentativas, onde os valores dos parâmetros são sorteados aleatoriamente

dentro da região de interesse. O Algoritmo Genético procura simular a evolução dos

seres vivos, onde os indivíduos evoluem ao longo das iterações mediante alguma função

que mede a aptidão de cada indivíduo (neste caso, a própria função objetivo). O

Page 10: Marcio Schwa Ab

3

Recozimento Simulado faz uma analogia ao recozimento de metais onde, a partir de

uma temperatura inicial, o metal é resfriado e os seus átomos são rearranjados de forma

a minimizar a energia de estrutura cristalina. O Enxame de Partículas faz uma analogia

com o comportamento gregário de animais (pássaros, peixes, abelhas, etc.), onde as

partículas do enxame trocam informações entre si para encontrar o ótimo global da

função objetivo.

Uma desvantagem marcante destes métodos heurísticos é a necessidade de um

grande número de avaliações da função objetivo, o que aumenta muito o custo

computacional. Por outro lado, este grande número de avaliações da função objetivo

pode ser utilizado para a avaliação estatística dos resultados finais, permitindo, por

exemplo, determinar a região de confiança dos parâmetros sem a necessidade de

aproximações.

O Enxame de Partículas, em particular, se destaca por sua simplicidade,

eficiência e robustez. Uma característica interessante deste método é marcada por uma

busca de caráter global no início do procedimento que, ao longo das iterações, torna-se

local, quando ocorre a convergência final das partículas. Assim, no início é realizada

uma busca global, permitindo a localização de possíveis mínimos globais; em seguida, o

caráter local da busca permite aumentar a precisão do valor obtido. Esta característica,

além de aumentar a probabilidade de se encontrar o mínimo global, garante uma boa

precisão do valor obtido e uma boa exploração da região próxima ao mínimo,

possibilitando uma boa representação da região de confiança dos parâmetros através da

utilização das avaliações da função objetivo realizadas pelo método durante a

minimização.

O objetivo principal desta dissertação é investigar a utilização dos métodos

heurísticos de otimização, em particular do Enxame de Partículas, na estimação de

parâmetros de modelos comuns na engenharia química. Serão avaliados os resultados

obtidos na estimação de parâmetros de problemas de grande dimensão e com alta

correlação entre os parâmetros, além da determinação da região de confiança dos

parâmetros sem a necessidade de aproximações.

No Capítulo 2 é feita uma revisão da literatura sobre o procedimento de

estimação de parâmetros, a síntese da função de máxima verossimilhança, os principais

métodos de otimização, incluindo os métodos heurísticos, e a análise estatística dos

parâmetros estimados.

Page 11: Marcio Schwa Ab

4

No Capítulo 3 são apresentados os algoritmos dos métodos heurísticos: Monte

Carlo, Algoritmo Genético, Recozimento Simulado e Enxame de Partículas. Os

algoritmos são validados através de problemas padrões de otimização e os resultados

obtidos por cada um são avaliados e comparados. Também é avaliada a capacidade de

cada algoritmo em descrever a região próxima ao mínimo da função, característica

importante para a determinação da região de confiança dos parâmetros.

No Capítulo 4 o Enxame de Partículas é avaliado quanto à capacidade de

descrever a região de confiança dos parâmetros em problemas de estimação com

modelo não lineares de pequena dimensão. Em seguida são resolvidos dois problemas

reais de estimação de parâmetros: o primeiro consiste na estimação dos parâmetros

cinéticos da reação de polimerização de etileno; o segundo consiste na estimação dos

parâmetros cinéticos da reação de conversão do metano em gás de síntese através da

reação de reforma do metano com dióxido de carbono combinada com a oxidação

parcial.

No Capítulo 5 são apresentadas as conclusões finais sobres os resultados obtidos

por cada método heurístico de otimização, tanto para a minimização das funções teste

quanto para a estimação dos parâmetros dos problemas reais. É feita uma avaliação final

da região de confiança dos parâmetros obtida através da utilização dos pontos avaliados

pelos métodos heurísticos durante a minimização da função objetivo. Por fim, são feitas

sugestões para trabalhos futuros.

Page 12: Marcio Schwa Ab

5

CAPÍTULO 2: REVISÃO BIBLIOGRÁFICA

Os primeiros trabalhos onde o procedimento de estimação de parâmetros é

discutido foram apresentados por Legendre, em 1806, e Gauss, em 1809. Apesar de

Legendre ter sido o primeiro a publicar o uso do método dos mínimos quadrados, Gauss

reivindicou o uso deste método já em 1795, o que dá a Gauss o reconhecimento de ter

sido o primeiro a usar esta importante ferramenta para a estimação de parâmetros

(BECK e ARNOLD, 1977).

Apesar de Gauss ter dado fundamentos estatísticos ao método de mínimos

quadrados, mostrando que as estimativas obtidas por este método maximizam a

densidade de probabilidade para uma distribuição normal dos erros, antecipando desta

forma o método de máxima verossimilhança, Gauss e seus contemporâneos parecem

preferir justificativas puramente heurísticas para utilização do referido método (BARD,

1974).

Os trabalhos seguintes, apresentados ainda no século XIX, realizados pelo

próprio Gauss e ainda por Cauchy, Bienaymé, Chebyshev, Gram, Schmidt e outros, se

concentraram nos aspectos computacionais do ajuste linear por mínimos quadrados. No

inicio do século XX, o desenvolvimento de métodos estatísticos de estimação foi

objetivo de estudo de Karl Pearson e R. A. Fisher. É atribuído a R. A. Fisher, apesar de

baseado nas idéias de Gauss, o desenvolvimento do método de máxima

verossimilhança, incluindo a análise das propriedades da estimação como consistência,

eficiência e suficiência (SEAL, 1967; HIMMELBLAU, 1970; BARD, 1974).

O advento dos computadores possibilitou o uso de métodos numéricos para

estimação de parâmetros, já que este procedimento envolve a busca do mínimo de uma

função não linear. Muitos métodos que levam os nomes de Newton, Gauss e Cauchy, e

que eram conhecidos há muito tempo, puderam finalmente ser aplicados de forma

intensiva aos mais diferentes problemas práticos de estimação, nas mais diferentes áreas

(BARD, 1974).

Um fato que mostra o grande impacto do desenvolvimento computacional na

área de estimação de parâmetros é o aumento do número de artigos publicados, como

foi observado pelo número de trabalhos no banco de dados do portal “Web of Science”

contendo as palavras “parameter” e “estimation”. Há um aumento exponencial do

número de artigos na área de estimação de parâmetros, mesmo considerando o número

Page 13: Marcio Schwa Ab

6

relativo de artigos em relação ao total de artigos em cada período de 10 anos, conforme

a Tabela 1.

Tabela 1: Evolução dos artigos publicados com estimação de parâmetros.

Período Artigos selecionados

Artigos totais

Percentual dos Artigos, %

1945-1954 7 721487 0,001 1955-1964 40 1526202 0,003 1965-1974 284 3894294 0,007 1975-1984 653 7308761 0,009 1985-1994 4018 9149986 0,044 1995-2004 13283 11757040 0,113

A seguir serão abordados os aspectos relacionados à função objetivo, sua síntese

e características, os métodos tradicionais e heurísticos de minimização e os aspectos

relacionados à análise estatística dos resultados.

2.1 A Função Objetivo

Uma definição estatisticamente rigorosa da função objetivo é o ponto de

fundamental importância no procedimento de estimação de parâmetros. É a partir da

definição da função objetivo que as análises estatísticas dos resultados podem ser feitas

com o rigor necessário. Uma forma de introduzir uma interpretação estatística rigorosa

no procedimento de estimação de parâmetros é através da definição da função de

máxima verossimilhança.

Considerando que as variáveis experimentais ez (independentes e dependentes)

estão sujeitas a erros, estas podem ser vistas como variáveis aleatórias descritas por uma

certa função de densidade de probabilidade

( ); ,eeP ∗z z V (2.1)

que fornece a probabilidade de obterem-se os valores experimentais ze, dados os valores

reais desconhecidos z* e uma medida do erro experimental Ve (matriz de covariância

dos desvios experimentais). O método de máxima verossimilhança consiste em

maximizar esta probabilidade, dadas as restrições do modelo

( )=f z,θ 0 (2.2)

onde z é um vetor que contém variáveis independentes e dependentes, f é um vetor com

as equações do modelo e θ é um vetor dos parâmetros do modelo.

Page 14: Marcio Schwa Ab

7

Admitindo que o modelo é perfeito e que os experimentos são bem feitos, de

modo que seja razoável admitir que os resultados experimentalmente obtidos sejam os

mais prováveis, parece razoável variarem-se os parâmetros de forma a maximizar a

probabilidade de observarem-se os resultados experimentais obtidos (BARD, 1974).

Assim, os valores reais e desconhecidos z* podem ser considerados iguais aos valores

calculados a partir do modelo, já que este é considerado perfeito, representados por zm.

Por outro lado, a definição da função de densidade de probabilidade não é uma

tarefa trivial, pois para isso um conjunto muito grande de experimentos seria necessário.

Como a realização de experimentos é freqüentemente custosa e pode levar muito tempo,

a opção mais comum é admitir que os dados experimentais seguem uma distribuição

normal de probabilidade, haja visto que também é muito difícil mostrar que um

conjunto de dados é ou não normalmente distribuído (RATKOWSKY, 1990). Além

disto, esta hipótese é geralmente razoável por causa das características da distribuição

normal, que são (segundo BARD, 1974):

a) geralmente apresenta um comportamento próximo a muitas medidas

experimentais;

b) com o aumento do número de medidas, muitas distribuições aproximam-se da

distribuição normal;

c) dadas a média e a variância de um conjunto de medidas e usando os conceitos do

cálculo variacional, é possível mostrar que a distribuição normal é a que insere a

menor quantidade de informações extras ao problema;

d) a facilidade no tratamento matemático, que permite a definição das distribuições

χ2 e F, de fundamental importância para o estabelecimento de intervalos de

confiança dos parâmetros e para o teste da qualidade do ajuste do modelo aos

dados.

Assim, admitindo que os desvios experimentais são normalmente distribuídos, a

função de máxima verossimilhança pode ser escrita como:

( ) ( )( )

( ) ( )2

-12 1exp2det

NTe m e m

ee

Lπ −

= − − −

θ z z V z zV

(2.3)

onde N é o número total de variáveis. Maximizar a função acima é o mesmo que

minimizar a função

( ) ( ) ( )-1Te m e meS = − −θ z z V z z (2.4)

Page 15: Marcio Schwa Ab

8

Se os experimentos são realizados de forma independente, os termos da matriz

de covariância Ve que correlacionam as variáveis de diferentes experimentos são nulos

e, assim, pode-se reescrever a função objetivo como:

( ) ( ) ( )-1

1i

NE Te m e mi i e i i

i

S=

= − −∑θ z z V z z (2.5)

onde i denota o experimento e NE o número total de experimentos. Admitindo-se que os

desvios nas variáveis independentes e dependentes não estão correlacionados, pode-se

escrever a função objetivo da seguinte forma:

( ) ( )( ) ( )( ) ( ) ( )-1 -1

1 1

, ,i i

NE NET Te m m e m m e m e mi i i ey i i i i i ex i i

i i

S= =

= − − + − −∑ ∑θ y y x θ V y y x θ x x V x x (2.6)

onde xe e xm são vetores com os valores das variáveis independentes medidos

experimentalmente e otimizados, respectivamente, ye e ym são vetores com os valores

das variáveis dependentes medidos experimentalmente e calculados pelo modelo,

respectivamente, θ o vetor dos parâmetros do modelo e Vex e Vey a matriz de

covariância dos desvios experimentais relacionados, respectivamente, às variáveis

independentes e às variáveis dependentes.

É interessante observar que a Equação (2.6) mostra de forma clara que os

desvios em ambas variáveis podem ser considerados durante a estimação de parâmetros.

Entretanto, deve-se observar que os valores de ym são calculados a partir de um modelo,

dados os parâmetros θ e as variáveis independentes xm. Geralmente, os desvios nas

variáveis independentes são descartados, já que estas variáveis podem ser controladas

com certa precisão. Entretanto, em problemas de reconciliação de dados, principalmente

quando são utilizados dados de plantas industriais, as variáveis independentes podem

apresentar desvios da mesma ordem de grandeza presentes em de uma variável

dependente (como as vazões de correntes de entrada e saída de um reator). Assim, o

procedimento de estimação de parâmetros que leva em consideração os desvios nas

variáveis independentes, assume um papel fundamental na reconciliação de dados (KIM

et al., 1991a, 1991b). DOVI et al. (1993) desenvolveram um procedimento para o

planejamento seqüencial de experimentos especificamente para modelos onde as

variáveis independentes também estão sujeitas a desvios.

Page 16: Marcio Schwa Ab

9

Entretanto, em escala laboratorial é muito comum considerar que as variáveis

independentes são conhecidas com grande precisão. Assim, o segundo termo da direita

da Equação (2.6) é nulo, resultando na seguinte expressão:

( ) ( )( ) ( )( )-1

1

, ,i

NE Te m m e m mi i i ey i i i

i

S=

= − −∑θ y y x θ V y y x θ (2.7)

Por outro lado, se as medidas das variáveis dependentes são realizadas de forma

a não haver correlação entre elas; isto é, se os desvios entre as variáveis dependentes

não estão correlacionados, a função objetivo pode ser reescrita em uma forma conhecida

como função de mínimos quadrados ponderados:

( )( )( )22

1 1

,e m mNE NY ij ij i

i j ij

y yS

σ= =

−=∑∑

x θθ (2.8)

onde σ2 é a variância associada à medida da variável j no experimento i.

Admitindo que o erro das variáveis dependentes é constante; ou seja, que é o

mesmo para todas as variáveis dependentes e em todos os experimentos, chega-se à

conhecida função de mínimos quadrados:

( ) ( )( )21 1

,NE NY

e m mij ij i

i j

S y y= =

= −∑∑θ x θ (2.9)

Este resultado assegura uma interpretação estatística ao método de mínimos

quadrados, conforme descrito por Gauss. Entretanto, ao longo do desenvolvimento

acima realizado foram realizadas inúmeras hipóteses, as quais devem ser satisfeitas para

que o método de mínimos quadrados tenha tal interpretação estatística. Além disso, em

modelos com mais de uma saída (por exemplo, um modelo que prevê a conversão e a

temperatura na saída de um reator) a função de mínimos quadrados é inadequada devido

às diferentes magnitudes das respostas do modelo.

É importante ressaltar a importância da matriz de covariância das medidas

experimentais. O desenvolvimento das equações feito acima foi sempre baseado em

hipóteses com relação à estrutura desta matriz. Assim, a escolha adequada da função

objetivo passa pelo conhecimento prévio da estrutura dos erros experimentais do

problema em questão. Além disso, em certos casos a matriz de erros pode trazer

importantes informações sobre o problema em questão. Por exemplo, no trabalho de

Page 17: Marcio Schwa Ab

10

LARENTIS et al. (2003) a análise da matriz dos erros experimentais possibilitou a

obtenção de informações sobre o mecanismo da reação estudada.

Freqüentemente a matriz de covariância dos erros experimentais é totalmente ou

parcialmente desconhecida, o que é particularmente verdade quando os dados são

correlacionados. O método da máxima verossimilhança permite que estes componentes

desconhecidos sejam considerados como parâmetros adicionais e estimados em

conjunto com os parâmetros do modelo. Entretanto, quando o número de variáveis é

grande e, conseqüentemente, o número de componentes desconhecidos da matriz de

covariância também é grande, a dimensão do problema de estimação aumenta

consideravelmente, tornando esta prática não usual devido a dados experimentais

limitados, erros sistemáticos de modelagem e outras violações de hipótese assumidas

para a análise de máxima-verossimilhança (RICKER, 1984). Além disso, SANTOS e

PINTO (1998) mostraram que a estimação simultânea dos parâmetros e da correlação

entre as variáveis pode levar a parâmetros sem quaisquer significados físico e/ou

matemático, além de inserir dificuldades adicionais para minimização da função

objetivo.

Um ponto importante para a estimação de parâmetros é a classificação das

variáveis em independentes e dependentes. Dependendo da estrutura matemática do

modelo, pode ser mais conveniente obter x como função de y; ou seja, calcular a

variável independente como função da variável dependente, fato comum em modelos

cinéticos. Assim, alguns autores estimam os parâmetros de tal modelo minimizando os

desvios em x, dados os valores de y (FROMENT e MEZAKI, 1970). Entretanto,

MEZAKI et al. (1973) mostram que a troca de y por x na minimização dos desvios

quadráticos para estimação dos parâmetros, além de ser teoricamente errado, pode levar

a uma estimativa ruim dos parâmetros. Isso ocorre porque a variável x é geralmente

conhecida com grande precisão, já que esta é a condição que foi fixada para realização

do experimento, o que não ocorre com a variável y, que é o resultado do experimento e

que pode depender de uma série de fatores nem sempre controláveis.

Outro ponto muito importante e geralmente ignorado é a hipótese de que o

modelo é perfeito. Os modelos utilizados na engenharia química, mesmo quando

baseados em teorias, estão sujeitos a muitas simplificações, pois é virtualmente

impossível desenvolver um modelo matemático onde todos os fenômenos que ocorrem

em um processo químico sejam considerados simultaneamente. Por outro lado, mesmo

se a estrutura de um modelo é perfeita, os parâmetros deste modelo são estimados a

Page 18: Marcio Schwa Ab

11

partir de dados experimentais incertos, e, assim, as predições do modelo também são

incertas, de modo que o modelo não é perfeito. Por este motivo, os desvios entre as

predições do modelo e os dados experimentais devem ser descritos pela distribuição do

erro de predição, ou melhor, pela matriz de covariância de predição, que considera tanto

os erros experimentais quanto os erros de modelagem.

Baseado nestes fatos, SANTOS e PINTO (1998) propuseram um procedimento

de estimação onde a matriz de covariância dos desvios experimentais é substituída pela

matriz de covariância de predição. Como esta matriz não é conhecida a principio, já que

depende dos parâmetros estimados, um procedimento iterativo deve ser utilizado. A

matriz de covariância de predição é calculada ao final da minimização e a função

objetivo, com a nova matriz de predição, é minimizada novamente. Este processo

iterativo é repetido até que a matriz de predição atinja a convergência.

Além do método de máxima verossimilhança, existem outros métodos de

estimação que diferem na formulação da função objetivo, como a estimação Bayesiana

(BOX e DRAPER, 1965), o método dos momentos (GALLANT, 1987), além de outros

métodos, mas não serão considerados neste trabalho.

2.2 A Minimização

Uma vez definida a função objetivo, a etapa seguinte é a minimização desta

função. A minimização consiste na busca de um conjunto de parâmetros onde a função

atinja o menor valor possível. Nos casos onde as variáveis independentes também estão

sujeitas a erros, a minimização é realizada de forma a encontrar os valores ótimos dos

parâmetros e das variáveis independentes.

A princípio, a minimização da função objetivo pode ser feita por qualquer

método de otimização e o seu resultado não deveria influenciar no procedimento de

estimação, já que o mínimo da função objetivo independe do método usado para

encontrá-lo. Porém, fatos como a freqüente não linearidade dos modelos, a presença de

mínimos locais e a alta correlação entre os parâmetros podem tornam a minimização da

função objetivo uma tarefa nada fácil. Assim, a escolha adequada do método pode ser

determinante para o sucesso do procedimento de estimação de parâmetros.

Existem diversos métodos de otimização, os quais podem ser classificados em

(HIMMELBLAU, 1972):

a) Métodos Analíticos: consistem no uso de técnicas clássicas do cálculo diferencial e

cálculo das variações, buscando os pontos onde as derivadas são nulas; entretanto,

Page 19: Marcio Schwa Ab

12

para sistemas não lineares e de grande dimensão, estes métodos são

insatisfatórios;

b) Métodos Numéricos: consistem no uso de informações prévias para a otimização

do problema através de procedimentos iterativos; são métodos capazes de resolver

problemas não lineares e de grande dimensão;

c) Métodos Gráficos: consistem no uso de uma representação gráfica da função e da

determinação do ponto ótimo por inspeção direta; é um método elementar e tem

como vantagem a possibilidade de mostrar se o ponto ótimo existe ou não;

entretanto, é restrito a problemas de uma ou duas dimensões;

d) Métodos Experimentais: consistem na busca de um extremo através da

experimentação direta sobre as variáveis do processo ao invés da manipulação

matemática da função que descreve tal processo; este método consegue tratar de

problemas onde a representação matemática de determinado processo é difícil,

sendo estas técnicas conhecidas como EVOP (Evolutionary Optimization);

e) Métodos de Estudo de Caso: consistem na avaliação de um número representativo

de soluções; a melhor solução encontrada é admitida como o ótimo.

Na estimação de parâmetros, a implementação de métodos analíticos é viável no

caso onde a função objetivo não considera os desvios nas variáveis independentes e

quando o modelo é linear nos parâmetros. Um modelo é dito linear nos parâmetros

quando pode ser escrito na forma

( ) ( ) ( )1

,NP

Ti i

i

f θ=

= ⋅ = ⋅∑y x θ f x θ x (2.10)

Assim, considerando a função objetivo descrita pela Equação (2.7), o conjunto

de parâmetros que minimiza esta função objetivo, dado o modelo linear nos parâmetros

(Equação (2.10)), pode ser calculado pela seguinte equação:

( ) ( ) ( )1

-1 -1

1 1i i

NE NET T e

i ey i i ey ii i

= =

= ∑ ∑θ f x V f x f x V y (2.11)

Apesar de fortes restrições impostas ao resultado acima, um número muito

grande de trabalhos é devotado a este problema clássico de modelagem. Por isso,

mesmo que de modo aproximado, grande parte da análise estatística dos resultados que

são encontrados na literatura são baseados na aproximação linear dos modelos.

Page 20: Marcio Schwa Ab

13

DRAPER e SMITH (1988) fornecem uma descrição detalhada sobre a estimação de

parâmetros em modelos lineares.

Entretanto, os modelos geralmente usados na engenharia química são não-

lineares. Por isso, além do fato da minimização da função objetivo não ter uma solução

analítica, a análise estatística dos resultados baseada na aproximação linear dos modelos

está também sujeita à qualidade desta aproximação, o que pode levar a conclusões

equivocadas, conforme será discutido no próximo item.

Os métodos gráficos são muito restritos e os métodos experimentais são focados

para a busca de ótimos operacionais. Os métodos de estudo de caso, por sua vez, são

computacionalmente intensivos; contudo, com o aumento da capacidade de cálculo dos

computadores modernos, esses métodos constituem alternativas interessantes. A

descrição destes métodos remete aos métodos heurísticos de otimização, apesar de não

serem muito difundidos (ou mesmo nem existirem) naquela época.

Os métodos numéricos utilizam procedimentos iterativos para a otimização e são

tradicionalmente utilizados para a minimização da função objetivo da estimação de

parâmetros. Estes métodos podem ser divididos em duas classes distintas: a primeira

consiste nos métodos que utilizam o cálculo das derivadas da função objetivo para a

minimização; e a segunda, as técnicas em que a busca é feita sem o cálculo de

derivadas.

Dentre os métodos que usam derivadas, o método do gradiente é o mais simples

e utiliza somente o vetor gradiente da função objetivo (∇S); ou seja, as derivadas de

primeira ordem (HIMMELBLAU, 1972). A partir de um ponto inicial é definida a

direção de busca, já que o vetor gradiente indica a direção ao longo da qual a função

tem a maior variação. A equação a seguir mostra a forma iterativa com que a busca é

feita:

k k kλ+ = − ∇ θθ θ S (2.12)

onde k indica a iteração, θ é o vetor dos parâmetros, ∇S é o vetor gradiente e λ é o passo

que é dado na direção definida pelo vetor gradiente. λ pode ser um valor fixo ou

variável ao longo das iterações ou até mesmo ser otimizado em cada iteração. A busca

segue até que a variação nos valores dos parâmetros seja pequena ou quando o vetor

gradiente seja praticamente nulo (de acordo com a precisão desejada). A eficiência deste

método está intimamente ligada ao controle do passo λ e com a otimização deste valor

Page 21: Marcio Schwa Ab

14

em cada iteração, apesar de interessante, pode aumentar muito o custo computacional

deste método.

Já o método de Newton utiliza, além do vetor gradiente, a matriz Hessiana da

função objetivo (H); ou seja, as derivadas de segunda ordem (HIMMELBLAU, 1972).

Este método é baseado na aproximação quadrática da função objetivo e o procedimento

iterativo é descrito pela seguinte equação:

1 1ˆ ˆ

k k kλ+ −= − ∇θ θθ θ H S (2.13)

Na Equação (2.13), o produto entre a inversa da matriz Hessiana e o vetor

gradiente fornece a direção e o tamanho do passo que deve ser dado nesta direção,

baseado na aproximação quadrática da função objetivo. O passo λ é mantido na equação

recursiva para acelerar ou retardar a busca, de acordo com a qualidade da aproximação.

A grande vantagem destes métodos é a rapidez na convergência, já que uma

grande quantidade de informação sobre a função é utilizada. Por outro lado, obter

expressões analíticas do gradiente e da matriz Hessiana pode tornar-se inviável em

problemas de grande dimensão. Por outro lado, o cálculo das derivadas através de

aproximações numéricas envolve um grande número de avaliações da função, além de

não serem os valores exatos das derivadas.

Gauss propôs uma aproximação para o cálculo da matriz Hessiana para

problemas de estimação de parâmetros. Considerando como função objetivo a Equação

(2.7), ou seja, sem erros na variável independente e sem correlação entre experimentos,

a matriz Hessiana pode ser descrita por:

( )( )22

1 1

1 1

2 , 2Tm m mNE NETe m m i i i

i i i ey eyi i

S y y yH y y x θ V Vθ θ θ θ θ θ

− −

= =

∂ ∂ ∂∂ = =− − + ∂ ∂ ∂ ∂ ∂ ∂ ∑ ∑ (2.14)

A aproximação consiste em desprezar o primeiro termo da direita da Equação

(2.14), considerando a hipótese de que os desvios entre os valores experimentais e a

predição do modelo são pequenos. Assim, a matriz Hessiana pode ser aproximada por:

1

1

2Tm mNE

i iey

i

y yH Vθ θ

=

∂ ∂ ≈ ∂ ∂ ∑ (2.15)

Geralmente, quando a aproximação de Gauss é usada num procedimento de

estimação de parâmetros, o método é chamado de Gauss-Newton. Um aspecto

importante da aproximação de Gauss é a possibilidade de calcular a matriz Hessiana

Page 22: Marcio Schwa Ab

15

somente com valores das derivadas de primeira ordem do modelo, as quais já são

utilizadas para o cálculo do vetor gradiente da função objetivo. É interessante observar

que, para modelos lineares, a aproximação de Gauss é exata, já que o primeiro termo da

direita da Equação (2.14) é nulo, pois a derivada segunda de um modelo linear em

relação a seus parâmetros é nula. Devido a esta característica, KOSTINA (2004)

compara o método de Gauss-Newton a um método de programação linear seqüencial, já

que o termo relativo às derivadas de segunda ordem do modelo é ignorado.

Outro ponto importante na aplicação dos métodos de Newton é a inversão da

matriz Hessiana. Além de acarretar um custo computacional relativamente alto, a sua

inversão pode ser difícil ou mesmo não ser possível em problemas onde existam

parâmetros pouco significativos. Assim, no caso de uma minimização, é importante

garantir que a matriz Hessiana seja positiva definida. O método de Levenberg-

Marquardt (LEVENBERG, 1944; MARQUARDT, 1963) propõe uma forma prática e

eficiente de garantir a inversão da matriz Hessiana, de forma que a equação recursiva do

método de Newton (Equação (2.13)) seja reescrita como:

11

ˆ ˆk k k kλ β

−+ = − + ∇ θ θθ θ H I S (2.16)

onde I é a matriz identidade e βk é uma constante definida de forma a garantir que a

soma das matrizes seja positiva definida.

Uma forma interessante de manipular o valor de β é diminuir seu valor ao longo

das iterações, já que próximo ao mínimo espera-se que a matriz Hessiana seja positiva

definida. Assim, o procedimento se inicia com um comportamento próximo ao do

método do gradiente e, ao longo das iterações, aproxima-se do comportamento de um

método de Newton. Comparando diversos métodos que utilizam derivadas, BARD

(1970) mostra que o método de Levenberg-Marquardt é inferior aos métodos de Gauss-

Newton. Mais tarde, porém, BARD (1974) conclui que em problemas onde não se

dispõe de boas estimativas iniciais, o método de Levenberg-Marquardt pode ser mais

eficiente.

Uma forma interessante de evitar o cômputo da matriz Hessiana é proposta pelo

método Davidon-Fletcher-Powell (FLETCHER e POWELL, 1963), onde uma

aproximação da inversa da matriz Hessiana na iteração k+1 é calculada através das

informações coletadas até a iteração k.

Page 23: Marcio Schwa Ab

16

BRITT e LUECKE (1973) desenvolvem um método de estimação de parâmetros

incluindo problemas onde ambas variáveis, independentes e dependentes, estão sujeitas

a desvios, de acordo com a função de máxima verossimilhança. ANDERSON et al.

(1978) chegam a um método muito similar ao proposto por BRITT e LUECKE (1973),

já que ambos fazem uma aproximação linear do modelo em cada iteração (ambos os

métodos são do tipo de Gauss-Newton). Entretanto, ANDERSON et al. (1978) aplicam

um método descrito por LAW e BAILEY (1963), que limita a busca sempre que a

aproximação linear do modelo não é satisfatória. Isso torna o procedimento de

estimação mais robusto com relação à estimativa inicial dos parâmetros.

Como foi discutido acima, as vantagens dos métodos de Newton são a rapidez

na busca pelo mínimo, devido à convergência quadrática destes métodos, e a garantia de

que um ponto mínimo é encontrado, já que a busca é terminada quando o vetor

gradiente é nulo. Por outro lado, a necessidade do cálculo das derivadas (analíticas ou

numéricas), a inversão da matriz Hessiana e as dificuldades associadas à necessidade de

uma estimativa inicial dos parâmetros podem tornar inviável a utilização destes

métodos.

Uma forma de evitar algumas destas dificuldades pode ser feita através da

utilização dos métodos de busca direta; isto é, sem a necessidade do cálculo das

derivadas. Estes métodos são simples, sendo a direção e o sentido da busca determinada

através da inferência direta dos valores da função. Assim, a partir de um ponto inicial,

são feitas tentativas nas vizinhanças deste ponto. De acordo com os valores obtidos,

pode-se determinar uma direção para onde o ponto ótimo deve estar localizado. Uma

vez determinada a direção de busca, podem ser usados diversos critérios para que a

busca avance nesta direção, como dar passos nesta direção até que ocorra uma falha;

isto é, que um valor indesejável seja encontrado. Uma vez feito isso, uma nova direção

de busca deve ser determinada e o procedimento segue até que o ponto ótimo ou algum

critério de parada seja satisfeito.

A principal vantagem destes métodos é o fato de não utilizar o cômputo das

derivadas da função, tornando estes métodos extremamente simples e com iterações

muito rápidas. Com relação aos métodos de Newton, os métodos de busca direta são

mais robustos e o sucesso da busca é menos dependente das estimativas iniciais.

Entretanto, o fato de não usar as derivadas da função faz com que a busca seja baseada

em uma quantidade pequena de informações sobre a função, podendo dificultar a busca

do ponto ótimo. Por isso, geralmente, são necessárias muitas iterações para que o ponto

Page 24: Marcio Schwa Ab

17

ótimo seja alcançado. Os métodos que utilizam este tipo de procedimento são os

métodos de Hooke-Jeeves (HOOKE e JEEVES, 1961), Rosenbrock, (ROSENBROCK,

1960), Powell (POWELL, 1964 e 1965) e Simplex, também conhecido como Poliedros

Flexíveis (NELDER e MEAD, 1965). BARD (1974) conclui que, apesar de ser atraente

a idéia de realizar uma busca direta sem o cômputo de derivadas, a eficiência destes

métodos é baixa, quando comparada a dos métodos que usam o cômputo de derivadas.

Entretanto, devido à freqüente não linearidade dos modelos na área da

engenharia química, a função objetivo usada na estimação de parâmetros pode

apresentar mínimos locais e os métodos de otimização descritos acima são incapazes de

distinguir um ótimo local de um ótimo global. Considerando ainda as dificuldades

associadas ao cálculo das derivadas e inversão da matriz Hessiana, a utilização destes

métodos pode apresentar sérias dificuldades, muitas vezes, impossibilitando sua

aplicação.

As dificuldades acima estão intimamente relacionadas com a definição

apropriada de uma boa estimativa inicial dos parâmetros. Porém, nem sempre se dispõe

de uma boa estimativa inicial e em problemas de grande dimensão, onde a correlação

entre os parâmetros é geralmente alta, é virtualmente impossível dispor de uma boa

estimativa inicial.

Uma dificuldade adicional que a estimação de parâmetros pode apresentar é a

existência de parâmetros não significativos. Este problema provoca dificuldades

adicionais ao procedimento de estimação de parâmetros e pode ocorrer por causa da

superparametrização do modelo, à própria estrutura do modelo ou aos pontos

experimentais disponíveis, os quais podem não permitir a identificação de um efeito

específico.

Tendo em vista as dificuldades acima relatadas, os algoritmos heurísticos

surgem como uma solução eficiente, robusta e de simples implementação

computacional, tornando estes métodos muito atraentes. Se a classificação dos métodos

de otimização descrita por HIMMELBLAU (1972) for considerada, os métodos

heurísticos de otimização se enquadrariam nos métodos de Estudos de Caso. Entretanto,

esta definição não representa de forma correta estes métodos e os termos heurísticos ou

estocásticos são mais apropriados.

Os algoritmos heurísticos de otimização são caracterizados pela realização de

um grande número de avaliações da função objetivo em toda a região de busca, de

forma a aumentar a probabilidade de encontrar o ótimo global da função objetivo. Além

Page 25: Marcio Schwa Ab

18

disso, o caráter aleatório do procedimento de busca é elevado, para evitar que a busca

fique presa a um ótimo local. Estes métodos não precisam de uma estimativa inicial e

não utilizam as derivadas para chegar ao ponto ótimo, evitando assim muitas das

dificuldades associadas aos métodos tradicionais.

Dentre os métodos heurísticos, destacam-se o método de Monte Carlo, o

Algoritmo Genético (Genetic Algorithm), o Recozimento Simulado (Simulated

Annealing) e o Enxame de Partículas (Particle Swarm). Ainda existem outros métodos

heurísticos, como o método da Busca Tabu (Tabu Search) e o método da Colônia de

Formigas (Ant Colony), mas estes não serão considerados neste trabalho.

O método de Monte Carlo, também chamado de Busca Aleatória Pura (DIXON

e SZEGO, 1975), consiste em um procedimento totalmente aleatório de busca pelo

ponto ótimo de uma função qualquer. A busca consiste em sortear aleatoriamente um

número elevado de pontos na região de busca e considerar o melhor ponto encontrado

como o ponto ótimo, sendo esta a forma mais simples do método de Monte Carlo. É um

método extremamente simples e robusto, pois é facilmente aplicado e sempre funciona,

independente de chute inicial, do cálculo de derivadas e da dimensão do sistema. É claro

que a grande limitação deste método é a necessidade de um número muito grande de

avaliações da função objetivo; não existindo garantia formal de que o ponto ótimo tenha

sido encontrado.

Uma forma de tentar melhorar a eficiência deste método é identificar as regiões

mais prováveis em conter o ótimo global e concentrar a busca nestas regiões, sendo esta

forma conhecida como por Busca Adaptativa Pura (PATEL et al., 1988; ZABINSKY e

SMITH, 1992). PRICE (1976) propôs um algoritmo híbrido que combina a Busca

Aleatória Pura com o método Simplex, sendo este novo algoritmo chamado de Busca

Aleatória Controlada. KRIVY e TVRDÍK (1995) inseriram uma pequena modificação

no algoritmo de Price e aplicaram este método em problemas de regressão. KLEPPER e

HENDRIX (1994a) propuseram um novo método de busca aleatória chamado de UCPR

(“uniform covering by probabilistic rejection”) e aplicaram este algoritmo na estimação

e determinação da região de confiança de parâmetros de modelos não lineares.

BROOKS (1958) discutiu a aplicação destes métodos aleatórios no planejamento

de experimentos para busca de pontos ótimos, alegando que planejamentos fatoriais em

problemas com muitas variáveis levam a um número muito grande de experimentos.

Além disso, Brooks discutiu duas formas alternativas dos métodos aleatórios: a primeira

consiste em dividir a região de busca em sub-regiões e sortear um ponto aleatório em

Page 26: Marcio Schwa Ab

19

cada sub-região; a segunda alternativa consiste em realizar tentativas aleatórias em

torno de uma estimativa inicial do ponto ótimo, sorteando, por exemplo, pontos

normalmente distribuídos em torno desta estimativa inicial. De forma surpreendente,

Brooks observou que existem analogias intrigantes entre a segunda forma dos métodos

aleatórios e a evolução, aparentemente antecipando já em 1958 os algoritmos evolutivos

de otimização.

A grande aplicação do método de Monte Carlo está relacionada com a geração

de um grande número de resultados estocásticos para posterior análise estatística. Um

exemplo é o trabalho de DONALDSON e SCHNABEL (1987), onde o método de

Monte Carlo é aplicado para a geração de um grande conjunto de dados para posterior

determinação da região de confiança de parâmetros.

O Algoritmo Genético é um método de otimização baseado na evolução dos

seres vivos, onde os indivíduos mais aptos tendem a sobreviver (teoria da seleção

natural de Darwin) e uma herança genética é passada através do cruzamento dos

indivíduos (princípios da herança genética de Mendel). Este método foi desenvolvido

por John Holland e seus colegas da Universidade de Michigan e o primeiro trabalho

neste tópico á a monografia de John Holland de 1975, como descrito por GOLDBERG

(1989).

A analogia com um procedimento de otimização é baseada na avaliação da

aptidão de um individuo através de uma função objetivo. Cada indivíduo comporta um

conjunto de valores das variáveis de otimização que é considerado a carga genética

deste indivíduo. Durante o cruzamento entre indivíduos pré-selecionados, ocorre a troca

de informações, isto é, a carga genética dos indivíduos é combinada. Assim, novos

indivíduos são gerados a partir do cruzamento entre os indivíduos selecionados do

grupo, de acordo com o valor da função objetivo; isto é, indivíduos mais aptos têm

maior chance de serem selecionados e realizarem o cruzamento. Por outro lado,

ocasionalmente alguns indivíduos sofrem mutações; ou seja, sua carga genética é

alterada aleatoriamente, gerando um novo indivíduo com características diferentes. A

mutação ocorre com igual probabilidade em qualquer individuo do grupo, independente

de sua aptidão, e possibilita uma maior diversificação no processo de busca.

Tradicionalmente, a evolução da busca era creditada ao cruzamento, enquanto a

mutação evitava a convergência prematura. Porém, estudos recentes defendem a

mutação com sendo o verdadeiro responsável pela evolução de todo grupo (HIBBERT,

1993a; BALLAND et al., 2000).

Page 27: Marcio Schwa Ab

20

Na sua forma tradicional, o Algoritmo Genético é usado com codificação

binária; porém, em diversas aplicações esta codificação não vem mais sendo usada,

evitando assim a necessidade de converter os valores das variáveis de otimização da

codificação binária para codificação decimal, sempre que é necessário avaliar a função

objetivo (BALLAND et al., 2000).

Uma desvantagem deste método é que, conforme se aproxima do valor ótimo, a

convergência do procedimento se torna muito lento. Uma alternativa interessante é

utilizar o Algoritmo Genético para gerar uma solução aproximada e refinar esta solução

com um método de busca local, os chamados algoritmos híbridos. HIBBERT (1993b)

avaliou a utilização do Algoritmo Genético seguido de um método de gradiente,

comparando algumas formas de usar os resultados obtidos pelo Algoritmo Genético

como ponto de partida do método de gradiente. MOROS et al. (1996) e PARK e

FROMENT (1998) utilizaram o Algoritmo Genético para geração de estimativas iniciais

para o método de Levenberg-Marquardt em problemas de estimação de parâmetros de

modelos cinéticos. BALLAND et al. (2000) utilizaram o Algoritmo Genético seguido

do método de Rosenbrock, também para estimação de parâmetros. Em todos estes casos

a aplicação de um algoritmo híbrido permitiu a solução de problemas de grande

dimensão, garantindo tanto a precisão do valor ótimo encontrado como a realização de

uma busca global.

O método do Recozimento Simulado, proposto por KIRKPATRICK et al.

(1983) para problemas de otimização combinatória, baseia-se na minimização da

energia da estrutura cristalina sofrida por um corpo metálico durante o processo de

recozimento. Este processo consiste no aquecimento de um corpo metálico até uma

temperatura suficientemente alta. Em seguida, é resfriado lentamente, possibilitando a

reorganização dos átomos do corpo metálico em uma estrutura com menor energia.

O procedimento de otimização de uma função objetivo consiste em explorar a

região de busca através de transições sucessivas (análogas às mudanças da estrutura

cristalina). As transições consistem em pequenas mudanças das variáveis de otimização.

De acordo com o valor da função objetivo neste novo ponto, a transição pode ou não ser

aceita: se a transição para o novo ponto leva a um menor valor da função objetivo, a

mudança é aceita; se o valor da função objetivo no novo ponto for maior, a transição

pode ser aceita ou não, de acordo com o cálculo de uma probabilidade de transição.

A forma tradicional do cálculo da probabilidade de transição é feita de acordo

com a estatística de Boltzmann, que corresponde à seguinte equação:

Page 28: Marcio Schwa Ab

21

expN

k

F FpT

− = − (2.17)

onde p é a probabilidade de transição, F é o valor da função objetivo, FN é o valor da

função objetivo no novo ponto e Tk é a k-ésima temperatura. O valor calculado de p é

comparado com um número aleatório com distribuição uniforme no intervalo [0, 1]. Se

p é maior que este número aleatório, a transição é aceita; caso contrário, a transição não

ocorre. O valor de T deve ser alto no início, aumentando a probabilidade de ocorrerem

transições para pontos onde a função objetivo aumenta, possibilitando que ocorra uma

boa exploração do espaço de busca e levando o algoritmo a encontrar o mínimo global.

KIRKPATRICK et al. (1983) diminuem exponencialmente o valor de T, de acordo com

a seguinte equação:

1k kT Tα+ = ⋅ (2.18)

sendo α um valor positivo menor que 1.

TSALLIS e STARIOLO (1996) sugeriram um algoritmo onde a probabilidade

de transição é calculada de acordo com a estatística de Tsallis, definida pela seguinte

equação:

( )( )

11

1 1aN q

ak

F Fp qT

− − = − − (2.19)

onde qa é um parâmetro adicional. Neste caso, a diminuição do parâmetro temperatura é

feita da seguinte forma:

( )

1

1 12 1

1 1

v

v

q

k k qT Tt

+ −−= ⋅

+ − (2.20)

onde qv é outro parâmetro adicional e t é o tempo discreto que corresponde ao número

de pontos avaliados até o momento.

AHON et al. (2000) compararam diferentes algoritmos do Recozimento

Simulado. Com relação à forma de calcular a probabilidade de transição, concluíram

que, em problemas de grande dimensão, o algoritmo baseado na estatística de Tsallis

converge para o mínimo global com maior freqüência que o algoritmo baseado na

estatística de Boltzmann. Para isso o algoritmo utiliza um número maior de avaliações,

Page 29: Marcio Schwa Ab

22

tornando-se assim mais lento. Uma dificuldade adicional do uso da estatística de Tsallis

é a determinação dos valores dos parâmetros adicionais qa e qv.

Outro ponto importante para a aplicação do método do Recozimento Simulado é

a definição do valor inicial do parâmetro temperatura. DAS et al. (1990) determinaram

o parâmetro temperatura da seguinte forma:

0 lnFTp

∗∆=− (2.21)

onde ∆F* é a diferença entre o maior e o menor valor da função objetivo em um

pequeno número de avaliações da função objetivo e p é a probabilidade das transições

ocorrerem, que deve ser alta para que no início da busca muitas transições sejam aceitas,

possibilitando a exploração da região de busca.

Um algoritmo híbrido proposto por KVASNICKA e POSPÍCHAL (1997) utiliza

o Recozimento Simulado para tornar aleatória uma busca realizada com base no método

Simplex. COSTA et al. (2000) aplicaram o método do Recozimeto Simulado em um

problema de estimação de parâmetros de modelos termodinâmicos e mostraram que os

resultados obtidos pelo método do Recozimento Simulado são superiores aos resultados

obtidos pelo método de Powell, já que o método do Recozimento Simulado foi capaz de

encontrar mínimos menores que o método de Powell, que é um método de busca local.

EFTAXIAS et al. (2002) compararam o desempenho do Recozimento Simulado com o

método de Levenberg-Marquardt em problemas de estimação de parâmetros de modelos

cinéticos. Mostraram que, em problemas de grande dimensão, somente o Recozimento

Simulado consegue estimar os parâmetros, compensando o alto custo computacional

requerido por este método.

O método do Enxame de Partículas consiste em um algoritmo de otimização

heurística inspirado no comportamento gregário de animais (peixes, pássaros, etc.).

Proposto por KENNEDY e EBERHART (1995), este método consiste na otimização de

uma função objetivo através da troca de informações entre os elementos (partículas) do

grupo, resultando em um algoritmo eficiente, robusto e de simples implementação

computacional.

O movimento de cada partícula em cada iteração corresponde à soma de três

termos distintos: o primeiro é um termo relativo à inércia da partícula e que traduz o

modo com que a partícula vem se movendo; o segundo é um termo relativo à atração da

partícula ao melhor ponto que já encontrou; e o terceiro termo é relativo à atração da

Page 30: Marcio Schwa Ab

23

partícula ao melhor ponto que todo o grupo (ou uma parte do grupo) já encontrou. Desta

forma, KENNEDY e EBERHART (1995) propuseram as seguintes equações:

( ) ( )1, , 1 , , 2 , ,k k k k k ki d i d 1 i d i d 2 global d i dv v c r p x c r p x+ = + ⋅ ⋅ − + ⋅ ⋅ − (2.22a)

1 1, , ,k k ki d i d i dx x v+ += + (2.22b)

onde os índices k, i e d denotam, respectivamente, a iteração, a partícula e a direção de

busca; v é a velocidade e x a posição no espaço de busca; c1 e c2 são duas constantes

positivas, chamadas respectivamente de parâmetro cognitivo e social; r1 e r2 são dois

números aleatórios com distribuição uniforme no intervalo [0, 1] e são sempre

diferentes para cada direção, partículas e iteração; pi é o melhor ponto encontrado pela

partícula i e pglobal é o melhor valor encontrado por todo enxame.

O termo Enxame foi utilizado em acordo com o trabalho de Milonas

(KENNEDY e EBERHART, 1995) que desenvolveu modelos para aplicações em vida

artificial e articulou cinco princípios básicos da Inteligência de Enxames, sobre os quais

o método do Enxame de Partículas é fundamentado. Os princípios são:

a) Proximidade: o enxame deve ser capaz de realizar cálculos simples de tempo e

espaço. O método do Enxame de Partículas realiza uma série de cálculos no

espaço n-dimensional em muitos intervalos de tempo.

b) Qualidade: o enxame deve ser capaz de responder a fatores de qualidade do

ambiente. O método responde a fatores de qualidade definidos pela melhor

posição encontrada pelo enxame e por cada melhor posição encontrada por cada

partícula.

c) Respostas Diversas: o enxame não deve submeter sua atividade em meios

excessivamente limitados. A definição do movimento de cada partícula do

enxame em relação às melhores posições individuais e à melhor posição do

enxame garante a diversidade das respostas.

d) Estabilidade: o enxame não deve mudar seu comportamento a todo momento que

o ambiente se altera. O comportamento do procedimento só é alterado quando

ocorre mudança nos melhores valores encontrados.

e) Adaptabilidade: o enxame deve ser capaz de alterar seu comportamento, quando o

custo computacional não for proibitivo. Assim que os melhores valores são

alterados, o comportamento do enxame se adapta aos novos valores.

Page 31: Marcio Schwa Ab

24

É interessante observar que os dois últimos princípios são opostos, mas o

método do Enxame assegura que ambos sejam satisfeitos, já que o melhor valor

encontrado não é, necessariamente, alterado a todo instante. Quando este é alterado, o

enxame se adapta a este novo valor, de forma a assegurar estabilidade e adaptabilidade

ao método.

A versão proposta inicialmente sofreu algumas pequenas alterações. A primeira

alteração foi proposta pelos próprios “inventores” do método e consistiu em considerar

que cada partícula, além da influência do melhor ponto encontrado por ela mesma, sofre

a influência do melhor ponto encontrado pelas partículas vizinhas a ela, e não pelo

melhor valor encontrado por todo o enxame. Esta versão “local” do algoritmo promove

uma maior exploração do espaço de busca, mas requer em média mais iterações para

que o enxame convirja (EBERHART e KENNEDY, 1995).

A segunda alteração, proposta por SHI e EBERHART (1998a), consistiu na

introdução de um novo parâmetro, chamado de peso de inércia (inertia weight) ou fator

de inércia, o qual pondera o termo relativo à velocidade prévia da partícula, de acordo

com as seguintes equações:

( ) ( )1, , 1 , , 2 , ,k k k k k ki d i d 1 i d i d 2 global d i dv w v c r p x c r p x+ = ⋅ + ⋅ ⋅ − + ⋅ ⋅ − (2.23a)

1 1, , ,k k ki d i d i dx x v+ += + (2.23b)

O papel fundamental do peso de inércia w é balancear o caráter global e local da

busca. Este peso pode ser uma constante positiva ou mesmo uma função do tempo

(iterações) positiva linear ou não linear. SHI e EBERHART (1998a; 1998b) observaram

que o número de iterações necessárias para atingir o valor mínimo aumenta quando se

aumenta o valor de w. Por outro lado, o número de falhas aumenta quando se diminui o

valor de w. Isto mostra que valores maiores de w aumentam o caráter exploratório e

global da busca, aumentando as chances de encontrar o valor ótimo, enquanto valores

pequenos de w aumentam o caráter local da busca, fazendo com que as partículas

convirjam rapidamente, para o ponto ótimo ou não, sem que haja uma boa exploração

do espaço de busca.

Desta forma, uma alternativa interessante é iniciar a busca com um valor alto de

w e diminuir este valor ao longo das iterações, como proposto por SHI e EBERHART

(1998a; 1998b, 1999), onde propuseram que o valor de w deve diminuir linearmente

entre 0,9 e 0,4, aproximadamente. Valores inferiores a 0,4 “congelam” o conjunto de

Page 32: Marcio Schwa Ab

25

partículas muito rapidamente, enquanto valores superiores a 0,9 podem causar a

“explosão” das partículas. Uma forma de evitar a explosão do conjunto de partículas é

limitar o valor da velocidade de cada partícula. Entretanto, definir este limite máximo

da velocidade, além de implicar na definição de um parâmetro adicional, depende do

problema que se está resolvendo. Assim, uma forma simples de controlar a velocidade

das partículas, evitando a “explosão”, é evitar que uma partícula ultrapasse os limites do

espaço de busca zerando a sua velocidade.

TRELEA (2003) avaliou a convergência das partículas através de resultados da

teoria de sistemas dinâmicos e mostrou que o conjunto de partículas converge sempre

que o valor absoluto de w é menor que 1, a despeito dos valores de c1 e c2, que também

influenciam a convergência das partículas. Assim, a princípio não é necessário diminuir

o valor de w ao longo das iterações, pois mesmo com um valor de w igual a 0,9, por

exemplo, a busca tem um caráter global. Ao longo das iterações, à medida que as

partículas convergem, o caráter da busca passa a ser local. Porém, neste caso, a

velocidade com que as partículas convergem depende do problema que está sendo

resolvido. Assim, pode ser mais fácil controlar o número de iterações forçando a

convergência através da diminuição do valor de w ao longo das iterações.

Uma terceira versão do método do Enxame de Partículas foi proposta por

CLERC (1999), através do uso de um fator de constrição K, o qual pode ser necessário

para assegurar a convergência do algoritmo. As equações desta nova versão são:

( ) ( )1, , 1 , , 2 , ,k k k k k ki d i d 1 i d i d 2 global d i dv K v c r p x c r p x+ = ⋅ + ⋅ ⋅ − + ⋅ ⋅ − (2.24a)

1 1, , ,k k ki d i d i dx x v+ += + (2.24b)

2

2

2 4K

ϕ ϕ ϕ=

− − − ⋅ (2.24c)

onde ϕ = c1 + c2 e ϕ > 4. Entretanto, conforme observado por EBERHART e SHI

(2000), a forma do algoritmo descrita pela Equações (2.24) é idêntica à descrita na

Equação (2.23), sendo a única diferença a definição dos valores dos parâmetros de

busca. Uma vantagem do algoritmo proposto por CLERC (1999) é a autocorreção dos

valores dos parâmetros de busca. Se os valores de c1 e c2 são altos, o que pode provocar

a “explosão” do enxame, estes valores são multiplicados pelo fator de constrição, que é

menor que 1, impedindo a “explosão” do enxame. CLERC e KENNEDY (2002)

Page 33: Marcio Schwa Ab

26

analisaram a estabilidade e convergência do método do Enxame de Partículas, chegando

a uma forma generalizada do algoritmo, contendo um conjunto de coeficientes que

controlam a tendência de convergência do método.

Como já foi citado acima, TRELEA (2003) avaliou a convergência das

partículas através de resultados da teoria de sistemas dinâmicos e definiu uma região

que delimita os valores dos parâmetros w e c (onde c é o valor médio de c1 e c2) onde a

convergência do algoritmo é garantida. Mostrou também os diferentes comportamentos

oscilatórios que podem ocorrer durante movimento das partículas do enxame, de acordo

com os valores dos parâmetros w e c.

Apesar de ser um método proposto recentemente, o Enxame de Partículas tem

apresentado bons resultados em diferentes aplicações. PARSOPOULOS e VRAHATIS

(2002a) avaliam a eficiência do Enxame de Partículas em diferentes problemas de

otimização, e propuseram uma forma muito interessante de resolver problemas de

otimização multiobjetivo com este método, onde enxames paralelos trocam informações

entre si, sendo cada enxame avaliado por uma função objetivo específica. Em outro

trabalho, PARSOPOULOS e VRAHATIS (2002b) fizeram a inicialização das

partículas, que é geralmente feita de forma aleatória, utilizando o método Simplex,

procurando detectar mais rapidamente regiões promissoras e acelerar a busca pelo

mínimo global.

OURIQUE et al. (2002) fizeram uma aplicação muito interessante do método do

Enxame de Partículas para a análise dinâmica não linear de processos químicos.

COSTA et al. (2003) aplicaram o Enxame de Partículas na otimização da polimerização

de estireno em reatores tubulares.

Como o método do Enxame de Partículas não necessita da diferenciação da

função objetivo e é capaz de resolver problemas de grande dimensão, PARSOPOULOS

et al. (2001) aplicaram o método a problemas de estimação de parâmetros onde todas as

variáveis são consideradas na estimação (variáveis independentes e dependentes). Os

resultados mostraram que este método é capaz de resolver estes problemas com grande

eficiência.

Recentemente, BISCAIA et al. (2004) propuseram uma nova forma do Enxame

de Partículas baseada na dinâmica de sistemas dinâmicos lineares de segunda ordem

subamortecidos. As características principais da forma tradicional do método do

Enxame de Partículas são mantidas, sendo a principal diferença relacionada à definição

dos parâmetros de busca, já que no novo algoritmo estes parâmetros estão relacionados

Page 34: Marcio Schwa Ab

27

ao amortecimento e à freqüência das oscilações, facilitando a escolha destes parâmetros

de busca.

A maior dificuldade associada ao uso dos métodos heurísticos de otimização está

relacionada à necessidade de um número muito grande de avaliações da função

objetivo. Entretanto, este grande número de avaliações da função objetivo pode ser

usado para a determinação da região de confiança dos parâmetros, a qual pode então ser

obtida sem a necessidade de aproximações. Desta forma, a principal desvantagem dos

algoritmos heurísticos é transformada em uma vantagem, já que permite uma análise

estatística mais rigorosa dos resultados, como será discutido no próximo item.

É importante ressaltar que comparações entre os métodos tradicionais de

otimização e os métodos heurísticos são sempre difíceis e até mesmo impróprias. O que

deve estar em questão é o tipo de problema que está sendo resolvido. Em problemas que

não envolvam muitas variáveis, onde não haja dificuldades relacionadas ao cômputo das

derivadas e a escolha de uma estimativa inicial, os métodos tradicionais devem ser

aplicados, pois são muito mais rápidos que os métodos heurísticos, além de garantirem

que chegam a um ponto ótimo (as derivadas de primeira ordem são nulas). Por outro

lado, em problemas de grande dimensão, onde não se dispõe de uma boa estimativa

inicial, onde existem mínimos locais ou existem problemas relacionados às derivadas,

como a impossibilidade de inverter a matriz Hessiana, a aplicação dos métodos

tradicionais é inviável e os métodos heurísticos constituem uma ferramenta eficiente e

robusta para solução destes problemas.

2.3 A Interpretação Estatística dos Resultados

É importante ressaltar aqui que a interpretação estatística dos resultados da

estimação de parâmetros se inicia durante a formulação de uma função objetivo

consistente com as condições experimentais. Mais especificamente, a estrutura da

matriz de variância dos desvios experimentais deve ser conhecida e utilizada no

procedimento de estimação de parâmetros. A interpretação estatística dos resultados

pode incluir diversos aspectos, tais como o valor da função objetivo, a análise dos

resíduos, a correlação entre os dados experimentais e os valores preditos pelo modelo,

entre outros. Entretanto, este trabalho se restringe à análise da região de confiança dos

parâmetros e ao uso dos algoritmos heurísticos de otimização para, além de minimizar a

função objetivo, determinar a região de confiança dos parâmetros.

Page 35: Marcio Schwa Ab

28

Admitindo-se que os desvios entre os dados experimentais e as predições do

modelo têm distribuição normal, a função de máxima verossimilhança pode ser

interpretada como uma soma do quadrado de variáveis com distribuição normal.

Portanto, a distribuição de probabilidades associada ao valor da função objetivo é a

distribuição chi-quadrado (χ2) com N-NP graus de liberdade, onde N = NE.NY. Então é

possível determinar, com um nível de confiança α o intervalo esperado para a função

objetivo e, desta forma, avaliar o ajuste do modelo aos dados experimentais.

Como os valores experimentais utilizados para a estimação dos parâmetros são

variáveis que seguem uma determinada distribuição de probabilidade, ou seja, não são

valores exatos, os parâmetros estimados também devem ser vistos desta forma.Logo, a

caracterização dos parâmetros estimados envolve tanto a definição do valor que

minimiza a função objetivo como uma medida da incerteza deste valor. Esta medida

caracteriza uma região em torno do ponto mínimo onde o valor da função objetivo não

aumenta significativamente com a mudança nos valores dos parâmetros. Esta região em

torno do ponto mínimo é chamada de região de confiança dos parâmetros.

A forma tradicional para a determinação da região de confiança consiste na

aproximação quadrática da função objetivo em torno do ponto mínimo, como na

equação abaixo:

( ) ( ) ( ) ( ) ( )1ˆ ˆ2

ˆ ˆ ˆ ˆTS S= + − ∇ + − −θ θθ θ θ θ S θ θ H θ θ (2.31)

O termo onde aparece o gradiente é descartado, já que no ponto mínimo o

gradiente é nulo. A matriz Hessiana é calculada pela Equação (2.15). Definindo a matriz

de covariância dos parâmetros como:

1

1

1

ˆ ˆTNEi i

eyi

θy yV Vθ θ

=

∂ ∂ = ∂ ∂ ∑ (2.32)

onde y corresponde à variável dependente calculada com os valores estimados dos

parâmetros. A Equação (2.31) pode ser reescrita como:

( ) ( ) ( ) ( )1ˆ ˆ ˆTS S −− = − −θθ θ θ θ V θ θ (2.33)

Para modelos lineares nos parâmetros a aproximação quadrática da função

objetivo é exata. Deste modo, se os desvios entre os dados experimentais e as predições

do modelo têm distribuição normal, os parâmetros também são normalmente

Page 36: Marcio Schwa Ab

29

distribuídos e cada um dos lados da Equação (2.33) segue a distribuição χ2 com NP

graus de liberdade. Assim, para modelos lineares nos parâmetros a região de confiança

tem a forma elíptica, definida pelo lado direito da Equação (2.33), com o limite superior

definido pela distribuição χ2 com NP graus de liberdade e um nível de confiança de α,

como na equação abaixo:

( ) ( ) ( )1 2ˆ ˆ ,T

NPχ α−− − ≤θθ θ V θ θ (2.34)

Entretanto, rigorosamente, a matriz de covariância dos parâmetros não é

conhecida com exatidão, já que é determinada por um certo número limitado de

experimentos que representam apenas uma amostra de todos os possíveis experimentos.

Como foi discutido acima, a função objetivo tem distribuição χ2 com N-NP graus de

liberdade e a razão entre o lado direito da Equação (2.33) e a função objetivo, cada qual

ponderado pelos seus respectivos graus de liberdade, tem distribuição F de Fisher com

NP e N-NP graus de liberdade, respectivamente, como na seguinte equação:

( ) ( )

( ) ( )( )

1ˆ ˆ, ,ˆ

TNP

F NP N NPS N NP

α− − − ≤ −

θθ θ V θ θ

θ (2.35)

que pode ser reescrita como:

( ) ( ) ( ) ( )1ˆ ˆ ˆ , ,T NPS F NP N NP

N NPα−− − ≤ −

−θθ θ V θ θ θ (2.36)

Para modelos não lineares nos parâmetros, a região obtida pela Equação (2.36)

ou (2.34) não é exata e a qualidade da aproximação da região de confiança é dependente

da qualidade da aproximação linear. A vantagem de utilizarem-se regiões de confiança

obtidas por aproximação linear é que o ponto mínimo é o centro da região; ou seja, a

região de confiança é simétrica, a distribuição dos parâmetros é normal e o ponto de

mínimo pode ser visto como o ponto médio. Além disso, a forma elíptica da região

facilita a interpretação geométrica. Adicionalmente, como toda região é definida pelo

ponto mínimo e pela matriz de covariância dos parâmetros, a utilização desta região

para procedimentos subseqüentes ao da estimação, como o planejamento de

experimentos, são matematicamente facilitados (BARD, 1974).

Quando a aproximação linear do modelo e, conseqüentemente, a aproximação

quadrática da função objetivo não apresentam boa qualidade, os resultados obtidos pela

Page 37: Marcio Schwa Ab

30

Equação (2.36) podem levar a conclusões equivocadas sobre os parâmetros e sobre a

região de confiança. Uma forma de melhorar os resultados obtidos na determinação da

região de confiança através da aproximação linear é fazer o uso de transformações nos

parâmetros, que podem diminuir a correlação entre os parâmetros e melhorar a

qualidade da região de confiança (BATES e WATTS, 1981).

Outra forma de se obter a região de confiança é, ao invés de se utilizar o lado

direito da Equação (2.33), utilizar o lado esquerdo desta equação. Assim, a região de

confiança pode ser definida pelos parâmetros que satisfazem a seguinte equação:

( ) ( ) ( )2ˆ ,S S NPχ α− ≤θ θ (2.37)

Como discutido anteriormente, rigorosamente a matriz de covariância dos

desvios entre os dados experimentais e as predições do modelo não é conhecida com

exatidão, sendo apenas baseada em um conjunto limitado de experimentos. A razão

entre o lado esquerdo da Equação (2.33) e a função objetivo, cada qual ponderado pelos

seus respectivos graus de liberdade, tem distribuição F de Fisher com NP e N-NP graus

de liberdade. Então, a seguinte equação pode ser obtida:

( ) ( )( )( ) ( )

( )ˆ

, ,ˆ

S S NPF NP N NP

S N NPα

−≤ −

θ θ

θ (2.38)

A razão presente na Equação (2.38) é conhecida como razão de verossimilhança

e, por isto, a região de confiança obtida por esta via é chamada de região de

verossimilhança (DONALDSON e SCHNABEL, 1987). A Equação (2.38) pode ser

reescrita como:

( ) ( ) ( )ˆ 1 , ,NPS S F NP N NPN NP

α ≤ + − −

θ θ (2.39)

Apesar da região de confiança obtida pelas Equações (2.37) ou (2.39) só ser

exata para modelos lineares nos parâmetros, o erro fica restrito ao fato de que a Equação

(2.37) não tem distribuição χ2, já que para modelos não-lineares nos parâmetros a

distribuição de probabilidade dos parâmetros não é normal mesmo que os desvios entre

os valores experimentais e os valores calculados pelo modelo tenham distribuição

normal. Assim, o erro associado à definição da região de confiança a partir da Equação

(2.39) está ligado ao nível de confiança, ou seja, o valor da distribuição F não condiz

com a realidade. Entretanto, a forma da região de confiança é muito próxima da real

Page 38: Marcio Schwa Ab

31

(DONALDSON e SCHNABEL, 1987), o que não acontece na região determinada pelas

Equações (2.34) ou (2.36), onde a forma da região de confiança é sempre elíptica..

Para a determinação exata da região de confiança, WILLIAMS (1962),

HALPERIN (1962) e HARTLEY (1964) desenvolveram um método onde a soma dos

quadrados dos desvios é dividida em duas partes independentes, sendo estas partes

variáveis aleatórias com distribuição χ2. Assim, a razão entre estas duas partes tem

distribuição F e a determinação da região de confiança consiste nos valores de θ que

satisfazem a seguinte expressão:

( ) ( ) ( )( ) ( ) ( )

( ), ,

, ,, ,

Te m m e m m

Te m m e m m

NP F NP N NPN NP

α − − ≤ −

− − − −

y y x θ P θ y y x θ

y y x θ I P θ y y x θ (2.40)

onde I é a matriz identidade e P(θ) e a matriz definida como:

( ) ( ) ( ) ( ) ( )1T T− = P θ J θ J θ J θ J θ (2.41)

sendo J(θ) a matriz Jacobiana (N x NP) definida da seguinte forma:

mi

ijj

∂=∂

J (2.42)

É interessante observar que este método (chamado de “Lack-of-Fit Method”) não

necessita “a priori” da solução do procedimento de estimação, ou seja, a determinação

da região de confiança é independente do ponto mínimo da função objetivo

(DONALDSON e SCHNABEL, 1987). RATKOWSKY (1990) adverte que, além das

dificuldades computacionais associadas à determinação da região de confiança exata, já

que é necessário o cálculo da matriz Jacobiana para cada conjunto de parâmetros, esta

pode apresentar características indesejáveis, como consistir de sub-regiões

desconectadas. Por outro lado, no desenvolvimento deste método a função de mínimos

quadrados foi considerada. Deve-se lembrar que esta função objetivo implica na

hipótese de erros com distribuição normal e variância constante. Nos casos onde estas

hipóteses não podem ser aplicadas, a determinação da região de confiança por este

método fica comprometida.

DONALDSON e SCHNABEL (1987) compararam os resultados obtidos pelos

três métodos descritos acima (aproximação quadrática, razão de verossimilhança e

“Lack-of-Fit Method”) para determinar a região de confiança e concluíram que, nos

Page 39: Marcio Schwa Ab

32

problemas estudados, não há diferenças significativas entre as regiões obtidas pelo

método exato e pelo de verossimilhança. Já as regiões obtidas pela linearização do

modelo podem levar a resultado inverossímeis.

É importante ressaltar que os critérios definidos acima para obtenção da região

de confiança são baseados ou em aproximações quadráticas da função objetivo

(Equações (2.34), (2.36), (2.37) e (2.39)) ou são desenvolvidos a partir de hipóteses que

nem sempre são apropriadas, como a Equação (2.40), que considera a função de

mínimos quadrados para estimação dos parâmetros. Uma forma simples e exata de se

determinar a região de confiança dos parâmetros é através da própria definição da

função objetivo como uma variável estatística. Assim, determina-se apenas um limite

superior para esta variável. De uma forma simples, a região de confiança poderia ser

definida pelos valores dos parâmetros que satisfazem à seguinte expressão:

( ) ( )ˆS cS≤θ θ (2.43)

onde c é uma constante (c > 1) definida em termos estatísticos de forma a delimitar o

limite superior de confiança para a função objetivo. Esta região poder ser considerada

exata, já que não é baseada em nenhuma aproximação; porém, é difícil definir o valor

de c com algum significado estatístico (MARSILI-LIBELLI et al., 2003). A Equação

(2.39) tem a forma da Equação (2.43) e, quando o valor de N é grande, a região de

confiança definida pela Equação (2.39) assume um nível de confiança assintótico; ou

seja, conforme N aumenta a região se aproxima da região verdadeira. É interessante

observar que, apesar do nível de confiança definido pela constante c não ser exato, a

forma da região de confiança é muito próxima da real, o que não ocorre quando são

definidas regiões de confiança elípticas (Equações (2.34) e (2.36)).

Apesar de ser possível obterem-se regiões de confiança dos parâmetros muito

próximas da real, a utilização de equações do tipo da Equação (2.43) é pouco difundida.

Isso ocorre por causa das dificuldades computacionais associadas à obtenção da região,

já que a solução destas equações não é simples. Alguns autores sugerem a utilização de

extensivas simulações de Monte Carlo para obtenção da região de confiança

(DONALDSON e SCHNABEL, 1987; VANROLLEGHEM e KEESMAN, 1996).

Uma outra forma de obtenção dos intervalos de confiança dos parâmetros de

modelos não lineares consiste na construção dos gráficos de perfil t (“profile t plots”),

como sugerido por BATES e WATTS (1988). Os gráficos obtidos fornecem

informações sobre a característica não linear dos parâmetros estimados. BATES e

Page 40: Marcio Schwa Ab

33

WATTS (1988) usaram a função de mínimos quadrados para definir a seguinte

estatística t para θi:

( ) ( )ˆsign( )i i i iFτ θ θ θ θ= − (2.44)

sendo F(θi) definido abaixo:

( )( ) ( )( ) ( )

ˆ,ˆi i

i

S SF

S N NP

θθ

− −=

Θ θ

θ (2.45)

onde ( ),i iS θ −Θ o valor mínimo da função objetivo com o parâmetro θi fixo e re-

estimando os parâmetros restantes, que fazem parte do vetor i−Θ . Comparando τ(θi)

com a distribuição tN-NP é possível obter o intervalo de confiança do parâmetro θi. Já o

gráfico de τ(θi) contra θi é chamado de gráfico de perfil t e pode ser usado para uma

análise da não linearidade do modelo, já que para modelos lineares este gráfico é uma

reta. SULIEMAN et al. (2001, 2004) aplicam este método para análise de sensitividade

de modelos de regressão não linear.

Este método apresenta bons resultados para determinação dos intervalos de

confiança dos parâmetros e para análise da não linearidade dos modelos. Porém, este

método não fornece de forma direta a região de confiança dos parâmetros. Por outro

lado, a necessidade de inúmeras minimizações da função objetivo, sempre fixando um

parâmetro e re-estimando os demais, inviabiliza a aplicação deste método em problemas

com muitos parâmetros.

KLEPPER e HENDRIX (1994a, 1994b) e KLEPPER e BEDAUX (1997a,

1997b) utilizaram um método aleatório de otimização e incorporaram a Equação (2.43)

no algoritmo de otimização, procurando cobrir uniformemente a região de confiança dos

parâmetros através de extensivas avaliações da função objetivo. A região obtida é

comparada com a região de confiança obtida pela aproximação quadrática da função

objetivo, mostrando que regiões diferentes são obtidas, fortalecendo a necessidade de

uma análise estatística mais rigorosa em problemas de estimação de parâmetros de

modelos não lineares.

MARSEGUERRA et al. (2003) usaram o Algoritmo Genético para estimação de

parâmetros de um modelo de reator nuclear e mostraram que os pontos avaliados pelo

algoritmo durante a solução podem ser usados para obter informações sobre a

convergência e estabilização dos diferentes parâmetros do problema de estimação.

Page 41: Marcio Schwa Ab

34

Assim, a utilização dos algoritmos heurísticos de otimização para a estimação de

parâmetros, além das vantagens relacionadas ao procedimento de minimização, tem a

vantagem de possibilitar a determinação da região de confiança sem representar custos

adicionais de computação, já que os próprios pontos avaliados durante a minimização

podem ser utilizados para descrever a região de confiança, de acordo com algum critério

semelhante à Equação (2.43). Desta forma, como já foi dito acima, a principal

desvantagem dos métodos heurísticos é convertida uma importante vantagem,

permitindo uma avaliação estatística rigorosa dos parâmetros estimados, sem a

necessidade de aproximações que podem levar a regiões que, além de não representar a

forma real da região, podem conduzir a conclusões equivocadas sobre os níveis de

confiança dos parâmetros.

Page 42: Marcio Schwa Ab

35

CAPÍTULO 3: OS ALGORITMOS HEURÍSTICOS

Neste capítulo são apresentados os detalhes de implementação dos algoritmos

heurísticos. Em seguida, os algoritmos são avaliados e comparados na minimização de

funções conhecidas e que apresentam características específicas, como a presença de

mínimos locais, funções com um grande número de variáveis e a alta correlação entre as

variáveis. Os métodos são avaliados quanto a sua eficiência, relacionada ao número de

avaliações da função objetivo, e quanto a sua robustez, relacionada ao número de

sucessos em encontrar o mínimo global.

São utilizadas duas funções para a avaliação dos métodos de otimização. A

primeira função é a função de Rosenbrock (ROSENBROCK, 1960; TRELEA, 2003).

Esta é uma função de n variáveis, conforme a Equação (3.1). Na Figura 1 pode ser

observada gráficos desta função para duas variáveis, onde fica visível a forma em curva

da região próxima ao mínimo e o rápido aumento da função na medida que se afasta do

mínimo, dificultando muito o procedimento de minimização desta função.

( ) ( ) ( )1 2 22

11

100 1n

i i ii

f x x x−

+=

= − + −∑x (3.1)

Figura 1: Representação da Função de Rosenbrock para duas variáveis.

Neste trabalho, a avaliação dos métodos será feita usando a função de

Rosenbrock com dimensão n igual a 30 e delimitando o espaço de busca no intervalo

dado por x = [-30, 30]n. Esta função apresenta um mínimo global em xi = 1, i = 1 ... n,

onde a função é igual a 0. Devido ao grande número de dimensões usadas e à

Mínimo Global

Page 43: Marcio Schwa Ab

36

dificuldade em minimizar esta função, sempre que a busca resultar em um valor inferior

a 100 a busca será admitida como um sucesso, como feito por TRELEA (2003).

A segunda função utilizada na avaliação dos métodos é a função de Levy Nº 5

(PARSOPOULOS e VRAHATIS, 2002a) ou, simplesmente, função de Levy. Esta

função de duas variáveis apresenta no espaço de busca delimitado por x = [-10, 10]2

cerca de 760 mínimos locais e um mínimo global localizado em x1 = -1,3068 e x2 =

-1,4248, onde a função assume o valor aproximado de -176,14. Desta forma, a busca foi

limitada ao intervalo [-10, 10] para cada variável. Sempre que a minimização encontrar

um valor menor que -176, definido de forma arbitrária, a busca será considerada como

sucesso.

Devido ao grande número de mínimos locais e de somente um mínimo global, a

minimização desta função avalia a capacidade dos métodos de encontrar o mínimo

global de uma função, ou seja, dentre todos os mínimos da função encontrar e convergir

para o mínimo global. A Equação (3.2) mostra a forma matemática desta função. A

Figura 2 mostra o grande número de mínimos e máximos locais que esta função

apresenta.

( ) ( )( ) ( )( )

( ) ( )

5 5

1 2 1 21 1

2 21 2

, cos 1 cos 1

1,42513 0,80032

i j

f x x i i x i j j x j

x x

= =

= ⋅ − ⋅ + ⋅ ⋅ + ⋅ + +

+ + + +

∑ ∑ (3.2)

Figura 2: Função de Levy em 3 dimensões e em gráfico de contorno.

O desempenho de cada algoritmo, em cada uma das duas funções apresentadas

acima, foi avaliada através do número de sucessos que cada um obteve em 50 corridas,

pela média aritmética (Iter) e pelo desvio padrão (Dpiter) das iterações das corridas que

Mínimo Global

Page 44: Marcio Schwa Ab

37

obtiveram sucesso e pelo valor esperado de avaliações da função objetivo, sendo este

último valor calculado pela seguinte equação (TRELEA, 2003):

NF Npt Iter PS= ⋅ (3.3)

onde NF é o número esperado de avaliações da função objetivo, Npt é o número de

pontos avaliados por iteração, Iter é a média aritmética das iterações das corridas que

resultaram em sucesso e PS é o percentual de sucesso. Assim, o número esperado de

avaliações da função objetivo leva em conta as corridas que não conseguem atingir o

mínimo, aumentando o número esperado de avaliações da função objetivo na medida

que a porcentagem de sucessos diminui.

Os métodos heurísticos de otimização também serão avaliados quanto à

qualidade da região de confiança obtida para estimação de parâmetros de um modelo

linear, de acordo com a seguinte equação:

1 2y xθ θ= ⋅ + (3.4)

onde y é a variável dependente, x é a variável independente e θ1 e θ2 são os parâmetros,

cujos valores são ajustados para a minimização da função objetivo. Neste caso, a

variável independente é considerada exata e os erros de medida da variável y têm

distribuição de probabilidade normal com variância constante. Portanto, a função

objetivo usada para a estimação dos parâmetros é a função de mínimos quadrados

(Equação (2.9)). Os dados foram obtidos através da Equação (3.4), com os valores de θ1

e θ2 ambos iguais a 5,0, perturbando-se o valor de y obtido com um erro aleatório com

distribuição normal, com média 0,0 e desvio padrão 2,0. Na Tabela 2 são reportados os

dados experimentais simulados usados para estimação dos parâmetros.

Tabela 2: Dados experimentais para estimação dos parâmetros do modelo linear.

x y 1,0 9,92 2,0 16,89 3,0 17,12 4,0 26,03 5,0 30,71 6,0 33,28 7,0 39,83 8,0 42,44 9,0 50,44 10,0 53,63

Page 45: Marcio Schwa Ab

38

Como o modelo cujos parâmetros estão sendo estimados é um modelo linear, a

região de confiança dos parâmetros é uma elipse descrita pela Equação (2.36). Assim, é

possível avaliar a qualidade da região obtida por cada método de otimização, já que a

região de confiança exata é conhecida. A região de confiança obtida pelos métodos de

otimização é obtida pelos pontos avaliados pelo método que satisfazem a Equação

(2.39).

A região elíptica é obtida pela minimização da função objetivo através de um

algoritmo chamado Estima desenvolvido por NORONHA et al.(1993), o qual utiliza um

método do tipo Gauss-Newton. No caso de modelo linear, obtém a solução exata em

apenas uma iteração e fornece a matriz de covariância dos parâmetros (Equação (2.32)),

que, no caso da função objetivo ser de mínimos quadrados, utiliza a matriz identidade

no lugar da matriz de covariância dos erros experimentais (Vey).

A estimação dos parâmetros com este método obtém os valores estimados dos

parâmetros θ1 e θ2 iguais a 4,84 e 5,40, respectivamente, valores estes onde a função

objetivo tem o valor mínimo igual a 20,78. A matriz de covariância dos parâmetros

obtida é:

2 1

165 1571

15 15

− =

θV (3.5)

Desta forma, substituindo os valores obtidos na estimação na Equação (2.36) e

considerando um nível de confiança de 95%, a região de confiança dos parâmetros do

problema definido acima é descrita pela seguinte equação:

1 1

2 2

4,84 385 55 4,8423,16

5,40 55 10 5,40

Tθ θθ θ

− ⋅ ⋅ − ≤ (3.6)

onde a matriz de covariância dos parâmetros já está invertida. Após a multiplicação

matricial fornece a seguinte equação algébrica:

2 21 2 1 2 1 2385,0 +10,0 +110,0 -4322,0 -640,6 +12169.2 0θ θ θ θ θ θ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ≤ (3.7)

A Equação (3.7) define uma elipse e pode ser vista na Figura 3, onde também

aparece o ponto mínimo, no centro da elipse.

Page 46: Marcio Schwa Ab

39

Figura 3: Região de confiança e ponto ótimo dos parâmetros estimados do modelo linear.

A região de confiança definida pela Equação (2.39) consiste de todos os

conjuntos de parâmetros em que a função objetivo é menor que o valor de 43,95. Assim,

após a minimização realizada pelos métodos heurísticos, todos os pontos avaliados

durante a busca são comparados com o valor máximo de 43,95. Os pontos onde o valor

da função objetivo é menor que este valor máximo são salvos e usados para a

construção da região de confiança.

A seguir o algoritmo de cada um dos métodos heurísticos de otimização será

apresentado e os métodos serão avaliados na minimização das duas funções descritas

acima e no problema de estimação dos parâmetros de um modelo linear, para avaliação

da qualidade da região de confiança obtida por cada método.

3.1. Monte Carlo

O método de Monte Carlo é um algoritmo heurístico que realiza uma busca

totalmente aleatória. Um número significativo de pontos é sorteado aleatoriamente na

região de busca e o melhor valor encontrado é considerado o valor mínimo. Para que um

resultado razoável seja encontrado, é necessário um número muito grande de avaliações

Page 47: Marcio Schwa Ab

40

da função objetivo, tornando este método computacionalmente intensivo, o que também

diminui a eficiência deste método. Apesar disso, o método de Monte Carlo é robusto, já

que não depende de nenhum parâmetro de busca e que, a despeito da necessidade do

grande número de avaliações da função objetivo, sempre funciona. Os pontos são

sorteados aleatoriamente na região de busca segundo a Equação (3.8), a qual é usada

para sortear cada ponto em todas as direções de busca e em cada iteração:

( )L H L= + ⊗ −x x r x x (3.8)

onde x é um vetor que representa o ponto sorteado, xL e xH são, respectivamente,

vetores com o limite inferior e superior da região de busca, r é um vetor contendo

números aleatórios com distribuição uniforme entre 0 e 1, de forma que a probabilidade

de um ponto ser escolhido é a mesma em toda a região de busca. O símbolo ⊗ indica

multiplicação de elemento por elemento dos vetores. Todos estes vetores têm dimensão

n, que corresponde ao número de variáveis que serão otimizadas; isto é, a dimensão do

problema.

Uma forma de aumentar a eficiência deste método é utilizar um procedimento

iterativo onde em cada iteração um número de pontos é avaliado (Npt) e a região de

busca é diminuída em torno do melhor ponto encontrado até então. Para tanto, é

necessário fornecer a taxa de redução da região de busca. Dados os limites iniciais da

região de busca (xL0: limite inferior inicial; xH0: limite superior inicial), os limites em

cada iteração podem ser calculados através das seguintes equações:

( ) ( )0 1 kL L otm otmTR= − ⋅ − +x x x x (3.9a)

( ) ( )0 1 kH H otm otmTR= − ⋅ − +x x x x (3.9b)

onde xotm é o vetor de dimensão n contendo o melhor ponto encontrado até a k-ésima

iteração e TR é a taxa de redução.

Quando a redução da região de busca é aplicada durante a busca, é possível

correlacionar o número máximo de iterações (Niter) com a taxa de redução, já que a

busca pode ser interrompida assim que a região de busca seja reduzida a uma pequena

faixa. A região de busca pode ser calculada em função do número de iterações através

da diferença entre as Equações (3.9a) e (3.9b), resultando na seguinte equação:

( )0 1 NiterTR∆ =∆ ⋅ −x x (3.10)

Page 48: Marcio Schwa Ab

41

onde ∆x e ∆x0 são, respectivamente, vetores contendo o intervalo de busca reduzido e

inicial em cada uma das n dimensões. Assim, fazendo a razão entre o intervalo de busca

reduzido e o inicial, tem-se uma idéia da precisão do valor encontrado, já que o menor

valor encontrado está dentro do intervalo de busca reduzido. Assim, dado um certo

número máximo de iterações (Niter) e uma direção de busca d qualquer, é possível

definir um valor apropriado da taxa de redução de acordo com a seguinte equação:

( )101 Niterd dTR x x= − ∆ ∆ (3.11)

observando que a razão entre a faixa de busca atual e a inicial é a mesma para todas as

direções.

Porém, é preciso tomar um certo cuidado com o valor de TR, pois um valor

muito alto pode diminuir a região de busca prematuramente, prejudicando a exploração

de toda a região de busca e possibilitando que um mínimo global seja perdido. Assim,

ao invés de calcular o valor de TR, pode ser mais interessante fornecer um valor de TR e

realizar a busca até que um critério de parada seja satisfeito, o qual pode ser dado pela

razão entre intervalo de busca reduzido e o inicial, um certo valor da função objetivo, ou

ainda, o número de iterações, entre outros.

Em cada iteração, o melhor ponto sorteado é selecionado e comparado com o

ponto ótimo até a iteração anterior. Caso um ponto melhor seja encontrado, o ponto

ótimo é atualizado e o procedimento iterativo continua.

A eficiência do método de Monte Carlo é influenciada de forma direta pelo

número de pontos avaliados da função objetivo e pela taxa de redução da região de

busca em cada iteração. Assim, o efeito da taxa de redução e o número de pontos

avaliados por iteração foram avaliados na minimização das funções de Rosenbrock e de

Levy (Equações (3.1) e (3.2), respectivamente), observando o número médio e o desvio

padrão das iterações das corridas que resultaram em sucesso e o número esperado de

avaliações da função objetivo, calculado pela Equação (3.3).

Foram usados três critérios de parada: o primeiro é aplicado no momento em que

a função objetivo atinge o limite proposto para cada função (100 para a função de

Rosenbrock e -176 para a função de Levy); o segundo consiste em limitar a busca a no

máximo 50000 iterações; o terceiro critério é aplicado com relação à redução da região

de busca sempre que o valor da razão ∆xd/∆x0d for inferior a 1.10-5 a busca é

interrompida. O primeiro critério de parada é relacionado a um sucesso da busca e os

outros dois são relacionados à insucessos. Geralmente, quando a taxa de redução é nula

Page 49: Marcio Schwa Ab

42

ou muito pequena os insucessos estão ligados ao número máximo de iterações e quando

a taxa de redução é alta os insucessos estão ligados ao valor da razão ∆xd/∆x0d inferior a

1.10-5. Na Figura 4 é apresentado um fluxograma básico deste método.

Figura 4: Fluxograma do método de Monte Carlo.

Page 50: Marcio Schwa Ab

43

Os resultados obtidos nas 50 corridas realizadas na minimização da função de

Rosenbrock são apresentados a seguir na Tabela 3.

Tabela 3: Resultados obtidos com o método de Monte Carlo na minimização da função de Rosenbrock.

Npt TR [%] Sucessos Iter DPIter NF 0,00 0 ----- ----- ----- 0,01 50 44617 676 446170 0,05 50 9260 184 92602 0,10 50 4777 274 47766 0,50 50 1019 66 10186 1,00 50 522 50 5221 2,00 49 278 29 2842 5,00 34 130 18 1916

10

10,00 0 ----- ----- ----- 0,00 0 ----- ----- ----- 0,01 50 43949 581 1318460 0,05 50 9031 192 270917 0,10 50 4572 85 137145 0,50 50 985 141 29537 1,00 50 506 55 15179 2,00 50 261 26 7834 5,00 48 109 10 3405

30

10,00 19 67 12 5265 0,00 0 ----- ----- ----- 0,01 50 43668 619 2183380 0,05 50 9023 432 451168 0,10 50 4532 72 226615 0,50 50 938 32 46924 1,00 50 481 31 24042 2,00 50 248 17 12418 5,00 48 106 12 5512

50

10,00 36 63 11 4365

Os resultados obtidos na minimização da função de Rosenbrock mostram que a

redução da região de busca é fundamental para que o mínimo da função seja atingido, já

que sem a redução (TR igual a 0), nenhuma das 50 corridas conseguiu obter sucesso,

sendo que o número máximo de iterações foi alcançado e a busca interrompida. Nas

corridas em que a redução da região de busca foi realizada, verifica-se uma diminuição

no número médio de iterações e do número esperado de avaliações da função objetivo

com o aumento da taxa de redução, já que a região de busca é reduzida mais

rapidamente. Entretanto, esta queda nas iterações e no número de avaliações da função

objetivo também é acompanhada da queda do número de sucessos; isto é, do número de

corridas que alcançam valores da função objetivo inferiores a 100. Para a minimização

Page 51: Marcio Schwa Ab

44

desta função uma taxa de redução da ordem de 1% parece ser um valor adequado, já que

se consegue obter 100% de sucessos sem ser necessário um número muito grande de

iterações.

Um resultado interessante foi que o número médio de iterações não foi alterado

significativamente, à medida que o número de avaliações por iteração (Npt) aumentava.

Conseqüentemente, com um número menor de pontos avaliados por iteração, o valor

esperado de avaliações da função objetivo é menor para valores menores de Npt, já que

mesmo com um valor baixo de Npt o número de sucessos é de 100%, exceto para as

buscas realizadas com uma taxa de redução muito alta (5 e 10%) ou sem a taxa de

redução (0%).

Em seguida foi avaliada a eficiência do método de Monte Carlo para a

minimização da função de Rosenbrock com os mesmos valores de TR e Npt usados

acima, mas definindo a região de busca no intervalo [-5, 30]n; ou seja, fazendo com que

o mínimo da função fique mais próximo dos limites da região de busca. Nenhum

sucesso foi obtido em qualquer uma das configurações, mostrando de forma clara que a

redução de busca pode levar a resultados muito ruins quando o mínimo da função

objetivo fica próximo dos limites da região de busca, já que este valor pode ser excluído

da região de busca com a redução.

A próxima etapa consistiu na minimização da função de Levy, onde também

foram realizadas 50 corridas, com os critérios de parada que já foram definidos

anteriormente e cujos resultados são apresentados na Tabela 4.

Na minimização da função de Levy, o método de Monte Carlo conseguiu

encontrar o mínimo global mesmo sem a redução da região de busca, provavelmente

devido ao pequeno número de dimensões desta função quando comparada à função de

Rosenbrock. Nas corridas em que a redução da região de busca foi utilizada, os

resultados obtidos mostram que o aumento da taxa de redução provoca uma queda do

número de sucessos; isto é, o número de corridas que alcançam valores da função

objetivo inferiores a -176,0 diminui. Entretanto, o número médio de iterações

necessárias para alcançar este objetivo e o número esperado de iterações necessárias

para a minimização (NF) também diminuem.

Conforme foi observado na minimização da função de Rosenbrock, não é

observada uma queda significativa do número médio de iterações com o aumento do

número de avaliações por iteração, para um mesmo valor da taxa de redução.

Entretanto, na minimização da função de Levy, um valor adequado da taxa de redução

Page 52: Marcio Schwa Ab

45

deve ser de 0,1%, já que valores acima deste apresentam uma queda no número de

sucessos. Na minimização da função de Rosenbrock um valor de 1% foi considerado

adequado. Esta diferença deve ter ocorrido devido à grande quantidade de mínimos

locais que a função de Levy apresenta e uma taxa de redução elevada pode levar a

região a convergir em torno de um mínimo local.

Tabela 4: Resultados obtidos com o método de Monte Carlo na minimização da função de Levy.

Npt TR [%] Sucessos Iter DPIter NF 0,00 18 26149 15210 726367 0,01 50 13751 5513 137513 0,05 50 4295 1129 42948 0,10 50 2344 675 23439 0,50 43 662 143 7699 1,00 41 356 56 4337 2,00 38 199 33 2615 5,00 23 86 9 1873

10

10,00 29 47 4 808 0,00 35 22483 14166 963540 0,01 50 10585 3636 317548 0,05 50 3207 1062 96220 0,10 50 1999 618 59955 0,50 50 500 155 14990 1,00 47 304 59 9708 2,00 43 168 22 5868 5,00 38 76 16 3017

30

10,00 29 40 6 2067 0,00 43 14341 13011 833798 0,01 50 7845 3848 392255 0,05 50 3078 961 153886 0,10 50 1745 597 87249 0,50 50 467 150 23334 1,00 50 276 56 13776 2,00 49 160 31 8153 5,00 36 70 11 4878

50

10,00 35 37 6 2676

De uma forma geral, os resultados indicam que uma busca com poucas

avaliações por iteração e uma pequena taxa de redução é a melhor opção para este

método, já que dessa forma o número de avaliações da função objetivo é menor e a

probabilidade de sucesso é maior. Entretanto, é preciso tomar cuidado com a redução da

região de busca, principalmente se a função apresenta o mínimo global próximo dos

limites da região de busca ou muitos mínimos locais, já que este ponto pode ser

Page 53: Marcio Schwa Ab

46

excluído da região de busca, levando a busca para um mínimo local ou outro ponto

qualquer.

O método de Monte Carlo também foi avaliado com relação à região de

confiança obtida na solução do problema de estimação de parâmetros de um modelo

linear proposto acima. Foram usadas três configurações de busca, onde os parâmetros de

busca foram alterados para avaliar seus efeitos na região de confiança obtida. Como a

minimização desta função objetivo é muito mais fácil que a minimização das funções de

Resenbrock de Levy, foram usadas poucas iterações e avaliações por iteração, de modo

a dificultar a obtenção do mínimo e da região de confiança e, assim, permitir uma

melhor discriminação entre os resultados obtidos.

Na Tabela 5 as configurações e os resultados obtidos são relatados de forma

sucinta.

Tabela 5: Condições de busca e resultados obtidos pelo método de Monte Carlo na estimação de parâmetros de um modelo linear.

Condição MC01 MC02 MC03 Niter 100 100 200 Npt 10 10 05

TR [%] 5,0 2,0 1,0 S(θ) 20,78 20,79 20,79 θ1 4,84 4,84 4,83 θ2 5,40 5,39 5,45

PRC 596 137 136

Nas três condições o valor da função objetivo alcançado e dos parâmetros

estimados é muito próximo do valor encontrado pelo ESTIMA. Entretanto, o número de

pontos na região de confiança (PRC) obtido na condição MC01 é muito maior que nas

outras duas condições, apesar do número de avaliações da função objetivo ser o mesmo.

Como pode ser observado na Figura 5, onde a região obtida pelo método de Monte

Carlo com a configuração MC01 é comparada à região elíptica (Figura 3), o valor de

5% da taxa de redução concentra os pontos em uma região muito próxima ao mínimo,

prejudicando a qualidade da região de confiança obtida.

Diminuindo a taxa de redução para 2% (configuração MC02), apesar do número

de pontos na região de confiança diminuir, já que a região onde os pontos são sorteados

tem uma redução menor, a qualidade desta é melhor, já que estes pontos se encontram

mais distribuídos, como mostra a Figura 6.

Page 54: Marcio Schwa Ab

47

Com o objetivo de tentar aumentar a qualidade da região de confiança obtida,

sem, no entanto, aumentar o número de avaliações da função objetivo, foi usada uma

terceira configuração de busca (MC03), onde o número de iterações foi dobrado e o

número de avaliações por iteração diminuído à metade. Como o número de iterações é

maior, a taxa de redução foi diminuída para 1%, para evitar uma alta concentração de

pontos próxima do mínimo. A região obtida é apresentada na Figura 7. Como pode ser

visto, não há melhora significativa com relação à região apresentada na Figura 6.

Para aumentar a qualidade da região de confiança obtida seriam necessárias mais

avaliações da função objetivo, de forma a aumentar a cobertura da região de confiança

pelos pontos. Entretanto, o objetivo desta avaliação é observar os efeitos que os

parâmetros de busca têm na região de confiança obtida pelo método de otimização. Foi

possível observar que pequenas taxas de redução evitam uma concentração muito alta

de pontos nas proximidades do ponto mínimo.

Figura 5: Comparação da região de confiança obtida pelo método de Monte Carlo (MC01) com a região elíptica.

Page 55: Marcio Schwa Ab

48

Figura 6: Comparação da região de confiança obtida pelo método de Monte Carlo (MC02) com a região elíptica.

Figura 7: Comparação da região de confiança obtida pelo método de Monte Carlo (MC03) com a região elíptica.

Page 56: Marcio Schwa Ab

49

Assim sendo, tanto para avaliação da minimização das funções teste, como para

obtenção da região de confiança, as condições mais favoráveis consistem na utilização

de pequenas taxas de redução. Com relação ao número de pontos avaliados por iteração,

para a minimização a melhor condição consiste em poucas avaliações por iteração, de

modo a diminuir o número de avaliações da função objetivo para o mínimo ser

alcançado. Entretanto, para obtenção da região de confiança, quanto maior o número de

avaliações da função objetivo, melhor será a região de confiança obtida. Assim, quando

se deseja obter tanto o ponto mínimo quanto a região de confiança, a melhor

configuração parece ser a combinação de pequenas taxas de redução e o maior número

de avaliações por iteração, desde que a busca não seja prejudicada, pois um número

muito grande de avaliações por iteração impede que a busca progrida, tornando o

processo de minimização muito lento.

Enfim, o método de Monte Carlo conseguiu obter bons resultados tanto na

minimização quanto na obtenção da região de confiança dos parâmetros do modelo

linear, apesar de que cuidados devem ser tomados com relação à redução da região de

busca.

3.2. Algoritmo genético

O Algoritmo Genético é um método de otimização baseado na evolução dos

seres vivos. Durante o procedimento de otimização são realizadas operações que

procuram simular as condições de evolução dos seres vivos, como a seleção natural, o

cruzamento entre os indivíduos e a mutação. Estas operações são realizadas sobre um

conjunto de pontos distribuídos no espaço de busca, geralmente chamados de

indivíduos, procurando fazer com que estes pontos encontrem valor ótimo da função

objetivo.

Assim, cada indivíduo contém um conjunto de informações que descrevem sua

localização no espaço de busca; ou seja, cada indivíduo consiste em um vetor contendo

os valores das variáveis de busca que correspondem a uma determinada posição. Desta

forma, a codificação dos indivíduos é o próprio vetor das variáveis de busca de cada um

dos Npt pontos usados na busca.

A primeira geração de indivíduos (ou seja, o conjunto inicial de pontos) é gerada

de forma aleatória dentro da região de busca, de acordo com a Equação (3.8), dados os

limites da região de busca. Este conjunto inicial de indivíduos é avaliado segundo uma

função objetivo, cujo mínimo é procurado. O valor da função objetivo corresponde a

Page 57: Marcio Schwa Ab

50

uma medida da aptidão de cada indivíduo. No caso de uma minimização, quanto menor

o valor da função objetivo, maior é a aptidão do indivíduo.

Após a avaliação da aptidão de todos os indivíduos, é realizado o cruzamento

entre os indivíduos mais aptos. Para a seleção dos pares de indivíduos é utilizada a

seleção por torneio, onde dois indivíduos são sorteados aleatoriamente e aquele que tem

a maior aptidão (neste caso, o menor valor da função objetivo) é escolhido. Este

procedimento é repetido para a formação a escolha do segundo indivíduo do par. Os

dois indivíduos do par são combinados para a geração de um novo indivíduo.

Entretanto, a combinação dos indivíduos do par só ocorre de fato se a probabilidade do

cruzamento (Pcruz) for maior que um número aleatório com distribuição uniforme no

intervalo [0, 1], sorteado para cada par que irá ou não sofrer o cruzamento. Caso o

cruzamento não ocorra, um dos indivíduos selecionados, definido neste trabalho como o

primeiro indivíduo selecionado, é incorporado na nova geração sem sofrer nenhuma

alteração.

Caso o cruzamento ocorra, a combinação consiste em sortear um ponto de

quebra para os dois indivíduos. O novo indivíduo gerado é composto do segmento

inicial de um dos indivíduos selecionados e pelo segmento final do outro indivíduo

selecionado (BALLAND et al., 2000), de acordo com a Figura 8, onde o procedimento

de cruzamento entre dois indivíduos é ilustrado:

Figura 8: Operação de cruzamento entre dois indivíduos para

geração de um indivíduo novo.

Este procedimento de seleção e cruzamento é repetido até que seja formada uma

nova geração de indivíduos.

Page 58: Marcio Schwa Ab

51

Em seguida é realizada a mutação, que consiste na seleção de indivíduos de

acordo com uma probabilidade de mutação (Pmut) e modificá-los de forma aleatória,

garantindo uma maior exploração da região de busca. A forma de mutação utilizada

neste trabalho consiste em perturbar cada uma das direções de busca que definem o

indivíduo em torno dos valores atuais, de acordo com a seguinte equação:

( )2 1N = + ⋅ − ⊗∆x x r x (3.12)

onde xN e x são os vetores correspondentes ao novo indivíduo e ao indivíduo atual,

respectivamente, r é um vetor de números aleatórios com distribuição uniforme no

intervalo [0, 1] e ∆x é um vetor com o tamanho máximo da perturbação. Assim, a

perturbação em x é um valor entre -∆x e ∆x, de acordo com os valores aleatórios

contidos no vetor r, os quais são sorteados para cada mutação realizada. Como na

Equação (3.8), o símbolo ⊗ indica multiplicação elemento por elemento dos vetores. Na

Figura 9 é apresentada uma ilustração da operação de mutação para uma busca

bidimensional, onde aparece no centro da figura o ponto e a região em torno deste ponto

onde o novo ponto será gerado pela mutação do ponto atual.

Figura 9: Operação de mutação de um indivíduo,

mostrando a região onde o novo ponto será gerado.

Uma outra prática comum é a realização do elitismo, que consiste em manter o

melhor indivíduo da geração atual na geração seguinte, evitando que a melhor solução

encontrada seja perdida devido ao caráter aleatório do algoritmo.

Um fluxograma básico do algoritmo é apresentado na Figura 10.

Page 59: Marcio Schwa Ab

52

Figura 10: Fluxograma do método do Algoritmo Genético.

O Algoritmo Genético foi avaliado na minimização das funções objetivo

definidas anteriormente, observando as influências do número de pontos avaliados por

iteração, da probabilidade de mutação e do tamanho da perturbação realizada durante a

mutação, sendo que a probabilidade de cruzamento foi fixada em 70% para todas as

minimizações, já que este deve ser um valor alto (GOLDBERG, 1989), além de evitar a

Page 60: Marcio Schwa Ab

53

necessidade de se realizar um número muito grande de simulações para se observar a

influência deste parâmetro de busca.

Os critérios de parada consistem no valor limite para o valor da função objetivo

(100 para a função de Rosenbrock e -176 para a função de Levy) e o número máximo de

50000 iterações.

Nas Tabelas 6, 7, 8 e 9 são apresentados os resultados da minimização da função

de Rosenbrock, onde a perturbação realizada nas mutações corresponde a valores de até

1, 2, 5 e 10% da faixa total de cada direção de busca.

Os resultados obtidos e apresentados nas Tabelas 6, 7, 8 e 9 mostram que as

minimizações realizadas com uma probabilidade de mutação da ordem de 30 a 40%

obtêm os melhores resultados, alcançando um número de sucessos próximos a 50 (ou

seja, próximo a 100%) e com um baixo número médio de iterações, fazendo com que o

número esperado de avaliações da função objetivo seja menor que as minimizações com

probabilidade de mutação de 1, 10, 20 e 50%. Este comportamento foi observado para

todos os níveis de perturbações realizadas pela mutação.

Tabela 6: Resultados obtidos com o método do Algoritmo Genético na minimização da função de Rosenbrock, com perturbação de 1% nas mutações.

Npt Pmut [%] Sucessos Iter DPIter NF 1 0 ----- ----- ----- 10 30 23863 14927 397716 20 40 9205 9090 115061 30 42 5404 3994 64336 40 47 4735 4540 50371

10

50 41 6975 6395 85063 1 0 ----- ----- ----- 10 36 13103 12688 545956 20 44 5270 7719 179646 30 48 4488 7194 140264 40 48 2993 3841 93519

30

50 50 4665 5778 139963 1 3 22068 20529 18390278 10 41 7678 8485 468169 20 42 6064 12168 360923 30 48 3299 5554 171828 40 49 1953 3420 99633

50

50 45 7331 10091 407298

Page 61: Marcio Schwa Ab

54

Tabela 7: Resultados obtidos com o método do Algoritmo Genético na minimização da função de Rosenbrock, com perturbação de 2% nas mutações.

Npt Pmut [%] Sucessos Iter DPIter NF 1 0 ----- ----- ----- 10 8 39896 4149 2493469 20 35 21262 10323 303736 30 39 16424 10486 210564 40 40 15478 9603 193471

10

50 40 19435 9365 242936 1 0 ----- ----- ----- 10 26 29680 10585 1712328 20 45 9264 10407 308799 30 43 5126 7649 178828 40 45 4461 5304 148689

30

50 45 14459 8543 481967 1 0 ----- ----- ----- 10 34 22891 11325 1683155 20 46 6072 7111 330020 30 50 4109 7828 205440 40 47 4499 5962 239299

50

50 47 15214 9566 809238

Tabela 8: Resultados obtidos com o método do Algoritmo Genético na minimização da função de Rosenbrock, com perturbação de 5% nas mutações.

Npt Pmut [%] Sucessos Iter DPIter NF 1 0 ----- ----- ----- 10 0 ----- ----- ----- 20 0 ----- ----- ----- 30 16 37600 6650 1175006 40 8 37704 9020 2356492

10

50 0 ----- ----- ----- 1 0 ----- ----- ----- 10 0 ----- ----- ----- 20 24 31553 11321 1972060 30 41 16796 10623 614490 40 39 21665 11732 833253

30

50 15 39900 5647 3989993 1 0 ----- ----- ----- 10 0 ----- ----- ----- 20 41 25795 11292 1572893 30 48 14277 11431 743615 40 43 18234 10226 1060105

50

50 7 35929 7839 12831888

Page 62: Marcio Schwa Ab

55

Tabela 9: Resultados obtidos com o método do Algoritmo Genético na minimização da função de Rosenbrock, com perturbação de 10% nas mutações.

Npt Pmut [%] Sucessos Iter DPIter NF 1 0 ----- ----- ----- 10 0 ----- ----- ----- 20 0 ----- ----- ----- 30 0 ----- ----- ----- 40 0 ----- ----- -----

10

50 0 ----- ----- ----- 1 0 ----- ----- ----- 10 0 ----- ----- ----- 20 1 43547 ----- 65320500 30 23 31841 10267 2076578 40 19 36054 7398 2846356

30

50 0 ----- ----- ----- 1 0 ----- ----- ----- 10 0 ----- ----- ----- 20 7 30183 6767 10779745 30 30 26753 12017 2229417 40 19 29514 10313 3883470

50

50 0 ----- ----- -----

Os resultados também mostram que o aumento da perturbação da mutação

provoca uma queda do número de sucessos, assim como aumenta o número de iterações

necessárias para um valor da função objetivo inferior a 100 ser alcançado. Além disso,

quando a perturbação é de 5 e 10%, o número de sucessos é zero para grande parte das

configurações utilizadas, como mostram as Tabelas 8 e 9. Dessa forma, valores

pequenos de perturbação devem ser utilizados, já que estes valores levam a resultados

melhores.

Com relação ao número de pontos avaliados em cada iteração, ou melhor, ao

número de indivíduos em cada geração, vemos que o seu aumento provoca um aumento

do número de sucessos e uma queda do número de iterações. Porém, as minimizações

com uma perturbação de 1% na mutação apresentam um aumento do número esperado

de avaliações da função objetivo com o aumento do número de indivíduos, mesmo

tendo um aumento no número de sucessos e uma queda no número de iterações. Este

comportamento também foi observado nas minimizações realizadas pelo método de

Monte Carlo. Como os melhores resultados foram obtidos com perturbações pequenas,

realizar buscas com um pequeno número de pontos avaliados por iteração aparenta ser a

melhor alternativa, apesar do menor número de sucessos. Mas se o custo computacional

associado à avaliação da função objetivo não prejudicar a busca pelo mínimo, um

Page 63: Marcio Schwa Ab

56

número maior de pontos avaliados por iteração (ou de indivíduos) deve ser utilizado, já

que o número de sucessos é maior.

Ainda é importante observar os valores altos do desvio padrão do número de

iterações das corridas que resultaram em sucesso, o que salienta o caráter aleatório do

método. Esta aleatoriedade pode dificultar o controle da busca, já que é difícil

determinar quando uma busca pode ser encerrada. Uma forma de se conseguir isto é

diminuir a probabilidade de mutação ao longo das iterações (HIBBERT, 1993b),

fazendo com que o conjunto de pontos convirja, já que é a mutação que mantém uma

certa diversidade no conjunto de pontos.

A avaliação do método do Algoritmo Genético segue com a minimização da

função de Levy. Os resultados obtidos são apresentados nas Tabelas 10, 11, 12 e 13.

Os resultados obtidos pelo Algoritmo Genético na minimização da função de

Levy mostram que o número de sucessos aumenta com o aumento da probabilidade de

mutação, independentemente do tamanho da perturbação realizada nas mutações.

Entretanto, quando a perturbação realizada pela mutação tem um tamanho de 1 e 2% da

faixa total de busca (Tabelas 10 e 11, respectivamente), o número de sucessos é baixo.

Somente com uma perturbação correspondente a 5 e 10% da faixa total de busca é

possível atingir uma taxa de sucessos igual ou próxima a 100%.

Tabela 10: Resultados obtidos com o método do Algoritmo Genético na minimização da função de Levy, com perturbação de 1% nas mutações.

Npt Pmut [%] Sucessos Iter DPIter NF 1 2 2591 2597 647625 10 3 87 40 14444 20 4 1165 1298 145563 30 11 6610 7212 300475 40 15 11799 13703 393313

10

50 28 13021 13632 232520 1 8 719 645 134742 10 7 70 22 14969 20 12 3218 11028 402208 30 13 4306 8148 496828 40 25 9702 14306 582134

30

50 29 2965 5281 153355 1 10 801 639 200275 10 17 36 27 5234 20 10 5684 14425 1421050 30 15 4736 12018 789356 40 29 3285 7039 283231

50

50 33 1018 2297 77105

Page 64: Marcio Schwa Ab

57

Tabela 11: Resultados obtidos com o método do Algoritmo Genético na minimização da função de Levy, com perturbação de 2% nas mutações.

Npt Pmut [%] Sucessos Iter DPIter NF 1 0 ----- ----- ----- 10 16 9751 12664 304727 20 29 3948 4674 68076 30 30 868 995 14473 40 30 316 290 5262

10

50 29 218 177 3750 1 10 4369 10773 655275 10 26 8416 12534 485547 20 32 900 1181 42193 30 33 260 482 11826 40 29 111 138 5752

30

50 36 61 39 2530 1 10 1271 1504 317675 10 34 3768 6554 277087 20 38 335 501 22060 30 34 111 110 8188 40 40 62 50 3864

50

50 42 57 28 3374

Tabela 12: Resultados obtidos com o método do Algoritmo Genético na minimização da função de Levy, com perturbação de 5% nas mutações.

Npt Pmut [%] Sucessos Iter DPIter NF 1 34 18302 12854 269152 10 28 1767 1081 31561 20 31 653 436 10531 30 33 1006 3437 15244 40 47 5318 8732 56571

10

50 49 1463 1629 14933 1 20 7009 4160 525709 10 29 501 294 25926 20 35 150 109 6426 30 46 3638 9797 118615 40 49 135 105 4125

30

50 47 159 107 5068 1 31 5982 7695 482393 10 39 432 306 27668 20 39 122 94 7845 30 42 783 4633 46593 40 46 93 47 5065

50

50 48 120 69 6275

Page 65: Marcio Schwa Ab

58

Tabela 13: Resultados obtidos com o método do Algoritmo Genético na minimização da função de Levy, com perturbação de 10% nas mutações.

Npt Pmut [%] Sucessos Iter DPIter NF 1 26 22145 12858 425865 10 29 3162 4798 54510 20 50 5787 10908 57869 30 50 1732 2529 17318 40 50 607 456 6070

10

50 50 743 513 7431 1 21 18019 13099 1287071 10 32 1118 768 52403 20 49 2213 6851 67733 30 50 237 184 7124 40 50 220 146 6600

30

50 50 290 235 8698 1 27 10713 8020 991948 10 38 1416 4880 93161 20 50 1477 3750 73868 30 50 122 76 6094 40 48 145 100 7576

50

50 50 281 189 14053

Já o número médio de iterações diminuiu com o aumento da probabilidade de

mutação, com exceção dos resultados obtidos com 1% de perturbação, onde o número

de iterações apresentou um mínimo em 10% de probabilidade de mutação.

Como o aumento da probabilidade de mutação provoca um aumento no número

de sucessos e, simultaneamente, uma queda no número de iterações, o número esperado

de avaliações da função objetivo também cai com o aumento da probabilidade de

mutação.

Com relação ao número de indivíduos (número de pontos avaliados por

iteração), é observado que um aumento do número de indivíduos leva a um aumento do

número de sucessos e a uma queda do número médio de iterações. Além disso, o

número esperado de avaliações da função objetivo também cai em alguns casos;

enquanto em outros casos, fica constante ou muda pouco, o que indica a possibilidade

de se utilizar um número maior de indivíduos sem, no entanto, levar a um esforço

computacional excessivo.

É interessante observar que, no caso da minimização da função de Levy, uma

configuração mais aleatória foi a que resulta nos melhores resultados, indicando que

funções com muitos mínimos locais, um caráter mais aleatório é necessário para

garantir a realização de uma busca global e conseguir determinar o mínimo global. Já a

Page 66: Marcio Schwa Ab

59

função de Rosenbrock não resulta em bons resultados quando o caráter aleatório da

busca é predominante, já que os melhores resultados foram obtidos com uma pequena

perturbação realizada pela mutação.

Esta comparação mostra que a definição dos parâmetros ótimos de busca é

dependente do problema que se deseja resolver. Desta forma, quando se tem em mãos

um problema cuja solução não é conhecida, a melhor forma de resolvê-lo consiste em

realizar diversas buscas com valores diferentes dos parâmetros de busca, de forma a

aumentar a probabilidade de encontrar o mínimo global deste problema.

O método do Algoritmo Genético foi também avaliado na determinação da

região de confiança dos parâmetros de um modelo linear. Foram avaliadas as

influências do tamanho da perturbação realizada pela mutação e do número de

avaliações por iteração. A probabilidade de cruzamento foi mantida igual a 70%, do

mesmo modo que foi feito na avaliação da minimização das funções-teste, e a

probabilidade de mutação foi mantida igual a 40%, já que esta valor resultou sempre em

bons resultados nas minimizações das funções teste. Foram usadas três configurações e

os resultados obtidos estão apresentados na Tabela 14.

Tabela 14: Condições de busca e resultados obtidos pelo método do Algoritmo Genético na estimação de parâmetros de um modelo linear.

Condição AG01 AG02 AG03 Niter 100 100 200 Npt 10 10 05

∆x [%] 1,0 5,0 5,0 S(θ) 20,79 20,78 20,79 θ1 4,85 4,38 4,83 θ2 5,38 5,43 5,41

PRC 924 523 570

As três configurações acima apresentadas alcançaram o ponto mínimo, já que os

valores obtidos são muito próximos dos valores obtidos pelo ESTIMA. A configuração

onde a perturbação da mutação foi de 1,0% (configuração AG01) obteve um número

alto de pontos na região de confiança (PRC), mas estes pontos ficaram concentrados no

cento da região de confiança, como mostra a Figura 11. O aumento da perturbação para

5,0% (configuração AG02) provocou uma queda no número de pontos na região de

confiança, mas como os pontos se encontraram mais dispersos, a região de confiança foi

descrita com uma qualidade superior, como mostra a Figura 12. Entretanto, os pontos

Page 67: Marcio Schwa Ab

60

estão dispersos na parte superior da região de confiança e a parte inferior ficou sem

pontos.

Uma terceira configuração foi avaliada (configuração AG03), onde foram

avaliados 5 pontos por iteração e o número de iterações foi 200, mantendo-se constante

o número de avaliações da função objetivo. A perturbação da mutação foi mantida igual

a 5%. A região de confiança obtida (Figura 13), apesar de apresentar os pontos mais

distribuídos, ainda obteve mais pontos na região superior da região de confiança.

Os resultados obtidos na determinação da região de confiança com o método do

Algoritmo Genético mostram que um caráter mais aleatório é indicado para uma boa

descrição da região de confiança, apesar de que a região obtida pode apresentar

problemas como os apontados acima. Como no método de Monte Carlo, um maior

número de iterações e de avaliações por iteração seriam necessários para uma boa

descrição da região de confiança. Por outro lado, devido ao caráter aleatório deste e dos

demais algoritmos, um número maior de resultados deve ser analisado em conjunto para

ser possível obterem-se conclusões mais consistentes. Entretanto, mesmo com este

pequeno número de casos avaliados, é possível verificar-se a influência de alguns

parâmetros e mostrar que é possível obter-se uma descrição aproximada da região de

confiança, mesmo com um número pequeno de avaliações da função objetivo.

Figura 11: Comparação da região de confiança obtida pelo método do Algoritmo Genético (AG01) com a região elíptica.

Page 68: Marcio Schwa Ab

61

Figura 12: Comparação da região de confiança obtida pelo método do Algoritmo Genético (AG02) com a região elíptica.

Figura 13: Comparação da região de confiança obtida pelo método do Algoritmo Genético (AG03) com a região elíptica.

Page 69: Marcio Schwa Ab

62

Assim, o método do Algoritmo Genético se mostrou eficiente para minimização

das funções-teste e para obtenção da região de confiança de parâmetros. Além disso,

ficou claro que a escolha adequada dos parâmetros de busca tem grande influência nos

resultados obtidos. Assim, como foi discutido acima, a realização de diversas buscas

com diferentes configurações é necessária para aumentar a probabilidade de se ter

encontrado o mínimo global e para a obtenção de uma boa descrição da região de

confiança.

3.3. Recozimento Simulado

O método do Recozimento Simulado consiste em um procedimento de

minimização (ou maximização) exaustiva, sendo a busca realizada através de pequenas

transições aleatórias que exploram toda a região de busca durante a procura do mínimo

global. No início da minimização são realizadas transições para valores maiores da

função objetivo, permitindo que a busca “escape” de mínimos locais. Ao longo das

iterações, as transições passam gradativamente a serem aceitas somente quando a

função objetivo diminui, para que a busca convirja para o mínimo.

Neste trabalho será utilizada a forma tradicional para o cálculo da probabilidade

de transição; ou seja, aquela baseada na estatística de Boltzmann. A queda da

temperatura será feita exponencialmente (KIRKPATRICK et al., 1983) e o valor inicial

da temperatura será definido de tal forma que inicialmente 95% das transições para

pontos com valor maior da função objetivo sejam aceitas (DAS et al., 1990). Assim, as

equações para o cálculo da probabilidade de transição, para a redução da temperatura e

para a definição da temperatura inicial são, respectivamente:

expN

k

F FpT

− = − (3.13)

1k kT α T+ = ⋅ (3.14)

( )0 ln 0,95

FT∗∆=− (3.15)

onde p é a probabilidade de transição, F é o valor atual da função objetivo, FN é o valor

da função objetivo no novo ponto, Tk é a k-ésima temperatura, α é o fator de redução de

temperatura, T0 é a temperatura inicial, e ∆F* é a máxima variação da função objetivo

em torno do ponto inicial, calculada a partir de um certo número de transições.

Page 70: Marcio Schwa Ab

63

O procedimento iterativo consiste na realização de um certo número de

tentativas de transições Npt, que podem ou não ser aceitas, em uma mesma temperatura

Tk, de acordo com a Equação (3.13). Em seguida a temperatura é reduzida com o auxílio

da Equação (3.14) e são realizadas outras Npt tentativas de transições, reduz-se a

temperatura e assim por diante, até que um critério de parada seja satisfeito.

Apesar de AHON et al. (2000) mostrarem que com a estatística de Tsallis o

método encontra o mínimo global com maior freqüência, esta forma requer um maior

número de avaliações da função objetivo, além da necessidade da definição dos

parâmetros adicionais qv e qa. Assim, a aplicação da estatística de Boltzmann é uma

forma mais simples, já que não depende de muitos parâmetros. Como a temperatura

inicial é definida pela Equação (3.15), os parâmetros que ainda devem ser definidos são

o número de reduções da temperatura (NT) e o número de tentativas de transições feitas

em cada temperatura Npt. Uma alternativa é reduzir a temperatura até que um número

máximo de transições não-aceitas consecutivamente seja atingido, indicando que o

sistema está “congelado”.

O valor da temperatura inicial é determinado a partir da Equação (3.15), sendo o

valor de ∆F* determinado pela diferença entre o maior e o menor valor da função

objetivo encontrados em um número significativo de tentativas de transição a partir de

um ponto inicial sorteado de forma aleatória. Neste trabalho, foi arbitrado um valor de

50 tentativas de transição para a definição da temperatura inicial.

O último ponto importante é a questão da “vizinhança”, a qual define o conjunto

de pontos vizinhos ao ponto atual que podem vir a ser o novo ponto. De uma forma

simples, é definido um intervalo em torno do ponto atual em cada direção de busca, e

um novo ponto é sorteado nesta região. Neste trabalho a região em torno do ponto atual

foi definida como qualquer posição ∆x maior ou menor que a posição atual, sendo ∆x

um vetor contendo as variações em cada direção para definição de uma região em torno

do ponto atual onde o novo ponto será gerado usando a Equação (3.12), também

utilizada pelo método do Algoritmo Genético, reescrita abaixo:

( )2 1N = + ⋅ − ⊗∆x x r x (3.12)

onde xN e x são os vetores correspondentes ao novo ponto e ao ponto atual,

respectivamente e r é um vetor de números aleatórios com distribuição uniforme no

intervalo [0, 1].

Na Figura 14 é apresentado um fluxograma básico do método.

Page 71: Marcio Schwa Ab

64

Figura 14: Fluxograma do método do Recozimento Simulado.

Foram analisadas as influências do número de avaliações em cada temperatura

(Npt), do valor de α, usado na Equação (3.14) para reduzir o valor da temperatura, e da

perturbação realizada para o cálculo do novo ponto, que corresponde a um percentual da

faixa total de cada direção de busca. Os critérios de parada foram os limites da função

objetivo propostos para cada função (100 para a função de Rosenbrock e -176 para a

Page 72: Marcio Schwa Ab

65

função de Levy). A busca era interrompida sempre que após cinco reduções sucessivas

de temperatura nenhuma transição fosse aceita, indicando que a busca “congelou”.

Nas Tabelas 15 e 16 estão reportados os resultados obtidos na minimização da

função de Rosenbrock, com uma perturbação para o cálculo do novo ponto de 0,1 e

0,5% da faixa de busca em cada direção, respectivamente. Uma perturbação de 1%

também foi usada, mas como não foi obtido nenhum sucesso, os dados não foram

reportados.

Os resultados apresentados nas Tabelas 15 e 16 mostram que o método do

Recozimento Simulado alcançou bons resultados na minimização da função de

Rosenbrock, obtendo um número de sucessos alto, particularmente para uma

perturbação de 0,1%.

Os resultado também mostram que o aumento do número de avaliações em uma

mesma temperatura (Npt) provocou um aumento do número de sucessos e uma queda

do número de reduções de temperatura (NT). O número esperado de avaliações da

função objetivo, no entanto, aumenta muito.

Tabela 15: Resultados obtidos com o método do Recozimento Simulado na minimização da função de Rosenbrock, com perturbação de 0,1%.

alfa Npt Sucessos NT DPIter NF 100 41 98 24 11984 500 42 26 5 15334 1000 44 23 2 26679 3000 43 22 1 77637 5000 44 22 1 126808

0,5

10000 47 22 1 234948 100 34 104 28 15281 500 38 46 3 30419 1000 42 44 1 51899 3000 46 44 2 141848 5000 46 44 3 239249

0,7

10000 49 43 2 440650 100 34 162 18 23750 500 40 151 8 94594 1000 49 151 16 153769 3000 47 149 8 475396 5000 46 145 5 789816

0,9

10000 49 146 5 1488755

Page 73: Marcio Schwa Ab

66

Tabela 16: Resultados obtidos com o método do Recozimento Simulado na minimização da função de Rosenbrock, com perturbação de 0,5%.

alfa Npt Sucessos NT DPIter NF 100 7 40 8 28367 500 28 32 8 28476 1000 28 28 3 50829 3000 34 26 1 113668 5000 38 25 1 164647

0,5

10000 41 25 1 302499 100 3 61 1 101667 500 30 55 5 46028 1000 28 53 3 94069 3000 38 50 2 196641 5000 45 49 1 273333

0,7

10000 43 49 1 564359 100 8 177 3 110391 500 34 172 6 126125 1000 40 169 4 211250 3000 39 166 4 639250 5000 40 164 3 1027969

0,9

10000 46 163 3 1770794

O aumento de α não provocou mudanças significativas no número de sucessos,

mas provocou um aumento do número de reduções de temperatura e um conseqüente

aumento do número esperado de avaliações da função objetivo.

O aumento do tamanho da perturbação para as transições provocou uma ligeira

queda, quando os resultados das Tabelas 15 e 16 são comparados. Lembrando que

quando a perturbação foi de 1% nenhum sucesso foi obtido, chega-se à conclusão de

que o aumento da perturbação provoca uma queda do número de sucessos na

minimização da função de Rosenbrock.

Nas Tabelas 17 e 18 são reportados os resultados obtidos na minimização da

função de Levy com 0,5 e 1% de perturbação para o cálculo do novo ponto. O valor de

0,1% para perturbação também foi utilizado, mas, como o número de sucessos foi muito

baixo, os dados não foram reportados.

Desta vez, os resultados apresentados nas Tabelas 17 e 18 mostram que o

método do Recozimento Simulado encontrou uma certa dificuldade na minimização da

função de Levy, já que o número de sucessos não foi alto para a maioria das condições

de busca.

Page 74: Marcio Schwa Ab

67

Tabela 17: Resultados obtidos com o método do Recozimento Simulado na minimização da função de Levy, com perturbação de 0,5%.

alfa Npt Sucessos NT DPIter NF 100 0 ----- ----- ----- 500 0 ----- ----- ----- 1000 1 3 ----- 150000 3000 5 3 1 84000 5000 4 2 1 125000

0,5

10000 12 2 1 86806 100 1 9 ----- 45000 500 2 6 1 68750 1000 2 9 4 225000 3000 6 6 3 141667 5000 14 5 3 95663

0,7

10000 25 4 2 84000 100 0 ----- ----- ----- 500 3 22 8 183333 1000 10 19 10 95000 3000 20 17 9 128625 5000 30 16 9 134444

0,9

10000 42 13 9 149093

Tabela 18: Resultados obtidos com o método do Recozimento Simulado na minimização da função de Levy, com perturbação de 1,0%.

alfa Npt Sucessos NT DPIter NF 100 0 ----- ----- ----- 500 1 3 ----- 75000 1000 4 3 2 37500 3000 13 4 1 43491 5000 19 3 1 44321

0,5

10000 25 3 1 64800 100 0 ----- ----- ----- 500 3 8 2 63889 1000 7 6 5 43878 3000 21 6 2 43537 5000 24 7 3 71615

0,7

10000 34 6 3 90398 100 3 27 11 44444 500 13 30 11 58136 1000 22 24 9 54545 3000 30 19 8 96833 5000 39 21 10 132479

0,9

10000 45 17 9 193333

Os resultados mostram que o aumento do número de avaliações em uma mesma

temperatura (NT) leva a um aumento do número de sucessos e a uma ligeira queda do

Page 75: Marcio Schwa Ab

68

número de reduções de temperatura, mas com um aumento do número esperado de

avaliações da função objetivo, da mesma forma que foi observado para a minimização

da função de Rosenbrock.

O aumento de α provocou um aumento do número de sucessos, mas sendo

necessárias mais reduções de temperatura. Já o aumento do tamanho da perturbação

para as transições provocou um aumento do número de sucessos, enquanto o número de

reduções de temperatura não foi significativamente afetado. Este comportamento foi o

oposto do observado na minimização da função de Rosenbrock.

Os resultados obtidos pelo método do Recozimento Simulado indicam uma certa

dependência dos parâmetros de busca com o tipo de problema. Assim, para uma função

que apresenta diversos mínimos locais, como a função de Levy, é necessário um caráter

mais global da busca, obtido com uma queda lenta da temperatura e um grande número

de avaliações em cada temperatura. Além disso, uma perturbação relativamente grande

para as transições também leva a melhores resultados. Já para a minimização da função

de Rosenbrock, o aumento do caráter global da busca não aumenta muito o número de

sucessos, apesar do aumento do número de reduções de temperatura e do maior número

de avaliações da função objetivo.

É importante ressaltar que estas conclusões são restritas às faixas de parâmetros

de busca avaliadas aqui e que outros comportamentos podem ser observados se outros

parâmetros forem avaliados. Um exemplo é o fato do número de sucessos aumentar com

o aumento da perturbação para o cálculo de um novo ponto. É bem provável que este

aumento esteja restrito à faixa avaliada e que, se este parâmetro continuar a ser

aumentado, o número de sucessos venha a cair.

A avaliação da obtenção da região de confiança dos parâmetros de um modelo

linear pelo método do Recozimento Simulado foi feita a partir de três estimações, onde

os parâmetros de busca foram alterados conforme mostra a Tabela 19, onde também

estão apresentados os resultados obtidos em cada estimação. Os valores altos da

perturbação ∆x foram assim definidos para aumentar o caráter exploratório da busca, e

procurar obter uma boa descrição da região de confiança.

A configuração RS01, onde o fator de queda da temperatura foi 0,9 e a

perturbação foi de 5,0%, obteve uma boa descrição da região de confiança, apesar da

parte superior da região conter poucos pontos, como mostra a Figura 15.

Page 76: Marcio Schwa Ab

69

Tabela 19: Condições de busca e resultados obtidos pelo método do Recozimento Simulado na estimação de parâmetros de um modelo linear.

Condição RS01 RS02 RS03 Niter 100 100 200 Npt 10 10 05 α 0,9 0,7 0,9

∆x [%] 5,0 5,0 10,0 S(θ) 20,80 20,78 20,85 θ1 4,86 4,38 4,84 θ2 5,34 5,42 5,47

PRC 314 387 141

Figura 15: Comparação da região de confiança obtida pelo método do Recozimento Simulado (RS01) com a região elíptica.

A diminuição do fator de redução da temperatura para 0,7 (configuração RS02)

concentra os pontos no centro da região de confiança. A rápida redução da temperatura

concentra a busca próxima ao ponto mínimo, sem que seja realizada uma boa varredura

de toda a região de busca, de forma que a parte inferior da mesma não contém nenhum

ponto, como apresentado na Figura 16.

Page 77: Marcio Schwa Ab

70

Figura 16: Comparação da região de confiança obtida pelo método do Recozimento Simulado (RS02) com a região elíptica.

Figura 17: Comparação da região de confiança obtida pelo método do Recozimento Simulado (RS03) com a região elíptica.

Page 78: Marcio Schwa Ab

71

Assim, a configuração RS03 manteve o valor do fator de queda da temperatura

em 0,9 e aumentou o número de temperaturas para 200 e diminuiu o número de

avaliações por temperatura para 5. Entretanto, o valor da perturbação teve de ser

aumentado, pois como o número de quedas na temperatura é maior, é necessária uma

maior perturbação para evitar que os pontos fiquem muito concentrados no centro da

região de confiança. A Figura 17 mostra a região de confiança obtida. Apesar do

pequeno número de pontos, a descrição da região é boa, já que os pontos se encontram

bem espalhados.

Os resultados obtidos na determinação da região de confiança pelo método do

Recozimento Simulado mostraram que é importante a definição adequada dos

parâmetros de busca para a obtenção de regiões de confiança com boa qualidade. Em

particular, uma queda lenta na temperatura tende a gerar bons resultados.

Desta forma, o método do Recozimento Simulado resultou em bons resultados

para a minimização das funções teste e para definição da região de confiança de

parâmetros de um modelo linear. Entretanto, na etapa da avaliação da minimização de

funções-teste ficou clara a dependência dos parâmetros de busca com o tipo de

problema. Sendo assim, é necessária a realização de diversas buscas, alterando-se o

valor dos parâmetros, para que se aumente a probabilidade de ser encontrado o mínimo

global e se obtenha uma região de confiança com uma boa qualidade.

3.4. Enxame de Partículas

O Método do Enxame de Partículas proposto por KENNEDY e EBERHART

(1995) consiste em um método de otimização simples e robusto onde um conjunto de

pontos, denominados partículas, trocam informações de forma a levar o conjunto de

partículas para o ponto ótimo de uma função qualquer. A partir deste primeiro trabalho,

algumas modificações foram propostas, como a inserção do peso de inércia (SHI e

EBERHART, 1998a) e a inserção do fator de constrição (CLERC, 1999), além de uma

nova formulação proposta recentemente por BISCAIA et al. (2004), conforme foi

discutido no capítulo anterior. Neste trabalho foi utilizado o método do Enxame de

Partículas com a inserção do peso de inércia, de acordo com as seguintes equações:

( ) ( )11 2

k k k k k ki i 1 i i 2 global iw c c+ = ⋅ + ⋅ ⊗ − + ⋅ ⊗ −v v r p x r p x (3.13a)

1 1k k ki i i+ += +x x v (3.13b)

Page 79: Marcio Schwa Ab

72

onde xi e vi são, respectivamente, vetores da posição e velocidade da partícula i, w é o

peso de inércia, c1 e c2 são duas constantes, r1 e r2 são dois vetores contendo números

aleatórios com distribuição uniforme no intervalo [0, 1], pi é a posição com o menor

valor da função objetivo que a partícula i já encontrou e pglobal é o posição do menor

valor encontrado por todo o conjunto de partículas. O símbolo ⊗ indica multiplicação

de elemento por elemento dos vetores.

O peso de inércia w, segundo TRELEA (2003), deve ser um valor menor que 1

para que a busca convirja para um ponto, a despeito dos valores de c1 e c2. Outra forma

de realizar a busca é diminuir o peso de inércia ao longo das iterações, iniciando a busca

com um valor alto, aumentando o caráter exploratório da busca, e finalizando com um

valor baixo, aumentando o caráter local da busca e permitindo o refinamento da solução

encontrada. Uma forma simples de se fazer isto consiste em reduzir o valor de w

linearmente de um valor inicial wi até um valor final wf, de acordo com a seguinte

equação:

( ) ( )( )

11i i f

kw w w w

Niter−

= + − ⋅−

(3.14)

onde wi e wf são, respectivamente, os valores do peso de inércia no início e na última

iteração Niter, sendo Niter o número máximo de iterações. O valor de Niter também

pode ser definido como a iteração a partir da qual o valor do peso de inércia permaneça

constante em wf. Definindo wi igual a wf, a busca é feita com um valor constante do peso

de inércia w. Para ajudar a evitar a “explosão” do conjunto de partículas, foi definida

uma velocidade máxima em cada direção de busca, sendo esta velocidade igual a

metade do tamanho total de cada direção de busca.

( )max 2H Lv x x= − (3.15)

Além disso, sempre que uma partícula ultrapassar algum limite da região de

busca, esta é recolocada no limite da região e a sua velocidade é definida como a

metade de sua velocidade atual com o sentido invertido; ou seja, sua velocidade é

multiplicada por -1/2. A posição inicial de cada partícula é calculada utilizando a

Equação (3.8) (ou seja, de forma aleatória em toda região de busca) e a velocidade

inicial é calculada de forma aleatória segundo a seguinte equação:

( )0max2 1i = ⋅ − ⋅v r v (3.16)

Page 80: Marcio Schwa Ab

73

Na Figura 18 o fluxograma do método do Enxame de Partículas é apresentado.

Figura 18: Fluxograma do método do Enxame de Partículas.

Como os demais métodos, este foi avaliado na minimização das funções de

Rosenbrock e de Levy, onde alguns valores de w, c1 e c2 foram avaliados, verificando a

influência destes no desempenho do método. Os critérios de parada foram o valor

definido como mínimo para cada função (100 para a função de Rosenbrock e -176 para

a função de Levy), o limite máximo de 50000 iterações e uma análise da convergência

das partículas, de acordo com a seguinte equação:

Page 81: Marcio Schwa Ab

74

( )1

1 Npt

i otmi

F FNpt

ε=

− <∑ (3.17)

onde Fi é o valor da função objetivo da partícula i, Fotm é o melhor valor encontrado e ε

é uma tolerância fixada em 10-5. Sempre que a Equação (3.17) é satisfeita, considera-se

que o conjunto de partículas convergiu e a busca é encerrada.

Nas Tabelas 20, 21, 22 e 23 são apresentados os resultados obtidos com os

valores de c1 e c2 iguais a 0,5 1,0 1,5 e 2,0, respectivamente, para minimização da

função de Rosenbrock. Os valores de c1 e c2 usados neste trabalho sempre são iguais

entre si, para manter um equilíbrio entre a atração na direção do melhor ponto

encontrado pela partícula e a atração na direção do melhor ponto encontrado por todo

enxame. Em cada tabela foram usados diferentes valores de w, ou mantidos fixos

durante a busca, ou variando de 1,2 a 0,5, conforme a Equação (3.14). Também foi

avaliado o efeito do número de partículas, sendo realizadas buscas com 10, 30 e 50

partículas.

Tabela 20: Resultados obtidos com o método do Enxame de Partículas na minimização da função de Rosenbrock, com os valores de c1 e c2 iguais a 0,5.

Npt w Sucessos Iter DPIter NF 0,8 0 ----- ----- ----- 0,9 0 ----- ----- ----- 10

1,2 - 0,5 0 ----- ----- ----- 0,8 0 ----- ----- ----- 0,9 39 464 240 17841 30

1,2 - 0,5 0 ----- ----- ----- 0,8 0 ----- ----- ----- 0,9 49 325 222 16560 50

1,2 - 0,5 6 550 26 229167

Tabela 21: Resultados obtidos com o método do Enxame de Partículas na minimização da função de Rosenbrock, com os valores de c1 e c2 iguais a 1,0.

Npt w Sucessos Iter DPIter NF 0,8 0 ----- ----- ----- 0,9 50 3850 5766 38501 10

1,2 - 0,5 0 ----- ----- ----- 0,8 31 417 216 20165 0,9 50 2998 4500 89953 30

1,2 - 0,5 1 633 ----- 949500 0,8 43 315 220 18295 0,9 50 1289 846 64461 50

1,2 - 0,5 23 647 23 70317

Page 82: Marcio Schwa Ab

75

Tabela 22: Resultados obtidos com o método do Enxame de Partículas na minimização da função de Rosenbrock, com os valores de c1 e c2 iguais a 1,5.

Npt w Sucessos Iter DPIter NF 0,6 0 ----- ----- ----- 0,7 1 910 ----- 455000 0,8 50 2274 3778 22741 0,9 0 ----- ----- -----

10

1,2 - 0,5 0 ----- ----- ----- 0,6 1 196 ----- 294000 0,7 50 410 262 12295 0,8 50 1010 546 30293 0,9 0 ----- ----- -----

30

1,2 - 0,5 25 829 34 49726 0,6 13 209 130 40163 0,7 50 327 286 16355 0,8 50 1312 1570 65602 0,9 0 ----- ----- -----

50

1,2 - 0,5 36 823 75 57176

Tabela 23: Resultados obtidos com o método do Enxame de Partículas na minimização da função de Rosenbrock, com os valores de c1 e c2 iguais a 2,0.

Npt w Sucessos Iter DPIter NF 0,1 0 ----- ----- ----- 0,2 0 ----- ----- ----- 0,3 0 ----- ----- ----- 0,4 10 620 192 30975 0,5 50 1837 2947 18374 0,6 50 2681 2960 26815 0,7 48 11124 6945 115870

10

1,2 - 0,5 49 2810 2667 28669 0,1 5 580 321 174000 0,2 5 369 102 110700 0,3 35 346 122 14831 0,4 50 478 393 14340 0,5 50 827 1135 24800 0,6 50 1347 1220 40420 0,7 49 7069 5622 216409

30

1,2 - 0,5 50 1632 1107 48968 0,1 21 434 178 51718 0,2 25 459 711 45920 0,3 48 397 555 20654 0,4 50 507 555 25364 0,5 50 701 952 35073 0,6 50 1004 755 50180 0,7 50 6658 3598 332920

50

1,2 - 0,5 50 1508 1040 75412

Page 83: Marcio Schwa Ab

76

Os resultados obtidos com o valor de c1 e c2 igual a 0,5, apresentados na Tabela

20, mostram que somente foram obtidos sucessos em três configurações, sendo que as

configurações que não obtiveram sucesso algum foram excluídas da tabela. Com 30

partículas e w igual a 0,9 o número de sucessos não foi alto, mas o número médio das

iterações foi baixo, de forma que o número esperado de avaliações da função objetivo se

manteve próximo aos menores valores obtidos. Já com 50 partículas e w igual a 0,9, o

número de sucessos foi alto e houve uma queda do número de iterações com relação à

busca com 30 partículas, mantendo baixo o número esperado de avaliações da função

objetivo. Já a busca realizada com 50 partículas e w variando de 1,2 a 0,5 obteve um

número de sucessos muito baixo.

Com o aumento de valor de c1 e c2 para 1,0 a quantidade de configurações que

obtiveram sucesso foi maior, conforme a Tabela 21, e sempre com um valor alto de w.

Com o valor de 0,9 para w foram obtidos 50 sucessos, ou seja, 100% de sucesso, com

10, 30 e 50 partículas. Entretanto o número médio de iterações foi alto, de forma a

aumentar o número esperado de avaliações da função objetivo. Com o valor de w igual a

0,8 o número de sucessos foi menor e foram obtidos com 30 e 50 partículas. Mas o

número médio de iterações foi bem menor que as buscas realizadas com w igual a 0,9,

fazendo com que o número esperado de avaliações da função objetivo fosse menor,

mesmo com um menor número de sucessos.

Na Tabela 22, onde o valor de c1 e c2 é igual a 1,5, o número de configurações

que alcançam algum sucesso aumenta, mas a faixa de w que consegue obter uma taxa de

sucessos próximo a 100% ainda é estreita e se limita aos valores de w iguais a 0,7 e 0,8.

As buscas realizadas com w igual a 0,9, que obtiverem bons resultados com valores

menores de c1 e c2, não obtiveram sucesso algum, já que atingiram o limite máximo de

iterações antes de ser alcançado o sucesso. Com 30 e 50 partículas foram obtidos os

melhores resultados com o valor de w igual a 0,7, já que tem uma taxa de sucessos de

100% e um número médio de iterações menor que as buscas com w igual a 0,8 e que

também obtiveram 100% de sucesso.

Quando o valor de c1 e c2 é igual a 2,0, o número de configurações que alcançam

sucessos aumenta ainda mais, conforme a Tabela 23. Porém, agora a faixa dos valores

de w que alcançam os melhores resultados é de 0,4 a 0,7. As corridas com w igual a 0,8

e 0,9 atingiram o limite máximo de iterações e as corridas com valores baixos de w têm

pouco ou nenhum sucesso devido à convergência prematura do enxame.

Page 84: Marcio Schwa Ab

77

Desta forma, ficou claro que para a minimização da função de Rosenbrock é

importante uma definição adequada dos parâmetros de busca. Uma vez que estes

parâmetros sejam escolhidos corretamente, os resultados obtidos pelo método do

Enxame de Partículas são muito bons, alcançando uma taxa de sucesso de 100% e com

um número médio de iterações relativamente baixo.

Nas Tabelas 24, 25, 26 e 27 apresentadas a seguir, o método do Enxame de

Partículas é avaliado na minimização da função de Levy, a fim de verificar como este

método se comporta na minimização de uma função com diversos mínimos locais.

Tabela 24: Resultados obtidos com o método do Enxame de Partículas na minimização da função de Levy, com os valores de c1 e c2 iguais a 0,5.

Npt w Sucessos Iter DPIter NF 0,1 7 32 16 2286 0,2 4 25 12 3094 0,3 9 49 44 2728 0,4 17 42 32 1237 0,5 25 18 5 370 0,6 26 32 38 615 0,7 25 34 25 671 0,8 32 57 51 897 0,9 42 105 53 1250

10

1,2 – 0,5 49 468 24 4776 0,1 15 35 20 3453 0,2 20 26 17 1980 0,3 21 18 12 1259 0,4 26 16 15 908 0,5 22 15 3 1026 0,6 32 25 33 1150 0,7 38 53 60 2095 0,8 48 59 78 1842 0,9 50 72 32 2150

30

1,2 – 0,5 50 445 18 13363 0,1 27 25 18 2342 0,2 27 17 7 1571 0,3 29 14 6 1186 0,4 34 14 8 1062 0,5 33 14 3 1072 0,6 34 22 21 1581 0,7 45 44 48 2430 0,8 49 48 48 2436 0,9 50 55 16 2750

50

1,2 – 0,5 50 434 25 21690

Page 85: Marcio Schwa Ab

78

Na Tabela 24, onde c1 e c2 são iguais a 0,5, o número de sucessos só alcança

valores altos quando os valores de w também são altos. Para valores baixos de w o

número de sucessos é baixo, mas o número médio de iterações é baixo, de forma que

nessas condições o número esperado de avaliações da função objetivo é baixo. Assim,

de acordo com este dado, valores de w próximos a 0,5 resultam nos melhores resultados.

No caso do custo computacional na avaliação da função objetivo não limitar a busca,

valores de w iguais a 0,9 devem ser utilizados, já que resultam num maior número de

sucessos.

Tabela 25: Resultados obtidos com o método do Enxame de Partículas na minimização da função de Levy, com os valores de c1 e c2 iguais a 1,0.

Npt w Sucessos Iter DPIter NF 0,1 5 17 4 1680 0,2 18 20 11 546 0,3 18 14 6 401 0,4 24 16 5 332 0,5 28 23 10 405 0,6 27 27 18 508 0,7 32 45 35 698 0,8 45 77 47 859 0,9 49 169 114 1726

10

1,2 – 0,5 50 511 35 5115 0,1 31 11 7 535 0,2 29 12 8 626 0,3 33 17 13 763 0,4 32 16 10 749 0,5 43 24 21 840 0,6 42 33 30 1195 0,7 47 32 15 1027 0,8 46 50 22 1632 0,9 50 94 39 2824

30

1,2 – 0,5 50 465 42 13957 0,1 33 12 11 939 0,2 35 13 15 914 0,3 35 19 18 1351 0,4 42 16 13 979 0,5 39 19 13 1200 0,6 49 34 33 1739 0,7 48 25 8 1319 0,8 50 40 14 1976 0,9 50 76 26 3809

50

1,2 – 0,5 50 463 18 23153

Page 86: Marcio Schwa Ab

79

Com o aumento dos valores de c1 e c2 para 1,0 (Tabela 25), o número de

configurações que resultam num alto número de sucessos aumenta. Como o número

médio de iterações cai um pouco em relação aos resultados obtidos com c1 e c2 iguais a

0,5, o número esperado de avaliações da função objetivo é menor. Aumentando ainda

mais os valores de c1 e c2 para 1,5 (Tabela 26) e 2,0 (Tabela 27), o número de sucessos

aumenta ainda mais, mostrando que o aumento de c1 e c2 aumenta o caráter exploratório

da busca, mas com um maior custo computacional, já que o número esperado de

avaliações da função objetivo aumenta.

Tabela 26: Resultados obtidos com o método do Enxame de Partículas na minimização da função de Levy, com os valores de c1 e c2 iguais a 1,5.

Npt w Sucessos Iter DPIter NF 0,1 13 20 11 751 0,2 28 23 16 415 0,3 25 28 35 561 0,4 29 32 27 552 0,5 35 32 19 452 0,6 39 58 50 746 0,7 43 61 39 712 0,8 50 117 159 1168 0,9 50 332 168 3325

10

1,2 – 0,5 50 570 157 5703 0,1 32 24 39 1124 0,2 46 37 61 1194 0,3 40 29 36 1083 0,4 43 30 33 1054 0,5 45 28 22 942 0,6 48 32 20 992 0,7 50 37 16 1119 0,8 50 64 29 1920 0,9 50 167 75 5017

30

1,2 – 0,5 50 509 30 15273 0,1 46 23 27 1249 0,2 45 20 22 1089 0,3 47 24 31 1297 0,4 50 20 15 1004 0,5 50 25 16 1262 0,6 50 25 12 1226 0,7 50 33 12 1656 0,8 50 51 18 2559 0,9 50 133 58 6627

50

1,2 – 0,5 50 493 40 24650

Page 87: Marcio Schwa Ab

80

Tabela 27: Resultados obtidos com o método do Enxame de Partículas na minimização da função de Levy, com os valores de c1 e c2 iguais a 2,0.

Npt w Sucessos Iter DPIter NF 0,1 37 44 52 592 0,2 45 47 57 524 0,3 40 52 41 649 0,4 45 51 50 565 0,5 44 71 76 806 0,6 46 75 93 819 0,7 47 105 50 1116 0,8 50 168 80 1677 0,9 50 3253 2605 32528

10

1,2 – 0,5 50 640 63 6403 0,1 50 33 41 980 0,2 50 29 26 873 0,3 50 25 14 762 0,4 50 32 26 955 0,5 50 36 19 1085 0,6 50 41 14 1218 0,7 50 56 26 1680 0,8 50 98 43 2927 0,9 50 1300 1141 39007

30

1,2 – 0,5 50 588 35 17645 0,1 50 28 35 1379 0,2 50 22 18 1090 0,3 50 21 17 1045 0,4 50 23 10 1174 0,5 50 32 17 1589 0,6 50 35 10 1770 0,7 50 48 17 2375 0,8 50 86 33 4306 0,9 50 786 670 39301

50

1,2 – 0,5 50 562 40 28117

De uma forma geral, o método do Enxame de Partículas se mostrou eficiente e

robusto na minimização da função de Levy, conseguindo obter uma alta taxa de

sucessos com um número relativamente baixo de iterações. Ainda é importante ressaltar

que no caso da função de Levy os resultados não se mostraram tão dependentes da

definição dos parâmetros de busca como na minimização da função de Rosenbrock,

provavelmente devido ao grande número de variáveis de busca presentes na

minimização da função de Rosenbrock.

Entretanto, apesar da maior dependência da definição dos parâmetros de busca

que o método apresentou na minimização de função de Rosenbrock, comportamentos

idênticos foram observados. O primeiro deles é que o aumento de w provoca um

Page 88: Marcio Schwa Ab

81

aumento do número de sucessos e do número médio de iterações. Isso ocorre devido à

mudança do caráter da busca, já que valores baixos de w provocam uma convergência

mais rápida, diminuindo o número de iterações, mas podendo levar o enxame a

convergir em um ponto ainda distante do mínimo global, o que diminui o número de

sucessos. Já com valores altos de w, a busca tem uma convergência mais lenta. Como

são feitas mais avaliações da função objetivo, a probabilidade de ser encontrado o

mínimo aumenta, o que leva a um maior número de sucessos, mas com um maior

número de iterações. Algumas vezes, a convergência é tão lenta que o limite máximo de

iterações é alcançado e a busca é terminada. Este efeito do parâmetro de busca w, o peso

de inércia, já foi descrito pela literatura (SHI e EBERHART, 1998a,b) e fica claro na

minimização da função de Levy. Na minimização da função de Rosenbrock este efeito

não é tão claro devido à estreita faixa de valores de w que levaram a algum sucesso, mas

também está presente, já que as corridas com baixos valores de w não resultam em

sucessos devido à rápida convergência. As corridas com valores altos de w não

conseguem sucesso devido ao limite máximo de iterações que é alcançado.

Também é importante observar a correlação que existe entre os valores de c1 e c2

e o valor de w. Quando os valores de c1 e c2 são pequenos, valores altos de w fornecem

os melhores resultados. Já com valores altos de c1 e c2, valores pequenos de w tendem a

fornecer resultados melhores. Assim, além da influência de w na convergência do

enxame, os valores de c1 e c2 também modificam o comportamento do movimento das

partículas. Esta correlação é praticamente a mesma para as duas funções usadas na

avaliação do algoritmo. Assim, a definição dos parâmetros não parece ser dependente

do tipo de problema que se está resolvendo. Esta é uma característica importante e

facilita a definição dos parâmetros de busca.

Um ponto que pode trazer alguma dificuldade para a definição dos parâmetros é

se deve ser dada maior importância para o número de sucessos, o número médio de

iterações ou ainda o número esperado de avaliações da função objetivo, que de certa

forma usa as informações do número de sucessos e o número de iterações

simultaneamente. Quando o tempo computacional está sendo considerado, o número de

iterações e o número esperado de avaliações da função objetivo devem ser analisados, já

que estes valores estão diretamente ligados ao tempo que será necessário para solução

do problema. Desta forma, os valores escolhidos para os parâmetros de busca irão levar

a uma rápida convergência, pois nessas condições o número médio de iterações ou o

valor esperado de avaliações da função objetivo é menor. Entretanto, nessas mesmas

Page 89: Marcio Schwa Ab

82

condições a taxa de sucessos geralmente não é alta e devem ser realizadas algumas

corridas para que se tenha maior confiança nos resultados obtidos.

Por outro lado, se o tempo computacional não está sendo considerado ou não

limita a busca, o número (ou a taxa) de sucessos é que deve ser levado em consideração.

Deste modo, condições onde ocorre uma grande exploração da região de busca devem

ser usadas. Um número menor de corridas é necessário, já que existe uma maior

exploração da região de busca em cada corrida e o número de sucessos nestas condições

é alto.

Além disso, um dos objetivos deste trabalho é a solução de problemas de

estimação de parâmetros, onde se planeja determinar a região de confiança em torno do

ponto mínimo utilizando o próprio método de otimização. Para que esta região em torno

do mínimo seja bem determinada, é necessário que as partículas explorem bem a região

de busca, o que também aumenta a taxa de sucessos.

Assim, na avaliação da determinação da região de confiança dos parâmetros de

um modelo linear, como proposto no início deste capítulo, foram usadas três

configurações de busca onde o parâmetro de busca w foi mantido igual a 0,9,

aumentando o caráter exploratório da busca. Os parâmetros c1 e c2 foram alterados, bem

como o número de pontos avaliados por iteração, como mostra a Tabela 28, onde

também são apresentados os resultados de cada estimação.

Tabela 28: Condições de busca e resultados obtidos pelo método do Enxame de Partículas na estimação de parâmetros de um modelo linear.

Condição EP01 EP02 EP03 Niter 100 100 200 Npt 10 10 05

c1, c2 1,0 2,0 2,0 S(θ) 20,78 20,79 20,80 θ1 4,85 4,84 4,83 θ2 5,38 5,42 5,43

PRC 245 125 155

A configuração EP01 obteve uma boa descrição da função objetivo, apesar dos

pontos se concentrarem um pouco no centro da região de confiança, como mostra a

Figura 19. O aumento do valor de c1 e c2 para 2,0 (configuração EP02) aumenta o

caráter exploratório da busca e, como mostra a Figura 20, a região de confiança obtida

apresenta os pontos um pouco mais espalhados por toda a região de confiança.

Page 90: Marcio Schwa Ab

83

Entretanto, o número de pontos é pequeno e isto diminui a qualidade da região de

confiança obtida.

Figura 19: Comparação da região de confiança obtida pelo método do Enxame de Partículas (EP01) com a região elíptica.

Figura 20: Comparação da região de confiança obtida pelo método do Enxame de Partículas (EP02) com a região elíptica.

Page 91: Marcio Schwa Ab

84

Figura 21: Comparação da região de confiança obtida pelo método do Enxame de Partículas (EP03) com a região elíptica.

A terceira configuração usada manteve os valores de c1 e c2 iguais a 2,0, mas

aumentou o número de iterações para 200 e diminuiu o número de pontos avaliados por

iteração para 5. A qualidade da região de confiança obtida é boa, já que os pontos se

encontram bem distribuídos por toda a região de confiança, como mostra a Figura 21.

Como nas três configurações utilizadas acima o valor do peso de inércia w foi

sempre igual a 0,9, as regiões de confiança sempre apresentaram uma boa dispersão, já

que um valor alto para este parâmetro garante um boa exploração da região de busca. É

importante observar que as minimizações com valores altos de w resultam em maior

número de iterações para que o mínimo seja encontrado, isto devido à maior exploração

da região de busca e, conseqüentemente, à maior demora na convergência dos pontos.

Por outro lado, valores altos de w aumentam a probabilidade de ser encontrado o

mínimo global, como pode ser observado pelo alto número de sucessos obtidos com

estes valores de w.

Assim, na resolução de um problema qualquer, cuja solução não é conhecida, a

configuração que deve ser adotada é realizar buscas com valores altos de w, aumentando

a probabilidade de ser encontrado o mínimo global e possibilitando a obtenção de

regiões de confiança com boa qualidade.

Page 92: Marcio Schwa Ab

85

Por fim, uma importante qualidade do método do Enxame de Partículas que

pôde ser observada é a facilidade em alterar a característica da busca. O próprio método

já apresenta uma característica global no inicio da busca, explorando a região de busca,

e ao longo das iterações convergindo as partículas e tornando a busca com característica

local. Além desta característica intrínseca do método, a manipulação dos parâmetros de

busca permite acelerar ou retardar a convergência do algoritmo, permitindo adaptar o

método a uma variedade muito grande de problemas.

3.5. Conclusões

Os quatro métodos heurísticos de otimização considerados neste trabalho foram

avaliados na minimização de duas funções conhecidas: a função de Rosenbrock e a

função de Levy. Todo os métodos foram capazes de solucionar os problemas e foi

possível identificar os pontos fortes e fracos de cada método, bem como a influência dos

parâmetros de busca. Também foram avaliadas as regiões de confiança obtidas pelos

métodos de otimização na estimação dos parâmetros de um modelo linear, mostrando

que é possível obter-se uma boa descrição da região de confiança utilizando os próprios

pontos que foram usados para a minimização da função objetivo.

O método de Monte Carlo obteve resultados bons, alcançando o número máximo

de sucessos sem, no entanto, ser necessário um número muito grande de iterações.

Porém, quando não é aplicada a redução da região de busca, o número de sucessos é

nulo. Por outro lado, quando o ponto mínimo se encontra próximo dos limites da região

de busca, a redução da região de busca acaba excluindo o ponto mínimo,

impossibilitando que a busca alcance o sucesso. Quando a redução da região de busca é

muito rápida, a região de confiança obtida apresenta os pontos muito concentrados em

torno do mínimo, prejudicando a qualidade da região obtida. Deste modo, em problemas

cuja solução não é conhecida, a aplicação deste método pode levar a resultados

indesejados e um cuidado especial deve ser tomado para a definição da taxa de redução

da região de busca.

O método do Algoritmo Genético também obteve resultados bons, alcançando

um alto número de sucessos para uma grande variedade de parâmetros de busca.

Entretanto, os parâmetros de busca se mostraram dependentes do tipo de problema que

se está resolvendo e este ponto pode prejudicar a definição destes valores. É importante

observar que as probabilidades de mutação usadas neste trabalho e que apresentaram os

melhores resultados são muito maiores que os valores encontrados na literatura,

Page 93: Marcio Schwa Ab

86

provavelmente devido à forma com que o cruzamento é realizado, pois no cruzamento

os pontos apenas trocam os valores dos parâmetros e é com a mutação que novos

valores são gerados. Por outro lado, o número de avaliações da função objetivo

realizadas pelo método do Algoritmo Genético é, em geral, maior que o número de

avaliações realizadas pelo método de Monte Carlo. Entretanto, como o método do

Algoritmo Genético não utiliza a redução da região de busca, este não apresenta as

mesmas dificuldades encontradas pelo método de Monte Carlo quando o mínimo se

encontra próximo dos limites da região de busca, o que o torna um método mais

confiável. Com relação à determinação da região de confiança de parâmetros, foi

possível observar que a qualidade da região depende da definição apropriada dos

parâmetros de busca e um número maior de avaliações da função objetivo deve ser

usado para a obtenção de uma região de confiança com maior qualidade.

Já o método do Recozimento Simulado apresentou dificuldades na minimização

da função de Levy, onde os resultados bons se restringiram a algumas condições de

busca. Por outro lado, na minimização da função de Rosenbrock os resultados obtidos

foram muito bons. Os parâmetros de busca que levaram aos melhores resultados se

mostraram dependentes do tipo de problema, como foi observado para o método do

Algoritmo Genético. Com relação ao número de avaliações da função objetivo, o

método do Recozimento Simulado realiza, em geral, um número de avaliações maior

que o Algoritmo Genético. Já as regiões de confiança obtidas, apesar de dependerem da

escolha adequada dos parâmetros de busca, apresentaram uma boa qualidade, já que os

pontos se encontraram bem dispersos por toda a região de confiança.

O método do Enxame de Partículas obteve os melhores resultados em relação

aos demais métodos, já que consegue alcançar o número máximo de sucessos com um

número de iterações e avaliações da função objetivo inferior aos demais métodos, sendo

assim um método mais robusto e eficiente. Entretanto, é importante definir os

parâmetros de busca de forma a obter uma maior eficiência deste método. É interessante

observar que a definição dos parâmetros de busca não aparenta ser dependente do tipo

do problema que se está resolvendo. Entretanto, os parâmetros parecem estar

correlacionados, já que quanto maior o valor dos parâmetros c1 e c2, menor é o valor de

w que alcança os melhores resultados. Com relação à determinação da região de

confiança dos parâmetros do modelo linear, foi observado que as regiões de confiança

apresentam uma boa qualidade, com os pontos bem distribuídos por toda a região.

Page 94: Marcio Schwa Ab

87

Entretanto, é necessário garantir que a definição dos parâmetros de busca seja traduzida

em uma busca com um caráter global, explorando bem toda a região de busca.

Desta forma, os resultados apresentados acima apontam para o método do

Enxame de Partículas como o mais indicado para a solução de problemas de otimização

de grande dimensão e que apresentam mínimos locais e alta correlação entre as

variáveis, problemas geralmente encontrados em estimação de parâmetros. Além disso,

a boa qualidade das regiões de confiança obtidas com este método reforça a indicação

deste método como o mais adequado para problemas de estimação de parâmetros.

Page 95: Marcio Schwa Ab

88

CAPÍTULO 4: PROBLEMAS DE ESTIMAÇÃO DE PARÂMETROS

Neste capítulo são apresentados os resultados obtidos na estimação de

parâmetros de dois problemas reais. O primeiro consiste na estimação dos parâmetros

para o ajuste de um modelo da taxa de consumo de eteno em uma reação de

polimerização de eteno. No segundo problema é feita a estimação de parâmetros de

modelos cinéticos da conversão de gás natural em gás de síntese através da combinação

das reações de reforma com CO2 e oxidação parcial.

Entretanto, antes de partir para a resolução dos problemas reais de estimação de

parâmetros, serão resolvidos problemas de estimação de parâmetros de modelos não-

lineares simples, cujas soluções já foram obtidas por outros métodos. Estes modelos

apresentam regiões de confiança onde a aproximação elíptica é ruim, reforçando a

necessidade de métodos que possibilitam a obtenção da região de confiança verdadeira.

Além disso, neste capítulo o método do Enxame de Partículas será

preferencialmente usado, pois obteve até aqui os melhores resultados conforme o

capitulo anterior.

4.1. Determinação da Região de Confiança em Modelos Não-Lineares

Como foi discutida no Capítulo 2, a determinação da região de confiança dos

parâmetros como uma região elíptica em torno dos parâmetros está associada à

linearidade do modelo, cujos parâmetros estão sendo estimados. Porém, em modelos

não-lineares a região de confiança dos parâmetros pode ter uma forma muito diferente

de uma elipse.

Desta forma, a obtenção da região de confiança dos parâmetros de um modelo

não-linear deve ser feita através de outros métodos. Um método que apresenta bons

resultados é o da região de verossimilhança (Equação (2.39)). Apesar da aparente

dificuldade na construção desta região, a utilização dos métodos heurísticos facilita a

obtenção da mesma, já que os pontos utilizados para a minimização da função objetivo

podem ser empregados na construção da região de confiança, sem implicar em um custo

computacional adicional.

A seguir, serão apresentados os resultados obtidos na determinação da região de

confiança dos parâmetros de modelos não-lineares simples, mostrando as diferenças

ente a região de confiança elíptica e a de verossimilhança, além da facilidade na

Page 96: Marcio Schwa Ab

89

obtenção das regiões de confiança de verossimilhança quando associadas a um método

heurístico de otimização.

O método utilizado para solução destes problemas é o método do Enxame de

Partículas, devido à melhor qualidade dos resultados obtidos, conforme discutido no

capítulo anterior. É importante ressaltar que o nível de confiança usado para a

construção das regiões de confiança foi sempre igual a 95% em todos os problemas

resolvidos neste capítulo e que as regiões elípticas foram construídas da mesma forma

como foram obtidas para o modelo linear no capítulo anterior, onde este procedimento

foi descrito com detalhes.

4.1.1. Modelo de Michaelis-Menten

O modelo de Michaelis-Menten é muito aplicado para representação matemática

de fenômenos bioquímicos (por exemplo, na descrição da velocidade de reações

enzimáticas). Apesar de ser um modelo simples contendo apenas dois parâmetros,

conforme mostra a Equação (4.1), a não-linearidade do modelo pode prejudicar a

aproximação elíptica da região de confiança.

1

2

xyx

θθ

⋅=+

(4.1)

No modelo definido na equação acima, y é a variável independente, x é a

variável dependente e θ1 e θ2 são os dois parâmetros do modelo. Os dados

experimentais utilizados para estimação dos parâmetros estão descritos na Tabela 29.

Estes dados referem-se a uma reação enzimática onde a variável independente

corresponde à concentração de substrato em ppm e a variável dependente a velocidade

da reação, medida em contagens / min2.

Tabela 29: Dados experimentais usados na estimação dos parâmetros do modelo de Michaelis-Menten (BATES e WATTS, 1988).

x y x y 0,02 76 0,22 159 0,02 47 0,22 152 0,06 97 0,56 191 0,06 107 0,56 201 0,11 123 1,10 207 0,11 139 1,10 200

A função objetivo utilizada para a estimação foi a de mínimo quadrados e foi

utilizado o algoritmo Estima (NORONHA et al., 1993) para a estimação dos parâmetros

Page 97: Marcio Schwa Ab

90

e obtenção da matriz de covariância dos parâmetros e construção da região de confiança

elíptica. Nestes modelos simples e quando se dispõe de uma boa estimativa inicial dos

parâmetros, métodos de Newton são muito superiores aos demais e são usados aqui para

obtenção de um resultado preciso, o qual será usado para verificar se o método do

Enxame de Partículas alcançou o mínimo.

Os parâmetros estimados com o Estima foram 212,7 e 0,06412 para θ1 e θ2,

respectivamente, onde a função objetivo assume o valor mínimo de 1195,45. A matriz

de covariância dos parâmetros é

3

3 6

0, 4037 0,3682.100,3682.10 0,5736.10

− −

=

θV (4.2)

indicando uma correlação de 0,7651 entre os parâmetros. A região de confiança elíptica

foi construída de acordo com a Equação (2.36) para comparação com a região obtida

através da Equação (2.39), utilizando o método do Enxame de Partículas.

A estimação de parâmetros foi realizada com o método do Enxame de Partículas

na região de busca limitada por [0, 500] para o parâmetro θ1 e [0, 1] para o parâmetro

θ2. A busca consistiu na realização de 1000 iterações com 40 partículas, c1 e c2 iguais a

2,0 e w decrescendo linearmente de 1,2 a 0,8. Os parâmetros estimados e o valor

mínimo da função objetivo foram os mesmos alcançados com o Estima e a região de

confiança dos parâmetros pode ser descrita com um total de 1240 pontos avaliados pelo

método do Enxame de Partículas e cuja função objetivo é inferior ao limite de 2176,40,

calculado através da Equação (2.39). Na Figura 22 é apresentado o ajuste do modelo aos

dados experimentais, com os parâmetros estimados pelo Enxame de Partículas. Na

Figura 23 as regiões de confiança obtidas por aproximação elíptica, utilizando a matriz

de covariância fornecida pelo Estima, e através do Enxame de Partículas são

comparadas.

A Figura 23 mostra que as regiões de confiança obtidas pelos dois métodos são

muito próximas e estão praticamente sobrepostas. É possível observar que o contorno da

região obtida com os pontos avaliados pelo método do Enxame de Partículas é bem

definido, mostrando a boa qualidade das regiões obtidas desta forma. Este resultado

indica que apesar da não-linearidade do modelo, a aproximação linear que resulta na

região de confiança elíptica tem uma boa qualidade. Porém, isto pode ter ocorrido

devido à boa qualidade do ajuste obtido, pois se o ajuste é bom, a região de confiança

dos parâmetros é pequena. Quanto mais próximo do ponto mínimo, melhor é a

Page 98: Marcio Schwa Ab

91

aproximação quadrática da função objetivo e, conseqüentemente, a região de confiança

tem a forma próxima de uma elipse.

Figura 22: Ajuste do modelo de Michaelis-Menten aos dados experimentais.

Figura 23: Comparação entre a região de confiança elíptica e a obtida com o Enxame de Partículas na estimação dos parâmetros do modelo de Michaelis-Menten.

Page 99: Marcio Schwa Ab

92

Para verificar este fato, dados experimentais foram simulados usando o modelo

de Michaelis-Menten com os parâmetros θ1 e θ2 iguais a 100 e 0,5 e adicionando aos

valores calculados pelo modelo um erro aleatório com distribuição normal com média

zero e desvio padrão igual a 10, ou seja, simulando o erro experimental. Os dados

simulados são apresentados na Tabela 30.

Tabela 30: Dados simulados usados na estimação dos parâmetros do modelo de Michaelis-Menten.

x y 0,2 20,41 0,2 43,01 0,5 55,31 0,5 72,80 1,0 56,91 1,0 76,57 2,5 100,86 2,5 73,98 5,0 84,86 5,0 78,21

Os parâmetros do modelo foram ajustados aos dados experimentais com a

minimização da função de mínimos quadrados e os resultados obtidos com o Estima

foram os valores de 91,39 e 0,2972, respectivamente, para os parâmetros θ1 e θ2, onde o

valor da função objetivo é 1265,13. A matriz de covariância dos parâmetros foi igual a:

2

2 4

0, 4395 0,4728.100,4728.10 0,8487.10

− −

=

θV (4.3)

de onde é calculado o valor de 0,7742 para a correlação entre os parâmetros.

A minimização com o Enxame de Partículas foi realizada na região definida pelo

intervalo [0, 500] para o parâmetro θ1 e [0, 10] para o parâmetro θ2. A busca consistiu

na realização de 1000 iterações com 40 partículas, c1 e c2 iguais a 2,0 e w decrescendo

linearmente de 1,2 a 0,8. Os valores encontrados foram 91,57 e 0,3008,

respectivamente, para os parâmetros θ1 e θ2, com o valor mínimo da função objetivo

igual a 1265,30, sendo o limite máximo da função objetivo para a construção da região

de confiança igual a 2675,77, calculado pela Equação (2.39), de forma que o número de

pontos na região de confiança foi igual a 1854. Como pode ser observado, o valor

mínimo obtido é muito próximo do valor obtido pelo Estima, indicado que o método do

Enxame de Partículas alcançou o mínimo da função. Na Figura 24 é apresentado o

Page 100: Marcio Schwa Ab

93

ajuste do modelo aos dados simulados, com os parâmetros estimados pelo Enxame de

Partículas. Em seguida foram construídas as regiões de confiança, a elíptica e a obtida

com o método do Enxame de Partículas, as quais são apresentadas na Figura 25.

Figura 24: Ajuste do modelo de Michaelis-Menten aos dados simulados.

Figura 25: Comparação entre a região de confiança elíptica e a obtida com o Enxame de Partículas na estimação dos parâmetros do modelo de Michaelis-Menten com os

dados simulados.

Page 101: Marcio Schwa Ab

94

Diferentemente do que foi observado na Figura 23, onde as duas regiões de

confiança são próximas, na Figura 25 as duas regiões não estão tão próximas. É possível

observar que a região de confiança obtida pelo método do Enxame de Partículas não é

simétrica em relação ao ponto mínimo, indicando que é falha a hipótese de que a

distribuição de probabilidade dos erros paramétricos seja normal. Além disso, é possível

observar que a região elíptica contém valores negativos do parâmetro θ2, o que na

realidade não ocorre como mostra a região obtida pelo método do Enxame de Partículas.

Assim, a partir do resultado obtido pela aproximação elíptica da região de confiança a

conclusão obtida seria que o parâmetro θ2 não é importante para o modelo, pois pode

assumir o valor zero. Mas o que ocorre de fato, é que a aproximação elíptica é ruim e a

região de confiança verdadeira obtida pelo método do Enxame de Partículas não contém

valores negativos do parâmetro θ2, permitindo assim a análise correta dos resultados

obtidos.

4.1.2. Modelo Cinético de 1ª Ordem

Modelos cinéticos de 1ª ordem são modelos onde a velocidade de uma reação é

considerada proporcional à concentração de um determinado reagente. Assim, dada uma

reação da forma A → B, pode-se definir o consumo de A ao longo do tempo como:

AA

dC k Cdt=− ⋅ (4.4)

Sendo CA0 a concentração do reagente A no tempo inicial 0, a solução da

Equação (4.4) pode ser escrita como:

( )0 expA AC C k t= ⋅ − ⋅ (4.5)

Considerando a relação estequiométrica, onde 1 mol de A gera 1 mol de B, e que

no tempo inicial a concentração de B é nula, a Equação (4.5) pode ser reescrita de forma

a mostrar a concentração do produto B ao longo do tempo, como segue abaixo:

( )( )0 1 expB AC C k t= ⋅ − − ⋅ (4.6)

Para manter a nomenclatura que vem sendo usada, a Equação (4.6) pode ser

escrita da seguinte forma:

( )( )1 21 expy xθ θ= ⋅ − − ⋅ (4.7)

Page 102: Marcio Schwa Ab

95

onde y é a variável dependente e corresponde à concentração de B, x é a variável

independente e corresponde ao tempo de reação e θ1 e θ2 são os parâmetros do modelo

que serão ajustados. Apesar de não ser comum, considerar a concentração inicial de um

componente da reação como um parâmetro a ser estimado facilita o procedimento de

estimação. Além disso, em alguns casos a medida deste valor inicial, particularmente

em processos industriais, é difícil ou inviável (PRATA, 2005).

Assim, o modelo definido na Equação (4.7) foi ajustado aos dados experimentais

apresentados na Tabela 31, sendo minimizada a função de mínimos quadrados. Os

dados apresentados referem-se a medidas da demanda bioquímica de oxigênio em mg/L

(variável dependente) em função do tempo de amostragem em dias (variável

independente).

Tabela 31: Dados experimentais usados na estimação dos parâmetros do modelo cinético de 1ª ordem (BATES e WATTS, 1988).

x y 1,0 8,3 2,0 10,3 3,0 19,0 4,0 16,0 5,0 15,6 7,0 19,8

O resultado obtido com o Estima foi um valor mínimo de 25,99 da função

objetivo, com os parâmetros, θ1 e θ2 iguais a 19,14 e 0,5311, respectivamente. A matriz

de covariância dos parâmetros obtida foi:

1

1 2

0,9588 0,6653.100,6653.10 0,6347.10

− −

− = − θV (4.8)

sendo a correlação entre os parâmetros igual a -0,8528.

A minimização realizada pelo método do Enxame de Partículas foi limitada na

região de busca definida por [0, 100] para ambos parâmetros, realizando um total de

1000 iterações com 40 partículas, c1 e c2 iguais a 2,0 e w decrescendo linearmente de

1,2 a 0,8. Os resultados obtidos foram os valores de 19,15 e 0,5273 para os parâmetros

θ1 e θ2, respectivamente, e o valor de 26,00 para o mínimo da função objetivo,

mostrando que o mínimo foi alcançado. O número de pontos na região de confiança que

o método do Enxame de Partículas obteve foi de 2722, sendo o limite máximo da

função objetivo de 116,3 para os pontos na região de confiança, de acordo com a

Page 103: Marcio Schwa Ab

96

Equação (2.39). Na Figura 26 o ajuste do modelo aos dados experimentais é

apresentado e na Figura 27 as regiões de confiança dos parâmetros são comparadas.

Figura 26: Ajuste do modelo cinético de 1ª ordem aos dados experimentais.

Figura 27: Comparação entre a região de confiança elíptica e a obtida com o Enxame de Partículas na estimação dos parâmetros do modelo cinético de 1ª ordem.

Page 104: Marcio Schwa Ab

97

Como pode ser observado na Figura 27, a região de confiança obtida com o

método do Enxame de Partículas apresenta uma forma não convexa, mostrando que a

aproximação elíptica da região de confiança leva a resultados totalmente equivocados.

A aproximação elíptica mostra que o parâmetro θ2 tem um erro pequeno e que a maior

variação ocorre com relação ao parâmetro θ1. Entretanto, o que ocorre na realidade é

que ambos os parâmetros têm uma incerteza muito grande, já que fixando um dos

parâmetros é possível ajustar o outro de forma a conseguir um ajuste com uma

qualidade dentro do nível de confiança estabelecido. Além disso, devido à forma

matemática do modelo, quando θ2 tende ao infinito, o modelo fica reduzido a y = θ1, ou

seja, uma reta paralela ao eixo x na Figura 26. Como os dados experimentais contêm

uma incerteza muito grande, o ajuste de uma reta paralela ao eixo x não é

significativamente pior que o melhor ajuste do modelo e, assim, a região de confiança

dos parâmetros é aberta, conforme discutido por BATES e WATTS (1988).

Neste problema, também fica claro que quando a incerteza é grande, já que são

poucos pontos experimentais e o espalhamento destes é relativamente grande, a

aproximação elíptica é inviável. Assim, para evidenciar este fato, novos dados foram

simulados usando mais pontos experimentais e com a inserção de um erro experimental

pequeno. Os dados simulados, apresentados na Tabela 32, foram gerados através da

Equação (4.7) com θ1 e θ2 iguais a 20,0 e 0,5, respectivamente, e adicionando ao valor

calculado um erro aleatório com distribuição normal com média zero e desvio padrão

igual a 1,0.

Tabela 32: Dados simulados usados na estimação dos parâmetros do modelo cinético de 1ª ordem.

x y 1 8,79 2 11,21 3 16,27 4 16,27 5 18,63 6 20,51 7 20,32 8 18,91 9 19,13 10 19,74

Com os dados da Tabela 32 os parâmetros do modelo cinético de 1ª ordem

foram estimados com o Estima, sendo que os valores encontrados para θ1 e θ2 foram

Page 105: Marcio Schwa Ab

98

20,06 e 0,4961, respectivamente, onde a função objetivo assume o valor mínimo de

8,637. A matriz de covariância dos parâmetros obtida foi:

1

1 2

0, 2991 0,2129.100,2129.10 0,2735.10

− −

− = − θV (4.9)

de onde é calculada uma correlação de -0,7443 entre os parâmetros.

A minimização realizada com o método do Enxame de Partículas foi limitada no

intervalo de [0, 50] para θ1 e [0, 1] para θ2. Foram utilizadas 1000 iterações com 40

partículas, c1 e c2 iguais a 2,0 e w decrescendo linearmente de 1,2 a 0,8. O valor mínimo

da função objetivo foi de 8,637 e os valores para os parâmetros θ1 e θ2 foram 20,06 e

0,4960, mostrando que o método alcançou resultado idêntico ao obtido pelo Estima.

Além disso, número de pontos na região de confiança que o método do Enxame de

Partículas avaliou e que têm o valor da função objetivo inferior ao limite de 18,26

(calculado pela Equação (2.39)) foi de 3333. Na Figura 28 é apresentado o modelo

ajustado aos dados simulados e na Figura 29 a região de confiança dos parâmetros

obtida pelo Enxame de Partículas é comparada com a região elíptica.

Figura 28: Ajuste do modelo cinético de 1ª ordem aos dados simulados.

Page 106: Marcio Schwa Ab

99

Figura 29: Comparação entre a região de confiança elíptica e a obtida com o Enxame de Partículas na estimação dos parâmetros do modelo cinético de 1ª ordem com os

dados simulados.

Como pode ser observado na Figura 29, a qualidade da aproximação elíptica é

muito superior à que foi observada na Figura 27, mostrando que quanto maior a

precisão dos parâmetros obtidos, que é função do número de dados experimentais e do

erro experimental, melhor é a aproximação elíptica da região de confiança. Por outro

lado, a região de confiança obtida pelo método do Enxame de Partículas apresenta uma

boa qualidade, já que uma linha que delimita a região de confiança pode ser facilmente

traçada contornando os pontos da Figura 29.

4.1.3. Modelo de Exponencial Dupla

Modelos de Exponencial Dupla, isto é, constituídos por uma soma de duas

exponenciais, como mostra a Equação (4.10), são muito comuns em sistemas onde

ocorrem transformações sucessivas, como reações em série.

( ) ( )1 3 2 4exp expy x xθ θ θ θ= ⋅ − ⋅ − ⋅ − ⋅ (4.10)

Este modelo apresenta 4 parâmetros para serem ajustados, os quais apresentam

uma alta correlação entre si, sendo esta uma característica de modelos que apresentam a

função exponencial como na Equação (4.10).

Page 107: Marcio Schwa Ab

100

Os dados foram simulados com os parâmetros θ1, θ2, θ3 e θ4 respectivamente

iguais a 200, 100, 0,1 e 0,8 e foi adicionado ao valor calculado pelo modelo um erro

aleatório com distribuição de probabilidade normal com média zero e desvio padrão 5.

Na Tabela 33 os dados simulados e usados para a estimação dos parâmetros do modelo

de Exponencial Dupla são apresentados.

Tabela 33: Dados simulados usados na estimação dos parâmetros do modelo de Exponencial Dupla.

x y 0,0 97,56 0,5 117,10 1,0 138,93 1,5 148,13 2,0 150,00 2,5 139,25 3,0 138,91 3,5 136,89 4,0 134,71 4,5 124,53 5,0 125,35 5,5 115,33 6,0 105,78 6,5 105,95 7,0 95,10 7,5 98,82 8,0 82,67 8,5 93,83 9,0 79,79 9,5 75,92 10,0 74,52

A minimização realizada com o Estima alcançou um valor mínimo de 318,51

para a função objetivo, sendo que os parâmetros estimados foram :

θ1 = 203,70; θ2 =108,12; θ3 = 0,1020; e θ4 = 0,8440.

A minimização realizada pelo Enxame de Partículas foi realizada na região de

busca delimitada por [0, 1000] para os parâmetros θ1 e θ2, e por [0, 5] para os

parâmetros θ3 e θ4. Foram realizadas 5000 iterações, com 20 partículas, c1 e c2 iguais a

1,5 e w igual a 0,7. Porém, como este modelo apresenta um número maior de

parâmetros e uma certa dificuldade na minimização devido à alta correlação entre os

parâmetros, a minimização foi realizada de forma que sempre que a convergência das

partículas foi detectada, a busca era totalmente reiniciada. Esta forma de busca dá maior

Page 108: Marcio Schwa Ab

101

garantia de que o mínimo encontrado seja o mínimo global e aumenta a exploração da

região de busca em torno do ponto mínimo.

O valor mínimo da função objetivo e dos parâmetros estimados encontrado pelo

método do Enxame de Partículas foi o mesmo valor encontrado acima pelo Estima. Na

Figura 30 o ajuste do modelo aos dados simulados é apresentado.

Figura 30: Ajuste do modelo de Exponencial Dupla aos dados simulados.

A região de confiança dos parâmetros estimados foi construída com os dados

avaliados pelo método do Enxame de Partículas, com um total de 53426 pontos cujas

funções objetivo têm valor inferior ao limite de 540,69. As regiões são apresentadas

para cada par de parâmetros, conforme ilustrado na Figura 31. Além disso, o ponto onde

a função tem o valor mínimo aparece destacado para cada par de parâmetros.

Neste caso, a região elíptica não é apresentada, já que a determinação desta

região quando o número de parâmetros é maior que dois começa a apresentar

dificuldades adicionais relacionadas às projeções das elipses aos planos de

representação. Já com a utilização do método do Enxame de Partículas, a definição das

regiões de confiança não apresenta maiores dificuldades, já que os pontos utilizados

para isso são os mesmos usados pelo método de otimização para encontrar o mínimo da

função objetivo, sem a necessidade de um esforço computacional adicional.

Page 109: Marcio Schwa Ab

102

Figura 31: Regiões de confiança dos parâmetros do modelo de Exponencial Dupla obtidas pelo método do Enxame de Partículas.

Como pode ser observado na Figura 31, as regiões de confiança são curvas e não

simétricas com relação ao ponto mínimo, já que este se encontra deslocado para uma

das extremidades da região, mostrando que a aproximação elíptica, neste caso, além de

ser de difícil obtenção não representa a forma real das regiões de confiança.

Os três problemas de estimação de parâmetros em modelos não-lineares

resolvidos até aqui evidenciaram que a aproximação elíptica das regiões de confiança

não consegue representar com qualidade as regiões de confiança verdadeiras. Por outro

lado, as regiões verdadeiras podem ser bem definidas com a utilização dos pontos

avaliados pelo método do Enxame de Partículas, os quais são selecionados pelo critério

descrito pela Equação (2.39) que define a região de verossimilhança.

Page 110: Marcio Schwa Ab

103

Além disso, foi possível observar que o aumento da incerteza paramétrica

prejudica ainda mais a qualidade da aproximação elíptica. A incerteza paramétrica é

tanto maior quanto maiores são os erros experimentais e quanto menor é o número de

dados experimentais. Este comportamento ocorre devido ao fato da qualidade da

aproximação quadrática da função objetivo melhorar à medida que se aproxima do

ponto mínimo, o que acontece quando a região de confiança é pequena; ou seja, quando

os parâmetros são precisos.

Por fim, foi possível mostrar a facilidade e a boa qualidade das regiões de

confiança obtidas com o método do Enxame de Partículas. Além disso, a construção da

região de confiança não implica em custo computacional adicional, já que as regiões de

confiança são construídas a partir dos pontos que o método do Enxame de Partículas já

havia utilizado para a minimização da função objetivo, sendo necessária apenas a

seleção dos pontos que fazem parte da região de confiança, de acordo com algum

critério, como o definido pela Equação (2.39).

4.2. Exemplo: Cinética da Reação de Polimerização de Etileno

O estudo dos processos de polimerização tem uma importância muito grande,

devido à grande aplicabilidade que os materiais poliméricos têm atualmente. O processo

de polimerização do etileno com complexos de níquel contendo isotiocianato, em

particular, pode ser suscetível a mudanças da atividade do catalisador devido a

alterações em sua estrutura, que alteram as características do sítio ativo onde a reação de

polimerização ocorre e que, assim, promovem alterações de atividade do catalisador,

conforme apontado por CROSSETTI et al. (2004).

O catalisador presente no reator sofre inicialmente uma ativação, ou seja, passa

de uma espécie inativa para uma espécie ativa, capaz de promover as reações de

polimerização. Em seguida, esta espécie ativa pode sofrer alterações gerando,

gradativamente, novas espécies ativas, cujas atividades são diferentes da primeira. O

número de alterações que o catalisador pode sofrer não é conhecido a priori, a menos

que avaliações mais rigorosas sejam feitas, como uma análise química detalhada do

catalisador ao longo da reação de polimerização. Entretanto, admite-se que devem

ocorrer somente duas ou três transformações sucessivas, podendo então estar presentes

duas ou três diferentes espécies ativas do catalisador, já que a primeira espécie é inativa.

O conjunto de transformações sucessivas que ocorrem com o catalisador pode

ser visto como um conjunto de reações sucessivas descritas abaixo:

Page 111: Marcio Schwa Ab

104

0

1

1

0 1

1 2

1

1

n

m

k

k

kn n

km m

C C

C C

C C

C C−

+

(4.11)

onde Cn corresponde à concentração presente da n-ésima forma do catalisador, kn é a

constante de velocidade de transformação da forma n na forma posterior n+1 e m é o

número de transformações que ocorrem e também corresponde ao número de espécies

ativas.

Assim, a Equação (4.11) denota um conjunto de transformações sucessivas a

partir de uma espécie inicial C0 até uma espécie final Cm. Considerando que as

velocidades das transformações ocorrem na forma de taxas de primeira ordem em

relação a cada uma das espécies e que a atividade de cada espécie do catalisador é

proporcional a sua concentração, é possível escrever um conjunto de equações de

balanço cuja solução (descrita com maiores detalhes no Apêndice) é dada pela seguinte

equação:

( )0

i

m mk tn

n ii n i

Rp t Kp A e−

= =

=

∑ ∑ (4.12)

onde Kpn é a constante de polimerização da n-ésima espécie do catalisador. Os demais

coeficientes são calculados de forma recursiva, como descrito abaixo:

( )

00

1 1

1

0

1

; 0.. -1 ; 0n n ni i

n i

nn nn i

i

AkA A i n n

k k

A A

− −

=

=

= = >−

=−∑

(4.13)

É importante observar que o valor de km é nulo, já que m é a última espécie do

catalisador e não sofre nenhuma transformação. Além disso, o valor de Kp0 também é

nulo, já que a espécie inicial de catalisador é inativa.

Assim, o número de constantes cinéticas que devem ser estimadas a partir de

dados experimentais é igual a 2m, sendo m constantes cinéticas de transformação de

Page 112: Marcio Schwa Ab

105

espécie do catalisador (k0, k1, ..., km-1) e m constantes cinéticas de polimerização (Kp1,

Kp2, ..., Kpm).

Um caso especial que pode ser considerado é a hipótese de que as

transformações sucessivas das espécies têm uma mesma constante de velocidade, ou

seja, k0 = k1 = ... = km-1 = k. Neste caso, a solução do sistema de equações de balanço

(descrita com maiores detalhes no Apêndice) tem a seguinte forma:

( ) ( )1

0 !

nmkt

n m mn

kRp t Kp Kp e Kpn

−−

=

= − + ∑ (4.14)

Neste modelo, onde as constantes de transformação das espécies do catalisador

são iguais, o número de parâmetros que são estimados é igual a m+1, sendo 1 constante

cinética de transformação de espécie do catalisador e m constantes cinéticas de

polimerização (Kp1, Kp2, ..., Kpm), lembrando que Kp0 é igual a zero.

Os dados experimentais da polimerização de etileno usados para estimação dos

parâmetros foram obtidos por SILVA (2004). A reação foi realizada em um reator semi-

batelada, onde o etileno era alimentado ao reator de modo a manter a pressão do sistema

constante, já que à medida que o etileno é consumido pela reação, a pressão tende a cair.

Assim, a taxa de polimerização pode ser acompanhada pelo consumo de etileno durante

o decorrer da reação. Estes dados serão usados para a estimação dos 2m parâmetros

cinéticos do modelo desenvolvido acima. Maiores detalhes do sistema reacional podem

ser obtidos em SILVA (2004).

Foram usados os dados de consumo de etileno da corrida NCS07, que

corresponde a uma reação realizada na temperatura de 0 ºC, pressão de 3 atm e o

catalisador com a razão de Al/Ni igual a 460. Foram medidos 402 valores de consumo

de etileno em aproximadamente 22 minutos de reação.

Inicialmente foram consideradas duas transformações sucessivas do catalisador

(ou seja, m = 2) e foram ajustados 4 parâmetros, de acordo com o modelo definido pelas

Equações (4.12) e (4.13). A função objetivo usada para a estimação dos parâmetros foi a

função de mínimos quadrados, cuja minimização foi realizada pelo método do Enxame

de Partículas. Foram realizadas diversas buscas com diferentes configurações para

aumentar a confiança de que o mínimo global foi encontrado. Os resultados

apresentados a seguir foram obtidos, em particular, com uma busca consistindo de um

total de 15000 iterações e 25 partículas. Os valores de c1 e c2 usados foram iguais a 1,5 e

w foi mantido constante e igual a 0,8. Além disso, sempre que a convergência das

Page 113: Marcio Schwa Ab

106

partículas era detectada, as posições e as velocidades das partículas eram reiniciadas e

os melhores pontos encontrados, tanto o individual quanto o global, eram “apagados da

memória das partículas”; ou seja, uma nova busca era iniciada. A convergência das

partículas foi definida pela diferença média das funções objetivo de todo o conjunto de

partículas com relação ao melhor ponto encontrado até o momento, conforme a Equação

(3.17). Neste caso, sempre que a diferença média for menor que 50, a busca é reiniciada,

de forma a serem realizadas diversas novas buscas dentro do total de 15000 iterações,

aumentando o caráter exploratório da busca.

É importante observar que, por se tratar de um problema real de estimação de

parâmetros, cuja solução não é conhecida a priori, a despeito da eficiência do método

que já foi comprovada nos exemplos anteriores, um número grande de avaliações é

necessário para garantir que a solução encontrada é realmente a melhor possível. Por

isso, o número de iterações deve ser grande.

A faixa de busca dos parâmetros foi definida após algumas buscas inicias, onde

se procurou encontrar faixas mais restritas para a busca pelo mínimo, já que

inicialmente a única informação disponível é que os parâmetros são positivos e finitos.

Além disso, uma minimização em uma região de busca mais restrita possibilita uma

melhor descrição da região de confiança dos parâmetros. As faixas definidas para os

parâmetros correspondem aos seguintes intervalos: para k0 [0, 100]; para k1 [0, 50]; para

Kp1 [0, 1000]; e para Kp2 [0, 10].

O valor mínimo da função objetivo encontrado pelo Enxame de Partículas ao

final da busca foi de 1001,21, sendo os parâmetros estimados apresentados na Tabela

34.

Tabela 34: Parâmetros estimados considerando 2 transformações no catalisador.

Parâmetro Valor estimado k0 6,3805 k1 6,3807

Kp1 150,54 Kp2 6,2781

É interessante observar que as constantes de transformação do catalisador são

praticamente as mesmas, indicando que o modelo definido pela Equação (4.14) deve ser

mais apropriado para este caso. Além disso, vê-se que a primeira espécie ativa tem uma

atividade de polimerização muito superior à da segunda espécie ativa.

Page 114: Marcio Schwa Ab

107

Na Figura 32 é apresentado o ajuste do modelo aos dados experimentais, bem

como a taxa correspondente a cada espécie do catalisador (sendo Rpi a taxa de

polimerização da espécie i). Já a Figura 33 apresenta a variação da fração de cada

espécie de catalisador ao longo do tempo.

Figura 32: Ajuste do modelo (linha preta) com 2 transformações no catalisador aos dados experimentais (pontos) e taxa correspondente a cada espécie do catalisador

(Rp1: linha amarela; Rp2: linha verde).

Figura 33: Fração de cada espécie de catalisador ao longo do tempo de reação. C0: linha vermelha; C1: linha amarela; C2: linha verde.

Page 115: Marcio Schwa Ab

108

As Figuras 32 e 33 mostram que as transformações no catalisador são

extremamente rápidas e em apenas 1 minuto praticamente todo o catalisador já está

transformado na espécie final. A partir daí, a atividade total de polimerização é

constante e igual à atividade da última espécie do catalisador.

A região de confiança dos parâmetros foi determinada com um nível de

confiança de 95 % e foi usado o critério da razão de verossimilhança (descrito pela

Equação (2.39)) em conjunto com os pontos avaliados pelo método do Enxame de

Partículas durante a minimização da função objetivo. A região de confiança dos

parâmetros é apresentada na Figura 34, onde a região de confiança é mostrada para cada

par de parâmetros, para permitir a avaliação do efeito da correlação entre cada par dos

parâmetros.

A constante de transformação da espécie inicial e a constante de polimerização

da primeira espécie ativa, k0 e Kp1 respectivamente, são muito pouco significativas, já

que os intervalos de confiança destes parâmetros são muito amplos. A constante de

transformação da primeira espécie ativa na ultima espécie do catalisador, k1, também

apresenta um intervalo de confiança grande, mas não tanto como os anteriores. Este fato

ocorre devido às rápidas transformações que ocorrem logo no início da reação, já que

todo o catalisador é convertido na sua espécie final em apenas 1 minuto. Assim, não faz

muita diferença se a constante de transformação da espécie inicial é igual a 6 ou 100, já

que de uma forma ou outra esta espécie é convertida rapidamente. O mesmo ocorre com

a primeira espécie ativa que, além de ter uma constante de transformação cujo limite do

intervalo de confiança vai de 4,0 a 36,0, tem uma constante de polimerização que pode

assumir valores desde 100 a 1000, pois mais uma vez esta espécie é consumida

rapidamente e seu efeito no ajuste é pequeno.

Já a constante de polimerização da última espécie ativa, Kp2, tem um intervalo

de confiança muito estreito. Isto ocorre porque, após o minuto inicial, este é o único

parâmetro que tem algum efeito no ajuste do modelo e pequenas mudanças no seu valor

têm um grande efeito no valor da função objetivo, já que pode afastar o modelo de um

número muito grande de dados experimentais.

É interessante observar que as regiões de confiança não são simétricas em

relação ao ponto mínimo. Além disso, no caso dos parâmetros k0, k1 e Kp1 as regiões de

confiança apresentam formas nada convencionais, deixando claro que a aproximação

elíptica da região de confiança não deve ser utilizada.

Page 116: Marcio Schwa Ab

109

Figura 34: Regiões de confiança dos parâmetros do modelo de polimerização com 2 transformações no catalisador.

Como foi comentado acima, os valores das constantes de transformação do

catalisador são praticamente iguais e, por isso, foi usado o modelo definido pela

Equação (2.14) que, no caso de duas transformações sucessivas, tem 3 parâmetros. É

importante observar que quando as constantes de transformação do catalisador são

iguais, o modelo definido pela Equação (2.12) apresenta uma descontinuidade, mas os

valores calculados tendem no limite aos valores calculados pelo modelo de k’s iguais;

ou seja, se os k’s são iguais, ambos modelos fornecem os mesmos resultados. A

diferença fica no menor número de parâmetros que o modelo de k’s iguais apresenta.

Apesar de isso não alterar o ajuste do modelo,essa modificação altera as regiões de

confiança dos parâmetros de forma muito significativa, conforme mostra a Figura 35.

Page 117: Marcio Schwa Ab

110

A estimação dos parâmetros do modelo com k’s iguais foi realizada com as

mesmas condições descritas anteriormente para o modelo com constantes de

transformação diferentes. O valor mínimo da função objetivo e os parâmetros estimados

foram os mesmos obtidos antes e apresentados na Tabela 34, lembrando que para este

modelo k0 e k1 são definidos por k, já que têm o mesmo valor.

Figura 35: Regiões de confiança dos parâmetros do modelo de polimerização com 2 transformações no catalisador com constantes de transformação iguais.

As regiões de confiança apresentadas na Figura 35 mostram que os parâmetros

estimados com o segundo modelo apresentam uma maior precisão, já que, a despeito da

curvatura da região de confiança para os parâmetros k e Kp1, as regiões são bem

definidas, sendo regiões fechadas e menores que as obtidas para o caso anterior.

Page 118: Marcio Schwa Ab

111

O parâmetro Kp1 ainda apresenta um intervalo de confiança relativamente

grande, mas é mais restrito que o obtido para o caso anterior. Os parâmetros k e Kp2

apresentam intervalos de confiança bem mais restritos, indicando uma boa precisão

nestes valores.

Apesar da melhor qualidade dos parâmetros estimados, o ajuste do modelo aos

dados experimentais não apresenta uma boa qualidade, já que após o primeiro minuto de

reação a taxa de polimerização permanece constante, comportamento que não é

observado nos dados experimentais.

Desta forma, foi utilizado um modelo com 3 transformações sucessivas no

catalisador, considerando que as constantes de transformação podem ser diferentes, de

acordo com a Equação (4.12). Neste caso, são 6 parâmetros que devem ser estimados e

o método do Enxame de Partículas foi usado para a minimização da função objetivo

com um total de 10000 iterações, com 50 partículas, c1 e c2 iguais a 1,5 e w constante e

igual a 0,8. Como no caso anterior, a busca era reiniciada sempre que a convergência

das partículas era detectada.

Durante as buscas preliminares para a definição da região de busca dos

parâmetros, foi observado que este problema apresenta dois mínimos onde a função

objetivo assume o valor de 106,39, ou seja, dois mínimos globais. Na Tabela 35 os

valores dos parâmetros em cada mínimo são apresentados, bem como os limites de

busca para cada parâmetro.

Tabela 35: Limite de busca e parâmetros estimados para cada mínimo e com o modelo considerando 3 transformações no catalisador.

Parâmetro Limite de busca Mínimo A Mínimo B k0 [0, 100] 10,31 2,313 k1 [0, 50] 2,313 10,30 k2 [0, 1] 0,232 0,232

Kp1 [0, 1000] 55,98 249,6 Kp2 [0, 1] 0,000 0,000 Kp3 [0, 10] 7,587 7,587

A Tabela 35 mostra que de um mínimo para outro, apenas os valores de k0 e k1

são trocados entre si e o valor de Kp1 muda, sendo que os demais valores são iguais nos

dois mínimos. O que ocorre é que apesar do ajuste da taxa global de polimerização ser

igual para os dois mínimos (Figura 36), as quantidades presentes das espécies inicial e

da segunda espécie do catalisador são diferentes (Figura 37) e a atividade da primeira

espécie ativa é diferente.

Page 119: Marcio Schwa Ab

112

Figura 36: Ajuste do modelo (linha preta) com 3 transformações no catalisador aos dados experimentais (pontos) e taxa correspondente a cada espécie do catalisador

(Rp1: linha amarela; Rp3: linha azul).

Figura 37: Fração de cada espécie de catalisador ao longo do tempo de reação para o mínimo A (linha tracejada) e mínimo B (linha pontilhada).

C0: linha vermelha; C1: linha amarela; C2: linha verde; C3: linha azul.

A Figura 36 mostra o ajuste do modelo aos dados experimentais e as taxas das

espécies 1 e 3, já que a espécie 2 tem atividade nula. Apesar da constante de

Page 120: Marcio Schwa Ab

113

polimerização da espécie 1 mudar de um mínimo para outro, a sua concentração

também muda, como mostra a Figura 37, de forma que a atividade desta espécie é a

mesma nos dois casos. Na Figura 37 também é possível observar que as concentrações

das espécies 2 e 3 são as mesmas para os dois mínimos. Já as espécies 0 e 1 têm valores

diferentes para cada mínimo.

Para o ponto de mínimo A, a espécie inicial é rapidamente convertida na espécie

1 e a conversão desta espécie na espécie 2 é mais lenta, aumentando rapidamente a

concentração da espécie 1. Para manter a taxa de polimerização em um nível adequado,

a constante de polimerização da espécie 1 é menor que no caso do mínimo B.

Para o ponto de mínimo B a constante de transformação da espécie 0 na espécie

1 é mais lenta que a constante de transformação da espécie 1 na espécie 2, mantendo

baixa a concentração da espécie 1. Neste caso, a constante de polimerização da espécie

1 é maior para compensar a baixa concentração e manter a taxa global em um nível

adequado para o bom ajuste do modelo.

Na Tabela 35 são apresentados os intervalos obtidos para cada mínimo na

estimação dos parâmetros. No mínimo A o parâmetro k0 alcança o limite máximo da sua

faixa de busca, enquanto no mínimo B é o parâmetro Kp1 que alcança o limite máximo

da sua faixa de busca. Além disso, o parâmetro Kp2 atinge o limite mínimo da sua faixa

de busca nos dois mínimos e é este limite que é o valor ótimo deste parâmetro. Os

intervalos de confiança dos parâmetros k2, Kp2 e Kp3 são aproximadamente os mesmos

nos dois mínimos. As regiões de confiança obtidas dos parâmetros estimados são

apresentadas para cada par de parâmetros na Figura 38.

Tabela 35: Intervalos de confiança para cada mínimo do ajuste do modelo com 3 transformações no catalisador.

Parâmetro Mínimo A Mínimo B k0 6,186 100,0 2,137 2,525 k1 2,150 2,515 6,467 35,87 k2 0,2165 0.2484 0,2157 0,2453

Kp1 49,39 75,58 139,6 1000 Kp2 0,000 0,1128 0,0000 0,1228 Kp3 7,444 7,746 7,467 7,741

Page 121: Marcio Schwa Ab

114

Figura 38a: Regiões de confiança dos parâmetros do modelo de polimerização com 3 transformações no catalisador.

Page 122: Marcio Schwa Ab

115

Figura 38b: Regiões de confiança dos parâmetros do modelo de polimerização com 3 transformações no catalisador.

Page 123: Marcio Schwa Ab

116

Figura 38c: Regiões de confiança dos parâmetros do modelo de polimerização com 3 transformações no catalisador.

Como pode ser observado na Figura 38, as regiões de confiança que envolvem

os parâmetros k0, k1 e Kp1 apresentam regiões de confiança desconectadas, devido ao

fato da função objetivo possuir dois mínimos globais. É interessante observar a alta

correlação entre os parâmetros k1 e Kp1 em torno do mínimo B. Isto ocorre porque k1 é a

constante de velocidade com que a espécie 1 é convertida. Se a espécie 1 desaparece

com maior velocidade, a constante de taxa de polimerização desta espécie, Kp1, tende a

aumentar para compensar a menor concentração desta espécie. Um fato idêntico, mas

não tão pronunciado, ocorre entre os parâmetros k0 e Kp1 para o mínimo A. As demais

regiões, bem como as discutidas acima, não são simétricas com relação aos valores

estimados o que, em conjunto com as outras características citadas acima, tornam a

aproximação elíptica destas regiões totalmente inapropriada.

Page 124: Marcio Schwa Ab

117

Como esperado, o ajuste do modelo considerando três transformações sucessivas

do catalisador apresentou uma qualidade bem superior ao ajuste obtido com o modelo

com somente duas transformações sucessivas. Por outro lado, a existência de dois

mínimos locais e a conseqüente formação de regiões de confiança desconectadas reforça

a necessidade da utilização de métodos heurísticos de otimização para estimação de

parâmetros.

É importante informar que a aplicação dos algoritmos baseados no método de

Gauss-Newton, como o Estima, não obteve sucesso, isto é, a busca não convergiu. Isto

ocorreu devido à presença de parâmetros pouco significativos e de parâmetros cujo

valor ótimo se encontram no limite da região de busca. Além disso, mesmo quando os

parâmetros eram significativos, uma estimativa inicial muito próxima do ponto mínimo

era necessária. Já no caso onde a função objetivo apresentou dois mínimos globais, o

que foi prontamente apontado pelos resultados obtidos com o Enxame de Partículas, a

utilização de um método de Gauss-Newton apresentaria dificuldades adicionais, já que

diversas estimativas iniciais deveriam ser fornecidas e a chance de serem observados os

dois mínimos é menor do que quando é utilizado o método do Enxame de Partículas,

cuja busca é baseada em um grande conjunto de pontos, o que aumenta a exploração da

região de busca, possibilitando a observação da região de confiança verdadeira, mesmo

quando esta é composta por regiões disjuntas.

4.3. Exemplo: Cinética da Reação de Produção de Gás de Síntese

A produção do gás de síntese é uma importante etapa para a conversão do gás

natural em combustíveis mais práticos e outros produtos químicos de interesse. As

reações envolvidas na produção do gás de síntese fazem parte de um sistema reacional

muito complexo (LARENTIS, 2001) e incluem, além de outras que não serão

consideradas, as seguintes reações:

4 2 2CH +H O CO+3H reforma com vapor (1)

4 2 2CH +CO 2CO+2H reforma com CO2 (2)

14 2 22CH + O CO+2H oxidação parcial (3)

4 2 2 2CH +2O CO +2H O oxidação total(4)

2 2 2CO+H O CO +H deslocamento água-gás(5)

Page 125: Marcio Schwa Ab

118

Considerando as cinco reações reversíveis acima, mesmo considerando um

modelo cinético simplificado, onde as taxas das reações são consideradas de primeira

ordem em relação a cada reagente e as constantes cinéticas modeladas de acordo com a

equação de Arrhenius, o número de parâmetros é igual a 20.

Assim, a utilização de métodos do tipo Gauss-Newton torna-se difícil, pois além

das dificuldades associadas ao grande número de parâmetros, a correlação entre estes

parâmetros e a presença de parâmetros não significativos torna o procedimento de

estimação um verdadeiro desafio. Um exemplo disso foram as dificuldades encontradas

no trabalho de LARENTIS et al. (2001), onde para ser possível a estimação, diversas

hipóteses tiveram de ser realizadas para a simplificação do sistema de equações,

permitindo a estimação de um conjunto menor de parâmetros. Estas hipóteses nem

sempre podem ser feitas de forma segura, o que pode dar uma desconfiança nos

resultados obtidos. Esta dificuldade motiva ainda mais a utilização dos algoritmos

heurísticos em problemas de estimação de parâmetros.

Para a estimação dos parâmetros cinéticos deste sistema reacional serão

utilizados os dados experimentais obtidos por LARENTIS (2000). Estes dados foram

obtidos na reação de reforma com dióxido de carbono combinada com a oxidação

parcial do gás natural em um reator tubular de bancada contendo catalisador 1,12% Pt/γ-

Al2O3. A alimentação consistia de uma mistura de gás natural (79% CH4, 17% C2H6 e

4% C3H8), ar seco e comprimido (20% O2, 79% N2 e 1% Ar) e dióxido de carbono em

diferentes quantidades. A análise dos gases foi feita por cromatografia gasosa

(Chrompack CP 9001), usando uma coluna recheada (HAYESEP D) e detector de

condutividade térmica. Maiores detalhes podem ser encontrados em LARENTIS (2000).

O reator foi modelado como um reator tubular de fluxo empistonado e foram

usadas as equações de balanço de massa individual para cada espécie i

catii cat i

dC dvv C m Rdx dx

+ = (4.15)

onde v é a vazão volumétrica (m3/s), Ci é a concentração da espécie i (gmol/m3), mcat é a

massa de catalisador (g), Ricat é taxa de reação do componente i em relação à massa de

catalisador (gmol/gcat/s) e x é o comprimento adimensional do reator. Também foi usada

a equação do balanço de massa total

1

NCcat

cat ii

dv RT m Rdx P =

= ∑ (4.16)

Page 126: Marcio Schwa Ab

119

onde R é a constante dos gases ideais (J/gmol/K), T é a temperatura (K) e P é a pressão

(Pa), de forma a considerar a variação na vazão volumétrica devido à mudança do

número de moles nas reações. As condições iniciais das equações diferenciais são as

concentrações medidas na entrada do reator e a taxa de alimentação volumétrica,

também medida na entrada do reator. O número de espécies envolvidas na reação é

igual a 7: CH4, O2, CO, H2, CO2, H2O e N2, sendo o nitrogênio um composto inerte

neste sistema. As concentrações de cada espécie foram calculadas com a equação dos

gases ideais, como uma função das frações molares (yi) obtidas por cromatografia

gasosa, como segue:

ii

y PCRT

= (4.17)

O sistema de equações diferenciais é integrado com o auxílio do procedimento

numérico Dassl (PETZOLD, 1989).

As taxas de reação foram descritas com base na massa de catalisador e são

definidas como:

cati i catR R ρ= (4.18)

onde ρcat é a densidade do catalisador e Ri é a taxa em termos volumétricos (gmol/m3/s).

A taxa de reação individual de uma espécie i em particular consiste na soma das

taxas de cada reação em que a espécie i participa ponderada pelo coeficiente

estequiométrico, como a equação abaixo:

,i i j jj

R q r= ⋅∑ (4.19)

onde i indica o componente, j indica a reação, q indica o coeficiente estequiométrico do

componente i na reação j e o termo rj representa a taxa de cada reação (gmol/m3/s).

Inicialmente será considerado um modelo simples, já citado acima, onde as taxas

das cinco reações apresentadas acima são consideradas de primeira ordem em relação a

cada componente, como mostram as equações a seguir:

4 2 21 1 1CH H O r H COr k C C k C C= − (4.20a)

4 2 22 2 2CH CO r H COr k C C k C C= − (4.20b)

4 2 23 3 3CH O r H COr k C C k C C= − (4.20c)

4 2 2 24 4 4CH O r H O COr k C C k C C= − (4.20d)

2 2 25 5 5CO H O r H COr k C C k C C= − (4.20e)

Page 127: Marcio Schwa Ab

120

onde kj são as constantes cinéticas (m3/gmol/s) calculadas pela equação de Arrhenius:

( )( )exp 1j j j refk A B T T= − + − (4.21)

onde Tref é uma temperatura de referência (K) e Aj e Bj são os parâmetros que serão

estimados. Esta forma da equação de Arrhenius apresenta uma menor correlação entre

os parâmetros e faz com que a faixa de busca dos parâmetros Aj e Bj seja mais limitada

(KITTRELL, 1970). As relações dos parâmetros Aj e Bj com a energia de ativação e

fator pré-exponencial são dadas pelas seguintes equações:

j j refE B RT= (4.22)

( )expoj j jk A B= − + (4.23)

Como foi discutido acima, o número de parâmetros que devem ser estimados é

igual a 20 e a faixa de busca é definida pelo intervalo [0, 50] para todos os parâmetros.

Esta faixa corresponde a uma faixa em termos de energia de ativação de [0, 508]

kJ/gmol e em termos de do fator pré-exponencial de [0, 5,2.1021] m3/gmol/s.

Este modelo simplificado será usado inicialmente para uma comparação entre os

métodos heurísticos de otimização considerados neste trabalho. A função objetivo é a

função de mínimos quadrados ponderados (Equação 2.08), sendo a variância de cada

medida, estimada a partir de réplicas, igual a 8.10-3 para todas as concentrações

(LARENTIS, 2000).

Foi definido um valor máximo de 100 mil avaliações da função objetivo e foram

usadas três configurações diferentes para cada método.

O método de Monte Carlo foi usado com o valor da taxa de redução calculado

pela Equação (3.11) de tal forma que, ao final das iterações, a faixa de busca

correspondia a 1 % da faixa inicial; ou seja, a faixa final tinha um tamanho de 0,5, já

que a faixa inicial tem o tamanho de 50. Na Tabela 36 são apresentadas as

configurações de busca e o valor mínimo da função objetivo encontrado em cada busca.

A busca MC2 alcançou o valor da função objetivo de 9220,4 e foi o menor valor

encontrado pelo método de Monte Carlo, enquanto a busca MC5 obteve o valor de

14016,9, a busca com o maior valor final. Já o valor médio da função objetivo foi igual

a 11640,8.

Page 128: Marcio Schwa Ab

121

Tabela 36: Configurações de busca e valor mínimo da função objetivo encontrado pelo método de Monte Carlo.

Nome Npt Niter TR [%] Fmin MC1 25 4000 0,115 10516,0 MC2 25 4000 0,115 9220,4 MC3 25 4000 0,115 10561,4 MC4 50 2000 0,230 10085,4 MC5 50 2000 0,230 14016,9 MC6 50 2000 0,230 12869,2 MC7 100 1000 0,460 11031,2 MC8 100 1000 0,460 13764,0 MC9 100 1000 0,460 12702,8

O método do Algoritmo Genético foi usado com uma probabilidade de

cruzamento igual a 70 %, valor utilizado na minimização das funções teste no Capítulo

3, uma probabilidade de mutação de 40 %, valor que obteve os melhores resultados na

minimização das funções teste, e uma perturbação relativa nos parâmetros realizada pela

mutação igual a 1 %. Foram realizadas 9 minimizações com três configurações

diferentes, as quais em conjunto com o valor mínimo encontrado são apresentadas na

Tabela 37.

Tabela 37: Configurações de busca e valor mínimo da função objetivo encontrado pelo método do Algoritmo Genético.

Nome Npt Niter Fmin AG1 25 4000 9630,6 AG2 25 4000 10321,7 AG3 25 4000 11672,2 AG4 50 2000 11120,0 AG5 50 2000 8438,6 AG6 50 2000 11107,9 AG7 100 1000 9827,5 AG8 100 1000 10264,5 AG9 100 1000 10269,3

O menor valor da função objetivo encontrado pelo método do Algoritmo

Genético foi de 8438,6 na busca AG5 e é inferior ao menor valor encontrado pelo

método de Monte Carlo. A corrida AG3 teve como valor mínimo 11672,2 e foi o pior

resultado obtido pelo Algoritmo Genético. O valor médio da nove minimizações foi

igual a 10294,7.

Para a utilização do método do Recozimento Simulado a temperatura inicial e

final foram definidas “a priori” como auxilio da Equação (2.21), sendo o valor do fator

Page 129: Marcio Schwa Ab

122

de redução de temperatura calculado de tal forma a satisfazer estas temperaturas. A

temperatura inicial foi definida como 20000, de forma que com uma probabilidade de

95%, uma transição para um valor maior da função objetivo de até 1026,0 fosse aceita.

Já a temperatura final foi definida como 1, sendo que com uma probabilidade de apenas

10%, uma transição para um valor maior da função objetivo só será aceita se o aumento

for menor que 2,3, de forma a busca para algum mínimo próximo de onde a ponto de

busca se encontra. O valor das perturbações para a realização das transições é de 5 % da

faixa total. Na Tabela 38 são apresentadas as configurações de cada busca e o valor

mínimo encontrado para a função objetivo em cada busca.

Tabela 38: Configurações de busca e valor mínimo da função objetivo encontrado pelo método do Recozimento Simulado.

Nome Npt Niter α Fmin RS1 250 400 0,976 6622,5 RS2 250 400 0,976 7265,3 RS3 250 400 0,976 9489,0 RS4 500 200 0,952 9398,3 RS5 500 200 0,952 9642,6 RS6 500 200 0,952 11181,7 RS7 1000 100 0,906 7368,8 RS8 1000 100 0,906 6280,8 RS9 1000 100 0,906 10145,3

Os resultados obtidos pelo método de Recozimento simulado são superiores ao

encontrados pelo método de Monte Carlo e pelo método do Algoritmo Genético. O

menor valor foi igual a 6280,8 na busca RS7, e o maior valor foi 11181,7 na busca RS6.

O valor médio encontrado foi de 8599,3, um valor menor que os valores médios obtidos

pelo método de Monte Carlo e do Algoritmo Genético.

A boa qualidade dos resultados obtidos nas minimizações realizadas pelo

Recozimento Simulado, superiores aos obtidos pelo Algoritmo Genético e Monte Carlo,

deve ter sido fruto da definição adequada das temperaturas inicial e final, de forma a

garantir uma boa exploração no início da busca, seguida de uma convergência para o

ponto mínimo na medida que a temperatura diminui. Antes da realização das corridas

acima, pode ser observado que o valor da temperatura inicial determinada

automaticamente, com a utilização da Equação (3.15), era muito diferente em cada

busca realizada, já que este valor depende do ponto inicial, determinado aleatoriamente,

e das transições iniciais realizadas em torno deste ponto para a determinação do valor de

∆F* usado na Equação (3.15). Desta forma, parece ser mais adequada uma análise

Page 130: Marcio Schwa Ab

123

preliminar do problema, particularmente da grandeza nos valores da função objetivo,

para a definição adequada dos valores das temperaturas inicial e final.

O método do Enxame de Partículas foi usado com os parâmetros c1 e c2 iguais a

1,5 e o valor de w decrescendo ao longo das iterações de 1,0 até 0,7. O valor máximo de

avaliações da função objetivo foi fixado em 100 mil avaliações e o número de partículas

e de iterações ajustados de forma a satisfazer este valor. Na Tabela 39 são apresentadas

as condições de cada busca e o valor mínimo alcançado em cada busca.

Tabela 39: Configurações de busca e valor mínimo da função objetivo encontrado pelo método do Enxame de Partículas.

Nome Npt Niter Fmin EP1 25 4000 6253,6 EP2 25 4000 7163,6 EP3 25 4000 7087,7 EP4 50 2000 7821,1 EP5 50 2000 6295,8 EP6 50 2000 6275,8 EP7 100 1000 6328,9 EP8 100 1000 7302,9 EP9 100 1000 9316,0

O menor valor encontrado pelo método do Enxame de Partículas foi igual a

6253,6 na busca EP1, e este valor é o menor encontrado por todos os algoritmos usados.

Além disso, as corridas EP5, EP6, e EP7 alcançaram valores muito próximos ao menor

valor obtido. A busca com o pior resultado foi a EP9 e a função objetivo foi igual a

9316,0. O valor médio alcançado pelo método do Enxame de Partículas foi de 7039,9, o

menor valor médio dentre os métodos usados. Este resultado mostra mais uma vez que o

método do Enxame de Partículas tem uma eficiência superior aos demais métodos

avaliados neste trabalho.

Como o método do Enxame de Partículas obteve os melhores resultados, este foi

utilizado para a definição dos intervalos de confiança dos parâmetros, através da seleção

dos pontos avaliados cuja função objetivo satisfez a razão de verossimilhança, definida

pela Equação (2.39), e em seguida, observando os valores mínimo e máximo de cada

parâmetro neste conjunto de dados selecionados. Para isso, o método do Enxame de

Partículas foi utilizado com 50 partículas, 5000 iterações, w com o valor constante de

0,7 e c1 e c2 iguais a 1, 5. Para aumentar a exploração da região de busca, sempre que

após 500 iterações sucessivas sem que o melhor ponto encontrado tivesse uma queda

Page 131: Marcio Schwa Ab

124

maior que 1 a velocidade de cada partícula era reiniciada, de forma a aumentar a

exploração da região de busca e evitar uma convergência prematura.

Ao final desta busca o melhor valor encontrado foi de 6227,8 e na Tabela 40

podem ser observados os valores dos parâmetros estimados e os seus respectivos limites

inferiores e superiores e o valor das constantes cinéticas na temperatura de referência.

Tabela 40: Resultados da estimação dos parâmetros do modelo simplificado da reação.

Parâmetro Limite inferior Valor estimado Limite superior ki (em 1222 K) A1 10,80 42,93 50,00 B1 0,00 13,36 50,00 1,44.10-13

A1r 13,74 15,04 50,00 B1r 0,00 12,95 22,86 1,24.10-01

A2 10,60 49,57 50,00 B2 0,00 24,92 50,00 1,99.10-11

A2r 14,82 42,03 50,00 B2r 0,00 0,00 3,51 5,59.10-19

A3 4,19 5,04 6,62 B3 0,94 3,15 8,65 1,51.10-01

A3r 13,70 48,18 50,00 B3r 0,00 0,99 30,82 3,20.10-21

A4 5,59 6,06 9,92 B4 0,00 0,14 3,72 2,69.10-03

A4r 6,99 8,19 9,17 B4r 2,34 7,44 12,03 4,69.10-01

A5 13,59 16,98 50,00 B5 29,75 49,99 50,00 2,16.10+14

A5r 11,01 11,45 11,84 B5r 2,89 4,76 7,35 1,25.10-03

A Tabela 40 mostra que são poucos os parâmetros que têm um intervalo de

confiança pequeno. A maioria dos parâmetros relacionados com as reações de reforma

com vapor (reação 1), reforma com CO2 (reação 2) e deslocamento água-gás (reação 5)

são pouco significativos.

Este resultado pode ter ocorrido devido à má qualidade do modelo utilizado,

onde as equações das taxas de reação foram consideradas de primeira ordem (Equações

4.20), sendo assim um modelo simplificado e que não leva em consideração a adsorção

dos compostos presentes. Assim, com o objetivo de melhorar o ajuste do modelo, foi

considerada a adsorção dos compostos presentes, com exceção do oxigênio. A adsorção

do oxigênio não foi considerada, uma vez que como em nenhuma condição foi

observado oxigênio na saída do reator, indicando que seu consumo é rápido, acreditou-

Page 132: Marcio Schwa Ab

125

se que não seria necessário considerar a sua adsorção. As equações das taxas são então

definidas como:

4 2 21 11 1

CH H O r H CO

i ii

k C C k C Cr

K C−

=+∑

(4.24a)

4 2 22 22 1

CH CO r H CO

i ii

k C C k C Cr

K C−

=+∑

(4.24b)

4 2 23 33 1

CH O r H CO

i ii

k C C k C Cr

K C−

=+∑

(4.24c)

4 2 2 24 44 1

CH O r H O CO

i ii

k C C k C Cr

K C−

=+∑

(4.24d)

2 2 25 55 1

CO H O r H CO

i ii

k C C k C Cr

K C−

=+∑

(4.24e)

onde i = CH4, CO, H2, CO2 e H2O.

A dependência das constantes de adsorção com a temperatura é definida como

( )( )exp 1i i i refK A B T T= − + − (4.25)

onde i define um componente específico e Ai e Bi são os parâmetros a serem estimados.

Na expressão da constante de adsorção o parâmetro Bi está relacionado com a entalpia

de adsorção que pode ser tanto positiva quanto negativa, dependendo se a adsorção é

exotérmica ou endotérmica. Assim, para sua estimação, os 20 parâmetros relacionados

com as constantes cinéticas tiveram seus limites de busca definidos pelo intervalo [0,

50] e os 10 parâmetros relacionados com as constantes de adsorção tiveram o limite de

busca definido pelo intervalo [-30, 30], tanto para os Ai como para os Bi.

As estimações realizadas, mais uma vez, somente com o método do Enxame de

Partículas foram realizadas com um total de 20000 iterações, 50 partículas, w fixo e

igual a 0,7 e c1 e c2 iguais a 1,5. Mais uma vez, sempre que após 500 iterações

consecutivas a queda do melhor valor encontrado da função objetivo for menor que 1 a

velocidade das partículas eram reiniciadas, impedindo a convergência prematura e

aumentando a exploração da região de busca.

O valor mínimo da função objetivo encontrado foi igual a 5036,96, mostrando

uma queda significativa com relação ao modelo simplificado onde as adsorções não

Page 133: Marcio Schwa Ab

126

foram consideradas. Na Tabela 41 são apresentados os parâmetros estimados, os limites

de confiança e o valor das constantes na temperatura de referência.

Tabela 41: Resultados da estimação dos parâmetros do modelo considerando as adsorções dos componentes da reação.

Parâmetro Limite inferior Valor estimado Limite superior ki (em 1222 K) A1 0,39 2,68 3,38 B1 15,95 17,18 22,67 1,98.1006

A1r 5,41 7,21 8,22 B1r 9,86 16,27 26,77 8,56.1003

A2 1,71 3,63 5,37 B2 25,88 30,15 45,53 3,29.1011

A2r 7,93 8,63 16,06 B2r 0,00 6,44 11,02 1,12.10-01

A3 0,00 0,00 0,87 B3 9,00 11,39 19,29 8,82.1004

A3r 24,87 45,28 50,00 B3r 0,00 0,06 50,00 2,29.10-20

A4 0,68 0,86 4,53 B4 0,00 0,00 0,07 4,24.10-01

A4r 0,00 0,00 2,33 B4r 12,24 14,93 47,48 3,05.1006

A5 3,18 4,41 5,03 B5 9,08 10,47 14,77 4,30.1002

A5r 2,17 3,20 3,54 B5r 13,03 14,17 17,79 5,80.1004

ACH4 -8,58 -7,98 -7,82 BCH4 14,80 16,04 18,67 2,68.1010

ACO -5,98 29,15 30,00 BCO -30,00 -29,97 30,00 2,11.10-26

AH2 -5,23 29,64 30,00 BH2 -30,00 21,04 30,00 1,84.10-04

ACO2 -6,91 29,44 30,00 BCO2 -21,56 -6,25 30,00 3,17.10-16

AH2O -6,43 27,26 30,00 BH2O -30,00 -29,98 30,00 1,38.10-25

Como pode ser visto na Tabela 41, com exceção da constante de adsorção do

metano, as demais constantes de adsorções têm valores muito pequenos, como pode ser

visto pelo valor destes constantes na temperatura de referência. Além disso, os

parâmetros relacionadas a estas constantes de adsorção têm intervalos de confiança

muito grandes, mostrando que são parâmetros pouco significativos.

Assim, uma nova estimação foi realizada com um modelo onde somente a

adsorção do metano foi considerada, conforme as equações abaixo.

Page 134: Marcio Schwa Ab

127

4 2 2

4 4

1 11 1

CH H O r H CO

CH CH

k C C k C Cr

K C−

=+

(4.26a)

4 2 2

4 4

2 22 1

CH CO r H CO

CH CH

k C C k C Cr

K C−

=+

(4.26b)

4 2 2

4 4

3 33 1

CH O r H CO

CH CH

k C C k C Cr

K C−

=+

(4.26c)

4 2 2 2

4 4

4 44 1

CH O r H O CO

CH CH

k C C k C Cr

K C−

=+

(4.26d)

2 2 2

4 4

5 55 1

CO H O r H CO

CH CH

k C C k C Cr

K C−

=+

(4.26e)

Na Tabela 42 são apresentados os valores estimados dos parâmetros, os limites

inferior e superior de confiança e o valor das constantes na temperatura de referência.

Tabela 42: Resultados da estimação dos parâmetros do modelo considerando somente a adsorção do metano.

Parâmetro Limite inferior Valor estimado Limite superior ki (em 1222 K) A1 3,95 6,37 13,50 B1 0,00 8,22 12,29 6,36.10+00

A1r 9,23 11,16 19,66 B1r 0,00 10,38 50,00 4,56.10-01

A2 4,39 6,68 50,00 B2 0,00 22,30 36,37 6,05.10+06

A2r 11,07 22,03 50,00 B2r 0,00 0,01 50,00 2,73.10-10

A3 2,47 3,42 6,09 B3 0,00 0,70 7,59 6,58.10-02

A3r 10,14 10,72 50,00 B3r 0,00 0,00 50,00 2,21.10-05

A4 0,48 0,77 5,29 B4 0,00 0,00 1,19 4,63.10-01

A4r 0,00 0,00 5,78 B4r 6,54 15,54 23,63 5,62.10+06

A5 6,80 7,99 15,40 B5 0,00 0,47 2,94 5,39.10-04

A5r 5,85 6,76 9,83 B5r 2,10 4,72 6,96 1,308.10-01

ACH4 -5,57 -4,61 -0,94 BCH4 2,74 7,25 9,60 1,41.10+05

Com este modelo, onde somente a adsorção do metano foi considerada

(Equações 4.26), o valor mínimo da função objetivo alcançado foi 5033,50. Este valor é

Page 135: Marcio Schwa Ab

128

até um pouco inferior ao encontrado com o modelo onde todas as adsorções foram

consideradas, indicando que naquele caso a busca não alcançou um mínimo,

provavelmente devido ao grande número de parâmetros do modelo.

Ainda procurando encontrar uma forma mais simples do modelo, mas mantendo

a qualidade do ajuste, foi realizada uma simplificação relacionada ao fato de que em

todas as condições experimentais observadas a concentração de oxigênio na saída do

reator era nula, ou seja, as reações de oxidação são rápidas ou praticamente instantâneas

quando comparadas às demais reações. Assim, como feito por LARENTIS et al. (2001),

o consumo de oxigênio pela reação de oxidação parcial (Reação 3) é definido pelo

parâmetro λ e o consumo de oxigênio pela reação de oxidação total (Reação 4) é

definido pelo parâmetro ξ. As modificações causadas por estas reações instantâneas são

inseridas no problema através da mudança das concentrações da alimentação do reator,

como segue abaixo

2

1,22 O inCξ λ+ = (4.27a)

4 4,0 ,CH CH inC C ξ λ= − − (4.27b)

2 ,0 0OC = (4.27c)

,0 ,CO CO inC C λ= + (4.27d)

2 2,0 , 2H H inC C λ= + (4.27e)

2 2,0 ,CO CO inC C ξ= + (4.27f)

2 2,0 , 2H O H O inC C ξ= + (4.27g)

2 2,0 ,N N inC C= (4.27h)

Como a concentração de alimentação de oxigênio é conhecida, é apenas

necessário estimar um dos parâmetros λ ou ξ, já que o outro pode ser calculado através

da Equação (4.27a). Assim, o parâmetro estimado foi λ e o seu intervalo de busca foi de

2 ,0, 2 O inC . Aqui os parâmetros λ e ξ foram considerados independentes da

temperatura, diferente do trabalho de LARENTIS et al. (2001) onde estes parâmetros

foram considerados como funções da temperatura.

Com a simplificação, as equações do modelo foram reduzidas às equações a

seguir:

4 2 2

4 4

1 11 1

CH H O r H CO

CH CH

k C C k C Cr

K C−

=+

(4.28a)

Page 136: Marcio Schwa Ab

129

4 2 2

4 4

2 22 1

CH CO r H CO

CH CH

k C C k C Cr

K C−

=+

(4.28b)

2 2 2

4 4

5 55 1

CO H O r H CO

CH CH

k C C k C Cr

K C−

=+

(4.28e)

sendo apenas excluídas as equações de taxa relacionadas às reações de oxidação parcial

e total.

Desta forma, o número de parâmetros a serem estimados foi reduzido a 14,

sendo 12 parâmetros das constantes cinéticas, 2 parâmetros de adsorção do metano e 1

parâmetro relacionado às reações de oxidação.

O valor mínimo encontrado foi igual a 5664,83, indicando que considerar

instantâneas as reações de oxidação prejudicou a qualidade do ajuste. Na Tabela 43 são

apresentados os parâmetros estimados, os limites inferior e superior de confiança e o

valor das constantes na temperatura de referência.

Tabela 43: Resultados da estimação dos parâmetros do modelo considerando somente a adsorção do metano e considerando instantâneas as reações de oxidação.

Parâmetro Limite inferior Valor estimado Limite superior ki (em 1222 K) A1 2,06 2,42 6,76 B1 5,54 13,03 13,57 4,06.10+04

A1r 9,20 9,90 13,41 B1r 0,00 0,00 5,68 5,03.10-05

A2 5,65 11,43 50,00 B2 0,00 49,87 50,00 4,98.10+16

A2r 7,76 8,51 19,58 B2r 0,00 15,28 34,32 8,77.10+02

A5 4,57 5,08 9,53 B5 0,00 4,79 6,93 7,52.10-01

A5r 3,78 4,11 8,85 B5r 0,14 8,08 9,53 5,31.10+01

ACH4 -6,43 -6,23 -1,51 BCH4 -0,70 7,58 8,06 9,96.10+05

λ 0,00 0,00 0,00 ---

O parâmetro λ, correspondente à reação de oxidação parcial, teve seu valor

estimado igual a 0 e, mesmo os limites inferior e superior de confiança também iguais a

0, mostrando que todo o oxigênio é consumido na reação de oxidação total, de acordo

com o resultado obtido por LARENTIS et al. (2001).

Na Figura 39 os valores das concentrações preditas pelo modelo são comparados

aos valores experimentais, considerando os resultados obtidos com o modelo

Page 137: Marcio Schwa Ab

130

representado pelas Equações (4.26), onde somente a adsorção do metano foi

considerada e não foram feitas simplificações com relação às reações de oxidação,

sendo este o melhor resultado obtido.

Como pode ser observado na Figura 39, a concordância entre os valores preditos

e experimentais é de boa qualidade. Entretanto, apesar do modelo ter uma relação com

os modelos de Langmuir-Hinshelwood, este modelo não pode ser considerado um

modelo fenomenológico, já que as reações, a despeito da adsorção do metano, são

reações de primeira ordem em relação a cada componente, sendo definidas assim

unicamente para simplificar o modelo e diminuir o número de parâmetros a serem

estimados.

Contudo, neste trabalho foi mostrado que os algoritmos heurísticos, em

particular o Enxame de Partículas, são capazes de estimar parâmetros, mesmo em

problemas com muitos parâmetros. Desta forma, a obtenção de modelos baseados no

mecanismo das reações que ocorrem na superfície do catalisador, mesmo resultando em

modelos com muitos parâmetros, não deve apresentar dificuldades adicionais, já que o

próprio método de otimização, no caso o Enxame de Partículas, ajuda na determinação

de parâmetros pouco significativos, indicando onde modificações podem ser realizadas.

Neste exemplo, ficou claro que a adsorção dos componentes era desnecessária,

com exceção da adsorção do metano, tanto que o modelo que só considerou a adsorção

do metano alcançou-se um resultado similar ao modelo com as adsorções de todos os

componentes, mesmo com 8 parâmetros a menos.

Por outro lado, os resultados não indicavam que as reações de oxidação

poderiam ser simplificadas, já que os parâmetros ligados a estas reações sempre

apresentaram intervalos de confiança pequenos. Assim, após fazer a hipótese de que

estas reações poderiam ser consideradas instantâneas, a qualidade do ajuste do modelo

foi prejudicada, como já poderia ter sido esperado.

Enfim, este exemplo mostra que, além de possibilitar a estimação de parâmetros

em problemas de grande dimensão, os algoritmos heurísticos, em particular o Enxame

de Partículas, permitem uma avaliação do modelo e indicam possíveis simplificações,

ajudando no desenvolvimento dos modelos.

Page 138: Marcio Schwa Ab

131

Figura 39: Comparação dos valores preditos e observados das concentrações de cada componente. (a) metano; (b) monóxido de carbono; (c) hidrogênio; (d) dióxido de

carbono; (e) água; (f) nitrogênio.

(a)

(c) (d)

(e) (f)

(b)

Page 139: Marcio Schwa Ab

132

CAPÍTULO 5: CONCLUSÕES E SUGESTÕES

Neste trabalho foi avaliada a utilização de algoritmos heurísticos de otimização

em problemas de estimação de parâmetros. Foram considerados quatro algoritmos

heurísticos: Monte Carlo, Algoritmo Genético, Recozimento Simulado e Enxame de

Partículas. Inicialmente foram avaliados os efeitos dos parâmetros de busca de cada

método na minimização de funções teste e na construção da região de confiança dos

parâmetros de um modelo. Em seguida, o método do Enxame de Partículas foi aplicado

na determinação de regiões de confiança de parâmetros de modelos não-lineares. O

método do Enxame de Partículas também foi aplicado para a estimação de parâmetros

de dois problemas cinéticos: um da reação de polimerização de etileno e outro da reação

de produção de gás de síntese.

Na minimização das funções teste, o método de Monte Carlo obteve uma alta

taxa de sucessos, mas somente nas buscas com redução na região de busca. Apesar dos

bons resultados, a redução da região de busca pode trazer certos inconvenientes, como

excluir o ponto ótimo da região de busca e concentrar os pontos avaliados em uma

região muito próxima ao mínimo, prejudicando a representação da região de confiança

dos parâmetros.

O método do Algoritmo Genético alcançou um bom número de sucessos.

Entretanto, a definição dos parâmetros parece depender muito da forma da função

objetivo, como foi observado para as funções de Rosenbrock e função de Levy. Apesar

do número de iterações necessárias ser maior que no método de Monte Carlo, o

Algoritmo Genético não apresenta os inconvenientes da utilização da redução da região

de busca.

O método do Recozimento Simulado apresentou dificuldades na minimização

das funções teste, devido provavelmente ao cálculo da temperatura inicial, o qual é feito

com um certo grau de aleatoriedade e, dependendo dos valores calculados, pode

prejudicar muito a minimização por este método. Por outro lado, as regiões de confiança

obtidas, apesar de dependerem da escolha adequada dos parâmetros de busca,

apresentaram uma boa qualidade.

O método do Enxame de Partículas obteve os melhores resultados nas

minimizações das funções teste. Entretanto, a definição adequada dos parâmetros de

busca é fundamental para garantir o sucesso da minimização. Também foi observada

Page 140: Marcio Schwa Ab

133

uma forte correlação entre os valores de c1 e c2 e o valor de w, aparentemente

independente do problema que se está resolvendo. Com relação às regiões de confiança

dos parâmetros, as regiões obtidas com o método do Enxame de Partículas apresentaram

uma boa distribuição dos pontos, possibilitando uma boa descrição da região de

confiança. Outra vantagem do método do Enxame de Partículas é o controle da

exploração da região de busca através da definição do parâmetro w, que pode ser um

valor fixo ou um valor decrescente para promover uma busca com caráter global nas

primeiras iterações e forçar a convergência ao final da busca.

Na obtenção das regiões de confiança dos modelos não-lineares foi possível

mostrar casos onde a aproximação elíptica não representa a verdadeira região de busca e

reforçar a necessidade da utilização dos métodos heurísticos de otimização, permitindo

uma análise estatística rigorosa dos parâmetros estimados. Além disso, a determinação

da região de confiança a partir dos pontos avaliados pelo método do Enxame (ou por

outro método heurístico) não representa custo computacional adicional e transforma

uma desvantagem do método em uma importante vantagem sobre os métodos

tradicionais de estimação.

Na estimação dos parâmetros cinéticos da reação de polimerização de etileno

foram obtidas regiões de confiança dos parâmetros, cuja forma em nada se assemelha

com elipses. Além disso, este problema de estimação apresenta uma alta correlação

entre os parâmetros, parâmetros pouco significativos e a surpreendente existência de

dois mínimos globais (no caso de três transformações no catalisador), e mesmo assim, o

método do Enxame de Partículas foi capaz de obter a solução e de descrever com boa

qualidade as regiões de confiança dos parâmetros.

Na estimação dos parâmetros cinéticos da reação de produção de gás de síntese,

os algoritmos heurísticos foram comparados novamente, com um modelo contendo 20

parâmetros cinéticos. Mais uma vez, o método do Enxame de Partículas obteve os

melhores resultados. Entretanto, o método do Recozimento Simulado obteve resultados

superiores aos obtidos pelo Algoritmo Genético e pelo Monte Carlo, o que não foi

observado na minimização das funções teste. Este melhor desempenho do método do

Recozimento Simulado ocorreu devido à uma definição adequada da temperatura inicial

e final da busca, sendo o valor da taxa de redução de temperatura α calculado de forma

a satisfazer os valores definidos para a temperatura inicial e final e para o número de

temperaturas utilizadas. Pode ser observado que o cálculo da temperatura inicial estava

levando a valore muito diferentes em cada inicialização da minimização, devido ao alto

Page 141: Marcio Schwa Ab

134

caráter aleatório do procedimento de cálculo da temperatura inicial. Assim, a definição

adequada destes valores permitiu um melhor desempenho do método do Recozimento

Simulado.

Mas como o método do Enxame de Partículas ainda obteve os melhores

resultados, este foi usado em novas estimações, onde foram consideradas as adsorções

dos componentes, com exceção do oxigênio, totalizando 30 parâmetros. Os resultados

obtidos permitem, através da análise dos intervalos de confiança de cada parâmetro,

identificar as reações que têm maior importância no sistema e permite a realização de

simplificações no modelo de forma prática e segura.

Desta forma, foi possível mostra as vantagens da utilização dos métodos

heurísticos de otimização, em particular do método do Enxame de Partículas, na

resolução de problemas de estimação de parâmetros, garantindo uma busca de caráter

global, sem necessidade de estimativas iniciais e do cálculo de derivadas e permitindo

uma análise estatística rigorosa dos parâmetros estimados, através da construção de uma

região de confiança sem a necessidade de aproximações elípticas da região de

confiança.

Uma sugestão para os usuários destes métodos é a necessidade da realização de

diversas minimizações para garantir que o menor valor encontrado realmente é o

mínimo, já que estes métodos têm um alto caráter aleatório e, assim, a solução

encontrada pode nem sempre ser o ponto ótimo. Com a realização de diversas buscas,

alterando os parâmetros de busca para alterar o caráter exploratório da busca, é possível

garantir a definição do ponto ótimo.

Algumas sugestões para trabalhos futuros são:

a) a necessidade de um número maior de avaliações de cada um dos métodos, de

forma a identificar de forma precisa o efeito dos parâmetros de busca nos

resultados obtidos por cada um dos métodos heurísticos;

b) procurar formas automáticas para a definição apropriada dos parâmetros de

busca, sem a necessidade do usuário do método precisar definir estes parâmetros

de busca, já que para a definição apropriada destes valores é preciso um bom

conhecimento do método que se está usando;

c) utilizar minimizações paralelas, com ou sem troca de informação, onde cada

busca pode ter uma característica diferente, explorando mais a região de busca

ou convergindo mais rapidamente, evitando a necessidade da realização de

diversas buscas individuais;

Page 142: Marcio Schwa Ab

135

d) utilizar técnicas de programação paralela para aumentar a velocidade dos

cálculos realizados e diminuir o tempo necessário para a estimação de

parâmetros com os métodos heurísticos, sem sombra de dúvida a maior

dificuldade para a utilização rotineira destas técnicas;

e) construção de algoritmos híbridos, onde a busca inicia com um método

heurístico e na medida que se aproxima da solução, um método do tipo de

Newton é utilizado para acelerar a convergência e garantir que o ponto

encontrado é um mínimo;

Por fim, a utilização dos métodos heurísticos na otimização e, em particular, na

estimação de parâmetros tem se apresentado como uma excelente alternativa, pois além

dos resultados promissores obtidos por estes algoritmos, são métodos extremamente

simples e fáceis de serem implementados. São métodos robustos e conseguem obter

boas soluções mesmo em problemas de difícil solução. Assim, a utilização destes

métodos deve ser considerada como uma alternativa viável e eficiente, pois apesar do

custo computacional para a avaliação de um grande número de pontos, estes podem ser

usados tanto na definição da região de confiança dos parâmetros, quanto para a análise

de sensitividade em problemas de otimização, aumentando assim a qualidade dos

resultados obtidos.

Page 143: Marcio Schwa Ab

136

REFERÊNCIAS

AHON, V.R., TAVARES, F.W., CASTIER, M., 2000, “A Comparison of Simulated

Annealing Algorithms in the Scheduling of Multiproduct Serial Batch Plants”,

Brazilian Journal of Chemical Engineering, v. 17, pp. 199-209.

ANDERSON, T.F., ABRAMS. D.S., GRENS II, E.A., 1978, “Evaluation of Parameters

for Nonlinear Thermodynamic Models”, AIChE Journal, v. 24, pp. 20-29.

BALLAND, L., ESTEL, L., COSMAO, J.M., MOUHAB, N., 2000, “A Genetic

Algorithm with Decimal Coding for the Estimation of Kinetic and Energetic

Parameters”, Chemometrics and Intelligent Laboratory Systems, v. 50, pp. 121-

135.

BARD, Y., 1970, “Comparison of Gradient Methods for the Solution of Nonlinear

Parameter Estimation Problems”, SIAM Journal on Numerical Analysis, v. 7, pp.

157-186.

BARD, Y., 1974, Nonlinear Parameter Estimation. New York, Academic Press.

BATES, D.M., WATTS, D.G., 1981, “Parameter Transformations for Improved

Approximate Confidence Regions in Nonlinear Least Squares”, Annals of

Statistics, v. 9, pp. 1152-1167.

BATES, D.M., WATTS, D.G., 1988, Nonlinear Regression Analysis and Its

Applications. New York, John Wiley & Sons.

BECK, J.V., ARNOLD, K.J., 1977, Parameter Estimation in Engineering and Science.

New York, John Wiley & Sons.

BISCAIA JR., E.C., SCHWAAB, M., PINTO, J.C., 2004, “Um Novo Enfoque do

Método do Enxame de Partículas”. Em: Nanobio: Workshop em Nanotecnologia e

Computação Inspirada em Biologia, Rio de Janeiro, RJ, Brasil, 23-24 Março.

BOX, G.E.P., DRAPER, N.R., 1965, “The Bayesian Estimation of Common Parameters

from Several Responses”, Biometrika, v. 52, pp. 355-365.

BRITT, H.I., LUECKE, R.H., 1973, “The Estimation of Parameters in Nonlinear,

Implicit Models”, Technometrics, v. 15, pp. 233-247.

BROOKS, S.H., 1958, “A Discussion of Random Methods for Seeking Maxima”,

Operations Research, v. 6, pp. 244-251.

Page 144: Marcio Schwa Ab

137

CLERC, M., 1999, “The Swarm and the Queen: Towards a Deterministic and Adaptive

Particle Swarm Optimization”. In: Proc. International Conference on

Evolutionary Computation, Washington, USA, pp. 1951-1957.

CLERC, M., KENNEDY, J., 2002, “The Particle Swarm - Explosion, Stability, and

Convergence in a Multidimensional Complex Space”, IEEE Transactions on

Evolutionary Computation, v. 6, pp. 58-73.

COSTA, A.L.H., SILVA, F.P.T., PESSOA, F.L.P., 2000, “Parameter Estimation of

Thermodynamic Models for High-Pressure Systems Employing a Stochastic

Method of Global Optimization”, Brazilian Journal of Chemical Engineering, v.

17, pp. 349-353.

COSTA, E.F., LAGE, P.L. C., BISCAIA JR., E.C., 2003, “On the Numerical Solution

and Optimization of Styrene Polymerization in Tubular Reactors”, Computers and

Chemical Engineering, v. 27, pp. 1591-1604.

CROSSETTI, G.L., DIAS, M.L., QUEIROZ, B.T., SILVA, L.P., ZIGLIO, C.M.,

BOMFIM, J.A.S., FILGUEIRAS, C.A.L., 2004, “Ethylene Polymerization with

Imine and Phosphine Nickel Complexes Containing Isothiocyanate”, Applied

Organometallic Chemistry, v. 18, pp. 331-336.

DAS, H., CUMMINGS, P.T., LEVAN, M.D., 1990, “Scheduling of Serial Multiproduct

Batch Processes via Simulated Annealing”, Computers and Chemical

Engineering, v. 14, pp. 1351-1362.

DIXON, L.C.W., SZEGO, G.P., 1975, Towards Global Optimization, Amsterdam,

North-Holland Publishing Company.

DONALDSON, J.R., SCHNABEL, R.B., 1987, “Computational Experience with

Confidence Regions and Confidence Intervals for Nonlinear Least Squares”,

Technometrics, v.29, pp. 67-82.

DOVI, V.G., REVERBERI, A.P., MAGA, L., 1993, “Optimal Design of Sequential

Experiments for Error-in-Variables Models”, Computers and Chemical

Engineering, v. 17, pp. 111-115.

DRAPER, N.R., SMITH, H., 1988, Applied Regression Analysis, 3a ed. New York, John

Wiley & Sons.

EBERHART, R., KENNEDY, J., 1995, “A New Optimizer Using Particle Swarm

Theory”. In: Proc. Sixth Symposium on Micro Machine and Human Science,

Nagoya, Japan, pp. 39-43.

Page 145: Marcio Schwa Ab

138

EBERHART, R., SHI, Y., 2000, “Comparing Inertia Weights and Constrictions Factors

in Particle Swarm Optimization”. In: Proc. Congress on Evolutionary

Computation, San Diego, USA, pp. 84-88.

EFTAXIAS, A., FONT., J., FORTUNY, A., FABREGAT, A., STÜBER, F., 2002,

“Nonlinear Kinetic Parameter Estimation Using Simulated Annealing”,

Computers and Chemical Engineering, v. 26, pp. 1725-1733.

FLETCHER, R e POWELL, M.J.F., 1963 “A Rapidly Convergent Descent Method for

Minimization”, Computer Journal, v. 6, pp. 163-168.

FROMENT, G.F., MEZAKI, R., 1970, “Sequential Discrimination and Estimation

Procedures for Rate Modeling in Heterogeneous Catalysis”, Chemical

Engineering Science, v. 25, pp. 293-301.

GALLANT, A.R., 1987, Nonlinear Statistical Models. New York, John Wiley & Sons.

GOLDBERG, D.E., 1989, Genetic algorithms in search, optimization and machine

learning. 1 ed. Boston, Addison Wesley Longman, Inc.

HALPERIN, M., 1963, “Confidence Interval Estimation in Non-Linear Regression”,

Journal of the Royal Statistical Society B, v. 25, pp. 330-333.

HARTLEY, H.O., 1964, “Exact Confidence Regions for the Parameters in Non-Linear

Regression laws”, Biometrika, v. 51, pp. 347-353.

HIBBERT, D.B., 1993a, “Genetic Algorithms in Chemistry”, Chemometrics and

Intelligent Laboratory Systems, v. 19, pp. 277-293.

HIBBERT, D.B., 1993b, “A Hybrid Genetic Algorithm for the Estimation of Kinetic

Parameters”, Chemometrics and Intelligent Laboratory Systems, v. 19, pp. 319-

329.

HIMMELBLAU, D.M., 1970, Process Analysis by Statistical Methods. New York, John

Wiley & Sons.

HIMMELBLAU, D.M., 1972, Applied Nonlinear Programming. New York, McGraw-

Hill.

HOOKE, R., JEEVES, T.A., 1961, “Direct Search Solution of Numerical and Statistical

Problems”, Journal of ACM., v. 8, pp. 212-229.

KENNEDY, J., EBERHART, R., 1995, “Particle Swarm Optimization”. In: Proc. IEEE

International Conference on Neural Networks, Perth, Australia, pp. 1942-1948.

KIM, I.W., EDGAR, T.F., BELL, N.H., 1991b, “Parameter Estimation for a Laboratory

Water-Gas-Shift Reactor using a Nonlinear Error-in-Variables Method”,

Computers and Chemical Engineering, v. 15, pp. 361-367.

Page 146: Marcio Schwa Ab

139

KIM, I.W., LIEBMAN, M.J., EDGAR, T.F., 1991a, “A Sequential Error-in-Variables

Method for Nonlinear Dynamic Systems”, Computers and Chemical Engineering,

v. 15, pp. 663-670.

KIRKPATRICK, S., GELLAT JR., C.D., VECCHI, M.P., 1983, “Optimization by

Simulated Annealing”, Science, v. 220, pp. 671-680.

KITTRELL, J.R., 1970, “Mathematical modeling of chemical reactions”, Advances in

Chemical Engineering, v. 8, pp. 97-183.

KLEPPER, O., HENDRIX, E.M.T., 1994a, “A Comparison of Algorithms for Global

Characterization of Confidence Regions for Nonlinear Models”, Environmental

Toxicology and. Chemistry, v. 13, pp. 1887-1899.

KLEPPER, O., HENDRIX, E.M.T., 1994b, “A Method for Robust calibration of

Ecological Models under Different Types of Uncertainty”, Ecological Modelling,

v. 74, pp. 161-182.

KLEPPER, O., BEDAUX, J.J.M., 1997a, “A Robust Method for Nonlinear Parameter

Estimation Illustrated on a Toxicological Model”, Nonlinear Analysis, v. 30, pp.

1677-1686.

KLEPPER, O., BEDAUX, J.J.M., 1997b, “Nonlinear Parameter Estimation for

Toxicological Threshold Models”, Ecological Modelling, v. 102, pp. 315-324.

KOSTINA, E., 2004, “Robust Parameter Estimation in Dynamic Systems”,

Optimization and Engineering, v. 5, pp. 461-484.

KRIVY, I., TVRDIK, J., 1995, “The Controlled Random Search Algorithm in

Optimizing Regression Models”, Computational Statistics and Data Analysis, v.

20, pp. 229-234.

KVASNICKA, V., POSPÍCHAL, J., 1997, “A Hybrid of Simplex Method and

Simulated Annealing”, Chemometrics and Intelligent Laboratory Systems, v. 39,

pp. 161-173.

LARENTIS, A.L., 2000, Modelagem do Processo Combinado de Reforma com Dióxido

de Carbono e Oxidação Parcial de Metano. Dissertação de M.Sc., PEQ/COPPE,

Universidade Federal do Rio de Janeiro, Rio de Janeiro, RJ, Brasil.

LARENTIS, A.L., RESENDE, N.S., SALIM, V.M.M., PINTO, J.C., 2001, “Modeling

and Optimization of the Combined Carbon Dioxide Reforming and Partial

Oxidation of Natural Gas”, Applied Catalysis A, v. 215, pp. 211-224.

Page 147: Marcio Schwa Ab

140

LARENTIS, A.L., BENTES JR., A.M.P., RESENDE, N.S., SALIM, V.M.M., PINTO,

J.C., 2003, “Analysis of Experimental Errors in Catalytic Tests for Production of

Synthesis Gas”, Applied Catalysis A, v. 242, pp. 365-379.

LAW, V.J., BAILEY, R.V., 1963, “A Method for the Determination of Approximate

Systems Transfer Functions”, Chemical Engineering Science, v. 18, pp. 189-202.

LEVENBERG, K., 1944, “A method for the solution of certain non-linear problems in

least squares”, Quarterly of Applied Mathematics, v. 2, pp. 164-168.

MARQUARDT, D.W., 1963, “An Algorithm for Least-Squares Estimation of

Nonlinear Parameters”, SIAM Journal, v. 11, pp. 431-441.

MARSEGUERRA, M., ZIO, E., PODOFILLINI, L., 2003, “Model Parameters

Estimation and Sensitivity by Genetic Algorithms”, Annals of Nuclear Energy, v.

30, pp. 1437-1456.

MARSILI-LIBELLI, S., GUERRIZIO, S., CHECCHI, N., 2003, “Confidence Regions

of Estimated Parameters for Ecological Systems”, Ecological Modelling, v. 165,

pp. 127-146.

MEZAKI, R., DRAPER, N.R., JOHNSON, R.A., 1973, “On the Violations of

Assumptions in Nonlinear Least Squares by Interchange of Response and

Predictor Variables”, Industrial and Engineering Chemistry Fundamentals, v. 12,

pp. 251-254.

MOROS, R., KALIES, H., REX, H.G., SCHAFFARCZYK, S., 1996, “A Genetic

Algorithm for Generating Initial Parameter Estimations for Kinetic Models of

Catalytic Processes”, Computers and Chemical Engineering, v. 20, pp. 1257-

1270.

NELDER, J.A., MEAD, R., 1965, A Simplex Method for Function Minimization”,

Computer Journal, v. 7, pp. 308-313.

NORONHA, F.B., PINTO, J.C., MONTEIRO, J.L., LOBÃO, M.W., SANTOS, T.J.,

1993, Um Pacote Computacional para Estimação de Parâmetros e Projeto de

Experimentos, Relatório Técnico PEQ/COPPE, Universidade Federal do Rio de

Janeiro, Rio de Janeiro, RJ, Brasil.

OURIQUE, C.O., BISCAIA JR., E.C., PINTO, J.C., 2002. “The Use of Particle Swarm

Optimization for Dynamical Analysis in Chemical Processes”, Computers and

Chemical Engineering, v. 26, pp. 1783-1793.

Page 148: Marcio Schwa Ab

141

PARK, T.Y., FROMENT, G.F., 1998, “A Hybrid Genetic Algorithm for the Estimation

of Parameters in Detailed Kinetic Models”, Computers and Chemical

Engineering, v. 22, pp. S103-S110.

PARSOPOULOS, K.E., LASKARI, E.C., VRAHATIS, M.N., 2001, “Solving L1 Norm

Errors-In-Variables Problems Using Particle Swarm Optimization”, In: Artificial

Intelligence and Applications, Anahein, USA, pp. 185-190.

PARSOPOULOS, K.E., VRAHATIS, M.N., 2002a, “Recent Approaches to Global

Optimization Problems through Particle Swarm Optimization”, Natural

Computing, v. 1, pp. 235-306.

PARSOPOULOS, K.E., VRAHATIS, M.N., 2002b, “Initializing the Particle Swarm

Optimizer Using the Nonlinear Simplex Method”, In: Advances in Intelligent

Systems, Fuzzy Systems, Evolutionary Computation, pp. 216-221, WSEAS Press.

PATEL, N.R., SMITH, R.L., ZABINSKY, Z.B., 1988, “Pure Adaptive Search in Monte

Carlo Optimization”, Mathematical Programming, v. 43, pp. 317-328.

PETZOLD, L.R., 1989, Dassl Code (Differential Algebraic System Solver), Computing

and Mathematics Research Division, Lawrence Livermore National Laboratory,

Livermore.

POWELL, M.J.D., 1964, “An Efficient Method for Finding the Minimum of a Function

of Several Variables without Calculating Derivatives”, Computer Journal, v. 7,

pp. 155-162.

POWELL, M.J.D., 1965, “A Method for Minimizing a Sum of Squares of Non-Linear

Functions without Calculating Derivatives”, Computer Journal, v. 7, pp. 303-307.

PRATA, D.M., 2005, Reconciliação de Dados em um Reator de Polimerização,

Dissertação de M.Sc., PEQ/COPPE, Universidade Federal do Rio de Janeiro, Rio

de Janeiro, RJ, Brasil

PRICE, W.L., 1976, “A Controlled Random Search Procedure for Global

Optimization”, Computer Journal, v. 20, pp. 367-370.

RATKOWSKY, A.D., 1990, Handbook of Nonlinear Regression Models. New York,

Marcel Dekker.

RICKER, N.L., 1984, “Comparison of Methods for Nonlinear Parameter Estimation”,

Industrial and Engineering Chemistry Process Design and Development, v. 23,

pp. 283-286

ROSENBROCK, H.H., 1960, “An Automatic Method for Finding the Greatest or Least

Value of a Function”, Computer Journal, v. 3, pp. 175-184.

Page 149: Marcio Schwa Ab

142

SANTOS, T.J., PINTO, J.C., 1998, “Taking Variable Correlation into Consideration

During Parameter Estimation”, Brazilian Journal of Chemical Engineering, v. 15,

pp. 1-20.

SEAL, H.L., 1967, “Studies in the History of Probability and Statistics. XV: The

Historical Development of the Gauss Linear Model”, Biometrika, v. 54, pp. 1-24.

SHI, Y., EBERHART, R., 1998a, “A Modified Particle Swarm Optimizer”. In: Proc.

Conference on Evolutionary Computation, Anchorage, Alaska, pp. 69-73.

SHI, Y., EBERHART, R., 1998b, “Parameter Selection in Particle Swarm

Optimization”. In: Proc. Annual Conference on Evolutionary Programming, New

York, EUA, pp. 591-600.

SHI, Y., EBERHART, R., 1999, “Empirical Study of Particle Swarm Optimization”. In:

Proc. Congress on Evolutionary Computation, pp. 1945-1950.

SILVA, L.P., 2003, Síntese e caracterização de poliolefinas obtidas com complexos

diimínicos de níquel contendo pseudo-haletos. Dissertação de M.Sc., IMA,

Universidade Federal do Rio de Janeiro, Rio de Janeiro, RJ, Brasil.

SULIEMAN, H., MCLELLAN, P.J., BACON, D.W., 2001, “A Profile-Base Approach

to Parametric Sensitivity Analysis of Nonlinear Regression Models”,

Technometrics, v. 43, pp. 425-433.

SULIEMAN, H., MCLELLAN, P.J., BACON, D.W., 2001, “A Profile-Base Approach

to Parametric Sensitivity in Multiresponse Regression Models”, Computational

Statistics & Data Analysis, v. 45, pp. 721-740.

TRELEA, I. C., 2003, “The Particle Swarm Optimization Algorithm: Convergence

Analysis and Parameter Selection”, Information Processing Letters, v. 85, pp.

317-325.

TSALLIS, C., STARIOLO, D.A., 1996, “Generalized Simulated Annealing”, Physica

A, v. 233, pp. 395-406.

VANROLLEGHEM, P.A., KEESMAN, K.J., 1996, “Identification of Biodegradation

Models Under Model and Data Uncertainty”, Water Science and Technology, v.

33, pp. 91-105.

WILLIAMS, E.J., 1962, “Exact Fiducial Limits in Non-Linear Estimation”, Journal of

the Royal Statistical Society B, v. 24, pp. 125-139.

ZABINSKY, Z.B., SMITH, R.L., 1992, “Pure Adaptive Search in Global

Optimization”, Mathematical Programming, v. 53, pp. 323-338.

Page 150: Marcio Schwa Ab

143

APÊNDICE

Considerando o conjunto de transformações sucessivas descritas abaixo:

0

1

1

0 1

1 2

1

1

n

m

k

k

kn n

km m

C C

C C

C C

C C−

+

(A.01)

onde Cn corresponde a concentração presente da n-ésima forma do catalisador, kn a

constante de velocidade específica de transformação da forma n na forma posterior n+1

e m é o número de transformações que ocorrem e também corresponde ao número de

espécies ativas, é possível escrever o seguinte conjunto de equações de balanço:

( )00 0 00 ; 0 1dC k C C

dt+ = = (A.02)

( )11 1 0 0 1 ; 0 0dC k C k C C

dt+ = = (A.03)

( )1 1 ; 0 0nn n n n n

dC k C k C Cdt − −+ = = (A.04)

( )1 1 ; 0 0mm m m

dC k C Cdt − −= = (A.05)

As soluções das Equações (A.02) e (A.03) são respectivamente:

( ) 00

k tC t e−= (A.06)

( ) ( )0 101

1 0

k t k tkC t e ek k

− −= −−

(A.07)

A solução da Equação (A.04) pode ser obtida observando-se que esta é uma

equação diferencial ordinária, linear, não homogênea e com coeficientes constantes.

Neste caso, a solução pode ser colocada na forma:

Page 151: Marcio Schwa Ab

144

( ) * nk tn n nC t C e Sp−= + (A.08)

onde En é uma constante de integração e Spn é uma solução particular, obtida a partir da

solução de Cn-1(t). Observe que a solução Cn(t) tem necessariamente a forma de uma

soma de exponenciais, em função da natureza recursiva do conjunto de equações e em

função da forma das Equações (A.06) e (A.07). Assim, supondo que a solução Cn(t)

possa ser escrita na forma:

( )0

i

nk tn

n i ni

C t A e b−

== +∑ (A.09)

de modo que:

( )

0

i

nn k tn

i ii

dC tA k e

dt−

==−∑ (A.10)

Substituindo as Equações (A.09) e (A.10) na Equação (A.04) obtém-se:

1

11

1 1 10 0 0

i i i

n n nk t k t k tn n n

i i n i n n n i n ni i i

A k e k A e k b k A e k b−

−− − −−

− − −= = =

− + + = +∑ ∑ ∑ (A.11)

e agrupando os dois primeiros termos da equação acima tem-se:

( ) 1

11

1 1 10 0

i i

n nk t k tn n

i n i n n i n n ni i

A k k e k b A k e k b−

−− −−

− − −= =

− + = +∑ ∑ (A.12)

Para que a Equação (A.12) seja satisfeita, as seguintes igualdades devem ser

observadas:

1 1n n n nk b k b− −= (A.13)

( ) 11

n ni n i i nA k k A k−

−− = (A.14)

A Equação (A.14) pode ser escrita na seguinte forma recursiva:

( )

1 1n n ni i

n i

kA Ak k

− −=−

(A.15)

para i = 0 ... n-1. As soluções para C0(t) e C1(t) (Equações (A.06) e (A.07)) indicam que

b0 e b1 são iguais a 0 e assim, de acordo com a Equação (A.13), todos os bi’s são nulos.

A Equação (A.06) ainda mostra que 00A = 1.

Page 152: Marcio Schwa Ab

145

Além disso, a condição inicial da Equação (A.04) é Cn(0) = 0 e substituindo na

Equação (A.09), lembrando que bn =0, obtém-se:

0

0n

ni

i

A=

=∑ (A.16)

ou ainda

1

0

nn nn i

i

A A−

==−∑ (A.17)

Deste modo, a concentração de uma forma n qualquer ao longo do tempo pode

ser obtida pela seguinte equação:

( )0

i

nk tn

n ii

C t A e−

==∑ (A.18)

sendo os coeficientes calculados através das seguintes equações:

( )

00

1 1

1

0

1

; 0.. -1 ; 0n n ni i

n i

nn nn i

i

AkA A i n n

k k

A A

− −

=

=

= = >−

=−∑

(A.19)

Finalmente, a solução da Equação (A.05) pode ser obtida pela substituição da

solução de Cm-1(t) e posterior por integração, como segue:

1

11

0

i

mk tmm

m ii

dC k A edt

−−−

−=

= ∑ (A.20)

1

110 0 0

mi

mC t kmm i m

i

d A k e dτς τ−

−−−

== ∑∫ ∫ (A.21)

( )1 1

1 11 1

0 0

i

m mk tm mm m

m i ii ii i

k kC t A e Ak k

− −−− −− −

= == −

− −∑ ∑ (A.22)

A Equação (A.22) pode ser reescrita com a inserção da constante cinética de

transformação km, lembrando que esta constante é igual a 0, já que a m-ésima forma não

sofre transformações.

( )1 1

1 11 1

0 0

i

m mk tm mm m

m i ii im i m i

k kC t A e Ak k k k

− −−− −− −

= == −

− −∑ ∑ (A.23)

Page 153: Marcio Schwa Ab

146

O segundo termo do lado direito da Equação (A.23) é, de acordo com as

Equações (A.15) e (A.17), igual a mmA . Utilizando ainda a Equação (A.15) o primeiro

temor do lado direito também pode ser simplificado e, assim, a Equação (A.23) pode ser

reescrita com a seguinte forma:

( )1

1

0

i

mk tm m

m i ii

C t A e A−

− −

== +∑ (A.24)

ou ainda:

( )0

i

mk tm

m ii

C t A e−

==∑ (A.25)

que uma equação com a mesma forma da solução genérica para uma forma

intermediária n qualquer. Assim, a concentração de qualquer forma do catalisador pode

ser representada pela Equação (A.18), sempre lembrando que para a última forma o

valor de km é nulo.

Sendo a taxa de polimerização catalisada por uma forma n qualquer do

catalisador proporcional à concentração desta forma, a taxa global de polimerização

pode ser definida como:

( ) ( )0 0 0

i

m m nk tn

n n n in n i

Rp t Kp C t Kp A e−

= = == =∑ ∑ ∑ (A.26)

onde Kpn é a constante da taxa de polimerização associada a n-ésima forma do

catalisador. A Equação (A.26) ainda pode ser reescrita em uma forma mais concisa,

como segue:

( )0

i

m mk tn

n ii n i

R t Kp A e−

= =

=

∑ ∑ (A.27)

Desta forma a taxa global de polimerização é descrita pela Equação (A.27),

sendo os coeficientes calculados pelas Equações (A.19) e lembrando que km = 0. Além

disso, o valor de Kp0 também deve ser nulo, já que o catalisador em sua forma inicial

não é ativo. Assim, o número de parâmetros que devem ser estimados é igual a 2m,

sendo m constantes cinéticas de transformação de espécie do catalisador (k0, k1, ..., km-1)

e m constantes cinéticas de polimerização (Kp1, Kp2, ..., Kpm).

Considerando um caso muito particular em que k0 = k1 = ... = kn = ... = km-1 = k, a

Equação (A.04) tem solução genérica:

Page 154: Marcio Schwa Ab

147

( ) * ktn n nC t C e Sp−= + (A.28)

onde Cn* é uma constante de integração e Spn é uma solução particular. Como a

constante de transformação k é a mesma para todas as espécies, a solução particular tem

a forma

0

ni kt

n ii

Sp At e−

=

=∑ (A.29)

gerando a seguinte solução genérica:

( ) ( )

0

nn i kt

n ii

C t A t e−

=

=∑ (A.30)

Substituindo a equação (A.30) na equação (A.04) temos:

( ) ( ) ( ) ( ) ( )1

11 1

0 0 0 0

n n n nn n n ni kt i kt i kt i kt

i i i ii i i i

A it e A t e k kA t e kA t e−

−− − − − − −

= = = =

+ − + =∑ ∑ ∑ ∑ (A.31)

Agrupando os termos de igual potência de t temos:

( ) ( ) ( ) ( ) ( ) ( ) ( )( )1

11

01 0

nn n n n n ni kt n kt

i i i i n ni

A i kA kA kA t e kA kA t e−

− − −+

=

+ − + − + − + = ∑ (A.32)

de onde é possível obter:

( ) ( )11 1

n ni i

kA Ai

−+ =

+ (A.33)

Da condição inicial

( ) ( ) ( )00

00 0 0 0

nn ni k

n ii

C A e A−

=

= = → =∑ (A.34)

Sabendo que ( )00 1A = (Equação (A.06)), conclui-se que

( ) ( )

( ) ( ) ( )

( ) ( )

1 10 1

22 2 1

0 1 1

0

0,

0, 0, 2

0, , !

nn n

n

A A k

kA A A

kA i n An

= =

= = =

= ≠ =

(A.35)

de maneira que

Page 155: Marcio Schwa Ab

148

( ) , 1 1!

nn kt

nkC t t e n mn

−= ≤ ≤ − (A.36)

que é a solução de Poisson.

Para n = m, a Equação (A.05) leva a

( ) ( )1 ; 0 0

-1 !

mm ktm

mdC k t e Cdt m

− −= = (A.37)

Postulando a solução

( ) ( )

0

mm i kt

m i mi

C t A t e b−

=

= +∑ (A.38)

e substituindo a Equação (A.38) na Equação (A.37) temos:

( ) ( ) ( ) ( )1 1

0 0 -1 !

mm mm mi kt i kt m kt

i ii i

kA it e A t e k t em

− − − − −

= =

+ − =∑ ∑ (A.39)

( ) ( ) ( ) ( )

( )1 1

11

0 01

-1 !

mm mm m mi kt i kt m kt m kt

i i mi i

kA i t e kA t e kA t e t em

− −− − − − −

+= =

+ − − =∑ ∑ (A.40)

de onde se conclui que

( )

( ) ( )

( )

( )

1

1

1

0;

, 1... 2;1

1 !

mm

m mi i

mm

m

AkA A i m

ikA

m

+

=

= = −+

= −−

(A.41)

Da equação (A.41)

( )

!

im

ikAi

= − (A.42)

de maneira que, da condição inicial, temos bm = 1 e a solução é:

( )1

01

!

imi kt

mi

kC t t ei

−−

=

= −∑ (A.43)

Dessa forma, se a taxa de polimerização ocorre com velocidade proporcional às

quantidades de cada espécie:

Page 156: Marcio Schwa Ab

149

1 1

0 0 01

! !

n nm m mn kt n kt

n n n mn n n

k kR Kp C Kp t e Kp t en n

− −− −

= = =

= = + −

∑ ∑ ∑ (A.44)

( )1

0 !

nmn kt

n m mn

kR Kp Kp t e Kpn

−−

=

= − +∑ (A.45)

Neste modelo o número de parâmetros que são estimados é igual a m+1, sendo 1

constante cinética de transformação de espécie do catalisador e m constantes cinéticas

de polimerização (Kp1, Kp2, ..., Kpm), lembrando que Kp0 é igual a zero.