【研赛C题成品论文】24华为杯数学建模研赛C题成品论文+可运行代码丨免费分享

2024华为杯研究生数学建模竞赛C题成品论文已出！

C题 数据驱动下磁性元件的磁芯损耗建模

一、问题一模型建立与求解——基于SVM的分类模型

1.1 SVM算法简介

支持向量机（Support Vector Machine，SVM）是一种经典的监督学习算法，用于解决二分类和多分类问题。其核心思想是通过在特征空间中找到一个最优的超平面来进行分类，并且间隔最大。

SVM能够执行线性或非线性分类、回归，甚至是异常值检测任务。它是机器学习领域最受欢迎的模型之一。SVM特别适用于中小型复杂数据集的分类。

超平面最大间隔介绍：下左图显示了三种可能的线性分类器的决策边界；右图中的实线代表SVM分类器的决策边界，不仅分离了两个类别，且尽可能远离最近的训练实例。

虚线所代表的模型表现非常糟糕，甚至都无法正确实现分类。其余两个模型在这个训练集上表现堪称完美，但是它们的决策边界与实例过于接近，导致在面对新实例时，表现可能不会太好。

SVM是通过优化一个凸二次规划问题来求解最佳的超平面，其中包括最小化模型的复杂度（即最小化权重的平方和），同时限制训练样本的误分类情况。这个优化问题可以使用拉格朗日乘子法来求解。对于非线性可分的情况，SVM可以通过核函数（Kernel Function）将输入特征映射到高维空间，使得原本线性不可分的数据在高维空间中变得线性可分。常用的核函数包括线性核、多项式核、高斯核等。

以上决策函数就称为线性可分支持向量机。

Φ(x)  是某个确定的特征空间转换函数, 它的作用是将x映射到更高的维度, 它有一个以后我们经常会见到的专有称号"核函数"。

以上就是线性可分支持向量机的模型表达式。我们要去求出这样一个模型，或者说这样一个超平面y(x) ，它能够最优地分离两个集合。

其实也就是我们要去求一组参数(w,b) ，使其构建的超平面函数能够最优地分离两个集合。如下就是一个最优超平面：

再比如下图的阴影部分是一个“过渡带”，“过渡带”的边界是集合中离超平面最近的样本点落在的地方：

1.2 特征提取与分析

数据分析：从给定数据集中，磁通密度的波形数据（第5到1029列）包含了不同励磁波形（正弦波、三角波、梯形波）的特征。需要提取反映这些波形特征的变量。

特征选择：

·波形的基本统计特征：均值、标准差、最大值、最小值、峰值峰值差等。

·波形的形状特征：

·峰度和偏度：描述波形的对称性和尖锐程度。

 （1）均值（Mean, μ ）：

（2）标准差（Standard Deviation, σ ）

（3）最大值（Maximum, Bmax ）与最小值（Minimum, Bmin ）

（4）峰值峰值差（Peak-to-Peak）

（5）偏度（Skewness）

（6）上升沿斜率（Slope of Rising Edges）

上升沿斜率描述波形从低谷上升到峰值时的变化速率，即波形上升过程中每个采样点的斜率，取斜率的最大值。

（7）下降沿斜率（Slope of Falling Edges）

下降沿斜率描述波形从峰值下降到低谷时的变化速率，即波形下降过程中每个采样点的斜率，取斜率的最小值。

（8）零交点数（Number of Zero Crossings）

零交点数表示波形与零点的交点个数，即波形在一个周期内穿过零点的次数。该特征能够反映波形的周期性和对称性。

sign(Bi) 表示第 i 个采样点的符号（正数为1，负数为-1）。

如果相邻两个采样点的符号相反，则波形在这两个采样点之间穿过零点。

1.3 模型建立与求解

基于以上特征值训练SVM模型，模型准确率接近100%。

特征分布直方图：展示各个特征的分布情况，有助于理解特征的集中趋势及离散情况。

特征相关性热力图：可视化特征之间的相关性，帮助理解特征的关系。

       最终求得预测结果准确率达到100%：

具体结果如下表所示：

序号	附件二（80个样品）励磁波形分类结果
0	2
1	2
2	2
3	2
4	2
5	3
6	3
7	3
8	3
9	3
10	1
11	1
12	1
13	1
14	1
15	2
16	2
17	2
18	2
……	……

统计出附件二中三种波形的各自数量，如下图所示：

可以看出，其中三角波占比最多，大于其他两类的总和。

1..4参考代码

import pandas as pd
import numpy as np
from scipy.stats import skew, kurtosis
from scipy.fft import fft  # 导入傅里叶变换函数

# 读取磁通密度数据（假设已经存储在DataFrame中）
magnetic_flux_density = pd.read_excel('附件一（训练集）.xlsx', sheet_name='材料1').iloc[:, 5:1029]  # 假设数据位于第5列到第1029列
label = pd.read_excel('附件一（训练集）.xlsx', sheet_name='材料1').iloc[:, 3]  # label在第4列
# 其他几个sheet的数据读取方式类似，只需修改sheet_name即可，放到一个循环中处理更好
for i in range(2, 5):
    magnetic_flux_density = pd.concat([magnetic_flux_density, pd.read_excel('附件一（训练集）.xlsx', sheet_name=f'材料{i}').iloc[:, 5:1029]], axis=0)
    # label在第4列, 读取标签数据，放到一个列表中
    label = pd.concat([label, pd.read_excel('附件一（训练集）.xlsx', sheet_name=f'材料{i}').iloc[:, 3]], axis=0)

    print(f"第{i}个材料的数据读取完毕")

# 检查是否存在缺失值
missing_count = magnetic_flux_density.isnull().sum().sum()
print(f"缺失值总数：{missing_count}")

# 填补缺失值（以列的均值进行填补，也可以用其他值替换，如0）
magnetic_flux_density_filled = magnetic_flux_density.fillna(magnetic_flux_density.mean())

# 计算基本统计特征
mean_val = magnetic_flux_density_filled.mean(axis=1)  # 每行均值
std_val = magnetic_flux_density_filled.std(axis=1)    # 每行标准差
max_val = magnetic_flux_density_filled.max(axis=1)    # 每行最大值
min_val = magnetic_flux_density_filled.min(axis=1)    # 每行最小值
peak_to_peak = max_val - min_val                      # 每行峰峰值
skewness = magnetic_flux_density_filled.apply(lambda row: skew(row), axis=1)  # 每行偏度
kurt = magnetic_flux_density_filled.apply(lambda row: kurtosis(row), axis=1)  # 每行峰度

# 构建特征DataFrame
features = pd.DataFrame({
    'Mean': mean_val,
    'Std': std_val,
    'Max': max_val,
    'Min': min_val,
    'Peak-to-Peak': peak_to_peak,
    'Skewness': skewness,
    'Kurtosis': kurt
})

label

# 输出magnetic_flux_density、feature和label的shape
print(f"磁通密度数据的shape：{magnetic_flux_density_filled.shape}")
print(f"特征数据的shape：{features.shape}")
print(f"label的shape：{len(label)}")

# 计算上升沿和下降沿的斜率
slope_rising = magnetic_flux_density.diff(axis=1).apply(lambda row: np.max(row), axis=1)
slope_falling = magnetic_flux_density.diff(axis=1).apply(lambda row: np.min(row), axis=1)

# 计算零交点数
zero_crossings = magnetic_flux_density.apply(lambda row: np.sum(np.diff(np.sign(row)) != 0), axis=1)

# 添加到特征表中
features['Slope_Rising'] = slope_rising
features['Slope_Falling'] = slope_falling
features['Zero_Crossings'] = zero_crossings


# 可视化和验证特征
import matplotlib.pyplot as plt
import seaborn as sns

# # 可视化不同类别的特征分布（假设已经有标签列 'Label'）
# sns.pairplot(features, hue='励磁波形')  

# 分类模型（示例：SVM）
from sklearn.model_selection import train_test_split
from sklearn.svm import SVC
from sklearn.metrics import accuracy_score

# 划分训练集和测试集
X_train, X_test, y_train, y_test = train_test_split(features, label, test_size=0.3, ran-dom_state=42)

# 训练SVM模型
clf = SVC()
clf.fit(X_train, y_train)

# 预测
y_pred = clf.predict(X_test)

# 计算准确率
accuracy = accuracy_score(y_test, y_pred)
print(f"准确率：{accuracy}")

# 基于附件二中的数据，假设数据位于第5列到第1029列，分类出波形类型，输出在附件四的第二列上
magnetic_flux_density_test = pd.read_excel('附件二（测试集）.xlsx').iloc[:, 5:1029]
# 填补缺失值（以列的均值进行填补，也可以用其他值替换，如0）
magnetic_flux_density_test_filled = magnet-ic_flux_density_test.fillna(magnetic_flux_density_test.mean())

# 计算基本统计特征
mean_val_test = magnetic_flux_density_test_filled.mean(axis=1)  # 每行均值
std_val_test = magnetic_flux_density_test_filled.std(axis=1)    # 每行标准差
max_val_test = magnetic_flux_density_test_filled.max(axis=1)    # 每行最大值
min_val_test = magnetic_flux_density_test_filled.min(axis=1)    # 每行最小值
peak_to_peak_test = max_val_test - min_val_test                      # 每行峰峰值
skewness_test = magnetic_flux_density_test_filled.apply(lambda row: skew(row), axis=1)  # 每行偏度
kurt_test = magnetic_flux_density_test_filled.apply(lambda row: kurtosis(row), axis=1)  # 每行峰度

# 构建特征DataFrame
features_test = pd.DataFrame({
    'Mean': mean_val_test,
    'Std': std_val_test,
    'Max': max_val_test,
    'Min': min_val_test,
    'Peak-to-Peak': peak_to_peak_test,
    'Skewness': skewness_test,
    'Kurtosis': kurt_test
})

# 计算上升沿和下降沿的斜率
slope_rising_test = magnetic_flux_density_test.diff(axis=1).apply(lambda row: np.max(row), axis=1)
slope_falling_test = magnetic_flux_density_test.diff(axis=1).apply(lambda row: np.min(row), axis=1)

# 计算零交点数
zero_crossings_test = magnetic_flux_density_test.apply(lambda row: np.sum(np.diff(np.sign(row)) != 0), axis=1)

# 添加到特征表中
features_test['Slope_Rising'] = slope_rising_test
features_test['Slope_Falling'] = slope_falling_test
features_test['Zero_Crossings'] = zero_crossings_test

# 预测
y_pred_test = clf.predict(features_test)

# 统计出我们预测结果中三种波形的各自数量，呈现出来，字体都是中文
import matplotlib.pyplot as plt

# 统计波形数量
waveform_counts = pd.Series(y_pred_test).value_counts()

# 设置中文显示
plt.rcParams['font.sans-serif'] = ['SimHei']

# 绘制柱状图
plt.bar(waveform_counts.index, waveform_counts.values)

# 把数量显示在柱状图上
for x, y in zip(waveform_counts.index, waveform_counts.values):
    plt.text(x, y, f'{y}', ha='center', va='bottom')

# 添加标题和标签
plt.title('预测结果中三种波形的数量统计')
plt.xlabel('波形类型')
plt.ylabel('数量')

# 显示图形
plt.show()

# 特别把预测结果样本序号为：1、5、15、25、35、45、55、65、75、80的分类结果，以表格形式呈现在论文正文中。
# 由于这里的数据是从0开始的，所以需要减1

# 选取样本序号
sample_indices = [1, 5, 15, 25, 35, 45, 55, 65, 75, 80]
# 选取样本
sample_indices = np.array(sample_indices)
selected_samples = pd.DataFrame(y_pred_test).iloc[sample_indices - 1]
# 索引用sample_indices
selected_samples.index = sample_indices
# 修改列名
selected_samples.columns = ['波形类型']
# 输出
print(selected_samples)

# 保存到Excel
selected_samples.to_excel('特殊预测结果样本序号.xlsx', index=False)


# 把结果替代成[1,2,3]，1表示正弦波，2表示三角波，3表示梯形波，比如：附件二中第1个样品分类结果是三角波，在第2列样本序号为1对应行就填数字2；
y_pred_test = pd.Series(y_pred_test).map({'正弦波': 1, '三角波': 2, '梯形波': 3})

# 输出预测结果在附件四（Excel表）.xlsx，只从第二行第二列开始写入，其余别覆盖
y_pred_test.to_excel('附件四（Excel表）.xlsx', index=True, header=['附件二（80个样品）励磁波形分类结果'])


# 绘制饼图
plt.figure(figsize=(6, 6))
plt.pie(waveform_counts.values, labels=waveform_counts.index, autopct='%1.1f%%', star-tangle=90)
plt.title('预测结果中三种波形的比例')
plt.show()
# 绘制特征的分布直方图
features_test.hist(figsize=(12, 10), bins=30)
plt.suptitle('特征分布直方图', size=16)
plt.show()


# 计算特征之间的相关性矩阵
correlation_matrix = features_test.corr()

# 绘制热力图
plt.figure(figsize=(10, 8))
sns.heatmap(correlation_matrix, annot=True, cmap='coolwarm', linewidths=0.5)
plt.title('特征相关性热力图')
plt.show()

二、 问题二模型建立与求解

2.1 问题二求解思路

磁芯损耗是磁性材料在交变磁场中的一个重要性能指标，对电力电子设备的效率和热管理有重要影响。常见的磁芯损耗模型包括原始斯坦麦茨方程（Original Steinmetz Equation, OSE），但因为其虽作为经典模型被广泛应用，却显著受限于其特定的适用条件，如：该方程主要针对正弦波形设计；对于不同种类的磁芯材料及工作温度的变化，SE方程会造成较大的误差，这在实际工程应用中带来了诸多不便与复杂性。所以需要考虑修正斯坦麦茨方程（Modified Steinmetz Equation, MSE），即在原斯坦麦茨方程基础上，增加温度这个因素，以适应不同温度变化，使磁芯损耗预测效果更好。

2.2问题二模型建立

（1）原始斯坦麦茨方程

原始斯坦麦茨方程是基于实验数据的经验公式，用于描述磁性材料在交变磁场下的损耗。其数学形式为：

其中：

- P  是单位体积的磁芯损耗（单位：W/m³）。

-  f 是 外加磁场的频率（单位：Hz）。

-Bmax 是磁通密度的峰值（单位：T）。

- k、α、β 是与材料特性相关的常数。

 k  是比例常数，反映了材料的固有损耗特性。 α 反映了频率对损耗的影响，通常在1到2之间。 β 反映了磁通密度峰值对损耗的影响，通常大于2。

尽管OSE在一定范围内能够有效地描述磁芯损耗，但它忽略了温度对磁芯损耗的影响。温度对磁性材料的磁滞损耗、涡流损耗以及剩磁损耗都有显著影响，因此需要对OSE进行修正。

（2）修正斯坦麦茨方程（MSE）

为了更精确地描述磁芯损耗，特别是在考虑温度变化的情况下，引入了修正斯坦麦茨方程。其数学形式为：

其中：

-  T  是温度（单位：°C）。

- γ、δ 是与温度相关的修正系数。

 γ  表示温度一阶对磁芯损耗的修正量。δ  表示温度二阶对磁芯损耗的修正量。这些温度修正项使得模型能够更好地适应实际工作环境中温度变化对磁芯损耗的影响。

（3）修正模型的理论基础

1.磁滞损耗

磁滞损耗是由于磁化过程中磁畴壁运动引起的内耗，随着温度的升高，磁畴壁的移动变得更加容易，从而影响磁滞回线的面积，进而影响总损耗。

2.涡流损耗

涡流损耗由交变磁场在磁性材料内部感应的电流引起。温度升高会改变材料的电导率，影响涡流的大小和分布，进而影响涡流损耗。

3.剩磁损耗

剩磁损耗是材料在磁化后保持残余磁化强度所引起的损耗。温度变化会影响剩磁强度和矫顽力，从而影响这部分损耗。

2.3 问题二的求解与分析

为了确定修正斯坦麦茨方程中的参数 k 、 α、 β、 γ、 δ， 需要使用实验数据进行非线性拟合。

1. 数据准备

实验数据包括在不同温度、频率和磁通密度峰值下测量的实际磁芯损耗 Pactual 。这些数据可以表示为：

2. 拟合过程

通过最小化模型预测值与实际值之间的误差，可以获得模型参数。具体方法是最小化以下目标函数：

3. 模型评估

为了评估模型的性能，可以使用以下指标：

- 均方误差（Mean Squared Error, MSE）：

- 平均绝对误差（Mean Absolute Error, MAE）：

- 拟合优度（R²）：

其中，Ppred,i  为模型预测值，Pactual 为实际值的平均值。

4. 结果与分析

在对损耗模型进行了参数拟合之后，我们得到了一组优化后的参数以及相关的评价指标，以下是详细的结果展示和分析：

（1）拟合参数和评价指标

拟合参数:

 拟合优度:

 均方误差:

MSE=134806625.10564047

 平均绝对误差:

MAE≈0.2219

（2）原斯坦麦茨方程和修正后的斯坦麦茨方程对比

我们对原斯坦麦茨方程和修正后的斯坦麦茨方程进行了比较。原斯坦麦茨方程的初始参数估计为：

原斯坦麦茨方程的评价指标

 拟合优度:

 平均绝对误差:

 修正斯坦麦茨方程的评价指标

如上文所述，修正后的斯坦麦茨方程拟合参数和评价指标为：

 拟合优度:

 均方误差:

平均绝对误差:

（3）结果分析

1. 拟合参数的含义:

     k 、 α  和  β 是斯坦麦茨方程中的核心参数，分别代表了损耗与频率、最大磁通密度的关系。

     γ  和  δ 是修正项中的参数，考虑了温度对损耗的影响。

2. 拟合优度 :

    修正斯坦麦茨方程的R2 值为 0.9955，远高于原斯坦麦茨方程的 0.0577，表明修正模型能更好地解释数据的变异性。

3. 误差分析:

    修正方程的均方误差MSE 和平均绝对误差MAE 都显著低于原方程，表明修正模型在预测性能上有显著提升。

MAE=0.2219

可视化展示：

为了直观展示模型的预测性能，上图呈现了修正Steinmetz方程与原始Steinmetz方程的预测结果与实际磁芯损耗之间的对比。横轴表示实际测量的磁芯损耗，纵轴表示模型预测的磁芯损耗。蓝色点代表修正模型的预测结果，红色点代表原始Steinmetz模型的预测结果，黑色虚线则表示理想情况下的y=x参考线，即模型预测值与实际值完全一致时的位置。

从图中可以看出，修正模型的预测结果与实际损耗值更为接近，点云分布更加接近y=x参考线。而原始Steinmetz模型的预测结果则存在更大的偏差，表明其在某些数据点上可能存在系统性误差。这表明修正后的Steinmetz方程通过引入温度等影响因素，提升了模型的拟合能力。

图2展示了原始Steinmetz方程与修正Steinmetz方程在预测磁芯损耗时的相对误差分布情况。相对误差通过预测值与实际值的比值计算而得，横轴表示误差大小，纵轴表示误差的出现频率。蓝色柱状图对应修正Steinmetz模型的误差分布，红色柱状图则对应原始Steinmetz模型的误差分布。

从误差分布的对比中可以看出，修正模型的误差集中分布在较低的范围内，峰值更高且分布更窄，说明大多数预测点的误差较小。而原始Steinmetz模型的误差分布较为分散，误差相对较大。这进一步验证了修正模型在处理不同温度和磁通密度条件下的预测性能优于原始模型。

2.4 结论

通过参数拟合和模型修正，我们显著提高了损耗模型的预测精度。修正后的斯坦麦茨方程不仅拟合优度更高，而且在误差指标上也优于原始模型。这表明引入温度影响的修正模型更能准确地反映实际情况，对工程应用具有更高的实用价值。

2.5 参考代码

import pandas as pd
import numpy as np
from scipy.optimize import curve_fit

# 读取材料1中的正弦波形数据
data = pd.read_excel('附件一（训练集）.xlsx', sheet_name='材料1')

# 筛选正弦波形的数据
sine_data = data[data['励磁波形'] == '正弦波']  # 假设波形类型列名为'波形类型'
temperature = sine_data['温度，oC']  # 温度列
frequency = sine_data['频率，Hz']  # 频率列
P_actual = sine_data['磁芯损耗，w/m3']  # 实际磁芯损耗列
# 磁通密度峰值列，从第五列开始到底为止，每行的最大值
B_max = sine_data.iloc[:, 4:].max(axis=1)

# 温度、频率、磁通密度数据预处理
X = np.vstack([frequency, B_max, temperature]).T  # 特征矩阵
y = P_actual  # 实际损耗
X

# 定义修正后的损耗模型
def loss_model(X, k, alpha, beta, gamma, delta):
    f, B_max, T = X.T
    return k * (f ** alpha) * (B_max ** beta) * (1 + gamma * T + delta * (T ** 2))

# 初始参数估计
initial_params = [1e-4, 1.5, 2.0, 0.001, 0.0001]

# 拟合修正方程参数
params_opt, params_cov = curve_fit(loss_model, X, y, p0=initial_params)

# 输出拟合的参数
k_opt, alpha_opt, beta_opt, gamma_opt, delta_opt = params_opt
print(f"拟合参数：k = {k_opt}, alpha = {alpha_opt}, beta = {beta_opt}, gamma = {gam-ma_opt}, delta = {delta_opt}")

# 拟合优度、MSE等评价指标
P_pred_corrected = loss_model(X, *params_opt)
r2 = 1 - np.sum((y - P_pred_corrected) ** 2) / np.sum((y - np.mean(y)) ** 2)
print(f"拟合优度R^2 = {r2}")
mse = np.mean((y - P_pred_corrected) ** 2)
print(f"均方误差MSE = {mse}")
mae = np.mean(np.abs(y - P_pred_corrected))

# 定义原始Steinmetz方程
def original_steinmetz(X, k, alpha, beta):
    f, B_max = X.T
    return k * (f ** alpha) * (B_max ** beta)
    
# 拟合原始Steinmetz方程

# 初始参数估计 (这里是一个合理的估计，可以根据具体问题调整)
initial_params = [1e-4, 1.5, 2.0]

# 拟合原始Steinmetz方程参数
params_opt, params_cov = curve_fit(original_steinmetz, X[:, :2], y, p0=initial_params)

# 输出拟合的参数
k0_opt, alpha0_opt, beta0_opt = params_opt
print(f'k0={k0_opt}, alpha0={alpha0_opt}, beta0={beta0_opt}')

# 使用拟合参数预测损耗
P_pred_original = original_steinmetz(X[:, :2], *params_opt)

# 定义原始斯坦麦茨方程
def original_steinmetz(f, B_max, k, alpha, beta):
    return k * (f ** alpha) * (B_max ** beta)

# 计算原斯坦麦茨方程的预测损耗
P_pred_original = original_steinmetz(frequency, B_max, k_opt, alpha_opt, beta_opt)

# 计算修正后的损耗预测
P_pred_corrected = loss_model(X, k_opt, alpha_opt, beta_opt, gamma_opt, delta_opt)

# 计算两种模型的误差
error_original = np.abs(P_actual - P_pred_original) / P_actual
error_corrected = np.abs(P_actual - P_pred_corrected) / P_actual

# 计算平均误差
mean_error_original = np.mean(error_original)
mean_error_corrected = np.mean(error_corrected)

print(f"原斯坦麦茨方程的平均误差：{mean_error_original:.4f}")
print(f"修正斯坦麦茨方程的平均误差：{mean_error_corrected:.4f}")

# 比较拟合优度等评价指标

# 计算原斯坦麦茨方程的拟合优度
r2_original = 1 - np.sum((P_actual - P_pred_original) ** 2) / np.sum((P_actual - np.mean(P_actual)) ** 2)

# 计算修正后的损耗预测的拟合优度
r2_corrected = 1 - np.sum((P_actual - P_pred_corrected) ** 2) / np.sum((P_actual - np.mean(P_actual)) ** 2)

print(f"原斯坦麦茨方程的拟合优度：{r2_original:.4f}")
print(f"修正斯坦麦茨方程的拟合优度：{r2_corrected:.4f}")

import matplotlib.pyplot as plt
# 计算原斯坦麦茨方程的预测损耗
P_pred_original = original_steinmetz(frequency, B_max, k_opt, alpha_opt, beta_opt)

# 计算修正后的损耗预测
P_pred_corrected = loss_model(X, k_opt, alpha_opt, beta_opt, gamma_opt, delta_opt)

# 计算两种模型的误差
error_original = np.abs(P_actual - P_pred_original) / P_actual
error_corrected = np.abs(P_actual - P_pred_corrected) / P_actual

# 计算平均误差
mean_error_original = np.mean(error_original)
mean_error_corrected = np.mean(error_corrected)

# 实际损耗与修正模型预测损耗的对比
plt.figure(figsize=(12, 6))
plt.scatter(y, P_pred_corrected, label='修正模型预测', color='blue', alpha=0.6)
plt.scatter(y, P_pred_original, label='原始模型预测', color='red', alpha=0.6)
plt.plot([y.min(), y.max()], [y.min(), y.max()], 'k--', lw=2)  # y=x参考线
plt.xlabel('实际损耗 (W/m³)')
plt.ylabel('预测损耗 (W/m³)')
plt.title('实际损耗与模型预测的对比')
plt.legend()
plt.grid()
plt.show()

# 误差分布
error_corrected = np.abs(P_actual - P_pred_corrected) / P_actual
error_original = np.abs(P_actual - P_pred_original) / P_actual

plt.figure(figsize=(12, 6))
plt.hist(error_corrected, bins=20, alpha=0.5, label='修正模型误差', color='blue')
plt.hist(error_original, bins=20, alpha=0.5, label='原始模型误差', color='red')
plt.xlabel('相对误差')
plt.ylabel('频率')
plt.title('模型预测误差分布')
plt.legend()
plt.grid()
plt.show()