الگوریتم جدیدی برای حل مسأله مسیریابی-موجودی با ارسال مستقیم

1- مقدمه

مسأله مسیریابی-موجودی[1] بسط مهمی از مسأله مسیریابی وسیله نقلیه[2] است که در آن تصمیمات کنترل موجودی و مسیریابی در هم ادغام می­شوند (کوردیو و همکاران[3]، 2007). مسأله مسیریابی-موجودی بیشتر در سیستم­های مدیریت موجودی توسط فروشنده[4] (VMI) کاربرد دارد. در سیستم­های مدیریت موجودی توسط فروشنده، فروشنده قادر است تا زمانبندی و اندازه تحویل محصول به خرده­فروشان را کنترل نماید. در قبال این آزادی عمل، فروشنده تضمین می­کند که مشتریان با کمبود مواجه نمی­شوند. در روابط سنتی­تر میان فروشنده و مشتری که در آن مشتریان درخواست سفارش محصولات را به فروشنده می­دادند، به دلیل زمانبندی سفارشات مشتریان، ممکن است کارایی به شدت کاهش و به نوبه آن هزینه­های موجودی و توزیع به شدت افزایش یابد. با وجود این، تحقق کاهش هزینه­های ناشی از به کارگیری سیستم­های VMI در عمل ساده نیست به ویژه با افزایش تعداد و تنوع مشتریان این امر دشوارتر نیز می­شود. با استفاده از مسأله مسیریابی-موجودی دستیابی به این هدف امکانپذیر است. در مسأله مسیریابی-موجودی با تعیین برنامه توزیع بهینه­ای که مجموع هزینه­ها را کمینه سازد، می­توان به این هدف دست یافت (کوردیو و همکاران، 2007). مسأله مسیریابی-موجودی در پژوهش­های متعددی بررسی شده است که مرور جامعی از پژوهش­های پیشین توسط اندرسون و همکاران[5] (2010) ارائه شده است. نویسندگان با بررسی ابعاد صنعتی مسأله، طبقه­بندی و مرور جامعی از پژوهش­های موجود ارائه داده­اند.

از جمله استراتژی­های توزیعی که می­توان در مسأله مسیریابی-موجودی به کار گرفت، استراتژی توزیع ارسال مستقیم[6] است. ارسال مستقیم به حالتی اشاره می­کند که در آن در هر دوره هر وسیله حمل تنها به یک خرده­فروش محصول تحویل می­دهد. به دلیل سادگی پیاده­سازی این استراتژی توزیع در عمل، طبیعی است که در مسأله مسیریابی-موجودی ابتدا چنین استراتژیی بررسی شود. از این رو، استراتژی ارسال مستقیم در مسأله مسیریابی-موجودی، بطور مکرر در سیستم­های توزیع صنعتی استفاده می­شود. کارایی استراتژی توزیع ارسال مستقیم حداقل به بزرگی جذر کوچکترین نرخ بهره­برداری از وسایل حمل است. به بیان دیگر، در صورتی که نرخ تقاضای هر خرده­فروش 90 درصد ظرفیت وسایل حمل ضرب در کران بالای فرکانس تحویل به خرده­فروشان باشد، آنگاه کارایی استراتژی ارسال مستقیم در حدود 86/94 درصد خواهد بود. کارایی ارسال مستقیم در شرایط خاص، حتی به 100 درصد نیز می­رسد. از این رو، در بسیاری از شرایط استفاده از استراتژی ارسال مستقیم توجیه منطقی دارد (لی و همکاران[7]، 2010). مقالات مختلفی با تمرکز بر استراتژی ارسال مستقیم، به بررسی مسأله مسیریابی-موجودی پرداخته­اند که از آن جمله می­توان به موارد زیر اشاره کرد. لی و همکاران (2010) به بررسی تحلیلی استراتژی ارسال مستقیم پرداخته­اند و کارایی ارسال مستقیم را در شرایط مختلف محاسبه کرده­اند. لی و همکاران (2008) مسأله مسیریابی-موجودی را در حالتی در نظر گرفته­اند که تأمین­کننده تنها یک وسیله حمل در اختیار داشته و در هر دوره تنها می­تواند برای یک مشتری موجودی ارسال کند. نویسندگان الگوریتمی ابتکاری برای یافتن توالی شدنی ارائه کرده­اند. کمبل و هاردین[8] (2005) با بررسی حداقل تعداد وسایل حمل مورد نیاز در مسأله مسیریابی-موجودی با ارسال مستقیم، برای حل مسأله، الگوریتمی حریصانه[9] ارائه داده­اند. چنگ و دوران[10] (2004) مدلی برای مسأله مسیریابی-موجودی در زنجیره تأمین جهانی نفت خام ارائه کرده­اند که در آن تقاضای مشتریان و طول زمان سفر، غیر قطعی است. همچنین، تقاضای مشتریان متغیر با زمان (پویا) است. در این مقاله یک سیستم پشتیبانی تصمیم­گیری برای طراحی و کنترل سیستم موجودی و حمل توسعه داده شده است. کلیوگت و همکاران[11] (2002) با بررسی مسأله مسیریابی-موجودی تصادفی با ارسال­های مستقیم و مدلسازی آن به صورت فرآیند تصمیم­گیری مارکوفی زمان گسسته، روشی تخمینی بر مبنای برنامه­ریزی پویا برای حل ارائه کرده­اند. بارنز شوستر و بسک[12] (1997) به ارزیابی کارایی استراتژی ارسال مستقیم در مسأله مسیریابی-موجودی با افق زمانی نامحدود پرداخته­اند که در آن تقاضای خرده­فروشان احتمالی با تابع توزیع معین و کمبود مجاز است.

مسأله مسیریابی وسایل نقلیه در زمره مسایل با پیچیدگی سخت[13] قرار دارد (لنسترا و رینوی[14]، 1981)؛ در نتیجه، مسأله مسیریابی-موجودی که توسعه­ای از مسأله مسیریابی وسایل نقلیه است، یک مسأله با پیچیدگی سخت است. بنابراین، در مقالات بسیاری، با استفاده از الگوریتم­های فراابتکاری، روش­های حل تقریبی برای مسأله توسعه داده شده است. عزیز و معین[15] (2007) مسأله را در حالت چند محصولی، چند دوره­ای با چند تأمین­کننده و یک کارخانه مونتاژ و با هدف کمینه­سازی مجموع هزینه‌های حمل و نگهداری موجودی بررسی کرده­اند. آنان با بررسی دو نحوه نمایش مختلف جواب، یک الگوریتم ژنتیک ترکیبی بر اساس رویکرد نخست تخصیص و سپس مسیریابی ارائه کرده­اند. معین و همکاران[16] (2010) در تحقیق مشابه­ای، الگوریتم ژنتیک ترکیبی بهبودیافته­ای ارائه داده­اند. بودیا و همکاران[17] (2007) در تحقیق خود با بررسی یک مسأله مسیریابی-موجودی که در آن برنامه­ریزی تولید نیز لحاظ شده است، الگوریتم­های حل مبتنی بر رویه جستجوی انطباقی حریصانه تصادفی[18] توسعه داده­اند. بودیا و پرینز[19] (2009) در تحقیق مشابه­ای، ساختاری مبتنی بر الگوریتم ممتیک[20] با مدیریت جمعیت جواب­ها توسعه داده­اند. ژائو و همکاران[21] (2008) مدل یکپارچه­ای برای مسأله مسیریابی-موجودی در زنجیره تأمین سه رده­ای ارائه کرده و برای حل آن الگوریتم جستجوی همسایگی بزرگ متغیری^[22] توسعه داده­اند. ژائو و همکاران (2007) رویکرد حل جدیدی بر اساس الگوریتم فراابتکاری جستجوی ممنوع^[23] برای حل مسأله مسیریابی-موجودی در یک زنجیره تأمین دو رده­ای توسعه داده­اند. اسپارچی-الکازار و همکاران[24] (2007) با پیاده­سازی الگوریتم ژنتیک برای حل مسأله مسیریابی-موجودی با چند محصول، تأثیر مقادیر مختلف پارامترهای ورودی الگوریتم ژنتیک را ارزیابی کرده­اند تا بهترین مجموع مقادیر پارامترهای الگوریتم بدست آید. راسدیانسیا و سائو[25] (2005) با توسعه مدلی برای مسأله مسیریابی-موجودی ماشین­های سکه­ای، به ارائه الگوریتمی دومرحله­ای مبتنی بر الگوریتم­های الحاق[26] و جستجوی ممنوع پرداخته­اند.

این مقاله به بررسی مسأله مسیریابی-موجودی در حالت چند دوره­ای چند محصولی و با هدف کمینه­سازی مجموع هزینه­های سیستم می­پردازد. هزینه­های سیستم شامل هزینه­های راه­اندازی، توزیع و نگهداری موجودی است. علاوه بر مدیریت موجودی و برنامه­ریزی توزیع، برنامه­ریزی تولید نیز در مسأله، ادغام شده است. در این مقاله، مسأله مسیریابی-موجودی در یک زنجیره تأمین دوبخشی متشکل از یک تولیدکننده و مجموعه­ای از خرده‌فروشان بررسی شده است که در آن، در هر دوره مقدار مشخصی از هریک از محصولات با استفاده از ناوگان همسانی از وسایل حمل با ظرفیت محدود و تحت استراتژی ارسال مستقیم، بین خرده­فروشان توزیع می­گردد. خدمت­رسانی به هریک از خرده­فروشان توسط وسایل حمل، یک دوره زمانی طول می­کشد. بنابراین، هر وسیله حمل در هر دوره تنها می­تواند به یک خرده­فروش محصول تحویل دهد. همچنین فرض شده است که محصول ارسالی به هریک از خرده­فروشان در هر دوره تنها باید طی یک ارسال به خرده­فروش تحویل داده شود. ظرفیت انباشت و تولید محدود است و کمبود مجاز نیست. مسأله مسیریابی-موجودی چند دوره­ای چند محصولی با ارسال مستقیم تحت مفروضات فوق، حالت توسعه یافته­ای از مسأله تعیین اندازه انباشته با ظرفیت محدود[27] در حالت چند محصولی، با ظرفیت تولید و هزینه­های نگهداری ثابت است. در پژوهش‌های پیشین نشان داده شده است که مسأله تعیین اندازه انباشته در شرایط فوق، یک مسأله با پیچیدگی سخت است (بیترن و یاناس[28]، 1981). بنابراین، می­توان نتیجه گرفت که مسأله بررسی شده در این مقاله نیز، با پیچیدگی سخت است. در نتیجه، حل مسأله مذکور با استفاده از روش­های حل دقیق[29] به ویژه برای مسایل با ابعاد بزرگ، در زمان محاسباتی معقول امکانپذیر نیست. در این صورت، می­توان با استفاده از الگوریتم­های فراابتکاری به جواب­های تقریبی مناسب در زمان محاسباتی معقول دست یافت. بنا به دانش نویسندگان، در پژوهش­های پیشین الگوریتم ابتکاری یا فراابتکاری برای حل مسأله مسیریابی-موجودی با مفروضات فوق، توسعه داده نشده است. از این رو، برای حل مسأله، الگوریتم بهینه­سازی گروه ذرات بهبودیافته[30] (IPSO) کارایی توسعه داده شده است. کارایی الگوریتم پیشنهادی با استفاده از مسایل نمونه تصادفی متعددی با الگوریتم­های ژنتیک و بهینه­سازی گروه ذرات[31] (PSO) مقایسه شده است.

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

2- مدل ریاضی

مسأله مسیریابی-موجودی چند دوره­­ای چند محصولی با ارسال مستقیم را می­توان در قالب مدل برنامه­ریزی عدد صحیح آمیخته[32] فرموله کرد. پیش از بیان مدل، ابتدا نمادهای استفاده شده برای توصیف مدل ارائه می­شود:

p : نشانگر محصول

j : نشانگر گره­ها (گره 0 نشانگر تولیدکننده و گره­های 1 تا N نشانگر خرده­فروشان)

t: نشانگر دوره‌های برنامه‌ریزی

N: تعداد خرده­فروشان

: هزینه ثابت راه­اندازی تولید برای تولیدکننده در دوره t.

: میزان موجودی محصول نوع p در پایان دوره t در مرکز j.

: هزینه حمل محصول از تولیدکننده به خرده­فروشی j .

: ضریب مصرف ظرفیت انباشت توسط محصول نوع p از فضای انبار یا وسیله حمل.

: ضریب مصرف ظرفیت تولید توسط محصول نوع p.

: ظرفیت انباشت در مرکز j .

: مقدار تقاضای محصول نوع p در خرده‌فروشی j در دوره t.

: ظرفیت تولیدکننده برای تولید محصول.

: ظرفیت وسیله حمل و نقل.

: هزینه نگهداری یک واحد محصول نوع p در مرکز j.

: میزان محصول نوع p ارسال شده از تولیدکننده به مرکز توزیع j در طی دوره t.

: میزان تولید محصول نوع p در طی دوره t توسط تولیدکننده.

: متغیر صفر و یک نشانگر ارسال محصول به خرده­فروش j در دوره t.

: متغیر صفر و یک نشانگر تولید محصول در دوره t.

مدل ریاضی مسأله، با استفاده از معادله‌های (1) الی (8) ارائه شده است:

(1) 

(2)  

(3)           

(4)                          

(5)               

(6)                        

(7)                               

(8)                        

معادله (1) بیانگر تابع هدف مدل پیشنهادی است که شامل هزینه­های ثابت راه­اندازی، هزینه ارسال محصول به خرده­فروشان، هزینه­های نگهداری محصول توسط تولیدکننده و خرده­فروشان است. معادله (2) بیانگر توزان موجودی بین تقاضا برای محصول نوع p در خرده‌فروشی j در دوره t و مجموع محصول نوع p ارسال شده از تولیدکننده به خرده‌فروشی l در دوره t است. توازن میان مجموع محصول نوع p ارسال شده از تولیدکننده به تمامی خرده‌فروشی‌ها در دوره t و مجموع محصول نوع p که در دوره t ام تولید شده است، با استفاده از معادله (3) بیان می­شود. معادله (4) محدودیت ظرفیت تولید را تضمین می­کند. محدودیت ظرفیت وسایل حمل و نقل با استفاده از معادله (5) ارائه شده است. معادله‌های (4) و (5) بطور ضمنی نشان می‌دهند که آیا در دوره t به ترتیب، توسط تولیدکننده تولیدی صورت گرفته است یا محصولی برای خرده­فروشی j ام ارسال شده است. معادله (6) تضمین­کننده محدودیت ظرفیت نگهداری برای تولیدکننده و خرده‌فروشی­ها است. محدودیت‌های فنی روی متغیرهای تصمیم از طریق معادله‌های (7) و (8) برقرار شده است.

3- الگوریتم بهینه­سازی گروه ذرات پیشنهادی

الگوریتم فراابتکاری بهینه­سازی گروه ذرات روش محاسباتی تکاملی مبتنی بر جمعیت جواب­ها است. مانند سایر الگوریتم­های فراابتکاری، الگوریتم مذکور ابزار بهینه­سازیی است که می­تواند برای حل انواع مختلفی از مسایل بهینه­سازی به­کار گرفته شود. این الگوریتم از جدیدترین روش­های فراابتکاری است که با الهام­گیری از رفتار اجتماعی گروهی از پرندگان مهاجر که در تلاش برای دستیابی به مقصد ناشناخته­ای هستند، توسط کندی و ابرهارت[33] (1995) در سال 1995 میلادی توسعه داده شده است. در الگوریتم PSO، جمعیت جواب­ها، گروه^[34] نامیده می­شود و هر جواب مانند یک پرنده در گروهی از پرندگان است و ذره^[35] نام دارد و شبیه کرموزوم در الگوریتم ژنتیک است. تمامی ذرات دارای مقدار شایستگی^[36] هستند که با استفاده از تابع شایستگی^[37] محاسبه می­گردند و تابع شایستگی ذرات باید بهینه گردد. جهت حرکت هر ذره توسط بردار سرعت^[38] آن ذره معین می­شود. برخلاف الگوریتم ژنتیک، در فرآیند تکاملی الگوریتم مذکور، پرندگان جدیدی از نسل قبل (جواب­های جدید از جواب‌های قبلی) ایجاد نمی­گردد، بلکه هر پرنده رفتار اجتماعی خود را با توجه به تجربیاتش و رفتار سایر پرندگان گروه تکامل بخشیده و مطابق آن حرکت خود را به سوی مقصد بهبود می­دهد. به عبارت دیگر، در این الگوریتم عملگرهای تکاملی چون تقاطع^[39] و جهش^[40] وجود ندارد (هو و همکاران[41]، 2004؛ جائو و همکاران[42]، 2006). ساختار کلی الگوریتم بهینه­سازی گروه ذرات در شکل 1 ارائه شده است.

شکل 1. ساختار الگوریتم بهینه­سازی گروه ذرات

الگوریتم PSO با گروهی از جواب­های (ذرات) تصادفی آغاز و سپس با به هنگام­سازی ذرات در هر تکرار به دنبال جواب بهینه می­گردد (هو و همکاران، 2004). اگر متغیرهای تصمیم و به نوبه آن موقعیت ذرات، از نوع صفر و یک باشند؛ بردارهای سرعت و موقعیت هریک از ذرات در هر تکرار الگوریتم، طبق روابط (9) الی (12) محاسبه می­شوند (انگلبرکت[43]، 2005):

(9)(10)                              

(11)                                   

(12)                          

طبق رابطه (9)بردار سرعت جدید هر ذره بر اساس سرعت قبلی خود ذره (V_i(t-1))، بهترین موقعیتی که ذره تاکنون به آن دست یافته است (pBest_i) و موقعیت بهترین ذره در همسایگی ذره که تابحال بدست آمده است (nBest_i)، محاسبه می‌گردد. در صورتی که همسایگی هر ذره شامل تمام ذرات گروه باشد، آنگاه nBest_i بیانگر موقعیت بهترین ذره در میان گروه است که با gBest به آن اشاره می­شود. r₁ و r₂ دو عدد تصادفی (با توزیع یکنواخت بین [0,1]) هستند که مستقل از یکدیگر تولید می­شوند. c₁ و c₂ که با نام ضرایب یادگیری به آنان اشاره شده است، تأثیر pBest و nBest را بر فرآیند جستجو کنترل می­نمایند. w بیانگر ضریب وزنی اینرسی است. بردار سرعت ذرات با مقدار V_max محدود شده است. V_max به عنوان محدودیتی است که قابلیت جستجوی جهانی گروه ذرات را کنترل می­کند. با استفاده از رابطه (11) بردار سرعت هریک از ذرات به بردار احتمال تغییر، تبدیل می­شود. در رابطه فوق، s_i بیانگر احتمال آن است که x_it  برابر با 1 شود. سپس با استفاده از رابطه (12) بردار موقعیت هریک از ذرات بهنگام می­گردد. در رابطه فوق،عددی تصادفی با توزیع یکنواخت بین صفر و یک است.

مسأله مسیریابی-موجودی چند دوره­ای چند محصولی با ارسال مستقیم تحت مفروضات مورد بررسی در این مقاله، از جمله مسایل با درجه پیچیدگی سخت است. در نتیجه در این مقاله، الگوریتم بهینه­سازی گروه ذرات بهبودیافته جدیدی برای حل مسأله فوق توسعه داده شده است.

با توجه به معادلات (1) الی (8) مسأله مسیریابی-موجودی، یک مسأله برنامه­ریزی عدد صحیح آمیخته است که در آن متغیرهای تصمیم به دو دسته متغیرهای پیوسته (w_pjt، I_pjt و P_pt) و متغیرهای صفر و یک (Z_t و X_jt) قابل تفکیک هستند. از آنجایی که دشواری حل مسأله مذکور ناشی از متغیرهای صفر و یک است؛ بنابراین، با تفکیک روی متغیرهای تصمیم، می­توان راه­حلی دو بخشی برای حل مسأله توسعه داد. در راه­حل پیشنهادی ابتدا با استفاده از الگوریتم بهینه­سازی گروه ذرات بهبودیافته، مقادیر متغیرهای صفر و یک (Z_t و X_jt) تعیین می­شود، سپس با توجه به مشخص بودن مقادیر متغیرهای صفر و یک، مسأله مسیریابی-موجودی، به یک مسأله برنامه­ریزی خطی برای تعیین مقادیر متغیرهای پیوسته (w_pjt، I_pjt و P_pt) با هدف کمینه­سازی مجموع هزینه­ها تبدیل می­شود که به سادگی می­توان مسأله برنامه­ریزی خطی حاصل را با روش­های موجود برای حل مسایل برنامه­ریزی خطی حل کرد و سرانجام، به جواب نهایی مسأله اصلی دست یافت.

ساختار کلی الگوریتم پیشنهادی در شکل 2 ارائه شده است. همانگونه که در شکل نشان داده شده است، ساختار پیشنهادی مبتنی بر الگوریتم بهینه­سازی گروه ذرات، طراحی شده است. مطابق با شکل 2، الگوریتم پیشنهادی از دو مرحله تشکیل شده است که مرحله اول خود متشکل از دو بخش مجزاست. در بخش نخست، ابتدا به تعداد SSize ذره با استفاده از روشی شبه تصادفی به عنوان ذرات اولیه، ایجاد می­شوند. پس از آن، همسایگی تصادفی به اندازه NSize ذره، حول هر جواب اولیه، شکل می­گیرد. در طی گام­های الگوریتم، کیفیت و پراکندگی جواب­ها بهبود می­یابد و حداکثر به تعداد SSize ذره از جواب­های خوب و غیرتکراری انتخاب و در مجموعه­ای با نام مجموعه مرجع -برای استفاده در مراحل بعدی الگوریتم- قرار داده می­شود. در پایان بخش نخست، جواب­های موجود در مجموعه مرجع به عنوان ذرات اولیه برای شروع بخش دوم در نظر گرفته می­شوند. در صورتی که تعداد ذرات موجود در مجموعه مرجع کمتر از SSize ذره باشد، ذرات باقیمانده موردنیاز به صورت تصادفی ایجاد و افزوده می­شوند. در بخش دوم نیز فرآیند بهبود جواب­ها با مکانیزمی متفاوت ادامه می­یابد و در هر تکرار الگوریتم جواب­های خوب و غیرتکراری در مجموعه مرجع قرار داده می­شوند. در پایان بخش دوم، جواب­های موجود در مجموعه مرجع به عنوان ذرات اولیه برای شروع بخش نخست در نظر گرفته می­شوند. در این بخش نیز، در صورتی که تعداد ذرات موجود در مجموعه مرجع کمتر از SSize ذره باشد، ذرات باقیمانده موردنیاز به صورت تصادفی تولید و اضافه می­شوند. بخش اول و دوم، هر کدام به ترتیب پس از P1S1 و P1S2 مرتبه تکرار یا در صورتی که به ترتیب پس از طی CR11 و CR12 مرتبه تکرار در مقدار شایستگی بهترین جواب تغییری حاصل نشود، خاتمه می­یابند. همچنین، مرحله اول الگوریتم -رویه تکرار پیاپی بخش اول و دوم- پس از طی P1 مرتبه تکرار پایان می­پذیرد. در پایان مرحله نخست الگوریتم، مجموعه مرجع تشکیل شده در طی این مرحله، به عنوان جواب­های آغازین مرحله دوم الگوریتم، در نظر گرفته و الگوریتم وارد مرحله دوم می­شود. در مرحله دوم، ابتدا با استفاده از مدل برنامه­ریزی خطی، مقادیر بهینه متغیرهای تصمیم پیوسته با توجه به مقادیر متغیرهای صفر و یک (Z_t و X_jt) تعیین و از بین جواب­ها، بهترین جواب انتخاب می­شود. پس از آن با استفاده از رویه جستجوی محلی کیفیت جواب منتخب بهبود می­یابد. مرحله دوم الگوریتم، بعد از P2 مرتبه تکرار یا در صورتی که پس از CR2 مرتبه تکرار تغییری در مقدار بهترین ذره حاصل نگردد، خاتمه می­یابد.

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

در مرحله اول الگوریتم، بردارهای سرعت و موقعیت هریک از ذرات در هر تکرار الگوریتم، طبق روابط (9) الی (12) محاسبه می­شوند. در بخش دوم از مرحله نخست الگوریتم، در رابطه (9) مقدار nBest_i با مقدار بهترین ذره در میان تمام ذرات (gBest) جایگزین می­شود.

شکل 2. ساختار الگوریتم بهینه­سازی گروه ذرات پیشنهادی

3-1- نحوه نمایش ذرات

یکی از مهمترین تصمیماتی که در زمان طراحی یک الگوریتم فراابتکاری می­بایستی اتخاذ شود، نحوه نمایش جواب­ها[44] و نیز چگونگی برقراری ارتباطی مؤثر، یکتا و قابل شناسایی میان جواب­ها و فضای جستجوی مسأله مورد بررسی است. در این مقاله برای نمایش ذرات، با توسعه روش عزیز و معین (2007)، از شیوه "نمایش زمان ارسال" برای نمایش بردار موقعیت ذرات استفاده شده است. در این شیوه نمایش از مقادیر صفر و یک برای نمایش زمان ارسال محصولات به خرده­فروشان مختلف و زمان تولید استفاده شده است (مقادیر متغیرهای Z_t و X_jt). هر 1 نشاندهنده تولید (1=Z_t) یا ارسال محصول به خرده­فروش در آن دوره (1=X_jt) و مقدار صفر بیانگر عدم تولید (0=Z_t) یا ارسال محصول در آن دوره (0=X_jt) است. نمایی کلی از بردار موقعیت ذرات در نحوه نمایش زمان ارسال برای حل مسأله مسیریابی-موجودی با ارسال مستقیم برای حالتی با 6 دوره برنامه­ریزی و 3 خرده­فروش، در شکل 3 ارائه شده است. موجودی خرده­فروشان و نیز تولیدکننده در ابتدای دوره برنامه­ریزی صفر فرض شده است، بنابراین، در نخستین دوره برنامه­ریزی تمام مقادیر مربوط به محصولات و خرده­فروشان برابر با 1 است.

شکل 3. نمونه­ای از نحوه نمایش ذرات

3-2- تولید جواب­های اولیه

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

1-    تمام مؤلفه­های ذره نخست را برابر 1 قرار دهید. این ذره، ذره­ای با بیشترین تعداد دفعات ارسال و همواره شدنی است. به عبارت دیگر، ذره نخست بیانگر حالتی است که در آن، در تمام دوره­ها تولید و در تمام دوره­ها و برای تمام خرده­فروشان ارسال داشته باشیم.

2-    با استفاده از رویه ارائه شده در شکل 4، ذره دوم را ایجاد کنید. ذره دوم، ذره­ای با کمترین تعداد دفعات ارسال است. برای تولید ذره دوم، محدودیت انباشت که با استفاده از رابطه (6) بیان می­شود، آزاد شده است.

شکل 4. شبه کد تولید ذره با کمترین تعداد دفعات تولید و ارسال

3-    ذره اول و دوم را به عنوان والد اول و دوم قرار دهید. پس از آن، با استفاده از عملگر تقاطع پراکنده[45] سایر ذرات موردنیاز را ایجاد کنید. برای انجام عملگر تقاطع پراکنده، ابتدا یک آرایه تصادفی صفر و یک با اندازه برابر با اندازه ذرات ایجاد می­گردد. مؤلفه i ام ذره فرزند (ذره جدید) اگر مؤلفه i ام آرایه تصادفی 1 بود، برابر با مؤلفه i ام والد 1 (ذره اول) و در غیر اینصورت برابر با مؤلفه i ام والد 2 (ذره دوم) خواهد بود. نحوه انجام عملگر تقاطع پراکنده در شکل 5 ارائه شده است.

شکل 5. نحوه انجام عملگر تقاطع پراکنده

3-3. محاسبه مقادیر شایستگی

به طور معمول، در الگوریتم­های فراابتکاری هر جواب در هر مرحله با استفاده از تابع شایستگی ارزیابی و مقدار شایستگی آن محاسبه می­شود. در این مقاله نیز مانند سایر مسایل تک هدفه از تابع هدف مسأله که توسط رابطه (1) بیان شده است، به عنوان تابع شایستگی ذرات استفاده می­شود. در ساختار الگوریتم پیشنهادی، موقعیت هر ذره، تنها بیانگر متغیرهای Z_t و X_jtاست. بنابراین، برای محاسبه تابع هدف لازم است، w_pjt، I_pjt و P_pt و در پی آن مقدار تابع هدف بر اساس Z_t و X_jtمحاسبه شود. به ازای مقادیر ثابتی از متغیرهای Z_t و X_jt، می‌توان مقادیر مختلفی برای متغیرهای w_pjt، I_pjt و P_pt بدست آورد که از میان این مقادیر، w_pjt، I_pjt و P_pt را باید به گونه­ای تعیین کرد که ضمن برقراری محدودیت­های مسأله، کمترین هزینه­ها را نیز دربرداشته باشد. در الگوریتم پیشنهادی برای یافتن مقادیر بهینه مذکور، از یک مدل برنامه­ریزی خطی استفاده می­شود. بدین منظور، با توجه به روابط (2)و (3) و با توجه به اینکهخواهیم داشت:

(13)  

(14)          

با جایگذاریاز روابط (13)و (14) در رابطه (1) و با مشخص بودن زمان­های تولید و ارسال، پس از ساده­سازی خواهیم داشت:

(15)

همچنین با جایگذاریاز روابط (13)و (14) در رابطه (6) خواهیم داشت:

(16)      

(17)  

همچنین با توجه به روابط (13)، (14) و چون داریم ، می­توان محدودیت­های (2) و (3) را با محدودیت­های (18) و (19) جایگزین کرد:

(18)           

(19)                 

بنابراین از روابط فوق و با توجه به معادلات (1)الی(8)، مدل ساده شده به صورت زیر است:

(20)

(21)           

(22)         

(23)                                  

(24)    

(25)  

(26)                 

(27)                                  

و معادلات (16) و (19).

که در رابطه (22)،
  است. مدل ساده شده فوق به تعداد P(N+1)T متغیر و حداکثرN(P+1) محدودیت کمتر از مدل اصلی دارد. حال در هر مرحله از الگوریتم، با مشخص بودن زمان­های تولید و ارسال (و) و با استفاده از مدل مذکور، می­توان متغیرهای w_pjt و P_pt و نیز مقادیر شایستگی هریک از ذرات را محاسبه و با توجه به مقادیر حاصل از حل مدل برنامه­ریزی خطی برای متغیرهای و، در صورت نیاز، اصلاح لازم در زمان­های تولید و ارسال را نیز اعمال کرد -زیرا ممکن است با حل مدل برنامه­ریزی خطی مقدار بهینه برخی از مقادیر تولید یا ارسال صفر بدست آید که در این صورت زمان ارسال متناظر با آن نیز صفر خواهد شد.

محاسبه مقدار شایستگی ذرات با استفاده از حل مدل برنامه­ریزی خطی تا حدی به زمان محاسباتی بالایی نیاز دارد، بنابراین برای افزایش سرعت الگوریتم، فقط در مرحله دوم الگوریتم از این روش استفاده می­شود و برای محاسبه مقدار شایستگی در مرحله نخست الگوریتم، روشی تقریبی که در ادامه تشریح خواهد شد، توسعه داده شده است:

1-    با استفاده از رابطه (28) مقدار اولیه w_pjt را محاسبه کنید:

(28)               

که در رابطه مذکور، r نشانگر دوره بعدی است که ارسال بعدی در آن صورت می­گیرد.

2-    با استفاده از روابط (29) تا (31) مقدار بدست آمده برای w_pjt را تعدیل کنید به گونه­ای که تمام مقادیر w_pjt شدنی شود ():

(29)   

(30)                

(31)

3-    مقدار اولیه P_pt را با استفاده از رابطه (32) محاسبه کنید:

(32)                

که در رابطه مذکور، r نشانگر دوره بعدی است که تولید بعدی در آن صورت می­گیرد.

4-    با استفاده از روابط (33) تا (35) مقدار بدست آمده برای P_pt را تعدیل کنید به گونه­ای که تمام مقادیر P_pt شدنی شود ():

(33)               

(34)   

(35)

5-    با توجه به روابط (36) و (37) مقدار I_p0t و I_pjt را محاسبه کنید:

(36)        

(37)           

6-    با به­کارگیری روابط (38) تا (40) میزان نشدنی بودن جواب­های حاصل را بدست آورید:

(38)                          

(39)                    

(40)

که در روابط بالا،  ضریب جریمه برای تخطی از محدودیت k ام است.

7-    با استفاده از رابطه (41) مقدار شایستگی ذره را محاسبه کنید:

(41)                

که در رابطه فوق، Z(X) بیانگر تابع هدف حاصل از رابطه (1) است.

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

تفاوت سرعت الگوریتم پیشنهادی در دو حالت محاسبه مقدار شایستگی ذرات با استفاده از حل مدل برنامه­ریزی خطی در هر دو مرحله (LP-IPSO) و محاسبه مقدار شایستگی ذرات با استفاده از رویه تقریبی فوق در مرحله نخست و استفاده از مدل برنامه­ریزی خطی فقط در مرحله دوم الگوریتم (IPSO) را می­توان با استفاده از چند مثال بررسی کرد. میانگین نتایج حاصل از 10 بار اجرای الگوریتم پیشنهادی در دو حالت مذکور برای حل چند مسأله نمونه (مسایل نمونه در بخش 4 تشریح خواهند شد)، در جدول 1 ارائه شده است.

جدول 1. نتایج حاصل از حل مسایل نمونه با استفاده از LP-IPSO و IPSO

شکل 6. میزان اختلاف زمان محاسباتی موردنیاز برای حل مسایل نمونه با استفاده از LP-IPSO و IPSO

همانگونه که در جدول 1 و شکل 6 ارائه شده است، هرچند LP-IPSO به مقادیر تابع هدف بهتری نسبت به IPSO دست یافته است، اما به زمان محاسباتی به مراتب بیشتری نیاز دارد. با افزایش ابعاد مسأله، تفاوت سرعت الگوریتم در دو حالت افزایش می­یابد.

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

(42)

در رابطه بالا،ضریب جریمه تخطی از محدودیت k ام در تکرار t ام،میزان افزایش در ضریب جریمه و infeasibility(X) میزان نشدنی بودن بهترین ذره است. ضریب جریمه در طی گام­های مرحله دوم الگوریتم ثابت باقی می­ماند.

3-4. ایجاد همسایگی تصادفی ذرات

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

1-    تعداد NSize-1 کپی از ذره ایجاد می­شود.

2-    به تعداد NSize-1 آرایه تصادفی صفر و یک، با تعداد مؤلفه­های برابر با مؤلفه­های یک ذره، ایجاد می­گردد. HDRate درصد از مؤلفه­های آرایه تولیدشده، به صورت تصادفی برابر 1 قرار داده می­شود.

3-    مطابق با مؤلفه­های با مقدار 1 در آرایه تصادفی، مؤلفه متناظر در ذره کپی معکوس می­گردد؛ به عبارت دیگر، مقدار 0 به 1 و بالعکس مقدار 1 به 0 تبدیل می­شود.

نمایی کلی از نحوه انجام عملگر جهش در شکل 7 ارائه شده است.

شکل 7. نحوه انجام عملگر جهش

3-5. بهبود همسایگی ذرات

پس از ایجاد همسایگی ذرات، طی گام­های الگوریتم، جواب­های موجود در هر همسایگی مستقل از ذرات سایر همسایگی­ها بهبود می­یابند. اما به منظور افزایش میزان بهبود در کیفیت جواب­ها طی گام­های الگوریتم، لازم است میان همسایگی­های مختلف ذرات، به گونه­ای تبادل اطلاعات صورت گیرد. بدین منظور، پس از بهنگام­سازی مقادیر nBest در هر تکرار الگوریتم، با به­کارگیری روش انتخاب مسابقه[46] به تعداد 2SSize ذره از میان ذرات nBest به عنوان والدین[47] انتخاب می­شوند. پس از آن با بهره­گیری از عملگر تقاطع پراکنده که در بخش 3-2 تشریح شد، SSize ذره جدید (فرزند[48]) ایجاد می­شود. ذرات جدید، جایگزین بدترین ذره در هر همسایگی می­شوند.

3-6. تشکیل و بهنگام­سازی مجموعه مرجع (RSet)

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

3-7. متنوع­سازی ذرات

در ساختار الگوریتم پیشنهادی، به منظور جلوگیری از همگرایی زودهنگام و ایجاد جواب­های اولیه با پراکندگی و کیفیت خوب برای مراحل بعدی الگوریتم، مکانیزم متنوع­سازی ذرات انجام می­گیرد. به همین منظور، در بخش نخست مرحله اول الگوریتم، بهترین ذره در میان مجموعه مرجع انتخاب و به تعداد MRate درصد از ذرات گروه، از ذره منتخب کپی ایجاد می­شود. ذرات کپی با به­کارگیری عملگر جهش که در بخش 3-4 تشریح شد، تغییر می­یابند و جایگزین بدترین ذرات گروه می­شوند. با استفاده از این روش، با اعمال تغییر تصادفی محدود در بهترین ذره گروه، با افزایش مؤلفه­های خوب در میان جواب­ها، فرآیند جستجو بهبود می­یابد.

3-8. جستجوی محلی

در مرحله دوم از ساختار الگوریتم پیشنهادی، پس از انتخاب بهترین ذره، کیفیت ذره منتخب با استفاده از جستجوی محلی بهبود می­یابد. برای این هدف، یکی از مؤلفه­های ذره به صورت تصادفی انتخاب و مقدار آن مؤلفه معکوس می­شود؛ به عبارت دیگر، اگر مؤلفه برابر 1 باشد به 0 و بالعکس اگر 0 باشد به 1 تبدیل می­شود. بدین ترتیب ذره جدیدی تنها با یک مؤلفه متفاوت نسبت به ذره اولیه بدست می­آید. با استفاده از مدل برنامه­ریزی خطی که در بخش 3-3 ارائه شد، مقدار شایستگی ذره جدید محاسبه می­گردد. اگر مقدار شایستگی ذره جدید از ذره اولیه بهتر باشد، ذره جدید جایگزین ذره اولیه می­شود و در غیر اینصورت ذره جدید حذف و ذره اولیه بدون تغییر باقی می­ماند. رویه مذکور حداکثر پس از P2 مرتبه تکرار یا در صورتی که در بهترین مقدار شایستگی که ذره بدان دست یافته است، طی CR2 تکرار متوالی بهبودی حاصل نشود، الگوریتم متوقف می­گردد.

4- نتایج محاسباتی

با استفاده از مسایل نمونه تصادفی متعددی، عملکرد الگوریتم پیشنهادی در مقایسه با الگوریتم بهینه­سازی گروه ذرات، الگوریتم ژنتیک و نرم­افزار Lingo (تنها برای مسایل با ابعاد کوچک) ارزیابی می­شود. با مقایسه عملکرد الگوریتم پیشنهادی با الگوریتم بهینه­سازی گروه ذرات، می­توان به میزان بهبود حاصل از ساختار جدید طراحی شده برای الگوریتم پیشنهادی نسبت به ساختار پایه­ای الگوریتم بهینه­سازی گروه ذرات پی برد. الگوریتم ژنتیک نیز، از جمله معروفترین و پرکاربردترین الگوریتم­های فراابتکاری است که توسط پ‍ژوهش­های متعددی در حل مسایل بهینه­سازی مختلفی به­کار گرفته شده است. به همین دلیل، برای ارزیابی عملکرد الگوریتم پیشنهادی، از دو الگوریتم فوق به عنوان الگوریتم­های معیار استفاده شده است.

الگوریتم­های پیشنهادی و بهینه­سازی گروه ذرات، در محیط Matlab 7.5 برنامه­نویسی و الگوریتم ژنتیک با استفاده از جعبه ابزار نرم­افزار Matlab 7.5 پیاده­سازی شده است. تمامی الگوریتم­های مذکور در کامپیوتری با مشخصات پردازنده AMD، 64 بیتی، 3 گیگاهرتز، ویندوز XP و 2 گیگابایت رم، اجرا شده­اند.

4-1- چگونگی ایجاد مسایل نمونه

مسایل نمونه در دو گروه با ابعاد کوچک و بزرگ به صورت تصادفی ایجاد شده­اند. ابعاد مسایل، تعداد دوره­های برنامه­ریزی و تعداد محصولات در جدول 2 ارائه شده است. همانگونه که در جدول مذکور نشان داده شده است، در تولید مسایل نمونه، هزینه حمل محصولات از مسایل معیار مسأله مسیریابی وسیله نقلیه با ظرفیت محدود[49] استخراج شده است (مجموعه داده­های مسیریابی وسیله نقلیه[50]). اطلاعات لازم برای تولید مسایل نمونه در جدول 3 ارائه شده است. در این مقاله، همانند بودیا و پرینز (2009) ظرفیت تولید و وسایل حمل متناسب با تقاضا، هزینه ثابت راه­اندازی تولید و همچنین ظرفیت انبار تولیدکننده متناسب با ظرفیت تولید در نظر گرفته شده است. ظرفیت انبار خرده‌فروشان نیز با ظرفیت وسایل حمل متناسب است. بدین منظور، ظرفیت وسایل حمل 2 برابر () میانگین تقاضای یک خرده­فروش و ظرفیت تولید 5/3 برابر میانگین تقاضای روزانه تمام خرده­فروشان () تعیین شده است. هزینه ثابت راه­اندازی تولید 5/1 برابر () و ظرفیت انبار تولیدکننده برابر با () ظرفیت تولید در نظر گرفته شده است. ظرفیت انبار تمام خرده­فروشان یکسان و برابر با ظرفیت وسیله حمل () تعیین شده است. در جدول 3، d_p بیانگر میانگین تقاضای محصول p ام در یک دوره است. در تولید مسایل نمونه با 5 محصول، 5 سطح متفاوت تقاضا شامل تقاضای خیلی کم (d_1jt)، کم (d_2jt)، متوسط (d_3jt)، زیاد (d_4jt) و خیلی زیاد (d_5jt) مدنظر قرار گرفته است که مقدار تقاضای هر محصول در هر دوره به صورت یکنواخت در بازه ارائه شده در جدول 3، ایجاد می­شود. برای مسایل با 3 محصول، 3 سطح تقاضای خیلی کم (d_1jt)، متوسط (d_3jt) و خیلی زیاد (d_5jt) لحاظ شده است.

جدول 2. چگونگی ایجاد مسایل نمونه

جدول 3. نحوه تولید پارامترهای مسایل نمونه

4-2- مفروضات و پارامترهای الگوریتم­ها

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

الگوریتم ژنتیک: تعداد تکرارها، اندازه جمعیت، تعداد جواب نخبه[51] و ضریب جریمه برای مسایل نمونه کوچک به ترتیب برابر با 200، 100، 10 و 100 و برای مسایل نمونه بزرگ به ترتیب برابر با 1500، 350، 50 و 5^10 قرار داده شده است. سایر پارامترها همانند گزینه­های پیش فرض جعبه­ابزار Matlab در نظر گرفته شده است.

الگوریتم بهینه­سازی گروه ذرات: تعداد تکرارها، تعداد ذرات، ضرایب یادگیری (c₁ و c₂)، V_max و ضریب جریمه برای مسایل نمونه کوچک به ترتیب برابر با 200، 100، 2، 2، 6 و 100 و برای مسایل نمونه بزرگ به ترتیب برابر با 1500، 350، 2، 2، 6 و 5^10 قرار داده شده است. ضریب اینرسی به صورت خطی از مقدار 1.4 به 0.9 در طول اجرای الگوریتم کاهش می­یابد.

الگوریتم بهینه­سازی گروه ذرات بهبودیافته: تعداد تکرارهای بخش اول (P1S2) و دوم (P1S1)، مرحله اول (P1) و دوم (P2)، حداکثر تکرار بدون بهبود در بهترین جواب برای بخش اول (C11) و دوم (C12) و مرحله دوم (C2)، تعداد ذرات (SSize)، اندازه همسایگی ذرات (NSize)، ضرایب یادگیری (c₁ و c₂)، V_max در بخش اول و دوم، HDRate، MRate، ضریب جریمه اولیه برای بخش اول () و دوم () و مرحله دوم () و میزان افزایش در ضریب جریمه در هر گام () برای مسایل نمونه کوچک به ترتیب برابر است با 50، 50، 10، 250، 25، 25، 50، 20، 10، 2، 2، 3، 6، 0.07، 0.1، 10، 10، 75 و 0.1 و برای مسایل نمونه بزرگ به ترتیب برابر با 75، 75، 10، 350، 35، 35، 75، 30، 10، 2، 2، 6، 6، 0.1، 0.1، 100، 100، 750 و 0.5 قرار داده شده است. ضریب اینرسی به صورت خطی از مقدار 1.4 به 0.9 در طول اجرای الگوریتم کاهش می­یابد.

4-3- نتایج عددی

هریک از الگوریتم­ها برای حل هر مسأله نمونه، 10 بار اجرا شده­اند که با تحلیل نتایج حاصل، از دیدگاه 3 معیار عملکردی[52] به ارزیابی عملکرد الگوریتم پیشنهادی پرداخته شده است: 1) میانگین مقادیر تابع هدف بدست آمده، 2) میانگین زمان محاسباتی و 3) میزان استواری[53] نتایج.

هر اندازه تغییرات مقادیر تابع هدف در اجراهای مختلف الگوریتمی مشخص، کمتر باشد؛ میزان استواری الگوریتم بیشتر خواهد بود. معیار عملکردی میزان استواری نتایج با استفاده از رابطه (43) محاسبه می­گردد (انگلبرکت، 2005):

(43)             

در رابطه فوق، میانگین و واریانس مقادیر تابع هدف بدست آمده در اجراهای مختلف الگوریتم است. هر اندازه واریانس نتایج و به نوبه آن بازه ارائه شده در رابطه (43) کوچکتر باشد، میزان استواری الگوریتم بیشتر است.

در جداول 4 و 5، میانگین مقادیر تابع هدف بدست آمده و زمان­های محاسباتی صرف شده توسط سه الگوریتم برای حل هر مسأله، ارائه شده است. در برخی از موارد، در حل یک مسأله مشخص، پس از اختتام الگوریتم، جوابی نشدنی حاصل شده است. در این صورت، میزان تخطی از محدودیت­ها در هر اجرای الگوریتم برای آن مسأله تعیین و ذخیره و در پایان میانگین میزان تخطی از محدودیت­ها در تمام اجراهای الگوریتم برای آن مسأله محاسبه می­شود. در جداول 4 و 5، میانگین مقدار تخطی از محدودیت­ها در حل یک مسأله، در ستون نشدنی بودن ارائه شده است. با توجه به جداول مذکور، الگوریتم بهینه­سازی گروه ذرات بهبودیافته پیشنهادی در تمام مسایل نمونه نسبت به الگوریتم­های معیار، به نتایج بهتری دست یافته است. در مسایل نمونه با ابعاد کوچک، الگوریتم پیشنهادی توانسته است به مقادیری بسیار نزدیک به بهینه و در مسأله نمونه 1 به مقدار بهینه دست یابد. همچنین، الگوریتم پیشنهادی در تمام مسایل نمونه به جوابی شدنی رسیده است، در حالی که الگوریتم بهینه­سازی گروه ذرات در برخی از موارد به نتایج نشدنی منتج شده است.

جدول 4. میانگین مقدار تابع هدف و زمان محاسباتی حل مسایل نمونه تصادفی با ابعاد کوچک

- بهترین مقدار بدست آمده برای هر مسأله نمونه پررنگ شده است.

*  نرم­افزار به جواب بهینه نرسیده و متوقف شده است، بهترین جواب بدست آمده گزارش شده است.

جدول 5. میانگین مقدار تابع هدف و زمان محاسباتی حل مسایل نمونه تصادفی با ابعاد بزرگ

- بهترین مقدار بدست آمده برای هر مسأله نمونه پررنگ شده است.

برای آزمودن وجود تفاوتی معنی­دار میان عملکرد الگوریتم پیشنهادی با الگوریتم­های معیار از تحلیل واریانس دو طرفه[54] بر روی مقادیر تابع هدف بدست آمده توسط سه الگوریتم در حل مسایل با ابعاد مختلف استفاده می­کنیم. بدین منظور الگوریتم­ها به عنوان تیمار[55] و مسایل با ابعاد مختلف به عنوان بلوک[56] در نظر گرفته شده­اند. تحلیل واریانس با استفاده از نرم­افزار Minitab 14 انجام شده است که نتایج تحلیل واریانس مسایل با ابعاد کوچک و بزرگ، به ترتیب در شکل­های 8 و 9 ارائه شده است. با توجه به شکل­های 8 و 9، در هر دو گروه مسایل نمونه با ابعاد کوچک و بزرگ در سطح معنی­دار 01/0 () می­توان یکسان نبودن عملکرد الگوریتم­ها را نتیجه گرفت.

شکل 8. تحلیل واریانس دو طرفه بر روی نتایج حاصل از حل مسایل با ابعاد کوچک

شکل 9. تحلیل واریانس دو طرفه بر روی نتایج حاصل از حل مسایل با ابعاد بزرگ

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

شکل 10. مقایسه زمان محاسباتی الگوریتم پیشنهادی و الگوریتم­های معیار در حل مسایل با ابعاد کوچک

شکل 11. مقایسه زمان محاسباتی الگوریتم پیشنهادی و الگوریتم­های معیار در حل مسایل با ابعاد بزرگ

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

میزان استواری هر سه الگوریتم در حل مسایل کوچک و بزرگ، به ترتیب در جداول 6 و 7 ارائه شده است. مطابق با جداول مذکور، به طور میانگین الگوریتم پیشنهادی در مقایسه با دو الگوریتم معیار در هر دو گروه مسایل کوچک و بزرگ از تغییرپذیری کمتر و به نوبه آن از استواری بیشتری برخوردار است.

عملکرد بهتر الگوریتم پیشنهادی را می­توان ناشی از ساختار جدید طراحی شده برای آن دانست. این مطلب از چهار بعد قابل بررسی است: 1) نحوه تولید جواب­های اولیه در الگوریتم پیشنهادی به گونه­ای است که جواب­هایی خوب برای مراحل بعدی فراهم می­گردد، 2) با تشکیل مجموعه مرجع در مرحله نخست الگوریتم، جواب­های با کیفیت و منحصربفردی که در هر تکرار تولید شده­اند، حفظ می­شوند. بدین ترتیب جواب­هایی خوب و متنوع برای گام­های بعدی بدست می­آید و جواب­های نامناسب از چرخه الگوریتم حذف می­شوند، 3) با تشکیل همسایگی حول ذرات در بخش نخست الگوریتم، همسایگی هر ذره با دقت بیشتری جستجو و امکان بهبود در ذرات بررسی می­شود، 4) استفاده از رویه تقریبی در مرحله نخست و مدل برنامه­ریزی خطی در مرحله دوم الگوریتم پیشنهادی برای محاسبه تابع هدف ذرات، به شدت بر سرعت الگوریتم می­افزاید. بنابراین، می­توان نتیجه گرفت که الگوریتم پیشنهادی با سرعت و ثبات بیشتری به نتایج بهتری دست می­یابد.

4-4- تحلیل حساسیت

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

جدول 6. میزان استواری الگوریتم­ها در حل مسایل نمونه تصادفی با ابعاد کوچک

جدول 7. میزان استواری الگوریتم­ها در حل مسایل نمونه تصادفی با ابعاد بزرگ

حمل به تقاضا، ارسال و تولید محصولات در فواصل

زمانی کوتاه­تر رخ می­دهد. بنابراین، برای ارزیابی جامع­تر عملکرد الگوریتم پیشنهادی، نسبت به مقادیر پارامترهای ، ، ،  و  تحلیل حساسیت صورت می­گیرد. بدین منظور، با در نظر گرفتن ،  و ، مطابق با جدول 8، 5 مقدار متفاوت به ترتیب برای پارامترهای ، ، ،  و  تعیین می­شود.

جدول 8. مقادیر پارامترهای مسأله در حالات مختلف برای تحلیل حساسیت

هریک از الگوریتم­ها برای حل هر حالت، 10 بار اجرا شده­اند که میانگین مقادیر تابع هدف بدست آمده و نیز درصد اختلاف میان میانگین مقادیر تابع هدف بدست آمده توسط هریک از الگوریتم­های معیار با الگوریتم پیشنهادی، در نمودارهای 12 الی 16 نشان داده شده است. این نکته قابل ذکر است که هر سه الگوریتم در تمامی حالت­ها به جوابی شدنی دست یافته­اند.

با توجه به نمودارهای مذکور، الگوریتم پیشنهادی تنها در یک مورد نسبت به الگوریتم ژنتیک از عملکرد اندکی ضعیف­تر برخوردار است (در تحلیل حساسیت نسبت به پارامتر  در حالت نخست) و در سایر موارد نسبت به دو الگوریتم معیار از عملکرد مناسب­تری برخوردار است. بنابراین، در کل می­توان عملکرد بهتر الگوریتم پیشنهادی را در ساختارهای مختلف مسایل نمونه نتیجه گرفت.

شکل 12. تحلیل حساسیت نسبت به پارامتر Q

شکل 13. تحلیل حساسیت نسبت به پارامتر

 شکل 14. تحلیل حساسیت نسبت به پارامتر

شکل 15. تحلیل حساسیت نسبت به پارامتر

شکل 16. تحلیل حساسیت نسبت به پارامتر