Phương pháp thiết kế phối hợp có tâm

Trong nhóm các phương pháp thiết kế thí nghiệm theo hướng xây dựng bề mặt đáp ứng, phương pháp phối hợp có tâm (central composite design - CCD) được dùng phổ biến nhất. Một thí nghiệm khảo sát ảnh hưởng của `k` yếu tố theo phương pháp này là sự phối hợp của ba nhóm nghiệm thức sau (minh họa trên Hình 1 và Hình 2):

Các nghiệm thức theo cách thiết kế kết hợp yếu tố (có thể là kết hợp đủ `2^k` hay kết hợp giảm `2^(k-p)`, trong đó mỗi yếu tố được thử nghiệm ở hai mức, mức thấp và mức cao. Như vậy nếu có `k` yếu tố kết hợp đủ sẽ có `2^k` nghiệm thức thuộc nhóm này (hình vuông xanh trên Hình 1 và Hình 2).
Các nghiệm thức “sao” (hay nghiệm thức trục) : các nghiệm thức trong nhóm này, thường được ký hiệu là "*", có giá trị của `(k – 1)` yếu tố là giá trị tâm, yếu tố còn lại có giá trị đã được mã hóa là `–alpha` hay `alpha` với `alpha >= 1`. Như vậy sẽ có `2k` nghiệm thức thuộc nhóm này (vòng tròn đỏ trên Hình 1 và Hình 2). Giá trị `alpha` còn được gọi là giá trị trục, giá trị “sao”, hay cánh tay đòn.
Nghiệm thức tâm : giá trị của tất cả các yếu tố đều ở mức tâm (vòng tròn xanh lá cây trên Hình 1 và Hình 2).
Trong khi thực hiện thí nghiệm, các nghiệm thức thuộc nhóm 1 và 2 đều không cần lặp lại, riêng nghiệm thức tâm được lặp lại `n` lần.

Hình 1 Minh họa thiết kế thí nghiệm phối hợp có tâm hai yếu tố

Hình 2 Minh họa thiết kế thí nghiệm phối hợp có tâm ba yếu tố

Như vậy một thí nghiệm thiết kế theo phương pháp phối hợp có tâm để khảo sát ảnh hưởng của `k` yếu tố sẽ gồm `2^k+2k+1` nghiệm thức, và nếu chỉ lặp ở nghiệm thức tâm sẽ gồm `2^k+2k+n` đơn vị thí nghiệm (khi dùng kết hợp đủ).

Tùy theo sự lụa chọn các thông số, thí nghiệm phối hợp có tâm có thể có các tính chất sau:

Nhìn chung mỗi yếu tố đều được khảo sát ở 5 mức. Riêng khi chọn `alpha=1`, mỗi yếu tố chỉ còn được khảo sát ở 3 mức.
Tính trực giao (orthogonality) : khi các vectơ cột của ma trận yếu tố mã hóa trực giao với nhau từng đôi một.
Tính quay được (rotatability) : Thuật ngữ này do Box và Hunter đề xuất vào 1957. Theo đó thiết kế có tính quay được khi phương sai theo mọi phương đều như nhau, nghĩa là nếu hai nghiệm thức có cùng khoảng cách đến tâm (trong không gian các yếu tố đã được mã hóa) thì phương sai của đáp ứng `Y` của hai nghiệm thức đó giống nhau. Để thỏa mãn tính chất này ta cần có:
`alpha=2^(k//4)`(1)
Tính cầu (sphericity) : Thiết kế CCD được xem là có tính cầu khi tất cả các nghiệm thức (trừ nghiệm thức tâm) đều nằm trên một mặt cầu trong khống gian yếu tố đã được mã hóa. Để đáp ứng điều này, ta chọn:
`alpha=sqrt k`(2)

Một trong các lợi thế của thiết kế phối hợp có tâm là ta có thể chia quá trình thí nghiệm làm một số giai đoạn. Sau mỗi giai đoạn, xử lý dữ liệu thu được và dựa vào kết quả xử lý để quyết định hành động tiếp theo. Thí dụ như tiến trình sau cho thí nghiệm khảo sát tác động của `k` yếu tố:

Giai đoạn 1 : thực hiện `2^k` (hay `2^(k - p)`) nghiệm thức nhóm kết hợp yếu tố.
Xử lý dữ liệu để đánh giá mức độ tác động của các yếu tố. Dựa vào kết quả xử lý, có thể loại bỏ một số yếu tố.
Giai đoạn 2 : thực hiện nghiệm thức tâm.
Xử lý dữ liệu để đánh giá tính phi tuyến của tác động của các yếu tố đến đáp ứng. Nếu tác động chỉ là tuyến tính, ta có thể dừng thí nghiệm ở bước này.
Giai đoạn 3 : thực hiện các nghiệm thức trục.
Xử lý dữ liệu để thu được phương trình hồi quy, xác định được các thông số tối ưu.

Qua đó, ta thấy giai đoạn một có thể dùng như một thí nghiệm sàng lọc để chỉ giữ lại các yếu tố có tác động đáng kể đến đáp ứng. Tùy theo kết quả xử lý số liệu mà ta có thể tiếp tục thí nghiệm một cách hợp lý và hiệu quả hơn hay thậm chí ngừng tiến trình thí nghiệm.

Trong trường hợp đã biết các yếu tố khảo sát thực sự có ảnh hưởng đến đáp ứng, ta có thể chia quá trình thí nghiệm làm hai giai đoạn:

Giai đoạn 1 : thực hiện các nghiệm thức thuộc nhóm kết hợp yếu tố và một số đơn vị của nghiệm thức tâm.
Xử lý dữ liệu thu được để đánh giá tính phi tuyến. Nếu tác động chỉ là tuyến tính, có thể dừng thí nghiệm.
Giai đoạn 2 : thực hiện các nghiệm thức trục và một số đơn vị của nghiệm thức tâm.

Như ta thấy, cách thực hiện qua một số giai đoạn giúp tiến trình thí nghiệm được định hướng tốt hơn, và trong một số trường hợp có thể tiết kiệm được nguồn lực và thời gian. Vì thế cách thực hiện này được áp dụng rộng rãi khi có nhiều yếu tố, đặc biệt là trong công nghiệp.

Thí nghiệm nhằm khảo sát tác động của nhiệt độ, pH và tỷ lệ dung môi/nguyên liệu (DN) lên tỷ lệ thu hồi hoạt chất B từ nguyên liệu A. Qua khảo sát ban đầu ta chọn mức thấp và mức cao của nhiệt độ là 60 và 80°C, pH là 5 và 6, tỷ lệ DN là 10 và 20.

Trong thí nghiệm này, ta sử dụng phụ kiện rsm của R để thiết kế thí nghiệm và xử lý số liệu thu được.

Thiết kế thí nghiệm

Ta dùng hàm ccd trong phụ kiện rsm để thiết kế thí nghiệm với các câu lệnh sau:.

library(rsm)
ccd(basis = ~ nhd + ph + tldgmngl, n0 = c(4, 2), alpha = 1.6, randomize = FALSE)
    coding = list(nhd ~ (NhietDo-70)/10, ph ~ (pH-5.5)/0.5, tldgmngl ~ (DN-15)/5))

Trong các câu lệnh trên :

library(rsm) dùng để kích hoạt phụ kiện rsm (giả sử phụ kiện này đã được tải về và cài đặt vào máy tính).
ccd( . . . ) là hàm dùng để thiết kế thí nghiệm phối hợp có tâm.
basis là đối số đầu tiên dùng để khai báo tên của yếu tố mã hóa hay số yếu tố. Trong thí dụ này ta khai báo 3 biến mã hóa: nhd cho nhiệt độ, ph cho pH, tldgmngl cho tỷ lệ dung môi / nguyên liệu.
coding = list( . . . ) dùng để khai báo tương quan giữa biến mã hóa (thí dụ nhd) và biến chưa mã hóa (thí dụ NhietDo)
alpha = 1,6 dùng để khai báo giá trị trục. Giá trị mặc định là "orthogonal", trong trường hợp 3 yếu tố là 1,8257... Ở đây ta chọn là 1,6 để giá trị các yếu tố thí nghiệm tương đối chẵn, dễ thực hiện hơn.

Các đối số khác sẽ được giải thích trong phần sau.

Cho thực hiện các câu lệnh trên, ta thu được kết quả sau:

> ccd(basis = ~ nhd + ph + tldgmngl, n0 = c(4, 2), alpha = 1.6, randomize = FALSE,
+    coding = list(nhd ~ (NhietDo - 70)/10, ph ~ (pH - 5.5)/0.5, tldgmngl ~ (DN - 15)/5))
   run.order std.order NhietDo  pH DN Block
1          1         1      60 5.0 10     1
2          2         2      80 5.0 10     1
3          3         3      60 6.0 10     1
4          4         4      80 6.0 10     1
5          5         5      60 5.0 20     1
6          6         6      80 5.0 20     1
7          7         7      60 6.0 20     1
8          8         8      80 6.0 20     1
9          9         9      70 5.5 15     1
10        10        10      70 5.5 15     1
11        11        11      70 5.5 15     1
12        12        12      70 5.5 15     1
13         1         1      54 5.5 15     2
14         2         2      86 5.5 15     2
15         3         3      70 4.7 15     2
16         4         4      70 6.3 15     2
17         5         5      70 5.5  7     2
18         6         6      70 5.5 23     2
19         7         7      70 5.5 15     2
20         8         8      70 5.5 15     2

Data are stored in coded form using these coding formulas ...
nhd ~ (NhietDo - 70)/10
ph ~ (pH - 5.5)/0.5
tldgmngl ~ (DN - 15)/5

Trong kết quả trên, ta thấy một bảng với các chi tiết sau:

Bảng có 7 cột : cột đầu tiên dùng để đánh số thứ tự các dòng của bảng, cột run.order cho biết thứ tự thực hiện các đơn vị thí nghiệm, cột std.order cho biết thứ tự sắp xếp chuẩn các đơn vị ấy dựa vào giá trị của các biến, 3 cột tương ứng với 3 yếu tố chưa mã hóa, và cột Block.
Bảng có 21 dòng : dòng đầu tiên là dòng tiêu đề, 20 dòng còn lại tương ứng với 20 đơn vị thí nghiệm.
Cột Block cho ta biết thí nghiệm được chia làm hai khối, ký hiệu Block 1 và Block 2. Các đơn vị thí nghiệm của hai khối này được đánh số thứ tự riêng.
Khối đầu tiên (Block 1) có 12 đơn vị thí nghiệm gồm 8 thí nghiệm tương ứng với 8 nghiệm thức của nhóm kết hợp yếu tố đủ (số thứ tự 1 đến 8) và 4 thí nghiệm của nghiệm thức tâm (số thứ tự 9 đến 12). Trong R, nhóm này được gọi là nhóm "cube".
Khối thứ hai (Block 2) có 8 đơn vị thí nghiệm gồm 6 thí nghiệm tương ứng với 6 nghiệm thức trục (số thứ tự 1 đến 6) và 2 thí nghiệm của nghiệm thức tâm (số thứ tự 7 và 8). Trong R, nhóm này được gọi là nhóm "star".
Do đối số randomize của hàm ccd được chọn là FALSE nên thứ tự thực hiện các đơn vị trùng với thứ tự sắp xếp chuẩn, dễ theo dõi hơn.
Khai báo n0 = c(4, 2) trong hàm ccd có nghĩa là nghiệm thức tâm gồm tổng cộng 6 đơn vị thí nghiệm: 4 đơn vị cho khối 1 và 2 đơn vị cho khối 2.
Thiết kế này cho phép ta chia thí nghiệm làm hai giai đoạn. Giai đoạn đầu thực hiện khối 1, xử lý dữ liệu, đánh giá tác động của các yếu tố và tính phi tuyến. Trên cơ sở đó quyết định có thực hiện tiếp khối 2 hay không.

Xử lý số liệu

Sau khi thực hiện thí nghiệm trên, ta thu được kết quả của tỷ lệ thu hồi và được trình bày ở Bảng 1 sau:

Bảng 1 Tỷ lệ thu hồi thu được từ thí nghiệm
Khối	Đơn vị thí nghiệm	Nhiệt độ (°C)	pH	Tỷ lệ DN	Tỷ lệ thu hồi (%)
1	1	60	5,0	10	45
1	2	80	5,0	10	93
1	3	60	6,0	10	65
1	4	80	6,0	10	95
1	5	60	5,0	20	34
1	6	80	5,0	20	83
1	7	60	6,0	20	42
1	8	80	6,0	20	88
1	9	70	5,5	15	90
1	10	70	5,5	15	85
1	11	70	5,5	15	84
1	12	70	5,5	15	90
2	1	54	5,5	15	42
2	2	86	5,5	15	85
2	3	70	4,7	15	53
2	4	70	6,3	15	80
2	5	70	5,5	7	94
2	6	70	5,5	23	51
2	7	70	5,5	15	90
2	8	70	5,5	15	92

Để xử lý kết quả này, ta chuẩn bị tập tin TrichLy.csv, chuyển vào R, đặt tên cho bảng dữ liệu là trly và xử lý bằng hàm rsm của phụ kiện rsm với các câu lệnh sau:

kq <- rsm(ThuHoi ~ SO(NhietDo, pH, DN), data = trly)
summary(kq)

Trong đoạn lệnh trên, ta dùng hàm rsm để xử lý bảng dũ liệu trly và lưu kết quả xử lý vào biến kq. Sau đó dùng lệnh summary để xem kết quả xử lý. Trong hàm rsm, ta dùng phương pháp xử lý SO (second order) để thu được phương trình bậc 2. Kết quả thu được như sau:

> kq <- rsm(ThuHoi ~ SO(NhietDo, pH, DN), data = trly)
> summary(kq)

Call:
rsm(formula = ThuHoi ~ SO(NhietDo, pH, DN), data = trly)

               Estimate  Std. Error t value  Pr(>|t|)    
(Intercept) -1.6634e+03  3.5925e+02 -4.6304 0.0009356 ***
NhietDo      1.6319e+01  4.1859e+00  3.8987 0.0029673 ** 
pH           3.8812e+02  1.0050e+02  3.8618 0.0031508 ** 
DN           3.9349e+00  7.2739e+00  0.5410 0.6003769    
NhietDo:pH  -5.2500e-01  5.2255e-01 -1.0047 0.3387396    
NhietDo:DN   4.2500e-02  5.2255e-02  0.8133 0.4349736    
pH:DN       -4.5000e-01  1.0451e+00 -0.4306 0.6759068    
NhietDo^2   -8.7332e-02  2.0954e-02 -4.1679 0.0019245 ** 
pH^2        -3.0245e+01  8.3814e+00 -3.6086 0.0047792 ** 
DN^2        -2.0870e-01  8.3814e-02 -2.4901 0.0319869 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Multiple R-squared:  0.9357,	Adjusted R-squared:  0.8778 
F-statistic: 16.16 on 9 and 10 DF,  p-value: 7.769e-05

Analysis of Variance Table

Response: ThuHoi
                     Df Sum Sq Mean Sq F value    Pr(>F)
FO(NhietDo, pH, DN)   3 6016.4 2005.45 36.7217 1.042e-05
TWI(NhietDo, pH, DN)  3  101.4   33.79  0.6188  0.618580
PQ(NhietDo, pH, DN)   3 1825.1  608.37 11.1398  0.001577
Residuals            10  546.1   54.61                  
Lack of fit           5  494.6   98.92  9.6043  0.013364
Pure error            5   51.5   10.30                  

Stationary point of response surface:
  NhietDo        pH        DN 
79.247142  5.643505 11.411664 

Eigenanalysis:
$values
[1]  -0.08078122  -0.21128688 -30.24927489

$vectors
               [,1]         [,2]         [,3]
NhietDo  0.98347062  0.180858768 -0.008697175
pH      -0.00990701  0.005787113 -0.999934178
DN       0.18079653 -0.983492052 -0.007483226

Ta có thể điểm qua một số nét chính của kết quả này như sau:

Phần đầu của kết quả là bảng cho ta hệ số của phương trình hồi quy (cột Estimate) và giá trị p của các hệ số ấy (cột Pr(>|t|). Qua đó ta có thể kết luận rằng phương trình hồi quy có bậc hai, giữa các yếu tố không có tương tác, ...
Phương trình hồi quy có `R^2=0,9357` và `R_(hc)^2=0,8778` tương đối cao. Ngoài ra sự khác biệt giữa hai đại lượng này không quá lớn. Vì thế ta có thể kết luận phương trình hồi quy tương thích tốt với số liệu thí nghiệm.
Phân tích phương sai cũng được thực hiện cho các phương pháp xử lý FO (first order), TWI (two way interaction) và PQ (pure quadratic).
Cực trị của tỷ lệ thu hồi đạt được với nhiệt độ 79,247142°C; pH 5,643505; tỷ lệ dung môi / nguyên liệu 14,411664.

Phương pháp thiết kế phối hợp có tâm

Mô tả chung

Một số tính chất của thí nghiệm phối hợp có tâm

Thực hiện nhiều giai đoạn

Thí dụ

Thiết kế thí nghiệm

Xử lý số liệu

Khối	Đơn vị thí nghiệm	Nhiệt độ (°C)	pH	Tỷ lệ DN	Tỷ lệ thu hồi (%)
1	1	60	5,0	10	45
1	2	80	5,0	10	93
1	3	60	6,0	10	65
1	4	80	6,0	10	95
1	5	60	5,0	20	34
1	6	80	5,0	20	83
1	7	60	6,0	20	42
1	8	80	6,0	20	88
1	9	70	5,5	15	90
1	10	70	5,5	15	85
1	11	70	5,5	15	84
1	12	70	5,5	15	90
2	1	54	5,5	15	42
2	2	86	5,5	15	85
2	3	70	4,7	15	53
2	4	70	6,3	15	80
2	5	70	5,5	7	94
2	6	70	5,5	23	51
2	7	70	5,5	15	90
2	8	70	5,5	15	92

Khối	Đơn vị thí nghiệm	Nhiệt độ (°C)	pH	Tỷ lệ DN	Tỷ lệ thu hồi (%)
1	1	60	5,0	10	45
1	2	80	5,0	10	93
1	3	60	6,0	10	65
1	4	80	6,0	10	95
1	5	60	5,0	20	34
1	6	80	5,0	20	83
1	7	60	6,0	20	42
1	8	80	6,0	20	88
1	9	70	5,5	15	90
1	10	70	5,5	15	85
1	11	70	5,5	15	84
1	12	70	5,5	15	90
2	1	54	5,5	15	42
2	2	86	5,5	15	85
2	3	70	4,7	15	53
2	4	70	6,3	15	80
2	5	70	5,5	7	94
2	6	70	5,5	23	51
2	7	70	5,5	15	90
2	8	70	5,5	15	92

Khối	Đơn vị thí nghiệm	Nhiệt độ (°C)	pH	Tỷ lệ DN	Tỷ lệ thu hồi (%)
1	1	60	5,0	10	45
1	2	80	5,0	10	93
1	3	60	6,0	10	65
1	4	80	6,0	10	95
1	5	60	5,0	20	34
1	6	80	5,0	20	83
1	7	60	6,0	20	42
1	8	80	6,0	20	88
1	9	70	5,5	15	90
1	10	70	5,5	15	85
1	11	70	5,5	15	84
1	12	70	5,5	15	90
2	1	54	5,5	15	42
2	2	86	5,5	15	85
2	3	70	4,7	15	53
2	4	70	6,3	15	80
2	5	70	5,5	7	94
2	6	70	5,5	23	51
2	7	70	5,5	15	90
2	8	70	5,5	15	92