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  14080
x3  324
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]
51
 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 B1b
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' IB1b
2
4
1
ˆ
ˆy  0  b' B1b
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
Download

Superficie de Resposta