基于粒子群优化（PSO）的电力系统经济调度程序实现

一、问题建模与目标函数

电力系统经济调度的核心目标是在满足负荷需求和约束条件下，最小化发电总成本。典型数学模型如下：

1. 目标函数（发电成本最小化）

\(minF=i=1∑N(aiPGi2+biPGi+ci+di⋅ui+ei⋅(PGi−PGi,min)2)\)

• \(N\)：发电机组数量
• \(PGi\)：第i台机组出力（MW）
• \(ai,bi,ci\)：燃料成本系数（二次/一次/固定成本）
• \(di\)：启动成本系数（\(ui=1\)表示启动，0表示停机）
• \(ei\)：爬坡惩罚系数（约束出力变化率）

2. 约束条件

• 功率平衡：\(∑i=1NPGi=PD+PL\)（\(PD\)为负荷需求，\(PL\)为网损）
• 出力上下限：\(PGi,min≤PGi≤PGi,max\)
• 爬坡率约束：\(∣PGit−PGit−1∣≤ΔPGimax\)（\(t\)为时间步长）
• 最小运行时间：\(TGion≥TGi,min\)（连续运行时间约束）

---

二、粒子群算法设计

1. 粒子编码与初始化

• 编码方式：实数编码，粒子位置向量 \(X=[PG1,PG2,...,PGN]\)
• 速度向量：\(V=[v1,v2,...,vN]\)（初始速度随机生成）
• 初始化范围：\(PGi∈[PGi,min,PGi,max]\)（基于机组参数）

2. 适应度函数设计

直接采用目标函数作为适应度函数，同时处理约束：

function fitness = evaluate_fitness(X, a, b, c, d, e, P_D, P_L, Pmin, Pmax, dPmax)
    N = length(X);
    F = 0;
    % 计算燃料成本与启动成本
    for i = 1:N
        F = F + a(i)*X(i)^2 + b(i)*X(i) + c(i);
        if X(i) > Pmin(i)  % 启动成本（假设停机时X(i)=Pmin）
            F = F + d(i);
        end
    end
    % 爬坡惩罚项
    if exist('prev_X', 'var')
        for i = 1:N
            F = F + e(i)*(X(i) - prev_X(i))^2;
        end
    end
    % 功率平衡约束惩罚（软约束）
    balance_error = abs(sum(X) - (P_D + P_L));
    F = F + 1e6 * balance_error;  % 大权重惩罚不平衡
end

3. 粒子更新策略（改进PSO）

• 惯性权重自适应：\(ω(t)=ωmax−(ωmax−ωmin)⋅t/Tmax\)

• 学习因子动态调整：\(c1(t)=c1max−(c1max−c1min)⋅t/Tmax，c2(t)\) 同理

• 混沌扰动：对全局最优粒子施加Tent映射扰动，避免早熟：

  function perturbed_pos = chaos_perturb(pos, scale)
      tent_map = @(x) 2*x .* (x < 0.5) + 2*(1-x) .* (x >= 0.5);
      perturbed_pos = pos + scale * (tent_map(rand(size(pos))) - 0.5);
  end
  

---

三、MATLAB程序实现

1. 主程序框架

% 参数设置
N = 5;                  % 发电机组数量
T = 24;                 % 调度时长（小时）
P_D = 1000;             % 负荷需求（MW）
P_L = 50;               % 网损（MW）

% 机组参数（示例）
a = [0.001, 0.0015, 0.002, 0.0012, 0.0018];  % 二次成本系数
b = [10, 12, 8, 15, 9];                     % 一次成本系数
c = [500, 600, 400, 700, 550];              % 固定成本
d = [1000, 1200, 800, 1500, 900];           % 启动成本
e = [5, 4, 6, 3, 7];                        % 爬坡惩罚系数
Pmin = [100, 80, 120, 90, 110];             % 最小出力
Pmax = [400, 350, 500, 450, 300];           % 最大出力
dPmax = [50, 40, 60, 30, 70];               % 爬坡率限制(MW/h)

% PSO参数
n_particles = 30;       % 粒子数量
max_iter = 100;         % 最大迭代次数
w_max = 0.9; w_min = 0.4; % 惯性权重范围
c1_max = 2.0; c1_min = 0.5; % 个体学习因子
c2_max = 2.0; c2_min = 0.5; % 群体学习因子

% 初始化粒子群
[positions, velocities] = initialize_particles(N, n_particles, Pmin, Pmax);
pbest_pos = positions;    % 个体最优位置
pbest_fit = inf(n_particles, 1); % 个体最优适应度
gbest_pos = zeros(1, N);  % 全局最优位置
gbest_fit = inf;          % 全局最优适应度

% 主循环
for iter = 1:max_iter
    % 更新惯性权重与学习因子
    w = w_max - (w_max - w_min) * iter/max_iter;
    c1 = c1_max - (c1_max - c1_min) * iter/max_iter;
    c2 = c2_min + (c2_max - c2_min) * iter/max_iter;
    
    % 评估适应度
    for i = 1:n_particles
        current_fit = evaluate_fitness(positions(i,:), a, b, c, d, e, P_D, P_L, Pmin, Pmax, dPmax, positions(i,:), iter>1);
        if current_fit < pbest_fit(i)
            pbest_fit(i) = current_fit;
            pbest_pos(i,:) = positions(i,:);
        end
        if current_fit < gbest_fit
            gbest_fit = current_fit;
            gbest_pos = positions(i,:);
        end
    end
    
    % 更新粒子速度与位置
    for i = 1:n_particles
        r1 = rand(1, N);
        r2 = rand(1, N);
        velocities(i,:) = w*velocities(i,:) + ...
                          c1*r1.*(pbest_pos(i,:) - positions(i,:)) + ...
                          c2*r2.*(gbest_pos - positions(i,:));
        positions(i,:) = positions(i,:) + velocities(i,:);
        
        % 边界处理（越界则重置）
        positions(i,:) = max(positions(i,:), Pmin);
        positions(i,:) = min(positions(i,:), Pmax);
    end
    
    % 混沌扰动（每10代）
    if mod(iter,10) == 0
        gbest_pos = chaos_perturb(gbest_pos, 0.1*(Pmax-Pmin));
        gbest_pos = max(gbest_pos, Pmin);
        gbest_pos = min(gbest_pos, Pmax);
    end
end

% 输出结果
disp(['最优总成本: ', num2str(gbest_fit)]);
disp('各机组最优出力:');
disp(gbest_pos');

2. 辅助函数（初始化与适应度计算）

function [positions, velocities] = initialize_particles(N, n_particles, Pmin, Pmax)
    positions = zeros(n_particles, N);
    velocities = zeros(n_particles, N);
    for i = 1:n_particles
        positions(i,:) = Pmin + (Pmax - Pmin).*rand(1, N);
        velocities(i,:) = 0.1*(Pmax - Pmin).*randn(1, N); % 初始速度
    end
end

function fitness = evaluate_fitness(X, a, b, c, d, e, P_D, P_L, Pmin, Pmax, dPmax, prev_X, is_iter>1)
    N = length(X);
    F = 0;
    % 燃料与固定成本
    for i = 1:N
        F = F + a(i)*X(i)^2 + b(i)*X(i) + c(i);
        if X(i) > Pmin(i)  % 启动成本（假设停机时X(i)=Pmin）
            F = F + d(i);
        end
    end
    % 爬坡惩罚（仅迭代后期）
    if is_iter>1
        for i = 1:N
            F = F + e(i)*(X(i) - prev_X(i))^2;
        end
    end
    % 功率平衡惩罚（软约束）
    balance_error = abs(sum(X) - (P_D + P_L));
    F = F + 1e6 * balance_error;  % 大权重确保平衡
end

参考代码 粒子群进行电力系统经济调度程序 www.youwenfan.com/contentcnd/82865.html

粒子群算法可有效解决电力系统经济调度问题，在保证约束的同时显著降低发电成本。实际应用中需根据具体电网结构调整参数，并结合实时数据进行在线优化。