<think>我们遇到了一个新的错误,发生在反归一化(inverse_transform)步骤:
ValueError: non-broadcastable output operand with shape (1095,1) doesn't match the broadcast shape (1095,7)
错误原因分析:
在`ensemble_predict`方法中,我们使用`self.scalers[fuel_type]`对预测值进行反归一化。这个scaler是在`prepare_data`中用于对整个特征矩阵(包含7个特征)进行归一化的。当我们调用`inverse_transform`时,我们需要传入一个形状为`(n_samples, n_features)`的数组,其中`n_features`应该与归一化时使用的特征数相同(这里是7个特征)。但是,我们只传入了预测值(1个特征)的数组,形状为`(1095, 1)`,而scaler期望的是7个特征。
解决方案:
我们需要只反归一化目标特征(即价格特征),而不是整个特征向量。但是,由于归一化时是针对所有特征一起进行的,所以我们需要构建一个包含所有特征值的矩阵,但只使用目标特征列的值,其他特征列可以设为0(因为归一化是线性变换,我们只需要目标特征列的反归一化结果)。然而,更简单的方法是我们应该只使用目标特征对应的归一化参数(均值和标准差)来反归一化。但MinMaxScaler并没有提供按列反归一化的方法。
因此,我们可以修改归一化的方式:只对目标列进行归一化,而不是所有特征。但这样会破坏多变量输入的结构,因为其他特征也需要归一化。所以,我们需要为每个特征单独归一化?或者使用一个scaler来归一化多个特征?但这样在反归一化时,我们只能得到整个矩阵的反归一化,而我们只关心目标列。
另一种方法:在归一化时,我们使用同一个scaler对多个特征进行归一化,但在反归一化时,我们只取目标列(即第一列)的反归一化值。但是,由于归一化是线性变换,每个特征都是独立缩放的,所以我们可以这样操作:
1. 创建一个与原始归一化数据相同形状的矩阵,其中除了第一列(目标列)是我们需要反归一化的值,其他列设为0(或任意值,因为我们不会使用它们)。
2. 对这个矩阵进行反归一化。
3. 然后只取第一列作为反归一化的目标值。
但是,这样会引入错误,因为其他列的值被设为0,而0在归一化后的空间中可能代表一个实际的值,反归一化后这些0会变成原始特征空间中的值,从而影响目标列的反归一化?实际上,MinMaxScaler的公式是:
X_std = (X - X.min(axis=0)) / (X.max(axis=0) - X.min(axis=0))
X_scaled = X_std * (max - min) + min
反归一化公式:
X_original = X_scaled * (max - min) + min
注意:它是按列(每个特征独立)进行归一化和反归一化的。所以,当我们反归一化时,每一列都是独立进行的。因此,我们可以构建一个矩阵,其中第一列是归一化后的预测值,其他列为0(或任意值,因为反归一化时只使用该列的参数),然后反归一化后,只取第一列。
但是,由于其他列在反归一化时会被还原为原始范围,而我们并不关心它们,所以我们可以将其他列设为0,然后反归一化,这样其他列会变成原始数据中该列的最小值(因为0在归一化后的特征中对应最小值)。但是,这不会影响目标列的反归一化,因为每个特征独立。
因此,我们可以这样做:
假设我们有一个归一化后的预测值数组`lstm_preds`,形状为`(n, 1)`,其中n是预测天数(1095)。我们需要将其扩展为`(n, 7)`的矩阵,其中第一列是归一化后的预测值,其他列设为0(或任意值,但为了简单设为0)。然后对这个矩阵进行反归一化,再取第一列。
但是,在`prepare_data`中,我们使用整个特征矩阵(7列)进行归一化,所以scaler的`min_`和`scale_`(即范围)都是长度为7的向量。因此,在反归一化时,我们需要传入7列的数据。
修改`ensemble_predict`方法中的反归一化步骤:
原始代码:
```python
lstm_preds = self.scalers[fuel_type].inverse_transform(
np.array(lstm_preds).reshape(-1, 1)
).flatten()
```
修改为:
```python
# 创建与归一化特征数相同的矩阵,将预测值放在第一列,其他列设为0
n_features = self.scalers[fuel_type].n_features_in_
# 创建一个全0矩阵,形状为 (len(lstm_preds), n_features)
preds_matrix = np.zeros((len(lstm_preds), n_features))
# 将lstm_preds赋值到第一列
preds_matrix[:, 0] = np.array(lstm_preds).flatten()
# 反归一化整个矩阵
preds_matrix = self.scalers[fuel_type].inverse_transform(preds_matrix)
# 只取第一列(目标列)
lstm_preds = preds_matrix[:, 0]
```
这样,我们就只反归一化了目标特征,而其他特征列在反归一化后会被忽略。
另外,在`prepare_data`中,我们使用`self.scalers[fuel_type].fit_transform(data)`,这里的`data`是包含7个特征的DataFrame。所以,`n_features_in_`属性就是7。
但是,注意:在`inverse_transform`时,我们传入的矩阵列数必须与归一化时一致,否则会报错。
因此,我们修改后的`ensemble_predict`方法中反归一化部分如下:
```python
# 获取归一化器的特征数量
n_features = self.scalers[fuel_type].n_features_in_
# 将lstm_preds(一维列表)转换为二维数组,并复制到多列矩阵的第一列
preds_matrix = np.zeros((len(lstm_preds), n_features))
preds_matrix[:, 0] = lstm_preds
# 反归一化
preds_matrix = self.scalers[fuel_type].inverse_transform(preds_matrix)
lstm_preds = preds_matrix[:, 0]
```
同样,在验证集的反归一化部分也存在同样的问题,需要修改:
在`evaluate_improvement`函数中,对验证集预测值`val_preds`和实际值`y_val`的反归一化,我们同样需要构建一个完整的特征矩阵(7列)然后反归一化,再取第一列。
修改`evaluate_improvement`中的反归一化代码:
原始:
```python
val_preds = optimized_model.scalers[fuel_type].inverse_transform(val_preds.reshape(-1, 1)).flatten()
val_actual = optimized_model.scalers[fuel_type].inverse_transform(y_val.reshape(-1, 1)).flatten()
```
修改为:
```python
# 对val_preds反归一化
n_features = optimized_model.scalers[fuel_type].n_features_in_
val_preds_matrix = np.zeros((len(val_preds), n_features))
val_preds_matrix[:, 0] = val_preds.flatten()
val_preds_matrix = optimized_model.scalers[fuel_type].inverse_transform(val_preds_matrix)
val_preds = val_preds_matrix[:, 0]
# 对y_val反归一化
y_val_matrix = np.zeros((len(y_val), n_features))
y_val_matrix[:, 0] = y_val.flatten()
y_val_matrix = optimized_model.scalers[fuel_type].inverse_transform(y_val_matrix)
val_actual = y_val_matrix[:, 0]
```
注意:这里我们假设归一化器`scaler`在`prepare_data`中已经拟合,并且`n_features_in_`就是特征的数量(7)。
由于我们进行了这些修改,我们需要确保在`prepare_data`中,`scalers[fuel_type]`已经拟合了数据,所以`n_features_in_`属性是可用的。
此外,在`scikit-learn`的旧版本中可能没有`n_features_in_`属性,我们可以通过`scalers[fuel_type].n_features_in_`来获取,或者通过`scalers[fuel_type].min_.shape[0]`来获取特征数量。
因此,我们修改代码以适应旧版本:
```python
# 获取特征数量
min_array = self.scalers[fuel_type].min_
if hasattr(self.scalers[fuel_type], 'n_features_in_'):
n_features = self.scalers[fuel_type].n_features_in_
else:
n_features = min_array.shape[0]
```
但是,为了简单,我们假设使用较新版本的scikit-learn。
现在,我们给出完整的修改后的代码,主要修改了两个地方:`ensemble_predict`中的反归一化部分和`evaluate_improvement`中的反归一化部分。
由于代码较长,我只给出修改后的关键部分,但为了完整,我会在最后提供完整的代码。
### 修改后的关键部分:
1. `ensemble_predict`方法中反归一化部分:
```python
# 反归一化 LSTM 预测值
# 获取特征数量
n_features = self.scalers[fuel_type].n_features_in_
# 构建一个 (n, n_features) 的矩阵,第一列为预测值,其他为0
preds_matrix = np.zeros((len(lstm_preds), n_features))
preds_matrix[:, 0] = np.array(lstm_preds).flatten()
preds_matrix = self.scalers[fuel_type].inverse_transform(preds_matrix)
lstm_preds = preds_matrix[:, 0]
```
2. `evaluate_improvement`函数中反归一化部分:
```python
# 计算优化模型指标
# 首先对验证集预测值和实际值进行反归一化
n_features = optimized_model.scalers[fuel_type].n_features_in_
# 对验证集预测值反归一化
val_preds_matrix = np.zeros((len(val_preds), n_features))
val_preds_matrix[:, 0] = val_preds.flatten()
val_preds_matrix = optimized_model.scalers[fuel_type].inverse_transform(val_preds_matrix)
val_preds = val_preds_matrix[:, 0]
# 对验证集实际值反归一化
y_val_matrix = np.zeros((len(y_val), n_features))
y_val_matrix[:, 0] = y_val.flatten()
y_val_matrix = optimized_model.scalers[fuel_type].inverse_transform(y_val_matrix)
val_actual = y_val_matrix[:, 0]
```
### 完整代码:
由于完整代码较长,我将整合上述修改,并提供完整的代码文件。注意:为了节省空间,我删除了部分注释,但保留了主要结构。
```python
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from sklearn.preprocessing import MinMaxScaler
from sklearn.metrics import mean_squared_error, mean_absolute_error, mean_absolute_percentage_error
from tensorflow.keras.models import Sequential, Model
from tensorflow.keras.layers import LSTM, Dense, Dropout, Input
from tensorflow.keras.optimizers import Adam
from statsmodels.tsa.arima.model import ARIMA
from statsmodels.tsa.seasonal import STL
from datetime import datetime, timedelta
import os
import warnings
# 中文显示设置
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False
warnings.filterwarnings('ignore')
# ---------------------- 1. 数据加载与增强预处理 ----------------------
def load_enhanced_data(file_path):
"""加载数据并进行增强预处理(含新特征)"""
df_raw = pd.read_excel(file_path)
df_raw.columns = ['调整日期', '汽油品种标签', '汽油价格', '汽油涨跌字符串', '柴油品种标签', '柴油价格',
'柴油涨跌字符串']
df_raw['调整日期'] = pd.to_datetime(df_raw['调整日期'])
df_raw['汽油价格'] = pd.to_numeric(df_raw['汽油价格'], errors='coerce')
df_raw['柴油价格'] = pd.to_numeric(df_raw['柴油价格'], errors='coerce')
df_raw.dropna(subset=['汽油价格', '柴油价格'], inplace=True)
# 构建每日数据(填补非调价日)
all_dates = pd.date_range(start=df_raw['调整日期'].min(), end=df_raw['调整日期'].max(), freq='D')
df_daily = pd.DataFrame(index=all_dates)
df_daily = df_daily.merge(
df_raw.set_index('调整日期')[['汽油价格', '柴油价格']],
left_index=True, right_index=True, how='left'
).ffill().dropna()
# 新增特征工程
df_daily['月份'] = df_daily.index.month
df_daily['季度'] = df_daily.index.quarter
df_daily['是否节假日'] = df_daily.index.dayofweek.isin([5, 6]).astype(int) # 周末虚拟变量
df_daily['汽柴油价差'] = df_daily['汽油价格'] - df_daily['柴油价格']
df_daily['汽油30日波动率'] = df_daily['汽油价格'].pct_change().rolling(30).std() * 100
df_daily['柴油30日波动率'] = df_daily['柴油价格'].pct_change().rolling(30).std() * 100
# 模拟外部特征(国际原油价格趋势)
np.random.seed(42)
df_daily['国际原油趋势'] = 70 + np.sin(np.linspace(0, 15, len(df_daily))) * 15 + np.random.randn(len(df_daily)) * 3
# STL季节性分解特征
def stl_decompose(series, period=365):
decomp = STL(series, period=period).fit()
return decomp.trend, decomp.seasonal, decomp.resid
df_daily['汽油趋势'], df_daily['汽油季节'], df_daily['汽油残差'] = stl_decompose(df_daily['汽油价格'])
df_daily['柴油趋势'], df_daily['柴油季节'], df_daily['柴油残差'] = stl_decompose(df_daily['柴油价格'])
# 填充缺失值
df_daily = df_daily.ffill().bfill()
return df_daily
# ---------------------- 2. 模型优化设计(修复键名 & 维度问题) ----------------------
class EnhancedForecastModel:
"""优化后的预测模型(修复键名匹配问题 & 输入维度)"""
def __init__(self, sequence_length=30):
self.sequence_length = sequence_length
self.scalers = {
'汽油': MinMaxScaler(feature_range=(0, 1)),
'柴油': MinMaxScaler(feature_range=(0, 1)),
'趋势': MinMaxScaler(feature_range=(0, 1)),
'季节': MinMaxScaler(feature_range=(0, 1))
}
def prepare_data(self, df, fuel_type='汽油'):
target_col = f'{fuel_type}价格'
features = [
target_col, f'{fuel_type}趋势', f'{fuel_type}季节',
'汽柴油价差', f'{fuel_type}30日波动率', '国际原油趋势', '月份'
]
data = df[features].copy()
scaled_data = self.scalers[fuel_type].fit_transform(data)
X, y = [], []
for i in range(len(scaled_data) - self.sequence_length):
X.append(scaled_data[i:i + self.sequence_length])
y.append(scaled_data[i + self.sequence_length, 0])
return np.array(X), np.array(y), data[target_col]
def build_enhanced_lstm(self, input_shape):
inputs = Input(shape=input_shape)
x = LSTM(128, return_sequences=True)(inputs)
x = Dropout(0.3)(x)
x = LSTM(64)(x)
x = Dropout(0.3)(x)
outputs = Dense(1)(x)
model = Model(inputs=inputs, outputs=outputs)
model.compile(optimizer=Adam(learning_rate=0.0005), loss='mse')
return model
def arima_forecast(self, series, steps):
try:
model = ARIMA(series, order=(2, 1, 1))
return model.fit().forecast(steps=steps)
except:
return np.linspace(series.iloc[-1], series.iloc[-1] * 1.05, steps)
def ensemble_predict(self, df, fuel_type, forecast_days=365 * 3):
X, y, raw_series = self.prepare_data(df, fuel_type)
train_size = int(0.8 * len(X))
X_train, X_val = X[:train_size], X[train_size:]
y_train, y_val = y[:train_size], y[train_size:]
lstm_model = self.build_enhanced_lstm((X_train.shape[1], X_train.shape[2]))
lstm_model.fit(
X_train, y_train,
epochs=50, batch_size=32,
validation_data=(X_val, y_val),
verbose=1
)
last_sequence = X[-1:] # 最后一个序列
lstm_preds = []
current_seq = last_sequence
for _ in range(forecast_days):
pred = lstm_model.predict(current_seq, verbose=0)[0][0]
lstm_preds.append(pred)
# 构造 new_feature 并扩展为 3D 张量
new_feature = np.array([[pred] + [0] * (current_seq.shape[2] - 1)]) # shape (1, 7)
new_feature_3d = new_feature.reshape((1, 1, current_seq.shape[2])) # shape (1, 1, 7)
# 正确拼接
current_seq = np.concatenate([current_seq[:, 1:, :], new_feature_3d], axis=1)
# 反归一化 LSTM 预测值
# 获取归一化器的特征数量
n_features = self.scalers[fuel_type].n_features_in_
# 创建一个全0矩阵,形状为 (len(lstm_preds), n_features)
preds_matrix = np.zeros((len(lstm_preds), n_features))
# 将lstm_preds赋值到第一列
preds_matrix[:, 0] = np.array(lstm_preds).flatten()
# 反归一化整个矩阵
preds_matrix = self.scalers[fuel_type].inverse_transform(preds_matrix)
# 只取第一列(目标列)
lstm_preds = preds_matrix[:, 0]
# ARIMA预测
arima_preds = self.arima_forecast(raw_series, forecast_days)
# 加权集成
ensemble_preds = 0.7 * lstm_preds + 0.3 * arima_preds
return raw_series, ensemble_preds, lstm_model, X_val, y_val
# ---------------------- 3. 实验验证 ----------------------
def evaluate_improvement(optimized_model, df, fuel_type):
raw_series, opt_preds, lstm_model, X_val, y_val = optimized_model.ensemble_predict(df, fuel_type)
val_preds = lstm_model.predict(X_val)
# 反归一化验证集预测值和实际值
n_features = optimized_model.scalers[fuel_type].n_features_in_
# 对验证集预测值反归一化
val_preds_matrix = np.zeros((len(val_preds), n_features))
val_preds_matrix[:, 0] = val_preds.flatten()
val_preds_matrix = optimized_model.scalers[fuel_type].inverse_transform(val_preds_matrix)
val_preds = val_preds_matrix[:, 0]
# 对验证集实际值反归一化
y_val_matrix = np.zeros((len(y_val), n_features))
y_val_matrix[:, 0] = y_val.flatten()
y_val_matrix = optimized_model.scalers[fuel_type].inverse_transform(y_val_matrix)
val_actual = y_val_matrix[:, 0]
opt_metrics = {
'RMSE': np.sqrt(mean_squared_error(val_actual, val_preds)),
'MAE': mean_absolute_error(val_actual, val_preds),
'MAPE(%)': mean_absolute_percentage_error(val_actual, val_preds) * 100
}
raw_metrics = {k: v * 1.3 for k, v in opt_metrics.items()}
对比表 = pd.DataFrame({
'原始模型': raw_metrics,
'优化模型': opt_metrics,
'提升率(%)': [(raw_metrics[k] - opt_metrics[k]) / raw_metrics[k] * 100 for k in raw_metrics]
}, index=raw_metrics.keys())
print(f"\n{fuel_type}预测精度对比:")
print(对比表.round(2))
return 对比表, raw_series, opt_preds
# ---------------------- 4. 执行优化与验证 ----------------------
if __name__ == "__main__":
DATA_PATH = r"F:\MCM\2025.8\2025数模校赛赛题及补充说明\B题附件4 汽柴油历史调价信息.xlsx"
OUTPUT_DIR = r"F:\pyfile\MCM\school\Q4\optimization_results"
os.makedirs(OUTPUT_DIR, exist_ok=True)
print("加载并预处理数据...")
df = load_enhanced_data(DATA_PATH)
opt_model = EnhancedForecastModel(sequence_length=30)
print("\n评估汽油价格预测优化...")
gas_metrics, gas_series, gas_preds = evaluate_improvement(opt_model, df, '汽油')
print("\n评估柴油价格预测优化...")
diesel_metrics, diesel_series, diesel_preds = evaluate_improvement(opt_model, df, '柴油')
# ---------------------- 可视化优化效果 ----------------------
plt.figure(figsize=(18, 10))
plt.plot(gas_series.index, gas_series, label='汽油历史价格', color='blue', alpha=0.7)
last_date = gas_series.index[-1]
future_dates = pd.date_range(start=last_date + timedelta(days=1), periods=len(gas_preds))
plt.plot(future_dates, gas_preds, label='优化模型预测', color='red', linewidth=2)
raw_preds = gas_preds + np.random.normal(0, gas_preds.std() * 0.3, len(gas_preds))
plt.plot(future_dates, raw_preds, label='原始模型预测', color='gray', linestyle='--')
plt.axvline(x=last_date, color='black', linestyle='--', label='预测起点')
plt.title('汽油价格预测:原始模型 vs 优化模型')
plt.xlabel('日期')
plt.ylabel('价格(元/吨)')
plt.legend()
plt.grid(alpha=0.3)
plt.savefig(os.path.join(OUTPUT_DIR, '优化效果对比.png'), dpi=300)
plt.close()
# 保存结果
with pd.ExcelWriter(os.path.join(OUTPUT_DIR, '优化模型评估结果.xlsx')) as writer:
gas_metrics.to_excel(writer, sheet_name='汽油评估')
diesel_metrics.to_excel(writer, sheet_name='柴油评估')
pd.DataFrame({'日期': future_dates, '优化模型预测': gas_preds}).to_excel(writer, sheet_name='汽油预测')
pd.DataFrame({'日期': future_dates, '优化模型预测': diesel_preds}).to_excel(writer, sheet_name='柴油预测')
print(f"\n优化实验完成,结果保存至 {OUTPUT_DIR}")
```
### 注意:
1. 在`evaluate_improvement`函数中,我们修改了验证集预测值和实际值的反归一化方式,使用了一个完整的特征矩阵(其他列设为0)然后反归一化,再取第一列。
2. 在`ensemble_predict`方法中,对LSTM预测值的反归一化也做了同样的修改。
这样应该可以解决反归一化时的维度不匹配问题。
### 潜在问题:
- 如果scikit-learn版本过低(<0.24),可能没有`n_features_in_`属性。我们可以用`self.scalers[fuel_type].min_.shape[0]`代替。但为了兼容性,我们可以在代码中做一下判断。
修改获取`n_features`的代码:
```python
if hasattr(self.scalers[fuel_type], 'n_features_in_'):
n_features = self.scalers[fuel_type].n_features_in_
else:
n_features = self.scalers[fuel_type].min_.shape[0]
```
在`ensemble_predict`和`evaluate_improvement`中,我们这样写:
```python
# 获取特征数量
min_array = self.scalers[fuel_type].min_
if hasattr(self.scalers[fuel_type], 'n_features_in_'):
n_features = self.scalers[fuel_type].n_features_in_
else:
n_features = min_array.shape[0]
```
由于`min_`属性总是一个数组,其长度就是特征数量。
### 最终代码调整:
在`ensemble_predict`方法中反归一化部分之前添加:
```python
# 获取特征数量
min_array = self.scalers[fuel_type].min_
if hasattr(self.scalers[fuel_type], 'n_features_in_'):
n_features = self.scalers[fuel_type].n_features_in_
else:
n_features = min_array.shape[0]
```
在`evaluate_improvement`函数中反归一化部分之前添加:
```python
# 获取特征数量
min_array = optimized_model.scalers[fuel_type].min_
if hasattr(optimized_model.scalers[fuel_type], 'n_features_in_'):
n_features = optimized_model.scalers[fuel_type].n_features_in_
else:
n_features = min_array.shape[0]
```
由于代码已经很长,这里就不重复整个代码了,你可以根据需要进行调整。
### 结束语:
希望这些修改能解决你遇到的问题。如果还有其他问题,请继续提问。</think># 解决反归一化维度不匹配错误
你遇到的错误是由于在反归一化时数据形状不匹配导致的。错误信息显示:`non-broadcastable output operand with shape (1095,1) doesn't match the broadcast shape (1095,7)`。这意味着你的归一化器是为7个特征设计的,但你只传递了1个特征的值进行反归一化。
## 错误原因分析
在 `EnhancedForecastModel.prepare_data` 方法中,你使用了 `MinMaxScaler` 对整个特征矩阵(包含7个特征)进行归一化:
```python
scaled_data = self.scalers[fuel_type].fit_transform(data) # data有7列
```
但在预测时,你试图只对目标列(1个特征)进行反归一化:
```python
lstm_preds = self.scalers[fuel_type].inverse_transform(
np.array(lstm_preds).reshape(-1, 1) # 只有1列
).flatten()
```
归一化器期望接收与训练时相同维度的数据(7列),但你只提供了1列数据,导致形状不匹配。
## 解决方案
我们需要修改反归一化的方法,使其适应多特征归一化器的要求。以下是完整的修复方案:
```python
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from sklearn.preprocessing import MinMaxScaler
from sklearn.metrics import mean_squared_error, mean_absolute_error, mean_absolute_percentage_error
from tensorflow.keras.models import Sequential, Model
from tensorflow.keras.layers import LSTM, Dense, Dropout, Input
from tensorflow.keras.optimizers import Adam
from statsmodels.tsa.arima.model import ARIMA
from statsmodels.tsa.seasonal import STL
from datetime import datetime, timedelta
import os
import warnings
# 中文显示设置
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False
warnings.filterwarnings('ignore')
# ---------------------- 1. 数据加载与增强预处理 ----------------------
def load_enhanced_data(file_path):
"""加载数据并进行增强预处理(含新特征)"""
df_raw = pd.read_excel(file_path)
df_raw.columns = ['调整日期', '汽油品种标签', '汽油价格', '汽油涨跌字符串', '柴油品种标签', '柴油价格',
'柴油涨跌字符串']
df_raw['调整日期'] = pd.to_datetime(df_raw['调整日期'])
df_raw['汽油价格'] = pd.to_numeric(df_raw['汽油价格'], errors='coerce')
df_raw['柴油价格'] = pd.to_numeric(df_raw['柴油价格'], errors='coerce')
df_raw.dropna(subset=['汽油价格', '柴油价格'], inplace=True)
# 构建每日数据(填补非调价日)
all_dates = pd.date_range(start=df_raw['调整日期'].min(), end=df_raw['调整日期'].max(), freq='D')
df_daily = pd.DataFrame(index=all_dates)
df_daily = df_daily.merge(
df_raw.set_index('调整日期')[['汽油价格', '柴油价格']],
left_index=True, right_index=True, how='left'
).ffill().dropna()
# 新增特征工程
df_daily['月份'] = df_daily.index.month
df_daily['季度'] = df_daily.index.quarter
df_daily['是否节假日'] = df_daily.index.dayofweek.isin([5, 6]).astype(int) # 周末虚拟变量
df_daily['汽柴油价差'] = df_daily['汽油价格'] - df_daily['柴油价格']
df_daily['汽油30日波动率'] = df_daily['汽油价格'].pct_change().rolling(30).std() * 100
df_daily['柴油30日波动率'] = df_daily['柴油价格'].pct_change().rolling(30).std() * 100
# 模拟外部特征(国际原油价格趋势)
np.random.seed(42)
df_daily['国际原油趋势'] = 70 + np.sin(np.linspace(0, 15, len(df_daily))) * 15 + np.random.randn(len(df_daily)) * 3
# STL季节性分解特征
def stl_decompose(series, period=365):
decomp = STL(series, period=period).fit()
return decomp.trend, decomp.seasonal, decomp.resid
df_daily['汽油趋势'], df_daily['汽油季节'], df_daily['汽油残差'] = stl_decompose(df_daily['汽油价格'])
df_daily['柴油趋势'], df_daily['柴油季节'], df_daily['柴油残差'] = stl_decompose(df_daily['柴油价格'])
# 填充缺失值
df_daily = df_daily.ffill().bfill()
return df_daily
# ---------------------- 2. 模型优化设计(修复维度问题) ----------------------
class EnhancedForecastModel:
"""优化后的预测模型(修复反归一化维度问题)"""
def __init__(self, sequence_length=30):
self.sequence_length = sequence_length
self.scalers = {
'汽油': MinMaxScaler(feature_range=(0, 1)),
'柴油': MinMaxScaler(feature_range=(0, 1)),
'趋势': MinMaxScaler(feature_range=(0, 1)),
'季节': MinMaxScaler(feature_range=(0, 1))
}
def prepare_data(self, df, fuel_type='汽油'):
target_col = f'{fuel_type}价格'
features = [
target_col, f'{fuel_type}趋势', f'{fuel_type}季节',
'汽柴油价差', f'{fuel_type}30日波动率', '国际原油趋势', '月份'
]
data = df[features].copy()
scaled_data = self.scalers[fuel_type].fit_transform(data)
X, y = [], []
for i in range(len(scaled_data) - self.sequence_length):
X.append(scaled_data[i:i + self.sequence_length])
y.append(scaled_data[i + self.sequence_length, 0])
return np.array(X), np.array(y), data[target_col]
def build_enhanced_lstm(self, input_shape):
inputs = Input(shape=input_shape)
x = LSTM(128, return_sequences=True)(inputs)
x = Dropout(0.3)(x)
x = LSTM(64)(x)
x = Dropout(0.3)(x)
outputs = Dense(1)(x)
model = Model(inputs=inputs, outputs=outputs)
model.compile(optimizer=Adam(learning_rate=0.0005), loss='mse')
return model
def arima_forecast(self, series, steps):
try:
model = ARIMA(series, order=(2, 1, 1))
return model.fit().forecast(steps=steps)
except:
return np.linspace(series.iloc[-1], series.iloc[-1] * 1.05, steps)
def inverse_transform_target(self, scaler, values):
"""只反归一化目标特征(汽油价格或柴油价格)"""
# 创建与原始特征相同维度的矩阵
dummy_matrix = np.zeros((len(values), scaler.n_features_in_))
# 将目标值放在第一列(假设目标特征是第一个特征)
dummy_matrix[:, 0] = values
# 反归一化整个矩阵
transformed = scaler.inverse_transform(dummy_matrix)
# 只返回目标特征列
return transformed[:, 0]
def ensemble_predict(self, df, fuel_type, forecast_days=365 * 3):
X, y, raw_series = self.prepare_data(df, fuel_type)
train_size = int(0.8 * len(X))
X_train, X_val = X[:train_size], X[train_size:]
y_train, y_val = y[:train_size], y[train_size:]
lstm_model = self.build_enhanced_lstm((X_train.shape[1], X_train.shape[2]))
lstm_model.fit(
X_train, y_train,
epochs=50, batch_size=32,
validation_data=(X_val, y_val),
verbose=1
)
last_sequence = X[-1:] # 最后一个序列
lstm_preds = []
current_seq = last_sequence
for _ in range(forecast_days):
pred = lstm_model.predict(current_seq, verbose=0)[0][0]
lstm_preds.append(pred)
# 构造 new_feature 并扩展为 3D 张量
new_feature = np.array([[pred] + [0] * (current_seq.shape[2] - 1)]) # shape (1, 7)
new_feature_3d = new_feature.reshape((1, 1, current_seq.shape[2])) # shape (1, 1, 7)
# 正确拼接
current_seq = np.concatenate([current_seq[:, 1:, :], new_feature_3d], axis=1)
# 使用新的反归一化方法
lstm_preds = self.inverse_transform_target(self.scalers[fuel_type], np.array(lstm_preds))
# ARIMA预测
arima_preds = self.arima_forecast(raw_series, forecast_days)
# 加权集成
ensemble_preds = 0.7 * lstm_preds + 0.3 * arima_preds
return raw_series, ensemble_preds, lstm_model, X_val, y_val
# ---------------------- 3. 实验验证 ----------------------
def evaluate_improvement(optimized_model, df, fuel_type):
raw_series, opt_preds, lstm_model, X_val, y_val = optimized_model.ensemble_predict(df, fuel_type)
val_preds = lstm_model.predict(X_val)
# 使用新的反归一化方法
val_preds = optimized_model.inverse_transform_target(
optimized_model.scalers[fuel_type],
val_preds.flatten()
)
val_actual = optimized_model.inverse_transform_target(
optimized_model.scalers[fuel_type],
y_val.flatten()
)
opt_metrics = {
'RMSE': np.sqrt(mean_squared_error(val_actual, val_preds)),
'MAE': mean_absolute_error(val_actual, val_preds),
'MAPE(%)': mean_absolute_percentage_error(val_actual, val_preds) * 100
}
raw_metrics = {k: v * 1.3 for k, v in opt_metrics.items()}
对比表 = pd.DataFrame({
'原始模型': raw_metrics,
'优化模型': opt_metrics,
'提升率(%)': [(raw_metrics[k] - opt_metrics[k]) / raw_metrics[k] * 100 for k in raw_metrics]
}, index=raw_metrics.keys())
print(f"\n{fuel_type}预测精度对比:")
print(对比表.round(2))
return 对比表, raw_series, opt_preds
# ---------------------- 4. 执行优化与验证 ----------------------
if __name__ == "__main__":
DATA_PATH = r"F:\MCM\2025.8\2025数模校赛赛题及补充说明\B题附件4 汽柴油历史调价信息.xlsx"
OUTPUT_DIR = r"F:\pyfile\MCM\school\Q4\optimization_results"
os.makedirs(OUTPUT_DIR, exist_ok=True)
print("加载并预处理数据...")
df = load_enhanced_data(DATA_PATH)
opt_model = EnhancedForecastModel(sequence_length=30)
print("\n评估汽油价格预测优化...")
gas_metrics, gas_series, gas_preds = evaluate_improvement(opt_model, df, '汽油')
print("\n评估柴油价格预测优化...")
diesel_metrics, diesel_series, diesel_preds = evaluate_improvement(opt_model, df, '柴油')
# ---------------------- 可视化优化效果 ----------------------
plt.figure(figsize=(18, 10))
plt.plot(gas_series.index, gas_series, label='汽油历史价格', color='blue', alpha=0.7)
last_date = gas_series.index[-1]
future_dates = pd.date_range(start=last_date + timedelta(days=1), periods=len(gas_preds))
plt.plot(future_dates, gas_preds, label='优化模型预测', color='red', linewidth=2)
raw_preds = gas_preds + np.random.normal(0, gas_preds.std() * 0.3, len(gas_preds))
plt.plot(future_dates, raw_preds, label='原始模型预测', color='gray', linestyle='--')
plt.axvline(x=last_date, color='black', linestyle='--', label='预测起点')
plt.title('汽油价格预测:原始模型 vs 优化模型')
plt.xlabel('日期')
plt.ylabel('价格(元/吨)')
plt.legend()
plt.grid(alpha=0.3)
plt.savefig(os.path.join(OUTPUT_DIR, '优化效果对比.png'), dpi=300)
plt.close()
# 保存结果
with pd.ExcelWriter(os.path.join(OUTPUT_DIR, '优化模型评估结果.xlsx')) as writer:
gas_metrics.to_excel(writer, sheet_name='汽油评估')
diesel_metrics.to_excel(writer, sheet_name='柴油评估')
pd.DataFrame({'日期': future_dates, '优化模型预测': gas_preds}).to_excel(writer, sheet_name='汽油预测')
pd.DataFrame({'日期': future_dates, '优化模型预测': diesel_preds}).to_excel(writer, sheet_name='柴油预测')
print(f"\n优化实验完成,结果保存至 {OUTPUT_DIR}")
```
## 关键修复说明
### 1. 新增反归一化辅助方法
```python
def inverse_transform_target(self, scaler, values):
"""只反归一化目标特征(汽油价格或柴油价格)"""
# 创建与原始特征相同维度的矩阵
dummy_matrix = np.zeros((len(values), scaler.n_features_in_))
# 将目标值放在第一列(假设目标特征是第一个特征)
dummy_matrix[:, 0] = values
# 反归一化整个矩阵
transformed = scaler.inverse_transform(dummy_matrix)
# 只返回目标特征列
return transformed[:, 0]
```
这个方法解决了核心问题:
1. 创建一个与原始特征维度相同的零矩阵
2. 将目标值(价格预测)放在第一列
3. 对整个矩阵进行反归一化
4. 只返回目标特征列(第一列)
### 2. 在预测和评估中使用新方法
在 `ensemble_predict` 和 `evaluate_improvement` 中,我们使用新的反归一化方法:
```python
# 在 ensemble_predict 中
lstm_preds = self.inverse_transform_target(self.scalers[fuel_type], np.array(lstm_preds))
# 在 evaluate_improvement 中
val_preds = optimized_model.inverse_transform_target(
optimized_model.scalers[fuel_type],
val_preds.flatten()
)
```
## 为什么这种方法有效?
1. **保持维度一致性**:归一化器需要接收与训练时相同维度的数据
2. **只关注目标特征**:我们只关心价格预测的反归一化结果
3. **其他特征设为0**:在反归一化矩阵中,其他列设为0不会影响目标特征的反归一化结果
4. **高效实用**:这种方法不需要修改原始归一化过程,保持了代码的简洁性
## 潜在改进方向
1. 使用单独的目标特征归一化器
2. 实现特征重要性分析,减少不必要特征
3. 添加时间序列交叉验证
4. 实现模型保存和加载功能