Experimento Catapulta: RSM BOX-BEHNKEN DESIGN: Modelagem e Otimização pela Linguagem R
library(rsm)
library(car)
library(lmtest)
library(nortest)
library(Rsolnp)
library(NlcOptim)
DOE<-bbd(k=~x1+x2+x3,
block = F,
n0=5,
randomize = F,
coding = list(x1~(ang-4)/2,
x2~(compr-4)/1,
x3~(alt-3)/1))
DOE run.order std.order ang compr alt
1 1 1 2 3 3
2 2 2 6 3 3
3 3 3 2 5 3
4 4 4 6 5 3
5 5 5 2 4 2
6 6 6 6 4 2
7 7 7 2 4 4
8 8 8 6 4 4
9 9 9 4 3 2
10 10 10 4 5 2
11 11 11 4 3 4
12 12 12 4 5 4
13 13 13 4 4 3
14 14 14 4 4 3
15 15 15 4 4 3
16 16 16 4 4 3
17 17 17 4 4 3
Data are stored in coded form using these coding formulas ...
x1 ~ (ang - 4)/2
x2 ~ (compr - 4)/1
x3 ~ (alt - 3)/1DOE$y<-(c(2.71,2.75,4.17,4.12,2.42,2.61,4.36,
4.82,2.65,3.43,3.68,5.62,3.43,3.45,
3.82,3.46,3.68))
DOE run.order std.order ang compr alt y
1 1 1 2 3 3 2.71
2 2 2 6 3 3 2.75
3 3 3 2 5 3 4.17
4 4 4 6 5 3 4.12
5 5 5 2 4 2 2.42
6 6 6 6 4 2 2.61
7 7 7 2 4 4 4.36
8 8 8 6 4 4 4.82
9 9 9 4 3 2 2.65
10 10 10 4 5 2 3.43
11 11 11 4 3 4 3.68
12 12 12 4 5 4 5.62
13 13 13 4 4 3 3.43
14 14 14 4 4 3 3.45
15 15 15 4 4 3 3.82
16 16 16 4 4 3 3.46
17 17 17 4 4 3 3.68
Data are stored in coded form using these coding formulas ...
x1 ~ (ang - 4)/2
x2 ~ (compr - 4)/1
x3 ~ (alt - 3)/1funcao<- rsm(y~ SO(x1,x2,x3),data = DOE)
summary(funcao)
Call:
rsm(formula = y ~ SO(x1, x2, x3), data = DOE)
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.568000 0.090199 39.5570 1.719e-09 ***
x1 0.080000 0.071309 1.1219 0.29892
x2 0.693750 0.071309 9.7288 2.564e-05 ***
x3 0.921250 0.071309 12.9192 3.869e-06 ***
x1:x2 -0.022500 0.100846 -0.2231 0.82982
x1:x3 0.067500 0.100846 0.6693 0.52473
x2:x3 0.290000 0.100846 2.8757 0.02380 *
x1^2 -0.211500 0.098292 -2.1518 0.06844 .
x2^2 0.081000 0.098292 0.8241 0.43707
x3^2 0.196000 0.098292 1.9941 0.08637 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Multiple R-squared: 0.9756, Adjusted R-squared: 0.9443
F-statistic: 31.16 on 9 and 7 DF, p-value: 7.854e-05
Analysis of Variance Table
Response: y
Df Sum Sq Mean Sq F value Pr(>F)
FO(x1, x2, x3) 3 10.6911 3.5637 87.6050 6.482e-06
TWI(x1, x2, x3) 3 0.3566 0.1189 2.9225 0.1096
PQ(x1, x2, x3) 3 0.3598 0.1199 2.9486 0.1079
Residuals 7 0.2848 0.0407
Lack of fit 3 0.1641 0.0547 1.8128 0.2846
Pure error 4 0.1207 0.0302
Stationary point of response surface:
x1 x2 x3
-0.2851951 0.6256368 -2.7638622
Stationary point in original units:
ang compr alt
3.4296098 4.6256368 0.2361378
Eigenanalysis:
eigen() decomposition
$values
[1] 0.29540968 -0.01347974 -0.21642994
$vectors
[,1] [,2] [,3]
x1 0.04278373 0.1402676 0.98918884
x2 0.55814419 -0.8245405 0.09277986
x3 0.82864022 0.5481405 -0.11356645funcao<-rsm(y~ FO(x1,x2,x3)+
TWI(x2,x3)+PQ(x1,x3),data = DOE)
summary(funcao)
Call:
rsm(formula = y ~ FO(x1, x2, x3) + TWI(x2, x3) + PQ(x1, x3),
data = DOE)
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.602105 0.072471 49.7040 2.625e-13 ***
x1 0.080000 0.064482 1.2407 0.243048
x2 0.693750 0.064482 10.7589 8.099e-07 ***
x3 0.921250 0.064482 14.2870 5.578e-08 ***
x2:x3 0.290000 0.091191 3.1801 0.009817 **
x1^2 -0.207237 0.088759 -2.3348 0.041705 *
x3^2 0.200263 0.088759 2.2563 0.047670 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Multiple R-squared: 0.9716, Adjusted R-squared: 0.9545
F-statistic: 56.92 on 6 and 10 DF, p-value: 3.73e-07
Analysis of Variance Table
Response: y
Df Sum Sq Mean Sq F value Pr(>F)
FO(x1, x2, x3) 3 10.6911 3.5637 107.1372 6.685e-08
TWI(x2, x3) 1 0.3364 0.3364 10.1133 0.009817
PQ(x1, x3) 2 0.3322 0.1661 4.9939 0.031346
Residuals 10 0.3326 0.0333
Lack of fit 6 0.2120 0.0353 1.1709 0.459702
Pure error 4 0.1207 0.0302
Stationary point of response surface:
x1 x2 x3
0.1930159 0.1272608 -2.3922414
Stationary point in original units:
ang compr alt
4.3860317 4.1272608 0.6077586
Eigenanalysis:
eigen() decomposition
$values
[1] 0.27634546 -0.07608231 -0.20723684
$vectors
[,1] [,2] [,3]
x1 0.0000000 0.0000000 1
x2 0.4646295 0.8855052 0
x3 0.8855052 -0.4646295 0#======== Análise dos Resíduos ===============
curve(dnorm(x,mean(funcao$residuals),sd(funcao$residuals)),
from = min(funcao$residuals),to=max(funcao$residuals),
col="blue",lwd=4,main="Curva Normal dos Resíduos",
xlab = "Observações",ylab = "Densidade")
qqnorm(funcao$residuals,lwd=3,col="blue")
qqline(funcao$residuals,lwd=3)
library(nortest)
ad.test(funcao$residuals)
Anderson-Darling normality test
data: funcao$residuals
A = 0.6899, p-value = 0.05845residualPlot(funcao)
bptest(funcao)
studentized Breusch-Pagan test
data: funcao
BP = 1.3968, df = 6, p-value = 0.9661durbinWatsonTest(funcao) lag Autocorrelation D-W Statistic p-value
1 0.1141991 1.729611 0.43
Alternative hypothesis: rho != 0dwt(funcao) lag Autocorrelation D-W Statistic p-value
1 0.1141991 1.729611 0.444
Alternative hypothesis: rho != 0#Otimização
obj=function(x){a<-predict(funcao,newdata = data.frame(x1=x[1],
x2=x[2],
x3=x[3]))
return(-a)}
x0<-c(0,0,0)
optim(par = x0,fn=obj,lower = c(-1,-1,-1),upper = c(1,1,1),
method = c("L-BFGS-B"))$par
[1] 0.1930159 1.0000000 1.0000000
$value
[1] -5.715089
$counts
function gradient
7 7
$convergence
[1] 0
$message
[1] "CONVERGENCE: NORM OF PROJECTED GRADIENT <= PGTOL"library(NlcOptim)
solnl(X=x0,objfun = obj,lb=c(-1,-1,-1),ub=c(1,1,1))$par
[,1]
[1,] 0.1930158
[2,] 1.0000000
[3,] 1.0000000
$fn
1
-5.715089
$counts
nfval ngval
[1,] 21 6
$lambda
$lambda$lower
[,1]
[1,] 0
[2,] 0
[3,] 0
$lambda$upper
[,1]
[1,] 0.000000
[2,] 0.983750
[3,] 1.611776
$grad
[,1]
[1,] 2.145767e-06
[2,] -9.837500e-01
[3,] -1.611776e+00
$hessian
[,1] [,2] [,3]
[1,] 4.144752e-01 -7.325474e-07 -7.325635e-07
[2,] -7.325474e-07 5.011736e-03 -1.321530e-02
[3,] -7.325635e-07 -1.321530e-02 3.485172e-02DOE run.order std.order ang compr alt y
1 1 1 2 3 3 2.71
2 2 2 6 3 3 2.75
3 3 3 2 5 3 4.17
4 4 4 6 5 3 4.12
5 5 5 2 4 2 2.42
6 6 6 6 4 2 2.61
7 7 7 2 4 4 4.36
8 8 8 6 4 4 4.82
9 9 9 4 3 2 2.65
10 10 10 4 5 2 3.43
11 11 11 4 3 4 3.68
12 12 12 4 5 4 5.62
13 13 13 4 4 3 3.43
14 14 14 4 4 3 3.45
15 15 15 4 4 3 3.82
16 16 16 4 4 3 3.46
17 17 17 4 4 3 3.68
Data are stored in coded form using these coding formulas ...
x1 ~ (ang - 4)/2
x2 ~ (compr - 4)/1
x3 ~ (alt - 3)/1#Gráficos
persp(funcao,~x1+x2,zlab = "y",col="blue",
contours = "colors",theta = 30,phi = 30)
persp(funcao,~x1+x3,zlab = "y",col="blue",
contours = "colors",theta = 30,phi = 30)persp(funcao,~x2+x3,zlab = "y",col="blue",
contours = "colors",theta = 30,phi = 30)
contour(funcao,~x1+x2,image = T)points(3.6140,5,col="blue",lwd=3)
contour(funcao,~x1+x3,col="black")points(3.6140,4,col="blue",lwd=3)
contour(funcao,~x2+x3,col=rainbow(10))points(5,4,col="blue",lwd=3)
predict(funcao,newdata = data.frame(x1=0,x2=1,x3=1),
se.fit = T,interval = "conf")$fit
fit lwr upr
1 5.707368 5.37123 6.043507
$se.fit
[1] 0.1508606
$df
[1] 10
$residual.scale
[1] 0.1823815#Validação
DOE run.order std.order ang compr alt y
1 1 1 2 3 3 2.71
2 2 2 6 3 3 2.75
3 3 3 2 5 3 4.17
4 4 4 6 5 3 4.12
5 5 5 2 4 2 2.42
6 6 6 6 4 2 2.61
7 7 7 2 4 4 4.36
8 8 8 6 4 4 4.82
9 9 9 4 3 2 2.65
10 10 10 4 5 2 3.43
11 11 11 4 3 4 3.68
12 12 12 4 5 4 5.62
13 13 13 4 4 3 3.43
14 14 14 4 4 3 3.45
15 15 15 4 4 3 3.82
16 16 16 4 4 3 3.46
17 17 17 4 4 3 3.68
Data are stored in coded form using these coding formulas ...
x1 ~ (ang - 4)/2
x2 ~ (compr - 4)/1
x3 ~ (alt - 3)/1predict(funcao,newdata = data.frame(x1=0,x2=1,x3=1),
se.fit = T,interval = "conf")$fit
fit lwr upr
1 5.707368 5.37123 6.043507
$se.fit
[1] 0.1508606
$df
[1] 10
$residual.scale
[1] 0.1823815predict(funcao,newdata = data.frame(x1=0,x2=1,x3=1),
se.fit = T,interval = "pred")$fit
fit lwr upr
1 5.707368 5.179991 6.234746
$se.fit
[1] 0.1508606
$df
[1] 10
$residual.scale
[1] 0.1823815LS0tDQp0aXRsZTogIkV4cGVyaW1lbnRvIENhdGFwdWx0YTogUlNNIEJPWC1CRUhOS0VOIERFU0lHTjogTW9kZWxhZ2VtIGUgT3RpbWl6YcOnw6NvIHBlbGEgTGluZ3VhZ2VtIFIiDQpvdXRwdXQ6IGh0bWxfbm90ZWJvb2sNCi0tLQ0KIA0KDQogDQoNCmBgYHtyfQ0KbGlicmFyeShyc20pDQpsaWJyYXJ5KGNhcikNCmxpYnJhcnkobG10ZXN0KQ0KbGlicmFyeShub3J0ZXN0KQ0KbGlicmFyeShSc29sbnApDQpsaWJyYXJ5KE5sY09wdGltKQ0KRE9FPC1iYmQoaz1+eDEreDIreDMsDQogICAgICAgICBibG9jayA9IEYsDQogICAgICAgICBuMD01LA0KICAgICAgICAgcmFuZG9taXplID0gRiwNCiAgICAgICAgIGNvZGluZyA9IGxpc3QoeDF+KGFuZy00KS8yLA0KICAgICAgICAgICAgICAgICAgICAgICB4Mn4oY29tcHItNCkvMSwNCiAgICAgICAgICAgICAgICAgICAgICAgeDN+KGFsdC0zKS8xKSkNCkRPRQ0KDQpET0UkeTwtKGMoMi43MSwyLjc1LDQuMTcsNC4xMiwyLjQyLDIuNjEsNC4zNiwNCiAgICAgICAgICAgNC44MiwyLjY1LDMuNDMsMy42OCw1LjYyLDMuNDMsMy40NSwNCiAgICAgICAgICAgMy44MiwzLjQ2LDMuNjgpKQ0KRE9FDQoNCmZ1bmNhbzwtIHJzbSh5fiBTTyh4MSx4Mix4MyksZGF0YSA9IERPRSkNCnN1bW1hcnkoZnVuY2FvKQ0KZnVuY2FvPC1yc20oeX4gRk8oeDEseDIseDMpKw0KICAgICAgICAgICAgIFRXSSh4Mix4MykrUFEoeDEseDMpLGRhdGEgPSBET0UpDQpzdW1tYXJ5KGZ1bmNhbykNCmBgYA0KDQpgYGB7cn0NCiM9PT09PT09PSBBbsOhbGlzZSBkb3MgUmVzw61kdW9zID09PT09PT09PT09PT09PQ0KY3VydmUoZG5vcm0oeCxtZWFuKGZ1bmNhbyRyZXNpZHVhbHMpLHNkKGZ1bmNhbyRyZXNpZHVhbHMpKSwNCiAgICAgIGZyb20gPSBtaW4oZnVuY2FvJHJlc2lkdWFscyksdG89bWF4KGZ1bmNhbyRyZXNpZHVhbHMpLA0KICAgICAgY29sPSJibHVlIixsd2Q9NCxtYWluPSJDdXJ2YSBOb3JtYWwgZG9zIFJlc8OtZHVvcyIsDQogICAgICB4bGFiID0gIk9ic2VydmHDp8O1ZXMiLHlsYWIgPSAiRGVuc2lkYWRlIikNCnFxbm9ybShmdW5jYW8kcmVzaWR1YWxzLGx3ZD0zLGNvbD0iYmx1ZSIpDQpxcWxpbmUoZnVuY2FvJHJlc2lkdWFscyxsd2Q9MykNCmxpYnJhcnkobm9ydGVzdCkNCmFkLnRlc3QoZnVuY2FvJHJlc2lkdWFscykNCnJlc2lkdWFsUGxvdChmdW5jYW8pDQpicHRlc3QoZnVuY2FvKQ0KZHVyYmluV2F0c29uVGVzdChmdW5jYW8pDQpkd3QoZnVuY2FvKQ0KDQpgYGANCg0KDQpgYGB7cn0NCiNPdGltaXphw6fDo28NCm9iaj1mdW5jdGlvbih4KXthPC1wcmVkaWN0KGZ1bmNhbyxuZXdkYXRhID0gZGF0YS5mcmFtZSh4MT14WzFdLA0KICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIHgyPXhbMl0sDQogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgeDM9eFszXSkpDQpyZXR1cm4oLWEpfQ0KDQp4MDwtYygwLDAsMCkNCm9wdGltKHBhciA9IHgwLGZuPW9iaixsb3dlciA9IGMoLTEsLTEsLTEpLHVwcGVyID0gYygxLDEsMSksDQogICAgICBtZXRob2QgPSBjKCJMLUJGR1MtQiIpKQ0KDQpsaWJyYXJ5KE5sY09wdGltKQ0Kc29sbmwoWD14MCxvYmpmdW4gPSBvYmosbGI9YygtMSwtMSwtMSksdWI9YygxLDEsMSkpDQpET0UNCg0KI0dyw6FmaWNvcw0KcGVyc3AoZnVuY2FvLH54MSt4Mix6bGFiID0gInkiLGNvbD0iYmx1ZSIsDQogICAgICBjb250b3VycyA9ICJjb2xvcnMiLHRoZXRhID0gMzAscGhpID0gMzApDQpwZXJzcChmdW5jYW8sfngxK3gzLHpsYWIgPSAieSIsY29sPSJibHVlIiwNCiAgICAgIGNvbnRvdXJzID0gImNvbG9ycyIsdGhldGEgPSAzMCxwaGkgPSAzMCkNCnBlcnNwKGZ1bmNhbyx+eDIreDMsemxhYiA9ICJ5Iixjb2w9ImJsdWUiLA0KICAgICAgY29udG91cnMgPSAiY29sb3JzIix0aGV0YSA9IDMwLHBoaSA9IDMwKQ0KY29udG91cihmdW5jYW8sfngxK3gyLGltYWdlID0gVCkNCnBvaW50cygzLjYxNDAsNSxjb2w9ImJsdWUiLGx3ZD0zKQ0KDQpjb250b3VyKGZ1bmNhbyx+eDEreDMsY29sPSJibGFjayIpDQpwb2ludHMoMy42MTQwLDQsY29sPSJibHVlIixsd2Q9MykNCg0KY29udG91cihmdW5jYW8sfngyK3gzLGNvbD1yYWluYm93KDEwKSkNCnBvaW50cyg1LDQsY29sPSJibHVlIixsd2Q9MykNCg0KcHJlZGljdChmdW5jYW8sbmV3ZGF0YSA9IGRhdGEuZnJhbWUoeDE9MCx4Mj0xLHgzPTEpLA0KICAgICAgICBzZS5maXQgPSBULGludGVydmFsID0gImNvbmYiKQ0KDQpgYGANCg0KDQpgYGB7cn0NCiNWYWxpZGHDp8Ojbw0KRE9FDQpwcmVkaWN0KGZ1bmNhbyxuZXdkYXRhID0gZGF0YS5mcmFtZSh4MT0wLHgyPTEseDM9MSksDQogICAgICAgIHNlLmZpdCA9IFQsaW50ZXJ2YWwgPSAiY29uZiIpDQpwcmVkaWN0KGZ1bmNhbyxuZXdkYXRhID0gZGF0YS5mcmFtZSh4MT0wLHgyPTEseDM9MSksDQogICAgICAgIHNlLmZpdCA9IFQsaW50ZXJ2YWwgPSAicHJlZCIpDQpgYGANCg0KDQo=