【数学建模】基于matlab改进量子行为的粒子群算法机组燃烧控制系统建模【含Matlab源码 1609期】

原创

佛怒火莲 2021-12-28 10:19:47 博主文章分类：matlab完整代码 ©著作权

文章标签 matlab 算法开发语言迭代初始化 文章分类 数据结构与算法人工智能

©著作权归作者所有：来自51CTO博客作者佛怒火莲的原创作品，请联系作者获取转载授权，否则将追究法律责任

一、获取代码方式

二、部分源代码

%*****************应用QPSO算法进行500MW--600MW燃料量对主蒸汽压力系统辨识——QPSO_526_B2PT_new.m
clear;  clc;    tic;
load 526select;                 %加载数据：526select表示500-600MW挑选了250组升负荷过程数据
NUM=250;                        %加载数据有多少组,记得改mdl里仿真时间
dimension = 5;                  %问题维数--即辨识多少个参数
popsize = 20;                   %群体规模            
MAXITER = 40;                   %最大迭代次数              

xmax = [10,300,300,300,100];    %搜索范围上界，注意寻优算法默认全是正值     
xmin = [-10,0.0,0.0,0.0,0.01];  %搜索范围下界   
M = (xmax-xmin)/2;              %搜索范围的中值
data1 = zeros(MAXITER);         %记录每一迭代步的最好适应值

%*****************初始化每个粒子最好位置、全局最好位置
x = rand(popsize,dimension,1);  %初始化粒子：popsize行dimension列的随机化 
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;                      %初始化每个粒子最好位置、popsize行dimension列
gbest = zeros(1,dimension);     %初始化全局最好位置为0、1行dimension列

%*****************计算初始适应值   
for i = 1:popsize               %计算初始时粒子的适应值
%_______________________________要改参数的地方_____________________________%   
    k = x(i,1);                 %将初始粒子的值代入mdl模型
    t1 = x(i,2);
    t2 = x(i,3);    
    t3 = x(i,4);
    tdelay = x(i,5);

    sim('B2PT');                                                    %进行仿真
    yt4to5 = PTout526.signals.values;                               %yt4to5为模型的输出
    f_pbest(i) = (PT526lvbo-yt4to5)'*(PT526lvbo-yt4to5);            %求均方差
end
    g = find(   f_pbest==min(f_pbest)   );                          %初始时全局最好粒子在第几行
    gbest = pbest(g,:);                                             %初始时全局最好粒子位置    
    f_gbest = f_pbest(g);                                           %记录全局最好位置的适应值    
    MINIMUM = f_pbest(g);                                           %记录算法每次迭代找到的最好适应值

%*****************算法的迭代开始     


    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);    %计算随机点吸引子、没看到2和2.1作用      
        b = alpha*abs(mbest-x(i,:));
        v = -log(u);
        x(i,:) = p+((-1).^ceil(0.5+rand(1,dimension))).*b.*v;       %粒子位置的更新           

        z = x(i,:) - (xmax+xmin)/2;                                 %以下3行将粒子位置限制在搜索范围
        z = sign(z) .* min(abs(z),M);
        x(i,:) = z + (xmax+xmin)/2;
%_______________________________要改参数的地方_____________________________%         
        k = x(i,1);
        t1 = x(i,2);
        t2 = x(i,3);    
        t3 = x(i,4);
        tdelay = x(i,5);                          
%*****************计算新粒子的适应值,以确定每次迭代时的最优粒子
        sim('B2PT');                                                %更新了粒子，进行仿真
        yt4to5 = PTout526.signals.values;                           % yt4to5为模型的输出
        f_x(i) = (PT526lvbo-yt4to5)'*(PT526lvbo-yt4to5);            % 求均方差

        if (f_x(i)<f_pbest(i))                                      %更新粒子个体最好位置
            pbest(i,:) = x(i,:);
            f_pbest(i) = f_x(i);
        end           
        if f_pbest(i)<f_gbest                                       %更新粒子群全局最好位置
            gg = [t,i]                                              %记录第t次迭代中最优的粒子在种群中的位置，即第几个。
            gbest = pbest(i,:)                                      %gbest是一个1*dimension的行向量
            f_gbest = f_pbest(i);                                   %记录全局最好位置适应值，在每次迭代过程中，随着粒子的循环而变化，粒子循环完成之后，得到的是本次迭代的种群最小适应值
        end            
    end                                                             %结束了一个粒子 popsize+1
    MSE = sqrt(f_gbest/NUM)                                         %每迭代一次 ，就更新一次全局最优粒子的标准差
    data1(t) = f_gbest;                                             %记录本次迭代找到的最好适应值
end                                                                 %结束了一次迭代 t+1  
gbest                                                              

%*****************把全局最优粒子仿真一遍，画出曲线对比
%_______________________________要改参数的地方_____________________________%
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');                                      %实际输出，蓝色实线
    hold on;
    plot(t_526,yt4to5,':r');                                        %辨识输出，红色虚线
    legend('实际主蒸汽压力曲线','辨识主蒸汽压力曲线'); 
    grid on;
    xlabel('时间/s');  
    ylabel('主蒸汽压力/MPa');    
toc;                                                                %记录本次迭代找到的最好适应值