Introdução à Metodologia de Superfície de Resposta Uso: modelagem e análise de problemas nos quais a variável de resposta de interesse é influenciada por diversas variáveis independentes ou fatores e cujo objetivo é otimizar a variável resposta. Exemplo: o pesquisador está interessado em encontrar os níveis de temperatura, tempo e pH que maximizam a produção de um processo. A produção é um função da temperatura, tempo e do pH, ou seja, y f ( x1, x2 , x3 ) Onde, representa o ruído branco ou os erros aleatórios observados na reposta y. O valor esperado da resposta, E(y), é dado por: E( y) f ( x1, x2 , x3 ) Então, a superfície representada por: f ( x1, x2 , x3 ) é chamada de superfície de resposta. 1 Representação de uma superfície de resposta * Geralmente é feita graficamente: SAUPERFÍCIE DE RESPOSTA E GRÁFICO DE CONTORNOS Exemplo: Schneider e Stockett (1963) citados por JOHN (1971) realizaram um experimento com o objetivo de verificar a influência dos fatores Temperatura (x1), taxa gás/líquido (x2) e altura da embalagem (x3), na redução do odor desagradável de um produto químico que está sendo estocado para uso residencial. (Arquivo SAS: reducaoodorsuperficieresposta.sas) 1 40 120 40 120 40 120 40 120 80 80 80 80 80 80 80 Variáveis naturais 2 0,3 0,3 0,7 0,7 0,5 0,5 0,5 0,5 0,3 0,7 0,3 0,7 0,5 0,5 0,5 3 4 4 4 4 2 2 6 6 2 2 6 6 4 4 4 x1 -1 1 -1 1 -1 1 -1 1 0 0 0 0 0 0 0 Variáveis codificadas x2 -1 -1 1 1 0 0 0 0 -1 1 -1 1 0 0 0 Resposta x3 0 0 0 0 -1 -1 1 1 -1 -1 1 1 0 0 0 66 39 43 49 58 17 -5 -40 65 7 43 -22 -31 -35 -26 2 A codificação utilizada foi: x1 14080 x3 324 x2 20,02,5 Projeto experimental: DELINEAMENTO BOX-BEHNKEN (delineamnetos com 3 níveis). Geometricamente o delineamento fica: ° ° ° +1 ° ° ° x3 ° ° ° +1 ° ° -1 ° -1 ° -1 x1 x2 +1 3 Superfície de resposta Gráfico de contornos. Cada contorno corresponde a uma particular altura da superfície de resposta. Em cada linha a resposta é constante. 4 Na maioria dos problemas em superfície de resposta, a forma do relacionamento entre as variáveis dependentes e independentes é desconhecida. Assim, o primeiro passo é encontrar uma aproximação para o verdadeiro relacionamento entre a variável resposta (y) e as variáveis independentes (fatores).Geralmente utiliza-se de uma regressão polinomial de baixo grau em alguma região das variáveis independentes. Exemplo: o modelo de regressão polinomial de primeiro grau é dado por: y 0 1x1 2 x2 ... k xk O modelo de segunda ordem é dado por: k k k 1 k i 1 i 1 i 1 j i y 0 i xi ii xi2 ij xi x j Efeito linear Efeito quadrático Efeito da interação Estes modelos geralmente funcionam bem para uma região relativamente pequena do espaça dos fatores. 5 Os parâmetros do modelo são mais adequadamente estimados se forem utilizados planos adequados para a coleta dos dados. Os planos para ajsutar superfícies de resposta são denominados de delineamentos para superfície de resposta. Estes serão discutidos mais adiante. A metodologia de superfície de resposta é um procedimento sequêncial. Quando estamos num ponto da superfície de resposta que está longe do ótimo, como na condição operacional atual da figura, há pouca curvatura no sistema e o modelo de 1ª ordem será apropriado. • Região do processo •• • • • • •• • • • • • • Região de ótimo Caminho para a região de ótimo Condição operacional atual (está longe do ótimo 6 O objetivo é auxiliar o pesquisador, de forma rápida e eficiente, a encontrar a região de ótimo, isto é, determinar a melhor região de estudo. Encontrada a região de ótimo, um modelo mais elaborado, por exemplo, um modelo de segunda ordem, pode ser empregado, e uma análise pode ser feita para localizar o ponto de máximo ou de mínimo (ponto ótimo). Um outro objetivo da MSR é determinar as condições de operação ótima para o sistema, ou determinar uma região do espaço dos fatores no qual as especificações (requerimentos) de operação são satisfeitas. Método da inclinação ascendente (steepest ascent) Frequentemente, as condições iniciais (os pontos iniciais, região inicial) estão afastadas daqueles que otimizam a resposta. Em tais condições, o objetivo é mover o experimento rapidamente para a vizinhança geral do ótimo utilizando um procedimento experimental, simples, rápido, econômico e eficiente. Quando se está distante do ótimo, vamos assumir um modelo de primeira ordem como aproximação da verdadeira superfície de resposta em uma pequena região das variáveis independentes (xi) . Se estamos buscando o máximo incremento na resposta temos o método da máxima inclinação ascendente (steepest ascent); se estamos buscando um ponto de mínimo o método chama-se máxima inclinação descendente (steepest descent). 7 O modelo ajustado de primeira ordem é: k yˆ ˆ0 ˆi xi i 1 O gráfico de contornos dos valores preditos da variável resposta ( y chapéu), é uma série de linhas paralelas, como na figura, Caminho da inclinação ascendente (É a direção em que os valores ajustados aumentam mais rapidamente) x2 Região dos valores preditos pelo modelo de primeira ordem yˆ 10 yˆ 20 yˆ 30 yˆ 40 yˆ 50 x1 Figura. Superfície de resposta de primeira ordem e o caminho da inclinação ascendente 8 Os passos ao longo do caminho são proporcionais aos sinais e grandezas dos coeficientes de regressão { ˆ i }. O tamanho real do passo é determinado pelo pesquisador, baseado em considerações práticas ou conhecimento do processo. Experimentos (tratamentos) são conduzidos ao longo do caminho da inclinação ascendente até que não ocorre mais acréscimos na resposta. Um novo modelo de primeira ordem pode ser ajustado, um novo caminho de inclinação ascendente determinado, e o processo continuado. Eventualmente, o pesquisador pode chegar na vizinhança do ponto ótimo. Isto é indicado pela falta de ajuste do modelo de primeira ordem. Neste momento, experimentos adicionais (tratamentos) são realizados para obter uma estimativa mais precisa do ótimo. 9 Exemplo 1: Um engenheiro químico está interessado em determinar os níveis de tempo e temperatura de reação que maximizam a produção de um processo. Normalmente, opera-se com um tempo de 35 minutos e temperatura de 155 oF, que resulta numa produção de 40% aproximadamente. Como esta região provavelmente não contém o ótimo, ajusta-se um modelo de primeira ordem e aplica-se o Método da Máxima Inclinação Ascendente. O engenheiro decide que a região experimental será: Tempo: 1=30 2=40 Temperatura: 1=150 2=160 10 Codificação (para simplicação dos cálculos: x1 1 35 5 x2 2 155 5 Delineamento experimental: o delineamento experimental utilizado é um fatorial 22 aumentado de 5 pontos centrais. As repetições no ponto central são utilizadas para estimar o erro experimental e para checar o ajuste do modelo de primeira ordem. Os pontos centrais do delineamento são os correspondentes às condições de operação atual. Arquivo SAS: producaoexemplosteepestascent.sas) Fatorial Pontos centrais Variáveis originais 1 2 30 150 30 160 40 150 40 160 35 155 35 155 35 155 35 155 35 155 Variáveis naturais x1 x2 -1 -1 -1 1 1 -1 1 1 0 0 0 0 0 0 0 0 0 0 Resposta y 39.3 40.0 40.9 41.5 40.3 40.5 40.7 40.2 40.6 Modelo ajustado de 1a ordem: yˆ 40,44 0,775x1 0,325x2 ( R2 94,1%) Falta de ajuste: F 47 ,82 valor p 0,0002 modeloajustado 11 Diagnóstico do modelo de 1a ordem: • Obter uma estimativa do erro •Checar se interações devem ser incluídas no modelo •Checar se termos quadráticos devem ser incluídos no modelo 12 Estimativa do erro experimental [com os pontos centrais (repetição)] ˆ 2 ( 40, 3)2 ( 40,5)2 ( 40, 7 )2 ( 40, 2 )2 ( 40, 6 )2 [( 202 , 3)2 / 5] 51 0,0430 Estimativa de mínimos quadrados do coeficiente da interação ˆ12 14 1(39,3) 1(41,5) 1(40,0) 1(40,9) 0,025 Soma de quadrados da interação, com 1 grau de liberdade: SQint eração ( 0 ,1) 2 4 0,0025 Teste F para interação: 0025 F 00,,0430 0,058NS Estudo do efeito quadrático: a análise de variância completa obtemos através do programa estatístico SAS. 13 Response Surface for Variable YIELD Regression Degrees of Freedom Linear Quadratic Crossproduct Total Regress Residual Type I Sum of Squares 2 1 1 4 2.825000 0.002722 0.002500 2.830222 Degrees of Freedom Lack of Fit Pure Error Total Error Sum of Squares 0 4 4 0 0.172000 0.172000 R-Square F-Ratio Prob > F 0.9410 0.0009 0.0008 0.9427 32.849 0.0633 0.0581 16.455 0.0033 0.8137 0.8213 0.0095 Mean Square F-Ratio . 0.043000 0.043000 . ** N.S. N.S. ** Prob > F . Os efeitos da interação e da curvatura não foram significativos, ao passo que o efeito linear é significativo, portanto, o modelo de 1a ordem está ajustado. O erro padrão dos coeficientes lineares vale: ep ( ˆi ) QM erro 4 ˆ 2 4 0, 043 4 0,10 i 1,2 14 Ambos os coeficientes de regressão são maiores do que os seus erros padrões (2 x erro padrão). Podemos considerar que o modelo de 1a ordem está ajustado aos dados. Direção da máxima inclinação ascendente: Para andar (mover-se) do centro do delineamento (x1=0 e x2=0) no caminho da inclinação ascendente, deveríamos mover 0,775 unidades na direção x1 para cada 0,325 unidades na direção de x2. Assim, a direção da inclinação ascendente passa pelo ponto central (x1=0 e x2=0) e tem inclinação 0,325/0,775=0,42. O engenheiro decide usar um tempo de reação de 5 minutos como tamanho do passo inicial. Usando a relação entre 1 e x1, vimos que 5 minutos no tempo de reação corresponde a um intervalo (passo), na variável codificada x1, de x1=1. Os passos no caminho da inclinação ascendente são: x1=1 x2=(0,325/0,775) x1=0,42. 15 Os pontos experimentais são obtidos e a produção para estes pontos observados até que se perceba um decréscimo na produção. Os resultaods são mostrados na tabela a seguir: Para converter o tamanho dos passos codificados (x1=1 e x2=(0,325/0,775) x1=0,42) para as unidades originais de tempo e temperatura, usamos as relações: x1 1 5 x2 2 5 1 x1 (5) 1.0(5) 5 2 x2 (5) 0,42(5) 2 16 Passos Origem Origem + Origem + 2 Origem +3 Origem +4 Origem +5 Origem +6 Origem + 7 Origem + 8 Origem + 9 Origem + 10 Origem + 11 Origem + 12 Variáveis codificadas x1 x2 0 0 1,00 0,42 1,00 0,42 2,00 0,84 3,00 1,26 4,00 1,68 5,00 2,10 6,00 2,52 7,00 2,94 8,00 3,36 9,00 3,78 10,00 4,20 11,00 4,62 12,00 5,04 Variáveis originais 1 2 35 155 5 2 40 157 45 159 50 161 55 163 60 165 65 167 70 169 75 171 80 173 85 175 90 179 95 181 Resposta y Não faz 41,0 42,9 47,1 49,7 53,8 59,9 65,0 70,4 77,6 80,3 76,2 75,1 17 Observa-se um acréscimo na resposta até o passo 10; entretanto, depois deste passo têm-se um decréscimo na resposta. Portanto, outro modelo de primeira ordem deve ser ajustado na vizinhança do ponto (1 =85 e 2 =175). A região de 1 é [80, 90] e para 2 é [170, 180]. O delineamento experimental e os resultados são apresentados na tabela a seguir. Novamente usou-se um fatorial 22 completo com cinco pontos centrais. 18 Variáveis codificadas x1 x2 -1 -1 -1 1 1 -1 1 1 0 0 0 0 0 0 0 0 0 0 x1 1 85 5 Variáveis originais 1 2 80 170 80 180 90 170 90 180 85 175 85 175 85 175 85 175 85 175 e x2 Resposta y 76.5 77.0 78.0 79.5 79.9 80.3 80.0 79.7 79.8 2 175 5 19 O modelo de primeira ordem ajustado aos dados é: yˆ 78,97 1,00x1 0,50x2 A análise de variância para este modelo, incluindo os termos da interação e a regressão quadrática fica: (saída do SAS) Regression Degrees of Freedom Linear Quadratic Crossproduct Total Regress Residual Lack of Fit Pure Error Total Error 2 1 1 4 Degrees of Freedom 0 4 4 Type I Sum of Squares 5.000000 10.658000 0.250000 15.908000 Sum of Squares 0 0.212000 0.212000 R-Square F-Ratio 0.3102 0.6612 0.0155 0.9868 47.170 201.1 4.717 75.038 Prob > F 0.0017 0.0001** 0.0956 NS 0.0005 Mean Square F-Ratio . 0.053000 0.053000 . Prob > F . 20 O efeito da interação não é significativo, porém, o efeito da regressão do segundo grau é significativo. Este efeito significativo da regressão quadrática indica que nós estamos próximos do ponto de ótimo. Neste ponto, uma análise adicional deve ser feita para localizar o ótimo com mais precisão. Exemplo 2: Achar condições em x1, x2, x3, x4 que maximizam uma resposta. Considerar um planejamento fatorial fracionário 24-1 onde os 4 fatores tem os seguintes níveis: Fator A Fator B Fator C Fator D -1 10 1 25 75 1 2 3 4 +1 15 2 35 85 Considerar um planejamento fatorial 24-1 com gerador I=ABCD, de resolução IV. Os pontos experimentais obtidos foram (saída do SAS): x1 1 12 ,5 2 ,5 3 30 x3 5 x2 2 1,5 0 ,5 4 80 x4 5 21 OBS 1 2 3 4 5 6 7 8 A B C D -1 -1 -1 -1 1 1 1 1 -1 -1 1 1 -1 -1 1 1 -1 1 -1 1 -1 1 -1 1 -1 1 1 -1 1 -1 -1 1 Parameter Estimate INTERCEPT A B C D 63.43750000 1.96250000 2.11250000 -0.31250000 -1.61250000 T for H0: Parameter=0 255.14 7.89 8.50 -1.26 -6.49 y (1) cd bd bc ad ac ab abcd 62.0 57.0 62.2 64.7 61.8 64.5 69.0 66.3 Pr > |T| 0.0001 0.0042 0.0034 0.2978 0.0074 Std Error of Estimate 0.24864215 0.24864215 0.24864215 0.24864215 0.24864215 A equação de regressão fica: yˆ 63.44 1.9625x1 2.1125x2 0.3125x3 1.6125x4 22 O centro do planejamento é dado por: (xi=0, x2=0, x3=0 e x4=0). A partir do centro do planejamento, avançar R unidades na superfície de uma hiperesfera na direção do máximo. Observação: não foi feito o diagnóstico do modelo Algoritmo geral para determinar as coordenadas de um ponto no caminho da máxima inclinação ascendente: assumir que o ponto x1=0, x2=0, ...,xk=0 é a base ou origem. 1. Escolha um tamanho de em uma das variáveis independentes, por exemplo, xj. Geralmente, selecionamos a variável que temos maior conhecimento, ou aquela que tem maior coeficiente de regressão em módulo | ˆ | . j 2. O passo nas demais variáveis é: xi ˆi i 1,2,...,k; i j ˆ j / x j 3. Converter xi das variáveis codificadas para as variáveis naturais. O pesquisador decide usar para o fator 1, 2,5 unidades como tamanho do passo inicial. Usando a relação entre 1 e x1, vimos que 2,5 unidades no fator 1 corresponde a um intervalo (passo), na variável codificada x1 de x1=1. Os passos ao longo da direção da inclinação ascendente são: x2=2,1125/(1,9625. (1,0))=1,0764 x3=-0,3125/(1,9625(1,0))=-0,1592 x4=-1,6125/(1,9625(1,0))=-0,8217. 23 Conversão dos passos nas variáveis codificadas para as variáveis naturais, usamos as relações: x1 1 2 ,5 1 1,0( 2,5) 2,5 x2 2 0,5 2 1,0764(0,5) 0,5382 x3 3 5, 0 3 0,1592(5,0) 0,7962 x4 4 5, 0 4 0,8217(5,0) 4,1085 Ponto Variáveis codificadas experimental X1 X2 X3 X4 Base 0 0 0 0 1,0 1,076 -0,159 -0,822 1,00 1,076 -0,159 -0,822 Base+ Base+2 Base+3 Base+4 1 12,5 2,5 15,0 17,5 20,0 22,5 Variáveis naturais Y 2 3 4 1,5 30 80 --0,538 -0,796 -4,108 2,038 29,204 75,892 --2,576 28,408 71,784 74,0 3,114 27,612 67,676 79,0 3,652 26,816 63,568 77,0 A região de máximo deve estar em torno de 1=20,0 2=3,114 3=27,612 4=67,676 24 * Os pontos em torno de (20,0; 3,114; 27,612; 67,676) decrescem. * Deve ajustar outro modelo de primeira ordem na vizinhança de (20,0; 3,114; 27,612; 67,676) * Fazer análise de variância (o efeito da curvatura é significativo?) 25 Análise de uma superfície de resposta de segunda ordem Quando o pesquisador está próximos da região de ótimo, um modelo que incorpora o efeito de curvatura é indicado. O modelo de segunda ordem é dado por: k k i 1 i 1 y 0 i xi ii xi2 ij xi x j • Como encontrar o ponto ótimo? • Qual a natureza da superfície de reposta? Localização do ponto estacionário Desejamos encontrar os níveis de x1, x2, ...,xk, que maximizam a resposta estimada (predita). Este ponto, se existir, será um conjunto de x1, x2, ...,xk para o qual as derivadas parcias são iguais a zero: yˆ / x1 yˆ / x2 .... yˆ / xk 0. Este ponto, x1.S, x2.S, ...,xk.S é chamado de PONTO ESTACIONÁRIO. Este ponto pode representar um MÁXIMO, MÍNIMO ou PONTO DE SELA. 26 Ponto de máximo (xS) Ponto de mínimo (xS) • • x2 x2 80 60 70 70 60 x1 80 x1 27 28 Determinação do ponto estacionário: uma solução matemática geral. O modelo de segunda ordem escrita na forma matricial fica: yˆ ˆ0 x'b x' Bx Onde: x1 x 2 . x . . xk ˆ1 ˆ 2 . b . . ˆ k ˆ11 ˆ12 / 2, ˆ22 B sim . . . . . . . . ˆ1k / 2 ˆ2k / 2 . . ˆkk b é um vetor (k x 1) dos coeficientes de regressão de primeira ordem e B é uma matriz simétrica (k x k) onde na diagonal têm-se os coeficientes de regressão de segunda ordem e fora da diagonal os coeficientes da interação. As derivadas parciais dos valores preditos da resposta ( y chapéu) com relação aos elementos de x e colocadas iguais a zero são dadas por: yˆ x b 2Bx 0 29 O ponto estacionário é a solução das equações, ou seja, xs 12 B1b O valor predito da variável resposta no ponto estacionário é: yˆ s ˆ0 12 xs' b Demonstração no próximo slide Natureza da superfície de resposta Desejamos saber se o ponto estacionário é um ponto de máximo, mínimo ou ponto de cela. A forma mais direta de se fazer isso é através do gráfico de contornos do modelo de regressão ajustado aos dados. Entretanto, mesmo com poucas variáveis independentes, uma análise mais formal, denominada de Análise Canônica, pode ser útil. Análise canônica (Facilitar a interpretação dos resultados) Considerar uma translação (novo sistema de coordenadas) da superfície de respostas da origem (x1, x2,...,xk)=(0, 0,...,0) para o ponto estacionário xS e então rotacionar os eixos desse sistema até que eles fiquem paralelos aos eixos principais da superfície de resposta ajustada. Veja figura no próximo slide. 30 Opcional: ˆy ˆ 0 x' b x' Bx 1 ' 1 1 ' 1 1 ˆ yˆ 0 b B b b B BB b 2 4 1 1 ' 1 ˆ ˆy 0 b B b b' IB1b 2 4 1 ˆ ˆy 0 b' B1b 2 1 ˆ ˆy 0 x' b 2 31 x2 70 w1 75 80 x1,S xs w2 x1,S x1 A função de respostas em termos das novas variáveis w1, w2,...,wk (forma canônica) é dada por: yˆ yˆ s 1w12 2 w22 ... k wk2 FORMA CANÔNICA DO MODELO Onde os (wi) são as variáveis independentes transformadas e os (i) são constantes. O ys chapéu é a resposta estimada no ponto estacionário. Os (i) são os autovalores ou raízes características da matriz B. 32 OPCIONAL: Redução de uma forma quadrática para uma forma canônica No estudo da forma de uma superfície de resposta e localização das regiões de condições ótimas, é de grande utilidade reduzir uma forma quadrática para uma forma canônica. Resultado: se 1, 2,..., k são raízes características (todas reais) da matriz simétrica real A, então existe uma transformação ortogonal X=Pw, tal que a forma quadrática real Q=X’AX é transformada para a forma canônica 1w12, 2w22,..., kwk2. Isto é, a forma quadrática Q é transformada para uma forma com uma matriz diagonal, onde seus elementos diagonais são as raízes características da matriz A. P : as suas colunas são os autovetores normalizados da matriz A. Q=x’Ax=w’P’ APw=w’w=iwi2 33 Estudo da natureza da superfície de resposta Este estudo pode ser feito considerando o ponto estacionário e os sinais e magnitudes dos (i). Suponha que o ponto estacionário esteja dentro da região de estudo na qual foi ajustado o modelo de segunda ordem. Se todos os valores de (i) são positivos, então, xs é um ponto de resposta mínima; se os (i) são todos negativos, então, xs é um ponto de resposta máxima;se os valores de (i) tem sinais positivos e negativos, então, xs é um ponto de sela. Além disso, a superfície tem inclinação na direção de wi para o qual o valor de |i| é maior. Por exemplo, na figura anterior, xs é um ponto de máximo (todos os (i) são negativos) e |1|> |2|. 34 Exemplo: vamos continuar com a análise do processo químico do exemplo 1 (segunda fase do estudo). Para ajustar um modelo de segunda ordem, o pesquisador decide aumentar o delineamento com pontos adicionais ( o engenheiro usou 4 observações adicionais mais ou menos no mesmo tempo em que executou os 9 tratamentos anteriores. Se passou muito tempo entre as 2 realizações dos tratamentos, então, deve-se usar blocos). Os 4 tratamentos adicionais foram: x1=0, x2=± 1,414 x1=± 1,414, x2=0 O delineamento completo é mostrado na tabela e figura a seguir. Este delineamento denomina-se de DELINEAMENTO CENTRAL COMPOSTO. * Arquivo SAS: chemicalprocesssuperficieresposta2aordem.sas 35 Variáveis naturais 1 2 80 80 90 90 85 85 85 85 85 92,07 77,93 85 85 170 180 170 180 175 175 175 175 175 175 175 182,07 167,93 Delineamento central composto Variáveis codificadas Respostas x1 x2 y1 y2 y3 (produção) (viscosidade) (peso molecular) -1 -1 76.5 62 2940 -1 1 77.0 60 3470 1 -1 78.0 66 3680 1 1 79.5 59 3890 0 0 79.9 72 3480 0 0 80.3 69 3200 0 0 80.0 68 3410 0 0 79.7 70 3290 0 0 79.8 71 3500 1,414 0 78.4 68 3360 -1,414 0 75.6 71 3020 0 1,414 78.5 58 3630 0 -1,414 77.0 57 3150 36 x2 •(0,1,414) (-1,1) • • (-1,414,0) •(1,1) • • (1,414,0) (0,0) (-1,-1) • x1 •(1,-1) •(0,-1,414) 37 Response Surface for Variable YIELD Response Mean Root MSE R-Square Coef. of Variation Regression Linear Quadratic Crossproduct Total Regress Residual Lack of Fit Pure Error Total Error Degrees of Freedom 2 2 1 5 Degrees of Freedom 3 4 7 Type I Sum of Squares 10.042955 17.953749 0.250000 28.246703 Sum of Squares 0.284373 0.212000 0.496373 78.476923 0.266290 0.9827 0.3393 R-Square F-Ratio Prob > F 0.3494 0.6246 0.0087 0.9827 70.814 126.6 3.526 79.669 0.0000 0.0000** 0.1025NS 0.0000 Mean Square F-Ratio Prob > F 0.094791 0.053000 0.070910 1.789 0.2886NS 38 O termo quadrático, comp valor p de 0,0000, foi significativo,portanto, decidimos usar (ajustar) um modelo de segunda ordem para a resposta. Parameter INTERCEPT TIME TEMPERA TIME*TIME TEMPERA*TIME TEMPERA*TEMPERA Degrees of Freedom 1 1 1 1 1 1 Parameter Estimate 79.939955 0.995050 0.515203 -1.376449 0.250000 -1.001336 Standard Error 0.119089 0.094155 0.094155 0.100984 0.133145 0.100984 T for H0: Parameter=0 Prob > |T| 671.3 10.568 5.472 -13.630 1.878 -9.916 0.0000 0.0000 0.0009 0.0000 0.1025 0.0000 Gráfico de contornos e superfície de resposta tridimensional para a variável resposta produção em função do tempo e da temperatura. Observa-se visualmente que o ótimo está próximo de 175oF e 85 minutos e que a resposta neste ponto é um ponto de máximo. Examinando o gráfico de contornos, observa-se que o processo é mais sensível (levemente) à mudanças no tempo de reação do 39 que mudanças na temperatura. 40 41 Determinação da localização do ponto estacionário (máximo). Temos que: 1,376 0,1250 B 0 , 1250 1 , 001 0,995 b 0 , 515 O ponto estacionário é dado por: X1,s x s 12 B 1b 0,7345 0,0917 0,995 0,389 - . 0,0917 1,0096 0,515 0,306 1 2 X2,s Em termos das variáveis naturais, o ponto estacionário é dado por: 0 ,389 1 585 1 86 ,95 87 0 ,306 2 5175 2 176 ,5 42 O valor da resposta estimada no ponto estacionário é: yˆ s ˆ0 12 x 'sb 0,995 ˆys 79.94 12 0,389 0,306 80,21. 0 , 515 ANÁLISE CANÔNICA: Objetivo: caracterizar a superfície de resposta Vamos expressar o modelo ajustado na forma canônica. Primeiro precisamos encontrar os autovalores, 1 e 2. Os autovalores são as raízes do determinante da equação: B I 0 A equação fica: 1,3770 0 ,1250 0 ,1250 1,0013 0 2 2.377786 1.36266393 7 0 As raízes desta equação de segundo grau são: 1=-0,9635 e 2=-1,4143. A forma canônica do modelo ajustado fica: yˆ 80 ,21 0 ,9635w12 1,4143w22 negativos Ponto de máximo 43 Saída do SAS: Canonical Analysis of Response Surface (based on coded data) Critical Value Factor TIME TEMPERA 0.389230 0.305847 Predicted value at stationary point Eigenvalues -0.963498 -1.414287 ** ** 80.212393 Eigenvectors TIME TEMPERA 0.289717 0.957112 Stationary point is a maximum. 0.957112 -0.289717 44 Variáveis canônicas (wi) e as covariáveis (xi) Em muitos problemas é necessário encontrar a relação entre estas duas variáveis. Exemplo: o custo para fazer o experimento no ponto estacionário (1=87 min e 2=176,5°F ) pode ser muito alto (inviabiliza), aí torna-se necessário encontrar um outro ponto (com menor custo) que não signifique muita perda na produção. Para explorar a forma canônica, necessita-se converter os pontos no espaço (w1,w2) para os pontos no espaço (x1,x2). 45 As variáveis x e w são relacionadas por: w M (x x s ) ' onde M é uma matriz ortogonal de dimensão (k x k). As colunas de M são os autovetores normalizados associados com cada (i). Isto é, se mi é a i-ésima coluna de M, então mié a solução para k para o qual : (B i I)mi 0 2 m ji 1. (Normalizado) j 1 Encontrar as colunas de M Exemplo: continuação do exemplo com os fatores tempo e temperatura na produção. 46 Para 1=-0,963499, temos: ( B i I )mi 0 0 ,1250 1,37645 0 ,963499 m11 0 0 ,1250 1,001336 0 ,963499 m21 0 Fazendo-se as operações matriciais, chegamos ao sistema de equações 0 ,412951m11 0 ,1250m21 0 0 ,1250m11 0 ,037837m21 0 Daí, achar m11 e m21, tal que sejam normalizados, isto é, 2 2 2 2 m m m ji 11 21 1. j 1 Estas equações não tem solução única, portanto, vamos dar um valor arbitrário para uma delas, por exemplo, m*21=1, resolver o sistema e, então, normalizar a solução. Seja m*21=1, obtemos m*11=0,302696. Para normalizar esta solução, devemos dividir m*21 e m*11 por * 2 * ( m11 ) ( m21 )2 ( 0 ,302696 )2 ( 1 )2 1,0448085 47 O vetor normalizado fica (a primeira coluna de M): * m11 m11 1,0448085 0 ,2897 m m* 0 ,9571 21 12 1,0448085 Aplicando o mesmo procedimento, agora para 2=-1,414287 obtemos a segunda coluna de M: m12 0 ,9571 m 0 ,2897 22 A relação entre w e x é dada por: w M' ( X X0 ) w1 0 ,2897( x1 0 ,3890 ) 0 ,9571( x2 0 ,3056 ) w2 0 ,9571( x1 0 ,3890 ) 0 ,2897( x2 0 ,3056 ) Se desejamos explorar a superfície de resposta na vizinhança do ponto estacionário, podemos determinar pontos no espaço de (w1, w2) e usar a relação acima para converter estes pontos no espaço de (x1, x2). 48 Exemplo: verificar os efeitos de 4 fatores numa reação química e encontrar as condições que maximizam a resposta. Os fatores são: NH3 (1): amônia (gramas) T (2): temperatura (oC) H2O (3): água (gramas) P (4): pressão do hidrogênio (Psi) Os níveis dos fatores são dados por: Fator Amônia Temperatura Água Pressão -1,4 30,6 222 20 360 -1,0 51 230 100 500 0,0 102 250 300 850 1,0 153 270 500 1200 1,4 173,4 278 580 1340 Os níveis dos fatores, na forma codificada, são dadas por: x1 (1 102) / 51; x2 ( 2 250) / 20; x3 ( 3 300) / 200; x4 ( 4 850) / 350. 49 A matriz de planejamento (Delineamento Central Composto) e os pontos experimentais obtidos são: f a t o r i a l central axial x1 -1 +1 -1 +1 -1 +1 -1 +1 -1 +1 -1 +1 -1 +1 -1 +1 0 -1,4 1,4 0 0 0 0 0 0 x2 -1 -1 +1 +1 -1 -1 +1 +1 -1 -1 +1 +1 -1 -1 +1 +1 0 0 0 -1,4 1,4 0 0 0 0 x3 -1 -1 -1 -1 +1 +1 +1 +1 -1 -1 -1 -1 +1 +1 +1 +1 0 0 0 0 0 -1,4 1,4 0 0 x4 -1 -1 -1 -1 -1 -1 -1 -1 +1 +1 +1 +1 +1 +1 +1 +1 0 0 0 0 0 0 0 -1,4 1,4 y 58.2 23.4 21.9 21.8 14.3 6.3 4.5 21.8 46.7 53.2 23.7 40.3 7.5 13.3 49.3 20.1 32.8 31.1 28.1 17.5 49.7 49.9 34.2 31.1 43.1 4 4 1,414 24=16 2x4=8 50 As estimativas dos parâmetros do modelo se segunda ordem ajustado aos dados valem: (Arquivo SAS: estudocondicoesotimasexperimentosuperficieresposta.sas) Parameter Degrees of Freedom INTERCEPT AMONIA TEMPERA AGUA PRESSAO AMONIA*AMONIA TEMPERA*AMONIA TEMPERA*TEMPERA AGUA*AMONIA AGUA*TEMPERA AGUA*AGUA PRESSAO*AMONIA PRESSAO*TEMPERA PRESSAO*AGUA PRESSAO*PRESSAO 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 Parameter Estimate Standard Error 40.198215 -1.511044 1.284137 -8.738956 4.954819 -6.332399 2.193750 -4.291583 -0.143750 8.006250 0.019642 1.581250 2.806250 0.293750 -2.505869 T for H0: Parameter=0 8.321708 3.151968 3.151968 3.151968 3.151968 5.035510 3.516952 5.035510 3.516952 3.516952 5.035510 3.516952 3.516952 3.516952 5.035510 4.831 -0.479 0.407 -2.773 1.572 -1.258 0.624 -0.852 -0.0409 2.276 0.0039 0.450 0.798 0.0835 -0.498 Prob > |T| 0.0007 0.6420 0.6923 0.0197 0.1470 0.2371 0.5467 0.4140 0.9682 0.0461 0.9970 0.6626 0.4435 0.9351 0.6295 Superfície ajustada 51 Canonical Analysis of Response Surface (based on coded data) Critical Value Factor AMONIA TEMPERA AGUA PRESSAO Ponto estacionário 0.264687 1.033646 0.290578 1.667961 predicted value at stationary point 43.524455 Valor no ponto estacionário Os coeficientes canônicos são as raízes características 1, 2, 3 e 4 da matriz: -6.332 1.097000000 -.07200000000 .7905000000 1.097000000 -4.292 4.003000000 1.403000000 m atrix1 := .020 .1470000000 -.07200000000 4.003000000 .7905000000 1.403000000 .1470000000 -2.506 Matriz B 52 Os coeficientes canônicos (raízes características), valem:-7,55, -6,01, -2,16 e 2,60. Assim, a forma canônica da superfície de resposta é: yˆ 43,53 7,55w12 6,01w22 2,16w32 2,60w42 A relação entre as variáveis naturais e canônicas é dada por: w M ' ( X X0 ) Desenvolvendo, chegamos: w1 0,5977 0,7025 0,3756 w 0,7688 0,4568 0,2858 2 w3 0,2151 0,1374 0,3071 w4 0,0741 0,5282 0,8264 0,0908 x1 0,265 0,3445 x2 1,034 0,9168 x3 0,291 0,1803 x4 1,668 Na equação canônica, observamos que a resposta decresce ao movimentarmos na direção dos eixos w1, w2 , w3,,mas cresce na direção de w4. Se quiser aumento fazemos w1=0, w2=0 , w3=0 Objetivo: encontrar condições operacionais cuja resposta seja alta na região do experimento planejado. 53 Com x3 -1,5 (pois, x3< -1,5 é impossível, sistema sem água), devemos encontrar condições em x1, x2, x3 e x4 que dão valor zero para w1, w2, w3 e vários valores para w4. Resultados na tabela a seguir. Tabela: valores para x1,x2,x3 e x4 com w1=w2=w3=0 e vários valores para w4 w4 x1 x2 x3 x4 1,0 0,339 1,562 1,117 1,848 1,5 0,376 1,826 1,531 1,938 2,0 0,413 2,090 1,944 2,028 -1,0 0,191 0,506 -0,535 1,488 -1,5 0,154 0,242 -0,949 1,398 Sai da faixa experimental. Parar. Por exemplo, para obtermos a primeira coluna da tabela, fazemos: 0,5977x1 0,7025x2 0,7688x 0,4568x 1 2 0,1374x2 0,2151x1 0,5282x2 0,0741x1 0,3756x3 0,0908x4 0,2858x3 0,3445x4 0,3071x3 0,9168x4 0,8264x3 0,1803x4 -2,0 0,117 -0,022 -1,362 1,307 -2,5 0,080 -0,287 -1,775 1,217 Fora de x3-1,5 0,3072 0 0,0182 0 1,6389 0 1,1070 1 Para w4=1 Mx M*(-0,265 -1,034 -0,291 -1,668)’ x= M-1y x=(0,3391 1,562 1,1174 1,8482)’ y=(-0,3072 -0,0182 1,6389 2,10701)’ 54 Interpretações: 1) para valores positivos de w4 os valores de x estão fora da região experimental; 2) condição mais adequada (com w4-2,2): x3=1,50 (sem água), x1 =0,102, x2=0,128 e x4=1,27. 3) qualquer ponto onde w1=0, w2=0, w3=0 e |w4|>2,2 está fora da região experimental pois x3 < -1,5 (valor impossível). 55