الگوریتم بهینه سازی ازدحام ذرات – کد الگوریتم PSO | راهنمای جامع

در این مطلب، الگوریتم بهینه سازی ازدحام ذرات (Particle Swarm Optimization | PSO) به طور کامل و همراه با مثال مورد بررسی قرار گرفته و پیاده‌سازی الگوریتم PSO در پایتون، متلب و جاوا انجام شده است. شایان توجه است که به منظور تشریح محاسبات ریاضی نهفته در پس الگوریتم بهینه سازی ازدحام ذرات یا همان الگوریتم PSO از نسخه کلاسیک این الگوریتم استفاده خواهد شد. در مطلب «الگوریتم بهینه سازی ازدحام ذرات | کد الگوریتم PSO در پایتون ، متلب و جاوا | راهنمای جامع» ابتدا به مفهوم بهینه‌سازی پرداخته شده و سپس، الگوریتم بهینه سازی ازدحام ذرات به طور جامع و کامل مورد بررسی قرار گرفته است. در ادامه مطلب، انواع الگوریتم بهینه سازی ازدحام ذرات تشریح می‌شود. همچنین، روش‌های ترکیبی موجود با بهره‌گیری از الگوریتم بهینه سازی ازدحام ذرات که ترکیبی از روش‌های بهینه‌سازی هیوریستیک و قطعی هستند نیز معرفی می‌شوند.

علاوه بر پرداختن به مباحث بیان شده، در مطلب «الگوریتم بهینه سازی ازدحام ذرات | کد الگوریتم PSO در پایتون ، متلب و جاوا | راهنمای جامع» چالش‌های اساسی که کاربر ضمن استفاده از الگوریتم PSO با آن‌ها مواجه می‌شود نیز مورد بررسی قرار گرفته‌اند. یک بررسی موردی (Case Study) نیز با بهره‌گیری از الگوریتم بهینه سازی ازدحام ذرات یا همان الگوریتم PSO انجام شده است که به درک بهتر مبحث کمک می‌کند. این مثال پیرامون بهینه‌سازی تابع هزینه برای سیستم تولید مثل با استفاده از الگوریتم PSO در بهینه‌سازی ترکیبی است. در نهایت، پیاده‌سازی الگوریتم بهینه سازی ازدحام ذرات در پایتون، متلب و جاوا انجام شده است.

مقدمه‌ای بر بهینه‌سازی و الگوریتم‌های آن

«بیشینه» (Maximizing) کردن «سود» یا «کمینه» (Minimizing) کردن «زیان» (Loss) از جمله مسائل بسیار حائز اهمیت در زمینه‌های گوناگون از جمله حوزه‌های فنی و مهندسی است. در یک تعریف ساده و کوتاه می‌توان گفت که مسائلی که در آن‌ها هدف بیشینه یا کمینه کردن یک تابع است را «مسئله بهینه‌سازی» (Optimization Problem) می‌گویند. برای مطالعه بیشتر پیرامون بهینه‌سازی، مطالعه مطالب «بهینه سازی (ریاضیاتی) چیست؟ — راهنمای جامع» و «بهینه سازی چند هدفه چیست؟ — راهنمای جامع» پیشنهاد می‌شود.

با توسعه فناوری، تعداد و پیچیدگی مسائل بهینه‌سازی نیز در زمینه‌های علمی گوناگون افزایش پیدا کرده است. از متداول‌ترین مسائل موجود در حوزه‌های مهندسی که نیاز به استفاده از روش‌های بهینه‌سازی برای حل آن‌ها وجود دارد می‌توان به تبدیل و توزیع انرژی، لجستیک (Logistics | آمادگاری) و بارگذاری مجدد رآکتورهای هسته‌ای اشاره کرد. مسائل بهینه‌سازی در دیگر زمینه‌ها از جمله هندسه و اقتصاد نیز کاربرد دارند. از دیگر زمینه‌های اصلی کاربرد بهینه‌سازی می‌توان به «هوش مصنوعی» (Artificial Intelligence | AI) و یادگیری ماشین «Machine Learning» اشاره کرد.

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

برای آشنایی بیشتر با روش‌های انتخاب بهترین الگوریتم برای یک مسئله خاص، مطالعه مطلب «روش انتخاب الگوریتم داده کاوی — راهنمای کاربردی» پیشنهاد می‌شود. یکی از مهم‌ترین موضوعات در انتخاب روش مناسب برای یک مسئله بهینه‌سازی، آشنایی کارشناس با انواع الگوریتم‌ها است تا بتواند مناسب‌ترین الگوریتم برای یک مسئله بهینه‌سازی را انتخاب کند. در این مطلب، یکی از الگوریتم‌های محبوب بهینه‌سازی یعنی الگوریتم بهینه سازی ازدحام ذرات مورد بررسی قرار گرفته و پیاده‌سازی الگوریتم PSO در متلب ، پایتون و جاوا انجام شده است. برای آشنایی با دیگر الگوریتم‌های بهینه‌سازی، مطالعه مطالب زیر پیشنهاد می‌شود.

• رویکرد هوش ازدحامی با استفاده از کلونی زنبور عسل مصنوعی برای حل مسائل بهینه‌سازی
• حل مسائل خوشه‌بندی با استفاده از الگوریتم کلونی زنبور عسل مصنوعی
• الگوریتم بهینه سازی فاخته – از صفر تا صد
• الگوریتم کرم شب تاب — از صفر تا صد
• الگوریتم ژنتیک – از صفر تا صد
• گرادیان کاهشی (Gradient Descent) و پیاده سازی آن در پایتون — راهنمای کاربردی
• الگوریتم کلونی مورچگان — از صفر تا صد
• الگوریتم بهینه سازی کلونی مورچگان در جاوا — راهنمای کاربردی
• شبیه سازی تبرید (Simulated Annealing) – به زبان ساده
• بهینه سازی نسبت طلایی — از صفر تا صد (+ دانلود فیلم آموزش رایگان)
• مهم‌ترین الگوریتم‌های یادگیری ماشین (به همراه کدهای پایتون و R) — بخش یازدهم و پایانی: الگوریتم‌های ارتقای گرادیان

الگوریتم بهینه سازی ازدحام ذرات (PSO)

در اوایل سال ۱۹۹۰ میلادی، پژوهش‌های گوناگونی پیرامون رفتار اجتماعی گروه‌های حیوانات انجام شد. این پژوهش‌ها حاکی از آن بودند که برخی از حیوانات که به یک گروه خاص متعلق هستند، مانند پرندگان، ماهی‌ها و دیگر موارد، قادر به آن هستند که اطلاعات را در گروه‌های (دسته‌های | گله‌های) خودشان به اشتراک بگذارند و چنین قابلیتی به این حیوانات مزایای قابل توجهی برای بقا اعطا می‌کرد.

با الهام گرفتن از این مطالعات، «کندی» (Kennedy) و «ابِرهارت» (Eberhart) در سال ۱۹۹۵ الگوریتم بهینه سازی ازدحام ذرات (Particle Swarm Optimization | PSO) یا الگوریتم PSO را در یک مقاله معرفی کردند. الگوریتم بهینه سازی ازدحام ذرات یا الگوریتم PSO یک الگوریتم «فراابتکاری» (Metaheuristic) است که برای بهینه‌سازی توابع پیوسته غیر خطی مناسب محسوب می‌شود. نویسندگان مقاله مذکور، الگوریتم بهینه سازی ازدحام ذرات یا الگوریتم PSO را از مفهوم هوش ذرات (Swarm Intelligence) که معمولا در گروه‌های حیوانات مانند گله‌ها و دسته‌های حیوانات وجود دارد الهام گرفته و ساخته‌اند.

برای شفاف شدن هر چه بیشتر ساز و کار کلی الگوریتم بهینه سازی ازدحام ذرات و دیگر الگوریتم‌هایی که از رفتار گروهی حیوانات الهام گرفته شده‌اند، توضیحاتی پیرامون رفتار گروهی (گله‌ای) حیوانات ارائه می‌شود. این توضیحات می‌تواند به درک چگونگی ساخت الگوریتم بهینه سازی ازدحام ذرات (و دیگر الگوریتم‌های دارای رویکرد مشابه) برای حل مسائل پیچیده ریاضی کمک کند.

الگوریتم بهینه سازی ازدحام ذرات و رفتار گروهی حیوانات

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

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

پژوهش‌هایی که از سال ۱۹۹۰ پیرامون رفتار پرندگان انجام شد، حاکی از آن است که همه پرندگان یک ازدحام (گروه | دسته) که به دنبال نقطه خوبی برای فرود هستند، قادر به آن هستند که از بهترین نقطه برای فرود در هنگامی که آن نقطه توسط یکی از اعضای ازدحام پیدا شد، آگاه شوند. با استفاده از این آگاهی، هر یک از اعضای این ازدحام، تجربه دانش شخصی و ازدحامی خود را متوازن می‌کنند که با عنوان «دانش اجتماعی» (Social Knowledge) شناخته شده است.

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

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

کندی و اِبِرهارت، از رفتار جمعی پرندگان الهام گرفتند؛ رفتاری که مزایای بقای قابل توجهی را برای پرندگان در هنگام جستجو برای یک نقطه امن برای فرود تضمین می‌کرد. آن‌ها بر همین اساس، الگوریتمی را ارائه کردند که الگوریتم ازدحام ذرات (Particle Swarm Optimization) نامیده می‌شود. الگوریتم PSO می‌تواند رفتاری به مثابه آنچه برای دسته پرندگان گفته شد را تقلید کند.

الگوریتم بهینه سازی ازدحام ذرات کلاسیک

نسخه اولیه الگوریتم ازدحام ذرات یا الگوریتم PSO که با عنوان نسخه کلاسیک این الگوریتم نیز شناخته شده است، در سال ۱۹۹۵ ارائه شد. از آن زمان تاکنون، انواع دیگری از این الگوریتم به عنوان نسخه‌های دیگر الگوریتم کلاسیک ارائه شده‌اند که از جمله آن‌ها می‌توان به «کاهش خطی وزن اینرسی» (Linear-Decreasing Inertia Weight)، «وزن عامل انقباض» (The Constriction Factor Weight) و «اینرسی پویا» (Dynamic Inertia) در کنار مدل‌های ترکیبی یا حتی روش‌های بهینه‌سازی الهام گرفته شده از کوانتوم که روی الگوریتم PSO اعمال شده‌اند اشاره کرد.

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

هدف از مسائل بهینه‌سازی، تعیین متغیری است که با بردار X=[x₁x₂x₃…x_n‎]‎ نشان داده می‌شود و بسته به فرمول بهینه‌سازی ارائه شده توسط تابع f(X)‎، بیشینه یا کمینه می‌شود. بردار متغیر X، به عنوان یک بردار مثبت شناخته شده است. این بردار، یک مدل متغیر و بردار n بُعدی آن را نمایش می‌دهد که در آن، n نشانگر تعداد متغیرهایی است که ممکن است در مسئله تعیین شوند. n در مسئله پیدا کردن بهترین نقطه برای فرود دسته پرندگان، طول و عرض جغرافیایی است.

از سوی دیگر، تابع f(X)‎ تابع برازش (Fitness Function) یا تابع هدف (Objective Function) نامیده می‌شود و تابعی است که میزان خوب یا بد بودن یک موقعیت X را ارزیابی می‌کند. این تابع برای مسئله دسته پرندگان، میزان خوب بودن یک نقطه برای فرود است که پرنده پس از پیدا کردن یک نقطه به آن فکر می‌کند. چنین ارزیابی برای مسئله فرود گروه پرندگان، براساس معیارهای بقای گوناگون انجام می‌شود. اکنون، ازدحامی با P ذره در نظر گرفته می‌شود؛ یک بردار مکان $$X_t^i$$ و یک بردار سرعت $$V_t^i$$ در هر تکرار برای هر یک از i ذره‌ای این سرعت را ایجاد می‌کنند، به صورت زیر وجود دارد:

$$X_i^t = (x_{i1}x_{i2}x_{i3} ... x_{in})^{T}$$

$$V_i^t = (v_{i1}v_{i2}v_{i3} ... v_{in})^{T}$$

این بردارها بر اساس بُعد j مطابق با معادله‌ای که در ادامه آمده است، به روز رسانی می‌شوند:

$$V_{i j}^{t+1}=w V_{i j}^{t}+c_{1} r_{1}^{t}\left(p b e s t_{i j}-X_{i j}^{t}\right)+c_{2} r_{2}^{t}\left(g b e s t_{j}-X_{i j}^{t}\right)$$

و

$$X_{i j}^{t+1}=X_{i j}^{t}+V_{i j}^{t+1}$$

که در آن‌ها، داریم:

i = 1, 2, …, P و j = 1, 2, …, n.

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

سرعتی که اولین عبارت در معادله را به روز رسانی می‌کند، ضرب داخلی پارامتر w و سرعت پیشین ذره است. به همین دلیل است که حرکت پیشین ذره به حرکت کنونی نمایش داده می‌شود. از همین رو، برای مثال، اگر w = 1 بود، حرکت ذره به طور کامل به وسیله حرکت قبلی خودش تحت تاثیر قرار گرفته است؛ بنابراین، ذره ممکن است به حرکت خود در همان جهت ادامه دهد.

از سوی دیگر، اگر $$0≤w<1$$، این تاثیر کاهش پیدا می‌کند و این یعنی ذرات به منطقه دیگری در ناحیه جستجو می‌روند. بنابراین، با توجه به کاهش پارامتر وزن اینرسی، ازدحام (گروه | دسته) ممکن است نواحی بیشتری را در ناحیه جستجو مورد اکتشاف قرار دهد و این یعنی شانس پیدا کردن بهینه سراسری افزایش پیدا می‌کند. اگرچه، در حالاتی که از مقادیر w کم‌تر استفاده می‌شود نیز هزینه‌ای وجود دارد که شبیه‌سازی‌ها را زمان‌برتر خواهد کرد.

عبارت درک فردی که دومین عبارت در معادله یک است، به وسیله تفاضل بین بهترین موقعیت خود ذره، برای مثال $$pbestij_{ij}$$ و موقعیت کنونی آن $$X_{ij}^t$$ محاسبه می‌شود. شایان توجه است که ایده نهفته در پس این ایده آن است که هر چه فعالیت‌ها فاصله بیشتری از موقعیت $$pbestij_{ij}$$ بگیرند، تفاضل $$(pbestij_{ij}-X_{ij}^t)$$ باید افزایش پیدا کند. بنابراین، این عبارت افزایش پیدا کرده و ذره را به بهترین موقعیت آن جذب می‌کند. پارامتر $$c_{1}$$ که به صورت حاصل‌ضرب در این رابطه وجود دارد، یک ثابت مثبت و یک پارامتر شناخت فردی محسوب می‌شود و به اهمیت تجربیات پیشین خود ذره وزن می‌دهد.

دیگر پارامتری که ضرب عبارت دوم را شکل می‌دهد، عبارت $$r_{1}$$ است. $$r_{1}$$ یک پارامتر مقدار تصادفی با طیف [0,1] است. این پارامتر تصادفی، نقش مهمی را بازی می‌کند، زیرا از همگرایی پارامترها ممانعت و بهینه سراسری احتمالی را بیشینه می‌کند. در نهایت، سومین عبارت مربوط به یادگیری اجتماعی است. به دلیل وجود این پارامتر، همه ذرات در ازدحام قادر به آن هستند که اطلاعات پیرامون بهترین نقطه به دست آمده را صرف نظر از اینکه کدام ذره آن را پیدا کرده است، با یکدیگر به اشتراک بگذارند؛ برای مثال $$gbestij_{ij}$$. فرمت این عبارت نیز درست مانند دومین عبارت است که مربوط به یادگیری فردی می‌شود. بنابراین، تفاضل $$(gbestij_{ij}-X_{ij}^t)$$ مانند یک جاذبه برای ذرات برای بهترین نقطه تا هنگام پیدا شدن نقطه در تکرار t عمل می‌کند. به طور مشابه، $$c_{2}$$ پارامتر یادگیری اجتماعی و وزن آن، اهمیت یادگیری سراسری ذرات است. همچنین، $$r_{2}$$ نیز نقشی مشابه با $$r_{1}$$ دارد.

در ادامه، الگوریتم PSO ارائه شده است و افراد ممکن است متوجه منطق بهینه‌سازی موجود در جستجوهای آن برای کمینه‌ها شوند و همه بردارهای مکانی که توسط تابع f(X)‎ ارزیابی می‌شوند. تابع f(X)‎ با عنوان «تابع برازش» (Fitness Function) شناخته شده است. در تصاویر ۲ و ۳ نیز به روز رسانی‌هایی در سرعت ذرات و موقعیت آن در تکرار t با در نظر داشتن مسئله دوبُعدی با متغیرهای $$x_{1}$$ و $$x_{2}$$ انجام شده است.

1. مقداردهی اولیه
   1. برای هر i در جمعیت ازدحام با اندازه p:
      1. $$X_{i}$$ را به طور تصادفی مقداردهی اولیه کن.
      2. $$$$x_{i}$V_{i}$$را به طور تصادفی مقداردهی اولیه کن. تابع برازش$$f(X_{i})$$را ارزیابی کن.$$pbestij_{ij}$$را با یک کپی از$$X_{i}$$مقداردهی اولیه کن. gbest را با یک نسخه از$$X_{i}$$با بهترین برازش مقداردهی اولیه کن. مراحل را تا هنگامی که یک معیار توقف ارضا شود، تکرار کن: برای هر ذره i:$$V_i^t$$و$$X_i^t$$را مطابق با معادلات ۱ و ۲ مقداردهی اولیه کن. تابع برازش$$f(X_i^t)$$را ارزیابی کن.$$pbest_{i} \leftarrow X_i^t$$اگر f(pbest_{i})$$gbest \leftarrow X_i^t$$h'v اگر f(gbest) بردار سرعت در تکرار t به صورتی که به وسیله دو مولفه ترکیب شده با ارجاع به یک مسئله دوبُعدی است. این بردار مکانی در تکرار t به روز رسانی شده، در حالیکه به وسیله دو مولفه با ارجاع به مسئله دوبُعدی ترکیب شده است. <b>برای مشاهده مجموعه فیلم‌های الگوریتم‌های بهینه‌سازی هوشمند</b> <b>+ اینجا کلیک کنید</b>. ترکیب الگوریتم ازدحام ذرات با روش‌های قطعی به طور کلی، روش‌های بهینه‌سازی به دو دسته قطعی (Deterministic) و هیوریستیک (Heuristic) تقسیم می‌شوند. هدف روش‌های قطعی پیدا کردن یک فرایند تکرار شونده شامل گرادیان (شیب) است که پس از تعداد مشخصی از تکرارها، به کمینه تابع هدف همگرا می‌شود. روال تکرار شونده این نوع از روش را می‌توان به صورت زیر نوشت:$$x^{k+1} = x^{k}+a^{k}d^{k}$$x بردار متغیر، α اندازه گام، d جهت کاهش و k تعداد تکرار است. بهترین ویژگی که از هر روش گرادیان قطعی انتظار می‌رود، همگرایی آن به یک نقطه ثابت است که معمولا کمینه محلی محسوب می‌شود. روش‌های هیوریستیک برخلاف روش‌های قطعی، از تابع هدف گرادیان در جهت رو به پایین استفاده نمی‌کنند. هدف آن‌ها تقلید طبیعت به منظور پیدا کردن بیشینه یا کمینه تابع هدف با انتخاب کردن نقاطی که این تابع محاسبه خواهد کرد، به شیوه‌ای ظریف و سازمان یافته است.. روش‌های ترکیبی یا ترکیبیاتی، در واقع ترکیبی از روش‌های قطعی و هیوریستیک را به منظور استفاده از مزایای هر دو روش ارائه می‌کنند. روش‌های ترکیبی معمولا از یک روش هیوریستیک برای پیدا کردن محتمل‌ترین مناطقی که کمینه محلی در آن‌ها قرار دارد استفاده می‌کنند. هنگامی که این منطقه تعیین شد، الگوریتم دارای فرمول ترکیبی به یک روش قطعی مبدل می‌شود تا به نقطه کمینه، نزدیک‌تر و سریع‌تر شود. معمولا، متداول‌ترین رویکرد مورد استفاده برای این فرمول‌سازی، استفاده از روش هیوریستیک برای تولید کاندیداهای خوب برای یک راه حل بهینه و سپس، استفاده از بهترین نقطه پیدا شده به عنوان یک نقطه شروع برای روش‌های قطعی به منظور همگرایی به کمینه محلی است. مقالات گوناگونی طی این سال‌ها ارائه شده‌اند که حاکی از کارایی و اثربخشی فرمول‌های ترکیبی هستند. همچنین، طی یک دهه اخیر تعداد رو به رشدی مقاله پیرامون روش‌های ترکیبی برای بهینه‌سازی ارائه شده است. در این زمینه، الگوریتم PSO را می‌توان با روش‌های قطعی ترکیب کرد و بدین شکل، شانس پیدا کردن محتمل‌ترین بهینه سراسری را افزایش داد. در ادامه، سه روش قطعی که الگوریتم PSO با آن‌ها ترکیب شده است معرفی می‌شوند. این روش‌ها عبارتند از: «روش گرادیان مزدوج» (Conjugate Gradient Method)، «روش نیوتون» (Newton’s Method) و «روش شبه‌نیوتون» (Quasi-Newton Method | BFGS). فرمول‌های هر یک از این الگوریتم‌ها در ادامه مطلب «الگوریتم بهینه سازی ازدحام ذرات | پیاده سازی الگوریتم PSO در پایتون ، متلب و جاوا | راهنمای جامع» به طور خلاصه مورد بررسی قرار گرفته‌اند. الگوریتم گرادیان مزدوج روش گرادیان مزدوج نرخ همگرایی «روش شدیدترین کاهش» (Steepest Descent Method) را با انتخاب جهت‌های کاهشی فراهم می‌کند که ترکیب خطی از جهت گرادیان با جهت کاهش تکرار پیشین هستند. بنابراین، روابط آن‌ها به صورت زیر است.$$x^{k+1} = x^{k}+a^{k}d^{k}$$

         $$d^{k} = -\triangledown(x^{k})+\gamma^{k}d^{k-1}$$در روابط بالا،$$\gamma$$ضریب مزدوجی است که با تنظیم اندازه بردارها کار می‌کند. در نسخه «فلچر-ریوز» (Fletcher-Reeves)، ضریب مزدوجی (هم‌یوغی) بر اساس رابطه زیر به دست می‌آید.$$\gamma^{k}=\frac{\parallel-\triangledown(x^{k})\parallel^{2}}{\parallel-\triangledown(x^{k-1})\parallel^{2}}$$روش نیوتون در حالی که روش‌های شدیدترین کاهش و گرادیان مزدوج از اطلاعات مشتق مرتبه اول استفاده می‌کنند، روش نیوتون (Newton’s Method) از اطلاعات مشتق مرتبه دوم برای شتاب‌دهی به هم‌گرایی فرایند تکرار شونده استفاده می‌کند. الگوریتم مورد استفاده در این روش، در ادامه ارائه شده است.$$x^{k+1}=x^{k}+\alpha^{k}d^{k}$$

         $$d^{k}=-\mid{H(x)}\mid^{-1}\triangledown{U(x^{k})}$$در رابطه بالا، H(x)‎، «ماتریس هسین» (Hessian Matrix) تابع است. به طور کلی، این روش نیازمند تکرارهای کمی برای همگرا شدن است. اگرچه، این روش نیاز به یک ماتریس دارد که با اندازه مسئله، رشد کند. اگر تخمین با کمینه تفاوت زیادی داشته باشد، ماتریس هسین ممکن است به طور ضعیفی در شرایط صدق کند. علاوه بر آن، این مورد شامل کشف ماتریسی است که موجب می‌شود روش‌ها به لحاظ کامپیوتری حتی پرهزینه‌تر هم باشند. روش شبه نیوتنی (BFGS) BFGS گونه‌ای از روش شبه نیوتنی (Quasi-Newton) است. این روش به دنبال تخمین معکوس ماتریس هسین با استفاده از اطلاعات گرادیان تابع است. این تخمین به گونه‌ای است که شامل مشتق مرتبه دوم نمی‌شود. بدین ترتیب، این روش دارای نرخ همگرایی کمتری نسبت به روش نیوتون است؛ هرچند که به لحاظ محاسباتی، سریع‌تر از روش نیوتون عمل می‌کند. الگوریتم روش شبه نیوتونی در ادامه ارائه شده است.$$x^{k+1} = x^{k} + a^{k} d^{k}$$

         $$d^{k} = -H^{k} \triangledown{U} (x^{k})$$

         $$H^{k} = H^{k-1} + M^{k-1} + N^{k-1}$$

         $$
         M^{k-1}=\left[1+\frac{\left(Y^{k-1}\right)^{T} \cdot H^{k-1} \cdot Y^{k-1}}{\left(Y^{k-1}\right)^{T} \cdot d^{k-1}}\right] \frac{d^{k-1} \cdot\left(d^{k-1}\right)^{T}}{\left(d^{k-1}\right)^{T} \cdot Y^{k-1}}
         $$

         $$N^{k-1} = -\frac{d^{k-1} (Y^{k-1})^T H^{k-1} + H^{k-1} Y^{k-1} (d^{k-1})^T)}{(d^{k-1})^T}$$

         $$Y^{k-1} = \triangledown{U} (x^{k}) - \triangledown{U} (x^{k-1})$$کاربردهای الگوریتم PSO و چالش‌های آن الگوریتم PSO قابل اعمال بر انواع مسائل در زمینه‌های گوناگون علمی است. به عنوان مثال، از الگوریتم PSO در حوزه بهداشت و درمان به منظور تشخیص بیماری سرطان خون از طریق تصویربرداری میکروسکوپی استفاده شده است. در علم اقتصاد، الگوریتم PSO برای تست پورتفولیوی ریسک سرمایه‌گذاری محدود و نامحدود به منظور به دست آوردن پرتفولیو ریسک بهینه استفاده شده است. در حوزه مهندسی، کاربردهای الگوریتم ازدحام ذرات یا همان الگوریتم PSO بسیار گوناگون است. مسائل بهینه‌سازی از جمله PSO را می‌توان در ادبیات پژوهش‌های سیستم‌های انتقال حرارت و الگوریتم‌های پیشی‌بینی ضریب انتقال حرارت نیز پیدا کرد. در این زمینه از ترمودینامیک، می‌توان به مقالات بهینه‌سازی شامل سیستم‌های حرارت مانند «سیکل رانکین آلی-موتور دیزل» (Diesel Engine–Organic Rankine Cycle)، سیستم دیزل ترکیبی ORC/فتوولتاییک (Hybrid Diesel-ORC/Photovoltaic) و «نیروگاه سیکل ترکیبی خورشیدی» (Solar Combined Cycle Power Plants | ISCC) اشاره کرد. همچنین از الگوریتم PSO برای مسائل بهینه‌سازی جغرافیایی به منظور پیدا کردن بهترین پیکربندی سیستم استفاده می‌شود که به بهترین شکل محدودیت‌های طراحی را ارضا می‌کند. در این زمینه، می‌توان به مطالعاتی اشاره کرد که شامل «بهینه‌سازی نوری-هندسی» (Optical-Geometric Optimization) متمرکز کننده‌های تابش خورشیدی و بهینه‌سازی جغرافیای برای «محوطه‌های تابشی» (Radiative Enclosures) می شود که توزیع دمایی و جریان گرما را ارضا می‌کند. پس از ارائه انواع گوناگونی از الگوریتم PSO، همچون مواردی که به آن‌ها در بخش‌های پیشین اشاره شد، PSO قادر به سر و کار داشتن با طیف وسیعی از مسائل، از مسائل موجود با تعداد کمی هدف و متغیرهای پیوسته گرفته تا دیگر مسائل چند هدفه چالش برانگیز با تعداد زیادی متغیر گسسته و یا پیوسته است. علاوه بر توانایی بالقوه الگوریتم PSO کاربر باید آگاه باشد که الگوریتم بهینه سازی ازدحام ذرات (الگوریتم PSO) تنها زمانی به نتایج مطلوب دست پیدا می‌کند که فرد یک تابع هدف را پیاده‌ئازی کند که قادر به منعکس کردن همه اهداف به صورت یکباره است. استخراج چنین تابعی ممکن است کار چالش برانگیزی باشد که نیاز به درک خوبی از مسئله فیزیکی برای حل شدن و توانایی استخراج ایده‌ها در یک معادله ریاضی دارد. مسائلی که در بخش پیشین این پژوهش ارائه شده‌اند، مثال‌هایی را از تابع هدف فراهم می‌کنند که می‌توانند این نقش را ایفا کنند. چالش دیگر برای افرادی که از الگوریتم PSO ‌استفاده می‌کنند، چگونگی مدیریت مرزهای فضای جستجو است که ذرات در آن حرکت می‌کنند. بسیاری از استراتژی‌های متداولی که در حال حاضر بر اساس الگوریتم PSO نسخه کلاسیک ارائه شده‌اند در مقالات گوناگون مورد بررسی قرار گرفته‌اند. مثال از الگوریتم PSO: هزینه یک سیستم تولید هم‌زمان این مسئله، یک مسئله کمینه‌سازی تابع است که کل هزینه عملیات یک سیستم تولید هم‌زمان به نام CGAM را نشان می‌دهد. این سیستم بر اساس نام سازندگان آن یعنی A. Valero ،G. Tsatsaronis ،C. Frangopoulos و M. von Spakovsky نام‌گذاری شده است که تصمیم به استفاده از سیستم مشابه برای مقایسه پاسخ‌های مسئله بهینه‌سازی با روش‌های گوناگون داشتند. در تصویر زیر، سیستم CGAM قابل مشاهده است. سیستم CGAM یک سیستم تولید هم‌زمان شامل «کمپرسور هوا» (Air Compressor | AC)، «محفظه احتراق» (Combustion Chamber | CC)، «توربین گاز» (Gas Turbine | GT)، «پیش‌گرم‌کن هوا» (Air Preheater | APH) و یک «دیگ‌بخار بازیافت حرارت» (Heat Recovery Steam Generator | HRSG) hsj که شامل یک اکونومایزر برای پیش‌گرم کردن آب و تبخیر کننده می‌شود. هدف از این چرخه، تولید ۳۰ مگاوات الکتریسیته و ۱۱ کیلوگرم بخار اشباع شده در فشار ۲۰ بار است. توضیح اقتصادی این سیستم، در این مطلب کاملا مشابه با مقاله اصلی است و هزینه سوخت سالیانه و هزینه‌های سالیانه مرتبط با تحصیل و عملیات هر تجهیزاتی می‌شود. معادلات برای هر مولفه در ادامه آمده است. کمپرسور هوا:$$Z_{AC} = (\frac{C_{11} \dot{m}_a}{C_{12} - \eta_{AC}})(\frac{P_{2}} {P_{1}}) ln(\frac{P_{2}} {P_{1}})$$محفظه احتراق:$$Z_{cc} = (\frac{C_{21}\dot{m}_a} {{C}_{22}-\frac{P_4}{P_3}})$$توربین:$$Z_{GT} = (\frac{C_{31}\dot{m}_{g}}{C_{32}-\eta_{GT}})ln(\frac{P_4}{P_5} )[1+exp(C_{33}T_{4}-C_{34})]$$پیش‌گرم‌کن:$$Z_{APH} = C_{41}(\frac{\dot{m}(h_5-h_6)}{(U)(\triangledown{TLM})})^{0.6}$$دیگ بخار بازیافت حرارت:$$Z_{HRSG} = C_{51}((\frac{Q_{PH}}{(\triangledown{TLM})_{PH}})^{0.8}+(\frac{Q_{PH}}{(\triangledown{TLM})_{PH}})^{0.8})+C_{52}\dot{m}_{st}+C_{53}\dot{m}_g^{1.2}$$</p> <p>عبارت کلی برای نرخ هزینه مربوط به سرمایه‌گذاری (S/$) برای هر مولفه در معادله زیر داده شده است.</p> <p dir=

         معادله زیر نشان‌گر هزینه کلی نرخ عملیات است.

         $$" />F = c_{1}\dot{m}_1PCI + \dot{Z}_{AC}+\dot{Z}_{APH}+\dot{Z}_{GT}+\dot{Z}_{HRSG}$$به منظور انجام بهینه‌سازی روی معادله بالا، سه متغیر تصمیم پذیرفته شده در تعریف مسئله اصلی در نظر گرفته می‌شوند. این متغیرها عبارتند از: نرخ فشرده‌سازی ($$\frac{P_{2}}{P_{1}}$$) بازدهی هم‌آنتروپی کمپرسور$$(\eta_{CA})$$بازدهی هم‌آنتروپی توربین$$(\eta_{GT})$$دمای هوا در خروجی پیش‌گرم کن$$(T_{3})$$دمای گاز سوخت در ورودی توربین$$(T_{4})$$

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

      	هیوریستیک	قطعی
      ترکیب ۱	ازدحام ذرات	گرادیان هم‌مزدوج
      ترکیب ۲	ازدحام ذرات	شبه نیوتون
      ترکیب ۳	ازدحام ذرات	نیوتون

      برای حل معادله ترمودینامیکی این مسئله، شبیه‌سازی تخصصی فرایند IPSEpro®‎ نسخه ۶.۰ مورد استفاده قرار گرفته است. IPSEpro®‎ یک شبیه‌سازی فرایند است که برای مدل‌سازی و شبیه‌سازی سیستم‌های حرارتی مختلف از طریق معادلات ترمودینامیکی آن‌ها استفاده می‌شود. این برنامه به وسیله «سیم‌تک» (SimTech) توسعه پیدا کرده و دارای یک رابط کاربرپسند و همچنین، طیف وسیعی از مولفه‌ها است که به کاربر امکان مدل‌سازی و شبیه‌سازی نیروگاه‌های متداول، سیستم‌های تولید هم‌زمان، چرخه‌های خنک کننده، چرخه‌های ترکیبی و بسیاری از دیگر موارد را می‌دهد.

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

      

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

      

      نتایج بهینه‌سازی در نمودارهای زیر ارائه شده است.

      

      

      

      پیاده‌سازی الگوریتم بهینه سازی ازدحام ذرات در پایتون

      همانطور که پیش از این بیان شد، اساس الگوریتم ازدحام ذرات یا الگوریتم PSO بر این موضوع بنا شده است که با هم کار کردن برای رسیدن به هدف، موثرتر از کار نکردن به صورت تیمی است. در الگوریتم PSO این رفتار اجتماعی با به اشتراک‌گذاری اطلاعات بین یک مجموعه از راه حل‌های ممکن که هر یک «ذره» (Particle) نامیده می‌شوند، انجام شده است.

      
      فیلم آموزش الگوریتم PSO در متلب MATLAB در فرادرس
      کلیک کنید

      یک ذره (چنانکه پیش از این نیز اشاره شد) به صورت زیر تعریف می‌شود:

      • مکان: برای مثال، هنگامی که متغیر طراحی به دو (صفحه) محدود شده است، یک ذره بر اساس مختصات آن تعریف می‌شود (x,y).
      • سرعت: این مقداری است که برای حرکت دادن محل ذره  به سمت بهترین راه حل انجام می‌شود.

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

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

      

      در انیمیشن زیر ساز و کار الگوریتم بهینه سازی ازدحام ذرات یا الگوریتم PSO قابل مشاهده است.

      

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

      
      فیلم آموزش حل مسائل عدد صحیح با استفاده از الگوریتم PSO (رایگان) در فرادرس
      کلیک کنید

      در ادامه، شبه کد الگوریتم PSO سراسری آمده است.

      For each particle
          Randomly initialize its position and velocity
      END
      Do
        For each particle
            Calculate its fitness value (measure of success)
            If the fitness value is better than the best fitness value (pBest) in history
                  Set current value as the new pBest
        END
      Choose the particle with the best fitness value of all the particles as the gBest
        For each particle
            Calculate particle velocity
            Update particle position using the computed velocity
        End
      While maximum iterations or minimum error criteria is not attained

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

      
      فیلم آموزش الگوریتم PSO در متلب MATLAB در فرادرس
      کلیک کنید
      • ذرات به طور تصادفی و با استفاده از «توزیع یکنواخت» (Uniform Distribution) مقداردهی اولیه می‌شوند.
      • به روز رسانی‌ها با استفاده از نسخه اندکی ویرایش شده الگوریتم اصلی (ویرایش توسط نویسندگان مقاله اصلی انجام شده) انجام می‌شود. در ادامه، کد پایتون برای پیاده‌سازی الگوریتم PSO ساده ارائه شده است.

      1import numpy as np
      2
      3class PSO(object):
      4  """
      5    Class implementing PSO algorithm.
      6  """
      7  def __init__(self, func, init_pos, n_particles):
      8    """
      9      Initialize the key variables.
      10      Args:
      11        func (function): the fitness function to optimize.
      12        init_pos (array-like): the initial position to kick off the
      13                               optimization process.
      14        n_particles (int): the number of particles of the swarm.
      15    """
      16    self.func = func
      17    self.n_particles = n_particles
      18    self.init_pos = np.array(init_pos)
      19    self.particle_dim = len(init_pos)
      20    # Initialize particle positions using a uniform distribution
      21    self.particles_pos = np.random.uniform(size=(n_particles, self.particle_dim)) \
      22                        * self.init_pos
      23    # Initialize particle velocities using a uniform distribution
      24    self.velocities = np.random.uniform(size=(n_particles, self.particle_dim))
      25
      26    # Initialize the best positions
      27    self.g_best = init_pos
      28    self.p_best = self.particles_pos
      29
      30  def update_position(self, x, v):
      31    """
      32      Update particle position.
      33      Args:
      34        x (array-like): particle current position.
      35        v (array-like): particle current velocity.
      36      Returns:
      37        The updated position (array-like).
      38    """
      39    x = np.array(x)
      40    v = np.array(v)
      41    new_x = x + v
      42    return new_x
      43
      44  def update_velocity(self, x, v, p_best, g_best, c0=0.5, c1=1.5, w=0.75):
      45    """
      46      Update particle velocity.
      47      Args:
      48        x (array-like): particle current position.
      49        v (array-like): particle current velocity.
      50        p_best (array-like): the best position found so far for a particle.
      51        g_best (array-like): the best position regarding
      52                             all the particles found so far.
      53        c0 (float): the cognitive scaling constant.
      54        c1 (float): the social scaling constant.
      55        w (float): the inertia weight
      56      Returns:
      57        The updated velocity (array-like).
      58    """
      59    x = np.array(x)
      60    v = np.array(v)
      61    assert x.shape == v.shape, 'Position and velocity must have same shape'
      62    # a random number between 0 and 1.
      63    r = np.random.uniform()
      64    p_best = np.array(p_best)
      65    g_best = np.array(g_best)
      66
      67    new_v = w*v + c0 * r * (p_best - x) + c1 * r * (g_best - x)
      68    return new_v
      69
      70  def optimize(self, maxiter=200):
      71    """
      72      Run the PSO optimization process untill the stoping criteria is met.
      73      Case for minimization. The aim is to minimize the cost function.
      74      Args:
      75          maxiter (int): the maximum number of iterations before stopping
      76                         the optimization.
      77      Returns:
      78          The best solution found (array-like).
      79    """
      80    for _ in range(maxiter):
      81      for i in range(self.n_particles):
      82          x = self.particles_pos[i]
      83          v = self.velocities[i]
      84          p_best = self.p_best[i]
      85          self.velocities[i] = self.update_velocity(x, v, p_best, self.g_best)
      86          self.particles_pos[i] = self.update_position(x, v)
      87          # Update the best position for particle i
      88          if self.func(self.particles_pos[i]) < self.func(p_best):
      89              self.p_best[i] = self.particles_pos[i]
      90          # Update the best position overall
      91          if self.func(self.particles_pos[i]) < self.func(self.g_best):
      92              self.g_best = self.particles_pos[i]
      93    return self.g_best, self.func(self.g_best)
      94  
      95# Example of the sphere function
      96def sphere(x):
      97  """
      98    In 3D: f(x,y,z) = x² + y² + z²
      99  """
      100  return np.sum(np.square(x))
      101
      102if __name__ == '__main__':
      103  init_pos = [1,1,1]
      104  PSO_s = PSO(func=sphere, init_pos=init_pos, n_particles=50)
      105  res_s = PSO_s.optimize()
      106  print("Sphere function")
      107  print(f'x = {res_s[0]}') # x = [-0.00025538 -0.00137996  0.00248555]
      108  print(f'f = {res_s[1]}') # f = 8.14748063004205e-06

      مشاهده کامل کدها

      برخی از نکاتی که باید پیرامون پیاده سازی الگوریتم بهینه سازی ازدحام ذرات در پایتون به آن‌ها توجه داشت در ادامه بیان شده‌اند.

      1. نقطه اولیه اهمیت زیادی در کارایی الگوریتم بهینه‌سازی دارد. به عنوان یک قاعده سر انگشتی، هرچقدر که نقطه اولیه از بهینه دور باشد، همانقدر طول می‌کشد تا الگوریتم همگرا شود. (این مورد برای همه الگوریتم‌های بهینه‌سازی صحت دارد.)
      2. PSO یک الگوریتم تصادفی است و به روز رسانی‌ها با استفاده از یک فرایند تصادفی انجام می‌شوند. این امر موجب می‌شود که پیگیری راه حل برای اجرای چند فرایند بهینه‌سازی دشوار باشد. اگرچه، همیشه می‌توان از دانه‌دهی (Seed) برای اعداد تصادفی استفاده کرد.
      3. این الگوریتم همیشه راه حل بهینه سراسری را پیدا نمی‌کند، اما در پیدا کردن بهینه‌ای که بسیار نزدیک به بهینه سراسری است، عملکرد بسیار خوبی دارد. برای تابع sphere بهینه سراسری در (0, 0, 0) است. کد الگوریتم PSO در پایتون که در اینجا ارائه شده است، نقطه دیگری را پیدا می‌کند که خیلی بد نیست.
      4. بسته به تعداد ذرات، همگرایی ممکن است زمان بیشتری طول بکشد. به طور کلی، بهتر است که این کار بیش از ۵۰ تکرار طول نکشد.
      5. برای توابع دشوار، ممکن است نیاز به تکرارهای بیشتری پیش از پیدا کردن بهترین راه حل باشد. برای مثال، می‌توان بین ۵۰۰ تا ۱۰۰۰ تکرار داشت.
      
      فیلم مجموعه آموزش الگوریتم های بهینه سازی هوشمند – مقدماتی تا پیشرفته در فرادرس
      کلیک کنید

      کد الگوریتم بهینه‌سازی ازدحام ذرات در متلب

      در ادامه، یک پیاده‌سازی ساده از الگوریتم بهینه‌سازی ازدحام ذرات در متلب نیز ارائه شده است.

      
      فیلم آموزش الگوریتم PSO در متلب MATLAB در فرادرس
      کلیک کنید

      Main

      1% Written by Dr. Seyedali Mirjalili
      2% ************************************************************************************************************************************************* 
      3clear 
      4close all
      5clc
      6% Problem preparation 
      7problem.nVar = 2;
      8problem.ub = 50 * ones(1, 2);
      9problem.lb = -50 * ones(1, 2);
      10problem.fobj = @ObjectiveFunction;
      11% PSO parameters 
      12noP = 4;
      13maxIter = 500;
      14visFlag = 1; % set this to 0 if you do not want visualization
      15RunNo  = 30; 
      16BestSolutions_PSO = zeros(1 , RunNo);
      17 [ GBEST  , cgcurve ] = PSO( noP , maxIter, problem , visFlag ) ;
      18 
      19 disp('Best solution found')
      20 GBEST.X
      21 disp('Best objective value')
      22 GBEST.O

      مشاهده کامل کدها

      ObjectiveFunction

      1% Written by Dr. Seyedali Mirjalili
      2% ************************************************************************************************************************************************* 
      3
      4
      5function [ GBEST ,  cgCurve ] = PSO ( noP, maxIter,  problem, dataVis )
      6% Define the details of the objective function
      7nVar = problem.nVar;
      8ub = problem.ub;
      9lb = problem.lb;
      10fobj = problem.fobj;
      11% Extra variables for data visualization
      12average_objective = zeros(1, maxIter);
      13cgCurve = zeros(1, maxIter);
      14FirstP_D1 = zeros(1 , maxIter);
      15position_history = zeros(noP , maxIter , nVar );
      16% Define the PSO's paramters
      17wMax = 0.9;
      18wMin = 0.2;
      19c1 = 2;
      20c2 = 2;
      21vMax = (ub - lb) .* 0.2;
      22vMin  = -vMax;
      23% The PSO algorithm
      24% Initialize the particles
      25for k = 1 : noP
      26    Swarm.Particles(k).X = (ub-lb) .* rand(1,nVar) + lb;
      27    Swarm.Particles(k).V = zeros(1, nVar);
      28    Swarm.Particles(k).PBEST.X = zeros(1,nVar);
      29    Swarm.Particles(k).PBEST.O = inf;
      30    
      31    Swarm.GBEST.X = zeros(1,nVar);
      32    Swarm.GBEST.O = inf;
      33end
      34% Main loop
      35for t = 1 : maxIter
      36    
      37    % Calcualte the objective value
      38    for k = 1 : noP
      39        
      40        currentX = Swarm.Particles(k).X;
      41        position_history(k , t , : ) = currentX;
      42        
      43        
      44        Swarm.Particles(k).O = fobj(currentX);
      45        average_objective(t) =  average_objective(t)  + Swarm.Particles(k).O;
      46        
      47        % Update the PBEST
      48        if Swarm.Particles(k).O < Swarm.Particles(k).PBEST.O
      49            Swarm.Particles(k).PBEST.X = currentX;
      50            Swarm.Particles(k).PBEST.O = Swarm.Particles(k).O;
      51        end
      52        
      53        % Update the GBEST
      54        if Swarm.Particles(k).O < Swarm.GBEST.O
      55            Swarm.GBEST.X = currentX;
      56            Swarm.GBEST.O = Swarm.Particles(k).O;
      57        end
      58    end
      59    
      60    % Update the X and V vectors
      61    w = wMax - t .* ((wMax - wMin) / maxIter);
      62    
      63    FirstP_D1(t) = Swarm.Particles(1).X(1);
      64    
      65    for k = 1 : noP
      66        Swarm.Particles(k).V = w .* Swarm.Particles(k).V + c1 .* rand(1,nVar) .* (Swarm.Particles(k).PBEST.X - Swarm.Particles(k).X) ...
      67            + c2 .* rand(1,nVar) .* (Swarm.GBEST.X - Swarm.Particles(k).X);
      68        
      69        
      70        % Check velocities
      71        index1 = find(Swarm.Particles(k).V > vMax);
      72        index2 = find(Swarm.Particles(k).V < vMin);
      73        
      74        Swarm.Particles(k).V(index1) = vMax(index1);
      75        Swarm.Particles(k).V(index2) = vMin(index2);
      76        
      77        Swarm.Particles(k).X = Swarm.Particles(k).X + Swarm.Particles(k).V;
      78        
      79        % Check positions
      80        index1 = find(Swarm.Particles(k).X > ub);
      81        index2 = find(Swarm.Particles(k).X < lb);
      82        
      83        Swarm.Particles(k).X(index1) = ub(index1);
      84        Swarm.Particles(k).X(index2) = lb(index2);
      85        
      86    end
      87    
      88    if dataVis == 1
      89        outmsg = ['Iteration# ', num2str(t) , ' Swarm.GBEST.O = ' , num2str(Swarm.GBEST.O)];
      90        disp(outmsg);
      91    end
      92    
      93    cgCurve(t) = Swarm.GBEST.O;
      94    average_objective(t) = average_objective(t) / noP;
      95    
      96    fileName = ['Resluts after iteration # ' , num2str(t)];
      97    save( fileName)
      98end
      99GBEST = Swarm.GBEST;
      100if dataVis == 1
      101    iterations = 1: maxIter;
      102    
      103%% Draw the landscape 
      104    figure
      105    
      106    x = -50 : 1 : 50;
      107    y = -50 : 1 : 50;
      108    
      109    [x_new , y_new] = meshgrid(x,y);
      110    
      111    for k1 = 1: size(x_new, 1)
      112        for k2 = 1 : size(x_new , 2)
      113            X = [ x_new(k1,k2) , y_new(k1, k2) ];
      114            z(k1,k2) = ObjectiveFunction( X );
      115        end
      116    end
      117    
      118    subplot(1,5,1)
      119    surfc(x_new , y_new , z);
      120    title('Search landscape')
      121    xlabel('x_1')
      122    ylabel('x_2')
      123    zlabel('Objective value')
      124    shading interp
      125    camproj perspective
      126    box on
      127    set(gca,'FontName','Times')
      128    
      129    
      130%% Visualize the cgcurve    
      131    subplot(1,5,2);
      132    semilogy(iterations , cgCurve, 'r');
      133    title('Convergence curve')
      134    xlabel('Iteration#')
      135    ylabel('Weight')
      136    
      137%% Visualize the average objectives
      138    subplot(1,5,3)
      139    semilogy(iterations , average_objective , 'g')
      140    title('Average objecitves of all particles')
      141    xlabel('Iteration#')
      142    ylabel('Average objective')
      143    
      144 %% Visualize the fluctuations 
      145    subplot(1,5,4)   
      146    plot(iterations , FirstP_D1, 'k');
      147    title('First dimention in first Particle')
      148    xlabel('Iteration#')
      149    ylabel('Value of the first dimension')
      150    
      151 %% Visualize the search history
      152    subplot(1,5,5)
      153    hold on
      154    for p = 1 : noP
      155        for t = 1 : maxIter
      156            x = position_history(p, t , 1);
      157            y = position_history(p, t , 2);
      158            myColor = [0+t/maxIter  0 1-t/maxIter ];
      159            plot(x , y , '.' , 'color' , myColor );
      160        end
      161    end
      162    contour(x_new , y_new , z);
      163    plot(Swarm.GBEST.X(1) , Swarm.GBEST.X(2) , 'og');
      164    xlim([lb(1) , ub(1)])
      165    ylim([lb(2) , ub(2) ])
      166    title('search history')
      167    xlabel('x')
      168    ylabel('y')
      169    box on
      170    
      171    set(gcf , 'position' , [128         372        1634         259])
      172    
      173    
      174    
      175    
      176end

      مشاهده کامل کدها

      PSO

      1% Written by Dr. Seyedali Mirjalili
      2% ************************************************************************************************************************************************* 
      3
      4function [ GBEST ,  cgCurve ] = PSO ( noP, maxIter,  problem, dataVis )
      5% Define the details of the objective function
      6nVar = problem.nVar;
      7ub = problem.ub;
      8lb = problem.lb;
      9fobj = problem.fobj;
      10% Extra variables for data visualization
      11average_objective = zeros(1, maxIter);
      12cgCurve = zeros(1, maxIter);
      13FirstP_D1 = zeros(1 , maxIter);
      14position_history = zeros(noP , maxIter , nVar );
      15% Define the PSO's paramters
      16wMax = 0.9;
      17wMin = 0.2;
      18c1 = 2;
      19c2 = 2;
      20vMax = (ub - lb) .* 0.2;
      21vMin  = -vMax;
      22% The PSO algorithm
      23% Initialize the particles
      24for k = 1 : noP
      25    Swarm.Particles(k).X = (ub-lb) .* rand(1,nVar) + lb;
      26    Swarm.Particles(k).V = zeros(1, nVar);
      27    Swarm.Particles(k).PBEST.X = zeros(1,nVar);
      28    Swarm.Particles(k).PBEST.O = inf;
      29    
      30    Swarm.GBEST.X = zeros(1,nVar);
      31    Swarm.GBEST.O = inf;
      32end
      33% Main loop
      34for t = 1 : maxIter
      35    
      36    % Calcualte the objective value
      37    for k = 1 : noP
      38        
      39        currentX = Swarm.Particles(k).X;
      40        position_history(k , t , : ) = currentX;
      41        
      42        
      43        Swarm.Particles(k).O = fobj(currentX);
      44        average_objective(t) =  average_objective(t)  + Swarm.Particles(k).O;
      45        
      46        % Update the PBEST
      47        if Swarm.Particles(k).O < Swarm.Particles(k).PBEST.O
      48            Swarm.Particles(k).PBEST.X = currentX;
      49            Swarm.Particles(k).PBEST.O = Swarm.Particles(k).O;
      50        end
      51        
      52        % Update the GBEST
      53        if Swarm.Particles(k).O < Swarm.GBEST.O
      54            Swarm.GBEST.X = currentX;
      55            Swarm.GBEST.O = Swarm.Particles(k).O;
      56        end
      57    end
      58    
      59    % Update the X and V vectors
      60    w = wMax - t .* ((wMax - wMin) / maxIter);
      61    
      62    FirstP_D1(t) = Swarm.Particles(1).X(1);
      63    
      64    for k = 1 : noP
      65        Swarm.Particles(k).V = w .* Swarm.Particles(k).V + c1 .* rand(1,nVar) .* (Swarm.Particles(k).PBEST.X - Swarm.Particles(k).X) ...
      66            + c2 .* rand(1,nVar) .* (Swarm.GBEST.X - Swarm.Particles(k).X);
      67        
      68        
      69        % Check velocities
      70        index1 = find(Swarm.Particles(k).V > vMax);
      71        index2 = find(Swarm.Particles(k).V < vMin);
      72        
      73        Swarm.Particles(k).V(index1) = vMax(index1);
      74        Swarm.Particles(k).V(index2) = vMin(index2);
      75        
      76        Swarm.Particles(k).X = Swarm.Particles(k).X + Swarm.Particles(k).V;
      77        
      78        % Check positions
      79        index1 = find(Swarm.Particles(k).X > ub);
      80        index2 = find(Swarm.Particles(k).X < lb);
      81        
      82        Swarm.Particles(k).X(index1) = ub(index1);
      83        Swarm.Particles(k).X(index2) = lb(index2);
      84        
      85    end
      86    
      87    if dataVis == 1
      88        outmsg = ['Iteration# ', num2str(t) , ' Swarm.GBEST.O = ' , num2str(Swarm.GBEST.O)];
      89        disp(outmsg);
      90    end
      91    
      92    cgCurve(t) = Swarm.GBEST.O;
      93    average_objective(t) = average_objective(t) / noP;
      94    
      95    fileName = ['Resluts after iteration # ' , num2str(t)];
      96    save( fileName)
      97end
      98GBEST = Swarm.GBEST;
      99if dataVis == 1
      100    iterations = 1: maxIter;
      101    
      102%% Draw the landscape 
      103    figure
      104    
      105    x = -50 : 1 : 50;
      106    y = -50 : 1 : 50;
      107    
      108    [x_new , y_new] = meshgrid(x,y);
      109    
      110    for k1 = 1: size(x_new, 1)
      111        for k2 = 1 : size(x_new , 2)
      112            X = [ x_new(k1,k2) , y_new(k1, k2) ];
      113            z(k1,k2) = ObjectiveFunction( X );
      114        end
      115    end
      116    
      117    subplot(1,5,1)
      118    surfc(x_new , y_new , z);
      119    title('Search landscape')
      120    xlabel('x_1')
      121    ylabel('x_2')
      122    zlabel('Objective value')
      123    shading interp
      124    camproj perspective
      125    box on
      126    set(gca,'FontName','Times')
      127    
      128    
      129%% Visualize the cgcurve    
      130    subplot(1,5,2);
      131    semilogy(iterations , cgCurve, 'r');
      132    title('Convergence curve')
      133    xlabel('Iteration#')
      134    ylabel('Weight')
      135    
      136%% Visualize the average objectives
      137    subplot(1,5,3)
      138    semilogy(iterations , average_objective , 'g')
      139    title('Average objecitves of all particles')
      140    xlabel('Iteration#')
      141    ylabel('Average objective')
      142    
      143 %% Visualize the fluctuations 
      144    subplot(1,5,4)   
      145    plot(iterations , FirstP_D1, 'k');
      146    title('First dimention in first Particle')
      147    xlabel('Iteration#')
      148    ylabel('Value of the first dimension')
      149    
      150 %% Visualize the search history
      151    subplot(1,5,5)
      152    hold on
      153    for p = 1 : noP
      154        for t = 1 : maxIter
      155            x = position_history(p, t , 1);
      156            y = position_history(p, t , 2);
      157            myColor = [0+t/maxIter  0 1-t/maxIter ];
      158            plot(x , y , '.' , 'color' , myColor );
      159        end
      160    end
      161    contour(x_new , y_new , z);
      162    plot(Swarm.GBEST.X(1) , Swarm.GBEST.X(2) , 'og');
      163    xlim([lb(1) , ub(1)])
      164    ylim([lb(2) , ub(2) ])
      165    title('search history')
      166    xlabel('x')
      167    ylabel('y')
      168    box on
      169    
      170    set(gcf , 'position' , [128         372        1634         259])
      171    
      172    
      173    
      174    
      175end

      مشاهده کامل کدها
      
      فیلم مجموعه آموزش الگوریتم های بهینه سازی هوشمند – مقدماتی تا پیشرفته در فرادرس
      کلیک کنید

      کد الگوریتم بهینه سازی ازدحام ذرات در جاوا

      در ادامه، کد الگوریتم PSO در جاوا انجام شده است.

      1import java.util.Arrays;
      2import java.util.Objects;
      3import java.util.Random;
      4import java.util.function.Function;
      5 
      6public class App {
      7    static class Parameters {
      8        double omega;
      9        double phip;
      10        double phig;
      11 
      12        Parameters(double omega, double phip, double phig) {
      13            this.omega = omega;
      14            this.phip = phip;
      15            this.phig = phig;
      16        }
      17    }
      18 
      19    static class State {
      20        int iter;
      21        double[] gbpos;
      22        double gbval;
      23        double[] min;
      24        double[] max;
      25        Parameters parameters;
      26        double[][] pos;
      27        double[][] vel;
      28        double[][] bpos;
      29        double[] bval;
      30        int nParticles;
      31        int nDims;
      32 
      33        State(int iter, double[] gbpos, double gbval, double[] min, double[] max, Parameters parameters, double[][] pos, double[][] vel, double[][] bpos, double[] bval, int nParticles, int nDims) {
      34            this.iter = iter;
      35            this.gbpos = gbpos;
      36            this.gbval = gbval;
      37            this.min = min;
      38            this.max = max;
      39            this.parameters = parameters;
      40            this.pos = pos;
      41            this.vel = vel;
      42            this.bpos = bpos;
      43            this.bval = bval;
      44            this.nParticles = nParticles;
      45            this.nDims = nDims;
      46        }
      47 
      48        void report(String testfunc) {
      49            System.out.printf("Test Function        : %s\n", testfunc);
      50            System.out.printf("Iterations           : %d\n", iter);
      51            System.out.printf("Global Best Position : %s\n", Arrays.toString(gbpos));
      52            System.out.printf("Global Best value    : %.15f\n", gbval);
      53        }
      54    }
      55 
      56    private static State psoInit(double[] min, double[] max, Parameters parameters, int nParticles) {
      57        int nDims = min.length;
      58        double[][] pos = new double[nParticles][];
      59        for (int i = 0; i < nParticles; ++i) {
      60            pos[i] = min.clone();
      61        }
      62        double[][] vel = new double[nParticles][nDims];
      63        double[][] bpos = new double[nParticles][];
      64        for (int i = 0; i < nParticles; ++i) {
      65            bpos[i] = min.clone();
      66        }
      67        double[] bval = new double[nParticles];
      68        for (int i = 0; i < bval.length; ++i) {
      69            bval[i] = Double.POSITIVE_INFINITY;
      70        }
      71        int iter = 0;
      72        double[] gbpos = new double[nDims];
      73        for (int i = 0; i < gbpos.length; ++i) {
      74            gbpos[i] = Double.POSITIVE_INFINITY;
      75        }
      76        double gbval = Double.POSITIVE_INFINITY;
      77        return new State(iter, gbpos, gbval, min, max, parameters, pos, vel, bpos, bval, nParticles, nDims);
      78    }
      79 
      80    private static Random r = new Random();
      81 
      82    private static State pso(Function<double[], Double> fn, State y) {
      83        Parameters p = y.parameters;
      84        double[] v = new double[y.nParticles];
      85        double[][] bpos = new double[y.nParticles][];
      86        for (int i = 0; i < y.nParticles; ++i) {
      87            bpos[i] = y.min.clone();
      88        }
      89        double[] bval = new double[y.nParticles];
      90        double[] gbpos = new double[y.nDims];
      91        double gbval = Double.POSITIVE_INFINITY;
      92        for (int j = 0; j < y.nParticles; ++j) {
      93            // evaluate
      94            v[j] = fn.apply(y.pos[j]);
      95            // update
      96            if (v[j] < y.bval[j]) {
      97                bpos[j] = y.pos[j];
      98                bval[j] = v[j];
      99            } else {
      100                bpos[j] = y.bpos[j];
      101                bval[j] = y.bval[j];
      102            }
      103            if (bval[j] < gbval) {
      104                gbval = bval[j];
      105                gbpos = bpos[j];
      106            }
      107        }
      108        double rg = r.nextDouble();
      109        double[][] pos = new double[y.nParticles][y.nDims];
      110        double[][] vel = new double[y.nParticles][y.nDims];
      111        for (int j = 0; j < y.nParticles; ++j) {
      112            // migrate
      113            double rp = r.nextDouble();
      114            boolean ok = true;
      115            Arrays.fill(vel[j], 0.0);
      116            Arrays.fill(pos[j], 0.0);
      117            for (int k = 0; k < y.nDims; ++k) {
      118                vel[j][k] = p.omega * y.vel[j][k] +
      119                    p.phip * rp * (bpos[j][k] - y.pos[j][k]) +
      120                    p.phig * rg * (gbpos[k] - y.pos[j][k]);
      121                pos[j][k] = y.pos[j][k] + vel[j][k];
      122                ok = ok && y.min[k] < pos[j][k] && y.max[k] > pos[j][k];
      123            }
      124            if (!ok) {
      125                for (int k = 0; k < y.nDims; ++k) {
      126                    pos[j][k] = y.min[k] + (y.max[k] - y.min[k]) * r.nextDouble();
      127                }
      128            }
      129        }
      130        int iter = 1 + y.iter;
      131        return new State(
      132            iter, gbpos, gbval, y.min, y.max, y.parameters,
      133            pos, vel, bpos, bval, y.nParticles, y.nDims
      134        );
      135    }
      136 
      137    private static State iterate(Function<double[], Double> fn, int n, State y) {
      138        State r = y;
      139        if (n == Integer.MAX_VALUE) {
      140            State old = y;
      141            while (true) {
      142                r = pso(fn, r);
      143                if (Objects.equals(r, old)) break;
      144                old = r;
      145            }
      146        } else {
      147            for (int i = 0; i < n; ++i) {
      148                r = pso(fn, r);
      149            }
      150        }
      151        return r;
      152    }
      153 
      154    private static double mccormick(double[] x) {
      155        double a = x[0];
      156        double b = x[1];
      157        return Math.sin(a + b) + (a - b) * (a - b) + 1.0 + 2.5 * b - 1.5 * a;
      158    }
      159 
      160    private static double michalewicz(double[] x) {
      161        int m = 10;
      162        int d = x.length;
      163        double sum = 0.0;
      164        for (int i = 1; i < d; ++i) {
      165            double j = x[i - 1];
      166            double k = Math.sin(i * j * j / Math.PI);
      167            sum += Math.sin(j) * Math.pow(k, 2.0 * m);
      168        }
      169        return -sum;
      170    }
      171 
      172    public static void main(String[] args) {
      173        State state = psoInit(
      174            new double[]{-1.5, -3.0},
      175            new double[]{4.0, 4.0},
      176            new Parameters(0.0, 0.6, 0.3),
      177            100
      178        );
      179        state = iterate(App::mccormick, 40, state);
      180        state.report("McCormick");
      181        System.out.printf("f(-.54719, -1.54719) : %.15f\n", mccormick(new double[]{-.54719, -1.54719}));
      182        System.out.println();
      183 
      184        state = psoInit(
      185            new double[]{0.0, 0.0},
      186            new double[]{Math.PI, Math.PI},
      187            new Parameters(0.3, 3.0, 0.3),
      188            1000
      189        );
      190        state = iterate(App::michalewicz, 30, state);
      191        state.report("Michalewicz (2D)");
      192        System.out.printf("f(2.20, 1.57)        : %.15f\n", michalewicz(new double[]{2.20, 1.57}));
      193    }
      194}

      مشاهده کامل کدها

      نتیجه‌گیری

      در مطلب «الگوریتم بهینه سازی ازدحام ذرات | کد الگوریتم PSO در پایتون ، متلب و جاوا | راهنمای جامع»، مباحث پایه‌ای روش PSO بیان، الگوریتم PSO تشریح و مزایا و معایب این روش مورد بررسی قرار گرفت. همچنین، روش‌های ترکیبی با الگوریتم PSO که در واقع ترکیبی از روش‌های قطعی و هیوریستیک هستند نیز به منظور بررسی مزایا و معایب آن‌ها مورد بررسی قرار گرفتند. سپس، مثال‌هایی پیرامون الگوریتم PSO ارائه و پیاده‌سازی الگوریتم ازدحام ذرات در زبان‌های برنامه‌نویسی پایتون، متلب و جاوا انجام شد.

      
      فیلم مجموعه آموزش الگوریتم های بهینه سازی هوشمند – مقدماتی تا پیشرفته در فرادرس
      کلیک کنید

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

      یکی از این استراتژی‌ها، مقایسه نتایج به دست آمده به وسیله دیگر الگوریتم‌های بهینه‌سازی است. در غیاب داده‌های بهینه در دسترس، به دلیل محدودیت‌های محاسباتی و یا حتی فقدان نتایج یک موضوع، این امکان وجود دارد که به عنوان استراتژی، مقایسه اطلاعات از مدل‌های فیزیکی واقعی انجام شود که از طریق الگوریتم‌های بهینه‌سازی قابل حصول نیستند. ولیکن منجر به فعالیت‌های مهندسی خوب و مقایسه نتایج از طریق تجربیات فنی می‌شوند. همچنین، این امکان وجود دارد که الگوریتم PSO روی مسائل مهندسی گوناگون اعمال شود که در مطلب «الگوریتم بهینه سازی ازدحام ذرات | کد الگوریتم PSO در پایتون ، متلب و جاوا | راهنمای جامع» نمونه‌هایی از آن مورد بررسی قرار گرفت.