Chapter 50 Response Surfaces Designs
50.1 Choosing Designs
The rsm
packages has several functions that help you
pick/generate central composite designs.
Generally the assumption is that the design will
be run in two blocks: (1) factorial points and some center points,
and (2) axial points and some center points.
The ccd()
function generates a CCD design. By default,
it creates two blocks and chooses alpha for
orthogonal blocking. Here we say not
to randomize the design so that we can
see the design more clearly.
> library(cfcdae)
Registered S3 method overwritten by 'DoE.base':
method from
factorize.factor conf.design
> library(rsm)
> ccd(3,randomize=FALSE)
run.order std.order x1.as.is x2.as.is x3.as.is Block
1 1 1 -1.000000 -1.000000 -1.000000 1
2 2 2 1.000000 -1.000000 -1.000000 1
3 3 3 -1.000000 1.000000 -1.000000 1
4 4 4 1.000000 1.000000 -1.000000 1
5 5 5 -1.000000 -1.000000 1.000000 1
6 6 6 1.000000 -1.000000 1.000000 1
7 7 7 -1.000000 1.000000 1.000000 1
8 8 8 1.000000 1.000000 1.000000 1
9 9 9 0.000000 0.000000 0.000000 1
10 10 10 0.000000 0.000000 0.000000 1
11 11 11 0.000000 0.000000 0.000000 1
12 12 12 0.000000 0.000000 0.000000 1
13 1 1 -1.825742 0.000000 0.000000 2
14 2 2 1.825742 0.000000 0.000000 2
15 3 3 0.000000 -1.825742 0.000000 2
16 4 4 0.000000 1.825742 0.000000 2
17 5 5 0.000000 0.000000 -1.825742 2
18 6 6 0.000000 0.000000 1.825742 2
19 7 7 0.000000 0.000000 0.000000 2
20 8 8 0.000000 0.000000 0.000000 2
21 9 9 0.000000 0.000000 0.000000 2
22 10 10 0.000000 0.000000 0.000000 2
Data are stored in coded form using these coding formulas ...
x1 ~ x1.as.is
x2 ~ x2.as.is
x3 ~ x3.as.is
You could tell it one block with a rotatable alpha.
> ccd(3,alpha="rotatable",oneblock=TRUE,randomize=FALSE)
run.order std.order x1.as.is x2.as.is x3.as.is
1 1 1 -1.000000 -1.000000 -1.000000
2 2 2 1.000000 -1.000000 -1.000000
3 3 3 -1.000000 1.000000 -1.000000
4 4 4 1.000000 1.000000 -1.000000
5 5 5 -1.000000 -1.000000 1.000000
6 6 6 1.000000 -1.000000 1.000000
7 7 7 -1.000000 1.000000 1.000000
8 8 8 1.000000 1.000000 1.000000
9 9 9 0.000000 0.000000 0.000000
10 10 10 0.000000 0.000000 0.000000
11 11 11 0.000000 0.000000 0.000000
12 12 12 0.000000 0.000000 0.000000
13 1 1 -1.681793 0.000000 0.000000
14 2 2 1.681793 0.000000 0.000000
15 3 3 0.000000 -1.681793 0.000000
16 4 4 0.000000 1.681793 0.000000
17 5 5 0.000000 0.000000 -1.681793
18 6 6 0.000000 0.000000 1.681793
19 7 7 0.000000 0.000000 0.000000
20 8 8 0.000000 0.000000 0.000000
21 9 9 0.000000 0.000000 0.000000
22 10 10 0.000000 0.000000 0.000000
Data are stored in coded form using these coding formulas ...
x1 ~ x1.as.is
x2 ~ x2.as.is
x3 ~ x3.as.is
You can specify a fractional factorial for the factorial part, in this case with the factorial part broken into two blocks.
> ccd(~A+B+C+D, generators=E~A*B*C*D, block = Day~A*B*C,randomize=FALSE)
run.order std.order A.as.is B.as.is C.as.is D.as.is E.as.is Day
1 1 1 -1.000000 -1.000000 -1.000000 -1.000000 1.000000 1
4 2 2 1.000000 1.000000 -1.000000 -1.000000 1.000000 1
6 3 3 1.000000 -1.000000 1.000000 -1.000000 1.000000 1
7 4 4 -1.000000 1.000000 1.000000 -1.000000 1.000000 1
9 5 5 -1.000000 -1.000000 -1.000000 1.000000 -1.000000 1
12 6 6 1.000000 1.000000 -1.000000 1.000000 -1.000000 1
14 7 7 1.000000 -1.000000 1.000000 1.000000 -1.000000 1
15 8 8 -1.000000 1.000000 1.000000 1.000000 -1.000000 1
11 9 9 0.000000 0.000000 0.000000 0.000000 0.000000 1
2 10 10 0.000000 0.000000 0.000000 0.000000 0.000000 1
3 11 11 0.000000 0.000000 0.000000 0.000000 0.000000 1
41 12 12 0.000000 0.000000 0.000000 0.000000 0.000000 1
13 1 1 1.000000 -1.000000 -1.000000 -1.000000 -1.000000 2
42 2 2 -1.000000 1.000000 -1.000000 -1.000000 -1.000000 2
61 3 3 -1.000000 -1.000000 1.000000 -1.000000 -1.000000 2
71 4 4 1.000000 1.000000 1.000000 -1.000000 -1.000000 2
91 5 5 1.000000 -1.000000 -1.000000 1.000000 1.000000 2
121 6 6 -1.000000 1.000000 -1.000000 1.000000 1.000000 2
141 7 7 -1.000000 -1.000000 1.000000 1.000000 1.000000 2
151 8 8 1.000000 1.000000 1.000000 1.000000 1.000000 2
111 9 9 0.000000 0.000000 0.000000 0.000000 0.000000 2
21 10 10 0.000000 0.000000 0.000000 0.000000 0.000000 2
31 11 11 0.000000 0.000000 0.000000 0.000000 0.000000 2
411 12 12 0.000000 0.000000 0.000000 0.000000 0.000000 2
16 1 1 -2.160247 0.000000 0.000000 0.000000 0.000000 3
22 2 2 2.160247 0.000000 0.000000 0.000000 0.000000 3
32 3 3 0.000000 -2.160247 0.000000 0.000000 0.000000 3
43 4 4 0.000000 2.160247 0.000000 0.000000 0.000000 3
5 5 5 0.000000 0.000000 -2.160247 0.000000 0.000000 3
62 6 6 0.000000 0.000000 2.160247 0.000000 0.000000 3
72 7 7 0.000000 0.000000 0.000000 -2.160247 0.000000 3
8 8 8 0.000000 0.000000 0.000000 2.160247 0.000000 3
92 9 9 0.000000 0.000000 0.000000 0.000000 -2.160247 3
10 10 10 0.000000 0.000000 0.000000 0.000000 2.160247 3
112 11 11 0.000000 0.000000 0.000000 0.000000 0.000000 3
122 12 12 0.000000 0.000000 0.000000 0.000000 0.000000 3
131 13 13 0.000000 0.000000 0.000000 0.000000 0.000000 3
142 14 14 0.000000 0.000000 0.000000 0.000000 0.000000 3
Data are stored in coded form using these coding formulas ...
A ~ A.as.is
B ~ B.as.is
C ~ C.as.is
D ~ D.as.is
E ~ E.as.is
Sometimes it’s difficult to pick the numbers of center points. We really want orthogonal blocking, but rotatable is also nice. Can we get an orthogonal design that is nearly rotatable?
ccd.pick()
evaluates different numbers of center
points for different blocks, computes the alpha
for orthogonal, and tries to find combinations
that are nearly orthogonal.
Here we just ask for three factors in cube and axial blocks. None of the designs is both rotatable and orthogonal, but the one that is closest has 9 center points in the cube block and 6 in the axial block. That could be too many, or the block size could be too large. The next one has only three total center points (2 and 1), which is probably too few.
> ccd.pick(3)
n.c n0.c blks.c n.s n0.s bbr.c wbr.s bbr.s N alpha.rot alpha.orth
1 8 9 1 6 6 1 1 1 29 1.681793 1.680336
2 8 2 1 6 1 1 1 1 17 1.681793 1.673320
3 8 6 1 6 4 1 1 1 24 1.681793 1.690309
4 8 5 1 6 3 1 1 1 22 1.681793 1.664101
5 8 10 1 6 7 1 1 1 31 1.681793 1.699673
6 8 8 1 6 5 1 1 1 27 1.681793 1.658312
7 8 3 1 6 2 1 1 1 19 1.681793 1.705606
8 8 7 1 6 5 1 1 1 26 1.681793 1.712698
9 8 4 1 6 2 1 1 1 20 1.681793 1.632993
10 8 4 1 6 3 1 1 1 21 1.681793 1.732051
You can get fancier and ask for the cube portion to be blocked further. Here was ask for two blocks of size 4 in the cube portion.
> ccd.pick(3,n.c=4,blks.c=2)
n.c n0.c blks.c n.s n0.s bbr.c wbr.s bbr.s N alpha.rot alpha.orth
1 4 1 2 6 1 1 1 1 17 1.681793 1.673320
2 4 6 2 6 8 1 1 1 34 1.681793 1.673320
3 4 3 2 6 4 1 1 1 24 1.681793 1.690309
4 4 5 2 6 7 1 1 1 31 1.681793 1.699673
5 4 4 2 6 5 1 1 1 27 1.681793 1.658312
6 4 7 2 6 10 1 1 1 38 1.681793 1.705606
7 4 7 2 6 9 1 1 1 37 1.681793 1.651446
8 4 2 2 6 2 1 1 1 20 1.681793 1.632993
9 4 2 2 6 3 1 1 1 21 1.681793 1.732051
10 4 4 2 6 6 1 1 1 28 1.681793 1.732051
rsm
also has a bbd()
function that
will generate Box-Behnken designs.
> bbd(4,block=FALSE,randomize=FALSE)
run.order std.order x1.as.is x2.as.is x3.as.is x4.as.is
1 1 1 -1 -1 0 0
2 2 2 1 -1 0 0
3 3 3 -1 1 0 0
4 4 4 1 1 0 0
5 5 5 0 0 -1 -1
6 6 6 0 0 1 -1
7 7 7 0 0 -1 1
8 8 8 0 0 1 1
9 9 9 -1 0 0 -1
10 10 10 1 0 0 -1
11 11 11 -1 0 0 1
12 12 12 1 0 0 1
13 13 13 0 -1 -1 0
14 14 14 0 1 -1 0
15 15 15 0 -1 1 0
16 16 16 0 1 1 0
17 17 17 -1 0 -1 0
18 18 18 1 0 -1 0
19 19 19 -1 0 1 0
20 20 20 1 0 1 0
21 21 21 0 -1 0 -1
22 22 22 0 1 0 -1
23 23 23 0 -1 0 1
24 24 24 0 1 0 1
25 25 25 0 0 0 0
26 26 26 0 0 0 0
27 27 27 0 0 0 0
28 28 28 0 0 0 0
Data are stored in coded form using these coding formulas ...
x1 ~ x1.as.is
x2 ~ x2.as.is
x3 ~ x3.as.is
x4 ~ x4.as.is
50.2 Potatoes data
From Eren and Kaymak-Erktekin (2007) via Lye (2019). We want optimum conditions to dehydrate potatoes via osmotic dehydration. This is a rotatable central composite design. There are several response variables, but we would like to maximize water loss (WL).
Variable | Unit | -2 | -1 | 0 | 1 | 2 |
---|---|---|---|---|---|---|
Temperature | \(^o\)C | 20 | 30 | 40 | 50 | 60 |
Sucrose conc. | % | 40 | 45 | 50 | 55 | 60 |
Salt conc. | % | 0 | 3.75 | 7.5 | 11.25 | 15 |
Time | min | 29.5 | 142 | 254.5 | 367 | 479.5 |
Read the data in. There are 16 factorial points, 8 axial points, and 7 center points for a total of 31 runs. They are entered in that order in the data set, although they would have been randomized when the experiment was run.
Then use coded.data()
to create
coded versions of the predictor variables within
the data frame.
> library(emmeans)
> potatoes <- read.csv("Potatoes.csv",header=TRUE)
> cpotatoes <- coded.data(potatoes,cTemp~(Temp-40)/10,cSucrose~(Sucrose-50)/5,cSalt~(Salt-7.5)/3.75,cTime~(Time-254.4)/112.5)
> cpotatoes
Temp Sucrose Salt Time WL SG WR aw
1 30 45 3.75 142.0 40.0 3.6 36.4 0.954
2 50 45 3.75 142.0 46.9 4.5 42.5 0.931
3 30 55 3.75 142.0 46.2 4.0 42.2 0.942
4 50 55 3.75 142.0 54.6 5.5 49.0 0.919
5 30 45 11.25 142.0 48.6 5.0 43.5 0.878
6 50 45 11.25 142.0 56.0 6.6 49.4 0.855
7 30 55 11.25 142.0 54.2 5.9 48.3 0.861
8 50 55 11.25 142.0 60.5 7.1 53.4 0.828
9 30 45 3.75 367.0 48.9 5.8 43.1 0.929
10 50 45 3.75 367.0 52.0 7.4 44.6 0.919
11 30 55 3.75 367.0 55.9 6.5 49.5 0.911
12 50 55 3.75 367.0 60.5 8.0 52.5 0.896
13 30 45 11.25 367.0 56.9 7.0 49.8 0.849
14 50 45 11.25 367.0 58.6 8.2 50.4 0.838
15 30 55 11.25 367.0 61.4 7.4 54.0 0.816
16 50 55 11.25 367.0 64.8 8.6 56.2 0.798
17 20 50 7.50 254.5 54.6 4.3 50.3 0.897
18 60 50 7.50 254.5 62.2 7.6 54.6 0.864
19 40 40 7.50 254.5 50.3 7.0 43.3 0.891
20 40 60 7.50 254.5 62.8 8.1 54.7 0.846
21 40 50 0.00 254.5 43.4 5.1 38.3 0.957
22 40 50 15.00 254.5 61.9 6.7 55.1 0.778
23 40 50 7.50 29.5 40.9 3.7 37.2 0.941
24 40 50 7.50 479.5 60.6 9.3 51.3 0.871
25 40 50 7.50 254.5 60.5 7.0 53.4 0.878
26 40 50 7.50 254.5 60.5 6.5 53.9 0.869
27 40 50 7.50 254.5 61.5 7.1 54.4 0.874
28 40 50 7.50 254.5 60.6 6.8 53.7 0.876
29 40 50 7.50 254.5 61.9 6.4 55.5 0.878
30 40 50 7.50 254.5 59.1 6.4 52.8 0.880
31 40 50 7.50 254.5 63.7 6.8 56.9 0.881
Data are stored in coded form using these coding formulas ...
cTemp ~ (Temp - 40)/10
cSucrose ~ (Sucrose - 50)/5
cSalt ~ (Salt - 7.5)/3.75
cTime ~ (Time - 254.4)/112.5
This is a central composite design, but let’s pretend for a moment that we only have the factorial points and some of the center points. That would be a good design for fitting a first order model.
Fit the first order model to this subset of the data.
The FO(a,b,c)
creates the first order model for
the variables a
, b
, and c
.
The lack of fit test is significant, meaning that there is enough curvature that the first order model is not sufficient. There are 12 df for lack of fit because there are 11 for interactions from the factorial portion of the data, plus 1 df for the contrast between the center points and the factorial points.
> use1 <- rep(TRUE,31)
> use1[17:24]<-FALSE # don't use axial points
> use1[29:31]<-FALSE # don't use last three center points
> WL.fo <- rsm(WL~FO(cTemp,cSucrose,cSalt,cTime),
+ data=cpotatoes,subset=use1)
> summary(WL.fo)
Call:
rsm(formula = WL ~ FO(cTemp, cSucrose, cSalt, cTime), data = cpotatoes,
subset = use1)
Estimate Std. Error t value Pr(>|t|)
(Intercept) 55.45211 0.75314 73.6282 < 2.2e-16 ***
cTemp 2.61250 0.84203 3.1026 0.0072783 **
cSucrose 3.13750 0.84203 3.7261 0.0020278 **
cSalt 3.50000 0.84203 4.1566 0.0008438 ***
cTime 3.25000 0.84203 3.8597 0.0015430 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Multiple R-squared: 0.7878, Adjusted R-squared: 0.7312
F-statistic: 13.92 on 4 and 15 DF, p-value: 6.168e-05
Analysis of Variance Table
Response: WL
Df Sum Sq Mean Sq F value Pr(>F)
FO(cTemp, cSucrose, cSalt, cTime) 4 631.71 157.926 13.921 6.168e-05
Residuals 15 170.16 11.344
Lack of fit 12 169.46 14.121 59.879 0.003106
Pure error 3 0.71 0.236
Direction of steepest ascent (at radius 1):
cTemp cSucrose cSalt cTime
0.4157757 0.4993287 0.5570201 0.5172330
Corresponding increment in original units:
Temp Sucrose Salt Time
4.157757 2.496644 2.088825 58.188708
Now we’ve added more runs to fill out the central composite design, so we can fit the second order model (OK, the data were already there, we’re just using all of them now).
The SO(a,b,c)
form gives the full second order model:
linear effects, pure quadratic effects, and two way interaction
effects (the cross products). An equivalent, but less compact,
form of the same model is FO(a,b,c)+TWI(a,b,c)+PQ(a,b,c)
yielding first order, two way interaction, and pure quadratic
effects.
> WL.so <- rsm(WL~SO(cTemp,cSucrose,cSalt,cTime),data=cpotatoes)
The summary()
gives us tons of information. The
first block is the usual summary for regression terms.
The first order terms are all highly significant.
One of the pure quadratics is marginally significant,
but the others are highly significant. Only
one of the two way interactions (cross product)
terms is even marginally significant.
The next block of output groups the terms into FO, TWI, and PQ, testing them as groups. It also does the lack of fit test, which is not significant for these data.
The third block prints the stationary point for the model surface in both coded and original variables. This point is reasonably within the experimental region, so it’s not an extrapolation.
The final block of output is the canonical analysis. The eigenvalues here are all negative, indicating a maximum. The eigenvectors (canonical directions) are each dominated by a single variable. This means that the canonical directions are fairly similar to the variable axes, which happens when the TWI terms are relatively small.
> summary(WL.so)
Call:
rsm(formula = WL ~ SO(cTemp, cSucrose, cSalt, cTime), data = cpotatoes)
Estimate Std. Error t value Pr(>|t|)
(Intercept) 61.11090 0.59429 102.8303 < 2.2e-16 ***
cTemp 2.37590 0.32095 7.4026 1.495e-06 ***
cSucrose 3.13321 0.32095 9.7622 3.835e-08 ***
cSalt 3.87540 0.32095 12.0747 1.881e-09 ***
cTime 3.81308 0.32095 11.8805 2.378e-09 ***
cTemp:cSucrose 0.22500 0.39309 0.5724 0.5750098
cTemp:cSalt -0.26250 0.39309 -0.6678 0.5137829
cTemp:cTime -1.01250 0.39309 -2.5758 0.0203150 *
cSucrose:cSalt -0.53750 0.39309 -1.3674 0.1904068
cSucrose:cTime 0.13750 0.39309 0.3498 0.7310529
cSalt:cTime -0.45000 0.39309 -1.1448 0.2691328
cTemp^2 -0.75565 0.29403 -2.5700 0.0205556 *
cSucrose^2 -1.21815 0.29403 -4.1429 0.0007646 ***
cSalt^2 -2.19315 0.29403 -7.4589 1.359e-06 ***
cTime^2 -2.66815 0.29403 -9.0743 1.042e-07 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Multiple R-squared: 0.9731, Adjusted R-squared: 0.9496
F-statistic: 41.4 on 14 and 16 DF, p-value: 7.048e-10
Analysis of Variance Table
Response: WL
Df Sum Sq Mean Sq F value Pr(>F)
FO(cTemp, cSucrose, cSalt, cTime) 4 1079.46 269.865 109.1573 2.125e-11
TWI(cTemp, cSucrose, cSalt, cTime) 6 26.48 4.413 1.7851 0.1657
PQ(cTemp, cSucrose, cSalt, cTime) 4 326.92 81.730 33.0589 1.499e-07
Residuals 16 39.56 2.472
Lack of fit 10 27.03 2.703 1.2944 0.3911
Pure error 6 12.53 2.088
Stationary point of response surface:
cTemp cSucrose cSalt cTime
1.3693434 1.3056026 0.5966421 0.4380645
Stationary point in original units:
Temp Sucrose Salt Time
53.693434 56.528013 9.737408 303.682260
Eigenanalysis:
eigen() decomposition
$values
[1] -0.6064319 -1.1566542 -2.1874299 -2.8846030
$vectors
[,1] [,2] [,3] [,4]
cTemp 0.95448498 0.1724550 0.0542147 0.237230786
cSucrose 0.18569220 -0.9420455 -0.2794036 0.001550546
cSalt -0.07927487 0.2524180 -0.9018507 0.341563955
cTime -0.21952749 -0.1381839 0.3250533 0.909424661
We can do predictions using emmeans()
. For example,
here are the predictions for various salt and sucrose
coded values when time and temperature are at the
level of the stationary point.
> emmeans(WL.so,~cSalt+cSucrose,at=list(cSalt=c(-1,0,1),cSucrose=c(-1,0,1),cTime=1.37,cTemp=.44))
NOTE: Results may be misleading due to involvement in interactions
cSalt cSucrose emmean SE df lower.CL upper.CL
-1 -1 51.1 1.218 16 48.5 53.7
0 -1 57.0 0.970 16 54.9 59.0
1 -1 58.5 1.218 16 55.9 61.0
-1 0 56.3 0.970 16 54.2 58.3
0 0 61.6 0.755 16 60.0 63.2
1 0 62.6 0.970 16 60.5 64.6
-1 1 59.0 1.218 16 56.4 61.6
0 1 63.8 0.970 16 61.8 65.9
1 1 64.2 1.218 16 61.6 66.8
Confidence level used: 0.95
Most of the two way interactions are not significant. Does it make sense to fit a reduced model including only the significant terms? Yes, because including non-significant terms in a model makes predictions worse. If you want more precise prediction, you don’t want the noise terms in your model.
See below how we can select just the terms we want.
If you look at the help for rsm()
you’ll see that
the TWI()
form allows a lot of flexibility (that we
don’t need here with just a single cross product
term).
> WL.red <- rsm(WL~FO(cTemp,cSucrose,cSalt,cTime)+PQ(cTemp,cSucrose,cSalt,cTime)+TWI(cTime,cTemp),data=cpotatoes)
> summary(WL.red)
Call:
rsm(formula = WL ~ FO(cTemp, cSucrose, cSalt, cTime) + PQ(cTemp,
cSucrose, cSalt, cTime) + TWI(cTime, cTemp), data = cpotatoes)
Estimate Std. Error t value Pr(>|t|)
(Intercept) 61.11090 0.58107 105.1695 < 2.2e-16 ***
cTemp 2.37590 0.31381 7.5710 1.970e-07 ***
cSucrose 3.13333 0.31381 9.9847 1.990e-09 ***
cSalt 3.87500 0.31381 12.3481 4.290e-11 ***
cTime 3.81308 0.31381 12.1507 5.782e-11 ***
cTemp^2 -0.75565 0.28749 -2.6284 0.0157053 *
cSucrose^2 -1.21815 0.28749 -4.2372 0.0003685 ***
cSalt^2 -2.19315 0.28749 -7.6285 1.750e-07 ***
cTime^2 -2.66815 0.28749 -9.2808 7.042e-09 ***
cTime:cTemp -1.01250 0.38434 -2.6344 0.0155014 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Multiple R-squared: 0.9663, Adjusted R-squared: 0.9518
F-statistic: 66.89 on 9 and 21 DF, p-value: 1.964e-13
Analysis of Variance Table
Response: WL
Df Sum Sq Mean Sq F value Pr(>F)
FO(cTemp, cSucrose, cSalt, cTime) 4 1079.46 269.865 114.1799 6.236e-14
PQ(cTemp, cSucrose, cSalt, cTime) 4 326.92 81.730 34.5800 5.814e-09
TWI(cTime, cTemp) 1 16.40 16.402 6.9399 0.0155
Residuals 21 49.63 2.364
Lack of fit 15 37.10 2.474 1.1847 0.4444
Pure error 6 12.53 2.088
Stationary point of response surface:
cTemp cSucrose cSalt cTime
1.2525887 1.2860982 0.8834306 0.4768896
Stationary point in original units:
Temp Sucrose Salt Time
52.52589 56.43049 10.81286 308.05008
Eigenanalysis:
eigen() decomposition
$values
[1] -0.6299144 -1.2181548 -2.1931548 -2.7938951
$vectors
[,1] [,2] [,3] [,4]
cTemp 0.9705122 0 0 0.2410519
cSucrose 0.0000000 -1 0 0.0000000
cSalt 0.0000000 0 -1 0.0000000
cTime -0.2410519 0 0 0.9705122
The results are very similar (the same coefficients, for example) in the reduced model. But the proof of the pudding is in the predictions. If you compare the predictions from the reduced model, they are somewhat different and have smaller standard errors.
> emmeans(WL.red,~cSalt+cSucrose,at=list(cSalt=c(-1,0,1),cSucrose=c(-1,0,1),cTime=1.37,cTemp=.44))
NOTE: Results may be misleading due to involvement in interactions
cSalt cSucrose emmean SE df lower.CL upper.CL
-1 -1 51.2 0.813 21 49.5 52.9
0 -1 57.3 0.771 21 55.7 58.9
1 -1 58.9 0.813 21 57.3 60.6
-1 0 55.5 0.771 21 53.9 57.2
0 0 61.6 0.739 21 60.1 63.2
1 0 63.3 0.771 21 61.7 64.9
-1 1 57.5 0.813 21 55.8 59.2
0 1 63.5 0.771 21 61.9 65.1
1 1 65.2 0.813 21 63.5 66.9
Confidence level used: 0.95
The lack of fit test says everything is OK, but what haven’t we done? We haven’t checked residuals!
So we do, and it looks bad. We are over estimating the high and low values and under estimating the middle values. (You can see a much weaker form of the same pattern in the full model residuals.)
> plot(WL.red,which=1)

Figure 50.1: Residual plot for reducted water loss model
Box-Cox picks the power 2, but to my great surprise the power 1 is just inside the interval. Things are a little better on the squared scale, but we’re going to continue with the original data to control the length of the example.
50.2.1 But wait! There’s more!
I hadn’t given you the full problem. The full problem is to maximize water loss while minimizing solid gain (SL). Let’s interpret “minimizing sold gain” to mean “keep solid gain under 6.5”.
We need to model solid gain. First fit the full model. Then fit a reduced model including only terms that are reasonably significant.
> SG.so <- rsm(SG~SO(cTemp,cSucrose,cSalt,cTime),data=cpotatoes)
> summary(SG.so)
Call:
rsm(formula = SG ~ SO(cTemp, cSucrose, cSalt, cTime), data = cpotatoes)
Estimate Std. Error t value Pr(>|t|)
(Intercept) 6.713252 0.133778 50.1822 < 2.2e-16 ***
cTemp 0.720817 0.072248 9.9769 2.836e-08 ***
cSucrose 0.295872 0.072248 4.0952 0.000845 ***
cSalt 0.571028 0.072248 7.9037 6.487e-07 ***
cTime 1.162642 0.072248 16.0923 2.651e-11 ***
cTemp:cSucrose 0.006250 0.088486 0.0706 0.944565
cTemp:cSalt -0.018750 0.088486 -0.2119 0.834862
cTemp:cTime 0.018750 0.088486 0.2119 0.834862
cSucrose:cSalt -0.031250 0.088486 -0.3532 0.728575
cSucrose:cTime -0.043750 0.088486 -0.4944 0.627722
cSalt:cTime -0.218750 0.088486 -2.4722 0.025036 *
cTemp^2 -0.217113 0.066188 -3.2802 0.004711 **
cSucrose^2 0.182887 0.066188 2.7631 0.013855 *
cSalt^2 -0.229613 0.066188 -3.4691 0.003164 **
cTime^2 -0.079613 0.066188 -1.2028 0.246542
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Multiple R-squared: 0.9676, Adjusted R-squared: 0.9392
F-statistic: 34.1 on 14 and 16 DF, p-value: 3.079e-09
Analysis of Variance Table
Response: SG
Df Sum Sq Mean Sq F value Pr(>F)
FO(cTemp, cSucrose, cSalt, cTime) 4 54.825 13.7062 109.4090 2.088e-11
TWI(cTemp, cSucrose, cSalt, cTime) 6 0.824 0.1373 1.0959 0.4064144
PQ(cTemp, cSucrose, cSalt, cTime) 4 4.164 1.0411 8.3102 0.0007954
Residuals 16 2.004 0.1253
Lack of fit 10 1.516 0.1516 1.8615 0.2307223
Pure error 6 0.489 0.0814
Stationary point of response surface:
cTemp cSucrose cSalt cTime
2.7205417 0.5983824 -7.1207671 17.2405071
Stationary point in original units:
Temp Sucrose Salt Time
67.20542 52.99191 -19.20288 2193.95705
Eigenanalysis:
eigen() decomposition
$values
[1] 0.18483027 -0.02185009 -0.21777830 -0.28865425
$vectors
[,1] [,2] [,3] [,4]
cTemp -0.00642643 -0.06564489 0.99643763 -0.05254991
cSucrose -0.99701518 -0.05701335 -0.01284677 -0.05044965
cSalt 0.01799251 0.46973620 -0.01547838 -0.88248773
cTime 0.07480440 -0.87851473 -0.08189887 -0.46465982
> SG.red <- rsm(SG~FO(cTemp,cSucrose,cSalt,cTime)+PQ(cTemp,cSucrose,cSalt)+TWI(cSalt,cTime),data=cpotatoes)
> summary(SG.red)
Call:
rsm(formula = SG ~ FO(cTemp, cSucrose, cSalt, cTime) + PQ(cTemp,
cSucrose, cSalt) + TWI(cSalt, cTime), data = cpotatoes)
Estimate Std. Error t value Pr(>|t|)
(Intercept) 6.631945 0.104163 63.6688 < 2.2e-16 ***
cTemp 0.720833 0.065189 11.0576 1.880e-10 ***
cSucrose 0.295833 0.065189 4.5381 0.0001621 ***
cSalt 0.571028 0.065189 8.7596 1.268e-08 ***
cTime 1.162500 0.065189 17.8328 1.443e-14 ***
cTemp^2 -0.208644 0.059382 -3.5136 0.0019592 **
cSucrose^2 0.191356 0.059382 3.2225 0.0039190 **
cSalt^2 -0.221144 0.059382 -3.7241 0.0011793 **
cSalt:cTime -0.218750 0.079840 -2.7399 0.0119567 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Multiple R-squared: 0.9637, Adjusted R-squared: 0.9505
F-statistic: 73.01 on 8 and 22 DF, p-value: 4.738e-14
Analysis of Variance Table
Response: SG
Df Sum Sq Mean Sq F value Pr(>F)
FO(cTemp, cSucrose, cSalt, cTime) 4 54.825 13.7062 134.3884 4.015e-15
PQ(cTemp, cSucrose, cSalt) 3 3.983 1.3277 13.0177 4.210e-05
TWI(cSalt, cTime) 1 0.766 0.7656 7.5069 0.01196
Residuals 22 2.244 0.1020
Lack of fit 16 1.755 0.1097 1.3472 0.37622
Pure error 6 0.489 0.0814
Stationary point of response surface:
cTemp cSucrose cSalt cTime
1.7274272 -0.7729905 5.3142857 -8.1344592
Stationary point in original units:
Temp Sucrose Salt Time
57.27427 46.13505 27.42857 -660.72666
Eigenanalysis:
eigen() decomposition
$values
[1] 0.19135638 0.04495637 -0.20864362 -0.26609999
$vectors
[,1] [,2] [,3] [,4]
cTemp 0 0.0000000 1 0.0000000
cSucrose 1 0.0000000 0 0.0000000
cSalt 0 -0.3801685 0 0.9249173
cTime 0 0.9249173 0 0.3801685
What we’re going to do is make some contour plots of our WL predictions and overlay contour plots of our SG predictions. Then we look at how well we can do. If we’re lucky, high values of WL are in the same regions as low values of SG, but we’re not lucky.
In this first plot we look at the two variables as time and temperature vary, with sucrose and salt set to the stationary point values for WL. The black countours are WL; the red contours are SG. The SL=6.5 contour curves from the southeast up to the northwest. At one point it gets between the 60 and 65 contours of WL. The best we can do in this slice is about 61 for WL.
> contour(WL.red,cTime~cTemp,at=c(cSalt=.88,cSucrose=1.29))
> contour(SG.red,cTime~cTemp,at=c(cSalt=.88,cSucrose=1.29),add=TRUE,col="red")

Figure 50.2: Contour plot of SG overlaid on WL; by Time and Temperature
Now flip around and look at the results by salt and sucrose with time and temperature set to the stationary point levels.
We can clearly see the saddle point nature of SG in this slice. The 6.5 contour for SG is way over in the west southwest. It does get up between 50 and 55 on WL, but no higher. Clearly, the other point we found is better.
> contour(WL.red,cSucrose~cSalt,at=c(cTime=.48,cTemp=1.25))
> contour(SG.red,cSucrose~cSalt,at=c(cTime=.48,cTemp=1.25),add=TRUE,col="red")

Figure 50.3: Contour plot of SG overlaid on WL; by Sucrose and Salt
50.3 Conversion data
These are data taken from the examples in the Design-Expert manual. We are making a product that should have an activity of 63; we could accept the range 60 to 66, but 63 is more desirable. There are three process variables: time (minutes), temperature (degrees C), and catalyst (percent). When we do our experiment, we have two responses: the conversion factor (essentially the percent yield) and the activity.
Our goal is to get process variables so that we hit an activity of 63 (or close enough) and maximize the conversion. In addition, we want to operate within the “cube” part of the experimental region.
Read in the data and covert to coded units.
> conversion <- read.csv("Conversion.csv",header=TRUE)
> conversion
Time Temperature Catalyst Conversion Activity
1 45.00 85.00 2.50 75 60.4
2 45.00 85.00 2.50 83 60.6
3 45.00 85.00 2.50 76 59.1
4 45.00 85.00 2.50 81 59.2
5 45.00 85.00 2.50 91 58.9
6 45.00 85.00 2.50 80 60.8
7 53.41 85.00 2.50 79 65.9
8 36.59 85.00 2.50 76 53.6
9 45.00 85.00 1.66 55 57.4
10 45.00 85.00 3.34 81 63.2
11 45.00 93.41 2.50 97 60.7
12 45.00 76.59 2.50 85 60.0
13 40.00 80.00 2.00 74 53.2
14 40.00 90.00 3.00 66 59.8
15 50.00 90.00 2.00 70 62.6
16 40.00 90.00 2.00 88 53.4
17 50.00 90.00 3.00 97 67.8
18 50.00 80.00 3.00 90 67.9
19 40.00 80.00 3.00 71 57.3
20 50.00 80.00 2.00 51 62.9
> c.conversion <- coded.data(conversion,cTime~(Time-45)/5,cTemp~(Temperature-85)/5,
+ cCat~(Catalyst-2.5)/.5)
Let’s fit conversion starting with a second order model. Looking at the summary, one pure quadratic is not even marginally significant; one cross product is not marginally significant; and one first order term is not marginally significant. However, all variables need to be retained due to hierarchy.
> conversion.so <- rsm(Conversion~SO(cTime,cTemp,cCat),data=c.conversion)
> summary(conversion.so)
Call:
rsm(formula = Conversion ~ SO(cTime, cTemp, cCat), data = c.conversion)
Estimate Std. Error t value Pr(>|t|)
(Intercept) 81.0917 1.9233 42.1632 1.351e-12 ***
cTime 1.0284 1.2760 0.8059 0.4390217
cTemp 4.0403 1.2760 3.1664 0.0100492 *
cCat 6.2060 1.2766 4.8612 0.0006602 ***
cTime:cTemp 2.1250 1.6673 1.2745 0.2312930
cTime:cCat 11.3750 1.6673 6.8225 4.611e-05 ***
cTemp:cCat -3.8750 1.6673 -2.3242 0.0424712 *
cTime^2 -1.8311 1.2419 -1.4744 0.1711487
cTemp^2 2.9407 1.2419 2.3678 0.0394226 *
cCat^2 -5.2027 1.2443 -4.1811 0.0018843 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Multiple R-squared: 0.92, Adjusted R-squared: 0.8479
F-statistic: 12.77 on 9 and 10 DF, p-value: 0.0002207
Analysis of Variance Table
Response: Conversion
Df Sum Sq Mean Sq F value Pr(>F)
FO(cTime, cTemp, cCat) 3 762.93 254.31 11.4357 0.001428
TWI(cTime, cTemp, cCat) 3 1191.37 397.12 17.8576 0.000243
PQ(cTime, cTemp, cCat) 3 601.51 200.50 9.0160 0.003414
Residuals 10 222.38 22.24
Lack of fit 5 56.38 11.28 0.3397 0.869479
Pure error 5 166.00 33.20
Stationary point of response surface:
cTime cTemp cCat
-1.0179845 -0.5294920 -0.3192329
Stationary point in original units:
Time Temperature Catalyst
39.910077 82.352540 2.340384
Eigenanalysis:
eigen() decomposition
$values
[1] 3.408437 2.322758 -9.824322
$vectors
[,1] [,2] [,3]
cTime 0.1380685 0.7982920 0.5862312
cTemp -0.9428308 0.2872127 -0.1690534
cCat 0.3033270 0.5293758 -0.7923093
Refit omitting the noise terms. Notice how you can add
just certain two way interactions in the TWI()
form.
Residuals seem OK.
> conversion.red <- rsm(Conversion~FO(cTime,cTemp,cCat)+PQ(cTemp,cCat)+TWI(formula=~cTime*cCat+cTemp*cCat),data=c.conversion)
> summary(conversion.red)
Call:
rsm(formula = Conversion ~ FO(cTime, cTemp, cCat) + PQ(cTemp,
cCat) + TWI(formula = ~cTime * cCat + cTemp * cCat), data = c.conversion)
Estimate Std. Error t value Pr(>|t|)
(Intercept) 79.5933 1.7509 45.4589 8.375e-15 ***
cTime 1.0284 1.3683 0.7516 0.4667826
cTemp 4.0403 1.3683 2.9529 0.0120800 *
cCat 6.2060 1.3690 4.5334 0.0006857 ***
cTemp^2 3.1226 1.3252 2.3564 0.0362905 *
cCat^2 -5.0213 1.3278 -3.7817 0.0026153 **
cTime:cCat 11.3750 1.7878 6.3624 3.597e-05 ***
cCat:cTemp -3.8750 1.7878 -2.1674 0.0510298 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Multiple R-squared: 0.8896, Adjusted R-squared: 0.8251
F-statistic: 13.81 on 7 and 12 DF, p-value: 7.061e-05
Analysis of Variance Table
Response: Conversion
Df Sum Sq Mean Sq F value
FO(cTime, cTemp, cCat) 3 762.93 254.31 9.9453
PQ(cTemp, cCat) 2 553.16 276.58 10.8163
TWI(formula = ~cTime * cCat + cTemp * cCat) 2 1155.25 577.62 22.5891
Residuals 12 306.85 25.57
Lack of fit 7 140.85 20.12 0.6061
Pure error 5 166.00 33.20
Pr(>F)
FO(cTime, cTemp, cCat) 0.001418
PQ(cTemp, cCat) 0.002063
TWI(formula = ~cTime * cCat + cTemp * cCat) 8.545e-05
Residuals
Lack of fit 0.736262
Pure error
Stationary point of response surface:
cTime cTemp cCat
-0.86490534 -0.70305722 -0.09040788
Stationary point in original units:
Time Temperature Catalyst
40.675473 81.484714 2.454796
Eigenanalysis:
eigen() decomposition
$values
[1] 4.588664 2.460161 -8.947536
$vectors
[,1] [,2] [,3]
cTime 0.5989293 -0.5989112 0.5315910
cTemp -0.6385820 -0.7577546 -0.1342427
cCat 0.4832150 -0.2590625 -0.8362953
> plot(conversion.red,which=1)

Figure 50.4: Residual plot of reduced conversion model
Now fit activity. The pure quadratics and interactions are not significant.
> activity.so <- rsm(Activity~SO(cTime,cTemp,cCat),data=c.conversion)
> summary(activity.so)
Call:
rsm(formula = Activity ~ SO(cTime, cTemp, cCat), data = c.conversion)
Estimate Std. Error t value Pr(>|t|)
(Intercept) 59.81833 0.42544 140.6029 < 2.2e-16 ***
cTime 4.26033 0.28226 15.0935 3.293e-08 ***
cTemp 0.25460 0.28226 0.9020 0.3883
cCat 2.23118 0.28240 7.9008 1.313e-05 ***
cTime:cTemp -0.38750 0.36881 -1.0507 0.3181
cTime:cCat -0.03750 0.36881 -0.1017 0.9210
cTemp:cCat 0.31250 0.36881 0.8473 0.4166
cTime^2 0.06769 0.27472 0.2464 0.8104
cTemp^2 0.27977 0.27472 1.0184 0.3325
cCat^2 0.26294 0.27525 0.9553 0.3620
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Multiple R-squared: 0.9672, Adjusted R-squared: 0.9376
F-statistic: 32.74 on 9 and 10 DF, p-value: 2.955e-06
Analysis of Variance Table
Response: Activity
Df Sum Sq Mean Sq F value Pr(>F)
FO(cTime, cTemp, cCat) 3 316.71 105.571 97.0169 1.079e-07
TWI(cTime, cTemp, cCat) 3 1.99 0.665 0.6107 0.6232
PQ(cTime, cTemp, cCat) 3 1.94 0.645 0.5931 0.6336
Residuals 10 10.88 1.088
Lack of fit 5 7.23 1.446 1.9786 0.2359
Pure error 5 3.65 0.731
Stationary point of response surface:
cTime cTemp cCat
15.05168 17.57008 -13.61031
Stationary point in original units:
Time Temperature Catalyst
120.258424 172.850398 -4.305155
Eigenanalysis:
eigen() decomposition
$values
[1] 0.48931661 0.18268588 -0.06160297
$vectors
[,1] [,2] [,3]
cTime 0.3691679 0.4710760 0.8011258
cTemp -0.7502855 -0.3576299 0.5560328
cCat -0.5484402 0.8063425 -0.2214160
Refit with the first order and also omitting unneeded first order term. The reduced model is acceptable from an ANOVA perspective, and residuals look OK.
> activity.fo <- rsm(Activity~FO(cTime,cTemp,cCat),data=c.conversion)
> activity.red <- rsm(Activity~FO(cTime,cCat),data=c.conversion)
> anova(activity.red,activity.fo,activity.so)
Analysis of Variance Table
Model 1: Activity ~ FO(cTime, cCat)
Model 2: Activity ~ FO(cTime, cTemp, cCat)
Model 3: Activity ~ FO(cTime, cTemp, cCat) + TWI(cTime, cTemp, cCat) +
PQ(cTime, cTemp, cCat)
Res.Df RSS Df Sum of Sq F Pr(>F)
1 17 15.697
2 16 14.812 1 0.8853 0.8136 0.3883
3 10 10.882 6 3.9299 0.6019 0.7240
> plot(activity.red,which=1)

Figure 50.5: Residual plot for first order activity model
We have our models, now we work on the experimental goal. We are going to approach this graphically. For each level of cTemp from -1 to 1 by increments of .25 we will:
- Make a contour plot of activity with contours at levels 60, 63, and 66 (the minimum, desired, and maximum levels);
- Overlay a contour plot of conversion (red, to make it visible);
- Overlay a box (in green) to show the experimental area.
The goal is to walk through the images and find where we want to operate the system. This will involve some tradeoff between being as close as possible to 63 for activity and having the highest conversion.
Note: you can do all these within a loop, but I’ve kept them separate to get better alt text for the graphics.
> contour(activity.red,cTime~cCat,levels=c(60,63,66))
> contour(conversion.red,cTime~cCat,at=list(cTemp=-1),add=TRUE,col="red")
> lines(c(3,3,2,2,3),c(40,50,50,40,40),col="green")

Figure 50.6: Overlaid contours of activity and conversion with coded temperature at -1.
> contour(activity.red,cTime~cCat,levels=c(60,63,66))
> contour(conversion.red,cTime~cCat,at=list(cTemp=-.75),add=TRUE,col="red")
> lines(c(3,3,2,2,3),c(40,50,50,40,40),col="green")

Figure 50.7: Overlaid contours of activity and conversion with coded temperature at -.75.
> contour(activity.red,cTime~cCat,levels=c(60,63,66))
> contour(conversion.red,cTime~cCat,at=list(cTemp=-.5),add=TRUE,col="red")
> lines(c(3,3,2,2,3),c(40,50,50,40,40),col="green")

Figure 50.8: Overlaid contours of activity and conversion with coded temperature at -.5.
> contour(activity.red,cTime~cCat,levels=c(60,63,66))
> contour(conversion.red,cTime~cCat,at=list(cTemp=-.25),add=TRUE,col="red")
> lines(c(3,3,2,2,3),c(40,50,50,40,40),col="green")

Figure 50.9: Overlaid contours of activity and conversion with coded temperature at -.25.
> contour(activity.red,cTime~cCat,levels=c(60,63,66))
> contour(conversion.red,cTime~cCat,at=list(cTemp=0),add=TRUE,col="red")
> lines(c(3,3,2,2,3),c(40,50,50,40,40),col="green")

Figure 50.10: Overlaid contours of activity and conversion with coded temperature at 0.
> contour(activity.red,cTime~cCat,levels=c(60,63,66))
> contour(conversion.red,cTime~cCat,at=list(cTemp=.25),add=TRUE,col="red")
> lines(c(3,3,2,2,3),c(40,50,50,40,40),col="green")

Figure 50.11: Overlaid contours of activity and conversion with coded temperature at .25.
> contour(activity.red,cTime~cCat,levels=c(60,63,66))
> contour(conversion.red,cTime~cCat,at=list(cTemp=.5),add=TRUE,col="red")
> lines(c(3,3,2,2,3),c(40,50,50,40,40),col="green")

Figure 50.12: Overlaid contours of activity and conversion with coded temperature at .5.
> contour(activity.red,cTime~cCat,levels=c(60,63,66))
> contour(conversion.red,cTime~cCat,at=list(cTemp=.75),add=TRUE,col="red")
> lines(c(3,3,2,2,3),c(40,50,50,40,40),col="green")

Figure 50.13: Overlaid contours of activity and conversion with coded temperature at .75.
> contour(activity.red,cTime~cCat,levels=c(60,63,66))
> contour(conversion.red,cTime~cCat,at=list(cTemp=1),add=TRUE,col="red")
> lines(c(3,3,2,2,3),c(40,50,50,40,40),col="green")

Figure 50.14: Overlaid contours of activity and conversion with coded temperature at 1.
In general, we can achieve higher conversions and stay within our activity bounds as we move to higher levels of temperature. If we stick to activity right at 63, we can almost hit conversion of 90. If we allow ourselves an activity of 66, we can hit a conversion of about 95.
50.4 Constrained maximization
Optimizing a quadratic function subject to linear equality and/or inequality constraints is called quadratic programming. If you want to maximize the function, use -2 times matrix of second order terms and the linear terms as the first two arguments. If you want to minimize, use 2B and -b instead.
The constraints are “greater than or equal to” constraints. The third argument has a column of coefficients for every constraint, and the four argument has the lower bound. Here we’re saying that the first argument must be at least -100, which seems safe enough. The constrained solution is the same as the unconstrained solution, because the unconstrained solution meets the constraint.
Since we are maximizing, the value at the optimum
is the Intercept from the model minus the
value
reported here, in this case,
61.1- (-6.1) = 67.2.
> library(quadprog)
> solve.QP(-2*WL.red$B,WL.red$b,matrix(c(1,0,0,0),ncol=1),-100)
$solution
[1] 1.2525887 1.2860982 0.8834306 0.4768896
$value
[1] -6.123755
$unconstrained.solution
[1] 1.2525887 1.2860982 0.8834306 0.4768896
$iterations
[1] 1 0
$Lagrangian
[1] 0
$iact
[1] 0
OK, let’s make it bite a little bit. Let’s constrain the third variable to be at least 1, and the first variable to be no more than one (its negative no more than -1).
> solve.QP(-2*WL.red$B,WL.red$b,matrix(c(0,0,1,0,-1,0,0,0),ncol=2),c(1,-1))
$solution
[1] 1.0000000 1.2860982 1.0000000 0.5248153
$value
[1] -6.05187
$unconstrained.solution
[1] 1.2525887 1.2860982 0.8834306 0.4768896
$iterations
[1] 3 0
$Lagrangian
[1] 0.5113095 0.3332150
$iact
[1] 2 1
The first meq
of the constraints are equality constraints. It defaults to zero, but here we say that the first variable
must equal 1.1.
> solve.QP(-2*WL.red$B,WL.red$b,matrix(c(1,0,0,0),ncol=1),1.1,meq=1)
$solution
[1] 1.1000000 1.2860982 0.8834306 0.5058415
$value
[1] -6.108397
$unconstrained.solution
[1] 1.2525887 1.2860982 0.8834306 0.4768896
$iterations
[1] 2 0
$Lagrangian
[1] 0.201295
$iact
[1] 1
You can use more interesting constraints like the first variable minus the second variable must be at least .5.
> solve.QP(-2*WL.red$B,WL.red$b,matrix(c(1,-1,0,0),ncol=1),.5)
$solution
[1] 1.5986920 1.0986920 0.8834306 0.4112207
$value
[1] -6.00196
$unconstrained.solution
[1] 1.2525887 1.2860982 0.8834306 0.4768896
$iterations
[1] 2 0
$Lagrangian
[1] 0.4565795
$iact
[1] 1