目录
%% 清除变量,导入数据
clear, clc;
load('Data_Create_For_Test.mat');
Data_real = Data_Create;
% 如果想导入EXCEL的数据，请注释掉上面两段代码，使用下面的一段代码
% Data_real = xlsread('Data_Create_For_Test.xlsx'); %导入EXCEL表格的真实数据
%% 初始化参数
Pc = 0.8; %交叉概率
Pm = 0.01; %变异概率
iter_max = 200; %最大进化代数，即最大迭代次数
iter = 1; %初始化进化代数，即当前迭代次数
N_chrom = 20; %种群染色体的数目，解的个数
N_gen = 3; %染色体基因的数目，自变量的个数（表示论文里需反演的三个参数Alpha, Beta, S）
Up_gen = [0.5; 0.01; 0.01]; %三个基因的上界取值范围（三个参数的取值范围, 根据实际情况设定）
Down_gen = [0; 0; 0]; %三个基因的下界取值范围
%% 生成初始解
chrom = zeros(N_chrom, N_gen); %初始化染色体矩阵（解集）
for i = 1 : N_chrom %初始化N_chrom组的参数值
chrom(i, 1) = Down_gen(1) + rand(1) * (Up_gen(1) - Down_gen(1)); %初始化Alpha(i)
chrom(i, 2) = Down_gen(2) + rand(1) * (Up_gen(2) - Down_gen(2)); %初始化Beta(i)
chrom(i, 3) = Down_gen(3) + rand(1) * (Up_gen(3) - Down_gen(3)); %初始化S(i)
end
fitness = zeros(N_chrom, 1); %初始化适应度矩阵（每一个解对应的适应度大小）
for i = 1 : N_chrom
fitness(i) = CalFitness(Data_real, chrom(i, : )); %计算每一个初始解对应的适应度大小
end
fitness_best = zeros(iter_max, 1); %初始化每一代的最优适应度大小
for i = 1 : iter_max
fitness_best(i) = inf; %初始化为无穷大
end
chrom_best = zeros(iter_max, N_gen); %初始化每一代的最优染色体
%% 进化迭代
[row, line] = size(chrom); %初始化选择操作中父本和母本的对数
dad_mom = zeros(2 * row, line); %初始化父本和母本,这里奇数行号认为是父本，偶数行号认为是母本
New_chrom = zeros(2 * row, line); %初始化交叉操作生成的新种群
for iter = 1 : iter_max
%% 选择操作（论文中的竞争选择法，锦标赛选择法）
dad_mom = Selection(chrom, fitness);
%% 交叉操作（论文中的离散交叉法）
New_chrom = Cross(dad_mom, Pc);
%% 变异操作（论文中的非均匀变异法）
New_chrom = Mutation(New_chrom, Up_gen, Down_gen, iter, iter_max, Pm);
%% 计算新种群的适应度大小
[row, line] = size(New_chrom);
New_fitness = zeros(row, line);
for i = 1 : row
New_fitness(i) = CalFitness(Data_real, New_chrom(i, : )); %计算每一个新染色体的适应度大小
if New_fitness(i) < fitness_best(iter) %判断当前新染色体的适应度，是否优于本迭代次数的最优适应度
fitness_best(iter) = New_fitness(i); %判断结果为是，则更新当前迭代次数的最优适应度
chrom_best(iter, : ) = New_chrom(i, : ); %同时更新最优染色体
end
end
if ( (iter ~= 1) && ( fitness_best(iter - 1) < fitness_best(iter) ) )
fitness_best(iter) = fitness_best(iter - 1); %记录截止到当前迭代次数，最优的适应度大小
chrom_best(iter, : ) = chrom_best(iter - 1, : ); %记录截止到当前迭代次数，最优的适应度对应的最优染色体
end
%% 形成新的种群，用于下一次迭代（自然选择）
Sort_fitness = sort(New_fitness); %排序新种群的适应度（注意这里适应度值的个数是N_chrom的两倍）
for i = 1 : (row/2)
index = find(New_fitness == Sort_fitness(i), 1);
chrom(i, : ) = New_chrom(index, : ); %取适应度排名前50%的染色体，作为下一次迭代的新种群，这里可以理解为自然选择操作
end
end
%% 绘图
figure(1)
plot(1:iter_max, fitness_best, 'r');
title('迭代图')
xlabel('迭代次数')
ylabel('适应度大小')
%% 输出结果
disp(['最优适应度:' num2str( fitness_best(iter) )])
disp(['Alpha最优参数值:' num2str( chrom_best(iter, 1) )])
disp(['Beta最优参数值:' num2str( chrom_best(iter, 2) )])
disp(['S最优参数值:' num2str( chrom_best(iter, 3) )]) %% 计算适应度大小
function Fitness_Chrom = CalFitness(Data_real, chrom)
%输入：真实观测数据Data_real（两列），包括时间t（第一列数据）和真实水位抬升值s（第二列数据）
% 染色体解集chrom（三列，分别代表三个参数alpha，beta，S）
%输出：适应度大小
M = 60; %含水层厚度
Q = 58; %定流量
r = 0.1; %井径
A = Q / (4 * pi * M); %论文中的斜率A
B = A * log( 2.25*M / (r*r) ); %论文中的截距B
Data_t = Data_real(:, 1); %取第一列数据，即时间t
Data_s = Data_real(:, 2); %取第二列数据，即观测到的真实水位抬升值s
Alpha = chrom(:, 1); %取第一列数据，即参数alpha
Beta = chrom(:, 2); %取第二列数据，即参数Beta
S = chrom(:, 3); %取第三列数据，即参数S
Data_num = size(Data_s, 1); %计算真实观测数据的数目
Error_sum = 0; %总误差的平方和，即优化目标
X = zeros(Data_num, 1);
for j = 1 : Data_num
X(j) = -1*Beta*Data_t(j) + log( Alpha*Data_t(j) ) - log( S ); %计算时刻t(i)对应的X(i)
error = Data_s(j) - ( A*X(j)+B ) / ( Alpha * exp( -1*Beta*Data_t(j) ) ); %计算时刻t(i)对应的误差（s-s*)
Error_sum = Error_sum + abs(error * error); %计算总误差的平方和
end
Fitness_Chrom = Error_sum; %记录染色体（一组参数）对应的优化目标值function dad_mom = Selection(chrom, fitness)
%% 选择操作（竞争选择法，锦标赛选择法）
% 输入：染色体种群， 适应度矩阵
% 输出：选择操作得到的父本和母本（用于交叉操作）
sn = 3; %从种群中选择的个体数目，进行竞争
[row, line] = size(chrom); %选择父本和母本的对数
dad_mom = zeros(2 * row, line); %初始化父本和母本,这里奇数行号认为是父本，偶数行号认为是母本
for i = 1 : 2 : 2*row %间距为2
r = randi([1 row], sn, 1); %从[1,sum]中随机生成sn个整数,即随机抽出来进行比赛的个体编号
index = find(fitness == min(fitness(r)), 1); %在进行比赛的个体中，找到适应度最优的个体编号，即在原种群chrom中对应的编号
dad_mom(i, : ) = chrom(index, : ); %将比赛胜出者作为父本
dad_mom(i+1, : ) = chrom(randi([1, row]), : ); %随机选择一个个体作为母本
end
function New_chrom = Cross(dad_mom, pc)
%% 交叉操作（离散交叉法）
%输入：选择操作得到的父本和母本；交叉概率
%输出：交叉操作后的新种群(子本)
[row, line] = size(dad_mom);
New_chrom = zeros(row, line); %初始化交叉操作的新种群
for i = 1 : 2 : row
if rand < pc %生成一个随机数，判断是否进行交叉操作
for j = 1 : line %进行交叉操作
if randi([0, 1]) == 1 %随机数为1，则对应基因不进行交叉操作
New_chrom(i, j) = dad_mom(i, j);
New_chrom(i+1, j) = dad_mom(i+1, j);
else %随机数为0，则相邻两个个体(即父本和母本)，进行相应基因的交叉操作（互换基因）
New_chrom(i, j) = dad_mom(i+1, j);
New_chrom(i+1, j) = dad_mom(i, j);
end
end
else %不进行交叉操作
New_chrom(i, : ) = dad_mom(i, : ); %子代保持不变
New_chrom(i, : ) = dad_mom(i+1, :); %子代保持不变
end
end
function New_chrom = Mutation(Cross_chrom, Up_gen, Down_gen, iter, iter_max, pm)
%% 变异操作（非均匀变异法）
%输入：交叉操作得到的新种群；基因的上界；基因的下界；当前迭代次数；最大迭代次数；变异的概率
%输出：变异操作后的新种群
[row, line] = size(Cross_chrom);
New_chrom = zeros(row, line); %初始化变异操作得到的新种群
b = 2; %非均匀变异公式中的常数
a = 1 - iter/iter_max; %此常数随当前迭代次数的变化而变化，具体见非均匀变异中的公式
for i = 1 : row %遍历每一个染色体，即每一个解
for j = 1 : line %遍历每一个基因，即要求解的参数(这里可以根据个人情况修改代码，如果你只想变异一个特定的基因，只需只定j为特定常数)
r = rand; %生成一个随机数，表示非均匀变异公式中随机数
if rand < pm %判断是否发生变异操作（此处为发生变异操作）
if randi([0, 1]) == 1
New_chrom(i, j) = Cross_chrom(i, j) + ( Up_gen(j) - Cross_chrom(i, j) ) * ( 1 - r^(a*b) ); %公式1
else
New_chrom(i, j) = Cross_chrom(i, j) + ( Cross_chrom(i, j) - Down_gen(j) ) * ( 1 - r^(a*b) ); %公式2
end
else
New_chrom(i, j) = Cross_chrom(i, j);
end
end
end
%% 生成测试数据
clear, clc;
M = 60; %含水层厚度
Q = 58; %定流量
r = 0.1; %井径
A = Q / (4 * pi * M); %论文中的斜率A
B = A * log( 2.25*M / (r*r) ); %论文中的截距B
Alpha = 0.346; %参数1
Beta = 0.0044; %参数2
S = 0.00192; %参数3
Sum = 15;
Data_Create = zeros(Sum, 2);
X = zeros(Sum, 1);
for i = 1 : 15
Data_Create(i, 1) = i;
X(i) = -1*Beta*Data_Create(i, 1) + log( Alpha * Data_Create(i, 1) ) - log( S );
Data_Create(i, 2) = ( A*X(i)+B ) / ( Alpha * exp( -1*Beta*Data_Create(i, 1) ) );
end
save Data_Create_For_Test.mat Data_Create;
