Chương 6: MÔ HÌNH HÓA VÀ TỐI ƯU HÓA THỰC NGHIỆM 6.1. Các phương pháp qui hoạch thực nghiệm Trong công nghệ hóa học, chúng ta nghiên cứu một đối tượng công nghệ thường phụ thuộc đồng thời vào nhiều yếu tố mà bản chất qui luật của quá trình xảy ra bên trong đối tượng chưa được biết rõ. Dựa vào những hiểu biết ban đầu về đối tượng, trước khi tiến hành thực nghiệm chúng ta cần xác định sơ bộ mô hình toán học của đối tượng, cần giải thích những yếu tố nào phải thay đổi trong quá trình làm thí nghiệm, những yếu tố nào giữ ở mức cố định và mục tiêu cần đạt được tối ưu. Phương pháp mô hình hóa toán học là phương pháp tính toán và phân tích các quá trình kỹ thuật, là vấn đề chọn công thức thực nghiệm và ước lượng các tham số của công thức đó. Để ước lượng giá trị thực của các đại lượng được khảo sát và độ chính xác của các ước lượng cần phải tiến hành xử lý các số liệu thực nghiệm. Phương pháp xử lý số liệu được sử dụng là phương pháp phân tích hồi qui. Đối tượng trong công nghệ hóa học phụ thuộc vào các yếu tố công nghệ có thể điều chỉnh như nồng độ, áp suất, nhiệt độ, độ pH,… và các yếu tố của biến ngẫu nhiên không điều khiển được. Hàm mục tiêu của đối tượng là hiệu suất, chất lượng sản phẩm, chi phí sản xuất,… Để tìm mối quan hệ giữa hàm mục tiêu y và các yếu tố đầu vào xi bằng thực nghiệm, chúng ta cho thay đổi tác động các biến đầu vào xi (đại lượng có thể đo và điều khiển được) và đo hàm đáp ứng đầu ra y (đại lượng đo được nhưng không điều khiển được) theo mô hình thực nghiệm như hình 6.1. x
X1 X2
y
Xi
xi – biến đầu vào (1 ≤ i ≤ k). y – thông số đầu ra (biến bị điều khiển). x - biến ngẫu nhiên không điều khiển được. Hình 6.1. Mô hình nghiên cứu thực nghiệm các đối tượng công nghệ. Ta cần thiết lập quan hệ: y = f(x1, x2,…,xk) + x Hay: y = f(X) + x
(6.1) (6.2)
Khi giả thiết biến ngẫu nhiên có phương sai D(x) = 2 và có kỳ vọng toán học E(x) = 0 thì đại lượng ngẫu nhiên này tuân theo một luật phân bố nào đó (như phân bố chuẩn) dạng: x = N(0, 2)
(6.3)
1
Nếu loại trừ được ảnh hưởng của nhiễu (x = 0) thì ta thu được các dạng mô hình thống kê như sau: - Nếu mô hình thống kê với các biến đầu vào chứa yếu tố thời gian thì gọi là mô hình động, không chứa yếu tố thời gian là mô hình tĩnh. - Nếu mô hình chỉ áp dụng cho một miền giới hạn nhất định của biến (điều kiện ràng buộc của biến) gọi là mô hình cục bộ (địa phương). Với mô hình địa phương, triển khai hàm f(x1, x2,…,xk) dưới dạng chuỗi Taylor: ∑
∑
∑
Với 1 ≤ i ≤ j ≤ k Các hệ số β0, βj, βjj,… được xác định từ số liệu thực nghiệm thì phương trình thu được gọi là phương trình hồi qui thực nghiệm của hệ thống, khi đó: ̂
∑
Phương trình hồi qui thực nghiệm (6.5) phụ thuộc vào bộ n thí nghiệm và phương pháp xứ lý số liệu thực nghiệm. Quá trình tiến hành thí nghiệm sao cho số lần thí nghiệm ít nhất, tính toán đơn giản và thu được kết quả chính xác nhưng vẫn khảo sát được ảnh hưởng đồng thời của nhiều thông số tác động đến đối tượng công nghệ để phản ánh chính xác bản chất của quá trình và giảm chi phí tiến hành thí nghiệm. Do đó, cần tiến hành thí nghiệm theo kế hoạch định trước sao cho có tính trực giao, nghĩa là bố trí thực nghiệm theo ma trận biến đầu vào có dạng: (ma trận n dòng, (k + 1) cột) [
Ma trận cột có dạng:
]
(6.6)
[ ]
(6.7)
[
Ma trận các hệ số hồi qui tuyến tính có dạng:
]
(6.8)
Ma trận chuyển vị của ma trận X có dạng: (ma trận (k + 1) dòng, n cột)
2
[
]
(6.9)
Nếu các kết quả thực nghiệm được biểu diễn theo phương trình hồi qui tuyến tính bằng phương pháp bình phương cực tiểu thì ta có dạng ma trận của hệ phương trình chuẩn: (
)
Suy ra:
(6.10) (
)
(6.11)
Với (XTX)-1 là ma trận nghịch của ma trận XTX khi định thức của ma trận XTX khác không hay ma trận XTX không suy biến. Ma trận trực giao X có những tính chất sau: - Tính trực giao: tích vô hướng của hai vectơ cột bất kỳ của X bằng 0. n
với j, m= ̅̅̅̅̅ 0, k
∑ xim xij =0
(6.12)
i=1
- Tính chất đối xứng: tổng các phần tử trong một cột bất kỳ đều bằng 0. n
∑ xij =0
với j
0.
(6.13)
i=1
Áp dụng tính chất của ma trân trực giao, chúng ta có ma trận XTX trở thành ma trận đường chéo:
∑
∑
[
]
Khi đó, ma trận nghịch đảo (XTX)-1 trong qui hoạch trực giao có dạng:
[
] 3
Với: n
∑
với j= ̅̅̅̅̅ 0, k
(6.16)
i=1
Ta tính ma trận XTY: ∑
[
][
]
∑
∑ [
]
Khi đó ma trận các hệ số hồi qui được tính theo công thức (6.11). Trường hợp tổng quát, nếu các số liệu thực nghiệm được biểu diễn bằng đa thức bậc α thì số hệ số hồi qui trong phương trình hồi qui bằng
. Sử dụng một số phép biến đổi,
chúng ta biến đổi phương trình hồi qui đa thức bậc α về phương trình hồi qui tuyến tính bằng cách thay các số hạng phi tuyến của thức bằng các số hạng tuyến tính. Phương trình này được gọi là phương trình tuyến tính hóa theo thông số. Nếu bậc của đa thức chưa biết trước thì việc tính toán thường phải tiến hành vài lần, tăng dần bậc của đa thức đến khi nào phương trình hồi qui nhận được tương thích với thực nghiệm. Mỗi lần tăng bậc đa thức, toàn bộ việc tính toán liên quan đến phân tích hồi qui đều phải tính toán lại. Sự thay đổi bậc của đa thức hoặc loại bỏ một phần các số hạng trong phương trình hồi qui dẫn đến thay đổi các giá trị tính toán các hệ số còn lại trong đa thức. Khi đó xuất hiện sự không xác định trong việc ước lượng các hệ số hồi qui, gây khó khăn cho việc giải thích ảnh hưởng của các biến độc lập vào đại lượng nghiên cứu. Để loại bỏ các khó khăn nói trên, chúng ta tiến hành thực nghiệm theo qui hoạch trực giao và ma trận thực nghiệm trực giao phải thỏa mãn các điều kiện sau: - Công thức tính các hệ số phương trình hồi qui bj. - Hệ số bj là ước lượng trúng của các hệ số βj. - Phương trình hồi qui ̂ là ước lượng trúng của y. 6.1.1. Qui hoạch thực nghiệm yếu tố toàn phần 2k Qui hoạch thực nghiệm yếu toàn phần là thực nghiệm mà mọi tổ hợp các mức của các yếu tố đều được thực hiện để nghiên cứu. Bố trí thực nghiệm yếu tố toàn phần theo ma trận trực giao và phân bố đối xứng các biến độc lập với tâm đối xứng. Ma trận X thỏa mãn điều kiện trực giao và có thêm tính chất chuẩn hóa: Tổng bình phương các phần tử của một cột bằng số thí nghiệm. 4
n
∑
với j= ̅̅̅̅̅ 0, k
=N
(6.18)
i=1
Với k là các yếu tố, n là số mức thì số thí nghiệm N = nk . Nếu các thí nghiệm chỉ thực hiện ở hai mức, thường là hai giá trị biến của mỗi yếu tố khảo sát thì N = 2k. Giả sử khảo sát biến Zj với hai giá trị biên là aj < bj (1 ≤ j ≤ k), chúng ta thực hiện một số phép biến đổi sau: Đặt: Khi đó ta đổi biến số theo biểu thức: Suy ra: Zj [aj, bj] thì xj [-1, 1]. Khi đó ta áp dụng công thức (6.18) vào công thức (6.14) ta có: [
] = N.I
Suy ra: [
]
Ma trận hệ số hồi qui được tính theo công thức (6.11): ∑
(
∑
)
[
∑
]
Hay: ∑
∑
Khi đó phương trình hồi qui có dạng: ̂
∑
Trong trường hợp phương trình hồi qui được mô tả dưới dạng: ̂
∑
∑ 5
Các hệ số bij được xác định theo công thức sau: ∑ Với số lượng hệ số bij được xác định theo công thức:
6.1.2. Qui hoạch thực nghiệm yếu tố từng phần 2k - p Khi số biến k trong mô hình qui hoạch thực nghiệm yếu tố toàn phần 2k càng lớn thì làm cho số thí nghiệm N càng lớn. Điều này làm cho qui hoạch trở nên cồng kềnh, chi phí lớn, kém hiệu quả. Để khắc phục điều đó, người ta dùng qui hoạch thực nghiệm yếu tố từng phần N = 2k – p, với p là giá trị đặc trưng cho độ từng phần, là số hiệu ứng tương tác được thay bằng số hiệu ứng tuyến tính. Thực chất đây là qui hoạch thực nghiệm yếu toàn phần bớt đi p cột của k thông số độc lập, số thí nghiệm giảm đi 2p lần nhưng vẫn đảm bảo tính trực giao của ma trận X. Quá trình qui hoạch thực nghiệm yếu tố từng phần theo các bước sau: Bước 1: chọn ra r thông số chính ảnh hưởng đến hàm mục tiêu trong k thông số đầu vào: r = k – p. Lập qui hoạch thực nghiệm yếu toàn phần 2r với số thí nghiệm là N = 2r. Tuy nhiên, khi lựa chọn giá trị p phải đảm bảo điều kiện sau: k + 1 ≤ N = 2r = 2k-p ≤ 2k
(6.24)
Bước 2: tiến hành thiết lập các biểu thức tương quan sinh biểu diễn các mối tương quan giữa mỗi thông số p với một tích các thông số trong r thông số chính. Các biểu thức tương quan sinh có thể là tích của các thông số trong r thông số chính mang dấu dương hay âm. Bước 3: kiểm tra tính tiện lợi của mô hình đã được lập: nếu ma trận X không có các cột giống nhau hoặc ngược dấu nhau thì vẫn đảm bảo tính trực giao, qui hoạch thực nghiệm đạt yêu cầu. Bước 4: tiến hành xác định và kiểm tra ý nghĩa của các hệ số hồi qui bj, kiểm tra sự tương thích của phương trình hồi qui thu được. 6.1.3. Qui hoạch trực giao cấp 2 Qui hoạch trực giao cấp 2 là qui hoạch thực nghiệm và xử lý số liệu thực nghiệm theo phương pháp xây dựng mô hình hồi qui cấp 2 với các điều kiện tương tự như qui hoạch thực nghiệm yếu tố toàn phần (qui hoạch trực giao cấp 1). Phương trình hồi qui bậc 2 đầy đủ có dạng: ̂
∑
∑
∑ 6
Xây dựng ma trận trực giao X bao gồm ba loại thí nghiệm: - Phần cơ sở gồm n = 2k thí nghiệm theo qui hoạch thực nghiệm yếu tố toàn phần. - Phần điểm “*” gồm nk = 2k điểm nằm trên các trục tọa độ của không gian k yếu tố và cách tâm phương án khoảng cách α > 0. - Phần tâm gồm n0 (n0 ≥ 1) thí nghiệm ở tâm phương án dùng để xác định phương sai tái hiện trong công thức kiểm tra ý nghĩa của các hệ số hồi qui. Tổng số thí nghiệm trong phương án là N = 2k + 2k + n0 Tuy nhiên, những phương án cấu trúc có tâm không trực giao do
> 0 nên:
N
∑
với j= ̅̅̅̅̅ 1, k
0
(6.26)
i=1 N
∑
với j= ̅̅̅̅̅ 1, k , j
0
u
(6.27)
i=1
Vì vậy, khi xây dựng ma trận trực giao X cần chọn α sao cho trực giao hóa công thức (6.26), (6.27). Giả sử xét qui hoạch thực nghiệm k = 2 yếu tố và n0 = 1 thì ma trận X có dạng:
(6.28)
[ Đặt
]
(i = 1, 2) thì ta có:
(6.29)
[
]
Sử dụng các điều kiện trực giao của ma trận X: Tích vô hướng của vectơ cột 1 và cột 6 bằng 0, tức là ta có đẳng thức sau: 22(1 - ) + 2α2 – 2.2
-
=0
Hay tổng quát với k yếu tố: 2k(1 - ) + 2α2 – 2k 7
- n0. = 0
Suy ra:
(6.30)
Tích vô hướng của vectơ cột 5 và cột 6 bằng 0, tức là ta có đẳng thức sau: 22(1 - )2 - 2.2 ( α2- ) +
2
=0
Hay tổng quát với k yếu tố: 2k(1 - )2 - 2k ( α2- ) + n0 .
2
=0
√√
Suy ra:
(6.31)
Thế công thức (6.30), (6.31) vào để tính được ma trận X. Khi đó các hệ số hồi qui được tình theo công thức (6.11):
(
)
6.1.4. Qui hoạch thực nghiệm để tìm cực trị Các phương trình hồi qui thu được là biểu diễn gần đúng đối tượng công nghệ mà ta đang khảo sát. Vấn đề đặt ra là tìm
,
,…,
sao cho: (6.32)
Khi hàm hồi qui đạt đủ độ chính xác cần thiết, ta dùng qui hoạch toán học để tìm cực trị: - Nếu hàm tuyến tính ta dùng phương pháp qui hoạch tuyến tính; - Nếu hàm phi tuyến ta áp dụng phương pháp qui hoạch phi tuyến. Nếu độ chính xác chưa đạt yêu cầu, kết quả còn thô, kém tin cậy thì trước khi sử dụng phương pháp qui hoạch toán học cần thu hẹp vùng chứa điểm cực trị, tức là có thể tìm vùng cực trị bẳng qui hoạch thực nghiệm. Phương pháp qui hoạch thục nghiệm để tìm cực trị được chia làm hai giai đoạn: - Tìm vùng chứa điểm cực trị bằng qui hoạch trực giao cấp 1. - Tìm phương trình hồi qui cấp 2 bằng qui hoạch trực giao cấp 2 và cuối cùng là dùng phương pháp qui hoạch phi tuyến để tìm cực trị: qui hoạch lồi hoặc qui hoạch toàn phương. Vectơ gradient của hàm y = f(x) tại x* (ký hiệu là grad[f(x*)]) là vectơ có chiều biểu thị sự biến thiên nhanh nhất của hàm y tại x*, giá trị của grad[f(x*)] thay đổi từ điểm này sang điểm khác trong không gian yếu tố. Với mô hình tuyến tính k yếu tố: [ Với:
+
]
(6.33)
: đạo hàm riêng của hàm f(x) theo biến xj tại x*.
+ ij (j = 1, 2, …,k) : vectơ đơn vị theo các trục tọa độ. Xét một miền con D0 có tâm Z0 ứng với x0, giả sử phương trình hồi qui có dạng:
8
̂
∑
Kiểm định sự tương hợp của ̂, nếu tương hợp thì có nghĩa là mặt cong được xấp xỉ bằng mặt phẳng thì D0 không chứa điểm cực trị. Chuyển sang vùng D1 theo hướng gradient [f(x0)] cho đến khi y không tăng được nữa. Lặp lại quá trình này đến khi hàm ̂ không tương hợp thì đã chuyển sang vùng chứa điểm cực trị, cần chuyển sang bước 2. Theo thuật toán leo Box – Wilson như sau: Tính các thành phần của gradient theo triển khai Taylor: ∑
∑∑
∑
Nếu ta lấy gần đúng đến số hạng bậc nhất và đặt: ∑
∑
Lượng mà hàm f(x) đã tăng là: ∑ Bước tiến nhanh nhất tương ứng với số hạng làm (6.37) tăng nhanh nhất j* là: |
|
j = ̅̅̅̅̅
(6.38)
Độ dài các bước hj của các yếu tố được tính theo j*:
Chú ý: độ dài các bước không nên quá lớn hoặc quá nhỏ. Chuyển động theo grad phải bắt đầu từ điểm 0. Mức cơ sở của mỗi yếu tố và ngừng lại nếu tìm được điểm tối ưu hoặc nếu những hạn chế đặt vào các yếu tố làm cho chuyển động tiếp tục theo grad không hợp lý nữa. Cách tiến hành: gọi Z0 là tâm miền D0 trong tọa độ biến thật ký hiệu là M0. Điểm M1 có ̅̅̅̅̅ . Tại đây ta làm thí nghiệm xác định được các tọa độ xác định theo: y2. Lặp lại ta sẽ thu được dãy: y0, y1, …, yn. Khi thực hiện: y0 < y1 tiếp tục làm y2. y2 < y3 tiếp tục làm y3. y3 < y4 tiếp tục làm y4. 9
yp-1 > yp mặt cong bắt đầu. Dừng lại ở Mp-1 với tâm mới Mp-1 = Zp-1. Tại Xp-1 lại tiếp tục quá trình xấp xỉ bằng mặt phẳng. Kiểm định sự phù hợp của mô hình bậc nhất, nếu thỏa mãn thì tiếp tục tìm độ dài bước làm thí nghiệm. Tiến hành cho đến khi mô hình bậc nhất không phù hợp thì đó là vùng cực trị. 6.2. Xác định các hệ số của phương trình hồi qui Chọn phương án thí nghiệm với k biến độc lập, thực hiện n thí nghiệm, không có thí nghiệm lập lại và bố trí thực nghiệm sao cho ma trận X có tính trực giao. 6.2.1. Tính hệ số hồi qui bj của phương trình hồi qui tuyến tính Áp dụng công thức (6.11), (6.15), (6.16), (6.17) ta có: ∑
(
∑
)
] ∑ [
[
]
Suy ra: ∑ ∑
[
∑
]
6.2.2. Tính hệ số hồi qui bj của phương trình hồi qui bậc 2 Tiến hành qui hoạch thực nghiệm theo mô hình qui hoạch trực giao cấp 2 thì các hệ số hồi qui được xác định theo công thức (6.11) sau khi lập được ma trận trực giao X, hay xác định theo các công thức sau: bj =
∑N i=1 xij y
i
2 ∑N i=1 xij
∑ ∑ ∑ ∑ 10
(6.42)
Khi đó phương sai
được xác định theo công thức: ∑
∑
∑ 6.3. Kiểm tra sự tương hợp của phương trình hồi qui Sau khi tính toán được các hệ số của phương trình hồi qui, chúng ta cần tiến hành kiểm định sự tương hợp của phương trình hồi qui với số liệu thực nghiệm. Quá trình kiểm định được tiến hành theo hai bước sau: Bước 1: kiểm tra sự có nghĩa của các hệ số phương trình hồi qui bằng tiêu chuẩn thống kê (chuẩn số) Student tα (với α: mức ý nghĩa). Bản chất của phương pháp là kiểm tra các hệ số bj = 0 hay không, hoặc kiểm định xem thực chất có bao nhiêu yếu tố ảnh hưởng đến hàm mục tiêu. Chọn thống kê:
| |
(6.46)
sbj: độ lệch quân phương của hệ số thứ i. Phươn sai
được xác định theo công thức: ∑
Với phương sai tái hiện
được tính theo số thí nghiệm lặp ở tâm n0 theo công thức: ∑
̅
Trong đó: - N: số thí nghiệm. - n0: số thí nghiệm lặp ở tâm. - fth = n0 – 1 : bậc tự do tái hiện. Thông thường chọn mức ý nghĩa α = 0,05. Tra theo chuẩn Student ta có được giá trị t α(fth). - Nếu tbi > tα(fth) thì hệ số bi được giữ lại trong phương trình hồi qui, ảnh hưởng của yếu tố xi có ý nghĩa đối với việc thay đổi thông số tối ưu y. - Ngược lại, nếu tbi < tα(fth) thì hệ số bi bị loại khỏi phương trình hồi qui. Bước 2: kiểm tra sự tương thích của phương trình hồi qui theo tiêu chuẩn Fisher:
Phươn sai
được xác định theo công thức: 11
∑
̅
Trong đó: - L: số hệ số có ý nghĩa trong phương trình hồi qui. Tra theo chuẩn Fisher ta có được giá trị Fα(α, ftt, fth), với ftt = N – L , fth = n0 – 1. - Nếu F < Fα thì mô hình thống kê phù hợp với số liệu thực nghiệm. - Nếu F ≥ Fα thì mô hình thống kê không phù hợp với số liệu thực nghiệm. 6.4. Ứng dụng Excel để tối ưu hóa thực nghiệm Để giải bài toán tối ưu hóa thực nghiệm chúng ta cần tiến hành các bước sau: - Bước 1: chọn phương án tiến hành thí nghiệm. - Bước 2: lập ma trận thực nghiệm X. - Bước 3: tiến hành thí nghiệm để xác định giá trị biến đầu ra Y. - Bước 4: xác định các hệ số trong phương trình hồi qui. - Bước 5: đánh giá phương trình hồi qui thu được. - Bước 6: xác định chế độ thực nghiệm tối ưu. Mỗi bước có thể thực hiện phép lặp nếu chưa đạt yêu cầu về mục tiêu của công nghệ, sai số…Các bước 1, 2 và 4, 5 đã được trình bày ở phần trên về phương pháp thực hiện. Bước 3 được thực hiện ở phòng thí nghiệm, phụ thuộc vào trang thiết bị, dụng cụ thí nghiệm và tay nghề của người làm thí nghiệm. Bước 6 được tối ưu hóa theo mô hình tìm ra. Khảo sát phương pháp thực hiện bước 4 và 5 bằng phần mềm Microsoft Excel để sử dụng trong tính toán bằng các hàm: - Nhân hai ma trận: MMULT(array1, array2). - Tính định thức của ma trận: MDETERM(array). - Tính ma trận nghịch đảo: MINVERSE(array). - Tính ma trận chuyển vị: TRANPOSE(array). - Tính giá trị trung bình các số hạng: AVERAGE(number1, number2,…). - Tính tổng bình phương các số hạng ∑ - Tính tổng bình phương độ lệch ∑ - Tính độ lệch chuẩn của mẫu √
∑
: SUMSQ(number1, number2,…). ̅ : SUMXMY2(array_x, array_y).
̅
: STDEV(number1, number2,…).
- Tra chuẩn số Student: TINV(p1, p2). 12
- Tra chuẩn số Fisher: FINV(α, p1, p2). Lưu ý: để thực hiện các phép toán ma trận trong Excel, chúng ta phải chọn các phần tử của ma trận cần tính (quét khối kích thước ma trận) và ấn giữ đồng thời ba phím Ctrl + Shift + Enter. Ví dụ 1: Tiến hành qui hoạch thực nghiệm nghiên cứu ảnh hưởng của 3 yếu tố Z1, Z2, Z3 lên hàm mục tiêu y với số liệu thu được như sau: Biến thực
n=2
k
n0 = 3
1
2
3
4
5
6
7
8
9
10
11
Z1
150
300
150
300
150
300
150
300
225
225
225
Z2
30
30
90
90
30
30
90
90
60
60
60
Z3
15
15
15
15
45
45
45
45
30
30
30
y
3,0
6,0
10,0
12,0
15,0
23,0
12,0
18,0
12,0
13,8
13,2
Hãy tìm mối quan hệ giữa y và các biến Z1, Z2, Z3 theo mô hình trực giao cấp 1? Giải: Chọn phương trình hồi qui biểu diễn mối quan hệ giữa y và các biến Z 1, Z2, Z3 có dạng theo (6.21) thì qui trình đề nghị giải theo các bước sau: Bước 1: lập 1 bảng tính excel với các thông tin sau:
Trong đó: - Giá trị mức thấp, mức cao của các biến thực lấy từ số liệu đã cho. - Zj0 : mức cơ sở, là trung bình cộng của mức thấp và mức cao nên ta nhập công thức ở cell D2 = AVERAGE(B2:C2)
Sau đó dùng draf kéo chuột xuống các hạng trong cột Zj0. 13
Khoảng biến thiên Zj được tính theo công thức trong cell E2 = (C2-B2)/2
Sau đó dùng draf kéo chuột xuống các hạng trong cột Zj. Bước 2: lập bảng tính chuyển đổi các biến thực sang biến mã hóa:
Nhập giá trị các ô trong cột x0 đều bằng 1. nên ở cell G8 ta nhập: =(C8-$D$2)/$E$2, sau đó dùng draf kéo chuột xuống các hạng trong cột x1.
nên ở cell H8 ta nhập: =(D8-$D$3)/$E$3, sau đó dùng draf kéo chuột xuống các hạng trong cột x2. 14
nên ở cell I8 ta nhập: =(E8-$D$4)/$E$4, sau đó dùng draf kéo chuột xuống các hạng trong cột x3.
Bước 3: xác định các hệ số phương trình hồi qui bằng phương pháp ma trận: (
)
, ta lập bảng excel như sau:
Ma trận hệ số X: do có 3 thí nghiệm được thực hiện ở tâm nên ma trận X là ma trận 8 dòng, 4 cột. Để xác định giá trị các phần tử trong ma trận X ta copy giá trị các ô F8:I15 vào các ô A21:D28. Ma trận các biến đầu ra y: tương tự như ma trận X, ta copy giá trị các ô J8:J15 vào các ô F21:F28. Ma trận chuyển vị XT: là ma trận có kích thước 4 dòng, 8 cột. Để xác định giá trị các phần tử trong ma trận XT ta làm như sau: sau khi quét khối các ô H21:O24, ta nhập công thức: =TRANSPOSE(A21:D28) và nhấn đồng thời ba phím Ctrl + Shift + Enter.
15
Ma trận XTX: là ma trận có kích thước 4 dòng, 4 cột. Nhân hai ma trận XT với ma trận X: ta quét khối các ô H27:K30 và nhập công thức: =MMULT(H21:O24,A21:D28) và nhấn đồng thời ba phím Ctrl + Shift + Enter.
Kiểm tra định thức det (XTX): ta nhập công thức =MDETERM(H27:K30). Ma
trận
nghịch
đảo
(XTX)-1:
quét
khối
các
ô
M27:P30,
nhập
công
thức:
=MINVERSE(H27:K30) và nhấn đồng thời ba phím Ctrl + Shift + Enter.
Ma trận XTY: là ma trận có kích thước 4 dòng, 1 cột. Ta nhập công thức: =MMULT(H21:O24,F21:F28). Ma trận hệ số hồi qui B: ta quét khối các ô D31:D34, nhập công thức = MMULT(M27:P30,A31:A34) nhấn đồng thời ba phím Ctrl + Shift + Enter.
16
Ta thu được các hệ số hồi qui: b0= 12,375; b1= 2,375; b2= 0,625; b3= 4,625. Bước 4: tiến hành lập bảng tính để kiểm tra ý nghĩa của các hệ số hồi qui:
Phương sai tái hiện
: được tính theo số thí nghiệm lặp ở tâm n0=3, theo công thức (6.48)
nên tại ô O8 ta nhập công thức: =STDEV(J16:J18)^2. Độ lệch quân phương Sbj: được tính theo công thức (6.47) nên tại ô P8 ta nhập công thức: =SQRT(O8/B15). Các hệ số tbj: tính theo công thức
| |
nên ta nhập công thức tại ô Q8: =Abs(N8)/$P$8
và sau đó dùng draf kéo chuột xuống các hạng trong cột tbj.
Chuẩn số student tα(fth): chọn mức ý nghĩa α = 0,05, bậc tự do tái hiện f th = n0 – 1 = 2 nên ta nhập công thức tại ô O12: =TINV(0.05,2). Kết quả chọn hệ số hồi qui: tại ô R8 ta nhập công thức: =IF(Q8<$O$12,”loại”,”chọn”) và kéo draf chuột xuống đến hết cột.
17
Kết quả ta loại bỏ hệ số b2 nên phương trình hồi qui có dạng: ̂ Bước 5: lập bảng tính kiểm tra sự tương hợp của phương trình hồi qui thu được: Tính giá trị của cột ̂ : sử dụng phương trình hồi qui thu được để tính các giá trị ̂, tại ô K8 nhập công thức: =12.375 + 2.375*G8+4.625*I8 . Sau đó kéo draf chuột xuống đến hết cột.
Phương sai
: xác định theo công thức (6.50) nên ta nhập công thức tính tại ô I1 như sau:
=SUMX2MY2(J8:J15,K8:K15)/(B15-3), với L = 3.
Tính hệ số
: ta nhập công thức tại ô I2: =I1/O8. 18
Tra chuẩn số Fisher Fα(α, ftt, fth): với ftt = N – L = 5 , fth = n0 – 1 = 2 nên ta nhập công thức tại ô I3: =FINV(0.05,5,2). Kiểm tra tính tương thích: tại ô I4 ta nhập công thức: =IF(I2
Kết quả so sánh ta thấy phương trình hồi qui thu được tương thích với thực nghiệm. Ví dụ 2: Tìm mối quan hệ hồi qui theo qui hoạch cấp 2 theo 3 biến với số liệu thực nghiệm như sau: N
1
2
3
4
5
6
7
8
9
y
13,9
18,5
2,0
3,0
16,0
18,5
9,0
12,0
15,0
N
10
11
12
13
14
15
16
17
18
y
8,0
7,5
15,8
11,5
5,0
10,1
11,2
9,9
12,3
Giải: Chọn phương trình hồi qui bậc 2 đầy đủ theo công thức (6.25) thì qui trình đề nghị giải theo các bước sau: Bước 1: Xác định các thông số α, các biến phụ: Với mô hình trực giao cấp 2 với k = 3 yếu tố thì tổng số thí nghiệm phải thực hiện là N = 2k + 2k + n0 = 23 + 2 . 3 + n0 = 18, nên số thí nghiệm lặp lại ở tâm n0 = 4. Thông số α được tính theo công thức (6.31): √√
= √√
√
Các biến phụ được tính theo công thức (6.30): = Bước 2: lập bảng tính ma trận thực nghiệm X với các biến đã mã hóa: Các giá trị cột x0, x1, x2, x3: được mã hóa từ các giá trị biên của 3 biến thực và phân bố thí nghiệm ở các điểm “*” và ở tâm phương án. 19
Giá trị cột x1x2, x1x3, x2x3: là tích của hai cột x1 và x2, x2 và x3, x2 và x3. Ta lần lượt nhập công thức ở ô G3=D3*E3, tại ô H3=D3*F3, ô I3=E3*F3. Sau đó kéo draf chuột xuống đến hết từng cột. Giá trị cột x21 – 2/3, x22 – 2/3, x23 – 2/3: lần lượt nhập công thức cho các cột tương ứng như sau: ô J3=D3*D3-2/3, ô K3=E3*E3-2/3, ô L3=F3*F3-2/3 và lần lượt kéo draf chuột đến hết từng cột. Kết quả ta có bảng tính sau:
Bước 3: xác định các hệ số phương trình hồi qui bằng phương pháp ma trận: (
)
Ma trận X: có kích thước là 14 dòng, 10 cột, giá trị các phần tử của ma trận X là giá trị của các ô C3:L16. Từ ma trận X, ta xác định được các ma trận XT, XTX, (XTX)-1, XTY và các ma trận này có kích thước như sau:
n dòng x m cột
XT
XTX
(XTX)-1
Y
XTY
B
10 x 14
10 x 10
10 x 10
14 x 1
10 x 1
10 x 1
Ngoài phương pháp tính ma trận B theo các ma trận trên, ta có thể tính trực tiếp ma trận B trực tiếp thông qua ma trận X như công thức: ta quét khối các ô C23:L23, nhập công thức: =TRANSPOSE(MMULT(MINVERSE(MMULT(
TRANSPOSE
(C3:L16),
C3:L16)),MMULT(TRANSPOSE(C3:L16),M3:M16))) và nhấn đồng thời ba phím Ctrl + Shift + Enter.
20
Kết quả ta thu được ma trận B như sau:
Bước 4: tiến hành lập bảng tính để kiểm tra ý nghĩa của các hệ số hồi qui:
Phương sai tái hiện
: được tính theo số thí nghiệm lặp ở tâm n0=4, theo công thức (6.48)
nên tại ô N22 ta nhập công thức: =STDEV(M17:M20)^2. Phương sai các hệ số S2bj: được tính theo công thức (6.45) nên tại ô C24 ta nhập công thức: =$N$22/SUMSQ(C3:C20). Sau đó kéo draf chuột đến hết dòng S2bj.
Các hệ số tbj: tính theo công thức
| |
nên ta nhập công thức tại ô C25: =
ABS(C23)/SQRT(C24) và sau đó dùng draf kéo chuột đến hết dòng tbj.
Chuẩn số student tα(fth): chọn mức ý nghĩa α = 0,05, bậc tự do tái hiện f th = n0 – 1 = 3 nên ta nhập công thức tại ô N23: =TINV(0.05,3). 21
Kết quả chọn hệ số hồi qui: tại ô C26 ta nhập công thức: =IF(C25<$N$23,"Loại","Chọn") và kéo draf chuột đến hết dòng.
Kết quả ta loại bỏ hệ số b12, b13, b33 nên phương trình hồi qui có dạng: ̂ = 10,4667 + 1,75 x1 – 4,3865 x2 + 2,2744 x3 + 1,7375 x2x3 + 1,6625 (x21 – ) + 1,7375 (x22 – ) Bước 5: lập bảng tính kiểm tra sự tương hợp của phương trình hồi qui thu được: Tính giá trị của cột ̂ : sử dụng phương trình hồi qui thu được để tính các giá trị ̂, tại ô N3 nhập
công
thức:
=$C$23*C3
+$D$23*D3+$E$23*E3+$F$23*F3+$I$23*I3+$J$23*J3+$K$23*K3. Sau đó kéo draf chuột xuống đến hết cột.
Phương sai
: xác định theo công thức (6.50) nên ta nhập công thức tính tại ô N24 như
sau: =SUMX2MY2(M3:M16,N3:N16)/(B20-7), với L = 7.
Tính hệ số
: ta nhập công thức tại ô N25: =N24/N22.
Tra chuẩn số Fisher Fα(α, ftt, fth): với ftt = 18 – 7 = 11 , fth = n0 – 1 = 3 nên ta nhập công thức tại ô N26: =FINV(0.05,11,3). Kiểm tra tính tương thích: tại ô O23 ta nhập công thức: =IF(N25
22
Kết quả so sánh ta thấy phương trình hồi qui thu được tương thích với thực nghiệm. 6.5. Bài tập: Bài tập 1: Lập phương trình hồi qui mô tả ảnh hưởng của: hàm lượng kim loại trong tổng đioxit Zr(Hf)O2, g/l; nồng độ acid , đlg/ l; nồng độ tributyl photphat trong O – xylen, phần (V); tỷ số thể tích 2 pha lên quá trình tách hafini va ziricon bằng tributyl photphat theo số liệu thực nghiệm sau (một thí nghiệm ở tâm có y0 = 14,1): n = 2k = 24 = 16
N
1
2
3
4
5
6
7
8
9
10
11
12
y
11,5
8,2
8,5
14,0
12,7
10,0
12,5
13,5
10,58
9,48
14,1
9,65
n = 2k = 16
2k = 2.4 = 8
N
13
14
15
16
17
18
19
20
21
22
23
24
y
11,38
11,17
11,8
13,3
13,8
12,95
7,9
13,05
12,9
10,3
12,5
13,9
Bài tập 2: Khi khảo sát hiệu suất tạo sản phẩm phụ thuộc vào nồng độ đầu các chất (mol/l) CA0 , CB0, CC0 và nhiệt độ T (0C) thu được bảng số liệu sau: Nồng độ (mol/l).10 N
CA
CB
CC
2
T =80 ph
=160 ph
=320 ph
=640 ph
=1280 ph
1
-
-
-
-
3,17
5,39
8,66
15,90
22,60
2
+
-
-
-
14,8
23,4
34,30
34,60
20,30
3
-
+
-
-
4,80
10,80
22,50
34,60
42,00
4
+
+
-
-
23,20
39,00
55,60
63,40
41,60
23
5
-
-
+
-
3,72
3,81
17,20
20,00
23,90
6
+
-
+
-
17,90
28,30
40,50
34,20
21,60
7
-
+
+
-
8,60
13,30
25,90
39,80
50,80
8
+
+
+
-
30,90
51,40
72,20
76,40
38,90
9
-
-
-
+
7,48
9,93
20,00
30,90
24,90
10
+
-
-
+
25,30
35,30
39,10
28,40
7,50
11
-
+
-
+
13,50
27,10
43,00
58,00
49,40
12
+
+
-
+
50,80
75,60
84,20
57,00
11,50
13
-
-
+
+
9,15
15,80
27,50
33,90
23,00
14
+
-
+
+
30,80
44,40
46,70
24,90
2,94
15
-
+
+
+
22,80
37,20
57,90
69,10
53,90
16
+
+
+
+
62,60
88,00
89,50
43,40
5,80
Các phản ứng xảy ra theo cơ chế sau: A+BF A + F G, F là sản phẩm chính. Hãy xác định hằng số tốc độ các phản ứng trên? Bài tập 3: Xác lập phương trình hồi qui khi khảo sát quá trình biến tính nhôm bằng molipden (Mo) phụ thuộc vào : lượng Mo đưa vào (%); nhiệt độ nung ( 0C); thời gian nung (ph) và tốc độ nguội có tính định tính (nhanh – grafit; chậm – samot). Hàm mục tiêu là số hạt nhôm trên 1cm2. Các yếu tố ảnh hưởng: Yếu tố
Z1
Z2
Z3
Z4
Cận trên
0,55
940
120
Grafit
Cận dưới
0,25
740
0
Samot
Số liệu thực nghiệm thu được: 24
N
1
2
3
4
5
6
7
8
9
10
11
y
64
90
69
130
36
95
81
100
80
82
78
25
