Matlab实现综合能源系统两阶段随机优化

1. 项目概述：能源系统优化的新挑战

在能源转型的大背景下，综合能源系统(IES)正成为提高能源利用效率的关键解决方案。这个项目聚焦于一个核心难题：如何在源(发电侧)和荷(用电侧)双重不确定性的情况下，实现综合能源生产单元的最优运行调度与容量配置。

我从事能源系统优化研究已有八年，处理过数十个类似项目。这个课题的特殊之处在于它同时考虑了运行调度(短期决策)和容量配置(长期规划)两个时间尺度的优化问题，并且直面可再生能源出力波动和负荷需求变化这两类不确定性因素。用Matlab实现这类复杂优化问题，需要巧妙地将数学建模、算法设计和工程实践相结合。

2. 核心问题拆解

2.1 源荷不确定性的数学表征

处理不确定性的第一步是建立合适的数学模型。在电力系统中，风电、光伏等可再生能源的出力波动通常采用概率分布或场景生成方法来描述：

matlab复制% 风速Weibull分布参数拟合示例
wind_data = [3.2, 4.1, 5.6, ...]; % 实际风速数据
parmhat = wblfit(wind_data); % 参数估计
k = parmhat(1); % 形状参数
lambda = parmhat(2); % 尺度参数

负荷侧的不确定性则常表现为时间序列的随机波动。我推荐采用拉丁超立方抽样(LHS)生成典型场景：

matlab复制% 拉丁超立方抽样示例
num_samples = 100; % 场景数量
load_mean = [50, 60, 70, ...]; % 各时段均值
load_std = [5, 6, 7, ...]; % 各时段标准差
samples = lhsnorm(load_mean, diag(load_std.^2), num_samples);

2.2 两阶段优化框架设计

本项目采用两阶段随机规划框架：

1. 第一阶段：容量配置决策(长期投资)
2. 第二阶段：运行调度决策(短期调整)

在Matlab中实现时，我习惯先构建主问题(容量配置)，再通过回调函数处理子问题(运行调度)：

matlab复制function [total_cost, capacity] = main_problem()
    % 初始化容量决策变量
    capacity = optimvar('capacity', 'LowerBound', 0);
    
    % 调用子问题计算期望运行成本
    [operational_cost] = sub_problem(capacity);
    
    % 总成本=投资成本+运行成本
    total_cost = capital_cost(capacity) + operational_cost;
    
    % 求解优化问题
    prob = optimproblem('Objective', total_cost);
    [sol, fval] = solve(prob);
end

3. 关键技术实现细节

3.1 随机场景缩减技术

当场景数量过多时，计算量会急剧增加。我通常采用场景缩减技术保留最具代表性的场景：

matlab复制function [reduced_scenarios, weights] = scenario_reduction(original_scenarios, target_num)
    % 使用k-medoids聚类
    [idx, C] = kmedoids(original_scenarios, target_num);
    
    % 计算每个簇的权重
    weights = zeros(target_num, 1);
    for i = 1:target_num
        weights(i) = sum(idx == i) / length(idx);
    end
    
    reduced_scenarios = C;
end

提示：场景缩减后务必检查概率权重之和是否为1，这是很多初学者容易忽略的细节。

3.2 混合整数规划建模

设备启停决策需要引入二元变量，形成混合整数线性规划(MILP)问题。Matlab的intlinprog函数对此有良好支持：

matlab复制% 设备启停决策示例
units = 5; % 机组数量
hours = 24; % 时间跨度

% 定义变量：x(启停状态), y(出力)
x = optimvar('x', units, hours, 'Type', 'integer', 'LowerBound', 0, 'UpperBound', 1);
y = optimvar('y', units, hours, 'LowerBound', 0);

% 添加最小启停时间约束
for u = 1:units
    for t = 2:hours
        constraints = [constraints, 
            x(u,t) - x(u,t-1) <= x(u,min(t+min_up(u)-1,hours))];
    end
end

4. 完整实现流程

4.1 数据准备阶段

1. 收集历史源荷数据
2. 建立不确定性模型
3. 生成典型场景集
4. 场景缩减处理

matlab复制% 数据预处理完整流程
raw_data = readtable('energy_data.csv');
[wind_model, pv_model] = fit_distribution(raw_data);
scenarios = generate_scenarios(wind_model, pv_model, 1000);
[reduced_scenarios, weights] = scenario_reduction(scenarios, 10);

4.2 模型构建阶段

1. 定义决策变量
2. 建立目标函数
3. 添加系统约束
4. 设置求解器参数

matlab复制% 完整模型构建示例
prob = optimproblem;

% 决策变量
capacity = optimvar('capacity', 'LowerBound', 0);
dispatch = optimvar('dispatch', [units, hours, size(scenarios,3)], 'LowerBound', 0);

% 目标函数
prob.Objective = capital_cost(capacity) + sum(operational_cost(dispatch).*weights);

% 添加约束
prob.Constraints.power_balance = power_balance_constraint(dispatch, scenarios);
prob.Constraints.ramping = ramping_constraint(dispatch);

% 求解器配置
options = optimoptions('intlinprog', 'Display', 'iter', 'Heuristics', 'advanced');

4.3 结果分析与可视化

1. 求解优化问题
2. 提取关键结果
3. 生成分析图表

matlab复制% 结果分析示例
[sol, fval] = solve(prob, 'Options', options);

% 绘制容量配置结果
figure;
bar(sol.capacity);
title('Optimal Capacity Allocation');
xlabel('Device Type');
ylabel('Capacity (MW)');

% 绘制典型场景下的调度计划
figure;
for s = 1:3 % 展示前三个场景
    subplot(3,1,s);
    plot(squeeze(sol.dispatch(:,:,s))');
    title(['Dispatch Plan - Scenario ' num2str(s)]);
end

5. 实战经验与避坑指南

5.1 计算效率优化技巧

1. 并行计算加速：利用parfor并行处理不同场景

matlab复制% 启用并行池
if isempty(gcp('nocreate'))
    parpool('local', 4); % 使用4个核心
end

% 并行计算各场景结果
parfor s = 1:num_scenarios
    results(s) = solve_scenario(scenarios(:,:,s));
end

2. 热启动策略：用上一个场景的解作为初始点

matlab复制x0 = []; % 初始为空
for s = 1:num_scenarios
    [sol, fval] = solve(prob, 'Options', options, 'InitialPoint', x0);
    x0 = sol; % 更新初始点
end

5.2 常见问题排查

1. 求解器无法收敛：

   • 检查约束条件的可行性
   • 放宽整数变量的容差(IntTol)
   • 尝试不同的启发式策略

2. 内存不足错误：

   • 采用稀疏矩阵存储
   • 分块处理大规模场景
   • 使用问题分解算法

3. 结果不合理的可能原因：

   • 目标函数权重设置不当
   • 场景概率未归一化
   • 忽略了重要物理约束

6. 扩展应用与进阶方向

基于这个基础框架，可以考虑以下扩展：

1. 考虑碳排放约束：在目标函数中加入碳成本项

matlab复制prob.Objective = ... + carbon_price*sum(emission(dispatch).*weights);

2. 引入鲁棒优化：采用区间不确定性模型

matlab复制% 定义不确定参数范围
wind_uncertainty = interval(0.8*wind_nominal, 1.2*wind_nominal);

3. 结合机器学习：用LSTM预测源荷不确定性

matlab复制% 训练预测模型
net = trainLSTMNetwork(trainingData, responses);
prediction = predict(net, newData);

在实际项目中，我发现将随机优化与鲁棒优化结合，既能考虑概率信息，又能防范极端情况，是当前比较前沿的解决方案。