Embed Size (px)
Citation preview
UM ALGORITMO PARA OTIMIZAÇÃO RESTRITA COM APROXIMAÇÃODE DERIVADAS
Henry Octavio Cortés Ramos
Tese de Doutorado apresentada ao Programade Pós-graduação em Engenharia Mecânica,COPPE, da Universidade Federal do Riode Janeiro, como parte dos requisitosnecessários à obtenção do título de Doutorem Engenharia Mecânica.
Orientadores: José Herskovits NormanAurélio Lima Araújo
Rio de JaneiroNovembro de 2011
UM ALGORITMO PARA OTIMIZAÇÃO RESTRITA COM APROXIMAÇÃODE DERIVADAS
Henry Octavio Cortés Ramos
TESE SUBMETIDA AO CORPO DOCENTE DO INSTITUTO ALBERTO LUIZCOIMBRA DE PÓS-GRADUAÇÃO E PESQUISA DE ENGENHARIA (COPPE)DA UNIVERSIDADE FEDERAL DO RIO DE JANEIRO COMO PARTE DOSREQUISITOS NECESSÁRIOS PARA A OBTENÇÃO DO GRAU DE DOUTOREM CIÊNCIAS EM ENGENHARIA MECÂNICA.
Examinada por:
Prof. José Herskovits Norman, D.Ing.
Prof. Marcelo José Colaço, D.Sc.
Prof. Anatoli Leontiev, D.Sc.
Prof. Cristóvão Manuel Mota Soares, D.Sc.
Prof. Antonio André Novotny, D.Sc.
RIO DE JANEIRO, RJ – BRASILNOVEMBRO DE 2011
Ramos, Henry Octavio CortésUm algoritmo para otimização restrita com aproximação
de derivadas/Henry Octavio Cortés Ramos. – Rio deJaneiro: UFRJ/COPPE, 2011.
XII, 118 p.: il.; 29, 7cm.Orientadores: José Herskovits Norman
Aurélio Lima AraújoTese (doutorado) – UFRJ/COPPE/Programa de
Engenharia Mecânica, 2011.Referências Bibliográficas: p. 111 – 118.1. otimização não linear. 2. aproximação das
derivadas. 3. modelos aproximados. I. Norman,José Herskovits et al. II. Universidade Federal do Rio deJaneiro, COPPE, Programa de Engenharia Mecânica. III.Título.
iii
A Deus. A minha filha Isabela.A meus pais, irmãos e a minha
esposa Andreia.
iv
Agradecimentos
Aos professores José Herskovits Norman e Aurélio Lima Araújo pela orientação,amizade, apoio e conhecimentos transmitidos ao longo da realização deste trabalho.
Ao Prof. Cristóvão Manuel Mota Soares (Universidade Técnica de Lisboa), porsua cooperação constante.
Ao Prof. Jean R. Roche (Nancy Université, França) pela amizade.Aos docentes e pessoal administrativo, do Programa de Engenharia Mecânica
PEM-COOPE/UFRJ, sempre dispostos a ajudar e por proporcionar um excelenteambiente de convívio.
A Andreia, e a meus pais: Carlos e Myriam, pelo carinho, incentivo e apoioconstantes.
A todos os colegas e amigos do Laboratório e da universidade: Eduardo SousaMota, Alfredo Canelas, Mario Tanaka, Gabriel Guerra, Miguel Aroztegui, Elmer,Jorge Zerpa, Pavel e Helena, Sandro e aos outros com os quais trabalhei junto.
Agradeço ao Conselho Nacional de Desenvolvimento Científico e Tecnológico(CNPq), a Fundação Carlos Chagas Filho de Amparo à Pesquisa do Estado do Riode Janeiro (FAPERJ), e a Coordenação de Aperfeiçoamento de Pessoal de NívelSuperior (CAPES) pelo suporte financeiro.
v
Resumo da Tese apresentada à COPPE/UFRJ como parte dos requisitos necessáriospara a obtenção do grau de Doutor em Ciências (D.Sc.)
UM ALGORITMO PARA OTIMIZAÇÃO RESTRITA COM APROXIMAÇÃODE DERIVADAS
Henry Octavio Cortés Ramos
Novembro/2011
Orientadores: José Herskovits NormanAurélio Lima Araújo
Programa: Engenharia Mecânica
Problemas de otimização em que não se dispõe das derivadas são cada vez maisfrequentes na prática da engenharia. As duas principais estratégias para trabalharcom estes problemas consistem no uso de algoritmos sem derivadas ou na adequaçãode algoritmos convencionais baseados no gradiente através da aproximação dasderivadas. O baixo desempenho dos algoritmos sem derivadas e sua dificuldade emrespeitar as restrições de projeto desestimula seu uso na comunidade científica. Poroutra parte, a técnica mais usada para capacitar algoritmos convencionais baseadosno gradiente na solução deste tipo de problemas é sem dúvida o uso de diferençasfinitas, mas ainda esta técnica apresenta um alto custo computacional e em algunscasos não pode ser aplicada, como acontece na presença de ruido externo. Comisso em vista, neste trabalho é apresentado um método para otimização restritaque é baseado num algoritmo de direções viáveis (FDIPA) e na aproximação dasderivadas. O mesmo, aqui chamado FDIPA-AG, permite obter soluções viáveismantendo a eficiência dos algoritmos baseados no gradiente. Sua variante FDIPA-AG-N é apresentada para o caso onde há ruido moderado.
As técnicas propostas são usadas na solução de problemas de otimização deestruturas reticuladas e são comparadas com dois dos algoritmos mais conhecidos.Duas aplicações na área de estruturas feitas com compósitos laminados sãoabordadas, consistindo na identificação de propriedades mecânicas e na distribuiçãoótima de atuadores piezoelétricos.
vi
Abstract of Thesis presented to COPPE/UFRJ as a partial fulfillment of therequirements for the degree of Doctor of Science (D.Sc.)
AN ALGORITHM FOR CONSTRAINED OPTIMIZATION WITHDERIVATIVES APPROXIMATION
Henry Octavio Cortés Ramos
November/2011
Advisors: José Herskovits NormanAurélio Lima Araújo
Department: Mechanical Engineering
Constrained optimization problems with no derivatives available are increasinglycommon in engineering practice. The two main strategies for working with suchproblems are the use of derivative free algorithms or the adaptation of gradientbased conventional algorithms by derivatives approximation. The poor performanceof derivatives free algorithms and their difficulty to respect the design constraintsdiscourages their use in the scientific community. On the other hand, the more usedtechnique to enable conventional algorithms based on the gradient in solving thiskind of problems is certainly the use of finite differences, but this technique has a highcomputational cost and in some cases can not be applied, as happen in the presenceof external noise. With this in mind, this work presents a method for constrainedoptimization that is based on a feasible directions algorithm (FDIPA) and derivativesapproximation. So, here called FDIPA-AG, allows get feasible solutions keeping theperformance of gradient based algorithms. Its variant FDIPA-AG-N is presentedfor the case where there is mild noise.
The proposed techniques are used to solve truss structures optimization problemsand are compared with two of the most well known algorithms. Two applications instructures made from laminated composite are addressed, consisting of mechanicalproperties identification and optimal distribution of piezoelectric actuators.
vii
Sumário
Lista de Figuras x
Lista de Tabelas xii
1 Introdução 1
2 Otimização sem derivadas 32.1 Limitações da otimização sem derivadas . . . . . . . . . . . . . . . . 42.2 Classificação das técnicas de otimização sem derivadas . . . . . . . . 5
3 Modelos aproximados e amostragem 83.1 Modelos lineares . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
3.1.1 Modelos lineares de interpolação e regressão . . . . . . . . . . 83.1.2 Limites do erro para interpolação e regressão linear . . . . . . 11
3.2 Modelos não lineares . . . . . . . . . . . . . . . . . . . . . . . . . . . 153.2.1 Interpolação de modelos não lineares . . . . . . . . . . . . . . 163.2.2 Regressão de modelos não lineares . . . . . . . . . . . . . . . . 21
3.3 Modelos não polinomiais . . . . . . . . . . . . . . . . . . . . . . . . . 253.3.1 Classificação dos modelos aproximados . . . . . . . . . . . . . 253.3.2 Funções de base radial . . . . . . . . . . . . . . . . . . . . . . 273.3.3 Modelos de Kriging . . . . . . . . . . . . . . . . . . . . . . . . 29
3.4 Amostragem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 303.4.1 Projeto de experimentos . . . . . . . . . . . . . . . . . . . . . 31
4 Estimação das derivadas e controle do erro 384.1 Cálculo das derivadas aproximadas . . . . . . . . . . . . . . . . . . . 38
4.1.1 Estimativa das derivadas usando interpolação linear . . . . . . 384.1.2 Estimativa das derivadas usando regressão linear . . . . . . . . 39
4.2 Controle do erro na aproximação das derivadas . . . . . . . . . . . . . 42
viii
5 Algoritmo de pontos interiores e direções viáveis com aproximaçãodas derivadas 505.1 Conceitos preliminares . . . . . . . . . . . . . . . . . . . . . . . . . . 50
5.1.1 Problema de otimização não linear . . . . . . . . . . . . . . . 505.1.2 Definições . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 515.1.3 Condições de otimalidade . . . . . . . . . . . . . . . . . . . . 53
5.2 Algoritmo FDIPA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 545.3 Algoritmo FDIPA-WLS-AG . . . . . . . . . . . . . . . . . . . . . . . 56
5.3.1 Descrição do algoritmo . . . . . . . . . . . . . . . . . . . . . . 575.4 Algoritmo FDIPA-AG . . . . . . . . . . . . . . . . . . . . . . . . . . 61
5.4.1 Descrição do algoritmo . . . . . . . . . . . . . . . . . . . . . . 62
6 Resultados numéricos 666.1 Exemplo 1 - Problema Algébrico de Barnes . . . . . . . . . . . . . . . 666.2 Exemplo 2 - Treliça de 10 barras . . . . . . . . . . . . . . . . . . . . 706.3 Exemplo 3 - Domo de 52 barras . . . . . . . . . . . . . . . . . . . . . 72
7 Tolerância ao ruído externo 767.1 Tamanho mínimo admissível para a região de interpolação . . . . . . 777.2 Direção de busca conservativa . . . . . . . . . . . . . . . . . . . . . . 807.3 Algoritmo FDIPA-AG-N . . . . . . . . . . . . . . . . . . . . . . . . . 82
7.3.1 Descrição do algoritmo . . . . . . . . . . . . . . . . . . . . . . 827.4 Resultados numéricos . . . . . . . . . . . . . . . . . . . . . . . . . . . 84
7.4.1 Treliça de 10 barras . . . . . . . . . . . . . . . . . . . . . . . . 857.4.2 Domo de 52 barras . . . . . . . . . . . . . . . . . . . . . . . . 87
8 Otimização de estruturas feitas com compósitos laminados 898.1 Materiais compósitos laminados . . . . . . . . . . . . . . . . . . . . . 90
8.1.1 Classificação dos materiais compósitos . . . . . . . . . . . . . 908.2 Identificação e otimização de estruturas compósitas laminadas . . . . 92
8.2.1 Modelo numérico da placa . . . . . . . . . . . . . . . . . . . . 938.3 Identificação das propriedades elásticas do laminado . . . . . . . . . . 988.4 Distribuição ótima dos arranjos sensor-atuador piezoelétricos . . . . . 102
9 Considerações finais 1099.1 Conclusões . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1099.2 Sugestões para trabalhos futuros . . . . . . . . . . . . . . . . . . . . . 110
Referências Bibliográficas 111
ix
Lista de Figuras
3.1 Classificação geral dos modelos de aproximação . . . . . . . . . . . . 263.2 Projeto 33 fatorial completo (27 pontos). . . . . . . . . . . . . . . . . 333.3 Projeto composto central de três variáveis (15 pontos). . . . . . . . . 343.4 Projeto D-ótimo de três variáveis (11 pontos). . . . . . . . . . . . . . 353.5 Um projeto LHS com n = 2, ne = 5 para x = (x1, x2) uniformemente
distribuída num quadrado unitário. . . . . . . . . . . . . . . . . . . . 37
4.1 Função de ponderação wi para WLS. . . . . . . . . . . . . . . . . . . 414.2 Representação da matriz M0 dada pela Eq. (4.14). . . . . . . . . . . . 444.3 Curvas de nível para cond(M0) (preto) e para
∥∥∥(M0)−1∥∥∥ (cinza). . . . 45
4.4 Iso-superfícies para o número de condição de M0, κ =∥∥∥M0
∥∥∥ ∥∥∥(M0)−1∥∥∥ 45
4.5 Iso-superfícies para χ =∥∥∥(M0)−1
∥∥∥ . . . . . . . . . . . . . . . . . . . . 46
4.6 Curvas de nível para |v1Tv2|
‖v1‖‖v2‖ . . . . . . . . . . . . . . . . . . . . . . . 474.7 Iso-superfícies para ς: Cosseno Máximo entre as linhas da matriz M0. 474.8 Curvas de nível para
∥∥∥(M0)−1∥∥∥
estim.. . . . . . . . . . . . . . . . . . . . 48
4.9 Iso-superfícies para χ :=∥∥∥(M0)−1
∥∥∥estim.
. . . . . . . . . . . . . . . . . 49
5.1 Direção de descida viável d no ponto x. Espaço tangente T (x) . . . . 535.2 Condições de KKT . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
6.1 Desempenho dos quatro algoritmos: FDIPA (azul), FDIPA-FD(preto), FDIPA-AG (vermelho) e ALGA (verde). . . . . . . . . . . . 68
6.2 Curvas de nível da função de Barnes e trajeto dos quatro algoritmos. 686.3 Convergência dos quatro algoritmos no problema de Barnes. . . . . . 696.4 Estrutura reticulada de 10 barras: L=9.144 m (360 in). . . . . . . . . 706.5 Estrutura ótima obtida pelos algoritmos FDIPA-FD e FDIPA-AG
para a treliça de 10 barras. . . . . . . . . . . . . . . . . . . . . . . . . 716.6 Convergência dos algoritmos FDIPA-FD e FDIPA-AG para a treliça
de 10 barras. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 726.7 Projeto inicial do domo. . . . . . . . . . . . . . . . . . . . . . . . . . 73
x
6.8 Convergência dos algoritmos FDIPA-FD, FDIPA-AG e ALGA noproblema do domo. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75
7.1 Erro no gradiente aproximado sob ruído moderado, i.e., ‖φ‖Y ≤ 1%. . 797.2 Direção de busca do FDIPA. . . . . . . . . . . . . . . . . . . . . . . . 817.3 Direção de busca conservativa (azul) proposta na presença de
restrições com ruído. . . . . . . . . . . . . . . . . . . . . . . . . . . . 817.4 Convergência dos algoritmos FDIPA-FD, FDIPA-AG e FDIPA-AG-N
para a treliça de 10 barras no caso de ruído nas restrições. . . . . . . 867.5 Convergência dos algoritmos FDIPA-FD, FDIPA-AG e FDIPA-AG-N
para o domo de 52 barras no caso de ruído nas restrições. . . . . . . . 88
8.1 Materiais compósitos com diferentes formas dos constituintes. . . . . 928.2 Camadas da placa sanduíche. . . . . . . . . . . . . . . . . . . . . . . 948.3 Lâmina de material ortotrópico com direções principais associadas ao
referencial material (x1, x2, x3) e sistema de coordenadas global (x, y, z). 948.4 Convergência dos algoritmos FDIPA-FD e FDIPA-AG no problema
de identificação de propriedades. . . . . . . . . . . . . . . . . . . . . . 1028.5 Placa CFRP sem (esquerda) e com (direita) arranjos piezoelétricos. . 1058.6 Distribuções ótimas dos arranjos piezoelétricos (sensor/atuador). . . . 1078.7 Convergência dos algoritmos FDIPA, FDIPA-FD e FDIPA-AG
no problema de distribuição ótima de arranjos sensor/atuadorpiezoelétricos. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108
xi
Lista de Tabelas
4.1 Resumo dos resultados teóricos para o erro na interpolação polinomial. 42
6.1 Resultados estatísticos problema de Barnes, para 10 execuçõesindependentes com o algoritmo ALGA. . . . . . . . . . . . . . . . . . 69
6.2 Resultados do número da avaliações da função objetivo necessáriaspara solucionar o problema de Barnes, partindo de diferentes pontosiniciais. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
6.3 Estrutura ótima obtida para a treliça de 10 barras (seções transversaisem in2). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
6.4 Classificação das barras para o problema do domo. . . . . . . . . . . 736.5 Resultado ótimo para o domo, coordenadas (m) e seção transversal
(m2). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 746.6 Resultado ótimo para o domo, frequência angular (rad/s). . . . . . . 746.7 Resultados estatísticos domo, para 10 execuções independentes com
o algoritmo ALGA. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75
7.1 Estrutura obtida para a treliça de 10 barras sob ruído moderado nastensões e deslocamentos (seções transversais em in2). . . . . . . . . . 85
7.2 Estrutura obtida para o domo sob ruído moderado nas frequênciasnaturais, coordenadas (m) e seção transversal (cm2). . . . . . . . . . 87
8.1 Dimensões da placa CFRP. . . . . . . . . . . . . . . . . . . . . . . . 1018.2 Resultado ótimo para a placa de laminado, parâmetros adimensionais. 1018.3 Resultado ótimo da distribuição dos arranjos de sensor/atuador
piezoelétricos. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106
xii
Capítulo 1
Introdução
Os algoritmos de otimização em engenharia mais confiáveis e eficientes são osbaseados no gradiente. É bem sabido que podemos encontrar muita informação útilcontida nas derivadas de uma função que queremos otimizar. Afinal de contas, acaracterização matemática padrão de um mínimo local, fornecida pelas condiçõesnecessárias de primeira ordem, requer, para funções continuamente diferenciáveis,que as derivadas de primeira ordem sejam zero. No entanto, por vários motivossempre houve muitos casos onde não dispomos de (ao menos algumas) derivadasou estas não são confiáveis. Não obstante, ainda sob essas circunstâncias podeser desejável realizar otimização. Consequentemente, sempre tem sido necessária aincorporação de diferenças finitas ou uma classe de técnicas de otimização conhecidacomo métodos de otimização sem derivadas (DFO, Derivative-Free Optimization).De fato, a otimização sem derivadas é uma das áreas mais importante, aberta,e desafiante na engenharia e nas ciências da computação, com enorme potencialprático. A razão de ser um desafio é que, desde o ponto de vista de otimização,desistimos de muita informação ao não dispor das derivadas. Sua importânciaprática atual vem da crescente necessidade de resolver problemas de otimizaçãodefinidos por funções para as quais não existem derivadas ou elas são disponíveisa um custo proibitivo. A complexidade crescente na modelagem matemática,alta sofisticação da computação científica, e a abundância de códigos privados sãoalgumas das razões pelas quais a otimização sem derivadas é atualmente uma áreade grande demanda.
A maioria das aplicações de otimização em engenharia trata de problemas nãolineares com restrições. Na área de otimização estrutural, entre outras mais, existemrestrições que precisam ser respeitadas durante todo o processo de otimização. Daquia importância de desenvolver algoritmos que tenham esta característica.
Até a presente data, na comunidade de otimização, não se dispõe de uma técnicade otimização eficiente sem derivadas que tenha sido projetada para gerar sempresoluções viáveis (i.e. que cumpram com as restrições não lineares impostas). A
1
maioria das técnicas mais conhecidas para DFO foram concebidas para problemasirrestritos. As poucas técnicas que tentam adaptar algoritmos para problemasrestritos fazem uso de alguma estratégia de penalização ou a resolução aproximadade subproblemas de otimização. Ditas estratégias possibilitam o trabalho comrestrições, mas de uma maneira precária e não garantem a geração de pontos viáveis.Tudo isto motivou a realização do presente trabalho, onde se apresentam novastécnicas de otimização restrita que aproximam as derivadas para suprir a falênciamencionada anteriormente. Opta-se por ter como ferramenta base os algoritmosde direções viáveis baseados no gradiente, pela sua eficiência e robustez. Assimsendo, resultados matemáticos da teoria de aproximação estudados neste trabalhopossibilitaram o desenvolvimento de um novo algoritmo de otimização não linearonde não se precisa do fornecimento das derivadas.
No Capítulo 2 apresenta-se uma revisão não exaustiva das técnicas estabelecidasde otimização sem derivadas, de maneira a ilustrar as características principais decada uma delas.
No Capítulo 3 se apresenta o estudo teórico com as características que tornamas aproximações polinomiais interpolantes a nossa principal eleição para a obtençãode gradientes aproximados a serem integrados com a técnica de direções viáveistomada como base. A estimação das derivadas levando em conta a qualidade dasaproximações é apresentada no Capítulo 4.
Posteriormente, no Capítulo 5, Seção 5.2, será explicada uma versão dométodo de direções viáveis desenvolvido por HERSKOVITS [1] que constitui partefundamental da nossa nova técnica e lhe proporciona características importantestais como obter soluções viáveis em cada iteração e a convergência garantida aum mínimo local para problemas sob a hipótese de regularidade. No restante docapítulo se apresentará o novo algoritmo de otimização restrita com aproximação dederivadas.
No Capítulo 6 se apresentam exemplos numéricos resolvidos com o algoritmoproposto e comparado com outros algoritmos. Dentro dos exemplos tem-se umexemplo algébrico muito conhecido e diversos problemas de otimização estrutural.
Um dos objetivos principais desta tese é que a técnica proposta possa abordarproblemas de otimização onde há ruido moderado nas funções. Com este propósito,no Capítulo 7 se propõem modificações, consolidadas numa nova versão específicado algoritmo para a abordagem deste tipo de problema.
Finalmente, no Capítulo 8 se abordam duas aplicações do método proposto emestruturas de materiais compósitos laminados. No primeiro problema se realizaa identificação de propriedades mecânicas do compósito laminado. No segundoproblema, se apresenta uma nova formulação da distribuição ótima de atuadorespiezoelétricos para a dissipação de vibrações.
2
Capítulo 2
Otimização sem derivadas
Com o aumento na escala e na dificuldade das aplicações em engenharia, métodossofisticados de otimização baseados em derivadas tornaram-se necessários. Como crescimento e desenvolvimento de métodos de otimização não linear baseadosem derivadas tornou-se evidente que os problemas de grande porte podem serresolvidos eficientemente, mas somente se há disponibilidade de informação precisadas derivadas. Os usuários começaram a providenciar tais derivadas (programandoelas ou aplicando ferramentas de derivação automática em seus programas) ouestimando as derivadas por diferenças finitas.
No entanto, há situações em que não funciona nenhuma dessas abordagens paraa obtenção das derivadas. Por exemplo, no caso de códigos registrados, reescrevertais códigos, ou adicionar a estes o necessário para prover derivadas de primeiraordem, pode ser uma tarefa que demanda muito tempo. As técnicas de derivaçãoautomática (ver, e.g., [2, 3]) também não podem ser aplicadas em todos os casos.Em particular, se a função objetivo é calculada usando um pacote de simulação tipocaixa preta, é tipicamente impossível aplicar derivação automática, e ainda no casoem que o cômputo é mais acessível, assuntos de registro ou propriedade intelectualpodem fazer tal método inaceitável.
Há também duas situações onde aplicar aproximações das derivadas pordiferenças finitas é inapropriado: quando as avaliações da função têm alto custocomputacional e quando elas contêm ruido. No primeiro caso, realizar o númeronecessário de avaliações da função para prover uma estimativa do gradiente podeser proibitivo. No segundo caso, o gradiente estimado pode não servir para nada.Ironicamente, com a crescente sofisticação de hardware e algoritmos matemáticos,situações como essa são cada vez mais frequentes. A razão é simplesmenteque enquanto, antes, a simulação de sistemas complexos foi uma tarefa difícil edispendiosa e não fornecia um ambiente bom o suficiente para otimização, agoratais simulações estão se tornando mais habituais e também mais precisas; portanto,a otimização de sistemas complexos está se tornando uma possibilidade razoável. A
3
demanda crescente de métodos sofisticados de otimização sem derivadas aumentousua visibilidade na comunidade de otimização não linear.
Não é de se surpreender, como já foi indicado, que há desvantagens consideráveisao não ter informação da derivada, por isso não se pode esperar que o desempenhodos métodos sem derivadas seja comparável ao dos métodos baseados nas derivadas.Em particular, a dimensão dos problemas que atualmente podem ser resolvidoseficientemente por métodos sem derivadas ainda é relativamente pequeno. Oscritérios de parada são também um desafio na ausência de derivadas, quando asavaliações da função contém ruido e/ou são dispendiosas. Portanto, uma soluçãoquase ótima obtida através de um método sem derivadas é muitas vezes menosprecisa do que a obtida por um método com derivadas, assumindo que a informaçãoda derivada está disponível. Mesmo assim, para muitas das aplicações que existemhoje estas limitações são aceitáveis.
2.1 Limitações da otimização sem derivadas
Uma das principais limitações dos métodos sem derivadas é que, numcomputador com processamento em série, geralmente não é razoável tentar otimizarproblemas com dezenas de variáveis, embora algumas das técnicas mais recentespermitam trabalhar com problemas sem restrições com centenas de variáveis(ver, e.g., [4]). Além disso, mesmo em problemas relativamente simples e bemcondicionados normalmente não é realista esperar soluções precisas. A convergênciaé tipicamente mais lenta. Tipicamente pode esperar-se uma velocidade deconvergência local que seja mais próxima a ser linear do que quadrática, podendo-sepreferir uma terminação antecipada. Os métodos de otimização sem derivadas, porestarem capacitados para abordar unicamente problemas de relativamente poucasvariáveis (menos de 50), o esquema habitual para abordar problemas com maisvariáveis consiste em usar métodos estatísticos como análise de variância (ver, e.g.,[5]) para determinar as variáveis mais críticas e otimizar apenas sobre elas. Tambémpode ser razoável tirar vantagem da simplicidade relativa de alguns algoritmos semderivadas, como os métodos de busca direta direcional, e usar processamento emparalelo.
É geralmente aceito que os métodos de otimização sem derivadas têm a habilidadede encontrar ótimos locais “bons” no seguinte sentido. Se tiver uma função comum aparente número elevado de ótimos locais, talvez por causa do ruido, então asestrategias sem derivadas, pela sua crueza relativa, têm a tendência a inicialmenteir a regiões com menor função objetivo (por causa de sua cegueira) e posteriormentetendem a suavizar a função, onde um método mais sofisticado pode achar o ótimolocal mais próximo. A tendência a “suavizar” as funções é a razão pela qual são
4
efetivos para funções com ruido moderado. Há muitas situações onde os métodos semderivadas são as únicas estratégias adequadas, capazes de ter melhores resultadosdo que os algoritmos heurísticos e provendo suporte teórico para a convergência.
Há, no entanto, classes de problemas para os quais o uso de métodos livres dederivadas baseados na matemática não são adequados. Tipicamente, os métodosrigorosos para tais problemas precisariam de uma desmesurada quantidade detrabalho que cresce exponencialmente com o tamanho do problema. Esta categoriade problemas incluem problemas de otimização global de meio e grande porte comou sem derivadas, problemas que além de não ter derivadas disponíveis também nãosão suaves, problemas gerais não lineares de grande porte com variáveis discretasincluindo muitos outros de otimização combinatória (chamados problemas NPdifíceis), e problemas de otimização estocástica.
Para alguns dos casos citados anteriormente, costuma-se usar algoritmosheurísticos, tais como recozimento simulado [6], algoritmos genéticos e outrosalgoritmos evolutivos [7, 8], redes neurais artificiais [9], métodos de busca tabu[10], e métodos de enxame de partículas ou métodos baseados em população [11],incluindo técnicas de enumeração. A comunidade científica considera eles comométodos de último recurso; ou seja, aplicáveis a problemas onde o espaço de buscaé necessariamente grande, complexo, ou pobremente compreendido.
2.2 Classificação das técnicas de otimização semderivadas
Os métodos usados em otimização sem derivadas podem ser classificados emtrês categorias: métodos de busca direta, métodos baseados em modelos, e métodosheurísticos.
Os métodos de busca direta exploram o domínio usando regras sistemáticas.Seguindo a LEWIS et al. [12] estes podem se classificar em três subcategorias:
? Métodos de busca por padrão: caracterizados por uma sequência demovimentos exploratórios que considera o comportamento da função objetivono padrão de pontos que ficam numa grade regular. Os movimentosexploratórios consistem numa estratégia sistemática para visitar os pontosna grade na vizinhança imediata da iteração atual. O algoritmo de buscapor padrão (Pattern Search Algorithm) de HOOKE e JEEVES [13] éprovavelmente o mais famoso. TORCZON [14] desenvolveu uma teoria deconvergência global para os métodos de busca por padrão, sob condiçõesadequadas.
5
? Métodos de busca por simplex: caracterizados por um esquema simples queguia a busca. O primeiro método de busca por simplex é devido a SPENDLEYet al. [15]. Tal método é baseado na observação de que para identificar umadireção de descida (ou acenso) são necessárias não mais do que n+1 avaliaçõesda função objetivo, sendo n o número de variáveis de projeto. A ideia destesmétodos é atualizar um simplex substituindo o pior vértice por um novo. Esteúltimo pode ser obtido por transformações geométricas do simplex; tais comoexpansões, reflexões ou contrações. O algoritmo proposto por NELDER eMEAD [16] é o mais famoso e o mais popular. No entanto MCKINNON [17]mostrou, através de um contra-exemplo simples dado por uma função objetivosuave e convexa em duas dimensões, que o algoritmo de Nelder-Mead nãoconverge.
? Métodos com conjuntos adaptativos de direções de busca: propostos a saberpor ROSENBROCK [18] e POWELL [19]. Eles tentam acelerar a busca aoconstruir direções projetadas para usar informação da curvatura da funçãoobjetivo recolhida durante a busca.
Na literatura sobressaem duas abordagens baseadas em modelos:
? Métodos de região de confiança (Trust Region) baseados em modelosquadráticos que interpolam a função objetivo: estes métodos são muitoefetivos. Em cada iteração constroem modelos quadráticos em torno daiteração atual. Se assume que este modelo aproxima suficientemente bema função objetivo numa região chamada região de confiança. O raio daregião de confiança precisa ser atualizado em cada iteração e um minimizadordo modelo quadrático dentro da região de confiança é considerado comocandidato para a próxima iteração. Se a redução da função prevista pelomodelo corresponde à redução efetiva da função objetivo suficientemente bem,o candidato é aceito e o raio da região de confiança é possivelmente aumentado.Caso contrário, o candidato é rejeitado e o raio da região de confiança édiminuído. O proceso começa novamente até satisfazer um certo critério deconvergência. DFO (Derivative-Free Optimization, ver CONN e TOINT [20]),UOBYQA (Unconstrained Optimization BY Quadratic Approximation, verPOWELL [21]) e NEWUOA (NEW Unconstrained Optimization Algorithm,ver POWELL [22]) podem ser considerados como algoritmos de referência.
? A plataforma de gestão de modelos substitutos (SMF, Surrogate ManagementFramework), proposta por BOOKER et al. [23], é baseada numa buscapor padrão. Propõe uma plataforma geral para gerar uma sequênciade aproximações da função objetivo e gerir o uso dessas aproximações
6
como substitutos no processo de otimização. O algoritmo converge a umminimizador de uma função objetivo cara sujeita a restrições simples.
Os métodos heurísticos em sua grande maioria são inspirados na natureza. Hojeo termo “heurística” é usado para descrever um método “que, baseado na experiênciaou julgamento, parece conduzir a uma boa solução de um problema, mas quenão garante produzir uma solução ótima”. Esse termo pode ser considerado comoassociado a um conhecimento circunstancial, não verificável matematicamente. Noentanto, os métodos heurísticos são usados em algumas situações extremas pelasua simplicidade, capacidade de trabalhar com funções de qualquer tipo, e suatendência a obter ótimos globais. A maioria dos métodos heurísticos são inspiradosna natureza. Diz-se que um algoritmo heurístico é bio-inspirado quando suas regrasde busca tentam simular alguns aspectos do comportamento de seres vivos ou defenômenos naturais como, por exemplo, os já citados algoritmos genéticos [7, 8],o algoritmo das formigas [24, 25], o algoritmo dos pássaros [26], o algoritmo dasabelhas [27] ou dos vermes luminosos [28].
Os métodos heurísticos também podem se classificar em duas subcategorias:
? Nos métodos construtivos ou de trajetória parte-se de um conjunto soluçãovazio e segue-se acrescentando elementos a esse conjunto até obter uma soluçãoviável para o problema. Um exemplo de método construtivo é o algoritmo dasformigas, pois este descreve uma trajetória no espaço de busca durante a suaexecução.
? Os métodos populacionais partem de um conjunto de soluções iniciais (apopulação inicial) e tentam encontrar uma solução melhor modificando apouplação inicial. Exemplos de métodos baseados em população são osalgoritmos genéticos (GA, Genetic Algorithms) e o recozimento simulado (SA,Simulated Annealing).
7
Capítulo 3
Modelos aproximados eamostragem
Nesta tese o foco é sobre modelos polinomiais a fim de obter aproximações dogradiente e desta maneira capacitar algoritmos convencionais robustos a trataremproblemas onde as derivadas não estão disponíveis.
Os modelos estudados aqui vão ser polinômios de baixa ordem. Basicamenteserão estudados os modelos lineares em detalhe, em menos detalhe se apresentarãoas propriedades dos modelos quadráticos e dos modelos não polinomiais.
3.1 Modelos lineares
Os modelos lineares podem dispor da precisão necessária para realizaraproximações razoáveis do gradiente e sua obtenção é relativamente mais barata(requer menos amostras, i.e. menos avaliações da função) do que no caso de modelosde ordem superior.
3.1.1 Modelos lineares de interpolação e regressão
O objetivo desta seção é estabelecer os fundamentos teóricos para usar modelosde interpolação e regressão como aproximações das funções (objetivo e restrições) oude suas derivadas num algoritmo de otimização. Para garantir a convergência globalde um algoritmo de otimização que usa o modelo das funções ou aproximações desuas derivadas é necessário garantir uma certa qualidade desse modelo. Quando ummodelo consiste numa expansão em série de Taylor truncada de primeira ou segundaordem, podemos obter facilmente uma medida da qualidade do modelo a partirdos limites do erro típicos de uma expansão em série de Taylor. No caso de umainterpolação ou regressão linear existem limites semelhantes, mas, ao contrário doslimites da expansão em série de Taylor, eles dependem não só do ponto de expansão
8
e da função que está sendo aproximada como também dos pontos do conjuntode amostras (como será mostrado ao longo desta seção). Assim, para manter aqualidade do modelo de interpolação ou regressão é necessário entender e mantera qualidade do conjunto de pontos amostrados (para interpolação ou regressão).Nesta seção, serão examinadas vários indicadores que caracterizam a qualidade deum conjunto de amostras. Serão estudadas as relações entre essas constantes e opapel que elas têm no limite do erro entre o polinômio de aproximação e a funçãoreal (e entre suas derivadas).
O conjunto de amostras considerado aqui não necessariamente será formadopor um conjunto de direções predefinidas, ditas direções sendo usadas comumenteem algoritmos de busca direta direcional (directional direct-search). Considere oconjunto de amostras Y = y0,y1, . . . ,yp em Rn. O modelo mais simples quepodemos pensar baseado em n + 1 pontos de amostragem (p = n) é um modelo deinterpolação.
Interpolação linear
Seja m(x) um polinômio de grau d = 1 que interpola f nos pontos em Y , i.e.,satisfazendo as condições de interpolação
m(yi) = f(yi), i = 0, . . . , n. (3.1)
Pode se expressar m(x) na forma
m(x) = α0 + α1x1 + · · ·+ αnxn,
usando as bases polinomiais φ = 1, x1, . . . , xn como bases para o espaço P1n de
polinômios lineares de grau 1. Pode se reescrever (3.1) como
1 y01 · · · y0
n
1 y11 · · · y1
n... ... ... ...1 yn1 · · · ynn
α0
α1...αn
=
f(y0)f(y1)
...f(yn)
. (3.2)
A matriz deste sistema linear é denotada por
M = M(φ, Y ) =
1 y0
1 · · · y0n
1 y11 · · · y1
n... ... ... ...1 yn1 · · · ynn
. (3.3)
9
Neste documento se escreve M(φ, Y ) para destacar que M depende das bases φ edo conjunto de amostras Y .
Definição 3.1.1. O conjunto Y = y0,y1, . . . ,yn é bem posto para interpolaçãolinear em Rn se a matriz M(φ, Y ) correspondente é não singular.
A Definição 3.1.1 é independente das bases escolhidas. Em outras palavras, seY é bem posto para uma base φ então ele é bem posto para qualquer outra base emP1n. A definição de m(x) também é independente da base escolhida. É fácil ver que
o conjunto de amostras Y é bem posto para interpolação linear se o polinômio deinterpolação linearm(x) = α0+α1x1+· · ·+αnxn está definido de maneira única [29].
Regressão linear
Quando o número de pontos p + 1 do conjunto de amostras excede em maisde 1 a dimensão n do espaço de amostragem, poderia não ser possível ajustar umpolinômio linear. Neste caso uma opção é usar regressão linear e obter os coeficientesda regressão (por mínimos quadrados) polinomial m(x) = α0 + α1x1 + · · · + αnxn
através da solução do sistema de mínimos quadrados
1 y01 · · · y0
n
1 y11 · · · y1
n... ... ... ...1 yp1 · · · ypn
α0
α1...αn
=
f(y0)f(y1)
...f(yp)
. (3.4)
De novo denotamos por M(φ, Y ) a matriz deste sistema (possivelmentesobredeterminado) de equações.
Generalizando a Definição 3.1.1 obtém-se a definição de conjunto bem posto pararegressão linear.
Definição 3.1.2. O conjunto Y = y0,y1, . . . ,yp, com p > n, é bem posto pararegressão linear em Rn se a matriz M(φ, Y ) correspondente tem posto por colunasigual a n, i.e. se todas suas colunas são linearmente independentes.
No caso de regressão linear também pode-se provar que se o conjunto Y é bemposto para uma base φ também é bem posto para qualquer outra base em P1
n, e quem(x) é independente da base escolhida [29]. Finalmente, pode-se ressaltar que oconjunto de amostras é bem posto para regressão linear se o polinômio de regressãolinear m(x) está definido de maneira única.
10
3.1.2 Limites do erro para interpolação e regressão linear
Limites do erro para interpolação linear
Fazendo α0 = c e αi = gi, i = 1, . . . , n pode-se reescrever o polinômio deinterpolação linear na forma
m(y) = c+ gTy.
Dependendo da finalidade da interpolação o interesse pode ser na medida daqualidade de m(y) como uma aproximação para f(y) ou na medida da qualidade de∇m(y) como uma aproximação para ∇f(y). Primeiro se abordará a qualidade dogradiente g = ∇m(y) do modelo, como uma aproximação para ∇f(y). Considerar-se-á que os pontos de interpolação y0,y1, . . . ,yn estão numa bola de raio ∆ comcentro em y0. Na prática, pode se escolher
∆ = ∆(Y ) = max1≤i≤n
‖yi − y0‖.
O interesse é na qualidade de ∇m(y) na bola de raio ∆ centrada em y0. Assuposições necessárias para esse resultado estão resumidas a seguir.
Suposição 3.1.1. Assume-se que Y = y0,y1, . . . ,yn ⊂ Rn é um conjunto bemposto de pontos de amostragem (no sentido de interpolação linear) contido na bolaB(y0; ∆(Y )) de radio ∆ = ∆(Y ).
Além disso, assume-se que a função f é continuamente diferenciável num dominoaberto Ω contendo B(y0; ∆) e que ∇f é Lipschitz contínuo em Ω com constantev > 0.
A obtenção dos limites de erro é baseada na aplicação de um passo de eliminaçãoGaussiana ao sistema M(φ, Y )α = f(Y ) em (3.1). Após realizar tal passo chega-seao sistema
1 y01 · · · y0
n
0 y11 − y0
1 · · · y1n − y0
n... ... ... ...0 yn1 − y0
1 · · · ynn − y0n
α0
α1...αn
=
f(y0)
f(y1)− f(y0)...
f(yn)− f(y0)
, (3.5)
que pode ser representado por blocos como1 (y0)T
0 M0
α0
∇f
= f(y0)δ(f, Y )
(3.6)
11
com
M0 =
y1
1 − y01 · · · y1
n − y0n
... ... ...yn1 − y0
1 · · · ynn − y0n
,∇f = [α1 α2 . . . αn]T,
δ(f, Y ) = [f(y1)− f(y0) . . . f(yn)− f(y0)]T.
É evidente que M0 é não singular se M é não singular, já que det(M0) = det(M).Cabe notar que os pontos aparecem listados em M0 por linhas, o que favorecefatorações por linhas.
Os limites de erro para a aproximação a se obter estão em termos da matrizescalada
M0 = 1∆M0 =
y1
1−y01
∆ · · · y1n−y0
n
∆... ... ...yn1−y
01
∆ · · · ynn−y0n
∆
.Esta matriz M0 corresponde ao conjunto de amostras escalado
Y = y0/∆,y1/∆, . . . ,yn/∆ ⊂ B(y0/∆; 1),
que está contido na bola de raio 1 centrada em y0/∆.De [29] temos o seguinte teorema.
Teorema 3.1.1. Mantendo a Suposição 3.1.1. O gradiente do modelo deinterpolação linear satisfaz, para todos os pontos y em B(y0; ∆), um limite no erroda forma
‖∇f(y)−∇m(y)‖ ≤ κeg∆, (3.7)
onde κeg = v(1 + n12
∥∥∥(M0)−1∥∥∥ /2) e M0 = M0/∆.
Demonstração. Se o conjunto Y é bem posto, então a matriz M = M(φ, Y ) detamanho (n + 1) × (n + 1) é não singular, e por tanto também a matriz M0 detamanho n× n.
Primeiro se considera o gradiente de f no ponto y0. Subtraindo a primeiracondição de interpolação das outras n condições, obtém-se
(yi − y0)Tg = f(yi)− f(y0), i = 1, . . . , n. (3.8)
Usando a forma integral do teorema do valor médio
f(yi)− f(y0) =∫ 1
0(yi − y0)T∇f(y0 + t(yi − y0))dt, (3.9)
12
de (3.8) e substraindo (yi − y0)T∇f(y0) de ambos os lados de (3.9) tem-se
(yi−y0)Tg− (yi−y0)T∇f(y0) =∫ 1
0(yi−y0)T
[∇f(y0 + t(yi − y0))−∇f(y0)
]dt.
Usando o fato que∥∥∥∫ 1
0 ∇f(t)dt∥∥∥ ≤ ∫ 1
0 ‖∇f(t)‖ dt, e a desigualdade de cauchy-swartz|a · b| < ‖a‖ ‖b‖ pode-se escrever:
∣∣∣(yi − y0)T(∇f(y0)− g)∣∣∣ ≤ ∫ 1
0
∥∥∥yi − y0∥∥∥ ∥∥∥∇f(y0 + t(yi − y0))−∇f(y0)
∥∥∥ dt.Usando a propriedade de ∇f ser Lipchitz contínua, i.e. ‖∇f(y1)−∇f(y0)‖ ≤v ‖yi − y0‖ , tem-se
∣∣∣(yi − y0)T(∇f(y0)− g)∣∣∣ ≤ ∫ 1
0v∥∥∥yi − y0
∥∥∥ ∥∥∥t(yi − y0)∥∥∥ dt
≤ v∥∥∥yi − y0
∥∥∥2 ∫ 1
0tdt
≤ v
2∥∥∥yi − y0
∥∥∥2, i = 1, . . . , n.
Então, destas últimas n desigualdades, obtém-se∥∥∥M0(∇f(y0)− g)
∥∥∥ ≤ (n 12v/2)∆2,
daqui se conclui que∥∥∥∇f(y0)− g
∥∥∥ ≤ (n 12∥∥∥(M0)−1
∥∥∥ v/2)∆.
O limite do erro para qualquer ponto y na bola B(y0; ∆) é obtido facilmente dacontinuidade de Lipchitz do gradiente de f :
‖∇f(y)− g‖ ≤∥∥∥∇f(y)−∇f(y0)
∥∥∥+∥∥∥∇f(y0)− g
∥∥∥ ≤ (v + n12∥∥∥(M0)−1
∥∥∥ v/2)∆.
Assumindo um limite uniforme em∥∥∥(M0)−1
∥∥∥ independente de ∆, o erro nogradiente é linear em ∆. Também se pode ver que o erro no modelo de interpolaçãom(x) é quadrático em ∆ [29].
Teorema 3.1.2. Mantendo a Suposição 3.1.1. O modelo de interpolação linearsatisfaz, para todos os pontos y em B(y0; ∆), um limite no erro da forma
|f(y)−m(y)| ≤ κef∆2, (3.10)
onde κef = κeg + v2 e κeg é fornecido no Teorema 3.1.1.
13
Demonstração. Usando os mesmos argumentos como na demonstração do Teorema3.1.1 obtém-se
f(y)− f(y0) ≤ ∇f(y0)T(y− y0) + v
2∥∥∥y− y0
∥∥∥2.
Daqui se obtém
f(y)− f(y0)− gT(y− y0) ≤ (∇f(y0)− g)T(y− y0) + v
2∥∥∥y− y0
∥∥∥2.
O limite de erro (3.10) vem da combinação desta desigualdade com (3.7) e denotar que o termo constante no modelo pode ser escrito como c = f(y0)−gTy0.
Limites do erro para regressão linear
No caso de regressão está se considerando um conjunto de amostras Y =y0,y1, . . . ,yp com mais de n+ 1 pontos contidos na bola B(y0; ∆(Y )) de raio
∆ = ∆(Y ) = max1≤i≤p
‖yi − y0‖.
Reescreve-se o polinômio de regressão linear na forma
m(y) = c+ gTy,
onde c = α0 e gi = αi, i = 1, . . . , n, são as componentes da solução por mínimosquadrados do sistema (3.4).
Suposição 3.1.2. Assume-se que Y = y0,y1, . . . ,yp ⊂ Rn, com p > n, é umconjunto bem posto de pontos de amostragem (no sentido de regressão linear)contido na bola B(y0; ∆(Y )) de radio ∆ = ∆(Y ).
Além disso, assume-se que a função f é continuamente diferenciável num dominoaberto Ω contendo B(y0; ∆) e que ∇f é Lipschitz continuo em Ω com constantev > 0.
A obtenção dos limites no erro para a aproximação também é feita em termosda matriz escalada
M0 = 1∆M0 = 1
∆[y1 − y0 · · ·yp − y0
]T,
que corresponde ao conjunto de amostras escalado contido numa bola de raio 1centrado em y0/∆, i.e.,
Y = y0/∆,y1/∆, . . . ,yp/∆ ⊂ B(y0/∆; 1).
14
A prova dos limites segue exatamente os mesmos passos da prova feita para ocaso de interpolação linear. Por exemplo, para a aproximação do gradiente, umavez chegado ao ponto na prova onde
∥∥∥M0(∇f(y0)− g)∥∥∥ ≤ (p 1
2v/2)∆2,
ou, equivalentemente,∥∥∥M0(∇f(y0)− g)
∥∥∥ ≤ (p 12v/2)∆,
passa-se M0 ao lado direito através da sua inversa esquerda (M0)†1. Podem seestabelecer os limites no seguinte formato.
Teorema 3.1.3. Mantendo a Suposição 3.1.2. O gradiente do modelo de regressãolinear satisfaz, para todos os pontos y em B(y0; ∆), um limite no erro da forma
‖∇f(y)−∇m(y)‖ ≤ κeg∆,
onde κeg = v(1 + p12
∥∥∥(M0)†∥∥∥ /2) e M0 = M0/∆.
O modelo de regressão linear satisfaz, para todos os pontos y em B(y0; ∆), umlimite no erro da forma
|f(y)−m(y)| ≤ κef∆2,
onde κef = κeg + v2 .
3.2 Modelos não lineares
Os modelos lineares discutidos na seção anterior são atraentes pela simplicidadena sua construção e analise. Mas como acontece com qualquer modelo linear elesnão capturam informação da curvatura da função que estão aproximando a menosque seja plana). Para conseguir melhores velocidades de convergência local em geralé importante considerar modelos não lineares.
O modelo polinomial quadrático pode ser considerado o modelo não linearmais simples, no entanto, é muitas vezes o mais eficiente. Nesta seção serãoreportadas as propriedades dos modelos de interpolação e regressão quadráticos.Modelos não polinomiais, tais como funções de base radial, também podem serusados efetivamente em otimização sem derivadas. Esses modelos serão cobertosbrevemente na Seção 3.3.
1A†. denota a inversa generalizada de Moore-Penrose da matriz A, que pode ser expressa peladecomposição em valor singular de A para qualquer matriz A real ou complexa. No contextoatual, onde M0 tem posto completo por colunas, tem-se (M0)† = ((M0)TM0)−1(M0)T.
15
3.2.1 Interpolação de modelos não lineares
A abordagem feita nesta seção aplica-se à interpolação de polinômios dequalquer ordem. Mas o caso quadrático é o mais representativo para os propósitosdesta pesquisa.
Conceitos básicos em interpolação
Bases polinomiaisSeja Pdn, o espaço de polinômios de grau menor ou igual a d em Rn. Seja p1 = p+1
a dimensão deste espaço. É bem sabido que para d = 1, p1 = n+1 e que para d = 2,p1 = (n+ 1)(n+ 2)/2.
Uma base φ = φ0(x), φ1(x), . . . , φp(x) de Pdn é um conjunto de p1 polinômiosde grau menor ou igual a d que gera Pdn. Já que φ é uma base em Pdn, então qualquerpolinômio m(x) ∈ Pdn pode ser escrito como m(x) = ∑p
j=0 αjφj(x), onde os αj sãoalguns coeficientes reais.
Há várias bases polinomiais interessantes a serem consideradas para váriasaplicações. Serão focadas somente aquelas bases que são de interesse no contextodesta tese. A base polinomial mais simples e comum é a base de monômios,conhecida como a base natural.
A base natural φ pode ser descrita convenientemente usando multi-índices. Sejao vetor αi = (αi1, . . . , αin) ∈ Nn
0 um multi-índice, e para qualquer x ∈ Rn, seja xαi
definido como
xαi =n∏j=1
xαij
j .
Também, define-se
∣∣∣αi∣∣∣ =∑j=1
αij e (αi)! =∏j=1
(αij!).
Então os elementos da base natural são
φi(x) = 1(αi)!x
αi , i = 0, . . . , p,∣∣∣αi∣∣∣ ≤ d.
As bases naturais podem ser escritas como:
φ = 1, x1, x2, . . . , xn, x21/2, x1x2, . . . , x
d−1n−1xn/(d− 1)!, xdn/d!.
Por exemplo, quando n = 3 e d = 2, tem-se
φ = 1, x1, x2, x3, x21/2, x1x2, x
22/2, x1x3, x2x3, x
23/2.
16
As bases naturais são também as bases dos polinômios que resultam de umaexpansão de Taylor. Por exemplo, assumindo que f(x) seja suave, o modelo deTaylor de ordem d = 2 em R3, centrado no ponto y, é o seguinte polinômio emz1, z2, e z3 (os elementos da base natural foram escritos dentro de parêntesesquadrados para uma melhor identificação):
f(y)[1] + ∂f
∂x1(y)[z1] + ∂f
∂x2(y)[z2] + ∂f
∂x3(y)[z3]
+∂2f
∂x21(y)[z2
1/2] + ∂2f
∂x1x2(y)[z1z2] + ∂2f
∂x22(y)[z2
2/2]
+ ∂2f
∂x1x3(y)[z1z3] + ∂2f
∂x2x3(y)[z2z3] + ∂2f
∂x23(y)[z2
3/2].
Quando d = 1, obtém-se φ = 1, x1, . . . , xn, que já apareceu na Seção 3.1.1.
Interpolação polinomialPode-se dizer que o polinômio m(x) interpola a função f(x) num ponto dado y
se m(y) = f(y). Assumindo que foi fornecido o conjunto Y = y0,y1, . . . ,yp ⊂ Rn
de pontos de interpolação, e seja m(x) o polinômio de grau menor ou igual a d
que interpola a função dada f(x) nos pontos em Y . Determinando os coeficientesα0, . . . , αp se obtém o polinômio de interpolação m(x). Os coeficientes α0, . . . , αp
podem ser determinados através das condições de interpolação
(m(yi)
)=
p∑j=0
αjφj(yi) = f(yi), i = 0, . . . , p. (3.11)
As condições (3.11) formam um sistema linear em termos dos coeficientes deinterpolação, que pode ser escrito na forma matricial como
M(φ, Y )αφ = f(Y ),
onde
M(φ, Y ) =
φ0(y0) φ1(y0) · · · φp(y0)φ0(y1) φ1(y1) · · · φp(y1)
... ... ... ...φ0(yp) φ1(yp) · · · φp(yp)
,
αφ =
α0
α1...αp
, e f(Y ) =
f(y0)f(y1)
...f(yp)
.
17
Por exemplo, quando n = d = 2 e p = 5, a matriz M(φ, Y ) fica
1 y01 y0
2 (y01)2/2 y0
1y02 (y0
2)2/21 y1
1 y12 (y1
1)2/2 y11y
12 (y1
2)2/21 y2
1 y22 (y2
1)2/2 y21y
22 (y2
2)2/21 y3
1 y32 (y3
1)2/2 y31y
32 (y3
2)2/21 y4
1 y42 (y4
1)2/2 y41y
42 (y4
2)2/21 y5
1 y52 (y5
1)2/2 y51y
52 (y5
2)2/2
.
Para que o sistema anterior tenha solução única a matriz M(φ, Y ) tem que sernão singular.
Definição 3.2.1. O conjunto Y = y0,y1, . . . ,yp é bem posto para interpolaçãopolinomial em Rn se a matriz M(φ, Y ) correspondente é não singular para algumabase φ em Pdn.
É fácil mostrar que se M(φ, Y ) é não singular para alguma base φ, então énão singular para qualquer base de Pdn. Sob essas circunstâncias, pode-se mostrartambém que o polinômio de interpolação m(x) existe e é único [29].
Dado um espaço de polinômios Pdn e uma base φ, seja
φ(x) = [φ0(x), φ1(x), . . . , φp(x)]T
um vetor em Rp1 formado pelos valores dos elementos da base polinomial em x (φpode ser visto como um mapeamento de Rn para Rp1).
Lema 3.2.1. Dada uma função f : Rn → Rp1 e um conjunto bem posto Y ∈ Rn, opolinômio de interpolação m(x) existe e é único.
Demonstração. Já que Y é bem posto, então é evidente que m(x) existe e é únicopara uma base φ(x). Agora é necessário mostrar quem(x) é independente da escolhade φ(x). Considere outra base ψ(x) = PTφ(x), onde P é uma matriz p1 × p1 nãosingular. Claramente, tem-se que M(ψ, Y ) = M(φ, Y )P. Seja αφ uma solução de(3.11) para a base φ, e seja αψ uma solução de (3.11) para a base ψ. Então, paraqualquer lado direito f(Y ), tem-se que αψ = P−1(M(φ, Y ))−1f(Y ) e
αTψψ(x) = f(Y )T(M(φ, Y ))−TP−TPTφ(x)
= f(Y )T(M(φ, Y ))−Tφ(x) = αTφφ(x).
Conclui-se então que m(x) = αTψψ(x) = αT
φφ(x) para todo x.
18
A matriz M(φ, Y ) é singular se existe γ ∈ Rp1 tal que γ 6= 0 e M(φ, Y )γ = 0 eisso implica que existe um polinômio, de grau menor ou igual a d, expressado como
m(x) =p∑j=0
γjφj(x),
tal que m(y) = 0 para todo y ∈ Y. Em outras palavras, M(φ, Y ) é singular se ospontos de Y fazem parte de uma “variedade polinomial” de grau menor ou iguala d. Por exemplo, seis pontos numa curva de segunda ordem em R2 formam umconjunto não bem posto para interpolação quadrática2.
Agora que foram estabelecidos critérios para saber se um conjunto deinterpolação é bem posto ou não, é natural se perguntar se essas condições podemser estendidas. Por exemplo, já que a não singularidade de M(φ, Y ) é um indicadorde ser bem posto, poderia o número de condição de M(φ, Y ) ser um indicador dequão bem posto é o conjunto?. Em termos gerais, a resposta, é “não” já que onúmero de condição de M(φ, Y ) depende da escolha de φ. Além disso, para umconjunto de interpolação Y dado, pode-se escolher a base φ de maneira que onúmero de condição de M(φ, Y ) possa igualar qualquer número entre 1 e +∞.Também, para qualquer escolha fixa de φ, o número de condição de M(φ, Y )depende do escalado de Y . Portanto, o número de condição de M(φ, Y ) não é umaboa caracterização do nível em que o conjunto de pontos é bem posto. No entanto,pode ser demonstrado que para uma escolha específica de φ, chamada de basenatural, e para Y , uma versão escalada de Y , o número de condição de M(φ, Y ) éuma medida significativa de quão bem posto é Y .
Obtenção dos limites no erro em termos do número de condiçãoA seguir serão mostrados limites no erro entre o polinômio de interpolação
e a função sendo interpolada, entre seus gradientes e, possivelmente, entre suasHessianas, expressados em termos do número de condição de M(φ, Y ).
O número de condição de M é denotado por cond(M) = ‖M‖∥∥∥M−1
∥∥∥. Pode-sedemonstrar que ‖M‖ está limitado em termos da constante p1 (no caso quadrático1 ≤ ‖M‖ ≤ p
321 ). Então para obter um limite em termos de cond(M) é suficiente
obter ele em termos de∥∥∥M−1
∥∥∥.Considere-se a interpolação de uma função f(y) através do polinômio m(y), que
pode ser escrito usando a base natural como:
m(y) =p∑
k=0γkφk(y).
2O fato de que os pontos fazem parte de uma curva de segunda ordem implica que as colunasde M(φ, Y ) podem ser combinadas com alguns coeficientes, não todos zero, para obter um vetorzero.
19
O conjunto de interpolação satisfaz Y = y0,y1, . . . ,yp ⊂ B(y0; ∆(Y )), ondeB(y0; ∆(Y )) é a menor bola fechada, com centro em y0 e raio ∆ = ∆(Y ), contentoY . Os coeficientes γk, k = 0, . . . , p, estão definidos pelo sistema linear resultantedas condições (3.11).
Será considerado em detalhe apenas o caso quadrático (d = 2). O caso linear(d = 1) já foi visto na Seção 3.1.2. Os limites no erro para o caso quadrático sãoestimados sob as seguintes suposições.
Suposição 3.2.1. Assume-se que Y = y0,y1, . . . ,yp ⊂ Rn, com p1 = p + 1 =(n+ 1)(n+ 2)/2, é um conjunto de amostras bem posto (no sentido de interpolaçãoquadrática, d = 2) contido na bola B(y0; ∆(Y )) de raio ∆ = ∆(Y ).
Adicionalmente, assume-se que a função f é duas vezes continuamentediferenciável num domínio aberto Ω contendo B(y0; ∆(Y )) e que ∇2f é Lipschitz-contínua em Ω com constante v2 > 0.
Sob essas suposições é possível obter o modelo de interpolação quadrática
m(y) = c+ gTy + 12yTHy = c+
n∑j=1
gjyj +∑
1≤k,l≤nhklykyl, (3.12)
onde H é uma matriz simétrica de ordem n. Os coeficientes desconhecidosc, g1, . . . , gn, e hkl, 1 ≤ l ≤ k ≤ n, estão definidos de maneira única pelas condiçõesde interpolação (3.11) já que o conjunto de amostras é bem posto.
Seja Q a matriz formada pelas ultimas p linhas e colunas de M(φ, Y ), e suaversão escalada Q dada por.
Q = Q
D−1∆ 00 D−1
∆2
, (3.13)
onde D∆ é uma matriz diagonal de dimensão n com ∆ na diagonal e D∆2 é umamatriz diagonal de dimensão p−n com ∆2 em sua diagonal. Esta matriz escalada éo mesmo que a matriz Q correspondente ao conjunto escalado Y = Y/∆ ⊂ B(0; 1).
O seguinte teorema declara os limites no erro no caso quadrático [29]. Como eraesperado, o limites no erro são lineares em ∆ para as derivadas segundas, quadráticosem ∆ para las primeiras derivadas, e cúbicos em ∆ para os valores da função, onde∆ = ∆(Y ) é o raio da menor bola que contém Y . Quanto menor for
∥∥∥Q−1∥∥∥ o
conjunto Y vira bem posto e melhor será a estimativa dos erros.
Teorema 3.2.2. Mantendo a Suposição 3.2.1. Então, para todos os pontos y emB(y0; ∆(Y )), tem-se que
? o erro entre a Hessiana do modelo de interpolação quadrática e a Hessiana dafunção satisfaz ∥∥∥∇2f(y)−∇2m(y)
∥∥∥ ≤ κeh∆,
20
? o erro entre o gradiente do modelo de interpolação quadrática e o gradiente dafunção satisfaz
‖∇f(y)−∇m(y)‖ ≤ κeg∆2,
? o erro entre o modelo de interpolação quadrática e a função satisfaz
|f(y)−m(y)| ≤ κef∆3,
onde κeh, κeg, e κef são dadas por
κeh = 3√
2p 12v2
∥∥∥Q−1∥∥∥ /2
κeg = 3(1 +√
2)p 12v2
∥∥∥Q−1∥∥∥ /2
κef = (6 + 9√
2)p 12v2
∥∥∥Q−1∥∥∥ /4 + v2/6.
Os limites no erro obtidos aqui, para interpolação quadrática, podem sergeneralizados para interpolação de ordem superior.
3.2.2 Regressão de modelos não lineares
Considerando novamente Pdn, o espaço de polinômios de grau menor ouigual a d em Rn. Seja q1 = q + 1 a dimensão deste espaço, e seja φ =φ0(x), φ1(x), . . . , φq(x) uma base para Pdn. Dada uma base polinomial φ, define-se φ(x) = [φ0(x), φ1(x), . . . , φq(x)]T sendo um vetor em Rq1 . Como antes, Y =y0,y1, . . . ,yp ⊂ Rn denota o conjunto de p1 = p + 1 pontos amostrados. Sejaf(Y ) o vetor com elementos dados por f(yi), i = 0, . . . , p. Será considerado ocontexto onde a dimensão q1 do espaço polinomial é fixa mas o número p1 de pontosamostrados pode mudar.
É bem sabido que algumas das aplicações de optimização sem derivadas contémfunções com ruído. A maioria dos métodos sem derivadas são baseados naamostragem da função objetivo em vários pontos, suficientemente espaciados. Esseesquema permite que eles sejam razoavelmente tolerantes ao ruido. No entanto,pode se preferir levar em consideração a presença de ruído quando se projeta umalgoritmo.
A maneira mais natural de manipular dados com ruido é substituir a interpolaçãoda função objetivo pela regressão por mínimos quadrados. Neste caso, as condiçõesde interpolação descritas previamente,
M(φ, Y )α = f(Y ),
21
são resolvidos no sentido de mínimos quadrados, significando que a solução α
procurada é tal que ‖M(φ, Y )α− f(Y )‖ é minimizada.Se as avaliações da função são relativamente baratas, mas ainda com ruído, então
pode ser efetivo amostrar a função objetivo em mais pontos locais (i.e., em pontospróximos ao “centro” do modelo) do que os necessários para interpolação. Nessecaso, p1 > q1, o que significa que o número de linhas número de linhas na matrizM(φ, Y ) é maior do que o número de colunas e assim o sistema de interpolação ésobre determinado. Na medida em que o número de pontos amostrados aumenta,a solução do problema com ruído por mínimos quadrados converge (sob algumascondições) à solução por mínimos quadrados da função verdadeira subjacente.
Especificamente, seja a função com ruído f(x) = fsuave(x) + ε, onde fsuave(x) é afunção verdadeira, e ε uma variável aleatória, independente de x e obtida a partirde uma distribuição com média zero. Assumindo que a base polinomial φ é fixa, econsiderando um conjunto de amostras Y p com tamanho |Y p| = p+ 1 variável. Sejao vetor aleatório αp a solução por mínimos quadrados do sistema
M(φ, Y p)α = f(Y p)
e o vetor αpsuave a solução por mínimos quadrados do sistema
M(φ, Y p)α = fsuave(Y p).
Agora, assumindo que o tamanho do conjunto de amostras Y p está indo parainfinito e que as seguintes condições são válidas:
lim infp→∞
λmin
(1
p+ 1M(φ, Y p)TM(φ, Y p))> 0, (3.14)
onde λmin denota o menor autovalor de uma matriz. Com estas condições pode-seprovar o seguinte resultado de consistência [30].
Teorema 3.2.3.E(αp) = αp
suave ∀p ≥ q
elimp→∞
Var(αp −αpsuave) = 0.
A condição (3.14) imposta na sequência de conjuntos Y p significa que os dadossão amostrados numa maneira uniforme gerando conjuntos bem postos.
Quando as avaliações da função são caras, é quase impensável acumularlocalmente os pontos suficientes para superar o número de pontos necessário parauma interpolação quadrática, a saber (n+ 1)(n+ 2)/2. Neste caso tanto a regressão
22
por mínimos quadrados como a interpolação vão produzir resultados idênticos. Mas,mesmo assim, a regressão pode ser útil para relaxar as condições de interpolação.Se a função objetivo contém muito ruido, então o modelo de interpolação pode vira ser desnecessariamente “não suave”. Pode-se então usar regressão regularizada,para tentar equilibrar entre a otimização do ajuste por mínimos quadrados e a“suavização” da interpolação polinomial ao mesmo tempo.
A seguir a análises serão feitas para o caso quando p > q. Primeiro se discutirãoas propriedades dos modelos de mínimos quadrados e por último se apresentarão oslimites no erro.
Conceitos básicos em regressão por mínimos quadrados
Seja m(x) o polinômio de grau menor ou igual a d que aproxima uma dadafunção f(x) nos pontos em Y através de regressão por mínimos quadrados. Já queφ é uma base em Pdn, então m(x) = ∑q
k=0 αkφk(x), onde os escalares αk representamos coeficientes desconhecidos. Os coeficientes α podem ser determinados a partirdas condições de regressão por mínimos quadrados
M(φ, Y )α l.s.= f(Y ), (3.15)
o que equivale a resolver
minα‖M(φ, Y )α− f(Y )‖2 .
O sistema anterior têm solução única se a matriz
M(φ, Y ) =
α0(y0) α1(y0) · · · αq(y0)α0(y1) α1(y1) · · · αq(y1)
... ... ... ...
... ... ... ...α0(yp) α1(yp) · · · αq(yp)
(3.16)
tiver posto completo por colunas.É fácil ver que se M(φ, Y ) é quadrada e não singular, então o problema anterior
vira um problema de interpolação linear. Igual que no caso de interpolação, seM(φ, Y ) tem posto completo por colunas para uma escolha de φ, então esta tambémtêm posto completo por colunas para qualquer base de Pdn. Então, a seguintedefinição para conjunto bem posto para regressão polinomial por mínimos quadradosé independente da base escolhida.
23
Definição 3.2.2. O conjunto Y = y0,y1, . . . ,yp é bem posto para regressãopolinomial por mínimos quadrados em Rn se a matriz correspondente M(φ, Y ) têmposto completo por colunas para alguma base φ em Pdn.
Agora será mostrado que se Y é bem posto, então sua regressão polinomial pormínimos quadrados é única e independente da escolha de φ [29].
Lema 3.2.4. Dada uma função f : Rn → R e um conjunto bem posto Y ∈ Rn emrelação a regressão polinomial por mínimos quadrados, a regressão polinomial pormínimos quadrados m(x) existe e é única.
Omitimos a apresentação da demonstração já que é feita da mesma maneira àdo Lema 3.2.1.
Obtenção dos limites no erro em termos do número de condicionamentoSerão apresentados os limites no erro para o caso de regressão por mínimos
quadrados. Considerando a decomposição em valor singular de M = M(φ, Y ) =UΣVT (significando que U é uma matriz p1×q1 com colunas ortonormais, Σ é umamatriz q1 × q1 de valores singulares diferentes de zero, e V é uma matriz q1 × q1
ortonormal). Seja σ1 o maior valor singular de M e σq1 o menor. Então ‖M‖ =‖Σ‖ = σ1 e
∥∥∥Σ−1∥∥∥ = 1
σq1. Note que cond(M) = σ1
σq1. Podem se obter diretamente
limites no erro, usando Σ−1, como feito na Seção 3.2.1. Aqui se apresentam oslimites para o caso quadrático. Cabe notar que o caso linear já foi apresentado naSeção 3.1.2.
Suposição 3.2.2. Assume-se que Y = y0,y1, . . . ,yp ⊂ Rn, com p1 = p + 1 >
(n + 1)(n + 2)/2, é um conjunto de amostras bem posto (no sentido de regressãoquadrática, d = 2) contido na bola B(y0; ∆(Y )) de raio ∆ = ∆(Y ). Adicionalmente,assume-se que a função f é duas vezes continuamente diferenciável num domínioaberto Ω contendo B(y0; ∆(Y )) e que ∇2f é Lipschitz-continua em Ω com constantev2 > 0.
Teorema 3.2.5. Mantendo a Suposição 3.2.2. Então, para todos os pontos y emB(y0; ∆(Y )), tem-se que
? o erro entre a Hessiana do modelo de regressão quadrática e a Hessiana dafunção satisfaz ∥∥∥∇2f(y)−∇2m(y)
∥∥∥ ≤ κeh∆,
? o erro entre o gradiente do modelo de regressão quadrática e o gradiente dafunção satisfaz
‖∇f(y)−∇m(y)‖ ≤ κeg∆2,
24
? o erro entre o modelo de regressão quadrática e a função satisfaz
‖f(y)−m(y)‖ ≤ κef∆3,
onde κeh, κeg, e κef são dadas por
κeh = 3√
2p 12v2
∥∥∥Σ−1∥∥∥ /2
κeg = 3(1 +√
2)p 12v2
∥∥∥Σ−1∥∥∥ /2
κef = (6 + 9√
2)p 12v2
∥∥∥Σ−1∥∥∥ /4 + v2/6.
3.3 Modelos não polinomiais
Nas aplicações de modelagem em engenharia é frequente o caso em que afunção a ser otimizada é dispendiosa de ser avaliada. Para resolver o problema(aproximadamente) pode ser necessário o uso extensivo de simulação de sistemas deequações diferenciais, possivelmente associados com diferentes disciplinas, ou podeenvolver outros cômputos numéricos demorados. Na notação deste documento, afunção real é denotada por f(x), onde f : Rn → R.
O modelo de aproximação pode servir para vários propósitos; em particular estepode ser usado unicamente para modelamento e análise, a fim de obter informaçõessobre as características do problema e comportamento, sem avaliações muitas caras.No contexto de otimização, com frequência o uso dado aos modelos de aproximaçãoé de substituir a função real com o propósito de realizar otimização. Algunsdos métodos para otimização sem derivadas são baseados na seleção de modelosadequados, em alguns aspectos essencialmente simples, como modelos substitutosda função real. Um modelo substituto é tipicamente menos preciso ou têm menosqualidade do que a função verdadeira, mas é mais barato para avaliar ou consomepoucos recursos computacionais. Várias avaliações do modelo aproximado podemser ainda mais baratas do que uma avaliação da função real (modelo real).
A seguir uma classificação geral dos modelos precede uma breve apresentaçãodos principais modelos não polinomiais.
3.3.1 Classificação dos modelos aproximados
Desde um ponto de vista geral e simplificado os modelos de aproximação podemser classificados em dois tipos básicos: soluções numéricas das equações governantesdos sistemas físicos; e aproximações funcionais das soluções das equações construídassem recorrer ao conhecimento do sistema físico, ou seja, usando somente os valoresda função. Por conveniência denotaremos esses tipos como modelos físicos e
25
funcionais, respectivamente. A principal característica que os distingue é que oprimeiro incorpora conhecimento do sistema físico, enquanto o último é puramenteuma construção matemática. O modelo funcional só incorpora conhecimento docomportamento da função ao aproximar nos pontos para os quais se forneceramvalores, enquanto o modelo físico incorpora conhecimento do comportamento dosistema em todos os pontos.
A distinção entre esses dois tipos de modelos não é absoluta. Há algunsmodelos que podem ser definidos como pertencendo aos dois tipos. Por exemplo,os modelos em séries de Taylor são modelos funcionais já que a diferenciação dasequações governantes é um processo puramente matemático que pode ser feito semconhecimento algum do sistema. Pode-se dizer que este também é um modelo físicojá que as derivadas das equações governantes descrevem o comportamento do mesmosistema que as equações governantes descrevem.
Outra diferença é que para construir um modelo do tipo funcional tipicamente serequer de algum número de avaliações das funções reais, enquanto que um modelofísico pode ser usado sem tal custo inicial. Visto que estamos assumindo que emgeral as avaliações das funções reais são custosas, o preço para construir um modelofuncional pode ser um fator significativo na escolha de uma estratégia de modelagempara uma aplicação particular.
Dentro dos modelos funcionais, se distinguem dois tipos de modelos: modelosinterpolantes e modelos não interpolantes. Os modelos interpolantes passam pelosdados fornecidos enquanto os modelos não interpolantes não cumprem com isto masincorporam certo grau de flexibilidade (se incorporar mais dados) na aproximação.
A Figura 3.1 ilustra a classificação descrita.
Modelos Físicos, fenomenológicos, ou Conduzidos por
Modelos
CAIXA BRANCA Ex.: Equações de
Euler
Modelos Físicos e Funcionais
(Conduzidos por Modelos e Dados)
CAIXA CINZA
Ex.: Expansão em Série de Taylor.
Modelos Funcionais ou
Conduzidos por Dados
CAIXA PRETA
Ex.: Superfícies de Resposta, Redes Neurais Artificiais
MODELOS INTERPOLANTES
Ex: - Polinômios de Lagrange - Kriging
MODELOS NÃO INTERPOLANTES
Ex: - Ajustes por Mínimos Quadrados - Splines
Figura 3.1: Classificação geral dos modelos de aproximação
Exemplos de modelos do tipo de equações governantes são o uso de fluxopotencial ou equações de Euler como modelos aproximados das equações de Navier-Stokes. A fidelidade de um modelo físico pode variar dependendo de algo tão simplescomo um parâmetro de espaçamento da malha, ou tão complexa como a formulaçãode um modelo de turbulência num solucionador de equações de Navier-Stokes.
26
Exemplos de modelos do tipo funcional são interpolações polinomiais, splines,redes neurais e outros tipos de ajustes de curvas, particularmente os modelos desuperfícies de resposta, e modelos de mínimos quadrados.
Ao longo desta tese, o termo “modelo” será usado para se referir a qualquerfunção que seja usada como uma aproximação de outra função com o propósito dereduzir o custo de avaliação ou de melhorar o comportamento da função. Ainda quea discussão se concentrará na modelagem da função objetivo, também se aplica àmodelagem das restrições.
3.3.2 Funções de base radial
Até aqui, foi tratado com detalhe somente modelos polinomiais (linear equadrático) de interpolação e regressão. Agora se abordarão alguns detalhesrelacionados a uma classe de modelos dados por funções de base radial, os quaistêm sido usados com frequência em otimização sem derivadas [31, 32].
Uma das características atraentes das funções de base radial é sua habilidade paramodelar bem a curvatura da função subjacente. Outra característica importante éque a matriz de coeficientes definida pelas condições de interpolação é não singularsob condições relativamente fracas. De fato, pode-se garantir limites úteis naestimação da norma e número de condição da matriz de interpolação (ver, e.g.,[33]). Os sistemas lineares de interpolação obtidos ao modelar com funções de baseradial são tipicamente densos. No entanto, isso é relevante unicamente quando háum número elevado de pontos de interpolação.
Para interpolar uma função f cujos valores em um conjunto Y = y0, . . . ,yp ⊂Rn são conhecidos, pode-se usar um modelo funcional de base radial na forma
mrb(x) =p∑i=0
λiφ(∥∥∥x− yi
∥∥∥), (3.17)
onde φ : R+ → R e λ0, . . . , λp ∈ R. O termo base radial vem do fato que φ(‖x‖)é constante em qualquer esfera centrada na origem em Rn. Para que mrb(x) sejacontinuamente diferenciável duas vezes, a função de base radial φ(x) além de sercontinuamente diferenciável duas vezes deve ter derivada nula na origem.
Algumas das funções de base radial mais populares são as seguintes:
? cúbica φ(r) = r3,
? Gaussiana φ(r) = e− r
2ρ2 ,
? multi-quadrática φ(r) = (r2 + ρ2) 32 ,
? inversa multi-quadrática φ(r) = (r2 + ρ2)− 12 ,
27
onde ρ2 é qualquer constante positiva.Já que em muitas aplicações é desejável que o espaço linear gerado pela função
de base radial inclua funções constantes (e as bases radiais padrão não fornecemesta propriedade quando n é finito), resulta útil aumentar a função de base radialem (3.17) adicionando um termo constante. Similarmente, se for desejável incluirfunções lineares (e/ou outras funções polinomiais de baixa ordem), adiciona-se umpolinômio de baixa ordem de grau d−1 que pode ser expressado como ∑q
j=0 γjpj(x),onde pj ∈ Pd−1
n , j = 0, . . . , q, são as funções base para o polinômio e γ0, . . . , γq ∈ Rn
Com isto o novo modelo de interpolação fica da seguinte forma:
mrb(x) =p∑i=0
λiφ(∥∥∥x− yi
∥∥∥) +q∑j=0
γjpj(x). (3.18)
Além disso, os coeficientes λi precisam satisfazer
p∑i=0
λipj(yi) = 0, j = 0, . . . , q.
Esta ultima expressão junto com as condições de interpolação mrb(yi) = f(yi), i =0, . . . , p, resultam no sistema linear Φ P
PT 0
λ
γ
=f(Y )
0
, (3.19)
onde Φij = φ(‖yi − yj‖) e Pij = pj(yi) para i, j ∈ 0, . . . , p e f(Y ) é o vetorformado pelos valores f(y0), . . . , f(yp).
A matriz simétrica (3.19) tem solução única para as bases φ fornecidas acima,desde que P tenha posto completo e d ≥ 2.
A abordagem de OEUVRAY [31] para otimização sem derivadas usa funções debase radial cúbicas e polinômios lineares:
mbr(x) =p∑i=0
λi∥∥∥x− yi
∥∥∥3+ c+ gTx. (3.20)
Neste caso, se o número de pontos de interpolação é p + 1, então o modelo temp + n + 2 parâmetros, p + 1 para os termos da base radial e n + 1 para os termosdo polinômio linear. No entanto, quando o número de pontos é n + 1 (ou menos),a solução do sistema de interpolação fornece um polinômio linear, já que todos osparâmetros λi, i = 0, . . . , p, são nulos. Consequentemente, o modelo mais simplesnão linear mrb(x) da forma (3.20) é baseado em n+ 2 pontos de interpolação e têm2n+ 3 parâmetros.
28
Pode-se mostrar que o modelo (3.20), construído num conjunto de amostras bemposto contido numa bola de raio ∆, provê um erro nos valores da função da ordemde ∆2 e no gradiente da ordem de ∆, como acontece para interpolação ou regressãolineares.
3.3.3 Modelos de Kriging
Os modelos funcionais podem incorporar um componente estocástico, o que podetorná-los mais adequados para fins de otimização global. Um exemplo popular éKriging. Falando em termos gerais, os modelos de Kriging estão formados por doiscomponentes. Tipicamente o primeiro componente é um modelo simples destinadoa captar a tendência dos dados. O outro componente mede o desvio entre o modelosimples e a função real.
Para facilitar a explicação, no exemplo seguinte o modelo simples escolhido éuma constante. Assume-se que a função real é da forma
f(x) = β + z(x),
onde z(x) segue uma distribuição Gaussiana com média 0 e variância σ2. Acovariância entre dois pontos amostrais y e w é modelada por:
R(y,w) = Cov(z(y), z(w)),
onde R(y,w) é tal que R(Y, Y ) é definida positiva para qualquer conjunto Y depontos amostrais distintos. Uma escolha muito comum é modelar a covariança porfunções de base radial R(y,w) = φ(‖y−w‖).
Dado um conjunto Y de pontos amostrais diferentes, o modelo de Krigingmkr(x)é definido como o valor esperado da função real dados os valores observados f(Y ):
mkr(x) = E(f(x)|Y ).
Pode-se provar que mkr(x) é da forma
mkr(x) = β + R(Y,x)TR(Y, Y )−1(f(Y )− βe), (3.21)
onde e é um vetor de uns e
β = eTR(Y, Y )−1f(Y )eTR(Y, Y )−1e
.
Pode-se interpretar a β como uma aproximação de β e aR(Y,x)TR(Y, Y )−1(f(Y )− βe) como uma aproximação média de z(x).
29
Os modelos de Kriging foram desenvolvidos originalmente na comunidade degeo-estatística (ver MATHERON [34]) e usados pela primeira vez por CURRINet. al. [35] para experimentos computacionais. Um artigo com uma boa revisão,referenciado por muitos pesquisadores é [5].
3.4 Amostragem
A construção de um modelo aproximado requer da obtenção de informaçãoda função f a ser aproximada. O processo comum consiste na obtenção depontos (amostras) onde se avalia f para depois obter o modelo aproximado sejapor interpolação ou por ajuste. Conhece-se como amostragem ao processo queseleciona os pontos. Na amostragem se precisa selecionar cuidadosamente os pontospara cumprir com os critérios requeridos pelo modelo, e.g., quantidade de pontosrequeridos e regularidade na distribuição (conjunto bem posto). As estrategias deamostragem geralmente visam a distribuição uniforme dos pontos com a finalidadede que eles obtenham a melhor caracterização da função na região de aproximação.
Os métodos mais conhecidos para amostragem são o projeto de experimentos(DOE, Design Of Experiments) e a construção de padrões regulares ou simplex.
Método SIMPLEX para geração de um novo conjunto bem posto depontos amostrais
Definição 3.4.1. Um simplex S em Rn é o envoltório convexo de n + 1 pontos,yjnj=0, onde yj representa cada um dos vértices de S.
Neste método, todo novo conjunto bem posto de pontos amostrais Y =y0,y1, . . . ,yn ⊂ B(x(k); ∆(k)) é um simplex construído segundo o seguinteprocedimento:
? Primeiro ponto: y0 = x(k). i.e., o ponto base é o resultado da iteração atual.
? Do ponto y1 ao yn:
yi = x(k) ±∆(k)ei, i = 1, . . . , n, (3.22)
com o sinal e ∆ tais que yi ∈ Ωbc :=x : xl < x < xu
onde ei denota o i-ésimo vetor canônico.
30
3.4.1 Projeto de experimentos
A seleção de pontos no espaço de projeto onde a resposta tem que seravaliada é freqüentemente chamada projeto de experimentos. Assim, um projetode experimentos representa a seqüencia de experimentos (análises) a ser realizada,expressada em termos de fatores (variáveis de projeto) fixos em níveis especificados(valores predefinidos). Um projeto experimental é comumente representado pelamatriz X onde as linhas denotam execuções de experimentos, e as colunas denotamas configurações de um dado fator. Assim, X pode ser representada como
X =
x
(1)1 x
(1)2 · · · x(1)
n... ... . . . ...
x(ne)1 x
(ne)2 · · · x(ne)
n
. (3.23)
Na maioria das aplicações onde assume-se que não é conhecido o comportamentoda função, e os dados são colecionados a partir de simulações computacionaisdeterminísticas, o principal interesse é minimizar o erro de desvio (bias error) já queo erro aleatório é pequeno. Para ter um ponto de vista estatístico de comparaçãodos diferentes projetos de experimentos explica-se a seguir o significado dos termosdesvio e variância.desvio (bias): quantifica a magnitude em que a saída do modelo aproximado(resposta), ou seja f(x) difere do valor real f(x). O desvio se calcula como a médiade todos os conjuntos de dados possíveis D.variância (var): mede quão sensível é o modelo aproximado f(x) ao conjuntoparticular de dados D. Cada conjunto de dados D corresponde a uma amostraaleatória da função real.Para uma formulação de erro médio as expressões para desvio e variância se mostramnas Eqs. (3.24) e (3.25) respectivamente. Nas duas expressões ETCD denota o valoresperado considerando todos os conjuntos de dados possíveis.
Ebias2(x) = ETCD[f(x)]− f(x)2, (3.24)
Evar(x) = ETCD[f(x)− ETCD[f(x)]]2. (3.25)
Existe uma relação entre o desvio e a variância. Um modelo aproximado que ajustaestreitamente um conjunto particular de dados (pouco desvio) tenderá a prover umagrande variância. Podemos diminuir a variância ao suavizar o modelo aproximadomas ao simplificar demais o modelo o erro de desvio fica significativamente grande.Em princípio, podemos reduzir tanto o desvio (pode-se escolher modelos maiscomplexos) como a variância (cada modelo mais fortemente restringido pelos dados)
31
ao aumentar o número de pontos de maneira que o aumento na restrição feita pelosdados seja maior que o aumento na complexidade do modelo.Na prática, o número de pontos no conjunto de dados é fortemente limitado devidoao custo computacional, por isto durante a construção do modelo aproximadofreqüentemente procura-se um equilíbrio entre os erros originados por desvio e oscausados por variância. Este equilíbrio pode ser atingido reduzindo o erro de desvioenquanto se impõem penalidades na complexidade do modelo.
A seguir mostra-se uma breve explicação dos projetos de experimentos maisusados nas técnicas de otimização baseadas em metamodelos.
Projeto fatorial completoAntes de criar um projeto de experimentos, define-se o alcance permitido para cadauma das n variáveis mediante limites inferior e superior. Depois se discretiza oalcance em níveis igualmente espaçados. Para conseguir estabilidade numérica epara facilitar a notação escala-se o alcance de cada variável para cobrir [−1, 1] [36].A região confinada pelos limites inferior e superior nas variáveis é denominada espaçode projeto, e seus vértices determinam um cubo n-dimensional. Assim um projeto2n fatorial completo é criado se cada uma das variáveis é especificada somente noslimites inferior e superior (dois níveis). Similarmente um projeto 3n fatorial completoé criado ao especificar os limites inferior, intermédio e superior (três níveis: [ -1, 0,1 ]) para cada uma das n variáveis. Na Figura 3.2 se ilustra um projeto 33 fatorialcompleto.
32
Figura 3.2: Projeto 33 fatorial completo (27 pontos).
O projeto fatorial completo permite a estimação dos efeitos principais(x1, x2, . . . , xn) e das interações por pares dos fatores (xixj).
Os projetos 2n e 3n fatoriais completos tornam-se impraticáveis na medida emque aumenta o número de variáveis de projeto (n > 10). Por exemplo para n = 29usando o projeto 2n fatorial completo precisaríamos de 229 ≈ 5.4 · 108 avaliações(análises). Por isto, o projeto fatorial completo é usado tipicamente para dez oumenos variáveis. Para mais de dez variáveis podem ser usados outros projetos deexperimentos. Entre eles se encontram o projeto fracionário fatorial e o projetocomposto pequeno [36].
Projeto composto centralPara reduzir o número de experimentos necessários pode ser usado o projeto deexperimentos conhecido como o projeto composto central (CCD, Central CompositeDesign) [36]. O CCD é formado por um projeto 2n fatorial completo com 2n pontos‘estrela’ e um ponto adicional ‘centro’. Na Figura 3.3 ilustra-se um projeto CCDde três variáveis. Neste projeto de experimentos os pontos ‘estrela’ ficam por forado contorno criado pelo experimento 2n fatorial completo. A distância desde ospontos estrela ao centro do CCD tipicamente variam de 1.0 a
√n [36]. O projeto
CCD proporciona 2n + 2n + 1 experimentos. Como no casso do projeto fatorialcompleto, o número de experimentos requerido pelo projeto CCD também torna-se
33
impraticável quando n é muito grande (n > 10).
- Pontos fatoriais ∗∗∗∗ - Pontos estrela - Ponto central
Figura 3.3: Projeto composto central de três variáveis (15 pontos).
Projeto D-ótimoOs projetos de experimentos do tipo factorial completo e composto central foramdesenvolvidos para serem usados com espaços de projeto retangulares. No entanto,neste trabalho trata-se com problemas não lineares e, por tanto, com espaços deprojeto de forma irregular. Um projeto de experimentos do tipo D-ótimo forneceum método atrativo para criar projetos de experimentos dentro de um espaço deprojeto com forma irregular [36]. Além disso, o projeto de experimentos D-ótimopode ser construído com menos experimentos (avaliações) do que o requerido peloprojeto composto central. Na Figura 3.4 mostra-se um exemplo de projeto D-ótimo.
34
Figura 3.4: Projeto D-ótimo de três variáveis (11 pontos).
Um projeto de experimentos D-ótimo é uma coleção de lugares de amostragemno espaço de projeto para os quais se maximiza a quantidade |XTX| sobre todasas localizações possíveis. Esta escolha minimiza a variância das estimações dosparâmetros. As propriedades mais relevantes do projeto de experimentos D-ótimosão:
1. Minimiza a variância (incerteza) nos coeficientes estimados.
2. Minimiza a máxima variância de qualquer valor previsto y(x).
3. A quantidade |XTX| é invariante ao escalado de x.
Para criar um projeto de experimentos D-ótimo usam-se métodos de otimização.Enquanto o problema de otimização pode ser representado tanto de forma contínuaquanto em forma discreta, prefere-se a formulação discreta quando se consideramespaços de projeto que tenham contornos com forma irregular para os quais não sedisponha de descrição analítica. Na formulação discreta, selecionam-se um conjuntode ne experimentos no espaço de projeto a partir de um grupo de nc candidatos (nc ≥ne), onde o grupo dos nc candidatos é definido a priori pelo usuário. Assim, paraachar os ne experimentos que maximizam |XTX|,emprega-se um método iterativode otimização.
35
A seleção de ne experimentos de um grupo de nc candidatos é um problemacombinatório que tem
(ncne
)= nc!/(ne!(nc − ne)!) combinações possíveis. Este
problema combinatório cresce rapidamente como se ilustra num exemplo queenvolve duas variáveis (n=2), onde o espaço de projeto bi-dimensional se discretizausando uma malha de 11 × 11 (nc = 121), para obter 25 experimentos (ne = 25).Então neste caso o problema combinatório é
(12125
)para o qual existem 5 · 26 · 1025
combinações possíveis, onde um ou mais experimentos são D-ótimos. Para este casoos métodos mais usados de otimização combinatória são baseados em algoritmosgenéticos (GA, Genetic Algorithms) e por limitações em custo computacional sãousados em projeto de experimentos do tipo D-ótimo com até 25 variáveis, mas serecomenda seu uso para problemas com até cinco variáveis.
Projeto de preenchimento do espaçoMuitos pesquisadores defendem o uso de projetos de preenchimento do espaço(space filling) para amostrar experimentos computacionais determinísticos já queeles tratam uniformemente todas as regiões do espaço de projeto [37]. A propriedadede uniformidade nos projetos de experimentos é procurada ao maximizar a distanciamínima entre os pontos de projeto, ou minimizando medidas de correlação dentrodos dados amostrais. Um exemplo de implementação prática destas estratégias é aamostragem do hiper-cubo de Latin (LHS, Latin Hypercube sampling).Uma amostragem estratificada garante que todas as porções de uma dada regiãosejam amostradas. LHS [38] é uma estratégia de amostragem estratificada coma restrição de que cada uma das variáveis de entrada (xi) têm todas as porçõesda sua distribuição representada pelos valores de entrada. Seguindo a notação deMcKay et al., uma amostra de tamanho ne pode-se construir ao dividir o alcancede cada variável de entrada em ne estratos de igual probabilidade marginal e pegaruma amostra de cada estrato. Denotamos esta amostra por x(p)
i , i = 1, . . . , n; p =1, . . . , ne. A amostra é feita de componentes de cada uma das (xi), combinadasaleatoriamente. A Figura 3.5 ilustra um projeto LHS para obter 5 experimentos(ne = 5) a partir de duas variáveis de projeto (n = 2).
36
0 0.2 0.4 0.6 0.8 10
0.2
0.4
0.6
0.8
1
x1
x2
Figura 3.5: Um projeto LHS com n = 2, ne = 5 para x = (x1, x2) uniformementedistribuída num quadrado unitário.
37
Capítulo 4
Estimação das derivadas e controledo erro
Neste capítulo se apresentam as estratégias que vão ser usadas para o cálculo dasderivadas de maneira aproximada e para estimar e controlar seu correspondente erro.A maioria destas estratégias são baseadas nos conceitos apresentados no capítulo 3.
4.1 Cálculo das derivadas aproximadas
As derivadas da função objetivo (∇f(x)) e das restrições (∇gj(x), j = 1, . . . ,m)serão obtidas através de modelos polinomiais (f(x), gj(x), j = 1, . . . ,m). Se assumeque avaliar as funções têm um custo computacional alto e por isto procura-se usara menor quantidade de avaliações possível.
A obtenção das derivadas envolve a criação de um conjunto de pontos Y
através de amostragem e posteriormente a solução de um sistema linear para ter oscoeficientes do polinômio que formam os componentes do gradiente correspondente.O procedimento é explicado a seguir para o caso de interpolação linear e tambémpara o caso de mínimos quadrados ponderados.
4.1.1 Estimativa das derivadas usando interpolação linear
Primeiro se cria um conjunto de pontos de interpolação Y bem posto (nosentido de interpolação linear) na bola B(y0; ∆), centrada em y0, de radio ∆ =∆(Y ) = max
1≤i≤p‖yi − y0‖, através de um dos procedimentos explicados na Seção 3.4,
especificamente o que cria um simplex a partir de n+1 pontos. O outro caso é quandoo conjunto já está criado, nesse caso o conjunto é atualizado pela substituição de umdos seus pontos ou pode ser totalmente recriado para diminuir o erro da aproximaçãocomo será explicado com mais detalhe na Seção 4.2.
38
Seja cada função fj (função objetivo em restrições de desigualdade) representadapela sua correspondente aproximação linear fj(x) que interpola fj nos pontos emY = y0,y1, . . . ,yp, i.e., satisfazendo as condições de interpolação
fj(yi) = fj(yi), j = 1, . . . ,m+ 1, i = 0, . . . , n. (4.1)
Pode-se expressar fj(x) na forma
fj(x) = αj0 + αj1x1 + · · ·+ αjnxn, j = 1, . . . ,m+ 1
Assim, pode se reescrever (4.1) como
1 y01 · · · y0
n
1 y11 · · · y1
n... ... ... ...1 yn1 · · · ynn
αj0
αj1...αjn
=
fj(y0)fj(y1)
...fj(yn)
. (4.2)
Sendo o sistema (4.2) representado por Mαj = fj, da Seção 3.1.2 temos que oselementos de ∇fj(x) são proporcionados pela solução do sistema
y1
1 − y01 · · · y1
n − y0n
... ... ...yn1 − y0
1 · · · ynn − y0n
αj1...αjn
=
fj(y1)− fj(y0)
...fj(yn)− fj(y0)
, (4.3)
representado por M0∇fj = δ(fj, Y ) Ou seja,
∇fj(x) = [αj1 αj2 . . . αjn]T , j = 1, . . . ,m+ 1. (4.4)
4.1.2 Estimativa das derivadas usando regressão linear
Nesta seção se explica uma estratégia proposta para a obtenção das derivadasaproveitando toda a informação disponível em relação às avaliações das funções reais(i.e., usando todas as avaliações feitas nos pontos obtidos ao longo do processo deotimização).
Para isto primeiro se cria um conjunto de pontos de interpolação Y bem posto(no sentido de regressão linear) na bola B(y0; ∆), centrada em y0, de radio ∆ =∆(Y ) = max
1≤i≤p‖yi − y0‖, através de um dos procedimentos explicados na Seção 3.4,
especificamente o que realiza um projeto de experimentos (DOE) para obter q pontos(q > n + 1). Esta amostragem obtida pelo DOE escolhido é feita apenas uma vez.Após, o conjunto criado é atualizado pela adição de novos pontos obtidos ao longodo processo de otimização.
39
Seja cada função fj (função objetivo em restrições de desigualdade) representadapela sua correspondente aproximação linear fj(x) que ajusta fj aos dados em Y pormínimos quadrados.
Sejam os pontos de amostragem Y = y0,y1, . . . ,yp, p ≥ n representados comolinhas dentro da matriz M do sistema de regressão linear Mαj = fj dado por
1 y0
1 · · · y0n
1 y11 · · · y1
n... ... ... ...1 yp1 · · · ypn
αj0
αj1...αjp
=
fj(y0)fj(y1)
...fj(yp)
, (4.5)
onde a matriz M é comum a todas as funções (função objetivo e restrições).Um quesito importante que deve ser imposto nos modelos aproximados fj é o
de consistência. Este quesito é definido na Eq. (4.6), e determina que as funçõesaproximadas precisam ter o mesmo valor das funções originais quando avaliadas noponto de iteração atual xk.
fj(xk) = fj(xk), j = 1, . . . ,m+ 1. (4.6)
A maneira mais simples de atingir a condição de consistência (4.6) é através deuma mudança de variáveis. Após a mudança a nova matriz representada por χ édada por (4.7).
χ = M− e (xk)T. (4.7)
Da mesma maneira, o vetor de resposta após a mudança é dado por
yyyj = fj − f(xk) e. (4.8)
Para obter os parâmetros αj do modelo de aproximação usamos mínimosquadrados ponderados (WLS,Weighted Least Squares), dando maior peso aos pontospróximos ao ponto de iteração atual xk, deixando a aproximação localmente maisprecisa.
Seja W ∈ R(p+1)×(p+1) a matriz de ponderação, onde p+ 1 é o número de pontosenvolvidos na construção do modelo; i.e.,
W def= diagw1, w2, . . . , w(p+1).
Os parâmetros αj estimados por WLS são
αj = [χTWχ]−1χTWyyyj, (4.9)
40
O uso de WLS envolve armazenar p + 1 valores adicionais (i.e. o vetor deponderação w). É fácil ver que multiplicando cada linha de χ e cada elemento deyyyj por
√wi o problema de WLS pode ser resolvido através da solução do problema
convencional de mínimos quadrados (LS, Least Squares). Assim, definindo yyy∗j e χ∗
tal que [yyy∗j ]i :=√
[w]i [yyyj]i e [χ∗]il :=√
[w]i [χ]il, (i = 1, . . . , p + 1, l = 1, . . . , n)temos a seguinte solução para WLS.
αj = [χ∗Tχ∗]−1χ∗Tyyy∗j . (4.10)
Seja yi o i-ésimo ponto amostrado. O peso dado para cada ponto é fornecidopela função de ponderação wi : Rn → R definida pela spline cúbica
wi(s) def= 1− 3s2 + 2s3, (4.11)
onde, s = dist(yi,xk)/ ir, ir é um raio que define a região de influência daponderação (i.e., para incluir todos os pontos ir := maxdist(yi,xk)), e dist(yi,xk)é uma distância (neste caso: norma Euclidiana ‖•‖2) entre yi e xk, i = 1, . . . , p+ 1.
Como ilustrado na Figura 4.1, a função de ponderação proposta fornece umaponderação maior para pontos bem próximos ao ponto atual xk e uma ponderaçãomenor para pontos distantes de xk.
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 10
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
s
wi(s)
Figura 4.1: Função de ponderação wi para WLS.
A solução αj obtida de (4.9) ou de (4.10) corresponde diretamente aos elementosde ∇fj(x). Ou seja,
∇fj(x) = [αj1 αj2 . . . αjn]T , j = 1, . . . ,m+ 1. (4.12)
41
4.2 Controle do erro na aproximação dasderivadas
Um dos principais objetivos na obtenção das derivadas de maneira aproximada éfazer com que o erro em ditas aproximações seja diminuído sucessivamente atingindoum valor mínimo desejado após de um número finito de avaliações. Isto equivalea dizer que as aproximações sejam consistentes e convergentes. Para começarnossa análise é necessário retomar os resultados teóricos obtidos no Capítulo 3,que caracterizam o erro nas aproximações. Na Tabela 4.1 se resumem os resultados(Teoremas 3.1.1, 3.1.2 e 3.2.2) para aproximações polinomiais. Cabe destacar queas expressões para regressão são as mesmas que para o caso de interpolação, porisso são omitidas as primeiras para não redundar.
Tabela 4.1: Resumo dos resultados teóricos para o erro na interpolação polinomial.
Interpolação linear Interpolação quadrática
Erro no gradiente ‖∇f(y)−∇m(y)‖ ≤ κeg∆ ‖∇f(y)−∇m(y)‖ ≤ κ2eg∆2
Erro na função |f(y)−m(y)| ≤ κef∆2 |f(y)−m(y)| ≤ κ2ef∆3
Pontos necessários n+ 1 (n+ 1)(n+ 2)/2
Lembra-se que estes resultados se referem a aproximações do gradientee da função obtidas a partir de um conjunto de pontos de amostragemY = y0,y1, . . . ,yp ⊂ Rn contido na bola B(y0; ∆(Y )), centrada em y0, deradio ∆ = ∆(Y ) = max
1≤i≤p‖yi − y0‖.
Sendo M0 a matriz resultante de um passo de eliminação Gaussiana na matrizde interpolação M(φ, Y ), M0 = M0/∆. κeg e κef são dadas por
κeg = v(1 + n12∥∥∥(M0)−1
∥∥∥ /2),
κef = κeg + v
2 .
Sendo Q a matriz formada pelas últimas p linhas e colunas de M(φ, Y ), e Q suaversão escalada. κ2eg, e κ2ef são dadas por
κ2eg = 3(1 +√
2)p 12v2
∥∥∥Q−1∥∥∥ /2,
κ2ef = (6 + 9√
2)p 12v2
∥∥∥Q−1∥∥∥ /4 + v2/6.
42
Observa-se que as constantes κ que limitam o erro são de característicasemelhante e da mesma ordem de magnitude. O outro fator que limita o erro éo tamanho ∆ da região em que estão contidos os pontos.
Pode-se ver claramente que para um mesmo tipo de interpolação, por exemplopara o caso de interpolação linear, o erro no gradiente é maior do que o erro nafunção, para um conjunto de pontos que se encontrem muito próximos (∆ < 1).
Comparando a interpolação quadrática com a linear, observa-se que o errona interpolação quadrática é ∆ vezes o erro na interpolação linear. Embora ainterpolação quadrática tenha menos erro para pontos que estejam muito próximos(∆ < 1), a construção de um modelo quadrático requer da ordem de n2 avaliaçõesque são mais do que as requeridas para um modelo linear (da ordem de n). Comoo objetivo do novo algoritmo é apenas obter uma direção de busca aproximadacom o menor custo computacional possível (i.e., menor quantidade de avaliações dasfunções reais) se escolhe o uso de interpolações lineares para aproximar os gradientesda função objetivo e das restrições.
Feita esta escolha, no que segue, vamos trabalhar com as Eqs. (4.13) paracontrolar o erro na aproximação do gradiente obtido por interpolação linear.
‖∇f(y)−∇m(y)‖ ≤ κeg∆, (4.13a)κeg = v(1 + n
12∥∥∥(M0)−1
∥∥∥ /2), (4.13b)
∆ = ∆(Y ) = max1≤i≤n
‖yi − y0‖. (4.13c)
Das Eqs. (4.13) temos que dois fatores delimitam o máximo erro na aproximaçãodo gradiente. O primeiro fator (κeg) depende da norma
∥∥∥(M0)−1∥∥∥, que está
diretamente relacionada com o número de condição da matriz (M0), representandoassim a independência linear dos vetores yi − y0. Por isso vamos nos referir a esseprimeiro fator como fator geométrico. O segundo fator (∆) é a máxima normadesses vetores, que representa o tamanho da região que contém os pontos deinterpolação. A estratégia a ser utilizada no novo algoritmo para ir diminuindoo erro é manter uma adequada geometria nos pontos de interpolação enquantopropositalmente vamos diminuindo o tamanho da região de interpolação. Com istovamos fazer tender a zero o erro na aproximação dos gradientes, obtendo assimuma aproximação consistente e convergente aos valores reais.
Fator geométrico
Sabendo a importância do fator geométrico na delimitação do erro precisa-seobter uma estimativa dele. Isto equivale a estimar a norma da inversa de uma
43
matriz de tamanho n × n. O cálculo da inversa requer da ordem de n3 operações,que é um alto custo computacional. Outra maneira de calcular dita norma é atravésdo número de condição de M0 que é denotado por cond(M0) =
∥∥∥M0∥∥∥ ∥∥∥(M0)−1
∥∥∥,no entanto o procedimento padrão para obter ele requer fazer uma decomposiçãoem valor singular que também têm um custo computacional alto (O(n3)). A idéiaadotada neste trabalho é usar uma medida de baixo custo do fator geométrico, quepermita comparar diferentes configurações de pontos de amostragem com o objetivode obter o melhor conjunto possível (menor valor do fator geométrico) evitandoenvolver avaliações adicionais das funções originais. Para não prejudicar a eficiênciado novo algoritmo se precisa de algum procedimento para ter uma estimativa destefator, mesmo que seja aproximada.
Para estudar as estimativas do fator geométrico vamos visualizar cada uma nosbaseando em modificações da matriz identidade para conter um novo ponto x emR2 e R3. Isto é, usando as matrizes fornecidas pela Eq. (4.14).
M0 = 1 0x1 x2
parax ∈ R2, M0 =
1 0 00 1 0x1 x2 x3
parax ∈ R3. (4.14)
A Figura 4.2 ilustra estas duas matrizes, onde cada linha da matriz representaum vetor contendo a posição de um ponto do conjunto de pontos amostrais.
(x1,x
2)
(1,0)
x1
x2
(a) M0, x ∈ R2.
00.2
0.40.6
0.81
0
0.2
0.4
0.6
0.8
10
0.2
0.4
0.6
0.8
1
x 3
x1
x2
(x1,x
2,x
3)
(1,0,0)(0,1,0)
(b) M0, x ∈ R3.
Figura 4.2: Representação da matriz M0 dada pela Eq. (4.14).
Então, fazendo variar a posição do ponto x dentro da região delimitada na Figura4.2 podemos visualizar a magnitude de
∥∥∥(M0)−1∥∥∥, cond(M0) e suas estimativas
44
aproximadas correspondentes, através de curvas de nível (para x ∈ R2) e de iso-superfícies (para x ∈ R3). Ver Figuras 4.3, 4.4 e 4.5.
1.1
1.5
1.5
2
2
2
33
3
3
55
5
5
1010
10
100 100 100
1.1 1.5
1.5
2
2
22
33
3
3
5 55
5
10 10 10
100 100 100
x1
x 2
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 10
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Figura 4.3: Curvas de nível para cond(M0) (preto) e para∥∥∥(M0)−1
∥∥∥ (cinza).
← κ = 1.1
← κ = 1.5
← κ = 2 κ = 5 →
κ = 10 →
κ = 50 →
Figura 4.4: Iso-superfícies para o número de condição de M0, κ =∥∥∥M0
∥∥∥ ∥∥∥(M0)−1∥∥∥
45
← χ = 1.1
← χ = 1.5
← χ = 2χ = 5 →
χ = 10 →χ = 50 →
Figura 4.5: Iso-superfícies para χ =∥∥∥(M0)−1
∥∥∥Através das ilustrações apresentadas nas Figuras 4.3-4.5 fica evidente a afirmação
demonstrada no Capítulo 3, que estabelece uma relação diretamente proporcionalentre
∥∥∥(M0)−1∥∥∥ e cond(M0). Graças a isto, a partir deste momento, pode-se
estimar o fator geométrico através de qualquer uma destas duas quantidades,indistintamente.
Numa primeira tentativa de quantificar o fator geométrico de maneiraaproximada poderia se pensar na máxima colinearidade entre qualquer par depontos, como uma medida da dependência linear das linhas de M0. Isto é: sejaV = y1−y0, . . . ,yn−y0 o conjunto das linhas da matriz M0, então uma medidaaproximada de cond(M0) estaria dada por
cond(M0) ≈ cn maxvi,vj∈V,i 6=j
∣∣∣viTvj∣∣∣
‖vi‖‖vj‖,
onde cn é uma constante que depende apenas de n.No entanto a colinearidade não é o único quesito a se levar em conta na
independência linear. De fato, um conjunto de vetores é linearmente independentese eles não são colineares e se nenhum deles pode ser obtido a partir da combinaçãolinear dos outros. A Figuras 4.6 e 4.7 ilustram claramente que o segundo quesito éimportante e não pode ser omitido (comparar com as Figuras 4.3 e 4.4).
46
0.01
0.01
0.1
0.1
0.1
0.2
0.2
0.2
0.3
0.3
0.3
0.5
0.5
0.5
0.7
0.7
0.7
0.7
0.9
0.9
0.9
0.9
0.990.99
0.99
x1
x 2
0 0.1 0,2 0.3 0,4 0,5 0.6 0.7 0,8 0.9 10
0.1
0,2
0,3
0,4
0,5
0,6
0,7
0,8
0,9
1
Figura 4.6: Curvas de nível para |v1Tv2|
‖v1‖‖v2‖ .
← ς = 0.05← ς = 0.15← ς = 0.25
← ς = 0.35
ς = 0.55 →
ς = 0.8 →
ς = 0.95 →
Figura 4.7: Iso-superfícies para ς: Cosseno Máximo entre as linhas da matriz M0.
Por outro lado, das propriedades das normas, sabe-se que se z é a solução paraAz = y, então
‖z‖‖y‖≤∥∥∥A−1
∥∥∥ ,
47
e o limite é atingido para algum vetor y escolhido segundo algum critério deotimalidade. Pretende-se então escolher um vetor y de maneira que a proporção‖z‖ / ‖y‖ seja tão grande quanto for possível e portanto uma estimativa razoávelpara ‖A−1‖ . Encontrar o y ótimo pode ter um custo proibitivo, mas pode se obteruma aproximação útil a um custo baixo. Uma abordagem (ver [39]) é escolher ycomo a solução do sistema ATy = c, onde c é um vetor com componentes ±1,com os sinais alternados para que o y resultante seja o maior possível. O custocomputacional desta metodologia é o da solução dos dois sistemas lineares, i.e., daordem de n2 operações.
A motivação para esta aproximação não parece ser obvia, mas basicamenteé equivalente a um passo na iteração inversa para o cálculo do vetor singularcorrespondente ao menor valor singular de A.
As Figuras 4.8 e 4.9 ilustram os resultados deste último procedimento para obteruma estimativa aproximada de
∥∥∥(M0)−1∥∥∥. Ao se compararem com os resultados
exatos das Figuras 4.3 e 4.5 se observa que a estimativa é boa o bastante para nossosfins. Por tanto, esta estrategia é a escolhida para ser usada no novo algoritmo parapoder quantificar o fator geométrico.
1.1
1.5
1.5
1.5
2
2
22
3 33
3
5 55
5
10 10 10
100 100 100
x1
x 2
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 10
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Figura 4.8: Curvas de nível para∥∥∥(M0)−1
∥∥∥estim.
.
48
← χ = 1.1
← χ = 1.5
← χ = 2χ = 5 →χ = 10 →χ = 50 →
Figura 4.9: Iso-superfícies para χ :=∥∥∥(M0)−1
∥∥∥estim.
49
Capítulo 5
Algoritmo de pontos interiores edireções viáveis com aproximaçãodas derivadas
5.1 Conceitos preliminares
O algoritmo desenvolvido neste trabalho resolve um problema de otimização nãolinear com restrições. Ele é baseado numa técnica de programação matemáticae usa modelos aproximados do gradiente. A seguir descreve-se matematicamenteo problema de otimização não linear e algumas definições importantes para oentendimento do algoritmo proposto neste trabalho (FDIPA-AG).
5.1.1 Problema de otimização não linear
Um problema de otimização consiste na minimização ou maximização deuma determinada função (função objetivo) sujeito a restrições nas suas variáveis.Matematicamente o problema geral de otimização não linear com restrições podeser formulado na seguinte forma:
minimizarx∈Rn
f(x),
sujeito a: gi(x) ≤ 0, i = 1, . . . ,m.(5.1)
A formulação (5.1) tem a seguinte interpretação: procurar na região viável
Ω := x ∈ Rn : gi(x) ≤ 0, i = 1, . . . ,m (5.2)
o vetor x∗ tal quef(x∗) ≤ f(x), para todo x ∈ Ω. (5.3)
50
O vetor x ∈ Rn é o vetor de variáveis do problema, ou variáveis de projeto,f : Rn → R é a função custo ou função objetivo. As funções gi : Rn → R sãodenominadas funções de restrição de desigualdade. As funções f e gi devem serfunções contínuas com derivadas contínuas.
Um ponto x∗ que verifica a condição (5.3) é denominado mínimo global doProblema (5.1). Infelizmente, a condição (5.3) não é fácil de verificar em geral. Omínimo global somente pode ser caracterizado em situações especiais, por exemplo,quando as funções f e gi são convexas.
As técnicas gerais de programação matemática procuram pontos denominadosmínimos locais.
Para o desenvolvimento do algoritmo proposto consideraram-se as seguintesdefinições relacionadas com o problema de otimização não linear.
5.1.2 Definições
Definição 5.1.1. Dado um ponto x ∈ Ω, o conjunto de restrições de desigualdadeativas é: A(x) := i : 1 ≤ i ≤ m, gi(x) = 0.
Definição 5.1.2. Um vetor d ∈ Rn é uma direção viável em x ∈ Ω se existe θ > 0tal que x + td ∈ Ω para todo t ∈ [0, θ].
Observações:1) Se x pertence ao interior de Ω, então toda direção d ∈ Rn é viável.2) Se i ∈ A(x) e dT∇gi(x) < 0, i = 1, . . . ,m, então d é uma direção viável em x.
Definição 5.1.3. Um vetor d ∈ Rn é uma direção de descida para f em x ∈ Rn seexiste δ > 0 tal que f(x + td) < f(x) para todo t ∈ (0, δ).
Observação:Se dT∇f(x) < 0, então d é uma direção de descida para f em x.
Definição 5.1.4. O ponto x∗ ∈ Ω é um mínimo local do Problema (5.1) se existeum ε > 0 tal que tal que f(x∗) ≤ f(x) para todo x ∈ Ω ∩ B(x∗, ε). Quandof(x∗) < f(x), x∗ é dito mínimo local estrito.
Se f e gj são da classe C1, f e gj possuem as n derivadas parciais
∂f
∂xi(x), ∂gj
∂xi(x), x ∈ Rn, j = 1, . . . ,m, i = 1, . . . , n
e estas derivadas parciais são funções contínuas. O gradiente de f é denotado por
∇f(x) =[∂f
∂x1(x) · · · ∂f
∂xn(x)
]T
∈ Rn. (5.4)
51
As derivadas parciais de gi em relação a x(k) se agrupam na seguinte matriz∇g(x) definida como
∇g(x) = [∇g1(x) · · · ∇gm(x)] ∈ Rn×m. (5.5)
Definição 5.1.5. Define-se a função de Lagrange (ou Lagrangiana) L : Rn×Rm → Rassociada ao Problema (5.1):
L(x,λ) = f(x) + g(x)Tλ = f(x) +m∑i=1
λigi(x), (5.6)
onde g(x)T = [g1(x), . . . , gm(x)].
Em vista de (5.5) e (5.6), o gradiente da função Lagrangeana em relação a xrescreve-se:
∇xL(x,λ) = ∇f(x) +∇g(x)λ
= ∇f(x) +m∑i=1
λi∇gi(x). (5.7)
Definição 5.1.6. A Hessiana da função Lagrangiana para o Problema (5.1) é afunção H : Rn × Rm → Rn×n definida por:
H(x,λ) = ∇2xxL(x,λ) = ∇2f(x) +
m∑i=1
λi∇2gi(x). (5.8)
Definição 5.1.7. Dado um ponto x ∈ Ω, o Requisito de Independência Linear dasRestrições (RILR), conhecido também como condição de regularidade, é satisfeitoem x se os gradientes das restrições ativas são linearmente independentes, ou seja,se o conjunto ∇gi(x) : i ∈ A(x) é linearmente independente. Para estes pontos oespaço tangente se expressa como:
T (x) :=d : ∇gT
i (x)d = 0 ∀i ∈ A(x). (5.9)
Neste trabalho será sempre assumido que a condição de regularidade definidaacima é válida. Outras condições mais fracas, como por exemplo o Requisito deMangasarian-Fromovitz [40], não serão consideradas aqui.
A Figura 5.1 exemplifica os conceitos de direção de descida viável e espaçotangente.
52
ΩΩΩΩ x g2=0 T(x)
∇g2(x) ∇f(x) d curvas de nível de f Figura 5.1: Direção de descida viável d no ponto x. Espaço tangente T (x)
5.1.3 Condições de otimalidade
A seguir estão descritas as condições necessárias e suficientes de otimalidade quecaracterizam um mínimo local do Problema (5.1).
Teorema 5.1.1. (Condições necessárias) Se x∗ ∈ Rn é um mínimo local e pontoregular das restrições do Problema (5.1), então existe um único vetor λ∗ ∈ Rm talque
∇xL(x∗,λ∗) = 0 (5.10)gi(x∗)λ∗i = 0, i = 1, . . . ,m (5.11)
gi(x∗) ≤ 0 (5.12)λ∗i ≥ 0 (5.13)
dTH(x∗,λ∗)d ≥ 0 ∀d ∈ T (x∗). (5.14)
Demonstração. Ver [41].
As condições (5.10)-(5.13) são ditas condições necessárias de primeira ordem deKarush-Kuhn-Tucker (KKT). O ponto (x∗,λ∗) que verificar estas condições seráchamado de ponto de KKT. A Figura 5.2 mostra uma interpretação geométricadestas condições.
53
ΩΩΩΩ curvas de nível de f g1=0
λ2*∇g2(x* ) ∇f(x* ) g3=0 g2=0
λ1*∇g1(x* ) -∇f(x* ) x* Figura 5.2: Condições de KKT
Teorema 5.1.2. (Condições suficientes) Seja x∗ ∈ Rn um ponto de KKT e sejaλ∗ ∈ Rm o vetor de multiplicadores de Lagrange correspondente. Seja H(x∗,λ) talque satisfaz:
dTH(x∗,λ∗)d > 0 ∀d ∈ T (x∗). (5.15)
Então x∗ é um mínimo local estrito do Problema (5.1).
Demonstração. Ver [41].
5.2 Algoritmo FDIPAO algoritmo de ponto interior e direções viáveis (FDIPA, Feasible Directions
Interior Point Algorithm) [1] é um método iterativo que serve para resolver oProblema (5.1). Este algoritmo requer um ponto inicial x(0) pertencente ao interiorda região viável Ω.
A partir do ponto inicial o algoritmo FDIPA gera uma seqüência x(k)k∈Ntotalmente contida no interior de Ω e que converge a um ponto de KKT do Problema(5.1). A seqüência gerada pelo algoritmo reduz monotonamente a função f(x).
A iteração do algoritmo FDIPA é dada pela regra x(k+1) = x(k) + t(k)d, ondex(k) é o ponto da iteração atual, d é a direção de busca e t(k) é o passo da iteraçãok definido por um procedimento chamado busca linear. A definição da direção debusca do algoritmo FDIPA é inspirada na aplicação do método de Newton no sistemanão linear dado pelas Eqs. (5.10)-(5.11). A iteração de Newton para este sistema de
54
equações é H(x(k),λ(k)) ∇g(x(k))Λ(k)∇g(x(k))T G(x(k))
x(k+1) − x(k)
λ(k+1) − λ(k)
= −∇xL(x(k),λ(k))
G(x(k))λ(k)
,onde H é a matriz hessiana definida na Eq. (5.8) , L é definida na Eq. (5.6), Λ(k) =diag(λ(k)) e G(x(k)) = diag(g(x(k))). Definindo d = x(k+1) − x(k), e desenvolvendoo sistema linear anterior obtém-se
H(x(k),λ(k)) ∇g(x(k))Λ(k)∇g(x(k))T G(x(k))
dλ(k+1)
= −∇f(x(k))
0
. (5.16)
Na referência [1] é demonstrado que a direção de Newton pode não seruniformemente viável para a região Ω nem de descida para a função f(x). Por estarazão a direção de busca do algoritmo FDIPA é diferente da fornecida pela Eq.(5.16). O pseudo-código do algoritmo FDIPA é o seguinte:
Algoritmo: Algoritmo de Ponto Interior e Direções Viáveis (FDIPA)Parâmetros: α ∈ (0, 1), η ∈ (0, 1), ϕ > 0, e ν ∈ (0, 1).Dados: x(0) ∈ Ω, λ(0) ∈ Rmpositivo, B(0) ∈ Rn×n simétrica e definida positiva.Passo 1: Teste de convergência.Passo 2: Cálculo da direção de busca.
Passo 2a: Calcule (d0,λ0) e (d1,λ1) como a solução do sistema: B(k) ∇g(x(k))Λ(k)∇g(x(k))T G(x(k))
d0 d1
λ0 λ1
= −∇f(x(k)) 0
0 λ(k)
. (5.17)
Passo 2b: Se dT1∇f(x(k)) > 0, faça
ρ = minϕ‖d0‖2
2; (α− 1)dT0∇f(xk)
dT1∇f(xk)
.
Caso contrario, façaρ = ϕ‖d0‖2
2.
Passo 2c: Calcule a direção de busca d e o vetor λ:
d = d0 + ρd1,
λ = λ0 + ρλ1.
55
Passo 3: Busca linear. Achar um cumprimento de passo t(k) como o primeironúmero da seqüência 1, ν, ν2, ν3, . . . que satisfaz:
f(x(k) + t(k)d) ≤ f(x(k)) + t(k)η∇f(x(k))d,
gi(x(k) + t(k)d) < 0 se λi ≥ 0, ou
gi(x(k) + t(k)d) ≤ gi(x(k)) se λi < 0, i = 1, . . . ,m.
Passo 4: Atualização.
Passo 4a: Obtenha o novo ponto como x(k+1) = x(k) + t(k)d, e defina novosvalores: λ(k+1) > 0, e B(k+1) definida positiva.
Passo 4b: Retorne ao Passo 1.
5.3 Algoritmo FDIPA-WLS-AGComo foi afirmado no Capítulo 2, dispõe-se de diversos tipos de algoritmos
para resolver um problema de otimização quando não temos as derivadas. Todoseles possuem algumas características úteis, no entanto eles apresentam grandesdesvantagens, que em parte inspiraram a realização deste trabalho para desenvolveruma nova técnica nesta área. As principais desvantagens desses métodos são:
? Os métodos de busca direta: requerem muitas avaliações da função objetivona medida em que vai se aproximando da solução, prejudicando assim suaeficiência.
? Os métodos baseados em modelos: em sua maioria usam modelos quadráticosou não polinomiais que são usados para substituir as funções reais e obteruma solução aproximada. A construção desse tipo de modelo requer muitasavaliações da função objetivo.
? Os métodos heurísticos: em sua maioria carecem de convergência teórica queé observada também na prática. Geralmente têm diversos parâmetros cujoajuste depende do tipo de problema. Estes algoritmos levam muitas iterações(avaliações da função objetivo) até chegarem a uma solução razoável.
Além das desvantagens anteriores, todos eles foram desenvolvidos para resolverproblemas de otimização sem restrições. A adaptação deles para o caso comrestrições não é direta, e geralmente requer a formulação de subproblemas deotimização envolvendo alguma função de penalidade. É bem sabido que o uso defunções de penalidade não garantem a viabilidade da solução em cada iteração, o que
56
acaba prejudicando seu uso numa grande quantidade de aplicações de otimização emengenharia, como é o caso de otimização estrutural. Observa-se também que nenhumdesses métodos proporciona uma trajetória eficiente e rápida como a proporcionadaatravés dos gradientes.
Para superar os problemas citados acima foi desenvolvida uma primeira versão donovo algoritmo apresentado neste trabalho. Esta versão denominada FDIPA-WLS-AG (FDIPA with Weighted Least Squares and Approximated Gradients) resolve umproblema de otimização não linear com restrições e é baseada no algoritmo FDIPA,proporcionando soluções viáveis em cada iteração onde a informação do gradienteé obtida através de modelos aproximados da função objetivo e das restrições.Particularmente os modelos escolhidos para esta versão são polinômios linearesnão interpolantes; que permitem uma aproximação do gradiente, o suficientementerazoável, a um baixo custo computacional. As aproximações são consistentes eprocuram aproveitar a maior quantidade de informação disponível ao longo doprocesso de otimização. As ideias do algoritmo FDIPA-WLS-AG foram apresentadasem [42], onde se observam resultados promissores em otimização de forma estrutural.
5.3.1 Descrição do algoritmo
O algoritmo FDIPA-WLS-AG por ser baseado na adaptação de um algoritmode direções viáveis que usa o gradiente procura manter as características deconvergência deste na resolução de um problema não linear com restrições. Amodificação consiste basicamente na obtenção dos gradientes através da construção eatualização de modelos lineares não interpolantes usando os dados do modelo iniciale todos os pontos x(k) que vão sendo obtidos após cada iteração do algoritmo. A idéiade usar os pontos que formam o histórico do processo de otimização tem o objetivode evitar fazer avaliações desnecessárias das funções originais (da função objetivoe restrições) já que no caso geral assumimos que são computacionalmente caras.Com isto, espera-se uma maior eficiência em relação com o custo computacional daobtenção das derivadas por diferenças finitas, já que este último precisa realizar navaliações de cada função por iteração enquanto esta versão do algoritmo propostoem princípio usa apenas uma avaliação de cada função por iteração que se somamaos dados das n + 1 avaliações feitas na primeira iteração. A necessidade deincorporar mais dados do que incógnitas no modelo de aproximação fez com quea obtenção das aproximações seja feita por ajuste (mínimos quadrados) e não porinterpolação. Sendo assim, em vez de aproximar as derivadas usando uma espéciede secante aproximada, estamos usando sua tendência. Em princípio as tendênciasconseguidas por ajustes de mais pontos do que os necessários são menos sensíveisao ruído externo como confirmado no Teorema 3.2.3. No entanto isto requer que os
57
pontos sejam amostrados de maneira uniforme conformando um conjunto bem postopara regressão. Esta última condição em princípio não é garantida nesta versão doalgoritmo, já que os pontos obtidos no final do processo de otimização tendem a ficarmuito próximos e alinhados. Para minorar este problema optou-se pela ponderaçãoda informação usando mínimos quadrados ponderados WLS e pelo uso de diferençasfinitas quando necessário.
O algoritmo começa criando os modelos de aproximação do gradiente. Naprimeira iteração (k = 0) a informação do gradiente é adquirida por diferençasfinitas no ponto atual xk, isto é com o propósito de obter uma direção precisa inicialque seja de decida e viável e para ter uma estimativa sobre o tamanho da região deprojeto, isto ultimo sendo útil para a construção do modelo aproximado na primeiraiteração (k = 1).
Assim, uma vez aplicado diferenças finitas no ponto atual xk temos a seguinteestimativa do tamanho da região de amostragem
Definição 5.3.1. Seja Ωsdef= x ∈ Rn | x ≥ xk − re, x ≤ xk + re, x ≥ xl, x ≤
xu, e = [1, . . . , 1]T a região de amostragem, onde xl, xu são as restrições tipo caixaincluidas na função de restrições de desigualdade g(x).
A metade r do tamanho de Ωs é obtida de:
rdef= ‖x1 − x0‖2, (5.18)
isto é, o primeiro cumprimento de passo (magnitude do progresso) é uma boamedida do raio r da região Ωs para ser explorada na amostragem por projeto deexperimentos.
Um projeto de experimentos (DOE) é realizado somente uma vez para a primeiraconstrução (segunda iteração, k = 1) do modelo de aproximação. Para esta tarefapode se usar qualquer um dos procedimentos de amostragem por DOE contempladosna Seção 3.4.
Como o modelo de aproximação usado aqui é um polinômio linear nãointerpolante LHS produz N (N > n) novos pontos dentro de Ωs para formar amatriz M ∈ R(p+1)×(p+1). Cada linha de M é uma amostra independente adquiridapor LHS. A matriz M é comum para os modelos aproximados da função objetivo edas restrições.
Um importante quesito para tentar manter a convergência do algoritmo emelhorar a qualidade da aproximação é o requisito de consistência. Aqui, somente éimposto o requisito de consistência nas funções. Como definido por (4.6), os modelosde aproximação devem ter o mesmo valor que as funções originais quando avaliadasna iteração atual xk. Para isto é realizada uma mudança de variáveis (ver Eqs. (4.7)e (4.8)). Portanto, os coeficientes do modelo de aproximação são obtidos por WLS
58
através de (4.9) ou (4.10), onde a ponderação dada pela Eq. (4.11) outorga maiorpeso (importância) aos pontos mais próximos ao ponto atual. Tudo isto faz que asaproximações sejam localmente mais precisas.
Uma vez construídos os modelos aproximados, as derivadas aproximadas seobtém facilmente da Eq. (4.12) como explicado na Seção 4.1.2. A atualização dosmodelos aproximados nas outras iterações consiste na inclusão do ponto atual x(k)
no conjunto de amostragem, sem requirer outro procedimento de amostragem, e noposterior ajuste da mesma forma como foi descrito anteriormente (i.e., shifting eWLS).
Os gradientes aproximados obtidos são usados no cálculo da direção de busca,onde é efetuada uma busca linear inexata usando as funções originais. Quandoo passo obtido na busca linear for muito pequeno, isto pode indicar uma pobreaproximação possivelmente pelo mal condicionamento da matriz de ajuste M devidoa uma distribuição pobre dos pontos não compensada pela ponderação. Quandoacontece isto se realiza uma vez diferenças finitas para obter um melhor gradiente,avançar e obter um novo ponto x(k+1). Este último procedimento é apenas usadoquando necessário auxiliando a convergência de um algoritmo que sempre obtémum decremento da função objetivo, por isto foi denominado aqui de passo auxiliar(spacer step, ver LUENBERGER [41]). O procedimento continua até que um critériode convergência seja atingido.
A primeira versão do novo algoritmo proposto, i.e., Algoritmo de DireçõesViáveis e Ponto Interior com Mínimos Quadrados Ponderados e GradientesAproximados (FDIPA-WLS-AG), é enunciado a seguir.
Algoritmo: FDIPA com Gradiente Aproximado usando MínimosQuadrados Ponderados (FDIPA-WLS-AG)
Parâmetros: α ∈ (0, 1), ϕ > 0, η ∈ (0, 1) e ν ∈ (0, 1).Dados: Valores iniciais para x(0) ∈ Ω, λ(0) ∈ Rm, positivo, B(0) ∈ Rn×n simétrica edefinida positiva.
Passo 0. Teste do critério de convergência.Passo 1. Obter os modelos de aproximação e suas derivadas correspondentes.
(i) Se for a primeira iteração principal (k=0), então faça diferenças finitas paraobter os gradientes aproximados: ∇f , ∇gi, e retorne ao Passo 2.
(ii) Se estiver na segunda iteração principal (k=1), então:Calcule r da Eq. (5.18) e obtenha M(k) através de amostragem por LHS emΩs.
59
(iii) Mudança de variável (Shifting): após a mudança de variáveis obtenha a matrizde ajuste χ da Eq. (4.7) e o vetor de resposta f para a função objetivo e asrestrições a partir da Eq. (4.8).
(iv) WLS : obtenha o vetor de ponderação w aplicando a Eq. (4.11) nos p pontosde amostragem e resolva o problema de WLS para os parâmetros α aplicandoa Eq. (4.10) às respostas da função objetivo e das restrições.
(v) Obtenha as derivadas aproximadas: ∇f e ∇g. Isto é da Eq. (4.12) nosparâmetros α do modelo correspondente.
Passo 2. Cálculo da direção de decida viável aproximada d.
(i) Ache (d0, λ0) resolvendo o sistema linear: S ∇gΛ∇gT G
d0
λ0
= −∇f
0
, (5.19)
Se ‖d0‖2 = 0, então Pare.
(ii) Encontre (d1, λ1) resolvendo o sistema linear: S ∇gΛ∇gT G
d1
λ1
= 0−λ
. (5.20)
(iii) Seja a função potencialφ(x) := f(x). (5.21)
(iv) Se dT1∇φ(x(k)) > 0, faça
ρ = minϕ‖d0‖2
2; (α− 1) dT0∇φ(x(k))dT1∇φ(x(k))
(5.22)
Caso contrário, façaρ = ϕ‖d0‖2
2. (5.23)
(v) Calcule a direção viável: d = d0 + ρd1.
Também: ˆλ = λ0 + ρλ1.
Passo 3. Critério de Spacer step.Seja δ > 0 um valor muito pequeno.Avalie f(x(k) + δd),g(x(k) + δd).Se f(x(k) + δd) > f(x(k)) ou gi(x(k) + δd) > 0 para algum i, então:Spacer step: realize diferenças finitas para obter os gradientes aproximados: ∇f ,
60
∇gi, e retorne ao Passo 2.Passo 4. Busca linear aproximada.Encontre um cumprimento de passo t(k) que satisfaça um dado critério de busca
linear restrita, na função potencial φ(x(k) + t(k)d) e tal que
gi(x(k) + t(k)d) < 0 se λi ≥ 0, (5.24)
ougi(x(k) + t(k)d) ≤ gi(x(k)) (5.25)
caso contrario.Passo 5. Atualizações.
(i) Atualize a matriz do modelo:
M(k+1) :=M(k)
x(k)T
.(ii) Determine o novo ponto:
x(k+1) := x(k) + t(k)d.
(iii) Atualize o conjunto de pontos de amostragem: Faça Y (k+1) = Y (k) ∪x(k+1)
(iv) Defina novos valores para λ(k+1)
i > 0 e B(k+1) simétrica e definida positiva.
(v) Retorne ao Passo 0.
5.4 Algoritmo FDIPA-AGO algoritmo descrito na seção anterior pode ter problemas de perda de precisão
ou eficiência quando os pontos obtidos no processo de otimização vão ficandodistribuídos de uma maneira pouco uniforme levando o sistema linear que representaos modelos aproximados a números de condição muito altos, o que se manifestaem erros excessivos na aproximação do gradiente. Quando isto acontece o númerode chamadas de passos auxiliares, i.e., a quantidade de cálculos necessários pordiferenças finitas faz com que a economia em custo computacional seja perdida.
Para superar os problemas da primeira versão foi criada uma segunda versãodo novo algoritmo. Esta versão também é baseada no algoritmo FDIPA,proporcionando soluções viáveis em cada iteração onde a informação do gradienteé obtida através de modelos aproximados da função objetivo e das restrições.Particularmente, os modelos escolhidos para esta versão são polinômios lineares
61
interpolantes. O erro na aproximação foi caracterizado teoricamente na Seção 3.1.2,podendo assim usar tal caracterização como descrito na Seção 4.2 para fazer asaproximações consistentes e convergentes aos valores reais, isto com a intenção deque o novo algoritmo seja robusto e procure manter a convergência global.
De novo, a busca linear é feita numa direção de descida viável aproximada que éobtida resolvendo os sistemas lineares do FDIPA usando os gradientes aproximadosobtidos pelo procedimento descrito na Seção 4.1.1.
Os modelos aproximados que fornecem a informação do gradiente são obtidos porum processo de amostragem e interpolação linear, onde a amostragem consideradaaqui é o procedimento simplex descrito na Seção 3.4. Esta amostragem fornece umanova distribuição uniforme de pontos (ou seja, um conjunto bem posto), e ela é feitana primeira iteração e também quando a direção de busca não seja de descida ou nãoseja viável. Nos outros casos (quando a direção é de descida e viável) se mantém oconjunto de pontos atual apenas substituindo o pior dos pontos y(k),s (i.e., aquele quepiora o condicionamento do sistema linear de interpolação) pelo novo ponto obtidox(k+1). A seleção do ponto a ser substituído para diminuir o erro na aproximaçãodos gradientes é feita tendo em conta o fator geométrico e o tamanho do conjuntode interpolação, procurando que ambas quantidades sejam o mínimo possível. Parafazer esta escolha se multiplica a estimativa do fator geométrico
∥∥∥(M0 (k)\i )−1
∥∥∥ como tamanho (∆(Y )) do conjunto de pontos amostrados. O ponto que tiver o menorvalor relativo nesta conta é o ponto que ao ser substituído diminuirá o erro naaproximação dos próximos gradientes.
A segunda versão do novo algoritmo apresentado, algoritmo de direçõesviáveis com derivadas aproximadas (FDIPA-AG, Feasible Directions Interior PointAlgorithm with Approximated Gradient), é descrito na Seção 5.4.1.
5.4.1 Descrição do algoritmo
Algoritmo: Algoritmo de Ponto Interior e Direções Viáveis comDerivadas Aproximadas (FDIPA-AG)
Parâmetros: α ∈ (0, 1), ϕ > 0, η ∈ (0, 1) e ν ∈ (0, 1).∆(0) > 0, τ ∈ (0, 1), 0 < ε∗f << 1, 0 < ε∗L < ε
(0)L << 1,
0 < δ << 1.Dados: x(0) ∈ Ω, λ(0) ∈ Rmpositivo, B(0) ∈ Rn×n simétrica e definida positiva.
Passo 0: Inicialização.Selecione um conjunto de interpolação bem posto Y (0) = y(0),0,y(0),1, . . . ,y(0),ncom y(0),0 = x(0). Partindo do conjunto de interpolação obtenha os modelos lineares
62
de interpolação f (0) e g(0)i , i = 1, . . . ,m, em torno de x(0). Faça λ = λ(0), k = 0 e
c = 0.
Passo 1: Teste de ponto crítico.
Passo 1a: Seja f (c) = f (k),
g(c)i = g
(k)i .
Seja ∇xL(c)(x,λ) = ∇xf
(c)(x) +∑mi=1 λi∇xg
(c)i (x) o gradiente do Lagrangiano
aproximado correspondente.
Passo 1b: Se k > 0 e ‖∇xL(c)(x(k))‖ < ε∗L ou
∣∣∣f(x(k))− f(x(k−1))∣∣∣ < ε∗f Pare.
Se ‖∇xL(c)(x(k))‖ < ε
(c)L , faça ε(c+1)
L = τ‖∇xL(c)(x(k))‖ construa os modelos
f (c+1) e g(c+1)i bem postos em B(x(k), ε
(c+1)L ), faça c = c + 1 e comece o Passo
1b de novo.
Passo 1c: Faça f (k) = f (c),
g(k)i = g
(c)i .
Passo 2: Cálculo da direção de busca d.
Passo 2a: Calcule (d0, λ0) e (d1, λ1) resolvendo o sistema: B(k) ∇g(x(k))Λ(k)∇g(x(k))T G(x(k))
d0 d1
λ0 λ1
= −∇f(x(k)) 0
0 λ(k)
. (5.26)
Passo 2b: Se dT1∇f(x(k)) > 0, faça
ρ = minϕ‖d0‖2
2; (α− 1) dT0∇f(xk)
dT1∇f(xk)
.Caso contrario, faça
ρ = ϕ‖d0‖22.
Passo 2c: Calcule a direção de busca d e o vetor ˆλ:
d = d0 + ρd1,
ˆλ = λ0 + ρλ1.
Passo 3: Teste da direção d.
Passo 3a: Avalie f(x(k) + δd),g(x(k) + δd).
63
Passo 3b: Se f(x(k) + δd) > f(x(k)) ou gi(x(k) + δd) > 0 para algum i, então:Seja ∆(k) = ∆(Y (k)) = max
1≤i≤n‖y(k),i − y(k),0‖.
Seja a matriz M0 (k) = M0 (k)
∆(k) , onde
M0 (k) =
y
(k),11 − y(k),0
1 · · · y(k),1n − y(k),0
n... ... ...
y(k),n1 − y(k),0
1 · · · y(k),nn − y(k),0
n
.
Seja K(k) = K(Y (k)) =∥∥∥(M0 (k))−1
∥∥∥estim.
.Faça ∆(k)
(new) < ∆(k) e obtenha um conjunto de interpolação Y (k)(new)
bem posto em B(x(k),∆(k)(new)) com K
(k)(new) ≤ K(k).
Obtenha os modelos de interpolação f (k)(new) e g
(k)(new),i.
Faça Y (k) = Y(k)
(new), f (k) = f(k)(new), g
(k)i = g
(k)(new),i.
Retorne ao Passo 2.
Passo 4: Busca linear. Achar um cumprimento de passo t(k) como o primeironúmero da seqüência 1, ν, ν2, ν3, . . . que satisfaz:
f(x(k) + t(k)d) ≤ f(x(k)) + t(k)η∇f(x(k))d,
gi(x(k) + t(k)d) < 0 se ˆλi ≥ 0, ou
gi(x(k) + t(k)d) ≤ gi(x(k)) se ˆλi < 0, i = 1, . . . ,m.
Passo 5: Substituição. Faça x(k+1) = x(k) + td.Seja M0 (k)
\i a matriz M0 (k) após de tirar y(k),i e incluir x(k+1), onde agoray(k),0 = x(k+1). Faça Y (k+1) = Y (k) \
y(k),s
∪x(k+1)
para
y(k),s = arg maxy(k),i∈Y (k)
∥∥∥(M0 (k)\i )−1
∥∥∥estim.
‖y(k),i − x(k+1)‖.
Passo 6: Atualização. Se Y (k+1) 6= Y (k), obtenha os modelos de interpolaçãof (k+1) e g(k+1)
i , i = 1, . . . ,m, em torno de x(k+1) usando Y (k+1). Defina novos valorespara λ(k+1)
i > 0 e B(k+1) simétrica e definida positiva. Faça k = k + 1 e retorne aoPasso 1.
O algoritmo parte de um ponto viável x(0), gera um conjunto de interpolaçãobem posto Y (0) e obtém os modelos lineares de interpolação iniciais para a funçãoobjetivo f (0) e para as restrições g(0)
i , i = 1, . . . ,m.O passo 1 testa o ponto atual x(k) para verificar se verdadeiramente é um ponto
crítico, em cujo caso para e retorna o ponto atual como solução.
64
Se o ponto atual não for ponto crítico continua com o passo 2 onde se procurauma direção de descida viável aproximada d. Como d pode não ser viável ou dedescida devido a seu carácter aproximado é necessário fazer um teste para verificaristo, que é o que se realiza no passo 3. Se a direção for de descida e viável continuacom o passo 4, caso contrário, se geram novos conjuntos amostrais tendo em contaa geometria para melhorar a qualidade da aproximação.
Se a direção de busca d passa pelo teste do passo 3 se procede a realizar umabusca linear restrita com o critério de Armijo no passo 4.
No passo 5 se atualiza o ponto atual x(k) e o conjunto de amostragem,substituindo o ponto que mais empobrece a qualidade da aproximação pelo novoponto x(k+1).
No passo 6 se atualizam os modelos, os multiplicadores de Lagrange e a Hessianaaproximada (por BFGS). Após disto se retorna ao passo 1 para testar se foi obtidoum ponto crítico e se continua assim até satisfazer o criterio de parada.
65
Capítulo 6
Resultados numéricos
Na prática da engenharia, os métodos de otimização com derivadas aproximadas,e também os que não usam derivadas se aplicam geralmente na resolução deproblemas de otimização com poucas variáveis de projeto (n < 50). Isto se deveao número requerido de avaliações das funções, ou seja, da função objetivo e dasrestrições, que cresce drasticamente com o número de variáveis de projeto, comoassinalado na Seção 2.1. No entanto com poucas variáveis de projeto os métodoscom derivadas aproximadas são mais robustos e eficientes, permitindo um maiorcampo de ação e uma maior flexibilidade ao usuário.
A seguir tratar-se-á problemas de otimização que podem ser representadoscom poucas variáveis de projeto. Especificamente apresenta-se um problemateste algébrico e vários exemplos de otimização estrutural aplicados a estruturasreticuladas. No Capítulo 8 se apresentarão resultados de vários problemas deotimização com materiais compósitos laminados.
Os exemplos foram realizados usando o código do algoritmo FDIPA-AGdesenvolvido neste trabalho. Para efeito de comparação, cada exemplo tambémfoi resolvido usando o algoritmo de ponto interior de direções viáveis FDIPA [1],calculando as derivadas por diferenças finitas (FDIPA-FD); e o algoritmo genéticode lagrangiano aumentado (ALGA, Augmented Lagrangian Genetic Algorithm)implementado no Optimization Toolbox de MATLAB.Todos os algoritmos foram implementados em linguagem MATLAB.
6.1 Exemplo 1 - Problema Algébrico de Barnes
O problema de BARNES [43] é um problema com duas variáveis de projetocontinuas (x1, x2) e três restrições não lineares de desigualdade. Este problemaé relativamente fácil de resolver e comumente utilizado para testar algoritmos de
66
otimização baseados em métodos aproximados [44]. A seguir se apresentam a funçãoobjetivo e as restrições para o problema de Barnes [45].
Função objetivo:
f(x) = 75.196− 3.8112x1 − 0.0020567x31 + 1.0345× 10−5x4
1
− 6.8306x2 + 0.030234x1x2 − 1.28134× 10−3x2x21
− 2.266× 10−7x41x2 + 0.25645x2
2 − 0.0034604x32 + 1.3514× 10−5x4
2
− 28.106/(x2 + 1)− 5.2375× 10−6x21x
22 − 6.3× 10−8x3
1x22
+ 7× 10−10x31x
32 + 3.405× 10−4x1x
22 − 1.6638× 10−6x1x
32
− 2.8673e0.0005x1x2 + 3.5256× 10−5x31x2 + 0.12694x2
1
Restrições:
−x1x2 + 700 ≤ 0−x2 + x2
1/125 ≤ 0−(x2 − 50)2 + 5(x1 − 55) ≤ 0
Restrições de caixa0 ≤ x1 ≤ 75
0 ≤ x2 ≤ 65
Começando desde o ponto viável (x1, x2)=(30,40) os algoritmos considerados(FDIPA, FDIPA-FD, FDIPA-AG e ALGA) convergem ao mínimo local(x∗1, x∗2)=(50,20) após 27/29, 30/101, 29/77 e 7/3660 iterações/avaliações dafunção objetivo respectivamente, onde o algoritmo FDIPA usou derivadas exataspara a função objetivo e as restrições, e FDIPA-FD é sua versão com derivadascalculadas por diferenças finitas. Isto ilustra que o algoritmo FDIPA foi maiseficiente, mas requer das derivadas exatas. O algoritmo FDIPA-AG mostroumelhores resultados do que os obtidos pela versão do FDIPA com diferenças finitas(FDIPA-FD). Já o algoritmo ALGA (com 10 individuos por geração) teve o piordos desempenhos avaliado em termos do número de avaliações da função objetivo.
A Figura 6.1 ilustra o comportamento dos quatro algoritmos. A Figura 6.2ilustra as restrições em cor azul, e as curvas de nível da função de BARNES. AFigura 6.3 mostra a convergência dos algoritmos ao mínimo local com decrementoda função objetivo em cada iteração. Podemos dizer que o comportamento da famíliade algoritmos de direções viáveis foi similar enquanto o algoritmo genético teve umconsiderável custo computacional adicional para atingir a mesma solução. Cabe
67
notar que o ALGA mostrou uma alta confiabilidade como pode se observar nosresultados estatísticos na Tabela 6.1.
Figura 6.1: Desempenho dos quatro algoritmos: FDIPA (azul), FDIPA-FD (preto),FDIPA-AG (vermelho) e ALGA (verde).
x1
x 2
x0
0 10 20 30 40 50 60 70 800
10
20
30
40
50
60
70
80
−120
−100
−80
−60
−40
−20
0
20
40
Figura 6.2: Curvas de nível da função de Barnes e trajeto dos quatro algoritmos.
68
0 1000 2000 3000−35
−30
−25
−20
−15
−10
−5
0
NF: avaliações da função objetivo
f(x)
: val
or d
a fu
nção
obj
etiv
o
FDIPAFDIPA−FDFDIPA−AGALGA
Detalhe para NF no intervalo [0 105]
Figura 6.3: Convergência dos quatro algoritmos no problema de Barnes.
Tabela 6.1: Resultados estatísticos problema de Barnes, para 10 execuçõesindependentes com o algoritmo ALGA.
Função objetivo no ótimo Avaliações da função objetivo
Média Desvio padrão Média Desvio padrão
-31.5583 0.0674 4763 431.948Nota: foram usados os parâmetros padrão do ALGA com tamanho da população 10.
Para uma avaliação mais completa do desempenho dos algoritmos (FDIPA,FDIPA-FD, FDIPA-AG e ALGA) se procedeu a rodar o mesmo problema partindode outros pontos iniciais. Os algoritmos chegaram ao mesmo mínimo local(x∗1, x∗2)=(50,20). Os resultados estão resumidos na Tabela 6.2, onde se apresentamo número de avaliações da função objetivo requeridas por cada algoritmo atéchegar à solução. Observando os resultados podemos afirmar que as conclusõesfeitas anteriormente (resultados com x0 = [30 40]T) são mantidas e podem sergeneralizadas para este problema.
69
Tabela 6.2: Resultados do número da avaliações da função objetivo necessárias parasolucionar o problema de Barnes, partindo de diferentes pontos iniciais.
Ponto Inicial FDIPA FDIPA-FD FDIPA-AG ALGA
x0 = [30 40]T 29 101 77 3660x0 = [40 45]T 26 69 55 4690x0 = [55 40]T 50 127 118 3130
6.2 Exemplo 2 - Treliça de 10 barras
Este é um exemplo de otimização de uma estrutura reticulada formada por 10barras onde se procura minimizar o peso da estrutura com restrições estáticas (natensão e no deslocamento) e dinâmicas (na frequência natural). A Figura 6.4 mostraa geometria da estrutura e as condições de carga, onde as variáveis de projeto são asáreas da seção transversal das barras. A estrutura é feita de alumínio com módulo deYoung E = 1× 107 psi (6.98× 1010 Pa) e densidade ρ = 0.1 lbm/in3 (2770 Kg/m3).Foram impostas restrições de tensão de 25 ksi (172368.925 kPa) para cada membroe restrições de deslocamento de 2.0 in (5.08×10−2 m) na direção vertical em cada nólivre. O limite inferior para todas as variáveis de projeto é de 0.1in2 (0.645×10−4m2).As condições de carga são: 50 kips (222.4111 kN) para cima aplicada nos dois nóssuperiores livres, e 150 kips (667.2333kN) para baixo aplicada nos dois nós inferioreslivres. Adicionalmente impõe-se um limite inferior de 22 Hz na frequência natural.
L
L L
1 2
3 4
5 6
7
8
9
10
Figura 6.4: Estrutura reticulada de 10 barras: L=9.144 m (360 in).
70
As soluções ótimas e o peso de cada estrutura obtida por vários métodos estãoresumidas na Tabela 6.3. O método ALGA não conseguiu resolver o problema(convergência prematura por não atingir ponto viável) enquanto os métodos FDIPA-FD e FDIPA-AG obtiveram a mesma solução (ilustrada na Figura 6.5) após de 285e 358 análises respectivamente como se ilustra na Figura 6.6.
Tabela 6.3: Estrutura ótima obtida para a treliça de 10 barras (seções transversaisem in2).
Elemento no. FDIPA-FD FDIPA-AG
1 24.876 24.8742 0.100 0.1003 25.986 25.9874 13.125 13.1255 0.100 0.1006 1.973 1.9737 13.198 13.1988 15.342 15.3429 17.507 17.50710 0.100 0.100
Peso (lbf) 4731.13 4731.13NF 285 358
Figura 6.5: Estrutura ótima obtida pelos algoritmos FDIPA-FD e FDIPA-AG paraa treliça de 10 barras.
71
0 100 200 3004000
5000
6000
7000
8000
9000
10000
11000
NF: avaliações da função objetivo
f(x)
: val
or d
a fu
nção
obj
etiv
o
FDIPA−FDFDIPA−AG
Figura 6.6: Convergência dos algoritmos FDIPA-FD e FDIPA-AG para a treliça de10 barras.
6.3 Exemplo 3 - Domo de 52 barras
Este problema trata de uma estrutura espacial hemisférica (domo) com 52 barras,como ilustrado na Figura 6.7. Uma massa não estrutural de 50 kg é atada a cadanó livre. A área da seção de cada barra é inicialmente 0.0002 m2, e se permite quevarie entre 0.0001 m2 e 0.001 m2. Como mostrado na tabela 6.4, as 52 barras estãoagrupadas em 8 classes. Todos os membros do domo são modelados como elementostipo barra. Adicionalmente, as três coordenadas de cada nó livre são consideradascomo variáveis independentes. A extensão do movimento de cada coordenada é±2m. A simetria estrutural deve ser mantida durante o processo de otimização.Então, há em total 14 variáveis de projeto independentes. A primeira frequênciaangular não pode ser maior que ω1 = 100rad/s, e a segunda não pode ser menorque ω2 = 180rad/s. O módulo de Young é E = 2.1 × 1011 Pa e a densidade éρ = 7800kg/m3.
72
Figura 6.7: Projeto inicial do domo.
Tabela 6.4: Classificação das barras para o problema do domo.
Grupo no. Barras
1 A-B A-C A-D A-E2 B-F C-H D-J E-L3 B-G C-G C-I D-I D-K EK E-M B-M4 B-C C-D D-E B-E5 F-G G-H H-I I-J J-K K-L L-M F-M6 F-N G-O H-P I-Q J-R K-S L-T M-U7 F-O H-O H-Q J-Q J-S L-S L-U F-U8 G-N G-P I-P I-R K-R K-T M-T M-N
O exemplo constitui um problema de otimização dinâmica altamente não linear,com restrições na frequência numa faixa proibitiva. As Tabelas 6.5 e 6.6 mostram acomparação dos resultados ótimos dos algoritmos FDIPA-FD, FDIPA-AG e ALGA,respectivamente. Pode se ver claramente na Figura 6.8 que o algoritmo com melhordesempenho (menor custo computacional) foi o FDIPA-AG. Os resultados do ALGAsão ligeiramente melhores do que o resultado obtido pelos outros dois algoritmos.No entanto, as restrições na frequência não foram satisfeitas completamente comovisto na Tabela 6.6. Cabe notar que o ALGA mostrou uma baixa confiabilidadecomo pode se observar nos resultados estatísticos na Tabela 6.7.
73
Tabela 6.5: Resultado ótimo para o domo, coordenadas (m) e seção transversal (m2).
Variáveis de projeto FDIPA-FD FDIPA-AG ALGA
ZA 4.2056 4.2025 4.0000XB 2.9391 2.9260 0.7546ZB 3.7000 3.7000 4.1127XF 4.2270 4.2229 4.1909ZF 2.5000 2.5000 2.5370XG 2.8667 2.8628 0.8284A1 1.0000× 10−4 1.0000× 10−4 1.0002× 10−4
A2 1.1057× 10−4 1.1096× 10−4 1.1079× 10−4
A3 1.3153× 10−4 1.3138× 10−4 1.4291× 10−4
A4 1.1343× 10−4 1.1266× 10−4 1.0002× 10−4
A5 1.0000× 10−4 1.0000× 10−4 1.0038× 10−4
A6 1.0000× 10−4 1.0000× 10−4 1.5233× 10−4
A7 1.1442× 10−4 1.1456× 10−4 1.0000× 10−4
A8 1.7998× 10−4 1.8056× 10−4 1.4545× 10−4
Peso (Kgf) 183.293 183.285 178.061NF 1519 673 3316Nota: XB , XF , ZB , ZF são as coordenadas Z e X dos nós B, F;
ZA é a coordenada Z do nó A; XG é a coordenada X do nó G;
Aj é a seção transversal do j-ésimo grupo de barras, j = 1, . . . , 8.
Tabela 6.6: Resultado ótimo para o domo, frequência angular (rad/s).
Frequência angular no. FDIPA-FD FDIPA-AG ALGA
1 98.9854 98.9711 102.32672 190.9019 190.8949 186.00743 190.9019 190.8949 188.53434 192.2024 192.2203 188.53435 195.6094 195.6144 209.0310
74
Tabela 6.7: Resultados estatísticos domo, para 10 execuções independentes com oalgoritmo ALGA.
Função objetivo no ótimo Avaliações da função objetivo
Média Desvio padrão Média Desvio padrão
236.503 34.7428 3780.5 1034.82Nota: foram usados os parâmetros padrão do ALGA com tamanho da população 20.
0 500 1000 1500 2000 2500 3000150
200
250
300
350
400
450
500
550
NF: avaliações da função objetivo
f(x)
: val
or d
a fu
nção
obj
etiv
o
FDIPA−FDFDIPA−AGALGA
Figura 6.8: Convergência dos algoritmos FDIPA-FD, FDIPA-AG e ALGA noproblema do domo.
75
Capítulo 7
Tolerância ao ruído externo
Uma das características mais importantes que se procura num algoritmo deotimização sem derivadas é a capacidade de tolerar ruído externo moderado.Conceder esta característica a um algoritmo baseado no gradiente é um dos objetivosprincipais desta tese.
Considera-se uma função (objetivo e/ou restrição) f que é o resultado da adesãode uma função de ruido φ a uma função suave fs
f(x) = fs(x) + φ(x). (7.1)
Pequenas oscilações em φ podem levar a que f tenha vários mínimos locais quepoderiam emboscar um algoritmo convencional baseado no gradiente, mesmo usandodiferenças finitas. De fato, o esquema de diferenças finitas envolve incrementos tãopequenos nas variáveis de projeto que acabam sempre sendo de ordem igual oumenor à ordem do ruido, o que impossibilita seu uso neste tipo de problema. Aperturbação φ pode, em geral, ser aleatória, determinística ou baseada no resultadode um experimento [46]. Sendo assim, φ poderia retornar ou não o mesmo valorquando chamada duas vezes com o mesmo argumento. Portanto, φ não precisamesmo ser uma função. Para simplificar a apresentação dos resultados neste trabalhoassumimos que φ é definida em todo lugar e limitada.
Dependendo do tipo de ruído podemos propor estratégias diferentes. Porexemplo, se o ruído for sempre aleatório podemos usar um esquema de amostrageme regressão baseado nos resultados da Seção 3.2.2 através do Teorema 3.2.3. Noentanto, no caso geral, o ruído pode não ser aleatório. Assim, no caso geral devemoscontinuar com o processo original de aproximação das derivadas por modelos deinterpolação linear, mas tendo controle no tamanho mínimo admissível da região deamostragem, de maneira que obtenha a maior precisão possível para o nível de ruídopresente. Também no caso geral de ruído moderado podemos fazer uma deflexãona direção de busca com o fim de que esta direção seja mais conservativa (com
76
maior probabilidade de que seja viável em relação às funções sem ruído) elevandoa eficiência dos algoritmos baseados em direções viáveis. Nas Seções 7.1 e 7.2se apresenta uma breve descrição das duas últimas estrategias mencionadas parapoder resolver problemas contendo funções com ruído moderado. Na Seção 7.3 seapresenta o algoritmo FDIPA-AG-N (Feasible Directions Interior Point Algorithmwith Approximated Gradient for Noise Functions) que consiste no algoritmo FDIPA-AG modificado para poder resolver problemas de otimização restrita com ruído nasfunções. Finalmente, na Seção 7.4 são apresentados alguns resultados numéricos napresença de ruído moderado.
7.1 Tamanho mínimo admissível para a região deinterpolação
Nesta seção são apresentados resultados necessários para descrever e analisaro comportamento e a precisão do algoritmo de otimização proposto, na presençade ruído moderado. A idéia fundamental é que o algoritmo FDIPA-AG, por estarbaseado em aproximações do gradiente, já possui uma tolerância ao ruído moderado.No entanto, para problemas com funções da forma (7.1) deve se impor restrições notamanho mínimo permitido à região de interpolação, isto para evitar o inútil esforçode tentar derivar o ruído.
As idéias apresentadas nesta seção são baseadas em resultados da Seção 3.2.2e em algumas características do gradiente aproximado obtido com amostragem porsimplex [47].
Da Seção 4.1.1, partindo de um simples não singular Y com vértices yjnj=0 osistema 4.3 fornece o gradiente por simplex dado por
∇m = (M0)−1δ(f, Y ).
Procura-se uma estimativa do erro no gradiente aproximado, similar à doTeorema 3.1.1, mas que leva em conta o ruído.
Precisaremos de medir o nível de ruído em cada simplex. Para este fim definimospara qualquer conjunto T
‖φ‖T = supy∈T‖φ(y)‖ .
Para uma função f que satisfaz (7.1) temos o seguinte resultado
Teorema 7.1.1. Seja Y um conjunto bem posto para interpolação linear em Rn,contido na bola B(y0; ∆(Y )) com centro em y0 e radio ∆(Y ). Seja f satisfazendo(7.1). Seja fs continuamente diferenciável num domínio aberto Ω contendo B(y0; ∆)e seja ∇fs Lipschitz continuo em Ω com constante vs > 0. O gradiente do modelo
77
de interpolação linear satisfaz, para todos os pontos y em B(y0; ∆), um limite noerro da forma
‖∇fs(y)−∇m(y)‖ ≤ κegs∆ + κegφ‖φ‖Y
∆ , (7.2)
onde κegs = vs(1 + n12
∥∥∥(M0)−1∥∥∥ /2), M0 = M0/∆; e
κegφ = 2n 12κ(M). Sendo κ(M) = ‖M‖ ‖M−1‖ .
Demonstração. O Teorema 3.1.1 aplicado a fs implica
‖∇fs(y)−∇ms(y)‖ ≤ κegs∆, (7.3)
onde κegs = vs(1 + n12
∥∥∥(M0)−1∥∥∥ /2), e M0 = M0/∆.
Agora, visto que ‖δ(φ, Y )‖ ≤ 2n 12 ‖φ‖Y , e ∆(Y ) ≤ ‖M‖ ,
‖∇m(y)−∇ms(y)‖ ≤∥∥∥M−1
∥∥∥ ‖δ(f, Y )− δ(fs, Y )‖ =∥∥∥M−1
∥∥∥ ‖δ(φ, Y )‖
≤ 2n 12∥∥∥M−1
∥∥∥ ‖φ‖Y≤ 2n 1
2κ(M)‖φ‖Y∆ . (7.4)
Finalmente, de (7.3) e (7.4) temos
‖∇fs(y)−∇ms(y)‖ ≤ ‖∇fs(y)−∇ms(y)‖+ ‖− (∇m(y)−∇ms(y))‖
≤ κegs∆ + 2n 12κ(M)‖φ‖Y∆ . (7.5)
Comparando o Teorema 3.1.1 com o Teorema 7.2 temos que este último temuma componente adicional 2n 1
2κ(M)‖φ‖Y∆ que pode chegar a valores consideráveisna medida em que o tamanho ∆ do simplex for muito menor do que o nível do ruído‖φ‖Y . Isto nos leva a modificar a estratégia original de controle do erro no algoritmoFDIPA-AG, que consistia em diminuir ∆ quanto for necessário de maneira irrestritamantendo uma boa geometria. Agora a componente adicional no erro do gradienteaproximado nos impõe uma limitante no tamanho mínimo permitido à região deinterpolação.
Seja a função εgn o erro no gradiente aproximado quando há presença de ruído,definida como
εgn = κegs∆ + κegφ‖φ‖Y
∆ . (7.6)
A maneira de ilustração a Figura 7.1 apresenta εgn para ruído moderado (valorespequenos de ‖φ‖Y ). Tomou-se como referência κegs = κegφ = 1.
78
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 10
1
2
3
4
5
6
7
∆
ε gn
φmax
=1/100
φmax
=1/1000
(a) εgn para ∆ ∈ (0, 1).
0 0.05 0.1 0.15 0.2
0.5
1
1.5
2
∆
ε gn
φmax
=1/100
φmax
=1/1000
(b) εgn para ∆ ∈ (0, 0.2).
Figura 7.1: Erro no gradiente aproximado sob ruído moderado, i.e., ‖φ‖Y ≤ 1%.
Da Figura 7.1 observa-se que para os dois casos de ruído moderado há um valorpara ∆ (∆∗) tal que para valores maiores a primeira componente do erro na Eq. (7.6)é a dominante (εgn ∝ ∆), enquanto para valores menores a segunda componente doerro domina (εgn ∝ 1/∆) chegando rapidamente a valores de erro muito altos (i.e.,εgn →∞ quando ∆→ 0).
Para obter uma estimativa do mínimo valor admissível para ∆, que também é oque gera o menor erro no gradiente, derivamos εgn para ter o ∆ crítico.
∂εgn(∆)∂∆ = κegs − κegφ
‖φ‖Y∆2 .
Então, fazendo ∂εgn(∆)∂∆ = 0 temos
∆∗ =√κegφκegs‖φ‖Y .
Agora da condição de otimalidade de segunda ordem temos que
∂2εgn(∆∗)∂∆2 = 2κegφ
‖φ‖Y(∆∗)3 > 0,
portanto, verifica-se que ∆∗ é um mínimo local.Como κegφ e κegs são da mesma ordem de magnitude podemos adotar (7.7), sem
perda de generalidade, como valor crítico de ∆
∆∗ =√‖φ‖Y . (7.7)
79
7.2 Direção de busca conservativa
Quando há ruído nas restrições podemos levar em conta o nível de ruído para teruma direção de busca mais conservativa, no sentido de ter maior probabilidade deque seja viável. Sendo ε a magnitude do maior erro presente nas restrições, propõe-se uma deflexão da direção de busca para dentro da região viável. Especificamenteesta deflexão é obtida modificando o primeiro sistema de (5.17). Sendo o sistemaoriginal B(k) ∇g(x(k))
Λ(k)∇g(x(k))T G(x(k))
d0
λ0
= −∇f(x(k))
0
,o novo sistema (7.8) deflete d0 para o interior da região viável, numa quantidadeque é proporcional a ε com o fim de compensar o erro no pior dos casos, i.e. parafazer d viável com relação a gsuave(x).
B(k) ∇g(x(k))Λ(k)∇g(x(k))T G(x(k))
d0
λ0
= −∇f(x(k))
ελ(k)
. (7.8)
Lembrando que o segundo sistema é B(k) ∇g(x(k))Λ(k)∇g(x(k))T G(x(k))
d1
λ1
= − 0λ(k)
, (7.9)
e que a direção de busca d é dada por
d = d0 + ρd1,
as Figuras 7.2 e 7.3 ilustram a direção de busca original do FDIPA e a nova direçãoconservativa, respectivamente. A região cinza da Figura 7.3 representa a fronteirada região viável correspondente à restrição g(x) = gsuave(x) + φ(x); que é agoradifusa pela presença do ruído φ, onde φ tem amplitude máxima ‖φ‖Y = ε.
80
Figura 7.2: Direção de busca do FDIPA.
d0 rd1
d
Figura 7.3: Direção de busca conservativa (azul) proposta na presença de restriçõescom ruído.
81
7.3 Algoritmo FDIPA-AG-NA incorporação das estratégias propostas neste Capítulo quando há ruído nas
funções requer pequenas modificações do algoritmo FDIPA-AG, apresentado naSeção 5.4. O algoritmo já modificado foi denominado FDIPA-AG-N (FDIPA withApproximated Gradients for Noise Functions) e é apresentado na Seção 7.3.1.
7.3.1 Descrição do algoritmo
O tamanho mínimo admissível para a região de interpolação ∆∗ =√‖φ‖Y =
√ε,
foi levado em conta nos Passos 1b e 3b do algoritmo FDIPA-AG-N, enquanto aincorporação de uma direção de busca conservativa levou à inclusão da Eq. (7.8) noprimeiro sistema do Passo 2a.
Algoritmo: Algoritmo de Ponto Interior e Direções Viáveis comDerivadas Aproximadas para Funções com Ruido Moderado (FDIPA-AG-N )
Parâmetros: α ∈ (0, 1), ϕ > 0, η ∈ (0, 1) e ν ∈ (0, 1).∆(0) > 0, τ ∈ (0, 1), 0 < ε∗f << 1, 0 < ε∗L < ε
(0)L << 1,
0 < δ << 1.Dados: x(0) ∈ Ω, λ(0) ∈ Rmpositivo, B(0) ∈ Rn×n simétrica e definida positiva.
Passo 0: Inicialização.Selecione um conjunto de interpolação bem posto Y (0) = y(0),0,y(0),1, . . . ,y(0),ncom y(0),0 = x(0). Partindo do conjunto de interpolação obtenha os modelos linearesde interpolação f (0) e g(0)
i , i = 1, . . . ,m, em torno de x(0). Faça λ = λ(0), k = 0 ec = 0.
Passo 1: Teste de ponto crítico.
Passo 1a: Seja f (c) = f (k)
g(c)i = g
(k)i .
Seja ∇xL(c)(x,λ) = ∇xf
(c)(x) +∑mi=1 λi∇xg
(c)i (x) o gradiente do Lagrangiano
aproximado correspondente.
Passo 1b: Se k > 0 e ‖∇xL(c)(x(k))‖ < ε∗L ou
∣∣∣f(x(k))− f(x(k−1))∣∣∣ < ε∗f Pare.
Se ‖∇xL(c)(x(k))‖ < ε
(c)L e ε(c)
L > ∆∗ , faça ε(c+1)L = max∆∗; τ‖∇xL
(c)(x(k))‖construa os modelos f (c+1) e g(c+1)
i bem postos em B(x(k), ε(c+1)L ), faça c = c+1
e comece o Passo 1b de novo.
82
Passo 1c: Faça f (k) = f (c)
g(k)i = g
(c)i .
Passo 2: Cálculo da direção de busca d.
Passo 2a: Calcule (d0, λ0) e (d1, λ1) resolvendo o sistema: B(k) ∇g(x(k))Λ(k)∇g(x(k))T G(x(k))
d0 d1
λ0 λ1
= −∇f(x(k)) 0
ελ(k) λ(k)
. (7.10)
Passo 2b: Se dT1∇f(x(k)) > 0, faça
ρ = minϕ‖d0‖2
2; (α− 1) dT0∇f(xk)
dT1∇f(xk)
Caso contrario, faça
ρ = ϕ‖d0‖22.
Passo 2c: Calcule a direção de busca d e o vetor ˆλ:
d = d0 + ρd1.
ˆλ = λ0 + ρλ1
.
Passo 3: Teste da direção d.
Passo 3a: Avalie f(x(k) + δd),g(x(k) + δd).
Passo 3b: Se f(x(k) + δd) > f(x(k)) ou gi(x(k) + δd) > 0 para algum i, então:Seja ∆(k) = ∆(Y (k)) = max
1≤i≤n‖y(k),i − y(k),0‖.
Se ∆(k) ≤ ∆∗, Pare.Seja a matriz M0 (k) = M0 (k)
∆(k) , onde
M0 (k) =
y
(k),11 − y(k),0
1 · · · y(k),1n − y(k),0
n... ... ...
y(k),n1 − y(k),0
1 · · · y(k),nn − y(k),0
n
.
Seja K(k) = K(Y (k)) =∥∥∥(M0 (k))−1
∥∥∥estim.
.Faça ∆(k)
(new) = max∆∗; 12∆(k) e obtenha um conjunto de interpolação
Y(k)
(new)
bem posto em B(x(k),∆(k)(new)) com K
(k)(new) ≤ K(k).
83
Obtenha os modelos de interpolação f (k)(new) e g
(k)(new),i.
Faça Y (k) = Y(k)
(new), f (k) = f(k)(new), g
(k)i = g
(k)(new),i.
Retorne ao Passo 2.
Passo 4: Busca linear. Achar um comprimento de passo t(k) como o primeironúmero da seqüência 1, ν, ν2, ν3, . . . que satisfaz:
f(x(k) + t(k)d) ≤ f(x(k)) + t(k)η∇f(x(k))d,
gi(x(k) + t(k)d) < 0 se ˆλi ≥ 0, ou
gi(x(k) + t(k)d) ≤ gi(x(k)) se ˆλi < 0, i = 1, . . . ,m.
Passo 5: Substituição. Faça x(k+1) = x(k) + td.Seja M0 (k)
\i a matriz M0 (k) após de tirar y(k),i e incluir x(k+1), onde agoray(k),0 = x(k+1). Faça Y (k+1) = Y (k) \
y(k),s
∪x(k+1)
para
y(k),s = arg maxy(k),i∈Y (k)
∥∥∥(M0 (k)\i )−1
∥∥∥estim.
‖y(k),i − x(k+1)‖
Passo 6: Atualização. Se Y (k+1) 6= Y (k), obtenha os modelos de interpolaçãof (k+1) e g(k+1)
i , i = 1, . . . ,m, em torno de x(k+1) usando Y (k+1). Defina novos valorespara λ(k+1)
i > 0 e B(k+1) simétrica e definida positiva. Faça k = k + 1 e retorne aoPasso 1.
7.4 Resultados numéricos
Nesta seção adaptaremos alguns dos problemas apresentados no Capítulo 6 paramostrar a capacidade do algoritmo FDIPA-AG-N de tolerar ruído moderado. Amodificação consiste em adicionar a cada valor estimado (tensão, deslocamento,frequência natural) ruído aleatório proporcional à magnitude do valor calculado. Porexemplo, o erro a ser adicionado ao funcional fsuave = fsuave(x(k)) vai ser um valorpseudoaleatório uniformemente distribuído no intervalo
[− fac
100 |fsuave| , + fac100 |fsuave|
],
sendo “fac” o valor absoluto do fator (em [%]) de erro adicionado, o que é equivalentea substituir fsuave por
fperturbado = fsuave
(1 +X
fac100
), (7.11)
onde X é um valor pseudoaleatório uniformemente distribuído no intervalo [−1, 1]3.3Em MATLAB a rotina RAND retorna um valor pseudoaleatório de uma distribuição uniforme
no intervalo [0, 1] . Para obter X em [−1, 1] tomou-se: X = −1 + 2 · RAND.
84
7.4.1 Treliça de 10 barras
Tomando o problema da Seção 6.2 sem a restrição na frequência natural,adicionamos ruído aleatório nas restrições de tensão e deslocamento substituindofsuave por fperturbado conforme (7.11), onde fsuave está representando a tensão e odeslocamento.
As soluções (x), peso e número de avaliações (NF) obtidas pelo algoritmo FDIPA-AG e sua versão adaptada para o ruído FDIPA-AG-N estão resumidas na Tabela7.1 para diversos níveis de ruído moderado (fac = 0.1% e fac = 1%). O métodoFDIPA-FD, i.e. FDIPA com diferenças finitas, não conseguiu avançar na presençade ruído externo nem mesmo para o nível mais baixo de ruído estudado, no entantoo método FDIPA-FD teve um desempenho ligeiramente melhor do que o FDIPA-AG na ausência de ruído. Também observa-se que embora o algoritmo FDIPA-AGconsiga uma solução na presença de ruído as modificações específicas para o ruídofizeram que o algoritmo FDIPA-AG-N tivesse um desempenho melhor, diferençaque se faz mais notória na medida em que o nível do ruído for maior. A Figura 7.4mostra a convergência destes algoritmos. Como se esperava, pequenas imprecisõesnos resultados estão diretamente relacionadas ao nível de ruído adicionado, entremaior nível de ruído mais imprecisa fica a solução.
Tabela 7.1: Estrutura obtida para a treliça de 10 barras sob ruído moderado nastensões e deslocamentos (seções transversais em in2).
SEM RUÍDO EXTERNO RUÍDO EXTERNO: 0.1% RUÍDO EXTERNO: 1%
Elem. no. FDIPA-FD FDIPA-AG FDIPA-AG FDIPA-AG-N FDIPA-AG FDIPA-AG-N
1 23.5308 24.3699 23.6153 25.5014 29.2092 27.74522 25.2851 25.0905 26.2749 23.6582 26.5784 24.08553 0.1000 0.1000 0.1032 0.1001 0.1144 0.10774 14.3745 14.1133 15.1628 13.7865 11.3119 12.62075 0.1000 0.1000 0.1015 0.1011 0.1071 0.10506 1.9697 2.0425 1.9695 1.9878 3.2108 3.59977 12.8277 13.3911 14.4397 13.3633 11.8428 12.51078 0.1000 0.1001 0.1019 0.1009 0.1221 0.10699 12.3906 12.7012 12.8072 13.0429 12.7949 13.351410 20.3286 19.2961 17.7463 19.6924 19.9253 19.7653
Peso (lbf) 4676.92 4685.26 4716.05 4697.04 4794.14 4785.9NF 300 487 514 496 635 536
85
0 100 200 300 400 500 6004000
5000
6000
7000
8000
9000
10000
11000
NF: avaliações da função objetivo
f(x)
: val
or d
a fu
nção
obj
etiv
o
FDIPA−FDFDIPA−AGFDIPA−AG, fac=0.1%FDIPA−AG−N, fac=0.1%FDIPA−AG, fac=1%FDIPA−AG−N, fac=1%
(a) Peso vs. número de avaliações. Intervalo completo
100 200 300 400 5004500
5000
5500
6000
6500
NF: avaliações da função objetivo
f(x)
: val
or d
a fu
nção
obj
etiv
o
FDIPA−FDFDIPA−AGFDIPA−AG, fac=0.1%FDIPA−AG−N, fac=0.1%FDIPA−AG, fac=1%FDIPA−AG−N, fac=1%
(b) Peso vs. número de avaliações. Intervalo: NF em (50, 500)
Figura 7.4: Convergência dos algoritmos FDIPA-FD, FDIPA-AG e FDIPA-AG-Npara a treliça de 10 barras no caso de ruído nas restrições.
86
7.4.2 Domo de 52 barras
Tomando o problema da Seção 6.3, adicionamos ruído aleatório nas restriçõesda frequência natural substituindo fsuave por fperturbado conforme (7.11), onde fsuave
está representando a frequência natural da estrutura.As soluções (x), peso e número de avaliações (NF) obtidas pelo algoritmo FDIPA-
AG e sua versão adaptada para o ruído FDIPA-AG-N estão resumidas na Tabela 7.2para diversos níveis de ruído moderado (fac = 0.1% e fac = 1%). Como esperado,o método FDIPA-FD, i.e. FDIPA com diferenças finitas, não conseguiu avançar napresença de ruído externo nem mesmo para o nível mais baixo de ruído estudado. Jáo método FDIPA-AG teve um desempenho notoriamente melhor do que o FDIPA-FD na ausência de ruído, manifestado em que este último precisou de mais dodobro de avaliações da função objetivo para convergir à mesma solução. Tambémobserva-se que embora o algoritmo FDIPA-AG consiga uma solução na presença deruído o algoritmo FDIPA-AG-N teve um desempenho melhor, diferença que se fazmais notória na medida em que o nível do ruído for maior. A Figura 7.5 mostra aconvergência destes algoritmos. De novo se observa a relação direta entre o nível deruído adicionado e as imprecisões nos resultados, entre maior nível de ruído maismais imprecisa fica a solução.
Tabela 7.2: Estrutura obtida para o domo sob ruído moderado nas frequênciasnaturais, coordenadas (m) e seção transversal (cm2).
SEM RUÍDO EXTERNO RUÍDO EXTERNO: 0.1% RUÍDO EXTERNO: 1%
Variável FDIPA-FD FDIPA-AG FDIPA-AG FDIPA-AG-N FDIPA-AG FDIPA-AG-N
ZA 4.2056 4.2025 4.0888 6.7458 6.1918 4.1013XB 2.9391 2.9260 1.9056 2.4588 1.2577 1.2989ZB 3.7000 3.7000 3.8546 3.8518 4.4909 4.2000XF 4.2270 4.2229 3.3463 3.8664 3.0595 2.9446ZF 2.5000 2.5000 2.8298 2.6126 3.2015 3.1153XG 2.8667 2.8628 2.3353 2.7451 2.2068 2.1194A1 1.0000 1.0000 1.1721 1.1985 1.8997 1.2024A2 1.1057 1.1096 1.2410 1.5599 2.1367 1.1004A3 1.3153 1.3138 1.5323 1.2852 2.0080 2.0519A4 1.1343 1.1266 1.2263 1.0080 1.9307 1.1046A5 1.0000 1.0000 1.4551 1.4404 1.9929 1.9562A6 1.0000 1.0000 1.3688 1.0130 1.4720 1.7049A7 1.1442 1.1456 1.6017 1.4656 1.5489 1.8049A8 1.7998 1.8056 1.6403 1.3180 1.5668 1.6945
Peso (Kgf) 183.293 183.285 206.782 190.619 243.780 238.506NF 1519 673 932 856 1556 1369
87
0 500 1000 1500150
200
250
300
350
400
450
500
550
NF: avaliações da função objetivo
f(x)
: val
or d
a fu
nção
obj
etiv
o
FDIPA−FDFDIPA−AGFDIPA−AG, fac=0.1%FDIPA−AG−N, fac=0.1%FDIPA−AG, fac=1%FDIPA−AG−N, fac=1%
Figura 7.5: Convergência dos algoritmos FDIPA-FD, FDIPA-AG e FDIPA-AG-Npara o domo de 52 barras no caso de ruído nas restrições.
88
Capítulo 8
Otimização de estruturas feitascom compósitos laminados
O uso de materiais compósitos em projeto estrutural tem se incrementadonas últimas três décadas devido a várias vantagens que esses materiais oferecemcomparados com os materiais tradicionais, tais como aço, alumínio, e várias ligas.Uma das principais razões para sua popularidade é seu baixo peso. Adicionalmentea sua vantagem em peso, alguns compósitos possuem maior rigidez e propriedades deresistência ao serem comparados com metais. A vantagem em peso pode também serefletir nas propriedades de desempenho descritas acima ao introduzir as definiçõesde rigidez específica e resistência específica. Estas últimas são definidas pela divisãoda propriedade em questão pela densidade do material. Portanto, ainda que a rigidezou resistência de um material compósito seja comparável às duma liga convencional,as vantagens de uma alta rigidez específica e/ou resistência específica fazem oscompósitos mais atraentes do que ligas metálicas.
Embora os materiais compósitos tenham mecanismos de falha desconhecidos eseu processo de fabricação tenha um custo elevado; foram essas vantagens, citadasno parágrafo anterior, que inicialmente estimularam a utilização generalizada decompósitos em projeto estrutural, especialmente para aplicações onde o peso écrítico tais como as encontradas na indústria aeroespacial; ver e.g. ZAGAINOV& LOZINO-LOZINSKI [48]. Além de sua vantagem em peso, o uso de materiaiscompósitos em projeto estrutural teve implicações significantes em termos daformulação do problema de projeto, particularmente em termos de projeto ótimo.Com a introdução de materiais compósitos na formulação de um problema deprojeto, os projetistas ganharam uma nova flexibilidade no uso de variáveis quemudam diretamente as propriedades do material. Com o propósito de melhoraro desempenho estrutural e atingir os requerimentos de uma situação específica deprojeto, agora é possível ajustar as propriedades do material além das dimensões daestrutura.
89
8.1 Materiais compósitos laminados
Os materiais compósitos podem ser definidos de maneira geral como materiaisque são feitos de dois ou mais constituintes, e cujas características são derivadasdos componentes individuais. O material compósito pode ter as característicascombinadas dos constituintes ou pode ter características muito diferentes,dependendo na maneira em que os componentes são combinados. De fato, às vezesas propriedades dos materiais compósitos podem até exceder as dos constituintes[49].
8.1.1 Classificação dos materiais compósitos
Várias abordagens podem ser usadas para classificar materiais compósitos. Umadelas é baseada na natureza dos materiais constituintes, classificando-os em duascategorias: orgânicos ou inorgânicos. A madeira é historicamente o materialestrutural mais comum de natureza orgânica. Ele é um material compósito porser composto de dois constituintes distintos: fibras duras embebidas numa estruturade células mais suaves (material matriz).
O uso de materiais puramente orgânicos tais como madeira não é apropriadoem muitas situações por causa da baixa rigidez e propriedades de resistência,sensibilidade a efeitos ambientais tais como temperatura e umidade, e baixaresistência a substâncias químicas e solventes. Em compósitos inorgânicos osconstituintes podem ser classificados como metálicos ou não metálicos. A principaldesvantagem dos constituintes metálicos é o alto peso específico, que é crítico emaplicações tais como as encontradas na indústria aeroespacial. No entanto, eles sãoum dos melhores constituintes em se tratando de resistência a alta temperatura.
Os materiais compósitos usados com maior frequência são os feitos de todosos constituintes não metálicos, tais como constituintes poliméricos assim comosubstâncias inorgânicas similares ao vidro. Materiais de matriz polimérica podemser classificados como termoestáveis ou termoplásticos, ambos tendo diferentescaracterísticas como materiais da matriz. Termoestáveis pré-impregnados podemser usados para fabricar formas complexas já que eles são maleáveis, e requerembaixa temperatura de processamento embora de precisar tempos longos para curado,que é um processo irreversível. Termoplásticos, por outro lado não requeremcurado e podem ser remodelados e reutilizados, mas requerem altas temperaturas deprocessamento e não são maleáveis o que torna difícil seu uso na fabricação de formascomplexas. Termoplásticos típicos são o polietileno, poliestireno e poliéster (PEEK,polyetheretherketone). Um aspecto importante quando comparamos materiais dematriz termoplásticos e termoestáveis é sua temperatura operacional, a qual paratermoplásticos pode alcançar até 260 °C enquanto os materiais termoestáveis podem
90
somente operar até 120 °C. Os termoplásticos também oferecem em geral altaresistência à umidade.
Outra abordagem possível para uma classificação geral dos materiais compósitosé baseada em sua forma. A Figura 8.1 apresenta vários exemplos de materiaiscompósitos com diferentes formas. Quando constituintes em partículas durasestão aleatoriamente dispersas numa matriz suave, como no caso do concreto,o material compósito entra na categoria de compósito em partículas. Quandoflocos finos de constituintes em partículas são adicionados (aleatoriamente oualinhados) ao material matriz, o material compósito é chamado de compósito emflocos. Compósitos em fibras reforçadas são os materiais compósitos mais usados.Geralmente fibras fortes e rígidas, que provem capacidade de suportar cargas, sãoembebidas numa matriz mais suave que é responsável por manter as fibras juntase redistribuir carga de uma fibra para suas adjacentes, no caso de falha da fibra.Finalmente, os compósitos laminados são obtidos por empilhamento de camadasfinas de material (coladas juntas). As camadas podem ser feitas de uma amplavariedade de materiais e podem até mesmo ser compósitos de fibra reforçada, quepodem ser usadas em estruturas de alto desempenho pois isso permite a adaptaçãodas propriedades globais às necessidades específicas.
Dependendo da disposição das fibras, uma camada individual do materialcompósito laminado pode assumir diferentes formas. As fibras podem ser curtase alinhadas numa direção dada ou podem estar distribuídas aleatoriamente nomaterial matriz. As camadas de compósito de fibra reforçada contínua estão feitastipicamente de fibras cilíndricas de diâmetro pequeno, cujo rádio é da ordem de0.005 mm para fibras de carbono até cerca de 0.05 mm para o caso de fibras de boro.Comercialmente, estes materiais estão disponíveis na forma de fita unidirecional eas fibras estão pre-impregnadas com o material matriz (prepregs). Outra forma decamadas compósitas de fibra reforçada é a de tipo tecido, onde um pacote grandede tipicamente 10000 ou mais fibras são tecidas em duas ou mais direções, comdiferentes características de tecido possíveis.
91
(a) Compósito em partículas.
(b) Compósito em flocos.
(c) Compósito em fibras reforçadas.
(d) Compósito laminado.
Figura 8.1: Materiais compósitos com diferentes formas dos constituintes.
8.2 Identificação e otimização de estruturascompósitas laminadas
A redução de vibração em estruturas têm sido o objetivo de intensivas pesquisasna área de materiais compósitos e estruturas adaptativas. Dentro deste marco, umtópico importante é o número ótimo de dispositivos de controle ativo, seu tamanhoe distribuição na estrutura, que pode ser fundamental para o funcionamentorobusto de sistemas de controle ativo. Trabalhos pioneiros sobre o desenvolvimentode metodologias para determinar a distribuição de sensores e atuadores usandotécnicas de otimização heurísticas têm sido conduzidos por FRANCO CORREIAet al. [50], e MOITA et al. [51], usando algoritmos de recozimento simulado.Previamente, SULEMAN & GONÇALVES [52] consideraram varias funções objetivosimultaneamente, tais como maximizar o deslocamento vertical médio de uma viga;minimizar a masa dos atuadores, e minimizar a voltagem de acionamento, ondeas variáveis de projeto eram as coordenadas dos pares de atuadores e o tamanhodos arranjos de atuadores. Também foram consideradas restrições geométricas elimites (superiores e inferiores) nas variáveis de projeto. ADALI et al. [53] tambémconsiderou o problema de uma viga onde procurou minimizar a máxima flambagemvertical de uma viga laminada, usando um par de atuadores, através de uma
92
técnica de projeto robusto. Uma revisão em otimização de estruturas inteligentes eatuadores pode ser encontrada em FRECKER [54].
Neste trabalho serão abordados dois problemas. O primeiro consiste naidentificação de propriedades elásticas dos materiais constituintes de uma placacompósita laminada, onde se procura minimizar uma função de erro que expressao desvio entre frequências de ressonância numéricas e experimentais. Os dadosexperimentais e a formulação do primeiro problema foram obtidos por ARAÚJOet al. [55] usando um vibrômetro laser para medir a resposta das placas, e oestimulo das amostras foi realizado por um martelo de baixo impacto. O segundoproblema consiste no posicionamento ótimo de arranjos de piezoelétricos colados nasuperfície de placas sanduíche com núcleo viscoelástico e camadas de laminado. Oobjetivo é minimizar a massa (e por tanto a quantidade) dos atuadores restringindoos fatores de amortecimento usando um esquema baseado em otimização topológica.Os arranjos piezoelétricos são montados na superfície e se introduz amortecimentoativo através de uma lei de controle proporcional ao deslocamento e velocidademedidas nos sensores.
Para resolver estes dois problemas usaremos o modelo da placa por elementosfinitos desenvolvido recentemente por ARAÚJO et al. [56–59]. O modelo porelementos finitos é apresentado a seguir.
8.2.1 Modelo numérico da placa
Este modelo é usado para obter a resposta dinâmica da placa sanduíche daFigura 8.2. O modelo é baseado numa abordagem mista por camadas onde o núcleoviscoelástico (v) é modelado de acordo com uma teoria de deformação cortantede alta ordem com expansão cúbica dos componentes de deslocamento no plano nacoordenada da espessura. As camadas de compósito laminado (e1, e2) são modeladasde acordo a uma teoria de deformação cortante de primeira ordem para placaslaminadas, e os arranjos de sensor (s) e atuador (a) piezoelétricos também sãomodelados segundo esta teoria. Todos os materiais são lineares, homogêneos eortotrópicos.
93
Figura 8.2: Camadas da placa sanduíche.
Figura 8.3: Lâmina de material ortotrópico com direções principais associadas aoreferencial material (x1, x2, x3) e sistema de coordenadas global (x, y, z).
Para cada lâmina (Figura 8.3) no laminado, as equações constitutivas sãoexpressadas nas direções principais materiais, e para a situação de tensão normal
94
transversal nula podem ser escritas como:
σ11
σ22
σ23
σ13
σ12
=
Q∗11 Q∗12 0 0 0Q∗12 Q∗22 0 0 00 0 Q44 0 00 0 0 Q55 00 0 0 0 Q66
E
ε11
ε22
γ23
γ13
γ12
−
0 0 e∗31
0 0 e∗32
0 e24 0e15 0 00 0 0
E1
E2
E3
,
D1
D2
D3
=
0 0 0 e15 00 0 e24 0 0e∗31 e∗32 0 0 0
ε11
ε22
γ23
γ13
γ12
+
ε11 0 00 ε22 00 0 ε∗33
ε E1
E2
E3
,
(8.1)
onde σij são os componentes da tensão, εij e γij são os componentes dadeformação, Ei e Di são os componentes do campo elétrico e deslocamento elétricorespectivamente, QE
ij e Q∗Eij são os coeficientes de rigidez e de rigidez reduzidos emcampo elétrico constante, eij e e∗ij são as constantes piezoelétricas e piezoelétricasreduzidas respectivamente, e εεij e ε∗ε33 são as constantes dielétricas e dielétricasreduzidas, medidas a deformação constante. O conjunto completo de equaçõesem 8.1 é usado apenas para as camadas de sensor e atuador, enquanto para ascamadas elásticas restantes se considera apenas a primeira equação sem a partepiezoelétrica. As constantes modificadas (reduzidas) seguem a hipótese de tensãonormal transversal nula:
Q∗Eij = QEij −
QEi3Q
Ej3
QE33
, i, j = 1, 2
e∗3i = e3i − e33QEi3
QE33, i = 1, 2
ε∗S33 = εS33 + e233QE
33.
(8.2)
Os coeficientes elásticos em tensão normal transversal nula da Eq. (8.2) podemser expressados em termos de parâmetros adimensionais do material, como indicadona Eq. (8.3):
Q∗E = E1
8α0
8 α4 − α3 0 0 0α4 − α3 8− 2α2 0 0 0
0 0 α8 − α9 0 00 0 0 α8 + α9 00 0 0 0 1
2(8− α2 − 3α3 − α4)
. (8.3)
95
Os parâmetros adimensionais estão definidos em termos das constantes deengenharia como mostrado na Eq. (8.4), onde α0 = 1− ν2
12E2/E1.
α2 = 4− 4E2/E1
α3 = 1 + (1− 2ν12)E2/E1 − 4α0G12/E1
α4 = 1 + (1 + 6ν12)E2/E1 − 4α0G12/E1
α8 = 4(G13 +G23)α0/E1
α9 = 4(G13 −G23)α0/E1
(8.4)
As constantes de engenharia E1 e E2 na Eq. (8.4) são o módulo de elasticidade nasdireções principais materiais x1 e x2 respectivamente, enquanto G12, G13 e G23 são omódulo cisalhante transversal nos planos (x1, x2), (x1, x3) e (x2, x3) respectivamentee ν12 é o principal coeficiente de Poisson. Estas constantes de engenharia podemser recuperadas a partir dos parâmetros adimensionais usando as seguintes relaçõesinversas:
E2/E1 = (4− α2)/4G12/E1 = (8− α2 − 3α3 − α4)/(16α0)
ν12 = (α4 − α3)/(8− 2α2)G13/E1 = (α8 + α9)/(8α0)G23/E1 = (α8 − α9)/(8α0)
α0 = 1− (α4 − α3)2
16(4− α2)
(8.5)
A restrição termodinâmica nos valores das constantes elásticas para materiaisem engenharia é equivalente à afirmação de que a soma do trabalho realizado portodos os componentes da tensão deve ser positiva para evitar a criação de energia.Esta restrição foi generalizada por LEMPRIERE [60] para materiais ortotrópicos,requirindo que as matrizes constitutivas (rigidez e complacência) sejam definidaspositivas, levando aos seguintes limites nos parâmetros de engenharia [61]:
E1, E2, G23, G13, G12 > 0
|ν12| <√E1
E2
(8.6)
As direções principais de ortotropia frequentemente não coincidem com asdireções do sistema de coordenadas que é usado na solução de problemas comlaminados. Então, se faz necessário aplicar transformação de coordenadas paratensão e deformação entre as direções principais materiais (x1, x2, x3) e o sistema decoordenadas geometricamente natural (x, y, z), ilustrado na Figura 8.3. Aplicando
96
uma transformação rotacional partindo dos eixos principais materiais (x1, x2, x3)para o sistema global da placa, de coordenadas retangulares (x, y, z), são obtidas asequações constitutivas neste último sistema coordenado [56–59]:
σxx
σyy
σyz
σxz
σxy
=
Q11 Q12 0 0 Q16
Q12 Q22 0 0 Q26
0 0 Q44 Q45 00 0 Q45 Q55 0Q16 Q26 0 0 Q66
E
εxx
εyy
γyz
γxz
γxy
−
0 0 e31
0 0 e32
e14 e24 0e15 e25 00 0 e36
Ex
Ey
Ez
,
Dx
Dy
Dz
=
0 0 e14 e15 00 0 e24 e25 0e31 e32 0 0 e36
εxx
εyy
γyz
γxz
γxy
+
εxx εxy 0εxy εyy 00 0 εzz
ε Ex
Ey
Ez
,
(8.7)
onde os coeficientes de rigidez reduzidos Qij estão dados por (8.8), e asconstantes piezoelétricas e dielétricas transformadas estão dados por (8.9) e (8.10)respectivamente.
Q11 = Q∗E11 cos4 θ + 2(Q∗E12 + 2QE66) sin2 θ cos2 θ +Q22 sin4 θ
Q12 = (Q∗E11 +Q∗E22 − 4QE66) sin2 θ cos2 θ +Q∗E12 (sin4 θ + cos4 θ)
Q22 = Q∗E11 sin4 θ + 2(Q∗E12 + 2QE66) sin2 θ cos2 θ +Q∗E22 cos4 θ
Q16 = (Q∗E11 −Q∗E12 − 2QE66) sin θ cos3 θ + (Q∗E12 −Q∗E22 + 2QE
66) sin3 θ cos θQ26 = (Q∗E11 −Q∗E12 − 2QE
66) sin3 θ cos θ + (Q∗E12 −Q∗E22 + 2QE66) sin θ cos3 θ
Q66 = (Q∗E11 +Q∗E22 − 2Q∗E12 − 2QE66) sin2 θ cos2 θ +QE
66(sin4 θ + cos4θ)Q44 = QE
44 cos2 θ +QE55 sin2 θ
Q45 = (QE55 −QE
44) cos θ sin θQ55 = QE
55 cos2 θ +QE44 sin2 θ
(8.8)
e31 = e∗31 cos2 θ + e∗31 sin2 θ
e32 = e∗31 sin2 θ + e∗32 cos2 θ
e36 = (e∗31 − e∗32) sin θ cos θe14 = (e15 − e24) sin θ cos θe24 = e24 cos2 θ + e15 sin2 θ
e15 = e15 cos2 θ + e24 sin2 θ
e25 = (e15 − e24) sin θ cos θ
(8.9)
97
εεxx = εε11 sin2 θ + εε22 cos2 θ
εεxx = (εε22 − εε11) sin θ cos θεεzz = ε∗ε33
(8.10)
As equações de movimento para a placa são obtidas aplicando o principio deHamilton, usando um elemento placa serendipity de oito nós com 17 graus deliberdade mecânicos por nó, e um grau de liberdade de potencial elétrico porcamada de piezoelétrico, assumindo que o potencial varia linearmente na direçãoda espessura. Assumindo vibrações harmônicas, as equações de equilíbrio dinâmicopara as vibrações livres podem ser escritas na seguinte forma:
[K∗(ω)− λ∗mM] um = 0, (8.11)
onde M é a matriz de massa, um é o autovetor, λ∗m = λm(1 + iηm) é o autovalorcomplexo, λm = ω2
m e ηm é o fator de amortecimento modal. A matriz de rigidezcondensada é escrita como:
K∗(ω) = K(ω)−[(Gd + iωGv) Ka
uφ + Ksuφ
](Ks
φφ)−1(Ksuφ)T , (8.12)
onde K(ω) é a matriz complexa de rigidez mecânica dependente da frequência,Kuφ é a matriz de rigidez piezoelétrica, Kφφ é a matriz de rigidez dielétrica (ossubíndices a e s se referem ao atuador e sensor respectivamente); Gd e Gv são osganhos de controle do deslocamento e velocidade respectivamente. Assume-se que onúcleo viscoelástico é homogêneo e isotrópico com um módulo cisalhante complexoe dependente da frequência, de acordo a um modelo de derivada fracionária domaterial viscoelástico. Maiores detalhes sobre o modelo podem ser encontrados emARAÚJO et al. [62].
8.3 Identificação das propriedades elásticas dolaminado
O problema inverso da estimação das propriedades elásticas do material, emestruturas adaptativas de compósito laminado pode ser resolvido de diferentesmaneiras. A abordagem usada neste trabalho consiste na minimização de umafunção do erro, que expressa o desvio entre as frequências naturais em vibraçãolivre experimentais e as obtidas com o modelo numérico. Esta expressão do erroentre o modelo e os dados experimentais está relacionada em forma implícita com
98
os parâmetros elásticos e piezoelétricos do material. A expressão proposta para oerro estimado é do tipo mínimos quadrados ponderados:
Φ =I∑
m=1wm
(1− λm
λm
)2
, (8.13)
onde λm são a parte real dos autovalores experimentais, λm são a parte real dosautovalores previsto pelo modelo numérico, wm são os pesos usados para expressaro nível de confiança em cada autovalor experimental, e I é o número total deautovalores experimentais em consideração.
As variáveis de projeto são os coeficientes elásticos adimensionais da relaçãoconstitutiva (8.1) para cada material elástico. A experiência prática em resolvereste tipo de problema inverso mostra que é recomendado o uso de parâmetrosadimensionais para as variáveis de projeto elásticas. Assim o vetor de variáveisde projeto é dado por:
b =[α1, . . . ,αmat
]T, (8.14)
onde mat é o número de materiais diferentes no laminado. Então, para um materiali em particular temos
αi =[Ei
1/0Ei
1, αi2, α
i3, α
i4, α
i8, α
i9
]T, (8.15)
onde 0Ei1 representa o estimativo inicial para Ei
1.
O problema de otimização pode ser formulado como a minimização restrita daestimativa do erro na Eq. (8.13), como segue:
minimizarb
Φ(b) ≥ 0
sujeito a: gj(b) ≤ 0, j = 1, . . . , ngbli ≤ bi ≤ bui , i =, . . . , ndv,
(8.16)
onde bli e bui são restrições laterais nas variáveis de projeto, ndv é o número devariáveis de projeto e ng é o número de restrições. As restrições gj são necessáriaspara garantir que as matrizes de elasticidade sejam definidas positivas para todos osmateriais presentes na estrutura, como descrito na Seção 8.2.1 através da Eq. (8.6).Essas restrições só têm significado para variáveis de projeto elásticas e podem ser
99
escritas na seguinte forma para cada material elástico i:
gi1 ≡ −αi1 ≤ 0
gi2 ≡αi2
αi2 − 4 ≤ 0
gi3 ≡8− αi2 − 3αi3 − αi4
−16αi0≤ 0
gi4 ≡∣∣∣∣∣αi4 − αi38− 2αi2
∣∣∣∣∣−√
44− αi2
≤ 0
gi5 ≡αi5 − αi6−8αi0
≤ 0
gi6 ≡αi5 + αi6−8αi0
≤ 0
(8.17)
ou, em termos das constantes elásticas de engenharia:
gi1 ≡ −Ei
1
Ei(0)1≤ 0
gi2 ≡ Ei2 ≤ Ei
1
gi3 ≡ −Gi
12Ei
1≤ 0
gi4 ≡∣∣∣νi12
∣∣∣−√√√√Ei
1Ei
2≤ 0
gi5 ≡ −Gi
23Ei
1≤ 0
gi6 ≡ −Gi
13Ei
1≤ 0
(8.18)
Para a obtenção das propriedades piezoelétricas, devido à diferença na ordemde magnitude das derivadas das frequências para os diferentes tipos de variáveisde projeto possíveis (elásticas e piezoelétricas) ARAÚJO et al. [63] propõe umprocedimento de identificação em três etapas:
Etapa I: Estimação das propriedades elásticas do material ortotrópico do compósitolaminado base.
Etapa II: Colagem superficial dos sensores/atuadores piezoelétricos; estimação daspropriedades elásticas do piezoelétrico (condição de circuito fechado).
Etapa III: Estimação das propriedades piezoelétricas (condição de circuito aberto).
Nesta seção será abordada apenas a Etapa I, onde a placa considerada éuma placa compósita laminada reforçada por fibra de carbono (CFRP, compositelaminated carbon fibre reinforced plate), numa resina epóxica, com 16 camadas,
100
sendo 200mm × 300mm suas dimensões no plano, espessura total de 4.13mm,e um esquema de laminado simétrico [90°/+45°/0°/-45°]2s, com massa específica eespessura apresentadas na Tabela 8.1.
Tabela 8.1: Dimensões da placa CFRP.
Espessura media Peso Volume Massa específica[mm] [g] [mm3] [kg/dm3]
4.13 382 24750 1.543
Devido à sequência simétrica de apilhamento, espera-se um comportamentoequivalente a um material de camada única quase-isotrópico.
Para a análise numérica considera-se que a placa está livre no espaço.Similarmente, para a análise modal experimental a placa foi suspensa com elásticos.
Na medição das frequências naturais do laminado se usou um vibrômetro laser,enquanto o estímulo foi dado através de um martelo de baixo impacto.
ARAÚJO et al. [55] realizaram uma análise modal experimental para identificaraté 20 frequências das quais 13 serão usadas na identificação de parâmetros.
Os resultados obtidos com os algoritmos FDIPA-FD e FDIPA-AG sãoapresentados na Tabela 8.2 e na Figura 8.4. Observa-se que o desempenho dosdois algoritmos foi similar, chegando à mesma solução.
Tabela 8.2: Resultado ótimo para a placa de laminado, parâmetros adimensionais.
Parâmetro FDIPA-FD FDIPA-AG
α0 0.8107 0.8207α2 3.6697 3.6585α3 0.8401 0.8402α4 1.0988 1.0972α8 0.2232 0.2241α9 0.0326 0.0324
Desvio Φ 9.0396× 10−5 9.0401× 10−5
NF 130 128
101
0 20 40 60 80 100 1200
0.05
0.1
0.15
0.2
0.25
0.3
0.35
NF: avaliações da função objetivo
f(x)
: val
or d
a fu
nção
obj
etiv
o
FDIPA−FDFDIPA−AG
Figura 8.4: Convergência dos algoritmos FDIPA-FD e FDIPA-AG no problema deidentificação de propriedades.
8.4 Distribuição ótima dos arranjos sensor-atuador piezoelétricos
Os materiais piezoelétricos possuem como característica fundamental acapacidade de conversão de energia elétrica em energia mecânica e vice-versa. As estruturas multilaminadas do tipo placa ou casca que integrammateriais piezoelétricos têm sido objeto de um crescente interesse em engenharia,especialmente no domínio das estruturas adaptativas, onde têm sido usadas comosensores e atuadores e, mais recentemente, em dispositivos de geração de energiaeléctrica [64–68].
As aplicações de dispositivos do tipo placa e casca multilaminados com materiaispiezoelétricos incluem o controle de deformação de asas de avião [64], estruturasde automóveis, antenas refletoras [65], espelhos deformáveis e atuadores [66].Adicionalmente, o uso de materiais piezoelétricos embebidos em estruturas com oobjetivo de produzir energia elétrica tem recebido cada vez mais atenção devidoà necessidade de eliminar ou reduzir fontes de energia externa ou baterias emdispositivos operados de forma remota [67, 68].
A distribuição do material piezoelétrico numa determinada camada afeta odesempenho da estrutura. Como consequência a quantidade, forma, tamanho,posição e polarização do material deve ser tido em conta simultaneamente numproblema de otimização. Neste sentido, a otimização de topologia tem sido aplicadaao projeto de transdutores piezoelétricos laminados do tipo placa e casca em busca
102
da topologia ótima. Assim, KOGL & SILVA [69] aplicaram a otimização de topologiaao projeto de atuadores piezoelétricos do tipo placa e casca. DONOSO & SIGMUND[70] usaram otimização de forma para determinar a espessura e forma de camadaspiezoelétricas numa viga, de modo a minimizar as cargas dinâmicas. DONOSO etal. [71] apresentaram métodos numéricos e analíticos para o projeto de sensorese atuadores modais em painéis cilíndricos. KANG & TONG [72, 73] usaramotimização de topologia para encontrar o empilhamento de camadas piezoelétricase estruturais, assim como as diferenças de potencial para atuação, com vistaao controle do campo de deslocamentos de placas multicamada piezoelétricas.No que respeita ao projeto de dispositivos de geração de energia baseados emestruturas do tipo placa e casca multilaminadas e usando otimização de topologia,ZHENG et al. [74] maximizaram a conversão de energia mecânica em energiaelétrica num dispositivo estático, enquanto RUPP [75, 76] considerou o projetodinâmico desses dispositivos. O problema de encontrar a distribuição ótima deum determinado número de pares piezoelétricos sensor/atuador pode ser formuladocomo um problema de otimização topológica. O problema de otimização consiste emdistribuir os sensores/actuadores piezoeléctricos na superfície da placa de modo aminimizar a massa da estrutura, especificando os níveis mínimos de amortecimentomodal para a gama de frequências de interesse. O amortecimento é introduzido naestrutura através do controle de velocidade, especificando um valor de ganho Gv
negativo na Eq. (8.12). A formulação usada neste trabalho é baseada no modeloPEMAP (piezoelectric material with penalization) [77], o qual é uma extensão domodelo SIMP (solid isotropic material with penalization), onde a pseudo-densidadedo material ρe de cada elemento de arranjo piezoelétrico é usada para atingir ummodelo de material continuo [69]:
K(e)uu = ρpke K0
(e)uu (8.19a)
K(e)uφ = ρpke K0
(e)uφ (8.19b)
M(e)uu = ρpme M0
(e)uu . (8.19c)
onde K0(e)uu , K0
(e)uφ e M0
(e)uu são respetivamente as matrizes de rigidez mecânica e
piezoelétrica e a matriz de massas de cada elemento na malha de elementos finitose K(e)
uu , K(e)uφ e M(e)
uu são as correspondentes matrizes penalizadas.Os fatores de penalidade pk e pm são usados para penalizar densidades
intermédias, facilitando assim a tarefa de obter uma solução discreta com valores0-1. Neste trabalho foram usados fatores de penalidade pk = 3 e pm = 1 para amatriz de rigidez e de massa respectivamente dos arranjos piezoelétricos. Note-se
103
que todas as matrizes elementares na Eq. (8.19) se referem unicamente aos arranjospiezoelétricos: a placa base nunca é afetada pelas variáveis de pseudo-densidade.
Uma primeira formulação apresentada em [78] tratou do posicionamento portopologia, de atuadores/sensores para a maximização dos fatores de amortecimentomodal, de uma placa sanduíche multilaminada, simplesmente apoiada, usando oalgoritmo FAIPA (Feasible Arc Interior Point Algorithm) [79, 80] e alternativamenteum algoritmo genético para efeitos comparativos. Embora os resultados foramsimilares, a abordagem por topologia usando FAIPA foi mais eficiente do que ouso de Algoritmos Genéticos com variáveis discretas.
A nova formulação proposta nesta tese, para o problema de otimização a serresolvido, consiste na minimização da massa dos sensores/atuadores com restriçõesnos fatores de amortecimento. Este problema pode ser formulado da seguinte forma:
minimizarρe
Ne∑e=1
ρe
sujeito a: ηm ≥ ηmin, m = 1, . . . , Nmodos
ρmin ≤ ρe ≤ 1,
(8.20)
onde Ne é o número de elementos da malha, ρe é a pseudo-densidade de cadaelemento piezoelétrico, e ηm são os fatores de amortecimento modal. Os fatoresde amortecimento modal ηm são obtidos resolvendo (8.11) iterativamente.
No problema que é apresentado, é usada uma placa de 300 mm × 200 mm,discretizada com uma malha de Ne = 24 elementos (6 × 4). A placa tem espessurade 2 mm e, por uma questão de simplicidade, é isotrópica (G = 80 GPa, ν = 0.3,ρ = 7860 kg/m3). Os elementos piezoelétricos têm uma espessura de 0.5 mm e terãoas dimensões de um dos elementos da malha (50 mm × 50 mm), e têm as seguintespropriedades: E1 = E2 = 47 GPa, G12 = 16.3 GPa, G23 = G13 = 17.4 GPa,ν = 0.33, ρ = 8036 kg/m3, e∗31 = e∗32 = −14.7 N/Vm e ε∗S33 = 2.12 × 10−8 F/m. Ovalor do ganho de controle de deslocamento é de Gd = 0, e do ganho de controle develocidade é Gv = −0.01. Considera-se que a placa está livre no espaço. A Figura8.5, (ver [55]), ilustra a placa sem piezoelétrico suspensa por elástico e tambémcom um arranjo de piezoelétrico, obtido após optimizar sua localização usando umAlgoritmo Genético.
104
Figura 8.5: Placa CFRP sem (esquerda) e com (direita) arranjos piezoelétricos.
São impostas restrições nos primeiros 6 fatores de amortecimento modais comηmin = 0.15 para os 6 modos. No ponto inicial as pseudo-densidades são ρ0
e = 0.999e toma-se ρmin = 0.1 por questões de estabilidade numérica.
Os resultados para a distribuição dos arranjos (sensor/atuador) piezoelétricosestão registrados na Tabela 8.3 e a distribuição correspondente se ilustra na Figura8.6. A Figura 8.7 mostra que o melhor desempenho se obtém usando derivadasanalíticas (FDIPA) seguido pelo desempenho obtido com aproximação das derivadasusando o algoritmo proposto (FDIPA-AG) e o pior desempenho foi obtido usandodiferenças finitas (FDIPA-FD).
Os três algoritmos chegaram a distribuições diferentes de 4 arranjos(sensor/atuador) piezoelétricos. O que indica a existência de mínimos locais nesteproblema.
105
Tabela 8.3: Resultado ótimo da distribuição dos arranjos de sensor/atuadorpiezoelétricos.
Densidade ρ FDIPA FDIPA-FD FDIPA-AG
ρ1 0.1000 0.1000 0.1003ρ2 0.1000 0.1000 0.1032ρ3 0.1000 0.1000 0.1003ρ4 0.1000 0.1000 0.1003ρ5 0.1000 0.1000 0.9300ρ6 0.1000 0.1000 0.1003ρ7 0.1138 0.1000 0.1001ρ8 0.1000 0.9999 0.1002ρ9 0.1000 0.1000 0.9978ρ10 0.9920 0.1000 0.9405ρ11 0.1000 0.8787 0.1002ρ12 0.9999 0.1000 0.1002ρ13 0.1000 0.1000 0.1003ρ14 0.9999 0.9694 0.9980ρ15 0.1000 0.1000 0.1003ρ16 0.9497 0.9999 0.1004ρ17 0.1000 0.1000 0.1037ρ18 0.1000 0.1000 0.1003ρ19 0.1000 0.1000 0.1003ρ20 0.1000 0.1000 0.1003ρ21 0.1000 0.1000 0.1003ρ22 0.1000 0.1000 0.1003ρ23 0.1000 0.1000 0.1003ρ24 0.1000 0.1000 0.100324∑i=1
ρi 5.9556 5.8482 5.8786
NF 108 1500 1288
106
(a) Resultados obtidos com FDIPA.
(b) Resultados obtidos com FDIPA-FD
(c) Resultados obtidos com FDIPA-AG
Figura 8.6: Distribuções ótimas dos arranjos piezoelétricos (sensor/atuador).
107
0 500 1000 15005
10
15
20
25
NF: avaliações da função objetivo
f(x)
: val
or d
a fu
nção
obj
etiv
o
FDIPAFDIPA−FDFDIPA−AG
Figura 8.7: Convergência dos algoritmos FDIPA, FDIPA-FD e FDIPA-AG noproblema de distribuição ótima de arranjos sensor/atuador piezoelétricos.
108
Capítulo 9
Considerações finais
9.1 Conclusões
Neste trabalho foi desenvolvido uma nova técnica de otimização restrita comaproximação de derivadas. Este método foi baseado numa família de algoritmosconhecida como algoritmos de direções viáveis, especificamente num algoritmo tipoquase-Newton de direções viáveis e ponto interior (FDIPA). Como foi mostrado,por causa das características da construção das aproximações, o método propostoé especialmente indicado para problemas de otimização de pequeno porte emEngenharia. Para ilustrar o desempenho do método de otimização proposto foramestudados os problemas de otimização de estruturas reticuladas e de otimização deestruturas em materiais compósitos laminados.
No Capítulo 5, Seção 5.3, foi apresentada uma primeira versão do algoritmoproposto, aqui chamado FDIPA-WLS-AG, que baseia-se no algoritmo de otimizaçãoFDIPA. O algoritmo FDIPA-WLS-AG incorpora aproximações não interpolantesnum esquema acumulativo de pontos visando o aproveitamento dos dados obtidosao longo do processo de otimização. Nesta versão não se exerce controle nocondicionamento do conjunto de dados amostrais, o que leva o algoritmo a realizarpassos auxiliares usando diferenças finitas toda vez que a qualidade da aproximaçãoseja prejudicada.
Posteriormente na Seção 5.4 foi apresentada uma versão mais robustado algoritmo proposto, denominado aqui FDIPA-AG, que usa aproximaçõesinterpolantes para a obtenção dos gradientes aproximados. Tais aproximaçõessão atualizadas quando necessário visando a redução do custo computacional. Arobustez vem da convergência garantida do algoritmo base (FDIPA) junto com umcuidadoso esquema de controle da qualidade nas aproximações do gradiente.
O algoritmo FDIPA-AG mostrou-se eficaz na solução de problemas de otimizaçãonão linear com restrições. Foi estudado seu desempenho em aplicações de otimização
109
estrutural sob restrições mecânicas (estáticas e dinâmicas). Seu desempenho erobustez foi maior do que a obtida por um dos algoritmos mais usados na área(algoritmos genéticos) e resultou comparável ao conhecido método de diferençasfinitas.
No Capítulo 7 foram apresentadas duas estratégias desenvolvidas paraproporcionar uma maior tolerância ao ruído externo, sendo esta uma característicaimportante necessária nas aplicações práticas da otimização na ausência dederivadas. A modificação do algoritmo proposto para o caso de ruído nas funçõesficou consolidada na Seção 7.3 sob o nome de FDIPA-AG-N, a qual mostrouresultados satisfatórios e um melhor desempenho na presença de ruído moderado.Cabe destacar que na presença de ruído moderado o esquema de diferenças finitasse mostrou inaplicável.
Por último, duas aplicações do método proposto em estruturas de materiaiscompósitos laminados foram apresentadas no Capítulo 8. No primeiro problemase realiza a identificação de propriedades mecânicas do compósito laminado. Nosegundo problema, se apresenta uma nova formulação da distribuição ótima deatuadores piezoelétricos para a dissipação de vibrações. Nos dois problemaso algoritmo FDIPA-AG se mostrou competitivo em comparação com o uso dediferenças finitas FDIPA-FD.
9.2 Sugestões para trabalhos futuros
Propõem-se como parte de estudos futuros os seguintes tópicos:.
? Estudar estratégias para abordar problemas de maior porte, e.g. filtragem derestrições.
? Nenhum estudo sobre análise de complexidade foi realizado, o que pode sugerirfuturas melhoras no desempenho do algoritmo.
? Explorar mais problemas que possam ser resolvidos com o algoritmo aquiapresentado, principalmente aqueles que apresentam ruído moderado, e.g.(CFD, Computational Fluid Dynamics).
? Estudar estratégias para obter um mínimo global no problema deposicionamento ótimo de piezoelétricos para a dissipação de vibrações.Propõe-se o uso do algoritmo apresentado neste trabalho, partindo de pontosiniciais viáveis diferentes que representem toda a região viável.
110
Referências Bibliográficas
[1] HERSKOVITS, J. “Feasible direction interior-point technique for nonlinearoptimization”, J. Optim. Theory Appl., v. 99, n. 1, pp. 121–146, 1998.
[2] GRIEWANK, A. Evaluating Derivatives: Principles and Techniques ofAlgorithmic Differentiation. Philadelphia, SIAM, 2000.
[3] GRIEWANK, A., CORLISS, G. F. Automatic Differentiation of Algorithms:Theory, Implementation, and Application. Philadelphia, SIAM, 1991.
[4] POWELL, M. J. D. The NEWUOA software for unconstrained optimizationwithout derivatives. DAMTP 2004/NA08, University of Cambridge, 2004.
[5] SACKS, J., WELCH, W. J., MITCHELL, T. J., et al. “Design and analysis ofcomputer experiments”, Statist. Science, v. 4, n. 4, pp. 409–423, 1989.
[6] KIRKPATRICK, S., GELATT, C. D., VECCHI, M. P. “Optimization bysimulated annealing”, Science, v. 220, pp. 671–680, 1983.
[7] GOLDBERG, D. E. Genetic Algorithms in Search, Optimization and MachineLearning. Boston, Addison-Wesley Longman, 1989.
[8] HOLLAND, J. H. Adaptation in Natural and Artificial Systems: IntroductoryAnalysis with Applications to Biology, Control, and Artificial Intelligence.Cambridge, MIT Press, 1992.
[9] HAYKIN, S. Neural Networks: A Comprehensive Foundation. Upper SaddleRiver, Prentice-Hall, 1994.
[10] GLOVER, F. “Tabu search – Part I”, ORSA J. Comput., v. 1, pp. 190–206,1989.
[11] KENNEDY, J., EBERHART, R. “Particle swarm optimization”. In:Proceedings of the 1995 IEEE International Conference on NeuralNetworks, pp. 1942–1948, Perth, 1995. IEEE Service Center.
111
[12] LEWIS, R. M., TORCZON, V., TROSSET, M. W. “Direct Search Methods:Then and Now”, Journal of Computational and Applied Mathematics,v. 124, pp. 191–207, 2000.
[13] HOOKE, R., JEEVES, T. A. “Direct search solution of numerical and statisticalproblems”, Journal of the ACM, v. 8, pp. 212–229, 1961.
[14] TORCZON, V. “On the convergence of pattern search algorithms”, SIAMJournal on Optimization, v. 7, n. 1, pp. 1–25, 1997.
[15] SPENDLEY, W., HEXT, G. R., HIMSWORTH, F. R. “Sequential applicationof simplex designs in optimization and evolutionary operation”,Technometrics 4, v. 4, n. 4, pp. 441–461, 1962.
[16] NELDER, J. A., MEAD, R. “A simplex method for function minimization”,Computer Journal, v. 7, pp. 308–313, 1965.
[17] MCKINNON, K. I. “Convergence of the Nelder-Mead simplex method to anonstationary point”, SIOPT, v. 9, n. 1, pp. 14–158, 1998.
[18] ROSENBROCK, H. “An automatic method for finding the greatest or leastvalue of a function”, Computer Journal, v. 3, pp. 175–184, 1960.
[19] POWELL, M. J. D. “An efficient method for finding the minimum of a functionof several variables without calculating derivatives”, Computer Journal,v. 17, pp. 155–162, 1964.
[20] CONN, A. R., TOINT, P. L. “An algorithm using quadratic interpolation forunconstrained derivative free optimization”. In: Pillo, G. D., Gianessi, F.(Eds.), Nonlinear Optimization and Applications, Plenum Publishing, pp.27–47, New York, 1996.
[21] POWELL, M. J. D. Uobyqa: unconstrained optimization by quadraticapproximation. Relatório Técnico DAMTP NA14, Department of AppliedMathematics and Theoretical Physics, Cambridge University, Cambridge,2000.
[22] POWELL, M. J. D. On the use of quadratic models in unconstrainedminimzation without derivatives. Relatório Técnico DAMTP NA03,Cambridge University, Cambridge, 2003.
[23] BOOKER, A. J., DENNIS, J. E., FRANK, P. D., et al. “A rigourous frameworkfor optimization of expensive functios by surrogates”, Structural andMultidisciplinary Optimization, v. 17, n. 1, pp. 1–13, 1999.
112
[24] BONABEAU, E., DORIGO, M., THERAULAZ, G. Swarm intelligence: fromnatural to artificial systems. New York, Oxford University Press, 1999.
[25] DORIGO, M., MANIEZZO, V., COLORNI, A. “The ant system: optimizationby a colony of cooperating agents”, IEEE T. on Syst. Man Cy. B, v. 26,n. 2, pp. 29–41, 1996.
[26] EBERHART, R., SHI, Y. “Particle swarm optimization: developments,applications and resources”. In: Proceedings of Congress on EvolutionaryComputation, pp. 81–86, Seoul, 2001. IEEE.
[27] ABBASS, H. A. “MBO: Marriage in Honey Bees Optimization, AHaplometrosis Polygynous Swarming Approach”. In: Proceedings of theIEEE Congress on Evolutionary Computation, pp. 207–214. IEEE, 2001.
[28] KRISHNANAND, K. N., GHOSE, D. “Glowworm swarm optimization forsimultaneous capture of multiple local optima of multimodal functions”,Swarm Intelligence, v. 3, n. 2, pp. 87–124, 2009.
[29] CONN, A. R., SCHEINBERG, K., VICENTE, L. N. “Global convergence ofgeneral derivative-free trust-region algorithms to first and second ordercritical points”, SIAM Journal on Optimization, v. 20, pp. 387–415, 2009.
[30] GUJARATI, D. N. Basic Econometrics. 3 ed. Singapore, McGraw-HillInternational Editions, 1995.
[31] OEUVRAY, R. Trust-Region Methods Based on Radial Basis Functionswith Application to Biomedical Imaging. Tese de Doutorado, ÉcolePolytechnique Fédérale de Lausanne, 2005.
[32] WILD, S. M., REGIS, R. G., SHOEMAKER, C. A. “ORBIT: Optimizationby radial basis function interpolation in trust-regions”, SIAM J. Sci.Comput., v. 30, n. 6, pp. 3197–3219, 2008.
[33] BUHMANN, M. D. Radial Basis Functions. Cambridge, Cambridge UniversityPress, 2003.
[34] MATHERON, G. “Principles of geostatistics”, Econ. Geol., v. 58, pp. 1246–1266, 1963.
[35] CURRIN, C., MITCHELL, T., MORRIS, M., et al. A Bayesian approachto the design and analysis of computer experiments. Relatório TécnicoORNL-6498, Oak Ridge National Laboratory, Oak Ridge, 1988.
113
[36] MYERS, R. H., MONTGOMERY, D. C. Response Surface Methodology:Process and Product Optimization Usign Designed Experiments. NewYork, John Wiley, 1995.
[37] SIMPSON, T. W. A Concept Exploration Method for Product FamilyDesign. Tese de Doutorado, George W. Woodruff School of MechanicalEngineering, Georgia Institute of Technology, Atlanta, GA, 1998.
[38] MCKAY, M., CONOVER, W., BECKMAN, R. “A comparison of threemethods for selecting values of input variables in the analysis of outputfrom a computer code”’, Technometrics, v. 21, pp. 239–245, 1979.
[39] HEATH, M. T. Scientific Computing: An Introductory Survey. McGraw-Hill,1997.
[40] BAZARAA, M. S., SHERALI, H. D., SHETTY, C. M. Nonlinear programming.3 ed. Hoboken, Wiley-Interscience, 2006.
[41] LUENBERGER, D. G., YE, Y. Linear and Nonlinear Programming. Springer,2008.
[42] CORTÉS, H., HERSKOVITS, J. “A Globally Convergent Technique forNonlinear Constrained Optimization With Surrogates”. In: Proceedings ofthe 8th World Congress on Structural and Multidisciplinary Optimization,June 1-5 2009.
[43] WUJEK, B. A., RENAUD, J. E. “New adaptive move-limit managementstrategy for approximate optimization”, AIAA J., , n. 36, pp. 1911–1934,1998.
[44] ELDRED, M. S., DUNLAVY, D. M. “Formulations for Surrogate-BasedOptimization with Data Fit, Multifidelity, and Reduced-Order Models”.In: Proc. 11th AIAA/ISSMO Multidisciplinary Analysis and OptimizationConf., AIAA Paper 2006-7117, Portsmouth, Virginia, september 2006.AIAA.
[45] SCHITTKOWSKI, K. Test Examples for Nonlinear Programming Codes.Springer-Verlag, 1980.
[46] STONEKING, D., BILBRO, G., TREW, R., et al. “Yield optimization using aGaAs process simulator coupled to a physical device model”, IEEE Trans.Microwave Theory and Techniques, v. 40, pp. 1353–1363, 1992.
114
[47] BORTZ, D. M., KELLEY, C. T. “The simplex gradient and noisy optimizationproblems”. In: Borggaard, J. T., Burns, J., Cliff, E., et al. (Eds.),Computational Methods in Optimal Design and Control, Birkhauser, pp.77–90, Boston, 1998.
[48] ZAGAINOV, G. I., LOZINO-LOZINSKI, G. E. Composite materials inaerospace design. Chapman & Hall, 1996.
[49] GÜRDAL, Z., HAFTKA, R. T., HAJELA, P. Design and optimization oflaminated composite materials. John Wiley & Sons, 1999.
[50] CORREIA, V. M. F., SOARES, C. M. M., SOARES, C. A. M. “Bucklingoptimization of composite laminated adaptive structures”, CompositeStructures, v. 62, pp. 315–321, 2003.
[51] MOITA, J. M., CORREIA, V. M. F., MARTINS, P. G., et al. “Optimal designin vibration control of adaptive structures using a simulated annealingalgorithm”, Composite Structures, v. 75, pp. 79–87, 2006.
[52] SULEMAN, A., GONÇALVES, M. A. “Multi-objective optimization of anadaptive composite beam using the physical programming approach”,Journal of Intelligent Material Systems and Structures, v. 10, pp. 56–70,1999.
[53] ADALI, S. J., BRUCH, J. C., SADEK, I. S., et al. “Robust shape controlof beams with load uncertainties by optimally placed piezo actuators”,Structural and Multidisciplinary Optimization, v. 19, pp. 274–281, 2000.
[54] FRECKER, M. I. “Recent advances in optimization of smart structures andactuators”, Journal of Intelligent Material Systems and Structures, v. 14,pp. 207–216, 2003.
[55] ARAÚJO, A. L., SOARES, C. M. M., FRIEDMANN, H., et al. “Optimallocation of piezoelectric patches and identification of material properties inlaminated composite structures”. In: International Conference on SmartStructures and Materials, Porto, July 13-15 2009.
[56] ARAÚJO, A. L., SOARES, C. M. M., SOARES, C. A. M., et al. “Developmentof a finite element model for the identification of mechanical andpiezoelectric properties through gradient optimisation and experimentalvibration data”, Composite Structures, v. 58, pp. 307–318, 2002.
115
[57] ARAÚJO, A. L. Identification of Elastic and Piezoelectric Properties ofAdaptive Materials. Tese de Doutorado, Universidade Técnica de Lisboa,2007.
[58] ARAÚJO, A. L., SOARES, C. M. M., SOARES, C. A. M. “Finite elementmodel for hybrid active-passive damping analysis of anisotropic laminatedsandwich structures”, Journal of Sandwich Structures and Materials,v. 12, n. 4, pp. 397–419, 2010.
[59] ARAÚJO, A. L., SOARES, C. M. M., SOARES, C. A. M. “A ViscoelasticSandwich Finite Element Model for the Analysis of Passive, Active andHybrid Structures”, Applied Composite Materials, v. 17, pp. 529–542,2010.
[60] LEMPRIERE, B. M. “Poissonś ratio in orthotropic materials”, AIAA Journal,pp. 2226–2227, 1968.
[61] JONES, R. M. Mechanics of Composite Materials. New York, McGraw-Hill,1975.
[62] ARAÚJO, A. L., SOARES, C. M. M., SOARES, C. A. M., et al.“Characterisation by inverse Techniques of Elastic, Viscoelastic andPiezoelectric Properties of Anisotropic Sandwich Adaptive Structures”,Applied Composite Materials, v. 17, pp. 543–556, 2010.
[63] ARAÚJO, A. L., LOPES, H. M. R., VAZ, M. A. P., et al. “Parameter estimationin active plate structures”, Computers and Structures, v. 84, pp. 1471–1479, 2006.
[64] TZOU, H. S., TSENG, C. I. “Distributed piezoelectric sensor/actuator designfor dynamic measurement/control of distributed parameter systems: Apiezoelectric finite element approach”, Journal of Sound and Vibration,v. 138, n. 1, pp. 17–34, 1990.
[65] AGRAWAL, B. N., TREANOR, K. E. “Shape control of a beam usingpiezoelectric actuators”, Smart Materials and Structures, v. 8, n. 6,pp. 729–739, 1999.
[66] MUKHERJEE, A., JOSHI, S. “Piezoelectric sensor and actuator spatial designfor shape control of piezolaminated plates”, AIAA Journal, v. 40, n. 6,pp. 1204–1210, 2002.
116
[67] SODANO, H. A., INMAN, D. J., PARK, G. “Comparison of piezoelectricenergy harvesting devices for recharging batteries”, Journal of IntelligentMaterial Systems and Structures, v. 16, n. 10, pp. 799–807, 2005.
[68] PARK, K. I., XU, S., LIU, Y., et al. “Piezoelectric BaTiO3 Thin FilmNanogenerator on Plastic Substrates”, Nano Letters, v. 10, n. 12,pp. 4939–4943, 2010.
[69] KÖGL, M., SILVA, E. C. N. “Topology optimization of smart structures:design of piezoelectric plates and shell actuators”, Smart Materials andStructures, v. 14, pp. 387–399, 2005.
[70] DONOSO, A., SIGMUND, O. “Optimization of piezoelectric bimorph actuatorswith active damping for static and dynamic loads”, Structural andMultidisciplinary Optimization, v. 38, n. 2, pp. 171–183, 2009.
[71] DONOSO, A., BELLIDO, J., CHACN, J. “Numerical and analytical methodfor the design of piezoelectric modal sensors/actuators for shell-typestructures”, International Journal for Numerical Methods in Engineering,v. 81, pp. 1700–1712, 2010.
[72] KANG, Z., TONG, L. Y. “Topology optimization-based distribution designof actuation voltage in static shape control of plates”, Computers andStructures, v. 86, n. 19-20, pp. 1885–1893, 2008.
[73] KANG, Z., TONG, L. Y. “Integrated optimization of material layout andcontrol voltage for piezoelectric laminated plates”, Journal of IntelligentMaterial Systems and Structures, v. 19, n. 8, pp. 889–904, 2008.
[74] ZHENG, B., CHANG, C., GEA, H. “Topology optimization of energyharvesting devices using piezoelectric materials”, Structural andMultidisciplinary Optimization, v. 38, pp. 17–23, 2009.
[75] RUPP, C. J., EVGRAFOV, A., MAUTE, K., et al. “Design of PiezoelectricEnergy Harvesting Systems: A Topology Optimization Approach Basedon Multilayer Plates and Shells”, Journal of Intelligent Material Systemsand Structures, v. 20, n. 16, pp. 1923–1939, 2009.
[76] RUPP, C., DUNN, M., MAUTE, K. “Analysis of Piezoelectric EnergyHarvesting Systems with Non-linear Circuits Using the Harmonic BalanceMethod”, Journal of Intelligent Material Systems and Structures, v. 21,pp. 1383–1396, 2010.
117
[77] SILVA, E. C. N., KIKUCHI, N. “Design of piezoelectric transducers usingtopology optimization”, Smart Materials and Structures, v. 8, pp. 350–364, 1999.
[78] ARAÚJO, A. L., SOARES, C. M. M., SOARES, C. A. M., et al. “Optimalpositioning of piezoelectric patches in sandwich structures for maximumdamping”. In: Proceedings of the 5th ECCOMAS Thematic Conference onSmart Structures and Materials SMART11, July 6-8 2011.
[79] HERSKOVITS, J., SANTOS, G. “Feasible arc interior point algorithm fornonlinear optimization”. In: Proceedings of the 4th World Congress onComputational Mechanics - New Trends and Applications, Barcelona,1998. CIMNE.
[80] HERSKOVITS, J., MAPPA, P., GOULART, E., et al. “Mathematicalprogramming models and algorithms for engineering design optimization”,Computer Methods in Applied Mechanics and Engineering, v. 194,pp. 3244–3268, 2005.
118
![[Alexandre] 8. Não Linear Restrita](https://img.document.onl/doc/110x75/557604d9d8b42a81158b518a/alexandre-8-nao-linear-restrita.jpg)
