湖南人文科技学院袁野

1、湖南人文科技学院机电系机械优化设计班级：机械三班姓名：袁野学号： 11428326教师：禹翼题目 给药方案设计需要依据药物 吸收与排除过程的原理。药物进入 机体后随血液输送到全身，不断地 被吸收、分布、代谢，最终排出体 外。药物在血液中的浓度，即单位 体积血液中的药物含量，称血药浓 度。在最简单的一室模型中，将整个机体看作一个房室，称中心室，室内的血药浓度是均匀的。 这里我们用一室模型，讨论在口服给药方式下血药浓度的变化规律，及根据实验数据拟合参数的方法。口服给药方式相当于先有一个将药从肠胃吸收入血液的过程， 这个过程可简化为在药物 进入中心室之前有一个吸收室（如图） ，记中心室和吸收室的容积

2、分别为 V ，V1 ，而 t 时刻 的血药浓度分别为 c（ t）， c1（ t）；中心室的排除速率为 k，吸收速率为 k1（这里 k和k1 分。设 t=0 时刻口服剂量为 d 的别是中心室和吸收室血药浓度变化率与浓度本身的比例系数）药物，容易写出吸收室的血药浓度 c1（ t）的微分方程为dc1dtk1c1,c1(0)dV1中心室血药浓度 c（ t）的变化率由两部分组成：与c 成正比的排除（比例系数） ；与成正比的吸收（比例系数）再考虑到中心室和吸收室的容积分别为，得到（）的微分方程为dcdtkc V1 k1c1,c(0) 0.由以上两个微分方程不难解出中心室血药浓度c(t)d k1ktV k1

3、 k (ek1t).在制定给药方案时必须知道这种药物的3 个参数 k，k1，b（=d/v） ，实际中通常通过实验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数据确定。设 t=0 时刻口服一定剂量的药物，下表是实验数据c（t），请由此确定 k，k1， b。【模型建立】问题可以转换为利用数据拟合方程k1 ktk1tc(t) b (e e 1 ) k1 kb=d/V 这个方程是非线性的。由于这是一个非线性最小二乘拟合的问题，故可以选用 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);【输出结果】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);输出结果】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 命令所计算得到的结果基 本上是一致的。【结果分析】( 1) Isqnonlin 和 Isqcurvefit 的对比可以分别使用 Isqnonlin 和 Isqcurvefit 求解非线性最小二乘拟合。实验结果发现结果基 本一致。这是由两者所用的算法完全一样所决定的。此外，还可以对比分析采用分析导数(雅克比矩阵) 、对比分析步长一维搜索算法(混 合二三次插值和三次插值) 等对 Matlab 的非线性二乘拟合算法进行探索， 所得结果都是一致 的。(2)LM 法和 GN法的对比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法运行结果如下：>> 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 函数调用时默认使用的下降方向算法为LM法。可以通过更改控制参数中的LevenbergMarquard 而使用 GN算法求解非线性最小二乘法。结果发现，两种下降方向算法 所得结果是完全一致的。【实验结论】利用lsqnonlin命令和 lsqcurvefit命令作数据拟合， 得到的结果基本一致。 对于本题， 可以 得到b =46.8275 ,k1=3.6212 ,k=0.2803 。用Isqnonlin命令所得到的结果要