【运动学】基于改进的遗传和粒子群算法高斯烟羽模型模拟气体扩散matlab源码...-CSDN博客

本文链接：https://blog.csdn.net/Matlab_dashi/article/details/140115973

一、简介

1 高斯羽烟模型方程【运动学】基于改进的遗传和粒子群算法高斯烟羽模型模拟气体扩散matlab源码_运动学
式中c为污染物浓度(单位：kg/m3)
Q为源强(单位：kg/s)
u为泄漏高度的平均风速(单位：m/s)
y、z分别用浓度标准偏差表示的y轴及z轴上的扩散参数
H为泄漏有效高度(单位：m)

2 扩散系数系数【运动学】基于改进的遗传和粒子群算法高斯烟羽模型模拟气体扩散matlab源码_运动学_02

二、源代码

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;%高度