ĐẠI HỌC THÁI NGUYÊN
TRƯỜNG ĐẠI HỌC KHOA HỌC -------------------------------
NGUYỄN HỮU ĐẠT
MỘT SỐ THUẬT TOÁN GIẢI SỐ BÀI TOÁN TỐI ƯU PHI TUYẾN.
LUẬN VĂN THẠC SĨ TOÁN HỌC
THÁI NGUYÊN - 2019
ĐẠI HỌC THÁI NGUYÊN
TRƯỜNG ĐẠI HỌC KHOA HỌC -------------------------------
NGUYỄN HỮU ĐẠT
MỘT SỐ THUẬT TOÁN GIẢI SỐ BÀI TOÁN TỐI ƯU PHI TUYẾN.
Chuyên ngành: Toán ứng dụng Mã số : 8 46 01 12
LUẬN VĂN THẠC SĨ TOÁN HỌC NGƯỜI HƯỚNG DẪN KHOA HỌC
TS. Vũ Vinh Quang
THÁI NGUYÊN - 2019
i
Lời cảm ơn
Trước hết, em xin bày tỏ lòng kính trọng và lòng biết ơn sâu sắc tới TS.
Vũ Vinh Quang, người thầy đã tận tình hướng dẫn, chỉ bảo và cung cấp
những tài liệu rất hữu ích để em có thể hoàn thành luận văn.
Xin cảm ơn lãnh đạo Trường Đại học Khoa học - Đại học Thái Nguyên
đã tạo điều kiện giúp đỡ tôi về mọi mặt trong suốt quá trình học tập và
thực hiện luận văn.
Em xin bày tỏ lòng biết ơn tới các thầy, cô giáo giảng dạy lớp K11C
đã truyền đạt kiến thức và phương pháp nghiên cứu khoa học trong suốt
những năm học vừa qua.
Xin chân thành cảm ơn anh chị em học viên cao học K11C và bạn bè
đồng nghiệp đã động viên và khích lệ tôi trong quá trình học tập, nghiên
cứu và làm luận văn.
Tôi xin bày tỏ lòng biết ơn sâu sắc đến gia đình, người thân, những
người luôn động viên, khuyến khích và giúp đỡ về mọi mặt để tôi có thể
hoàn thành công việc nghiên cứu.
Thái Nguyên, tháng 4 năm 2019
Tác giả luận văn
Nguyễn Hữu Đạt
ii
Lời cam đoan
Tôi xin cam đoan:
Những nội dung trong luận văn này là do tôi thực hiện dưới sự hướng
dẫn trực tiếp của thầy giáo hướng dẫn TS. Vũ Vinh Quang.
Mọi tham khảo dùng trong luận văn đều được trích dẫn rõ ràng tác giả,
tên công trình, thời gian, địa điểm công bố.
Tôi xin chịu trách nhiệm với lời cam đoan của mình.
Thái Nguyên, tháng 4 năm 2019
Tác giả luận văn
Nguyễn Hữu Đạt
iii
Mục lục
Lời cảm ơn i
Lời cam đoan ii
Bảng ký hiệu v
Mở đầu 1
3 1 Một số kiến thức cơ bản
3 1.1. Mô hình tổng quát của bài toán tối ưu hóa . . . . . . . . .
4 1.2. Phân loại bài toán tối ưu . . . . . . . . . . . . . . . . . . .
5 1.3. Một số phương pháp giải cơ bản bài toán tuyến tính . . .
5 1.3.1. Thuật toán hình học . . . . . . . . . . . . . . . . .
6 1.4. Mô hình bài toán quy hoạch lồi tổng quát . . . . . . . . .
6 1.4.1. Khái niệm tập lồi, hàm lồi . . . . . . . . . . . . . .
8 1.4.2. Khái niệm về Gradient và đạo hàm hướng . . . . .
8 1.4.3. Bài toán quy hoạch lồi tổng quát, điều kiện tối ưu .
1.4.4. Cực tiểu hàm lồi một biến . . . . . . . . . . . . . . 10
1.5. Phương pháp giải bài toán quy hoạch tuyến tính tổng quát
trên phần mềm MATLAB . . . . . . . . . . . . . . . . . . 15
2 Một số thuật toán giải số bài toán tối ưu phi tuyến không
16 ràng buộc
iv
2.1. Một số kiến thức cơ bản . . . . . . . . . . . . . . . . . . . 16
2.1.1. Định nghĩa . . . . . . . . . . . . . . . . . . . . . . 16
2.1.2. Điều kiện tối ưu . . . . . . . . . . . . . . . . . . . . 17
2.2. Các thuật toán sử dụng đạo hàm . . . . . . . . . . . . . . 18
2.2.1. Thuật toán Gradient . . . . . . . . . . . . . . . . . 18
2.2.2. Thuật toán đường dốc nhất . . . . . . . . . . . . . 20
2.2.3. Thuật toán Newton . . . . . . . . . . . . . . . . . . 23
2.3. Các thuật toán không sử dụng đạo hàm . . . . . . . . . . 26
2.3.1. Phương pháp tìm trực tiếp (Direct search) . . . . . 26
2.3.2. Phương pháp Powell . . . . . . . . . . . . . . . . . 27
2.3.3. Phương pháp Nelder và Mead . . . . . . . . . . . . 28
3 Một số thuật toán giải số bài toán tối ưu phi tuyến có ràng
buộc 32
3.1. Một số kiến thức cơ bản . . . . . . . . . . . . . . . . . . . 32
3.1.1. Hàm Lagrange . . . . . . . . . . . . . . . . . . . . 32
3.1.2. Thiết lập điều kiện tối ưu Kuhn - Tucker . . . . . . 33
3.2. Một số thuật toán . . . . . . . . . . . . . . . . . . . . . . . 35
3.2.1. Thuật toán Gradient . . . . . . . . . . . . . . . . . 35
3.2.2. Phương pháp hàm phạt (Penalty function method) 38
Kết luận 52
Tài liệu tham khảo 53
v
Bảng ký hiệu
QHTT Quy hoạch tuyến tính;
tọa độ thứ i của x;
xi xT vectơ hàng (chuyển vị của x );
||x|| = chuẩn Euclide của x;
∂f (x) dưới vi phân của f tại x;
∇f (x) hoặc đạo hàm của f tại x; f (cid:48)(x) đạo hàm f tại x
1
Mở đầu
Mô hình bài toán quy hoạch phi tuyến nói chung là một mô hình quan
trọng trong lớp các bài toán tối ưu hóa. Mô hình này có rất nhiều ứng
dụng trong các bài toán cơ học và vật lý, kinh tế và thương mại. Về mặt lý
thuyết, đã có rất nhiều các tài liệu trình bày các thuật toán lý thuyết giải
mô hình các bài toán này trên mô hình tổng quát. Tuy nhiên việc nghiên
cứu và cài đặt chi tiết các thuật toán và ứng dụng vào một số mô hình đối
với một số bài toán cụ thể trong cơ học và vật lý là chưa nhiều người đề
cập đến.
Nội dung chính của luận văn là nghiên cứu cơ sở toán học của các thuật
toán cơ bản giải bài toán quy hoạch phi tuyến tính không ràng buộc và
có ràng buộc, tìm hiểu chi tiết các bước mô tả thuật toán, xây dựng sơ đồ
khối và cài đặt các thuật toán trên ngôn ngữ lập trình cụ thể.
Nội dung luận văn gồm ba chương, phần phụ lục được cấu trúc như sau:
Chương 1. Một số kiến thức cơ bản
Trình bày mô hình tổng quát của bài toán tối ưu hóa, phân loại bài toán tối
ưu, các phương pháp biến đổi cơ bản, một số thuật toán giải bài toán tối ưu
hàm lồi một biến, giải bài toán quy hoạch tuyến tính trên MATLAB. Các
kết quả là những kiến thức quan trọng được ứng dụng trong các chương
sau của luận văn.
Chương 2. Mô hình bài toán tối ưu hóa phi tuyến
Trong chương này, trình bày một số thuật toán giải số bài toán tối ưu phi
2
tuyến không ràng buộc. Cụ thể của chương nêu mô hình tổng quát, điều
kiện tối ưu, các thuật toán sử dụng đạo hàm như thuật toán Gradient,
thuật toán đường dốc nhất, thuật toán Newton, thuật toán Gradient liên
hợp, các thuật toán không sử dụng đạo hàm như thuật toán tìm trực tiếp,
thuật toán Powell, thuật toán Nelder và Mead.
Chương 3. Một số thuật toán giải số bài toán tối ưu phi tuyến có ràng
buộc
Nội dung của chương là tìm hiểu một số thuật toán giải số bài toán tối ưu
phi tuyến có ràng buộc như khái niệm hàm Lagrange, phương pháp hàm
phạt, các thuật toán tìm nghiệm xấp xỉ. Các thuật toán đã được cài đặt
trên môi trường MATLAB version 7.0.
3
Chương 1
Một số kiến thức cơ bản
Nội dung chính của chương là trình bày mô hình tổng quát của bài toán
tối ưu hóa, phân loại bài toán tối ưu, các phương pháp biến đổi cơ bản,
một số thuật toán giải bài toán tối ưu hàm lồi một biến, giải bài toán quy
hoạch tuyến tính trên MATLAB. Các kết quả là các kiến thức quan trọng
được ứng dụng trong các chương sau của luận văn. Các kiến thức được
1.1. Mô hình tổng quát của bài toán tối ưu hóa
tham khảo trong các tài liệu [1], [2], [4].
Tối ưu hóa là một trong những lĩnh vực quan trọng của bài toán có ảnh
hưởng đến hầu hết các lĩnh vực khoa học, công nghệ, kinh tế và xã hội.
Việc tìm giải pháp tối ưu cho một bài toán thực tế nào đó chiếm một vai
trò hết sức quan trọng như việc tiến hành lập kế hoạch sản xuất hay thiết
kế hệ thống điều khiển các quá trình... Nếu sử dụng các kiến thức trên
nền tảng của toán học để giải quyết các bài toán cực trị, người ta sẽ đạt
được hiệu quả kinh tế rất cao. Điều này phù hợp với mục đích của các bài
toán đặt ra trong thực tế hiện nay.
Mô hình bài toán tối ưu tổng quát được phát biểu như sau:
4
Cực đại hóa (cực tiểu hóa) hàm:
f (X) → max/min
Với các điều kiện:
(1.1) gi(X) = bi, i ∈ J1
(1.2) gj(X) ≤ bj, j ∈ J2
(1.3) gk(X) ≥ bk, k ∈ J3
(1.4) x1, x2, ..., xn ≥ 0.
Trong đó f (X) được gọi là hàm mục tiêu, các điều kiện (1.1) được
gọi là ràng buộc đẳng thức. Các điều kiện (1.2), (1.3) được gọi là ràng
buộc bất đẳng thức. Các điều kiện (1.4) được gọi là ràng buộc về dấu. X = (x1, x2, ..., xn)T là vectơ thuộc không gian Rn. Tập các vectơ X thỏa mãn hệ ràng buộc lập nên một miền D được gọi là miền phương án (hay
miền chấp nhận được), mỗi điểm X ∈ D gọi là một phương án. Một phương án X ∗ ∈ D làm cho hàm mục tiêu f (X) đạt cực đại hoặc cực tiểu
1.2. Phân loại bài toán tối ưu
được gọi là phương án tối ưu.
Dựa trên mô hình tổng quát, người ta thường phân loại lớp các bài toán
tối ưu như sau:
- Quy hoạch tuyến tính: Là những bài toán mà hàm mục tiêu f (X)
và tất cả các hàm ràng buộc gi(X), gj(X), gk(X) là tuyến tính.
- Quy hoạch phi tuyến: Là những bài toán một trong hàm mục tiêu
f (X) hoặc các hàm ràng buộc gi(X), gj(X), gk(X) là phi tuyến.
- Quy hoạch lồi: Là các bài toán quy hoạch mà các hàm mục tiêu
f (X) là lồi trên tập các ràng buộc D lồi.
5
- Quy hoạch lõm: Là các bài toán quy hoạch mà các hàm mục tiêu
f (X) là lõm trên tập các ràng buộc D lõm.
- Quy hoạch rời rạc: Bài toán tối ưu được gọi là quy hoạch rời rạc
nếu miền ràng buộc D là tập hợp rời rạc. Trong trường hợp riêng khi các
biến chỉ nhận giá trị nguyên thì ta có quy hoạch nguyên.
- Quy hoạch đa mục tiêu: Nếu trên cùng một miền ràng buộc ta
xét đồng thời các hàm mục tiêu khác nhau.
Trong các lĩnh vực kinh tế kỹ thuật thì quy hoạch phi tuyến, quy hoạch
1.3. Một số phương pháp giải cơ bản bài toán tuyến tính
tuyến tính là những bài toán thường gặp.
1.3.1. Thuật toán hình học
Xét bài toán:
f (x1, x2) = c1x1 + c2x2 → M ax;
a11x1 + a12x2 ≤ b1;
a21x1 + a22x2 ≤ b2;
...
an1x1 + an2x2 ≤ bn;
x1 ≥ 0; x2 ≥ 0.
Nhận xét
1. Vì các ràng buộc của bài toán luôn luôn là các nửa mặt phẳng, do
đó miền phương án luôn luôn là một đa giác lồi (là giao của các nửa
mặt phẳng).
2. Xét đường thẳng f = m được gọi là đường mức. Hiển nhiên khi
đường mức chuyển động song song trong miền phương án thì điểm chạm
6
cuối cùng của đường mức với miền luôn luôn là một trong các đỉnh của đa
giác (hoặc một cạnh của đa giác). Đấy chính là phương án tối ưu cần tìm.
Xuất phát từ nhận xét trên, chúng ta có thuật toán hình học gồm các
bước như sau:
Thuật toán:
Bước 1: Vẽ miền phương án D là đa giác lồi bằng cách xác định miền
giao của các nửa mặt phẳng trong hệ ràng buộc.
Bước 2 : Xác định tọa độ của các đỉnh đa giác: Giả sử là các điểm
A1, A2, ..., Ak.
Bước 3: Xác định phương án tối ưu
fmax = max(f (A1), f (A2), ..., f (Ak)).
Chú ý
Trong trường hợp khi miền phương án không phải là miền kín thì tùy
thuộc vào hướng di chuyển của đường mức, chúng ta sẽ xác định được
1.4. Mô hình bài toán quy hoạch lồi tổng quát
phương án tối ưu của bài toán.
1.4.1. Khái niệm tập lồi, hàm lồi
a. Tập lồi
Định nghĩa: Tập C ⊂ Rn được gọi là tập lồi nếu x, y ∈ C ⇒ λx + (1 − λ)y ∈ C,
∀λ ∈ [0; 1]. Nghĩa là nếu x, y ∈ C thì đoạn thẳng [x, y] ∈ C
b. Hàm lồi
Định nghĩa:
Hàm số f (x) là lồi (convex function) trên tập C nếu với mọi cặp điểm
(x1, x2) thuộc C và mọi số λ ∈ [0, 1], ta có:
f [λx1 + (1 − λ)x2] ≤ λf (x1) + (1 − λ)f (x2)
7
có nghĩa là điểm x = λx1 + (1 − λ)x2 trong [x1, x2] thì mọi điểm của đồ thị luôn nằm dưới M1M2
Điểm A: x = λx1 + (1 − λ)x2 Điểm B: f (x) = λf (x1) + (1 − λ)f (x2) Điểm C: f (x) = f (λx1 + (1 − λ)x2) Một số điều kiện
- Hàm f (x) là lồi, nếu đối với hai điểm x1, x2 thỏa mãn điều kiện f (x2) ≥ f (x1) + (cid:53)f (cid:48)(x1).(x2 − x1) - Hàm f (x) là hàm lồi, nếu ma trận Hesian H(x) = [∂2f (x)/∂x2∂x2] là bán xác định dương. Khi H(x) xác định dương thì hàm f (x) gọi là hàm
lồi chặt.
Cực trị của hàm lồi
Bất cứ cực tiểu địa phương nào của hàm lồi trên tập lồi cũng là cực tiểu
của hàm trên tập đó.
Ta sẽ chứng minh tính chất này bằng phản chứng:
(cid:48)
Giả thiết hàm f (x) có hai điểm cực tiểu là x1 và x2. Vì f (x) là lồi nên:
f (x2) − f (x1) ≥ (cid:53)f (x2)(x2 − x1)
Hay (cid:53)f (cid:48)(xs).s ≤ 0 Trong đó s = (x2 − x1) là vectơ nối điểm x1 và x2. Theo (*), hàm f (x) giảm khi di chuyển theo hướng s khi xuất phát từ điểm x1. Điều này trái
8
với giả thiết x1 là cực tiểu. Do đó f (x) chỉ có một cực tiểu duy nhất.
Như vậy trong quy hoạch lồi thì giá trị tối ưu địa phương cũng là giá
trị tối ưu toàn cục.
1.4.2. Khái niệm về Gradient và đạo hàm hướng
+ Gradient của f (x) là một vectơ có các thành phần là đạo hàm riêng
∂f (x)/∂x1 (cid:21)T , , ..., (cid:53)f (x) = (cid:20) ∂f ∂x1 ∂f ∂x2 ∂f ∂xn
+ Vectơ (cid:53)f (x0) vuông góc với đường mức của f (x) tại x0, tốc độ hàm f (x) tăng nhanh nhất theo hướng của (cid:53)f (x0). Nếu đi theo hướng − (cid:53) f (x0) thì f (x) giảm nhanh nhất.
(cid:48)
+ Đạo hàm theo hướng z của hàm f (x) tại điểm x0:
z(x0) = (cid:104)(cid:53)f (x0), z(cid:105) = | (cid:53) f (x0)|.|z|.cos((cid:53)f (x0), z)
f
Đó là hình chiếu của vectơ (cid:53)f (x0) lên hướng z. + Ma trận Hessian H(x) là ma trận có các thành phần là Gradient cấp
hai của f (x)
...
∂2f ∂x1∂x1 ∂2f ∂x2∂x1 ... ∂2f ∂xn∂x1
∂2f ∂x1∂x2 ∂2f ∂x2∂x2 ... ∂2f ∂xn∂x2
∂2f ∂x1∂xn ∂2f ∂x2∂xn ... ∂2f ∂xn∂xn
... H(x) = ∇2f (x) = ... ...
1.4.3. Bài toán quy hoạch lồi tổng quát, điều kiện tối ưu
a. Phát biểu bài toán
Tìm x sao cho hàm mục tiêu
f (x) → min
Các ràng buộc:
x ∈ C : gi(x) ≤ 0; i = 1, 2, ..., m.
9
trong đó C là tập lồi, f , gi là các hàm lồi trên C. b. Điều kiện tối ưu
+ Miền nghiệm chấp nhận được:
D = {x ∈ C : gi(x) ≤ 0} : i = 1, 2, ..., m.
+ Khái niệm về điểm yên ngựa (saddle point)
m (cid:80) i=1
Ký hiệu hàm Lagrange L(x, λ) = f (x) + λigi(x) = f (x) + λT g(x)
Khi đó điểm yên ngựa của hàm L(x, λ) là điểm (x∗, λ∗) với x∗ ∈ D; λ∗ ≥ 0 sao cho: L(x, λ∗) ≤ L(x∗, λ∗) ≤ L(x∗, λ). Các thành phần λi của vectơ λ = [λ1, ..., λm]T được gọi là các nhân tử Lagrange.
Khi λ = λ∗ thì điểm (x∗, λ∗) là điểm cao nhất của L(x, λ). Khi x = x∗ thì điểm (x∗, λ∗) là điểm thấp nhất của L(x, λ).
Phương pháp xác định điểm yên ngựa: Điểm (x∗, λ∗) là điểm yên ngựa của hàm L(x∗, λ∗) khi và chỉ khi:
+(cid:79)L(x, λ) = 0
+gi(x) ≤ 0; i = 1, 2, ..., m
+λigi(x) ≤ 0; i = 1, 2, ..., m
+ Điều kiện cần và đủ của tối ưu
Định lí 1.1
Điểm x∗ là tối ưu khi và chỉ khi
fz(x∗) = (cid:104)(cid:53)f (x∗), z(cid:105) ≥ 0; ∀z ∈ D(x∗)
Tức là nếu xuất phát từ x∗ theo hướng bất kỳ z mà f (x) tăng thì f (x) đạt giá trị min tại x∗.
Định lí 1.2 (Định lý Kuhn - Tucker, phát biểu năm 1951)
Giả sử bài toán quy hoạch lồi thỏa mãn điều kiện Slater:
∃x0 ∈ C : gi(x0) ≤ 0; i = 1, 2, ..., m
10
2, ..., λ∗
1, λ∗ ngựa của hàm Lagrange L(x, λ).
Khi đó điều cần và đủ để x∗ trở thành nghiệm tối ưu là tồn tại một vectơ m]T sao cho cặp (x∗, λ∗) là điểm yên m chiều, không âm λ∗ = [λ∗
Chú ý Điều kiện Slater không được thỏa mãn thì có thể không tồn tại điểm yên ngựa của hàm L(x, λ) loại (x∗, λ∗).
Ví dụ 1 Tìm min của f (x) = −x với ràng buộc g = x2 < 0 Ta có x∗ = 0. Hàm Lagrange L(x, λ) = −x + λ.x2
Xác định điểm dừng:
= −1 + 2λ.x = 0 ∂L ∂x
Điều kiện Slater không thỏa mãn khi x = 0, do đó không tồn tại λ và hàm
L(x, λ) không có điểm yên ngựa loại (0, λ).
1.4.4. Cực tiểu hàm lồi một biến
Thuật toán chia đôi
Cho hàm số f (x) xác định trên đoạn [a, b] với điều kiện f (x) lồi trên
[a, b]. Cần xác định điểm xopt để hàm f (x) đạt min.
Tư tưởng: Thu hẹp dần miền tìm kiếm để xác định được xopt
Bước 1: Chọn hai điểm x1, x2 ∈ (a, b) bất kỳ Bước 2: Nếu f (x1) < f (x2) thì xopt ∈ [a, x2] Nếu f (x1) > f (x2) thì xopt ∈ [x1, b] Nếu f (x1) = f (x2) thì xopt ∈ [x1, x2]
Thuật toán chia đôi:
Bước 1: Đặt L = (b − a)
Bước 2: While L > ε
11
4
4 , x2 = b − L
x1 = a + L f (x1) < f (x2) thì b := x2 f (x1) > f (x2) thì a := x1 f (x1) = f (x2) thì a := x1, b := x2 L = b − a;
Lặp lại cho đến khi L < ε
Nhận xét
+ Thuật toán chắc chắn hội tụ sau hữu hạn bước lặp, tốc độ hội tụ phụ
thuộc vào việc chọn vị trí của hai điểm x1, x2 sau từng bước lặp.
+ Phương pháp chia đôi dễ lập trình nhưng không tối ưu về tính toán
vì việc chọn x1, x2 là tùy ý.
Thuật toán được mô phỏng dưới ngôn ngữ lập trình MATLAB.
function cd=chia_doi(a, b, epxilon)
L=b-a;
count = 0;
while L > epxilon
count=count + 1
x1=a + L/4;
x2=b - L/4;
if f(x1) < f(x2)
b = x2;
elseif f(x1) > f(x2)
a = x1;
elseif f(x1) == f(x2)
a = x1;
b = x2;
end;
L= b - a;
end;
12
cd = x1;
Sau đây là kết quả chạy thuật toán giải bài toán
f (x) = x2 + , x ∈ [−10, 10].
54 x Hiển nhiên hàm f (x) là hàm lồi trên đoạn [1, 10], cực tiểu đạt tại x=3.
Kết quả chạy thuật toán chia đôi cho trong bảng 2.1
Bảng 2.1: Nghiệm xấp xỉ tối ưu sau các bước lặp
Số bước lặp Nghiệm tối ưu Số bước lặp Nghiệm tối ưu
5 3.6738 35 2.9938
10 3.3165 40 2.9999
15 3.1148 45 2.9982
20 3.0113 50 3.0001
25 3.0275 55 3.0004
30 3.0003 60 3.0000
Nhận xét Thuật toán hội tụ tốt sau 60 bước lặp đạt đến nghiệm
tối ưu.
Thuật toán mặt cắt vàng
Phương pháp này sử dụng tính chất của dãy Fibonacci {Fn}.
F1 = F2 = 1 (n=3, 4,...) Fn = Fn−1 + Fn−2
Dãy số Fibonacci có tính chất đặc biệt đáng chú ý là: Tỷ số giữa hai số kế
tiếp nhau của dãy tiến tới tỷ số vàng: √
= lim n→∞ 5 − 1 2 Fn−1 Fn
+ Tư tưởng phương pháp
Đặt bk − ak = 1 và chia [ak, bk] theo tỷ lệ:
= ⇔ p2 + p − 1 = 0 1 − p p 1 − 2(1 − p) 1 − p
13
5
√ 2 ≈ 0, 61803
Phương trình p2 + p − 1 = 0 có nghiệm dương: p = −1+
Hai số p và (1-p)=0,38197 được gọi là hằng số Fibonacci
+ Trong đoạn [ak, bk] lấy hai điểm αk, βk theo tỷ lệ:
= = p (*) βk − ak bk − ak bk − αk bk − ak
Nếu chia [ak, bk] thành hai đoạn thỏa mãn (*) thì ở mỗi phép lặp chỉ cần tính một giá trị của f(x). Do đó việc tìm giá trị tối ưu x∗ sẽ nhanh hơn so
với việc chia đều [ak, bk] thành những đoạn bằng nhau.
Các mặt cắt [αk, βk] thỏa mãn (*) được gọi là mặt cắt vàng và p là
số vàng.
Thuật toán lát cắt vàng
+ Chọn αk, βk theo tỷ lệ vàng (*). + Tính giá trị hàm f tại αk và βk và rút ngắn đoạn [ak, bk] Nếu f (αk) < f (βk) thì lấy giá trị [ak+1, bk+1] = [ak, bk] Nếu f (αk) > f (βk) thì lấy giá trị [ak+1, bk+1] = [ak, bk] Nếu f (αk) = f (βk) thì lấy giá trị [ak+1, bk+1] = [ak, bk] Giá trị tối ưu x∗ đạt được khi độ dài đoạn còn lại nhỏ hơn sai số
cho trước.
Nhận xét
a/ Trong cả hai trường hợp ta đều có:
bk+1 − ak+1 := p(bk − ak);
b/ Một trong hai điểm chia ở bước sau trùng với điểm chia ở bước trước.
Do đó thuật toán này cho phép giảm số phép tính.
Thuật toán mặt cắt vàng được mô phỏng bằng ngôn ngữ MATLAB
function mc=mat_cat_vang_1(a, b, epxilon)
L=b-a;
p=(-1 + sqrt(5))/2; count=0;
while L > epxilon
14
count=count +1
alpha=b - p*L;
beta=a + p*L;
if f(alpha) <= f(beta)
b = beta;
end;
if f(alpha) >f(beta)
a = alpha;
end;
L= b - a;
end;
mc = a;
Sau đây là kết quả chạy thuật toán giải bài toán
, x ∈ [0, 10]. f (x) = x2 + 54 x
Bảng 1.3: Nghiệm xấp xỉ tối ưu sau các bước lặp
Số bước lặp Nghiệm tối ưu Số bước lặp Nghiệm tối ưu
5 3.2204 20 2.9999
10 3.0148 25 3.0000
15 2.9998 30 3.0000
Nhận xét Thuật toán hội tụ tốt sau 30 bước lặp đạt đến nghiệm tối ưu,
tốc độ hội tụ của thuật toán mặt cắt vàng là nhanh hơn thuật toán chia đôi.
15
1.5. Phương pháp giải bài toán quy hoạch tuyến tính tổng
quát trên phần mềm MATLAB
Trong phần mềm MATLAB, bài toán Quy hoạch tuyến tính có dạng
mặc định là: Hàm mục tiêu: f = C T X → min
Các ràng buộc:
AX < b
AeqX = beq
Các giới hạn biên của nghiệm lb ≤ X ≤ ub (lb: lower bounds, ub: upper
bounds) Khi gặp bài toán f = C T X → max, chỉ cần đặt maxf = −min(−f )
Lệnh thường dùng để giải QHTT là:
[X, f val, exitf lag, output] = linprog(C, A, b, A ∈ q, beq, lb, ub)
[X, f val, exitf lag, output] = intprog(C, A, b, Aeq, beq)
Trong đó:
Lệnh linprog để lấy các nghiệm không âm
Lệnh bintprog để lấy các nghiệm nguyên có giá trị 1 hoặc 0.
Trong dấu () là các vectơ và ma trận đã cho của hàm mục tiêu và các ràng
buộc.
Trong dấu [ ] là các đại lượng cần tính:
X - giá trị tối ưu của nghiệm,
f val - giá trị min của hàm mục tiêu.
Exitf lag - số nguyên thông báo kết thúc tính toán. Các kết quả tính khi
exitf lag = 1, được coi là thành công tốt đẹp, nghĩa là hàm số hội tụ về
một nghiệm. Các kết quả tính tương ứng exitf lag ≤ 0 được coi là không
thành công với các giải thích tương ứng.
Output- cho các thông tin về phép tính đã thực hiện.
16
Chương 2
Một số thuật toán giải số bài toán
tối ưu phi tuyến không ràng buộc
Nội dung của chương 2 là trình bày một số thuật toán giải số bài toán
tối ưu phi tuyến không ràng buộc: mô hình tổng quát, điều kiện tối ưu, các
thuật toán sử dụng đạo hàm như thuật toán Gradient, thuật toán đường
dốc nhất, thuật toán Newton, thuật toán Gradient liên hợp; các thuật toán
không sử dụng đạo hàm như thuật toán tìm trực tiếp, thuật toán Powell,
thuật toán Nelder và Mead. Các kết quả được tham khảo trong các tài
2.1. Một số kiến thức cơ bản
liệu [2], [3], [6], [7], [8], [9].
2.1.1. Định nghĩa
Cho dãy {xk} hội tụ về dãy x∗, ta nói:
- Dãy {xk} hội tụ với tốc độ một cấp số nhân công bội q nếu tồn tại k0, 0 < q < 1, C < +∞ sao cho ∀k > k0: Khi đó ta có các đánh giá ||xk+1 − x∗|| ≤ q||xk − x∗|| hoặc ||xk − x∗|| ≤ qk.C - Dãy {xk} hội tụ với tốc độ trên cấp số nhân nếu ta có: ||xk+1 − x∗|| ≤ qk||xk − x∗|| hoặc ||xk − x∗|| ≤ q1.q2...qk.C với dãy ||qk|| → 0
17
khi k → +∞, C < +∞.
- Dãy {xk} hội tụ với tốc độ bậc hai nếu tồn tại C>0, tồn tại k0 sao cho: ||xk+1 − x0|| ≤ ||xk − x0||2, ∀k > k0.
2.1.2. Điều kiện tối ưu
Xét bài toán: Tìm x∗ để f (x) → min, x ∈ En trong đó x = [x1, x2, ..., xn]T
+ Điều kiện cần của tối ưu địa phương
1/ f (x) khả vi tại x∗. 2/ (cid:53)f (x∗) = 0 nghĩa là x∗ là điểm dừng
Chú ý Khi thỏa mãn điều kiện cần thì x∗ có thể là cực tiểu địa phương,
cực đại địa phương hoặc điểm yên ngựa.
+ Điều kiện đủ của cực tiểu địa phương H(x) = (cid:53)2f (x) > 0, nghĩa là ma trận Hessian xác định dương.
1 + x3
2 + 2x2
1 + 4x2
2 + 6
Ví dụ 2 Xác định các điểm dừng của hàm mục tiêu: f (x) = x3
1 + 4x1 = 0
= 3x2 ∂f ∂x1
+ Xác định điểm dừng: (cid:53)f (x) = 0
Nghiệm của hệ phương trình là các điểm: (0, 0); (0, -8/3); (-4/3, 0); (-4/3,
-8/3) + Xác định loại điểm dừng. Tính (cid:53)2f (x) :
= 0 = 6x1 + 4; = 6x2 + 8; ∂2f ∂x1∂x2 ∂2f ∂x2 1 ∂2f ∂x2 2
Ma trận Hessian:
(cid:34) (cid:35) 0 H(x) = 6x1 + 4 0 6x2 + 8
18
Các định thức:
(cid:34) (cid:35) 0 D1 = |6x1 + 4|; D2 = 6x1 + 4 0 6x2 + 8
Giá trị D1, D2 và loại điểm dừng cho ở bảng sau:
Điểm x∗ H(x) Loại điểm dừng f (x)
(0,0) D1 D2 4 32 Xác định dương Min địa phương 6
(0, -8/3) 4 -32 Không xác định Điểm yên ngựa 418/27
(-4/3, 0) -4 -32 Không xác định Điểm yên ngựa 194/27
-4 (-4/3, -8/3) 32 Xác định âm Max địa phương 50/3
Nhận xét Trong trường hợp tổng quát khi hàm mục tiêu là phi tuyến
thì việc xác định các điểm dừng là không thực hiện được. Do đó chúng ta
2.2. Các thuật toán sử dụng đạo hàm
phải nghiên cứu các thuật toán xác định nghiệm xấp xỉ.
2.2.1. Thuật toán Gradient
Đây là phương pháp phổ biến nhất, luôn luôn hội tụ.
Thuật toán sử dụng sơ đồ lặp
x(k+1) = x(k) − λk (cid:53) f (x(k))
Trong đó λk là hệ số xác định độ dài của bước đi theo hướng Gradient. Ta có thể chọn bước đi bằng hằng số cho cả quá trình hoặc tính giá trị bước
đi tối ưu theo phương pháp tối ưu một tham số. Điều kiện dừng lặp ||x(k+1) − x(k)|| ≤ ε.
Quy tắc Armijo xác định bước đi tối ưu tổng quát
Bước 1: Chọn λ tùy ý như nhau với mọi bước lặp (chẳng hạn λk = 1), Xác định điểm x = x(k) − λ (cid:53) f (x(k)) là một điểm trên hướng giảm. Bước 2: Tính f (x) = f (x(k) − λ (cid:53) f (x(k)))
19
Bước 3: Kiểm tra
f (x(k)) − f (x) ≥ ελ|| (cid:53) f (x(k))||2
Thì lấy λk = λ, ngược lại lấy λk = αλ, α ∈ (0, 1) Sử dụng quy tắc Armijo, thuật toán Gardient được mô tả bằng ngôn ngữ
MATLAB như sau:
function gd=gradient_amo(epxilon, k)
clc;
X=[-1; -2]; Xluu=[1; 1]; count =0; ss=10; ep=10∧(-4); alpha=1/2;
while and(ss > epxilon, count < k)
count=count + 1;
grad=[dh1x(X); dh1y(X)];
dk=grad;
%Quytac Armijo;
lamda=1;
while f(X)-f(X- lamda*grad) < ep*lamda * (grad.*dk)
lamda=lamda*alpha;
end;
lamdak = lamda;
X= X - lamdak*dk;
ss=chuan1(X - Xluu, 2);
Xluu=X;
end;
count
gd=X
Các kết quả thực hiện thuật toán được cho trong bảng dữ liệu sau đây.
20
Bảng 2.2: Nghiệm xấp xỉ tối ưu sau các bước lặp
f (x, y) = x2 + y2 + x − y
Số bước lặp Sai số Nghiệm xấp xỉ
1 [1; 1]
2 1.5 10×e − 10 [-0.5; 0.5]
Nhận xét Thuật toán hội tụ tốt sau hai bước lặp đạt đến nghiệm
tối ưu.
2.2.2. Thuật toán đường dốc nhất
Cơ sở phương pháp
Phương pháp đường dốc nhất (The steepest descent method) là một
trong các phương pháp cổ điển thông dụng nhất giải bài toán quy hoạch
phi tuyến không ràng buộc nhiều biến.
Xét bài toán quy hoạch phi tuyến không ràng buộc tổng quát:
Minf (x), x = (x1, x2, ..., xn) ∈ Rn
Ta gọi vectơ d là hướng giảm của hàm f : Rn → R tại x nếu ∃δ > 0 sao
cho f (x + λd) < f (x), ∀x ∈ (0, δ).
Giả sử hàm f khả vi tại x. Ngoài ra giả sử rằng (cid:53)f (x) (cid:54)= 0. Lúc đó, có thể chứng minh được ¯d = − (cid:53) f (x)/|| (cid:53) f (x)|| là hướng giảm nhanh nhất, tức là ¯d là lời giải của bài toán Minf (cid:48)(x, d), trong đó f (cid:48)(x, d) là đạo hàm theo hướng d tại x, với điều kiện ||d|| ≤ 1.
Thật vậy, do f khả vi tại x nên:
f (x + λd) = f (x) + (cid:53)f (x)T (λd) + ||λd||α(x, λd) (2.1)
(cid:48)
Với α(x, λd) = 0. Vậy đạo hàm theo hướng d tại x chính là lim λ→0−
= (cid:53)f (x)T d. f (x, d) = lim λ→0 f (x + λd) − f (x) λ
21
Do ||d|| ≤ 1 nên theo bất đẳng thức schwartz ta có
(cid:53)f (x)T d ≥ −|| (cid:53) f (x)|| ||d|| ≥ −|| (cid:53) f (x)||.
Với ¯d = − (cid:53) f (x)/|| (cid:53) f (x)|| ta có (cid:53)f (x)T = −|| (cid:53) f (x)||, nên ¯d là hướng giảm nhanh nhất của hàm f tại x. Nếu biểu thức ||λd||α(x, λd) được coi là
bằng 0 trong công thức (2.1) thì với một giá trị λ > 0 cố định và với điều kiện ||d|| ≤ 1, f (x + λd) đạt giá trị cực tiểu tại ¯d = − (cid:53) f (x)/|| (cid:53) f (x)||. Tuy nhiên, biểu thức ||λd||α(x, λd) không nhất thiết phải bằng 0, nên sau khi hướng giảm nhanh nhất ¯d đã được chọn, cần xác định λ ≥ 0 để cực tiểu hóa f (x + λ ¯d).
Sau đây là thuật toán của phương pháp đường dốc nhất. Dựa trên lý
thuyết về ánh xạ co, có thể chứng minh thuật toán này hội tụ tới điểm ¯x có (cid:53)f (¯x) = 0 với điều kiện dãy {xk} được phát sinh trong thuật toán đều
nằm trong một tập giới nội. Nếu hàm f (x) là hàm lồi thì ¯x sẽ là phương
án tối ưu toàn cục của bài toán quy hoạch phi tuyến không ràng buộc
đã cho.
Thuật toán đường dốc nhất
Bước khởi tạo: Cho ε > 0 làm sai số kết thúc. Lấy một điểm xuất phát
x1, đặt k:=1 và chuyển sang các bước lặp.
Các bước lặp (bước lặp thứ k) Bước 1: Nếu || (cid:53) f (xk)|| > ε thì đặt dk = − (cid:53) f (xk) và chuyển sang
Bước 2.
Bước 2: Tìm λk là phương án tối ưu của bài toán cực tiểu hóa hàm một biến f (xk + λdk) (phụ thuộc vào λ ≥ 0). Đặt xk+1 = xk + λkdk, k := k + 1
và chuyển về Bước 1.
Bước kết thúc: Nếu || (cid:53) f (xk)|| ≤ ε thì dừng.
Thuật toán có thể dùng quy tắc Armijo để chọn bước đi tối ưu. Thuật
toán đường dốc nhất được mô tả bằng ngôn ngữ MATLAB như sau:
function dg=dg(epxilon, k)
22
clc;
X=[1; 3]; ss=10; ep=1/2; lamda =1/3; count =0;
while and(ss > epxilon, count < k)
count=count + 1;
grad=[dfx(X); dfy(X)];
dk=grad;
%Quytac Armijo;
alpha=1;
while f(X+ alpha*dk) - f(X) > ep*alpha * grad.*dk
alpha=alpha*lamda;
grad=[dfx(X); dfy(X)];
dk = -grad;
end;
lamdak=alpha;
X= X+lamdak*dk;
ss=chuan1(grad, 2);
end;
count
dg=X;
Các kết quả thực hiện thuật toán được cho trong bảng dữ liệu sau đây
Bảng 2.3: Nghiệm xấp xỉ tối ưu sau các bước lặp
f (x, y) = x2 + y2 + x − y
23
Số bước lặp Sai số Nghiệm xấp xỉ
1 5 [0; 1.3]
2 1 [-0.3; 0.7]
3 0.5 [-0.4; 0.5]
4 0.1 [-0.48; 0.53]
5 0.06 [-0.49; 0.51]
6 0.02 [-0.497; 0.50]
7 [-0.4999; 0.501]
8 [-0.4999; 0.5001]
9 0.001 10×e − 5 10×e − 7 [-0.5; 0.5]
Nhận xét Thuật toán hội tụ tốt sau 9 bước lặp đạt đến nghiệm tối ưu.
2.2.3. Thuật toán Newton
Trong phương pháp đường dốc nhất, quy tắc dịch chuyển cho bởi xk+1 = xk + λkdk với dk = − (cid:53) f (xk). Trong phương pháp Newton, ta cũng có quy tắc dịch chuyển tương tự với λk được thay thế bởi H −1(xk), trong đó H(xk) là ma trận Hessian được tính tại điểm xk với điều kiện ma trận này khả nghịch. Giả sử rằng dãy {xk} hội tụ tới ¯x với (cid:53)f (¯x) = 0 và H(¯x) xác định dương, trong đó f (x) là hàm khả vi cấp hai. Lúc đó với các điểm xk khả sát ¯x. H(xk) cùng xác định dương nên là ma trận khả nghịch.
Sau đây, chúng ta giải thích ý nghĩa của quy tắc dịch chuyển:
xk+1 = xk − H −1(xk). (cid:53) f (xk) trong phương pháp Newton. Đối với hàm
khả vi cấp hai chúng ta có thể viết:
f (x) = f (xk)+(cid:53)f (xk)T (x−xk)+ (x−xk)T H(xk)(x−xk)+||x−xk||2α(x, x−xk) 1 2
trong đó, α(xk, x − xk) = 0. Bởi vậy, có thể xấp xỉ f (x) bởi: lim x→xk
q(x) = f (xk) + (cid:53)f (xk)T (x − xk) + (x − xk)T H(xk)(x − xk) ≈ f (x) 1 2
24
Ngoài ra, dễ thấy điều kiện cần để q(x) đạt giá trị cực tiểu là:
(cid:53)q(x) = 0 ⇔ (cid:53)f (xk) + H(xk)(x − xk) = 0
Giả sử ma trận H(xk) khả nghịch thì điểm tiếp theo nên xem xét chính là điểm xk+1 = xk − H −1(xk) (cid:53) f (xk).
Có thể chứng minh được phương pháp Newton hội tụ (khá nhanh) với điều kiện điểm xuất phát x1 nằm sát gần ¯x với (cid:53)f (¯x) = 0 và ma trận
H(¯x) là khả nghịch. Để khắc phục điều kiện ngặt nghèo này, phương pháp
Newton cải biên đã được đề xuất. Tuy nhiên đây là giải thuật phức tạp.
xkf (xk)
(cid:53)f (xk)
H(xk)
H(xk)−1
−H(xk)−1 (cid:53) f (xk)
xk+1
Lặp k
(0; 3)
50 −4
4
(-44; 24)
(0,67; -2,67)
(0,67; 0,33)
1
1 384
52
−4
8
8 4
50
(0, 67; 0, 33)
23, 23 −4
4
2
(-9,39;-0,04)
(0,44; 0,23)
(1,11; 0,56)
1 169,84
3, 33
−4
8
8 4
23, 23
11, 5 −4
4
(1, 11; 0, 56)
3
(-0,84;-0,04)
(0,3; 0,14)
(1,41;0,7)
1 76
−4
8
11, 5
0, 63
6, 18 −4
4
(1, 41; 0, 7)
4
(-0,8;-0,04)
(0,2; 0,1)
(1,61; 0,80)
1 33,4
−4
8
6, 18
0, 12
(1, 61; 0, 8)
3, 83 −4
4
5
(-0,22;-0,04)
(0,13; 0,07)
(1,74; 0,87)
1 16,64
0, 02
−4
8
8 4 8 4 8 4
3, 88
(1, 74; 0, 87)
2, 81 −4
4
6
(-0,07; 0)
(0,09; 0,04)
(1,83; 0,91)
1 6,48
0, 05
−4
8
8 4
2, 81
(1, 83; 0, 91)
7
(0,0003;-0,04)
0, 0009
Ví dụ 3 Giải bài toán Minf (x) = (x1 −2)4 +(x1 −2x2)2 bằng phương pháp Newton. Quá trình giải được minh họa trong bảng tóm tắt:
Thuật toán có thể dùng quy tắc Armijo để chọn bước đi tối ưu. Thuật
25
toán Newton được mô tả bằng ngôn ngữ MATLAB như sau:
function nw=newton_amo(epxilon, k)
format short e
clc;
X=[1; 3.5]; a=0; b=1;
Xluu=[1; 1]; count =0; ss=10; ep=1/2;
while and(ss > epxilon, count < k)
count=count + 1;
grad=[dh1x(X); dh1y(X)];
H=[dh2x(X) dh2xy(X); dh2xy(X) dh2y(X)];
%Quytac Armijo;
alpha=1;
dk=-inv(H)*grad;
while f(X+ alpha*dk)-f(X) > ep*alpha*grad.*dk
alpha=alpha*lamda;
dk=-inv(H)*grad;
end;
lamdak = alpha;
X= X + lamdak*dk;
ss=chuan11(X - Xluu, 2);
Xluu=X;
end;
count
X
ss
26
2.3. Các thuật toán không sử dụng đạo hàm
2.3.1. Phương pháp tìm trực tiếp (Direct search)
Xét bài toán: Tìm x∗ sao cho f (x) → max; x ∈ Rn
Tư tưởng chính của thuật toán: ở mỗi bước của thuật toán chỉ biến đổi
một thành phần xi của x còn các thành phần khác giữ nguyên và cứ làm như vậy cho tới khi nhận được cực đại.
n ]T và số gia của biến số (cid:52)x((cid:52)x1, ..., (cid:52)xn)
1 , ..., x(0)
Các bước thuật toán:
Bước 1: Thăm dò bước 1 * Cho các giá trị của x(0) = [x(0) * Tính giá trị f (x(0)).
2 , ..., x(0) n )
1 , x(0)
1 + (cid:52)x1 → tính α − f (x(0)
1 = x(0) x(1)
t = x(0)
t − (cid:52)x(0)
1 → Tính α
1 + (cid:52)x1 và x(0)
1 − (cid:52)x1)
2 = x(0)
2 ± (cid:52)x2 ... cho đến khi tất cả các biến
* Cho một thành phần của x biến đổi, các thành phần khác giữ nguyên:
+ Nếu α ≤ f (x(0)) thì lấy: x(1) + Nếu f (x) không cải tiến được cả hai phía (với x(0) thì cố định x(0) 1 không cho biến đổi nữa. + Tiếp tục tiến hành với x(1) đều được biến đổi.
Như vậy trong Bước 1, tại mỗi bước dịch chuyển theo một biến độc lập,
giá trị hàm mục tiêu được so sánh với giá trị của nó tại điểm trước. Nếu
hàm mục tiêu được cải tiến tại một bước nào đó thì giá trị cũ được thay
bằng giá trị mới trong những so sánh sau đó. Nếu hàm mục tiêu được cải
tiến thì giữ nguyên giá trị cũ.
Bước 2: Tìm theo mẫu (Pattern search)
Khi kết thúc Bước 1, ta xác định được x(k). Ở bước tiếp theo ta lấy
theo mẫu.
Tức là:
x(k+1) = tx(k) − x(b)
27
i − x(b)
i
x(k+1) = tx(k)
Trong đó: x(b) là điểm cơ sở (base point), ở lần lặp đầu x(b) = x(0)
t là số biến thăm dò cần thiết. Đối với f (x1, x2) ta có t = 2. Bước 3:
* Xuất phát từ điểm x(k+1) tính f (x(k+1)) và so với giá trị f (x) ở bước
tìm theo mẫu để xem Bước 2 có kết quả không.
* Nếu Bước 2 có kết quả tại x(k+2) thì tìm theo mẫu được coi là có kết quả nếu: f (x(k+2)) ≥ f (x(k)) và quá trình được lặp lại với điểm xuất phát x(b) = x(k).
* Nếu Bước 2 không kết quả thì kết luận rằng tìm theo mẫu đã thất
bại và lại phải thực hiện lại Bước 1 sao cho việc xác định hướng mới có
hiệu quả.
* Nếu Bước 1 liên tiếp không cho hướng mới hiệu quả thì liên tiếp giảm
(cid:52)x cho tới khi hoặc xác định được hướng mới có hiệu quả, hoặc (cid:52)x đã
nhỏ hơn sai số cho phép.
Việc không có khả năng tăng f (x) khi (cid:52)x đã khá bé tức là tối ưu cục bộ
đã đạt được.
Nhận xét Phương pháp trực tiếp nói chung chậm hội tụ hơn so với
phương pháp dùng đạo hàm nhưng sử dụng lại tiện hơn vì không phải
tính đạo hàm.
Thuật toán có thể cài đặt được bằng ngôn ngữ lập trình MATLAB.
2.3.2. Phương pháp Powell
1 , ..., s(k) n i ) → min để xác
i+1 + λis(k)
Bước 1: Chọn điểm xuất phát x(k) 0
(điểm dừng x0, bước lặp k) Khi k = 0, chọn các hướng song song với các trục tọa độ của Rn: s(k) Dùng phương pháp tối ưu một tham số: f (x(k) định λi.
28
i
Khi đó xác định được λ(k) và do đó xác định được:
i+1 = x(k) x(k)
i + λ(k)
i s(k)
i
n − x(k) 0
n+1 = 2x(k) Bước 2: Chọn điểm tiếp theo: x(k) Tiếp tục tìm theo hướng Gradient liên hợp.
; i = 1, 2, ..., n
2.3.3. Phương pháp Nelder và Mead
J.A Nelder và R.Mead đã dùng đơn hình là những đa diện biến dạng
(flexible polyhedron) nhờ ba phép biến đổi: ánh xạ gương, phép co và
phép dãn.
Thuật toán được thực hiện theo các bước: Bước 1: Trong không gian Rn, chọn một đơn hình có đỉnh tại các điểm
(cid:80) xi; i (cid:54)= 1. Tính f (x0).
(x1, x2, ..., xn+1). Tính giá trị các hàm mục tiêu f (x) ở các đỉnh của đơn hình fi = f1(x), ..., fn+1(x) = f (xn+1). Bước 2: Xác định các giá trị của hàm f : giá trị lớn nhất f1, giá trị tiếp theo fg, giá trị nhỏ nhất fb; các đỉnh tương ứng là: x1, x2, xb. Bước 3: Xác định trọng tâm của hình (trừ đỉnh x1) theo công thức: x0 = 1 n Bước 4: Thực hiện ánh xạ gương x1 qua x0 nhận được xr. Tính f (xr) = fr. Nếu hệ số ánh xạ α > 0 thì vị trí xr được xác định theo công thức:
xr − x0 = α(x0 − x1)
xr = (1 + α)x0 − αx1
α = |xr − x0|/|x0 − x1|
Bước 5: So sánh giá trị fr và fb
1/ Nếu fr < fb thì ta nhận được giá trị nhỏ nhất của hàm. Hướng từ x0 đến xr là hướng tốt nhất cho dịch chuyển, ta kéo dãn theo hướng đó và
29
tìm điểm xp. Hệ số dãn γ > 1 được tìm từ biểu thức:
xp − x0 = γ(xr − x0)
xp = γxr + (1 − γ)x0
γ = |xp − x0|/|xr − x0|
Sau đó tính f (xp) = fp
a) Nếu fp < fb thì thay x1 bởi điểm xp và kiểm tra điểm thứ (n+1) của đơn hình về tính hội tụ đến cực tiểu. Nếu hội tụ thì ngừng tính, nếu chưa
thì quay lại Bước 2.
b) Nếu fp > fb thì bỏ điểm xp vì đã đi quá xa. Thay x1 bởi xr, kiểm tra
sự hội tụ. Nếu vẫn chưa hội tụ thì quay về Bước 2.
2/ Nếu fr > fb nhưng fr < fg thì xr là điểm tốt hơn so với hai điểm khác của đơn hình. Thay x1 bằng xr. Kiểm tra nếu chưa hội tụ thì quay về Bước 2.
3) Nếu fr > fb nhưng fr > fg thì chuyển sang Bước 6.
Bước 6: So sánh các giá trị fr và f1
1/ Nếu fr > f1 thì chuyển về phép co (Bước 6, đoạn 2): Nếu fr < f1 thì thay x1 bởi điểm xr và thay f1 bởi fr. Nếu fr > fg chuyển sang phép co. 2/ Thực hiện phép co:
- Trường hợp fr > f1 ta đã dịch chuyển quá xa theo xu hướng từ x1 đến x0 do đó cần sửa lại nhờ phép co để tìm x0 như sau:
xc − x0 = β(x1 − x0) với 0 < β < 1; β : hệ số co
xc = βx1 + (1 − β)x0
- Trường hợp fr < f1 ta thay x1 bằng xr và sau đó thực hiện phép co để tìm xc như sau:
xc − x0 = β(xr − x0)
xc = βxr + (1 − β)x0
30
Bước 7: So sánh các giá trị fc và f1. Tính fc
1/ Nếu fc
đạt được cần sang Bước 8.
Bước 8: Thu nhỏ kích thước của đơn hình tại một nửa, lấy xb làm chuẩn. Như vậy điểm xi được thay bằng điểm:
(xi − xb) = (xi + xb) xi − 1 2 1 2
Tính các giá trị fi(i = 1, 2, ..., n + 1). Kiểm tra sự hội tụ nếu chưa hội tụ thì quay lại Bước 2.
n+1 (cid:88)
Bước 9: Tính
i=1
σ2 = (fi − f )2/(n + 1); f = fi/(n + 1)
Nếu σ nhỏ hơn độ chính xác (cid:15) thì giá trị các hàm rất gần nhau và do đó
1 + x2
2 − 40x1 − 12x2 + 136 = 4(x1 − 5)2 + (x2 − 6)2 → min
điểm cực tiểu gần là xb. Ví dụ 4 f (x) = 4x2
1 = [8; 9]T ; x(0) x(0)
2 = [10; 11]T ; x(0)
3 = [8; 11]T .
* Chọn đơn hình là tam giác với ba đỉnh ban dầu là:
* Bước lặp k=0
f (8; 9) = 45; f (10; 11) = 12; f (8; 11) = 65
3
Xác định trọng tâm của x(0)
1 và x(0) 1 2
[(8 + 10 + 8) − 10] = 8
[(9 + 11 + 11) − 11] = 10 x(0) 4.2 = x(0) 4.1 = 1 2
4 = [8; 10]T
Ta nhận được x(0)
31
Ánh xạ qua x(0) 4 :
x(0) 5.1 = 8 + 1(8 − 10) = 6
5 = [6; 9]T và f(6; 9)=13. Vì f(6; 9)=13 < f(8; 9) =45 ta thực
x(0) 5.2 = 10 + 1(10 − 11) = 9
Ta được x(0) hiện phép dãn:
x(0) 6.1 = 8 + 2(6 − 8) = 4
6 = [4; 8]T và f(4; 8)=8.
x(0) 6.2 = 10 + 2(9 − 10) = 8
Ta nhận được x(0) Cuối cùng ta nhận được x∗ = [5; 6]T và f (x∗) = 0.
32
Chương 3
Một số thuật toán giải số bài toán
tối ưu phi tuyến có ràng buộc
Nội dung chương 3 là tìm hiểu một số thuật toán giải số bài toán tối
ưu phi tuyến có ràng buộc như khái niệm hàm Lagrange, phương pháp
hàm phạt, các thuật toán tìm nghiệm xấp xỉ. Các kết quả được tham khảo
3.1. Một số kiến thức cơ bản
trong các tài liệu [2], [3], [9].
3.1.1. Hàm Lagrange
Xét bài toán quy hoạch phi tuyến có ràng buộc tổng quát:
(3.1) Min(Max)f (x), với x ∈ D = {x ∈ Rn : gi(x) ≤ 0, ∀i = 1, m}.
Lúc đó, hàm Lagrange tương ứng với bài toán trên có dạng sau:
F (x, λ) = f (x) + λ1g1(x) + ... + λmgm(x), với điều kiện λi ≥ 0, ∀i = 1, m (các số λi ≥ 0, ∀i = 1, m được gọi là các nhân tử Lagrange).
33
Ký hiệu
i , hàm Lagrange được định nghĩa trên đây được viết lại dưới
λ = và G(x) = thì F (x, λ) = f (x) + λT G(x) λ1 λ2 ... g1(x) g2(x) ... λm gm(x)
m).
2, ..., s2
1, s2
m (cid:80) i=1
i gi(x), với s2 = (s2 s2 Chúng ta gọi các điểm (x, λ) = (x, s2) là điểm dừng của hàm Lagrange nếu điểm (x, s) ∈ Rn+m thỏa mãn hệ điều kiện sau đây:
Đặt λi = s2 dạng F (x, s2) = f (x) +
∂f ∂xj
∂gi(x) ∂xj
∂F ∂xj
m (cid:80) i=1
∂F ∂si
+ , ∀j = 1, n = 0, ∀j = 1, n s2 i ⇔ (3.1) = 0, ∀i = 1, m sigi(x) = 0, ∀i = 1, m
3.1.2. Thiết lập điều kiện tối ưu Kuhn - Tucker
Trong mục này, với mục đích với tìm hiểu bước đầu, chúng ta sẽ nghiên
cứu cách thiết lập điều kiện tối ưu Kuhn - Tucker đối với bài toán quy
hoạch phi tuyến có ràng buộc và xem xét nó qua một số ví dụ cụ thể mà
không đi sâu vào việc chứng minh các điều kiện này một cách chặt chẽ.
Có thể nói rằng điều kiện Kuhn - Tucker là điều kiện cơ bản nhất trong
lý thuyết tối ưu phi tuyến và là nhiều cơ sở cho các phương pháp tối ưu
phi tuyến cổ điển.
Xét hệ điều kiện bao gồm điều kiện dừng của hàm Lagrange và điều
kiện ràng buộc của Bài toán quy hoạch phi tuyến có ràng buộc (3.1)
∂f ∂xj
∂gi(x) ∂xj
m (cid:80) i=1
m (cid:80) i=1 λT G(x) = 0
= 0, ∀i = 1, n + (cid:53)f (x) + λi (cid:53) gi(x) = 0 λi
λigi(x), ∀i = 1, m ⇔
G(x) ≤ 0 gi(x) ≤ 0, ∀i = 1, m
λ ≥ 0 λi ≥ 0, ∀i = 1, m
34
Hệ điều kiện trên đây được gọi là điều kiện Kuhn - Tucker của Bài toán
quy hoạch phi tuyến có ràng buộc (3.1)
Ví dụ 5
Thiết lập điều kiện Kuhn - Tucker cho Bài toán quy hoạch phi tuyến có
ràng buộc sau:
Minf (x) = (x1 + 1)2 + (x2 − 1)2
Với điều kiện x = (x1, x2) ∈ D là miền ràng buộc được xác định bởi
g1(x) = x1 − 2 ≤ 0 x1 − 2 ≤ 0 g2(x) = x2 − 1 ≤ 0 ⇔ x2 − 1 ≤ 0 g3 = −x1 ≤ 0 x1, x2 ≥ 0 g4 = −x1 ≥ 0.
Có thể kiểm nghiệm được rằng trong ví dụ này chúng ta có Bài toán quy
hoạch phi tuyến có ràng buộc với hàm Lagrange:
F (x, λ) = (x1+1)2+λ1(x1−2)+λ2(x2−1)−λ3x1−λ4x2, (λi ≥ 0, ∀i = 1, 4).
Điều kiện Kuhn - Tucker của bài toán được viết như sau: (3.2)
(3.3)
(3.4)
(3.5)
(3.6)
(3.7)
(3.8)
(3.9)
(3.10)
(3.11)
(3.12) 2(x1 + 1) + λ1 − λ3 = 0 2(x2 − 1) + λ2 − λ4 = 0 λ1(x1 − 2) = 0 λ2(x2 − 1) = 0 λ3(−x1) = 0 λ4(−x2) = 0 x1 − 2 ≤ 0 x2 − 1 ≤ 0 −x1 ≤ 0 −x2 ≤ 0 λ1, λ2, λ3, λ4 ≥ 0
35
Từ (3.2) và (3.6) suy ra: x1[2(x1 + 1) + λ1] = 0 ⇒ x1 = 0 ⇒ theo (3.4) ta có λ1 = 0 Từ (3.3) và (3.7) suy ra:
(cid:34)
x2[2(x1 − 1) + λ2] = 0 ⇒ x2 = 0 ⇒ λ2 = 0 2(x2 − 1) + λ2 = 0 ⇒ x2 = 1, λ2 = 0
Với điều kiện λ1 = λ2 = 0, ta thấy trong hai điểm (x1 = 0, x2 = 0) và (x1 = 0, x2 = 1) chỉ có điểm (x1 = 0, x2 = 1) với (λ1 = λ2 = 0, λ3 = 3, λ4 = 0) thỏa mãn điều kiện dừng của hàm Lagrange.
Vậy phương án tối ưu toàn cục là (x1 = 0, x2 = 1) với fmin = 1.
Nhận xét Bằng cách đưa vào hàm Lagrange, bài toán quy hoạch phi
tuyến có ràng buộc luôn được đưa về bài toán không ràng buộc bằng cách
thay hàm mục tiêu bằng hàm Lagrange và thêm vào các biến là các nhân
tử. Từ đó chúng ta luôn luôn có thể sử dụng các thuật toán trong chương
2 để giải bài toán này. Tuy nhiên chúng ta có thể xây dựng các thuật toán
giải trực tiếp không cần sử dụng hàm Lagrange. Sau đây chúng ta sẽ xét
3.2. Một số thuật toán
một số thuật toán.
3.2.1. Thuật toán Gradient
Xét bài toán
f (X) → min
gi(X) ≤ 0, i = 1, 2, ..., m
Thuật toán Gradient
Định nghĩa: Nón các hướng chấp nhận được tại a ∈ D là tập:
P (a) = {X| (cid:104)(cid:53)gi(a), (x − a)(cid:105) ≤ −gi(a); i = 1, 2, ..., m}
36
Hướng chấp nhận tại a: z = (x−a)+ε(x(0) −a); ε > 0. x(0) điểm xuất phát. Bước 1: Tìm điểm xuất phát x(0) ∈ D
Bước 2: ∀k = 0, 1, ...
+ Xây dựng nón P (x(k)) + Giải bài toán quy hoạch tuyến tính phụ (cid:10)(cid:53)f (x(k)), (x − x(k))(cid:11) →
min với x ∈ P (x(k))
+ Gọi lời giải là x(∗)
Bước 3: Kiểm tra điều kiện tối ưu: Đặt σk = (cid:10)(cid:53)f (x(k)), (x∗ − x(k))(cid:11)
+ Nếu σk ≥ 0 ⇒ Xopt = x(∗) ⇒ Dừng + Nếu σk < 0 chưa tối ưu. Khi đó hướng chấp nhận được là:
z(k) = (x(∗) − x(k)) + εk(x(0) − x(k))
(cid:111) (cid:110) 0 < σ < 1 trong đó 1, σ σk |δk|
Để thuật toán hội tụ, ta sẽ chọn εk = min ta thường lấy (cid:69) (cid:68) (cid:53)f (x(k)), (x(0) − x(k)) σ = δk = 1 2
Bước 4: Xác định phương án tiếp sau:
(cid:48)
+ Xác định λ(k) lớn nhất: Đặt ϕ(λ) = f (x(k) + λz(k))
λ(k) = max{λ > 0, x(k) + λz(k) ∈ D, ϕ (λ) ≤ 0}
+ Khi đó x(k+1) = x(k) + λ(k)z(k)
Chú ý λ(k) thỏa mãn f (x(k) + λz(k)) → min|λ > 0, x(k) + λz(k) ∈ D
Nhận xét
+ Phương pháp Gardient đã đưa việc giải bài toán phi tuyến về một
dãy các bài toán QHTT phụ. Các bài toán này dễ dàng tìm được phương
án tối ưu trên phần mềm MATLAB bằng thủ tục linprog(...).
+ Việc xác định bước đi tối ưu luôn thực hiện được bằng các thuật toán
chia đôi hoặc lát cắt vàng vì f là hàm lồi.
Thuật toán có thể mô phỏng bằng các ngôn ngữ lập trình nếu chúng ta
mô tả được thủ tục xác định nón P (a) tùy theo từng bài toán cụ thể.
37
Ví dụ 6 Hàm mục tiêu
1 + x2
2 + 18 → min
f (x) = −6x1 − 4x2 + x2
Các ràng buộc:
2 < 6
4x1 + 3x2 ≤ 24
xj ≥ 0, j = 1, 2
1 − 5x1 + x2 x2 Tìm x∗ sao cho thỏa mãn mục tiêu và các ràng buộc.
Giải:
Lấy một điểm bất kỳ trong miền ràng buộc: x(0) = (1; 5; 3)
Nghiệm lại: 4 x 1,5 + 3 x 3 =15 <24
(1, 5)2 − 5.1, 5 + 32 = 3, 75 < 6
Tại điểm này f (x(0)) = 8, 25
Đối Gradient − (cid:53) f (x):
− = 6 − 2.1, 5 = 3 = 6 − 2x1 → −
= 4 − 2.3 = −2 − = 4 − 2x2 → − ∂f ∂x1 ∂f ∂x2 ∂f (x(0)) ∂x1 ∂f (x(0)) ∂x2
− (cid:53)f (x(0)) = (3; −2)
− (cid:53)f (x(0)) (cid:54)= 0 nên x(0) không phải điểm cực trị.
Lấy x(1) = x(0) − λ (cid:53) f (x(0)) = (1, 5 + 3λ; 3 − 2λ)
Tính bước đi λ:
= [− (cid:53) f (x(1))].[− (cid:53) f (x(0))] = 0
d (cid:52) f dλ = [6 − 2(1, 5 + 3λ); 4 − 2(3 − 2λ)].[3; −2] = 0
dλ2 = −26 < 0 nên λ = 0, 5 đại lượng (cid:52)f đạt giá trị lớn nhất.
13 − 26λ = 0 ⇒ λ = 0, 5
Vì d2(cid:52)f Chọn điểm x(1) = (1, 5 + 3.0, 5; 3 − 2.0, 5) = (3; 2)
38
Thử lại xem x(1) có trong miền ràng buộc không?
4.3 + 3.2 = 18 < 24
32 − 5.3 + 22 = −2 < 6
Tính − (cid:53) f (x(1)) = (6 − 2.3; 4 − 2.2) = (0; 0) Nghĩa là không có cách di chuyển nào từ x(1) làm giảm f nữa, vậy x(1) = x(∗)
x(∗) = [3; 2]T ; f (x(∗)) = 5
3.2.2. Phương pháp hàm phạt (Penalty function method)
Bài toán cơ bản của quy hoạch phi tuyến tính bị ràng buộc tổng quát
có dạng:
f (x) → min
gj(x) ≤ 0; j = 1, m
Tìm x(∗) sao cho mục tiêu và các ràng buộc được thỏa mãn.
Tư tưởng chính:
Mục đích của phương pháp hàm phạt là đưa việc giải bài toán bị ràng
buộc về việc giải một chuỗi các bài toán không bị ràng buộc.
Trong phương pháp hàm phạt, ta thay thế hàm mục tiêu ban đầu f (x)
m (cid:88)
bởi hàm mục tiêu mở rộng chứa thông số dương rk; k = 1, 2, 3, ...
j=1
Pk = P (x, rk) = f (x) + rk Gj[gj(x)]
Hàm Pk được gọi là hàm phạt, tham số rk dương được gọi là tham số phạt. Hàm Gj của ràng buộc gj(x) được xây dựng sao cho:
- Bên ngoài tập chấp nhận được thì P (x, rk) khác với f (x): - Khi giải liên tiếp các bài toán không ràng buộc
P (x, rk) → min; k = 1, 2, ... thì Gj → 0, do đó P (x∗, r) → f (x∗)
39
Các phương pháp hàm phạt được chia thành hai nhóm: phương pháp
hàm phạt trong, phương pháp hàm phạt ngoài.
a) Phương pháp hàm phạt trong (interior penalty function method),
hàm Gj thường được chọn dưới dạng:
Gj = − 1 gj(x)
Gj = log[−gj(x)]
Cực tiểu hàm Pk nằm trong miền nghiệm chấp nhận được và hội tụ về
nghiệm của bài toán cơ bản f (x) khi thông số rk biến đổi đều đặn.
b) Phương pháp hàm phạt ngoài (exterior penalty function method),
hàm GJ thường được chọn dưới dạng:
Gj = max[0, gj(x)] hoặc Gj = {max[0, −gj(x)]}2
3(x1 + 1)3 + x2 → min
Khi đó cực tiểu của hàm Pk nằm ngoài miền nghiệm chấp nhận được và hội tụ về nghiệm của bài toán cơ bản f (x) khi thông số rk biến đổi đều. Ví dụ 7 Hàm mục tiêu: f (x) = 1
Ràng buộc:
g1(x) = −x1 + 1 ≤ 0
g2(x) = −x2 ≤ 0
Giải:
Theo phương pháp hàm phạt trong, hàm mục tiêu mở rộng là:
(cid:18) (cid:19) − → min P (x, r) = (x1 + 1)3 + x2 − r 1 3 1 −x1 + 1 1 x2
Để giải bài toán không ràng buộc của P (x, r) ta dùng điều kiện cần:
1 − 1)2 = r
r = (x1 + 1)2 − ∂P ∂x1 (1 − x1)2 = 0 → (x2
2 = r
= 1 − → x2 ∂P ∂x2 r x2 2
40
1
Các phương trình này cho:
1 2
2 + 1
2 ; x∗
2(r) = r
(cid:16) (cid:17) 1 r x∗ 1(r) =
1
2 + 1
2 + 1
2
1 r −
3 2
(cid:21)3 (cid:20)(cid:16) (cid:17) 1 r − Pmin(r) = (cid:17) 1 1 3 + 1 r2 1 (cid:16) 1 r
Để nhận được nghiệm của bài toán gốc, ta thấy:
Pmin(r) x∗ 1(r)
x∗ 2(r) fmin = lim r→0 x∗ 1 = lim r→0 x∗ 2 = lim r→0
1, x∗
2 được cho trong bảng sau:
Các giá trị của f, x∗
Bảng 3.1:
Giá trị r f (r)
1000 x∗ 1(r) 5,71164 x∗ 2(r) 31,62278 Pmin(r) 376,2636 132,4003
100 3,31662 10,00000 89,9772 36,8109
10 2,04017 3,16228 25,3048 12,5286
1 1,41421 1,00000 9,1046 5,6904
0,1 1,14727 0,31623 4,6117 3,6164
0,01 1,04881 0,10000 3,2716 2,9667
0,001 1,01569 0,03162 2,8569 2,7615
0,0001 1,00499 0,01000 2,7267 2,6967
0,00001 1,00158 0,00316 2,6856 2,6762
0,000001 1,00050 0,00100 2,6727 2,6697
0 8/3 8/3 Nghiệm chính xác 1
Sau đây chúng ta sẽ nghiên cứu chi tiết một số thuật toán sử dụng tư
tưởng hàm phạt
A. Một số phương pháp hàm phạt điểm trong
Hàm phạt p(x) thỏa mãn tính chất:
41
i/ không âm và liên tục trên tập intD={x ∈ R : gi(x) < 0, i = 1, ..., m} ii/ p(x) → +∞ khi gi(x) → 0−
Hai hàm phạt điểm trong được sử dụng nhiều, do Fiacco và McCormick
m (cid:88)
m (cid:88)
đưa ra, là:
i=1
i=1
p(x) = − (ln(−gi(x)) và p(x) = − 1 gi(x)
Phương pháp hàm phạt điểm trong xuất phát từ một điểm trong x1 của tập chấp nhận được D, giải một dãy các bài toán tối ưu không ràng buộc minψ(x, αk) với điều kiện x ∈ Rn, trong đó ψ(x, αk) = f (x) + αkp(x), {αk} là dãy tham số dương, giảm đơn điệu về 0.
Thuật toán 3.1
Bước chuẩn bị: Cho số ε > 0 đủ bé (để kiểm tra điều kiện dừng của thuật
toán). Chọn một điểm x1 ∈ D thỏa mãn gi(x1) < 0, i = 1, ..., m. Chọn một tham số phạt α1 > 0 và một số µ ∈ (0, 1). Đặt k = 1 Bước lặp k: (k=1, 2, ...)
Bước k1: Xuất phát từ xk giải bài toán tối ưu không ràng buộc minψ(x, αk) với điều kiện x ∈ Rn nhận được nghiệm xk+1, Bước k2:
if αkp(xk+1) < ε Then Dừng thuật toán (lấy xk+1 là nghiệm tối ưu của
bài toán (P)
Else chuyển Bước k3;
Bước k3; Đặt tham số phạt mới αk+1 := µαk. Đặt k := k + 1. Chuyển về Bước lặp k.
Định lí 3.1
Thuật toán 3.1 có các tính chất sau:
i/ ψ(x, α) ≥ f (x) với mọi α > 0, x chấp nhận được;
ii// gi(xk) < 0, i = 1, ..., m iii/ (cid:8)ψ(xk, αk)(cid:9) hội tụ đến giá trị tối ưu của bài toán (P) khi {αk} → 0
và mọi điểm tụ của {xk} là nghiệm tối ưu của (P ).
42
Chứng minh
i/ Hiển nhiên
ψ(xk, αk) = f (x∗) ii/ Nếu gi(xk)=0 , ta suy ra ψ(xk, αk) không hữu hạn. iii/ Nếu x∗ là nghiệm của (P), ta cần chứng minh lim k→∞
Vì αk (cid:38) 0, {ψ(xk, αk)} là dãy giảm và ψ bị chặn dưới (do f liên tục trên tập compac D nên bị chặn dưới) và αkp(x) ≥ 0), nên ta suy ra ψ(xk, αk) hội tụ với ψ0 ≥ f (x∗). Thật vậy:
ψ(xk, αk) ≥ f (xk) ≥ f (x∗)
ψ(xk, αk) ≥ f (x∗) ψ0 = lim k→+∞
Nếu ε = ψ0 − f (x∗) > 0, do f liên tục nên ta có thể chọn ¯x sao cho
f (¯x) < f (x∗) +
ε 2 gi(¯x) < 0, i = 1, ..., m
Lấy k đủ lớn sao cho
m lim i=1
< −αk ε 2 1 gi(¯x)
và
ψ(xk, αk) − ψ0 < 3 2
Ta có
ψ0 ≤ ψ(xk, αk) (vì {ψ(xk, αk)} là dãy giảm hội tụ với ψ0)
≤ ψ(¯x, αk)
i=1
(vì ψ(xk, αk) đạt cực tiểu tại xk) m (cid:88) ≤ f (¯x) − αk 1 gi(¯x)
< f (¯x) +
ε 2 < f (x∗) + ε
= ψ0
43
Điều này mâu thuẫn. Do đó, ε = ψ0 − f (x∗) = 0, tức là ψ0 = f (x∗) Tiếp theo, nếu ¯x là điểm hội tụ của {xk}, tức ¯x = limxk1, theo i/
ψ(xki, αki) ≥ f (xki)
Từ đây ta suy ra
ψ0 = limψ(xki, αki) ≥ f (¯x)
(cid:3) Điều này có nghĩa là ¯x là cực tiểu.
Chú ý Điểm cực tiểu của ψ(x, α) đều nằm trong tập chấp nhận được.
Mệnh đề 3.1 Cho hàm lồi f : Rn → R và tập lồi khác rỗng D ⊂ Rn. Xét
bài toán min{f (x)|x ∈ D}. Khi đó:
i/ nếu x∗ là một nghiệm tối ưu địa phương của bài toán này thì x∗ cũng
là nghiệm tối ưu toàn cục;
ii/ Nếu x∗ là nghiệm tối ưu địa phương chặt hoặc f là hàm lồi chặt thì
x∗ là nghiệm tối ưu toàn cục duy nhất của bài toán.
Định lí 3.2
Giả sử f là hàm lồi khả vi trên Rn. Khí đó, x∗ ∈ Rn là nghiệm cực tiểu toàn cục của bài toán min{f (x)|x ∈ Rn}, f : Rn → R là hàm phi tuyến khi
và chỉ khi (cid:53)f (x∗) = 0.
a1. Thuật toán Fiacco - Cormick
Hàm mục tiêu: f (x) → min(max)
Các ràng buộc:
gi(x) ≥ 0; i = 1, 2, ..., I
hj(x) = 0; j = 1, 2, ..., J
J (cid:88)
I (cid:88)
Hàm mục tiêu mở rộng được xây dựng ở dạng hàm phạt trong:
j=1
i=1
± . P (x, rk) = f (x) ± rk. h2 j(x) 1 gi(x) 1 √ rk
44
rk là dãy số thực, đơn điệu giảm trong quá trình tính: r0 > r1 > r2 > ... > 0.
Ba cách chọn giá trị đầu r0 khi tính P (x, rk) → min a) rk = 1
1 gi(x)
I (cid:80) i=1
b) Đặt p(x) =
Lấy x(0) là một điểm trong của miền chấp nhận được,
I (cid:88)
i=1
(cid:16) x(0)(cid:17) R = 1 (cid:0)x(0)(cid:1)
r0 = gi − (cid:53) f (cid:0)x(0)(cid:1) . (cid:53) f (cid:0)x(0)(cid:1) (cid:2)(cid:53)p (cid:0)x(0)(cid:1)(cid:3)2
c) Chọn
(cid:115) (cid:16) x(0)(cid:17) (cid:54)= 0 r0 = (cid:53)T f (cid:0)x(0)(cid:1) .H −1. (cid:53) f (cid:0)x(0)(cid:1) (cid:53)T p (cid:0)x(0)(cid:1) .H −1. (cid:53) p (cid:0)x(0)(cid:1) ; (cid:53)p
Trong đó H là ma trận Hessian của hàm p(x) tại x(0).
Thuật toán được mô tả bằng ngôn ngữ lập trình MATLAB như sau:
function mfc=method_fiacco_cormick(epxilon, k)
format short e;
clc;
r=1;
muy=1/2;
X=[0; 1/2]; ss1=10;
while ss1> epxilon
count=0; ss=10; ep=1/2; lamdak =2/3; Xluu=[1; 1];
while and(ss > epxilon, count < k)
count=count + 1;
grad=[dh1px(X, r); dh1py(X, r)];
alpha=1;
dk=-grad;
45
%+ xac dinh buoc di toi uu lamdak
X= X+lamdak*dk;
ss=chuan1(X - Xluu, 2);
Xluu = X;
end;
ss1=abs(r*g(X));
r = r*muy;
end;
mcf=X;
count;
Các kết quả thực hiện thuật toán được cho trong bảng dữ liệu sau đây
Bảng 3.3: Nghiệm xấp xỉ tối ưu sau các bước lặp
f (x, y) = x2 + y2
g(x, y) = x + y − 2
Số bước lặp Sai số Nghiệm xấp xỉ
2 1.5 [-2xe+1; -2xe+1]
4 0.9 [-033; -033]
6 0.15 [-2xe-2; -2xe-2]
8 [-6xe-2; -6xe-2]
10 [-1 x e-3; -1 x e-3]
12 [-4 x e-4; -4 x e-4]
14 [-6 x e-5; -6 x e-5]
16 [-7 x e-6; -7 x e-6]
18 [-1 x e-6; -1 x e-6]
20 0.01 3×e − 3 6×e − 4 10×e − 5 3×e − 5 9×e − 6 9×e − 7 [-5 x e-8; -5 x e-8]
a2. Thuật toán Carroll
Hàm mục tiêu: f (x) → min(max)
46
m (cid:88)
Các ràng buộc: gi(x) ≥ 0; i = 1, 2, ..., m Hàm mục tiêu mở rộng được xây dựng ở dạng hàm phạt trong:
j=1
Pk = P (x, rk) = f (x) ± rk. wj gi(x)
Dấu + khi tìm minf(x), dấu - khi tìm maxf(x);
rk là nhân tử ở bước lặp thứ k; wj là trọng số (thường chọn wj = 1) Thuật toán: Bước 0: Khởi động x0 ∈ D; r0 = 1; µ ∈ (0, 1); ε = 10(−10) Bước lặp: Mọi k=0, 1, 2, ...
wj gi(x) → min
m (cid:80) j=1
+ Giải bài toán không ràng buộc Pk = P (x, rk) = f (x)±rk.
(Sử dụng thuật toán Gradient và quy tắc Armijo)
wj gi(x)
m (cid:80) j=1
+ Kiểm tra: Nếu ≤ ε. (Trong cài đặt thường chọn điều kiện rk. (cid:12) (cid:12) (cid:12) (cid:12) (cid:12)
(cid:12) (cid:12) (cid:12) (cid:12) (cid:12) rk < ε) thì thuật toán dừng. Ngược lại hiệu chỉnh rk+1 = µrk Thuật toán Carroll được mô tả bằng ngôn ngữ MATLAB như sau:
function mc=method_carroll(epxilon, k)
format short e
clc;
r=1;
X=[2; 2]; while r>10∧ -10
count=0; ss=10; ep=1/2; lamda =2/3; Xluu=[1; 1];
while and(ss > epxilon, count < k)
count=count + 1;
grad=[dh1px(X, r); dh1py(X, r)];
alpha=1;
dk=-grad;
while p(X+ alpha*dk, r) - p(X, r) > ep*alpha*(grad.*dk)
47
alpha=alpha*lamda;
grad=[dh1px(X, r); dh1py(X, r)];
dk = -grad;
end;
lamdak=alpha;
X= X+lamdak*dk;
ss=chuan1(X - Xluu, 2);
Xluu = X;
end;
r=r/2
end;
mcf=X;
Các kết quả thực hiện thuật toán được cho trong bảng dữ liệu sau đây
Bảng 3.4: Nghiệm xấp xỉ tối ưu sau các bước lặp
f (x, y) = x2 + y2
g(x, y) = x + y − 2
Số bước lặp Sai số Nghiệm xấp xỉ
1 1.5 [-0.59; -0.59]
2 0.9 [033; 033]
3 0.15 [0.48; 0.48]
4 [0.49; 0.49]
5 [0.499; 0.499]
6 [0.5; 0.5]
9 0.01 3×e − 4 6×e − 8 10×e − 7 [-0,5; 0.5]
B. Phương pháp hàm phạt điểm ngoài
48
m (cid:88)
Hàm phạt p(x) định nghĩa bởi
i=1
p(x) = θ(gi(x))
trong đó θ là một hàm một biến liên tục và thỏa mãn
θ(y) = 0 nếu y ≤ 0 và θ(y) > 0 nếu y > 0.
Thông thường, hàm phạt p(x) có dạng
m (cid:80) i=1
m (cid:80) i=1
p(x) := max{0, gi(x)} hoặc p(x) := [max {0, gi(x)}]2
Hàm mục tiêu của dãy bài toán tối ưu không ràng buộc tương ứng với bài
toán (P) là ψ(x, αk) = f (x) + αkp(x) Trong đó dãy tham số {αk} là dãy số dương, đơn điệu tăng đến +∞. Đại lượng αkp(x) là lượng phạt. Dễ thấy rằng khi x là phương án chấp nhận được của bài toán (P) thì nó không bị phạt (lượng phạt bằng 0),
nhưng khi x không phải phương án chấp nhận được thì nó phải chịu một
lượng phạt là αkp(x). Thuật toán 3.3
Bước chuẩn bị: Cho số ε > 0 đủ bé (để kiểm tra điều kiện dừng của thuật toán). Chọn một điểm x1 ∈ Rn. Chọn một tham số phạt α1 > 0 và một số µ ∈ (0, 1). Đặt k = 1;
Bước lặp k: (k=1, 2, ...)
Bước k1: Giải bài toán tối ưu không ràng buộc minψ(x, αk) với điều kiện x ∈ Rn
Gọi nghiệm của bài toán này là xk+1. Bước k2:
if αkp(xk+1) < ε Then Dừng thuật toán (lấy xk+1 là nghiệm tối ưu của
bài toán (P)
Else chuyển Bước k3;
Bước k3: Đặt tham số phạt mới αk+1 := µαk. Đặt k := k + 1. Chuyển về Bước k.
49
Định lí 3.3
Thuật toán 3.3 có các tính chất sau: i/ ψ(x, α) ≥ f (x) với mọi α > 0, x ∈ Rn;
ii// {ψ(xk, αk)} tăng nếu αk tăng; iii/ {ψ(xk, αk)} hội tụ tới giá trị tối ưu của bài toán (P), khi αk → +∞
và mọi điểm tụ của xk là nghiệm tối ưu của (P ).
Phương pháp Pietrzykowski
Hàm mục tiêu: f (x) → min(max)
Các ràng buộc: gi(x) ≥ 0; i = 1, 2, ..., I
hj(x) = 0; j = 1, 2, ..., J
I (cid:88)
J (cid:88)
Hàm mục tiêu mở rộng được xây dựng ở dạng hàm phạt ngoài:
i (x) −
i=1
j=1
P (x, r) = r.f (x) − wig2 h2 j(x)
Trong đó: wi = 1 đối với gi(x) ≥ 0 wi = 0 đối với gi(x) < 0
Thuật toán được mô tả bằng ngôn ngữ MATLAB như sau
function mp=pietrzykowski_chuan(epxilon, k)
format short
clc;
t=4/3;
X=[0; 0]; ss1=p(X,t)-f(X); Xluu1=[10; 10]; muy=2;
while ss1 > epxilon
count=0; ss=10; ep=1/2; lamda =2/3; Xluu=[1; 1];;
while and(ss > epxilon, count < k)
count=count + 1;
grad=[dh1px(X, t); dh1py(X, t)];
alpha=1;
50
dk=-grad;
while p(X+ alpha*dk, t) - p(X, t) > ep*alpha*(grad.*dk)
alpha=alpha*lamda;
grad=[dh1px(X, t); dh1py(X, t)];
dk = -grad;
end;
lamdak=alpha;
X= X+lamdak*dk;
ss=chuan1(X - Xluu, 2);
Xluu = X;
end;
ss2=chuan1(X-Xluu1, 2)
Xluu1 = X;
ss1=p(X, t)-f(X)
t=t*muy
end;
mcf=X
Các kết quả thực hiện thuật toán được cho trong bảng dữ liệu sau đây
Bảng 3.5: Nghiệm xấp xỉ tối ưu sau các bước lặp
f (x, y) = x2 + y2
g(x, y) = x + y − 3
51
Số bước lặp Sai số Nghiệm xấp xỉ
1.5 [1.0; 1.0] 2
0.9 [1.33; 1.33] 4
0.15 [1.45; 1.45] 5
[1.47; 1.47] 6
[1.48; 1.48] 10
[1.4996; 1.4996] 12
[1.4999; 1.4999] 14
[1.5; 1.5] 16
[1.5; 1.5] 18
0.01 2×e − 2 1×e − 2 3×e − 3 1×e − 4 3×e − 6 3×e − 7 [1.5; 1.5] 20
52
Kết luận
Nội dung chính của luận văn là nghiên cứu các thuật toán giải gần đúng
bài toán quy hoạch phi tuyến có ràng buộc, các kết quả chính luận văn đã
đạt được:
1. Trình bày các thuật toán cơ bản giải mô hình bài toán quy hoạch
tuyến tính bao gồm: Thuật toán hình học, thuật toán đơn hình.
2. Nghiên cứu mô hình bài toán quy hoạch lồi tổng quát.
3. Nghiên cứu các thuật toán tìm nghiệm tối ưu của hàm lồi một biến
số, cài đặt các thuật toán trên ngôn ngữ lập trình MATLAB.
4. Nghiên cứu mô hình bài toán quy hoạch phi tuyến tổng quát không
ràng buộc. Một số thuật toán tổng quat có sử dụng đạo hàm (Gradient,
đường dốc nhất, Newton); các thuật toán không dùng đạo hàm (phương
pháp tìm trực tiếp, phương pháp Powell, phương pháp Nelder và Mead, ...)
5. Nghiên cứu mô hình bài toán quy hoạch phi tuyến tổng quát, có ràng
buộc. Khái niệm hàm Lagrange, các thuật toán thuộc lớp hàm phạt.
Các thuật toán đưa trong luận văn đều đã được nghiên cứu và cài đặt
chi tiết trên ngôn ngữ lập trình MATLAB. Các số liệu thực nghiệm chứng
tỏ các thuật toán là chính xác, tốc độ hội tụ nhanh.
Hướng phát triển tiếp theo của luận văn là nghiên cứu các mô hình bài
toán phi tuyến tính trong các lĩnh vực tự động hóa, tin học, sinh học để
có thể ứng dụng các thuật toán đã nghiên cứu.
53
Tài liệu tham khảo
Tiếng Việt
[1] Nguyễn Nhật Lệ, Các bài toán cơ bản của tối ưu hóa và điều khiển tối
ưu, NXB Khoa học và kỹ thuật, 2009.
[2] Nguyễn Đức Nghĩa, Tối ưu hóa, Nhà xuất bản Giáo dục, 2002.
[3] Bùi Thế Tâm, Trần Vũ Thiệu Các phương pháp tối ưu hóa, Nhà xuất
bản Giao thông Vận tải, 1998.
[4] Nguyễn Hải Thanh, Tối ưu hóa (Giáo trình dành cho ngành Tin hoc
và Công nghệ thông tin), Nhà xuất bản Bách Khoa Hà Nội, năm 2006.
[5] Bùi Minh Trí, Quy hoạch toán học, Nhà xuất bản Khoa học và kỹ
thuật, 1999.
Tiếng Anh
[6] J.Cea, Lectures on Optimization - Theory and Algorithms Tata Insti-
tute of Fuzdamental Ressarch, 1978
[7] R, Fletchers, Practical Methods of Optimization, Wiley, 2000
[8] C.T. Kelley, Iterative Methods of Optimization, SIAM, 1999
[9] J. Nocedal, S. J. Wright, Numerical Optimization, Springer, 1999
