Tối ưu hoá điều độ phát điện

BÁO CÁO ĐỀ TÀI Tối ưu hoá điều độ phát điện LỜI NÓI ĐẦU Sau thời gian 4,5 năm học tập và nghiên cứu tại Trường Đại Học Điện Lực, chúng em đã được các thầy, cô truyền đạt cho những kiến thức cả về lý thuyết và thực hành, để chúng em áp dụng những kiến thức đó vào thực tế và làm quen với công việc của người kỹ sư trong tương lai. Để chuẩn bị cho kỳ tốt nghiệp kết thúc khóa học 2008-2013, em đã được nhận đề tài tốt nghiệp đó là: "Tối ưu hóa điều độ phát điện" do thầy giáo TS Trần Thanh Sơn - giảng viên Bộ môn mạng và hệ thống điện trực tiếp hướng dẫn em làm đồ án này. Đựơc sự giúp đỡ tận tình của thầy cùng các thầy, cô giáo trong khoa, trong trường, với sự lỗ lực của bản thân đến nay em đã hoàn thành đề tài tốt nghiệp của mình. Em xin gửi lời cảm ơn sâu sắc nhất đến thầy cùng toàn thể các thầy cô đã tạo mọi điều kiện giúp em hoàn thành đồ án một cách tốt nhất. Cuối cùng, em xin gửi lời cảm ơn chân thành đến gia đình, bạn bè- những người đã luôn bên cạnh em, tiếp thêm nguồn động lực cho em trong suốt những năm học qua. Đặc biệt, cảm ơn bạn cùng nhóm đề tài với em, chúng em đã cùng nhau nghiên cứu, trau dồi, thảo luận kiến thức để hoàn thành tốt đề tài này. Tuy nhiên, trong quá trình thực hiện đề tài, em còn nhiều bỡ ngỡ, do chưa có kinh nghiệm thực tiễn nên không tránh khỏi những sai sót, nhầm lẫn. Vì vậy, em rất mong nhận được sự góp ý của các thầy, cô giáo để hoàn thành tốt đồ án tốt nghiệp và nhiệm vụ học tập tại trường với kết quả cao. Hà Nội, tháng 1 năm 2013 Sinh viên thực hiện: Nguyễn Quang TùngMỤC LỤC DANH SÁCH BẢNG BIỂU Bảng 2.1: Bảng thông số của các đường dây trong ví dụ 7.9 44 Bảng 2.2: Bảng kết quả modun và góc pha điện áp 45 Bảng 2.3: Bảng kết quả công suất nguồn tại các nút 45 Bảng 2.4: Bảng kết quả công suất tải tại các nút 46 Bảng 2.5: Bảng kết quả sau khi giải tích lưới điện trong ví dụ tài liệu [1] 46 Bảng 2.6: Bảng kết quả các hệ số tổn thất trên lưới 46 Bảng 2.7: Bảng kết quả các hệ sô tổn thất trong ví dụ tài liệu [1] 47 Bảng 3.1: Dữ liệu các nút 49 Bảng 3.2: Dữ liệu máy biến áp 50 Bảng 3.3: Dữ liệu tụ bù 50 Bảng 3.4: Giới hạn công suất tác dụng tại các nút 51 Bảng 3.5: Thông số của các đường dây 51 Bảng 3.6: Kết quả mô đun và góc pha điện áp 53 Bảng 3.7: Công suất phát trước khi tối ưu 54 Bảng 3.8: Hệ số tổn thất trên lưới 54 Bảng 3.9: Bảng tổng kết kết quả so sánh chi phí nhiên liệu 56 DANH SÁCH HÌNH VẼ Hình 1.1: Sơ đồ thuật toán bài toán tối ưu hóa điều độ phát điện 18 Hình 1.2 :Sơ đồ thuật toán chương trình tính toán các hệ số tổn thất B 29 Hình 2.1 : Tổ chức chương trình lập trình trong Matlab 32 Hình 2.2 :Sơ đồ thuật toán chương trình tính toán các hệ số tổn thất B 41 Hình 2.3 : Sơ đồ thuật toán chương trình tối ưu hóa điều độ phát điện 43 Hình 2.4 : Sơ đồ hệ thống điện 5 nút 45 Hình 3.1: Sơ đồ lưới điện 26 nút 49 KÝ HIỆU CÁC CỤM TỪ VIẾT TẮT OPF Optimal power flow HTĐ Hệ thống điện Nút SL Nút cân bằng công suất Nút PV Nút giữ điện áp Nút PQ Nút phụ tải MBA Máy biến áp pu Đơn vị tương đối cb Đơn vị cơ bản GIỚI THIỆU CHUNG Hệ thống điện Việt Nam là một hệ thống điện (HTĐ) hợp nhất. Cùng với sự phát triển của kinh tế, xã hội dẫn đến sự phát triển không ngừng của phụ tải làm cho hệ thống điện vận hành trong tình trạng quá tải. Nguồn điện trong hệ thống điện Việt Nam có rất nhiều loại: thủy điện, nhiệt điện than, nhiệt điện dầu, tuabin khí chu trình đơn và chu trình hỗn hợp, điện diesel,...với các đặc tính vận hành rất khác nhau nên các phương án phân bố công suất cho các nhà máy khác nhau dẫn đến chi phí cho các nhà máy cũng khác nhau đáng kể. Mặt khác các nguồn điện phân bố không đều: Miền Bắc chủ yếu là các nhà máy thủy điện và nhiệt điện than; Miền Nam: bao gồm các nhà máy nhiệt điện dầu, tuabin khí và tuabin khí hỗn hợp; Miền Trung: cho đến nay có ít nhà máy điện. Lưới cung cấp điện còn tồn tại nhiều vấn đề bất cập: có quá nhiều cấp điện áp (110, 66, 35, 22, 15, 10, 6 kV), cũ nát, phi tiêu chuẩn (thiếu bù và các phương tiện điều chỉnh điện áp), tổn thất truyền tải cao, dung lượng bù công suất phản kháng rất thiếu, ... Do sự tăng lên về quy mô và sự phức tạp đó trong hệ thống điện nên vấn đề tối ưu trong quy hoạch và khai thác hệ thống trở nên thiết yếu. Đường dây liên kết với giả thiết có công suất truyền tải khác nhau ảnh hưởng rất nhiều đến chi phí vận hành của toàn hệ thống. Khi đó bài toán tối ưu hóa là rất cần thiết nhằm mục đích đạt được hiệu quả cao trong vận hành. Các bài toán tối ưu là một công cụ hữu hiệu giúp chúng ta có những giải pháp đơn giản nhất để giải quyết một vấn đề dù đơn giản hay phức tạp. Bài toán tối ưu mà bản chất là bài toán giải tìm cực trị của một hàm dưới những ràng buộc nào đó nên có rất nhiều thuật toán giải. Ngày nay, với sự phát triển của khoa học và kĩ thuật, phạm vi ứng dụng của tối ưu hóa ngày càng được mở rộng, các bài toán tối ưu được giải quyết nhanh và chính xác hơn. Ngành Hệ thống điện là một trong những lĩnh vực mà bài toán tối ưu hóa được ứng dụng rất nhiều như: tối ưu hóa chi phí nhiên liệu, tối ưu hóa tổn thất công suất tác dụng, tối ưu hóa tổn thất công suất phản kháng, tối ưu hóa công suất tác dụng của máy phát ở nút cân bằng, tối ưu hóa công suất phản kháng của máy phát ở nút cân bằng, ... Khoảng cách giữa các nhà máy điện đến phụ tải khác nhau nên tổn thất trong quá trình truyền tải là khác nhau và phụ thuộc vào công suất của các nhà máy. Hơn nữa chi phí nhiện liệu cho các nhà máy cũng khác nhau. Do đó vấn đề phân bố công suất tác dụng và công suất phản kháng phát của mỗi nhà máy nhằm cực tiểu hóa hàm chi phí nhiên liệu của toàn bộ hệ thống được đặt ra. Có nghĩa rằng, ta phải tính công suất tác dụng và công suất phản kháng nằm trong giới hạn cho trước mà vẫn đảm bảo nhu cầu phụ tải với chi phí nhiên liệu nhỏ nhất. Bài toán này được gọi là bài toán tối ưu hóa điều độ phát điện. Trong phạm vi nghiên cứu của đề tài này, chúng ta phân tích bài toán tối ưu hóa điều độ công suất tác dụng phát. Bài toán này phục vụ tính phân bố công suất cho các nhà máy sao cho tổng chi phí nhiên liệu của hệ thống đạt giá trị nhỏ nhất với điều kiện có tính đến tổn thất công suất tác dụng và giới hạn công suất tác dụng phát. Đề tài tốt nghiệp gồm 3 chương như sau: Chương 1: Tối ưu hóa điều độ phát điện. Chương này giới thiệu bài toán tối ưu hóa điều độ phát điện trong ngành HTĐ. Mô hình toán học của bài toán được đưa ra sau đó bài toán được giải bằng phương pháp Lagrange. Một ví dụ với các số liệu cụ thể để có sự hình dung rõ hơn về vấn đề được giới thiệu. Phần cuối của chương giới thiệu phương pháp tính toán các hệ số tổn thất để phục vụ cho bài toán tối ưu hóa điều độ phát điện. Chương 2: Chương trình tính toán tối ưu hóa điều độ phát điện. Chương này gồm 2 phần: phần đầu chương giới thiệu về lập trình các thuật toán ở trên trong Matlab, phần thứ 2 đưa ra một hệ thống điện đơn giản để kiểm tra thuật toán và chương trình đã lập trình. Chương 3: Ứng dụng chương trình. Chương này sử dụng chương trình để tính toán cho trường hợp phức tạp hơn với bài toán quy mô lớn hơn. Và cuối cùng là một số những kết luận chung. Tối ưu hóa điều độ phát điện Mô hình bài toán điều độ phát điện Hệ thống điện bao gồm nhiều nhà máy điện và các phụ tải được nối với nhau thông qua các đường dây truyền tải điện. Ở mỗi chế độ làm việc, các phụ tải sẽ yêu cầu cung cấp một lượng công suất nhất định. Việc phân chia lượng công suất này cho các nhà máy điện trong hệ thống để đạt được chi phí sản xuất nhỏ nhất là một bài toán lớn, vì: Mỗi nhà máy điện có một hàm chi phí nhiên liệu riêng và giới hạn công suất phát khác nhau . Tổn thất công suất từ các nhà máy đến các phụ tải cũng rất khác nhau do khoảng cách và công suất tiêu thụ của các phụ tải là khác nhau. Trong phần này chúng ta sẽ xét mô hình toán học của bài toán tối ưu hóa điều độ phát điện khi tính đến tổn thất công suất tác dụng và giới hạn công suất phát. Khi truyền tải một lượng công suất lớn mà khoảng cách truyền tải ngắn thì tổn thất truyền tải có thể được bỏ qua và việc tối ưu hóa công suất phát đạt được với tất cả các nhà máy được thể hiện ở phương trình chi phí sản xuất. Tuy nhiên, trong một lưới điện liên kết lớn nơi mà năng lượng điện được truyền tải trên một khoảng cách dài với khu vực mật độ tải thấp thì tổn thất truyền tải là một yếu tố quyết định và có ảnh hưởng đến việc tối ưu hóa điều độ phát điện của các nhà máy. Để giải quyết bài toán đặt ra ta coi tổng tổn thất truyền tải như một hàm bậc hai của các công suất phát của các máy phát điện. Dạng phương trình đơn giản nhất là: ( 0.1) Trong đó ng là tổng số máy phát trong hệ thống xét. Khai triển công thức thành các thành phần tuyến tính và hằng số gọi là công thức tổn thất Kron được trình bày ở mục 1.4 có dạng: ( 0.2) Trong đó Bij: hệ số tổn thất giữa các nút i và nút j còn gọi là hệ số B , có thể coi là không đổi trong quá trình tính toán B0i: hệ số tổn thất ứng với nút thứ i B00: hệ số tổn thất cố định ng – số máy phát phân bố tối ưu công suất trong tổng số n máy phát Pi , Pj - công suất ứng với nút thứ i và nút thứ j. Có nhiều phương pháp để xây dựng phương trình tổn thất. Một phương pháp tính toán các hệ số B được trình bày trong mục 1.4. Đặc tính chi phí nhiên liệu sản xuất điện trong 1 giờ của nhà máy nhiệt điện thứ i có thể biểu thị dưới dạng hàm bậc 2 theo công suất tác dụng của nguồn phát như sau: ( 0.3) Trong đó là các hệ số hồi quy được tính toán từ đặc tính thực nghiệm suất tiêu hao nhiên liệu của nguồn thứ i Pi là công suất phát của nguồn thứ i. Bài toán điều độ phát điện là tìm các giá trị Pi vừa thỏa mãn các điều kiện rằng buộc vừa cực tiểu hóa được hàm chi phí nhiên liệu Ct của hệ thống. Các điều kiện rằng buộc là tổng công suất phát cân bằng với với tổng công suất tiêu thụ của phụ tải cộng với tổn thất: ( 0.4) và điều kiện về giới hạn công suất phát của mỗi nhà máy trong hệ thống xét: ( 0.5) Trong đó: Pi(max); Pi(min) là các giới hạn công suất của nhà máy thứ i. Phương pháp hệ số Lagrange Một trong những phương pháp thông dụng và phổ biến để giải bài toán tối ưu hóa có xét đến các điều kiện rằng buộc là phương pháp hệ số Lagrange. Ta xét bài toán được ra như sau: Xác định giá trị nhỏ nhất của hàm sau: ( 0.6) với điều kiện các biến x1, x2, x3, . . . xn phải thỏa mãn các điều kiện rằng buộc dạng phương trình: ( 0.7) hoặc các rằng buộc dạng bất phương trình: ( 0.8) Bài toán trên được giải bằng phương pháp hệ số Lagrance như sau: Ta sẽ lập hàm chi phí mới có tính đến các điều kiện rằng buộc. Hàm Lagrange được thành lập như sau: ( 0.9) Trong đó và là các hằng số Tính đạo hàm riêng của hàm Lagrance ở trên theo các biến và cho triệt tiêu ta nhận được hệ phương trình sau: ( 0.10) Giải hệ phương trình trên ta nhận được các giá trị x1, x2 , x3…, xn và () Thay lại hàm ban đầu các giá trị vừa tìm được rồi kết luận. Chú ý rằng phương trình là rằng buộc đơn giản ban đầu. Giả sử (x1,x2,…, xn) là giá trị nhỏ nhất thỏa mãn. Bất phương trình rằng buôc không xảy ra nếu bất phương trình tại điểm (x1,x2,…, xn) có =0. Mặt khác, khi phương trình cân bằng được xác lập, rằng buộc thỏa mãn ở điểm (x1,x2,…, xn) đó (ví dụ nếu rằng buộc và > 0) . Đây chính là điều kiện Kuhn-Tucker. Áp dụng cho bài toán tối ưu hóa điều độ phát điện Trong phần này ta xét ứng dụng phương pháp hệ số Lagrange để giải bài toán tối ưu hóa điều độ ở phần 1.1. Thêm các điều kiện rằng buộc vào hàm mục tiêu trong phương trình (1.3) ta nhận được hàm Lagrance có dạng: ( 0.11) Các ràng buộc được hiểu theo nghĩa là: ( 0.12) Giá trị nhỏ nhất của các hàm ràng buộc này được tìm thấy tại điểm mà ở đó đạo hàm riêng của hàm theo các biến của nó bằng 0. ( 0.13) ( 0.14) ( 0.15) ( 0.16) Trong công thức (1.15) và (1.16) Pi không được phép vượt giới hạn của nó và khi Pi tiến dần đến giới hạn thì: ( 0.17) và hàm Kuhn-Tucker trở thành một hàm Lagrange. Điều kiện thứ nhất đưa bởi (1.13) cho ta kết quả sau: ( 0.18) Ta có: ( 0.19) Suy ra: ( 0.20) Vì vậy, điều kiện cho tối ưu hóa điều độ là: ( 0.21) Trong đó được gọi độ gia tăng tổn thất truyền tải. Điều kiện thứ 2 được đưa bởi (1.14), kết quả là: ( 0.22) Phương trình (1.22) chính là phương trình điều kiện rằng buộc của bài toán Biến đổi phương trình (1.21) dưới dạng như sau: ( 0.23) Hoặc: ( 0.24) Trong đó Li là hệ số phạt của máy phát thứ i, được cho bởi công thức: ( 0.25) Phương trình (1.24) chỉ ra rằng chi phí nhỏ nhất đạt được khi độ gia tăng chi phí của mỗi nhà máy theo cấp số nhân bởi hệ số phạt với các nhà máy thì như nhau. Độ gia tăng chi phí sản xuất được đưa bởi công thức: ( 0.26) và gia tăng tổn thất công suất truyền tải thu được từ công thức tính tổn thất của Kron chính là lợi nhuận (sản lượng): ( 0.27) Thay thế kết quả vào công thức tính độ gia tăng chi phí sản xuất (1.26) và gia tăng tổn thất công suất truyền tải (1.21) ta có kết quả sau: ( 0.28) Hoặc: ( 0.29) Áp dụng (1.29) cho tất cả các nhà máy, ta nhận được 1 hệ phương trình tuyến tính ở dạng ma trận như sau: ( 0.30) Hay dưới dạng ngắn gọn: ( 0.31) Trong đó Để tìm tối ưu hóa điều độ phát điện ta chọn một giá trị ước lượng ban đầu của λ(1). Công thức (1.29) để tính Pi ở vòng lặp thứ k được xác định như sau: ( 0.32) Thay thế Pi từ (1.32) vào (1.22) ta nhận được: ( 0.33) Hoặc: ( 0.34) Khai triển vế bên phải của phương trình trên theo khai triển Taylor tại một điểm λ(k) và bỏ qua các thành phần bậc cao thì được kết quả là: ( 0.35) Hoặc: ( 0.36) Trong đó: ( 0.37) ( 0.38) Ta được giá trị mới của λ là: ( 0.39) Quá trình lặp cứ tiếp tục cho đến khi ∆P(k) nhỏ hơn một giá trị xác định. Nếu công thức tổn thất được lấy xấp xỉ theo công thức sau: ( 0.40) Khi đó: Bij=0, B00=0 và kết quả của phương trình cho bởi (1.32) được xác định đơn giản như sau: ( 0.41) Và cho bởi (1.37) là: ( 0.42) Dưới đây là sơ đồ thuật toán của bài toán trên Hình 0.1: Sơ đồ thuật toán bài toán tối ưu hóa điều độ phát điện Để làm rõ hơn ta tính toán cho 1 bài toán được trích từ trang 284 – tài liệu [1]. Bài toán: Chi phí nhiên liệu của 3 nhà máy nhiệt điện của một hệ thống điện như sau: Công suất phát của các nhà máy được giới hạn sau: Giả sử tổn thất điện năng được cho bởi biểu thức đơn giản: Lấy Scb=100 MVA. Xác định điều độ tối ưu của các nhà máy khi tổng phụ tải hệ thống là 150MW. MW MW Cho các giải pháp số sử dụng phương pháp Gradient, giả sử giá trị ban đầu từ phương trình phối hợp (1.41), ta tính được P1, P2, P3 là: Tổn thất công suất Với PD = 150 MW. Sai lệch theo công thức (7.68) là: Từ (1.42) ta suy ra Từ (1.36) ta tính được: Khi đó, giá trị mới của là Tiếp tục như trên với bước 2 ra có Tổn thất công suất: Với PD = 150 MW. Sai lệch là: Từ (1.42) Từ (1.36) Khi đó, giá trị mới của là Lặp với bước 3 ta có Tổn thất công suất: Với PD= 150 MW. Sai lệch theo công thức (7.68) là Từ (1.42) Từ (1.36) Khi đó, giá trị mới của là Với , đó là điều kiện rằng buộc nhỏ được đáp ứng qua 4 bước lặp, và điều độ tối ưu với Tổn thất công suất Và tổng chi phí nhiên liệu là: Ct = 200 + 7.(35,0907) + 0,008(35,0907)2 + 180 + 6,3(64,1317) + + 0,009(64,1317)2 + 140 + 6,8(52,4767) + 0,007(52,4767)2 = 1592,65 $/h Phương pháp tính hệ số tổn thất Một trong những thông số phục vụ bài toán tính toán tối ưu điều độ công suất phát của nhà máy điện là xác định được hàm tổn thất công suất trên hệ thống theo các biến là công suất phát của các nhà máy. Có một vài phương pháp xây dựng công thức tính tổn thất. Một phương pháp được phát triển bởi Kron và bổ sung bởi Kirchmayer là phương pháp hệ số tổn thất B. Nội dung phương pháp như sau: Công suất nút bơm vào tại một nút i, kí hiêụ là Si được tính theo điện áp dây và dòng nút Ii như sau: ( 0.43) Tổng tổn thất công suất trên hệ thống điện xét được tính bằng tổng công suất nút của tất cả các nút trong hệ thống: ( 0.44) PL và QL là tổn thất công suất tác dụng và công suất phản kháng của hệ thống. Ubus là véc tơ cột của điện áp nút và Ibus là véc tơ cột của dòng điện bơm vào nút. Mặt khác công thức tính dòng điện bơm vào nút : ( 0.45) Ybus là ma trận tổng dẫn nút của hệ thống điện. Từ công thức trên ta suy ra: ( 0.46) Ma trận Zbus được gọi là ma trận tổng trở nút của hệ thống Rút Vbus từ (1.46) thay vào (1.44) ta được kết quả sau: ( 0.47) Zbus là ma trận đối xứng, do đó: và tổng tổn thất công suất là: ( 0.48) Công thức tính (1.48) có thể được viết dưới dạng sau: ( 0.49) Do ma trận tổng trở nút của hệ thống đối xứng nên đẳng thức trên được viết lại là : ( 0.50) Tách phần thực và phần ảo, như vậy tổn thất công suất là: ( 0.51) Trong đó: Rij và Xij là các thông số của ma trận tổng trở nút. Như vậy từ thì phương trình tính tổn thất công suất tác dụng đổi lại là: ( 0.52) Viết dưới dạng ma trận là: ( 0.53) Trong đó : Rbus là phần thực của ma trận tổng trở nút Để đưa ra công thức tính tổng tổn thất công suất theo công suất phát của các nhà máy điện, ta đi xác định tổng dòng tải ID từ từng dòng tải ILi: ( 0.54) Trong đó: n là số nút, là dòng tải tại nút i. ( 0.55) Để đơn giản hóa ta đặt ( 0.56) Coi nút 1 là nút slack bus, khai triển công thức ta có: Thay tất cả các chữ ( 0.57) Trong đó : tổng số nút có máy phát và ( 0.58) Đặt: ( 0.59) Suy ra: ( 0.60) ( 0.61) Đặt : ( 0.62) Suy ra : ( 0.63) Dòng nút : ( 0.64) Với :dòng nút ; : dòng nguồn; : dòng tải Viết dưới dạng ma trận ta có: ( 0.65) Hay: ( 0.66) ( 0.67) ( 0.68) ( 0.69) Rút Ibus thay vào ta có: ( 0.70) Nếu là công suất phức nguồn tại nút i thì dòng điện nguồn tương ứng của các nguồn là: ( 0.71) Hay: ( 0.72) Viết dưới dạng ma trận: ( 0.73) Hay : ( 0.74) Thay thế ta có : ( 0.75) ( 0.76) Suy ra: ( 0.77) Các thông số trong ma trận H là phức hợp. Phần thực được dùng để tính tổn thất công suất. Đó gọi là ma trận Hermitian, =*. Vì vậy, phần thực của được xác định : ( 0.78) Các thành phần của ma trận trên như sau : ( 0.79) Thay thế vào công thức PL ta được: ( 0.80) Hay: ( 0.81) Đổi các hệ số về dạng đơn vị có tên theo công thức sau: ( 0.82) Với SB là công suất cơ bản có đơn vị là MVA Hình 1.2 :Sơ đồ thuật toán chương trình tính toán các hệ số tổn thất B Chương trình tính tối ưu hóa điều độ phát điện Giới thiệu về chương trình Matlab Matlab là một trong những hệ thống được soạn thảo hết sức chu đáo, mạch lạc với việc áp dụng nhiều thuật toán phong phú, đặc biệt như chính tên của nó – “Matric Laboratory’’, thuật toán ma trận được coi là hạt nhân quan trọng trong lập trình Matlab. Cú pháp ngôn ngữ của chương trình Matlab được thiết lập tinh tế đến mức người sử dụng không có cảm giác là đang tiếp xúc với các phép tính phức tạp của ma trận. Ma trận được áp dụng rộng rãi trong các bài toán kỹ thuật điện, điện tử, cơ học như: toán tử phức, ma trận, véc tơ, đa thức, xử lý tín hiệu... Matlab có thể làm việc ở hai chế độ: tương tác và lập trình. Trong tương tác, Matlab thực hiện từng lệnh được gõ trong cửa sổ lệnh sau dấu nhắc lệnh (>>). Trong chế độ lập trình, một tập lệnh được soạn thảo và ghi lại thành một tệp đuôi .m (m-file). Để gọi chương trình chỉ cần gọi tên m- file trong cửa sổ lệnh. Hệ thống Matlab cho phép người sử dụng giải các bài toán chuyên sâu như: bài toán tối ưu, hệ thống điều khiển, mạng noron, logic mở, phân tích tài chính, mô phỏng dạng sóng...với độ chính xác cao và tốc độ nhanh. Ngôn ngữ của hệ thống Matlab phong phú hơn bất kỳ một ngôn ngữ lập trình nào. Nó gần như sử dụng tất cả các phương tiện lập trình đã biết, kể cả mô phỏng (phương tiện Simulink), đồ họa, hoạt hình...Hệ thống Matlab còn chứa một lượng khổng lồ các phép tính, thuật giải và hàm cho phép giải rất nhiều bài toán phức tạp như: lấy nghịch đảo ma trận, đạo hàm, vi phân...Ngôn ngữ đầu vào của Matlab gần giống như Basic ( hỗn hợp của Fortran và Pascal). Vì vậy, rất quen thuộc và dễ dàng cho người sử dụng. Thêm vào đó, hệ thống còn có khả năng hiệu chỉnh chương trình, khi người sử dụng không đúng thì máy sẽ báo lỗi giúp cho việc tìm lỗi dễ dàng hơn. Do đó, Matlab là một phần mềm hết sức thuận tiện và hiệu quả. Là một hệ thống mở, Matlab kết hợp rất nhiều phương pháp tính mà có thể áp dụng thuận tiện cho bất kì người sử dụng nào. Matlab có đặc điểm linh hoạt và dễ thích nghi. Vì vậy, cho dù là người mới bắt đầu hay là một chuyên gia sành sỏi, vẫn có thể sử dụng Matlab một cách thoải mái và cuốn hút. Matlab rất ưu ái với người sử dụng, bạn có thể tạo ra chương trình chuyên môn của mình mà Matlab sẵn sang đáp ứng. Matlab cho phép tiếp cận và áp dụng dễ dàng các hàm có sẵn để giải các bài toán cần thiết và đồng thời có thể sang tạo ra các m.file mà khi được lưu giữ với thủ tục phù hợp, Matlab coi như là hàm của chính hệ thống, điều đó cho phép mở rộng khả năng vô hạn của Matlab. Đó cũng chính là ưu điểm nổi bật của Matlab mà không chương trình nào có được. Tổ chức chương trình Dựa trên các thuật toán đã giới thiệu ở chương 1, một chương trình tính toán tối ưu hóa điều độ phát điện đã được mã hóa trong môi trường Matlab. Chương trình được tổ chức như trên hình 2.1. Hình 0.1 : Tổ chức chương trình lập trình trong Matlab Để thực hiện tính toán, chúng ta phải thực hiện 2 bước sau: Thứ nhất là nhập dữ liệu vào các file. Các file đã được dịnh dạng sẵn, người dùng phải nhập các dữ liệu theo các định dạng này. Nếu không chương trình sẽ báo lỗi hoặc cho kết quả sai. Trong tính toán tối ưu hóa điều độ phát điện, chúng ta cần nhập các loại dữ liệu như sau: Các thông số của nút Các thông số của đường dây và máy biến áp Các hệ số của hàm chi phí nhiên liệu Thứ 2 là cho chương trình chính chạy. Trong quá trình chạy, chương trình chính sẽ thức hiện các chức năng sau: Gọi các chương trình con để đọc dữ liệu nút, dữ liệu nhánh và các hệ số của hàm chi phí nhiên liệu rồi gán cho các biến tương ứng. Thực hiện tính toán chế độ xác lập lưới điện bằng phương pháp Newton-Raphson hoặc Gauss-Seidel Gọi chương trình con tính các hệ số tổn thất Gọi chương trình tính tối ưu hóa điều độ phát điện. Bước thứ 2 này được tổ chức trong testMainFunc.m, nội dung file như sau: clc; clear all; % Vi du 7.9 --- % Doc du lieu nut [tongSoNut, nut, Udm, Ucb, Scb, code, modunU, phaU, Sng, Stai, Qmin, Qmax, Ys] ... = readBusData('filetest_BusData.xls') maxNut = max(nut); % Doc du lieu duong day [tongSoDuongDay, soDuongDay, id, nutDau, nutCuoi, UcbDZ, ScbDZ, chieuDai, ... Z0, Y0ngang, heSoPhanAp]= readLineData(maxNut, 'filetest_LineData.xls'); % Doc du lieu ham chi phi va cong suat gioi han [alpha,beta,gamma,Pmax,Pmin] ... = readData('filetest_Data.xls'); % Tinh ma tran tong dan Ybus [Ybus] = ybus(tongSoDuongDay, soDuongDay, id, nutDau, nutCuoi, UcbDZ, ... ScbDZ, chieuDai, Z0, Y0ngang, heSoPhanAp); %================================================= % Tinh Toan GS % Giai tich luoi dien bang GS epsilon = 1e-10; MaxNbrIte = 1000; accel = 1.8; [magU, theta, nbrIte, doLechCongSuat] = plfNR(tongSoNut, nut, Udm, Ucb, Scb, code, modunU, phaU, Sng, Stai, Qmin, Qmax, Ybus, ... epsilon, MaxNbrIte); U = magU(nbrIte,:).*exp(1i*theta(nbrIte,:)); % Tinh cong suat nut [Snut] = calculateSbus(tongSoNut, code, Sng, Stai, Scb, U... , Ybus); %Tinh ma tran ton that bloss [B, B0, B00, nbrGen ] = Bloss(Snut, Stai, Ybus, code, U); % Tinh toi uu hoa phan bo cong suat cho cac nha may epsilon1 = 0.5e-5; [P_matrix, c_matrix,lamda,P_D,P_L,PL_matrix,Pk,C] = mainfunc(epsilon1, tongSoNut, nbrGen, Snut, Stai, code, alpha, beta, gamma, Pmin, Pmax, B, B0, B00); Các thông số, kí hiệu được sử dụng trong các hàm của chương trình: Nhập các dữ liệu chương trình trong 3 file: filetest_BusData.xls Số nút Giá trị: Udm: điện áp định mức của mạng điện(kV) Ucb: điện áp cơ bản dạng có tên(kV) Scb: công suất cơ bản dạng có tên(MVA) Loại nút (code) : dựa vào code để phân loại các nút (nút SL code=0, nút PV code=1, nút PQ code=2) Mo dun U, góc pha U : giá trị điện áp và góc pha điện áp ở dạng đơn vị tương đối(pu) Giá trị: Pnguon, Qnguon : các công suất của nguồn(MW, MVAr) Ptai, Qtai : các công suất của tải (MW, MVAr) Qmax, Qmin: các giới hạn công suất phản kháng tại nút (MVAr) Gshunt, Bshunt : thành phần phức và ảo của tổng dẫn nút filetest_LineData.xls Thứ tự nhánh Dòng ID biểu thị số nhánh nối vào 1 nút. Nút đầu: tên nút đầu đường dây Nút cuối: tên nút cuối đường dây Giá trị: Ucb: điện áp cơ bản(kV) Scb: công suất cơ bản (MVA) Chiều dài đường dây: (km) Giá trị điện trở R0, điện kháng X0, điện dẫn B0, dung dẫn G0: các thông số cơ bản của đường dây Hệ số phân áp filetest_Data.xls Nút: tên nút Giá trị: Pmax, Pmin- các giới hạn công suất tác dụng tại nút (MW) Giá trị: alpha, beta, gama: các hệ số của hàm chi phí Các chương trình lập trình Matlab như sau: readBusData.m Những chỗ này bạn thêm dòng như ở dưới cho tôi [tongSoNut, nut, Udm, Ucb, Scb, code, modunU, phaU, Sng, Stai, Qmin, Qmax, Ys] = readBusData(filename) Đọc dữ liệu từ file filetest_BusData.xls và gán cho các biến: tongSoNut, nut, Udm, Ucb, Scb, code, modunU, phaU, Sng, Stai, Qmin, Qmax, Ys. Trong đó: Input: filename: Tên file (dạng excel) và đường dẫn chứa dữ liệu nút Output: tongSoNut: Tổng số nút của lưới điện mô phỏng nut: Bảng đánh số nút của lưới điện mô phỏng Udm: Bảng chứa điện áp định mức của các nút Ucb: Bảng chứa điện áp cơ bản của các nút Scb: Bảng chứa công suất cơ bản của các nút code: Bảng đánh code(loại nút) của nút 0 : nút SL 1 : nút PV 2 : nút PQ modunU: Bảng chứa mô đun điện áp của nút phaU: Bảng chứa góc pha của điện áp Sng: Bảng chứa công suất của nguồn tại nút Stai: Bảng chứa công suất của tải tại nút Qmin: Bảng chứa giới hạn dưới của công suất phản kháng Qmax: Bảng chứa giới hạn trên của công suất phản kháng Ys: Bảng chứa tổng dẫn Shunt của nút readLineData.m [tongSoDuongDay, soDuongDay, id, nutDau, nutCuoi, UcbDZ, ScbDZ,chieuDai,Z0,Y0ngang, heSoPhanAp] = readLineData(maxNut, filename) Đọc dữ liệu từ file filetest_LineData.xls và gán cho các biến: tongSoDuongDay, soDuongDay, id, nutDau, nutCuoi, UcbDZ, ScbDZ, chieuDai, Z0, Y0ngang, heSoPhanAp. Trong đó: Input: maxNut: Số nút lớn nhất filename: Tên file (dạng excel) và đường dẫn chứa dữ liệu đường dây Output: tongSoDuongDay: Tổng số đường dây của lưới điện mô phỏng soDuongdayk: Bảng ghi đánh số đường dây id: Bảng ghi id cả đường dây nutDau: Bảng ghi nút đầu của đường dây nutCuoi: Bảng ghi nút cuối của đường dây UCbDZ: Bảng ghi Ucb của đường dây ScbDZ: Bảng ghi Scb của đường dây chieuDai: Bảng ghi chiều dài cơ bản của đường dây Z0: Bảng ghi tổng trở dọc đơn vị của đường dây Y0ngang: Bảng ghi tổng dẫn ngang cơ bản của đường dây heSoPhanAp: Hệ số phân áp, đối với đường dây có giá trị bằng 1, đối với MBA thì bằng tỷ số đầu phân áp readData.m [alpha,beta,gamma,Pmax,Pmin] = readData(filename) Đọc dữ liệu từ file filetest_Data.xls và gán cho các biến: alpha, beta, gamma, Pmax, Pmin. Trong đó: Input: filename: Ten file (dang excel) va duong dan chua du lieu nut Output: alpha, beta, gamma: các hệ số của hàm chi phí Pmax, Pmin: các giới hạn công suất tác dụng tại nút (MW) ybus.m [Y] = ybus(tongSoDuongDay, soDuongDay, id, nutDau, nutCuoi, UcbDZ, ScbDZ, chieuDai, Z0, Y0ngang, heSoPhanAp) Chương trình tính toán ma trận tổng dẫn từ ma trận tổng trở và được gán cho các biến: tongSoDuongDay, soDuongDay, id, nutDau, nutCuoi, UcbDZ, ScbDZ, chieuDai, Z0, Y0ngang, heSoPhanAp. Trong đó: Input: tongSoDuongDay: Tổng số đường dây của lưới điện mô phỏng soDuongdayk: Bảng ghi đánh số đường dây id: Bảng ghi ID đường dây nutDau: Bảng ghi nút đầu của đường dây nutCuoi: Bảng ghi nút cuối của đường dây UCbDZ: Bảng ghi Ucb của đường dây ScbDZ: Bảng ghi Scb của đường dây chieuDai: Bảng ghi chiều dài của đường dây Z0: Bảng ghi tổng trở dọc đơn vị của đường dây Y0ngang: Bảng ghi tổng dẫn ngang đơn vị của đường dây Output: Y: ma trận tổng dẫn plfNR.m [magU, theta, tongSoBuocLap, doLechCongSuat] = plfNR(n, nut, Udm, Ucb, Scb, code, modunU, phaU, Sng, Stai, Qmin, Qmax, Ybus, ... epsilon, soBuocLapLonNhat) Chương trình giải tích lưới điện bằng phương pháp Newton- Raphson và được gán cho các biến magU, theta, tongSoBuocLap, doLechCongSuat. Trong đó: Input: n: Tổng số nút của lưới điện mô phỏng nut: Bảng đánh số nút của lưới điện mô phỏng Udm: Bảng chứa điệ áp định mức của các nút Ucb: Bảng chứa điệ áp cơ bản của các nút Scb:Bảng chứa công suất cơ bản của các nút code: Bảng đánh code (loại nút) của nút 0 : nút SL 1 : nút PV 2 : nút PQ modunU: Bảng chứa mô đun điện áp của nút phaU: Bảng chứa góc pha điện áp Sng: Bảng chứa công suất của nguồn tại nút Stai: Bảng chứa công suất của tải tại nút Qmin : Bảng chứa giới hạn dưới của công suất phản kháng Qmax: Bảng chứa giới hạn trên của công suất phản kháng Ybus: ma trận tổng dẫn nút epsilon: sai số cho phép soBuocLapLonNhat: số bước lặp lớn nhất Output: magU(k,i): giá trị modun điện áp của nút i tại bước lặp thứ k theta(k,i): góc pha điện áp của nút i tại bước lặp thứ k tongSoBuocLap: Tổng số bước lặp doLechCongSuat: Độ lệch công suất lớn nhất sau khi lặp calculateSnut.m [Snut] = calculateSbus(tongSoNut, code, Sng, Stai, Scb, U, Ybus) Chương trình tính toán công suất tại các nút và được gán cho biến Snut. Trong đó: Input: tongSoNut: tổng số nút code: Bảng đánh code (loại nút) của nút 0 : nút SL 1 : nút PV 2 : nút PQ Sng: Bảng chứa công suất của nguồn tại nút Stai: Bảng chứa công suất của tải tại nút Scb:Bảng chứa công suất cơ bản của các nút U: Bẳng chứa điện áp tại nút Ybus: Ma trận tổng dẫn nút Output: Snut: Công suất nút dạng có tên bloss.m [B B0 B00 nbrGen] = Bloss(Snut, Stai, Ybus, code, Unut) Chương trình tính toán hệ số tổn thất của lưới điện. Input: Snut : Công suất nút dạng có tên Stai : Công suất tải dạng có tên Ybus : ma trận tổng dẫn nút code : Bảng đánh code phân loại nút Unut :Điện áp nút Output: B B0 B00: các hệ số tổn thất nbrGen: Tổng số nút nguồn Hình 0.2 :Sơ đồ thuật toán chương trình tính toán các hệ số tổn thất B mainfunc.m [P_matrix,c_matrix,lamda,P_D,P_L,PL_matrix,Pk,C] = mainfunc(epsilon,tongSoNut, nbrGen, Snut, Stai, code, alpha, beta, gamma, ...Pmin, Pmax, B, B0, B00) Chương trình chính tính toán tối ưu hóa điều độ lưới điện và được gán các biến P_matrix, c_matrix,lamda,P_D,P_L,PL_matrix,Pk,C. Trong đó: Input: epsilon:sai số cho phép tongSoNut:tổng số nút nbrGen:số nút nguồn Snut, Stai:công suất tại nút, tại nút tải code :bảng code phân loại nút alpha, beta, gamma :các giá trị hệ số của hàm chi phí Pmin, Pmax :giới hạn dưới và trên của công suất tác dụng B, B0, B00 : các hệ số tỏn thất tính được từ chương trình bloss.m Output: P_matrix: ma trận lưu giữ giá trị P cho tất cả các lần c_matrix:ma trận chi phí của từng nhà máy lamda: chi phí trên một đơn vị điện năng P_D :tổng côg suất tác dụng của các nút tải P_L : lượng tổn thất công suất tác dụng Pk: giá trị P đã tối ưu hóa C : tổng chi phí của toàn hệ thống Hình 0.3 : Sơ đồ thuật toán chương trình tối ưu hóa điều độ phát điện Kiểm tra chương trình Để kiểm tra chương trình ta sử dụng lưới điện được trích từ ví dụ 7.9 – trang 295-296 – tài liệu [1]. Sơ đồ lưới điên như hình vẽ 2.4. Lưới điệnĐường dây gồm 5 nút với 3 máy phát được đặt tại các nút 1, 2, 3. Nút 1 là nút cân bằng công suất, có điện áp pu. Nút 2 và 3 là nút giữ điện áp có mMô đun điện áp và công suất tác dụng ở nút 2 và 3 lần lượt tương ứng là: 1,045 pu, 40MW và 1,03 pu, 30MW. Công suất cơ bản lấy bằng 100 MVA. Thông số của các đường dây cho trong bảng sau: Bảng 0.1: Bảng thông số của các đường dây trong ví dụ 7.9 Nhánh Nút đầu Nút cuối R () X () B() 1 1 2 0,02 0,06 0,06 2 1 3 0,08 0,24 0,05 3 2 3 0,06 0,18 0,04 4 2 4 0,06 0,18 0,04 5 2 5 0,04 0,12 0,03 6 3 4 0,01 0,03 0,02 7 4 5 0,08 0,24 0,05 Giới hạn công suất tác dụng tại các nút 1, 2 và 3 của mỗi máy là: Giới hạn công suất phản kháng tại các nút 1, 2 và 3 là: Hàm chi phí cho 3 máy phát tương ứng như sau: Tính toán các dòng công suất và sử dụng chương trình bloss để tìm các hệ số tổn thất ở dạng đơn vị tương đối. Tìm kết quả của bài toán tối ưu hóa điều độ phát điện. Hình 0.4 : Sơ đồ hệ thống điện 5 nút Ta thực hiện mô phỏng bằng chương trình đã viết và so sánh với kết quả trong tài liệu. Kết quả thu được sau khi giải tích lưới điện bằng phương pháp Newton-Raphson Bảng 0.2: Bảng kết quả modun và góc pha điện áp Nút ModunMô đun Góc pha (0) đơn vị là gì? 1 1,06 0 2 1,045 -1,782435777649380 3 1,03 -2,663957385290456 4 1,018631056611613 -3,243097920284439 5 0,990100467545556 -4,405014491141186 Bảng 0.3: Bảng kết quả công suất nguồn tại các nút Nút P (MW) Q (Mvar) 1 83,051200146755550 7,271372039145252 2 40 41,810740606005000 3 30 24,147606040691870 4 0 0 5 0 0 Bảng 0.4: Bảng kết quả công suất tải tại các nút Nút P (MW) Q (Mvar) 1 0 0 2 20 10 3 20 15 4 50 30 5 60 40 Đáp án của ví dụ trong tài liệu [1] như sau: Bảng 0.5: Bảng kết quả sau khi giải tích lưới điện trong ví dụ tài liệu [1] Kết quả thu được sau khi tính toán các hệ số tổn thất của lưới điện Bảng 0.6: Bảng kết quả các hệ số tổn thất trên lưới B = 0,021800675639716 0,009279886585469 24,147606040691870 0,009279886585469 0,022825143554939 0,001661695963918 0,002802801931712 0,001661695963918 0,017917992374288 B0 = 0,000333555810561 0.003069049858572 0.001507022700345 B00 = 3,052281917197353e-04 Đáp án của ví dụ trong tài liệu [1] như sau: Bảng 0.7: Bảng kết quả các hệ sô tổn thất trong ví dụ tài liệu [1] Kết quả thu được sau khi chạy chương trình chính tính toán tối ưu hóa phát điện Công suất phát tối ưu của các nhà máy P1 = 33,455812813762690 MW P2 = 64,110142224041010 MW P3 = 55,100493405916800 MW Tổn thất công suất tác dụng P_L  = 2,666445472622878 MW Tổng chi phí của toàn hệ thống C = 1599,965715795177 $/h Chi phí trên một đơn vị điện năng λ = 7,767607813261342 $/MWh Đáp án của ví dụ trong tài liệu [1] như sau: Kết luận: So sánh các kết quả sau khi chạy chương trình với đáp án ta có thể đưa ra được kết luận là chương trình chạy đúng. Ứng dụng chương trình Ở các chương trước chúng ta đã nghiên cứu lý thuyết, xây dựng nên chương trình tính toán tối ưu hóa phát điện và áp dụng chương trình để tính toán một lưới điện nhỏ, đơn giản. Trong chương này chúng ta sẽ xét ứng dụng của chương trình để mô phỏng lưới điện phức tạp hơn và so sánh hàm chi phí trước và sau khi đã tối ưu. Ứng dụng chương trình Một lưới điện gồm 26 nút như hình vẽ Hình 0.1: Sơ đồ lưới điện 26 nút Nút 1 là nút cân bằng công suất với điện áp pu. Các bảng dữ liệu như sau: Bảng 0.1: Dữ liệu các nút Tên nút Loại nút Giá trị điện áp P nguồn (MW) Q nguồn (MVAr) P tải (MW) Q tải (MVAr) Qmin, MVAr QMax, MVAr 1 SL 1,025 0 0 51 37 40 250 2 PV 1,02 79 0 22 15 40 150 3 PV 1,025 20 0 64 50 40 80 4 PV 1,05 100 0 25 8 40 160 5 PV 1,045 300 0 50 25 15 50 6 PQ 1 0 0 76 27 40 250 7 PQ 1 0 0 0 0 40 150 8 PQ 1 0 0 0 0 40 80 9 PQ 1 0 0 89 47 40 160 10 PQ 1 0 0 0 0 15 50 11 PQ 1 0 0 25 13,5 40 250 12 PQ 1 0 0 89 46 40 150 13 PQ 1 0 0 31 15 40 80 14 PQ 1 0 0 24 12 40 160 15 PQ 1 0 0 70 30,5 15 50 16 PQ 1 0 0 55 27 40 250 17 PQ 1 0 0 78 38 40 150 18 PQ 1 0 0 153 67 40 80 19 PQ 1 0 0 75 10 40 160 20 PQ 1 0 0 48 27 15 50 21 PQ 1 0 0 46 23 40 250 22 PQ 1 0 0 45 22 40 150 23 PQ 1 0 0 25 12 40 80 24 PQ 1 0 0 54 27 40 160 25 PQ 1 0 0 28 13 15 50 26 PV 1,015 60 0 40 20 15 50 Bảng 0.2: Dữ liệu máy biến áp Dữ liệu MBA Nhánh dây có MBA Tỉ số phân áp 2-3 0,96 2-13 0,96 3-13 1,017 4-8 1,05 4-12 1,05 6-19 0,95 7-9 0,95 Bảng 0.3: Dữ liệu tụ bù Dữ liệu tụ bù Tên nút bù Dung lượng bù (MVAr) 1 4,0 4 2,0 5 5,0 6 2,0 9 3,0 11 1,5 12 2,0 15 0,5 19 5,0 Các hàm phương trình chi phí tại các nút có nguồncủa hệ thống : Bảng 0.4: Giới hạn công suất tác dụng tại các nút Giới hạn công suất tác dụng của máy phát Nút Min, MW Max, MW 1 100 500 2 50 200 3 80 300 4 50 150 5 50 200 6 50 120 Bảng 0.5: Thông số của các đường dây Nhánh Nút đầu Nút cuối R (Ω) X (Ω) B (S) 1 1 2 0,00055 0,0048 0,06 2 1 18 0,0013 0,0115 0,05 3 2 3 0,00146 0,0513 0,04 4 2 7 0,0103 0,0586 0,04 5 2 8 0,0074 0,0321 0,03 6 2 13 0,00357 0,0967 0,02 7 2 26 0,0323 0,1967 0,05 8 3 13 0,0007 0,00548 0,06 9 4 8 0,0008 0,024 0,05 10 4 12 0,0016 0,0207 0,04 11 4 6 0,0069 0,03 0,04 12 5 7 0,00535 0,0306 0,03 13 6 11 0,0097 0,057 0,02 14 6 18 0,00374 0,0222 0,05 15 6 19 0,0035 0,066 0,06 16 6 21 0,005 0,09 0,05 17 7 8 0,0012 0,00693 0,04 18 7 9 0,00095 0,0429 0,04 19 8 12 0,002 0,018 0,03 20 9 10 0,00104 0,0493 0,02 21 10 12 0,00247 0,0132 0,05 22 10 19 0,0547 0,236 0,06 23 10 20 0,0066 0,016 0,05 24 10 22 0,0069 0,0298 0,04 25 11 25 0,096 0,27 0,04 26 11 26 0,0165 0,097 0,03 27 12 14 0,0327 0,0802 0,02 28 12 15 0,0018 0,0598 0,05 29 13 14 0,0046 0,0271 0,06 30 13 15 0,0116 0,061 0,05 31 13 16 0,01793 0,0888 0,04 32 14 15 0,0069 0,0382 0,04 33 15 16 0,0209 0,0512 0,03 34 16 17 0,099 0,06 0,02 35 16 20 0,0239 0,0585 0,05 36 17 18 0,0032 0,06 0,06 37 17 21 0,229 0,445 0,05 38 19 23 0,03 0,131 0,04 39 19 24 0,03 0,125 0,04 40 19 25 0,119 0,2249 0,03 41 20 21 0,0657 0,157 0,02 42 20 22 0,015 0,0366 0,05 43 21 24 0,0476 0,151 0,06 44 22 23 0,029 0,099 0,05 45 22 24 0,031 0,088 0,04 46 23 25 0,0987 0,1168 0,04 Giải bài toán tìm công suất nút và tối ưu hóa điều độ phát điện. Kết quả thu được sau khi giải tích lưới điện bằng phương pháp Newton-Raphson Bảng 0.6: Kết quả mô đun và góc pha điện áp Nút Modun Góc pha(0) 1 1,025 0 2 1,02 -0,71731 3 1,025 -3,88424 4 1,05 -3,41132 5 1,045 2,148035 6 1,006713 -4,22773 7 1,005018 -1,73444 8 1,002967 -2,36884 9 1,016792 -4,41207 10 0,99342 -5,22779 11 1,005716 -4,5216 12 0,995268 -4,23472 13 1,005916 -4,10678 14 0,996307 -4,68162 15 0,990114 -5,20451 16 0,98261 -5,67823 17 0,986468 -5,36076 18 1,009554 -2,39561 19 1,010932 -7,29483 20 0,984245 -5,83542 21 0,981155 -6,62352 22 0,982561 -6,35568 23 0,985197 -7,55865 24 0,973396 -7,80414 25 0,987526 -7,75229 26 1,015 -2,55509 Bảng 0.7: Công suất phát trước khi tối ưu Nút P (MW) Q (Mvar) 1 669,0314 205,4919 2 79 117,8168 3 20 2,796625 4 100 137,4694 5 300 125,5916 26 60 24,1115 Từ kết quả trên ta có chi phí của từng nhà máy là: Vậy tổng chi phí của toàn hệ thống trước tối ưu C = 17916,33054987297 Ta chạy chương trình tối ưu hóa điều độ phát điện Bảng 0.8: Hệ số tổn thất trên lưới B = 0,0014 0,0014 0,0006 -0,0002 -0,0003 -0,0001 0,0014 0,0039 0,0011 -0,0000 -0,0003 -0,0003 0,0006 0,0011 0,0028 -0,0001 -0,0001 -0,0005 -0,0002 -0,0000 0,0001 0,0048 0,0010 -0,0002 -0,0003 -0,0003 -0,0001 0,0010 0,0073 -0,0006 -0,0001 -0,0003 -0,0005 -0,0002 -0,0076 0,0019 B0 = -0,0004 -0,0014 0,0003 0,0021 -0,0008 -0,0004 B00 = 0,0142 Kết quả thu được sau khi chạy chương trình chính tính toán tối ưu hóa phát điện Công suất phát tối ưu của các nhà máy P1 = 447.6008 MW P2 = 164.7873 MW P3 = 261.6360 MW P4 = 127.7905 MW P5 = 170.5165 MW P26 = 104.0022 MW Tổn thất công suất tác dụng P_L  = 13.3333 MW Chi phí của từng nhà máy là: Tổng chi phí của toàn hệ thống C = 17456,31034124316 $/h Chi phí trên một đơn vị điện năng λ = 13.5140 $/MWh Nhận xét: Ta có bảng so sánh chi phí nhiên liệu trong các trường hợp Bảng 0.9: Bảng tổng kết kết quả so sánh chi phí nhiên liệu trong các trường hợp C ($/h) Trước tối ưu 17916,33054987297 Sau tối ưu 17456,31034124316 Như vậy so với trường hợp chạy ở chế độ xác lập thì chi phí sau khi tính toán tối ưu hóa đã giảm là: 17916,33054987297 - 17456,31034124316 = 460,020208629798 ($/h) Từ các bảng kết quả ở trên ta thấy sau khi sử dụng chương trình tối ưu hóa điều độ phát điện thì có sự giảm về chi phí nhiên liệu từ 17916,33054987297 $/h xuống 17456,31034124316 $/h. Như vậy, chi phí nhiên liệu tiết kiệm được có thể chiếm một vài phần trăm của chi phí cho việc phát điện. Như vậy, tối ưu hóa điều độ phát điện đã mang lại hiệu quả đáng kể cho ngành HTĐ. Kết luận chung Như vậy, qua các bài toán ví dụ ở trên đã mô tả rõ nét về bài toán tối ưu hóa điều độ phát điện trong HTĐ. Cách giải quyết bài toán bằng phương pháp lập trình sử dụng ngôn ngữ Matlab, trong đó có sử dụng các thuật toán Newton-Raphson giải tích lưới điện ở chế độ xác lập và phương pháp hệ số Lagrange để giải bài toán tối ưu hóa cho ta kết quả tối ưu nhất đáp ứng các yêu cầu đặt ra. Chương trình có nhiều ưu điểm: Tính ứng dụng cao trong thực tế cho lưới điện từ đơn giản đến phức tạp mang lại nhiều lợi ích về kinh tế - kỹ thuật cho ngành điện cũng như các ngành kỹ thuật khác. Mang lại hiệu quả cao với độ chính xác tin cậy Đảm bảo điện năng cung cấp cho các phụ tải theo yêu cầu đặt ra. Luôn cho kết quả tối ưu nhất. Dễ sử dụng Hiện tại chương trình còn có nhiều giới hạn như hàm chi phí đưa ra mới chỉ xét đến chi phí nhiên liệu mà chưa tính đến các chi phí khác. Mới chỉ giải quyết được bài toán tối ưu theo chi phí nhiên liệu mà chưa xét đến các chỉ tiêu khác như: tối ưu về tổn thất công suất, tối ưu điện áp…Việc kết hợp các bài toán tối ưu đó đòi hỏi chương trình phức tạp hơn. TÀI LIỆU THAM KHẢO Power System Analysis – Hadi Saadat Phân tích chế độ xác lập hệ thống điện – PGS.TS Phạm Văn Hòa – Nhà xuất bản Bách Khoa-Hà Nội Lưới điện & Hệ thống điện (tập 1) – Trần Bách – Nhà xuất bản khoa học và kỹ thuật Matlab ứng dụng – Trần Quang Khánh – Nhà xuất bản khoa học và kỹ thuật PHỤ LỤC : Chương trình tính toán các hệ số tổn thất trên lưới bloss function [B B0 B00 nbrGen] = Bloss(Snut, Stai, Ybus, code, Unut); %format long Snut = Snut/100; %matra (1xn) Stai = Stai/100; %matra (1xn) Sng = Snut + Stai;%matra (1xn) Itai=-conj(Stai./Unut); %tinh Itai = - *Stai/*U; %matra (1xn) chu y Unut chinh la U Id = sum(Itai); %Tinh Id = Tong Itai l_k = Itai/Id; %l_k =Itai/Id %matra (1xn) Zbus = inv(Ybus); %Zbus = -1/Ybus Rbus = real(Zbus); %Rbus = phan thuc Zbus T=0; for i=1:length(l_k) j=i; lkz = l_k(j)*Zbus(1,j); T = T + lkz; end rho = l_k./(-1*T); nbrGen = 0; nbrLoad = 0; for i = 1:length(code) if code(i)==2 nbrLoad = nbrLoad +1; %neu code =2 nut nbrLoad else nbrGen = nbrGen +1 end end %xu ly Sng Sngmoi = zeros(0,nbrGen); tt=0; for i=1:nbrGen+nbrLoad if code(i)==0 tt=tt+1; Sngmoi(tt) = Sng(i); end if code(i)==1 tt=tt+1; Sngmoi(tt) = Sng(i); end end C = zeros(length(Snut), nbrGen +1); Cio = rho*Zbus(1,1); %tinh Ci0 = rho(i)*Zbus(1,1) for i=1:length(Snut) C(i,(nbrGen+1))= Cio(i); end %------------------------------------------- for i=1:length(Snut) for k=1:nbrGen if i==k C(i,k) = 1 + rho(i)*Zbus(1,i) else C(i,k) = rho(i)*Zbus(1,k) end end end %------------- for i=1:nbrGen C(i,i) = 1 + rho(i)*Zbus(1,i); end %================ %--------------------------------------------- Psi = zeros((nbrGen+1),(nbrGen+1)); Psi((nbrGen+1),(nbrGen+1)) = -1*Unut(1,1)/Zbus(1,1); for m =1:nbrGen tgt = 1 -1i*imag(Sngmoi(m))/real(Sngmoi(m)); Psi(m,m) = tgt/conj(Unut(m)); end %-------------------------------------------- H = transpose(Psi)*transpose(C)*Rbus*conj(C)*conj(Psi); ReH = (H + conj(H))/2 %------------------------------------------- B = zeros(nbrGen,nbrGen); for i=1:nbrGen for j=1:nbrGen B(i,j) = ReH(i,j) end end %-------------------------------------------- B00 = ReH(nbrGen+1, nbrGen+1) %------------------------------------------- B0 = zeros(1,nbrGen); for i=1:nbrGen B0(i) = 2*ReH(nbrGen+1,i) end %=====================The end =========== : Chương trình chính tính toán tối ưu hóa sau khi đã giải tích lưới điện ở chế độ xác lập function [P_matrix, c_matrix,lamda,P_D,P_L,PL_matrix,Pk,C] = mainfunc(epsilon,tongSoNut, nbrGen, Snut, Stai, code, alpha, beta, gamma, ... Pmin, Pmax, B, B0, B00) P_D = sum(real(Stai)); P=[]; for i = 1:tongSoNut if(code(i)~=2) P=[P, real(Snut(i))]; end end %======================================= % Tinh toan lambda for i = 1:nbrGen tmp1 = beta(i)/(2*gamma(i)); tmp2 = 1/(2*gamma(i)); end lamda = (P_D+tmp1)/tmp2; %================================================= % Can bang gia tri B Sb = 100; B = B/Sb; B00 = B00*Sb; %============================================================= % Cac gia tri Pk: gia tri P tinh cho tung vong lap % P_matrix: matran luu gia tri P cho tat ca cac lan % P_final: ma tran kiem tra ket qua, khoi tao = [0 0...0] % neu gia tri Pi nao thoa man, vi tri tuong ung % chuyen thanh 1. khi = [1 1 ...] ket thuc main Pk = []; % matrix for Pk_i, ma tran tinh toan cho tung vong lap Pk1= P; % gia tri P tai buoc lap dau tien P_matrix = []; PL_matrix = []; P_final = zeros(1, nbrGen); % final P % So vong lap toi da maxLoop = 100; nbrLoop = 0; % bien dem so lan lap flag =0; % co xac dinh trang thai lap % while nbrLoop < maxLoop while(flag == 0 && nbrLoop < maxLoop) % while(flag == 0) tmp = 0; nbrLoop = nbrLoop +1 for i = 1:nbrGen tmp = 0; for j = 1:nbrGen if j ~= i tmp = tmp + B(i,j)*Pk1(j); end end tuso = lamda*(1- B0(i)) - beta(i)-2*lamda*tmp; mauso = 2*(gamma(i) + lamda*B(i,i)); var = tuso/mauso; if(P_final(i) == 0) Pk(i) = var; end; end; % Tinh toan cac gia tri P_L loopStatus = 0; while (loopStatus == 0) %Cal culate P_L P_L = 0; for i = 1:nbrGen P_L = P_L + B0(i)*Pk(i); end for i = 1:nbrGen for j = 1:nbrGen P_L = P_L + Pk(i)*B(i,j)*Pk(j); end end P_L = P_L + B00; delta_P = P_D + P_L - sumCal(Pk) % So sanh delta_P and epsilon if (abs(delta_P) <=epsilon) for i = 1:nbrGen % Vong lap tinh toan cac phan tu Pk(i) if(P_final(i)==0) % neu gia tri Pk(i) chua chinh xac, tinh lai if (Pmin(i) thoa man P_final(i) = 1; % gan P_final(i) = 1 loopStatus = 1; % ngung vong lap nay else if (Pk(i)<Pmin(i)) % neu ko, gan gia tri = min hoac max Pk(i)=Pmin(i); P_final(i) = 1; end; if (Pk(i)>Pmax(i)) Pk(i)=Pmax(i); P_final(i) = 1; break; end; end; end; end; else loopStatus = 1; % tinh toan gia tri delta_lamda %Tinh dao ham rieng Pk theo lamda Atmp = 0; if (P_final(i)== 0) for i = 1:nbrGen tmp = 0; for j = 1:nbrGen if j ~= i tmp = tmp + B(i,j)*Pk(j); end end Atmp = Atmp + (gamma(i)*(1-B0(i)) + B(i, i)*beta(i) - 2*gamma(i)*tmp)/(2*(gamma(i) + lamda*B(i,i))^2); end end % Tinh Delta lamda delta_lamda = delta_P/Atmp lamda = lamda + delta_lamda; end end; % Luu ket qua vao mang P_matrix = [P_matrix;Pk]; Pk1= Pk; % kiem tra xem tinh duoc chua if(sum(P_final)== nbrGen) flag=1; break; end; end %Cal Cost value c_matrix =[]; for i = 1:nbrGen tmp = alpha(i)+beta(i)*Pk(i) + gamma(i)*(Pk(i)^2); c_matrix = [c_matrix, tmp]; end C = sum(c_matrix); : Chương trình tổng quát tính toán tối ưu hóa điều độ phát điện clc; clear all; % Vi du 7.9 --- % Doc du lieu nut [tongSoNut, nut, Udm, Ucb, Scb, code, modunU, phaU, Sng, Stai, Qmin, Qmax, Ys] ... = readBusData('filetest_BusData.xls') maxNut = max(nut); % Doc du lieu duong day [tongSoDuongDay, soDuongDay, id, nutDau, nutCuoi, UcbDZ, ScbDZ, chieuDai, ... Z0, Y0ngang, heSoPhanAp]= readLineData(maxNut, 'filetest_LineData.xls'); % doc du lieu ham chi phi va cong suat gioi han [alpha,beta,gamma,Pmax,Pmin] ... = readData('filetest_Data.xls'); % Tinh ma tran tong dan Ybus [Ybus] = ybus(tongSoDuongDay, soDuongDay, id, nutDau, nutCuoi, UcbDZ, ... ScbDZ, chieuDai, Z0, Y0ngang, heSoPhanAp); %================================================= % Tinh Toan GS % Giai tich luoi dien bang NR epsilon = 1e-10; MaxNbrIte = 1000; accel = 1.8; [magU, theta, nbrIte, doLechCongSuat] = plfNR(tongSoNut, nut, Udm, Ucb, Scb, code, modunU, phaU, Sng, Stai, Qmin, Qmax, Ybus, ... epsilon, MaxNbrIte); U = magU(nbrIte,:).*exp(1i*theta(nbrIte,:)); % Tinh cong suat nut [Snut] = calculateSbus(tongSoNut, code, Sng, Stai, Scb, U... , Ybus); %tinh ma tran ton that bloss [B, B0, B00, nbrGen ] = Bloss(Snut, Stai, Ybus, code, U); % chuong trinh chinh epsilon1 = 0.5e-5; [P_matrix, c_matrix,lamda,P_D,P_L,PL_matrix,Pk,C] = mainfunc(epsilon1, tongSoNut, nbrGen, Snut, Stai, code, alpha, beta, gamma, ... Pmin, Pmax, B, B0, B00);