محاكاة مونت كارلو

تجيب محاكاة مونت كارلو (Monte Carlo simulation) عن سؤال يخص نظاما غير مؤكد بأن تسحب أرقاما عشوائية وتقرأ كل رقم منها على أنه حدث: طلب يوم أو زمن توريد أو علاج أو زمن خدمة. ونحل في هذا القسم أمثلة التمرين 5 الخمسة كلها خطوة بخطوة: التحكم في المخزون وعيادة طبيب الأسنان والتحويل العكسي للتوزيع الأسي وطابور بأزمنة أسية ومنشأة تصنيع.

الأهداف

• تحويل جدول التكرار إلى فترات أرقام عشوائية (random number intervals) من خانتين.
• حساب متوسط الطلب وزمن التوريد وانحرافهما المعياري ثم كمية الطلب الاقتصادية (EOQ) ومستوى إعادة الطلب (ROL).
• بناء جدول محاكاة المخزون وتحديد موعد الطلبية.
• محاكاة طابور بوصول مجدول وإيجاد متوسط زمن الانتظار.
• توليد متغيرات عشوائية أسية (exponential random variates) باستخدام التحويل العكسي (inverse transform).
• استخدام المولد المتطابق المختلط (mixed congruential generator) لمحاكاة منشأة تصنيع.

1. المثال 1: التحكم في المخزون

لنفترض أننا نعرف الطلب (demand) في آخر 100 يوم وزمن التوريد (lead time) الذي وصلت فيه الكمية المطلوبة. ونحتاج إلى معرفة الطلب في الأيام القادمة وزمن التوريد المحتمل وموعد الطلبية.

الطلب في آخر 100 يوم:

زمن التوريد:

قد يكون لزمن التوريد توزيع تكراري قياسي بدلا من التكرارات المعطاة لإيجاد الأرقام العشوائية المقابلة لزمن توريد معين.

ولنفترض أن التحكم في المخزون له هذه التكاليف:

عند حدوث عجز (shortage) يجب إحضار الأصناف من السوق المحلية بسعر أعلى.

الخطوة 1: المتوسط والانحراف المعياري

نحصل على مستوى إعادة الطلب والكميات الأخرى بحسابات إحصائية على البيانات المتاحة. ومجموع قيم الطلب المئة هو $\sum D_i = 4523$ ومجموع مربعاتها هو $\sum D_i^2 = 214077$.

متوسط الطلب وانحرافه المعياري (standard deviation) هما:

$$\begin{aligned} \mu_D &= \dfrac{\sum D_i}{n} = \dfrac{4523}{100} \\ &= 45.23 \text{ unit/day} \end{aligned}$$ $$\begin{aligned} \sigma_D &= \sqrt{\dfrac{\sum D_i^2 - n (\mu_D)^2}{n - 1}} \\ &= \sqrt{\dfrac{214077 - 100 (45.23)^2}{99}} \\ &= 9.8 \end{aligned}$$

ومتوسط زمن التوريد وانحرافه المعياري هما:

$$\begin{aligned} \mu_L &= \dfrac{(28 \times 3) + (44 \times 4) + (28 \times 5)}{100} \\ &= \dfrac{400}{100} \\ &= 4 \text{ days} \end{aligned}$$ $$\begin{aligned} & 28 (3)^2 + 44 (4)^2 + 28 (5)^2 - 100 (4)^2 \\ &= 252 + 704 + 700 - 1600 \\ &= 56 \end{aligned}$$ $$\begin{aligned} \sigma_L &= \sqrt{\dfrac{56}{99}} \\ &= 0.75 \end{aligned}$$

ملاحظة على ورقة التمرين. تكتب الورقة متوسط زمن التوريد هكذا: 400/4 = 4. والقسمة على 100: أي 400 / 100 = 4. والناتج 4 أيام لا يتغير.

الخطوة 2: الأرقام العشوائية للطلب

نوجد المدرج التكراري (histogram) للطلب ونحسب التكرار (frequency) والتكرار التراكمي (cumulative frequency) ثم نخصص لكل فترة أرقاما عشوائية من خانتين (من 00 إلى 99) بحسب تكرارها. فالفترة التي تكرارها f تأخذ f من الأرقام وتنتهي عند تكرارها التراكمي ناقص 1. وعندما نتوقع الطلب نأخذ القيمة المتوسطة للطلب في كل فترة.

ملاحظة على ورقة التمرين. تقول الورقة إن الطلب يمتد من 20 إلى 70. والبيانات تمتد من 22 إلى 72. والجدول صحيح: ففتراته الإحدى عشرة من 20 - 24 إلى 70 - 74 تضم القيم المئة كلها. وإعادة عد القيم تعطي التكرارات نفسها.

الخطوة 3: الأرقام العشوائية لزمن التوريد

ينطبق الحساب نفسه على زمن التوريد:

الخطوة 4: كمية الطلب الاقتصادية

$$\begin{aligned} EOQ &= \sqrt{\dfrac{2 \mu_D C_0}{C_i}} \\ &= \sqrt{\dfrac{2 \times 45.23 \times 100}{0.1}} \\ &\approx 300 \text{ units} \end{aligned}$$

كلما هبط المخزون تحت حد معين هو مستوى إعادة الطلب (reorder level) نطلب هذه الكمية نفسها وهي 300 وحدة لأن هذا نظام كمية ثابتة (fixed quantity system).

الخطوة 5: مستوى إعادة الطلب

$$\begin{aligned} ROL &= \mu_D \mu_L \\ &\quad + K \sqrt{\mu_L \sigma_D^2 + \mu_D^2 \sigma_L^2} \end{aligned}$$

K ثابت يحدد احتمال حدوث عجز في المخزون. فعند K = 1.65 يكون احتمال العجز 5% أي أننا متأكدون بنسبة 95% من عدم حدوث عجز. وعند K = 2.33 لا يتجاوز احتمال العجز 1%.

عند K = 1.65:

$$\sqrt{4 (9.8)^2 + (45.23)^2 (0.75)^2} = 39.18$$ $$\begin{aligned} ROL &= (45.23 \times 4) + 1.65 \times 39.18 \\ &= 180.92 + 1.65 \times 39.18 \\ &\approx 245.6 \end{aligned}$$

ملاحظة على ورقة التمرين. تكتب الورقة ROL = 245. والقيمة الدقيقة نحو 245.6 فالرقم 245 مقرب إلى الأسفل. ونبقي على 245 لأن جدول المحاكاة يستخدمها. ومع 246 تصدر الطلبيات في الأيام نفسها لأن أقل كمية لا تصدر عندها طلبية هي 248.

الخطوة 6: بدء المحاكاة

نبدأ المحاكاة باستخدام EOQ = 300 مع ROL = 245 وحدة. ونفترض أن المخزون الابتدائي 300 وحدة مساء الجمعة فتبدأ المحاكاة من يوم الاثنين.

يمكن توليد الأرقام العشوائية من 00 إلى 99 في Excel أو باستخدام مولدات أرقام عشوائية أخرى.

ملاحظة على ورقة التمرين. تطلب الورقة أن نضع 100 داخل الدالة rand(). لكن RAND() في Excel لا تأخذ وسيطا وتعطي رقما أقل من 1. فاضرب في 100 بدلا من ذلك: INT(RAND()*100) تعطي رقما من 00 إلى 99.

تأخذ الورقة أرقامها العشوائية من جدول مطبوع لأرقام من خانتين:

• للطلب من يسار الصف الأول: 81, 78, 89, 59, 40, 33, 60, 67, 19, 21.
• لزمن التوريد من يمين الصف الأخير: 65, 67, 22. ولا نسحب منها إلا عند إصدار طلبية.

ملاحظة على ورقة التمرين. يذكر النص أرقام الطلب هكذا: 81, 78, 59 فيسقط منها رقم. أما الصف الأول وجدول المحاكاة ففيهما 81, 78, 89, 59. ويوم الأربعاء يستخدم 89 فيكون طلبه 57.

الخطوة 7: جدول المحاكاة

يغطي الجدول أسبوعي عمل (10 أيام عمل) والمخزون الابتدائي 300 وحدة وكمية الطلب الاقتصادية EOQ تساوي 300 وحدة ومستوى إعادة الطلب ROL يساوي 245 وحدة:

القواعد التي يبنى عليها كل صف:

1. المخزون الابتدائي هو المخزون النهائي السابق.
2. نحصل على الطلب من الرقم العشوائي عبر جدول الخطوة 2: يقع 81 في 69 - 81 فيكون الطلب 52.
3. المخزون النهائي = المخزون الابتدائي - الطلب اليومي + الكمية المستلمة.
4. الكمية المتاحة + الطلبية = المخزون النهائي + الطلبية التي في الطريق.
5. إذا هبطت هذه الكمية تحت ROL فإننا نطلب EOQ ونسحب زمن توريد من جدول الخطوة 3.

والآن نتتبع الأيام الأولى:

• الاثنين الأول: 300 - 52 = 248 وحدة. وهذا فوق 245 فلا تصدر طلبية.
• الثلاثاء: المخزون الابتدائي 248 والمخزون النهائي 248 - 52 = 196. وهذا تحت 245 فيجب إصدار طلبية. ويقع الرقم العشوائي 65 في 28 - 71 فيكون زمن التوريد 4 أيام.
• من الأربعاء إلى الجمعة: المخزون النهائي تحت 245 لكن الوحدات الـ300 التي في الطريق تحتسب: 139 + 300 = 439 ثم 92 + 300 = 392 ثم 50 + 300 = 350. فلا تصدر طلبية جديدة.
• الأربعاء الثاني: هبط المخزون إلى 214 وليس في الطريق شيء فنصدر طلبية أخرى. ويعطي الرقم العشوائي 67 زمن توريد قدره 4 أيام.

ملاحظة على ورقة التمرين. تعرف الورقة الكمية المتاحة + الطلبية بأنها المخزون النهائي + الكمية المطلوبة. وهذا يعطي في يوم الطلبية 196 + 300 = 496 ثم 214 + 300 = 514 لكن الجدول يعرض 196 ثم 214. ففي الجدول تدخل الطلبية هذا العمود من اليوم التالي. ونبقي على الجدول وأيام الطلبيات واحدة في الحالتين.

الخطوة 8: وصول الطلبيات والتكاليف

عندما يكون زمن التوريد 4 أيام عمل من الثلاثاء تصل الكمية المطلوبة يوم الاثنين التالي ويكون المخزون النهائي:

$$(50 - 42) + 300 = 308 \text{ units}$$

وتصل الكمية الثانية يوم الثلاثاء التالي بعد الأيام العشرة في الجدول.

ويمكننا أيضا حساب مؤشرات مثل التكلفة في مدة معينة:

• تكلفة الاحتفاظ بالمخزون = تكلفة الاحتفاظ بالوحدة × المخزون النهائي.
• تكلفة إصدار الطلبيات = عدد الطلبيات × تكلفة الطلبية.

2. المثال 2: عيادة طبيب الأسنان

في عيادة طبيب أسنان طبق محاكاة مونت كارلو لإيجاد متوسط زمن انتظار العملاء. تبدأ المحاكاة عند الساعة 8 صباحا والوصول مجدول (وقد يكون في حالات أخرى عشوائيا بحسب توزيع احتمالي). يصل العميل الأول عند الساعة 8 صباحا ثم يأتي عميل كل 30 دقيقة من أجل واحد من خمسة علاجات مختلفة. والخدمات العلاجية الخمس هي الحشو والتاج والتنظيف والخلع والفحص ولكل منها احتمال وزمن خدمة.

الخطوة 1: فترات الأرقام العشوائية

لا نعرف الخدمة العلاجية التي يأتي من أجلها العميل. ونخمنها بأرقام عشوائية من أي مولد عبر هذا الجدول:

الخطوة 2: أي خدمة؟

تحدد الأرقام العشوائية 40, 82, 11, 34, 25, 66, 17, 79 الخدمة لـ8 مرضى:

لاحظ أن 40 تاج وليس حشوا: ففترة الحشو تنتهي عند 39.

الخطوة 3: جدول المحاكاة

تبدأ الخدمة عند الأحدث بين الوصول ونهاية الخدمة السابقة. وزمن الانتظار هو بداية الخدمة مطروحا منها الوصول.

فمثلا يصل المريض 2 عند 8:30 والطبيب مشغول حتى 9:00 فيكون الانتظار 30 دقيقة. وتنتهي هذه الخدمة عند 9:15 بعد 15 دقيقة فينتظر المريض 3 الذي يصل عند 9:00 مدة 15 دقيقة.

متوسط زمن الانتظار هو:

$$\begin{aligned} & 0 + 30 + 15 + 25 + 35 \\ & \quad + 45 + 30 + 40 = 220 \end{aligned}$$ $$\dfrac{220}{8} = 27.5 \text{ minutes}$$

ملاحظة على ورقة التمرين. تطلب الورقة متوسط زمن الانتظار لكنها لا تذكره. ويمتد جدولها أيضا خارج الصفحة: فالمطبوع منه 6 من أزمنة الانتظار الثمانية ومعها 7 من نهايات الخدمة الثماني. أما القيم الأخيرة (النهاية 12:55 وزمنا الانتظار 30 ثم 40) والمتوسط 27.5 دقيقة فناتجة من القاعدة نفسها.

3. المثال 3: التحويل العكسي للتوزيع الأسي

ليكن للمتغير $X$ التوزيع الأسي (exponential distribution) بمتوسط $\beta$. دالة التوزيع هي:

$$F(x) = \begin{cases} 1 - e^{-x/\beta} & \text{if } x \ge 0 \\ 0 & \text{otherwise} \end{cases}$$

لإيجاد $F^{-1}$ نضع $u = F(x)$ ونحل بالنسبة إلى $x$:

$$\begin{aligned} u &= 1 - e^{-x/\beta} \\ e^{-x/\beta} &= 1 - u \\ F^{-1}(u) = x &= -\beta \ln (1 - u) \end{aligned}$$

ولتوليد المتغير العشوائي (random variate) المطلوب نولد أولا $U \sim U(0, 1)$ ثم نجعل:

$$X = -\beta \ln U$$

ويمكن في هذه الحالة استخدام $U$ بدلا من $1 - U$ لأن $1 - U$ لها التوزيع $U(0, 1)$ نفسه الذي للمتغير $U$. وهذا يوفر عملية طرح.

4. المثال 4: طابور بأزمنة أسية

لنفترض أن زمن ما بين الوصول (inter-arrival time) في نظام طابور موزع توزيعا أسيا بمتوسط $\beta_1 = 2$ وأن زمن الخدمة (service time) موزع توزيعا أسيا أيضا بمتوسط $\beta_2 = 1.33$. من جدول المحاكاة لأزمنة الوصول والخدمة والمغادرة لـ8 عملاء احسب متوسط زمن انتظار العميل.

الخطوة 1: المتغيرات العشوائية

نحصل على زمن ما بين الوصول $X_i$ وزمن الخدمة $Y_i$ من معادلتي التحويل العكسي:

$$\begin{aligned} X_i &= -\beta_1 \ln U_1 \\ Y_i &= -\beta_2 \ln U_2 \end{aligned}$$

وتقع $U_1$ مع $U_2$ في الفترة من 0 إلى 1. وإذا أعطينا متتاليتين من الأرقام العشوائية ذات الخانتين $u_1, u_2$ من 00 إلى 99 نحصل على $U_1, U_2$ بقسمة كل رقم على 100. والمتتاليتان هما:

• u1 = 20, 23, 86, 09, 92, 35, 38, 01
• u2 = 21, 44, 27, 70, 73, 36, 59, 42

للعميل الأول:

$$\begin{aligned} U_1 &= 20 / 100 = 0.20 \\ X_1 &= -2 \ln (0.20) = 3.22 \\ U_2 &= 21 / 100 = 0.21 \\ Y_1 &= -1.33 \ln (0.21) = 2.08 \end{aligned}$$

الخطوة 2: جدول المحاكاة

يقرب كل متغير إلى خانتين عشريتين ثم يجمع كما تفعل الورقة. وزمن الوصول هو الوصول السابق زائد $X_i$. وتبدأ الخدمة عند الأحدث بين الوصول والمغادرة السابقة ويغادر العميل بعد $Y_i$. وزمن الانتظار هو بداية الخدمة مطروحا منها الوصول.

لا ينتظر إلا عميلان:

• العميل 3 يصل عند 6.46 ويغادر العميل 2 عند 7.25: فالانتظار 7.25 - 6.46 = 0.79 والمغادرة 7.25 + 1.74 = 8.99.
• العميل 5 يصل عند 11.45 ويغادر العميل 4 عند 11.75: فالانتظار 11.75 - 11.45 = 0.30 والمغادرة 11.75 + 0.42 = 12.17.

متوسط زمن الانتظار هو:

$$\dfrac{0.79 + 0.30}{8} = \dfrac{1.09}{8} \approx 0.14$$

ويمكننا حساب مقاييس أداء أخرى من الجدول نفسه.

ملاحظة على ورقة التمرين. عند u2 = 27, 70, 42 تعطي الورقة أزمنة الخدمة 1.75, 0.48, 1.16. ومع القيمة المعطاة $\beta_2 = 1.33$ تكون 1.74, 0.47, 1.15 وهي المستخدمة في الجدول أعلاه. أما عمود الورقة كله فيطابق $\beta_2 = 4/3$. ولذلك تصبح أزمنة المغادرة في الورقة 9.00, 11.76, 12.18, 25.86 هي 8.99, 11.75, 12.17, 25.85 ويصبح الانتظار الثاني 0.31 هو 0.30. أما متوسط الورقة (0.79 + 0.31) / 8 = 0.14 فيقرب إلى القيمة نفسها.

5. المثال 5: التصنيع

لدينا الوحدات المنتجة في منشأة تصنيع وعدد المركبات المتاحة للنقل في محطة معينة بهذه الاحتمالات:

الخطوة 1: المولد المتطابق المختلط

يولد المولد المتطابق المختلط (mixed congruential generator) بمعاملاته الثلاثة a ثم c ثم m أرقاما عشوائية من خانتين من 00 إلى 99:

$$X_i = (a X_{i-1} + c) \bmod m$$

• للمتتالية $u_1$: a = 21 مع c = 3 مع m = 100 والبذرة (seed) $X_0 = 20$.
• للمتتالية $u_2$: a = 41 مع c = 7 مع m = 100 والبذرة $X_0 = 1$.

ملاحظة على ورقة التمرين. تكتب الورقة $X_i = a X_{i-1} + c \bmod (m)$ دون القوسين الخارجيين. وخطواتها تحسب الصيغة ذات القوسين مثل 423 mod 100 = 23. وكل قيم المتتاليتين تتبع هذه الصيغة.

للمتتالية $u_1$:

$$\begin{aligned} X_1 &= (21 \times 20 + 3) \bmod 100 \\ &= 423 \bmod 100 = 23 \\ X_2 &= (21 \times 23 + 3) \bmod 100 \\ &= 486 \bmod 100 = 86 \end{aligned}$$

للمتتالية $u_2$:

$$\begin{aligned} X_1 &= (41 \times 1 + 7) \bmod 100 \\ &= 48 \bmod 100 = 48 \\ X_2 &= (41 \times 48 + 7) \bmod 100 \\ &= 1975 \bmod 100 = 75 \end{aligned}$$

وهكذا حتى $X_9$:

الخطوة 2: جدول المحاكاة لعشرة أيام

تستخدم الأيام العشرة القيم من $X_0$ إلى $X_9$ فيستخدم اليوم 1 البذرة نفسها. ونقرأ كل رقم عبر جدول الفترات أعلاه.

متوسط عدد الوحدات المنتجة في اليوم:

$$\dfrac{6000}{10} = 600 \text{ units/day}$$

متوسط عدد المركبات المتاحة في اليوم:

$$\dfrac{64}{10} = 6.4 \text{ vehicle/day}$$

6. الأخطاء الشائعة

1. استخدام حدي الفترة طلبا. الطلب الناتج من فترة هو متوسطها: يقع 81 في 69 - 81 وهذا الصف يعطي 52.
2. إغفال الحدود. 40 تاج لأن فترة الحشو تنتهي عند 39. أما 79 فهو خلع لأن فترته 70 - 79.
3. مقارنة المخزون النهائي وحده بمستوى إعادة الطلب. في الأربعاء الأول يكون المخزون النهائي 139 لكن الوحدات الـ300 التي في الطريق تعطي 439 فلا تصدر طلبية جديدة.
4. حساب زمن التوريد بالأيام التقويمية. الحساب بأيام العمل: 4 أيام عمل من الثلاثاء هي الاثنين التالي.
5. بدء كل خدمة عند الوصول. تبدأ الخدمة عند الأحدث بين الوصول ونهاية الخدمة السابقة. ويقاس زمن الانتظار من الوصول.
6. نسيان القسمة على 100 أو إسقاط إشارة السالب. $U = u / 100$ ثم $X = -\beta \ln U$. ولأن $\ln U$ سالب فإن إشارة السالب تجعل $X$ موجبا.

7. تدريب

حل كل مسألة ثم راجع الإجابات في نهاية هذه الصفحة.

1. باستخدام جدولي الطلب وزمن التوريد في المثال 1 أوجد الطلب عند الأرقام العشوائية 94, 51, 37 وزمن التوريد عند 20, 71, 72.
2. استخدم صيغتي التكلفة في المثال 1 لإيجاد تكلفة الاحتفاظ بالمخزون وتكلفة إصدار الطلبيات للأيام العشرة في جدول المحاكاة.
3. في عيادة طبيب الأسنان في المثال 2 يصل 4 مرضى كل 30 دقيقة من 8:00 بالأرقام العشوائية 39, 55, 80, 70. أوجد أزمنة الانتظار ومتوسطها.
4. عند $\beta_1 = 2$ مع $\beta_2 = 1.33$ كما في المثال 4 أجر محاكاة لـ3 عملاء عند u1 = 50, 90, 40 مع u2 = 30, 60, 80. أوجد متوسط زمن الانتظار.
5. استخدم مولد $u_1$ في المثال 5 مع البذرة $X_0 = 5$ للحصول على 4 أرقام عشوائية (يستخدم اليوم 1 البذرة) ثم أوجد الوحدات المنتجة في كل يوم ومتوسطها.

أهم النقاط

• تأخذ الفترة من الأرقام ذات الخانتين بقدر تكرارها وتنتهي عند تكرارها التراكمي ناقص 1.
• تحدد EOQ مقدار الطلبية ويحدد ROL موعدها.
• تحتسب الطلبية من اليوم التالي لإصدارها وتدخل المخزون يوم وصولها.
• في الطابور يكون زمن الانتظار هو بداية الخدمة مطروحا منها الوصول.
• يعطي التحويل العكسي المتغير الأسي $X = -\beta \ln U$.
• يحسب المولد المتطابق المختلط $X_i = (a X_{i-1} + c) \bmod m$.

الإجابات

التدريب 1

يقع 94 في 94 - 95 فيكون الطلب 62. ويقع 51 في 26 - 51 فيكون الطلب 42. ويقع 37 أيضا في 26 - 51 فيكون الطلب 42. وفي زمن التوريد يقع 20 في 00 - 27 فيعطي 3 أيام ويقع 71 في 28 - 71 فيعطي 4 أيام ويقع 72 في 72 - 99 فيعطي 5 أيام.

التدريب 2

مجموع المخزون النهائي هو 248 + 196 + 139 + 92 + 50 + 308 + 261 + 214 + 177 + 140 = 1825 وحدة يوم. فتكلفة الاحتفاظ بالمخزون هي 0.1 × 1825 = 182.5 جنيه. وفي الجدول طلبيتان فتكلفة إصدار الطلبيات هي 2 × 100 = 200 جنيه. ولا يوجد عجز فلا توجد تكلفة عجز.

التدريب 3

39 حشو (40 دقيقة) ثم 55 تنظيف (15 دقيقة) ثم 80 فحص (15 دقيقة) ثم 70 خلع (45 دقيقة).

متوسط زمن الانتظار هو 10 / 4 = 2.5 دقيقة.

التدريب 4

متوسط زمن الانتظار هو (1.39 + 0.24) / 3 = 0.54.

التدريب 5

$$\begin{aligned} X_1 &= (21 \times 5 + 3) \bmod 100 \\ &= 108 \bmod 100 = 8 \\ X_2 &= (21 \times 8 + 3) \bmod 100 \\ &= 171 \bmod 100 = 71 \\ X_3 &= (21 \times 71 + 3) \bmod 100 \\ &= 1494 \bmod 100 = 94 \end{aligned}$$

وتعطي الأرقام 05, 08, 71, 94 الوحدات 500 ثم 550 ثم 650 ثم 700. والمجموع 2400 فيكون المتوسط 2400 / 4 = 600 وحدة في اليوم.