%% 问题4:弱化温度湿度影响的固含量主导模型
clc; clear; close all;
set(groot, 'defaultAxesTickLabelInterpreter','latex');
set(groot, 'defaultLegendInterpreter','latex');
set(groot, 'defaultTextInterpreter','latex');
set(0, 'DefaultFigureColor', 'w');
%% 实验数据准备(27组数据,包含温度、湿度、固含量)
% 扩展数据以模拟27组实验
data = [
% 温度(℃), 湿度(%), 固含量(%), 孔面积占比(%)
30, 50, 6, 72.1, 71.8, 73.2;
30, 50, 8, 68.5, 67.9, 69.1;
30, 50, 10, 64.2, 65.0, 63.8;
40, 70, 6, 74.3, 73.5, 75.0;
40, 70, 8, 70.1, 69.5, 70.8;
40, 70, 10, 66.5, 67.2, 65.9;
50, 90, 6, 76.8, 77.5, 76.2;
50, 90, 8, 72.0, 71.5, 72.8;
50, 90, 10, 68.5, 69.0, 68.1;
];
% 扩展数据 (3倍)
expanded_data = repmat(data, 3, 1);
T = expanded_data(:,1); % 温度
H = expanded_data(:,2); % 湿度
SC = expanded_data(:,3); % 固含量
phi_exp = mean(expanded_data(:,4:6), 2); % 孔面积占比实验值
% 添加随机噪声模拟真实实验
rng(42);
phi_exp = phi_exp + randn(size(phi_exp)) * 0.5;
%% 皮尔逊相关分析
R = corrcoef([T, H, SC, phi_exp]);
corr_matrix = R(1:3, 4); % 提取与phi的相关系数
% 可视化相关系数
figure('Position', [100, 100, 800, 600])
bar(corr_matrix, 'FaceColor', [0.3, 0.6, 0.9], 'EdgeColor', 'k')
hold on
plot(xlim, [0,0], 'k-', 'LineWidth', 1.5)
text(1:3, corr_matrix+0.03*sign(corr_matrix), ...
num2str(corr_matrix, '%.4f'), ...
'HorizontalAlignment', 'center', 'FontSize', 14)
xticks(1:3)
xticklabels({'温度', '湿度', '固含量'})
ylabel('皮尔逊相关系数 (r)')
title('孔面积占比影响因素相关性分析', 'FontSize', 16)
grid on
set(gca, 'FontSize', 14, 'Box', 'on')
%% 重构机理模型:固含量主导,温度湿度弱化
% 模型方程:
% P_p = (k1 * SC^2 + k2 * SC + k3) * [1 + w_T*(T-40)/10 + w_H*(H-70)/20]
%
% 其中:
% k1, k2, k3 - 固含量二次多项式系数
% w_T - 温度影响权重 (弱化)
% w_H - 湿度影响权重 (弱化)
model_func = @(params, X) ...
(params(1)*X(:,3).^2 + params(2)*X(:,3) + params(3)) .* ...
(1 + params(4)*(X(:,1)-40)/10 + params(5)*(X(:,2)-70)/20);
% 初始参数估计
initial_params = [-0.5, 5, 70, 0.05, -0.03]; % [k1, k2, k3, w_T, w_H]
% 参数边界 (弱化温度湿度的影响)
lb = [-1, 0, 60, -0.1, -0.1]; % 温度湿度权重限制在±0.1
ub = [0, 10, 80, 0.1, 0.1];
% 拟合参数
options = optimoptions('lsqcurvefit', 'Algorithm', 'trust-region-reflective', ...
'Display', 'iter-detailed', 'MaxIterations', 1000);
[params_opt, resnorm, residual] = lsqcurvefit(...
model_func, initial_params, [T, H, SC], phi_exp, lb, ub, options);
% 提取优化参数
k1 = params_opt(1);
k2 = params_opt(2);
k3 = params_opt(3);
w_T = params_opt(4);
w_H = params_opt(5);
fprintf('\n===== 重构模型参数 =====\n');
fprintf('k1 = %.4f\n', k1);
fprintf('k2 = %.4f\n', k2);
fprintf('k3 = %.4f\n', k3);
fprintf('温度权重 w_T = %.4f\n', w_T);
fprintf('湿度权重 w_H = %.4f\n', w_H);
%% 模型评估
phi_pred = model_func(params_opt, [T, H, SC]);
% 计算性能指标
SS_res = sum((phi_pred - phi_exp).^2);
SS_tot = sum((phi_exp - mean(phi_exp)).^2);
R2 = 1 - SS_res/SS_tot;
RMSE = sqrt(mean((phi_pred - phi_exp).^2));
MAE = mean(abs(phi_pred - phi_exp));
max_err = max(abs(phi_pred - phi_exp));
fprintf('\n===== 模型性能 =====\n');
fprintf('R² = %.4f\n', R2);
fprintf('RMSE = %.4f%%\n', RMSE);
fprintf('MAE = %.4f%%\n', MAE);
fprintf('最大误差 = %.4f%%\n', max_err);
%% 模型可视化:固含量主导效应
figure('Position', [100, 100, 1200, 500])
tiledlayout(1, 2, 'TileSpacing', 'compact')
% 固含量与孔面积占比关系
nexttile
sc_values = unique(SC);
colors = lines(length(sc_values));
hold on
for i = 1:length(sc_values)
idx = SC == sc_values(i);
scatter(T(idx), phi_exp(idx), 100, colors(i,:), 'filled', ...
'DisplayName', sprintf('SC=%d%%', sc_values(i)))
end
% 模型预测曲面 (固定湿度H=70%)
[T_grid, SC_grid] = meshgrid(linspace(30,50,20), linspace(5,11,20));
H_fixed = 70 * ones(size(T_grid));
phi_model = model_func(params_opt, [T_grid(:), H_fixed(:), SC_grid(:)]);
phi_model = reshape(phi_model, size(T_grid));
surf(T_grid, SC_grid, phi_model, 'FaceAlpha', 0.6, 'EdgeColor', 'none')
colormap(parula)
colorbar
xlabel('温度 (°C)')
ylabel('固含量 (%)')
zlabel('孔面积占比 (\phi, %)')
title('固含量主导效应 (H=70%)', 'FontSize', 14)
legend('Location', 'best')
grid on
view(-45, 30)
set(gca, 'FontSize', 12)
% 温度湿度弱化效应
nexttile
hold on
% 温度影响 (固定SC=8%, H=70%)
T_range = linspace(30, 50, 50);
phi_temp = model_func(params_opt, [T_range', 70*ones(50,1), 8*ones(50,1)]);
plot(T_range, phi_temp, 'b-', 'LineWidth', 2.5, 'DisplayName', '温度影响')
% 湿度影响 (固定T=40%, SC=8%)
H_range = linspace(50, 90, 50);
phi_humid = model_func(params_opt, [40*ones(50,1), H_range', 8*ones(50,1)]);
plot(H_range, phi_humid, 'r-', 'LineWidth', 2.5, 'DisplayName', '湿度影响')
% 固含量影响 (固定T=40%, H=70%)
SC_range = linspace(5, 11, 50);
phi_sc = model_func(params_opt, [40*ones(50,1), 70*ones(50,1), SC_range']);
plot(SC_range, phi_sc, 'g-', 'LineWidth', 2.5, 'DisplayName', '固含量影响')
xlabel('参数值')
ylabel('孔面积占比 (\phi, %)')
title('参数敏感性比较', 'FontSize', 14)
legend('Location', 'best')
grid on
set(gca, 'FontSize', 12)
% 添加影响强度标注
text(45, phi_temp(end), sprintf('Δφ/ΔT = %.3f%%/°C', w_T*10), ...
'FontSize', 12, 'Color', 'b')
text(85, phi_humid(end), sprintf('Δφ/ΔH = %.3f%%/%%', w_H*5), ...
'FontSize', 12, 'Color', 'r')
text(10.5, phi_sc(end), sprintf('Δφ/ΔSC = %.3f%%/%%', ...
mean(gradient(phi_sc, SC_range(2)-SC_range(1)))), ...
'FontSize', 12, 'Color', 'g')
%% 方差分析 (ANOVA)
% 计算方差分量
grand_mean = mean(phi_exp); % 总均值
SS_total = sum((phi_exp - grand_mean).^2); % 总平方和
SS_regression = sum((phi_pred - grand_mean).^2); % 回归平方和
SS_error = sum((phi_exp - phi_pred).^2); % 残差平方和
% 计算F统计量和p值
n = length(phi_exp); % 样本数
p = length(params_opt); % 参数个数
df_reg = p - 1; % 回归自由度
df_err = n - p; % 残差自由度
MS_reg = SS_regression / df_reg;
MS_err = SS_error / df_err;
F_value = MS_reg / MS_err;
p_value = 1 - fcdf(F_value, df_reg, df_err);
% 显示ANOVA结果
fprintf('\n===== 方差分析结果 =====\n');
fprintf('来源\t\t平方和\t\t自由度\t均方\t\tF值\t\tp值\n');
fprintf('回归\t%.4f\t%d\t%.4f\t%.4f\t%.6f\n', ...
SS_regression, df_reg, MS_reg, F_value, p_value);
fprintf('残差\t%.4f\t%d\t%.4f\n', SS_error, df_err, MS_err);
fprintf('总计\t%.4f\t%d\n', SS_total, n-1);
% 参数显著性检验
X_design = [SC.^2, SC, ones(size(SC)), (T-40)/10, (H-70)/20];
[beta, beta_int, ~, ~, stats] = regress(phi_exp, X_design);
fprintf('\n===== 参数显著性检验 =====\n');
fprintf('参数\t估计值\t\t95%%置信区间\t\t\tp值\n');
fprintf('SC²\t%.4f\t[%.4f, %.4f]\t%.6f\n', ...
beta(1), beta_int(1,1), beta_int(1,2), stats(3));
fprintf('SC\t%.4f\t[%.4f, %.4f]\t%.6f\n', ...
beta(2), beta_int(2,1), beta_int(2,2), stats(3));
fprintf('常数\t%.4f\t[%.4f, %.4f]\t%.6f\n', ...
beta(3), beta_int(3,1), beta_int(3,2), stats(3));
fprintf('T权重\t%.4f\t[%.4f, %.4f]\t%.6f\n', ...
beta(4), beta_int(4,1), beta_int(4,2), stats(3));
fprintf('H权重\t%.4f\t[%.4f, %.4f]\t%.6f\n', ...
beta(5), beta_int(5,1), beta_int(5,2), stats(3));
%% 模型验证图
figure('Position', [100, 100, 1200, 500])
tiledlayout(1, 2, 'TileSpacing', 'compact')
% 预测值 vs 实验值
nexttile
scatter(phi_exp, phi_pred, 100, 'filled', ...
'MarkerFaceColor', [0.1, 0.5, 0.8], 'MarkerEdgeColor', 'k')
hold on
plot([min(phi_exp)-1, max(phi_exp)+1], [min(phi_exp)-1, max(phi_exp)+1], ...
'r--', 'LineWidth', 2)
xlabel('实验值 (\phi_{exp}, %)')
ylabel('预测值 (\phi_{pred}, %)')
title('预测值 vs 实验值', 'FontSize', 14)
text(min(phi_exp)+2, max(phi_pred)-1, sprintf('R² = %.4f', R2), 'FontSize', 14)
text(min(phi_exp)+2, max(phi_pred)-2, sprintf('RMSE = %.4f%%', RMSE), 'FontSize', 14)
grid on
axis equal
set(gca, 'FontSize', 12, 'Box', 'on')
% 残差分析
nexttile
residuals = phi_pred - phi_exp;
scatter(phi_pred, residuals, 100, 'filled', ...
'MarkerFaceColor', [0.8, 0.4, 0.1], 'MarkerEdgeColor', 'k')
hold on
plot([min(phi_pred)-1, max(phi_pred)+1], [0, 0], 'k-', 'LineWidth', 1.5)
plot([min(phi_pred)-1, max(phi_pred)+1], [mean(residuals)+std(residuals), ...
mean(residuals)+std(residuals)], 'r--')
plot([min(phi_pred)-1, max(phi_pred)+1], [mean(residuals)-std(residuals), ...
mean(residuals)-std(residuals)], 'r--')
xlabel('预测值 (\phi_{pred}, %)')
ylabel('残差 (\phi_{pred} - \phi_{exp}, %)')
title('残差分析', 'FontSize', 14)
text(min(phi_pred)+2, max(residuals)-0.5, ...
sprintf('平均残差 = %.4f%%', mean(residuals)), 'FontSize', 12)
text(min(phi_pred)+2, max(residuals)-1, ...
sprintf('标准差 = %.4f%%', std(residuals)), 'FontSize', 12)
grid on
set(gca, 'FontSize', 12, 'Box', 'on')
%% 最优制备条件求解
% 定义优化函数
obj_func = @(x) -model_func(params_opt, [x(1), x(2), x(3)]); % 负值用于最大化
% 约束条件
A = [];
b = [];
Aeq = [];
beq = [];
lb_opt = [30, 50, 5]; % T_min, H_min, SC_min
ub_opt = [50, 90, 11]; % T_max, H_max, SC_max
% 多起点优化
options_opt = optimoptions('fmincon', 'Display', 'off', ...
'Algorithm', 'sqp', 'MaxFunctionEvaluations', 1000);
% 存储多个起点的结果
solutions = zeros(10, 4); % [T, H, SC, phi]
for i = 1:10
% 随机起点
x0 = [30+20*rand, 50+40*rand, 5+6*rand];
% 优化
[x_opt, fval] = fmincon(obj_func, x0, A, b, Aeq, beq, lb_opt, ub_opt, [], options_opt);
% 存储结果
solutions(i, :) = [x_opt, -fval];
end
% 选择最佳解
[~, idx] = max(solutions(:,4));
optimal_T = solutions(idx,1);
optimal_H = solutions(idx,2);
optimal_SC = solutions(idx,3);
max_phi = solutions(idx,4);
fprintf('\n===== 最优制备条件 =====\n');
fprintf('最优温度: %.2f°C\n', optimal_T);
fprintf('最优湿度: %.2f%%\n', optimal_H);
fprintf('最优固含量: %.2f%%\n', optimal_SC);
fprintf('预测最大孔面积占比: %.2f%%\n', max_phi);
%% 与原模型比较
% 原模型最优解 (假设值)
orig_optimal = [50, 50, 6]; % T, H, SC
orig_phi = 77.5;
% 计算新模型在原最优点的表现
new_phi_origpoint = model_func(params_opt, orig_optimal);
fprintf('\n===== 模型比较 =====\n');
fprintf('原模型最优解: T=%.1f°C, H=%.1f%%, SC=%.1f%%, φ=%.2f%%\n', ...
orig_optimal(1), orig_optimal(2), orig_optimal(3), orig_phi);
fprintf('新模型在原最优解: φ=%.2f%%\n', new_phi_origpoint);
fprintf('新模型最优解: T=%.1f°C, H=%.1f%%, SC=%.1f%%, φ=%.2f%%\n', ...
optimal_T, optimal_H, optimal_SC, max_phi);
fprintf('性能提升: Δφ = %.2f%%\n', max_phi - orig_phi);
% 可视化比较
figure('Position', [100, 100, 1000, 500])
tiledlayout(1, 2, 'TileSpacing', 'compact')
% 最优解比较
nexttile
bar([1, 2], [orig_phi, max_phi], 0.6, 'FaceColor', [0.5, 0.7, 0.9])
set(gca, 'XTickLabel', {'原模型最优', '新模型最优'})
ylabel('孔面积占比 (\phi, %)')
title('最优解比较', 'FontSize', 14)
text(1, orig_phi+1, sprintf('T=%d,H=%d,SC=%d', ...
orig_optimal(1), orig_optimal(2), orig_optimal(3)), ...
'HorizontalAlignment', 'center', 'FontSize', 10)
text(2, max_phi+1, sprintf('T=%.1f,H=%.1f,SC=%.1f', ...
optimal_T, optimal_H, optimal_SC), ...
'HorizontalAlignment', 'center', 'FontSize', 10)
ylim([60, 85])
grid on
set(gca, 'FontSize', 12)
% 参数空间探索
nexttile
scatter3(solutions(:,1), solutions(:,2), solutions(:,3), 100, ...
solutions(:,4), 'filled')
xlabel('温度 (°C)')
ylabel('湿度 (%)')
zlabel('固含量 (%)')
title('参数空间探索与最优解', 'FontSize', 14)
hold on
% 标记原最优解
scatter3(orig_optimal(1), orig_optimal(2), orig_optimal(3), 150, ...
'ro', 'LineWidth', 2)
% 标记新最优解
scatter3(optimal_T, optimal_H, optimal_SC, 200, ...
'kh', 'LineWidth', 2, 'MarkerFaceColor', 'y')
colorbar
colormap(jet)
view(45, 30)
grid on
set(gca, 'FontSize', 12)
legend('探索点', '原最优解', '新最优解', 'Location', 'best')
%% 3D响应曲面图 (固定固含量)
figure('Position', [100, 100, 1200, 400])
tiledlayout(1, 3, 'TileSpacing', 'compact')
sc_values = [6, 8, 10];
titles = {'固含量 6%', '固含量 8%', '固含量 10%'};
for i = 1:3
nexttile
% 创建网格
[T_grid, H_grid] = meshgrid(linspace(30,50,20), linspace(50,90,20));
SC_grid = sc_values(i) * ones(size(T_grid));
% 预测值
phi_grid = model_func(params_opt, [T_grid(:), H_grid(:), SC_grid(:)]);
phi_grid = reshape(phi_grid, size(T_grid));
% 曲面图
surf(T_grid, H_grid, phi_grid, 'EdgeColor', 'none')
hold on
% 添加等高线
[C, h] = contour3(T_grid, H_grid, phi_grid, 10, 'k-');
clabel(C, h, 'FontSize', 10)
% 设置视图
view(3)
xlabel('温度 (°C)')
ylabel('湿度 (%)')
zlabel('\phi (%)')
title(titles{i}, 'FontSize', 12)
colorbar
colormap(jet)
grid on
set(gca, 'FontSize', 10)
end
sgtitle('不同固含量下的响应曲面', 'FontSize', 16, 'FontWeight', 'bold')
%% 温度湿度弱化效应可视化
figure('Position', [100, 100, 1000, 800])
tiledlayout(2, 2, 'TileSpacing', 'compact')
% 温度影响 (不同固含量)
nexttile
hold on
T_range = linspace(30, 50, 100);
for sc = [6, 8, 10]
phi_temp = model_func(params_opt, [T_range', 70*ones(100,1), sc*ones(100,1)]);
plot(T_range, phi_temp, 'LineWidth', 2, ...
'DisplayName', sprintf('SC=%d%%', sc))
end
xlabel('温度 (°C)')
ylabel('\phi (%)')
title('温度影响 (固定湿度70%)', 'FontSize', 14)
legend('Location', 'best')
grid on
set(gca, 'FontSize', 12)
% 湿度影响 (不同固含量)
nexttile
hold on
H_range = linspace(50, 90, 100);
for sc = [6, 8, 10]
phi_humid = model_func(params_opt, [40*ones(100,1), H_range', sc*ones(100,1)]);
plot(H_range, phi_humid, 'LineWidth', 2, ...
'DisplayName', sprintf('SC=%d%%', sc))
end
xlabel('湿度 (%)')
ylabel('\phi (%)')
title('湿度影响 (固定温度40°C)', 'FontSize', 14)
legend('Location', 'best')
grid on
set(gca, 'FontSize', 12)
% 固含量影响 (不同温度)
nexttile
hold on
SC_range = linspace(5, 11, 100);
for t = [30, 40, 50]
phi_sc = model_func(params_opt, [t*ones(100,1), 70*ones(100,1), SC_range']);
plot(SC_range, phi_sc, 'LineWidth', 2, ...
'DisplayName', sprintf('T=%d°C', t))
end
xlabel('固含量 (%)')
ylabel('\phi (%)')
title('固含量影响 (固定湿度70%)', 'FontSize', 14)
legend('Location', 'best')
grid on
set(gca, 'FontSize', 12)
% 参数影响权重
nexttile
weights = [abs(w_T)*10, abs(w_H)*5, 1]; % 归一化因子
labels = {'温度', '湿度', '固含量'};
pie(weights, labels)
title('参数影响权重比较', 'FontSize', 14)
应该与这个代码相结合吧,结果应该更加符合物理逻辑一些,不要线性