黄河小浪底调水调沙问题

黄河小浪底调水调沙问题

摘要

为了确定排沙量与时间、排沙量与水流量的函数关系，我们可以用SAS软件做线性回归得到排沙量与时间的函数关系式，再利用所求函数在区间[0,24]上进行积分得到总排沙量1.93962亿吨。对于排沙量与水流量之间的关系，按时间分为两段进行拟合，最终用MATHLAB软件来画出图像，确定排沙量与排水量之间的函数关系式。

关键词：调水调沙实验，排沙量，matlab，拟合

1

目 录

摘要 ............................... 1

问题的提出 .......................... 3

问题分析............................ 3

模型假设............................ 3

模型的建立与求解 .................... 4

结论以及分析检验 ................... 10

讨论与推广 ......................... 10

参考文献........................... 11

2

一．问题的提出

2004年6月至7月黄河进行了第三次调水调沙试验，特别是首次由小浪底、三门峡和万家寨三大水库联合调度，采用接力式防洪预泄放水，形成人造洪峰进行调沙试验获得成功．整个试验期为20多天，小浪底从6月19日开始预泄放水，直到7月13日恢复正常供水结束．小浪底水利工程按设计拦沙量为75.5亿立方米，在这之前，小浪底共积泥沙达14.15亿吨．这次调水调试验一个重要目的就是由小浪底上游的三门峡和万家寨水库泄洪，在小浪底形成人造洪峰，冲刷小浪底库区沉积的泥沙．在小浪底水库开闸泄洪以后，从6月27日开始三门峡水库和万家寨水库陆续开闸放水，人造洪峰于29日先后到达小浪底，7月3日达到最大流量2700立方米/每秒，使小浪底水库的排沙量也不断地增加．下面是由小浪底观测站从6月29日到7月10日检测到的试验数据：

表1: 试验观测数据 单位：水流为立方米 / 秒，含沙量为公斤 / 立方米 日期 6.29 6.30 7.1 7.2 时间 8:020:8:020:8:020:8:020:0 00 0 00 0 00 0 00 水流180190210220230240250260量 0 0 0 0 0 0 0 0 含沙32 60 75 85 90 98 100 102 量 日期 7.5 7.6 7.7 7.8 时间 8:020:8:020:8:020:8:020:0 00 0 00 0 00 0 00 水流260250230220200185182180量 0 0 0 0 0 0 0 0 含沙118 120 118 105 80 60 50 30 量 现在，根据试验数据建立数学模型研究下面的问题： (1) 给出估算任意时刻的排沙量及总排沙量的方法； (2) 确定排沙量与水流量的变化关系。

7.3 7.4 8:020:8:020:0 00 0 00 2652702722650 0 0 0 108 112 115 116 7.9 8:020:0 00 1751500 0 26 20 7.10 8:020:0 00 100900 0 8 5 二．问题分析

1、对于问题一，所给数据中水流量x和含沙量h的乘积即为该时刻的排沙量y即：y=hx。

2、对于问题二，研究排沙量与排水量的关系，从实验数据中可以看出，开始排沙量随水量增加而增加，而后随水流量的增加而减少，显然变化关系并非线性的关系，为此，把问题分为两部分，从水流量增加到最大值为第一阶段，从水流量最大值到结束为第二阶段，分别来研究水流量与排沙量之间的函数关系。

三．模型假设

1、水流量和排沙量都是连续的，不考虑上游泄洪所带来的含沙量和外界带来的含沙量。

2、时间是连续变化的，所取时间点依次为1,2,3，…,24,单位时间为12h.

四．模型的建立与求解

3

已知给定的观测时刻是等间距的，以6月29日零时刻开始计时，则各次观测时刻分别为ti3600(12i4),i1,2,...,24,

其中：计时单位为s。第一次观察的时刻t128800，最后一次观察的时刻

t241022400。记第i（i1,2,...,24）次观测时水流量为vi，含沙量为ci，则第

i次观测时的排沙量为yicivi。相关数据见表2 日期 表2 插值数据对应关系 单位：排沙量为kg/s 6.29 6.30 7.1 7.2 7.3 8:00 20:00 3 4 8:00 20:00 5 6 8:00 20:00 7 8 7.4 时8:00 20:0间 0 节1 2 点 8:00 20:08:00 20:00 0 9 10 11 12 排 11401575187020702352250026522862302431283074沙57600 00 00 00 00 00 00 00 00 00 00 00 量 日7.5 7.6 7.7 7.8 7.9 7.10 期 时8:00 20:08:00 20:08:00 20:08:00 20:08:00 20:08:00 20:0间 0 0 0 0 0 0 节13 14 15 16 17 18 19 20 21 22 23 24 点 排306830002714231016001110910000455030008000 4500 沙00 00 00 00 00 00 0 0 0 0 量 对于问题（1），根据所给问题的实验数据，要计算任意时刻的排沙量，就要确定出排沙量随时间变化的规律，可以通过插值来实现。考虑到实际中的排沙量应该是时间的连续函数，为了提高模型的精度，采用三次样条函数进行插值。

利用MATLAB函数，求出三次样条函数，得到排沙量y=y（t）与时间的关系，

然后进行积分，就可以得到总得排沙量zy(t)dt。最后求得总得排沙量为

t1t21.844109t，计算的MATLAB程序如下：

排沙量对时间的关系图像的MATLAB程序 t=1:1:24;

y=[57600,114000 ,157500,187000 ,207000,235200,250000,265200,286200,302400,312800,307400,306800,300000,271400,231000,160000,111000,91000,000,45500,30000,8000,4500]; plot(t,p,'b')

4

3.5x 10532.521.510.500510152025

对于问题（2），研究排沙量与水流量的关系，从实验数据可以看出，开始排沙量是随着水流量的增加而增长，而后是随着水流量的减少而减少。显然，变化规律并非是线性的关系，为此，把问题分为两部分，从开始水流量增加到最大值2720m3/s（即增长的过程）为第一阶段，从水流量的最大值到结束为第二阶段，分别来研究水流量与排水量的关系。画出排沙量与水流量的散点图 两个阶段的数据如表3表4所示 表3：第一阶段试验数据 序号 1 2 3 4 5 6 7 水流量 1800 1900 2100 2200 2300 2400 2500 含沙量 32 60 75 85 90 98 100 表4：第二阶段的试验观测数据 序号 1 2 3 4 5 6 7 水流量 2650 2600 2500 2300 2200 2000 1850 含沙量 116 118 120 118 105 80 60 排沙量与时间的关系图像的MATLAB程序： >>t=1:1:24;

y=[57600,114000,157500,187000,207000,235200,250000,265200,2862000,2400,312800,307400,306800,300000,271400,231000,160000,111000,91000,000,45500,30000,8000,4500]; plot(t,y,'r')

8 2600 102 8 1820 50 5

3x 1062.521.510.500510152025

第一阶段的排沙量与水流量之间的关系MATLAB程序：

>>x=[1800,1900,2100,2200,2300,2400,2500,2600,2650,2700,2720]; h=[32,60,75,85,90,98,100,102,108,112,115];

x1=[2650,2600,2500,2300,2200,2000,1850,1820,1800,1750,1500,1000,900]; h1=[116,118,120,118,105,80,60,50,40,32,20,8,5]; plot(x,h,'r:')

6

1201101009080706050403018001900200021002200230024002500260027002800

第一阶段三次多项式拟合函数以及拟合效果程序与结果：

三次多项式拟合由MATLAB拟合函数求解出a0=a1=0,a2=0.0032,a3=-2.4929.则拟合函数

h=0.0032x^2-2.4929x^3,

>> A1=polyfit(x,h,3) z1=polyval(A1,x);

plot(x,h,'k+',x,z1,'r')

>Warning: Polynomial is badly conditioned. Add points with distinct X values, reduce the degree of the polynomial, or try centering and scaling as described in HELP POLYFIT. > In polyfit at 76 A1 =

1.0e+003 *

0.0000 -0.0000 0.0032 -2.4929

7

1201101009080706050403018001900200021002200230024002500260027002800

第一阶段四次多项式拟合函数以及拟合效果程序与结果：

四次多项式拟合由MATLAB拟合函数求解出a0=a1=a2=0,a3=0.0121,a4= -7.4347则拟合函数

h=0.0121x^3-7.4347x^4 >>A2=polyfit(x,h,4) z2=polyval(A2,x);

plot(x,h,'*',x,z2,'r')

>Warning: Polynomial is badly conditioned. Add points with distinct X values, reduce the degree of the polynomial, or try centering and scaling as described in HELP POLYFIT. > In polyfit at 80 A2 =

-0.0000 0.0000 -0.0000 0.0121 -7.4347

8

1201101009080706050403018001900200021002200230024002500260027002800

第二阶段三次多项式拟合函数以及拟合效果程序与结果：

三次多项式拟合由MATLAB拟合函数求解出a0=a1=0,a2=-0.9475,a3= 4.9601.则拟合函数

h=-0.9475x^2+4.9601x^3 >>A3=polyfit(x1,h1,3) z3=polyval(A3,x1);

plot(x,h,'*',x1,z3,'b')

>Warning: Polynomial is badly conditioned. Add points with distinct X values, reduce the degree of the polynomial, or try centering and scaling as described in HELP POLYFIT. > In polyfit at 80 A3 =

-0.0000 0.0006 -0.9475 4.9601

9

1401201008060402008001000120014001600180020002200240026002800

第二阶段四次多项式拟合函数以及拟合效果程序与结果：

四次多项式拟合由MATLAB拟合函数求解出a0=a1=0，a2=-0.0013，a3= 1.1219 a4=-3.5952则拟合函数

h=-0.0013x^2+1.1219x^3-3.5952x^4 >>A4=polyfit(x1,h1,4) z4=polyval(A4,x1);

plot(x1,h1,'k*',x1,z4,'r:')

>Warning: Polynomial is badly conditioned. Add points with distinct X values, reduce the degree of the polynomial, or try centering and scaling as described in HELP POLYFIT. > In polyfit at 80 A4 =

-0.0000 0.0000 -0.0013 1.1219 -3.5952

五．结论以及分析检验

对于第一阶段三次多项式拟合由MATLAB拟合函数求解

a0=a1=0,a2=0.0032,a3=-2.4929.则拟合函数h=0.0032x^2-2.4929x^3 对于第一阶段四次多项式拟合由MATLAB拟合函数求解出

a0=a1=a2=0,a3=0.0121,a4= -7.4347则拟合函数h=0.0121x^3-7.4347x^4 对于第二阶段三次多项式拟合由MATLAB拟合函数求解出

a0=a1=0,a2=-0.9475,a3= 4.9601.则拟合函数h=-0.9475x^2+4.9601x^3, 对于第二阶段四次多项式拟合由MATLAB拟合函数求解出a0=a1=0，a2=-0.0013，a3= 1.1219 a4=-3.5952则拟合函数h=-0.0013x^2+1.1219x^3-3.5952x^4

六．讨论与推广

10

对于第二个问题，由于MATLAB软件的计算可能有些偏差导致拟合的函数关系可能与实际有稍微偏差，此外，还可以进行高次的拟合。

七．参考文献

[1]. [2]. [3].

司守奎，《数学建模算法与应用》，国防工业出版社，2012. 姜启源，《数学模型》（第四版），北京：高等教育出版社，2011。 韩中庚，《数学建模方法及其应用》（第二版），北京：高等教育出版,2009

11