当前位置: 首页 > 计算模拟 >最优化方法求教!!!

最优化方法求教!!!

作者 vvaa
来源: 小木虫 750 15 举报帖子
+关注

我的课题要做一个动力学模型,目前打算用集总的方法,建立一个微分方程组,比如
dc1/dt=-k1c1+k2c2
dc2/dt=k1c1-(k2+k3)c2+k4c3
dc3/dt=k3c2-(k4+k5)c3+k6c4
dc4/dt=k5c3-(k6+k7)c4
dc5/dt=k7c4
因为方程之间耦联,不能得到解析解。用bfgs法优化求解参数时,导函数不知道怎么求解。本人小白,没有最优化基础,希望大神能指导一下,怎么求dc/k,或者有没有比较好的不用导函数的优化方法推荐?

 返回小木虫查看更多

今日热帖
  • 精华评论
  • huab1984666

    function dy=myfunx(t,y,Parameters)

    %%
    k1=Parameters(1);
    k2=Parameters(2);
    k3=Parameters(3);
    k4=Parameters(4);
    k5=Parameters(5);
    k6=Parameters(6);
    k7=Parameters(7);
    %%
    dy=zeros(5,1);
    dy(1)=-k1*y(1)+k2*y(2);
    dy(2)=k1*y(1)-(k2+k3)*y(2)+k4*y(2);
    dy(3)=k3*y(2)-(k4+k5)*y(3)+k6*y(4);
    dy(4)=k5*y(5)-(k6+k7)*y(4);
    dy(5)=k7*y(4);
    %%%%


    function f=Kinetic(beta)
    global c;
    global t;
    Parameters=beta;
    t0=t(1);
    tn=t(end);
    %%
    options=odeset('RelTol',1e-4;'AbsTol',[1e-4 1e-4 1e-4 1e-4 1e-5]);
    [tx,cx]=ode15s(@myfunx,[t0 tn],[xx 0 0 0 0],options,Parameters);%%xx是t=0的时候,c的值。
    c_x=interp1(tx,cx,t,'spline');
    c1=c_x(:,1);c2=c_x(:,2);c3=c_x(:,3);c4=c_x(:,4);c5=c_x(:,5);
    ce=c1+c2+c3+c4+c5;
    N=length(t);
    for i=1:1:N
        rx(i)=(ce(i)-c(i)).^2;
    end
    f=sum(rx);


    global c; c=[- - - - -];%实验值
    global t; t=[- - - - -];实验值
    lb=[0 0 0 0 0 0 0];
    ub=[2 2 2 2 2 2 2];
    options =psoptimset('Display','Iter','MaxIter',2000,'TolFun',1e-8,'CompleteSearch','on');
    beta0=[0.7 0.01 0.01 0.0176 0.02 0.03 0.05];%% k1-7初始值
    [beta,fval,exitflag,output] = patternsearch(@Kinetic,beta0,[],[],[],[],lb,ub,options );
    %%
    Parameters=beta;
    t0=t(1);
    tn=t(end);
    %%
    Options=odeset('RelTol',1e-4,'AbsTol',[1e-4 1e-4 1e-4 1e-4 1e-5]);
    [tx,cx]=ode15s(@myfunx,[t0 tn],[1 0 0 0 0],Options,Parameters);
    c_x=interp1(tx,cx,t,'spline');
    c1=c_x(:,1);c2=c_x(:,2);c3=c_x(:,3);c4=c_x(:,4);c5=c_x(:,5);
    ce=c1+c2+c3+c4+c5;

    plot(t,c,'o');hold on;
    plot(t,ce,'r-');,

  • huab1984666

    引用回帖:
    2楼: Originally posted by huab1984666 at 2017-06-27 20:53:54
    function dy=myfunx(t,y,Parameters)

    %%
    k1=Parameters(1);
    k2=Parameters(2);
    k3=Parameters(3);
    k4=Parameters(4);
    k5=Parameters(5);
    k6=Parameters(6);
    k7=Parameters(7);
    %%
    dy=zeros(5,1);
    dy(1 ...

    global c; c=[- - - - -];%实验值
    global t; t=[- - - - -];实验值

    c,t都是列向量

  • huab1984666

    直接搜索,patternsearch

  • huab1984666

    引用回帖:
    7楼: Originally posted by vvaa at 2017-06-28 09:26:25
    嗯嗯,谢谢,我学习学习模式搜索法。还有一个问题,我每个实验条件下要做几个时间点才能求出来7个k值呢?
    ...

    这个看你自己啊。我觉得试验点越多越好

猜你喜欢
下载小木虫APP
与700万科研达人随时交流
  • 二维码
  • IOS
  • 安卓