SUPERFÍCIE DE RESPOSTA

Superfície de resposta

Objetivo do estudo

Otimização do conteúdo dos componentes de um catalisador

X1: Massas de Cobalto (g)

X2: Massas de tungstênio (g)

X3: Massas de Cerio (g)

E com isso deseja-se otimizar o rendimento de Monóciddo de carbono, que será nossa variáve resposta

y: rendimento de Monóciddo de carbono

Para criação desse planejamento irei utilizar para criar e analisar o Pacote “RSM”.

Chamando nossos dados

## Warning: package 'rsm' was built under R version 4.3.3

##    run.order std.order        Co         W       Ce Block
## 1          1         1  8.000000 1.0000000 3.000000     1
## 2          2         2 12.000000 1.0000000 3.000000     1
## 3          3         3  8.000000 2.0000000 3.000000     1
## 4          4         4 12.000000 2.0000000 3.000000     1
## 5          5         5  8.000000 1.0000000 5.000000     1
## 6          6         6 12.000000 1.0000000 5.000000     1
## 7          7         7  8.000000 2.0000000 5.000000     1
## 8          8         8 12.000000 2.0000000 5.000000     1
## 9          1         1  6.636414 1.5000000 4.000000     2
## 10         2         2 13.363586 1.5000000 4.000000     2
## 11         3         3 10.000000 0.6591036 4.000000     2
## 12         4         4 10.000000 2.3408964 4.000000     2
## 13         5         5 10.000000 1.5000000 2.318207     2
## 14         6         6 10.000000 1.5000000 5.681793     2
## 15         7         7 10.000000 1.5000000 4.000000     2
## 16         8         8 10.000000 1.5000000 4.000000     2
## 17         9         9 10.000000 1.5000000 4.000000     2
## 
## Data are stored in coded form using these coding formulas ...
## x1 ~ (Co - 10)/2
## x2 ~ (W - 1.5)/0.5
## x3 ~ (Ce - 4)/1

Adicionando a respostas

O planejamneto codificado, obteve um alfa de 1,68, e teve 6 pontos axiais que são os experimentos, 9 ao 14, oito pontos fatorias e 3 pontos centrais.

Modelagem

modelo de primeira ordem, completo via função rsm

## 
## Call:
## rsm(formula = y ~ FO(x1, x2, x3), data = plan)
## 
##             Estimate Std. Error  t value Pr(>|t|)    
## (Intercept) 91.39882    0.67302 135.8037   <2e-16 ***
## x1           1.03757    0.75089   1.3818   0.1903    
## x2          -0.74374    0.75089  -0.9905   0.3400    
## x3          -0.21178    0.75089  -0.2820   0.7824    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Multiple R-squared:  0.186,  Adjusted R-squared:  -0.001885 
## F-statistic:  0.99 on 3 and 13 DF,  p-value: 0.428
## 
## Analysis of Variance Table
## 
## Response: y
##                Df  Sum Sq Mean Sq F value   Pr(>F)
## FO(x1, x2, x3)  3  22.869  7.6230    0.99 0.428000
## Residuals      13 100.104  7.7003                 
## Lack of fit    11  99.999  9.0908  173.16 0.005755
## Pure error      2   0.105  0.0525                 
## 
## Direction of steepest ascent (at radius 1):
##         x1         x2         x3 
##  0.8018018 -0.5747435 -0.1636573 
## 
## Corresponding increment in original units:
##         Co          W         Ce 
##  1.6036037 -0.2873718 -0.1636573

Devido à falta de ajuste significativa (p-valor = 0,005), deve-se então utilizar um modelo de ordem maior. O modelo ajustado não apresenta caracteristicas que nos permita estimar valores para a variável resposta. Como não há termos quadráticos e de interação no ajuste do modelo, nota-se que as variáveis x1 x2 e x3 não são significativas.

modelo de segunda ordem, completo via função rsm

## 
## Call:
## rsm(formula = y ~ SO(x1, x2, x3), data = plan)
## 
##              Estimate Std. Error  t value  Pr(>|t|)    
## (Intercept) 94.341870   0.164148 574.7379 < 2.2e-16 ***
## x1           1.037567   0.077085  13.4600 2.932e-06 ***
## x2          -0.743744   0.077085  -9.6483 2.708e-05 ***
## x3          -0.211780   0.077085  -2.7474  0.028613 *  
## x1:x2       -0.460000   0.100717  -4.5673  0.002582 ** 
## x1:x3       -0.975000   0.100717  -9.6806 2.649e-05 ***
## x2:x3       -0.530000   0.100717  -5.2623  0.001170 ** 
## x1^2        -2.331911   0.084843 -27.4849 2.166e-08 ***
## x2^2        -1.578842   0.084843 -18.6089 3.211e-07 ***
## x3^2         0.247261   0.084843   2.9143  0.022521 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Multiple R-squared:  0.9954, Adjusted R-squared:  0.9894 
## F-statistic: 167.6 on 9 and 7 DF,  p-value: 2.464e-07
## 
## Analysis of Variance Table
## 
## Response: y
##                 Df Sum Sq Mean Sq F value    Pr(>F)
## FO(x1, x2, x3)   3 22.869  7.6230  93.937 5.112e-06
## TWI(x1, x2, x3)  3 11.545  3.8483  47.422 5.104e-05
## PQ(x1, x2, x3)   3 87.991 29.3302 361.429 4.909e-08
## Residuals        7  0.568  0.0812                  
## Lack of fit      5  0.463  0.0926   1.764    0.4001
## Pure error       2  0.105  0.0525                  
## 
## Stationary point of response surface:
##         x1         x2         x3 
##  0.1693026 -0.3289331  0.4095174 
## 
## Stationary point in original units:
##        Co         W        Ce 
## 10.338605  1.335533  4.409517 
## 
## Eigenanalysis:
## eigen() decomposition
## $values
## [1]  0.3614754 -1.5213599 -2.5036078
## 
## $vectors
##          [,1]        [,2]      [,3]
## x1  0.1675212  0.30486530 0.9375467
## x2  0.1138866 -0.95060326 0.2887616
## x3 -0.9792683 -0.05840032 0.1939663

Obtém-se agora falta de ajuste não significativa (p-valor = 0,40), ou seja, o modelo de segunda ordem se ajustou bem aos dados. Observa-se também que o ponto de estacionariedade do modelo ajustado está em (0.1693026, -0.3289331, 0.4095174 ), que é ponto de máximo.

Para construir as curvas de níveis e a superfície de resposta, pode ser empregada a função lm (linear model) ou a função rsm (response surface methodology).

condição de otimalidade - Análise canonica

## $xs
##         x1         x2         x3 
##  0.1693026 -0.3289331  0.4095174 
## 
## $eigen
## eigen() decomposition
## $values
## [1]  0.3614754 -1.5213599 -2.5036078
## 
## $vectors
##          [,1]        [,2]      [,3]
## x1  0.1675212  0.30486530 0.9375467
## x2  0.1138866 -0.95060326 0.2887616
## x3 -0.9792683 -0.05840032 0.1939663

Checando os Pressupostos Da Metodologia de superficie de resposta, que são eles

1° Normalidade dos Resíduos

2° indepêndência dos Resíduos

3° Variância constante

Agora iremos fazer esse checagem

Teste de normalidade dos Resíduos

## 
##  Shapiro-Wilk normality test
## 
## data:  res.rsm1$residuals
## W = 0.94687, p-value = 0.4089

indepêndência dos Resíduos

## Warning: package 'lmtest' was built under R version 4.3.3

## Carregando pacotes exigidos: zoo

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

## 
##  Breusch-Godfrey test for serial correlation of order up to 6
## 
## data:  res.rsm1
## LM test = 12.222, df = 6, p-value = 0.0572

Variância constante - Homodasticidade

## 
##  studentized Breusch-Pagan test
## 
## data:  res.rsm1
## BP = 11.145, df = 9, p-value = 0.2659

Graficamente

Gráficos de contorno e Gráficos de superficie

Visualizando os Gráficos separadamente

Vendo x1(Coblato) e x2(tungstênio), possivelmente teria um ponto estacinário de máximo uma superfície com convexidade para baixo os alto valores negativos.

Mas como temos uma variável a mais na hora que foi plotado x1, x3 e x2 e x3 podemos ver que a superfície parece ter aparência de um sela de cavalo, então o ponto estacionário proxímo ao centro, mas isso na verdade aqui corte da superfície de resposta, e não conseguimos enxergar no \(R^4\), esse ponto ele não maximiza a função a gente vê que máximo está indo pra região branca, como vemos nos gráficos de contorno. O ponto estacionário é simplesmente um ponto onde você tem um ponto de mudança de convexidade.

pontos estacionários e de Previsão

Matrizes b e B

##         x1         x2         x3 
##  1.0375673 -0.7437437 -0.2117799

##           x1        x2         x3
## x1 -2.331911 -0.230000 -0.4875000
## x2 -0.230000 -1.578842 -0.2650000
## x3 -0.487500 -0.265000  0.2472609

Ponto estacionarios

##         x1         x2         x3 
##  0.1693026 -0.3289331  0.4095174

##          [,1]
## x1  0.1693026
## x2 -0.3289331
## x3  0.4095174

Deconposição em valores singulares

## eigen() decomposition
## $values
## [1]  0.3614754 -1.5213599 -2.5036078
## 
## $vectors
##            [,1]        [,2]      [,3]
## [1,]  0.1675212  0.30486530 0.9375467
## [2,]  0.1138866 -0.95060326 0.2887616
## [3,] -0.9792683 -0.05840032 0.1939663

## eigen() decomposition
## $values
## [1]  0.3614754 -1.5213599 -2.5036078
## 
## $vectors
##          [,1]        [,2]      [,3]
## x1  0.1675212  0.30486530 0.9375467
## x2  0.1138866 -0.95060326 0.2887616
## x3 -0.9792683 -0.05840032 0.1939663

Com isso vamos fazer a otimização não linear ou restrita ou análise rígida o pacote “RSM”, tm esse comando “steepest”, que é do método de escalada acentuada, e esse comando providência o passo de escalada acentuada ou descida acentuada para um modelo de superfície de resposta , e ele é utilizado para modelos de primeira e segunda ordem.

Realizando a otimização

Análise rígida

## Path of steepest ascent from ridge analysis:

##    dist    x1     x2     x3 |     Co      W    Ce |   yhat
## 1   0.0 0.000  0.000  0.000 | 10.000 1.5000 4.000 | 94.342
## 2   0.1 0.076 -0.059 -0.027 | 10.152 1.4705 3.973 | 94.455
## 3   0.2 0.140 -0.111 -0.090 | 10.280 1.4445 3.910 | 94.540
## 4   0.3 0.188 -0.136 -0.190 | 10.376 1.4320 3.810 | 94.609
## 5   0.4 0.221 -0.142 -0.302 | 10.442 1.4290 3.698 | 94.674
## 6   0.5 0.247 -0.140 -0.412 | 10.494 1.4300 3.588 | 94.743
## 7   0.6 0.270 -0.134 -0.519 | 10.540 1.4330 3.481 | 94.816
## 8   0.7 0.291 -0.126 -0.624 | 10.582 1.4370 3.376 | 94.896
## 9   0.8 0.311 -0.117 -0.727 | 10.622 1.4415 3.273 | 94.981
## 10  0.9 0.331 -0.108 -0.830 | 10.662 1.4460 3.170 | 95.075
## 11  1.0 0.350 -0.098 -0.932 | 10.700 1.4510 3.068 | 95.175
## 12  1.1 0.368 -0.088 -1.033 | 10.736 1.4560 2.967 | 95.281
## 13  1.2 0.386 -0.078 -1.134 | 10.772 1.4610 2.866 | 95.395
## 14  1.3 0.404 -0.067 -1.234 | 10.808 1.4665 2.766 | 95.516
## 15  1.4 0.422 -0.056 -1.334 | 10.844 1.4720 2.666 | 95.644
## 16  1.5 0.440 -0.046 -1.433 | 10.880 1.4770 2.567 | 95.778
## 17  1.6 0.457 -0.035 -1.533 | 10.914 1.4825 2.467 | 95.921

Com isso temos os níveis codificados da minha resposta no raio igual Raiz 4 de \(2^k\), e os níveis decodificados e o ótimo restrito que foi \(95,921\), que seria o nosso ótimo restrito, então quando tenho cobalto igual a \(10,91\), tungstênio igual a \(1,48\) cerio igual a \(2,47\), tudo isso em gramas vamos ter o rendimento de monóxido de carbono igual a \(95,92%\) aproximadamente

Graficamente

Aqui podemos ver gráficos de contorno mostrando os pontos do planejamento, e podemos ver neles os pontos ótimos restritos. em coloração magenta

## Warning in xy.coords(x, y): NAs introduzidos por coerção

## Warning in xy.coords(x, y): NAs introduzidos por coerção

## Warning in xy.coords(x, y): NAs introduzidos por coerção

Conclusão

Com base nos fatos apresentados, conclui-se que o modelo linear não apresentou o melhor desempenho em comparação ao modelo de segunda ordem. Esse resultado sugere que os dados podem não possuir uma relação linear. Assim, para uma modelagem e otimização mais robusta, é necessário recorrer a um modelo de segunda ordem.

O modelo de segunda ordem, por sua vez, é mais adequado para os dados. Isso é evidenciado pelo fato de o \(R^2\) ajustado ter superado o valor de 0,98. o modelo foi capaz de identificar pontos ótimos que maximizam a função, cumprindo seu objetivo principal. Dessa forma, ele se apresenta como a melhor opção para futuras predições.

Com isso temos que a otimização restrita é muito importante, pois em trabalhos que temos \(k\) a partir de 3 fatores, porque nesses casos começaremos a ter os altovalores mistos, o que leva a ter superfície de resposta no formato de sela como visto aqui, e nesses casos temos que usar a otimização restrita.