数模培训数据拟合方法

1、数模培训数据拟合方法 设设 R=at+ba,b为待定系数为待定系数求电阻求电阻R R随温度随温度t t的变化规律的变化规律。热敏电阻数据：热敏电阻数据：t(0C) 20.5 32.7 51.0 73.0 95.7电阻电阻R( ) 765 826 873 942 1032引例引例1 1：热敏电阻电阻值的变化规律：热敏电阻电阻值的变化规律 t (h) 0.25 0.5 1 1.5 2 3 4 6 8c ( g/ml) 19.21 18.15 15.36 14.10 12.89 9.32 7.45 5.24 3.01 对某人用快速静脉注射方式一次性注射某种药物对某人用快速静脉注射方式一次性注射某种药

2、物300mg后，后，经过时间经过时间t采集血样，测得血药浓度采集血样，测得血药浓度c如下表：如下表：求血药浓度随时间的变化规律求血药浓度随时间的变化规律c(t).半对数坐标系半对数坐标系(semilogy)(semilogy)下的图形下的图形为待定系数kcectckt,)(0Log10c(t)=a t + b引例引例2 2：血药浓度的变化规律：血药浓度的变化规律数据拟合问题的提法数据拟合问题的提法数据拟合问题：一维二维，数据拟合问题：一维二维，数据，即平面上的数据，即平面上的n n个点个点(xi(xi，yi)yi)，i=1,2,ni=1,2,n，xixi互不相同，寻求一个函数互不相同，寻求一个

3、函数( (曲线曲线)y=f(x)y=f(x)，使使f(x)f(x)在某种准那么下与所有数据点最为接近，即曲线拟合的在某种准那么下与所有数据点最为接近，即曲线拟合的最好，如以下图所示最好，如以下图所示( (图中图中ii为为(xi(xi，yi)yi)与与y=f(x)y=f(x)的距离的距离) )。 Oxyi(xi，yi)数据拟合问题的求解思路数据拟合问题的求解思路线性最小二乘法是解决数据拟合最常用的方法。线性最小二乘法是解决数据拟合最常用的方法。 根本思路：根本思路： 令令 f(x)=a1r1(x)+a2r2(x)+amrm(x) (1) 其中其中rk(x)是事先选定的一组函数，是事先选定的一组函

4、数，ak是待定系数是待定系数 (k=1,2,m，m0), 称排除速率称排除速率.b) 中心室血液容积为常数中心室血液容积为常数 V, t=0 瞬时注入药物的瞬时注入药物的剂量为剂量为 d, 血药浓度立即为血药浓度立即为.dV由假设由假设 a), 中心室的血药浓度中心室的血药浓度 c(t)应满足微分方程应满足微分方程dckcdt 由假设由假设 b), 方程的初始条件为方程的初始条件为: 0.dcV 求解得求解得: .ktdc teV 即血药浓度即血药浓度c(t)按指数规律下降按指数规律下降.2. 给药方案设计给药方案设计简单实用的给药方案是简单实用的给药方案是:每隔一定时间每隔一定时间 , 重复

5、注入固定剂量重复注入固定剂量 D, 使血药浓使血药浓度度 c(t) 呈周期性变化呈周期性变化, 并保持在并保持在 c1-c2 之间之间. x0yc1c2 2. 给药方案设计给药方案设计简单实用的给药方案是简单实用的给药方案是:每隔一定时间每隔一定时间 , 重复注入固定剂量重复注入固定剂量 D, 使血药浓使血药浓度度 c(t) 呈周期性变化呈周期性变化, 并保持在并保持在 c1-c2 之间之间. 为此为此, 初次剂量需加大到初次剂量需加大到 D0.由式由式 得到得到: ktdc teV 2022111,lncDVcDV cckc 显然显然, 当当 c1, c2 给定后给定后, 要确定给药方案要确

6、定给药方案 0,DD 必须知道参数必须知道参数 V 和和 k.2. 由实验数据作曲线拟合以确定参数由实验数据作曲线拟合以确定参数lnlndcktV问题化为由数据问题化为由数据 ti , yi ( i=1,8 ) 拟合直线拟合直线12ya ta记记为了用线性最小二乘法拟合为了用线性最小二乘法拟合 的系数的系数 V 和和 k, 先取对数得先取对数得 ktdc teV 12ln ,lndyc akaV 用用Matlab作线性最小二乘法拟合作线性最小二乘法拟合, 得到得到120.2347,2.9943.aa lnlndcktV问题化为由数据问题化为由数据 ti , yi ( i=1,8 ) 拟合直线拟

7、合直线12ya ta记记为了用线性最小二乘法拟合为了用线性最小二乘法拟合 的系数的系数 V 和和 k, 先取对数得先取对数得 ktdc teV 12ln ,lndyc akaV 用用Matlab作线性最小二乘法拟合作线性最小二乘法拟合, 得到得到120.2347,2.9943.aa 由实验数据由实验数据 d=300 (mg) 算出算出:0.2347,15.02.kV拟合曲线为拟合曲线为: 0.234730015.02tc te 3. 结论结论将将 k, V 和给出的和给出的 c1=10, c2=25 代入代入lnlndcktV得得:D0=375.5, D=225.3, 3.9. 给药方案不妨定

8、为给药方案不妨定为:D0=375 mg , D=225 mg , 4 小时小时.范例：薄膜渗透率的测定范例：薄膜渗透率的测定 一、问题：一、问题： 某种医用薄膜，具有从高浓度的溶液向低浓度某种医用薄膜，具有从高浓度的溶液向低浓度的溶液扩散的功能，在试制时需测定薄膜被物质分的溶液扩散的功能，在试制时需测定薄膜被物质分子穿透的能力。子穿透的能力。 测定方法：用面积为测定方法：用面积为S S的薄膜将容器分成体积分的薄膜将容器分成体积分别为别为 的两部份，在两局部中分别注满该物质的两部份，在两局部中分别注满该物质的两种不同浓度的溶液。此时该物质分子就会从高浓的两种不同浓度的溶液。此时该物质分子就会从高

9、浓度溶液穿过薄膜向低浓度溶液中扩散。平均每单位时度溶液穿过薄膜向低浓度溶液中扩散。平均每单位时间通过单位面积薄膜的物质分子量与膜两侧溶液的浓间通过单位面积薄膜的物质分子量与膜两侧溶液的浓度差成正比，比例系数度差成正比，比例系数K K表征了薄膜被该物质分子穿表征了薄膜被该物质分子穿透的能力，称为渗透率。定时测量容器中薄膜某一侧透的能力，称为渗透率。定时测量容器中薄膜某一侧的溶液浓度，以此确定的溶液浓度，以此确定K K。BAVV 、VAVBS二、问题分析二、问题分析 考察时段考察时段tt，t+tt+t薄膜两侧容器中该物质质量薄膜两侧容器中该物质质量的变化。的变化。 jt(s)10020030040

10、05006007008009001000jC4.544.995.355.655.906.106.266.396.506.59设设 ，对容器的，对容器的B B部分溶液部分溶液浓度的测试结果如下表：（浓度单位浓度的测试结果如下表：（浓度单位 ）23BAcm10S,cm1000VV3cm/mg 1 在容器的一侧，物质质量的增加是由于另一侧的物质向该侧渗透的结果，因此物质质量的增量应等于另一侧的该物质向这侧的渗透量。) t (CV) tt (CVAAAA 以容器以容器A A侧为例，在时段侧为例，在时段tt，t+tt+t物质质量的增量物质质量的增量为：为：) t (C) t (CBA、分别表示在时刻分别

11、表示在时刻t t膜两侧溶液膜两侧溶液设设的浓度，浓度单位的浓度，浓度单位: :3cm/mg 由于平均每单位时间通过单位面积薄膜的物由于平均每单位时间通过单位面积薄膜的物质分子量与膜两侧溶液的浓度差成正比，比例系质分子量与膜两侧溶液的浓度差成正比，比例系数为数为K K。 因此，在时段因此，在时段tt，t+tt+t，从，从B B侧渗透至侧渗透至A A侧侧的该物质的质量为：的该物质的质量为：tS)CC(KAB于是有：于是有：) t (CV) tt (CVAAAAtS)CC(KAB两边除以两边除以tt，并令，并令t0t0取极限再稍加整理即得：取极限再稍加整理即得：)CC(VSKdtdCABAABA 、

12、分别表示在初始时刻两侧溶液的浓度分别表示在初始时刻两侧溶液的浓度其中其中12) 注意到整个容器的溶液中含有该物质的质量不变注意到整个容器的溶液中含有该物质的质量不变, ,与初与初始时刻该物质的含量相同，因此始时刻该物质的含量相同，因此 BBAABBAAVV) t (CV) t (CV从而：从而：加上初值条件：加上初值条件：.)0(CBB)VV(SKC)V1V1(SKdtdC) t (CVVVV) t (CABBABBABBABBABAA代入式（代入式（1）得：）得：便可得出便可得出CB(t)的变化规律，从而根据实验数据进行的变化规律，从而根据实验数据进行拟合，估计出参数拟合，估计出参数K, 。

13、BA 、三、数学模型三、数学模型假设：假设：1 1薄膜两侧的溶液始终是均匀的；薄膜两侧的溶液始终是均匀的；2 2平均每单位时间通过单位面积薄膜的物质分平均每单位时间通过单位面积薄膜的物质分 子量与膜两侧溶液的浓度差成正比。子量与膜两侧溶液的浓度差成正比。3 3薄膜是双向同性的即物质从膜的任何一侧向薄膜是双向同性的即物质从膜的任何一侧向 另一侧渗透的性能是相同的。另一侧渗透的性能是相同的。基于假设和前面的分析，基于假设和前面的分析，B B侧的浓度侧的浓度CB(t)CB(t)应满足如下应满足如下微分方程和初始条件：微分方程和初始条件：BBABBABBAB)0(C)VV(SKC)V1V1(SKdtd

14、C四、求解方法：四、求解方法：1. 函数拟合法函数拟合法t )V1V1(SKBAABABABBAABBAeVV)(VVVVV) t (C前面得到的模型是一个带初值的一阶线性微分方程，解前面得到的模型是一个带初值的一阶线性微分方程，解之得：之得：问题归结为利用问题归结为利用C CB B在时刻在时刻t tj j的测量数据的测量数据C Cj j(j=1,2,.,N)(j=1,2,.,N)来辨识来辨识 K K 和和 。BA,BA 、引入引入BAABABABBAAVV)(VbVVVVa，t )V1V1(SKBBAbea) t (C从而从而Kt02. 0Bbea) t (C 用函数用函数CB(t)来拟合所

15、给的实验数据，从而来拟合所给的实验数据，从而估计出其中的参数估计出其中的参数a,b,K。23BAcm10S,cm1000VV将将代入上式有：代入上式有：用用MATLABMATLAB软件进行计算软件进行计算. .1 1编写函数编写函数M-M-文件文件 nongdu.m nongdu.mfunction f = nongdu(x,tdata)function f = nongdu(x,tdata)f = x(1)+x(2)f = x(1)+x(2)* *exp(-0.02exp(-0.02* *x(3)x(3)* *tdata);tdata);其中其中 x(1) = a;x(2) = b;x(3)

16、 = k; x(1) = a;x(2) = b;x(3) = k;2) 2) 在工作空间中执行以下命令在工作空间中执行以下命令(test1.m)(test1.m) tdata = linspace(100,1000,10); tdata = linspace(100,1000,10); cdata =4.54 4.99 5.35 5.65 5.90 6.10 . cdata =4.54 4.99 5.35 5.65 5.90 6.10 . 6.26 6.39 6.50 6.59; 6.26 6.39 6.50 6.59; x0 = 0.2,0.05,0.05; x0 = 0.2,0.05,0.

17、05; x=lsqcurvefit ( x=lsqcurvefit (nongdunongdu,x0,tdata,cdata),x0,tdata,cdata)3) 3) 输出结果输出结果: x = 0.007 -0.003 0.1012: x = 0.007 -0.003 0.1012 即即 k = 0.1012, a = 0.007, b = -0.003, k = 0.1012, a = 0.007, b = -0.003, nongdutest1进一步求得：进一步求得：)cm/mg(01. 0),cm/mg(004. 03A3B2. 非线性规划法非线性规划法 利用利用C CB B在时刻在

18、时刻t tj j的测量数据的测量数据C Cj j(j=1,2,.,N)(j=1,2,.,N)来辨识来辨识 K K 和和 。BA,N1j2jjBBAmin)C)t (C(),K(E问题可转化为求函数问题可转化为求函数202. 0),(jKtCbeabaKEj即求函数即求函数的最小值点的最小值点K K，a a，b b。3. 导函数拟合法导函数拟合法前面得到的微分方程为：前面得到的微分方程为：C02. 0)(01. 0kdtdCBBAB令令BAh上式变为：上式变为：)C02. 0h01. 0(kdtdCBB这可以看作这可以看作dtdCB随随CB的变化规律的变化规律(j=1, 2,.,N)若知道一组数

19、据若知道一组数据)dtdC( ,CjBj则可用最小二乘拟合的方法来求出函数则可用最小二乘拟合的方法来求出函数dtdCB中的未知参数中的未知参数K和和h。即为求参数即为求参数K, a使以下误差函数到达最小：使以下误差函数到达最小：2jjBN1j)C02. 0a01. 0(K)dtdC()a,K( J 该问题等价于用函该问题等价于用函 数数 f(K,a,CBB)来拟合数据来拟合数据)dtdC( ,CjBj(j=1, 2,.,N)用用MATLABMATLAB软件进行计算软件进行计算. .)dtdC( ,CjBj求数据点求数据点(j=1, 2,.,N)tdata = linspace(100,1000

20、,10);cdata = 1e-05.*454 499 535 565 590. 610 626 639 650 659;d,ifail=e01bef(tdata,cdata);cj,dcj=e01bgf(tdata,cdata,d,tdata);1 1编写函数编写函数M-M-文件文件 baomof.m baomof.mfunction f=baomof(x,cdata)function f=baomof(x,cdata)f=x(1)f=x(1)* *(0.01(0.01* *x(2)-0.02x(2)-0.02* *cdata)cdata)其中其中 x(1) = K; x(2) =h x(1

21、) = K; x(2) =h2) 2) 编写命令编写命令M M文件文件(baomo21.m)(baomo21.m)3) 3) 输出结果输出结果: x = 0.1009 0.014: x = 0.1009 0.014 即即 k = 0.1009, h = 0.014k = 0.1009, h = 0.014作函数拟合作函数拟合x0=0.2,0.1;x=lsqcurvefit (baomof,x0,cdata,dcj)4.线性化迭代法线性化迭代法Kt02. 0Bbea) t (C前面带初始条件的一阶线性微分方程的解为前面带初始条件的一阶线性微分方程的解为其中：其中：)(5 . 0b)(5 . 0a

22、ABBA， 如果得到了参数K的一个较好的近似值K*，那么将CB(t)关于K在K*处展开，略去K的二次及以上的项得CB(t)的一个近似式teKbbeatCtKtKB*02. 002. 002. 0)(通过极小化KbdtedbeaCKbaEjtKtKnjjjj02. 0),(202. 002. 01*其中确定确定a, b, d, 再由再由 K=d/0.02b得到得到K*的修正值的修正值 K。K*K*- K, 得到得到K的一个新的近似值，用同样的方的一个新的近似值，用同样的方法再求新的修正值法再求新的修正值 K。这个过程可以不断重复，。这个过程可以不断重复，直到修正值足够小为止。直到修正值足够小为止

23、。1当当K的初值取为的初值取为k=0.3时，出现奇异情况，迭代不时，出现奇异情况，迭代不收敛；收敛；2当当K的初值取为的初值取为k=0.2时，经四次迭代，已经收敛到时，经四次迭代，已经收敛到一个很好的解。迭代结果如下表。一个很好的解。迭代结果如下表。迭代次数abk误差 E一6.7312-2.54410.06292.3676二6.5933-2.57820.11230.2963三6.9750-2.98180.10124.5779e-004四6.9850-2.99410.10125.6531e-005五、结果及误差分析五、结果及误差分析 几种方法得出的结果及相应的误差总结于下表，几种方法得出的结果及

24、相应的误差总结于下表，误差为计算数据与实验数据之差的平方和。误差为计算数据与实验数据之差的平方和。注：导函数拟合法得出的参数值精度有限，线性化注：导函数拟合法得出的参数值精度有限，线性化迭代法要求参数的初值比较接近精确值。因此可将迭代法要求参数的初值比较接近精确值。因此可将导函数拟合法和线性化迭代法结合起来使用，把前导函数拟合法和线性化迭代法结合起来使用，把前者得到的参数者得到的参数K的值作为迭代法中的值作为迭代法中K的初值，这样可的初值，这样可使迭代法收敛或收敛更快。使迭代法收敛或收敛更快。3取取K的初值为的初值为k=0.1009,只一次迭代就得到只一次迭代就得到2中的最中的最后结果。后结果。方法abk误差 E函数拟合6.9852-2.99410.10125.6541e-005非线性规划6.9143-2.99010.10962.0248e-003导函数拟合6.98480.1009线性化迭代法6.9850-2.99410.10025.6531e-005函数拟合法的拟合效果函数拟合法的拟合效果求解参数辨识模型的方法：求解参数辨识模型的方法： 函数拟合；函数拟合； 非线性规划；非线性规划； 导函数拟合；导函数拟合； 线性化迭代；线性化迭代； 其它方法其它方法 。 布置布置“函数拟合实验函数拟合实验目的目的 1 1. 掌握用掌握用MATLABMATLAB计算函数拟合