应《网络安全法》要求，自2017年10月1日起，未进行实名认证将不得使用互联网跟帖服务。为保障您的帐号能够正常使用，请尽快对帐号进行手机号验证，感谢您的理解与支持！

24小时热门版块排行榜    

小木虫论坛-学术科研互动平台 » 计算模拟区 » 计算模拟 » 求助常微分方程组的参数优化问题，有没有师兄可以用1stopt帮忙试试？
51/1返回列表
查看: 1411  |  回复: 5
只看楼主 @他人 存档 新回复提醒 (忽略) 收藏    在APP中查看
当前只显示满足指定条件的回帖，点击这里查看本话题的所有回帖

快乐的小尼玛

新虫 (初入文坛)

[求助] 求助常微分方程组的参数优化问题，有没有师兄可以用1stopt帮忙试试？ 已有1人参与

师兄师姐们好！

  我最近遇到一个常微分方程组的参数优化问题。方程如下：
dT <- -beta*T*V-h*F*T+rho*R
dI <- beta*T*V -delta*I-k*I*F
dR <- h*F*T-rho*R
dV <- p*I-c*V
dF <- q*I-d*F
其中：
delta = 2*exp(de*(t-7))
beta,p,c,k,h,rho,q,d,de是待优化参数
初值：
T = 2e7
I = 0
R = 0
V = 3
F = 1
需要拟合F和V的实验值：
t        V        F
1        3.5        1.63
2        4.5        2.67
3        5.5        5.35
4        6.5        1.09
5        4.5        1.32
6        1.5        0
7        1.5        0
8        0        0
我尝试用r的minpack.lm包求解，但是试了两万多组初始值，始终达不到好的效果。
目前拟合最好的一组参数如下：
beta =1.7e-07
p = 40
c = 17
k = 0.85
h = 0.2
rho = 1e-07
q = 3e-06
d = 0.74501237318802
de = 1
图在附件：

本人数学功底一般，刚接触这方面话题，有什么说的不清楚的地方请大家指出，我再补充。
平时潜水较多，只有10个金币啦，希望各位师兄师姐帮忙！谢谢！
下面是我的r代码，我有说不清的大家可以看代码（解ode并画图的，不是优化的）。
#Function
my_ode <- function(V_input,F_input,beta = beta,
        de = de,
        p = p,
        c = c,
        k = k,
        h = h,
        rho = rho,
        q = q,
        d = d,
        miu = miu,
        state)
{

        library(deSolve)
        library(minpack.lm)
        library(rlist)

        parameters <- c(beta = beta,
        de = de,
        p = p,
        c = c,
        k = k,
        h = h,
        rho = rho,
        q = q,
        d = d,
        miu=miu)
        Target<-function(t, state, parameters) {
        with(as.list(c(state, parameters)),{
                        delta = 2*exp(de*(t-miu))
                        dT <- -beta*T*V-h*F*T+rho*R
                        dI <- beta*T*V -delta*I-k*I*F
                        dR <- h*F*T-rho*R
                        dV <- p*I-c*V
                        dF <- q*I-d*F
                        list(c(dT, dI, dR, dV, dF))
                })
        }
        times <- seq(0,8,by = 0.01)
        out=ode(y=state,times=times,func=Target,parms=parameters)
        residue = c(log10(out[100,5])-3.5,log10(out[200,5])-4.5,log10(out[300,5])-5.5,log10(out[400,5])-6.5,log10(out[500,5])-4.5,log10(out[600,5])-1.5,log10(out[700,5])-1.5,out[100, 6]-1.64,out[200,6]-2.676,out[300,6]-5.35,out[400,6]-1.09,out[500,6]-1.33)
        residue = sqrt(sum(residue[1:7]^2)/7 + sum(residue[8:12]^2)/5)

        V_curve = log10(out[,5])
        F_curve = out[,6]



        SOUTH<-1; WEST<-2; NORTH<-3; EAST<-4;
         
        GenericFigure <- function(ID, size1, size2, point, line)
        {
          plot(0:8, 0:8, type="n", xlab="X", ylab="Y"
          lines(times,line, col = 2, cex = 2)
          points(point,col=3, pch = 15)
          #text(5,5, ID, col="red", cex=size1)
          box("plot", col="red"
          #mtext(paste("cex",size2,sep="", SOUTH, line=3, adj=1.0, cex=size2, col="blue"
          title(paste("Figure",ID,"(",residue,"",sep=" ")
        }
         
        MultipleFigures <- function()
        {
          GenericFigure("V", 3, 1, V_input, V_curve)
          box("figure", lty="dotted", col="blue"
         
          GenericFigure("F", 3, 1, F_input, F_curve)
          box("figure", lty="dotted", col="blue"
        }
        par(mfrow=c(1,2))
         
        MultipleFigures()
}
#Experimental data

#Initial values
state <- c(T = 2e7,
        I = 0,
        R = 0,
        V = 3,
        F = 1)
#Parameters

leiode<- function(beta,p,c,k,h,rho,q,d,de,miu)
{
        my_ode(V_input,F_input,
        beta = beta,
        de = de,
        p = p,
        c = c,
        k = k,
        h = h,
        rho = rho,
        q = q,
        d = d,
        miu = miu,
        state)
}

#1
V_input = c(3.5,4.5,5.5,6.5,4.5,1.5,1.5,0,0,0)
F_input = c(1.63,2.67,5.35,1.09,1.32)

leiode(1.7e-07, 40, 17, 0.85, 0.2, 1e-07, 3e-06, 0.74501237318802,1,7)

求助常微分方程组的参数优化问题，有没有师兄可以用1stopt帮忙试试？
rplot.png
回复此楼

» 猜你喜欢
1楼 2016-02-26 00:32:17
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

快乐的小尼玛

新虫 (初入文坛)
引用回帖:
2楼: Originally posted by dingd at 2016-02-26 16:12:27
1stOpt求这种微分方程拟合问题代码要简单的多，效果也好。

ConstStr delta = 2*exp(de*(tt-7));
InitialODEValue tt=0,T = 2e7,I = 0,R = 0,V = 3,F = 1;
Variable tt,V,F;
ODEFunction T' = -beta1*T*V-h*F* ...

谢谢师兄！
回复此楼
3楼2016-02-26 22:31:45
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
查看全部 6 个回答

dingd

铁杆木虫 (职业作家)

【答案】应助回帖

★ ★ ★ ★ ★ ★ ★ ★ ★ ★
感谢参与，应助指数 +1
快乐的小尼玛: 金币+10, ★★★★★最佳答案, 非常感谢！ 2016-02-26 22:30:05
1stOpt求这种微分方程拟合问题代码要简单的多，效果也好。
CODE:
ConstStr delta = 2*exp(de*(tt-7));
InitialODEValue tt=0,T = 2e7,I = 0,R = 0,V = 3,F = 1;
Variable tt,V,F;
ODEFunction T' = -beta1*T*V-h*F*T+rho*R;
            I' = beta1*T*V -delta*I-k*I*F;
            R' = h*F*T-rho*R;
            V' = p*I-c*V;
            F' = q*I-d*F;
Data;
tt        V        F
1        3.5        1.63
2        4.5        2.67
3        5.5        5.35
4        6.5        1.09
5        4.5        1.32
6        1.5        0
7        1.5        0
8        0        0

均方差(RMSE):0.85269322610406
残差平方和(SSE):11.6333718055
相关系数(R): 0.877968714324806
相关系数之平方(R^2): 0.770829063333154
修正R平方(Adj. R^2): 0.60992322741393
确定系数(DC): 0.771292669112812
F统计(F-Statistic): -1.4490058555454

参数                  最佳估算
--------------------        -------------
beta1        -6.71463644100751E-8
h        -7.08544427462871
rho        23.3711437688932
de        -0.897079268067582
k        -5.49659440033554
p        0.354150320382865
c        -0.265869901817656
q        0.225238201093483
d        -0.399871335562109
求助常微分方程组的参数优化问题，有没有师兄可以用1stopt帮忙试试？-1
c1.jpg


求助常微分方程组的参数优化问题，有没有师兄可以用1stopt帮忙试试？-2
c2.jpg
赞一下 回复此楼
2楼2016-02-26 16:12:27
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

快乐的小尼玛

新虫 (初入文坛)
引用回帖:
2楼: Originally posted by dingd at 2016-02-26 16:12:27
1stOpt求这种微分方程拟合问题代码要简单的多，效果也好。

ConstStr delta = 2*exp(de*(tt-7));
InitialODEValue tt=0,T = 2e7,I = 0,R = 0,V = 3,F = 1;
Variable tt,V,F;
ODEFunction T' = -beta1*T*V-h*F* ...

师兄抱歉，我发帖时候遗漏了对v的描述。
v的实际值我写错了，应该是：
tt        V        F
1        1e3.5        1.63
2        1e4.5        2.67
3        1e5.5        5.35
4        1e6.5        1.09
5        1e4.5        1.32
6        1e1.5        0
7        1e1.5        0
8        0        0
其中v的初始值也是不确定的，我随便选了个3.
我画线时候把v取了10的对数。抱歉师兄，麻烦您能帮我再试试吗？
赞一下 回复此楼
4楼2016-02-26 23:38:29
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

dingd

铁杆木虫 (职业作家)
1e3.5
这是什么数值表示法？是指10^3.5？
赞一下 回复此楼
5楼2016-03-02 17:15:43
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
查看全部 6 个回答
普通表情 龙 兔 虎 猫
信息提示
关闭
请填处理意见
关闭