24小时热门版块排行榜    

查看: 2523  |  回复: 12
【悬赏金币】回答本帖问题,作者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客户端
回复此楼
已阅   回复此楼   关注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

» 本帖已获得的红花(最新10朵)

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,操作更方便哦

» 本帖已获得的红花(最新10朵)

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

新虫 (小有名气)

送红花一朵
引用回帖:
7楼: Originally posted by johnliu1983 at 2024-10-09 20:17:28
要不你把数据和你的m文件发给我,我有空帮你看看。 johnliu1983@126.com...

太感谢您了,邮件已发,请查收

发自小木虫Android客户端
10楼2024-10-10 11:14:27
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 ZXR319 的主题更新
不应助 确定回帖应助 (注意:应助才可能被奖励,但不允许灌水,必须填写15个字符以上)
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[考博] 26申博求助 +3 跳跃饼干 2026-03-10 4/200 2026-03-10 21:15 by Tntcnn
[考研] 求调剂 +4 yfihxh 2026-03-09 4/200 2026-03-10 16:50 by Pegasus0429
[考研] 材料与化工求调剂 +11 与冬清宁 2026-03-07 14/700 2026-03-10 14:23 by Linda Hu
[基金申请] PDF在线压缩软件 10+3 haxia 2026-03-08 4/200 2026-03-10 13:44 by lfy8008
[考研] 0856求调剂 +3 squirtle11 2026-03-07 3/150 2026-03-09 09:54 by @飒飒飒飒
[考研] 一志愿武理314求调剂 +4 ( ̄~ ̄;) 2026-03-08 5/250 2026-03-08 23:12 by L135790
[考研] 347求调剂 +4 浮云满足 2026-03-07 4/200 2026-03-08 16:46 by 星空星月
[考研] 085701环境工程317分求调剂 +9 6汆尼9 2026-03-07 9/450 2026-03-08 06:41 by 刘兵
[考研] 081700学硕一志愿北京化工大学数二英一过六级有竞赛求调剂 +5 galaxary 2026-03-07 7/350 2026-03-08 04:37 by wang_dand
[考研] 301求调剂 +5 一二LV 2026-03-07 5/250 2026-03-07 22:20 by 18137688336
[考研] 2026调剂】考试A区0703化学类323分 诚求接收 +3 卷柏卷柏 2026-03-05 4/200 2026-03-07 17:00 by yuandd_2001
[考研] 复试调剂 +7 呼呼?~+123456 2026-03-05 10/500 2026-03-06 16:07 by @飒飒飒飒
[考研] 085600材料与化工(高分子)290分求调剂 +9 wengyujian 2026-03-04 9/450 2026-03-06 15:59 by 聪明的大松鼠
[考研] 289求调剂 +3 驰光绝景 2026-03-04 4/200 2026-03-06 12:34 by 驰光绝景
[考研] 085600,一志愿郑州大学,280分求调剂 +7 Wuqi725 2026-03-05 7/350 2026-03-06 08:12 by 杨杨杨紫
[考研] 求调剂 +3 泡了个椒 2026-03-04 4/200 2026-03-05 14:37 by 泡了个椒
[考研] 376求调剂 +3 王浩然sam 2026-03-04 3/150 2026-03-05 07:48 by bxbo
[考研] 材料学硕080500复试调剂294 +3 四叶zjz 2026-03-04 3/150 2026-03-05 07:09 by kunm555
[考研] 286 +6 ksncj 2026-03-04 6/300 2026-03-04 20:49 by lature00
[考研] 0703化学求调剂 +3 将未有期 2026-03-04 3/150 2026-03-04 16:19 by sslc1985
信息提示
请填处理意见