数学建模算法与应用：预测算法（6）预测习题练习

一，水塔总水量以及流速预测问题

1.1、题目

1.2、建立模型

1.3、用MATLAB计算，将“-”替换为-1。

1.4、拟合法

二、预测产值问题

2.1、题目

2.2、建立模型

一，水塔总水量以及流速预测问题

1.1、题目

某地区用水管理机构需要对居民的用水速度（单位时间的用水量）和日总用水量进行估计。现有一居民区，其自来水是由一个圆柱形水塔提供，水塔高12.2m，塔的直径为17.4m。水塔是由水泵根据水塔中的水位自动加水。按照设计，当水塔中的水位降至最低水位，约8.2m时，水泵自动启动加水；当水泵升高到最高水位，约10.8m时，水泵停止工作。

表1给出的是28个时刻的数据，但由于水泵正向水塔供水，有四个时刻无法测到水位（-）。

水塔中水位原始数据
时刻（t）/h	0	0.92	1.84	2.95	3.87	4.98	5.90
水位/m	9.68	9.48	9.31	9.13	8.98	8.81	8.69
时刻（t）/h	7.01	7.93	8.97	9.98	10.92	10.95	12.03
水位/m	8.52	8.39	8.22	-	-	10.82	10.5
时刻（t）/h	12.95	13.88	14.98	15.9	16.83	17.93	19.04
水位/m	10.21	9.94	9.65	9.41	9.18	8.92	8.66
时刻（t）/h	19.96	20.84	22.01	22.96	23.88	24.99	25.91
水位/m	8.43	8.22	-	-	10.59	10.35	10.18

试建立数学模型，来估计居民的用水速度和日总用水量。

1.2、建立模型

解：（1）插值法。要估计在任意时刻（包括水泵灌水期间）t 居民的用水速度和日总用水量，分如下三步。

①水塔中水的体积的计算。计算水的流量，首先需要计算出水塔中水的体积，即：

$V=\frac{\pi }{4}D^{2}h$ ，式中：D为水塔的直径；h为水塔中水位高度。

②水塔中水流速度的估计。居民的用水速度就是水塔中的水流速度，水流速度应该是水塔中水的体积对时间的导数，但由于没有每一时刻水体积的具体数学表达式，只能用差商近似导数。

由于在两个时段，水泵向水塔供水，无法确定水位的高度，因此在计算水塔中水流速度时要分三段计算。第一时段为0~8.97h，第二时段为10.95~20.84h，第三时段为23.88~25.91h。

上面计算仅给出流速的离散值，如果需要得到流速的连续型曲线，需要做插值处理，这里可以使用三次样条插值。

③日总用水量的计算。日用水量是对水流速度做积分，其积分区间是[0,24]，可以采用数值积分的方法计算。

1.3、用MATLAB计算，将“-”替换为-1。

计算的MATLAB程序如下：

clc, clear, close all
a=load('bitter_tea_seeds.txt');
t0=a([1:2:end],:); t0=t0'; t0=t0(:); %提出时间数据，并展开成列向量
h0=a([2:2:end],:); h0=h0'; h0=h0(:); %提出高度数据，并展开成列向量
D=17.4;
V=pi/4*D^2*h0; %计算各时刻的体积
dv=gradient(V,t0); %计算各时刻的数值导数（导数近似值）
no1=find(h0==-1) %找出原始无效数据的地址
no2=[no1(1)-1:no1(2)+1,no1(3)-1:no1(4)+1] %找出导数数据的无效地址
t=t0; t(no2)=[]; %删除导数数据无效地址对应的时间
dv2=-dv; dv2(no2)=[]; %给出各时刻的流速
plot(t,dv2,'*'),title('流速图') %画出流速的散点图
pp=csape(t,dv2); %对流速进行插值
tt=0:0.1:t(end); %给出插值点
fdv=fnval(pp,tt); %计算各插值点的流速值
hold on, plot(tt,fdv) %画出插值曲线
I=trapz(tt(1:241),fdv(1:241)) %计算24小时内总流量的数值积分

结果展示


no1 =11       12       24       25       no2 =1 至 5 列10             11             12             13             23       6 至 8 列24             25             26       I =21221/17

1.4、拟合法

解：（2）拟合法。要估计在任意时刻（包括水泵灌水期间）t 居民的用水速度和日总用水量，分如下三步。

①水塔中水的体积的计算。计算水的流量，首先需要计算出水塔中水的体积，即：

$V=\frac{\pi }{4}D^{2}h$ ，式中：D为水塔的直径；h为水塔中水位高度。

上面的计算仅给出流速的离散值，流速的散点图如图所示中的 “ * ” 点。如果需要得到流速的连续型曲线，可以拟合多项式曲线，原始数据总共有28个观测值，其中4个无效数据。

图中总共有20个数据点，这里我们分三段进行三次多项式拟合，应用前6个数据点拟合三次多项式，即在时间区间[0，4.98]上拟合三次多项式；应用第6个数据点到第10个数据点，即在时间区间[4.98，12.03]，拟合第二个三次多项式；应用第10个数据点到第20个数据点，共11个数据点，即在时间区间[12.03，25.91]，拟合第三个三次多项式。

拟合得到的分段三次多项式曲线如上图所示。

③日总用水量的计算。日用水量是对水流速度做积分，其积分区间是[0,24]，可以采用数值积分的方法计算。

MATLAB程序如下：

clc, clear, close all
a=load('bitter_tea_seeds.txt');
t0=a([1:2:end],:); t0=t0'; t0=t0(:); %提出时间数据，并展开成列向量
h0=a([2:2:end],:); h0=h0'; h0=h0(:); %提出高度数据，并展开成列向量
D=17.4;
V=pi/4*D^2*h0; %计算各时刻的体积
dv=gradient(V,t0); %计算各时刻的数值导数（导数近似值）
no1=find(h0==-1) %找出原始无效数据的地址
no2=[no1(1)-1:no1(2)+1,no1(3)-1:no1(4)+1] %找出导数数据的无效地址
t=t0; t(no2)=[]; %删除导数数据无效地址对应的时间
dv2=-dv; dv2(no2)=[]; %给出各时刻的流速
hold on, plot(t,dv2,'*') %画出流速的散点图
a1=polyfit(t(1:6),dv2(1:6),3); %拟合第一个多项式的系数
a2=polyfit(t(6:10),dv2(6:10),3); %拟合第二个多项式的系数
a3=polyfit(t(10:20),dv2(10:20),3); %拟合第三个多项式的系数
dvf1=polyval(a1,[t(1):0.1:t(6)]); %计算第一个多项式的函数值
dvf2=polyval(a2,[t(6):0.1:t(10)]); %计算第二个多项式的函数值
dvf3=polyval(a3,[t(10):0.1:t(end)]); %计算第三个多项式的函数值
tt=t(1):0.1:t(end); dvf=[dvf1,dvf2,dvf3];
plot(tt,dvf) %画出拟合的三个分段多项式曲线
I=trapz(tt(1:241),dvf(1:241)) %计算24小时内总流量的数值积分

二、预测产值问题

2.1、题目

某大型企业1997---2000年的产值资料如表所示，试建立GM（1，1）预测模型，预测该企业2001---2005年的产值。

年份	1997	1998	1999	2000
产值/万元	27260	29547	32411	35388

2.2、建立模型（笔记打印，有曲面）

解：1.级比检验

建立原始序列数据 $x^{(0)}=(x^{(0)}(1),x^{(0)}(2),x^{(0)}(3),x^{(0)}(4))=(27260,29547,32411,35388)$

2.GM（1，1）建模

3.模型检验

模型的各种检验指标值的计算结果见表2。经检验，该模型的精度较高，可进行预测和预报。

年份	原始值	预测值	残差	相对误差	级比偏差
1997	27260	27260	0	0
1998	29547	29553.4421	-6.4421	0.0002	-0.0095
1999	32411	32336.4602	74.5398	0.0023	0.0025
2000	35388	35381.5524	6.4476	0.0002	-0.0022

4.预测值

2001---2005年预测值见表2。

年份	2001	2002	2003	2004	2005
产值/万元	38713.3978	42358.9998	46347.9045	50712.4404	55487.9803

MATLAB程序设计代码，如下：

clc,clear, format long g
x0=[27260  29547   32411   35388]'; %注意这里为列向量
n=length(x0);
lamda=x0(1:n-1)./x0(2:n)  %计算级比
range=minmax(lamda')  %计算级比的范围
theta=[exp(-2/(n+1)),exp(2/(n+1))]  % 计算级比的容许区间
x1=cumsum(x0)  %累加运算
B=[-0.5*(x1(1:n-1)+x1(2:n)),ones(n-1,1)];
Y=x0(2:n);
u=B\Y  %拟合参数u(1)=a,u(2)=b
syms x(t)
x=dsolve(diff(x)+u(1)*x==u(2),x(0)==x0(1)); %求微分方程的符号解
xt=vpa(x,6) %以小数格式显示微分方程的解
yuce1=subs(x,t,[0:n+4]); %求已知数据和未来5期的预测值
yuce1=double(yuce1); %符号数转换成数值类型，否则无法作差分运算
yuce=[x0(1),diff(yuce1)]  %差分运算，还原数据
epsilon=x0'-yuce(1:n)    %计算已知数据预测的残差
delta=abs(epsilon./x0')  %计算相对误差
rho=1-(1-0.5*u(1))/(1+0.5*u(1))*lamda'  %计算级比偏差值，u(1)=a
yhat=yuce(n+1:end)  %提取未来5期的预测值

结果展示：


lamda =0.9225978948793450.9116349387553610.915875438001582range =0.911634938755361         0.922597894879345theta =0.670320046035639          1.49182469764127x1 =272605680789218124606u =-0.089995172386872625790.2838424515xt =313834.0*exp(0.0899952*t) - 286574.0yuce =1 至 3 列27260          29553.4420565748          32336.46018497974 至 6 列35381.5523516033          38713.3978069354          42358.99982184127 至 9 列46347.9045382397          50712.4404287318          55487.9803059021epsilon =1 至 3 列0         -6.44205657481507          74.53981502031094 列6.44764839668642delta =1 至 3 列0      0.000218027433404917       0.002299830767958754 列0.000182198722637233rho =-0.00953940978346313       0.00245664647930999      -0.00218345852198909yhat =1 至 3 列38713.3978069354          42358.9998218412          46347.90453823974 至 5 列50712.4404287318          55487.9803059021