优化调度基于matlab粒子群算法求解水火电调度优化问题含matlab源码1181期(代码片段)

一、简介

粒子群算法源于复杂适应系统（Complex Adaptive System,CAS）。CAS理论于1994年正式提出，CAS中的成员称为主体。比如研究鸟群系统，每个鸟在这个系统中就称为主体。主体有适应性，它能够与环境及其他的主体进行交流，并且根据交流的过程“学习”或“积累经验”改变自身结构与行为。整个系统的演变或进化包括：新层次的产生（小鸟的出生）；分化和多样性的出现（鸟群中的鸟分成许多小的群）；新的主题的出现（鸟寻找食物过程中，不断发现新的食物）。
所以CAS系统中的主体具有4个基本特点（这些特点是粒子群算法发展变化的依据）：
首先，主体是主动的、活动的。
主体与环境及其他主体是相互影响、相互作用的，这种影响是系统发展变化的主要动力。
环境的影响是宏观的，主体之间的影响是微观的，宏观与微观要有机结合。
最后，整个系统可能还要受一些随机因素的影响。
粒子群算法就是对一个CAS系统－－－鸟群社会系统的研究得出的。
粒子群算法（ Particle Swarm Optimization, PSO）最早是由Eberhart和Kennedy于1995年提出，它的基本概念源于对鸟群觅食行为的研究。设想这样一个场景:一群鸟在随机搜寻食物，在这个区域里只有一块食物，所有的鸟都不知道食物在哪里，但是它们知道当前的位置离食物还有多远。那么找到食物的最优策略是什么呢?最简单有效的就是搜寻目前离食物最近的鸟的周围区域。
PSO算法就从这种生物种群行为特性中得到启发并用于求解优化问题。在PSO中，每个优化问题的潜在解都可以想象成d维搜索空间上的一个点，我们称之为“粒子”（Particle），所有的粒子都有一个被目标函数决定的适应值(Fitness Value )，每个粒子还有一个速度决定他们飞翔的方向和距离，然后粒子们就追随当前的最优粒子在解空间中搜索。Reynolds对鸟群飞行的研究发现。鸟仅仅是追踪它有限数量的邻居但最终的整体结果是整个鸟群好像在一个中心的控制之下.即复杂的全局行为是由简单规则的相互作用引起的。

二、源代码

clear all
clc

%opt计算相关量
M=20;%粒子个数
I=4;%水电站个数
T=24;%时间
K=400;%迭代次数
lamda=100;%罚函数惩罚因子
error=0.01;%误差精度
wini=0.9;wend=0.4;c1=2;c2=2;%粒子群算法惯性因子和学习因子
x=zeros(I,T,M);%粒子位置
v=zeros(I,T,M);%粒子速度
vmax=zeros(I,1);%粒子最大速度
alpha=60;beta=-1.355;gama=0.0105;yeta=0.4968;deta=0.01925;%目标函数e的系数
e=zeros(M,1);%粒子适应值
pbest=zeros(I,T,M);%个体最好位置
epbest=zeros(M,1);%个体最好适应值
gbest=ones(1,1);%全局最好序号
as=5000;bs=19.2;cs=0.002;%煤耗成本系数
f=0;%煤耗成本

%负载相关量
Pd=[1370 1390 1360 1290 1290 1410 1650 2000 2240 2320 2230 2310 2230 2200 2130 2070 2130 2140 2240 2280 2240 2120 1850 1590];%24小时的负载功率

%火电站相关量
Ps=zeros(1,T);%火电输出功率
Psmax=2500;%火电输出最大功率
Psmin=500;%火电输出最小功率
Psbest=zeros(1,T);%火电输出功率最优解

%水电站相关量
Vh=zeros(I,T);%水电站库容量
Vhmax=[150 120 240 160];%水电站库容量最大值
Vhmin=[80 60 100 70];%水电站库容量最小值
Vhini=[100 80 170 120];%水电站库容量初始值
Vhend=[120 70 170 140];%水电站库容量最终值
Vhx=ones(I,1);%调整发电流量Qh时，水电站库容量计算值与最终值的差值
Vhbest=zeros(I,T);%水电站库容量最优解
Qh=zeros(I,T);%发电流量
Qhmax=[15 15 30 25];%发电流量最大值
Qhmin=[5 6 10 13];%发电流量最小值
Qhbest=zeros(I,T);%发电流量最优解
Ph=zeros(I,T);%水电输出功率
Phsum=zeros(1,T);%水电输出功率之和
Phmax=500;%水电输出功率最大值
Phmin=0;%水电输出功率最小值
Phbest=zeros(I,T);%水电输出功率最优解
Phsumbest=zeros(1,T);%水电输出功率最优解之和
Ih=[10   9   8   7   6   7   8   9   10  11  12  10  11  12  11  10   9   8   7   6   7   8   9   10;
     8   8   9   9   8   7   6   7    8   9   9   8   8   9   9   8   7   6   7   8   9   9   8    8;
    8.1 8.2  4   2   3   4   3   2    1   1   1   2   4   3   3   2   2   2   1   1   2   2   1    0;
    2.8 2.4 1.6  0   0   0   0   0    0   0   0   0   0   0   0   0   0   0   0   0   0   0   0    0];%水电站径流量
c=[-0.0042 -0.42 0.030 0.90 10.0 -50;
     -0.0040 -0.30 0.015 1.14 9.5 -70;
     -0.0016 -0.30 0.014 0.55 5.5 -40;
     -0.0030 -0.31 0.027 1.44 14.0 -90];%水电站输出功率系数
tao=[2 3 4 0];%水流时滞系数

%/**********************初始化**********************/%
for m=1:M
    %初始化Qh
    for i=1:I
        for t=1:T
            Qh(i,t)=Qhmin(i)+rand()*(Qhmax(i)-Qhmin(i));
        end
    end
    [Qh,Ps,Vh,Ph]=adjust(Qh);
    %初始化x，v，vmax
    for i=1:I
        x(i,1:T,m)=Qh(i,1:T);
        vmax(i,1)=(Qhmax(i)-Qhmin(i))/5;
    end
    for i=1:I
        for t=1:T
            v(i,t,m)=rand()*vmax(i,1);
            if rand()<0.5
                v(i,t,m)=-1*v(i,t,m);
            end
        end
    end
    [e(m,1)]=aime(Ps,Vh,Ph);
end
%初始化个体最好
pbest=x;
%初始化个体最好适应值
epbest=e;
%min (x-3)^2 
%s.t. 1<x<4.
clear all
clc
x=zeros(20,1);%粒子位置
v=zeros(20,1);%粒子速度
vmax=(4-1)./5;%最大速度
fx=zeros(20,1);%粒子适应值
pbest=zeros(20,1);%个体最好位置
fpbest=zeros(20,1);%个体最好适应值
gbest=ones(1,1);%全局最好序号
w=0.9;%惯性权重
c1=2.0;%学习因子
c2=2.0;
gmax=200;%最大迭代次数

%初始化
for i=1:20
    x(i,1)=1+rand().*(4-1);%初始化粒子位置
    tmp = rand().*vmax;%初始化速度
    v(i,1)=tmp;
    if rand()<0.5
        v(i,1)=-1.*tmp;
    end
    fx(i)=(x(i,1)-3)^2;%每个个体目标函数值
    pbest(i,1)=x(i,1);%初始化个体最好
    fpbest(i,1)=fx(i);%初始化个体最好适应值
end
t=1;
while t<gmax
    w=0.9-0.5.*t/gmax;%更新惯性权重
    for i=1:20
        v(i,1)=w.*v(i,1)+c1.*rand().*(pbest(i,1)-x(i,1))+c2.*rand().*(pbest(gbest,1)-x(i,1));%速度更新
        if v(i,1)>vmax%速度越界判断
            v(i,1)=vmax;
        end
        if v(i,1)<-1.*vmax
            v(i,1)=-1*vmax;
        end
        x(i,1)=x(i,1)+v(i,1);%位置更新
        if x(i,1)>4%限幅
            x(i,1)=4;
        end
        if x(i,1)<1
            x(i,1)=1;
        end
function [Qhbest,Phbest,Vhbest,Psbest,Phsumbest]=best(x,gbest)
%opt计算相关量
I=4;%水电站个数
T=24;%时间

%负载相关量
Pd=[750 780 700 650 670 800 950 1010 1090 1080 1100 1150 1110 1030 1010 1060 1050 1120 1070 1050 910 860 850 800];%24小时的负载功率

%火电站相关量
Psbest=zeros(3,T);%火电输出功率最优解
Psum=zeros(1,T);
Psmax=[175 300 500];%火电输出最大功率
Psmin=[20 40 50];%火电输出最小功率
Pssumbest=zeros(1,T);
%水电站相关量
Vhini=[100 80 170 120];%水电站库容量初始值
Vhbest=zeros(I,T);%水电站库容量最优解
Qhbest=zeros(I,T);%发电流量最优解
Phmin=0;%水电输出功率最小值
Phbest=zeros(I,T);%水电输出功率最优解
Phsumbest=zeros(1,T);%水电输出功率最优解之和
Ih=[10   9   8   7   6   7   8   9   10  11  12  10  11  12  11  10   9   8   7   6   7   8   9   10;
     8   8   9   9   8   7   6   7    8   9   9   8   8   9   9   8   7   6   7   8   9   9   8    8;
    8.1 8.2  4   2   3   4   3   2    1   1   1   2   4   3   3   2   2   2   1   1   2   2   1    0;
    2.8 2.4 1.6  0   0   0   0   0    0   0   0   0   0   0   0   0   0   0   0   0   0   0   0    0];%水电站径流量
c=[-0.0042 -0.42 0.030 0.90 10.0 -50;
     -0.0040 -0.30 0.015 1.14 9.5 -70;
     -0.0016 -0.30 0.014 0.55 5.5 -40;
     -0.0030 -0.31 0.027 1.44 14.0 -90];%水电站输出功率系数
tao=[2 3 4 0];%水流时滞系数
%计算Qh最优解
for i=1:I
    Qhbest(i,1:T)=x(i,1:T,gbest);
end
%计算Vh最优解
for i=1:I
    for t=1:T
        if t==1
            Vhbest(i,t)=Vhini(i)+Ih(i,t)-Qhbest(i,t);
        else
            Vhbest(i,t)=Vhbest(i,t-1)+Ih(i,t)-Qhbest(i,t);
        end
        if i==3 
            if t>tao(1)
                Vhbest(i,t)=Vhbest(i,t)+Qhbest(1,t-tao(1));
            end
            if t>tao(2)
                Vhbest(i,t)=Vhbest(i,t)+Qhbest(2,t-tao(2));
            end
        end
        if i==4
            if t>tao(3)
                Vhbest(i,t)=Vhbest(i,t)+Qhbest(3,t-tao(3));
            end
        end
    end
end
%计算Ph最优解
 for i=1:I
        for t=1:T
            Phbest(i,t)=c(i,1)*Vhbest(i,t)*Vhbest(i,t)+c(i,2)*Qhbest(i,t)*Qhbest(i,t)+c(i,3)*Vhbest(i,t)*Qhbest(i,t)+c(i,4)*Vhbest(i,t)+c(i,5)*Qhbest(i,t)+c(i,6);
            if Phbest(i,t)<Phmin
                Phbest(i,t)=Phmin;
            end
            if Phbest(i,t)<0
                Phbest(i,t)=100;
            end    
        end
 end
 %检验功率平衡
for t=1:T
    for i=1:I
        Phsumbest(1,t)=Phsumbest(1,t)+Phbest(i,t);
    end
end
%计算Ps
 for t=1:T 
         for i=1:2
            Ps(t)=Pd(t)-Phsumbest(t); 
            Psbest(i,t)=Psmin(i)+rand()*(Psmax(i)-Psmin(i));
            Psum(1,t)=Psum(1,t)+Psbest(i,t);
         end
         Psbest(3,t)=Ps(t)-Psum(1,t);  
         if Psbest<0
             Psbest(3,t)=0;
         end
            if Psbest(3,t)<Psmin(3)
                 Psbest(3,t)=0;
                 Psbest(1,t)=Psbest(1,t)+ Psbest(3,t)./2;
                 Psbest(2,t)=Psbest(2,t)+ Psbest(3,t)./2;
             end
    end
      %检验功率平衡
for t=1:T
    for i=1:3
        Pssumbest(1,t)=Pssumbest(1,t)+Psbest(i,t);
    end
end
    

三、运行结果

四、备注

版本：2014a