Seleção de Modelos - Comparativo de Modelos Matemáticos para Velocidade de Paraquedista

Exemplo Seleção de modelos O teste e seleção de modelos são atividades comuns e extremamente importantes em todos os campos de engenharia Considerase aqui um modelo matemático teórico para a velocidade de um paraquedista é dado pela equação abaixo vt gm 1 ecmt onde v é a velocidade ms g é a constante gravitacional m é a massa do paraquedista igual a 681 kg e c é o coeficiente de arrasto de 125 kgs O modelo prediz que a velocidade do paraquedista como função do tempo Um modelo empírico alternativo para a velocidade do paraquedista é dado por vt gmc t 375 t 616 Supondo que desejase comparar a adequação desses dois modelos matemáticos Isso pode ser feito medindo a velocidade real do paraquedista em valores conhecidos de tempo e comparando as velocidades preditas por cada modelo Tabela 61 Velocidades de queda de um paraquedista medidas e calculadas Tempo s Velocidade Medida ms Velocidade Calculada ms Modelo 1 Velocidade Calculada ms Modelo 2 1 1000 8953 11240 2 1630 16405 18570 3 23 22607 23729 4 2750 27769 27556 5 3100 32605 30509 6 3560 35641 32855 7 3900 38617 34766 8 4150 41095 36351 9 4290 43156 37687 10 4500 44872 38829 11 4600 48479 40678 12 4550 47409 40678 13 4600 48479 41437 14 4900 49303 42110 15 5000 49988 42712 Solução A adequação dos modelos pode ser testada plotando a velocidade calculada pelo modelo contra a velocidade medida empiricamente Uma regressão linear pode ser usada para calcular o coeficiente angular e o intercepto do gráfico Essa linha então deve ter coeficiente angular 1 e intercepto 0 e um r2 1 se o modelo se ajusta perfeitamente aos dados Um desvio significativo desses valores pode ser usado como indicador de inadequação do modelo As figuras 61a e 61b são gráficos dos dados e da linha de regressão das colunas b e c respectivamente versus a coluna a Para o primeiro modelo vmodelo 0859 1032vmedida 617 106 E para o segundo modelo vmodelo 5776 0752vmedida 618 Esses gráficos indicam que a regressão linear entre os dados experimentais e cada um dos modelos é altamente significante Ambos os modelos ajustam os dados com um coeficiente de correlação maior do que 099 Figura 61 Adequação de Modelos No entanto primeiro modelo descrito confirma a hipótese do critério de teste muito melhor do que o segundo modelo já que o o coeficiente angular e o intercepto estão mais próximos de 1 e 0 Dessa forma apesar de cada gráfico ser bem representado por uma linha reta o primeiro modelo parece ser um melhor modelo do que o segundo Partidose desse resultado é possível fazer um estudo mais aprofundado do modelo que melhor ajustou os dados Nesse caso é necessário recalcular os dados utilizando a abordagem matricial para estimar os erros para os parâmetros e então empregar esses erros no desenvolvimento de intervalos de confiança e usalos para tirar conclusões estatísticas sobre a qualidade do ajuste Esses dados podem ser então escritos na forma de matriz para uma regressão linear Z 1 10 1 163 1 23 1 2750 1 31 1 356 1 39 1 415 1 429 1 45 1 46 1 455 1 46 1 49 1 50 Y 8953 16405 22607 27769 32065 35641 38617 41095 43156 44872 46301 47490 48476 49303 49988 A transposição e multiplicação das matrizes pode ser então utilizada para gerar as equações normais tal como 107 A ZT Z1 ZT Y 0689414 001701 001701 0000465 552741 2242143 085872 1031592 Portanto o intercepto e o coeficiente angular são determinados como a0 085872 e a1 1031592 respectivamente Esses valores por sua vez podem ser usados para calcular o erro padrão da estimativa como syx 08634403 Esse valor pode ser usado juntamente com os elementos da diagonal da matriz inversa para calcular o erro padrão dos coeficientes sa0 sqrtz111 syx2 sqrt0688414 08634032 0716372 sa0 sqrtz111 syx2 sqrt0688414 08634032 0716372 A estatística tα2n1 necessária para um intervalo de confiança de 95 com n2 13 graus de liberdade pode ser determinado através de uma tabela estatística ou utilizando o Excel função INVT ou INVT cuidado pois essas possuem inputs diferentes Isso resulta no valor de t2513 2160368 então as equações 614 e 615 podem ser utilizadas para calcular os intervalos de confiança a0 085872 2160368 0716372 085872 1547627 240634 0688912 a1 1031592 2160368 0018625 1031592 0040237 0991355 0071828 Notase que os valores desejados 0 para o intercepto e 1 para o coeficiente angular ficam nesses intervalos Com base nessa análise é possível fazer uma afirmação quanto ao coeficiente angulares Temse grandes evidências para acreditar que o coeficiente angular da regressão linear para a linha real tem seu valor no intervalo de 0991355 a 1071828 Por esse motivo temse também fortes evidências para acreditar que o resultado sustenta a relação entre as medições e o modelo Devido ao fato de que o intervalo no qual se encontra o intercepto contém o valor de zero uma conclusão similar pode ser tirada Como mencionado anteriormente as equações normais são geralmente mal condicionadas Portanto se resolvidas através de técnicas convencionais tal como decomposição LU os coeficientes computados podem se altamente suscetíveis a erros de arredondamento Como consequência algoritmos de ortogonalização tal como fatorização QR podem ser utilizados para contornar o problema no entanto isso vai além do escopo desse texto 63 Regressão nãolinear Há muitos casos em engenharia onde modelos não lineares devem ser ajustados a um certo conjunto de dados No presente contexto esses modelos são definidos como aqueles que possuem uma dependência não linear com seus parâmetros Por exemplo fx a0 1 ea1x e 619 108 Essa equação não pode ser manipulada para se adequar à forma geral da equação 61 Tal como mínimos quadrados no caso linear a regressão nãolinear é baseada em determinar os valores dos parâmetros que minimizam a soma dos quadrados dos resíduos No entanto para o caso nãolinear a solução deve proceder de maneira iterativa O método de GaussNewton é um algoritmo para minimizar a soma dos quadrados dos resíduos entre dados e as equações não lineares O conceito chave que fundamenta a técnica é que uma expansão em série de Taylor é usado para expressar a equação não linear original de maneira aproximada ou seja há a linearização da função Então a teoria dos mínimos quadrados pode ser utilizada para obter novas estimativas de parâmetros que se movem em direção de minimizar o resíduo Para ilustrar como isso é feito primeiro a relação entre a equação nãolinear e os dados podem ser expressos genericamente por yi fxi a0 a1 am ei 620 onde yi é um valor medido da variável dependente fxi a0 a1 am é a equação que é uma função da variável independente xi e uma função nãolinear dos parâmetros a0 a1 am e ei é o erro aleatório Por conveniência esse modelo pode ser expresso de maneira abreviada omitindose os parâmetros yi fxi ei 621 O modelo nãolinear pode ser expandido em uma série de Taylor em torno dos valores dos parâmetros e truncados depois da primeira derivada Por exemplo para um caso de dois parâmetros fxij1 fxij fracpartial fxijpartial a0 Delta a0 fracpartial fxijpartial a1 Delta a1 622 onde j é o chute inicial j 1 é a predição Delta a0 a0j1 a0j e Delta a1 a1j1 a1j Dessa forma o modelo geral está linearizado com respeito a seus parâmetros e a equação A equação 622 pode ser então substituída por 621 para resultar em yi fxij fracpartial fxijpartial a0 Delta a0 fracpartial fxijpartial a1 Delta a1 ei 623 ou na forma matricial D ZjDelta A E 624 onde Zj é a matriz das derivadas parciais da função avaliada no chute inicial j Zj beginbmatrix fracpartial f1partial a0 fracpartial f1partial a1 fracpartial f2partial a0 fracpartial f2partial a1 vdots vdots fracpartial fnpartial a0 fracpartial fnpartial a1 endbmatrix 625 onde n é o número de dados e partial fi partial ak é a derivada parcial da função com respeito ao késimo parâmetro avaliado no iésimo ponto experimental O vetor D contém as diferenças entre as medições e os valores das funções D beginbmatrix y1 fx1 y2 fx2 vdots yn fxn endbmatrix 626 e o vetor Delta A contém as mudanças nos valores dos parâmetros A beginbmatrix Delta a0 Delta a1 vdots Delta am endbmatrix 627 Aplicando a teoria dos mínimos quadrados linear à equação 624 resulta nas seguintes equações normais ZjTZj Delta A ZjTD 628 Assim a abordagem consiste em resolver a equação 628 para A a qual pode ser empregada para calcular valores melhorados para os parâmetros tal como em a0j1 a0j Delta a0 629 e a1j1 a1j Delta a1 630 Esse procedimento é repetido até que a solução convirja ou seja até que varepsilonak fracakj1 akjakj1 100 631 atinja o valor de parada estabelecido Exemplo Método de GaussNewton Desejase ajustar a função fx a0 a1 a01 ea1 x aos dados x 025 075 125 175 225 y 028 057 068 074 079 Utilizando chutes iniciais de a0 10 e a1 10 para os parâmetros Note que para esses chutes a soma inicial dos resíduos é 00248 Solução As derivadas parciais da função com respeito aos parâmetros são fracpartial fpartial a0 1 ea1 x fracpartial fpartial a0 a0 x ea1 x As equações acima podem ser utilizadas para avalia a matriz Z0 beginbmatrix 02212 01947 05276 03543 07135 03581 08262 03041 08946 02371 endbmatrix Essa matriz multiplicada por sua transposta resulta em Z0T Z01 beginbmatrix 36397 78421 78421 191678 endbmatrix 632 O vetor D consiste de diferenças entre as medidas e as predições do modelo D beginbmatrix 028 02212 057 05276 068 07135 074 08262 079 08946 endbmatrix beginbmatrix 00588 00424 00335 00862 01046 endbmatrix 633 É multiplicado por Z0T para resultar em Z0T D beginbmatrix 01533 00365 endbmatrix 634 O vetor Delta A é então calculado resolvendo a Equação 628 para Delta A beginbmatrix 02714 05019 endbmatrix 635 que pode ser adicionado os chutes iniciais para resultar em beginbmatrix a0 a1 endbmatrix beginbmatrix 10 10 endbmatrix beginbmatrix 02714 05019 endbmatrix beginbmatrix 07286 15019 endbmatrix 636 Então as estimativas dos parâmetros são a0 07286 e a1 15019 os novos parâmetros resultam em uma soma dos quadrados dos resíduos igual a 00242 A Equação 631 pode ser utilizada para computar varepsilon0 e varepsilon1 igual a 37 e 33 porcento respectivamente O cálculo pode ser várias vezes repetido até que esses valores fiquem abaixo do critério de parada estabelecido O resultado final é a0 079186 e a1 16751 Esses coeficientes resultam em uma soma dos quadrados dos resíduos de 0000662 Um problema em potencial com o método de GaussNewton tal como desenvolvido até aqui é que as derivadas parciais da função em questão podem ser difíceis de de avaliar ΦA 0 646 Como Y é nãolinear com respeito aos parâmetros a equação 646 resultará em uma equação nãolinear que pode ser de difícil resolução para A Esse problema foi amenizado por Gauss que determinou que o ajuste de funções por mínimos quadrados pode ser alcançado por um método iterativo envolvendo uma série de aproximações lineares Em cada estágio de iteração a teoria dos quadrados lineares pode ser usada para obter a próxima aproximação Esse método conhecido como o método de GaussNewton converte o problema nãolinear em um problema linear aproximando a função Y por uma expansão de Taylor em torno de um valor estimado do vetor de parâmetros A FxA FxAk ΔA FxAk FAbk ΔA Fx Zj ΔA 647 Onde a série de Taylor é truncada depois do segundo termo A Equação 647 é linear em ΔA Dessa maneira o problema se torna de encontrar A para encontrar a correção de A isto é ΔA que deve ser adicionado à estimativa de A para minimizar a soma dos quadrados dos resíduos Para fazer isso Y é substituído na Equação 641 com o lado direito da Equação 647 para fornecer Φ Y Fx Zj ΔAT Y Fx Zj ΔA 648 Tomando a derivada parcial de Φ com respeito a ΔA igualando a zero e resolvendo para ΔA obtémse ΔA ZjT Zj1 ZjT Y Fx 649 O método de GaussNewton se aplica a ambos problemas modelos uni e multivariável O algoritmo do método GaussNewton envolve os seguintes passos 1 Assuma um chute inicial para o vetor de parâmetros A 2 Se o modelo está na forma de equações diferenciais use o vetor A e as condições de contorno para integrar as equações para obter os perfis de Y Se o modelo está na forma de equações algébricas então use o vetor A para avaliar Y a partir das equações 3 Calcule a matriz Jacobiana Zj das equações do modelo 4 Utilize a Equação 649 para ober a correção do vetor ΔA 5 Calcule a nova estimativa do vetor dos parâmetros a partir da equação 643 Ak1 Ak ΔA 6 Também é possível aplicar o fator de relaxação para evitar que os cálculos divirjam 7 Repita os passos 25 até que uma ou ambas das condições a seguir sejam satisfeitas Φ não muda significativamente ΔA se torne muito pequeno O método de GaussNewton é baseado na linearização de um modelo nãolinear dessa forma é esperado que esse método funcione bem nos casos em que o modelo não é altamente nãolinear ou se a estimativa inicial do vetor de parâmetros esteja próxima da soma mínima dos quadrados dos resíduos Os contornos de Φ constante no espaço dos parâmetros de um modelo linear são elipsoides Para um modelo nãolinear esses contornos são distorcidos mas na vizinhança do mínimo Φ os contornos estão muito próximos de serem elípticos Portanto podese dizer que o método de GaussNewton é bem efetivo se o ponto de início para a busca do mínimo está próxima à região elíptica Por outro lado esse método pode divergir se o ponto de início das iterações está numa região muito distorcida do hiperespaço dos parâmetros 643 Método de Newton A Equação 646 representa um conjunto de equações nãolineares assim o método de Newton pode ser aplicado para resolver esse conjunto de equações nãolineares Primeiramente é possível expandir Φ por série de Taylor até o terceiro termo ΦxA ΦxAk ΦAk ΔA 12 ΔAT 2ΦA2k ΔA 650 Tomandose a derivada parcial de ambos os lados da equação 650 com respeito a A resulta em ΦA ΦAk 2ΦA2k ΔA 651 A primeira derivada de Φ com respeito a A pode ser calculada derivandose a Equação 641 ΦAk 2ZjT Y Fx 652 E a segunda derivada de Φ com respeito a A é chamada matriz Hessiana das derivadas de segunda ordem de Φ com respeito a A avaliadas em todos os n pontos onde as observações experimentais estão disponíveis H 2Φa12 2Φa1a2 2Φa1am 2Φa2a1 2Φb22 2Φb1bm 2Φama1 2Φbmb2 2Φb22m 653 Aplicandose a condição necessária de se ter um mínimo local de Φ Equação 646 na Equação 651 e combinando com as equações 652 e 653 podese calcular o vetor de correção Δb ΔA 2H1ZjT Y Fx 654 É interessante notar que no caso da regressão de um único parâmetro a Equação 646 se torna dΦda Φ 0 655 E a equação 654 é simplificada para Δa Φk Φk 656 que é a solução da equação não linear Φ 0 através do método de NewtonRaphson O procedimento de cálculo para o método de Newton é quase o mesmo que o método de GaussNewton com exceção de que o vetor de correções dos parâmetros é calculado da equação 654 Se Φ é quadrático com respeito a A isto é regressão linear o método de Newton converge em somente um passo Tal como todos os outros métodos aplicandose a técnica de Newton para a solução do conjunto de equações nãolineares um fator de relaxação pode ser usado juntamente com a Equação 654 na correção dos parâmetros 644 Método de Marquardt Marquardt desenvolveu uma técnica de interpolação entre os métodos GaussNewton e descida mais ingrime Essa interpolação é atingida pela adição da matriz diagonal λI à matriz ZjT Zj na equação 649 ΔA ZjT Zj λI1 ZjT Y Fx 657 O valor de λ é escolhido em cada iteração de forma que o parâmetro corrigido resulte na menor soma dos quadrados na iteração seguinte Pode ser facilmente notado que quando o valor de λ é pequeno em comparação com os elementos da matriz ZjT Zj o método de Marquardt se aproxima do método de GaussNewton quando λ é muito grande esse método se torna idêntico ao método da descida mais ingrime com a diferença de que o fator de escala não afeta a direção do parâmetro de correção mas gera um passo pequeno De acordo com Marquardt é desejável minimizar Φ na vizinhança do máximo sobre o qual a função linearizada fornece uma representação adequada da função nãolinear Portanto o método para escolher λ deve fornecer valores pequenos de λ em situações onde o método de GaussNewton convergiria mais eficientemente e grandes valores de λ em situações onde o método de descida mais ingrime é mais adequado O método de Marquardt pode ser aplicado da mesma forma ao método de Newton Nesse caso a matriz diagonal λI é adicionada à matriz Hessiana na equação 654 Δb 2H λI1 Zj Y Fx 658 O método de Marquardt consiste nos seguintes passos 1 Assuma um chute inicial para o vetor de parâmetros A 2 Assuma um valor por exemplo 1000 para λ Isso significa que na primeira iteração o método de descida ingrime é predominante e deveria assegurar que o método está movendo em direção da menor soma dos quadrados dos resíduos 3 Calcule a matriz Zj a partir das equações do modelo Também avalie a matriz Hessiana H utilizando o método de Newton 4 Utilize a Equação 657 ou a Equação 658 para obter o vetor de correção ΔA 5 Calcule a nova estimativa do vetor de parâmetros a partir da Equação 643 Ak1 Ak ΔA 659 6 calcule o novo valor de Φ Se Φk1 Φk reduza λ por um fator de 4 por exemplo Se Φk1 Φm mantenha os parâmetros anteriores e aumente o valor de λ por um fator de 2 por exemplo 7 Repita os passos 36 até que um ou ambas das seguintes condições seja satisfeita Φ não muda significativamente ΔA se torna muito pequeno 65 Regressão nãolinear múltipla Anteriormente discutindo sobre regressão nãolinear a soma dos quadrados dos resíduos que foi minimizada foi aquela dada pela equação Φ DTD Y FxTY Fx 660 Essa era a soma dos quadrados dos resíduos determinada a partir do ajuste de uma equação às medições de uma variável Contudo a maioria dos modelos matemáticos pode envolver equações simultâneas com múltiplas variáveis dependentes Para tais casos quando mais do que uma equação é ajustada a dados multiresposta onde há ν variáveis dependentes no modelo a soma ponderada dos quadrados dos resíduos é dada por Φ k1ν wkDTD k1ν wkφk k1ν wkYFxTYFx 661 onde wk fator de ponderação correspondente à késima variável dependente φk soma do quadrado dos resíduos correspondente à késima variável Para minimizar Φ pelo método de GaussNewton primeiramente se lineariza os modelos utilizando a equação YxA YxAm ΔA YxAm YA ΔA Y ZjΔA 662 E combinandoa com a Equação 661 para obter Φ k1ν wkYk Fxk Zjk ΔAT Yk Fxk Zjk ΔA 663 Tomando a derivada parcial de Φ com respeito a ΔA igualando a zero e resolvendo para ΔA obtémse ΔA k1ν wkZjkT Zjk1 k1ν wkZjkT Yk Fxk 664 A Equação 665 dá a correção do vetor de parâmetros quando ajustando múltiplas variáveis dependentes simultaneamente A Equação 665 tornase idêntica à 649 quando ν 1 ou seja quando somente uma variável dependente é ajustada Utilizando o método de Marquardt a correção do vetor de parâmetros é calculada a partir de ΔA λI k1ν wkZjkT Zjk1 k1ν wkZjkT Yk Fxk 665 Os fatores de ponderação wk são determinados como se segue A premissa básica na derivação do algoritmo de regressão foi de que a variância da distribuição do erro nas medições era constante por todo o perfil de uma única variável dependente No entanto no caso de uma regressão múltipla é muito incomum que a variância σ2 de todas as curvas seja a mesma Portanto com a finalidade de se construir uma soma dos quadrados dos resíduos imparcial as somas dos quadrados dos resíduos individuais devem ser multiplicadas por um fator de ponderação que seja proporcional a aσ2 A equação para calcular os fatores de ponderação é dada por wk 1σ2 1i1ν ni l1ni 1σ2 666 onde σj2 ou σi2 variância para cada curva ni número de pontos experimentais disponíveis para cada curva ν número de variáveis sendo ajustadas O denominador da Equação 666 leva em consideração a possibilidade de que cada curva possa ter um número de pontos experimentais diferentes ni e os pondera de acordo Se a hipótese de que σ2 é constante dentro de uma das curvas não se aplica então a Equação 666 pode ser expandida para que o fator de ponderação possa ser calculado em cada ponto com o valor apropriado de σ2 Na maioria dos casos os valores de σ2 podem não ser conhecidos no entanto as estimativas dessas variâncias sk2 podem ser obtidas de experimentos repetidos e os valores de sk2 são então utilizados na Equação 666 para calcular os fatores de ponderação No pior dos casos onde não há a repetição de experimentos e não há conhecimento à priori de σ2 então os valores de wk devem ser supostos Caso contrário o algoritmo de regressão nãolinear pode introduzir uma tendência na direção de ajustar mais satisfatoriamente a curva com φk mais elevado e ignorar as curvas com baixo φk A regressão nãolinear pode ser também estendida para ajustar valores experimentais múltiplos da variável dependente em cada valor da variável independente Isso pode ser feito alterandose a Equação 661 para que os quadrados dos resíduos sejam somados também dentro de cada grupo de pontos Finalmente se o valor da variância do erro é proporcional ao valor da variável dependente o resíduo na soma dos quadrados deve ser dividido pelo valor teórico calculado da variável dependente em cada ponto do cálculo Bibliografia 1 A Constantinides and N Mostoufi Numerical Methods for Chemical Engineers with MATLAB Applications Prentice Hall PTR 1999 2 S C Chapra and R P Canale Numerical Methods for Engineers McGrawHill 2015 Modelagem e Simulação de Processos II Prof Dr Júlio Cesar de Carvalho Miranda