水流量估计课件.doc

对估计水塔流量的探究

摘要
本文是一个关于投资收益和风险权衡规划问题的解答，问题中给出四种（或多种）投资方案供投资者选择，每种方案都有相应的收益率和风险率，我们的目的是

关键字：

一．问题重述
某居民区有一供居民用水的圆柱形水塔，一般可以通过测量其水位来估计水的流量，但面临的困难是，当水塔水位下降到设定的最低水位时，水泵自动启动向水塔供水，到设定的最高水位时停止供水，这段时间无法测量水塔的水位和水泵的供水量．通常水泵每天供水一两次，每次约两小时.
水塔是一个高12.2m，直径17.4m的正圆柱．按照设计，水塔水位降至约8.2m时，水泵自动启动，水位升到约10.8m时水泵停止工作．
表1 是某一天的水位测量记录，试估计任何时刻（包括水泵正供水时）从水塔流出的水流量，及一天的总用水量．

时刻(h)
水位(cm)
0 0.92 1.84 2.95 3.87 4.98 5.90 7.01 7.93 8.97
968 948 931 913 898 881 869 852 839 822
时刻(h)
水位(cm)
9.98 10.92 10.95 12.03 12.95 13.88 14.98 15.90 16.83 17.93
// // 1082 1050 1021 994 965 941 918 892
时刻(h)
水位(cm)
19.04 19.96 20.84 22.01 22.96 23.88 24.99 25.91
866 843 822 // // 1059 1035 1018
表1 水位测量记录（符号//表示水泵启动）

二．模型的假设
1) 假设储水器的空气压力强稳定；
2) 假设题目涉及的水泵功率恒定；
3) 假设一天中温度不变或者影响不大；
4) 居民的用水量没有突变点。

三． 符号说明
表1未完善
符号
符号说明

赤纬角

方位角

太阳高度角

太阳时角(地球每个小时自转15°称之为时角)
W
地理纬度
J
经度
L
影长
D
日期
H
杆子高度
n
某点的地方时间
L
杆影长度

四．模型的建立与求解
4.1问题一
4.1.1问题1分析
题目要求描建立每个时刻水塔流出用水流量变化的数学模型，分析水箱水位关于各个时刻的变化规律。首先通过查阅相关数学建模文献，找寻建模中可能遇到的概念、方法，再挖掘建立模型相关的方法，为之后的建模奠定坚实的基础。模型建立第一步,理清自变量（时间），因变量（流量）。先预处理数据，将水位转化为水流量，再处理流量—时间散点图，分段多项式拟合构造对应函数，使之符合题意，同时根据相应的标准进行统计、分析和构建相关的数学模型。

4.1.2数据的预处理
将数据导进MATLAB中并将水塔水位转化为水塔中体积，公式如下
①
得出的体积与时间关系如下
表二未完善
时刻(h)
体积(cm)
0 0.92 1.84 2.95 3.87 4.98 5.90 7.01 7.93 8.97

时刻(h)
体积(cm)
9.98 10.92 10.95 12.03 12.95 13.88 14.98 15.90 16.83 17.93

时刻(h)
体积cm)
19.05 19.96 20.84 22.01 22.96 23.88 24.99 25.91

4.1.2模型1的建立
通过对题目的探讨，我们先利用MATLAB作出水体积与时间的散点图。
体积—时间散点图

图一
我们发现，图中散点可以分为五个阶段：未供水第一段，供水第一段，未供水第二阶段，供水第二阶段，未供水第三阶段。
由于供水阶段没有数据，所以先根据未供水段的散点图，利用插值拟合出未供水段的体积关于时间函数方程。根据流体力学知识可得流量是指单位时间内流经封闭管道或明渠有效截面的流体量，既是
②
其中Q是流量，S为截面面积,v为水流速度，L为体积，t为时间；
③
根据以上的方程式可以知道，先要拟合出流量—时间方程，再对方程求一次微分（导数），就可以得出流量—时间方程。

4.1.3步骤一：拟合方式的选择

由于供水阶段没有数据，所以先根据未供水段的散点图，利用插值拟合出未供水段的函数方程。
我们优先选择多项式拟合的方式，下面对多项式原理进行建模分析：
假设给定数据点 (i=0,1,…,m)，为所有次数不超过的多项式构成的函数类，现求一,使得
; ④
当拟合函数为多项式时，称为多项式拟合，满足式（1）的称为最小二乘拟合多项式。特别地，当n=1时，称为线性拟合或直线拟合。
显然
; ⑤
为的多元函数，因此上述问题即为求的极值 问题。由多元函数求极值的必要条件，得
; ⑥
即
 ; ⑦
（3）是关于的线性方程组，用矩阵表示为
; ⑧
式（3）或式（4）称为正规方程组或法方程组。
可以证明，方程组（4）的系数矩阵是一个对称正定矩阵，故存在唯一解。从式（4）中解出 (k=0,1,…，n)，从而可得多项式
; ⑨
可以证明，式（5）中的满足式（1），即为所求的拟合多项式。我们把称为最小二乘拟合多项式的平方误差，记作
; ⑩
由式(2)可得
⑪
多项式拟合的一般方法可归纳为以下几步：
(1) 由已知数据画出函数粗略的图形——散点图，确定拟合多项式的次数n；
(2) 列表计算和；⑫
(3) 写出正规方程组，求出；
(4) 写出拟合多项式。 ⑬
由于要求拟合函数比较圆滑，所以多项式的最高次幂在三到六之间比较合适;

4.1.4步骤二：利用MATLAB引用拟合工具箱对未供水段进行三次多项式拟合
将未供水段一二数据利用三次多项式拟合（polyfit();polyval();）过所有散点拟合出两个时间段的曲线图如下：

图二
拟合出曲线后就可以利用polysym2()语句直接得出曲线方程如下

；①
L（t）求导之后可得，既是为未供水段流量方程；
②

图三

4.1.5步骤三：拟合出供水段的流量—时间曲线

由于在步骤二已经求出未供水的流量方程以及流量数据，我们采用拟合预测的方法，利用已有的未供水段数据，适当取点，做时间序列或者多项式拟合预测。
在未供水第一段取三个采样点，以及第二段取三个取样点，将它们的数据（t,:）中的t结合成一个一个6X1矩阵，（：Q）同理；再利用多项式拟合出相应过这六点的拟合曲线，优化取点，再写出供水段的方程。拟合出的曲线如下:

图四
可见，我们拟合出来的曲线不够圆滑，证明拟合预测具有相应的误差值，我们会在以下的步骤进行修正；利用poly2sym()语句得出拟合所得的方程。
；③
4.1.6步骤四：差分法求最后一段未供水数据流量方程
原始数据可知未供水第三段只有三个数据，直接进行拟合会产生很大的误差，所以不能采用数据拟合的方式得出流量方程。
三个散点数据点 ，因此我们采用差分法求流量方程。
差分法是微分方程的一种近似数值解法。具体地讲，差分法就是把微分用有限差分代替，把导数用有限差商代替，从而把基本方程和边界条件（一般均为微分方程）近似地改用差分方程（代数方程）来表示，把求解微分方程的问题改换成为求解代数方程的问题。
对三个点的横坐标进行两两相减，纵坐标同理，再求微分，微分之后再求商，就是最后所求的流量值。利用MATLAB进行操作后得出的流量数值为
；
去除偏离值得到的流量值为
；
在未供水第二段取两个个采样点，以及未供水第三段取两个个取样点，将它们的数据（t,:）中的t结合成一个一个4X1矩阵，（：Q）同理；再利用多项式拟合出相应过这六点的拟合曲线，优化取点，再拟合写出供水段的方程。
曲线如下：

图五
拟合出的方程为：
.
以上得出的图形和方程既是所求的流量曲线图，根据方程和曲线可以求解任何时候从水塔中流出的流量！
综上所述：
流量方程为：

4.1.7模型的建立与检验
模型建立之后，我们需要检验它的准确程度，因为拟合数据具有相当大的数据误差，所以我们需要对模型进行检验以及修正，根据已有数据，我们采用拟合度分析—— 拟合优度检验。
根据已有数据，我们代进拟合方程求出系列数据，原始数据的水位值与拟合曲线求出的水位值进行regress拟合度分析如下：
我们得出它们的残差分析图：
我们发现在11到12时刻，有一个很大的偏离值，所以在这个位置需要优化这个位置既是之前提到的：流量函数的斜率突变点。最后得出的拟合度为R=？？？

4.2问题二，居民用水量的求解
根据题一建立的模型以及得出的流量—曲线方程，我们可以通过曲线对时间求定积分就可以得出居民的用水量；
未供水第一段：

供水第一段：

未供水第二段：

供水第二段：

未供水第三段：

居民总用水：

利用MATLAB求得：居民一天用水量为:65682m3;

4.2.1问题二的分析

4.2.2模型的建立与检验：

五、结论

六、模型的改进与推广

参考文献：
[1] 司守奎，孙玺菁:数学建模算法与应用[S]: 国防工业出版社,2013.
[2]何银涛，张梅黄华.垂直单轴跟踪光伏支架方阵间距计算界面设计.太阳能报刊.2015.
[3]郑鹏飞,林大钧,刘小羊,吴志庭.基于影子轨迹线反求采光效果的技术研究[].广东:华南理工大学20023 7.
[4] 林根石,利用太阳视坐标的计算进行物高测量与定位[A].南京大学学报,1991.
[5] 吴济廉,影端轨迹周年变化的实践与分析.地理教学报，2013,50(5):9-20.
[6]张丽梅,赵建立,乔立山.用二值矩阵表示研究格矩阵的{1}-广义逆和{1,2}-广义逆.中国知网.
未修改

附录：
%%
x=[0,0.92,1.84,2.95,3.87,4.98,5.90,7.01,7.93,8.97,9.98,10.92,10.95,12.03,12.95,13.88,14.98,15.90,16.83,17.93,19.04,19.96,20.84,22.01,22.96,23.88,24.99,25.91];
y_y=[968,948,931,913,898,881,869,852,839,822,0,0,1082,1050,1021,994,965,941,918,892,866,843,822,0,0,1059,1035,1018];
y=pi.*(17.4.^2./4)*y_y;
xlswrite('高程表.xlsx',y);
%导入数据，时间以及水位
figure(1)
plot(x,y,'*');%做出水位与时间的关系图
% text(7.93,852,' \leftarrow 水位一阶段');
% text(pi/2,0,' \leftarrow 水位二阶段');
grid on;
set(gca,'xtick',[0:2:28]);
title('体积—时间曲线');
xlabel('时间/h');
ylabel('水体积/m^3');

nihezhiduoxiangshi=polyfit(x(1:10),y(1:10),3);
x_x=0:0.1:9;
y_y_2=polyval(nihezhiduoxiangshi,x_x);
nihezhiduoxiangshi=polyfit(x(13:23),y(13:23),3);
x_x1=10.95:0.1:20.84;
y_y_3=polyval(nihezhiduoxiangshi,x_x1);
figure(4)
plot(x_x,y_y_2);
hold on;
plot(x_x1,y_y_3);
set(gca, 'XLim',[0,30]);
set(gca, 'YLim',[0,3*1.0e+05]);
grid on;
set(gca,'xtick',[0:2:28]);
title('体积—时间拟合曲线');
xlabel('时间/h');
ylabel('水体积/m^3');

poly2sym(nihezhiduoxiangshi);
poly2sym(nihezhiduoxiangshi_3);
%%
%拟合供水一段流量
temp_1=polyval(weifen,[8,9]);
temp_2=polyval(weifen_2,[11,12]);
temp_12=[temp_1,temp_2];
x_4=[8,9,11,12];
b_4=polyfit(x_4,temp_12,3)
x_5=[9:11];
y_4=abs(polyval(b_4,x_5));
length(y_4);
plot(x_5,y_4,'RED');
hold on;
%%差分法求流量
%拟合供水第二段
dt3=diff(x(26:28));% 最后3个时刻的两两之差 >>
dh3=diff(y(26:28));% 最后3个水位的两两之差 >>
dht3=-dh3./dt3;
t3=[20 20.8 x(26) x(27)];
xx3=[abs(polyval(weifen_2,t3(1:2))),dht3];
c3=polyfit(t3,xx3,3)
tp3=20.8:0.1:24;
x3=polyval(c3,tp3);
plot(tp3,x3,'RED');
%积分求用水量
y_M=0.1*trapz(abs(weifen))
y_M1=0.1*trapz(abs(weifen_2))
y_M2=0.1*trapz(abs(y_4)) %gongshui1
y_M3=0.1*trapz(abs(c3)) %gongshui2
y_M+y_M1+y_M2+y_M3

14