A list,
1. Gaussian plume model equation
Where C is pollutant concentration (unit: kg/m3)
Q is source strength (unit: kg/s)
U is the average wind speed at the leakage height (unit: m/s)
The diffusion parameters on y axis and z axis expressed by concentration standard deviation respectively
H is the effective height of leakage (unit: m)
2 diffusion coefficient coefficient
Ii. Source code
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.to (popsize, lenchrom maxgen, popmax, popmin) % -- -- -- -- -- -- -- -- -- -- -- -- -- genetic algorithm solution leakage source points, Best effect %% Input parameter % POPSize % population size % Lenchrom % Variable string length % maxGen % number of evolutions % POPMax % population maximum % POPmin % population minimum %% Output parameter % xm Leak source coordinates. For leak source coordinates, the paper only considers ground coordinates, namely xM (1) = x and xm (2)=y % XV leak source strength estimate, main test method validity % accuracy method accuracy (%) bound=[POPmin popmax;popmin popmax]; The % variable range %% produces the initial particle and speedfor i=1GApop(I,:)=Code(lenchrom,bound); Fitvalue (I)=fitness(GApop(I,:)); End % Find the best chromosome [bestfitness,bestindex]=min(fitValue); xm=GApop(bestindex,:); % global gBest =GApop; % Individual fitnessgBest = fitValue; % Individual optimum fitness xv=bestfitness; % global optimum fitness value %% iterative optimizationfor i=1GApop=Select(GApop, fitValue,popsize); % Select % GA PC = I /maxgen; % maxgen evolutionary times GApop = Cross (PC, lenchrom GApop, popsize, bound). % where Cross is the crossover operator function % mutation operator GA mutation PM = I /maxgen; GApop=Mutation(pm,lenchrom,GApop,popsize,[i maxgen],bound); pop=GApop;for j=1: POPsize % fitness value is the variable constraintif 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 % Individually optimal updateif fitvalue(j) < fitnessgbest(j)
gbest(j,:) = pop(j,:); fitnessgbest(j) = fitvalue(j); End % group optimal updateif fitvalue(j) < xv
xm = pop(j,:);
xv = fitvalue(j);
end
end
yy(i)=xv; End %% graph analysis Q =10000.5; % Leakage source strong accuracy =abs(Q-xv)/Q; figure; plot(yy,'linewidth'.2);hold on
xlabel('Evolutionary algebra'); ylabel('Fitness');
title('Convergence Curve of improved GA Algorithm');
legend('Weight adaptive GA algorithm') grid on [x,y,C] = point; figure; mesh(x,y,C); xlabel('X downwind distance (m)'); ylabel('Y axial distance (m)'); zlabel('Diffusion concentration of gas') figure; surf(x,y,0*C,C,'edgecolor'.'none'.'facecolor'.'interp'); hold on; plot(xm(1),xm(2),'r.'.'MarkerSize'.10)
xlabel('X downwind distance (m)'); ylabel('Y axial distance (m)'); Function V=fitness(V) %%1); % represents the x value of the leakage source y = v(:,2); % represents the y value of the leak source %% The output parameter % V represents the best fitness value, and represents the strong value c = [of the leak source in the iterative process.112.3840.99.0546.51.0527.8.7177.19.9691.257.4304. You could have just four; u =5; % Wind speed Hr =10; % Effective height of leak point z =0; % height Q =10000.5; % leakage source strength % atmospheric stability 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; % Select atmospheric stability w ='F'; % if w ='A'; For the case where atmospheric stability is selected as A, similarly, w ='B'; And so on.switchW % Atmospheric stability 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; % Select atmospheric stability w = F; % if w ='A'; For the case where atmospheric stability is selected as A, similarly, w ='B'; And so on.switch w
case 'A'
Ty = ty1*x.*(1+ty2*x).^0.5; % horizontal diffusion coefficient Tz = tz1*x; % vertical diffusion coefficientcase 'B'ty1=B1; ty2=B2; tz1=B3; Ty = ty1*x.*(1+ty2*x).^0.5; % horizontal diffusion coefficient Tz = tz1*x; % vertical diffusion coefficientcase 'C'ty1=C1; ty2=C2; tz1=C3; tz2=C4; Ty = ty1*x.*(1+ty2*x).^0.5; % horizontal diffusion coefficient Tz = tz1*x.*(1+tz2*x).^0.5; % vertical diffusion coefficientcase 'D'ty1=D1; ty2=D2; tz1=D3; tz2=D4; Ty = ty1*x.*(1+ty2*x).^0.5; % horizontal diffusion coefficient Tz = tz1*x.*(1+tz2*x).^0.5; % vertical diffusion coefficientcase 'E'ty1=E1; ty2=E2; tz1=E3; tz2=E4; Ty = ty1*x.*(1+ty2*x).^0.5; % horizontal diffusion coefficient Tz = tz1*x.*(1+tz2*x).^0.5; % vertical diffusion coefficientcase 'F'ty1=F1; ty2=F2; tz1=F3; tz2=F4; Ty = ty1*x.*(1+ty2*x).^0.5; % horizontal diffusion coefficient Tz = tz1*x.*(1+tz2*x).^0.5; Function Q=fit(u1,hr,data,F) x = data(:,2); % represents the x value of the leak source y = data(:,3); % represents the y value of the leak source %% The output parameter % V represents the best fitness value, and represents the leak source strong value % c = [in the iterative process112.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]; % of the concentration detected, we used17Of course, you can only use four; c=data(:,1); u = u1; % Wind speed Hr = Hr; % Effective height of leak point z =0; % heightCopy the code
3. Operation results
Fourth, note
Version: 2014 a