软测量作业1-4题答案

工业过程建模与软测量课程作业1

作业1：假设已知Q i 、R 1和R 2、A 1和A 2，试根据物料平衡方程建立机理模型，求液位h 1和h 2。

解：假设忽略两个储罐的蒸发量，根据物料平衡方程，可列出如下的微分方程 1

1

12A i dh Q Q dt =- （1）

2212o dh

A Q Q dt =- （2）

根据流体运动方程可得：

12

121

h h Q R -=

（3）

2

2

o h Q R =

（4）

于是，串接液位储罐过程机理模型如下：

1

12112122212

12

1212

2--i o dh h h A Q dt R dh h h h A dt R R h h Q R h Q R ?=??

-?=-???

-?=??

?=??

（5）

作业2：

MATLAB 回归分析程序如下

x=[0.5;0.7;1.0;2.0;4.0;5.0;7.0;8.0;10.0;12.0;15.0]; %输入数据

X=[ones(11,1),x];

Y=[0.139;0.196;0.272;0.545;1.033;1.272;1.671;1.807;1.991;2.088;2.121];

plot(x,Y); %输入图像

[b,bint,r,rint,stats]=regress (Y,X,0.02); %回归分析 b,bint,stats

rcoplot(r,rint); %作残差图 z=b(1)+b(2)*x; %预测及作图 plot(x,Y,'k+',x,z,'r'); grid on;

结果如下 b = 0.2855 0.1533 bint =

-0.0731 0.6440 0.1059 0.2007 stats =

0.9023 83.1657 0.0000 0.0685

从MATLAB 计算结果中，我们可以得出，0β =0.2855，1β =0.1533，0β的置信区间为[-0.0731，0.6440]，1β的置信区间为[0.1059,0.2007],2

0.9023r

=，0.0000p =。

0.02p <，则回归模型0.28550.1533y x =+成立。

图一为原始数据的拟合曲线，从曲线中可以看出，x 与Y 之间存在线性关系。

图一 x 与Y 拟合曲线

图二为残差图，从残差图中可以卡出，除最后一个数据外，其余数据的残杀离零

点均较近，且残差的置信区间均包含零点，这说明回归模型0.28550.1533y x =+能较好的符合原始数据，而最后一个数据可视为异常点。

图二 残差图

图三为拟合直线与散点图

图三 回归直线

作业3：

（1）建模方法如下

采用多元线性回归分析来建立模型，根据所给出的数据，系统有两个输出变量1y （PH 值）和2y （亚硫酸浓度），五个输入变量1x （甲基丙稀磺酸钠浓度）、2x （氯酸钠浓度）

、3x （焦亚硫酸钠浓度）、4x （β羟基乙硫醇浓度）、5x （聚合釜反应温度），且共有34组测量数据。在这里，我们选用前32组数据进行回归分析，最后两组数据进行验证。在建模过程中研究12y y 、与(5)m m =个自变量

12345x x x x x 、、、、之间的线性相关关系，12y y 、为随机变量，12345x x x x x 、、、、为一般变量，i x 相互之间不相关，i ε为随机误差。

假定因变量12y y 、与自变量12345x x x x x 、、、、线性相关。以1y 为例，收集到

(32)n n =组数据125(,,,,)(1,2,,32)i i i y x x x i = 满足以下关系：

011223344552

()0()cov(,)0(..1,2,,)i i j y x x x x x E i Var i j i j n ββββββεεσεε=++++++??=?

?=?

?=≠=?

111532,132,51(1),1n x x C x x x ??

??==??????

01112

2213324325,,y y y y ββεβεβεββεβ??

??

????????????????===??????????????????

??????

则所建回归模型的矩阵形式为：

2

~(0,1)

n n y C N βε

εσ=+??? 以上为经典多远回归模型，其中y 是可观测的随机向量，ε是不可观测的随机向量，C 是已知矩阵，,βσ是未知参数，并且n m >，且()1rank C m =+。 根据最小二乘可得β的最佳线性无偏估计为：

-11?=[C'C]C y β

- 由此解得多元线性回归方程的回归系数β，构建多元线性回归模型。写成标量形式：

01155

???y x x βββ=+++ （2）仿真与验证

MATLAB 回归分析程序如下：

x=[0.050 0.169 0.450 0.077 56.0;0.055 0.169 0.450 0.077 56.0; 0.065 0.169 0.450 0.077 56.0;0.070 0.169 0.450 0.077 56.0; 0.080 0.169 0.450 0.077 56.0;0.085 0.169 0.450 0.077 56.0; 0.070 0.150 0.450 0.077 56.0;0.070 0.155 0.450 0.077 56.0; 0.070 0.160 0.450 0.077 56.0;0.070 0.165 0.450 0.077 56.0;

0.070 0.170 0.450 0.077 56.0;0.070 0.175 0.450 0.077 56.0;

0.070 0.180 0.450 0.077 56.0;0.070 0.169 0.435 0.077 56.0;

0.070 0.169 0.440 0.077 56.0;0.070 0.169 0.445 0.077 56.0;

0.070 0.169 0.450 0.077 56.0;0.070 0.169 0.455 0.077 56.0;

0.070 0.169 0.460 0.077 56.0;0.070 0.169 0.465 0.077 56.0;

0.070 0.169 0.470 0.077 56.0;0.070 0.169 0.450 0.060 56.0;

0.070 0.169 0.450 0.065 56.0;0.070 0.169 0.450 0.070 56.0;

0.070 0.169 0.450 0.075 56.0;0.070 0.169 0.450 0.080 56.0;

0.070 0.169 0.450 0.085 56.0;0.070 0.169 0.450 0.090 56.0;

0.070 0.169 0.450 0.077 50.0;0.070 0.169 0.450 0.077 52.0;

0.070 0.169 0.450 0.077 54.0;0.070 0.169 0.450 0.077 56.0];

%输入数据

X=[ones(32,1),x];

y1=[2.0;2.0;2.0;2.0;2.0;1.95;2.0;2.0;2.0;2.0;2.0;2.0;2.0;1.9;1.9;1.9;

2.0;

2.0;1.9;1.9;1.9;2.2;2.0;2.0;2.0;1.9;1.85;1.85;1.9;1.9;2.0;2.0];

y2=[0.00154;0.00144;0.00133;0.00123;0.00133;0.00143;0.00297;0.00277;

0.00266;0.00143;0.00113;0.00113;0.00113;0.00267;0.00215;0.00144;

0.00123;0.00133;0.00195;0.00205;0.00215;0.00164;0.00154;0.00154;

0.00143;0.00102;0.00092;0.00062;0.00133;0.00133;0.00123;0.00123]; [b1,bint1,r1,rint1,stats1]=regress(y1,X); %回归分析及检验

b1,bint1,stats1

[b2,bint2,r2,rint2,stats2]=regress(y2,X);

b2,bint2,stats2

rcoplot(r1,rint1); %作残差图

rcoplot(r2,rint2);

结果如下

b1 =

2.7783

-1.2997

-1.1891

-0.7367

-10.0725

0.0105

bint1 =

1.0917 4.4650

-4.4686 1.8691

-4.7070 2.3289

-3.7272 2.2538

-13.7423 -6.4027

-0.0033 0.0242

stats1 =

0.5823 7.2495 0.0002 0.0023

b2 = 0.0140 -0.0017 -0.0768 0.0009 -0.0280 0.0000 bint2 =

-0.0008 0.0288 -0.0295 0.0260 -0.1076 -0.0459 -0.0253 0.0271 -0.0602 0.0041 -0.0001 0.0002 stats2 =

0.5372 6.0369 0.0008 0.0000

从MATLAB 的运算结果可得出

112345

2

123452.7783 1.2997 1.18910.736710.07250.01050.01400.00170.07680.00090.02800.0000y x x x x x y x x x x x =----+??

=--+-+? 图四为1y 的残差图，图五为2y 的残差图，从图中可以看出，回归模型能较好的符合测量数据。

验证如下：

1123452123452.7783 1.2997 1.18910.736710.07250.01050.01400.00170.07680.00090.02800.0000y x x x x x y x x x x x =----+??

=--+-+? 得

00166

.02 03027.2100166.02 00927.21-==-==y y y y

作业4：

MATLAB程序如下：

f11=fopen('sampledata.txt','r');

[r01,c]=fread(f11,[100,4],'float');

fclose(f11);

X=r01(:,1:3);

Y=r01(:,4);

pz=[X,Y];

% function y=pls(pz);

[row,col]=size(pz);

[row,col]=size(pz);

aver=mean(pz);

stdcov=std(pz); %求均值和标准差

rr=corrcoef(pz); %求相关系数矩阵

data=zscore(pz); %数据标准化

stdarr=(pz-aver(ones(row,1),:))./stdcov(ones(row,1),:); %标准化数据结果与zscore()一致

x0=pz(:,1:col-1);

y0=pz(:,end); %提取原始的自变量、因变量数据

e0=stdarr(:,1:col-1);

f0=stdarr(:,end); %提取标准化后的自变量、因变量数据

num=size(e0,1); %求样本点的个数

temp=eye(col-1); %对角阵

%以下计算w,w*和t的得分向量

for i=1:col-1

p=e0'*f0;

q=norm(p);

w(:,i)=p/q;

t(:,i)=e0*w(:,i); %计算成分ti的得分

alpha(:,i)=e0'*t(:,i)/(t(:,i)'*t(:,i)) %计算alpha_i,其中

(t(:,i)'*t(:,i))等价于norm(t(:,i))^2

e=e0-t(:,i)*alpha(:,i)' %计算残差矩阵

e0=e;

%计算w*矩阵

if i==1

w_star(:,i)=w(:,i);

else

for j=1:i-1

temp=temp*(eye(col-1)-w(:,j)*alpha(:,j)');

end

w_star(:,i)=temp*w(:,i);

end

%以下计算ss(i)的值

beta=[t(:,1:i),ones(num,1)]\f0 %求回归方程的系数

beta(end,:)=[]; %删除回归分析的常数项

cancha=f0-t(:,1:i)*beta; %求残差矩阵

ss(i)=sum(sum(cancha.^2)); %求误差平方和

%以下计算press(i)

for j=1:num

t1=t(:,1:i);

f1=f0;

she_t=t1(j,:);

she_f=f1(j,:); %把舍去的第j个样本点保存起来

t1(j,:)=[];

f1(j,:)=[]; %删除第j个观测值

beta1=[t1,ones(num-1,1)]\f1; %求回归分析的系数

beta1(end,:)=[]; %删除回归分析的常数项

cancha=she_f-she_t*beta1; %求残差向量

press_i(j)=sum(cancha.^2);

end

press(i)=sum(press_i)

if i>1

Q_h2(i)=1-press(i)/ss(i-1)

else

Q_h2(1)=1

end

if Q_h2(i)<0.0985

fprintf('提出的成分个数r=%d',i);

r=i;

break

end

end

beta_z=[t,ones(num,1)]\f0; %求标准化Y关于主成分得分向量t的回归系数

beta_z(end,:)=[]; %删除常数项

xishu=w_star*beta_z; %求标准化Y关于X的回归系数，且是针对标准数据的回归系数，每一列是一个回归方程

mu_x=aver(1:col-1);

mu_y=aver(end);

sig_x=stdcov(1:col-1);

sig_y=stdcov(end);

ch0=mu_y-mu_x./sig_x*sig_y*xishu; %计算原始数据的回归方程的常数项

xish=xishu'./sig_x*sig_y; %计算原始数据的回归方程的系数，每一列是一

个回归方程

Rc=corrcoef(x0*xish'+ch0,y0)

sol=[ch0;xish'] %显示回归方程的系数，每一列是一个方程，每一列的第一个数是常数

y=sol;

结果如下：

sol =

1.0e+05 *

-0.0000

1.2931

-0.6463

-0.6463

因此，偏最小二乘回归方程如下：

x

y-

-

=

x

1

3

64634

2

64634

129310x