粒子群算法求解电力系统经济调度问题
使用粒子群优化算法(PSO)求解电力系统经济调度问题的完整MATLAB实现。程序考虑了发电机组的运行成本、功率平衡约束和发电机组出力限制。
% 粒子群算法求解电力系统经济调度问题
function pso_economic_dispatch()
% 清空环境
clc;
clear;
close all;
% 设置随机数种子
rng(0, 'twister');
% 发电机参数（成本系数：a, b, c；出力范围：Pmin, Pmax）
% 6台发电机参数
gen_data = [
0.001562, 7.92, 561, 150, 600; % 机组1
0.001940, 7.85, 310, 100, 400; % 机组2
0.004820, 7.97, 78, 50, 200; % 机组3
0.002110, 8.03, 152, 80, 300; % 机组4
0.003980, 7.87, 113, 70, 250; % 机组5
0.007120, 8.05, 95, 60, 180 % 机组6
];
% 系统参数
Pd = 1800; % 总负荷需求 (MW)
loss_coeff = 1e-4; % 网损系数
% PSO参数
n_particles = 50; % 粒子数量
max_iter = 200; % 最大迭代次数
w = 0.9; % 初始惯性权重
w_min = 0.4; % 最小惯性权重
c1 = 2.0; % 个体学习因子
c2 = 2.0; % 群体学习因子
n_gens = size(gen_data, 1); % 发电机组数量
% 初始化粒子群
particles = struct();
for i = 1:n_particles
% 随机初始化位置（机组出力）
particles(i).position = zeros(1, n_gens);
for j = 1:n_gens
Pmin = gen_data(j, 4);
Pmax = gen_data(j, 5);
particles(i).position(j) = Pmin + (Pmax - Pmin) * rand();
end
% 调整出力以满足负荷需求
particles(i).position = adjust_generation(particles(i).position, Pd, gen_data);
% 初始化速度
particles(i).velocity = zeros(1, n_gens);
% 计算初始成本
[cost, penalty] = calculate_cost(particles(i).position, Pd, gen_data, loss_coeff);
particles(i).cost = cost + penalty;
% 初始化个体最优
particles(i).best_position = particles(i).position;
particles(i).best_cost = particles(i).cost;
end
% 寻找全局最优粒子
[global_best_cost, idx] = min([particles.best_cost]);
global_best_position = particles(idx).best_position;
% 记录迭代过程
best_costs = zeros(max_iter, 1);
avg_costs = zeros(max_iter, 1);
convergence_curve = zeros(max_iter, 1);
best_penalties = zeros(max_iter, 1);
% PSO主循环
for iter = 1:max_iter
% 更新惯性权重（线性递减）
w = w_min + (w - w_min) * (max_iter - iter) / max_iter;
total_cost = 0;
total_penalty = 0;
for i = 1:n_particles
% 更新速度
r1 = rand(1, n_gens);
r2 = rand(1, n_gens);
cognitive = c1 * r1 .* (particles(i).best_position - particles(i).position);
social = c2 * r2 .* (global_best_position - particles(i).position);
particles(i).velocity = w * particles(i).velocity + cognitive + social;
% 限制速度范围
max_velocity = 0.1 * ([gen_data(:,5) - gen_data(:,4)])';
particles(i).velocity = min(max(particles(i).velocity, -max_velocity), max_velocity);
% 更新位置
particles(i).position = particles(i).position + particles(i).velocity;
% 确保位置在可行范围内
for j = 1:n_gens
Pmin = gen_data(j, 4);
Pmax = gen_data(j, 5);
particles(i).position(j) = max(min(particles(i).position(j), Pmax), Pmin);
end
% 调整出力以满足负荷需求
particles(i).position = adjust_generation(particles(i).position, Pd, gen_data);
% 计算新成本
[cost, penalty] = calculate_cost(particles(i).position, Pd, gen_data, loss_coeff);
particles(i).cost = cost + penalty;
% 更新个体最优
if particles(i).cost < particles(i).best_cost
particles(i).best_position = particles(i).position;
particles(i).best_cost = particles(i).cost;
end
% 更新全局最优
if particles(i).best_cost < global_best_cost
global_best_cost = particles(i).best_cost;
global_best_position = particles(i).best_position;
end
total_cost = total_cost + cost;
total_penalty = total_penalty + penalty;
end
% 记录统计信息
best_costs(iter) = min([particles.cost]);
avg_costs(iter) = total_cost / n_particles;
convergence_curve(iter) = global_best_cost;
best_penalties(iter) = min([particles.cost] - [particles.cost] + [particles.cost]); % 占位符
% 显示迭代信息
if mod(iter, 10) == 0
fprintf('迭代 %d: 最优成本 = %.2f $/h\n', iter, global_best_cost);
end
end
% 计算网损
total_gen = sum(global_best_position);
Ploss = loss_coeff * total_gen^2;
% 显示最终结果
fprintf('\n======= 经济调度最优解 =======\n');
fprintf('总负荷需求: %.2f MW\n', Pd);
fprintf('总发电量: %.2f MW\n', total_gen);
fprintf('网损: %.4f MW\n', Ploss);
fprintf('总发电成本: %.2f $/h\n\n', global_best_cost);
fprintf('机组\t出力(MW)\t成本($/h)\t运行状态\n');
fprintf('------------------------------------\n');
for i = 1:n_gens
P = global_best_position(i);
a = gen_data(i, 1);
b = gen_data(i, 2);
c = gen_data(i, 3);
cost = a*P^2 + b*P + c;
status = '运行';
if P < gen_data(i, 4) + 0.1
status = '停机';
end
fprintf('%d\t%.2f\t\t%.2f\t\t%s\n', i, P, cost, status);
end
% 可视化结果
visualize_results(gen_data, global_best_position, best_costs, avg_costs, convergence_curve, max_iter);
end
% 调整发电量以满足负荷需求
function new_position = adjust_generation(position, Pd, gen_data)
n_gens = length(position);
Pmin = gen_data(:, 4)';
Pmax = gen_data(:, 5)';
% 计算当前总发电量
total_gen = sum(position);
% 计算发电量偏差
deviation = Pd - total_gen;
% 如果偏差很小，直接返回
if abs(deviation) < 1e-3
new_position = position;
return;
end
% 按比例调整发电量（考虑机组出力限制）
if deviation > 0 % 发电不足，需要增加出力
available_headroom = Pmax - position;
total_headroom = sum(available_headroom);
if total_headroom < deviation
% 如果可用裕度不足，将所有机组设为最大出力
new_position = Pmax;
else
% 按裕度比例分配增加量
scaling_factors = available_headroom / total_headroom;
new_position = position + scaling_factors * deviation;
end
else % 发电过剩，需要减少出力
available_rampdown = position - Pmin;
total_rampdown = sum(available_rampdown);
if total_rampdown < abs(deviation)
% 如果可用减量不足，将所有机组设为最小出力
new_position = Pmin;
else
% 按减量比例分配减少量
scaling_factors = available_rampdown / total_rampdown;
new_position = position + scaling_factors * deviation;
end
end
% 确保调整后的出力在限制范围内
new_position = max(min(new_position, Pmax), Pmin);
end
% 计算总成本和罚值
function [cost, penalty] = calculate_cost(position, Pd, gen_data, loss_coeff)
n_gens = length(position);
cost = 0;
% 计算发电成本（二次成本函数）
for i = 1:n_gens
P = position(i);
a = gen_data(i, 1);
b = gen_data(i, 2);
c = gen_data(i, 3);
cost = cost + a*P^2 + b*P + c;
end
% 计算网损
total_gen = sum(position);
Ploss = loss_coeff * total_gen^2;
% 计算功率平衡偏差
power_balance = abs(total_gen - Ploss - Pd);
% 计算罚值
penalty_factor = 1000; % 罚值系数
penalty = penalty_factor * power_balance^2;
% 增加机组出力越限罚值
for i = 1:n_gens
P = position(i);
Pmin = gen_data(i, 4);
Pmax = gen_data(i, 5);
if P < Pmin
penalty = penalty + penalty_factor * (Pmin - P)^2;
elseif P > Pmax
penalty = penalty + penalty_factor * (P - Pmax)^2;
end
end
end
% 结果可视化函数
function visualize_results(gen_data, best_position, best_costs, avg_costs, convergence_curve, max_iter)
% 设置图形
figure('Position', [100, 100, 1200, 800], 'Color', 'w');
% 机组出力分布
subplot(2, 2, 1);
n_gens = size(gen_data, 1);
bar(best_position, 'FaceColor', [0.2, 0.6, 0.9]);
hold on;
% 添加上下限线
for i = 1:n_gens
line([i-0.4, i+0.4], [gen_data(i,4), gen_data(i,4)], 'Color', 'r', 'LineWidth', 1.5, 'LineStyle', '--');
line([i-0.4, i+0.4], [gen_data(i,5), gen_data(i,5)], 'Color', 'r', 'LineWidth', 1.5, 'LineStyle', '--');
end
title('机组最优出力分配', 'FontSize', 12);
xlabel('机组编号', 'FontSize', 10);
ylabel('出力 (MW)', 'FontSize', 10);
grid on;
set(gca, 'FontSize', 9);
% 成本曲线
subplot(2, 2, 2);
plot(1:max_iter, best_costs, 'b-', 'LineWidth', 2);
hold on;
plot(1:max_iter, avg_costs, 'r--', 'LineWidth', 1.5);
plot(1:max_iter, convergence_curve, 'g-', 'LineWidth', 2);
title('算法收敛曲线', 'FontSize', 12);
xlabel('迭代次数', 'FontSize', 10);
ylabel('成本 ($/h)', 'FontSize', 10);
legend('当前最优', '种群平均', '全局最优', 'Location', 'best');
grid on;
set(gca, 'FontSize', 9);
% 成本函数构成
subplot(2, 2, 3);
costs = zeros(1, n_gens);
for i = 1:n_gens
P = best_position(i);
a = gen_data(i, 1);
b = gen_data(i, 2);
c = gen_data(i, 3);
costs(i) = a*P^2 + b*P + c;
end
pie(costs, ones(1, n_gens));
title('各机组成本占比', 'FontSize', 12);
legend_labels = arrayfun(@(i) sprintf('机组%d', i), 1:n_gens, 'UniformOutput', false);
legend(legend_labels, 'Location', 'eastoutside');
set(gca, 'FontSize', 9);
% 机组成本特性曲线
subplot(2, 2, 4);
colors = lines(n_gens);
for i = 1:n_gens
Pmin = gen_data(i, 4);
Pmax = gen_data(i, 5);
P_range = linspace(Pmin, Pmax, 100);
cost_curve = gen_data(i,1)*P_range.^2 + gen_data(i,2)*P_range + gen_data(i,3);
plot(P_range, cost_curve, 'Color', colors(i,:), 'LineWidth', 1.5);
hold on;
% 标记最优运行点
P_opt = best_position(i);
cost_opt = gen_data(i,1)*P_opt^2 + gen_data(i,2)*P_opt + gen_data(i,3);
plot(P_opt, cost_opt, 'o', 'MarkerSize', 8, 'MarkerFaceColor', colors(i,:), 'MarkerEdgeColor', 'k');
end
title('机组成本特性曲线', 'FontSize', 12);
xlabel('出力 (MW)', 'FontSize', 10);
ylabel('成本 ($/h)', 'FontSize', 10);
grid on;
legend_labels_opt = arrayfun(@(i) sprintf('机组%d (最优点)', i), 1:n_gens, 'UniformOutput', false);
legend(legend_labels_opt, 'Location', 'eastoutside');
set(gca, 'FontSize', 9);
% 添加标题
sgtitle('粒子群算法求解经济调度问题', 'FontSize', 16, 'FontWeight', 'bold');
% 保存结果
saveas(gcf, 'economic_dispatch_result.png');
end
% 运行程序
pso_economic_dispatch();
说明
MATLAB程序使用粒子群优化算法(PSO)解决电力系统经济调度问题，主要功能包括：
1. 问题建模
- 发电机组成本函数：二次成本函数 \(F_i(P_i) = a_i P_i^2 + b_i P_i + c_i\)
- 系统约束：
- 功率平衡约束：\(\sum P_i = P_d + P_{loss}\)
- 机组出力限制：\(P_{i,min} \leq P_i \leq P_{i,max}\)
- 网损模型：\(P_{loss} = K \times (\sum P_i)^2\)
2. 粒子群算法实现
- 粒子表示：每个粒子代表一个发电计划（各机组出力）
- 约束处理：
- 使用罚函数法处理约束违反
- 专门设计
adjust_generation函数确保功率平衡
- 参数设置：
- 惯性权重线性递减策略
- 速度限制防止粒子移动过快
- 边界检查确保机组出力在可行范围内
参考代码 基于matlab语言，粒子群算法pso，求解经济调度问题 youwenfan.com/contentcnb/79575.html
3. 结果分析
程序提供以下分析结果：
- 各机组最优出力分配
- 算法收敛曲线（全局最优、当前最优和种群平均成本）
- 各机组成本占比饼图
- 机组成本特性曲线及最优运行点
4. 可视化功能
- 机组出力分配柱状图（显示上下限）
- 成本收敛曲线
- 成本构成饼图
- 机组成本特性曲线
算法步骤详解
-
初始化粒子群：
- 随机生成粒子位置（机组出力）
- 调整出力以满足负荷需求
- 初始化粒子速度和成本
-
迭代优化：
- 更新惯性权重（线性递减）
- 更新粒子速度
- 更新粒子位置
- 边界检查确保出力在允许范围内
- 调整出力以满足功率平衡
- 计算成本（包含罚值）
- 更新个体最优和全局最优
-
结果输出：
- 最优发电计划
- 总发电成本
- 网损计算
- 各机组运行状态
-
结果可视化：
- 多种图形展示优化结果
- 保存结果图像
程序输出示例
======= 经济调度最优解 =======
总负荷需求: 1800.00 MW
总发电量: 1803.72 MW
网损: 0.3253 MW
总发电成本: 15450.25 $/h
机组 出力(MW) 成本($/h) 运行状态
------------------------------------
1 436.28 4276.27 运行
2 398.44 3920.19 运行
3 130.00 1322.90 运行
4 240.00 2371.20 运行
5 200.00 1979.00 运行
6 180.00 1580.70 运行
浙公网安备 33010602011771号