BỘ GIÁO DỤC VÀ ĐÀO TẠO
TRƯỜNG ĐẠI HỌC CÔNG NGHỆ TP. HCM
---------------
VÕ KHÁNH DƯƠNG
ÁP DỤNG THUẬT TOÁN GRAVITATIONAL
SEARCH ALGORITHM TÍNH TOÁN PHÂN BỐ
CÔNG SUẤT TỐI ƯU TRONG HỆ THỐNG ĐIỆN
LUẬN VĂN THẠC SỸ
Chuyên ngành: KỸ THUẬT ĐIỆN
Mã số ngành: 60520202
TP. HỒ CHÍ MINH, Tháng 7 năm 2016
BỘ GIÁO DỤC VÀ ĐÀO TẠO
TRƯỜNG ĐẠI HỌC CÔNG NGHỆ TP. HCM
---------------
VÕ KHÁNH DƯƠNG
ÁP DỤNG THUẬT TOÁN GRAVITATIONAL
SEARCH ALGORITHM TÍNH TOÁN PHÂN BỐ
CÔNG SUẤT TỐI ƯU TRONG HỆ THỐNG ĐIỆN
LUẬN VĂN THẠC SỸ
Chuyên ngành: KỸ THUẬT ĐIỆN
Mã số ngành: 60520202
CÁN BỘ HƯỚNG DẪN KHOA HỌC: PGS.TS.NGÔ CAO CƯỜNG
THS.LÊ ĐÌNH LƯƠNG
TP. HỒ CHÍ MINH, Tháng 7 năm 2016
CÔNG TRÌNH ĐƯỢC HOÀN THÀNH TẠI
TRƯỜNG ĐẠI HỌC CÔNG NGHỆ TP. HCM
Cán bộ hướng dẫn khoa học : PGS.TS.NGÔ CAO CƯỜNG
Luận văn Thạc sĩ được bảo vệ tại Trường Đại học Công nghệ TP. HCM
ngày 25 tháng 09 năm 2016.
Thành phần Hội đồng đánh giá Luận văn Thạc sĩ gồm:
TT Họ và tên Chức danh Hội đồng
1 PGS.TS.Quyền Huy Ánh Chủ tịch
2 TS.Nguyễn Hùng Phản biện 1
3 TS.Võ Công Phương Phản biện 2
4 PGS.TS.Đồng Văn Hướng Ủy viên
5 TS.Võ Viết Cường Ủy viên, Thư ký
Xác nhận của Chủ tịch Hội đồng đánh giá Luận sau khi Luận văn đã được
sửa chữa (nếu có).
Chủ tịch Hội đồng đánh giá Luận văn
TRƯỜNG ĐH CÔNG NGHỆ TP. HCM CỘNG HÒA XÃ HỘI CHỦ NGHĨA VIỆT NAM
PHÒNG QLKH – ĐTSĐH Độc lập – Tự do – Hạnh phúc
TP. HCM, ngày 31 tháng 07 năm 2016
NHIỆM VỤ LUẬN VĂN THẠC SĨ
Họ tên học viên: Võ Khánh Dương. Giới tính: Nam.
Ngày, tháng, năm sinh: 16/10/1991. Nơi sinh: Quảng Ngãi.
Chuyên ngành: Kỹ Thuật Điện. MSHV:1441830033.
I- Tên đề tài:
Áp dụng thuật toán Gravitational Search Algorithm tính toán phân bố công suất tối
ưu trong hệ thống điện.
II- Nhiệm vụ và nội dung:
Nghiên cứu lý thuyết về thuật toán GSA từ những công trình đã công bố
trước đây trên tạp chí khoa học thế giới.
Nghiên cứu cách áp dụng thuật toán GSA vào tính toán trong hệ thống
điện.
Xây dựng giải thuật GSA giải bài toán điều phối kinh tế công suất ED.
Lập trình tính toán điều phối tối ưu trong mạng điện.
Nhận xét, đánh giá kết quả thu được, so sánh với kết quả dùng các giải thuật
khác đã công bố trên tạp chí khoa học trên thế giới.
III- Ngày giao nhiệm vụ : 23/01/2016.
IV- Ngày hoàn thành nhiệm vụ : 31/07/2016.
V- Cán bộ hướng dẫn : PGS.TS. Ngô Cao Cường.
Thạc sĩ Lê Đình Lương
CÁN BỘ HƯỚNG DẪN KHOA QUẢN LÝ CHUYÊN NGÀNH
BỘ GIÁO DỤC VÀ ĐÀO TẠO TRƯỜNG ĐH CÔNG NGHỆ TP. HCM CỘNG HÒA XÃ HỘI CHỦ NGHĨA VIỆT NAM Độc lập – Tự do – Hạnh phúc
TP. Hồ Chí Minh, ngày 28 tháng 10 năm 2016.
BẢN CAM ĐOAN
Họ và tên học viên: VÕ KHÁNH DƯƠNG.
Ngày sinh: 16/10/1991. Nơi sinh: Quảng Ngãi.
Trúng tuyển đầu vào năm: 2014
Là tác giả luận văn: Áp dụng thuật toán Gravitational Search Algorithm tính
toán phân bố công suất tối ưu trong hệ thống điện.
Chuyên ngành: Kỹ Thuật Điện. Mã ngành: 60520202.
Bảo vệ ngày: 25 tháng 09 năm 2016.
Điểm bảo vệ luận văn: 8,7.
Tôi cam đoan chỉnh sửa nội dung luận văn thạc sĩ với đề tài trên theo góp ý của Hội
đồng đánh giá luận văn Thạc sĩ. Các nội dung đã chỉnh sửa:
- Chỉnh sửa bố cục lời mở đầu và chương 1 của đề tài không còn bị trùng lắp về mặt nội dung.
- Bổ sung sơ đồ nhất thứ tổng quan của các mạng điện 3 nút, 13 nút và 40 nút.
Người cam đoan Cán bộ Hướng dẫn (Ký, ghi rõ họ tên) (Ký, ghi rõ họ tên)
i
LỜI CAM ĐOAN
Tôi xin cam đoan Luận văn tốt nghiệp: “Áp dụng thuật toán Gravitational
Search Algorithm tính toán phân bố công suất tối ưu trong hệ thống điện” là đề
tài nghiên cứu do bản thân tôi tự thực hiện, không sao chép dưới bất kỳ hình thức
nào.
Những kết quả và các số liệu trong Luận văn đều được lấy từ những nguồn
tài liệu chính thống và có uy tín đã được tôi trích dẫn đầy đủ và ghi chép rõ ràng
trong phần tài liệu tham khảo. Tôi hoàn toàn chịu trách nhiệm trước nhà trường về
sự cam đoan này.
TP.HCM, Ngày 31 tháng 07 năm 2016
Học viên thực hiện
VÕ KHÁNH DƯƠNG
ii
LỜI CẢM ƠN
Lời đầu tiên, Em xin được cảm ơn Quý Thầy Cô trường Đại học Công Nghệ
Thành phố Hồ Chí Minh, những người đã tận tình truyền đạt kiến thức chuyên
ngành cho em, cũng như những kiến thức trong xã hội trong suốt thời gian em học
tập tại trường.
Em cũng xin chân thành cảm ơn Thầy Ngô Cao Cường và Thầy Lê Đình
Lương, những người đã tận tình dìu dắt, bổ sung những kiến thức thực tế, nhiệt tình
giúp đỡ và hướng dẫn em trong suốt thời gian làm Luận văn tốt nghiệp.
Cuối cùng, Em xin chúc các Thầy Cô giáo ở trường lời chúc sức khỏe và công
tác tốt. Chúc các thầy, cô trong Phòng Quản lý khoa học và Đào tạo sau đại học lời
chúc tốt đẹp nhất.
TP.HCM, Ngày 31 tháng 07 năm 2016
Học viên thực hiện
VÕ KHÁNH DƯƠNG
iii
TÓM TẮT
Trong những năm gần đây, sự phát triển của khoa học công nghệ khiến cho
nhu cầu năng lượng trong các hệ thống điện tăng cao và làm cho các hệ thống này
ngày càng phức tạp. Điều phối kinh tế (ED), một trong những vấn đề tối ưu hóa phi
tuyến trong các hệ thống năng lượng, có một vị trí quan trọng trong hoạt động kinh
tế của hệ thống điện. Trong việc giải quyết các vấn đề của ED, mục tiêu là để giảm
thiểu tổng chi phí nhiên liệu, trong khi vẫn đáp ứng các yêu cầu về vật lý và các
ràng buộc khác nhau.
Luận văn này trình bày lý thuyết về thuật toán Gravitational Search
Algorrithm (GSA) và nêu tổng quan các thuật toán tối ưu khác đã được áp dụng để
giải bài toán ED. Trên cơ sở đó, áp dụng thuật toán Gravitational Search Algorithm
(GSA) giải các bài toán điều phối kinh tế trong hệ thống điện như sau:
Bài toán ED hệ thống 3 nút xét đến điểm van công suất có tổng công suất
nhu cầu của phụ tải PD = 850 (MW), tổn thất PL = 0 (MW).
Bài toán ED hệ thống 13 nút xét đến điểm van công suất có tổng công suất
nhu cầu của phụ tải PD = 1800 (MW), tổn thất PL = 0 (MW).
Bài toán ED hệ thống 40 nút xét đến điểm van công suất có tổng công suất
nhu cầu của phụ tải PD = 10500 (MW), tổn thất PL = 0 (MW).
Tiến hành so sánh kết quả tính toán của thuật toán GSA với kết quả của các
thuật toán khác rút ra nhận xét và kinh nghiệm lập trình để áp dụng cho phù hợp
với các bài toán khác.
Vạch ra những hướng phát triển nghiên cứu, hướng tiếp cận mới để tiếp tục
cải thiện thuật toán GSA đưa đến kết quả tốt hơn, cũng như ứng dụng vào các bài
toán có quy mô lớn và phức tạp hơn trong hệ thống điện.
iv
ABSTRACT
In recent years, the developments in technology cause more energy demand
in power systems and make these systems more complicated. The Economic
Dispatch (ED) problem, one of the nonlinear optimization problems in electrical
power systems, has an important role in the economical operation of the power
system. For solving the ED problem, the objective is to minimize the total fuel cost,
while satisfying the various physical and operational constraints.
This thesis presents the theory of algorithms Gravitational Search Algorrithm
(GSA) and provides an overview of different optimization algorithms which have
been applied to solve the ED problem. Based on which, apply the Gravitational
Search Algorithm algorithm (GSA) to solve the problems of economic coordination
(ED) in the electrical system as follows:
The problem of ED 3-units system considering the capacity valve points with
the total capacity of the load demand PD = 850 (MW), losses PL = 0 (MW).
The problem of ED 13-units system considering the capacity valve points
with the total capacity of the load demand PD = 1800 (MW), losses PL = 0
(MW).
The problem of ED 40-units system considering the capacity valve points
with the total capacity of the load demand PD = 10500 (MW), losses PL = 0
(MW).
Conducted comparing the results of the algorithm calculates the GSA with the results of other algorithms draws comment and programming experience to apply for matching other problems.
Outlines the development of research, new approaches to further improve
GSA algorithm leads to better results, as well as apply to the more complex problem
in the power system.
v
MỤC LỤC
LỜI CAM ĐOAN .......................................................................................................i
LỜI CẢM ƠN ........................................................................................................... ii
TÓM TẮT ................................................................................................................ iii
ABSTRACT...............................................................................................................iv
MỤC LỤC .................................................................................................................. v
CHỮ VIẾT TẮT TRONG LUẬN VĂN .............................................................. viii
DANH MỤC CÁC BẢNG BIỂU ............................................................................ix
DANH MỤC CÁC HÌNH ......................................................................................... x
CHƯƠNG 1 TỔNG QUAN ĐỀ TÀI ....................................................................... 1
1.1 ĐẶT VẤN ĐỀ ............................................................................................... 1
1.2 TÓM TẮT MỘT SỐ TÀI LIỆU VÀ BÀI BÁO LIÊN QUAN ................ 3
1.3 HƯỚNG TIẾP CẬN ĐỀ TÀI ..................................................................... 4
1.4 MỤC TIÊU CỦA ĐỀ TÀI .......................................................................... 5
1.5 ĐIỂM MỚI CỦA ĐỀ TÀI .......................................................................... 5
CHƯƠNG 2 GIỚI THIỆU BÀI TOÁN ĐIỀU PHỐI KINH TẾ TRONG HỆ
THỐNG ĐIỆN ........................................................................................................... 6
2.1 THUẬT TOÁN TỐI ƯU: .............................................................................. 6
Hình 2.1: Minh hoạ tối ưu toàn cục hàm Peaks. .................................................... 7
2.2 BÀI TOÁN ĐIỀU PHỐI KINH TẾ CỔ ĐIỂN: ........................................ 7
2.2.1 Ràng buộc đẳng thức: .............................................................................. 9
2.2.2 Ràng buộc bất đẳng thức: ..................................................................... 10
2.2.2.1 Giới hạn công suất thực phát ra: ....................................................... 10
2.2.2.2 Giới hạn tốc độ: ................................................................................. 10
2.2.2.3 Giới hạn về vùng cấm vận hành: ....................................................... 10
2.2.3 Ràng buộc về công suất truyền tải: ......................................................... 11
2.2.4 Bài toán điều phối kinh tế với hàm chi phí nhiên liệu không trơn: .. 11
2.2.4.1 Bài toán điều phối kinh tế có điểm van công suất: ....................... 12
2.2.4.2 Biểu thức điều phối kinh tế với điểm van công suất: ................... 12
vi
2.3 TỔNG QUAN CÁC PHƯƠNG PHÁP ĐÃ ĐƯỢC ÁP DỤNG ĐỂ GIẢI BÀI TOÁN ĐIỀU PHỐI KINH TẾ (ED): ........................................................ 13
2.3.1 DE ( Differential Evolution ) [22][23] ............................................... 13
2.3.2 PSO ( Particle Swarm Optimization ) [24][25][26] .......................... 13
2.3.3 ABC (Thuật toán Artificial Bee Colony) [27][28][29] ..................... 13
2.3.4 HNN ( Hopfield Neuron Network ) [30][31] ..................................... 14
2.3.5 ELANN ( Enhanced Lagrangian Artificial Neural Network) [32]. 14
2.3.6 HS (Harmony Search) [33] ................................................................ 14
2.3.7 CS (Cuckoo Search) [34][35][36][37] ................................................ 15
CHƯƠNG 3 THUẬT TOÁN GRAVITAIONAL SEARCH ALGORITHM ... 16
3.1 THUẬT TOÁN GSA CỔ ĐIỂN: ................................................................ 16
3.2 CÁC BƯỚC TRONG THUẬT TOÁN GSA: .......................................... 21
3.3 CÁC ĐẶC TÍNH CỦA GIẢI THUẬT GSA ........................................... 23
3.4 ẢNH HƯỞNG CỦA CÁC THAM SỐ TRONG GSA ............................... 23
3.4.1 Hệ số gia tốc trọng trường G0: .......................................................... 24
3.4.2 Hệ số suy giảm α : ............................................................................... 24
3.4.3 Số cá thể trong tập hợp N: ................................................................. 24
3.5 KẾT LUẬN ................................................................................................. 24
CHƯƠNG 4 ỨNG DỤNG THUẬT TOÁN GSA GIẢI BÀI TOÁN ĐIỀU PHỐI
KINH TẾ TRONG HỆ THỐNG ĐIỆN ................................................................ 26
4.1 GIẢI BÀI TOÁN ĐIỀU PHỐI KINH TẾ (ED) VỚI HÀM CHI PHÍ CÓ XÉT ĐIỂM VAN CÔNG SUẤT: ................................................................ 26
4.1.1 Các bước áp dụng thuật toán GSA giải bài toán ED: ..................... 26
4.1.2 Lưu đồ giải thuật của thuật toán GSA: ............................................ 28
4.1.3 Hàm mục tiêu của bài toán ED: ........................................................ 30
4.2 GIẢI BÀI TOÁN ED MẠNG 3 NÚT: ........................................................ 31
4.3 GIẢI BÀI TOÁN ED MẠNG 13 NÚT: ...................................................... 34
4.4 GIẢI BÀI TOÁN ED 40 MẠNG NÚT: ..................................................... 38
4.5 KẾT LUẬN: ............................................................................................... 43
CHƯƠNG 5 KẾT LUẬN VÀ HƯỚNG PHÁT TRIỂN ....................................... 45
5.1 TỔNG KẾT .................................................................................................... 45
vii
5.2 HƯỚNG PHÁT TRIỂN ĐỀ TÀI ............................................................. 46
5.3 LỜI KẾT .................................................................................................... 46
TÀI LIỆU THAM KHẢO ........................................................................................ 47
viii
CHỮ VIẾT TẮT TRONG LUẬN VĂN
Harmony Search HS
Economic Load Dispatch ELD
Nonlinear programming NLP
Linear programming LP
Evolutionary programming EP
Genetic algorithm GA
Improved Particle Swarm Optimization IPSO
Quadratic programming QP
Simulated Annealing- Particle Swarm Optimization SA-PSO
Interior Point Methods IP
Different Evolution DE
Dynamic Programming DP
Artificial Bee Colony ABC
Ant Colony Optimization ACO
Harmony Memory HM
Genetic Algorithm GA
Cuckoo Search CS
Gravitational Search Algorithm GSA
Hopfield Neuron Network HNN
Enhanced Lagrangian Artificial Neural Network ELNN
ix
DANH MỤC CÁC BẢNG BIỂU
Bảng 4.1: Các thông số của bài toán ED mạng 3 nút................................................ 31
Bảng 4.2: Thông số của thuật toán GSA áp dụng giải bài toán ED mạng 3 nút. ...... 31
Bảng 4.3: Kết quả tính toán bài toán ED mạng 3 nút. .............................................. 31
Bảng 4.4: So sánh kết quả tính toán bài toán ED mạng 3 nút. .................................. 32
Bảng 4.5: Các thông số của bài toán ED mạng 13 nút.............................................. 34
Bảng 4.6: Thông số của thuật toán GSA áp dụng giải bài toán ED mạng 13 nút. .... 35
Bảng 4.7: Kết quả tính toán bài toán ED mạng 13 nút. ............................................ 35
Bảng 4.8: So sánh kết quả tính toán bài toán ED mạng 13 nút. ................................ 35
Bảng 4.9: Các thông số của bài toán ED mạng 40 nút.............................................. 39
Bảng 4.10: Thông số của thuật toán GSA áp dụng giải bài toán ED mạng 40 nút. .. 39
Bảng 4.11: Kết quả tính toán bài toán ED mạng 40 nút. .......................................... 39
Bảng 4.12: So sánh kết quả tính toán bài toán ED mạng 40 nút. .............................. 41
x
DANH MỤC CÁC HÌNH
Hình 2.1: Minh hoạ tối ưu toàn cục hàm Peaks. ......................................................... 7
Hình 2.2: Đường cong chi phí phổ biến của nhà máy nhiệt điện................................ 9
Hình 2.3: Đồ thị biểu diễn vùng cấm của nút nhiệt cơ bản. ...................................... 11
Hình 2.4 Hàm chi phí nhiên liệu của nhà máy nhiệt điện với điểm van công suất. . 12
Hình 3.1 Lực và gia tốc tương tác lên vật thể 1 do các vật thể khác sinh ra [10]. .... 19
Hình 3.2 Nguyên lý của giải thuật GSA [10]. ........................................................... 21
Hình 4.1: Lưu đồ giải thuật GSA cho bài toán ED. .................................................. 29
Hình 4.2: Sơ đồ mạng nhất thứ các nhà máy điện phân phối công suất đến phụ tải
thông qua hệ thống truyền tải điện............................................................................30
Hình 4.3: Đồ thị giá trị hàm chi phí của bài toán ED mạng 3 nút. ........................... 32
Hình 4.4: Đồ thị phân bố công suất tại các nút của bài toán ED mạng 3 nút. .......... 33
Hình 4.5: Giao diện tính toán của bài toán ED mạng 3 nút. ..................................... 33
Hình 4.6: Đồ thị giá trị hàm chi phí của bài toán ED mạng 13 nút. ......................... 36
Hình 4.7: Đồ thị phân bố công suất tại các nút của bài toán ED mạng 13 nút. ........ 37
Hình 4.8: Giao diện tính toán của bài toán ED mạng 13 nút. ................................... 37
Hình 4.9: Đồ thị giá trị hàm chi phí của bài toán ED mạng 40 nút. ......................... 42
Hình 4.10: Đồ thị phân bố công suất tại các nút của bài toán ED mạng 40 nút. ...... 42
Hình 4.11: Giao diện tính toán của bài toán ED mạng 40 nút. ................................. 43
1
CHƯƠNG 1 TỔNG QUAN ĐỀ TÀI
1.1 ĐẶT VẤN ĐỀ
Trong những thập kỷ qua, nhu cầu năng lượng điện trên toàn thế giới đã đột
ngột tăng do tăng trưởng kinh tế. Bên cạnh đó các nguồn năng lượng hóa thạch vốn
là nguyên liệu chính để sản xuất diện đang đứng trên nguy cơ cạn kiệt, nguồn cung
không ổn định giá cả biến động. Mục đích của hệ thống điện là nâng cao hiệu quả
sử dụng, chất lượng điện năng, độ tin cậy cấp điện đồng thời giảm chi phí đầu tư,
vận hành và bảo trì. Bài toán điều độ kinh tế được đặt ra bức thiết hơn bao giờ
hết[1].
Phân bố công suất tối ưu trong hệ thống điện là một trong những vấn đề quan
trọng nhất được những người vận hành hệ thống công suất sử dụng phổ biến mỗi
ngày. Công việc này nhằm tìm kiếm, phân phối thông số trạng thái vận hành tối ưu
công suất thực và công suất phản kháng làm giảm chi phí và cải thiện hiệu suất cho
toàn hệ thống. Vấn đề này có thể được trình bày rõ ràng và giải quyết thành hai vấn
đề riêng lẻ. Thứ nhất, điều phối kinh tế làm giảm chi phí hệ thống bằng cách phân
phối hợp lý công suất thực cho các máy phát hiện hành. Thứ hai, điều phối công
suất phản kháng làm giảm tổn thất và cải thiện điện áp hệ thống bằng cách phân
phối hợp lý công suất phản kháng.
Lưới điện bao gồm máy phát điện, máy biến áp, đường dây truyền tải, thiết
bị đóng cắt, rơ le bảo vệ và thiết bị bù công suất phản kháng…. Các hệ thống truyền
tải được sử dụng để truyền công suất đi xa. Việc kiểm soát các mục tiêu khác nhau
như vận hành và thiết kế trong những hệ thống như vậy đòi hỏi về tối ưu hóa. Đối
với hệ thống phi tuyến, giải pháp cho vấn đề tối ưu hóa lại càng cần thiết. Hơn thế
nữa, cần phải chú ý các vấn đề sau:
- Các kỹ thuật tối ưu hóa được lựa chọn nên phù hợp với vấn đề nghiên
cứu.
- Tất cả các khía cạnh khác nhau của vấn đề phải được kể đến.
- Tất cả các ràng buộc của hệ thống phải được trình bày chính xác.
- Phải xác định được hàm mục tiêu.
Tối ưu hóa chiếm một vị trí quan trọng trong hệ thống năng lượng và là một
kỹ thuật thường được sử dụng trong vận hành hệ thống điện. Tối ưu hóa tìm cách
2
phân bố lại công suất thực và công suất phản kháng nhằm làm giảm chi phí nhiên
liệu, giảm lượng khí thải gây ảnh hưởng trực tiếp đến môi trường xung quanh và cải
thiện hiệu quả toàn bộ hệ thống. Việc mô hình hóa bài toán điều phối công suất phát
là chỉ tiêu để đạt tới được kết quả tối ưu. Trong bài điều phối kinh tế, công thức cổ
điển thể hiện các khiếm khuyết do mô hình quá đơn giản. Trong bài công suất phản
kháng, phương pháp thông thường là mô hình hóa các máy biến thế và dàn tụ điện
thành các biến liên tục thay vì các biến rời rạc.Để cải thiện vấn đề này, các mô hình
mới vẫn liên tục được phát triển cho việc vận hành hệ thống thêm hiệu quả. Mức độ
phức tạp của bài tối ưu cũng được nâng lên do liên kết các ràng buộc trở nên phi
tuyến tính.
Các biện pháp để giảm chi phí nhiên liệu trong vận hành là:
- Tăng lượng công suất phát ra của các nhà máy nhiệt điện gần phụ tải
nhằm giảm tổn hao truyền tải, do đó làm giảm chi phí tiêu hao nhiên liệu trên
toàn bộ hệ thống.
- Tăng lượng công suất phát tại các nhà máy nhiệt điện có đặc tính tiêu hao
nhiên liệu thấp.
- Phối hợp giữa các nhà máy nhiệt điện với nhau sao cho chi phí sản xuất
điện năng là nhỏ nhất.
Vì vậy người ta đặt ra bài toán điều phối tối ưu để nâng cao khả năng tận
dụng hệ thống điện hiện có. Đây là bài toán mà ngành điện lực phải tìm cách giải
quyết từ rất lâu, đã dùng nhiều thuật toán cổ điển như: Linear Programming (Lập
trình tuyến tính) [2, 3], Nonlinear Programming (Lập trình phi tuyến), Newton-
Raphson [4]. Những cải tiến gần đây trong việc giải quyết các bài toán tối ưu phức
tạp với kết quả chính xác hơn đã làm phát triển các kỹ thuật mới mang tên Thuật
toán tiến hóa. Thuật toán Tiến hóa là kỹ thuật tối ưu dựa trên nền tảng tìm kiếm đáp
án ngẫu nhiên bằng việc sử dụng mô hình được đơn giản hóa trong tiến trình tiến
hóa, cho ra được kết quả tối ưu toàn bộ, đặc biệt trong các khoảng không đáp án
không liên tục, không lồi và phi tuyến tính cao. Nó dựa theo quần thể, thám hiểm
khoảng không đáp án ngẫu nhiên bằng cách sử dụng một vài đáp án đề cử thay cho
cách ước lượng đáp án đơn lẻ được sử dụng trong nhiều kỹ thuật cổ điển. Sự thành
công của thuật toán này nằm ở khả năng tìm kiếm đáp án theo cách thám hiểm ngẫu
3
nhiên trong khu vực khả thi chứ không phải thám hiểm toàn bộ khu vực. Kết quả
tìm được theo cách này nhanh hơn, tiêu tốn ít tài nguyên phần cứng máy tính hơn
mà vẫn cho được khả năng tối ưu toàn bộ. Một vài kỹ thuật tiến hóa đã được phát
triển trong lĩnh vực tính toán tối ưu hóa mà phổ biến nhất là các kỹ thuật: Genetic
Algorithm (Thuật toán Di truyền) [5, 6], Differential Evolution (Thuật toán tiến
hóa) [7], Ant Colony Optimization (Tối ưu hóa đàn kiến) [8], Interior Point
Methods (Phương pháp điểm nội) [9], Particle Swarm Optimization (Thuật toán bầy
đàn) [13] … Trong sự phát triển của trí tuệ nhân tạo, gần đây trong lĩnh vực công
nghệ thông tin xuất hiện thuật toán Gravitational Search Algorithm [10], đây là
thuật toán có tuổi đời khá non trẻ nhưng đã được ứng dụng vào một số các lĩnh vực
nghiên cứu, một trong những lĩnh vực ứng dụng của GSA là lĩnh vực hệ thống điện.
Một số nhà khoa học trên thế giới đã triển khai đưa thuật toán GSA vào ứng dụng
giải bài toán điều phối kinh tế trong hệ thống điện và đã cho ra những kết quả tốt.
1.2 TÓM TẮT MỘT SỐ TÀI LIỆU VÀ BÀI BÁO LIÊN QUAN
Nhóm Esmat Rashedi
- Bài báo: “GSA: A Gravitational Search Algorithm”. Tác giả Esmat
Rashedi, Hossein Nezamabadi-pour, Saeid Saryazdi. [10]
Trong bài báo này, một thuật toán tối ưu hóa mới dựa trên định luật hấp dẫn
và khối lượng tương tác được giới thiệu. Trong các thuật toán đề xuất, các đại lý tìm
kiếm là một tập hợp quần chúng tương tác với nhau dựa trên các lực hấp dẫn
Newton và định luật chuyển động. Các phương pháp được đề xuất đã được so sánh
với một số phương pháp tìm kiếm heuristic nổi tiếng. Các kết quả thu được xác
nhận hiệu suất cao của phương pháp được đề xuất trong việc giải quyết các chức
năng phi tuyến khác nhau.
Nhóm Norlina Mohd Sabri, Mazidah Puteh, and Mohamad Rusop
Mahmood
- Bài báo: “A Review of Gravitational Search Algorithm”. Tác giả
Norlina Mohd Sabri, Mazidah Puteh, and Mohamad Rusop Mahmood. [11]
Bài viết này nhằm mục đích để khám phá thuật toán GSA và xác định thuật
toán đã được cải tiến như thế nào cho đến nay, các nghiên cứu và phát triển đã được
thực hiện kể từ khi có sự ra đời của thuật toán. Mục tiêu của bài viết này là để phân
4
tích các công trình liên quan đến GSA, để xem xét tiến bộ GSA và hiệu suất của
thuật toán này, xem xét các ứng dụng và để đưa ra những thách thức và tính khả thi
trong tương lai.
- Bài báo: “Gravitational Search Algorithm for Optimal Economics
Dispatch”. Tác giả P.K.Swain, N.C.Sahu, P.K.Hota. [12]
Bài báo này trình bày một phương pháp tối ưu hóa mới để tìm kiếm lời giải
tối ưu cho bài toán điều độ kinh tế (ED) bằng phương pháp sử dụng thuật toán tìm
kiếm hấp dẫn (GSA). Điều phối kinh tế xác định năng lượng điện được tạo ra bởi
các đơn vị phát điện đã cam kết trong một hệ thống điện do đó chi phí trong hệ
thống được giảm thiểu trong khi đáp ứng nhu cầu các phụ tải. Bài viết này trình bày
một thuật toán mới dựa trên định luật hấp dẫn và tương tác khối lượng để giải quyết
kinh tế vấn đề tải công văn (ED) bởi một thuật toán tối ưu hóa mới gọi là thuật toán
Tìm kiếm hấp dẫn (GSA). Các kết quả mô phỏng cho thấy kỹ thuật này được thực
hiện dễ dàng, hội tụ với thời gian thực hiện ít hơn và giải pháp rất tối ưu cho bài
toán điều độ kinh tế với chi phí tối thiểu mà hệ thống có thể đạt được. Mô phỏng kết
quả đã được thực hiện trên các hệ thống điện khác nhau với số lượng nguồn phát
khác nhau và so sánh với các phương pháp tiếp cận phổ biến khác. Những phát hiện
này đã khẳng định sự vững mạnh, hội tụ nhanh so với phương pháp đề xuất trên các
kỹ thuật hiện có khác.
1.3 HƯỚNG TIẾP CẬN ĐỀ TÀI
Thực tế, khi giải quyết bất kỳ bài toán nào người giải đều mong muốn có
phương án tốt nhất theo một hoặc một vài tiêu chí nào đó như tiết kiệm thời gian
nhất, chi phí nhỏ nhất, năng suất lớn nhất, quãng đường đi ngắn nhất, thiết kế kết
cấu với trọng lượng vật liệu nhỏ nhất…. Tuy nhiên trong phần lớn các bài toán tối
ưu, người sử dụng thường có băn khoăn đó là: Kết quả nhận được từ quá trình tính
toán đã thật sự là phương án tốt nhất chưa? Vì vậy, việc phát triển các thuật toán tối
ưu đủ mạnh luôn được người làm kỹ thuật quan tâm. Các phương pháp cổ điển
trước đây chỉ thích hợp được với những bài toán tìm kiếm tối ưu trong không gian
tìm kiếm nhỏ nên để tìm được nghiệm tối ưu trong không gian tìm kiếm lớn cần
phải có những phương pháp tìm kiếm đặc biệt hơn mà nổi bật trong số đó là thuật
toán Gravitational Search Algorithm (GSA). Bằng việc sử dụng phần mềm
5
MATLAB, Thuật toán GSA vẫn đang được quan tâm, phát triển và ứng dụng vào
giải các bài toán tối ưu hóa hiện nay.
Bài toán điều phối kinh tế tối ưu trong hệ thống điện ED đã được giải bằng
nhiều phương pháp khác nhau như đã nêu trong những bài báo ở trên. Qua một số
bài báo ở trên cho thấy các nhà khoa học trên thế giới đã ứng dụng thuật toán GSA
vào trong bài toán hệ thống điện đơn giản và đã cho ra những kết quả đạt hiệu suất
cao. Thuật toán GSA có tuổi đời còn khá non trẻ: ra đời vào năm 2009 [10] hiện
được áp dụng trong một vài lĩnh vực như: Công nghệ thông tin, Hệ thống điện,… vì
vậy thuật toán này đang được nhiều nhà khoa học quan tâm nghiên cứu đi sâu.
Giải thuật GSA có ưu điểm là đơn giản, ổn định và có khả năng thích ứng
nên có thể ứng dụng trong nhiều bài toán tối ưu. Khả năng tìm kiếm toàn cục của
GSA cũng tốt hơn các giải thuật nổi tiếng khác trong hầu hết các trường hợp
1.4 MỤC TIÊU CỦA ĐỀ TÀI
Ứng dụng thuật toán Gravitional Search Algorithrm vào giải bài toán điều
phối kinh tế công suất (ED) giữa các nhà máy điện.
Áp dụng thuật toán GSA giải mạng điện hàm chi phí nhiên liệu có xét ảnh
hưởng của điểm van công suất nhằm tạo ra chương trình có thể tính toán tốt hơn,
nhanh hơn khi so sánh với kết quả tính toán của các thuật toán khác. Đối chiếu và
đánh giá kết quả tính toán được bằng thuật toán GSA với một số thuật toán khác.
1.5 ĐIỂM MỚI CỦA ĐỀ TÀI
Xây dựng giải thuật GSA giải bài toán điều phối kinh tế công suất ED.
Lập trình tính toán điều phối tối ưu trong mạng điện.
Nhận xét, đánh giá kết quả thu được, so sánh với kết quả dùng các giải thuật
khác đã công bố trên tạp chí khoa học trên thế giới.
1.6 KẾT CẤU ĐỀ TÀI
Tổng quan
Giới thiệu bài toán phân bố công suất tối ưu ED.
Giới thiệu về thuật toán GSA.
Ứng dụng thuật toán GSA vào giải bài toán phân bố công suất tối ưu trong hệ
thống điện.
Tổng kết và hướng phát triển đề tài.
6
CHƯƠNG 2 GIỚI THIỆU BÀI TOÁN ĐIỀU PHỐI KINH TẾ
TRONG HỆ THỐNG ĐIỆN
2.1 THUẬT TOÁN TỐI ƯU:
Trong kỹ thuật, khi giải quyết bất kỳ nhiệm vụ nào chúng ta đều mong muốn
có phương án tốt nhất theo một hoặc một vài tiêu chí nào đó. Có thể liệt kê rất nhiều
những ví dụ cụ thể như: tiết kiệm thời gian nhất, chi phí nhỏ nhất, năng suất lớn
nhất, quãng đường đi ngắn nhất, thiết kế kết cấu với trọng lượng vật liệu nhỏ nhất…
Để giải được những bài toán này, toán học đã cho ra đời một ngành là “Quy hoạch
toán học” hay “tối ưu hóa” [14], [15].
Bài toán tối ưu nói chung được viết dưới dạng toán học như sau:
n
f(x) min(max); x R
Tìm giá trị cực tiểu (hoặc cực đại) hàm:
(2.1)
Với các điều kiện:
gi(x) ≥ 0; i = 1,2,... , m
(2.2)
hi(x) = 0; i = 1,2,... , l (2.3)
Bài toán đặt ra yêu cầu là tìm tập hợp các biến xi, i = 1, … ,n thoả mãn các
điều kiện ràng buộc đồng thời hàm f(x) đạt giá trị cực tiểu (hoặc cực đại).
Hàm f(x) trong biểu thức (2.1) được gọi là hàm mục tiêu hoặc tiêu chuẩn tối
ưu, biểu diễn mối quan hệ giữa tiêu chuẩn chất lượng của quá trình khảo sát và các
biến độc lập x.
Các hàm số gi(x), hi(x) là các điều kiện ràng buộc của bài toán tối ưu dưới
dạng đẳng thức và bất đẳng thức. Trong không gian các biến, các hàm số này tạo ra
miền giới hạn D các khả năng cho phép của hàm f(x).
Nếu như D Rn (với R là số chiều của hàm mục tiêu), có nghĩa là không tồn
tại bất kỳ một điều kiện giới hạn nào ta nói rằng bài toán quy hoạch phi tuyến không
có điều kiện ràng buộc.
Tuy nhiên trong phần lớn các bài toán tối ưu, người sử dụng thường có băn
khoăn đó là: kết quả nhận được từ quá trình tính đã thật sự là phương án tốt nhất
7
chưa. Để minh họa vấn đề này ta có thể xét ví dụ như hàm Peaks (2.4) - hai biến là
1
2
2 (-x -(x +1) ) 2
2 1
2 2 (-x -x ) 1 2
2 2 (-(x +1) -x ) 2
1
f(x) = 3.(1- x ) .e
-10.(
- x - x ).e
-
.e
1
3 1
5 2
x 1 5
3
hàm đơn điệu đa cực trị, được biểu diễn bằng đồ thị như trên Hình 2.1.
(2.4)
Hình 2.1: Minh hoạ tối ưu toàn cục hàm Peaks.
Như trên hình 2.1, xung quanh phương án tốt nhất (ở đây chọn là điểm thấp
nhất - điểm A) còn có một điểm đạt cực trị địa phương là điểm B và một số điểm
nghi ngờ có cực tiểu khác. Trong quá trình giải, rất có khả năng kết quả giải bài
toán tối ưu của chúng ta bị "kẹt" tại một cực trị nào đó (không phải điểm A) và
không thoát ra được. Vì vậy, việc phát triển các thuật toán đủ mạnh tiệm cận được
giá trị tối ưu luôn được người làm kỹ thuật quan tâm.
2.2 BÀI TOÁN ĐIỀU PHỐI KINH TẾ CỔ ĐIỂN:
Để hệ thống điện hoạt động hiệu quả và tin cậy thì một số kỹ thuật đã được
phát triển để tính toán xác định dự báo công suất và mức công suất phát. Điều phối
công suất là một trong các kỹ thuật trên để điều chỉnh biến điều khiển và phân phối
công suất cho hệ thống điện hoạt động tối ưu. Điều phối công suất có hai cách: điều
phối công suất thực và điều phối công suất phản kháng. Bài toán điều phối kinh tế
tìm điểm hoạt động tối ưu để phân phối công suất thực giữa các nhà máy nhằm
giảm thấp nhất chi phí sản xuất. Điều phối công suất phản kháng dùng để cực tiểu
tổn thất hệ thống, nâng cao hiệu suất và khả năng tận dụng nguồn.
8
Bài toán điều phối công suất làm cải thiện việc hoạt động ổn định của hệ
thống điện. Thường làm giảm mô hình hệ thống điện, làm đơn giản các giải pháp
chi phí về chất lượng. Việc sử dụng đúng đắn và chính xác hơn các mô hình sản
lượng điện làm cho lời giải bài toán tốt hơn nhưng vấn đề khó khăn cũng tăng lên
đáng kể.
Mô hình phổ biến cải tiến bài toán điều phối kinh tế bao gồm: hàm chi phí có
xét ảnh hưởng của điểm van công suất, vùng hoạt động không liên tục và sự chuyển
đổi các loại nhiên liệu; các loại ràng buộc an ninh hệ thống điện như giới hạn dòng
công suất, dự trữ công suất máy phát và cấu hình điện áp. Trong chương này chúng
tôi trình bày hệ thống các biểu thức của bài toán điều phối kinh tế với hàm chi phí
trơn dạng bậc hai cổ điển và hàm chi phí có xét ảnh hưởng của điểm van công suất.
Hàm mục tiêu của bài toán điều phối kinh tế cổ điển là cực tiểu tổng chi phí
hệ thống điện (1) bằng cách hiệu chỉnh công suất phát của mỗi nhà máy kết nối với
lưới điện. Tổng chi phí được biểu diễn bằng hàm tổng của các chi phí ở mỗi nhà
máy.
Xét một hệ thống có N nhà máy, mỗi nhà máy đảm nhận một lượng công
suất Pi MW. Các nhà máy nên phát công suất sao cho tổng chi phí nhiên liệu F là
N
G
min
)
nhỏ nhất.
F P ( i G i
i
1
(
)
(2.5)
GF P là hàm chi phí của nhà máy thứ i,
i
iGP là công suất thực
i
Trong đó
phát ra của nhà máy thứ i và NG là tổng số lượng các nhà máy kết nối với hệ thống
điện.
9
Hình 2.2: Đường cong chi phí phổ biến của nhà máy nhiệt điện.
Mỗi hàm chi phí của nhà máy thiết lập mối quan hệ giữa nhà máy và hệ
thống thông qua khả năng phát công suất với chi phí phát của nhà máy. Thông
thường các nhà máy được mô hình bằng hàm chi phí trơn như trong (2.6) để đơn
)
a
giản bài toán tối ưu và khả năng ứng dụng các kỹ thuật truyền thống để tính toán.
i
F P ( i G i
b P i G i
2 c P i G i
(2.6)
Trong đó ai, bi, ci là hệ số chi phí của hàm chi phí nhà máy thứ i.
2.2.1 Ràng buộc đẳng thức:
Ràng buộc cân bằng công suất: Ràng buộc cân bằng công suất là ràng buộc
đẳng thức dùng để giảm bớt công suất hệ thống dựa trên nguyên lý cơ bản cân bằng
giữa tổng công suất nhà máy phát với tổng tải của hệ thống. Cân bằng chỉ xảy ra khi
iGP bằng với tổng tải trong hệ thống PD cộng thêm
tổng công suất nhà máy phát
N
G
một lượng tổn hao PL được biểu diễn như trong (2.7).
P D
P L
P G i
i
1
(2.7)
Tổn thất trong hệ thống có thể xác định một cách chính xác nhờ phương
pháp phân luồng công suất. Một cách điển hình để ước lượng tổn thất bằng cách mô
hình chúng dạng hàm của hệ thống nhà máy phát sử dụng công thức tổn thất của
Kron (2.8). Một số cách khác để mô hình hóa tổn thất là sử dụng hệ số phạt hoặc
xem tổn thất là hằng số.
N
N
N
G
G
G
B
P L
0
00
P B P G ij G i i
P B G i i
j
1
i
1
j
1
10
(2.8)
Trong đó Bij, Bi0, B00 gọi là tổn thất hay hệ số B.
2.2.2 Ràng buộc bất đẳng thức:
min
2.2.2.1 Giới hạn công suất thực phát ra:
iGP và
max
Giới hạn công suất thực phát ra: Mỗi nhà máy có giới hạn thấp nhất
iGP phát công suất vì nó phụ thuộc vào cấu trúc của máy phát.
giới hạn cao nhất
Các giới hạn trên được định nghĩa bằng một cặp của ràng buộc bất đẳng thức (2.9).
min P G i
P G i
max P G i
(2.9) , i = 1,…, NG
2.2.2.2 Giới hạn tốc độ:
Công suất của máy phát không thể tăng hay giảm đến giá trị bất kỳ tức thời.
Tầm hoạt động của tất cả các máy phát được giới hạn bởi độ dốc tốc độ được trình
0
DR
bày như sau:
P P i i 0 P P UR i i i
(2.10)
0
iP : công suất phát trước đó của tổ máy thứ i
Trong đó:
DRi và URi : là giới hạn độ dốc đi xuống và đi lên
2.2.2.3 Giới hạn về vùng cấm vận hành:
Vùng cấm của một nút có thể được trình bày qua công thức sau:
(2.11)
i kP ,
L ,
U ,
i kP là giá trị giới hạn dưới và giới hạn trên của vùng cấm
Trong đó:
thứ k của nút thứ i ; k là hệ số của vùng cấm.
11
Hình 2.3: Đồ thị biểu diễn vùng cấm của nút nhiệt cơ bản.
2.2.3 Ràng buộc về công suất truyền tải:
P
k
1, 2,...,
L
Ràng buộc về công suất truyền tải của đường dây được trình bày như sau:
Lf k ,
max P Lf k
, ;
(2.12)
PLf,k : công suất thực của đường dây thứ k và L là số lượng đường dây truyền
tải.
2.2.4 Bài toán điều phối kinh tế với hàm chi phí nhiên liệu không trơn:
Các nhà máy phát thường được mô hình hóa sử sụng hàm chi phí trơn để
biểu diễn mối quan hệ giữa công suất phát ra và chi phí sản xuất. Hàm chi phí loại
này có ưu điểm là làm đơn giản bài toán điều phối kinh tế và khả năng sử dụng
nhiều kỹ thuật áp dụng vào để giải bài toán này. Trong một số trường hợp, biểu diễn
dưới dạng bậc hai không mô hình hết được đặc điểm của nhà máy điện, do đó cần
mô hình chính xác hơn để cho kết quả tốt hơn trong việc giải bài toán điều phối
kinh tế. Mô hình chính xác hơn thường có dạng hàm phi tuyến hơn, không trơn và
nằm trong miền lõm. Một số ví dụ của hàm chi phí không trơn là: hàm chi phí có
xét ảnh hưởng điểm van công suất, hàm bậc hai liên tục từng khúc gồm hàm có
nhiều loại nhiên liệu và hàm có vùng hoạt động không liên tục. Trong đó, hàm chi
phí nhiên liệu có xét ảnh hưởng điểm van công suất được sử dụng khá phổ biến trên
thế giới.
12
2.2.4.1 Bài toán điều phối kinh tế có điểm van công suất:
Nhà máy điện thường sử dụng nhiều van để điều khiển công suất phát của
nhà máy [16-21]. Trong giai đoạn đầu khi van nạp hơi nước được mở trong nhà
máy nhiệt điện, chi phí do tổn hao gia tăng một cách đột ngột làm cho hàm chi phí
có độ nhấp nhô như hình 2. Hiệu ứng này được gọi là điểm van công suất. Loại bài
toán này vô cùng khó giải quyết với những kỹ thuật thông thường bởi vì tồn tại sự
thay đổi đột ngột và không liên tục trong sự gia tăng của hàm chi phí.
Hình 2.4 Hàm chi phí nhiên liệu của nhà máy nhiệt điện với điểm van công suất.
2.2.4.2 Biểu thức điều phối kinh tế với điểm van công suất:
Điều phối kinh tế với điểm van công suất dùng để cực tiểu chi phí hệ thống
(2.8) dựa trên hàm chi phí có xét ảnh hưởng của vị trí van. Vị trí van công suất
thường được mô hình bằng cách thêm hàm sin vào hàm chi phí bậc hai cổ điển
)
d
sin
a i
i
F P ( i G i
b P i G i
2 c P i G i
P G i
min e P i G i
(2.13).
(2.13)
Trong đó ai, bi, ci, di và ei là hệ số chi phí của nhà máy thứ i.
Biểu thức cơ bản của bài toán này là các vấn đề ràng buộc cân bằng công
suất (2.9) và giới hạn máy phát (2.10). Những ràng buộc khác có thể thêm vào tùy
thuộc vào mô hình yêu cầu. Bài toán điều phối kinh tế với điểm van công suất đã
được một số nhà khoa học nghiên cứu. Sheblé và Walters [18] sử dụng GA để giải
bài toán này. Ngoài ra, K. Wong và Y. Wong đã đề xuất cách giải bài toán điều phối
13
kinh tế với điểm van công suất sử dụng GA và giải thuật luyện kim SA (Simulated
Annealing). K. Wong cùng với B. Lau và A. Fry [20] đã trình bày phương pháp
dùng mạng noron giải bài toán điều phối kinh tế có xét ảnh hưởng của điểm van
công suất. Nguyên lý cơ bản của điểm van công suất và hàm chi phí của nó được
trình bày trong [16, 17].
2.3 TỔNG QUAN CÁC PHƯƠNG PHÁP ĐÃ ĐƯỢC ÁP DỤNG ĐỂ GIẢI BÀI
TOÁN ĐIỀU PHỐI KINH TẾ (ED):
2.3.1 DE ( Differential Evolution ) [22][23]
DE là thuật toán tiến hóa được đề xuất bởi Storn và Price. Phương pháp này
rất hiệu quả trong những bài toán tối ưu không tuyến tính với nhiều ràng buộc. Ưu
điểm của DE so với các phương pháp tiến hóa khác là cấu trúc đơn giản và gọn, ít
thông số điều khiển, điểm hội tụ cao. Do khả năng dò tìm tin cậy và mạnh nên nó
được ứng dụng trong nhiều lĩnh vực khác nhau như cơ khí, công nghệ sinh học, điện
tử… Đặc biệt là ứng dụng để giải bài toán ED trong các trường hợp khác nhau. Với
bài toán ED có xét đến điểm van công suất, phương pháp này cho các kết quả tốt
hơn so với các phương pháp khác từ trước đến giờ. Sự mạnh mẽ của phương pháp
được chứng minh bởi việc thay đổi các nhu cầu phụ tải của bài toán. Với bài toán
ED có hàm chi phí bậc hai theo từng đoạn, DE cho ra kết quả tốt hơn các phương
pháp truyền thống khác. Đối với các bài toán có xét đến tổn hao thì DE cho ra kết
quả tốt hơn nhiều so với các phương pháp truyền thống khác.
2.3.2 PSO ( Particle Swarm Optimization ) [24][25][26]
Thuật toán PSO được giới thiệu đầu tiên bởi Eberhart và Kennedy. Thuật
toán này phỏng đoán hành vi của đàn chim hay đàn cá tìm thức ăn. Mỗi cá thể thay
đổi vận tốc và vị trí dựa trên kinh nghiệm của bản thân và kinh nghiệm của cả quần
thể. Thuật toán PSO được áp dụng cho nhiều bài toán ELD khác nhau như bài toán
điều độ kinh tế đa vùng, bài toán điều độ kinh tế với hàm chi phí gồm nhiều đoạn
bậc thang… Phương pháp PSO hội tụ tới điểm tối ưu toàn cục hay gần toàn cục bất
kể hình dạng của hàm chi phí như hàm chi phí không liên tục, hàm chi phí lồi…
Hiệu quả tính toán và đặc tính hội tụ của phương pháp PSO rất tốt và nó được áp
dụng rộng rãi trong nhiều bài toán tối ưu.
2.3.3 ABC (Thuật toán Artificial Bee Colony) [27][28][29]
14
ABC là kỹ thuật thông minh lấy cảm hứng từ hành vi tìm kiếm thức ăn của
loài ong mật. Lý thuyết này lần đầu tiên được phát triển bởi Karaboga năm 2005.
Trong tự nhiên bầy ong mật tìm kiếm thức ăn theo quy trình sau: Đầu tiên các ong
do thám sẽ được cử đi thăm dò các vùng thức ăn tiềm năng; các ong do thám sẽ di
chuyển ngẫu nhiên từ bụi hoa này sang bụi hoa khác; sau đó chúng quay về tổ và
thông tin cho cả bầy về hướng di chuyển, khoảng cách từ tổ đến vùng thức ăn và
chất lượng của các vùng thức ăn; những thông tin này sẽ làm cho bầy ong bay đến
các vùng thức ăn nhanh chóng và chính xác hơn. Một vấn đề tự nhiên nhưng là
điểm trọng tâm cho ý tưởng thuật toán bầy ong là: Nơi nào có thức ăn dồi dào hơn
thì nơi đó sẽ có nhiều ong được cử đến hơn.
2.3.4 HNN ( Hopfield Neuron Network ) [30][31]
HNN là một phương pháp dựa trên mạng Neuron Hopfield để giải bài toán
ED. Phương pháp này rất thành công trong việc giải bài toán ED với hàm chi phí
bậc hai theo từng đoạn. HNN có thể áp dụng dễ dàng trong trường hợp số tổ máy
lớn. Tuy nhiên, hạn chế của HNN là tốc độ hội tụ của nó rất chậm. Ngoài ra, số
vòng lặp lớn và sự dao động là nhược điểm của HNN khi giải các bài toán tối ưu.
2.3.5 ELANN ( Enhanced Lagrangian Artificial Neural Network) [32]
ELANN được áp dụng để giải các bài toán ED với hàm chi phí gồm nhiều
đoạn bậc hai và các ràng buộc không tuyến tính. Tốc độ hội tụ của ELANN tăng lên
bởi việc áp dụng phương pháp động lượng và việc cung cấp các tiêu chuẫn cho sự
lựa chọn các thông số tốc độ. Nhược điểm của ELANN là có số vòng lặp lớn và
thường dao động trong suốt quá trình quá độ.
2.3.6 HS (Harmony Search) [33]
HS được đề xuất bởi Zong Woo Geem, lấy cảm hứng từ việc nghe nhạc Jazz
để sáng tạo nên một thuật toán tối ưu. Việc nhạc sĩ ngẫu hứng chơi một khúc nhạc,
bằng cảm quan âm nhạc nhạc sĩ ấy tìm thấy sự hòa hợp của các nốt nhạc và sáng tác
nên một khác nhạc cũng giống như ngẫu nhiên tạo ra một phương án cho một bài
toán rồi tìm kiếm trạng thái tốt nhất trong số các phương án đó. Phương pháp HS đã
tìm ra mối tương đồng thú vị và đáng kinh ngạc này.
15
2.3.7 CS (Cuckoo Search) [34][35][36][37]
Bài báo trình bày việc áp dụng thuật toán Cuckoo Search vào giải bài toán
xác định vị trí đặt máy phát phân tán DG trong mạng điện để cải thiện chất lượng
điện áp và giảm tổn thất công suất trong mạng điện phân phối. Kết quả tính toán của
thuật toán CS so sánh với các thuật toán tiến hóa khác như: GA, PSO trong nhiều
trường hợp vận hành của Hệ thống lưới điện. Kết quả tính toán cho thấy thuật toán
CS có nhiều ưu điểm như: số biến điều khiển ít hơn các thuật toán khác, kết quả tính
toán tốt hơn, đặc tính hội tụ nhanh hơn và vững mạnh hơn các thuật toán khác.
16
CHƯƠNG 3 THUẬT TOÁN GRAVITAIONAL SEARCH
ALGORITHM
Trong những năm gần đây nhiều phương pháp tối ưu được nghiên cứu, phát
triển. Trong đó, có nhóm phương pháp lấy cảm hứng từ hành vi bầy đàn trong tự
nhiên. Các phương pháp này là khá phù hợp để giải các bài toán tối ưu phức tạp mà
các phương pháp tối ưu truyền thống gặp nhiều khó khăn vì số chiều của không
gian bài toán lớn, hoặc vùng xác định không liên tục, hàm mục tiêu không liên tục,
có nhiều điểm cực trị địa phương... Trong các phương pháp thuộc nhóm này thì
Gravitational Search Algorithm (GSA) là một trong những giải thuật tối ưu mới
nhất được Rashedi đề xuất năm 2009. Như trong [10], các tác giả đã khẳng định,
không một giải thuật tối ưu nào là tốt hơn hẳn các phương pháp khác, nên trong
chương này sẽ tìm hiểu giải thuật GSA và GSA cải tiến để áp dụng giải điều phối
kinh tế.
Người ta có thể nhận ra hai qui luật phổ biến trong thuật toán: thăm dò và
khai thác. Các thăm dò là khả năng mở rộng không gian tìm kiếm, nơi khai thác là
khả năng của việc tìm kiếm một giải pháp tốt. Trong lần lặp đầu, một thuật toán tìm
kiếm tiến hóa khám phá không gian tìm kiếm để tìm giải pháp mới. Để tránh bẫy
trong một tối ưu hóa địa phương, các thuật toán phải sử dụng các thăm dò trong vài
lần lặp đầu tiên. Do đó, việc thăm dò là một vấn đề quan trọng trong một thuật toán
heuristic dựa vào cá thể. Để có một tìm kiếm hiệu suất cao, một chìa khóa quan
trọng là một sự cân bằng thích hợp giữa thăm dò và khai thác. Nói cách khác, tất cả
các thuật toán tìm kiếm có một khuôn khổ chung [10].
3.1 THUẬT TOÁN GSA CỔ ĐIỂN:
Bài toán điều phối kinh tế (ED) liên quan đến việc lập kế hoạch phát tối ưu
của các nhà máy phát điện trong một hệ thống điện nhằm mục đích giảm thiểu tổng
chi phí nhiên liệu trong khi đáp ứng nhu cầu phụ tải và hạn chế hoạt động. ED đóng
một vai trò quan trọng trong việc lập kế hoạch và kiểm soát việc hoạt động của hệ
thống điện hiện đại Trong vài năm qua, một số phương pháp lập trình toán học cổ
điển đã được phát triển để giải quyết ED. Các phương pháp tối ưu cổ điển rất nhạy
cảm với điểm khởi đầu và thường xuyên hội tụ về với giải pháp tối ưu hóa địa
phương hoặc phân tách hoàn toàn. Phương pháp lập trình tuyến tính nhanh chóng và
17
đáng tin cậy nhưng có nhược điểm là cho ra các chi phí gần đúng từng phần theo
tuyến tính [17]. Phương pháp lập trình phi tuyến liên quan đến vấn đề hội tụ và độ
phức tạp của thuật toán. Gần đây, để các phương pháp thuận tiện hơn cho việc giải
quyết các bài toán tối ưu hóa hiện đại , bài toán ED được xem như một bài toán tối
ưu hóa không trơn. Một thuật toán tìm kiếm heuristic mới, cụ thể là thuật toán tìm
kiếm hấp dẫn (GSA) dựa trên luật hấp dẫn và các định luật về chuyển động đã được
đề xuất bởi Rashedi [10]. Mục tiêu chính của nghiên cứu này nhằm để trình bày
việc sử dụng các kỹ thuật tối ưu hóa GSA cho hoạt động kinh tế của hệ thống điện.
Giải thuật tối ưu GSA lần đầu tiên được Rashedi đề xuất năm 2009 [10], giải
thuật này lấy cảm hứng từ luật trọng trường của Newton. Trong giải thuật này, các
cá thể có đặc tính của chúng được đo bởi các khối lượng. Tất cả các vật thể hấp dẫn
nhau bởi lực trọng trường, lực này tạo ra chuyển động của tất cả vật thể và các vật
thể có xu hướng dịch chuyển đến các vật thể có khối lượng nặng hơn. Các vị trí của
vật thể có khối lượng nặng được xem là những lời giải tốt của bài toán [12].
Trong GSA, mỗi vật thể được đặc trưng bởi 4 thông số: vị trí, khối lượng
quán tính, khối lượng trọng trường chủ động, khối lượng trọng trường thụ động.
Trong đó, vị trí của vật thể tương ứng với lời giải của bài toán, khối lượng quán tính
và trọng trường được xác định thông qua tính toán hàm mục tiêu của vật thể đó. Nói
cách khác là mỗi vật thể đại diện cho một lời giải và sau mỗi lần lặp giải thuật sẽ
điều chỉnh khối lượng của các vật thể, các vật thể sẽ bị hấp dẫn bởi vật có khối
lượng nặng nhất. Vật thể nặng nhất này sẽ biểu thị lời giải tối ưu trong không gian
tìm kiếm [10].
Trong giải thuật GSA, các vật thể được xem như trong một hệ cách ly, trong
đó các vật thể chuyển động theo luật trọng trường và luật chuyển động của Newton.
Cụ thể:
+ Luật trọng trường: Mỗi phần tử hấp dẫn các phần tử khác và lực trọng
trường giữa chúng là tỷ lệ thuận với tích của các khối lượng và tỷ lệ nghịch với khoảng cách giữa chúng,R. Các tác giả đề nghị sử dụng ở đây là R thay vì R2, bởi vì theo những kết quả thí nghiệm, R tạo ra kết quả tốt hơn R2 trong tất cả các trường
hợp thí nghiệm [10].
18
+Luật chuyển động: Vận tốc tức thời của bất cứ vật thể nào bằng tổng của
vận tốc trước đó và sự thay đổi vận tốc. Sự thay đổi vận tốc hay gia tốc của vật thể
bằng lực tác động lên hệ chia cho khối lượng quán tính [10].
Bây giờ, xét một hệ có N phần tử (agents). Chúng ta định nghĩa vị trí
X
của phần tử thứ i là:
i
n (x , x ..., x ,..., x ) i
2 i
d i
1 i
d
với i = 1,2,…,N (3.1)
ix là vị trí của vật thể i theo chiều d.
Trong đó
M (t)xM (t)
(x (t) x (t))
Tại thời điểm t, ta định nghĩa lực tác dụng lên vật thể i từ vật thể j là:
d F (t) G(t). ij
d i
d j
aj
pi R (t) ij
(3.2)
Trong đó:
+ Maj là khối lượng hấp dẫn chủ động của vật thể j
+ Mpi là khối lượng hấp dẫn thụ động của vật thể i
+ G(t) là hằng số trọng trường tại thời điểm t
+ ε là hằng số
+ Rij(t) là khoảng cách Euclide của vật thể i và j:
R ( t) ij
|| X ( t), X ( t) || 2
i
j
(3.3)
Để cho giải thuật có đặc tính ngẫu nhiên, giả thiết rằng lực tổng hợp (total
force) tương tác lên vật thể i theo chiều d là một tổng ngẫu nhiên của các thành
N
phần lực theo chiều d của các vật thể:
d F (t) i
d ij
j 1, j
i
(3.4) rand F (t) j
Trong đó: randj là một số ngẫu nhiên trong đoạn [0, 1].
d
Theo luật chuyển động, gia tốc của vật thể i, theo hướng d, ở thời điểm t là
d a (t) i
ia ( t ) :
d F (t) i M (t) ii
(3.5)
Trong đó: Mii là khối lượng quán tính (inertial mass) của vật thể i.
Ngoài ra, vận tốc thời điểm tiếp theo của một vật thể là một phần của vận tốc
hiện tại cộng thêm gia tốc. Cho nên, vị trí và vận tốc của vật thể ở thời điểm tiếp
theo được tính như sau:
v ( t 1)
19
d i
rand .v (t) i
d i
d a (t) i
x (t 1)
v (t 1)
(3.6)
d i
d x (t) i
d i
(3.7)
Trong đó: randi là một biến ngẫu nhiên chuẩn (uniform random variable)
trong đoạn [0,1].
Hình 3.1 Lực và gia tốc tương tác lên vật thể 1 do các vật thể khác sinh ra [10].
Thuật toán tìm kiếm dựa vào số lượng cá thể phải vượt qua ba bước trong
mỗi lần lặp để nhận ra các khái niệm về thăm dò và khai thác: tự thích ứng, hợp tác
và cạnh tranh. Trong tự thích ứng bước, mỗi thành viên (agent) cải thiện hiệu quả
của nó. Trong bước hợp tác, các cá thể cộng tác với nhau bằng cách truyền tải thông
tin. Cuối cùng, ở bước cạnh tranh, các thành viên cạnh tranh để tồn tại. Những bước
này có các hình thức thường ngẫu nhiên, và có thể được thực hiện theo những cách
khác nhau [10].
Hằng số trọng trường G, được khởi tạo lúc bắt đầu và sẽ bị giảm dần theo
thời gian để điều khiển độ chính xác tìm kiếm (search accurary). Nói cách khác, G
là một hàm của giá trị ban đầu G0 và thời gian t :
(3.8) G(t) = G(G0, t).
Khối lượng quán tính và trọng trường được tính toán đơn giản nhờ đánh giá
giá trị hàm mục tiêu (fitness evaluation). Khối lượng càng nặng nghĩa là vật thể
càng hiệu quả trong việc làm lời giải. Điều này nghĩa là khối lượng tốt hơn có lực
hấp dẫn các vật khác lớn hơn và di chuyển chậm hơn. Giả thiết các khối lượng quán
20
tính và trọng trường bằng nhau, giá trị của các khối lượng được tính thông qua giá
trị hàm mục tiêu. Chúng ta cập nhật các khối lượng quán tính và trọng trường bởi
phương trình sau:
i
Mai=Mpi=Mii=Mi, i=1,2,…N (3.9)
m (t) i
fit (t) worst(t) best(t) worst(t)
(3.10)
M (t) i
m (t) i N
m (t) j
j 1
(3.11)
Trong đó:
+ fiti(t) là giá trị fitness của vật thể i tại thời điểm t.
best(t) min fit (t), j 1..N
+ worst(t), best(t) được định nghĩa (trong bài toán tìm Min) như sau :
j
worst(t) maxfit (t), j 1..N
j
(3.12)
(3.13)
Để tránh việc kẹt ở cực trị địa phương thì ban đầu giải thuật sẽ thực hiện việc
mở rộng không gian tìm kiếm (exploration). Dần dần theo các lần lặp về sau giải
thuật giảm khả năng exploration và tăng khả năng phân tích - exploitation [10].
Để cải thiện hiệu suất của GSA nhờ điều chỉnh không gian tìm kiếm và tốc
độ tìm kiếm chỉ có những vật thể trong tập Kbest mới hấp dẫn những vật thể khác.
Số lượng các vật thể trong Kbest một phụ thuộc vào thời gian (số lần lặp), ban đầu
có giá trị khởi tạo K0 và giảm dần theo thời gian. Ví dụ, lúc đầu tất cả các vật thể
đều tác dụng lực lên nhau, và sau đó thì Kbest giảm dần dần, đến lần lặp cuối còn
lại một số vật thể nặng nhất tác dụng lực đến các vật thể khác [10]. Cho nên phương
N
trình (3.4) có thể được chỉnh sửa lại là:
d F (t) i
rand F (t) j
d ij
j Kbest, j i
(3.14)
Trong đó Kbest là tập của K phần tử đầu tiên có giá trị fitness tốt nhất và
khối lượng lớn nhất.
21
Hình 3.2 Nguyên lý của giải thuật GSA [10].
3.2 CÁC BƯỚC TRONG THUẬT TOÁN GSA:
Bước 1- Khởi tạo vị trí các vật thể
Khởi tạo một tập hợp các vật thể có vị trí ngẫu nhiên trong không gian bài
toán D chiều sử dụng hàm phân bố xác suất. Để làm điều này GSA xem hệ thống N
2 n X (x , x ..., x ,..., x ) i i
1 i
d i
i
phần tử , vị trí thứ i được định nghĩa như (3.15) sau:
với i = 1,2,…,N.
d là vị trí của vật thể i theo chiều d và D là số chiều không gian.
(3.15)
Ở đây xi
Bước 2: Đánh giá hàm mục tiêu của các vật thể.
Bước 3: Cập nhật G(t), best(t), worst(t), Mi(t) đối với i=1...N:
G(t) G(G , t)
Trong GSA, hằng số trọng trường G, lúc đầy lấy một giá trị khởi tạo G0, sau
0
đó sẽ được giảm dần theo thời gian:
(3.16)
Khối lượng của vật thể được tính toán sau khi tính được hàm mục tiêu của
tập hợp hiện tại:
i
22
m (t) i
fit (t) worst(t) best(t) worst(t)
M ( t ) i
m ( t ) i N
(3.17)
m ( t ) j
j 1
(3.18)
Ở đây, Mi(t) và fiti(t) là khối lượng và giá trị hàm mục tiêu của vật thể i tại
thời điểm t và
best(t) max fit (t), j 1..N
worst(t) và best(t) là giá trị kém nhất và tốt nhất (đối với bài toán maximization) thì:
j
worst(t) min fit (t), j 1..N
(3.19)
j
(3.20)
Bước 4: Tính toán lực tổng hợp theo các hướng khác nhau:
Lực tổng hợp của tập các vật thể có khối lượng nặng hơn sẽ tác dụng lên vật
M (t)xM (t)
(x (t) x (t))
d F (t) i
rand G(t). j
d j
d i
thể i được tính dựa vào lực trọng trường:
(3.21)
i
j Kbest, j
i
j R (t) ij
Trong đó:
+ randj là giá trị ngẫu nhiên từ [0,1]
+ ε là một số nhỏ và Rij(t) là khoảng cách Euclidian giữa 2 vật thể i và j
+ Kbest là tập hợp các phần tử có khả năng tác dụng lực lên các phần tử còn
lại.
Bước 5: Tính toán gia tốc và vận tốc:
(x (t) x (t))
Gia tốc của vật thể được tính nhờ luật chuyển động
d a (t) i
rand G(t). j
d j
d i
j Kbest, j
i
d F (t) i M (t) i
M (t) i R (t) ij
(3.22)
Trong đó randj là số ngẫu nhiên trong đoạn [0,1]
v (t 1)
Vận tốc mới của vật thể bằng vận tốc cũ của vật thể cộng gia tốc:
d i
rand .v (t) i
d i
d a (t) i
(3.23)
x (t 1)
v (t 1)
Bước 6-Cập nhật vị trí các vật thể:
d i
d x (t) i
d i
(3.24)
Bước 7-Lặp lại chu trình trên
23
Lặp lại từ bước 2 đến bước 6 đến khi gặp điều kiện dừng, thường điều kiện
dừng là hàm mục tiêu đủ tốt hoặc số lần lặp đến tối đa.
3.3 CÁC ĐẶC TÍNH CỦA GIẢI THUẬT GSA
+ Do mỗi phần tử có thể quan sát đặc tính của các phần tử khác, cho nên lực
trọng trường là công cụ truyền thông tin.
+ Do lực này tác dụng lên mỗi phần tử từ các phần tử lân cận, nên nó có thể
quan sát được không gian xung quanh nó.
+ Một vật thể nặng có bán kính tương tác hiệu quả lớn vì thế nó có cường độ
hấp dẫn mạnh. Cho nên các vật thể có hiệu quả cao hơn thì có khối lượng trọng
trường lớn hơn. Kết quả là các vật thể có khuynh hướng di chuyển về phần tử tốt
nhất.
+ Khối lượng quán tính của vật thể đặc trưng cho sự chống lại sự chuyển
động và làm cho sự chuyển động của vật thể chậm. Vì vậy, vật thể có khối lượng
quán tính lớn sẽ chuyển chậm lại giúp cho việc tìm kiếm địa phương tốt hơn.
+ Hằng số trọng trường hiệu chỉnh độ chính xác tìm kiếm, cho nên nó được
giảm theo thời gian.
+ GSA là giải thuật không cần nhớ, tức là nó cập nhật vận tốc, vị trí dựa trên
thông tin của hiện tại của vật thể. Điều này khác với PSO là cần thông tin quá khứ
của cá thể là Pbesti, Gbest. Tuy nhiên GSA lại làm việc hiệu quả giống như giải thuật
cần nhớ.
+ Trong giải thuật này, tác giả giả sử khối lượng quán tính và khối lượng
trọng trường là một. Tuy nhiên, ở một số bài toán khác nhau thì chúng có thể được
sử dụng. + Khối lượng quán tính lớn sẽ làm vật thể di chuyển chậm cho nên
tăng độ chính xác tìm kiếm, ngược lại, khối lượng trọng trường lớn sẽ hấp dẫn các
vật thể khác mạnh khiến cho tốc độ hội tụ nhanh.
3.4 ẢNH HƯỞNG CỦA CÁC THAM SỐ TRONG GSA
Ta biết rằng, một giải thuật stochastic search optimization (tối ưu hóa tìm
kiếm ngẫu nhiên) thì giai đoạn đầu phải exploration (mở rộng không gian tìm kiếm)
còn giai đoạn sau là exploitation (tìm kiếm chi tiết). Việc lựa chọn cá tham số điều
khiển của giải thuật tác động trực tiếp đến các quá trình này, cho nên cần lựa chọn
giá trị tham số điều khiển phù hợp.
24
3.4.1 Hệ số gia tốc trọng trường G0:
Hằng số này sẽ tác động trực tiếp đến vận tốc của các vật thể, do đó giai
đoạn đầu để có thể exploration thì vận tốc vật thể nhanh để có thể khám phá được
hết các khu vực trong không gian tìm kiếm, đến giai đoạn sau thì khả năng
exploitation tăng lên. Theo số lần lặp thì khả năng exploration giảm dần còn khả
năng exploitation tăng dần. Ngay từ đầu giá trị G0 phải chọn trong một vùng phù
hợp tùy theo từng bài toán. Ví dụ, khi vùng tìm kiếm và giá trị số biến nhiều thì G0
cũng phải chọn khá cao để đảm bảo exploration được hiệu quả, tuy nhiên, nếu G lớn
quá thì exploitation không hiệu quả vào giai đoạn cuối, ngược lại, G0 nhỏ quá thì
vận tốc di chuyển của vật thể quá nhỏ trong giai đoạn đầu khiến tốc độ hội tụ giảm,
và có thể kẹt ở cực trị địa phương.
3.4.2 Hệ số suy giảm α :
Hằng số trọng trường G0 quyết định việc exploration, nhưng giai đoạn này
kéo dài lâu hay mau thì phụ thuộc vào hệ số α. Khi α càng lớn thì tốc độ suy giảm G
càng lớn, điều này có nghĩa là giai đoạn exploration nhanh chóng, còn Explotation
dài hơn, và ngược lại khi α nhỏ thì tốc độ suy giảm G chậm, exploration chậm, còn
exploration nhanh. Chọn một giá trị α phù hợp sẽ quyết định đến chất lượng tìm
kiếm.
3.4.3 Số cá thể trong tập hợp N:
Hai thông số G0, α ảnh hưởng đến việc exploration, exploitation theo cách
khác với tham số N ảnh hưởng. Việc phân bố N phần tử ban đầu của việc tìm kiếm
được nhiều bài báo nghiên cứu, nó ảnh hưởng đến chất lượng tìm kiếm và nhất là
thời gian tìm kiếm. Thời gian tìm kiếm phụ thuộc vào thời gian tính toán, khi số cá
thể tăng lên thì số lượng phép tính tăng lên nhiều lần. Đối với một giải thuật thì việc
tăng số lượng cá thể lên cũng không cải thiện được nhiều đối với chất lượng lời giải
trong khi đó thời gian tính toán tăng lên nhiều.
3.5 KẾT LUẬN
Mặc dù thực tế rằng GSA vẫn là một thuật toán gần đây, sự phát triển của
các nghiên cứu liên quan GSA đã được hứa hẹn. Tại thời điểm viết bài, có những
biến thể khác nhau của GSA mà đã được phát triển và các thuật toán đã được áp
dụng trong việc giải quyết vấn đề khác nhau như trong xử lý ảnh, phân loại, phân
25
nhóm, tối ưu hóa mạng đa mục tiêu, mô hình bộ lọc, thiết kế bộ điều khiển,.... Dựa
trên đánh giá, có thể thấy rằng GSA sẽ được ứng dụng rộng rãi trong việc giải quyết
các vấn đề kỹ thuật đặc biệt là các vấn đề liên quan đến hệ thống điện và thiết kế bộ
điều khiển.
Các biến thể GSA khác nhau đã được đề xuất để khắc phục những điểm yếu
của GSA. Dựa trên các tài liệu, các GSA ban đầu có một số điểm yếu như sử dụng
các toán phức tạp và thời gian tính toán lâu dài [38]. GSA cũng bị từ tốc độ tìm
kiếm chậm trong các lần lặp cuối cùng [39]. Một vấn đề khác là khó khăn cho việc
lựa chọn thích hợp của tham số hằng số hấp dẫn, G. Tham số điều khiển tìm kiếm
chính xác và không đảm bảo một giải pháp hoàn hảo ở tất cả các mốc thời gian [40].
Mặc dù có những điểm yếu, GSA đã được chấp nhận rộng rãi do nó dễ dàng
thực hiện và khả năng giải quyết vấn đề tối ưu hóa của hệ thống kỹ thuật phi tuyến
phức tạp [41]. Các nghiên cứu đã chỉ ra rằng GSA có thể thực hiện hiệu quả về thời
gian tính toán của CPU và cho ra kết quả phù hợp hơn với độ chính xác cao hơn
[42]. GSA đã được chứng minh để làm tốt hơn tối ưu hóa khác các thuật toán như
PSO về tốc độ hội tụ và tránh được các cực tiểu địa phương [39] và có thể tạo ra
giải pháp chất lượng tốt hơn trong thời gian tính toán ngắn và đặc tuyến hội tụ ổn
định hơn so với PSO và GA [40].
Xét về phạm vi của các ứng dụng, vẫn còn nhiều lĩnh vực có thể ứng dụng
được GSA. Nó được quan sát thấy rằng ứng dụng GSA chưa thâm nhập nhất định
vào các lĩnh vực như tài chính, kinh tế, quân sự và các lĩnh vực y tế. nhiều nghiên
cứu nên được thực hiện bởi cộng đồng khoa học để kiểm tra khả năng GSA trong
các lĩnh vực. Do đó, trong những năm tới, dự kiến sẽ có nhiều chủ đề được nghiên
cứu dựa trên GSA trong nhiều lĩnh vực đa dạng hơn.
26
CHƯƠNG 4 ỨNG DỤNG THUẬT TOÁN GSA GIẢI BÀI
TOÁN ĐIỀU PHỐI KINH TẾ TRONG HỆ THỐNG ĐIỆN
4.1 GIẢI BÀI TOÁN ĐIỀU PHỐI KINH TẾ (ED) VỚI HÀM CHI PHÍ CÓ
XÉT ĐIỂM VAN CÔNG SUẤT:
Mục tiêu chính của bài toán điều phối kinh tế là để giảm thiểu chi phí vận hành
của các nhà máy điện trong khi vẫn thõa mãn các ràng buộc và đáp ứng nhu cầu của
các phụ tải trong một hệ thống điện. Trong chương này tôi xét loại bài toán điều
phối kinh tế công suất phát của các nhà máy điện có xét ảnh hưởng của điểm van
công suất.
4.1.1 Các bước áp dụng thuật toán GSA giải bài toán ED:
Bước 1: Khởi tạo vị trí các vật thể.
Khởi tạo một tập hợp các vật thể có vị trí ngẫu nhiên trong không gian bài
toán D chiều sử dụng hàm phân bố xác suất. Để làm điều này GSA xem hệ thống có
N phần tử, vị trí thứ i được định nghĩa như sau:
= ( , ,…. ,…, ) với i = 1,2,…,N. (5.1)
d là công suất phát của nhà máy thứ i theo chiều d thỏa mãn biểu d được đưa ra như dưới đây để đáp ứng hạn chế:
Ở đây Pi
d =
thức (5.13). Giá trị của Pi
d) của các nhà máy đã được
+ rand ( - ) (5.2) Pi
Bước 2: Đánh giá hàm chi phí nhiên liệu F(PGi
khởi tạo với i = 1,2,…,N.
Bước 3: Cập nhật G(t), best(t), worst(t), Mi(t) đối với i=1...N:
Trong GSA, hằng số trọng trường G, lúc đầu lấy một giá trị khởi tạo G0, sau
đó sẽ được giảm dần theo thời gian: G(t) = G0 exp(-α(t/T)) (5.3)
Ở đây, G0 là giá trị ban đầu của các hằng số hấp dẫn được lựa chọn ngẫu
nhiên, α là một hằng số, t là lần lặp hiện tại và T là tổng số lần lặp.
Đánh giá sự phù hợp của tất cả các vật thể trong không gian tìm kiếm.
cho thấy các giá trị mục tiêu của các vật thể thứ i theo thời gian t, giá trị
worst(t) và best(t) được định nghĩa như sau:
Best(t) = min( ) (5.4)
27
Worst(t) = max( ) (5.5)
Với best(t) và tồi worst(t) là giá trị mục tiêu tốt nhất và tệ nhất của các vật
thể tương ứng. Khối lượng của vật thể được tính toán sau khi tính được hàm mục
tiêu của tập hợp hiện tại:
i
(5.6) Mai=Mpi=Mii=Mi, i=1,2,…N
m (t) i
fit (t) worst(t) best(t) worst(t)
(5.7)
M (t) i
m (t) i N
m (t) j
j 1
(5.8)
Bước 4: Tính toán lực tổng hợp theo các hướng khác nhau.
Lực tổng hợp của tập các vật thể có khối lượng nặng hơn sẽ tác dụng lên vật
M (t)xM (t)
(x (t) x (t))
thể i được tính dựa vào lực trọng trường:
d F (t) i
rand G(t). j
d i
d j
i
j Kbest, j
i
j R (t) ij
(5.9)
Trong đó:
+ randj là giá trị ngẫu nhiên từ [0,1]
+ ε là một số nhỏ và Rij(t) là khoảng cách Euclidian giữa 2 vật thể i và j.
+ Kbest là tập hợp các phần tử có khả năng tác dụng lực lên các phần tử còn
lại.
Bước 5: Tính toán gia tốc và vận tốc của các vật thể.
(x (t) x (t))
Gia tốc của vật thể được tính nhờ luật chuyển động:
d a (t) i
rand G(t). j
d j
d i
j Kbest, j
i
d F (t) i M (t) i
M (t) i R (t) ij
(5.10)
Trong đó randj là số ngẫu nhiên trong đoạn [0,1]
v (t 1)
Vận tốc mới của vật thể bằng vận tốc cũ của vật thể cộng gia tốc:
d i
rand .v (t) i
d i
d a (t) i
(5.11)
x (t 1)
v (t 1)
Bước 6: Cập nhật vị trí các vật thể. Khởi tạo các ma trận X_max, X_min và X.
d i
d x (t) i
d i
(5.12)
Tính toán giá trị [best best_X] = min(fitness); [best best_X] = max(fitness).
Fbest = best; Lbest = X(best_X, :)
28
Bước 7: Lặp lại chu trình trên.
Lặp lại từ bước 2 đến bước 6 đến khi gặp điều kiện dừng là hàm giá trị chi phí
nhiên liệu của tất cả các nhà máy đủ tốt Fbest(t) ≥ best(t) thì X = Xmin, Fbest(t) <
best(t) thì X = Xmax, hoặc số lần lặp đến tối đa.
4.1.2 Lưu đồ giải thuật của thuật toán GSA:
Lưu đồ thể hiện các bước trong lập trình và tính toán mô phỏng bằng phần
mềm Matlab của giải thuật GSA, phần này tôi xin trình bày tóm tắt nhằm nêu lên
được ý chính của giải thuật.
29
Bước 1: Khởi tạo các thông số của thuật toán: f(d)= sum(c + b*X(d,:) + a.*(X(d,:).^2)) + abs(e.*sin(f.*(Xmin - X(d,:))),
với N, a , b , c , e , f , PD, PL
Bước 2: Đánh giá hàm mục tiêu của bài toán ED với các thông số đầu vào đã xác định
Bước 3: Cập nhật các thông số G(t), best(t), worst(t), Mi(t),...đối với i = 1:N
i
m (t) i
M ( t ) i
N
m ( t ) j
j 1
best(t) max fit (t), j 1..N
j
j
fit (t) worst(t) best(t) worst(t) m ( t ) i
worst(t) min fit (t), j 1..N
Bước 4: Tính toán lực tổng hợp theo các hướng khác nhau:
x (t))
d F ( t) i
rand G ( t). j
d (x (t) j
d i
j Kbest , j
i
M (t)xM ( t) j R (t) ij
i
Bước 5: Tính toán gia tốc và vận tốc:
(x (t) x (t))
d a (t) i
rand G(t). j
d i
d j
j Kbest, j
i
M (t) i R (t) ij
d F (t) i M (t) i
v (t 1)
d i
rand .v (t) i
d i
d a (t) i
Bước 6: Cập nhật vị trí các vật thể
x (t 1)
d i
d x (t) i
d v (t 1) i
Bước 7: Lặp lại từ bước 2 đến bước 6 đến khi gặp
điều kiện dừng là giá trị hàm chi phí đủ tốt hoặc
số lần lặp đến tối đa.
Hình 4.1: Lưu đồ giải thuật GSA cho bài toán ED.
30
4.1.3 Hàm mục tiêu của bài toán ED:
Hàm chi phí trong bài toán điều phối kinh tế là cực tiểu tổng chi phí hệ thống
N G
)
)
F P ( G
F P ( i G i
i
1
các nhà máy phát điện.
(5.13)
Ở bài toán ED này tôi xét hàm chi phí nhiên liệu có tính đến ảnh hưởng của
điểm van điều chỉnh công suất phát trong nhà máy, thường được mô hình bằng cách
sin
)
thêm hàm sin vào hàm chi phí bậc hai cổ điển như bên dưới.
e i
a i
2 c P i G i
P G i
F P ( i G i
b P i G i
min f P i G i
(5.14)
N G
P P P G D L i
Ràng buộc cân bằng:
1 i
(5.15)
min P G i
P G i
max P G i
Ràng buộc bất cân bằng:
(5.16) , i = 1,…, NG
Trong đó:
iGP
PD là tổng công suất yêu cầu của hệ thống, PL là tổn thất trong hệ thống,
là công suất phát của nhà máy thứ i.
Hình 4.2: Sơ đồ mạng nhất thứ các nhà máy điện phân phối công suất đến phụ tải
thông qua hệ thống truyền tải điện.
31
4.2 GIẢI BÀI TOÁN ED MẠNG 3 NÚT:
Số liệu của bài toán:
Các thông số của hệ thống này được nêu trong Bảng 4.1. Trong trường hợp
này, tổng nhu cầu phụ tải sẽ được xác định là PD = 850 (MW), tổn thất PL = 0
(MW) và các thông số khác của các máy phát điện được lấy từ [43].
Bảng 4.1: Các thông số của bài toán ED mạng 3 nút. min
max
Nút bi ai ci Pi ei Pi fi
100 100 50 600 400 200 561 310 78 7.92 7.85 7.97 0.001562 0.001940 0.004820 300 200 150 0.0315 0.042 0.063 1 2 3
Các thông số thiết lập cho thuật toán GSA được đề xuất trong bảng 4.2 như sau:
Bảng 4.2: Thông số của thuật toán GSA áp dụng giải bài toán ED mạng 3 nút.
THÔNG SỐ GIÁ TRỊ
Hằng số trọng trường (G0) Hằng số suy giảm (α) Khoảng cách giữa các vật thể (D) Số vòng lặp (Int) 100 10 10 100
Kết quả tính toán của bài toán:
Bảng 4.3: Kết quả tính toán bài toán ED mạng 3 nút.
NÚT Pmin Pmax CÔNG SUẤT
100 100 50 600 400 200 393.1697 334.6039 122.2264
1 2 3 Tổng công suất phát (MW) 850
Giá trị chi phí tốt nhất ($/h) 8194.3561
Số vòng lặp (vòng) 100
Thời gian / số vòng lặp 0.044942
32
So sánh kết quả tính toán với các thuật toán khác:
Bảng 4.4: So sánh kết quả tính toán bài toán ED mạng 3 nút.
PHƯƠNG PHÁP GIÁ TRỊ CHI PHÍ ($/h) Thời gian / số vòng lặp
8196.14 0.356 MRPSO [37]
8194.48 0.091 CS [37]
8197.16 0.3561 CPSO [37]
8194.58 0.0986 DE [44]
8194.59 0.0989 BBO [45]
8194.35 GSA 0.044942
Đồ thị kết quả tính toán:
Hình 4.3: Đồ thị giá trị hàm chi phí của bài toán ED mạng 3 nút.
33
Hình 4.4: Đồ thị phân bố công suất tại các nút của bài toán ED mạng 3 nút.
Hình 4.5: Giao diện tính toán của bài toán ED mạng 3 nút.
34
Nhận xét :
Áp dụng thuật toán GSA giải bài toán ED mạng 3 nút (bằng phần mềm
Matlab phiên bản 8.1.0.604 phát hành năm 2013 chạy trên máy tính hiệu DELL cấu
hình: Pentium(R) Dual-Core, CPU E5500, 2.80GHz, RAM 3GB) cho kết quả tốt
hơn 0.00159% và tốc độ tính toán nhanh hơn 50.61% so với thuật toán CS [37] đã
nêu ở bảng 4.4, đây cũng là tiền đề để áp dung thuật toán GSA cho các bài toán hệ
thống điện lớn hơn.
4.3 GIẢI BÀI TOÁN ED MẠNG 13 NÚT:
Các thông số đầu vào của bài toán:
Các thông số của hệ thống này được nêu trong Bảng 4.5. Trong trường hợp
này, tổng nhu cầu phụ tải sẽ được xác định là PD = 1800 (MW), tổn thất PL = 0
(MW) và các thông số khác của các máy phát điện được lấy từ [43].
Bảng 4.5: Các thông số của bài toán ED mạng 13 nút. min
max
Nút bi ai ei ci Pi Pi fi
0 680 550 8.1 0.00028 300 0.035 1
0 360 309 8.1 0.00056 200 0.042 2
0 360 307 8.1 0.00056 200 0.042 3
60 180 240 7.74 0.00324 150 0.063 4
60 180 240 7.74 0.00324 150 0.063 5
60 180 240 7.74 0.00324 150 0.063 6
60 180 240 7.74 0.00324 150 0.063 7
60 180 240 7.74 0.00324 150 0.063 8
60 180 240 7.74 0.00324 150 0.063 9
40 120 126 8.6 0.00284 100 0.084 10
40 120 126 8.6 0.00284 100 0.084 11
55 120 126 8.6 0.00284 100 0.084 12
55 120 126 8.6 0.00284 100 0.084 13
Các thông số thiết lập cho thuật toán GSA được đề xuất trong bảng 4.6 như sau:
35
Bảng 4.6: Thông số của thuật toán GSA áp dụng giải bài toán ED mạng 13 nút.
THÔNG SỐ GIÁ TRỊ
Hằng số trọng trường (G0) 100
10 Hằng số suy giảm (α)
Khoảng cách giữa các vật thể (D) 50
Số vòng lặp (Int) 1000
Kết quả tính toán của bài toán:
Bảng 4.7: Kết quả tính toán bài toán ED mạng 13 nút.
NÚT Pmin Pmax CÔNG SUẤT (MW)
1 0 680 482.3320
2 0 360 243.6327
3 0 360 248.4931
4 60 180 109.4314
5 60 180 112.4325
6 60 180 108.3449
7 60 180 104.4846
8 60 180 100.7960
9 60 180 98.2176
10 40 120 40.2990
11 40 120 40.2351
12 55 120 55.6431
13 55 120 55.6579
Tổng công suất phát (MW) 1800
Giá trị chi phí tốt nhất ($/h) 17934.8774
Số vòng lặp (vòng) 1000
Thời gian / số vòng lặp 0.067971
So sánh kết quả tính toán với các thuật toán khác:
Bảng 4.8: So sánh kết quả tính toán bài toán ED mạng 13 nút.
PHƯƠNG PHÁP GIÁ TRỊ HÀM CHI PHÍ ($/h) Thời gian / số vòng lặp
18048.21 0.0978 CEP [46]
36
18018.00 N/A FEP [46]
18028.09 0.1053 MFEP [12]
17994.07 N/A IFEP [12]
17934.87 GSA 0.0679
Đồ thị kết quả tính toán:
Hình 4.6: Đồ thị giá trị hàm chi phí của bài toán ED mạng 13 nút.
37
Hình 4.7: Đồ thị phân bố công suất tại các nút của bài toán ED mạng 13 nút.
Hình 4.8: Giao diện tính toán của bài toán ED mạng 13 nút.
Nhận xét :
38
Áp dụng thuật toán GSA giải bài toán ED mạng 13 nút (bằng phần mềm
Matlab phiên bản 8.1.0.604 phát hành năm 2013 chạy trên máy tính hiệu DELL cấu
hình: Pentium(R) Dual-Core, CPU E5500, 2.80GHz, RAM 3GB) cho kết quả tốt
hơn 0.33% so với thuật toán IFEP [13] và tốc độ tính toán nhanh hơn 44.03% so với
thuật toán CEP [46] đã nêu ở bảng 4.8, đây cũng là tiền đề để áp dung thuật toán
GSA cho các bài toán hệ thống điện lớn hơn.
4.4 GIẢI BÀI TOÁN ED 40 MẠNG NÚT:
1/Các thông số đầu vào của bài toán:
Các thông số của hệ thống này được nêu trong Bảng 4.9. Trong trường hợp
này, tổng nhu cầu phụ tải sẽ được xác định là PD = 10500 (MW), tổn thất PL = 0
(MW) và các thông số khác của các máy phát điện được lấy từ [43].
39
Bảng 4.9: Các thông số của bài toán ED mạng 40 nút.
max
max
Nút Pi
min Pi
ai
bi
ci
ei
fi
Nút Pi
min Pi
ai
bi
ci
ei
fi
36
114
94.705 6.73
0.0069
100 0.084
254
550
785.96 6.63 0.00298 300 0.035
21
1
36
114
94.705 6.73
0.0069
100 0.084
254
550
785.96 6.63 0.00298 300 0.035
22
2
60
120
309.54 7.07 0.02028 100 0.084
254
550
794.53 6.66 0.00284 300 0.035
23
3
80
190
369.03 8.18 0.00942 150 0.063
254
550
794.53 6.66 0.00284 300 0.035
24
4
47
97
148.89 5.35
0.0114
120 0.077
254
550
801.32
7.1
0.00277 300 0.035
25
5
68
140
222.33 8.05 0.01142 100 0.084
254
550
801.32
7.1
0.00277 300 0.035
26
6
110
300
287.71 8.03 0.00357 200 0.042
150
1055.1 3.33 0.52124 120 0.077
10
27
7
135
300
391.98 6.99 0.00492 200 0.042
150
1055.1 3.33 0.52124 120 0.077
10
28
8
135
300
455.76
6.6
0.00573 200 0.042
150
1055.1 3.33 0.52124 120 0.077
10
29
9
130
300
722.82 12.9 0.00605 200 0.042
97
148.89 5.35
0.0114
120 0.077
47
30
10
94
375
635.2
12.9 0.00515 200 0.042
190
222.92 6.43
0.0016
150 0.063
60
31
11
94
375
654.69 12.8 0.00569 200 0.042
190
222.92 6.43
0.0016
150 0.063
60
32
12
125
500
913.4
12.5 0.00421 300 0.035
190
222.92 6.43
0.0016
150 0.063
60
33
13
125
500
1760.4 8.84 0.00752 300 0.035
200
107.87 8.95
0.0001
200 0.042
90
34
14
125
500
1728.3 9.15 0.00708 300 0.035
200
116.58 8.62
0.0001
200 0.042
90
35
15
125
500
1728.3 9.15 0.00708 300 0.035
200
116.58 8.62
0.0001
200 0.042
90
36
16
220
500
647.85 7.97 0.00313 300 0.035
110
307.45 5.88
0.0161
80
0.098
25
37
17
220
500
649.69 7.95 0.00313 300 0.035
110
307.45 5.88
0.0161
80
0.098
25
38
18
242
550
647.83 7.97 0.00313 300 0.035
110
307.45 5.88
0.0161
80
0.098
25
39
19
242
550
647.81 7.97 0.00313 300 0.035
242
550
647.83 7.97 0.00313 300 0.035
40
20
Các thông số thiết lập cho thuật toán GSA được đề xuất trong bảng 4.10 như sau:
Bảng 4.10: Thông số của thuật toán GSA áp dụng giải bài toán ED mạng 40 nút.
THÔNG SỐ GIÁ TRỊ
Hằng số trọng trường (G0) 100
10 Hằng số suy giảm (α)
Khoảng cách giữa các vật thể (D) 50
Số vòng lặp (Int) 2000
Kết quả tính toán của bài toán:
Bảng 4.11: Kết quả tính toán bài toán ED mạng 40 nút.
CÔNG SUẤT STT Pmin Pmax (MW)
36 114 113.9969 1
36 114 113.9999 2
40
60 120 119.9999 3
80 190 190.0000 4
47 97 97.0000 5
68 140 140.0000 6
110 300 300.0000 7
135 300 300.0000 8
135 300 299.9665 9
130 300 246.3067 10
94 375 235.9074 11
94 375 214.0734 12
125 500 280.4458 13
125 500 275.1626 14
125 500 294.3407 15
125 500 313.8960 16
220 500 363.6258 17
220 500 406.0826 18
242 550 436.9626 19
242 550 473.3370 20
254 550 494.2423 21
254 550 492.1875 22
254 550 510.4019 23
254 550 539.2708 24
254 550 543.0784 25
254 550 544.5812 26
10 150 10.2811 27
10 150 10.4759 28
10 150 10.8424 29
47 97 96.6642 30
60 190 189.9986 31
60 190 189.8529 32
41
60 190 189.9947 33
90 200 198.8544 34
90 200 195.7301 35
90 200 199.0430 36
25 110 108.1421 37
25 110 109.4256 38
25 110 108.1535 39
242 550 543.6758 40
Tổng công suất phát 10500 (MW)
Giá trị chi phí tốt nhất 121044.4873 ($/h)
Số vòng lặp (vòng) 2000
Thời gian / số vòng lặp 0.076212
So sánh kết quả tính toán với các thuật toán khác:
Bảng 4.12: So sánh kết quả tính toán bài toán ED mạng 40 nút.
PHƯƠNG PHÁP GIÁ TRỊ HÀM CHI PHÍ ($/h) Thời gian / số vòng lặp
121416.29 0.615 DE [47]
121663.52 N/A APSO [48]
121468.82 0.318 PSO-RDL [49]
121430.00 N/A SA-PSO [50]
121044.48 0.0762 GSA
Đồ thị kết quả tính toán:
42
Hình 4.9: Đồ thị giá trị hàm chi phí của bài toán ED mạng 40 nút.
Hình 4.10: Đồ thị phân bố công suất tại các nút của bài toán ED mạng 40 nút.
43
Hình 4.11: Giao diện tính toán của bài toán ED mạng 40 nút.
Nhận xét:
Áp dụng thuật toán GSA giải bài toán ED mạng 40 nút (bằng phần mềm
Matlab phiên bản 8.1.0.604 phát hành năm 2013 chạy trên máy tính hiệu DELL cấu
hình: Pentium(R) Dual-Core, CPU E5500, 2.80GHz, RAM 3GB) cho kết quả tốt
hơn 0.306% so với thuật toán DE [47] và tốc độ tính toán nhanh hơn 76.03% so
thuật toán PSO-RDL [49] đã nêu ở bảng 4.12, đây cũng là tiền đề để áp dung thuật
toán GSA cho các bài toán hệ thống điện lớn hơn và có phức tạp hơn.
4.5 KẾT LUẬN:
Phương pháp tiếp cận sử dụng giải thuật GSA có thể giải quyết tốt bài toán
ED có xét điểm van công suất. Kết quả tính toán từ các phương pháp trên có thể so
sánh với các kết quả tốt nhất được tìm thấy đến thời điểm hiện tại của bài toán ED
với mạng 3 nút, 13 nút và 40 nút.
44
Đối với kết quả tính toán cho mạng 3 nút, 13 nút và 40 nút, giải thuật GSA
cho kết quả tốt hơn và thời gian hội tụ nhanh so với các giải thuật khác, điều này
cho thấy rằng việc kết hợp giữa thuật toán GSA có thể áp dụng đc cho các mạng có
số lượng nút lớn hơn nữa. Đối với luận văn này, vẫn còn khó khăn trong quá trình
lập trình mô phỏng các hệ thống mạng điện do có sự hạn chế kinh nghiệm lựa chọn
thông số cài đặt tính toán cho thuật toán của người lập trình.
Trong khuôn khổ luận văn này, thuật toán GSA sẽ là hướng đi mới để tôi
hướng đến việc nghiên cứu áp dụng tính toán thực tế vào các mạng điện sau này
dựa trên nền tảng của bài toán kinh điển ED đã được trình bày.
45
CHƯƠNG 5 KẾT LUẬN VÀ HƯỚNG PHÁT TRIỂN
5.1 TỔNG KẾT
Từ những kết quả nghiên cứu ta rút ra được những kết luận sau:
Ưu điểm của Thuật toán GSA:
- Thuật toán GSA tìm kiếm trong tất cả các không gian bài toán chứ
không phải riêng từng điểm.
- Tìm được điểm tối ưu toàn cục, bởi vì việc tính toán dựa trên các cá thể
riêng biệt, GSA thừa hưởng khả năng tính toán đồng thời.
- GSA sử dụng các hàm mục tiêu và hàm tính toán độ phù hợp để trả về
trực tiếp kết quả. GSA thích hợp với những hàm mục tiêu không liên tục,
không khả vi tồn tại phổ biến trong hệ thống điện.
- GSA có khả năng tìm kiếm trong những vùng không gian phức tạp,
không chắc chắn để tìm ra lời giải tối ưu toàn cục do đó nó linh hoạt và
tốt hơn các phương pháp truyền thống khác.
- Thuật toán có tính chất mở, thời gian phát triển không quá dài nên việc
cải tiến thuật toán cũng như ứng dụng thuật toán vào các lĩnh vực mới
thực sự hấp dẫn và đem lại nhiều kỳ vọng.
- Trong nội dung luận văn đã áp dụng thuật toán GSA vào giải bài toán
ED, từ những kết quả thu được có thể khẳng định GSA giải quyết tốt bài
toán này, thực sự có thể so sánh với các thuật toán thông minh khác.
Nhược điểm của GSA:
- Thuật toán GSA sử dụng các toán tử phức tạp và thời gian tính toán lâu
dài.
- Kết quả tính toán phụ thuộc nhiều vào các thông số cài đặt ban đầu, buộc
người giải phải thử qua nhiều phương án khác nhau, vận dụng yếu tố
kinh nghiệm. Đây là điểm yếu cố hữu của các thuật toán từ trước đến
nay, gần đây các kỹ thuật đáp ứng thích nghi ra đời đã hứa hẹn giải
quyết được vấn đề này.
46
5.2 HƯỚNG PHÁT TRIỂN ĐỀ TÀI
Từ những phân tích trên tôi xin đưa ra hướng phát triển tiếp theo của đề tài:
- Giải bài toán phân bố công suất tối ưu trong mạng điện có số lượng nút
lớn hơn.
- Tiếp tục nghiên cứu khai thác phương án GSA để giải các bài toán khác,
trên cơ sở đó cải thiện thêm cho phương pháp này.
- Nghiên cứu các thuật toán khác để đưa đến những ý tưởng mới cùng với
kinh nghiệm từ luận văn này để kết hợp tạo ra thuật toán lai giữa thuật
toán GSA với các thuật toán phù hợp nhăm đạt được kết quả tính toán
tốt hơn.
- Tiến tới áp dụng thuật toán tối ưu cho các bài toán khác như: bài toán
quy hoạch và phát triển hệ thống điện, bài toán phân bố công suất tối ưu
trong thị trường điện cạnh tranh, bài toán tối ưu phân bố công suất trong
hệ thống Thuỷ điện – Nhiệt điện...
- Tìm hiểu bài toán phân bố công suất tối ưu OPF, áp dụng thuật toán
GSA, và các phương án cải tiến để giải quyết bài toán này.
Thời gian thực hiện luận văn này không phải là quá ít cho một đề tài nghiên
5.3 LỜI KẾT
cứu nhưng cũng không đủ dài để có thể hoàn chỉnh một công trình khoa học nhất là
trong điều kiện vẫn còn tồn tại những khó khăn nhất định về kiến thức, kinh nghiệm
cũng như việc tìm kiếm tài liệu tham khảo chính thống và sử dụng sao cho đạt hiệu
quả cao nhất.
Luận văn tốt nghiệp có thể xem là công trình khoa học đầu tiên và quan
trọng đối với tôi. Trong quá trình tìm hiểu và nghiên cứu thuật toán tối ưu GSA, tôi
đã tìm thấy được sự thích thú, đam mê. Tôi thực sự cố gắng để hoàn thành tốt luân
văn này.
Tuy nhiên cũng không thể tránh được những sai sót, nhầm lẫn rất mong sự
chia sẻ đóng góp của hội đồng cũng như của người đọc để tôi bổ sung hoàn chỉnh
trong tương lai.
Xin chân thành cảm ơn quí thầy (cô) !
47
TÀI LIỆU THAM KHẢO
[1] Nguyễn Hồng Việt Phương, Ngô Văn Dưỡng, “A calculation for building a
simulation model to research the operation systems of UPFC,” Đại học Đà Nẵng.
[2] Scott B, Marinho JL, “Linear programming for power system network security
applications,” IEEE Trans.on PAS, vol 98, no 3, pp 837-848, 1979.
[3] Alsc O, Bright J, Praise M, Scott B, “Further developments in LP – basesd
optimal power flow,” IEEE Transactions on Power System, vol5, no3, pp 697-711,
1990.
[4] C.R. Fuerte-Esquivel and E. Acha. “A Newton-type Algorithm for the control of
power flow in electrical power network,” IEEE Transactions on Power Systems,
Vol. 12, No.4, November 1997. pp.1471-1480.
[5] Phạm Việt Cường, “Ứng dụng thuật toán di truyền phân bố công suất tối ưu
trong hệ thống điện,” Đại học Bách Khoa Tp Hồ Chí Minh, Luận văn Thạc sĩ,
07/2003, 700998, Thư viện Đại Học Bách Khoa Tp.HCM.
[6] Tarek Bouktir, Linda Slimani, M. Belkacemi, “A Genetic Algorithm for Solving
the Optimal Power Flow Problem,” Leonardo Journal of Sciences, Issue 4, January-
June 2004, p. 44-58.
[7] K.Vaisakh, L.R.Srinivas, “Differential Evolution Approach for Optimal Power
Flow Solutions,” Journal of Theoretical and Applied Information Technology, 2005
- 2008 JATIT, p. 261-268.
[8] Boumediène ALLAOUA, Abdellah LAOUFI, “Optimal Power Flow Solution
Using Ant Manners for Electrical Network,” Advances in Electrical and Computer
Engineering, Volume 9, Number 1, 2009, p. 34-40.
[9] Ding Xiaoying, Wang Xifan and Liu Lin, “Interior Point Cutting Plane Method
for Discrete Decoupled Optimal Power Flow,” IAENG International Journal of
Applied Mathematics, 37:1, IJAM 37_1_8.
[10] E. Rashedi, H. Nezamabadi-pour, and S. Saryazdi, “GSA: A Gravitational
Search Algorithm,” Information Sciences, vol. 179, no. 13, (2009), pp. 2232–2248.
[11] Norlina Mohd Sabri, Mazidah Puteh, and Mohamad Rusop Mahmood “A
Review of Gravitational Search Algorithm” Int. J. Advance. Soft Comput. Appl.,
Vol. 5, No. 3, November 2013 ISSN 2074-8523.
48
[12] P.K.Swain, N.C.Sahu, P.K.Hota “Gravitational Search Algorithm for Optimal Economics Dispatch” 2nd International Conference on Communication, Computing
& Security [ICCCS – 2012], Procedia Technology 6 (2012), pp 441-419.
[13] An Improved Particle Swarm Optimization for Nonconvex Economic Dispatch
Problems Jong-Bae Park, Member, IEEE, Yun-Won Jeong, Joong-Rin Shin, Senior
Member, IEEE, and Kwang Y. Lee, Fellow, IEEE.
[14] Hoàng Kiếm, Thuật giải di truyền, NXB Giáo dục, Hà nội, 2000.
[15] Bùi Minh Trí, Bài tập tối ưu hoá, NXB Khoa học Kỹ thuật, Hà nội, 2002.
[16] J. S. Alsumait, J. K. Sykulski, and A. K. Al-Othman, “A hybrid GA– PS–SQP
method to solve power system valve-point economic dispatch problems,” Applied
Energy, vol. 87, no. 5, pp. 1773-1781, May 2010.
[17] A. Wood, B. Wollenberg, “Power generation, operation and control,” New
York: Wiley, 1996.
[18] D.C Walters, G. B. Sheblé, “Genetic Algorithm Solution of Economic
Dispatch with Valve Point Loading,” IEEE Trans. Power Systems, Vol. 8, No. 3,
pp. 1325-1332, August 1993.
[19] K. Wong, Y. Wong, “Genetic and genetic/simulated-annealing approaches to
economic dispatch,” IEEE Proceedings Gener, Trans and Distr, Vol. 141, No. 5,
pp. 507-513, Sep 1994.
[20] K. Wong, B. Lau, A. Fry, “Modelling Generator Input-Output Characteristics with Valve-Point Loading Using Neural Networks,” IEEE 2nd International
Conference on Advances in Power System Control Operation and Management, pp.
843-848, 7-10 Dec 1993.
[21] H. Yang, P. Yang, C Huang, “Evolutionary Programming Based Economic
Dispatch for Units with Non-Smooth Fuel Cost Functions,” IEEE Trans. Power
System, Vol. 11, No. 1, pp. 112-118, February 1996.
[22] Storn R, Price K. Differential evolution – a simple and efficient heuristic for
global optimization over continuous spaces. J Global Optim 1997;11:341–59.
[23] Price K, Storn RM, Lampinen JA. Differential evolution: a practical approach
to global optimization. New York: Springer; 2005.
49
[24] K. S. Swarup, “Swarm intelligence approach to the solution of optimal power
flow”, Indian Institute of Science, Oct. 2006, 86, 439–455.
[25] Hirotaka Yoshida, Kenichi Kawata, Yoshikazu Fukuyama, Yosuke Nakanishi,
“A Particle Swarm Optimization for Reactive Power and Voltage control
considering Voltage stability,” IEEE International Conference on Intelligent System
Applications to Power Systems (ISAP'99), Rio de Janeiro, April 4-8, 1999.
[26] An Improved Particle Swarm Optimization for Nonconvex Economic Dispatch
Problems Jong-Bae Park, Member, IEEE, Yun-Won Jeong, Joong-Rin Shin, Senior
Member, IEEE, and Kwang Y. Lee, Fellow, IEEE.
[27] D. Karaboga, B. Basturk, A powerful and efficient algorithm for numerical
function optimization: artificial bee colony algorithm, J. Global Optim. 3(2007)
9459–9471.
[28] D. Karaboga, B. Basturk, On the performance of artificial bee colony
algorithm, Appl. Soft Comput. 8 (2008) 687–697.
[29] D. Karaboga, Akay Bahriye, A comparative study of artificial bee colony
algorithm, Appl. Math. Comput. 214 (2009) 108–132].
[30] T. D. King, M. E. El-Hawary, and F. El-Hawary, “Optimal environmental
dispatching of electric power system via an improved hopfield neural network
model”, IEEE Transactions on Power Systems, Vol. 1, No. 3, pp. 1559-1565,
August 1995.
[31] K. Y. Lee, A. Sode-Yome, and J. H. Park, ‘Adaptive Hopfield Neural
Networks for Economic Load Dispatch’, IEEE Trans. on Power Systems, Vol 13
No 2 (May 1998) pp 519-526.
[32] Lee, S.C., Kim, Y.H., “An Enhanced Lagrangian Neural Network for the ELD
Problems with Piecewise Quadratic Cost Functions and Nonlinear Constraints”,
Elect. Power Systems Research, vol. 60, pp.167-177, 2002.
[33] A new meta-heuristic algorithm for continuous engineering optimization:
harmony search theory and practice Kang Seok Lee, Zong Woo Geem, Comput.
Methods Appl. Mech. Engrg. 194 (2005) 3902–3933.
[34] Xin-She Yang, Suash Deb, “Cuckoo Search via L´evy Flights”, IEEE , 2009.
50
[35] Ramin Rajabioun, “Cuckoo Optimization Algorithm,” Elsevier, Applied
Soft Computing 11 (2011) 5508–5518, 2011.
[36] Zahra Moravej, Amir Akhlaghi, “A novel approach based on cuckoo search for
DG allocation in distribution network”, Elsevier, Electrical Power and Energy
Systems 44 (2013) 672–679.
[37] Ramandeep Kaur, Gurpreet Kaur Gill, "Solution to Economic Load Dispatch
Problem Using Cuckoo Search Algorithm", International Journal of Electrical and
Electronics Research, April - June 2015, Vol. 3, Issue 2, pp: (362-369).
[38] Z. Beheshti, S. M. Shamsuddin, “A Review of Population-based MetaHeuristic
Algorithm,” International Journal of Advance Soft Computing, vol. 5, no. 1, (2013),
pp. 1–35.
[39] S. Mirjalili, S. Z. Mohd Hashim, and H. Moradian Sardroudi, “Training
feedforward neural networks using hybrid particle swarm optimization and
gravitational search algorithm,” Applied Mathematics and Computation, vol. 218,
no. 22, (2012), pp. 11125–11137.
[40] J. Vijaya Kumar, D. M. Vinod Kumar, and K. Edukondalu, “Strategic bidding
using fuzzy adaptive gravitational search algorithm in a pool based electricity
market,” Applied Soft Computing, vol. 13, (2012), pp. 2445- 2455.
[41] R. K. Khadanga and S. Panda, “Gravitational search algorithm for Unified
Power Flow Controller based damping controller design,” 2011 International
Conference on Energy, Automation and Signal, (2011), pp. 1–6.
[42] E. Rashedi, H. Nezamabadi-pour, and S. Saryazdi, “Filter modeling using
gravitational search algorithm,” Engineering Applications of Artificial Intelligence,
vol. 24, no. 1, (2011), pp. 117–122.
[43] T. K. A. Rahman, Z. M. Yasin, W. N. W. Abdullah, "Artificial-immune-based
for solving economic dispatch in power system", National Power & Energy
Conference (PECon), Kuala Lumpur, Malaysia, 2004, pp. 31-35.
[44] Abido M. A, 2003 “Environmental/economic power dispatch using
multiobjective evolutionary algorithms” IEEE Trans. Power System, Vol. 18,No. 4,
pp. 1529–1537.
51
[45] Zhang Y., Gong D.W., Ding Z., 2012, “A bare-bones multi-objective particle
swarm optimization algorithm for environmental/economic dispatch” Information
Sciences, Vol.192, pp.213–22.
[46]. Sinha N, Chakrabarti R, Chattopadhyay PK. Evolutionary programming
techniques for economic load dispatch. IEEE Trans Evolut Comput 2003;7(1):83
94.
[47] N. Noman and H. Iba, “Differential evolution for economic load dis-patch
problems,” Electric Power Systems Research, vol. 78, no. 8, pp. 1322-1331, Aug.
2008.
[48] A. I. Selvakumar and K. Thanushkodi, “Anti-predatory particle swarm
optimization: Solution to nonconvex economic dispatch problems,” Electric Power
Systems Research, vol. 78, no. 1, pp. 2-10, Jan. 2008.
[49] Y.-P. Chen, W.-C. Peng, and M.-C. Jian, “Particle swarm optimization with
recombination and dynamic linkage discovery,” IEEE Trans. Power Systems, Man,
and Cybernetics - Part B: Cybernetics, vol. 37, no. 6, pp. 1460-1470, Dec. 2007.
[50] C.-C. Kuo, “A novel coding scheme for practical economic dispatch by
modified particle swarm approach,” IEEE Trans. Power Systems, vol. 23, no. 4, pp.
1825-1835, Nov. 2008.
Phụ Lục
CODE THUẬT TOÁN GRAVITATIONAL
SEARCH ALGORITHM
Code File Main:
% GSA code v1.1.
% Generated by Esmat Rashedi, 2010.
% " E. Rashedi, H. Nezamabadi-pour and S. Saryazdi,
% “GSA: A Gravitational Search Algorithm”, Information sciences, vol. 179,
% no. 13, pp. 2232-2248, 2009."
% Main function for using GSA algorithm.
clear all;clc
% inputs:
% N: Number of agents.
% max_it: Maximum number of iterations (T).
% ElitistCheck: If ElitistCheck=1, algorithm runs with eq.21 and if =0, runs
with eq.9.
% Rpower: power of 'R' in eq.7.
% F_index: The index of the test function. See tables 1,2,3 of the mentioned
article.Insert your own objective function with a new F_index in
'test_functions.m' and 'test _functions _ range.m'.
% outputs:
% Fbest: Best result.
% Lbest: Best solution. The location of Fbest in search space.
% BestChart: The best so far Chart over iterations.
% MeanChart: The average fitnesses Chart over iterations.
fig
Code File Initialization:
% GSA code v1.1.
% Generated by Esmat Rashedi, 2010.
% " E. Rashedi, H. Nezamabadi-pour and S. Saryazdi,
%“GSA: A Gravitational Search Algorithm”, Information sciences, vol. 179,
%no. 13, pp. 2232-2248, 2009."
%This function initializes the position of the agents in the search space,
randomly. function [X]=initialization(sumfixed,Xmax,Xmin,N,dim)
% if size(Xmax,2)==1
% X=rand(N,dim).*(Xmax-Xmin)+Xmin;
% end
% if size(Xmax,2)>1
% for i=1:dim
% high=Xmax(i);low=Xmin(i);
% X(:,i)=rand(N,1).*(high-low)+low;
% end
% end
X = zeros(N,dim);
% for d = 1:1:N
% for j = 1:1:dim
% while(1)
% j
% if j==dim
% X(d,j) = sumfixed - sum(X(d,1:dim-1));
% break;
% end
% X(d,j) = Xmin(j) + rand * (Xmax(j) - Xmin(j));
% sumd = sum(Xmin(j+1:dim));
% sumu = sum(Xmax(j+1:dim));
% sumt = sumfixed - sum(X(d,1:j));
% if sumt >= sumd && sumt <= sumu
% break;
% end
% end
% end
% end
for d = 1:N
total = sumfixed;
xmin = Xmin;
xmax = Xmax;
for j = 1:dim
if(j == dim)
X(d,j) = total;
continue;
end
l = max(xmin(1),total - sum(xmax(2:end)));
p = min(xmax(1), total - sum(xmin(2:end)));
X(d,j) = l + rand * (p - l);
total = total - X(d,j);
xmin = xmin(2:end);
xmax = xmax(2:end);
end
end
Code File G_Constant:
% This function calculates Gravitational constant. eq.13.
function G=Gconstant(iteration,max_it,alfa,G0)
%%%here, make your own function of 'G'
G=G0*exp(-alfa*iteration/max_it); %eq. 28.
Code File MassCaculation:
%This function calculates the mass of each agent. eq.14-20
function [M]=massCalculation(fit,min_flag);
%%%%here, make your own function of 'mass calculation'
Fmax=max(fit); Fmin=min(fit); Fmean=mean(fit);
[i N]=size(fit);
if Fmax==Fmin
M=ones(N,1);
else
if min_flag==1 %for minimization
best=Fmin;worst=Fmax; %eq.17-18.
else %for maximization
best=Fmax;worst=Fmin; %eq.19-20.
end
M=(fit-worst)./(best-worst); %eq.15,
end
M=M./sum(M); %eq. 16.
Code File Space_bound:
%This function checks the search space boundaries for agents.
function X=space_bound(sumfixed,Xmax,Xmin,X);
N = size(X,1);
dim = size(X,2);
%[N,dim]=size(X);
% phao modified
%for i=1:N
% %%Agents that go out of the search space, are reinitialized randomly .
%
Tp=X(i,:)>Xmax;Tm=X(i,:) Xmax-low)+low).*(Tp+Tm)); % %%Agents that go out of the search space, are returned to the boundaries. % Tp=X(i,:)>Xmax;Tm=X(i,:) *Tm; %end % for d = 1:1:N % err = (X(d,:) > Xmax) + (X(d,:) < Xmin) % if sum(err)==0 && sum(X(d,:))==sumfixed % continue; % end % while(1) % sumd = sum(Xmin.*err); % sumu = sum(Xmax.*err); % sumcheck = sumfixed - sum(X(d,:).*(~err)); % check = (sumd <= sumcheck) + (sumcheck <= sumu); % % if check == 2 % break; % else % f0 = find(err==0); % err(f0(1)) = 1; % end % end % while(1) % sumrem = sumfixed - sum(X(d,:).*(~err)); % if nnz(err)==1 % f = find(err==1); % X(d,f(1)) = sumrem; % break; % else % f1 = find(err==1); % err(f1(1)) = 0; % while(1) % X(d,f1(1)) = Xmin(f1(1)) + rand * (Xmax(f1(1)) - Xmin(f1(1))); % sumd = sum(Xmin.*err); % sumu = sum(Xmax.*err); % sumt = sumrem - X(d,f1(1)); % if sumt >= sumd && sumt <= sumu % break; % end % end % end % end % end for d = 1:N err = (X(d,:) > Xmax) + (X(d,:) < Xmin); if sum(err)==0 && sum(X(d,:))==sumfixed continue; end while(1) sumd = sum(Xmin.*err); sumu = sum(Xmax.*err); sumcheck = sumfixed - sum(X(d,:).*(~err)); check = (sumd <= sumcheck) + (sumcheck <= sumu); if check == 2 break; else f0 = find(err==0); err(f0(1)) = 1; end end pos = find(err); total = sumfixed - sum(X(d,:).*(~err)); xmin = Xmin(pos); xmax = Xmax(pos); for j = 1:size(pos,2) if(j == size(pos,2)) X(d,pos(j)) = total; continue; end l = max(xmin(1),total - sum(xmax(2:end))); p = min(xmax(1), total - sum(xmin(2:end))); X(d,pos(j)) = l + rand * (p - l); total = total - X(d,pos(j)); xmin = xmin(2:end); xmax = xmax(2:end); end end Code File Test_functions: function fit=test_functions(L,F_index,dim) %Insert your own objective function with a new F_index. if F_index==1 fit=sum(L.^2); end if F_index==2 fit=sum(abs(L))+prod(abs(L)); end if F_index==3 fit=0; for i=1:dim fit=fit+sum(L(1:i))^2; end end if F_index==4 fit=max(abs(L)); end if F_index==5 fit=sum(100*(L(2:dim)-(L(1:dim-1).^2)).^2+(L(1:dim-1)-1).^2); end if F_index==6 fit=sum(floor((L+.5)).^2); end if F_index==7 fit=sum([1:dim].*(L.^4))+rand; end if F_index==8 fit=sum(-L.*sin(sqrt(abs(L)))); end if F_index==9 fit=sum(L.^2-10*cos(2*pi.*L))+10*dim; end if F_index==10 fit=-20*exp(-.2*sqrt(sum(L.^2)/dim))- exp(sum(cos(2*pi.*L))/dim)+20+exp(1); end if F_index==11 fit=sum(L.^2)/4000-prod(cos(L./sqrt([1:dim])))+1; end if F_index==12 fit=(pi/dim)*(10*((sin(pi*(1+(L(1)+1)/4)))^2)+sum((((L(1:dim- 1)+1)./4).^2).*... (1+10.*((sin(pi.*(1+(L(2:dim)+1)./4)))).^2))+((L(dim)+1)/4)^2)+sum(Ufun(L,1 0,100,4); end if F_index==13 fit=.1*((sin(3*pi*L(1)))^2+sum((L(1:dim-1)- 1).^2.*(1+(sin(3.*pi.*L(2:dim))).^2))+... ((L(dim)-1)^2)*(1+(sin(2*pi*L(dim)))^2))+sum(Ufun(L,5,100,4)); end if F_index==14 aS=[-32 -16 0 16 32 -32 -16 0 16 32 -32 -16 0 16 32 -32 -16 0 16 32 -32 -16 0 16 32;,... -32 -32 -32 -32 -32 -16 -16 -16 -16 -16 0 0 0 0 0 16 16 16 16 16 32 32 32 32 32]; for j=1:25 bS(j)=sum((L'-aS(:,j)).^6); end fit=(1/500+sum(1./([1:25]+bS))).^(-1); end if F_index==15 aK=[.1957 .1947 .1735 .16 .0844 .0627 .0456 .0342 .0323 .0235 .0246]; bK=[.25 .5 1 2 4 6 8 10 12 14 16];bK=1./bK; fit=sum((aK-((L(1).*(bK.^2+L(2).*bK))./(bK.^2+L(3).*bK+L(4)))).^2); end if F_index==16 fit=4*(L(1)^2)-2.1*(L(1)^4)+(L(1)^6)/3+L(1)*L(2)-4*(L(2)^2)+4*(L(2)^4); end if F_index==17 fit=(L(2)-(L(1)^2)*5.1/(4*(pi^2))+5/pi*L(1)-6)^2+10*(1- 1/(8*pi))*cos(L(1))+10; end if F_index==18 fit=(1+(L(1)+L(2)+1)^2*(19-14*L(1)+3*(L(1)^2)- 14*L(2)+6*L(1)*L(2)+3*L(2)^2))*... (30+(2*L(1)-3*L(2))^2*(18-32*L(1)+12*(L(1)^2)+48*L(2)- 36*L(1)*L(2)+27*(L(2)^2))); end if F_index==19 aH=[3 10 30;.1 10 35;3 10 30;.1 10 35];cH=[1 1.2 3 3.2]; pH=[.3689 .117 .2673;.4699 .4387 .747;.1091 .8732 .5547;.03815 .5743 .8828]; fit=0; for i=1:4 fit=fit-cH(i)*exp(-(sum(aH(i,:).*((L-pH(i,:)).^2)))); end end if F_index==20 aH=[10 3 17 3.5 1.7 8;.05 10 17 .1 8 14;3 3.5 1.7 10 17 8;17 8 .05 10 .1 14]; cH=[1 1.2 3 3.2]; pH=[.1312 .1696 .5569 .0124 .8283 .5886;.2329 .4135 .8307 .3736 .1004 .9991;... .2348 .1415 .3522 .2883 .3047 .6650;.4047 .8828 .8732 .5743 .1091 .0381]; fit=0; for i=1:4 fit=fit-cH(i)*exp(-(sum(aH(i,:).*((L-pH(i,:)).^2)))); end end aSH=[4 4 4 4;1 1 1 1;8 8 8 8;6 6 6 6;3 7 3 7;2 9 2 9;5 5 3 3;8 1 8 1;6 2 6 2;7 3.6 7 3.6]; cSH=[.1 .2 .2 .4 .4 .6 .3 .7 .5 .5]; if F_index==21 fit=0; for i=1:5 fit=fit-((L-aSH(i,:))*(L-aSH(i,:))'+cSH(i))^(-1); end end if F_index==22 fit=0; for i=1:7 fit=fit-((L-aSH(i,:))*(L-aSH(i,:))'+cSH(i))^(-1); end end if F_index==23 fit=0; for i=1:10 fit=fit-((L-aSH(i,:))*(L-aSH(i,:))'+cSH(i))^(-1); end end function y=Ufun(x,a,k,m) y=k.*((x-a).^m).*(x>a)+k.*((-x-a).^m).*(x<(-a)); return Code File Gfiled: %This function calculates the accelaration of each agent in gravitational field. eq.7-10,21. function a=Gfield(M,X,G,Rnorm,Rpower,ElitistCheck,iteration,max_it); [N,dim]=size(X); final_per=2; %In the last iteration, only 2 percent of agents apply force to the others. %%%%total force calculation if ElitistCheck==1 kbest=final_per+(1-iteration/max_it)*(100-final_per); %kbest in eq. 21. kbest=round(N*kbest/100); else kbest=N; %eq.9. end [Ms ds]=sort(M,'descend'); for i=1:N E(i,:)=zeros(1,dim); for ii=1:kbest j=ds(ii); if j~=i R=norm(X(i,:)-X(j,:),Rnorm); %Euclidian distanse. for k=1:dim E(i,k)=E(i,k)+rand*(M(j))*((X(j,k)-X(i,k))/(R^Rpower+eps)); %note that Mp(i)/Mi(i)=1 end end end end %%acceleration a=E.*G; %note that Mp(i)/Mi(i)=1 Code File EvaluateF: %This function Evaluates the agents. function fitness=evaluateF(X, a, b, c, e, f, Xmin); % phao modified %for i=1:N %L is the location of agent number 'i' %L=X(i,:); %calculation of objective function for agent number 'i' %fitness(i)=test_functions(L,F_index,dim); %end N = size(X,1); for d = 1:1:N % fitness(d) = sum(c + b.*X(d,:) + a.*(X(d,:).^2)+ abs(e.*sin(f.*(Xmin - X(d,:))))); fitness(d) = sum(c + b.*X(d,:) + a.*(X(d,:).^2)) ;%+abs(e.*sin(f.*(Xmin - X(d,:)))); end Code File GSA.m: % Gravitational Search Algorithm. function [Fbest,Lbest,BestChart,MeanChart]=GSA(options, AxesPower, AxesBestCost) %V: Velocity. %a: Acceleration. %M: Mass. Ma=Mp=Mi=M; %dim: Dimension of the test function. %N: Number of agents. %X: Position of agents. dim-by-N matrix. %R: Distance between agents in search space. %[low-up]: Allowable range for search space. %Rnorm: Norm in eq.8. %Rpower: Power of R in eq.7. %parse data Rnorm = options.Vars.Rnorm; Rpower = options.Vars.Rpower; min_flag = options.Vars.MinFlag; ElitistCheck = options.Vars.ElitistCheck; N = options.Vars.NbOfAgents;% dimension alfa = options.Vars.Alpha; G0 = options.Vars.G0; PD = options.Vars.PD; dim = options.Vars.Dim; max_it = options.Vars.Iterations; Xmin = options.Parameters.Pmin; Xmax = options.Parameters.Pmax; ci = options.Parameters.c; bi = options.Parameters.b; ai = options.Parameters.a; ei = options.Parameters.e; fi = options.Parameters.f; %get allowable range and dimension of the test function. % phao modified %[low,up,dim]=test_functions_range(F_index); % ndata = [ 100 600 561 7.92 0.001562 300 0.0315 % 100 400 310 7.85 0.001940 200 0.042 % 50 200 78 7.97 0.004820 150 0.063] % ndata = [0 680 550 8.10 0.00028 300 0.035 % 0 360 309 8.10 0.00056 200 0.042 % 0 360 307 8.10 0.00056 200 0.042 % 60 180 240 7.74 0.00324 150 0.063 % 60 180 240 7.74 0.00324 150 0.063 % 60 180 240 7.74 0.00324 150 0.063 % 60 180 240 7.74 0.00324 150 0.063 % 60 180 240 7.74 0.00324 150 0.063 % 60 180 240 7.74 0.00324 150 0.063 % 40 120 126 8.6 0.00284 100 0.084 % 40 120 126 8.6 0.00284 100 0.084 % 55 120 126 8.6 0.00284 100 0.084 % 55 120 126 8.6 0.00284 100 0.084] % PD = 1800; % ndata = [36 114 94.705 6.7300 0.00690 100 0.0840; % 36 114 94.705 6.7300 0.00690 100 0.0840; % 60 120 309.540 7.0700 0.02028 100 0.0840; % 80 190 369.030 8.1800 0.00942 150 0.0630; % 47 97 148.890 5.3500 0.01140 120 0.0770; % 68 140 222.330 8.0500 0.01142 100 0.0840; % 110 300 287.710 8.0300 0.00357 200 0.0420; % 135 300 391.980 6.9900 0.00492 200 0.0420; % 135 300 455.760 6.6000 0.00573 200 0.0420; % 130 300 722.820 12.900 0.00605 200 0.0420; % 94 375 635.200 12.900 0.00515 200 0.0420; % 94 375 654.690 12.800 0.00569 200 0.0420; % 125 500 913.400 12.500 0.00421 300 0.0350; % 125 500 1760.400 8.8400 0.00752 300 0.0350; % 125 500 1728.300 9.1500 0.00708 300 0.0350; % 125 500 1728.300 9.1500 0.00708 300 0.0350; % 220 500 647.850 7.9700 0.00313 300 0.0350; % 220 500 649.690 7.9500 0.00313 300 0.0350; % 242 550 647.830 7.9700 0.00313 300 0.0350; % 242 550 647.810 7.9700 0.00313 300 0.0350; % 254 550 785.960 6.6300 0.00298 300 0.0350; % 254 550 785.960 6.6300 0.00298 300 0.0350; % 254 550 794.530 6.6600 0.00284 300 0.0350; % 254 550 794.530 6.6600 0.00284 300 0.0350; % 254 550 801.320 7.1000 0.00277 300 0.0350; % 254 550 801.320 7.1000 0.00277 300 0.0350; % 10 150 1055.100 3.3300 0.52124 120 0.0770; % 10 150 1055.100 3.3300 0.52124 120 0.0770; % 10 150 1055.100 3.3300 0.52124 120 0.0770; % 47 97 148.890 5.3500 0.01140 120 0.0770; % 60 190 222.920 6.4300 0.00160 150 0.0630; % 60 190 222.920 6.4300 0.00160 150 0.0630; % 60 190 222.920 6.4300 0.00160 150 0.0630; % 90 200 107.870 8.9500 0.00010 200 0.0420; % 90 200 116.580 8.6200 0.00010 200 0.0420; % 90 200 116.580 8.6200 0.00010 200 0.0420; % 25 110 307.450 5.8800 0.01610 80 0.0980; % 25 110 307.450 5.8800 0.01610 80 0.0980; % 25 110 307.450 5.8800 0.01610 80 0.0980; % 242 550 647.830 7.9700 0.00313 300 0.0350] % PD = 10500; % Xmin = ndata(:,1).'; % Xmax = ndata(:,2).'; % dim = size(Xmin, 2); % N = 50;% dimension % ci = ndata(:,3).' % bi = ndata(:,4).' % ai = ndata(:,5).' % ei = ndata(:,6).' % fi = ndata(:,7).' % alfa=10;G0=100; %random initialization for agents. X=initialization(PD,Xmax,Xmin,N, dim); %create the best so far chart and average fitnesses chart. BestChart=[];MeanChart=[]; V=zeros(N,dim); Fbest = 0; tPosMin = 1; tElapsed = 0; tStart = tic; for iteration=1:max_it % iteration %Checking allowable range. X=space_bound(PD,Xmax,Xmin,X); %Evaluation of agents. fitness=evaluateF(X, ai, bi, ci, ei, fi, Xmin); if min_flag==1 [best best_X]=min(fitness); %minimization. else [best best_X]=max(fitness); %maximization. end if iteration==1 Fbest=best;Lbest=X(best_X,:); end if min_flag==1 if best Fbest=best;Lbest=X(best_X,:); tPosMin = iteration; tElapsed = toc(tStart) end else if best>Fbest %maximization Fbest=best;Lbest=X(best_X,:); end end BestChart=[BestChart Fbest]; MeanChart=[MeanChart mean(fitness)]; %Calculation of M. eq.14-20 [M]=massCalculation(fitness,min_flag); %Calculation of Gravitational constant. eq.13. G=Gconstant(iteration,max_it,alfa,G0); %Calculation of accelaration in gravitational field. eq.7-10,21. a=Gfield(M,X,G,Rnorm,Rpower,ElitistCheck,iteration,max_it); %draw DrawGSA(Lbest, BestChart, tPosMin, tElapsed, AxesPower, AxesBestCost); %Agent movement. eq.11-12 [X,V]= move(X,a,V); end %iteration Code File Move: %This function updates the velocity and position of agents. function [X,V]=move(X,a,V) %movement. [N,dim]=size(X); V=rand(N,dim).*V+a; %eq. 11. X=X+V; %eq. 12. Code File Test_function_range: function [down,up,dim]=test_functions_range(F_index) %If lower bounds of dimensions are the same, then 'down' is a value. %Otherwise, 'down' is a vector that shows the lower bound of each dimension. %This is also true for upper bounds of dimensions. %Insert your own boundaries with a new F_index. dim=30; if F_index==1 down=-100;up=100; end if F_index==2 down=-10;up=10; end if F_index==3 down=-100;up=100; end if F_index==4 down=-100;up=100; end if F_index==5 down=-30;up=30; end if F_index==6 down=-100;up=100; end if F_index==7 down=-1.28;up=1.28; end if F_index==8 down=-500;up=500; end if F_index==9 down=-5.12;up=5.12; end if F_index==10 down=-32;up=32; end if F_index==11 down=-600;up=600; end if F_index==12 down=-50;up=50; end if F_index==13 down=-50;up=50; end if F_index==14 down=-65.536;up=65.536;dim=2; end if F_index==15 down=-5;up=5;dim=4; end if F_index==16 down=-5;up=5;dim=2; end if F_index==17 down=[-5 0];up=[10 15];dim=2; end if F_index==18 down=-2;up=2;dim=2; end if F_index==19 down=0;up=1;dim=3; end if F_index==20 down=0;up=1;dim=6; end if F_index==21 down=0;up=10;dim=4; end if F_index==22 down=0;up=10;dim=4; end if F_index==23 down=0;up=10;dim=4; end % Duong add if F_index==24 down = [0 0 0 60 60 60 60 60 60 40 40 55 55 ].'; up = [680 360 360 180 180 180 180 180 180 120 120 120 120].'; dim = 1; end Code File getOptions: function options = getOptions() options = struct(... ... %variables 'Vars', struct(... 'Dim', 0,...%dim: Dimension of the test function. 'NbOfAgents', 0,...%N: Number of agents. 'Rnorm', 2,...%Rnorm: Norm in eq.8. 'Rpower', 1,...%Rpower: Power of R in eq.7. 'MinFlag', 1,... 'ElitistCheck', 0,... 'Alpha', 10,... 'G0', 100,... 'Iterations', 1000,... 'PD', 0),... 'Parameters', struct()... ); Code File Fig: function varargout = fig(varargin) % FIG MATLAB code for fig.fig % FIG, by itself, creates a new FIG or raises the existing % singleton*. % H = FIG returns the handle to a new FIG or the handle to % the existing singleton*. % FIG('CALLBACK',hObject,eventData,handles,...) calls the local % function named CALLBACK in FIG.M with the given input arguments. % FIG('Property','Value',...) creates a new FIG or raises the % existing singleton*. Starting from the left, property value pairs are % applied to the GUI before fig_OpeningFcn gets called. An % unrecognized property name or invalid value makes property application % stop. All inputs are passed to fig_OpeningFcn via varargin. % *See GUI Options on GUIDE's Tools menu. Choose "GUI allows only one % instance to run (singleton)". % See also: GUIDE, GUIDATA, GUIHANDLES % Edit the above text to modify the response to help fig % Last Modified by GUIDE v2.5 19-Jul-2016 13:37:48 % Begin initialization code - DO NOT EDIT gui_Singleton = 1; gui_State = struct('gui_Name', mfilename, ... 'gui_Singleton', gui_Singleton, ... 'gui_OpeningFcn', @fig_OpeningFcn, ... 'gui_OutputFcn', @fig_OutputFcn, ... 'gui_LayoutFcn', [] , ... 'gui_Callback', []); if nargin && ischar(varargin{1}) gui_State.gui_Callback = str2func(varargin{1}); end if nargout [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:}); else gui_mainfcn(gui_State, varargin{:}); end % End initialization code - DO NOT EDIT % --- Executes just before fig is made visible. function fig_OpeningFcn(hObject, eventdata, handles, varargin) % This function has no output args, see OutputFcn. % hObject handle to figure % eventdata reserved - to be defined in a future version of MATLAB % handles structure with handles and user data (see GUIDATA) % varargin command line arguments to fig (see VARARGIN) % Choose default command line output for fig initialize_gui(hObject,handles,true); handles = guidata(hObject); handles.output = hObject; bTitle = get(handles.axespower,'Title'); set(bTitle, 'String', 'Power'); bxLabel = get(handles.axespower,'xLabel'); set(bxLabel, 'String', 'Power Units'); byLabel = get(handles.axespower,'yLabel'); set(byLabel, 'String', 'Power'); pTitle = get(handles.axesbestcost,'Title'); set(pTitle, 'String', 'Best Fuel Cost(GSA)'); pxLabel = get(handles.axesbestcost,'xLabel'); set(pxLabel, 'String', 'Iteration'); pyLabel = get(handles.axesbestcost,'yLabel'); set(pyLabel, 'String', 'Fuel Cost($/h)'); % Update handles structure guidata(hObject, handles); % UIWAIT makes fig wait for user response (see UIRESUME) % uiwait(handles.figure1); % --- Outputs from this function are returned to the command line. function varargout = fig_OutputFcn(hObject, eventdata, handles) % varargout cell array for returning output args (see VARARGOUT); % hObject handle to figure % eventdata reserved - to be defined in a future version of MATLAB % handles structure with handles and user data (see GUIDATA) % Get default command line output from handles structure varargout{1} = handles.output; function update_display(hObject,handles, data) handles.options.Parameters = data.Parameters; handles.options.Vars.Iterations = data.Iterations; handles.options.Vars.Alpha = data.Alpha; handles.options.Vars.NbOfAgents = data.NbOfAgents; handles.options.Vars.G0 = data.G0; handles.options.Vars.PD = data.PD; handles.options.Vars.Dim = size(handles.options.Parameters.Pmin, 2); set(handles.iteration, 'String', handles.options.Vars.Iterations); set(handles.dimension, 'String', handles.options.Vars.NbOfAgents); set(handles.alpha, 'String', handles.options.Vars.Alpha); set(handles.G0, 'String', handles.options.Vars.G0); set(handles.PD, 'String', handles.options.Vars.PD); data = [[1:handles.options.Vars.Dim]',... handles.options.Parameters.Pmin',... handles.options.Parameters.Pmax',... handles.options.Parameters.a',... handles.options.Parameters.b',... handles.options.Parameters.c',... handles.options.Parameters.e',... handles.options.Parameters.f']; cName = {'units','Pmin', 'Pmax', 'a', 'b', 'c', 'e', 'f'}; set(handles.uitable1,'ColumnWidth',{60},'data', data, 'ColumnName', cName) guidata(hObject, handles); % -------------------------------------------------------------------- function initialize_gui(hObject, handles, isreset) % If the metricdata field is present and the reset flag is false, it means % we are we are just re-initializing a GUI by calling it from the cmd line % while it is up. So, bail out as we dont want to reset the data. handles = guidata(hObject); if isfield(handles, 'options') && ~isreset return; end handles.options = getOptions(); data = load('13_units_data_test.mat'); set(handles.units, 'SelectedObject', handles.thirteenunits); update_display(hObject,handles, data.data); handles = guidata(hObject); % Update handles structure guidata(hObject, handles); function iteration_Callback(hObject, eventdata, handles) % hObject handle to iteration (see GCBO) % eventdata reserved - to be defined in a future version of MATLAB % handles structure with handles and user data (see GUIDATA) % Hints: get(hObject,'String') returns contents of iteration as text % str2double(get(hObject,'String')) returns contents of iteration as a double iterations = str2double(get(hObject, 'String')); if isnan(iterations) set(hObject, 'String', 0); errordlg('Input must be a number','Error'); end handles = guidata(hObject); % Save the new density value handles.options.Vars.Iterations = iterations; guidata(hObject,handles); % --- Executes during object creation, after setting all properties. function iteration_CreateFcn(hObject, eventdata, handles) % hObject handle to iteration (see GCBO) % eventdata reserved - to be defined in a future version of MATLAB % handles empty - handles not created until after all CreateFcns called % Hint: edit controls usually have a white background on Windows. % See ISPC and COMPUTER. if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor')) set(hObject,'BackgroundColor','white'); end function dimension_Callback(hObject, eventdata, handles) % hObject handle to dimension (see GCBO) % eventdata reserved - to be defined in a future version of MATLAB % handles structure with handles and user data (see GUIDATA) % Hints: get(hObject,'String') returns contents of dimension as text % str2double(get(hObject,'String')) returns contents of dimension as a double dimension = str2double(get(hObject, 'String')); if isnan(dimension) set(hObject, 'String', 0); errordlg('Input must be a number','Error'); end handles = guidata(hObject); % Save the new density value handles.options.Vars.NbOfAgents = dimension; guidata(hObject,handles) % --- Executes during object creation, after setting all properties. function dimension_CreateFcn(hObject, eventdata, handles) % hObject handle to dimension (see GCBO) % eventdata reserved - to be defined in a future version of MATLAB % handles empty - handles not created until after all CreateFcns called % Hint: edit controls usually have a white background on Windows. % See ISPC and COMPUTER. if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor')) set(hObject,'BackgroundColor','white'); end function alpha_Callback(hObject, eventdata, handles) % hObject handle to alpha (see GCBO) % eventdata reserved - to be defined in a future version of MATLAB % handles structure with handles and user data (see GUIDATA) % Hints: get(hObject,'String') returns contents of alpha as text % str2double(get(hObject,'String')) returns contents of alpha as a double alpha = str2double(get(hObject, 'String')); if isnan(alpha) set(hObject, 'String', 0); errordlg('Input must be a number','Error'); end handles = guidata(hObject); % Save the new density value handles.options.Vars.Alpha = alpha; guidata(hObject,handles) % --- Executes during object creation, after setting all properties. function alpha_CreateFcn(hObject, eventdata, handles) % hObject handle to alpha (see GCBO) % eventdata reserved - to be defined in a future version of MATLAB % handles empty - handles not created until after all CreateFcns called % Hint: edit controls usually have a white background on Windows. % See ISPC and COMPUTER. if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor')) set(hObject,'BackgroundColor','white'); end function G0_Callback(hObject, eventdata, handles) % hObject handle to G0 (see GCBO) % eventdata reserved - to be defined in a future version of MATLAB % handles structure with handles and user data (see GUIDATA) % Hints: get(hObject,'String') returns contents of G0 as text % str2double(get(hObject,'String')) returns contents of G0 as a double G0 = str2double(get(hObject, 'String')); if isnan(G0) set(hObject, 'String', 0); errordlg('Input must be a number','Error'); end handles = guidata(hObject); % Save the new density value handles.options.Vars.G0 = G0; guidata(hObject,handles) % --- Executes during object creation, after setting all properties. function G0_CreateFcn(hObject, eventdata, handles) % hObject handle to G0 (see GCBO) % eventdata reserved - to be defined in a future version of MATLAB % handles empty - handles not created until after all CreateFcns called % Hint: edit controls usually have a white background on Windows. % See ISPC and COMPUTER. if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor')) set(hObject,'BackgroundColor','white'); end function PD_Callback(hObject, eventdata, handles) % hObject handle to PD (see GCBO) % eventdata reserved - to be defined in a future version of MATLAB % handles structure with handles and user data (see GUIDATA) % Hints: get(hObject,'String') returns contents of PD as text % str2double(get(hObject,'String')) returns contents of PD as a double PD = str2double(get(hObject, 'String')); if isnan(PD) set(hObject, 'String', 0); errordlg('Input must be a number','Error'); end handles = guidata(hObject); % Save the new density value handles.options.Vars.PD = PD; guidata(hObject,handles) % --- Executes during object creation, after setting all properties. function PD_CreateFcn(hObject, eventdata, handles) % hObject handle to PD (see GCBO) % eventdata reserved - to be defined in a future version of MATLAB % handles empty - handles not created until after all CreateFcns called % Hint: edit controls usually have a white background on Windows. % See ISPC and COMPUTER. if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor')) set(hObject,'BackgroundColor','white'); end % --- Executes on button press in runstop. function runstop_Callback(hObject, eventdata, handles) % hObject handle to runstop (see GCBO) % eventdata reserved - to be defined in a future version of MATLAB % handles structure with handles and user data (see GUIDATA) %data = get(handles.figure1,'UserData'); handles = guidata(hObject); [Fbest,Lbest] = GSA(handles.options, handles.axespower,handles.axesbestcost); cName = 'Power'; f1 = @(x) sprintf('P %d',x); rNames = cellfun(f1, num2cell(1:size(Lbest,2)),'UniformOutput', false); rNames{size(Lbest,2)+1} = 'Total'; data = [Lbest sum(Lbest)]; set(handles.uitable2,'ColumnWidth',{75},'data',data','ColumnName',cName,... 'RowName', rNames); %legend(handles.axesbestcost,sprintf('Best Cost = %0.3f',Fbest)) % --- Executes during object creation, after setting all properties. function units_CreateFcn(hObject, eventdata, handles) % hObject handle to units (see GCBO) % eventdata reserved - to be defined in a future version of MATLAB % handles empty - handles not created until after all CreateFcns called % --- Executes when selected object is changed in units. function units_SelectionChangeFcn(hObject, eventdata, handles) % hObject handle to the selected object in units % eventdata structure with the following fields (see UIBUTTONGROUP) % EventName: string 'SelectionChanged' (read only) % OldValue: handle of the previously selected object or empty if none was selected % NewValue: handle of the currently selected object % handles structure with handles and user data (see GUIDATA) handles = guidata(hObject); if (hObject == handles.thirteenunits) data = load('13_units_data_test.mat') update_display(hObject,handles, data.data); % Update handles structure handles = guidata(hObject); guidata(hObject, handles); elseif (hObject == handles.fourtyunits) data = load('40_units_data_test.mat') update_display(hObject,handles, data.data); % Update handles structure handles = guidata(hObject); guidata(hObject, handles); else data = load('3_units_data_test.mat'); update_display(hObject,handles, data.data); % Update handles structure handles = guidata(hObject); guidata(hObject, handles); end % -------------------------------------------------------------------- function Untitled_1_Callback(hObject, eventdata, handles) % hObject handle to Untitled_1 (see GCBO) % eventdata reserved - to be defined in a future version of MATLAB % handles structure with handles and user data (see GUIDATA) % -------------------------------------------------------------------- function Untitled_2_Callback(hObject, eventdata, handles) % hObject handle to Untitled_2 (see GCBO) % eventdata reserved - to be defined in a future version of MATLAB % handles structure with handles and user data (see GUIDATA) % -------------------------------------------------------------------- function Untitled_3_Callback(hObject, eventdata, handles) % hObject handle to Untitled_3 (see GCBO) % eventdata reserved - to be defined in a future version of MATLAB % handles structure with handles and user data (see GUIDATA) Code File Draw GSA: %DrawSwarm >> Internal function of psotoolbox. % Purpose: To draw a visual display of the Swarm. % You shouldn't need to mess around with this fn. if u don't wanna change the visualization. % see also: pso.m function DrawGSA(X, Fbest, tPosMin, tElapsed, AxesPower, AxesBestCost) bar(AxesPower,X); plot(AxesBestCost,Fbest); xl = xlim(AxesBestCost); yl = ylim(AxesBestCost); x = xl(1) + (xl(2) - xl(1))/2; y = yl(1) + (yl(2) - yl(1))/2; strtelapsed = ['time/iter=', num2str(tElapsed), ' / ', num2str(tPosMin), '=', num2str(tElapsed/tPosMin)]; strbestcost = ['best cost=', num2str(Fbest(tPosMin))]; str = {strtelapsed,strbestcost}; text(x, y, str,'VerticalAlignment', 'middle', 'HorizontalAlignment','center'); % bTitle = get(AxesPower,'Title'); % set(bTitle, 'String', 'Power'); % bxLabel = get(AxesPower,'xLabel'); % set(bxLabel, 'String', 'Power xLabel'); % byLabel = get(AxesPower,'yLabel'); % set(byLabel, 'String', 'Power yLabel'); % pTitle = get(AxesBestCost,'Title'); % set(pTitle, 'String', 'Best Cost'); % pxLabel = get(AxesBestCost,'xLabel'); % set(pxLabel, 'String', 'Best Cost xLabel'); % pyLabel = get(AxesBestCost,'yLabel'); % set(pyLabel, 'String', 'Best Cost yLabel'); drawnow;
