| 查看: 2074 | 回复: 12 | ||
| 【悬赏金币】回答本帖问题，作者ZXR319将赠送您 50 个金币 | ||
| 当前只显示满足指定条件的回帖，点击这里查看本话题的所有回帖 | ||
[求助]
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，那这五个参数初值应该一样吗？  发自小木虫Android客户端 |
回复此楼» 猜你喜欢
比利时鲁文大学固体物理系招收CSC博士
已经有0人回复
-
🌟 比利时新鲁汶大学（UCLouvain）诚邀CSC博士加入Pascal Gehring教授团队
已经有0人回复
-
物理学I论文润色/翻译怎么收费?
已经有202人回复
-
第一性原理计算方向2026级博士申请 PRB*1，四级484
已经有1人回复
-
求助VISSIM破解版软件
已经有0人回复
-
求2026年在台湾举行的物理和材料领域国际学术会议信息
已经有0人回复
-
求国际会议网站
已经有1人回复
-
求取一些关于纳米材料和纳米技术相关的英文PPT。
已经有0人回复
-
【复旦大学】二维材料方向招收2026年博士研究生1名
已经有0人回复
-
北京纳米能源与系统研究所 王中林院士/曹南颖研究员课题组2026级硕/博/博后招生
已经有10人回复
送红花一朵
11楼2024-10-10 11:16:44
xrs333
木虫 (正式写手)
- 应助: 151 (高中生)
- 金币: 5117.2
- 红花: 23
- 帖子: 390
- 在线: 239.5小时
- 虫号: 1867708
- 注册: 2012-06-23
- 性别: GG
- 专业: 内流流体力学
2楼2024-10-07 09:10:05
johnliu1983
至尊木虫 (著名写手)
- 应助: 43 (小学生)
- 金币: 15079.3
- 红花: 8
- 帖子: 2424
- 在线: 149.7小时
- 虫号: 370817
- 注册: 2007-05-14
- 性别: GG
- 专业: 物理学II
3楼2024-10-08 18:39:53
|
谢谢您的回复！关键问题是我不知道合适的参数初值应该是多少?我做的东西别人没做过，类似的文献中的参数带进去拟合效果又不好，所以我应该考虑全局最优化算法吗？ 发自小木虫Android客户端 |
4楼2024-10-09 11:47:21