Dissertação_Otimização Multiobjetivo Robusta de Estruturas

UNIVERSIDADE FEDERAL DE PERNAMBUCO CENTRO DE TECNOLOGIA E GEOCIÊNCIAS DEPARTAMENTO DE ENGENHARIA CIVIL

OTIMIZAÇÃO ROBUSTA DE ESTRUTURAS UTILIZANDO O MÉTODO DA BASE REDUZIDA

Renato de Siqueira Motta

Dissertação submetida ao Corpo Docente do Curso de Pós-Graduação da Universidade Federal de Pernambuco, como parte dos requisitos necessários à obtenção do Grau de Mestre em Ciências em Engenharia Civil.

Orientador: Silvana Maria Bastos Afonso da Silva Co-orientador: Paulo Roberto Maciel Lyra

Recife, Pernambuco – Brasil Novembro de 2009.

OTIMIZAÇÃO ROBUSTA DE ESTRUTURAS UTILIZANDO O MÉTODO DA BASE REDUZIDA

Renato de Siqueira Motta

Dissertação submetida ao Corpo Docente do Curso de Pós-Graduação da Universidade Federal de Pernambuco, como parte dos requisitos necessários à obtenção do Grau de Mestre em Ciências em Engenharia Civil.

Aprovada por: _____________________________________________ Prof. Silvana Maria Bastos Afonso da Silva, Ph.D. Orientador - UFPE

_____________________________________________ Prof. Paulo Roberto Maciel Lyra, Ph.D. Co-orientador - UFPE

_____________________________________________ Prof. Bernardo Horowitz, Ph.D. Examinador Interno - UFPE

_____________________________________________ Prof. Ramiro Brito Willmersdorf, Ph.D. Examinador Externo - UFPE

_____________________________________________ Prof. Luis Eloy Vaz, Dr.-Ing. Examinador Externo - UFRJ

Recife, Pernambuco – Brasil Novembro de 2009.

“No processo de evolução das espécies, somos um ótimo local” iii   

MÉTODO DA BASE REDUZIDA PARA OTIMIZAÇÃO ROBUSTA DE ESTRUTURAS Renato de Siqueira Motta Resumo Com o rápido aumento da capacidade computacional, o tema otimização avançou de maneira notável nos últimos anos. Atualmente inúmeras aplicações de projetos ótimos em diferentes especialidades, como mecânica estrutural, custos de produção, escoamento de fluidos, acústica, etc. têm sido descritas na literatura. Entretanto, na maioria das aplicações da engenharia, a abordagem tradicional é considerar modelos e parâmetros determinísticos. Infelizmente a abordagem determinística pode levar a soluções cujo desempenho pode cair significativamente devido às perturbações decorrentes das incertezas. Nestas circunstâncias, um objetivo melhor seria um projeto ótimo que tenha um alto grau de robustez. O processo de encontrar este ótimo é chamado Otimização Robusta (OR). Aqui, abordaremos duas técnicas para a análise de propagação de incerteza, não intrusivas, que utiliza modelos computacionais determinísticos: o método de Monte Carlo (MC) e o método da Colocação Probabilística (“Probabilistic Collocation Method”) (PCM). A análise de propagação de incerteza essencialmente envolve o cálculo de momentos estatísticos da função de interesse. Várias medidas de robustez têm sido propostas na literatura, em particular, o valor médio e o desvio padrão da função envolvida no problema de otimização serão considerados aqui. Quando estas medidas de robustez são usadas combinadas, a procura de projetos ótimos robustos surge como um problema de Otimização Multiobjetivo Robusta (OMR). Técnicas de Otimização Multiobjetiva permitem o projetista modelar um problema específico considerando um comportamento mais realista, o qual comumente envolve o atendimento de vários objetivos simultaneamente. O procedimento adequado, quando um problema multiobjetivo precisa ser resolvido, é determinar a fronteira de Pareto. Nos últimos 15 anos, distribuições eficientes de pontos de Pareto têm sido obtidas através de novos algoritmos como o NBI (Normal-Boundary Intersection) e o NNC (Normalized Normal-Constraint). Estas estratégias são implementadas aqui, junto com outras abordagens comumente utilizadas na literatura, como o método da soma ponderada e o método Min-Max. Como a geração de pontos de Pareto e a análise de incerteza podem ser muito custosas, técnicas de aproximação, baseada no uso do Método da Base Reduzida (MBR), são incorporadas ao nosso procedimento. O propósito do método é obter um modelo de alta fidelidade com custo computacional aceitável. Além disto, uma estratégia de separabilidade com uma decomposição afim, permite o desenvolvimento de uma estratégia eficiente de cálculo “off-line/on-line”, para a implementação computacional do MBR. Problemas contínuos em duas dimensões submetidos a carregamentos estáticos e térmicos são as aplicações consideradas neste trabalho, os desempenhos das diferentes estratégias examinadas são comparadas. A combinação das várias técnicas de aproximação descritas permitiu a obtenção das soluções OMR em pouco tempo computacional. Palavras chaves: Otimização robusta, método da Colocação Probabilística, Otimização Multiobjetivo, Método da base reduzida. iv   

REDUCED-BASIS METHOD FOR STRUCTURAL ROBUST OPTIMIZATION

Renato de Siqueira Motta Abstract The topic of optimization has advanced in a remarkable manner in recent years, with the rapid growth of computational power. Nowadays a number of applications on design optimization from different disciplines such as structural mechanics, acoustics, etc. have been already reported in literature. However, in most engineering applications, the traditional approach is to consider deterministic models and parameters. Unfortunately, the deterministic approach generally leads to a final design whose performance may degrade significantly because of perturbations arising from uncertainties. In this scenario a better target that provides an optimal design is one that gives a high degree of robustness. The process to find such optimal is referred to as robust optimization (RO). Here, we discuss two nonintrusive uncertainty propagation analysis techniques that exploit deterministic computer models: Monte Carlo (MC) method and Probabilistic Collocation Method (PCM). The uncertainty propagation analysis essentially involves computing the statistical moments of the output. Several robustness measures have been proposed in the literature, in particular, the expected value and standard deviation of the function involved in the optimization problem are considered here. When using these robustness measures combined the search of optimal robust design appears as a robust multiobjetive optimization (RMO) problem. Multiobjective optimization techniques allow a designer to model a specific problem considering a more realistic behavior, which commonly involves the satisfaction of several targets simultaneously. The computation of the Pareto front solutions is the adequate procedure when a multiobjective problem has to be solved. In the last 15 years efficient Pareto point distribution has been obtained through the use of new algorithms such as NBI (Normal-Boundary Intersection), and NNC (Normalized Normal-Constraint). These two strategies are implemented in this work together with other commonly considered approaches in literature such as weighted sum method and min-max method. As the generation of Pareto points and the uncertainty analysis could be very costly, approximation techniques based on the use of Reduced Basis Method (RBM) are also incorporated in our procedure. The purpose of such scheme is to get high fidelity model information with acceptable computational expense. Moreover, a parameter separability strategy together with the affine decomposition allows the development of an efficient off-line/on-line calculation strategy for the computational implementation of the RBM. Two dimensional continua problems under static loads are the applications addressed in this work and the performance of the different strategies discussed are compared.  The combination of all the approximate methodologies described in this work allows the computations of RMO solutions, with very low computational time. Keywords: Robust optimization, Probabilistic Collocation Method, Multiobjetive optimization, Reduced basis method v   

SUMÁRIO LISTA DE SÍMBOLOS 1. INTRODUÇÃO

ix 1

1.1. MOTIVAÇÃO E CONSIDERAÇÕES INICIAIS

1

1.2. OBJETIVOS

3

1.3. ORGANIZAÇÃO DA DISSERTAÇÃO

5

1.4. REFERÊNCIAS

5

2. PROBLEMAS BIDIMENSIONAIS 2.1. FORMULAÇÃO DE PROBLEMAS TÉRMICOS 2.1.1. Método dos Elementos Finitos aplicado à problemas térmicos 2.2. FORMULAÇÕES DOS PROBLEMAS ELÁSTICOS

8 8 9 11

2.2.1. Método dos Elementos Finitos aplicado à elasticidade bidimensional 13 2.2.2. Acoplamento termo-elástico 2.3. MAPEAMENTO DAS EQUAÇÕES DO MEF

14 15

2.3.1. Mapeamento da equação térmica

16

2.3.2. Mapeamento da equação elástica

18

2.4. MÉTODO DA BASE REDUZIDA 2.4.1. Procedimento Computacional Off-Line/On-Line 2.5. EXEMPLOS

20 23 24

2.5.1. Placa quadrada com orifício central

24

2.5.2. Exemplo termo-elástico acoplado

28

2.6. REFERÊNCIAS 3. OTIMIZAÇÃO ESCALAR E ANÁLISE DE SENSIBILIDADE

32

34

3.1. INTRODUÇÃO

34

3.2. FORMULAÇÃO MATEMÁTICA

35

3.3. PROGRAMAÇÃO MATEMÁTICA

35

3.3.1. Condições de ótimo para problemas com restrições

36

3.3.2. Programação Quadrática Sequencial - SQP

38

3.4. ANÁLISE DE SENSIBILIDADES

40

3.4.1. Método das diferenças finitas globais

41

3.4.2. Método direto

41

3.5. INTEGRAÇÃO ANÁLISE/OTIMIZAÇÃO 3.5.1. O método da base reduzida no procedimento de otimização (a) Análise de sensibilidade através do RBM

42 43 43

3.6. EXEMPLOS

45

3.7. REFERÊNCIAS

48

4. OTIMIZAÇÃO MULTIOBJETIVO

52

4.1. DEFINIÇÃO DO PROBLEMA

52

4.2. CONCEITO DE ÓTIMO DE PARETO

53

4.3. MÉTODOS PARA GERAÇÃO DE PONTOS DE PARETO

54

4.3.1. Método da Soma Ponderada dos Objetivos

54

4.3.2. Método Min-Max

56

4.3.3. Método da Interseção Contorno-Normal (NBI)

56

4.3.4. Método da Restrição Normal Normalizada (NNC)

59

4.3.5. Diferenças entre o NBI e o NNC

61

4.4. EXEMPLOS

62

4.4.1. Treliça de 3 barras

62

4.4.2. Placa quadrada com um orifício central

63

4.4.3. Placa engastada - problema aclopado

66

4.5. PROBLEMAS COM MAIS DE DUAS FUNÇÕES OBJETIVO

68

4.5.1. Exemplo geométrico com 3 funções objetivo

72

4.5.2. Exemplo analítico com 3 funções objetivo

73

4.5.3. Placa sob ação termo-estrutural acoplada com 3 funções objetivo

75

4.6. REFERÊNCIAS 5. OTIMIZAÇÃO CONSIDERANDO INCERTEZAS

77

80

5.1. INTRODUÇÃO

80

5.2. TEORIA PROBABILÍSTICA E ESTATÍSTICA

81

5.2.1. Variáveis Aleatórias

82

5.2.2. Distribuições de probabilidade

83

5.3. CÁLCULO DAS ESTATÍSTICAS

85

5.3.1. Método de Monte Carlo

86

5.3.2. Técnicas de amostragem

89

(a) Exemplo - MC por diferentes amostragens

90

5.3.3. Método da Colocação Probabilística (PCM)

92

(a) Polinômios ortogonais

92

(b) Quadratura de Gauss

93

(c) Aplicando Quadratura de Gauss à estatística - PCM

95

(d) Implementação Computacional

99

5.3.4. Exemplos

101

(a) Verificação - PCM

101

(b) Função Periódica

103

(c) Função com singularidade

104

5.4. OTIMIZAÇÃO ROBUSTA

109

5.4.1. Medidas de Robustez

109

5.4.2. Exemplo: Placa quadrada com orifício central

112

(a) Determinação das amostras

112

(b) Resultado da otimização

114

5.5. REFERÊNCIAS 6. CONCLUSÃO E TRABALHOS FUTUROS

115

118

6.1. PRINCIPAIS DESENVOLVIMENTOS

118

6.2. CONCLUSÕES DOS RESULTADOS OBTIDOS

118

6.3. TRABALHOS FUTUROS

120

LISTA DE SÍMBOLOS ABREVIATURAS E SIGLAS Capítulo 2 MEF

Método dos Elementos Finitos.

EF

Elementos Finitos.

Q4

Elemento retangular isoparamétrico linear de quatro nós.

T3

Elemento triangular linear de três nós (vide CST).

CST

“Constant Strain Triangle” – Elemento triangular de deformação constante.

MBR

Método da Base Reduzida.

NGL

Número de Graus de Liberdade

Le

Tamanho médio dos elementos

Capítulo 3 SSO

“Structural Shape Optimization” - Otimização estrutural de forma.

SQP

“Sequential Quadratic Programming” - Programação sequencial quadrática.

KKT

Karush-Kuhn-Tucker – Relativo às condições de ótimo.

BFGS

Broyden-Fletcher-Goldfarb-Shann - Métodos para aproximação da matriz Hessiana

MDF

Método das Diferenças Finitas

Capítulo 4 OM

Otimização Multiobjetivo

POM

Problema de Otimização Multiobjetivo

WS

“Weight Sum” – Método da soma ponderada

NBI

“Normal-Boundary Intersection” - Interseção Contorno-Normal

NBIm

NBI modificado (vide NBI)

ECMI

Envoltória Convexa de Mínimos Individuais

NNC

“Normalized Normal Constraint” - Método da Restrição Normal Normalizada

NC

“Normal Constraint” - Restrição Normal ix 

 

NU

Normal à linha utópica

Capítulo 5  OMR

Otimização Multiobjetivo Robusta

OR

Otimização Robusta

PDF

“Probability Density Function” - Função densidade de probabilidade

CDF

“Cumulative Distribution Function” - Função de distribuição acumulada

MC

Método de Monte Carlo

DoE

“Design of experiments” – Plano de amostragem

LHS

“Latin Hipercube Sampling” - Método para gerar amostras aleatórias

MT

“Mersenne Twister” - método para a geração de amostras pseudo-aleatória

PCM

“Probabilistic Collocation Method” - Método da colocação probabilística

ME-PCM “Multi-Element Probabilistic Collocation Method”

LISTA DE SÍMBOLOS ROMANOS  

Escalares c1

Constante de Lamé

c2

Constante de Lamé

Dijkl

Componente do tensor de elasticidade

E Fobj(x) gi(x) hi(x) kij N

Módulo de elasticidade do material Função objetiva Restrição de desigualdade Restrição de igualdade Termos do tensor de condutividade térmica Tamanho da Base

Ni

Função de forma do nó i

f (μ )

Flexibilidade

f N (μ)

Flexibilidade aproximada

L(x,α,β,λ) ti Ti

xri x kl

Função Lagrangeana Largura da região Temperatura no nó i Coordenada de um nó i na sub-região r Limite inferior da variável de projeto x 

 

x ku

Limite superior da variável de projeto

Vetores b br(μ)

de ε fn

Forças de volume agindo sobre o domínio V* Forças de volume agindo sobre o domínio computacional para uma subregião r Vetor de deslocamentos nodais de um elemento Vetor de deformação Forças normais aplicadas no contorno Γ*

ft

Forças tangenciais e tangentes aplicadas no contorno Γ*

f (μ)

Forças normais aplicadas no contorno Γ para uma sub-região r

f tr (μ)

Forças tangenciais aplicadas no contorno Γ para uma sub-região r

r FTE

Carregamento devido à deformação térmica

F F*

Vetor de forças, Vetor das funções objetivo Vetor do mínimo individual das funções objetivo

F N (μ)

Vetor de carregamentos transformado

r n

*r

F F jr

Vetor de carregamentos do domínio real de uma sub-região r

Fˆ

Vetor de Pseudo-força Fluxo de calor no contorno Calor interno gerado sobre o domínio V*,

q qv

S TN

Vetor de carregamento independente de μ para uma sub-região r

Vetor da direção de busca

u

Vetor de temperaturas aproximado Vetor de deslocamento

u N (μ )

Deslocamento aproximado

T0

Temperatura inicial

Matrizes B *e

C(μ) D D* G Gr

Matriz que relaciona as componentes da derivada de deslocamento com deslocamentos Operador simétrico usado para relaxação Matriz de elasticidade Pseudo-matriz de elasticidade no domínio real Matriz de transformação

H I

Matriz de transformação de uma sub-região r Hessiana da função Lagrangeana Matriz identidade

Ks

Matriz de rigidez

K *s

Matriz de rigidez no domínio real xi 

 

K *s N K

*r s

Matriz de rigidez no domínio real no espaço reduzido Matriz de rigidez do domínio real de uma sub-região r

K rs j

Matriz de rigidez independente de μ para uma sub-região r

KT

Matriz de “rigidez térmica” ou condutividade térmica generalizada

K *T

Matriz de “rigidez térmica” transformada

K

*N T

Matriz de “rigidez térmica” transformada no espaço reduzido

K rk j

Matriz de condutividade independente de μ

K Γr c

Matriz de convecção no contorno de referência

K

Matriz de convecção no contorno de referência no espaço reduzido

Nr Γc

K Ωr c

Matriz de convecção no domínio de referência

K ΩNrc

Matriz de convecção no domínio de referência no espaço reduzido

Kk N

Matriz de condutividade térmica Matriz de função de forma Conjunto de amostras para campo de solução Espaço da Base Reduzida Matriz compostas por soluções base do espaço reduzido

SN WN

Z

GREGOS Escalares α

Passo da busca linear

αΩ

Coeficiente de convecção no domínio

αΩ

Coeficiente de convecção no contorno

βi

Componentes do vetor com os pesos para os problemas multiobjetivo

β sr ( μ ) j βTr ( μ ) j

Termos que dependem do parâmetro μ da transformação da Matriz de rigidez Termos que dependem do parâmetro μ da transformação da Matriz de condutividade térmica

Δxi

Pertubação da variável

λ, β, φ Γ μ

Multiplicadores de Lagrange Contorno Parâmetros variáveis

μ low

Limite inferior do espaço de projeto

μ

Limite superior do espaço de projeto

up

(μ ) ϕ (μ ) ϕ

r Ω

Termos da transformação no domínio que dependem de μ

r Γ

Termos da transformação no contorno que dependem de μ xii 

 

Ω Ω*

Domínio Domínio real (transformado)

σ ij ν γ

Tensão Coeficiente de Poisson Coeficiente de expansão térmica

Vetores α(μ)

ε

* 0

Coeficiente linear da aproximação no espaço reduzido Deformações iniciais (térmicas)

μ

Vetor dos Multiplicadores de Lagrange das funções restrições para um x qualquer Vetor dos Multiplicadores de Lagrange das funções restrições para x* ótimo Vetor de parâmetros variáveis

ζi

Conjunto de soluções do espaço

∇N i

Derivadas da função de forma para um nó i

ϕ, λ ϕ *, λ*

Matrizes Φ

σ

Matriz auxiliar na definição dos pontos que compõem a ECMI *

Tensor de tensões

xiii   

CAPÍTULO 1 INTRODUÇÃO  1 INTRODUÇÃO

1.1 MOTIVAÇÃO E CONSIDERAÇÕES INICIAIS Técnicas de otimização têm sido extensivamente usadas para obter projetos viáveis e econômicos nos mais variados campos das Engenharias. Atualmente as abordagens têm se tornado cada vez mais realistas, tendo sido comumente empregadas na solução de problemas não triviais da engenharia prática, incluindo o projeto de estruturas ótimas através de simulação computacional. No entanto, dois pontos apenas recentemente têm sido abordados de forma mais incisiva. A primeira questão diz respeito a atender simultaneamente várias metas (funções objetivo) em geral conflitantes, além das várias restrições a serem satisfeitas, quase sempre envolvidas no desenvolvimento de projetos ótimos de problemas reais. Otimizadores de propósito geral não resolvem tais problemas. Uma classe de estratégias baseadas no denominado conceito de Pareto (ARORA et al., 2007), constitui a abordagem adequada quando problemas de otimização multiobjetivo (OM) devem ser resolvidos. Nos últimos 15 anos distribuições eficientes de pontos de Pareto têm sido obtidas graças ao desenvolvimento de algoritmos eficientes tais como o NBI (“NormalBoundary Intersection”) (DAS e DENNIS, 1996) e o NNC (“Normalized Normal Constraint”) (MESSAC et al, 2003). Essas estratégias juntamente com outras abordagens de mais fácil implementação (Método da soma ponderada e Método min-max) (ARORA et al., 2007) são aqui implementados e analisadas. Outra questão importante, são as incertezas embutidas de várias formas no problema de otimização. Na maioria das aplicações na engenharia, a abordagem tradicional é considerar modelos determinísticos. Porém, algum grau de incerteza ou variação nas características de qualquer sistema estrutural é inevitável. Infelizmente, a abordagem determinística pode levar a soluções cujo desempenho pode cair significativamente devido às perturbações decorrentes de incertezas. A Otimização Robusta (OR) leva em consideração as variáveis incertas (aleatórias) e suas probabilidades de ocorrência, de modo a encontrar um ótimo menos vulnerável a variação dos parâmetros incertos. Nesta dissertação, serão examinadas algumas abordagens para a consideração das incertezas no processo de otimização e assim obter projetos robustos. As medidas de robustez utilizadas aqui foram: a esperança e a variância da função de interesse. Quando usa-se estas medidas, a busca por um projeto robusto ótimo, surge com um problema de decisão com múltiplos critérios (otimização multi-objetivo robusta - OMR). 1   

Neste trabalho, uma ferramenta para obter projetos ótimos robustos, sobre vários critérios, é descrita, implementada e aplicada no âmbito estrutural. Para o cálculo dos parâmetros estatísticos serão empregadas duas técnicas não intrusivas de análise de propagação de incerteza, o método de Monte Carlo (MC) e o método da colocação probabilística (“Probabilistic Collocation Method” - PCM) (RAMAMURTHY, 2005). O MC é um dos métodos mais tradicionais para este tipo de análise, porém é inviável para ser aplicado diretamente em modelos complexos de alta fidelidade, devido ao grande número de análises necessárias. O PCM é uma ferramenta desenvolvida visando uma análise de incerteza eficiente, mesmo em modelos complexos e computacionalmente custosos. A idéia básica do PCM é aproximar a resposta do modelo em função das variáveis aleatórias, por funções polinomiais, e estimar os parâmetros estatísticos por integração numérica. Foi ainda estudado o comportamento do MC para o cálculo das estatísticas de funções especificas, através do uso de diferentes técnicas de amostragem. Os problemas de otimização estrutural aqui considerados, abordam modelos elásticos lineares e/ou térmicos bidimensionais, além de treliças planas. Para sua solução, foi utilizado um algoritmo padrão de otimização de forma (“Structural Shape Optimization” - SSO) integrando aos procedimentos de otimização, a modelagem geométrica, a análise estrutural (ou térmica) e a análise de sensibilidade, para o cálculo dos gradientes requeridos pelo otimizador. Este algoritmo usa a programação sequencial quadrática (“Sequential Quadratic Programming” - SQP) (VANDERPLAATS, 2004) como otimizador e na versão inicial do algoritmo SSO, cada análise estrutural (e/ou térmica) é calculada através de simulação computacional utilizando o Método dos Elementos Finitos (MEF) (COOK et al, 2003). O algoritmo SSO será aproveitado e adaptado para os métodos de obtenção de pontos de Pareto, e nele também serão introduzidos os procedimentos para os cálculos dos parâmetros estatísticos das funções de interesse. A depender da aplicação, o custo das análises dos modelos, bem como da análise de sensibilidade via MEF, pode ser elevado. No entanto, as obtenções dos pontos de Pareto para os problemas de OM e principalmente nos problemas de OMR estão associados a um grande número de avaliações de funções e cálculos de gradientes, assim sendo o tempo total de CPU no procedimento de otimização pode ser bastante elevado. Para amenizar esta dificuldade, técnicas de aproximação baseadas no método das bases reduzidas (MBR) (AFONSO, 2003; ALBUQUERQUE, 2005; AFONSO et al., 2009) são inseridas na metodologia, para as etapas de análise e de análise de sensibilidade requeridas pelo algoritmo, durante o processo de otimização, produzindo um grande ganho computacional, sem grande perda de fidelidade do modelo. O MBR aqui implementado, é uma projeção do tipo Galerkin em um espaço de aproximação de baixa ordem, formado por soluções do modelo de alta fidelidade, para o problema de interesse em pontos selecionados do espaço de projeto. O procedimento combina dois aspectos, a saber: uma aproximação com precisão e uma melhoria na eficiência computacional. No MBR utiliza-se do mapeamento da geometria das regiões do modelo para aplicar uma estratégia de separabilidade de parâmetros (variáveis aleatórias, variáveis de projeto). Todos os cálculos são conduzidos em um domínio fixo denominado domínio 2   

de referência. Transformações geométricas entre o domínio real e computacional são conduzidas e inseridas nas equações governantes de cada problema específico. Isto, junto com a decomposição dos termos da rigidez e de força permite o desenvolvimento de um algoritmo para implementação computacional do método dividido em dois estágios denominados “off-line” e “on-line”. Os cálculos “off-line” são conduzidos somente uma vez e usados subseqüentemente no estágio “on-line” para cada novo parâmetro desejado. Com isso, para cada novo projeto, funções e gradientes são obtidos muito rapidamente. A conseqüência de avaliações de funções (e seus gradientes) precisas e de baixo custo computacional torna este procedimento bem atrativo para propósitos de OMR. Vários exemplos utilizando as várias estratégias serão investigados com o objetivo de se verificar a efetividade e robustez dos procedimentos desenvolvidos. Será demonstrado que a união de estratégias eficientes de OM e OMR e procedimentos de aproximação constitui uma ferramenta efetiva e confiável para se obter os pontos de Pareto nos problemas aqui estudados e por seguinte dá indícios de sucesso para outros problemas mais complexos.

1.2 OBJETIVOS O objetivo final desta dissertação é investigar, desenvolver e implementar procedimentos gerais para a otimização multiobjetivo convencional e robusta, incorporando um procedimento de aproximação (MBR), aplicado na análise estrutural e/ou térmica e análise de sensibilidade. Na Figura 1.1 é apresentado o fluxograma do processo de otimização para as diferentes configurações possíveis implementadas. Os principais aspectos relacionados ao desenvolvimento deste trabalho são: •

Adaptar e implementar rotinas para análise estrutural e/ou térmica e análise de sensibilidade via Método dos Elementos Finitos (aqui denominado modelo real);

•

Automatizar, desenvolver e implementar os procedimentos para análise estrutural e/ou térmica e o cálculo das sensibilidades no contexto do MBR (aqui denominado modelo aproximado);

•

Realizar um estudo comparativo dos dois procedimentos descritos nos itens anteriores, através de exemplos envolvendo análise estrutural e/ou térmica;

•

Desenvolver e implementar rotinas que solucionem o problema OM para ambos os modelos (real e aproximado) a fim de obter os pontos de Pareto;

•

Acoplar aos algoritmos de otimização multiobjetivo, as rotinas de otimização estrutural;

•

Realizar um estudo comparativo através exemplos envolvendo otimização multiobjetivo;

3   

•

Implementar rotina para o cálculo dos parâmetros estatísticos da resposta de interesse, além de seus gradientes, através do método de MC, com diferentes técnicas de amostragem;

•

Implementar o procedimento para o desenvolvimento e aplicação do PCM para variáveis aleatórias considerando as distribuições de probabilidade;

•

Incorporar o PCM, na rotina para o cálculo dos parâmetros estatísticos e seus gradientes;

•

Realizar um estudo comparativo entre o MC e o PCM;

•

Adaptar o MBR para considerar as variáveis aleatórias no espaço da base reduzida;

•

Considerar os parâmetros estatísticos no processo de otimização uni e multiobjetivo, tanto utilizando apenas o MEF como via MBR;

•

Obter resultados finais, considerando as diferentes metodologias para a otimização multiobjetivo, para a análise, e para o cálculo dos parâmetros estatísticos.

4   

Início

Caso Estocástico: Definir variáveis aleatórias e suas distribuições

Definir: • Geometria do problema • Condições de contorno • Variáveis de projeto • Funções objetivo e restrição

Caso MBR: Definir amostras 

Caso MBR: Gerar soluções via MEF e a estruturas de dados off-line

Geração do projeto inicial  

  Análise Estrutural/Térmica  

   

Análise de Sensibilidade

Otimização Escalar

 

 

MEF

Determinístico

ou

ou

MBR

Estocástico

Otimizador SQP

ou Otimização Multiobjetivo (neste caso, repetir o processo para cada ponto de Pareto)

 

Geração do novo projeto

Não Convergiu?

Sim Fim

Figura 1.1 Algoritmo para as diferentes opções de otimização consideradas.

5   

1.3 ORGANIZAÇÃO DA DISSERTAÇÃO Esta dissertação consiste de seis capítulos organizados conforme será descrito a seguir. Após este breve capítulo introdutório, é apresentado no segundo capítulo a formulação do problema bidimensional térmico e estrutural, o MEF aplicado a estes, a separabilidade das variáveis do problema e a aproximação pelo MBR. Este capítulo é uma continuação da dissertação de (ALBUQUERQUE, 2005), com a inclusão do uso do elemento retangular bilinear (Q4) para o cálculo via MEF (o elemento triangular linear (CST) já havia sido considerado). No terceiro capítulo, é explanado o problema de otimização escalar, seus procedimentos de solução e as diferentes técnicas para a análise de sensibilidade, aplicadas via o MEF e o MBR. Em seguida são confrontados os resultados obtidos por ambas as abordagens. No capítulo quarto é apresentada a formulação do problema de otimização multiobjetivo, o conceito de pontos de Pareto e os métodos para encontrá-los. Alguns exemplos são resolvidos para a comparação das técnicas. Este capítulo é baseado na dissertação de (MACEDO, 2002) aplicado para otimização de treliças planas. As extensões aqui conduzidas, foram a incorporação de ferramentas para a solução de problemas 2D contínuos, a inclusão do método NNC (MESSAC et al., 2003) e de uma modificação da estratégia NBI (DAS e DENNIS, 1996), proposta neste trabalho, para problemas com mais de duas funções objetivo. O capítulo cinco é onde as incertezas passam a ser consideradas no problema de otimização. Nele é feita uma breve introdução sobre probabilidade e conceitos estatísticos para, a seguir, serem apresentados e investigados os métodos MC e o PCM para o cálculo das estatísticas de interesse. Posteriormente, é feita uma introdução à otimização robusta, para só então, serem apresentados os exemplos utilizando este conceito de otimização. No sexto capítulo serão feitas algumas conclusões e sugestões de trabalhos futuros. 1.4 REFERÊNCIAS AFONSO, S.M.B, LYRA, P.R.M, ALBUQUERQUE, T.M. M., R. S., MOTTA. “Structural Analysis and Optimization in the Framework of Reduced-Basis Method”. Structural and Multidisciplinary Optimization, Springer Berlin / Heidelberg, DOI: 10.1007/s00158-008-0350-4, Published online, January 2009. AFONSO, S. M. B. e PATERA, A. T. “Structural Optimization in Framework of Reduced Basis Method”. in XXIV CILAMCE, Congresso Ibero Latino Americano de Métodos Computacionais na Engenharia, 2003. ARORA J. S.; MESSAC, A.; MULLUR, A. A. “OPTIMIZATION OF STRUCTURAL AND MECHANICAL SYSTEMS. Chapter 4 - Multiobjective Optimization: Concepts and Methods”. Jasbir S Arora, University of Iowa, USA, 2007. 6   

ALBUQUERQUE T. M. M.. “Análise e Otimização de Problemas Térmicos e Estruturais Bidimensionais Através do Método das Bases Reduzidas”. Dissertação de mestrado, Dept. de Engenharia Civil, Universidade Federal de Pernambuco, Recife-PE, Brasil, 2005. COOK, R. D., MALKUS, D. S., PLESHA, M. E. e WITT R. J. “Concepts and Applications of Finite Element Analysis”. Fourth edition, Wiley, 2003. DAS, I. e DENNIS, J.E. “Normal Boundary Intersection: A New Method for Generating Pareto Surface in Nonlinear Multicriteria Optimization Problems”. SIAM J. Otimization, Vol. 8, No. 3, pp. 631-657, 1996. MACEDO, C. M. H. “Otimização de Treliças Planas sob Várias Solicitações com Ênfase a Problemas Multiobjetivos”. Dissertação de Mestrado, Dept. de Engenharia Civil, Universidade Federal de Pernambuco. Recife-PE, Brasil, 2002. MESSAC, A., ISMAIL-YAHAYA A. e MATTSON C. A. “The Normalized Normal Constraint Method for Generating the Pareto Frontier”. Structural Optimization, Vol. 25, No. 2, pp. 86-98, 2003. RAMAMURTHY, D. “Smart simulation techniques for the evaluation of parametric uncertainties in black box systems”. Thesis (M.S. in computer science). Washington State University. 2005. VANDERPLAATS, G. N. “Numerical Optimization Techniques for Engineering Design - with Applications”. 4rt Edition, Vanderplaats Research & Development, Inc., Colorado Springs, CO, 2004.

7   

CAPÍTULO 2 PROBLEMAS BIDIMENSIONAIS 2 PROBLEMAS BIDIMENSIONAIS 2.1 FORMULAÇÃO DE PROBLEMAS TÉRMICOS Nesta seção estão apresentadas as equações governantes dos problemas térmicos no regime estacionário (ou permanente). Inicialmente deve-se levar em conta que as equações governantes são descritas para um material homogêneo e ortotrópico, onde este corpo Ω possui contorno Γ. A variável T, que representa neste caso a temperatura, satisfaz a equação diferencial parcial nd

(2.1) ∑ i =1

∂qi = Q − α Ω (T − T0 ) em Ω ∂xi

(2.1)

onde qi representa a componente i do vetor de fluxo de calor, nd é o número de dimensões do problema, Q representa o termo de fonte (ou sumidouro) de calor no domínio, α Ω é o coeficiente de convecção no domínio, T0 é a temperatura ambiente e Ω representa o domínio do problema. O fluxo de calor se relaciona com o gradiente de temperatura através da lei de Fourier da condução, nd

(2.2) qi = ∑ −kij j =1

∂T ∂x j

(2.2)

onde kij representa os termos do tensor da condutividade térmica. A Equação (2.1) pode então ser expressa como nd

(2.3) −∑ i =1

∂ ∂xi

⎛ nd ∂T ⎞ ⎜⎜ ∑ kij ⎟⎟ + α Ω (T − T0 ) = Q em Ω ⎝ j =1 ∂x j ⎠

(2.3)

Além disso, as seguintes condições de contorno são consideradas

(2.4)

T =T

em ΓT

q n = q n − α Γ (T − T0

)

em Γ q

(2.4)

8

onde ΓT representa a porção de contorno com temperaturas prescritas T (condição de contorno de Dirichlet), Γ q representa a parcela do contorno sujeita a fluxo e/ou a convecção (condição de contorno de Cauchy), q n é o fluxo na direção n ortogonal à Γ q ,

α Γ é o coeficiente de convecção no contorno e T0 é a temperatura ambiente. As expressões foram escritas até então, utilizando notação indicial. Na forma matricial, as equações governantes no domínio Ω (2.3) e no contorno Γ q (2.4) podem ser expressas como

(2.5)

−∇T k∇T + α ΩT = α ΩT0 + Q em Ω k∇Tn + α ΓT = α ΓT0 + q n

(2.5)

em Γ eq

Para o caso bidimensional ortotrópico, ∇ = [ ∂ /∂x, ∂ /∂y ] e k é a matriz constitutiva T

que contém as condutividades térmicas do material nas direções x e y, dado por ⎡kx (2.6) k = ⎢ ⎣0

0⎤ k y ⎥⎦

(2.6)

2.1.1 Método dos Elementos Finitos aplicado à problemas térmicos

Na aproximação via elementos finitos é feita uma discretização do domínio em subdomínios chamados de elementos, onde serão efetuadas as aproximações. Para os problemas térmicos em regime estacionário, onde todos os graus de liberdade estão associados à temperatura em determinado nó, a temperatura Te(ξ,η) no elemento “e”, em qualquer ponto (ξ,η) deste elemento, é aproximada da seguinte forma (2.7) T e (ξ ,η )

n

∑ N (ξ ,η ) T i =1

e i

e

i

= N e Te

(2.7)

onde Ti e é a temperatura no nó i e N ie é a função de forma associado ao nó i. Considerando (ξ j ,η j ) as coordenadas locais do nó j as funções de forma são definidas de tal

forma que (ZIEKIEWICZ e TAYLOR, 2000) ⎧1, para j = i (2.8) N i (ξ j ,η j ) = ⎨ ⎩0, para j ≠ i

i, j = 1...n

(2.8)

9

Definimos, então, os vetores N e e Te , onde N e = ⎡⎣ N1e ,… , N ie ,… , N ne ⎤⎦ é o vetor das T

funções de forma do elemento “e”, Te = ⎡⎣T1e ,… , Ti e ,… , Tne ⎤⎦ é o vetor das temperaturas nodais deste mesmo elemento e n é o número de nós do elemento. Utilizando a aproximação dada em (2.7), as equações governantes (2.5) para o domínio Ω e para o contorno Γ q , podem ser aproximadas como

( −∇ k∇N + α N ) T = α (2.9) T

Ω

T + Q em Ω

Ω 0

( k∇Nn + α Γ N ) T = α ΓT0 + qn

em Γ q

(2.9)

onde N é o vetor contendo as funções de formas relacionadas aos nós, e T é o vetor contendo as temperaturas nodais dos nós. Para a solução do sistema de equações diferenciais (2.9), será usado o Método dos Resíduos Ponderados. A partir da parametrização da aproximação da solução do problema, este método consiste em encontrar os parâmetros referentes à aproximação considerada (no nosso caso T), que anule a integral em todo o domínio, do resíduo das equações ponderado por diferentes funções peso Wl , para l =1..nl, onde nl é o número de funções de peso utilizadas. As equações do problema estacionário de condução de calor devem satisfazer, então, a seguinte equação ⎡ ⎤ ⎢ ∫ ( −∇T k∇N + α Ω N )Wl d Ω + ∫ ( k∇Nn + α Γ N ) Wl d Γ q ⎥ × T = ⎢ ⎥ Γq ⎦ (2.10) ⎣ Ω

∫ (α

Ω

T0 + Q ) Wl d Ω +

Ω

∫ (α

Γ

T0 + q n ) Wl d Γ q

(2.10)

Γq

para l =1..nl Considerando uma discretização em elementos finitos do domínio em estudo, aplicando a aproximação de Galerkin do Método dos Resíduos Ponderados, onde as funções de peso são as próprias funções de forma (utilizadas para a aproximação da solução), e integrando por partes o termo com derivada de maior ordem, chega-se a forma fraca do problema estacionário de condução de calor (ALMEIDA, 2001) ⎡ ⎤ ⎢ ∫ ( ∇NT k∇N + α Ω NT N ) d Ω + ∫ α Γ NT Nd Γ q ⎥ × T = ⎢ ⎥ Γq ⎦ (2.11) ⎣ Ω

∫ (α

Ω

T0 + Q ) N d Ω + T

Ω

∫ (α

T0 + q n ) N d Γ q

(2.11)

T

Γ

Γq

As integrais que aparecem na Equação (2.11), podem ser avaliadas, somando as contribuições de cada elemento. A equação governante do problema térmico na forma discreta, pode ser expressa como 10

(2.12) K T T = FT

(2.12)

onde KT é a matriz de “rigidez térmica”, FT é o vetor com termos de cargas térmicas total e T é o vetor das temperaturas nodais incógnitos. O vetor FT é dado por (2.13) FT = FΩ + FΓ

(2.13)

onde FΩ representa as cargas no domínio e FΓ as cargas no contorno, dado por (2.14) FΩ = ∫ (α Ω T0 + Q ) NT d Ω ;

FΓ =

Ω

∫ (α T Γ

0

+ q n ) NT dΓ q

(2.14)

Γq

Para o problema descrito na Equação (2.12) a matriz de “rigidez térmica” apresenta a contribuição das parcelas relativas à condutividade térmica Kk, convecção no domínio K Ωc e convecção contorno K Γc , i.e. (2.15) K T = K k + K Ωc + K Γc

(2.15)

onde (2.16) K k = ∫ ∇NT k∇Nd Ω; Ω

K Ωc = ∫ α Ω NT Nd Ω; Ω

K Γc = ∫ α Γ NT Nd Γ q

(2.16)

Γq

2.2 FORMULAÇÕES DOS PROBLEMAS ELÁSTICOS Nesta seção serão apresentadas as equações governantes da elasticidade linear nas formas forte. Para tornar mais prática a apresentação da formulação, considera-se a deformação de um corpo homogêneo Ω com contorno Γ sujeito à forças volumétricas e carregamentos distribuídos. Além disso, considera-se que o gradiente do deslocamento é suficientemente pequeno, de forma que o modelo de elasticidade linear descreva adequadamente a deformação. Desta forma as equações de equilíbrio são satisfeitas na forma nd ⎛ ∂σ ij (2.17) ∑ ⎜ ⎜ i =1 ⎝ ∂x j

⎞ ⎟⎟ + b j = 0 em Ω ⎠

(2.17)

onde bj é a força volumétrica agindo sobre o domínio Ω na direção j e o tensor de tensões σ ij está relacionado linearmente com o vetor de deformação. Na forma matricial, as equações de equilíbrio para o caso bidimensional, podem ser expressas da seguinte forma

11

(2.18) ∇ σ + b = 0 T

(2.18)

onde ∇ é um operador diferencial, σ é o vetor das tensões e b o vetor das forças de volume, definidos respectivamente por ⎡∂ ⎢ ∂x (2.19) ∇T = ⎢ ⎢ ⎢0 ⎣

0 ∂ ∂y

∂ ∂y 0

⎤ 0⎥ ⎥, ∂⎥ ∂x ⎥⎦

⎡σ xx ⎤ ⎢σ ⎥ σ = ⎢⎢σ yy ⎥⎥ , xy ⎢ ⎥ ⎢⎣σ yx ⎥⎦

⎡ bx ⎤ b=⎢ ⎥ ⎣ by ⎦

(2.19)

A equação constitutiva, para a elasticidade linear na forma matricial, é apresentada através da seguinte relação (2.20) σ = Dε

(2.20)

onde ε o vetor de deformações dado por ε = ∇u onde u é o vetor dos deslocamentos, que no caso de problemas de elasticidade plana é dado por

⎡u ( x, y ) ⎤ (2.21) u = ⎢ ⎥ ⎣ v ( x, y ) ⎦

(2.21)

Para materiais isotrópicos com módulo de elasticidade E e coeficiente de Poisson υ , a matriz de elasticidade D (Equação (2.20)) para o estado plano de deformação é dada por (VILLAÇA e GARCIA, 2000). ⎡ 2c2 + c1 ⎢ c1 (2.22) D = ⎢ ⎢ 0 ⎢ ⎢⎣ 0

c1 2c2 + c1 0 0

0 0 c2 c2

0⎤ ⎥ 0⎥ c2 ⎥ ⎥ c2 ⎥⎦

(2.22)

onde, c1 e c2 são as constantes de Lamé definidas por (2.23) c1 =

Eυ E , c2 = (1 + υ )(1 − 2υ ) 2(1 + υ )

(2.23)

A equação de equilíbrio para elasticidade plana pode então ser expressa, em função dos deslocamentos, da seguinte forma 12

T (2.24) ∇ D∇u + b = 0

(2.24)

Além das equações no domínio, têm-se as condições de contorno de deslocamento em Γ D e de tensão prescrita em Γ N , ou seja,

(2.25)

u = u em Γ D

(2.25)

σ = σ em Γ N

2.2.1 Método dos Elementos Finitos aplicado à elasticidade bidimensional Na aproximação via elementos finitos para os problemas bi-dimensionais onde todos os graus de liberdade estão associados aos deslocamentos nodais. Os deslocamentos ue(ξ,η) e ve (ξ,η) em qualquer ponto (ξ,η) no elemento “e”, onde u e v são, respectivamente, as componentes horizontais e verticais do deslocamento, são aproximados da seguinte forma u e ( ξ ,η )

(2.26)

n

∑ N (ξ ,η ) u i =1

v (ξ ,η ) e

n

e i

e i

(2.26)

∑ N (ξ ,η ) v i =1

e i

e i

onde uie e vie são as componentes do deslocamento no nó i do elemento e e N i é a função de forma associado ao nó i, da mesma forma como foi definido para o caso térmico na seção 2.1.1. Com isso pode-se reescrever a equação que relaciona deformação e deslocamento para um elemento genérico e da seguinte maneira (2.27) ε e

n

B ed e = ∑ Biedie

(2.27)

i=1

onde Be é matriz que relaciona as deformações com deslocamentos, die = [uie

T

vie ] e n é

o número de nós por elemento. Para um nó específico i, Bie pode ser escrita como (2.28) Bie = ∇N ie

(2.28)

Partindo da formulação forte e aplicando a aproximação de Galerkin em Elementos Finitos (ZIEKIEWICZ e TAYLOR, 2000) a equação diferencial governante para problemas de elasticidade pode ser escrita como

13

⎡ ⎤ (2.29) ⎢ ∫ BT DBd Ω ⎥ u = ∫ NT bd Ω + ∫ NT f n d Γ + ∫ NT ft d Γ ⎣Ω ⎦ Ω Γn Γt

(2.29)

Na Equação (2.29) b é o vetor das forças de volume agindo sobre o domínio Ω , f n e ft são respectivamente os vetores das forças normais e tangentes aplicadas no contorno Γ n e Γt respectivamente, N é a matriz relacionada com as funções de forma utilizadas pelo Método dos Elementos Finitos em modelos bi-dimensional e B a matriz que relaciona as componentes das deformações com deslocamentos, apresentadas anteriormente. A Equação (2.29), pode ser reescrita como (2.30) K s u = Fs

(2.30)

onde Ks é a matriz de rigidez da estrutura e Fs é o vetor com termos de carregamento e condições de contorno.

2.2.2 Acoplamento termo-elástico A equação governante para problemas termo-elásticos constitui uma extensão da apresentada na Equação (2.29) com a inclusão de um termo de carregamento adicional relacionado com as deformações térmicas. Assim, tem-se ⎡ ⎤ (2.31) ⎢ ∫ BT DBd Ω ⎥ u = ∫ NT bd Ω + ∫ NT f n d Γ + ∫ NT ft d Γ + ∫ BT Dε 0 d Ω Ω Γ Γ Ω ⎣Ω ⎦

(2.31)

Na equação anterior, ε 0 é o vetor das deformações iniciais devido à variação térmica. Assim, o último termo da Equação (2.31) é o responsável pelo acoplamento entre os problemas térmico e estrutural. As deformações iniciais são causadas por variação térmica do material em uma estrutura e devido a restrições oriundas das condições de apoio. Considerando um material isotrópico com coeficiente de expansão térmica γ, coeficiente de Poisson υ e uma variação da temperatura em relação a uma temperatura de referência ΔT, tem-se para as condições de Estado Plano de Deformação que

⎡(1 + ν ) γ ΔT ⎤ ⎢ ⎥ 1 + ν ) γ ΔT ⎥ (2.32) ε 0 = ⎢( ⎢ ⎥ 0 ⎢ ⎥ 0 ⎣ ⎦

(2.32)

14

As deformações térmicas produzem apenas uma variação no volume da estrutura, sem variar sua forma. Por isso, as componentes da deformação de cisalhamento são iguais à zero. Como se pode observar, o efeito dos carregamentos térmicos é introduzido na análise estrutural transferindo-se o campo de temperaturas do modelo térmico para o modelo elástico. Este campo de temperaturas é utilizado para computar as deformações iniciais térmicas, que introduzem carregamentos mecânicos na estruturas através do último termo da Equação (2.31).

2.3 MAPEAMENTO DAS EQUAÇÕES DO MEF Os elementos utilizados neste trabalho serão: o quadrilateral isoparamétrico bilinear de quatro nós (Q4) e o triângular linear de três nós (T3) conhecido como CST – “Constant Strain Triangle”, pois como suas funções de forma são lineares, a deformação é constante no elemento. A fim de evitar a necessidade de uma nova discretização, bem como reduzir o custo computacional relativo ao processamento da solução via MEF a cada mudança de parâmetro que altere a geometria da estrutura, pode-se realizar os cálculos sobre um domínio de referência Ω mantido constante durante todo o procedimento e aplicar transformações a fim de obter as equações em outro domínio Ω *, chamado domínio real. Para tal utiliza-se inicialmente o conceito de separabilidade (ALBUQUERQUE, 2005), escrevendo os termos da matriz de rigidez e vetor de carga em forma decomposta. Assim, os domínios bidimensionais serão divididos em várias regiões r = 1… R , de acordo com a dependência dos parâmetros variáveis para o problema em particular. Considerando a contribuição de cada sub-região Ω r a matriz de rigidez e o vetor dos termos de carregamento definidos nas Equações (2.12) e (2.30) podem ser reescritos como R

(2.33) K = ∑ K r r =1

R

F = ∑ Fr

(2.33)

r =1

Separando os termos (regiões) que dependem do parâmetro variável, um novo conjunto de parâmetros constituirá somente um novo dado a ser ajustado na forma decomposta, não se tornando necessário outro processo de discretização de um novo domínio, por conseguinte, contribuindo para a melhoria da eficiência computacional. Este procedimento se mostrará ainda mais interessante mais adiante, ao se utilizar o MBR (Método da Base Reduzida) (AFONSO e PATERA, 2003). Cada novo domínio real Ω* bidimensional é dividido em um número de regiões R, onde cada região r possui uma transformação diferente, entre o domínio real Ω*r e o domínio de referência Ω r , a depender dos parâmetros ajustáveis de cada região. Isto significa que para cada região, há uma associação linear entre os domínios de referência e o real, definido da seguinte forma (2.34) xi*r = G ijr x rj

ou x*r = G r xr

com i, j = 1, 2, e r = 1,..., R

(2.34)

15

onde x1 e x2 correspondem às coordenadas nas direções x e y respectivamente, e a matriz de mapeamento (ou transformação) G é dada por ⎡∂ x1* ⎢ ∂x r (2.35) G (μ) = ⎢ 1* ⎢∂ x2 ⎢ ⎣∂ x1

∂ x1* ⎤ ⎥ ∂ x2 ⎥ ∂ x2* ⎥ ⎥ ∂ x2 ⎦

(2.35)

Esta matriz jacobiana (2.35) é escrita em função dos parâmetros variáveis μ , para maiores informações vide (ALBUQUERQUE, 2005). O mapeamento desta transformação permitirá obter as equações para novos domínios, sem a necessidade de gerar nova malha nem construir as matrizes ou vetores do novo sistema de equação. Considerando a transformação apresentada, pode-se então obter as equações governantes para um novo domínio qualquer Ω*r em função dos parâmetros μ , a partir das equações no domínio de referência Ω r . Para uma região específica r, as transformações dos infinitesimais do domínio e do contorno podem ser escritas como (2.36) d Ω*r = ϕΩr (μ)d Ω r e d Γ*r = ϕΓr (μ)d Γ r

(2.36)

onde (2.37) ϕΩr (μ) = det[G r (μ)]−1 e ϕΓr (μ) = [G r (μ)]−1 n r t onde n*r t é o vetor unitário tangente ao contorno Γr de uma região particular r e

(2.37) sig-

nifica a norma euclidiana do vetor.

2.3.1 Mapeamento da equação térmica Aplicando os conceitos de separabilidade, juntamente com as transformações aplicada á matriz ∇N e aos infinitesimais d Ω e d Γ , tem-se que as matrizes de “rigidez térmica” e o vetor de carregamento térmico decompostos, podem ser transformados do domínio de referência para o domínio real como K *kr ( μ ) =

∫ ∇N

Ωr

(2.38) K *Ωrc ( μ ) =

∫α

Ω

r Ω

T

⎡⎣G rT (μ)k r G r (μ) ⎤⎦ ∇N det[G r (μ)]−1 d Ω r

NT N det[G r (μ )]−1 d Ω r

(2.38)

r

K *Γrc ( μ ) = ∫ α Γr NT N [G r (μ)]−1 n r t d Γ rq Γ rq

16

∫ (α

FΩ*r = (2.39)

r Ω

T0 + Q r ) NT det[G r (μ )]−1 d Ω r

Ωr

FΓ*r =

∫ (α T r Γ

0

− q rn ) NT [G r (μ)]−1 n r t d Γ rq

(2.39)

Γ rq

Na presente aplicação, o único termo dependente de μ é G r (μ) . Estas transformações, nas aplicações aqui analisadas, são constantes no domínio de cada região, sendo possível “retirá-las” das integrais e desta forma reescrever as matrizes e vetores decompostos apresentados nas Equações (2.38) e (2.39) com o auxilio da Equação (2.37), da seguinte maneira nt

(2.40) K *kr ( μ ) = ∑ βTr ( μ ) j K rk j ; K *Ωrc ( μ ) = ϕΩr (μ)K Ωr c ; K *Γrc ( μ ) = ϕΓr (μ)K Γr c

(2.40)

(2.41) FT*r ( μ ) = ϕΩr (μ)FΩr + ϕΓr (μ)FΓr

(2.41)

j =1

onde os termos independentes do parâmetro µ são dados pelas Equações (2.14) e (2.16) aplicada a cada região. As matrizes K rk j são formadas a partir da sub-decomposição da matriz K *kr de condutividade térmica. Nesta decomposição nt e βTr ( μ ) j variam de acordo com a transformação específica da região. As matrizes K rk j são calculadas da seguinte forma (2.42) K rk j =

∫ ∇N k ∇Nd Ω T

Ω

r j

r

,

j = 1...nt

(2.42)

r

onde as sub-matrizes constitutivas k rj são obtidas através da decomposição de k *r (μ) , que é definida a partir do seu valor de referência como (2.43) k *r (μ) = ⎡⎣G rT (μ)k r G r (μ) ⎤⎦ det[G r (μ )]−1

(2.43)

A partir desta transformação (Equação (2.43)) são retiradas não só as matrizes k rj como também o parâmetro β kr ( μ ) j que contêm os fatores de dependência da variável µ, de tal modo que nt

(2.44) k *r (μ) = ∑ βTr ( μ ) j k rj

(2.44)

j =1

Para maiores detalhes vide (ALBUQUERQUE, 2005).

17

A partir da Equação (2.33) e das matrizes e vetores decompostos mostrados nas Equações (2.40) e (2.41), respectivamente, obtém-se as equações governantes, no domínio real, para valores de μ quaisquer, da mesma forma como foi feito nas Equações (2.15) e (2.16). Temos assim (2.45) K *T T* = FT*

(2.45)

Resolve-se este novo sistema para se obter as temperaturas atualizadas T* .

2.3.2 Mapeamento da equação elástica Utilizando as transformações de coordenadas das Equações (2.34) e (2.35), aplicadas à matriz B, as relações infinitesimais da Equação (2.36) e o conceito da separabilidade, tem-se as novas sub-matrizes de rigidez K *er e os termos de carregamento Fe*r no domínio real de cada região (Equações (2.29) e (2.30)), da seguinte forma (2.46) K *s r ( μ ) =

∫ (G ( μ ) B rT

T

)Dr (G r ( μ ) B) det[G r (μ)]−1 d Ω r

(2.46)

Ωr

Fs*r ( μ ) =

∫N b T

r

det[G r (μ)]−1 d Ω r +

Ωr

(2.47)

∫

+

∫Nf T

r n

[G r (μ)]−1 n nr t d Γ nr

Γ nr

(2.47)

NT ftr [G r (μ)]−1 ntr t d Γtr

Γtr

Assim como para o caso térmico, as componentes que dependem de μ podem ser retirados das integrais e desta forma obter nt

(2.48) K *s r ( μ ) = ∑ β sr ( μ ) j K rj

(2.48)

(2.49) Fs*r ( μ ) = ϕΩr (μ)Fbr + ϕΓr (μ)F fr

(2.49)

j =1

onde

∫B

(2.50) K rj =

T

Drj BdV r

(2.50)

Vr

(2.51) Fbr =

∫ N b dV T

V

r

r

r

, F fr =

∫Nf T

Γ

nr

r n

d Γr +

∫N f T

Γ

r t

d Γr

(2.51)

tr

18

Nas equações anteriores, cada componente das sub-matrizes de elasticidade Drj e

β sr ( μ ) j são definidos de maneira análoga à decomposição da matriz de condutividade térmica k *r (Equação (2.43) e (2.44)), de tal modo que nt

(2.52) D*r (μ) = ∑ β sr ( μ ) j Drj = ⎡⎣G rT (μ) Dr G r (μ) ⎤⎦ det[G r (μ)]−1

(2.52)

j =1

Como se pode verificar na Equação (2.52), dependendo da transformação a ser aplicada em cada região, a pseudo-matriz de elasticidade D*r depende diretamente de Gr e, portanto, da geometria/transformação da região. O vetor que acopla os efeitos térmicos ao carregamento elástico é calculado novamente para cada variação do parâmetro μ (estagio on-line do MBR que será tratado mais adiante), pois os valores das deformações variam em função da temperatura, que não pode ser mapeada no domínio da região r. Desta forma, para cada valor de μ se tem *r (2.53) FTE (μ ) =

∫B

*T

Dε*o (μ)dV *r

(2.53)

V *r

O vetor carregamento que acopla os efeitos térmicos, pode ser decomposto de forma linear em relação às temperaturas, porém tal decomposição não foi aplicada. O vetor de carregamento total considerando os efeitos térmicos nas deformações iniciais de uma região r para um valor de μ qualquer, pode ser calculado como *r (2.54) Fs*r ( μ ) = ϕΩr (μ)Fbr + ϕΓr (μ)F fr + FTE (μ )

(2.54)

A partir da Equação (2.33) as equações governantes são obtidas para valores de μ quaisquer. Para se obter os deslocamentos atualizados u* , no domínio real, resolve-se o novo sistema. (2.55) K *s u* = Fs*

(2.55)

2.4 MÉTODO DA BASE REDUZIDA O MBR consiste numa projeção do tipo Galerkin em um espaço de ordem reduzida que contém soluções (base) para o problema de interesse em pontos selecionados do espaço de projeto. O campo de soluções u(μ) (onde, no presente trabalho, u pode representar deslocamentos ou temperaturas) faz parte do espaço de solução de dimensões infinitas Y , e claramente os possíveis valores de u(μ) não cobrem inteiramente o espaço Y . Considerando Y como um espaço tri-dimensional, observa-se que u em função de μ pode ser representado por uma curva ou superfície conforme a Figura 2.1 (a). A partir disto, po19

de-se afirmar que o espaço das soluções não se trata de um espaço de dimensão infinita N Y e sim um espaço de dimensão reduzida Y , onde é observado a dependência do parâmetro variável μ (PRUD` HOMME et al, 2002), como ilustrado na Figura 2.1 (b).

u(μnovo) u(μ) u(μ1)

u(μN) u(μ2) Y

(a)

(b)

Figura 2.1 Redução de dimensões do espaço solução - (a) Variação da solução com o parâmetro μ; e (b) aproximação da solução μnovo através de uma combinação linear de soluções u(μi) pré-calculadas. Conforme comentado acima, não é mais preciso representar todas as possíveis funções em Y para aproximar u(μ), tornando-se necessário aproximar somente aquelas funções no espaço reduzido YN. Com isso, pode-se simplesmente calcular a solução u em vários pontos do conjunto correspondente a diferentes valores de μ e para qualquer novo parâmetro μnovo pode-se aproximar u(μnovo) através de uma combinação linear de soluções conhecidas u(μi). A ilustração desta combinação linear encontra-se na Figura 2.1 (b). Para se construir uma aproximação para o campo de solução, inicialmente deve-se construir a amostra de pontos sob a qual a base será montada. Desta forma tem-se a definição do conjunto de amostras representada da seguinte forma

{

(2.56) S N = (μ1 ,…, μ R ) ,…, (μ1 ,…, μ R ) 1

N

}

(2.56)

onde cada (μ 1 , … , μ R ) está contido no espaço Dμ, isto significa que i

(2.57) μ low ≤ μ ri ≤ μ up

(2.57)

20

no qual μ low e μ up são os limites inferior e superior, respectivamente, do espaço de projeto Dμ e N é o número de amostras. Para cada componente de SN calcula-se a solução u(μ) através do Método dos Elementos Finitos, utilizando as Equações (2.12) e/ou (2.30). Isto nos permite construir um espaço de Base Reduzida WN tal que (2.58) W N = span{u[ (μ1 ,…, μ R ) ],…,u[ (μ1 , …, μ R ) 1

N

]}

(2.58)

Para simplificar a notação defini-se ζ i ∈ R n como

(

(2.59) ζ i = u (μ1 ,..., μ R )

i

)

(2.59)

e assim WN pode ser escrito como (2.60) W N = span{ζ i } com i = 1,…, N

(2.60)

A equação acima significa que qualquer vetor do espaço reduzido WN pode ser expresso como uma combinação linear das soluções ζ i . Vale salientar que as soluções nos pontos amostrais devem ser linearmente independentes, para que ζ seja base de WN. Caso haja soluções redundantes elas devem ser desconsideradas. Isto nos leva à definição da aproximação do MBR para a solução como N

(2.61) u N (μ) = ∑ α ( μ ) ζ j j

α j ∈ Rn

(2.61)

j =1

ou na forma matricial (2.62) u N (μ) = Zα (μ )

(2.62)

no qual (2.63) Z = [ζ1 ,..., ζ N ] ; αT ( μ ) = [α ( μ ) ,..., α ( μ ) ] 1

N

(2.63)

Partindo da equação fundamental do MEF (Equações (2.12) e (2.30)) e substituindo o vetor solução u (ou T) por uN tem-se (2.64) K (μ )Zα (μ ) = F (μ )

(2.64)

Pre-multiplicando ambos os lados de (2.64) por ZT , ou alternativamente, aplicando mínimos quadrados em (2.62), obtemos 21

(2.65) Z T K (μ )Zα(μ ) = Z T F(μ )

(2.65)

Com isso, os coeficientes lineares α(μ) são obtidos resolvendo o problema no espaço WN dado da seguinte forma (2.66) K N (μ)α (μ) = F N (μ)

(2.66)

no qual (2.67) K N (μ ) = Z T K (μ)Z ∈ R N × N , F N (μ) = Z T F(μ) ∈ R N

(2.67)

Para ganhar ainda mais agilidade, pode-se utilizar o conceito da separabilidade, já discutido, onde a matriz de rigidez e o vetor de carregamento no espaço WN podem ser escrito como R

R

(2.68) K = ∑ K , F = ∑ F Nr N

Nr

N

r =1

(2.68)

r =1

com K Nr e F Nr representando, respectivamente, a contribuição de cada sub-região na matriz de rigidez e no vetor de carregamento no espaço WN, i.e. (2.69) K N r (μ ) = Z T K r (μ)Z , F N r (μ) = Z T F r (μ)

(2.69)

Nas equações acima os termos K r (μ ) e F r (μ) podem ser mapeados como foi visto nas seções 2.3.1 e 2.3.2, e assim reescrevê-los, para se obter os mesmos correspondentes ao domínio real de cada região r . Por exemplo, para o caso elástico, resulta em nt

(2.70) K *s Nr ( μ ) = ∑ βsr ( μ ) j K sNr j ; Fs* Nr ( μ ) = ϕΩr (μ)FbNr + ϕΓr (μ )F fNr

(2.70)

j =1

r

onde apenas os termos βs j , ϕbr , ϕnr e ϕtr dependem do parâmetro variável, K sNr j e F jNr são constantes, calculados para o domínio de referência e projetados no espaço de base reduzida. Além disso, nt é o número de termos da matriz e tal como anteriormente mencionado varia com a transformação específica da região correspondente r. Vale salientar que estes parâmetros foram discutidos com mais profundidade na seção 2.3. Para o caso térmico, as matrizes e vetores de cada região podem ser escritos da seguinte forma nt

r Nr (2.71) K *k Nr ( μ ) = ∑ βTr ( μ ) j K kNr j ; K *ΩNrc ( μ ) = ϕΩr (μ )K ΩNrc ; K *ΓNr c ( μ ) = ϕ Γ (μ )K Γc

(2.71)

j =1

22

(2.72) F* Nr ( μ ) = ϕdr (μ)FdNr + ϕcr (μ)FcNr

(2.72)

2.4.1 Procedimento Computacional Off-Line/On-Line No mapeamento das equações governantes da base reduzida, os termos que são dependentes/não dependentes do parâmetro μ são claramente distinguidos. Como uma conseqüência desta subdivisão, a implementação computacional para cálculo das grandezas pelo Método da Base Reduzida é conduzida por um algoritmo “off-line” (independente de μ)/ “on-line” (dependente de μ). A idéia deste algoritmo consiste no fato que a parte “off-line” só é executado uma vez, gerando variáveis contendo as matrizes K iNr e vetores de força FiNr . Consequentemente, no estágio on-line utiliza-se estes dados gerados anteriormente para executar uma resposta em tempo real para um novo μ. No caso de problemas de termo-elasticidade, o procedimento do Método da Base Reduzida é implementado de acordo com o algoritmo mencionado na Tabela 2.1. Para os problemas térmicos o algoritmo é bastante similar. Para mais detalhes sobre o procedimentos do MBR utilizado, bem como suas aplicações e resultados, vide (AFONSO e PATERA, 2003; ALBUQUERQUE, 2005; MOTTA et al., 2007; AFONSO et al., 2009) Tabela 2.1 Algoritmo do Método da Base Reduzida: OFF-LINE/ON-LINE. (Problemas acoplados) OFF-LINE - μ independent: 1- Problema Térmico 1.1 Escolher a amostra: S N =

{( μ ,..., μ ) ,..., ( μ ,..., μ ) } ; 1

1

R

N

1

R

1.2 Construir a matriz de soluções térmicas via MEF: ΖT = [ζ T1 ,..., ζ TN ] ; 1.3 Construir as componentes matrizes de “rigidez térmica” no espaço reduzido: K kNr j = ZTT K rk j ZT ; K ΩNrc = ZTT K Ωr c ZT ; K ΓNrc = ZTT K Γr c ZT ;

1.4 Constrói as componentes dos vetores de cargas térmicas no espaço reduzido: FΩNr = ZTT FΩr ; FΓNr = ZTT FΓr . 2- Elástico Acoplado 2.1 Construir a matriz de soluções dos deslocamentos via MEF (a partir da amostra e das soluções térmicas): Ζ s = [ζ s1 ,..., ζ sN ] ; 2.3 Construir as componentes da matriz de rigidez no espaço reduzido:

K sNr j = ZTs K rs j Z s ;

2.4 Construir as componentes dos vetores de carga estrutural no espaço reduzido: FbNr = ZTs Fbr ; F fNr = ZTs F fr .

23

ON-LINE – para um novo vetor μ : 1- Problema Térmico 1.1 Formar a matriz de “rigidez térmica” no espaço reduzido: R nt ⎛ r ⎞ Nr Nr r *N K T ( μ ) = ∑ ⎜ ϕΩ (μ)K Ωc + ϕΓ (μ )K Γc + ∑ βTr ( μ ) j K kNr j ⎟ ; r =1 ⎝ j =1 ⎠ 1.2 Formar o vetor de cargas térmicas no espaço reduzido: R

FT* N ( μ ) = ∑ (ϕΩr (μ)FΩNr + ϕΓr (μ)FΓNr ) ; r =1

1.3 Resolver: K T* N (μ)αT (μ) = FT* N (μ) ; 1.4 Calcular: T N (μ) = ZT αT (μ) . 2- Elástico Acoplado R

nt

2.1. Formar a matriz de rigidez no espaço reduzido: K *s N (μ) = ∑∑ β sr (μ ) j K sNr j ; r =1 j =1

2.2 Formar o vetor de carga estrutural no espaço reduzido: R

* Nr Fs* N (μ) = ∑ ⎡⎣ϕΩr (μ)FΩNr + ϕΓr (μ)F fNr + FTE ( μ )⎤⎦ ; r =1

2.3 Resolver: K *s N ( μ)αs ( μ) = Fs* N ( μ) ; 2.4 Calcular: u N ( μ) = Z s αs ( μ) ;

2.5 EXEMPLOS 2.5.1 Placa quadrada com orifício central O primeiro exemplo a ser considerado será uma placa quadrada com um orifício central cujo comportamento estrutural pode ser modelado como estado plano de tensões. Devido à simetria do problema apenas um quarto da placa será modelada. A geometria, as condições de contorno e as variáveis de projeto estão apresentadas na Figura 2.1. As variáveis do problema são as dimensões μ1 e μ2 representadas na figura. 24

μ1 μ2

1 2

3

Figura 2.1 Um quarto da placa quadrada – Definição do problema. As propriedades do material e as dimensões do modelo são as seguintes: •

O Módulo de elasticidade das regiões 1 e 2 é E = 105 MPa,

•

O Módulo de elasticidade da região 3 é E3 = 5x104 MPa,

•

O Coeficiente de Poisson υ = 0.3,

•

Espessura t = 1mm,

•

Comprimento lateral é 100mm,

•

Carregamento distribuído p = 1 N/mm.

Para as dimensões do orifício os limites inferior e superior impostos são 25mm e 75mm, respectivamente. Neste exemplo, é analisada a estrutura citada para μ1 = 45 mm e μ2 = 55 mm. Um teste de convergência da malha foi realizado para dois tipos de elementos utilizados no modelo de Elementos Finitos (EF), variando seus tamanhos médios (Le). A convergência foi estudada em relação à flexibilidade total da estrutura (FTu). Os elementos utilizados foram o quadrilateral isoparamétrico bilinear de quatro nós (Q4) e o triângulo linear de três nós (T3). A Figura 2.2 apresenta um gráfico do valor da flexibilidade total para diferentes refinamentos de malha. Na Figura 2.2 (a) o eixo x é o tamanho dos elementos (Le), enquanto que na Figura 2.2 (b) o eixo x é o número de graus de liberdade. As malhas utilizadas são estruturadas e uniformes, os elementos quadrilaterais possuem lados iguais (quando possível) e as malhas de elementos triangulares se sobrepõem as malhas de elementos quadrilaterais, ou seja, possuem os mesmos nós. A densidade (refinamento) da malha é definida pelo tamanho médio dos elementos (Le).

25

a)

b)

Figura 2.2 Um quarto da placa quadrada – Convergência da malha de EF em relação à flexibilidade total (10-3 J): a) Le – tamanho médio do elemento, b) Número de graus de liberdade. Apesar da convergência mais eficiente do elemento quadrilateral, a formulação do sistema de equações do elemento triangular é mais rápida devido à obtenção da integração analítica das equações, sem a necessidade de utilizar a integração numérica. Para maiores detalhes vide (COOK et al, 2003). Agora, será examinado o procedimento de análise via MBR. Para tal, as malhas de elementos finitos consideradas, foram as que apresentaram um erro (em relação à flexibilidade utilizando a malha mais refinada) próximo de 0,1 %, portanto: •

Para o elemento Q4 adotou-se elementos de dimensão média de 2mm. O modelo possui 3.952 graus de liberdade (erro 0,12%).

•

Para o elemento T3 adotou-se elementos de dimensão média de 1mm. O modelo possui 15.402 graus de liberdade (erro 0,10%).

Os resultados para estas malhas estão indicados com um “X” na Figura 2.2 (a e b). O valor absoluto dos deslocamentos máximos para o modelo utilizando o elemento Q4 de 2mm e o elemento T3 de 1mm foram, respectivamente, 6.5819 x10-3 mm e 6.5828x10-3 mm, as flexibilidade foram, respectivamente, 0.42714 x10-3 J e 0.42722 x10-3 J e o tempo de CPU foram 4.41s e 9.83s, respectivamente. A Figura 2.3 apresenta a distribuição das tensões de von Mises do modelo deformado, utilizando as discretizações mencionadas, distribuições similares podem ser observadas. Na Figura 2.3 (a) foi utilizada uma malha de elementos Q4 de tamanho médio de 2mm e na Figura 2.3 (b) foi utilizada uma malha de elementos T3 de tamanho médio de 1mm.

26

a) Q4

b) T3

Figura 2.3 Um quarto da placa quadrada – Contornos das tensões von Mises para a configuração deformada calculadas considerando-se os elementos Q4 e T3. Para a aproximação via MBR, 3 regiões (indicadas na Figura 2.1) são definidas. A base reduzida foi construída no espaço viável das variáveis geométricas, D ={[25, 75]²}. O número de pontos amostrais N foi variado de 4 à 10, e as análises foram feitas considerando as malhas (modelos computacionais) utilizadas anteriormente. Para verificar a convergência do MBR, construímos o gráfico da flexibilidade aproximada (FNTα) com o número de pontos amostrais N, apresentado na Figura 2.4. No gráfico pode-se verificar a convergência para poucas soluções base (tamanho da amostra).

Figura 2.4 Um quarto da placa quadrada – Convergência do MBR com malhas de EF utilizando elementos Q4 e T3. A Tabela 2.2 ilustra os tempos de CPU do MBR “OFF-Line” e “ON-Line” (para uma reavaliação), para diferentes tamanhos de amostras, utilizando o elemento Q4. Fica evidenciada a eficácia no processo de reavaliação de funções, a ser utilizado para procedimentos de otimização e cálculos estocásticos, os quais serão considerados nesta dissertação.

27

Tabela 2.2 Desempenho do MBR utilizando o elemento Q4. N Tempo “OFF-Line” (s) Tempo “ON-Line” (s) 10 54.9937 0.0035 9 50.4739 0.0031 8 45.9483 0.0018 7 41.3479 0.0020 6 36.8836 0.0019 5 32.4563 0.0017 4 27.9801 0.0017

2.5.2 Exemplo termo-elástico acoplado Como segundo exemplo, será considerada uma estrutura plana bidimensional indicada na Figura 2.5 (a) e (b) com o total de treze (13) regiões, especificadas para o presente problema, numeradas conforme indicado na Figura 2.5 (a) em função dos parâmetros variantes do problema em questão (μ1 à μ4 indicados na Figura 2.5(a)). Todas as dimensões do problema estão especificadas na Figura 2.5 (a). Condições de contorno para o problema elástico e de transferência de calor estão indicadas na Figura 2.5 (b). A estrutura é engastada na lateral esquerda e esta isolada termicamente (não possui fluxo de calor) nos vazios no centro da estrutura e na parte superior.

a)

28

b) Figura 2.5 Definição do problema: a) Geometria, b) Condições de contorno. As propriedades do material são: E= 104 MPa (1000 kN/cm²), o coeficiente de Poisson υ = 0.2, coeficiente de convecção h = 10 W/m² ºC, a condutividade térmica κx = κy = 1 W/(mK), o coeficiente de dilatação térmica α = 2x10-4 ºC-1 e as condições de convecção são para temperatura ambiente Ta = 20 ºC na base inferior e 100 ºC na lateral direita da estrutura (o sentido da convecção indicado na figura é meramente ilustrativo). Os quatro parâmetros considerados são: as espessuras t1, t2 e as distâncias horizontais L1 e L2, indicados na Figura 2.5 (a). A aproximação foi construída de forma que os parâ2 metros μ=( t1, t2, L1, L2 ) pertençam ao espaço de projeto Dμ = [ 0,5;5] × [10;21] ×

[28;38] cm4. A discretização de um exemplo de domínio real, onde μ1 = 2, μ2 = 2, μ3 = 14, μ4 = 31 cm e o domínio de referência se encontram respectivamente na Figura 2.6 (a) e (b).

a)

b)

Figura 2.6 Malha de elementos finitos: a) domínio real (projeto inicial) e b) domínio computacional.

29

Para este exemplo de discretização foram usados elementos triangulares T3, totalizando 1976 graus de liberdade para o problema elástico e 1023 para o térmico. Um teste de convergência, similar ao realizado no exemplo anterior, foi conduzido para esta estrutura. As malhas foram refinadas, variando o tamanho médio (Le) dos elementos. Os elementos utilizados foram os mesmos do exemplo anterior: o quadrilateral (Q4) e o triangular (T3). A Figura 2.7 apresenta um gráfico do valor da flexibilidade total para diferentes refinamentos de malha. Na Figura 2.7 (a) o eixo da abscissa representa o tamanho dos elementos (Le), enquanto que na Figura 2.7 (a) o eixo da abscissa indica o número de graus de liberdade.

Figura 2.7 Problema acoplado – Convergência da malha de EF Foi então examinado o procedimento de análise via MBR. Para tal, as malhas de elementos finitos consideradas, foram as que apresentaram um erro (em relação à flexibilidade utilizando a malha mais refinada) próximo de 0,1 %, portanto: •

Para o elemento Q4 adotou-se elementos de dimensão média de 0.5 cm. O modelo possui 2.046 graus de liberdade (erro 0,11%).

•

Para o elemento T3 adotou-se elementos de dimensão média de 0.25 cm. O modelo possui 7.342 graus de liberdade (erro 0,10%).

A Figura 2.8 apresenta a distribuição das temperaturas e o contorno das tensões de von Mises no modelo deformado, utilizando as malhas descritas anteriormente. Na Figura 2.8 (a) e (b) são apresentados o resultado da análise para uma malha de elementos Q4 de tamanho médio de 0.5 cm e na Figura 2.8 (c) e (d) o resultado da análise para uma malha de elementos T3 de tamanho médio de 0.25cm.

a) Térmico – Q4 (Cº)

b) Acoplado – Q4 (kN/cm²) 30

c) Térmico – T3 (Cº)

d) Acoplado – T3 (kN/cm²)

Figura 2.8 Problema acoplado – Contornos das temperaturas e das tensões von Mises para as malhas consideradas. O valor absoluto dos deslocamentos máximos para o modelo utilizando o elemento Q4 e para o elemento T3 foram, respectivamente, 0.4141 cm e 0.4046 cm, já as flexibilidades foram, respectivamente, 4,6828 kJ e 4,6835 kJ enquanto que o tempo de CPU foram 3.11 s e 3.56 s, respectivamente. Na Figura 2.9 (a) é apresentado o gráfico da convergência da flexibilidade calculada pelo MBR em função do tamanho da base N utilizada, para as malhas citadas anteriormente. Na Figura 2.9 (b) é apresentado o tempo de CPU off-line, para a construção das bases.

Figura 2.9 Problema acoplado – Convergência do MBR e tempo off-line. O tempo “on-line” não variou significativamente com o número de pontos da amostra, sendo o tempo on-line médio para a malha de elementos T3 = 0.8s e para a malha de elementos Q4 = 1.3s. O MBR para o problema acoplado se mostra menos eficiente que nos outros tipos de solicitações (puramente elástico ou puramente térmico) devido à atualização do vetor de carga gerado pela deformação térmica que é desenvolvido durante o estágio “on-line”.

31

2.6 REFERÊNCIAS

AFONSO, S.M.B, LYRA, P.R.M, ALBUQUERQUE, T.M. M., R. S., MOTTA, 2009, “Structural Analysis and Optimization in the Framework of Reduced-Basis Method. Structural and Multidisciplinary Optimization”, Springer Berlin / Heidelberg, DOI: 10.1007/s00158-008-0350-4, Published online, January 2009. AFONSO, S. M. B. e PATERA, A. T. “Structural Optimization in Framework of Reduced Basis Method”. in XXIV CILAMCE, Congresso Ibero Latino Americano de Métodos Computacionais na Engenharia, 2003. ALBUQUERQUE T. M. M. “Análise e Otimização de Problemas Térmicos e Estruturais Bidimensionais Através do Método das Bases Reduzidas”. Msc. Thesis ( in Portuguese), Dept. de Engenharia Civil, Universidade Federal de Pernambuco, Recife-PE, Brasil, 2005. ALMEIDA, F. P. C. “OTIMIZAÇÃO DE FORMA EM PROBLEMAS TERMOESTRUTURAIS ACOPLADOS”. Dissertação de mestrado, Dept. de Engenharia Mecânica, Universidade Federal de Pernambuco, Recife-PE, Brasil, 2001. COOK, R. D., MALKUS, D. S., PLESHA, M. E. e WITT R. J. “Concepts and Applications of Finite Element Analysis”. Fourth edition, Wiley, 2003. MOTTA, R. S.; AFONSO, S. M. B.; LYRA, P. R. M. “Multiobjective Solutions for 2d Continua Problems Considering Reduced-Basis Approximations”. CMNE/CILAMCE 2007. Porto, Portugal, 2007. PRUD’HOMME, C., ROVAS, D. V., VEROY, K., MACHIELS, L., MADAY, Y., PATERA, A.T. e TURICINI, G. “Reliable Real-time Solution of Parameterized Partial Differential Equations: Reduced-basis Output Bound Method”. J. of Fluids Eng., Vol.124, pp.70-79, 2002. VILLAÇA, S. F. e GARCIA, L. F. T. “Introdução à Teoria da Elasticidade”. 4ª edição, COPPE/UFRJ, Rio de. Janeiro, RJ, Brasil, 2000. VEROY, K. “Reduced-Basis Methods Applied to Problems in Elasticity: Analysis and Applications” PhD Thesis. Department of Civil and Environmental Engineering, Massachusetts Institute of Technology, 2003. ZIENKIEWICZ O C & TAYLOR R L. “The finite element method. Vol. I. Basic formulations and linear problems”. London: McGraw-Hill. 648 p. Vol. 2, 1989.

32

33

CAPÍTULO 3 OTIMIZAÇÃO ESCALAR E ANÁLISE DE SENSIBILIDADE 3 OTIMIZAÇÃO ESCALAR E ANÁLISE DE SENSIBILIDADE 3.1 INTRODUÇÃO Os procedimentos de otimização de forma consistem numa abordagem geral para projetar estruturas dos mais variados campos da engenharia tais como Civil, Mecânicas, Aeroespaciais, Marítimas, etc. Através da variação de parâmetros de forma e/ou dimensões que definam a estrutura, um dado projeto inicial é melhorado com relação a um conjunto de funções objetivo e restrições pré-definidas. A pesquisa científica no campo da otimização (uni ou multiobjetivo) foi substancialmente aumentada durante as últimas décadas, e um progresso considerável foi atingido. Este desenvolvimento foi devido principalmente ao progresso de ferramentas confiáveis para a análise em geral, tais como o método dos elementos finitos, métodos de análise de sensibilidades e métodos de programação matemática. Tal procedimento foi fortemente alavancado pelo crescimento acentuado das velocidades e capacidades dos computadores digitais. Uma idéia geral das atividades em tal campo nas últimas décadas pode ser encontrada na literatura (BATES, 2003; VENKATARAMAN e HAFTKA, 2004; VANDERPLAATS, 2006). Com o objetivo de dimensionar formas estruturais eficientes, surge o processo de Otimização Estrutural de Forma (“Structural Shape Optimization” - SSO) que inclui, além do algoritmo de otimização adequado, a análise estrutural e a análise de sensibilidade, quando o método de otimização requer a determinação dos gradientes das funções objetivo e restrição. Desta forma, as alterações no processo de dimensionamento são feitas de forma automática, baseadas em informações sobre o comportamento do projeto, obtidas ao longo do processo de otimização. Para isso, essas técnicas de programação matemática podem ser integradas ao procedimento de análise por Elementos Finitos, resultando em uma excelente ferramenta automatizada para a otimização do projeto estrutural. Estes procedimentos podem ser associados também, a métodos aproximados como o Método da Base Reduzida, permitindo assim, uma grande economia computacional como será visto posteriormente. Para a obtenção de projetos ótimos através da otimização de forma será utilizado um algoritmo SSO (ALBUQUERQUE, 2005), o qual permite tanto a otimização seguindo o procedimento convencional (MEF) quanto através das aproximações via o Método da Base Reduzida. Este procedimento iterativo envolve vários processos de análise antes de atingir a solução ótima.

34

3.2 FORMULAÇÃO MATEMÁTICA Quando os problemas de otimização envolvem uma única função objetivo ele é denominado problema de otimização escalar (HAFTKA e GURDAL, 1993). Sua representação genérica é dada por: (3.1)

(3.1) min f ( x ) x∈

nvp

Sujeito às seguintes condições: hk ( x ) = 0

(3.2) g ( x ) ≤ 0 i

xlj ≤ x j ≤ xuj

k = 1,… , ne i = 1,… , ni

(3.2)

j = 1,… , npv

sendo: f(x) → função objetivo; h(x) → função restrição de igualdade; g(x) → função restrição de desigualdade; x → vetor das variáveis de projeto; ne → número de funções restrições de igualdade; ni → número de funções restrições de desigualdade; nvp → número de variáveis de projeto; nvp

→ espaço das variáveis de projeto (nvp dimensional).

As funções restrições de igualdade, h(x), e de desigualdade, g(x) definem o espaço viável ou admissível para a(s) variável(eis) de projeto. As restrições, em certo ponto x* viável, são ditas inativas, caso g i ( x ) < 0 e ativas caso g i ( x ) = 0 , todas as restrições de igualdade são consideradas ativas.

3.3 PROGRAMAÇÃO MATEMÁTICA A Programação Matemática pode ser considerada como a primeira linha de métodos para resolução de problemas de otimização. Ela trata o problema de forma iterativa e determinística, isto é, através de gradientes, funcionais, operações matriciais, etc., para encontrar o ponto ótimo. Na maioria dos casos de otimização, supõe-se a continuidade das funções, assim como de suas derivadas, o que nem sempre ocorre acarretando na principal dificuldade encontrada por estes tipos de métodos. O processo de otimização é abordado no espaço ℜnvp. A diferenciabilidade e a convexidade do problema influem, fundamentalmente, sobre a natureza das condições de ótimo.

35

Os métodos mais empregados para solucionar problemas de otimização nãolineares com restrições, são os métodos diretos (VANDERPLAATS, 2004). Estes métodos definem um subproblema para a determinação de uma direção de busca S e em seguida realizam uma busca linear para se determinar o tamanho do passo α, nessa direção. A forma clássica da atualização das variáveis, por este procedimento iterativo, é dada por: (3.3)xk+1 = xk + αSk+1

(3.3)

onde S é o vetor da direção de busca e α representa o tamanho do passo na direção de S. A seleção do algoritmo de otimização depende fundamentalmente do problema envolvido. Isto é importante para se obter uma otimização confiável e um alto nível de eficiência computacional (tempo de computação, taxa de convergência, memória requerida). Vários são os algoritmos utilizados na resolução dos problemas de otimização nãolineares com restrições baseadas no método direto (KLEIBER, 2005). A resolução dos problemas de otimização não-linear com restrições pode ser feita por meio dos chamados métodos indiretos. Estes transformam o problema original definido pela Equações (3.1) e (3.2), num problema sem restrições através de algum artifício como: penalidade, dualidade, lagrangeano aumentado, etc... (BATES, 2003; HAFTKA e GURDAL, 1993). O algoritmo usado neste trabalho, na busca da solução ótima para os problemas de otimização, é o da programação quadrática seqüencial, mais conhecido como SQP (“Sequential Quadratic Programming”). As etapas e as equações básicas deste método estão descritas na seção 2.4.2.

3.3.1 Condições de ótimo para problemas com restrições Um ponto x é dito regular, caso os gradientes das funções restrição ativas sejam linearmente independentes. Neste caso o número de restrições ativas deve ser menor ou igual ao número de variáveis de projeto. Um ponto regular viável x, para ser dito um ponto de mínimo local (x*), deve atender às condições necessárias de primeira ordem, dadas pelas condições de KarushKuhn-Tucker (condições KKT) (KARUSH, 1939; KUHN e TUCKER, 1950). Neste ponto o gradiente da função objetivo deve ser capaz de se anular (mesma direção e sentido opostos), como uma combinação linear dos gradientes das restrições ativas, e os coeficiente lineares relacionados às restrições de desigualdade ativas devem ser não negativos. As condições de KKT podem ser escritas como: ne

ni

k =1

i =1

(3.4) ∇f (x* ) + ∑ λk*∇hk (x* ) + ∑ φi*∇gi (x* ) = 0

(3.4)

(3.5) φi* ≥ 0, i = 1..ni

(3.5)

(3.6) φi* gi (x* ) = 0, i = 1..ni

(3.6) 36

onde λk e φi, são os multiplicadores de Lagrange associados respectivamente às restrições hk, gi, no ponto x e formam as variáveis duais do problema. Os coeficientes φi, também são conhecidos como os multiplicadores de KKT. Os asteriscos indicam se tratar de um ponto de ótimo, incluindo os multiplicadores de Lagrange ótimos. Nas restrições gi estão incluídas também, as 2npv restrições de lateralidade definidas na Equação (3.2). A Equação (3.4) é chamada de condição de estacionaridade. A Equação (3.5) é conhecida como as condições de viabilidade dual e impõe que os multiplicadores de Lagrange associados às restrições de desigualdade (multiplicadores de KKT) ativas em x* sejam positivos (ou nulas para restrições redundantes). Isto porque as restrições são escritas na forma gi ≤ 0. A Equação (3.6) é chamada de condição de complementaridade, ela implica que uma restrição inativa tem um multiplicador nulo e que uma restrição ativa pode ter multiplicador diferente de zero. É conveniente utilizar a função Lagrangeana, para o desenvolvimento metodológico do processo de otimização. A função Lagrangeana é definida por (CHONG, 2001): ne

ni

k =1

i =1

(3.7) L ( x, λ , φ ) = f (x) + ∑ λk hk (x) + ∑ φi gi (x)

(3.7)

A Equação (3.4) representa, então, o gradiente da função Lagrangeana em relação a x num ponto de ótimo, ou seja, ne

ni

k =1

i =1

(3.8) ∇L ( x, λ , φ ) = ∇f (x) + ∑ λk ∇hk (x) + ∑ φi ∇gi (x) = 0

(3.8)

Para as condições de ótimo de segunda ordem, é necessário definir S que representa uma direção de decréscimo da função objetivo pertencente ao subespaço tangente das direções viáveis T, onde pequenas perturbações nestas direções mantêm a viabilidade do problema. O subespaço tangente T é definido como: (3.9)T = N (∇h(x)) ∪ N (∇g j (x)) , j ∈ J (x)

(3.9)

onde N (⋅) representa o núcleo e J (x) é o conjunto dos índices das restrições de desigualdade ativas, ou seja, (3.10) J (x) = { j : g j (x) = 0}

(3.11)

(3.10)

N (∇h(x)) = núcleo(∇h(x)) = {S | S ≠ 0 e ST ∇h(x) = 0} N (∇g j (x)) = núcleo(∇g j (x)) = {S | S ≠ 0 e ST ∇g j (x) ≤ 0} ,

j ∈ J ( x)

(3.11)

37

Assim, as direções viáveis estacionárias S pertencem a N (∇h(x)) e N (∇g j (x)) , que são os espaços tangentes das superfícies definidas pelas restrições de igualdade e pelas restrições de desigualdade ativas ( g j (x) = 0 ), respectivamente. A condição necessária de segunda ordem pode ser enunciada da seguinte maneira (CHONG, 2001): Se existirem os vetores λ, φ tais que as restrições (3.4) a (3.6) sejam satisfeitas, então a condição dada na equação abaixo é necessária para x* ser um mínimo local. (3.12) ST HS ≥ 0,

∀S ∈ T

(3.12)

ou ainda, podemos ter a condição suficiente para x* ser um mínimo local estrito. (3.13) ST HS > 0,

∀S ∈ T

(3.13)

sendo S qualquer direção de decréscimo da função objetivo pertencente ao subespaço tangente das direções viáveis T e H representa a matriz Hessiana da função Lagrangeana em relação a x: ne

ni

k =1

i =1

(3.14) H = ∇ 2 L ( x, λ , φ ) = ∇ 2 f (x) + ∑ λk ∇ 2 hk (x) + ∑ φi ∇ 2 gi (x)

(3.14)

3.3.2 Programação Quadrática Sequencial - SQP No método da programação quadrática sequencial, o problema é resolvido por meio de uma seqüência de subproblemas quadráticos convexos aproximados (POWEL, 1978). Partindo-se da solução do passo anterior, para se obter um novo incremento das variáveis de projeto, resolve-se uma seqüência de problemas quadráticos aproximados. Essa seqüência deve convergir para um ponto de mínimo x*. Assim, dado x, procura-se o valor de S que resolve o problema original: (3.15) min f ( x + S ) S∈

(3.15)

nvp

Sujeito as seguintes condições: (3.16) hk ( x + S ) = 0 gi (x + S) ≤ 0

k = 1,… , ne i = 1,… , n i

(3.16)

A função Lagrangeana associada às Equações (3.15) e (3.16) é:

38

(3.17) L ( x + S, λ , φ ) = f (x + S) + λ T h(x + S) + φT g (x + S)

(3.17)

Sendo h = (h1,..., hne)T e g = (g1,..., gni)T vetores que contêm as restrições e λ = (λ1,..., λne)T e φ = (φ1,..., φni)T os vetores que contêm os multiplicadores de Lagrange de dimensão ne e ni, respectivamente. As condições de ótimo de primeira ordem aplicadas ao problema (3.15) e (3.16) fornecem as equações:

(3.18) ∇L ( x + S, λ , φ ) = ∇f (x + S) + λ T ∇h(x + S) + φT ∇g (x + S) = 0

(3.18)

(3.19) h(x + S) = 0

(3.19)

(3.20) g(x + S) ≤ 0

(3.20)

(3.21) φT g(x + S) = 0

(3.21)

(3.22) φi ≥ 0, i = 1..ni

(3.22) (3.23)

onde ∇g e ∇h são matrizes Jacobianas cujas linhas são os gradientes de ∇gi e ∇hk, respectivamente. Aproximando as Equações de (3.18) até (3.21) por uma série de Taylor até termos de primeira ordem em x resulta: (3.24) ∇L(x) + ∇ 2 L(x)S = 0

(3.23)

(3.25) h(x) + ∇h(x)S = 0

(3.24)

(3.26) g (x) + ∇g (x)S ≤ 0

(3.25)

(3.27) φT [ g (x) + ∇g (x)S ] = 0

(3.26)

onde ∇L(x) é o gradiente e ∇ 2 L(x) é a matriz Hessiana da função Lagrangeana, definidos por: (3.28) ∇L(x) = ∇f (x) + λ T ∇h(x) + φT ∇g(x)

(3.27)

(3.29) ∇ 2 L(x) = ∇ 2 f (x) + ⎡⎣ λ∇ 2h(x) ⎤⎦ + ⎡⎣φ∇ 2 g (x) ⎤⎦

(3.28)

ne

ni

i =1

i =1

Utilizando a seguinte notação: ⎡⎣ λ∇ 2h(x) ⎤⎦ = ∑ λi ∇ 2 hi (x) e ⎡⎣ φ∇ 2 g(x) ⎤⎦ = ∑ φi ∇ 2 gi (x) . As Equações (3.22) e de (3.23) a (3.26) são as condições de otimalidade do subproblema de determinação da direção de busca S dado por:

39

1 (3.30) min ∇L(x)T S + ST ∇ 2L(x)S ndv S∈ 2

(3.29)

Sujeito as seguintes restrições (3.31)

h(x) + ∇h(x)S = 0 g(x) + ∇g(x)S ≤ 0

(3.30)

Este é o esquema básico de linearização que, demonstra-se, é localmente de segunda ordem, desde que satisfeitas algumas hipóteses, evidentemente, uma delas a condição de ótimo de segunda ordem. As dificuldades na solução deste subproblema são a determinação dos multiplicadores de Lagrange e o cálculo a cada iteração da Hessiana ∇2L que depende dos multiplicadores, sendo geralmente, um processo caro. Assim, surgiram os métodos Quase-Newton, nos quais uma aproximação da matriz Hessiana da função Lagrangeana ∇2L é construída a partir dos gradientes ∇L da Lagrangeana ao longo das iterações. Estes métodos possuem convergência superlinear e são largamente utilizados na solução dos problemas de otimização. Entre os métodos Quase-Newton, o Broyden-Fletcher-Goldfarb-Shann (BFGS) é o utilizado neste trabalho. Entretanto, a questão dos multiplicadores permanece, e é necessário utilizar alguma técnica para estimá-los (VANDERPLAATS, 2004). Após a obtenção da direção de busca S, é necessário calcular o tamanho do passo nessa direção, a fim de se obter o novo vetor das variáveis de projeto. A forma como estes dois componentes são determinados reflete, fundamentalmente, na eficiência e confiabilidade do algoritmo de programação para cada caso. O parâmetro tamanho do passo é determinado de modo a produzir o maior decréscimo na função mérito. A função mérito usada aqui, é a própria função Lagrangeana, porém os multiplicadores de Lagrange são parâmetros de penalidade (POWELL, 1978). Todo o processo de otimização (obter direção de busca, busca linear, atualização da matriz hessiana, atualização das variáveis, etc..) é feito automaticamente, através de ferramentas do ambiente computacional MATLAB (MATHWORKS, 2007), utilizadas neste trabalho.

3.4 ANÁLISE DE SENSIBILIDADES A análise de sensibilidades constitui uma fase importante de qualquer algoritmo de otimização que necessita do cálculo dos gradientes. A análise de sensibilidade avalia a resposta da função objetivo ou das restrições às alterações das variáveis de projeto e normalmente está relacionada com procedimentos numéricos para o cálculo das derivadas destas funções e restrições em relação às variáveis de projeto. Vários são os procedimentos para o cálculo de sensibilidades (CHOI et al, 2005; KEULEN et al, 2005). Na elaboração deste trabalho foram utilizados dois métodos, são eles o método das diferenças finitas (MDF) e método direto. A seguir, os métodos usados serão discutidos e em seguida aplicados na solução de exemplos-padrão. O método das diferenças finitas apresenta a vantagem de ser de fácil implementação, porém podem induzir a erros elevados. Já os métodos direto e adjunto, que não fazem uso de cálculos via diferenças finitas, apresentam como mérito principal a preci40

são na solução e a eficiência computacional, mas sua implementação é mais elaborada e consequentemente mais propícia a erros.

3.4.1 Método das diferenças finitas globais É um dos mais antigos métodos usados para a análise de sensibilidade. Possui como característica principal a simplicidade de implementação do algoritmo. Apresenta um inconveniente quando o número de variáveis de projeto é grande, pois o custo computacional passa a ser alto. Isto ocorre porque para cada variável de projeto perturbada é necessária a resolução das equações governantes. Um outro problema é encontrar um valor do tamanho da perturbação que modifique as variáveis de projeto de tal forma que forneça uma boa precisão para os gradientes (HAFTKA e GURDAL, 1993; CHOI et al, 2005). Sua formulação vem da representação da função em análise na forma da série de Taylor, cuja forma é:

(3.32) f ( μ + ei Δμi ) = f ( μ ) +

∂ 2 f ( μ + eiδΔμi ) Δμi2 ∂f ( μ ) Δμi + ∂μi ∂μi 2 2!

(3.31)

para algum δ entre 0 e 1. Truncando-se os termos de ordem superior ( Δμi2 ) e realizando-se algumas operações matemáticas, tem-se:

(3.33)

∂f ∂μi

f ( μi + ei Δμi ) − f ( μi ) Δf = , para i = 1,..., m Δμi Δμi

(3.32)

O principal fator que influencia a precisão deste método é o tamanho do passo Δμi usado. Quando as sensibilidades são obtidas usando o método das diferenças finitas, o dilema em relação ao tamanho da perturbação surge (HAFTA e GURDAL, 1993; KEULEN et al, 2005), pois grandes perturbações geram erros de truncamento (estão relacionados aos termos negligenciados no truncamento da série de Taylor) e pequenas amplificam os erros de arredondamento. No caso deste método o tamanho do passo em torno de 10-4 μ i a 10-6 μ i é tido como satisfatório (AFONSO, 1995).

3.4.2 Método direto Este método possibilita a obtenção de gradientes com custo computacional mais baixo quando comparado ao MDF, anteriormente apresentado. O gradiente é determinado aplicando-se a derivada de primeira ordem à função em análise com relação a uma dada variável de projeto. Apesar deste método não ser tão simples quanto o anterior é computacionalmente mais eficiente. Para demonstrar este método parte-se da forma discretizada da equação governante (CHOI et al, 2005). No caso de problemas elásticos, 41

(3.34) Ku = F

(3.33)

onde K é a matriz de rigidez, u é o vetor de deslocamento e F é o vetor de carregamento. Derivando a equação em relação à variável de projeto x tem-se

(3.35) K

∂u ∂K ∂F + u= ∂μk ∂μk ∂μk

(3.34)

o gradiente dos deslocamentos pode ser calculado, então, da seguinte forma:

(3.36)

⎛ ∂F ∂K ⎞ ∂u u⎟ = K −1 ⎜ − ∂μk ⎝ ∂μ k ∂μ k ⎠

(3.35)

∂K ∂F e sejam calculadas indiretamente, como por diferenças ∂x ∂x finitas, tem-se as versões do método direto, denominada semi-analítica (convencional ou exato) (AFONSO, 1995). Porém, através da decomposição por região e do mapeamento aplicado ao domínio de referência, abordado no seção 2.3.2, pode-se obter as derivadas analíticas das matriz de rigidez e vetor carregamento de cada região. Para isto, basta derivar as Equações (2.48) e (2.49) em relação a k-esima variável de projeto, logo Caso as derivadas

r ∂K *s r ( μ ) nt ∂β s ( μ ) j r (3.37) =∑ Kj ; ∂μk ∂μk j =1

∂Fs*r ( μ ) ∂ϕΩr (μ) r ∂ϕΓr (μ) r Ff = Fb + ∂μk ∂μk ∂μk

(3.36)

Os gradientes da matriz de rigidez e vetor carregamento são calculados somando-se a contribuição de cada região ∂K *s nr ∂K s ( μ ) ∂Fs* nr ∂F ( μ ) ; =∑ =∑ ∂μk r =1 ∂μk ∂μk r =1 ∂μk *r

(3.38)

*r

(3.37)

Através deste método podemos usar a matriz de rigidez já decomposta e para cada novo x calculamos apenas as derivadas dos parâmetros variáveis (dependentes de x), estas derivadas são analiticamente conhecidas, obtendo assim um ganho computacional. Apesar das equações para a análise de sensibilidade apresentadas, estarem relacionadas com as equações do problema de elasticidade, estas são estendidas para problemas de condução de calor tratados aqui (regime estacionário), bem como para o caso acoplado.

3.5 INTEGRAÇÃO ANÁLISE/OTIMIZAÇÃO

O sistema computacional desenvolvido neste trabalho incorpora todos os elementos necessários para a otimização de forma/dimensões, incluindo: modelagem geométrica,

42

geração de malhas, análise estrutural e/ou térmica, análise de sensibilidades e o módulo que integra o procedimento completo com um algoritmo de programação matemática. Este sistema permite executar uma otimização de projetos estruturais envolvendo efeitos mecânicos, termo-mecânicos acoplados em regime estático, e otimização térmica em regime permanente. Na Tabela 3.1 é mostrado o Algoritmo SSO, que resume o procedimento utilizado na otimização de forma. Tabela 3.1 Algoritmo SSO: procedimento de otimização de forma. SSO1. Definir o problema de otimização de forma; SSO2. Gerar a geometria e discretizar o problema; SSO3. Efetuar a análise estrutural, térmica ou termo-estrutural conforme o caso; SSO4. Efetuar a análise de sensibilidades (via Diferenças Finitas ou Método Direto

Analítico) ; SSO5. Obter um novo projeto utilizando o procedimento de programação matemática; SSO6. Verificar o critério de convergência para a otimização de forma: SSO6.1 parar se o novo projeto satisfaz ao critério; SSO6.2 caso contrário atualizar o projeto e retornar ao passo SSO2.

3.5.1 O método da base reduzida no procedimento de otimização

Nesta seção serão descritas algumas peculiaridades relacionadas ao MBR no procedimento de otimização. No algoritmo SSO apresentado na Tabela 3.1, esta técnica se insere no módulo de análise (avaliação de funções) e no módulo de análise de sensibilidades (avaliação dos gradientes das funções). O procedimento aplicado para análise via MBR foram discutidos nos capítulo 2. No que se segue, apresenta-se a análise de sensibilidades via MBR, onde as equações foram particularizadas para problemas de elasticidade.

(a) Análise de sensibilidade através do MBR

O uso das matrizes decompostas para as matrizes de rigidez e vetores de carregamento de cada região, além do mapeamento e projeção das grandezas, permite que o Método da Base Reduzida aplique o Método Direto Analítico na análise de sensibilidade, de forma simplificada. Estas matrizes e vetores são independentes do valor da variável de projeto vide Equações (3.36) e (3.37). Com isso, para cada novo parâmetro μ executamos a análise de sensibilidade com um menor custo computacional que os métodos baseados nas equações tradicionais provenientes do MEF.

43

No caso do cálculo das sensibilidades utilizando o método das diferenças finitas, o MBR também é vantajoso pois as análises são obtidas de maneira mais rápida que o método convencional (totalmente via MEF). No caso de problemas de elasticidade por exemplo, a análise de sensibilidade via o método direto é conduzida conforme descrito a seguir (CHOI et al, 2005). Utilizando a equação dos deslocamentos aproximados Equação (2.62), e aplicando diferenciação direta tem-se

(3.39)

∂u N ∂α N =Z ∂μk ∂μk

(3.38)

onde a derivada dos coeficientes lineares α , pode ser obtida derivando-se diretamente a Equação (2.66),

(3.40)

∂K N ∂α ∂F N α + KN = ∂μk ∂μk ∂μk

⇒

N −1 ⎛ ∂F ∂α ∂K N ⎞ α⎟ = ⎡⎣K N ⎤⎦ ⎜ − ∂μk ⎠ ∂μk ⎝ ∂μk

(3.39)

ou na forma compacta (KEULEN et al, 2005)

(3.41)

−1 ∂α = ⎡⎣K N ⎤⎦ Fˆ N ∂μk

(3.40)

onde Fˆ N é o pseudo vetor de força da base reduzida, definido como ⎛ ∂F N ∂K N ⎞ N ˆ (3.42) F = ⎜ − α⎟ ∂μk ⎠ ⎝ ∂μk

(3.41)

De forma análoga à diferenciação no MEF, através da decomposição por região, do mapeamento aplicado ao domínio de referencia abordado no seção 2.3.2, e aplicando as projeções na base reduzida (seção 2.4), pode-se então derivar as Equações (2.70) em relação a k-esima variável de projeto, para se obter, no domínio real, as derivadas das matriz de rigidez e vetor carregamento de cada região no espaço de base reduzida WN. Somando-se as contribuição de cada região, os gradientes no domínio real, das matriz de rigidez e vetor carregamento totais, na base reduzida, são obtidos através de r ∂K *s N ( μ ) nr nt ∂βs ( μ ) j Nr = ∑∑ Ks j ∂μk ∂μk r =1 j =1

(3.43)

(3.42) ∂F ( μ ) nr ⎡ ∂ϕ (μ) Nr ∂ϕ (μ) Nr ⎤ = ∑⎢ Fb + Ff ⎥ ∂μk ∂μk r =1 ⎣ ∂μ k ⎦ *N s

r Ω

r Γ

44

As Equações (3.38), (3.40), (3.41) e (3.42) são desenvolvidas para cada nova variável de projeto ( μ ). O algoritmo completo para a implementação computacional via o Método da Base Reduzida Tabela 2.1, ficará acrescido de um item na etapa “on-line” para a análise de sensibilidade, após o cálculo do vetor uN , o cálculo do gradiente de uN através da Equação (3.38). Para comprovar a eficácia do método, são apresentados na seção seguinte alguns exemplos envolvendo análise de sensibilidade e otimização.

3.6 EXEMPLOS

Neste exemplo iremos realizar via Método dos Elementos Finitos e Método da Base Reduzida a análise de sensibilidade e otimização do problema apresentado na seção 2.5.1 e que se refere à análise do estado plano de tensões em uma placa quadrada com um orifício central. Neste exemplo prático as dimensões do orifício são as variáveis do projeto selecionadas para a otimização. Os valores iniciais das variáveis de projeto são μ1 = 60 mm e μ2 = 43 mm, e os limites inferior e superior impostos são 25mm e 75mm, respectivamente. O objetivo considerado é minimizar a flexibilidade. Ao volume, é imposto ser inferior ou igual ao seu valor inicial. O problema de otimização pode ser formulado como: Minimizar: f ( μ ) sujeito a V ( μ ) ≤ V0 25mm ≤ μk ≤ 75mm k = 1,...ndv

Onde f ( μ ) é a flexibilidade, V ( μ ) o volume da estrutura, V0 o volume inicial e

ndv o número de variáveis de projeto. Serão considerados vários modelos de elementos finitos, onde o número de graus de liberdade será variado através da alteração na dimensão média dos elementos, na configuração de referência onde μ1 = μ2 = 50mm. Uma análise de sensibilidade foi realizada, considerando as malhas utilizadas para o estudo do comportamento do MBR para o cálculo da flexibilidade, definidas na seção 2.5.1. O gradiente da flexibilidade foi calculado via MEF e via MBR, o número de amostras utilizada neste caso foi 8 e os resultados obtidos estão ilustrados na Tabela 3.2. Pode-se notar que para poucas reavaliações da sensibilidade, o tempo total via MEF superara o tempo total via MBR, pois as reavaliações via MBR são feitas todas no estagio on-line (o tempo da avaliação é o tempo “on-line”).

45

Tabela 3.2 Resultados da análise de sensibilidade. ∂f ( μ ) ∂μ1 MEF – T3 -0.0085236 MBR – T3 -0.0085230 MEF – Q4 -0.0085210 MBR – Q4 -0.0085206

∂f ( μ ) Tempo (s) da Tempo (s) avaliação “off-line” ∂μ2 -0.015176 17.05 -0.015177 0.064 75.5 -0.015172 8.75 -0.015174 0.054 45.1

Para o processo de otimização via MEF, o valor da dimensão média dos elementos foi variada de 10 mm à 1 mm, resultando em modelos de EF (Elementos Finitos) que variam de malhas grosseiras até as mais refinadas. Para cada refinamento de malha será feita uma otimização. Analisou-se a função objetivo e as variáveis de projeto ótimas resultantes das otimizações. A Figura 3.1 ilustra o resultado de cada otimização variando com o número de graus de liberdade da malha utilizada para a otimização. Na Figura 3.1 (a) é apresentado a convergência do valor ótimo de X1 (X1 ótimo), onde X1= 50 − μ1 se refere à primeira variável de projeto, com a variação do número de graus de liberdade (NGL). Já na Figura 3.1 (b), vemos o gráfico da função objetivo (flexibilidade) no ponto ótimo versus o NGL.

a) variável X1

b) função objetivo f

Figura 3.1 Um quarto da placa quadrada – Variação do ótimo com a variação do número de graus de liberdade: a) valor da primeira variável no ponto ótimo, b) valor da função objetivo. Esses gráficos ilustram a influência da aproximação do MEF no resultado final da otimização, ressaltando a importância de uma boa modelagem. Também se pode imaginar um processo mais eficiente, utilizando recursivamente como ponto inicial, um ponto ótimo obtido com uma malha mais grosseira, para obter um ponto de ótimo de um modelo mais refinado, porém este procedimento não foi testado aqui. Foi examinado o procedimento de otimização, realizando a etapa das análises via MBR. Para tal, as malhas de elementos finitos consideradas foram a que apresentaram resultados satisfatórios para a função objetivo no ponto ótimo, são elas: 46

•

Para o elemento Q4 adotou-se elementos de dimensão média de 10/3 mm. O modelo possui 1.472 graus de liberdade.

•

Para o elemento T3 adotou-se elementos de dimensão média de 2 mm. O modelo possui 3.952 graus de liberdade.

A base foi construída no espaço das variáveis de projeto, D ={[25, 75]²}. Foram realizados vários processos de otimização para o problema considerado, variando o número de pontos amostrados para a base reduzida N de 5 à 16. Para verificar a convergência do MBR, construímos o gráfico apresentado na Figura 3.2 (a) do valor referente à primeira variável no ponto ótimo (X1 ótimo) variando com número de pontos amostrados N, e o gráfico do objetivo (flexibilidade) ótimo variando com N, Figura 3.2 (b). No gráfico podemos verificar a convergência para poucas soluções base. A Figura 3.3 ilustra os tempos de CPU para o procedimento “off-line” e para o processo de otimização (“on-line”) via MBR.

a) variável X1

b) função objetivo f

Figura 3.2 Um quarto da placa quadrada – Convergência do ponto ótimo com o tamanho da base.

a) Tempo “off-line”

b) Tempo “on-line”

Figura 3.3 Um quarto da placa quadrada: a) Tempo “off-line”, b) tempo de otimização “on-line”. 47

A Tabela 3.3 resume o desempenho de otimizações via MEF e via MBR, onde X1 e X2 são respectivamente 50 − μ1 e 50 − μ2 . As malhas utilizada foram as mesmas consideradas para o estudo da otimização via MBR e o número de pontos amostrados para a aproximação do MBR foi 11. Tabela 3.3 Comparativo dos resultados das otimizações: MEF x MBR X1 MEF – T3

X2

Flexibilidade

Tempo

Tempo

Tempo

“on-line”(s) “off-line” (s)

Total

-21.5248 13.9286

0.37320

-

-

67.3

MBR – T3 -21.5410 13.9368 MEF – Q4 -21.5051 13.9187 MBR – Q4 -21.5149 13.9236

0.37317

0.6

14.7

15.3

0.37360

-

-

118.6

0.37356

0.3

25.7

26.0

Os resultados comprovam a eficiência do método, principalmente para problemas de otimização complexos onde são exigidos muitos cálculos das funções de interesse. Além disso, o estágio “off-line” pode ser facilmente paralelizável diminuindo significativamente o tempo “off-line”, consequentemente o tempo total da otimização.

3.7 REFERÊNCIAS

AFONSO, S. M. B. “Shape Optimization of Shells Under Static and Free Vibrations Conditions”. Tese de Doutorado, University of Swansea, Wales, U.K. 1995. AFONSO, S. M. B.; MACEDO, C. M. H.; OLIVEIRA, D. A. P. “Structural Shape Optimization under Multicriteria Conditions In: V World Congress on Computational Mechanics, Viena. 2002. ALBUQUERQUE T. M. M. “Análise e Otimização de Problemas Térmicos e Estruturais Bidimensionais Através do Método das Bases Reduzidas”. Msc. Thesis ( in Portuguese), Dept. de Engenharia Civil, Universidade Federal de Pernambuco, Recife-PE, Brasil, 2005. BATES, S. “Development of Robust Simulation, Design and Optimization Techniques for Engineering Applications”. PhD Thesis. School of Engineering, University of Wales, Swansea. 2003. CHOI K.K.; KIM N. H. “Structural Sensitivity Analysis and Optimization 1 – Linear Systems”. Springer Mechanical Engineering Series, Springer New York, 2005.

48

CHONG, E. K. P.; ZAK, S. H. “An Introduction to Optimization”, Second edition, John Wiley & Sons, 2001. HAFTKA, R.T.; GURDAL, Z. “Elements of Structural Optimization”. Third Revised and Expanded Edition, Kluwer Academic Publishers, 1993. KARUSH, W. “Minima of Functions of Several Variables with Inequalities as Side Constraints”. M.Sc. Dissertation. Dept. of Mathematics, Univ. of Chicago, Chicago, Illinois, 1939. KEULEN F. V.; HAFTKA, R. T.; KIM N. H. “Review of options for structural design sensitivity analysis. Part 1: Linear systems”. Comput. Methods Appl. Mech. Engrg. 194, pp. 3213–3243, 2005. KIRSCH, U. “Structural Optimization Fundamentals and Applications”. Springer Verlag, 1993. KLEIBER M. “Parameter Sensitivity in Nonlinear Mechanics – Theory and Finite Element Computations”. John Wiley & Sons, 1997. KUHN, H. W.; TUCKER, A. W. “Nonlinear Programming”. In: Procedeering 2nd Berkeley Symposium Mathematics Statistics and Probability, 481, 1950. MACEDO, C. M. H. “Otimização de Treliças Planas sob Várias Solicitações com Ênfase a Problemas Multiobjetivos”. dissertação de Mestrado, Universidade Federal de Pernambuco. Recife-PE, Brasil, 2002. MATHWORKS. “MATLAB User’s Guide”. Mathworks Inc., Natacki, 2007. POWEL, M. J. D. “Algorithms for Nonlinear Constraints that use Lagrangian Function”. Math. Programming, vol 14, pp. 224-248, 1978. ROCKAFELLAR, R. T. “LAGRANGE MULTIPLIERS AND OPTIMALITY”. Vol. 35, No. 2, pp. 183-238, June 1993. VANDERPLAATS, G. N. “Numerical Optimization Techniques for Engineering Design - with Applications”. 4rt Edition, Vanderplaats Research & Development, Inc., Colorado Springs, CO, 2004.

49

VANDERPLAATS, G. N. “Structural Optimization for Statics, Dynamics, and Beyond”. Journal of the Brazilian Society of Mechanical Sciences and Engineering. vol. 28, no.3, p.316-322. ISSN 1678-5878, September 2006. VENKATARAMAN, S.; HAFTKA, R.T. “Structural optimization complexity: what has Moore’s law done for us?”. Structural and Multidisciplinary Optimization, 28(6), pp. 375-387, 2004.

50

CAPÍTULO 4 OTIMIZAÇÃO MULTIOBJETIVO 4 OTIMIZAÇÃO MULTIOBJETIVO O projeto ótimo de problemas reais quase sempre está envolvido com várias metas (funções objetivo) a serem aprimoradas e várias restrições a serem satisfeitas. Os algoritmos de otimização de propósito geral (existentes na literatura) são capazes de resolver problemas que envolvam uma única função objetivo. A maioria dos otimizadores disponíveis na literatura não lida, portanto com várias funções objetivo simultaneamente. No entanto, o problema de otimização multiobjetivo (POM) pode ser resolvido empregando-se uma das estratégias: •

O uso do conceito de Pareto

•

O uso da Programação Hierárquica

Na presente pesquisa, serão consideradas metodologias baseadas no uso do conceito de Pareto. Nos últimos 15 anos distribuições eficientes de pontos de Pareto têm sido obtidas graças ao desenvolvimento de algoritmos eficientes tais como o NBI (Normal-Boundary Intersection) (DAS e DENNIS, 1996) e o NNC (Normalized Normal Constraint) (MESSAC et al, 2003 ). Essas estratégias juntamente com outras abordagens de mais fácil implementação (Método da soma ponderada e Método min-max) são aqui implementadas.

4.1 DEFINIÇÃO DO PROBLEMA O problema multiobjetivo pode ser expresso como: (4.1) min F ( x ) = ⎡ f1 ( x ), f 2 ( x ), f 3 ( x ),… , f nobj ( x ) ⎤ , nobj ≥ 2 … ⎣ ⎦ x

(POM)

(4.1)

Sujeito às seguintes condições: hk ( x ) = 0

(4.2) g ( x ) ≤ 0 i

xlj ≤ x j ≤ xuj

k = 1,… , ne i = 1,… , ni

(4.2)

j = 1,… , npv

onde as restrições são as mesmas definidas no capítulo 3, porém os objetivos agora formam um vetor de nobj funções objetivo, as quais precisam ser minimizadas.

52

4.2 CONCEITO DE ÓTIMO DE PARETO Nos problemas de otimização multiobjetivo encontrar um x* que minimize várias funções objetivo simultaneamente é uma tarefa muito árdua, senão impossível para a quase totalidade dos problemas. Uma forma de determinar um x que satisfaça em parte os POM (Equações (4.1) e (4.2)) está contida na definição de Otimalidade de Pareto (ARORA et al., 2007). Pontos de Pareto são pontos xP tais que não exista nenhum ponto x o qual: i) f k ( x) ≤ f k ( x P ) para k = 1, … , n

ii) f k ( x) < f k ( x P )

para uma função objetivo ao menos

Os pontos de Pareto apresentam a propriedade de que quando se movem na direção decrescente de uma das funções, pelo menos uma das outras funções restantes tem seu valor aumentado. Vê-se isso na Figura 4.1, onde o ponto ótimo de Pareto é qualquer ponto no intervalo x1 ≤ x ≤ x2. Também, devido às restrições, o ponto ótimo de Pareto pode estar localizado ao longo do contorno da região viável. f

f2

f1

f1*

f2* x2

x1

x

Figura 4.1 Problema de otimização com uma variável e duas funções objetivo. Em problemas irrestritos, os ótimos de Pareto são pontos onde os gradientes das funções objetivo se anulem através de coeficientes lineares não-negativos. No caso biobjetivo, os gradientes possuem direções opostas (Figura 4.2). f2

x1

xP f1

x

x2 Figura 4.2 Pareto mínima com duas variáveis de projeto e duas funções objetivo.

53

Em problemas de otimização multiobjetivo é muito importante formular o problema no espaço das funções objetivo. Isto pode ser feito usando-se um sistema de equações geradas pelas funções objetivo e conjuntos das restrições ativas. Para cada projeto viável, haverá correspondentes valores das funções objetivo que definirão o espaço viável das funções objetivo. Sobre seu contorno se localizam os pontos ótimos de Pareto. Na Figura 4.3, tem-se o exemplo de um problema com duas variáveis de projeto e duas funções objetivo. Em ambas as Figura 4.3 (a) e (b), a linha tracejada representa os pontos ótimos de Pareto. x2

f2 x

F(x)

P

x

F(xP) f1

x1 a)

b)

Figura 4.3 Região viável e pontos de Pareto: a) no espaço das variáveis de projeto, b) no espaço das funções objetivo Conforme já mencionado, em problemas multiobjetivo, o interesse do projetista é encontrar um vetor das variáveis de projeto x* tal que as Equações (4.1) e (4.2) sejam satisfeitas. Usualmente, não existe tal x* devido ao aspecto de conflito comum entre as funções objetivo. Usando o conceito de Pareto, o projetista tem que encontrar tantos pontos quanto possíveis. A partir destes pontos, será escolhido o projeto o qual irá satisfazer, mais adequadamente, cada função objetivo.

4.3 MÉTODOS PARA GERAÇÃO DE PONTOS DE PARETO

Existem várias técnicas para se obter o chamado conjunto de mínimo de Pareto (DAS e DENNIS, 1997; BATES, 2003; COLLETTE, 2004). Entre eles, métodos baseados em metodologias metaheurísticas têm sido usados com algum sucesso (LALONDE et al 2009). Porém, neste trabalho apenas serão considerados procedimentos que fazem uso de programação matemática. Os métodos considerados serão os seguintes: Método da Soma Ponderada (STEUER, 1985), Método Min-Max (HWANG et all, 1980), Método NBI (DAS e DENNIS, 1997) e o Método da Restrição Normal Normalizada (MESSAC et al, 2003).

4.3.1 Método da Soma Ponderada dos Objetivos

Dentre os métodos desenvolvidos para otimização multiobjetivo, no qual se substitui as funções objetivo por uma única função, denominada de função substituta, o mais 54

empregado e de uso mais simples é o método da soma ponderada (“Weight Sum method” – WS) (KOSKI, 1985; AFONSO, 1997; AFONSO e SIENZ, 1999; AFONSO et al, 2002). Sua técnica baseia-se em atribuir um vetor de coeficientes de ponderação β j , às funções objetivo normalizadas, combinando-as linearmente, ou seja, transformandoas em uma única função objetivo. Sua representação algébrica é dada da seguinte forma

f f nobj = ∑ β j ,k k (4.3) F = β f0k f0 k =1

(4.3)

T j

onde os elementos de β j ,k são normalizados, da seguinte maneira: nobj

(4.4)

∑β k =1

j ,k

= 1 , 0 ≤ β j ,k ≤ 1

(4.4)

e f 0 k é a função objetiva k no projeto inicial x0. O algoritmo desse método pode ser representado pelos seguintes passos: i) Definir o número de subconjuntos β; ii) Normalizar as funções objetivo; iii) Para cada β j faça: iii.a) Obter a função objetivo substituta usando a Equação (4.3); iii.b) Otimizar a função substituta e encontrar o ponto xj*; iii.c) Substituir o xj* nas funções objetivo e obter os seus valores; Problemas na obtenção de pontos de Pareto via WS poderão surgir quando o contorno da região viável no espaço das funções objetivos for não-convexa, como mostra a Figura 4.4. Neste caso, não existirá nenhum β j capaz de fornecer uma solução que esteja na parte não-convexa. f2 f1*

f2* f1 Figura 4.4 Região viável não-convexa no espaço das funções objetivo. Em geral, quando se utiliza essa metodologia, ocorre que uma distribuição uniforme dos pesos β não fornece uma distribuição uniforme de pontos de Pareto. 55

4.3.2 Método Min-Max

Este método foi desenvolvido com o propósito de sanar os problemas listados anteriormente no método da soma ponderada dos objetivos, o que é parcialmente alcançado no que se refere a obter pontos uniformemente distribuídos. Neste procedimento, as funções objetivo são normalizadas, de forma diferente daquela do método da soma ponderada dos objetivos, sendo para isso acrescido dois parâmetros: max fk e min fk. Estes parâmetros são obtidos das soluções das otimizações escalares, isto é, que envolvem a consideração individual das funções objetivo isoladas. Obtém-se, então, um conjunto de variáveis xk*, resultante de cada otimização k isolada. Aplica-se esse conjunto a cada função objetivo e daí atribui-se o máximo valor da função a max fk e o mínimo a min fk. As novas funções objetivo normalizadas assumirão a forma: (4.5) f k =

f k − min f k max f k − min f k

, k = 1,..., nobj

(4.5)

O novo vetor das funções objetivo é (4.6) F = ⎡ f1 , f 2 ,..., f k ,..., f nobj ⎤ ⎣ ⎦

(4.6)

Então, um coeficiente βk é associado a cada nova função objetivo e o seguinte problema é proposto: (4.7) min ( γ )

(4.7)

onde (4.8) γ = max( β k f k ),

k = 1,..., n obj

(4.8)

sujeita às restrições definidas na Equação (4.2) e às restrições adicionais (4.9) β k f k ≤ γ

k = 1,..., nobj

(4.9)

para vários conjuntos de vetores β, pois para cada vetor β diferente um novo subproblema de otimização é formulado e um novo ponto de Pareto é obtido.

4.3.3 Método da Interseção Contorno-Normal (NBI)

O método da Interseção Contorno-Normal ou Normal-Boundary Intersection (NBI) (DAS e DENNIS, 1996) é uma técnica criada para encontrar pontos eficientes (ou pontos NBI) do contorno do espaço viável gerado pelos vetores objetivos alcançáveis, {F(x): x∈C}, que possibilitem a construção de uma curva suave, de forma que o projetista possa definir em qual daqueles pontos será considerada a solução compromisso para o problema multiobjetivo. Quando os pontos eficientes estiverem sobre uma parte 56

do contorno suficientemente convexa daquele espaço viável, esses são definidos como pontos de Pareto. Isto acontece para a grande maioria dos casos estudados na engenharia. Porém, se aqueles pontos estiverem na parte côncava do contorno, não há a garantia de que eles sejam pontos de Pareto. Apesar disso, esses pontos contribuem para que a curva Pareto seja definida. Este método é inovador com relação à obtenção dos pontos eficientes sobre a superfície Pareto, permitindo que se obtenha uma distribuição uniforme daqueles pontos até mesmo para um pequeno conjunto de vetores do parâmetro β, independentemente do número de funções objetivos. Este método difere daqueles descritos anteriormente por não ocorrer a transformação do vetor das funções objetivo em uma única função objetivo, através de uma combinação linear. A idéia central do NBI é encontrar uma porção do contorno do denominado espaço das funções objetivo (DAS e DENNIS, 1996), o qual contém os pontos ótimos de Pareto. Tais pontos podem ser encontrados resolvendo-se um problema de otimização. No que se segue, apresentam-se, inicialmente, algumas terminologias específicas do método para em seguida detalharmos a metodologia. Define-se F* como sendo o vetor do mínimo local das funções objetivo, denominado de Ponto Utópico (Shadow Minima ou Utopia Point) (DAS e DENNIS, 1996), representado por: (4.10) F* = ⎡⎣ f1* , f 2* , f 3* ,… , f nobj * ⎤⎦

(4.10)

T

onde cada f i* representa um mínimo local individual. Sendo o vetor xi* a solução ótima de f i , temos que f i * = f i ( xi* ) . Define-se a envoltória convexa do mínimo individual (ECMI) como: n ⎧ ⎫ (4.11) ⎨Φβ :β ∈ ℜ n , ∑ β i =1, β i ≥ 0⎬ i =1 ⎩ ⎭

(4.11)

onde (4.12) Φ (i, j ) = fi ( x*j ) − f i * , i = 1, … , nobj;

j = 1, … , nobj

(4.12)

Assim, os pontos pertencentes a ECMI são definidos por um conjunto de pontos do ℜn, que são definidos pelas combinações convexas de {F(xi*)- F*}, armazenados sob a forma de matriz, Φ, denominada de "pay-off" (DAS e DENNIS, 1996). Um exemplo da representação gráfica da ECMI é ilustrada na Figura 4.5. Nesta figura é considerado que na origem esteja o ponto de utopia F* e, dessa forma, todas as funções objetivos são não-negativas, isto é, F(x) é substituída por F que é definido da seguinte forma: F (x) = F(x) – F*

57

f2

A

C

B

f1

0 Figura 4.5 Representação gráfica da ECMI num espaço bidimensional. Com esta redefinição, observa-se na Figura 4.5 que o ponto A é F (x1*) e o ponto B é F ( x2*), 0 é a origem e ao mesmo tempo o ponto Utópico F*, o segmento tracejado é a ECMI, enquanto que o arco ACB é a fronteira de Pareto no espaço das funções objetivo. O conjunto das funções objetivo no espaço viável {F(x): x∈C} será denominado por ƒ e o seu contorno por ∂ƒ. A técnica NBI possui o objetivo de encontrar parte do contorno ∂ƒ que contém os pontos ótimos de Pareto. A idéia geométrica associada ao método é que tais pontos de Pareto são encontrados a partir da interseção da reta quase-normal à ECMI, apontada para a origem, e o contorno ∂ƒ como ilustrado na Figura 4.6. Nesta, observa-se que a família dos vetores quase-normais uniformemente espaçados intercepta os pontos igualmente espaçados sobre o contorno. f2

F ( x1*)

ECMI

∂ƒ

F ( x2 *)

f1

Figura 4.6 A imagem do conjunto viável sobre o mapeamento de ƒ no espaço das funções objetivo. Matematicamente, tais pontos podem ser encontrados resolvendo-se um problema de otimização descrito a seguir. Dados os parâmetros β, Φβ representa pontos sobre a ECMI. Seja nˆ o vetor unitário quase-normal à ECMI, i.e. direção que liga o ponto média da ECMI e o ponto Utópico F*. Então, Φβ + tnˆ , com t∈ℜ, representa o conjunto de pontos sobre nˆ , que formam uma reta quase-normal à ECMI. A interseção entre a

58

reta quase-normal a ECMI, e o contorno que define o espaço {F(x) | x∈C}, ∂ƒ, mais próximo da origem é a solução global do seguinte problema de programação não-linear: (4.13) max t x,t

(4.13)

sujeita às restrições definidas na Equação (4.2) e às restrições adicionais

ˆ = F ( x) (4.14) Φβ + tn

(4.14)

sendo esta equação de restrição a garantia do mapeamento de x por F ( x ) sobre a reta quase-normal. O problema apresentado nas Equações (4.13), (4.14) e (4.2) passa a ser definido com um subproblema NBI, representado por NBIβ, considerando que β seja o parâmetro que caracteriza o subproblema. Cada vetor parâmetro β, resultará em uma solução pretendente à Pareto. Resolvendo esse subproblema para um conjunto de parâmetros β, encontra-se um conjunto de pontos sobre ∂ƒ, que poderão fornecer uma curva suavizada. Esses pontos serão pontos de Pareto de ótimo caso estejam numa porção convexa do ∂ƒ. Em caso contrário, eles poderão não ser ponto de Pareto de ótimo. Em não sendo, poderão ser úteis na suavização da curva Pareto, além da garantia de que todo ∂ƒ foi explorado, o que não ocorre para os métodos descritos anteriormente. Para maiores detalhes, ver referências (DAS e DENNIS, 1996; MACEDO, 2002).

4.3.4 Método da Restrição Normal Normalizada (NNC)

Este método foi introduzido por (MESSAC et al, 2003) e é uma melhoria sobre o método da Restrição Normal (NC) de Yahaya e Messac removendo problemas numéricos de escala com a normalização das funções objetivo. Este método mostrou ser capaz de gerar uma propagação uniforme de pontos Pareto. O NNC trabalha de forma similar ao método da Interseção Contorno-Normal (NBI), (discutido anteriormente), e sua representação gráfica, tirada da referência (MESSAC et al, 2003), é mostrada a seguir. No espaço das funções objetivo normalizadas (mesma normalização usada no método Min-Max), todos os pontos mínimos locais das funções objetivo, estão a uma unidade de distância do Ponto Utópico. E o Ponto Utópico está na origem – por definição Na Figura 4.7 nós observamos o espaço viável e a fronteira de Pareto correspondente. Nós também notamos os dois pontos de mínimo das funções objetivo que são obtidos pela minimização separada do primeiro e do segundo objetivo do projeto. Uma linha ligando os dois pontos é chamada de Linha Utópica, ela é equivalente à envoltória convexa do mínimo individual (ECMI) no método NBI. Ela é dividida em m-1 segmentos, resultando em m pontos.

59

f2

F( x1* )

1

Nk m −1 F ( x2* )

0

1

0

f1

Figura 4.7 Um conjunto de pontos espaçados na linha Utópica para um problema bi-objetivo. Na Figura 4.8 um ponto genérico X interceptando a NU é utilizado para definir sua reta normal. Está linha normal (à NU) diminui a região viável, como indicado na Figura 4.8. Como pode ser visto, minimizando-se f 2 , o resultado ótimo será o ponto f ∗ . Transladando o ponto genérico X , pela Linha Utópica, podemos ver que um conjunto de Pontos Pareto correspondente será obtido.

f2

1

f 1∗ Região X pj

0

Linha NU

Viavel

f 0

∗

Região Inviavel

f

2∗

1

Plano Utópico

f1

Figura 4.8 Representação gráfica do método da NNC para um problema bi-objetivo. No que se segue, a formulação matemática do Método da Restrição Normal Normalizada para a otimização multiobjetivo em geral é apresentada. O vetor x define as variáveis de projeto e fi (f ∈ Rn ) a i-esima função objetivo. x* ( xi* ∈ R nx ) é o ponto ótimo da função fi . Sendo nobj o número de funções objetivo de um problema genérico, os * vetores do Plano Utópico N k são definidos como a direção dos F( xk* ) até F( xnobj ) , para k=1..nobj-1, pela equação:

* * (4.15) N k = F( xnobj ) − F( xk ),

k = 1, 2,...nobj − 1

(4.15)

60

ma.

Um conjunto de np pontos distribuídos no Plano Utópico é gerado da seguinte for-

nobj

(4.16) X j = ∑ β j ,k F ( xk* ) , k =1

j = 1, 2,...n p

(4.16)

onde (4.17) 0 ≤ β j ,k < 1 e

nobj

∑ β j ,k = 1,

k =1

j = 1, 2,..., n p

(4.17)

Agora nós geramos pontos distribuídos no espaço das funções objetivas normalizadas. Para cada ponto X gerados anteriormente é obtido um Ponto Pareto correspondente, solucionando o problema abaixo:

(4.18)

Para j = 1, 2,..., np

(4.18)

f nobj min x

sujeita às restrições definidas na Equação (4.2) e às restrições adicionais: (4.19) N k ( F − Xj ) ≤ 0, T

k = 1, 2,..., n − 1

(4.19)

4.3.5 Diferenças entre o NBI e o NNC As diferenças entre o método da Interseção Contorno-Normal e o método da Restrição Normal Normalizada são três: 1. Enquanto o NNC faz a normalização das funções objetivo, o NBI faz apenas uma translação. 2. O NNC diminui a região viável, do problema original, à uma região limitada (na direção de fn ) pela reta NU (normal ao Plano Utópico). Enquanto que no NBI a região viável se restringe estritamente à reta (quase) normal á ECMI. 3. Com relação à reta normal à ECMI (no NBI) ou ao Plano Utópico (no NNC), enquanto no NNC ela é na direção perpendicular, no NBI ela é quase normal (direção de um ponto médio da ECMI ao Ponto utópico). Sem estas diferenças, os resultados seriam absolutamente os mesmos.

61

4.4 EXEMPLOS 4.4.1 Treliça de 3 barras Neste primeiro exemplo consideramos a treliça de 3 barras de Koski (1985), ilustrada na Figura 4.9, na qual L = 1 m. Para este exemplo em particular, a carga aplicada é F = 20 KN e o módulo de Young é E = 200 GPa.

Figura 4.9 Treliça de três barras: definição do problema. As variáveis de projeto são as áreas das seções transversais das barras e os limites superiores e inferiores são respectivamente XL = 0.1 cm2 e XU = 2 cm2. Neste problema, soluções multiobjetivo serão encontradas considerando a minimização de duas funções objetivas: (1) Volume total da estrutura, (2) Combinação linear dos deslocamentos: dp = 0.25δ1 + 0.75δ2. Além das restrições das variáveis de projeto, são impostas restrições de tensão para todas as barras σall < 200 Mpa (para compressão e tração). Na seqüência 100 soluções de Pareto para o problema são apresentadas, considerando os métodos WS, Min-Max, NBI e NNC, respectivamente nas Figura 4.10 (a), (b), (c) e (d). Os resultados estão em acordo com os reportados na literatura (KOSKI, 1985, MESSAC et al., 2003). Como esperado, os métodos NBI e NNC produziram pontos de Pareto melhores distribuídos. Além disso existe uma região côncava na fronteira de Pareto como indicado na Figura 4.10 (a). É importante enfatizar que as técnicas Min-Max, NBI e NNC são capazes de calcular pontos nesta área ao contrário do método WS, que assume uma combinação convexa entre as funções. Na região côncava existe alguns pontos não Pareto, o método NBI encontrou esses pontos, como apresenta a Figura 4.10 (c). Apesar disso, (DAS e DENNIS, 1996) argumentam que esses pontos contribuem para que a curva Pareto seja definida. O método Min-Max definiu bem a curva de Pareto, desconsiderando corretamente a parcela não 62

Pareto da região convexa, Figura 4.10 (b). A distribuição de pontos do NNC excluiu alguns deles, como observado na Figura 4.10 (d). Porém todos os pontos não Pareto podem ser facilmente separados através de um filtro nos pontos.

a) WS

b) Min-Max

c) NBI

d) NNC

Figura 4.10 Treliças de três barras – Distribuição de Pareto: (a) WS, (b) Min-Max, (c) NBI e (d) NNC.

4.4.2 Placa quadrada com um orifício central O segundo exemplo se refere ao já apresentado na seção 2.5.1 e considerado também, na seção 3.6. Nele foi considerada a análise do estado plano de tensão em uma placa quadrada com um orifício central. O módulo de elasticidade é considerado o mesmo para todas as regiões: E = 105 N/mm. As dimensões do orifício são as variáveis do projeto selecionadas para a otimização. Os valores iniciais das variáveis de projeto são μ1 = μ2 = 50mm, e os limites infe63

rior e superior impostos são 25mm e 75mm, respectivamente. Dois objetivos são considerados, são eles: minimizar o volume total e minimizar a flexibilidade total da placa. As soluções MO serão obtidas considerando 10 pontos de Pareto. Para a aproximação via MBR, três regiões são definidas. A base reduzida será construída no espaço viável do projeto D = [25, 75]² e o número de amostras analisadas foi N = 10, baseada no estudo de convergência do MBR, já realizado para este caso, na seção 2.5.1. O domínio de referência (o qual neste caso é igual ao projeto inicial) com sua respectiva malha de EF para elementos triangulares CST, com número total de graus de liberdade – ndof = 6.732, é mostrado na Figura 4.11 (a). A distribuição das tensões de von Mises e a configuração deformada para o projeto inicial é mostrada na Figura 4.11 (b).

100

80

60

40

20

0

0

20

40

60

80

100

Figura 4.11 Projeto inicial (domínio de referência). a) malha EF (ndof = 6.732), b) tensão de von Mises (N/mm²). A Figura 4.12 apresenta as distribuições dos pontos Paretos obtidos para os métodos aqui considerados. As soluções são obtidas considerando o modelo de MEF e o modelo aproximado via MBR. Verifica-se que as curvas de Pareto via MEF e MBR estão de acordo. Como pode ser observado, soluções via o método da soma ponderada (WS) são bastante pobres com pontos aglutinados (cluster) e regiões vazias sobre a fronteira de Pareto esperada. Soluções via Min-Max se mostraram melhores. Entretanto os métodos NBI e NNC são os que conseguem obter pontos uniformemente espaçados em todas as partes da fronteira de Pareto.

64

3

3 FEM RBM

FEM RBM

2.5 Strain Energy (N.mm)

Strain Energy (N.mm)

2.5 2 1.5 1 0.5

2 1.5 1 0.5

0 4000

5000

6000 7000 8000 Volume (mm³)

9000

0 4000

10000

5000

6000 7000 8000 Volume (mm³)

a)WS

10000

b)Min-Max

3

3 FEM RBM

FEM RBM

2.5 Strain Energy (N.mm)

2.5 Strain Energy (N.mm)

9000

2 1.5 1 0.5

2 1.5 1 0.5

0 4000

5000

6000 7000 8000 Volume (mm³)

9000

0 4000

10000

5000

6000 7000 8000 Volume (mm³)

c)NBI

9000

10000

d)NNC

Figura 4.12 Placa quadrada com um orifício central – Pontos Paretos: a) WS, b) Min-Max, c) NBI, d) NNC. A Tabela 4.1 sumariza as performances de cada método investigado. Soluções via MBR são mais que duas ordens de magnitude mais rápidas em comparação com os resultados totalmente obtidos via MEF. Para a obtenção de pontos de Pareto pelos procedimento considerados, utilizando o MBR, o estágio “off-line”, foi avaliado apenas uma vez, totalizando 37.4 s. Além disso, o processo de obtenção dos pontos de Pareto (estágio “on-line”) pelos quatro métodos, foram feitas em paralelo em um computador dotado de um processador quad core, através de ferramentas do ambiente computacional MATLAB. Um ganho computacional ainda maior seria alcançado paralelizando, também, o estágio off-line. Tabela 4.1 Placa quadrada com um orifício central – Desempenho dos algoritmos. WS FEM RBM (37.4 s)

Min-Max

NBI

NNC

4902.0s

3846.2s

3444.3s

4613.7s

21.5s

5.2s

7.3s

7.8s

Na Figura 4.13 os projetos ótimos (na configuração deformada) para os pontos de Pareto 1, 3, 5, 7, 9 e 10 são ilustrados, para dar uma idéia dos pontos de Pareto obtidos pelo método NBI.

65

Ponto de Pareto - 1

Ponto de Pareto - 3

Ponto de Pareto - 5

Ponto de Pareto - 7

Ponto de Pareto - 9

Ponto de Pareto - 10

Figura 4.13 Placa quadrada com um orifício central – Tensões von Mises (N/mm²) na configuração deformada, dos projetos para pontos de Pareto ótimos obtidos pelo método NBI.

4.4.3 Placa engastada - problema aclopado Finalizando as aplicações considerando duas funções objetivo, será estudada uma estrutura plana bidimensional, a mesma analisada na seção 2.5.2 e indicada na Figura 4.14 (a), porém com condições de contorno diferentes, indicada na Figura 4.14 (b). Tem-se um total de treze (13) regiões especificadas para o presente problema, numeradas conforme indicado na Figura 4.14 (a) em função dos parâmetros variantes do problema em questão. Condições de contorno para o problema estático e de transferência de calor estão indicadas na Figura 4.14 (b), onde a lateral esquerda está engastada e a parte superior está submetida a um carregamento distribuído de 1MPa, ambas regiões, assim como o interior vazio da estrutura, encontram-se isoladas termicamente. Todas as dimensões do problema estão especificadas na figura. p=1 MPa

q = 1 W/m² Ta = 20 ºC Figura 4.14 Definição do problema: a) Geometria , b) Condições de contorno. 66

As propriedades do material são E= 10³ MPa, o coeficiente de Poisson υ = 0.2, coeficiente de convecção h = 10 W/m² ºC, a condutividade térmica κx = κy = 1 W/(mK), o coeficiente de dilatação térmica α = 10-4 ºC-1 e as condições de convecção são para temperatura ambiente Ta = 20ºC. Quatro parâmetros (variáveis de projeto), indicados na Figura 4.14 (a) serão considerados na otimização. São eles: as espessuras t1, t2 e as distâncias horizontais L1, L2 indicadas na figura acima. A aproximação foi construída de forma que os parâmetros μ=( t1, t2, L1, L2 ) pertençam ao espaço de projeto Dμ = [0.5;5]2 × [10,21]× [28,38] . Para o presente problema o número de amostras consideradas foi N=30, baseado no estudo realizado no capítulo 2, seção 2.5.2. A discretização do domínio real (projeto inicial) onde μ1 = 2, μ2 = 2, μ3 = 19, μ4 = 31 e o domínio de referência (1976 graus de liberdade para o problema elástico e 1023 para o térmico) se encontram, respectivamente, na Figura 4.15 (a) e (b). A otimização será conduzida considerando os seguintes objetivos: Minimização da temperatura máxima e minimização da tensão Von Mises máxima. Além dos limites das variáveis, o volume inicial deve ser mantido constante. As soluções OM serão obtidas considerando 15 pontos de Pareto.

a) domínio real

b) domínio de referência

Figura 4.15 Malha de elementos finitos: a) domínio real (projeto inicial), b) domínio de referência. A Figura 4.16 apresenta a distribuição de Pontos Pareto para os métodos aqui considerados. As curvas de Pareto obtidas pelo MEF e MBR estão de acordo. Como podese observar, as soluções via o método WS e o Min-Max apresentam regiões vazias na fronteira de Pareto esperada, o WS ainda apresenta pontos sobrepostos. Já os métodos NBI e NNC conseguiram obter pontos bem distribuídos em todas as partes da fronteira de Pareto.

67

a)WS

b)Min-Max

c)NBI

d)NNC

Figura 4.16 Pontos de Pareto: a) Soma Ponderada, b) Min-Max, c) NBI, d) NNC. A Tabela 2 resume as performances de cada procedimento investigado. Soluções via MBR foram mais rápidas comparadas com o procedimento baseado apenas no MEF. O ganho computacional do MBR não foi tão grande quanto o obtido no exemplo anterior, devido à computação do termo de acoplamento, realizado na etapa on-line. Tabela 2. Desempenho dos algoritmos. WS

Min-Max

NBI

NNC

FEM

777,12s

826,92s

1239,30s

963,00s

RBM

340,53s

387,46s

353,25s

323,51s

4.5 PROBLEMAS COM MAIS DE DUAS FUNÇÕES OBJETIVO Pode-se obter qualquer ponto de Pareto de um problema qualquer, tanto pelo NBI como pelo NNC, resolvendo um subproblema com um determinado vetor β associado (Equações (4.14), (4.16) e (4.19)). Porém no caso de problemas com mais de duas funções objetivo, as componentes do vetor β não são necessariamente não-nulas. Encon68

trar pontos sobre toda a fronteira de Pareto em problemas com muitas funções objetivo ainda é um problema em aberto (ARORA, MESSAC e MULLUR, 2007). Como exemplo, considere a Figura 4.17 (a) (DAS e DENNIS, 1996) onde é apresentado o espaço de três funções objetivo de um problema qualquer. Caso o vetor β só tenha componentes positivas, as Equações (4.11) e (4.16) indicam que apenas é possível encontrar pontos partindo do interior do triângulo formado por F(X1*), F(X2*) e F(X3*) (vide Figura 4.17 (b)), desprezando deste modo, a região da fronteira de Pareto adjacente, incluindo os arcos F(X1*)F12*F(X2*), F(X1*)F13*F(X3*) e F(X2*)F23*F(X3*). Além disto, os pontos podem recair sobre uma região fora da superfície de Pareto, é o que ocorre neste caso. Como podemos observar na Figura 4.17, uma parte do segmento de reta F(X1* ) F(X3* ) , se projeta no exterior da região de Pareto. Neste caso o NBI irá encontrar um ponto não-Pareto, devido a sua restrição rígida (de igualdade). Já o ótimo encontrado através do NNC, pode se deslocar para a região de Pareto, devido à suas restrições mais flexíveis (de desigualdade). A Figura 4.17 (b), apresenta os pontos θβ (Equação (4.11)) na ECMI (para o NBI) ou os pontos X j (Equação (4.16)) no plano utópico (para o NNC), para coeficientes β normalizados igualmente espaçados. Nota-se que os métodos só se preocupam com a região formada pelos pontos F(X1*), F(X2*) e F(X3*).

F12* F13*

F23*

a) Frente de Pareto

69

b) Pontos base para os sub-problemas Figura 4.17 Problema com 3 funções objetivo (DAS e DENNIS, 1996), a) Frente de Pareto no espaço de três funções objetivo, b) Pontos base para a definição dos sub-problemas de otimização dos métodos NBI e NNC. Para encontrar toda a fronteira de Pareto, (MESSAC e MATTSON, 2004) apresentam uma modificação no método NNC original, que visa solucionar o problema das regiões inexploradas pelo método original, através da ampliação do Plano Utópico (ou a ECMI) de forma que ele contenha toda a superfície de Pareto. Deste novo Plano Utópico são desconsideradas algumas partes onde há menos probabilidade em haver pontos Pareto projetados. Porém o problema das otimizações em regiões onde não há pontos de Pareto projetados continua e ainda é amplificado, podendo ainda os pontos se projetarem em regiões inviáveis, gerando um problema de otimização mal formulado. Quanto mais pontos não-Paretos encontrados, menos eficiente se torna o procedimento, mesmo retirando-os por um processo de filtragem, pois é uma otimização desperdiçada. Outra possível solução, proposta no presente trabalho, é definir inicialmente os contornos da fronteira de Pareto e encontrar os pontos no seu interior. Para isso é preciso resolver o problema multiobjetivo por pares de funções objetivo e assim definir cada contorno da fronteira de Pareto, ou seja, resolvendo o POM considerando apenas os objetivos f1 e f2, obtemos a curva de Pareto F12*, resolvendo o POM para f2 e f3, encontramos a curva de Pareto F23*, já desconsiderando a função f2, alcançamos a curva de Pareto F13*. Na Figura 4.18 vemos um exemplo de pontos de Pareto encontrados para cada par de funções objetivo marcados em preto, a partir destes pontos de Pareto são distribuídos os restantes dos pontos centrais, marcados de azul.

70

Figura 4.18 Problema com 3 funções objetivo – Distribuição dos pontos. Depois de encontrados os pontos de Pareto que formam o contorno da fronteira total (pontos pretos na Figura 4.18), projeta-se as curvas formadas por estes pontos, no plano formado por F(X1*), F(X2*) e F(X3*). Estas curvas projetadas delimitam uma região a qual é utilizada (ao invés do triângulo original) para encontrarmos o restante dos pontos centrais. As curvas F12*, F23* e F13* são consideradas as soluções de Pareto dos subproblemas onde o coeficiente β da função desconsiderada é nulo. Um exemplo desta metodologia é tratada a seguir a partir da seção 4.5.2 onde será utilizado o NBI para se encontrar as curvas de Pareto do contorno da superfície de Pareto, bem como para as otimizações dos subproblemas relacionados aos pontos centrais, este procedimento proposto, será chamado NBIm. Para problemas com mais de 3 funções objetivo a idéia básica é a mesma, 1. Resolver os problemas de otimização uni-objetivo para cada objetivo. 2. Em posse dos mínimos individuais, que limitam as curvas de Pareto, resolver os problemas bi-objetivos para todas as possíveis permutações de objetivos. 3. Em posse das curvas de Pareto para os pares de objetivos, que delimitam as superfícies de Pareto, resolver os problemas com três objetivos para todas as possíveis permutações. 4. Em posse das superfícies de Pareto para os trios de objetivos, que delimitam as hiper-superfícies de Pareto, resolver os problemas com quatro objetivos para todas as possíveis permutações. E assim sucessivamente até que se chegue ao número total de objetivos.

71

4.5.1 Exemplo geométrico com 3 funções objetivo Considere o problema definido na Figura 4.19, a qual ilustra o espaço das (duas) variáveis de projeto. Serão considerados três objetivos, são eles minimizar a distância do ponto definido pelas variáveis, aos pontos A, B e C. Tem-se ainda a restrição que o espaço viável exclui o circulo preto indicado na Figura 4.19, ou seja, a distância do ponto definido pelas variáveis, ao ponto D deve ser maior que 0.5.

Figura 4.19 Definição do problema. A Figura 4.20 (a) ilustra a região ótima exata dos pontos de Pareto, indicada pela cor azul. As figuras seguintes apresentam os resultados obtidos para as diferentes técnicas, onde cada ponto azul indica o valor da variável de projeto para cada subproblema otimizado. Na legenda de cada figura, a partir da letra (b), se encontram o método utilizado e, entre parênteses, o número total de funções avaliadas, utilizadas pelo método. Nota-se, novamente, os espaçamentos vazios entre as soluções do método da soma ponderada e o grande numero de funções calculadas. O Min-Max já apresenta um resultado melhor apresentando uma boa definição dos contornos da região de ótimos e um numero bem menor de funções avaliadas, porém com poucos resultados centrais. O NNC e o NBI, apresentam pontos bem distribuídos, porém apresentando soluções nãoPareto, principalmente o NBI, por motivos já citados. Ambos os métodos também realizaram um grande número de avaliações de funções, isto ocorreu principalmente, devido a procura de pontos viáveis em sub-problemas mal formulados sem região viável. O NBIm, foi o que melhor cobriu a região de ótimos, com pontos relativamente bem distribuídos e contorno definido, foi também, o que utilizou o menor número de avaliações de funções entre os métodos investigados.

a) Solução exata do problema

b) WS (9263) 72

c) Min-Max (3802)

e) NNC (6819)

d) NBI (7553)

f) NBIm (2048)

Figura 4.20 Soluções de Pareto: a) Solução exata, b) Soma Ponderada (9263), c) Min-Max (3802), d) NBI (7553), e) NNC (6819), f) NBIm (2048).

4.5.2 Exemplo analítico com 3 funções objetivo Considere o seguinte problema de otimização multiobjetivo: min f ( x)i = xi , i = 1, 2,3 x

Sujeito às seguintes restrições : x1 ≥

1 1 1 1 1 1 + ; x2 ≥ + ; x3 ≥ + x2 x3 x1 x3 x1 x2

0.2 ≤ xi ≤ 10, i = 1, 2,3 A Figura 4.21 apresenta os pontos de Pareto para 15 valores de βi , i = 1..3 uniformemente distribuídos, que combinados totalizam 120 vetores β diferentes, para mais detalhes sobre as combinações dos parâmetros β vide (DAS e DENNIS, 1996). O gráfico da esquerda apresenta uma vista de topo ( f ( x)1 x f ( x) 2 ) do espaço das funções objetivo, onde a escala de cor se refere ao valor da terceira função objetivo f ( x)3 . O gráfico da direita apresenta uma vista isométrica do espaço das funções objetivo. Na legenda de cada par de figura, se encontram o método utilizado e , entre parênteses, o número total de funções avaliadas.

73

a) WS (7092)

b) Min-Max (6090)

c) NBI (3281)

74

d) NNC (5617)

e) NBIm (2748) Figura 4.21 Pontos de Pareto: a) Soma Ponderada (7092), b) Min-Max (6090), c) NBI (3281), d) NNC (5617), e) NBIm (2748). Como se pode observar, os métodos da Soma Ponderada e o Min-Max apresentaram concentrações de pontos em alguns lugares, espaços vazios em outros e os maiores números de funções calculadas. Já os métodos NNC e o NBI obtiveram pontos igualmente espaçados, porém não foram capazes de cobrir toda a superfície de Pareto, apenas as regiões onde se projeta a ECMI (triangulo formado pelos mínimos individuais). Entre estes, o NBI apresentou um numero menor de funções avaliadas. Já com o NBIm, através da estratégia de definir os contornos da região de Pareto, foram obtidos pontos bem distribuídos sobre toda a região de Pareto, apresentando ainda o menor número de funções avaliadas entre os métodos investigados.

4.5.3 Placa sob ação termo-estrutural acoplada com 3 funções objetivo Como exemplo, foi considerada a placa engastada do exemplo 4.4.3 com condições de contorno, para o problema estático e de transferência de calor, diferentes, mostradas na Figura 4.22. 75

Figura 4.22 Definição do problema - Condições de contorno. As propriedades do material são: E= 104 MPa, o coeficiente de Poisson υ = 0.2, coeficiente de convecção h = 10 W/m² ºC, a condutividade térmica κx = κy = k = 1 W/(mK), o coeficiente de dilatação térmica α = 2x10-4 ºC-1 e as condições de convecção são para temperatura ambiente Ta = 20 ºC e uma temperatura externa de 100 ºC na lateral direita da estrutura (o sentido da convecção indicado na figura é meramente ilustrativo). Os mesmos quatro parâmetros (variáveis de projeto), indicados na Figura 4.14 (a) serão considerados na otimização. Para o presente problema o número de amostras consideradas foi N=30, vide estudo de convergência seção 2.5.2. A otimização será conduzida considerando três objetivos: Minimizar o volume da estrutura, minimizar a flexibilidade total e minimizar o deslocamento máximo. Além dos limites das variáveis, a tensão Von Misses foi limitada a 124 kN/cm². As soluções OM serão obtidas considerando 55 pontos de Pareto. A Figura 4.23 apresenta a distribuição de Pontos Pareto para os métodos MO aqui considerados obtidas pelo MBR. Nos resultados obtidos, foram filtrados os pontos dominados (não Pareto), esse filtro verifica para toda solução se existe algum outra que a domina. Os “x” marcados de vermelho na Figura 4.23, são pontos não Paretos, pois foram dominados por alguma solução de algum dos métodos. Como podemos observar na Figura 4.23, a solução via o método WS e o Min-Max são ineficientes, pois apresentam pontos sobrepostos e regiões vazias na superfície de Pareto. Assim como, os métodos NBI e NNC apresentaram resultados bem diferentes, como o NBI não tem flexibilidade nas restrições adicionais, ele não foi capaz de encontrar a verdadeira superfície de Pareto, apenas alguns pontos, pois ela não se encontra totalmente na projeção da ECMI. O NNC se comportou pouco melhor, porém não conseguiu definir toda a superfície de Pareto e não consegui obter um resultado melhor que o WS. Já o NBIm foi o que conseguiu definir melhor todas as partes da fronteira de Pareto, seus resultados podiam ser ainda melhor através de melhorias no procedimento de distribuição dos pontos centrais (estudo ainda não concluído).

76

a) WS

b) Min-Max

c) NBI

d) NNC

e) NBIm Figura 4.23 Pontos de Pareto: a) Soma Ponderada, b) Min-Max, c) NBI, d) NNC, e) NBIm.

4.6 REFERÊNCIAS AFONSO, S. M. B.; MACEDO, C. M. H.; OLIVEIRA, D. A. P. “Structural Shape Optimization under Multicriteria Conditions In: V World Congress on Computational Mechanics, Viena. 2002.

77

ARORA J. S.; MESSAC, A.; MULLUR, A. A. “OPTIMIZATION OF STRUCTURAL AND MECHANICAL SYSTEMS. Chapter 4 - Multiobjective Optimization: Concepts and Methods”. Jasbir S Arora, University of Iowa, USA, 2007. BATES, S. “Development of Robust Simulation, Design and Optimization Techniques for Engineering Applications”. PhD Thesis. School of Engineering, University of Wales, Swansea. 2003. COLLETTE, Y., SIARRY, P. “Multiobjective Optimization: Principles and Case Studies”, Springer, 2004. DAS, I.; DENNIS, J.E. Normal Boundary Intersection: A New Method for Generating Pareto Surface in Nonlinear Multicriteria Optimization Problems. SIAM J. Optimization, Vol. 8, No. 3, pp. 631-657, 1996. HWANG, C. L.; PAIDY, S. R.; YOON, K. e MASUD, A. S. M., 1980. Mathematical Programing whith Multiple Objectives: A Tutorial, Comput. and Ops. Res., Vol. 7, pp. 5-31. LALONDE, N.; KIM, I. Y.; WECK, O. “A Comprehensive Comparison between Deterministic and Probabilistic Multiobjective Optimization Algorithms with Mathematical and Practical Applications”. 8º World Congress on Structural and Multidisciplinary Optimization. Lisbon, Portugal, 2009. MACEDO, C. M. H. “Otimização de Treliças Planas sob Várias Solicitações com Ênfase a Problemas Multiobjetivos”. dissertação de Mestrado, Universidade Federal de Pernambuco. Recife-PE, Brasil, 2002. MESSAC, A., ISMAIL-YAHAYA, A. e MATTSON C. A. “The Normalized Normal Constraint Method for Generating the Pareto Frontier”. Structural Optimization, Vol. 25, No. 2, pp. 86-98, 2003. MESSAC, A.; MATTSON C.A. “Normal constraint method with guarantee of even representation of complete Pareto frontier”, 45th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics & Material Conference, Palm Springs, CA, 2004. STEUER, R. E. “Multicriteria optimization – theory, computation and application” . John Wiley & Sons, 1985.

78

79

CAPÍTULO 5 OTIMIZAÇÃO CONSIDERANDO INCERTEZAS 5 OTIMIZAÇÃO CONSIDERANDO INCERTEZAS

5.1 INTRODUÇÃO Este capítulo se concentra no projeto de sistemas estruturais na presença de incertezas. A motivação para tal vem do fato de que algum grau de incerteza ou variação nas características de qualquer sistema estrutural é inevitável. Quando se analisa e projeta-se um sistema estrutural, a caracterização determinística de suas propriedades e do ambiente que o cerca pode não ser oportuno por várias razões, incluindo incertezas nas propriedades dos materiais, variação na geometria devido à tolerância na fabricação, incerteza no carregamento devido à natureza não determinística do ambiente operacional, incertezas devido a degradação do sistema em serviço, e assim por diante. Na prática da engenharia é comum assumir valores nominais para os parâmetros incertos ao realizar estudos de projeto, de fato, a consideração determinística tem sido feita implicitamente em toda esta dissertação até o presente. Infelizmente, a abordagem determinística leva a um projeto final cujo desempenho pode cair significativamente devido a perturbações decorrentes de incertezas. Este problema é acentuado quando lida-se com projetos que foram severamente otimizados, pois projetos ótimos tendem a se localizar nos extremos da função objetivo ou no contorno das restrições (pequenas perturbações podem levar a queda de desempenho e/ou violação das restrições do projeto). Como conseqüência, um ótimo determinístico pode potencialmente ser uma solução de alto risco (BEYER et al, 2007). Marczyk escreveu “Otimização, na verdade é o oposto de robustez” (MARCZYK, 2000). Apesar de não ser de todo certo, há alguma verdade nisso, porém pode-se utilizar os atuais algoritmos de otimização, de forma apropriada para atingir projetos robustos. Uma maneira de garantir segurança contra incertezas é aplicar restrições mais rígidas que as idealmente impostas. Por exemplo, quando se projeta um sistema estrutural, para garantir que ele não falhe geralmente são impostas restrições de projetos nas tensões da forma σ (x) < σ max , onde σ max é a tensão de falha do material e x é o conjunto das variáveis de projeto. Para considerar as incertezas já citadas, a restrição pode ser escrita como FSσ (x) < σ max , onde Fs é chamado o fator de segurança. Normalmente o valor do fator de segurança varia entre 1.2 e 3.0. Este valor muitas vezes é definido com base na experiência a priori do material usado e de projetos com conceitos similares. Claramente, o novo projeto ótimo será mais conservador quanto mais se aumenta o valor de Fs, haja visto que o ótimo se distancia do contorno da restrição original σ (x) − σ max = 0 . Para o uso de novos materiais e projetos com novos conceitos, para os quais não há a priori experiência no projeto e limitada informação experimental, não é 80

trivial decidir um valor apropriado para os fatores de segurança (KEANE e NAIR, 2005). Neste capítulo, serão examinadas algumas abordagens para a consideração das incertezas no processo de otimização e assim obter projetos robustos. Isto contrasta com as abordagens que utilizam o fator de segurança, nos quais os parâmetros incertos não são explicitamente incorporados na formulação do projeto. Projetos robustos requerem que seu desempenho pouco se altere quando exposto à incertezas. Várias medidas de robustez foram propostas, incluído medidas de esperança e de dispersão. Em particular, a esperança e a variância são as medidas básicas usadas para otimização robusta (DOLTSINIS e KANG, 2004; BEYER et al, 2007; SCHUËLLER e JENSEN, 2008). Quando usa-se estas medidas, à busca por um projeto robusto ótimo, surge com um problema de decisão com múltiplos critérios. Por exemplo, ao otimizar a esperança pode se encontrar um valor alto de dispersão ou variância, o que em geral não é desejável. Nestes casos há uma “negociação” (trade-off) entre o valor médio e a variância, então uma solução combinada deve ser encontrada (otimização multiobjetivo robusta - OMR). Alternativamente, pode ser usada alguma das metodologias descritas no Cap. 4 para geração de um conjunto de soluções ótimas de Pareto que são possíveis candidatas à solução.

5.2 TEORIA PROBABILÍSTICA E ESTATÍSTICA Esta seção introduz conceitos elementares da teoria da probabilidade e conceitos estatísticos. Intuitivamente, a probabilidade é uma medida da freqüência de ocorrência de um determinado evento. Em um experimento isto pode ser medido dividindo o número de situações favoráveis pelo número de eventos possíveis. Entretanto, uma definição matematicamente mais rigorosa pode ser formulada a partir de três axiomas básicos. Considerando eventos relacionados aos subconjuntos A, B, C.... contidos num conjunto Ω , que é o conjunto de todos os eventos possíveis, a função Prob (i.e. probabilidade) é definida através de três axiomas:

(5.1) 0 ≤ Prob [ A] ≤ 1

(5.1)

(5.2) Prob [ Ω ] = 1

(5.2)

(5.3) Prob [ A ∪ B ] = Prob [ A] + Prob [ B ] − Prob [ A ∩ B ]

(5.3)

A probabilidade de um evento A condicionado a ocorrência de um evento B, é a probabilidade de ocorrer A dado que o evento B ocorreu, e pode ser definida como: (5.4) Prob ⎡⎣ A B ⎤⎦ =

Prob [ A ∩ B ] Prob [ B ]

(5.4)

81

5.2.1 Variáveis Aleatórias Variáveis aleatórias são funções de um espaço amostral que assumem valores numéricos através de um mapeamento do fenômeno aleatório associado. Elas podem ser de dois tipos: discretas ou contínuas. Variáveis aleatórias discretas podem assumir apenas um número finito de valores. Variáveis aleatórias contínuas podem assumir infinitos valores distintos, em geral são valores reais de alguma medida. Neste trabalho serão usadas apenas variáveis aleatórias contínuas reais. Para toda variável aleatória real, existe uma função chamada de Função Densidade de Probabilidade (“Probability Density Function” - PDF) P (ξ ) que define a distribuição de ocorrências de ξ ∈ associada a um fenômeno aleatório, tal que (MEYER, 1983): b

(5.5) Prob(a ≤ ξ ≤ b) = ∫ P(ξ )dξ a

∀a, b ∈

ea≤b

(5.5)

Um evento aleatório A pode ser definido como a ocorrência de uma determinada variável aleatória real ξ ser menor que um valor determinístico prescrito x. Assim, A = {ξ ξ < x} . A probabilidade Prob[A] associado com este evento obviamente depen-

de do valor prescrito x, i.e. Prob[A] = Fξ ( x) , esta função Fξ ( x) é chamada função de distribuição acumulada (“Cumulative Distribution Function” - CDF), calculada como: x

(5.6) Fξ ( x) = ∫ P(ξ )d ξ

∀x ∈

−∞

(5.6)

Se ξ é uma variável aleatória, então uma função qualquer f (ξ ) também será, porém com uma PDF própria associada. As PDFs são dependentes de parâmetros, os quais podem possuir interpretações práticas, tais como o seu valor médio ou média f (ξ ) e a sua variância υ. Matematicamente, a média em relação a uma função aleatória é chamada de esperança da função E [ f (ξ ) ] , calculado como (5.7) f = E [ f (ξ ) ] = ∫

∞

−∞

f (ξ ) P(ξ )d ξ

(5.7)

A variância ( σ 2f = σ [ f (ξ )] ) pode ser calculada como: 2

(5.8) σ 2f = E ⎡⎢( f (ξ ) − f ⎣

) ⎤⎦⎥ = ∫ ( f (ξ ) − f ) ∞

2

2

−∞

P(ξ )dξ

(5.8)

ou (5.9) σ 2f = E ⎡⎣ f (ξ ) 2 ⎤⎦ − ( f

)

2

= f 2 −( f

)

2

=∫

∞

−∞

( f (ξ ) ) P(ξ )dξ − ( f ) 2

2

(5.9)

82

onde σ f é a raiz quadrada da variância de f (ξ ) , conhecida como desvio padrão. A descrição das variáveis aleatórias em termos da média e do desvio padrão é chamada “representação de segundo momento”. Uma generalização desta representação é dada pela definição do momento central de k-ésima ordem (BUCHER, 2009). (5.10) μk = E ⎡( f (ξ ) − f ⎢⎣

) ⎤⎥⎦ = ∫ ( f (ξ ) − f ) ∞

k

k

−∞

P(ξ )dξ , k ≥ 2

(5.10)

Estes momentos serão úteis aqui, para o cálculo de duas outras estatísticas: a obliqüidade ("skewness") e a curtose (“kurtosis”). A obliqüidade é uma medida da assimetria de uma determinada distribuição, enquanto que a curtose é uma medida de dispersão que caracteriza o "achatamento" da curva da função de distribuição e são definidas respectivamente por:

(5.11) s =

μ3 σ3

e κ=

μ4 −3 σ4

(5.11)

5.2.2 Distribuições de probabilidade Neste trabalho serão utilizadas apenas dois tipos de distribuição, a distribuição Normal ou Gaussiana e a distribuição Lognormal. A distribuição Normal, também conhecida como Gaussiana, é uma das distribuições fundamentais da teoria estatística. O teorema do limite-central, diz que o somatório de variáveis aleatórias de distribuição quaisquer, tende a uma distribuição Normal (PAPOULIS, 1965). Ela descreve bem, vários fenômenos aleatórios “naturais”. É também bastante simples, pois é definida com apenas dois parâmetros básicos, a média e o desvio padrão. Se uma variável ξ tem uma distribuição Normal, escrevesse simbolicamente ξ ~ N ( μ , σ ) , onde μ é a média e σ o desvio padrão da variável. A PDF normal para uma variável aleatória ξ é expressa por:

⎛ ξ −μ ⎞ ⎡ ⎜ ⎟ (5.12) P (ξ ) = ⎢ 2πσ e⎝ σ ⎠ ⎢⎣

2

⎤ ⎥ ⎥⎦

−1/2

(5.12)

Sua formula padrão N ( 0,1) considera média igual a zero e desvio padrão unitário:

(

(5.13) P (ξ ) = 2π eξ

2

)

−1/2

(5.13)

83

A Figura 5.1 apresenta o gráfico da PDF normal padrão.

Função Densidade de Probabilidade 0.4 0.35 0.3

P(x)

0.25 0.2 0.15 0.1 0.05 0 -3

-2

-1

0 x

1

2

3

Figura 5.1 Função densidade de probabilidade (PDF) Normal. A distribuição Normal e a Lognormal estão fortemente relacionadas. Se uma variável aleatória ξ tem distribuição Normal com média μ e variância v , então exp( ξ ) tem distribuição Lognormal, com média μlog e variância vlog , respectivamente, iguais a: (5.14) μlog = e μ + v 2

vlog = e(

2 μ +v)

(e

v

− 1)

(5.14)

Esta distribuição é muito usada para variáveis aleatórias que só aceitam valores positivos, pois ela é definida apenas para valores reais maiores do que zero. Quando sua média se afasta de zero, com um valor fixo do desvio padrão, a distribuição Lognormal se aproxima de uma distribuição Normal. A PDF Lognormal é definida como: ⎛ − ln (ξ a )2 ⎞ 1 exp ⎜ (5.15) P (ξ ) = ⎟= ⎜ ⎟ 2 b ξ b2π ⎝ ⎠

onde b = ln(c 2 + 1) , a =

μ

ln (ξ / a ) ⎞ 1⎛ b 2 π ξ a b ( ) ⎜ ⎟ x⎝ ⎠

e c=

−1 2

(5.15)

σ

é o coeficiente variacional definido apenas X c +1 para variáveis aleatórias com médias não-zero. A Figura 5.2 apresenta o gráfico da PDF Lognormal com variância unitária para diferentes valores de média. 2

84

Pdf Lognormal com variância 1 1 Media 1 Media 2 Media 6

0.9 0.8 0.7

P(x)

0.6 0.5 0.4 0.3 0.2 0.1 0 0

2

4

6

8

10

x

Figura 5.2 Função densidade de probabilidade Lognormal para valores distintos de média.

5.3 CÁLCULO DAS ESTATÍSTICAS

Sendo ξ uma variável aleatória com distribuição de probabilidade conhecida, pode-se calcular diretamente a distribuição de probabilidade para uma função f ( X ) , caso f seja uma função monotônica, através da sua função inversa (BUCHER, 2009). Foi apresentado também como encontrar diretamente parâmetros estocásticos de f utilizando a distribuição de probabilidade da variável aleatória ξ , através de Equações (5.7) à (5.10). Porém, na grande maioria dos problemas envolvendo simulação numérica, estes procedimentos analíticos não são viáveis ou possíveis. Neste caso é necessário o uso de métodos aproximados. Há duais classes básicas de métodos usados no cálculo das estatísticas da função f (ξ ) , os métodos não-intrusivos (“Black-box”) e os métodos intrusivos (“Physicsbased”). Na análise de incertezas o objetivo central é ter um método geral que possa ser aplicado a sistemas complexos com parâmetros incertos. Grandes esforços têm sido aplicados no desenvolvimento de um método não-intrusivo para análise de incerteza, que explora modelos computacionais determinísticos existentes. Estas abordagens consideram o sistema computacional como uma caixa-preta (“Black-box”) que retornam o valor da função para um dado vetor de entradas. Neste trabalho serão utilizados dois métodos que satisfazem esta condição, o método de Monte Carlo (MC) e o Método da Colocação Probabilística (“Probabilistic Collocation Method” - PCM) (RAMAMURTHY, 2005; WEBSTER). Os problemas aqui tratados, essencialmente envolveram o cálculo apenas de alguns momentos estatísticos sem a necessidade de se encontrar toda a distribuição de probabilidade da função aleatória. Considere as Equações (5.7) e (5.9) a idéia básica é aproximar numericamente estas integrais, através de simulações determinísticas. Ambas as metodologias serão tratadas nas seções subsequentes.

85

5.3.1 Método de Monte Carlo

O método de Monte Carlo é o mais popular método não-intrusivo e pode ser utilizado para qualquer problema de propagação de incerteza (KEANE e NAIR, 2005). Dado uma distribuição conjunta de probabilidade das variáveis aleatórias envolvidas, o método de MC pode ser aplicado para o cálculo aproximado da estatística da resposta de interesse, incluindo sua distribuição, com um nível de erro arbitrário, desde que seja fornecido um número suficiente de amostras. Esta abordagem é utilizada também para validar outras técnicas, no cálculo das estatísticas. No método de MC as funções f (ξ ) de interesse, são calculadas em m pontos ξ (i ) , i = 1… m gerados aleatoriamente a partir de suas distribuições P (ξ ) , então as integrais das Equações (5.7) e (5.9) são aproximadas, respectivamente, por:

(5.16) f

f MC =

1 m ∑ f (ξ (i ) ) m i =1

(5.17) σ [ f (ξ ) ] = σ 2f 2

σˆ 2f =

(5.16)

2⎤ 1 ⎡m f (ξ (i ) ) 2 ) − m ( f MC ) ⎥ ( ∑ ⎢ m − 1 ⎣ i =1 ⎦

(5.17)

onde m é o número de pontos amostrados, f MC é a aproximação de MC da média de f (ξ ) e σˆ f é a aproximação de MC do desvio padrão de f (ξ) . Se f é integrável em ξ , então f MC → f a medida em que m → ∞ . O cálculo da variância pode ser usada para verificar a aproximação f MC através da Equação (5.9), e então pode-se estimar o erro dado por σ f m (KEANE e NAIR, 2005 ), o qual independe do número de variáveis aleatórias do problema. O maior problema do método de MC é sua baixa taxa de convergência, já que o seu erro estimado é de ordem O(1 m ) . Por exemplo, para melhorar em um décimo a aproximação, é necessário uma amostra 100 vezes maior. Porém o método é facilmente paralelizável, pois o cálculo de cada f (ξ (i ) ) pode ser feito de forma independente. No caso das avaliações de f (ξ (i ) ) serem caras computacionalmente, recomenda-se o emprego de modelos substitutos. Pode-se calcular a partir das Equações (5.16) e (5.17) os gradientes da média e do desvio padrão em relação a variáveis quaisquer (aleatórias ou determinísticas). Estes gradientes serão requeridos durante o processo de otimização, e podem ser calculados derivando as Equações (5.16) e (5.17) em relação às variáveis de projeto x:

(5.18)

∂f (x, ξ ) ∂x

∂f MC 1 m ∂f (x, ξ (i ) ) = ∑ ∂x m i =1 ∂x

(5.18)

86

(5.19)

∂ (σ 2f )

∂ (σˆ 2f )

∂x

∂x

derivando σˆ f =

=

∂f (x, ξ (i ) ) ⎞ ∂f MC ⎤ 1 ⎡m ⎛ ⎢ ∑ ⎜ 2 f ( x, ξ ( i ) ) ⎥ ⎟ − 2m ( f MC ) ∂x m − 1 ⎣ i =1 ⎝ ∂x ⎦ ⎠

(5.19)

(σˆ ) , tem-se que 2 f

2 ⎡m ⎛ ∂f (x, ξ (i ) ) ⎞ ∂f ⎤ 1 2 −1/2 ∂ (σˆ f ) 1 = (σˆ f ) = − m ( f MC ) MC ⎥ (5.20) ⎢ ∑ ⎜ f ( x, ξ ( i ) ) ⎟ 2 σˆ f ( m − 1) ⎣ i =1 ⎝ ∂x ∂x ∂x ∂x ⎦ ⎠

∂σˆ f

(5.20) Um ponto muito importante, quando se utiliza o método de MC para um processo de otimização iterativo, é usar uma mesma semente (“seed”), para a geração de amostras parentes, para as várias etapas do processo de otimização. Caso a PDF das variáveis aleatórias não dependam das variáveis de projeto, as amostras geradas devem ser sempre iguais. Caso contrário, as amostras devem sofrer alterações em função do valor da variável de projeto correspondente. Este requisito é fundamental, principalmente, para alguns métodos que iterativamente aproximam gradientes ( ∇f ) ou hessianas ( ∇ 2 f ) de uma função f calculada a partir de MC, pois a utilização de amostras independentes (aleatórias) para cada iteração gera um erro aleatório O(1 m ) dificultando (ou até impossibilitando) a convergência de tais métodos. Nesse caso são necessários outros métodos de otimização que lidam diretamente com as incertezas, é o caso de métodos de otimização estocástica como algoritmos evolucionários, métodos de aproximação estocástica (“stochastic approximation methods”), entre outros (ANDRADÓTTIR, 1998; SCHUËLLER e JENSEN, 2008). Esta dificuldade no uso de amostras independentes é ressaltada quando se considera o cálculo da função f no mesmo ponto em duas iterações sucessivas, o que leva a um gradiente (ou hessiana) indeterminado, pois a função f terá dois valores (ou gradientes) distintos para o mesmo ponto, causado pelo uso de amostras diferentes. (a) Exemplo Para ilustrar a conseqüência do uso de amostras aleatórias e a importância de se fix2 (1 + x − x 2 ) xar a semente (“seed”), considere a função f ( x, ξ ) = ξ 2 − ξ , onde 10 3 ξ ~ N ( 0,1) , e pretende-se encontrar x* que minimize a média da função f ( x, ξ ) em relação à ξ . A média de f ( x, ξ ) pode ser analiticamente encontrada através das Equações (5.7) e (5.13), obtém-se f ( x) = x 2 10 , logo x* = 0 . Foi realizada uma aproximação da média da função de interesse por MC ( f MC ( x) ) utilizando uma amostra ξ com 30 pontos, para vários valores de x, variando-o de -1 a 1 com incrementos de 0.01. Porém duas estratégias diferentes para a construção das amostras foram consideradas: 87

mantendo a mesma amostra ξ para todos os pontos de x (“seed” fixo) e gerando uma nova amostra aleatória para cada ponto de x (“seed” aleatório). Na Figura 5.3 (a) é apresentado o gráfico de f ( x) (média de f ( x, ξ ) ) em função de x, calculada por quatro procedimentos diferentes: 1. Cálculo analítico já mencionado, 2. Cálcular via PCM com 3 pontos de integração, 3. Aproximar f MC ( x) utilizando a mesma amostra para todos os pontos (“seed” fixo) 4. Aproxima f MC ( x) utilizando uma amostra diferente para cada ponto de x (“seed” aleatório). Nos itens 3 e 4, o MC foi feito com 30 pontos. Para a função utilizada neste exemplo, o PCM com 3 pontos encontra a resposta exata da média e do desvio padrão (este método será tratado adiante), por isso as curvas da solução analítica e da solução por PCM estão sobrepostas (exato* (PCM)).

a)

b)

Figura 5.3 Monte Carlo: a) Gráfico da média de f ( x, ξ ) e b) Gráfico do erro f MC ( x) − f ( x) .

Na Figura 5.3, nota-se claramente a dificuldade em otimizar a função f MC ( x) usando MC com amostras aleatórias (MC - seed aleatório) indicada em vermelho. Já o MC com o uso da mesma amostra (MC - seed fixo), mesmo apresentando um erro médio considerável, apresenta uma curva tão suave quanto a função f ( x, ξi ) para um ξi qualquer. Isto ocorre, pois se f MC ( x) trata da média das funções fi ( x) = f ( x, ξi ) , i = 1..m, onde m é o tamanho da amostra. Assim, a função f MC ( x) pode ser otimizada e se encon* = trar ótimos adequados. Neste caso o ótimo de f MC ( x) será xˆMC

ξ 0.6ξ 2 + 2ξ

para qual-

quer amostra (fixa) de ξ , pois

88

(5.21) f MC ( x) =

df MC x2 2 (1 + x − x 2 ) 1 x ξ −ξ e ( x, ξ ) = (3ξ 2 + 10ξ ) − ξ dx 10 3 15 3

(5.21)

onde ξ e ξ 2 são respectivamente a média da amostra e a média dos quadrados da amostra. À medida que a amostra aumenta (melhora) ξ 2 − (ξ

)

2

→ σ ξ2 (Equação (5.9)),

com isso ξ e ξ 2 tendem, respectivamente, a zero e a um, consequentemente, * xˆMC → x* . Para a amostra de 30 pontos utilizada no gráfico, ξ = -0.0275 e * xˆMC = -0.0469 .

5.3.2 Técnicas de amostragem A geração das amostras pode ser feita de maneira totalmente aleatória (ou pseudoaleatória), ou utilizando técnicas mais eficientes para o plano de amostragem (“Design of experiments” – DoE) (GIUNTA et al, 2002) tal como o método LHS (“Latin Hipercube Sampling”) que melhora a distribuição dos pontos da amostra e conseqüentemente aumenta a convergência do MC. As amostras com distribuição Lognormal são geradas a partir de amostras com distribuição Normal, estas são geradas a partir de amostras com distribuição Uniforme. A geração computacional desta ultima é feita através de algoritmos determinísticos, capazes de gerar recursivamente uma sequência finita de números inteiros ou de ponto flutuante, com um determinado período, sendo por isso chamados de números pseudo-aleatórios. No presente trabalho as metodologias utilizadas para a geração das amostras serão uma técnica pseudo-aleatória e o LHS. Na primeira abordagem será considerado o algoritmo “Mersenne Twister” (MT) (MATSUMOTO, 1998) para a geração de amostras com distribuição uniforme, seguido por um algoritmo polar para obter as amostras com distribuição Normal. Ambos os algoritmos usados são do ambiente MATLAB 7.5 (MATHWORKS, 2007). O LHS é um método utilizado para a geração de uma amostra que cubra mais eficientemente o espaço das variáveis aleatórias, para um determinado número de pontos. A sua idéia básica é dividir o intervalo de cada uma das n dimensões da amostra pelo número de pontos pretendidos N, onde cada subintervalo tem a mesma probabilidade de ocorrência. Os N pontos são dispostos de forma que cada subintervalo de cada uma das n variáveis tenha apenas um ponto, para mais detalhes vide (BARÓN, et al , 1999). As amostras LHS são geradas a partir de uma amostra aleatória com distribuição normal, a qual é ajustada para que as distribuições marginais de cada variável se aproxime da sua distribuição de probabilidade teórica (STEIN, 1987). Quando se faz o cálculo das estatísticas de uma amostra de variáveis aleatórias geradas por alguma das técnicas mencionadas, em geral esses parâmetros não apresentam os mesmos valores originais das variáveis aleatórias, calculados através das PDF das variáveis. Estes parâmetros são desejados, pois, como se pode notar na aproximação de MC do exemplo anterior, eles estão diretamente ligados à convergência da aproxima-

89

ção. Pode-se ver no exemplo anterior, que quanto mais esses valores se aproximam de seus valores originais, melhor será a aproximação do MC para esta amostra. Deste modo, como será utilizada a mesma amostra-base, durante todo o processo de otimização, foi criada uma estratégia para a seleção de uma amostra, escolhida a partir de um conjunto de amostras. A amostra selecionada será a que apresentar as estatísticas mais próximas da desejada, a partir de uma determinada norma. Para a distinção das amostras foi proposta uma função que faz uma ponderação do erro (de forma empírica) das estatísticas da amostra. As amostras são geradas considerando uma distribuição Normal padrão ∼ N ( 0,1) . Essa equação empírica, que qualifica as amostras, faz a ponderação dos logaritmos dos erros quadráticos dos momentos centrais estatísticos das amostras, dando maior peso aos momentos de maior ordem. Esta função é definida como:

(5.22) FErr = ln ( μ

2

)+

ln ( (σ x − 1) 2 ) 2

+

ln ( s 2 ) 6

+

ln (κ 2 ) 24

(5.22)

onde μ é a média da amostra, σ x o seu desvio padrão, s a obliqüidade e κ a curtose.

Para uma amostra com distribuição Normal padrão N ( 0,1) : μ , s e κ devem ser iguais

a zero, enquanto σ x deve se aproximar de um. A amostra selecionada é a que apresenta o menor valor de FErr . Este procedimento não pretende escolher a amostra que conseguirá a melhor aproximação por MC do parâmetro de interesse, apenas se quer evitar que ocasionalmente se utilize uma amostra deficiente que empobreça a aproximação. Vale salientar que o mesmo foi criado de maneira intuitiva e adaptado em função de experimentos numéricos realizados. Outras duas funções “teste”, para se diferenciar as amostras, foram aplicadas: o teste de Kolmogorov-Smirnov e o teste de Anderson-Darling, (NIST/SEMATECH, 2009). Elas testam o quanto a CDF de uma dada amostra se aproxima da CDF de uma dada distribuição. Os resultados de propagação de incerteza obtidos com amostras LHS selecionadas por estes critérios não apresentaram grande avanço em comparação com as amostras LHS originais (isto é, sem critério de seleção), para o exemplo considerado (especificado a seguir).

(a) Exemplo - MC por diferentes amostragens Para exemplificar o uso da amostra selecionada, considere a função f ( X , Y ) = sen( X )cos(Y ) , onde X e Y são variáveis aleatórias ∼ N (1, 0.1) . Foi realizada uma aproximação da média f MC e do desvio padrão σˆ f de f ( X , Y ) por MC, utilizando três técnicas de amostragem diferentes: totalmente aleatória, LHS padrão e LHS selecionada, para este último, foram analisados conjuntos de cinco mil amostras. Variou-se o tamanho da amostra de 64 até 16.384 para cada tipo de amostra. Este procedimento foi repetido 50 vezes, pois seus resultados são aleatórios, então foi calculado o erro médio no cálculo de f MC de cada tipo de amostra. 90

No caso das amostras LHS selecionadas, para cada uma das 50 repetições feitas, um conjunto diferente de cinco mil amostras foi analisado, consequentemente 50 amostras selecionadas diferentes foram obtidas (totalizando 250.000 amostras analisadas e 50 selecionadas), para cada tamanho de amostra diferente. Na Figura 5.4 é apresentado o gráfico da do erro médio de f MC , para os diferentes tamanhos de amostras (de 64 até 16384) pelos três tipos de amostragem. O erro foi medido em relação à média exata f calculada simbolicamente pelo do MATLAB, através das Equações (5.7), (5.8) e (5.13).

a)

b)

Figura 5.4 Erro médio na aproximação por MC utilizando três técnicas de amostragem diferentes: totalmente aleatória, LHS padrão e LHS selecionada: a) Erro da média e b) Erro do desvio padrão. Como se pode observar na Figura 5.4 (a) a amostra LHS Selecionada apresentou melhor resultado no cálculo da média para tamanhos de amostras menores, porém à medida que o número de pontos aumenta a diferença entre as amostras LHS diminui, sendo necessário um conjunto maior de amostras para extrair uma que se sobressaia significativamente. Os resultados da LHS Selecionada chega a apresentar um erro médio maior que o LHS simples para amostras maiores. Já para o cálculo do desvio padrão os resultados não foram muitos diferentes, quando se compara o LHS com o LHS Selecionado. No gráfico nota-se também, a convergência mais acentuada do LHS em relação à amostra totalmente aleatória*, para o cálculo da média e do desvio padrão. Para este caso simples o custo computacional do LHS Selecionado não justificaria o seu uso, pois o tempo de CPU pra se gerar uma amostra LHS é maior que a propia avaliação da função f ( X , Y ) . Porém, para problemas em que o custo computacional, de se gerar e avaliar uma amostra LHS, é irrelevante em comparação ao custo do cálculo da função de interesse, o processo de seleção amostral mostrou-se adequado para os casos aqui considerados (seção 5.4.2). Assim sendo, nos exemplo analíticos posteriores serão usadas amostras LHS padrão, enquanto que nos exemplos práticos o processo de seleção amostral será empregado.

91

5.3.3 Método da Colocação Probabilística (PCM) Métodos tradicionais como o MC, mesmo com técnicas de amostragem que melhoram sua eficiência, são inviáveis para serem aplicadas diretamente em modelos complexos de alta fidelidade. O Método da Colocação Probabilística (“Probabilistic Collocation Method” - PCM) (RAMAMURTHY, 2005) é uma ferramenta desenvolvida visando uma análise de incerteza eficiente em modelos complexos e computacionalmente custosos. A idéia básica do PCM é aproximar a resposta do modelo em função das variáveis aleatórias, por funções polinomiais, e estimar as integrais das Equações (5.7) e (5.9) por Quadratura de Gauss (STOER e BULIRSCH, 1991). Em particular, um PCM de grau n1, aproxima a resposta do modelo por uma função polinomial de grau 2n-1, através de n pontos calculados, obtendo de forma exata as estatísticas de funções polinomiais de grau menor ou igual a 2n-1. Este método é indicado principalmente para problemas onde as funções de interesse são suaves, pois as aproximações polinomiais podem apresentar dificuldades para funções oscilatórias ou com singularidades. Além disso, o número de variáveis aleatórias consideráveis deve ser pequeno, pois o número de pontos necessários, para a aproximação de um mesmo grau, aumenta exponencialmente com o número de variáveis aleatórias. Esta dificuldade pode ser superada utilizando técnicas de integração por grades esparsas (“sparse grids) (HEISS e WINSCHEL, 2008). O PCM se baseia nos conceitos de Quadratura de Gauss e de polinômios ortonormais. A idéia básica do método é aproximar a função de interesse por funções polinomiais e calcular as integrais (5.7) e (5.9) por quadratura de Gauss. Então, antes de explanar sobre os pormenores da estrutura do PCM, é necessário apresentar os assuntos supracitados.

(a) Polinômios ortogonais Polinômios são ditos ortogonais entre si com relação a um produto interno relacionado a um espaço, se este for nulo. Todo produto interno em um dado espaço F , de funções reais f, g e h pertencentes a F , deve satisfazer as seguintes condições (APOSTOL, 1967): f , g + h = f , g + f ,h

(5.23)

α f , g = α f , g = f , α g , onde α é um escalar f , g = g, f

(5.23)

f , f > 0, se f ≠ 0

Dado um espaço linear de funções reais F e considerando duas funções polinomiais f ( x), g ( x) ∈ F , o produto interno aqui tratado é definido como:

92

(5.24) f ( x ) , g ( x ) = ∫ f ( x ) g ( x ) P( x)dx

(5.24)

F

onde P(x) é uma função de ponderação não negativa definida no espaço F . Este produto interno forma a base da integração por Quadratura de Gauss e para o PCM. Como já mencionado, as funções polinomiais f ( x), g ( x) ∈ F são ortogonais se seu produto interno for nulo. Um conjunto de polinômios hi ( x) pertencentes ao espaço polinomial H , pode ser definido como: i

(5.25) hi ( x) = ai ,0 + ai ,1 x + ai ,2 x 2 ... + ai ,i x i = ∑ ai , j x j

(5.25)

j =0

Estes polinômios serão ortonormais em relação a uma função de ponderação P(x), se a seguinte relação existe para todos hi ( x) , i=0,1..n, onde o índice de hi indica o grau do polinômio: ⎧1, para i = j (5.26) hi ( x), h j ( x) = ⎨ ⎩0, para i ≠ j

(5.26)

Através destas relações são encontrados os coeficientes ai , j que definem os polinômios ortonormais. Estes polinômios são únicos para cada função de ponderação dada, e formam uma base para H . Todas as raízes x*j , j = 1..i de um polinômio hi ( x) , estão contidas no espaço real, ou seja, hi ( x*j ) = 0, x*j ∈ F , j = 1..i , e dependem apenas da função de ponderação P(x). As raízes de hi ( x) formam os pontos de colocação da quadratura de Gauss. Para mais detalhes consulte (GAUTSCHI, 2005). Note que h0 ( x ) é uma constante, h0 ( x ) = a0,0 , consequentemente propriedades importantes, a partir da Equação (5.26), podem ser obtidas, as quais serão usadas adiante.

(5.27)

h0 ( x ) , hi ( x ) = 0, h0 ( x ) , h0 ( x ) = 1,

∫ h ( x ) P( x)dx = 0, a ∫ P ( x)dx = 1 F

i

i>0

(5.27)

2

0,0

F

(b) Quadratura de Gauss Na integração numérica via quadratura de Gauss para integrais da seguinte forma:

93

(5.28)

∫ f ( x) P ( x ) dx

(5.28)

F

Aproxima-se a função f ( x) , por um polinômio de grau 2n-1 , a partir da base ortonormal do espaço H , em relação à função de ponderação P( x) , onde n é o número de pontos de integração, tal como segue

(5.29) f ( x)

⎛ n −1 ⎞ ⎛ n −1 ⎞ fˆ ( x) = ⎜ ∑ bi hi ( x) ⎟ + hn ( x) ⎜ ∑ ci hi ( x) ⎟ ⎝ i =0 ⎠ ⎝ i =0 ⎠

(5.29)

Na integral da Equação (5.28) aproximada por quadratura de Gauss, o segundo termo da aproximação se cancela (por ortogonalidade) e lembrando as propriedades (5.27) os termos para i = 1..n − 1 do primeiro somatório da aproximação também é cancelado. A integral desejada (5.28), aproximada por quadratura de Gauss, pode ser expressa então da seguinte forma:

(5.30)

∫ f ( x ) P ( x ) dx F

b0 h0 ∫ P ( x ) dx

(5.30)

F

Para encontrar os coeficientes bi e ci da aproximação (Equação (5.29)), seria necessário o cálculo da função f ( x) em 2n pontos. Porém, como a integral não depende dos coeficientes ci, pode-se calcular a função f ( x) nas n raízes x* de hn ( x) , cancelando assim o segundo termo da aproximação apresentada na Equação (5.29), pois hn ( xi* ) = 0, i = 1..n . Os coeficientes ai , j que definem os polinômios ortonormais são calculados através das relações (5.26), como já mencionado. Desta forma, os coeficientes bi são encontrados resolvendo o sistema: ⎡ f ( x1* ) ⎤ ⎡ h0 n −1 ⎛ ⎞ ⎢ ⎥ ⎢ (5.31) f ( xi* ) = ⎜ ∑ b j h j ( xi* ) ⎟ ; ou ⎢ ⎥=⎢ * ⎥ ⎝ j =0 ⎠ ⎢ ⎢ ⎣ f ( xn ) ⎦ ⎣ h0

hn −1 ( x1* ) ⎤ ⎡ b0 ⎤ ⎥⎢ ⎥ ⎥⎢ ⎥ hn −1 ( xn* ) ⎥⎦ ⎢⎣bn −1 ⎥⎦

(5.31)

Definindo-se a matriz h j ( xi* ) inversa, como ⎡ h0 ⎢ (5.32) Ap = ⎢ ⎢ h0 ⎣

hn −1 ( x1* ) ⎤ ⎥ ⎥ hn −1 ( xn* ) ⎥⎦

−1

(5.32)

os coeficientes b0 podem ser calculados como:

94

⎡ f ( x1* ) ⎤ ⎡ b0 ⎤ n ⎢ ⎥ (5.33) ⎢⎢ bi ⎥⎥ = Ap ⎢ f ( xi* ) ⎥ , logo b0 = ∑ Ap1,i f ( xi* ) i =1 ⎢ f ( xn* ) ⎥ ⎢⎣bn −1 ⎥⎦ ⎣ ⎦

(5.33)

Para o cálculo da integral, basta conhecer a primeira linha da matriz Ap e ponderar as respostas da função nos pontos de integração (raízes de hn ), então, o vetor dos pesos P é definido como Pi = Ap1,i . Definindo-se (5.34) C0 = h0 ∫ P( x)dx

(5.34)

F

o valor da integral (5.30) é aproximado por b0C0 . Os passos para a determinação dos parâmetros necessários a aplicação da metodologia são: 1. A partir de uma função de ponderação qualquer P( x) (não negativa), encontra-se os polinômios hi ( x) , i=1..n (cujos coeficientes ai , j são obtidos utilizando (5.26)). 2. Calcula-se as raízes xi* de hn ( x) . 3. Calcula-se a matriz Ap . 4. Computa-se o vetor de ponderação P, onde Pi = Ap1,i . 5. Calcula-se o valor da integral utilizada na Equação (5.30) e defini-se uma constante C0 Equação (5.34). Restando apenas avaliar a resposta f ( xi* ) em cada ponto de integração, para então calcular o valor desejado: (5.35)

∫

F

n

f ( x) P( x)dx

C0 ∑ Pi f ( xi* )

(5.35)

i =1

Vale ressaltar que não é necessário calcular estes parâmetros para as funções de distribuição conhecidas, pois estes valores são conhecidos e podem ser encontrados na literatura, bem como em bibliotecas computacionais da maioria dos programas.

(c) Aplicando Quadratura de Gauss à estatística - PCM A avaliação das estatísticas definidas nas Equações (5.7) e (5.9) considerando o PCM (“Probabilistic Collocation Method”) é uma aplicação direta da quadratura de Gauss considerando o espaço das variáveis aleatórias ξ e sua PDF como função de 95

ponderação. Portanto tem-se que

∫

F

P(ξ )dξ = 1 e pela segunda propriedade em (5.27)

h0 = 1 , consequentemente a constante C0 = 1 (definida pela Equação (5.34)). Com isso o valor da integral (5.30) é aproximado apenas por b0 . Os polinômios ortonormais são definidos para cada PDF. A média e o desvio padrão de uma resposta de interesse serão aproximados pelo PCM da seguinte forma: n

f PC = ∑ Pi f (ξ*(i ) ) i =1

(5.36)

σˆ

2 PC

n

= ∑ Pi f (ξ ) − f PC * 2 (i )

i =1

(5.36) 2

onde ξ*(i ) , i = 1..n são as raízes do polinômio ortonormal. Pode-se calcular a partir da Equação (5.36) os gradientes da média e do desvio padrão em relação a variáveis quaisquer, aleatórias ou determinísticas, da mesma forma como foi feito para o MC, derivando a Equação (5.36) em relação às variáveis de projeto x

(5.37)

(5.38)

∂f (x, ξ ) ∂x

m ∂f (x, ξ (i ) ) ∂f PC = ∑ Pi ∂x ∂x i =1

∂ (σ 2f )

2 ∂ (σˆ PC )

∂x

∂x

derivando σˆ PC =

m ∂f (x, ξ (i ) ) ⎞ ⎛ ∂f PC = ∑ ⎜ 2 Pi f (x, ξ (i ) ) ⎟ − 2 f PC ∂x ∂x i =1 ⎝ ⎠

(5.37)

(5.38)

(σˆ ) , tem-se que 2 PC

2 ∂f (x, ξ (i ) ) ⎞ ∂σˆ PC 1 2 −1/2 ∂ (σˆ PC ) ∂f ⎤ 1 ⎡m ⎛ = (σˆ PC ) = − f PC PC ⎥ (5.39) ⎢ ∑ ⎜ f ( x, ξ ( i ) ) ⎟ 2 ∂x ∂x ∂x ∂x ⎦ σˆ PC ⎣ i =1 ⎝ ⎠

(5.39)

Uma dificuldade da integração por Quadratura de Gauss e que também padece o PCM é a chamada maldição dimensional (“dimensional curse”), pois o número de pontos de integração (para um mesmo grau de aproximação) cresce exponencialmente com o número de dimensões do problema. Ou seja, no caso do PCM, o número de variáveis aleatórias a ser considerada não deve ser elevado. Por exemplo, para um problema de 10 variáveis aleatórias, caso seja utilizado o PCM com 3 pontos de colocação para todas as variáveis (aproximação de 5º grau), serão necessários 59.049 (310) cálculos da função de interesse.

96

Para problemas com muitas variáveis aleatórias pode-se utilizar técnicas de grades esparsas (“sparse grids”) baseadas nas regras de Smolyak (1963) para integração multivariável, diminuindo consideravelmente o número de pontos de colocação, para mais detalhes vide (HEISS e WINSCHEL, 2008). Será mostrado aqui, como encontrar cada valor ξ (*i ) , i = 1..n , raízes do polinômio hn (ξ ) e os coeficientes de ponderação Pi , porém estes valores apenas devem ser calculados uma única vez para cada distribuição e seus valores devem ser armazenados para serem usados em problemas diversos. No nosso caso, será tratada apenas a distribuição normal padrão ∼ N ( 0,1) , pois os mesmos valores podem ser aplicados a uma distribuição normal com parâmetros diferentes, bem como a uma distribuição Lognormal, através de uma transformação de variáveis. Os polinômios ortogonais para uma distribuição gaussiana são conhecidos como polinômios Hermitianos. Como exemplo, será calculado os coeficientes ai ,k , k = 0..i dos polinômios ortonormais hi (ξ ), i = 0..n , os quais foram definidos na Equação (5.25), usando como função de ponderação a PDF Normal padrão (Equação (5.13)). Isso será feito até ser encontrado o polinômio ortonormal do primeiro grau, ou seja, i = 0,1 , através das relações descritas na Equação (5.26) e sabendo que para uma distribuição normal padrão

(5.40)

∫

∞

−∞

ξ P(ξ )d ξ = 0 ,

∫

∞

−∞

ξ 2 P(ξ )dξ = 1 e ξ ∈ , obtém-se:

h0 (ξ ), h0 (ξ ) = 1

(5.40)

∞

a0,0 2 ∫ P(ξ )dξ = 1 −∞

logo, a0,0 = 1 . E para h0 (ξ ), h1 (ξ ) = 0

(5.41)

∫ ( a )( a + a ξ ) P(ξ )dξ = 0 a ∫ P(ξ )dξ + a ∫ ξ P(ξ )dξ = 0 ∞

0,0

−∞

1,0

(5.41)

1,1

∞

∞

1,0 −∞

1,1 −∞

logo, a1,0 = 0 . E ainda, para h1 (ξ ), h1 (ξ ) = 1 (5.42)

∫ (a a ∫ ∞

1,0

−∞

2

1,0

∞

−∞

+ a1,1ξ )( a1,0 + a1,1ξ ) P(ξ )d ξ = 1

(5.42)

∞

∞

−∞

−∞

P(ξ )d ξ + 2a1,1a1,0 ∫ ξ P (ξ )d ξ + a1,12 ∫ ξ 2 P (ξ )d ξ = 1

logo, a1,1 = 1 .

97

Portanto, os dois primeiros polinômios ortonormais são: h0 (ξ ) = 1 e h1 (ξ ) = ξ . Para se encontrar um polinômio de grau j, a partir dos polinômios de menor grau já calculados (de 0 à j-1), resolve-se um sistema de j+1 equações não lineares definidos na Equação (5.26), onde i = 0...j . Assim o seguinte sistema é resolvido:

(5.43)

∞ ⎛ i ⎞ ⎧0, i ≠ j ⎞⎛ j hi (ξ ), h j (ξ ) = ∫ ⎜ ∑ ai ,k ξ k ⎟ ⎜ ∑ a j ,k ξ k ⎟ P(ξ )d ξ = ⎨ −∞ ⎝ k =0 ⎠ ⎝ k =0 ⎠ ⎩1, i = j

(5.43)

para i = 0… j . Onde todos os coeficientes ai ,k dos polinômios hi (ξ ) , para i = 0… j − 1 , já são conhecido. O sistema é resolvido considerando os coeficientes a j ,k do polinômio h j (ξ ) , como incógnitas, onde k = 0… j . Esse procedimento é repetido recursivamente variando o j até atingir o polinômio de grau n ( hn (ξ ) ) desejado. Depois de encontrados os n polinômios ortonormais, continua-se analogamente os mesmos passos de 1 a 5 utilizados na quadratura de Gauss, Item 5.3.3 (b)., i.e: 1. Calcula-se as n raízes ξ j * , j = 1..n de hn (ξ ) , 2. Cria-se a matriz de hi (ξ j * ) avaliando os polinômios hi (ξ ), i = 0… n − 1 em cada ξ j * , j = 1..n . 3. Calcula-se sua inversa: Ap 4. Extraí-se da 1ª linha de Ap os pesos Pi = Ap1,i . Os valores de Pj e ξ j * , j = 1..n são armazenados e utilizados para obter as estatísticas de qualquer função em estudo, a partir da Equação (5.36). Para poder aplicar os valores de ξ* , obtidos para uma distribuição normal padrão ∼ N ( 0,1) , a uma distribuição normal com parâmetros diferentes ∼ N ( μ , σ ) , basta realizar um simples transformação (WEBSTER et al, 1996). (5.44) ξ ' = ξ *σ + μ

(5.44)

onde ξ ' será o vetor das raízes do polinômio ortonormal com relação à uma distribuição normal ∼ N ( μ , σ ) . Os pesos P são os mesmos para este caso, bem como para o caso de uma distribuição lognormal. Através de uma transformação de variáveis inversa à citada na Equação (5.14), se obtém as raízes ξ 'log do polinômio ortonormal para uma distribuição lognormal. Esta transformação é realizada da seguinte forma: (5.45) ξ 'log = exp (ξ *θ + η )

(5.45) 98

onde ⎛ μ2 (5.46) η = ln ⎜ ⎜ μ2 +σ 2 ⎝

⎞ ⎛σ 2 ⎞ ⎟ e θ = ln ⎜ 2 + 1⎟ ⎟ ⎝μ ⎠ ⎠

(5.46)

(d) Implementação Computacional Para o desenvolvimento da metodologia, foi implementado no MATLAB uma função para o cálculo dos coeficientes ai , j de cada polinômio ortonormal. Esta função simbolicamente calcula as integrais de

∫

b

a

ξ i P(ξ )dξ para i = 0...2n − 1 , onde a função de

ponderação P (ξ ) , os limites [a,b] onde ela é definida e o valor de n são dados de entrada. A partir dos valores das integrais obtidos, cria-se um sistema de equações não lineares, apresentado na Equação (5.43), o qual também é resolvido simbolicamente, para se obter os coeficientes ai , j de cada polinômio ortonormal. Como exemplo da aplicação da metodologia citada, os coeficientes que formam os dez primeiros polinômios ortonormais (até o polinômio de nono grau), obtidos para uma distribuição Normal padrão com limites [a,b]=[-∞, ∞], estão indicados abaixo: ai , j = [

1]

[

0,

1]

[

-1,

0,

[

0,

-3,

0,

[

3,

0,

-6,

0,

[

0,

15,

0,

-10,

[ -15,

0,

45,

0,-105,

0,

[

[ 105, [

1]/2*2^(1/2)

0,-420,

0, 945,

1]/6*6^(1/2) 1]/12*6^(1/2) 0,

0, -15, 105,

0,

0, -21,

0, 210,

0,-1260,

1]/60*30^(1/2) 1]/60*5^(1/2) 0,

0, -28,

0, 378,

0,

1]/420*35^(1/2) 0, -36,

1]/1680*70^(1/2) 0,

1]/5040*70^(1/2)

para i = 0… n, e j = 0… i . Estes polinômios ortonormais estão de acordo com os polinômios ortogonais encontrados na literatuda (WEBSTER et al, 1996), após normalização. Depois de encontrado os polinômios ortonormais, pode-se encontrar os pontos de colocação calculando suas raízes através da função “roots” do próprio MATLAB (MATHWORKS, 2007). Para o polinômio ortonormal de nona ordem ( h9 (ξ ) ), encontra-se:

ξ* = 0

99

-4.512745863399775 4.512745863399775 -3.205429002856465 3.205429002856465 -2.076847978677832 -1.023255663789131 2.076847978677832 1.023255663789131

Os resultados acima obtidos estão em concordância com os apresentados na referência (GREENWOOD e MILLER, 1948), onde são obtidas as raízes de polinômios ortogonais Hermitianos, sendo necessário a multiplicação pela raiz de 2, para se obter as raízes dos polinômios ortonormais em relação à função de peso Gaussiana utilizada. Estes nove pontos de colocação são utilizados para uma aproximação de 17ª ordem. Como já foi mencionado, a partir destes pontos, cria-se a matriz de hi (ξ j * ) avaliando os polinômios hi (ξ ), i = 0… n − 1 em cada ξ j * , j = 1..n . Só restando então calcular Ap, invertendo a matriz criada, e extrair da 1ª linha os pesos P, onde Pi = Ap1,i . Os pesos obtidos para os nove pontos de colocação foram:

P= 0.406349206349205 2.2345844007746e-5 2.2345844007746e-5 2.789141321231750e-3 2.789141321231750e-3 4.9916406765217772e-2 0.244097502894940 4.9916406765217772e-2 0.244097502894940

Os resultados acima obtidos estão em concordância com os apresentados na referência (GREENWOOD e MILLER, 1948), onde são obtidas os coeficientes de ponderação para os polinômios ortogonais Hermitianos, sendo necessário dividí-los por π , para se obter os coeficientes de ponderação dos polinômios ortonormais em relação à função de peso Gaussiana utilizada. Vale notar que a partir da matriz Ap pode-se encontrar os coeficientes da aproximação polinomial para posterior verificação, como será feito mais adiante. Os valores de P e ξ* são armazenados e utilizados para obter as estatísticas da função em estudo a partir da Equação (5.36).

100

5.3.4 Exemplos (a) Verificação - PCM Além da verificação, já efetuada no item anterior, relacionada com a obtenção dos parâmetros do método: pontos de colocação e coeficientes de ponderação, a verificação geral da metodologia irá ser feita aqui, através do cálculo da média e do desvio padrão para uma função polinomial de grau N P = 8 de uma variável aleatória ξ ∼ N (0,1) , para este caso o PCM deve obter uma solução exata para um determinado número de pontos. A função polinomial que será utilizada é:

f (ξ ) = (5.47)

1 NP 1 ξ − i (−1)i ) , ou f (ξ ) = (ξ + 1)(ξ − 2 )(ξ + 3) ... (ξ − 8 ) ( ∏ C i =1 C NP

(5.47)

onde, C = ( N P !) ∏ −(−1) = 1(−2)3...7(−8) i

i =1

Vale salientar que para o cálculo exato da média de um polinômio de 8º grau, são necessários 5 pontos de colocação, pois 2n-1 = 9. Já para o cálculo do desvio padrão, a função aproximada pelo PCM é a f (ξ ) 2 , o que leva a um polinômio de 16º grau, sendo necessários 9 pontos de colocação (2n-1 = 17) para a convergência do PCM. A Tabela 5.1 ilustra o erro no cálculo da média e no desvio padrão pelo PCM, para aproximações que utilizam de 3 até 9 pontos, onde nota-se a convergência esperada do método. Tabela 5.1 Aproximação da média e do desvio padrão via PCM de uma função polinomial de oitavo grau. PCM n

2n-1

Erro na média

Erro no desvio padrão

3 

5 

0.012053546253266     0.055207201802620 

4 

7 

0.000595258325395     0.029324964183213 

5 

9 

0.000000000000000     0.009362017482057 

6 

11 

0.000000000000001     0.002193953762619 

7 

13 

0 

8 

15 

0.000000000000000     0.000022214207163 

9 

17 

0.000000000000000     0 

   0.000319231598934 

101

Como comparativo, a Tabela 5.2 ilustra a aproximação do mesmo caso por MC utilizando amostras do LHS, onde se percebe a grande vantagem do PCM para este tipo de caso. Tabela 5.2 Aproximação da média e do desvio padrão via. MC Nº pontos

Erro da média

Erro do desvio padrão

103 

   0.001008014052612 

   0.000064040096951 

209 

   0.000299157888075 

   0.001628113576030 

327 

   0.000301517289170 

   0.000156040536578 

481 

   0.000531314411038 

   0.000984779588322 

743 

   0.000688695761748 

   0.001169150370338 

1329 

   0.000185486317499 

   0.000273969426433 

2887 

   0.000072895050160 

   0.000098277073118 

7361 

   0.000036488736102 

   0.000024487079049 

20583 

   0.000009992521449 

   0.000003760397742 

60049 

   0.000002035240651 

   0.000005664316955 

A Figura 5.5 sumariza os resultados obtidos pelos dois métodos o MC utilizando amostras LHS e o PCM, onde o erro mínimo foi fixado em 10-10. Percebe-se a grande vantagem do PCM para este caso.

Figura 5.5 Teste de convergência dos métodos MC e PCM, para um caso polinomial. (erro mínimo fixado em 10-10).

102

(b) Exemplo - Função Periódica Neste exemplo será comparado o PCM com o LHS-MC considerando-se novamente a função f ( X , Y ) = sen( X )cos (Y ) , onde X e Y ∼ N (1, σ ) são variáveis aleatórias independentes. Para este exemplo o valor do desvio padrão σ será variado para que se tenham diferentes “regiões de abrangência” (regiões onde há probabilidade de ocorrência significante) das variáveis aleatórias. O MC foi testado com amostra LHS de 50.000 pontos, enquanto que para o PCM utilizou-se 7 pontos de colocação para cada variável, totalizando 49 pontos de colocação analisados. Assim um estudo paramétrico sobre o desvio padrão σ das variáveis aleatórias foi conduzido, variando-o de 0.2 até 2.6. A Figura 5.6 mostra para cada valor de σ , no intervalo (de 0.2 até 2.6), o valor das estatísticas da função, média (Figura 5.6(a)) e desvio padrão (Figura 5.6(b)), bem como os erros de cada método no cálculo da média (Figura 5.6(c)) e no cálculo do desvio padrão (Figura 5.6(d)).

(a) média

(c) erro da média

(b) desvio padrão

(d) erro no desvio padrão

Figura 5.6 Estudo paramétrico da função f ( X , Y ) variando σ das variáveis aleatórias: (a) média, (b) desvio padrão, erros de cada método no cálculo (c) da média e (d) do desvio padrão. 103

Note que, quanto menor o desvio padrão das variáveis aleatórias, mais “simples” será sua região de influência na função f ( X , Y ) (menos oscilações, i.e. picos e vales). Como pode-se observar o PCM com 49 pontos se mostrou melhor para valores menores de σ (menor região de abrangência), para todos as aproximações da média e do desvio padrão. Já o MC apresentou um desempenho quase constante, não apresentando grande variação na ordem de grandeza do erro com o aumento do desvio padrão. Vale frisar que com os parâmetros utilizados neste exemplo, a solução via o PCM é três ordens de grandeza mais rápido que o LHS-MC, devido a diferença no número de pontos em que a função é avaliada por cada método (49 pontos o PCM e 50x10³ pontos o MC). Na Figura 5.7 é apresentada a superfície da função f ( X , Y ) , os pontos de colocação do PCM e os pontos de integração do MC, para o valor inicial de σ = 0.2 , o que da uma idéia da região de abrangência inicial. Os resultados deste caso correspondem aos indicados pelos inícios das curvas apresentadas na Figura 5.6.

Figura 5.7 Superfície da função f ( X , Y ) , pontos de colocação do PCM, e amostra do MC para o valor inicial de σ = 0.2 . Na Figura 5.8 é apresentada a superfície da função f ( X , Y ) e os pontos de colocação do PCM, para o valor de σ = 1.6 , o que já seria um valor “elevado” de desvio padrão, porém como pode-se ver nas Figura 5.6 (c) e (d), através dos pontos circulados no gráfico, o PCM com 49 pontos ainda apresenta um resultado melhor que o MC com 50.000 pontos, tanto para a média quanto para o desvio padrão, este ultimo com pequena diferença. Só a partir deste ponto, com o aumento do desvio padrão da variável aleatória ( σ > 1.6 ), o MC com amostra LHS de 50.000 pontos começa a apresentar resultados melhores que o PCM no cálculo de desvio padrão da função. Para o cálculo da média o PCM obteve erros menores para todo o intervalo do σ analisado. 104

(a)

(b) Figura 5.8 Valores de σ = 1.6 (ou183°): (a) Pontos de colocação do PCM e (b) Pontos de integração do MC. À medida que se aumenta o valor de σ , aparecem regiões de mínimos e máximos e um comportamento multimodal é verificado para a função f ( X , Y ) . Quando o número de regiões de mínimos e máximos torna-se elevado para um mesmo número de pontos de colocação, a aproximação do PCM apresenta dificuldades. Uma alternativa seria usar métodos que utilizam aproximações discretas, dividindo o domínio das variáveis aleatórias, como o “Multi-Element Probabilistic Collocation Method” (ME-PCM) (FOO, WAN e KARNIADAKIS, 2008). 105

(c) Exemplo – Função com singularidade Neste exemplo serão testados o PCM e o MC para o caso de uma função não continuamente diferenciável, confrontando-os com a resposta exata. Considere a função f ( x, ξ ) = x − ξ − ξ (1 + x − x 2 ) , onde ξ ∼ N ( 0,1) é a variável aleatória. A Figura 5.9 ilustra o gráfico de f ( x) (a média de f ( x, ξ ) em função de x), com x variando de -1 à 1. Três procedimentos para o seu cálculo foram usados: •

O primeiro (MC - 30) aproxima f MC ( x) por MC utilizando uma amostra aleatória fixa de 30 pontos para ξ (para todos os valores de x);

•

O segundo faz o cálculo exato, através da Equação (5.7) (integrando simbolicamente pelo MATLAB);

•

Do terceiro (PCM - 3) ao sexto (PCM - 8) foi aproximada a média da função por PCM utilizando 3, 5, 6 e oito pontos de colocação respectivamente.

Como a função f ( x, ξ ) não é continuamente diferenciável (singular em x = ξ ), a aproximação polinomial “global” apresenta dificuldades. Como se pode notar na Figura 5.9 o PCM não consegue capturar a função de modo satisfatório. Já a aproximação por MC, mesmo para esta função singular, apresentou um comportamento semelhante ao obtido no item 5.3.1(a).

Figura 5.9 Gráfico de f ( x) (média da função f ( x, ξ ) em relação à ξ ), calculadas analiticamente, por MC e por PCM para diferentes graus de aproximação. Como para cada valor ξi da amostra, f ( x, ξi ) apresenta uma singularidade em xi = ξi , (i = 1..n) , então a função f MC ( x) terá n pontos de singularidade (onde n é o tamanho da amostras). Isto resulta em singularidades do gráfico f MC ( x) , que são mais perceptíveis ao se observar o gráfico do erro, mostrado na Figura 5.10. A Figura 5.10 mostra o erro da aproximação por PCM com 8 pontos de colocação e o erro de f MC ( x) 106

calculado por MC com 30 pontos. Para este caso o PCM apresentou erros maiores que os obtidos usando MC.

Figura 5.10 Gráfico do erro no cálculo de f ( x) (média da função f ( x, ξ ) em relação à ξ ), calculadas por MC e por PCM. Para melhor ilustrar as dificuldades enfrentadas pelo PCM, considere um caso particular da função do exemplo anterior, onde x=0: f (0, ξ ) = f (ξ ) = ξ − ξ onde

ξ ∼ N (1, 0 ) . Utilizando o PCM, aproximou-se a função f (ξ ) através da aproximação apresentada na Equação (5.29) e verificou-se o seu comportamento com a variação no grau dos polinômios utilizado. A aproximação completa apresentada na Equação (5.29), uma vez determinado os ξ tal que hn (ξ * ) = 0 , não necessita ser calculada, pois para qualquer valor de c a integral desejada (Equação (5.30) não tem seu valor alterado. Pode-se então definir qualquer valor aos coeficientes c. Neste estudo em particular foi analisada a forma de fˆ (ξ ) para ci=0 e apenas os coeficientes b, foram encontrados. *

Para diferentes graus aproximação, ou seja, para diferentes números de pontos de colocação n, as várias aproximações obtidas estão apresentadas na Figura 5.11(a). Na Figura 5.11(b) ilustra-se o gráfico de fˆ (ξ ) P (ξ ) , usada para a obtenção da média (Equação (5.7)), para os diferentes graus de aproximação. Em ambas as figuras a função analítica (exata) é também apresentada. Observa-se a oscilação da aproximação polinomial do PCM com o aumento do número de pontos, mostrando a dificuldade na convergência do método para casos onde a função a ser aproximada não é continuamente diferenciavel.

107

ξ a)

ξ b) Figura 5.11 PCM para diferentes graus de aproximação, aplicadas para a: a) função f (ξ ) e b) função f (ξ ) P(ξ ) . A Figura 5.12 ilustra o erro no cálculo da média via PCM para vários graus de aproximação, onde observa-se a convergência não monotônica do PCM para este caso. Onde o valor da média f (ξ ) , foi calculada analiticamente a partir da função f (ξ ) e da PDF P (ξ ) e vale 0,79788456 (Equação (5.7)).

108

Figura 5.12 Variação do erro no cálculo da média via PCM com o número de pontos. Como já mencionado, uma alternativa seria utilizar o “multi-element probabilistic collocation method” (ME-PCM) (FOO, WAN e KARNIADAKIS, 2008). Este método é uma generalização do PCM onde a aproximação é feita por subdomínios (elementos), a partir da discretização do domínio das variáveis aleatórias.

5.4 OTIMIZAÇÃO ROBUSTA Nos capítulos 3 e 4 foram apresentadas a formulação geral de um problemas de otimização uni e multiobjetivo, respectivamente. As soluções deste problema podem ser muito sensíveis à variação de parâmetros incertos, introduzidos no modelo como determinísticos. Sob estas circunstâncias os valores da(s) função(ões) objetivo e das restrições podem possuir uma grande probabilidade de sofrer variações inesperadas prejudiciais ao projeto. A Otimização Robusta leva em consideração as variáveis incertas (aleatórias) e suas probabilidades de ocorrência (PDF) de modo a encontrar um ótimo menos vulnerável a variabilidade dos parâmetros incertos.

5.4.1 Medidas de Robustez Considere um processo de otimização com incertezas, onde a função de interesse é f (x, ξ ) , x ∈ n são as variáveis de projeto e ξ ∈ E as variáveis incertas, onde E é o espaço amostral das variáveis. Um objetivo robusto e utópico seria encontrar um x* , de tal forma que para todo ξ ∈ E (5.48) f (x* , ξ ) ≤ f (x, ξ ) ∀x ∈

n

(5.48)

109

Achar tal x* que simultaneamente minimize f (x, ξ ) , para todos os ξ ∈ E , é o problema central da teoria da decisão estatística (KEANE e NAIR, 2005 Apud PRATT et al, 1996; BARROS, 2009). Porém, esta condição não pode ser atingida na maioria dos problemas, sendo necessário definir alguma medida de robustez que viabilize o problema, varias abordagens são usadas neste contexto. A primeira visa encontrar um x* que minimiza o pior caso ( ξ ∈ E mais desfavorável para f (x, ξ ) ), ou seja,

(5.49)

F (x) = max ( f (x, ξ ) ) ξ∈E

minn [ F (x) ]

(5.49)

x∈

Para problemas cujas variáveis aleatórias são continuas e ilimitadas e ainda não se pode determinar o máximo da função avaliada em todos os ξ ∈ E , outras abordagens podem ser consideradas. Visando diminuir o conservadorismo, pode-se considerar uma região menor A(ξ, ε ) , ao invés de todo o domínio das variáveis aleatórias. Ou seja, a otimização é realizada sobre o pior caso dentro da região A(ξ, ε ) , neste caso (5.50) F (x, ε ) = max

ξ∈ A ( ξ ,ε )

( f (x, ξ) )

(5.50)

onde ε é o parâmetro regularizador que define o tamanho da região, tal que: (5.51) lim ( F (x, ε ) ) → f (x, ξ 0 ) ε →0

(5.51)

onde ξ 0 é um determinado valor central da região A(ξ, ε ) e o parâmetro ε regula o conservadorismo da otimização. Esta abordagem é conhecida por regulação robusta (“robust regularization approach”) (BEYER et al, 2007). Outra abordagem semelhante é ver a função f como uma variável aleatória, onde em muitos casos, sua distribuição e seus limites são desconhecidos. Imaginando a função f como tendo, por exemplo, uma distribuição normal, pode-se admitir o pior caso com uma probabilidade p ( w) ( F (x, w) ), como sendo o valor de f onde há uma probabilidade p ( w) de f apresentar um valor menor, assim (5.52) F (x, w) = f (x) + wσ f (x)

(5.52)

onde f e σ f são a média e o desvio padrão, respectivamente, da função de interesse em relação a variável aleatória ξ ∈ E . O parâmetro w controla o conservadorismo da solução, por exemplo, para w=1, 2 e 3, tem-se o pior caso com probabilidades 84.13%, 97.72% e 99.87% respectivamente. Já para w=0 ( p ( w) = 50% ) apenas a média seria o110

timizada, não sendo considerada esta, uma abordagem de pior caso, pois se teria o pior caso com uma probabilidade de 50%. Quando w → ∞ tende-se a um pior caso com uma probabilidade de 100% e apenas o desvio padrão passa a ser considerado, assim (5.53) minn ⎡ lim F (x, w) ⎤ → minn ⎡⎣σ f (x) ⎤⎦ ⎣ w→∞ ⎦ x∈ x∈

(5.53)

Neste caso pode ser considerada a otimização do desvio padrão da função. Assim a solução será o ponto onde há a menor variação da resposta de interesse, esta constitui a idéia central da otimização robusta. A segunda abordagem, motivada pelo princípio de Bayes, tenta encontrar x* , tal que (5.54) f (x* , ξ ) ≤ f (x, ξ ) ∀x ∈

n

(5.54)

onde f é a média da função de interesse em relação a variável aleatória ξ ∈ E . Nota-se que a primeira abordagem é mais conservativa, já que se trata da otimização do pior caso. Por outro lado o principio de Bayes se preocupa com caso geral, já que apenas o valor da média é minimizado (KEANE e NAIR, 2005). Neste trabalho serão utilizados dois controles principais do objetivo: a média e o desvio padrão de uma dada função. Esta consideração tem sido adotada em várias referências no contexto da otimização robusta (DOLTSINIS e KANG, 2004; BEYER et al, 2007; SCHUËLLER e JENSEN, 2008). Otimizando a média encontra-se um ótimo menos conservador, pois pode existir uma probabilidade razoável do desempenho ser bem pior (ou melhor) que o valor encontrado. Enquanto que quando um processo de otimização é realizado sobre o desvio padrão se está sendo conservativo, e encontra-se o ponto onde há a menor variação da função de interesse, sendo esta uma das principais medidas de robustez. Como já mencionado este tipo de consideração leva a um problema de otimização multiobjetivo, para se encontrar soluções intermediárias de interesse (otimização multiobjetivo robusta - OMR). O problema de OMR mencionado anteriormente pode ser formulado como: Minimizar: ⎡⎣ E ( f (x) ) , σ ( f (x) ) ⎤⎦ gi ( x ) ≤ 0

i = 1,...m

hj ( x) = 0

j = 1,...

xk ≤ xk ≤ xku

k = 1,...ndv

onde E(.) é a esperança, definida na Equação (5.7), σ (.) o desvio padrão definido na Equação (5.8), f é a função de interesse e x o vetor das variáveis de projeto.

111

Para a solução deste problema, as técnicas de MO apresentadas no Capítulo 5 serão utilizadas.

5.4.2 Exemplo: Placa quadrada com orifício central O primeiro exemplo prático a ser considerado o problema apresentado na seção 2.5.1 e que se refere à análise do estado plano de tensões em uma placa quadrada com um orifício central. O Módulo de elasticidade da região 3 é uma variável aleatória com distribuição lognormal de média 5x104 MPa e desvio padrão 104 MPa. As dimensões do orifício são as variáveis do projeto selecionadas para a otimização. Os valores iniciais das variáveis de projeto são μ1 = μ2 = 50mm, e os limites inferior e superior impostos são 25mm e 75mm, respectivamente. Tal como definido, dois objetivos estocásticos são considerados são eles: minimizar a média e o desvio padrão da flexibilidade. Ao volume, é imposto ser inferior ou igual ao seu valor inicial. Além desta restrição, a média da tensão mais três vezes seu desvio padrão é restrito a magnitude de 7.0 N/mm². As soluções OMR serão obtidas considerando 15 pontos de Pareto. O problema de otimização pode ser formulado como: Minimizar: ⎡⎣ E ( f ( μ , E 3 ) ) , σ ( f ( μ , E 3 ) ) ⎤⎦ sujeito à E (τ eq ( i ) ( μ , E 3 ) ) + 3σ (τ eq (i ) ( μ , E 3 ) ) ≤ 7 MPa

i = 1,...nel

V ( μ ) ≤ V0 25mm ≤ μk ≤ 75mm k = 1,...ndv Onde f ( μ , E 3 ) é a flexibilidade, τ eq (i ) ( μ , E 3 ) é a tensão Von Misses no elemento i, V ( μ ) o volume da estrutura, V0 o volume inicial, nel é o número de elementos e ndv o

número de variáveis de projeto. O modelo de elementos finitos considerado possui 3900 graus de liberdade com elementos de dimensão média de 2mm, na configuração de referência onde μ1 = μ2 = 50mm. Para a aproximação via MBR, 3 regiões (indicadas na figura) são definidas. A base reduzida será construída no espaço viável das variáveis de projeto, bem como das variáveis aleatórias, D ={ [1, 9]x104,[25, 75]²} e o número de amostras analisadas foi N = 16.

(a) Determinação das amostras Para o cálculo das estatísticas do problema, serão usados os dois métodos já mencionados o MC e o PCM. Para definir o número de pontos (tamanho da amostras) utili112

zados por cada método, foi realizado um teste de convergência apresentado na Figura 5.13, Tabela 5.3 e Tabela 5.4. Os resultados do MC apresentam uma variabilidade, devido à aleatoriedade das amostras. Para quantificar esta flutuação nos resultados do MC todo o processo de cálculo das estatísticas da flexibilidade, foi repetido 100 vezes, para cada tamanho de amostra. E a partir destes dados foi calculado o desvio padrão dos resultados (“D. P.”), o que permite quantificar de modo aproximado a variação dos parâmetros calculados pelo MC para cada tamanho de amostra. A Figura 5.13 apresenta os resultados do MC para a média (Figura 5.13(a)) e o desvio padrão (Figura 5.13 (b)), da flexibilidade. O tamanho das amostras (“Número de pontos”) criadas por LHS foram variadas exponencialmente de 127 à 16255. As linhas pretas indicam o intervalo de flutuação dos resultados do MC, ou seja os resultados do MC variaram entre a linha preta superior e a linha preta inferior com uma certa freqüência. Considerando que um parâmetro estatístico calculado via MC segue uma distribuição normal, o intervalo entre as linhas pretas ( μ ± σ ) compreende aproximadamente 58% das ocorrências da aproximação. Vale salientar que este estudo só foi possível graças ao desempenho computacional conseguido através do MBR. Também foram analisados, para cada tamanho da amostra (“Número de pontos”), os resultados das estatísticas com uma amostra selecionada, indicado pela curva azul, amostra esta, utilizada no processo de otimização com MC utilizando uma amostra de 5000 pontos.

a) Média de f

b) Desvio padrão de f

Figura 5.13 Placa quadrada com um orifício central – Convergência do MC para o cálculo das estatísticas da flexibilidade: a) média e b) desvio padrão. A Tabela 5.3 mostra os resultados do MC com amostra LHS selecionada, ilustrado no gráfico da Figura 5.13 com a linha azul, para alguns tamanhos de amostras.

113

Tabela 5.3 Placa quadrada com um orifício central – Cálculo via MC com amostra LHS selecionada, para a média e desvio padrão. Tamanho

Média

da amostra

Desvio Padrão

202

0.308530959713673

0.035326555620201

403

0.308550213518812

0.035410907482470

806

0.308556674596121

0.035455905892811

1613

0.308558482597648

0.035455246318295

3225

0.308559106575475

0.035453460837231

6451

0.308559457041847

0.035451940103915

12902

0.308559963766802

0.035455220139521

A Tabela 5.4 mostra os resultados da média e do desvio padrão calculados pelo PCM, para diferentes números de pontos de colocação. Observa-se a convergência acentuada para poucos pontos, devido à suavidade da função flexibilidade em função do módulo de elasticidade. Tabela 5.4 Placa quadrada com um orifício central – Resultados do PCM Tamanho

Média

da amostra

Desvio Padrão

2

0.308541286854411

0.034854261923426

3

0.308559984101546

0.035445314900078

4

0.308560062653515

0.035455152887810

5

0.308560062770421

0.035455264870394

6

0.308560062772797

0.035455265822598

7

0.308560062777446

0.035455265989171

8

0.308560062777447

0.035455265989205

(b) Resultado da otimização A Figura 5.14 apresenta as distribuições dos pontos Paretos obtidos para os métodos aqui considerados para a otimização multiobjetiva, comparando os métodos de Monte Carlo (MC) e o Método da Colocação Probabilística (PCM). Tendo em vista os resultados obtidos na seção anterior, para o cálculo das estatísticas foram utilizados 5000 pontos para o MC e 5 pontos de colocação para o PCM. As curvas de Pareto via MC e PC estão de acordo, mesmo com a grande diferença no número de pontos calcula-

114

dos. Como pode ser observado, soluções pelos métodos NBI e NNC são as que conseguem obter pontos uniformemente espaçados em todas as partes da fronteira de Pareto.

a)WS

b)Min-Max

c)NBI

d)NNC

Figura 5.14 Placa quadrada com um orifício central – Pontos de Pareto: a) WS, b) MinMax, c) NBI, d) NNC. A Tabela 5.5 sumariza os desempenhos de cada método investigado em segundos. Soluções via PC são três ordens de magnitude mais rápidos em comparação com os resultados obtidos via MC. Tabela 5.5 Placa quadrada com um orifício central – desempenhos dos algoritmos. (segundos) MC 5000 pontos PC 5 pontos

WS

Min-Max

NBI

NNC

38 623

26 486

19 480

18 633

60

43

31

28

REFERÊNCIAS 115

APOSTOL T. M. “Calculus, Volume 1, One-Variable Calculus with an Introduction to Linear Algebra”. 2nd Edition, John Wiley & Sons, 1967. PAPOULIS, A. “Probability, Random Variables, and Stochastic Processes”. McGraw– Hill Kogakusha. 1965. MEYER, P. L. Probabilidade: aplicações à estatística, 2a edição, Rio de Janeiro, LTC. 1983. MARCZYK, J. “Stochastic multidisciplinary improvement: beyond optimization, American Institute of Aeronautics and Astronautics”, AIAA-2000-4929, 2000. BEYER, H. G.; SENDHOFF, B., “Robust optimization – A comprehensive survey”. Computational Methods and Applications in Mechanical Engineering. 196 (2007). SCHUËLLER, G.I.; JENSEN, H.A. Computational methods in optimization considering uncertainties – An overview. Computational Methods and Applications in Mechanical Engineering, 2008. MATSUMOTO, M.; NISHIMURA, T. “Mersenne Twister: A 623-Dimensionally Equidistributed Uniform Pseudorandom Number Generator”. ACM Transactions on Modeling and Computer Simulation, (1998), 8(1):3-30. ANDRADÓTTIR, S. “A review of simulation optimization techniques”. in: D. Medeiros, E. Watson, J. Carson, M. Manivannan (Eds.), Proceedings of the 30th Conference on Winter simulation, IEEE, Piscataway, NJ, pp. 151– 158, 1998. DOLTSINIS, I.; KANG, Z. Robust design of structures using optimization methods. Computational Methods and Applications in Mechanical Engineering, 194. 2004 KEANE, A. J.; NAIR P. B. Computational Approaches for Aerospace Design: The Pursuit of Excellence. John-Wiley and Sons. 602 p., August 2005. NIST/SEMATECH, e-Handbook of Statistical Methods. Disponível em: http://www.itl.nist.gov/div898/handbook, acessado em: set. 2009. BARÓN, J.H.; MAC NÚÑEZ LEOD, J.E. “SCALABILITY ON LHS SAMPLES FOR USE IN UNCERTAINTY ANALYSIS OF LARGE NUMERICAL MODELS”. International Conference on Safety and Reliability – ESREL ´99, Munich, Alemanha, 1999.

116

GAUTSCHI, W. “Orthogonal polynomials (in Matlab)”. Journal of Computational and Applied Mathematics, Vol. 178, 2005, pp. 215-234, 2005. MATHWORKS. “MATLAB User’s Guide”. Mathworks Inc., Natacki, 2007. STOER, J.; BULIRSCH, R. “Introduction to Numerical Analysis - Second Edition”. Springer-Verlag, Heidelberg, Berlin. p. 150-166, 1991. WEBSTER, M.; TATANG, M. A.; MCRAE, G. J. “Application of the Probabilistic Collocation Method for an Uncertainty Analysis of a Simple Ocean Model”. MIT Joint Program on the Science and Policy of Global Change, 1996. HEISS, F.; WINSCHEL, V., “Likelihood approximation by numerical integration on sparse grids”, Jornal of Econometrics, 144, p. 62-80. 2008. GREENWOOD, R. E.; MILLER, J. J. “Zeros of the Hermite Polynomials and Weights for Gauss' Mechanical Quadrature Formula”. THE UNIVERSITY OF TEXAS. Bull. Amer. Math. Soc. 54, MathSciNet, 1948. LAGAROS, N., PLEVRIS, V., PAPADRAKAKIS, M. “Multi-objective design optimization using cascade evolutionary computation”, Comput. Methods Appl. Mech. Engrg. 194 (2005) 3496–3515. FOO, J., WAN, X., KARNIADAKIS, G. E. “The multi-element probabilistic collocation method (ME-PCM): Error analysis and applications”. Journal of Computational Physics 227. 2008. BARROS, M., 2009. “Teoria da Decisao: NEYMAN-PEARSON e BAYES”, disponível em, http://www.mbarros.com/sitebuildercontent/sitebuilderfiles/Teoria_Decisao.pdf. Acesso em: 09 de set. de 2009. STEIN, M. “Large Sample Properties of Simulations Using Latin. Hypercube Sampling”. Technometrics, vol 29, no.2, May 1987. WAN, X.; KARNIADAKIS, G. E. “Multi-Element generalized Polynomial Chaos for arbitrary probability measures”, SIAM J. Sci. Comput. 28 (3) (2006) 901–928. RAMAMURTHY, D. “Smart simulation techniques for the evaluation of parametric uncertainties in black box systems”. Thesis (M.S. in computer science). Washington State University. 2005.

117

CAPÍTULO 6 CONCLUSÃO E TRABALHOS FUTUROS 6 CONCLUSÃO E TRABALHOS FUTUROS 6.1 PRINCIPAIS DESENVOLVIMENTOS Nesta dissertação procedimentos para a solução das diversas configurações de otimização apresentadas, foram desenvolvidos, implementados e verificados. Foram considerados incertezas e/ou múltiplos objetivos, permitindo a utilização de dois métodos para a avaliação dos parâmetro estatísticos e diversos esquemas para a obtenção da frente de Pareto, bem como diferentes formas para avaliação das funções relacionadas com a análise estrutural e/ou térmica e análise de sensibilidade. Todas as opções desenvolvidas foram incorporadas para análise de estruturas, considerando modelos elásticos lineares e térmicos, bidimensionais. Na etapa de análise, foi incluído a opção do uso do elemento quadrilateral bilinear (Q4) para o cálculo via MEF (o elemento triangular linear (CST) já havia sido implementado) e um gerador de malhas estruturadas uniformes. Além da implementação do MBR para o elemento quadrilateral, todo o procedimento de atualização para as diversas formas de mapeamento consideradas foi automatizado. Foram implementados os métodos para a obtenção de pontos de Pareto para problemas OM para a solução de problemas 2D contínuos. Os esquemas considerados foram: WS, Min-Max, NBI e o NNC. Também foi proposto uma modificação da estratégia NBI para problemas com mais de duas funções objetivo. Foram introduzidas as incertezas no problema de otimização, através dos métodos MC e o PCM, implementados e investigados para o cálculo das estatísticas de interesse, e assim obter resultados ótimos robustos sob múltiplos critérios. E por fim, foram apresentados os exemplos utilizando este conceito, para problemas de otimização estrutural multiobjetivo robusto. Os aspectos desenvolvidos que destacamos são: • Foram adaptadas e implementadas rotinas para análise estrutural e/ou térmica e análise de sensibilidade via Método dos Elementos Finitos, para a consideração do elemento Q4, além do CST já implementado; • Foram automatizados, desenvolvidos e implementados os procedimentos para análise estrutural e/ou térmica e o cálculo das sensibilidades utilizando o Método da Base Reduzida; • Foi realizado um estudo comparativo dos métodos através de exemplos envolvendo análise estrutural e/ou térmica. 118   

• Foram implementados os métodos: WS, Min-Max, NBI e NNC, além de uma modificação no NBI para problemas com mais de duas funções objetivo, a fim de obter os pontos de Pareto para problemas de otimização multiobjetivo, e seus resultados mostraram concordância com os da literatura; • Foram acoplados aos algoritmos de otimização multiobjetivo, as rotinas de otimização estrutural através de ambos os modelos (real e aproximado); • Uma rotina para o cálculo dos parâmetros estatísticos da resposta de interesse, além de seus gradientes, através do método de MC, foi implementada, com diferentes técnicas de amostragem (totalmente aleatória e LHS); • Foi implementado um procedimento para o desenvolvimento e aplicação do PCM para variáveis aleatórias de distribuições de probabilidade específica; • Foram feitos estudos em funções analíticas, com o objetivo de avaliar o comportamento do método MC e PCM para algumas situações específicas. • O MBR foi adaptado para considerar as variáveis aleatórias no espaço da base reduzida; • Além do MC, o PCM foi incorporado na rotina para o cálculo dos parâmetros estatísticos e seus gradientes, para o processo de otimização estrutural através de ambos os modelos (real e aproximado); • Soluções robustas e/ou multiobjetivo foram obtidas para problemas contínuos 2D. Diferentes esquemas para a obtenção dos parâmetros estatísticos de interesse, bem como, para a obtenção da frente de Pareto, foram incorporados em um algoritmo de SSO, de modo eficiente.

6.2 CONCLUSÕES DOS RESULTADOS OBTIDOS As principais conclusões extraídas, através dos resultados dos exemplos analisados, foram: • Os resultados da análise comparativa entre os elementos Q4 e o CST não indicam claramente o mais eficaz. Os testes de convergência indicam que a implementação do elemento quadrilateral foi adequada. • Do estudo comparativo realizado, entre o procedimento de análise via MBR e via MEF, foi mostrada uma convergência do MBR para poucos pontos amostrados, em relação ao número de variáveis de projeto. Devido à própria metodologia do MBR, o método não é indicado para problemas onde o número de variáveis de projeto se aproxima do número de graus de liberdade do modelo à ser aproximado. •

As técnicas NBI and NNC demonstraram ser as mais eficientes para problemas com duas funções objetivo.

•

Foram confrontados os resultados de problemas com mais de duas funções objetivo, obtidos através das quatro técnicas, bem como pela modificação proposta do NBI, que mostrou ser eficaz para os exemplos estudados. 119 

 

•

A importância de se incorporar técnicas de aproximação no procedimento de OM e OMR foi destacada.

•

Os resultados foram comparados com aqueles obtidos quando do uso de soluções que se baseiam integralmente no MEF. Conforme esperado, as soluções via o MBR foram acompanhadas por uma grande redução de tempo computacional, sem comprometimento substancial dos resultados das análises dos modelos.

•

Duas metodologias para o cálculo das estatísticas de várias respostas foram empregadas de modo satisfatório: o método de Monte Carlo (MC) e o método da colocação probabilística (PCM). Para a maioria dos casos, o número de ponto de integração necessários para uma aproximação com erros de mesma magnitude, foi mais de três ordens de grandeza menor via PCM, quando comparado ao MC. No entanto o PCM pode requerer um número elevado de pontos de colocação, em casos em que a função de interesse apresentar muitas variáveis aleatórias. Além disso, foi observado que o PCM pode não obter soluções satisfatórias em casos de funções com pontos de singularidades ou comportamento multimodal.

•

Os resultados obtidos mostram a grande vantagem em usar o PCM para os problemas de otimização na engenharia aqui considerados, i.e. para funções suaves e com poucas variáveis aleatórias.

•

A combinação das várias técnicas de aproximação descritas, além de métodos eficientes para lidar com problemas multiobjetivo, permitiram a obtenção das soluções OMR em pouco tempo computacional. Em alguns casos a eficiência alcançada pelo uso do modelo aproximado foi de um tempo computacional 100 vezes menor que sem o seu uso.

6.3 TRABALHOS FUTUROS Baseado nos estudos/implementações aqui conduzidos e resultados obtidos, apresenta-se as seguintes sugestões de continuidade do presente trabalho: •

Abordar outros tipos de problemas, como problemas dinâmicos transiente e/ou interação fluido-estrutura e/ou escoamentos em meios porosos, etc.

•

Incorporar processos de adaptação de malhas, para, não apenas garantir a acurácia da solução e diminuir o tempo da análise via MEF, como o tempo do estágio “off-line” no MBR e o tempo do procedimento “on-line” para variáveis de projeto não mapeáveis. Lembrando que o procedimento “on-line” do MBR independe da malha de EF considerada.

•

Adicionar a opção de elementos de maior ordem, lembrando que o procedimento “on-line” do MBR também independe tipo de elemento considerado.

•

Adaptar o MBR para análises que considerem as não linearidades físicas e geométricas do problema.

•

Adaptar o MBR para otimização de forma. 120 

 

•

Implementar o NNC com a expansão do Plano Utópico, para encontrar frentes de Pareto com mais de dois objetivos.

•

Investigar problemas com mais de 3 objetivos.

•

Aplicar o PCM para variáveis aleatórias com pdfs diferentes.

•

Implementar o ME-PCM (“Multi-Element Probabilistic Collocation Method”), objetivando tratar os casos de funções multimodal ou com singularidades.

•

Implementar a integração do PCM por grades esparsas (“sparse grids”) para diminuir o número de pontos de integração em problemas com muitas variáveis aleatórias.

•

Explorar novos problemas de OMR. Considerando outras estruturas, funções objetivo, restrições, etc.

•

Considerar restrições probabilísticas pela probabilidade de falha, através de análise de confiabilidade estrutural.

•

Utilizar paradigmas da computação paralela, que pode ser eficientemente incorporada em várias etapas do processo, como: para conduzir o estágio “off-line” do MBR, para a geração dos pontos de Pareto e para o cálculo das grandezas estatísticas.

•

Automatizar, por meio de análise de erro, o cálculo do tamanho das amostras para a base do MBR, bem como o número de pontos de integração para o cálculo das estatísticas via MC e PCM.

•

Investigar procedimento para determinar o número adequado de pontos de Pareto, necessários para a definição da fronteira de Pareto.

121