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

24小时热门版块排行榜    

小木虫论坛-学术科研互动平台 » 计算模拟区 » 计算模拟 » 程序代码 » Matlab 拟合动力学参数时设置的初值不同，输出的结果不同怎么办？
91/1返回列表
查看: 144  |  回复: 8
只看楼主@他人 存档 新回复提醒 (忽略) 收藏    在APP中查看
【悬赏金币】回答本帖问题，作者ZXR319将赠送您 50 个金币

ZXR319

新虫 (小有名气)

[求助] Matlab 拟合动力学参数时设置的初值不同，输出的结果不同怎么办？已有3人参与

我在做 LH 机理动力学参数估计，发现设置的参数初值不同输出的结果就不同，我通过调整参数初值能使 lnk 和1/T 呈线性相关，但 lnK 和1/T 就不成线性了怎么办？我该怎么调整参数初值？以下是我的代码:
% 主脚本
% 步骤1: 读取Excel文件中的数据
filename = 'D:\桌面文件\新60.xlsx';
data = readtable(filename);
time = data{:, 't'}; % 时间数据
CA_exp_333K = data{:, 'CA'}; % 333K下反应物A的浓度数据
CB_exp_333K = data{:, 'CB'}; % 333K下生成物B的浓度数据
% 步骤2: 初始参数估计
params_initial = [0.5,0.5,0.5,0.5,0.5]; % 需要根据实际情况进行调整
lb=[0 0 0 0 0];
ub=[100 100 100 100 100];
% 使用非线性最小二乘法拟合参数
options = optimoptions('lsqnonlin', 'Algorithm', 'levenberg-marquardt', 'Display', 'iter');
params_est = lsqnonlin(@(params) ObjFunc(params, time, CA_exp_333K, CB_exp_333K), params_initial, [], [], options);
% 使用拟合参数重新计算浓度
options = odeset('RelTol', 1e-6, 'AbsTol', 1e-9);
[t, C] = ode45(@(t, C) KineticModel2(t, C, params_est), time, [CA_exp_333K(1); CB_exp_333K(1)], options);
CA_pred = C(:, 1);
CB_pred = C(:, 2);
% 步骤3: 计算总方差和 (SST) 并检查
residuals = ObjFunc(params_est, time, CA_exp_333K, CB_exp_333K);
SST = sum(residuals.^2); % 计算总方差和
% 步骤4: 根据SST值绘制图形或重新设定初值
if SST < 1
% 如果SST小于1，绘制图形并输出参数
fprintf('拟合的参数值为:\n');
fprintf('k: %f\n', params_est(1));
fprintf('K_A: %f\n', params_est(2));
fprintf('K_B: %f\n', params_est(3));
fprintf('K_H: %f\n', params_est(4));
fprintf('K_bc: %f\n', params_est(5));
fprintf('总方差和 (SST): %f\n', SST); % 输出总方差和
% 绘制实验浓度和计算浓度相关图
figure;
plot(time, CA_exp_333K, 'b-', time, CA_pred, 'r--');
hold on;
plot(time, CB_exp_333K, 'g-', time, CB_pred, 'k--');
legend('实验CA', '拟合CA', '实验CB', '拟合CB');
xlabel('时间 (min)');
ylabel('浓度 (mol/L)');
title('实验浓度与计算浓度对比');
hold off;
else
% 如果SST大于或等于1，提示需要重新设定初值
fprintf('SST值过大，需要重新设定参数初值。当前SST: %f\n', SST);
% 这里可以添加代码来调整初始参数或采取其他措施
end
function dCdt = KineticModel2(t, C, params)
k = params(1);
K_A = params(2);
K_B = params(3);
K_H = params(4);
K_bc = params(5);
P_H2 = 3e3; % 氢气压力，单位为KPa
c_A = C(1);
c_B = C(2);
numerator = k * K_A * c_A * sqrt((K_H * P_H2) / K_bc);
denominator = (1 + K_A * c_A + K_B * c_B+ sqrt((K_H * P_H2) / K_bc))^2;
r = numerator / denominator;
dC_A_dt = -r;
dC_B_dt = r;
dCdt = [dC_A_dt; dC_B_dt];
end
function res = ObjFunc(params, time, CA_exp, CB_exp)
options = odeset('RelTol', 1e-6, 'AbsTol', 1e-9);
[t, C] = ode45(@(t, C) KineticModel2(t, C, params), time, [CA_exp(1); CB_exp(1)], options);
CA_pred = C(:, 1);
CB_pred = C(:, 2);
% 计算总方差和
res = (CA_exp - CA_pred).^2 + (CB_exp - CB_pred).^2;
end
以下图片中是我的数据，跪求大佬指导?
参数初值通过查文献给出的方法不现实，因为我做这个没有人做过，类似的文献初值差距也很大，我觉得没有什么参考价值???而且我不太明白，比如我有五个参数 k、KA、KB、KH、Kbc，那这五个参数初值应该一样吗？

Matlab 拟合动力学参数时设置的初值不同，输出的结果不同怎么办？


发自小木虫Android客户端
回复此楼

» 猜你喜欢
1楼2024-10-06 15:38:04
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

xrs333

新虫 (正式写手)

【答案】应助回帖

感谢参与，应助指数 +1
你的方程有多个根，比如，0.1 - x + 2*x^3 = 0，有三个根。使用不同的迭代算法以及给定不同预估值，可能会收敛到不同的根，需要提供适当的预估值（迭代初值）才能迭代求得所要求的根。
赞一下 回复此楼
2楼2024-10-07 09:10:05
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

johnliu1983

至尊木虫 (著名写手)

【答案】应助回帖

感谢参与，应助指数 +1
matlab在拟合是时候，你应该可以限制参数范围，先找到你认为比较合理的参数范围，然后把设置参数的范围设置好就行。
赞一下 回复此楼
3楼2024-10-08 18:39:53
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

ZXR319

新虫 (小有名气)
引用回帖:
2楼: Originally posted by xrs333 at 2024-10-07 09:10:05
你的方程有多个根，比如，0.1 - x + 2*x^3 = 0，有三个根。使用不同的迭代算法以及给定不同预估值，可能会收敛到不同的根，需要提供适当的预估值（迭代初值）才能迭代求得所要求的根。

谢谢您的回复！关键问题是我不知道合适的参数初值应该是多少?我做的东西别人没做过，类似的文献中的参数带进去拟合效果又不好，所以我应该考虑全局最优化算法吗？

发自小木虫Android客户端
赞一下 回复此楼
4楼2024-10-09 11:47:21
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

ZXR319

新虫 (小有名气)
引用回帖:
3楼: Originally posted by johnliu1983 at 2024-10-08 18:39:53
matlab在拟合是时候，你应该可以限制参数范围，先找到你认为比较合理的参数范围，然后把设置参数的范围设置好就行。

谢谢您的回复！你说的意思我明白，可我参数初值代0.5或1的时候得到的结果都不一样，这个参数范围实在是不好确定，而且就算确定了范围，给出一个初值，拟合效果也不好，不知道为什么，真的愁死了

发自小木虫Android客户端
赞一下 回复此楼
5楼2024-10-09 11:50:23
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

xrs333

新虫 (正式写手)
引用回帖:
4楼: Originally posted by ZXR319 at 2024-10-09 11:47:21
谢谢您的回复！关键问题是我不知道合适的参数初值应该是多少?我做的东西别人没做过，类似的文献中的参数带进去拟合效果又不好，所以我应该考虑全局最优化算法吗？
...

可以尝试，把方程增加一个变量b，然后，在给定状态下，计算一系列待定系数时该变量的值，找到使变量b为零（或接近为零）的待定系数。
赞一下 回复此楼
6楼2024-10-09 18:12:01
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

johnliu1983

至尊木虫 (著名写手)

【答案】应助回帖

引用回帖:
5楼: Originally posted by ZXR319 at 2024-10-09 11:50:23
谢谢您的回复！你说的意思我明白，可我参数初值代0.5或1的时候得到的结果都不一样，这个参数范围实在是不好确定，而且就算确定了范围，给出一个初值，拟合效果也不好，不知道为什么，真的愁死了
...

要不你把数据和你的m文件发给我，我有空帮你看看。 johnliu1983@126.com
赞一下 回复此楼
7楼2024-10-09 20:17:28
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

hzlhm

至尊木虫 (著名写手)

【答案】应助回帖

1、先用rand函数作为初值
2、然后修正初值
3、最后得到结果
链接: https://pan.baidu.com/s/1X3gWCp4EEbFsWrwCKLyeEA?pwd=wur3 提取码: wur3 复制这段内容后打开百度网盘手机App，操作更方便哦
赞一下 回复此楼
QQ：2120156492
8楼2024-10-09 21:36:21
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

hzlhm

至尊木虫 (著名写手)

【答案】应助回帖

CA
R^2=  0.98593573
CB
R^2=  0.98593573
333℃时的系数：
k：  0.45591792
K_A：  0.45589602
K_B：  0.47448748
K_H：  0.50329217
K_bc：  0.45984830
赞一下 回复此楼
QQ：2120156492
9楼2024-10-09 21:37:30
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 ZXR319 的主题更新
91/1返回列表
不应助 确定回帖应助 (注意：应助才可能被奖励，但不允许灌水，必须填写15个字符以上)
普通表情 龙 兔 虎 猫
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[教师之家] 女博士寻找有缘人 (EPI+-1)(金币-50) +7 AlishaChen 2024-10-09 8/400 2024-10-09 23:34 by mch20120521
[考博] 环境或者土木专业想读博的学生可以私聊我 +5 zmr8 2024-10-08 9/450 2024-10-09 22:44 by wdcy
[论文投稿] 论文 +3 1821588220 2024-10-08 6/300 2024-10-09 21:56 by 1821588220
[基金申请] 为啥诺贝尔奖不要填表盖章，不要答辩，不公布评审人名单 +17 babu2015 2024-10-08 17/850 2024-10-09 19:39 by asdew124
[基金申请] 这几天，有院士等瑞典的电话吗？ +12 babu2015 2024-10-07 13/650 2024-10-09 19:10 by zjctec
[基金申请] 博后面上，求中！！！ +9 近未来0814 2024-10-08 11/550 2024-10-09 16:47 by 近未来0814
[教师之家] 严重违反诚信行为！复旦大学公布一批学生名单 +6 苏东坡二世 2024-10-09 6/300 2024-10-09 16:18 by ma0526
[论文投稿] science china material投稿 50+5 bingbinglove 2024-10-05 8/400 2024-10-09 12:14 by ￥笑傲江湖￥
[教师之家] 老师上课迟到5分钟算事故，领导让一群下属等半小时以上 +12 河西夜郎 2024-10-04 14/700 2024-10-09 11:50 by 小钱景K
[基金申请] 青年基金中了没立项，能否再次申请？ +10 leicahfs 2024-10-07 17/850 2024-10-09 10:52 by caiqingchong
[基金申请] 谢美蓝与沈磊离婚，和路杰结婚，终成人生赢家 +9 苏东坡二世 2024-10-05 9/450 2024-10-09 10:41 by forever411
[有机交流] 酸溶于水萃取不出来 5+3 遥遥万里. 2024-10-08 3/150 2024-10-09 10:26 by 张国馗1
[论文投稿] JACS投稿经验求助 20+3 mmjun1015 2024-10-06 10/500 2024-10-08 15:19 by HaLaFulan
[论文投稿] 投稿求助 5+3 科研@ 2024-10-03 5/250 2024-10-08 14:59 by 西楚淮左鲁北
[论文投稿] 投稿英文期刊中有中文 10+3 wxj824 2024-10-05 5/250 2024-10-08 11:56 by 小李123123
[海外博后] 已入职编制副教授，可以出国博后吗 +3 blake1111 2024-10-05 4/200 2024-10-06 05:52 by zhangysbad
[论文投稿] 三审论文被拒，已申诉，等待结果中。 +7 sadhioaw 2024-10-04 9/450 2024-10-05 16:21 by iwdw2012
[论文投稿] ECM能源期刊 5+4 ！学习爱我！ 2024-10-04 4/200 2024-10-05 16:14 by babero
[论文投稿] 关于爱思唯尔proof central 中 corresponding footnote and email address 的问题？ 20+3 maoxiao 2024-10-05 6/300 2024-10-05 16:09 by babero
[教师之家] 国自然结题问题咨询 +3 lyz2003cd 2024-10-05 4/200 2024-10-05 10:21 by lyz2003cd
信息提示
关闭
请填处理意见
关闭