ارائۀ دو روش فراابتکاری برای حل مسئله مکان‌یابی هاب مرکز ظرفیت‌دار

نوع مقاله : مقاله پژوهشی

نویسندگان

1 استادیار گروه مهندسی صنایع، دانشکده مهندسی، دانشگاه کردستان، سنندج، ایران

2 دانشیار گروه مهندسی صنایع، دانشکده مهندسی، دانشگاه کردستان، سنندج، ایران

3 دانش‌آموخته کارشناسی ارشد مهندسی صنایع، دانشکده مهندسی، دانشگاه کردستان، سنندج، ایران

چکیده

مسئلۀ مکان‌یابی هاب یکی از موضوعات جدید در حوزۀ مسائل مکان‌یابی است. این دسته از مسائل، کاربردهای فراوانی در سیستم‌های حمل‌ونقل دارند. در این پژوهش، مسئلۀ مکان‌یابی تخصیص هاب مرکز با درنظرگرفتن محدودیت ظرفیت، بررسی می‌شود. هدف از مدل ارائه‌شده، یافتن مکان هاب‌ها و مشخص‌کردن نحوۀ تخصیص گره‌های غیرهاب به هاب است به‌گونه‌ای که بیشینه زمان سفر بین جفت گره‌های مبدأ-مقصد، کمینه شود. از آنجایی که مسئلۀ تحت بررسی از نوع مسائل ناچندجمله‌ای سخت محسوب می‌شود، در این تحقیق دو الگوریتم فراابتکاری شامل الگوریتم‌های بازپخت شبیه‌سازی شده و اجتماع مورچگان، توسعه داده می‌شود. کارایی الگوریتم‌ها از طریق حل تعدادی مثال عددی که از مجموعه مسائل شناخته‌شدۀ پست استرالیایی (AP) برگرفته شده، ارزیابی می‌شود و نتایج به‌دست‌آمده با راه‌حل‌های نرم‌افزار Lingo مقایسه می‌شوند. نتایج مثال‌ها حاکی از کارایی مناسب الگوریتم‌های توسعه ‌داده‌ شده است.

کلیدواژه‌ها

موضوعات


عنوان مقاله [English]

Two Metaheuristic Approaches for p-Hub Center Location Problem under Capacity Constraints

نویسندگان [English]

  • Alireza Eydi 1
  • Jamal Arkat 2
  • Ehsan Parhizkar Mehrabadi 3
1 Assistant Professor in Industrial Engineering, University of Kurdistan, Sanandaj, Iran
2 Associate Professor in Industrial Engineering, University of Kurdistan, Sanandaj, Iran
3 M.Sc. in Industrial Engineering, University of Kurdistan, Sanandaj, Iran
چکیده [English]

Hub location problem is one of the new issues in location problems. This kind of location problem is widely used in transportation systems. In this paper, we investigate p-hub center location allocation problem under capacity constraint. The aim of the proposed model is to determine the location of hub nodes and the allocation of non-hub nodes to the hub in such a way that the maximum traveling time is minimized. In addition, since the problem is an NP hard problem, two metaheuristic algorithms including simulated annealing algorithm and ant colony are developed for solving large size real world problems. The performances of the proposed algorithms are examined via some numerical examples taken from known related benchmark sets (AP dataset). The best solutions found using metaheuristic algorithms are also compare to the results achieved using Lingo software. The results demonstrate that the proposed algorithms are able to find optimum or near optimum solutions in acceptable run times.

کلیدواژه‌ها [English]

  • Simulated Annealing
  • Ant colony optimization
  • Capacity constraint
  • p-hub center location

مقدمه

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

در مسائل مکان‌یابی هاب، مکان استقرار هاب‌ها و نحوۀ تخصیص گره‌های شبکه به هاب‌ها به‌گونه‌ای تعیین می‌شوند که مجموع هزینه‌های استقرار هاب‌ها و هزینه‌های حمل‌ونقل در کل شبکه کمینه شوند. مسئلۀ مکان‌یابی هاب، کاربردهای فراوانی در مسائل دنیای واقعی همچون خطوط هوایی، سیستم‌های ترافیکی و حمل‌ونقل، خدمات پستی و شبکه‌های مخابراتی دارد. در تحقیقات اوکلی2 و میلر3 (1994)، کمپل4 و همکاران (2002) و فراهانی و همکاران (2013) بررسی جامعی در خصوص ادبیات این حوزه از مسائل مکان‌یابی انجام شده است. به‌طور کلی مسائل مکان‌یابی هاب را می‌توان به چهار نوع شامل مسئلۀ مکان‌یابی هاب میانه، مسئلۀ مکان‌یابی هاب با هزینۀ ثابت، مسئلۀ مکان‌یابی هاب مرکز و مسئلۀ مکان‌یابی هاب پوششی تقسیم کرد (الومار5 و کارا6، (2008)). از آنجایی که تحقیق حاضر، مسئلۀ مکان‌یابی هاب مرکز را بررسی می‌کند، در ادامه برخی از جدیدترین تحقیقات مرتبط با مسئلۀ هاب مرکز مورد بررسی بیشتر قرار می‌گیرند.

مسائل هاب مرکز از جمله مسائل کمینه‌سازی بیشینه (minimax) به شمار می‌روند. چنین مسائلی اغلب در سیستم‌های حمل‌ونقلی کاربرد دارند که در آنها زمان‌های پاسخ‌دهی به تقاضای مشتریان از اهمیت و حساسیت بالایی برخوردار است. مکان‌یابی تسهیلات اضطراری (اورژانسی) و مسیریابی حمل برای کالاهای فاسدشدنی و یا حساس به زمان نظیر مواد خوراکی یا دارویی از جمله چنین کاربردهایی هستند. چنین مسئله‌ای حتی در مواردی که بحث خدمت‌رسانی اضطراری مطرح نیست (همانند کمینه‌کردن بیشینه نارضایتی مسافران در سفرهای هوایی) نیز کاربرد دارد. هدف از مسئلۀ هاب مرکز، یافتن بهترین مکان استقرار هاب‌ها در شبکه و همچنین بهترین نحوۀ تخصیص گره‌ها به هاب‌ها است به‌گونه‌ای که بیشینه هزینه (زمان) سفر بین هر جفت از گره‌های مبدأ–مقصد، کمینه شود. اولین مدل‌های ریاضی برای مسئلۀ هاب مرکز را کمپل (1994) ارائه کرد. در این تحقیق سه نوع مدل ریاضی شامل کمینه‌کردن بیشینه هزینۀ بین هر زوج مبدأ-مقصد، کمینه‌کردن بیشینه هزینۀ جابه‌جایی بر هر یال (مبدأ به هاب، هاب به هاب و هاب به مقصد) و کمینه‌کردن بیشینه هزینۀ جابه‌جایی بین هاب و یک مبدأ یا یک مقصد، ارائه شده‌اند. کارا و تانسل7 (2000) یک مدل خطی عدد صحیح برای مسئلۀ هاب مرکز ارائه کردند و اثبات کردند که این مسئله ناچندجمله‌ای کامل8 است. مؤلفان همچنین نشان داده‌اند که مدل ارائه‌شده نسبت به نسخه‌های خطی دیگر برای مسئلۀ هاب مرکز از کارایی بالاتری برخوردار است. کمپل و همکاران (2007) مسئلۀ تخصیص هاب مرکز9 را بررسی کرده‌اند. مسئلۀ تخصیص هاب حالت خاصی از مسئلۀ مکان‌یابی هاب است که در آن، مکان‌های استقرار هاب‌ها از قبل مشخص هستند و صرفاً نحوۀ تخصیص گره‌های غیرهاب به گره‌های هاب باید مشخص شود. محققان چند مدل عدد صحیح برای این مسئله ارائه کرده‌اند و پیچیدگی محاسباتی هر یک از مدل‌ها را بررسی کرده‌اند.

کاربردهای متعدد مسئلۀ هاب مرکز در حوزه‌های مختلف باعث شده است الگوریتم‌های دقیق و ابتکاری متعددی برای این مسئله توسعه داده شود. به‌طور متعارف، برای حل مسائل با اندازۀ کوچک از روش‌های دقیق بهینه‌سازی مانند شاخه و کرانه، شاخه و برش و تجزیۀ بندرز10 استفاده می‌شود. این در حالی است که به‌دلیل ناچندجمله‌ای سخت‌بودن مسائل مکان‌یابی هاب، برای مسائل با ابعاد بزرگ، از روش‌های ابتکاری و فراابتکاری مختلفی نظیر جستجوی حریص11، جستجوی ممنوعه12 (TS)، بازپخت شبیه‌سازی‌شده13 (SA) و الگوریتم ژنتیک14استفاده شده‌اند (فراهانی و حکمت‌فر 2009).

نخستین روش حل ابتکاری برای مسئلۀ تخصیص هاب مرکز را پاموک15 و سپیل16 (2001) ارائه کرده‌اند. روش ابتکاری ارائه‌شده ترکیبی از الگوریتم جستجوی ممنوعه و یک الگوریتم جستجوی حریص است و در رویه‌ای ترتیبی، دو مسئلۀ مکان‌یابی و تخصیص حل می‌شوند؛ بدین معنی که نخست مکان‌های هاب‌ها با استفاده از رویکردهای ابتکاری تعیین می‌شوند، سپس با توجه به مکان‌های هاب‌ها، نحوۀ تخصیص گره‌های غیرهاب به هاب‌ها مشخص می‌شوند. در گام بعد، الگوریتم حریص با توجه به تخصیص به‌دست‌آمده، مکان‌های استقرار هاب‌ها را تغییر می‌دهد و سپس مجدداً مسئلۀ تخصیص با توجه به راه‌حل به‌دست‌آمده برای مکان‌یابی، اصلاح می‌شود. رویۀ ترتیبی حل مسئله مکان‌یابی مسئلۀ تخصیص تا رسیدن به یک راه‌حل مطلوب تکرار می‌شود. الگوریتم حلی که می‌یر17 و همکاران (2009) ارائه کرده‌اند یک الگوریتم ترکیبی است. این الگوریتم برای حل مسئلۀ هاب مرکز در حالت تخصیص تکی توسعه داده شده است و مشتمل بر دو الگوریتم شامل الگوریتم کوتاه‌ترین مسیر (فاز اول) و الگوریتم بهینه‌سازی اجتماع مورچگان18 (فاز دوم) است. شایان ذکر است که محققان در فاز اول الگوریتم توسعه‌داده‌شده از یک رویۀ شاخه و کرانه استفاده کرده‌اند.

ازجمله مطالعاتی که درصدد توسعۀ الگوریتم‌های دقیق برای حل مسئلۀ هاب مرکز بوده‌اند می‌توان به بامگارتنر19 (2003)، همچیر20 و می‌یر (2006) و ارنست21 و همکاران (2009) اشاره کرد. بومگارتنر (2003) مدل‌های ریاضی مختلفی را که برای مسئلۀ هاب مرکز بدون محدودیت ظرفیت ارائه شده‌اند، مرور کرده است و سپس با معرفی حدود بالا و پایین کارا، یک الگوریتم شاخه و کرانه را توسعه داده و آن را با نسخه‌های متعارف الگوریتم شاخه و کرانه مقایسه کرده است. همچیر و مایر (2006) نشان داده‌اند که الگوریتم ترکیبی جستجوی دودویی22 که کارایی آن برای مسائل متعارف مکان‌یابی مرکز اثبات شده است در مسائل مکان‌یابی هاب مرکز نیز به‌خوبی عمل می‌کند. ارنست و همکاران (2009) مسئلۀ هاب مرکز بدون ظرفیت را برای دو حالت تخصیص تکی و تخصیص چندگانه بررسی کرده‌اند. مؤلفان نشان داده‌اند که مسئله در هر دو حالت ناچندجمله‌ای سخت است و حتی در حالتی که مکان استقرار هاب‌ها نیز مشخص باشد و هدف صرفاً یافتن نحوۀ تخصیص گره‌های غیرهاب باشد، مسئله همچنان ناچندجمله‌ای سخت باقی می‌ماند. محققان برای حالت تخصیص چندگانه، یک الگوریتم شاخه و کرانۀ کارا توسعه داده‌اند.

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

در پژوهش حاضر، یک مدل برنامه‌ریزی ریاضی عدد صحیح برای مسئلۀ مکان‌یابی هاب مرکز با ظرفیت محدود ارائه می‌شود. هدف مدل ریاضی ارائه‌شده، تعیین مکان‌های استقرار هاب‌ها و نحوۀ تخصیص گره‌های غیرهاب است؛ به‌گونه‌ای که بیشینه زمان سفر بین هر جفت گره مبدأ–مقصد، کمینه شود. از آنجایی که مسئلۀ تحت بررسی، ناچندجمله‌ای سخت است، برای حل آن، دو الگوریتم فراابتکاری شامل الگوریتم بازپخت شبیه‌سازی‌شده و بهینه‌سازی اجتماع مورچگان، ارائه می‌شوند. برای ارزیابی کارایی الگوریتم‌های ارائه‌شده، تعدادی از مسائل مربوط به مجموعه مسائل شناخته‌شدۀ پست استرالیا23 (AP) توسط دو الگوریتم حل می‌شوند و نتایج حل، بررسی می‌شوند. در ادامه، ساختار مقاله بدین صورت سازمان‌دهی شده است:

نخست، مسئلۀ تحت بررسی، تشریح می‌شود و در قالب یک مدل عدد صحیح، مدل‌سازی می‌شود. در بخش بعد، ساختار الگوریتم‌های حل پیشنهادی تشریح می‌شود. سپس مثال‌های عددی و بررسی نتایج ارائه می‌شود. در آخر، جمع‌بندی، نتیجه‌گیری و ارائۀ پیشنهادها برای مطالعات آتی در بخش نهایی ارائه می‌شوند.

 

2- بیان مسئله و ارائۀ مدل ریاضی

شکل کلی مسئله تحت بررسی بدین صورت است: مجموعۀ ، شبکه‌ای با  گره را نشان می‌دهد، مجموعۀ  نیز مجموعه هاب‌های بالقوه را که شامل  گره است نشان می‌دهد. زیرشبکۀ هاب‌ها یک شبکۀ کامل است؛ بدین معنی که هر جفت هاب به‌صورت مستقیم به یکدیگر متصل هستند. ارتباط جفت گره‌های غیرهاب صرفاً از طریق هاب‌ها امکان‌پذیر است، بدین معنی که بین گره‌های غیرهاب ارتباط مستقیم وجود ندارد. هدف از مسئلۀ تحت بررسی، انتخاب تعداد مشخصی از گره‌های بالقوه برای احداث هاب و تخصیص گره‌های غیرهاب به هاب‌ها است؛ به‌نحوی که بیشینه زمان‌های سفر بین هر جفت مبدأ–مقصد، کمینه شود. برای مدل‌سازی این مسئله از مفهوم جدید شعاع هاب24(می‌یر و همکاران، 2009) استفاده می‌شود. برای روشن‌شدن این مفهوم، موقعیت نسبی دو هاب و تعدادی از گره‌های غیرهاب تخصیص‌یافته به آنها در شکل (1) نشان داده شده است. در این شکل، هاب‌ها با مربع و گره‌های غیرهاب با دایره مشخص شده‌اند. برای به‌دست‌آوردن مدت زمان سفر بین یک مبدأ در شعاع پوشش یکی از هاب‌ها و یک مقصد در شعاع پوشش هاب دیگر، مجموع شعاع‌های پوشش دو هاب و ضریبی از فاصلۀ زمانی بین دو هاب در نظر گرفته می‌شود. فرض بر این است که سرعت انتقال مستقیم بین دو هاب بیشتر از سرعت انتقال از یک گره غیرهاب به یک هاب یا برعکس است.

 شکل 1- مثالی از مفهوم شعاع هاب

 

مفروضات زیر در مدل‌سازی مدنظر قرار می‌گیرند:

  • هر گره غیرهاب به یک هاب تخصیص می‌یابد.
  • به‌دلیل صرفۀ اقتصادی در ارتباط مستقیم دو هاب از فاکتور تنزیل (اوکلی، 1987) استفاده می‌شود.
  • میزان جریان عبوری از هر هاب، محدود است.
  • یال‌های شبکه فاقد محدودیت ظرفیت هستند.
  • گراف مربوط به گره‌های غیرهاب و هاب‌ها، متقارن است و نامعادلۀ مثلثی برای فواصل زمانی در آن صدق می‌کند.
  • کلیۀ پارامترهای مسئله، قطعی هستند.

نمادگذاری مدل ریاضی به شرح زیر است:

مجموعۀ اندیس‌ها

: مجموعۀ گره‌ها (  و : اندیس گره)

: مجموعۀ هاب‌ها (  و : اندیس گره)

 

پارامترها

: شعاع پوشش هاب

: هزینۀ (زمان یا مسافت) سفر از گره  به گره

: ظرفیت هاب

: مقدار تقاضای (جریان خروجی از) گره

: ضریب تنزیل

 

متغیرهای تصمیم

: متغیر دودویی که مقدار یک می‌گیرد اگر یک هاب در گره  مستقر شود.

: متغیر دودویی که مقدار یک می‌گیرد اگر گره  به هاب مستقر در گره  تخصیص داده شود.

با درنظرگرفتن مفروضات عنوان‌شده و نمادگذاری فوق، مدل ریاضی مسئله به‌شکل زیر ارائه می‌شود:

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

تابع هدف مدل ریاضی ارائه‌شده، غیرخطی است. برای خطی‌کردن این عبارت، پس از تعریف متغیر تصمیم جدید ، تابع هدف به‌صورت رابطۀ (9) بازنویسی و محدودیت (10) به مدل اضافه می‌شود.

 

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

 

3- الگوریتم‌های حل

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

 

3-1- الگوریتم بازپخت شبیه‌سازی‌شده

الگوریتم بازپخت شبیه‌سازی‌شده از فرایند فیزیکی بازپخت یا تبرید الهام گرفته شده است. در فرایند بازپخت، فلز یا آلیاژ مذاب به‌آرامی و با کاهش تدریجی دما، سرد می‌شود تا ساختار ملکولی مناسبی برای محصول نهایی به دست آید. هدف بازپخت شبیه‌سازی‌شده، یافتن ساختار یا چیدمان بهینه (حالتی با کمترین انرژی) برای یک مسئلۀ پیچیده است (بلوم25 و رولی26، 2003 و گلور27 و کوچنبرگر28، 2003). نقطۀ شروع این الگوریتم برخلاف الگوریتم‌های تکاملی، یک راه‌حل اولیه است. راه‌حل اولیه می‌تواند به‌صورت تصادفی یا براساس رویه‌ای نظام‌مند تولید شود. SA دارای دو حلقۀ اصلی شامل حلقۀ درونی و حلقۀ بیرونی است. در هر تکرار حلقۀ درونی، یک راه‌حل همسایه برای راه‌حل فعلی تولید می‌شود. در صورتی که راه‌حل همسایه بهتر از راه‌حل فعلی باشد، جایگزین آن می‌شود. حتی در صورتی که راه‌حل همسایه بهتر از راه‌حل کنونی نباشد، با احتمال مشخصی، شانس جایگزین‌شدن با راه‌حل فعلی را دارد. این ویژگی باعث می‌شود که الگوریتم تا حدود زیادی بتواند از دام تله‌های بهینگی محلی بگریزد. حلقۀ درونی تا یافتن تعداد مشخصی راه‌حل بهبوددهنده یا تعداد مشخصی تکرار (طول زنجیره مارکوف29) تکرار می‌شود. پس از آنکه شرط خروج از حلقۀ درونی برقرار شد، دما در حلقۀ بیرونی کاهش می‌یابد و تکرارهای حلقۀ درونی مجدداً تکرار می‌شوند. فرایند کاهش دما در حلقۀ بیرونی باعث کاهش احتمال پذیرش راه‌حل‌های غیربهبوددهنده می‌شود. الگوریتم پس از رسیدن به دمای مشخصی که دمای انجماد نامیده می‌شود، متوقف می‌شود. عناصر اصلی SA در ادامه تشریح می‌شوند.

3-1-1- ساختار راه‌حل و تولید راه‌حل اولیه

اولین مرحله در هر الگوریتم فراابتکاری، طراحی ساختاری مناسب برای ارائۀ راه‌حل است. ساختار درنظرگرفته‌شده، یک بردار دوبخشی است. بخش اول که اندازۀ آن برابر با تعداد هاب‌ها است، شمارۀ گره‌هایی را نشان می‌دهد که برای استقرار هاب‌ها، انتخاب شده‌اند. اندازۀ بخش دوم برابر با تعداد گره‌های شبکه است و اعداد آن، شمارۀ هاب‌هایی را نشان می‌دهند که گره‌ها به آنها تخصیص یافته‌اند. شکل (2) نمونه‌ای از چنین ساختاری را برای یک مسئله با 10 گره و سه هاب نشان می‌دهد. در راه‌حلی که این بردار نشان می‌دهد، هاب‌های 1،‌ 2 و 3 به‌ترتیب در گره‌های 6، 3 و 7 مستقر هستند. پیش از آنکه نحوۀ تخصیص مشخص شود، لازم است بردار تخصیص (بخش دوم راه‌حل) در صورت نیاز اصلاح شود. همان‌گونه که در این مثال دیده می‌شود اگرچه هاب‌های شماره 1 و 3 به‌ترتیب در گره‌های 6 و 7 مستقر شده‌اند، گره 6 به هاب شماره 2 و گره 7 به هاب شماره 1 تخصیص یافته است. با این اوصاف باید این راه‌حل اصلاح شود. بردار اصلاح‌شدۀ شکل(2) در شکل (3) نشان داده شده است. با توجه به بردار اصلاح‌شده نحوۀ تخصیص بدین صورت است: گره‌های 1، 5، 6، 8 و 9 به هاب شماره 1 (مستقر در گره 6)، گره‌های 2، 3، 4 و 10 به هاب شماره 2 (مستقر در گره 3) و گره 7 به هاب شماره 3 (مستقر در گره 7) تخصیص داده می‌شوند.

 

2

1

1

1

2

1

2

2

2

1

7

3

6

شکل 2- نمونه‌ای از بردار راه‌حل (پیش از اصلاح)

 

2

1

1

3

1

1

2

2

2

1

7

3

6

شکل 3- راه‌حل اصلاح‌شدۀ متناظر بردار شکل 2

راه‌حل اولیۀ SA به‌صورت تصادفی تولید می‌شود بدین صورت که اعداد مربوط به بخش اول به‌صورت تصادفی از بین اعداد 1 تا تعداد گره‌ها انتخاب می‌شوند و اعداد بخش دو نیز به‌صورت تصادفی از بین اعداد 1 تا تعداد هاب‌ها انتخاب می‌شوند. در صورت نشدنی‌بودن، براساس شیوه‌ای که گفته شد بخش دوم راه‌حل باید اصلاح شود.

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

 

3-1-2- تولید راه‌حل همسایه

شیوه‌ای که SA برای تولید راه‌حل همسایه از راه‌حل فعلی به کار می‌گیرد در سرعت همگرایی و کیفیت راه‌حل نهایی الگوریتم تأثیر فراوانی دارد. در الگوریتم ارائه‌شده، از دو رویۀ متفاوت با احتمالات مساوی برای تولید همسایه استفاده می‌شود. رویۀ نخست برای بخش اول بردار راه‌حل (انتخاب هاب‌ها) و رویۀ دوم برای بخش دوم (تخصیص گره‌ها به هاب‌ها) استفاده می‌شود. در رویۀ نخست، یکی از درایه‌های بخش اول بردار راه‌حل به‌صورت تصادفی انتخاب می‌شود و مقدار آن، به‌صورت تصادفی به مقداری دیگر تغییر داده می‌شود. چنین تغییری معادل انتخاب گرهی دیگر برای استقرار هاب است. در رویۀ دوم، از تعویض‌های دوتایی استفاده می‌شود بدین معنی که دو درایۀ متناظر با گره‌های غیرهاب، از بخش دوم بردار راه‌حل به‌صورت تصادفی انتخاب می‌شوند و مقادیر آنها (در صورتی که متفاوت باشند) معاوضه می‌شوند. در صورتی که مقادیر درایه‌های انتخاب‌شده یکسان باشند، دو درایۀ دیگر به‌صورت تصادفی انتخاب می‌شوند.

 

3-1-3- تنظیم پارامترهای الگوریتم

تنظیم پارامترها یکی از مراحل مهم در به‌کارگیری یک الگوریتم فراابتکاری است. در الگوریتم SA دمای اولیه، دمای انجماد، ضریب کاهش دما30 و طول زنجیرۀ مارکوف، مهم‌ترین پارامترهای الگوریتم محسوب می‌شوند. دمای اولیه باید به‌گونه‌ای تنظیم شود که در تکرارهای نخست الگوریتم، احتمال پذیرفته‌شدن یک راه‌حل غیربهبوددهنده، نزدیک به یک (معمولاً حداقل 90 درصد) باشد. در این تحقیق، احتمال پذیرش یک راه‌حل همسایه غیربهبوددهنده از رابطۀ (11) به دست می‌آید:

 

که در آن،  میزان اختلاف بین راه‌حل همسایه و راه‌حل فعلی را نشان می‌دهد و مقدار دمای فعلی است. با درنظرگرفتن رابطۀ فوق، می‌توان دمای اولیه را براساس شرط نشان‌داده‌شده در رابطۀ (12) تنظیم کرد:

 

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

مقدار دما در ابتدای هر حلقۀ بیرونی از طریق ضرب‌کردن دمای فعلی در فاکتور کاهش دما، کاهش داده می‌شود. در صورتی که دما به‌آرامی کاهش یابد، احتمال به‌دست‌آوردن راه‌حل‌هایی با کیفیت بالاتر، بیشتر خواهد شد؛ اما در مقابل زمان محاسبات نیز بیشتر خواهد شد. از سویی، اگر دما به‌سرعت کاهش یابد، ممکن است الگوریتم در تلۀ بهینگی محلی قرار گیرد. در این تحقیق از مقادیر 8/0 تا 99/0 به‌عنوان ضریب کاهش دما استفاده شده است.

تعداد تکرارهای حلقۀ درونی به طول زنجیره مارکوف اشاره دارند. واضح است که هرچه طول زنجیرۀ مارکوف بیشتر باشد، جستجوی محلی بیشتر انجام می‌گیرد و شانس دستیابی به راه‌حل‌های مرغوب‌تر بیشتر است؛ اما در مقابل، زمان حل نیز افزایش خواهد یافت. طول زنجیرۀ مارکوف با توجه به نوع مسئله و تعداد متغیرهای تصمیم مسئله، تعیین می‌شود.

 

3-2- الگوریتم کلونی مورچگان

بهینه‌سازی کلونی مورچگان یک روش احتمالی برای حل مسائل بهینه‌سازی ترکیباتی است که براساس مشاهدۀ رفتار مورچه‌ها توسعه یافته است. بهینه‌سازی دستۀ مورچگان یک روش جمعیت محور است که می‌تواند راه‌حل‌های مناسبی را برای مسائل پیچیدۀ بهینه‌سازی به دست آورد (بلوم ورولی،2003). مطالعات زیست‌شناسی نشان داده است که مورچه‌ها ماده‌ای به نام فرومون31 از خود ترشح می‌کنند که به‌وسیلۀ آن می‌توانند با یکدیگر ارتباط برقرار کنند و اطلاعات مسیر را به‌هنگام جستجوی غذا در اختیار سایرین قرار دهند (گلوورو کوچنبرگ،2003). به‌طور کلی، ساختار الگوریتم ACO را می‌توان در چهار مرحلۀ اصلی خلاصه کرد. در مرحلۀ نخست، پارامترهای اولیۀ الگوریتم (شامل تعداد مورچه‌ها، سرعت تبخیر و ...) تنظیم می‌شوند. در مرحلۀ دوم هر مورچه، یک راه‌حل اولیه (شدنی) را تولید می‌کند. مراحل سوم (جستجوی محلی) و چهارم (به‌هنگام‌کردن) در یک حلقه تکرارشونده، انجام می‌شوند. در مرحلۀ سوم، با استفاده از یک رویۀ ابتکاری، راه‌حل‌های فعلی که توسط مورچه‌ها تولید شده‌اند، بهبود داده می‌شوند. در مرحلۀ چهارم، مقادیر فرومون‌های ترشح‌شده در هر مسیر براساس توابعی مشخص، به‌هنگام می‌شوند. دو مرحلۀ جستجوی محلی و به‌هنگام‌سازی تا جایی تکرار می‌شوند که الگوریتم به شرط توقف برسد. جزئیات این مراحل در ادامه تشریح شده‌اند.

 

3-2-1- مقداردهی اولیۀ فرومون‌ها

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

در ادامه توضیح داده خواهد شد که چگونه مقادیر فرومون‌ها چنین انتخاب‌هایی را کنترل می‌کنند و چگونه فرومون‌ها در طی مراحل الگوریتم به‌هنگام می‌شوند.

 

3-2-2- حلقۀ اصلی الگوریتم

در این مرحله تعداد تکرارها و تعداد مورچه‌ها تعیین می‌شود. تعداد  مورچه برای شروع الگوریتم تعیین می‌شوند و هر کدام در طی الگوریتم، راه‌حل جزئی را کامل می‌کنند. ساختن یک راه‌حل شامل دو مرحله است: در مرحلۀ نخست، هر مورچه ابتدا مکان هاب‌ها را مشخص می‌کند و در مرحلۀ دوم، سایر گره‌ها به هاب‌های انتخاب‌شده تخصیص می‌یابند. در مرحلۀ اول تعداد  هاب انتخاب می‌شود. یک گره به‌صورت تصادفی با احتمال زیر به‌عنوان هاب انتخاب می‌شود:

 

اگر مقدار  از مقدار  (یک عدد تصادفی انتخاب‌شده از بازه صفر تا یک) بزرگ‌تر باشد آنگاه گره  به‌عنوان یک هاب انتخاب می‌شود و در غیر این صورت  به‌عنوان یک هاب جدید انتخاب می‌شود. هر زمان یک هاب انتخاب می‌شود، عمل تبخیر موضعی32 صورت می‌گیرد و مقدار فرومون براساس رابطۀ (14) به‌هنگام می‌شود:

 

که در آن ضریب  نشان‌دهندۀ نرخ عمل تبخیر است. عمل تبخیر باعث می‌شود که در راه‌حل‌های بعدی، احتمال انتخاب هاب‌های گفته‌شده کمتر شود و بنابراین تنوع‌پذیری33 الگوریتم (جستجوی فضایی گسترده‌تر) افزایش می‌یابد. پس از تعیین هاب‌ها، در مرحلۀ دوم ساختن راه‌حل، سایر گره‌های باقیمانده به هاب‌ها تخصیص می‌یابند. برای انتخاب دورترین گره برای تخصیص به یک هاب، از روشی مشابه با آنچه در بخش قبل گفته شد، استفاده می‌شود با این تفاوت در این مرحله از فاکتوری تحت عنوان تمایل ابتکاری34 استفاده می‌شود. نحوۀ تخصیص بدین صورت است که صرفاً گره‌هایی تخصیص می‌یابند که منجر به افزایش مقدار تابع هدف نشوند. مورچه‌ها گره موردنظر را براساس رابطۀ احتمالی (15) انتخاب می‌کنند:

 

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

 

در رابطۀ فوق پارامتر  تأثیر  را محدود می‌کند. مقدار تمایل ابتکاری ( ) از رابطۀ (17) به دست می‌آید:

 

این ضریب در واقع نشان می‌دهد که گره‌هایی که در فاصلۀ دورتری نسبت به هاب  قرار دارند به این هاب تخصیص نمی‌یابند تا زمانی که گره‌های نزدیک‌تر به هاب  تخصیص یابند. واضح است که متوسط فاصلۀ یک گره غیرهاب از یک گره هاب با افزایش تعداد هاب‌ها کاهش می‌یابد. هنگامی که یک گره غیرهاب انتخاب می‌شود تبخیر موضعی فرومون براساس رابطۀ (18) انجام می‌گیرد:

 

اگر گرهی پس از مرحلۀ قبل، تخصیص نیابد با احتمالی که توسط فرومون  تعیین می‌شود، هاب  انتخاب می‌شود و گره  به آن تخصیص می‌یابد. در چنین شرایطی، دسته‌ای از پارامترهای ابتکاری پویا  که بستگی به مجموعۀ فعلی هاب‌ها دارد براساس رابطۀ (19) انتخاب می‌شود:

 

احتمال انتخاب هاب  و تخصیص گره  به آن براساس رابطۀ (20) محاسبه می‌شود:

 

در الگوریتم ارائه‌شده مقدار 5/0 برای  در نظر گرفته شده است. صرف‌نظر از اینکه  چگونه انتخاب شود، تبخیر موضعی مطابق معادلۀ (18) انجام می‌گیرد. همچنین اگر مشخص شود که جواب فعلی نشدنی است هاب‌های انتخاب‌شده نیز به‌صورت موضعی و با مقدار ثابت 5/0 تبخیر می‌شوند.

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

 

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

 

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

در این بخش، به‌منظور بررسی کارایی الگوریتم‌های ارائه‌شده در حل مسائل بزرگ (دنیای واقعی) تعدادی مثال عددی ارائه می‌شود. مثال‌های ارائه‌شده از مجموعه مسائل نمونه AP (شرکت پست استرالیا) استخراج شده‌اند. این مجموعه از داده‌ها از مطالعۀ سیستم ارسال مرسولات کشور استرالیا به دست آمده و برای نخستین بار در تحقیقات ارنست و کریشنامورتی36 (1997) استفاده شده است. به‌منظور انجام مقایسات، مدل ریاضی با استفاده از نرم‌افزار Lingo 8 و بر رایانۀ شخصی با مشخصات Pentium 4, 2.53GHz و 4GB RAM اجرا شده است. همچنین الگوریتم‌های حل معرفی‌شده با استفاده از نرم‌افزار برنامه‌نویسی C# کدنویسی شده‌اند. برای اجرای مدل ریاضی بر نرم‌افزار Lingo سقف دو ساعت در نظر گرفته شده است به عبارتی در صورتی که نرم‌افزار نتواند ظرف دو ساعت، راه‌حل بهینه را تولید کند، بهترین جواب به‌دست‌آمده، گزارش می‌شود. در این تحقیق نرخ تنزیل ( ) برابر با 75/0 در نظر گرفته شده است و نمونه‌های مختلفی با شبکه‌های دارای 20، 25، 40 و 50 گره و 2، 3، 4 و 5 هاب بررسی و تحلیل شده‌اند.

تنظیم پارامترها در کیفیت راه‌حل نهایی الگوریتم‌های فراابتکاری تأثیر قابل‌ملاحظه‌ای دارد. بدین جهت پیش از آنکه از الگوریتم‌های ارائه‌شده در حل مثال‌های عددی استفاده شود، لازم است پارامترهای آنها تنظیم شود. در الگوریتم SA دو پارامتر نرخ انجماد و طول زنجیره مارکوف، پارامترهای مهم به شمار می‌روند. نرخ انجماد در الگوریتم ارائه‌شده برای مثال‌های مختلف بین اعداد 90/0 تا 95/0 انتخاب شده است. همچنین از آنجایی که نتایج محاسباتی حاکی از آن است که افزایش طول زنجیره مارکوف، تأثیری بر کیفیت راه‌حل نهایی ندارد و صرفاً باعث افزایش زمان حل الگوریتم می‌شود، این پارامتر برابر عدد یک در نظر گرفته شده است. به عبارتی در هر سطح دما با دستیابی به یک راه‌حل همسایۀ بهبوددهنده، حلقۀ درونی به اتمام می‌رسد و دما کاهش می‌یابد.

برای تنظیم پارامترهای الگوریتم ACO نیز پس از انجام تعدادی آزمایش، این مقادیر به‌عنوان مقادیر مناسب انتخاب شده‌اند:

 

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

جدول 1- نتایج حل مثال‌های عددی

اندازه مسئله

Lingo

ACO

SA

BS

t (s)

BS

t (s)

BS

t (s)

(2) 20

9/45

40

9/45

4

9/45

0

(3) 20

4/43

31

4/43

4

4/43

0

(4) 20

6/38

11

6/38

5

6/38

1

(5) 20

9/37

15

9/37

5

9/37

1

(2) 25

1/57

41

1/57

9

1/57

2

(3) 25

6/46

110

6/46

10

6/46

2

(4) 25

5/45

46

5/45

11

5/45

2

(5) 25

5/45

23

5/45

10

5/45

3

(2) 40

4/65

366

2/66

91

2/66

7

(3) 40

7/58

1180

2/62

43

4/60

7

(4) 40

-

-

1/55

47

1/58

7

(5) 40

7/49

710

8/55

30

6/52

7

(2) 50

4/75

4614

6/87

273

3/91

10

(3) 50

-

-

5/80

290

6/89

10

(4) 50

-

-

8/80

192

9/88

11

(5) 50

-

-

9/81

161

9/87

11

                 

 

در جدول فوق علامت «-» بدین معنی است که نرم‌افزار Lingo نتوانسته است در مدت‌زمان سه ساعت، جواب بهینه را برای مسئلۀ متناظر به دست آورد. ستون اول جدول، معرف اندازۀ مسئله (تعداد هاب.گره) و ستون دوم، راه‌حل به‌دست‌آمده (بهترین جواب و زمان حل) توسط Lingo است. ستون‌های بعدی به‌ترتیب جواب‌های به‌دست‌آمده از روش‌های ACO و SA هستند.

  

شکل 4- نتایج حل مثال‌های عددی

به‌منظور تفسیر شکلی نتایج، مقایسه این سه روش از نظر مقدار تابع هدف، در شکل (4) نشان داده شده است. همان‌گونه که مشاهده می‌شود برای مسائل کوچک با اندازه‌های 20 و 25، راه‌حل‌های به‌دست‌آمده توسط هر سه روش، بهینه هستند و برای مسائل با اندازۀ بزرگ‌تر، کیفیت راه‌حل‌های به‌دست‌آمده از ACO از SA بهتر بوده است. همچنین برای مسائل با اندازۀ 50 گره و تعداد سه هاب یا بیشتر نرم‌افزارLingo نتوانسته است در مدت‌زمان سه ساعت به راه‌حل بهینه دست یابد. با بزرگ‌شدن ابعاد مسئله، زمان حل الگوریتم ACO بیشتر از الگوریتم SA شده است. نتایج فوق حاکی از کارایی مناسب دو الگوریتم فراابتکاری است که در این تحقیق ارائه شدند. همان‌گونه که جدول (1) نشان می‌دهد افزایش تعداد هاب‌ها در تعداد ثابتی از گره‌ها، مقدار تابع هدف را کاهش می‌دهد و تأثیر محدودیت ظرفیت با افزایش تعداد هاب‌ها کمتر می‌شود.

 

5- نتیجه‌گیری و ارائۀ پیشنهادها

در این مقاله، مسئلۀ مکان‌یابی هاب مرکز با درنظرگرفتن محدودیت ظرفیت بررسی شد. به‌منظور لحاظ‌کردن خدمت‌دهی به دورترین نقاط تقاضا در کمترین زمان خدمت، برای مسئلۀ تحت بررسی یک مدل ریاضی کمینه‌سازی بیشینه ارائه شد. با توجه به آنکه مسئله تحت بررسی از مسائل ناچندجمله‌ای سخت به شمار می‌رود، برای حل مسئله در اندازه‌های بزرگ، دو الگوریتم فراابتکاری شامل یک الگوریتم بازپخت شبیه‌سازی‌شده و یک الگوریتم بهینه‌سازی اجتماع مورچگان، توسعه داده شدند. با استفاده از تعدادی از مثال‌های عددی به‌دست‌آمده از مجموعه مسائل معتبر AP، کارایی الگوریتم‌های توسعه‌داده‌شده در مقابل راه‌حل‌های به‌دست‌آمده از Lingo به‌عنوان نرم‌افزار بهینه‌ساز، ارزیابی شد. نتایج حاصل از مثال‌های حل‌شده حاکی از آن است که الگوریتم‌های توسعه‌داده‌شده می‌توانند با صرف زمان‌های اجرای بسیار کم، راه‌حل‌هایی با کیفیت مناسب تولید کنند. در مقایسۀ دو الگوریتم ACO و SA از آنجایی که مسائل مکان‌یابی جزء مسائل طراحی محسوب می‌شوند و میزان اهمیت کیفیت راه‌حل نهایی به‌مراتب از میزان زمان اجرای صرف‌شده بیشتر است، الگوریتم ACO کارایی بالاتری نسبت به الگوریتم SA داشته است.

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

Alumur, S., & Kara, B. Y. (2008). "Network hub location problems: The state of the art". European Journal of Operational Research190(1), 1-21.
Bashiri, M., &Mirzaei, M. (2008). "Hybrid Fuzzy Capacitated Hub Center Allocation Problem with bothqualitative and quantitative variables". World Applied Sciences Journal5(4), 507-516.
Baumgartner, S. (2003). Polyhedral analysis of hub center problems (Doctoral dissertation, TechnischeUniversität Kaiserslautern).
Blum, C., &Roli, A. (2003). "Metaheuristics in combinatorial optimization: Overview and conceptual comparison". ACM Computing Surveys (CSUR)35(3), 268-308.
Campbell, J. F. (1994). "Integer programming formulations of discrete hub location problems". European Journal of Operational Research72(2), 387-405
Campbell, J.F., Ernst, A.T.,Krishnamoorthy, M. (2002) Hub location problems, Drezner, Z. andHamacher,H.W., Eds., "Facility Location: Applications and Theory", New York,Springer-Verlag, 373-408.
Campbell, A. M., Lowe, T. J., & Zhang, L. (2007). "The p-hub center allocation problem". European Journal of Operational Research176(2), 819-835.
Ernst, A. T., Hamacher, H., Jiang, H., Krishnamoorthy, M., &Woeginger, G. (2009). "Uncapacitated single and multiple allocation p-hub center problems". Computers & Operations Research36(7), 2230-2241.
Ernst, A. T., &Krishnamoorthy, M. (1996). "Efficient algorithms for the uncapacitated single allocation p-hub median problem". Location science4(3), 139-154.
Farahani, R. Z., Hekmatfar, M., Arabani, A. B., &Nikbakhsh, E. (2013). "Hub location problems: A review of models, classification, solution techniques, and application". Computers& Industrial Engineering64(4), 1096-1109.
Farahani,R.Z., Hekmatfar,M., (2009). Facility Location: Concepts, Models, Algorithms and Case Studies, Germany, physica-verlag, pp.243-270.
Glover, F., Kochenberger, G.A. (2003). Handbook of metaheuristics, Springer.
Hamacher, H. W., & Meyer, T. (2006). Hub cover and hub center problems. TechnischeUniversität Kaiserslautern, FachbereichMathematik.
Kara, B. Y., &Tansel, B. C. (2000). "On the single - assignment p-hub center problem". European Journal of Operational Research, 125(3), 648-655.
Meyer, T., Ernst, A. T., &Krishnamoorthy, M. (2009). "A 2-phase algorithm for solving the single allocation p -hub center problem". Computers & Operations Research36(12), 3143-3151.
 
O'kelly, M. E. (1987). "A quadratic integer program for the location of interacting hub facilities". European Journal of Operational Research32(3), 393-404.
O'Kelly, M. E., & Miller, H. J. (1994). "The hub network design problem: a review and synthesis". Journal of Transport Geography, 2(1), 31-40.
Pamuk, F. S., & Sepil, C. (2001). A solution to the hub center problem via a single-relocation algorithm with tabu search. IIE Transactions, 33(5), 399-411.