改进粒子群算法求解电力系统经济调度问题（Matlab实现）

最新推荐文章于 2024-03-01 12:52:47 发布

程序辅导帮

最新推荐文章于 2024-03-01 12:52:47 发布

阅读量373

点赞数

文章标签： matlab 算法开发语言

本文链接：https://blog.csdn.net/2201_75454341/article/details/128091719

版权

1 相关知识点

2 Matlab完整代码实现

3 结果及可视化

1 相关知识点

这里写图片描述

2 Matlab代码实现

知识点讲解完毕，下面就是Matlab部分代码：

clc;
clear;
close all;

%% 经济调度问题

extmodel=CreateModel();

CostFunction=@(x) MyCostExt(x,extmodel); % 成本函数（目标函数）

nVar=extmodel.nPlant;             % 发电机台数（决策变量的个数）

VarSize=[1 nVar];   % 决策变量矩阵的大小

VarMin=0;         % 变量下限
VarMax=1;         % 变量上限


%% 粒子群算法相关参数

MaxIt=100;      % 最大迭代次数

nPop=10;        % 总群数量

% w=1;            % 惯性权重
% wdamp=0.99;     % 惯性重量阻尼比
% c1=2;           % 个体学习系数
% c2=2;           % 种群学习系数

%% 约束系数
phi1=2.05;
phi2=2.05;
phi=phi1+phi2;
chi=2/(phi-2+sqrt(phi^2-4*phi));
w=chi;          % 惯性权重
wdamp=1;        % 惯性重量阻尼比
c1=chi*phi1;    % 个体学习系数
c2=chi*phi2;    % 种群学习系数

%% 飞行速度限制
VelMax=0.1*(VarMax-VarMin);
VelMin=-VelMax;

%% 初始化

empty_particle.Position=[];
empty_particle.Cost=[];
empty_particle.Out=[];
empty_particle.Velocity=[];
empty_particle.Best.Position=[];
empty_particle.Best.Cost=[];
empty_particle.Best.Out=[];

particle=repmat(empty_particle,nPop,1);

BestSol.Cost=inf;

for i=1:nPop
    
    %=====初始化粒子群位置===============
    particle(i).Position=unifrnd(VarMin,VarMax,VarSize);
    
    %=====初始化速度======
    particle(i).Velocity=zeros(VarSize);
    
    %=====目标函数计算===========
    [particle(i).Cost, particle(i).Out]=CostFunction(particle(i).Position);
    
    %====更新粒子个体最优=====
    particle(i).Best.Position=particle(i).Position;
    particle(i).Best.Cost=particle(i).Cost;
    particle(i).Best.Out=particle(i).Out;
    
    %====更新粒子群全局最优========
    if particle(i).Best.Cost<BestSol.Cost
        
        BestSol=particle(i).Best;
        
    end
    
end

BestCost=zeros(MaxIt,1);


%% PSO 主循环

for it=1:MaxIt
    
    for i=1:nPop
        
        %============更新速度（跟着公式写就可以啦）===============
        particle(i).Velocity = w*particle(i).Velocity ...
            +c1*rand(VarSize).*(particle(i).Best.Position-particle(i).Position) ...
            +c2*rand(VarSize).*(BestSol.Position-particle(i).Position);
        
        %============适用速度限制============
        particle(i).Velocity = max(particle(i).Velocity,VelMin);
        particle(i).Velocity = min(particle(i).Velocity,VelMax);
        
        %============更新位置================
        particle(i).Position = particle(i).Position + particle(i).Velocity;
        
        IsOutside=(particle(i).Position<VarMin | particle(i).Position>VarMax);
        particle(i).Velocity(IsOutside)=-particle(i).Velocity(IsOutside);
        
        %============适用位置限制============
        particle(i).Position = max(particle(i).Position,VarMin);
        particle(i).Position = min(particle(i).Position,VarMax);
        
        %============计算目标函数===========
        [particle(i).Cost, particle(i).Out] = CostFunction(particle(i).Position);
        
        %==========更新个体最优==========
        if particle(i).Cost<particle(i).Best.Cost
            
            particle(i).Best.Position=particle(i).Position;
            particle(i).Best.Cost=particle(i).Cost;
            particle(i).Best.Out=particle(i).Out;
            
            %=======更新全局最优============
            if particle(i).Best.Cost<BestSol.Cost
                
                BestSol=particle(i).Best;
                
            end
            
        end
        
    end
    
    BestCost(it)=BestSol.Cost;
    
    disp(['迭代次数' num2str(it) ': 最优解为 = ' num2str(BestCost(it))]);
    
    w=w*wdamp;
    
end

%% 结果

figure;
plot(BestCost,'LineWidth',2);
xlabel('迭代次数');
ylabel('最优解');

%% 目标函数
function [z, out]=MyCost(x,model)

    P=ParseSolution(x,model);
    
    out=ModelCalculations(P,model);

    z=out.z;

end

%% 约束条件：功率平衡
function out=ModelCalculations(P,model)

    alpha=model.Plants.alpha;
    beta=model.Plants.beta;
    gamma=model.Plants.gamma;
    
    C=alpha+beta.*P+gamma.*P.*P;
    
    CTotal=sum(C);
    
    B=model.B;
    B0=model.B0;
    B00=model.B00;
    
    PL=P*B*P'+B0*P'+B00;
    
    PTotal=sum(P);
    
    PD=model.PD;
    
    PowerBalanceViolation=max(1-(PTotal-PL)/PD,0);
    
    q=1000000;
    z=CTotal+q*PowerBalanceViolation;
    
    out.P=P;
    out.PTotal=PTotal;
    out.C=C;
    out.CTotoal=CTotal;
    out.PL=PL;
    out.PowerBalanceViolation=PowerBalanceViolation;
    out.z=z;

end

%% 约束条件：功率平衡
function out=ModelCalculations(P,model)

alpha=model.Plants.alpha;
beta=model.Plants.beta;
gamma=model.Plants.gamma;

C=alpha+beta.*P+gamma.*P.*P;

CTotal=sum(C);

B=model.B;
B0=model.B0;
B00=model.B00;

PL=P*B*P'+B0*P'+B00;

PTotal=sum(P);

PD=model.PD;

PowerBalanceViolation=max(1-(PTotal-PL)/PD,0);

q=1000000;
z=CTotal+q*PowerBalanceViolation;

out.P=P;
out.PTotal=PTotal;
out.C=C;
out.CTotoal=CTotal;
out.PL=PL;
out.PowerBalanceViolation=PowerBalanceViolation;
out.z=z;

end

3 结果及可视化

function intmodel=CreateInternalModel(k,extmodel)

intmodel.PD=extmodel.PD;

nPlant=numel(k);

intmodel.nPlant=nPlant;

Plants=extmodel.Plants;
EmptyPlant=extmodel.EmptyPlant;

intmodel.Plants.Pmin=zeros(1,nPlant);
intmodel.Plants.Pmax=zeros(1,nPlant);
intmodel.Plants.alpha=zeros(1,nPlant);
intmodel.Plants.beta=zeros(1,nPlant);
intmodel.Plants.gamma=zeros(1,nPlant);
intmodel.Plants.PZ=cell(1,nPlant);

for i=1:nPlant
ki=k(i);
if ki==0
% Set Plant #i to be Empty
intmodel.Plants.Pmin(i)=EmptyPlant.Pmin;
intmodel.Plants.Pmax(i)=EmptyPlant.Pmax;
intmodel.Plants.alpha(i)=EmptyPlant.alpha;
intmodel.Plants.beta(i)=EmptyPlant.beta;
intmodel.Plants.gamma(i)=EmptyPlant.gamma;
intmodel.Plants.PZ{i}=EmptyPlant.PZ;
else

intmodel.Plants.Pmin(i)=Plants.Pmin(ki);
intmodel.Plants.Pmax(i)=Plants.Pmax(ki);
intmodel.Plants.alpha(i)=Plants.alpha(ki);
intmodel.Plants.beta(i)=Plants.beta(ki);
intmodel.Plants.gamma(i)=Plants.gamma(ki);
intmodel.Plants.PZ{i}=Plants.PZ{ki};
end
end

intmodel.B=extmodel.B;
intmodel.B0=extmodel.B0;
intmodel.B00=extmodel.B00;

end