{"id":2800,"date":"2024-07-20T03:16:22","date_gmt":"2024-07-20T10:16:22","guid":{"rendered":"https:\/\/gantovnik.com\/bio-tips\/?p=2800"},"modified":"2024-07-20T05:07:50","modified_gmt":"2024-07-20T12:07:50","slug":"437-particle-swarm-optimization-using-python","status":"publish","type":"post","link":"https:\/\/gantovnik.com\/bio-tips\/2024\/07\/437-particle-swarm-optimization-using-python\/","title":{"rendered":"#437 Particle swarm optimization using python"},"content":{"rendered":"

\"\"<\/a><\/p>\n

\nimport numpy as np\nimport matplotlib.pyplot as plt\nfrom matplotlib.animation import FuncAnimation\n\ndef f(x,y):\n    return (x-3.14)**2 + (y-2.72)**2 + np.sin(3*x+1.41) + np.sin(4*y-1.73)\n    \n# Compute and plot the function in 3D within [0,5]x[0,5]\nx, y = np.array(np.meshgrid(np.linspace(0,5,100), np.linspace(0,5,100)))\nz = f(x, y)\n\n# Find the global minimum\nx_min = x.ravel()[z.argmin()]\ny_min = y.ravel()[z.argmin()]\n\n# Hyper-parameter of the algorithm\nc1 = c2 = 0.1\nw = 0.8\n\n# Create particles\nn_particles = 10\nnp.random.seed(100)\nX = np.random.rand(2, n_particles)*5\nV = np.random.randn(2, n_particles)*0.1\n\n# Initialize data\npbest = X\npbest_obj = f(X[0], X[1])\ngbest = pbest[:, pbest_obj.argmin()]\ngbest_obj = pbest_obj.min()\n\ndef update():\n    # Function to do one iteration of particle swarm optimization\n    global V, X, pbest, pbest_obj, gbest, gbest_obj\n    r1, r2 = np.random.rand(2)\n    V = w*V + c1*r1*(pbest - X) + c2*r2*(gbest.reshape(-1,1)-X)\n    X = X + V\n    obj = f(X[0], X[1])\n    pbest[:, (pbest_obj >= obj)] = X[:, (pbest_obj >= obj)]\n    pbest_obj = np.array([pbest_obj, obj]).min(axis=0)\n    gbest = pbest[:, pbest_obj.argmin()]\n    gbest_obj = pbest_obj.min()\n\n# Set up base figure: The contour map\nfig, ax = plt.subplots(figsize=(8,6))\nfig.set_tight_layout(True)\nimg = ax.imshow(z, extent=[0,5,0,5], origin='lower', cmap='viridis', alpha=0.5)\nfig.colorbar(img, ax=ax)\nax.plot([x_min], [y_min], marker='x', markersize=5, color=\"white\")\ncontours = ax.contour(x, y, z, 10, colors='black', alpha=0.4)\nax.clabel(contours, inline=True, fontsize=8, fmt=\"%.0f\")\npbest_plot = ax.scatter(pbest[0], pbest[1], marker='.', color='red', alpha=0.5)\np_plot = ax.scatter(X[0], X[1], marker='.', color='blue', alpha=0.5)\np_arrow = ax.quiver(X[0], X[1], V[0], V[1], color='blue', width=0.002,\n                    angles='xy', scale_units='xy', scale=1)\ngbest_plot = plt.scatter([gbest[0]], [gbest[1]], marker='*', s=100, color='black', alpha=0.4)\nax.set_xlim([0,5])\nax.set_ylim([0,5])\n\ndef animate(i):\n    # Steps of PSO: algorithm update and show in plot\n    title = 'Iteration {:02d}'.format(i)\n    update()\n    ax.set_title(title)\n    pbest_plot.set_offsets(pbest.T)\n    p_plot.set_offsets(X.T)\n    p_arrow.set_offsets(X.T)\n    p_arrow.set_UVC(V[0], V[1])\n    gbest_plot.set_offsets(gbest.reshape(1,-1))\n    return ax, pbest_plot, p_plot, p_arrow, gbest_plot\n\nanim = FuncAnimation(fig, animate, frames=list(range(1,50)),\n                     interval=3000, blit=False, repeat=True)\nanim.save(\"ex437.gif\", dpi=120, writer=\"imagemagick\")\n\nprint(\"PSO found best solution at f({})={}\".format(gbest, gbest_obj))\nprint(\"Global optimal at f({})={}\".format([x_min,y_min], f(x_min,y_min)))\n<\/pre>\n
Output:<\/p>\n
PSO found best solution at f([3.18522118 3.13002418])=-1.8083515794181544\nGlobal optimal at f([3.1818181818181817, 3.131313131313131])=-1.8082706615747688<\/pre>\n","protected":false},"excerpt":{"rendered":"
import numpy as np import matplotlib.pyplot as plt from matplotlib.animation import FuncAnimation def f(x,y): return (x-3.14)**2 + (y-2.72)**2 + np.sin(3*x+1.41) + np.sin(4*y-1.73) # Compute and plot the function in 3D within [0,5]x[0,5] x, y = np.array(np.meshgrid(np.linspace(0,5,100), np.linspace(0,5,100))) z = f(x, y) # Find the global minimum x_min = x.ravel()[z.argmin()] y_min = y.ravel()[z.argmin()] # Hyper-parameter of […]<\/p>\n","protected":false},"author":1,"featured_media":0,"comment_status":"open","ping_status":"open","sticky":false,"template":"","format":"standard","meta":{"nf_dc_page":"","_et_pb_use_builder":"","_et_pb_old_content":"","_et_gb_content_width":"","_lmt_disableupdate":"yes","_lmt_disable":"","jetpack_post_was_ever_published":false,"_jetpack_newsletter_access":"","_jetpack_dont_email_post_to_subs":false,"_jetpack_newsletter_tier_id":0,"_jetpack_memberships_contains_paywalled_content":false,"_jetpack_memberships_contains_paid_content":false,"footnotes":""},"categories":[89,9,88,2],"tags":[105,129,128,3],"class_list":["post-2800","post","type-post","status-publish","format-standard","hentry","category-animation","category-optimization","category-pso","category-python","tag-animation","tag-imagemagick","tag-pso","tag-python"],"modified_by":"gantovnik","jetpack_featured_media_url":"","jetpack_sharing_enabled":true,"jetpack_shortlink":"https:\/\/wp.me\/p8bH0k-Ja","jetpack_likes_enabled":true,"jetpack-related-posts":[{"id":2059,"url":"https:\/\/gantovnik.com\/bio-tips\/2024\/01\/407-multiple-axes-animation-using-python\/","url_meta":{"origin":2800,"position":0},"title":"#407 Multiple axes animation using python","author":"gantovnik","date":"2024-01-14","format":false,"excerpt":"y=sin(x)","rel":"","context":"In "animation"","block_context":{"text":"animation","link":"https:\/\/gantovnik.com\/bio-tips\/category\/animation\/"},"img":{"alt_text":"","src":"https:\/\/i0.wp.com\/gantovnik.com\/bio-tips\/wp-content\/uploads\/2024\/01\/ex407.gif?resize=350%2C200&ssl=1","width":350,"height":200,"srcset":"https:\/\/i0.wp.com\/gantovnik.com\/bio-tips\/wp-content\/uploads\/2024\/01\/ex407.gif?resize=350%2C200&ssl=1 1x, https:\/\/i0.wp.com\/gantovnik.com\/bio-tips\/wp-content\/uploads\/2024\/01\/ex407.gif?resize=525%2C300&ssl=1 1.5x"},"classes":[]},{"id":3022,"url":"https:\/\/gantovnik.com\/bio-tips\/2024\/07\/442-the-partial-sums-of-the-fourier-series-for-the-square-wave-function-using-python\/","url_meta":{"origin":2800,"position":1},"title":"#442 The partial sums of the Fourier series for the square-wave function using python","author":"gantovnik","date":"2024-07-25","format":false,"excerpt":"import numpy as np import matplotlib.pyplot as plt from matplotlib.widgets import Slider nmax = 5 pi = np.pi x = np.linspace(-2*pi, 2*pi, 1001) def f(xarray): y = np.zeros_like(xarray) for ind, x in enumerate(xarray): xmod = x%(2*pi) if xmod