Aplicações de um Modelo Substituto de Ordem Reduzida

Universidade Federal de Pernambuco Centro de Tecnologia ecnolog ia e Geociˆ Geo ciências encias Departamento de Engenharia Civil

Aplica¸c˜ coes o ˜es de um modelo substituto de ordem reduzida a estudos de gerenciamento de reservat´ orios orios de petr´ oleo oleo

Manuel Fragoso Machado Junior

Recife/2014

Universidade Federal de Pernambuco Centro de Tecnologia ecnolog ia e Geociˆ Geo ciências encias Departamento de Engenharia Civil

Manuel Fragoso Machado Junior

Aplica¸c˜ coes o ˜es de um modelo substituto de ordem reduzida a estudos de gerenciamento de reservat´ orios orios de petr´ oleo oleo

Disserta¸ c˜ cao a õ de Mestrado apresentada apresentada ao Programa de P´ os-gradua¸ os-gradua¸ c˜ c˜ ao ao do Departamento Departa mento de Engenharia Civil da Universidade Federal de Pernambuco como requisito requisito parcial parcial para obten¸ obten¸ c˜ cao a õ do grau de Mestre na ´ area area de Geren Gerencia ciamen mento to de Reserv Reservat´ at´ orio o rioss de Petr´ oleo. oleo.

Orientado Orientador: r: Prof. Dr. Bernardo Bernardo Horowitz Horowitz Co-Orientador: Co-Orientador: Dr. Jos´ e Roberto Pereira Rodrigues

Recife/2014

Catalogação na fonte Bibliotecária Margareth Malta, CRB-4 / 1198

M149a

Machado Junior, Manuel Fragoso. Aplicações de um modelo substituto de ordem reduzida a estudos de gerenciame gerenciamento nto de reservatór reservatórios ios de petróleo petróleo / Manuel Manuel Fragoso Fragoso Machado Machado Junior. - Recife: O Autor, 2014. 2014. 135 folhas, folhas, il., gráfs., gráfs., tabs. tabs. Orientador: Prof. Dr. Bernardo Horowitz. Coorientador: Prof. Dr. José Roberto Pereira Rodrigues. Dissertaç Dissertação ão (Mestrado (Mestrado)) – Universida Universidade de Federal Federal de Pernambu Pernambuco. co. CTG. CTG. Programa de Pós-Graduação em Engenharia Civil, 2014. Inclui Referências. Engenharia Civil. 2. Modelo Modelo de ordem reduzida. reduzida. 3. Trajectory Piecewise Linearization Linearization.. 4. Simulação de reservatórios. 5. Otimização sequencia sequenciall aproximad aproximadaa multifide multifidelidad lidade. e. 6. Control variates. variates. 7. Propagação de Incertezas. I. Horowitz, Bernardo. (Orientador). II. Rodrigues, José Roberto Pereira. (Coorientador). III. Título. 1.

624 CDD (22. ed.)

UFPE BCTG/2012-304

UNIVERSIDADE FEDERAL DE PERNAMBUCO ´ PROGRAMA DE POS-GRADUAC ¸ ˜ AO EM ENGENHARIA CIVIL A comissão examinadora da defesa da disserta¸caõ de mestrado APLICAC ¸ ˜ OES DE UM MODELO SUBSTITUTO DE ORDEM REDUZIDA A ´ ´ ESTUDOS DE GERENCIAMENTO DE RESERVATORIOS DE PETROLEO

defendida por Manuel

Fragoso Machado Junior

considera o candidato APROVADO.

Recife, 28 de maio de 2014

Orientadores: Prof. Dr. Bernardo Horowitz (Orientador) Dr. José Roberto Pereira Rodrigues (Co-Orientador) Banca Examinadora: Prof. Dr. Paulo Goldfeld — UFRJ (Examinador Externo) Dr. Marco Antônio Cardoso — Petrobras (Examinador Externo) Prof. Dr. Ramiro Willmersdorf - UFPE (Examinador Interno) Profa. Dra. Silvana Maria Bastos da Silva (Examinadora Interna)

´ DEDICATORIA

Dedico este trabalho aos meus queridos amigos e familiares e em especial à minha noiva, Elaine, e à minha mãe, Maria Luiza.

AGRADECIMENTOS Muitas foram as pessoas que contribuiram para a realiza¸cão deste trabalho. São tantas que tenho receio de me esquecer de alguém. Mas vou iniciar minha longa lista de agradecimentos. ´ Ele quem me permite acordar todos os Antes de tudo quero agradecer a Deus. E dias, seguir minha jornada, fazer minhas escolhas. Sem Ele não teria for¸cas para realizar este trabalho. Quero agradecer a meu orientador Prof. Bernardo Horowitz por toda a sua boa vontade em orientar, suas ideias valiosas e as conversas durante o almo¸co ou café. Desde o primeiro dia em que cheguei a Recife para este trabalho me senti direcionado a desenvolver um trabalho interessante e desafiante. ´ Agrade¸co também ao Prof. Ezio Araujo por todas as conversas técnicas ou não, pelo acolhimento no grupo de pesquisa e por emprestar-me um carro que foi útil em diversos momentos. Também não posso me esquecer do Prof. Leonardo Guimar˜ aes, pelas conversas cheias de ideias, pelas pedaladas no fim de semana, pelas excelentes aulas e pela amizade que foi cultivada. Pelas ótimas aulas e trocas de ideias devo agradecer tamb´ em ao Prof. Ramiro Wilmersdorf e Profa. Silvana Bastos. Dois dos cap´ıtulos desta disserta¸caõ surgiram em trabalhos de suas disciplinas. Também não posso esquecer dos Profs. Darlan Karlo e Antônio pelas aulas e conversas de corredor sobre assuntos do mundo do petróleo ou não. Entre os alunos do Programa de Pós-gradua¸caõ em Engenharia Civil não posso deixar de mencionar alguns nomes: Liliane Fonseca, pelo aux´ılio na minha chegada a Recife; Jefferson Wellano e Julio pelos bons momentos de descontra¸cão no laboratório de informática; Roberto Pareja e Rodrigo pelo companheirismo nas disciplinas que cursamos juntos.

Muito obrigado também a`s funcionárias do Programa Andrea e Rose Mary por resolverem todos os tipos de problemas que tive com a burocracia da universidade. Também tenho muito a agradecer à empresa em que trabalho, a PETROBRAS, por me proporcionar a oportunidade de realizar esta pesquisa e pelo apoio financeiro, pagando o meu salá rio. Em especial, não posso esquecer de minha gerente, Flávia Pacheco, pela permissão e est´ımulo. Tamb´ em agrade¸co aos Engos . Régis Kruel, por indicar-me a UFPE, Marco Antônio e José Roberto, por responderem às questões sobre TPWL e SIMPAR com toda a paciência e ao Leonardo Cabral, que fazendo doutorado na mesma universidade, me proporcionou o´timas conversas. Agrade¸co também aos meus amigos de sala no CENPES, Engos . Guilherme Teixeira e Lu´ıs Carlos e Rodrigo Cárpio pelas conversas produtivas. Além das pessoas da universidade, muitos foram aqueles que me proporcionaram ótimos momentos de descontra¸caõ durante a estada em Recife. Por causa deles, o trabalho árduo e a distância da fam´ılia não foram grandes problemas. Cito aqui os principais nomes que deixarão saudades, mas a lista completa é bem extensa. Agrade¸co então a todos os meus amigos que moravam na Rep´ ublica das Palmeiras em especial a Ana Lúcia, Ciancony Rocha, Cristina Nepomuceno, Daniel Lazo, Elisa Griti, Fagunes Ferreira, Manuela, Mar´ılia Santos, Murilo Ara´ ujo, Pedro Souza, Ranielder e Vanessa Pedrosa. Tenho muito a agradecer também aos meus familiares que me apoiaram, mesmo que a distância. Agrade¸co a meu irmão, Carlos Augusto, que passou comigo ótimos momentos em Recife e sempre me motivou na busca pelo conhecimento, a meu irmão Marcelo Machado e meu pai, Manuel Fragoso, pela ajuda constante com todo tipo de problema que pudesse ter em minha terra natal, aos meus padrinhos, Heloisa Moreira e Jorge Guilherme por tudo que fizeram por mim ao longo de minha vida, a Sr a . Maria Jos´ e, por ter sido sempre a minha segunda m˜ ae e a meus cumpadres, Renato Cabral e Luciana Ferrari, por me darem a honra de batizar sua filha, Larinha. Tamb´ em agrade¸co a minha sobrinha e afilhada Ma´ıssa Machado, por todo o carinho que me dispunha em todas as vezes que voltei ao Rio. Finalmente agrade¸co muito à minha noiva Elaine Andrade, pelo amor, carinho e paciência que a mim ela dedica. Também à minha mãe, que não só me deu a vida, mas me dá todo o apoio em todas as minhas escolhas.

“A imagina¸caõ é mais importante que a ciência, porque a ciência é limitada, ao passo que a imagina¸caõ abrange o mundo inteiro.” (Albert Einstein)

RESUMO Esta disserta¸caõ tem como objetivo a aplica¸caõ de esquemas multifidelidade de otimiza¸caõ e propaga¸cão de incertezas ao gerenciamento de reservatórios de petróleo. A metodologia Trajectory Piecewise Linearization -TPWL em combina¸caõ com Proper Orthogonal Decomposition - POD foi utilizada como modelo substituto de ordem reduzida para o simulador de fluxo. Estas técnicas reduzem a complexidade e a dimensão do problema atrav´ es da lineariza¸caõ de suas equa¸co˜es governantes em torno de estados convergidos e armazenados durante uma simula¸cão de treinamento. O método mostra-se acurado nas vizinhan¸cas da trajetória de treinamento. Ganho de tempo consider´ avel foi obtido pela sua utiliza¸cão, o que ficou demonstrado nas aplica¸co˜es realizadas com modelos de dois reservatórios diferentes. O simulador de reservatórios utilizado foi desenvolvido pela Petrobras e cedido às universidades participantes da rede SIGER para realizarem seus estudos de simula¸caõ. Consiste em um simulador black-oil , mas sua formula¸cão difere da normalmente utilizada em softwares comerciais, pois em vez de utilizar as satura¸co˜es como variáveis primárias utiliza fra¸cões mássicas. Seu código foi alterado de modo a realizar a exporta¸caõ dos dados necessários a` constru¸caõ do modelo substituto. Para comunica¸caõ entre o simulador e as rotinas de treinamento do modelo de ordem reduzida, utilizou-se um padrão hierárquico e genérico de arquivos binários chamado Hierachical Data Format . O problema do retreinamento do modelo de ordem reduzida foi estudado em mais detalhe, propondo-se um método para agregar os estados e derivadas de uma dada simula¸caõ ao conjunto já utilizado na constru¸caõ inicial. Além disso, estudou-se a necessidade de retreinamento durante o processo de otimiza¸cão, propondo-se um critério para realizá-lo a` medida que o processo de otimiza¸caõ avan¸ca, passando por pontos que produzem trajetórias cada vez mais distantes da utilizada no treinamento. O trabalho mostra que o problema com a instabilidade do TPWL/POD reportado na literatura é mitigado trocando-se a proje¸caõ de Bubnov-Galerkin, utilizada nos primeiros trabalhos com este método, pela proje¸caõ de Petrov-Galerkin. Isto já havia sido mencionado na literatura como mecanismo de estabiliza¸caõ. O estudo de otimiza¸caõ visou a maximiza¸caõ do valor presente

l´ıquido considerando pressões de fundo dos po¸cos como variáveis de projeto. Foi desenvolvido um algoritmo de Otimiza¸caõ Sequencial Aproximada com Multifidelidade (MSAO), que converte o problema de otimiza¸caõ em uma sequência de problemas aproximados, v´ alidos dentro de uma região de confian¸ca. Este algoritmo utiliza uma combina¸caõ entre o modelo de ordem reduzida aqui estudado e um modelo de krigagem, baseado em algumas simula¸co˜es realizadas no interior de cada região de confian¸ca. Bons resultados foram obtidos, tendo o algoritmo em alguns casos superado a referência que consiste na otimiza¸cão através da aplia¸caõ direta do algoritmo Sequential Quadratic Programing - SQP sobre o modelo de alta fidelidade. O estudo de incertezas consistiu no cálculo da média e do desvio padrão do valor presente l´ıquido considerando-se distribui¸cões de probabilidades para as pressões de fundo dos po¸cos. Um algoritmo Monte Carlo Multifidelidade foi desenvolvido baseado na técnica estat´ıstica das variáveis de controle. Este método tira proveito da correla¸cão entre os modelos de alta e baixa fidelidade, lan¸cando muitas simula¸co˜es do modelo de ordem reduzida, corrigindo seus resultados com uma quantidade menor de execu¸cões do simulador. Com base nesta correla¸caõ e na razão do custo computacional para simula¸caõ de alta e baixa fidelidade, o algoritmo determina a razã o do n´ umero de execu¸co˜es dos dois modelos.

Palavras-Chave: Modelo de ordem reduzida. Trajectory piecewise linearization . Simula¸caõ de reservatórios. Otimiza¸caõ sequencial aproximada multifidelidade. Control variates . Propaga¸caõ de incertezas.

ABSTRACT The goal of this dissertation is the application of optimization and uncertainty propagation based on multifidelity frameworks to oil reservoir managment. The reduced order surrogate model used to replace the flow simulator is based on the Trajectory Piecewise Linearization (TPWL) methodology combined with Proper Orthogonal Decomposition (POD). These techniques reduce the numerical complexity and dimension of the problem by performing the linearization of governing equations around converged states stored during a training simulation. The method is shown to be accurate in the neighbourhood of the training trajectory. Good speedups were achieved by its application to two different reservoir models. The reservoir simulator used in this work was developed by Petrobras and was given to Brazilian universities members of SIGER Research Network to perform simulation studies. It is a blackoil simulator but its formulation is different from what is used in commercial software. It uses mass fractions instead of saturations. The simulator code, in general, was changed to export the necessary data to build the reduced order model. Hierarchical data format was used to link the simulator to reduced order model. The TPWL/POD retraining process was studied in detail. It consists of a method to add new simulation states and derivatives to the training set built initially. We proposed a criterion to determine when retraining is needed during the optimization process. This work shows that the problem with instability of TPWL/POD reported in the literature is avoided by replacing the Bubnov-Galerkin projection with the Petrov-Galerkin one. This is not something new, but is confirmed by the results obtained. The optimization study developed in this work aims to maximize the net present value (NPV) using bottom hole pressures (BHP’s) as design variables. We developed a Multifidelity Sequential Approximate Optimization algorithm (MSAO) which transforms the real problem into a sequence of approximate problems, accurately represented inside a trust region. This approximation relies on a metamodel which is a combination of TPWL/POD and a kriging correction model. Good results were obtained, sometimes better than traditional optimization based on high fidelity model. The uncertainty propagation study aims to cal-

culate the average and standard deviation of NPV considering some probability distribution for the BHP’s. A Multifidelity Monte Carlo Algorithm was developed using Control Variates Tecnhnique. It takes advantage of the correlation between high and low fidelity models and launch many low fidelity runs, correcting the results with fewer high fidelity runs. Based on that correlation and computational cost rate, the algorithm determines the optimal number of high and low fidelity model evaluations.

Key-Words: Reduced order model. Trajectory piecewise linearization. Reservoir simulation. Multifidelity sequential approximation optimization. Control variates. Uncertainty propagation.

Lista de Figuras 2.1

Volume de controle e seus vizinhos . . . . . . . . . . . . . . . . . . . . . . . .

2.2

Modelo sintético (baseado no caso SPE10) com 24000 células, 5 po¸cos produtores e 2 injetores. Mapa de permeabilidades na dire¸cão x. . . . . . . . . . . . . . .

31 44

2.3

Compara¸cão de resultados do SIMPAR com IMEX - Modelo SPE10 - Modificado 45

2.4

Modelo de Brugge Modificado. Mapa de permeabilidades na dire¸caõ x. . . . .

47

2.5

Compara¸cã o de resultados do SIMPAR com IMEX . . . . . . . . . . . . . . .

48

3.1

Região Validade do TPWL em torno da trajetória de treinamento . . . . . . .

55

3.2

Funcionamento do POD - Consideram-se somente as dire¸cões de maior variabilidade . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3.3

63

Sequência de controles de treinamento (pressões de fundo) aplicados aos po¸cos produtores . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

74

3.4

Valores singulares dispostos em ordem decrescente - Modelo SPE10-Modificado

75

3.5

Compara¸cão do VPL calculado pelo SIMPAR e pelo TPWL/POD para o modelo SPE10 - Modificado sem compressibilidades e densidades iguais.

3.6

. . . . . . .

Compara¸cão do VPL calculado pelo SIMPAR e pelo TPWL/POD para o modelo SPE10 - Modificado sem compressibilidades e densidades diferentes.

3.7

. . . . . . . . . .

80

Vazã o de óleo dos po¸cos - SPE10 - Modificado compress´ıvel e com densidades de fluidos diferentes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4.1

79

Compara¸cão do VPL calculado pelo SIMPAR e pelo TPWL/POD para o modelo SPE10 - Modificado com compressibilidades e densidades diferentes.

3.8

78

83

Interpola¸caõ por krigagem - regressão polinomial adicionada a deforma¸co˜es locais causadas por uma fun¸cão aleatória com correla¸cã o espacial definida . . . .

92

4.2

Hipercubo latino 2D com 4 amostras . . . . . . . . . . . . . . . . . . . . . . .

93

4.3

Esquema de amostragem e atualiza¸caõ da região de confian¸ca no algoritmo MSAO 96

4.4

Evolu¸caõ do processo de otimiza¸caõ nas itera¸co˜es do algoritmo - Modelo SPE10 - Modificado . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4.5

105

Vazões totais do campo obtidas no processo de otimiza¸caõ — SPE10 - Modificado106

4.6

Pressões de fundo dos po¸cos produtores otimizadas - Modelo SPE10 - Modificado107

4.7

Compara¸caõ do VPL obtido durante o processo de otimiza¸ca˜ o - 1a Itera¸caõ SPE10 - Modificado . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4.8


4.9

108 108


109

4.10

Valores singulares dispostos em ordem decrescente - Modelo Brugge-Modificado 110

4.11

Evolu¸ caõ do processo de otimiza¸caõ nas itera¸co˜es do algoritmo - Modelo SPE10 - Modificado . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

112

4.12

Vazões totais do campo obtidas no processo de otimiza¸caõ — Brugge - Modificado113

5.1

Compara¸cão entre MC Multifidelidade e Tradicional - Reservatório Incompress´ıvel e Densidades dos Fluidos Iguais . . . . . . . . . . . . . . . . . . . . . . . . . .

5.2

Evolu¸caõ dos Parâmetros do MC Multifidelidade para Estima¸caõ da Média e da Variância - Reservatório Incompress´ıvel e Densidades dos Fluidos Iguais . . . .

5.3

124

Compara¸cão entre MC Multifidelidade e Tradicional - Reservatório Compress´ıvel e Densidades dos Fluidos Diferentes . . . . . . . . . . . . . . . . . . . . . . . .

5.4

123

125

Evolu¸caõ dos Parâmetros do MC Multifidelidade para Estima¸caõ da Média e da Variância - Reservatório Compress´ıvel e Densidades dos Fluidos Diferentes . .

126

Lista de Tabelas 2.1

Compara¸cão de volumes in place calculados pelos simuladores . . . . . . . . .

49

3.1

Sequências de controles utilizadas no estudo . . . . . . . . . . . . . . . . . . .

77

3.2

RMSD percentual para os dois casos considerados . . . . . . . . . . . . . . . .

81

4.1

Otimiza¸caõ do VPL utilizando diferentes esquemas - SPE10 com 4 ciclos de controle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4.2

103

Otimiza¸caõ do VPL em diferentes esquemas - Brugge-Modificado com 3 ciclos de controle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

111

Sum´ ario 1.

Introdu¸c˜ ao e Revis˜ ao Bibliogr´ afica

17

1.1

Revisão Bibliogr´ afica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

19

1.2

Ob jetivos do Trabalho . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

22

1.3

Organiza¸cão da Disserta¸caõ . . . . . . . . . . . . . . . . . . . . . . . . . . . .

23

1.4

Lista de publica¸co˜es . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

26

2.

Conceitos de Simula¸ c˜ ao de Reservat´ orios

28

2.1

Equa¸cões de fluxo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

30

2.2

Discretiza¸cão das equa¸cões de fluxo . . . . . . . . . . . . . . . . . . . . . . . .

37

2.3

Cálculo das vazões dos po¸cos . . . . . . . . . . . . . . . . . . . . . . . . . . . .

40

2.4

Cá lculo do Valor Presente L´ıquido . . . . . . . . . . . . . . . . . . . . . . . . .

42

2.5

Compara¸cão com simulador comercial . . . . . . . . . . . . . . . . . . . . . . .

43

3.

Modelo de Ordem Reduzida

50

3.1

Lineariza¸cão utilizando TPWL . . . . . . . . . . . . . . . . . . . . . . . . . .

51

3.2

Exporta¸caõ de mapas de estados e derivadas . . . . . . . . . . . . . . . . . . .

56

3.3

Redu¸caõ de dimensã o - Proper Orthogonal Decomposition (POD) . . . . . . .

59

3.3.1 Descri¸cão do POD . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

60

3.3.2 Decomposi¸caõ em Valores Singulares . . . . . . . . . . . . . . . . . . . . . . .

63

3.4

Modelo de Ordem Reduzida - TPWL/POD . . . . . . . . . . . . . . . . . . .

64

3.5

Detalhes de Implementa¸caõ . . . . . . . . . . . . . . . . . . . . . . . . . . . .

69

3.6

Compara¸cão com o modelo de alta fidelidade . . . . . . . . . . . . . . . . . . .

73

4.

Aplica¸ca õ em Otimiza¸c˜ ao

85

4.1

Algoritmos de Otimiza¸ca˜ o baseados em Modelos Substitutos . . . . . . . . . .

87

4.2

Modelos de Kriging . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

88

4.3

Amostragem por Hipercubo Latino . . . . . . . . . . . . . . . . . . . . . . . .

92

4.4

Algoritmo SAO Multifidelidade . . . . . . . . . . . . . . . . . . . . . . . . . .

93

4.5

Detalhes de Implementa¸caõ . . . . . . . . . . . . . . . . . . . . . . . . . . . .

99

4.6

Desempenho do algoritmo MSAO . . . . . . . . . . . . . . . . . . . . . . . . .

100

4.6.1 Resultados com modelo SPE10 - Modificado . . . . . . . . . . . . . . . . . . .

101

4.6.2 Resultados com modelo de Brugge - Modificado . . . . . . . . . . . . . . . . .

109

5.

Aplica¸ca õ em Propaga¸c˜ ao de Incertezas

114

5.1

M´ etodo das Variáveis de Controle . . . . . . . . . . . . . . . . . . . . . . . . .

115

5.2

Monte Carlo Multifidelidade . . . . . . . . . . . . . . . . . . . . . . . . . . . .

116

5.3

Implementa¸caõ do algoritmo . . . . . . . . . . . . . . . . . . . . . . . . . . . .

119

5.4

Resultados obtidos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

121

6.

Conclus˜ oes e Trabalhos Futuros

127

6.1

Resultados e conclus˜ oes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127

6.2

Sugestões de Trabalhos Futuros . . . . . . . . . . . . . . . . . . . . . . . . . .

Referˆ encias

130

131

˜ E REVISAO ˜ 1. INTRODUC ¸ AO ´ BIBLIOGRAFICA O avan¸co dos computadores, com o aumento da sua capacidade de processamento e a paraleliza¸cão, tem permitido à simula¸caõ de reservatórios, através de suas previs˜ oes, passar a ter papel preponderante nas tomadas de decisão sobre o gerenciamento dos campos. Modelos cada vez mais detalhados, com um número de blocos cada vez maior e capazes de representar as diferentes f´ısicas, têm sido desenvolvidos considerando grandes e complexos sistemas de equa¸cões que consomem enormes recursos computacionais. Por outro lado, o significativo n´ıvel de incertezas nos parâmetros do reservatório (geologia, fluido, etc) e a necessidade de otimiza¸caõ dos recursos utilizados para a produ¸caõ (como equipamentos e investimentos) têm fomentado a concep¸cão de estudos cada vez mais complexos e que requerem grande número de simula¸cões. Estudos como otimiza¸cão dos controles dos po¸c os e análise de incertezas, consideram todo o tempo de concessão dentro da simula¸cão de reservatórios e podem possuir um número grande de variáveis. Sendo assim, o uso de modelos simplificados (baixa fidelidade), que aproximem o modelo de reservatórios (alta fidelidade) com boa acurácia e grande economia de tempo de simula¸caõ, tem se disseminado, permitindo o desenvolvimento de algoritmos multifidelidade para estudos de otimiza¸caõ e incertezas. Neste contexto, o presente trabalho teve como objetivo a utiliza¸cão de um modelo simplificado de ordem reduzida para a diminui¸caõ do custo do processo de simula¸cã o de reservatórios e aplica¸caõ a dois tipos de estudo de gerenciamento de reservatórios de petróleo: otimiza¸caõ dos controles dos po¸cos e propaga¸caõ de incertezas dos controles. Tal modelo de ordem reduzida possui um custo para seu treinamento, mas apresenta grande ganho de tempo e boa acurácia, podendo ser retreinado quando necessário. A técnica escolhida para a constru¸cão do modelo simplificado foi o TPWL (Tra-

17

18 jectory Piecewise Linearization ) combinado ao POD (Proper Orthogonal Decomposition ). A primeira faz basicamente a lineariza¸caõ das equa¸cões de fluxo resolvidas no simulador em torno de estados convergidos e gravados durante uma ou mais simula¸cões de treinamento e é responsável pela redu¸cão da complexidade numérica. A segunda consiste na redu¸cã o da dimens˜ ao do problema através da proje¸caõ dos estados em um espa¸co reduzido gerado pelas dire¸co˜es principais de maior variabilidade. O TPWL elimina a necessidade de resolver as equa¸cões de forma iterativa já que as lineariza, transformando a simula¸caõ em uma sequência de solu¸co˜es de sistemas lineares. Entretanto, essencialmente não há redu¸caõ da dimensão do problema, isto é, a mesma quantidade de equa¸co˜es deve ser resolvida. Desta forma, com o objetivo de reduzir drasticamente o tempo de simula¸caõ do metamodelo quando comparado ao tempo gasto pelo simulador, aplicou-se uma técnica conhecida como POD (Proper Orthogonal Decomposition ) ou PCA (Principal Component Analisys ). Esta técnica baseia-se na existência de correla¸cão entre os mapas de satura¸cões e fra¸cões mássicas gravados, de forma que estes possam ser representados em um espa¸co de dimensão bem menor com pouca perda de informa¸cão. Calcula-se uma base com os mapas gravados que será aplicada a`s equa¸cões do TPWL reduzindo drasticamente sua dimensão. Para a constru¸caõ do modelo de ordem reduzida, al´ em de armazenarem-se os estados, faz-se necessário a grava¸caõ de algumas derivadas que aparecem na expressã o da lineariza¸caõ. Sendo assim, o código do simulador precisou ser adaptado para exportar esta gama de dados (mapas e derivadas) de forma a permitir a constru¸cão do modelo TPWL/POD. Isto é conhecido como metodologia semi-intrusiva, pois requer conhecimento do código do simulador, mas as altera¸co˜es realizadas não atingem os cálculos. A escolha do simulador a ser utilizado foi diretamente influenciada por essa necessidade de acesso ao código. O simulador utilizado tem o nome de SIMPAR e foi desenvolvido na década de 90 pelo CENPES/Petrobras e está em processo de revitaliza¸caõ. Sua formula¸cão é um pouco diferente do usual na indústria. As variáveis primárias que descrevem os estados ao inv´ es ´ um simulador black-oil, mas o de pressão e satura¸co˜es são pressões e fra¸cões mássicas. E problema que se quer atacar é bifásico e portanto utiliza uma op¸cão espec´ıfica para modelar isso. O metamodelo constitu´ıdo pelo TPWL/POD foi constru´ıdo utilizando-se o Matlab, software amplamente utilizado para constru¸caõ de protótipos de algoritmos numéricos. Sendo assim, avaliou-se qual seria a melhor forma de exportar as matrizes de mapas e deri-

1.1. Revis˜ ao Bibliogr´ afica

19

vadas calculadas pela simula¸cão de treinamento de tal maneira que fosse fácil de escrever a partir do código em FORTRAN do SIMPAR e de ler a partir do Matlab. Concluiu-se que a melhor maneira seria utilizar uma biblioteca para manipula¸caõ de grandes massas de dados, denominado HDF5 (Hierarchical Data Format Versão 5). O seu uso permitiu que os tempos de escrita e leitura dos dados não tornassem o treinamento do metamodelo muito caro. No estudo de otimiza¸cão, considerou-se o problema de maximizar o valor presente l´ıquido para todo o tempo de concess˜ ao. O algoritmo utilizado é uma adapt¸ caõ do SAO (Sequential Approximation Optimization ) que considera dois n´ıveis de fidelidade. O metamodelo considerado é o TPWL/POD corrigido por uma superf´ıcie de resposta com krigagem. Esta é constru´ıda no interior da região de confian¸ca através da simula¸caõ de alta (SIMPAR) e baixa (TPWL/POD) fidelidades em alguns pontos e interpola o erro entre os dois modelos. A regi˜ ao de confian¸ca (dimensões e centro) é atualizada de acordo com a qualidade da aproxima¸caõ corrigida. Também foi proposto um critério para o retreinamento do TPWL/POD, adicionando os estados calculados na simula¸caõ do centro da região de confian¸ca ao conjunto de estados inicialmente considerados na constru¸cão do modelo. Esta abordagem baseia-se no fato de que a acurácia do TPWL é garantida nas proximidades da trajetória da simula¸caõ de treinamento e que à medida que a qualidade da aproxima¸caõ diminua, faz-se necessário adicionar uma nova trajetória mais próxima aos u ´ ltimos estados calculados. No estudo de incertezas, considerou-se um algoritmo Monte Carlo multifidelidade que tira proveito da correla¸cão entre os modelos de alta e baixa fidelidade para reduzir o número de execu¸co˜es do simulador de reservatórios. O esquema estudado inclusive calcula qual deve ser a propor¸caõ o´tima entre o número de simula¸co˜es de alta e baixa fidelidade baseado na correla¸caõ entre elas e na razão entre os custos computacionais.

1.1

Revis˜ ao Bibliogr´ afica Na a´rea de gerenciamento de reservatórios dois problemas de grande interesse são

a otimiza¸caõ dos controles dos po¸cos e a propaga¸cão de incertezas. A estratégia de produ¸caõ mais comum dos campos de petróleo no Brasil e no mundo utiliza-se da inje¸ca˜ o de água, sendo a mesma utilizada nos casos aqui estudados. A fun¸caõ objetivo a ser avaliada, tanto nos casos de otimiza¸caõ como para a análise de incertezas, é o Valor Presente L´ıquido, ou simplesmente, VPL, que mede o retorno econômico do campo ao fim do tempo de concessão.


20

Os controles aqui considerados são as pressões de fundo dos po¸cos. Isto, na pr´ atica, com o advento das completa¸cões inteligentes, é uma simplifica¸cão grosseira e deveria ser substitu´ıdo por modelos integrados que simulariam válvulas e facilidades, transferindo os controles do fundo dos po¸cos para os verdadeiros pontos de controle da opera¸caõ. Existe uma extensa literatura sobre otimiza¸cão de controles de po¸cos. Pode-se classificar os métodos de acordo com o grau de intrus˜ a o no código do simulador de reservatórios. Os métodos mais intrusivos utilizam-se da técnica adjunta para calcular o gradiente da fun¸cão objetivo com rela¸cão a`s variáveis de controle (Jansen, 2011) e estão entre os métodos mais eficientes (Brower and Jansen, 2004; Sarma et al., 2008; Chen et al., 2010, 2012). Os métodos adjuntos necessitam de grande esfor¸co de programa¸caõ para serem implementados e portanto não estão dispon´ıveis nos simuladores comerciais até agora. Os métodos semi-intrusivos utilizam-se de modelos de ordem reduzida como Trajectory Piecewise Linearization - TPWL (Cardoso, 2009; He, 2010; Cardoso and Durlofsky, 2010; He et al., 2011; He and Durlofsky, 2013) ou Discrete Empirical Interpolation - DEIM (Gildin et al., 2013; Chaturantabut and Sorensen, Dez/2009) ou ainda podem utilizar a simula¸caõ de linhas de fluxo (Alhuthali et al., 2009). Os métodos não intrusivos usam o simulador como uma “caixa preta” e são puramente baseados nos dados das simula¸co˜es anteriores. Alguns destes algoritmos usam técnicas evolucionárias (Oliveira, 2006; Almeida et al., 2010; Souza et al., 2010), métodos de busca de padrões (Asadollahi and Naevdal, 2010) e métodos baseados em modelos substitutos (Queipo et al., 2002; CMOST, 2012). Existe uma outra classe de algoritmos considerados derivative free (sem derivadas) que aproximam os gradientes da fun¸caõ objetivo utilizando metodologias estocásticas (Wang et al., 2009) e métodos baseados em conjunto (ensemble-based ) (Chen and Oliver, 2010) que podem ser corrigidos por cálculos de diferen¸cas finitas (Xia and Reynolds, 2013) ou serem incorporados em um modelo de interpola¸caõ quadrático (Zhao et al., 2011). Uma discussão mais completa é encontrada em (Conn et al., 2009). No estudo de otimiza¸cão da produ¸cão realizado neste trabalho, o método considerado é do tipo semi-intrusivo conhecido como Trajectory Piecewise Linearization - TPWL (Rewienski, 2003). Este foi inicialmente aplicado ao problema de simula¸caõ de reservatórios em (Cardoso, 2009; Cardoso and Durlofsky, 2010) e aperfei¸coado em (He, 2010; He et al., 2011). Este é um modelo de ordem reduzida f´ısico que se baseia na lineariza¸caõ das equa¸cões de fluxo em torno de estados previamente convergidos e gravados durante uma ou mais ´ utilizado como redutor da complexidade numérica na aprosimula¸cões de treinamento. E xima¸cão (Markovinovic, 2003). Este método é normalmente combinado com a técnica Proper


21

Orthogonal Decomposition - POD para reduzir a dimensão do problema. Existem alguns relatos de instabilidade na literatura acerca deste método (He, 2010; Gildin et al., 2013). Entretanto, o problema de instabilidade parece ter sido resolvido ao trocar-se a proje¸cão de Bubinov-Galerkin por uma proje¸caõ de Petrov-Galerkin (Carlberg et al., 2009; He and Durlofsky, 2013). Os desenvolvimentos do TPWL/POD em (Cardoso, 2009; He, 2010) foram aplicados ao problema de simula¸cão bifásica (óleo e água) com formula¸cão baseada em satura¸co˜es. Já em (He and Durlofsky, 2013) aplicou-se a técnica de redu¸caõ a uma formula¸caõ composicional baseada em fra¸cões molares. No presente trabalho a formula¸caõ considerada é a´gua-óleo baseada em fra¸co˜es mássicas. Os aspectos te´ oricos e as vantagens desta formula¸caõ podem ser obtidas em (Maliska et al., 1997). Como a simula¸cão de reservatórios tem um elevado custo computacional, muitas vezes não é poss´ıvel acoplá-la diretamente ao algoritmo de otimiza¸cão. Sendo assim, utilizase um modelo substituto composto pelo TPWL/POD corrigido por uma aproxima¸cã o por krigagem de dados obtidos em simula¸co˜es no interior de uma região de confian¸ca (Giunta, 2002; Forrester et al., 2008). Efetua-se então a simula¸caõ do modelo de reservatórios (alta fidelidade) e do TPWL/POD (baixa fidelidade) em pontos amostrados no interior da região de confian¸ca utilizando-se a técnica do hipercubo latino. Um preditor de krigagem é constru´ıdo com a diferen¸ca dos dois e ser´ a utilizado na corre¸cão do TPWL/POD durante o processo de otimiza¸caõ no interior desta região de confian¸ca utilizando-se o algoritmo de otimiza¸caõ sequencial quadrática (SQP) (Powel, 1978). O algoritmo baseado em um modelo substituto válido no interior de uma região de confian¸ca, tranformando o problema de otimiza¸caõ em uma sequência de problemas aproximados, é conhecido na literatura como Sequential Aproximation Optimization Algorithm - SAO (Alexandrov et al., 1997). Tal algoritmo, em sua versão clássica considerando somente a krigagem como modelo substituto, foi aplicado ao problema de otimiza¸cão do gerenciamento de campos de petróleo em (Horowitz et al., 2013). No caso deste trabalho, como utilizam-se dois n´ıveis de fidelidade, dá-se o nome de Multifidelity SAO - MSAO. Com rela¸caõ a aplica¸ca˜ o de m´ ultiplas fidelidades à quantifica¸caõ de incertezas, existem alguns esfor¸cos na literatura reportados em (Ng et al., 2013). Em análise de confiabilidade, uma abordagem muito parecida com a otimiza¸cão multifidelidade consiste em executar um número grande de simula¸cões do modelo de baixa fidelidade de forma a aproximar a fronteira de estados limite e corrigir esta estimativa com poucas simula¸cões de alta

1.2. Objetivos do Trabalho

22

fidelidade (Li and Xiu, 2010). Aproxima¸ co˜es polinomiais (e. g. expansão em caos polinomial) das sa´ıdas do modelo de alta fidelidade podem ser constru´ıdas combinando-se uma aproxima¸caõ de baixa fidelidade em um gird mais fino e uma corre¸caõ que utiliza o modelo de alta fidelidade num grid mais grosseiro (Ng and Eldred, 2012). Tamb´ em existem abordagens para amostragem multifidelidade para propaga¸cão de incertezas que se baseiam numa regressão Bayesiana entre as sa´ıdas dos modelos de alta e baixa fidelidade. A aplica¸caõ de quantifica¸caõ de incerteza aqui estudada baseia-se na técnica de vari´ aveis de controle (control variate ) aplicada à simula¸cão Monte Carlo (Glynn and Szechtman, 2000; Nelson, Nov-Dez/1990) e foi desenvolvida para o problema de estima¸cã o de parâmetros estat´ısticos da sa´ıda do modelo de alta fidelidade como, por exemplo, média e variˆ ancia dadas as distribui¸co˜es das variáveis de entrada.

1.2

Objetivos do Trabalho O desenvolvimento de modelos multifidelidade em geral e aplicados a simula¸caõ

de reservatórios tem crescido bastante nos últimos anos. Estudos de gerenciamento de reservatórios como otimiza¸cão de controles ou posi¸caõ dos po¸cos, análise de incertezas ou a combina¸ca˜ o dos dois u ´ ltimos, resultando na otimiza¸cão sob incertezas necessitam de um número cada vez maior de simula¸cões. Por outro lado a tendência é que estes modelos de simula¸cão se tornem cada vez mais complexos e portanto, mesmo com o desenvolvimento de computadores mais rápidos e com a tendência de paraleliza¸caõ, consumam ainda muito tempo. Sendo assim, se torna imperativo o uso de modelos substitutos simplificados que aproximem o resultado das simula¸co˜es mas consumam muito menos tempo. Esta disserta¸cão não está dedicada ao desenvolvimento de uma nova técnica para concep¸caõ de modelos substitutos, mas sim à utiliza¸caõ de um tipo de modelo de ordem reduzida em estudos de otimiza¸caõ e de propaga¸cão de incertezas utilizando o recurso da simula¸cão de reservatórios. A ideia é construir um modelo baseado na combina¸caõ de duas técnicas: Trajectory Piecewise Linearization - TPWL e Proper Orthogonal Decomposition POD como redutores de complexidade numérica e dimensão, respectivamente. Este modelo precisa se acoplar a um simulador de reservat´ orios e, por necessitar de mais dados do que o normalemente exportado, necessita de conhecimento do código-fonte do simulador. Finalmente, dois algoritmos, um para otimiza¸caõ e outro para propaga¸cão de incertezas, são constru´ıdos e adaptados para utilizarem um esquema multifidelidade que considere o modelo

1.3. Organiza¸ cão da Disserta¸caõ

23

de ordem reduzida em questão. São considerados objetivos espec´ıficos deste trabalho:

• Estudar em detalhe a formula¸cão de fluxo e o código-fonte do simulador SIMPAR, desenvolvido pelo CENPES/Petrobras para poder utilizá-lo em algoritmos semi-intrusivos

de otimiza¸cão. Contribuir com a proposta de revitaliza¸cão do simulador para uso nas pesquisas desempenhadas nas universidades brasileiras inseridas na rede SIGER.

• Estudar em detalhes a técnica e desenvolver um código para o modelo de ordem reduzida baseado no TPWL/POD que possa ser acoplado ao simulador SIMPAR.

• Avaliar a proje¸caõ de Petrov-Galerkin em substitui¸caõ de Bubnov-Galerkin na estabiliza¸caõ do TPWL.

• Estudar a metodologia e desenvolver um código para o algoritmo de otimiza¸caõ multifidelidade baseado em modelos substitutos que se acople ao simulador em questão e ao modelo de ordem reduzida estudado considerando uma única fun¸cão objetivo.

• Propor uma metodologia para retreinamento do TPWL ao longo do processo de oti-

miza¸caõ, com a defini¸cão de um critério para decidir quando o retreinamento deve ocorrer.

• Estudar o arcabou¸co matemático e estat´ıstico e desenvolver um código para o algoritmo Monte Carlo Multifidelidade acoplado ao mesmo simulador e ao mesmo modelo de ordem reduzida e considerando incertezas nos controles dos po¸cos.

1.3

Organiza¸ c˜ ao da Disserta¸c˜ ao A disserta¸caõ foi dividida em quatro partes que correspondem exatamente aos

passos transcorridos na pesquisa. O Cap´ıtulo 2 contem todo o conte´ udo de simula¸cã o de reservatórios necessário para o desenvolvimento do restante da pesquisa. O software SIMPAR, desenvolvido na década de 90 pelo CENPES/Petrobras, foi escolhido como simulador de reservatórios a ser utilizado neste trabalho. A escolha se deve à necessidade de acesso irrestrito ao código do simulador. No entanto, apesar de ser considerado um simulador black-oil , possui uma formula¸cão distinta da maioria dos softwares comerciais. Suas vari´ aveis primárias são pressão


24

e fra¸co˜es mássicas, diferentemente da formula¸cão por pressões e satura¸co˜es usualmente utilizada nos simuladores comerciais. Esta formula¸caõ foi escolhida no desenvolvimento desse simulador por ser mais facilmente extens´ıvel ao caso composicional e por eliminar o principal inconveniente da formula¸caõ por balan¸co de fases: o desaparecimento da fase gasosa em caso de subsatura¸caõ. Sendo assim, um cap´ıtulo foi dedicado a` dedu¸caõ da equa¸caõ de fluxo e à sua discretiza¸caõ e a outros aspectos importantes, como modelo de po¸co e as defini¸cões de grandezas básicas na simula¸caõ, sempre atendo-se às peculiaridades do simulador que será utilizado. Isto é uma tentativa de fazer da disserta¸caõ um documento autocontido, permitindo a um leitor leigo obter os fundamentos básicos da simula¸cão de reservatórios que terá papel preponderante nos próximos cap´ıtulos. Este mesmo cap´ıtulo possui alguns resultados de valida¸cão do SIMPAR, comparando resultados obtidos com ele e com um simulado comercial, o IMEX da CMG. Tais resultados foram gerados com dois modelos de reservatórios. Um deles é um corte do modelo SPE10 e possui 24.000 células e o outro é uma adapta¸caõ com malha cartesiana para o modelo de Brugge amplamente divulgado na literatura e contendo em torno de 60 .000 células. Ambos consideram somente o fluxo de o´leo e água. O Cap´ıtulo 3 dedica-se a estudar o modelo de ordem reduzida escolhido para ser aplicado. Este modelo utiliza a técnica TPWL (Trajectory Piecewise Linearization ) como redutor da complexidade numérica e o POD (Proper Orthogonal Decomposition ) como redutor da dimensão do problema. Os aspectos teóricos sobre ambas as técnicas foram expostos, com toda a formula¸caõ aplicada ao problema de simula¸caõ de fluxo em reservatórios. Entretanto, este modelo reduzido não é meramente baseado em resultados de simula¸cõ es. O TPWL concebe um modelo substituto f´ısico através da lineariza¸caõ da equa¸cão de fluxo em torno de estados convergidos em uma simula¸caõ de treinamento, transformando a simula¸caõ de reservatórios em uma sequência de solu¸co˜es de sistemas lineares. Portanto, depende da exporta¸cão não só dos mapas de estados mas também de diversas derivadas dos res´ıduos. Os detalhes da implementa¸caõ, como, por exemplo, a maneira de exporta¸cão dos dados do simulador e as ferramentas escolhidas para a constru¸caõ do modelo estão descritos. Também é descrito o processo de retreinamento em torno de uma nova simula¸caõ e que virá a ser utilizado no algoritmo de otimiza¸cão. Alguns resultados comparativos da simula¸ cão do SIMPAR (alta fidelidade) com o TPWL (baixa fidelidade) são apresentados com uma u ´ nica simula¸cão de treinamento e considerando retreinamento em torno de diferentes simula¸cões de modo a avaliar a acurácia da aproxima¸caõ de ordem reduzida. Para tal estudo é utilizado o


25

modelo SPE10 apresentado no cap´ıtulo 2 considerando três versões: uma totalmente incompress´ıvel e com densidades de ambos os fluidos idênticas (caso mais estudado na literatura até então), uma considerando as densidades diferentes, mas ainda total incompressibilidade (o que aumenta enormemente a não linearidade do problema) e uma tamb´ em com densidades diferentes e considerando compressibilidade de rocha e para os fluidos. Prova-se que a acurácia do modelo reduzido diminui consideravelmente com a considera¸cão das densidades diferentes e se reduz um pouco mais ao considerar-se as compressibilidades. O Cap´ıtulo 4 ataca o problema da otimiza¸caõ dos controles dos po¸c os em um campo de petróleo. A fun¸cão objetivo a ser maximizada é o Valor Presente L´ıquido - VPL e as variáveis de projeto são os controles dos po¸cos, que neste trabalho são as pressões de fundo de po¸c o. Sendo assim, propõe-se uma estratégia multifidelidade na qual o modelo de ordem reduzida TPWL é válido no interior de uma região de confian¸ca e corrigido por um modelo de krigagem. Para tanto, faz-se necess´ ario treinar o TPWL considerando uma simula¸cão controlada por uma sequência aleatória de pressões de fundo e a cada itera¸cão devese amostrar o espa¸co contido no interior da região de confian¸ca e simular os pontos amostrados com o modelo de alta fidelidade (simulador SIMPAR) e baixa fidelidade (TPWL). Controi-se então uma superf´ıcie de resposta para o erro entre eles que deve ser adicionada ao TPWL como corre¸caõ. Desta forma, a cada itera¸cão o centro e as dimensões da região de confian¸ca são atualizados de acordo com a qualidade da aproxima¸caõ. A otimiza¸caõ realizada no interior desta região baseia-se no algoritmo Sequential Quadratic Programming e é aplicada ao modelo aproximado com corre¸caõ. Isto é o que se conhece como algoritmo de otimiza¸cão aproximada sequencial multifidelidade (MSAO - Multifidelity Sequential Approximation Optmization ), o qual transforma a otimiza¸cão do problema com a fun¸caõ objetivo real em uma sequência de otimiza¸co˜es aproximadas. Eventualmente espera-se que haja a necessidade de retreinamento do TPWL. Isto consiste em adicionar ao conjunto de estados já considerados no primeiro treinamento, um conjunto de estados e derivadas calculados na simula¸caõ do centro da última região de confian¸ca, o que aumenta a acurácia do metamodelo nas proximidades deste ponto. Alguns resultados foram gerados com os modelos SPE10 e Brugge modificados, em ambos os casos considerando densidades diferentes para o´leo e água e compressibilidades não-nulas. Várias situa¸co˜es foram testadas variando-se o crit´ erio de retreinamento, podendose assim avaliar o desempenho desta técnica. Como compara¸caõ o otimizador foi conectado diretamente ao simulador e ao TPWL para avaliar os ganhos de utilizar-se o MSAO.

1.4. Lista de publica¸co˜es

26

A u ´ltima parte da pesquisa, relatada no cap´ıtulo 5, consisitu na utiliza¸ca˜ o do TPWL em um algoritmo para propaga¸caõ de incertezas. Este é um tipo de Monte Carlo Multifidelidade, baseado na metodologia das variáveis de controle (Control Variates Method ) e se vale da correla¸cão entre os modelos de alta e baixa fidelidade para reduzir o número de simula¸co˜es utilizando o simulador, substituindo-as de maneira ótima por simula¸co˜es de baixa fidelidade. O arcabou¸co estat´ıstico e matemático bem como a dedu¸caõ do método são apresentados, construindo-se um algoritmo que gradativamente aumenta o número de avalia¸cões da fun¸cão VPL de alta fidelidade e é capaz de calcular a média e o desvio padrão desta variável com menos simula¸co˜es do que o método de Monte Carlo convencional. Os resultados obtidos utilizam o modelo baseado no SPE10 na versão incompress´ıvel e com densidades iguais e compress´ıvel com densidades iguais. Ao final do trabalho algumas conclusões acerca do uso do TPWL em estudos de otimiza¸caõ e de propaga¸cão de incertezas são apresentadas no cap´ıtulo 6. Algumas sugest˜ oes de trabalhos futuros são apresentadas, considerando tanto novos estudos que precisariam ser realizados como melhorias dos algoritmos aqui utilizados.

1.4

Lista de publica¸ co ˜es Esta disserta¸caõ culminou na publica¸caõ dos quatro resumos e trabalhos completos

em anais de congressos internacionais listados abaixo.

• Fragoso, Manuel Jr.; Horowitz, Bernardo; Rodrigues, José Roberto Pereira. A Re-

duced Order Surrogate Model for Optimal Reservoir Management. 14th European Conference on Mathematics of Oil Recovery, ECMOR. Set, 2014. Catania, Sicilia, It´ alia.

• Fragoso, Manuel Jr.; Horowitz, Bernardo; Rodrigues, José Roberto Pereira. A Multi-

fidelity Approach to Waterflooding Optimization. 4th International Conference on Engineering Optimization, EngOpt . Set, 2014. Lisboa, Portugal.

• Fragoso, Manuel Jr.; Horowitz, Bernardo; Rodrigues, José Roberto Pereira. A Redu-

ced Order Method to Waterflooding Optimization. 35th Iberian Latin American Congress on Computational Methods in Engineering, CILAMCE . Nov, 2014 Fortaleza, Brasil.

• Fragoso, Manuel Jr.; Horowitz, Bernardo; Rodrigues, José Roberto Pereira. Retrai-

1.4. Lista de publica¸co˜es

27

ning Criteria for TPWL/POD Surrogate Based Waterflodding Optimization. Reservoir Simulation Simposium, RSS . Fev, 2015. Houston, EUA. Também resultou na apresenta¸caõ oral do seguinte trabalho.

• Horowitz, Bernardo; Afonso, Silvana Maria Bastos; Fragoso, Manuel Jr. Surrogate Based Waterflooding Optimization. 1st International Symposium on Energy Chalenges and Mechanics, ECM . Jul, 2014. Aberdeen, Escócia, Reino Unido. Finalmente, culminou na aceita¸caõ para a publica¸caõ do seguinte artigo de revista que se encontra em processo de revisão.

• Fragoso, Manuel Jr.; Horowitz, Bernardo; Rodrigues, José Roberto Pereira. A Re-

duced Order Surrogate Model for Optimal Reservoir Management . Selected contributions from the 14th European Conference on the Mathematics of Oil Recovery (ECMOR XIV). Computational Geosciences, Springer.

˜ 2. CONCEITOS DE SIMULAC ¸ AO ´ DE RESERVATORIOS Com o aumento da demanda global por hidrocarbonetos e a consequente alta no pre¸co do petróleo, novas fronteiras de explora¸cão de petróleo cada vez mais complexas têm se mostrado viáveis. Com isso, engenheiros de reservatórios se deparam com desafios crescentes para gerar previsões de produ¸cão e estimativas de reserva. Para realizar tais previsões de produ¸caõ, são constru´ıdos modelos matemáticos que representam o que ocorre nas rochas que servem de acumula¸cão para os hidrocarbonetos em subsuperf´ıcie. O problema básico a ser resolvido é o cálculo do fluxo através do reservatório e das vazões de produ¸cão e inje¸cão em cada po¸co, dados os controles aplicados ao longo do tempo. Note que, em geral, ocorre a produ¸caõ de três fluidos distintos: a´gua, o´leo e gás. A composi¸cão das fases óleo e gás pode ser bastante complexa, o que pode levar a considera¸caõ de diferentes componentes presentes nestas fases e à concep¸cão de modelos denominados composicionais. Esta considera¸cão aumenta consideravelmente a complexidade do problema e não será considerada neste trabalho. O modelo aqui considerado se limitar´ a ao fluxo de três fluidos (água, o´leo e gás) que fluem independentemente de sua composi¸caõ. ´ o que se dá o nome de modelo black-oil . E Por outro lado, há muitas incertezas nos parâmetros dos modelos de reservatório que são constru´ıdos de modo que nenhum estudo de reservatórios pode negligencia-las. Isto complica ainda mais o processo, demandando o uso de técnicas estat´ısticas espec´ıficas para lidar com incertezas, além de consumir muito recurso computacional com a necessidade de um n´ umero grande de simula¸co˜es. Todavia, o interesse final do processo de previsão de produ¸cão é a defini¸caõ da estratégia ótima para a produ¸caõ de uma acumula¸caõ de petróleo, considerando as muitas maneiras de posicionar ou controlar os po¸cos. Para tanto, demanda-se a utiliza¸caõ de técnicas 28

29 de otimiza¸caõ e aumenta-se ainda mais o n´ umero de simula¸co˜es a serem executadas. Todo este processo significa em suma estimar a produ¸caõ em m´ ultiplos cenários de forma a quantificar riscos e melhorar os resultados. Tudo isto, além da necessidade de modelos cada vez mais detalhados, motiva o desenvolvimento de simuladores de reservatório mais eficientes e acurados. Nesse cap´ıtulo serão abordados os aspectos do modelo matemático para fluxo em meios porosos utilizado na simula¸caõ de reservatórios de petróleo e as principais caracter´ısticas do software (simulador de fluxo) utilizado em todos os estudos deste trabalho. Utilizou-se um simulador desenvolvido por uma equipe do CENPES/Petrobras no in´ıcio da ´ um simulador década de noventa e conhecido como SIMPAR (Rodrigues and Bonet, 1991). E do tipo black-oil (Ertekin et al., 2011), isto é, que considera o fluxo de três fases (aquosa, oleosa e gasosa) e três pseudo-compononentes (óleo, á gua e gás). Na verdade, este simulador é capaz de considerar um quarto componente chamado de pol´ımero o que possibilita a simula¸caõ de estratégias de produ¸caõ que considerem a inje¸caõ de pol´ımeros voltada, por exemplo, para o controle da produ¸caõ de água. Como o objetivo do presente trabalho reside no problema de otimiza¸cão da inje¸caõ de água em um reservatório subsaturado, utlizar-se-á uma simplifica¸caõ do modelo blackoil , conhecida como modelo água-óleo, que considera somente a existência de duas fases no reservatório (á gua e óleo). Quando configurado para simular este modelo, o SIMPAR considera que o componente gás existe somente em solu¸cã o na fase óleo e, portanto, caso o usuário determine uma razão de solubilidade, haverá produ¸caõ de gás, mas esse gás não estará fluindo livremente no reservatório. No presente estudo, considera-se a razão de solubilidade igual a zero somente por simplicidade. Uma importante caracter´ıstica do SIMPAR que o diferencia da maioria dos simuladores black-oil comerciais é a sua formula¸caõ por fra¸co˜es mássicas. Ao invés de utilizar-se das equa¸cões diferenciais expressas com base em satura¸cões volumétricas das fases como é comum entre os outros softwares de simula¸caõ blackoil , este simulador utiliza equa¸cões baseadas em fra¸cões mássicas dos componentes. Esta formula¸cão possui vantagens e caracter´ısticas próprias como descrito na se¸caõ 2.1. O cap´ıtulo tratará inicialmente da formula¸caõ matemática do problema, deduzindo muito rapidamente na primeira se¸caõ as equa¸cões diferencias de fluxo para o modelo água-oleo na sua formula¸cão por fra¸cões mássicas. Na segunda se¸caõ, descreve-se a discretiza¸caõ em diferen¸cas finitas do sistema de equa¸co˜es diferenciais, a qual é utilizadoa para

2.1. Equa¸cões de fluxo

30

a solu¸caõ num´ erica realizada no simulador. A terceira se¸cão estuda o tratamento dado aos termos de produ¸cão/inje¸caõ (termos de po¸co) pelo SIMPAR. Isto é devido à importância que o modelo de po¸co tem no processo de cálculo das vazões produzidas e injetadas, as quais são as entradas da fun¸caõ objetivo utilizada ao longo deste trabalho. A quarta se¸cão contem os resultados da compara¸cão dos resultados obtidos com o SIMPAR e com um simulador comer´ uma tentativa de valida¸cão do cial, desenvolvido pela Computer Modeling Group, o IMEX. E primeiro considerando os resultados do outro, amplamente utilizado.

2.1

Equa¸c˜ oes de fluxo Do ponto de vista matemático, a fun¸cão básica de um simulador de reservatórios

é resolver numericamente um sistema de equa¸co˜es diferenciais que modela o fluxo em um meio poroso de fases distintas contendo um ou mais componentes. Em um modelo conhecido como black-oil , são três as fases — oleosa, gasosa e aquosa — e três os componentes — óleo, á gua e gás. Além disso, admite-se que o gás pode estar livre ou dissolvido no óleo. Neste trabalho, só serão considerados dois destes componentes e fases: óleo e a fase oleosa e água e a fase aquosa. Isto é que se chama de modelo água-óleo. Classicamente, a maioria dos simuladores considera em sua formula¸cão o balan¸co de massa entre as fases e, para representar a parcela de cada fase, consideram as satura¸cões volumétricas de cada fase em rela¸cão ao total (Ertekin et al., 2011). Entretanto, o simulador utilizado neste estudo possui uma formula¸caõ distinta, considerando o balan¸co de massa dos componentes. Para tanto, contabiliza-se a participa¸cão de cada um dos componentes através de sua fra¸cão mássica em rela¸cão à massa total. Desta forma, utiliza como variáveis primárias as fra¸co˜es mássicas ao invés das satura¸co˜es (Maliska et al., 1997). Tal formula¸caõ possui algumas vantagens como descritas em (Maliska et al., 1997) e (Rodrigues and Bonet, 1991). Basicamente a escolha se deve a dois fatores: 1. não é necessário nenhum tipo de tratamento especial para o gás. Quando se consideram as satura¸cõ es das fases, se todo o gás estiver dissolvido, a fase g´ a s deixa de estar presente, isto é, sua satura¸caõ é zero, o que exige um tratamento especial para o desaparecimento da fase, por exemplo, mudan¸ca de variáveis ou manuten¸caõ de uma satura¸caõ residual; 2. esta é a formula¸caõ natural para a constru¸caõ de um modelo composicional, o que


31

facilitaria bastante se houvesse o interesse de extender o simulador para simula¸cão composicional. Neste trabalho, como j´ a mencionado, só será considerado o problema bifásico água-óleo. Entretanto, o SIMPAR é capaz de simular um modelo black-oil completo com a presen¸ca de pol´ımeros que seriam injetados. O tipo de modelo a ser simulado é uma das entradas de seu arquivo de entrada conforme seu manual (CENPES, 1995). Abaixo são deduzidas as equa¸co˜es de fluxo para esse modelo de duas fases e dois componentes. Entretanto, a dedu¸cão é extens´ıvel para o problema trifásico, adicionandose uma equa¸caõ de balan¸co e considerando-se as caracte´ısticas espec´ıficas do componente gás, como compressibilidade e solubilidade no óleo. Inicialmente considera-se um modelo unidimensional e um fluido, embora a extensão para três dimensões e dois fluidos seja simples e será considerada ao final desta se¸caõ. Sendo assim, para in´ıcio da dedu¸caõ, considere um volume de controle e seus vizinhos como mostra a Figura 2.1.

Fig. 2.1: Volume de controle e seus vizinhos Fazendo-se o balan¸co de massa no interior desse volume de controle obtem-se a Equa¸caõ 2.1, onde M in é a massa que entra no volume de controle, M out é a massa que sai do volume de controle, M s é a massa injetada/produzida a partir de um po¸co completado no volume de controle e M acc é a massa acumulada no volume de controle.

M in

− M

out + M s = M acc

(2.1)

A vazão mássica m que ˙ entra ou sai do volume de controle é calculada pelo produto da velocidade do fluido v x por sua densidade ρ e pela área da interface A x , como na Equa¸caõ 2.2. A massa que entra ou sai é a integral de m. ˙

· ·

˙ m = A x ρ vx

(2.2)


32

A massa acumulada é obtida na Equa¸cão 2.3. Ela é calculada pela diferen¸ca entre a massa contida no interior do volume de controle no instante t + δt e contida no instante anterior t e depende da porosidade no volume de controle φ e do volume de controle V b .

·

M acc = (ρ φ V b )t+∆t

− (ρ · φ V )

b t

(2.3)

onde Aplicando as Equa¸cões 2.2 e 2.3 à Equa¸cão 2.1, obtem-se a Equa¸cão 2.4.

· ·

(Ax ρ vx,i

1/2

−

−A ·ρ·v

x,i+1/2 +

x

·

· ·

˙ s ) ∆t = (ρ φ V b )t+∆t m

− (ρ · φ · V )

b t

(2.4)

·

Manipulando a equa¸caõ anterior, dividindo-se ambos os membros por V b ∆t, obtem-se a Equa¸caõ 2.5.

·

ρ vx,i

1/2

−

−ρ·v

x,i+1/2

∆x

·

− ·

˙s (ρ φ)t+∆t (ρ φ)t m = ∆t V b

+

(2.5)

Supondo-se que o volume de controle e o intervalo de tempo considerados sejam pequenos o suficiente para uma aproxima¸caõ cont´ınua, a Equa¸cão 2.5 pode ser reescrita em termos da equa¸caõ diferencial 2.6

·

·

˙s ∂ (ρ vx ) m ∂ (ρ φ) + = ∂x V b ∂t

(2.6)

Considerando-se que a rocha reservatório é um meio poroso e que as velocidades dos fluidos em seu interior são pequenas, pode-se calcular tal velocidade vx como fun¸caõ do potencial hidráulico Φ através da lei de Darcy generalizada, conforme a Equa¸caõ 2.7 onde k é a permeabilidade ao fluido (caso haja mais de um fluido, é o produto da permeabilidade absoluta pela relativa) e µ é a viscosidade do fluido.

vx =

k ∂ Φ µ ∂x

·

(2.7)

Aqui define-se o potencial hidráulico conforme a Equa¸caõ 2.8 em fun¸cão da profundidade D em rela¸cão a` uma referência, à pressão do fluido. Nesta equa¸caõ, γ representa o peso espec´ıfico do fluido.

Φ = p

− γD

(2.8)


33

Substituido a expressão para vx da Equa¸caõ 2.7 na equa¸cão de balan¸co 2.6 obtemse 2.9. ∂ ∂x

� · · � ρ

k ∂ Φ µ ∂x

˙s m ∂ = (ρ φ) V b ∂t

·

+

(2.9)

Esta equa¸cão pode ser extendida para três dimensões somando-se as contribui¸co˜es das três dire¸cões conforme Equa¸caõ 2.10.

∂ ∂x

�·

k x ∂ Φ ρ µ ∂x

·

� �· ∂ + ∂y

k y ∂ Φ ρ µ ∂y

·

� �·

∂ k z ∂ Φ + ρ ∂z µ ∂z

·

�

+

˙s m ∂ = (ρ φ) V b ∂t

·

(2.10)

Pode-se simplificar esta equa¸caõ utilizando-se a nota¸caõ vetorial. Para tanto, de forma a contemplar os casos anisotrópicos (permeabilidades diferentes nas três dire¸cões), é necess´ ario definir o tensor de permeabilidades K, matriz cuja diagonal é preenchida com as permeabilidades nas dire¸co˜es x, y e z , como mostrada na Equa¸cão 2.11.

K =

 

kx

0

0

0

ky

0

0

0

kz

 

(2.11)

A Equa¸cão 2.12 é a forma simplificada pela nota¸ cão vetorial da Equa¸cão 2.10. Também considera-se que m ˙ s já representa a vazão injetada por unidade de volume e introduzse uma nova grandeza denominda mobilidade do fluido, cujo tensor é Λ = ρ K/µ. Note que au ´ nica variável primária que representa a fun¸cão a ser calculada nesta equa¸caõ é a pressão que está contida no potencial hidráulico representado na Equa¸cão 2.8.

∇ · (Λ∇Φ) + m˙ = ∂t∂ (ρ · φ) s

(2.12)

Entretanto, existe mais de um fluido percolando o meio poroso. Considerando que cada fluido possui sua própria mobilidade, o modelo deve possuir então duas equa¸co˜es de balan¸co que interagem entre si de modo a representar o escoamento bifásico. Portanto, duas serão as variáveis primárias, em rela¸cão a`s quais o problema será resolvido. Na formula¸caõ aqui considerada, as variáveis primárias escolhidas são a pressã o do o´leo e a fra¸cão mássica de água. As equa¸cões diferenciais que são resolvidas pelo simulador estão expressas nas eqs.


34

2.13 e 2.14 e representam, repectivamente, o balan¸co de massa total e do componente água. Aqui os ´ındices w e o representam respectivamente água e óleo. ∂ (φ ρm ) ∂t

·

− ∇ · (Λ ∇Φ

∂ (φ z w ρm ) ∂t

· ·

Nestas equa¸cões

w +

w

∇

Λo Φo ) + m˙w + m˙ o = 0

− ∇ · (Λ ∇Φ

w) +

w

(2.13)

m˙w = 0

(2.14)

∇Φ representa o gradiente do potencial de cada fase e tem origem

na lei de Darcy, onde a velocidade é proporcional ao gradiente do potencial. Desta forma, ambos os gradientes podem ser expressos nas equa¸cões 2.15 para a fase água e 2.16, para o óleo, respectivamente. Note que, no caso do fluxo de duas fases concorrentes, o gradiente, antes representado na Equa¸cão 2.8, deve contemplar uma parcela que representa a pressão capilar entre as fases que é denotada por P cow .

∇Φ

w

∇ − ∇P − γ ∇D

= po

cow

∇Φ = ∇ p − γ ∇D o

o

w

o

(2.15)

(2.16)

De maneira sint´ etica, podem-se definir as grandezas que aparecem nas Equa¸cões 2.13 a 2.16 na lista a seguir, onde o subscrito l, representa as fases oleosa e aquosa e portanto assume os valores l = o, w e o subscrito c representa os componentes óleo e água e também assume os valores l = o, w.

• z - fra¸cão mássica global do componente c; • ρ - massa espec´ıfica da fase l; • ρ - massa espec´ıfica da mistura, ρ = S ρ • S - satura¸caõ da fase l; • λ =- mobilidade mássica da fase l; • Φ - Potencial da fase l; • p - pressão da fase óleo; • P - pressão capilar no sistema água-óleo; c l

m l

l

l

o

cow

m

w w + S o ρo ;


35

• γ - peso espec´ıfico da fase l, γ = ρ g; • g - acelera¸caõ da gravidade; • D - profundidade tomada em rela¸caõ a uma referência; • m˙ = vazão mássica injetada/produzida da fase l por unidade de volume. l

l

l

l

A mobilidade da fase l é expressa na Equa¸cão 2.17.

λl = k

ρl krl µl

(2.17)

Vale ressaltar que as satura¸cões das fases S l se relacionam com as fra¸co˜es mássicas dos componentes z c através da Equa¸cão 2.18.

S l =

z l /ρl , l = o, w. z w /ρw + z o /ρo

(2.18)

As satura¸co˜es são utilizadas na determina¸cão da permeabilidade relativa k rl , uma vez que suas curvas são expressas em fun¸cão destas variáveis. Tais curvas de permeabilidade relativa podem ser fornecidas através de tabelas contendo alguns pontos que são interpolados ou podem utilizar o modelo de Corey dado pelas Equa¸co˜es 2.19 e 2.20 que representam, respectivamente, as permeabilidades relativas do óleo k ro e da água k rw em fun¸cão da satura¸caõ de água S w . 0 kro = k ro

0 krw = k rw

�− − �− 1 1

1 1

S w S wc

− S − S

or or

− S − − S S w S wc

wc

� �

a

(2.19) b

(2.20)

or

As curvas de permeabilidade relativa no modelo de Corey são calibradas através dos seguintes parâmetros determinados a partir de ensaios:

• S

Satura¸cão de óleo residual que representa a m´ınima satura¸caõ de óleo que se pode

• S

Satura¸caõ de água conata que representa a m´ınima satura¸ca˜ o de água que nor-

or -

obter; wc -

malmente existia no reservatório antes da inje¸cão de água; 0 ro -

•k

ponto terminal para permeabilidade relativa do óleo, medido na satura¸caõ S wc ;


36

0 rw

• k - ponto terminal para permeabilidade relativa da água, medido na satura¸cão 1 − S ; • a - expoente de Corey para a curva de permeabilidade relativa de óleo e • b - expoente de Corey para a curva de permeabilidade relativa de água. or

Para cálculo das mobilidades das fases, λl faz-se necessário determinar suas densidades e viscosidades em fun¸cão da pressão. As densidades podem ser calculada utizando-se os fatores volume de forma¸caõ Bl , l = o, w fornecidos dentre as informa¸co˜es PVT conforme a Equa¸caõ 2.21. 1 ρl,st , Bl onde ρ l,st é a densidade da fase em condi¸caõ padrão. ρl =

(2.21)

As viscosidades µl são calculadas através da Equa¸cão 2.22, onde µl,b é a viscosidade tomada na pressão de bolha,

∂µ é ∂p

a viscosibilidade, p o é a pressão do óleo e p b é a pressão de

bolha.

µl = µ l,b +

∂ µ ( po ∂p

− p ), b

(2.22)

A porosidade φ nas Equa¸cões 2.13 e 2.14 é calculada considerando-se o modelo pouco compress´ıvel com compressibilidade da rocha constante. Desta considera¸cão, pode-se calcular a porosidade através da Equa¸cão 2.23, onde φ l,ref é a porosidade tomada na pressão de referência, c r é a compressibilidade da rocha, po é a pressã o do o´leo e pref é a pressão de referência para o cálculo da porosidade.

φ = φ l,ref + cr ( po

− p

ref ),

(2.23)

Por fim, vale ressaltar que para simplificar o problema, a razão de solubilidade Rs do gá s no óleo foi desconsiderada em todos os estudos. Entretanto, o simulador permite

̸

que se considere Rs = 0 no modelo água-óleo. Neste caso, não há gás livre, mas existe componente gás fluindo junto com o óleo na fase oleosa. Sendo assim, a quantidade de gás produzida seria completamente proporcional ao óleo produzido e seria necessário conderar-se a fra¸caõ mássica do componente gás na fase oleosa X go e do componente óleo, X oo . Isto adiciona alguma complexidade as Equa¸co˜es 2.18 e 2.21 e pode ser encontrado em (Rodrigues and Bonet, 1991).

2.2. Discretiza¸caõ das equa¸cões de fluxo

2.2

37

Discretiza¸c˜ ao das equa¸ c˜ oes de fluxo Para a solu¸cão numérica das Equa¸co˜es 2.13 e 2.14 faz-se necessário discretizá-las.

Existem diversos métodos para discretiza¸caõ: elementos finitos, volumes finitos, métodos multiponto, etc. Entretanto, classicamente, no caso do problema de fluxo, utiliza-se o método das diferen¸cas finitas, aproximando as derivadas temporais e espaciais por diferen¸cas divididas (Ertekin et al., 2011). Este é o método utilizado na formula¸caõ do simulador aqui considerado. A discretiza¸caõ tamb´ em deve considerar uma malha que particiona o espa¸co em pequenos blocos para aproxima¸caõ dos gradientes. Quando se utiliza o método das diferen¸cas finitas, é mais comum utilizarem-se malhas estruturadas (Ertekin et al., 2011), nas quais é poss´ıvel numerar os nós considerando a rela¸caõ com os vizinhos. Apesar de existirem hoje em dia tipos mais flex´ıveis, aqui só serão consideradas malhas regulares com topo variável, isto é, todos os blocos são paralelep´ıpedos retângulos, uma vez que o simulador SIMPAR só é capaz de lidar com este tipo de configura¸caõ. Para discretiza¸cão do espa¸co, considera-se o método por diferen¸cas finitas centradas na interface, também conhecido por CVFD - control volume finite diference ou diferen¸cas finitas em volume de controle (Ertekin et al., 2011). Este método é uma pequena varia¸caõ do método das diferen¸cas finitas centradas, uma vez que ao discretizar-se a lei de ( Darcy), as propriedades são consideradas nas interfaces de cada bloco. A Equa¸caõ 2.24 ilustra a aproxima¸caõ por diferen¸cas finitas em uma única dimensão. A extensão de tal aproxima¸caõ para as três dimensões é imediata. Nesta equa¸cão só estão ilustrados os termos de fluxo e os coeficientes T l representam as transmissibilidades de cada fase na interface entre os blocos e mais adiante serão relacionados a` mobilidade.

� �

∂ ∂ Φl λl ∂x ∂x

n+1

≈

�

1 n+1 (T l )n+1 1 (Φi 1 i 2 V −

−

−

Φn+1 ) + (T l )n+1 (Φn+1 i i+1 i+ 1 2

−

�

Φn+1 ) i

(2.24)

Já no caso da discretiza¸caõ do tempo, a aproxima¸caõ por diferen¸cas finitas corresponde à Equa¸caõ 2.25. Nesta equa¸cão só está ilustrado o termo dependente do tempo, também conhecido como termo de acumula¸caõ. ∂ (φ z l ρm ) ∂t

· ·

≈ ∆t1



(φz l ρm )n+1 i,j,k

n l m )i,j,k

− (φz ρ



(2.25)

Note que nas Equa¸co˜es 2.24 e 2.25, todos os termos são tomados no passo de


38

tempo n + 1 com exce¸caõ de uma parcela do termo de acumula¸cão. Isso é chamado de formula¸cão totalmente impl´ıcita e demandará a solu¸cão de um sistema de Equa¸co˜es em cada passo de tempo. Sendo assim, ap´ os aplica¸caõ do método das diferen¸cas finitas às Equa¸cões 2.13 e 2.14 e considerando o tratamento do tempo de forma impl´ıcita resultam as Equa¸co˜es 2.26 e 2.27.



−

V i,j,k (φρm )n+1 (φρm )ni,j,k i,j,k ∆tn+1 n+1 (T w ∆Φw )n+1 i+1/2,j,k + (T w ∆Φw )i−1/2,j,k n+1 (T o ∆Φo )n+1 i+1/2,j,k + (T o ∆Φo )i−1/2,j,k

−

n+1 (T w ∆Φw )n+1 i,j+1/2,k + (T w ∆Φw )i,j 1/2,k −

− − − − −

n+1 (T o ∆Φo )n+1 i,j+1/2,k + (T o ∆Φo )i,j −1/2,k n+1 (T w ∆Φw )n+1 i,j,k+1/2 + (T w ∆Φw )i,j,k−1/2 n+1 (T o ∆Φo )n+1 i,j,k+1/2 + (T o ∆Φo )i,j,k−1/2 + (m ˙ o )n+1 ˙ w )n+1 i,j,k + (m i,j,k =

V i,j,k ∆tn+1



(φz w ρm )n+1 i,j,k

− (φz ρ w

n m )i,j,k

(2.26)

0

−

n+1 (T w ∆Φw )n+1 i+1/2,j,k + (T w ∆Φw )i−1/2,j,k n+1 (T w ∆Φw )n+1 i,j+1/2,k + (T w ∆Φw )i,j −1/2,k n+1 (T w ∆Φw )n+1 i,j,k+1/2 + (T w ∆Φw )i,j,k−1/2 + (m ˙ w )n+1 i,j,k =

− −

(2.27) 0

Esta discretiza¸caõ considera um bloco de coordenadas ijk e utiliza o balan¸co entre os fluxos através das faces do bloco considerando os volumes acumulados e os produzidos/injetados neste bloco por unidade de tempo. O tratamento da dependência temporal é totalmente impl´ıcito, como já mencionado. O sobrescrito n ou n + 1 indica o passo de tempo onde é calculada a propriedade e os subscritos i, j, k indica as coordenadas espaciais do bloco ou interface. Já, V i,j,k é o volume do bloco de coordenadas i, j, k e é dado na Equa¸cão 2.28.

V i,j,k = ∆xi ∆y j ∆kijk

(2.28)

Nas Equa¸cões 2.26 e 2.27, ∆Φl são as diferen¸cas de potencial das fases l = o, w entre os blocos adjacentes. ∆t é o tamanho do passo de tempo, (m ˙ l )n+1 e a vazão mássica da i,j,k ´ fase l = o, w injetada/produzida no bloco ijk ao final do passo de tempo tn+1 e T l representa a


39

transmissibilidade da fase l, calculada na Equa¸cão 2.29, onde F g representa o fator geométrico na interface e λ l =

ρl krl , µl

a mobilidade da fase l na interface.

T l = F g

ρl krl , µl

(2.29)

O fator geométrico em cada dire¸caõ pode ser calculado através da média harmônica entre as permeabilidades das células adjacentes ponderadas pelos seus comprimentos e multiplicadas pela área da se¸caõ de fluxo. Por simplicidade, considere-se o problema unidimensional cujos blocos possuem comprimentos iguais ∆x e área de se¸cão A. No caso tridimensional, o mesmo racioc´ınio pode ser aplicado. Neste caso, o fator geométrico est´ a expresso na Equa¸caõ 2.30.

(F g )i+ 1 = 2

ki+ 1 Ai+ 1 2

2

∆xi+ 1

,

(2.30)

2

Considerando por simplicidade as dimensões das células iguais, a permeabilidade k da interface pode ser calculada através da simples média harmônica entre as permeabilidades das células adjacentes conforme a Equa¸cão 2.31. Entretanto, em caso de dimensões diferentes, esta média deve ser ponderada pelas dimensões. Para mais detalhes, ver (Ertekin et al., 2011; Rodrigues and Bonet, 1991).

ki+ 1 = 2

2ki ki+1 ki + ki+1

(2.31)

Quanto a mobilidade, esta é calculada utilizando o que a literatura chama de upwind , que consiste em tomar ρ, k r e µ do bloco de maior potencial seguindo-se a dire¸cão do fluxo. A Equa¸caõ 2.32 representa este racioc´ınio na dire¸cão X, e deve se repetir nas dire¸cões Y e Z.

λl,i+ 1 ,j,k = 2

 

k

+1,j,k ρl,i+1,j,k µrl,i , l,i+1,j,k

k

ρl,i,j,k µrl,i,j,k , l,i,j,k

≥ 0

se ∆Φl,i+ 1 ,j,k 2

se ∆Φl,i+ 1 ,j,k < 0 2

(2.32)

2.3. Cálculo das vazões dos po¸cos

2.3

40

C´ alculo das vaz˜ oes dos po¸cos Os termos de produ¸cão nas equa¸co˜es de res´ıduos 2.13 e 2.14 representam a vazão

de ó leo e água produzidos ou de a´gua injetada. Para calcul´ a-las, faz-se necessário que se estabele¸ca uma rela¸caõ entre a pressão de fundo do po¸co, a pressão do bloco onde se encontra cada completa¸caõ e a vazão produzida/injetada de cada componente. Além disso, como se admite m´ ultiplas completa¸cões, é necessário relacionar a pressã o de fundo de po¸c o BHP à pressão no interior do po¸co em cada uma das completa¸cões. A rela¸caõ básica entre as pressões no interior do po¸co p wf ijk e no bloco onde este se ˙ l,ijk é que recebe o nome modelo de po¸co, onde localiza p oijk e a vazão produzida/injetada M o subscrito l = o, w representa a fase e ijk representa as coordenadas do bloco do po¸c o. O modelo de po¸co aqui considerado é descrito em (Peaceman, 1983) e necessita de um fator de proporcionalidade conhecido como ´ındice de produtividade P I l,ijk e está exposto na Equa¸caõ 2.33. ˙ ijk = P I l,ijk ρn+1 λn+1 ( pn+1 M oijk lijk ijk

n+1 wf ijk )

− p

(2.33)

A mobilidade do fluido no po¸co λ n+1 e calculada na Equa¸caõ 2.34 e depende se o ijk ´ po¸co é produtor ou injetor.

   ∑

krlijk µlijk

λijk =

f =o,w

para po¸co produtor; (2.34)

krf ijk µf ijk

para po¸co injetor.

No caso de po¸c os de m´ ultiplas completa¸co˜es, faz-se necessário calcular, a partir da pressão de fundo de po¸co medida na primeira célula do po¸co, as pressões no interior do po¸co nas demais completa¸cões. Por simplicidade, aqui só serão considerados po¸cos verticais, mas o exposto se aplica a qualquer tipo de po¸cos. Sendo assim, considere pwf ijN z a pressão de fundo de po¸co (primeiro bloco completado de baixo para cima ou, como considerado no código do simulador, N z é a camada da u ´ltima completa¸caõ) e p wf ijk a pressão no interior da completa¸caõ na camada k. A pressão no interior do po¸co nas completa¸co˜es de cima pode ser obtida adicionando-se o peso da coluna de fluido como na Equa¸caõ 2.35 .

2.3. 2.3. Calculo a´lculo das vazões oes dos po¸cos

41

N z −1 n+1 n+1 = pwf + pwf ijk ijk ijN ijN z

�

ρnr+ 1 g (hi,j,r 2

r =k

−h

i,j,r+1 i,j,r+1 ),

(2.35)

edia dos fluidos presentes no interior do po¸co co na interface ρnr+ 1 é a densidade média 2

entre os blocos i blocos i,, j, r e i, edia das densidades dos i, j, r + 1 e pode ser calculada fazendo-se a média fluidos que adentram o po¸co co desde o fundo até a célula elula considerada consider ada como na Equa¸c˜ cao aõ 2.36.

ρnr+ 1 = 2

   ∑

l=o,w

ρnwr

para po¸co co injetor; (2.36)

n co produtor. produtor. f l,r+ ρn para po¸co l,r+ 1 lr 2

Note que, no caso de po¸cos cos produtores, as densidades das fases são ao ponderadas pelo fluxo fracionário ario de cada fase no interior do po¸co. co. Define-se Define-se por fluxo fracion´ fracionário ario a razão ao entre a vazão ao volum´ volumétrica etrica de cada fase e a vazão ao volum´ volumétrica etrica total considerando-se todo o volume volume que adentrou adentrou o po¸co co desd d esdee o fundo fund o at´ a té a célula elul a cons c onsider iderada ada.. Isso Iss o é matem m atematic aticame amente nte expresso através es da Equa¸c˜ cao aõ 2.37, onde Qlk é a vaz˜ vazão ao volumétrica etri ca da fase fas e l na camada k camada k.. Nz

n = f l,r+ l,r + 1 2

∑

k=r+1 Nz

∑

k=r+1

n + Qw k

Qnlk Nz

∑

k=r+1

,

(2.37)

Qonk

Entretanto, já que as vazões oes volumétricas etri cas são ao diretamente proporcionais às as mobilidades de cada fase, sendo mais fácil acil calcular as mobilidades do que as vazões oes — principalmente quando há gás as dissolvido — o fluxo fracionário ario é aproximado aproximad o pela razão ao entre as mobilidades como exposto na Equa¸c˜ cao aõ 2.38, onde λ onde λ lk a mobilidade da fase l na camada k camada k.. Nz

n = f l,r+ l,r + 1 2

∑

k=r+1 Nz

∑

k=r+1

n + λw k

λnlk Nz

∑

k=r+1

,

(2.38)

λonk

Note que, apesar de a formula¸c˜ cão ao aqui considerada conside rada tratar tra tar a dependˆ d ependência encia tempora tem porall de forma impl´ impl´ıcita, para cálculo alculo da pressão ao no interior do po¸co co em cada completa¸c˜ cao, aõ, como na Equa¸c˜ cão ao 2.35, necessitanecessita-se se da densidade densidade média edia dos fluidos no interior interior do po¸co. c o. Para ara

2.4. 2.4. Ca´lculo do Valor Presente L´ıquido

42

calcul´ a-la a- la,, é necess nec ess´ário ario utilizar propriedades calculadas no passo de tempo anterior o que, de certa forma, altera o esquema de cálculo alcu lo totalm to talmente ente impl´ i mpl´ıcito. ıcit o. Por fim, o controle de cada po¸co co pode assumir duas formas. formas. A primeira, primeira, por vaz˜ ao ao especificada, fixa a vazão a o do po¸co co (soma de todas as completa¸c˜ coes) o˜es) e portanto tranforma a pressão ao em cada completa¸c˜ cao ão em variáveis, aveis, adicionando ao sistema de equa¸c˜ cões oes uma nova equa¸c˜ cao aõ baseada na Equa¸c˜ cao aõ 2.33 para cada uma das completa¸c˜ coes o˜es de cada po¸co co com vazão ao especificada especificada.. A segunda, por press˜ ao especificada, fixa a pressão ao ao na completa¸c˜ cao ão mais profunda p funda p wf ijN e calcula da forma fo rma descrita d escrita na Equa¸ Equ a¸c˜ cão ao 2.35 a pressão ao nas outras completa¸ completa ¸c˜ coes o˜es ijN z acima. acima. Sendo Sendo assim, assim, n˜ ao a o há necessidade de adi¸c˜ cão ao de equa¸c˜ coes o˜e s e não ao há vari´ var iáveis aveis extras. A vazão ao passa a ser calculada pela Equa¸c˜ cao aõ 2.33 onde a unica u ´ nica variável avel advinda do sistema de equa¸c˜ coes o˜es de fluxo é a pressão ao na célula elula de cada completa¸c˜ cão. O unico uńico tipo de controle considerado consider ado neste trabalho trabalh o é este est e ultimo, u ´ ltimo, por pressão ao especificada.

2.4

C´ alculo alcu lo do Valor alo r Prese Pr esente nte L´ L´ıquido ıqu ido Por fim, o simulador simulador ser´ a utilizado utilizado para a previs˜ previs˜ ao ao do Valor Presente L´ıquido

(VPL) obtidos ob tidos ao a o aplicar-se aplic ar-se uma um a estratégia egia de produ¸ prod u¸c˜ cao aõ baseada exclusivamente exclusivamente na inje¸c˜ cao aõ de agua, água, onde os controle controless dos po¸cos cos são ao as pressões oes de fundo aplicadas em cada um deles, tanto produtores como injetores. Sendo assim, nos cap´ cap´ıtulos 3, 4 e 5 o simulador será considerado como uma fun¸c˜ cao V aõ V P L = f = f ((p), onde p é o vetor de controles contendo todas to das as pressões oes de fundo. fundo. Desta forma, forma, considere considere um campo com nw po¸cos cos simulados nos passos de tempo τ t , p ode ser calculado atrav´ através es da Equa¸c˜ cao aõ 2.39. t = 1 . . . nt , o VPL pode V P L =

nt

nw

∑∑

t=1 w=1

1 (1+ j) j )τ t

× − −

·

[P Oil Oil q Oil,w Oil,w ( pw )

·

C WProd W Prod q WProd,w W Prod,w ( pw )

·

(2.39)

C WInj W Inj q WInj,w W Inj,w ( pw )]

´ importante ressaltar que a taxa de desconto j deve relacionar-se ao mesmo E per´ per´ıodo a que se relacionam vaz˜ vazões oes q w . Po Porr exem exempl plo, o, se j for uma taxa mensal, q w deve ser vazão ao mensal, isto é, e, volume produzido/injetado pro duzido/injetado em um mês. es. O pre¸co c o do óleo é P Oil Oil e os custos de produ¸c˜ cao aõ e inje¸c˜ cao a˜ o de água agua são ao respectivamente respec tivamente C WProd odos s˜ sao aõ W Prod e C WInj W Inj . Todos considerad considerados os constantes. constantes. A vaz˜ ao a o do po¸co co q w é determinada pelo simulador como a soma das vazões oes em todas as completa¸c˜ cões. oes. Para tanto as pressões oes em cada completa¸c˜ cao aõ deve ser

2.5. 2.5. Compar Compara¸ a¸c˜ caõ com simulador comercial

43

calcula calc ulada da através es da Equa¸ Equa ¸c˜ cão ao 2.35 para aplicar-se o modelo de po¸co co da Equa¸c˜ cao a˜ o 2.33 2.33.. A rotina de cálculo alculo do Valor Presente L´ıquido é externa ao simulador simulado r e aplica a Equa¸c˜ cão ao 2.39 após os receber a vazão ao de cada po¸co co em cada passo de tempo.

2.5

Compara¸ Compar a¸c˜ cao a õ com c om simulador s imulador comercial comerci al Como uma tentativa de valida¸c˜ cao aõ do simulador, buscou-se comparar os resultados

da simula¸c˜ cao aõ de um modelo de reservatórios orios obtida com o uso do SIMPAR ao resultado obtido com algum simulador comercial já consagrado pela ind´ ustria ustria.. Sendo Sendo assim, assim, por ser mais comumente utilizado na Petrobras, escolheu-se o IMEX da CMG como simulador a ser comparado. O primeiro caso teste proposto nesse estudo comparativo, e utilizado também em nos Cap´ıtulos ıtu los 3, 4 e 5, é uma por¸ po r¸c˜ cão ao modificada do modelo SPE10, um reservatório orio sintético etic o paralelepip´ edico edico (“caixa de sapato”) com canais proposto inicialmente em (Christie and Blunt, 2001). Esta modifica¸ modifica¸c˜ cao, aõ, considerando a mesma por¸c˜ cão ao aqui utilizada, foi determinada em (Cardoso, 2009) como um dos modelos para aplica¸c˜ cao aõ da técnica. ecni ca. O mod m odelo elo é tridi tr idimens mension ional al

× 80 × 5 células, elulas, isto é, e, possui pos sui 24000 células. elulas . Todas oda s as células elulas possuem poss uem dimens˜ oes o es de 50ft × 70ft × 20ft. Há 4 po¸cos cos produtores (nomeados P1-P4) e dois injetores de dimens˜ dim ensões oes 60

(I1 e I2). O valor médio edio para a permeabilidad permeabilidadee na dire¸c˜ cao aõ x é kx = 418 mD. mD. A permea permeabi bi--

lidade vertical é definida como kz = 0,3kx nos canais e kz = 10 3 kx na por¸c˜ cao aõ restante restante.. As −

satura¸c˜ coes o˜es iniciais de óleo o leo e água agua s˜ ao ao respectivamente S oi oi = 0,8 e S wc wc = 0,2 e as residuais, oleo oleo tem-se ρo = 45lb/ 45lb/ft3 e µo = 3 cp e para a água agua ρw = 60lb/ 60lb/ft3 S wr wr = S or or = 0,2. Para o ´ e µw = 0.3 cp. O sistema é pouco p ouco compress´ıvel ıvel com co m co = cw = cr = 10

6

−

psi

1

−

(compressi-

biliad biliadade adess dos fluidos fluidos e da rocha iguais). iguais). O efeito efeito da press˜ press˜ ao ao capilar é desconsiderado. As permeabilidades relativas são ao calculadas atrav´ através es do modelo de Corey (Equa¸c˜ coes ões 2.19 e 2.20 0 0 com os seguintes parâmetros: ametros: kro = krw = 1 e a = b = 2. A figura 2.2 repres represen enta ta o mapa mapa

de permeabilidades na dire¸c˜ cao aõ x, kx com os po¸cos cos posicionados. posicionados. Neste Neste teste comparativo comparativo,, todos os po¸cos cos foram controlados por pressão ao de fundo especificada, sendo para os produtores para os injetor injetores es P 11000 psi. psi. O tempo da simula simula¸c˜ çao aõ é de 3000 dias. P w = 4000 psi e para P w = 11000 Os resultados obtidos para este modelo SPE10 - Adaptado se encontram na Figura 2.3. 2.3 . A compa co mpara¸ ra¸c˜ cão ao foi feita considerando as vazões oes de óleo oleo produzido (Figura 2.3(a)) a vazão ao de água agua produzida (Figura 2.3(b)) e a vazão a o de água agua injetada (Figura 2.3(c)) por po¸co. co. Estes mostraram quase perfeita aderência encia entre as vaz˜ vazões oes calculadas pelos dois simuladores

2.5. Compara¸caõ com simulador comercial

44

Fig. 2.2: Modelo sintético (baseado no caso SPE10) com 24000 células, 5 po¸cos produtores e 2 injetores. Mapa de permeabilidades na dire¸cão x. para todos os 3000 dias de simula¸caõ. Portanto, nem a formula¸caõ por fra¸co˜es mássicas nem qualquer outra diferen¸ca que por ventura ambos os simuladores possuam, tiveram efeito significativo nos resultados. Este ent˜ ao foi um dos modelos utilizados nos estudos contidos nos Cap´ıtulos 3, 4 e 5 posteriores. Apesar disso, o desempenho do SIMPAR foi bastante inferior ao do IMEX em termos de tempo de simula¸caõ. Utilizando um computador com processador Intel Core 2 Duo 3 GHz com 4 GB de RAM e sistema operacional Windows 7, enquanto o IMEX utilizou 194 s, o SIMPAR utilizou 2175 s na simula¸caõ. Observou-se que o u ´ ltimo cortou mais os passos de tempo (em torno de duas vezes) e fez em torno de 4 vezes o número de itera¸co˜es lineares total. Isto se deve a`s configura¸co˜es numéricas, mas também a` qualidade do solver antigo e menos eficiente que é utilizado no SIMPAR. Na verdade, o solver desenvolvido no CENPES no fim da década de 80 é do tipo orthomin (Vinsome, 1976), mas seu pré-condicionador está muito ultrapassado se comparado ao estado da arte atual neste área. Este modelo SPE10 - Modificado foi bastante utilizado nos estudos de acurácia do TPWL, otimiza¸caõ e incertezas. Entretanto, além de ser um modelo pequeno em termos de n´ umero de c´ elulas, possui poucos controles. São apenas 6 po¸cos onde somente os 4 produtores serão controlados nos estudos de otimiza¸caõ. Considerando-se o n´ u mero de 6 ciclos de controle, máximo utilizado, obtém-se um problema com 24 vari´ aveis de projeto para o problema de otimiza¸caõ.


45

´ (a) Vazão de Oleo Produzido por Po¸co

´ (b) Vazão de Agua Produzida por Po¸co

´ (c) Vazão de Agua Injetada por Po¸co

Fig. 2.3: Compara¸caõ de resultados do SIMPAR com IMEX - Modelo SPE10 - Modificado


46

Por outro lado, juntamente com a ideia de revitaliza¸cão do simulador SIMPAR, surgiu uma demanda por constru¸caõ de modelos para testes. Um campo que tem sido utilizado frequentemente na literatura e tem um caráter geológico mais realista é o campo de Brugge (Peters et al., 2010), proposto como caso benchmark para problemas de otimiza¸caõ ´ um campo sint´ e ajuste ao histórico. E etico, cuja estrutura consiste em um domo alongado de leste para oeste com uma falha de borda ao norte e uma falha interna com um rejeito

× 3 km.

modesto. As dimens˜ oes do campo são de aproximadamente 10 km

A sequência

vertical das forma¸co˜es diz respeito a estratigrafia BRENT (forma¸cões Broom-Ranoch-EtiveNess-Tarbert). O modelo de reservatórios proposto para um estudo comparativo em (Peters et al., 2010) possui células de dimensões aproximadas de 75 m

× 75 m × 2, 5 m e um total de 60048

blocos, sendo 44550 ativos. Contém 20 po¸cos produtores e 10 po¸cos injetores, todos verti-

cais. Está subsaturado o que permite a simula¸cão bifásica e a constru¸cão do modelo de ordem reduzida TPWL. A figura 2.4 apresenta o mapa de permeabilidades na dire¸cão x. As propriedades de fluido, rocha e rocha-fluido estão descritas em (Peters et al., 2010). Na descri¸caõ do modelo foram fornecidas 104 realiza¸co˜es das propriedades de reservatório (permeabilidades, porosidade, razão net-to-gross , etc), tendo sido utilizada neste estudo a primeira delas. A dura¸cão da simula¸caõ considerada aqui diz respeito somente ao per´ıodo relativo ao histórico de produ¸caõ proposto no estudo comparativo: 10 anos. O SIMPAR, por ser um simulador antigo, possui algumas limita¸cões em rela¸caõ aos mais modernos. Sendo assim, para comparar seus resultados a respeito deste modelo com um outro simulador (IMEX), fez-se necessário realizar algumas modifica¸co˜es que deveriam ser replicadas nos modelos para os dois sumuladores. A primeira diz respeito à descri¸caõ da malha de simula¸caõ, não possibilitando a utiliza¸cão de malhas mais genéricas, ficando limitado ao tipo cartesiano com topo variável. O caso em questão possui uma malha do tipo corner point que precisou ser convertida em cartesiano, utilizando-se a ferramenta Cyclope desenvolvida em parceria Petrobras-ESSS (CENPES and ESSS, 2012). Esta conversão (juntamente com as outras modifica¸co˜es descritas a seguir) alterou consideravelmente o volume poroso (redu¸caõ de 16 %) e os volumes de óleo in place (aumento de 5%) e água in place (redu¸cão de 20 %), pois causa uma deforma¸cão da malha nas partes periféricas, onde as células da malha original são muito diferentes de paralelep´ıpedos. Os volumes de rocha e fluido do modelo antes da conversão encontram-se na primeira coluna da tabela 2.1 para compara¸cão com o caso modificado. Esta diferen¸ca se deve ao fato de


47

a conversão simplesmente fazer um mapeamento um para um das c´ elulas das duas malhas (original e convertida) mantendo assim a quantidade de células. O mais correto seria fazer um recobrimento do espa¸co representado no modelo com a nova malha cartesiana minimizando as diferen¸cas de volumes, tarefa que pode vir a ser feita no futuro.

Fig. 2.4: Modelo de Brugge Modificado. Mapa de permeabilidades na dire¸cão x. Além disso, houve a necessidade de se redefinir o desenho dos po¸cos. Todos no modelo original são verticais, mas devido ao grid não regular não possuem todas as completa¸co˜es com as mesmas coordenadas x e y. Entretanto, com a “regulariza¸caõ” da malha, os po¸cos deixariam de ser verticais se fossem mantidas as coordenadas originais para as completa¸cões. Sendo assim, alteraram-se as posi¸cões das completa¸cões e desprezou-se a descri¸caõ exata da trajetória (este u ´ ltimo recurso não existe no SIMPAR). Outras modifica¸co˜es foram realizadas tais como: (1) simplifica¸caõ do modelo de fluido, pois o SIMPAR não admite tabelas PVT para um modelo bi-fásico; (2) raz˜ a o de solubilidade anulada, devido a diferen¸ca de formula¸caõ; (3) reposicionamento da profundidade de referência, pois o SIMPAR só permite que esteja no contato; (4) pressão capilar desprezada devido a diferen¸cas na forma de inicializa¸caõ que geravam discrepâncias entre os volumes iniciais de óleo e água calculados pelos simuladores; (5) considera¸caõ de uma u ńica curva de


48

permeabilidade relativa para todo o campo devido a complexidade de convers˜ ao do esquema de re-escalamento das curvas utilizado no modelo IMEX. Devido a todas estas altera¸cões no modelo, este caso é chamadado daqui por diante de Brugge - Modificado. Todos os po¸cos produtores foram controlados por vazã o máxima de l´ıquido de 320m3 /d, exceto os po¸c os P10 e P15, onde a vazã o máxima de l´ıquido foi alterada para 220m3 /d e o po¸co P9, que foi controlado por pressão m´ınima de fundo de 51 kgf /cm2 . Os po¸cos injetores foram controlados por vaz˜ ao máxima de água injetada de 636 m3 /d. Desta forma, todas estas modifica¸ co˜es foram aplicadas tanto ao modelo simulado no IMEX quanto ao simulado pelo SIMPAR, possibilitando uma compara¸caõ justa dos resultados para fins de valida¸caõ. A Figura 2.5 contém as vazões de óleo, água e l´ıquido produzidos e vazão de água injetada total do campo obtidas pelos dois simuladores, utilizando-se a primeira realiza¸caõ do modelo de Brugge, com todas as modifica¸cões descritas.

´ (a) Vazão Total de Oleo Produzido

(c) Vazão Total de L´ıquido Produzido

´ (b) Vazão Total de Agua Produzida

´ (d) Vazão Total de Agua Injetada

Fig. 2.5: Compara¸caõ de resultados do SIMPAR com IMEX As vazões totais produzidas e injetadas se mostraram praticamente idênticas nos resultados obtidos pelos dois simuladores. Houve uma pequena diferen¸ca na produ¸cã o de


49

óleo e de água devida à pequenas discrepâncias após a erup¸cão de água em alguns dos po¸cos e a uma discrepância razoável nas vazões calculadas para o po¸co controlado por pressão de fundo. Em todos os po¸cos produtores, exceto neste último, a vazão de l´ıquido máxima é atingida durante todo o tempo. As vaz˜ oes de inje¸cã o de todos os po¸cos se mativeram no valor máximo nos dois resultados. Outro importante resultado a ser comparado, calculado na inicializa¸caõ do simulador, são os volumes de ó leo e água in place e o volume poroso. Como ambos os modelos possuem a mesma malha e mesmos mapas espera-se que estes coincidam. A Tabela 2.1 apresenta os resultados. Pode-se observar que os volumes coincidem, mas isto só pode ser obtido desprezando-se a pressão capilar. Diferen¸cas no cálculo da inicializa¸caõ (posicionamento dos contatos) dos dois simuladores causavam discrepâncias nos volumes de óleo e água quando a curva de pressão capilar foi considerada. Tab. 2.1: Compara¸ cão

de volumes in place calculados pelos simuladores

IMEX - Inicial Volume Poroso(m3 ) Volume de óleo (Sm3 ) Volume de água (Sm3 ) Número de células ativas

8

× 10 1,318 × 10 7,589 × 10 8,875

8 8

44550

Simpar

IMEX 8

× 10 1,387 × 10 6,028 × 10 7,380

8 8

44550

8

× 10 1,388 × 10 6,027 × 10 7,380

8 8

44550

Apesar da proximidade dos resultados, assim como com o modelo SPE10 - Modificado, o desempenho do SIMPAR com o Brugge - Modificado foi bastante inferior ao do IMEX em termos de tempo de simula¸cão. Utilizando um computador com processador Intel Core 2 Duo 3GHz com 4GB de RAM e sistema operacional Windows 7, enquanto o IMEX utilizou 247s, o SIMPAR utilizou 1477s na simula¸caõ. Observou-se que o u ´ ltimo fez em torno de 1,4 vezes o número de itera¸co˜es lineares total. Entretanto, o acréscimo de tempo foi muito mais expressivo. Comparando com o primeiro modelo, demonstra-se que o solver do SIMPAR não é bem escalável com o n´ umero de células da simula¸caõ. O modelo de Brugge - Modificado foi utilizado nos estudos de otimiza¸cã o do cap´ıtulo 4 com mais uma modifica¸cão: os controles utilizados foram somente pressõ es de fundo de po¸co para que fosse compat´ıvel com o modelo de ordem reduzida.

3. MODELO DE ORDEM REDUZIDA Com o advento de modelos de simula¸caõ cada vez maiores e mais detalhados e com a necessidade de estudos que demandam cada vez mais simula¸co˜es destes modelos, como otimiza¸caõ de controles e posi¸cão de po¸cos, análise de incertezas ou, ainda mais, otimiza¸caõ sob incerteza, cresce a necessidade de constru¸cão de modelos substitutos de ordem reduzida que sejam capazes de obter resultados próximos aos dos simuladores com custo computacional bem inferior. Este cap´ıtulo descreve os aspectos teóricos acerca do modelo substituto de ordem reduzida utilizado neste trabalho, bem como os detalhes para a sua implementa¸cão. A técnica escolhida para a redu¸cão da complexidade numérica é conhecida como Trajectory Piecewise Linearization - TPWL e consiste basicamente na lineariza¸caõ da equa¸caõ do res´ıduo em torno de estados salvos durante uma simula¸caõ previamente executada, também chamada de simula¸cão de treinamento. Tal técnica foi proposta pela primeira vez em (Rewienski, 2003) num contexto de microeletrônica, mas foi aplicada inicialmente à problemas de simula¸caõ de reservatórios por Cardoso (Cardoso, 2009). A simples aplica¸caõ desta técnica transforma o problema não linear da simula¸caõ de escoamento em uma sequência de solu¸co˜es de sistemas lineares. A dimensão deste problema ainda é a mesma do problema original. Portanto, faz-se necess´ ario uma técnica para redu¸caõ da dimensão do problema. A técnica de redu¸caõ escolhida é conhecida como Proper Orthogonal Linearization - POD ou Principal Component Analysis - PCA amplamente aplicada em diversas áreas da engenharia, tendo sido combinada com o TPWL já em (Cardoso, 2009). Por outro lado, há um grande número de trabalhos que mencionam um problema de estabilidade com a técnica (He, 2010; Gildin et al., 2013). Estudos mais recentes mostram

50

3.1. Lineariza¸cão utilizando TPWL

51

que a técnica pode ser estabilizada pela substitui¸cã o do tipo de proje¸cã o que era usado inicialmente (Petrov-Galerkin) pela proje¸caõ de Bubnov-Galerkin (Carlberg et al., 2009; He and Durlofsky, 2013). Isto foi utilizado em todos os estudos aqui apresentados, eliminando problemas com estabilidade do método. A primeira parte deste cap´ıtulo consiste em explica¸co˜es sobre a matemática e o funcionamento algor´ıtmico do TPWL. A se¸caõ 3.1 apresenta os detalhes do método de lineariza¸caõ aplicado à equa¸caõ de fluxo. Em seguida, a se¸caõ 3.2 apresenta em detalhes a forma das matrizes de derivadas que precisarão ser exportadas a partir do código do simulador de reservatórios. A se¸caõ 3.3 detalha metodologia de redu¸caõ de dimensão por POD (Proper Orthogonal Decomposition ). A u ´ ltima se¸cão desta primeira parte, se¸caõ 3.4, descreve o modelo de ordem reduzida obtido pela combina¸caõ da lineariza¸cão com o TPWL e da redu¸caõ de dimensão com o POD. Na segunda parte são expostos a implementa¸caõ do método acoplado ao simulador SIMPAR, detalhando a constru¸cão de rotinas de exporta¸cão, a forma de comunica¸ca˜ o e o desenvolvimento utilizando o Matlab na se¸cão 3.5. Por fim, apresentam-se alguns resultados de medida de acurácia do método na se¸caõ 3.6, comparando valores obtidos para o Valor Presente L´ıquido obtidos com o simulador e com o modelo de ordem reduzida. Isto ilustra a sua qualidade e evidencia algumas das limita¸co˜es.

3.1

Lineariza¸ c˜ ao utilizando TPWL

Para construir o modelo reduzido, primeiro é necessário analisar as equa¸cões de ´ poss´ıvel reescrevê-las, dividindo-as em três partes: termos de fluxo descritas em 2.26 e 2.27. E acumula¸caõ, que representa toda a massa acumulada no bloco e depende da compressibilidade dos fluidos, termo de fluxo, que representa a massa que flui de um bloco para outro e depende das transmissibilidades das interfaces, e termo de fonte/sumidouro, que representa a massa produzida/injetada em eventuais po¸cos completados nos blocos e depende do modelo de po¸co. A Equa¸cão 3.1 representa estas três parcelas que somadas resultam em zero, fechando assim o balan¸co de massa.

g(xn+1 , xn , un+1 ) = A(xn+1 , xn ) + F(xn+1 ) + Q(xn+1 , un+1 ) = 0

(3.1)


52

Na Equa¸caõ 3.1, g representa o vetor de res´ıduo para as duas equa¸co˜es que representam o modelo bifásico. Como a equa¸caõ é não linear, utiliza-se o método de NewtonRaphson para a solu¸caõ. Os ´ındices n e n + 1 designam o passo de tempo da simula¸ caõ e A ,

F e Q são os termos de acumula¸caõ, fluxo e fonte/sumidouro respectivamente aqui representados em negrito por representarem vetores com uma componente para cada equa¸caõ, sendo duas equa¸cões para cada bloco. A técnica conhecida como Trajectory Piecewise Linearization - TPWL se baseia na lineariza¸cão da equa¸caõ do res´ıduo em torno de estados salvos durante uma simula¸caõ previamente executada, tamb´ em chamada de simula¸caõ de treinamento. A cada passo de tempo lineariza-se o modelo em torno de um ponto de alguma trajetória de treinamento, cujos estados foram previamente convergidos e armazenados (Cardoso, 2009; He, 2010). Analisando-se a equa¸cão 3.1 sabe-se que o res´ıduo g é fun¸caõ do vetor de estados do passo atual xn+1 , do vetor de estados do passo anterior xn e do vetor de controles dos po¸cos un+1 . Sendo assim, considere os vetores de estados xi+1 e xi convergidos e memorizados durante uma simula¸caõ de treinamento utilizando o vetor de controles u i+1 . Seja a expansão em série de Taylor do res´ıduo g n+1 representado na equa¸caõ 3.1 em torno de x i+1 , x i e u i+1 . Podemos aproximar o res´ıduo em uma simula¸caõ g n+1 a partir do estado calculado no passo anterior xn , utilizando os dados memorizados pelo truncamento da expansão, representado na equa¸cão 3.2.

gn+1 = g i+1 +

∂ gi+1 n+1 (x ∂ xi+1

−x

i+1

∂ gi+1 n )+ (x ∂ xi

i+1

− x ) + ∂ ∂ ug i

i+1

(un+1

i+1

−u

),

(3.2)

´ importante ressaltar que na Equa¸cão 3.2, o sobrescrito i representa um passo E de tempo de uma simula¸caõ de treinamento, onde as derivadas, estados e controles utilizados foram convergidos e memorizados e o sobrescrito n representa os passos de tempo de uma nova simula¸cão, considerando o mesmo estado inicial, mas outros controles. Além disso, s˜ ao feitas as seguintes simplifica¸cões de nota¸caõ g n+1 = g (xn+1 , xn , un+1 ) e g i+1 = g (xi+1 , xi , ui+1 ). Por outro lado, gi+1 = 0, pois representa o res´ıduo previamente convergido durante a simula¸cão de treinamento. Por defini¸cão a matriz jacobiana é a derivada do res´ıduo em rela¸cão ao estado. Desta forma a primeira derivada parcial que aparece na Equa¸caõ 3.2 pode ser renomeada conforme Equa¸caõ 3.3.


53

i+1

J

∂ gi+1 = ∂ xi+1

(3.3)

As outras duas derivadas parciais que aparecem na equa¸cão do res´ıduo podem ser reescritas. Isto se deve ao fato de que o res´ıduo é a soma de fatores que não dependem todos da tripla x i+1 , xi , ui+1 . Expandindo o res´ıduo conforme Equa¸cão 3.2 no cálculo de sua derivada em rela¸cão ao estado x i , resulta a Equa¸caõ 3.4. ∂ gi+1 ∂ = A(xi+1 , xi ) + F(xi+1 ) + Q(xi+1 , ui+1 ) i i ∂ x ∂ x





(3.4)

Como somente o termo de acumula¸cão A (xi+1 , xi ) possui dependência em rela¸caõ ao estado no passo de tempo anterior, a Equa¸caõ 3.4 pode ser simplificada como na Equa¸caõ 3.5. ∂ gi+1 ∂ Ai+1 = ∂ xi ∂ xi

(3.5)

` primeira vista, a derivada na Equa¸caõ 3.5 parece não estar dispon´ıvel facilmente A em um simulador. Entretanto, ao observar-se as Equa¸co˜es 2.26 e 2.27, nota-se que o termo de acumula¸caõ é composto por uma diferen¸ca entre duas parcelas quase idênticas. A primeira tomada no passo atual (i + 1) e a segunda no passo anterior (i). Isto posto, pode-se expressar uma rela¸cão entre a derivada que aqui se deseja calcular

∂ Ai+1 ∂ xi

e a derivada do termo de

acumula¸caõ no passo anterior em rela¸caõ ao vetor de estados também do passo anterior Para tanto, deve-se fazer uma corre¸cão conforme se expressa na Equa¸caõ 3.6, onde ∆t

∂ Ai . ∂ xi i+1

e

∆ti representam a dura¸caõ do passo de tempo atual e anterior respectivamente. ∂ Ai+1 = ∂ xi

i

∂ Ai i+1 ∂ xi

− ∆t∆t

(3.6)

A mesma simplifica¸cão realizada na derivada do res´ıduo em rela¸cão ao vetor de estados x i e que resulta na equa¸cão 3.5 pode ser repetida na derivada em rela¸cão ao vetor de controles u i+1 . Como somente o termo de fonte/sumidouro Q(xi+1 , ui+1 ) possui dependência em rela¸caõ ao controle, a derivada

∂ gi+1 pode ∂ ui+1

ser simplificada como na Equa¸caõ 3.7.

∂ gi+1 ∂ Qi+1 = . ∂ ui+1 ∂ ui+1

(3.7)


54

Como o res´ıduo gn+1 = 0, a Equa¸caõ 3.2 pode ser reescrita de forma mais compacta na Equa¸caõ 3.8. Esta é a forma da equa¸cão motriz do TPWL expressa em (He, 2010).

Ji+1 (xn+1

i+1

−x

)=

�−

∂ Ai+1 n (x ∂ xi

−

∂ Qi+1 n+1 (u xi ) + ∂ ui+1

i+1

−u

�

)

(3.8)

Observando a Equa¸cão 3.8 percebe-se que esta é uma forma de relacionar o desvio entre os vetores de estado simulado pelo TPWL e gravado na simula¸cão de treinamento (xn+1

i+1

−x

) com os desvios entres os estados simulados e gravados para o passo de tempo

anterior (xn xi ) e os desvios entre os controles aplicados na simula¸caõ TPWL e na simula¸caõ

−

de treinamento ( un+1

i+1

−u

).

Além disso, observa-se que o TPWL aproxima a simula¸caõ eliminando a necessidade de um processo iterativo de solu¸cão. Isto fica mais claro ao explicitar-se o estado xn+1 na Equa¸caõ 3.8 resultando na Equa¸caõ 3.9.

n+1

x

= x

i+1

(−  � J

i+1

1

−

∂ Ai+1 n (x ∂ xi

−

∂ Qi+1 n+1 (u x )+ ∂ ui+1 i

−u

i+1

)

�

(3.9)

Mas ainda há um ponto em aberto na Equa¸caõ 3.9: como escolher os estados gravados durante a simula¸cão de treinamento xi e xi+1 e as derivadas Ji+1 ,

∂ Ai+1 ∂ xi

e

∂ Qi+1 ∂ ui+1

a

serem utilizados a cada passo de tempo da simula¸caõ TPWL? Na verdade, basta determinar-se o ´ındice i do estado xi . A partir dele, determinase o estado consecutivo na simula¸caõ de treino que será xi+1 e as derivadas correspondentes a este u ´ ltimo. Mas como determinar xi ? Este será o estado “mais próximo”, segundo alguma norma, ao estado calculado pelo TPWL no passo de tempo anterior xn . Existem várias maneiras de se medir tal proximidade. A maneira escolhida é parecida com a definida em (Cardoso, 2009) e também utilizada em (He, 2010). Seja zn o mapa de fra¸co˜es mássicas no passo de tempo n da simula¸caõ TPWL, representado por somente as linhas de x n correspondentes às fra¸cões mássicas. O estado gravado xi será aquele cujo mapa de satura¸co˜es zi for mais próximo de zn medindo-se a distância pela norma euclidiana. Sendo assim, determina-se o ´ındice i do estado gravado a ser utilizado no passo de tempo n minimizando-se a distância no espa¸co das fra¸co˜es mássicas como mostrado na Equa¸cão 3.10.

i = arg min zk k

n

| −z |

(3.10)


55

Esta lineariza¸caõ é melhor se o estado simulado for mais próximo do estado gravado mais próximo. Além disso é dependente da trajetória, isto é, dos estados pelos quais passou anteriormente. Sendo assim, existe em torno de cada estado da trajetória de treinamento, uma região de validade para o modelo linearizado dentro da qual a aproxima¸caõ é boa. Se a trajetória simulada passar por dentro de tais regiões ent˜ ao a aproxima¸caõ será melhor. A figura 3.1 foi extra´ıda de (Albunni et al., 2008) esquematiza este racioc´ınio.

Fig. 3.1: Região Validade do TPWL em torno da trajetória de treinamento Desta forma, considere uma sequência de controles (pressões de fundo de po¸co) t N c t=1 ,

{BHP }

onde N c é o n´ umero de ciclos de controle ao longo do tempo de simula¸caõ

T para os quais se desejam calcular as vazõ es dos po¸cos numa simula¸cã o do modelo de ordem reduzida TPWL. Considere uma simula¸cão de reservatórios chamada de simula¸caõ de treinamento. Esta simula¸cão passa por N t passos de tempo. Armazenam-se e exportam-se todas as seguintes informa¸co˜es ao final de todos os passos de tempo da simula¸cã o após a convergência do método de Newton: i N t i=1

• os incrementos de tempo da simula¸caõ {∆t } ; • os controles utilizados {u } ; • os estados convergidos {x } ; • as derivadas convergidas {J } , { } e { i N t i=1

i N t i=1

∂ Ai+1 N t i=1 ∂ xi

i+1 N t i=1

∂ Qi+1 N t ; ∂ ui+1 i=1

}

O algoritmo executado na simula¸caõ TPWL consiste nos seguintes passos: 1. Faz-se n = i = 1, tn = 0 e x n = x i . 2. Determina-se o vetor de controles do passo de tempo de un+1 a partir da lista de controles a serem simulados BHPt

{

N c t=1 .

}

Para po¸cos com múltipla completa¸caõ, faz-se

necess´ ario calcular a pressão nas células superiores a` de fundo utilizando o modelo de

3.2. Exporta¸caõ de mapas de estados e derivadas

56

po¸co da Equa¸cão 2.35. 3. Busca-se entre todos os estados gravados xi

N t i=1 qual

{ }

é mais próximo de xn conside-

rando somente as fra¸cões mássicas. Determina-se i e selecionam-se os estados gravados

xi e x i+1 e o intervalo de tempo entre os passos de tempo i e i + 1, ∆ti . 4. Selecionam-se também as derivadas Ji+1 ,

∂ Ai+1 ∂ xi

e

∂ Qi+1 . ∂ ui+1

5. Calcula-se o vetor de estados x n+1 Resolve-se a Equa¸cão 3.9 com os dados selecionados e o estado anterior xn . 6. Determina-se o tempo transcorrido até então na simula¸caõ por tn+1 = t n + ∆ti . Além disso, atualiza-se o vetor de estados, x n

←x

n+1

.

7. Calculam-se as vaz˜ oes nos po¸cos aplicando-se as pressões das células dos po¸cos contidas em x n+1 ao modelo de po¸co da Equa¸caõ 2.33. 8. se t n+1 < T retorna-se ao passo 2, caso contrário a simula¸caõ termina. ´ importante ressaltar que até agora não houve nenhuma redu¸caõ de dimensão E do problema. Tal redu¸caõ será realizada pela aplica¸caõ de uma técnica descrita na se¸caõ 3.3 conhecida como Proper Orthogonal Decomposition e modifica minimamente o algoritmo do TPWL.

3.2

Exporta¸ ca õ de mapas de estados e derivadas A Equa¸caõ 3.9 demonstra que o TPWL basicamente transforma a solu¸caõ de um

sistema de equa¸cões não lineares que modela o problema de escoamento em uma sequência de sistemas lineares. Entretanto, para que o método funcione, é necessário conhecer os mapas de pressã o e fra¸cõ es mássicas em cada passo de tempo da simula¸caõ de treinamento que constituirão juntos os vetores de estado xi . Além disso, também é necessário conhecer as matrizes de derivadas convergidas em todos os passos de tempo, Ji+1 ,

∂ Ai+1 ∂ xi

e

∂ Qi+1 . ∂ xi+1

Portanto, fez-se necessário alterar o código do SIMPAR de forma a que ele fosse capaz de exportar todos estes dados enquanto executava a simula¸cão. Algumas rotinas foram criadas, de uma forma menos intrusiva poss´ıvel, para que, a cada passo de tempo da simula¸caõ de fluxo, as matrizes de derivadas e mapas fossem exportadas para arquivo texto para serem lidos pelo programa TPWL. Nenhuma parte do cálculo do simulador foi alterada, mas o código precisou ser compreendido em um n´ıvel razoalvelmente profundo para acessarem-se


57

todas as variáveis necessárias e replicarem-se alguns trechos na exporta¸cão e na posterior importa¸cão no Matlab para treinamento. Os mapas de pressão e fra¸co˜es mássicas já são calculados a toda itera¸cão do método de Newton-Raphson. Inclusive já existe a op¸cão de sua exporta¸cão para pós-processamento. Sendo assim, bastou fazer a escrita de todos os valores em ordem predefinida (sequência de células definida pelo simulador) após a convergência em cada passo de tempo. A exporta¸cão dos mapas de controle possui uma complica¸caõ adicional. Basicamente há um controle, ou seja, uma pressão fixada, em cada bloco com completa¸caõ de algum po¸co. Entretanto, o que se controla efetivamente é a pressão no fundo do po¸co (BHP - Bottom Hole Pressure ). Portanto, exporta-se o valor do BHP na c´ elula mais profunda correspondente a cada po¸co, mas deve-se também exportar a pressão em fluxo, calculada na equa¸caõ 2.35, para as células superiores de cada po¸co. Ent˜ ao, aproveitou-se todo o cálculo do peso da coluna de fuido já realizado para obter-se as pressões em fluxo de todos os blocos dos po¸cos e exportá-las como controles. Em se tratando das derivadas, a mais simples de exportar foi a matriz jacobiana, uma vez que ela já deve ser montada para a realiza¸caõ da simula¸caõ. Para sua exporta¸caõ, bastou copiar a rotina de montagem do jacobiano trocando a escrita na memó ria por um arquivo de sa´ıda. Esta rotina, como todas as outras de exporta¸cão de derivadas, só pode ser chamada ap´ os a convergência do método de Newton-Raphson. ´ importante ressaltar que devido a natureza do método de diferen¸cas finitas e E da malha estruturada, esta matriz, no caso tridimensional, é praticamente heptadiagonal por blocos (a excessão ocorre quando há conexão entre células que não são vizinhas naturalmente e isto acontece em casos de pinch-out em que uma camada desaparece, havendo a necessidade de conectar duas camadas não consecutivas ou no caso de falhas com rejeito, onde parte do reservatório se desloca para cima ou para baixo). Sendo assim, é uma matriz esparsa e portanto não há a necessidade de exportar todos os seus elementos. Além disso, como no problema aqui em questão, só se admite controle dos po¸cos por pressão, esta matriz não possui as linhas e colunas correspondentes às equa¸cões de po¸co e seus elementos não-nulos são sempre os mesmos ao longo do tempo. Portanto, exportam-se as coordenadas desses elementos uma u ´ nica vez, só sendo necessário exportar os valores em cada passo de tempo. Esse mesmo racioc´ınio se aplica a`s outras derivadas. As outras derivadas normalmente não precisariam ser calculadas durante a simula¸cão. Sendo assim, fez-se necessário a contru¸cão de rotinas para cálculo e exporta¸caõ.


58

Portanto, antes de tentar implementa-las, foi necessário calcular qual seria a forma dessas derivadas baseados em propriedades e derivadas que já estivessem dispon´ıveis no decorrer da simula¸caõ. A Equa¸caõ 3.11 contém um bloco da derivada do termo de acumula¸cão em ´ importante ressaltar que este rela¸caõ ao estado anterior para uma célula N , no passo i. E bloco possui duas linhas, uma para cada equa¸caõ de balan¸co e duas colunas, uma para cada componente (ó leo e água). A matriz

∂ Ai+1 ∂ xi

para a formula¸caõ totalmente impl´ıcita é bloco-

diagonal. Isto se deve ao fato de o termo de acumula¸cão da célula N só depender de termos relacionados a` própria célula, reduzindo o número de elementos a serem exportados.

� � ∂ Ai+1 ∂ xi

= N

− ∆tV

N i+1

O cálculo desta derivada

 

∂φ ρm ∂p o

m + φ ∂ρ ∂p o

m φ ∂ρ ∂z w

∂φ m m + φz w ∂ρ + φρm z w ρm ∂p φz w ∂ρ ∂p o ∂z w o

∂ Ai+1 ∂ xi

 

i

(3.11) N

já estava em parte implementado no SIMPAR,

pois é um dos componentes necessários para o método adjunto, parcialmente implementado no simulador. Note que, assim como exposto na Equa¸caõ 3.6, tal derivada depende de propriedades calculadas no passo de tempo n e n + 1. Portanto, optou-se por exportá-la no passo de tempo n +1 junto com as outras matrizes e sem a divisão por ∆tn+1 . O seu uso deve ser então atrasado no tempo o que implica que o passo de tempo utilizado no denominador seria o imediatamente posterior. Por fim, uma outra rotina foi desenvolvida para exportarem-se as derivadas do termo de fonte/sumidouro. A partir das Equa¸cões 2.26 e 2.27 tais derivadas puderam ser calculadas. Seu c´ alculo baseia-se no modelo de po¸co utilizado exposto na Equa¸caõ 2.33. A Equa¸caõ 3.12 contém um bloco da matriz para a célula N que contenha um po¸co produtor no passo de tempo i. Nesta equa¸caõ, o termo F g é o fator geométrico do po¸co (Peaceman, 1983). Note que um bloco desta matriz possui duas linhas, uma para cada equa¸caõ, representando, a primeira, a vazão total de fluido produzido ou injetado e, a segunda, a vazã o de á gua. H´ a uma u ´ nica coluna, pois só h´ a um controle no máximo por célula. Os u ´ nicos blocos não nulos são os correspondentes às células com po¸cos completados, portanto, esta matriz é bloco-diagonal, mas com pouqu´ıssimos termos n˜ ao nulos.

3.3. Redu¸caõ de dimensão - Proper Orthogonal Decomposition (POD)

 � � −  ∂ Qi+1 ∂ ui+1

˙o ∂ M ∂p wf

+

˙w ∂ M ∂p wf

=

N

˙w ∂ M ∂p wf

 

i+1

F g

=

N

 � −  � � ρo kµroo

59

+ ρw kµrw w

F g ρw kµrw w

� 

i+1

(3.12) N

A equa¸caõ 3.13 cont´ em um bloco da matriz para a célula N que contenha um po¸co injetor no passo de tempo i. Note que neste caso as linhas correspondentes às duas equa¸co˜es são idênticas, já que a mobilidade da água num po¸co injetor é a total para os dois fluidos, como definido na Equa¸caõ 2.34.

 � � −  ∂ Qi+1 ∂ ui+1

3.3

˙w ∂ M ∂p wf

=

N

˙w ∂ M ∂p wf

 

i+1

F g ρw

= N

 � −  � F g ρw

kro µo

kro µo

+ +

krw µw

krw µw

�  �

i+1

(3.13) N

Redu¸c˜ ao de dimens˜ ao - Proper Orthogonal Decomposition (POD) A metodologia TPWL consiste em um redutor de complexidade numérica pois

lineariza o problema, eliminando a necessidade de um processo iterativo como NewtonRaphson para a solu¸cão das equa¸cões em cada passo de tempo. Entretanto, a quantidade de variáveis e equa¸cões permanece a mesma, não havendo portanto redu¸cão da dimensão do problema. Sendo assim, faz-se necessário a combina¸cão do TPWL com alguma técnica eficiente para a redu¸caõ de dimensão de modo a obter-se um modelo de ordem reduzida com consequente redu¸cão significativa no tempo computacional. A técnica assim escolhida para promover a redu¸cão de dimensão é conhecida como Proper Orthogonal Decomposition - POD ou Principal Component Analysis - PCA (Jolliffe, 2002). Ela tira proveito da correla¸cão existente entre elementos de um mesmo conjunto de dados. Pode ser aplicada em diversos ramos da engenharia. Esta técnica é amplamente utilizada em problemas de compressão, reconstru¸caõ e reconhecimento de faces devido à correla¸caõ existente entre as faces humanas (Muller et al., 2004). Por exemplo, no contexto de identifica¸cão automática de imagens, ao considerar-se um conjunto de rostos de pessoas, apesar de todas diferirem entre si, há caracter´ısticas presentes


60

em todas elas. Todas as pessoas possuem dois olhos, duas orelhas, uma testa, uma boca, etc, mais ou menos nas mesmas posi¸co˜es. Sendo assim, as imagens possuem muita correla¸caõ umas com as outras, o que sugere que o “espa¸co vetorial” gerado pelo conjunto de imagens de faces humanas tem uma dimensão inferior ao número de imagens consideradas, , se um pequeno erro for tolerado. Correla¸cões similares ocorrem, por exemplo, entre os mapas de pressão e satura¸caõ resultantes de uma simula¸caõ de reservatórios. Isto se deve ao fato de formarem-se frentes de avan¸co que seguem uma tendência fortemente dependente das caracter´ısticas do meio simulado. Este trabalho visa avaliar a reconstru¸caõ de tais mapas utilizando-se uma quantidade menor de informa¸caõ. A ideia é descobrir as dire¸co˜es principais do espa¸co gerado pelos mapas conhecidos, mas só considerar as mais significativas com base em um crit´ erio de energia. A ideia é baseada na análise espectral da matriz de covariância cujos autovalores indicam a variabilidade em cada uma das dire¸co˜es principais, que são os seus autovetores. Dire¸cões principais com baixa variabilidade (pequenos autovalores) podem ser desprezadas. Entretanto, nã o há a necessidade do cálculo da matriz de covariância; isto seria bastante custoso. A teoria demonstra a equivalência entre esse conceito e a decomposi¸cão em valores singulares, fatora¸cão de matrizes bastante conhecida pela álgebra linear.

3.3.1 Descri¸caõ do POD Considere um sistema de dimensão N , no caso deste trabalho, representando as equa¸co˜es resolvidas pelo simulador de reservatórios e descritas na Equa¸caõ 2.26 e 2.27 e um M i=1 ,

{ }

conjunto de M solu¸cões desse sistema, isto é, vetores de dimensão N , xi

chamados de

snapshots e que foram gerados por simula¸caõ num´ erica. Monta-se a matriz X da Equa¸caõ

× M em que as colunas representam os M snapshots .

3.14 de dimensão N

X = [x1 , x2 ,

·· · , x

M ]

(3.14)

Como o modelo é uma aproxima¸cão discreta de uma equa¸cão diferencial parcial, a dimensão do snapshot é normalmente muito maior que a sua quantidade (N >> M ). O υ i=1 de

{ }

objetivo do POD é determinar υ(υ < M << N ) vetores próprios ortonormais φi

tal

modo que a distância entre os snapshots e sua proje¸caõ no subespa¸co gerado pelos vetores próprios apresentada na Equa¸caõ 3.15 seja minimizada para um dado υ (Markovinovic, 2003).


1 Q(Φ) = M

M

�

xi

i=1

T

− ΦΦ

xi

dimens˜ ao n

× υ apresentada na Equa¸caõ 3.16. Φ=

�

φ1 φ2

·· ·

φυ

υ e i=1 ´

{ }

O subespa¸co gerado pelos vetores próprios φi



61

(3.15)

definido pela matriz Φ de

�

(3.16)

considerando υ autovetores da matriz n

υ i=1

{ }

Pode-se mostrar que a solu¸caõ deste problema de minimiza¸cão é dada por φi

× n dada pela Equa¸caõ 3.17 e que representa a

covariˆ ancia entre os snapshots . Estes autovetores correspondem aos υ maiores autovalores υ i i=1 .

{λ }

R =

1 XXT M

(3.17)

Em alguns casos, ao inv´ es de utilizar-se a matriz de snapshots X diretamente no ˜ conforme cálculo da matriz de covariância R definida na Equa¸caõ 3.17, calcula-se a matriz R Equa¸caõ 3.18. ˜ = 1 X ˜X ˜ T = 1 R M M

M

� i=1

(xi

T

− x¯)(x − x¯ ) i

(3.18)

˜ é apresentada na Equa¸cão 3.19, onde x ¯ representa o snapshot médio A matriz X ¯ = (1/M ) x

∑

M i=1

xi . ˜ = [x1 X

− x¯, x − x¯, ··· , x − x¯ ] , 2

M

(3.19)

Geometricamente, a subtra¸caõ da média de um conjunto de pontos de cada um dos pontos do conjunto é equivalente a mover o centro de massa do conjunto para a origem. Nas aplica¸co˜es de constru¸cão de modelos de ordem reduzida determin´ısticos, esta transla¸caõ dos pontos transformando-os em um conjunto de média zero não é de grande importância (Markovinovic, 2003) e portanto não será considerada neste trabalho. Determinar os autovetores da matriz de covariância significa descobrir as dire¸cões ortogonais nas quais os dados apresentem maior variabilidade. Ao desprezarem-se algumas


62

dire¸co˜es, as dire¸co˜es com pouca variabilidade, obtem-se uma aproxima¸caõ Xυ conforme a Equa¸caõ 3.20.

Xυ = ΦΦ T X

(3.20)

O erro nesta aproxima¸cão, então, medido pela norma-2 da diferen¸ca entre X e X υ é no máximo igual ao maior valor singular desprezado, conforme representado na Equa¸caõ 3.21.

∥X − X ∥ ≤ σ υ

υ+1

(3.21)

Pode-se mostrar que X υ é a melhor aproxima¸cão para X (no sentido da norma-2) entre todas as matrizes com posto r

≤ υ.

Se συ+1 é suficiente pequeno, manter somente υ

valores singulares resulta em uma boa aproxima¸caõ e, portanto, se diz que o posto efetivo de

X é υ. A redu¸caõ de dimensão ocorre através da proje¸caõ dos estados contidos na matriz

X de dimensão n em um subespa¸co de dimensão υ muito menor. Tal proje¸caõ é representada M i=1 .

{ }

a xi

M i=1 correspondentes

{ }

pela matriz Φ, chamada de base do POD, e gera estados reduzidos zi

Esta transforma¸caõ pode ser aplicada a qualquer estado x no espa¸co de dimensão

M para o cálculo de um estado reduzido z como mostrado na Equa¸caõ 3.22.

z = Φ T x

(3.22)

Por fim, para ilustrar o funcionamento do POD, considere o exemplo ilustrado na Figura 3.2. Suponha que os estados de uma sistema sejam representados em um espa¸co tridimensional dentro de um elipsóide como na Figura 3.2(a). As dire¸co˜es principais deste conjunto são os autovetores da matriz de covariância e estão representadas na Figura 3.2(b) por PC1, PC2 e PC3; o módulo dos vetores é proporcional aos autovalores e portanto a` variabilidade naquela dire¸caõ. Ao desprezar-se a dire¸caõ PC3 com menor autovalor, isto é, com menor variabilidade, os pontos passam a ser representados de maneira aproximada em um plano no interior de uma elipse como na Figura 3.2(c). Isto só é poss´ıvel pois a variabilidade na dire¸caõ principal desprezada é bastante pequena.


63

Fig. 3.2: Funcionamento do POD - Consideram-se somente as dire¸co˜es de maior variabilidade

3.3.2 Decomposi¸caõ em Valores Singulares Como descrito anteriormente, o POD necessita do cálculo dos autovalores e autovetores da matriz de covariância dos dados. Entretanto, só o cálculo dessa matriz, por si só, já seria custoso. Além disso, a determina¸caõ de autovetores é uma opera¸cão cara compu´ nesse sentido que se aplica a decomposi¸cão em valores singulares descrita tacionalemente. E nesta subse¸caõ.

× M de snapshots de simula¸caõ na qual será aplicada

Considere a matriz N

o POD. A sua decomposi¸cão em valores singulares (sigla em inglês, SVD – Singular Value Decomposition) consiste na sua fatora¸cão em um produto de três matrizes U, V e Σ conforme representado na Equa¸caõ 3.23.

X = U Σ VT

· ·

(3.23)

× N e M × M , respectivamente.

U e V são ortonormais e possuem dimensão N

× M , contendo os valores singulares não-negativos σ , j

Σ é uma matriz diagonal N

j

=

1, . . . , m i n(N, M ), dispostos em ordem não-crescente ao longo da diagonal. As colunas de U e V são denotadas pelos vetores u j , j = 1, . . . , n e v j , j = 1, . . . , M . Duas importantes expressões que relacionam U e V a Σ estão nas Equa¸co˜es 3.24 e 3.25 e podem ser deduzidas diretamente da Equa¸caõ 3.23.

X XT U = U Σ ΣT

·

·

· ·

(3.24)

3.4. Modelo de Ordem Reduzida - TPWL/POD

64

XT X V = V ΣT Σ

· ·

·

·

(3.25)

Das Equa¸co˜es 3.24 e 3.25 conclui-se que as colunas de U são autovalores de X XT

·

e as colunas de V são autovalores de XT X. Além disso, tanto os autovalores de X XT ,

·

·

encontrados na diagonal da matriz Σ ΣT , quanto os de X T X, na diagonal de Σ T Σ, são

·

·

·

iguais aos quadrados dos valores singulares de X. Portanto, estas expressões permitem que se utilize a SVD para o cálculo de autovalores e autovetores da matriz de covariância. ´ poss´ıvel demonstrar que o posto de X é igual ao número de valores singulares E diferentes de zero. Portanto, se rank(X) = r, é poss´ıvel reescrever a SVD em sua forma reduzida. Na Equa¸caõ 3.26 , U + e V + são as primeiras r colunas de U e V , respectivamente.

Σ+ é uma matriz diagonal r x r. Esta forma para a decomposi¸caõ reduz o espa¸co em memória necess´ ario na opera¸caõ. T X = U+ Σ+ V+

·

·

(3.26)

A computa¸cão das dire¸co˜es principais da matriz de covariância, que são a solu¸caõ do problema de minimiza¸cão da Equa¸cão 3.15 objetivo do POD, é feita utilizando-se a SVD descrita na Equa¸caõ 3.26. A matriz U + , de acordo com a Equa¸cão 3.24, é a matriz de todos os autovalores de X XT .

·

Sendo υ, υ < M << N , o n´ umero de dire¸co˜es principais utilizadas na aproxima¸caõ do POD, monta-se a matriz Uυ , considerando-se somente as υ primeiras colunas de U + . As colunas estão ordenadas em ordem decrescente dos autovalores associados e, portanto, correspondem aos maiores autovalores, representando os modos de maior energia. Desta forma, Uυ é a base do POD representada na Equa¸caõ 3.16. Isto está exposto na Equa¸caõ 3.27.

Uυ = Φ =

3.4

�

φ1 φ2

···

φυ

�

(3.27)

Modelo de Ordem Reduzida - TPWL/POD A lineariza¸caõ em torno de estados convergidos e gravados reduz a complexidade

do modelo, eliminando a necessidade de métodos iterativos para solu¸cão do sistema, transformando o problema de simula¸caõ em uma sequência de sistemas lineares. Neste trabalho,


65

a técnica POD (Proper Orthogonal Decomposition ) foi combinada ao TPWL (Trajectory Piecewise Linearization ) com o objetivo de reduzir a dimensão do problema a ser resolvido. O primeiro trabalho a aplicar este esquema à simula¸caõ de reservatórios foi (Cardoso, 2009). Faz-se necessário então gerar a base do POD, matriz que projeta vetores de um espa¸co de estados bastate grande (e. g., mapas de pressão e fra¸cões mássicas gerados em uma simula¸cão de reservatórios) em um espa¸co de dimensão reduzida que contém a maior parte da informa¸cão necessária. Para tanto, efetua-se uma ou mais simula¸co˜es do reservatório estudado, armazenando os estados (ou snapshots ) pelos quais a simula¸caõ passou a cada passo de tempo simulado. Entende-se por snapshot um par de vetores x p e xz contendo, respectivamente, os valores de pressã o e fra¸cã o mássica de todos as N c células do modelo. Considere então que X p e X z são matrizes que contém em suas colunas todos os snapshots calculados em todos os N passos de tempo, como mostra a Equa¸caõ 3.28.



X p = x p1 , x p2 ,



N p

··· , x



, Xz = x1z , x2z ,

· ·· , x



N z

(3.28)

Note que, o TPWL/POD já havia sido aplicado a` formula¸cão o´leo-água que considera as variáveis de estado: press˜ ao e satura¸cão (Cardoso, 2009; He, 2010). Portanto, o modelo de ordem reduzida aqui desenvolvido possui esta diferen¸ca. Uma aplica¸cão parecida foi realizada em (He and Durlofsky, 2013), que considerou um modelo composicional por fra¸co˜es molares. ´ importante ressaltar que, por causa da diferen¸ca de magnitude entre os valores E de pressão e fra¸co˜es mássicas, com o objetivo de evitar problemas num´ ericos no algoritmo, os estados são normalizados conforme é mostrado na Equa¸cão 3.29. x pi =

− −

− −

( pi pmin ) (z i z min ) , xzi = ( pmax pmin ) (z max z min )

(3.29)

A base do POD utilizada para redu¸cão da dimensão dos mapas de pressã o e satura¸caõ é a matriz U da decomposi¸cão em valores singulares da matriz de snapshots . São realizadas duas decomposi¸co˜es em separado com os mapas normalizados de pressões e os de fra¸co˜es mássicas, obtendo-se as matrizes Φ p e Φz . Isto permite a escolha de um número de autovalores diferentes para as grandezas que têm naturezas bastante distintas. Ao combinaremse essas duas matrizes, obtem-se a base Φ que será usada nos cálculos de redu¸caõ da dimensão do TPWL como proje¸cão dos estados no espa¸co reduzido representada na Equa¸caõ 3.30.


66

·

x = Φ z

(3.30)

Em suma, Φ p foi determinada a partir da decomposi¸caõ da matriz de snapshots de pressão X p e Φ z e foi calculada a partir da decomposi¸caõ da matriz de snapshots de fra¸cões mássicas X z . Os snapshots de dimensão reduzida têm duas partes também. A maneira como as bases e snapshots foram combinados está descrita na Equa¸caõ 3.31.

    

     ≈     

x p1 xz1 x p2 xz2 .. . x pN c xzN c

φ p11 0 φ p12 0 .. . φ p1N c 0

·· · ·· · ·· · ·· ·

φ p1p

υ

0

0

φ1z1

φ p2p

υ

0

0 .. .

φ1z2 .. .

·· · ·· ·

p φ pN c

υ

0

0

φ1zN c

...

··· ··· ··· ···

...

··· ···

0 φυz1z 0 φυz2z .. . 0 z φυzN c

    ·    

   

z p1 .. . z pυp z z1 .. . z zυz

(3.31)

Finalmente, a redu¸caõ de dimensã o pode ser aplicada ao TPWL. A Equa¸caõ 3.32 mostra como fica a expressão do TPWL que era representada pela Equa¸cão 3.9 após a aplica¸cão do POD para redu¸caõ da dimensão através da pro je¸cão da Equa¸cão 3.30. Note que os snapshots x de dimensão completa se transformam em z de dimensão reduzida através da transforma¸caõ representada na Equa¸cão 3.30.

n+1

z

= z

(−  ��

i+1

i+1

J

1

−

r

∂ Ai+1 ∂ xi

n

(z r

i

−z )+

� � ∂ Qi+1 ∂ ui+1

n+1

(u r

i+1

−u

)

�

(3.32)

Por outro lado, há a necessidade de calcular também as derivadas com dimensão reduzida (Ji+1 )r ,

∂ Ai+1 ∂ xi

r

,

∂ Qi+1 ∂ ui+1

r

. Coerentemente com a Equa¸caõ 3.30, faz-se necessário

pós-multiplicar todas as derivadas em rela¸caõ aos estados pela base do POD. Entretanto, isto faz com que o sistema fique com mais equa¸cões do que incognitas. Existem diferentes estrat´ egias para resolver este sistema. A primeira maneira, já utilizada em um contexto de elementos finitos, é conhecida como proje¸caõ de BubnovGalerkin e representa uma proje¸caõ ortogonal no espa¸co de estados reduzido. Foi aplicada


67

a este problema de constru¸caõ de modelo de ordem reduzida em simula¸cão de reservatórios em (Cardoso, 2009; He, 2010). Consiste em pré-multiplicar as derivadas pela transposta da base, Φ T , e, portanto, é uma proje¸caõ ortogonal no espa¸co reduzido. Entretanto, problemas de estabilidade do método foram relatados em (He, 2010; Gildin et al., 2013) o que levou a busca de outras formas. A proje¸cão escolhida neste trabalho é a de Petrov-Galerkin. Foi relatada como respons´ avel pela estabiliza¸caõ em (Carlberg et al., 2009; He and Durlofsky, 2013) e é equivalente a resolver o sistema por m´ınimos quadrados, j´ a que consiste em pré-multiplicar as derivadas T ´ uma proje¸caõ obl´ıqua no espa¸co de dimensão reduzida. A Equa¸caõ 3.33 por ΦT (Ji+1 ) . E mostra todas as três derivadas com dimensões reduzidas a serem utilizadas em 3.32.

  � �  � � 

T

(Ji+1 )r = Φ T (Ji+1 ) Ji+1 Φ ∂ Ai+1 ∂ xi

∂ Qi+1 ∂ ui+1

T

r

i+1 T

= Φ (J

)

T

r

= Φ T (Ji+1 )

� � � � ∂ Ai+1 ∂ xi

Φ

(3.33)

∂ Qi+1 ∂ ui+1

Por fim, a constru¸caõ do modelo de ordem reduzida é o que se chama aqui de treinamento e precisa calcular todos os objetos utilizadas na solu¸caõ sequencial da Equa¸caõ 3.32. O processo de treinamento do TPWL/POD engloba a execu¸caõ de uma ou mais simula¸cões do SIMPAR armazenando os mapas e derivadas convergidos, seguida do cálculo da base e dos mapas reduzidos e finalmente o cálculo das derivadas reduzidas. Neste trabalho o treinamento foi realizado com uma única simula¸caõ do SIMPAR com controles de pressão de fundo variando aleatoriamente entre os limites máximo e m´ınimo estipulados. Estes controles variam com uma frequência maior do que será utilizado na simula¸cão posterior para estimular a variabilidade dos estados pelos quais a trajetória de treinamento passa. A quantidade de autovalores (υ p e υz ) considerados na redu¸caõ da dimens˜ ao é definida empiricamente através de um simples estudo da qualidade da aproxima¸caõ. O algoritmo de treinamento da maneira como foi implementado está representado abaixo: 1. Seleciona-se uma sequência aleatória de controles de pressão de fundo a serem aplicadas na simula¸caõ de treinamento. 2. Realiza-se a simula¸caõ de treinamento. Armazenam-se os mapas, derivadas, controles utilizados, informa¸co˜es de po¸cos, etc. em um arquivo binário.


68

3. Leem-se os mapas e montam-se as matrizes normalizadas de snapshots de pressão X p e fra¸cões mássicas X z . 4. Calculam-se as bases Φ p e Φ z através da SVD considerando as quantidades de autovalores υ p e υz respectivamente. Combinam-se as bases, obtendo-se Φ. 5. Calculam-se os mapas reduzidos z i através da Equa¸caõ 3.30, mantendo-os na memória para uso futuro. 6. Leem-se as derivadas (J

i+1

treinamento.

),

� �� ∂ Ai+1 ∂ xi

,

∂ Qi+1 ∂ ui+1

7. Calculam-se as derivadas reduzidas (Ji+1 )r , mantendo-as na memória.

exportadas durante a simula¸ca˜ o de

∂ Ai+1 ∂ xi

r

,

∂ Qi+1 ∂ ui+1

r

conforme Equa¸caõ 3.33,

8. Realiza-se a simula¸caõ TPWL/POD utilizando-se a Equa¸caõ 3.32 e o algoritmo de simula¸caõ TPWL descrito na se¸caõ 3.1 com os controles utilizados para o treinamento com o objetivo de avaliar a qualidade da aproxima¸caõ. Tamb´ em foi desenvolvida uma metodologia para retreinamento. Isto consiste na adi¸caõ dos estados e derivadas calculados em uma nova simula¸caõ ao conjunto de estados considerados no treinamento inicial. A base precisa ser recalculada considerando os estados de todas as simula¸cões. A justificativa para o retreinamento é extrapolar a capacidade de representa¸cão do TPWL/POD em uma trajetória próxima a um conjunto de controles espec´ıfico e será u ´ til durante a otimiza¸cão na medida em que o processo avance e a qualidade da aproxima¸caõ diminua. No cap´ıtulo 4 ser´ a desenvolvido um crit´ erio para iniciar o retreinamento. O algoritmo de retreinamento pressupõe que um conjunto de simula¸cõ es com o SIMPAR já foi realizado e que uma delas teve os mapas e derivadas exportadas para serem considerados na reconstru¸caõ do TPWL/POD. Pode considerar uma ou mais simula¸co˜es no conjunto para retreinamento, dependendo da configura¸cã o escolhida. Note que os mapas e derivadas reduzidos de todas as simula¸co˜es ficam na memória, dispon´ıveis para serem utilizados e cuidados devem ser tomados para não perder a sua correspondência temporal. Os passos do retreinamento, após a simula¸cão do SIMPAR, são: 1. Realiza-se a leitura dos mapas de pressã o e fra¸co˜es mássicas da primeira simula¸caõ de treinamento e de todas aquelas simula¸cões que se deseja utilizar no retreinamento. Montam-se as matrizes normalizadas de snapshots de pressão X p e fra¸co˜es mássicas X z

3.5. Detalhes de Implementa¸caõ

69

com os mapas de todas as simula¸co˜es concatenados. 2. Executam-se os passos 4 a 7 do algoritmo de treinamento obtendo-se os mapas redui

i+1

zidos z e as derivadas reduzidas ( J

)r ,

(treinamento e retreinamento).

� �� ∂ Ai+1 ∂ xi

r

,

∂ Qi+1 ∂ ui+1

r

de todas as simula¸cões

3. Realizam-se a simula¸co˜es TPWL/POD utilizando-se a Equa¸caõ 3.32 e o algoritmo descrito na se¸caõ 3.1 com os controles utilizados para o retreinamento e outras sequências de controles que já tenham sido simuladas com o objetivo de corrigir o resultado. Foi desenvolvida uma maneira de medir a qualidade da aproxima¸caõ realizada pelo modelo de ordem reduzida baseada nos estados calculados na simula¸cão. Isto é necessário para o algoritmo de otimiza¸caõ estudado no Cap´ıtulo 4 como um mecanismo para avaliar a necessidade de retreinamento. Na equa¸caõ 3.32, observa-se que o modelo linear que constitui o TPWL/POD é baseado na distˆ ancia entre os estados e controles gravados e calculados durante a simula¸caõ. O estado reduzido gravado zi a ser utilizado no passo de tempo n é escolhido de forma a minimizar a distância ao u ´ ltimo estado zn gravado na simula¸cão pelo TPWL/POD. Esta distância considerada diz respeito somente à parte do vetor de estados correspondente às fra¸cõ es mássicas znz

i z

∥ − z ∥ =

√ − (znz

ziz )T (znz

i z

· − z ).

Sendo assim, a

medida de erro na estimativa de estados realizada pelo modelo E tpwl foi definida como o valor médio normalizado desta distância ao longo do tempo de simula¸caõ T , conforme equa¸caõ 3.34.

E tpwl =

3.5

N

1

� ( ∥ − × √

T

υz

znz

n=1

ziz

∥ × ∆t

n



(3.34)

Detalhes de Implementa¸ c˜ ao A implementa¸caõ do algoritmo TPWL foi totalmente realizada utilizando-se o

Matlab. Entretanto, foi necessário alterar o código do SIMPAR de forma a proporcionar a exporta¸caõ dos dados necessários para a constru¸cão do modelo de ordem reduzida. Este simulador é em sua maior parte desenvolvido em FORTRAN 77, contendo algumas rotinas e módulos em FORTRAN 90. Desta forma, foram desenvolvidas algumas rotinas dentro do simulador que acessassem as variáveis e matrizes a serem exportadas e escrevessem em um arquivo de interface. Isto faz com que o TPWL seja considerado uma técnica semi-intrusiva, uma vez que não


70

se precisa interferir nos cálculos do simulador, mas necessita acessar variáveis internas que normalmente não estariam dispon´ıveis. A u ńica interferência nos cálculos do simulador diz respeito ao critério de convergência do método de Newton-Raphson. Há dois crit´ erios para que o método pare de

∥ ∥ < ϵ) e o segundo

iterar: o primeiro pelo valor do res´ıduo, que deve tender a zero ( g

∥ ∥

por varia¸caõ do vetor de estados que também tende a zero ∆x < ϵ. Entretanto, caso a convergência ocorra devido à varia¸cão dos estados, o simulador não tem a necessidade de

atualizar as derivadas para o cálculo do res´ıduo. Sendo assim, os estados ficam “defasados” em rela¸cão a`s derivadas, o que gera uma incoerência e traz problemas ao funcionamento do TPWL. Portanto, há duas formas de resolver o problema. A mais robusta e mais intrusiva seria alterar o loop do método de Newton de modo a obrigar a atualiza¸caõ das derivadas em qualquer caso. Entretanto, a forma escolhida neste trabalho, consiste em desabilitar a convergência por varia¸caõ dos estados, só permitindo que o processo iterativo pare com um res´ıduo suficientemente pequeno e com as derivadas calculadas coerentemente com os mapas de estados. Outro aspecto importante que refor¸ca o aspecto semi-intrusivo do metamodelo ´ necessário conhecer precisamente os detalhes de está relacionado ao modelo de po¸co. E implementa¸cão do modelo de po¸co utilizado pelo simulador. O TPWL só é capaz de estimar as variáveis de estado (pressã o e fra¸cões mássicas), mas o que de fato se está interessado em calcular são as vazões de fluido produzido/injetado. Este cálculo é realizado a partir das vari´ aveis de estado nas células dos po¸cos pelo modelo de po¸co do simulador, que deve estar inteiramente reproduzido dentro do metamodelo. Ainda tratando dos po¸cos, o TPWL/POD precisa fazer o cálculo das variáveis de estado nas células dos po¸cos para que o modelo de po¸co possa ser aplicado e permitindo o cálculo de vazão. Entretanto, o restante é inteiramente computado no espa¸co reduzido, onde o vetor de estado não tem significado f´ısico direto. Desta forma, utilizam-se somente as linhas correspondentes às células dos po¸cos na matriz Φ para reverter-se a redu¸cão de dimensão. Isto reduz o produto de matrizes a ser efetuado e o custo computacional. Um problema relacionado à aproxima¸caõ realizada pelo TPWL, e que ainda não foi abordado na literatura, é o fato de ele poder calcular valores não-f´ısicos para os estados como, por exemplo, pressões menores que zero e fra¸cões mássicas fora da faixa entre zero e um. No u ´ ltimo caso, propaga-se o problema para as satura¸co˜es obtidas pela Equa¸caõ 2.18, o que leva a um erro no cálculo da permeabilidade relativa e mobilidade de cada fluido no


71

po¸co (ver Equa¸caõ 2.34) e a uma estima¸caõ errônea da vazão produzida/injetada. Este problema foi de certa forma evitado, impondo-se limites f´ısicos às fra¸cões mássicas nos po¸cos antes dos cálculos do modelo de po¸co. Entretanto, é poss´ıvel que as vari´ aveis de estado calculadas pelo metamodelo possam assumir valores não-f´ısicos em muitas células e estes valores sejam propagados ao longo do tempo. Como o TPWL é inteiramente computado no espa¸co reduzido, a única forma de avaliar a viabilidade das variáveis de estado e aplicar os limites f´ısicos seria reverter a redu¸caõ de dimensão em todas as células, o que poderia tornar o cálculo bastante caro. Testes estão sendo realizados nesse sentido. Um ponto cr´ıtico desta implementa¸caõ está exatamente na troca de informa¸cões entre o simulador e a ferramenta desenvolvida no Matlab. De fato, esta exporta¸caõ das derivadas e mapas precisa de uma boa quantidade de espa¸co em disco e consome tempo de escrita/leitura dos arquivos. Portanto, a escolha do formato deste arquivo tem grande influência no desempenho da ferramenta. Inicialmente, o formato escolhido para a troca de informa¸cões era basicamente um conjunto de arquivos-texto. Entretanto, o volume de dados escrito era bastante grande e, por causa da lentidão da escrita, o tempo de simula¸caõ era multiplicado por três. Além disso, a leitura destes dados fazia com que o tempo para treinamento se comparasse ao de duas simula¸co˜es com o SIMPAR. Com este tempo de treinamento de cinco a sete vezes o de uma simula¸caõ, tornar-se-ia proibitivo realizar um retreinamento ao longo de um processo de otimiza¸caõ caso fosse necessário. Sendo assim, iniciou-se a busca de um formato de arquivo que fosse mais eficiente. O ideal é que tal formato fosse binário, para economizar espa¸co e reduzir o tempo de acesso através de acesso direto. Entretanto, o formato binário do FORTRAN padrão não se mostrou compat´ıvel com o Matlab, o que demandaria a constru¸caõ de rotinas complexas para leitura e a concep¸caõ de um padrão de ponteiros e organiza¸caõ do arquivo. Decidiu-se então utilizar o formato HDF - Hierarchical Data Format , que está em sua quinta versão: HDF5. (HDFGroup, 2014) Tal formato é binário e genérico, onde as informa¸co˜es podem ser organizadas de maneira hierárquica, como numa árvore de diretórios, mas são escritas de forma sequencial. Para acessá-las existe um sistema de ponteiros que permitem a leitura de maneira direta. Outra importante vantagem reside no fato de tudo isto ser totalmente transparente ao usuário. Existe uma biblioteca para diversas linguagens de programa¸cão como C, C++ e Fortran, o que permite utilizá-lo de maneira bastante simples. Além disso, o Matlab já possui fun¸co˜es relativamente simples para leitura. Ao optar-se por


72

esta forma de exporta¸caõ reduziu-se significamente o tempo de escrita e de leitura, bem como o volume de dados gerados. Todas as informa¸cões exportadas, vazões, mapas, derivadas, controles, etc, são encapsuladas em um arquivo HDF que contém uma estrutura dividida por categorias. Sendo assim, deixou-se de exportar um conjunto de arquivos-texto para exportar-se um único arquivo binário com toda a informa¸cão necessária. Com o uso do HDF o tempo de treinamento, considerando a simula¸caõ com exporta¸cão de dados, somado ao tempo de leitura dos mapas e derivadas, constru¸caõ da base e cálculo das derivadas reduzidas, equivale ao tempo de uma e meia simula¸cão simples. Outro aspecto importante da implementa¸caõ é que o processo de treinamento foi dividido em duas etapas. A primeira parte lˆ e todos os mapas e realiza o c´ alculo da base utilizando SVD e o cálculo dos mapas reduzidos. A segunda parte faz a leitura das derivadas e controles utilizados e efetua a redu¸caõ da dimensão, efetuando o produto da base pelas matrizes de cada passo de tempo. Estas rotinas foram separadas para permitir que a base fosse calculada a partir de uma simula¸caõ, e as derivadas, caso se fizesse necessário, a partir de uma simula¸caõ diferente. A terceira rotina implementada no MATLAB, e que integra a base do modelo substituto, é responsável pela simula¸caõ, isto é, a sequência de solu¸co˜es de sistemas lineares que representa a simplifica¸cão do problema de simula¸caõ. Esta é a rotina que deve ser chamada ao longo do processo de otimiza¸caõ como substituta do SIMPAR. Vale ressaltar que inicialmente os mapas reduzidos, as respectivas derivadas reduzidas, os controles e todos os outros dados necessários eram mantidos em arquivos binários do MATLAB. Entretanto, isto se mostrou ineficiente, uma vez que para cada simula¸cão do TPWL era necessário a leitura de arquivos. Isto posto, uma vez que a rotina de otimiza¸cão que chamaria o modelo substituto ser´ a também desenvolvida em MATLAB e executada pela mesma instância, optou-se por manter todas aquelas variáveis na memória, como variáveis globais, pagando-se o pre¸co de ocupar uma quantidade significativa de memória RAM. Também foram constru´ıdas rotinas para encapsular a execu¸caõ do TPWL, bem como a execu¸caõ do simulador. Na verdade, o papel dela é permitir que o otimizador em processo de desenvolvimento seja capaz de executar o SIMPAR e o TPWL com controles normalizados entre zero e um, onde zero significa o valor m´ınimo da pressão de fundo e um, o valor máximo.

3.6. Compara¸caõ com o modelo de alta fidelidade

73

Por fim, foi desenvolvida uma rotina para retreinamento. Esta basicamente inclui entre os mapas e derivadas reduzidas já calculados um novo conjunto obtido com outra simula¸cão. Para tanto, precisa recalcular a base considerando todos os mapas. Al´ em disso, como a base é recalculada, é necessário recalcular os mapas reduzidos da simula¸caõ que já haviam sido calculados no treinamento.

3.6

Compara¸ca õ com o modelo de alta fidelidade Com o objetivo de avaliar a qualidade da aproxima¸caõ realizada pelo modelo de

ordem reduzida aqui estudado, foram executados alguns estudos comparativos. O u ńico modelo de reservatórios aqui considerado é o SPE10 - Modificado descrito no cap´ıtulo 2, considerando três varia¸co˜es. A primeira, mais simples, considera que as compressibilidades do óleo, água e rocha são nulas e que as densidades dos fluidos são idênticas. A segunda também é incompress´ıvel, mas considera uma diferen¸ca de densidades entre água e óleo de 5 lbm/ft3 . A terceira é um modelo compress´ıvel com a mesma diferen¸ca de densidades entre fluidos; considera compressibilidades do óleo de co = 3 da rocha, cr = 1

× 10

7

−

× 10

6

−

psi 1 , da água cw = 1 −

× 10

6

−

psi

1

−

e

psi 1 . A primeira varia¸caõ se deve ao fato de o problema ser menos −

não-linear e portanto mais reprodut´ıvel pelo modelo de ordem reduzido aqui considerado. ´ importante ressaltar que na literatura ainda é bastante comum considerar-se os casos E incompress´ıveis para avaliar o uso dessa técnica (Cardoso, 2009; He, 2010; He et al., 2013). A compara¸cão de resultados é feita em termos do VPL (Valor Presente L´ıquido) que representa a fun¸caõ objetivo a ser otimizada. O pre¸co do óleo é 100 $/BBL e os custos de inje¸cão e produ¸caõ de água são de 10 $/BBL e a taxa de desconto é de 10 % a. a. Há quatro po¸cos produtores e dois injetores instalados, todos abrindo no in´ıcio da simula¸caõ. A pressão de inje¸cão é fixada em 12000 psi no fundo do po¸co. A pressão de fundo em fluxo (BHP) nos po¸cos produtores varia entre 2000 psi e 6000 psi. O per´ıodo de simula¸caõ aqui considerado é de 3000 dias. O reservatório (mapa de permeabilidade) está representado na Figura 2.2. A sequência de controles para o primeiro treinamento foi determinada aleatoriamente através de um hipercubo latino considerando-se os controles dos quatro po¸c os a cada 200 dias como variáveis. Portanto há 15 ciclos de controles. As sequências utilizadas em todos os po¸cos estão representadas na Figura 3.3. O que se espera com isto é uma maior variabilidade dos estados gravados com o incremento da capacidade de aproxima¸caõ do TPWL/POD. O processo de retreinamento que será aplicado a este teste considera os


74

estados gravados inicialmente adicionando os estados de mais uma simula¸ cão ao conjunto. 6000 5500 ) i 5000 s p ( o d4500 n u F 4000 e d o ã3500 s s e r P3000

2500 2000 0 W1

500 W2

1000 W3

1500 W4

2000 2500 3000 Tempo (Dias)

Fig. 3.3: Sequência de controles de treinamento (pressões de fundo) aplicados aos po¸cos produtores O número de autovalores considerados na redu¸caõ de dimensã o do POD foi de ´ importante ressaltar que υ p = 128, para a pressão e υz = 256, para as fra¸co˜es mássicas. E a dura¸cã o máxima dos passos de tempo neste problema foi de 30 dias, mas o simulador decidiu cortá-los bastante o que gerou estados suficientes para utilizar tamanha quantidade autovalores. A escolha do n´ umero de autovalores foi feita após um conjunto rápido de testes que mostrou que estas quantidades eram suficientes para uma boa aproxima¸caõ. Os u ´ltimos valores singulares considerados estão entre a ordem de 10

3

−

e 10 4 . A figura 3.4 apresenta −

um gráfico dos valores singulares para a SVD da matriz de snapshots de pressões (3.4(a)) e de fra¸co˜es mássicas (3.4(b)) na simula¸caõ de treinamento dispostos em ordem decrescente com a indica¸caõ dos respectivos pontos de corte escolhidos. Portanto, em resumo, o teste aqui realizado para avalia¸caõ da acurácia do modelo reduzido pode ser descrito pelos passos a seguir. 1. Escolhem-se aleatoriamente um conjunto de 15 controles entre 2000 psi e 6000 psi para cada po¸co produtor utilizando-se hipercubo latino. Estes controles serão aplicados aos po¸cos, mudando-se o valor a cada 200 dias, totalizando 3000 dias de simula¸caõ. 2. Utilizando-se o SIMPAR, realiza-se a simula¸cão do modelo aplicando-se os controles escolhidos.


5

5

10

10

Incompressível − Dens. Iguais Incompressível − Dens. Dif. Compressível − Dens. Dif.

0

0

10

10

−5

−5

10

10

−10

−10

10

75

0

100

200

300

400

(a) Pressão

500

600

700

10

0

100

200

300

400

500

600

700

(b) Fra¸cão mássica

Fig. 3.4: Valores singulares dispostos em ordem decrescente - Modelo SPE10-Modificado 3. Executam-se em sequência as rotinas de cálculo da base e mapas reduzidos (POD) e de cálculo das derivadas reduzidas, realizando-se assim o treinamento do modelo de ordem reduzida. 4. Executam-se algumas simula¸cões com o SIMPAR com dura¸cão de 3000 dias, considerando-se que em cada haverá 3 ciclos de controle de 1000 dias. Anotam-se para cada uma delas o valor de VPL calculado. Exportam-se todos os dados de mapas e derivadas, para que se possa vir a fazer retreinamento do TPWL no futuro. Considera-se que os controles de pressão aplicados em todos os po¸cos são iguais, permitindo-se a varia¸caõ dos controles de um ciclo para outro. São três os valores que serão avaliados em cada ciclo: 2000 psi (limite inferior), 4000 psi (meio da escala) e 6000 psi (limite superior). Como os controles são iguais para todos os po¸cos, tem-se três variáveis, uma para cada ciclo. Como são três ciclos, com a possibilidade de três valores para cada, serão assim 27 simula¸co˜es no total. 5. Executam-se as mesmas simula¸cões com o TPWL treinado no passo 3, anotando-se o seu resultado de VPL. 6. Efetua-se o retreinamento do TPWL, adicionando-se os mapas e derivadas de uma simula¸caõ escolhida entre todas as calculadas no passo 4. Utiliza-se para isso o algoritmo de retreinamento descrito na se¸caõ 3.4, que sempre mantém as informa¸co˜es da simula¸caõ realizada no passo 2 no conjunto de dados do modelo de ordem reduzida. Adicionam-se as informa¸co˜es de uma das simula¸cões efetuadas em 4 eliminando-se as informa¸cões do


76

u ´ltimo retreinamento. 7. Executam-se as mesmas 27 simula¸cões com o TPWL retreinado, anotando o VPL. 8. Retorna-se ao passo 6, de modo a retreinar o TPWL com uma outra simula¸caõ. Repetese este ciclo até esgotar as simula¸co˜es escolhidas para retreinamento. A tabela 3.1 contém as 27 sequências dos controles aplicados igualmente aos quatro po¸cos nos trˆ es ciclos de controle. Nesta tabela, a palvra MIN, siginifica valor m´ınimo do controle, 2000 psi, a palavra MED representa valor médio, 4000 psi e MAX, o máximo, 6000 psi. As simula¸co˜es escolhidas para retreinamento do TPWL no passo 6 do teste são as obtidas com as sequências de controles de número 1, 5, 10, 20, 24 e 27. Estas foram escolhidas de forma a que haja a mesma quantidade de simula¸cões de retreinamento com os três n´ıveis de controle no primeiro ciclo. Tal escolha se deve à forte dependência que se espera haver entre o VPL e o valor inicial dos controles devido a taxa de desconto. Sendo assim, entre as simula¸cões utilizadas para retreinamento há duas com valor máximo de pressão de fundo no primeiro ciclo, duas com valor médio e duas com valor m´ınimo. Para avaliar a qualidade da aproxima¸caõ de maneira rápida e visualmente simples, foram constru´ıdos crossplots nos quais o eixo x representa o VPL calculado pelo SIMPAR, que seria a medida real, e o eixo y representa a aproxima¸cão pelo modelo de ordem reduzida. Sendo assim, se o TPWL fosse perfeito, todos os pontos cairiam sobre a reta identidade em azul. H´ a um conjunto de crossplots para o caso incompress´ıvel com densidades iguais, um para o caso incompress´ıvel com densidades diferentes e um para o compress´ıvel. Em cada conjunto foi feito um gráfico para o primeiro treinamento e um para cada retreinamento considerado no teste, como mostram as Figuras 3.5, 3.6 e 3.7. Na Figura 3.5 é considerado o modelo SPE10 - modificado sem compressibilidades e densidades iguais para á gua e óleo. Em cada crossplot o eixo x representa o VPL obtido na simula¸cão do SIMPAR e o eixo y o VPL calculado com o TPWL. As nuvens representam os seguintes casos: (a) sem retreinamento, (b) com retreinamento na simula¸caõ 1, (c) com retreinamento na simula¸caõ 5, (d) com retreinamento na simula¸cão 10, (e) com retreinamento na simula¸caõ 20, (f) com retreinamento na simula¸cão 24 e (g) com retreinamento na simula¸caõ 27. Neste caso, obtem-se boa aderência entre simulador e modelo de ordem reduzida devido ao baixo n´ıvel de não linearidades proporcionado pela ausência de diferen¸cas em densidades dos fluidos.


77

1o Ciclo 2o Ciclo 3o Ciclo 1

MIN

MIN

MIN

2

MED

MIN

MIN

3

MAX

MIN

MIN

4

MIN

MED

MIN

5

MED

MED

MIN

6

MAX

MED

MIN

7

MIN

MAX

MIN

8

MED

MAX

MIN

9

MAX

MAX

MIN

10

MIN

MIN

MED

11

MED

MIN

MED

12

MAX

MIN

MED

13

MIN

MED

MED

14

MED

MED

MED

15

MAX

MED

MED

16

MIN

MAX

MED

17

MED

MAX

MED

18

MAX

MAX

MED

19

MIN

MIN

MAX

20

MED

MIN

MAX

21

MAX

MIN

MAX

22

MIN

MED

MAX

23

MED

MED

MAX

24

MAX

MED

MAX

25

MIN

MAX

MAX

26

MED

MAX

MAX

27

MAX

MAX

MAX

encias Tab. 3.1: Sequˆ

de controles utilizadas no estudo

3.6. 3.6. Compar Compara¸ a¸c˜ caõ com o mo delo de alta fidelidade

78

8

8

8

7.5

7.5

7.5

7

7

7

6.5

6.5

6.5

6

6

6

5.5

5.5

5.5

5

5

5

4.5

4.5

4.5

4 4

5

6

7

8

(a) Somente Treinamento

4 4

5

6

7

8

(b) Retreinamento - Sim. 01

4 4

8

8

8

7.5

7.5

7

7

7

6.5

6.5

6.5

6

6

6

5.5

5.5

5.5

5

5

5

4.5

4.5

4.5

5

6

7

8

(d) Retreinamento - Sim. 10

4 4

5

6

7

6

7

8

(c) Retreinamento - Sim. 05

7.5

4 4

5

8

(e) Retreinamento - Sim. 20

4 4

5

6

7

8

(f) Retreinamento - Sim. 24

8 7.5 7 6.5 6 5.5 5 4.5 4 4

5

6

7

8

(g) Retreinamento - Sim. 27

cao aõ do VPL calculado pelo SIMPAR e pelo TPWL/POD para o modelo Fig. 3.5: Compara¸c˜ SPE10 - Modificado sem compressibilidades e densidades iguais. Já na Figura 3.6 é considerado o modelo mo delo SPE10 - modificado mo dificado sem compressibilidades e densidades diferentes para água agua e óleo. oleo. A ordem dos casos de retreinamento retreinam ento é idêntica entica a` Figura Figura 3.5. Nota-se Nota-se uma piora significativ significativaa na qualidade da aproxima¸ aproxima¸ c˜ cao aõ com grandes discrepâncias ancias devido ao aumento da não-linearid ao-linearidade ade por p or causa da diferen¸ diferen¸ca ca entre as densidades dos fluidos. Finalmente, a Figura 3.7 mostra o resultado para o modelo SPE10 - modificado com compressibilidades e densidades diferentes para água a gua e óleo o leo.. A orde ordem m dos dos caso casoss de retrei ret reina namen mento to é idêntica enti ca às as Figuras 3.5 e 3.6. Ocorrem Ocorrem resultados resultados ainda mais discrepan discrepantes tes entre TPWL e SIMPAR devido à atribui¸c˜ cao aõ de compressibili compres sibilidade dade acarretand ac arretandoo a dependˆ d ependência encia


79

8

8

8

7.5

7.5

7.5

7

7

7

6.5

6.5

6.5

6

6

6

5.5

5.5

5.5

5

5

5

4.5

4.5

4.5

4 4

5

6

7

8


4 4

5

6

7

8


4 4

8

8

8

7.5

7.5

7

7

7

6.5

6.5

6.5

6

6

6

5.5

5.5

5.5

5

5

5

4.5

4.5

4.5

5

6

7

8


4 4

5

6

7


6

7

8


7.5

4 4

5

8

4 4

5

6

7

8


8 7.5 7 6.5 6 5.5 5 4.5 4 4

5

6

7

8


cao aõ do VPL calculado pelo SIMPAR e pelo TPWL/POD para o modelo Fig. 3.6: Compara¸c˜ SPE10 - Modificado sem compressibilidades e densidades diferentes. temporal temporal dos resultados. resultados. Em todos os três es casos estudados observou-se uma dependência encia muito forte da qualidade dos resultados do TPWL com a escolha do primeiro controle da simula¸c˜ cao aõ utilizada para retreinamento. Nos crossplots Nos crossplots apresentados apresentados nas Figuras 3.5, 3.6 e 3.7, os pontos verdes representam simula¸c˜ cões oes cujo primeir pri meiroo controle contr ole é o valor máximo axim o para a press˜ pres são ao de fundo, os pontos po ntos preto p retoss têm em o prim p rimeiro eiro controle contr ole no n o valor valo r m´ m´ınimo ıni mo e os vermelh ver melhos os o têm em primei pr imeiro ro control co ntrolee no valor intermediário. ario. O pon ponto to marcado marcado com um diaman diamante te rosa repres represen enta ta a simula simula¸c˜ çao aõ utilizada para retreinamento e por isso recai sobre a reta, já que seu resultado passa a ser praticamente exato.


80

8

8

8

7.5

7.5

7.5

7

7

7

6.5

6.5

6.5

6

6

6

5.5

5.5

5.5

5

5

5

4.5

4.5

4.5

4 4

5

6

7

8


4 4

5

6

7

8


4 4

8

8

8

7.5

7.5

7

7

7

6.5

6.5

6.5

6

6

6

5.5

5.5

5.5

5

5

5

4.5

4.5

4.5

5

6

7

8


4 4

5

6

7


6

7

8


7.5

4 4

5

8

4 4

5

6

7

8


8 7.5 7 6.5 6 5.5 5 4.5 4 4

5

6

7

8


cao aõ do VPL calculado pelo SIMPAR e pelo TPWL/POD para o modelo Fig. 3.7: Compara¸c˜ SPE10 - Modificado com compressibilidades e densidades diferentes. Observa-se Observa-se uma certa estratifica¸c˜ cão, ao, onde o valor utilizado para o primeiro controle divide divide a nuvem nuvem em trˆ es es grupos de ponto p ontoss quase colineares. colineares. Isto se deve deve à grande for¸ca c a do primeiro ciclo de controle na forma¸c˜ cao aõ do VPL por causa da taxa de desconto, que maximiza a influência encia de erros no primeiro ciclo no erro total. Entretanto, Entretanto, o mais interessante interessante é que o retreinamento se mostra menos efetivo quando a simula¸c˜ cão ao utilizada possui o primeiro ciclo controlado pelo valor de pressão ao intermedi´ ario ario (pontos (pontos pretos) pretos).. Uma explica explica¸¸c˜ cao aõ para este fato é que a sequência encia de treinamento treiname nto utilizada utiliza da tem média edia próxima oxima do valor intermediário ario e, portanto, a adi¸c˜ cão ao de uma simula¸ simu la¸c˜ cao aõ de retreinamento com primeiro ciclo controlado no valor valo r médio edi o não ao agrega muita informa¸c˜ cao aõ para o modelo TPWL/POD.


81

Outra forma escolhida para medir a qualidade da aproxima¸c˜ cao aõ é atrav at rav´és es do erro er ro médio edi o quadr´ qua drático atico considerando-se os resultados de todas as simula¸c˜ coes. o˜ es. Para ara o caso caso inincompress´ compress´ıvel e o compress´ compress´ıvel foram tabelados os valores valores de VPL obtidos com a simula¸c˜ cao aõ do SIMPAR SIMPAR e do TPWL sem retreinam retreinament entoo e com retreinamen retreinamento. to. Considerand Considerandoo V P LREF e V P LROM respctivamente o valor calculado pelo SIMPAR SIMPAR (referˆ (referência) encia) e o calculado pelo modelo de ordem reduzida, calculam-se os desvios médios edios quadráticos aticos pela Equa¸c˜ cão ao 3.35 levando-se em conta as N as N = = 27 simula¸c˜ coes o˜es em cada um dos casos considerados.

  ∑ �  N

RMSD = RMSD =

i=1

V P LROM −V P LREF V P LREF

N

�

2

∗ 100%

(3.35)

Os valores obtidos para o RMSD estão ao compilados compilados na Tabela 3.2, de onde é poss´ poss´ıvel observar que os casos com densidades diferentes, devido a sua não ao linearidade maior apresent apresentam am desvios maiores. Note que a presen¸ presen¸ca ca de compress c ompressibilida ibilidade de tamb´ t ambém em aumenta a umenta o desvi des vioo médio ed io.. Tab. 3.2: RMSD

percentual para os dois casos considerados

Incompr Inco mpress´ ess´ıvel ıvel

Incompr Inco mpress´ ess´ıvel ıvel

Compres Com presss´ıvel

Dens. Igu Iguais

Dens. Dif.

Sem Retreinamento

3, 06

8, 83

11, 17

Retreinamento na 1 a sim.

3, 18

3, 59

4, 90


2, 33

8, 93

9, 49


3, 53

3, 38

4, 97


2, 62

6, 99

7, 72


1, 92

4, 60

4, 65


1, 98

5, 11

5, 72

Outro importante resultado diz respeito aos tempos de simula¸c˜ cao. aõ. Vale ale a pena pena utilizar um modelo de ordem reduzida que reproduz os resultados do modelo de alta fidelidade consumindo menos tempo de simula¸c˜ cao. a˜ o. D´ a-se a-se o nome de speedup à razão ao entre os tempos de simula¸c˜ cão ao de alta e baixa fidelidade. O modelo TPWL/POD TPWL/POD foi desenvol desenvolvido vido em Matlab, Matlab, equanto equanto que o SIMPAR SIMPAR é escrito em FORTRAN. Isto pode limitar um pouco o speedup, speedup, pois o Matlab, apesar de possuir muito boas bibliotecas, com fun¸c˜ coes o˜es e métodos etodos desenvolvidos desenvolvidos para serem rápidos, apidos,


82

acaba se mostrando mostrando lento por ser interpret interpretado. ado. O computador em que estes estudos foram foram realizados consiste em um Intel Xeon X5680 3, 3 ,33GHz, com 24GB de RAM, 12 núcleos ucleos e sistema sistema operacional operacional Windows Windows 7. Uma simula¸c˜ cao aõ do modelo de reservatórios orios aqui estudado estudad o para o caso compress´ıvel ıvel com densidades diferentes diferentes utilizando o SIMPAR levou em média edia 1873 s, considerando-se as 27 realizadas no estudo. Por outro lado, as simula¸c˜ coes o˜es correspondentes do TPWL/POD custaram em média edi a 4,8 s, o que resulta resulta em um speedup um speedup de de 390. 3 90. Além em disso, d isso, após os o retreinamento, o custo do TPWL/POD TPWL/POD passou a ser em média edia 6, 6 ,4 s, resu result ltan ando do em um speedup de 292. Este aumento no custo computacional do TPWL/POD após os o retreinamento era esperado e é devido ao maior número umero de estados gravados, o que complica o processo de busca do mais próximo. oximo. Entretanto, note que as avalia¸c˜ cões oes do SIMPAR podem ser realizadas em paralelo (a rotina já contempla isso), mas as do TPWL/POD devem ser feitas sequencialmente devido a falta de licen¸cas cas do módulo odulo de processamento paralelo do Matlab. Outro tempo importante a ser avaliado é o tempo de treinamento. Considerando o mesmo caso compress´ compress´ıvel com densidades diferentes, diferentes, a simula¸c˜ cão ao do SIMPAR com controles aleatórios orios para para treina treinamen mento to custou custou 1830 s. O process processoo de c´ alculo alculo da base base do POD e dos estado estadoss reduzi reduzidos dos levou levou 44 s consid considera erando ndo a leitur leituraa dos mapas, o process processoo de c´ alculo alculo das derivadas derivadas reduzidas levou 580 s, tempo consumido em grande parte pela leitura das derivadas e a simula¸c˜ cão ao do TPWL/POD da mesma sequência encia de controles de treinamento utilizou 5 s. Sendo assim, o processo de treinamento custou mais ou menos 1, 1,4 simul si mula¸ a¸c˜ cões. oes. O processo pro cesso de retreinamento retreinamento é um pouco mais caro. Não ao há a necessidade de uma simula¸c˜ cão ao do SIMPAR exclusiva para o retreino, já que a ideia é utilizar uma das simula¸c˜ coes o˜es já executadas durante o processo de otimiza¸c˜ cão. ao. En Entre tretan tanto, to, há a necessidade de realizar o POD considerando os mapas utilizados no primeiro treinamento e os da nova simula¸c˜ cao. aõ. Desta forma consumiu-se nessa etapa um tempo de em torno de 80s 80 s. O cálculo alculo das derivadas reduzidas durante o retreinamento retreinamento é bem b em mais caro, consumindo em torno de 1200s 1200 s devido a uma maior quantidade de derivadas derivadas a serem lidas e reduzidas. Além em disto, é necessário ario recalcular o VPL utilizando o TPWL/POD retreinado em todos os pontos nos quais já havia sido calculado anteriormente, o que pode custar bastante, dependendo da quantidade de avalia¸c˜ coes o˜es a serem realizadas, já que as avalia¸c˜ coes o˜es do modelo reduzido devem ser feitas sequencialmente. A ultima u ´ ltima análise alise realizada neste estudo diz respeito ao efeito de estabiliza¸cão ao causado pelo uso da proje¸c˜ cao aõ de Petrov-Galerkin no POD. Inicialmente, foi utilizada a proje¸c˜ cao aõ


83

de Bubnov-Galerkin, como na maioria dos trabalhos sobre o assunto (Cardoso, 2009; He, 2010). Ent˜ ao, observou-se uma tendˆ encia muito forte à instabilidade. Entretanto, alguns trabalhos demonstraram que o uso da proje¸caõ de Petrov-Galerkin (equivalente à solu¸caõ por m´ınimos quadrados) leva à estabiliza¸caõ do processo (Carlberg et al., 2009; He and Durlofsky, 2013). Ao utilizar-se tal recurso não se experimentou em nenhuma simula¸caõ qualquer problema de instabilidade. A Figura 3.8 mostra sobrepostos trˆ es gr´ aficos. Todos representam a vaz˜ a o de óleo dos quatro po¸cos no modelo SPE10 Modificado com compressibilidades e diferen¸ca de densidades e com os controles da simula¸caõ 1 utilizados no estudo de acurácia. As linhas cheias representam as vazões de óleo dos po¸cos produtores calculada pelo SIMPAR, as cruzes (x) são as vazões calculadas pelo TPWL/POD estabilizado pelo uso da pro je¸caõ de PetrovGalerkin (pré-multiplica¸caõ das derivadas por ΦT (J i+1 )T ) e os c´ırculos (o) representam as vazões calculadas pelo TPWL/POD com proje¸cão de Bubnov-Galerkin (pré-multiplica¸caõ por ΦT ).

6000 W1(SIMPAR) W2(SIMPAR) W3(SIMPAR) W4(SIMPAR) W1(TPWL−Petrov) W2(TPWL−Petrov) W3(TPWL−Petrov) W4(TPWL−Petrov) W1(TPWL−Bubnov) W2(TPWL−Bubnov) W3(TPWL−Bubnov) W4(TPWL−Bubnov)

5000

4000

3000

2000

1000

0 0

500

1000

1500

2000

2500

3000

3500

Fig. 3.8: Vazã o de óleo dos po¸cos - SPE10 - Modificado compress´ıvel e com densidades de fluidos diferentes


84

Note que a vazão calculada pelo TPWL/POD com a proje¸caõ de Petrov-Galerkin mostra alguns desvios em rela¸caõ a` obtida pelo simulador. De fato, o erro no VPL para este caso foi o maior entre todas as 27 simula¸cões sem retreinamento, atingindo 21 %. Entretanto, percebe-se que o resultado é estável e os erros podem ser atribu´ıdos às varia¸cões frequentes dos controle na simula¸cão de treinamento. Por outro lado, a vazão calculada utilizando a outra proje¸caõ se mostra verdadeiramente instá vel. O VPL calculado n˜ ao faz nenhum sentido e quando se observam os mapas de pressão e satura¸cões, muitos valores negativos ou acima dos limites podem ser encontrados. Este problema de instabilidade foi reportado na literatura (Gildin et al., 2013) e depende do n´ umero de autovalores escolhidos (ver (He, 2010), onde um estudo mais detalhado foi realizado). Por outro lado, a solu¸caõ de utilizar-se a proje¸cão de Petrov-Galerkin também foi apontada anteriormente (Carlberg et al., 2009; He and Durlofsky, 2013) e se mostrou bastante eficaz. Não foram observados problemas de instabilidade em nenhuma execu¸caõ do TPWL/POD com qualquer modelo quando utilizou-se a proje¸caõ de Petrov-Galerkin. ´ Isto não pode ser conclusivo, mas é um bom indicativo do papel estabilizador da técnica. E necess´ ario fazer um estudo mais cuidadoso sobre a estabilidade do método, aplicando técnicas apropriadas de análise numérica.

˜ EM OTIMIZAC ˜ 4. APLICAC ¸ AO ¸ AO Em diversos ramos da engenharia os estudos de otimiza¸caõ têm se tornado indispensá veis. Na a´rea de gerenciamento de reservatórios os estudos de otimiza¸cão vêm se tornando cada vez mais comuns. Os objetivos podem ser variados, como a determina¸caõ das posi¸cões para perfura¸cão de novos po¸cos, ou a defini¸cão da sequência o´tima de controles a serem aplicados aos po¸cos. A fun¸cão-objetivo utilizada tamb´ em pode variar, como, por exemplo, minimizar a produ¸cã o de água instantânea ou maximizar o valor presente l´ıquido de todo o tempo de vida do empreendimento de produ¸cão do o´leo do reservatório. Entretanto, os modelos de simula¸caõ utilizados nesta área possuem elevado custo computacional. Ultimamente, cada vez mais detalhes são inclu´ıdos nos modelos e, apesar dos avan¸cos dos computadores, o tempo de avalia¸caõ de um modelo destes é razoavelmente alto. Por outro lado, a quantidade de avalia¸cões que devem ser realizadas durante o processo de otimiza¸caõ, por mais eficiente que seja o algoritmo, pode ser bastante grande. Sendo assim, fazem-se necessários esfor¸cos para reduzir o custo computacional total do processo de otimiza¸caõ, reduzindo-se assim o tempo para a obten¸cão dos resultados. Uma alternativa que vem sendo empregada é a utiliza¸caõ de modelos de ordem reduzida substitutos da simula¸caõ que aproximam o seu resultado com boa acurácia, mas são executados em um tempo bem menor. O algoritmo de otimiza¸caõ sequencial aproximada (Sequential Approximate Optimization - SAO) é um daqueles baseados em modelos substitutos (Surrogate Based Optimization – SBO) e realiza a otimiza¸cão de uma fun¸caõ cuja avalia¸cão é custosa através de uma sequência de otimiza¸co˜es de modelos substitutos que têm validade em uma região limitada e que vão sendo melhorados à medida que o processo avan¸ca (Giunta and Eldred, 2000). Uma aplica¸cão deste algoritmo ao problema de otimiza¸caõ do gerenciamento de reservatórios de petróleo pode ser encontrada em (Horowitz et al., 2013). Sua versão tradicional utiliza como modelo de ordem reduzida uma abordagem

85

86 baseada em processos estocásticos (Jones et al., 1998). Tal abordagem é conhecida por kriging (Forrester et al., 2008) ou krigagem e utiliza pontos que tenham sido previamente amostrados e avaliados no modelo de alta fidelidade. Este modelo é exclusivamente baseado em dados e, portanto, é dito data driven . Entretanto, o algoritmo aqui desenvolvido permite utilizar um segundo modelo substituto com base na f´ısica do problema. Sendo assim, a krigagem é utilizada somente para corrigir poss´ıveis erros que este modelo f´ısico possa vir a cometer na estima¸cão do modelo de alta fidelidade. Este é chamado de modelo de baixa fidelidade. Portanto, o metamodelo ou modelo substituto que é de fato utilizado pelo algoritmo SAO passa a ser a soma do modelo baseado no kriging com o modelo de baixa fidelidade. A rotina desenvolvida permite a utiliza¸cão dos dois modelos em conjunto — o que é conhecido como multifidelidade — mas também permite que o usuário opte por utilizar SAO tradicional, caso um modelo f´ısico de baixa fidelidade não esteja dispon´ıvel. Todo o algoritmo foi desenvolvido em Matlab (Mathworks, 2012), considerando-se como fun¸caõ de alta fidelidade o simulador de reservatórios SIMPAR, descrito no Cap´ıtulo 2 e como fun¸cão de baixa fidelidade o modelo de ordem reduzida TPWL/POD desenvolvido no Cap´ıtulo 3. Em ambos os modelos a sa´ıda utilizada como fun¸cão-objetivo no processo de otimiza¸caõ é o Valor Presente L´ıquido definido na se¸cão 2.4 e os controles são as pressões de fundo dos po¸cos produtores. O algoritmo aqui desenvolvido recebe o nome de SAO Multifidelidade (MSAO é a sigla em inglês). Foi concebido de forma a permitir o acoplamento com qualquer tipo de modelo e, portanto, para a integra¸caõ destes modelos, foram constru´ıdas rotinas que recebem como parâmetros as variáveis de projeto da otimiza¸caõ normalizadas, aplicando-as aos simuladores e recebem os resultados de VPL das simula¸co˜es, normalizando-os para que sejam avaliados pelo otimizador. O presente cap´ıtulo come¸ca fazendo primeiramente algumas considera¸co˜es teóricas sobre as técnicas que são utilizadas pelo algoritmo desenvolvido e sobre as diretrizes nele aplicadas. Logo em seguida descreve-se o algoritmo MSAO e o seu desenvolvimento no Matlab. Por fim são apresentados alguns resultados obtidos com os dois diferentes modelos de reservatórios descritos no Cap´ıtulo 2.

4.1. Algoritmos de Otimiza¸caõ baseados em Modelos Substitutos

4.1

87

Algoritmos de Otimiza¸c˜ ao baseados em Modelos Substitutos O algoritmo de otimiza¸caõ aqui desenvolvido é uma adapta¸caõ do algoritmo de

otimiza¸caõ sequencial aproximada tradicional. Como em qualquer algoritmo baseado em modelos substitutos ou proxy models, ocorrer´ a a aproxima¸caõ através de modelos simplificados de um problema de elevado custo computacional, representado em um primeiro momento por um modelo de alta fidelidade. Este modelo de alta fidelidade pode ser, por exemplo, um modelo de simula¸caõ de reservatórios e pode representar a fun¸caõ-ob jetivo ou as restri¸cões do problema de otimiza¸cã o. Já os modelos substitutos podem ser dos seguintes tipos, assim classificados em (Afonso et al., 2011): 1. Aproxima¸caõ funcional – baseada na interpola¸caõ de respostas obtidas a partir do modelo de alta fidelidade através de uma fun¸caõ especificada; esta aproxima¸caõ também é dita baseada em dados ou data driven . Exemplos deste tipo de aproxima¸co˜es são as redes neurais, fun¸cões de base radial (radial basis fuctions ) a krigagem ou kriging . 2. Aproxima¸caõ f´ısica – baseada na f´ısica do problema real, ou seja, o modelo de ordem reduzida é obtido a partir de alguma simplifica¸caõ no equacionamento do problema de alta ordem. Um exemplo para este tipo é a combina¸caõ da técnica TPWL (Trajectory Piecewise Linearization ) com POD (Proper Orthogonal Decomposition ) . A primeira lineariza o problema real em torno de estados gravados a partir da simula¸caõ de alta fidelidade enquanto a segunda reduz a dimensão do problema através da decomposi¸caõ em dire¸cões principais e elimina¸cão das dire¸cões menos importantes. 3. Aproxima¸caõ mista – é um misto dos dois tipos anteriores, onde se combinam técnicas ´ este tipo de aproxima¸caõ que foi considerada neste f´ısicas com técnicas funcionais. E trabalho. Existe uma classe de algoritmos utilizados para otimiza¸caõ que são ditos baseados em modelos substitutos (Surrogate Based Optimization - SBO). Estes algoritmos se baseiam na aproxima¸cão do modelo de alta fidelidade, complexo e de elevado custo computacional, de alguma das três formas descritas acima. Um destes algoritmos é o SAO (Sequential Approximate Optimization ) (Afonso et al., 2011). A ideia básica é que essa aproxima¸caõ é confiável dentro de uma regiã o de confian¸ca cujo raio será avaliado a cada itera¸caõ. Sendo assim, a otimiza¸caõ é realizada sobre

4.2. Modelos de Kriging

88

a fun¸cão aproximada, simples e de baixo custo computacional. O ponto ótimo obtido com o modelo aproximado é aplicado ao modelo real e os resultados são comparados segundo algum crit´ erio. Dependendo da proximidade dos dois resultados, o ponto ótimo é aceito como centro da região de confian¸ca da próxima itera¸cão do algoritmo ou não e o raio da nova região é atualizado, podendo ser reduzido, aumentado ou mantido. Um novo modelo substituto é constru´ıdo para a nova regi˜ ao de confian¸ca e a otimiza¸cão é novamente realizada. Desta sequência de otimiza¸co˜es advém o nome do algoritmo. Note que no caso deste trabalho, não haverá outras restri¸co˜es além dos limites das vari´ aveis de projeto, e portanto só se faz necess´ ario a aproxima¸caõ do modelo de simula¸cão utilizado como fun¸caõ-objetivo. Mas num caso geral, onde existam restri¸co˜es não-lineares baseadas em simula¸cão, modelos substitutos podem ser utilizados tamb´ em para as restri¸cões. No caso do SAO comumente encontrado na literatura, um modelo reduzido comumente utilizado é o kriging (ou krigagem em português), um tipo de aproxima¸cão funcional que interpola resultados obtidos a partir do modelo de alta fidelidade. Portanto, uma desvantagem direta é que este modelo real deve ser simulado algumas vezes para obten¸caõ de pontos que serão utilizados na constru¸cão do modelo substituto, fazendo-o sofrer da “maldi¸cã o de dimensionalidade”, pois quanto mais mais variáveis o problema considerar, mais simula¸co˜es de alta fidelidade devem ser realizadas para constru¸cão da aproxima¸caõ. O algoritmo MSAO aqui desenvolvido utilizará uma forma de aproxima¸cão mista como modelo substituto, a qual é resultado da composi¸caõ de um modelo f´ısico, o TPWL, desenvolvido no cap´ıtulo 3 e uma corre¸cão baseada em krigagem. O MSAO, assim como o SAO tradicional, é baseado em regiões de confian¸ca e possui as mesmas regras para atualiaza¸cão do centro e do tamanho da região de confian¸ca, necessitando fazer algumas amostras no interior desta região. Foi nomeado Multifidelity Sequential Aproximate Optmization por se basear em dois n´ıveis de fidelidade.

4.2

Modelos de Kriging A krigagem, do nome do criador Krige (Krige, 1951), faz a interpola¸cão de dados

através da associa¸cão de um polinômio (regressão polinomial) a uma fun¸cão erro cuja matriz de covariância é conhecida. Nesta aproxima¸caõ, assume-se que os dados possuem correla¸caõ entre si e que esta correla¸caõ depende da distância (Forrester et al., 2008). Este tipo de modelo é bastante utilizado nas áreas da estat´ıstica e geo-estat´ıstica


89

e foi inicialmente introduzido por (Matheron, 1963) baseado na disserta¸cão de mestrado (Krige, 1951). O modelo considera a covariância entre os pontos da amostra usando uma correla¸caõ espacial através de uma distribui¸cão gaussiana em cada ponto. Devido a isto, métodos que utilizam esta técnica são facilmente influenciados pelo arranjo dos pontos das amostras (da Silva, 2009). A krigagem, diferentemente dos modelos de regressão tradicionais (que utilizam polinômios), em geral são capazes de prover aproxima¸cões globais a` resposta e também capturar tendências oscilatórias (Afonso et al., 2011). O modelo de krigagem é um interpolador que, de forma geral, pode ser expresso segundo a Equa¸caõ 4.1. k

f (x) =

�

β j N j (x) + Z (x)

(4.1)

j=1

Na defini¸cão do modelo de krigagem, N j (x) representam fun¸co˜es polinomiais ponderadas por coeficientes β j ajustados por regressão e Z (x) é um desvio localizado considerando uma fun¸caõ aleatória que normalmente segue uma distribui¸caõ gaussiana. Suponha que exitem m observa¸co˜es da fun¸caõ a ser a justada. Além disso o vetor de entradas do modelo tem dimensão n. O termo aleat´ orio Z (x) é constru´ıdo assumindo média igual a zero e covariância dada pelo variograma da Equa¸caõ 4.2, onde os ´ındices i e j representam as diferentes observa¸cões e portanto 1 < i, j < m..





( 

cov Z (xi ), Z (x j ) = σ 2 R θ, xi , x j A matriz de correla¸caõ R tem dimensão m

(4.2)

× m e é definida na Equa¸cão 4.3 e

correlaciona xi e x j , i, j = 1 . . . m. Esta utiliza um tipo de fun¸cão de correla¸caõ dependente dos parâmetros θk , k = 1 . . . n inicialmente desconhecidos. n

ij

R

( ∏ ( i

Rl θl , xil , x jl

j

θ, x , x =

l=1



(4.3)

Existem várias formas para as fun¸co˜es de correla¸caõ R l que podem ser utilizadas na montagem de R . A fun¸cão de correla¸caõ aqui considerada é conhecida como exponencial gaussiana e tem a forma apresentada na Equa¸caõ 4.4.


90

Rk = exp

−

j 2 l

θk (xil

−x )



(4.4)

ˆ x), Uma estimativa do modelo a ser aproximado em qualquer ponto do espa¸co, f ( considerando as respostas já obtidas da fun¸cão real f (x) pode ser calculada como na Equa¸caõ 4.5. ˆ x) = E [[f (x) f (x1 ), . . . , f ( xm )]] f (

|

(4.5)

´ Para avaliar a qualidade desta aproxima¸caõ faz necessário uma medida de erro. E escolhido então o Erro Médio Quadrático (Mean Square Error (MSE)), que mede a diferen¸ca entre a fun¸cão aproximada pelo modelo de krigagem e a fun¸caõ real em todo o espa¸c o de projeto conforme exposto na Equa¸caõ 4.6 (Horowitz et al., 2012).

MSE = E

��

f ˆk (x)

−

f (x)

�� 2

(4.6)

O estimador f ˆk (x) que minimiza o MSE pode ser calculado conforme a Equa¸caõ 4.7. ˆT + r(x)R f ˆk (x) = N (x)β

1

−

�

f S

� − N ˆT β

(4.7)

Para compreensão deste resultado fazem-se necessárias as defini¸co˜es listadas abaixo.

• N(x) = [N (x) . . . N (x)] é um vetor contendo os polinômios utilizados na regressão; • β ˆ é um vetor de coeficientes do modelo de regressão; • f = [f . . . f ] é um vetor com os valores do modelo real avaliado nos pontos amos1

S

1

k

m

T

trados;

• N é uma matriz m × k que contém os polinômios escolhidos para a cria¸caõ do modelo

de regressão num´ erica aplicados nas amostras com o seguinte arranjo apresentado na Equa¸caõ 4.8.

N =

 

N 1 (x1 ) .. .

···

N k (x1 ) .. .

N 1 (xm )

···

N k (xm )

 

(4.8)


91

• r(x) = [R(θ, x, x ) . . . R(θ, x, x 1

m

)]T é o vetor de correla¸cão entre os pontos da amostra

à qualquer ponto dentro de espa¸co de projeto. ˆ e a covariância σ 2 são obtidos usando o método dos Os coeficientes da regressão β

m´ınimos quadrados segundo a Equa¸caõ 4.9.

( N N  N � − N � �

ˆ = β σ2 =

1 m

T

R

f S

1

−

1

−

ˆT β

T

T

R 1 f S

R

1

−

f S

−

� − N

(4.9)

ˆT β

Note que r e R dependem do vetor de parâmetros θ. Para determiná-los utilizase a estimativa de máxima verossimilhan¸ca, obtida através da solu¸caõ de um problema de minimiza¸caõ definido na Equa¸cão 4.10. θ = arg min 12 [m ln σ2 + ln R ] s. a. θ

≥0

| |

(4.10)

Ao resolver este problema de otimiza¸caõ que está acoplado a Equa¸cã o 4.9, o modelo kriging está constru´ıdo. Resumindo, a krigagem funciona como mostrado na figura 4.1. Na figura 4.1(a) é poss´ıvel observar uma aproxima¸cão de um conjunto de pontos por regressão de um polinômio de segundo grau, uma parábola. Note que a curva não honra nenhum deles, possuindo erro maior que zero para todos, já que estes pontos não recaem sobre uma curva polinomial. Nas figuras 4.1(b) e 4.1(c) observa-se a aproxima¸caõ por krigagem. Note que a par´ abola foi combinada a um desvio local associado a` parte estocástica do kriging que deforma a curva polinomial (pontilhada) de tal forma a fazê-la honrar todos os pontos. O efeito de tal deforma¸caõ depende do raio de correla¸caõ utilizado no kriging . Na figura 4.1(b) utilizou-se um raio de correla¸caõ menor enquanto na figura 4.1(c) utilizou-se um raio maior. Observa-se uma diferen¸ca na distância em que a krigagem atua. Observa-se também que, com um raio mais longo, o método altera a extrapola¸cão polinomial, n˜ ao honrando mais a parábola além dos limites dos pontos conhecidos.

4.3. Amostragem por Hipercubo Latino

Y

92

10

10

5

5

0

0

Y

−5

−5

−10 −20

0

20

40

60

80

100

−10 −20

120

0

20

40

X

60

80

100

120

X

(a) Regressão polinomial

(b) Krigagem - correla¸cão de curto alcance

10

5

Y

0

−5

−10 −20

0

20

40

60

80

100

120

X

(c) Krigagem - correla¸caõ de longo alcance

Fig. 4.1: Interpola¸caõ por krigagem - regressão polinomial adicionada a deforma¸co˜es locais causadas por uma fun¸caõ aleatória com correla¸caõ espacial definida

4.3

Amostragem por Hipercubo Latino Para realizar a krigagem e, portanto, construir o interpolador, como mencionado

anteriormente, faz-se necessário avaliar-se o modelo de alta fidelidade em diferentes pontos do espa¸co de projeto dentro da região de confian¸ca, de tal forma que precisa-se amostrar tal espa¸co. Para tanto existem alguns métodos de amostragem (Forrester et al., 2008; Afonso et al., 2011), podendo-se entre eles destacar: Quasi–Monte Carlo (QMC), Centroidal Voronoi Tesselation (CVT) e Latin Hypercube Sampling (LHS). A técnica aqui escolhida foi a LHS (Giunta, 2002). Para obter o conjunto de amostras pelo método LHS, o intervalo de cada variável deve ser dividido em m subintervalos equiprováveis, onde m é o n´ umero de amostras. Se

4.4. Algoritmo SAO Multifidelidade

93

considerarmos n variáveis, esta divisão resulta em mn subdivisões no espa¸co de projeto. Sendo assim, m amostras são selecionadas aleatoriamente dentro do espa¸co de projeto satisfazendo as seguintes restri¸co˜es: cada amostra é aleatoriamente posicionada dentro de uma parti¸caõ do dom´ınio e para toda proje¸caõ unidimensional das amostras e parti¸cões, haver´ a uma e somente uma amostra em cada parti¸caõ (Horowitz et al., 2012).

Fig. 4.2: Hipercubo latino 2D com 4 amostras A figura 4.3 mostra como exemplo o que seria um hipercubo latino em 2D, isto é, para duas variáveis (x1 e x 2 ). Neste caso o número de amostras m = 4 e portanto 4 parti¸cões são constru´ıdas nos dois eixos. Isso resulta em 16 “subespa¸cos” dos quais 4 são escolhidos de forma a satisfazer às duas restri¸cões. As amostras são representadas pelos pontos pretos posicionados no interior de cada quadrado e foram sorteados aleatoriamente.

4.4

Algoritmo SAO Multifidelidade O algoritmo de otimiza¸caõ sequencial aproximada - SAO aqui desenvolvido possui

uma diferen¸ca básica com rela¸cão ao tradicionalmente explorado na literatura. Este não utiliza a krigagem dos dados diretamente como modelo substituto. Um modelo de ordem reduzida, que no caso desta disserta¸caõ é representado pelo TPWL apresentado no cap´ıtulo 3, é utilizado e, portanto, a krigagem é aplicada à diferen¸ca entre as respostas do modelo real e deste modelo substituto f´ısico. Sendo assim, existem duas fun¸cões que devem ser avaliadas: o modelo real, chamado de modelo de alta fidelidade e o modelo reduzido, chamado de modelo de baixa fidelidade. O proxy ou metamodelo utilizado pelo algoritmo de otimiza¸caõ é, portanto, a soma do modelo de baixa fidelidade com a corre¸caõ por kriging . A rotina


94

desenvolvida tamb´ em permite o uso do algoritmo SAO tradicional caso n˜ ao haja modelo de baixa fidelidade dispon´ıvel. O objetivo é resolver o problema de otimiza¸caõ descrito na Equa¸caõ 4.11, onde f (x) representa o modelo alta fidelidade que no caso em questão é o Valor Presente L´ıquido calculado por uma simula¸caõ de reservatórios utilizando o SIMPAR. x é o vetor de variáveis de projeto, isto é, a pressão de fundo de po¸co a ser aplicada como controle dos po¸cos em cada ciclo de controle. xl , xu são os limites superiores e inferiores do espa¸co de projeto, isto é, os valores máximos e m´ınimos para as pressões de fundo em cada po¸co e cada ciclo. f (x)

max s.a. xl

(4.11)

≤ x ≤ x

u

Expressando matematicamente, a ideia por trás do SAO multifidelidade é converter o problema de otimiza¸caõ em uma sequência de problemas aproximados da forma expressada na Equa¸cão 4.12. Note que a fun¸cão de alta fidelidade f (x) é substitu´ıda pelo ˆ x) resultante da composi¸caõ de um modelo de ordem metamodelo ou modelo aproximado f ( reduzida e um modelo de corre¸caõ e é válido no interior de uma região de confian¸ca. F k (x) é um modelo de baixa fidelidade que aproxima a f´ısica de f (x) com ordem reduzida. No caso deste trabalho F k (x) representa o VPL calculado pela simula¸cão do TPWL e pode ser atualizado ou retreinado caso seja necessário. δ k (x) é um modelo de krigagem constru´ıdo a partir da diferen¸ca entre o de alta fidelidade f (x) e o de baixa F k (x) servindo como elemento de corre¸caõ ao de baixa fidelidade. xkl , xku são respectivamente os limites superiores e inferiores da região de confian¸ca da itera¸caõ e devem ser atualizados a cada itera¸caõ do algoritmo. k e

xkc , ∆k são respectivamente o centro e o raio da região de confian¸ca da itera¸caõ k e também são recalculados a cada passo. max f ˆk (x) = F k (x) + δ k (x) s.a. xl

k l k k c l xku = x kc

k u

≤ x ≤ x ≤ x ≤ x x = x − ∆ k

u

(4.12)

+ ∆k

k = 0, 1, 2, . . . O algoritmo utilizado para a otimiza¸caõ em cada um dos problemas da sequência é o de programa¸cão quadrática sequencial (em inglês, Sequential Quadratic Programing -


95

SQP) (Biggs, 1979), já implementado na fun¸caõ fmincon do Matlab. Entretanto, qualquer algoritmo de otimiza¸caõ pode ser empregado. A região de confian¸ca consiste em um paralelep´ıpedo multidimensional no espa¸co de projeto cujas dimensões podem ser iguais ou diferentes dentro do qual se supõe boa a aproxima¸caõ feita pelo metamodelo. Genericamente seu tamanho (dado por uma ou mais medidas se as dimensões forem diferentes) é conhecido como raio. Além do raio, a região de confian¸ca possui um centro. O primeiro raio e centro são pré-definidos e podem ter influência no desempenho do algoritmo. ` medida que o algoritmo avan¸ca, a cada itera¸caõ, é necessário atualizar o centro e A do raio região de confian¸ca. Para tanto, é necessário medir a qualidade da aproxima¸caõ que o metamodelo está realizando. Sendo assim, o ponto ótimo obtido com o metamodelo é aplicado ao modelo de alta fidelidade. De posse das respostas dos dois modelos, calcula-se o parˆ ametro de qualidade da aproxima¸caõ ρk através da Equa¸caõ 4.13 onde f (xk ), f (xk ), f ˆk (xk ), f ˆk (xk ) f

c

∗

c

∗

são respectivamente o modelo de alta fidelidade avaliado no centro e no ponto ótimo calculado na itera¸cão k e o modelo aproximado (composto) avaliado no centro e no ponto ótimo calculado na itera¸caõ k.

ρkf

f (xkc ) = f ˆk (xkc )

k

− f (x ) − f ˆ (x )

(4.13)

∗

k

k ∗

A regra aplicada pelo algoritmo para atualiza¸caõ do raio da região de confian¸ca é de fato a “arte” do método. Neste trabalho utilizou-se a mesma regra prosposta em (Giunta and Eldred, 2000) exposta na equa¸caõ 4.14.

k+1

∆

=

  

0, 25∆k , 0, 25∆k ,

ρk 0

∆k , 0, 25 2∆k ,

≤ ≤

ρk ρk ρk

≤ ≤ ≤ ≥

0 0, 25 0, 75

(4.14)

0, 75

A regra para atualiza¸caõ do centro da região é também proposta em (Giunta and Eldred, 2000) e está descrita na Equa¸cão 4.15 onde x k+1 é o centro da próxima itera¸cão, x kc c é o centro da itera¸caõ corrente e xk é o ponto ótimo calculado utilizando o metamodelo na ∗

itera¸cão corrente.


xk+1 c

96

=

{

xk , ρ > 0 ∗

xkc , ρ < 0

(4.15)

Fig. 4.3: Esquema de amostragem e atualiza¸caõ da região de confian¸ca no algoritmo MSAO A figura 4.4 esquematiza como o algoritmo MSAO funciona na busca pelo ótimo. Nela está representado um problema de duas variáveis de projeto (fun¸cão banana de Rosenbrok). A estrela representa o ponto ótimo. A cada itera¸cão o algoritmo considera uma região de confian¸ca representada pelos retângulos. Faz-se a amostragem de um conjunto de pontos no interior desta região, avaliando-se os modelos de alta e baixa fidelidade (SIMPAR e TPWL/POD) em todos os pontos amostrados. Constr´ oi-se o modelo de corre¸caõ por krigagem utilizando-se a diferen¸ca entre os resultados nos dois n´ıveis de fidelidade. Otimiza-se no interior da região de confian¸ca utilizando-se o metamodelo representado pelo modelo de baixa fidelidade corrigido. O resultado ao fim de cada itera¸caõ é um ponto candidato a o´timo que deve ser aplicado ao modelo de alta fidelidade e ao de baixa fidelidade corrigido. Com base na metodologia aqui exposta, representada pelas Equa¸cões 4.14 e 4.15, o algoritmo calcula as dimens˜ oes e o centro da próxima região de confian¸ca, repetindo todo o processo na próxima itera¸cão. Os critérios de parada que são testados pelo algoritmo e que levam à parada do processo estão listados a seguir e, para que qualquer um dos critérios acima leve à parada do processo, é necessário que a situa¸cão ocorra em um n´ umero de itera¸co˜es consecutivas definido pelo usuário. 1. Distância m´ınima entre os pontos de ótimo atingida;


97

2. Varia¸caõ m´ınima da fun¸caõ-objetivo atingida; 3. Raio m´ınimo da região de confian¸ca atingido; 4. N´ umero máximo de itera¸cões atingido. O algoritmo MSAO também foi concebido de forma a permitir que o modelo de ordem reduzida seja retreinado ao longo do processo de otimiza¸caõ. Na realidade, observouse uma tendência de redu¸cão da qualidade do metamodelo ao se aproximar do ponto ótimo. Isto pode ser atribu´ıdo ao pequeno raio da regi˜ ao de confian¸ca que pode levar a varia¸co˜es com sinais diferentes do VPL calculado pelo simulador e o calculado pelo metamodelo. Neste caso, não se deseja o retreinamento, mas sim parar o algoritmo. Entretanto, o parâmetro de qualidade ρkf pode ficar com sinal negativo, o que pode ser confundido com qualidade ruim do metamodelo e portanto gerar a necessidade de retreinar o modelo de ordem reduzida. Com o objetivo de avaliar se o defeito do metamodelo é o TPWL/POD, foi proposto o parâmetro de qualidade E T P W L da Equa¸caõ 3.34. Este representa uma soma acumulada das distâncias entre os estados calculados pelo modelo reduzido e os estados gravados utilizados na lineariza¸caõ, sendo considerado uma medida da distância da trajetória simulada pelo TPWL/POD e a trajetória em torno da qual a lineariza¸caõ foi realizada. Devido ao custo computacional elevado do retreinamento, espera-se que este só seja necessário quando esta medida for grande. Portanto, só se cogita retreinar quando ρ kf < ρ min e E T P W L > E max , onde ρmin e E max são definidos pelo usuário devido a experiência com o problema a ser otimizado. Desta forma, o algoritmo SAO multifidelidade aqui desenvolvido pode ser resumidamente descrito como a seguir. Considere x0c , a coordenada do centro e ∆0 o raio da primeira região de confian¸ca. Considere tamb´ em que já foi executada uma simula¸cã o de treinamento do SIMPAR com uma sequência aleatória de controles. Esta foi utilizada para a constru¸caõ do modelo TPWL/POD inicial, F 0 . 1. Inicializa o contador de itera¸cões k = 0. 2. Calcula o modelo de alta e baixa fidelidade (VPL atrav´ es da simula¸caõ SIMPAR e TPWL/POD respectivamente) no centro, f (x0c ) e F 0 (x0c ). 3. Calcula os limites da regi˜ ao de confian¸ca xkl = x kc

k

− ∆

e xku = x kc + ∆k . Observe que

os limites em cada itera¸caõ não devem ultrapassar os limites do espa¸co de projeto x l e

xu .


98

4. Amostra N pontos dentro desta região de confian¸ca xki , 1

≤ i ≤ N (quantidade definida

pelo usuário – recomendado N = 2N var + 1 para a primeira itera¸caõ, onde N var é o número de variáveis e N = N proc para as demais, onde N proc equivale ao número de processadores). Estas amostras são adicionadas ao conjunto de todas as amostras obtidas nas itera¸cões anteriores. 5. Calcula o modelo de alta e baixa fidelidade para todos os pontos amostrados na itera¸caõ k, f (xki ) e F k (xki ), 1

≤ i ≤ N . Acumulam-se estas avalia¸cões das fun¸cões ao conjunto

de todas as anteriores.

6. Calcula a diferen¸ca entre as respostas de alta e baixa fidelidade δ ik = f (xki )

k

k i

− F (x ).

7. Faz o kriging considerando todos os pontos do conjunto de amostras no interior da região de confian¸ca da itera¸caõ k (amostras x diferen¸ca das respostas). Caso não haja pontos suficientes no interior da região de confian¸ca, são utilizados os pontos mais pr´ oximos de todos os já calculados até que se complete a quantidade necessária. Constroi-se assim o modelo δ k (x) e o metamodelo f ˆk (x) = F k (x) + δ k (x). 8. Otimiza o metamodelo f ˆk (x) dentro da região de confian¸ca (problema mostrado na Equa¸caõ 4.16) utilizando o algoritmo de programa¸caõ quadrática sequencial (SQP) determinando o ponto candidato a ótimo x k . ∗

k

x = ∗

{

max

f ˆk (x)

s.a. xkl

≤x≤x

(4.16)

k u

9. Verifica convergência ou algum dos critérios de parada definidos na lista anterior. Caso algum dos crit´ erios seja atingido, incrementa o contador correspondente. 10. Calcula modelo de alta fidelidade no ponto de ótimo obtido com o metamodelo f (xk ). ∗

Al´ em do VPL, anota o erro do TPWL/POD calculado conforme Equa¸cã o 3.34. Se algum dos contadores de convergência atinge o valor definido pelo usuário, o algoritmo pára. 11. Caso o algoritmo não tenha atingido a convergˆ encia, simultaneamente com o passo 10, amostram-se N proc pontos x jk , 1

≤ j ≤ N

proc dentro

de uma região de confian¸ca

centrada no u ´ ltimo candidato a ótimo xk com o mesmo raio anterior ∆k . ∗

12. Avaliam-se a fun¸caõ de alta e baixa fidelidade f (x jk ) e F k (x jk ), 1

≤ i ≤ N

proc nos

pontos

amostrados simultaneamente com a avalia¸cão do candidato a o´timo f (xk ) do passo 10. ∗


99

Aplica-se a corre¸caõ δ k (x) ao modelo de baixa fidelidade. Acumulam-se estas avalia¸co˜es ao conjunto de todas as avalia¸co˜es anteriores. 13. Calcula parˆ ametro de qualidade do metamodelo ρ kf definido na Equa¸caõ 4.13. 14. Aplica a regra do algoritmo definida nas Equa¸co˜es 4.14 e 4.15 a partir do valor de ρkf , calculando novo centro x k+1 e raio ∆ k+1 da nova região de confian¸ca. c 15. Se ρkf < ρmin , onde ρmin é configurado pelo usuário, ativa-se o alerta para retreinamento que obrigar´ a a exporta¸cão dos mapas e derivadas na simula¸caõ do próximo candidato a o´timo xk+1 . ∗

16. Se ρkf < ρmin e o alerta de retreinamento já se encontrava ligado, realiza-se o processo de retreinamento com a simula¸caõ do candidato a o´timo atual xk . Desliga-se o alerta ∗

de retreinamento.

← k + 1. Retorna a 3, iniciando a próxima itera¸caõ.

17. k

4.5

Detalhes de Implementa¸ c˜ ao Todo o desenvolvimento foi realizado utilizando-se o Matlab. Para tanto dois

pacotes se fizeram necessários. O primeiro, utilizado para a otimiza¸caõ propriamente dita, é o Optimization Toolbox (Mathworks, 2012) e o segundo, responsável por realizar o processo de krigagem, é o DACE (Design and Analysis of Computer Experiments ) (Lophaven et al., 2002). Para otimiza¸caõ foi utilizada a fun¸caõ fmincon , aplicada ao metamodelo. Esta fun¸caõ permite o uso de restri¸co˜es lineares e não lineares e pode resolver o problema utilizando os algoritmos mais avan¸cados da atualidade. Entretanto, nenhuma restri¸cão, além dos limites da região de confian¸ca de cada itera¸caõ do SAO, foi aplicada e o algoritmo para otimiza¸caõ foi o SQP (Sequential Quadratic Programming ). Outro aspecto importante é que esta fun¸caõ tem por objetivo resolver um problema de minimiza¸cão. Portanto, como o problema em quest˜ ao é de maximiza¸caõ do Valor Presente L´ıquido, fez-se necessário considerar a fun¸caõobjetivo com sinal trocado. Esta é uma estrat´ egia muito comum para o uso de pacotes de otimiza¸caõ. Na constru¸caõ do modelo kriging utilizou-se a fun¸caõ dacefit do pacote DACE que recebe como parâmetros de entrada os pontos amostrados e suas respectivas avalia¸co˜es no

4.6. Desempenho do algoritmo MSAO

100

modelo real (neste caso diferen¸ca entre a avalia¸cão no modelo de alta e baixa fidelidade). Esta fun¸caõ é capaz de utilizar diferentes tipos de correla¸cão, mas neste trabalho só foi utilizada a de Gauss. Tamb´ em pode utilizar polinômio regressor constante (kriging ordinário) ou de primeiro ou segundo graus. Tanto a fun¸cão de correla¸cã o como o grau do polinˆ omio de regressão são escolhidos entre os parâmetros de entrada. A principal sa´ıda da fun¸caõ dacefit é uma estrutura com todos os parˆ ametros necess´ arios para reconstruir o modelo kriging que se ajusta aos pontos. Basicamente são os coeficientes do polinômio e os parâmetros da fun¸caõ de correla¸cão. A outra fun¸caõ do pacote DACE utilizada no trabalho é a fun¸caõ predictor . Esta fun¸caõ avalia o modelo de krigagem em um ponto dado como entrada. O modelo utilizado para avalia¸caõ tem que ter sido previamente ajustado e é representado por uma estrutura sa´ıda da fun¸caõ dacefit que deve ser passada como entrada para predictor . Para realizar a amostragem por hipercubo latino também foi utilizada uma fun¸caõ do pacote DACE chamada lhsamp. Entretanto, esta fun¸caõ somente amostra pontos em um hipercubo de aresta unitária. Portanto, fez-se necessário escalarem-se estas amostras para as dimens˜ oes da região de confian¸ca. Neste estudo a fun¸caõ de alta fidelidade é representada por uma simula¸caõ de reservatórios utilizando o simulador SIMPAR. A entrada do modelo sã o os controles e a sa´ıda, o Valor Presente L´ıquido (VPL). O modelo de baixa fidelidade é uma simula¸cã o do modelo de ordem reduzida TPWL/POD que foi previamente treinado utilizando-se uma sequência aleatória de controles. Este modelo pode ser retreinado ao longo do processo da maneira como descrito no cap´ıtulo 3. Portanto, foram constru´ıdas duas fun¸co˜es de interface do otimizador com os simuladores (SIMPAR e TPWL). Como o algoritmo foi concebido de maneira genérica, as variáveis de projeto utilizadas e a fun¸cão-objetivo devem ser normalizadas. Logo, estas fun¸co˜es fazem a adequa¸cão das variáveis aos limites determinados pelo usuário para o espa¸co de projeto e à ordem de grandeza do VPL.

4.6

Desempenho do algoritmo MSAO Com o objetivo de avaliar-se o algoritmo de otimiza¸cão MSAO desenvolvido neste

cap´ıtulo foram propostos alguns testes utilizando-se os modelos de reservatório apresentados


101

no Cap´ıtulo 2 com algumas modifica¸co˜es. Os modelos s˜ ao: SPE10 - Modificado, que é um corte de 24000 células do modelo SPE10 e Brugge - Modificado, que é uma adapta¸caõ ao SIMPAR do modelo de Brugge (benchmark apresentado em (Peters et al., 2010)) na qual a malha é convertida para cartesiana. A fun¸caõ-ob jetivo aqui considerada para ambos os casos é o Valor Presente L´ıquido - VPL com os mesmos parâmetros econômicos utilizados na avalia¸cão do TPWL/POD no Cap´ıtulo 3. O pre¸co do óleo é de 100 $/BBL e os custos de inje¸caõ e produ¸cão de água são de 10 $/BBL.

4.6.1 Resultados com modelo SPE10 - Modificado O primeiro problema que foi atacado baseou-se no modelo SPE10 exposto no Cap´ıtulo 2 considerando as densidades diferentes e fluido e rocha compress´ıveis como no u ´ ltimo caso estudado no Cap´ıtulo 3. O tempo de concessão do reservatório é de 3000 dias e o controle é feito sobre a pressão de fundo dos 4 po¸cos produtores em 4 ciclos de controle de 750 dias de dura¸caõ cada um. Desta forma o problema de otimiza¸cão possui 16 variáveis de projeto. A sequência de controles de treinamento foi determinada aleatoriamente da mesma forma que no teste do Cap´ıtulo 3, variando-se os controles dos quatro po¸cos entre 2000 psi e 6000 psi a cada 200 dias, isto é, considerando-se 15 ciclos de controles. A justificativa para isto é que se espera uma maior variabilidade dos estados gravados, o que incrementa a capacidade de aproxima¸caõ do TPWL/POD. A quantidade de autovalores foi mantida idêntica aos testes ´ importante ressaltar que a dura¸cão máxima dos passos de acurácia: υ p = 128 e υz = 256. E de tempo neste problema foi de 30 dias como no estudo de acurácia. O simulador decidiu cortá-los gerando estados suficientes para utilizar tamanha quantidade de autovalores. A quantidade de pontos avaliados no interior da primeira região de confian¸ca foi de 2N var +1 = 33, onde N var = 16 é o n´ umero de variáveis. Nesta quantidade já está inclu´ıdo o centro da primeira região. A justificativa para esta escolha é que só com esta quantidade de pontos seria poss´ıvel o ajuste de um polinômio quadrático sem os termos cruzados. O computador utilizado neste processo de otimiza¸caõ possui 12 processadores, sendo que um deles fica dedicado ao proceso do Matlab. Sendo assim, somente N proc = 11 é o n´ umero de processadores livres para simula¸caõ do SIMPAR. Em todas as avalia¸co˜es do modelo de alta fidelidade, são efetuadas N proc simula¸co˜es simultaneamente.


102

Desta forma, a quantidade de pontos avaliados em cada região de confian¸ca a partir da segunda itera¸cão foi fixada em N proc , o que permite que todas sejam feitas simultaneamente. A quantidade de avalia¸co˜es realizadas em conjunto com o candidato a ótimo foi fixada em N proc

− 1 pelo mesmo motivo.

Foram avaliadas duas possibilidades para o grau do polinômio interpolador utilizado na krigagem. O processo de otimiza¸caõ foi executado utilizando-se polinômio de grau zero (krigagem ordinária) e grau um (polinômio linear). O estudo comparativo entre os modelos de alta e baixa fidelidade realizado no Cap´ıtulo 3, indica que uma corre¸caõ linear pode ser mais eficiente. Também foram avaliados os parâmetros de controle de retreinamento. O valor limite para E T P W L (distância entre estados gravados no treinamento e os calculados na simula¸cão TPWL/POD) foi alterado de forma a aumentar ou diminuir o número de retreinamentos e é representado por E max . Note que quanto menor este limiar, maior será a tendência ao retreinamento. Também avaliou-se o limiar m´ınimo do parâmetro de qualidade do metamodelo ρ f para o qual deve-se cogitar o retreinamento. Este limiar é designado por ρmin . Como nos casos simulados neste estudo nunca se observou ρ f > 1, utilizou-se ρ min = 2 como uma maneira de obrigar o retreinamento o máximo n´ umero de vezes. O valor inicial para as pressões de fundo em todos os po¸cos e todos os ciclos de controle foi de 4000 psi, podendo variar de 2000 psi a 6000 psi que s˜ ao os limites aplicados durante o treinamento. O tamanho inicial da região de confian¸ca foi de 20 % do espa¸co de projeto baseado em (Horowitz et al., 2013). Além de aplicar-se o algoritmo MSAO variando-se alguns de seus parâmetros de configura¸caõ, foi realizada a otimiza¸caõ do VPL através da aplica¸caõ direta do algoritmo SQP (Programa¸cão Quadrática Sequencial) sobre o modelo de alta fidelidade (simulador SIMPAR) e sobre o modelo de baixa fidelidade (TPWL/POD). No último caso, apesar da otimiza¸caõ se realizar sobre o modelo reduzido, o resultado ao final é medido utilizando-se o simulador. Este é o método mais investigado na literatura para o uso de modelos reduzidos em otimiza¸cão e necessita de modelos substitutos quase perfeitos. No caso aqui estudado, onde se consideram não-linearidades bastante grandes, o TPWL/POD não é um substituto perfeito e portanto tal estratégia não parece ser aplicável. Os resultados obtidos estão apresentados na Tabela 4.1, onde SQP(HF) significa a aplica¸caõ do SQP diretamente sobre o SIMPAR, SQP(LF) significa aplicar o SQP sobre o TPWL/POD e medir o resultado final do VPL com o SIMPAR e MSAO significa a aplica¸caõ


103

do algoritmo MSAO proposto neste cap´ıtulo. Nesta tabela E max representa a distância máxima entre os estados gravados e calculados no TPWL/POD para que o retreinamento ocorra no caso de ρ f > ρ min . Também foram anotados os n´ umeros de itera¸co˜es (# It.), de avalia¸co˜es da fun¸cão de alta fidelidade, o SIMPAR (# HF), de avalia¸cões da fun¸caõ de baixa fidelidade, o TPWL/POD (# LF) e de retreinamentos. cão Tab. 4.1: Otimiza¸

do VPL utilizando diferentes esquemas - SPE10 com 4 ciclos de controle

Algorithm Polin. E max

ρmin

# It. # HF # LF #Retrein. NPV (108 $) ∗

SQP(HF)

-

-

-

5

88

-

-

6, 9229

SQP(LF)

-

-

-

11

-

187

-

6, 4826

MSAO

0

10

2

0,00

7

164

676

0

6,5424

MSAO

0

10

2

0,25

7

165

777

1

6,5318

MSAO

0

10

2

2,00

8

186

1027

2

6,5626

MSAO

1

10

2

0,25

5

120

715

0

7,2000

MSAO

1

10

2

2,00

11

254

1446

2

7,2238

MSAO

1

10

3

2,00

11

254

2040

5

7,2238

−

−

−

−

−

−

´ importante recordar que, conforme análise do speedup realizada no Cap´ıtulo E 3, uma avalia¸caõ de fun¸cão de alta fidelidade (SIMPAR) gasta 400 vezes o tempo de uma de baixa fidelidade (TPWL/POD). Esta propor¸caõ diminui para 300 após o retreinamento. Portanto, um alto número de avalia¸co˜es de baixa fidelidade não impacta tanto no tempo total do processo de otimiza¸cão. A referência da compara¸caõ de resultados é o algoritmo SQP acoplado diretamente ao simulador SIMPAR que está representado na primeira linha da tabela. Este é um algoritmo de otimiza¸caõ restrita de aplica¸cão geral comumente utilizado. Precisa dos gradientes para funcionar, que neste caso são calculados numericamente através de perturba¸co˜es nas variáveis e utilizando diferen¸cas finitas. Isto aumenta bastante o número de simula¸co˜es. O resultado utilizando-se tal algoritmo está na primeira linha da tabela. Nota-se que a simples substitui¸caõ do modelo de alta fidelidade pelo reduzido pode não levar a um bom resultado na otimiza¸cão. A segunda linha da tabela contém o resultado onde a fun¸cão-objetivo acoplada ao SQP foi o VPL calculado pelo TPWL/POD. Ao final do processo de otimiza¸cão calcula-se o VPL com o simulador. Este foi o pior resultado entre


104

todos, o que demonstra que o modelo de ordem reduzida precisa de corre¸cã o ao longo do processo de otimiza¸cão. Além disso, observa-se, ao utilizar-se uma krigagem ordin´ aria (polinômio de grau 0) consegue-se melhorar o resultado da otimiza¸cão. Entretanto, estes ainda n˜ ao foram os melhores resultados. Por outro lado, retomando-se a análise de acurácia do TPWL/POD apresentada no Cap´ıtulo 3, observa-se que o erro deste modelo reduzido não é constante para qualquer sequência de controles e depende muito do primeiro controle. Portanto, o uso de um polinômio regressor constante pode realmente não ser u ´ til, já que a parte estocástica do modelo kriging tem influência local. Esta foi a justificativa para a utiliza¸caõ de polinômios lineares. Observa-se na linha 6 da tabela de resultados que com o uso de um interpolador linear na krigagem, não houve necessidade de retreinamento pois o parâmetro de qualidade ρf do SAO se manteve superior a 0,25. Melhor ainda, mesmo sem retreinamento, o resultado conseguiu superar aquele obtido com o SQP. Isso é explicado pelo fato de este último algoritmo sofrer com certas descontinuidades do modelo de simula¸cão que são mitigadas pelas aproxima¸co˜es e interpola¸co˜es realizadas pela krigagem. Al´ em disso, observa-se uma ligeira melhora no resultado quando se aumenta o valor de ρf a partir do qual se avalia a necessidade de retreinamento. Na linha 7 tem-se um resultado um pouco maior onde foram realizados dois retreinamentos. A necessidade de retreinamento é avaliada pela distância entre os estados gravados e simulados E tpwl (calculada conforme Equa¸cão 3.34) que deve estar abaixo de 10

2

−

neste caso.

No entanto, quando este patamar máximo de erro é reduzido em uma ordem de grandeza, obtem-se o resultado da linha 8. Note que não h´ a melhora alguma apesar do aumento expressivo do n´ umero de retreinamentos. Tudo isto demonstra que retreinar nem sempre é indicado. Primeiramente, sabese que a qualidade do metamodelo medida por ρf fica comprometida na medida em que o processo se aproxima do ótimo devido a pequenas varia¸co˜es da fun¸cão-objetivo, fazendo com que um erro de aproxima¸cão, mesmo que pequeno, leve a valores baixos ou at´ e negativos para ρ f . Este problema tende a ser mitigado pela aplica¸caõ do critério sobre E tpwl que pode avaliar se a “culpa” pelos erros é da má qualidade do modelo reduzido. A Figura 4.4 mostra a evolu¸cã o de ρf , E tpwl e do valor ótimo ao longo das itera¸co˜es do MSAO nos casos com


105

polinômio interpolador linear. 0.07 1

0.06

0.8

0.05

0.6

L0.04 W P T

f

ρ

E

0.4

0.03

0.2

0.02

0

0.01

−0.2

2

4 6 # Iteração do MSAO

8

10

(a) Parâmetro de qualidade do MSAO ( ρf )

0 0

2

4 6 8 # Iteração do MSAO

10

12

(b) Indicador de erro do TPWL/POD ( E tpwl )

7.6 7.4 7.2 ) $ M M (

7

L P6.8 V

6.6 ro=0.25,E=1E−2 ro=2.00,E=1E−2 ro=2.00,E=1E−3

6.4 6.2 0

2


10

12

(c) Valor presente l´ıquido - VPL(10 8 $)

Fig. 4.4: Evolu¸caõ do processo de otimiza¸cão nas itera¸cões do algoritmo - Modelo SPE10 Modificado Observe primeiramente as curvas em vermelho que representam o primeiro resultado onde não houve retreinamento. O valor ótimo obtido para o VPL é ligeiramente menor do que os outros casos e isto se deve ao retreinamento que ocorreu na quarta itera¸cão. Esta atualiza¸caõ do modelo reduzido levou a um VPL um pouco menor, mas, que devido à melhor aproxima¸caõ, veio a ficar melhor na próxima itera¸caõ. Note também que o valor de E tpwl cai a um patamar de 10

2

−

depois do segundo

retreinamento na itera¸cão 6 (curva azul). Mesmo assim, ρf é próximo de zero (em algumas itera¸cões, negativo), o que indica que o metamodelo não é razoável. Entretanto, a culpa não


106

é do TPWL/POD já que E tpwl se mantém controlado. Isto se deve à proximidade do ótimo, onde as varia¸cões da fun¸caõ-objetivo são da ordem do erro da aproxima¸caõ. Pode-se confirmar esta ideia ao “apertar-se” o crit´ erio do erro para E tpwl

≤ 10

3

−

como nas curvas em preto. Não há nenhum ganho na fun¸caõ-objetivo em rela¸cão a curva azul. Tamb´ em não há melhora de ρf e observa-se um aumento de E tpwl apó s um dos retreinamentos ocorrerem. Isto demonstra que fazer E tpwl

≤ 10

3

−

é desnecessário e só leva a

retreinamentos excessivos, que podem vir a piorar o metamodelo. Na prática, o valor ideal para este parâmetro dependerá da experiência do profissional responsável pela otimiza¸cão. 4

16000 ) d / l b b ( o d i z u d o r p l a t o t o e l ó e d o ã z a V

3.5

Base Ótimo

14000

) d / l b 3 b ( a d i 2.5 z u d o r p 2 l a t o t

12000 10000 8000

1.5

a u g á e 1 d o ã z0.5 a V

6000 4000 2000 0

x 10

500

1000 1500 2000 Tempo (Dias)

2500

3000

0 0

(a) Vazão de óleo total

500

1000 1500 2000 Tempo (Dias)

2500

3000

(b) Vazão de água produzida total

4

3.5

x 10

) d / l b 3 b ( a d a2.5 t e j n i l a t o 2 t a u g á1.5 e d o ã z 1 a V

0.5 0

500

1000 1500 2000 Tempo (Dias)

2500

3000

(c) Vazão de água injetada total

Fig. 4.5: Vazões totais do campo obtidas no processo de otimiza¸cão — SPE10 - Modificado Analisando-se os resultados do melhor caso (Retreino com ρf

≤ 2,00 e E ≤ tpwl

10 2 ), obtém-se as vaz˜ oes totais do campo apresentadas na Figura 4.5. Nesta figura, os −


107

gráficos em vermelho representam as vazões para configura¸cão de controles inicial do processo de otimiza¸cão, os gráficos em preto, a configura¸caõ o´tima e os azuis, são as itera¸cões intermediárias. Nota-se destes resultados que a estrat´ egia do otimizador consiste em antecipar produ¸caõ de óleo e reduzir a produ¸caõ e inje¸cã o de água drasticamente em quase todo o tempo simulado. Também pode-se perceber que há grandes ganhos nas primeiras quatro itera¸cões do algoritmo, o que cessa nas itera¸cões seguintes. Os controles de pressão de fundo otimizados para os quatro po¸cos produtores estão representados na Figura 4.6 onde se observa que a decisão tomada pelo otimizador é de aumentar a press˜ a o de fundo nos u ´ ltimos ciclos de controle o que reduz a produ¸caõ de água, já que os po¸cos já produzem muita água nesta época.

6000 5500 ) i s5000 p ( o d4500 n u F 4000 e d o ã3500 s s e r P3000

W1 W2 W3 W4

2500 2000 0

500

1000 1500 2000 Tempo (Dias)

2500

3000

Fig. 4.6: Pressões de fundo dos po¸cos produtores otimizadas - Modelo SPE10 - Modificado Finalmente é importante avaliar a qualidade da aproxima¸caõ do metamodelo corrigido e do modelo reduzido (TPWL/POD). Para tanto, com o mesmo resultado utilizado na Figura 4.5, foram constru´ıdos, para cada itera¸caõ do algoritmo, crossplots com os pontos avaliados durante a itera¸caõ, onde o eixo x representa o resultado (VPL) obtido com o modelo de alta fidelidade (SIMPAR) e o eixo y, o resultado (VPL) em baixa fidelidade. Quanto melhor a aproxima¸caõ, mais próximos da reta de 45o estarão os pontos. As Figuras 4.7 a 4.9 apresentam as nuvens de pontos para as três primeiras itera¸cões. Os pontos apresentados são somente os considerados em cada uma das itera¸co˜es.


108

7.5

7.5

7

7

6.5

6.5

6 6

6.5

7

7.5

6 6

(a) Simpar x TPWL/POD

6.5

7

7.5

(b) Simpar x TPWL/POD+kriging

Fig. 4.7: Compara¸caõ do VPL obtido durante o processo de otimiza¸caõ - 1a Itera¸caõ - SPE10 - Modificado 7.6

7.5

7.4 7.2 7

7 6.8 6.6

6.5

6.4 6.2 6 6

6.5

7


7.5

6 6

6.5

7

7.5


Fig. 4.8: Compara¸caõ do VPL obtido durante o processo de otimiza¸caõ - 2a Itera¸caõ - SPE10 - Modificado Observando-se os crossplots da direita (compara¸cão entre VPL calculado pelo SIMPAR e pelo metamodelo com corre¸cão) em todos as figuras, observam-se dois grupos de pontos. O primeiro grupo recai exatamente sobre a reta, pois são os pontos considerados para a constru¸cão do modelo de krigagem e portanto são perfeitamente corrigidos. O segundo grupo representa o candidato a ótimo da itera¸caõ e aqueles que são avaliados dentro da candidata a próxima região de confian¸ca. Note que estes ultimos não recaem sobre a reta pois não foram utilizados na corre¸caõ do kriging . O valor do VPL obtido no pela simula¸caõ no centro da região de confian¸ca é marcado com um x vermelho e o candidato a ótimo obtido


109

nesta região é marcado com um diamante vermelho. Nestas itera¸co˜es a corre¸caõ por krigagem se mostrou bastante efetiva, melhorando razoavelmente os resultados da aproxima¸caõ (comparem-se os crossplots com os resultados do TPWL/POD e do modelo reduzido com a corre¸cão por krigagem). Isto pode também ser medido pelo valor de ρf próximo de 1, como mostra a Figura 4.4(a). Entretanto, desta mesma Figura 4.4(a), observa-se que, a partir da 5 a itera¸caõ, o metamodelo mostra resultados piores com a queda drástica do ρ f , o que pode ser devido à proximidade do ótimo. 7.5

7.5

7

7

6.5

6.5

6 6

6.5

7


7.5

6 6

6.5

7

7.5


Fig. 4.9: Compara¸caõ do VPL obtido durante o processo de otimiza¸caõ - 3a Itera¸caõ - SPE10 - Modificado

4.6.2 Resultados com modelo de Brugge - Modificado O segundo problema atacado baseou-se no modelo Brugge - Modificado exposto no Cap´ıtulo 2 considerando as densidades da a´gua e óleo diferentes e fluido e rocha compress´ıveis conforme sua descri¸caõ em (Peters et al., 2010). O tempo de concessão do reservatório é de 3600 dias e o controle é feito sobre a pressão de fundo dos 20 po¸cos produtores, variando de 40Kgf /cm2 a 140 Kgf /cm2 , em 3 ciclos de controle de 1200 dias de dura¸cão cada um. Os 10 po¸cos injetores são mantidos com 300 Kgf /cm2 constante de pressão de fundo. Desta forma o problema de otimiza¸cão possui 60 variáveis de projeto, transformando-o em um problema razoavelmente grande. A sequência de controles de treinamento foi determinada aleatoriamente trocandose os controles dos 20 po¸cos independentemente a cada 200 dias, isto é, considerando-se 18


110

ciclos de controles. O valor destes controles tem os mesmos limites da defini¸caõ do problema de otimiza¸ca˜ o: 40Kgf /cm2 a 140 Kgf /cm2 . A quantidade de autovalores foi υ p = 80 e υz = 100. Esta quantidade de autovalores foi determinada através de alguns testes simples. Entretanto, para gerar número de estados suficiente para isto, utilizou-se uma dura¸cã o máxima para o passo de tempo de 10 dias. Os u ´ltimos valores singulares considerados estão entre a ordem de 10

1

−

e 10 2 . A −

figura 4.10 apresenta um gráfico dos valores singulares para a SVD da matriz de snapshots de pressões (3.4(a)) e de fra¸co˜es mássicas (3.4(b)) para a simula¸caõ de treinamento dispostos em ordem decrescente com a indica¸caõ dos respectivos pontos de corte escolhidos. 4

4

10

10

2

2

10

10

0

0

10

10

−2

−2

10

10

−4

10

0

−4

50

100

(a) Pressão

150

200

10

0

50

100

150

200

(b) Fra¸cão mássica

Fig. 4.10: Valores singulares dispostos em ordem decrescente - Modelo Brugge-Modificado O número de pontos avaliados no interior da primeira região de confian¸ca foi de 2N var + 1 = 121, onde N var = 60 é o número de variáveis. Nesta quantidade já está inclu´ıdo o centro da primeira regiã o. O polinˆ omio interpolador utilizado na krigagem é de primeiro grau. Os valores dos parâmetros de controle do retreinamento foram fixados: E max = 10

2

−

e ρmin = 0,25. O valor inicial para todas as press˜ oes de fundo controladas, centro da

primeira região de confian¸ca, foi fixado no meio da escala, 90 kgf /cm2 . O número de processadores livres para simula¸cão do SIMPAR é N proc = 11 e o número de simula¸co˜es executadas em paralelo é sempre igual a este número. Portanto, a quantidade de pontos simulados em cada região de confian¸ca para constru¸caõ do modelo


111

de corre¸cão a partir da segunda itera¸cão e a quantidade de pontos simulados em conjunto com o candidato a ótimo foram parâmetros avaliados, variando em múltiplos de N proc e são designadas pela palavra #Amostras. Tamb´ em foram avaliados diferentes valores para o raio da região de confian¸ca inicial, designado neste estudo por Rinit e dado por uma fra¸caõ das dimens˜ oes do espa¸co de projeto. A Tabela 4.2 reune os resultados obtidos neste estudo. Note que a referência para compara¸caõ está na primeira linha e consiste na aplica¸cão do algoritmo SQP diretamente ao modelo de alta fidelidade. Entretanto, por necessidade de economia de tempo, só se permitiu que tal algoritmo realizasse duas itera¸cões, não tendo neste momento atingido a convergência. Tab. 4.2: Otimiza¸ cão

do VPL em diferentes esquemas - Brugge-Modificado com 3 ciclos de controle

Algorithm #Samples R Init. # It. # HF # LF #Retrein. NPV (109 $) ∗

SQP(HF)

-

-

2

122

-

-

28,01

MSAO

11

0,25

9

298

4099

1

26,79

MSAO

33

0,25

7

518

3467

1

26,35

MSAO

11

0,50

11

342

5942

3

27,31

MSAO

33

1,00

13

913

9661

2

27,54

Com este modelo de Brugge houve uma redu¸caõ do speedup do TPWL/POD em rela¸caõ ao SIMPAR se comparado ao speedup obtido com o modelo SPE10. O computador em que estes estudos foram realizados é o mesmo utilizado nos estudos do Cap´ıtulo 3 e consiste em um Intel Xeon X5680 3,33 GHz, com 24 GB de RAM, 12 n´ ucleos e sistema operacional Windows 7. Enquanto uma simula¸cão de alta fidelidade custou em m´ edia 1700 s, uma de baixa fidelidade custou 10 s, levando a um speedup de 170, menor do que o do caso SPE10Modificado, apresentado no Cap´ıtulo 3. O tempo da simula¸cão de treinamento foi de 1600 s, o tempo para o cálculo da base foi 28 s e para o cálculo das derivadas reduzidas, 234 s, custando o treinamento um total de 1,2 simula¸cões. Além disso, o tempo de retreinamento foi de 42s para constru¸caõ da base e 470s para cálculo de todas as derivadas reduzidas considerando a nova base. Para entender o efeito dos retreinamentos do TPWL/POD no processo de otimiza¸caõ, foi feita uma análise similar à realizada com o primeiro caso SPE10-Modificado. A Figura 4.11 mostra a evolu¸caõ de ρf , E tpwl e do valor ótimo ao longo das itera¸cões do


112

MSAO nos dois melhores casos apresentados nas duas u ´ ltimas linhas da Tabela 4.2. Note que em ambos os casos houve redu¸caõ do erro E tpwl , mas que o parâmetro de qualidade do metamodelo ρf não se mostrou muito satisfatório. No caso com 33 amostras no interior da região de confian¸ca houve momentos no final da otimiza¸caõ em que ρ f ficou próximo de 1,0, mas em geral seu valor se mostrou baixo. Entretanto, como E tpwl pôde ser controlado com retreinamento, conclui-se que neste caso, a “culpa” pelo problema é do kriging . Na verdade, utilizaram-se poucas amostras no interior da regi˜ ao para o n´ umero de variáveis considerado. Portanto, houve a necessidade de reutiliza¸caõ de um grande número de pontos das regiões anteriores, o que pode causar problemas na aproxima¸caõ por krigagem. 1.2

0.16

1

0.14 0.12

0.8

0.1 0.6

l w p t

f

ρ

E

0.4

0.08 0.06

0.2

0.04

0 −0.2 0

0.02 2


10

12

(a) Parâmetro de qualidade do MSAO ( ρf )

0 0

2

4 6 8 10 # Iteração do MSAO

12

(b) Indicador de erro do TPWL/POD ( E tpwl )

2.8 2.75 2.7

) $ M M ( 2.65 L P V

2.6

2.55

11 amostras, ρ = 0.5 0

33 amostras, ρ = 1.0 2.5 0

0

2

4 6 8 10 # Iteração do MSAO

12

(c) Valor presente l´ıquido - VPL(10 8 $)

Fig. 4.11: Evolu¸caõ do processo de otimiza¸caõ nas itera¸cões do algoritmo - Modelo SPE10 - Modificado


113

5

8

5

x 10

Base Ótimo

) d / l b b (

7 6

o d i z u d o r p l a t o t

5 4

o e l ó e d o ã z a V

3 2 1 0 0

1000

2000 Tempo (Dias)

4.5

x 10

) d / 4 l b b ( a3.5 d i z u 3 d o r p2.5 a u g á 2 e d l 1.5 a t o t o 1 ã z a0.5 V

0 0

3000

(a) Vazão de óleo total

1000

2000 Tempo (Dias)

3000

(b) Vazão de água produzida total

6

2.5

x 10

) d / l b b ( 2 a d a t e j n1.5 i l a t o t a u 1 g á e d o ã0.5 z a V

0 0

1000

2000 Tempo (Dias)

3000

(c) Vazão de água injetada total

Fig. 4.12: Vazões totais do campo obtidas no processo de otimiza¸caõ — Brugge - Modificado Como u ´ ltima análise, observe-se a estrat´ egia do otimizador para maximizar o VPL. A Figura 4.12 mostra em vermelho as vazões totais iniciais (caso base), em preto as otimizadas e em azul, as itera¸cões intermediárias. Assim como no caso anterior do SPE10Modificado, a estrat´ egia consistiu em antecipa¸ca˜ o de óleo e redu¸cão da produ¸caõ e inje¸caõ de água. Entretanto, a diferen¸ca entre caso base e otimizado só foi expressiva na produ¸caõ de água.

˜ EM 5. APLICAC ¸ AO ˜ DE INCERTEZAS PROPAGAC ¸ AO Não há como se pensar em gerenciar um campo de petróleo sem levar em conta as incertezas envolvidas. Estas podem ser de natureza geológica, por não se saber exatamente as caracter´ısticas de subsuperf´ıcie; de natureza operacional, por não se ter pleno controle ou não se confiar totalmente em um equipamento; econômica, devido às flutua¸cões do câmbio, dos custos dos equipamentos, do pre¸co do petróleo, etc. Entretanto, levar-se em conta as incertezas no processo de gerenciamento do reservatório, aplicando técnicas estat´ısticas capazes de medir os riscos, é, por si só, muito custoso computacionalmente por causa da necessidade de lidar com muitos modelos que representam os diferentes cenários incertos e simula¸co˜es cada vez mais caras. ´ portanto de grande utilidade a propaga¸caõ da incerteza nos parâmetros do E modelo de simula¸caõ (mapas geológicos, controles dos po¸cos, capacidade dos equipamentos, pre¸co do petróleo, etc) para a sua sa´ıda (e.g. Valor Presente L´ıquido - VPL). Isto demanda a execu¸cão de muitas simula¸co˜es do modelo considerando múltiplos cenários. Além disso, a tendência atual de criar modelos cada vez mais detalhados e complexos aumenta ainda mais o custo dos estudos de incerteza. Isso fica ainda pior se a propaga¸cão de incertezas está combinada a um outro processo, como por exemplo otimiza¸caõ ou ajuste de histórico sob incerteza. Este cap´ıtulo descreve uma técnica que tem por objetivo baixar o custo computacional da propaga¸caõ de incerteza em um modelo de simula¸caõ de alta fidelidade pela utiliza¸caõ de um metamodelo que consegue prever com um boa acurácia os resultados da sa´ıda de alta fidelidade. O método aqui utilizado est´ a descrito em (Ng et al., 2013) sendo que o modelo de alta fidelidade considerado será o simulador de reservatórios SIMPAR e o modelo de baixa fidelidade, o TPWL. A ideia é propagar a incerteza nos controles (no caso pressões de fundo de po¸co) para a sa´ıda do modelo que neste caso é o valor presente l´ıquido 114

5.1. Método das Vari´ aveis de Controle

115

obtido no fim do tempo de concessão. ´ importante ressaltar que a incerteza nos controles é pequena se comparada aos E parâmetros geológicos. Entretanto, ainda não está dispon´ıvel na implementa¸caõ do TPWL aqui implementado, a considera¸caõ de tais parâmetros. Um esfor¸co neste sentido é encontrado em (He et al., 2013). Al´ em disso, não é de todo artificial considerar incertezas nos controles, uma vez que existem erros grandes nas medidas de pressão de fundo (as medi¸cões são normalmente indiretas) e além disso existem grandes dificuldades na aplica¸caõ destes controles, o que traz um certo n´ıvel de desconhecimento sobre quais são as verdadeiras pressões de fundo aplicadas num dado instante de tempo. O cap´ıtulo está dividido em duas partes. A primeira descreve em detalhes a matemática da técnica e o algoritmo para a propaga¸cão da incerteza na média e no desvio padr˜ ao. A segunda parte apresenta resultados obtidos com o problema de simula¸ca˜ o de reservatórios controlado pela pressão de fundo de po¸co.

5.1

M´ etodo das Vari´ aveis de Controle Considere o problema de estima¸cão do valor esperado de uma variável aleatória

X , x ref = E [X ]. Tome um conjunto de amostras independentes e identicamente distribu´ıdas de X , x1 , x2 ,

··· , x . O estimador Monte Carlo para x

e ref ´

n

descrito na Equa¸caõ 5.1 e dado

pela média aritmética dos valores amostrados.

1 x¯n = n

n

�

xi

(5.1)

i=1

Considerando σ x2 a variância da variável aleatória X , pode-se calcular a variância no estimador de x ref , conforme Equa¸cão 5.2 como medida do erro para este estimador. 1 V ar[¯xn ] = 2 V ar n

�� n

i=1

σx2 xi = n

(5.2)

A ideia da técnica aqui descrita é utilizar a correla¸cão entre X (que representará o modelo de alta fidelidade) e uma variável aleatória auxiliar Y (que representará o modelo de baixa fidelidade) de forma a reduzir o número de amostras de X necessárias para atingir uma variância aceitável para seu estimador.

5.2. Monte Carlo Multifidelidade

116

Sendo assim, sejam y 1 , y2 ,

··· , y um conjunto de amostras independentes e idenn

ticamente distribu´ıdas da variável aleatória auxiliar Y cuja média é conhecida yref = E [Y ]. Como

− Y ] = 0 o estimador de variáveis de controle (control variate estimator ) para

E [yref

xref é definido como na Equa¸caõ 5.3, onde β é um parâmetro de controle.

− y¯ )

¯n + β (yref xˆn = x

(5.3)

n

A variância deste estimador, ap´ os algumas simplifica¸co˜es, está expressa na Equa¸caõ 5.4.

V ar[ˆxn ] =

1 2 (σX + β 2 σY 2 + 2βρ XY σX σY ) n

(5.4)

Na Equa¸cão 5.4, ρXY é o coeficiente de correla¸cão entre X e Y . O que se quer é o estimador com menor variância e, ao minimizar-se V ar[ˆxn ] com rela¸caõ a β obtem-se o valores ótimos β e V ar[ˆxn ] das Equa¸cões 5.5 e 5.6. ∗

∗

β = ρ xy ∗

V ar[ˆ xn ] = ∗

1 (1 n

σX σY

−ρ

2 2 XY )σX .

(5.5)

(5.6)

Desta forma observa-se uma redu¸caõ da variância do estimador para xref se Y ´ nisso que se baseia o m´ for bastante correlacionado a X . E etodo de Monte Carlo com multifidelidade que é utilizado neste cap´ıtulo.

5.2

Monte Carlo Multifidelidade Daqui por diante as variáveis aleatórias X e Y serão substitu´ıdas respectivamente

por avalia¸co˜es de um modelo de alta e um de baixa fidelidade. Desta forma, seja a fun¸caõ hhigh (z) um resultado da simula¸caõ do modelo de alta fidelidade que toma como entrada o vetor z e retorna a variável x como sa´ıda. No problema da simula¸ cão de reservatórios aqui considerado, isto representa o VPL calculado pelo SIMPAR considerando-se uma sequência de BHP’s para todos os po¸c os. Quando se fala de propaga¸caõ de incertezas, o vetor de entrada representa um vetor aleatório Z e o objetivo é estimar x ref = E [hhigh (z)].


117

Por outro lado, seja a fun¸caõ hlow (z) um resultado da simula¸caõ do modelo de baixa fidelidade que toma como entrada o mesmo vetor z aplicado ao de alta fidelidade, mas retorna a variável y como sa´ıda. No problema em quest˜ ao, representaria o VPL calculado através da simula¸caõ do modelo de ordem reduzida TPWL/POD treinado anteriormente, considerando uma sequência de pressões de fundo para os po¸cos. Este modelo é u ´til pois é altamente correlacionado ao de alta fidelidade. O procedimento aqui aplicado, descrito em detalhes em (Ng et al., 2013), é resultado da adapta¸caõ da técnica de variáveis de controle a um contexto onde existem modelos de múltiplas fidelidades x i = h high(zi ) e y i = hlow (zi ), sendo z i amostras i.i.d. da vetor aleatório

Z que representa os controles dos po¸cos neste caso. Entretanto, yref =

E [hlow (zi )]

reque-

rido pela Equa¸cã o 5.3 não é conhecido e será aproximado por uma média de m simula¸cões y¯m , considerando m >> n. Note que isto faz com que sejam realizadas m

− n simula¸cões

adicionais do modelo de baixa fidelidade. Assim, a Equa¸caõ 5.3 se transforma na Equa¸caõ 5.7.

xñ = x¯n + β (y¯m

− y¯ ) = β ¯y n

¯n m + ( x

− β ¯y )

(5.7)

n

A segunda igualdade mostra que o maior trabalho é realizado sobre o modelo de baixa fidelidade cujo resultado deve ser corrigido por um termo de corre¸cão β . A variˆ ancia deste estimador multifidelidade pode ser obtida pela Equa¸caõ 5.8.

V ar[˜xn ] =

1 2 2 (σX + β 2 σY + 2βρ XY σX σY ) n

− m1 (β σ

2 2 Y

+ 2βρ XY σX σY )

(5.8)

Observe que o processo de cálculo deste estimador multifidelidade x ˜ precisa de n simula¸cões de alta fidelidade e m de baixa. Portanto, além de determinar o valor ótimo para β é poss´ıvel distribuir o recurso computacional entre as simula¸cões de alta e baixa fidelidade de maneira ótima (Ng et al., 2013). Com o objetivo de comparar com o método de Monte Carlo tradicional define-se uma unidade de esfor¸co computacional p contando-se o n´ umero equivalente de simula¸cões de alta fidelidade como mostrado na Equa¸cão 5.9.

p

≡

� �

m r = n 1 + n + w w

(5.9)

Na Equa¸cão 5.9, w é a razão entre o tempo médio de computa¸caõ do modelo de alta fidelidade e o de baixa fidelidade. No caso do SIMPAR e do TPWL/POD, esta razão


118

é aproximadamente 400. Além disso, r = m/n é a razã o entre os números de simula¸co˜es de baixa e alta fidelidade. Pode-se re-escrever a variância do estimador da Equa¸caõ 5.8 em termos de p, r e w conforme Equa¸caõ 5.10.

� �� − �   

1 r 1+ V ar[˜xn ] = p w

1 r

2 + 1 σX

2 (β 2 σY + 2βρ XY σX σY )

�

(5.10)

1

n

Minimizando V ar[˜xn ] em rela¸cão a β e r e considerando fixo o recurso computacional p, obtem-se como solu¸caõ o´tima β e r como exposto na Equa¸caõ 5.11. ∗

β = ρ xy ∗

∗

σX , σY

r = ∗

�

wρ2XY 1 ρ2XY

(5.11)

−

A Equa¸cão 5.12 apresenta o valor ótimo para o erro, V ar[˜xn ] considerando β e ∗

∗

r . ∗

� �� − � − � �

1 r 1+ V ar[˜xn ] = p w

∗

∗

1

1

1 r

∗

2 ρ2XY σX .

(5.12)

Com esta solu¸cão, o que se percebe é que r aloca a maior parte do esfor¸co ∗

computacional para a simula¸caõ do modelo de baixa fidelidade se este é barato (w é grande) e/ou acurado (ρXY é próximo de 1). At´ e este ponto desenvolveu-se uma metodologia para estimar-se a esperan¸ca xref = ´ poss´ıvel aplicar a mesma simula¸cão Monte Carlo Multifidelidade de forma a E [hhigh (z)]. E 2 estimar outros parâmetros estat´ısticos como, por exemplo, V ar[X ] = σX = V ar[hhigh(z)],

desde que tomando-se os devidos cuidados para manter a eficiência do método. Sendo assim, como utilizar o que já foi desenvolvido para a média no cálculo da variância? Seria imediato utilizar-se a rela¸caõ V ar[X ] = E [X 2 ]

2

− (E [X ]) . Entretanto, apesar de simples, segundo (Ng

et al., 2013), esta abordagem é ingênua por produzir resultados pouco satisfat´ orios. Isto se deve ao fato de a correla¸caõ ρX 2 Y 2 se mostrar pior do que ρXY . Portanto é melhor utilizar o algoritmo de cálculo da variância em um passo e descrito em (Knuth, 1997). Para tanto, definam-se η i e ν i como exposto nas Equa¸cões 5.13 e 5.14.

5.3. Implementa¸caõ do algoritmo

ηi =

n

119

� · � · � · �

hhigh (zi )

n−1

ν i =

m m−1

·

−

−

hhigh (zi )

hlow (zi )

n

hhigh (z j )

j=1

1 n−1

1 m

−

−

hlow (zi )

1 n

∑ �· ∑ � ∑ �· ∑ � n−1

(5.13)

(5.14)

hhigh (z j )

j=1

m

hlow (z j )

j=1

1 m−1

m−1

hlow (z j )

j=1

Usando essa defini¸cão para as amostras ηi e ν i obtem-se que suas médias representam as variâncias desejadas conforme Equa¸cões 5.15 e 5.16.

V ar[hhigh (z)] = η¯n = =

V ar[hlow (z)] = ν¯ m = =

n

1 n

∑ � ∑ ∑ � ∑ ηi =

i=1

n

1

n−1

1 m

hhigh(zi )

i=1

−

1 n

n

∑

�

2

(5.15)

hhigh(z j )

j=1

m

ν i =

i=1

1 m−1

m

hlow (zi )

i=1

−

1 m

m

∑

j=1

�

2

hlow (z j )

(5.16)

Portanto, pode-se aplicar, a primeira vista, a mesma metodologia desenvolvida anteriormente para estimar a esperan¸ca de variáveis definidas em 5.13 e 5.14, o que significaria a sua variância. Segundo (Ng et al., 2013), o coeficiente de correla¸caõ entre estas variáveis se mostra mais alto do que se fosse aplicada a rela¸caõ com os quadrados das variáveis aleatórias. Entretanto, η i e ν i não são mais independentes, o que viola uma hip´ otese utilizada na teoria. Por outro lado, os resultados de (Ng et al., 2013) e os obtidos com o modelo reduzido aqui estudado mostram que essa abordagem ainda se mostra efetiva.

5.3

Implementa¸c˜ ao do algoritmo O algoritmo para a simula¸caõ Monte Carlo Multifidelidade é executado de forma

iterativa, incrementando-se o número de amostras de alta fidelidade e calculando-se os valores

5.3. Implementa¸caõ do algoritmo

120

de β , r a cada incremento. Na pr´ atica, como as variâncias σx e σy e o coeficiente de ñ e correla¸caõ ρ XY são desconhecidos, o melhor que se pode fazer é calcular os estimadores β ∗

∗

rñ baseados nas n primeiras amostras dos modelos de alta fidelidade — variável aleatória X — e baixa fidelidade — vari´ avel Y , conforme Equa¸co˜es 5.17, 5.18 e 5.19.

∑∑ − − − � − ∑ ∑ −∑−

ñ = β

n x¯n )(yi i=1 (xi n y¯n )2 i=1 (yi

ρñ =

[

−

(5.17)

(5.18) 2

n i=1 (xi

n i=1 (xi

wρ˜2n 1 ˜ ρ2n

rñ =

[

y¯n )

x¯n )(yi y¯n )] x¯n )2 ] [ ni=1 (yi y¯n )2 ]

−

(5.19)

O processo de simula¸cão Monte Carlo Multifidelidade precisa da escolha de alguns parâmetros: uma raz˜ ao média entre os tempos de computa¸caõ dos modelos de alta e baixa fidelidade, w, um n´ umero de amostras iniciais, n0 , um valor para o incremento no número de amostras, n∆ , um modelo de alta fidelidade hhigh(z), que neste caso é o simulador de reservatórios SIMPAR descrito no cap´ıtulo 2, um modelo de baixa fidelidade hlow (z), representado pelo modelo de ordem reduzida TPWL descrito no cap´ıtulo 3. Também é sorteada uma sequência pseudo-aleatória de vetores de entrada z i sobre a distribui¸cão de Z , que é a vari´ avel que possui incerteza e que no caso deste trabalho representa os valores de controle dos po¸cos, pressões de fundo em fluxo. O algoritmo como proposto em (Ng et al., 2013) e utilizado nesta disserta¸caõ pode ser descrito pela sequˆ encia de passos abaixo. 1. Sejam n old = 0, m old = 0, l = 0 e n = n 0 . 2. Avaliam-se as amostras de alta fidelidade (simula¸caõ do SIMPAR) para o estimador da média x i = h high (zi ), i = nold + 1, . . . , n. 3. Se l < n, avaliam-se as amostras de baixa fidelidade (simula¸caõ TPWL) para o estimador da média yi = hlow (zi ), i = l + 1, . . . , n e faz-se l = n. 4. Calculam-se as médias x¯n e y¯n , médias de alta e baixa fidelidade respectivamente. 5. Avaliam-se as amostras para o estimador de variˆ ancia ηi e ν i , i = n old +1, . . . , n conforme Equa¸co˜es 5.13 e 5.14. 6. Calculam-se η¯n e ν¯n , representando as variâncias, conforme Equa¸co˜es 5.15 e 5.16.

5.4. Resultados obtidos

121

˜xn , ρ˜2 , r˜xn usando xi , yi 7. Computam-se β xn

{

n i=1

}

˜ηn , ρ˜2 , r˜ηn utilizando ηi , ν i e β ηn

{

n i=1 .

}

8. Calcula-se o número de avalia¸cões de baixa fidelidade m = n max(˜ rxn , r˜ηn ). 9. Se l < m, avaliam-se as amostras de baixa fidelidade que faltam para o estimador da média y i = h low (zi ), i = l + 1, . . . , m. Faz-se l = m. 10. Calcula-se a média y¯m . 11. Avaliam-se as amostras de baixa fidelidade que faltam para o estimador da variˆ ancia ν i , i = l + 1, . . . , m de acordo com a Equa¸caõ 5.14. 12. Calcula-se a média ν¯ m. 13. Computam-se os estimadores multi-fidelidade xñ e ηñ para média e variância respectivamente utilizando a Equa¸cão 5.7. 14. Estima-se o erro V ar[˜xn ] pela equa¸cão 5.12. ∗

15. Se o erro é grande, faz-se n old = n, m old = m, n <

5.4

−n + n

∆ e

volta-se ao passo 2.

Resultados obtidos Esta disserta¸cão tem como u ´ltimo objetivo avaliar o método Monte Carlo Multifi-

delidade para propaga¸caõ de incertezas no contexto de simula¸caõ de reservatórios de petróleo utilizando como modelo de ordem reduzida o TPWL/POD desenvolvido no Cap´ıtulo 3 e como simulador o SIMPAR descrito no Cap´ıtulo 2. Os resultados aqui apresentados utilizaram o modelo SPE10 apresentado no Cap´ıtulo 2 tamb´ em utilizado em um processo de otimiza¸caõ no Cap´ıtulo 4. Foram utilizados duas variantes do modelo: a primeira desconsiderando as compressibilidades e com as densidades iguais para ambos os fluidos; a segunda considera as compressibilidades dos fluidos e rocha e as densidades diferentes. As incertezas aqui consideradas estão associadas a` pressão de fundo de po¸co. Estas não são as mais cr´ıticas no gerenciamento de reservatórios uma vez que a grande fonte de incertezas nesta área é a geologia. Entretanto, ainda não está dispon´ıvel (por não fazer parte do escopo deste trabalho) um modelo de ordem reduzida capaz de avaliar a influência de parâmetros geológicos como porosidade, permeabilidade, etc. no VPL (He (2010) - Apêndice 1 e He et al. (2013) cont´ em um exemplo TPWL considerando lineariza¸caõ em torno das


122

transmissibilidades em todas as células). Por outro lado, existem sim incertezas associadas à pressão de fundo de po¸co uma vez que os modelos de escoamento nã o são perfeitamente precisos. Al´ em disso, existem problemas de confiabilidade de equipamentos (e. g. bombas de fundo, gas lift ) que quando falham causam varia¸cões da pressão de fundo e têm influência na produ¸cão. Sendo assim, tais estudos possuem boa utilidade. Neste estudo, as press˜ o es de fundo dos po¸cos produtores são independentes e uniformemente distribu´ıdas entre 2000 psi e 6000 psi. Os po¸cos injetores ficam com pressão de fundo constante em 12000 psi. O tempo de concess˜ ao considerado é de 3000 dias. Os parâmetros para cá lculo do VPL são os mesmos utilizados nos Cap´ıtulos 3 e 4: taxa de retorno de 10 % a. a., pre¸co do óleo de $100 /BBL e custos de inje¸cão e produ¸cão de água de $10/BBL. O modelo de ordem reduzida foi treinado uma única vez para cada um dos casos considerados com uma sequência de 15 controles aleatórios da mesma forma que nos estudos de acurácia e de otimiza¸caõ dos Cap´ıtulos 3 e 4. Com rela¸cão a` configura¸caõ inicial do algoritmo Monte Carlo Multifidelidade, utilizou-se um n´ umero inicial de amostras n0 = 11 e um incremento do número de simula¸cões n∆ = 11 já que havia 11 processadores dispon´ıveis para a simula¸caõ do SIMPAR. Foram realizada 40 itera¸co˜es do algoritmo, considerando-se assim 440 simula¸co˜es de alta fidelidade. A razão entre os tempos de simula¸caõ de alta e baixa fidelidade w foi calculada pelo tempo obtido no treinamento e neste caso foi igual a aproximadamente 350 para os dois casos considerados. A Figura 5.1 mostra a compara¸caõ de resultados entre o método de Monte Carlo Tradicional e Multifidelidade para o caso incompress´ıvel com densidades iguais. Em cada gráfico o eixo x representa o esfor¸co computacional pn em simula¸co˜es de alta fidelidade equivalentes, aumentando na medida em que se aumenta o número de simula¸cões de alta fidelidade. O erro médio quadrático (Root Mean Square Error - RMSE) representado pela variˆ ancia do estimador da esperan¸ca nos dois casos é apresentado na Figura 5.1(a). Nota-se que esta cai mais rapidamente quando se utiliza o esquema multifidelidade. A estimativa para a média obtida com os dois métodos está representada na Figura 5.1(b). Observa-se a convergência entre os dois resultados por causa da interse¸cão entre os intervalos de confian¸ca (para um desvio padrão) representados por x e +.


123

−1

10

6.36 6.34 ) 6.32 $ M 6.3 M ( a d6.28 a m i t 6.26 s E a i d6.24 é M

E S −2 M10 R

6.22

MC Multifidelidade MC

−3

10

1

2

10

10 Esforço Computacional (p)

(a) RMSE (V ar(˜ xn ))

6.2 3

10

6.18 1 10

MC Multifidelidade MC 2

10 Simulações HF

3

10

(b) Valor Esperado x ñ

Fig. 5.1: Compara¸cão entre MC Multifidelidade e Tradicional - Reservatório Incompress´ıvel e Densidades dos Fluidos Iguais Os parâmetros de controle do Monte Carlo Multifidelidade estão apresentados na Figura 5.2 para as duas variáveis estimadas: média e variância. A evolu¸cão da correla¸caõ com o esfor¸co computacional está na Figura 5.2(a). Observa-se uma alta correla¸caõ entre os modelos de alta e baixa fidelidade (amostras xi e yi ), o que demonstra que o TPWL/POD é uma boa representa¸caõ aproximada para a simula¸caõ deste reservatório com o SIMPAR. Entretanto, a correla¸caõ para a estimativa da variância (amostras η i e ν i ) é menor. A razão r entre o número de simula¸cões de alta e baixa fidelidades é mais alta para a estimativa média do que para a variˆ ancia devido à correla¸caõ mais alta. O valor utilizado deve ser o máximo entre os dois e flutuou em torno de 40, conforme mostrado na Figura 5.2(b). Confirma-se também a redu¸cã o no parâmetro de corre¸caõ a medida que o

→ ∞, tem-se que

processo evolui já que, como se espera, quando o número de amostras n y¯n

→ y ≈ y¯ ref

m.

Por outro lado, os resultados mudam razoavelmente ao considerar-se o reservatório compress´ıvel e as densidades dos fluidos diferentes. Isto aumenta bastante a não linearidade do problema, interferindo na qualidade dos resultados. Além disso, conforme visto nos resultados de acurácia do TPWL/POD no Cap´ıtulo 3, o erro da aproxima¸caõ pode, nestes casos, depender muito da sequência de treinamento, levando ao aparecimento de algum viés.


124 65

1

Média Variância

Média Variância

60 55

0.95

50 0.9

45 r

ρ

40 0.85

35 30

0.8

25 0.75 1 10

2

10 Simulações HF

20 1 10

3

10

2

3

10 Simulações HF

(a) Correla¸cão (ρXY )

10

(b) Razão de no. de Sim’s (r )

0.9

Média Variância

0.8

0.7 β

0.6

0.5

0.4 1 10

2

3

10 Simulações HF

10

(c) Corre¸cão β

Fig. 5.2: Evolu¸caõ dos Parâmetros do MC Multifidelidade para Estima¸caõ da Média e da Variância - Reservatório Incompress´ıvel e Densidades dos Fluidos Iguais A Figura 5.3 mostra a compara¸caõ de resultados entre o método de Monte Carlo Tradicional e Multifidelidade para o caso compress´ıvel com densidades diferentes. Para obterem-se tais resultados, utilizou-se os mesmos valores de compressibilidades utilizados nos Cap´ıtulos 3 e 4: compressibilidade do o´leo co = 3 e da rocha co = 1

× 10

7

−

× 10

6

−

psi 1 , da água c w = 1 −

× 10

6

−

psi

1

−

psi 1 . A diferen¸ca nas densidades dos fluidos foi de 5 lbm/ft3 . −

Observa-se que o ganho na velocidade de convergência (redu¸caõ do erro) para a estimativa da média ao comparar-se com o Monte Carlo tradicional é menor do que no caso incompress´ıvel. Isto ocorre devido à menor correla¸cão entre o modelo de ordem reduzida e o real. A evolu¸caõ do erro com o custo computacional se apresenta na Figura 5.3(a).


125

Por outro lado, observa-se na Figura 5.3(b) que com este número de simula¸co˜es, as estimativas para a média calculadas pelos dois métodos ainda n˜ ao convergiram para um desvio padrão (observando-se os intervalos marcados po x e + no fim do gráfico). Isto pode se dever ao viés criado pela dependência dos controles de treinamento ou ao número insuficiente de simula¸cões. 10

10

6.5

−1.1


6.45 ) $ 6.4 M M (

−1.3

a d6.35 a m i t s 6.3 E a i d é6.25 M

E −1.5 S10 M R

10

−1.7

6.2 10


−1.9

10

1

2

10 Esforço Computacional (p)

(a) RMSE (V ar(˜ xn ))

3

10

6.15 1 10

2

10 Simulações HF

3

10

(b) Valor Esperado x ñ

Fig. 5.3: Compara¸caõ entre MC Multifidelidade e Tradicional - Reservatório Compress´ıvel e Densidades dos Fluidos Diferentes Os parâmetros de controle do método multifidelidade para o caso compress´ıvel estão apresentados na Figura 5.4. Observa-se de fato na Figura 5.4(a) um expressiva redu¸cão da correla¸caõ entre o modelo de alta e baixa fidelidade ρxy quando comparadas com o caso incompress´ıvel. Isto demonstra que a aproxima¸caõ TPWL/POD perde qualidade devido à inclusão das não linearidades no problema. Mas o que chama a aten¸caõ é que a correla¸caõ ρην entre as amostras para estima¸cão da variância é bem reduzida. Devido a esta redu¸cão na correla¸caõ, observa-se também um redu¸caõ da razão r entre o n´ umero de simula¸co˜es de alta e baixa fidelidade que flutuou em torno de 16 conforme Figura 5.4(b). Isto mostra que no caso incompress´ıvel com densidades diferentes a redu¸caõ do custo computacional para estima¸caõ dos parâmetros é menor do que no caso incompress´ıvel com densidades iguais devido a menor qualidade da aproxima¸cão de ordem reduzida.


126 25

0.8

Média Variância

Média Variância

0.7

20

0.6

15

0.5 r

ρ

0.4

10

0.3

5 0.2 0.1 1 10

2

10 Simulações HF

0 1 10

3

10

(a) Correla¸cão (ρXY )

2

10 Simulações HF

3

10

(b) Razão de no. de Sim’s (r )

0.4

Média Variância

0.35 0.3 0.25 β

0.2 0.15 0.1 0.05 1 10

2

10 Simulações HF

3

10

(c) Corre¸cão β

Fig. 5.4: Evolu¸caõ dos Parâmetros do MC Multifidelidade para Estima¸caõ da Média e da Variância - Reservatório Compress´ıvel e Densidades dos Fluidos Diferentes Finalmente, o que se observa na Figura 5.4(c) é que o parâmetro de corre¸caõ β é pequeno pois, com a qualidade menor da aproxima¸caõ, a influência do modelo de baixa fidelidade no estimador deve ser menor. Nota-se que para o estimador da variância, esta influência será ainda menor devido a seus baixos n´ıveis de correla¸cão.

˜ 6. CONCLUSOES E TRABALHOS FUTUROS Este trabalho de mestrado estudou a técnica Trajectory Piecewise Linearization combinada a Proper Orthogonal Decomposition como modelo de ordem reduzida aplicado ao gerenciamento de reservatórios de petróleo. O método foi aplicado aos processos de otimiza¸caõ de controles de po¸cos e de propaga¸caõ de incertezas nestes controles. Neste cap´ıtulo, é apresentado um breve resumo sobre os resultados e algumas conclus˜ oes sobre o uso e a robustez das aplica¸co˜es serão apresentadas. Além disso, algumas sugestões de trabalhos futuros serão mostradas com base nos problemas enfrentados, nas limita¸cões percebidas e nas ideias levantadas durante o desenvolvimento do trabalho. Todas as técnicas aqui combinadas tiveram origem em trablhos de outros contex´ importante portanto ressaltar que o presente trabalho pode ser aplicado a outras áreas tos. E do conhecimento, tanto na engenharia de reservatórios como em outros ramos da engenharia.

6.1

Resultados e conclus˜ oes Os principais resultados est˜ ao resumidos na lista a seguir.

• Estudou-se uma formula¸caõ black-oil ligeiramente diferente daquela implantada nos

simuladores comerciais. Esta é baseada em fra¸co˜es mássicas e está incorporada no simulador de reservatórios utilizado neste trabalho. Tal simulador é chamado SIMPAR, foi desenvolvido pela PETROBRAS na década de 90 e teve seu código aberto para identifica¸caõ dos pontos necessários para a exporta¸cão de dados de mapas de estados e derivadas.

• Foram constru´ıdos dois modelos de reservatórios para os testes das aplica¸cões: um ba127

6.1. Resultados e conclusões

128

seado em corte do SPE10 e outro derivado da conversão da malha do modelo de Brugge para cartesiana. Os resultados deste simulador foram comparados com os obtidos com o IMEX da CMG para valida¸caõ em dois modelos diferentes.

• Foi desenvolvido um modelo substituto de ordem reduzida baseado nas técnicas de Trajectory Piecewise Linearization (TPWL) combinada ao Proper Orthogonal Decom-

position (POD) em Matlab. Tal modelo foi utilizado para reduzir o tempo de simula¸cão de reservatórios utilizando o SIMPAR, conseguindo-se speedups da ordem de 350.

• A técnica TPWL/POD já havia sido aplicada ao problema de simula¸caõ o´leo-água de reservatórios em sua formula¸caõ por satura¸cões. O fato de o simulador utilizado consi-

derar fra¸cões mássicas como variáveis primárias não alterou a qualidade dos resultados.

• Também utilizou-se neste desenvolvimento a proje¸cão de Petrov-Galerkin em lugar de

Bubnov-Galerkin. Sendo assim, para solucionar-se o sistema de equa¸co˜es de ordem reduzida, pré-multiplicam-se as derivadas reduzidas por ΦT (J i+1 )T em lugar de Φ T , onde Φ é a base do POD e J i+1 é o jacobiano do estado gravado mais pr´ oximo. O uso desta proje¸caõ transforma a solu¸caõ por m´ınimos quadrados e foi capaz de estabilizar o TPWL/POD. De fato, depois da decisão de utilizá-la, nenhum problema de instabilidade foi experimentado. Até então, a literatura apontava como principal problema da técnica uma certa instabilidade, que parece ser resolvido pela utiliza¸cão da pro je¸caõ obl´ıqua em lugar da ortogonal.

• Foi proposto que o fluxo de trabalho de treinamento do modelo reduzido somente

contivesse uma u ´ nica simula¸caõ com controles escolhidos aleatoriamente e ciclos de controle mais curtos que o problema a ser resolvido. A troca aleatória e frequente dos controles proporciona uma boa explora¸caõ do espa¸co de estados e permitiu atingir-se um n´ıvel de erro aceitável para a aproxima¸cão.

• A escolha de uma u´ nica simula¸caõ para o treinamento se deve ao interesse de que não se espera que o modelo reduzido seja perfeito, já que poderá ser retreinado quando

necess´ ario, acrescentando-se assim uma nova trajet´ oria ao conjunto de estados. Isto permite uma razoável redu¸caõ no tempo de simula¸caõ para o treinamento, mas pode vir a acrescentar um tempo de retreinamento durante um processo de otimiza¸cão.

• Observou-se que, ao considerarem-se fluidos e rochas incompress´ıveis e nenhuma dife-

ren¸ca nas densidades do fluidos, o problema é pouco não-linear, proporcionando que os resultados obtidos com o modelo substituto TPWL/POD sejam praticamente perfei-

6.1. Resultados e conclusões

129

tos quando comparados aos obtidos com o simulador com erros da ordem de 3%, sem nenhum retreinamento.

• Ao acrescentarem-se primeiro a diferen¸ca de densidades dos fluidos e depois as compres-

sibilidades, fazendo do modelo mais realista, observa-se uma degrada¸caõ da qualidade dos resultados obtidos pelo modelo reduzido, atingindo-se um patamar de erro de 10%, mesmo após retreinamento. Isto se deve à forte não-linearidade do problema.

• O parâmetro de erro do TPWL/POD representa uma tentativa de medir a distância no espa¸co reduzido da trajetória simulada e a trajetória gravada utilizada. O critério de retreinamento baseado num valor máximo para esta distância evitou retreinamento excessivo, em situa¸co˜es onde o parâmetro de controle do MSAO se mostrava baixo por causa da proximidade do ponto de ótimo.

• Foi constru´ıdo um algoritmo multifidelidade de otimiza¸caõ sequencial aproximada

(MSAO) que foi aplicado ao problema de otimiza¸cão do VPL controlando-se as pressões de fundo dos po¸c os. Este algoritmo resolve o problema de otimiza¸caõ como uma sequência de problemas aproximados resolvidos no interior de uma região de confian¸ca que é atualizada a cada itera¸caõ. A aproxima¸cão utilizada é o modelo TPWL/POD corrigido por uma krigagem obtida com algumas amostras dos controles na região de confian¸ca.

• O algoritmo MSAO aqui desenvolvido obteve resultado com maior VPL do que o obtido

com o acoplamento do simulador diretamente ao algoritmo de programa¸caõ sequencial quadrática (SQP) em um dos problemas estudados. Isto se deve a alguma descontinuidade existente na fun¸cão-ob jetivo que é suavizada devido à aplica¸caõ da krigagem.

• A corre¸caõ aplicada ao modelo de ordem reduzida TPWL/POD que obteve melhor

resposta foi o kriging com polinômio interpolador de primeiro grau. Isto se deve ao fato de haver uma certa dependência linear entre o erro da aproxima¸caõ e a sequência de controles utilizada.

• Foi constru´ıdo também um algoritmo Monte Carlo Multifidelidade baseado na técnica

estat´ıstica das variáveis de controle. Este foi utilizado para a estima¸caõ da média e desvio padrão do VPL considerando-se incertas as pressõ es de fundo dos po¸cos. O modelo reduzido foi utilizado para reduzir-se o número de execu¸cões do simulador. O fator de redu¸caõ no n´ umero de simula¸co˜es é calculado baseado na correla¸caõ entre os resultados de alta e baixa fidelidade e a razão entre os custos computacionais.

6.2. Sugestões de Trabalhos Futuros

130

• O algoritmo Monte Carlo Multifidelidade foi aplicado ao modelo baseado no SPE10

e conseguiu obter melhores estimativas para a média e desvio padrão com um menor número de execu¸co˜es do SIMPAR quando comparado ao método Monte Carlo clássico. Entretanto, os resultados foram melhores com a versão sem compressibilidades de fluidos e rocha e com densidades iguais. Ao considerarem-se as compressibilidades e a diferen¸ca nas densidades, a correla¸caõ do TPWL/POD com o SIMPAR se reduz bastante o que faz com que o algoritmo necessite de mais execu¸co˜es do simulador.

6.2

Sugest˜ oes de Trabalhos Futuros Algumas sugestões de trabalhos futuros para dar continuidade a esta pesquisa e

acrescentar mais conhecimento a esta área estão listadas abaixo.

• Analisar em mais detalhes a estabilidade o método TPWL/POD, dando um enfoque mais matemático e resultando em um crit´ erio de estabilidade.

• Fazer a análise de erro do TPWL/POD de modo a obter um critério de retreinamento mais robusto.

• Implementar a possibilidade de controlar os po¸cos por vazão no modelo TPWL/POD. • Utilizar o TPWL/POD para representar de maneira aproximada as restri¸cões de pressão m´ınima de fundo de po¸co em um problema de otimiza¸caõ do VPL com controles de vazão.

• Aplicar o TPWL/POD ao cálculo de gradientes aproximados em algoritmos de otimiza¸caõ baseados em derivadas.

• Substituir a krigagem por uma rede neural no algoritmo MSAO. • Considerar, além dos controles de po¸cos, a varia¸cão de parâmetros geológicos na formula¸caõ do TPWL/POD (em (He, 2010; He et al., 2013) encontra-se um estudo para

transmissibilidades), implementando assim um modelo susbtituto a ser utilizado em estudos de análise de incertezas geológicas e ajuste de histórico.

• Utilizar a técnica das variáveis de controle e a ideia anterior para propaga¸cão de incertezas geológicas e otimiza¸caõ sob incertezas.

ˆ REFERENCIAS S. M. B. Afonso, B. Horowitz, and M. F. da Silva. Performance of a physical surrogate in the sequential approximate optimization of space trusses. International Journal for ComputerAided Engineering and Software , 2011. M. N. Albunni, V. Rischmuller, T. Fritzsche, and B. Lohmann. Model-order reduction of moving nonlinear electromagnetic devices. IEEE Transactions on Magnetics , 2008. N. Alexandrov, J. E. D. Jr, R. M. Lewisand, and V. Torczon. A trust region framework for managing the use of approximation models in optimization. Technical report, NASA, Hampton, VA, USA, 1997. A. H. Alhuthali, A. Datta-Gupta, B. Yuen, and J. P. Fontanilla. Field applications of waterflooding optimization via rate control with smart wells. In SPE Reservoir Simulation Symposium , 2009. SPE 118948-PA. L. F. Almeida, M. M. B. R. Vellasco, and M. A. C. Pacheco. Optimization system for valve control in intelligent wells under uncertainties. Journal of Petroleum Science and Engineering , 2010. M. Asadollahi and G. Naevdal. Selection of decision variables for large-scale production optimization problems applied to the brugge field. In SPE Russian Oil and Gas Conference and Exhibition , 2010. SPE 136378-MS. M. C. Biggs. Constrained minimization using recursive quadratic programming. Towards Global Optimization (L.C.W. Dixon and G.P. Szergo, eds.), 1979. D. R. Brower and J. D. Jansen. Dynamic optimization of waterflooding with smart wells using optimal control theory. SPE Journal , 2004. SPE 78278-PA. M. A. Cardoso. Development and application of reduced-order modeling procedures for reservoir simulation . PhD thesis, Stanford University, 2009.

ˆ REFERENCIAS

132

M. A. Cardoso and L. J. Durlofsky. Linearized reduced-order models for subsurface flow simulation. Journal of Computational Physics , Vol. 229(No. 3):pp. 681–700, 2010. K. Carlberg, C. Bou-Mosleh, and C. Farhat. Effcient nonlinear model reduction via a leastsquares petrov-galerkin projection and compressive tensor approximations. International Journal for Numerical Methods in Engineering , 2009. CENPES. Simpar 2.0 - Manual do Usu´ ario. Petrobras/CENPES, 1995. CENPES and ESSS. Cyclope - Manual do usuário. Petrobras/CENPES e ESSS, 2012. S. Chaturantabut and D. C. Sorensen. Discrete empirical interpolation for nonlinear model reduction. In 48th IEEE Conference , Dez/2009. C. Chen, Y. Wang, G. Li, and A. C. Reynolds. Closed-loop reservoir management on the brugge test case. Computational Geosciences , 2010. C. Chen, G. Li, and A. C. Reynolds. Robust constrained optimization of short and long-term net present value for closed-loop reservoir management. SPE Journal , 2012. SPE 141314-MS. Y. Chen and D. S. Oliver. Ensemble-based closed-loop optimization applied to brugge field. SPE Reservoir Evaluation Engineering , 2010. SPE 118926-MS. M. A. Christie and M. J. Blunt. Tenth spe comparative solution project: A comparison of upscaling techniques. SPE Reservoir Evaluation and Engineering , 2001. SPE 66599-MS. CMOST. CMOST - User’s Guide . Computer Modeling Group Ltd., Calgary, Alberta, Canada, 2012. A. R. Conn, K. Scheinberg, and L. N. Vicente. Introduction to derivative-free optimization. In SIAM , Philadelphia,USA, 2009. M. F. da Silva. Estrat´ egias de aproxima¸caõ para a otimiza¸caõ estrutural. Master’s thesis, Universidade Federal de Pernambuco, 2009. T. Ertekin, J. H. Abou-Kassem, and G. R. King. Basic Applied Reservoir Simulation . Society of Petroleum Engineering - SPE, 2011. A. I. J. Forrester, A. Sóbester, and A. J. Keane. Engineering Design via Surrogate Modelling: a Practical Guide . Wiley, 2008.

ˆ REFERENCIAS

133

E. Gildin, M. Ghasemi, A. Protasov, and Y. Efendiev. Nonlinear complexity reduction for fast simulation of flow in heterogeneous porous media. In SPE Reservoir Simulation Symposium , 2013. SPE 163618-MS. A. A. Giunta. Use of data sampling, surrogate models, and numerical optimization in engineering design. In 40th AIAA Aerospace Sciences Meeting and Exhibit , 2002. AIAA 2002-0538. A. A. Giunta and M. S. Eldred. Implementation of a trust region model management strategy in the dakota optimization toolkit. In 8th AIAA/USAF/NASA/ISSMO Symposium on Multidisciplinary Analysis and Optimization , 2000. AIAA 2000-4935. P. W. Glynn and R. Szechtman. Some New Perspectives on the Method of Control Variates, in Monte Carlo and Quasi-Monte Carlo Methods . Springer-Verlag, Berlin, Germany, 2000. HDFGroup. Hdf5 tutorial. http://www.hdfgroup.org/HDF5/Tutor/, 2014. Accessed: 201403. J. He. Enhanced linearized reduced-order models for subsurface flow simulation. Master’s thesis, Stanford University, 2010. J. He and L. J. Durlofsky. Reduced-order modelling for compositional simulation using trajectory piecewise linearization. In Reservoir Simulation Symposium , The Woodlands, Texas, USA, 2013. SPE 163634-MS. J. He, J. Saetrom, and L. J. Durlofsky. Enhanced linearized reduced-order models for subsurface flow simulation. Journal of Computational Physics , Vol. 230(No. 23):pp. 8313–8341, 2011. J. He, P. Sarma, and L. Durlofsky. Reduced-order modeling and geological parametrization for ensemble-based data assimilation. Computers and Geosciences , Vol. 55, 2013. B. Horowitz, S. M. B. Afonso, C. V. P. Mendon¸ca, and R. B. Willmersdorf. Rate control optimization of waterflooding management. International Journal of Modeling and Simulation for the Petroleum Industry , 2012. B. Horowitz, S. M. B. Afonso, and C. V. P. de Mendon¸ca. Surrogate based optimal water flooding management. Journal of Petroleum Science and Engineering , 2013. J. D. Jansen. Adjoint-based optimization of multiphase flow through porous media - a review. Computational Fluids , 2011.

ˆ REFERENCIAS

134

I. T. Jolliffe. Principal Component Analysis . Springer-Verlag, New York, 2 edition, 2002. D. R. Jones, M. Schonlau, and W. J. Welch. Efficient global optmization of expensive blackbox functions. Journal of Global Optimization , 1998. D. E. Knuth. The Art of Computer Programming , volume 2. Addison-Wesley, Massachussets, 3 edition, 1997. D. G. Krige. A statistical approach to some mine valuations and allied problems at the witwatersrand. Master’s thesis, University of Witwatersrand, 1951. J. Li and D. Xiu. Evaluation of failure probability via surrogate models. Journal of Computational Physics , 2010. S. N. Lophaven, H. B. Nielsen, and J. Søndergaard. DACE - A MATLAB Kriging Toolbox . IMM Technical University of Denmark, 2002. C. R. Maliska, A. F. C. da Silva, A. O. Czesnat, R. M. Lucianetti, and C. R. M. Junior. Three-dimensional multiphase flow simulation in petroleum reservoirs using the mass fractions as dependent variables. In Fifth Latin American and Caribean Petroleum Engineering Conference and Exhibition . Society of Petroleum Engineering - SPE, 1997. SPE 39067-MS. R. Markovinovic. System-theoretical model reduction for reservoir simulation and optimization . PhD thesis, Technische Universiteit Delft, 2003. G. Matheron. Principles of geostatistics. Economic Geology , 1963. Mathworks. Optimization Toolbox TM . User’s Guide . Mathworks, 2012.

†

N. Muller, L. Magaia, and B. M. Herbst . Singular value decomposition, eigenfaces, and 3d reconstructions. SIAM Review , 2004.

B. L. Nelson. Control variate remedies. Operations Research , vol 38(no. 6), Nov-Dez/1990. L. W. T. Ng and M. S. Eldred. Multifidelity uncertainty quantification using non-intrusive polynomial chaos and stochastic collocation. In 53rd AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics and Materials Conference , Honolulu, HI, 2012. L. W. T. Ng, D. B. P. Huynh, and K. E. Willcox. Multifidelity uncertainty propagation for optimization under uncertainty. In SIAM/ASA Journal on Uncertainty Quantification , 2013.

ˆ REFERENCIAS

135

D. F. Oliveira. Production optimization techniques for petroleum reservoirs: Derivative free approaches to dynamic rate allocation for injection and production. Master’s thesis, Universidade Federal de Pernambuco, 2006. D. W. Peaceman. Interpretation of well-block pressures in numerical reservoir simulation with nonsquare grid blocks and anisotropic permeability. SPE Journal , 1983. SPE 10528-PA. E. Peters, R. J. Arts, G. K. Brower, C. R. Geel, S. Cullick, R. J. Lorentzen, Y. Chen, K. N. B. Dunlop, F. C. Vossepoel, R. Xu, P. Sarma, A. H. Alhutali, and A. C. Reynolds. Results of the brugge benchmark study for flooding optmization and history matching. In SPE Reservoir Simulation Symposium , 2010. SPE 119094-MS. M. J. D. Powel. Algorithms for nonlinear constraints that use lagrangian functions. Mathematical Programming , 14, 1978. N. V. Queipo, J. V. Goicochea, and P. Pintos. Surrogate modeling-based optimization of sagd process. Journal of Petroleum Science and Engineering , 2002. M. J. Rewienski. A trajectory piecewise-linear approach to model order reduction of nonlinear dynamical systems . PhD thesis, Massachusets Institute of Technology, 2003. J. R. P. Rodrigues and L. Bonet. Desenvolvimento de um simulador black-oil. Technical report, Petrobras/CENPES, 1991. P. Sarma, W. H. Chen, L. J. Durlofsky, and K. Aziz. Production optimization with adjoint models under nonlinear control-state path inequality constraints. SPE Reservoir Evaluation Engineering , 2008. SPE 99959-PA. S. A. Souza, S. M. B. Afonso, and B. Horowitz. Optimal management of oil production using the particle swarm algorithm and adaptive surrogate models. In XXXI Iberian LatinAmerican Congress on Computational Methods in Engineering , Buenos Aires, 2010. P. K. W. Vinsome. Orthomin, an iterative method for solving sparse sets of simultaneous linear equations. In Fourth Symposium of Numerical Simulation of Reservoir Performance , Los Angeles, CA, USA, 1976. SPE 5729. C. Wang, G. Li, and A. C. Reynolds. Production optimization in closed-loop reservoir management. SPE Journal , 2009. SPE 109805-PA.

Aplicações de um Modelo Substituto de Ordem Reduzida

Recommend Documents