RSM-CCD-ARRANJO DE FACE: EXPERIMENTO DA CATAPULTA
Experimento da catapulta RSM-CCD-FACE. Experimento real da catapulta realizado pelo primeiro ano de Engenharia de Produção.
library(rsm)
library(car)
library(lmtest)
library(rio)
library(FrF2)
library(nortest)
library(NlcOptim)
library(broom)
DOE<-ccd(basis = ~x1+x2+x3,
oneblock = T,
alpha = "face",
randomize = F,
n0=c(3,3),
coding = list(x1~(ang-4)/2,
x2~(compr-3)/1,
x3~(alt-3)/1))
DOE run.order std.order ang compr alt
1 1 1 2 2 2
2 2 2 6 2 2
3 3 3 2 4 2
4 4 4 6 4 2
5 5 5 2 2 4
6 6 6 6 2 4
7 7 7 2 4 4
8 8 8 6 4 4
9 9 9 4 3 3
10 10 10 4 3 3
11 11 11 4 3 3
12 1 1 2 3 3
13 2 2 6 3 3
14 3 3 4 2 3
15 4 4 4 4 3
16 5 5 4 3 2
17 6 6 4 3 4
18 7 7 4 3 3
19 8 8 4 3 3
20 9 9 4 3 3
Data are stored in coded form using these coding formulas ...
x1 ~ (ang - 4)/2
x2 ~ (compr - 3)/1
x3 ~ (alt - 3)/1DOE<-djoin(DOE,DOE) # Blocagem
DOE run.order std.order ang compr alt Block
1 1 1 2 2 2 1
2 2 2 6 2 2 1
3 3 3 2 4 2 1
4 4 4 6 4 2 1
5 5 5 2 2 4 1
6 6 6 6 2 4 1
7 7 7 2 4 4 1
8 8 8 6 4 4 1
9 9 9 4 3 3 1
10 10 10 4 3 3 1
11 11 11 4 3 3 1
12 1 1 2 3 3 1
13 2 2 6 3 3 1
14 3 3 4 2 3 1
15 4 4 4 4 3 1
16 5 5 4 3 2 1
17 6 6 4 3 4 1
18 7 7 4 3 3 1
19 8 8 4 3 3 1
20 9 9 4 3 3 1
21 1 1 2 2 2 2
22 2 2 6 2 2 2
23 3 3 2 4 2 2
24 4 4 6 4 2 2
25 5 5 2 2 4 2
26 6 6 6 2 4 2
27 7 7 2 4 4 2
28 8 8 6 4 4 2
29 9 9 4 3 3 2
30 10 10 4 3 3 2
31 11 11 4 3 3 2
32 1 1 2 3 3 2
33 2 2 6 3 3 2
34 3 3 4 2 3 2
35 4 4 4 4 3 2
36 5 5 4 3 2 2
37 6 6 4 3 4 2
38 7 7 4 3 3 2
39 8 8 4 3 3 2
40 9 9 4 3 3 2
Data are stored in coded form using these coding formulas ...
x1 ~ (ang - 4)/2
x2 ~ (compr - 3)/1
x3 ~ (alt - 3)/1#set.seed(30);DOE<-DOE[sample(nrow(DOE)),]
#DOE
DOE$y<-c(0.10,0.10,2.40,2.35,3.36,3.10,5.66,
5.84, 2.44, 2.51, 2.58, 2.71, 2.50, 1.25,3.88,
0.90, 4.12, 2.55,2.35,2.40,
0.12,0.12,2.40,2.45,3.31,3.10,5.68,
5.80,2.58, 2.67, 2.59,2.77,2.62,1.38,3.95,
1.05,4.34,2.42,2.70,2.72)
DOE run.order std.order ang compr alt Block y
1 1 1 2 2 2 1 0.10
2 2 2 6 2 2 1 0.10
3 3 3 2 4 2 1 2.40
4 4 4 6 4 2 1 2.35
5 5 5 2 2 4 1 3.36
6 6 6 6 2 4 1 3.10
7 7 7 2 4 4 1 5.66
8 8 8 6 4 4 1 5.84
9 9 9 4 3 3 1 2.44
10 10 10 4 3 3 1 2.51
11 11 11 4 3 3 1 2.58
12 1 1 2 3 3 1 2.71
13 2 2 6 3 3 1 2.50
14 3 3 4 2 3 1 1.25
15 4 4 4 4 3 1 3.88
16 5 5 4 3 2 1 0.90
17 6 6 4 3 4 1 4.12
18 7 7 4 3 3 1 2.55
19 8 8 4 3 3 1 2.35
20 9 9 4 3 3 1 2.40
21 1 1 2 2 2 2 0.12
22 2 2 6 2 2 2 0.12
23 3 3 2 4 2 2 2.40
24 4 4 6 4 2 2 2.45
25 5 5 2 2 4 2 3.31
26 6 6 6 2 4 2 3.10
27 7 7 2 4 4 2 5.68
28 8 8 6 4 4 2 5.80
29 9 9 4 3 3 2 2.58
30 10 10 4 3 3 2 2.67
31 11 11 4 3 3 2 2.59
32 1 1 2 3 3 2 2.77
33 2 2 6 3 3 2 2.62
34 3 3 4 2 3 2 1.38
35 4 4 4 4 3 2 3.95
36 5 5 4 3 2 2 1.05
37 6 6 4 3 4 2 4.34
38 7 7 4 3 3 2 2.42
39 8 8 4 3 3 2 2.70
40 9 9 4 3 3 2 2.72
Data are stored in coded form using these coding formulas ...
x1 ~ (ang - 4)/2
x2 ~ (compr - 3)/1
x3 ~ (alt - 3)/1funcao<-rsm(y~SO(x1,x2,x3)+Block,data = DOE)
summary(funcao)
Call:
rsm(formula = y ~ SO(x1, x2, x3) + Block, data = DOE)
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.488295 0.026840 92.7075 < 2.2e-16 ***
x1 -0.026500 0.020697 -1.2804 0.210549
x2 1.223500 0.020697 59.1162 < 2.2e-16 ***
x3 1.616000 0.020697 78.0807 < 2.2e-16 ***
x1:x2 0.048125 0.023139 2.0798 0.046492 *
x1:x3 -0.010625 0.023139 -0.4592 0.649533
x2:x3 0.059375 0.023139 2.5660 0.015720 *
x1^2 0.138636 0.039467 3.5127 0.001474 **
x2^2 0.103636 0.039467 2.6259 0.013654 *
x3^2 0.091136 0.039467 2.3092 0.028254 *
Block2 0.083500 0.029269 2.8528 0.007912 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Multiple R-squared: 0.997, Adjusted R-squared: 0.996
F-statistic: 970.8 on 10 and 29 DF, p-value: < 2.2e-16
Analysis of Variance Table
Response: y
Df Sum Sq Mean Sq F value Pr(>F)
FO(x1, x2, x3) 3 82.182 27.3941 3197.6524 < 2.2e-16
TWI(x1, x2, x3) 3 0.095 0.0318 3.7068 0.022611
PQ(x1, x2, x3) 3 0.820 0.2733 31.9070 2.538e-09
Block 1 0.070 0.0697 8.1386 0.007912
Residuals 29 0.248 0.0086
Lack of fit 5 0.094 0.0187 2.8985 0.034755
Pure error 24 0.155 0.0065
Stationary point of response surface:
x1 x2 x3
0.4710516 -3.8385847 -7.5879652
Stationary point in original units:
ang compr alt
4.9421032 -0.8385847 -4.5879652
Eigenanalysis:
eigen() decomposition
$values
[1] 0.15248949 0.11896004 0.06195956
$vectors
[,1] [,2] [,3]
x1 0.8349161 -0.4879181 0.2546589
x2 0.5202954 0.5488353 -0.6542726
x3 0.1794656 0.6787605 0.7120929funcao<-rsm(y~ FO(x1,x2,x3)+PQ(x1,x2,x3)+
TWI(x1,x2)+TWI(x2,x3)+Block,data = DOE)
summary(funcao)
Call:
rsm(formula = y ~ FO(x1, x2, x3) + PQ(x1, x2, x3) + TWI(x1, x2) +
TWI(x2, x3) + Block, data = DOE)
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.488295 0.026485 93.9514 < 2.2e-16 ***
x1 -0.026500 0.020423 -1.2976 0.204318
x2 1.223500 0.020423 59.9094 < 2.2e-16 ***
x3 1.616000 0.020423 79.1284 < 2.2e-16 ***
x1^2 0.138636 0.038944 3.5599 0.001259 **
x2^2 0.103636 0.038944 2.6612 0.012387 *
x3^2 0.091136 0.038944 2.3402 0.026118 *
x1:x2 0.048125 0.022833 2.1077 0.043522 *
x2:x3 0.059375 0.022833 2.6004 0.014314 *
Block2 0.083500 0.028882 2.8911 0.007074 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Multiple R-squared: 0.997, Adjusted R-squared: 0.9961
F-statistic: 1108 on 9 and 30 DF, p-value: < 2.2e-16
Analysis of Variance Table
Response: y
Df Sum Sq Mean Sq F value Pr(>F)
FO(x1, x2, x3) 3 82.182 27.3941 3284.0402 < 2.2e-16
PQ(x1, x2, x3) 3 0.820 0.2733 32.7690 1.35e-09
TWI(x1, x2) 1 0.037 0.0371 4.4424 0.043522
TWI(x2, x3) 1 0.056 0.0564 6.7621 0.014314
Block 1 0.070 0.0697 8.3584 0.007074
Residuals 30 0.250 0.0083
Lack of fit 6 0.095 0.0159 2.4621 0.053459
Pure error 24 0.155 0.0065
Stationary point of response surface:
x1 x2 x3
0.7737689 -3.9074286 -7.5929979
Stationary point in original units:
ang compr alt
5.5475377 -0.9074286 -4.5929979
Eigenanalysis:
eigen() decomposition
$values
[1] 0.15437132 0.11529979 0.06373797
$vectors
[,1] [,2] [,3]
x1 0.8106107 0.5455157 0.2128918
x2 0.5300749 -0.5290582 -0.6626598
x3 0.2488591 -0.6500077 0.7180244curve(dnorm(x,mean(funcao$residuals),sd(funcao$residuals)),
from = min(funcao$residuals),to=max(funcao$residuals),
col="purple",lwd=3,xlab = "Observações",ylab = "Densidade",
main="Análise dos ResÃduos")
qqnorm(funcao$residuals,lwd=2)
qqline(funcao$residuals,col="red",lwd=2)
shapiro.test(funcao$residuals)
Shapiro-Wilk normality test
data: funcao$residuals
W = 0.97466, p-value = 0.4985ad.test(funcao$residuals)
Anderson-Darling normality test
data: funcao$residuals
A = 0.28226, p-value = 0.6189residualPlot(funcao)
bptest(funcao) # Homocedasticidade
studentized Breusch-Pagan test
data: funcao
BP = 5.7539, df = 9, p-value = 0.7643durbinWatsonTest(funcao) lag Autocorrelation D-W Statistic p-value
1 -0.01287799 1.936935 0.466
Alternative hypothesis: rho != 0dwt(funcao) lag Autocorrelation D-W Statistic p-value
1 -0.01287799 1.936935 0.512
Alternative hypothesis: rho != 0# Otimizar o modelo
obj=function(x){a<-predict(funcao,newdata = data.frame(x1=x[1],
x2=x[2],
x3=x[3],
Block="2"))
return(-a)}
x0<-c(0,0,0)
optim(par = x0,fn=obj,lower = c(-1,-1,-1),upper = c(1,1,1),
method = "L-BFGS-B")$par
[1] 1 1 1
$value
[1] -5.825705
$counts
function gradient
6 6
$convergence
[1] 0
$message
[1] "CONVERGENCE: NORM OF PROJECTED GRADIENT <= PGTOL"persp(funcao,~x1+x2,zlab = "Dist",col="blue",
contours = "colors")
persp(funcao,~x1+x3,zlab = "Dist",col="blue",
contours = "colors")
persp(funcao,~x2+x3,zlab = "Dist",col="blue",
contours = "colors")
contour(funcao,~x1+x2,col=topo.colors(10))
points(6,4,col="red",lwd=4)
contour(funcao,~x1+x3,col=topo.colors(5))
points(6,4,col="red",lwd=4)
contour(funcao,~x2+x3,col="blue")
points(4,4,col="red",lwd=4)
DOE run.order std.order ang compr alt Block y
1 1 1 2 2 2 1 0.10
2 2 2 6 2 2 1 0.10
3 3 3 2 4 2 1 2.40
4 4 4 6 4 2 1 2.35
5 5 5 2 2 4 1 3.36
6 6 6 6 2 4 1 3.10
7 7 7 2 4 4 1 5.66
8 8 8 6 4 4 1 5.84
9 9 9 4 3 3 1 2.44
10 10 10 4 3 3 1 2.51
11 11 11 4 3 3 1 2.58
12 1 1 2 3 3 1 2.71
13 2 2 6 3 3 1 2.50
14 3 3 4 2 3 1 1.25
15 4 4 4 4 3 1 3.88
16 5 5 4 3 2 1 0.90
17 6 6 4 3 4 1 4.12
18 7 7 4 3 3 1 2.55
19 8 8 4 3 3 1 2.35
20 9 9 4 3 3 1 2.40
21 1 1 2 2 2 2 0.12
22 2 2 6 2 2 2 0.12
23 3 3 2 4 2 2 2.40
24 4 4 6 4 2 2 2.45
25 5 5 2 2 4 2 3.31
26 6 6 6 2 4 2 3.10
27 7 7 2 4 4 2 5.68
28 8 8 6 4 4 2 5.80
29 9 9 4 3 3 2 2.58
30 10 10 4 3 3 2 2.67
31 11 11 4 3 3 2 2.59
32 1 1 2 3 3 2 2.77
33 2 2 6 3 3 2 2.62
34 3 3 4 2 3 2 1.38
35 4 4 4 4 3 2 3.95
36 5 5 4 3 2 2 1.05
37 6 6 4 3 4 2 4.34
38 7 7 4 3 3 2 2.42
39 8 8 4 3 3 2 2.70
40 9 9 4 3 3 2 2.72
Data are stored in coded form using these coding formulas ...
x1 ~ (ang - 4)/2
x2 ~ (compr - 3)/1
x3 ~ (alt - 3)/1predict(funcao,newdata = data.frame(x1=1,x2=1,x3=1,Block="2"),
se.fit = T,interval = "pred")$fit
fit lwr upr
1 5.825705 5.608253 6.043156
$se.fit
[1] 0.0547301
$df
[1] 30
$residual.scale
[1] 0.09133223predict(funcao,newdata = data.frame(x1=1,x2=1,x3=1,Block="1"),
se.fit = T,interval = "pred")$fit
fit lwr upr
1 5.742205 5.524753 5.959656
$se.fit
[1] 0.0547301
$df
[1] 30
$residual.scale
[1] 0.09133223library(NlcOptim)
solnl(X=x0,objfun = obj, lb=c(-1,-1,-1),ub=c(1,1,1))$par
[,1]
[1,] 1
[2,] 1
[3,] 1
$fn
1
-5.825705
$counts
nfval ngval
[1,] 31 6
$lambda
$lambda$lower
[,1]
[1,] 0
[2,] 0
[3,] 0
$lambda$upper
[,1]
[1,] 0.2988977
[2,] 1.5382727
[3,] 1.8576477
$grad
[,1]
[1,] -0.2988977
[2,] -1.5382727
[3,] -1.8576477
$hessian
[,1] [,2] [,3]
[1,] 1.993231e-02 -0.04812498 1.759114e-15
[2,] -4.812498e-02 0.61619424 -4.999991e-01
[3,] 1.759114e-15 -0.49999914 5.000027e-01LS0tDQp0aXRsZTogIlJTTS1DQ0QtQVJSQU5KTyBERSBGQUNFOiBFWFBFUklNRU5UTyBEQSBDQVRBUFVMVEEiDQpvdXRwdXQ6IGh0bWxfbm90ZWJvb2sNCi0tLQ0KDQpFeHBlcmltZW50byBkYSBjYXRhcHVsdGEgW1JTTS1DQ0QtRkFDRV0oaHR0cDovL3JtYXJrZG93bi5yc3R1ZGlvLmNvbSkuIEV4cGVyaW1lbnRvIHJlYWwgZGEgY2F0YXB1bHRhIHJlYWxpemFkbyBwZWxvIHByaW1laXJvIGFubyBkZSBFbmdlbmhhcmlhIGRlIFByb2R1w6fDo28uDQpgYGB7cn0NCmxpYnJhcnkocnNtKQ0KbGlicmFyeShjYXIpDQpsaWJyYXJ5KGxtdGVzdCkNCmxpYnJhcnkocmlvKQ0KbGlicmFyeShGckYyKQ0KbGlicmFyeShub3J0ZXN0KQ0KbGlicmFyeShObGNPcHRpbSkNCmxpYnJhcnkoYnJvb20pDQoNCkRPRTwtY2NkKGJhc2lzID0gfngxK3gyK3gzLA0KICAgICAgICAgb25lYmxvY2sgPSBULA0KICAgICAgICAgYWxwaGEgPSAiZmFjZSIsDQogICAgICAgICByYW5kb21pemUgPSBGLA0KICAgICAgICAgbjA9YygzLDMpLA0KICAgICAgICAgY29kaW5nID0gbGlzdCh4MX4oYW5nLTQpLzIsDQogICAgICAgICAgICAgICAgICAgICAgIHgyfihjb21wci0zKS8xLA0KICAgICAgICAgICAgICAgICAgICAgICB4M34oYWx0LTMpLzEpKQ0KRE9FDQoNCkRPRTwtZGpvaW4oRE9FLERPRSkgIyBCbG9jYWdlbQ0KRE9FDQoNCiNzZXQuc2VlZCgzMCk7RE9FPC1ET0Vbc2FtcGxlKG5yb3coRE9FKSksXQ0KI0RPRQ0KDQpET0UkeTwtYygwLjEwLDAuMTAsMi40MCwyLjM1LDMuMzYsMy4xMCw1LjY2LCANCjUuODQsIDIuNDQsIDIuNTEsIDIuNTgsIDIuNzEsIDIuNTAsIDEuMjUsMy44OCwNCjAuOTAsIDQuMTIsIDIuNTUsMi4zNSwyLjQwLA0KICANCjAuMTIsMC4xMiwyLjQwLDIuNDUsMy4zMSwzLjEwLDUuNjgsDQo1LjgwLDIuNTgsIDIuNjcsIDIuNTksMi43NywyLjYyLDEuMzgsMy45NSwNCjEuMDUsNC4zNCwyLjQyLDIuNzAsMi43MikNCg0KRE9FDQpgYGANCg0KYGBge3J9DQpmdW5jYW88LXJzbSh5flNPKHgxLHgyLHgzKStCbG9jayxkYXRhID0gRE9FKQ0Kc3VtbWFyeShmdW5jYW8pDQoNCmZ1bmNhbzwtcnNtKHl+IEZPKHgxLHgyLHgzKStQUSh4MSx4Mix4MykrDQogICAgICAgICAgICAgIFRXSSh4MSx4MikrVFdJKHgyLHgzKStCbG9jayxkYXRhID0gRE9FKQ0Kc3VtbWFyeShmdW5jYW8pDQoNCmN1cnZlKGRub3JtKHgsbWVhbihmdW5jYW8kcmVzaWR1YWxzKSxzZChmdW5jYW8kcmVzaWR1YWxzKSksDQogICAgICBmcm9tID0gbWluKGZ1bmNhbyRyZXNpZHVhbHMpLHRvPW1heChmdW5jYW8kcmVzaWR1YWxzKSwNCiAgICAgIGNvbD0icHVycGxlIixsd2Q9Myx4bGFiID0gIk9ic2VydmHDp8O1ZXMiLHlsYWIgPSAiRGVuc2lkYWRlIiwNCiAgICAgIG1haW49IkFuw6FsaXNlIGRvcyBSZXPDrWR1b3MiKQ0KcXFub3JtKGZ1bmNhbyRyZXNpZHVhbHMsbHdkPTIpDQpxcWxpbmUoZnVuY2FvJHJlc2lkdWFscyxjb2w9InJlZCIsbHdkPTIpDQpzaGFwaXJvLnRlc3QoZnVuY2FvJHJlc2lkdWFscykNCmFkLnRlc3QoZnVuY2FvJHJlc2lkdWFscykNCg0KcmVzaWR1YWxQbG90KGZ1bmNhbykNCmJwdGVzdChmdW5jYW8pICMgSG9tb2NlZGFzdGljaWRhZGUNCmR1cmJpbldhdHNvblRlc3QoZnVuY2FvKQ0KZHd0KGZ1bmNhbykNCg0KYGBgDQoNCmBgYHtyfQ0KIyBPdGltaXphciBvIG1vZGVsbw0Kb2JqPWZ1bmN0aW9uKHgpe2E8LXByZWRpY3QoZnVuY2FvLG5ld2RhdGEgPSBkYXRhLmZyYW1lKHgxPXhbMV0sDQogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgeDI9eFsyXSwNCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICB4Mz14WzNdLA0KICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIEJsb2NrPSIyIikpDQpyZXR1cm4oLWEpfQ0KDQp4MDwtYygwLDAsMCkNCm9wdGltKHBhciA9IHgwLGZuPW9iaixsb3dlciA9IGMoLTEsLTEsLTEpLHVwcGVyID0gYygxLDEsMSksDQogICAgICBtZXRob2QgPSAiTC1CRkdTLUIiKQ0KDQoNCnBlcnNwKGZ1bmNhbyx+eDEreDIsemxhYiA9ICJEaXN0Iixjb2w9ImJsdWUiLA0KICAgICAgY29udG91cnMgPSAiY29sb3JzIikNCnBlcnNwKGZ1bmNhbyx+eDEreDMsemxhYiA9ICJEaXN0Iixjb2w9ImJsdWUiLA0KICAgICAgY29udG91cnMgPSAiY29sb3JzIikNCg0KcGVyc3AoZnVuY2FvLH54Mit4Myx6bGFiID0gIkRpc3QiLGNvbD0iYmx1ZSIsDQogICAgICBjb250b3VycyA9ICJjb2xvcnMiKQ0KY29udG91cihmdW5jYW8sfngxK3gyLGNvbD10b3BvLmNvbG9ycygxMCkpDQpwb2ludHMoNiw0LGNvbD0icmVkIixsd2Q9NCkNCmNvbnRvdXIoZnVuY2FvLH54MSt4Myxjb2w9dG9wby5jb2xvcnMoNSkpDQpwb2ludHMoNiw0LGNvbD0icmVkIixsd2Q9NCkNCmNvbnRvdXIoZnVuY2FvLH54Mit4Myxjb2w9ImJsdWUiKQ0KcG9pbnRzKDQsNCxjb2w9InJlZCIsbHdkPTQpDQoNCkRPRQ0KDQpwcmVkaWN0KGZ1bmNhbyxuZXdkYXRhID0gZGF0YS5mcmFtZSh4MT0xLHgyPTEseDM9MSxCbG9jaz0iMiIpLA0KICAgICAgICBzZS5maXQgPSBULGludGVydmFsID0gInByZWQiKQ0KDQpwcmVkaWN0KGZ1bmNhbyxuZXdkYXRhID0gZGF0YS5mcmFtZSh4MT0xLHgyPTEseDM9MSxCbG9jaz0iMSIpLA0KICAgICAgICBzZS5maXQgPSBULGludGVydmFsID0gInByZWQiKQ0KDQpsaWJyYXJ5KE5sY09wdGltKQ0Kc29sbmwoWD14MCxvYmpmdW4gPSBvYmosIGxiPWMoLTEsLTEsLTEpLHViPWMoMSwxLDEpKQ0KDQpgYGANCg0K