The algorithm is primarily inspired by three concepts in cosmology: white holes, black holes, and wormholes. The mathematical models of these three concepts are used for exploration, development, and local search, respectively. In 2015, it was published in the authoritative SCI journal Neural Computing and Applications. With nearly 2000+ citations, it still has a lot of potential. The multiverse optimization algorithm utilizes the concepts of white holes and wormholes to explore the search space, on the contrary, wormholes assist in the development of the search space in the multiverse optimization algorithm. Suppose that each solution is a universe and that the variables in each solution are objects in a universe. In the process of solving the optimization problem, the following rules are applied to the multiverse optimization algorithm: objects can move between different universes through white holes and black hole tunnels. When two universes are established between a white hole and a black hole tunnel, the universe has an expansion rate, and the universe with a higher expansion rate has a white hole, and the universe with a lower expansion rate is considered to have a black hole, and the particles in the universe are searched for a simulation by the principle of the transfer of particles from the white hole to the black hole position through the wormhole. This mechanism allows the universe to easily exchange objects with each other, regardless of the expansion rate of matter in the universe, all matter in the universe will randomly move through the wormhole to reach the optimal position of the universe. Among them, the number of iterations throughout the universe, can ensure an improvement in the average inflation rate. In each iteration, according to the classification of their inflation rate, the roulette mechanism is relied on to create a white hole. The composition of the algorithm individuals is as follows:
where d is the number of variables and n is the number of universes (candidate solutions).
step1.Initialization control parameters: number of universes n, parameters wep, tdr, threshold h and maximum number of iterations itermax;
step2.Initialize the cosmic group, find the current cosmic black hole and its location, and calculate its fitness value;
step3.According to the calculated value obtained in step 2, the optimal value and the optimal recording point are initialized.
step4.All individuals are transferred from a white hole to a black hole. Taking the current optimal cosmic location point as the center, the new optimal location point is generated according to the transfer aggregation rule, and the original point is replaced.
step5.The fitness of each point was calculated, and the optimal value and optimal cosmic point were recorded.
step6.Black hole individual update. Set the WEP, if the 2R WEP and the updated individual fitness value of the black hole is better than the original, then follow the updated formula (23) Update the individual, otherwise it will not be updated.
step7.Record the optimal values and the best values.
step8.Evaluate the entire cosmic community and analyze the results.
step9.If the number of iterations is less than the specified number of iterations, go to step4, otherwise go to step10.
step10.If the end prerequisites of the algorithm are met, the optimal solution is output and the end is completed, otherwise step 4 is continued.
In order to verify the results, the MVO algorithm is compared with the gray wolf optimization algorithm GWO, the particle swarm optimization algorithm PSO, the genetic algorithm GA and the gravitational search algorithm GSA. The results prove that the algorithm is able to provide very competitive results and outperforms the best algorithms in the literature on most test platforms. The results of practical case studies also demonstrate the potential of MVO in solving real-world problems with unknown search spaces. 
function[best_universe_inflation_rate,best_universe,convergence_curve]=mvo(n,max_time,lb,ub,dim,fobj)%two variables for s**ing the position and inflation rate (fitness) of the best universebest_universe=zeros(1,dim);best_universe_inflation_rate=inf;%initialize the positions of universesuniverses=initialization(n,dim,ub,lb);%minimum and maximum of wormhole existence probability (min and max in% eq.(3.3) in the **wep_max=1;wep_min=0.2;convergence_curve=zeros(1,max_time);%iteration(time) countertime=1;%main loopwhile time %eq. (3.3) in the **wep=wep_min+time*((wep_max-wep_min)/max_time); tr**elling distance rate (formula): eq. (3.4) in the **tdr=1-((time)^(1/6)/(max_time)^(1/6));inflation rates (i) (fitness values) inflation_rates=zeros(1,size(universes,1));for i=1:size(universes,1) %boundary checking (to bring back the universes inside search % space if they go beyoud the boundaries flag4ub=universes(i,:)ub; flag4lb=universes(i,:)universes(i,:)=(universes(i,:)flag4ub+flag4lb)))ub.*flag4ub+lb.*flag4lb; %calculate the inflation rate (fitness) of universes inflation_rates(1,i)=fobj(universes(i,:)elitism if inflation_rates(1,i) best_universe_inflation_rate=inflation_rates(1,i); best_universe=universes(i,:)end end [sorted_inflation_rates,sorted_indexes]=sort(inflation_rates); for newindex=1:n sorted_universes(newindex,:)=universes(sorted_indexes(newindex),:end %normaized inflation rates (ni in eq. (3.1) in the **normalized_sorted_inflation_rates=normr(sorted_inflation_rates); universes(1,:)= sorted_universes(1,:)update the position of universes for i=2:size(universes,1)%starting from 2 since the firt one is the elite back_hole_index=i; for j=1:size(universes,2) r1=rand();if r1 white_hole_index=roulettewheelselection(-sorted_inflation_rates);%for maximization problem -sorted_inflation_rates should be written as sorted_inflation_rates if white_hole_index==-1 white_hole_index=1; end %eq. (3.1) in the **universes(back_hole_index,j)=sorted_universes(white_hole_index,j); end if (size(lb,2)==1) %eq. (3.2) in the **if the boundaries are all the same r2=rand();if r2 r3=rand();if r3<0.5 universes(i,j)=best_universe(1,j)+tdr*((ub-lb)*rand+lb); end if r3>0.5 universes(i,j)=best_universe(1,j)-tdr*((ub-lb)*rand+lb); end end end if (size(lb,2)~=1) %eq. (3.2) in the **if the upper and lower bounds are %different for each variables r2=rand();if r2 r3=rand();if r3<0.5 universes(i,j)=best_universe(1,j)+tdr*((ub(j)-lb(j))*rand+lb(j));end if r3>0.5 universes(i,j)=best_universe(1,j)-tdr*((ub(j)-lb(j))*rand+lb(j));end end end end end %update the convergence curve convergence_curve(time)=best_universe_inflation_rate; %print the best universe details after every 50 iterations% if mod(time,50)==0% display(['at iteration ', num2str(time), ' the best universes fitness is ', num2str(best_universe_inflation_rate)])end time=time+1;endendfunction choice = roulettewheelselection(weights) accumulation = cumsum(weights); p = rand() accumulation(end); chosen_index = -1; for index = 1 : length(accumulation) if (accumulation(index) >p) chosen_index = index; break; end end choice = chosen_index;end