A brief introduction to workshop scheduling

1. Definition of Shop Scheduling Shop scheduling refers to the allocation of processing shop sequence according to the reasonable demand of product manufacturing, so as to achieve the rational use of product manufacturing resources and improve the economic benefits of enterprises. The shop floor scheduling problem can be described mathematically as having n parts to be processed on M machines. The conditions to be met include that each process of each part shall use each machine no more than once and each part shall be processed in a certain order.

Traditional job shop schedulingTraditional job shop with scheduling instancesThere are a number of workpiece, each workpiece has a number of processes, there are multiple processing machines, but each process can only be processed on one machine. Corresponding to the example in the above table, there are two workpieces, and the workpiece J1 has three processes. The process Q11 can only be processed on M3, and the processing time is 5 hours. A constraint is that the relative order of operations cannot be changed for a workpiece. 11 – > O12 – > O13. Each workpiece can be processed on only one machine at a time; There can be only one work piece per machine. The task of scheduling is to arrange the processing order of the working procedure, and the processing order is determined, because there is only one machine available for each working procedure, and the processing machine is determined. The goal of scheduling is to minimize the total completion time (or any other goal). For example, after determining the processing order of O21->O22->O11->O23->O12->O13, we can calculate the total processing time according to the constraints of the processing machine. M2 processing O21 takes 6 hours, the current processing time of workpiece J2 is 6 hours. M1 processing O22 takes 9 hours, the current processing time of workpiece J2 is 6+9=15 hours. It takes 5 hours for M3 to process O11 and 5 hours for workpiece J1. M4 O23 processing cost 7 hours, workpiece J2 processing time 15+7=22 hours. M1 processing O12 takes 11 hours, but O12 can be processed only after M1 processes O22, so the current processing time of workpiece J1 is Max (5,9)+11=20 hours. M5 processing O13 consumes 8 hours, the workpiece J2 processing time 20+8=28 hours. The total completion time is Max (22,28)=28 hours.

Flexible job shop schedulingFlexible Job Shop with Scheduling Example (reference from Gao Liang’s paper “Improved Genetic Algorithm for Flexible Job Shop Scheduling Problem” — Journal of Mechanical Engineering)Compared with traditional job shop scheduling, flexible job shop scheduling relaxes the constraints on processing machines and is more in line with the reality of production. The optional processing machines in each process become multiple, which can be processed by one of multiple processing machines. For example, as shown in the above table, THE O12 process of J1 can be processed by M2 and M4, and the processing time is 8 hours and 4 hours respectively. However, M4 process is not always selected, and the final total completion time is shorter. Therefore, scheduling algorithm is needed to solve and optimize.

Compared with traditional job shop, the scheduling task of FLEXIBLE job shop scheduling not only needs to determine the processing sequence of each process, but also needs to determine the machine allocation of each process. For example, determine the processing order of O21->O22->O11->O23->O12->O13, we can not process the machine corresponding to the process, Therefore, the corresponding machine combination of [M1, M3, M5]->[M1, M2, M3]->[M1, M2, M3, M4, M5]->[M2, M3, M4, M5]->[M2, M4]->[M1, M3, M4, M5] should also be determined. The goal of scheduling is still the shortest total completion time (it can also be other goals, such as the shortest maximum machine load, the shortest total machine load).

Introduction to particle swarm optimization

The group behavior of birds and fish in nature has always been the research interest of scientists. Biologist Craig Reynolds proposed a highly influential model of bird clustering in 1987. In his simulation, each individual follows the following rules: avoid collisions with neighboring individuals; match the speed of neighboring individuals; Fly to the center of the flock, and the whole flock flies to the target. Using only the above three simple rules, we can simulate the phenomenon of bird flock flight very closely. In 1990, biologist Frank Heppner also proposed a bird model with a difference: birds are attracted to habitats. In the simulation, each bird has no specific flight target at the beginning, but only uses simple rules to determine its flight direction and speed. When one bird flies to the habitat, the birds around it will also fly to the habitat, and eventually the whole flock will fall into the habitat. In 1995, American social psychologist James Kennedy and electrical engineer Russell Berhart jointly proposed ParticleS Warm Optimization (PSO), The proposed algorithm is inspired by the research results of modeling and simulating the group behavior of birds. Their model and simulation algorithm mainly modified Frank Heppner’s model to make particles fly to the solution space and land at the optimal solution. Particle swarm optimization (PSO) has attracted extensive attention and become a research hotspot in the field of evolutionary computing due to its simplicity and easy implementation. Swarm Intelligence, co-authored by J.Kennedy and R.Eberhart in 2001, further expanded the influence of swarm intelligence [], followed by a large number of research reports and results on particle swarm optimization algorithm, which set off a research boom at home and abroad [2-7]. Particle swarm optimization (PSO) algorithm is derived from the regularity of bird swarm activity, and a simplified model is established by using swarm intelligence. It simulates the foraging behavior of birds, compares the search space for solving the problem to the flight space of birds, abstracts each bird into a particle without mass and volume, and uses it to represent a possible solution of the problem. The process of finding the optimal solution of the problem is regarded as the process of birds looking for food, and then solves the complex optimization problem. Particle swarm optimization (PSO), like other evolutionary algorithms, is also based on the concept of “population” and “evolution”, and achieves the search of optimal solution in complex space through cooperation and competition among individuals. At the same time, it is not like other evolutionary algorithms for individual evolution operator operation such as crossover, mutation, selection, but the individual in the group as there is no quality and volume in l d search space of particles, each particle moves in a certain speed in the solution space, and to its own historical best position P best and neighborhood historical best position g best gathered, Implement the evolution of candidate solutions. Particle swarm optimization (PSO) has a good biological social background and is easy to understand. Due to its few parameters and easy implementation, pSO has strong global search ability for nonlinear and multi-peak problems, and has been widely concerned in scientific research and engineering practice. At present, the algorithm has been widely used in function optimization, neural network training, pattern classification, fuzzy control and other fields.

2 Theory of particle swarm Algorithm 2.1 Particle swarm algorithm describes birds in the process of predation, members of a flock can obtain the discovery and flight experience of other members through information exchange and sharing between individuals. In a world where food sources are sporadic and unpredictable, the advantages of such collaboration are decisive, far outweighing the disadvantages of competition for food. Particle swarm optimization (pso) and inspired by the behavior of birds feed on to mimic this behavior, will the searching space of optimization problems in birds flying space, will every bird abstracted as a particle, particle no quality, no volume, used to characterize the a feasible solution of the problem, optimization problem to search for the optimal solution is equal to birds for food source. Particle swarm algorithm formulates simple behavior rules for each particle similar to bird movement, so that the movement of the whole particle swarm shows characteristics similar to bird predation, so as to solve complex optimization problems. The information sharing mechanism of particle swarm optimization algorithm can be interpreted as a symbiotic cooperative behavior, that is, each particle is constantly searching, and its searching behavior is affected by other individuals in the group to varying degrees [8]. Meanwhile, these particles also have the ability to remember the best position they have experienced. That is, their search behavior is not only influenced by other individuals but also guided by their own experience. Based on the unique search mechanism, pSO firstly generates the initial population, that is, randomly initializes the velocity and position of particles in the feasible solution space and velocity space, in which the position of particles is used to represent the feasible solution of the problem, and then solves the optimization problem through the cooperation and competition of individual particles among the population. Particle swarm optimization (PSO) algorithm is derived from the research on the predation behavior of birds: a group of birds randomly search for food in the area, and all birds know how far away they are from the food, so the simplest and most effective search strategy is to search the area around the nearest bird to the food. Particle swarm optimization (PSO) is applied to solve optimization problems. In particle swarm optimization, the potential solution to each optimization problem is a bird in the search space, called a particle. All particles have a fitness value determined by the optimized function, and each particle has a velocity that determines the direction and distance it flies. The particles then follow the current optimal particle to search the solution space [9].

Particle swarm optimization (PSO) firstly initializes the particle swarm randomly in a given solution space, and the dimension of the knowledge space is determined by the number of variables of the problem to be optimized. Each particle has its initial position and velocity, and then iteratively optimizes it. In each iteration, each particle by tracking two “extreme” to update their spatial location in the solution space and speed: a extreme value is a single particles themselves to find the optimal solution in the process of iterative particle, the particle called the individual extremum: another extreme value is population all particles in the process of iteration to find the optimal solution of the particle, the particle is the global extremum. The method described above is called global particle swarm optimization. If we do not use all the particles of the population but only use some of them as the neighbor particles of the particle, then the extreme value in all the neighbor particles is the local extreme value, which is called the local particle swarm optimization algorithm.

2.3 Characteristics of PARTICLE Swarm Optimization Algorithm Particle swarm optimization algorithm is a random search algorithm, it is a new intelligent optimization technology. The algorithm can converge to the global optimal solution with high probability. It is proved by practice that it is suitable for dynamic and multi-objective optimization environment, and has faster calculation speed and better global search ability compared with the traditional optimization algorithm. (1) Particle swarm optimization algorithm is an optimization algorithm based on swarm intelligence theory. It guides optimization search through swarm intelligence generated by cooperation and competition among particles in the population. Compared with other algorithms, particle swarm optimization is an efficient parallel search algorithm. (2) Both particle swarm optimization and genetic algorithm randomly initialize the population, using the adaptive value to evaluate the merits of the individual and carry out a certain random search. However, particle swarm optimization (PSO) decides the search according to its own speed without the crossover and variation of genetic algorithm. Compared with the evolutionary algorithm, particle swarm optimization retains the global search strategy based on population, but its speed-displacement model is simple to operate and avoids complex genetic operation. (3) Since each particle still maintains its individual extreme value at the end of the algorithm, that is, particle swarm optimization algorithm can not only find the optimal solution of the problem, but also obtain several better sub-optimal solutions, so the application of particle swarm optimization algorithm to scheduling and decision problems can give a variety of meaningful solutions. (4) The unique memory of particle swarm optimization enables it to dynamically track the current search situation and adjust its search strategy. In addition, particle swarm optimization is not sensitive to the size of the population.

3 types of particle swarm optimization 3.1 Elementary particle swarm optimization 3.2 Standard particle swarm optimizationTwo concepts commonly used in particle swarm optimization are introduced: first, “exploration”, which means that particles leave the original search track to a certain extent and search in a new direction, reflecting the ability to explore unknown areas, similar to global search; The second is “development”, which means that particles continue to search in a more detailed step on the original search track to a certain extent, mainly referring to the further search of the region searched in the process of exploration. Exploration is to deviate from the original search trajectory to find a better solution, and exploration ability is the global search ability of an algorithm. Development is to use a good solution, continue the original search path to search for a better solution, it is the local search ability of the algorithm. How to determine the proportion of local search ability and global search ability is very important for solving a problem. In 1998, Shi Yuhui et al. proposed an improved particle swarm optimization algorithm with inertial weight [10], which is regarded as the standard particle swarm optimization algorithm by default because it can guarantee good convergence effect. Its evolutionary process is as follows:In Equation (6.7), the first part represents the previous velocity of the particle to ensure the global convergence performance of the algorithm. The second part and the third part make the algorithm have local convergence ability. It can be seen that the inertia weight w in Equation (6.7) indicates to what extent the original speed should be retained: larger w means stronger global convergence and weaker local convergence. If W is small, the local convergence is strong, but the global convergence is weak. When w=1, Equation (6.7) is exactly the same as Equation (6.5), indicating that particle swarm optimization algorithm with inertial weight is an extension of elementary particle swarm optimization algorithm. The experimental results show that w is 0.8When between 1.2, particle swarm optimization algorithm has faster convergence speed; When w>1.2, the algorithm tends to fall into local extremum. In addition, w can be dynamically adjusted in the search process: at the beginning of the algorithm, w can be given a large positive value, and with the progress of the search, W can be linearly reduced, which can ensure that at the beginning of the algorithm, each particle can explore a better region in the global range with a large speed step. In the later stage of the search, a smaller value of W ensures that particles can conduct a fine search around the extreme point, so that the algorithm has a greater probability of convergence to the global optimal solution position. Adjusting W allows you to balance global and local search capabilities. At present, the linear decreasing weight strategy proposed by Shi adopts more dynamic inertia weight values, and its expression is as follows: 3.3 Compression factor particle swarm optimizationClerc et al. proposed the use of constraint factors to control the final convergence of system behavior [11]. This method can effectively search different regions and obtain high-quality solutions. The speed updating formula of compression factor method is:The experimental results show that compared with the particle swarm optimization algorithm with inertia weight, the particle swarm optimization algorithm with constraint factor has faster convergence speed.3.4 Discrete particle swarm optimizationThe basic particle swarm optimization (PSO) algorithm is a powerful tool for searching the extreme value of a function in continuous domain. After elementary particle swarm optimization, Kennedy and Eberhart proposed a discrete binary version of particle swarm optimization [12]. In this discrete particle swarm optimization method, the discrete problem space is mapped to the continuous particle motion space, and the particle swarm optimization algorithm is modified appropriately to solve the problem. The values and variations of particles in the state space are limited to 0 and 1, and each dimension of the velocity, vi and y, represents the probability that each xi of position is equal to 1. Therefore, the viJ updating formula in the continuous particle swarm remains unchanged, but P best and: best are only valued within [0, 1]. Its position update equation is expressed as follows: 4 process of particle swarm optimizationBased on the concepts of “population” and “evolution”, particle swarm optimization algorithm realizes the search for optimal solution in complex space through cooperation and competition among individuals [13]. Its process is as follows: (1) Initialization of particle swarm, including population size N and position X of each particle; And the velocity Vio (2) calculates the fitness value of each particle fit[I]. (3) For each particle, compare its fitness value fit[gate and individual extremum P best(I). If fit[I] <P best(I), replace P best(I) with FIT [I]. (4) For each particle, compare its fitness value fit[I] with the global extreme value G Best. If Fit [I] < 8 Best, replace G Best with Fit [I]. (5) Iteratively update the particle velocity V; and position Xj. (6) perform boundary condition processing. (7) Determine whether the termination conditions of the algorithm are met: if so, terminate the algorithm and output the optimization results; otherwise, return to Step (2). 5 Key ParametersIn particle swarm optimization, the choice of control parameters can affect the performance and efficiency of the algorithm. It is a complex optimization problem how to choose appropriate control parameters to optimize the performance of the algorithm. In practical optimization problems, the control parameters are usually selected according to the user’s experience. The control parameters of pSO mainly include particle population size N, inertia weight W, acceleration coefficient C and C, maximum velocity Via x, stop criterion, neighborhood structure setting, boundary condition processing strategy, etc. [14].Particle population size NThe selection of particle population size depends on the specific problem, but generally the number of particles is set to 2050. For most problems 10 particles is a good result: but for more difficult problems or for certain types of problems, the number of particles can go up to 100 or 200. In addition, the larger the number of particles is, the larger the search space of the algorithm is, and it is easier to find the global optimal solution. Of course, the longer the algorithm runs.Inertial weight WThe inertial weight W is a very important control parameter in the standard particle swarm optimization algorithm, which can be used to control the development and exploration ability of the algorithm. The magnitude of the inertial weight indicates how much is inherited from the particle’s current velocity. When the inertia weight value is large, the global optimization ability is strong, and the local optimization ability is weak: When the inertia weight value is small, the global optimization ability is weak, and the local optimization ability is strong. The choice of inertia weight usually includes fixed weight and time varying weight. Fixed weight means that constant is selected as the inertia weight value, and its value remains unchanged in the evolution process, with a general value of [0.8, 1.2]. Time-varying weight is to set a certain change interval and gradually reduce the inertia weight in a certain way in the evolution process. The selection of time – varying weights includes range of change and decline rate. A fixed inertial weight allows the particle to maintain the same exploration and exploitation capacity, while a time-varying weight allows the particle to have different exploration and exploitation capacity at different stages of evolution.The acceleration constants c1 and C2The acceleration constants C and C2 adjust the maximum stride length to P best and G best, respectively. They determine the influence of individual particle experience and group experience on particle trajectory, respectively, and reflect the information exchange between particle swarm. If Cr =c2=0, the particle will fly to the boundary at its current flight speed. At this point, the particle can only search a limited area, so it is difficult to find the optimal solution. If q=0, it is a “social” model, in which the particle lacks cognitive ability and only has group experience. Its convergence speed is fast, but it is easy to fall into local optimum. If OY =0, it is a “cognitive” mode. There is no shared information of the society and no information interaction between individuals, so the probability of finding the optimal solution is small. A group with a size of D is equivalent to running N independent particles. So c1=C2, usually c1= CG =1.5. In this way, individual experience and group experience have the same important influence, making the final optimal solution more accurate.The maximum velocity of the particle is vmaxThe particle speed in each dimension of space has a maximum speed limit value vd Max, used to control the particle speed, so that the speed is controlled within the range [-vimax, + VA Max], which determines the strength of the problem space search, the value is generally set by the user. Vmax is a very important parameter. If it is too large, particles may pass through the good region; if it is too small, particles may not be able to detect the region beyond the local optimal region adequately. They may get stuck in local optimality and not be able to move far enough to jump out of local optimality to a better position in space. The researchers point out that setting Vmax and adjusting inertia weights are equivalent, so! Max is generally used to set the initialization of the population, that is, vmax is set as the variation range of each dimension variable, instead of careful selection and adjustment of the maximum speed.Stopping criteriaThe maximum number of iterations, computational accuracy, or the maximum number of stagnant steps for the optimal solution ▲ T (or acceptable satisfactory solution) is usually considered to be the stop criterion, the termination condition of the algorithm. According to the specific optimization problem, the setting of the stop criterion should take into account the performance of the algorithm, such as the solution time, optimization quality and search efficiency.Setting up the neighborhood structureThe global pSO takes the whole population as the neighborhood of particles and has the advantage of fast convergence, but sometimes the algorithm will fall into local optimization. The local version of particle swarm optimization algorithm takes the individuals with similar positions as the neighborhood of particles, and the convergence speed is slow and it is not easy to fall into the local optimal value. In practical applications, the global particle swarm optimization algorithm can be used to find the direction of the optimal solution, that is, to get a rough result, and then the local particle swarm optimization algorithm can be used to conduct fine search near the optimal point.Boundary condition processingWhen the position or velocity of one or several dimensions exceeds the set value, the boundary condition processing strategy can limit the position of particles in the feasible search space, so as to avoid the expansion and divergence of the population, and avoid the blind search of particles in a large range, thus improving the search efficiency. There are many ways to do this, such as setting the maximum position limit Xmax and the maximum speed limit Vmax. When the maximum position or speed is exceeded, a numerical value is randomly generated in the range to replace it, or it is set to the maximum value, that is, boundary absorption.

Case and part of the source code

1 caseA garment enterprise decides to process 9 batches of its brand clothing. Due to seasonal changes, each batch of products has a delivery deadline. If the products are completed before this deadline, they can be sold at a higher price; after the deadline, they will face more fierce price competition and reduce production efficiency. It is assumed that the processing time, delivery time and profit of each batch of clothing are shown in Table 1, and the processing process of each batch of products is not allowed to be interrupted, that is, the processing of one batch of products cannot be inserted into the processing of other batches of products, so as to find the processing sequence with the maximum total profit. Table 1 Processing time and delivery deadline of the brand products of the enterprise 2 Part source code

clear all;
close all;
clc;
N=100; % Number of population particles T=20; % Maximum number of iterations c1=2; % Individual optimal learning factor c2=2; % global optimal learning factor Vmax=4; % Maximum speed Vmin=4 -; % Minimum speed Wmax=0.9; % Maximum inertia weight Wmin=0.4; % Minimum inertia weight RRR Channel=9; % Sets the number of idle channels Counter=0; % number of updates Flag=zeros(T,1); % overMatrix=zeros(T,N); % overdue time matrix minTbest=zeros(20.1); % overdue time matrix results=zeros(20.11); % Result matrixfor i=1:20
    minTbest(i,1) =1000;
end
step=4; D=step.*Channel; % rbest=ones(20,Channel);

% vacant_time=[10.60.25.170.83.5.54.155]; % Set the average idle time of each channel (calculate the mean) time=[3 4 1 2 6 1 4 7 5]; % time matrix deadline=[5 9 3 12 10 24 5 6 6];
    profilein=[750 1200 800 900 2500 500 3000 5600 4500];
    profileout=[500 900 400 750 1800 300 1500 4000 2000];
    %time=[0 3 6 5 1 3 1 5 6 3 2]; % time matrix %deadline=[2 4 5 5 5 15 15 20 24 35 35];
    %profilein=[1250 1600 1800 2500 500 1700 650 1500 1800 1200 800];
    %profileout=[1000 1280 1440 200 0 1190 450 1050 1260 660 500]; % % % % % % % % % % % % % % % initialization population individuals (limited speed and position) % % % % % % % % % % % % % % % % % % % % % % % x = randint (N, D); % Initial population v=rand(N,D)*(vmax-vmin)+Vmin; The initialization speed % % % % % % % % % % % % % % % % initialization individual optimal location and the optimal value % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % p = x; pbest=ones(N,1);
tbest=ones(N,1);
for i=1:N [pbest(i),tbest(i)]=fitness(x(i,:),step,Channel,time,deadline, profilein,profileout); End % % % % % % % % % % % % % % % % initializing global optimal location and the optimal value % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % g = ones (1,D);
gbest=inf;
for i=1:N
    if(pbest(i)>gbest)
        g=p(i,:);
        gbest=pbest(i);
    end
end
gb=ones(1,T); % % % % % % % % % % % % in accordance with the formula iteration in order to meet the precision and the optimal value % % % % % % % % % % % % % % % % % % % % %for i=1:T
    for j=1:N %%%%%%%%% Updates individual optimal location and optimal value %%%%%%%%%%%%%%% %if (fitness(x(j,:),step,Channel,time,deadline, profilein,profileout))
        p(j,:)=x(j,:);
        [pbest(j),tbest(j)]=fitness(x(j,:),step,Channel,time,deadline, profilein,profileout);
        overTime = overdueTime(time,deadline);
        totalProcessTime = 0; %end %%%%%%%%% Updates the global optimal location and value %%%%%%%%%%%%%%%%%%%if(pbest (j) < gbest) % % % % % % % % % % % update total overdue time % % % % % % % % % % % % % % %if tbest(j) < minTbest(i,1)
                minTbest(i,1) = tbest(j);
                g=p(j,:);
                gbest=pbest(j);
                Flag(i,1) =1; End end % % % % % % % % % % % % update overdue time global optimal location and the optimal value % % % % % % % % % % % % % % %if tbest(j) < minTbest(i,1)
        %    minTbest(i,1) = tbest(j); % g=p(j,:); % gbest=pbest(j); End % % % % % % % % % % % % to calculate dynamic inertia weight value % % % % % % % % % % % % % % % w = Wmax - (Wmax - Wmin) * I/T; % % % % % % % % % % % update position and velocity value % % % % % % % % % % % % % % % % % % v (j) = v (j:) + c1 * w * rand * (p (j:) - x (j:))... +c2*rand*(g-x(j,:)); %%%%%%%%%%%% boundary condition processing %%%%%%%%%%%%%%%%%%%%%%%for ii=1:D
            if(v(j,ii)>Vmax)|(v(j,ii)<Vmin)
                v(j,ii)=rand*(Vmax-Vmin)+Vmin;
            end
        end
        vx(j,:)=1./ (1+exp(-v(j,:)));
        for jj=1:D
            if vx(j,jj)>rand
                x(j,jj)=1;
            else
                x(j,jj)=0; End end end % % % % % % % % % % % % % % % % % % % % % % fitness function % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % function [result,totalProcessTime]=fitness(x,step,Channel,time,deadline, Profilein, profileOut) %singleDec: channel sequence represented by a single particle R %vacant_time: channel idle time %fitness: expectation as fitness singleDec=bin2decStep(x,step,Channel); singleDec=singleDec+1;
M=length(singleDec);
sortorder=singleDec;
totalProcessTime=0;
% for i=1:M
%     if(((i- 1)*T+Th)<permChannel(i))
%         pFail(i)=((i- 1)*T+Th)./permChannel(i);
%     else
%         pFail(i)=1;
%     end
% end

% for i=1:M
%     s1=s1*pFail(i);
% end
    time=time(sortorder);
    deadline=deadline(sortorder);
    profilein=profilein(sortorder);
    profileout=profileout(sortorder);
    processTime=accumulate(time);
    overTime=overdueTime(time,deadline);
    index1= deadline>=processTime;
    index2= deadline<processTime;
    %index3=find(overTime);
    profilein(index2)=0;
    profileout(index1)=0;
    profile=profilein+profileout;
    money=sum(profile);
    for i=1:M
        totalProcessTime = totalProcessTime + overTime(i);
   
end

Copy the code

4. Operation results

Matlab version and references

1 matlab version 2014A

[1] Yang Baoyang, YU Jizhou, Yang Shan. Intelligent Optimization Algorithm and Its MATLAB Example (2nd Edition) [M]. Publishing House of Electronics Industry, 2016. [2] ZHANG Yan, WU Shuigen. MATLAB Optimization Algorithm source code [M]. Tsinghua University Press, 2017.