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

24小时热门版块排行榜    

小木虫论坛-学术科研互动平台 » 计算模拟区 » 计算模拟 » 理论与计算方法 » 如何使用1stOpt进行三自变量三因变量的非线性拟合
71/1返回列表
查看: 1972  |  回复: 6
只看楼主 @他人 存档 新回复提醒 (忽略) 收藏    在APP中查看

mushixilin

铜虫 (初入文坛)

[求助] 如何使用1stOpt进行三自变量三因变量的非线性拟合 已有2人参与

正在学习使用1stOpt，见到的示例都是只有单个因变量的非线性拟合。想要使用1stopt对三元液液相平衡的NRTL热力学数据进行拟合，由于出现三个因变量，不太确定如何编程。编写的程序如下，如不使用sharedmodel，可以得到结果，但不理想；如使用sharedmodel，软件无法计算。还望前辈可以给予一定指导。

Title "NRTL simulation";
//sharedmodel;
Parameters alpha[0,1],tau12,tau13,tau21,tau23,tau31,tau32;
Variables x1,x2,x3,y1,y2,y3;

ConstStr F12=exp(-alpha*tau12);
ConstStr F13=exp(-alpha*tau13);
ConstStr F21=exp(-alpha*tau21);
ConstStr F23=exp(-alpha*tau23);
ConstStr F31=exp(-alpha*tau31);
ConstStr F32=exp(-alpha*tau32);

ConstStr lngammax1=(x1*0*1+x2*tau21*F21+x3*tau31*F31)/(x1*1+x2*F21+x3*F31)+
         (x1*1/(x1*1+x2*F21+x3*F31)*(0-(x1*0*1+x2*tau21*F21+x3*tau31*F31)/(x1*1+x2*F21+x3*F31)))+
         (x2*F12/(x1*F12+x2*1+x3*F32)*(tau12-(x1*tau12*F12+x2*0*1+x3*tau32*F32)/(x1*F12+x2*0+x3*F32)))+
         (x3*F13/(x1*F12+x2*F23+x3*1)*(tau13-(x1*tau13*F13+x2*tau23*F23+x3*0*1)/(x1*F13+x2*F23+x3*1)));
ConstStr lngammax2=(x1*tau12*F12+x2*0*1+x3*tau32*F32)/(x1*F12+x2*0+x3*F32)+
         (x1*F21/(x1*1+x2*F21+x3*F31)*(tau21-(x1*0*1+x2*tau21*F21+x3*tau31*F31)/(x1*1+x2*F21+x3*F31)))+
         (x2*1/(x1*F12+x2*1+x3*F32)*(0-(x1*tau12*F12+x2*0*1+x3*tau32*F32)/(x1*F12+x2*0+x3*F32)))+
         (x3*F23/(x1*F12+x2*F23+x3*1)*(tau23-(x1*tau13*F13+x2*tau23*F23+x3*0*1)/(x1*F13+x2*F23+x3*1)));
ConstStr lngammax3=(x1*tau13*F13+x2*tau23*F23+x3*0*1)/(x1*F13+x2*F23+x3*1)+
         (x1*F31/(x1*1+x2*F21+x3*F31)*(tau31-(x1*0*1+x2*tau21*F21+x3*tau31*F31)/(x1*1+x2*F21+x3*F31)))+
         (x2*F32/(x1*F12+x2*1+x3*F32)*(tau32-(x1*tau12*F12+x2*0*1+x3*tau32*F32)/(x1*F12+x2*0+x3*F32)))+
         (x3*1/(x1*F12+x2*F23+x3*1)*(0-(x1*tau13*F13+x2*tau23*F23+x3*0*1)/(x1*F13+x2*F23+x3*1)));
         
ConstStr gammax1=exp(lngammax1);
ConstStr gammax2=exp(lngammax2);
ConstStr gammax3=exp(lngammax3);

ConstStr lngammay1=(y1*0*1+y2*tau21*F21+y3*tau31*F31)/(y1*1+y2*F21+y3*F31)+
         (y1*1/(y1*1+y2*F21+y3*F31)*(0-(y1*0*1+y2*tau21*F21+y3*tau31*F31)/(y1*1+y2*F21+y3*F31)))+
         (y2*F12/(y1*F12+y2*1+y3*F32)*(tau12-(y1*tau12*F12+y2*0*1+y3*tau32*F32)/(y1*F12+y2*0+y3*F32)))+
         (y3*F13/(y1*F12+y2*F23+y3*1)*(tau13-(y1*tau13*F13+y2*tau23*F23+y3*0*1)/(y1*F13+y2*F23+y3*1)));
ConstStr lngammay2=(y1*tau12*F12+y2*0*1+y3*tau32*F32)/(y1*F12+y2*0+y3*F32)+
         (y1*F21/(y1*1+y2*F21+y3*F31)*(tau21-(y1*0*1+y2*tau21*F21+y3*tau31*F31)/(y1*1+y2*F21+y3*F31)))+
         (y2*1/(y1*F12+y2*1+y3*F32)*(0-(y1*tau12*F12+y2*0*1+y3*tau32*F32)/(y1*F12+y2*0+y3*F32)))+
         (y3*F23/(y1*F12+y2*F23+y3*1)*(tau23-(y1*tau13*F13+y2*tau23*F23+y3*0*1)/(y1*F13+y2*F23+y3*1)));
ConstStr lngammay3=(y1*tau13*F13+y2*tau23*F23+y3*0*1)/(y1*F13+y2*F23+y3*1)+
         (y1*F31/(y1*1+y2*F21+y3*F31)*(tau31-(y1*0*1+y2*tau21*F21+y3*tau31*F31)/(y1*1+y2*F21+y3*F31)))+
         (y2*F32/(y1*F12+y2*1+y3*F32)*(tau32-(y1*tau12*F12+y2*0*1+y3*tau32*F32)/(y1*F12+y2*0+y3*F32)))+
         (y3*1/(y1*F12+y2*F23+y3*1)*(0-(y1*tau13*F13+y2*tau23*F23+y3*0*1)/(y1*F13+y2*F23+y3*1)));
         
ConstStr gammay1=exp(lngammay1);
ConstStr gammay2=exp(lngammay2);
ConstStr gammay3=exp(lngammay3);

Function y1=x1*gammax1/gammay1;
                y2=x2*gammax2/gammay2;
                y3=x3*gammax3/gammay3;

Data;
0.9561 0.0289 0.0150 0.0010 0.0011 0.9979
0.9101 0.0739 0.0160 0.0013 0.0028 0.9959
0.7982 0.1675 0.0343 0.0010 0.0071 0.9919
0.6980 0.2558 0.0462 0.0026 0.0164 0.9810
0.6461 0.2941 0.0598 0.0018 0.0253 0.9729
0.5686 0.3488 0.0826 0.0011 0.0404 0.9585
0.5194 0.3921 0.0885 0.0037 0.0525 0.9438
0.4919 0.4101 0.0980 0.0045 0.0583 0.9372
0.4664 0.4319 0.1017 0.0079 0.0688 0.9233
0.4497 0.4410 0.1093 0.0111 0.0736 0.9153
0.4374 0.4498 0.1128 0.0125 0.0884 0.8991
0.4189 0.4515 0.1296 0.0186 0.0942 0.8872
回复此楼

» 猜你喜欢

» 本主题相关价值贴推荐，对您同样有帮助:
1楼 2017-02-17 15:15:08
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

dingd

铁杆木虫 (职业作家)

【答案】应助回帖

★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★
感谢参与，应助指数 +1
mushixilin: 金币+15, ★★★很有帮助 2017-02-18 09:50:46
sharedmodel;可以实现多因子参数共享拟合，1.5的估计不行。
赞一下 回复此楼

» 本帖已获得的红花（最新10朵）

mushixilin
2楼2017-02-17 17:15:32
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

mushixilin

铜虫 (初入文坛)
送红花一朵
引用回帖:
2楼: Originally posted by dingd at 2017-02-17 17:15:32
sharedmodel;可以实现多因子参数共享拟合，1.5的估计不行。

谢谢，得到大神的解答。在小木虫上看了很多类似的帖子，有提到1.5不可以使用sharedmodel，也有虫友说lib文件夹可以用，下载后使用了一下，发现1.5可以使用sharedmodel功能，可以计算较为简单的sharedmodel问题。但是对于我这个问题，可能由于代码太长，无法计算，提示如下：Fatal: Line too long (more than 1023 characters)，Compile failed, check your program codes please!。
赞一下 回复此楼
3楼2017-02-17 22:53:23
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

mushixilin

铜虫 (初入文坛)
将代码进行了一定处理，仍然提示Fatal: Line too long (more than 1023 characters)，Compile failed, check your program codes please!。是因为最后的公式过于复杂吗？

Title "NRTL simulation";
sharedmodel;
Parameters alpha[0,1],t12,t13,t21,t23,t31,t32;
Variables x1,x2,x3,y1,y2,y3;
ConstStr G12=exp(-alpha*t12),
         G13=exp(-alpha*t13),
         G21=exp(-alpha*t21),
         G23=exp(-alpha*t23),
         G31=exp(-alpha*t31),
         G32=exp(-alpha*t32),
         a=(x2*t21*G21+x3*t31*G31)/(x1+x2*G21+x3*G31),
         b=(x1*t12*G12+x3*t32*G32)/(x1*G12+x2+x3*G32),
         c=(x1*t13*G13+x2*t23*G23)/(x1*G13+x2*G23+x3),
         d=x1/(x1+x2*G21+x3*G31),
         e=x2/(x1*G12+x2+x3*G32),
         f=x3/(x1*G13+x2*G23+x3),
         j=(y2*t21*G21+y3*t31*G31)/(y1+y2*G21+y3*G31),
         k=(y1*t12*G12+y3*t32*G32)/(y1*G12+y2+y3*G32),
         l=(y1*t13*G13+y2*t23*G23)/(y1*G13+y2*G23+y3),
         m=y1/(y1+y2*G21+y3*G31),
         n=y2/(y1*G12+y2+y3*G32),
         o=y3/(y1*G13+y2*G23+y3);
Conststr lngammax1=a-d*a+e*G12*(t12-b)+f*G13*(t13-c),
         lngammax2=b+d*G12*(t21-a)-e*b+f*G23*(t23-c),
         lngammax3=c+d*G31*(t31-a)+e*G32*(t32-b)-f*c,
         lngammay1=j-m*j+n*G12*(t12-k)+o*G13*(t13-l),
         lngammay2=k+m*G12*(t21-j)-n*k+o*G23*(t23-l),
         lngammay3=l+m*G31*(t31-j)+n*G32*(t32-k)-o*l;
Conststr gammax1=exp(lngammax1),
         gammax2=exp(lngammax2),
         gammax3=exp(lngammax3),
         gammay1=exp(lngammay1),
         gammay2=exp(lngammay2),
         gammay3=exp(lngammay3);
Function y1=x1*gammax1/gammay1;
         y2=x2*gammax2/gammay2;
         y3=x3*gammax3/gammay3;
Data;
0.9561 0.0289 0.0150 0.0010 0.0011 0.9979
0.9101 0.0739 0.0160 0.0013 0.0028 0.9959
0.7982 0.1675 0.0343 0.0010 0.0071 0.9919
0.6980 0.2558 0.0462 0.0026 0.0164 0.9810
0.6461 0.2941 0.0598 0.0018 0.0253 0.9729
0.5686 0.3488 0.0826 0.0011 0.0404 0.9585
0.5194 0.3921 0.0885 0.0037 0.0525 0.9438
0.4919 0.4101 0.0980 0.0045 0.0583 0.9372
0.4664 0.4319 0.1017 0.0079 0.0688 0.9233
0.4497 0.4410 0.1093 0.0111 0.0736 0.9153
0.4374 0.4498 0.1128 0.0125 0.0884 0.8991
0.4189 0.4515 0.1296 0.0186 0.0942 0.8872
赞一下 回复此楼
4楼2017-02-18 16:20:57
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

lipenggg

专家顾问 (正式写手)

【答案】应助回帖

★ ★ ★ ★ ★
感谢参与，应助指数 +1
mushixilin: 金币+5, ★★★很有帮助 2017-02-18 22:05:54
1.5版本太低了，即使能计算共享参数拟合，你的代码有点复杂，软件功能不够。
需要正版高版本软件来计算。

[ 发自手机版 http://muchong.com/3g ]
赞一下 回复此楼

» 本帖已获得的红花（最新10朵）

mushixilin
一团和气
5楼2017-02-18 18:09:04
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

mushixilin

铜虫 (初入文坛)
送红花一朵
引用回帖:
5楼: Originally posted by lipenggg at 2017-02-18 18:09:04
1.5版本太低了，即使能计算共享参数拟合，你的代码有点复杂，软件功能不够。
需要正版高版本软件来计算。

谢谢前辈解答。
试了试，字符数多但拟合公式简单的sharedmodel问题可以解决，我这个公式确实有点儿复杂，算不了。
之前在用matlab，可惜对初值要求太大，专程学了1stOpt，没想到还是解不了。我再想想办法。
赞一下 回复此楼
6楼2017-02-18 22:11:33
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

小灵鱼儿

金虫 (小有名气)
引用回帖:
6楼: Originally posted by mushixilin at 2017-02-18 22:11:33
谢谢前辈解答。
试了试，字符数多但拟合公式简单的sharedmodel问题可以解决，我这个公式确实有点儿复杂，算不了。
之前在用matlab，可惜对初值要求太大，专程学了1stOpt，没想到还是解不了。我再想想办法。...

请问你的问题解决了吗？
赞一下 回复此楼
不打无准备之战
7楼2017-03-19 17:47:48
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 mushixilin 的主题更新
71/1返回列表
普通表情 龙 兔 虎 猫
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[基金申请] 欢迎发来filecode的Mz6后的代码验证其规律 +55 医学老男孩 2026-08-13 123/6150 2026-08-19 23:44 by stan_et
[基金申请] 今天基金会出结果吗？20260819 +13 kkkl_v 2026-08-19 13/650 2026-08-19 22:37 by 宝贝虫子
[基金申请] 时间戳变了，能看出什么问题？ +12 基诺咪客 2026-08-17 16/800 2026-08-19 22:04 by kk 的小木虫
[基金申请] filecode，4个jtjc了 +11 ziyangfang 2026-08-19 13/650 2026-08-19 21:28 by aasahr
[教师之家] 为什么余额宝的年化利率越来越低？主要原因有哪些？ +5 瞬息宇宙 2026-08-15 5/250 2026-08-19 20:23 by super2002521
[基金申请] 今天放榜没戏了吧 +4 yuleib84 2026-08-19 5/250 2026-08-19 19:52 by hhs666
[基金申请] 今天放榜吗？ +14 布布和一二 2026-08-19 15/750 2026-08-19 18:07 by gltch
[基金申请] filecode=后面第一个是大写字母 +10 wangze12014 2026-08-14 12/600 2026-08-19 16:56 by 苦难博士
[基金申请] 2027广东省杰青 +4 奶牛小黑 2026-08-15 10/500 2026-08-19 11:02 by wanfengnew
[基金申请] 朋友圈看到的 +6 wangzilk 2026-08-18 8/400 2026-08-19 10:55 by Haru815
[基金申请] 快农历七夕节了，轻松一下，男人悄悄话，女施主请不要进来。 +6 Tide man 2026-08-14 7/350 2026-08-19 09:56 by ZJTJZ
[基金申请] 什么时候开奖？ +6 CrisMessi 2026-08-18 6/300 2026-08-19 08:06 by Equinoxhua
[基金申请] 明天放榜？ +5 Shxjjxjkx 2026-08-18 5/250 2026-08-18 18:14 by -大大大大大-
[论文投稿] 投稿咨询 +4 wwm09 2026-08-17 6/300 2026-08-18 15:36 by wwm09
[基金申请] 今天维护系统维护 祝所有人 高中 +8 gjjjzhong 2026-08-18 9/450 2026-08-18 13:01 by 家与远方
[基金申请] 哪位老哥知道今年的国自然具体哪一天放榜？ +13 Ldrop2023 2026-08-13 16/800 2026-08-18 12:25 by 淀粉搬运工
[基金申请] 93BebMhtakh前后11位开头都是大写 +4 且听虎啸 2026-08-17 5/250 2026-08-18 00:49 by 蔡棒棒菂
[基金申请] 感觉是下周放榜了 +6 angus9576 2026-08-17 11/550 2026-08-17 23:57 by angus9576
[基金申请] 今天系统多次维护，明天很可能放榜！ +8 zju2000 2026-08-16 9/450 2026-08-17 12:20 by lmz0216
[精细化工] 招聘 金属平磨液，抛光液研发工程师 +3 小天0311 2026-08-14 3/150 2026-08-16 07:31 by H9PLUS
信息提示
关闭
请填处理意见
关闭