Bài giảng Tối ưu hóa trong thiết kế cơ khí - Chương 11: Ứng dụng Matlab giải các bài toán tối ưu hóa
Giới thiệu về Optimization Toolbox 10
STT Loại bài toán Hàm sử
dụng
1 Cực tiểu hóa hàm 1 biến số
(Scalar Minimization)
fminbnd
2 Cực tiểu hóa hàm nhiều biến số không có ràng buộc
(Unconstrained Minimization)
fminunc
fminsearch
3 Quy hoạch tuyến tính
(Linear Programming)
linprog
4 Quy hoạch bậc hai
(Quadratic Programming)
quadprog
5
Cực tiểu hóa hàm phi tuyến với các ràng buộc
tuyến tính và phi tuyến
(Constrained Minimization)
fmincon
6 Cực tiểu hóa nửa vô hạn
(Semi-Infinite Minimization) fseminf
Bạn đang xem 20 trang mẫu của tài liệu "Bài giảng Tối ưu hóa trong thiết kế cơ khí - Chương 11: Ứng dụng Matlab giải các bài toán tối ưu hóa", để tải tài liệu gốc về máy hãy click vào nút Download ở trên
Tóm tắt nội dung tài liệu: Bài giảng Tối ưu hóa trong thiết kế cơ khí - Chương 11: Ứng dụng Matlab giải các bài toán tối ưu hóa
Trường Đại học Công nghiệp thành phố Hồ Chí Minh Khoa Công nghệ Cơ khí CHƯƠNG 11: ỨNG DỤNG MATLAB GIẢI CÁC BÀI TOÁN TỐI ƯU HÓA Thời lượng: 3 tiết 2 Thống nhất phiên bản MATLAB Trường câu lệnh Các biến tính toán xong Trường lịch sử các câu lệnh đã dùng Chọn thư mục làm việc 3 Làm quen với MATLAB Dấu >> Ký tự nhắc mặc định trong MATLAB, đầu câu lệnh >> Dấu ; Dấu chấm phẩy ở cuối dòng sẽ tránh việc in kết quả ra trường câu lệnh (khi ta không cần in kết quả ra cho cửa sổ ngắn gọn) >> a=2; Dấu Dấu ba chấm ở cuối dòng cho phép tiếp tục code ở dòng tiếp theo >> a= 2 help tên_câu_lệnh Hiển thị các thông tin chi tiết về câu lệnh mà người dùng cần sử dụng >>help linprog Chữ cái viết thường và viết hoa được phân biệt khác nhau trong MATLAB >> a=2 >>A=2 >>A+a MATLAB coi tất cả các biến đều ở dạng mảng (arrays) Tên biến Trong MATLAB tên biến được bắt đầu từ chữ cái và có chiều dài tối đa 31 ký tự bao gồm chữ cái (in hoa và viết thường là khác nhau), số và dấu gạch dưới >>A_a_bc_9=12 Các phép toán thông dụng + - * / ^ >>2+3*6^2/4-7 Tránh trùng với các tên biến tích hợp của hệ thống, như các hằng số, tên hàm pi sin cos v.v.. >>pi 4 Ma trận trong MATLAB Đối tượng Cách thực hiện Ví dụ Véctơ hàng Ngoặc vuông, Dấu cách A=[1 2 3 4] Véctơ cột Cách 1: Ngoặc vuông và xuống dòng B=[5 2 ] Cách 2: Ngoặc vuông và dấu chấm phẩy B=[5;2;3] Cách 3: Đảo véc tơ hàng thành cột bằng dấu ‘ B=[5 2 3]' Ma trận [mxn] Cách 1: áp dụng véc tơ hàng và cột: dấu cách và xuống dòng A=[1 2 3 4 5 6 7 8 9] Cách 2: Dấu cách và dấu chấm phẩy A=[1 2 3;4 5 6; 7 8 9] Ma trận 1 đơn vị ở đường chéo - eye(m,n) – ma trận mxn - eye(n) – ma trận vuông kích thước n eye(4) eye(3,4) Ma trận toàn 1 ones(m,n) hoặc ones(n) ones(3) Ma trận toàn 0 zeros(m,b) hoặc zeros(n) zeros(4) Tạo ra 1 dãy cấp số cộng M:icr:N M:N – khi icr=1 theo mặc định 100:-7:50 50:100 5 Ma trận trong MATLAB (tiếp) Đối tượng Cách thực hiện Ví dụ Phần tử trong ma trận Dấu ngoặc tròn: - A(i,j) chọn phần tử hàng i, cột j của ma trận A - A(m1:m2,n1:n2) Chọn các phần tử từ hàng m1 đến hàng m2, cột n1 đến cột n2 A(2,3) A(2:3,1:3) 6 Kịch bản (Scripts) Là một dạng M-file đơn giản nhất, không có biến vào và biến ra. Nó chỉ gồm một chuỗi các trình tự câu lệnh. 7 Hàm số (Function) function [O_1,O_2,,O_m] = function_name (I_1,I_2,,I_n) %----------------------------------------------------------------------- % Ở dưới ghi function body end 8 function [ p ] = func1( x ) %func1 of this function goes here %Detailed explanation goes here p=x^3-2*x+cos(x); end 9 function [x1,x2] = PTB2(a,b,c) % Ham PTB2 dung de giai phuong trinh bac 2 co dang a*x^2+b*x+c=0 % Tham bien dau vao la 3 he so a, b, c trong do a !=0 Delta = b^2-4*a*c; x1 = (-b+sqrt(Delta))/(2*a); x2 = (-b-sqrt(Delta))/(2*a); [x1 x2] end [x1, x2] = feval('PTB2',1,2,-5) 10 Giới thiệu về Optimization Toolbox STT Loại bài toán Hàm sử dụng 1 Cực tiểu hóa hàm 1 biến số (Scalar Minimization) fminbnd 2 Cực tiểu hóa hàm nhiều biến số không có ràng buộc (Unconstrained Minimization) fminunc fminsearch 3 Quy hoạch tuyến tính (Linear Programming) linprog 4 Quy hoạch bậc hai (Quadratic Programming) quadprog 5 Cực tiểu hóa hàm phi tuyến với các ràng buộc tuyến tính và phi tuyến (Constrained Minimization) fmincon 6 Cực tiểu hóa nửa vô hạn (Semi-Infinite Minimization) fseminf 11 Các thuật toán của các công cụ 12 Các thuật toán của các công cụ 13 CỰC TIỂU HÓA HÀM MỘT BIẾN SỐ min;f x a x b [x,fval,exitflag,output] = fminbnd(@Objfun,a,b,options) x – xuất ra giá trị x làm cho hàm mục tiêu đạt cực tiểu fval – xuất ra giá trị hàm mục tiêu tại điểm cực tiểu x exitflag – xuất ra giá trị để xác định điều kiện dừng tính toán, cụ thể là: • exitflag=1: có nghĩa là hàm đã hội tụ tại điểm lời giải x nếu Stopping Criteria = 'TolX‘ • exitflag=0: có nghĩa là số lượng tính hàm mục tiêu hoặc số lượng vòng lặp đã đạt ngưỡng cho phép nếu Stopping Criteria = 'MaxIter' hoặc 'MaxFunEvals' • exitflag=-1: có nghĩa là thuật toán bị dừng vì hàm đầu ra • exitflag=-2: có nghĩa là khoảng giá trị bị sai (a>b) output – xuất ra các thông tin về số vòng lặp tính toán, số lần tính hàm số, các thuật toán tại các bước tính và thông báo cuối cùng Objfun – tên của M-file xác định hàm mục tiêu a,b – giá trị 2 biên của biến x options – các thuộc tính cần thiết cho việc giải bài toán bằng hàm fminbnd, được xác định trước dòng cú pháp 14 CỰC TIỂU HÓA HÀM MỘT BIẾN SỐ Tìm cực tiểu hàm số 2 0.75 1 0.65 0.65 arctan ; 0;0.5 1 f x x x x x 1) Bước 1: Tạo 1 thư mục cho bài toán, ví dụ fminbnd1 Vào trong thư mục 15 2) Bước 2: Tạo M-file hàm mục tiêu, ví dụ Objfun.m function f = Objfun( x ) %OBJFUN Summary of this function goes here % Detailed explanation goes here f= 0.65 - (0.75/(1+x^2))- 0.65*x*atan(1/x); end Sửa lại code của hàm 16 3) Bước 3: Tạo M-file Script lời giải, ví dụ Solve.m clear;clc;format long;warning('off'); % Nhap mien xac dinh cua bien x: a<=x<=b a=0; b=0.5; % Chon Algorithm (Chon mot trong so duoi day) %Alg='active-set'; %Alg='trust-region-reflective'; Alg='interior-point'; %Alg='levenberg-marquardt'; %Alg='trust-region-dogleg'; %Alg='lm-line-search'; % Chon Stopping Criteria GTN=1e-4;GTL=1e4; % Neu Stopping Criteria la sai so cua tham bien StCr='TolX'; % Neu Stopping Criteria la so luong Iterations %StCr='MaxIter'; % Neu Stopping Criteria la so luong tinh cac ham so %StCr='MaxFunEvals'; if strcmp(StCr,'MaxIter')==1 || strcmp(StCr,'MaxFunEvals')==1 GT=GTL; elseif strcmp(StCr,'TolX')==1 GT=GTN; end options = optimset('Algorithm',Alg,'Display','iter',StCr,GT,'PlotFcns',@optimplotfval); [x,fval,exitflag,output]=fminbnd(@Objfun,a,b,options) 17 0 1 2 3 4 5 6 7 8 9 -0.32 -0.31 -0.3 -0.29 -0.28 -0.27 -0.26 -0.25 -0.24 -0.23 Iteration F u n c ti o n v a lu e Current Function Value: -0.31002 18 CỰC TIỂU HÓA HÀM NHIỀU BIẾN SỐ KHÔNG RÀNG BUỘC 1 2min; , , , nf x x x x x [x,fval,exitflag,output, grad,hessian] = fminunc(@Objfun,x0,options) x – xuất ra véctơ tham biến x làm cho hàm mục tiêu đạt cực tiểu fval – xuất ra giá trị hàm mục tiêu tại điểm cực tiểu x exitflag – xuất ra giá trị để xác định điều kiện dừng tính toán, cụ thể là: • exitflag>0: có nghĩa là hàm đã hội tụ tại điểm lời giải x • exitflag=0: có nghĩa là số lượng tính hàm mục tiêu hoặc số lượng vòng lặp đã đạt ngưỡng cho phép nếu Stopping Criteria = 'MaxIter' hoặc 'MaxFunEvals' • exitflag<0: Hàm không hội tụ tại điểm lời giải output – xuất ra các thông tin về số vòng lặp tính toán, số lần tính hàm số, các thuật toán tại các bước tính, v.v grad – xuất ra véctơ Gradient của hàm số tại điểm lời giải hessian – xuất ra ma trận Hessian của hàm số tại điểm lời giải x0 – véc tơ tham biến khởi đầu Objfun - tên của M-file xác định hàm mục tiêu 19 Tìm cực tiểu hàm số: 1 2 1 2 1 2 50 20 , minf f x x x x x x x 1) Bước 1: Tạo 1 thư mục cho bài toán, ví dụ MultiUncon1 2) Bước 2: Tạo M-file hàm mục tiêu, ví dụ Objfun.m function f = Objfun( x ) %OBJFUN Summary of this function goes here % Detailed explanation goes here f = x(1)*x(2) + 50/x(1) + 20/x(2); end 3) Bước 3: Tạo M-file Script lời giải, ví dụ Solvefminunc.m 20 clear;clc;format long;warning('off'); % Nhap vecto tham bien khoi dau: x0 = [1,1]; % Chon Algorithm (Chon mot trong so duoi day) %Alg='active-set'; %Alg='trust-region-reflective'; %Alg='interior-point'; %Alg='levenberg-marquardt'; Alg='trust-region-dogleg'; %Alg='lm-line-search'; % Chon Stopping Criteria GTN=1e-4;GTL=1e4; % Neu Stopping Criteria la sai so cua tham bien %StCr='TolX'; % Neu Stopping Criteria la sai so cua ham so StCr='TolFun'; % Neu Stopping Criteria la so luong Iterations %StCr='MaxIter'; % Neu Stopping Criteria la so luong tinh cac ham so %StCr='MaxFunEvals'; if strcmp(StCr,'MaxIter')==1 || strcmp(StCr,'MaxFunEvals')==1 GT=GTL; elseif strcmp(StCr,'TolX')==1 || strcmp(StCr,'TolFun')==1 GT=GTN; end options = optimset('Algorithm',Alg,'Display','iter',StCr,GT,'PlotFcns',@optimplotfval); [x,fval,exitflag,output, grad,hessian] = fminunc(@Objfun,x0,options) 21 22 0 2 4 6 8 10 12 14 30 35 40 45 50 55 60 65 70 75 Iteration F u n c ti o n v a lu e Current Function Value: 30 23 CỰC TIỂU HÓA HÀM NHIỀU BIẾN SỐ KHÔNG RÀNG BUỘC 1 2min; , , , nf x x x x x [x,fval,exitflag,output] = fminsearch(@Objfun,x0,options) x – xuất ra véctơ tham biến x làm cho hàm mục tiêu đạt cực tiểu fval – xuất ra giá trị hàm mục tiêu tại điểm cực tiểu x exitflag – xuất ra giá trị để xác định điều kiện dừng tính toán, cụ thể là: • exitflag>0: có nghĩa là hàm đã hội tụ tại điểm lời giải x • exitflag=0: có nghĩa là số lượng tính hàm mục tiêu hoặc số lượng vòng lặp đã đạt ngưỡng cho phép nếu Stopping Criteria = 'MaxIter' hoặc 'MaxFunEvals' • exitflag<0: Hàm không hội tụ tại điểm lời giải output – xuất ra các thông tin về số vòng lặp tính toán, số lần tính hàm số, các thuật toán tại các bước tính, v.v x0 – véc tơ tham biến khởi đầu Objfun - tên của M-file xác định hàm mục tiêu Cú pháp này sử dung pp Nelder-Mead và chỉ tìm được cực tiểu local 24 4) Bước 4: Tạo M-file Script lời giải, ví dụ Solvefminsearch.m clear;clc;format long;warning('off'); % Nhap vecto tham bien khoi dau: x0 = [1,1]; % Chon Algorithm (Chon mot trong so duoi day) %Alg='active-set'; %Alg='trust-region-reflective'; %Alg='interior-point'; %Alg='levenberg-marquardt'; Alg='trust-region-dogleg'; %Alg='lm-line-search'; % Chon Stopping Criteria GTN=1e-4;GTL=1e4; % Neu Stopping Criteria la sai so cua tham bien %StCr='TolX'; % Neu Stopping Criteria la sai so cua ham so StCr='TolFun'; % Neu Stopping Criteria la so luong Iterations %StCr='MaxIter'; % Neu Stopping Criteria la so luong tinh cac ham so %StCr='MaxFunEvals'; if strcmp(StCr,'MaxIter')==1 || strcmp(StCr,'MaxFunEvals')==1 GT=GTL; elseif strcmp(StCr,'TolX')==1 || strcmp(StCr,'TolFun')==1 GT=GTN; end options = optimset('Algorithm',Alg,'Display','iter',StCr,GT,'PlotFcns',@optimplotfval); [x,fval,exitflag,output] = fminsearch(@Objfun,x0,options) 25 Optimization terminated: the current x satisfies the termination criteria using OPTIONS.TolX of 1.000000e-004 and F(X) satisfies the convergence criteria using OPTIONS.TolFun of 1.000000e-004 26 0 5 10 15 20 25 30 35 40 45 30 35 40 45 50 55 60 65 70 75 Iteration F u n c ti o n v a lu e Current Function Value: 30 27 QUY HOẠCH TUYẾN TÍNH minT m n p n f x c x A x a B x b lb x ub 1 1 1 1 1 1 2 2 2 2 2 2 ; ; ; ; ; ; n n n n n n c x a b l u c x a b l u c x a b l u c x a b lb ub 11 12 1 21 22 2 1 2 11 12 1 21 22 2 1 2 n n m n m m mn n n p n p p pn a a a a a a a a a b b b b b b b b b A B 28 [x,fval,exitflag,output,lambda] = linprog(c,A,a,B,b,lb,ub,x0,options) x – xuất ra véctơ tham biến x làm cho hàm mục tiêu đạt cực tiểu fval – xuất ra giá trị hàm mục tiêu tại điểm cực tiểu x exitflag – xuất ra giá trị để xác định điều kiện dừng tính toán, cụ thể là: • exitflag=1: có nghĩa là hàm đã hội tụ tại điểm lời giải x • exitflag=0: có nghĩa là số lượng tính hàm mục tiêu hoặc số lượng vòng lặp đã đạt ngưỡng cho phép nếu Stopping Criteria = 'MaxIter' • exitflag=-2: Không tìm được một điểm hợp lệ nào • exitflag=-3: Miền ràng buộc là vô hạn • exitflag=-4: Giá trị NaN (không phải là số) xuất hiện trong quá trình tính toán • exitflag=-5: Cả bài toán gốc và đối ngẫu của nó đều vô nghiệm • exitflag=-7: Hướng tìm kiếm càng lúc càng nhỏ, không thể tìm được các điểm tốt hơn output – xuất ra các thông tin về số vòng lặp tính toán, số lần tính hàm số, các thuật toán tại các bước tính, v.v x0 – véc tơ tham biến khởi đầu c,A,a,B,b,lb,ub – lần lượt các véc tơ, ma trận dữ liệu bài toán 29 Giải bài toán quy hoạch tuyến tính sau: 1 2 8 2 0 ; ; 10 ; ; ; ; 3 0 3 x x c x a b lb ub 3 2 1 2 2 1 ; 0 1 A B Đổi lại yêu cầu của mục tiêu: 1 22 3 minf x x x Nhận diện được các véc tơ và ma trận dữ liệu đề bài: 1) Bước 1: Phân tích bài toán 2) Bước 2: Tạo M-file script cho bài toán, ví dụ LP1.m 30 clear;clc;format long;warning('off'); c=[-2,-3] A=[1,2; 2,1; 0,1] a=[8,10,3] B=[] b=[] lb=zeros(2,1) ub=[] x0=[] options = optimset('LargeScale', 'off', 'Simplex', 'on','Display','iter') [x,fval,exitflag,output,lambda] = linprog(c,A,a,B,b,lb,ub,x0,options) 31 Giải bài toán quy hoạch tuyến tính sau: 1 1 2 2 minTf x c x A x a A x a B x b lb x ub 1 2 3 1 2 4 5 6 1 2 7 2 5 0 2 3 7 5 10 1 2 ; ; ; ; ; ; ; 4 8 1 37 15 5 8 5 4 6 1 10 x x x x x x c x a a b lb ub 1 2 1 2 3 0 0 0 0 1 2 3 0 0 3 0 0 0 2 1 1 1 1 1 1 1 ; ; 0 0 1 2 3 0 0 4 0 2 0 3 5 0 3 0 1 0 0 0 0 1 2 3 A A B 1) Bước 1: Phân tích bài toán 32 2) Bước 2: Tạo M-file script cho bài toán, ví dụ LP2.m clear;clc;format long;warning('off'); c=[1,-2,3,-4,5,-6]' A1=[1,2,3,0,0,0; 0,1,2,3,0,0; 0,0,1,2,3,0; 0,0,0,1,2,3] A2=[3,0,0,0,-2,1; 0,4,0,-2,0,3] B=[1,1,1,1,1,1; 5,0,-3,0,1,0] a1=[5,7,8,8]' a2=[5,7]' b=[10,15]' A=[-A1; A2] a=[-a1; a2] lb=[-2,0,-1,-1,-5,1]' ub=[7,2,2,3,4,10]' x0=[] options = optimset('LargeScale', 'off', 'Simplex', 'on','Display','iter') [x,fval,exitflag,output,lambda] = linprog(c,A,a,B,b,lb,ub,x0,options) 33 QUY HOẠCH BẬC HAI 1 min 2 T T n n m n p n f x x H x c x A x a B x b lb x ub 11 12 1 21 22 2 1 2 11 12 1 21 22 2 1 2 n n m n m m mn n n p n p p pn a a a a a a a a a b b b b b b b b b A B 1 1 1 1 1 1 2 2 2 2 2 2 ; ; ; ; ; ; n n n n n n c x a b l u c x a b l u c x a b l u c x a b lb ub [H] chính là ma trận Hessian của hàm mục tiêu 34 [x,fval,exitflag,output,lambda] = quadprog(H,c,A,a,B,b,lb,ub,x0,options) x – xuất ra véctơ tham biến x làm cho hàm mục tiêu đạt cực tiểu fval – xuất ra giá trị hàm mục tiêu tại điểm cực tiểu x exitflag – xuất ra giá trị để xác định điều kiện dừng tính toán, cụ thể là: • exitflag=1: có nghĩa là hàm đã hội tụ tại điểm lời giải x • exitflag=3: có nghĩa là sự thay đổi giá trị hàm mục tiêu đã nhỏ hơn sai số cho phép • exitflag=4: có nghĩa là cực tiểu địa phương đã được tìm thấy • exitflag=0: có nghĩa là số lượng tính hàm mục tiêu hoặc số lượng vòng lặp đã đạt ngưỡng cho phép nếu Stopping Criteria = 'MaxIter' • exitflag=-2: Bài toán không giải được • exitflag=-3: Miền ràng buộc là vô hạn • exitflag=-4: Hướng tìm kiếm không phải là hướng tốt, không thể tìm được các điểm tốt hơn • exitflag=-7: Giá trị của hướng tìm kiếm quá nhỏ, không có tiến triển • output – xuất ra các thông tin về số vòng lặp tính toán, số lần tính hàm số, các thuật toán tại các bước tính, v.v x0 – véc tơ tham biến khởi đầu H,c,A,a,B,b,lb,ub – lần lượt các véc tơ, ma trận dữ liệu bài toán 35 2 1 430 2 2 2 1 2 4 1 2 1 3 2 3 2 4 3 4 0 4 8 72 5 93 6 x xx x x x x x x x x x x x x x 4 4 6 12 0 0 0 7 7 4 9 9 5 8 4 5 4 8 H Dựa vào mầu để nhận diện 36 Giải bài toán quy hoạch bậc 2 sau: 1) Bước 1: Phân tích bài toán 1 2 2 2 8 2 0 ; ; 10 ; ; ; ; 3 0 3 1 2 2 0 ; 2 1 ; 0 4 0 1 x x c x a b lb ub H A B 2) Bước 2: Tạo M-file script cho bài toán, ví dụ QP1.m 37 clear;clc;format long;warning('off'); H=[2,0; 0,4]; c=[2,3]'; A=[1,2; 2,1; 0,1]; a=[8,10,3]'; B=[]; b=[]'; lb=[0,0]'; ub=[inf,inf]'; x0=[]'; % Chon Algorithm (Chon mot trong so duoi day) %Alg='active-set'; %Alg='trust-region-reflective'; %Alg='interior-point'; %Alg='levenberg-marquardt'; Alg='trust-region-dogleg'; %Alg='lm-line-search'; % Chon Stopping Criteria GTN=1e-4;GTL=1e4; % Neu Stopping Criteria la sai so cua tham bien %StCr='TolX'; % Neu Stopping Criteria la sai so cua ham so StCr='TolFun'; % Neu Stopping Criteria la so luong Iterations %StCr='MaxIter'; if strcmp(StCr,'MaxIter')==1 GT=GTL; elseif strcmp(StCr,'TolX')==1 || strcmp(StCr,'TolFun')==1 GT=GTN; end options = optimset('Algorithm',Alg,'Display','iter',StCr,GT); [x,fval,exitflag,output,lambda]=quadprog(H,c,A,a,B,b,lb,ub,x0,options) % Ve do thi x1=[-4:0.1:4]; x2=[-4:0.1:4]; [X1,X2]=meshgrid(x1,x2); f=X1.^2 + 2*X2.^2 + 2*X1 + 3*X2; meshc(X1,X2,f); hold on plot(x(1),x(2),'r*') 38 Giải bài toán quy hoạch bậc 2 sau: 1) Bước 1: Phân tích bài toán 1 2 3 3 3 4 0 6 6 ; ; ; ; 0 ; ; 2 12 0 2 1 0 1 1 1 1 4 2 ; ; 1 1 2 0 2 4 x x x c x a b lb ub H A B 2) Bước 2: Tạo M-file script cho bài toán, ví dụ QP2.m 39 clear; clc; format long; H=[2,1,0; 1,4,2; 0,2,4]; c=[4,6,12]'; A=[1,1,1; -1,-1,2]; a=[6,2]'; B=[]; b=[]'; lb=[0,0,0]'; ub=[inf,inf,inf]'; x0=[]'; % Chon Algorithm (Chon mot trong so duoi day) %Alg='active-set'; %Alg='trust-region-reflective'; %Alg='interior-point'; %Alg='levenberg-marquardt'; Alg='trust-region-dogleg'; %Alg='lm-line-search'; % Chon Stopping Criteria GTN=1e-4;GTL=1e4; % Neu Stopping Criteria la sai so cua tham bien %StCr='TolX'; % Neu Stopping Criteria la sai so cua ham so StCr='TolFun'; % Neu Stopping Criteria la so luong Iterations %StCr='MaxIter'; if strcmp(StCr,'MaxIter')==1 GT=GTL; elseif strcmp(StCr,'TolX')==1 || strcmp(StCr,'TolFun')==1 GT=GTN; end options = optimset('Algorithm',Alg,'Display','iter',StCr,GT); [x,fval,exitflag,output,lambda]=quadprog(H,c,A,a,B,b,lb,ub,x0,options) %Neu tim max f: doi cua f se chinh la gia tri cuc dai %[x,fval,exitflag,output,lambda]=QUADPROG(H,c,A,a,B,b,lb,ub,x0,options) 40 Cực tiểu hóa hàm phi tuyến với các ràng buộc 1 2 min 0; 1.. 0; 1.. j k m n p n T n f g j r h k s x x x x x x A x a B x b lb x ub x 1 1 1 1 1 2 2 2 2 2 ; ; ; ; n n n n n x a b l u x a b l u x a b l u x a b lb ub 11 12 1 21 22 2 1 2 11 12 1 21 22 2 1 2 n n m n m m mn n n p n p p pn a a a a a a a a a b b b b b b b b b A B 41 [x,fval,exitflag,output,lambda,grad,hessian]= fmincon(@Objective,x0,A,a,B,b,lb,ub,@NonLinCon,options) x – xuất ra véctơ tham biến x làm cho hàm mục tiêu đạt cực tiểu fval – xuất ra giá trị hàm mục tiêu tại điểm cực tiểu x exitflag – xuất ra giá trị để xác định điều kiện dừng tính toán, cụ thể là: a) Đối với tất cả các thuật toán: • exitflag=1: Thước đo tối ưu bậc nhất (First-order optimality measure) nhỏ hơn giá trị thông số TolFun và vi phạm ràng buộc lớn nhất (maximum constraint violation) nhỏ hơn giá trị thông số TolCon • exitflag=0: có nghĩa là số lượng tính hàm mục tiêu đã đạt ngưỡng cho phép nếu Stopping Criteria = FunEvals hoặc số lượng vòng lặp đã đạt ngưỡng cho phép nếu Stopping Criteria = 'MaxIter' • exitflag=-1: Hàm đầu ra kết thúc thuật toán • exitflag=-2: Không tìm được lời giải hợp lệ nào b) Đối với 2 thuật toán là ‘trust-region-reflective’ và ‘interior-point’: • exitflag=2: Sự biến thiên của x nhỏ hơn giá trị thông số TolX và vi phạm ràng buộc lớn nhất nhỏ hơn giá trị thông số TolCon c) Đối với riêng thuật toán ‘trust-region-reflective’: • exitflag=3: có nghĩa là sự thay đổi giá trị hàm mục tiêu đã nhỏ hơn sai số cho phép TolFun và vi phạm ràng buộc lớn nhất nhỏ hơn giá trị thông số TolCon 42 d) Đối với riêng thuật toán ‘active-set’: • exitflag=4: giá trị của hướng tìm kiếm (search direction) đã nhỏ hơn 2 lần giá trị cho phép TolX và vi phạm ràng buộc lớn nhất nhỏ hơn giá trị thông số TolCon • exitflag=5: giá trị của hướng tìm kiếm (search direction) đã nhỏ hơn 2 lần giá trị cho phép TolFun và vi phạm ràng buộc lớn nhất nhỏ hơn giá trị thông số TolCon e) Đối với riêng thuật toán ‘interior-point’: • exitflag=-3: Điểm x đang xét đi xuống dưới giá trị ObjectiveLimit và vi phạm ràng buộc lớn nhất nhỏ hơn giá trị thông số TolCon lambda – Giá trị toán tử Lagrange tại điểm lời giải grad – xuất ra véctơ Gradient của hàm số tại điểm lời giải hessian – xuất ra ma trận Hessian của hàm số tại điểm lời giải output – xuất ra các thông tin về số vòng lặp tính toán, số lần tính hàm số, các thuật toán tại các bước tính, v.v x0 – véc tơ tham biến khởi đầu A,a,B,b,lb,ub – lần lượt các véc tơ, ma trận dữ liệu bài toán Objective - tên của M-file xác định hàm mục tiêu NonLinCon - tên của M-file xác định hàm ràng buộc phi tuyến đẳng thức và bất đẳng thức 43 Giải bài toán tối ưu hóa phi tuyến với các ràng buộc sau 1 2 1 2 1 1 2 2 2 2 1 2 1 2 3 1 2 2 2 1 2 1 1 2 1 2 1 2 1 2 , 9.82 2 min 2500 500 0 2500 0 0.5882 0.3 48.59 0 10 1.056 0 3 18 12 3 8 10 2.531 2 14;0.2 0.8 f f x x x x x x x x x x x x x x x x x x x x x x x x x x x 2 2 1 2 1 2 3 18 1 3 1 10 12 ; ; 2.531 ; 8 2 14 ; 0.2 0.8 x x A B x a b lb ub 44 1) Bước 1: Tạo 1 thư mục cho bài toán, ví dụ NonlinCon1 2) Bước 2: Tạo M-file hàm mục tiêu, ví dụ Objective.m function f = Objective( x ) %OBJECTIVE Summary of this function goes here % Detailed explanation goes here f= 9.82*x(1)*x(2)+2*x(1); end 3) Bước 3: Tạo M-file hàm các hàm ràng buộc phi tuyến đẳng thức và bất đẳng thức, ví dụ NonLinCon.m function [ g, h ] = NonLinCon( x ) %NONLINCON Summary of this function goes here % Detailed explanation goes here % Cac rang buoc bat dang thuc g = [2500/(pi*x(1)*x(2))- 500; 2500/(pi*x(1)*x(2))- (pi^2*(x(1)^2+x(2)^2))/0.5882]; % Cac rang buoc dang thuc h = [0.3*x(1)^3*x(2)-48.59*x(2); 10*x(2)/x(1) - x(1)*x(2) + 1.056]; end 4) Bước 4: Tạo M-file Script lời giải, ví dụ Solve.m 45 clear;clc;format long;warning off; A=[3,-18; 1,3]; a=[12,8]'; B=[1,-10]; b=[2.531]'; lb=[2,0.2]'; ub=[14,0.8]'; x0 = [lb(:)+(ub(:)-lb(:))*rand(1,1)]' % Diem khoi dau % Kiem tra gia tri cua ham muc tieu cung nhu cac rang buoc phi tuyen tai diem khoi dau fprintf ('Gia tri ham muc tieu va rang buoc phi tuyen tai diem khoi dau\n'); f=Objective(x0) [g, h] = NonLinCon(x0) % Chon Algorithm (Chon mot trong so duoi day) %Alg='active-set'; %Alg='trust-region-reflective'; Alg='interior-point'; % Chon Stopping Criteria GTN=1e-6;GTL=1e6; % Neu Stopping Criteria la sai so cua tham bien %StCr='TolX'; % Neu Stopping Criteria la sai so cua ham so %StCr='TolFun'; % Neu Stopping Criteria la so luong Iterations StCr='MaxIter'; % Neu Stopping Criteria la so luong tinh cac ham so %StCr='MaxFunEvals'; if strcmp(StCr,'MaxIter')==1 || strcmp(StCr,'MaxFunEvals')==1 GT=GTL; elseif strcmp(StCr,'TolX')==1 || strcmp(StCr,'TolFun')==1 GT=GTN; end options = optimset('Algorithm',Alg,'Display','iter-detailed',StCr,GT,'PlotFcns',@optimplotfval); [x,fval,exitflag,output,lambda,grad,hessian]=fmincon(@Objective,x0,A,a,B,b,lb,ub,@NonLinCon,options ) fprintf('Gia tri rang buoc phi tuyen tai diem cuc tri\n'); [g, h] = NonLinCon(x) 46 CỰC TIỂU HÓA NỬA VÔ HẠN (Semi-Infinite Minimization) 1 2 1 1 2 2 min 0; 1.. 0; 1.. , 0 , 0 , 0 j k m n p n T n q q f g j r h k s x x x K w K w K w x x x A x a B x b lb x ub x x x x 11 12 1 21 22 2 1 2 11 12 1 21 22 2 1 2 n n m n m m mn n n p n p p pn a a a a a a a a a b b b b b b b b b A B 1 1 1 1 1 2 2 2 2 2 ; ; ; ; n n n n n x a b l u x a b l u x a b l u x a b lb ub 47 [x,fval,exitflag,output,lambda]=fseminf(@Objective,x0,q, myinfcon,A,a,B,b,lb,ub,options) x – xuất ra véctơ tham biến x làm cho hàm mục tiêu đạt cực tiểu fval – xuất ra giá trị hàm mục tiêu tại điểm cực tiểu x exitflag – xuất ra giá trị để xác định điều kiện dừng tính toán, cụ thể là: • exitflag=1: Hàm hội tụ tại điểm lời giải x • exitflag=4: giá trị của hướng tìm kiếm (search direction) đã nhỏ hơn giá trị cho phép và vi phạm ràng buộc lớn nhất nhỏ hơn giá trị thông số TolCon • exitflag=5: giá trị của đạo hàm định hướng (directional derivative) đã nhỏ giá trị cho phép và vi phạm ràng buộc lớn nhất nhỏ hơn giá trị thông số TolCon • exitflag=0: có nghĩa là số lượng tính hàm mục tiêu đã đạt ngưỡng cho phép nếu Stopping Criteria = FunEvals hoặc số lượng vòng lặp đã đạt ngưỡng cho phép nếu Stopping Criteria = 'MaxIter' • exitflag=-1: Hàm đầu ra kết thúc thuật toán • exitflag=-2: Không tìm được lời giải hợp lệ nào x0 – véc tơ tham biến khởi đầu A,a,B,b,lb,ub – lần lượt các véc tơ, ma trận dữ liệu bài toán Objective - tên của M-file xác định hàm mục tiêu myinfcon - tên của M-file xác định hàm ràng buộc phi tuyến đẳng thức, bất đẳng thức và ràng buộc nửa vô hạn 48 myinfcon function [g,h,K1,K2,..,Kq,S] = myinfcon(x,S) % Initial sampling interval if isnan(S(1,1)), S = [s1 0; s2 0; ... sq 0]; % S co q hang and 2 columns end % Sample set w1 = w1_begin:S(1,1):w1_end; % Compute sample set w2 = w2_begin:S(2,1):w2_end; % Compute sample set ... wq = wq_begin:S(q,1):wq_end; % Compute sample set % Semi-infinite constraints K1 = ... % 1st semi-infinite constraint at x and w K2 = ... % 2nd semi-infinite constraint at x and w ... Kq = ...% Last semi-infinite constraint at x and w % Cac rang buoc bat dang thuc phi tuyen g = [ham_g1; ham_g2; ... ham_gr]; % Compute nonlinear inequalities at x % Cac rang buoc dang thuc phi tuyen h = [ham_h1; ham_h2; ... ham_hs]; % Compute nonlinear equalities at x % Plot a graph of semi-infinite constraints % plot(w1,K1,'-',w2,K2,':'),title('Semi-infinite constraints') % drawnow 49 Giải bài toán tối ưu hóa nửa vô hạn sau: min 1) Bước 1: Tạo 1 thư mục cho bài toán, ví dụ Seminf1 2) Bước 2: Tạo M-file hàm mục tiêu, ví dụ Objective.m function f = Objective(x,s) %OBJECTIVE Summary of this function goes here % Detailed explanation goes here f=sum((x-0.5).^2); %f=(x(1)-0.5)^2 + (x(2)-0.5)^2 + (x(3)-0.5)^2; end 3) Bước 3: Tạo M-file hàm các hàm ràng buộc phi tuyến đẳng thức, bất đẳng thức và nửa vô hạn, ví dụ myinfcon.m 50 function [g,h,K1,K2,s] = myinfcon(x,s) % Initial sampling interval if isnan(S(1,1)), s = [0.2 0; 0.2 0]; % S co q hang and 2 columns end % Sample set w1 = 0:s(1,1):100; % Compute sample set w2 = 0:s(2,1):100; % Compute sample set % Semi-infinite constraints K1 = sin(w1*x(1)).*cos(w1*x(2)) - 1/1000*(w1-50).^2 - sin(w1*x(3))-x(3)-1; % 1st semi-infinite constraint at x and w K2 = sin(w2*x(2)).*cos(w2*x(1)) - 1/1000*(w2-50).^2 - sin(w2*x(3))-x(3)-1; % 2nd semi-infinite constraint at x and w % Cac rang buoc bat dang thuc phi tuyen g = []; %Trong bai toan nay khong co 1 rang buoc phi tuyen bat dang thuc nao % Cac rang buoc dang thuc phi tuyen h = []; %Trong bai toan nay khong co 1 rang buoc phi tuyen dang thuc nao % Plot a graph of semi-infinite constraints % plot(w1,K1,'-',w2,K2,':'),title('Semi-infinite constraints') % drawnow 4) Bước 4: Tạo M-file Script lời giải, ví dụ Solve.m 51 clear;clc;format long;warning off; A=[]; %Trong bai tap nay khong co rang buoc tuyen tinh bat dang thuc a=[]'; %Trong bai tap nay khong co rang buoc tuyen tinh bat dang thuc B=[]; %Trong bai tap nay khong co rang buoc tuyen tinh dang thuc b=[]'; %Trong bai tap nay khong co rang buoc tuyen tinh dang thuc lb=[]';%Trong bai tap nay khong co rang buoc phia truoc cua tham bien ub=[]';%Trong bai tap nay khong co rang buoc phia sau cua tham bien x0 = [0.5;0.2;0.3]; % Diem khoi dau % Kiem tra gia tri cua ham muc tieu cung nhu cac rang buoc phi tuyen tai diem khoi dau fprintf ('Gia tri ham muc tieu va rang buoc phi tuyen tai diem khoi dau\n'); f=Objective(x0,NaN) [g,h] = myinfcon(x0,NaN) % Chon Algorithm (Chon mot trong so duoi day) %Alg='active-set'; %Alg='trust-region-reflective'; %Alg='interior-point'; %Alg='levenberg-marquardt'; %Alg='trust-region-dogleg'; Alg='lm-line-search'; % Chon Stopping Criteria GTN=1e-6;GTL=1e6; % Neu Stopping Criteria la sai so cua tham bien %StCr='TolX'; % Neu Stopping Criteria la sai so cua ham so %StCr='TolFun'; % Neu Stopping Criteria la so luong Iterations StCr='MaxIter'; % Neu Stopping Criteria la so luong tinh cac ham so %StCr='MaxFunEvals'; if strcmp(StCr,'MaxIter')==1 || strcmp(StCr,'MaxFunEvals')==1 GT=GTL; elseif strcmp(StCr,'TolX')==1 || strcmp(StCr,'TolFun')==1 GT=GTN; end options = optimset('Algorithm',Alg,'Display','iter-detailed',StCr,GT,'PlotFcns',@optimplotfval); [x,fval,exitflag,output,lambda]=fseminf(@Objective,x0,2,@myinfcon,A,a,B,b,lb,ub,options) fprintf('Gia tri rang buoc phi tuyen tai diem cuc tri\n'); [g, h] = myinfcon(x,NaN)
File đính kèm:
- bai_giang_toi_uu_hoa_trong_thiet_ke_co_khi_chuong_11_ung_dun.pdf