热敏电阻系统辨识

基于递推最小二乘法的热敏电阻测量温度的系统辨识、热敏电阻温度变化对其阻值的影响1. 先测量常温下热敏电阻R,调节输出电压给热敏电阻进行水浴加热，与此同时不断调节电阻箱的电阻值，使得指针总在零刻度线附近从25°C开始每隔5°C测 量一次Rt，直到85Co2. 原理图如图(a)所示：R 二 Ri RxR0若G中检流为0则B和D等势，故此时 R2 ，可在检流计的灵敏度范围内得到Rx的值6)匝理图x3. 实验数据：温 度 C2025303540455055606570758085阻 值Qi945i6i0i327ii30960829707602522450392340302269二、递推最小二乘法的原理与设计1、数学模型设时不变SISO动态过程的数学模型为1.1)A(z-i)z(k)二 B(z-i)u(k) + e(k)其中，u(k)为过程的输入量，z(k)为过程的输出量，e(k)是噪声，多项式A(z-1)和B(z-i)为：A(z-i)二 1 + a z-i + a z — + —f a z-na1 2 naB (z-i) = b z-i + b z -2 f f b za~nbi 2 n在本次辨识实验设计中，第一次选取n = n =3•即 abA(zt) = i + a zt + a z-2 + a z-3i 2 3B(z-i) =bz-ifb z-2 fb z-3i 2 3第二次选取n = n =i，即abA(z-i) =i f a z-iiB(z-i) =b z-ii将此模型写成最小其中，格k)=hT (k|^^e(k) hT (k)是可观测的数h(k) =[_z(k_1),_z(k_2),_z(k_3),u2、辨识原理利用数据序列{h(k)}0 =[ 1,a2,, b ,b11),u(k _小化准则皿J (0)考虑模型z(k)二b (k)0 + e(k)的辨识问题，其中 和估计参数，准则函数J (0)可写成二次型的形式1.2)) 是均值为零的随机噪声。

量；)都是可观测的数据J (0) = (z _ H 0 口 (z _ H 0)l l l l量模型输出与实际过程输出的接近情况极小化J(0)，求得参数0(14)作来衡1.5)0 = (HtH )_iHt zl l l l首先将1.5 式一次完成算法写成1.6)0 二(Ht H ) _1 Ht z 2 P(1 )Ht zl l l l l l=[工h(i)hT (i)]_1[]Lh(i)z(i)]i =1 i =1定义P _1(k)=》h(i)hT (i) = Ht H1.7)k kP _1 (k _ 1)=丈 h(i)hT (i) = Ht Hk _1 k _1hT (1)hT (1)H=khT (2)， Hk_1 =hT (2)hT (k)hT (k _ 1)由(1.7)式可得：其中：i=1P_1 (k) = t1 h(i)hT (i) + h(k)hT (k) = P_i(k _ 1) + h(k)d (k) (1.8)i=1令：zk_1 =[z⑴’Z⑵‘…z(k_ 1”1.9 )1.10)1.14)则：e(k -1)二(Ht H ) -1 Ht z 二 P(k -1)[艺h(i)z(i)]k -1 k-1 k-1 k-1i=1于是有 p -1(k - 1)e (k -1)=艺 h(i)z (i)i=1令 z = [ z (1), z (2), z (k )]tk利用(1.8)和(1.9)式，可得e(k) = (HtH )-1 Htz = P(k)[h h(i)z(i)]k k k ki=1=P(k )[P-1 (k — 1)e (k — 1) + h(k) z (k)]=P(k){[P-1 (k) — h(k)hT (k)]G(k — 1) + h(k)z(k)}=e(k -1) + p(k)h(k)[z(k) -hT (k)6(k -1)]引进增益矩阵K (k)，定义为K(k) = P(k)h(k ) (1.11)则(1.10)式写成e (k) = e (k -1)+k( k)[ z (k) - hT (k )e (k -1)] (1.12)进一步把(1.8)式写成P(k) = [P-1 (k -1) + h(k )hT (k )]-1 (1.13)为了避免矩阵求逆运算，利用矩阵反演公式可将(1.13)式演变成P(k) = [P-1 (k -1) + h(k )hT (k )]-1 = [I - K (k )h 工(k )]P(k -1)将(1.14)式代入(2.11)式，整理后有K (k) = P(k - 1)h(k )险(k )P(k - 1)h(k) +1]-1 (1.15)综合(1.12)、(1.13)、(1.14)式便得到最小二乘参数估计递推算法e (k) = e (k -1) + K (k)[ z (k) - hT (k )e (k -1)]P(k) = [P-1 (k -1) + h(k )hT (k )]-1 = [I - K (k )h 工(k )]P(k -1)K (k) = P(k - 1)h(k )[hT (k )P(k - 1)h(k) +1]-1利用上述公式即可求得参数的估计值e三、辨识结果1、当模型建立为三阶系统时程序运行结果为0. 001000-0. 1840-0. 6364-0. 6294-1. 1973-0. 11132. 2281-0. 7510-0.8183-0. 7692-0. 8258-0.82520. 001000-0. 2235-0.0103-0. 3401-0.4229-0.4762-0.4389-0.2081-0.2182-0. 2202-0. 2343-0. 23430. 001000-0. 2702-0.13190. 13930. 57110. 0953-1. 04130. 11010. 14830. 12670. 16420. 16400. 0010000. 00520. 17920.2597-0. 396159. 7080193.30211. 6427-2.76140. 0063-3.3864-3.36170. 0010000. 00450. 16390.2389-0. 3842-3. 0163-8. 9025-0. 04250.15500. 03330. 17780.17680. 0010000. 00380. 14870.2180-0. 3723-65. 7405-211.1070-1. 72783.07150.06043. 74203.7153000 -185.04462.4578-0.01100.9024-0. 9070-21. 0107-1.33710.0895-0. 06000. 0736-0. 0007000-224.5077-0.954132. 11830. 24350. 1260-0. 0782-0. 52600.04850.00960. 06370. 0003000-271.2221-0.5119-2.05583. 1011-0. 8331-11. 9228-1. 10580.3462-0. 14540. 2963-0. 0015000■4. 183433.56980.4491-2. 5256-151.73122. 2375-0.9915-2.6810-1.0023-537. 7521-0. 00730LI口3.486135. 54070.4571-2. 60876. 85011. 9515-0. 9952-4.6448-0. 78494.3326-0. 00570LI口2.788938.23710.4667-2. 7076175.55982. 2112-0. 9918-2. 7777-0. 980360. 9790-0. 00710= [a a a b b1 2 3 1 2b ] = [- 0.8252 - 0.234330.1640 - 3.36170.1768 3.7153]最小二乘各次递推参数估计值参数的误差收敛情况kS实际曲线与模拟曲线的对比2、三阶程序运行结果分析从辨识结果中可以看出此系统辨识的效果不好，误差较大，从实际曲线与模拟曲线的对 比中可以明确看出。

说明此辨识系统的模型不适合用于热敏电阻温度变化对阻值的影响这一 实验3、当模型建立为一阶系统时程序运行结果为0.001000-0.8511-0. 8569-0.8365-0. 84680. 0010000.0203-0. 23750.57810. 2022—00-852.08670. 0068-0.02380.012300019.2635-12. 7228-3.4335-0. 6502-0.8509-0. 8481-0. 8488-0. 8484-0. 8490-0. 8482-0. 84800.06740. 15150. 13310. 14400. 12790. 14640. 15100. 01：148-0.00320.1：11：108-0.1：11：1050.1：11：108-0.1：11：109-0. 0002-0.66661.2471-0. 12140.0818-0.11210. 14490.0313则Gk)= (a b) =( -0.8480 0.1510)最 小二乘 各次递 推参数 估计值050差误数参0050--1000 02468 10 12 1420001500* 1000500020 30 40 50 60 70 80 90u/ °C4、一阶程序运行结果分析 有上述运行结果可以看出一阶系统模型下辨识后的曲线和原曲线误差很小，说明其比三阶更 适合用来模拟该系统。

四、程序代码：clc;clear;close;u=[20 25 30 35 40 45 50 55 60 65 70 75 80 85];z=[1945 1610 1327 1130 960 829 707 602 522 450 392 340 302 269]; N=length(u);c0=[0.001,0.001]'; %直接给出被辨识参数的初始值，取一个充分小的实向量p0=10'7*eye(2,2); %初始状态P0也采用直接取方式，取一个充分大的实数单位矩阵E=0.000001; %相对误差E参考值取0.000001c=[c0,zeros(2,13)]; %被辨识参数矩阵的初始值及大小e=zeros(2,14); %相对误差的初始值及大小%开始递推运算for k=4:1:N;h1=[-z(k-1),u(k-1)]'; %求 h1x=h1'*p0*h1+1;x1=inv(x);k1=p0*h1*x1; %求 k1c1=c0+k1*[z(k)-h1'*c0]; %求 ce2=(c1-c0)./c0;%求参数的相对变化e(:,k)=e2; %把当前相对变化的列向量加入误差矩阵的最后一列c0=c1; %新获得的参数作为下一次递推的旧参数c(:,k)=c1; %把当前所辨识参数的 c1 列向量加入辨识参数矩阵的最后一列p1=p0-k1*k1'*[h1'*p0*h1+1]; %求 p1 值p0=p1; %把当前值给下次用if e2==E break; %若参数收敛满足要求，终止计算endendc %显示被辨识参数e%显示辨识结果的收敛情况%开始计算模拟曲线u(1)=20;y1(1)=1945;for k=2:14 y1(k)=0.1510*u(k-1)+0.8480*y1(k-1);y1(k-1)=y1(k); u(k-1)=u(k);enda=c(1,:); b=c(2,:); ea=e(1,:); eb=e(2,:); %分离参数 figure(1); %画第 1 个图形i=1:14; %横坐标从 1 到 15 plot(i,a,'r+-',i,b,'mx-.') %画出 a1，a2，b1，b2 的各次辨识结果 xlabel('k')ylabel('辨识参数')％标注纵轴变量title('最小二乘各次递推参数估计值')％图形标题figure(2); %画第 2 个图形 i=1:14; %横坐标从 1 到 15plot(i,ea,'r',i,eb,'g') %画出al, a2, bl, b2的各次辨识结果的收敛情况 xlabel('k') %标注横轴变量ylabel('参数误差')％标注纵轴变量tit le('参数的误差收敛情况')％图形标题figure(3);plo t(u,z,'r')%画出实际曲线hold onplo t(u,y1,'b')%画出模拟曲线legend('z','yl');xlabel('u/°C') %标注横轴变量ylabel('z/Q') %标注横轴变量title('s实际曲线与模拟曲线的对比')%图形标题。