Introdução à Metodologia de Superfície de Resposta

1. 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, Onde,  representa o ruído branco ou os erros aleatórios observados na reposta y. O valor esperado da resposta, E(y), é dado por: Então, a superfície representada por: é chamada de superfície de resposta.

2. 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)

3. A codificação utilizada foi: Projeto experimental: DELINEAMENTO BOX-BEHNKEN (delineamnetos com 3 níveis). Geometricamente o delineamento fica: ° ° ° ° +1 ° ° ° x3 ° ° +1 ° ° ° x2 -1 -1 ° -1 x1 +1

4. 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.

5. 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: O modelo de segunda ordem é dado por: Efeito da interação Efeito linear Efeito quadrático Estes modelos geralmente funcionam bem para uma região relativamente pequena do espaça dos fatores.

6. 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 de ótimo • • • • • • • • • • • • Região do processo Caminho para a região de ótimo Condição operacional atual (está longe do ótimo

7. 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).

8. O modelo ajustado de primeira ordem é: 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 x1 Figura. Superfície de resposta de primeira ordem e o caminho da inclinação ascendente

9.  Os passos ao longo do caminho são proporcionais aos sinais e grandezas dos coeficientes de regressão { }.  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.

10. 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

11. Codificação (para simplicação dos cálculos: 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 Modelo ajustado de 1a ordem:

12. 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

13. Estimativa do erro experimental [com os pontos centrais (repetição)] Estimativa de mínimos quadrados do coeficiente da interação Soma de quadrados da interação, com 1 grau de liberdade: Teste F para interação: Estudo do efeito quadrático: a análise de variância completa obtemos através do programa estatístico SAS.

14. Response Surface for Variable YIELD Degrees of Type I Sum Regression Freedom of Squares R-Square F-Ratio Prob > F Linear 2 2.825000 0.9410 32.849 0.0033 ** Quadratic 1 0.002722 0.0009 0.0633 0.8137 N.S. Crossproduct 1 0.002500 0.0008 0.0581 0.8213 N.S. Total Regress 4 2.830222 0.9427 16.455 0.0095 ** Degrees of Sum of Residual Freedom Squares Mean Square F-Ratio Prob > F Lack of Fit 0 0 . . . Pure Error 4 0.172000 0.043000 Total Error 4 0.172000 0.043000 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:

15. 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.

16. 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:

18. 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.

20. O modelo de primeira ordem ajustado aos dados é: 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) Degrees of Type I Sum Regression Freedom of Squares R-Square F-Ratio Prob > F Linear 2 5.000000 0.3102 47.170 0.0017 Quadratic 1 10.658000 0.6612 201.1 0.0001** Crossproduct 1 0.250000 0.0155 4.717 0.0956 NS Total Regress 4 15.908000 0.9868 75.038 0.0005 Degrees of Sum of Residual Freedom Squares Mean Square F-Ratio Prob > F Lack of Fit 0 0 . . . Pure Error 4 0.212000 0.053000 Total Error 4 0.212000 0.053000

21. 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, x4que maximizam uma resposta. Considerar um planejamento fatorial fracionário 24-1 onde os 4 fatores tem os seguintes níveis: Considerar um planejamento fatorial 24-1 com gerador I=ABCD, de resolução IV. Os pontos experimentais obtidos foram (saída do SAS):

22. OBS A B C D y 1 -1 -1 -1 -1 (1) 62.0 2 -1 -1 1 1 cd 57.0 3 -1 1 -1 1 bd 62.2 4 -1 1 1 -1 bc 64.7 5 1 -1 -1 1 ad 61.8 6 1 -1 1 -1 ac 64.5 7 1 1 -1 -1 ab 69.0 8 1 1 1 1 abcd 66.3 T for H0: Pr > |T| Std Error of Parameter Estimate Parameter=0 Estimate INTERCEPT 63.43750000 255.14 0.0001 0.24864215 A 1.96250000 7.89 0.0042 0.24864215 B 2.11250000 8.50 0.0034 0.24864215 C -0.31250000 -1.26 0.2978 0.24864215 D -1.61250000 -6.49 0.0074 0.24864215 A equação de regressão fica:

23. 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 . 2. O passo nas demais variáveis é: 3. Converter xidas 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.

24. Conversão dos passos nas variáveis codificadas para as variáveis naturais, usamos as relações: A região de máximo deve estar em torno de 1=20,0 2=3,114 3=27,612 4=67,676

25. * 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?)

26. 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: • 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, ...,xkpara o qual as derivadas parcias são iguais a zero: 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.

27. • • x2 x2 80 60 70 70 60 80 x1 x1 Ponto de máximo (xS) Ponto de mínimo (xS)

29. Determinação do ponto estacionário: uma solução matemática geral. O modelo de segunda ordem escrita na forma matricial fica: Onde: 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:

30. O ponto estacionário é a solução das equações, ou seja, O valor predito da variável resposta no ponto estacionário é: 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.

31. Opcional:

32. x2 w1 x1,S w2 x1 x1,S 70 75 80 xs A função de respostas em termos das novas variáveis w1, w2,...,wk(forma canônica) é dada por: 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.

33. 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

34. 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|.

35. 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

37. x2 • (0,1,414) • • (1,1) (-1,1) • • • x1 (-1,414,0) (1,414,0) (0,0) • • (-1,-1) (1,-1) (0,-1,414) •

38. Response Surface for Variable YIELD Response Mean 78.476923 Root MSE 0.266290 R-Square 0.9827 Coef. of Variation 0.3393 Degrees of Type I Sum Regression Freedom of Squares R-Square F-Ratio Prob > F Linear 2 10.042955 0.3494 70.814 0.0000 Quadratic 2 17.953749 0.6246 126.6 0.0000** Crossproduct 1 0.250000 0.0087 3.526 0.1025NS Total Regress 5 28.246703 0.9827 79.669 0.0000 Degrees of Sum of Residual Freedom Squares Mean Square F-Ratio Prob > F Lack of Fit 3 0.284373 0.094791 1.789 0.2886NS Pure Error 4 0.212000 0.053000 Total Error 7 0.496373 0.070910

39. O termo quadrático, comp valor p de 0,0000, foi significativo,portanto, decidimos usar (ajustar) um modelo de segunda ordem para a resposta. Degrees of Parameter Standard T for H0: Parameter Freedom Estimate Error Parameter=0 Prob > |T| INTERCEPT 1 79.939955 0.119089 671.3 0.0000 TIME 1 0.995050 0.094155 10.568 0.0000 TEMPERA 1 0.515203 0.094155 5.472 0.0009 TIME*TIME 1 -1.376449 0.100984 -13.630 0.0000 TEMPERA*TIME 1 0.250000 0.133145 1.878 0.1025 TEMPERA*TEMPERA 1 -1.001336 0.100984 -9.916 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 que mudanças na temperatura.

42. Determinação da localização do ponto estacionário (máximo). Temos que: O ponto estacionário é dado por: X1,s X2,s Em termos das variáveis naturais, o ponto estacionário é dado por:

43. O valor da resposta estimada no ponto estacionário é: 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: A equação fica: 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: Ponto de máximo negativos

44. Saída do SAS: Canonical Analysis of Response Surface (based on coded data) Critical Value Factor TIME 0.389230 TEMPERA 0.305847 Predicted value at stationary point 80.212393 Eigenvectors Eigenvalues TIME TEMPERA -0.963498 ** 0.289717 0.957112 -1.414287 ** 0.957112 -0.289717 Stationary point is a maximum.

45. 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).

46. para o qual : (Normalizado) As variáveis x e w são relacionadas por: 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 Encontrar as colunas de M Exemplo: continuação do exemplo com os fatores tempo e temperatura na produção.

47. Para 1=-0,963499, temos: Fazendo-se as operações matriciais, chegamos ao sistema de equações Daí, achar m11e m21, tal que sejam normalizados, isto é, 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*21e m*11 por

48. O vetor normalizado fica (a primeira coluna de M): Aplicando o mesmo procedimento, agora para 2=-1,414287 obtemos a segunda coluna de M: A relação entre w e x é dada por: 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).

49. 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: Os níveis dos fatores, na forma codificada, são dadas por:

50. A matriz de planejamento (Delineamento Central Composto) e os pontos experimentais obtidos são: fatorial 24=16 central 2x4=8 axial