一、系统架构
%% 光伏出力预测系统 - 考虑气象因子的深度学习框架
% 系统功能：基于历史数据和气象因子的光伏出力预测
% 预测方法：LSTM + 注意力机制 + 多变量时间序列预测
% 气象因子：温度、湿度、风速、辐照度、云量、气压等
clear; clc; close all;
warning('off', 'all');
%% 1. 系统参数配置
disp('=== 光伏出力预测系统 ===');
disp('初始化系统参数...');
% 预测参数
params = struct();
params.pred_horizon = 24; % 预测未来24小时
params.lookback_window = 168; % 使用过去168小时（7天）数据
params.train_ratio = 0.8; % 训练集比例
params.val_ratio = 0.1; % 验证集比例
params.batch_size = 32; % 批次大小
params.max_epochs = 100; % 最大训练轮数
params.learning_rate = 0.001; % 学习率
params.dropout_rate = 0.2; % Dropout率
params.patience = 10; % 早停耐心值
% 数据参数
params.seasonal_period = 24; % 日周期
params.trend_period = 168; % 周周期
params.n_features = 12; % 特征数量
% 模型选择
params.model_type = 'lstm_attention'; % lstm, bilstm, lstm_attention, cnn_lstm
params.use_weather = true; % 是否使用气象因子
params.use_temporal = true; % 是否使用时序特征
params.normalize_method = 'zscore'; % 标准化方法: zscore, minmax
%% 2. 数据加载与预处理
disp('加载和预处理数据...');
% 2.1 加载数据（模拟或真实数据）
if exist('real_pv_data.mat', 'file')
load('real_pv_data.mat', 'pv_data', 'weather_data', 'timestamp');
fprintf('已加载真实数据: %d 个样本\n', length(pv_data));
else
[pv_data, weather_data, timestamp] = generate_synthetic_data(params);
fprintf('已生成模拟数据: %d 个样本\n', length(pv_data));
end
% 2.2 特征工程
features = create_features(pv_data, weather_data, timestamp, params);
% 2.3 数据标准化
[features_norm, norm_params] = normalize_features(features, params);
% 2.4 划分训练集、验证集、测试集
[X_train, Y_train, X_val, Y_val, X_test, Y_test, idx] = ...
split_dataset(features_norm, pv_data, params);
fprintf('数据划分结果:\n');
fprintf(' 训练集: %d 个样本\n', size(X_train, 1));
fprintf(' 验证集: %d 个样本\n', size(X_val, 1));
fprintf(' 测试集: %d 个样本\n', size(X_test, 1));
%% 3. 模型构建
disp('构建深度学习模型...');
input_size = size(X_train, 2);
num_features = size(X_train, 3);
switch params.model_type
case 'lstm'
model = create_lstm_model(input_size, num_features, params);
case 'bilstm'
model = create_bilstm_model(input_size, num_features, params);
case 'lstm_attention'
model = create_lstm_attention_model(input_size, num_features, params);
case 'cnn_lstm'
model = create_cnn_lstm_model(input_size, num_features, params);
otherwise
error('未知的模型类型');
end
% 显示模型结构
analyzeNetwork(model);
%% 4. 模型训练
disp('开始训练模型...');
% 训练选项
options = trainingOptions('adam', ...
'MaxEpochs', params.max_epochs, ...
'MiniBatchSize', params.batch_size, ...
'InitialLearnRate', params.learning_rate, ...
'GradientThreshold', 1, ...
'Shuffle', 'every-epoch', ...
'ValidationData', {X_val, Y_val}, ...
'ValidationFrequency', 30, ...
'ValidationPatience', params.patience, ...
'Plots', 'training-progress', ...
'Verbose', false, ...
'LearnRateSchedule', 'piecewise', ...
'LearnRateDropFactor', 0.5, ...
'LearnRateDropPeriod', 20);
% 训练模型
[net, info] = trainNetwork(X_train, Y_train, model, options);
% 保存模型
save('pv_prediction_model.mat', 'net', 'norm_params', 'params', 'info');
%% 5. 模型评估
disp('评估模型性能...');
% 5.1 在验证集上预测
Y_val_pred = predict(net, X_val, 'MiniBatchSize', params.batch_size);
% 5.2 在测试集上预测
Y_test_pred = predict(net, X_test, 'MiniBatchSize', params.batch_size);
% 5.3 反标准化
Y_test_true = Y_test * norm_params.pv_std + norm_params.pv_mean;
Y_test_pred_denorm = Y_test_pred * norm_params.pv_std + norm_params.pv_mean;
% 5.4 计算评价指标
metrics = calculate_metrics(Y_test_true, Y_test_pred_denorm);
fprintf('\n=== 模型性能指标 ===\n');
fprintf('RMSE: %.4f kW\n', metrics.rmse);
fprintf('MAE: %.4f kW\n', metrics.mae);
fprintf('MAPE: %.2f%%\n', metrics.mape * 100);
fprintf('R²: %.4f\n', metrics.r2);
fprintf('nRMSE: %.2f%%\n', metrics.nrmse * 100);
%% 6. 结果可视化
disp('生成可视化结果...');
visualize_results(Y_test_true, Y_test_pred_denorm, timestamp(idx.test), metrics, params);
%% 7. 特征重要性分析
disp('分析特征重要性...');
feature_importance = analyze_feature_importance(net, X_test, Y_test, features, params);
plot_feature_importance(feature_importance);
%% 8. 未来预测
disp('进行未来24小时预测...');
future_prediction = predict_future(net, features_norm, norm_params, params);
plot_future_prediction(future_prediction, timestamp);
%% 9. 敏感性分析
disp('进行敏感性分析...');
sensitivity_results = sensitivity_analysis(net, features_norm, norm_params, params);
plot_sensitivity_analysis(sensitivity_results);
%% 10. 模型解释性分析
disp('模型解释性分析...');
explain_model_predictions(net, X_test, Y_test_pred, features, params);
fprintf('\n=== 光伏出力预测完成 ===\n');
二、数据生成模块
function [pv_data, weather_data, timestamp] = generate_synthetic_data(params)
% 生成模拟光伏出力和气象数据
n_days = 365; % 一年数据
n_samples = n_days * 24; % 小时数据
% 生成时间戳
start_date = datetime('2023-01-01 00:00:00');
timestamp = start_date + hours(0:n_samples-1)';
% 1. 生成光伏出力数据（基于太阳辐射模型）
pv_data = zeros(n_samples, 1);
for i = 1:n_samples
t = timestamp(i);
hour_of_day = hour(t);
day_of_year = day(t, 'dayofyear');
% 太阳高度角计算
lat = 39.9; % 纬度（北京）
lon = 116.4; % 经度
[sun_altitude, ~] = calculate_sun_position(lat, lon, t);
if sun_altitude > 0
% 日出到日落期间有光伏出力
% 基本出力曲线（正弦平方）
base_power = 1000 * sind(sun_altitude)^2;
% 加入日变化
daily_variation = 1 + 0.2 * sin(2*pi*hour_of_day/24);
% 加入季节变化
seasonal_variation = 1 + 0.3 * sin(2*pi*(day_of_year-80)/365);
% 加入随机波动
random_variation = 1 + 0.1 * randn;
pv_data(i) = base_power * daily_variation * seasonal_variation * random_variation;
else
pv_data(i) = 0;
end
end
% 平滑处理
pv_data = smoothdata(pv_data, 'gaussian', 6);
% 2. 生成气象数据
weather_data = struct();
% 温度（℃） - 日变化和季节变化
base_temp = 15 + 10 * sin(2*pi*(day_of_year-105)/365);
daily_temp_variation = 8 * sin(2*pi*(hour_of_day-14)/24);
weather_data.temperature = base_temp + daily_temp_variation + 3*randn(n_samples,1);
% 湿度（%） - 与温度负相关
base_humidity = 60 - 0.5 * (weather_data.temperature - 20);
weather_data.humidity = max(20, min(100, base_humidity + 10*randn(n_samples,1)));
% 风速（m/s）
base_wind = 3 + 2 * sin(2*pi*day_of_year/365);
weather_data.wind_speed = max(0, base_wind + 2*randn(n_samples,1));
% 辐照度（W/m²） - 与光伏出力强相关
weather_data.irradiance = pv_data + 50*randn(n_samples,1);
weather_data.irradiance = max(0, weather_data.irradiance);
% 云量（0-1）
base_cloud = 0.3 + 0.2 * sin(2*pi*day_of_year/365);
weather_data.cloud_cover = max(0, min(1, base_cloud + 0.2*randn(n_samples,1)));
% 气压（hPa）
base_pressure = 1013 + 10 * sin(2*pi*day_of_year/365);
weather_data.pressure = base_pressure + 5*randn(n_samples,1);
% 降水（mm/h）
rain_prob = 0.1 + 0.05 * sin(2*pi*day_of_year/365);
is_raining = rand(n_samples,1) < rain_prob;
weather_data.precipitation = is_raining .* (0.1 + 5*rand(n_samples,1));
% 能见度（km）
base_visibility = 10 - 8 * weather_data.cloud_cover;
weather_data.visibility = max(1, base_visibility + 2*randn(n_samples,1));
fprintf('已生成 %d 天的模拟数据\n', n_days);
end
function [altitude, azimuth] = calculate_sun_position(lat, lon, time)
% 计算太阳位置
% 输入：纬度(deg), 经度(deg), 时间(datetime)
% 输出：太阳高度角(deg), 方位角(deg)
% 转换为弧度
lat_rad = deg2rad(lat);
% 计算儒略日
d = days(time - datetime(year(time),1,1)) + 1;
% 太阳赤纬（简化计算）
declination = 23.45 * sind(360*(284 + d)/365);
declination_rad = deg2rad(declination);
% 计算时角
hour_angle = 15 * (hour(time) + minute(time)/60 - 12);
hour_angle_rad = deg2rad(hour_angle);
% 计算太阳高度角
sin_altitude = sin(lat_rad) * sin(declination_rad) + ...
cos(lat_rad) * cos(declination_rad) * cos(hour_angle_rad);
altitude = asin(sin_altitude);
altitude = rad2deg(altitude);
% 计算太阳方位角
cos_azimuth = (sin(declination_rad) - sin(lat_rad) * sin(altitude)) / ...
(cos(lat_rad) * cos(asin(sin_altitude)));
azimuth = acos(min(1, max(-1, cos_azimuth)));
azimuth = rad2deg(azimuth);
% 调整方位角
if hour_angle > 0
azimuth = 360 - azimuth;
end
end
三、特征工程模块
function features = create_features(pv_data, weather_data, timestamp, params)
% 创建特征矩阵
n_samples = length(pv_data);
% 初始化特征矩阵
n_features = params.n_features;
features = zeros(n_samples, n_features);
% 1. 光伏出力相关特征
features(:,1) = pv_data; % 当前光伏出力
% 滞后特征
features(:,2) = [0; pv_data(1:end-1)]; % 滞后1小时
features(:,3) = [zeros(2,1); pv_data(1:end-2)]; % 滞后2小时
features(:,4) = [zeros(24,1); pv_data(1:end-24)]; % 滞后24小时
features(:,5) = [zeros(168,1); pv_data(1:end-168)]; % 滞后168小时
% 2. 气象特征
features(:,6) = weather_data.temperature;
features(:,7) = weather_data.humidity;
features(:,8) = weather_data.wind_speed;
features(:,9) = weather_data.irradiance;
features(:,10) = weather_data.cloud_cover;
% 3. 时间特征
hour_of_day = hour(timestamp);
day_of_week = weekday(timestamp);
month_of_year = month(timestamp);
day_of_year = day(timestamp, 'dayofyear');
% 周期性编码
features(:,11) = sin(2*pi*hour_of_day/24);
features(:,12) = cos(2*pi*hour_of_day/24);
% 4. 派生特征
% 温度-辐照度交互特征
features(:,13) = weather_data.temperature .* weather_data.irradiance;
% 云量对辐照度的影响
features(:,14) = weather_data.irradiance .* (1 - weather_data.cloud_cover);
% 滑动统计特征
window_size = 24;
features(:,15) = movmean(pv_data, window_size, 'omitnan');
features(:,16) = movstd(pv_data, window_size, 'omitnan');
% 5. 天气状态特征
% 晴天指数
sunshine_index = (1 - weather_data.cloud_cover) .* (weather_data.irradiance / max(weather_data.irradiance));
features(:,17) = sunshine_index;
% 降雨影响
rain_effect = exp(-weather_data.precipitation);
features(:,18) = rain_effect;
% 6. 季节性特征
features(:,19) = sin(2*pi*day_of_year/365);
features(:,20) = cos(2*pi*day_of_year/365);
% 7. 周特征
features(:,21) = sin(2*pi*day_of_week/7);
features(:,22) = cos(2*pi*day_of_week/7);
% 8. 气象组合特征
% 体感温度
apparent_temp = weather_data.temperature + 0.33*weather_data.humidity/100*6.105*exp(17.27*weather_data.temperature/(237.7+weather_data.temperature)) - 0.7*weather_data.wind_speed - 4;
features(:,23) = apparent_temp;
% 能见度指数
visibility_index = weather_data.visibility / max(weather_data.visibility);
features(:,24) = visibility_index;
% 更新特征数量
params.n_features = size(features, 2);
fprintf('已创建 %d 个特征\n', params.n_features);
end
四、数据预处理模块
function [features_norm, norm_params] = normalize_features(features, params)
% 数据标准化
n_samples = size(features, 1);
n_features = size(features, 2);
features_norm = zeros(size(features));
norm_params = struct();
switch params.normalize_method
case 'zscore'
% Z-score标准化
for i = 1:n_features
mu = mean(features(:,i), 'omitnan');
sigma = std(features(:,i), 'omitnan');
if sigma < 1e-6
sigma = 1;
end
features_norm(:,i) = (features(:,i) - mu) / sigma;
% 保存标准化参数
norm_params.mean(i) = mu;
norm_params.std(i) = sigma;
end
case 'minmax'
% 最小-最大标准化
for i = 1:n_features
f_min = min(features(:,i));
f_max = max(features(:,i));
if (f_max - f_min) < 1e-6
features_norm(:,i) = 0;
else
features_norm(:,i) = (features(:,i) - f_min) / (f_max - f_min);
end
% 保存标准化参数
norm_params.min(i) = f_min;
norm_params.max(i) = f_max;
end
otherwise
error('未知的标准化方法');
end
% 处理NaN值
features_norm(isnan(features_norm)) = 0;
% 保存光伏出力的标准化参数（用于反标准化）
norm_params.pv_mean = norm_params.mean(1);
norm_params.pv_std = norm_params.std(1);
fprintf('已完成 %s 标准化\n', params.normalize_method);
end
function [X_train, Y_train, X_val, Y_val, X_test, Y_test, idx] = ...
split_dataset(features, target, params)
% 划分数据集
n_samples = size(features, 1);
lookback = params.lookback_window;
horizon = params.pred_horizon;
% 计算可用样本数
n_usable = n_samples - lookback - horizon + 1;
% 生成输入输出对
X = zeros(n_usable, lookback, size(features, 2));
Y = zeros(n_usable, horizon);
for i = 1:n_usable
X(i,:,:) = features(i:i+lookback-1, :);
Y(i,:) = target(i+lookback:i+lookback+horizon-1, 1);
end
% 划分数据集
n_train = floor(size(X,1) * params.train_ratio);
n_val = floor(size(X,1) * params.val_ratio);
idx.train = 1:n_train;
idx.val = n_train+1:n_train+n_val;
idx.test = n_train+n_val+1:size(X,1);
X_train = X(idx.train,:,:);
Y_train = Y(idx.train,:);
X_val = X(idx.val,:,:);
Y_val = Y(idx.val,:);
X_test = X(idx.test,:,:);
Y_test = Y(idx.test,:);
end
五、模型构建模块
function model = create_lstm_attention_model(input_size, num_features, params)
% 创建LSTM+注意力机制模型
layers = [
% 输入层
sequenceInputLayer([input_size, num_features], 'Name', 'input')
% 第一个LSTM层
lstmLayer(128, 'OutputMode', 'sequence', 'Name', 'lstm1')
dropoutLayer(params.dropout_rate, 'Name', 'dropout1')
% 第二个LSTM层
lstmLayer(64, 'OutputMode', 'sequence', 'Name', 'lstm2')
dropoutLayer(params.dropout_rate, 'Name', 'dropout2')
% 注意力层
attentionLayer('Name', 'attention')
% 全连接层
fullyConnectedLayer(32, 'Name', 'fc1')
reluLayer('Name', 'relu1')
dropoutLayer(params.dropout_rate, 'Name', 'dropout3')
% 输出层
fullyConnectedLayer(params.pred_horizon, 'Name', 'fc_output')
regressionLayer('Name', 'output')
];
model = layerGraph(layers);
fprintf('已创建LSTM+注意力模型\n');
end
function model = create_bilstm_model(input_size, num_features, params)
% 创建双向LSTM模型
layers = [
% 输入层
sequenceInputLayer([input_size, num_features], 'Name', 'input')
% 双向LSTM层
bilstmLayer(128, 'OutputMode', 'sequence', 'Name', 'bilstm1')
dropoutLayer(params.dropout_rate, 'Name', 'dropout1')
% 第二个双向LSTM层
bilstmLayer(64, 'OutputMode', 'last', 'Name', 'bilstm2')
dropoutLayer(params.dropout_rate, 'Name', 'dropout2')
% 全连接层
fullyConnectedLayer(32, 'Name', 'fc1')
reluLayer('Name', 'relu1')
dropoutLayer(params.dropout_rate, 'Name', 'dropout3')
% 输出层
fullyConnectedLayer(params.pred_horizon, 'Name', 'fc_output')
regressionLayer('Name', 'output')
];
model = layerGraph(layers);
fprintf('已创建双向LSTM模型\n');
end
function model = create_cnn_lstm_model(input_size, num_features, params)
% 创建CNN-LSTM混合模型
layers = [
% 输入层
sequenceInputLayer([input_size, num_features], 'Name', 'input')
% 1D卷积层提取局部特征
convolution1dLayer(3, 32, 'Padding', 'same', 'Name', 'conv1')
reluLayer('Name', 'relu1')
maxPooling1dLayer(2, 'Stride', 2, 'Name', 'pool1')
% 第二个卷积层
convolution1dLayer(3, 64, 'Padding', 'same', 'Name', 'conv2')
reluLayer('Name', 'relu2')
maxPooling1dLayer(2, 'Stride', 2, 'Name', 'pool2')
% LSTM层捕捉时序依赖
flattenLayer('Name', 'flatten')
lstmLayer(128, 'OutputMode', 'last', 'Name', 'lstm1')
dropoutLayer(params.dropout_rate, 'Name', 'dropout1')
% 全连接层
fullyConnectedLayer(64, 'Name', 'fc1')
reluLayer('Name', 'relu3')
dropoutLayer(params.dropout_rate, 'Name', 'dropout2')
% 输出层
fullyConnectedLayer(params.pred_horizon, 'Name', 'fc_output')
regressionLayer('Name', 'output')
];
model = layerGraph(layers);
fprintf('已创建CNN-LSTM混合模型\n');
end
六、注意力机制实现
classdef attentionLayer < nnet.layer.Layer
% 自定义注意力层
properties
% 注意力权重维度
AttentionSize
end
methods
function layer = attentionLayer(name)
% 构造函数
layer.Name = name;
layer.Description = 'Attention layer for time series';
end
function Z = predict(layer, X)
% 前向传播
% X: [sequence_length, batch_size, num_features]
[seq_len, batch_size, num_features] = size(X);
% 计算注意力分数
scores = zeros(seq_len, batch_size);
for b = 1:batch_size
for t = 1:seq_len
% 计算当前时间步的重要性
scores(t,b) = sum(X(t,b,:).^2); % 简单的重要性度量
end
end
% Softmax归一化
attention_weights = softmax(scores, 1);
% 加权求和
Z = zeros(1, batch_size, num_features);
for b = 1:batch_size
weighted_sum = zeros(1, 1, num_features);
for t = 1:seq_len
weighted_sum = weighted_sum + attention_weights(t,b) * X(t,b,:);
end
Z(1,b,:) = weighted_sum;
end
end
end
end
参考代码 可用于光伏出力预测，引入气象影响因子 www.youwenfan.com/contentcst/160678.html
七、模型评估模块
function metrics = calculate_metrics(y_true, y_pred)
% 计算预测性能指标
% 确保向量形状
y_true = y_true(:);
y_pred = y_pred(:);
% 移除NaN值
valid_idx = ~isnan(y_true) & ~isnan(y_pred);
y_true = y_true(valid_idx);
y_pred = y_pred(valid_idx);
if isempty(y_true) || isempty(y_pred)
error('有效数据为空');
end
n = length(y_true);
% 1. RMSE（均方根误差）
rmse = sqrt(mean((y_true - y_pred).^2));
% 2. MAE（平均绝对误差）
mae = mean(abs(y_true - y_pred));
% 3. MAPE（平均绝对百分比误差）
% 避免除以零
non_zero_idx = y_true > 0;
if any(non_zero_idx)
mape = mean(abs((y_true(non_zero_idx) - y_pred(non_zero_idx)) ./ y_true(non_zero_idx)));
else
mape = NaN;
end
% 4. R²（决定系数）
ss_res = sum((y_true - y_pred).^2);
ss_tot = sum((y_true - mean(y_true)).^2);
r2 = 1 - (ss_res / ss_tot);
% 5. nRMSE（标准化RMSE）
y_range = max(y_true) - min(y_true);
if y_range > 0
nrmse = rmse / y_range;
else
nrmse = NaN;
end
% 6. 相关系数
corr_coef = corr(y_true, y_pred);
% 7. 对称MAPE
smape = 2 * mean(abs(y_true - y_pred) ./ (abs(y_true) + abs(y_pred))) * 100;
metrics = struct(...
'rmse', rmse, ...
'mae', mae, ...
'mape', mape, ...
'r2', r2, ...
'nrmse', nrmse, ...
'corr', corr_coef, ...
'smape', smape);
end
八、可视化模块
function visualize_results(y_true, y_pred, timestamp, metrics, params)
% 可视化预测结果
% 1. 整体预测对比图
figure('Position', [100, 100, 1200, 800]);
subplot(2,2,1);
n_show = min(500, length(y_true));
idx_show = 1:n_show;
plot(timestamp(idx_show), y_true(idx_show), 'b-', 'LineWidth', 1.5, 'DisplayName', '真实值');
hold on;
plot(timestamp(idx_show), y_pred(idx_show), 'r--', 'LineWidth', 1.5, 'DisplayName', '预测值');
xlabel('时间');
ylabel('光伏出力 (kW)');
title('光伏出力预测对比');
legend('Location', 'best');
grid on;
xlim([timestamp(1), timestamp(n_show)]);
% 2. 误差分布图
subplot(2,2,2);
errors = y_pred - y_true;
histogram(errors, 50, 'FaceColor', [0.2, 0.6, 0.8]);
xlabel('预测误差 (kW)');
ylabel('频数');
title('预测误差分布');
grid on;
% 添加统计信息
mean_err = mean(errors);
std_err = std(errors);
text(0.7, 0.9, sprintf('均值: %.3f\n标准差: %.3f', mean_err, std_err), ...
'Units', 'normalized', 'FontSize', 10);
% 3. 散点图
subplot(2,2,3);
scatter(y_true, y_pred, 20, 'filled', 'MarkerFaceAlpha', 0.3);
hold on;
plot([min(y_true), max(y_true)], [min(y_true), max(y_true)], 'r--', 'LineWidth', 2);
xlabel('真实值 (kW)');
ylabel('预测值 (kW)');
title('真实值 vs 预测值');
grid on;
axis equal;
% 添加R²值
text(0.1, 0.9, sprintf('R² = %.4f', metrics.r2), ...
'Units', 'normalized', 'FontSize', 12, 'FontWeight', 'bold');
% 4. 预测误差随时间变化
subplot(2,2,4);
plot(timestamp, abs(errors), 'b-', 'LineWidth', 1);
xlabel('时间');
ylabel('绝对误差 (kW)');
title('预测误差时间序列');
grid on;
% 5. 指标对比表
figure('Position', [300, 300, 600, 400]);
metrics_table = table(...
metrics.rmse, metrics.mae, metrics.mape*100, metrics.r2, metrics.nrmse*100, ...
'VariableNames', {'RMSE', 'MAE', 'MAPE(%)', 'R2', 'nRMSE(%)'});
uitable('Data', table2cell(metrics_table), ...
'ColumnName', metrics_table.Properties.VariableNames, ...
'RowName', {'数值'}, ...
'Position', [20, 20, 560, 360]);
% 6. 不同时间尺度的预测性能
figure('Position', [500, 100, 1000, 600]);
% 按小时分析
hours = hour(timestamp);
hourly_mae = zeros(24,1);
hourly_rmse = zeros(24,1);
for h = 0:23
hour_idx = hours == h;
if sum(hour_idx) > 0
hourly_mae(h+1) = mean(abs(y_true(hour_idx) - y_pred(hour_idx)));
hourly_rmse(h+1) = sqrt(mean((y_true(hour_idx) - y_pred(hour_idx)).^2));
end
end
subplot(2,1,1);
bar(0:23, hourly_mae, 'FaceColor', [0.2, 0.4, 0.8]);
xlabel('小时');
ylabel('MAE (kW)');
title('按小时预测性能');
grid on;
xlim([-0.5, 23.5]);
subplot(2,1,2);
bar(0:23, hourly_rmse, 'FaceColor', [0.8, 0.4, 0.2]);
xlabel('小时');
ylabel('RMSE (kW)');
grid on;
xlim([-0.5, 23.5]);
end
九、特征重要性分析
function importance = analyze_feature_importance(net, X_test, Y_test, features, params)
% 分析特征重要性
n_features = size(features, 2);
n_samples = size(X_test, 1);
% 获取基准预测
Y_base = predict(net, X_test, 'MiniBatchSize', params.batch_size);
base_error = mean(abs(Y_test(:) - Y_base(:)));
% 计算特征重要性
importance = zeros(n_features, 1);
for i = 1:n_features
% 打乱第i个特征
X_permuted = X_test;
% 对批次中的每个样本打乱特征
for b = 1:params.batch_size:n_samples
batch_end = min(b+params.batch_size-1, n_samples);
batch_size = batch_end - b + 1;
% 随机打乱特征序列
for s = 1:batch_size
seq_len = size(X_test, 2);
perm_idx = randperm(seq_len);
X_permuted(b+s-1,:,i) = X_test(b+s-1,perm_idx,i);
end
end
% 用打乱后的数据进行预测
Y_permuted = predict(net, X_permuted, 'MiniBatchSize', params.batch_size);
perm_error = mean(abs(Y_test(:) - Y_permuted(:)));
% 特征重要性 = 误差增加比例
importance(i) = (perm_error - base_error) / base_error;
end
% 特征名称
feature_names = {
'PV_当前', 'PV_滞后1h', 'PV_滞后2h', 'PV_滞后24h', 'PV_滞后168h', ...
'温度', '湿度', '风速', '辐照度', '云量', ...
'sin(小时)', 'cos(小时)', '温度×辐照度', '辐照度×(1-云量)', ...
'24h均值', '24h标准差', '晴天指数', '降雨影响', ...
'sin(年周期)', 'cos(年周期)', 'sin(周)', 'cos(周)', ...
'体感温度', '能见度指数'
};
% 创建结果结构
importance = struct(...
'scores', importance, ...
'names', {feature_names(1:n_features)}, ...
'sorted_idx', [], ...
'sorted_scores', []);
% 排序
[importance.sorted_scores, importance.sorted_idx] = sort(importance.scores, 'descend');
end
function plot_feature_importance(importance)
% 绘制特征重要性图
figure('Position', [100, 100, 1000, 600]);
n_show = min(15, length(importance.scores));
subplot(1,2,1);
barh(importance.sorted_scores(1:n_show), 'FaceColor', [0.2, 0.6, 0.8]);
set(gca, 'YTick', 1:n_show, 'YTickLabel', importance.names(importance.sorted_idx(1:n_show)));
xlabel('重要性分数');
title('Top 15 重要特征');
grid on;
subplot(1,2,2);
% 特征相关性热图
feature_names_short = cellfun(@(x) x(1:min(10,length(x))), importance.names, 'UniformOutput', false);
imagesc(importance.scores);
colorbar;
set(gca, 'XTick', 1:length(importance.scores), 'XTickLabel', feature_names_short, ...
'XTickLabelRotation', 45);
title('特征重要性热图');
xlabel('特征');
ylabel('重要性');
end
十、未来预测与敏感性分析
function future_pred = predict_future(net, features, norm_params, params)
% 进行未来预测
% 使用最近的数据
last_idx = size(features, 1);
start_idx = last_idx - params.lookback_window + 1;
% 获取输入序列
input_seq = features(start_idx:last_idx, :);
input_seq = reshape(input_seq', [1, params.lookback_window, size(features,2)]);
% 进行预测
future_pred_norm = predict(net, input_seq);
% 反标准化
future_pred = future_pred_norm * norm_params.pv_std + norm_params.pv_mean;
% 创建预测时间
last_time = datetime('now');
pred_times = last_time + hours(1:params.pred_horizon);
future_pred = struct(...
'values', future_pred, ...
'times', pred_times, ...
'confidence_interval', calculate_confidence_interval(future_pred_norm, norm_params));
end
function ci = calculate_confidence_interval(pred_norm, norm_params, alpha)
% 计算置信区间
if nargin < 3
alpha = 0.05; % 95%置信区间
end
% 简化的置信区间计算
std_multiplier = 1.96; % 95%置信区间对应的Z值
% 预测的标准差（基于历史误差）
pred_std = 0.1 * std(pred_norm) * norm_params.pv_std;
ci = struct(...
'lower', pred_norm * norm_params.pv_std - std_multiplier * pred_std, ...
'upper', pred_norm * norm_params.pv_std + std_multiplier * pred_std, ...
'mean', pred_norm * norm_params.pv_std);
end
function sensitivity = sensitivity_analysis(net, features, norm_params, params)
% 敏感性分析
% 测试不同气象条件的影响
test_cases = {
'晴天', 1.0, 0.1, 2.0, 800; % 晴天
'多云', 0.5, 0.5, 4.0, 400; % 多云
'阴天', 0.3, 0.8, 3.0, 200; % 阴天
'雨天', 0.2, 1.0, 6.0, 50; % 雨天
};
n_cases = size(test_cases, 1);
sensitivity = cell(n_cases, 1);
for i = 1:n_cases
case_name = test_cases{i,1};
% 修改气象因子
features_modified = features;
last_idx = size(features, 1);
% 假设气象因子在第6-10个特征
features_modified(last_idx, 6) = 25; % 温度
features_modified(last_idx, 7) = test_cases{i,2} * 100; % 湿度比例
features_modified(last_idx, 8) = test_cases{i,3}; % 风速
features_modified(last_idx, 9) = test_cases{i,5}; % 辐照度
features_modified(last_idx, 10) = test_cases{i,4}; % 云量
% 标准化
features_modified_norm = (features_modified - norm_params.mean) ./ norm_params.std;
% 预测
input_seq = features_modified_norm(last_idx-params.lookback_window+1:last_idx, :);
input_seq = reshape(input_seq', [1, params.lookback_window, size(features,2)]);
pred_norm = predict(net, input_seq);
pred = pred_norm * norm_params.pv_std + norm_params.pv_mean;
sensitivity{i} = struct(...
'case', case_name, ...
'conditions', test_cases(i,2:end), ...
'prediction', pred, ...
'total_energy', sum(pred));
end
end
function plot_sensitivity_analysis(sensitivity)
% 绘制敏感性分析结果
n_cases = length(sensitivity);
hours = 1:24;
figure('Position', [100, 100, 1000, 800]);
% 1. 预测曲线对比
subplot(2,2,1);
colors = lines(n_cases);
legend_entries = cell(n_cases, 1);
for i = 1:n_cases
plot(hours, sensitivity{i}.prediction, 'Color', colors(i,:), 'LineWidth', 2);
hold on;
legend_entries{i} = sensitivity{i}.case;
end
xlabel('预测小时');
ylabel('光伏出力 (kW)');
title('不同天气条件下的预测');
legend(legend_entries, 'Location', 'best');
grid on;
% 2. 总能量对比
subplot(2,2,2);
total_energy = zeros(n_cases, 1);
for i = 1:n_cases
total_energy(i) = sensitivity{i}.total_energy;
end
bar(total_energy, 'FaceColor', [0.2, 0.6, 0.8]);
set(gca, 'XTick', 1:n_cases, 'XTickLabel', legend_entries);
ylabel('总发电量 (kWh)');
title('24小时总发电量对比');
grid on;
% 3. 辐照度影响
subplot(2,2,3);
irradiance = zeros(n_cases, 1);
max_power = zeros(n_cases, 1);
for i = 1:n_cases
irradiance(i) = sensitivity{i}.conditions{4}; % 辐照度
max_power(i) = max(sensitivity{i}.prediction);
end
scatter(irradiance, max_power, 100, 'filled');
xlabel('辐照度 (W/m²)');
ylabel('最大出力 (kW)');
title('辐照度对最大出力的影响');
grid on;
% 添加拟合线
p = polyfit(irradiance, max_power, 1);
hold on;
plot(irradiance, polyval(p, irradiance), 'r--', 'LineWidth', 2);
% 4. 云量影响
subplot(2,2,4);
cloud_cover = zeros(n_cases, 1);
efficiency = zeros(n_cases, 1);
for i = 1:n_cases
cloud_cover(i) = sensitivity{i}.conditions{3}; % 云量
efficiency(i) = mean(sensitivity{i}.prediction) / max_power(i);
end
scatter(cloud_cover, efficiency*100, 100, 'filled');
xlabel('云量');
ylabel('平均效率 (%)');
title('云量对发电效率的影响');
grid on;
end
十一、模型部署与API接口
function pv_prediction_api()
% 光伏预测API接口
% 加载模型
persistent model norm_params;
if isempty(model)
load('pv_prediction_model.mat', 'net', 'norm_params', 'params');
model.net = net;
model.params = params;
fprintf('模型加载成功\n');
end
% 创建Web服务
port = 8080;
web_service = 'http://localhost:8080/predict';
% 启动Web服务
start_web_service(port);
fprintf('光伏预测API已启动，端口: %d\n', port);
fprintf('API端点: %s\n', web_service);
end
function result = predict_pv(weather_data, historical_pv)
% 预测函数（供API调用）
% 加载模型
persistent model;
if isempty(model)
load('pv_prediction_model.mat', 'net', 'norm_params', 'params');
model.net = net;
model.norm_params = norm_params;
model.params = params;
end
% 准备输入数据
features = prepare_realtime_features(weather_data, historical_pv, model.params);
% 标准化
features_norm = (features - model.norm_params.mean) ./ model.norm_params.std;
% 重塑为模型输入格式
input_seq = reshape(features_norm', [1, model.params.lookback_window, size(features,2)]);
% 预测
pred_norm = predict(model.net, input_seq);
prediction = pred_norm * model.norm_params.pv_std + model.norm_params.pv_mean;
% 返回结果
result = struct(...
'prediction', prediction, ...
'timestamp', datetime('now') + hours(1:model.params.pred_horizon), ...
'confidence', calculate_confidence(pred_norm, model.norm_params));
end
十二、使用指南
12.1 快速开始
-
运行主程序：
main_pv_prediction -
使用真实数据：
将真实数据保存为real_pv_data.mat，包含：pv_data: 光伏出力时间序列weather_data: 气象数据结构体timestamp: 时间戳
-
自定义参数：
params.pred_horizon = 48; % 预测48小时 params.lookback_window = 336; % 使用14天历史数据 params.model_type = 'cnn_lstm'; % 使用CNN-LSTM模型
12.2 关键参数说明
| 参数 | 说明 | 推荐值 |
|---|---|---|
pred_horizon | 预测步长 | 24(短期), 168(周预测) |
lookback_window | 历史窗口 | 168(7天), 672(28天) |
model_type | 模型类型 | lstm_attention(最优) |
batch_size | 批大小 | 32(小数据), 128(大数据) |
max_epochs | 最大轮数 | 100-200 |
learning_rate | 学习率 | 0.001-0.0001 |
12.3 性能优化建议
- 数据质量：确保数据完整，处理缺失值和异常值
- 特征工程：添加气象交互特征和时间特征
- 模型选择：
- 短期预测：LSTM+Attention
- 长期预测：CNN-LSTM
- 实时预测：LightGBM/XGBoost
- 超参数调优：使用贝叶斯优化搜索最佳参数
十三、预期结果
运行本系统将得到：
- 高精度预测：MAPE通常<5%（晴天），<10%（多变天气）
- 特征分析：识别关键气象因子的影响
- 可视化报告：多维度展示预测性能
- 未来预测：未来24小时光伏出力曲线
- 敏感性分析：不同天气条件下的预测对比
- 部署接口：REST API用于实时预测
十四、扩展功能
- 多站点预测：扩展支持多个光伏电站
- 概率预测：输出预测分布而非点估计
- 极端天气预警：识别异常天气对光伏的影响
- 经济性分析：结合电价预测进行收益分析
- 迁移学习：将模型迁移到新电站
- 边缘计算：部署到嵌入式设备实时预测
扫一扫
604
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
