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計算函數擬合