Bài giảng Giới thiệu lí thuyết mô phỏng và mô hình hàng chờ

Nội dung text: Bài giảng Giới thiệu lí thuyết mô phỏng và mô hình hàng chờ

1. Chương IV GIỚI THIỆU LÍ THUYẾT MÔ PHỎNG VÀ MÔ HÌNH HÀNG CHỜ 1. MỤC ĐÍCH VÀ CÁC CÔNG CỤ CỦA MÔ PHỎNG 1.1. Khái niệm về mô phỏng ngẫu nhiên Mô phỏng (Simulation) được ứng dụng rộng rãi trong kinh tế, kĩ thuật và nhiều lĩnh vực khác. Theo Từ điển chính xác Oxford, bản 1976, "mô phỏng có nghĩa là giả cách, , làm ra vẻ như, hành động như, bắt chước giống với, mang hình thức của, giả bộ như , làm giả các điều kiện của tình huống nào đó thông qua một mô hình với mục đích huấn luyện hoặc tiện lợi". Mô phỏng được áp dụng nhằm khảo sát hành vi hay sự vận động của một hệ thống thông qua các quan hệ tương tác của các thành phần của hệ thống đó để tìm ra các giá trị phù hợp của các tham số giúp cho hệ thống hoạt động tốt hơn. Một cách tổng quát, mô phỏng (hay nói đúng hơn, phương pháp mô phỏng) hàm chứa việc áp dụng một mô hình nào đó để tạo ra các số liệu đầu ra như vậy, chứ không có nghĩa là thử nghiệm một hệ thống thực tế nào đó đang cần nghiên cứu hay khảo sát. Nếu mô hình có chứa các thành phần hay yếu tố ngẫu nhiên thì chúng ta có mô phỏng ngẫu nhiên. Mô phỏng ngẫu nhiên có thể được coi là một thí nghiệm thống kê bởi các số liệu đầu ra của mô phỏng phụ thuộc vào cách thức thực hiện mô phỏng cũng như cách thức khảo sát các quan hệ tương tác của các thành phần của hệ thống. Do đó các kết quả của mô phỏng ngẫu nhiên luôn đi kèm với các sai số thí nghiệm. Tuy nhiên, mô phỏng ngẫu nhiên khác với các thí nghiệm thông thường ở chỗ nó được tiến hành hoàn toàn trên hệ thống máy tính. Thuật ngữ “phương pháp Monte−Carlo” xuất hiện ngay từ thế chiến thứ hai khi tiến hành các mô phỏng ngẫu nhiên trong quá trình phát kiến bom nguyên tử. Ngày nay, thuật ngữ này đôi khi vẫn được dùng, như khi ta nói đến phương pháp Monte−Carlo tính tích phân chẳng hạn. Tuy nhiên, phần lớn các tài liệu chuyên ngành hiện tại đang sử dụng thuật ngữ “phương pháp mô phỏng ngẫu nhiên”. Phạm vi ứng dụng của mô phỏng ngẫu nhiên ngày càng trở nên rộng rãi. Về phương diện nghiên cứu khoa học cơ bản, mô phỏng ngẫu nhiên được áp dụng để tính tích phân nhiều lớp, giải hệ phương trình, tìm các phương án tối ưu, nghiên cứu về khuyếch tán hạt, mô phỏng các chuyển động hỗn loạn. Về phương diện thiết kế các hệ thống hợp lí giải quyết các bài toán thực tiễn, mô phỏng ngẫu nhiên được sử dụng để phân tích các bài toán công nghệ, quản lí, kinh doanh, quân sự Mô phỏng được xem xét trên hai khía cạnh: nghệ thuật và kĩ thuật, mà trong một số trường hợp rất khó phân định ranh
2. giới rạch ròi. Trong chương này chúng ta nghiên cứu mô phỏng ngẫu nhiên về phương diện một số kĩ thuật, công cụ thường được sử dụng. 1.2. Các công cụ chủ yếu của mô phỏng Nguồn ngẫu nhiên (Source of randomness) Để áp dụng mô phỏng ngẫu nhiên trước hết cần phải có được một nguồn các số ngẫu nhiên. Các số ngẫu nhiên như vậy có thể được tạo ra bởi các hàm sinh số ngẫu nhiên. Trong nhiều ngôn ngữ lập trình (như Visual C++ 6.0, hay Builder C++ 5.0, ), ta sẽ thấy có một cặp hàm dạng SRAND (seed) và RANDOM để phát sinh các số (được coi là) ngẫu nhiên. Hàm SRAND, có tham số là seed được gọi là hạt mầm ngẫu nhiên, đóng vai trò khởi tạo dãy số ngẫu nhiên. Còn hàm RANDOM là hàm sinh các số ngẫu nhiên sau khi có giá trị khởi tạo. Thông thường, các nguồn này được coi như tồn tại một cách đương nhiên. Câu hỏi đặt ra là chúng đã "đủ tốt" hay chưa? Trong giáo trình này chúng ta không đi sâu vào phân tích vấn đề trên. Một cách khái quát có thể nói rằng, các số được gọi là số ngẫu nhiên được tạo ra như vậy còn xa mới thực sự là ngẫu nhiên. Một cách chính xác hơn, chúng chỉ có thể gọi là các số giả ngẫu nhiên mà thôi. Chất lượng của nguồn ngẫu nhiên có thể ảnh hưởng rất lớn tới kết quả nghiên cứu khi sử dụng phương pháp mô phỏng ngẫu nhiên. Xét về thực chất, các số giả ngẫu nhiên là các số có tính chất tất định (deterministic), nhưng chúng có tính chất giống với một dãy các giá trị thể hiện của các biến ngẫu nhiên độc lập, có phân phối đều. Ví dụ, xét dãy số: 13, 8, 1, 2, 11, 14, 7, 12, 13, 12, 17, 2, 11, 10, 3, Dãy số này trông thì có vẻ ngẫu nhiên, nhưng thực chất là tuân theo một quy tắc (hãy phát hiện ra quy tắc này). Việc tìm kiếm các thuật giải (hay các quy tắc tất định) để phát sinh ra các số giả ngẫu nhiên đủ tốt là một lĩnh vực nghiên cứu chuyên sâu của Toán học và Tin học. Mặc dù trong thực tế, khi áp dụng mô phỏng ngẫu nhiên, người ta ít khi dùng các số ngẫu nhiên tuân theo luật phân phối xác suất đều U[0, 1) trên [0, 1), nhưng nguồn số ngẫu nhiên loại này chính là cơ sở để mô phỏng các phân phối xác suất khác (xem mục 1.3). Mô hình ngẫu nhiên Hai lí do chính cho việc áp dụng mô phỏng ngẫu nhiên là: − Tổng hợp dữ liệu theo sự phân loại nhất định. − Đưa ra các dự báo. Muốn áp dụng mô phỏng ngẫu nhiên cần phải có mô hình. Như vậy, mục đích của mô phỏng ngẫu nhiên cũng gần với mục đích của mô hình hóa (modelling). Có hai loại mô hình thường được áp dụng, đó là: mô hình cơ chế (mechanistic model) và mô hình Trường Đại học Nông nghiệp Hà Nội – Giáo trình Vận trù học 100
3. tiện dụng (convenient model). Cả hai loại này đều có thể được sử dụng để trợ giúp các công việc nghiên cứu, khảo sát nhằm gia tăng sự nhận biết và tìm kiếm tri thức, dự báo và hỗ trợ việc đưa ra các quyết định điều hành hệ thống một cách hợp lí. Để ứng dụng một mô hình, ta có các lựa chọn sau: − Tiến hành các phân tích về mặt toán học để tìm hiểu hành vi của mô hình. Vấn đề này nhiều khi trở nên rất phức tạp với các hệ phi tuyến nhiều biến, do đó chúng ta cần đặt ra thêm các giả thiết. Tuy nhiên những giả thiết "chặt chẽ quá" của toán học đôi khi trở nên "đáng nghi ngờ" trong thực tế. − Thí nghiệm với mô hình đang xem xét. Đối với các mô hình ngẫu nhiên các giá trị phản hồi (số liệu đầu ra) sẽ biến thiên. Vì vậy, với những bộ tham số khác nhau, chúng ta cần tạo ra hàng loạt các kịch bản hành vi của mô hình, để từ đó tìm ra các giá trị phù hợp của các tham số giúp cho hệ thống hoạt động tốt hơn. − Đôi khi cũng cần xem xét tới sự lựa chọn thứ ba, đó là tiếp cận lai (hybrid approach) của hai lựa chọn trên. 1.3. Mô phỏng một số phân phối xác suất Một số phân phối xác suất thường gặp Để áp dụng mô phỏng ngẫu nhiên cần biết một số kiến thức cơ bản của lí thuyết xác suất − thống kê toán học mà chúng ta sẽ nhắc lại ngay sau đây. Biến ngẫu nhiên là một khái niệm quan trọng trong lí thuyết xác suất thống kê. Một cách giản lược, biến ngẫu nhiên (random variable), còn gọi là đại lượng ngẫu nhiên, được hiểu là biến nhận giá trị tuỳ thuộc vào kết quả của phép thử (phép đo, quan sát, thí nghiệm) mà không thể đoán trước được. Biến ngẫu nhiên chia làm hai loại chính: rời rạc và liên tục. Biến rời rạc có thể nhận các giá trị từ một tập hợp (có lực lượng) hữu hạn hoặc đếm được. Biến liên tục là một khái niệm toán học về loại biến ngẫu nhiên có thể nhận các giá trị dày sát nhau trên một hoặc một số khoảng/đoạn số thực nào đó (để trình bày vấn đề đơn giản, ở đây chúng ta chỉ nói tới biến ngẫu nhiên nhận các giá trị là số thực). Trong thực tế, không có một đại lượng ngẫu nhiên nào là liên tục theo nghĩa tuyệt đối, chẳng qua là chúng ta không nhận biết được (một cách cố ý hay không cố ý) khoảng cách giữa các giá trị rất sát nhau của nó mà thôi. Phân phối xác suất của biến ngẫu nhiên rời rạc được minh hoạ qua ví dụ sau: Xét biến X có thể rơi vào một trong ba trạng thái được định lượng bởi các giá trị 6, 9, 12 với các xác suất tương ứng của các trạng thái là 0,3, 0,4 và 0,3. Chú ý rằng tổng các xác suất bằng 1 (100%) được phân phối vào các giá trị biến ngẫu nhiên X có thể lấy như trình bày trong bảng sau đây, được gọi là bảng phân phối xác suất. Các giá trị của X: xi 6 9 12 Xác suất tương ứng: pi 0,3 0,4 0,3
4. 3 Cần chú ý rằng: ∑ pi = 0,3 + 0,4 +0,3 = 1. i1= Một số phân phối xác suất thường dùng của biến ngẫu nhiên liên tục và rời rạc được liệt kê dưới đây. Phân phối đều trong [0,1): X nhận các giá trị thuộc nửa khoảng [0,1) với khả năng “như nhau”. Hàm mật độ xác suất f(x) của nó được biển diễn trên hình IV.1. f(x) 1 Hình IV.1. Đồ thị hàm mật độ phân phối đều Phân phối Poát−xông: Với một hệ thống hàng chờ một kênh (xem mục 3), số lượng X tín hiệu đến trong một khoảng thời gian là một biến ngẫu nhiên, X có thể nhận các giá trị nguyên không âm 0, 1, , k, Giả sử số tín hiệu đến trung bình trong một khoảng thời gian đã biết được (kí hiệu số đó là λ) thì với một số điều kiện nhất định có thể coi X tuân theo luật phân phối xác suất Poát−xông (Poisson) như sau: Các giá trị của X: x i 0 1 k +∞ Xác suất pi P(X = 0) P(X = 0) λke−λ tương ứng P(X = k) = k! +∞ 012 k −λ⎡λλλ λ ⎤ −λ λ Dễ thấy: ∑ pi = e⎢ +++++ ⎥ = e ×= e 1. i0= ⎣ 0! 1! 2! k! ⎦ Chú ý rằng số đặc trưng cho giá trị trung bình của biến ngẫu nhiên X được gọi là kì vọng. Trong phân phối Poát−xông, kì vọng của X là λ. Số đặc trưng cho độ phân tán các giá trị của X xung quanh giá trị kì vọng của nó được gọi là độ lệch chuẩn σ. Với phân phối Poát−xông thì σ2 = λ. Trường Đại học Nông nghiệp Hà Nội – Giáo trình Vận trù học 102
5. Phân phối mũ: Trên đây ta đã xét phân phối Poát−xông của số các tín hiệu đến trong một đơn vị thời gian. Một kiểu biến ngẫu nhiên thường xét là khoảng thời gian giữa hai tín hiệu liên tiếp sẽ tuân theo phân phối mũ. Đây là biến ngẫu nhiên liên tục chỉ nhận các giá trị không âm với hàm mật độ xác suất làf(τ=λ ) e−λτ . Kí hiệu biến ngẫu t nhiên đang xét là τ thì xác suất P(τ ≤ t) = ∫ λed−λτ τ có thể hiểu là xác suất cộng dồn 0 tt cho tới t. Do đó hàm phân phối xác suất của τ là: F(t) = f(ττ=λ )d e−λτ d τ=− e −λτ t = ∫∫ 0 00 1 − e−λt. Phân phối chuẩn tắc N(0, 1): Giả sử X là biến ngẫu nhiên có phân phối chuẩn tắc N(0,1). Lúc đó nó có kì vọng m = 0 và độ lệch chuẩn σ = 1. Hàm phân phối xác suất của X có dạng: x x F(x) = P (X≤ x) = ∫ f(x)dx = ∫ (1/ 2π− )exp( x2 / 2)dx . −∞ −∞ Cho X là biến ngẫu nhiên tuân theo luật phân phối chuẩn N(m, σ2) có kì vọng m, độ Xm− lệch chuẩn σ. Lúc đó, thực hiện phép đổi biến Z = thì Z là một biến ngẫu nhiên σ tuân theo luật phân phối chuẩn tắc N(0,1). Mô phỏng các phân phối xác suất Ví dụ 1: Mô phỏng phân phối đều trên [0, 1) Cách 1: Dùng bảng số ngẫu nhiên (xem phụ lục 2A và 2B). Đây là các bảng số ghi lại các số (giả) ngẫu nhiên được phát sinh nhờ các hàm sinh số ngẫu nhiên trong máy tính. Chẳng hạn, sử dụng phụ lục 2B chúng ta nhận được một dãy số ngẫu nhiên: 0,10; 0,09; 0,73; 0,25 Cách 2: Sử dụng các hàm sinh số ngẫu nhiên (Random number generator) đã được cài đặt trên máy tính. Dù dùng bảng số ngẫu nhiên hay sử dụng các hàm sinh số ngẫu nhiên trong máy tính, ta cũng lấy ra hoặc tính được liên tiếp các số ngẫu nhiên xi trong [0, 1) với i = 1, 2, , n. Tần số các giá trị này rơi vào k khoảng nhỏ với độ dài bằng nhau 1/k được chia ra từ [0, 1) là gần như nhau (≈ n/k). Với n lớn thì các tần số đó càng sát gần n/k. Vì vậy ta coi các giá trị phát sinh được là các thể hiện của biến ngẫu nhiên X tuân theo phân phối đều trên [0, 1). Trong trường hợp cần mô phỏng biến Y phân phối đều trên [a, b), ta chỉ việc tính yi = a + (b − a)xi. Chú ý rằng để phát sinh các số ngẫu nhiên nhận giá trị nguyên 0, 1, 2, , N, chỉ cần áp dụng công thức yi = [(N + 1)xi], trong đó vế phải là phần nguyên của (N + 1)xi. Một số bảng số ngẫu nhiên nguyên hay hàm sinh số ngẫu nhiên nguyên cài đặt sẵn trong các hệ máy tính cũng giúp giải quyết vấn đề này.
6. Ví dụ 2: Mô phỏng phân phối rời rạc với luật phân phối xác suất sau Các giá trị của X: xi 6 9 12 Xác suất pi 0,3 0,4 0,3 Muốn mô phỏng phân phối trên, trước hết cần tạo ra một dãy các chữ số ngẫu nhiên bằng cách tra bảng số ngẫu nhiên hay dùng hàm sinh số ngẫu nhiên đã được cài đặt trong máy tính. Chẳng hạn ta có thể chọn dãy sau 1009732533 7652013586 3467354876 lấy từ hàng đầu bảng số ngẫu nhiên trong phụ lục 2B. Ta quy định nếu các chữ số 0, 1, 2 xuất hiện thì coi X = 6, nếu 3, 4, 5, 6 xuất hiện thì coi X = 9, còn nếu có 7, 8, 9 xuất hiện thì coi X = 12. Lúc đó ứng với 10 chữ số đầu tiên của dãy trên a1a2 a10 = 1009732533 ta có bảng sau đây cho biết các giá trị của X có thể lấy: ai 1 0 0 9 7 3 2 5 3 3 Các giá trị của X: xi 6 6 6 12 12 9 6 9 9 9 Như vậy, đã có 10 giá trị (thể hiện) của X được tạo ra. Tương tự, có thể tạo ra các thể hiện khác của X. Do tần suất (hay xác suất thực nghiệm) của mỗi chữ số ngẫu nhiên từ 0 tới 9 trong bảng số ngẫu nhiên là khoảng 10% nên tần suất (xác suất thực nghiệm) X nhận giá trị 6, 9 và 12 theo thứ tự là 30%, 40% và 30%. Do đó có thể coi P(X = 6) = 30%, P(X = 9) = 40%, P(X = 12) = 30%. Vậy muốn mô phỏng phân phối của X phải phát sinh ra một loạt các giá trị (các thể hiện) xi của biến ngẫu nhiên X tuân theo quy luật phân phối đã cho. Ví dụ 3: Mô phỏng phân phối mũ. Giả sử biến ngẫu nhiên τ tuân theo phân phối mũ với hàm phân phối xác suất là F(t) = P(τ ≤ t) = 1e− −λt , với λ là tham số đã cho của phân phối mũ. F(t) chính là xác suất để τ nhận giá trị không lớn hơn một số t cho trước. Nếu r là biến ngẫu nhiên có phân phối đều trên [0, 1) thì P(r ≥ e−λt ) = 1 − e−λt = 1 P(τ ≤ t) (xem hình III.1). Do đó, P(lnr ≥ − λt) = P(− lnr ≤ t) = P(τ ≤ t). Vậy để phát λ sinh ra giá trị ngẫu nhiên τ của τ thì trước hết cần phát sinh ra giá trị ngẫu nhiên r của r 1 và tính τ = − lnr. Chẳng hạn, từ bảng số ngẫu nhiên (phụ lục 2B), nếu lấy λ r = 0,10 và λ = 5 thì τ = −0,2×lnr = −0,2×ln0,1 = 0,46. Tiếp theo, nếu lấy r = 0,09 thì τ = − 0,2×ln 0,09 = 0,482. Cứ như vậy ta thu được một dãy các thể hiện của τ. Trường Đại học Nông nghiệp Hà Nội – Giáo trình Vận trù học 104
7. 2. ÁP DỤNG MÔ PHỎNG NGẪU NHIÊN 2.1. Vai trò của phương pháp mô phỏng Nhiều bài toán thực tế chứa các yếu tố ngẫu nhiên, bất ổn định không giải được bằng các phương pháp giải tích. Nếu chúng ta áp dụng các phương pháp giải tích thì trong nhiều trường hợp buộc phải công nhận những giả thiết chặt chẽ không được thoả mãn trên thực tế và do đó lời giải tìm được cũng ít có giá trị thực tiễn. Phương pháp mô phỏng được dùng rộng rãi để giải các bài toán loại đó, nhất là những bài toán liên quan đến hệ thống lớn, bất ổn định, hàm chứa nhiều yếu tố ngẫu nhiên. Chúng ta cần áp dụng phương pháp mô phỏng trong các tình huống sau đây: − Khi không tìm được mô hình giải tích nào thích hợp. − Các hoạt động của hệ thống thường bị ngắt quãng, đứt đoạn không theo quy luật nào cả. − Mô phỏng là phương pháp duy nhất cho chi phí tiết kiệm và tốn ít thời gian. Tuy nhiên phương pháp mô phỏng có một số điểm hạn chế sau: − Không đưa ra được lời giải chính xác. − Khó xác định được sai số. − Mô phỏng chỉ sử dụng khi môi trường có tính bất ổn định. − Mô phỏng chỉ tạo ra các phương án đánh giá chứ không đưa ra được kĩ thuật tìm lời giải tối ưu. − Mô phỏng đôi khi rất đắt tiền. 2.2. Các bước cần tiến hành khi áp dụng mô phỏng − Xác định vấn đề hay hệ thống cần mô phỏng. − Xác định mô hình mô phỏng. − Đo và thu thập số liệu cần thiết cho mô hình. − Chạy mô phỏng. − Phân tích kết quả mô phỏng, nếu cần thì phải sửa lại phương án đã được đánh giá qua chạy mô phỏng. − Chạy mô phỏng để kiểm chứng phương án mới. − Kiểm tra tính đúng đắn của mọi kết luận về hệ thống thực tế được rút ra sau khi chạy mô phỏng. Trên đây là các bước cần làm khi áp dụng mô phỏng ngẫu nhiên để tìm ra các phương án hợp lí cho các bài toán thực tế. Ngoài ra, mô phỏng còn được áp dụng để giải quyết nhiều vấn đề khác.
8. 2.3. Một số ví dụ về áp dụng phương pháp mô phỏng Ví dụ 1: Cần lựa chọn một trong hai chiến lược để phát triển sản phẩm, với các số liệu thu thập được cho trong ba bảng IV.1, IV.2 và IV.3. Bảng IV.1. Xác suất thời gian phát triển sản phẩm Thời gian phát triển sản Xác suất phẩm Chiến lược I Chiến lược II 6 0,2 0,4 9 0,3 0,4 12 0,5 0,2 Bảng IV.2. Chi phí lợi nhuận Chi phí/giá bán Chiến lược I Chiến lược II Chi phí cố định 600.000 1.500.000 Chi phí biến thiên/đơn vị 7,5 6,75 Giá bán/đơn vị sản phẩm 10 10 Bảng IV.3. Doanh số phụ thuộc thời gian phát triển sản phẩm Xác suất Doanh số 6 tháng 9 tháng 12 tháng 1.000.000 0,2 0,4 0,5 1.500.000 0,8 0,6 0,5 Vấn đề đặt ra là áp dụng phương pháp mô phỏng để tính lợi nhuận trung bình của từng chiến lược, sau đó kiểm tra kết quả (so sánh với kết quả lí thuyết). Như vậy có năm phân phối xác suất cần mô phỏng ứng với năm biến ngẫu nhiên: X1 − thời gian phát triển sản phẩm (theo chiến lược) I, X2 − thời gian phát triển sản phẩm II, X3 − doanh số cho thời gian 6 tháng, X4 − doanh số cho thời gian 9 tháng và X5 − doanh số cho thời gian 12 tháng. Trong ví dụ này, để trình bày đơn giản về vấn đề mô phỏng các phân phối xác suất của các biến trên, ta dùng mười số ngẫu nhiên, mỗi số gồm mười chữ số ngẫu nhiên rút ra từ bảng số ngẫu nhiên − phụ lục 2A (vì vậy các chữ số 0, 1, 2, , 9 mỗi số chiếm khoảng 10%). a1 a2 a3 a4 a5 a6 a7 a8 a9 a10 1 5 8 1 9 2 2 3 9 6 2 0 6 8 5 7 7 9 8 4 8 2 6 2 1 3 0 8 9 2 Trường Đại học Nông nghiệp Hà Nội – Giáo trình Vận trù học 106
9. 8 3 7 4 8 3 6 0 4 9 4 6 3 7 5 6 7 4 8 8 0 9 2 8 1 0 5 5 8 2 7 2 9 5 0 8 8 5 7 9 9 5 8 6 1 1 1 6 3 2 7 0 5 5 5 0 8 7 6 7 6 4 7 2 3 8 2 9 3 4 Ta quy định a1 ứng với X1, a2 ứng với X2, a6 ứng với X3, a8 ứng với X4 và a10 ứng với X5. Ngoài ra cũng quy định: ⎡01, thì X1 = 6 tháng (thời gian phát triển sản phẩm I) ⎢ thì X = 9 tháng a1 = ⎢234,, 1 ⎣⎢56789,,,, thì X1 = 12 tháng ⎡0123,, , thì X2 = 6 tháng (thời gian phát triển sản phẩm II) ⎢ a2 = ⎢4567,,, thì X2 = 9 tháng ⎢ ⎣89, thì X2 = 12 tháng 6 ⎡0,1 thì X3 = 10 (doanh số 6 tháng phát triển sản phẩm) a6 = ⎢ 6 ⎣2,3, , 9 thì X3 = 1,5.10 6 ⎡0,1,2,3 thì X4 = 10 (doanh số 9 tháng phát triển sản phẩm) a8 = ⎢ 6 ⎣4,5, , 9 thì X4 = 1,5.10 6 ⎡0,1,2,3,4 thì X5 = 10 (doanh số 12 tháng phát triển sản phẩm) a10 = ⎢ 6 ⎣5,6, , 9 thì X5 = 1,5.10 Cần nhắc lại một số công thức trong lĩnh vực quản trị kinh doanh như sau: + Lợi nhuận = (Doanh số − Điểm hoà vốn) × (Lợi nhuận/đơn vị sản phẩm) + Điểm hoà vốn = (Chi phí cố định)/(Lợi nhuận/đơn vị sản phẩm) + Lợi nhuận/đơn vị sản phẩm = (Giá bán/đơn vị sản phẩm) - (chi phí/đơn vị sản phẩm) Các tính toán mô phỏng được tổng hợp trong bảng IV4. Bảng IV.4. Kết quả tính toán mô phỏng Số ngẫu nhiên Thời gian Doanh số Lợi nhuận a1 a2 a6 a8 a10 I II I II I II 1 5 8 1 9 2 2 3 9 6 6 9 1,5.106 106 3,15.106 1,75.106 2 0 6 8 5 7 7 9 8 4 9 6 1,5.106 1,5.106 3,15.106 3,38.106
10. 8 2 6 2 1 3 0 8 9 2 12 6 106 1,5.106 1,9.106 3,38.106 8 3 7 4 8 3 6 0 4 9 12 6 1,5.106 1,5.106 3,15.106 3,38.106 4 6 3 7 5 6 7 4 8 8 9 9 1,5.106 1,5.106 3,15.106 3,38.106 0 9 2 8 1 0 5 5 8 2 6 12 106 106 1,9.106 1,75.106 7 2 9 5 0 8 8 5 7 9 12 6 1,5.106 1,5.106 3,15.106 3,38.106 9 5 8 6 1 1 1 6 3 2 12 9 106 1,5.106 1,9.106 3,38.106 7 0 5 5 5 0 8 7 6 7 12 6 1,5.106 106 3,15.106 1,75.106 6 4 7 2 3 8 2 9 3 4 12 9 106 1,5.106 1,9.106 3,38.106 600. 000 Điểm hoà vốn của chiến lược I = = 240. 000 10− 7, 5 1.500.000 Điểm hoà vốn của chiến lược II = = 461.538 10− 6,75 Bảng IV.5. So sánh lợi nhuận giữa chiến lược I và II Chiến lược I Chiến lược II Tổng lợi nhuận 26,5 × 106 28,91×106 Lợi nhuận trung bình 2,65 × 106 2,891×106 (Σ lợi nhuận/10) Cần chú ý rằng trong bảng IV.5 là kết quả tính toán khi chạy mô phỏng 10 lượt ứng với 10 số đã chọn ra. Nếu ta lấy càng nhiều số ngẫu nhiên thì độ chính xác đạt được càng cao. Vì vậy, nếu việc tính toán trên đây được lập trình và chạy trên máy tính với hàng trăm, hàng ngàn lượt thì độ chính xác sẽ rất cao. Qua các phân tích trên ta thấy, để tiến hành mô phỏng cần phải có: − Cơ sở dữ liệu (DataBase) − Cơ sở tri thức (KnowledgeBase) Kiểm tra kết quả mô phỏng trên bằng cách so sánh với kết quả lí thuyết được thực hiện như sau: Doanh số chiến lược I = 0,2×(0,2×106 + 0,8×1,5×106) + 0,3×(0,4×106 + 0,6×1,5×106) + 0,5×(0,5×106 + 0,5×1,5×106) = 1,295×106. Lợi nhuận trung bình chiến lược I = (1,295 - 0,24)×2,5×106 = 2,637×106. Kết quả tính toán mô phỏng là 2,65×106 rất sát với kết quả này. Tương tự ta tính được doanh số và lợi nhuận trung bình cho chiến lược II (2,84×106) và rút ra được kết luận về độ chính xác của tính toán mô phỏng. Trường Đại học Nông nghiệp Hà Nội – Giáo trình Vận trù học 108
11. Ví dụ 2: Tìm xác suất p để bao lồi của 4 điểm lấy bất kì trong vòng tròn đơn vị là một hình tam giác (bài toán Sylvester). Có lẽ cách đơn giản nhất để giải bài toán này là áp dụng mô phỏng ngẫu nhiên theo các bước sau đây: i) Gán cho biến đếm Counter giá trị ban đầu bằng 0. ii) Tiến hành một đợt gieo ngẫu nhiên tám số thực 0 ≤ ri ≤ 1 và 0 ≤ ϕi ≤ 2π (để gieo ϕi ta lấy số ngẫu nhiên thuộc [0, 1) gieo được nhân thêm với 2π), i = 1, 2, 3, 4. Đặt xi = risinϕi, yi = ricosϕi, ta có 4 điểm nằm trong hình tròn đơn vị. Đặt A = (x1, y1), B = (x2, y2), C = (x3, y3), D = (x4, y4). iii) Ta tính diện tích 4 tam giác ABC, ABD, ACD và BCD. Nếu ta có diện tích của một tam giác bằng tổng diện tích ba tam giác còn lại thì ta được bao lồi của bốn điểm A, B, C và D là một tam giác. Ta tăng giá trị của biến đếm Counter lên thêm 1, nếu trái lại biến đếm giữ nguyên giá trị cũ và quay về bước ii). Quá trình cứ thế tiếp diễn cho tới khi số đợt gieo đạt tới một giá trị khá lớn được chọn từ trước (chẳng hạn 10.000 đợt hoặc 20.000 đợt, hoặc 100.000 đợt). Mỗi lần biến đếm Counter sẽ có giá trị kết thúc khác nhau. Lấy tỉ số của số đó và số đợt, ta có tần suất xuất hiện của sự kiện "bao lồi của 4 điểm là tam giác". Số tần suất này theo luật số lớn là giá trị gần đúng của xác suất cần tính. Theo các tài liệu chuyên khảo, lời giải đúng của bài toán là: p = 35/(12π2) ≈ 0,29552. Rõ ràng, trong trường hợp này, ta nên áp dụng mô phỏng ngẫu nhiên để tính ra tần suất (việc dễ thực hiện), thay thế cho việc tính xác suất theo lí thuyết (việc khó thực hiện). Sau đây là văn bản chương trình máy tính với ngôn ngữ lập trình C giải bài toán Sylvester. #include #include #include #include #define PI 3.14159265358979 const double esp =4.5e−12; struct diem {double x,y;}; /* Tao bo so ngau nhien bang cach tron − shuffling */ int rd(){return rand()%10;}
12. double radm() {return rand()%100+0.1*rd()+0.001*rd()+0.0001*rd();} /* Chuong trinh chinh */ void main() { clrscr(); long int count = 0, reps; diem d[4]; double r, goc; printf("\n Provide number of repetitions:"); scanf("%ld",&reps); printf("\n reps= %ld",reps); srand(19587); /* Gieo ngau nhien 4 diem va tinh dien tich bon tam giac */ for(long int i=0;i 0&&s123>0&&s134>0&&s234>0) { Trường Đại học Nông nghiệp Hà Nội – Giáo trình Vận trù học 110
13. s123=sqrt(s123);s124=sqrt(s124); s134=sqrt(s134);s234=sqrt(s234); if(fabs(s123−(s124+s134+s234))<esp) count++; else if(fabs(s124−(s123+s134+s234))<esp) count++; else if(fabs(s134−(s123+s124+s234))<esp) count++; else if(fabs(s234−(s123+s124+s134))<esp) count++; } else count++; } printf("\n Number of repetitions = %ld",reps); printf("\n Number of successes = %ld",count); printf("\n Probability to compute= %0.9f", count*1.0/reps); getch(); } Các kết quả chạy chương trình trong bốn lần như sau: esp = 4.5e−12 esp = 4.5e−12 esp = 4.5e−12 esp = 4.5e−12 số lần lặp 10000 số lần lặp 20000 số lần lặp 100000 số lần lặp 200000 số lần thành công số lần thành công số lần thành công số lần thành công 3050 5941 29594 58993 xác suất: 0.30500 xác suất: 0.29705 xác suất: 0.29594 xác suất: 0.294965 Ví dụ 3: Giải bài toán tối ưu toàn cục sau (xem lại mục 4.2 chương I về bài toán tối ưu toàn cục phi tuyến). 246 24 f(x)=− 4x1111222 2,1x + x /3 + x x −+ 4x 4x → Min với điều kiện ràng buộc (miền ràng buộc D): ⎧−≤≤2,5 x1 2,5 ⎨ ⎩−≤≤1, 5 x2 1, 5 Bài toán trên đây có 4 cực tiểu địa phương và 2 cực tiểu toàn cục, có tên gọi là “Bài toán lưng lạc đà”.
14. Ta tìm phương án tối ưu toàn cục bằng phương pháp mô phỏng tôi luyện SA (Simulated Annealing). Phương pháp (SA) mô phỏng quá trình một vật thể rắn sau khi bị nung nóng ở nhiệt độ rất cao được để nguội từ từ về một nhiệt độ rất thấp. Lúc đó hàm năng lượng của vật thể sẽ đạt mức thấp nhất. Thuật giải SA áp dụng mô phỏng ngẫu nhiên (bằng lí thuyết xích Markov, như sẽ trình bày trong chương IV, có thể chứng minh được quá trình SA sẽ hội tụ theo xác suất về lời giải tối ưu toàn cục) như sau: Bước khởi tạo Ta xuất phát từ một phương án X bất kì ban đầu thoả điều kiện ràng buộc. Lấy nhiệt độ T = Tban đầu khá cao (Tban đầu =10000, chẳng hạn). Các bước lặp Tại mỗi mức nhiệt độ T thực hiện các bước sau: i) Chọn X’ ∈ D và thuộc một lân cận đủ nhỏ của X. ii) Xét Δf = f(X’) - f(X). Nếu Δf 0 thì chấp nhận X:= X’ với xác suất p=−Δ× exp( f /(Kb T)) , trong đó Kb là hằng số Boltzmann 23 (Kb = 1,38.10 ), T là nhiệt độ hiện thời trong quá trình nguội. Quy trình i) và ii) lặp lại một số lần L đủ lớn (chẳng hạn L = 200, 300, ). Sau đó tính mức nhiệt độ mới theo công thức T: = αT (α ≈ 1, chẳng hạn như α = 0,95 hay 0,99 ). Thuật toán dừng khi T ≤ Tcuối (Tcuối là giá trị đã chọn trước ≈ 0). Sau đây là văn bản chương trình annealing.cpp: /* Su dung ky thuat simulated annealing − mo phong toi luyen giai bai toan toi uu toan cuc co rang buoc */ #include #include #include #include /* Tinh gia tri ham so can cuc tieu hoa */ float f(float x,float y) { float fg = 4*pow(x,2)−2.1*pow(x,4)+pow(x,6)/3; fg = fg +x*y −4*pow(y,2)+4*pow(y,4); Trường Đại học Nông nghiệp Hà Nội – Giáo trình Vận trù học 112
15. return fg; } /* Kiem tra cac dieu kien rang buoc */ int constraint(float x,float y) { float fg; fg = x + 2.5; if (fg<0) goto ZERO; fg = 2.5 − x; if (fg<0) goto ZERO; fg = y + 1.5; if (fg<0) goto ZERO; fg = 1.5 − y; if (fg<0) goto ZERO; return 1; ZERO: return 0; } /* Thu tuc tim diem thay the mo phong qua trinh annealing */ void sa(float x0,float y0,int k,float T,float Tlast,float alfa,float delta) { float x1,y1,deltaf,u,p,ux,uy; int l=1; srand(27556); printf("\n starting value of the function=%f",f(x0,y0)); getch(); printf("\n obtained at x=%f,y=%f",x0,y0);getch(); do { do { ux=float(random(10000))/10000−0.5; x1=x0−delta*ux; uy=float(random(10000))/10000−0.5; y1=y0−delta*uy; if(constraint(x1,y1)==1)
16. { deltaf=f(x1,y1)−f(x0,y0); if(deltaf Tlast); printf("\n Optimal value fMin=%f",f(x0,y0)); getch(); printf("\n obtained at x=%f,y=%f",x0,y0);getch(); } /* Chuong trinh chinh tim diem toi uu */ void main() { clrscr(); float x0,y0,T,Tlast,k,alfa,delta; printf("\n Within the range −2.5 to 2.5 provide value x0="); scanf("%f",&x0); printf("\n Within the range −1.5 to 1.5 provide value y0="); scanf("%f",&y0); printf("\n Specify number of iterations at each T level="); scanf("%d",&k); printf("\n Specify a very high value for Tstart="); scanf("%f",&T); printf("\n Specify a very low value for Tlast="); scanf("%f",&Tlast); printf("\n Specify reduction coeffient alfa="); scanf("%f",&alfa); Trường Đại học Nông nghiệp Hà Nội – Giáo trình Vận trù học 114
17. printf("\n Specify neighbourhood radius delta="); scanf("%f",&delta); sa(x0,y0,k,T,Tlast,alfa,delta); getch(); } Kết quả chạy chương trình máy tính với thuật giải SA là: Hạt mầm Phương án Phương án Giá trị (Seed) ban đầu tối ưu fMin 27556 (0, 0) (−0.0898613, 0.7124848) −1.0316283 19587 (0.1, 0.1) (0.0898837, −0.7125957) −1.0316284 (với alfa = 0.997, delta = 0.01, Tban đầu = 10000, L = 500, Tcuối = 0.0001). 3. MỘT SỐ VẤN ĐỀ VỀ MÔ HÌNH HÀNG CHỜ 3.1. Một số yếu tố cơ bản của hệ thống hàng chờ Như đã biết, trong nhiều hoạt động sản xuất kinh doanh cũng như trong đời sống chúng ta áp dụng các hệ dịch vụ đám đông hay hệ phục vụ công cộng. Chúng có tên gọi chung là hệ thống hàng chờ (Waiting Line System). Chẳng hạn các xí nghiệp sửa chữa máy móc, các cửa hàng, các bến xe, bến cảng, trạm tổng đài, các hệ thống điện tử viễn thông, dịch vụ Internet, là các ví dụ về hệ thống hàng chờ. Mô hình hàng chờ Trong các hệ thống hàng chờ thường xuyên diễn ra hai quá trình: quá trình nảy sinh các yêu cầu (một yêu cầu còn được coi là một tín hiệu cần được phục vụ) và quá trình phục vụ các yêu cầu ấy. Song trong quá trình phục vụ của các hệ thống, do nhiều nguyên nhân khác nhau, thường xảy ra các tình trạng sau: Trong nhiều trường hợp, quá trình phục vụ không đáp ứng các yêu cầu và do đó dẫn đến kết quả là nhiều yêu cầu phải chờ để được phục vụ. Ngược lại, trong một số tình huống khác, khả năng phục vụ của hệ thống vượt quá số yêu cầu cần được phục vụ, với kết quả là hệ thống không sử dụng hết phương tiện phục vụ. Vì vậy bài toán đặt ra là: - Phân tích bản chất của quá trình diễn ra trong các hệ thống hàng chờ và thiết lập các mối liên hệ về lượng giữa các đặc trưng của các quá trình ấy. Điều đó có nghĩa là cần thiết lập hay lựa chọn một mô hình hàng chờ (Waiting Line Model) phản ánh được bản chất của hệ thống. - Trên cơ sở các mối liên hệ đã được xây dựng và các số liệu thu được từ hệ thống, cần tính toán, phân tích và đưa ra các quyết định nhằm tìm ra các giá trị thích hợp cho
18. các tham số điều khiển/thiết kế của hệ thống để thiết kế hay điều khiển các hoạt động của hệ thống hoạt động một cách có hiệu quả hơn. Các phương pháp giải bài toán mô hình hàng chờ Để tìm lời giải cho một mô hình hàng chờ người ta thường sử dụng hai phương pháp: phương pháp giải tích và phương pháp mô phỏng trên máy tính. Phương pháp giải tích để giải mô hình hàng chờ gồm các bước sau: Bước 1: Phân tích hệ thống, chủ yếu là phân tích bản chất của dòng yêu cầu/tín hiệu đến và các trạng thái của hệ thống. Bước 2: Thiết lập hệ phương trình trạng thái cho các xác suất trạng thái (xác suất để hệ thống ở một trạng thái nào đó tại thời điểm t). Bước 3: Giải hệ phương trình để tìm các xác suất trạng thái. Từ đó thiết lập các mối quan hệ giữa các chỉ tiêu cần phân tích. Bước 4: Tính toán, phân tích các chỉ tiêu, trên cơ sở đó đưa ra các nhận xét và các quyết định. Phương pháp giải tích thường sử dụng các giả thiết rất chặt chẽ của Toán học về các đặc trưng của hệ thống, vì vậy nó có một số hạn chế nhất định khi giải các bài toán thực tế. Trong khi đó, phương pháp mô phỏng/mô phỏng ngẫu nhiên để giải mô hình hàng chờ được áp dụng cho các bài toán dịch vụ đám đông không giải được bằng công cụ giải tích, nhất là những bài toán liên quan đến hệ thống lớn, bất ổn định, hàm chứa nhiều yếu tố ngẫu nhiên, không tuân theo các giả thiết quá chặt chẽ của Toán học. Trong nhiều trường hợp phương pháp mô phỏng cho ta tiết kiệm được thời gian và chi phí nghiên cứu. Tuy phương pháp mô phỏng chỉ tạo ra các phương án đủ tốt để đánh giá hoạt động của hệ thống chứ không đưa ra được kĩ thuật tìm lời giải tốt nhất, nó tỏ ra rất thành công khi giải quyết nhiều bài toán hàng chờ nảy sinh từ thực tiễn. Các bước cần tiến hành khi áp dụng phương pháp mô phỏng bao gồm: Bước 1: Xác định bài toán hay hệ thống hàng chờ cần mô phỏng và mô hình mô phỏng. Bước 2: Đo và thu thập số liệu cần thiết cần thiết để khảo sát thống kê các số đặc trưng/các yếu tố cơ bản của mô hình. Bước 3: Chạy mô phỏng kiểm chứng (test simulation) mô hình và so sánh kết quả kiểm chứng với các kết quả đã biết được trong thực tế. Phân tích kết quả chạy mô phỏng kiểm chứng, nếu cần thì phải sửa lại phương án đã được đánh giá qua chạy mô phỏng. Bước 4: Chạy mô phỏng để kiểm chứng phương án cuối cùng và kiểm tra tính đúng đắn của mọi kết luận về hệ thống thực tế được rút ra sau khi chạy mô phỏng. Triển khai hoạt động của hệ thống hàng chờ dựa trên phương án tìm được. Trường Đại học Nông nghiệp Hà Nội – Giáo trình Vận trù học 116
19. Từ những phân tích trên đây có thể thấy Lí thuyết hàng chờ (Waiting Line Theory) còn gọi là Lí thuyết hệ phục vụ công cộng hay Lí thuyết hệ dịch vụ đám đông là lĩnh vực rất quan trọng của Toán ứng dụng/Vận trù học. Nhiều bài toán thực tế trong các lĩnh vực hệ thống dịch vụ, kĩ thuật, đã được giải quyết thành công nhờ áp dụng phương pháp mô phỏng mô hình hàng chờ. Các yếu tố cơ bản của hệ thống hàng chờ Hệ thống hàng chờ tổng quát được minh hoạ như trên hình IV.2. Input Output dòng tín hiệu đến dòng tín hiệu ra KÊNH PHỤC VỤ hàng chờ Hình IV.2. Hệ thống hàng chờ Các yếu tố cơ bản của hệ thống hàng chờ bao gồm: Bố trí vật lí của hệ thống Hệ thống hàng chờ có một số dạng bố trí vật lí (phisical layout) như minh hoạ trên hình IV.3. Trên hình IV.3, các kênh phục vụ được hiểu là những thiết bị kĩ thuật hoặc con người hoặc những tổ hợp các thiết bị kĩ thuật và con người được tổ chức quản lí một cách thích hợp nhằm phục vụ các yêu cầu/các tín hiệu đến hệ thống. Chẳng hạn, ở các trạm điện thoại tự động, kênh phục vụ là các đường dây liên lạc cùng các thiết bị kĩ thuật khác phục vụ cho việc đàm thoại. Single Channel - Single Server (Một kênh phục vụ, một loại dịch vụ) Single Channel - Multi Server (Một kênh phục vụ, nhiều loại dịch vụ) Dịch vụ 1 Dịch vụ 2 Dịch vụ 3 Multi Channel - Single Server (Nhiều kênh phục vụ, một loại dịch vụ)
20. Multi Channel - Multi Server (Nhiều kênh phục vụ, nhiều loại dịch vụ) Dịch vụ 1 Dịch vụ 2 Hình IV.3. Các dạng hệ thống hàng chờ Nguyên tắc phục vụ Nguyên tắc phục vụ (hay nội quy) của hệ thống là cách thức nhận các yêu cầu vào các kênh phục vụ. Nguyên tắc phục vụ cho biết trường hợp nào thì yêu cầu được nhận vào phục vụ và cách thức phân bố các yêu cầu vào các kênh như thế nào. Đồng thời nguyên tắc phục vụ cũng cho biết trong trường hợp nào yêu cầu bị từ chối hoặc phải chờ và giới hạn của thời gian chờ. Một số nguyên tắc phục vụ thường được áp dụng trong các hệ thống hàng chờ là FCFS (First come first served), LCFS (Last come first served), SIRO (service in random order), có ưu tiên, không ưu tiên, Các phân phối xác suất của các dòng tín hiệu, dòng phục vụ Số tín hiệu đến trong một khoảng thời gian cũng như thời gian phục vụ từng tín hiệu nói chung là những biến ngẫu nhiên và do đó, chúng tuân theo các quy luật phân phối xác suất. Các quy luật phân phối xác suất này được thiết lập căn cứ các số liệu thực nghiệm thu thập từ các quan sát, thí nghiệm, hay từ cơ sở dữ liệu sẵn có. Đối với dòng tín hiệu đầu vào, thông thường chúng ta giả sử rằng số tín hiệu đến trong vòng một khoảng thời gian nào đó được ấn định trước (1 phút, 3 phút, 5 phút, 30 phút, ) tuân theo luật phân phối Poát−xông P(λ). Ở đây, tham số λ đặc trưng cho số tín hiệu đến (trung bình) trong khoảng thời gian trên. Ví dụ, số khách vào siêu thị (trung bình) là 100 người trong 1 giờ. Có nghĩa là, số khách vào siêu thị là biến ngẫu nhiên X có phân phối Poát−xông với λ = 100. Hoặc, với số cuộc gọi (trung bình) đến tổng đài trong vòng 1 phút là 3 (tín hiệu) thì có X ∼ P(3). Một cách chính xác hơn, trong những trường hợp trên, ta có dòng tín hiệu đến là dòng Poát-xông dừng (còn gọi là dòng tối giản) với các tính chất trên sau: − Tính không hậu quả: Một dòng tín hiệu có tính không hậu quả nếu xác suất xuất hiện một số tín hiệu nào đó trong một khoảng thời gian nhất định không phụ thuộc vào Trường Đại học Nông nghiệp Hà Nội – Giáo trình Vận trù học 118
21. việc đã có bao nhiêu tín hiệu đã xuất hiện và xuất hiện như thế nào trước khoảng thời gian đó. − Tính đơn nhất: Dòng tín hiệu có tính đơn nhất nếu xét trong khoảng thời gian khá bé thì sự kiện “có nhiều hơn một tín hiệu xuất hiện” hầu như không xảy ra. Về mặt thời gian ta có thể xem dòng tín hiệu có tính đơn nhất nếu thời điểm xuất hiện các tín hiệu không trùng nhau. − Tính dừng: Dòng tín hiệu có tính dừng nếu xác suất xuất hiện một số tín hiệu nào đó trong khoảng thời gian τ chỉ phụ thuộc vào độ dài của τ chứ không phụ thuộc vào điểm khởi đầu của τ. 3.2. Các chỉ số cần khảo sát Đối với một hệ thống hàng chờ, cần tìm cách để đánh giá được các chỉ số sau: − A (Arrival rate): cường độ dòng tín hiệu đến hay số tín hiệu đến trung bình trong một khoảng thời gian. Ví dụ: A = 6 (6 khách hàng đến trong 2 tiếng); A = 20 (20 cú điện thoại đến tổng đài trong 1 phút). − S (Service rate): cường độ phục vụ hay số tín hiệu trung bình được phục vụ trên một đơn vị thời gian. Ví dụ: S = 7 (hệ thống có thể phục vụ 7 khách trong 1 giờ); S = 25 (tổng đài phục vụ được 25 cú điện thoại trong 2 phút). − Lq (Number in queue hay Length of queue): số tín hiệu trung bình trong hàng chờ. − Ls (Number in system hay Length of system): số tín hiệu trung bình trong toàn hệ thống (như vậy Ls ≥ Lq). − Wq (Waiting time in queue): thời gian chờ trung bình trong hàng chờ của một tín hiệu. − Ws (Waiting time in system): thời gian chờ trung bình trong hệ thống của một tín hiệu. − Pw (Probability the system is busy): xác suất hệ thống bận (đang hoạt động) hay còn gọi là hệ số (chỉ số) sử dụng của toàn hệ thống (Utilization factor). 3.3. Tính toán các chỉ số Với mục đích tìm hiểu bước đầu, sau đây chúng ta chỉ xét các hệ thống hàng chờ với một loại dịch vụ. Bằng phương pháp giải tích (xem mục 3.1), có thể tìm được công thức tính toán các chỉ số với điều kiện: các giả thiết của mô hình được thỏa mãn. Mô hình một kênh phục vụ thoả mãn: số tín hiệu đến có phân phối Poát−xông, thời gian phục vụ có phân phối mũ Các công thức (I) sau đây đã được chứng minh (bằng phương pháp giải tích):
22. A2 A A Lq = ; Ls = ; Lq= .Ls ; S(S− A) SA− S A 1 A Wq = ; Ws = ; Pw = . S(S− A) SA− S Chẳng hạn với A = 3, S = 4 thì: 9 3 Lq = = 2, 25 ; Ls = = 3 ; 4(4− 3) 1 3 3 Wq = = 0,75 ; Ws = 1; Pw = = 0,75 . 4 4 Mô hình một kênh phục vụ thoả mãn: số tín hiện đến là phân phối Poát−xông, thời gian phục vụ có phân phối bất kì Các công thức (II) sau đây đã được chứng minh: A(A/S)22σ+ 2 A Lq==+ ; Ls Lq ; 2(1− A / S) S Lq 1 A Wq==+= ; Ws Wq ; Pw . ASS Trong đó σ độ lệnh chuẩn thời gian phục vụ một tín hiệu. Chú ý rằng, nếu thời gian 1 phục vụ tuân theo phân phối mũ (tức là có hàm mật độ là f(t) = Se-St ) thì σ= và cũng S là thời gian trung bình phục vụ một tín hiệu. Có thể nhắc lại rằng: ∞ 1 ∞ 1 m = ∫ tf (t)dt = và σ =−∫ (t m)2 f(t)dt =.v.v 0 S 0 S Lúc này các công thức (II) trở về (I). Mô hình một kênh phục vụ thoả mãn: số tín hiệu đến có phân phối Poát−xông, thời gian phục vụ có phân phối mũ, hàng chờ có giới hạn số tín hiệu tối đa M Các công thức (III) sau đã được chứng minh: 1A/S− P = ; Pw = 1 − P0 ; 0 1(A/S)− M1+ với P0 là xác suất không có tín hiệu nào trong hệ thống (hệ số không sử dụng); M P(A/S)PM0= là tỉ lệ % tín hiệu không được phục vụ do hệ thống “đầy”; Trường Đại học Nông nghiệp Hà Nội – Giáo trình Vận trù học 120
23. Pw− M(A /S)P A(1− P ) Ls = M ; Lq=− Ls M ; 1(A/S)− S Ls 1 Ws = ; Wq= Ws − . A(1− PM ) S Chú ý rằng nếu M = +∞ thì (III) trở về (I). Mô hình nhiều kênh phục vụ thoả mãn: tín hiệu đến có phân phối Poát−xông, thời gian phục vụ là phân phối mũ P0 − xác suất tất cả các kênh phục vụ đều không có tín hiệu, tìm được bằng cách tra phụ lục 3 dựa trên tỉ số A/kS (k số kênh phục vụ) hoặc tính trực tiếp từ công thức sau: 1 P0 = nk với kS> A; k1− 1A⎛⎞ 1A ⎛⎞ kA ∑ ⎜⎟+ ⎜⎟ n0= n!⎝⎠ S k! ⎝⎠ S kA− S 1A kS Pw= ( )k P ; k! S kS− A 0 AS(A / S)k A A Ls Lq Ls=+=− P ; Lq Ls ; Ws== ; Wq . (k−− 1)!(kS A)2 0 S S AA Một số điểm hạn chế của các mô hình hàng chờ Các mô hình hàng chờ giới thiệu ở trên là những mô hình tiện lợi nhất được áp dụng khá rộng rãi. Tuy nhiên, do các mô hình này công nhận các giả thiết “quá chặt chẽ” ít xảy ra trên thực tế, nên các chuyên gia trong lĩnh vực Toán ứng dụng/Vận trù học/Khoa học quản lí cũng đã đề xuất xem xét nhiều mô hình khác. Đó là các mô hình với các giả thiết như: số tín hiệu cần phục vụ là hữu hạn, dòng tín hiệu đến không phải kiểu Poát−xông, cường độ phục vụ phụ thuộc vào số tín hiệu trong hàng chờ và việc giải quyết những mô hình như vậy cần tới sự trợ giúp của phương pháp mô phỏng ngẫu nhiên. Ngay cả khi các giả thiết khá chặt chẽ của bốn mô hình đã nêu trong mục này (cũng như một số mô hình tương tự khác) là hợp lí thì việc các mô hình hàng chờ đưa ra các lời giải cho trạng thái vững (steady state solutions) cũng ít có ý nghĩa thực tế. Trong nhiều ứng dụng thực tiễn, các hệ thống hàng chờ không bao giờ đạt tới các trạng thái vững. Chẳng hạn, trong một hệ thống hàng chờ, cường độ tín hiệu đến trung bình thay đổi nhiều lần trong ngày không cho phép hệ thống đạt được trạng thái vững. Do đó, để giải quyết nhiều bài toán hàng chờ trong lĩnh vực dịch vụ đám đông và các lĩnh vực khác, cần áp dụng phương pháp mô phỏng để tìm ra các lời giải có tính thực tiễn cho các mô hình hàng chờ khi hệ thống không thể đạt tới trạng thái vững hoặc khi không có các mô hình lí thuyết thích hợp. 3.4. Áp dụng mô phỏng cho một số hệ thống hàng chờ
24. Ví dụ 1: Bài toán hệ dịch vụ hàng chờ ba kênh với dòng tối giản có từ chối. Cho biết: dòng tín hiệu đến là dòng Poát−xông dừng (còn gọi là dòng tối giản). Giãn cách thời gian giữa thời điểm đến của hai nhu cầu (tín hiệu) liên tiếp có phân phối mũ với tham số μ = 5, tức là có hàm mật độ f(t) = 5e−5t. Nếu tín hiệu xuất hiện mà có ít nhất một trong ba kênh không bận (kênh số 1 hoặc kênh số 2 hoặc kênh số 3 không bận) thì tín hiệu được phục vụ tại kênh không bận với số thứ tự nhỏ nhất; nếu trái lại (khi cả ba kênh đều bận) thì tín hiệu bị từ chối. Biết thời gian phục vụ mỗi nhu cầu là 0,5 phút, hãy xác định kì vọng toán số nhu cầu được phục vụ trong khoảng thời gian 4 phút. Như vậy, cần áp dụng mô hình hàng chờ MultiChannel − SingleServer System (Hệ thống nhiều kênh phục vụ - một loại dịch vụ) theo quy tắc First in first out (FIFO: Tín hiệu đến trước được phục vụ xong trước). Thời gian giữa hai tín hiệu liên tiếp có phân phối mũ với hàm mật độ xác suất f(t) = 5e−5t. Trong bài toán này (nhằm đơn giản các bước tính toán) thời gian phục vụ mỗi tín hiệu được coi là không đổi và bằng 0,5 phút. Chúng ta sẽ áp dụng mô phỏng để xác định số nhu cầu trung bình cần được phục vụ trong khoảng thời gian 4 phút như trình bày sau đây. Kí hiệu Ti là thời điểm đến của tín hiệu thứ i, Tki là thời điểm kết thúc dịch vụ của tín hiệu thứ i (nếu có), tại kênh thứ k (k = 1, 2, 3). Thời điểm đến của nhu cầu tiếp theo −5t là Ti = Ti−1 +τi với τ tuân theo luật phân phối mũ có hàm mật độ f(t) = 5e và hàm phân phối là F(t) = 1 −e−5t = P(τ ≤ t). Lúc đó T1 = 0, T11 = T1 + 0,5. Kết quả này cho biết thời điểm đến của tín hiệu thứ nhất là T1 = 0 và được kênh 1 phục vụ. Kết thúc phục vụ tín hiệu 1 là thời điểm T11 = T1 + 0,5 = 0,5. Máy đếm ghi nhận 1 đơn vị là số tín hiệu đã được phục vụ. Để tìm T2 theo công thức T2 = T1 + τ2, ta phát sinh τ2 theo cách đã biết ở mục 1.3: Trước hết, phát sinh số ngẫu nhiên r2 có 2 chữ số sau dấu phẩy 0≤ ri ≤1 (theo bảng số 1 1 ngẫu nhiên - phụ lục 2B) ta có r2 = 0,10. Sau đó tính τ2 = − ln r và T2 = T1 − ln r = 5 2 5 2 0 - 0,2ln0,1 = 0,46. Vậy tín hiệu tiếp theo phải vào kênh 2 vì kênh 1 còn đang bận. Máy đếm ghi thêm 1 đơn vị thời điểm kết thúc phục vụ tín hiệu 2 là T22 = T2 + 0,5 = 0,46 + 0,5 = 0,96. Tiếp tục phát sinh r3 = 0,09, ta có τ3 = −0,2ln 0,09 = 0,482. Do đó thời điểm đến của tín hiệu 3 là T3 = T2 + τ3 = 0,46 + 0,482 = 0,942. Lúc này kênh 1 đã được giải phóng do đã phục vụ xong tín hiệu 1, nên tín hiệu 3 được tiếp nhận vào kênh 1. Tại thời điểm kết thúc phục vụ tín hiệu 3 là T13 = T3 + 0,5 = 0,942 + 0,5 = 1,442 máy đếm lại ghi tiếp 1 đơn vị. Thực hiện tính toán tương tự, kết quả tổng hợp được ghi trong bảng IV.6. Bảng IV.6. Tính toán mô phỏng tìm số nhu cầu được phục vụ Trường Đại học Nông nghiệp Hà Nội – Giáo trình Vận trù học 122
25. Thời điểm Tki kết thúc phục Đếm số tín Thứ tự Số ngẫu τi = Thời điểm −lnri vụ tại kênh k hiệu tín hiệu nhiên ri −1/5lnri đến Ti 1 2 3 nhận bỏ 1 0 0,5 1 2 0,10 2,30 0,46 0,46 0,96 1 3 0,09 2,44 0,482 0,942 1,442 1 4 0,73 0,32 0,064 1,006 1,506 1 5 0,25 1,39 0,278 1,284 1,784 1 6 0,33 1,11 0,222 1,506 2,006 1 7 0,76 0,27 0,054 1,560 2,060 1 8 0,52 0,65 0,13 1,690 1 9 0,01 4,6 0,92 2,61 3,11 1 10 0,35 1,05 0,21 2,82 3,32 1 11 0,86 0,15 0,03 2,85 3,35 1 12 0,34 1,08 0,216 3,066 1 13 0,67 0,40 0,08 3,146 3,646 1 14 0,35 1,05 0,21 3,356 3,856 1 15 0,48 0,73 0,146 3,502 4,022 1 16 0,76 0,27 0,054 3,556 1 17 0,80 0,22 0,044 3,600 1 18 0,95 0,05 0,01 3,61 1 19 0,9 0,10 0,02 3,63 1 20 0,91 0,09 0,018 3,648 4,148 1 21 0,17 1,77 0,354 4,002 1 14 6 Phân tích kết quả tính toán ta thấy trong 20 nhu cầu đến thì chỉ có 14 nhu cầu được phục vụ. Tính toán tương tự 6 lần nữa ta có kết quả: x1 x2 x3 x4 x5 x6 14 15 14 12 13 15 Vậy số nhu cầu trung bình được hệ phục vụ trong vòng 4 phút vào khoảng x = (14 + 15 + 14 + 12 + 13 + 15)/6 = 13,83. Giải bài toán dịch vụ ba kênh có từ chối trên máy tính Một phần mềm máy tính với tên gọi MOPHONG1 phiên bản 1.0 đã được thiết kế dựa trên ngôn ngữ Builder C++ 5.0 để giải bài toán một dịch vụ với nhiều kênh phục vụ có từ chối như đã trình bày trong ví dụ 1 với các điều kiện sau: − Số kênh phục vụ n có thể lấy giá trị từ 2 tới 10. − Dòng tín hiệu đến là dòng Poát−xông, thời gian dãn cách giữa hai tín hiệu liên tiếp tuân theo phân phối mũ với hàm mật độ xác suất là f(t) = 5e−5t. − Thời gian trung bình phục vụ (xử lí) một tín hiệu là 0,5 phút.
26. − Tính số tín hiệu được phục vụ trong số m tín hiệu đến (m có thể lấy giá trị từ 10 tới 100). − Thực hiện k lần mô phỏng (k lần thử, k có thể lấy giá trị từ 4 tới 10). Mục đích của việc xây dựng phần mềm này nhằm vào: phục vụ dạy và học môn Mô phỏng ngẫu nhiên, cũng như tiếp tục nâng cấp phiên bản 1.0 để mô phỏng được hệ nhiều kênh phục vụ − nhiều loại dịch vụ có từ chối trong trường hợp tổng quát khi dòng tín hiệu đến và thời gian phục vụ tín hiệu có phân phối bất kì. Để chạy phần mềm MOPHONG1 chúng ta cần cài đặt Builder C++ 5.0 vào máy tính. Sau khi kích chuột vào biểu tượng của phần mềm, chọn File > Open Project > Look in > Mophong1 (Thư mục lưu trữ phần mềm) > Mp1.bpr. Sau đó chọn Run trên thanh công cụ để chạy phần mềm (xem hình IV.4). Chúng ta nhập số kênh phục vụ n = 3 vào ô số (kênh) dịch vụ, số tín hiệu phát sinh m = 20 vào ô số tín hiệu, số lần mô phỏng k = 6 vào ô số lần thử. Ngoài ra ta phải chọn hạt mầm là một số nguyên (đủ lớn) nhằm khởi tạo hàm sinh số ngẫu nhiên có phân phối đều trong [0, 1), chẳng hạn số 123456 để ghi vào ô hạt mầm ngẫu nhiên, nhằm từ đó mô phỏng dòng Poát−xông các tín hiệu đến. Sau đó chúng ta kích chuột vào nút Chạy để chạy chương trình. Kết quả ta thấy tổng cộng số tín hiệu đến là 120, trong đó có 83 tín hiệu được phục vụ (do đó có 37 tín hiệu bị từ chối). Vậy trung bình trong 20 tín hiệu đến có 13,833 tín hiệu được phục vụ. Kết quả này so với kết quả tính toán trên giấy cho ví dụ trên là khá sát nhau. Hình IV.4. Chạy phần mềm MOPHONG 1 Trường Đại học Nông nghiệp Hà Nội – Giáo trình Vận trù học 124
27. Chú ý: Việc tính toán dựa vào bảng số ngẫu nhiên cho trong phụ lục hoặc dựa vào hàm sinh số ngẫu nhiên trong máy tính không cho kết quả hoàn toàn giống nhau trong các lần chạy mô phỏng khác nhau. Điều này xảy ra là vì các bộ số ngẫu nhiên tạo được không giống nhau. Với các hạt mầm khác nhau thì hàm sinh số ngẫu nhiên cũng cho các bộ số ngẫu nhiên khác nhau và do đó các kết quả cuối cùng cũng không trùng nhau. Muốn kết quả mô phỏng ổn định hơn cần chọn số tín hiệu đến m đủ lớn. Ngoài ra, chúng ta có thể chọn số kênh phục vụ tuỳ ý, chẳng hạn n = 5. Kết quả chạy phần mềm với hạt mầm 123456 cho biết trong số 20 tín hiệu đến trung bình có 18,167 tín hiệu được phục vụ. Ví dụ 2: Một công ti điều hành một kho nguyên liệu để cấp phát cho các đốc công của 10 phân xưởng. Hiện tại, hai nhân viên phục vụ đang được công ti giao cho nhiệm vụ cấp phát nguyên liệu. Bộ phận quản lí công ti muốn cân nhắc liệu có nên thêm một nhân viên phục vụ nữa hay không. Rõ ràng rằng số đốc công (tín hiệu cần phục vụ) là hữu hạn, phân phối của số tín hiệu đến trong một đơn vị thời gian cũng không theo kiểu Poát−xông và thời gian phục vụ tín hiệu cũng không tuân theo luật phân phối mũ. Do đó, không thể tìm ra được lời giải giải tích thông qua một mô hình nhiều kênh phục vụ với các giả thiết như vậy. Phương pháp “duy nhất” là tìm cách áp dụng mô phỏng. Số liệu thu thập được − Quan sát trong vòng một tháng vào các ngày làm việc, mỗi ngày một giờ vào các thời điểm ngẫu nhiên, các số liệu về thời gian phục vụ một tín hiệu và tần suất tương ứng đã được thu thập (bảng IV.7). − Tính thời gian phục vụ trung bình cho một tín hiệu (tính kì vọng): 8×0,1 + 9×0,2 + 10×0,3 + 11×0,4 = 10 (phút). − Ngoài ra cũng đã khảo sát được: giãn cách thời gian trung bình giữa hai tín hiệu liên tiếp là 5 phút và số lượng tín hiện trung bình đến trong một khoảng 5 phút là một tín hiệu. Bảng IV.7. Thời gian phục vụ một tín hiệu và tần suất Thời gian phục vụ một Số lượng Tần suất/Xác suất thực tín hiệu (phút) nghiệm 8 15 0,1 9 30 0,2 10 45 0,3 11 60 0,4 Σ = 150 Σ = 1
28. Cần chú ý rằng, dòng tín hiệu đến chưa chắc tuân theo phân phối Poát−xông và thời gian phục vụ một tín hiệu không nhất thiết tuân theo phân phối mũ. Do đó, không áp dụng được công thức của mục 3.3 mà phải dùng mô phỏng để giải quyết vấn đề: cần bố trí bao nhiêu kênh phục vụ (nhân viên phục vụ) trong kho cấp phát nguyên liệu là hợp lí nhất? Như vậy mô phỏng có khả năng xử lí các tình huống, sự kiện như trong thực tế xảy ra chứ không bắt chúng tuân theo các phân phối xác suất nhất định hay theo các hành vi gò ép. Mô phỏng hệ thống hàng chờ − Mô phỏng tín hiệu đến: trung bình 5 phút có 1 tín hiệu đến. Chúng ta dùng 24 số ngẫu nhiên sau lấy ra từ bảng số ngẫu nhiên (phụ lục 2A), mỗi số gồm 10 chữ số để mô phỏng 24 khoảng 5 phút (như vậy tổng cộng là 120 phút, mỗi số dùng để mô phỏng một khoảng 5 phút từ 9h − 11h): 1581922396, 2068577984, 8262130892, 8374856049, 4637567488, 0928105582, 7295088579, 9586111652, 7055508767, 6472382984, 4112077556, 3440672486, , 5973470495. Nếu chữ số 7 xuất hiện trong số 10 chữ số đã chọn, ta coi như 1 tín hiệu đến trong khoảng thời gian tương ứng (vì trung bình trong 5 phút có 1 tín hiệu đến cũng giống như trung bình trong một số có 10 chữ số có một chữ số 7). Chẳng hạn trong khoảng 5 phút đầu không có tín hiệu nào, khoảng 5 phút thứ hai có 2 tín hiệu đến. Ta thấy số tín hiệu đến chỉ có thể là 0, 1, 2, 3, 4 tín hiệu. Thời điểm đến của tín hiệu trong mỗi khoảng 5 phút được quy định như sau tuỳ theo số tín hiệu đến trong khoảng đó (bảng IV.8). Chẳng hạn, nếu có hai tín hiệu đến thì thời điểm đến là vào đầu phút thứ nhất và đầu phút thứ ba. Bảng IV.8. Quy định thời điểm đến của tín hiệu Phút 1 2 3 4 5 Nếu 1 tín hiệu đến * Nếu 2 tín hiệu * * * Nếu 3 tín hiệu * * * * * Nếu 4 tín hiệu * − Mô phỏng thời gian phục vụ X một tín hiệu: Ta mô phỏng phân phối xác suất ở bảng IV.7 theo cách đã biết. Trước hết lấy ba số ngẫu nhiên có 10 chữ số: 9846413446, 8306646692, 0661684251 (hàng 4 từ dưới lên, phụ lục 2A). X = 8 phút nếu xuất hiện chữ số 0; X = 9 phút nếu xuất hiện chữ số 1, 2; X = 10 phút nếu xuất hiện chữ số 3, 4, 5; X = 11 phút nếu xuất hiên chữ số 6, 7, 8 hoặc 9. Bảng IV.9 tổng hợp kết quả mô phỏng số tín hiệu đến và thời gian phục vụ các tín hiệu. Bảng IV. 9. Kết quả mô phỏng số tín hiệu đến và thời gian phục vụ tín hiệu Trường Đại học Nông nghiệp Hà Nội – Giáo trình Vận trù học 126
29. Chu Số tín Số ngẫu nhiên 10 chữ số Thời gian phục vụ kì hiệu đến 1 1 5 8 1 9 2 2 3 9 6 0 2 2 0 6 8 5 7 7 9 8 4 2 11, 11 3 8 2 6 2 1 3 0 8 9 2 0 4 8 3 7 4 8 5 6 0 4 9 1 10 5 4 6 3 7 5 6 7 4 8 8 2 11, 10 6 0 9 2 8 1 0 5 5 8 2 0 7 7 2 9 5 0 8 8 5 7 9 2 9, 10 8 9 5 8 6 1 1 1 6 5 2 0 9 7 0 5 5 5 0 8 7 6 7 3 10, 10, 11 10 6 4 7 2 3 8 2 9 3 4 1 11 11 4 1 1 2 0 7 7 5 5 6 2 10, 8 12 3 4 4 0 6 7 2 4 8 6 1 11 13 1 8 8 2 4 1 2 9 6 3 0 14 0 6 8 4 0 1 2 0 0 6 0 15 0 9 3 3 1 4 7 9 1 4 1 11 16 7 4 5 7 4 7 7 4 6 8 4 10, 11, 11, 11 17 5 4 3 5 8 8 0 7 8 8 1 9 18 9 6 7 0 8 5 2 9 1 3 1 8 19 1 2 9 1 2 6 5 7 3 0 1 11 20 4 8 9 0 0 3 1 3 0 5 0 21 0 0 9 9 5 2 0 8 5 8 0 22 3 0 9 0 9 0 8 8 7 2 1 11 23 2 0 3 9 5 9 3 1 8 1 0 24 5 9 7 3 4 7 0 4 9 5 2 9,11 Để tiến hành minh hoạ quá trình tính toán mô phỏng cho số tín hiệu đến trong từng khoảng thời gian 5 phút và thời gian phục vụ mỗi tín hiệu, chúng ta quy ước các kí hiệu sau: Tín hiệu đến Chờ Được phục vụ Hình IV.5 tổng hợp kết quả tính toán mô phỏng cho hệ thống chờ hai kênh phục vụ - một dịch vụ. Theo hình IV.5 có thể thấy, trong thời gian 120 phút có 25 tín hiệu đến (25 số 7 xuất hiện ở các chu kì 2, 4, 5, 7, 9, 10, 11, 12, 15, 16, 17, 18, 19, 22, 24). Tổng thời gian đợi 213 phút, vì vậy thời gian đợi trung bình = 213/25 =8,52 phút. Chúng ta phân tích ý nghĩa kinh tế của tính toán mô phỏng như sau: Mô hình hàng chờ trên có hai kênh phục vụ (hai nhân viên phục vụ). Giả sử rằng lương 7$/1 giờ/1 nhân viên phục vụ, lương của mỗi đốc công là 12$/1 giờ/1 đốc công. Theo dữ kiện của bài toán, thời gian trung bình giữa hai lần tín hiệu đến liên tiếp là 5 phút. Như vậy, trong
30. trong 8 giờ có 96 lần đi và phải đợi 8,52 phút mỗi lần. Do đó, tổng hao hụt/mỗi ngày làm việc do thời gian đốc công lãng phí do phải đợi là: (8,52 x 12$ x 96)/60 = 163,56$. Ngoài ra, tổng chi lương/ngày cho hai nhân viên phục vụ là: 7$ × 8 × 2 người = 112$. Từ đó, tổng chi phí/ngày cho hệ thống hàng chờ trên là $ 275,56. 7 5 6 4 2 3 1 9:00 9:05 9:10 9:15 9:20 9:25 9:30 9:35 9:40 19 18 17 16 15 14 13 11 12 10 9 8 9:40 9:45 9:50 9:55 10:00 10:05 10:10 10:15 10:20 25 24 23 22 21 20 10:55 11:00 10:20 10:25 10:30 10:35 10:40 10:45 10:50 Hình IV.5. Tổng hợp kết quả mô phỏng hệ thống chờ Tương tự, nếu bố trí ba kênh phục vụ (ba nhân viên phục vụ) thì có thể tính được thời gian trung bình chờ đợi là 1,88 phút và tổng chi phí cho hệ thống ba kênh là 204$. Nếu dùng bốn kênh phục vụ thì thời gian chờ đợi trung bình là 0 phút, tổng chi phí là 224$. Vậy hệ thống hàng chờ trên nên dùng ba kênh phục vụ là tốt nhất. Nói cách khác, ban điều hành công ti nên điều thêm một nhân viên phục vụ tới làm việc tại kho cấp phát nguyên liệu nhằm giảm bớt chi phí “cơ hội”. So sánh với lời giải lí thuyết So sánh lời giải dùng mô phỏng với lời giải dựa trên lí thuyết mô hình hàng chờ (lúc đó cần giả sử rằng: số tín hiệu đến tuân theo luật Poát−xông, thời gian phục vụ một tín hiệu tuân theo phân phối mũ). Khi tra phụ lục 3 tìm P0 dựa vào tỉ số A/(kS) với k = 3 là số kênh phục vụ, ta được P0 = 0,11435 (A== 12,S 6;A /(kS) = 0,66) ; Trường Đại học Nông nghiệp Hà Nội – Giáo trình Vận trù học 128
31. AS(A / S)k A A Ls = P+= 2,915; Lq = Ls − = 0,915; (k−− 1)!(kS A)2 0 S S Lq Wq = = 0,076 (giờ) = 4,5 phút. A Trong khi đó theo kết quả mô phỏng thì thời gian chờ trung bình là 1,88 phút. Tuy nhiên, nếu chạy mô phỏng cho thời gian dài hơn (10 tới 15 giờ) thông qua chương trình trên máy tính thì có kết quả thời gian chờ trung bình là 2,63 phút tương đối ổn định. Từ các phân tích trên có thể thấy, trong các tình huống chúng ta không xử lí được các điều kiện toán học phức tạp ẩn chứa trong bài toán hàng chờ, lựa chọn duy nhất là áp dụng mô phỏng để giải bài toán. Phần mềm mô phỏng hệ thống hàng chờ Phần mềm máy tính với tên gọi MOPHONG2 phiên bản 1.0 đã được thiết kế dựa trên ngôn ngữ Builder C++ 5.0 để mô phỏng hệ thống hàng chờ như đã trình bày trong ví dụ 2 với các điều kiện sau: − Số kênh phục vụ n có thể lấy giá trị từ 2 tới 10. − Dòng tín hiệu đến theo quy luật: cứ trung bình 5 phút có một tín hiệu đến, nhưng chưa biết đây có phải là phân phối Poát−xông hay không. Ngoài ra, thời điểm đến của các tín hiệu trong vòng 5 phút cũng chưa biết rõ, nên chúng được mô phỏng căn cứ bảng III.8. − Thời gian phục vụ (xử lí) một tín hiệu tuân theo phân phối rời rạc đã biết: Thời gian xử lí một 8 9 10 11 tín hiệu (phút) Xác suất 0,1 0,2 0,3 0,4 − Khoảng thời gian mô phỏng có thể chọn k chu kì, mỗi chu kì 5 phút (k có thể chọn giá trị 12, 24, 32, 48, ). − Cần tính thời gian chờ trung bình (phút) cho mỗi tín hiệu đi đến hệ thống dịch vụ trên đây. Mục đích của việc xây dựng phần mềm này là phục vụ dạy và học môn Mô phỏng ngẫu nhiên, cũng như tiếp tục nâng cấp phiên bản 1.0 để mô phỏng được hệ nhiều kênh phục vụ − nhiều loại dịch vụ không có từ chối trong trường hợp tổng quát khi dòng tín hiệu đến và thời gian phục vụ tín hiệu có phân phối bất kì. Để chạy phần mềm MOPHONG2, cũng tương tự như chạy phần mềm MOPHONG1 chúng ta cần cài đặt Builder C++ 5.0 vào máy tính. Sau khi kích chuột vào biểu tượng của phần mềm, chọn File > Open Project > Look in > Mophong2 (Thư mục lưu trữ phần mềm) > Mp2.bpr. Sau đó chọn Run trên thanh công cụ để chạy phần mềm (xem hình IV.6).
32. Chúng ta nhập số kênh phục vụ n = 3 vào ô số (kênh) dịch vụ, khoảng thời gian mô phỏng với m = 24 chu kì 5 phút vào ô khoảng thời gian phát sinh tín hiệu. Ngoài ra, ta phải chọn hạt mầm là một số nguyên (đủ lớn) nhằm khởi tạo hàm sinh số ngẫu nhiên có phân phối đều trong [0, 1), chẳng hạn số 345678 để ghi vào ô hạt mầm ngẫu nhiên. Sau đó chúng ta kích chuột vào nút Chạy để chạy chương trình. Kết quả ta thấy tổng cộng có 27 tín hiệu đi tới hệ thống dịch vụ, tổng thời gian chờ là 288 phút, nên trung bình một tín hiệu phải chờ là 10,667 phút. Kết quả này so với kết quả tính toán trên giấy sử dụng bảng số ngẫu nhiên (phụ lục 2A) là khá sát nhau. Hình IV.6. Chạy phần mềm MOPHONG2 Cần nhắc lại rằng việc tính toán dựa vào bảng số ngẫu nhiên cho trong phụ lục 2A và dựa vào hàm sinh số ngẫu nhiên trong máy tính không cho kết quả hoàn toàn giống nhau do các bộ số ngẫu nhiên tạo được không giống nhau. Với các hạt mầm khác nhau thì hàm sinh số ngẫu nhiên cũng cho các bộ số ngẫu nhiên khác nhau và do đó các kết quả cuối cùng cũng không trùng nhau. Chẳng hạn, nếu thay hạt mầm 345678 (cho kết quả thời gian chờ trung bình 10,667 phút) bằng hạt mầm 456789 thì thời gian chờ trung bình sẽ chỉ là 4,158 phút. Nếu lấy trung bình của hai kết quả này thì được con số sát hơn với kết quả tính toán đã biết dựa trên phụ lục 2A. Từ các phân tích trên ta thấy, muốn kết quả mô phỏng ổn định hơn cần chọn số chu kì mô phỏng m đủ lớn và chạy cho nhiều hạt mầm khác nhau để tính giá trị trung bình của kết quả các lần chạy đó. Vấn đề này thực chất là vấn đề bố trí thí nghiệm mô phỏng Trường Đại học Nông nghiệp Hà Nội – Giáo trình Vận trù học 130
33. trên máy tính, cũng là một hướng trong nghiên cứu về mô phỏng ngẫu nhiên hiện nay đang được tiến hành trong nhiều lĩnh vực. BÀI TẬP CHƯƠNG IV 1 −x2 1. Sử dụng mô phỏng ngẫu nhiên để tính tích phân ∫ 2dx. 0 2 Hướng dẫn. Giá trị của tích phân chính là diện tích nằm dưới đường cong 2−x ứng với 0 ≤ x ≤ 1. Mỗi lần cần gieo ngẫu nhiên (một cách độc lập) một giá trị x ∈ [0, 1) và 2 một giá trị y ∈ [0, 1). Lần gieo được tính là thành công nếu y ≤ 2−x . 2. Xét mô hình một kênh phục vụ thoả mãn: số tín hiệu đến có phân phối Poat−xông, thời gian phục vụ có phân phối mũ. Giả sử A = 3 tín hiệu đến trung bình trong một phút, S = 4 tín hiệu được phục vụ trung bình trong một phút. Hãy kiểm tra lại 9 các kết quả sau và phân tích ý nghĩa của chúng: Lq = = 2, 25 ; Ls = 3/1 = 3 ; 4(4− 3) Wq = 3/ 4 = 0,75; Ws = 1; Pw = 3/ 4 = 0,75. Tuy nhiên, với Pw = 0,75, ta thấy 75% số tín hiệu phải chờ trước khi được phục vụ. Điều này có nghĩa là cần tiếp tục cải thiện hệ thống hàng chờ để hiệu quả phục vụ được tốt hơn. Một trong các biện pháp để đem lại hiệu quả phục vụ được tốt hơn là nâng cao tốc độ phục vụ. Hãy lập bảng tính và so sánh khi các tham số khác của hệ thống cố định, riêng S nhận các giá trị khác nhau S = 4, 6, 8 và 10. 3. Xét mô hình nhiều kênh phục vụ với các giả thiết như trong bài tập trên (A = 3 và S = 4). Hãy tính các chỉ số của mô hình khi số kênh phục vụ k = 2, 3. Từ đó phân tích ý nghĩa của các kết quả đạt được và so sánh với kết quả khi sử dụng mô hình một kênh phục vụ. 4. Xét mô hình hàng chờ một kênh phục vụ với dòng tín hiệu đến là dòng Poát−xông có cường độ A = 3 tín hiệu đến/giờ và dòng thời gian phục vụ có phân phối mũ với tham số S = 8 tín hiệu/giờ. Có thể hiểu đây là một trạm dịch vụ rửa xe ô tô với một bệ rửa, số ô tô đến sử dụng dịch vụ tuân theo phân phối Poát−xông với trung bình 3 xe/giờ và thời gian rửa mỗi xe tuân theo phân phối mũ với thời gian trung bình rửa mỗi xelà 1/8 giờ. − Đặt ρ = A/S, hãy chứng minh công thức tính xác suất (% thời gian) hệ thống có n n tín hiệu là pn = (1 − ρ)ρ , để từ đó có bảng phân phối xác suất sau đây: n 0 1 2 3 4 5 6 7 ≥ 8 pn 0.625 0.234 0.088 0.033 0.012 0.005 0.002 0.001 0
34. Hướng dẫn: Kí hiệu pn(t) là xác suất hệ thống có n tín hiệu tại thời điểm t, pn(t + ∆t) là xác suất hệ thống có n tín hiệu tại thời điểm t + ∆t. Thừa nhận rằng xác suất có một tín hiệu đến trong khoảng thời gian ∆t là A∆t, do đó xác suất không có tín hiệu đến trong khoảng thời gian ∆t là 1 − A∆t; xác suất có một tín hiệu được phục vụ trong khoảng thời gian ∆t là S∆t, do đó xác suất không có tín hiệu nào được phục vụ trong khoảng thời gian ∆t là 1 − S∆t. Lúc đó ta có công thức sau đây đúng với mọi n > 0 (hãy tham khảo thêm mục 2.4, Chương IV): pn(t + ∆t) = pn(t)(1 − A∆t)(1 − S∆t) + pn−1(t)(A∆t)(1 − S∆t) + pn+1(t)(1 − A∆t)(S∆t). Với n = 0, có: p0(t + ∆t) = p0(t)(1 − A∆t) + p1(t)(1 − A∆t)(S∆t). Từ các công thức trên bằng cách chuyển vế thích hợp, chia cả hai vế cho ∆t và cho qua giới hạn, sẽ có các công thức sau: dpn(t)/dt = Apn−1(t) + Spn+1(t) − (A + S)pn(t), đúng ∀n > 0, dp0(t)/dt = −Ap0(t) + Sp1(t). Do đó lời giải cho trạng thái vững (steady state solutions) phải thoả mãn: n Apn−1+ Spn+1 − (A + S)pn = 0, ∀n > 0, −Ap0 + Sp1 = 0 với n = 0 ⇒ pn = (1 − ρ)ρ , ∀n. ∞ − Từ đó tìm giá trị của các chỉ số thích hợp theo các công thức Ls = ∑npn , Ws = n0= Ls/A, Wq = Ws − 1/S, Lq = AWq. Hãy kiểm tra lại các kết quả trên theo công thức (I) mục 3.3 đã biết. 5. Xét các điều kiện của bài tập 3. Hãy xác định cần bố trí bao nhiêu chỗ cho xe chờ trước khi vào rửa để khi một xe tới sử dụng dịch vụ có chỗ đỗ với xác suất 90%. 6. Xét bài tập 3. Giả sử trạm xe chỉ có 4 chỗ cho xe chờ trước khi sử dụng dịch vụ và do đó nếu hàng chờ đã có 4 xe thì một xe mới đến sẽ bỏ đi tới chỗ rửa xe khác. Hãy áp dụng công thức (IV) mục 3.3 để xác định % khách hàng bị mất và thời gian trung bình rửa xong một xe tính từ lúc vào hàng chờ. 7. Các bài tập 3, 4 và 5 có thể được giải bằng phương pháp mô phỏng như thế nào? 8. Xét một quầy bán hàng ăn nhanh với các số liệu sau: giãn cách thời gian giữa thời điểm hai khách hàng liên tiếp đến quầy tuân theo phân phối đều trong khoảng từ 1 tới 5 phút, thời gian phục vụ mỗi một khách hàng là đúng 2 phút. Hãy thực hiện mô phỏng ngẫu nhiên và cho biết: hệ số sử dụng của quầy và số khách hàng trung bình trong hàng chờ. 9. Giải lại bài tập trên, biết rằng khách hàng chia thành hai loại: loại được ưu tiên và loại bình thường (khách thuộc loại được ưu tiên khi đến quầy được xếp phía trên tất cả các khách loại bình thường). Ngoài ra cho biết tỉ lệ khách hàng ưu tiên so với bình thường là 1/2. Trường Đại học Nông nghiệp Hà Nội – Giáo trình Vận trù học 132
35. 10. Khảo sát 200 xung tín hiệu qua các van điện đến điều khiển cơ cấu chấp hành, người ta thấy trung bình 2 giây có một chuỗi xung. Số liệu đã khảo sát được về thời gian xung của các chuỗi xung được cho trong bảng. Bằng phương pháp mô phỏng ngẫu nhiên (nếu có thể, lập chương trình tính toán trên máy tính) hãy xác định số van điện (tối thiểu) cần mở sao cho việc điều hành cơ cấu chấp hành được liên tục (nói cách khác, các chuỗi xung luôn được phục vụ kịp thời). Thời gian xung (giây) Số lần Tần suất 3 50 0,25 4 40 0,20 5 50 0,25 6 60 0,30 11. Một trạm bưu điện viễn thông có 13 cổng. Thời gian phục vụ mỗi khách hàng trung bình là 5 phút. Kết quả khảo sát thống kê cho biết số lượng tín hiệu khách hàng trung bình đến trong một giờ, còn kết quả thu thập phiếu thăm dò ý kiến khách hàng cho biết thời gian trung bình (số phút) một khách hàng có thể đợi trước khi được phục vụ như sau: Các giai đoạn Số tín hiệu/giờ Thời gian có thể đợi tối đa Nhu cầu cao 120 5 Nhu cầu trung bình 60 5,5 Nhu cầu thấp 30 6 Sử dụng mô phỏng ngẫu nhiên, hãy xác định quy trình tính toán tìm số cổng tối thiểu cần mở trong mỗi giai đoạn để đáp ứng được yêu cầu của khách hàng (những giả thiết nào cần đề ra để giải quyết vấn đề này). 12. Một trạm rút tiền có hai máy tự động, tại mỗi máy có hàng chờ cho tối đa 4 khách hàng (kể cả người đang sử dụng dịch vụ). Các khách hàng đến trạm tuân theo luật Poát − xông với thời gian giãn cách trung bình giữa hai lần liên tiếp khách đến là 1 phút. Nếu cả hai hàng chờ đều đã có đủ số lượng tối đa người chờ thì khách hàng tới trạm sẽ bỏ đi. Trong trường hợp hàng chờ còn chỗ, khách mới tới sẽ xếp vào hàng chờ ít người hơn, nếu hai hàng chờ cùng còn chỗ, khách mới tới vào hàng bên phải. Trong trường hợp hai hàng chờ đều còn chỗ, khách hàng có thể chuyển từ hàng dài sang hàng ngắn hơn nếu thấy hàng đó có ít nhất ít hơn hai người so với hàng đang đứng. Giả sử thời gian rút tiền tại mỗi máy đều tuân theo luật mũ với kì vọng là 1,5 phút. Hãy sử dụng mô phỏng ngẫu nhiên (và viết chương trình máy tính phù hợp) để xác định các chỉ số sau: − Thời gian trung bình một khách hàng sử dụng dịch vụ (tính từ khi vào hàng chờ cho tới khi rút được tiền). − Thời gian giãn cách trung bình giữa hai khách hàng bỏ đi.
36. − Hệ số sử dụng của mỗi máy rút tiền. − Độ dài trung bình của mỗi hàng chờ (kể cả khách hàng đang rút tiền). Trường Đại học Nông nghiệp Hà Nội – Giáo trình Vận trù học 134
37. Chương V PHÂN TÍCH MARKOV VÀ ỨNG DỤNG 1. CÁC KHÁI NIỆM CƠ BẢN VỀ XÍCH MARKOV 1.1. Một số định nghĩa Nhiều mô hình ngẫu nhiên trong Vận trù học, Kinh tế, Kĩ thuật, Dân số học, Di truyền học, dựa trên cơ sở là quá trình Markov. Đặc biệt, Tin sinh học (Bioinformatics), một lĩnh vực liên ngành của Sinh học phân tử và Tin học, hiện tại đang ứng dụng rất mạnh các vấn đề của lí thuyết các quá trình Markov. Các nhiều chuyên gia ngành Điện tử viễn thông và Cơ điện cũng rất quan tâm tới quá trình Markov nói chung, cũng như các quá trình sinh−tử hay quá trình hồi phục nói riêng. Ví dụ 1: Xét một hệ thống vật lí tiến triển theo thời gian. Tại thời điểm t = 0, hệ thống có thể rơi vào một trong ba trạng thái (hay vị trí) 1, 2 hoặc 3 một cách ngẫu nhiên. Kí hiệu X(0) là vị trí của hệ thống tại thời điểm t = 0 thì X(0) là một biến ngẫu nhiên, có thể nhận các giá trị 1 hoặc 2 hoặc 3 với các xác suất nhất định. Giả sử rằng căn cứ vào các kết quả quan sát hay nghiên cứu, chúng ta có bảng phân phối xác suất sau cho X(0): Các giá trị của X(0) 1 2 3 Xác suất tương ứng 0,2 0,5 0,3 Tại các thời điểm tiếp theo, chẳng hạn, t = 1, 2, 3, vị trí của hệ thống sẽ được mô tả bởi các biến ngẫu nhiên X(1), X(2), X(3), với các bảng phân phối xác suất tương ứng. Dựa trên ví dụ này, chúng ta xét định nghĩa sau về quá trình ngẫu nhiên. Định nghĩa 1 Xét một hệ thống (có thể là hệ thống vật lí, hệ thống sinh thái hay hệ thống dịch vụ, ) tiến triển theo thời gian. Gọi X(t) là vị trí (trạng thái) của hệ tại thời điểm t. Như vậy ứng với mỗi thời điểm t, X(t) chính là một biến ngẫu nhiên mô tả vị trí (trạng thái) của hệ thống. Quá trình {X(t)}t≥0 được gọi là một quá trình ngẫu nhiên. Tập hợp các vị trí có thể có của hệ gọi là không gian trạng thái. Không gian trạng thái được kí hiệu là S. Trong ví dụ trên, nếu giả sử rằng X(t) chỉ có thể nhận một trong ba giá trị 1, 2, 3 ∀t thì S = {1, 2, 3}. Giả sử trước thời điểm s, hệ đã ở trạng thái nào đó, còn tại thời điểm s, hệ ở trạng thái i. Chúng ta muốn đánh giá xác suất để tại thời điểm t (t >s), hệ sẽ ở trạng thái j. Nếu xác suất này chỉ phụ thuộc vào bộ bốn (s, i, t, j), tức là P[X(t) = j/X(s) = i] = p(s, i, t, j) là đúng ∀i, ∀j, ∀s, ∀t thì điều này có nghĩa là, sự tiến triển của hệ trong tương lai chỉ
38. phụ thuộc vào hiện tại (trạng thái của hệ tại thời điểm s) và hoàn toàn độc lập với quá khứ (tính không nhớ). Đó chính là tính Markov. Lúc này quá trình ngẫu nhiên X(t) được gọi là quá trình Markov. Trong ví dụ trên P[X(1) = 2/X(0) = 1] là xác suất có điều kiện của sự kiện X(1) = 2 (tại thời điểm t =1, hệ thống nằm tại vị trí 2) với điều kiện X(0) = 1 (tại thời điểm t = 0, hệ thống nằm tại vị trí 1). Nếu quá trình ngẫu nhiên có tính Markov thì xác suất này chỉ phụ thuộc vào trạng thái của hệ tại thời điểm s = 0 và hoàn toàn độc lập với các trạng thái của hệ trong quá khứ (trước thời điểm s = 0). Định nghĩa 2 Nếu không gian trạng thái S gồm một số hữu hạn hoặc vô hạn đếm được các trạng thái thì quá trình Markov X(t) được gọi là xích Markov. Lúc này, có thể kí hiệu S = {1, 2, 3, }, tức là các trạng thái được đánh số. Hơn nữa, nếu tập các giá trị t không quá đếm được (chẳng hạn, t = 0, 1, 2, ) thì ta có xích Markov với thời gian rời rạc, hay xích Markov rời rạc. Nếu t∈[0, ∞) thì ta có xích Markov với thời gian liên tục, hay xích Markov liên tục. Định nghĩa 3 Xét một xích Markov. Nếu xác suất chuyển trạng thái p(s, i, t, j) = p(s+h, i, t+h, j),∀i, ∀j, ∀s, ∀t và ∀h > 0 thì ta nói rằng xích Markov thuần nhất theo thời gian. Đây là một khái niệm mới và sẽ được giải thích ngay sau đây trong mục 1.2. Ngoài ra với mục đích tìm hiểu bước đầu, trong các mục 1.2 và 1.3 chúng ta sẽ chỉ xét xích Markov rời rạc và thuần nhất theo thời gian. Ví dụ về xích Markov liên tục sẽ được xem xét trong mục 2.4 và 2.5. 1.2. Ma trận xác suất chuyển trạng thái và phân phối dừng Trong mục này chúng ta đưa ra khái niệm về ma trận xác suất chuyển trạng thái của một xích Markov rời rạc và thuần nhất theo thời gian với không gian trạng thái gồm N phần tử. Trong trường hợp xích Markov rời rạc và thuần nhất có không gian trạng thái với số phần tử vô hạn đếm được, khái niệm về ma trận xác suất chuyển trạng thái sẽ được xây dựng một cách tương tự. Ví dụ 2: Xét lại ví dụ 1 trong mục 1.1, nhưng với một cách minh họa khác trong lĩnh vực dịch vụ. Trong một khu phố 1000 dân (khách hàng) có 3 siêu thị là A, B và C (A, B, C được coi là các vị trí 1, 2, 3 của hệ thống siêu thị này). Giả sử rằng, trong từng tháng mỗi khách hàng luôn trung thành với một siêu thị. Ngoài ra, cũng giả sử rằng trong tháng đầu số khách vào các siêu thị lần lượt là 200, 500 và 300; tức là có 20% khách hàng vào siêu thị A, 50% vào B và 30% vào C. Như vậy, có thể dự đoán rằng một khách hàng vào A với xác suất 0,2; vào B với xác suất 0,5 và vào C với xác suất 0,3. Để mô tả tình trạng phân chia thị phần trong tháng đầu (tháng 0) của hệ thống siêu thị trên, Trường Đại học Nông nghiệp Hà Nội – Giáo trình Vận trù học 136
39. chúng ta thiết lập biến ngẫu nhiên X(0) với quy tắc: nếu khách hàng mua hàng ở siêu thị A thì đặt X(0)=1, ở siêu thị B thì đặt X(0) = 2, còn ở siêu thị C thì X(0) = 3. Lúc đó, X(0) có bảng phân phối xác suất sau: Các giá trị của X(0) 1 2 3 Xác suất tương ứng 0,2 0,5 0,3 (0) (0) (0) (0) Kí hiệu P[X(0) = 1] = π1 , P[X(0) = 2] = π2 , P[X(0) = 3] = π3 thì véc tơ Π = (0) (0) (0) [π1 , π2 , π3 ] = [0,2 0,5 0,3] được gọi là véc tơ phân phối xác suất tại thời điểm t = 0 hay véc tơ phân phối ban đầu. Các thành phần của Π(0) cho biết tỉ lệ phần trăm (%) khách hàng vào các siêu thị A, B và C. Những tháng sau, ta giả sử xác suất để một người khách, đã vào mua hàng ở siêu thị A tháng trước vào lại A trong tháng sau luôn là 0,8; chuyển sang mua hàng ở B luôn là 0,1 và chuyển sang C luôn là 0,1. Xác suất để một người khách, đã vào mua hàng ở siêu thị B tháng trước chuyển sang A luôn là 0,07; vào lại B luôn là 0,9 và chuyển sang C luôn là 0,03. Còn xác suất để một người khách, đã vào siêu thị C tháng trước chuyển sang A luôn là 0,083; chuyển sang B luôn là 0,067 và vào lại C luôn là 0,85. Lúc đó các xác suất chuyển của khách hàng được cho thông qua ma trận xác suất chuyển trạng thái P (còn gọi là ma trận chuyển sau một bước): ⎡0,8 0,1 0,1 ⎤ ⎢ ⎥ P = ⎢0,07 0,9 0,03⎥ = [pij]3×3. ⎣⎢0,083 0,067 0,85⎦⎥ Chúng ta có thể vẽ được sơ đồ chuyển trạng thái như trên hình V.1. 0,1 0,8 1 0,083 3 0,85 0,1 0,03 0,07 0,067 2 0,9 Hình V.1. Sơ đồ chuyển trạng thái
40. Để mô tả tình trạng phân chia thị phần trong tháng t (t = 1, 2, 3, ) của hệ thống siêu thị trên, có thể thiết lập biến ngẫu nhiên X(t) với quy tắc tương tự như khi thiết lập X(0): nếu khách hàng mua hàng ở siêu thị A thì đặt X(t) = 1, ở siêu thị B thì đặt X(t) = 2, còn ở siêu thị C thì X(t) = 3. Vấn đề đặt ra là X(t) có bảng phân phối xác suất như thế nào. Trước hết ta đi tìm bảng phân phối xác suất cho X(1). Xét p12 = P[(X(1) = 2/X(0) = 1] = 0,1 là xác suất để một người khách, đã vào mua hàng ở siêu thị A tháng 0 chuyển sang mua hàng ở siêu thị B trong tháng 1. Ngoài ra, P[X(t+1) = 2/X(t) = 1] = 0,1 ∀t là số tự nhiên, vì theo giả thiết của bài toán thì xác suất để một người khách, đã vào mua hàng ở siêu thị A tháng trước chuyển sang mua hàng ở B luôn là 0,1. Vậy p12 được gọi là xác (1) suất chuyển sau một bước từ vị trí 1 sang vị trí 2, bởi vậy có thể dùng kí hiệu p12 để chỉ rõ đây là xác suất chuyển sau một bước. Các phần tử pij ∀i = 1, 2, 3 và ∀j = 1, 2, 3 của ma trận P có ý nghĩa tương tự. Dễ thấy rằng trong tháng 1 số khách hàng mua hàng tại siêu thị A là 200 × 0,8 + 500 × 0,07 + 300 × 0,083 = 219,9 (≈ 220); số khách hàng mua hàng tại siêu thị B là 200 × 0,1 + 500 × 0,9 + 300 × 0,067 = 490,1 (≈ 490); còn số khách hàng mua hàng tại siêu thị C sẽ là 200 × 0,1 + 500 × 0,03 + 300 × 0,85 = 290. Do tổng số khách hàng là 1000, nên X(1) có bảng phân phối xác suất sau: Các giá trị của X(1) 1 2 3 Xác suất tương ứng 0,2199 0,4901 0,2900 (1) (1) (1) (1) Vậy véc tơ phân phối xác suất tại thời điểm t = 1 là Π = [π1 , π2 , π3 ] cho biết tỉ lệ phần trăm khách hàng vào các siêu thị A, B và C trong tháng 1. Bằng phép tính ma trận cũng tìm được Π(1) như sau: ⎡0,8 0,1 0,1 ⎤ (1) (0) ⎢ ⎥ Π = Π × P=[0,2 0,5 0,3]× ⎢0,07 0,9 0,03⎥ = [0,2199 0,4901 0,2900]. ⎣⎢0,083 0,067 0,85⎦⎥ Tương tự có thể tìm được Π(2): ⎡0,8 0,1 0,1 ⎤ (2) (1) ⎢ ⎥ Π = Π × P = [0,2199 0,4901 0,2900] × ⎢0,07 0,9 0,03⎥ ⎣⎢0,083 0,067 0,85⎦⎥ = [0,234297 0,48251 0,283193]. Trường Đại học Nông nghiệp Hà Nội – Giáo trình Vận trù học 138
41. (2) Sau đây ta đi tìm ma trận xác suất chuyển trạng thái sau hai bước. Kí hiệu p12 là xác suất chuyển từ vị trí 1 sang vị trí 2 sau hai bước. Theo công thức xác suất toàn phần ta có: (2) p12 = P[X(2) = 2/X(0) = 1] = P[X(1) = 1/X(0) = 1] × P[X(2) = 2/X(1) = 1] + P[X(1) = 2/X(0) = 1] × P[X(2) = 2/X(1) = 2] + P[X(1) =3/X(0) = 1] × P[X(2) = 2/X(1) = 3] (1) (1) (1) (1) (1) (1) = p11 p12 + p12 p22 + p13 p32 3 (1) (1) = ∑ pp1k k2 = 0,8 × 0,1 + 0,1 × 0,9 + 0,1 × 0,067 = 0,1767. k1= Một cách hoàn toàn tương tự, ta có xác suất chuyển từ vị trí i sang vị trí j sau hai 3 (2) (1) (1) (1) (1) (1) (1) (1) (1) bước là pij = pi1 p1j + pi2 p2j + pi3 p3j = ∑ppik kj . Vậy ta có ma trận chuyển k1= sau hai bước là: (2) (2) (1) (1) 2 P = [pij ]3×3 = P × P =P × P= P ⎡0,8 0,1 0,1 ⎤ ⎡0,8 0,1 0,1 ⎤ ⎢ ⎥ ⎢ ⎥ = ⎢0,07 0,9 0,03⎥ × ⎢0,07 0,9 0,03⎥ . ⎣⎢0,083 0,067 0,85⎦⎥ ⎣⎢0,083 0,067 0,85⎦⎥ Dễ thấy Π(2) = Π(1)×P=Π(0)×P2. Tương tự, có thể chứng minh được Π(n+m) = Π(n) × P(m), trong đó Π(n+m) và Π(n) là các véc tơ phân phối tại các thời điểm t = m + n và t = n, còn P(m) là ma trận xác suất chuyển trạng thái sau m bước. (m) (m) Do P = [pij ]3×3 nên P[X(m) = j/X(0) = i] = P[X(n + m) = j/X(n) = i] = P[X(n’ + m) = (m) j/X(n’) = i] = pij , là xác suất chuyển từ vị trí i sang vị trí j sau m bước. Đặt n = s, t = n+m và h = n’ - n thì có ngay P[X(t) = j/X(s) = i] = P[X(t + h) = j/X(s + h) = i], hay p(s, i, t, j) = p(s + h, i, t + h, j) luôn đúng ∀s, ∀t, ∀h. Từ các phân tích trên đây và đối chiếu với các định nghĩa 1, 2 và 3 mục 1.1, ta thấy quá trình ngẫu nhiên X(t) với t = 0, 1, 2, trong ví dụ này chính là một xích Markov rời rạc và thuần nhất theo thời gian. Để khái quát hóa các khái niệm đã trình bày, chúng ta xét xích Markov rời rạc và thuần nhất theo thời gian X(t), t = 0, 1, 2, với không gian trạng thái gồm N phần tử mà ta kí hiệu là S = {1, 2, , N}. Định nghĩa 4 Giả sử tại thời điểm t = n, X(n) cũng có thể nhận một trong N giá trị 1, 2, , N với (n) (n) (n) (n) (n) (n) các xác suất tương ứng là π1 , π2 , , πN (với π1 + π2 + + πN = 1) thì véc tơ
42. (n) (n) (n) (n) Π = [π1 , π2 , , πN ] được gọi là véc tơ phân phối tại thời điểm t = n. Với t = 0 ta (0) (0) (0) (0) có véc tơ phân phối ban đầu Π = [π1 , π2 , , πN ]. Ma trận P = [pij]N×N, trong đó pij = p(t, i, t + 1, j) = P[X(t + 1) = j/X(t) = i] ∀t là xác suất chuyển trạng thái từ vị trí i sang vị trí j sau một bước, ∀i = 1, 2, , N và ∀j = 1, 2, , N, được gọi là ma trận xác suất chuyển trạng thái hay ma trận chuyển sau một bước. Ví dụ 3: Tiếp tục xét ví dụ trên, trong đó đã tìm được Π(1) = [0,2199 0,4901 0,2900], Π(2) = =[0,234297 0,482510 0,283193]. Dễ thấy, các véc tơ phân phối xác suất Π(1), Π(2), Π(3), tại các thời điểm t = 1, 2, 3, được tính theo công thức: Π(1) = Π(0) × P, Π(2) = Π(1) × P = Π(0) × P2 và Π(n+1) = Π(n) × P = Π(0) × Pn+1, ∀ n. Sau 21 bước (21 tháng), ta có Π(21) = [0,272257 0,455523 0,272220]. Các véc tơ phân phối (hay tỉ lệ phần trăm khách hàng vào các siêu thị A, B, C) sau 1, 2, 3, , 21 tháng được cho trong bảng V.1. (n) Vấn đề đặt ra là liệu Π = limn→∞ Π có tồn tại không và nếu tồn tại thì được tìm bằng cách nào. Trong ví dụ này, chúng ta sẽ tìm được Π = [0,273 0,454 0,273], biểu thị cho tỉ lệ phần trăm cân bằng dừng (stationary equylibrium) số khách hàng vào các siêu thị A, B, C sau một thời gian đủ dài. Bảng V.1. Tỉ lệ phần trăm khách hàng vào các siêu thị Tháng A B C 1 0,2199 0,4901 0,29 2 0,234297 0,48251 0,283193 3 0,2447183 0,476662631 0,27861905 4 0,2522664 0,472135676 0,2755979 5 0,2577373 0,46861381 0,27364893 6 0,2617056 0,465860633 0,27243373 7 0,2645868 0,463698194 0,27171505 8 0,2666806 0,461991958 0,27132742 9 0,2682041 0,460639762 0,27115613 10 0,269314 0,459563657 0,27112231 11 0,2701238 0,45870389 0,27117228 12 0,2707156 0,458014426 0,27126994 13 0,2711489 0,457459633 0,27139144 14 0,2714668 0,457011789 0,27152141 Trường Đại học Nông nghiệp Hà Nội – Giáo trình Vận trù học 140
43. 15 0,2717005 0,456649225 0,27165023 16 0,2718729 0,456354922 0,27177223 17 0,2720002 0,456115454 0,27188433 18 0,2720947 0,455920181 0,27198516 19 0,2721649 0,455760634 0,27207446 20 0,2722173 0,45563005 0,2721526 21 0,2722566 0,455523004 0,27222035 Cách tính Π Xuất phát từ Π(n+1) = Π(n) × P, cho qua giới hạn cả hai vế khi n → ∞ ta có: Π = Π × P, hay Π ×(I - P) = 0. Do P là ma trận đặc biệt (ma trận chuyển) nên nó là ma trận suy biến. Khi viết lại dưới dạng hệ phương trình (3 ẩn, 3 phương trình) ta phải loại bớt một phương trình đi và thêm vào hệ thức π1+ π2+ π3= 1 và ràng buộc πk ≥ 0 (k = 1, 2, 3). Kí hiệu x = π1, y = π2 và z = π3, ta sẽ có hệ: ⎧0,2x−− 0,07y 0,083z = 0 ⎧x= 0,273 ⎪ ⎪ ⎨−+0,1x 0,1y − 0,067z = 0 ⇔=⎨y 0,454 ⎪ ⎪ ⎩xyz1++= ⎩z= 0,273 Vậy Π = [0,273 0,454 0,273]. Định nghĩa 5 Xét xích Markov rời rạc và thuần nhất với ma trận chuyển P = [pij]N×N. Lúc đó, véc tơ phân phối xác suất Π = [π1, π2, , πN] thỏa mãn điều kiện Π ×(I - P) = 0 được gọi là phân phối dừng của xích Markov đã cho. Có thể thấy ngay, phân phối dừng Π không phụ thuộc vào Π(0) mà chỉ phụ thuộc vào ma trận P. Một cách toán học, ta nói mô hình xích Markov rời rạc thuần nhất chính là bộ ba (X(tn), S/Π, P). Áp dụng mô hình xích Markov để phân tích một vấn đề nào đó trong Kinh tế, Kĩ thuật, Sinh học, được coi là việc ứng dụng phân tích Markov. 1.3. Các tính chất và định lí Xét xích Markov rời rạc và thuần nhất với ma trận chuyển P = [pij]N×N. Có thể chứng minh được các tính chất và định lí sau: Các tính chất
44. N (n+m) (n) (m) 1/p ij = ∑ p ik p kj (đây là phương trình Chapman-Kolmogorov). k =1 2/P(2) = P × P = P2, P(n) = Pn và P(n+m) = P(n) × P(m). 3/Π(n+m) = Π(n) × P(m). Định lí 1 1/Giả sử P là ma trận xác suất chuyển chính quy, tức là tồn tại chỉ số n0, sao cho (n0 ) ∀ i, j thì xác suất chuyển từ i đến j sau n0 bước là một số dương: pij > 0. Khi đó tồn tại (n) π1, π2, , πN > 0 và π1 + π2 + + πN = 1 để cho lim p ij = πj, không phụ thuộc vào i. n→∞ Các số π1, π2, , πN được tìm từ hệ phương trình N N xjkkj==∑ x p , j 1, 2, , N , xj ≥ 0 ∀j và ∑ x1j = . k1= j1= 2/Nếu có các số π1, π2, , πN thoả mãn điều kiện π1+ π2+ + πN = 1 và (n) lim p ij = πj, không phụ thuộc vào i thì ma trận P là ma trận chính quy. n→∞ Chú ý: (n) Phân phối [π1, π2, , πN] thoả mãn điều kiện π1 + π2 + + πN = 1 và lim p ij = πj, không n→∞ phụ thuộc vào i, được gọi là phân phối giới hạn. Ngoài ra, nếu điều kiện πj > 0, ∀j được thỏa mãn thì phân phối này được gọi là phân phối Ergodic. Có thể chứng minh được rằng, nếu phân phối giới hạn tồn tại thì đó là phân phối dừng (duy nhất). Tuy nhiên, điều ngược lại không luôn đúng. 2. MỘT SỐ ỨNG DỤNG CỦA PHÂN TÍCH MARKOV Phân tích Markov có nhiều ứng dụng trong Kinh tế, Quản trị kinh doanh, Kĩ thuật, Sinh học, Xã hội học Trong mục này, chúng ta sẽ xem xét các ứng dụng như tìm cân bằng thị phần, xác định chính sách thay thế vật tư thiết bị, dự báo thất thu cho các hợp đồng thực hiện trước, tìm phân phối giới hạn của một hệ thống kĩ thuật và một ứng dụng của quá trình sinh − tử cho hệ thống hàng chờ. 2.1. Tìm cân bằng thị phần Ta nhắc lại một cách vắn tắt bài toán cho ở mục 1.2: Trong một khu phố 1000 dân (khách hàng) có 3 siêu thị là A, B và C. Giả sử, trong tháng đầu, số khách vào các siêu thị lần lượt là 200, 500 và 300. Những tháng sau đó, ta giả sử xác suất để một khách hàng (đã vào siêu thị A lúc trước) vào lại A luôn là 0,8; chuyển sang B luôn là 0,1 và chuyển sang C luôn là 0,1 Các xác suất chuyển khác của khách hàng ("trụ lại" B, chuyển sang A, chuyển sang C ) được cho thông qua ma trận chuyển P Trường Đại học Nông nghiệp Hà Nội – Giáo trình Vận trù học 142
45. ⎡0,8 0,1 0,1 ⎤ ⎢ ⎥ P = ⎢0,07 0,9 0,03⎥ ⎣⎢0,083 0,067 0,85⎦⎥ Lúc đó, theo kết quả đã biết, tỉ lệ phần trăm cân bằng dừng (khi thời gian đủ dài) số khách hàng vào các siêu thị A, B, C là 27,3%, 45,4% và 27,3% có thể tìm được từ hệ Π ×(I - P) = 0. 2.2. Chính sách thay thế vật tư thiết bị Trong một hệ thống điện kĩ thuật, các thiết bị cùng một loại được phân ra các trạng thái sau đây: vừa mới thay, còn tốt, vẫn dùng được và đã bị hỏng. Theo số liệu thống kê được, ta có ma trận xác suất chuyển trạng thái như sau: ⎡0 0,8 0,2 0 ⎤ ⎥ ⎢0 0,6 0,4 0 P= ⎢ ⎥ , ⎢0 0 0,5 0,5⎥ ⎢ ⎥ ⎣1,0 0 0 0 ⎦ trong đó, sau mỗi tuần (xem hàng đầu của ma trận P) có 0%, 80%, 20% và 0% số các thiết bị mới thay chuyển sang trạng thái mới thay, còn tốt, vẫn dùng được và đã bị hỏng. Các hàng khác của ma trận P được giải thích một cách tương tự. Ta đi tìm phân phối dừng bằng phương pháp đã biết. Xuất phát từ Π(n+1) = Π(n) × P, cho qua giới hạn cả hai vế khi n→∞ ta có: Π = Π × P, hay Π ×(I - P) = 0. Do P là ma trận đặc biệt (ma trận chuyển xác suất) nên nó là ma trận suy biến. Khi viết lại dưới dạng hệ phương trình (4 ẩn, 4 phương trình) ta phải loại bớt một phương trình đi và thêm vào hệ thức π1+ π2+ π3 + π4 = 1 và ràng buộc πk ≥ 0 (k = 1, 2, 3, 4). Kí hiệu x1 = π1, x2 = π2, x3 = π3 và x4 = π4 ta sẽ có hệ: ⎧xx14−= 0 ⎪ −+0,8x 0, 4x = 0 ⎧ 1 ⎪ 12 xx14= = ⎪ ⎪ 6 ⎨−−0, 2x123 0, 4x + 0,5x = 0 ⇔ ⎨ 1 ⎪−+=0,5x x 0 ⎪xx= = . ⎪ 34 ⎩⎪ 233 ⎩⎪xxxx11234+++= Vậy phân phối dừng Π = [1/6 1/3 1/3 1/6]. Giả sử rằng chi phí thay mới một thiết bị là 25 nghìn (đồng) và thất thu khi mỗi một thiết bị hỏng là 18,5 nghìn thì mỗi tuần hệ thống trên phải chi trung bình trên một thiết bị số tiền là: (1/6)×25 + (1/6)×18,5 = 7,25 nghìn/thiết bị/tuần.
46. Ta xét phương án thứ hai cho việc thay thế vật tư thiết bị với ma trận xác suất chuyển trạng thái sau đây: ⎡0 0,8. 0,2⎤ ⎢ ⎥ P = ⎢0 0,6 0, 4⎥ . ⎣⎢1,0 0 0 ⎦⎥ Ma trận này tương ứng với chính sách mới về thay thế vật tư thiết bị là: thay thế mỗi thiết bị một khi kiểm tra và phát hiện thiết bị ở trạng thái vẫn dùng được. Điều này có thể dẫn tới việc giảm thiểu thất thu do thiết bị hỏng gây nên. Thật vậy, ứng với ma trận P trên đây, phân phối dừng Π = [1/4 1/2 1/4]. Lúc này, mỗi tuần hệ thống trên phải chi trung bình trên một thiết bị số tiền là: (1/4)×25 + (0)×18,5 = 6,25 nghìn/thiết bị/tuần. Như vậy hệ thống sẽ tiết kiệm được 1 nghìn/thiết bị/một tuần. Nếu hệ thống có 2000 thiết bị thì nhờ chính sách thay thế vật tư mới, mỗi tuần hệ thống sẽ tiết kiệm được 2 triệu (đồng). 2.3. Phân tích Markov trong dự báo thất thu cho các hợp đồng thực hiện trước Một công ti kinh doanh trong ngành điện chuyên về sửa chữa và thay thế phụ tùng đề ra chính sách tín dụng: đáp ứng yêu cầu của khách hàng trước, thanh toán sau. Phần nhiều hợp đồng sẽ được thanh toán đúng thời hạn, một tỉ lệ nhất định sẽ được công ti cho thanh toán chậm, còn một số ít không thanh toán được. Theo kinh nghiệm, sau hai hay ba hợp đồng thanh toán chậm của một khách hàng nào đó là hợp đồng không thanh toán được sau một thời gian dài, công ti coi đây là hợp đồng “xấu” và sẽ cắt bỏ chính sách tín dụng với khách hàng đó. Như vậy tại từng thời điểm các hợp đồng có thể rơi vào một trong các trạng thái sau: − S0: hợp đồng được thanh toán, − S1: hợp đồng không được thanh toán, − S2: hợp đồng sẽ được thanh toán đúng thời hạn, − S3: hợp đồng sẽ được thanh toán chậm. Sau đây là ma trận xác suất chuyển trạng thái (sau từng tháng): ⎡1 0 0 0 ⎤ ⎢0 1 0 0 ⎥ P = ⎢ ⎥ ⎢0,5 0 0,3 0,2⎥ ⎢ ⎥ ⎣0,4 0,3 0,2 0,1⎦ Hiện tại công ti có các hợp đồng phải thanh toán đúng hạn với tổng số 500 triệu và các hợp đồng cho thanh toán chậm với tổng số 100 triệu. Hãy xác định trong tổng trên có bao nhiêu sẽ được thanh toán, còn bao nhiêu sẽ là nợ “xấu” không đòi được. Trường Đại học Nông nghiệp Hà Nội – Giáo trình Vận trù học 144
47. Đây là bài toán khá phức tạp liên quan tới phân loại các trạng thái của xích Markov là vấn đề chúng ta sẽ không trình bày trong giáo trình này. Tuy nhiên, có thể thấy ngay rằng các trạng thái S0 và S1 là các trạng thái “hấp thụ” (absorbing state), tức là mọi hợp đồng dù hiện đang ở trạng thái nào thì cuối cùng sau một thời gian nhất định cũng sẽ rơi vào một trong hai trạng thái trên. Trong khi đó các trạng thái S2 và S3 được gọi là các trạng thái truyền ứng (hay các trạng thái di chuyển). Để tìm câu trả lời cho vấn đề đặt ra, chúng ta cần thực hiện các bước sau: Trước hết, ta chia ma trận P theo khối. ⎡J Ο ⎤ ⎡1 0⎤ ⎡0,5 0 ⎤ ⎡0 0⎤ ⎡0,3 0, 2⎤ P = ⎢ ⎥ với J = ⎢ ⎥ , K= ⎢ ⎥ , O = ⎢ ⎥ , M = ⎢ ⎥ . ⎣Κ Μ⎦ ⎣0 1⎦ ⎣0,4 0,3⎦ ⎣0 0⎦ ⎣0,2 0,1⎦ Sau đó, ta tìm ma trận R = I - M và ma trận nghịch đảo của nó R−1, ở đây I là ma trận đơn vị cùng cỡ với ma trận M. Ta có: −1 ⎡1,5254 0,3390⎤ R = ⎢ ⎥ , ⎣0,3390 1,1864 ⎦ và tính được: −1 ⎡0,8983 0,1017⎤ R ×K = ⎢ ⎥ . ⎣0,6441 0,3559⎦ Các phần tử trong ma trận trên có ý nghĩa đặc biệt. Trong số các hợp đồng hiện tại ở trạng thái S2 (phải thanh toán đúng kì hạn) cuối cùng sau một thời gian nhất định có 89,83% sẽ rơi vào trạng thái S0 (được thanh toán) và 10,17% sẽ rơi vào trạng thái S1 (không dược thanh toán). Còn trong số các hợp đồng hiện tại ở trạng thái S3 (thanh toán chậm) cuối cùng sau một thời gian nhất định có 64,41% sẽ rơi vào trạng thái S0 (được thanh toán) và 35,59% sẽ rơi vào trạng thái S1 (không được thanh toán). Thực hiện phép tính: ⎡0,8983 0,1017⎤ [500 100]× ⎢ ⎥ = [459,32 140,68], ⎣0,6441 0,3559⎦ ta thấy trong 500 triệu phải thanh toán đúng kì hạn và 100 triệu thanh toán chậm cuối cùng sẽ có 459,32 triệu được thanh toán và 140,68 triệu là nợ “xấu” không đòi được. Để cải thiện tình trạng này, công ti cần nghiên cứu tìm ra một chính sách tín dụng hợp lí hơn. Ngoài ra, ma trận R−1 còn cho biết các thông tin sau: − Tổng của các phần tử trên hàng thứ nhất là 1,8644 là thời gian trung bình (tháng) mà một hợp đồng dạng phải thanh toán đúng kì hạn sẽ trải qua trước khi rơi vào một trong các trạng thái hấp thụ, tức là trở thành hợp đồng thanh toán được hoặc hợp đồng “xấu”.
48. − Tổng các phần tử trên hàng thứ hai của R−1 cũng có ý nghĩa tương tự đối với các hợp đồng dạng thanh toán chậm. − Phần tử nằm trên hàng 1 và cột 1 của R−1 cho biết thời gian trung bình (tháng) mà một hợp đồng dạng phải thanh toán đúng hạn sẽ ở trong trạng thái S2 trước khi nó rơi vào một trong các trạng thái hấp thụ là 1,5254 tháng. Phần tử nằm trên hàng 1 và cột 2 cho biết thời gian trung bình (tháng) mà một hợp đồng dạng phải thanh toán đúng hạn sẽ ở trong trạng thái S3 trước khi nó rơi vào một trong các trạng thái hấp thụ là 0,3390 tháng. − Các phần tử nằm trên hàng 2 của ma trận R-1 có ý nghĩa tương tự đối với một hợp đồng dạng được thanh toán chậm. Sau đây, chúng ta sẽ đưa ra một số công thức giải thích các phân tích trên đây cho trường hợp ma trận xác suất chuyển trạng thái của xích Markov có dạng sau: ⎡1 0 0 0 ⎤ ⎡1 0 0 0 ⎤ ⎢0 1 0 0 ⎥ ⎢0 1 0 0 ⎥ P = ⎢ ⎥ = ⎢ ⎥ (như trong bài toán trên), ⎢p20 p21 p22 p23 ⎥ ⎢0,5 0 0,3 0,2⎥ ⎢ ⎥ ⎢ ⎥ ⎣⎢p30 p31 p32 p33 ⎦⎥ ⎣0, 4 0,3 0,2 0,1⎦ ⎡J Ο ⎤ ⎡1 0⎤ ⎡p20 p21 ⎤ ⎡0 0⎤ ⎡p22 p23 ⎤ = ⎢ ⎥ với J = ⎢ ⎥ , K= ⎢ ⎥ , O = ⎢ ⎥ , M= ⎢ ⎥ , ⎣Κ Μ⎦ ⎣0 1⎦ ⎣p30 p31 ⎦ ⎣0 0⎦ ⎣p32 p33 ⎦ trong đó pij là xác suất chuyển từ trạng thái Si sang trạng thái Sj sau một bước. Không gian trạng thái gồm bốn trạng thái S0, S1, S2 và S3; các trạng thái S0 và S1 là các trạng thái hấp thụ, còn S2 và S3 là các trạng thái truyền ứng. Chúng ta dùng các kí hiệu: ⎡u20 u21 ⎤ U = ⎢ ⎥ , ⎣u30 u31 ⎦ với uik là xác suất hấp thụ vào trạng thái Sk khi trạng thái ban đầu là Si, k = 0, 1, còn i = 2, 3. ⎡v2 ⎤ V = ⎢ ⎥ , ⎣v3 ⎦ với vi là thời gian trung bình cho tới khi rơi vào một trong các trạng thái hấp thụ nếu trạng thái ban đầu là Si, i = 2, 3. ⎡w 22 w 23 ⎤ W = ⎢ ⎥ , ⎣w32 w33 ⎦ với wij là thời gian trung bình xích Markov ở trong trạng thái Sj trước khi nó rơi vào một trong các trạng thái hấp thụ nếu trạng thái ban đầu là Si, i = 2, 3. Lúc đó: Trường Đại học Nông nghiệp Hà Nội – Giáo trình Vận trù học 146
49. −1 −1 ⎡1⎤ −1 ⎡1 0⎤ U = (I - M) , V = (I - M) × ⎢ ⎥ và W = (I - M) × ⎢ ⎥ . ⎣1⎦ ⎣0 1⎦ Chú ý: Việc chứng minh các công thức trên cho trường hợp tổng quát thực ra cũng không quá khó, có thể tìm thấy trong các sách tham khảo về quá trình Markov. Cách ứng dụng phân tích Markov như trong mục này còn có thể được áp dụng rộng rãi trong nhiều lĩnh vực khác như Sinh học, Xã hội học, Lí thuyết nhận dạng và Thiết kế các hệ thống kĩ thuật, trong đó có Kĩ thuật điện. 2.4. Tìm phân phối giới hạn cho một hệ thống kĩ thuật Một hệ thống kĩ thuật có hai chi tiết có thể bị hỏng ở bất kì thời điểm nào. Tại mỗi thời điểm hệ thống có thể rơi vào một trong những trạng thái sau (xem hình V.2): S0: cả 2 chi tiết tốt; S1: chi tiết 1 hỏng, chi tiết 2 bình thường; S2: chi tiết 1 bình thường, chi tiết 2 hỏng; S3: cả hai chi tiết đều hỏng. S0 S1 S2 S 3 Hình V.2. Sơ đồ các trạng thái Nói cách khác, tại mỗi thời điểm t, biến X(t) có thể rơi vào một trong các vị trí/trạng thái S0, S1, S2 và S3. Chú ý rằng lúc này ta có xích Markov (thời gian) liên tục với không gian trạng thái S ={S0, S1, S2, S3}. Sau đây, chúng ta sẽ tìm cách xác định phân phối giới hạn (long run distribution) của {X(t)}t≥0. Đây là một vấn đề khá phức tạp nên chúng ta chỉ có thể trình bày vấn đề một cách vắn tắt. Trước hết ta nhắc lại về phân phối Poát−xông và phân phối mũ. Giả sử dòng tín hiệu đến (hay xảy ra) tuân theo phân phối Poát−xông P (λ) với λ là số tín hiệu đến trung bình trong một khoảng thời gian nhất định (coi là một đơn vị thời gian), λ còn được gọi là cường độ của dòng tín hiệu đến. Lúc đó, trong khoảng thời gian như trên thì λke−λ số tín hiệu xảy ra sẽ nhận giá trị k với xác suất . Ta gọi phần tử xác suất P là xác k!
50. suất xuất hiện (ít nhất) một tín hiệu trong khoảng thời gian Δ t. Thế thì, do tính “đơn nhất” của quá trình Poát−xông, P cũng là xác suất xuất hiện đúng một tín hiệu trong khoảng thời gian Δ t. Theo công thức đã biết thì P = λΔt (chính xác tới vô cùng bé o( Δ t)). Chẳng hạn, nếu λ = 6 tín hiệu/1 phút và Δt = 2 giây, ta sẽ có P = λΔt = 6 × (1/30) = 1/5 = 0,2. Từ đó, ta thấy xác suất để có 1 tín hiệu đến trong khoảng thời gian 2 giây là 0,2. Xét biến ngẫu nhiên T (chẳng hạn thời gian phục vụ một tín hiệu trong một hệ dịch vụ), có phân phối mũ ε(μ) với hàm mật độ là f(τ) = μe−μτ. μ cũng được gọi là cường độ phục vụ hay cường độ của “dòng phục vụ”. Hàm phân phối xác suất của T sẽ là ττ F (τ) = P (T ≤ τ) = ∫∫f(t)dt=μ e−μ−μτt dt =− 1 e . 00 Còn kì vọng toán và độ lệch chuẩn của biến ngẫu nhiên T là +∞ +∞ 1 1 −μt σ = mT = ∫∫tf (t)dt= μ= te dt , T . 00 μ μ Ta nhận thấy ngay rằng: P (0≤≤Δ= T t) F(Δt) - F(0) = 1 − e −μΔt - [1 − e 0] = 1 − e −μΔt = μ Δt (chính xác tới vô cùng bé o( Δ t)). Chú ý: Nếu dòng tín hiệu đến có phân phối Poát−xông P (λ) thì thời gian giữa hai tín hiệu liên tiếp có phân phối mũ ε(λ). Chúng ta quay lại bài toán đang xét. Gọi λ1 số lần chi tiết 1 hỏng và λ2 số lần chi tiết 2 hỏng (tính trung bình) trên 1 đơn vị thời gian. Lúc đó, ta có thể coi dòng tín hiệu chi tiết 1 và 2 hỏng là dòng Poát−xông với các tham số λ1 và λ2. Gọi T1 và T2 là thời gian sửa chữa chi tiết 1 và 2, có phân phối mũ với các kì vọng tsc1 và tsc2 là thời gian sửa chữa (trung bình) chi tiết 1 và chi tiết 2. Vậy T1 và T2 có phân phối mũ ε(μ1) và ε(μ2), với μ1 = 1/tsc1 và μ2 = 1/tsc2. Tại thời điểm t ta có biến ngẫu nhiên X(t) = Xt với phân phối xác suất sau đây: Xt S0 S1 S2 S3 P π0(t) π1(t) π2(t) π3(t) Ta tính π0(t + Δt ) tại thời điểm tiếp theo (t + Δt ) trong hai trường hợp sau đây: − Trường hợp 1: Tại thời điểm t, hệ thống ở trạng thái S0 và tại thời điểm t + Δt , hệ thống vẫn ở trạng thái S0 (không hỏng). Trường Đại học Nông nghiệp Hà Nội – Giáo trình Vận trù học 148
51. − Trường hợp 2: Tại thời điểm t, hệ thống ở trạng thái S1 hoặc S2, còn tại thời điểm t + Δt hệ thống ở trạng thái S0. Do đó, π0(t + Δt) = π0(t) [1 − (λ1 + λ2)Δt] + π1(t) μ1 Δt + π2(t) μ2 Δt (*). Thật vậy, xác suất do trường hợp 2 gây nên là π1(t)μ1 Δt + π2(t)μ2 Δt, với μ1Δt = P(0 ≤≤ΔTt)1 là xác suất sửa chữa xong chi tiết 1 trong khoảng thời gian Δt và μ2Δt= P(0 ≤≤ΔTt)2 là xác suất sửa chữa xong chi tiết 2 trong khoảng thời gian Δt. Trong khi đó, xác suất do trường hợp 1 gây nên là π0(t)[1 − (λ1 + λ2)Δt], với λ1Δt: xác suất hỏng chi tiết 1 trong khoảng Δt, còn λ2Δt: xác suất hỏng chi tiết 2 trong khoảng Δt. Nói cách khác, chúng ta đã thực hiện công thức xác suất đầy đủ π0(t + Δt) = π0(t)p00 + π1(t)p10 + π2(t)p20, trong đó: pi0 là xác suất hệ ở trạng thái Si tại thời điểm t và chuyển sang trạng thái S0 tại thời điểm (t+ Δt). π+Δ−π(t t) (t) Từ (*) ta có: 00= πμ+πμ−πλ−πλ(t) (t) (t) (t) . Δt 11220102 Cho Δt → 0 (vế phải không liên quan với Δt) thì d(t)π 0 =π(t) μ +π (t) μ −π (t) λ −π (t) λ dt 11220102 Khi t rất lớn (hệ thống hoạt động trong một khoảng thời gian đủ dài) thì hệ thống dần ổn định với phân phối giới hạn có thể tìm được, tức là: [π0(t), π1(t), π2(t), π3(t)] → [π0, π1, π2, π3]. Vậy ta có: d(t)π0 π1μ1 + π1μ2 − π0λ1 − π0λ2 = 0 (vì = 0 khi t đủ lớn). dt Một cách tương tự, ta đi đến hệ phương trình: ⎧dπ 0 =μπ+μπ−λ+λπ()0 = ⎪ dt 11 2 2 1 2 0 ⎪ dπ ⎪ 1 =λπ +μ π −()0 λ +μ π = ⎪ dt 10 23 2 1 1 ⎨ dπ ⎪ 2 =λ π +μπ −()0 λ +μ π = ⎪ dt 20 13 1 2 2 ⎪dπ ⎪ 3 =λ π +λπ −()0 μ +μ π = ⎩ dt 21 12 1 2 3 Một cách tổng quát, phân phối giới hạn được tìm từ hệ phương trình: −πjjjqq =∑ π iij ij≠ π=q0 hay ∑ iij ,∀j ∈S ( ) và ∑πi =1, trong đó − qii là cường độ chuyển từ trạng thái iS∈ iS∈
52. i sang các trạng thái khác (không kể i), còn qij là cường độ chuyển từ trạng thái i sang trạng thái j, được định nghĩa như sau: qii = − limΔt→0 (P[X(t+ Δ≠ t) i/ X(t) = i]/ Δ t), qij = limΔt→0 (P[X(t+ Δ= t) j/ X(t) = i]/ Δ t) , Lúc đó, Q = [qij] được gọi là ma trận cường độ. Từ điều kiện ( ) ta thấy, để tìm T T phân phối giới hạn cần phải giải hệ [π0 π1 π2 π3]Q = 0 hay Q [π0 π1 π2 π3] = 0. Ví dụ 1: Cho λ1 = 1, λ2 = 2, μ1 = 2, μ2 = 3. Từ sơ đồ cường độ chuyển trạng thái cho trên hình V.3, có thể tìm được ma trận cường độ Q, với QT có dạng sau: ⎡⎤−3230 ⎢⎥1403− QT = ⎢⎥ ⎢⎥2042− ⎢⎥ ⎣⎦0215− S0 μ1 μ2 λ λ1 2 S1 S2 μ μ 2 1 S3 λ2 λ1 Hình V.3. Sơ đồ cường độ chuyển trạng thái Giải thích: q00 = − 3 do cường độ chuyển từ trạng thái S0 sang các trạng thái khác là λ1+ λ2 = 3, còn q10 = 2 là cường độ chuyển từ trạng thái S1 vào trạng thái S0. T T Giải hệ [π0 π1 π2 π3]Q = 0 hay Q [π0 π1 π2 π3] = 0 (với điều kiện bổ trợ π0 + π1 + π2 + π3 = 1) có kết quả: π=0 6/15 = 0,4; π1 ==3/15 0,2; π2 ==4/15 0,27; π=3 2/15 = 0,13. Cần chú ý rằng, hệ [π0 π1 π2 π3] Q = 0 theo một nghĩa nhất định là tương tự với hệ Π ×(I - P) = 0, như đã trình bày trong các mục 1.2 và 2.1. Giả sử lợi nhuận/1 đơn vị thời gian hệ thống mang lại trong các trường hợp có thể xảy ra như sau: nếu hệ thống trong trạng thái S0 thì lợi nhuận là 8 USD, tại S1 là 3 USD, tại S2 là 5 USD, tại S3 là 0 USD. Vậy lợi nhuận trung bình/1 đơn vị thời gian là 8 × 0,4 + 3 × 0,2 + 5 × 0,27 = 5,15 (USD). Qua ví dụ ta thấy π0, π1, π2, π3 được xác định căn cứ vào các giá trị đã biết λ1, λ2, μ1, μ2. Trường Đại học Nông nghiệp Hà Nội – Giáo trình Vận trù học 150
53. λ1, λ2: số lần chi tiết hỏng (tuỳ thuộc hệ thống cụ thể), μ1, μ2: các tham số sửa chữa cần đưa vào. Lợi nhuận cuối cùng của hệ thống phụ thuộc vào λ1, λ2, μ1, μ2 và được xác định bằng cách giải bài toán tối ưu sau: Lợi nhuận L = c0π0 + c1π1 + c2π2 → Max (c0, c1, c2: lợi nhuận từng trạng thái) với các ràng buộc: ⎧μπ11 +μπ 2 2 −()0 λ 1 +λ 2 π 0 = ⎪ ⎪λπ10 +μπ 23 −()0 λ 2 +μ 1 π 1 = ⎪ ⎪λπ +μπ −()0 λ +μ π = ⎪ 20 13 1 2 2 ⎨ ⎪λπ+λπ21 12 −μ+μ()0 1 2 π 3 = ⎪ ⎪π0123 +π +π +π =1 ⎪ ⎪ ⎩ππππ≥0123,,, 0;, μμ≥ 12 0 Lưu ý: Bài toán trên đây có 6 biến (λ1, λ2 đã biết). Ta phải tìm được μ1, μ2 từ bài toán để có phương hướng xây dựng hệ thống với lợi nhuận lớn nhất. 2.5. Một ứng dụng của quá trình sinh - tử cho hệ thống hàng chờ Quá trình sinh−tử được ứng dụng khá rộng rãi trong Lí thuyết độ tin cậy, là một môn học của ngành Điện/Điện tử và một số ngành khoa học kĩ thuật khác cũng như trong Quản trị kinh doanh và Vận trù học. Quá trình sinh − tử là trường hợp riêng của xích Markov thuần nhất thời gian liên tục, với không gian trạng thái S không quá đếm được S = {S0, S1, S2, , Sn, } và ma trận cường độ Q = [qij] có tính chất qij = 0 với ⏐i − j⏐≥ 2. Điều này có nghĩa là việc chuyển trạng thái trong quá trình sinh−tử chỉ có thể tới “1 đơn vị lên hoặc xuống” (xem hình V.3). λ0 λ1 λn-1 λn S0 S1 S2 Sn-1 Sn Sn+1 μ1 μ2 μn μn+1 Hình V.3. Sơ đồ chuyển trạng thái trong quá trình sinh− tử Từ trạng thái Sn tại thời điểm t hệ X(t) chỉ có thể chuyển tới một trong các trạng thái Sn+1, Sn hoặc Sn−1. Vì vậy chúng ta có các cường độ chuyển: μ0 = λ−1 = 0, q00 = − λ0, qn, n+1 = λn, qn, n−1 = μn và qn, n = − (λn +μn) ∀n.
54. Trong trường hợp λn, μn > 0, ∀n > 0, theo định lí đã được chứng minh, phân phối giới hạn có thể tìm được bằng cách giải hệ: [π0 π1 π2 π3 ]Q = 0, với ma trận cường độ Q đã biết. Ma trận chuyển vị của Q có dạng: q ⎡q00 q10 qn0 n1,0+ ⎤ ⎢ q q q q ⎥ ⎢ 01 11 n1 n1,1+ ⎥ QT = ⎢ ⎥ ⎢ ⎥ ⎢q0n q1n qnn qn1,n+ ⎥ ⎢ ⎥ ⎣ ⎦ T T Ta có [π0 π1 π2 π3 ]Q = 0 ⇔ Q [π0 π1 π2 π3 ] = 0 ⇔ ⎡⎤⎡⎤⎡⎤q00 q 10 q 20 π 0 0 ⎢⎥⎢⎥⎢⎥q q q π 0 ⎢⎥⎢⎥⎢⎥01 11 21×= 1 . ⎢⎥⎢⎥⎢⎥q02 q 12 q 22 π 2 0 ⎢⎥⎢⎥⎢⎥ ⎣⎦⎣⎦⎣⎦ hay: q00π0 + q10π1 + q20π2 + = 0, q01π0 + q11π1 + q21π2 + = 0, q02π0 + q12π1 + q22π2 + = 0, Do tính chất đặc biệt, như đã phân tích ở trên, của ma trận cường độ Q của quá trình sinh−tử, hệ trên được viết một cách tường minh hơn như sau: −λ0π0 + μ1π1 + = 0, λ0π0 − (λ1 + μ1)π1 + μ2π2 + = 0, λ1π1 − (λ2 + μ2)π2 + μ3π3 + = 0, Từ đây đễ dàng tìm được πn+1 = (λn/μn+1)πn, ∀n = 1, 2, 3, để đi tới công thức tính πi,∀i. π1 = (λ0/μ1)π0, π2 = (λ1/μ2)π1 = (λ1λ0/μ2μ1)π0, π3 = (λ2/μ3)π2 = (λ2λ1/μ3μ2)π1 = (λ2λ1λ0/μ3μ2μ1)π0, πn+1 = (λn/μn+1)πn = = (λnλn−1 λ0/μn+1μn μ1)π0, Trường Đại học Nông nghiệp Hà Nội – Giáo trình Vận trù học 152