【优化调度】杂交粒子群算法求解风光储微网日前优化调度问题【含Matlab源码 3872期】

✅博主简介：热爱科研的Matlab仿真开发者，修心和技术同步精进，Matlab项目合作可私信。
🍎个人主页：海神之光
🏆代码获取方式：
海神之光Matlab王者学习之路—代码获取方式
⛳️座右铭：行百里者，半于九十。

更多Matlab仿真内容点击👇
Matlab图像处理（进阶版）
路径规划（Matlab）
神经网络预测与分类（Matlab）
优化求解（Matlab）
语音处理（Matlab）
信号处理（Matlab）
车间调度（Matlab）

⛄一、杂交粒子群算法求解风光储微网日前优化调度问题简介

微网中包含：风电、光伏、储能、微型燃气轮机，以最小化电网购电成本、光伏风机的维护成本、蓄电池充放电维护成本、燃气轮机运行成本及污染气体治理成本为目标，综合考虑：功率平衡约束、燃气轮机爬坡约束、电网交换功率约束、储能装置约束、可控微电源出力约束等约束条件建立优化调度模型。
1.系统结构模型
1.1 光伏
光伏电池的空载电压与光谱辐照度有关，与电池面积无关。在100MW/cm2的太阳光谱辐照度下，单晶硅光伏电池的空载电压为450——600mV，最高可达690mV。当入射光谱辐照度变化时，光伏电池的空载电压与入射光谱辐照度的对数成正比。环境温度升高时，光伏电池的空载电压将下降，一般温度每升高1℃，Uoc约下降2——3mV。
由上述定义，可列出光伏电池等效电路中各变量的方程式如下：


式中，IO为光伏电池内部等效二极管PN结反向饱和电流，一般为常数；UD为等效二极管端电压；q为电子电荷；k为玻尔兹曼常量；T为热力学温度，A为PN结曲线常数。
光伏电池的伏安特性与光照强度和电池温度有关，可以认为光伏电池的数学物理模型方程以光照强度和电池温度作为参变量。以下是环境因素对光伏电池数学物理模型的修正。
当电池温度T、光照强度S不是参考光照强度和参考电池温度时，必须考虑它们的影响，并需要对光伏电池的数学物理模型进行修正。根据参考光照强度和参考电池温度下的Isc、Uoc、Im和Um,再将它们代入实用表达式，即可得到新的光照强度和新的电池温度下的伏安特性曲线。

1.2 风电
风力发电机的输出功率可以用下式表示：
式中：Pwt为风力发电机的输出功率，为空气密度，V为风速，R为风力发电机组的叶片半径，

Cp为风能利用系数。
风力发电机的出力和风速的关系为:

式中：Prated为额定功率；Vrated为额定风速，Vcutin为切入风速，Vcutout为切出风速

1.3 微型燃气轮机
选取燃氢微燃机作为氢电转换部分的能量耦合设备进行发电，该设备以氢⽓为燃料，燃烧产物只有水，实现零排放。

2 目标函数
目标函数为微网运行经济成本，主要包括：与主网交换功率成本，光伏、风机、蓄电池的维护成本，燃气轮机运行成本，污染气体治理成本，燃气轮机开机成本。

3 求解算法
采用PSO粒子群算法，用一种粒子来表示一个个体，每个粒子可视为N维搜索空间中的一个搜索个体，粒子的当前位置即为对应优化问题的一个候选解，粒子的飞行过程即为该个体的搜索过程．粒子的飞行速度可根据粒子历史最优位置和种群历史最优位置进行动态调整．粒子仅具有两个属性：速度和位置，速度代表移动的快慢，位置代表移动的方向。每个粒子单独搜寻的最优解叫做个体极值，粒子群中最优的个体极值作为当前全局最优解。不断迭代，更新速度和位置。最终得到满足终止条件的最优解。
算法流程如下：
①初始化
首先，我们设置最大迭代次数，目标函数的自变量个数，粒子的最大速度，位置信息为整个搜索空间，我们在速度区间和搜索空间上随机初始化速度和位置，设置粒子群规模，每个粒子随机初始化一个速度。
②个体极值与全局最优解定义适应度函数，个体极值为每个粒子找到的最优解，从这些最优解找到一个全局值，叫做本次全局最优解。与历史全局最优比较，进行更新。
③更新粒子的速度和位置

④终止条件
（1）达到设定迭代次数；（2）代数之间的差值满足最小界限。
本文在PSO粒子群算法的基础上进行改进，首先将传统固定惯性权重系数更改为根据全局最优点自适应调整惯性权重系数，即权重根据粒子的位置不同而动态变化。

⛄二、部分源代码

clear all;
close all;
clc;
%% INPUTS
pso_flag = 1;
swarm_size = 80;
BATTpmax = 400000;%电池出力上限400kw
BATTpmin = -150000;%电池出力下限-150KW
PVpmax =100000;%光伏出力上限100KW
WTpmax =100000;%风机出力上限100KW
MTpmax = 50000;%燃气轮机出力上限50KW
MTpmin=5000;%不停机运行的出力下限5kw
Putilmax=150000;%与主网交换功率最大值150kW
%% INIALIZE TIME STEPS AND AC BUS（初始化时间步长以及交流母线）
stp = 3600;
pertime = 24;
totalTime = pertime1; % total hours simulation
totalUnits = 5;%总共五个单元
p_bus = zeros(totalUnits,totalTime);%规格：（5，24）
index_GenMT = 1;
index_GenPV = 2;
index_GenWT = 3;
index_GenBTR = 4;
index_Load = 5;
index_Util = 6;
zeros24(1:pertime) = 0.5;
%% MICROURBINE PARAMETERS FOR EMS （微型燃气轮机参数）
%pMTmax = 2010^3;
cmt = 0.3571; % ($/m^3)fmt=0.0085;dt=1;koc=0.0587;$/Wh)Wh:瓦特小时
hot_startup = 30; %sec
cold_startup = 200; %sec
cooling_time = 520; %sec
mut = 600;% minimum up time (sec)
mdt = 300;% minimum down time (sec)
ton=0;% counting the times turned on
tonmin=2;
toff=0;% counting the times turned off
toffmin=2;
statusMT=ones(1,pertime);%先假设MT一直开机
Rup=75000;%爬坡上升速率限制 75kWh
Rdown=75000;%爬坡下降速率限制 75kWh

    %% WIND TURBINE PARAMETERS FOR EMS（风机参数）
Cwt = 0.0001;% $0.1/Wh
statusWT=ones(1,totalTime);
    %% BATTTERY PARAMETERS FOR EMS（蓄电池参数）
soc = zeros(1,totalTime);
soc(1,:) = 0.5; %initialize battery's sate of charge to 50%
statusBTR=ones(1,totalTime);
minchargpwr = 1000;%最小充电功率
maxchargpwr = 16000;%最大充电功率
effBTR = 0.9;%电池效率
    %% PV array PARAMETERS FOR EMS（光伏系统参数）
Cpv = 0.2;% $0.2/kWh
    %% LOAD PARAMETERS FOR EMS（负载侧参数）
powerDemandmin = 150*10^3;%负载需求最小值
powerDemandmax = 300*10^3;%负载需求最大值
   %% UTILITY
cutil = 1.05;
counterDisp=0;
util = zeros(1,totalTime);
stp_time =0;
irrad = zeros(totalTime);
tempr = zeros(totalTime);
vwind = zeros(totalTime);
p_min=zeros(3,pertime);
gb_mt=zeros(1,pertime);
gb_batt=zeros(1,pertime);
gb_util=zeros(1,pertime);
irrad=[0 0 0 0 0 0 0 0.2 0.3 0.4 0.6 0.8 0.8 0.8 0.8 0.7 0.6 0.4 0.3 0.1 0 0 0 0];%定义光照强度
tempr=[10 10 9 9 8 9 10 12 15 18 20 22 25 25 25 25 24 21 18 16 14 12 11 10];%定义温度
%vwind=[10 10 11 11 10 9 8 8 7 9 10 10 8 7 7 8 9 9 9 10 11 10 9 8];%定义风速，3-13
vwind=[10,11,12,12,12,12,12,13,15,16,17,17,18,19,19,20,21,22,23,24,25,26,27,28];%每小时的风速


 for i=1:totalTime
    %% LOAD DEMAND
    stp_time=stp_time+1;

    if(stp_time<=pertime)% stp_time<=24，分时负荷需求赋值
        if (stp_time==1 || stp_time==8 || stp_time==9 || stp_time==16 || stp_time==17 || stp_time==24||stp_time==23||stp_time==22)
            p_bus(index_Load,i) = 30000;
        elseif (stp_time==11 || stp_time==12 || stp_time==15 ||stp_time==19)
            p_bus(index_Load,i) = 60000;
        elseif (stp_time==10||stp_time==18)
            p_bus(index_Load,i) = 90000;
        elseif (stp_time==21)
            p_bus(index_Load,i) = 220000;
        elseif (stp_time==14 || stp_time==20)
            p_bus(index_Load,i) = 180000;
        elseif (stp_time==13)
            p_bus(index_Load,i) = 280000;
        else
            p_bus(index_Load,i) = 8000;
        end
    else %当运行时间大于24h
          if (mod(stp_time,pertime)==1 || mod(stp_time,pertime)==8 || mod(stp_time,pertime)==9 || mod(stp_time,pertime)==16 || mod(stp_time,pertime)==17 || mod(stp_time,pertime)==24)
            p_bus(index_Load,i) = 30000;
        elseif (mod(stp_time,pertime)==11 || mod(stp_time,pertime)==12 || mod(stp_time,pertime)==15 || mod(stp_time,pertime)==18 || mod(stp_time,pertime)==19 || mod(stp_time,pertime)==23)
            p_bus(index_Load,i) = 60000;
        elseif (mod(stp_time,pertime)==10)
            p_bus(index_Load,i) = 90000;
        elseif (mod(stp_time,pertime)==21 || mod(stp_time,pertime)==22)
            p_bus(index_Load,i) = 150000;
        elseif (mod(stp_time,pertime)==14 || mod(stp_time,pertime)==20)
            p_bus(index_Load,i) = 210000;
        elseif (mod(stp_time,pertime)==13)
            p_bus(index_Load,i) = 300000;
        else
            p_bus(index_Load,i) = 3000;
        end
    end
          %% CONSTRAINTS MT（微型燃气轮机的约束）***
   % if( ton >= mut/stp )%600/3600，将单位由s转换为h；MT打开时间大于上升时间约束
    %      statusMT(i) = 0;%标志MT关闭
    %  end

   %   if( toff >= mdt/stp )%MT关闭时间大于下降时间约束
   %       statusMT(i) = 1;%标志MT打开
   %   end

    %  if( statusMT(i) == 0 )
      %    ton=0;%开启时间记为0
     %     toff=toff+1;%关闭时间增加
      %    if(i<length(statusMT))
     %         statusMT(i+1)=0;
     %     end
    %  else
     %     toff = 0;
     %     ton = ton+1;
     %     if(i<length(statusMT))
      %        statusMT(i+1)=1;%当时间小于24h，标记下一个小时MT的状态为打开，顺延上一时刻的MT的状态
      %    end
     % end
            %% POWER GENERATORS (RES) 发出电能

% irrad(i) = 0.4;
% tempr(i) = 25 ;
% vwind(i) = 23;
% irrad(i) = rand();%当前光照强度，随机选取0-1
% tempr(i) = 25 + (75-25).rand();%当前温度**
% vwind(i) = 10 + (28-10).*rand();%当前风速
p_bus(index_GenWT,i) = windturbine(vwind(i),statusWT(i),WTpmax);%P_bus第三行记录风机各个小时的发电量
%if(mod(stp_time,pertime)<=19 && mod(stp_time,pertime)>=7)%一天中的7：00到19：00光伏出力
% p_bus(index_GenPV,i) = pv_array(irrad(i),tempr(i),PVpmax);%Pbus第二行记录光伏出力
%end
p_bus(index_GenPV,i) = pv_array(irrad(i),tempr(i),PVpmax);%Pbus第二行记录光伏出力

     if(mod(stp_time,pertime)==0)%到达第24h
        if(pso_flag==1)
            soc_tmp=soc;
            %% INITIALIZE POPULATION
            [ DV, putil, soc_ ,status_] = init_swarms(p_bus(index_GenWT,i-pertime+1:i),p_bus(index_GenPV,i-pertime+1:i),p_bus(index_Load,i-pertime+1:i),soc(i-pertime+1:i),statusMT(i-pertime+1:i), BATTpmax, BATTpmin, MTpmax,MTpmin,Putilmax,swarm_size);%初始化种群，因为此时i是24的整数倍，故输入参数是p_bus每行的第1到24列。
            %输出DV：蓄电池出力、燃气轮机出力；putil与电网交换功率；蓄电池soc值

⛄三、运行结果

⛄四、matlab版本及参考文献

1 matlab版本
2014a

2 参考文献
[1] 包子阳,余继周,杨杉.智能优化算法及其MATLAB实例（第2版）[M].电子工业出版社，2016.
[2]张岩,吴水根.MATLAB优化算法源代码[M].清华大学出版社，2017.

3 备注
简介此部分摘自互联网，仅供参考，若侵权，联系删除

🍅 仿真咨询
1 各类智能优化算法改进及应用
生产调度、经济调度、装配线调度、充电优化、车间调度、发车优化、水库调度、三维装箱、物流选址、货位优化、公交排班优化、充电桩布局优化、车间布局优化、集装箱船配载优化、水泵组合优化、解医疗资源分配优化、设施布局优化、可视域基站和无人机选址优化

2 机器学习和深度学习方面
卷积神经网络（CNN）、LSTM、支持向量机（SVM）、最小二乘支持向量机（LSSVM）、极限学习机（ELM）、核极限学习机（KELM）、BP、RBF、宽度学习、DBN、RF、RBF、DELM、XGBOOST、TCN实现风电预测、光伏预测、电池寿命预测、辐射源识别、交通流预测、负荷预测、股价预测、PM2.5浓度预测、电池健康状态预测、水体光学参数反演、NLOS信号识别、地铁停车精准预测、变压器故障诊断

3 图像处理方面
图像识别、图像分割、图像检测、图像隐藏、图像配准、图像拼接、图像融合、图像增强、图像压缩感知

4 路径规划方面
旅行商问题（TSP）、车辆路径问题（VRP、MVRP、CVRP、VRPTW等）、无人机三维路径规划、无人机协同、无人机编队、机器人路径规划、栅格地图路径规划、多式联运运输问题、车辆协同无人机路径规划、天线线性阵列分布优化、车间布局优化

5 无人机应用方面
无人机路径规划、无人机控制、无人机编队、无人机协同、无人机任务分配

6 无线传感器定位及布局方面
传感器部署优化、通信协议优化、路由优化、目标定位优化、Dv-Hop定位优化、Leach协议优化、WSN覆盖优化、组播优化、RSSI定位优化

7 信号处理方面
信号识别、信号加密、信号去噪、信号增强、雷达信号处理、信号水印嵌入提取、肌电信号、脑电信号、信号配时优化

8 电力系统方面
微电网优化、无功优化、配电网重构、储能配置

9 元胞自动机方面
交通流 人群疏散 病毒扩散 晶体生长

10 雷达方面
卡尔曼滤波跟踪、航迹关联、航迹融合