MATLAB考虑源荷不确定性的电力系统优化
一、源荷不确定性建模方法综述
在电力系统优化中，源（可再生能源出力）和荷（负荷需求）的不确定性是核心挑战。以下是主要建模方法：
| 方法 | 原理 | 优点 | 缺点 | 适用场景 |
|---|---|---|---|---|
| 确定性优化 | 使用预测均值 | 简单，计算快 | 忽略不确定性，结果保守 | 不确定性小，要求不高的场景 |
| 随机规划 | 基于概率分布，多场景 | 考虑多种可能，期望最优 | 场景数多时计算量大 | 可再生能源高渗透 |
| 鲁棒优化 | 最坏情况下最优 | 保证可行，抗干扰强 | 结果保守 | 安全要求高的场景 |
| 机会约束规划 | 允许小概率违反约束 | 平衡经济性与可靠性 | 求解复杂 | 可接受一定风险的场景 |
| 分布鲁棒优化 | 结合随机和鲁棒 | 平衡保守性和经济性 | 建模求解复杂 | 不确定性信息有限 |
本文将重点介绍基于场景的随机规划和鲁棒优化两种主流方法。
二、场景生成与削减技术
2.1 不确定性源建模
%% 1. 风电、光伏、负荷不确定性建模
clear; clc; close all;
%% 1.1 参数设置
n_hours = 24; % 24小时调度
n_scenarios = 100; % 初始场景数
n_clusters = 10; % 削减后场景数
time = 1:n_hours;
%% 1.2 风电出力不确定性建模
% 均值曲线 (典型日)
wind_mean = 50 + 20 * sin(2*pi*(time-6)/24); % MW
% 标准差 (白天波动大)
wind_std = 0.2 * wind_mean; % 20%的波动
% 生成场景 (正态分布)
wind_scenarios = zeros(n_scenarios, n_hours);
for h = 1:n_hours
wind_scenarios(:, h) = wind_mean(h) + wind_std(h) * randn(n_scenarios, 1);
% 截断在[0, 最大容量]之间
wind_scenarios(:, h) = max(0, min(wind_scenarios(:, h), 100));
end
%% 1.3 光伏出力不确定性建模
pv_mean = 60 * sin(pi*(time-6)/12); % 白天出力
pv_mean(pv_mean < 0) = 0;
pv_std = 0.15 * pv_mean;
pv_std(isnan(pv_std)) = 0.1; % 处理夜晚0出力的标准差
pv_scenarios = zeros(n_scenarios, n_hours);
for h = 1:n_hours
if pv_mean(h) > 0
pv_scenarios(:, h) = pv_mean(h) + pv_std(h) * randn(n_scenarios, 1);
pv_scenarios(:, h) = max(0, min(pv_scenarios(:, h), 80));
end
end
%% 1.4 负荷不确定性建模
load_mean = 150 + 50 * sin(2*pi*(time-14)/24); % MW
load_std = 0.1 * load_mean;
load_scenarios = zeros(n_scenarios, n_hours);
for h = 1:n_hours
load_scenarios(:, h) = load_mean(h) + load_std(h) * randn(n_scenarios, 1);
load_scenarios(:, h) = max(100, load_scenarios(:, h)); % 最小负荷
end
%% 1.5 场景概率 (初始等概率)
prob_init = ones(n_scenarios, 1) / n_scenarios;
%% 2. 场景削减 (基于K-means聚类)
function [scenarios_red, prob_red, idx] = scenario_reduction(original_scenarios, original_prob, n_clusters)
% 输入:
% original_scenarios: n_scenarios × n_hours 矩阵
% original_prob: n_scenarios × 1 概率向量
% n_clusters: 削减后场景数
% 输出:
% scenarios_red: n_clusters × n_hours 削减后场景
% prob_red: n_clusters × 1 削减后概率
% idx: 每个原始场景所属的聚类
[n_scenarios, n_hours] = size(original_scenarios);
% 使用K-means聚类
[idx, C] = kmeans(original_scenarios, n_clusters, ...
'MaxIter', 1000, 'Replicates', 10);
scenarios_red = C;
% 计算每个聚类的概率
prob_red = zeros(n_clusters, 1);
for k = 1:n_clusters
cluster_members = (idx == k);
prob_red(k) = sum(original_prob(cluster_members));
end
% 可视化聚类结果
figure;
for k = 1:min(4, n_clusters)
subplot(2,2,k);
cluster_data = original_scenarios(idx == k, :);
plot(time, cluster_data', 'Color', [0.8, 0.8, 0.8], 'LineWidth', 0.5);
hold on;
plot(time, scenarios_red(k, :), 'r-', 'LineWidth', 2);
title(sprintf('Cluster %d (Prob=%.3f)', k, prob_red(k)));
xlabel('Hour');
ylabel('Value');
grid on;
end
end
% 应用场景削减
[wind_red, prob_wind, idx_wind] = scenario_reduction(wind_scenarios, prob_init, n_clusters);
[pv_red, prob_pv, idx_pv] = scenario_reduction(pv_scenarios, prob_init, n_clusters);
[load_red, prob_load, idx_load] = scenario_reduction(load_scenarios, prob_init, n_clusters);
%% 3. 联合场景生成 (考虑相关性)
function [joint_scenarios, joint_prob] = generate_joint_scenarios(...
wind_scenarios, prob_wind, pv_scenarios, prob_pv, load_scenarios, prob_load)
n_wind = size(wind_scenarios, 1);
n_pv = size(pv_scenarios, 1);
n_load = size(load_scenarios, 1);
% 假设独立，生成所有组合
n_joint = n_wind * n_pv * n_load;
joint_scenarios = zeros(n_joint, 3, 24); % 维度: 场景 × 3个变量 × 24小时
joint_prob = zeros(n_joint, 1);
idx = 1;
for i = 1:n_wind
for j = 1:n_pv
for k = 1:n_load
joint_scenarios(idx, 1, :) = wind_scenarios(i, :);
joint_scenarios(idx, 2, :) = pv_scenarios(j, :);
joint_scenarios(idx, 3, :) = load_scenarios(k, :);
joint_prob(idx) = prob_wind(i) * prob_pv(j) * prob_load(k);
idx = idx + 1;
end
end
end
% 进一步削减联合场景
% 将三维数据展平为二维
joint_scenarios_flat = reshape(joint_scenarios, n_joint, []);
[joint_scenarios_red, joint_prob_red] = scenario_reduction(...
joint_scenarios_flat, joint_prob, 20);
% 恢复为三维
joint_scenarios = reshape(joint_scenarios_red, [], 3, 24);
joint_prob = joint_prob_red;
end
% 生成联合场景
[joint_scenarios, joint_prob] = generate_joint_scenarios(...
wind_red, prob_wind, pv_red, prob_pv, load_red, prob_load);
n_joint = size(joint_scenarios, 1);
fprintf('削减后联合场景数: %d\n', n_joint);
三、基于场景的随机规划模型
3.1 两阶段随机规划模型
%% 两阶段随机规划：考虑源荷不确定性的电力系统调度
clear; clc;
%% 1. 系统参数
n_hours = 24;
n_units = 3; % 常规机组数
n_wind = 1; % 风电场数
n_pv = 1; % 光伏电站数
n_scenarios = size(joint_scenarios, 1); % 场景数
%% 2. 机组参数
% 机组1: 基荷机组
units(1).P_min = 50; % MW
units(1).P_max = 200; % MW
units(1).a = 0.01; % 成本系数 a ($/MW^2)
units(1).b = 20; % 成本系数 b ($/MW)
units(1).c = 100; % 成本系数 c ($)
units(1).ramp_up = 100; % MW/h
units(1).ramp_down = 100; % MW/h
units(1).min_up = 4; % 最小开机时间 (h)
units(1).min_down = 4; % 最小停机时间 (h)
units(1).startup_cost = 500; % 启动成本 ($)
units(1).shutdown_cost = 100; % 停机成本 ($)
% 机组2: 腰荷机组
units(2).P_min = 20;
units(2).P_max = 100;
units(2).a = 0.02;
units(2).b = 30;
units(2).c = 50;
units(2).ramp_up = 80;
units(2).ramp_down = 80;
units(2).min_up = 2;
units(2).min_down = 2;
units(2).startup_cost = 300;
units(2).shutdown_cost = 50;
% 机组3: 峰荷机组
units(3).P_min = 0;
units(3).P_max = 50;
units(3).a = 0.05;
units(3).b = 50;
units(3).c = 0;
units(3).ramp_up = 50;
units(3).ramp_down = 50;
units(3).min_up = 1;
units(3).min_down = 1;
units(3).startup_cost = 200;
units(3).shutdown_cost = 0;
%% 3. 可再生能源参数
% 弃风/弃光惩罚
wind_penalty = 50; % $/MWh
pv_penalty = 40; % $/MWh
% 旋转备用要求
spinning_reserve = 0.1; % 10%的负荷作为旋转备用
%% 4. 构建两阶段随机规划模型
% 第一阶段变量: 机组启停状态 (二进制)
% 第二阶段变量: 机组出力、可再生能源实际出力、弃风弃光量 (连续)
% 使用YALMIP建模
addpath(genpath('yalmip')); % 添加YALMIP路径
addpath(genpath('gurobi')); % 添加Gurobi路径 (或其他求解器)
% 决策变量
u = binvar(n_units, n_hours, 'full'); % 机组启停状态
v = binvar(n_units, n_hours, 'full'); % 启动状态
w = binvar(n_units, n_hours, 'full'); % 停机状态
% 第二阶段变量 (场景相关)
P = sdpvar(n_units, n_hours, n_scenarios, 'full'); % 机组出力
P_wind_actual = sdpvar(n_hours, n_scenarios, 'full'); % 风电实际出力
P_pv_actual = sdpvar(n_hours, n_scenarios, 'full'); % 光伏实际出力
wind_curtail = sdpvar(n_hours, n_scenarios, 'full'); % 弃风量
pv_curtail = sdpvar(n_hours, n_scenarios, 'full'); % 弃光量
%% 5. 目标函数: 最小化期望总成本
total_cost = 0;
% 第一阶段成本: 启停成本
for i = 1:n_units
for t = 1:n_hours
total_cost = total_cost + units(i).startup_cost * v(i,t) + ...
units(i).shutdown_cost * w(i,t);
end
end
% 第二阶段成本期望值
for s = 1:n_scenarios
scenario_prob = joint_prob(s);
for t = 1:n_hours
% 从联合场景中获取数据
wind_available = joint_scenarios(s, 1, t);
pv_available = joint_scenarios(s, 2, t);
load_demand = joint_scenarios(s, 3, t);
% 发电成本
for i = 1:n_units
total_cost = total_cost + scenario_prob * (...
units(i).a * P(i,t,s)^2 + ...
units(i).b * P(i,t,s) + ...
units(i).c * u(i,t));
end
% 弃风弃光惩罚
total_cost = total_cost + scenario_prob * (...
wind_penalty * wind_curtail(t,s) + ...
pv_penalty * pv_curtail(t,s));
end
end
%% 6. 约束条件
constraints = [];
% 6.1 启停逻辑约束
for i = 1:n_units
for t = 1:n_hours
% 启动和停机指示器
if t > 1
constraints = [constraints, v(i,t) >= u(i,t) - u(i,t-1)];
constraints = [constraints, w(i,t) >= u(i,t-1) - u(i,t)];
else
constraints = [constraints, v(i,1) >= u(i,1) - 0]; % 初始状态假设为0
constraints = [constraints, w(i,1) >= 0 - u(i,1)];
end
constraints = [constraints, v(i,t) <= u(i,t)];
constraints = [constraints, w(i,t) <= 1 - u(i,t)];
constraints = [constraints, v(i,t) + w(i,t) <= 1];
end
end
% 6.2 最小启停时间约束
for i = 1:n_units
for t = 1:n_hours
% 最小开机时间
if t + units(i).min_up - 1 <= n_hours
constraints = [constraints, ...
sum(u(i, t:min(t+units(i).min_up-1, n_hours))) >= ...
units(i).min_up * v(i,t)];
end
% 最小停机时间
if t + units(i).min_down - 1 <= n_hours
constraints = [constraints, ...
sum(1 - u(i, t:min(t+units(i).min_down-1, n_hours))) >= ...
units(i).min_down * w(i,t)];
end
end
end
% 6.3 场景相关的约束
for s = 1:n_scenarios
for t = 1:n_hours
% 从联合场景中获取数据
wind_available = joint_scenarios(s, 1, t);
pv_available = joint_scenarios(s, 2, t);
load_demand = joint_scenarios(s, 3, t);
% 功率平衡约束
total_gen = sum(P(:,t,s)) + P_wind_actual(t,s) + P_pv_actual(t,s);
constraints = [constraints, total_gen == load_demand];
% 机组出力上下限约束
for i = 1:n_units
constraints = [constraints, ...
units(i).P_min * u(i,t) <= P(i,t,s) <= units(i).P_max * u(i,t)];
end
% 爬坡约束
if t > 1
for i = 1:n_units
constraints = [constraints, ...
P(i,t,s) - P(i,t-1,s) <= units(i).ramp_up];
constraints = [constraints, ...
P(i,t-1,s) - P(i,t,s) <= units(i).ramp_down];
end
end
% 可再生能源出力约束
constraints = [constraints, 0 <= P_wind_actual(t,s) <= wind_available];
constraints = [constraints, 0 <= P_pv_actual(t,s) <= pv_available];
% 弃风弃光定义
constraints = [constraints, ...
wind_curtail(t,s) == wind_available - P_wind_actual(t,s)];
constraints = [constraints, ...
pv_curtail(t,s) == pv_available - P_pv_actual(t,s)];
% 旋转备用约束
quick_reserve = 0;
for i = 1:n_units
quick_reserve = quick_reserve + min(units(i).ramp_up, ...
units(i).P_max * u(i,t) - P(i,t,s));
end
constraints = [constraints, quick_reserve >= spinning_reserve * load_demand];
end
end
%% 7. 求解优化问题
ops = sdpsettings('solver', 'gurobi', 'verbose', 1, 'showprogress', 1);
ops.gurobi.TimeLimit = 3600; % 1小时时间限制
ops.gurobi.MIPGap = 0.01; % 1%的MIP间隙
% 求解
tic;
diagnostics = optimize(constraints, total_cost, ops);
solve_time = toc;
if diagnostics.problem == 0
disp('=== 求解成功 ===');
fprintf('求解时间: %.2f 秒\n', solve_time);
fprintf('目标函数值: %.2f $\n', value(total_cost));
% 提取结果
u_opt = value(u);
P_opt = value(P);
P_wind_opt = value(P_wind_actual);
P_pv_opt = value(P_pv_actual);
wind_curtail_opt = value(wind_curtail);
pv_curtail_opt = value(pv_curtail);
% 计算期望值
P_expected = zeros(n_units, n_hours);
for i = 1:n_units
for t = 1:n_hours
for s = 1:n_scenarios
P_expected(i,t) = P_expected(i,t) + joint_prob(s) * P_opt(i,t,s);
end
end
end
else
disp('求解失败');
disp(yalmiperror(diagnostics.problem));
end
%% 8. 结果可视化
plot_results(u_opt, P_expected, P_wind_opt, P_pv_opt, ...
wind_curtail_opt, pv_curtail_opt, joint_scenarios, joint_prob);
3.2 结果可视化函数
function plot_results(u_opt, P_expected, P_wind_opt, P_pv_opt, ...
wind_curtail_opt, pv_curtail_opt, joint_scenarios, joint_prob)
[n_units, n_hours] = size(P_expected);
n_scenarios = size(joint_scenarios, 1);
time = 1:n_hours;
figure('Position', [100, 100, 1400, 900]);
% 子图1: 机组启停状态
subplot(3,3,1);
imagesc(time, 1:n_units, u_opt);
colormap(flipud(gray));
colorbar;
xlabel('时间 (小时)');
ylabel('机组编号');
title('机组启停状态');
set(gca, 'YTick', 1:n_units);
% 子图2: 期望出力
subplot(3,3,2);
area_plot = area(time, P_expected');
legend('机组1', '机组2', '机组3', 'Location', 'best');
xlabel('时间 (小时)');
ylabel('出力 (MW)');
title('机组期望出力');
grid on;
% 子图3: 可再生能源出力
subplot(3,3,3);
hold on;
% 计算期望的可再生能源出力
P_wind_expected = zeros(1, n_hours);
P_pv_expected = zeros(1, n_hours);
for t = 1:n_hours
for s = 1:n_scenarios
P_wind_expected(t) = P_wind_expected(t) + joint_prob(s) * P_wind_opt(t,s);
P_pv_expected(t) = P_pv_expected(t) + joint_prob(s) * P_pv_opt(t,s);
end
end
plot(time, P_wind_expected, 'b-', 'LineWidth', 2, 'DisplayName', '风电期望出力');
plot(time, P_pv_expected, 'g-', 'LineWidth', 2, 'DisplayName', '光伏期望出力');
% 计算可用可再生能源
wind_available_exp = zeros(1, n_hours);
pv_available_exp = zeros(1, n_hours);
for t = 1:n_hours
for s = 1:n_scenarios
wind_available_exp(t) = wind_available_exp(t) + joint_prob(s) * joint_scenarios(s,1,t);
pv_available_exp(t) = pv_available_exp(t) + joint_prob(s) * joint_scenarios(s,2,t);
end
end
plot(time, wind_available_exp, 'b--', 'LineWidth', 1, 'DisplayName', '风电可用');
plot(time, pv_available_exp, 'g--', 'LineWidth', 1, 'DisplayName', '光伏可用');
xlabel('时间 (小时)');
ylabel('出力 (MW)');
title('可再生能源出力');
legend('Location', 'best');
grid on;
% 子图4: 弃风弃光
subplot(3,3,4);
wind_curtail_exp = zeros(1, n_hours);
pv_curtail_exp = zeros(1, n_hours);
for t = 1:n_hours
for s = 1:n_scenarios
wind_curtail_exp(t) = wind_curtail_exp(t) + joint_prob(s) * wind_curtail_opt(t,s);
pv_curtail_exp(t) = pv_curtail_exp(t) + joint_prob(s) * pv_curtail_opt(t,s);
end
end
bar_data = [wind_curtail_exp; pv_curtail_exp]';
bar(time, bar_data, 'stacked');
xlabel('时间 (小时)');
ylabel('弃电量 (MWh)');
title('弃风弃光量');
legend('弃风', '弃光', 'Location', 'best');
grid on;
% 子图5: 成本分解
subplot(3,3,5);
% 计算各项成本
fuel_cost = 0;
startup_cost = 0;
shutdown_cost = 0;
curtail_cost = 0;
% 这里需要从优化结果中计算具体成本
% 简化显示
cost_labels = {'燃料成本', '启停成本', '弃电惩罚'};
cost_values = [70, 15, 5]; % 示例值，单位: 千美元
bar(cost_values);
set(gca, 'XTickLabel', cost_labels);
ylabel('成本 (千$)');
title('成本分解');
grid on;
% 子图6: 典型场景下的功率平衡
subplot(3,3,6);
% 选择概率最大的场景
[~, max_prob_idx] = max(joint_prob);
load_scenario = squeeze(joint_scenarios(max_prob_idx, 3, :))';
total_gen = sum(P_opt(:,:,max_prob_idx), 1) + ...
P_wind_opt(:,max_prob_idx)' + P_pv_opt(:,max_prob_idx)';
area_plot = area(time, [P_opt(:,:,max_prob_idx)', ...
P_wind_opt(:,max_prob_idx)', ...
P_pv_opt(:,max_prob_idx)']);
hold on;
plot(time, load_scenario, 'k-', 'LineWidth', 2, 'DisplayName', '负荷');
xlabel('时间 (小时)');
ylabel('功率 (MW)');
title('典型场景功率平衡');
legend('机组1', '机组2', '机组3', '风电', '光伏', '负荷', 'Location', 'best');
grid on;
% 子图7: 备用充足性
subplot(3,3,7);
reserve_available = zeros(1, n_hours);
reserve_required = zeros(1, n_hours);
for t = 1:n_hours
load_exp = 0;
for s = 1:n_scenarios
load_exp = load_exp + joint_prob(s) * joint_scenarios(s,3,t);
end
reserve_required(t) = 0.1 * load_exp; % 10%备用
% 计算可用备用
for s = 1:n_scenarios
for i = 1:n_units
reserve_available(t) = reserve_available(t) + joint_prob(s) * ...
min(units(i).ramp_up, units(i).P_max*u_opt(i,t)-P_opt(i,t,s));
end
end
end
bar(time, [reserve_available; reserve_required]', 'grouped');
xlabel('时间 (小时)');
ylabel('备用容量 (MW)');
title('旋转备用');
legend('可用备用', '要求备用', 'Location', 'best');
grid on;
% 子图8: 场景概率分布
subplot(3,3,8);
[sorted_prob, idx] = sort(joint_prob, 'descend');
bar(1:length(sorted_prob), sorted_prob);
xlabel('场景编号 (按概率排序)');
ylabel('概率');
title('场景概率分布');
grid on;
% 子图9: 不确定性带
subplot(3,3,9);
% 计算负荷的不确定性带
load_min = min(squeeze(joint_scenarios(:,3,:)), [], 1);
load_max = max(squeeze(joint_scenarios(:,3,:)), [], 1);
load_mean = mean(squeeze(joint_scenarios(:,3,:)), 1);
fill([time, fliplr(time)], [load_max, fliplr(load_min)], ...
[0.8, 0.8, 0.8], 'FaceAlpha', 0.3, 'EdgeColor', 'none');
hold on;
plot(time, load_mean, 'b-', 'LineWidth', 2);
xlabel('时间 (小时)');
ylabel('负荷 (MW)');
title('负荷不确定性带');
legend('不确定性范围', '期望值', 'Location', 'best');
grid on;
end
四、鲁棒优化模型
4.1 自适应鲁棒优化框架
%% 自适应鲁棒优化：考虑最坏情况下的源荷不确定性
clear; clc;
%% 1. 系统参数
n_hours = 24;
n_units = 3;
n_wind = 1;
n_pv = 1;
%% 2. 不确定性集合定义
% 区间不确定性模型
load_nominal = [120, 115, 110, 105, 100, 110, 130, 150, 170, 180, ...
175, 170, 165, 160, 155, 160, 170, 185, 190, 180, ...
170, 160, 140, 125]; % MW
wind_nominal = [45, 48, 50, 52, 55, 60, 65, 70, 75, 80, ...
78, 75, 72, 70, 68, 65, 62, 60, 58, 55, ...
52, 50, 48, 46]; % MW
pv_nominal = [0, 0, 0, 0, 0, 5, 15, 30, 45, 60, ...
70, 75, 80, 75, 70, 60, 45, 30, 15, 5, ...
0, 0, 0, 0]; % MW
% 不确定性偏差范围 (±百分比)
load_uncertainty = 0.15; % ±15%
wind_uncertainty = 0.30; % ±30%
pv_uncertainty = 0.25; % ±25%
% 预算不确定性集合参数
Gamma_load = 8; % 负荷不确定性预算
Gamma_wind = 6; % 风电不确定性预算
Gamma_pv = 4; % 光伏不确定性预算
%% 3. 鲁棒优化建模
% 使用YALMIP的鲁棒优化功能
yalmip('clear');
% 第一阶段变量: 机组启停
u = binvar(n_units, n_hours, 'full');
v = binvar(n_units, n_hours, 'full');
w = binvar(n_units, n_hours, 'full');
% 第二阶段变量: 机组出力 (作为仿射可调决策)
P = sdpvar(n_units, n_hours, 'full');
% 不确定变量
xi_load = sdpvar(1, n_hours, 'full'); % 负荷不确定性
xi_wind = sdpvar(1, n_hours, 'full'); % 风电不确定性
xi_pv = sdpvar(1, n_hours, 'full'); % 光伏不确定性
% 不确定性集合定义
% 标准化不确定性变量: ξ ∈ [-1, 1]
uncertainty_set = [];
for t = 1:n_hours
uncertainty_set = [uncertainty_set, -1 <= xi_load(t) <= 1];
uncertainty_set = [uncertainty_set, -1 <= xi_wind(t) <= 1];
uncertainty_set = [uncertainty_set, -1 <= xi_pv(t) <= 1];
end
% 预算约束
uncertainty_set = [uncertainty_set, ...
norm(xi_load, 1) <= Gamma_load, ...
norm(xi_wind, 1) <= Gamma_wind, ...
norm(xi_pv, 1) <= Gamma_pv];
% 实际值 = 标称值 + 不确定性
load_actual = load_nominal + load_uncertainty * load_nominal .* xi_load;
wind_actual = wind_nominal + wind_uncertainty * wind_nominal .* xi_wind;
pv_actual = pv_nominal + pv_uncertainty * pv_nominal .* xi_pv;
%% 4. 约束条件
constraints = [];
% 4.1 启停逻辑约束
for i = 1:n_units
for t = 1:n_hours
if t > 1
constraints = [constraints, v(i,t) >= u(i,t) - u(i,t-1)];
constraints = [constraints, w(i,t) >= u(i,t-1) - u(i,t)];
else
constraints = [constraints, v(i,1) >= u(i,1)];
constraints = [constraints, w(i,1) >= 0 - u(i,1)];
end
constraints = [constraints, 0 <= v(i,t) <= 1, 0 <= w(i,t) <= 1];
end
end
% 4.2 机组出力约束
for i = 1:n_units
for t = 1:n_hours
constraints = [constraints, ...
units(i).P_min * u(i,t) <= P(i,t) <= units(i).P_max * u(i,t)];
end
end
% 4.3 爬坡约束
for i = 1:n_units
for t = 2:n_hours
constraints = [constraints, ...
P(i,t) - P(i,t-1) <= units(i).ramp_up];
constraints = [constraints, ...
P(i,t-1) - P(i,t) <= units(i).ramp_down];
end
end
% 4.4 功率平衡约束 (鲁棒约束，对所有不确定性实现成立)
power_balance_constraint = [];
for t = 1:n_hours
total_gen = sum(P(:,t)) + wind_actual(t) + pv_actual(t);
power_balance_constraint = [power_balance_constraint, ...
total_gen == load_actual(t)];
end
constraints = [constraints, implies(uncertainty_set, power_balance_constraint)];
% 4.5 备用约束
reserve_constraint = [];
for t = 1:n_hours
% 计算可用备用
available_reserve = 0;
for i = 1:n_units
available_reserve = available_reserve + ...
min(units(i).ramp_up, units(i).P_max*u(i,t) - P(i,t));
end
% 考虑最坏情况下的备用需求
worst_case_load = load_nominal(t) * (1 + load_uncertainty);
reserve_requirement = 0.1 * worst_case_load;
reserve_constraint = [reserve_constraint, ...
available_reserve >= reserve_requirement];
end
constraints = [constraints, reserve_constraint];
%% 5. 目标函数
% 最小化最坏情况下的总成本
total_cost = 0;
% 启停成本
for i = 1:n_units
for t = 1:n_hours
total_cost = total_cost + units(i).startup_cost * v(i,t) + ...
units(i).shutdown_cost * w(i,t);
end
end
% 发电成本 (考虑最坏情况)
gen_cost = 0;
for i = 1:n_units
for t = 1:n_hours
gen_cost = gen_cost + units(i).a * P(i,t)^2 + ...
units(i).b * P(i,t) + units(i).c * u(i,t);
end
end
% 弃电惩罚 (在最坏情况下)
wind_curtail_penalty = 0;
pv_curtail_penalty = 0;
% 由于是鲁棒优化，需要考虑最坏情况下的弃电
% 简化处理: 假设弃电发生在可再生能源过剩时
for t = 1:n_hours
% 弃电惩罚项
wind_curtail_penalty = wind_curtail_penalty + wind_penalty * ...
max(0, wind_actual(t) - (load_actual(t) - sum(P(:,t)) - pv_actual(t)));
pv_curtail_penalty = pv_curtail_penalty + pv_penalty * ...
max(0, pv_actual(t) - (load_actual(t) - sum(P(:,t)) - wind_actual(t)));
end
total_cost = total_cost + gen_cost + wind_curtail_penalty + pv_curtail_penalty;
%% 6. 求解鲁棒优化问题
ops = sdpsettings('solver', 'gurobi', 'robust.lplp', 'duality');
ops.verbose = 1;
ops.gurobi.TimeLimit = 1800;
tic;
diagnostics = optimize([constraints, uncertainty_set], total_cost, ops);
solve_time = toc;
if diagnostics.problem == 0
disp('=== 鲁棒优化求解成功 ===');
fprintf('求解时间: %.2f 秒\n', solve_time);
fprintf('最坏情况总成本: %.2f $\n', value(total_cost));
% 提取结果
u_robust = value(u);
P_robust = value(P);
% 分析鲁棒性
analyze_robustness(u_robust, P_robust, load_nominal, wind_nominal, pv_nominal, ...
load_uncertainty, wind_uncertainty, pv_uncertainty, units);
else
disp('鲁棒优化求解失败');
disp(yalmiperror(diagnostics.problem));
end
4.2 鲁棒性分析函数
function analyze_robustness(u_opt, P_opt, load_nominal, wind_nominal, pv_nominal, ...
load_uncert, wind_uncert, pv_uncert, units)
n_hours = length(load_nominal);
n_units = size(P_opt, 1);
time = 1:n_hours;
% 生成最坏情况场景
n_monte_carlo = 1000;
feasibility_count = 0;
cost_distribution = zeros(n_monte_carlo, 1);
for mc = 1:n_monte_carlo
% 随机生成不确定性实现
xi_load = 2*rand(1, n_hours) - 1; % [-1, 1]
xi_wind = 2*rand(1, n_hours) - 1;
xi_pv = 2*rand(1, n_hours) - 1;
% 应用预算约束
if sum(abs(xi_load)) > 8
xi_load = xi_load * 8 / sum(abs(xi_load));
end
if sum(abs(xi_wind)) > 6
xi_wind = xi_wind * 6 / sum(abs(xi_wind));
end
if sum(abs(xi_pv)) > 4
xi_pv = xi_pv * 4 / sum(abs(xi_pv));
end
% 计算实际值
load_actual = load_nominal + load_uncert * load_nominal .* xi_load;
wind_actual = wind_nominal + wind_uncert * wind_nominal .* xi_wind;
pv_actual = pv_nominal + pv_uncert * pv_nominal .* xi_pv;
% 检查可行性
feasible = true;
for t = 1:n_hours
% 检查功率平衡
total_gen = sum(P_opt(:,t)) + wind_actual(t) + pv_actual(t);
if abs(total_gen - load_actual(t)) > 1e-3
feasible = false;
break;
end
% 检查机组约束
for i = 1:n_units
if P_opt(i,t) < units(i).P_min * u_opt(i,t) - 1e-3 || ...
P_opt(i,t) > units(i).P_max * u_opt(i,t) + 1e-3
feasible = false;
break;
end
end
if ~feasible, break; end
end
if feasible
feasibility_count = feasibility_count + 1;
% 计算该场景下的成本
scenario_cost = 0;
for i = 1:n_units
for t = 1:n_hours
scenario_cost = scenario_cost + ...
units(i).a * P_opt(i,t)^2 + ...
units(i).b * P_opt(i,t) + ...
units(i).c * u_opt(i,t);
end
end
cost_distribution(mc) = scenario_cost;
end
end
feasibility_rate = feasibility_count / n_monte_carlo;
fprintf('蒙特卡洛可行性检验: %.2f%%\n', feasibility_rate * 100);
% 可视化鲁棒性分析结果
figure('Position', [100, 100, 1200, 800]);
% 子图1: 成本分布
subplot(2,3,1);
histogram(cost_distribution(cost_distribution > 0), 20);
xlabel('总成本 ($)');
ylabel('频数');
title('不同场景下的成本分布');
grid on;
% 子图2: 最坏情况与标称情况对比
subplot(2,3,2);
% 计算标称情况下的出力
P_nominal = zeros(n_units, n_hours);
for t = 1:n_hours
% 简单分配: 按机组顺序满足净负荷
net_load = load_nominal(t) - wind_nominal(t) - pv_nominal(t);
remaining = net_load;
for i = 1:n_units
if remaining > 0
P_nominal(i,t) = min(remaining, units(i).P_max * u_opt(i,t));
remaining = remaining - P_nominal(i,t);
end
end
end
plot(time, sum(P_opt,1), 'r-', 'LineWidth', 2, 'DisplayName', '鲁棒调度');
hold on;
plot(time, sum(P_nominal,1), 'b--', 'LineWidth', 2, 'DisplayName', '标称调度');
xlabel('时间 (小时)');
ylabel('常规机组总出力 (MW)');
title('鲁棒 vs 标称调度');
legend('Location', 'best');
grid on;
% 子图3: 备用充足性
subplot(2,3,3);
reserve_available = zeros(1, n_hours);
worst_case_reserve_needed = zeros(1, n_hours);
for t = 1:n_hours
% 可用备用
for i = 1:n_units
reserve_available(t) = reserve_available(t) + ...
min(units(i).ramp_up, units(i).P_max*u_opt(i,t) - P_opt(i,t));
end
% 最坏情况备用需求
worst_load = load_nominal(t) * (1 + load_uncert);
worst_case_reserve_needed(t) = 0.1 * worst_load;
end
bar(time, [reserve_available; worst_case_reserve_needed]');
xlabel('时间 (小时)');
ylabel('备用容量 (MW)');
title('备用充足性分析');
legend('可用备用', '最坏情况需求', 'Location', 'best');
grid on;
% 子图4: 不确定性对功率平衡的影响
subplot(2,3,4);
% 生成几个极端场景
n_extreme = 5;
imbalance_scenarios = zeros(n_extreme, n_hours);
for s = 1:n_extreme
% 生成极端不确定性
xi_load_extreme = sign(randn(1, n_hours));
xi_wind_extreme = sign(randn(1, n_hours));
xi_pv_extreme = sign(randn(1, n_hours));
% 计算实际值
load_extreme = load_nominal + load_uncert * load_nominal .* xi_load_extreme;
wind_extreme = wind_nominal + wind_uncert * wind_nominal .* xi_wind_extreme;
pv_extreme = pv_nominal + pv_uncert * pv_nominal .* xi_pv_extreme;
% 计算功率不平衡
for t = 1:n_hours
total_gen = sum(P_opt(:,t)) + wind_extreme(t) + pv_extreme(t);
imbalance_scenarios(s,t) = total_gen - load_extreme(t);
end
end
plot(time, imbalance_scenarios', 'LineWidth', 1);
hold on;
plot(time, zeros(1,n_hours), 'k-', 'LineWidth', 2);
xlabel('时间 (小时)');
ylabel('功率不平衡 (MW)');
title('极端场景下的功率平衡');
grid on;
% 子图5: 成本保守性分析
subplot(2,3,5);
% 计算标称成本
nominal_cost = 0;
for i = 1:n_units
for t = 1:n_hours
nominal_cost = nominal_cost + ...
units(i).a * P_nominal(i,t)^2 + ...
units(i).b * P_nominal(i,t) + ...
units(i).c * u_opt(i,t);
end
end
% 鲁棒成本 (期望)
robust_cost_mean = mean(cost_distribution(cost_distribution > 0));
robust_cost_max = max(cost_distribution);
bar_data = [nominal_cost, robust_cost_mean, robust_cost_max];
bar(1:3, bar_data);
set(gca, 'XTickLabel', {'标称成本', '鲁棒期望成本', '最坏情况成本'});
ylabel('成本 ($)');
title('成本保守性分析');
grid on;
% 添加数值标签
for i = 1:3
text(i, bar_data(i), sprintf('%.0f', bar_data(i)), ...
'HorizontalAlignment', 'center', 'VerticalAlignment', 'bottom');
end
% 子图6: 可行性热图
subplot(2,3,6);
% 测试不同不确定性水平下的可行性
gamma_test = 0:2:10;
feasibility_rates = zeros(length(gamma_test), 1);
for g = 1:length(gamma_test)
gamma = gamma_test(g);
feasible_count = 0;
for mc = 1:100
xi_test = 2*rand(1, n_hours) - 1;
if sum(abs(xi_test)) > gamma
xi_test = xi_test * gamma / sum(abs(xi_test));
end
load_test = load_nominal + load_uncert * load_nominal .* xi_test;
% 简单可行性检查
feasible = true;
for t = 1:n_hours
if load_test(t) > sum(P_opt(:,t)) + wind_nominal(t) + pv_nominal(t) + 10
feasible = false;
break;
end
end
if feasible
feasible_count = feasible_count + 1;
end
end
feasibility_rates(g) = feasible_count / 100;
end
plot(gamma_test, feasibility_rates*100, 'b-o', 'LineWidth', 2);
xlabel('不确定性预算 \Gamma');
ylabel('可行性率 (%)');
title('不同不确定性水平下的可行性');
grid on;
ylim([0, 105]);
end
参考代码 考虑源荷不确定性-电力系统优化 www.youwenfan.com/contentcnt/160705.html
五、实际应用建议
5.1 方法选择指南
| 应用场景 | 推荐方法 | 原因 |
|---|---|---|
| 日前市场出清 | 随机规划 | 考虑多种可能，经济性较好 |
| 实时调度 | 鲁棒优化 | 快速响应，保证安全性 |
| 规划问题 | 随机规划+鲁棒优化 | 兼顾长期经济性和短期可靠性 |
| 高可再生能源渗透 | 随机规划+机会约束 | 平衡经济性与消纳 |
| 极端天气应对 | 鲁棒优化 | 抗干扰能力强 |
5.2 参数设置经验
- 场景数：50-200个场景通常足够，更多场景收益递减
- 不确定性预算Γ：通常设为不确定性变量数的20-40%
- 弃电惩罚：设为最高发电成本的1.2-1.5倍
- 备用要求：常规10-15%，高可再生能源时15-20%
5.3 计算效率提升
- 并行计算：使用
parfor并行处理场景 - 分解算法：Benders分解、场景分解
- 场景削减：保留概率大、差异大的场景
- 简化模型：线性化、分段线性化复杂函数
六、总结
本文提供了完整的MATLAB实现，涵盖：
- 不确定性建模：场景生成与削减
- 随机规划：两阶段随机优化模型
- 鲁棒优化：自适应鲁棒优化框架
- 可视化分析：全面的结果可视化
关键要点：
- 随机规划更适合追求经济性的场景
- 鲁棒优化更适合安全要求高的场景
- 实际应用常结合两种方法
- 场景质量和数量显著影响结果
后续工作方向：
- 集成机器学习进行不确定性预测
- 考虑网络约束的扩展
- 多时间尺度协调优化
- 分布式优化算法
代码已模块化设计，便于扩展和应用于实际电力系统。
