Chương 11: Tối ưu hóa và nội suy

Tài liệu được phát hành theo giấy phép văn bản tự do GNU (FDL). Bản gốc: Physical Modeling in MATLAB—Green Tea Press. Tác giả Allen B. Downey. Người dịch: Nguyễn Quang Chiến.

Các sự kiện dùng trong hàm ODE

Thường thì khi gọi ode45, bạn đều phải chỉ ra các thời điểm bắt đầu và kết thúc. Nhưng trong nhiều trường hợp, bạn không biết trước lúc nào việc mô phỏng cần kết thúc. Thật may là MATLAB cung cấp một cơ chế xử lý vấn đề này. Điều không hay là cơ chế này hơi lủng củng một chút. Sau đây là cách hoạt động:

1. Trước khi gọi ode45 bạn dùng odeset để tạo ra một đối tượng có tên options để chứa các giá trị quy định cách hoạt động của ode45:

   options = odeset('Events', @events);
   

   Trong trường hợp này, tên của tùy chọn (option) là Events còn giá trị là một chuôi hàm. Khi ode45 hoạt động, nó sẽ gọi events sau mỗi bước thời gian. Bạn có thể gọi hàm này bằng bất cứ tên gọi nào, nhưng cái tên events thường được chọn theo thông lệ.

2. Hàm mà bạn cung cấp phải nhận vào cùng các biến như hàm tốc độ nhận vào. Chẳng hạn, sau đây là một hàm sự kiện gắn với hàm projectile từ Mục [tính đường bay].

   function [value,isterminal,direction] = events(t,X)
       value = X(2);        % Extract the current height.
       isterminal = 1;      % Stop the integration if height crosses zero.
       direction = -1;      % But only if the height is decreasing.
   end
   

   events trả lại ba biến đầu ra:

   value quyết định xem khi nào một sự kiện xảy ra. Trong trường hợp này value nhận giá trị của phần tử thứ hai thuộc X, vốn được hiểu là chiều cao của vật bay. Một “sự kiện” (event) là một thời điểm khi mà giá trị chiều cao nói trên đi qua 0.

   direction quyết định liệu rằng một sự kiện có xảy ra khi value đang tăng (direction=1), giảm (direction=-1, hoặc cả hai direction=0.

   isterminal quyết định xem điều gì cần thực hiện khi sự kiện xảy ra. Nếu isterminal=1, sự kiện là “terminal” và quá trình mô phỏng kết thúc. Nếu isterminal=0, quá trình mô phỏng vẫn tiếp tục, nhưng ode45 thực hiện thêm một số việc nữa để đảm bảo rằng nghiệm trong vùng lân cận của sự kiện là đúng, và một trong các giá trị tính được phải rơi đúng vào lúc xảy ra sự kiện.

3. Khi gọi ode45, bạn chuyển options làm đối số thứ tư:

   ode45(@projectile, [0,10], [0, 3, 40, 30], options);
   

Phải sửa events như thế nào để nó dừng chạy khi chiều cao của quả bóng rớt xuống thấp hơn 3 m?

Tối ưu hóa

Trong Bài tập [bóng chày], bạn được yêu cầu tính góc ban đầu tối ưu sau khi bóng bị vụt. “Tối ưu” là cách nói hoa mỹ cho “tốt nhất;” nghĩa của nó lại tùy thuộc vào bài toán cụ thể. Với bài toán Green Monster—tìm ra góc đánh tối ưu cho một cú hết sân thì nghĩa của “tối ưu” không phải là hiển nhiên.

Rất dễ nghĩ rằng chọn góc mà cho tầm xa là lớn nhất (khoảng cách từ chỗ vụt đến chỗ chạm đất). Nhưng ở đây ta cần cố vượt qua một tường chắn cao 12 m, vì vậy có lẽ ta sẽ cần góc vụt sao cho tầm xa là lớn nhất khi quả bóng bay qua ngưỡng 12 m.

Mặc dù định nghĩa nào trong số trên cũng dùng được cho những mục đích nhất định, nhưng không có định nghĩa nào hoàn toàn đúng. Ở trường hợp này góc “tối ưu” là góc mà cho chiều cao lớn nhất tại điểm mà quả bóng chạm tường, vốn cách chỗ vụt 97 m.

Vì vậy bước đầu tiên của mọi bài toán tối ưu là định nghĩa xem “tối ưu” là gì. Bước thứ hai là định nghĩa khoảng giá trị mà bạn muốn dò tìm. Trong trường hợp này, khoảng giá trị cho phép là giữa 0 độ (song song với mặt đất) đến 90 độ (thẳng lên trời). Ta trông đợi góc tối ưu gần với 45 độ, nhưng không chắc rằng nó chệch với 45 độ là bao nhiêu. Để an toàn, ta có thể bắt đầu tìm trong khoảng rộng nhất có thể.

Cách đơn giản nhất để tìm giá trị tối ưu là chạy chương trình mô phỏng với một khoảng rộng các giá trị rồi lựa ra giá trị cho kết quả tốt nhất. Cách này không hiệu quả lắm, đặc biệt trong bài này khi mà việc tính tầm bay là rất tốn công.

Một thuật toán khá hơn là cách Tìm kiếm theo lát cắt vàng.

Tìm kiếm theo lát cắt vàng

Để trình bày cách tìm kiếm theo lát cắt vàng, tôi sẽ bắt đầu với một dạng đơn giản hơn mà tôi đặt tên là tìm kiếm theo lát cắt bạc. Ý tưởng cơ bản cũng giống như các phương pháp tìm nghiệm mà ta đã thấy ở Mục [tìm nghiệm]. Trong bài toán tìm nghiệm, ta có một bức tranh như sau:

Cho trước một hàm f có thể lượng giá được, và ta cần tìm một nghiệm của f; nghĩa là một giá trị của x sao cho f(x) = 0. Nếu ta có thể tìm được một giá trị, x₁, sao cho f(x₁) dương và một giá trị khác, x₂, sao cho f(x₂) âm, thì phải có một nghiệm trong khoảng giữa chúng (miễn là f liên tục). Lúc này ta nói x₁ và x₂ bao lấy nghiệm.

Thuật toán tiếp diễn với việc chọn một giá trị thứ ba, x₃, ở giữa x₁ và x₂ rồi tính y = f(x₃). Nếu y dương thì ta có thể lập ra một cặp mới, (x₃, x₂), bao lấy nghiệm. Nếu y âm thì cặp (x₁, x₃) bao lấy nghiệm. Bằng cách nào đi nữa thì khoảng bao cũng hẹp lại và ước đoán của ta về vị trí nghiệm trở nên chính xác hơn.

Đó là việc tìm nghiệm. Còn cách Tìm kiếm theo lát cắt vàng cũng tương tự, nhưng ta phải khởi đầu với ba giá trị, và bức tranh sẽ có dạng sau:

Biểu đồ này cho thấy rằng ta đã tính f tại ba vị trí, x₁, x₂ và x₃, rồi biết được x₂ cho giá trị lớn nhất. Nếu f liên tục, thì có ít nhất là một cực trị địa phương ở giữa x₁ và x₃, vì vậy ta sẽ nói rằng bộ ba (x₁, x₂, x₃) bao lấy một cực đại.

Bước tiếp theo là chọn điểm thứ tư, x₄, rồi tính f(x₄). Có hai kết quả có thể xảy ra, tùy theo là f(x₄) có lớn hơn f(x₂)’ hay không:

Nếu f(x₄) nhỏ hơn f(x₂) (hình bên trái), thì bộ ba mới (x₁, x₂, x₄) bao lấy cực đại. Nếu f(x₄) lớn hơn f(x₂) (hình bên phải), thì (x₂, x₄, x₃) bao lấy cực đại. Dù với cách nào đi nữa thì khoảng bao cũng thu hẹp lại và ước tính giá trị cực đại của x càng tốt hơn.

Phương pháp này áp dụng được cho hầu hết các giá trị của x₄, nhưng có cách chọn nhất định sẽ có hiệu quả hơn. Ở đây, tôi chọn cách chia đôi khoảng lớn hơn trong số hai khoảng (x₁, x₂) và (x₂, x₃).

Sau đây là đoạn chương trình trong MATLAB:

function res = optimize(V)
    x1 = V(1);
    x2 = V(2);
    x3 = V(3);

    fx1 = height_func(x1);
    fx2 = height_func(x2);
    fx3 = height_func(x3);

    for i=1:50
        if x3-x2 > x2-x1
            x4 = (x2+x3) / 2;
            fx4 = height_func(x4);
            if fx4 > fx2
                x1 = x2;  fx1 = fx2;
                x2 = x4;  fx2 = fx4;
            else
                x3 = x4;  fx3 = fx4;
            end
        else
            x4 = (x1+x2) / 2;
            fx4 = height_func(x4);
            if fx4 > fx2
                x3 = x2;  fx3 = fx2;
                x2 = x4;  fx2 = fx4;
            else
                x1 = x4;  fx1 = fx4;
            end
        end

        if abs(x3-x1) < 1e-2
            break
        end
    end
    res = [x1 x2 x3];
end

Biến đầu vào là một véc-tơ có chứa ba giá trị và bao lấy một cực đại; trong trường hợp này là các góc tính theo độ. optimize bắt đầu bằng việc ước lượng height_func cho ba giá trị này. Ta giả sử rằng height_func trả lại kết quả mà ta cần tối ưu hóa; trong bài toán Green Monster đó chính là chiều cao của quả bóng khi chạm tường.

Mỗi lượt lặp qua vòng for hàm sẽ chọn một giá trị của x4, tính height_func, rồi cập nhật bộ ba x1, x2 và x3 tùy theo kết quả thu được.

Sau khi cập nhật, nó tính độ dài khoảng bao nghiệm, x3-x1, rồi kiểm tra xem nó đã đủ ngắn chưa. Nếu được rồi, nó sẽ thoát khỏi vòng lặp và trả lại kết quả là bộ ba hiện thời. Trong trường hợp dở nhất, vòng lặp sẽ được thực hiện 50 lần.

Tôi gọi thuật toán này là Tìm kiếm theo lát cắt bạc vì nó tốt gần bằng cách Tìm kiếm theo lát cắt vàng. Hãy đọc Wikipedia về cách Tìm kiếm theo lát cắt vàng (http://en.wikipedia.org/wiki/Golden_section_search) rồi sửa lại mã lệnh trên để tính theo cách mới này.

Bạn có thể viết các hàm nhận vào chuôi của hàm khác, như ta đã thấy ở fzero và ode45. Chẳng hạn, handle_func nhận vào một chuôi hàm tên là func rồi gọi nó, truyền vào đối số là pi.

function res = handle_func(func)
    func(pi)
end

Bạn có thể gọi handle_func từ Command Window và truyền vào các chuôi hàm khác nhau làm đối số:

>> handle_func(@sin)

ans = 0

>> handle_func(@cos)

ans = -1

Hãy sửa lại optimize để cho nó nhận vào một chuôi hàm rồi lấy hàm này làm mục tiêu để tối ưu hóa.

Hàm fminsearch của MATLAB nhận vào một chuôi hàm và tìm cực tiểu địa phương của hàm này. Hãy đọc lời hướng dẫn cách dùng fminsearch rồi dùng nó để tìm góc đánh tối ưu của quả bóng chày ứng với một vận tốc cho trước.

Ánh xạ rời rạc và liên tục

Khi bạn giải PVT theo cách giải tích, kết qủa thu được là một hàm, và bạn có thể coi rằng đó là một phép ánh xạ liên tục. Khi bạn dùng một hàm giải PVT, bạn thu được hai véc-tơ (hoặc một véc-tơ và một ma trận), mà bạn có thể coi là một phép ánh xạ rời rạc.

Chẳng hạn, ở Mục [ode45], ta đã dùng hàm tốc độ sau để ước tính số con chuột như một hàm theo thời gian:

function res = rats(t, y)
    a = 0.01;
    omega = 2 * pi / 365;
    res = a * y * (1 + sin(omega * t));
end

Kết quả thu được từ ode45 là hai véc-tơ:

>> [T, Y] = ode45(@rats, [0, 365], 2);

T chứa các giá trị thời gian tại đó cần ước tính số chuột bằng ode45; Y chứa các giá trị ước tính.

Bây giờ ta hình dung như cần biết số chuột vào ngày thứ 180 của năm. Ta có thể tìm giá trị 180 trong T:

>> find(T==180)

ans = Empty matrix: 0-by-1

Nhưng không có bảo đảm gì rằng sẽ tồn tại một giá trị như vậy trong T. Ta có thể tìm ra chỉ số mà tại đó giá trị của T cắt qua 180:

>> I = find(T>180); I(1)

ans = 23

I sẽ nhận được tất cả những chỉ số ứng với các phần tử của T mà lớn hơn 180, vì vậy I(1) chính là chỉ số đầu tiên.

Sau đó ta tìm giá trị tương ứng của Y:

>> [T(23), Y(23)]

ans = 184.3451   40.3742

Cách này cho ta một ước tính thô sơ về số chuột vào ngày 180. Nếu cần tính kĩ hơn, ta có thể tìm thêm giá trị ngay trước ngày 180:

>> [T(22), Y(22)]

ans = 175.2201   36.6973

Như vậy số chuột vào ngày 180 sẽ nằm giữa 36.6973 và 40.3742.

Nhưng cụ thể con số nào trong khoảng này sẽ là ước lượng chính xác nhất? Một cách làm đơn giản là chọn ngay giá trị nào tương ứng thời gian gần với 180 hơn. Trong bài này, cách làm như vậy không hay vì giá trị thời gian mà ta cần xác định lại nằm ngay giữa.

Nội suy

Một cách làm hay hơn là vạch một đường thẳng nối hai điểm bao ngày 180 và dùng đường thẳng đó để ước tính giá trị giữa chúng. Quá trình này được gọi là nội suy tuyến tính, và MATLAB cung cấp một hàm có tên interp1 để đảm nhiệm việc này:

>> pop = interp1(T, Y, 180)

pop = 38.6233

Hai đối số đầu dùng để chỉ một phép ánh xạ rởi rạc từ các giá trị có trong T đến các giá trị trong Y. Đối số thứ ba là giá trị thời gian mà ta cần nội suy. Kết quả thu được giống như ta trông đợi, chừng ở chính giữa hai giá trị đầu khoảng bao.

interp1 cũng có thể nhận một đối số thứ tư để chỉ dạng nội suy mà bạn muốn. Mặc định là 'linear', vốn thực hiện nội suy tuyến tính. Các lựa chọn khác bao gồm 'spline' vốn dùng một đường cong spline để lượn qua bốn điểm, với hai điểm mỗi phía, và 'cubic', vốn dùng phép nội suy Hermit bậc ba cho từng đoạn.

>> pop = interp1(T, Y, 180, 'spline')

pop = 38.6486

>> pop = interp1(T, Y, 180, 'cubic')

pop = 38.6491

Ở trường hợp này ta trông đợi các kết quả thu được từ nội suy spline và bậc ba tốt hơn so với tuyến tính, vì chúng dùng nhiều điểm số liệu hơn, và ta biết rằng hàm không phải tuyến tính. Nhưng chẳng có lý do gì để ta có thể trông đợi rằng hàm spline cho kết quả tốt hơn hàm bậc ba, hoặc ngược lại. Thật may là, các kết quả không quá khác biệt.

Ta cũng có thể dùng interp1 để kéo dài đồ thị của số chuột ra ngoài khoảng các giá trị có trong T:

>> [T(end), Y(end)]

ans = 365.0000   76.9530

>> pop = interp1(T, Y, 370, 'cubic')

pop = 80.9971

Quá trình này được gọi là ngoại suy. Với những giá trị thời gian còn gần 365, ngoại suy vẫn có vẻ hợp lý, nhưng khi càng xa về phía “tương lai,” ta càng trông đợi ở đó sự kém chính xác. Chẳng hạn, sau đây là con số ước tính được khi ngoại suy cả một năm:

>> pop = interp1(T, Y, 365*2, 'cubic')

pop = -4.8879e+03

Kết quả là hoàn toàn sai.

Nội suy hàm ngược

Ta đã dùng interp1 để tìm số chuột như một hàm của thời gian; và bằng cách đảo ngược vai trò của T và Y, ta cũng có thể nội suy thời gian như một hàm của số chuột. Chẳng hạn, ta cần biết sau bao lâu thì số chuột sẽ đạt đến 20.

>> interp1(Y, T, 20)

ans = 133.4128

Cách dùng interp1 thế này có thể dễ nhầm lẫn nếu bạn nghĩ rằng các đối số như là x và y. Có thể hay hơn là bạn hình dung chúng như tập nguồn và tập đích trong một phép ánh xạ (còn đối số thứ ba là phần tử trong tập nguồn).

Hình vẽ dưới đây cho thấy f (Y vẽ theo T) và nghịch đảo của f (T vẽ theo Y).

Ở trường hợp này ta có thể dùng interp1 theo cách nào cũng được vì f có ánh xạ đơn trị, tức là mỗi giá trị của tập đích chỉ có một giá trị từ tập nguồn có ánh xạ đến nó.

Nếu ta giảm lượng cung cấp thức ăn sao cho số chuột giảm trong thời kì mùa đông thì có thể sẽ thấy kết quả sau:

Ta vẫn dùng được interp1 để ánh xạ từ T đến Y:

>> interp1(T, Y, 260)

ans = 15.0309

Như vậy là vào ngày 260, số chuột có khoảng 15 con, nhưng nếu ta cần biết vào ngày nào số chuột có 15 con thì sẽ tồn tại hai câu trả lời: 172,44 và 260,44. Nếu thử dùng interp1, ta sẽ nhận được kết quả sai:

>> interp1(Y, T, 15)         

ans = 196.3833                % SAI

Vào ngày 196, số chuột thực tế là 16,8; vì vậy interp1 còn không đạt gần con số đó! Vấn đề là ở chỗ T với vai trò hàm số của Y là một ánh xạ đa trị; với một giá trị nào đó ở tập nguồn, có hơn một giá trị ở tập đích. Điều này làm cho interp1 không tính đúng. Tôi không thể tìm thấy tài liệu nào viết về sự hạn chế nói trên, thật hơi tệ.

Chuột đồng

Như ta đã thấy, một công dụng của nội suy là để diễn giải kết quả của một bài toán số trị; song còn một công dụng khác là để lấp đầy những khoảng trống giữa các số liệu đo rời rạc.

Chẳng hạn¹, giả sử rằng số chuột đồng bị chi phối bởi phương trình tốc độ:

g(t, y) = ay - b(t)y^1. 7

trong đó t là thời gian tính theo tháng, y là số chuột, a là một tham số đặc trưng cho tốc độ tăng số chuột trong trường hợp không hạn chế, còn b là hàm số theo thời gian, đặc trưng cho ảnh hưởng của lương thực được cấp đến tốc độ chết.

Mặc dù b xuất hiện trong phương trình như một hàm liên tục, ta có thể sẽ không biết được b(t) với mọi t. Thay vào đó, ta có thể chỉ có các số liệu đo rời rạc sau đây:

t          b(t)
-          ----
0          0.0070
1          0.0036
2          0.0011
3          0.0001
4          0.0004
5          0.0013
6          0.0028
7          0.0043
8          0.0056

Nếu dùng ode45 để giải phương trình vi phân, thì ta sẽ không thể tiến gần đến những giá trị của t khi hàm tốc độ (và do đó cả b) được lượng giá. Ta cần cung cấp một hàm cho phép xác định b bất kì lúc nào:

function res = interpolate_b(t)
    T = 0:8;
    B = [70 36 11 1 4 13 28 43 56] * 1e-4;
    res = interp1(T, B, t);
end

Nhìn bao quát, hàm này dùng một phép ánh xạ rời rạc để lập ra một ánh xạ liên tục.

Hãy viết một hàm tốc độ trong đó dùng interpolate_b để lượng giá g rồi dùng ode45 để tính số chuột đồng từ t = 0 đến t = 8 với số chuột ban đầu bằng 100 và a = 0. 9.

Sau đó sửa lại interpolate_b để dùng nội suy spline và chạy lại ode45 để xem phép nội suy ảnh hưởng nhiều đến các kết quả hay không.

Thuật ngữ

nội suy:

Ước tính giá trị của một hàm dựa vào các giá trị đã biết ở hai phía.

ngoại suy:

Ước tính giá trị của một hàm dựa vào các giá trị đã biết nhưng không bao lấy giá trì cần tìm.

ánh xạ đơn trị:

Ánh xạ trong đó mỗi giá trị của tập nguồn chiếu đến một giá trị ở tập đích.

ánh xạ đa trị:

Ánh xạ trong đó có ít nhất một giá trị ở tập nguồn chiếu đến nhiều giá trị ở tập đích.

Bài tập

Một quả bóng golf² được đánh theo cú xoáy ngược sẽ phát sinh lực nâng, vốn có thể làm tăng tầm xa, nhưng năng lượng để tạo ra độ xoáy có thể sẽ làm giảm vận tốc ban đầu. Hãy viết một chương trình mô phỏng đường bay của một quả bóng golf rồi dùng nó để tính góc đánh và sự phân chia năng lượng để tạo độ xoáy và vận tốc ban đầu (từ một nguồn năng lượng nhất định) sao cho tầm xa theo phương ngang của quả bóng đạt cực đại.

Hiện tượng quả bóng xoáy bị nâng lên là do lực Magnus³; lực này vuông góc với trục quay và đường bay. Hệ số nâng tỉ lệ với tốc độ quay, và bằng khoảng 0,1 đối với quả bóng quay được 3000 vòng/phút. Hệ số cản của quả bóng vào khoảng 0,2 khi quả bóng bay nhanh hơn 20 m/s.

---

1. Ví dụ này được trích có sửa đổi từ Gerald and Wheatley, Applied Numerical Analysis, Fourth Edition, Addison-Wesley, 1989. ↩

2. Xem http://en.wikipedia.org/wiki/Golf_ball. ↩

3. Xem http://en.wikipedia.org/wiki/Magnus_effect. ↩