2.6.1 嵌入法（embedding method）
前面介绍的线性规划模型，可用如图2.16所示的示意图表示。为与后面嵌入法后的LP模型以示区别，这里线性规划模型用LP＇表示。
图2.16 线性规划模型结构
【例2.17】二单元潜水含水层稳定流水量分配。
有一矩形潜水含水层，三面为隔水边界，河水位视为定水头，地下水接受大气降水补给，向河流排泄，已知含水层的几何参数和导水系数T、降水入渗补给强度N。规划要在两个单元抽水，单位抽水费用为c1和c2。规划要求对二单元的稳定抽水量最优分配，使其满足总需水量D和二单元最低水位限制条件下，总的费用最小（雅·贝尔，1985；陈爱光等，1991；鲍新华等，2010），如图2.17所示。
图2.17 二单元潜水含水层规划问题
有关参数为
L=10000m,B=10000m
补给强度：N=0.36m/a
导水系数：T=3.6×106m2/a
水位限制（假设河水位为零）：h1min=2.5m，h2min=5m
需水量：D=4.5×107m3/a
抽水费用：c1、c2
解法一：不考虑水量均衡问题。
可给出该问题的线性规划模型为
给定c1=0.2元/m3、c2=0.1元/m3，代入相关数据整理后，采用前面介绍的Excel求解，结果如图2.18所示。
将c1、c2值交换下，又能求解出新的一组解。两个方案下的结果汇总见表2.6。
图2.18 二单元潜水含水层解法一Excel求解结果
表2.6 [例2.17]解法一二方案结果汇总表
解法二：考虑水量均衡问题。
设I 为水力梯度，根据达西定律，单元1与河流、单元1与单元2断面的流量应为TBI。取河水位为0，可列出单元1和单元2的均衡方程如下：
上述两式整理后作为约束条件，加入到模型LP＇中，有新的LP模型。
用Excel求解上面LP模型，解见表2.7。
表2.7 [例2.17]解法二二方案结果汇总表
两个线性规划模型均有解，也都符合最优决策总是趋向于抽水费用最小的单元抽水。哪个解法对呢？
如果对解法一从单元上及整个区域上，对解进行水均衡验证计算，就能找到问题答案了。
解法一错误：水量不均衡（以方案1为例）。
没有数值模型的优化模型，可能有解，也可能无解。但即使是有解，解在单元上和均衡区上均没达不均衡。不均衡怎么可能是稳定运动！
解法二则把每个单元的均衡方程作为约束条件加入（嵌入）到规划模型中，解一定满足每个单元上的水均衡。每个单元上均衡了，整个区域上肯定也是均衡的。
解法二加入的两个方程，实际上是反映了地下水内在运动规律的水位流量关系，实际上就是地下水的数值方程。
本例题中，单元上的入渗补给量与开采量之差，为含水层向河流的排泄量。按此思路，可以计算傍河水源地袭多河流的补给量，或者计算这部分水量占开采量的比例。
嵌入法：将地下水数值方程组作为一组约束条件加入线性规划模型中的建模方法称为嵌入法。
数值方程形式
[A]·H+[I]Q=R
式中[A]——线性方程组系数矩阵；
H——各节点各时段水位或降深的列向量；
[I]——单位矩阵；
Q——可控输入量；
R——常数项。
参见上面的LP模型，嵌入法可用图2.19表示。
图2.19 嵌入前后线性规划模型的变化
【例2.18】承压水二维稳定问题。
二维承压均质含水层，四周为定水头边界，设含水层隔水顶板底面为水位基准线，其剖面及平面差分网格如图2.20所示。要求规划布置4口井，其中3眼分别位于6、7、11节点上，另1眼井可布置在10、14、15、18、19任一节点上，每井抽水量不得小于P（P单位为强度，m/d），总需水量为4P。并要求各井抽水产生的降深值之和最小（即水位值之和最大），已知含水层导水系数T=1000m2/d，抽水前天然水位为水平，水位h=20m，差分剖分网格Δx=Δy=1000m（陈爱光等，1991，结果有更正）。
图2.20 承压含水层平面问题示意图
描述该问题的数学模型为
式中 Q——抽水强度；
Lx、Ly——含水层沿x、y方向的长度。
设决策变量Qi为稳定抽水量，状态变量hi为水位值。依题意有：
目标函数 max Z=h6+h7+h10+h11+h14+h15+h18+h19
可以用均衡法写出各节点的差分方程。节点均衡区是以节点为中心，四周各取单元间一半距离构成的。本例题共8个计算节点，需要列出每个点的均衡方程。如节点6的方程为
代入Δx=Δy，整理有
式中 h02、h05表示边界节点上对应点已知水位，本题中河流边界初始水位均为20m。
类似可列出其他节点的方程，可用矩阵表示整个8个节点的方程：
其他约束为
将以上数据代入Excel规划求解，设定3 种不同的P 值（P=0，不抽水；P=0.01m3/d；P=0.03965m3/d，抽水稳定时某井降深刚好达到20m，即水位到达承压水隔水顶板。P=0.03965m3/d数值是试算出来的），运行结果见表2.8。
表2.8 [例2.18]3种不同抽水强度下Excel优化结果
定性分析看，当点6、7、11布置3口抽水井后，为使区域总体水位降深小，另1口抽水井应该布置在离开6、7、11点尽可能远的地方。从图2.20分析看，18点最佳。从计算结果看，在抽水情况下的后两个方案，优化的抽水井或在18点，或在19点，基本与定性分析一致。
使用Surfer软件，输入坐标点及水位高程，可以方便地画出不同优化方案下的水位等值线图。[例2.18]最后一种情况下（表2.8）等值线图如图2.21所示。从图2.21中可以看到，6、7、11井形成了一个水位下降漏斗连片区，19 井局部形成单井水位下降漏斗区。
【例2.19】潜水含水层三单元二时段非稳定抽水最优分配问题。
设有一矩形潜水含水层（长3L、宽L，L=1000m），西侧以河流作定水头边界，年平均水位H0=40m，其他3 边界为隔水边界（图2.22），年平均降水入渗强度N=154mm/a（或0.4mm/d），含水层平均导水系数T=1000m2/d，给水度μ=0.3。该含水层分为3个正方形单元网格，各单元均可抽水，目前该含水层尚未开发，现要求在满足今后第一年和第二年分别提供水量为P1=5000m3/d和P2=7500m3/d的条件下，如何合理开采地下水，使得第二年末3个单元平均水位之和最大。同时，要求在第二年末单元1的平均水位不低于H1=39m（杨悦所等，1992；卢文喜，1999；鲍新华等，2010）。
(https://www.daowen.com)
图2.21 抽水强度P=0.03965m3/d时流场等水位线图（图及绘图数据）
图2.22 三单元潜水非稳定问题
解：以1年为1个管理时段（Δt=365d），共分两个时段，决策变量为各单元各时段抽水量
，状态变量为水位
1，2，3为单元序号；k=1，2为时段）。
依题意有：
（1）为建立数值约束方程，先要确定未开发（未抽水）前天然水位。
抽水前稳定情况下三单元的均衡方程为
整理式（2.17）有
代入基本数据（N 取0.4mm/d），有方程：
用Excel规划求解，解为
（2）再确定非稳定均衡方程（隐式法）。
1）第一时段末。
2）第二时段末。类似有
将已知数据L=1000m、N =0.4mm/d、T=1000m2/d、H0=40m、Δt=365d、
代入，整理有
将式（2.22）嵌入式（2.16）中，构成新的约束，Excel求解结果如图2.23所示。
图2.23 [例2.19]Excel求解结果
图2.24 二单元潜水非稳定问题
【例2.20】潜水含水层二单元三时段非稳定抽水最优分配
为进一步熟悉非稳定抽水下嵌入法求解过程，对[例2.19]稍加改造：设有一矩形潜水含水层（长2L、宽L，L=1000m），西侧以河流作定水头边界，年平均水位H0=40m，其他3边为隔水边界（图2.24），年平均降水入渗强度N=0.4mm/d，含水层平均导水系数T=1000m2/d，给水度μ=0.3。该含水层分为2个正方形单元网格，各单元均可抽水，目前该含水层尚未开发，现要求在满足今后第一年、第二年和第三年分别提供水量为P1=2000m3/d、P1=3000m3/d 和P3=5000m3/d的条件下，如何合理开采地下水，使得第三年末3个单元平均水位之和最大。同时，要求在第三年末单元1的平均水位不低于H1=39m。
解：以1年为1个管理时段（Δt=365d），共分3个时段，决策变量为各单元各时段抽水量
，状态变量为水位
（i=1，2为单元序号；k=1，2，3为时段）。
以题意，有
数值约束：
（1）先确定初始水位。抽水前二单元的均衡方程为
整理式（2.5）有
代入基本数据（N 取0.4mm/d），有方程：
用Excel规划求解，解为
（2）再确定非稳定均衡方程。
1）第一时段末。
2）第二时段末。类似有
3）第三时段末。类似有
将式（2.31）嵌入式（2.24）中，构成新的约束，Excel求解结果见表2.9。
表2.9 【例2.20】Excel规划求解结果
【例2.21】潜水含水层四单元二时段非稳定抽水最优分配（[例2.19]、[例2.20]稍变动）。
设有一矩形潜水含水层（长4L、宽L，L=1000m），西侧以河流作定水头边界，年平均水位H0=40m，其他3边为隔水边界（图2.25），年平均降水入渗强度N=154mm/a（或平均0.4mm/d），含水层平均导水系数T=1000m2/d，给水度μ=0.3。该合水层分为4个正方形单元网格，各单元均可抽水，目前该含水层尚未开发，现要求在满足今后第一年和第二年分别提供水量为P1=5000m3/d和P2=7500m3/d的条件下，如何合理开采地下水，使得第二年末3个单元平均水位之和最大。同时，要求在第二年末单元1的平均水位不低于H1=39m。（该例题可以简略讲，大部分安排学生完成）
图2.25 四单元潜水非稳定问题
解：以1年为1个管理时段（Δt=365d），共分两个时段，决策变量为各单元各时段抽水量
，状态变量为水位
（i=1，2，3，4为单元序号；k=1，2为时段）。
依题意有
数值约束：
（1）先确定初始水位。抽水前稳定三单元的均衡方程为
整理式（2.34）有
代入基本数据（N 取0.4mm/d），有方程：
用Excel规划求解，解为
。
（2）再确定非稳定均衡方程。
1）第一时段末。
2）第二时段末。类似有
将已知数据L=1000m、N=0.4mm/d、T=1000m2/d、H0=40m、Δt=365d、
代入，整理有
观察式（2.39）可发现，系数矩阵每行零元素很多，特别是节点方程越多，这种现象越突出。这正是地下水数值约束方程系数矩阵的一般特点。一般来说，无论是采用有限差或有限元方法建立的数值方程，一般都是高阶、稀疏、对角优势（一般是绝对对角优势）矩阵。且节点单元适当编号时，非零元素会尽可能靠近对角线所在的一个比较窄的带状范围内，该范围之外，都是零元素。如果不考虑水量变量的话，对有限差分法来说，每行的非零元素（水位或降深变量）最多是5个[式（2.39）中前一个系数矩阵]。对有限元来说，每行的非零元素（水位或降深变量）数一般在个位数以内（10以内）。其中的原因大家可以对节点方程的建立过程稍作研究，就会明白的。
将式（2.39）嵌入式（2.33）中，构成新的约束，Excel求解结果见表2.10。
表2.10 [例2.21]Excel规划求解结果
上面采用嵌入法建立地下水资源管理模型的方法比较直观，且同时可求解水位和流量两类变量。但对实际问题来说，嵌入的数值方程规模一般较大，规划模型嵌入数值方程后，规模急剧膨胀。但不嵌入又不满足地下水运动内在规律。有没有一种方法，不嵌入这一庞大的数值方程，还能保留地下水内在运动规律呢？ 这就是下面要介绍的响应矩阵法。
