الگوریتم بهینه سازی ازدحام ذرات — کد الگوریتم 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=[x1x2x3…xn] نشان داده میشود و بسته به فرمول بهینهسازی ارائه شده توسط تابع 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}$$ انجام شده است.
- مقداردهی اولیه
- برای هر i در جمعیت ازدحام با اندازه p:
- $$X_{i}$$ را به طور تصادفی مقداردهی اولیه کن.
- $$$$x_{i}$V_{i}$$ را به طور تصادفی مقداردهی اولیه کن.
- تابع برازش $$f(X_{i})$$ را ارزیابی کن.
- $$pbestij_{ij}$$ را با یک کپی از $$X_{i}$$ مقداردهی اولیه کن.
- gbest را با یک نسخه از $$X_{i}$$ با بهترین برازش مقداردهی اولیه کن.
- برای هر i در جمعیت ازدحام با اندازه p:
- مراحل را تا هنگامی که یک معیار توقف ارضا شود، تکرار کن:
- برای هر ذره i:
- $$V_i^t$$ و $$X_i^t$$ را مطابق با معادلات ۱ و ۲ مقداردهی اولیه کن.
- تابع برازش $$f(X_i^t)$$ را ارزیابی کن.
- $$pbest_{i} \leftarrow X_i^t$$ اگر f(pbest_{i})<f(X_i^t)
- $$gbest \leftarrow X_i^t$$ h'v اگر f(gbest)<f(X_i^t)
- برای هر ذره i:
- برای مشاهده مجموعه فیلمهای الگوریتمهای بهینهسازی هوشمند + اینجا کلیک کنید.
ترکیب الگوریتم ازدحام ذرات با روشهای قطعی
به طور کلی، روشهای بهینهسازی به دو دسته قطعی (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}$$
عبارت کلی برای نرخ هزینه مربوط به سرمایهگذاری (S/$) برای هر مولفه در معادله زیر داده شده است.
$$\dot{Z}_{i,invest} = \frac{{Z}_{iφ}CRF}{N.3600}$$
CRF فاکتور بازیابی کلی (۱۸/۲ درصد)، N تعداد ساعات کاری نیروگاه (۸۰۰۰ ساعت)، و φ فاکتور نگهداری (۱/۰۶) است. علاوه بر آن، Cf هزینه سوخت به ازای واحد انرژی (۰.۰۰۴ MJ/$) است. جدول زیر نشانگر ثابتهای هزینه پذیرفته شده برای هر مولفه است.
معادله زیر نشانگر هزینه کلی نرخ عملیات است.
$$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) نامیده میشوند، انجام شده است.
یک ذره (چنانکه پیش از این نیز اشاره شد) به صورت زیر تعریف میشود:
- مکان: برای مثال، هنگامی که متغیر طراحی به دو (صفحه) محدود شده است، یک ذره بر اساس مختصات آن تعریف میشود (x,y).
- سرعت: این مقداری است که برای حرکت دادن محل ذره به سمت بهترین راه حل انجام میشود.
به اشتراکگذاری اطلاعات در الگوریتم ازدحام ذرات محلی و الگوریتم ازدحام ذرات سراسری به دو صورت مختلف انجام میشود که در ادامه بیان شدهاند.
- در 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
نقاطی که در اینجا باید مورد بررسی قرار بگیرند، مقداردهی اولیه ذرات و به روز رسانیهای آنها هم برای مکان و هم سرعت است.
- ذرات به طور تصادفی و با استفاده از «توزیع یکنواخت» (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
برخی از نکاتی که باید پیرامون پیاده سازی الگوریتم بهینه سازی ازدحام ذرات در پایتون به آنها توجه داشت در ادامه بیان شدهاند.
- نقطه اولیه اهمیت زیادی در کارایی الگوریتم بهینهسازی دارد. به عنوان یک قاعده سر انگشتی، هرچقدر که نقطه اولیه از بهینه دور باشد، همانقدر طول میکشد تا الگوریتم همگرا شود. (این مورد برای همه الگوریتمهای بهینهسازی صحت دارد.)
- PSO یک الگوریتم تصادفی است و به روز رسانیها با استفاده از یک فرایند تصادفی انجام میشوند. این امر موجب میشود که پیگیری راه حل برای اجرای چند فرایند بهینهسازی دشوار باشد. اگرچه، همیشه میتوان از دانهدهی (Seed) برای اعداد تصادفی استفاده کرد.
- این الگوریتم همیشه راه حل بهینه سراسری را پیدا نمیکند، اما در پیدا کردن بهینهای که بسیار نزدیک به بهینه سراسری است، عملکرد بسیار خوبی دارد. برای تابع sphere بهینه سراسری در (0, 0, 0) است. کد الگوریتم PSO در پایتون که در اینجا ارائه شده است، نقطه دیگری را پیدا میکند که خیلی بد نیست.
- بسته به تعداد ذرات، همگرایی ممکن است زمان بیشتری طول بکشد. به طور کلی، بهتر است که این کار بیش از ۵۰ تکرار طول نکشد.
- برای توابع دشوار، ممکن است نیاز به تکرارهای بیشتری پیش از پیدا کردن بهترین راه حل باشد. برای مثال، میتوان بین ۵۰۰ تا ۱۰۰۰ تکرار داشت.
کد الگوریتم بهینهسازی ازدحام ذرات در متلب
در ادامه، یک پیادهسازی ساده از الگوریتم بهینهسازی ازدحام ذرات در متلب نیز ارائه شده است.
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 در پایتون ، متلب و جاوا | راهنمای جامع» نمونههایی از آن مورد بررسی قرار گرفت.
با سلام و احترام خدمت خانم مهندس حصارکی
و تشکر از مطالب بسیار خوب و مفید شما
سپاسگزارم موفق باشید