Embed Size (px)
Citation preview
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 Per-nambuco, como parte dos requisitos necessários à obtenção do Grau de Mestre em Ciências em En-genharia 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.
M921o Motta, Renato de Siqueira
Otimização robusta de estruturas utilizando o método da base reduzida / Renato de Siqueira Motta. - Recife: O Autor, 2009.
xiii, 121 f.; il., gráfs., tabs. Dissertação (Mestrado) – Universidade Federal de Pernambuco.
CTG. Programa de Pós-Graduação em Engenharia Civil, 2009. Inclui Referências bibliográficas. 1. Engenharia civil. 2. Otimização Robusta. 3. Método da
Colocação Probabilística. 4. Otimização Multiobjetivo. 5. Método da Base Reduzida. I. Título.
UFPE 624 CDD (22. ed.) BCTG/2010-043
iv
“No processo de evolução das espécies, somos um ótimo local”
v
OTIMIZAÇÃO ROBUSTA DE ESTRUTURAS UTILIZANDO O MÉTODO DA BASE REDUZIDA
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 proje-tos ó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 mo-delos 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 óti-mo é chamado Otimização Robusta (OR).
Aqui, abordaremos duas técnicas para a análise de propagação de incerteza, não in-trusivas, 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 mo-mentos estatísticos da função de interesse. Várias medidas de robustez têm sido propos-tas 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 abor-dagens 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 cus-tosas, 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 sepa-rabilidade 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 aproxi-mação descritas permitiu a obtenção das soluções OMR em pouco tempo computacio-nal.
Palavras chaves: Otimização robusta, método da Colocação Probabilística, Otimiza-
ção Multiobjetivo, Método da base reduzida.
vi
STRUCTURAL ROBUST OPTIMIZATION CONSIDERING REDUCED-BASIS METHOD
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 ad-dressed in this work and the performance of the different strategies discussed are com-pared. 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 opti-
mization, Reduced basis method
SUMÁRIO LISTA DE SÍMBOLOS ix 1. INTRODUÇÃO 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 8 2.1. FORMULAÇÃO DE PROBLEMAS TÉRMICOS 8
2.1.1. Método dos Elementos Finitos aplicado à problemas térmicos 9
2.2. FORMULAÇÕES DOS PROBLEMAS ELÁSTICOS 11
2.2.1. Método dos Elementos Finitos aplicado à elasticidade bidimensional 13
2.2.2. Acoplamento termo-elástico 14
2.3. MAPEAMENTO DAS EQUAÇÕES DO MEF 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 20
2.4.1. Procedimento Computacional Off-Line/On-Line 23
2.5. EXEMPLOS 24
2.5.1. Placa quadrada com orifício central 24
2.5.2. Exemplo termo-elástico acoplado 28
2.6. REFERÊNCIAS 32
3. OTIMIZAÇÃO ESCALAR E ANÁLISE DE SENSIBILIDADE 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 42
3.5.1. O método da base reduzida no procedimento de otimização 43
(a) Análise de sensibilidade através do RBM 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 77
5. OTIMIZAÇÃO CONSIDERANDO INCERTEZAS 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 115
6. CONCLUSÃO E TRABALHOS FUTUROS 118 6.1. PRINCIPAIS DESENVOLVIMENTOS 118
6.2. CONCLUSÕES DOS RESULTADOS OBTIDOS 118
6.3. TRABALHOS FUTUROS 120
x
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 constan-te.
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áti-ca.
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 Normali-zada
NC “Normal Constraint” - Restrição Normal
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
1c Constante de Lamé
2c Constante de Lamé
ijklD Componente do tensor de elasticidade
E Módulo de elasticidade do material Fobj(x) Função objetiva gi(x) Restrição de desigualdade hi(x) Restrição de igualdade kij Termos do tensor de condutividade térmica N Tamanho da Base
iN Função de forma do nó i
( )f μ Flexibilidade
( )Nf μ Flexibilidade aproximada
L(x,α,β,λ) Função Lagrangeana ti Largura da região Ti Temperatura no nó i
irx Coordenada de um nó i na sub-região r lkx Limite inferior da variável de projeto
xi
ukx Limite superior da variável de projeto
Vetores b Forças de volume agindo sobre o domínio V*
br(μ) Forças de volume agindo sobre o domínio computacional para uma sub-região r
ed Vetor de deslocamentos nodais de um elemento ε Vetor de deformação
nf Forças normais aplicadas no contorno Γ*
tf Forças tangenciais e tangentes aplicadas no contorno Γ* rnf (μ) Forças normais aplicadas no contorno Γ para uma sub-região r rtf (μ) Forças tangenciais aplicadas no contorno Γ para uma sub-região r r
TEF Carregamento devido à deformação térmica F Vetor de forças, Vetor das funções objetivo F* Vetor do mínimo individual das funções objetivo
)(μF N Vetor de carregamentos transformado r*F Vetor de carregamentos do domínio real de uma sub-região r
rjF Vetor de carregamento independente de μ para uma sub-região r
F̂ Vetor de Pseudo-força q Fluxo de calor no contorno qv Calor interno gerado sobre o domínio V*, S Vetor da direção de busca
NT Vetor de temperaturas aproximado u Vetor de deslocamento
)(μNu Deslocamento aproximado
0T Temperatura inicial
Matrizes
e*B Matriz que relaciona as componentes da derivada de deslocamento com deslocamentos
)(μC Operador simétrico usado para relaxação D Matriz de elasticidade
*D Pseudo-matriz de elasticidade no domínio real G Matriz de transformação
rG Matriz de transformação de uma sub-região r H Hessiana da função Lagrangeana I Matriz identidade
sK Matriz de rigidez *sK Matriz de rigidez no domínio real
xii
*NsK Matriz de rigidez no domínio real no espaço reduzido *rsK Matriz de rigidez do domínio real de uma sub-região r rs jK Matriz de rigidez independente de μ para uma sub-região r
TK Matriz de “rigidez térmica” ou condutividade térmica generalizada *TK Matriz de “rigidez térmica” transformada *NTK Matriz de “rigidez térmica” transformada no espaço reduzido rk jK Matriz de condutividade independente de μ r
cΓK Matriz de convecção no contorno de referência Nr
cΓK Matriz de convecção no contorno de referência no espaço reduzido r
cΩK Matriz de convecção no domínio de referência Nr
cΩK Matriz de convecção no domínio de referência no espaço reduzido Kk Matriz de condutividade térmica N Matriz de função de forma
NS Conjunto de amostras para campo de solução WN Espaço da Base Reduzida Z Matriz compostas por soluções base do espaço reduzido
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
( )rs j
β μ Termos que dependem do parâmetro μ da transformação da Matriz de ri-gidez
( )rT j
β μ Termos que dependem do parâmetro μ da transformação da Matriz de condutividade térmica
ixΔ Pertubação da variável
λ, β, φ Multiplicadores de Lagrange Γ Contorno μ Parâmetros variáveis
lowμ Limite inferior do espaço de projeto upμ Limite superior do espaço de projeto rϕΩ ( )μ Termos da transformação no domínio que dependem de μ rϕΓ ( )μ Termos da transformação no contorno que dependem de μ
xiii
Ω Domínio *Ω Domínio real (transformado)
ijσ Tensão ν Coeficiente de Poisson γ Coeficiente de expansão térmica Vetores α(μ) Coeficiente linear da aproximação no espaço reduzido
*0ε 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
i∇N 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
xiv
1
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. Otimiza-dores de propósito geral não resolvem tais problemas. Uma classe de estratégias basea-das 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 (“Normal-Boundary Intersection”) (DAS e DENNIS, 1996) e o NNC (“Normalized Normal Cons-traint”) (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 proble-ma 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 de-vido à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).
2
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 pro-pagação de incerteza, o método de Monte Carlo (MC) e o método da colocação probabi-lí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 computacional-mente 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ásti-cos 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 reque-ridos pelo otimizador. Este algoritmo usa a programação sequencial quadrática (“Se-quential 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 pon-tos 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 redu-zidas (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 requeri-das 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 efi-ciência computacional.
No MBR utiliza-se do mapeamento da geometria das regiões do modelo para apli-car 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
3
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 ra-pidamente. 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 procedi-mentos 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 procedi-mentos 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 dife-rentes 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 sen-sibilidade 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 anterio-res, 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 estru-tural;
• Realizar um estudo comparativo através exemplos envolvendo otimização multiob-jetivo;
4
• Implementar rotina para o cálculo dos parâmetros estatísticos da resposta de inte-resse, 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 va-riá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 gradi-entes;
• 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 multiobjeti-vo, 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.
5
Figura 1.1 Algoritmo para as diferentes opções de otimização consideradas.
Caso Estocástico:
Definir variáveis aleatórias e suas
distribuições
Otimização Escalar
ou
Otimização Multiobjetivo (neste caso,
repetir o pro-cesso para cada ponto de Pare-
to)
Determinístico
ou
Estocástico
MEF
ou
MBR
Caso MBR:
Definir amostras
Não
Sim
Fim
Análise Estru-tural/Térmica
Geração do projeto inicial
Início
Geração do novo projeto
Otimizador SQP
Análise de Sensibilidade
Convergiu?
Definir: • Geometria do problema • Condições de contorno • Variáveis de projeto • Funções objetivo e
restrição
Caso MBR: Gerar soluções via MEF e a estruturas de dados off-line
6
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 formu-lação do problema bidimensional térmico e estrutural, o MEF aplicado a estes, a separa-bilidade 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 procedi-mentos 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 a-bordagens.
No capítulo quarto é apresentada a formulação do problema de otimização multiob-jetivo, 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 con-duzidas, 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 es-traté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ísti-cos 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 oti-mizaçã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. “Struc-tural Analysis and Optimization in the Framework of Reduced-Basis Method”. Struc-tural 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 Re-duced 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.
7
ALBUQUERQUE T. M. M.. “Análise e Otimização de Problemas Térmicos e Estrutu-rais Bidimensionais Através do Método das Bases Reduzidas”. Dissertação de mestra-do, 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 Applica-tions of Finite Element Analysis”. Fourth edition, Wiley, 2003. DAS, I. e DENNIS, J.E. “Normal Boundary Intersection: A New Method for Gene-rating Pareto Surface in Nonlinear Multicriteria Optimization Problems”. SIAM J. Oti-mization, Vol. 8, No. 3, pp. 631-657, 1996. MACEDO, C. M. H. “Otimização de Treliças Planas sob Várias Solicitações com Ênfa-se 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 un-certainties in black box systems”. Thesis (M.S. in computer science). Washington State University. 2005. VANDERPLAATS, G. N. “Numerical Optimization Techniques for Engineering De-sign - with Applications”. 4rt Edition, Vanderplaats Research & Development, Inc., Colorado Springs, CO, 2004.
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, satis-faz a equação diferencial parcial
Ω
(2.1) ( )01
ndi
i i
q Q T T emx
αΩ=
∂= − − Ω
∂∑ (2.1)
onde representa a componente i do vetor de fluxo de calor, nd é o número de dimen-sões do problema, representa o termo de fonte (ou sumidouro) de calor no domínio,
iqQ
αΩ é o coeficiente de convecção no domínio, é a temperatura ambiente e Ω repre-senta o domínio do problema. O fluxo de calor se relaciona com o gradiente de tempera-tura através da lei de Fourier da condução,
0T
(2.2)1
nd
i ijj j
Tq kx=
∂= −
∂∑ (2.2)
onde kij representa os termos do tensor da condutividade térmica. A Equação (2.1) pode então ser expressa como
(2.3) ( )01 1
nd nd
iji ji j
Tk T T Q em Ω x x
αΩ= =
⎛ ⎞∂ ∂− + − =⎜ ⎟⎜ ⎟∂ ∂⎝ ⎠∑ ∑ (2.3)
Além disso, as seguintes condições de contorno são consideradas
(2.4) ( )0
em T
n n
T T
T T emαΓ
= Γ
= − − Γq q q
(2.4)
8
onde representa a porção de contorno com temperaturas prescritas TΓ T (condição de contorno de Dirichlet), representa a parcela do contorno sujeita a fluxo e/ou a con-vecção (condição de contorno de Cauchy), é o fluxo na direção n ortogonal à
qΓ
nq qΓ ,
αΓ é o coeficiente de convecção no contorno e é a temperatura ambiente. 0T
As expressões foram escritas até então, utilizando notação indicial. Na forma matri-cial, as equações governantes no domínio Ω (2.3) e no contorno qΓ (2.4) podem ser expressas como
(2.5) 0
0
T
en q
T T T Q em
T T T em
α α
α αΩ Ω
Γ Γ
−∇ ∇ + = + Ω
∇ + = + Γ
k
k n q (2.5)
Para o caso bidimensional ortotrópico, [ ]/ , / Tx y∇ = ∂ ∂ ∂ ∂ e k é a matriz constitutiva que contém as condutividades térmicas do material nas direções x e y, dado por
(2.6)0
0x ⎤
⎥ y
kk
⎡= ⎢⎣ ⎦
k (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 as-sociados à temperatura em determinado nó, a temperatura Te(ξ,η) no elemento “e”, em qualquer ponto (ξ,η) deste elemento, é aproximada da seguinte forma
(2.7) ( ) ( )1
, ,n
e e ei i
iT N Tξ η ξ η
=
=∑ N Te e (2.7)
onde é a temperatura no nó i e é a função de forma associado ao nó i. Conside-rando
eiT e
iN( , )j jξ η as coordenadas locais do nó j as funções de forma são definidas de tal
forma que (ZIEKIEWICZ e TAYLOR, 2000)
(2.8) ( ) 1, para , , 1...
0, para i j j
j iN
j iξ η
=⎧= ⎨ ≠⎩
i j n= (2.8)
9
Definimos, então, os vetores e , onde eN eT 1 , , , ,e e eiN N N e
n⎡ ⎤= ⎣ ⎦N … …
, , , ,Te e e
i nT T
é o vetor das
funções de forma do elemento “e”, 1e T⎡ ⎤= ⎣ ⎦… …T é 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 do-mínio Ω e para o contorno , podem ser aproximadas como qΓ
(2.9)( )( )
0
0
T
n q
T Q em
T em
α α
α αΩ Ω
Γ Γ
−∇ ∇ + = + Ω
∇ + = + Γ
k N N T
k Nn N T q (2.9)
onde N é o vetor contendo as funções de formas relacionadas aos nós, e T é o vetor con-tendo 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 pro-blema, este método consiste em encontrar os parâmetros referentes à aproximação con-siderada (no nosso caso T), que anule a integral em todo o domínio, do resíduo das e-quações ponderado por diferentes funções peso , 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
lW
(2.10)( ) ( )
( ) ( )0 0
q
q
Tl l
l n l q
W d W d
W d W d
α α
α α
Ω ΓΩ Γ
Ω ΓΩ Γ
⎡ ⎤⎢ ⎥−∇ ∇ + Ω+ ∇ + Γ × =⎢⎣
+ Ω+ + Γ
∫ ∫
∫ ∫
k N N k Nn N T
T Q T q
q
=
⎥⎦ (2.10)
para l =1..nl
Considerando uma discretização em elementos finitos do domínio em estudo, apli-cando 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)
(2.11)( )
( ) ( )0 0
q
q
T T Tq
T Tn q
d d
d d
α α
α α
Ω ΓΩ Γ
Ω ΓΩ Γ
⎡ ⎤⎢ ⎥∇ ∇ + Ω+ Γ ×⎢⎣
+ Ω+ + Γ
∫ ∫
∫ ∫
N k N N N N N T
T Q N T q N
⎥⎦ (2.11)
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) T T=K T F (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 é dado por TF
(2.13) T Ω= +F F FΓ (2.13)
onde representa as cargas no domínio e ΩF ΓF as cargas no contorno, dado por
(2.14) ( ) ( )0 0;
q
Tndα αΩ Ω Γ Γ
Ω Γ
= + Ω = +∫ ∫F T Q N F T q NTqdΓ (2.14)
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
e convecção contorno , i.e. cΩK cΓK
(2.15) T k c= + +K K K K cΩ Γ (2.15)
onde
(2.16) ; ;q
T Tk c cd dαΩ Ω Γ Γ
Ω Ω Γ
= ∇ ∇ Ω = Ω = Γ∫ ∫ ∫K N k N K N N K N NTqd α (2.16)
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 de-formaçã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 ade-quadamente a deformação. Desta forma as equações de equilíbrio são satisfeitas na for-ma
(2.17)1
0nd
ijj
i j
b emxσ
=
⎛ ⎞∂+ = Ω⎜ ⎟⎜ ⎟∂⎝ ⎠
∑ (2.17)
onde bj é a força volumétrica agindo sobre o domínio Ω na direção j e o tensor de ten-sõ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) 0 T∇ + =σ b (2.18)
onde é um operador diferencial, σ é o vetor das tensões e b o vetor das forças de volume, definidos respectivamente por
∇
(2.19)0 0
, ,0 0
xx
yy xT
xy y
yx
bx yb
y x
σσσσ
⎡ ⎤∂ ∂⎡ ⎤⎢ ⎥⎢ ⎥ ⎡ ⎤∂ ∂ ⎢ ⎥⎢ ⎥∇ = = = ⎢ ⎥⎢ ⎥∂ ∂⎢ ⎥ ⎣ ⎦⎢ ⎥⎢ ⎥∂ ∂⎣ ⎦ ⎢ ⎥⎣ ⎦
bσ (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
(2.21)( )( )
,,
u x yv x y⎡ ⎤
= ⎢⎣ ⎦
u ⎥ (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).
(2.22)
2 1 1
1 2 1
2 2
2 2
2 02 0
0 00 0
c c cc c c
c cc c
⎡ ⎤+⎢ ⎥+⎢=⎢⎢ ⎥⎢ ⎥⎣ ⎦
D
00 ⎥⎥
(2.22)
onde, e são as constantes de Lamé definidas por 1c 2c
(2.23) 1 2,(1 )(1 2 ) 2(1 )
E Ec υ cυ υ
= =+ − υ+
(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
(2.24) 0 T∇ ∇ + =D u b (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 , ou seja, NΓ
(2.25)emem
D
N
= Γ= Γ
u uσ σ
(2.25)
2.2.1 Método dos Elementos Finitos aplicado à elasticidade bidimensional
Na aproximação via elementos finitos para os problemas bi-dimensionais onde to-dos 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, respectiva-mente, as componentes horizontais e verticais do deslocamento, são aproximados da seguinte forma
(2.26)( ) ( )
( ) ( )
1
1
, ,
, ,
ne e
i iin
e ei i
i
u N
v N
ξ η ξ η
ξ η ξ η
=
=
∑
∑
e
e
u
v (2.26)
onde e são as componentes do deslocamento no nó i do elemento e e é a fun-ção de forma associado ao nó i, da mesma forma como foi definido para o caso térmico na seção
eiu e
iv iN
2.1.1. Com isso pode-se reescrever a equação que relaciona deformação e des-locamento para um elemento genérico e da seguinte maneira
(2.27)1
ne e e e
i==∑ε B d B de
i i (2.27)
onde Be é matriz que relaciona as deformações com deslocamentos, e n é o número de nós por elemento. Para um nó específico i, pode ser escrita como
[ ]Te e ei i iu v=d
eiB
(2.28) e e
i iN= ∇B (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 pro-blemas de elasticidade pode ser escrita como
13
(2.29)n t
T T Tnd d d
Ω Ω Γ Γ
⎡ ⎤Ω = Ω+ Γ + Γ⎢ ⎥
⎣ ⎦∫ ∫ ∫ ∫B DB u N b N f N fT
t d (2.29)
Na Equação (2.29) b é o vetor das forças de volume agindo sobre o domínio , e são respectivamente os vetores das forças normais e tangentes aplicadas no contor-no e respectivamente, N é a matriz relacionada com as funções de forma utiliza-das pelo Método dos Elementos Finitos em modelos bi-dimensional e B a matriz que relaciona as componentes das deformações com deslocamentos, apresentadas anterior-mente.
Ω nf
tfΓn tΓ
A Equação (2.29), pode ser reescrita como
(2.30) s s=K u F (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) 0d T T T T Tn td d d d
Ω Ω Γ Γ Ω
⎡ ⎤Ω = Ω+ Γ + Γ + Ω⎢ ⎥
⎣ ⎦∫ ∫ ∫ ∫ ∫B DB u N b N f N f B Dε (2.31)
Na equação anterior, é o vetor das deformações iniciais devido à variação térmi-ca. Assim, o último termo da Equação
0ε(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ér-mica do material em uma estrutura e devido a restrições oriundas das condições de a-poio.
Considerando um material isotrópico com coeficiente de expansão térmica γ, coefi-ciente 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
(2.32)
( )( )
0
11
00
TT
ν γν γ
+ Δ⎡ ⎤⎢ ⎥+ Δ⎢=⎢ ⎥⎢ ⎥⎣ ⎦
ε ⎥ (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 é introduzi-do 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 defor-mações iniciais térmicas, que introduzem carregamentos mecânicos na estruturas atra-vé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 – “Cons-tant 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), escre-vendo 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 1r 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 carre-gamento definidos nas Equações (2.12) e (2.30) podem ser reescritos como
(2.33)1 1
R Rr
r r= =
= =∑ ∑K K F rF (2.33)
Separando os termos (regiões) que dependem do parâmetro variável, um novo con-junto de parâmetros constituirá somente um novo dado a ser ajustado na forma decom-posta, 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 proce-dimento se mostrará ainda mais interessante mais adiante, ao se utilizar o MBR (Méto-do 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 e o domínio de referência , 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
*Ω*rΩ
rΩ
(2.34) * *ou com , 1, 2, e 1,...,r r r r r r
i ij jx G x i j r R= = =x G x = (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
(2.35)
* *1 1
1 2* *
x
2 2
1 2
( )r
xx x
x xx x
⎡ ⎤∂ ∂⎢ ⎥∂ ∂⎢=⎢∂ ∂⎢ ⎥∂ ∂⎣ ⎦
μ ⎥⎥
G (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 gover-nantes 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) * *( ) e ( )r r r r rd d dϕΩΩ = Ω Γ = Γμ rdΓ μ ϕ (2.36)
onde
(2.37) 1 1)] r t−( ) det[ ( )] e ( ) [ (r r r rϕ ϕ−
Ω Γ= =μ G μ μ G μ n (2.37)
onde é o vetor unitário tangente ao contorno Γr de uma região particular r e *r tn 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 apli-cada á matriz e aos infinitesimais ∇N e d 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
(2.38)
( )
( )
( )
* 1
* 1
* 1
( ) ( ) det[ ( )]
det[ ( )]
[ ( )]
r
r
rq
r T rT r r rk
r r T r rc
r r T r r t rc q
d
d
d
α
α
−
Ω
−Ω Ω
Ω
−Γ Γ
Γ
⎡ ⎤= ∇ ∇ Ω⎣ ⎦
= Ω
= Γ
∫
∫
∫
K μ N G μ k G μ N G μ
K μ N N G μ
K μ N N G μ n
r
(2.38)
16
(2.39)( )
( )
* 10
* 1 r t rΓ
0
det[ ( )]
[ ( )]
r
rq
r r r T r r
r r r T rn q
d
d
α
α
−Ω Ω
Ω
−Γ Γ
Γ
= + Ω
= −
∫
∫
F T Q N G μ
F T q N G μ n(2.39)
Na presente aplicação, o único termo dependente de μ é ( )rG μ .
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 matri-zes e vetores decompostos apresentados nas Equações (2.38) e (2.39) com o auxilio da Equação (2.37), da seguinte maneira
(2.40) ( ) ( ) ( ) ( )* * *
1
; ( ) ;nt
r r r r r r r rk T k j c c cj
j
β ϕ ( ) rcϕΩ Ω Ω Γ Γ
=
= =∑K Γ=μ μ K K μ μ K K μ μ K (2.40)
(2.41) ( )* ( ) ( )r r r r
T ϕ ϕΩ Ω Γ= +F μ μ F μ FrΓ (2.41)
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 são formadas a partir da sub-decomposição da
matriz de condutividade térmica. Nesta decomposição nt e
rk jK
*rkK ( )r
T jβ μ variam de a-
cordo com a transformação específica da região. As matrizes são calculadas da se-guinte forma
rk jK
(2.42) , 1...
r
r T r rk j j d j
Ω
= ∇ ∇ Ω =∫K N k N nt (2.42)
onde as sub-matrizes constitutivas são obtidas através da decomposição de rjk * ( )rk μ ,
que é definida a partir do seu valor de referência como
(2.43) * ( ) ( ) ( ) det[ ( )]r rT r r r 1−⎡ ⎤= ⎣ ⎦k μ G μ k G μ G μ (2.43)
A partir desta transformação (Equação (2.43)) são retiradas não só as matrizes como
também o parâmetro
rjk
( )rk j
β μ que contêm os fatores de dependência da variável µ, de tal
modo que
(2.44) ( ) rjj
*
1
( )nt
r rT
j
β=
=∑k μ μ k (2.44)
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) * * *
T T=K T F (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), aplica-das à matriz B, as relações infinitesimais da Equação (2.36) e o conceito da separabili-dade, tem-se as novas sub-matrizes de rigidez e os termos de carregamento no domínio real de cada região (Equações
*reK *r
eF(2.29) e (2.30)), da seguinte forma
(2.46) ( ) ( ) ( )* 1( ) ( )det[ ( )]
r
r rT T r r rs d−
Ω
= ∫K μ G μ B D G μ B G μ rΩ (2.46)
(2.47)
( )* 1
1
det[ ( )] [ ( )]
[ ( )]
r nr
tr
r T r r r T r r nrts n
T r r tr t trt
d d
d
− −
Ω Γ
−
Γ
= Ω +
+ Γ
∫ ∫
∫
F μ N b G μ N f G μ n
N f G μ n
1 nrΓ
(2.47)
Assim como para o caso térmico, as componentes que dependem de μ podem ser retirados das integrais e desta forma obter
(2.48) ( ) ( )*
1
ntr r r
s sjβ
=
=∑K jjμ μ K (2.48)
(2.49) ( )* ( ) ( )r r r r r
s bϕ ϕΩ Γ= +F fμ μ F μ F (2.49)
onde
(2.50)
r
r T rj j
V
dV= ∫K B D B r (2.50)
(2.51) ,
r nr tr
r T r r r T r r T rb f n
V
dV d dΓ Γ
= = Γ +∫ ∫ ∫F N b F N f N f rt Γ (2.51)
18
Nas equações anteriores, cada componente das sub-matrizes de elasticidade e rjD
( )rs j
β μ são definidos de maneira análoga à decomposição da matriz de condutividade
térmica (Equação *rk (2.43) e (2.44)), de tal modo que
(2.52) ( )* 1
1
( ) ( ) ( ) det[ ( )]nt
r r r rT r r rs jj
j
β −
=
⎡ ⎤= = ⎣ ⎦∑D μ μ D G μ D G μ G μ (2.52)
Como se pode verificar na Equação (2.52), dependendo da transformação a ser a-plicada 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 nova-mente 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
(2.53) ( )
*
* * *
r
r TTE o
V
dV= ∫F μ B Dε μ *( ) r (2.53)
O vetor carregamento que acopla os efeitos térmicos, pode ser decomposto de for-ma 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 i-niciais de uma região r para um valor de μ qualquer, pode ser calculado como
(2.54) ( ) ( )* ( ) ( )r r r r r r
s b fϕ ϕΩ Γ= + +F *TEμ μ F μ F F μ (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 , no domínio real, resolve-se o novo sistema.
*u
(2.55) * * *
s s=K u F (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 deslo-camentos ou temperaturas) faz parte do espaço de solução de dimensões infinitas , e claramente os possíveis valores de u(μ) não cobrem inteiramente o espaço . Conside-rando 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
YY
YFigura 2.1 (a). A partir disto, po-
19
de-se afirmar que o espaço das soluções não se trata de um espaço de dimensão infinita e sim um espaço de dimensão reduzida , onde é observado a dependência do pa-
râmetro variável μ (PRUD` HOMME et al, 2002), como ilustrado na Y NY
Figura 2.1 (b).
(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 de-finição do conjunto de amostras representada da seguinte forma
(2.56) ( ) ( ){ }NRR
NS μμμμ ,,,,,, 11
1 ………= (2.56)
onde cada está contido no espaço Dμ, isto significa que ( iRμμ ,,1 … )
(2.57) upir low μμμ ≤≤ (2.57)
Y
u(μ) u(μN)
u(μ1)
u(μnovo)
u(μ2)
20
no qual e são os limites inferior e superior, respectivamente, do espaço de projeto Dμ e N é o número de amostras.
lowμ upμ
Para cada componente de SN calcula-se a solução u(μ) através do Método dos Ele-mentos 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) NW = span{u[ ( )11 ,, Rμμ … ],…,u[ ( )N
Rμ ]} μ ,,1 … (2.58)
Para simplificar a notação defini-se como iζ nR∈
(2.59) ( )( )iR
i μμ ,...,1uζ = (2.59)
e assim WN pode ser escrito como
(2.60) { } NiomspanW iN ,,1c …== ζ (2.60)
A equação acima significa que qualquer vetor do espaço reduzido WN pode ser ex-presso como uma combinação linear das soluções . Vale salientar que as soluções nos pontos amostrais devem ser linearmente independentes, para que
iζζ 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
(2.61) ( )1
( )N
jN j
j
j nRα α=
=∑u μ μ ζ ∈ (2.61)
ou na forma matricial
(2.62) ( )μZα μu =)(N (2.62)
no qual
(2.63) ( ) ( ) ( )11[ ,..., ] ; [ ,..., ]NN T α α= =Z ζ ζ α μ μ μ (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) ( ) ( ) ( )μFμZαμK = (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) ( ) ( ) ( )μFZμZαμKZ TT = (2.65)
Com isso, os coeficientes lineares α(μ) são obtidos resolvendo o problema no espa-ço WN dado da seguinte forma
(2.66) )( )()( μFμαμK NN = (2.66)
no qual
(2.67) ( ) NTNNNTN RR ∈=∈= × )()(,)( μFZμFZμKZμK (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
(2.68)1 1
,R R
N Nr N
r r= =
= =∑ F FNr∑K K (2.68)
com e 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.
NrK NrF
(2.69) ( ) )()(,)( μFZμFZμKZμK rTrNrTrN == (2.69)
Nas equações acima os termos e podem ser mapeados como foi visto nas seções
)(μK r )(μFr
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
(2.70) ( ) ( )*
1
ntNr r Nr
s s sj
β=
= ∑K μ μ K jj; ( )* ( ) ( )Nr r Nr r Nr
s bϕ ϕΩ Γ= +F μ μ F μ Ff (2.70)
onde apenas os termos rs jβ , , e r r
b n trϕ ϕ ϕ dependem do parâmetro variável, Nr
s jK e 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 men-cionado 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
NrjF
2.3. Para o caso térmico, as matrizes e vetores de cada região podem ser escritos da seguinte forma
(2.71) ( ) ( ) ( ) ( )* * *
1; ( ) ;
ntNr r Nr Nr r Nr Nr r Nr
k T k j c c cjjβ ϕ ( ) cϕΩ Ω Ω Γ Γ
=
= =∑K μ μ K K μ μ K K μ μ KΓ = (2.71)
22
(2.72) ( )* ( ) ( )Nr r Nr r Nrd d c cϕ ϕ= +F μ μ F μ F (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 de-pendentes/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 gran-dezas pelo Método da Base Reduzida é conduzida por um algoritmo “off-line” (inde-pendente 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
e vetores de força . Consequentemente, no estágio on-line utiliza-se estes da-dos gerados anteriormente para executar uma resposta em tempo real para um novo μ.
NriK Nr
iF
No caso de problemas de termo-elasticidade, o procedimento do Método da Base Redu-zida é 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 procedi-mentos 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: ) (({ ) }1N1 1,..., ,..., ,..., N
R Rμ μ μ μ=S ;
1.2 Construir a matriz de soluções térmicas via MEF: 1 ; [ ,..., ]NT T T=Ζ ζ ζ
1.3 Construir as componentes matrizes de “rigidez térmica” no espaço reduzido:
; ;Nr T r Nr T r Nr T rk j T k j T c T c T c T c TΩ Ω Γ Γ= = =K Z K Z K Z K Z K Z K Z ;
1.4 Constrói as componentes dos vetores de cargas térmicas no espaço reduzido: Γ . ;Nr T r Nr T r
T TΩ Ω Γ= =F Z F F Z F
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): 1[ ,..., ]N
s s s=Ζ ζ ζ ;
2.3 Construir as componentes da matriz de rigidez no espaço reduzido: Nr T rs j s s j=K Z K sZ
f
;
2.4 Construir as componentes dos vetores de carga estrutural no espaço reduzido: . ;Nr T r Nr T r
b s b f s= =F Z F F Z F
23
ON-LINE – para um novo vetor μ :
1- Problema Térmico 1.1 Formar a matriz de “rigidez térmica” no espaço reduzido:
k jj; ( ) ( )*
1 1( ) ( )
R ntN r Nr Nr r Nr
T c c Tr j
ϕ ϕ βΩ Ω Γ Γ= =
⎛ ⎞= + +⎜ ⎟
⎝ ⎠∑ ∑K μ μ K μ K μ Kr
1.2 Formar o vetor de cargas térmicas no espaço reduzido: Nr( ) ( )*
1( ) ( )
RN r Nr r
Tr
ϕ ϕΩ Ω Γ Γ=
= +∑F μ μ F μ F ;
1.3 Resolver: * *( ) ( ) ( )N NT T T=K μ μ F μα ;
1.4 Calcular: ( ) ( )NT T=T μ Z μα .
2- Elástico Acoplado
2.1. Formar a matriz de rigidez no espaço reduzido: * ( ) ( )R nt
N r Nrs s j s
r 1 j 1β
= =
=∑∑K μ μ K j
TE
;
2.2 Formar o vetor de carga estrutural no espaço reduzido:
( )* *( ) ( ) ( )R
N r Nr r Nr Nrs f
r 1ϕ ϕΩ Ω Γ
=
⎡ ⎤= + +⎣ ⎦∑F μ μ F μ F F μ ;
2.3 Resolver: ; * *( ) ( ) ( )N Ns s s=K Fμ α μ μ
2.4 Calcular: ( ) ( )Ns s=u Zμ α μ ;
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ên-cia foi estudada em relação à flexibilidade total da estrutura (FTu). Os elementos utili-zados 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 ele-mentos (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 pos-suem lados iguais (quando possível) e as malhas de elementos triangulares se sobre-põem as malhas de elementos quadrilaterais, ou seja, possuem os mesmos nós. A densi-dade (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 à flexi-bilidade 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 defor-mado, utilizando as discretizações mencionadas, distribuições similares podem ser ob-servadas. 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 con-figuraçã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 consi-derando 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 proce-dimentos de otimização e cálculos estocásticos, os quais serão considerados nesta dis-sertaçã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 indi-cada na Figura 2.5 (a) e (b) com o total de treze (13) regiões, especificadas para o pre-sente problema, numeradas conforme indicado na Figura 2.5 (a) em função dos parâme-tros 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â-metros μ=( t1, t2, L1, L2 ) pertençam ao espaço de projeto Dμ = [ ] [ ]20,5;5 10;21× ×
cm4. [28;38]
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, totali-zando 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 e-lementos. Os elementos utilizados foram os mesmos do exemplo anterior: o quadrilate-ral (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 repre-senta 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 e-lementos finitos consideradas, foram as que apresentaram um erro (em relação à flexibi-lidade 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 elemen-tos 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 flexibi-lidades 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 calcula-da pelo MBR em função do tamanho da base N utilizada, para as malhas citadas anteri-ormente. 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 amos-tra, 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 du-rante 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 Re-duced 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 Estrutu-rais Bidimensionais Através do Método das Bases Reduzidas”. Msc. Thesis ( in Portu-guese), Dept. de Engenharia Civil, Universidade Federal de Pernambuco, Recife-PE, Brasil, 2005.
ALMEIDA, F. P. C. “OTIMIZAÇÃO DE FORMA EM PROBLEMAS TERMO-ESTRUTURAIS 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 Applica-tions 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., PA-TERA, 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, Mas-sachusetts Institute of Technology, 2003.
ZIENKIEWICZ O C & TAYLOR R L. “The finite element method. Vol. I. Basic for-mulations 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 dimen-sõ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 substanci-almente aumentada durante as últimas décadas, e um progresso considerável foi atingi-do. 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 sensibili-dade, 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 proje-to, 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 Fini-tos, 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 computa-cional 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 é de-nominado problema de otimização escalar (HAFTKA e GURDAL, 1993). Sua repre-sentação genérica é dada por:
(3.1) min ( )nvp
f∈x
x (3.1)
Sujeito às seguintes condições:
(3.2) ( ) 0 1, ,( ) 0 1, ,
1, ,
k
i
lj j uj
h k neg i nix x x j npv
= =
≤ =
≤ ≤ =
xx
………
(3.2)
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 ( ) 0ig <x e ativas caso ( ) 0ig =x , 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éto-dos 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ão-lineares 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 envolvi-do. Isto é importante para se obter uma otimização confiável e um alto nível de eficiên-cia 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ão-lineares com restrições baseadas no método direto (KLEIBER, 2005). A resolução des-tes problemas pode ser feita por meio dos chamados métodos indiretos. Estes transfor-mam o problema original definido pela Equações (3.1) e (3.2), num problema sem res-triçõ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 li-
nearmente 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 aten-der às condições necessárias de primeira ordem, dadas pelas condições de Karush-Kuhn-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 sen-tido 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:
(3.4) * * * * *
1 1( ) ( ) ( ) 0
ne ni
k k i ik i
f h gλ φ= =
∇ + ∇ + ∇ =∑ ∑x x x (3.4)
(3.5) * 0, 1..i i nφ ≥ = i (3.5)
(3.6) * *( ) 0, 1..i ig iφ = =x 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 res-triçõ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 por-que 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 mul-tiplicador 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):
(3.7) ( )1 1
, ( ) ( ) (ne ni
k k i ik i
L f hλ φ= =
= + +∑ ∑x,λ x xφ )g 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,
(3.8) ( )1 1
, ( ) ( ) ( )ne ni
k k i ik i
L f h gλ φ= =
∇ = ∇ + ∇ + ∇∑ ∑x,λ x xφ 0=x (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 = ( ( )) ∪ ( ( ))j , ( )N h∇ x N g∇ x j J∈ x (3.9) onde representa o núcleo e é o conjunto dos índices das restrições de desi-gualdade ativas, ou seja,
( )N ⋅ ( )J x
(3.10) ( ) { : ( ) 0}jJ j g=x x = (3.10)
(3.11){ }{ }
( ( )) ( ( )) | 0 e ( ) 0
( ( )) ( ( )) | 0 e ( ) 0 , (
T
Tj j j
N h núcleo h h
N g núcleo g g j J
∇ = ∇ = ≠ ∇ =
∇ = ∇ = ≠ ∇ ≤ ∈
x x S S S x
x x S S S x )x (3.11)
37
Assim, as direções viáveis estacionárias S pertencem a ( ( )N h )∇ x e , que são os espaços tangentes das superfícies definidas pelas restrições de igualdade e pelas restrições de desigualdade ativas (
( ( )jN g∇ x )
( ) 0jg =x ), 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) 0,T ≥ ∀ ∈S HS S T (3.12)
ou ainda, podemos ter a condição suficiente para x* ser um mínimo local estrito.
(3.13) 0,T > ∀ ∈S HS 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 Lagrangea-na em relação a x:
(3.14) ( )2 2 2
1 1, ( ) ( )
ne ni
k k i ik i
L f hλ φ= =
= ∇ = ∇ + ∇ + ∇∑ ∑H x,λ x xφ 2 ( )g 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 va-riá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 ( )nvp
f∈S
x + S (3.15)
Sujeito as seguintes condições:
(3.16) ne( ) 0 1, ,( ) 0 1, ,
k
i
h kg i n
= =
≤ =
x + Sx + S
…… i
(3.16)
A função Lagrangeana associada às Equações (3.15) e (3.16) é:
38
(3.17 )( ), ( ) ( ) (T TL f+ = + + + + +x S,λ x S λ h x S g x Sφ φ ) (3.17)
e ótimo de primeira ordem aplicadas ao problema (3.15) e (3.16) fornecem as equações:
(3.18) 0 (3.18)
(3.22) i
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 di-mensão ne e ni, respectivamente. As condições d
( ), ( ) ( ) ( )T TL f∇ + = ∇ + + ∇ + + ∇ + =x S,λ x S λ h x S g x Sφ φ
( )h x S (3.19) 0+ = (3.19)
( )+ ≤g x S 0 (3.20) (3.20)
( )T + =g x Sφ (3.21) 0 (3.21)
0, 1..i i nφ ≥ = (3.22)
(3.23) onde ∇g e ∇h são matrizes Jacobianas cujas linhas são os gradientes de ∇gi e ∇hk, res-pectivamente. 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.27
2( ) ( ) 0L L∇ +∇ =x x S (3.24) (3.23)
( ) ( ) 0+∇ =h x h x S (3.24) (3.25)
( ) ( ) 0+∇ ≤g x g x S (3.25) (3.26)
[ ]( ) ( ) 0T +∇ =g x g x Sφ ) (3.26)
é o gradiente e é a matriz Hessiana da função Lagrangeana, defi-nidos por:
(3.28)
onde (L∇ x) 2 ( )L∇ x
( ) ( ) ( ) ( )T TL f∇ =∇ + ∇ + ∇x x λ h x g xφ (3.27)
(3.29) ⎤⎦2 2 2 2( ) ( ) ( ) ( )L f ⎡ ⎤ ⎡∇ = ∇ + ∇ + ∇⎣ ⎦ ⎣x x λ h x g xφ (3.28)
Utilizando a seguinte notação: e )
2 2( (
ne
hλ⎡ ⎤∇ = ∇∑λ h x) x1
i ii=
⎣ ⎦1
i ii=
⎣ ⎦) 2 2( (ni
gφ⎡ ⎤∇ = ∇∑g x) xφ .
ções de otimalidade do subpro-blem inação da direção de busca S dado por:
As Equações (3.22) e de (3.23) a (3.26) são as condia de determ
39
(3.30) 21min ( ) ( )2ndv
T T
Sf
∈∇ + ∇x S S L x S (3.29)
Sujeito as seguintes restrições
(3.31)( ( )( ( )
+∇ =+∇ ≤
h x) h x Sg x) g x S
00
(3.30)
Este é o esquema básico de linearização que, demonstra-se, é localmente de segun-
da 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 deter-minaçã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, surgi-ram 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 utili-zados 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 es-timá-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 deriva-das destas funções e restrições em relação às variáveis de projeto. Vários são os proce-dimentos 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 preci-
40
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 co-mo 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 com-putacional 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)2 2
2
( )( )( ) ( )2!
ii i i
i ii i
fff f μ δ μμμ μ μ μμ μ
μ∂ + Δ∂+ Δ = + Δ +
∂ ∂ee Δ (3.31)
para algum δ entre 0 e 1.
Truncando-se os termos de ordem superior ( 2iμΔ ) e realizando-se algumas opera-
ções matemáticas, tem-se:
(3.33) ( ) ( )ifii i
i i i
ff f μ μμ μ μ
+ Δ∂ Δ=
∂ Δ Δe μ− , para i = 1,..., m (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 é determina-do 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) FKu = (3.33) onde K é a matriz de rigidez, u é o vetor de deslocamento e F é o vetor de carregamen-to. Derivando a equação em relação à variável de projeto x tem-se
(3.35)k k kμ μ μ
∂ ∂ ∂+ =
∂ ∂ ∂u K FK u (3.34)
o gradiente dos deslocamentos pode ser calculado, então, da seguinte forma:
(3.36) 1
k kμ μ μ− ⎛ ⎞∂ ∂ ∂
= −⎜k∂ ∂⎝ ⎠
u FK ⎟∂K u (3.35)
Caso as derivadas x∂
∂K e x∂
∂F sejam calculadas indiretamente, como por diferenças
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 mapea-mento 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
(3.37) ( ) ( )*
1
rr nts js r
jjk k
β
μ μ=
∂∂=
∂ ∂∑μK μ
K ; ( )* ( ) ( )r r rs r r
fϕμΓ μ F b
k k k
ϕμ μ
Ω∂ ∂ ∂= +
∂ ∂ ∂F μ μ F (3.36)
Os gradientes da matriz de rigidez e vetor carregamento são calculados somando-se a contribuição de cada região
(3.38) ( )**
1
rnrss
rk kμ μ=
∂∂=
∂ ∂∑K μK ; ( )**
1
rnrs
rk kμ μ=
∂∂=
∂ ∂∑F μF (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 relacio-nadas com as equações do problema de elasticidade, estas são estendidas para proble-mas 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 proce-dimento 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 proce-dimento 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 sensibili-dades (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 sensi-bilidades via MBR, onde as equações foram particularizadas para problemas de elastici-dade.
(a) Análise de sensibilidade através do MBR O uso das matrizes decompostas para as matrizes de rigidez e vetores de carrega-mento 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 sensibilida-de, 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)N N
k kμ μ∂ ∂
=∂ ∂u αZ (3.38)
onde a derivada dos coeficientes lineares α , pode ser obtida derivando-se diretamente a Equação (2.66),
(3.40)1
N N N
N N
k k k k kμ μ μ μ μ μ− ⎛ ⎞∂ ∂ ∂ ∂ ∂ ∂⎡ ⎤+ = ⇒ = −⎜
N
k⎣ ⎦∂ ∂ ∂ ∂ ∂⎝ ⎠
K α F α F Kα K K ⎟∂α (3.39)
ou na forma compacta (KEULEN et al, 2005)
(3.41)1 ˆN N−F
kμ∂ ⎡ ⎤= ⎣ ⎦∂α K (3.40)
onde é o pseudo vetor de força da base reduzida, definido como ˆ NF
(3.42) ˆN N
N
k kμ μ⎛ ⎞∂ ∂
= −⎜ ⎟F α ∂ ∂⎝ ⎠
F 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
(3.43)
( ) ( )
( )
*
1 1
*
1
( ) ( )
rN nr nts js Nr
s jr jk k
N r rnrs Nr Nr
b frk k k
β
μ μ
ϕ ϕμ μ μ
= =
Ω Γ
=
∂∂=
∂ ∂
∂ ⎡ ⎤∂ ∂= +⎢ ⎥∂ ∂ ∂⎣ ⎦
∑∑
∑
μK μK
F μ μ μF F
(3.42)
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 selecio-nadas 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, respectiva-mente.
O objetivo considerado é minimizar a flexibilidade. Ao volume, é imposto ser infe-rior ou igual ao seu valor inicial.
O problema de otimização pode ser formulado como:
Minimizar: ( )f μ
sujeito a
( ) 0
25mm 75mm 1,...k d
V Vk n
μμ
≤
≤ ≤ = v
Onde ( )f μ é a flexibilidade, ( )V μ o volume da estrutura, o volume inicial e
ndv o número de variáveis de projeto. 0V
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 a-mostras 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 esta-gio on-line (o tempo da avaliação é o tempo “on-line”).
45
Tabela 3.2 Resultados da análise de sensibilidade.
( )1
fμ
∂∂μ
( )2
fμ
∂∂μ
Tempo (s) daavaliação
Tempo (s) “off-line”
MEF – T3 -0.0085236 -0.015176 17.05 - MBR – T3 -0.0085230 -0.015177 0.064 75.5 MEF – Q4 -0.0085210 -0.015172 8.75 - MBR – Q4 -0.0085206 -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= 150 μ− 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úme-ro 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 imagi-nar 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 mo-delo 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 re-alizados vários processos de otimização para o problema considerado, variando o núme-ro 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 à pri-meira 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 pro-cesso 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 tama-nho 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 150 μ− e 250 μ− . As malhas utilizada foram as mesmas consi-deradas 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 X2
Flexibilidade Tempo
“on-line”(s)
Tempo
“off-line” (s)
Tempo
Total
MEF – T3 -21.5248 13.9286 0.37320 - - 67.3MBR – T3 -21.5410 13.9368 0.37317 0.6 14.7 15.3MEF – Q4 -21.5051 13.9187 0.37360 - - 118.6MBR – Q4 -21.5149 13.9236 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 significa-tivamente 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 Op-timization under Multicriteria Conditions In: V World Congress on Computational Me-chanics, Viena. 2002.
ALBUQUERQUE T. M. M. “Análise e Otimização de Problemas Térmicos e Estrutu-rais Bidimensionais Através do Método das Bases Reduzidas”. Msc. Thesis ( in Portu-guese), 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 Ver-lag, 1993.
KLEIBER M. “Parameter Sensitivity in Nonlinear Mechanics – Theory and Finite Ele-ment 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 Ênfa-se a Problemas Multiobjetivos”. dissertação de Mestrado, Universidade Federal de Per-nambuco. Recife-PE, Brasil, 2002.
MATHWORKS. “MATLAB User’s Guide”. Mathworks Inc., Natacki, 2007.
POWEL, M. J. D. “Algorithms for Nonlinear Constraints that use Lagrangian Func-tion”. 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 De-sign - with Applications”. 4rt Edition, Vanderplaats Research & Development, Inc., Colorado Springs, CO, 2004.
49
50
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.
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 algo-ritmos 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 dis-poní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 empre-gando-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 concei-to 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 imple-mentadas.
4.1 DEFINIÇÃO DO PROBLEMA
O problema multiobjetivo pode ser expresso como:
(4.1) 1 2 3min ( ) ( ), ( ), ( ), , ( ) , 2 (POM)nobjf f f f nobj⎡ ⎤= ⎣ ⎦x
F x x x x x… …≥
ne
(4.1)
Sujeito às seguintes condições:
(4.2) ( ) 0 1, ,( ) 0 1, ,
1, ,
k
i
lj j uj
h kg i nix x x j npv
= =
≤ =
≤ ≤ =
xx
………
(4.2)
onde as restrições são as mesmas definidas no capítulo 3, porém os objetivos agora for-mam 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) ( ) ( ) para 1, ,Pk kf x f x k≤ = … n
ii) ( ) ( ) para uma função objetivo ao menosPk kf x f x<
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*
x1 x2 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 bi-objetivo, os gradientes possuem direções opostas (Figura 4.2).
Figura 4.2 Pareto mínima com duas variáveis de projeto e duas funções objetivo.
x1
xP
x2
f1 x
f2
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 pon-tos ótimos de Pareto.
x2 f2
f1
F(x)
F(xP)
x
xP
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á satis-fazer, 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 basea-dos 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éto-do NBI (DAS e DENNIS, 1997) e o Método da Restrição Normal Normalizada (MES-SAC 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 substi-tui 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 me-thod” – 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, transformando-as em uma única função objetivo. Sua representação algébrica é dada da seguinte forma
(4.3) ,10 0
nobjT kj j k
k k
fFf
β=
= =∑ff
β (4.3)
onde os elementos de ,j kβ são normalizados, da seguinte maneira:
(4.4) , ,1
1 , 0 1nobj
j k j kk
β β=
= ≤ ≤∑ (4.4)
e é a função objetiva k no projeto inicial x0. k0f
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 con-torno 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 este-ja na parte não-convexa.
f1*
f1
f2*
f2
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 ante-riormente 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 da-quela do método da soma ponderada dos objetivos, sendo para isso acrescido dois pa-râ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) min
, 1,...,max min
k kk
k k
f ff k nobj
f f−
=−
= (4.5)
O novo vetor das funções objetivo é
(4.6) 1 2, ,..., ,...,k nobjf f f f⎡ ⎤⎦F = ⎣ (4.6)
Então, um coeficiente βk é associado a cada nova função objetivo e o seguinte pro-
blema é proposto:
(4.7) )(min γ (4.7) onde
(4.8) max( ), 1,...,k kf k nobjγ β= = (4.8) sujeita às restrições definidas na Equação (4.2) e às restrições adicionais
(4.9) k kfβ γ≤ k = 1,..., nobj (4.9) para vários conjuntos de vetores β, pois para cada vetor β diferente um novo subpro-blema 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 pon-tos 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 proje-tista 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 engenhari-a. 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 super-fí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 objeti-vo, 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 Pare-to. 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, denomina-do de Ponto Utópico (Shadow Minima ou Utopia Point) (DAS e DENNIS, 1996), repre-sentado por:
(4.10) ⎤⎦F * * * * *1 2 3, , , , T
nobjf f f f⎡= ⎣ … (4.10) onde cada *
if representa um mínimo local individual. Sendo o vetor xi* a solução ótima
de if , temos que * ( )i i i*f f x= . Define-se a envoltória convexa do mínimo individual
(ECMI) como:
(4.11) ⎭⎬⎫
⎩⎨⎧
≥β=βℜ∈ ∑=
0,1,: i
n
1ii
nβΦβ (4.11)
onde
(4.12) * *( , ) ( ) , 1, , ; 1, ,i j ii j f x f i nobj j 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
C
A
B 0
1f
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 encontra-dos 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.
∂ƒ ECMI
f2
f1
1( *)xF
2( *)xF
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 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,
n̂ˆt+Φβ n , com t∈ℜ, representa o conjunto
de pontos sobre , que formam uma reta quase-normal à ECMI. A interseção entre a n̂
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)
tmaxt,x (4.13)
sujeita às restrições definidas na Equação (4.2) e às restrições adicionais
(4.14) ˆ ( )t+ =n F xΦβ (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éri-cos 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 re-presentaçã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 uni-dade 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 corresponden-te. 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 con-vexa do mínimo individual (ECMI) no método NBI. Ela é dividida em m-1 segmentos, resultando em m pontos.
59
2f
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 pertencente à Linha Utópica é utilizado para definir uma reta normal (NU). Está linha normal (NU) diminui a região viável, como indicado na Figura 4.8. Como pode ser visto, minimizando-se 2f , 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.
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 Norma-lizada para a otimização multiobjetivo em geral é apresentada. O vetor x define as vari-áveis de projeto e a i-esima função objetivo. )( n
i Rf ∈f * *( xni )x R∈x é o ponto ótimo
da função if . Sendo nobj o número de funções objetivo de um problema genérico, os vetores do Plano Utópico kN são definidos como a direção dos *( )kxF até *( nobjxF ) , para k=1..nobj-1, pela equação:
(4.15) * *( ) ( ), 1,2,... 1k nobj kN x x k nobj= − = − F F (4.15)
1f
*1( )F x1
1kN
m−F0
0 1
*2( )x
2f
1f
∗1f
∗2f1
1
00
pjX
ViavelRegião
InviavelRegião
∗f
UtópicoPlano
NULinha
60
Um conjunto de np pontos distribuídos no Plano Utópico é gerado da seguinte for-ma.
(4.16) *,
1( ) , 1, 2,...
nobj
j j k k px j n k
β=
= =∑X F (4.16)
onde
(4.17) , ,1
0 1 e 1, 1, 2,...,nobj
j k j kk
jβ β=
≤ < = =∑ pn (4.17)
Agora nós geramos pontos distribuídos no espaço das funções objetivas normaliza-
das. Para cada ponto X gerados anteriormente é obtido um Ponto Pareto corresponden-te, solucionando o problema abaixo:
(4.18) Para 1, 2,..., min nobj
x
j npf
= (4.18)
sujeita às restrições definidas na Equação (4.2) e às restrições adicionais:
(4.19) ( ) 0, 1, 2,..., 1T
jkN k n − − ≤ =F X (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 nf ) 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 (di-reçã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), ilus-trada 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 minimi-zaçã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, conside-rando 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 Pa-reto 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) argu-mentam 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 tam-bé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 otimi-zação. Os valores iniciais das variáveis de projeto são μ1 = μ2 = 50mm, e os limites infe-
63
rior e superior impostos são 25mm e 75mm, respectivamente. Dois objetivos são consi-derados, 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á cons-truí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).
0 20 40 60 80 100
0
20
40
60
80
100
Figura 4.11 Projeto inicial (domínio de referência). a) malha EF (ndof = 6.732), b) ten-sã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
4000 5000 6000 7000 8000 9000 100000
0.5
1
1.5
2
2.5
3
Stra
in E
nerg
y (N
.mm
)
Volume (mm³)
FEMRBM
4000 5000 6000 7000 8000 9000 100000
0.5
1
1.5
2
2.5
3
Stra
in E
nerg
y (N
.mm
)
Volume (mm³)
FEMRBM
a)WS b)Min-Max
4000 5000 6000 7000 8000 9000 100000
0.5
1
1.5
2
2.5
3
Stra
in E
nerg
y (N
.mm
)
Volume (mm³)
FEMRBM
4000 5000 6000 7000 8000 9000 100000
0.5
1
1.5
2
2.5
3
Stra
in E
nerg
y (N
.mm
)
Volume (mm³)
FEMRBM
c)NBI 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 re-sultados totalmente obtidos via MEF. Para a obtenção de pontos de Pareto pelos proce-dimento 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 dota-do de um processador quad core, através de ferramentas do ambiente computacional MATLAB. Um ganho computacional ainda maior seria alcançado paralelizando, tam-bém, o estágio off-line.
Tabela 4.1 Placa quadrada com um orifício central – Desempenho dos algoritmos.
WS Min-Max NBI NNC
FEM 4902.0s 3846.2s 3444.3s 4613.7sRBM (37.4 s) 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, numera-das conforme indicado na Figura 4.14 (a) em função dos parâmetros variantes do pro-blema 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, co-eficiente 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 dis-tâ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μ =
. 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. [ ] [ ] ]38,28[21,105;5.0 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 con-siderados. As curvas de Pareto obtidas pelo MEF e MBR estão de acordo. Como pode-se 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 anteri-or, 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. Encon-
68
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 é apre-sentado 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 adjacen-te, 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 *(X1 ) (X3 )F F * , 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 jX (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) apresen-tam 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ópi-co 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 projeta-rem 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 con-tornos 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, en-contramos 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 to-
tal (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 procedi-mento 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, resol-ver os problemas bi-objetivos para todas as possíveis permutações de obje-tivos.
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 pon-to 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écni-cas, 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 utili-zado 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 resul-tado 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ão-Pareto, principalmente o NBI, por motivos já citados. Ambos os métodos também reali-zaram 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 dis-tribuí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) d) NBI (7553)
e) NNC (6819) 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 ( ) , 1, 2,3i ix
f x x i= =
Sujeito às seguintes restrições :
1 2 32 3 1 3 1
1 1 1 1 1 1; ;
0.2 10, 1, 2,3i
x x x2x x x x x
x i
≥ + ≥ + ≥ +
≤ ≤ =x
A Figura 4.21 apresenta os pontos de Pareto para 15 valores de , 1..3i iβ = unifor-memente 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 ( 1 2( ) x ( )f x f x ) do espaço das funções obje-tivo, onde a escala de cor se refere ao valor da terceira função objetivo 3( )f x . 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 apresenta-ram 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 igual-mente 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, atra-vé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 ilus-trativo). 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 do-minados (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 fo-ram 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 encon-trar 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 con-seguiu 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 Pa-reto, 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 Op-timization under Multicriteria Conditions In: V World Congress on Computational Me-chanics, 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 Stu-dies”, Springer, 2004.
DAS, I.; DENNIS, J.E. Normal Boundary Intersection: A New Method for Generating Pareto Surface in Nonlinear Multicriteria Optimization Problems. SIAM J. Optimiza-tion, 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 De-terministic and Probabilistic Multiobjective Optimization Algorithms with Mathemati-cal 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 Ênfa-se a Problemas Multiobjetivos”. dissertação de Mestrado, Universidade Federal de Per-nambuco. 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 Struc-tures, 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 incerte-zas. 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 ambien-te que o cerca pode não ser oportuno por várias razões, incluindo incertezas nas propri-edades dos materiais, variação na geometria devido à tolerância na fabricação, incerteza no carregamento devido à natureza não determinística do ambiente operacional, incerte-zas 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 in-certos ao realizar estudos de projeto, de fato, a consideração determinística tem sido fei-ta implicitamente em toda esta dissertação até o presente. Infelizmente, a abordagem de-terminística leva a um projeto final cujo desempenho pode cair significativamente devi-do a perturbações decorrentes de incertezas. Este problema é acentuado quando lida-se com projetos que foram severamente otimizados, pois projetos ótimos tendem a se loca-lizar 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). Co-mo conseqüência, um ótimo determinístico pode potencialmente ser uma solução de al-to 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 nis-so, 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ígi-das 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 ten-sões da forma max( )σ σ<x , 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 es-crita como max( )Sσ σ<F x
0
, onde Fs é chamado o fator de segurança. Normalmente o va-lor 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 va-lor de Fs, haja visto que o ótimo se distancia do contorno da restrição original
max( )σ σ−x = . 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 in-certezas 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 exem-plo, 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 meto-dologias 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 es-tatí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 ma-tematicamente mais rigorosa pode ser formulada a partir de três axiomas básicos.
Considerando eventos relacionados aos subconjuntos contidos num con-junto , que é o conjunto de todos os eventos possíveis, a função Prob (i.e. probabili-dade) é definida através de três axiomas:
, , ....A B CΩ
(5.1) [ ]0 1A ≤ ≤ Prob (5.1)
(5.2) [ ] 1Ω = Prob (5.2)
(5.3) [ ] [ ] [ ] [ ]A B A B A∪ = + − ∩Prob Prob Prob Prob 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) [ ][ ]
A BA B
B∩
⎡ ⎤ =⎣ ⎦Prob
ProbProb
(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 nu-mé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 ape-nas 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 u-sadas 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):
∈ξ
(5.5) b ( ) ( ) , e b
aa b P d a b aξ ξ ξ≤ ≤ = ∀ ∈ ≤∫Prob (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] = , esta função é chamada função de distribuição acumulada (“Cumulative Distribution Function” - CDF), calculada como:
( )F xξ ( )F xξ
(5.6) ∀ ∈ ( ) ( )−∞
= ∫x
F x P d xξ ξ ξ (5.6)
Se ξ é uma variável aleatória, então uma função qualquer ( )f ξ também será, po-ré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 é cha-mada de esperança da função [ ]( )E f ξ , calculado como
(5.7) [ ]( ) ( ) ( )f E f f P dξ ξ ξ ξ∞
−∞= = ∫ (5.7)
A variância ( [ ]22 ( )f fσ σ ξ= ) pode ser calculada como:
(5.8) ( ) ( )2 22 ( ) ( ) ( )f E f f f f P dσ ξ ξ ξ ξ ∞
−∞⎡ ⎤= − = −⎢ ⎥⎣ ⎦ ∫ (5.8)
ou
(5.9) ( ) ( ) ( ) ( )2 2 2−2 2 2 2( ) ( ) ( )
∞
−∞⎡ ⎤= − = − =⎣ ⎦ ∫f E f f f P df fσ ξ ξ ξ fξ (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 é cha-mada “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) ( ) ( )( ) ( ) ( ) , 2k k
k E f f f f P d kμ ξ ξ ξ ξ∞
−∞⎡ ⎤= − = − ≥ ⎢ ⎥⎣ ⎦ ∫ (5.10)
Estes momentos serão úteis aqui, para o cálculo de duas outras estatísticas: a obli-qüidade ("skewness") e a curtose (“kurtosis”). A obliqüidade é uma medida da assime-tria 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 res-pectivamente por:
(5.11) 3 43 4 3 e s μ μκ
σ σ= = − (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 (PA-POULIS, 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 des-vio padrão. Se uma variável ξ tem uma distribuição Normal, escrevesse simbolicamen-te (~ ,N )ξ μ σ , onde μ é a média e σ o desvio padrão da variável. A PDF normal pa-ra uma variável aleatória ξ é expressa por:
(5.12) 2 1/2
( ) 2P eξ μσξ πσ
−−⎛ ⎞
⎜ ⎟⎝ ⎠
⎡ ⎤⎢ ⎥ =⎢ ⎥⎣ ⎦
)
(5.12)
Sua formula padrão considera média igual a zero e desvio padrão unitário: (0,1N
(5.13) ( )2 1/2( ) 2P eξξ π
−= (5.13)
83
A Figura 5.1 apresenta o gráfico da PDF normal padrão.
-3 -2 -1 0 1 2 30
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4Função Densidade de Probabilidade
x
P(x)
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 , então exp(v ξ ) tem distribuição Lognormal, com média logμ e variância , respectivamente, iguais a: logv
(5.14) 2log
veμμ += ( ) (2log 1v vv e eμ+= )− (5.14)
Esta distribuição é muito usada para variáveis aleatórias que só aceitam valores po-sitivos, 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:
(5.15) ( ) ( )( ) 1 22 ln /ln1 1( ) exp 2
22
ab
aP b
b xb
ξξξ
ξ π
−⎛ ⎞− ⎛ ⎞= =⎜ ⎟ ⎜ aπ ξ
⎜ ⎟ ⎝ ⎠⎝ ⎠⎟ (5.15)
onde , 2ln( 1)b c= +2 1
acμ
=+
e cXσ
= é o coeficiente variacional definido apenas
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.
84
0 2 4 6 8 100
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1Pdf Lognormal com variância 1
x
P(x)
Media 1Media 2Media 6
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, po-de-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 a-presentado 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 (“Physics-
based”). 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 a-plicados no desenvolvimento de um método não-intrusivo para análise de incerteza, que explora modelos computacionais determinísticos existentes. Estas abordagens conside-ram 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 distri-buiçã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ísti-cas. 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 utili-zado para qualquer problema de propagação de incerteza (KEANE e NAIR, 2005). Da-do 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 ( ) , 1=ξ …i i m gerados aleatoriamente a partir de suas distribuições ( )P ξ , então as integrais das Equações (5.7) e (5.9) são aproxima-das, respectivamente, por:
(5.16) ( )1
1 ( )m
MCi
f fm =
= ∑ if ξ (5.16)
(5.17) [ ] ( ) ( )22 2 2 2( )
1
1ˆ( ) ( )1
m
f f i MCi
f f fm
σ σ σ=
⎡ m= = −⎢− ⎣ ⎦∑ξ ξ ⎤
⎥ (5.17)
onde m é o número de pontos amostrados, MCf é 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 MCf f→ a medida em que . O cálculo da variância pode ser usada para verificar a aproximação
m →∞
MCf 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 a-leatórias do problema. O maior problema do método de MC é sua baixa taxa de conver-gência, já que o seu erro estimado é de ordem (1O r 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
)m . Po
)iξ pode ser feito de forma inde-pendente. No caso das avaliações de ( )(f )iξ serem caras computacionalmente, recomen-da-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) ( )
1
( , )( , ) 1 miMC
i
fffm =
∂∂∂=
∂ ∂ ∂∑x ξx ξ
x x x (5.18)
86
(5.19) ( ) ( ) ( )
2 2( )
( )1
ˆ ( , )1 2 ( , ) 21
mf f i MC
i MCi
f ff m fm
σ σ
=
∂ ∂ ⎡ ∂ ⎤∂∂x
⎛ ⎞
= −⎢ ⎥⎜ ⎟∂ ∂ − ∂⎝ ⎠⎣ ⎦∑
x ξx ξ
x x x(5.19)
derivando ( )2ˆ ˆf fσ σ= , tem-se que
(5.20) ( ) ( )( ) ( )
21/2 ( )2
( )1
ˆˆ ( , )1 1ˆ ( , )ˆ2 1
mff i MC
f i MCif
f ff m fm
σσσ
σ−
=
∂∂ ⎡ ∂⎛ ⎞ ∂= = −
⎤⎢ ⎥⎜ ⎟∂ ∂ − ∂⎝ ⎠ ∂⎣ ⎦∑
x ξx ξ
x x x x (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 amos-tras 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 sem-pre 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 (a-leatórias) para cada iteração gera um erro aleatório (1 )O m dificultando (ou até impos-sibilitando) 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 otimi-zaçã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 gradien-tes) 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 fi-
xar a semente (“seed”), considere a função 2 2
2 (1 )( , )10 3x x xf x ξ ξ ξ + −
= − , onde
, e pretende-se encontrar (~ 0,1Nξ ) *x que minimize a média da função ( , )f x ξ em re-lação à ξ . A média de ( , )f x ξ pode ser analiticamente encontrada através das Equa-ções (5.7) e (5.13), obtém-se 2x=( ) 10f x , logo * 0x = . Foi realizada uma aproxima-ção da média da função de interesse por MC ( ( )MCf x ) utilizando uma amostra ξ com 30 pontos, para vários valores de x, variando-o de -1 a 1 com incrementos de 0.01. Po-ré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 ( )MCf x utilizando a mesma amostra para todos os pontos (“seed” fixo)
4. Aproxima ( )MCf 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 exem-plo, 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 ( ) ( )MCf x f x− .
Na Figura 5.3, nota-se claramente a dificuldade em otimizar a função ( )MCf x usan-do 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 ( , )if x ξ para um iξ qual-quer. Isto ocorre, pois se (MCf x) trata da média das funções ( )if x ( , )f x iξ= , i = 1..m, onde m é o tamanho da amostra. Assim, a função ( )MCf x pode ser otimizada e se encon-
trar ótimos adequados. Neste caso o ótimo de ( )MCf x será *ˆ2
x ξ20.6
MCξ ξ+
= para qual-
quer amostra (fixa) de ξ , pois
88
(5.21) 2 2
2 (1 )( )10 3MCx x xf x ξ ξ + −
= − e 2 1( , ) (3 10 )15 3
MCdf xxdx
ξ ξ ξ ξ− = + (5.21)
onde 2 e ξ ξ são respectivamente a média da amostra e a média dos quadrados da a-
mostra. À medida que a amostra aumenta (melhora) ( )22 − → 2ξξ σξ (Equação (5.9)),
com isso 2 e ξ ξ tendem, respectivamente, a zero e a um, consequentemente, *ˆMC
*x x→ . Para a amostra de 30 pontos utilizada no gráfico, -0.0275ξ = e . *ˆ -0MCx = .0469
5.3.2 Técnicas de amostragem
A geração das amostras pode ser feita de maneira totalmente aleatória (ou pseudo-aleató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 Hiper-cube 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 algorit-mos 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 algo-ritmo “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 efici-entemente 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 ge-radas 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 va-riá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 conside-rando uma distribuição Normal padrão ( )0,1N∼ . Essa equação empírica, que qualifica as amostras, faz a ponderação dos logaritmos dos erros quadráticos dos momentos cen-trais estatísticos das amostras, dando maior peso aos momentos de maior ordem. Esta função é definida como:
(5.22) ( ) ( ) ( ) ( )2 22
ln ( 1) ln lnln
2 6
2x
Err
sF
24σ κ
μ−
= + + + (5.22)
onde μ é a média da amostra, xσ o seu desvio padrão, a obliqüidade e a curtose. Para uma amostra com distribuição Normal padrão
s κ
( )10,N : μ , e s κ devem ser iguais a zero, enquanto xσ deve se aproximar de um. A amostra selecionada é a que apresenta o menor valor de . Este procedimento não pretende escolher a amostra que conse-guirá 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 ex-perimentos numéricos realizados.
ErF r
Outras duas funções “teste”, para se diferenciar as amostras, foram aplicadas: o tes-te 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 sele-cionadas 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 e X Y são variáveis aleatórias ( )1,0.1N∼ . Foi realizada
uma aproximação da média MCf 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 sele-cionada, 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 MCf 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 amos-tras 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 MCf , para os diferentes tamanhos de amostras (de 64 até 16384) pelos três tipos de amostragem. O erro foi me-dido 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 à me-dida 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 significa-tivamente. 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 resul-tados 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 à amos-tra 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 ava-liaçã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 melho-ram sua eficiência, são inviáveis para serem aplicadas diretamente em modelos comple-xos de alta fidelidade. O Método da Colocação Probabilística (“Probabilistic Collocati-on Method” - PCM) (RAMAMURTHY, 2005) é uma ferramenta desenvolvida visando uma análise de incerteza eficiente em modelos complexos e computacionalmente custo-sos. 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 n-1, 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 es-parsas (“sparse grids) (HEISS e WINSCHEL, 2008).
O PCM se baseia nos conceitos de Quadratura de Gauss e de polinômios ortonor-mais. A idéia básica do método é aproximar a função de interesse por funções polinomi-ais 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 supraci-tados.
(a) Polinômios ortogonais
Polinômios são ditos ortogonais entre si com relação a um produto interno relacio-nado a um espaço, se este for nulo. Todo produto interno em um dado espaço , de funções reais f, g e h pertencentes a F , deve satisfazer as seguintes condições (APOS-TOL, 1967):
F
(5.23)
, , ,
, , , , onde é um escalar
, ,
, 0, se 0
f g h f g f h
f g f g f g
f g g f
f f f
α α α α
+ = +
= =
=
> ≠
(5.23)
Dado um espaço linear de funções reais e considerando duas funções polinomi-ais , o produto interno aqui tratado é definido como:
F( ), ( )f x g x ∈F
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 pertencentes ao espaço polinomial H , pode ser definido como:
( )ih x
(5.25) j 2,0 ,1 ,2 , ,
0
( ) ...i
ii i i i i i i j
j
h x a a x a x a x a x=
= + + + =∑ (5.25)
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 , i=0,1..n, onde o índice de indica o grau do polinômio:
( )ih x ih
(5.26) 1, para
( ), ( )0, para i j
i jh x h x
i j≠
=⎧= ⎨⎩
(5.26)
Através destas relações são encontrados os coeficientes que definem os poli-nômios ortonormais. Estes polinômios são únicos para cada função de ponderação dada, e formam uma base para . Todas as raízes
,i ja
H * , 1..jx j =* * , jx j= ∈F
i
,
).
de um polinômio , estão
contidas no espaço real, ou seja 1..i , e dependem apenas da função de ponderação P(x). As raízes d ( )x formam os pontos de colocação da qua-dratura de Gauss. Para mais detalhes consulte (GAUTSCHI, 2005
( )ih x
( ) 0,i jh x =
e ih
Note que é uma constante, ( )0h x ( )0h x a= 0,0 , consequentemente propriedades im-portantes, a partir da Equação (5.26), podem ser obtidas, as quais serão usadas adiante.
(5.27) ( ) ( ) ( )
( ) ( )
0
20 0 0,0
, 0, ( ) 0,
, 1, ( ) 1
i iF
F
h x h x h x P x dx i
h x h x a P x dx
0= = >
= =
∫∫
(5.27)
(b) Quadratura de Gauss
Na integração numérica via quadratura de Gauss para integrais da seguinte forma:
93
(5.28) ( ) ( )F
f x P x dx∫ (5.28)
Aproxima-se a função ( )f x , por um polinômio de grau 2n-1 , a partir da base orto-normal do espaço H , em relação à função de ponderação , onde n é o número de pontos de integração, tal como segue
( )P x
(5.29) 1 1
0 0
ˆ( ) ( ) ( ) ( ) ( )n n
i i n i ii i
f x f x b h x h x c h x− −
= =
⎛ ⎞ ⎛= +⎜ ⎟ ⎜⎝ ⎠ ⎝ ⎠∑ ∑ ⎞
⎟ (5.29)
Na integral da Equação (5.28) aproximada por quadratura de Gauss, o segundo ter-mo da aproximação se cancela (por ortogonalidade) e lembrando as propriedades (5.27) os termos para do primeiro somatório da aproximação também é cancelado. A integral desejada
1.. 1i n= −(5.28), aproximada por quadratura de Gauss, pode ser expressa en-
tão da seguinte forma:
(5.30) ( ) ( ) ( )0 0F Ff x P x dx b h P x dx∫ ∫ (5.30)
Para encontrar os coeficientes bi e ci da aproximação (Equação (5.29)), seria neces-sá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 (nh ncelando as-sim o segundo termo da aproximação apresentada na Equação
de ca
0
1n
)x ,(5.29), pois
. Os coeficientes que definem os polinômios ortonormais são calculados através das relações
*( ) 0, 1..n ih x i n= = ,i ja(5.26), como já mencionado. Desta forma, os coeficien-
tes bi são encontrados resolvendo o sistema:
(5.31)
* *1 0 1 11
* *
0 * *0 1
( ) ( )( ) ( ) ; ou
( ) ( )
nn
i j j ij
n n n
f x h h x bf x b h x
f x h h x b
−−
=− −
⎡ ⎤ ⎡ ⎤ ⎡ ⎤⎛ ⎞ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥= =⎜ ⎟ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥⎝ ⎠ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥⎣ ⎦⎣ ⎦ ⎣ ⎦
∑ (5.31)
Definindo-se a matriz inversa, como *( )j ih x
(5.32) ⎥⎥
*0 1 1
*0 1
1( )
( )
n
n n
h h xAp
h h x
−
−
−⎡ ⎤⎢= ⎢⎢ ⎥⎣ ⎦
(5.32)
os coeficientes b0 podem ser calculados como:
94
(5.33) *i
*0 1
*0 1,
1*1
( )( ) , logo ( )( )
n
i i ii
n n
b f xb Ap f x b Ap f x
b f x =−
⎡ ⎤⎡ ⎤⎢ ⎥⎢ ⎥ = =⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎥⎣ ⎦ ⎣ ⎦
∑ (5.33)
Para o cálculo da integral, basta conhecer a primeira linha da matriz e ponderar
as respostas da função nos pontos de integração (raízes de ), então, o vetor dos pesos P é definido como .
Ap
nh
1,i iP Ap=
Definindo-se
(5.34) x 0 0 ( )F
C h P x d= ∫ (5.34)
o valor da integral (5.30) é aproximado por . 0 0b C
Os passos para a determinação dos parâmetros necessários a aplicação da metodo-logia são:
1. A partir de uma função de ponderação qualquer ( )P x (não negativa), encon-tra-se os polinômios ( )ih x , i=1..n (cujos coeficientes ,i ja são obtidos utili-zando (5.26)).
2. Calcula-se as raízes *ix de ( )nh x .
3. Calcula-se a matriz Ap .
4. Computa-se o vetor de ponderação P, onde 1,i iP Ap= .
5. Calcula-se o valor da integral utilizada na Equação (5.30) e defini-se uma constante 0C Equação (5.34).
Restando apenas avaliar a resposta *( )if x em cada ponto de integração, para então calcular o valor desejado:
(5.35) *0
1( ) ( ) ( )
n
i iFi
f x P x dx C P f x=∑∫ (5.35)
Vale ressaltar que não é necessário calcular estes parâmetros para as funções de dis-
tribuição conhecidas, pois estes valores são conhecidos e podem ser encontrados na lite-ratura, 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 ( ) 1F
P d =∫ ξ ξ e pela segunda propriedade em (5.27)
, consequentemente a constante 0 1h = 0 1C = (definida pela Equação (5.34)). Com isso o valor da integral (5.30) é aproximado apenas por . 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:
0b
*( )
1
2 * 2( )
1
( )
ˆ ( )
n
PC i ii
n(5.36)
2PC i i PC
i
f P f
P f fσ
=
=
=
= −
∑
∑
ξ
ξ
n
(5.36)
onde são as raízes do polinômio ortonormal. *( ) , 1..i i =ξ
Pode-se calcular a partir da Equação (5.36) os gradientes da média e do desvio pa-drã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 proje-to x
(5.37) ( )( , )i
∂
x ξ
1
( , ) mPC
ii
fff P=
∂∂∂=
∂ ∂ ∑x ξx x x
(5.37)
(5.38) ( ) ( )2 2
( )( )
( ,, ) i
1
ˆ )2 (
mf PC 2i i
f PCPC
i
fP fσ σ
=
∂ ∂ ⎛ ⎞f
∂= ⎜
x ξx ξ ∂
−⎟∂ ∂ ∂⎝ ⎠∑x x x ∂x
(5.38)
derivando ( )2ˆ ˆPC PCσ σ= , tem-se que
(5.39) ( ) ( )21/2 ( )2
( )1
1 1 ( ,ˆ
m
iPC
ˆ ( , )ˆ ˆ2
PC iPCPC i )
f PCPC
ff fσσ σ
− ∂
σ =
⎡ ∂ ⎤⎛ ⎞∂ ∂= = −⎢ ⎥∂⎜ ⎟∂ ∂ ∂⎝ ⎠⎣ ⎦
x ξξ
x x x∑ xx
(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 pon-tos 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 multi-variá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 , raízes do polinômio *( ) , 1..i iξ = n
( )nh ξ e os coeficientes de ponderação , porém estes valores apenas devem ser calcu-lados 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
iP
( )0,1N∼ , pois os mesmos valores podem ser aplicados a uma distribui-ção normal com parâmetros diferentes, bem como a uma distribuição Lognormal, atra-vé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 dos polinômios ortonormais
, , 0..i ka k = in( ), 0..ih iξ = , 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, , através das relações descritas na Equação
0,1i =(5.26) e sabendo que para uma distribuição
normal padrão , ( )P dξ ξ ξ∞
−∞=∫ 0 2 ( )P dξ ξ
∞
−∞1ξ =∫ e ξ ∈ , obtém-se:
(5.40) 0 0
20,0
( ), ( ) 1
( ) 1
h h
a P d
ξ ξ
ξ ξ∞
−∞
=
=∫ (5.40)
logo, . E para 0,0 1a =
(5.41) ( )( )0 1
0,0 1,0 1,1
1,0 1,1
( ), ( ) 0
( ) 0
( ) ( ) 0
h h
a a a P d
a P d a P d
ξ ξ
ξ ξ ξ
ξ ξ ξ ξ ξ
∞
−∞
∞ ∞
−∞ −∞
=
+ =
+ =
∫∫ ∫
(5.41)
logo, . E ainda, para 1,0 0a =
(5.42) ( )( )1 1
1,0 1,1 1,0 1,1
2 21,0 1,1 1,0 1,1
( ), ( ) 1
( ) 1
( ) 2 ( ) ( ) 1
h h
a a a a P d
a P d a a P d a P d
ξ ξ
ξ ξ ξ ξ
ξ ξ ξ ξ ξ ξ ξ ξ
∞
−∞
∞ ∞ ∞
−∞ −∞ −∞
=
+ + =
2+ + =
∫∫ ∫ ∫
(5.42)
logo, . 1,1 1a =
97
Portanto, os dois primeiros polinômios ortonormais são: 0 1( ) 1 e ( )h hξ ξ ξ= = . Para se encontrar um polinômio de grau j, a partir dos polinômios de menor grau já calcula-dos (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) , ,0 0
0,( ), ( ) ( ) =
1,
jik k
i j i k j kk k
i jh h a a P d
i jξ ξ ξ ξ ξ ξ
∞
−∞= =
≠⎛ ⎞ ⎧⎛ ⎞= ⎨⎜ ⎟⎜ ⎟ =⎝ ⎠⎝ ⎠ ⎩∑ ∑∫ (5.43)
para . Onde todos os coeficientes dos polinômios 0i = … j ,i ka ( )ih ξ , para , já são conhecido. O sistema é resolvido considerando os coeficientes do polinômio
0 1i j= −…
,j ka(jh )ξ , como incógnitas, onde . Esse procedimento é repetido recursivamente
variando o j até atingir o polinômio de grau n (0k = … j
h ( )n ξ ) 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 *, 1..j j nξ = de ( )nh ξ ,
2. Cria-se a matriz de *( )i jh ξ avaliando os polinômios ( ), 0 1ih i nξ = −… em
cada *, 1..j j nξ = .
3. Calcula-se sua inversa: Ap
4. Extraí-se da 1ª linha de Ap os pesos 1,i iP Ap= .
Os valores de e jP *, 1..j j nξ = são armazenados e utilizados para obter as estatísti-cas 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 , a uma distribuição normal com parâmetros diferentes (0,1N∼ ) )( ,N μ σ∼ , basta rea-
lizar 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
(5.46) 2
2 2ln μη
μ σ
⎛ ⎞⎜ ⎟=⎜ ⎟+⎝ ⎠
e 2
2ln 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 de cada polinômio ortonormal. Esta função sim-
bolicamente calcula as integrais de
,i jab
a( )i P dξ ξ ξ∫ para 0...2 1i n= − , onde a função de
ponderação ( )P ξ , os limites [a,b] onde ela é definida e o valor de n são dados de entra-da. A partir dos valores das integrais obtidos, cria-se um sistema de equações não linea-res, apresentado na Equação (5.43), o qual também é resolvido simbolicamente, para se obter os coeficientes de cada polinômio ortonormal. ,i ja
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:
,i ja =
[ 1] [ 0, 1] [ -1, 0, 1]/2*2^(1/2) [ 0, -3, 0, 1]/6*6^(1/2) [ 3, 0, -6, 0, 1]/12*6^(1/2) [ 0, 15, 0, -10, 0, 1]/60*30^(1/2)
[ -15, 0, 45, 0, -15, 0, 1]/60*5^(1/2) [ 0,-105, 0, 105, 0, -21, 0, 1]/420*35^(1/2)
[ 105, 0,-420, 0, 210, 0, -28, 0, 1]/1680*70^(1/2) [ 0, 945, 0,-1260, 0, 378, 0, -36, 0, 1]/5040*70^(1/2)
para . 0 , e 0i n j= =… …i
Estes polinômios ortonormais estão de acordo com os polinômios ortogonais encon-trados 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 ( 9 ( )h ξ ), 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 refe-rência (GREENWOOD e MILLER, 1948), onde são obtidas as raízes de polinômios or-togonais 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. Es-tes 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 *( )i jh ξ avalian-
do os polinômios ( ), 0 1ih i nξ = … − em cada *, 1..j j nξ = . Só restando então calcular Ap, invertendo a matriz criada, e extrair da 1ª linha os pesos P, onde i iPobtidos para os nove pontos de colocação foram:
1,Ap= . Os pesos
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 refe-rê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
e P *ξ 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 8PN = de uma variável aleatória (0,1)Nξ ∼ , 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 é:
(5.47) ( ) ( )( )( ) (
( )
1
1
1 1( ) ( 1) , ou ( ) 1 2 3 ... 8
onde, ! ( 1) 1( 2)3...7( 8)
P
P
Ni
i
Ni
Pi
f i fC C
C N
ξ ξ ξ ξ ξ ξ ξ=
=
= − − = + − +
= − − = − −
∏
∏
)− (5.47)
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 2( )f ξ , 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 0.000319231598934
8 15 0.000000000000000 0.000022214207163
9 17 0.000000000000000 0
101
Como comparativo, a Tabela 5.2 ilustra a aproximação do mesmo caso por MC uti-lizando 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 e X Y ( )1,N σ∼ são variáveis aleatórias in-dependentes. Para este exemplo o valor do desvio padrão σ será variado para que se te-nham 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 in-tervalo (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 (b) desvio padrão
(c) erro da média (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 des-
vio padrão.
103
Note que, quanto menor o desvio padrão das variáveis aleatórias, mais “simples” se-rá 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 pa-drã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 peque-na diferença. Só a partir deste ponto, com o aumento do desvio padrão da variável alea-tória ( 1.6σ > ), o MC com amostra LHS de 50.000 pontos começa a apresentar resulta-dos 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 conti-
nuamente diferenciável, confrontando-os com a resposta exata. Considere a função 2( , ) (1 )f x x x xξ ξ ξ= − − + − , onde ξ ( )0,1N∼ é 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 ( )MCf x por MC utilizando uma amostra a-leató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 simbo-licamente 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 ana-
liticamente, por MC e por PCM para diferentes graus de aproximação.
Como para cada valor iξ da amostra, ( , )if x ξ apresenta uma singularidade em , ( 1.. )i ix i nξ= = , então a função ( )MCf x terá n pontos de singularidade (onde n é o ta-
manho da amostras). Isto resulta em singularidades do gráfico ( )MCf 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 ( )MCf 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 par-ticular da função do exemplo anterior, onde x=0: (0, ) ( )f fξ ξ ξ ξ= = − onde
. Utilizando o PCM, aproximou-se a função (1,0Nξ ∼ ) ( )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 , não necessita ser calculada, pois para qualquer valor de c a inte-
gral desejada (Equação
*( ) 0nh ξ =(5.30) não tem seu valor alterado. Pode-se então definir qual-
quer 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 (E-quaçã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 a-proximaçã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 o-timizaçã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 deter-miní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 prejudi-ciais ao projeto. A Otimização Robusta leva em consideração as variáveis incertas (alea-tó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 ξ , são as variáveis de projeto e n∈x ∈ξ E as variáveis incertas, onde E é o
espaço amostral das variáveis. Um objetivo robusto e utópico seria encontrar um , de tal forma que para todo
*x∈ξ E
(5.48) *( , ) ( , ) nf f≤ ∀ ∈x ξ x ξ x (5.48)
109
Achar tal que simultaneamente minimize *x ( , )f x ξ , para todos os , é o pro-blema 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 proble-ma, varias abordagens são usadas neste contexto. A primeira visa encontrar um que minimiza o pior caso ( mais desfavorável para
∈ξ E
*x∈ξ E (f x, )ξ ), ou seja,
(5.49) ( )
[ ]
( ) max ( , )
min ( )n
F f
F∈
∈
=ξ
x
x x ξ
xE (5.49)
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 o-timização é realizada sobre o pior caso dentro da região ( ,A )εξ , neste caso
(5.50) ( )( , )
( , ) max ( , )A
F fε
ε∈
=ξ ξ
x x ξ (5.50)
onde ε é o parâmetro regularizador que define o tamanho da região, tal que:
(5.51) ( ) 00lim ( , ) ( , )F fε
ε→
→x x ξ (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 ( ), como sendo o valor de f onde há uma probabilida-de
( , )F wx( )p w de f apresentar um valor menor, assim
(5.52) ( , ) ( ) ( )fF w f wσ= +x x 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 . O parâmetro w controla o conservadorismo da so-luçã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 (
∈ξ E
( ) 50%p w = ) apenas a média seria o-
110
timizada, não sendo considerada esta, uma abordagem de pior caso, pois se teria o pior caso com uma probabilidade de 50%. Quando tende-se a um pior caso com uma probabilidade de 100% e apenas o desvio padrão passa a ser considerado, assim
w →∞
(5.53) ⎤⎦x min lim ( , ) min ( )n fw
F w σ→∞∈ ∈
⎡ ⎤ ⎡→ ⎣⎣ ⎦x xx
n(5.53)
Neste caso pode ser considerada a otimização do desvio padrão da função. Assim a so-lução será o ponto onde há a menor variação da resposta de interesse, esta constitui a i-déia central da otimização robusta.
A segunda abordagem, motivada pelo princípio de Bayes, tenta encontrar , tal que
*x
(5.54) *( , ) ( , ) nf f≤ ∀ ∈x ξ x ξ x (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 refe-rê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 me-nos 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 otimi-zação é realizado sobre o desvio padrão se está sendo conservativo, e encontra-se o pon-to onde há a menor variação da função de interesse, sendo esta uma das principais me-didas 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 (o-timização multiobjetivo robusta - OMR).
O problema de OMR mencionado anteriormente pode ser formulado como:
Minimizar: ( ) ( )( ) , ( )E f fσ⎡ ⎤⎣ ⎦x x
( )( )
0 1,...
0 1,...
1,...
i
j
uk k k d
g i
h j
v
m
x x x k n
≤ =
= =
≤ ≤ =
x
x
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 o vetor das variáveis de projeto. x
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 conside-rando 15 pontos de Pareto.
O problema de otimização pode ser formulado como:
Minimizar: ( ) ( )3 3( ,E ) , ( ,E )E f fμ σ μ⎡ ⎤⎣ ⎦
sujeito à
( )( ) ( )( )( )
( ) 3 ( ) 3
0
,E 3 ,E 7 MPa 1,...eq i eq i elE i
V V
τ μ σ τ μ
μ
+ ≤
≤
n=
v
25mm 75mm 1,...k dk nμ≤ ≤ =
Onde 3( ,E )f μ é a flexibilidade, ( )( ) 3, Eeq iτ μ é a tensão Von Misses no elemento i,
( )V μ o volume da estrutura, o volume inicial, nel é o número de elementos e ndv o número de variáveis de projeto.
0V
O modelo de elementos finitos considerado possui 3900 graus de liberdade com e-lementos 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á men-cionados o MC e o PCM. Para definir o número de pontos (tamanho da amostras) utili-
112
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ál-culo das estatísticas da flexibilidade, foi repetido 100 vezes, para cada tamanho de a-mostra. 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 des-vio 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üên-cia. 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ál-culo 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
da amostra Média Desvio Padrão
202
403
806
1613
3225
6451
12902
0.308530959713673
0.308550213518812
0.308556674596121
0.308558482597648
0.308559106575475
0.308559457041847
0.308559963766802
0.035326555620201
0.035410907482470
0.035455905892811
0.035455246318295
0.035453460837231
0.035451940103915
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 a-centuada 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
da amostra Média Desvio Padrão
2
3
4
5
6
7
8
0.308541286854411
0.308559984101546
0.308560062653515
0.308560062770421
0.308560062772797
0.308560062777446
0.308560062777447
0.034854261923426
0.035445314900078
0.035455152887810
0.035455264870394
0.035455265822598
0.035455265989171
0.035455265989205
(b) Resultado da otimização
A Figura 5.14 apresenta as distribuições dos pontos Paretos obtidos para os méto-dos 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 conse-guem 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) Min-Max, 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 re-sultados obtidos via MC.
Tabela 5.5 Placa quadrada com um orifício central – desempenhos dos algoritmos.
(segundos) WS Min-Max NBI NNC
MC 5000 pontos 38 623 26 486 19 480 18 633
PC 5 pontos 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 consider-ing uncertainties – An overview. Computational Methods and Applications in Mechani-cal 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. Medei-ros, 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 Pur-suit 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”. Inter-national Conference on Safety and Reliability – ESREL ´99, Munich, Alemanha, 1999.
116
117
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 optimi-zation using cascade evolutionary computation”, Comput. Methods Appl. Mech. Engrg. 194 (2005) 3496–3515.
FOO, J., WAN, X., KARNIADAKIS, G. E. “The multi-element probabilistic colloca-tion 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 Sam-pling”. 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 un-certainties in black box systems”. Thesis (M.S. in computer science). Washington State University. 2005.
118
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 oti-mização apresentadas, foram desenvolvidos, implementados e verificados. Foram con-siderados 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 fren-te 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 desenvol-vidas 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 imple-mentado) 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 diver-sas formas de mapeamento consideradas foi automatizado.
Foram implementados os métodos para a obtenção de pontos de Pareto para pro-blemas OM para a solução de problemas 2D contínuos. Os esquemas considerados fo-ram: 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 apre-sentados 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 envol-vendo análise estrutural e/ou térmica.
119
• 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 oti-mizaçã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 dife-rentes 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 com-portamento 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 es-tatí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 interes-se, 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 indi-cam claramente o mais eficaz. Os testes de convergência indicam que a imple-mentaçã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 amostra-dos, em relação ao número de variáveis de projeto. Devido à própria metodolo-gia 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 a-proximado.
• 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 obje-tivo, obtidos através das quatro técnicas, bem como pela modificação proposta do NBI, que mostrou ser eficaz para os exemplos estudados.
120
• 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 em-pregadas 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 pro-blemas 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 e-ficientes para lidar com problemas multiobjetivo, permitiram a obtenção das so-luções OMR em pouco tempo computacional. Em alguns casos a eficiência al-canç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, apre-senta-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 proje-to 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 geo-métricas do problema.
• Adaptar o MBR para otimização de forma.
121
• 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 di-minuir 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 ob-jetivo, 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 incor-porada 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 es-tatísticas.
• Automatizar, por meio de análise de erro, o cálculo do tamanho das amostras pa-ra 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 Pare-to, necessários para a definição da fronteira de Pareto.
