电力系统优化：二进制粒子群算法 (BPSO) 求解配电网重构

电力系统优化：二进制粒子群算法 (BPSO) 求解配电网重构

使用二进制粒子群算法(BPSO)解决电力系统中的配电网重构问题。该算法旨在通过优化开关状态来最小化电网功率损耗。

%% 电力系统配电网重构 - 二进制粒子群算法 (BPSO)
% 作者: MATLAB技术
% 描述: 使用BPSO算法优化配电网开关状态以最小化功率损耗
% 日期: 2023年11月15日

function BPSO_PowerSystem()
clc;
clear;
close all;
warning('off', 'all');

%% 参数设置
display('正在初始化参数...');
% 算法参数
swarmSize = 50;         % 粒子群大小
maxIter = 200;          % 最大迭代次数
c1 = 2.0;               % 个体学习因子
c2 = 2.0;               % 社会学习因子
inertia = 1.0;          % 初始惯性权重
inertiaMin = 0.4;       % 最小惯性权重
inertiaMax = 0.9;       % 最大惯性权重

% 电力系统参数
numSwitches = 20;       % 开关数量
baseLoad = 10;          % 基础负载 (MW)
voltageBase = 11;       % 基础电压 (kV)
resistancePerKm = 0.5;  % 线路电阻 (Ω/km)

%% 初始化电力系统
display('正在初始化电力系统模型...');
% 创建配电网拓扑结构
[network, lines] = createDistributionNetwork(numSwitches);

% 随机生成负载分布
loads = baseLoad * (0.8 + 0.4*rand(1, numSwitches+1)); % +1 为变电站

% 计算初始功率损耗
initialLoss = calculatePowerLoss(network, lines, loads, voltageBase, resistancePerKm);
fprintf('初始功率损耗: %.4f MW\n', initialLoss);

%% 初始化粒子群
display('正在初始化粒子群...');
% 初始化粒子位置（二进制开关状态）
particles = rand(swarmSize, numSwitches) > 0.5;

% 初始化粒子速度
velocities = 0.1 * randn(swarmSize, numSwitches);

% 初始化个体最优位置和适应度值
personalBestPos = particles;
personalBestFit = inf(swarmSize, 1);  % 初始化为最大值

% 初始化全局最优
globalBestFit = inf;
globalBestPos = zeros(1, numSwitches);

%% 计算初始适应度
display('正在计算初始适应度...');
for i = 1:swarmSize
    % 获取当前开关状态
    switchState = particles(i, :);
    
    % 更新网络拓扑
    updatedNetwork = updateNetworkTopology(network, switchState);
    
    % 计算适应度（功率损耗）
    [fitness, voltageProfile] = calculateFitness(updatedNetwork, lines, loads, ...
                                                voltageBase, resistancePerKm);
    
    % 更新个体最优
    if fitness < personalBestFit(i)
        personalBestFit(i) = fitness;
        personalBestPos(i, :) = switchState;
    end
    
    % 更新全局最优
    if fitness < globalBestFit
        globalBestFit = fitness;
        globalBestPos = switchState;
        bestVoltage = voltageProfile;
    end
end

%% BPSO 主循环
display('开始BPSO优化...');
% 存储收敛曲线
convergence = zeros(maxIter, 1);
convergence(1) = globalBestFit;

% 创建进度条
h = waitbar(0, 'BPSO优化进行中...', 'Name', '配电网重构');

for iter = 2:maxIter
    % 更新惯性权重（线性递减）
    inertia = inertiaMax - (inertiaMax - inertiaMin) * iter / maxIter;
    
    for i = 1:swarmSize
        % 更新速度
        r1 = rand(1, numSwitches);
        r2 = rand(1, numSwitches);
        
        cognitive = c1 * r1 .* (personalBestPos(i, :) - particles(i, :));
        social = c2 * r2 .* (globalBestPos - particles(i, :));
        velocities(i, :) = inertia * velocities(i, :) + cognitive + social;
        
        % 限制速度范围
        velocities(i, :) = min(max(velocities(i, :), -6), 6);
        
        % 更新位置（使用sigmoid函数）
        sigmoidV = 1 ./ (1 + exp(-velocities(i, :)));
        particles(i, :) = rand(1, numSwitches) < sigmoidV;
        
        % 计算适应度
        switchState = particles(i, :);
        updatedNetwork = updateNetworkTopology(network, switchState);
        [fitness, voltageProfile] = calculateFitness(updatedNetwork, lines, loads, ...
                                                    voltageBase, resistancePerKm);
        
        % 更新个体最优
        if fitness < personalBestFit(i)
            personalBestFit(i) = fitness;
            personalBestPos(i, :) = switchState;
        end
        
        % 更新全局最优
        if fitness < globalBestFit
            globalBestFit = fitness;
            globalBestPos = switchState;
            bestVoltage = voltageProfile;
        end
    end
    
    % 存储收敛值
    convergence(iter) = globalBestFit;
    
    % 更新进度条
    waitbar(iter/maxIter, h, sprintf('迭代: %d/%d, 最小损耗: %.4f MW', iter, maxIter, globalBestFit));
    
    % 显示每20次迭代的信息
    if mod(iter, 20) == 0
        fprintf('迭代 %d: 最小功率损耗 = %.4f MW\n', iter, globalBestFit);
    end
end

% 关闭进度条
close(h);

%% 结果显示
display('优化完成!');
fprintf('\n======== 优化结果 ========\n');
fprintf('初始功率损耗: %.4f MW\n', initialLoss);
fprintf('优化后功率损耗: %.4f MW\n', globalBestFit);
fprintf('功率损耗降低: %.2f%%\n', (initialLoss - globalBestFit)/initialLoss*100);
fprintf('最优开关状态: \n');
disp(globalBestPos);

% 绘制收敛曲线
figure('Position', [100, 100, 800, 400]);
plot(convergence, 'LineWidth', 2);
title('BPSO收敛曲线', 'FontSize', 14);
xlabel('迭代次数', 'FontSize', 12);
ylabel('功率损耗 (MW)', 'FontSize', 12);
grid on;
set(gca, 'FontSize', 10);

% 绘制电压分布图
figure('Position', [100, 100, 800, 400]);
plot(1:length(bestVoltage), bestVoltage, 'o-', 'LineWidth', 2, 'MarkerFaceColor', 'b');
title('节点电压分布', 'FontSize', 14);
xlabel('节点编号', 'FontSize', 12);
ylabel('电压 (p.u.)', 'FontSize', 12);
ylim([0.9, 1.05]);
grid on;
set(gca, 'FontSize', 10);

% 绘制网络拓扑
figure('Position', [100, 100, 800, 600]);
plotNetworkTopology(network, globalBestPos);
title('优化后网络拓扑', 'FontSize', 14);
end

%% 辅助函数：创建配电网拓扑
function [network, lines] = createDistributionNetwork(numSwitches)
% 创建基本拓扑结构
network = struct();
network.numNodes = numSwitches + 1; % 节点数 = 开关数 + 1（变电站）
network.substation = 1; % 变电站节点

% 初始化连接矩阵
network.connections = zeros(network.numNodes);

% 创建线路
lines = cell(numSwitches, 1);
for i = 1:numSwitches
    % 随机确定线路长度 (1-10 km)
    length = 1 + 9*rand();
    
    % 随机确定起点和终点
    fromNode = randi([1, network.numNodes]);
    toNode = randi([1, network.numNodes]);
    while toNode == fromNode % 确保起点和终点不同
        toNode = randi([1, network.numNodes]);
    end
    
    % 添加到连接矩阵
    network.connections(fromNode, toNode) = i;
    network.connections(toNode, fromNode) = i;
    
    % 存储线路信息
    lines{i} = struct(...
        'id', i, ...
        'from', fromNode, ...
        'to', toNode, ...
        'length', length, ...
        'status', true); % 初始状态为闭合
end
end

%% 辅助函数：更新网络拓扑
function updatedNetwork = updateNetworkTopology(network, switchState)
% 复制原始网络
updatedNetwork = network;

% 更新开关状态
for i = 1:length(switchState)
    if ~switchState(i) % 如果开关状态为0，断开该线路
        from = find(updatedNetwork.connections == i, 1, 'first');
        [row, col] = find(updatedNetwork.connections == i);
        if ~isempty(row)
            updatedNetwork.connections(row(1), col(1)) = 0;
            updatedNetwork.connections(col(1), row(1)) = 0;
        end
    end
end
end

%% 辅助函数：计算适应度（功率损耗）
function [powerLoss, voltageProfile] = calculateFitness(network, lines, loads, voltageBase, resistancePerKm)
% 简化的功率流计算
numNodes = network.numNodes;

% 初始化节点电压（标幺值）
voltageProfile = ones(1, numNodes);

% 计算每条线路的电流和功率损耗
totalLoss = 0;
for i = 1:length(lines)
    line = lines{i};
    
    % 检查线路是否连接
    if network.connections(line.from, line.to) == line.id
        % 计算线路电流（简化计算）
        current = (loads(line.from) + loads(line.to)) / (voltageBase * 1000); % kA
        
        % 计算线路电阻
        resistance = resistancePerKm * line.length;
        
        % 计算线路功率损耗
        lineLoss = current^2 * resistance;
        totalLoss = totalLoss + lineLoss;
        
        % 更新节点电压（简化计算）
        voltageDrop = current * resistance;
        voltageProfile(line.to) = voltageProfile(line.from) - voltageDrop / (voltageBase * 1000);
    end
end

% 添加电压约束惩罚
voltageViolation = sum(max(0, 0.95 - voltageProfile)) + sum(max(0, voltageProfile - 1.05));
powerLoss = totalLoss / 1000 + 10 * voltageViolation; % 转换为MW并添加惩罚项
end

%% 辅助函数：计算功率损耗
function powerLoss = calculatePowerLoss(network, lines, loads, voltageBase, resistancePerKm)
[powerLoss, ~] = calculateFitness(network, lines, loads, voltageBase, resistancePerKm);
end

%% 辅助函数：绘制网络拓扑
function plotNetworkTopology(network, switchState)
numNodes = network.numNodes;
connections = network.connections;

% 创建节点位置
theta = linspace(0, 2*pi, numNodes+1);
theta(end) = [];
x = cos(theta);
y = sin(theta);

% 绘制节点
scatter(x, y, 100, 'filled', 'MarkerFaceColor', [0.2 0.6 0.9]);
hold on;
text(x(1), y(1), '变电站', 'VerticalAlignment', 'bottom', 'HorizontalAlignment', 'center');

% 绘制节点标签
for i = 2:numNodes
    text(x(i), y(i), sprintf('节点%d', i), 'VerticalAlignment', 'bottom', 'HorizontalAlignment', 'center');
end

% 绘制线路
lineCount = 0;
for i = 1:numNodes
    for j = i+1:numNodes
        if connections(i, j) > 0
            lineId = connections(i, j);
            lineCount = lineCount + 1;
            
            % 根据开关状态选择线条样式
            if switchState(lineId)
                lineStyle = '-';
                lineColor = [0 0.5 0];
                lineWidth = 2;
            else
                lineStyle = ':';
                lineColor = [0.8 0.2 0.2];
                lineWidth = 1.5;
            end
            
            % 绘制线路
            plot([x(i), x(j)], [y(i), y(j)], lineStyle, ...
                'Color', lineColor, 'LineWidth', lineWidth);
            
            % 添加线路标签
            midX = (x(i) + x(j)) / 2;
            midY = (y(i) + y(j)) / 2;
            text(midX, midY, sprintf('S%d', lineId), ...
                'BackgroundColor', 'w', 'Margin', 1, ...
                'FontSize', 8, 'HorizontalAlignment', 'center');
        end
    end
end

% 添加图例
h = zeros(2, 1);
h(1) = plot(NaN, NaN, '-', 'Color', [0 0.5 0], 'LineWidth', 2);
h(2) = plot(NaN, NaN, ':', 'Color', [0.8 0.2 0.2], 'LineWidth', 1.5);
legend(h, {'闭合开关', '断开开关'}, 'Location', 'best');

axis equal;
axis off;
title(sprintf('网络拓扑 (节点: %d, 开关: %d)', numNodes, length(switchState)));
hold off;
end

算法说明

1. 二进制粒子群算法 (BPSO) 核心原理

BPSO算法在标准粒子群优化基础上进行了修改，适用于解决二进制优化问题：

• 位置更新：使用sigmoid函数将速度转换为概率值

  sigmoidV = 1 ./ (1 + exp(-velocities(i, :)));
  particles(i, :) = rand(1, numSwitches) < sigmoidV;
  

• 速度更新：与传统PSO类似，但位置是二进制值

  velocities(i, :) = inertia * velocities(i, :) + cognitive + social;
  

• 惯性权重：线性递减策略平衡全局和局部搜索

  inertia = inertiaMax - (inertiaMax - inertiaMin) * iter / maxIter;
  

2. 电力系统建模

配电网模型包含以下关键组件：

• 变电站：电力系统的起点
• 节点：电力分配点
• 开关：控制线路连接状态
• 线路：具有电阻属性的连接

3. 适应度函数

适应度函数计算功率损耗并考虑电压约束：

function [powerLoss, voltageProfile] = calculateFitness(...)
    % 计算线路电流和功率损耗
    current = (loads(line.from) + loads(line.to)) / (voltageBase * 1000);
    resistance = resistancePerKm * line.length;
    lineLoss = current^2 * resistance;
    
    % 电压约束惩罚
    voltageViolation = sum(max(0, 0.95 - voltageProfile)) + ...
                      sum(max(0, voltageProfile - 1.05));
    powerLoss = totalLoss / 1000 + 10 * voltageViolation;
end

4. 可视化功能

算法提供三种可视化效果：

1. 收敛曲线：展示优化过程中功率损耗的变化
2. 电压分布：显示优化后各节点的电压水平
3. 网络拓扑：直观展示开关状态和网络连接

运行结果示例

运行该算法后，您将看到：

1. 控制台显示初始和优化后的功率损耗
2. 收敛曲线展示优化过程
3. 电压分布图显示各节点电压
4. 网络拓扑图展示开关状态

参数调整建议

1. 粒子群参数：

   • 增加swarmSize可提高搜索能力但增加计算时间
   • 调整inertiaMin和inertiaMax可平衡探索与开发

2. 电力系统参数：

   • 修改baseLoad调整系统负载水平
   • 调整resistancePerKm改变线路特性

3. 约束条件：

   • 在calculateFitness函数中修改电压约束范围
   • 调整惩罚系数以改变约束严格程度

参考 用二进制粒子群算法求解电力系统的BPSO算法源代码