一种基于粒子群算法的渠系优化配水方法与流程

本发明属于智能农业技术领域，特别涉及一种基于粒子群算法的渠系优化配水方法。

背景技术：

我国作为一个农业大国，农业用水量已经占到全国用水总量的72％，其中大部分都用于灌溉，而灌溉水利用效率极低，很多灌区的灌溉水资源利用率仅仅为40％，这与我国当前努力建设可持续发展的节约型社会是不相符的。因此，在灌溉水资源供需矛盾日益严重的背景下，在农业生产中占据重要位置的灌溉水资源的合理优化配置的重要性和必要性日益显现，为了能更好地促进农业的可持续发展，必须提高灌溉水的利用效率。灌区渠系配水的优化属于其中一个重要方面。

当前，我国多数灌区在进行渠系配水采用的是经验方法，从而产生渠系输水流量小、流量不集中和输水损失大等问题。如何通过优化渠道的流量和配水时间来减小渠道系统的输配水和渗水损失，是灌区管理层面急需解决的问题。

卢麾基于对作物需水量和渠系输配水损失加以考虑，构建了配水优化模型，并使用了遗传算法进行模型求解。同时，配合该模型开发了有GIS功能的灌溉决策支持系统，该系统能实现在满足作物需水的条件下输配水损失达到最小。徐建新考虑了降水的随机特征和作物不同生育期内土壤最低含水率的变化，研制了灌区水量实时调配决策软件。孟爽在VC#.NET的环境下，基于ArcGIS Engine组件包开发了具有查询统计、求解计算、地图操作、水情监测和水资源管理等功能的灌区信息管理系统。佛罗里达大学研发了基于作物种类、作物种植面积、作物生产期、土壤含水率以及灌区灌溉制度的可用于计算作物不同生育时期需水量的AFSIRS系统。

当前在渠系配水模型与灌区灌溉管理方法方面的研究存在诸多问题：

(1)适用于不同类型、不同地域、不同灌区的个体模型构建与算法研究仍有待加强，需因地制宜，增强针对性。

(2)一些模型所需要的参数众多且获取难度大，使得模型的实用性和通用性较低，影响了进一步的推广。

(3)相配套的灌区用水管理方法的研究尚处于初步阶段，需加大力度进行改进和整合。

技术实现要素：

本发明的目的是提供一种适用于不同类型、不同地域、不同灌区的，实用性和通用性高的基于粒子群算法的渠系优化配水方法。

为了实现上述目的，本发明提供了如下技术方案：

本发明提供一种基于粒子群算法的渠系优化配水方法，包括如下步骤：

步骤一，根据上级渠道流量Q、渠道水利用系数η、下级渠道数量N和下级渠道流量q_j构建目标函数，并确定轮灌组数目范围以及粒子群算法运行所需的粒子群粒子总数和最大迭代次数；

步骤二，利用粒子群算法在可行域内产生粒子群；

步骤三，判断步骤二产生的粒子群中的粒子是否满足约束条件，如果粒子不满足约束条件，则摈弃，并返回步骤二继续产生新的下代粒子，直到满足约束条件的粒子数目达到步骤一中的粒子群粒子总数为止；

步骤四，记录所有满足约束条件的粒子的位置；

步骤五，计算粒子的适应度值，并依适应度值由大到小为粒子排序，将适应度值最大的粒子作为备选优化结果；

步骤六，判断粒子群算法是否达到最大迭代次数；如果粒子群算法达到了步骤一确定的最大迭代次数，则输出优化结果；否则，更新该粒子群的所有粒子的速度和位置信息，并重新进行步骤五和步骤六。

所述步骤一中，目标函数为：

min(ω₁Z+ω₂ΔT)

其中，

Z表示实际总轮灌组数Z与下级渠道数量N的比值；

ΔT表示各个轮灌组引水时间相差最大值ΔT与当实际总轮灌组数Z为下级渠道数N时各轮灌组引水时间相差最大值ΔT的取值m的比值；

ω₁、ω₂为权重系数，ω₁+ω₂＝1，ω₁和ω₂分别代表了该目标函数的两个目标所占的比重。

ω₁和ω₂的取值优选为ω₁＝0.4，ω₂＝0.6；

该目标函数有两个目标，一为使总轮灌组数最少，二为使各轮灌组间灌水时间相差最小，分别由如下公式确定：

公式中，Z表示实际总轮灌组数Z与下级渠道数量N的比值；

i和j分别表示轮灌组序号和下级渠道序号；

M为轮灌组数目；

C_i表示第i个轮灌组的损失，现假定各个轮灌组损失相同，因此取相同值为1；

f_i为轮灌组存在系数，当时，f_i＝1,否则，f_i＝0；

N为下级渠道数量；

X_ij为决策变量，是0，1变量，其中，

X_ij＝0表示上级渠道在i轮灌组向第j条下级渠道不配水；

X_ij＝1表示上级渠道在i轮灌组向第j条下级渠道配水；

ΔT表示各个轮灌组引水时间相差最大值ΔT与当实际总轮灌组数Z为下级渠道数N时各轮灌组引水时间相差最大值ΔT的取值m的比值；m为Z＝N时ΔT的取值；

ΔT＝max(T_i-T_k),1≤i≤M,1≤k≤M,i≠k

ΔT表示各轮灌组引水时间相差的最大值；

m为Z＝N时ΔT的取值；

i和k表示不同的轮灌组序号。

所述步骤一中，轮灌组数目M的取值范围由下式确定：

M≤Qη/q_j

其中，Q为上级渠道流量，m³/s；η为渠道水利用系数；则Qη为其净引水流量；q_j为下级渠道流量。

所述步骤一中，粒子群算法运行所需的粒子群粒子总数为100。

所述步骤一中，最大迭代次数为200。

所述步骤三中，约束条件包括：

(1)轮期约束：在任意一个轮灌组内，所有下级渠道的总配水时间小于等于灌水轮期T₀；

其中，T₀为灌水轮期，N为下级渠道数量；t_j为第j条下级渠道的配水所需时间，i和j分别表示轮灌组序号和下级渠道序号；X_ij为决策变量；

(2)下级渠道一次性引水约束：任意一条下级渠道在上级渠道的一个轮期内进行配水时，渠道闸门只被打开一次，并在引水时间内不间断地持续引水；

其中，i和j分别表示轮灌组序号和下级渠道序号；X_ij为决策变量；M为轮灌组数目；

(3)决策变量约束：

X_ij＝{0,1}

其中，X_ij＝0表示上级渠道在i轮灌组向第j条下级渠道不配水；X_ij＝1表示上级渠道在i轮灌组向第j条下级渠道配水；

(4)来水量约束：各个轮灌组总流量之和不超过上级渠道来水净流量；

其中，i和j分别表示轮灌组序号和下级渠道序号；M为轮灌组数目；N为下级渠道数量；q_i为第i个轮灌组的引水流量，m³/s；f_i为轮灌组存在系数，当时，f_i＝1,否则，f_i＝0；Q为上级渠道的引水流量，m³/s，η为渠道水利用系数，则Qη为其净引水流量。

所述步骤五中，粒子的适应度值通过如下适应度函数F计算得到：

所述步骤六中，粒子的位置和速度信息的更新公式分别为：

其中，

为第e次迭代进行中第i个粒子在n维空间中的位置；

为第e+1次迭代进行中第i个粒子在n维空间中的位置；

r为速度的约束因子，通常设置为1；

为第e次迭代进行中第i个粒子在n维空间中的速度；

为第e+1次迭代进行中第i个粒子在n维空间中的速度；

θ为惯性权重，用来表示保持原来速度的系数；

c₁，c₂分别为粒子跟踪个体最优的权重系数和跟踪群体最优的权重系数，通常都设置为2；

h₁,h₂为[0,1]区间内均匀分布的随机数；

为第e次迭代进行中第i个粒子在n维空间中迄今为止搜索到的最优位置；

为第e次迭代进行中粒子群体在n维空间中迄今为止搜索到的最优位置。

与现有技术相比，本发明的有益效果在于：

本发明运用了灌溉渠系配水优化数学模型和高效、稳健的模型求解算法，建立具有较好的通用性，适用于多数灌区的区域灌溉渠系的优化配水决策模型。

本发明由计算机程序实现，包括输入模块，功能模块，数据库管理模块，概化图产生模块，参数传递模块，求解模块，输出模块等，实现与MATLAB编写的渠系优化配水模型嵌套。同时软件力求有良好的用户界面，简洁友好的人工对话方式从而方便工作人员的操作。

附图说明

图1为本发明基于粒子群算法的渠系优化配水方法流程图；

图2为本发明实施例的配水过程图；

图3a为本发明实施例的第一次求解结果的轮灌组1的分组情况图；

图3b为本发明实施例的第一次求解结果的轮灌组2的分组情况图；

图3c为本发明实施例的第一次求解结果的轮灌组3的分组情况图；

图3d为本发明实施例的第一次求解结果的轮灌组4的分组情况图。

其中的附图标记为：

1-12分别为下级斗渠序号

具体实施方式

下面结合附图和实施例对本发明进行进一步说明。

如图1所示，本发明基于粒子群算法的渠系优化配水方法包括如下步骤：

步骤一，根据上级渠道流量Q、渠道水利用系数η、下级渠道数量N和下级渠道流量q_j构建目标函数，并确定轮灌组数目范围以及粒子群算法运行所需的粒子群粒子总数和最大迭代次数。

目标函数为：

min(ω₁Z+ω₂ΔT) (公式1)

其中，

Z表示实际总轮灌组数Z与下级渠道数量N的比值；

ΔT表示各个轮灌组引水时间相差最大值ΔT与当实际总轮灌组数Z为下级渠道数N时各轮灌组引水时间相差最大值ΔT的取值m的比值；

ω₁、ω₂为权重系数，ω₁+ω₂＝1，ω₁和ω₂分别代表了该目标函数的两个目标所占的比重，ω₁和ω₂的取值优选为ω₁＝0.4，ω₂＝0.6。

目标函数的两个目标：一为使总轮灌组数最少，二为使各轮灌组间灌水时间相差最小，分别由如下公式确定：

公式2中，Z表示实际总轮灌组数Z与下级渠道数量N的比值；

i和j分别表示轮灌组序号和下级渠道序号；

M为轮灌组数目；

C_i表示第i个轮灌组的损失，现假定各个轮灌组损失相同，因此取相同值为1；

f_i为轮灌组存在系数，当时，f_i＝1,否则，f_i＝0；

N为下级渠道数量；

X_ij为决策变量，是0，1变量，其中，

X_ij＝0表示上级渠道在i轮灌组向第j条下级渠道不配水；

X_ij＝1表示上级渠道在i轮灌组向第j条下级渠道配水。

公式3中，ΔT表示各个轮灌组引水时间相差最大值ΔT与当实际总轮灌组数Z为下级渠道数N时各轮灌组引水时间相差最大值ΔT的取值m的比值；因为实际总轮灌组数Z与各轮灌组引水时间差的最大值ΔT量纲不同，故模型的目标采用了无量纲化处理，用比值形式的Z和ΔT来替代。可以得到实际的最小轮灌组数和最小的各轮灌组引水时间的最大差值；(m为Z＝N时ΔT的取值)。

ΔT＝max(T_i-T_k)(1≤i≤M,1≤k≤M,i≠k)

ΔT表示各轮灌组引水时间相差的最大值

m为Z＝N时ΔT的取值；

i和k表示不同的轮灌组序号。

确定轮灌组数目范围：

轮灌组数目M的取值范围由下式确定：

M≤Qη/q_j (公式4)

其中，Q为上级渠道流量，m³/s，η为渠道水利用系数，则Qη为其净引水流量，q_j为下级渠道流量，具体为第i个轮灌组的引水流量(q₁＝q₂＝…＝q_N)。

确定粒子群算法运行所需的粒子群粒子总数：

粒子群算法中诸多参数的选择带有经验性。粒子群大小代表每一代粒子群当中所包含的粒子的数目。面对不同的优化问题时，最佳粒子群规模并不相同。当选取的值太小时，算法的运算速度较快，但受限于较差的种群多样性，容易产生过早收敛，陷入局部最优；而当选取的K值太大时，算法的运算速度和运行寻优效率会降低。在实际应用中，为了充分考虑种群的丰富多样性，算法的搜索效率和足够的解空间等种种因素，人们通常采用经验规则，选取的值为粒子个体维数的5～10倍，本发明选取粒子群粒子总数为100。

确定最大迭代次数：

本发明选取最大迭代次数(总迭代次数)为200，在用实例进行验证时，产生的非利解在误差范围内，算法的鲁棒性表现不错。

步骤二，利用粒子群算法在可行域内产生粒子群。

粒子群算法的思想是：在一个优化问题的求解过程中，每个潜在解都是搜索空间中的一只鸟，我们称之为粒子，并为每个粒子规定了相应的行为规则，每个粒子都由一个依据目标函数设计的适应度值引导飞行。与此同时，每个粒子还有一个速度来决定它们飞行的距离和方向。单个粒子通过追随最优粒子在解空间中进行搜索，最终整个群体通过相互之间的个体协作降落在最优解位置，从而实现对问题的求解。

用产生随机数的方法在解空间内产生初始粒子。本发明中设定粒子群数目为100。即产生100个M行N列的矩阵，即粒子群。

对于M＝7，N＝7，则生成100个7×7的矩阵。每个矩阵即为一个粒子。

比如7×7维的粒子，其表现形式如以下粒子1和粒子2所示，粒子1

粒子2

即每一列有且仅有1个元素为1。

步骤三，判断步骤二产生的粒子群中的粒子是否满足下述约束条件，如果粒子不满足约束条件，则摈弃，并返回步骤二继续产生新的下代粒子，直到满足约束条件的粒子数目达到步骤一中的粒子群粒子总数为止。

约束条件包括：

(1)轮期约束：在任意一个轮灌组内，所有下级渠道的总配水时间必须小于等于灌水轮期T₀；

其中，T₀为灌水轮期，N为下级渠道数量；t_j为第j条下级渠道的配水所需时间，i和j分别表示轮灌组序号和下级渠道序号；X_ij为决策变量。

(2)下级渠道一次性引水约束：任意一条下级渠道在上级渠道的一个轮期内进行配水时，渠道闸门只被打开一次，并在引水时间内不间断地持续引水；

其中，i和j分别表示轮灌组序号和下级渠道序号；X_ij为决策变量；M为轮灌组数目。

(3)决策变量约束：

X_ij＝{0,1} (公式7)

其中，X_ij＝0表示上级渠道在i轮灌组向第j条下级渠道不配水；X_ij＝1表示上级渠道在i轮灌组向第j条下级渠道配水。

(4)来水量约束：各个轮灌组总流量之和不超过上级渠道来水净流量

其中，i和j分别表示轮灌组序号和下级渠道序号；M为轮灌组数目；N为下级渠道数量；q_i为第i个轮灌组的引水流量，m³/s；f_i为轮灌组存在系数，当时，f_i＝1,否则，f_i＝0；Q为上级渠道的引水流量，m³/s，η为渠道水利用系数，则Qη为其净引水流量。

本发明在约束条件的处理上，合理地利用了约束条件，从而使得所建立的最优化模型解空间大大缩小。传统的渠系优化配水模型求解方法大都将其视为一个一般化的最优化问题，由于渠道数量较多，导致解空间是巨大的，因而寻优过程消耗的时间也是相当大的。而本发明从一次性引水约束出发，得出上级渠道所对应的各条下级渠道在某一个配水轮期内只能有一条配水，从而排除了很多潜在的解，使得离散PSO算法的寻优过程变得更为有效。对约束条件的处理上，轮灌组数目约束没放在约束条件里面，而是采用逐次遍历轮灌组数目M的取值的方法处理，即让M从1取到最大值。

步骤四，记录所有满足约束条件的粒子的位置

例如某一个7×7离散粒子X为

基于纵坐标的坐标赋予粒子位置，则X对应的位置(即粒子矩阵中每列元素1所在行的行号)为(2,2,5,4,3,6,5)。

步骤五，计算粒子的适应度值，并依适应度值由大到小为粒子排序，将适应度值最大的粒子作为备选优化结果。

通过适应度函数F(公式9)计算出该粒子的适应度值。

步骤六，判断粒子群算法是否达到最大迭代次数；如果粒子群算法达到了步骤一确定的最大迭代次数，则输出优化结果；否则，更新该粒子群的所有粒子的速度和位置信息，并重新进行步骤五和步骤六。

粒子的位置和速度信息的更新公式分别为：

其中，

为第e次迭代进行中第i个粒子在n维空间中的位置；

为第e+1次迭代进行中第i个粒子在n维空间中的位置；

r为速度的约束因子，通常设置为1；

为第e次迭代进行中第i个粒子在n维空间中的速度；

为第e+1次迭代进行中第i个粒子在n维空间中的速度；

θ为惯性权重，用来表示保持原来速度的系数；

c₁，c₂分别为粒子跟踪个体最优的权重系数和跟踪群体最优的权重系数，通常都设置为2；

h₁,h₂为[0,1]区间内均匀分布的随机数；

为第e次迭代进行中第i个粒子在n维空间中迄今为止搜索到的最优位置；

为第e次迭代进行中粒子群体在n维空间中迄今为止搜索到的最优位置；

例如：某一个7×7离散粒子X为

基于纵坐标的坐标赋予粒子位置，则X对应的连续位置(即粒子矩阵中每列元素1所在行的行号)为(2,2,5,4,3,6,5)。假设该粒子飞行速度V＝(-8,8,8,8,8,8,8)，则新的位置应该为Y＝X+V＝(-6,10,13,12,11,14,13)。再将每个元素对7取余数，使其回到先前的解空间进行搜索。即mod(-6,7)＝1，mod(10,7)＝3，……从而Y应为(1，3，6，5，4，0，6)，Y实际上应为(1，3，6，5，4，7，6)，再将其还原为如下所示的离散粒子形式即可。

基于粒子群算法的局限性，每次求解结果并不是唯一的，需进行多次求解，以得出较为合理的结果，避免获得的仅仅是局部最优解。比较多次求解结果，目标函数值越小，求解结果越好。轮灌组划分方案以轮灌组组数少，且各轮灌组组间灌水时间相差小为优。

实施例

本发明选用的实施例为新疆莫索湾灌区150团一营二支渠的灌区资料。

莫索湾灌区地处天山北坡经济开发区中心地带，灌区辖四个大型农牧团场。所选用的150团一营二支渠位于三支干渠上游，全长2600m，控制面积256.33hm²，控制一连和三连两个用水单位。该支渠设计流量Q为0.65m³/s，渠道水利用系数η为0.975，下有12条斗渠，即下级渠道数量N为12，灌水轮期T₀为144h，各条斗渠的资料信息见表1。

表1莫索湾灌区150团一营二支渠信息

步骤一，构建目标函数，并确定轮灌组数目范围、粒子群粒子总数和最大迭代次数。

权重系数ω₁＝0.4，ω₂＝0.6；

则目标函数为min(0.4Z+0.6ΔT)；

根据公式4确定轮灌组数目M的取值范围为1～11；逐次遍历轮灌组数目M的取值；

选定粒子群粒子总数为100；

最大迭代次数为200。

步骤二，产生初始粒子群。即产生100个M行N列的矩阵，即粒子群。

步骤三，判断各粒子是否满足约束条件，若不满足约束条件，则摈弃该粒子，重新生成新的粒子，直到满足约束条件的粒子数目达到100。

步骤四，记录粒子的位置信息，粒子对应的位置即为粒子矩阵中每列元素1所在行的行号。

步骤五，计算粒子的适应度，并依适应度值由大到小为粒子排序，将适应度值最大的粒子作为备选优化结果。

步骤六，判断粒子群算法迭代次数是否达到最大迭代次数200，若达到则输出优化结果，若未达到，则更新粒子的速度和位置信息，并返回步骤五重新计算粒子的适应度值。

基于粒子群算法的局限性，每次求解结果并不是唯一的，需进行多次求解，以得出较为合理的结果，避免获得的仅仅是局部最优解，同时验证并检验算法的鲁棒性。

求解三次的结果。目标函数值越小，求解结果越好。轮灌组划分方案以轮灌组组数少，且各轮灌组组间灌水时间相差小为优。

表2第一次求解结果

此时该模型求解结果对应的目标函数值的值为0.2784。

对应模型求解结果一的轮灌组划分情况如下表3。

表3轮灌组划分结果

根据轮灌组划分结果可得到对应的渠系配水过程图，如图2所示。

产生的轮灌组划分结果采用渠系平面布局概化图的方式呈现如图3a至图3d。

根据以上求解结果所示，12条下级渠道共被划分为4个轮灌组，同时也意味着配水支渠被划分为4个续灌组，并依据组内轮灌，组间续灌的规则进行灌水。同时考虑到为了充分减少灌溉过程中的渗漏损失和蒸发损失，并依据渠道用水均衡的要求，在同一个轮灌组内对下级斗渠采取从后向前的灌水顺序。分析各个轮灌组引水时间的差异，比较4个轮灌组所用时间，轮灌组4所需时间最长，轮灌组1所需时间最短，所以总的灌水时间即为轮灌组4所需时间130.44h，轮灌组引水时间最大相差值6.71h。配水支渠在长达123.73h(约占总配水时间130.44h的95％)内保持稳定，无需调节闸门，只需在最后的6.71h里调节两次闸门，最后关闭闸门即可。在灌水过程中，无需调节闸门的时间越长，需要调节闸门的时间越集中，能在最大程度上不影响输水渠道水流的平稳，同时提高了配水精度，节省了人工劳力，给控制水闸的技术人员提供了操作上的便利，提高了渠道管理水平和效率。经过优化后的灌水方式减少了灌水时间，合理划分了轮灌组，提高了灌水效率，提高了管理人员工作效率，最终提高了灌溉水资源利用价值。

表4模型第二次求解结果

此时该模型求解结果对应的目标函数值为的值为0.2535。

表5模型第三次求解结果

此时该模型求解结果对应的目标函数值为的值为0.2961。

若不对该渠系进行优化而是顺次编组，则轮灌组划分如下表6。

表6顺次编组轮灌结果

对未优化与优化的结果对比见下表7。

表7结果对比

minZ即为最少轮灌组组数，min△T为各轮灌组引水时间相差最大值的最小值，代表各轮灌组间引水时间差最小。

未经优化的min Z(最少总轮灌组数)＝4，min△T＝70.21，目标函数值为43.726。结果明显不如进行优化之后的结果。总引水时间为151.61h，远比优化之后的130.44h所需时间长，而且轮灌组引水持续时间均匀性相比优化结果差了很多。可以表明经过优化后的轮灌组划分使得灌水过程中配水流量均匀平稳，输水总时间短，输水损失小，配水质量高。