1stopt 7参数拟合求助
需要根据两个式子（一个微分方程，一个代数方程）拟合7个未知参数fac(1)~fac(7)，网上下载的软件有参数限制，想向有正版软件的朋友求助，帮忙拟合一下。写了两版，程序分别如下。完整程序和数据放在了文件里，不知为什么上传失败，只能贴一下网盘链接，麻烦大家了。
链接：https://pan.baidu.com/s/18yjg6HL9hhQ5hmjGdjFhEQ
提取码：8xmm
第一种：
//title "fit";
parameter fac(1:7);
variable t,lamb,pa,p,lambv;
odefunction lambv'=((1/3/fac(7))*(fac(1)*((lamb/lambv)^fac(2)-(lambv/lamb)^(0.5*fac(2)))+fac(3)*((lamb/lambv)^fac(4)-(lambv/lamb)^(0.5*fac(4)))+fac(5)*((lamb/lambv)^fac(6)-(lambv/lamb)^(0.5*fac(6)))));
p=(pa+fac(1)*((lamb/lambv)^fac(2)/lamb-(lambv/lamb)^(fac(2)*0.5)/lamb)+fac(3)*((lamb/lambv)^fac(4)/lamb-(lambv/lamb)^(fac(4)*0.5)/lamb)+fac(5)*((lamb/lambv).^fac(6)/lamb-(lambv/lamb)^(fac(6)*0.5)/lamb));
//根据上面两个式子拟合fac(1)~fac(7)，共计7个参数。lambv为中间变量
data;
//t,lamb,pa,p
第二种：
//title "fit";
parameter fac(1:7);
variable t,lamb,pa,p;
conststr lambvdif=((1/3/fac(7))*(fac(1)*((lamb/lambv)^fac(2)-(lambv/lamb)^(0.5*fac(2)))+fac(3)*((lamb/lambv)^fac(4)-(lambv/lamb)^(0.5*fac(4)))+fac(5)*((lamb/lambv)^fac(6)-(lambv/lamb)^(0.5*fac(6)))));
function p=(pa+fac(1)*((lamb/lambv)^fac(2)/lamb-(lambv/lamb)^(fac(2)*0.5)/lamb)+fac(3)*((lamb/lambv)^fac(4)/lamb-(lambv/lamb)^(fac(4)*0.5)/lamb)+fac(5)*((lamb/lambv).^fac(6)/lamb-(lambv/lamb)^(fac(6)*0.5)/lamb));
//根据上面两个式子拟合fac(1)~fac(7)，共计7个参数。lambv为中间变量
data;
//t,lamb,pa,p
返回小木虫查看更多
京公网安备 11010802022153号
图形参考：https://blog.csdn.net/wlfc/article/details/119539076
使用 3楼、4楼 的公式进行求解。
有两个难点：
1、微分方程求解时要传递中间变量lamb, pa。
2、lambv为未知中间变量，需要求解。
第一个难点对Lu脚本来说不是问题，gsl支持的Lu扩展数学库中求解微分方程的函数gsl_ode中提供了该功能。
进一步讨论：微分方程求解传递中间变量时，按时间t的改变量进行插值，可提高求解精度。例如，进行线性插值：lamb=lamb0+(lamb1-lamb0)*[(t-t0)/(t1-t0)], pa=pa0+(pa1-pa0)*[(t-t0)/(t1-t0)]。
第二个难点若使用6楼的办法，当数据量大时似乎难以实现。故在这里使用幂级数逼近lambv，进行拟合的方法，需增加k0,k1,k2,k3四个拟合参数；当然拟合参数多时精度高，但耗时长。
仍使用部分数据拟合，以减少时间。
Lu脚本代码：
!!!using["luopt","math"]; //使用命名空间
f(t,p,dp, i : lamb,pa, A,w,lambdif,lambvdif,lambv: t_A, t_lamb, t_pa, fac1,fac2,fac3,fac4,fac5,fac6,fac7,k0,k1,k2,k3,k4) =
{
lamb=t_lamb+(t_lamb[i+1]-t_lamb)*[(t-t_A)/(t_A[i+1]-t_A)], pa=t_pa+(t_pa[i+1]-t_pa)*[(t-t_A)/(t_A[i+1]-t_A)], //线性插值计算lamb,pa
lambv=k0+k1*lamb+k2*lamb*lamb+k3*lamb^3+k4*lamb^4, //lambv=f(lamb,k0,k1,k2,k3,k4)，使用幂级数逼近
A=0.02,w=2*pi*0.1,
lambdif=A*w*cos(w*t),
lambvdif=((1.0/3/fac7)*(fac1*((lamb/lambv)^fac2-(lambv/lamb)^(0.5*fac2))+fac3*((lamb/lambv)^fac4-(lambv/lamb)^(0.5*fac4))+fac5*((lamb/lambv)^fac6-(lambv/lamb)^(0.5*fac6)))),
dp=pa+fac1*((fac2-1)*lamb^(fac2-2)*lambv^(-fac2)*lambdif-fac2*lamb^(fac2-1)*lambv^(-fac2-1)*lambvdif
+(0.5*fac2+1)*lamb^(-0.5*fac2-2)*lambv^(0.5*fac2)*lambdif-0.5*fac2*lamb^(-0.5*fac2-1)*lambv^(0.5*fac2-1)*lambvdif)
+fac3*((fac4-1)*lamb^(fac4-2)*lambv^(-fac4)*lambdif-fac4*lamb^(fac4-1)*lambv^(-fac4-1)*lambvdif+(0.5*fac4+1)*lamb^(-0.5*fac4-2)*lambv^(0.5*fac4)*lambdif
-0.5*fac4*lamb^(-0.5*fac4-1)*lambv^(0.5*fac4-1)*lambvdif)+fac5*((fac6-1)*lamb^(fac6-2)*lambv^(-fac6)*lambdif
-fac6*lamb^(fac6-1)*lambv^(-fac6-1)*lambvdif+(0.5*fac6+1)*lamb^(-0.5*fac6-2)*lambv^(0.5*fac6)*lambdif
-0.5*fac6*lamb^(-0.5*fac6-1)*lambv^(0.5*fac6-1)*lambvdif),
0 //必须返回0
};
目标函数(_fac1,_fac2,_fac3,_fac4,_fac5,_fac6,_fac7,_k0,_k1,_k2,_k3,_k4 : i,s,tf: t_A, t_p, fac1,fac2,fac3,fac4,fac5,fac6,fac7,k0,k1,k2,k3,k4)=
{
fac1=_fac1, fac2=_fac2, fac3=_fac3, fac4=_fac4, fac5=_fac5, fac6=_fac6, fac7=_fac7, k0=_k0, k1=_k1, k2=_k2, k3=_k3, k4=_k4, //传递优化变量
//最后一个参数50表示gsl_ode函数在计算时,最多循环计算50次,这样可以提高速度
tf=gsl_ode[@f, nil, nil, t_A, ra1(-0.060966222), 1e-6, 1e-6, gsl_rkf45, 1e-6,50],
sum{[tf(all:1).reshape()-t_p].^2.0} //代码矢量化加速，特别适合大数据量
};
main(: tArray : t_A, t_lamb, t_pa, t_p)=
{
tArray=matrix{ //存放实验数据//t,lamb,pa,p
"0 0.959926667 -0.218715994 -0.060966222
0.25 0.963063333 -0.200840559 -0.046836278
0.5 0.96612 -0.183553472 -0.035042256
0.75 0.969016667 -0.167289626 -0.024438988
1 0.97171 -0.152268899 -0.01706765
1.25 0.97409 -0.139075651 -0.010939712
1.5 0.976123333 -0.127862697 -0.006795087
1.75 0.977786667 -0.118729759 -0.002919842
2 0.979000001 -0.112089924 -0.001981819
2.25 0.979743333 -0.108031322 -0.001163676
2.5 0.98001 -0.106577018 -0.003491456
2.75 0.979766667 -0.107904033 -0.009723697"
},
t_A=tArray(all:0).reshape(), t_lamb=tArray(all:1).reshape(), t_pa=tArray(all:2).reshape(), t_p=tArray(all:3).reshape(), //t_A等取矩阵的列，并转为一维数组
Opt1[@目标函数, optwaysimdeep, optwayconfra] //Opt1函数全局优化
};
结果（fac1～fac7，k0,k1,k2,k3,k4，最小值）：
17.38947868627771 0.3416680168065553 -31901.02726305072 1.871542771429994e-004 2.297270203567084e-003 8.207606707403809 -41.52220237484659 -1.705988163594953 -0.3001058349255158 2.1410790976942 4.419612473682794 -4.029380465220802 2.405714024195563e-006
绘图代码：
!!!using["luopt","math","win"]; //使用命名空间
f(t,p,dp, i : lamb,pa, A,w,lambdif,lambvdif,lambv: t_A, t_lamb, t_pa, fac1,fac2,fac3,fac4,fac5,fac6,fac7,k0,k1,k2,k3,k4) =
{
lamb=t_lamb+(t_lamb[i+1]-t_lamb)*[(t-t_A)/(t_A[i+1]-t_A)], pa=t_pa+(t_pa[i+1]-t_pa)*[(t-t_A)/(t_A[i+1]-t_A)], //线性插值计算lamb,pa
lambv=k0+k1*lamb+k2*lamb*lamb+k3*lamb^3+k4*lamb^4, //lambv=f(lamb,k0,k1,k2,k3,k4)，使用幂级数逼近
A=0.02,w=2*pi*0.1,
lambdif=A*w*cos(w*t),
lambvdif=((1.0/3/fac7)*(fac1*((lamb/lambv)^fac2-(lambv/lamb)^(0.5*fac2))+fac3*((lamb/lambv)^fac4-(lambv/lamb)^(0.5*fac4))+fac5*((lamb/lambv)^fac6-(lambv/lamb)^(0.5*fac6)))),
dp=pa+fac1*((fac2-1)*lamb^(fac2-2)*lambv^(-fac2)*lambdif-fac2*lamb^(fac2-1)*lambv^(-fac2-1)*lambvdif
+(0.5*fac2+1)*lamb^(-0.5*fac2-2)*lambv^(0.5*fac2)*lambdif-0.5*fac2*lamb^(-0.5*fac2-1)*lambv^(0.5*fac2-1)*lambvdif)
+fac3*((fac4-1)*lamb^(fac4-2)*lambv^(-fac4)*lambdif-fac4*lamb^(fac4-1)*lambv^(-fac4-1)*lambvdif+(0.5*fac4+1)*lamb^(-0.5*fac4-2)*lambv^(0.5*fac4)*lambdif
-0.5*fac4*lamb^(-0.5*fac4-1)*lambv^(0.5*fac4-1)*lambvdif)+fac5*((fac6-1)*lamb^(fac6-2)*lambv^(-fac6)*lambdif
-fac6*lamb^(fac6-1)*lambv^(-fac6-1)*lambvdif+(0.5*fac6+1)*lamb^(-0.5*fac6-2)*lambv^(0.5*fac6)*lambdif
-0.5*fac6*lamb^(-0.5*fac6-1)*lambv^(0.5*fac6-1)*lambvdif),
0 //必须返回0
};
set(_fac1,_fac2,_fac3,_fac4,_fac5,_fac6,_fac7,_k0,_k1,_k2,_k3,_k4 :: fac1,fac2,fac3,fac4,fac5,fac6,fac7,k0,k1,k2,k3,k4)=
{
fac1=_fac1, fac2=_fac2, fac3=_fac3, fac4=_fac4, fac5=_fac5, fac6=_fac6, fac7=_fac7, k0=_k0, k1=_k1, k2=_k2, k3=_k3, k4=_k4 //传递优化变量
};
init(init: tArray,t_lambv,i,tf,k,lamb : t_A, t_lamb, t_pa, t_p, max, fac1,fac2,fac3,fac4,fac5,fac6,fac7,k0,k1,k2,k3,k4)=
{
tArray=matrix{ //存放实验数据//t,lamb,pa,p
"0 0.959926667 -0.218715994 -0.060966222
... ... 数据省略
2.75 0.979766667 -0.107904033 -0.009723697"
},
len[tArray,0,&max], t_A=tArray(all:0).reshape(), t_lamb=tArray(all:1).reshape(), t_pa=tArray(all:2).reshape(), t_p=tArray(all:3).reshape(),
set[17.38947868627771, 0.3416680168065553, -31901.02726305072, 1.871542771429994e-004, 2.297270203567084e-003, 8.207606707403809, -41.52220237484659, -1.705988163594953, -0.3001058349255158, 2.1410790976942, 4.419612473682794, -4.029380465220802 ],
t_lambv=new[real_s,max],
i=-1, while{++i<max, lamb=t_lamb, t_lambv=k0+k1*lamb+k2*lamb*lamb+k3*lamb^3+k4*lamb^4},
tf=gsl_ode[@f, nil, nil, t_A, ra1(-0.060966222), 1e-6, 1e-6, gsl_rkf45, 1e-6,50],
cwAttach[typeSplit], cwResizePlots(1,2,2), //左二右二分裂
k=cwAddCurve{t_A, t_p, max, 0}, //给0子图添加曲线
cwSetScatter(k,0), //设置绘制点
cwSetDataLineSize(5, k, 0), //设置点的大小
cwAddCurve{t_A, tf(all:1).reshape(), max, 0}, //给0子图添加曲线
cwAddCurve{t_A, t_lamb, max, 1}, //给1子图添加曲线
cwAddCurve{t_A, t_pa, max, 2}, //给2子图添加曲线
cwAddCurve{t_A, t_lambv, max, 3} //给3子图添加曲线
};
ChartWnd[@init];
图形：
https://blog.csdn.net/wlfc/article/details/119607664
非常感谢，您博客里给出的图形拟合效果很好，我再仔细看一下您的答复
7楼、8楼、9楼的回复，因理解有误，故重新解答。
仍使用 1楼 给出的第一种公式。仍使用部分数据拟合。
分析：需在微分方程求解中传递中间变量lamb，gsl支持的Lu扩展数学库中求解微分方程的函数gsl_ode提供了传递中间变量的功能。
进一步讨论：微分方程求解传递中间变量时，按时间t的改变量进行插值，可提高求解精度。例如，进行线性插值：lamb=lamb0+(lamb1-lamb0)*[(t-t0)/(t1-t0)]。
另外，微分方程初值lambv未知，将其追加为拟合变量lambv0。
Lu代码：
结果：
7.148976005697928e-003 -14.5708961980305 -1.059296609366114 0.2342635320112431 -9.68049051324744e-005 -31.79361445870142 -2.026215042289879 1.184269268536115 3.55377979903292e-006
绘图代码：
图形参考：https://blog.csdn.net/wlfc/article/details/119830682
，