[mathematical modeling] particle swarm optimization based on MATLAB improved quantum behavior modeling of unit combustion control system [including Matlab source code 1609]

One 、 How to get the code

How to get the code 1：
The complete code has been uploaded to my resources ：【 mathematical modeling 】 be based on matlab Particle swarm optimization algorithm with improved quantum behavior modeling of unit combustion control system 【 contain Matlab Source code 1609 period 】

How to get the code 2：
By subscribing to Ziji Shenguang blog Paid column , With proof of payment , Private Blogger , This code is available .

remarks ：
Subscribe to Ziji Shenguang blog Paid column , Free access to 1 Copy code （ The period of validity From the Subscription Date , Valid for three days ）;

Two 、 Partial source code

%***************** application QPSO Algorithm 500MW--600MW Effect of fuel quantity on main steam pressure system identification ——QPSO_526_B2PT_new.m
clear;  clc;    tic;
load 526select;                 % Load data ：526select Express 500-600MW Select the 250 Group load rise process data 
NUM=250;                        % How many groups of loaded data , Remember to change mdl Internal simulation time 
dimension = 5;                  % Problem dimension -- That is, how many parameters are identified 
popsize = 20;                   % Group size             
MAXITER = 40;                   % Maximum number of iterations               

xmax = [10,300,300,300,100];    % Upper bound of search range , Note that the optimization algorithm defaults to all positive values      
xmin = [-10,0.0,0.0,0.0,0.01];  % Lower bound of search range    
M = (xmax-xmin)/2;              % The median value of the search range 
data1 = zeros(MAXITER);         % Record the best adaptation value of each iteration step 

%***************** Initialize the best position of each particle 、 Global best position 
x = rand(popsize,dimension,1);  % Initialize particles ：popsize That's ok dimension Randomization of columns  
for i = 1:popsize
    for j = 1:dimension
        x(i,j) = (xmax(1,j)-xmin(1,j)) *x(i,j) + xmin(1,j);
    end
end    
pbest = x;                      % Initialize the best position of each particle 、popsize That's ok dimension Column 
gbest = zeros(1,dimension);     % The best place to initialize the global is 0、1 That's ok dimension Column 

%***************** Calculate initial fitness    
for i = 1:popsize               % Calculate the fitness value of the particle at the beginning 
%_______________________________ Where parameters need to be changed _____________________________%   
    k = x(i,1);                 % Substitute the value of the initial particle into mdl Model 
    t1 = x(i,2);
    t2 = x(i,3);    
    t3 = x(i,4);
    tdelay = x(i,5);
                 
    sim('B2PT');                                                    % Conduct simulation 
    yt4to5 = PTout526.signals.values;                               %yt4to5 For the output of the model 
    f_pbest(i) = (PT526lvbo-yt4to5)'*(PT526lvbo-yt4to5);            % Find mean square deviation 
end
    g = find(   f_pbest==min(f_pbest)   );                          % At the beginning, the best global particles are in the first row 
    gbest = pbest(g,:);                                             % Initial global best particle position     
    f_gbest = f_pbest(g);                                           % Record the fitness value of the global best position     
    MINIMUM = f_pbest(g);                                           % Record the best fitness found in each iteration of the algorithm 

%***************** The iteration of the algorithm begins      

      
    for i = 1:popsize
        fi1 = rand(1,dimension);fi2 = rand(1,dimension);u = rand(1,dimension);
            
        p = (2*fi1.*pbest(i,:)+2.1*fi2.*gbest)./(2*fi1+2.1*fi2);    % Calculate the random point attractor 、 I can't see 2 and 2.1 effect       
        b = alpha*abs(mbest-x(i,:));
        v = -log(u);
        x(i,:) = p+((-1).^ceil(0.5+rand(1,dimension))).*b.*v;       % Update of particle position            
     
        z = x(i,:) - (xmax+xmin)/2;                                 % following 3 The row limits the particle position to the search range 
        z = sign(z) .* min(abs(z),M);
        x(i,:) = z + (xmax+xmin)/2;
%_______________________________ Where parameters need to be changed _____________________________%         
        k = x(i,1);
        t1 = x(i,2);
        t2 = x(i,3);    
        t3 = x(i,4);
        tdelay = x(i,5);                          
%***************** Calculate the fitness of the new particle , To determine the optimal particle at each iteration 
        sim('B2PT');                                                % Updated particle , Conduct simulation 
        yt4to5 = PTout526.signals.values;                           % yt4to5 For the output of the model 
        f_x(i) = (PT526lvbo-yt4to5)'*(PT526lvbo-yt4to5);            %  Find mean square deviation 
          
        if (f_x(i)<f_pbest(i))                                      % Update the best position of individual particles 
            pbest(i,:) = x(i,:);
            f_pbest(i) = f_x(i);
        end           
        if f_pbest(i)<f_gbest                                       % Update the global best position of particle swarm 
            gg = [t,i]                                              % Record No. t The position of the optimal particle in the population in the next iteration , The first few .
            gbest = pbest(i,:)                                      %gbest It's a 1*dimension The row vector 
            f_gbest = f_pbest(i);                                   % Record the global best position Fitness , In each iteration , With the circulation of particles , After the particle cycle is complete , The minimum fitness of the population in this iteration is obtained 
        end            
    end                                                             % End a particle  popsize+1
    MSE = sqrt(f_gbest/NUM)                                         % Every iteration  , Update the standard deviation of the global optimal particle once 
    data1(t) = f_gbest;                                             % Record the best fitness found in this iteration 
end                                                                 % End an iteration  t+1  
gbest                                                              
  
%***************** Simulate the global optimal particle again , Draw a curve to compare 
%_______________________________ Where parameters need to be changed _____________________________%
k = gbest(1,1);
t1 = gbest(1,2);
t2 = gbest(1,3);             
t3 = gbest(1,4);
tdelay = gbest(1,5);
   
sim('B2PT');   
yt4to5 = PTout526.signals.values;           
    figure(1);
    plot(t_526,PT526lvbo,'b');                                      % The actual output , Blue solid line 
    hold on;
    plot(t_526,yt4to5,':r');                                        % Identification output , Red dotted line 
    legend(' Actual main steam pressure curve ',' Identify the main steam pressure curve '); 
    grid on;
    xlabel(' Time /s');  
    ylabel(' Main steam pressure /MPa');    
toc;                                                                % Record the best fitness found in this iteration 

3、 ... and 、 Running results

Four 、matlab Edition and references

1 matlab edition
2014a

2 reference
[1] Li Xin .MATLAB mathematical modeling [M]. tsinghua university press .2017
[2] Wang Jian , Zhao Guosheng .MATLAB Mathematical modeling and simulation [M]. tsinghua university press .2016
[3] Yu Shengwei .MATLAB Mathematical modeling classic case practice [M]. tsinghua university press .2015