湖南人文科技學(xué)院袁野

1、湖南人文科技學(xué)院機(jī)電系機(jī)械優(yōu)化設(shè)計(jì)班級(jí)：機(jī)械三班姓名：袁野學(xué)號(hào)： 11428326教師：禹翼題目 給藥方案設(shè)計(jì)需要依據(jù)藥物 吸收與排除過程的原理。藥物進(jìn)入 機(jī)體后隨血液輸送到全身，不斷地 被吸收、分布、代謝，最終排出體 外。藥物在血液中的濃度，即單位 體積血液中的藥物含量，稱血藥濃 度。在最簡(jiǎn)單的一室模型中，將整個(gè)機(jī)體看作一個(gè)房室，稱中心室，室內(nèi)的血藥濃度是均勻的。 這里我們用一室模型，討論在口服給藥方式下血藥濃度的變化規(guī)律，及根據(jù)實(shí)驗(yàn)數(shù)據(jù)擬合參數(shù)的方法?？诜o藥方式相當(dāng)于先有一個(gè)將藥從腸胃吸收入血液的過程， 這個(gè)過程可簡(jiǎn)化為在藥物 進(jìn)入中心室之前有一個(gè)吸收室（如圖） ，記中心室和吸收室的容積

2、分別為 V ，V1 ，而 t 時(shí)刻 的血藥濃度分別為 c（ t）， c1（ t）；中心室的排除速率為 k，吸收速率為 k1（這里 k和k1 分。設(shè) t=0 時(shí)刻口服劑量為 d 的別是中心室和吸收室血藥濃度變化率與濃度本身的比例系數(shù)）藥物，容易寫出吸收室的血藥濃度 c1（ t）的微分方程為dc1dtk1c1,c1(0)dV1中心室血藥濃度 c（ t）的變化率由兩部分組成：與c 成正比的排除（比例系數(shù)） ；與成正比的吸收（比例系數(shù)）再考慮到中心室和吸收室的容積分別為，得到（）的微分方程為dcdtkc V1 k1c1,c(0) 0.由以上兩個(gè)微分方程不難解出中心室血藥濃度c(t)d k1ktV k1

3、 k (ek1t).在制定給藥方案時(shí)必須知道這種藥物的3 個(gè)參數(shù) k，k1，b（=d/v） ，實(shí)際中通常通過實(shí)驗(yàn)t0.0830.1670.250.50.7511.5c(t)10.921.127.336.435.538.434.8t2.2534681012c(t)24.223.615.78.28.32.21.8數(shù)據(jù)確定。設(shè) t=0 時(shí)刻口服一定劑量的藥物，下表是實(shí)驗(yàn)數(shù)據(jù)c（t），請(qǐng)由此確定 k，k1， b?！灸Ｐ徒ⅰ繂栴}可以轉(zhuǎn)換為利用數(shù)據(jù)擬合方程k1 ktk1tc(t) b (e e 1 ) k1 kb=d/V 這個(gè)方程是非線性的。由于這是一個(gè)非線性最小二乘擬合的問題，故可以選用 MA TLA

4、B 中 的 Isqnonlin 命令或 Isqcurvefit 命令。【程序】function f=medicine(x,t)%x(1)-b,x(2)-k1,x(3)-kf = x(1)*x(2)/(x(2) - x(3)*(exp(-x(3)*t)-exp(-x(2)*t);方法一 用 Isqcurvefit 命令x0=1,1,0;t = 0.083 0.167 0.25 0.50 0.75 1.0 1.5 2.25 3.0 4.0 6.0 8.0 10.0 12.0;c = 10.9 21.1 27.3 36.4 35.5 38.4 34.8 24.2 23.6 15.7 8.2 8.3

5、2.21.8 ;opt = optimset;x,norm,res,ef,out,lam = lsqcurvefit(medicine,x0,t,c,opt);【輸出結(jié)果】x =46.8275 3.6212 0.2803norm =34.2317res = 1.1090 -0.3885 -0.5052 -0.5820 2.2767 -1.4071 -1.68592.8020 -1.7058 0.8438 1.2452 -2.9075 0.8787 -0.0423ef =4out =iterations: 34funcCount: 255stepsize: 1cgiterations: firs

6、torderopt: algorithm: 'medium-scale: Levenberg-Marquardt, line-search'message: 1x57 charlam = upper: lower: 方法二 用 Isqnonlin 命令 x0 = 50,1,0;t = 0.083 0.167 0.25 0.50 0.75 1.0 1.5 2.25 3.0 4.0 6.0 8.0 10.0 12.0; c = 10.9 21.1 27.3 36.4 35.5 38.4 34.8 24.2 23.6 15.7 8.2 8.3 2.21.8 ;opt=optimset

7、( 'LargeScale' , 'on' , 'MaxFunEvals' ,1000, 'MaxIter' ,200); opt=optimset(opt, 'tolx' ,1e-16, 'tolf' ,1e-16);x,norm,res,ef,out,lam=lsqnonlin(medicine,x0,t,c);輸出結(jié)果】x =46.8217 3.6219 0.2802norm =34.2318res =1.1096 -0.3881 -0.5051 -0.5835 2.2743 -1.4099

8、-1.6881 2.8011 -1.70560.8449 1.2469 -2.9060 0.8799 -0.0414ef =3out = firstorderopt: 0.0058iterations: 18funcCount: 76 cgiterations: 18algorithm: 'large-scale: trust-region reflective Newton'message: 1x87 char lam =lower: 3x1 doubleupper: 3x1 doubled=46.8217,k=3.6219,k1=0.2802 ，norm =34.2318

9、。這與用 Isqcurvefit 命令所計(jì)算得到的結(jié)果基 本上是一致的。【結(jié)果分析】( 1) Isqnonlin 和 Isqcurvefit 的對(duì)比可以分別使用 Isqnonlin 和 Isqcurvefit 求解非線性最小二乘擬合。實(shí)驗(yàn)結(jié)果發(fā)現(xiàn)結(jié)果基 本一致。這是由兩者所用的算法完全一樣所決定的。此外，還可以對(duì)比分析采用分析導(dǎo)數(shù)(雅克比矩陣) 、對(duì)比分析步長(zhǎng)一維搜索算法(混 合二三次插值和三次插值) 等對(duì) Matlab 的非線性二乘擬合算法進(jìn)行探索， 所得結(jié)果都是一致 的。(2)LM 法和 GN法的對(duì)比x0=1,1,0;t = 0.083 0.167 0.25 0.50 0.75 1.0 1

10、.5 2.25 3.0 4.0 6.0 8.0 10.0 12.0; c = 10.9 21.1 27.3 36.4 35.5 38.4 34.8 24.2 23.6 15.7 8.2 8.3 2.21.8 ;opt1=optimset( 'MaxFunEvals' ,2000);x1,norm1,res1=lsqcurvefit(medicine,x0,t,c,opt1);% LM 法opt2=optimset(opt1, 'LevenbergMarquardt' , 'off' );x2,norm2,res2=lsqcurvefit(medi

11、cine,x0,t,c,opt2);% GN法運(yùn)行結(jié)果如下：>> x1 x1 =46.8275 3.6212 0.2803>> out1out1 =firstorderopt: 5.4286e-004 iterations: 20 funcCount: 84cgiterations: 0algorithm: 'large-scale: trust-region reflective Newton'message: 1x463 char>> x2x2 =46.8275 3.6212 0.2803>> out2out2 =firsto

12、rderopt: 5.4286e-004iterations: 20funcCount: 84cgiterations: 0algorithm: 'large-scale: trust-region reflective Newton'message: 1x463 charlsqcurvefit 函數(shù)調(diào)用時(shí)默認(rèn)使用的下降方向算法為L(zhǎng)M法。可以通過更改控制參數(shù)中的LevenbergMarquard 而使用 GN算法求解非線性最小二乘法。結(jié)果發(fā)現(xiàn)，兩種下降方向算法 所得結(jié)果是完全一致的。【實(shí)驗(yàn)結(jié)論】利用lsqnonlin命令和 lsqcurvefit命令作數(shù)據(jù)擬合， 得到的結(jié)果基本一致。 對(duì)于本題， 可以 得到b =46.8275 ,k1=3.6212 ,k=0.2803 。用Isqnonlin命令所得到的結(jié)果要