（论文复现，matlab代码分享，可运行）改进遗传算法求解农业水资源调度问题

问题简介

目前国内外现有的渠道优化配水模型都是在下级渠道配水流量相等这一假定条件的基础上建立的。这与
绝大多数渠系实际配水要求不相符合。针对这一问题，建立了下级渠道引水流量不等情况下的渠道优化配水模型，研究了模型的遗传算法的编码及适应度设计方法。结果表明，渠系优化配水过程较原配水过程斗渠配水时间搭配合理、支渠配水流量均匀‚符合优化配水原则，渠系渗透损失水量比原方案可减少。表明模型及其求解算法是可行的，可为下级引水渠道流量不等时的渠系优化配水决策提供理论和技术支持。

渠系分布示意图

目标函数：

模型以 N 条下级渠道在 T个时段的渠首配水流量为决策变量，将目标函数设置为下级用水户缺水量和渠道输水水量损失最小，即：

1)对农业水资源进行优化的目的是为农田补充水分以保证农作物对水的需求，以此获得较高的产量和较为理想的收益。由于本文研究的是实时来水量不满足该次灌溉需求时的水资源配置，因此，将水资源配置结果最能满足实时灌溉预报为目标建立模型如公式：

式中：Ws 为下级用水户各时段的缺水总量，m3； Qjtd 为下级渠道 Dj 渠首第 t
时段的输水流量，m3/s； Δtj 为为第j条下级渠道的配水时间，s； ,
t1j、t2j分别为第j条下级渠道的配水开始、配水结束时间；Wjt 为第 j
条下级渠道第 t 时段的需水量， m3。

2）渠道渗漏是导致渠系水利用效率低下的主要因素之一，要使农业灌溉的水量损失最小，就必须考虑渠道渗漏的影响，因此，建立干渠及支渠灌溉损失最小模型如公式：

式中：Wl为上级渠道输水的水量损失，m3；qu为上级渠道输水水量损失，m3/s ；
qd为下级渠道输水水量损失，m3/s ；qitu 为上级渠道 Ui 第 t
时段的损失流量，m3/s；qitd 为下级渠道 Di 第 t 时段的损失流量，m3/s。

式中： Aiu 、Aju、 miu 、 mju 分别为上、下级渠道 Ui
、Dj的渠床土壤透水系数和指数； Qitu 为上级渠道 Ui末端第
t时段的输水流量，m3/s； Liu 为上级渠道 Ui 的长度，km。

约束条件

1. 渠道输水能力约束：为了保证渠道安全稳定运行，渠道输水流量不得超过设计流量。

式中： Qius 、 Qjds 分别为上级渠道 Ui 、下级渠道 Dj 的设计流量，m3/s。

1. 水量约束：受来水限制，上级渠道渠首引水量不得超过该干渠分配水量
   
   （轮期内上级渠道来水量不超过50万立方米）

   每条下级渠道的配水流量与引水时间的乘积等于该渠道的配水量。
   

2. 渠段水量平衡约束：

   上级渠道第i渠段末端流量应等于该渠段末端所连的各下级渠道渠首流量与上级渠道第
   i + 1渠段渠首流量之和。
   

3. 非负约束：

1. 轮期约束：

   各下级渠道配水开始时间和结束时间应在轮期内，且开始时间大于零。（轮期为5天）

渠系分布图

所有数据展示

表1 西浚灌区西洞干渠下级渠道基本信息

序号	渠道名称	渠道类型	长度/km	控制灌溉面积/hm2	渠道设计流量/(m3·s-1)	子灌区需水量	下级渠首到上级渠首（km）
1	西洞干渠	干渠	10.23	1259	2.5
2	直属一斗	斗渠	1.80	46	0.6	19000	3.9
3	直属二斗	斗渠	4.20	146	1.0	60000	4.4
4	直属三斗	斗渠	5.80	248	1.0	100000	4.7
5	直属四斗	斗渠	1.25	65	0.6	26000	5
6	直属五斗	斗渠	1.20	76	0.6	32000	5.3
7	直属六斗	斗渠	0.88	55	0.5	22000	5.7
8	直属七斗	斗渠	1.03	41	0.5	17000	5.9
9	直属八斗	斗渠	1.40	117	0.6	48000	6.2
10	直属九斗	斗渠	0.90	230	0.8	95000	6.7
11	西洞支渠	支渠	6.17	233	1.5	96000	10.23
11	毛家湾支渠	支渠	1.80	53	0.8	22000	3

结合结合《农田水利学》，模型系数见表 2。

参数	参数取值
防渗措施折减系数β	0.5
渠床土壤透水系数 Ａ	3.4
渠床土壤透水指数 ｍ	0.5
上级渠道最小流量系数Ｊｄ	0.4
上级渠道加大流量系数Ｊｕ	1.3
下级渠道最小流量系数αｄ	0.6
下级渠道加大流量系数αｕ	1.3

遗传算法设计及代码实现

本次代码大大优化了执行效率，增加了完整注释，值得一读。使用结构数组封装种群，效率更高，可读性大大提升。（如需完整代码，私信留言邮箱）

运行结果

// 主程序入口
%% 清空环境
clc
clear
close all

%% 全局参数
global T Qmax Qmin
T = 3;
Qmax = [0.8,0.6,1,1,0.6,0.6,0.5,0.5,0.6,0.8,1.5];
Qmin = 0.6.*Qmax;

%% 遗传算法参数
maxgen=200;                         %进化代数
sizepop=300;                       %种群规模

%% 个体初始化
population = Init(sizepop);

trace_obj = zeros(1,maxgen);
trace_con = zeros(1,maxgen);

%% 进化开始
for i=1:maxgen
    % 交叉变异
    offspring = Mutate(population,i/maxgen);
    % 挑选新个体
    population = Select(population,offspring,sizepop);
    
    % 记录信息
    bestobj = population(1).obj;
    trace_obj(i) = bestobj;
    trace_con(i) = population(1).cons;
    
    if ~mod(i,10)
        cons = [population.cons];
        num = sum(cons==0);
        avgcons = mean(cons);
        disp(['第' num2str(i) '代，满足约束个体数量：' num2str(num), '，最佳个体：' num2str(bestobj)])
    end
end
%进化结束
bestsol = population(1);

DrawResult()

// 目标函数计算
function population = CalObj(population)
% population = Init(10);
global T
N = length(population);

beta = 0.5;
A = 3.4;
m = 0.5;
Jd = 0.4;
Ju = 1.3;
ad = 0.6;
au = 1.3;

L = [1.8,1.8,4.2,5.8,1.25,1.2,0.88,1.03,1.4,0.9,6.17]; %渠道长度
distance = [3,3.9,4.4,4.7,5,5.3,5.7,5.9,6.2,6.7,10.23]; %下级渠首到上级渠首
Qneed = [22,19,60,100,26,32,22,17,48,95,96].*1000; %子灌区需水量


for i=1:N
    tstart = population(i).tstart; % 解码
    tend = population(i).tend;
    tlast = tend - tstart + 1;
    q = population(i).q;
    
    % 是否处于供水状态
    isgongshui = zeros(T*24/4,11);
    for qu = 1:11
       isgongshui(tstart(qu):tend(qu),qu) = 1; 
    end
    
    % 下级渠道的配水水量
    Q = q.*tlast.*3600.*4;
    
    % 下级渠道的损失流量
    Q_down_loss = 0.01 .* A .* q.^(1-m) .* L .* 1000;
    
    % 用水户缺水量
    Ws = max(Qneed - Q, 0);
    
    % 计算各时间段渠首流量
    Q_t_shou = zeros(T*24/4,11);
    for t=1:T*24/4
        Q_t_shou(t,:) = q(11).*isgongshui(t,11);
        for qu = 10:-1:1
            Q_t_shou(t,qu) = Q_t_shou(t,qu+1) + q(qu)*isgongshui(t,qu);
        end
    end
    
    % 上级渠道的损失流量
    Q_up_loss = 0.01 .* A .* sum(Q_t_shou).^(1-m) .* distance .* 1000;
    
    %% 计算目标函数值
    Wl = sum(Q_down_loss + Q_up_loss);
    F = sum(Ws) + Wl;
    
    %% 检查约束条件
    cons = 0;
    % 水量约束
    if sum(Q) > 500000
       cons = cons + sum(Q)-500000;
    end
    
    % 渠道供水能力约束
    Q_ganqu_t = Q_t_shou(:,1);
    cons = cons + sum(max(Q_ganqu_t-2.5*Ju,0)) + sum(max(2.5*Jd-Q_ganqu_t,0));
    
    %% 封装
    population(i).obj = F;
    population(i).cons = cons;
    population(i).Q = Q;
    population(i).Wl = Wl;
    population(i).Ws = sum(Ws);
    population(i).Q_down_loss = Q_down_loss;
    population(i).Q_up_loss = Q_up_loss;
end

最后，博主专注于论文的复现工作，有兴趣的同学可以私信共同探讨。相关代码已经上传到资源共享，点击我的空间查看分享代码。