【气动学】改进的遗传和粒子群算法高斯烟羽模型模拟气体扩散【含Matlab源码 1061期】

最新推荐文章于 2024-01-25 10:02:52 发布

Matlab领域

最新推荐文章于 2024-01-25 10:02:52 发布

阅读量1.1k

点赞数 22

分类专栏： Matlab物理应用（初级版）文章标签： matlab

本文链接：https://blog.csdn.net/weixin_63266434/article/details/135819637

版权

Matlab物理应用（初级版）专栏收录该内容

293 篇文章 13 订阅

订阅专栏

在这里插入图片描述

⛄一、获取代码方式

获取代码方式1：
完整代码已上传我的资源：【气动学】基于matlab改进的遗传和粒子群算法高斯烟羽模型模拟气体扩散【含Matlab源码 1061期】
点击上面蓝色字体，直接付费下载，即可。

获取代码方式2：
付费专栏Matlab物理应用（初级版）

备注：
点击上面蓝色字体付费专栏Matlab物理应用（初级版），扫描上面二维码，付费29.9元订阅海神之光博客付费专栏Matlab物理应用（初级版），凭支付凭证，私信博主，可免费获得1份本博客上传CSDN资源代码（有效期为订阅日起，三天内有效）；
点击CSDN资源下载链接：1份本博客上传CSDN资源代码

⛄二、高斯烟羽模型简介

1 高斯羽烟模型方程
在这里插入图片描述
式中c为污染物浓度(单位：kg/m3)
Q为源强(单位：kg/s)
u为泄漏高度的平均风速(单位：m/s)
y、z分别用浓度标准偏差表示的y轴及z轴上的扩散参数
H为泄漏有效高度(单位：m)

2 扩散系数系数
在这里插入图片描述

⛄三、部分源代码

clc;clear;close all;
[xm1,xv1,a1] = mGA(100,2,1000,500,0);
[xm2,xv2,a2] = mPSO(1000,2,2,0.9,0.4,500,3);
[x_zuobiao,y_zuobiao]= gaosiyanyu(5,10,10000.5);
function [xm,xv,accuracy] = mGA(popsize,lenchrom,maxgen,popmax,popmin)
%-------------遗传算法解泄漏源点，效果最好
%% 输入参数
% popsize %种群规模
% lenchrom %变量字串长度
% maxgen % 进化次数
% popmax % 种群最大值
% popmin % 种群最小值
%% 输出参数
% xm 泄漏源坐标,针对泄漏源坐标，论文只考虑地面坐标即xm(1)=x和xm(2)=y
% xv泄漏源强估计值,主要检验方法的有效性
% accuracy方法精确度（%）
bound=[popmin popmax;popmin popmax]; %变量范围
%% 产生初始粒子和速度
for i=1:popsize
%随机产生一个种群
GApop(i,:)=Code(lenchrom,bound); %随机产生个体
%计算适应度
fitvalue(i)=fitness(GApop(i,:)); %染色体的适应度
end

%找最好的染色体
[bestfitness,bestindex]=min(fitvalue);
xm=GApop(bestindex,:); %全局最佳
gbest=GApop; %个体最佳
fitnessgbest=fitvalue; %个体最佳适应度值
xv=bestfitness; %全局最佳适应度值

%% 迭代寻优
for i=1:maxgen
%种群更新 GA选择更新
GApop=Select(GApop,fitvalue,popsize); % 其中Select为选择算子函数

    % 交叉操作 GA
    pc=i/maxgen; % maxgen  进化次数  
    GApop=Cross(pc,lenchrom,GApop,popsize,bound); % 其中Cross为交叉算子函数

    % 变异操作 GA变异
    pm=i/maxgen;
    GApop=Mutation(pm,lenchrom,GApop,popsize,[i maxgen],bound);

    pop=GApop;
    
  for j=1:popsize
    %适应度值即变量约束条件
    if 1*pop(j,1)+0*pop(j,2)<=5000
        if (0*pop(j,1)+1*pop(j,2)<=200) && (0*pop(j,1)-1*pop(j,2)>=-200)
                fitvalue(j)=fitness(pop(j,:));
        end
    end
    %个体最优更新
    if fitvalue(j) < fitnessgbest(j)
        gbest(j,:) = pop(j,:);
        fitnessgbest(j) = fitvalue(j);
    end
    
    %群体最优更新
    if fitvalue(j) < xv
        xm = pop(j,:);
        xv = fitvalue(j);
    end
  end

yy(i)=xv;

end

%% 作图分析
Q = 10000.5;%泄漏源强
accuracy = abs(Q-xv)/Q;
figure;plot(yy,‘linewidth’,2);hold on
xlabel(‘进化代数’);ylabel(‘适应度’);
title(‘改进GA算法收敛曲线’);
legend(‘权重自适应GA算法’)
grid on
[x,y,C] = point;
figure;mesh(x,y,C);
xlabel(‘x轴下风向距离(m)’);ylabel(‘y轴向距离(m)’);zlabel(‘气体扩散浓度’)
figure;surf(x,y,0*C,C,‘edgecolor’,‘none’,‘facecolor’,‘interp’);hold on;plot(xm(1),xm(2),‘r.’,‘MarkerSize’,10)
xlabel(‘x轴下风向距离(m)’);ylabel(‘y轴向距离(m)’);
end

%% 检测函数
function V=fitness(v)
%% 输入参数
x = v(:,1); % 表示泄漏源的x值
y = v(:,2); % 表示泄漏源的y值
%% 输出参数
% V 表示最佳适应度值，同时表示迭代过程中的泄露源强值
c = [112.3840,99.0546,51.0527,8.7177,19.9691,257.4304,…
当然也可以只采用四个；
u = 5;% 风速
Hr = 10; %泄漏点有效高度
z = 0;%高度
Q = 10000.5;%泄漏源强

% 大气稳定度
A
B1=0.16;B2=0.0001;B3=0.12;
C1=0.11;C2=0.0001;C3=0.08;C4=0.0002;
D1=0.08;D2=0.0001;D3=0.06;D4=0.0015;
E1=0.06;E2=0.0001;E3=0.03;E4=0.0003;
F1=0.04;F2=0.0001;F3=0.016;F4=0.0003;
% 选择大气稳定度
w = ‘F’; % 如果w = ‘A’;就表示选择大气稳定度为A的情况进行研究，类似地，w = ‘B’;等。
switch w
% 大气稳定度
A1=0.22;A2=0.0001;A3=0.20;
B1=0.16;B2=0.0001;B3=0.12;
C1=0.11;C2=0.0001;C3=0.08;C4=0.0002;
D1=0.08;D2=0.0001;D3=0.06;D4=0.0015;
E1=0.06;E2=0.0001;E3=0.03;E4=0.0003;
F1=0.04;F2=0.0001;F3=0.016;F4=0.0003;
% 选择大气稳定度
w = F; % 如果w = ‘A’;就表示选择大气稳定度为A的情况进行研究，类似地，w = ‘B’;等。
switch w
case ‘A’

    Ty = ty1*x.*(1+ty2*x).^-0.5;%水平扩散系数
    Tz = tz1*x;%垂直扩散系数
case 'B'
    ty1=B1;ty2=B2;tz1=B3;
    Ty = ty1*x.*(1+ty2*x).^-0.5;%水平扩散系数
    Tz = tz1*x;%垂直扩散系数
case 'C'
    ty1=C1;ty2=C2;tz1=C3;tz2=C4;
    Ty = ty1*x.*(1+ty2*x).^-0.5;%水平扩散系数
    Tz = tz1*x.*(1+tz2*x).^-0.5;%垂直扩散系数
case 'D'
    ty1=D1;ty2=D2;tz1=D3;tz2=D4;
    Ty = ty1*x.*(1+ty2*x).^-0.5;%水平扩散系数
    Tz = tz1*x.*(1+tz2*x).^-0.5;%垂直扩散系数
case 'E'
    ty1=E1;ty2=E2;tz1=E3;tz2=E4;
    Ty = ty1*x.*(1+ty2*x).^-0.5;%水平扩散系数
    Tz = tz1*x.*(1+tz2*x).^-0.5;%垂直扩散系数
case 'F'
    ty1=F1;ty2=F2;tz1=F3;tz2=F4;
    Ty = ty1*x.*(1+ty2*x).^-0.5;%水平扩散系数
    Tz = tz1*x.*(1+tz2*x).^-0.5;%垂直扩散系数

end
%% 检测函数
function Q=fit(u1,hr,data,F)
%% 输入参数
x = data(:,2); % 表示泄漏源的x值
y = data(:,3); % 表示泄漏源的y值
%% 输出参数
% V 表示最佳适应度值，同时表示迭代过程中的泄露源强值
% c = [112.3840,99.0546,51.0527,8.7177,19.9691,257.4304,…
% 53.4442,51.7675,43.7906,28.0270,34.9047,66.5592,…
% 27.5724,27.1807,25.2141,20.6378,22.8799]; % 监测到的浓度，我们用了17个，当然也可以只采用四个；
c=data(:,1);
u = u1;% 风速
Hr = hr; %泄漏点有效高度
z = 0;%高度