【微电网优化】基于量子行为粒子群算法机组燃烧控制系统建模含Matlab源码

Matlab科研辅导帮

于 2022-03-28 21:33:25 发布

阅读量15

点赞数

文章标签：算法 matlab 开发语言

某大厂资深算法工程师，从事Matlab算法仿真工作10年,擅长智能优化算法、神经网络预测、机器学习、信号处理、元胞自动机、图像处理、路径规划、无人机、无线传感器网络、车间调度、生产调度等多种领域的Matlab仿真，更多仿真源码、算法改进、Matlab项目和期刊发表可私信合作。

本文链接：https://blog.csdn.net/m0_60703264/article/details/140030734

版权

1 简介

能源问题与环境问题随着现代社会的快速发展已成为中国乃至全世界关注的焦点。就我国现状来说,由于独特的能源架构和社会形态,直接决定了我国的电力工业在当今乃至未来相当长的一段时期内将以燃煤火电机组为主。伴随着电力工业架构的调整,国内的电源项目集中建设大型火电机组,着重发展大容量、高效率、高参数的超临界火电机组。处于发展中、又可靠成熟的超临界火电机组已经被国际纳入洁净煤燃烧技术的轨道,相比于传统的亚临界火力发电机组,超临界机组的发电煤耗率大大减少,更加具有节约和环保意义。本文在研究了超临界机组燃烧系统的工作原理、控制任务以及生产工艺流程的基础之上,给出了西安某电厂600MW超临界机组各个设备的概述及运行规范。整个火电厂的能量来源是由超临界机组燃烧系统将燃料燃烧产生的热能转化的蒸汽动能。本文以超临界机组燃烧控制系统为研究对象,集中研究系统建模和优化控制两个方面。精准的系统模型赋予了超临界机组燃烧控制系统更高的控制水平,同时能正确反映各个热力设备的动态响应过程及相关参数,充分了解机组运行的动态特性。我们可以根据现场采集到的运行历史数据,应用智能算法对超临界机组燃烧控制系统进行模型辨识,进而得到符合现场要求的数学模型。

2 部分代码


          
          
           
           %*****************应用QPSO算法进行500MW--600MW燃料量对主蒸汽压力系统辨识——QPSO_526_B2PT_new.m
           
           
%___________________要改参数的地方3处:(含sim引用的模型)用此横线标出、前两处一致、最后是gbest
           
           
clear;  clc;    tic;
           
           
load 526select;                 %加载数据：526select表示500-600MW挑选了250组升负荷过程数据
           
           
NUM=250;                        %加载数据有多少组,记得改mdl里仿真时间
           
           
dimension = 15;                  %问题维数--即辨识多少个参数
           
           
popsize = 30;                   %群体规模            
           
           
MAXITER = 1000;                   %最大迭代次数              
           
           
xmax = [10,10,10,300,300,300,300,300,300,300,300,300,100,100,100];    %搜索范围上界，注意寻优算法默认全是正值     
           
           
xmin = [-10,-10,-10,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0001,0.0001,0.0001];  %搜索范围下界  
           
           
M = (xmax-xmin)/2;              %搜索范围的中值
           
           
data1 = zeros(MAXITER);         %记录每一迭代步的最好适应值
           
           
%*****************初始化每个粒子最好位置、全局最好位置
           
           
x = zeros(popsize,dimension);
           
           
r = zeros(1,dimension);
           
           
for j = 1:dimension
           
           
    r(j)=exp(j)-floor(exp(j));
           
           
end
           
           
for i = 1:popsize
           
           
    for j = 1:dimension
           
           
        x(i,j)= r(j)*i - floor(r(j)*i);
           
           
        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               %计算初始时粒子的适应值
           
           
%_______________________________要改参数的地方_____________________________%   
           
           
    k1 = x(i,1);                 %将初始粒子的值代入mdl模型
           
           
    k2 = x(i,2);
           
           
    k3 = x(i,3);
           
           
    t11 = x(i,4);
           
           
    t12 = x(i,5);    
           
           
    t13 = x(i,6);
           
           
    t21 = x(i,7);
           
           
    t22 = x(i,8);    
           
           
    t23 = x(i,9);
           
           
    t31 = x(i,10);
           
           
    t32 = x(i,11);    
           
           
    t33 = x(i,12);
           
           
    tdelay1 = x(i,13);
           
           
    tdelay2 = x(i,14);
           
           
    tdelay3 = x(i,15);
           
           
    sim('BFASA2PT');                                                    %进行仿真
           
           
    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 t = 1:MAXITER            
           
           
    t                                                               %在workspace显示当前运行至第几轮  
           
           
    alpha = (1.0-0.5)*(MAXITER-t)/MAXITER+0.5;                      %收缩－扩张系数计算       
           
           
    mbest = sum(pbest)/popsize;                                     %计算所有最优粒子的平均值
           
           
    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;
           
           
%_______________________________要改参数的地方_____________________________%         
           
           
          k1 = x(i,1);                 
           
           
          k2 = x(i,2);
           
           
          k3 = x(i,3);
           
           
          t11 = x(i,4);
           
           
          t12 = x(i,5);    
           
           
          t13 = x(i,6);
           
           
          t21 = x(i,7);
           
           
          t22 = x(i,8);    
           
           
          t23 = x(i,9);
           
           
          t31 = x(i,10);
           
           
          t32 = x(i,11);    
           
           
          t33 = x(i,12);
           
           
          tdelay1 = x(i,13);
           
           
          tdelay2 = x(i,14);
           
           
          tdelay3 = x(i,15);                        
           
           
%*****************计算新粒子的适应值,以确定每次迭代时的最优粒子
           
           
        sim('BFASA2PT');                                                %更新了粒子，进行仿真
           
           
        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                                                              
           
           
%*****************把全局最优粒子仿真一遍，画出曲线对比
           
           
%_______________________________要改参数的地方_____________________________%
           
           
    k1 = gbest(1,1);               
           
           
    k2 = gbest(1,2);
           
           
    k3 = gbest(1,3);
           
           
    t11 = gbest(1,4);
           
           
    t12 = gbest(1,5);
           
           
    t13 = gbest(1,6);
           
           
    t21 = gbest(1,7);
           
           
    t22 = gbest(1,8);
           
           
    t23 = gbest(1,9);
           
           
    t31 = gbest(1,10);
           
           
    t32 = gbest(1,11);
           
           
    t33 = gbest(1,12);
           
           
    tdelay1 = gbest(1,13);
           
           
    tdelay2 = gbest(1,14);
           
           
    tdelay3 = gbest(1,15);
           
           
sim('BFASA2PT');   
           
           
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;                                                                %记录本次迭代找到的最好适应值