Đề tài Xây dựng mô hình động học lò hơi bằng phương pháp nhận dạng vòng kín sử dụng matlap

ố thấp, nên ước lượng đáp ứng trên miền tần số cao không đảm bảo chính xác. Nhận xét: Với phương pháp phân tích phổ, để thu được mô hình đáp ứng tần số tốt ta cần sử dụng tín hiệu vào u có phổ tần trải rộng trên miền tần số quan tâm. Chính vì lý do trên phương pháp này thích hợp cho nhận dạng chủ động, vòng hở vì khả năng linh động trong chọn tín hiệu vào hệ thống. 2.2.1.3. Ước lượng hàm truyền đạt liên tục từ đáp ứng tần số Mô hình đáp ứng tần số cho ta hình ảnh trực quan về đáp ứng của hệ thống trên miền tần số, nhưng lại chưa đưa ra được hàm truyền đạt cụ thể để có thể xây dựng được đáp ứng trên miền thời gian của hệ thống. Từ mô hình tần số ta có được đồ thị bode và có thể dựa vào đó có thể xác định được hệ số khuếch đại tĩnh và bậc của mô hình. Hệ số khuếc đại tĩnh K là giá trị biên độ của đồ thị bode tại ω = 0; bậc của mô hình xác định từ độ dốc của đồ thị khi ω → ∞. Như vậy ta đã có thông tin ban đầu về bậc và hệ số khuếch đại tĩnh của hệ thống qua đó sẽ xây dựng mô hình hàm truyền với tham số là hằng số thời gian Ti. Sử dụng phương pháp bình phương sai lệch cực tiểu để xác định Ti. Thực tế trong miền tần số lớn mô hình đáp ứng tần số ta thu được từ ước lượng trên bị ảnh hưởng của nhiễu và phép chuyển đồi Fourier gián đoạn. Chính vì vậy trong vùng tần số này đáp ứng tần số của hệ thống thu được là không trung thực không nên sử dụng để ước lượng mô hình hàm truyền. Dải tần số phản ảnh đặc tính của hệ thống là dải cho đáp ứng pha trong khoảng (0, 2π) nên ta chỉ sử dụng thông tin về đáp ứng tần số trong dải này. Hàm ước lượng từ đáp ứng tần số của hệ thống sang hàm truyền đạt bằng nguyên lý bình phương cực tiểu hàm sai lệch được xây dưng thành hàm invfreqs trong Matlab identification toolbox . Cấu trúc lệnh: invfreqs(h,w,n,m) Với: w: là dãy tần số sử dụng để ước lượng h: dãy đáp ứng của hệ thống tương ứng với mỗi tần số trên n: bậc của đa thức tử trong mô hình hàm truyền m: bậc của đa thức mẫu trong mô hình hàm truyền. Ví dụ 2.2: Sử dụng mô hình đáp ứng tần số h12 để ươc lượng ngược lại mô hình của hệ thống với đáp ứng tần số ước lượng được có dạng: (phụ lục 2.2) ≈ 40dB/dec Hình 2.14. Đồ thị bode của mô hình đáp ứng tần số h13 Từ đáp ứng tần số ước lượng được nhờ phương pháp spafdr từ hệ thống (*) ta có được những thông tin sau về hệ thống: 1. Hệ số khuếch đại của hệ thống ≈ 100.56 2. Độ dốc biên độ ở dải tần số cao ≈ 40dB, như vậy hàm truyền đạt của hệ thống (*) có bậc tử lớn hơn bậc mẫu ít nhất là 2. Ước lượng mô hình hàm truyền đạt liên tục từ đáp ứng tần số trên. Ta sử dụng hàm invtf (phụ lục 2.2) g1 = invtf(h13,0,1,80); g2 = invtf(h13,0,2,80); g3 = invtf(h13,1,2,80); g4 = invtf(h13,1,3,80); g5 = invtf(h13,1,4,80); bode(g,'b',g1,'g--',g2,'g',g3,'y',g4,'b--',g5,'r') Hình 2.15. Mô hình hàm truyền ước lượng được từ dãy đáp ứng tần số của mô hình h13 Ta thấy mô hình g1, g3 biểu không biểu diễn tốt đặc tính tần số của hệ thống, g2 cho đáp ứng biểu diễn chính xác hệ (*) với dải tần rộng. Chứng tỏ hệ thống thực có thể xấp xỉ tốt mô hình bậc 2. Mô hình g5 biểu diễn tốt đáp ứng của hệ thống tại dải tần số thấp nhưng lại không tốt ở dải tần sô cao. Mô hình g4 có đồ thị bode bám sát hệ thống thực do có bậc mô hình chính xác với bậc của hệ thống thực. 2.2.2. Hệ hồi quy tuyến tính và phương pháp bình phương cực tiểu (LSE) Hệ hồi quy tuyến tính Giả sử quá trình được mô tả bởi một mô hình toán học như sau: (2.25) Trong đó là giá trị quan sát được tại thời điểm , là các hàm đã xác định trước là vector tham số của mô hình cần xác định: Vector hàm được gọi là vector hồi quy, các phần tử của nó được gọi là biến hồi quy. Thông thường, vector hồi quy biểu diễn trực tiếp các dữ liệu đầu vào hoặc đầu ra đã biết trước. Mô hình trên được gọi là mô hình hồi qui tuyến tính. Nhận dạng hệ hồi quy tuyến tính theo phương pháp bình phương cực tiểu. Bài toán nhận dạng được đưa về bài toán xác định tham số mô hình sao cho sai lệch giữa các giá trị quan sát thực và các giá trị tính toán theo mô hình ước lượng là nhỏ nhất. Như tên gọi của phương pháp, tiêu chuẩn để đánh giá mức độ sai lệch này dựa trên tổng bình phương của từng sai lệch. Ta có hàm mục tiêu sau: (2.26) Trong đó: là giá trị quan sát được tại thời điểm ti là giá trị tính toán được từ quá trình nhận dạng. Ta đặt: (2.27a) (2.27b) Lúc này hệ thống (2.25) trở thành: (2.28a) (2.28b) Ước lượng tham θ trở thành giải bài toán tìm nghiệm tối ưu để cực tiểu hóa hàm V(θ, tN) (2.29) Do Φ không phụ thuộc vào θ nên (2.29) có nghiệm chính xác khi N ≥ n và là ma trận đối xứng không âm, nó không suy biến khi rank() = n tức là có đủ hạng cột. Ta có nghiệm θ của là nghiệm của phương trình: (2.30) (2.31) (2.32) trong đó là ma trận giả nghịch đảo (pseudo-inverse) của . * Khi mô hình có nhiễu tác động : (2.33) Khi đó vecto tham số sẽ được ước lượng: (2.34) (2.35) Nếu nhiễu là nhiễu không có tương quan với biến đầu vào Φ* ta có Phép ước lượng là xác thực với kỳ vọng đúng bằng của θ là: Phương sai của θ: (2.36) λ là kỳ vọng của nhiễu, Nhược điểm của phương pháp LSE: Vấn đề là khi nhiễu không phải là ồn trắng mà là nhiễu màu thì ma trận Ф và có tương quan với nhau hoặc là trong vòng kín khi u(t) là φ(t) thì φ(t) sẽ phụ thuộc vào v(t). Do đó phương pháp Least Squares sẽ có sai số.Điều này làm ảnh hưởng đến tính nhất quán của phương pháp bình phương cực tiểu, làm sai lệch kết quả nhận dạng. Để khắc phục vấn đề này, người ta sử dụng biến công cụ (Instrumental variable), là ma trận thỏa mãn hai điều kiện: không suy biến. (2.37a) (2.37b) Lúc này: (2.38) Phương pháp bình phương tối thiểu đã được xây dựng thành hàm trong matlab xác định tham số cho mô hình ARX hoặc AR. Cấu trúc lệnh: M=arx(data, orders) M=arx(data, ‘na’, na, ‘nb’, nb, ‘nk’, nk) M=arx(data, orders, ‘Property1’, value1,…, ‘PropertyN’, valueN) Với order: là ma trận tham số [na nb nc] của mô hình arx Value là phần định nghĩa các tên gọi, thuộc tính tương ứng với Property trong mô hình tham số thu được. Nhận dạng mô hình AR có thể dùng hàm arx với nb=0 hoặc dùng hàm ar có cấu trúc lệnh: M = ar(data, n) Với n là bậc của đa thức A(q) 2.2.3. Phương pháp sai số dự báo (PEM) Các hệ tuyến tính có thể được biểu diễn như sau: (2.39) Với G là hàm truyền đạt với tham số θ, có tính nhân quả chặt (2.40a) H là khâu lọc được sử dụng để mô hình hóa nhiễu: (2.40b) Bộ dự báo tuyến tính của mô hình là: (2.41) Với: (2.42) Bộ dự báo cũng có thể coi là một bộ lọc phi tuyến trong đó ta có thể viết (2.43) Với bộ dữ liệu thu thập (2.44) Bài toán ước lượng tham số trở thành bài toán xác định ánh xạ từ tập ZN → θ € DM Sai số giữa đầu ra dự báo và đầu ra thực là sai số dự báo: (2.45) Dựa vào bộ dữ liệu Z đã thu thập được, ta có thể tính sai số ε(t, θ) với t=1, 2,…, N. Xác định vectơ tham số θ sao cho dãy sai số dự báo ε(t,θ) càng nhỏ càng tốt. Mỗi tiêu chuẩn xác định độ lớn của sai số dự báo ε ta sẽ tìm ra những mô hình khác nhau thỏa mãn độ lớn đó nhỏ nhất. Sai số dự báo có thể được xem xét trên không gian RN . kích thước của véc tơ này có thể được xác định thông qua chuẩn trong không gian RN , chuẩn bình phương hoặc không. Theo đó nó đưa ra một lớp các lựa chọn tương ứng. Sai số được đi qua bộ lọc tuyến tính ổn định L(q) (2.46) Sau đó sử dụng chuẩn sau: (2.47) Với ℓ(.) là hàm vô hướng của véc εF(t,θ). Ý nghĩa của hàm L(q) có thể thấy rõ khi xem xét trong miền tần số: tập chung vào ảnh hưởng của sai số dự báo trong miền tần số thấp. Hàm VN là hàm vô hướng trên không gian RN. Lúc này ước lượng vectơ tham số θN là giá trị làm cho hàm VN nhỏ nhất. (2.48) Phương pháp xác đinh θN như trên được goi là phương pháp sao số dự báo (PEM). Mỗi cách chọn L(q) và ℓ cho chúng ta một trường hợp đặc biệt của phương pháp PEM với tên gọi riêng. Trong miền tần số với L(q)=1; ta có thể xấp xỉ hàm chuẩn VN như sau: (2.49) Phương pháp PEM còn có thể áp dụng cho mô hình đa thức dạng tổng quát như sau: (2.50) Với: Như vậy, mô hình ARX chính là (2.50) với , mô hình ARMAX chính là (2.50) với , mô hình BJ chính là (2.50) với , mô hình OE chính là (2.50) với A=C=D=1. Ngoài ra PEM còn được sử dụng để ước lượng mô hình hàm truyền đạt và mô hình không gian trạng thái. Phương pháp dự báo lỗi được xây dựng thành hàm pem, armax, oe trong Matlab identification Cấu trúc lệnh: Lệnh ước lượng mô hình không gian trạng thái: m = pem(data,n): trả về mô hình không gian trạng thái bậc n m = pem (data, ‘nx’,[n1,n2,…nN]): tính toán ước lượng mô hình không gian trạng thái với những bậc khác nhau. m = pem(data) hoặc m = pem(data,’best’): tính toán ước lượng mô hình không gian trạng thái với bậc trong khoảng [1:10] và trả về mô hình tốt nhất m=pem(data,’best’,’nk;,[]): tương tự như lệnh trên nhưng mô hình có trễ Lệnh ước lương mô hình đa thức tổng quát: m = pem(data, order): hàm trả về mô hình đa thức với bậc được cho trong vectơ order: order = [na nb nc nd nf nk] m = pem(data,.. ‘na’, na, ‘nb’, nb): hàm trả về mô hình đa thức với bậc được cho trong các tham số na, nb… Lệnh ước lượng mô hình hàm truyền đạt m = pem(data, ‘M’): trả về hàm truyền đạt được mô ta trong M = PxIDZU với: P: kí hiệu đầu tiên xác định ước lượng mô hình hàm truyền đạt, x = 0, 1, 2 hoặc 3: xác định sô hằng số thời gian, hoặc điểm cực trong mô hình I: có nếu trả về mô hình có khâu tích phân D: có nếu trả về mô hình có trễ X: có nếu trả về mô hình có điểm không U có nếu trả về mô hình có điểm cực phức. Ví dụ 2.3. Áp dụng 2 phương pháp trên ước lượng tham số cho mô hình ARX Mô hình ARX có dang: (2.51) trong đó cho trước (dùng để mô tả độ trễ của quá trình) và là nhiễu ồn trắng (white noise) và: (2.52a) (2.52b) 1. Xem xét ước lượng theo nguyên lý bình phương cực tiểu Mô hình trên có thể biểu diễn dưới dạng hồi quy tuyến tính: (2.53) Vector tham số của mô hình cần xác định: Vector hồi quy: Lúc này: (2.54) Hay: Hàm mục tiêu: (2.55) Chọn thời gian quan sát từ m đến t (với ). Lúc này nghiệm sẽ tính theo công thức : (2.56a) Với: (2.56b) (2.56c) 2. Xem xét ước lượng mô hình theo phương pháp sai số dự báo: (2.57) Bộ dự báo của hệ (2.57) là: Sai số dự báo: (2.58) Chọn L=1, hàm phí tổn trở thành (2.59) Như vậy nếu chọn L=1; thì hàm mục tiêu của 2 phương pháp LSE và PEM là giống nhau. Sự mở rộng của phương pháp PEM so với LSE 1. Mô hình nhiễu: (2.60) Khắc phục nhược điểm về tính nhất quán của LSE khi nhiễu không phải là tạp trắng hay khi nhiễu có tương quan với tín hiệu vào. 2. Hàm phí tổn: tổng quát hơn. Cho phép lưa chọn theo nhiều phương án. 3. Ưu, nhược điểm của phương pháp PEM Phức tạp hơn: Trong thuật toán tìm bộ tham số tối ưu. Đơn giản hơn: Thiết lập bài toán theo mô hình dự báo chuẩn (không phải xác định vector hồi quy). 4. LSE là một trường hợp đặc biệt của PEM với hệ hồi quy tuyến tính và chọn L=1; trong hàm phí tổn. Ví dụ 2.4: Sử dụng 2 phương pháp PEM và LSE để ước lượng mô hình cho hệ thống (*) như ở ví dụ 2.1. Ta khảo sát hệ thống khi tín hiệu vào là u2( xem thêm ở ví dụ 1) lúc này tín hiệu ra tương ứng là y2 và thu được đối tượng nhận dạng được đóng gói trong m2. Trong nhận dạng mô hình đáp ứng tần số ta thấy hệ thống (*) có thể xấp xỉ tốt về mô hình bậc [tử mẫu] là [0 2] , [1 3] hoặc [1 4] trên miền thời gian liên tục.Khi gián đoạn hóa các mô liên tục trên ta được mô hình gián đoạn có trễ là 1 và bậc tương ứng với từng mô hình là: 2, 3, 4 1. Xét mô hình khi khong chịu tác động của nhiễu. (phụ lục 2.3) Sử dụng một nửa bộ dư liệu m2 dùng dể ước lương m2e: nửa còn lại dùng để kiểm chứng mô hình m2v Ước lượng (*) theo mô hình ARX các bộ số [na nk] chọn tương ứng với bậc hệ thống là 2, 3, 4: h1=arx(m2e,'na',2,'nb',2,'nk',0); h2=arx(m2e,'na',2,'nb',2,'nk',1); h3=arx(m2e,'na',3,'nb',3,'nk',1); h4=arx(m2e,'na',4,'nb',3,'nk',1); h5=arx(m2e,'na',6,'nb',5,'nk',0); compare(m2v,h1,h2,h3,h4,h5) Hình 2.16. Ước lương (*) theo mô hình ARX khi hệ thống không chịu ảnh hưởng của nhiễu Từ đồ thị ta thấy khi không có nhiễu tác động và đã có thông tin về bậc của mô hình, chất chất lương mô hình thu được là rất cao. 2. Khi mô hình chịu tác động của nhiễu υ tác động đến đầu ra ta thu được bộ dữ liệu đầu ra y3 (phu lục 2.4) Ước lượng mô hình ARX cho hệ thống sử dụng phương pháp LSE v = wgn(2000,1,0.1) y3=y2+v; m3=iddata(y3,u2,0.1); m3e=m3(1:1000); m3v=m3(1001:2000); h1=arx(m3e,'na',2,'nb',2,'nk',0); h2=arx(m3e,'na',2,'nb',2,'nk',1); h3=arx(m3e,'na',3,'nb',3,'nk',1); h4=arx(m3e,'na',4,'nb',3,'nk',1); h5=arx(m3e,'na',6,'nb',5,'nk',0); compare(m3v,h1,h2,h3,h4,h5) Hình 2.17. Kết quả ước lương mô hình ARX khi hệ thống chịu ảnh hưởng của nhiễu trắng tới đầu ra sử dụng phương pháp LSE Khi có nhiễu ảnh hưởng tới đầu ra, mô hình ước lượng được có độ chính xác không cao cho dù chọn được đúng bậc mô hình. Sử dụng phương pháp PEM ước lương mô hình ARX cho hệ thống h_1=pem(m3e,'na',2,'nb',2,'nk',0); h_2=pem(m3e,'na',2,'nb',2,'nk',1); h_3=pem(m3e,'na',3,'nb',3,'nk',1); h_4=pem(m3e,'na',4,'nb',3,'nk',1); h_5=pem(m3e,'na',6,'nb',5,'nk',0); compare(m3v,h_1,h_2,h_3,h_4,h_5) Hình 2.18. Kết quả ước lượng (*) theo mô hình ARX sử dung phương pháp PEM Ta thấy sử dụng phương pháp PEM thu được mô hình có độ fit cao hơn phương pháp LSE nhưng kết quả vẫn còn chưa được tốt, điều này chứng tỏ mô hình ARX chưa phù hợp để biểu diễn hệ thống. Ước lượng hệ thống (*) sử với mô hình armax sử dụng phương pháp PEM h_11=pem(m3e,'na',2,'nb',2,'nc',1,'nk',0); h_12=pem(m3e,'na',2,'nb',2,'nc',1,'nk',1); h_13=pem(m3e,'na',3,'nb',3,'nc',1,'nk',1); h_14=pem(m3e,'na',4,'nb',3,'nc',1,'nk',1); h_15=pem(m3e,'na',6,'nb',5,'nc',1,'nk',0); compare(m3v,h_11,h_12,h_13,h_14,h_15) Hình 2.19. Kết quả ước lương (*) theo mô hình ARMAX khi đầu ra chịu tác động của nhiễu Ta thấy với mô hình ARMAX chất lượng mô hình ước lương được tốt hơn so với mô hình ARX, chứng tỏ khi hệ thống chịu ảnh hưởng của nhiễu, mô hình ARMAX phù hợp hơn để biểu diễn hệ thống. Ước lượng mô hình hàm truyền đạt sử dụng phương pháp PEM h_21=pem(m3e,'P1'); h_22=pem(m3e,'P1D'); h_23=pem(m3e,'P2'); h_24=pem(m3e,'P2Z'); h_25=pem(m3e,'P3'); h_26=pem(m3e,'P3Z'); compare(m3v,h_21,h_22,h_23,h_24,h_25,h_26 Hình 2.20. Kết quả ước lương (*) theo mô hình hàm truyền đạt khi đầu ra chịu tác động của nhiễu Từ kết quả trên ta thấy khi ước lương hàm truyền đạt cho hệ thống, khi chọn xấp xỉ hệ thống với mô hình thấp hơn hệ thống thực, ta có thể sử dụng khâu trễ để tăng độ chính xác của mô hình. ( h21 à h22). Kết quả ước lượng cho hàm truyền đạt sử dụng phương pháp PEM là rất cao cho dù có nhiễu tác đông, đặc biệt khi biết được chính xác bậc của hệ thống. Từ ví dụ trên ta thấy khi hệ thống tuyến tính không chịu ảnh hưởng của nhiễu thì mô hình ARX có rât phù hợp để mô hình hệ thống, và phương pháp ước lương LSE cho kết quả rất tốt. Nhưng khi hệ thống chịu tác động của nhiễu thì phương pháp PEM cho kết quả tốt hơn LSE với cùng một mô hình ARX. Mặc dù vậy mô hình ARX tỏ ra không thích hợp để mô hình hóa hệ thống khi có nhiễu, các mô hình ARMAX, mô hình hàm truyền đạt có cho chất lượng tốt hơn rất nhiều đăc biết khi có thông tin về bậc của hệ thống. CHƯƠNG 3: ĐỐI TƯỢNG LÒ HƠI Chương này trình bày tóm tắt về quá trình công nghệ cũng như các sách lược điều khiển hiện tại trong lò hơi để làm rõ các biến vào ra của đối tượng phục vụ cho mục đích nhận dạng. Nội dung trong phần này chủ yếu được trích dẫn trong tài liệu [1]. 3.1. Quá trình công nghệ lò hơi Lò hơi là thiết bị tạo ra hơi nước bão hòa hoặc hơi nước quá nhiệt. Hơi nước quá nhiệt dùng để làm nguồn năng lượng cung cấp cho các thiết bị quay (rotate device) tại các nhà máy sản xuất công nghiệp như turbine truyền động bơm hoặc máy nén...hay dẫn động các turbine để quay các máy phát điện. Bên cạnh việc tạo ra động năng, hơi nước quá nhiệt này còn có thể sử dụng trong một vài ứng dụng khác như làm khô sản phẩm hay gia nhiệt chất xúc tác… Hình 3.1. Sơ đồ cấu tạo của lò hơi Lò hơi 10-B-8001 tại nhà máy đạm Phú Mỹ có thể cung cấp công suất tối đa là 140 tấn/h hơi quá nhiệt ở nhiệt độ 380 ± 5 oC và áp suất 39 ± 0.5 bar. Nhiên liệu là khí đốt thiên nhiên (natural gas) được cung cấp bởi trạm cung cấp khí GDC của tổng công ty khí Việt Nam (PVGAS). Chu trình hoạt động diễn ra bên trong lò hơi như sau: Nước khử khoáng sau khi qua bộ trao đổi nhiệt (Economiser) sẽ đi vào bao hơi, sau đấy nó sẽ trao đổi nhiệt tại buồng lửa thông qua bộ gia nhiệt (Evaporator) và trở về lại bao hơi. Lúc này trong bao hơi sẽ là hỗn hợp giữa nước và hơi nước. Do đó mức nước sẽ thay đổi phụ thuộc vào áp suất trong bao hơi. Nếu áp suất giảm thì mức sẽ tăng và ngược lại. Dòng hơi bão hòa ra khỏi bao hơi sẽ đi vào Superheater. Bộ Superheater này bao gồm hai dàn trao đổi nhiệt. Sau khi trao đổi nhiệt tại dàn trao đổi nhiệt thứ nhất của Superheater, hơi quá nhiệt sẽ được làm mát bằng một lượng nước làm mát (quench water) để điều hòa nhiệt độ và tiếp tục đi vào dàn trao đổi nhiệt thứ hai. Hơi nước ra khỏi dàn trao đổi nhiệt thứ hai này (tức là ra khỏi Superheater) chính là hơi nước quá nhiệt được sử dụng cho quá trình sản xuất. Lượng hơi này sẽ được đưa vào mạng hơi, tùy vào lượng hơi được sử dụng trong quá trình sản xuất (downstream) mà lò hơi thay đổi công suất cho phù hợp. Trong quá trình vận hành, ta cần phải duy trì mực nước trong bao hơi ở mức độ an toàn (thông thường là giữa bao hơi): không thấp quá để đảm bảo đủ nước cho quá trình tạo hơi và tránh gây cháy các ống trao đổi nhiệt mà cũng không cao quá để ổn định được áp suất trong bao hơi khi tải (nhu cầu tiêu thụ hơi nước quá nhiệt) thay đổi. Do vậy khi mực nước nằm trong giới hạn an toàn ta còn phải bơm liên tục một lượng nước vào bao hơi bằng với lượng hơi được tiêu thụ ở đầu ra của lò hơi để đảm bảo mực nước không thay đổi. Để gia nhiệt cho nước thành hơi nước quá nhiệt, nguồn nhiệt năng được tạo từ buồng lửa bằng việc đốt cháy nhiên liệu. Chất lượng quá trình cháy được kiểm soát bằng nồng độ Oxy ở khí thải. 3.2. Giải pháp điều khiển đang được sử dụng Ta nhận thấy, việc điều khiển lò hơi được chi phối bởi bốn vấn đề chính: i. Mức nước trong bao hơi. ii. Nhiệt độ hơi quá nhiệt ở ngõ ra. iii. Áp suất hơi quá nhiệt ở ngõ ra. iv. Chất lượng quá trình cháy trong buồng lửa. Trong khuôn khổ đồ án này sẽ tập trung vào điều khiển chất lượng hơi quá nhiệt tức là vòng điều khiển nhiệt độ và áp suất nhằm xây dựng mô hình động học của lò hơi như mong muốn. 3.2.1. Điều khiển mức nước trong bao hơi Việc điều khiển mức nước trong bao hơi được thực hiện thông qua bộ điều khiển LIC8250 với cảm biến mức dạng chênh áp LIT8250 và van điều khiển LV8250. Hình 3.2. Sơ đồ điều khiển mức trong bao hơi Mức nước trong bao hơi chủ yếu được duy trì bằng cách giữ cân bằng vật chất vào ra của lò hơi (tức là khối lượng nước vào và khối lượng hơi ra phải bằng nhau). Quá trình này được điều tiết bằng vòng điều khiển cascade bao gồm 2 bộ điều khiển FIC8251 và LIC8250A. Bộ điều khiển LIC8250 chỉ hoạt động khi mức nước trong bao hơi xuống quá thấp, lúc đấy LIC8250 sẽ mở hoàn toàn van LV8250 để đưa nước vào để mức nước trong bao hơi trở về giá trị phù hợp (tuy nhiên lúc này thì lò hơi không ở trong trạng thái hoạt động ổn định nữa, do vậy ta không xét đến trạng thái này của lò hơi). Bởi vì tính độc lập giữa việc duy trì mức nước trong bao hơi và chất lượng hơi quá nhiệt ở ngõ ra, ta không đưa vòng điều khiển này vào mục tiêu nhận dạng của đồ án này. 3.2.2. Kiểm soát nhiệt độ hơi nước quá nhiệt Dòng hơi nước ra khỏi bao hơi sẽ đi vào bộ trao đối nhiệt Superheater, sau khi qua bộ Superheater 1 (10-B-8001/SH1) sẽ đi qua bộ Desuperheater và đi vào Superheater 2 (10-B-8001/SH2). Hình 3.3. Sơ đồ điều khiển nhiệt độ của hơi quá nhiệt Nhiệt độ hơi quá nhiệt ở đầu ra qua bộ điều khiển TIC8253 cùng với lưu lượng nhiên liệu được đo từ bộ FI8201 tạo ra giá trị đặt cho bộ điều khiển lưu lượng FIC8252. Đồng thời lưu lượng nước làm mát được đo qua bộ FIT8252 sẽ được so sánh với giá trị đặt ở bộ điều khiển lưu lượng FIC8252 để tạo tín hiệu điều khiển van. 3.2.3. Áp suất hơi quá nhiệt ở đầu ra Việc duy trì áp suất hơi quá nhiệt ở đầu ra (khoảng 39 bar) được thực hiện thông qua việc duy trì nhiệt năng cung cấp cho lò hơi mà cụ thể chính là năng lượng từ quá trình đốt cháy nhiên liệu. Hình 3.4. Sơ đồ điều khiển áp suất hơi quá nhiệt đầu ra Bộ điều khiển PIC4048 có nhiệm vụ duy trì áp suất của hơi quá nhiệt ở đầu ra của lò hơi. Giá trị đặt cho bộ PIC4048 sẽ do người vận hành nhập vào cùng với lưu lượng hơi quá nhiệt từ FI8253 sẽ tạo ra giá trị đặt cho lưu lượng nhiên liệu đầu vào FIC8201. Lưu lượng nhiên liệu qua bộ FIT 8201 sẽ so sánh với giá trị đặt để tạo tín hiệu điều khiển van FV 8201. 3.2.4. Chất lượng quá trình cháy trong buồng lửa Quá trình cháy trong buồng lửa luôn được quan tâm nhiều trong vận hành lò hơi, nếu quá trình cháy tốt thì lò hơi sẽ tận dụng được nhiều năng lượng hơn còn nếu quá trình cháy không tốt thì lượng nhiên liệu tiêu hao rất đáng kể. Hình 3.5. Sơ đồ điều khiển chất lượng quá trình cháy Chất lượng của quá trình đốt cháy nhiên liệu được đánh giá thông qua nồng độ khí Ôxy bên trong khí thải (flue gas). Công việc kiểm soát nồng độ Oxy trong khí thải được thực hiện bằng bộ phân tích online AT8250. Bộ điều khiển nồng độ khí ôxy AIC8250 sẽ thực hiện điều tiết tỉ lệ lưu lượng không khí FIT8250 và lưu lượng nhiên liệu FIT8201 để quá trình cháy đạt được chất lượng tốt nhất bằng cách thay đổi hệ số nhân biến thiên từ 0.85 tới 1.15. Để duy trì lượng O2 trong khí thải bằng cách điều chỉnh lượng nhiên liệu vào, ta chỉ cần sử dụng một bộ điều khiển tỉ lệ sao cho nồng độ O2 đúng như mong đợi (hình 3.6). Đây là vòng điều khiển tương đối độc lập với bài toán thiết lập mô hình nhiệt độ và áp suất của lò hơi nên ta cũng không xét đến. Hình 3.6. Nồng độ Ôxy tối ưu đối với từng mức tải của lò hơi Tóm lại: Việc nhận dạng lò hơi chỉ xét đến hai vòng điều khiển nhiệt độ và áp suất của hơi quá nhiệt. Hình 3.7. Quá trình và các biến quá trình được sử dụng trong mô hình Lúc này hai biến đầu vào là lưu lượng nhiên liệu ở đầu vào buồng đốt FI8201( điều chỉnh áp suất) và lưu lượng nước làm mát FI8252 (điều chỉnh nhiệt độ), lưu lượng hơi quá nhiệt phụ thuộc vào yêu cầu của quá trình tiếp theo không điều khiển được nên ta coi là nhiễu tải. Ngoài ra nhiễu đo và các ảnh hưởng không mong muốn khác từ môi trường được coi như nhiễu ồn trắng. Hai biến ra là nhiệt độ TI8253 và áp suất của hơi quá nhiệt PI4048 ở đầu ra. 3.3. Thu thập số liệu Hệ thống lò hơi được duy trì trong trạng thái vận hành, các vòng điều khiển vẫn để ở chế độ tự động. Chính vì vậy mà dữ liệu sẽ được thu thập theo phương pháp nhận dạng trực tiếp trong vòng kín. Việc thu thập dữ liệu được bắt đầu bằng việc thay đổi điểm đặt của từng vòng điều khiển (ở chế độ tự động) trong khi vòng điều khiển còn lại bị khống chế ở chế độ vận hành bằng tay. Các thay đổi của các biến quá trình và biến điều khiển được ghi nhận với chu kỳ lấy mẫu là 1 giây. Như vậy với mỗi lần làm cho một vòng điều khiển ta thu được dữ liệu mô tả ảnh hưởng của các biến vào tới từng biến đầu ra. Các dữ liệu sau khi được thu thập cần phải loại bỏ giá trị trung bình để thực hiện nhận dạng. 3.3.1. Vòng điều khiển nhiệt độ Hình 3.8. Sơ đồ khối của vòng điều khiển nhiệt độ TIC8253 Ta thực hiện thay đổi giá trị đặt của bộ điều khiển nhiệt độ TIC8253, giữ cố định giá trị đặt cho PIC8201 để chỉ xét ảnh hưởng của các biến vào tới nhiệt độ hơi quá nhiệt ngõ ra. Các biến vào chính là FIC8201.PV (lưu lượng nhiên liệu), FIC8252.PV (lưu lượng nước làm mát), FI8253 (lưu lượng hơi quá nhiệt), biến ra là TIC8253.PV(nhiệt độ của hơi quá nhiệt). Vòng điều khiển thứ cấp FIC8252 là vòng điều khiển van cấp nước làm mát với biến ra là FIC8252.PV. Như vậy, trong mô hình nhiệt độ của lò hơi ta có thể bỏ qua vòng này khi xét ảnh hưởng của lưu lượng nước làm mát tới nhiệt độ hơi quá nhiệt. Nhận dạng mô hình nhiệt độ là xác định hàm truyền đạt G11 của lưu lượng nước làm mát, G21 của lưu lượng nhiên liêu, G31 của lưu lượng hơi quá nhiệt đến biến nhiệt độ. 3.3.2. Vòng điều khiển áp suất Hình 3.9. Sơ đồ khối của vòng điều khiển áp suất PIC4048 Ta thực hiện thay đổi điểm đặt của bộ điều khiển áp suất PIC4048. Trong lúc đấy, điểm đặt của TIC8253 được giữ cố định để chỉ xét ảnh hưởng của các biến vào tới áp suất hơi quá nhiệt đầu ra. Các biến vào chính là FIC8201.PV (lưu lượng nhiên liệu), FIC8252.PV (lưu lượng nước làm mát), FI8253 (lưu lượng hơi quá nhiệt), biến ra là PIC4048.PV(nhiệt độ của hơi quá nhiệt). Vòng điều khiển thứ cấp FIC8201 là vòng điều khiển van cấp lưu lượng nhiên liệu với biến ra là FIC8201.PV. Như vậy, trong mô hình áp suất của lò hơi ta có thể bỏ qua vòng này khi xét ảnh hưởng của lưu lượng nhiên liệu tới áp suất hơi quá nhiệt. Nhận dạng mô hình áp suất là xác định hàm truyền đạt G12 của lưu lượng nước làm mát, G22 của lưu lượng nhiên liêu, G32 của lưu lượng hơi quá nhiệt đến biến ra áp suất. CHƯƠNG 4: NHẬN DẠNG MÔ HÌNH LÒ HƠI 4.1. Xử lý số liệu trước khi nhận dạng Để cho kết quả nhận dạng được chính xác và phản ánh đúng đặc tính động học của đối tượng thì trước khi tiến hành nhận dạng ta cần xử lý số liệu từ bộ số liệu thô đã thu thập được. Gồm hai bước cơ bản sau: Loại bỏ giá trị trung bình Chia bộ số liệu làm hai phần: một phần dùng để ước lượng, một phần dùng để kiểm chứng. Các hàm hỗ trợ trong Matlab được sử dụng như sau: Iddata Hàm iddata được sử dụng để tạo ra đối tượng dữ liệu miền thời gian hoặc miền tần số dùng cho nhận dạng trong công cụ System Identification Toolbox của Matlab. Cấu trúc câu lệnh: data = iddata(y,u,Ts) Trong đó: y là ma trận đại diện cho biến ra của hệ thống. Nếu hệ thống chỉ có một biến ra thì y là một vector cột, nếu hệ thống có nhiều biến ra thì y là ma trận có số cột bằng số biến ra. u là ma trận đại diện cho các biến vào của hệ thống. Nếu hệ thống chỉ có một biến vào thì u là một ma trận cột, nếu hệ thống có nhiều biến vào thì u là ma trận có số cột bằng số biến vào. Ts là thời gian lấy mẫu của bộ dữ liệu vào ra [y , u] Detrend Hàm detrend được sử dụng để loại bỏ giá trị trung bình và làm cho bộ số liệu rời rạc theo thời gian thu được từ thực nghiệm có xu hướng tuyến tính. Cú pháp: Để loại bỏ giá trị trung bình từ bộ dữ liệu data được tạo ra từ câu lệnh iddata, sử dụng cấu trúc lệnh: data = detrend(data); Để trừ đi xu hướng tuyến tính bậc 1, ta sử dụng cú pháp: data = detrend(data,1); Delayest Hàm delayest thực hiện công việc dự đoán thời gian trễ giữa tín hiệu input và output trong đối tượng iddata sử dụng mô hình ARX Cú pháp: nk = delayest(Data) nk = delayest(Data,na,nb,nkmin,nkmax,maxtest) Trong đó: Data là đối tượng iddata, chứa các thông tin về các biến vào ra. na là số phần tử của đầu ra y được sử dụng trong mô hình ARX nb là số phần tử của đầu vào u được sử dụng trong mô hình ARX. Nếu có nhiều đầu vào nb là ma trận có số cột bằng số đầu vào. nkmin và nkmax xác định giới hạn thời gian trễ cần kiểm tra. maxtest xác định số mẫu được sử dụng để thực hiện dự đoán trễ. Compare Dùng để so sánh dữ liệu mô phỏng ở ngõ ra của mô hình với dữ liệu ngõ ra thu thập được. Cú pháp thường dùng: Compare(data,m) Trong đó: data là đối tượng iddata chứa thông tin vào/ra để nhận dạng. m chính là mô hình thu được từ việc nhận dạng data. Ta sử dụng bộ số liệu EXP1 RAW DATA để nhận dạng và kiểm chứng mô hình nhiệt độ. Bộ số liệu này gồm 5000 mẫu, ta chia làm hai phần, phần 1 từ [1: 2000] sẽ dùng để ước lượng. Từ [2001:5000] sẽ sử dụng để kiểm chứng. Việc đóng gói dữ liệu và loại bỏ giá trị trung bình cho các biến vào ra xem phụ lục kèm theo. Xử lý với đầu ra nhiệt độ ta có: >>T=expr1(:,5); >>plot(T) ; Hình 4.1a. Đồ thị nhiệt độ khi chưa loại bỏ giá trị trung bình >>T1= detrend (T,1); >>Plot(T1); Hình 4.1b. Đồ thị nhiệt độ sau khi loại bỏ giá trị trung bình Ta sử dụng bộ số liệu EXP2 RAW DATA để nhận dạng và kiểm chứng mô hình áp suất. Bộ số liệu này gồm 4200 mẫu, ta chia làm hai phần, phần 1 từ [1: 2000] sẽ dùng để ước lượng. Từ [2001:4200] sẽ sử dụng để kiểm chứng. Việc đóng gói dữ liệu và loại bỏ giá trị trung bình cho các biến vào ra xem phụ lục kèm theo. Xử lý với đầu ra áp suất ta có: >> P2=expr2(:,4); >>plot(P2) ; Hình 4.2a. Đồ thị áp suất khi chưa loại bỏ giá trị trung bình >> P2e=detrend(P2,1); >> plot(P2e); Hình 4.2b. Đồ thị áp suất hơi quá nhiệt sau khi loại bỏ giá trị trung bình Sau đây, ta sẽ tiến hành nhận dạng và kiểm chứng mô hình nhiệt độ và mô hình áp suất theo 3 phương pháp và với các dạng mô hình khác nhau. 4.2. Kết quả nhận dạng cho mô hình nhiệt độ Vòng điều khiển nhiệt độ là vòng kín có biến ra là nhiệt độ, biến vào là lưu lượng nước, lưu lượng nhiên liệu và nhiễu là lưu lượng hơi quá nhiệt. Lưu lượng hơi quá nhiệt Nhiệt độ hơi quá nhiệt Lưu lượng nhiên liệu Quá trình nhiệt Lưu lượng nước làm mát SP nhiệt độ -- Hình 4.3. Mô hình nhận dạng nhiệt độ 4.2.1. Nhận dạng mô hình nhiệt độ theo phương pháp dựa trên đáp ứng trên miền tần số Phân tích phổ tần số của các biến vào dùng phép biến đổi Fourier (phụ lục 4.1) Ta thu được phổ tần số như sau: Hình 4.4. Phổ tần số của các biến vào trong mô hình nhiệt độ Các biến vào của mô hình có dải tần số hẹp phân bố ở vùng tần số thấp. Với dải tần số trên, nhận dạng mô hình bằng phương pháp phân tích phổ sẽ không có kết quả chính xác ở vùng tần số cao, chỉ phản ánh được đặc tính của hệ thống ở vùng tần số thấp. Ước lượng theo các phương pháp dựa trên đáp ứng tần số >>T1= etfe(Tmhe); % uoc luong bang phuong phap etfe >>T2=spafdr(Tmhe);% uoc luong bang phuong phap spafdr >>T3=spa(Tmhe); % uoc luong bang phuong phap spa >>Bode(T1, 'r', T2,'b', T3,'g'); ≈20dB Hình 4.5a. Đáp ứng tần số của nhiệt độ ứng với biến vào lưu lượng nhiên liệu ≈ 40dB Hình 4.5b. Đáp ứng tần số của nhiệt độ với biến vào lưu lượng nước làm mát ≈ 60dB Hình 4.5c. Đáp ứng tần số của nhiệt độ với biến vào lưu lượng hơi quá nhiệt Nhận xét: Ta thấy cả 3 phương pháp đều cho mô hình không trung thực ở vùng tần số cao. Lí do là các biến vào là các tín hiệu không giàu tần số tức là chỉ tập trung ở dải tần số thấp. Thay đổi độ rộng của hàm cửa sổ trong phép ước lượng etfe và phương pháp spa cũng không cải thiện được chất lượng của mô hình ở vùng tần số ω > 0.1 rad. Từ mô hình đáp ứng tần số ước lượng được ở trên có thể ước lượng được bậc của mô hình nhiệt độ với biến vào là lưu lượng nhiên liệu ≥ 1; với biến vào là lưu lượng nước làm mát ≥ 2; với biến vào là lưu lượng hơi quá nhiệt ≥ 3. 4.2.2. Nhận dạng mô hình nhiệt độ theo phương pháp LSE Để có được một mô hình phù hợp với mục đích điều khiển ta bắt đầu từ mô hình ARX với bậc của mô hình thấp. >>M1=arx(Tmhe,'na',2,'nb',[2 2 2],'nk',[1 1 1]); >>compare(Tmhv,M1); grid; Kết quả thu được như sau: Hình 4.6a. Nhận dạng nhiệt độ theo phương pháp PEM cho mô hình ARX với bậc mô hình thấp Ta thấy tuy bậc mô hình thấp (bậc 2) nhưng chất lượng đạt được rất kém (30.39%). Ta tiếp tục tăng bậc của mô hình tới khi: >>M2=arx(Tmhe,'na',1,'nb',[85 85 85],'nk',[3 3 3]); >>compare(Tmhv,M2); grid; Ta thu được kết quả như sau: Hình 4.6b. Nhận dạng nhiệt độ theo phương pháp PEM cho mô hình ARX với bậc mô hình cao Nhận xét: Ta thấy chất lượng mô hình khá tốt (92.58%) tuy nhiên bậc của mô hình lại rất cao (bậc 85) vì vậy mô hình này không phù hợp để thiết kế sách lược điều khiển. Điều này có thể giải thích như sau: Mô hình nhận dạng ARX với phương pháp LSE chỉ phù hợp với các đối tượng không chịu ảnh hưởng của nhiễu, hoặc chỉ bị ảnh hưởng của nhiễu đo vào đầu ra. Đối tượng thực tế ngoài nhiễu do thiết bị đo còn có các nhiễu quá trình khác ảnh hưởng đến hệ thống như nhiệt độ môi trường, chất lượng của nhiên liệu… Chính vì vậy mô hình ARX không phù hợp để mô ta quá trình nhiệt của lò hơi. 4.2.3. Nhận dạng mô hình nhiệt độ theo phương pháp PEM cho các dạng mô hình khác nhau. Sử dụng phương pháp PEM để nhận dạng cho mô hình ARX >>M3=pem(Tmhe,'na',2,'nb',[2 2 2],'nk',[1 1 1]); >>compare(Tmhv,M3); grid; Ta thu được kết quả với độ fit là 30.39% giống như việc áp dụng phương pháp LSE cho mô hình ARX và chất lượng mô hình rất thấp. Lí do là mô hình ARX chưa xét đến mô hình của nhiễu. Sử dụng phương pháp PEM để nhận dạng cho mô hình ARMAX. Ta sử dụng các bộ tham số khác nhau như sau: M4=pem(Tmhe,'na',2,'nb',[2 2 2],'nc',2,'nk',[1 1 1]); M5=pem(Tmhe,'na',3,'nb',[3 3 3],'nc',2,'nk',[1 1 1]); compare(Tmhv,M4,M5); grid; Mô hình nhiệt độ thu được theo phương pháp PEM cho mô hình ARMAX ( mô hình M4) như sau: Discrete-time IDPOLY model: A(q)y(t) = B(q)u(t) + C(q)e(t) A(q) = 1 - 1.978 q^-1 + 0.9776 q^-2 B1(q) = -2.667e-005 q^-1 + 3.396e-005 q^-2 B2(q) = 0.005461 q^-1 - 0.007149 q^-2 B3(q) = -6.136e-006 q^-1 + 5.825e-006 q^-2 C(q) = 1 - 1.375 q^-1 + 0.3801 q^-2 Estimated using PEM using SearchMethod = Auto from data set z Loss function 0.00230978 and FPE 0.0023329 Sampling interval: 1 Hình 4.7. Nhận dạng nhiệt độ theo phương pháp PEM cho mô hình ARMAX Ta thấy, khi sử dụng phương pháp PEM cho mô hình ARMAX, với các bộ tham số khác nhau ta thu được chất lượng mô hình gần giống nhau. Ta tiếp tục sử dụng phương pháp PEM cho các dạng mô hình khác nhau như sau: % mo hinh quan tinh bac 2 co tre, co diem khong M6=pem(Tmhe,'p2dz','Td',{'max',2}); % mo hinh state space M7=pem(Tmhe,4); % mo hinh OE M8=oe(Tmhe,'nb',[1 1 1],'nf',[2 2 2],'nk',[1 1 0]); compare(Tmhv,M6,M7,M8); grid; Mô hình nhiệt độ ước lượng theo phương pháp PEM cho mô hình quán tính bậc 2 có trễ, có điểm không Process model with 3 inputs: y = G_1(s)u_1 + ...+ G_3(s)u_3 where 1+Tz*s G_1(s) = Kp * ------------------ * exp(-Td*s) (1+Tp1*s)(1+Tp2*s) with Kp = 0.056577 Tp1 = 79.347 Tp2 = 79.38 Td = 2 Tz = -1.5035 1+Tz*s G_2(s) = Kp * ------------------ * exp(-Td*s) (1+Tp1*s)(1+Tp2*s) with Kp = -14.177 Tp1 = 144.89 Tp2 = 52.242 Td = 0.10512 Tz = -8.734 1+Tz*s G_3(s) = Kp * ------------------ * exp(-Td*s) (1+Tp1*s)(1+Tp2*s) with Kp = -0.0033729 Tp1 = 362.88 Tp2 = 2.4704 Td = 0 Tz = -15.489 Estimated using PEM using SearchMethod = Auto from data set z Loss function 0.0107199 and FPE 0.0108807 Mô hình nhiệt độ ước lượng theo phương pháp PEM cho mô hình không gian trạng thái State-space model: x(t+Ts) = A x(t) + B u(t) + K e(t) y(t) = C x(t) + D u(t) + e(t) A = x1 x2 x3 x4 x1 0.9845 0.0038504 0.010116 -0.014505 x2 0.0021587 1.1149 0.25222 -0.42929 x3 -0.0097151 -0.17313 -0.14724 -0.32387 x4 0.0054164 0.037878 -0.87939 -0.24233 B = u1 u2 u3 x1 2.3716e-006 -0.00047634 -2.049e-007 x2 -0.00010391 0.026915 8.7695e-008 x3 6.5788e-005 0.0037466 -8.1436e-005 x4 0.00010974 0.0047669 -0.00010742 C = x1 x2 x3 x4 y1 94.921 0.79819 -0.12301 0.44307 D = u1 u2 u3 y1 0 0 0 K = y1 x1 0.0067272 x2 -0.030863 x3 -0.04828 x4 0.011882 x(0) = x1 -0.0096244 x2 0.069502 x3 -0.046093 x4 -0.0095877 Estimated using PEM using SearchMethod = Auto from data set z Loss function 0.00223677 and FPE 0.00229045 Sampling interval: 1 Mô hình nhiệt độ ước lượng theo phương pháp PEM cho mô hình OE Discrete-time IDPOLY model: y(t) = [B(q)/F(q)]u(t) + e(t) B1(q) = 0.0006304 q^-1 B2(q) = -0.001031 q^-1 B3(q) = -2.62e-006 F1(q) = 1 - 0.07277 q^-1 - 0.9167 q^-2 F2(q) = 1 - 1.987 q^-1 + 0.9872 q^-2 F3(q) = 1 - 1.807 q^-1 + 0.8083 q^-2 Estimated using OE from data set z Loss function 0.0370939 and FPE 0.0374281 Sampling interval: 1 Kết quả nhận được như sau: Hình 4.8. Nhận dạng nhiệt độ theo phương pháp PEM cho mô hình SOPDT, không gian trạng thái và OE Nhận xét: Mô hình OE theo phương pháp PEM thu được độ fit là 88.24%, bậc của mô hình cũng tương đối thấp. Mô hình hàm truyền bậc 2 có trễ, có điểm 0 và mô hình State Space thu được độ fit lần lượt là 84.81% và 83.18% . Ta thấy kết quả nhận dạng theo phương pháp PEM cho các dạng mô hình khác nhau cho kết quả tương đối chính xác so với phương pháp LSE bên trên. Lí do là: phương pháp sai số dự báo đã khắc phục nhược điểm về tính nhất quán của LSE khi nhiễu không phải là tạp trắng hay khi nhiễu có tương quan với tín hiệu vào. 4.3. Kết quả nhận dạng cho mô hình áp suất Vòng điều khiển áp suất hơi quá nhiệt là vòng kín có biến ra là áp suất, biến vào là lưu lượng nước, lưu lượng nhiên liệu và nhiễu là lưu lượng hơi quá nhiệt. Lưu lượng hơi quá nhiệt Áp suất hơi quá nhiệt Lưu lượng nước làm mát Quá trình áp suất SP áp suất Lưu lượng nhiên liệu - Hình 4.9. Mô hình nhận dạng áp suất 4.3.1. Nhận dạng mô hình áp suất theo phương pháp phân tích phổ Phổ tần số của các biến vào: Hình 4.10. Phổ tần số của các biến vào trong mô hình áp suất Các biến vào trong mô hình áp suất cũng có dải tần hẹp tập trung ở vùng tần số thấp chính vì vậy khi dùng để nhận dạng mô hình đáp ứng tần số sẽ không thu được kết quả trung thực ở vùng tần số cao. Ước lượng theo các phương pháp dựa trên đáp ứng tần số >>P1= etfe(Pmhe); % uoc luong bang phuong phap etfe >>P2=spafdr(Pmhe);% uoc luong bang phuong phap spafdr >>P3=spa(Pmhe); % uoc luong bang phuong phap spa >>Bode(P1, 'r', P2,'b', P3,'g') ≈ 20dB Hình 4.11a. Đáp ứng tần số của áp suất với biến vào lưu lượng nhiên liệu ≈ 40dB Hình 4.11b. Đáp ứng tần số của áp suất với biến vào lưu lượng nước làm mát ≈ 40dB Hình 4.11c. Đáp ứng tần số của áp suất với biến nhiễu tải lưu lượng hơi quá nhiệt Nhận xét: Cũng như mô hình nhiệt độ, cả 3 phương pháp đều không phản ánh trung thực mô hình tại vùng tần số cao. Lí do là các biến vào là các tín hiệu không giàu tần số tức là chỉ tập trung ở dải tần số thấp. Do đó ước lượng ở dải tần số cao không chính xác. Dựa vào đáp ứng theo mô hình ước lượng bởi phương pháp SPAFDR ta có thể ước lượng bậc của mô hình với các biến vào như sau: biến vào là lưu lượng nhiên liệu ≥ 1; với biến vào là lưu lượng nước làm mát ≥ 2; với biến vào là lưu lượng hơi quá nhiệt ≥ 3. 4.3.2. Nhận dạng mô hình áp suất theo phương pháp LSE với mô hình ARX Cũng tương tự như cách làm với mô hình nhiệt độ. Để có được một mô hình phù hợp với mục đích điều khiển ta bắt đầu từ mô hình ARX với bậc của mô hình thấp. m1=arx(Pmhe,'na',2,'nb',[3 3 3],'nk',[1 1 1]); compare(Pmhv,m1);grid; Hình 4.12a. Nhận dạng áp suất theo phương pháp LSE áp dụng mô hình ARX với bậc của mô hình thấp Ta thấy tuy bậc mô hình thấp (bậc 3) nhưng chất lượng đạt được chưa được tốt lắm (66%). Ta tiếp tục tăng bậc của mô hình tới khi: >>m2=arx(Pmhe,'na',2,'nb',[88 88 88],'nk',[3 3 3]); compare(Pmhv,m2);grid; Hình 4.12b. Nhận dạng áp suất theo phương pháp LSE áp dụng mô hình ARX với bậc của mô hình cao Ta thấy chất lượng mô hình khá tốt (86.84%) tuy nhiên bậc của mô hình lại rất cao (bậc 88) vì vậy mô hình này không phù hợp để thiết kế sách lược điều khiển. Điều này có thể giải thích như sau: nếu như coi đầu ra nhiệt độ hơi quá nhiệt không chịu tác động của thành phần nhiễu (hoặc nhiễu là ồn trắng) thì mô hình ARX cho kết quả khá chính xác. Tuy nhiên, thực tế đầu ra nhiệt độ bao gồm thành phần của nhiều do đó phương pháp LSE không còn tính nhất quán, làm sai lệch kết quả nhận dạng. 4.3.3. Nhận dạng mô hình áp suất theo phương pháp PEM cho các dạng mô hình khác nhau Sử dụng phương pháp PEM để nhận dạng cho mô hình ARX m3=pem(Pmhe,'na',2,'nb',[3 3 3],'nk',[1 1 1]); Ta thu được kết quả với độ fit là 66% giống như việc áp dụng phương pháp LSE cho mô hình ARX và chất lượng mô hình chưa cao. Sử dụng phương pháp PEM để nhận dạng cho mô hình ARMAX. Ta sử dụng các bộ tham số khác nhau như sau: >>m4=pem(Pmhe,'na',2,'nb',[3 3 3],'nc',2,'nk', [1 1 1]); m5=pem(Pmhe,'na',1,'nb',[2 2 2],'nc',2,'nk', [1 1 1]); compare(Pmhv,m4,m5); grid; Mô hình ARMAX bậc na = 2; nb = 3; nc = 2; nk = 1 Discrete-time IDPOLY model: A(q)y(t) = B(q)u(t) + C(q)e(t) A(q) = 1 - 1.96 q^-1 + 0.9606 q^-2 B1(q) = 8.901e-005 q^-1 - 0.0001463 q^-2 + 5.76e-005 q^-3 B2(q) = 0.003282 q^-1 + 0.001472 q^-2 - 0.005666 q^-3 B3(q) = -3.259e-005 q^-1 + 6.306e-005 q^-2 - 3.039e-005 q^-3 C(q) = 1 - 1.309 q^-1 + 0.3678 q^-2 Estimated using PEM using SearchMethod = Auto from data set z Loss function 2.41995e-005 and FPE 2.45146e-005 Sampling interval: 1 Mô hình ARMAX bậc na = 1; nb = 2; nc = 2; nk = 1 Discrete-time IDPOLY model: A(q)y(t) = B(q)u(t) + C(q)e(t) A(q) = 1 - 0.9968 q^-1 B1(q) = 0.0001523 q^-1 - 0.0001355 q^-2 B2(q) = -0.003138 q^-1 - 0.004175 q^-2 B3(q) = -5.886e-005 q^-1 + 5.93e-005 q^-2 C(q) = 1 - 0.2903 q^-1 + 0.06295 q^-2 Estimated using PEM using SearchMethod = Auto from data set z Loss function 2.55504e-005 and FPE 2.57806e-005 Sampling interval: 1 Hình 4.13. Nhận dạng áp suất theo phương pháp PEM áp dụng cho mô hình ARMAX Ta thấy, khi sử dụng phương pháp PEM cho mô hình ARMAX, với các bộ tham số khác nhau ta thu được chất lượng mô hình gần giống nhau. Ta tiếp tục sử dụng phương pháp PEM cho các dạng mô hình khác nhau: m7=Pem(Pmhe,'P2DU'); m8=Pem(Pmhe,'P3D'); m9=Pem(Pmhe,'P3DU'); compare(Pmhv,m6,m7,m8,m9); Mô hình hàm truyền liên tục bậc 2 có trễ, có điểm cực phức: Process model with 3 inputs: y = G_1(s)u_1 + ...+ G_3(s)u_3 where Kp G_1(s) = ---------------------- * exp(-Td*s) 1+2*Zeta*Tw*s+(Tw*s)^2 with Kp = 0.0033078 Tw = 58.759 Zeta = 2.4964 Td = 0 Kp G_2(s) = ---------------------- * exp(-Td*s) 1+2*Zeta*Tw*s+(Tw*s)^2 with Kp = -2.8342 Tw = 9.424 Zeta = 3.4828 Td = 18.152 Kp G_3(s) = ---------------------- * exp(-Td*s) 1+2*Zeta*Tw*s+(Tw*s)^2 with Kp = 0.00019435 Tw = 9.4162 Zeta = 8.4444 Td = 14.227 Estimated using PEM using SearchMethod = Auto from data set z Loss function 0.00212907 and FPE 0.00215462 Mô hình quán tính bậc ba có trễ có điểm cực thực Process model with 3 inputs: y = G_1(s)u_1 + ...+ G_3(s)u_3 where Kp G_1(s) = --------------------------- * exp(-Td*s) (1+Tp1*s)(1+Tp2*s)(1+Tp3*s) with Kp = 0.0041342 Tp1 = 443.09 Tp2 = 0.007212 Tp3 = 0.0044328 Td = 0.93008 Kp G_2(s) = --------------------------- * exp(-Td*s) (1+Tp1*s)(1+Tp2*s)(1+Tp3*s) with Kp = 87.393 Tp1 = 37.687 Tp2 = 91.311 Tp3 = 16295 Td = 0 Kp G_3(s) = --------------------------- * exp(-Td*s) (1+Tp1*s)(1+Tp2*s)(1+Tp3*s) with Kp = 0.00037693 Tp1 = 202.81 Tp2 = 1.2197 Tp3 = 0.01133 Td = 28.189 Estimated using PEM using SearchMethod = Auto from data set z Loss function 0.00187927 and FPE 0.00190746 Mô hình quán tính bậc ba có trễ có điểm cực phức Process model with 3 inputs: y = G_1(s)u_1 + ...+ G_3(s)u_3 where Kp G_1(s) = --------------------------------- * exp(-Td*s) (1+2*Zeta*Tw*s+(Tw*s)^2)(1+Tp3*s) with Kp = 0.0034318 Tw = 0.59971 Zeta = 0.001 Tp3 = 315.78 Td = 0 Kp G_2(s) = --------------------------------- * exp(-Td*s) (1+2*Zeta*Tw*s+(Tw*s)^2)(1+Tp3*s) with Kp = -0.62724 Tw = 16.527 Zeta = 0.0099794 Tp3 = 28.385 Td = 0 Kp G_3(s) = --------------------------------- * exp(-Td*s) (1+2*Zeta*Tw*s+(Tw*s)^2)(1+Tp3*s) with Kp = 0.00033122 Tw = 0.3862 Zeta = 0.001 Tp3 = 226.04 Td = 20.715 Estimated using PEM using SearchMethod = Auto from data set z Loss function 0.00148695 and FPE 0.00150925 Kết quả nhận được như sau: Hình 4.14. Nhận dạng áp suất theo phương pháp PEM cho các mô hình hàm truyền đạt liên tục Nhận xét: Mô hình ARMAX theo phương pháp PEM thu được độ fit tương đối đáng tin cậy (75.49%) đồng thời bậc của mô hình lại thấp do ta đã giải quyết được ảnh hưởng của nhiễu đến quá trình. Mô hình hàm truyền liên tục bậc 2, có trễ và có điểm cực phức có độ fit là 73.84%. Mô hình hàm truyền liên tục bậc 3, có trễ có độ fit là 78.47%, mô hình hàm truyền liên tục bậc 3, có trễ và có điểm cực phức có độ fit là 78.04%. Ta thấy kết quả nhận dạng theo phương pháp PEM cho các dạng mô hình khác nhau cho kết quả tương đối chính xác so với phương pháp LSE bên trên. Lí do là: phương pháp sai số dự báo đã khắc phục nhược điểm về tính nhất quán của LSE khi nhiễu không phải là tạp trắng hay khi nhiễu có tương quan với tín hiệu vào. 4.4. Nhận xét chung Với đối tượng nhiệt độ, sử dụng phương pháp PEM ước lượng mô hình ARMAX thu được độ fit cao nhất 92.4% với mô hình là: na = 2; nb = 2; nc = 2; nk = 1. Mô hình trên có bậc 2 thích hợp cho mục đích thiết kế, chỉnh định tham số cho bộ điều khiển. Với đối tượng áp suất, sử dụng phương pháp PEM ước lượng mô hình hàm truyền bậc ba có trễ thu được độ fit cao nhất là 78.47%. CHƯƠNG 5: KẾT LUẬN 5.1. Các phương pháp nhận dạng được sử dụng Phương pháp phân tích Fourier và phân tích phổ Phương pháp phân tích phổ để nhận dạng mô hình đáp ứng tần số cho kết quả nhận dạng không cao mặc dù phương pháp có tính loại bỏ ảnh hưởng của nhiễu tạp trắng trong giá trị tín hiệu đo được. Nguyên nhân là tín hiệu vào là tín hiệu không giàu tần số (có dải tần hẹp ở vùng tần số thấp) do đó ở dải tần số cao mô hình đáp ứng tần số không chính xác. Phương pháp LSE Phương pháp LSE là phương pháp có thuật toán đơn giản thích hợp cho nhận dạng mô hình hồi quy tuyến tính với nhiễu là tạp trắng. Nếu nhiễu là nhiễu màu hay có tương quan với biến vào thì ước lượng cho kết quả không sát thực. Phương pháp PEM Phương pháp PEM cho phép ước lượng tập các mô hình rộng hơn LSE và có khả năng ước lượng mô hình nhiễu. Tuy nhiên chưa đánh giá được tính hội tụ của bộ tham số ước lượng với bộ tham số thực dẫn đến không có phương pháp ước lượng được bậc của mô hình. Phương pháp PEM thu được kết quả tốt hơn cả so với các phương pháp trên, bậc của mô hình cũng thấp nên rất phù hợp cho mục đích của bài toán điều khiển. Do phương pháp PEM đã khắc phục nhược điểm về tính nhất quán của LSE khi nhiễu không phải là tạp trắng hay khi nhiễu có tương quan với tín hiệu vào. 5.2. Các dạng mô hình sử dụng Mô hình đáp ứng tần số cho thông tin về bậc của mô hình tốt nhưng chỉ ước lượng được một mô hình tốt khi có tín hiệu vào giàu tần số. Do phương pháp lấy số liệu là thay đổi giá trị đặt của vòng điều khiển nên ta không lựa chọn được tín hiệu vào có dải tần theo mong muốn dẫn đến chất lượng mô hình đáp ứng tần số thu được là không cao. Mô hình ARX có ưu điểm là thuật toán ước lượng đơn giản (LSE) nhưng chưa mô tả được mô hình của nhiễu do đó không phù hợp với đối tượng có ảnh hưởng bởi nhiễu màu. Đối tượng thực tế ngoài nhiễu đo, luôn chịu tác động của các yếu tố môi trường bất định vào hệ thống. Chính vì vậy với hệ thống chịu tác động mạnh của nhiễu thì mô hình ARX là không thích hợp cho dù dùng các phương pháp ước lượng khác nhau. Mô hình ARMAX, OE, mô hình hàm truyền đạt liên tục, mô hình không gian trạng thái, có xét đến ảnh hưởng của nhiễu vào mô hình (có mô hình nhiễu) nên có thể mô tả chính xác hơn hệ thống thực. 5.3. Chất lượng mô hình nhận được Chất lượng nhận được của mô hình nhiệt độ là từ 88% - 93% ; của mô hình áp suất là 78%- 86%. Kết quả trên có thể chấp nhận được vì tính chất phi tuyến của lò hơi và một số yếu tố ảnh hưởng đến đối tượng mà chúng ta đã bỏ qua: Ngoài nhiễu đo và nhiễu tải chúng ta còn các yếu tố bất định khác bên trong và từ môi trường như nhiệt độ môi trường thay đổi, chất lượng nhiên liệu… Quá trình nhận dạng là nhận dạng vòng kín, trực tiếp nên các nhiễu đo và nhiễu quá trình ảnh hưởng đến hệ thống và có tương quan mạnh với tín hiệu vào và ra của hệ thống. 5.4. Công cụ System Identification Toolbox Công cụ System Identification Toolbox là công cụ rất mạnh, có khả năng hỗ trợ mạnh mẽ cho việc nhận dạng hệ MISO, MIMO. Tuy nhiên chất lượng của mô hình thu được còn phụ thuộc vào kinh nghiệm của người sử dụng trong việc lựa chọn tham số của mô hình và sử dụng các câu lệnh. 5.4. Bài học thu được Trong quá trình thực hiện đồ án, chúng em gặp phải nhiều vấn đề trong việc áp dụng lý thuyết nhận dạng vào đối tượng cụ thể. Qua đó thấy cần phải có sự hiểu biết thêm về vận hành của hệ thống thực để có thể có thể nhận dạng được mô hình có tính trung thực hơn. Bên cạnh những kiến thức đã tìm hiểu để giải quyết bài toán xây dựng mô hình động học bậc thấp cho lò hơi, chúng em đã biết cách triển khai công việc theo từng bước cụ thể, cách đọc và tìm hiểu tài liệu, cách trình bày các vấn đề một cách khoa học…. để hoàn thành đồ án một cách tốt nhất. Về kĩ năng làm việc theo nhóm: Để hoàn thành đồ án một cách tốt đẹp trong khoảng thời gian ngắn, chúng em đã học được cách làm việc, cách phối hợp với nhau để tận dụng được ưu điểm của mỗi người sao cho đạt hiệu quả cao nhất. 5.5. Hướng phát triển và nghiên cứu của đề tài Độ chính xác của các phương pháp nhận dạng mô hình tham số phụ thuộc rất nhiều vào chọn bậc của mô hình. Mô hình đáp ứng tần số có thể cho ta thông tin để ước lượng được bậc của mô hình hỗ trợ cho ước lượng mô hình tham số được chính xác hơn và không phải thử nhiều lần. Tuy nhiên việc ước lượng mô hình đáp ứng tần số theo phân tích phổ chỉ cho mô hình tốt khi tín hiệu vào có dải tần trải rộng trong pham vi đáp ứng pha của hệ thống từ 0 à -2π. Để có thể có được mô thông tin về bậc của mô hình, nhóm sẽ tiếp tục tìm hiêu để xem xét khả năng thực hiện tín hiệu đầu vào thỏa mãn yêu cầu trên. Bên cạnh đấy, như đã nêu trong mục 3.3, trong khuôn khổ đồ án này mô hình thu được chỉ mô tả mối quan hệ mối quan hệ giữa các ngõ vào nhiên liệu FI8201, nước làm mát FI2852 với các ngõ ra áp suất, nhiệt độ hơi quá nhiệt. Do vậy để thiết lập một mô hình đầy đủ cho lò hơi ta còn phải thiết lập thêm mô hình cho các vòng điều khiển mức nước và quá trình cháy trong buồng lửa. Sau khi có được mô hình động học đầy đủ thì ta có thể can thiệp một cách đầy đủ và hợp lý hơn. TÀI LIỆU THAM KHẢO Nguyễn Hữu Quốc Đạt (2009), Luận văn thạc sĩ: Xây dựng mô hình động học cho lò hơi trong Nhà máy đạm Phú Mỹ, Đại học Bách Khoa Hà Nội. Hoàng Minh Sơn (2006), Cơ sở hệ thống điều khiển quá trình, NXB Bách khoa – Hà nội. Lennart Ljung (1999), System Identification, Theory for user, Second Edition, Prentice Hall. Nguyễn Doãn Phước, Phan Xuân Minh (2001), Nhận dạng hệ thống điều khiển, NXB Khoa học và Kĩ thuật. Lennart Ljung (2007), System Identification Toolbox – For use with MATLAB, version 7, MATLAB User Guide, Math Works.