حل یک مدل ریاضی جدید برای زمانبندی در شبکه های توزیع با بهینه سازی ذرات انبوه چند هدفه

1- مقدمه

1-1- بیان مسأله

مدیریت زنجیره تأمین از فلسفه­هایی است که در چند دهه اخیر به علت افزایش روزافزون رقابت‌پذیری و تلاش سازمان­ها برای بقا و با تکیه بر پیشرفت­های حاصل در تکنولوژی اطلاعات و نزدیک شدن ارتباطات مورد توجه سازمان­ها قرار گرفته است. رمز بقای سازمان‌ها در ارضای نیازهای مشتریان است و نگرش مدیریت زنجیره تامین[1] (SCM) نه فقط موجودیت نهایی در ارتباط با مشتری که محصول نهایی را تحویل می­دهد، بلکه یک توالی از تأمین‌کنندگان را در این راستا مؤثر می­داند و ”یکپارچه­سازی[2]“ سازمان­های درگیر و ”هماهنگ‌سازی[3] “ جریان­های مواد، اطلاعات و مالی را برای بهبود رقابت پذیری زنجیره تأمین بررسی می‌کند،. در بحث مدیریت زنجیره تأمین در سال‌های اخیر، توزیع کالا به مشتریان از طریق شبکه­های توزیع مورد توجه بسیاری از محققان بوده است. در هنگام تکمیل سفارش‌های کالاها دو فاکتور اساسی، هزینه و تحویل هستند.  

یکی از بحث­های مطرح در زمینه توزیع کارآمد کالاها، تحویل بموقع در شبکه­های توزیع است. از یک طرف، از آنجا که هزینه­های لجستیکی بخش عمده­ای از هزینه­های یک زنجیره تأمین را تشکیل می­دهند، هماهنگ سازی برنامه­ریزی توزیع و حمل و نقل و کنترل موجودی­ها به­عنوان فرآیندهای اصلی لجستیکی در این راستا می­تواند صرفه­جویی قابل توجهی را در کل هزینه­های زنجیره تأمین موجب شود؛ و از طرف دیگر، نقش بسزایی در تعیین سطح خدمت مشتری خواهد داشت.

این مسأله، تصمیمات توزیع را در یک زنجیره تأمین سه سطحی در افق برنامه­ریزی [1,T]  بررسی می کند. زنجیره تأمینی با چند عمده فروش و چند خرده فروش خواهیم داشت که چندین قلم کالای مورد نیاز مشتریان را تأمین می‌کنند. کالاها از چند تأمین کننده خارجی تأمین می‌شود و توسط دو کانال حمل مستقیم (ارسال مستقیم از تأمین­کنندگان به خرده­فروشان) و غیر مستقیم (ارسال از تأمین­کنندگان به عمده­فروشان و سپس به خرده­فروشان) توزیع می‌شود. در ضمن، هر تأمین­کننده می­تواند بیش از یک نوع قلم را تأمین کند. تقاضای کالاها در محل خرده فروشی­ها قطعی و پویا است. مدت زمان تحویل کالاها توسط تأمین کنندگان به خرده ­فروشان یا انبارهای عمده­فروشان پس از دریافت سفارش از آنها به صورت قطعی، پویا و مضرب صحیحی از طول یک دوره است. خرده فروشان می­توانند هر کالا را از چندین انبار تهیه کنند یا به چندین تأمین ­کننده سفارش بدهند. انبارهای عمده­ فروشی نیز می­توانند هر کالا را به چندین تأمین­کننده سفارش بدهند. به همین، ترتیب عناصر هر رده می­توانند به چندین عنصر رده پایین­تر خود سرویس بدهند. بعلاوه، خروج کالاها از انبارها براساس قاعده اولین ورود-اولین خروج[4]انجام می‌شود.این مقاله در قالب مقاله علمی- پژوهشی ارائه شده است.

1-2- اهداف تحقیق

اهداف این تحقیق را می­توان به شرح ذیل خلاصه نمود:

الف) طراحی یک مدل برنامه­ریزی ریاضی چند هدفه برای مسأله حداقل سازی هزینه­ها و زمانبندی آنها؛ به گونه‌ای که اهداف مطلوب حاصل گردند. بنابراین هدف از طرح این مسأله این است که چه مقدار از هر یک از کالاها باید در هر دوره از هر مبدأ به هر مقصد حمل شود تا با توجه به محدودیت­های مسأله توابع هدف کمینه شوند.

ب) تایید روش طراحی فوق از طریق حل مسأله با استفاده از روش ریاضی و طراحی روند مناسب برای به کارگیری الگوریتم بهینه­سازی ذرات انبوه چندهدفه برای کاربرد در مسایل دنیای واقعی.

1-3- مرور ادبیات موضوع

فلسفه مدیریت زنجیره تأمین از اوایل دهه 1980 مورد توجه قرار گرفت و مقالات و تحقیقات موجود به سال‌های 1990 و پس از آن بر می­گردد. در دهه‌های 1950 و 1960 بیشتر تولید­کنندگان برای کاهش هزینه تولید هر واحد محصول بیشتر روی تولید انبوه به­عنوان استراتژی عملیاتی اولیه تکیه می‌کردند و به انعطاف­پذیری محصول و فرآیند توجهی نداشتند. پیشرفت­ها در زمینه محصول بسیار کند انجام می­شد و منحصراً به تکنولوژی فردی خود سازمان و ظرفیت آن متکی بود. در دهه 1970 برنامه‌ریزی منابع تولید[5] (MRP) مطرح شد و مدیران به مفاهیم جدیدی در مدیریت مواد برای بهبود عملکرد سازمان متوسل شدند. رقابت جهانی در دهه 1980 باعث شد تا سازمان­های بزرگ قیمت­ها (هزینه­ها) را پایین بیاورند، کیفیت را بالا برده، محصولات مطمئن را با انعطاف­پذیری در طراحی ایجاد کنند. استفاده از فلسفه تولید بهنگام نیز برای بهبود کارایی تولید و زمان تولید متداول شد. در دهه 1990 با گسترش مدیریت مواد در بحث زنجیره ارزش[6]، تولید کنندگان فقط از تعداد معدودی تأمین کننده معتبر و تأیید شده کالای مورد نیاز خود را تأمین می‌کردند و این کار موجب افزایش کارایی در زنجیره ارزش شده است. تن[7] (2001) تکامل SCM را در دو مسیر خرید و عرضه از دید خریداران، و حمل و نقل و لجستیک از دید خرده فروشان / عمده فروشان تقسیم­بندی می­کند.

کو و همکاران[8] (1999) شبکه توزیعی شامل یک انبار مرکزی و تعدادی تأمین­کننده را کرده‌اند. میزان تقاضای انبار غیر قطعی است و از یک سیاست تغییر یافته از مرور دوره­ای برای کنترل موجودی‌ها استفاده می‌شود. کورپلا و لووسمرا[9] (1999) برای طراحی بهینه شبکه انبارها، ابتدا با استفاده از روش فرایند تحلیلی سلسله­ مراتبی[10] (AHP) انبارها را بر اساس سطح خدمت مورد نیاز مشتریان رتبه‌بندی کرده­اند. سپس مدل ریاضی طراحی کرده­اند که در آن مقادیر قطعی تقاضا به­عنوان ورودی هستند و محدودیت‌هایی چون حداقل موجودی مجاز برای انبارهای باز و تعداد انبارها وجود دارد. صبری و بیمون[11] (2000) به بررسی شبکه زنجیره تأمینی شامل تأمین­کنندگان، کارخانجات، مراکز توزیع و نواحی تقاضا پرداخته‌اند. در سطح استراتژیک ساختار بهینه شبکه تعیین می‌شود. مسأله چند محصوله، با تقاضای غیر قطعی و شامل چند ماده اولیه است. چان و همکاران[12] (2005) مسأله­ای را برای پیدا کردن تورهای جدا از هم برای وسایل نقلیه برای توزیع کالا از انبارها به نقاط تقاضا با مقادیر احتمالی طراحی کرده‌اند. پیرکول و جایرمن[13] (2001) زنجیره تأمینی شامل تأمین کنندگان، کارخانجات، مراکز توزیع و نواحی تقاضا را با چند محصول با تقاضای قطعی بررسی می­کنند و مدلی برای تصمیم گیری همزمان در دو سطح استراتژیک و عملیاتی برای راه‌اندازی تسهیلات، تعیین میزان تولید، تخصیص مراکز توزیع به نواحی تقاضا و میزان حمل ارائه داده‌اند.

دسکی و ورتر[14] (2001) به مقوله راه‌اندازی تعدادی تسهیلات جید در یک فضا می‌پردازند؛ به‌گونه­ای که هر یک از تسهیلات به یک محدوده تقاضا سرویس بدهد. نوزیک[15] (2001) مکان‌یابی تسهیلات را در بین نقاط تقاضا با مقادیر قطعی بررسی کرده است. وانگ[16] (2003) با در نظر گرفتن سطوح خدمت مورد نیاز، به طراحی یک سیستم لجستیکی شامل تعدادی مراکز تولید، انبارها یا مراکز توزیع و مشتریان با مقادیر تقاضای غیر قطعی می‌پردازد. ساریف و همکاران[17](2002) شبکه زنجیره تأمینی شامل تأمین­کنندگان، کارخانجات و مراکز توزیع را بررسی کرد‌ه‌اند. تصمیماتی که باید در این شبکه اتخاذ شوند، شامل راه‌اندازی کارخانجات و مراکز توزیع، میزان تولید و حمل محصولات است.

ژو و همکاران[18] (2002) شبکه زنجیره تأمینی را طراحی کرده­اند که هزینه حمل و نقل و سطح خدمت را به بهترین وجه متعادل کند، به گونه­ای که تا حد امکان به همه مراکز توزیع بار کاری یکسان داده شود و کل مسافت حمل شبکه کمینه شود. سیام[19] (2002) مدل زنجیره تأمینی شامل کارخانجات و انبارها با مقادیر تقاضای قطعی مشتریان از چند محصول را ارائه داده است. میراندا و گاریدو[20] (2004) تصمیمات کنترل موجودی مثل اندازه سفارش اقتصادی، نقطه سفارش و موجودی اطمینان را همزمان با مکان‌یابی تسهیلات در طراحی شبکه توزیع اتخاذ کرده‌اند.

1-4- مروری بر مقالات

در این بخش، یک دسته­بندی مقالات مورد بررسی به­صورت خلاصه در جدول (1) ارایه می­گردد. باید اشاره کرد که در مورد اهداف، خروجی­ها، محدودیت‌ها و مفروضات مدل‌های مربوطه در ادامه مباحث توضیحاتی کوتاه آورده شده است.

الف) مقالات مطالعه شده را می­توان از لحاظ پرداختن به سطوح برنامه­ریزی در دو سطح استراتژیک و عملیاتی قرار داد.

ب) از نظر ماهیت مقادیر تقاضا، مقالات ارایه شده را می­توان به دو دسته قطعی و غیر قطعی طبقه­بندی کرد که بیشتر آنها مقادیر قطعی تقاضا را بررسی کرده‌اند و تعداد کمتری به مقادیر غیر قطعی و احتمالی پرداخته­اند.

ج) از نظر محدودیت‌ها باید گفت از جمله محدودیت‌های عمده به کار رفته در مدل‌سازی مسایل توزیع در زنجیره تأمین، محدودیت ظرفیت و سطح خدمت انبارهاست.

د) توابع هدف زیر به صورت تکی یا چندگانه در مدل­سازی­ها استفاده شده‌اند. به طوریکه حداقل هزینه بسیار مورد توجه بوده است. برقراری تعادل بین رده‌ها در موارد کمی مورد توجه بوده است.

ه) رده‌های مورد توجه شامل تدارک، تولید و توزیع هستند. در اغلب موارد به دو تا از این سه رده توجه شده است؛ ضمناً تدارک کمتر مورد توجه قرار گرفته و بیشتر تولید - توزیع تحت مطالعه بوده است.

و) بیشتر کارهای انجام شده روی یک محصول بوده است تا چند محصول.

ز) تعداد تحقیقات در حالت تک دوره‌ای بیشتر است؛ لیکن به نظر می‌رسد روند حرکت به سوی بیشتر منظور نمودن چند دوره است.

ح) در زمینه ارایه مدل ظاهراً بیان این مسایل در درجه اول به عنوان یک مدل ریاضی همیشه مورد توجه بوده است. این مدل‌ها اغلب به صورت خطی ولی با متغیرهای مختلط هستند.

ط) اغلب این مدل­های دارای اندازه‌های بزرگ و غیرقابل حل به صورت بهینه خواهند بود. به همین دلیل، از نگرش سلسه مراتبی برای حل و گاه در مدل‌سازی استفاده می‌شود. در حل این مسایل اغلب از روش‌های ابتکاری یا فراابتکاری برای حل استفاده شده است.

جدول 1: دسته­ بندی مقالات

2- مدل ریاضی و روش تحقیق

وانگ و همکاران (2003) شبکه توزیعی متشکل از انبارها و خرده فروشان را در نظر گرفته­اند. در این شبکه کالاها قبل از رسیدن به دست مشتری نهایی به طور موقت در انبارها ذخیره می شوند. هر انبار می‌تواند بسته به نیاز خرده فروشان به چند خرده فروش سرویس بدهد. سیستمی که بتواند با فرض محدودیت ظرفیت انبارها در تأمین هر کالا در هر دوره، تقاضای مورد نیاز مشتری را به موقع تأمین کند، بسیاراهمیت خواهد داشت. درصورتی که تأمین نیاز زودتر یا دیرتر از زمان مقرر انجام گیرد هزینه‌هایی را متناسب با قیمت محصول تأمین شده از هر انبار (قیمت تولید در کارخانه) در بر خواهد داشت. خروجی مدل طراحی شده، مقدار موجودی هر محصول در هر دوره در انبار کارخانجات و تابع هدف، کمینه کردن هزینه‌ها است. این مدل به یک مدل ‌برنامه‌ریزی خطی تبدیل شده است که با روش سیمپلکس قابل حل خواهد بود. ایده اولیه ایجاد مدلی که در این مقاله ارائه شده بر مبنای مدل وانگ بوده است. شباهت این مدل با مدل ارائه شده، توسط وانگ در ایده تحویل به موقع است. تفاوت‌های عمده این دو مدل را می توان در موارد زیر خلاصه کرد:

الف) در مدل وانگ و همکاران (2003) تنها حداقل کردن هزینه­ها مورد توجه بوده است، در صورتیکه در این مدل بحث هزینه ها و تحویل بموقع به­عنوان دو مقوله متفاوت مد نظر قرار گرفته است.

ب) مدل وانگ و همکاران (2003) مربوط به یک زنجیره یک رده­ای است، ولی این مدل به صورت زنجیره تأمین مخلوط دو رده و سه رده­ای با دو کانال توزیع است.

شایان ذکر است که روش تحقیق به صورت توسعه الگوریتمی بوده و برای جمع­آوری داده­ها از اندازه­ها و داده­های مسائل حل شده توسط دیگران از جمله وانگ استفاده شده است.

2-1- ویژگی­های مدل

مدل مربوط ویژگی­های خاصی دارد که از جمله می­توان به دو هدفه بودن آن اشاره کرد. بحث حل مسائل چند هدفه از مباحث به­روز و البته، دشوار ریاضی است که در اکثر موارد به بررسی زیاد و استفاده از روش­های غیر ریاضی مانند الگوریتم­های فراابتکاری نیاز دارد که این موضوع به علت ساختار مسأله و افزایش زمان حل آن است. همچنین، در این مدل محدودیت­های زیادی مورد توجه قرار گرفته است که همگی سعی بر این دارند که کالاها در بهترین زمان انتقال داده شوند و در عین حال، حداقل کردن هزینه­ها نیز مورد توجه قرار بگیرد.

2-2- اندیس­ها، پارامترها و متغیرهای تصمیم مدل

در اینجا، ابتدا مدل ریاضی پیشنهادی و سپس توضیحی در مورد مشخصات آن ارایه می­شود. این مدل دارای پنج اندیس تأمین کنندگان (i)، عمده فروشان (j)، خرده فروشان (k)، دوره­ها (t) و کالاها (p) است.

متغیرهای مدل شامل موارد زیر است:

y_{pijt :} مقدار کالای p که در ابتدای دوره t از تأمین کننده i به عمده فروش j حمل می‌شود.

u_pkjt: مقدار کالای p که در ابتدای دوره t از عمده فروش j به خرده فروشk حمل می‌شود.

v_pikt: مقدار کالای p که در ابتدای دوره t از تأمین کننده i به خرده فروش k حمل می‌شود.

B_pkt: مقدار کمبود خرده فروش k از کالای p در دورهt.

In_pkt: مقدار اضافی تحویل داده شده از کالای p به خرده فروش k در دوره t .

v'_pkt: متغیر صفر و یک نشان­دهنده وجود کمبود یا اضافی در مقدار تحویل داده شده از کالای p  به خرده­فروش k در دورهt .

 پارامترهای مدل نیز شامل موارد زیر است:

a_ijp: مدت زمان تحویل کالای p از تأمین کننده i به عمده فروش j .

f_ikp: مدت زمان تحویل کالای p از تأمین کننده i به خرده فروش  k.

S_pit: ظرفیت تولید کالای p توسط تأمین کننده i در دوره t .

d_pkt: میزان تقاضای کالای p در خرده فروشی k در ابتدای دوره t .

c_p: هزینه حمل واحد کالای p به­ازای واحد مسافت.

cu_pi : قیمت‌خرید هر واحد کالای p از تأمین‌کنندهi.

bl_pkt : حداکثر میزان کمبود مجاز کالای p در انبار خرده فروش kدر دوره t.

m_ij: مسافت بین تأمین کننده i و عمده فروش j .

n_ik: مسافت بین تأمین کنندهi و خرده فروش k .

o_jk: مسافت بین عمده فروش j  و خرده فروش k .

h_pj: هزینه نگهداری هر واحد کالای p در انبار عمده فروش j در هر پریود.

h'_pk: هزینه نگهداری هر واحد کالای p در انبار خرده فروش k در هر پریود.

Q_pjt: ظرفیت نگهداری انبار عمده فروش j  برای کالای p در پریود t.

Q'_pkt: ظرفیت نگهداری انبار خرده فروش k برای کالای p در پریود t.

Ca_jt: ظرفیت تحویل­گیری کالا در انبار عمده فروش jدر دوره t.

Ca'_kt: ظرفیت تحویل­گیری کالا در انبار خرده فروش k در دوره t.

2-3- مدل ریاضی پیشنهادی

بر اساس پارامترهای فوق، مدل دو هدفه ذیل برای مسأله ارایه می‌شود:

(1--3)

رابطه (3-1) مجموع هزینه­های خرید عمده­فروشان و خرده­فروشان از تأمین­کنندگان و هزینه­های حمل کالاها و هزینه نگهداری کالاها در انبار عمده فروشان و خرده فروشان را حداقل می کند. رابطه (3-2) به عنوان تابع هدف تحویل بموقع، مجموع مدت زمان ضربدر مقادیر دیرکردها و زودکردهای تحویل کالاها به انبار خرده فروشان را در کل افق برنامه‌ریزی حداقل می­کند.

مجموعه محدودیت­های (3-3) تضمین می­کند که مقدار سفارش کالا از عمده فروشان و خرده فروشان به تأمین­کنندگان در هر دوره از ظرفیت تأمین کننده تجاوز نمی­کند. مجموعه محدودیت­های (3-4) تأمین کامل تقاضاهای مشتریان را نشان می­دهد، بدون اینکه میزان موجودی اضافی در انتهای افق برنامه­ریزی در انبار خرده فروش وجود داشته باشد؛ یعنی میزان کل ورودی به انبار خرده فروشان در طول افق برنامه‌ریزی باید برابر با کل تقاضای مشتری در طول افق برنامه‌ریزی باشد. مجموعه محدودیت­های (3-5) محدودیت ظرفیت نگهداری انبارهای عمده فروشان را در هر دوره بیان می­کند؛ یعنی تفاضل خروجی‌ها از ورودی­ها در هر دوره باید در حد ظرفیت انبار باشد.

مجموعه محدودیت­های (3-6) بیان می کند که میزان خروجی از انبار عمده فروشان در هر دوره نمی‌تواند از میزان موجودی انبار بیشتر باشد.  مجموعه محدودیت­های (3-7) برابری کل ورودی­ها به انبار عمده فروشان را با کل خروجی‌ها از آن در طول افق برنامه‌ریزی بیان می کند. یعنی در انتها موجودی انبارهای عمده فروشان باید صفر باشد. مجموعه محدودیت­های (3-8) و (3-9) بیانگر ظرفیت تحویل­گیری کالاها در انبار­های خرده فروشی و عمده فروشی هستند.

مجموعه روابط (3-10) میزان موجودی و یا کمبود کالا را در انبار خرده فروشی در هر دوره به دست می‌آورند. مجموعه محدودیت­های (3-11) و (3-2) نشان می‌دهند که در هر دوره در هر خرده فروشی برای هر کالا فقط یکی از حالت‌های کمبود یا موجودی اضافی می­تواند وجود داشته باشد. به علاوه؛ حداکثر میزان موجودی و کمبود مجاز هر دوره در انبارهای خرده فروشی را نشان می­دهند.

   مجموعه محدودیت­های (3-13) و (3-14) عنوان می­کند که به ازای اندیس­های خارج از افق برنامه‌ریزی مقادیر حمل صفر خواهد بود. مجموعه محدودیت­های (3-15) صفر و یک بودن متغیرهای نشانگر موجودی اضافی یا کمبود در تحویل را نشان می‌دهد. مجموعه محدودیت­های (3-16) مثبت بودن مقادیر حمل، کمبود و موجودی اضافی را نشان می‌دهد.

2-4- روش حل پیشنهادی

ابتدا توضیح مختصری در مورد الگوریتم بهینه­سازی ذرات انبوه[21] (PSO) که مورد استفاده در مدل پیشنهادی بوده است، ارایه می‌شود. این الگوریتم به خانواده جدیدی از الگوریتم­هایی که برای یافتن جواب‌های بهینه و نزدیک به بهینه در مسایل عددی به کار می­روند، اطلاق می­گردد. این روش یک الگوریتم بهینه­سازی احتمالی و بر اساس جمعیت بوده که پس از شبیه­سازی رفتار اجتماعی گروه پرندگان مدل­سازی گردید. PSO در آغاز توسط کندی و ابرهارت (1995) بر اساس الگوریتم­هایی که رفتار دسته جمعی در آنها مانند پرندگان مشاهده شده، توسعه داده شد. این الگوریتم از این جهت که یک الگوریتم بر اساس جمعیت بوده و برای ارزیابی هر جزو (جواب) از تابع برازندگی استفاده می­کند، شبیه به الگوریتم­های تکاملی[22] است. اما تفاوت عمده­ای که وجود دارد، این است که در الگوریتم بهینه­سازی ذرات انبوه، هر جواب از حافظه خود بهره­مند می‌شود و این در مورد الگوریتم­های تکاملی صادق نیست.

شایان ذکر است که برای تعیین میزان شایستگی[23] جواب‌ها از روشی بر پایه بهینگی پاراتو به صورت رده­بندی جواب‌های غیرمغلوب برای هدایت جامعه به سمت یک جبهه پاراتو استفاده می‌شود. یک بردار تصمیم  بهینه پاراتو است اگر و تنها اگر هیچ  یافت نشود که  بر غلبه کند. بدین ترتیب، با هر بار حل مسأله مجموعه‌ای از جواب‌های قابل قبول به دست می­آیند که هیچ یک به ­طور مطلق برتر از دیگری نیستند. اصطلاحاً این جواب‌ها را جواب‌های غیرمغلوب می‌گویند. انتخاب یکی از این جواب‌ها از میان سایرین به عواملی بستگی دارد که به ماهیت مسأله وابسته هستند. بنابراین، با تغییر شرایط و یا حتی از نظر تصمیم گیرندگان مختلف، جواب مؤثر انتخاب شده یکسان نخواهد بود.

در ادامه، پس از ارزیابی پارامترهای ذرات و تعیین مقادیر توابع هدف به ازای حل­های شدنی با استفاده از مفهوم بهینگی پاراتو ذرات جامعه را رده­بندی می‌کنیم؛ به این صورت که ابتدا از بین تمامی ذرات، ذره­های غیرمغلوب را در رده اول قرار می دهیم. با کنار گذاشتن ذرات رده اول، ذره­های غیر مغلوب را ازبین بقیه ذرات و با توجه به مخزن جدا کرده، در رده پسی قرار می دهیم. به همین ترتیب ادامه می‌دهیم تا تمام ذره­ها رده‌بندی شوند. در شکل (1) مثالی در سه سطح رده بندی شده اند. ذره­های سطح اول بیشترین مقدار شایستگی را دارند و به ترتیب از این میزان کاسته می‌شود. در واقع شایستگی در ازای نمودار تراز جبهه پاراتوی دو تابع هدف بررسی می‌شود که هرچه این خطوط تراز به مبدا نزدیک‌تر باشند؛ یعنی این­که دو جبهه به­هم نزدیک­تر و توزیع یکنواخت­تری را ارائه می­دهند.

شکل 1: مثالی از رده­بندی نقاط غیر مغلوب به­صورتی­

که: fitness₁=fitness₃>fitness₅>fitness₄

2-5- طراحی الگوریتم بهینه­سازی ذرات انبوه

در این بخش، رویکردی که برای بهینه­سازی ذرات انبوه طراحی شده ارایه می­گردد. در این مقاله، برای الگوریتم بهینه­سازی ذرات انبوه چندهدفه از روش کوئلو[24] (2004) و روش ابتکاری بالتر و فوتان[25] (2010) استفاده می­شود.

الف) رویکرد حل مسأله

مطابق با نظر کندی[26] (1999) در مدل g_best[27] هر یک از اجزا نسبت به بهترین موقعیت قبلی­اش و نسبت به بهترین موقعیتی که در کل دسته وجود دارد حرکت می­کند، اما در حالت در نظر گرفتن چند هدف بطور همزمان، چیزی به نام بهترین موقعیت وجود ندارد و مجموعه جواب‌های بهینه پاراتو به عنوان جواب‌های بهینه وجود دارند و باید از روش‌های بهینه­سازی چندهدفه استفاده نمود. باید توجه داشت، در الگوریتم بهینه­سازی ذرات انبوه استاندارد ابتدا مجموعه ذرات[28] تولید می‌شود. این آغازسازی شامل بردارهای مکان و سرعت خواهد بود. p_best برای هر ذره محاسبه و leader در آن مکان می­یابد. سپس برای تعداد مشخصی از تکرار، بردارهای سرعت و مکان بروز می­شوند که در نتیجه آن p_best[29] و پیشروها[30] نیز به­روز خواهند شد، اما با کمی تغییر برای حل مسایل چند هدفه با استفاده از الگوریتم بهینه­سازی ذرات انبوه باید پس از تولید جمعیت اولیه یا swarm ، یک مجموعه از leader ها از جمعیت مذکور که ذرات غیرمغلوب را ایجاد می‌کنند، تولید شوند. می­دانیم که این مجموعه در مخزن نگهداری می‌شود. در هر بار تولید یا به اصطلاح نسل برای هر ذره یک leader  انتخاب و یک حرکت یا پرواز صورت می­گیرد. اکثر الگوریتم‌های ذرات انبوه چند هدفه[31]از تنظیمات اپراتور جهش پس از انجام پرواز بهره می­برند. در ادامه، یک ذره مورد ارزیابی قرار می­گیرد و مطابق با p_best آن به­روز می­گردد، به­خصوص زمانی که این ذره مغلوب باشد. پس از اینکه همه ذرات بروز شدند، مجموعه leader ها نیز بروز می­شوند. این فرایند تا تکرار معینی ادامه پیدا می­کند. در ادامه، الگوریتم بهینه­سازی ذرات انبوه برای مدل پیشنهادی تشریح می­گردد.

ب) طراحی اجزا یا آرایه مکان هر پرنده

با توجه به مدل ریاضی ارایه شده، مشاهده می­گردد که متغیر اصلی که می­توان طراحی اجزا را بر اساس آن انجام داد، متغیرv'_pkt است. کلیه متغیرهای دیگر مدل بر اساس این متغیر به دست خواهند آمد. با توجه به مدل ریاضی مسأله، مقدار v'_pktدر صورتی می­تواند برابر 1 باشد که با توجه به محدودیت­های (3-11) متغیر In_pktبتواند مقدار مثبت بگیرد؛ بدین مفهوم که موجودی کالای p در خرده فروشی k در دوره t مقداری مثبت باشد. در صورتی می­توانیم موجودی داشته باشیم که با توجه به محدودیت‌های (3-13) و (3-14) حداقل به اندازه کمترین زمان تحویل ممکن از ابتدای افق برنامه‌ریزی گذشته باشد و به عبارتی:

بنابراین، متغیرها یا ذراتی می­توانند مقدار 1 بگیرند که در رابطه (1) صدق کنند.

برای ایجاد یک رابطه مستقیم میان مسأله و ذرات الگوریتم بهینه­سازی ذرات انبوه برای این مسأله، هر ذره به منزله یک جواب در نظر گرفته می‌شود و دارای p  پس است که به منزله در نظر گرفتن  pنوع کالاست. بدین ترتیب، هر بعد نشان دهنده یک کالاست و ذره به صورت  نمایش داده می‌شود. هربعد از هر ذره خود دارای kپس دیگر است که نشان­دهنده تعداد kخرده فروش برای کالای نوع p است و به صورت  نمایش داده می‌شود. بدین ترتیب، هر سطر ممکن است دارای تعداد ستون­های متفاوتی باشد و در کل هر ذره دارای kp بعد است.

ج) تعریف پارامترهای اولیه و اصطلاحات به­کار رفته

پارامترهای موجود در این الگوریتم به شرح زیر است:

T: این پارامتر تعیین کننده تعداد کل تکرار است.  NP: این پارامتر نشان دهنده تعداد پرندگان و یا اجزای موجود در هر دسته است که این مقدار بر اساس تحقیق تاسگترین[32] (2007) می­تواند دو برابر ابعاد هر ذره باشد. با توجه به اینکه در مسأله ما ابعاد هر ذره برابر با مجموع کالاهای مختلف برای خرده فروش­ها ست، مقدار NP به­ صورت  محاسبه می‌شود که در آن هر خرده فروش می­تواند هر تعداد از p کالا را دریافت نماید.

: حداقل مقداری ­است که مکان هر پرنده  آن را تجربه می­کند که به طور کلاسیک این مقدار برابر با صفر است.

: حداکثر مقداری ­است که مکان هر پرنده  آن را تجربه کند که به­طور کلاسیک این مقدار برابر با عدد 4 است.

: این مقدار که همان مقدار بیشینه سرعت       است، معمولاً برابر با است. همچنین، مقدار  .

: این پارامترها که به ضرایب شتاب دهنده معروف بوده، به ترتیب، پارامتر ادراکی و پارامتر اجتماعی نامیده می­شوند، مطابق با نظر کندی و ابرهارت می­توانند در بازه (4و1) قرار گیرند که معمولاً مقادیر  برابر 2 اختیار می­گردد (البته، نشان داده شده است که بهتر است 4c₁+c₂≤ ).

: این پارامتر ضریب اینرسی است که با توجه به فرمول زیر محاسبه می­گردد و با توجه به مقالات مختلف مقادیر مناسب برای  و برای  است. در فرمول زیر مقدار t  نشان­دهنده تعداد تکرارهایی است که تاکنون صورت پذیرفته است و مقدار T نیز نشان­دهنده تعداد کل تکرارها است.

بطور کلی با توجه به نظر کندی و ابرهارت (1995) مقادیر پارامترهای اولیه تعریف شده در بالا می­توانند در بازه­های زیر قرار گیرند:

د) تولید آرایه مکان و سرعت اولیه برای ذرات انبوه

شمارنده تکرار برابر صفر در نظر گرفته می‌شود یعنی. برای تولید آرایه اجزای اولیه، از یک رویکرد ساده تصادفی بر حسب توزیع یکنواخت استفاده شده است. به این منظور، یک جمعیت از ذرات به تعدادNP  آرایه مکان که نشان دهنده مجموعه جواب مسأله هستند، ایجاد شده، خانه­های آرایه اجزا با استفاده از فرمول زیر به صورت تصادفی با اعداد بین صفر و پر می­گردد:

 یک عدد تصادفی است که از توزیع یکنواخت (1و0) پیروی می­کند. بنابراین، به تعداد NP ذره {, i=1,2,...,NP} خواهیم داشت.

سرعت اولیه نیز از فرمول مشابهی برای کلیه مکان­ها و به صورت تصادفی از فرمول زیر ایجاد می‌شود:

 یک عدد تصادفی است که از توزیع یکنواخت (1و0) پیروی می­کند. بنابراین، به تعداد ذرات ایجاد شده، سرعت اولیه {, i=1,2,...,NP} خواهیم داشت.

ه) آزمون امکان­پذیری

از آنجایی که این مسأله دارای محدودیت­های خاصی است، بنابراین یک مرحله با عنوان آزمون امکان پذیری جواب به دست آمده از مرحله قبل در نظر گرفته شده است که محدودیت­های موجود از قبیل ظرفیت نگهداری انبارها و ورود و خروج از آنها و زمانبندی مربوط به زودکرد و دیرکرد تحویل کالاها در این مرحله بررسی می‌شوند. ذراتی که در این مرحله قابل قبول شناخته شوند، به مرحله پس و تعیین تابع برازندگی آنها خواهند رفت و ذراتی که در این آزمون قابل قبول شناخته نشوند، به مرحله قبل باز خواهند گشت تا مجدداً تولید شده و رویه را طی کنند. این فرآیند تا زمانی ادامه خواهد داشت که به تعداد NP ذره قابل قبول به وجود آید.

و) ارزیابی جواب‌ها و تابع برازش

در حالت تک هدفه بودن مدل، الگوریتم بهینه‌سازی ذرات انبوه پس از تعیین آرایه زمانبندی، مقدار آن تابع هدف برای همه اجزا محاسبه می­گردد و سپس با توجه به مدل g_best بهترین مکان برای هر جزو و برای کل اجزا در هر تکرار به دست می­آید.  با توجه به دوهدفه بودن مدل پیشنهادی که در قسمت قبل ارایه شده است، نمی­توان یک بهترین مکان را محاسبه نمود و در عوض، مجموعه جواب‌های بهینه پاراتو به عنوان جواب‌های بهینه وجود دارند، بدین ترتیب مطابق با روش کوئلو و همکاران (2004) عملکرد ذرات مختلف متناسب با مغلوبیت آنها سنجیده می‌شود و یک مخزن برای ذخیره ذرات غیرمغلوب در نظر گرفته شده است.

ز) به­روز رسانی رویه

به شمارنده تکرار یک عدد اضافه می‌شود و . ضریب اینرسی با توجه به فرمول گفته شده در قسمت ج بروز می‌شود. سرعت هر آرایه به صورت فرمول زیر به روز رسانی می‌شود:

 یک جوابی است که از مخزن خارجی در هر تکرار  انتخاب شده است. این انتخاب در حالتی انجام می‌شود که جواب‌های غیرمغلوب در مناطقی از فضای هدف که دارای تراکم بیشتر جمعیت است، واقع شده­اند و دارای احتمال کمی برای انتخاب هستند. مطابق با روش بالتر و فوتان (2010)، در فضای هدف، تراکم نقاط اطراف هر جواب ذخیره شده در مخزن را محاسبه می‌شود و چرخ رولت به عنوان تابع انتخاب جواب با تراکم جمعیت رابطه معکوس دارد. تراکم جمعیت از روی تعداد ذراتی که از لحاظ یکی از توابع هدف به یکدیگر نزدیک هستند، سنجیده می‌شود.

با توجه به سرعت به دست آمده، آرایه­های مکان به صورت فرمول زیر به روز می­شوند:

همانند رویه­های توضیح داده شده در قسمت‌های (ه- ز) آزمون امکان­پذیری انجام شده و ارزیابی جواب‌ها و تابع برازش به دست می­آید. اگر ذره یا ذراتی با توجه به آزمون امکان­پذیری قابل قبول نباشند سرعت آنها مجدداً بروز رسانی شده و آرایه مکان آنها مجدداً محاسبه خواهد شد و رویه تا زمانی ادامه پیدا می­کند که کلیه ذرات قابل قبول باشند. مطابق با روش بالتر و فوتان  در فضای هدف، تراکم نقاط اطراف هر جواب ذخیره شده در مخزن محاسبه می‌شود و چرخ رولت به­عنوان تابع انتخاب جواب با تراکم جمعیت رابطه معکوس دارد. تراکم جمعیت از روی تعداد ذراتی که از لحاظ یکی از توابع هدف به یکدیگر نزدیک هستند، سنجیده می‌شود.

در هر تکرار t ، موقعیت جدید هر ذره با تمام ذرات مقایسه می‌شود و غیرمغلوب­ها با کلیه جواب‌های ذخیره شده در مخزن مقایسه می­شوند. سپس مخزن به­هنگام می‌شود و جواب جدید غیرمغلوب را ذخیره کرده و جواب‌های قدیمی که اکنون مغلوب هستند، حذف می­شوند. هنگامی که مخزن پر باشد و جواب غیرمغلوب یافته شود، در این صورت این جواب جایگزین یکی از جواب‌های غیرمغلوب دیگر در مخزن می‌شود.

ح) معیار توقف

هنگامی که تعداد تکرارهای الگوریتم از مقدار تعیین شده فراتر رود، رویه الگوریتم متوقف می‌شود و جواب به دست آمده، بهترین جواب حاصل از اجرای الگوریتم در همه تکرارها و از همه ذرات و جواب نهایی است.

ط) جستجوی موضعی برای الگوریتم پیشنهادی

در این الگوریتم پس از به­دست آوردن بهترین مکان برای کل اجزا در هر تکرار از یک الگوریتم جستجوی موضعی برای یافتن بهترین مکان­های همسایگی نیز استفاده می­گردد. مکان­های همسایگی به صورت زیر تولید می­گردد:

1- جابه‌جایی دوره­های مربوط به هر آرایه.

2- استفاده از عملگر تقاطع ستونی که در آن ابتدا تصادفاً یک نقطه برش در برای ستونی به تصادف در آرایه انتخاب می‌شود و سپس مطابق شکل بخش­های متناظر آرایه با یکدیگر تعویض می­گردند.

پس از اجرای جستجوی موضعی، مقادیر جدیدی برای ذرات تولید خواهد شد که پس از تبدیل آنها به آرایه مورد نظر، این ذرات جدید نیز نیاز به عبور از آزمون امکان­پذیری­ دارند. چنانچه ذرات توسط این آزمون مورد قبول باشند، به مرحله تعیین تابع برازندگی و ادامه رویه می­روند؛ در غیر این صورت تقاطع ستونی دیگری برای جستجوی موضعی انجام می‌شود و این رویه تا زمانی ادامه خواهد داشت که ذرات حاصل از جستجوی موضعی همگی قابل قبول باشند.

3- تحلیل نتایج محاسباتی

به منظور حل مسأله مورد نظر، پس از تعریف آرایه مکان و پارامترها و الگوریتم جستجوی موضعی، در گام اول پارامترهای ورودی برای مسأله مورد نظر تعریف می­گردند. ورودی­های مسأله برای بکارگیری الگوریتم بهینه­سازی ذرات انبوه ارایه شده و روش شبیه­سازی به شرح زیر است:

 (پارامتر ادراکی و پارامتر اجتماعی) -   و (ضریب اینرسی) - (حداقل مقدار آرایه مکان هر جز) -( حداکثر مقدار آرایه مکان هر جز) - (مقدار بیشینه سرعت) - (تعداد حل­های شدنی ایجاد شده در هر گروه) - (تعداد تکرار برای رسیدن به حل بهتر) -  (تعداد انواع کالاها) -  (تعداد خرده­ فروش­ها) (تعداد انواع تأمین­کنندگان) -  (تعداد انواع عمده­ فروش­ها)

بر اساس اطلاعات ورودی، کد مورد نظر در برنامه ویژوال بیسیک  اجرا می‌شود که به شرح زیر است:

در ابتدا کلیه اطلاعات ورودی فرا خوانده می‌شود، آغازسازی اولیه مطابق با الگوریتم بهینه­سازی ذرات انبوه چندهدفه ایجاد خواهد شد. برای کل اجزای دسته، آرایه مکان و همچنین، آرایه مربوط به سرعت اولیه به صورت تصادفی تولید می­گردند. سپس آرایه‌های مکان ایجاد شده برای کل اجزای دسته مطابق با شرحی که در قسمت قبل داده شده، مورد آزمون امکان­پذیری قرار می­گیرند و اگر امکان­پذیر نباشند، فرآیند تولید ذرات مجدداً رخ می­دهد. سپس مطابق با فرمول­های ارایه شده در بخش قبل آرایه­های سرعت و مکان هر یک از اجزای دسته به­روز می‌گردد و رویه­ای که در بالا شرح داده شد، تا زمان برآورده شدن معیار توقف، که تعداد تکرارها به عنوان معیار توقف در نظر گرفته شده است، ادامه پیدا می‌کند؛ بدین ترتیب، یک حل اولیه خواهیم داشت. در نهایت، خروجی مورد نظر که به­شرح زیر خواهد بود، ارایه می‌شود.

برای مقایسه جواب‌های غیرمغلوب به­دست آمده توسط الگوریتم بهینه­سازی ذرات انبوه با نرم­افزار لینگو، مقدار تابع هدف Z₂ مربوط به هر جواب به صورت محدودیتی وارد مدل ریاضی اصلی شده و مقدار حداقل تابع هدف Z₁ در مدل جدید شامل محدودیت­های اصلی و محدودیت اضافه شده محاسبه می‌شود. همچنین، برای مقایسه کیفیت جواب‌ها در این دو روش نیز معیارهایی بر اساس اختلاف مقادیر متناظر تابع هدف Z₁ در دو روش به ازای Z₂ های یکسان تعریف شده است که در واقع نشان دهنده فاصله دو جبهه پاراتو به ­دست آمده از دو روش هستند. برای ارزیابی کیفیت الگوریتم طراحی شده، مسایل آزمایشی در ابعاد مختلف توسط الگوریتم حل شده‌اند. در ابعاد کوچک این مسایل توسط هر دو روش لینگو و الگوریتم بهینه­سازی ذرات انبوه طراحی شده حل شده و نتایج آنها با هم مقایسه شده است. ابعاد مسایل آزمایشی طراحی شده در جدول (2) آورده شده است.

شایان ذکر است به ­عنوان مرجعی برای تعیین اندازه مسائل آزمایشی، اندازه مسائل حل شده در ادبیات مسائل توزیع در مدیریت زنجیره تأمین از نظر پارامترهای مربوط به تعداد تسهیلات هر رده و یا تعداد انواع محصولات بررسی شد.

از آنجا که از یک طرف مقادیر پارامترهای مربوط به ظرفیت‌ها و مقادیر تقاضا در مدل در مشخصات جواب‌های به دست آمده شامل چگونگی انتخاب کانال‌ها و مقادیر توزیع و یا نگهداری کالاها در انبارها مؤثر است و از طرف دیگر، ترکیب مقادیر پارامترهای عنوان شده با پارامترهای مربوط به هزینه‌های شبکه توزیع در مقادیر توابع هدف نسبت به یکدیگر تأثیرگذار است، در هر اندازه، چندین مسأله با ترکیب‌های متفاوت از مقادیر این پارامترها مطابق جدول (3) و به­صورت تصادفی بین حدود بالا و پایین طراحی شده­اند تا به گونه­ای شرایط مختلف دنیای واقعی شبیه سازی شده باشد. نتایج این جدول نشان می­دهد که زمان حل الگوریتم بهینه­سازی ذرات انبوه با افزایش تعداد متغیرهای صفر و یک مدل افزایش می­یابد.

در شکل (2) متوسط زمان صرف شده توسط الگوریتم بهینه­سازی ذرات انبوه برای حل مسایل در ابعاد مختلف بسته به تعداد متغیرهای صفر و یک آنها آورده شده است. به علاوه، با افزایش ابعاد مسأله نسبت زمان حل الگوریتم بهینه­سازی ذرات انبوه به زمان حل لینگو کاهش می یابد که نشان می­دهد با افزایش ابعاد مسأله الگوریتم بهینه­سازی ذرات انبوه جواب‌هایی با کیفیت خوب را در مدت زمانی به مراتب کمتر از روش دقیق به دست می آورد. این موضوع در شکل (3) نشان داده شده است.

جدول 2: ابعاد مسایل آزمایشی طراحی شده

جدول 3: محدوده مقادیر پارامترهای مدل در مسایل آزمایشی

شکل (4) جبهه­های پاراتوی به دست آمده از دو روش را نشان می­دهد. این شکل به روشنتر شدن مفهوم معیارهای استفاده شده برای مقایسه دو روش شامل درصد خطای مقادیر Z₁ و درصد انطباق دو جبهه و همچنین، میزان پراکندگی نقاط روی آن کمک می کند. در این شکل نقاط متناظر روی دو لبه که مقادیر Z₁ یکسان دارند با عناوین یکسان و اندیس‌های متفاوت نشان داده شده‌اند. برای محاسبه معیار درصد خطای مقادیر Z₁ در جبهه پارتوی به دست آمده از الگوریتم بهینه­سازی ذرات انبوه، طبق جدول (1) درصد اختلاف مقادیر Z₁ نقاط A₁ تا I₁ از مقادیر Z₁ نقاط متناظر A₂ تا I₂ محاسبه می‌شود. میانگین، حداقل و حداکثر مقادیر محاسبه شده به عنوان معیارهای خطای جبهه پارتوی PSO نسبت به جبهه لینگو قلمداد می‌شوند. نقاطی که مقدار محاسبه شده برای آنها صفر باشد، بر نقاط متناظر خود روی جبهه پارتوی لینگو منطبق هستند. در این مثال، نقاط متناظر (D₁, D₂) و (E₁, E₂) با مقادیر درصد خطای 2/0% تقریباً بر هم منطبق هستند، که در این مثال به ترتیب 3069، 2161، 946، 2724،532، 924، 4647 و 2525 هستند.

ملاحظه می‌کنید که هر چه این اعداد در کل به هم نزدیک‌تر باشند، به این معنی است که نقاط به طور یکنواخت از لحاظ مقادیر Z₂ شان (با فواصل یکسان) در طول جبهه توزیع شده­اند و این بدان معنی است که انحراف معیار این اعداد نسبت به میانگین­شان کمتر خواهد بود. بنابراین، هر چه مقدار نسبت میانگین به انحراف معیار بیشتر باشد، توزیع نقاط یکنواخت‌تر است.

شکل 2: نمودار متوسط زمان حل الگوریتم بهینه‌سازی ذرات انبوه بر حسب تعداد متغیرهای صفر و یک

شکل 3: نمودار مقادیر میانگین، حداقل و حداکثر نسبت زمان حل  PSOبه Lingo در هر اندازه

شکل 4: نمودار دو جبهه به دست آمده از دو روش

 (خط پررنگ مربوط به PSO و خط کم رنگ مربوط به Lingo ارتباط دارند).