基于已知与未知变量求解幂律模型参数并优化拟合精度

需求说明 #

• 仅已知1个k值和3个E值，需猜测另外两个k值的上下界
• 求解幂律模型参数，要求已知k值与对应E值下模型计算的k值差值为0（当前差值为0.002）
• 通过类似Excel规划求解的逻辑，调整k的猜测值来消除差值

原代码（已翻译中文注释） #

% ------------------ 幂律模型分析 ------------------ %
E1 = Sheet.("ED (W/m3)");  % 自变量（已知）

% 提取E1列的前3个值
E = E1(1:3);

% 对应E的k值（已知或先前计算的）
% 替换为实际k值
k = [0.005, k_zero_order, 0.009]; 

% 对E和k取对数
logE = log(E);
logk = log(k);

% 线性回归求alpha和log(C)
p = polyfit(logE, logk, 1);
alpha = p(1);   % 拟合直线的斜率，即alpha
logC = p(2);    % 拟合直线的截距，即log(C)

% 由log(C)计算C
C = exp(logC);

% 显示幂律常数结果
 fprintf('回归分析结果:\n');
 fprintf('估计的幂律指数alpha: %.4f\n', alpha);
 fprintf('估计的系数C: %.4f\n', C);

% 用估计的幂律模型计算k值
k_estimated = C * E.^alpha;

% 用已知k值验证模型
% 替换为实际已知k值和对应E值
known_k_value = k_zero_order; % 用于验证的已知k值

% 对应已知k的E值
known_E_value = E(2); % 如果不同，替换为实际对应E值

% 计算已知E值对应的估计k值
k_est_known = C * known_E_value^alpha;

% 显示验证结果
fprintf('使用已知k = %.4f 和对应E = %.4f:\n', known_k_value, known_E_value);
fprintf('已知E对应的估计k值: %.4f\n', k_est_known);
fprintf('估计k值与已知k值的差值: %.4f\n', abs(k_est_known - known_k_value));

修改后的代码（实现参数优化） #

% ------------------ 幂律模型参数优化求解 ------------------ %
E1 = Sheet.("ED (W/m3)");  % 自变量（已知）
E = E1(1:3);               % 提取前3个E值
known_k_idx = 2;           % 已知k值对应的索引（此处为第2个）
known_k_value = k_zero_order; % 已知的k值

% 定义目标函数：输入猜测的k1和k3，返回已知k值的预测误差平方
function error_sq = objective_func(guess_k, E, known_k_idx, known_k_value)
    k = [guess_k(1), known_k_value, guess_k(2)];
    logE = log(E);
    logk = log(k);
    p = polyfit(logE, logk, 1);
    alpha = p(1);
    C = exp(p(2));
    % 计算已知E对应的预测k值
    k_pred = C * E(known_k_idx)^alpha;
    % 返回误差平方，用于最小化
    error_sq = (k_pred - known_k_value)^2;
end

% 初始猜测的两个未知k值（可根据实际情况调整）
initial_guess = [0.004, 0.010]; 
% 设置k值的上下界（根据物理意义调整）
lb = [0.001, 0.008]; % 下界
ub = [0.006, 0.012]; % 上界

% 使用带边界约束的优化函数fminsearchbnd（若无该函数，可替换为fminsearch）
optimized_k = fminsearchbnd(@(x) objective_func(x, E, known_k_idx, known_k_value), initial_guess, lb, ub);

% 获取优化后的完整k数组
k_opt = [optimized_k(1), known_k_value, optimized_k(2)];

% 重新拟合幂律模型
logE_opt = log(E);
logk_opt = log(k_opt);
p_opt = polyfit(logE_opt, logk_opt, 1);
alpha_opt = p_opt(1);
C_opt = exp(p_opt(2));

% 计算最终误差
k_pred_known = C_opt * E(known_k_idx)^alpha_opt;
error_final = abs(k_pred_known - known_k_value);

% 输出结果
fprintf('优化后的未知k值:\n');
fprintf('k1 = %.6f, k3 = %.6f\n', optimized_k(1), optimized_k(2));
fprintf('优化后的幂律参数:\n');
fprintf('alpha = %.6f, C = %.6f\n', alpha_opt, C_opt);
fprintf('已知k值与模型预测值的差值: %.6f\n', error_final);

% 输出所有E对应的估计k值
k_estimated_opt = C_opt * E.^alpha_opt;
fprintf('所有E对应的估计k值:\n');
for i = 1:length(E)
    fprintf('E%d = %.4f, 估计k值 = %.6f\n', i, E(i), k_estimated_opt(i));
end

代码说明 #

• 目标函数：以两个未知k值为输入，拟合幂律模型后计算已知k值的预测误差平方，将其作为优化目标（最小化该值）
• 优化逻辑：通过fminsearchbnd（带边界约束）或fminsearch（无约束）寻找最优未知k值，确保已知k的预测误差趋近于0
• 边界控制：可根据实际物理场景调整lb和ub，避免出现不符合逻辑的k值
• 结果验证：优化完成后重新拟合模型，输出最终参数与误差值，确认差值满足要求

内容的提问来源于stack exchange，提问作者Jacobus02