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

24小时热门版块排行榜    

小木虫论坛-学术科研互动平台 » 计算模拟区 » 计算模拟 » 程序代码 » Matlab 拟合动力学参数时设置的初值不同，输出的结果不同怎么办？
91/1返回列表
查看: 132  |  回复: 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个字符以上)
普通表情 龙 兔 虎 猫
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[论文投稿] 论文 +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
[论文投稿] 在线等，感谢小木虫，感谢虫友们了哈， +3 sdgame0 2024-10-07 10/500 2024-10-09 16:25 by nono2009
[教师之家] 严重违反诚信行为！复旦大学公布一批学生名单 +6 苏东坡二世 2024-10-09 6/300 2024-10-09 16:18 by ma0526
[考博] 课题组有学博名额！！ +4 宋宋要上岸。 2024-10-08 4/200 2024-10-09 11:05 by 岳牛仔
[基金申请] 等待基金委接受纸质 +9 阿呆不呆 2024-10-08 15/750 2024-10-09 11:03 by 阿呆不呆
[基金申请] 谢美蓝与沈磊离婚，和路杰结婚，终成人生赢家 +9 苏东坡二世 2024-10-05 9/450 2024-10-09 10:41 by forever411
[考博] 25申博 +5 sssvvv 2024-10-07 5/250 2024-10-09 10:27 by August1351
[有机交流] 酸溶于水萃取不出来 5+3 遥遥万里. 2024-10-08 3/150 2024-10-09 10:26 by 张国馗1
[硕博家园] 高校教职offer选择 +19 orangex2 2024-10-02 26/1300 2024-10-09 08:32 by 温柔的Rain
[考博] 东北大学新能源转换与存储团队诚招2025级化学、化工、材料与冶金方向博士 +3 实验一切顺利102 2024-10-06 3/150 2024-10-09 01:20 by 学员Eg9jRO
[基金申请] 在国家经费有限的情况下，之前申请面上两年停1年确实资助率不至于今年这么低 +6 苏东坡二世 2024-10-05 6/300 2024-10-09 01:03 by sujianquan
[论文投稿] Accept!!! +5 jibeidaxue 2024-10-05 6/300 2024-10-08 15:38 by hzwdzdf
[论文投稿] JACS投稿经验求助 20+3 mmjun1015 2024-10-06 10/500 2024-10-08 15:19 by HaLaFulan
[基金申请] 2024年广东省面上基金什么时候公布结果呀，区域联合基金还有么 +3 微风阵阵 2024-10-06 3/150 2024-10-07 22:39 by a441914426
[论文投稿] IOP系统投稿几天发现稿件ID无效了 是不是代表拒稿了 +3 Z780283717 2024-10-04 4/200 2024-10-07 12:57 by Z780283717
[论文投稿] Justification for Accelerated Publication期刊要求写这个 怎么填 5+4 tdrwx 2024-10-05 6/300 2024-10-06 22:31 by tdrwx
[论文投稿] 三审论文被拒，已申诉，等待结果中。 +7 sadhioaw 2024-10-04 9/450 2024-10-05 16:21 by iwdw2012
[论文投稿] 关于爱思唯尔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
信息提示
关闭
请填处理意见
关闭