| 查看: 3117 | 回复: 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客户端 |
回复此楼» 猜你喜欢
26申博推荐：南京航空航天大学国际前沿院光学方向招收博士生！
已经有1人回复
-
如何从铁电相到顺电相。
已经有1人回复
-
物理学I论文润色/翻译怎么收费?
已经有272人回复
-
各位大佬，求一份最新Wien2k版本作为学习用途
已经有0人回复
-
南科大活性流体和软物质课题组诚招2027级博士、硕士研究生和博后
已经有40人回复
-
瑞典林雪平大学博士后招聘｜PEC Water Splitting 方向
已经有24人回复
-
PRB投稿，Acknowledgment sent to author状态十天了
已经有3人回复
-
PVA溶解
已经有1人回复
-
基金申请
已经有44人回复
xrs333
木虫 (正式写手)
- 应助: 151 (高中生)
- 金币: 5117.2
- 红花: 23
- 帖子: 390
- 在线: 239.5小时
- 虫号: 1867708
- 注册: 2012-06-23
- 性别: GG
- 专业: 内流流体力学
2楼2024-10-07 09:10:05
johnliu1983
至尊木虫 (著名写手)
- 应助: 43 (小学生)
- 金币: 15180.3
- 红花: 8
- 帖子: 2459
- 在线: 150.5小时
- 虫号: 370817
- 注册: 2007-05-14
- 性别: GG
- 专业: 物理学II
3楼2024-10-08 18:39:53
|
谢谢您的回复！关键问题是我不知道合适的参数初值应该是多少?我做的东西别人没做过，类似的文献中的参数带进去拟合效果又不好，所以我应该考虑全局最优化算法吗？ 发自小木虫Android客户端 |
4楼2024-10-09 11:47:21
|
谢谢您的回复！你说的意思我明白，可我参数初值代0.5或1的时候得到的结果都不一样，这个参数范围实在是不好确定，而且就算确定了范围，给出一个初值，拟合效果也不好，不知道为什么，真的愁死了 发自小木虫Android客户端 |
5楼2024-10-09 11:50:23
xrs333
木虫 (正式写手)
- 应助: 151 (高中生)
- 金币: 5117.2
- 红花: 23
- 帖子: 390
- 在线: 239.5小时
- 虫号: 1867708
- 注册: 2012-06-23
- 性别: GG
- 专业: 内流流体力学
6楼2024-10-09 18:12:01
johnliu1983
至尊木虫 (著名写手)
- 应助: 43 (小学生)
- 金币: 15180.3
- 红花: 8
- 帖子: 2459
- 在线: 150.5小时
- 虫号: 370817
- 注册: 2007-05-14
- 性别: GG
- 专业: 物理学II
ZXR3197楼2024-10-09 20:17:28
【答案】应助回帖
|
1、先用rand函数作为初值 2、然后修正初值 3、最后得到结果 链接: https://pan.baidu.com/s/1X3gWCp4EEbFsWrwCKLyeEA?pwd=wur3 提取码: wur3 复制这段内容后打开百度网盘手机App，操作更方便哦 |
» 本帖已获得的红花（最新10朵）
8楼2024-10-09 21:36:21
9楼2024-10-09 21:37:30
10楼2024-10-10 11:14:27
15