Metodologia de Superfície de Resposta
Marília Canabarro Zordan
Sabrina Letícia Couto da Silva
MAT 02014 - Planejamento de Experimentos II
MÉTODOS DE SUPERFÍCIE DE
RESPOSTA
-Conjunto de técnicas estatísticas usadas
para analisar problemas com variáveis
independentes contínuas em relação a
variável aleatória.
- O objetivo é buscar a combinação dos
fatores que otimizam a resposta.
-É usado quando:
• Fixa-se um fator A e varia-se o outro B.
• Fixa-se B onde a resposta é máxima e varia-se
A.
• Variando os dois fatores conjuntamente.
- A modelagem é usada para estabelecer a
relação entre a variável resposta e a(s) variável
(is) independente(s), pois em geral esta não é
conhecida e precisa de aproximação.
O procedimento seqüencial
na MSR é
encontrar os resultados da ANOVA e verificar
se:
* São distante do ótimo, com pouca curvatura
→ Modelo de 1ª ordem
Leva de forma rápida e
vizinhança do valor ótimo.
eficiente
até
a
* A região ótima é encontrada → Modelo de 2ª
ordem e análise do ótimo
Função e Representação Gráfica
• y é a resposta (variável dependente, p.e, produção de
um objeto)
• x1 e x2 são os fatores (p.e, temperatura e pressão)
• y = f(x1,x2) + erro
• E(y) = η → η = f(x1,x2) = Superfície de Resposta
• y = η + erro
Uma possível representação gráfica para a
SR seria:
x1 e x2 aparecem no
plano
y aparece no eixo
perpendicular
Um representação alternativa é o gráfico de
contornos.
Desenham-se linhas de igual resposta de um
gráfico cujas coordenadas representam os níveis dos
fatores:
Para três ou mais fatores é mais comum a
representação gráfica da superfície através do
gráfico de contornos no espaço bidimensional,
fixando um ou mais fatores no nível ótimo.
Exemplo:
Y = produção de uma reação química
min)
A = Tempo de Reação (60,90,120,150,180
B = Temperatura (210,220,230,240,250°C)
Procedimento Clássico
• Fixando B = 225°C → 130 min e 75g
• Fixando A = 130min → 225°C e 75g
E como os fatores se comportam conjuntamente?
É preciso conduzir o experimento variando
ambos os fatores simultaneamente e então para
representação
gráfica
usar
um
gráfico
de
contornos.
Continuando o exemplo...
A = 65 minutos
B = 225°C
Y = 91g
Conclusão
do
Exemplo: O procedimento de
investigar uma variável independente (fator) de
cada vez falha, pois ele assume que o máximo
valor produzido por um fator é independente do
máximo produzido pelo outro fator, o que nem
sempre acontece.
Aproximações
Como nem sempre a relação entre fatores e
resposta é conhecida, o primeiro passo em MSR é
encontrar uma aproximação adequada para a
relação entre eles, dentro de uma faixa limitada do
“espaço de fatores”.
Se a resposta é bem modelada por uma função
linear das var. independentes:
y = β0 + β1x1 + ... + βkxk + ε
A análise da SR é então feita em termos da SR
ajustada.
Modelos de Primeira Ordem
Método da Máxima Inclinação Ascendente
1 - Delineamentos Exploratórios de Tratamentos
2 – Fatoriais da série 2k
Onde,
K = n° de fatores completos ou fracionados
2 = n° de níveis
3 – Verificar se há acréscimo (decréscimo) na
resposta ou não, com a presença dos fatores
• Condições iniciais estão afastadas daquelas que
otimizam a resposta → Modelo de 1ª Ordem
• O objetivo é mover o experimento rapidamente
para a vizinhança geral do ótimo utilizando um
procedimento
experimental
simples,
rápido,
econômico e eficiente.
• MMIA procura a MÁXIMA inclinação ascendente
• MMID procura a MÍNIMA inclinação descendente
• O gráfico de contornos da SR de 1ª ordem são uma
série de linhas paralelas. A direção da MIA é a
direção
na
qual
y
estimado
cresce
mais
rapidamente.
X2
X1
Usualmente, toma-se como caminho da MIA a
linha a partir do centro da região de interesse e os
“passos” ao longo do caminho são determinados
pela experiência do pesquisador.
Utiliza-se um conjunto de tratamentos em
torno do ponto inicial e estima-se por Mínimos
Quadrados as inclinações βi. A partir das magnitudes
e sinais destas inclinações, calcula-se a direção da
MIA.
• Experimentos são conduzidos ao longo do caminho
da MIA até que nenhum incremento na resposta seja
observado. E repetir a seqüência novamente.
•Eventualmente chega-se a vizinhança do valor
ótimo e isto será indicado pela falta de ajuste do
modelo de 1ª ordem.
• A aproximação por um plano se torna insatisfatória
pelo fato dos efeitos de ordens mais elevadas,
particularmente os de 2ª ordem (quadrático e de
interação linear), se tornarem relativamente mais
importantes. Nesse caso usa-se Modelo Quadráticos.
EXEMPLO (produção do reagente):
O processo normal é operado com um tempo
de 35min e uma temperatura de 155°F, que resulta
numa produção de 40%.
Como a região ótima é desconhecida ajustase um modelo de 1ª ordem por MMIA.
A região experimental será (30,40) min para
o
tempo
de
reação
e
(150,160)°F
para
temperatura.
As variáveis
recodificadas para
cálculos.
independentes podem ser
(-1,1) para simplificar os
Se ξ1 e ξ2 representam as variáveis naturais
para tempo e temperatura respectivamente, os
valores são codificados:
135
x1 
5
1155
x2 
5
“Tamanho do Passo”
Média do Intervalo
O delineamento a ser utilizado nesse
experimento consiste de um Fatorial 22 aumentado
por CINCO PONTOS (Tratamentos) CENTRAIS.
Com a utilização dos pontos centrais é
possível estimar o erro experimental e testar a
adequabilidade (“lack of fit”) do modelo de 1ª
ordem. Além disso, nesse caso, o tratamento
central representa as condições de operação
normalmente empregadas.
Deve-se ajustar aos dados um modelo de
primeira ordem, ou seja, a equação de regressão
yˆ  ˆ oˆ1x1  ˆ 2x 2
Dados
Variáveis Naturais
ξ1
30
30
40
40
35
35
35
35
35
ξ2
150
160
150
160
155
155
155
155
155
Variáveis Codificadas
x1
-1
-1
1
1
0
0
0
0
0
x2
-1
1
-1
1
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
A equação de regressão encontrada é:
yˆ  40,44  0,775x1  0,325x 2
A SQTotal é encontrada da mesma forma de
qualquer outra análise, mas é particionada em
SQRegressão e SQResíduos. A SQErro com 4GL
também é obtida de forma tradicional, porém só
com os valores dos pontos centrais.
Assim, o modelo de primeira ordem assume
que os fatores possuem um efeito aditivo sobre a
resposta.
A interação entre os fatores pode ser obtida
adicionando o termo x1x2 e é medida pelo
coeficiente β12. A estimativa é obtida (considerando
as variáveis codificadas) por:
ˆ 12 
1
( x1y1  x1y 2  x 2 y3  x 2 y 4 )
n de trat
SQInt 
( y extremos   y meio )
n de trat
2
Para o exemplo...
ˆ  1 1* 39,3  (1* 40)  (1* 40,9)  1* 41,55   0,025
12
4
[(39,3  41,5)  (40  40,5)]2
SQInt 
 0,0025
4
1 GL
A estatística da falta de ajuste é:
SQInt 0,0025
F

 0,058
QME 0,043
F(1,4)=7,71 a 5%
Logo, NS.
Outra verificação da adequabilidade do
modelo é obtida pela comparação da resposta
média dos quatro pontos do fatorial (40,425), com
a resposta média do centro do delineamento
(40,46).
Se β11 e β22 são os coeficientes dos termos
quadráticos x 2 e x 2 , então a diferença das médias
1
2
é uma estimativa de β11+ β22.
ˆ11ˆ12  y1 y2  40,425 40,46  0,035
n f n c ( y1  y 2 ) 2 (4 * 5)(0,035) 2
SQQuad 

 0,0027
nf  nc
45
SQQuadrático 0,0027
F

 0,0663
QME
0,0430
Efeito quadrático NS a 5% de significância.
Causas da Variação
SQ
GL
QM
F
Regressão
2,8250
2
1,4125
47,83*
Resíduo
0,1772
6
0,0295
Interação
0,0025
1
0,0025
0,058
Quadr. Puro
0,0027
1
0,027
0,063
Erro
0,1720
4
0,0430
Total
3,0022
8
* Significante a 1%
Assim, não existe nenhuma razão para
questionar o modelo de primeira ordem. Os
próximos passos da MIA devem seguir.
Os fatores devem variar nas proporções das
estimativas dos Betas, ou seja, para cada 0,775
unidades acrescidas em x1, x2 deve aumentar 0,325
unidades. Destransformando-a, direção da MIA é
definida por:
5x0,775 = 3,875 minutos → 5x0,325=1,625°
O caminho principal a partir da origem (0,0)
nesta direção pode ser obtido por um incremento
conveniente a um dos fatores (p.e, 5 minutos para
tempo). As mudanças proporcionais do outro fator:
0,325/0,775=0,42.
Assim,
estas
quantidades
devem
ser
acrescidas aos níveis base dando origem ao
caminho da MIA.
O engenheiro realizou os ensaios de acordo
com esse caminho e observou a produção nesses
pontos até que um decréscimo na resposta foi
notado.
Os resultados estão na tabela a seguir:
Variáveis Codificadas
Nível Base
X1
X2
ξ1
ξ2
0
0
35
155
5
5
3,875
1,625
Unidade
Unidade*
Variação do Nível
Variáveis Naturais
Resposta
y
ˆ i
1,00
0,42
5
2
Ensaio 1
1
0,42
40
157
41,0
Ensaio 2
2
0,84
45
159
42,9
Ensaio 3
3
1,26
50
151
47,1
Ensaio 4
4
1,68
55
163
49,7
Ensaio 5
5
2,10
60
165
53,8
Ensaio 6
6
2,52
65
167
59,9
Ensaio 7
7
2,96
70
169
65,0
Ensaio 8
8
3,36
75
171
70,4
Ensaio 9
9
3,78
80
173
77,6
Ensaio 10
10
4,20
85
175
80,3
Ensaio 11
11
4,62
90
177
76,2
Ensaio 12
12
5,04
95
179
75,1
Incrementos na resposta são observados até o
10° passo; depois há um decréscimo na produção.
Portanto, outro modelo de primeira ordem
pode ser ajustado em torno do ponto (85,175). A
região de exploração para tempo seria (80,90) e de
temperatura (170,180).
Codificam-se os níveis das variáveis tempo e
temperatura
novamente
como
(-1,1)
e
um
delineamento fatorial 22 acrescido de 5 pts centrais
será utilizado.
Assim repete-se todo o processo
resultados são analisados a seguir:
e
os
Dados para ajuste do segundo modelo de 1ª Ordem:
Variáveis
Naturais
Variáveis
Codificadas
1
2
X1
X2
y
80
170
-1
-1
76.5
80
180
-1
1
77.0
90
170
1
-1
78.0
90
180
1
1
79.5
85
175
0
0
79.9
85
175
0
0
80.3
85
175
0
0
80.0
85
175
0
0
79.7
85
175
0
0
79.8
x1 
1  85
5
e
Resposta
x2 
2  175
5
O modelo de 1ª Ordem ajustado aos dados codificados é
dado por:
Ŷ = 78.97 + 1.00 X1 + 0.50 X2
Causas da Variação
SQ
GL
QM
F
Regressão
5,0000
2
2,5000
47,17*
Resíduo
11,1200
6
Interação
0,2500
1
0,2500
4,72
Quadr. Puro
10,6580
1
10,6580
201,09 *
Erro Puro
0,2120
4
0,0530
Total
16,1200
8
* Significativo a 1%
- Pela tabela de ANOVA, o componente do termo quadrático
puro foi significativo, isso implica que o modelo de 1ª Ordem
não é uma aproximação adequada;
- Essa curvatura na real superfície pode indicar que se está
próximo do ótimo. Assim, análises adicionais devem ser feitas
para localizar o ótimo com mais precisão.
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 para 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
ˆ j / x j
i  1,2,...,k; i  j
3. Converter xi das variáveis codificadas para as variáveis
naturais.
Análise de Modelos de Segunda Ordem
- Quando o pesquisador está próximo da região de ótimo, um
modelo que incorpora o efeito de curvatura é indicado.
O modelo de 2ª Ordem é dado por:
k
k
y  0   i xi    x   ij xi x j  
i 1
i 1
2
ii i
• Como encontrar o ponto ótimo (estacionário)?
• 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 parciais são iguais a zero:
yˆ / x1  yˆ / x2  ....  yˆ / xk  0.
Esse ponto, x1.S, x2.S, ...,xk.S é o dito PONTO ESTACIONÁRIO.
Este ponto pode representar um MÁXIMO, MÍNIMO ou
PONTO DE SELA.
Ponto de máximo (xS)
x2
•
Ponto de mínimo (xS)
x2
80
•
60
70
70
60
x1
80
x1
Determinação do ponto estacionário:
Solução matemática geral. O modelo de 2ª Ordem escrito na forma matricial
fica:
yˆ  ˆ0  x'b  x' Bx
 x1 
x 
 2
.
x 
.
.
 
 xk 
 ˆ1 
ˆ
2 
 .
b 
 .
 .
 
 ˆk 
 ˆ11 ˆ12 / 2,

ˆ22


B



 sim .
.
.
.
.
.
.
.
ˆ1k / 2

ˆ2k / 2
.
.
ˆkk






Onde: b é um vetor (k x 1) dos coeficientes de regressão de 1ª ordem e B é
uma matriz simétrica (k x k) onde na diagonal têm-se os coeficientes de
regressão de 2ª ordem e fora da diagonal os coeficientes de 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
O ponto estacionário é a solução da equação anterior, ou seja,
xs   12 B1b
O valor predito da variável resposta no ponto estacionário é:
yˆ s  ˆ0  12 xs' b
Demonstração:
yˆ  ˆ 0  x' b  x' Bx
1
1
yˆ  ˆ 0  b' B1b  b' B1BB 1b
2
4
ˆy  ˆ 0  1 b' B1b  1 b' IB1b
2
4
1
yˆ  ˆ 0  b' B1b
2
1
yˆ  ˆ 0  x' b
2
Natureza da superfície de resposta
Desejamos saber se o ponto estacionário é um ponto de máximo,
mínimo ou ponto de cela.
-Forma mais direta: gráfico de contorno do modelo de regressão
ajustado aos dados.
-Mesmo com poucas variáveis independentes, uma análise mais
formal, denominada de Análise Canônica, pode ser útil.
Análise Canônica (Facilita a interpretação dos resultados!!!)
Considere uma translação (novo sistema de coordenadas)
da superfície de resposta da origem (x1, x2,...,xk)=(0, 0,...,0) para
o ponto estacionário xS e então rotacione os eixos desse sistema
até que eles fiquem paralelos aos eixos principais da superfície de
resposta ajustada.
Na figura no próximo slide isso pode ser visualizado...
x2
w1
x1,S
w2
x1,S
x1
A função de resposta 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 indep. transformadas e os i são
constantes. O ŷs é a resposta estimada no ponto estacionário xS. Os i
são os autovalores ou raízes características da matriz B.
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 na qual foi
ajustado o modelo de 2ª 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|.
Exemplo 2: vamos continuar com a análise do processo químico do
ex 1 (2ª fase do estudo).
Para ajustar um modelo de 2ª ordem, o pesquisador decide
aumentar o delineamento com pontos adicionais (o engenheiro usou 4
OBS adicionais mais ou menos no mesmo tempo em que executou os
9 tratamentos anteriores).
Os 4 tratamentos adicionais foram:
(x1=0; x2=± 1,414)
(x1=± 1,414; x2=0)
Este delineamento denomina-se de DELINEAMENTO
CENTRAL COMPOSTO.
VARIÁVEL RESPOSTA
Variáveis naturais
2
1
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
Respostas
Variáveis codificadas
y3
y2
y1
x2
x1
(produção) (viscosidade) (peso molecular)
2940
62
76.5
-1
-1
3470
60
77.0
1
-1
3680
66
78.0
-1
1
3890
59
79.5
1
1
3480
72
79.9
0
0
3200
69
80.3
0
0
3410
68
80.0
0
0
3290
70
79.7
0
0
3500
71
79.8
0
0
3360
68
78.4
0
1,414
3020
71
75.6
0
-1,414
3630
58
78.5
1,414
0
3150
57
77.0
-1,414
0
x2
•(0,1,414)
(-1,1) •
•
(-1,414,0)
•(1,1)
•
•
(1,414,0)
(0,0)
(-1,-1) •
•(1,-1)
•(0,-1,414)
Delineamento Central Composto para o exemplo 1 (processo químico)
x1
Response Surface for Variable PRODUCAO
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
0.3494
0.6246
0.0087
0.9827
70.814
126.6
3.526
79.669
Prob > F
Mean Square
F-Ratio
0.094791
0.053000
0.070910
1.789
0.0000
0.0000**
0.1025 NS
0.0000
Prob > F
0.2886
NS
Análise no DOE do SAS 8.0:
Efeito
Quadrático
significativo a 1%.
O modelo de 2ª
Ordem será ajustado!
Gráfico de Contorno
Pelo gráfico de contornos, observa-se que o processo é mais
sensível (levemente) à mudanças no tempo de reação do que para
mudanças na temperatura.
Superfície de Resposta
Ponto Ótimo  (85 min; 175º F)
Observa-se que o ótimo está próximo de 175oF e 85 min e que
a resposta neste ponto é um ponto de máximo.
Determinação da localização do ponto estacionário (máximo).
Temos que:
0,995
b

0
,
515


 1,376 0,1250
B

0
,
1250

1
,
001


O ponto estacionário é dado por:
1
xs   B b
1
2
X1,s
  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
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. 1º 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
Ambos negativos
Ponto de máximo