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

本文链接：https://blog.youkuaiyun.com/matlab_dingdang/article/details/123806885

1 简介

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

2 部分代码

%*****************应用QPSO算法进行500MW--600MW燃料量对主蒸汽压力系统辨识——QPSO_526_B2PT_new.m%___________________要改参数的地方3处:（含sim引用的模型）用此横线标出、前两处一致、最后是gbestclear;  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));endfor 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);    endendpbest = 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;                                                                %记录本次迭代找到的最好适应值