热电偶实验报告1
热电偶辨识实验
一 实验目的：
通过对热电偶的辨识，并对辨识结果进行动态误差修正，掌握系统辨识方法中的时域辨识方法和对测量结果的动态误差修正方法，了解动态误差修正在实际生活中的应用。
二 实验器材：
热电偶一个，应变放大器一台，桥盒一个，数采模块，PC机一台。
三 实验原理：
本实验是基于热电偶测温的工作原理所做，即：热电偶是由两种不同成份的材质导体组成闭合回路， 当两端存在温度梯度时，回路中就会有电流通过，此时两端之间就存在电动势——热电动势，这就是所谓的塞贝克效应(Seebeck effect）。两种不同成份的均质导体为热电极，温度较高的一端为工作端，温度较低的一端为自由端，自由端通常处于某个恒定的温度下。读出热端的电动势，然后根据热电动势与温度的函数关系可得出当前的温度值。当我们将热电偶放入热水中，由于温度的变化，产生一个阶跃信号，通过图形确定系统是几阶系统，然后对模型进行辨识，并对测量结果进行动态误差修正，将修正前后的响应特性曲线进行比较，对实验结果进行分析。
四 实验过程：
（1）将热电偶通过桥盒与应变放大器相连， 然后与PC机连接好，组成一个完整的传感器系统。按如图1所示方式将热电偶的两个接线端接入桥盒。
图1 热电偶与桥盒的连接
（2）PCI6013——AI接线分配如图2所示，我们这里选择的是第一通道，所以连接33号跟64号线。
图2 PCI6013——AI接线分配
（3）打开labview，单击启动采集按钮，将K型热电偶迅速放进热水瓶中，待输出稳定后保存数据然后取出热电偶冷却，然后重复多次试验，保存数据。
（4）利用所保存的数据进行系统辨识和误差修正。
五 实验数据分析
下面通过实验来进行系统辨识及其动态误差修正。
它利用不平衡电桥产生的热电势来补偿热电偶因冷端温度的变化而引起热电势的变化,经过设计,可使电桥的不平衡电压等于因冷端温度变化引起的热电势变化而实现的自动补偿。后接放大器来将热电偶输出的电压信号进行放大，经过数采卡进行数据采集，最后传到计算机处理。经实验测得几组数据，经过仿真得到如图3所示阶跃响应特性曲线。
图3 阶跃响应特性曲线
由于初始状态不为零，所以通过一定的计算使他的初始状态变为零。如图4所示。
图4 未加补偿前的阶跃响应图像
本实验采用动态测量系统的时域辨识，根据上述所获得的阶跃响应曲线可判断该系统为一阶或者二阶系统。本实验采用最小二乘法来辨识差分方程。
为了辨识差分方程的阶次，分别采用一阶和二阶进行辨识，通过比较拟合误差平方和来确定是几阶系统。调用[a,b,J]=DE_Id01(x,y,1)和[a,b,J]=DE_Id01(x,y,2)可得一阶差分方程的拟合误差为J =0.0016，二阶差分方程的拟合误差为J =0.0010，由于拟合误差差距不大且一阶系统辨识较为简单，故将待辨识差分方程的阶次设为1阶。
设离散传递函数为H1（z）=(b0+b1z-1)/(1+a1z-1)其中 a1, b1和b0为待辨识方程的系统参数。得到参数辨识结果为
Transfer function: 1.997e-005 - 1.649e-005 z^-1 ----------------------------
1 - 0.9993 z^-1 Sampling time: 0.001
这是一个一阶低通系统，其幅频特性和相频特性如图5
所示。
图5 原系统频响特性
由幅频特性曲线可见，原系统的工作频带是在低频段，公称增益应为增益下降到10%附近的公称增益的上限频率约为v10%=0.05。
因此，应将动态补偿合成系统HeD（Z
）设计成归一化频率为vc=v10%= 0.5的低通滤波器。于是，取vc=0.05，按2阶Butterworth低通滤波器设计HeD(Z),可得
Transfer function: 0.07296 + 0.07296 z^-1 ---------------------- 1 - 0.8541 z^-1 Sampling time: 0.001
该系统的极值点的位置如图6所示
图6 系统极值点
由此可知该补偿系统所有极值点均在单位圆内，该系统是一个稳定的可实现系统。
图7为系统动态特性补偿之前与之后的阶跃信号响应特性曲线对照示意图。
图7 补偿前后的阶跃信号响应对照图
图8 补偿综合系统的频响特性
由此可以看出补偿之后的动态特性明显得以改善。但是实验的结果仍有一定的误差。
误差产生的情况主要由以下几个方面构成：
（1）在本次实验中由于要做多次实验所以水的温度会不断降低，还有热电偶初始温度也因为试验次数的不同而不同，所以系统产生的阶跃信号不一样。
（2）外界环境（噪声、温度）的影响。
六 总结
通过这次大作业，对系统辨识及系统的动态误差修正有了更多的认识，也认识到了自己在很多方面的不足，最突出的问题就是不能够熟练的运用所学的知识来解决实际的问题，也不能够灵活的考虑问题。这就是以后的学习中所重点解决的问题。
附录
t=1:2000; L=length(t); for k=1:L if (k
x0(k)=18; end end
y0=importdata('100,2000.txt'); x=x0;
y=y0-0.01; plot(t,y);
[a,b,J]=DE_Id01(x,y,1); Ts=0.01;
H1=filt(b,a,Ts);
G1=d2c(H1,'tustin'); bode(G1);
title('幅频特性与相频特性'); [b,a]=butter(2,0.05/50); H3=filt(b,a,Ts); G3=d2c(H3); bode(G3);
H2=tf(1,[1 0],Ts); H4=H3/H1*H2;
[b,a]=tfdata(H4,'v') H4=filt(b,a,Ts); [z,p,k]=tf2zp(b,a); zplane(b,a);
title('数据零极点显示'); y1=filter(b,a,y);
plot(t,x,t,y*100,t,y1)
function [a,b,J]=DE_Id01(x,y,N) L=length(x); LN=L-N-1; N1=2*N+1; for m=1:N, m2=2*m; m1=m2-1;
for n=1:N, n2=2*n; n1=n2-1; c(m1,n1)=0; c(m1,n2)=0; c(m2,n1)=0; c(m2,n2)=0; for k=1:LN,
c(m1,n1)=c(m1,n1)+x(k+m)*x(k+n); c(m1,n2)=c(m1,n2)-x(k+m)*y(k+n); c(m2,n1)=c(m2,n1)+y(k+m)*x(k+n); c(m2,n2)=c(m2,n2)-y(k+m)*y(k+n); end end
c(m1,N1)=0; c(m2,N1)=0; d(m1)=0; d(m2)=0; for k=1:LN,
c(m1,N1)=c(m1,N1)+x(k+m)*x(k+N+1); d(m1)=d(m1)+x(k+m)*y(k+N+1);
c(m2,N1)=c(m2,N1)+y(k+m)*x(k+N+1); d(m2)=d(m2)+y(k+m)*y(k+N+1); end end
for n=1:N, n2=2*n; n1=n2-1; c(N1,n1)=0; c(N1,n2)=0; for k=1:LN,
c(N1,n1)=c(N1,n1)+x(k+N+1)*x(k+n); c(N1,n2)=c(N1,n2)-x(k+N+1)*y(k+n); end end
c(N1,N1)=0; d(N1)=0; for k=1:LN, c(N1,N1)=c(N1,N1)+x(k+N+1)^2; d(N1)=d(N1)+x(k+N+1)*y(k+N+1); end
bx=c\d'; %% =inv(c)*d' ÏßÐÔ·½³Ì×éÇó½â a(1)=1.0; b(1)=bx(N1); %% ²ÎÊý·ÖÀë for n=2:N+1, n1=2*(N-n+2);n2=n1-1; a(n)=bx(n1); b(n)=bx(n2); end
for k=1:LN, k1=k+N; e(k)=b(1)*x(k1)-y(k1); %% ¼ÆËãÄ£ÐÍÔ¤±¨Îó²î¼°Æäƽ·½ºÍ
for n=1:N, e(k)=e(k)+b(1+n)*x(k1-n)-a(1+n)*y(k1-n);end end
J=sum(e.^2); return